make time spacing finer in recovery

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