diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index 43d7fcf7a83cdcd94158dbf0f246cea8a6fbf5ec..1e4a6cec9c14b02a28455e99cceb0ea187ecca2e 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -1,5 +1,6 @@
from __future__ import division
import numpy as np
+import scipy.integrate as integrate
from scipy.interpolate import interp1d
try:
@@ -16,7 +17,6 @@ from .prior import BBHPriorDict
from .source import lal_binary_black_hole
from .utils import noise_weighted_inner_product, build_roq_weights, blockwise_dot_product
from .waveform_generator import WaveformGenerator
-from math import ceil
class GravitationalWaveTransient(likelihood.Likelihood):
@@ -439,11 +439,6 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
optimal_snr_squared = 0.
d_inner_h = 0.
- indices, in_bounds = self._closest_time_indices(
- self.parameters['geocent_time'] - self.interferometers.start_time)
- if not in_bounds:
- return np.nan_to_num(-np.inf)
-
waveform = self.waveform_generator.frequency_domain_strain(
self.parameters)
if waveform is None:
@@ -469,7 +464,8 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
h_plus_quadratic = f_plus * waveform['quadratic']['plus']
h_cross_quadratic = f_cross * waveform['quadratic']['cross']
- indices, in_bounds = self._closest_time_indices(ifo_time)
+ indices, in_bounds = self._closest_time_indices(
+ ifo_time, self.time_samples[ifo.name])
if not in_bounds:
return np.nan_to_num(-np.inf)
@@ -478,8 +474,8 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
self.weights[ifo.name + '_linear'][indices])
d_inner_h += interp1d(
- self.time_samples[indices],
- d_inner_h_tc_array, kind='quadratic')(ifo_time)
+ self.time_samples[ifo.name][indices],
+ d_inner_h_tc_array, kind='cubic')(ifo_time)
optimal_snr_squared += \
np.vdot(np.abs(h_plus_quadratic + h_cross_quadratic)**2,
@@ -497,25 +493,28 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
return log_l.real
- def _closest_time_indices(self, time):
+ @staticmethod
+ def _closest_time_indices(time, samples):
"""
- Get the closest an two neighbouring times
+ Get the closest five times
Parameters
----------
time: float
Time to check
+ samples: array-like
+ Available times
Returns
-------
indices: list
- Indices nearest to time.
+ Indices nearest to time
in_bounds: bool
- Whether the indices are for valid times.
+ Whether the indices are for valid times
"""
- closest = np.argmin(np.absolute(self.time_samples - time))
- indices = [closest + ii for ii in [-1, 0, 1]]
- in_bounds = (indices[0] >= 0) & (indices[2] < self.time_samples.size)
+ closest = np.argmin(abs(samples - time))
+ indices = [closest + ii for ii in [-2, -1, 0, 1, 2]]
+ in_bounds = (indices[0] >= 0) & (indices[-1] < samples.size)
return indices, in_bounds
def _set_weights(self):
@@ -526,6 +525,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
The times are chosen to allow all the merger times allows in the time
prior.
"""
+ self.time_samples = dict()
for ifo in self.interferometers:
# only get frequency components up to maximum_frequency
self.linear_matrix = \
@@ -535,30 +535,28 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
# array of relative time shifts to be applied to the data
# 0.045s comes from time for GW to traverse the Earth
- self.time_samples = np.linspace(
+ self.time_samples[ifo.name] = np.arange(
self.priors['geocent_time'].minimum - 0.045,
self.priors['geocent_time'].maximum + 0.045,
- int(ceil((self.priors['geocent_time'].maximum -
- self.priors['geocent_time'].minimum + 0.09) *
- ifo.strain_data.sampling_frequency)))
- self.time_samples -= ifo.strain_data.start_time
- time_space = self.time_samples[1] - self.time_samples[0]
+ self._get_time_resolution(ifo))
+ self.time_samples[ifo.name] -= ifo.strain_data.start_time
+ time_space = (self.time_samples[ifo.name][1] -
+ self.time_samples[ifo.name][0])
# array to be filled with data, shifted by discrete time_samples
tc_shifted_data = np.zeros([
- len(self.time_samples),
+ len(self.time_samples[ifo.name]),
len(ifo.frequency_array[ifo.frequency_mask])], dtype=complex)
- # shift data to beginning of the prior
- # increment by the time step
+ # shift data to beginning of the prior increment by the time step
shifted_data =\
ifo.frequency_domain_strain[ifo.frequency_mask] * \
np.exp(2j * np.pi * ifo.frequency_array[ifo.frequency_mask] *
- self.time_samples[0])
+ self.time_samples[ifo.name][0])
single_time_shift = np.exp(
2j * np.pi * ifo.frequency_array[ifo.frequency_mask] *
time_space)
- for j in range(len(self.time_samples)):
+ for j in range(len(self.time_samples[ifo.name])):
tc_shifted_data[j] = shifted_data
shifted_data *= single_time_shift
@@ -578,6 +576,64 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
1 / ifo.power_spectral_density_array[ifo.frequency_mask],
self.quadratic_matrix.real, 1 / ifo.strain_data.duration)
+ @staticmethod
+ def _get_time_resolution(ifo):
+ """
+ This method estimates the time resolution given the optimal SNR of the
+ signal in the detector. This is then used when constructing the weights
+ for the ROQ.
+
+ Parameters
+ ----------
+ ifo: bilby.gw.detector.Interferometer
+
+ Returns
+ -------
+ delta_t: float
+ Time resolution
+ """
+
+ def calc_fhigh(freq, psd, scaling=20.):
+ """
+
+ Parameters
+ ----------
+ freq: array-like
+ Frequency array
+ psd: array-like
+ Power spectral density
+ scaling: float
+ SNR dependent scaling factor
+
+ Returns
+ -------
+ f_high: float
+ The maximum frequency which must be considered
+ """
+ integrand1 = np.power(freq, -7. / 3) / psd
+ integral1 = integrate.simps(integrand1, freq)
+ integrand3 = np.power(freq, 2. / 3.) / (psd * integral1)
+ f_3_bar = integrate.simps(integrand3, freq)
+
+ f_high = scaling * f_3_bar**(1 / 3)
+
+ return f_high
+
+ def c_f_scaling(snr):
+ return (np.pi**2 * snr**2 / 6)**(1 / 3)
+
+ psd = ifo.power_spectral_density_array[ifo.frequency_mask]
+
+ freq = ifo.frequency_array[ifo.frequency_mask]
+
+ inj_snr = getattr(ifo.meta_data, 'optimal_SNR', 30)
+
+ fhigh = calc_fhigh(freq, psd, scaling=c_f_scaling(inj_snr))
+
+ delta_t = fhigh**-1
+
+ return delta_t
+
def get_binary_black_hole_likelihood(interferometers):
""" A rapper to quickly set up a likelihood for BBH parameter estimation
diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py
index 632aadac94536dcd425e7b2914ef1f63bcf17a94..0dbfd147eb7b452114ed62b02f7cfe4ccca7c4c7 100644
--- a/test/gw_likelihood_test.py
+++ b/test/gw_likelihood_test.py
@@ -473,7 +473,7 @@ class TestROQLikelihood(unittest.TestCase):
self.test_parameters = dict(
mass_1=36.0, mass_2=36.0, a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0,
- phi_12=1.7, phi_jl=0.3, luminosity_distance=5000., theta_jn=0.4,
+ phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4,
psi=0.659, phase=1.3, geocent_time=1.2, ra=1.3, dec=-1.2)
ifos = bilby.gw.detector.InterferometerList(['H1'])
@@ -481,7 +481,7 @@ class TestROQLikelihood(unittest.TestCase):
sampling_frequency=self.sampling_frequency, duration=self.duration)
self.priors = bilby.gw.prior.BBHPriorDict()
- self.priors['geocent_time'] = bilby.core.prior.Uniform(1.1, 1.3)
+ self.priors['geocent_time'] = bilby.core.prior.Uniform(1.19, 1.21)
non_roq_wfg = bilby.gw.WaveformGenerator(
duration=self.duration, sampling_frequency=self.sampling_frequency,
@@ -502,15 +502,19 @@ class TestROQLikelihood(unittest.TestCase):
reference_frequency=20., minimum_frequency=20.,
approximant='IMRPhenomPv2'))
- self.non_roq_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+ self.non_roq = bilby.gw.likelihood.GravitationalWaveTransient(
interferometers=ifos, waveform_generator=non_roq_wfg)
- self.roq_likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+ self.non_roq_phase = bilby.gw.likelihood.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=non_roq_wfg,
+ phase_marginalization=True, priors=self.priors.copy())
+
+ self.roq = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=ifos, waveform_generator=roq_wfg,
linear_matrix=linear_matrix_file,
quadratic_matrix=quadratic_matrix_file, priors=self.priors)
- self.roq_phase_like = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+ self.roq_phase = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=ifos, waveform_generator=roq_wfg,
linear_matrix=linear_matrix_file,
quadratic_matrix=quadratic_matrix_file,
@@ -520,31 +524,27 @@ class TestROQLikelihood(unittest.TestCase):
pass
def test_matches_non_roq(self):
- self.non_roq_likelihood.parameters.update(self.test_parameters)
- self.roq_likelihood.parameters.update(self.test_parameters)
- self.assertAlmostEqual(
- self.non_roq_likelihood.log_likelihood_ratio(),
- self.roq_likelihood.log_likelihood_ratio(), 0)
+ self.non_roq.parameters.update(self.test_parameters)
+ self.roq.parameters.update(self.test_parameters)
+ self.assertLess(
+ abs(self.non_roq.log_likelihood_ratio() -
+ self.roq.log_likelihood_ratio()) /
+ self.non_roq.log_likelihood_ratio(), 1e-3)
def test_time_prior_out_of_bounds_returns_zero(self):
- self.roq_likelihood.parameters.update(self.test_parameters)
- self.roq_likelihood.parameters['geocent_time'] = -5
+ self.roq.parameters.update(self.test_parameters)
+ self.roq.parameters['geocent_time'] = -5
self.assertEqual(
- self.roq_likelihood.log_likelihood_ratio(), np.nan_to_num(-np.inf))
+ self.roq.log_likelihood_ratio(), np.nan_to_num(-np.inf))
def test_phase_marginalisation_roq(self):
"""Test phase marginalised likelihood matches brute force version"""
- like = []
- self.roq_likelihood.parameters.update(self.test_parameters)
- phases = np.linspace(0, 2 * np.pi, 1000)
- for phase in phases:
- self.roq_likelihood.parameters['phase'] = phase
- like.append(np.exp(self.roq_likelihood.log_likelihood_ratio()))
-
- marg_like = np.log(np.trapz(like, phases) / (2 * np.pi))
- self.roq_phase_like.parameters = self.test_parameters.copy()
- self.assertAlmostEqual(
- marg_like, self.roq_phase_like.log_likelihood_ratio(), delta=0.5)
+ self.non_roq_phase.parameters = self.test_parameters.copy()
+ self.roq_phase.parameters = self.test_parameters.copy()
+ self.assertLess(
+ abs(self.non_roq_phase.log_likelihood_ratio() -
+ self.roq_phase.log_likelihood_ratio()) /
+ self.non_roq_phase.log_likelihood_ratio(), 1e-3)
class TestBBHLikelihoodSetUp(unittest.TestCase):