Commit ad26282f authored by Marc Arene's avatar Marc Arene
Browse files

bad commit.. ..but found the culprit for steppy logL !

parent 811e6c24
No preview for this file type
......@@ -174,6 +174,7 @@ def sigpar_to_injection_parameters(sigpar, start_time):
chi_2=0.0,
luminosity_distance=sigpar[3]/MPC,
iota=sigpar[0],
theta_jn=sigpar[0],
psi=sigpar[2],
# phase=0,
# phase=0.5 * (sigpar[1] - 2*np.pi*time_shift_H1/(sigpar[8] - time_shift_H1)),
......@@ -199,7 +200,10 @@ def injection_parameters_to_sigpar(injection_parameters, minimum_frequency):
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 180 / np.pi
inc = injection_parameters['iota'] * 180 / np.pi
if 'iota' in injection_parameters.keys():
inc = injection_parameters['iota'] * 180 / np.pi
else:
inc = injection_parameters['theta_jn'] * 180 / np.pi
psi = injection_parameters['psi'] * 180 / np.pi
phic = (injection_parameters['phase'] * 2 ) % (2*np.pi) * 180 / np.pi
dist = injection_parameters['luminosity_distance']
......@@ -219,7 +223,10 @@ def injection_parameters_to_binaryFile(injection_parameters, minimum_frequency,
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 180 / np.pi
inc = injection_parameters['iota'] * 180 / np.pi
if 'iota' in injection_parameters.keys():
inc = injection_parameters['iota'] * 180 / np.pi
else:
inc = injection_parameters['theta_jn'] * 180 / np.pi
psi = injection_parameters['psi'] * 180 / np.pi
phic = injection_parameters['phase'] * 2 * 180 / np.pi
dist = injection_parameters['luminosity_distance']
......@@ -242,7 +249,10 @@ def injection_parameters_to_binaryFile(injection_parameters, minimum_frequency,
def injection_parameters_to_q_pos(injection_parameters, start_time):
q_pos = np.zeros(13)
q_pos[0] = injection_parameters['iota']
if 'iota' in injection_parameters.keys():
q_pos[0] = injection_parameters['iota']
else:
q_pos[0] = injection_parameters['theta_jn']
# q_pos[1] = injection_parameters['phase'] * 2
q_pos[1] = (injection_parameters['phase'] * 2 ) % (2*np.pi)
q_pos[2] = injection_parameters['psi']
......@@ -266,6 +276,7 @@ def q_pos_to_dictionary_original(q_pos, start_time):
"""
dictionary = dict()
dictionary['iota'] = q_pos[0]
dictionary['theta_jn'] = q_pos[0]
dictionary['phase'] = q_pos[1] * 0.5
dictionary['psi'] = q_pos[2]
dictionary['luminosity_distance'] = q_pos[3] / MPC
......@@ -288,6 +299,7 @@ def q_pos_to_dictionary(q_pos, start_time, dictionary=None):
dictionary = set_inj.injection_parameters
if 0 in conf.param_indices:
dictionary['iota'] = q_pos[0]
dictionary['theta_jn'] = q_pos[0]
if 1 in conf.param_indices:
dictionary['phase'] = q_pos[1] * 0.5
if 2 in conf.param_indices:
......
......@@ -12,7 +12,7 @@ def WaveForm_ThreeDetectors(parameters, start_time, minimum_frequency, interfero
# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
waveform_arguments = dict(waveform_approximant='TaylorF2',
reference_frequency=0, minimum_frequency=minimum_frequency, f_max=0)
reference_frequency=0, minimum_frequency=minimum_frequency, maximum_frequency=0)
# Create the waveform_generator using a LAL Binary Neutron Star source function
waveform_generator = bilby.gw.WaveformGenerator(
......@@ -336,6 +336,9 @@ def dlogL_ThreeDetectors1(signal1, signal2, signal3, parameters, start_time, min
# Compute waveform and derivatives of waveform at initial position
h_wave_1, h_wave_2, h_wave_3, dh_Wave_1, dh_Wave_2, dh_Wave_3 = h_dh_Generation_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
# f_max = np.where(h_wave_1!=0)[0][-1] / parameters['meta']['duration']
# print("f_max = {}".format(f_max))
duration = interferometers.duration
# Compute gradient of log-likelihood. Cf equation (8.5) in Yann's
dlogL = np.zeros(9)
......
......@@ -29,16 +29,16 @@ param_inverse_functions = [np.arccos, identity, identity, np.exp, np.exp, np.exp
# param_offsets = [1e-7, 0, 1e-7, 0, 1e-7, 1e-7, 1e-7, 1e-2, 0] # 0 meaning analytical derivative
# param_offsets = [1e-6, 0, 1e-2, 0, 1e-5, 1e-6, 1e-3, 1e-2, 0] # 0 meaning analytical derivative
param_offsets = [1e-7, 1e-7, 1e-7, 0, 1e-7, 1e-7, 1e-5, 1e-7, 1e-7] # 0 meaning analytical derivative
param_offsets = [1e-7, 1e-7, 1e-7, 0, 1e-7, 1e-7, 1e-6, 1e-6, 1e-7] # 0 meaning analytical derivative
param_indices = [6]
param_indices = [3]
# param_indices = [0,1,2,3,4,5,8]
# param_indices = [0,1,2,3,4,5,6,7,8]
# param_indices = [7]
param_scales = [0,0,0,0,0,0,3.14e-03,0,0]
dlogL = 'dlogL2'
dlogL = 'dlogL1'
debug = False
use_conf_scales = False
......
......@@ -61,14 +61,17 @@ injection_parameters = dict(
chi_2=0.0,
luminosity_distance=40,
# iota=146*np.pi/180,
iota=2.81,
# iota=2.81,
theta_jn=2.81,
psi=2.21,
phase=5.497787143782138,
geocent_time=1187008882.43,
ra=3.44616, # value from fixed skyposition run
dec=-0.408084, # value from fixed skyposition run
lambda_1=0.,
lambda_2=0.
lambda_2=0.,
tilt_1=0,
tilt_2=0
)
injection_parameters = conv.generate_mass_parameters(injection_parameters)
......@@ -84,18 +87,22 @@ m2 = injection_parameters['mass_2']
tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
# We integrate the signal up to the frequency of the "Innermost stable circular orbit (ISCO)"
R_isco = 6. # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
m1_min = 1
m2_min = 1
f_high = Compute_fmax(R_isco, m1, m2)
f_high_run = Compute_fmax(R_isco, m1_min, m2_min)
# import IPython; IPython.embed();
duration = tc_3p5PN + 2
sampling_frequency = 4096
# sampling_frequency = 4096
# sampling_frequency = 4397
# sampling_frequency = 2 * fmax
sampling_frequency = 2 * f_high_run
start_time = injection_parameters['geocent_time'] - tc_3p5PN
meta_params = dict(
minimum_frequency = minimum_frequency,
tc_3p5PN = tc_3p5PN,
f_high = f_high,
f_high_run = f_high_run,
duration = duration,
sampling_frequency = sampling_frequency,
start_time = start_time
......
......@@ -748,3 +748,22 @@ Size of the scales:
- For each of the 9 parameters, plot
6. Prepare plots for bilby people at the LVC showing the steppiness of logL(sindec) and logL(ra)
7. Make merge request
Eric Thranes' input:
- artificially increase the time resolution by padding with zeros in the frequency domain for frequencies > ??
- set the noise to infinity for these frequencies, and signal to zero.
Eric Thranes' 2nd input:
- plot logL as a function of timeshift but keeping ra and dec constant...! This in order to see if the timescale with which logL changes on that plot is comparable to the timescale set by the sampling frequency
- Send Eric my code which plot the steppy logL(ra)
- plot with no noise: s = h
- keep all param constant except for `ln(D)`. Run a one dimensional leapfrog in `ln(D)`. 100 to 500 trajectories. Produce a ascii file
- float vs double precision in python: what is the default, which is best, why dtype=np.complex or np.cfloat ?
- from 1e-3 to 1e-5: gather informations on:
- time_shift
- ra and dec
- try using C-precision in python
No preview for this file type
......@@ -31,6 +31,7 @@ def main(signal1, signal2, signal3, injection_parameters, start_time, minimum_fr
sigpar, n = pu.injection_parameters_to_sigpar(injection_parameters, minimum_frequency)
n_freq = len(interferometers.frequency_array)
# n_freq = int(n/2)
# print("sigpar[14] = {}".format(sigpar[14]))
dlogL = tf2_for_bilby.dlogL_ThreeDetectors(signal1, signal2, signal3, sigpar, n_freq, interferometers)
prefix = '\ndlogL with hmc '
......@@ -40,9 +41,9 @@ def main(signal1, signal2, signal3, injection_parameters, start_time, minimum_fr
else:
info = prefix + 'using: dlogL = <s|dh> - <h|dh>'
print(info)
# print("{:15} | {:15}".format(" dlogL", "offsets"))
# for i in range(9):
# print("{:15.6e} | {:6.0e}\t".format(dlogL[i], conf.param_offsets[i]))
print("{:15} | {:15}".format(" dlogL", "offsets"))
for i in range(9):
print("{:15.6e} | {:6.0e}\t".format(dlogL[i], conf.param_offsets[i]))
return dlogL
......@@ -60,6 +61,9 @@ if __name__=='__main__':
parser.add_option("--new", default=False, action="store_true", help="""New noise realisation from PSDs is generated""", dest='new_noise_realisation_from_psd')
(opts,args)=parser.parse_args()
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# SET THE INJECTION PARAMETERS
injection_parameters = GW170817_parameters.main()
minimum_frequency = injection_parameters['meta']['minimum_frequency']
......
......@@ -5,6 +5,7 @@ import numpy as np
import bilby.gw.conversion as conv
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
from Library.Python_Utility import print_dict
import Library.set_injection_parameters as set_inj
# import sys
......@@ -57,7 +58,7 @@ def Compute_fmax(Rmin, m1, m2):
def main():
def main_old():
injection_parameters = dict(
mass_1=1.47,
mass_2=1.27,
......@@ -109,6 +110,8 @@ def main():
return injection_parameters
def main():
return set_inj.injection_parameters
if __name__=='__main__':
......
......@@ -36,7 +36,7 @@ injection_parameters = set_inj.injection_parameters
parameters = injection_parameters.copy()
nb_points = 100
nb_points = 1000
offset = conf.param_offsets[6]
sin_dec_inj = np.sin(injection_parameters['dec'])
......
# This script is meant to compare the influence of the offset used in the numerical derivative on trajectories that have revealed some instability.
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append('../')
import bilby
import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
import Library.config as conf
import Library.Python_Utility as pu
from Library.BNS_Utility_Functions import Mc_and_mu_to_M_eta_m1_and_m2
# import Library.set_injection_parameters as set_inj
import GW170817_parameters
import GW170817_set_psds
import GW170817_waveform_three_detectors
import GW170817_logL_snr
import GW170817_set_strain_data_from_psd
import GW170817_dlogL
def fparam_ranges(parameters, param_index):
"""
Returns the ranges for each parameter in which one can vary a parameter while keeping the other fixed at the position defined by the input `parameters`. When the bound is infinite we stop it at...
"""
Mc_min = parameters['reduced_mass'] / (0.2495**(2/5))
Mc_max = parameters['reduced_mass'] / (0.248**(2/5))
mu_max = 0.25**(2/5) * parameters['chirp_mass']
mu_min = 0.2435**(2/5) * parameters['chirp_mass']
ranges = [(-1,1), (0, 2*np.pi), (0, 2*np.pi), (54.55, 59), (np.log(Mc_min), np.log(Mc_max)), (np.log(mu_min), np.log(mu_max)), (-1, 0), (2, 4.5), (4.0413, 4.0418)]
return ranges[param_index]
def main(param_index, parameters, interferometers, signal1, signal2, signal3, start_time, minimum_frequency, generate_from_bilby, nb_points, whole_range, offsets):
offset = conf.param_offsets[param_index]
step = offset
if step == 0: # ie analytical derivative like for ln(D)
step = 1e-5
q_pos = pu.injection_parameters_to_q_pos(parameters, parameters['meta']['start_time'])
f_q_pos_pindex_inj = conf.param_functions[param_index](q_pos[param_index]) # does not work for phase and lntc...
# if param_index in [1,8]:
# print("PROBLEM: can't do plots for phic and lntc for the moment because I need to write special code for them....")
# sys.exit(1)
if whole_range:
f_q_pos_pindex_min, f_q_pos_pindex_max = fparam_ranges(parameters, param_index)
else:
f_q_pos_pindex_max = f_q_pos_pindex_inj + nb_points/2 * step
f_q_pos_pindex_min = f_q_pos_pindex_inj - nb_points/2 * step
f_q_pos_pindex_linspace = np.linspace(f_q_pos_pindex_min, f_q_pos_pindex_max, nb_points)
logL_list = []
for f_q_pos_pindex in f_q_pos_pindex_linspace:
q_pos[param_index] = conf.param_inverse_functions[param_index](f_q_pos_pindex)
parameters = pu.q_pos_to_dictionary(q_pos, start_time)
if param_index in [4,5]:
parameters['total_mass'], parameters['symmetric_mass_ratio'], parameters['mass_1'], parameters['mass_2'] = Mc_and_mu_to_M_eta_m1_and_m2(parameters['chirp_mass'], parameters['reduced_mass'])
# parameters.update({conf.param_dict_keys[param_index]: param_inverse_functions[param_index](f_q_pos_pindex)})
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, generate_from_bilby, False)
logL = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
logL_list.append(logL)
dlogL_list_offsets = []
for offset in offsets:
dlogL_list_offset = []
for f_q_pos_pindex in f_q_pos_pindex_linspace:
q_pos[param_index] = conf.param_inverse_functions[param_index](f_q_pos_pindex)
parameters = pu.q_pos_to_dictionary(q_pos, start_time)
if param_index in [4,5]:
parameters['total_mass'], parameters['symmetric_mass_ratio'], parameters['mass_1'], parameters['mass_2'] = Mc_and_mu_to_M_eta_m1_and_m2(parameters['chirp_mass'], parameters['reduced_mass'])
conf.param_offsets[param_index] = offset
dlogL = GW170817_dlogL.main(signal1, signal2, signal3, parameters, start_time, minimum_frequency, interferometers, generate_from_bilby, False)
dlogL_list_offset.append(dlogL[param_index])
dlogL_list_offsets.append(dlogL_list_offset)
return f_q_pos_pindex_linspace, logL_list, dlogL_list_offsets, f_q_pos_pindex_inj
if __name__=='__main__':
from optparse import OptionParser,OptionValueError
usage = """%prog [options]
Plots the logL and its derivative with respect to the chosen parameter."""
parser=OptionParser(usage)
parser.add_option("--psd", default=1, action="store", type="int", help="""1: sets the PSDs from official GWTC-1 open PSD data; 2: computes PSDs using Welch methods from gwosc open strain data; 3: uses bilby's analytical pre-stored PSD files.""", metavar="INT from {1,2,3}")
parser.add_option("--hmc", default=True, action="store_false", help="""Waveforms are generated using Ed and Yann's TF2 code, default is to generate them using bilby.lalsim""", dest='generate_from_bilby')
# parser.add_option("--offset", default=None,action="store", type="float", help="""offset used to compute dlogL but also stepsize between plotted points around the injeceted value.""")
parser.add_option("--new", default=False, action="store_true", help="""New noise realisation from PSDs is generated""", dest='new_noise_realisation_from_psd')
parser.add_option("--whole_range", default=False, action="store_true", help="""Plots on a predefined range that should give a good view of the logL shape overall""")
parser.add_option("-n", "--nb_points", default=500, action="store", type="int", help="""Number of points to plot. Default is 100""", metavar="INT")
parser.add_option("-i", "--p_index", default=None, action="store", type="int", help="""Index of the parameter to make the plot on. Default is the value stored in config.py""", metavar="INT")
parser.add_option("--no_seed", default=False, action="store_true", help="""New noise realisation from PSDs is generated""")
# parser.add_option("-p", "--param_index", default=6, action="store", type="int", help="""Index of the parameter which will be plotted""", metavar="INT from {0,1,2,3,4,5,6,7,8}")
(opts,args)=parser.parse_args()
if opts.p_index is None:
param_index = conf.param_indices[0]
conf.param_indices = [param_index]
else:
param_index = opts.p_index
conf.param_indices = [param_index]
if not(opts.no_seed):
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# offsets = [1e-6]
offsets = [1e-4, 1e-5, 1e-6, 1e-7]
# import IPython; IPython.embed(); sys.exit(1)
# SET THE INJECTION PARAMETERS
injection_parameters = GW170817_parameters.main()
minimum_frequency = injection_parameters['meta']['minimum_frequency']
duration = injection_parameters['meta']['duration']
sampling_frequency = injection_parameters['meta']['sampling_frequency']
f_high = injection_parameters['meta']['f_high']
start_time = injection_parameters['meta']['start_time']
# INITIALIZE THE THREE INTERFEROMETERS
interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
for ifo in interferometers:
ifo.minimum_frequency = minimum_frequency
# The following line is called automatically when setting the PSDs, which we donnot do here hence the line.
# Stangely however, in bilby, the `Interferometer` class does not have `duration` and `sampling_frequency` arguments when the `InterferometerStrainData` and the `InterferometerList` classes do. Hence one needs to set those arguments for the list of `InterferometerStrainData` so that the `InterferometerList` inherits from them. The values will be inherited from the first one in the list...
# Why do `interferometers.duration` and `interferometers.sampling_frequency` need to be set ? => Because I use them when creating my waveform_generator in bilby_wavefor
ifo.strain_data._set_time_and_frequency_array_parameters(duration, sampling_frequency, start_time)
# SET THEIR POWER SPECTRAL DENSITIES
# This function does not interpolate the psd
GW170817_set_psds.main(interferometers, opt_psd=opts.psd, sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
# SET NOISE STRAIN DATA FROM PSD
GW170817_set_strain_data_from_psd.main(interferometers, injection_parameters, sampling_frequency, duration, start_time, opts.psd, opts.new_noise_realisation_from_psd)
# WAVEFORM GENERATION
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(injection_parameters, start_time, minimum_frequency, interferometers, opts.generate_from_bilby)
# INJECT TEMPLATE INTO NOISE STRAIN DATA
# Add the template signal to each detector's strain data
interferometers[0].strain_data.frequency_domain_strain += h_wave_1
interferometers[1].strain_data.frequency_domain_strain += h_wave_2
interferometers[2].strain_data.frequency_domain_strain += h_wave_3
# COMPUTE AND PRINT SNR
logLBest, mf_snr_best, opt_snr_Best = GW170817_logL_snr.main(h_wave_1, h_wave_2, h_wave_3, interferometers, True)
# This copy is here just because for the moment I will keep manipulating 3 separate variables as originally done in the hmc
signal1 = interferometers[0].strain_data.frequency_domain_strain.copy()
signal2 = interferometers[1].strain_data.frequency_domain_strain.copy()
signal3 = interferometers[2].strain_data.frequency_domain_strain.copy()
# COMPUTE dlogL and print it to check with the plots
GW170817_dlogL.main(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers, opts.generate_from_bilby, True)
f_q_pos_pindex_linspace, logL_list, dlogL_list_offsets, f_q_pos_pindex_inj = main(param_index, injection_parameters, interferometers, signal1, signal2, signal3, start_time, minimum_frequency, opts.generate_from_bilby, opts.nb_points, opts.whole_range, offsets)
# import IPython; IPython.embed();
# PLOT
# offset = conf.param_offsets[param_index]
from matplotlib import rc
plt.rcParams['text.usetex'] = False
plt.rcParams['figure.titlesize'] = 18
plt.rcParams['figure.titleweight'] = 'bold'
plt.rcParams['axes.titlesize'] = 15
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['axes.titleweight'] = 'normal'
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['font.size'] = 10
if opts.generate_from_bilby:
title_prefix = '(bilby) '
else:
title_prefix = '(hmc) '
fig= plt.figure(figsize=(16,8))
# fig, axs = plt.subplots(2, 1, figsize=(16,8), gridspec_kw = {'height_ratios': [10, 1]})
# fig.suptitle('logL and dlogL computed with a step size of {:.0e}'.format(offset))
ax0 = plt.subplot(211)
ax0.plot(f_q_pos_pindex_linspace, logL_list, linewidth=0.5, color='green', label='logL')
# import IPython; IPython.embed()
ax0.axvline(x=f_q_pos_pindex_inj, linewidth=0.5, color='black', label='injected value')
ax0.set_xlabel('{}'.format(conf.fparam_names[param_index]))
ax0.set_ylabel('logL')
ax0.set_title(title_prefix + 'logL')
# ax0.legend(loc=(0.85,0.9))
ax0.legend()
ax1 = plt.subplot(212)
for i, offset in enumerate(offsets):
ax1.plot(f_q_pos_pindex_linspace, dlogL_list_offsets[i], linewidth=0.5, label='offset = {}'.format(offset))
ax1.axvline(x=f_q_pos_pindex_inj, linewidth=0.5, color='black', label='injected value')
ax1.set_xlabel('{}'.format(conf.fparam_names[param_index]))
ax1.set_ylabel('dlogL/d{}'.format(conf.fparam_names[param_index]))
ax1.set_title(title_prefix + 'dlogL/d{}'.format(conf.fparam_names[param_index]))
# ax1.tick_params('y', colors='red')
# ax1.legend(loc=(0.85,0.85))
ax1.legend()
# axs[0].plot(f_q_pos_pindex_linspace, logL_list, linewidth=1, color='green', label='logL')
# axs[0].axvline(x=f_q_pos_pindex_inj, linewidth=2, color='black')
# axs[0].set_xlabel('sin(dec)')
# axs[0].set_ylabel('logL')
# axs[0].set_title('logL and dlogL computed with a step size of {:.0e}'.format(offset))
#
# ax01 = axs[0].twinx()
# ax01.plot(f_q_pos_pindex_linspace, dlogL_list, linewidth=1, color='red', label='dlogL')
# ax01.set_ylabel('dlogL', color='red')
# ax01.tick_params('y', colors='red')
# axs[1].plot(f_q_pos_pindex_linspace, dlogL_list, linewidth=1, color='red', label='dlogL')
# axs[1].set_xlabel('sin(dec)')
# axs[1].set_title('dlogL computed with a step size of {:.0e}'.format(offset))
fig.tight_layout()
outdir = './.logL_and_dlogL_at_inj/'
if opts.generate_from_bilby:
suffix = 'bilby'
else:
suffix = 'hmc'
file_name = outdir + '{}_{}_{}'.format(param_index, conf.fparam_names[param_index], suffix)
fig.savefig(file_name)
# plt.show()
sys.exit(0)
# This script is meant to compare the influence of the offset used in the numerical derivative on trajectories that have revealed some instability.
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append('../')
import bilby
import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
import Library.config as conf
import Library.Python_Utility as pu
from Library.BNS_Utility_Functions import Mc_and_mu_to_M_eta_m1_and_m2
# import Library.set_injection_parameters as set_inj
import GW170817_parameters
import GW170817_set_psds
import GW170817_waveform_three_detectors
import GW170817_logL_snr
import GW170817_set_strain_data_from_psd
import GW170817_dlogL
def fparam_ranges(parameters, param_index):
"""
Returns the ranges for each parameter in which one can vary a parameter while keeping the other fixed at the position defined by the input `parameters`. When the bound is infinite we stop it at...
"""
Mc_min = parameters['reduced_mass'] / (0.249**(2/5))
Mc_max = parameters['reduced_mass'] / (0.245**(2/5))
mu_max = 0.249**(2/5) * parameters['chirp_mass']
mu_min = 0.240**(2/5) * parameters['chirp_mass']
ranges = [(-1,1), (0, 2*np.pi), (0, 2*np.pi), (54.55, 59), (np.log(Mc_min), np.log(Mc_max)), (np.log(mu_min), np.log(mu_max)), (-1, 1), (0, 2*np.pi), (4.041, 4.042)]
return ranges[param_index]
def main(param_index, parameters, interferometers, signal1, signal2, signal3, start_time, minimum_frequency, generate_from_bilby, nb_points, whole_range):
offset = conf.param_offsets[param_index]
step = offset
if step == 0:
step = 1e-5
q_pos = pu.injection_parameters_to_q_pos(parameters, parameters['meta']['start_time'])
f_q_pos_pindex_inj = conf.param_functions[param_index](q_pos[param_index]) # does not work for phase and lntc...
# if param_index in [1,8]:
# print("PROBLEM: can't do plots for phic and lntc for the moment because I need to write special code for them....")
# sys.exit(1)
if whole_range:
f_q_pos_pindex_min, f_q_pos_pindex_max = fparam_ranges(parameters, param_index)
else:
f_q_pos_pindex_max = f_q_pos_pindex_inj + nb_points/2 * step
f_q_pos_pindex_min = f_q_pos_pindex_inj - nb_points/2 * step
f_q_pos_pindex_linspace = np.linspace(f_q_pos_pindex_min, f_q_pos_pindex_max, nb_points)
logL_list = []
for f_q_pos_pindex in f_q_pos_pindex_linspace:
q_pos[param_index] = conf.param_inverse_functions[param_index](f_q_pos_pindex)
parameters = pu.q_pos_to_dictionary(q_pos, start_time)
if param_index in [4,5]:
parameters['total_mass'], parameters['symmetric_mass_ratio'], parameters['mass_1'], parameters['mass_2'] = Mc_and_mu_to_M_eta_m1_and_m2(parameters['chirp_mass'], parameters['reduced_mass'])
# parameters.update({conf.param_dict_keys[param_index]: param_inverse_functions[param_index](f_q_pos_pindex)})
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, generate_from_bilby, False)
logL = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
logL_list.append(logL)
return f_q_pos_pindex_linspace, logL_list, f_q_pos_pindex_inj
if __name__=='__main__':
from optparse import OptionParser,OptionValueError
usage = """%prog [options]
Plots the logL and its derivative with respect to the chosen parameter."""
parser=OptionParser(usage)
parser.add_option("--psd", default=1, action="store", type="int", help="""1: sets the PSDs from official GWTC-1 open PSD data; 2: computes PSDs using Welch methods from gwosc open strain data; 3: uses bilby's analytical pre-stored PSD files.""", metavar="INT from {1,2,3}")
parser.add_option("--hmc", default=True, action="store_false", help="""Waveforms are generated using Ed and Yann's TF2 code, default is to generate them using bilby.lalsim""", dest='generate_from_bilby')
# parser.add_option("--offset", default=None,action="store", type="float", help="""offset used to compute dlogL but also stepsize between plotted points around the injeceted value.""")
parser.add_option("--new", default=False, action="store_true", help="""New noise realisation from PSDs is generated""", dest='new_noise_realisation_from_psd')
parser.add_option("--whole_range", default=False, action="store_true", help="""Plots on a predefined range that should give a good view of the logL shape overall""")
parser.add_option("-n", "--nb_points", default=100, action="store", type="int", help="""Number of points to plot. Default is 100""", metavar="INT")
parser.add_option("--no_seed", default=False, action="store_true", help="""New noise realisation from PSDs is generated""")
# parser.add_option("-p", "--param_index", default=6, action="store", type="int", help="""Index of the parameter which will be plotted""", metavar="INT from {0,1,2,3,4,5,6,7,8}")
(opts,args)=parser.parse_args()
if len(conf.param_indices)!=1:
print("Error: please choose only one index to plot in `config.param_indices`")
sys.exit(1)
if not(opts.no_seed):
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# import IPython; IPython.embed(); sys.exit(1)
# SET THE INJECTION PARAMETERS
injection_parameters = GW170817_parameters.main()
minimum_frequency = injection_parameters['meta']['minimum_frequency']
duration = injection_parameters['meta']['duration']
sampling_frequency = injection_parameters['meta']['sampling_frequency']
f_high = injection_parameters['meta']['f_high']
start_time = injection_parameters['meta']['start_time']
# INITIALIZE THE THREE INTERFEROMETERS
interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
for ifo in interferometers:
ifo.minimum_frequency = minimum_frequency
# The following line is called automatically when setting the PSDs, which we donnot do here hence the line.
# Stangely however, in bilby, the `Interferometer` class does not have `duration` and `sampling_frequency` arguments when the `InterferometerStrainData` and the `InterferometerList` classes do. Hence one needs to set those arguments for the list of `InterferometerStrainData` so that the `InterferometerList` inherits from them. The values will be inherited from the first one in the list...
# Why do `interferometers.duration` and `interferometers.sampling_frequency` need to be set ? => Because I use them when creating my waveform_generator in bilby_wavefor
ifo.strain_data._set_time_and_frequency_array_parameters(duration, sampling_frequency, start_time)