window_factor discrepancy with bilby fixed

Codes/main.py

@@ -262,7 +262,6 @@ if __name__ == '__main__':
 
             ifo.psd_array = ifo.power_spectral_density_array
             ifo.inverse_psd_array = 1 / ifo.psd_array
-
         # import IPython; IPython.embed();sys.exit()
 
 
@@ -346,6 +345,7 @@ if __name__ == '__main__':
 
     starting_point_parameters['geocent_duration'] = starting_point_parameters['geocent_time'] - start_time
 
+    # Note: Calling the `ifo.strain_data.frequency_domain_strain` for the first time sets the `ifo.strain_data.window_factor` to a value potentially != 1.
     for ifo in interferometers:
         ifo.fd_strain = ifo.strain_data.frequency_domain_strain
 
@@ -354,7 +354,6 @@ if __name__ == '__main__':
 
 
 
-
     # SET THE LIKELIHOOD OBJECT, two different cases whether ROQ basis is used or not
     if config_dict['analysis']['roq']:
         if ext_analysis_dict['approximant'] != 'IMRPhenomPv2':

Library/likelihood_gradient.py

@@ -374,8 +374,10 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
             ifo.exp_two_pi_f_dt_at_point = ifo.get_fd_time_translation(self.likelihood.parameters)
             ifo.template_at_point = ifo.get_detector_response(waveform_polarizations, self.likelihood.parameters, exp_two_pi_f_dt=ifo.exp_two_pi_f_dt_at_point)
             ifo.dh = {}
-            ifo.s_inner_h_at_point = nwip(ifo.fd_strain, ifo.template_at_point, ifo.psd_array, ifo.strain_data.duration)
-            ifo.h_inner_h_at_point = np.real(nwip(ifo.template_at_point, ifo.template_at_point, ifo.psd_array, ifo.strain_data.duration))
+            # ifo.s_inner_h_at_point = nwip(ifo.fd_strain, ifo.template_at_point, ifo.psd_array, ifo.strain_data.duration)
+            # ifo.h_inner_h_at_point = np.real(nwip(ifo.template_at_point, ifo.template_at_point, ifo.psd_array, ifo.strain_data.duration))
+            ifo.s_inner_h_at_point = ifo.inner_product(ifo.template_at_point)
+            ifo.h_inner_h_at_point = np.real(ifo.optimal_snr_squared(ifo.template_at_point))
             ifo.matched_filter_snr_at_point = np.real(ifo.s_inner_h_at_point) / ifo.h_inner_h_at_point**0.5
             self.s_inner_h_at_point += ifo.s_inner_h_at_point
             self.h_inner_h_at_point += ifo.h_inner_h_at_point
@@ -502,9 +504,12 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
         dh = (template_plus - template_minus) / (2 * offset)
         ifo.dh[traj_param_key] = dh
 
-        ifo.s_inner_h_plus = nwip(ifo.fd_strain, template_plus, ifo.psd_array, ifo.strain_data.duration)
-        ifo.s_inner_h_minus = nwip(ifo.fd_strain, template_minus, ifo.psd_array, ifo.strain_data.duration)
-        ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.s_inner_h_plus = nwip(ifo.fd_strain, template_plus, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.s_inner_h_minus = nwip(ifo.fd_strain, template_minus, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        ifo.s_inner_h_plus = ifo.inner_product(template_plus)
+        ifo.s_inner_h_minus = ifo.inner_product(template_minus)
+        ifo.h_inner_dh = ifo.noise_weighted_inner_product(ifo.template_at_point, dh)
 
         return dh
 
@@ -550,8 +555,10 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
         dh = (template_plus - ifo.template_at_point) / offset
         ifo.dh[traj_param_key] = dh
 
-        ifo.s_inner_h_plus = nwip(ifo.fd_strain, template_plus, ifo.psd_array, ifo.strain_data.duration)
-        ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.s_inner_h_plus = nwip(ifo.fd_strain, template_plus, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        ifo.s_inner_h_plus = ifo.inner_product(template_plus)
+        ifo.h_inner_dh = ifo.noise_weighted_inner_product(ifo.template_at_point, dh)
 
         return dh
 
@@ -595,8 +602,10 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
         dh = (ifo.template_at_point - template_minus) / offset
         ifo.dh[traj_param_key] = dh
 
-        ifo.s_inner_h_minus = nwip(ifo.fd_strain, template_minus, ifo.psd_array, ifo.strain_data.duration)
-        ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.s_inner_h_minus = nwip(ifo.fd_strain, template_minus, ifo.psd_array, ifo.strain_data.duration)
+        # ifo.h_inner_dh = nwip(ifo.template_at_point, dh, ifo.psd_array, ifo.strain_data.duration)
+        ifo.s_inner_h_minus = ifo.inner_product(template_minus)
+        ifo.h_inner_dh = ifo.noise_weighted_inner_product(ifo.template_at_point, dh)
 
         return dh
 
@@ -705,7 +714,8 @@ class GWTransientLikelihoodGradient(LikelihoodGradient):
             for j, param_key_col in enumerate(parameters_keys):
                 key_col = self.dict_trajectory_parameters_keys[param_key_col]
                 for ifo in self.interferometers:
-                    fisher_matrix[i][j] += np.real(nwip(ifo.dh[key_row], ifo.dh[key_col], ifo.psd_array, ifo.strain_data.duration))
+                    # fisher_matrix[i][j] += np.real(nwip(ifo.dh[key_row], ifo.dh[key_col], ifo.psd_array, ifo.strain_data.duration))
+                    fisher_matrix[i][j] += np.real(ifo.noise_weighted_inner_product(ifo.dh[key_row], ifo.dh[key_col]))
 
         return fisher_matrix
 

examples/config_files/config_GW190814_IMRPhenomHM.ini

@@ -2,7 +2,7 @@
 # Interferometers to run the analysis on.
 ifos = H1,L1,V1
 approximant = IMRPhenomHM
-minimum_frequency_ifos = {H1: 20.0, L1: 30.0, V1: 20.0}
+minimum_frequency_ifos = {H1: 20.0, L1: 20.0, V1: 20.0}
 maximum_frequency = 1024
 reference_frequency = 20
 roq = False
@@ -18,9 +18,9 @@ prior_file=../gw/prior_files/GW190814_alignedspins_withphic.prior
 
 [hmc]
 # Total number of trajectories of HMC
-n_traj_hmc_tot = 1500
+n_traj_hmc_tot = 15
 # Number of trajectories using  numericalgradients in phase 1.
-n_traj_fit = 1500
+n_traj_fit = 15
 # Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
 n_traj_for_this_run = 1000000
 # Number of points used to sub-select points in Ordered Look-Up Tables for local fit bimodal parameters like (cos(inc), psi, logD)

gw/detector/interferometer.py

@@ -5,6 +5,7 @@ sys.path.append('../')
 import numpy as np
 
 from bilby.gw.detector.interferometer import Interferometer as BilbyInterferometer
+from bilby.gw.utils import noise_weighted_inner_product as gwutilsnwip
 
 import Library.param_utils as paru
 import Library.bilby_utils as bilby_utils
@@ -91,3 +92,13 @@ class Interferometer(BilbyInterferometer):
         # signal_ifo[self.strain_data.frequency_mask] *= self.calibration_model.get_calibration_factor(self.strain_data.frequency_array[self.strain_data.frequency_mask], prefix='recalib_{}_'.format(self.name), **parameters)
 
         return signal_ifo
+
+    def noise_weighted_inner_product(self, aa, bb):
+        """
+        Allows to compute the nwip of any two signals in the interferometers. Existing methods in Bilby allow for <s|h> and <h|h> but we also need <h|dh> in the HMC.
+        """
+        return gwutilsnwip(
+            aa=aa[self.strain_data.frequency_mask],
+            bb=bb[self.strain_data.frequency_mask],
+            power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask],
+            duration=self.strain_data.duration)