From ce8a510f65a50fca71d97133826801979401bb9e Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Wed, 27 Feb 2019 23:34:21 +1100
Subject: [PATCH] make time spacing finer in recovery
---
bilby/gw/conversion.py | 38 ++++++++++++++++++++++++--------------
1 file changed, 24 insertions(+), 14 deletions(-)
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index ffd10cc0d..1f4df1280 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -3,7 +3,7 @@ from __future__ import division
import numpy as np
from pandas import DataFrame
-from ..core.utils import logger, solar_mass
+from ..core.utils import logger, solar_mass, create_time_series
from ..core.prior import DeltaFunction, Interped
from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions
from .cosmology import get_cosmology
@@ -1043,6 +1043,8 @@ def generate_time_sample_from_marginalized_likelihood(
Generate a single sample from the posterior distribution for coalescence
time when using a likelihood which explicitly marginalises over time.
+ In order to resolve the posterior we artifically upsample to 16kHz.
+
See Eq. (C29-C32) of https://arxiv.org/abs/1809.02293
Parameters
@@ -1062,13 +1064,19 @@ def generate_time_sample_from_marginalized_likelihood(
if signal_polarizations is None:
signal_polarizations = \
likelihood.waveform_generator.frequency_domain_strain(sample)
- d_inner_h = np.zeros_like(likelihood.time_prior_array, dtype=np.complex)
+ n_time_steps = int(likelihood.waveform_generator.duration * 16384)
+ d_inner_h = np.zeros(n_time_steps, dtype=np.complex)
+ psd = np.ones(n_time_steps)
+ signal_long = np.zeros(n_time_steps, dtype=np.complex)
+ data = np.zeros(n_time_steps, dtype=np.complex)
rho_opt_sq = 0
for ifo in likelihood.interferometers:
+ ifo_length = len(ifo.frequency_domain_strain)
signal = ifo.get_detector_response(signal_polarizations, sample)
- d_inner_h += np.fft.fft(signal[0:-1] *
- ifo.frequency_domain_strain.conjugate()[0:-1] /
- ifo.power_spectral_density_array[0:-1])
+ signal_long[:ifo_length] = signal
+ data[:ifo_length] = np.conj(ifo.frequency_domain_strain)
+ psd[:ifo_length] = ifo.power_spectral_density_array
+ d_inner_h += np.fft.fft(signal_long * data / psd)
rho_opt_sq += ifo.optimal_snr_squared(signal=signal)
if likelihood.distance_marginalization:
@@ -1086,16 +1094,18 @@ def generate_time_sample_from_marginalized_likelihood(
else:
time_log_like = (d_inner_h.real - rho_opt_sq.real / 2)
- time_post = (np.exp(time_log_like - max(time_log_like)) *
- likelihood.time_prior_array)
+ times = create_time_series(
+ sampling_frequency=16384,
+ starting_time=likelihood.waveform_generator.start_time,
+ duration=likelihood.waveform_generator.duration)
- keep = (time_post > 1e-8)
- if sum(keep) == 1:
- sample['geocent_time'] = float(likelihood._times[keep])
- else:
- time_post = time_post[keep]
- times = likelihood._times[keep]
- sample['geocent_time'] = Interped(times, time_post).sample()
+ time_prior_array = likelihood.priors['geocent_time'].prob(times)
+ time_post = (np.exp(time_log_like - max(time_log_like)) * time_prior_array)
+
+ keep = (time_post > max(time_post) / 10)
+ time_post = time_post[keep]
+ times = times[keep]
+ sample['geocent_time'] = Interped(times, time_post).sample()
return sample
--
GitLab