Merge branch 'roq-frequency-matching' into 'master'

make frequency array checks more sensible in ROQ weight generation

See merge request !487

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)
 

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')

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):