Skip to content
Snippets Groups Projects
Commit 6faed1e1 authored by Moritz Huebner's avatar Moritz Huebner
Browse files

Removed some code redundancy for the time marginalisation

parent d26e9a0f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 46 additions and 53 deletions
tupak/likelihood.py
+ 46
− 53
View file @ 6faed1e1
from __future__ import division, print_function from __future__ import division, print_function
import numpy as np import numpy as np
try: try:
...@@ -7,7 +8,6 @@ except ImportError: ...@@ -7,7 +8,6 @@ except ImportError:
from scipy.misc import logsumexp from scipy.misc import logsumexp
from scipy.special import i0e from scipy.special import i0e
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
from scipy.integrate import quad as gaussian_quadrature
from scipy.interpolate import interp2d from scipy.interpolate import interp2d
import tupak import tupak
import logging import logging
...@@ -60,6 +60,7 @@ class GravitationalWaveTransient(Likelihood): ...@@ -60,6 +60,7 @@ class GravitationalWaveTransient(Likelihood):
some model parameters some model parameters
""" """
def __init__(self, interferometers, waveform_generator, time_marginalization=False, distance_marginalization=False, def __init__(self, interferometers, waveform_generator, time_marginalization=False, distance_marginalization=False,
phase_marginalization=False, prior=None): phase_marginalization=False, prior=None):
...@@ -77,7 +78,7 @@ class GravitationalWaveTransient(Likelihood): ...@@ -77,7 +78,7 @@ class GravitationalWaveTransient(Likelihood):
self.delta_distance = 0 self.delta_distance = 0
self.distance_prior_array = np.array([]) self.distance_prior_array = np.array([])
self.setup_distance_marginalization() 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: if self.phase_marginalization:
self.bessel_function_interped = None self.bessel_function_interped = None
...@@ -111,9 +112,10 @@ class GravitationalWaveTransient(Likelihood): ...@@ -111,9 +112,10 @@ class GravitationalWaveTransient(Likelihood):
matched_filter_snr_squared = 0 matched_filter_snr_squared = 0
optimal_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: for interferometer in self.interferometers:
duration = interferometer.duration
signal_ifo = interferometer.get_detector_response(waveform_polarizations, signal_ifo = interferometer.get_detector_response(waveform_polarizations,
self.waveform_generator.parameters) self.waveform_generator.parameters)
matched_filter_snr_squared += tupak.utils.matched_filter_snr_squared( matched_filter_snr_squared += tupak.utils.matched_filter_snr_squared(
...@@ -121,69 +123,61 @@ class GravitationalWaveTransient(Likelihood): ...@@ -121,69 +123,61 @@ class GravitationalWaveTransient(Likelihood):
optimal_snr_squared += tupak.utils.optimal_snr_squared( optimal_snr_squared += tupak.utils.optimal_snr_squared(
signal_ifo, interferometer, self.waveform_generator.time_duration) signal_ifo, interferometer, self.waveform_generator.time_duration)
if self.time_marginalization: if self.time_marginalization:
interferometer.time_marginalization = self.time_marginalization interferometer.time_marginalization = self.time_marginalization
matched_filter_snr_squared_tc_array += 4. * (1. / interferometer.duration) * np.fft.ifft(
if matched_filter_snr_squared_tc_array is None: 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])
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])
if self.time_marginalization: if self.time_marginalization:
delta_tc = 1./self.waveform_generator.sampling_frequency delta_tc = 1. / self.waveform_generator.sampling_frequency
tc_log_norm = np.log(interferometer.duration*delta_tc) tc_log_norm = np.log(duration * delta_tc)
if self.distance_marginalization: if self.distance_marginalization:
rho_opt_ref = optimal_snr_squared.real * \ rho_opt_ref = optimal_snr_squared.real * \
self.waveform_generator.parameters['luminosity_distance'] ** 2 \ self.waveform_generator.parameters['luminosity_distance'] ** 2 \
/ self.ref_dist**2. / self.ref_dist ** 2.
rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \ rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \
self.waveform_generator.parameters['luminosity_distance'] \ self.waveform_generator.parameters['luminosity_distance'] \
/ self.ref_dist / self.ref_dist
if self.phase_marginalization: if self.phase_marginalization:
phase_margd_rho_mf_tc_array = self.bessel_function_interped(abs(rho_mf_ref_tc_array)) 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 log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - tc_log_norm
else: 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: elif self.phase_marginalization:
log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\ log_l = logsumexp(self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc) \
- optimal_snr_squared/2. - tc_log_norm - optimal_snr_squared / 2. - tc_log_norm
else: 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: elif self.distance_marginalization:
rho_opt_ref = optimal_snr_squared.real * \ rho_opt_ref = optimal_snr_squared.real * \
self.waveform_generator.parameters['luminosity_distance'] ** 2 \ self.waveform_generator.parameters['luminosity_distance'] ** 2 \
/ self.ref_dist**2. / self.ref_dist ** 2.
rho_mf_ref = matched_filter_snr_squared.real * \ rho_mf_ref = matched_filter_snr_squared.real * \
self.waveform_generator.parameters['luminosity_distance'] \ self.waveform_generator.parameters['luminosity_distance'] \
/ self.ref_dist / 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: elif self.phase_marginalization:
matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared)) matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared))
...@@ -208,36 +202,35 @@ class GravitationalWaveTransient(Likelihood): ...@@ -208,36 +202,35 @@ class GravitationalWaveTransient(Likelihood):
for distance in self.distance_array]) for distance in self.distance_array])
### Make the lookup table ### ### Make the lookup table ###
self.ref_dist = 1000 # 1000 Mpc self.ref_dist = 1000 # 1000 Mpc
self.dist_margd_loglikelihood_array = np.zeros((400,800)) 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.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 ii, rho_opt_ref in enumerate(self.rho_opt_ref_array):
for jj, rho_mf_ref in enumerate(self.rho_mf_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 / self.distance_array ** 2
matched_filter_snr_squared_array = rho_mf_ref * self.ref_dist \ matched_filter_snr_squared_array = rho_mf_ref * self.ref_dist \
/ self.distance_array / self.distance_array
self.dist_margd_loglikelihood_array[ii][jj] = \ self.dist_margd_loglikelihood_array[ii][jj] = \
logsumexp(matched_filter_snr_squared_array - \ logsumexp(matched_filter_snr_squared_array - \
optimal_snr_squared_array / 2, \ optimal_snr_squared_array / 2, \
b= self.distance_prior_array * self.delta_distance) b=self.distance_prior_array * self.delta_distance)
log_norm = logsumexp(0. / self.distance_array - 0./self.distance_array**2.,\ log_norm = logsumexp(0. / self.distance_array - 0. / self.distance_array ** 2., \
b=self.distance_prior_array* self.delta_distance) b=self.distance_prior_array * self.delta_distance)
self.dist_margd_loglikelihood_array -= log_norm self.dist_margd_loglikelihood_array -= log_norm
self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array,
self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array, self.dist_margd_loglikelihood_array) self.dist_margd_loglikelihood_array)
def setup_phase_marginalization(self): def setup_phase_marginalization(self):
if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior): if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
...@@ -249,7 +242,6 @@ class GravitationalWaveTransient(Likelihood): ...@@ -249,7 +242,6 @@ class GravitationalWaveTransient(Likelihood):
bounds_error=False, fill_value=-np.inf) bounds_error=False, fill_value=-np.inf)
class BasicGravitationalWaveTransient(Likelihood): class BasicGravitationalWaveTransient(Likelihood):
""" A basic gravitational wave transient likelihood """ A basic gravitational wave transient likelihood
...@@ -272,6 +264,7 @@ class BasicGravitationalWaveTransient(Likelihood): ...@@ -272,6 +264,7 @@ class BasicGravitationalWaveTransient(Likelihood):
some model parameters some model parameters
""" """
def __init__(self, interferometers, waveform_generator): def __init__(self, interferometers, waveform_generator):
Likelihood.__init__(self, waveform_generator.parameters) Likelihood.__init__(self, waveform_generator.parameters)
self.interferometers = interferometers self.interferometers = interferometers
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment