diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index d6a5cde5ec0bdf3e48dd1ac4ed24c90ac7c5f24b..65d9b12d9294689e8b1047aa1fbc0dd572d74ec2 100644 --- a/bilby/gw/conversion.py +++ b/bilby/gw/conversion.py @@ -933,7 +933,7 @@ def compute_snrs(sample, likelihood): for ifo in likelihood.interferometers: signal = ifo.get_detector_response(signal_polarizations, sample) sample['{}_matched_filter_snr'.format(ifo.name)] =\ - ifo.matched_filter_snr_squared(signal=signal) ** 0.5 + ifo.matched_filter_snr(signal=signal) sample['{}_optimal_snr'.format(ifo.name)] = \ ifo.optimal_snr_squared(signal=signal) ** 0.5 else: @@ -950,7 +950,7 @@ def compute_snrs(sample, likelihood): signal = ifo.get_detector_response( signal_polarizations, sample.iloc[ii]) matched_filter_snrs[ifo.name].append( - ifo.matched_filter_snr_squared(signal=signal) ** 0.5) + ifo.matched_filter_snr(signal=signal)) optimal_snrs[ifo.name].append( ifo.optimal_snr_squared(signal=signal) ** 0.5) @@ -1023,20 +1023,20 @@ def _generate_distance_sample_from_marginalized_likelihood(sample, likelihood): """ signal_polarizations = \ likelihood.waveform_generator.frequency_domain_strain(sample) - rho_mf_sq = 0 + d_inner_h = 0 rho_opt_sq = 0 for ifo in likelihood.interferometers: signal = ifo.get_detector_response(signal_polarizations, sample) - rho_mf_sq += ifo.matched_filter_snr_squared(signal=signal) + d_inner_h += ifo.inner_product(signal=signal) rho_opt_sq += ifo.optimal_snr_squared(signal=signal) - rho_mf_sq_dist = (rho_mf_sq * sample['luminosity_distance'] / + d_inner_h_dist = (d_inner_h * sample['luminosity_distance'] / likelihood._distance_array) rho_opt_sq_dist = (rho_opt_sq * sample['luminosity_distance']**2 / likelihood._distance_array**2) - distance_log_like = (rho_mf_sq_dist.real - rho_opt_sq_dist.real / 2) + distance_log_like = (d_inner_h_dist.real - rho_opt_sq_dist.real / 2) distance_post = np.exp(distance_log_like - max(distance_log_like)) *\ likelihood.distance_prior_array diff --git a/bilby/gw/detector.py b/bilby/gw/detector.py index 22845a6d283bdb3a2c34bbb295ee16058f3f0c26..217fa143b3f2a361e789d10950f4b338548cec67 100644 --- a/bilby/gw/detector.py +++ b/bilby/gw/detector.py @@ -1357,7 +1357,7 @@ class Interferometer(object): self.meta_data['optimal_SNR'] = ( np.sqrt(self.optimal_snr_squared(signal=signal_ifo)).real) self.meta_data['matched_filter_SNR'] = ( - np.sqrt(self.matched_filter_snr_squared(signal=signal_ifo))) + self.matched_filter_snr(signal=signal_ifo)) self.meta_data['parameters'] = parameters logger.info("Injected signal in {}:".format(self.name)) @@ -1499,11 +1499,29 @@ class Interferometer(object): ------- float: The optimal signal to noise ratio possible squared """ - return gwutils.optimal_snr_squared(signal=signal, - power_spectral_density=self.power_spectral_density_array, - duration=self.strain_data.duration) + return gwutils.optimal_snr_squared( + signal=signal, + power_spectral_density=self.power_spectral_density_array, + duration=self.strain_data.duration) - def matched_filter_snr_squared(self, signal): + def inner_product(self, signal): + """ + + Parameters + ---------- + signal: array_like + Array containing the signal + + Returns + ------- + float: The optimal signal to noise ratio possible squared + """ + return gwutils.noise_weighted_inner_product( + aa=signal, bb=self.frequency_array, + power_spectral_density=self.power_spectral_density_array, + duration=self.strain_data.duration) + + def matched_filter_snr(self, signal): """ Parameters @@ -1516,10 +1534,10 @@ class Interferometer(object): float: The matched filter signal to noise ratio squared """ - return gwutils.matched_filter_snr_squared(signal=signal, - frequency_domain_strain=self.frequency_domain_strain, - power_spectral_density=self.power_spectral_density_array, - duration=self.strain_data.duration) + return gwutils.matched_filter_snr( + signal=signal, frequency_domain_strain=self.frequency_domain_strain, + power_spectral_density=self.power_spectral_density_array, + duration=self.strain_data.duration) @property def whitened_frequency_domain_strain(self): diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py index b56a377c603764b4e7853e9a34e438f39161f037..0e18210cb60d83447a946fb4a81883048dc763bd 100644 --- a/bilby/gw/likelihood.py +++ b/bilby/gw/likelihood.py @@ -156,19 +156,19 @@ class GravitationalWaveTransient(likelihood.Likelihood): if waveform_polarizations is None: return np.nan_to_num(-np.inf) - matched_filter_snr_squared = 0 + d_inner_h = 0 optimal_snr_squared = 0 - matched_filter_snr_squared_tc_array = np.zeros( + d_inner_h_squared_tc_array = np.zeros( self.interferometers.frequency_array[0:-1].shape, dtype=np.complex128) for interferometer in self.interferometers: signal_ifo = interferometer.get_detector_response( waveform_polarizations, self.parameters) - matched_filter_snr_squared += interferometer.matched_filter_snr_squared(signal=signal_ifo) + d_inner_h += interferometer.inner_product(signal=signal_ifo) optimal_snr_squared += interferometer.optimal_snr_squared(signal=signal_ifo) if self.time_marginalization: - matched_filter_snr_squared_tc_array +=\ + d_inner_h_squared_tc_array +=\ 4 / self.waveform_generator.duration * np.fft.fft( signal_ifo[0:-1] * interferometer.frequency_domain_strain.conjugate()[0:-1] / @@ -178,7 +178,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): if self.distance_marginalization: rho_mf_ref_tc_array, rho_opt_ref = self._setup_rho( - matched_filter_snr_squared_tc_array, optimal_snr_squared) + d_inner_h_squared_tc_array, optimal_snr_squared) if self.phase_marginalization: dist_marged_log_l_tc_array = self._interp_dist_margd_loglikelihood( abs(rho_mf_ref_tc_array), rho_opt_ref) @@ -191,34 +191,34 @@ class GravitationalWaveTransient(likelihood.Likelihood): b=self.time_prior_array) elif self.phase_marginalization: log_l = logsumexp(self._bessel_function_interped(abs( - matched_filter_snr_squared_tc_array)), + d_inner_h_squared_tc_array)), b=self.time_prior_array) - optimal_snr_squared / 2 else: log_l = logsumexp( - matched_filter_snr_squared_tc_array.real, + d_inner_h_squared_tc_array.real, b=self.time_prior_array) - optimal_snr_squared / 2 elif self.distance_marginalization: - rho_mf_ref, rho_opt_ref = self._setup_rho(matched_filter_snr_squared, optimal_snr_squared) + rho_mf_ref, rho_opt_ref = self._setup_rho(d_inner_h, optimal_snr_squared) if self.phase_marginalization: rho_mf_ref = abs(rho_mf_ref) log_l = self._interp_dist_margd_loglikelihood(rho_mf_ref.real, rho_opt_ref)[0] elif self.phase_marginalization: - matched_filter_snr_squared = self._bessel_function_interped(abs(matched_filter_snr_squared)) - log_l = matched_filter_snr_squared - optimal_snr_squared / 2 + d_inner_h = self._bessel_function_interped(abs(d_inner_h)) + log_l = d_inner_h - optimal_snr_squared / 2 else: - log_l = matched_filter_snr_squared.real - optimal_snr_squared / 2 + log_l = d_inner_h.real - optimal_snr_squared / 2 return log_l.real - def _setup_rho(self, matched_filter_snr_squared, optimal_snr_squared): + def _setup_rho(self, d_inner_h, optimal_snr_squared): rho_opt_ref = (optimal_snr_squared.real * self.parameters['luminosity_distance'] ** 2 / self._ref_dist ** 2.) - rho_mf_ref = matched_filter_snr_squared * \ - self.parameters['luminosity_distance'] / self._ref_dist + rho_mf_ref = (d_inner_h * self.parameters['luminosity_distance'] / + self._ref_dist) return rho_mf_ref, rho_opt_ref def log_likelihood(self): @@ -262,12 +262,12 @@ class GravitationalWaveTransient(likelihood.Likelihood): for ii, rho_opt_ref in enumerate(self._rho_opt_ref_array): for jj, rho_mf_ref in enumerate(self._rho_mf_ref_array): optimal_snr_squared_array = rho_opt_ref * self._ref_dist ** 2. / self._distance_array ** 2 - matched_filter_snr_squared_array = rho_mf_ref * self._ref_dist / self._distance_array + d_inner_h_array = rho_mf_ref * self._ref_dist / self._distance_array if self.phase_marginalization: - matched_filter_snr_squared_array =\ - self._bessel_function_interped(abs(matched_filter_snr_squared_array)) + d_inner_h_array =\ + self._bessel_function_interped(abs(d_inner_h_array)) self._dist_margd_loglikelihood_array[ii][jj] = \ - logsumexp(matched_filter_snr_squared_array - optimal_snr_squared_array / 2, + logsumexp(d_inner_h_array - optimal_snr_squared_array / 2, b=self.distance_prior_array * self._delta_distance) log_norm = logsumexp(0. / self._distance_array - 0. / self._distance_array ** 2., b=self.distance_prior_array * self._delta_distance)