From 61fb7ba8f41d7d18ad379a27531244285aa21e59 Mon Sep 17 00:00:00 2001
From: RorySmith <rory.smith@caltech.edu>
Date: Fri, 1 Jun 2018 13:51:27 +1000
Subject: [PATCH] working (dist,time,phase) marg
---
examples/injection_examples/basic_tutorial.py | 4 +-
tupak/detector.py | 14 ++++--
tupak/likelihood.py | 45 ++++++++++---------
tupak/prior.py | 6 +--
4 files changed, 39 insertions(+), 30 deletions(-)
diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py
index d0fdb9141..8a1be60e3 100644
--- a/examples/injection_examples/basic_tutorial.py
+++ b/examples/injection_examples/basic_tutorial.py
@@ -19,7 +19,7 @@ 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(2157034243)
+np.random.seed(88170235)
# 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),
@@ -55,7 +55,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=False, distance_marginalization=True, prior=priors)
+likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=False, phase_marginalization=False, distance_marginalization=False, prior=priors)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
diff --git a/tupak/detector.py b/tupak/detector.py
index f38f3aa34..6030523b8 100644
--- a/tupak/detector.py
+++ b/tupak/detector.py
@@ -69,6 +69,7 @@ class Interferometer(object):
self.frequency_array = []
self.sampling_frequency = None
self.duration = None
+ self.time_marginalization = False
@property
def minimum_frequency(self):
@@ -226,7 +227,7 @@ class Interferometer(object):
det_response = self.antenna_response(
parameters['ra'],
parameters['dec'],
- parameters['geocent_time'],
+ self.epoch,#parameters['geocent_time'],
parameters['psi'], mode)
signal[mode] = waveform_polarizations[mode] * det_response
@@ -237,9 +238,14 @@ class Interferometer(object):
time_shift = self.time_delay_from_geocenter(
parameters['ra'],
parameters['dec'],
- parameters['geocent_time'])
+ self.epoch)#parameters['geocent_time'])
+
+ if self.time_marginalization:
+ dt = time_shift # when marginalizing over time we only care about relative time shifts between detectors and marginalized over
+ # all candidate coalescence times
+ else:
+ dt = self.epoch - (parameters['geocent_time'] - time_shift)
- dt = self.epoch - (parameters['geocent_time'] - time_shift)
signal_ifo = signal_ifo * np.exp(
-1j * 2 * np.pi * dt * self.frequency_array)
@@ -688,7 +694,7 @@ def get_interferometer_with_fake_noise_and_injection(
interferometer = get_empty_interferometer(name)
interferometer.set_data(
sampling_frequency=sampling_frequency, duration=time_duration,
- from_power_spectral_density=True)
+ from_power_spectral_density=True, epoch=(injection_parameters['geocent_time']+2)-time_duration)
interferometer.inject_signal(
waveform_polarizations=injection_polarizations,
parameters=injection_parameters)
diff --git a/tupak/likelihood.py b/tupak/likelihood.py
index e3bd07b50..c983bb57d 100644
--- a/tupak/likelihood.py
+++ b/tupak/likelihood.py
@@ -72,9 +72,6 @@ class GravitationalWaveTransient(Likelihood):
self.phase_marginalization = phase_marginalization
self.prior = prior
- #if self.time_marginalization:
- # prior['geocent_time'] = 0 #FIXME: should be fixed to injection value; probably easiest to pass in a delta function prior
-
if self.distance_marginalization:
self.distance_array = np.array([])
self.delta_distance = 0
@@ -127,45 +124,57 @@ class GravitationalWaveTransient(Likelihood):
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]).real * len(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] ).real * len(interferometer.data[0:-1])
+ / interferometer.power_spectral_density_array[0:-1] ) * len(interferometer.data[0:-1])
if self.time_marginalization:
- delta_tc = self.waveform_generator.sampling_frequency
+ delta_tc = 1./self.waveform_generator.sampling_frequency
+ tc_log_norm = np.log(interferometer.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.
- rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array.real * \
+ rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \
self.waveform_generator.parameters['luminosity_distance'] \
/ self.ref_dist
- 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)
+ 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)
+ 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)
+
+ log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - dist_margd_loglikelihood_array
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.
+ - optimal_snr_squared/2. - tc_log_norm
else:
- log_l = logsumexp(matched_filter_snr_squared_tc_array, axis=0, b=delta_tc) - optimal_snr_squared / 2
+ 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:
- #TODO: order should be: (phase-> dist-> tc)
rho_opt_ref = optimal_snr_squared.real * \
self.waveform_generator.parameters['luminosity_distance'] ** 2 \
@@ -219,21 +228,15 @@ class GravitationalWaveTransient(Likelihood):
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)
+ 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)
+ 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)
- #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 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
diff --git a/tupak/prior.py b/tupak/prior.py
index 622dabad7..78e1f4286 100644
--- a/tupak/prior.py
+++ b/tupak/prior.py
@@ -418,9 +418,9 @@ def create_default_prior(name):
Default prior distribution for that parameter, if unknown None is returned.
"""
default_priors = {
- 'mass_1': Uniform(name=name, minimum=5, maximum=100),
- 'mass_2': Uniform(name=name, minimum=5, maximum=100),
- 'chirp_mass': Uniform(name=name, minimum=5, maximum=100),
+ 'mass_1': Uniform(name=name, minimum=20, maximum=100),
+ 'mass_2': Uniform(name=name, minimum=20, maximum=100),
+ 'chirp_mass': Uniform(name=name, minimum=25, maximum=100),
'total_mass': Uniform(name=name, minimum=10, maximum=200),
'mass_ratio': Uniform(name=name, minimum=0.125, maximum=1),
'symmetric_mass_ratio': Uniform(name=name, minimum=8 / 81, maximum=0.25),
--
GitLab