Commit f3afa08b authored by Colm Talbot's avatar Colm Talbot
Browse files

make marginalised likelihood (with no marginalisation) default

parent 7ef8da5b
......@@ -4,52 +4,20 @@ try:
from scipy.special import logsumexp
except ImportError:
from scipy.misc import logsumexp
from scipy.special import i0, i0e
from scipy.special import i0e
from scipy.interpolate import interp1d
import tupak
import logging
class Likelihood(object):
def __init__(self, interferometers, waveform_generator):
def __init__(self, interferometers, waveform_generator, distance_marginalization=False, phase_marginalization=False,
prior=None):
# Likelihood.__init__(self, interferometers, waveform_generator)
self.interferometers = interferometers
self.waveform_generator = waveform_generator
def noise_log_likelihood(self):
log_l = 0
for interferometer in self.interferometers:
log_l -= 2. / self.waveform_generator.time_duration * np.sum(
abs(interferometer.data) ** 2 / interferometer.power_spectral_density_array)
return log_l.real
def log_likelihood(self):
log_l = 0
waveform_polarizations = self.waveform_generator.frequency_domain_strain()
if waveform_polarizations is None:
return np.nan_to_num(-np.inf)
for interferometer in self.interferometers:
log_l += self.log_likelihood_interferometer(waveform_polarizations, interferometer)
return log_l.real
def log_likelihood_interferometer(self, waveform_polarizations, interferometer):
signal_ifo = interferometer.get_detector_response(waveform_polarizations, self.waveform_generator.parameters)
log_l = - 2. / self.waveform_generator.time_duration * np.vdot(interferometer.data - signal_ifo,
(interferometer.data - signal_ifo)
/ interferometer.power_spectral_density_array)
return log_l.real
def log_likelihood_ratio(self):
return self.log_likelihood() - self.noise_log_likelihood()
class MarginalizedLikelihood(Likelihood):
def __init__(self, interferometers, waveform_generator, distance_marginalization=False, phase_marginalization=False,
time_marginalization=False, prior=None):
Likelihood.__init__(self, interferometers, waveform_generator)
self.distance_marginalization = distance_marginalization
self.phase_marginalization = phase_marginalization
self.time_marginalization = time_marginalization
self.prior = prior
if self.distance_marginalization:
......@@ -64,29 +32,6 @@ class MarginalizedLikelihood(Likelihood):
self.setup_phase_marginalization()
prior['psi'] = 0
if self.time_marginalization:
logging.warning('Time marginalisation not yet implemented.')
self.time_marginalization = False
def setup_distance_marginalization(self):
if 'luminosity_distance' not in self.prior.keys():
logging.info('No prior provided for distance, using default prior.')
self.prior['luminosity_distance'] = tupak.prior.create_default_prior('luminosity_distance')
self.distance_array = np.linspace(self.prior['luminosity_distance'].minimum,
self.prior['luminosity_distance'].maximum, 1e4)
self.delta_distance = self.distance_array[1] - self.distance_array[0]
self.distance_prior_array = np.array([self.prior['luminosity_distance'].prob(distance)
for distance in self.distance_array])
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')
self.bessel_function_interped = interp1d(np.linspace(0, 1e6, 1e5),
np.log([i0e(snr) for snr in np.linspace(0, 1e6, 1e5)])
+ np.linspace(0, 1e6, 1e5),
bounds_error=False, fill_value=-np.inf)
def noise_log_likelihood(self):
log_l = 0
for interferometer in self.interferometers:
......@@ -101,10 +46,7 @@ class MarginalizedLikelihood(Likelihood):
if waveform_polarizations is None:
return np.nan_to_num(-np.inf)
if self.time_marginalization:
signal_times_data = 0 * 1j
else:
matched_filter_snr_squared = 0
matched_filter_snr_squared = 0
optimal_snr_squared = 0
for interferometer in self.interferometers:
......@@ -142,4 +84,54 @@ class MarginalizedLikelihood(Likelihood):
def log_likelihood(self):
return self.log_likelihood() + self.noise_log_likelihood()
def setup_distance_marginalization(self):
if 'luminosity_distance' not in self.prior.keys():
logging.info('No prior provided for distance, using default prior.')
self.prior['luminosity_distance'] = tupak.prior.create_default_prior('luminosity_distance')
self.distance_array = np.linspace(self.prior['luminosity_distance'].minimum,
self.prior['luminosity_distance'].maximum, 1e4)
self.delta_distance = self.distance_array[1] - self.distance_array[0]
self.distance_prior_array = np.array([self.prior['luminosity_distance'].prob(distance)
for distance in self.distance_array])
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')
self.bessel_function_interped = interp1d(np.linspace(0, 1e6, 1e5),
np.log([i0e(snr) for snr in np.linspace(0, 1e6, 1e5)])
+ np.linspace(0, 1e6, 1e5),
bounds_error=False, fill_value=-np.inf)
class BasicLikelihood(object):
def __init__(self, interferometers, waveform_generator):
self.interferometers = interferometers
self.waveform_generator = waveform_generator
def noise_log_likelihood(self):
log_l = 0
for interferometer in self.interferometers:
log_l -= 2. / self.waveform_generator.time_duration * np.sum(
abs(interferometer.data) ** 2 / interferometer.power_spectral_density_array)
return log_l.real
def log_likelihood(self):
log_l = 0
waveform_polarizations = self.waveform_generator.frequency_domain_strain()
if waveform_polarizations is None:
return np.nan_to_num(-np.inf)
for interferometer in self.interferometers:
log_l += self.log_likelihood_interferometer(waveform_polarizations, interferometer)
return log_l.real
def log_likelihood_interferometer(self, waveform_polarizations, interferometer):
signal_ifo = interferometer.get_detector_response(waveform_polarizations, self.waveform_generator.parameters)
log_l = - 2. / self.waveform_generator.time_duration * np.vdot(interferometer.data - signal_ifo,
(interferometer.data - signal_ifo)
/ interferometer.power_spectral_density_array)
return log_l.real
def log_likelihood_ratio(self):
return self.log_likelihood() - self.noise_log_likelihood()
Markdown is supported
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