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

define log_likelihood_ratio

parents 6adb8053 edf9b93c
Pipeline #15867 passed with stages
in 3 minutes and 3 seconds
......@@ -75,27 +75,27 @@ class Interferometer:
detector_tensor = 0.5 * (np.einsum('i,j->ij', self.x, self.x) - np.einsum('i,j->ij', self.y, self.y))
return detector_tensor
def inject_signal(self, source, params):
def inject_signal(self, source):
"""
Inject a signal into noise.
Adds the requested signal to self.data
:param source: source type
:param params: parameters
"""
signal = source.frequency_domain_strain(params)
signal = source.frequency_domain_strain()
for mode in signal.keys():
det_response = self.antenna_response(
params['ra'], params['dec'], params['geocent_time'],
params['psi'], mode)
source.ra, source.dec, source.geocent_time,
source.psi, mode)
signal[mode] *= det_response
signal_ifo = sum(signal.values())
time_shift = self.time_delay_from_geocenter(
params['ra'], params['dec'], params['geocent_time'])
source.ra, source.dec, source.geocent_time)
signal_ifo *= np.exp(-1j*2*np.pi*time_shift*source.frequency_array)
self.data += signal_ifo
......
import numpy as np
class likelihood:
class Likelihood:
def __init__(self, interferometers, source):
self.interferometers = interferometers
self.source = source
self.parameter_keys = set(self.source.parameter_keys +
['ra', 'dec', 'geocent_time', 'psi'])
self.noise_log_likelihood = 0
self.set_noise_log_likelihood()
def loglikelihood(self, parameters):
log_l = -self.noise_log_likelihood
waveform_polarizations = self.source.frequency_domain_strain(parameters)
for interferometer in self.interferometers:
h = []
for mode in waveform_polarizations:
det_response = interferometer.antenna_response(
parameters['ra'], parameters['dec'],
parameters['geocent_time'], parameters['psi'], mode)
def get_interferometer_signal(self, waveform_polarizations, interferometer):
h = []
for mode in waveform_polarizations:
det_response = interferometer.antenna_response(
self.source.ra, self.source.dec,
self.source.geocent_time, self.source.psi, mode)
h.append(waveform_polarizations[mode] * det_response)
signal = np.sum(h, axis=0)
h.append(waveform_polarizations[mode] * det_response)
time_shift = interferometer.time_delay_from_geocenter(
self.source.ra, self.source.dec,
self.source.geocent_time)
signal *= np.exp(-1j * 2 * np.pi * time_shift * self.source.frequency_array)
signal_ifo = np.sum(h, axis=0)
return signal
time_shift = interferometer.time_delay_from_geocenter(
parameters['ra'], parameters['dec'],
parameters['geocent_time'])
signal_ifo *= np.exp(-1j*2*np.pi*time_shift*self.source.frequency_array)
def log_likelihood(self):
log_l = 0
waveform_polarizations = self.source.frequency_domain_strain()
for interferometer in self.interferometers:
log_l += self.log_likelihood_interferometer(waveform_polarizations, interferometer)
return log_l.real
log_l -= 4. / self.source.time_duration * np.vdot(
interferometer.data - signal_ifo,
(interferometer.data - signal_ifo) / (
interferometer.power_spectral_density_array))
def log_likelihood_interferometer(self, waveform_polarizations, interferometer):
signal_ifo = self.get_interferometer_signal(waveform_polarizations, interferometer)
log_l = - 4. / self.source.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
def set_noise_log_likelihood(self):
log_l = 0
for interferometer in self.interferometers:
......@@ -44,40 +50,39 @@ class likelihood:
self.noise_log_likelihood = log_l.real
class likelihoodB(likelihood):
class LikelihoodB(Likelihood):
def __init__(self, interferometers, source):
likelihood.__init__(self, interferometers, source)
Likelihood.__init__(self, interferometers, source)
for interferometer in self.interferometers:
interferometer.whiten_data()
def loglikelihood(self, parameters):
def log_likelihood(self):
log_l = 0
waveform_polarizations = self.source.frequency_domain_strain(parameters)
waveform_polarizations = self.source.frequency_domain_strain()
for interferometer in self.interferometers:
for mode in waveform_polarizations.keys():
det_response = interferometer.antenna_response(
parameters['ra'], parameters['dec'],
parameters['geocent_time'], parameters['psi'], mode)
self.source.ra, self.source.dec,
self.source.geocent_time, self.source.psi, mode)
waveform_polarizations[mode] *= det_response
signal_ifo = np.sum(waveform_polarizations.values(), axis=0)
time_shift = interferometer.time_delay_from_geocenter(
parameters['ra'], parameters['dec'],
parameters['geocent_time'])
signal_ifo *= np.exp(-1j*2*np.pi*time_shift)
self.source.ra, self.source.dec,
self.source.geocent_time)
signal_ifo *= np.exp(-1j * 2 * np.pi * time_shift)
signal_ifo_whitened = signal_ifo / (
interferometer.amplitude_spectral_density_array)
interferometer.amplitude_spectral_density_array)
log_l -= 4. * self.source.sampling_frequency * (
np.real(sum(
(interferometer.whitened_data - signal_ifo_whitened)**2)))
(interferometer.whitened_data - signal_ifo_whitened) ** 2)))
return log_l
......@@ -9,8 +9,6 @@ import peyote
class Result(dict):
def __init__(self):
pass
def __getattr__(self, name):
try:
......@@ -34,8 +32,8 @@ class Result(dict):
if os.path.isfile(file_name):
logging.info(
'Renaming existing file {} to {}.old'
.format(file_name, file_name))
os.rename(file_name, file_name+'.old')
.format(file_name, file_name))
os.rename(file_name, file_name + '.old')
logging.info("Saving result to {}".format(file_name))
with open(file_name, 'w+') as f:
......@@ -71,12 +69,20 @@ class Sampler:
self.label = label
self.outdir = outdir
self.kwargs = kwargs
self.sampler_string = sampler_string
self.external_sampler = None
self.import_external_sampler()
self.active_parameter_values = self.prior.__dict__.copy()
self.search_parameter_keys = []
self.ndim = 0
self.initialise_parameters()
self.verify_prior()
self.result = Result()
self.add_initial_data_to_results()
self.set_kwargs()
......@@ -89,25 +95,22 @@ class Sampler:
pass
def add_initial_data_to_results(self):
self.result = Result()
self.result.search_parameter_keys = self.search_parameter_keys
self.result.labels = [
self.prior[k].latex_label for k in self.search_parameter_keys]
self.result.labels = [self.prior.__dict__[k].latex_label for k in self.search_parameter_keys]
def initialise_parameters(self):
self.fixed_parameters = self.prior.copy()
self.search_parameter_keys = []
for key in self.likelihood.parameter_keys:
if key in self.prior:
p = self.prior[key]
CA = isinstance(p, numbers.Real)
CB = hasattr(p, 'prior')
CC = getattr(p, 'is_fixed', False) is True
if CA is False and CB and CC is False:
for key in self.likelihood.source.__dict__:
if key in self.prior.__dict__:
p = self.prior.__dict__[key]
ca = isinstance(p, numbers.Real)
cb = hasattr(p, 'prior')
cc = getattr(p, 'is_fixed', False) is True
if ca is False and cb and cc is False:
self.search_parameter_keys.append(key)
self.fixed_parameters[key] = np.nan
elif CC:
self.fixed_parameters[key] = p.value
self.active_parameter_values[key] = np.nan
elif cc:
self.active_parameter_values[key] = p.value
else:
try:
self.prior[key] = getattr(peyote.parameter, key)
......@@ -119,31 +122,31 @@ class Sampler:
logging.info("Search parameters:")
for key in self.search_parameter_keys:
logging.info(' {} ~ {}'.format(key, self.prior[key].prior))
logging.info(' {} ~ {}'.format(key, self.prior.__dict__[key]))
def verify_prior(self):
required_keys = self.likelihood.parameter_keys
required_keys = self.likelihood.source.__dict__
unmatched_keys = [
r for r in required_keys if r not in self.prior]
r for r in required_keys if r not in self.prior.__dict__]
if len(unmatched_keys) > 0:
raise ValueError(
"Input prior is missing keys {}".format(unmatched_keys))
def prior_transform(self, theta):
return [self.prior[k].prior.rescale(t)
for k, t in zip(self.search_parameter_keys, theta)]
return [self.prior.__dict__[key].prior.rescale(t)
for key, t in zip(self.search_parameter_keys, theta)]
def loglikelihood(self, theta):
for i, k in enumerate(self.search_parameter_keys):
self.fixed_parameters[k] = theta[i]
return self.likelihood.loglikelihood(self.fixed_parameters)
self.likelihood.source.__dict__[k] = theta[i]
return self.likelihood.log_likelihood()
def run_sampler(self):
pass
def import_external_sampler(self):
try:
self.extenal_sampler = __import__(self.sampler_string)
self.external_sampler = __import__(self.sampler_string)
except ImportError:
raise ImportError(
"Sampler {} not installed on this system".format(
......@@ -162,7 +165,7 @@ class Nestle(Sampler):
self.kwargs = self.kwargs_defaults
def run_sampler(self):
nestle = self.extenal_sampler
nestle = self.external_sampler
if self.kwargs.get('verbose', True):
self.kwargs['callback'] = nestle.print_progress
......@@ -180,7 +183,7 @@ class Nestle(Sampler):
class Dynesty(Sampler):
def run_sampler(self):
dynesty = self.extenal_sampler
dynesty = self.external_sampler
nested_sampler = dynesty.NestedSampler(
loglikelihood=self.loglikelihood,
prior_transform=self.prior_transform,
......@@ -210,7 +213,7 @@ class Pymultinest(Sampler):
self.kwargs['outputfiles_basename'])
def run_sampler(self):
pymultinest = self.extenal_sampler
pymultinest = self.external_sampler
out = pymultinest.solve(
LogLikelihood=self.loglikelihood, Prior=self.prior_transform,
n_dims=self.ndim, **self.kwargs)
......@@ -226,13 +229,12 @@ class Pymultinest(Sampler):
def run_sampler(likelihood, prior, label='label', outdir='outdir',
sampler='nestle', **sampler_kwargs):
if hasattr(peyote.sampler, sampler.title()):
SamplerClass = getattr(peyote.sampler, sampler.title())
sampler = SamplerClass(likelihood, prior, sampler, outdir=outdir,
label=label, **sampler_kwargs)
sampler_class = getattr(peyote.sampler, sampler.title())
sampler = sampler_class(likelihood, prior, sampler, outdir=outdir,
label=label, **sampler_kwargs)
result = sampler.run_sampler()
result.save_to_file(outdir=outdir, label=label)
return result
else:
raise ValueError(
"Sampler {} not yet implemented".format(sampler))
from __future__ import division, print_function
import peyote
import numpy as np
import os.path
import os
import lal
import lalsimulation as lalsim
from astropy.table import Table
from peyote.utils import sampling_frequency, nfft
try:
import lal
import lalsimulation as lalsim
except ImportError:
print("lal is not installed")
import peyote
from peyote.utils import sampling_frequency, nfft
class Source:
......@@ -49,41 +47,59 @@ class BinaryBlackHole(Source):
"""
A way of getting a BBH waveform from lal
"""
parameter_keys = ['mass_1', 'mass_2', 'luminosity_distance', 'spin_1',
'spin_1', 'iota', 'phase', 'waveform_approximant',
'reference_frequency', 'ra', 'dec', 'geocent_time',
'psi']
def frequency_domain_strain(self, parameters):
luminosity_distance = parameters['luminosity_distance'] * 1e6 * lal.PC_SI
mass_1 = parameters['mass_1'] * lal.MSUN_SI
mass_2 = parameters['mass_2'] * lal.MSUN_SI
def __init__(self, name, sampling_frequency, time_duration, mass_1, mass_2, luminosity_distance, spin_1, spin_2,
iota, phase, waveform_approximant, reference_frequency, ra, dec, geocent_time, psi):
Source.__init__(self, name, sampling_frequency, time_duration)
self.mass_1 = mass_1
self.mass_2 = mass_2
self.luminosity_distance = luminosity_distance
self.spin_1 = spin_1
self.spin_2 = spin_2
self.iota = iota
self.phase = phase
self.waveform_approximant = waveform_approximant
self.reference_frequency = reference_frequency
self.ra = ra
self.dec = dec
self.geocent_time = geocent_time
self.psi = psi
def copy(self):
copy = BinaryBlackHole(self.name, self.sampling_frequency, self.time_duration, self.mass_1, self.mass_2,
self.luminosity_distance, self.spin_1, self.spin_2, self.iota, self.phase,
self.waveform_approximant, self.reference_frequency, self.ra, self.dec,
self.geocent_time, self.psi)
return copy
def frequency_domain_strain(self):
luminosity_distance = self.luminosity_distance * 1e6 * lal.PC_SI
mass_1 = self.mass_1 * lal.MSUN_SI
mass_2 = self.mass_2 * lal.MSUN_SI
longitude_ascending_nodes = 0.0
eccentricity = 0.0
meanPerAno = 0.0
mean_per_ano = 0.0
waveform_dictionary = lal.CreateDict()
approximant = lalsim.GetApproximantFromString(
parameters['waveform_approximant'])
approximant = lalsim.GetApproximantFromString(self.waveform_approximant)
frequency_minimum = 20
frequency_maximum = self.frequency_array[-1]
delta_frequency = self.frequency_array[1] - self.frequency_array[0]
hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
mass_1, mass_2, parameters['spin_1'][0], parameters['spin_1'][1],
parameters['spin_1'][2], parameters['spin_2'][0],
parameters['spin_2'][1], parameters['spin_2'][2],
luminosity_distance, parameters['iota'],
parameters['phase'], longitude_ascending_nodes,
eccentricity, meanPerAno, delta_frequency, frequency_minimum,
frequency_maximum, parameters['reference_frequency'],
h_plus, h_cross = lalsim.SimInspiralChooseFDWaveform(
mass_1, mass_2, self.spin_1[0], self.spin_1[1], self.spin_1[2],
self.spin_2[0], self.spin_2[1], self.spin_2[2],
luminosity_distance, self.iota,
self.phase, longitude_ascending_nodes,
eccentricity, mean_per_ano, delta_frequency, frequency_minimum,
frequency_maximum, self.reference_frequency,
waveform_dictionary, approximant)
h_plus = hplus.data.data
h_cross = hcross.data.data
h_plus = h_plus.data.data
h_cross = h_cross.data.data
return {'plus': h_plus, 'cross': h_cross}
......@@ -122,15 +138,6 @@ class Supernova(AstrophysicalSource):
AstrophysicalSource.__init__(self, name, right_ascension, declination, luminosity_distance)
# class BinaryBlackHole(CompactBinaryCoalescence):
# def __init__(self, name, right_ascension, declination, luminosity_distance, mass_1, mass_2, spin_1, spin_2,
# coalescence_time, inclination_angle, waveform_phase, polarisation_angle, eccentricity):
# CompactBinaryCoalescence.__init__(self, name, right_ascension, declination, luminosity_distance, mass_1,
# mass_2, spin_1, spin_2, coalescence_time, inclination_angle, waveform_phase,
# polarisation_angle, eccentricity)
class BinaryNeutronStar(CompactBinaryCoalescence):
def __init__(self, name, right_ascension, declination, luminosity_distance, mass_1, mass_2, spin_1, spin_2,
coalescence_time, inclination_angle, waveform_phase, polarisation_angle, eccentricity,
......@@ -171,10 +178,10 @@ class BinaryNeutronStarMergerNumericalRelativity(Source):
full_filename = '{}/{}'.format(directory_path, file_name)
if not os.path.isfile(full_filename):
print('{} does not exist'.format(full_filename)) # add exception
return(-1)
else: # ok file exists
print('{} does not exist'.format(full_filename)) # add exception
return (-1)
else: # ok file exists
strain_table = Table.read(full_filename)
Hplus, _ = nfft(strain_table["hplus"], sampling_frequency(strain_table['time']))
Hcross, frequency = nfft(strain_table["hcross"], sampling_frequency(strain_table['time']))
return(strain_table['time'],strain_table["hplus"],strain_table["hcross"],frequency,Hplus,Hcross)
return (strain_table['time'], strain_table["hplus"], strain_table["hcross"], frequency, Hplus, Hcross)
......@@ -45,12 +45,12 @@ class Test(unittest.TestCase):
self.assertAlmostEqual(all(self.msd['hf_signal_and_noise'] - hf_signal_and_noise_saved), 0.00000000, 5)
def test_recover_luminosity_distance(self):
likelihood = peyote.likelihood.likelihood(
likelihood = peyote.likelihood.Likelihood(
[self.msd['IFO']], self.msd['source'])
prior = self.msd['simulation_parameters'].copy()
dL = self.msd['simulation_parameters']['luminosity_distance']
prior['luminosity_distance'] = peyote.parameter.Parameter(
prior = self.msd['source'].copy()
dL = self.msd['source'].luminosity_distance
prior.luminosity_distance = peyote.parameter.Parameter(
'luminosity_distance',
prior=peyote.prior.Uniform(lower=dL-10, upper=dL+10))
......
import numpy as np
import pylab as plt
import dynesty.plotting as dyplot
import corner
import peyote
import dynesty.plotting as dyplot
import corner
peyote.setup_logger()
# peyote.setup_logging()
time_duration = 1.
sampling_frequency = 4096.
simulation_parameters = dict(
mass_1=36., mass_2=29.,
spin_1=[0, 0, 0],
spin_2=[0, 0, 0],
luminosity_distance=100.,
iota=0.4, #np.pi/2,
phase=1.3,
waveform_approximant='IMRPhenomPv2',
reference_frequency=50.,
ra=1.375,
dec=-1.2108,
geocent_time=1126259642.413,
psi=2.659
)
source = peyote.source.BinaryBlackHole('BBH', sampling_frequency, time_duration)
hf_signal = source.frequency_domain_strain(simulation_parameters)
source = peyote.source.BinaryBlackHole('BBH', sampling_frequency, time_duration, mass_1=36., mass_2=29.,
spin_1=[0, 0, 0], spin_2=[0, 0, 0], luminosity_distance=100., iota=0.4,
phase=1.3, waveform_approximant='IMRPhenomPv2', reference_frequency=50.,
ra=1.375, dec=-1.2108, geocent_time=1126259642.413, psi=2.659)
# source = peyote.source.BinaryBlackHole('BBH', sampling_frequency, time_duration)
hf_signal = source.frequency_domain_strain()
# Simulate the data in H1
H1 = peyote.detector.H1
H1_hf_noise, frequencies = H1.power_spectral_density.get_noise_realisation(
sampling_frequency, time_duration)
H1.set_data(sampling_frequency, time_duration, frequency_domain_strain=H1_hf_noise)
H1.inject_signal(source, simulation_parameters)
H1.inject_signal(source)
H1.set_spectral_densities()
# Simulate the data in L1
......@@ -41,7 +36,7 @@ L1 = peyote.detector.L1
L1_hf_noise, frequencies = L1.power_spectral_density.get_noise_realisation(
sampling_frequency, time_duration)
L1.set_data(sampling_frequency, time_duration, frequency_domain_strain=L1_hf_noise)
L1.inject_signal(source, simulation_parameters)
L1.inject_signal(source)
L1.set_spectral_densities()
IFOs = [H1, L1]
......@@ -57,32 +52,26 @@ plt.xlabel(r'frequency')
plt.ylabel(r'strain')
fig.savefig('data')
likelihood = peyote.likelihood.likelihood(IFOs, source)
likelihood = peyote.likelihood.Likelihood(IFOs, source)
prior = simulation_parameters.copy()
prior['mass_1'] = peyote.parameter.Parameter(
prior = source.copy()
prior.mass_1 = peyote.parameter.Parameter(
'mass_1', prior=peyote.prior.Uniform(lower=35, upper=37),
latex_label='$m_1$')
#prior['mass_2'] = peyote.parameter.Parameter(
# 'mass_2', prior=peyote.prior.Uniform(lower=27, upper=31),
# latex_label='$m_2$')
#prior['iota'] = peyote.parameter.Parameter(
# prior.mass_2 = peyote.parameter.Parameter(
# 'mass_2', prior=peyote.prior.Uniform(lower=27, upper=31),
# latex_label='$m_2$')
#prior.iota = peyote.parameter.Parameter(
# 'iota', prior=peyote.prior.Uniform(lower=0, upper=np.pi),
# latex_label='$iota$')
#prior['phase'] = peyote.parameter.Parameter(
# 'phase', prior=peyote.prior.Uniform(lower=0, upper=np.pi))
#prior['geocent_time'] = peyote.parameter.Parameter(
# 'geocent_time', prior=peyote.prior.Uniform(
# lower=simulation_parameters['geocent_time']-100,
# upper=simulation_parameters['geocent_time']+100))
prior['luminosity_distance'] = peyote.parameter.Parameter(
prior.luminosity_distance = peyote.parameter.Parameter(
'luminosity_distance', prior=peyote.prior.Uniform(lower=30, upper=200),
latex_label='$d_L$')
result = peyote.run_sampler(likelihood, prior, sampler='nestle',
n_live_points=200, verbose=True)
truths = [simulation_parameters[x] for x in result.search_parameter_keys]
truths = [source.__dict__[x] for x in result.search_parameter_keys]
fig = corner.corner(result.samples, truths=truths, labels=result.search_parameter_keys)
fig.savefig('corner')
......
# coding: utf-8
# # GW150914 analysis
......@@ -15,12 +14,12 @@ import peyote
import corner
import logging
logging.getLogger().addHandler(logging.StreamHandler())
logging.getLogger().setLevel('DEBUG')
import matplotlib.mlab as mlab
data_file = 'tutorial_data/GW150914_strain_data.npy'
time_series, strain_H1, strain_L1 = np.load(data_file)
time_duration = time_series[-1] - time_series[0]
......@@ -29,9 +28,8 @@ time_of_event = 1126259462.44