diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index dde655533ab342f58d0bd46feb8a4aa674b09a14..136d9567c832bcd9947421b08f65609d000706b0 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -17,7 +17,8 @@ from scipy.special import i0e
from ..core import likelihood
from ..core.utils import (
logger, UnsortedInterp2d, BilbyJsonEncoder, decode_bilby_json,
- create_frequency_series, create_time_series)
+ create_frequency_series, create_time_series, speed_of_light,
+ radius_of_earth)
from ..core.prior import Interped, Prior, Uniform
from .detector import InterferometerList
from .prior import BBHPriorDict
@@ -913,29 +914,38 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
return indices, in_bounds
def _set_weights(self, linear_matrix, quadratic_matrix):
- """
- Setup the time-dependent ROQ weights.
- This follows FIXME: Smith et al.
+ """ Setup the time-dependent ROQ weights.
+
+ Parameters
+ ----------
+ linear_matrix, quadratic_matrix: array_like
+ Arrays of the linear and quadratic basis
- The times are chosen to allow all the merger times allows in the time
- prior.
"""
- self.weights['time_samples'] = np.arange(
- self.priors['geocent_time'].minimum - 0.045,
- self.priors['geocent_time'].maximum + 0.045,
- self._get_time_resolution()) - self.interferometers.start_time
+ earth_light_crossing_time = radius_of_earth / speed_of_light
+ time_space = self._get_time_resolution()
+ delta_times = np.arange(
+ self.priors['geocent_time'].minimum - earth_light_crossing_time,
+ self.priors['geocent_time'].maximum + earth_light_crossing_time,
+ time_space)
+ time_samples = delta_times - self.interferometers.start_time
+ self.weights['time_samples'] = time_samples
+ logger.info("Using {} ROQ time samples".format(len(time_samples)))
for ifo in self.interferometers:
if self.roq_params is not None:
- frequencies = create_frequency_series(
- sampling_frequency=self.roq_params['fhigh'] * 2,
- duration=self.roq_params['seglen'])
- roq_mask = [frequencies >= self.roq_params['flow']]
- frequencies = frequencies[roq_mask]
+ if ifo.maximum_frequency > self.roq_params['fhigh']:
+ raise ValueError("Requested maximum frequency larger than ROQ basis fhigh")
+ # Generate frequencies for the ROQ
+ roq_frequencies = np.arange(
+ self.roq_params['flow'],
+ self.roq_params['fhigh'] + 1,
+ 1 / self.roq_params['seglen'])
overlap_frequencies, ifo_idxs, roq_idxs = np.intersect1d(
- ifo.frequency_array[ifo.frequency_mask], frequencies,
+ ifo.frequency_array[ifo.frequency_mask], roq_frequencies,
return_indices=True)
+
else:
overlap_frequencies = ifo.frequency_array[ifo.frequency_mask]
roq_idxs = np.arange(linear_matrix.shape[0], dtype=int)
@@ -950,11 +960,6 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
ifo.name, len(overlap_frequencies),
min(overlap_frequencies), max(overlap_frequencies)))
- # array of relative time shifts to be applied to the data
- # 0.045s comes from time for GW to traverse the Earth
- time_space = (self.weights['time_samples'][1] -
- self.weights['time_samples'][0])
-
# array to be filled with data, shifted by discrete time_samples
tc_shifted_data = np.zeros([
len(self.weights['time_samples']), len(overlap_frequencies)],
@@ -990,6 +995,8 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
quadratic_matrix[roq_idxs].real,
1 / ifo.strain_data.duration)
+ logger.info("Finished building weights for {}".format(ifo.name))
+
def save_weights(self, filename):
with open(filename, 'w') as file:
json.dump(self.weights, file, indent=2, cls=BilbyJsonEncoder)
@@ -1055,6 +1062,10 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
delta_t = fhigh**-1
+ # Apply a safety factor to ensure the time step is short enough
+ delta_t = delta_t / 5
+
+ logger.info("ROQ time-step = {}".format(delta_t))
return delta_t
def _rescale_signal(self, signal, new_distance):