Skip to content

Commit

Permalink
Merge pull request #10 from adivijaykumar/relbin_aditya
Browse files Browse the repository at this point in the history 
Changes to make max frequency consistent, remove zeros in likelihood evaluation, fix bug in summary data calculation
  • Loading branch information
@adivijaykumar
adivijaykumar authored Nov 25, 2020
2 parents e7bf340 + cf2afaa commit 645bf13
Showing 1 changed file with 82 additions and 75 deletions.
157 changes: 82 additions & 75 deletions bilby/gw/likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
import json
import copy
import sys
import matplotlib.pyplot as plt

import numpy as np
import scipy.integrate as integrate
Expand Down Expand Up @@ -1417,10 +1416,6 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
parameter_bounds: dict, optional
Dictionary of bounds (lists) for the initial parameters when finding
the initial maximum likelihood (fiducial) waveform.
min_bin_frequencty: int, optional
Minimum frequency of the bin range used.
max_bin_frequencty: int, optional
Maximum frequency of the bin range used.
chi: float, optional
Tunable parameter which limits the perturbation of alpha when setting
up the bin range. See https://arxiv.org/abs/1806.08792.
Expand All @@ -1437,8 +1432,7 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):

# Make sure that working with the individual polarizations still works...
def __init__(self, interferometers, waveform_generator,
initial_parameters={}, parameter_bounds={},
min_bin_frequency=20, max_bin_frequency=1000, chi=1,
initial_parameters={}, parameter_bounds={}, chi=1,
epsilon=.5, debug=False):
super(RelativeBinningGravitationalWaveTransient, self).__init__(
interferometers=interferometers,
Expand All @@ -1451,8 +1445,6 @@ def __init__(self, interferometers, waveform_generator,

self.initial_parameters = initial_parameters
self.parameter_bounds = parameter_bounds
self.min_bin_frequency = min_bin_frequency
self.max_bin_frequency = max_bin_frequency
self.chi = chi
self.epsilon = epsilon
self.gamma = np.array([-5 / 3, -2 / 3, 1, 5 / 3, 7 / 3])
Expand All @@ -1463,25 +1455,23 @@ def __init__(self, interferometers, waveform_generator,
self.fiducial_waveform_obtained = False
self.check_if_bins_are_setup = False
self.fiducial_polarizations = None
self.fiducial_polarizations_binned = None
self.per_detector_fiducial_waveforms = {}
self.bin_freqs = None
self.bin_inds = None
self.bin_freqs = dict()
self.bin_inds = dict()
self.initial_parameter_keys_sorted = sorted(self.initial_parameters)
self.maximum_likelihood_parameters = initial_parameters
self.setup_bins()
logger.info('Bin setup completed. Number of bins = {}'.format(len(self.bin_freqs) - 1))
self.set_fiducial_waveforms(self.initial_parameters)
logger.info("Initial fiducial waveforms set up")
self.setup_bins()
self.compute_summary_data()
logger.info("Summary Data Obtained")
#self.find_maximum_likelihood_waveform(self.initial_parameters, self.parameter_bounds, iterations=1)
# self.find_maximum_likelihood_waveform(self.initial_parameters, self.parameter_bounds, iterations=1)
# maxl_logl = self.log_likelihood_ratio_approx(None, parameter_dictionary=self.maximum_likelihood_parameters)
# maxl_logl = self.log_likelihood_ratio_approx()
# print(maxl_logl)

if debug:
print('maxl value = %s' % maxl_logl)
# print('maxl value = %s' % maxl_logl)
print('actual maxl value = %s' % self.log_likelihood_ratio_full(self.maximum_likelihood_parameters))

def __repr__(self):
Expand All @@ -1490,52 +1480,60 @@ def __repr__(self):

def setup_bins(self):
frequency_array = self.waveform_generator.frequency_array
frequency_array_useful = frequency_array[np.intersect1d(
np.where(frequency_array >= self.min_bin_frequency),
np.where(frequency_array <= self.max_bin_frequency))]
gamma = self.gamma

d_alpha = self.chi * 2 * np.pi / np.abs(
(self.min_bin_frequency ** gamma) * np.heaviside(
-gamma, 1) - (self.max_bin_frequency ** gamma) * np.heaviside(
gamma, 1))
d_phi = np.sum(np.array([np.sign(gamma[i]) * d_alpha[i] * (
frequency_array_useful ** gamma[i]) for i in range(len(gamma))]), axis=0)
d_phi_from_start = d_phi - d_phi[0]
num_bins = int(d_phi_from_start[-1] // self.epsilon)
self.bin_freqs = np.array([frequency_array_useful[np.where(d_phi_from_start >= (
(i / num_bins) * d_phi_from_start[-1]))[0][0]] for i in range(
num_bins + 1)])

self.bin_inds = np.array([np.where(frequency_array >= bin_freq)[0][0]
for bin_freq in self.bin_freqs])

self.waveform_generator.waveform_arguments['frequency_bin_edges'] = self.bin_freqs
for interferometer in self.interferometers:
frequency_array_useful = frequency_array[np.intersect1d(
np.where(frequency_array >= interferometer.minimum_frequency),
np.where(frequency_array <= interferometer.maximum_frequency))]

d_alpha = self.chi * 2 * np.pi / np.abs(
(interferometer.minimum_frequency ** gamma) * np.heaviside(
-gamma, 1) - (interferometer.maximum_frequency ** gamma) * np.heaviside(
gamma, 1))
d_phi = np.sum(np.array([np.sign(gamma[i]) * d_alpha[i] * (
frequency_array_useful ** gamma[i]) for i in range(len(gamma))]), axis=0)
d_phi_from_start = d_phi - d_phi[0]
self.number_of_bins = int(d_phi_from_start[-1] // self.epsilon)
self.bin_freqs[interferometer.name] = np.array([frequency_array_useful[np.where(d_phi_from_start >= (
(i / self.number_of_bins) * d_phi_from_start[-1]))[0][0]] for i in range(
self.number_of_bins + 1)])

self.bin_inds[interferometer.name] = np.array(
[np.where(frequency_array >= bin_freq)[0][0] for bin_freq in self.bin_freqs[interferometer.name]])
logger.info("Set up {} bins for {} between {} Hz and {} Hz".format(
self.number_of_bins, interferometer.name, interferometer.minimum_frequency,
interferometer.maximum_frequency))
self.waveform_generator.waveform_arguments["frequency_bin_edges"] = self.bin_freqs[interferometer.name]
return

def set_fiducial_waveforms(self, parameters):
self.waveform_generator.waveform_arguments["fiducial"] = True
self.fiducial_polarizations = self.waveform_generator.frequency_domain_strain(
parameters)

self.fiducial_polarizations_binned = {mode: (self.fiducial_polarizations[mode][self.bin_inds]
) for mode in (self.fiducial_polarizations.keys())}
maximum_nonzero_index = np.where(self.fiducial_polarizations["plus"] != 0j)[0][-1]
logger.info("Maximum Nonzero Index is {}".format(maximum_nonzero_index))
maximum_nonzero_frequency = self.waveform_generator.frequency_array[maximum_nonzero_index]
logger.info("Maximum Nonzero Frequency is {}".format(maximum_nonzero_frequency))

if self.fiducial_polarizations is None:
return np.nan_to_num(-np.inf)
plt.loglog(self.waveform_generator.frequency_array, self.fiducial_polarizations['plus'])

for interferometer in self.interferometers:
logger.info("Maximum Frequency is {}".format(interferometer.maximum_frequency))
if interferometer.maximum_frequency > maximum_nonzero_frequency:
interferometer.maximum_frequency = maximum_nonzero_frequency

self.per_detector_fiducial_waveforms[interferometer.name] = (
interferometer.get_detector_response(
self.fiducial_polarizations, parameters))
self.waveform_generator.waveform_arguments["fiducial"] = False
return

def log_likelihood(self):
return self.log_likelihood_ratio_relative_binning() + self.noise_log_likelihood()
return self.log_likelihood_ratio() + self.noise_log_likelihood()

def log_likelihood_ratio_relative_binning(self):
def log_likelihood_ratio(self):
d_inner_h = 0.
optimal_snr_squared = 0.
complex_matched_filter_snr = 0.
Expand All @@ -1556,7 +1554,6 @@ def log_likelihood_ratio_relative_binning(self):
# print('logl in inner calculation = ', log_l)
return float(log_l.real)


def find_maximum_likelihood_waveform(self, initial_parameter_guess,
parameter_bounds, iterations=10,
likelihood_threshold=1):
Expand Down Expand Up @@ -1586,14 +1583,13 @@ def get_best_fit_parameters(self, initial_parameter_bounds, maxiter=500,
atol=1e-10):
# Walk uphill using differential evolution from scipy.
print('computing maxL parameters...')
print('par bounds',initial_parameter_bounds)
print('par bounds', initial_parameter_bounds)
output = differential_evolution(
self.log_likelihood_ratio_relative_binning,
bounds=initial_parameter_bounds, atol=atol,
maxiter=maxiter)
best_fit = output['x']
log_likelihood = -output['fun']

# Output best-fit parameters if requested.
self.maximum_likelihood_parameters = (
self.get_parameter_dictionary_from_list(best_fit))
Expand All @@ -1620,35 +1616,46 @@ def compute_summary_data(self):
for interferometer in self.interferometers:
mask = interferometer.frequency_mask
masked_frequency_array = interferometer.frequency_array[mask]
maximum_bin_frequency_array = np.ones_like(masked_frequency_array)
start_index = 0
for edge in self.bin_freqs[1:]:
masked_bin_inds = []
for edge in self.bin_freqs[interferometer.name]:
index = np.where(masked_frequency_array == edge)[0][0]
maximum_bin_frequency_array[start_index:index +
1] = maximum_bin_frequency_array[start_index:index + 1] * edge
start_index = index + 1
factor = masked_frequency_array - maximum_bin_frequency_array

a0 = noise_weighted_inner_product(
interferometer.frequency_domain_strain[mask],
self.per_detector_fiducial_waveforms[interferometer.name][mask],
interferometer.power_spectral_density_array[mask],
self.waveform_generator.duration)
b0 = noise_weighted_inner_product(
self.per_detector_fiducial_waveforms[interferometer.name][mask],
self.per_detector_fiducial_waveforms[interferometer.name][mask],
interferometer.power_spectral_density_array[mask],
self.waveform_generator.duration)
a1 = noise_weighted_inner_product(
interferometer.frequency_domain_strain[mask] * factor,
self.per_detector_fiducial_waveforms[interferometer.name][mask],
interferometer.power_spectral_density_array[mask],
self.waveform_generator.duration)
b1 = noise_weighted_inner_product(
self.per_detector_fiducial_waveforms[interferometer.name][mask] * factor,
self.per_detector_fiducial_waveforms[interferometer.name][mask],
interferometer.power_spectral_density_array[mask],
self.waveform_generator.duration)
masked_bin_inds.append(index)
masked_strain = interferometer.frequency_domain_strain[mask]
masked_h0 = self.per_detector_fiducial_waveforms[interferometer.name][mask]
masked_psd = interferometer.power_spectral_density_array[mask]
a0, b0, a1, b1 = np.zeros((4, self.number_of_bins), dtype=np.complex)

for i in range(self.number_of_bins):

central_frequency_i = 0.5 * (masked_frequency_array[i] + masked_frequency_array[i + 1])
masked_strain_i = masked_strain[masked_bin_inds[i]:masked_bin_inds[i + 1]]
masked_h0_i = masked_h0[masked_bin_inds[i]:masked_bin_inds[i + 1]]
masked_psd_i = masked_psd[masked_bin_inds[i]:masked_bin_inds[i + 1]]
masked_frequency_i = masked_frequency_array[masked_bin_inds[i]:masked_bin_inds[i + 1]]

a0[i] = noise_weighted_inner_product(
masked_strain_i,
masked_h0_i,
masked_psd_i,
self.waveform_generator.duration)

b0[i] = noise_weighted_inner_product(
masked_h0_i,
masked_h0_i,
masked_psd_i,
self.waveform_generator.duration)

a1[i] = noise_weighted_inner_product(
masked_strain_i * (masked_frequency_i - central_frequency_i),
masked_h0_i,
masked_psd_i,
self.waveform_generator.duration)

b1[i] = noise_weighted_inner_product(
masked_h0_i * (masked_frequency_i - central_frequency_i),
masked_h0_i,
masked_psd_i,
self.waveform_generator.duration)

summary_data[interferometer.name] = dict(a0=a0, a1=a1, b0=b0, b1=b1)

Expand All @@ -1659,8 +1666,8 @@ def compute_relative_ratio(self, parameter_dictionary, interferometer):
self.waveform_generator.parameters = parameter_dictionary
new_polarizations = self.waveform_generator.frequency_domain_strain(parameter_dictionary)
h = interferometer.get_detector_response_relative_binning(
new_polarizations, parameter_dictionary, self.bin_freqs)
h0 = self.per_detector_fiducial_waveforms[interferometer.name][self.bin_inds]
new_polarizations, parameter_dictionary, self.bin_freqs[interferometer.name])
h0 = self.per_detector_fiducial_waveforms[interferometer.name][self.bin_inds[interferometer.name]]
waveform_ratio = h / h0

if (self.debug):
Expand All @@ -1676,7 +1683,7 @@ def compute_relative_ratio(self, parameter_dictionary, interferometer):

r0 = (waveform_ratio[1:] + waveform_ratio[:-1]) / 2
r1 = (waveform_ratio[1:] - waveform_ratio[:-1]) / (
self.bin_freqs[1:] - self.bin_freqs[:-1])
self.bin_freqs[interferometer.name][1:] - self.bin_freqs[interferometer.name][:-1])
return r0, r1

def calculate_snrs_from_summary_data(self, summary_data_per_interferometer, r0, r1):
Expand Down

0 comments on commit 645bf13

Please sign in to comment.