From e2ab870a0d178d26f7c2399a47607d8d92a35d58 Mon Sep 17 00:00:00 2001
From: RorySmith <rory.smith@caltech.edu>
Date: Thu, 31 May 2018 12:02:47 +1000
Subject: [PATCH] stable version of dist marg, and time-and-phase marg. Need to
fix dist and time marg before moving onto dist and tiem nad phase marg
---
examples/injection_examples/basic_tutorial.py | 4 +-
tupak/likelihood.py | 72 ++++++++++---------
2 files changed, 41 insertions(+), 35 deletions(-)
diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py
index 7fae64df1..d0fdb9141 100644
--- a/examples/injection_examples/basic_tutorial.py
+++ b/examples/injection_examples/basic_tutorial.py
@@ -19,13 +19,13 @@ label = 'basic_tutorial'
tupak.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
-#np.random.seed(170)
+np.random.seed(2157034243)
# We are going to inject a binary black hole waveform. We first establish a dictionary of parameters that
# includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2),
# spins of both black holes (a, tilt, phi), etc.
injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3,
- luminosity_distance=4000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
+ luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
# Create the waveform_generator using a LAL BinaryBlackHole source function
diff --git a/tupak/likelihood.py b/tupak/likelihood.py
index fea574fd2..e3bd07b50 100644
--- a/tupak/likelihood.py
+++ b/tupak/likelihood.py
@@ -85,7 +85,7 @@ class GravitationalWaveTransient(Likelihood):
if self.phase_marginalization:
self.bessel_function_interped = None
self.setup_phase_marginalization()
- prior['psi'] = 0
+ prior['phase'] = 0
@property
def prior(self):
@@ -141,25 +141,24 @@ class GravitationalWaveTransient(Likelihood):
if self.time_marginalization:
delta_tc = self.waveform_generator.sampling_frequency
- #print (1./delta_tc, (interferometer.duration) )
- if self.distance_marginalization:
- rho_opt_1 = optimal_snr_squared.real * \
- self.waveform_generator.parameters['luminosity_distance'] ** 2
- rho_mf_1_tc_array = matched_filter_snr_squared_tc_array.real * \
- self.waveform_generator.parameters['luminosity_distance']
+ if self.distance_marginalization:
- ## TRY time marg before dist marg:
- #time_marged_log_l_dist_array = logsumexp(rho_mf_1_tc_array, b=delta_tc)
- #log_l = self.interp_dist_margd_likelihood(rho_mf_1_tc_array, rho_opt_1)[0]
+ rho_opt_ref = optimal_snr_squared.real * \
+ self.waveform_generator.parameters['luminosity_distance'] ** 2 \
+ / self.ref_dist**2.
+ rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array.real * \
+ self.waveform_generator.parameters['luminosity_distance'] \
+ / self.ref_dist
- dist_marged_log_l_tc_array = self.interp_dist_margd_likelihood(rho_mf_1_tc_array, rho_opt_1)[0]
+ dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array, rho_opt_ref)[0]
log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc)
+
elif self.phase_marginalization:
- phase_marged_log_l_tc_array = self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array))
- log_l = logsumexp(phase_marged_log_l_tc_array, axis=0, b=delta_tc)
+ log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
+ - optimal_snr_squared/2.
else:
@@ -168,12 +167,14 @@ class GravitationalWaveTransient(Likelihood):
elif self.distance_marginalization:
#TODO: order should be: (phase-> dist-> tc)
- 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']
+ rho_opt_ref = optimal_snr_squared.real * \
+ 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
- log_l=self.interp_dist_margd_likelihood(rho_mf_1, rho_opt_1)[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))
@@ -198,48 +199,53 @@ 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.dist_margd_likelihood_array = np.zeros((1000,1000))
+ self.rho_opt_ref_array = np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[0]) # optimal filter snr at fiducial distance of 1 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 1 Mpc
- self.rho_opt_1_array = np.linspace(5e5,7e8,self.dist_margd_likelihood_array.shape[0]) # optimal filter snr at fiducial distance of 1 Mpc
- self.rho_mf_1_array = np.linspace(-45000,2e6,self.dist_margd_likelihood_array.shape[1]) # matched filter snr at fiducial distance of 1 Mpc
+ for ii, rho_opt_ref in enumerate(self.rho_opt_ref_array):
- for ii, rho_opt_1 in enumerate(self.rho_opt_1_array):
+ for jj, rho_mf_ref in enumerate(self.rho_mf_ref_array):
- for jj, rho_mf_1 in enumerate(self.rho_mf_1_array):
-
- optimal_snr_squared_array = rho_opt_1 \
+ optimal_snr_squared_array = rho_opt_ref * self.ref_dist**2. \
/ self.distance_array ** 2
- matched_filter_snr_squared_array = rho_mf_1 \
+ matched_filter_snr_squared_array = rho_mf_ref * self.ref_dist \
/ self.distance_array
- self.dist_margd_likelihood_array[ii][jj] = \
+ 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)
log_norm = logsumexp(0 / self.distance_array - 0/self.distance_array**2. / 2, b=self.distance_prior_array * self.delta_distance)
- self.dist_margd_likelihood_array -= log_norm
+ self.dist_margd_loglikelihood_array -= log_norm
- self.interp_dist_margd_likelihood = interp2d(self.rho_mf_1_array, self.rho_opt_1_array, self.dist_margd_likelihood_array)
+ self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array, self.dist_margd_loglikelihood_array)
- #np.save("dist_margd_likelihood_array",self.dist_margd_likelihood_array)
+ #np.save("dist_margd_loglikelihood_array",self.dist_margd_loglikelihood_array)
#sys.exit()
+ '''def setupt_time_marginalization(self):
+ if 'geocent_time' not in self.prior.keys() or not isinstance(self.prior['geocent_time'], tupak.prior.Prior):
+ logging.info('No prior provided for geocent_time, using default prior.')
+ self.prior['geocent_time'] ='''
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.')
- self.prior['psi'] = tupak.prior.create_default_prior('psi')
+ if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
+ logging.info('No prior provided for phase at coalescence, using default prior.')
+ self.prior['phase'] = tupak.prior.create_default_prior('phase')
self.bessel_function_interped = interp1d(np.linspace(0, 1e6, int(1e5)),
np.log([i0e(snr) for snr in np.linspace(0, 1e6, int(1e5))])
+ np.linspace(0, 1e6, int(1e5)),
bounds_error=False, fill_value=-np.inf)
+
class BasicGravitationalWaveTransient(Likelihood):
""" A basic gravitational wave transient likelihood
--
GitLab