From 6c03849f1734624ac01f99e029d6141f6320a332 Mon Sep 17 00:00:00 2001
From: Ethan Payne <ethan.payne@ligo.org>
Date: Fri, 12 Mar 2021 02:17:54 +0000
Subject: [PATCH] removed the likelihood matrix calculation within the
likelihood function, and instead use a loop in the post-processing script
---
bilby/core/result.py | 121 +++++--
bilby/gw/conversion.py | 17 +-
bilby/gw/detector/calibration.py | 82 +++++
bilby/gw/likelihood.py | 322 ++++++++++++++++--
containers/Dockerfile-test-suite-python37 | 76 +++++
.../calibration_marginalization_example.py | 136 ++++++++
requirements.txt | 1 +
test/gw/likelihood_test.py | 4 +-
8 files changed, 709 insertions(+), 50 deletions(-)
create mode 100644 containers/Dockerfile-test-suite-python37
create mode 100644 examples/gw_examples/injection_examples/calibration_marginalization_example.py
diff --git a/bilby/core/result.py b/bilby/core/result.py
index 7dcb4b04f..3e2ef48fe 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -5,6 +5,7 @@ from copy import copy
from distutils.version import LooseVersion
from importlib import import_module
from itertools import product
+from tqdm import tqdm
import corner
import h5py
@@ -103,7 +104,7 @@ def read_in_result(filename=None, outdir=None, label=None, extension='json', gzi
def get_weights_for_reweighting(
result, new_likelihood=None, new_prior=None, old_likelihood=None,
- old_prior=None):
+ old_prior=None, resume_file=None, n_checkpoint=5000):
""" Calculate the weights for reweight()
See bilby.core.result.reweight() for help with the inputs
@@ -116,17 +117,33 @@ def get_weights_for_reweighting(
An array of the natural-log likelihoods
new_log_prior_array: array
An array of the natural-log priors
+ resume_file: string
+ filepath for the resume file which stores the weights
+ n_checkpoint: int
+ Number of samples to reweight before writing a resume file
"""
+
nposterior = len(result.posterior)
- old_log_likelihood_array = np.zeros(nposterior)
- old_log_prior_array = np.zeros(nposterior)
- new_log_likelihood_array = np.zeros(nposterior)
- new_log_prior_array = np.zeros(nposterior)
- for ii, sample in result.posterior.iterrows():
+ if (resume_file is not None) and os.path.exists(resume_file):
+ old_log_likelihood_array, old_log_prior_array, new_log_likelihood_array, new_log_prior_array = \
+ np.genfromtxt(resume_file)
+
+ starting_index = np.argmin(np.abs(old_log_likelihood_array))
+ logger.info(f'Checkpoint resuming from {starting_index}.')
+
+ else:
+ old_log_likelihood_array = np.zeros(nposterior)
+ old_log_prior_array = np.zeros(nposterior)
+ new_log_likelihood_array = np.zeros(nposterior)
+ new_log_prior_array = np.zeros(nposterior)
+
+ starting_index = 0
+
+ for ii, sample in tqdm(result.posterior.iloc[starting_index:].iterrows()):
# Convert sample to dictionary
- par_sample = {key: sample[key] for key in result.search_parameter_keys}
+ par_sample = {key: sample[key] for key in result.posterior}
if old_likelihood is not None:
old_likelihood.parameters.update(par_sample)
@@ -152,10 +169,17 @@ def get_weights_for_reweighting(
# Don't perform prior reweighting (i.e. prior isn't updated)
new_log_prior_array[ii] = old_log_prior_array[ii]
+ if (ii % (n_checkpoint) == 0) and (resume_file is not None):
+ checkpointed_index = np.argmin(np.abs(old_log_likelihood_array))
+ logger.info(f'Checkpointing with {checkpointed_index} samples')
+ np.savetxt(
+ resume_file,
+ [old_log_likelihood_array, old_log_prior_array, new_log_likelihood_array, new_log_prior_array])
+
ln_weights = (
new_log_likelihood_array + new_log_prior_array - old_log_likelihood_array - old_log_prior_array)
- return ln_weights, new_log_likelihood_array, new_log_prior_array
+ return ln_weights, new_log_likelihood_array, new_log_prior_array, old_log_likelihood_array, old_log_prior_array
def rejection_sample(posterior, weights):
@@ -179,7 +203,9 @@ def rejection_sample(posterior, weights):
def reweight(result, label=None, new_likelihood=None, new_prior=None,
- old_likelihood=None, old_prior=None):
+ old_likelihood=None, old_prior=None, conversion_function=None, npool=1,
+ verbose_output=False, resume_file=None, n_checkpoint=5000,
+ use_nested_samples=False):
""" Reweight a result to a new likelihood/prior using rejection sampling
Parameters
@@ -198,6 +224,21 @@ def reweight(result, label=None, new_likelihood=None, new_prior=None,
old_prior: bilby.core.prior.PriorDict, (optional)
If given, calculate the old prior from this object. If not given,
the values stored in the posterior are used.
+ conversion_function: function, optional
+ Function which adds in extra parameters to the data frame,
+ should take the data_frame, likelihood and prior as arguments.
+ npool: int, optional
+ Number of threads with which to execute the conversion function
+ verbose_output: bool, optional
+ Flag determining whether the weight array and associated prior and
+ likelihood evaluations are output as well as the result file
+ resume_file: string, optional
+ filepath for the resume file which stores the weights
+ n_checkpoint: int, optional
+ Number of samples to reweight before writing a resume file
+ use_nested_samples: bool, optional
+ If true reweight the nested samples instead. This can greatly improve reweighting efficiency, especially if the
+ target distribution has support beyond the proposal posterior distribution.
Returns
=======
@@ -207,32 +248,56 @@ def reweight(result, label=None, new_likelihood=None, new_prior=None,
"""
result = copy(result)
+
+ if use_nested_samples:
+ result.posterior = result.nested_samples
+
nposterior = len(result.posterior)
logger.info("Reweighting posterior with {} samples".format(nposterior))
- ln_weights, new_log_likelihood_array, new_log_prior_array = get_weights_for_reweighting(
- result, new_likelihood=new_likelihood, new_prior=new_prior,
- old_likelihood=old_likelihood, old_prior=old_prior)
+ ln_weights, new_log_likelihood_array, new_log_prior_array, old_log_likelihood_array, old_log_prior_array =\
+ get_weights_for_reweighting(
+ result, new_likelihood=new_likelihood, new_prior=new_prior,
+ old_likelihood=old_likelihood, old_prior=old_prior,
+ resume_file=resume_file, n_checkpoint=n_checkpoint)
+
+ weights = np.exp(ln_weights)
+
+ if use_nested_samples:
+ weights *= result.posterior['weights']
# Overwrite the likelihood and prior evaluations
result.posterior["log_likelihood"] = new_log_likelihood_array
result.posterior["log_prior"] = new_log_prior_array
- weights = np.exp(ln_weights)
-
result.posterior = rejection_sample(result.posterior, weights=weights)
+ result.posterior = result.posterior.reset_index(drop=True)
logger.info("Rejection sampling resulted in {} samples".format(len(result.posterior)))
result.meta_data["reweighted_using_rejection_sampling"] = True
- result.log_evidence += logsumexp(ln_weights) - np.log(nposterior)
- result.priors = new_prior
+ if use_nested_samples:
+ result.log_evidence += np.log(np.sum(weights))
+ else:
+ result.log_evidence += logsumexp(ln_weights) - np.log(nposterior)
+
+ if conversion_function is not None:
+ data_frame = result.posterior
+ if "npool" in inspect.getargspec(conversion_function).args:
+ data_frame = conversion_function(data_frame, new_likelihood, new_prior, npool=npool)
+ else:
+ data_frame = conversion_function(data_frame, new_likelihood, new_prior)
+ result.posterior = data_frame
if label:
result.label = label
else:
result.label += "_reweighted"
- return result
+ if verbose_output:
+ return result, weights, new_log_likelihood_array, \
+ new_log_prior_array, old_log_likelihood_array, old_log_prior_array
+ else:
+ return result
class Result(object):
@@ -1398,7 +1463,7 @@ class Result(object):
if isinstance(self.injection_parameters.get(key, None), float)}
return credible_levels
- def get_injection_credible_level(self, parameter):
+ def get_injection_credible_level(self, parameter, weights=None):
"""
Get the credible level of the injected parameter
@@ -1415,11 +1480,15 @@ class Result(object):
if self.injection_parameters is None:
raise(TypeError, "Result object has no 'injection_parameters'. "
"Cannot copmute credible levels.")
+
+ if weights is None:
+ weights = np.ones(len(self.posterior))
+
if parameter in self.posterior and\
parameter in self.injection_parameters:
credible_level =\
- sum(self.posterior[parameter].values <
- self.injection_parameters[parameter]) / len(self.posterior)
+ sum(np.array(self.posterior[parameter].values <
+ self.injection_parameters[parameter]) * weights) / (sum(weights))
return credible_level
else:
return np.nan
@@ -1849,7 +1918,7 @@ def plot_multiple(results, filename=None, labels=None, colours=None,
@latex_plot_format
def make_pp_plot(results, filename=None, save=True, confidence_interval=[0.68, 0.95, 0.997],
lines=None, legend_fontsize='x-small', keys=None, title=True,
- confidence_interval_alpha=0.1,
+ confidence_interval_alpha=0.1, weight_list=None,
**kwargs):
"""
Make a P-P plot for a set of runs with injected signals.
@@ -1873,6 +1942,8 @@ def make_pp_plot(results, filename=None, save=True, confidence_interval=[0.68, 0
A list of keys to use, if None defaults to search_parameter_keys
confidence_interval_alpha: float, list, optional
The transparency for the background condifence interval
+ weight_list: list, optional
+ List of the weight arrays for each calculation.
kwargs:
Additional kwargs to pass to matplotlib.pyplot.plot
@@ -1886,10 +1957,14 @@ def make_pp_plot(results, filename=None, save=True, confidence_interval=[0.68, 0
if keys is None:
keys = results[0].search_parameter_keys
+ if weight_list is None:
+ weight_list = [None] * len(results)
+
credible_levels = pd.DataFrame()
- for result in results:
+ for i, result in enumerate(results):
credible_levels = credible_levels.append(
- result.get_all_injection_credible_levels(keys), ignore_index=True)
+ result.get_all_injection_credible_levels(keys, weights=weight_list[i]),
+ ignore_index=True)
if lines is None:
colors = ["C{}".format(i) for i in range(8)]
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index caf68fd26..f034dfe84 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -788,7 +788,8 @@ def _generate_all_cbc_parameters(sample, defaults, base_conversion,
if (
hasattr(likelihood, 'phase_marginalization') or
hasattr(likelihood, 'time_marginalization') or
- hasattr(likelihood, 'distance_marginalization')
+ hasattr(likelihood, 'distance_marginalization') or
+ hasattr(likelihood, 'calibration_marginalization')
):
try:
generate_posterior_samples_from_marginalized_likelihood(
@@ -1162,6 +1163,7 @@ def compute_snrs(sample, likelihood):
for ifo in likelihood.interferometers:
per_detector_snr = likelihood.calculate_snrs(
signal_polarizations, ifo)
+
matched_filter_snrs[ifo.name].append(
per_detector_snr.complex_matched_filter_snr)
optimal_snrs[ifo.name].append(
@@ -1201,7 +1203,8 @@ def generate_posterior_samples_from_marginalized_likelihood(
"""
if not any([likelihood.phase_marginalization,
likelihood.distance_marginalization,
- likelihood.time_marginalization]):
+ likelihood.time_marginalization,
+ likelihood.calibration_marginalization]):
return samples
# pass through a dictionary
@@ -1227,6 +1230,8 @@ def generate_posterior_samples_from_marginalized_likelihood(
samples['geocent_time'] = new_samples[:, 0]
samples['luminosity_distance'] = new_samples[:, 1]
samples['phase'] = new_samples[:, 2]
+ if likelihood.calibration_marginalization:
+ samples['recalib_index'] = new_samples[:, 3]
return samples
@@ -1254,4 +1259,10 @@ def fill_sample(args):
sample = dict(sample).copy()
likelihood.parameters.update(dict(sample).copy())
new_sample = likelihood.generate_posterior_sample_from_marginalized_likelihood()
- return new_sample["geocent_time"], new_sample["luminosity_distance"], new_sample["phase"]
+
+ if not likelihood.calibration_marginalization:
+ return new_sample["geocent_time"], new_sample["luminosity_distance"],\
+ new_sample["phase"]
+ else:
+ return new_sample["geocent_time"], new_sample["luminosity_distance"],\
+ new_sample["phase"], new_sample['recalib_index']
diff --git a/bilby/gw/detector/calibration.py b/bilby/gw/detector/calibration.py
index efb6336e4..88a5a664f 100644
--- a/bilby/gw/detector/calibration.py
+++ b/bilby/gw/detector/calibration.py
@@ -2,9 +2,91 @@
"""
import numpy as np
+import tables
from scipy.interpolate import interp1d
+def read_calibration_file(filename, frequency_array, number_of_response_curves, starting_index=0):
+ """
+ Function to read the hdf5 files from the calibration group containing the physical calibration response curves.
+
+ Parameters
+ ----------
+ filename: str
+ Location of the HDF5 file that contains the curves
+ frequency_array: array-like
+ The frequency values to calculate the calibration response curves at
+ number_of_response_curves: int
+ Number of random draws to use from the calibration file
+ starting_index: int
+ Index of the first curve to use within the array. This allows for segmenting the calibration curve array
+ into smaller pieces.
+
+ Returns
+ -------
+ calibration_draws: array-like
+ Array which contains the calibration responses as a function of the frequency array specified.
+ Shape is (number_of_response_curves x len(frequency_array))
+
+ """
+
+ calibration_file = tables.open_file(filename, 'r')
+ calibration_amplitude = \
+ calibration_file.root.deltaR.draws_amp_rel[starting_index:number_of_response_curves + starting_index]
+ calibration_phase = \
+ calibration_file.root.deltaR.draws_phase[starting_index:number_of_response_curves + starting_index]
+
+ calibration_frequencies = calibration_file.root.deltaR.freq[:]
+
+ calibration_file.close()
+
+ if len(calibration_amplitude.dtype) != 0: # handling if this is a calibration group hdf5 file
+ calibration_amplitude = calibration_amplitude.view(np.float64).reshape(calibration_amplitude.shape + (-1,))
+ calibration_phase = calibration_phase.view(np.float64).reshape(calibration_phase.shape + (-1,))
+ calibration_frequencies = calibration_frequencies.view(np.float64)
+
+ # interpolate to the frequency array (where if outside the range of the calibration uncertainty its fixed to 1)
+ calibration_draws = calibration_amplitude * np.exp(1j * calibration_phase)
+ calibration_draws = interp1d(
+ calibration_frequencies, calibration_draws, kind='cubic',
+ bounds_error=False, fill_value=1)(frequency_array)
+
+ return calibration_draws
+
+
+def write_calibration_file(filename, frequency_array, calibration_draws, calibration_parameter_draws=None):
+ """
+ Function to write the generated response curves to file
+
+ Parameters
+ ----------
+ filename: str
+ Location and filename to save the file
+ frequency_array: array-like
+ The frequency values where the calibration response was calculated
+ calibration_draws: array-like
+ Array which contains the calibration responses as a function of the frequency array specified.
+ Shape is (number_of_response_curves x len(frequency_array))
+ calibration_parameter_draws: data_frame
+ Parameters used to generate the random draws of the calibration response curves
+
+ """
+
+ calibration_file = tables.open_file(filename, 'w')
+ deltaR_group = calibration_file.create_group(calibration_file.root, 'deltaR')
+
+ # Save output
+ calibration_file.create_carray(deltaR_group, 'draws_amp_rel', obj=np.abs(calibration_draws))
+ calibration_file.create_carray(deltaR_group, 'draws_phase', obj=np.angle(calibration_draws))
+ calibration_file.create_carray(deltaR_group, 'freq', obj=frequency_array)
+
+ calibration_file.close()
+
+ # Save calibration parameter draws
+ if calibration_parameter_draws is not None:
+ calibration_parameter_draws.to_hdf(filename, key='CalParams', data_columns=True, format='table')
+
+
class Recalibrate(object):
name = 'none'
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index c33a67d23..0f1136839 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -7,6 +7,8 @@ import copy
import numpy as np
import scipy.integrate as integrate
from scipy.interpolate import interp1d
+import pandas as pd
+from tqdm import tqdm
try:
from scipy.special import logsumexp
@@ -19,8 +21,8 @@ from ..core.utils import BilbyJsonEncoder, decode_bilby_json
from ..core.utils import (
logger, UnsortedInterp2d, create_frequency_series, create_time_series,
speed_of_light, radius_of_earth)
-from ..core.prior import Interped, Prior, Uniform
-from .detector import InterferometerList, get_empty_interferometer
+from ..core.prior import Interped, Prior, Uniform, PriorDict, DeltaFunction
+from .detector import InterferometerList, get_empty_interferometer, calibration
from .prior import BBHPriorDict, CBCPriorDict, Cosmological
from .source import lal_binary_black_hole
from .utils import (
@@ -64,6 +66,9 @@ class GravitationalWaveTransient(Likelihood):
If true, marginalize over phase in the likelihood.
This is done analytically using a Bessel function.
The phase prior is set to be a delta function at phase=0.
+ calibration_marginalization: bool, optional
+ If true, marginalize over calibration response curves in the likelihood.
+ This is done numerically over a number of calibration response curve realizations.
priors: dict, optional
If given, used in the distance and phase marginalization.
distance_marginalization_lookup_table: (dict, str), optional
@@ -74,6 +79,19 @@ class GravitationalWaveTransient(Likelihood):
The lookup table is stored after construction in either the
provided string or a default location:
'.distance_marginalization_lookup_dmin{}_dmax{}_n{}.npz'
+ calibration_lookup_table: dict, optional
+ If a dict, contains the arrays over which to marginalize for each interferometer or the filepaths of the
+ calibration files.
+ If not provided, but calibration_marginalization is used, then the appropriate file is created to
+ contain the curves.
+ number_of_response_curves: int, optional
+ Number of curves from the calibration lookup table to use.
+ Default is 1000.
+ starting_index: int, optional
+ Sets the index for the first realization of the calibration curve to be considered.
+ This, coupled with number_of_response_curves, allows for restricting the set of curves used. This can be used
+ when dealing with large frequency arrays to split the calculation into sections.
+ Defaults to 0.
jitter_time: bool, optional
Whether to introduce a `time_jitter` parameter. This avoids either
missing the likelihood peak, or introducing biases in the
@@ -107,13 +125,16 @@ class GravitationalWaveTransient(Likelihood):
['d_inner_h',
'optimal_snr_squared',
'complex_matched_filter_snr',
+ 'd_inner_h_array',
+ 'optimal_snr_squared_array',
'd_inner_h_squared_tc_array'])
def __init__(
self, interferometers, waveform_generator, time_marginalization=False,
- distance_marginalization=False, phase_marginalization=False, priors=None,
- distance_marginalization_lookup_table=None, jitter_time=True,
- reference_frame="sky", time_reference="geocenter"
+ distance_marginalization=False, phase_marginalization=False, calibration_marginalization=False, priors=None,
+ distance_marginalization_lookup_table=None, calibration_lookup_table={},
+ number_of_response_curves=1000, starting_index=0, jitter_time=True, reference_frame="sky",
+ time_reference="geocenter"
):
self.waveform_generator = waveform_generator
@@ -122,6 +143,7 @@ class GravitationalWaveTransient(Likelihood):
self.time_marginalization = time_marginalization
self.distance_marginalization = distance_marginalization
self.phase_marginalization = phase_marginalization
+ self.calibration_marginalization = calibration_marginalization
self.priors = priors
self._check_set_duration_and_sampling_frequency_of_waveform_generator()
self.jitter_time = jitter_time
@@ -180,11 +202,19 @@ class GravitationalWaveTransient(Likelihood):
priors['luminosity_distance'] = float(self._ref_dist)
self._marginalized_parameters.append('luminosity_distance')
+ if self.calibration_marginalization:
+ self.number_of_response_curves = number_of_response_curves
+ self.starting_index = starting_index
+ self._setup_calibration_marginalization(calibration_lookup_table)
+ self._marginalized_parameters.append('recalib_index')
+
def __repr__(self):
return self.__class__.__name__ + '(interferometers={},\n\twaveform_generator={},\n\ttime_marginalization={}, ' \
- 'distance_marginalization={}, phase_marginalization={}, priors={})'\
+ 'distance_marginalization={}, phase_marginalization={}, '\
+ 'calibration_marginalization={}, priors={})'\
.format(self.interferometers, self.waveform_generator, self.time_marginalization,
- self.distance_marginalization, self.phase_marginalization, self.priors)
+ self.distance_marginalization, self.phase_marginalization, self.calibration_marginalization,
+ self.priors)
def _check_set_duration_and_sampling_frequency_of_waveform_generator(self):
""" Check the waveform_generator has the same duration and
@@ -221,23 +251,60 @@ class GravitationalWaveTransient(Likelihood):
"""
signal = interferometer.get_detector_response(
waveform_polarizations, self.parameters)
+ _mask = interferometer.frequency_mask
+
+ if 'recalib_index' in self.parameters:
+ signal[_mask] *= self.calibration_draws[interferometer.name][int(self.parameters['recalib_index'])]
+
d_inner_h = interferometer.inner_product(signal=signal)
optimal_snr_squared = interferometer.optimal_snr_squared(signal=signal)
complex_matched_filter_snr = d_inner_h / (optimal_snr_squared**0.5)
- if self.time_marginalization:
- d_inner_h_squared_tc_array =\
+ d_inner_h_array = None
+ optimal_snr_squared_array = None
+
+ if self.time_marginalization and self.calibration_marginalization:
+
+ d_inner_h_integrand = np.tile(
+ 4. / self.waveform_generator.duration *
+ interferometer.frequency_domain_strain.conjugate() * signal /
+ interferometer.power_spectral_density_array, (self.number_of_response_curves, 1)).T
+
+ d_inner_h_integrand[_mask] *= self.calibration_draws[interferometer.name].T
+
+ d_inner_h_array =\
+ 4 / self.waveform_generator.duration * np.fft.fft(
+ d_inner_h_integrand[0:-1], axis=0).T
+
+ optimal_snr_squared_integrand = 4. / self.waveform_generator.duration *\
+ np.abs(signal)**2 / interferometer.power_spectral_density_array
+ optimal_snr_squared_array = np.dot(optimal_snr_squared_integrand[_mask],
+ self.calibration_abs_draws[interferometer.name].T)
+
+ elif self.time_marginalization and not self.calibration_marginalization:
+ d_inner_h_array =\
4 / self.waveform_generator.duration * np.fft.fft(
signal[0:-1] *
interferometer.frequency_domain_strain.conjugate()[0:-1] /
interferometer.power_spectral_density_array[0:-1])
- else:
- d_inner_h_squared_tc_array = None
+
+ elif self.calibration_marginalization and ('recalib_index' not in self.parameters):
+ d_inner_h_integrand = 4. / self.waveform_generator.duration * \
+ interferometer.frequency_domain_strain.conjugate() * signal / \
+ interferometer.power_spectral_density_array
+ d_inner_h_array = np.dot(d_inner_h_integrand[_mask], self.calibration_draws[interferometer.name].T)
+
+ optimal_snr_squared_integrand = 4. / self.waveform_generator.duration *\
+ np.abs(signal)**2 / interferometer.power_spectral_density_array
+ optimal_snr_squared_array = np.dot(optimal_snr_squared_integrand[_mask],
+ self.calibration_abs_draws[interferometer.name].T)
return self._CalculatedSNRs(
d_inner_h=d_inner_h, optimal_snr_squared=optimal_snr_squared,
complex_matched_filter_snr=complex_matched_filter_snr,
- d_inner_h_squared_tc_array=d_inner_h_squared_tc_array)
+ d_inner_h_array=d_inner_h_array,
+ optimal_snr_squared_array=optimal_snr_squared_array,
+ d_inner_h_squared_tc_array=None)
def _check_marginalized_prior_is_set(self, key):
if key in self.priors and self.priors[key].is_fixed:
@@ -304,13 +371,27 @@ class GravitationalWaveTransient(Likelihood):
d_inner_h = 0.
optimal_snr_squared = 0.
complex_matched_filter_snr = 0.
- if self.time_marginalization:
+
+ if self.time_marginalization and self.calibration_marginalization:
+ if self.jitter_time:
+ self.parameters['geocent_time'] += self.parameters['time_jitter']
+
+ d_inner_h_array = np.zeros(
+ (self.number_of_response_curves, len(self.interferometers.frequency_array[0:-1])),
+ dtype=np.complex128)
+ optimal_snr_squared_array = np.zeros(self.number_of_response_curves, dtype=np.complex128)
+
+ elif self.time_marginalization:
if self.jitter_time:
self.parameters['geocent_time'] += self.parameters['time_jitter']
- d_inner_h_tc_array = np.zeros(
- self.interferometers.frequency_array[0:-1].shape,
+ d_inner_h_array = np.zeros(
+ len(self.interferometers.frequency_array[0:-1]),
dtype=np.complex128)
+ elif self.calibration_marginalization:
+ d_inner_h_array = np.zeros(self.number_of_response_curves, dtype=np.complex128)
+ optimal_snr_squared_array = np.zeros(self.number_of_response_curves, dtype=np.complex128)
+
for interferometer in self.interferometers:
per_detector_snr = self.calculate_snrs(
waveform_polarizations=waveform_polarizations,
@@ -320,12 +401,27 @@ class GravitationalWaveTransient(Likelihood):
optimal_snr_squared += np.real(per_detector_snr.optimal_snr_squared)
complex_matched_filter_snr += per_detector_snr.complex_matched_filter_snr
- if self.time_marginalization:
- d_inner_h_tc_array += per_detector_snr.d_inner_h_squared_tc_array
+ if self.time_marginalization or self.calibration_marginalization:
+ d_inner_h_array += per_detector_snr.d_inner_h_array
- if self.time_marginalization:
+ if self.calibration_marginalization:
+ optimal_snr_squared_array += per_detector_snr.optimal_snr_squared_array
+
+ if self.calibration_marginalization and self.time_marginalization:
+ log_l = self.time_and_calibration_marginalized_likelihood(
+ d_inner_h_array=d_inner_h_array,
+ h_inner_h=optimal_snr_squared_array)
+ if self.jitter_time:
+ self.parameters['geocent_time'] -= self.parameters['time_jitter']
+
+ elif self.calibration_marginalization:
+ log_l = self.calibration_marginalized_likelihood(
+ d_inner_h_calibration_array=d_inner_h_array,
+ h_inner_h=optimal_snr_squared_array)
+
+ elif self.time_marginalization:
log_l = self.time_marginalized_likelihood(
- d_inner_h_tc_array=d_inner_h_tc_array,
+ d_inner_h_tc_array=d_inner_h_array,
h_inner_h=optimal_snr_squared)
if self.jitter_time:
self.parameters['geocent_time'] -= self.parameters['time_jitter']
@@ -361,12 +457,22 @@ class GravitationalWaveTransient(Likelihood):
caching, as the signal is overwritten in place.
"""
if any([self.phase_marginalization, self.distance_marginalization,
- self.time_marginalization]):
+ self.time_marginalization, self.calibration_marginalization]):
signal_polarizations = copy.deepcopy(
self.waveform_generator.frequency_domain_strain(
self.parameters))
else:
return self.parameters
+
+ if self.calibration_marginalization and self.time_marginalization:
+ raise AttributeError(
+ "Cannot use time and calibration marginalization simultaneously for regeneration at the moment!"
+ "The matrix manipulation has not been tested.")
+
+ if self.calibration_marginalization:
+ new_calibration = self.generate_calibration_sample_from_marginalized_likelihood(
+ signal_polarizations=signal_polarizations)
+ self.parameters['recalib_index'] = new_calibration
if self.time_marginalization:
new_time = self.generate_time_sample_from_marginalized_likelihood(
signal_polarizations=signal_polarizations)
@@ -381,6 +487,38 @@ class GravitationalWaveTransient(Likelihood):
self.parameters['phase'] = new_phase
return self.parameters.copy()
+ def generate_calibration_sample_from_marginalized_likelihood(
+ self, signal_polarizations=None):
+ """
+ Generate a single sample from the posterior distribution for the set of calibration response curves when
+ explicitly marginalizing over the calibration uncertainty.
+
+ Parameters
+ ----------
+ signal_polarizations: dict, optional
+ Polarizations modes of the template.
+
+ Returns
+ -------
+ new_calibration: dict
+ Sample set from the calibration posterior
+ """
+ if 'recalib_index' in self.parameters:
+ self.parameters.pop('recalib_index')
+ self.parameters.update(self.get_sky_frame_parameters())
+ if signal_polarizations is None:
+ signal_polarizations = \
+ self.waveform_generator.frequency_domain_strain(self.parameters)
+
+ log_like = self.get_calibration_log_likelihoods(signal_polarizations=signal_polarizations)
+
+ calibration_post = np.exp(log_like - max(log_like))
+ calibration_post /= np.sum(calibration_post)
+
+ new_calibration = np.random.choice(self.number_of_response_curves, p=calibration_post)
+
+ return new_calibration
+
def generate_time_sample_from_marginalized_likelihood(
self, signal_polarizations=None):
"""
@@ -483,6 +621,7 @@ class GravitationalWaveTransient(Likelihood):
if signal_polarizations is None:
signal_polarizations = \
self.waveform_generator.frequency_domain_strain(self.parameters)
+
d_inner_h, h_inner_h = self._calculate_inner_products(signal_polarizations)
d_inner_h_dist = (
@@ -562,12 +701,17 @@ class GravitationalWaveTransient(Likelihood):
d_inner_h_ref = np.abs(d_inner_h_ref)
else:
d_inner_h_ref = np.real(d_inner_h_ref)
+
return self._interp_dist_margd_loglikelihood(
d_inner_h_ref, h_inner_h_ref)
def phase_marginalized_likelihood(self, d_inner_h, h_inner_h):
d_inner_h = self._bessel_function_interped(abs(d_inner_h))
- return d_inner_h - h_inner_h / 2
+
+ if self.calibration_marginalization and self.time_marginalization:
+ return d_inner_h - np.outer(h_inner_h, np.ones(np.shape(d_inner_h)[1])) / 2
+ else:
+ return d_inner_h - h_inner_h / 2
def time_marginalized_likelihood(self, d_inner_h_tc_array, h_inner_h):
if self.distance_marginalization:
@@ -585,6 +729,81 @@ class GravitationalWaveTransient(Likelihood):
time_prior_array = self.priors['geocent_time'].prob(times) * self._delta_tc
return logsumexp(log_l_tc_array, b=time_prior_array)
+ def time_and_calibration_marginalized_likelihood(self, d_inner_h_array, h_inner_h):
+ times = self._times
+ if self.jitter_time:
+ times = self._times + self.parameters['time_jitter']
+
+ _time_prior = self.priors['geocent_time']
+ time_mask = np.logical_and((times >= _time_prior.minimum), (times <= _time_prior.maximum))
+ times = times[time_mask]
+ time_probs = self.priors['geocent_time'].prob(times) * self._delta_tc
+
+ d_inner_h_array = d_inner_h_array[:, time_mask]
+ h_inner_h = h_inner_h
+
+ if self.distance_marginalization:
+ log_l_array = self.distance_marginalized_likelihood(
+ d_inner_h=d_inner_h_array, h_inner_h=h_inner_h)
+ elif self.phase_marginalization:
+ log_l_array = self.phase_marginalized_likelihood(
+ d_inner_h=d_inner_h_array,
+ h_inner_h=h_inner_h)
+ else:
+ log_l_array = np.real(d_inner_h_array) - np.outer(h_inner_h, np.ones(np.shape(d_inner_h_array)[1])) / 2
+
+ prior_array = np.outer(time_probs, 1. / self.number_of_response_curves * np.ones(len(h_inner_h))).T
+
+ return logsumexp(log_l_array, b=prior_array)
+
+ def get_calibration_log_likelihoods(self, signal_polarizations=None):
+ self.parameters.update(self.get_sky_frame_parameters())
+ if signal_polarizations is None:
+ signal_polarizations =\
+ self.waveform_generator.frequency_domain_strain(self.parameters)
+
+ d_inner_h = 0.
+ optimal_snr_squared = 0.
+ complex_matched_filter_snr = 0.
+ d_inner_h_array = np.zeros(self.number_of_response_curves, dtype=np.complex128)
+ optimal_snr_squared_array = np.zeros(self.number_of_response_curves, dtype=np.complex128)
+
+ for interferometer in self.interferometers:
+ per_detector_snr = self.calculate_snrs(
+ waveform_polarizations=signal_polarizations,
+ interferometer=interferometer)
+
+ d_inner_h += per_detector_snr.d_inner_h
+ optimal_snr_squared += np.real(per_detector_snr.optimal_snr_squared)
+ complex_matched_filter_snr += per_detector_snr.complex_matched_filter_snr
+ d_inner_h_array += per_detector_snr.d_inner_h_array
+ optimal_snr_squared_array += per_detector_snr.optimal_snr_squared_array
+
+ if self.distance_marginalization:
+ log_l_cal_array = self.distance_marginalized_likelihood(
+ d_inner_h=d_inner_h_array, h_inner_h=optimal_snr_squared_array)
+ elif self.phase_marginalization:
+ log_l_cal_array = self.phase_marginalized_likelihood(
+ d_inner_h=d_inner_h_array,
+ h_inner_h=optimal_snr_squared_array)
+ else:
+ log_l_cal_array = np.real(d_inner_h_array - optimal_snr_squared_array / 2)
+
+ return log_l_cal_array
+
+ def calibration_marginalized_likelihood(self, d_inner_h_calibration_array, h_inner_h):
+ if self.distance_marginalization:
+ log_l_cal_array = self.distance_marginalized_likelihood(
+ d_inner_h=d_inner_h_calibration_array, h_inner_h=h_inner_h)
+ elif self.phase_marginalization:
+ log_l_cal_array = self.phase_marginalized_likelihood(
+ d_inner_h=d_inner_h_calibration_array,
+ h_inner_h=h_inner_h)
+ else:
+ log_l_cal_array = np.real(d_inner_h_calibration_array - h_inner_h / 2)
+
+ return logsumexp(log_l_cal_array) - np.log(self.number_of_response_curves)
+
def _setup_rho(self, d_inner_h, optimal_snr_squared):
optimal_snr_squared_ref = (optimal_snr_squared.real *
self.parameters['luminosity_distance'] ** 2 /
@@ -735,6 +954,58 @@ class GravitationalWaveTransient(Likelihood):
self.time_prior_array = \
self.priors['geocent_time'].prob(self._times) * self._delta_tc
+ def _setup_calibration_marginalization(self, calibration_lookup_table):
+ self.calibration_draws = {}
+ self.calibration_abs_draws = {}
+ self.calibration_parameter_draws = {}
+ for interferometer in self.interferometers:
+
+ # Force the priors
+ calibration_priors = PriorDict()
+ for key in self.priors.keys():
+ if 'recalib' in key and interferometer.name in key:
+ calibration_priors[key] = copy.copy(self.priors[key])
+ self.priors[key] = DeltaFunction(0.0)
+
+ # If there is no entry in the lookup table, make an empty one
+ if interferometer.name not in calibration_lookup_table.keys():
+ calibration_lookup_table[interferometer.name] =\
+ f'{interferometer.name}_calibration_file.h5'
+
+ # If the interferometer lookup table file exists, generate the curves from it
+ if os.path.exists(calibration_lookup_table[interferometer.name]):
+ self.calibration_draws[interferometer.name] =\
+ calibration.read_calibration_file(
+ calibration_lookup_table[interferometer.name], self.interferometers.frequency_array,
+ self.number_of_response_curves, self.starting_index)
+
+ else: # generate the fake curves
+ self.calibration_parameter_draws[interferometer.name] =\
+ pd.DataFrame(calibration_priors.sample(self.number_of_response_curves))
+
+ self.calibration_draws[interferometer.name] = \
+ np.zeros((self.number_of_response_curves, len(interferometer.frequency_array)), dtype=complex)
+
+ for i in tqdm(range(self.number_of_response_curves)):
+ self.calibration_draws[interferometer.name][i, :] =\
+ interferometer.calibration_model.get_calibration_factor(
+ interferometer.frequency_array,
+ prefix='recalib_{}_'.format(interferometer.name),
+ **self.calibration_parameter_draws[interferometer.name].iloc[i])
+
+ calibration.write_calibration_file(
+ calibration_lookup_table[interferometer.name],
+ self.interferometers.frequency_array,
+ self.calibration_draws[interferometer.name],
+ self.calibration_parameter_draws[interferometer.name])
+
+ interferometer.calibration_model = calibration.Recalibrate()
+
+ _mask = interferometer.frequency_mask
+ self.calibration_draws[interferometer.name] = self.calibration_draws[interferometer.name][:, _mask]
+ self.calibration_abs_draws[interferometer.name] =\
+ np.abs(self.calibration_draws[interferometer.name])**2
+
@property
def interferometers(self):
return self._interferometers
@@ -821,6 +1092,7 @@ class GravitationalWaveTransient(Likelihood):
time_marginalization=self.time_marginalization,
phase_marginalization=self.phase_marginalization,
distance_marginalization=self.distance_marginalization,
+ calibration_marginalization=self.calibration_marginalization,
waveform_generator_class=self.waveform_generator.__class__,
waveform_arguments=self.waveform_generator.waveform_arguments,
frequency_domain_source_model=self.waveform_generator.frequency_domain_source_model,
@@ -1072,7 +1344,9 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
return self._CalculatedSNRs(
d_inner_h=np.nan_to_num(-np.inf), optimal_snr_squared=0,
complex_matched_filter_snr=np.nan_to_num(-np.inf),
- d_inner_h_squared_tc_array=None)
+ d_inner_h_squared_tc_array=None,
+ d_inner_h_array=None,
+ optimal_snr_squared_array=None)
d_inner_h_tc_array = np.einsum(
'i,ji->j', np.conjugate(h_plus_linear + h_cross_linear),
@@ -1092,7 +1366,9 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
return self._CalculatedSNRs(
d_inner_h=d_inner_h, optimal_snr_squared=optimal_snr_squared,
complex_matched_filter_snr=complex_matched_filter_snr,
- d_inner_h_squared_tc_array=d_inner_h_squared_tc_array)
+ d_inner_h_squared_tc_array=d_inner_h_squared_tc_array,
+ d_inner_h_array=None,
+ optimal_snr_squared_array=None)
@staticmethod
def _closest_time_indices(time, samples):
diff --git a/containers/Dockerfile-test-suite-python37 b/containers/Dockerfile-test-suite-python37
new file mode 100644
index 000000000..b1d3762e9
--- /dev/null
+++ b/containers/Dockerfile-test-suite-python37
@@ -0,0 +1,76 @@
+LABEL name="bilby Base Enterprise Linux 7" \
+maintainer="Gregory Ashton <gregory.ashton@ligo.org>" \
+date="20190104"
+
+ENV PATH /opt/conda/bin:$PATH
+
+# Install backend
+RUN apt-get update --fix-missing \
+&& apt-get install -y libglib2.0-0 libxext6 libsm6 libxrender1 libgl1-mesa-glx \
+dh-autoreconf build-essential libarchive-dev wget curl git libhdf5-serial-dev
+
+# Install python3.7
+RUN wget --quiet https://repo.anaconda.com/miniconda/Miniconda3-4.5.11-Linux-x86_64.sh -O ~/miniconda.sh && \
+/bin/bash ~/miniconda.sh -b -p /opt/conda && \
+rm ~/miniconda.sh && \
+/opt/conda/bin/conda clean -tipsy && \
+ln -s /opt/conda/etc/profile.d/conda.sh /etc/profile.d/conda.sh && \
+echo ". /opt/conda/etc/profile.d/conda.sh" >> ~/.bashrc && \
+echo "conda activate base" >> ~/.bashrc
+
+# Update pip and setuptools
+RUN pip install --upgrade pip \
+&& LC_ALL=C pip install --upgrade setuptools
+
+# Install conda-installable programs
+RUN conda install -y numpy scipy matplotlib pandas
+
+# Install conda-forge-installable programs
+RUN conda install -c conda-forge deepdish arviz
+
+# Install requirements
+RUN pip install future \
+pycondor>=0.5 \
+configargparse \
+flake8 \
+mock \
+pipenv \
+coverage \
+pytest-cov \
+coverage-badge
+
+# Install documentation requirements
+RUN pip install sphinx numpydoc nbsphinx sphinx_rtd_theme sphinx-tabs
+
+# Install additional bilby requirements
+RUN pip install corner lalsuite astropy gwpy theano healpy tables
+
+# Install samplers
+RUN pip install cpnest dynesty emcee nestle ptemcee pymc3 pymultinest kombine ultranest dnest4
+
+# Install pymultinest requirements
+RUN apt-get install -y libblas3 libblas-dev liblapack3 liblapack-dev \
+libatlas3-base libatlas-base-dev cmake build-essential gfortran
+
+# Install pymultinest
+RUN git clone https://github.com/farhanferoz/MultiNest.git \
+&& (cd MultiNest/MultiNest_v3.11_CMake/multinest && mkdir build && cd build && cmake .. && make)
+
+ENV LD_LIBRARY_PATH $HOME/MultiNest/MultiNest_v3.11_CMake/multinest/lib:
+
+# Install Polychord
+RUN git clone https://github.com/PolyChord/PolyChordLite.git \
+&& (cd PolyChordLite && python setup.py --no-mpi install)
+
+# Install PTMCMCSampler
+RUN git clone https://github.com/jellis18/PTMCMCSampler.git \
+&& (cd PTMCMCSampler && python setup.py install)
+
+# Add the ROQ data to the image
+RUN mkdir roq_basis \
+ && cd roq_basis \
+ && wget https://git.ligo.org/lscsoft/ROQ_data/raw/master/IMRPhenomPv2/4s/B_linear.npy \
+ && wget https://git.ligo.org/lscsoft/ROQ_data/raw/master/IMRPhenomPv2/4s/B_quadratic.npy \
+ && wget https://git.ligo.org/lscsoft/ROQ_data/raw/master/IMRPhenomPv2/4s/fnodes_linear.npy \
+ && wget https://git.ligo.org/lscsoft/ROQ_data/raw/master/IMRPhenomPv2/4s/fnodes_quadratic.npy \
+ && wget https://git.ligo.org/lscsoft/ROQ_data/raw/master/IMRPhenomPv2/4s/params.dat \
diff --git a/examples/gw_examples/injection_examples/calibration_marginalization_example.py b/examples/gw_examples/injection_examples/calibration_marginalization_example.py
new file mode 100644
index 000000000..73279aadb
--- /dev/null
+++ b/examples/gw_examples/injection_examples/calibration_marginalization_example.py
@@ -0,0 +1,136 @@
+#!/usr/bin/env python
+"""
+Tutorial to demonstrate running parameter estimation with calibration
+uncertainties marginalized over using a finite set of realizations.
+"""
+
+import numpy as np
+import bilby
+from copy import copy
+import matplotlib.pyplot as plt
+
+# Set the duration and sampling frequency of the data segment
+# that we're going to create and inject the signal into.
+duration = 4.
+sampling_frequency = 2048.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'calibration_marginalization'
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility. This is optional!
+np.random.seed(170817)
+
+# We are going to inject a binary black hole waveform. We first establish a
+# dictionary of parameters that includes all of the different waveform
+# parameters, including masses of the two black holes (mass_1, mass_2),
+# spins of both black holes (a, tilt, phi), etc.
+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=2000., theta_jn=0.4, psi=2.659,
+ phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+start_time = injection_parameters['geocent_time'] - duration + 2
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+ reference_frequency=50.)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+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)
+waveform_generator_rew = copy(waveform_generator)
+
+# Set up interferometers. In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
+# These default to their design sensitivity
+ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+for ifo in ifos:
+ injection_parameters.update({
+ 'recalib_{}_amplitude_{}'.format(ifo.name, ii): 0.0 for ii in range(10)})
+ injection_parameters.update({
+ 'recalib_{}_phase_{}'.format(ifo.name, ii): 0.0 for ii in range(10)})
+ ifo.calibration_model = bilby.gw.calibration.CubicSpline(
+ prefix='recalib_{}_'.format(ifo.name),
+ minimum_frequency=ifo.minimum_frequency,
+ maximum_frequency=ifo.maximum_frequency, n_points=10)
+ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
+ifos.inject_signal(parameters=injection_parameters,
+ waveform_generator=waveform_generator)
+ifos_rew = copy(ifos)
+
+# Set up prior, which is a dictionary
+# Here we fix the injected cbc parameters (except the distance)
+# to the injected values.
+priors = injection_parameters.copy()
+priors['luminosity_distance'] = bilby.prior.Uniform(
+ injection_parameters['luminosity_distance'] - 1000, injection_parameters['luminosity_distance'] + 1000,
+ name='luminosity_distance', latex_label='$d_L$')
+for key in injection_parameters:
+ if 'recalib' in key:
+ priors[key] = injection_parameters[key]
+
+# Convert to prior dictionary to replace the floats with delta function priors
+priors = bilby.core.prior.PriorDict(priors)
+priors_rew = copy(priors)
+
+# Initialise the likelihood by passing in the interferometer data (IFOs) and
+# the waveform generator. Here we assume there is no calibration uncertainty
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator, priors=priors)
+
+# Run sampler. In this case we're going to use the `dynesty` sampler
+result = bilby.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# Setting the log likelihood to actually be the log likelihood and not the log likelihood ratio...
+# This is used the for reweighting
+result.posterior['log_likelihood'] = result.posterior['log_likelihood'] + result.log_noise_evidence
+
+# Setting the priors we want on the calibration response curve parameters - as an example.
+for name in ['recalib_H1_amplitude_1', 'recalib_H1_amplitude_4']:
+ priors_rew[name] = bilby.prior.Gaussian(
+ mu=0, sigma=0.03, name=name, latex_label='H1 $A_{}$'.format(name[-1]))
+
+# Setting up the calibration marginalized likelihood.
+# We save the calibration response curve files into the output directory under {ifo.name}_calibration_file.h5
+cal_likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos_rew, waveform_generator=waveform_generator_rew,
+ calibration_marginalization=True, priors=priors_rew,
+ number_of_response_curves=100,
+ calibration_lookup_table={ifos[i].name: f'{outdir}/{ifos[i].name}_calibration_file.h5' for i in range(len(ifos))})
+
+# Plot the magnitude of the curves to be used in the marginalization
+plt.semilogx(ifos[0].frequency_array[ifos[0].frequency_mask][0:-1],
+ cal_likelihood.calibration_draws[ifos[0].name][:, 0:-1].T)
+plt.xlim(20, 1024)
+plt.ylabel('Magnitude')
+plt.xlabel('Frequency [Hz]')
+plt.savefig(f"{outdir}/calibration_draws.pdf")
+plt.clf()
+
+# Reweight the posterior samples from a distribution with no calibration uncertainty to one with uncertainty.
+# This method utilizes rejection sampling which can be inefficient at drawing samples at higher SNRs.
+result_rew = bilby.core.result.reweight(result,
+ new_likelihood=cal_likelihood,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+
+# Plot distance posterior with and without the calibration
+for res, label in zip([result, result_rew], ["No calibration uncertainty", "Calibration uncertainty"]):
+ plt.hist(res.posterior['luminosity_distance'], label=label, bins=50, histtype='step', density=True)
+plt.legend()
+plt.xlabel('Luminosity distance [Mpc]')
+plt.savefig(f"{outdir}/luminosity_distance_posterior.pdf")
+plt.clf()
+
+plt.hist(
+ result_rew.posterior['recalib_index'],
+ bins=np.linspace(0, cal_likelihood.number_of_response_curves - 1, cal_likelihood.number_of_response_curves),
+ density=True)
+plt.xlim(0, cal_likelihood.number_of_response_curves - 1)
+plt.xlabel('Calibration index')
+plt.savefig(f"{outdir}/calibration_index_histogram.pdf")
diff --git a/requirements.txt b/requirements.txt
index 50a9de636..05d86e2f5 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,3 +9,4 @@ mock
dill
tqdm
h5py
+tables
diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py
index 5e5c719e7..8f8bbcfc1 100644
--- a/test/gw/likelihood_test.py
+++ b/test/gw/likelihood_test.py
@@ -164,12 +164,13 @@ class TestGWTransient(unittest.TestCase):
expected = (
"GravitationalWaveTransient(interferometers={},\n\twaveform_generator={},\n\t"
"time_marginalization={}, distance_marginalization={}, phase_marginalization={}, "
- "priors={})".format(
+ "calibration_marginalization={}, priors={})".format(
self.interferometers,
self.waveform_generator,
False,
False,
False,
+ False,
self.prior,
)
)
@@ -211,6 +212,7 @@ class TestGWTransient(unittest.TestCase):
time_marginalization=False,
phase_marginalization=False,
distance_marginalization=False,
+ calibration_marginalization=False,
waveform_generator_class=self.waveform_generator.__class__,
waveform_arguments=self.waveform_generator.waveform_arguments,
frequency_domain_source_model=self.waveform_generator.frequency_domain_source_model,
--
GitLab