diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index 409ca05d1f125fdcdcd452f53d964d945c55bac5..72966d35f7380973beca2276fe6ca7f7ad9e2eab 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -17,7 +17,7 @@ from scipy.special import i0e
from ..core import likelihood
from ..core.utils import (
logger, UnsortedInterp2d, BilbyJsonEncoder, decode_bilby_json,
- create_time_series)
+ create_frequency_series, create_time_series)
from ..core.prior import Interped, Prior, Uniform
from .detector import InterferometerList
from .prior import BBHPriorDict
@@ -912,17 +912,27 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
for ifo in self.interferometers:
if self.roq_params is not None:
- frequencies = np.arange(
- self.roq_params['flow'],
- self.roq_params['fhigh'] + 1 / self.roq_params['seglen'],
- 1 / self.roq_params['seglen'])
+ 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]
overlap_frequencies, ifo_idxs, roq_idxs = np.intersect1d(
ifo.frequency_array[ifo.frequency_mask], frequencies,
return_indices=True)
else:
overlap_frequencies = ifo.frequency_array[ifo.frequency_mask]
roq_idxs = np.arange(linear_matrix.shape[0], dtype=int)
- ifo_idxs = ifo.frequency_mask
+ ifo_idxs = [True] * sum(ifo.frequency_mask)
+ if sum(ifo_idxs) != len(roq_idxs):
+ raise ValueError(
+ "Mismatch between ROQ basis and frequency array for "
+ "{}".format(ifo.name))
+ logger.info(
+ "Building ROQ weights for {} with {} frequencies between {} "
+ "and {}.".format(
+ 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
@@ -936,7 +946,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
# shift data to beginning of the prior increment by the time step
shifted_data =\
- ifo.frequency_domain_strain[ifo_idxs] * \
+ ifo.frequency_domain_strain[ifo.frequency_mask][ifo_idxs] * \
np.exp(2j * np.pi * overlap_frequencies *
self.weights['time_samples'][0])
single_time_shift = np.exp(
@@ -951,14 +961,16 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
max_elements = int((max_block_gigabytes * 2 ** 30) / 8)
self.weights[ifo.name + '_linear'] = blockwise_dot_product(
- tc_shifted_data / ifo.power_spectral_density_array[ifo_idxs],
+ tc_shifted_data /
+ ifo.power_spectral_density_array[ifo.frequency_mask][ifo_idxs],
linear_matrix[roq_idxs],
max_elements) * 4 / ifo.strain_data.duration
del tc_shifted_data
self.weights[ifo.name + '_quadratic'] = build_roq_weights(
- 1 / ifo.power_spectral_density_array[ifo_idxs],
+ 1 /
+ ifo.power_spectral_density_array[ifo.frequency_mask][ifo_idxs],
quadratic_matrix[roq_idxs].real,
1 / ifo.strain_data.duration)
diff --git a/examples/injection_examples/roq_example.py b/examples/injection_examples/roq_example.py
index 6d54aaa99459347cc5d791f8798d0f0450903ed7..b98f4f2bfb279ce08dff05dc63cfaa66358e8242 100644
--- a/examples/injection_examples/roq_example.py
+++ b/examples/injection_examples/roq_example.py
@@ -17,19 +17,32 @@ import bilby
outdir = 'outdir'
label = 'roq'
+# The ROQ bases can be given an overall scaling.
+# Note how this is applied to durations, frequencies and masses.
+scale_factor = 1.6
+
# Load in the pieces for the linear part of the ROQ. Note you will need to
# adjust the filenames here to the correct paths on your machine
basis_matrix_linear = np.load("B_linear.npy").T
-freq_nodes_linear = np.load("fnodes_linear.npy")
+freq_nodes_linear = np.load("fnodes_linear.npy") * scale_factor
# Load in the pieces for the quadratic part of the ROQ
basis_matrix_quadratic = np.load("B_quadratic.npy").T
-freq_nodes_quadratic = np.load("fnodes_quadratic.npy")
+freq_nodes_quadratic = np.load("fnodes_quadratic.npy") * scale_factor
+
+# Load the parameters describing the valid parameters for the basis.
+params = np.genfromtxt("params.dat", names=True)
+params['flow'] *= scale_factor
+params['fhigh'] *= scale_factor
+params['seglen'] /= scale_factor
+params['chirpmassmin'] /= scale_factor
+params['chirpmassmax'] /= scale_factor
+params['compmin'] /= scale_factor
np.random.seed(170808)
-duration = 4
-sampling_frequency = 2048
+duration = 4 / scale_factor
+sampling_frequency = 2048 * scale_factor
injection_parameters = dict(
mass_1=36.0, mass_2=29.0, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
@@ -37,7 +50,8 @@ injection_parameters = dict(
phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=20., minimum_frequency=20.)
+ reference_frequency=20. * scale_factor,
+ minimum_frequency=20. * scale_factor)
waveform_generator = bilby.gw.WaveformGenerator(
duration=duration, sampling_frequency=sampling_frequency,
@@ -58,17 +72,22 @@ search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
frequency_domain_source_model=bilby.gw.source.roq,
waveform_arguments=dict(frequency_nodes_linear=freq_nodes_linear,
frequency_nodes_quadratic=freq_nodes_quadratic,
- reference_frequency=20., minimum_frequency=20.,
+ reference_frequency=20. * scale_factor,
+ minimum_frequency=20. * scale_factor,
approximant='IMRPhenomPv2'),
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
-# Here we add constraints on chirp mass and mass ratio to the prior
+# Here we add constraints on chirp mass and mass ratio to the prior, these are
+# determined by the domain of validity of the ROQ basis.
priors = bilby.gw.prior.BBHPriorDict()
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra',
'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
priors[key] = injection_parameters[key]
+for key in ['mass_1', 'mass_2']:
+ priors[key].minimum = max(priors[key].minimum, params['compmin'])
priors['chirp_mass'] = bilby.core.prior.Constraint(
- name='chirp_mass', minimum=12.5, maximum=45)
+ name='chirp_mass', minimum=float(params['chirpmassmin']),
+ maximum=float(params['chirpmassmax']))
priors['mass_ratio'] = bilby.core.prior.Constraint(0.125, 1, name='mass_ratio')
priors['geocent_time'] = bilby.core.prior.Uniform(
injection_parameters['geocent_time'] - 0.1,
@@ -77,7 +96,7 @@ priors['geocent_time'] = bilby.core.prior.Uniform(
likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=ifos, waveform_generator=search_waveform_generator,
linear_matrix=basis_matrix_linear, quadratic_matrix=basis_matrix_quadratic,
- priors=priors)
+ priors=priors, roq_params=params)
# write the weights to file so they can be loaded multiple times
likelihood.save_weights('weights.json')
diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py
index a1584bdc7c45a902b6573bb14403f4f8a25806b9..8447b8fd6b0dcd57ba0df4dd76c2f28e6e83facc 100644
--- a/test/gw_likelihood_test.py
+++ b/test/gw_likelihood_test.py
@@ -497,6 +497,9 @@ class TestROQLikelihood(unittest.TestCase):
fnodes_linear = np.load(fnodes_linear_file).T
fnodes_quadratic_file = "{}/fnodes_quadratic.npy".format(roq_dir)
fnodes_quadratic = np.load(fnodes_quadratic_file).T
+ self.linear_matrix_file = "{}/B_linear.npy".format(roq_dir)
+ self.quadratic_matrix_file = "{}/B_quadratic.npy".format(roq_dir)
+ self.params_file = "{}/params.dat".format(roq_dir)
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,
@@ -520,6 +523,8 @@ class TestROQLikelihood(unittest.TestCase):
ifos.inject_signal(
parameters=self.test_parameters, waveform_generator=non_roq_wfg)
+ self.ifos = ifos
+
roq_wfg = bilby.gw.waveform_generator.WaveformGenerator(
duration=self.duration, sampling_frequency=self.sampling_frequency,
frequency_domain_source_model=bilby.gw.source.roq,
@@ -529,6 +534,8 @@ class TestROQLikelihood(unittest.TestCase):
reference_frequency=20., minimum_frequency=20.,
approximant='IMRPhenomPv2'))
+ self.roq_wfg = roq_wfg
+
self.non_roq = bilby.gw.likelihood.GravitationalWaveTransient(
interferometers=ifos, waveform_generator=non_roq_wfg)
@@ -548,7 +555,8 @@ class TestROQLikelihood(unittest.TestCase):
phase_marginalization=True, priors=self.priors.copy())
def tearDown(self):
- pass
+ del (self.roq, self.non_roq, self.non_roq_phase, self.roq_phase,
+ self.ifos, self.priors)
def test_matches_non_roq(self):
self.non_roq.parameters.update(self.test_parameters)
@@ -573,6 +581,31 @@ class TestROQLikelihood(unittest.TestCase):
self.roq_phase.log_likelihood_ratio()) /
self.non_roq_phase.log_likelihood_ratio(), 1e-3)
+ def test_create_roq_weights_with_params(self):
+ roq = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+ interferometers=self.ifos, waveform_generator=self.roq_wfg,
+ linear_matrix=self.linear_matrix_file, roq_params=self.params_file,
+ quadratic_matrix=self.quadratic_matrix_file, priors=self.priors)
+ roq.parameters.update(self.test_parameters)
+ self.roq.parameters.update(self.test_parameters)
+ self.assertEqual(
+ roq.log_likelihood_ratio(), self.roq.log_likelihood_ratio())
+
+ def test_create_roq_weights_frequency_mismatch_works_with_params(self):
+ self.ifos[0].maximum_frequency = self.ifos[0].maximum_frequency / 2
+ _ = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+ interferometers=self.ifos, waveform_generator=self.roq_wfg,
+ linear_matrix=self.linear_matrix_file, roq_params=self.params_file,
+ quadratic_matrix=self.quadratic_matrix_file, priors=self.priors)
+
+ def test_create_roq_weights_frequency_mismatch_fails_without_params(self):
+ self.ifos[0].maximum_frequency = self.ifos[0].maximum_frequency / 2
+ with self.assertRaises(ValueError):
+ _ = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+ interferometers=self.ifos, waveform_generator=self.roq_wfg,
+ linear_matrix=self.linear_matrix_file,
+ quadratic_matrix=self.quadratic_matrix_file, priors=self.priors)
+
class TestBBHLikelihoodSetUp(unittest.TestCase):