Merge branch 'fix_ROQ_scaling' into 'master'

Fix ROQ scaling checks

Closes bilby_pipe#147

See merge request !678

bilby/gw/likelihood.py

@@ -835,8 +835,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
         If true, run tests using the roq_params to check the prior and data are
         valid for the ROQ
     roq_scale_factor: float
-        The ROQ scale factor used. WARNING: this does not apply the scaling,
-        but is only used for checking that the ROQ basis is appropriate.
+        The ROQ scale factor used.
     priors: dict, bilby.prior.PriorDict
         A dictionary of priors containing at least the geocent_time prior
     distance_marginalization_lookup_table: (dict, str), optional
@@ -1005,22 +1004,22 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             logger.info(msg)
 
         roq_params = self.roq_params
-        roq_params['flow'] *= self.roq_scale_factor
-        roq_params['fhigh'] *= self.roq_scale_factor
-        roq_params['seglen'] /= self.roq_scale_factor
-        roq_params['chirpmassmin'] /= self.roq_scale_factor
-        roq_params['chirpmassmax'] /= self.roq_scale_factor
-        roq_params['compmin'] /= self.roq_scale_factor
-
-        if ifo.maximum_frequency > roq_params['fhigh']:
+        roq_minimum_frequency = roq_params['flow'] * self.roq_scale_factor
+        roq_maximum_frequency = roq_params['fhigh'] * self.roq_scale_factor
+        roq_segment_length = roq_params['seglen'] / self.roq_scale_factor
+        roq_minimum_chirp_mass = roq_params['chirpmassmin'] / self.roq_scale_factor
+        roq_maximum_chirp_mass = roq_params['chirpmassmax'] / self.roq_scale_factor
+        roq_minimum_component_mass = roq_params['compmin'] / self.roq_scale_factor
+
+        if ifo.maximum_frequency > roq_maximum_frequency:
             raise BilbyROQParamsRangeError(
                 "Requested maximum frequency {} larger than ROQ basis fhigh {}"
-                .format(ifo.maximum_frequency, roq_params['fhigh']))
-        if ifo.minimum_frequency < roq_params['flow']:
+                .format(ifo.maximum_frequency, roq_maximum_frequency))
+        if ifo.minimum_frequency < roq_minimum_frequency:
             raise BilbyROQParamsRangeError(
                 "Requested minimum frequency {} lower than ROQ basis flow {}"
-                .format(ifo.minimum_frequency, roq_params['flow']))
-        if ifo.strain_data.duration != roq_params['seglen']:
+                .format(ifo.minimum_frequency, roq_minimum_frequency))
+        if ifo.strain_data.duration != roq_segment_length:
             raise BilbyROQParamsRangeError(
                 "Requested duration differs from ROQ basis seglen")
 
@@ -1031,27 +1030,27 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
         if priors.minimum_chirp_mass is None:
             logger.warning("Unable to check minimum chirp mass ROQ bounds")
-        elif priors.minimum_chirp_mass < roq_params["chirpmassmin"]:
+        elif priors.minimum_chirp_mass < roq_minimum_chirp_mass:
             raise BilbyROQParamsRangeError(
                 "Prior minimum chirp mass {} less than ROQ basis bound {}"
                 .format(priors.minimum_chirp_mass,
-                        roq_params["chirpmassmin"]))
+                        roq_minimum_chirp_mass))
 
         if priors.maximum_chirp_mass is None:
             logger.warning("Unable to check maximum_chirp mass ROQ bounds")
-        elif priors.maximum_chirp_mass > roq_params["chirpmassmax"]:
+        elif priors.maximum_chirp_mass > roq_maximum_chirp_mass:
             raise BilbyROQParamsRangeError(
                 "Prior maximum chirp mass {} greater than ROQ basis bound {}"
                 .format(priors.maximum_chirp_mass,
-                        roq_params["chirpmassmax"]))
+                        roq_maximum_chirp_mass))
 
         if priors.minimum_component_mass is None:
             logger.warning("Unable to check minimum component mass ROQ bounds")
-        elif priors.minimum_component_mass < roq_params["compmin"]:
+        elif priors.minimum_component_mass < roq_minimum_component_mass:
             raise BilbyROQParamsRangeError(
                 "Prior minimum component mass {} less than ROQ basis bound {}"
                 .format(priors.minimum_component_mass,
-                        roq_params["compmin"]))
+                        roq_minimum_component_mass))
 
     def _set_weights(self, linear_matrix, quadratic_matrix):
         """ Setup the time-dependent ROQ weights.
@@ -1077,11 +1076,15 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
         for ifo in self.interferometers:
             if self.roq_params is not None:
                 self.perform_roq_params_check(ifo)
+                # Get scaled ROQ quantities
+                roq_scaled_minimum_frequency = self.roq_params['flow'] * self.roq_scale_factor
+                roq_scaled_maximum_frequency = self.roq_params['fhigh'] * self.roq_scale_factor
+                roq_scaled_segment_length = self.roq_params['seglen'] * self.roq_scale_factor
                 # Generate frequencies for the ROQ
                 roq_frequencies = create_frequency_series(
-                    sampling_frequency=self.roq_params['fhigh'] * 2,
-                    duration=self.roq_params['seglen'])
-                roq_mask = roq_frequencies >= self.roq_params['flow']
+                    sampling_frequency=roq_scaled_maximum_frequency * 2,
+                    duration=roq_scaled_segment_length)
+                roq_mask = roq_frequencies >= roq_scaled_minimum_frequency
                 roq_frequencies = roq_frequencies[roq_mask]
                 overlap_frequencies, ifo_idxs, roq_idxs = np.intersect1d(
                     ifo.frequency_array[ifo.frequency_mask], roq_frequencies,

examples/gw_examples/injection_examples/roq_example.py

@@ -32,12 +32,11 @@ 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
+
+# Get scaled ROQ quantities
+minimum_chirp_mass = params['chirpmassmin'] / scale_factor
+maximum_chirp_mass = params['chirpmassmax'] / scale_factor
+minimum_component_mass = params['compmin'] / scale_factor
 
 np.random.seed(170808)
 
@@ -84,10 +83,10 @@ 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[key].minimum = max(priors[key].minimum, minimum_component_mass)
 priors['chirp_mass'] = bilby.core.prior.Constraint(
-    name='chirp_mass', minimum=float(params['chirpmassmin']),
-    maximum=float(params['chirpmassmax']))
+    name='chirp_mass', minimum=float(minimum_chirp_mass),
+    maximum=float(maximum_chirp_mass))
 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,
@@ -96,7 +95,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, roq_params=params)
+    priors=priors, roq_params=params, roq_scale_factor=scale_factor)
 
 # write the weights to file so they can be loaded multiple times
 likelihood.save_weights('weights.npz')

test/gw_likelihood_test.py

@@ -772,12 +772,6 @@ class TestRescaledROQLikelihood(unittest.TestCase):
 
         scale_factor = 0.5
         params = np.genfromtxt(self.params_file, 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
 
         self.duration = 4 / scale_factor
         self.sampling_frequency = 2048 * scale_factor
@@ -807,8 +801,9 @@ class TestRescaledROQLikelihood(unittest.TestCase):
 
         self.roq = bilby.gw.likelihood.ROQGravitationalWaveTransient(
             interferometers=ifos, waveform_generator=self.roq_wfg,
-            linear_matrix=linear_matrix_file, roq_params=params,
-            quadratic_matrix=quadratic_matrix_file, priors=self.priors)
+            linear_matrix=linear_matrix_file, roq_params=params, 
+            roq_scale_factor=scale_factor, quadratic_matrix=quadratic_matrix_file,
+            priors=self.priors)
 
 
 class TestBBHLikelihoodSetUp(unittest.TestCase):