From 6faed1e16c3932e2c284333988faceabb962c298 Mon Sep 17 00:00:00 2001
From: Moritz Huebner <email@moritz-huebner.de>
Date: Mon, 4 Jun 2018 19:11:37 +1000
Subject: [PATCH] Removed some code redundancy for the time marginalisation
---
tupak/likelihood.py | 99 +++++++++++++++++++++------------------------
1 file changed, 46 insertions(+), 53 deletions(-)
diff --git a/tupak/likelihood.py b/tupak/likelihood.py
index 1e7922305..479c77471 100644
--- a/tupak/likelihood.py
+++ b/tupak/likelihood.py
@@ -1,4 +1,5 @@
from __future__ import division, print_function
+
import numpy as np
try:
@@ -7,7 +8,6 @@ except ImportError:
from scipy.misc import logsumexp
from scipy.special import i0e
from scipy.interpolate import interp1d
-from scipy.integrate import quad as gaussian_quadrature
from scipy.interpolate import interp2d
import tupak
import logging
@@ -60,6 +60,7 @@ class GravitationalWaveTransient(Likelihood):
some model parameters
"""
+
def __init__(self, interferometers, waveform_generator, time_marginalization=False, distance_marginalization=False,
phase_marginalization=False, prior=None):
@@ -77,7 +78,7 @@ class GravitationalWaveTransient(Likelihood):
self.delta_distance = 0
self.distance_prior_array = np.array([])
self.setup_distance_marginalization()
- prior['luminosity_distance'] = 1 # this means the prior is a delta function fixed at the RHS value
+ prior['luminosity_distance'] = 1 # this means the prior is a delta function fixed at the RHS value
if self.phase_marginalization:
self.bessel_function_interped = None
@@ -111,9 +112,10 @@ class GravitationalWaveTransient(Likelihood):
matched_filter_snr_squared = 0
optimal_snr_squared = 0
- matched_filter_snr_squared_tc_array = None
-
+ matched_filter_snr_squared_tc_array = np.zeros(self.interferometers[0].data.shape)
+ duration = 0
for interferometer in self.interferometers:
+ duration = interferometer.duration
signal_ifo = interferometer.get_detector_response(waveform_polarizations,
self.waveform_generator.parameters)
matched_filter_snr_squared += tupak.utils.matched_filter_snr_squared(
@@ -121,69 +123,61 @@ class GravitationalWaveTransient(Likelihood):
optimal_snr_squared += tupak.utils.optimal_snr_squared(
signal_ifo, interferometer, self.waveform_generator.time_duration)
-
if self.time_marginalization:
-
interferometer.time_marginalization = self.time_marginalization
-
- if matched_filter_snr_squared_tc_array is None:
-
- matched_filter_snr_squared_tc_array = 4. * (1./interferometer.duration) \
- * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
- / interferometer.power_spectral_density_array[0:-1]) * len(interferometer.data[0:-1])
- else:
-
- matched_filter_snr_squared_tc_array += 4. * (1./interferometer.duration) \
- * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
- / interferometer.power_spectral_density_array[0:-1] ) * len(interferometer.data[0:-1])
+ matched_filter_snr_squared_tc_array += 4. * (1. / interferometer.duration) * np.fft.ifft(
+ signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1]
+ / interferometer.power_spectral_density_array[0:-1]) * len(interferometer.data[0:-1])
if self.time_marginalization:
- delta_tc = 1./self.waveform_generator.sampling_frequency
- tc_log_norm = np.log(interferometer.duration*delta_tc)
+ delta_tc = 1. / self.waveform_generator.sampling_frequency
+ tc_log_norm = np.log(duration * delta_tc)
if self.distance_marginalization:
rho_opt_ref = optimal_snr_squared.real * \
- self.waveform_generator.parameters['luminosity_distance'] ** 2 \
- / self.ref_dist**2.
+ self.waveform_generator.parameters['luminosity_distance'] ** 2 \
+ / self.ref_dist ** 2.
rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \
- self.waveform_generator.parameters['luminosity_distance'] \
- / self.ref_dist
-
+ self.waveform_generator.parameters['luminosity_distance'] \
+ / self.ref_dist
if self.phase_marginalization:
phase_margd_rho_mf_tc_array = self.bessel_function_interped(abs(rho_mf_ref_tc_array))
- dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(phase_margd_rho_mf_tc_array, rho_opt_ref)
+ dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(phase_margd_rho_mf_tc_array,
+ rho_opt_ref)
log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - tc_log_norm
else:
- dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array.real, rho_opt_ref)
+ dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array.real,
+ rho_opt_ref)
- log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc)
+ log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc)
elif self.phase_marginalization:
- log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
- - optimal_snr_squared/2. - tc_log_norm
+ log_l = logsumexp(self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc) \
+ - optimal_snr_squared / 2. - tc_log_norm
else:
- log_l = logsumexp(matched_filter_snr_squared_tc_array.real, axis=0, b=delta_tc) - optimal_snr_squared / 2. - tc_log_norm
+ log_l = logsumexp(matched_filter_snr_squared_tc_array.real, axis=0,
+ b=delta_tc) - optimal_snr_squared / 2. - tc_log_norm
elif self.distance_marginalization:
rho_opt_ref = optimal_snr_squared.real * \
- self.waveform_generator.parameters['luminosity_distance'] ** 2 \
- / self.ref_dist**2.
+ self.waveform_generator.parameters['luminosity_distance'] ** 2 \
+ / self.ref_dist ** 2.
rho_mf_ref = matched_filter_snr_squared.real * \
- self.waveform_generator.parameters['luminosity_distance'] \
- / self.ref_dist
+ self.waveform_generator.parameters['luminosity_distance'] \
+ / self.ref_dist
- log_l=self.interp_dist_margd_loglikelihood(rho_mf_ref, rho_opt_ref)[0]
+ log_l = self.interp_dist_margd_loglikelihood(rho_mf_ref, rho_opt_ref)[0]
elif self.phase_marginalization:
matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared))
@@ -208,36 +202,35 @@ class GravitationalWaveTransient(Likelihood):
for distance in self.distance_array])
### Make the lookup table ###
- self.ref_dist = 1000 # 1000 Mpc
- self.dist_margd_loglikelihood_array = np.zeros((400,800))
-
- self.rho_opt_ref_array = np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[0]) # optimal filter snr at fiducial distance of ref_dist Mpc
- self.rho_mf_ref_array = np.hstack( (-np.logspace(2,-3,self.dist_margd_loglikelihood_array.shape[1]/2), \
- np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[1]/2)) ) # matched filter snr at fiducial distance of ref_dist Mpc
+ self.ref_dist = 1000 # 1000 Mpc
+ self.dist_margd_loglikelihood_array = np.zeros((400, 800))
+ self.rho_opt_ref_array = np.logspace(-3, 4, self.dist_margd_loglikelihood_array.shape[
+ 0]) # optimal filter snr at fiducial distance of ref_dist Mpc
+ self.rho_mf_ref_array = np.hstack((-np.logspace(2, -3, self.dist_margd_loglikelihood_array.shape[1] / 2), \
+ np.logspace(-3, 4, self.dist_margd_loglikelihood_array.shape[
+ 1] / 2))) # matched filter snr at fiducial distance of ref_dist Mpc
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. \
+ 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
+ / self.distance_array
self.dist_margd_loglikelihood_array[ii][jj] = \
- logsumexp(matched_filter_snr_squared_array - \
- optimal_snr_squared_array / 2, \
- b= self.distance_prior_array * self.delta_distance)
+ logsumexp(matched_filter_snr_squared_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)
+ log_norm = logsumexp(0. / self.distance_array - 0. / self.distance_array ** 2., \
+ b=self.distance_prior_array * self.delta_distance)
self.dist_margd_loglikelihood_array -= log_norm
-
- self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array, self.dist_margd_loglikelihood_array)
-
+ self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array,
+ self.dist_margd_loglikelihood_array)
def setup_phase_marginalization(self):
if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
@@ -249,7 +242,6 @@ class GravitationalWaveTransient(Likelihood):
bounds_error=False, fill_value=-np.inf)
-
class BasicGravitationalWaveTransient(Likelihood):
""" A basic gravitational wave transient likelihood
@@ -272,6 +264,7 @@ class BasicGravitationalWaveTransient(Likelihood):
some model parameters
"""
+
def __init__(self, interferometers, waveform_generator):
Likelihood.__init__(self, waveform_generator.parameters)
self.interferometers = interferometers
--
GitLab