diff --git a/examples/open_data_examples/GW150914.py b/examples/open_data_examples/GW150914.py
index 56ea1cde47d0274b68c2088ccd79ca571937d5a2..fbc858df47752084785d1e491efa8856be2c3c9a 100644
--- a/examples/open_data_examples/GW150914.py
+++ b/examples/open_data_examples/GW150914.py
@@ -37,6 +37,7 @@ prior['mass_1'] = tupak.prior.Uniform(30, 50, 'mass_1')
prior['mass_2'] = tupak.prior.Uniform(20, 40, 'mass_2')
prior['geocent_time'] = tupak.prior.Uniform(
time_of_event - 0.1, time_of_event + 0.1, name='geocent_time')
+#prior['geocent_time'] = time_of_event
prior['luminosity_distance'] = tupak.prior.PowerLaw(
alpha=2, minimum=100, maximum=1000)
@@ -53,7 +54,7 @@ waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=in
# In this step, we define the likelihood. Here we use the standard likelihood
# function, passing it the data and the waveform generator.
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
+likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers, waveform_generator, time_marginalization=False, distance_marginalization=True, prior=prior)
# Finally, we run the sampler. This function takes the likelihood and prio
# along with some options for how to do the sampling and how to save the data
diff --git a/tupak/likelihood.py b/tupak/likelihood.py
index bde687537a157609b86d65f524d72fe36db8613e..045a4d5cea0cf885cb6c4879e7c8996b503f4f0a 100644
--- a/tupak/likelihood.py
+++ b/tupak/likelihood.py
@@ -7,6 +7,8 @@ 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
@@ -72,19 +74,19 @@ class GravitationalWaveTransient(Likelihood):
if self.time_marginalization:
self.setup_time_marginalization()
- self.prior['geocent_time'] = 1
+ prior['geocent_time'] = 0 #FIXME: should be fixed to injection value
if self.distance_marginalization:
self.distance_array = np.array([])
self.delta_distance = 0
self.distance_prior_array = np.array([])
self.setup_distance_marginalization()
- self.prior['luminosity_distance'] = 1
+ 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
self.setup_phase_marginalization()
- self.prior['psi'] = 0
+ prior['psi'] = 0
@property
def prior(self):
@@ -137,7 +139,6 @@ class GravitationalWaveTransient(Likelihood):
if self.time_marginalization:
-
delta_tc = self.waveform_generator.sampling_frequency
#print (1./delta_tc, (interferometer.duration) )
if self.distance_marginalization: #FIXME: This doens't implement the lookup-table method
@@ -159,7 +160,7 @@ class GravitationalWaveTransient(Likelihood):
log_l = logsumexp(matched_filter_snr_squared_tc_array, axis=0, b=delta_tc) - optimal_snr_squared / 2
elif self.distance_marginalization: #FIXME: This doens't implement the lookup-table method
-
+ '''
optimal_snr_squared_array = optimal_snr_squared \
* self.waveform_generator.parameters['luminosity_distance'] ** 2 \
/ self.distance_array ** 2
@@ -174,6 +175,21 @@ class GravitationalWaveTransient(Likelihood):
log_l = logsumexp(matched_filter_snr_squared_array - optimal_snr_squared_array / 2,
b=self.distance_prior_array * self.delta_distance)
+ '''
+ rho_opt_1 = optimal_snr_squared.real * \
+ self.waveform_generator.parameters['luminosity_distance'] ** 2
+ rho_mf_1 = matched_filter_snr_squared.real * \
+ self.waveform_generator.parameters['luminosity_distance']
+
+ #nearest_rho_opt_1_idx = np.argmin(abs(self.rho_opt_1_array-rho_opt_1))
+ #nearest_rho_mf_1_idx = np.argmin(abs(self.rho_mf_1_array-rho_mf_1))
+
+ #log_l = self.dist_margd_likelihood_array[nearest_rho_opt_1_idx][nearest_rho_mf_1_idx]
+ #print("SDG", self.interp_dist_margd_likelihood(rho_opt_1, rho_mf_1))
+ #print (self.interp_dist_margd_likelihood(rho_opt_1, rho_mf_1)[0], self.interp_dist_margd_likelihood(rho_opt_1, rho_mf_1))
+
+ log_l=self.interp_dist_margd_likelihood(rho_mf_1, rho_opt_1)[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
@@ -201,6 +217,51 @@ class GravitationalWaveTransient(Likelihood):
self.distance_prior_array = np.array([self.prior['luminosity_distance'].prob(distance)
for distance in self.distance_array])
+ ### Make the lookup table ###
+
+ self.dist_margd_likelihood_array = np.zeros((1000,1000))
+
+ self.rho_opt_1_array = np.linspace(1e-10,100000000,self.dist_margd_likelihood_array.shape[0]) # optimal filter snr at fiducial distance of 1 Mpc
+ self.rho_mf_1_array = np.linspace(-100000000,100000000,self.dist_margd_likelihood_array.shape[1]) # matched filter snr at fiducial distance of 1 Mpc
+ #print(min(self.rho_opt_1_array)/max(self.distance_array ** 2), max(self.rho_opt_1_array)/min(self.distance_array ** 2),\
+ #cmin(self.rho_mf_1_array)/max(self.distance_array), max(self.rho_mf_1_array)/min(self.distance_array))
+
+ integrand = lambda D, a: np.exp(a[0]/D - a[1]/D**2. / 2. - a[2]) *\
+ self.prior['luminosity_distance'].prob(D)
+
+ for ii, rho_opt_1 in enumerate(self.rho_opt_1_array):
+
+ for jj, rho_mf_1 in enumerate(self.rho_mf_1_array):
+
+ optimal_snr_squared_array = rho_opt_1 \
+ / self.distance_array ** 2
+
+ matched_filter_snr_squared_array = rho_mf_1 \
+ / self.distance_array
+
+ #regularizer = np.min(unmarged_log_l)'''
+ self.dist_margd_likelihood_array[ii][jj] = \
+ logsumexp(matched_filter_snr_squared_array - \
+ optimal_snr_squared_array / 2, \
+ b=self.distance_prior_array * self.delta_distance)
+
+ #reg = np.max(rho_mf_1/(self.distance_array) - rho_opt_1/(self.distance_array)**2. / 2.)
+
+ #self.dist_margd_likelihood_array[ii][jj] = #np.log( gaussian_quadrature(integrand,
+ #self.prior['luminosity_distance'].minimum,self.prior['luminosity_distance'].maximum,\
+ #args=[rho_mf_1,rho_opt_1,reg])[0] ) + reg
+
+ #print(ii, self.dist_margd_likelihood_array[ii])
+ self.interp_dist_margd_likelihood = interp2d(self.rho_mf_1_array, self.rho_opt_1_array, self.dist_margd_likelihood_array)
+
+ #np.save("dist_margd_likelihood_array",self.dist_margd_likelihood_array)
+ #sys.exit()
+ '''self.dist_margd_likelihood_array[ii][jj] = \
+ logsumexp(matched_filter_snr_squared_array - \
+ optimal_snr_squared_array / 2, \
+ b=self.distance_prior_array * self.delta_distance)'''
+
+
def setup_phase_marginalization(self):
if 'psi' not in self.prior.keys() or not isinstance(self.prior['psi'], tupak.prior.Prior):
logging.info('No prior provided for polarization, using default prior.')