Skip to content
Snippets Groups Projects
Commit 38b24605 authored by Colm Talbot's avatar Colm Talbot Committed by Paul Lasky
Browse files

Roq set time step

parent 6e249e8a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bilby/gw/likelihood.py
+ 83
− 27
View file @ 38b24605
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
......
test/gw_likelihood_test.py
+ 24
− 24
View file @ 38b24605
......@@ -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):
......
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