Commit 28563df7 authored by Marc Arene's avatar Marc Arene
Browse files

Copied all the /Tests/GW170817_.py files in /Library.

parent dabc0d48
......@@ -28,13 +28,12 @@ import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
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
import GW170817_FIM
import Library.set_psds as set_psds
import Library.waveform_three_detectors as waveform_three_detectors
import Library.logL_snr as logL_snr
import Library.dlogL as dlogL
import Library.FIM as FIM
if __name__=='__main__':
# previously was : run BNS_HMC_9D.py input_file #n_traj_HMC_tot #n_traj_fit #n_1_fit #n_2_fit #source_id SkipPhase1=False
......@@ -199,17 +198,16 @@ if __name__=='__main__':
# 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_psds.main(interferometers, opt_psd=opts.psd, sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
# SET NOISE STRAIN DATA FROM PSD
interferometers.set_strain_data_from_power_spectral_densities(
sampling_frequency=sampling_frequency,
duration=duration,
start_time=start_time)
# GW170817_set_strain_data_from_psd.main(interferometers, injection_parameters, sampling_frequency, duration, start_time, opts.psd, False)
# WAVEFORM GENERATION
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(injection_parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby)
h_wave_1, h_wave_2, h_wave_3 = waveform_three_detectors.main(injection_parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby)
# INJECT TEMPLATE INTO NOISE STRAIN DATA
# Add the template signal to each detector's strain data
......@@ -218,7 +216,7 @@ if __name__=='__main__':
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)
logLBest, mf_snr_best, opt_snr_Best = logL_snr.main(h_wave_1, h_wave_2, h_wave_3, interferometers, True)
# snrBest = np.sqrt(2.0*logLBest)
logL = logLBest # Values of LogL for initial point
......@@ -230,7 +228,7 @@ if __name__=='__main__':
signal3 = interferometers[2].strain_data.frequency_domain_strain.copy()
# COMPUTE dlogL and print it to check with the plots
dlogL = GW170817_dlogL.main(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, True)
dlogL = dlogL.main(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, True)
# Define boundary conditions for HMC using values from input files
# The boundary conditions are given both for masses m1, m2 (solar masses) and for the distance (meters)
......@@ -252,7 +250,7 @@ if __name__=='__main__':
print("m2min = {}, m2max = {}".format(boundary[1][0], boundary[1][1]))
print("Dmin = {}, Dmax = {}".format(boundary[2][0]/MPC, boundary[2][1]/MPC))
FisherMatrix, FisherMatrix_hmc, scale, scale_SVD, scale_hmc = GW170817_FIM.main(injection_parameters, start_time, minimum_frequency, interferometers)
FisherMatrix, FisherMatrix_hmc, scale, scale_SVD, scale_hmc = FIM.main(injection_parameters, start_time, minimum_frequency, interferometers)
if not(conf.generate_from_bilby):
scale = scale_hmc
......
......@@ -16,8 +16,8 @@ import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
import Library.config as conf
import GW170817_waveform_three_detectors
import GW170817_dlogL
import Library.waveform_three_detectors as waveform_three_detectors
import Library.dlogL as dlogL
def FindClosestPoints(q_pos_fit, PointFit_Sort, n_pt_fit, n_param, n_1_fit, param_index):
"""
......@@ -248,7 +248,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, si
# QUESTION: should I pass as an input the injection_parameters dictionary and update it here with the following parameters dictionary or is it sufficient to just give as input to bilby/lal this parameters dictionary?
parameters = su.q_pos_to_dictionary(q_pos, start_time)
# h_wave_1, h_wave_2, h_wave_3 = bilby_wv.WaveForm_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
h_wave_1, h_wave_2, h_wave_3 = waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
# Calculate Hamiltonian at final position, and evaluate Hastings ratio with appropriate priors
# logL_final, SNR_final = buf.Log_Likelihood_SNR_ThreeDetectors(signal_1, signal_2, signal_3, h_wave_1, h_wave_2, h_wave_3, shf_L, shf_V, q_pos, n_time)
......@@ -363,7 +363,7 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
if conf.debug:
parameters = su.q_pos_to_dictionary(q_pos, start_time)
# h_wave_1, h_wave_2, h_wave_3 = bilby_wv.WaveForm_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
h_wave_1, h_wave_2, h_wave_3 = waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
logL = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
H = computeHamiltonian(p_mom,logL)
......@@ -476,7 +476,7 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
# QUESTION: should I pass as an input the injection_parameters dictionary and update it here with the following parameters dictionary or is it sufficient to just give as input to bilby/lal this parameters dictionary?
parameters = su.q_pos_to_dictionary(q_pos, start_time)
# dlogL_pos = bilby_wv.dlogL_ThreeDetectors(signal_1, signal_2, signal_3, parameters, start_time, minimum_frequency, interferometers, l, i)
dlogL_pos = GW170817_dlogL.main(signal_1, signal_2, signal_3, parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
dlogL_pos = dlogL.main(signal_1, signal_2, signal_3, parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
# dlogL_pos = tf2_for_bilby.dlogL_ThreeDetectors(signal_1, signal_2, signal_3, q_pos, n_freq, interferometers)
# 1/2 step in momentum (leapfrog algorithm)
......@@ -497,7 +497,7 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
if conf.debug:
# h_wave_1, h_wave_2, h_wave_3 = bilby_wv.WaveForm_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers)
h_wave_1, h_wave_2, h_wave_3 = GW170817_waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
h_wave_1, h_wave_2, h_wave_3 = waveform_three_detectors.main(parameters, start_time, minimum_frequency, interferometers, conf.generate_from_bilby, False)
logL = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
H = computeHamiltonian(p_mom,logL)
......
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import numpy as np
import bilby
import Library.bilby_waveform as bilby_wv
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.sigpar_utils as su
import Library.TaylorF2_for_bilby as tf2_for_bilby
import Library.config as conf
import Library.set_injection_parameters as set_inj
import Library.set_psds as set_psds
import Library.waveform_three_detectors as waveform_three_detectors
def main(injection_parameters, start_time, minimum_frequency, interferometers):
######################################
# FISHER MATRICES #
###################
# classic decomp
# Prepare scales used in HMC equations using the Fisher matrix depending on option taken:
# option = 1 : full inversion of Fisher matrix
# other option : inverse of square root of diagonal FIM terms
# Compute the Fisher Matrix
# FisherMatrix = tf2.FisherMatrix_BNSWaveform_ThreeDetectors(sigpar, n_freq, shf_L, shf_V)
FisherMatrix = bilby_wv.FisherMatrix_ThreeDetectors(injection_parameters, start_time, minimum_frequency, interferometers)
# If I want to run the hmc on less than 9 dimensions, my 9x9 FIM will be really singular, hence I should remove the rows and columns of the parameters which I keep fixed to get a non-singular sub_FIM which I can invert, otherwise just run the last commented two lines.
ixgrid = np.ix_(conf.param_indices, conf.param_indices)
if len(conf.param_indices)<9:
sub_FIM = FisherMatrix[ixgrid]
inverse_sub_FIM = np.linalg.inv(sub_FIM)
sub_scale = np.sqrt(inverse_sub_FIM.diagonal())
scale = np.zeros(9)
scale[conf.param_indices] = sub_scale
else:
InverseFisherMatrix = np.linalg.inv(FisherMatrix)
scale = np.sqrt(InverseFisherMatrix.diagonal())
# re-adjust scales if FIM is incorrect
# if scale[0] >1.0: scale[0] = 1.0 # assume 50% error in cosinc
# if scale[1] > PI: scale[1] = PI # assume 50% error in phi_c
# if scale[2]>PIb2: scale[2] = PIb2 # assume 50% error in psi
# if scale[3]> 0.5: scale[3] = 0.5 # assume 50% error in lnDL
# classic decomp
###################
###################
# SVD decomp of FIM
# If I want to run the hmc on less than 9 dimensions, my 9x9 FIM will be really singular, hence I should remove the rows and columns of the parameters which I keep fixed to get a non-singular sub_FIM which I can invert, otherwise just run the last commented four lines.
if len(conf.param_indices)<9:
u, s, vh = np.linalg.svd(sub_FIM)
inverse_sub_FIM_SVD = vh.T * 1/s @ u.T
sub_scale_SVD = np.sqrt(inverse_sub_FIM_SVD.diagonal())
scale_SVD = np.zeros(9)
scale_SVD[conf.param_indices] = sub_scale_SVD
else:
# u * s @ vh = u @ np.diag(s) @ vh = FisherMatrix
u, s, vh = np.linalg.svd(FisherMatrix)
InverseFisherMatrix_SVD = vh.T * 1/s @ u.T
scale_SVD = np.sqrt(InverseFisherMatrix_SVD.diagonal())
# re-adjust scales if FIM is incorrect
# if scale_SVD[0] >1.0: scale_SVD[0] = 1.0 # assume 50% error in cosinc
# if scale_SVD[1] > PI: scale_SVD[1] = PI # assume 50% error in phi_c
# if scale_SVD[2]>PIb2: scale_SVD[2] = PIb2 # assume 50% error in psi
# if scale_SVD[3]> 0.5: scale_SVD[3] = 0.5 # assume 50% error in lnDL
# SVD decomp of FIM
###################
###################
# tf2_for_bilby FIM
sigpar, n = su.injection_parameters_to_sigpar(injection_parameters, minimum_frequency)
n_freq = len(interferometers[0].power_spectral_density_array)
FisherMatrix_hmc = tf2_for_bilby.FisherMatrix_ThreeDetectors(sigpar, n_freq, interferometers)
if len(conf.param_indices)<9:
sub_FIM_hmc = FisherMatrix_hmc[ixgrid]
inverse_sub_FIM_hmc = np.linalg.inv(sub_FIM_hmc)
sub_scale_hmc = np.sqrt(inverse_sub_FIM_hmc.diagonal())
scale_hmc = np.zeros(9)
scale_hmc[conf.param_indices] = sub_scale_hmc
else:
InverseFisherMatrix_hmc = np.linalg.inv(FisherMatrix_hmc)
scale_hmc = np.sqrt(InverseFisherMatrix_hmc.diagonal())
# re-adjust scales if FIM is incorrect
# if scale_hmc[0] >1.0: scale_hmc[0] = 1.0 # assume 50% error in cosinc
# if scale_hmc[1] > PI: scale_hmc[1] = PI # assume 50% error in phi_c
# if scale_hmc[2]>PIb2: scale_hmc[2] = PIb2 # assume 50% error in psi
# if scale_hmc[3]> 0.5: scale_hmc[3] = 0.5 # assume 50% error in lnDL
# tf2_for_bilby FIM
###################
# np.savetxt('FIM_compare.txt', FisherMatrix, fmt='%14.6e', header='FIM from Bilby')
# np.savetxt('FIM_Yann.txt', FisherMatrix_hmc, fmt='%14.6e', header='FIM from TaylorF2_Waveform_Threedetectors()')
# np.savetxt('FIM_ratio.txt', FisherMatrix/FisherMatrix_hmc, fmt='%14.6e', header='FIM_bilby/FIM_Yann')
# FISHER MATRICES #
######################################
# import IPython; IPython.embed()
return FisherMatrix, FisherMatrix_hmc, scale, scale_SVD, scale_hmc
if __name__=='__main__':
from optparse import OptionParser
usage = """%prog [options]
Plots the three waveforms in the time domain"""
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')
(opts,args)=parser.parse_args()
# SET THE INJECTION PARAMETERS
injection_parameters = set_inj.ini_file_to_dict()
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. But the interpolation is done automatically if `ifo.power_spectral_density_array` is called at some point.
set_psds.main(interferometers, opt_psd=opts.psd, sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
# COMPUTE THE FISHER MATRICES
FisherMatrix, FisherMatrix_hmc, scale, scale_SVD, scale_hmc = main(injection_parameters, start_time, minimum_frequency, interferometers)
# print("\nValues of FIM:")
# print("FIM from Bilby")
# for i in range(9):
# for j in range(9):
# print('{:12.4e}'.format(FisherMatrix[i][j]), end=' ')
# print('')
#
# print("FIM from TaylorF2_Waveform_Threedetectors()")
# for i in range(9):
# for j in range(9):
# print('{:12.4e}'.format(FisherMatrix_hmc[i][j]), end=' ')
# print('')
ndscale = np.zeros(9)
ndscale[0] = 1.705314e+00
ndscale[1] = 3.188156e-01
ndscale[2] = 1.027065e+02
ndscale[3] = 1.561083e+00
ndscale[4] = 2.813453e-05
ndscale[5] = 9.993015e-04
ndscale[6] = 3.007305e-03
ndscale[7] = 8.457038e-03
ndscale[8] = 1.403338e-06
print("\nValue of scales")
print("{:20} | {:20} | {:20} | {:20} | {:20} | {:20} | {:20} | {:20} | {:20}".format("param name", "bilby linalg.svd()", "hmc in python", "hmc in C", "hmc_py / hmc_C", "bilby / hmc_py", "bilby / hmc_C", "offset python", "param name"))
for i in range(9):
if scale_hmc[i]!=0:
ratio_bilby_hmc = scale_SVD[i]/scale_hmc[i]
else:
ratio_bilby_hmc = np.NaN
print("{:20} | {:20.6e} | {:20.6e} | {:20.6e} | {:20.2f} | {:20.2f} | {:20.2f} | {:20.0e} | {:20}\t".format(conf.fparam_names[i], scale_SVD[i], scale_hmc[i], ndscale[i], scale_hmc[i]/ndscale[i], ratio_bilby_hmc, scale_SVD[i]/ndscale[i], conf.param_offsets[i], conf.fparam_names[i]))
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import numpy as np
import bilby
import Library.bilby_waveform as bilby_wv
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.sigpar_utils as su
import Library.TaylorF2_for_bilby as tf2_for_bilby
import Library.config as conf
import Library.set_injection_parameters as set_inj
import Library.set_psds as set_psds
import Library.waveform_three_detectors as waveform_three_detectors
import Library.logL_snr as logL_snr
def main(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers, generate_from_bilby, do_prints=False):
if generate_from_bilby:
dlogL = bilby_wv.dlogL_ThreeDetectors(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers)
prefix = '\ndlogL with bilby '
else:
# Template generation with Yann's tf2
sigpar, n = su.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 '
if do_prints:
if conf.dlogL == 'dlogL2' and generate_from_bilby:
info = prefix + 'using: dlogL = logL(p+offset) - logL(p-offset)'
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]))
return dlogL
if __name__=='__main__':
from optparse import OptionParser
usage = """%prog [options]
Plots the three waveforms in the time domain"""
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("--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 = set_inj.ini_file_to_dict()
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
set_psds.main(interferometers, opt_psd=opts.psd, sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
# SET NOISE STRAIN DATA FROM PSD
interferometers.set_strain_data_from_power_spectral_densities(
sampling_frequency=sampling_frequency,
duration=duration,
start_time=start_time)
# WAVEFORM GENERATION
h_wave_1, h_wave_2, h_wave_3 = 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 = 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
dlogL = main(signal1, signal2, signal3, injection_parameters, start_time, minimum_frequency, interferometers, opts.generate_from_bilby, True)
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import numpy as np
import bilby
import Library.bilby_utils as bilby_utils
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.set_injection_parameters as set_inj
import Library.set_psds as set_psds
import Library.waveform_three_detectors as waveform_three_detectors
def main(h_wave_1, h_wave_2, h_wave_3, interferometers, do_prints):
logLBest, mf_snr_best, opt_snr_Best = bilby_utils.loglikelihood_snr(h_wave_1, h_wave_2, h_wave_3, interferometers)
if do_prints:
print("")
print("Optimal logL = {:.2f}".format(logLBest))
print("Matched filter SNR = {:.2f}".format(mf_snr_best))
print("Optimal SNR = {:.2f}".format(opt_snr_Best))
print("")
return logLBest, mf_snr_best, opt_snr_Best
if __name__=='__main__':
from optparse import OptionParser
usage = """%prog [options]
Plots the three waveforms in the time domain"""
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("--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 = set_inj.ini_file_to_dict()
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
set_psds.main(interferometers, opt_psd=opts.psd, sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
# SET NOISE STRAIN DATA FROM PSD
interferometers.set_strain_data_from_power_spectral_densities(
sampling_frequency=sampling_frequency,
duration=duration,
start_time=start_time)
# WAVEFORM GENERATION
h_wave_1, h_wave_2, h_wave_3 = 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 = main(h_wave_1, h_wave_2, h_wave_3, interferometers, True)
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
......@@ -97,3 +96,12 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
injection_parameters['meta'] = meta_params
return injection_parameters
if __name__=='__main__':
injection_parameters = ini_file_to_dict()
pu.print_dict(injection_parameters, indent=3, align_keys=True)
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import numpy as np
import bilby
import Library.psd_utils as psd_utils
import Library.set_injection_parameters as set_inj
def psd_file_used(interferometer):
psd_file = interferometer.power_spectral_density.psd_file
psd_file = psd_file.split(sep='/').pop()
return psd_file
def print_psd_file_used(interferometers):
for ifo in interferometers:
psd_file = psd_file_used(ifo)
print("The PSD file used for {} is: {}".format(ifo.name, psd_file))
def main(interferometers, opt_psd=1, sampling_frequency=None, duration=None, start_time=None):
# Those two methods are going to store psd into ifo.power_spectral_density.psd_array and the corresponding frequencies in ifo.power_spectral_density.frequency_array
if opt_psd==1:
GWTC1_GW170817_PSDs_file = '../.input_data/GW170817/GWTC1_GW170817_PSDs.dat'
GWTC1_GW170817_PSDs = np.loadtxt(GWTC1_GW170817_PSDs_file)
for i, ifo in enumerate(interferometers):
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(psd_array=GWTC1_GW170817_PSDs[:, i+1], frequency_array=GWTC1_GW170817_PSDs[:, 0])
elif opt_psd==2:
hdf5_file_path_H1 = '../.input_data/GW170817/H-H1_GWOSC_4KHZ_R1-1187006835-4096.hdf5'
hdf5_file_path_L1 = '../.input_data/GW170817/L-L1_GWOSC_4KHZ_R1-1187006835-4096.hdf5'
hdf5_file_path_V1 = '../.input_data/GW170817/V-V1_GWOSC_4KHZ_R1-1187006835-4096.hdf5'
hdf5_file_path_list = [hdf5_file_path_H1, hdf5_file_path_L1, hdf5_file_path_V1]
# offset in seconds added to the time of coalescence to get to start_time of the strain which will be used to estimate the psd. This strain_psd should not overlap with the gw signal into the strain data otherwise the snr will be less than real.
# psd_offset = -1024
psd_offset = -2048
# duration of the strain which will be used to estimate the psd. The longer this duration, the smoother the psd will be but that shouldn't change the snr much
# psd_duration = 100
psd_duration = 1800
filter_freq = None
# NFFT input description from `mlab.psd(.., NFFT = )`: The number of data points used in each block for the FFT. A power 2 is most efficient. The default value is 256. This should *NOT* be used to get zero padding, or the scaling of the result will be incorrect.
# NFFT description in gwosc script `lotsofplots.py`: number of sample for the fast fourier transform
NFFT = 4 * sampling_frequency
for i, ifo in enumerate(interferometers):
psd_utils.set_interferometer_psd_from_hdf5_file(interferometer=ifo, hdf5_file_path=hdf5_file_path_list[i], signal_start_time=start_time, signal_duration=duration, psd_duration=psd_duration, psd_offset=psd_offset, filter_freq=filter_freq, NFFT=NFFT)
elif opt_psd==3:
# Change LIGO psd files if wanted:
path = bilby.__path__[0] + '/gw/noise_curves/'
new_LIGO_psd_file_name = 'aLIGO_early_high_psd.txt'
# new_LIGO_psd_file_name = 'aLIGO_mid_psd.txt'
# new_LIGO_psd_file_name = 'aLIGO_ZERO_DET_high_P_psd.txt'
interferometers[0].power_spectral_density.psd_file = path + new_LIGO_psd_file_name
interferometers[1].power_spectral_density.psd_file = path + new_LIGO_psd_file_name
else:
print("Error: the psd option you have set must be either 1, 2 or 3. Now it is set to {}".format(opt_psd))
sys.exit(1)
print("")
if opt_psd==1:
print("PSDs estimated from the official PSDs published with GWTC1 which file is locally stored here: \n {}".format(GWTC1_GW170817_PSDs_file))
title = GWTC1_GW170817_PSDs_file.split('/').pop()
elif opt_psd==2:
print("PSDs estimated using the Welch methods from the following strain files:")
for i in range(len(interferometers)):
hdf5_file_name = hdf5_file_path_list[i].split(sep='/').pop()
print("- {}".format(hdf5_file_name))
title = 'Welch methods on strain data files'
elif opt_psd==3:
print_psd_file_used(interferometers)
title = psd_file_used(interferometers[0]) + ' and ' + psd_file_used(interferometers[2])
print("")
return title
if __name__=='__main__':
from optparse import OptionParser
usage = """%prog [options]
Plots the PSDs for each interferometers."""
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}")
(opts,args)=parser.parse_args()
# set_psds_from_open_data = False
# SET THE INJECTION PARAMETERS
injection_parameters = set_inj.ini_file_to_dict()
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