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

`pu.XX` to `su.XX`

parent a38f9cec
......@@ -11,6 +11,7 @@ from optparse import OptionParser,OptionValueError
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.BNS_Utility_Functions as buf
import Library.Python_Utility as pu
import Library.sigpar_utils as su
import Library.BNS_HMC_Tools as bht
import Library.LVC_Noise_Functions as lvnf
import Library.TaylorF2 as tf2
......@@ -266,7 +267,7 @@ if __name__=='__main__':
# Print starting quantities of the HMC
print("\n\n********** MAIN code *********\n")
q_pos = pu.injection_parameters_to_q_pos(injection_parameters, start_time)
q_pos = su.injection_parameters_to_q_pos(injection_parameters, start_time)
print("{:11} | {:20} | {:12} | {:20} | {:20} | {:20} | {:14}".format("param name", "param value (q_pos)", "fparam_name", "fparam value", "scale", "dlogL", "offset python"))
for i in range(9):
......@@ -338,7 +339,7 @@ if __name__=='__main__':
for l in range(1,n_traj_fit+1):
# generation of epsilon according to a Gaussian law with limitations between 0.001 and 0.01
# I import Yann's random numbers generated here since I can't use the same random number generator as his. This way I should be able to reproduce exactly the same chain and compare my results with his precisely.
# rannum = pu.get_next_random_numbers(1)
# rannum = su.get_next_random_numbers(1)
rannum = randGen.normal(loc=0, scale=1.5)
epsilon = (epsilon0 + rannum) * 1.0e-3
if epsilon <= 0.001: epsilon = 0.001
......@@ -511,7 +512,7 @@ if __name__=='__main__':
# Draw epsilon according to a Gaussian law limited between 0.001 and 0.01
# rannum = pu.get_next_random_numbers(1)
# rannum = su.get_next_random_numbers(1)
rannum = randGen.normal(loc=0, scale=1.5)
epsilon = (epsilon0 + rannum) * 1.0e-3
if epsilon <= 0.001: epsilon = 0.001
......@@ -530,7 +531,7 @@ if __name__=='__main__':
# Main HMC routines. Depending on the value of acceptance rate, we either use the full analytcal trajectories,
# the hybrid trajectories or the full numerical trajectories.
if Acc >= 0.65:
# lengthTana = 50 + int(pu.get_next_random_numbers(1))
# lengthTana = 50 + int(su.get_next_random_numbers(1))
lengthTana = 50 + int(randGen.uniform() * 100)
# QUESTION: flagNumTable seems to be not used anymore
q_pos, logL, Snr, dlogL, Accept = bht.AnalyticalGradientLeapfrog_ThreeDetectors(q_pos, logL, Snr, dlogL, signal1, signal2, signal3, shf_L, shf_V, randGen, n, epsilon, lengthTana, scale, fit_coeff_good, PointFit_SortCosInc, PointFit_SortPsi, PointFit_SortLogD, boundary, n_param, Accept,flagNumTable, n_pt_fit_bad, n_fit_1, n_fit_2,SquaredScale)
......
......@@ -7,7 +7,8 @@ import numpy as np
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.BNS_Utility_Functions as buf
import Library.Python_Utility as pu
import Library.Python_Utility as pu
import Library.sigpar_utils as su
import Library.TaylorF2 as tf2
import Library.TaylorF2_for_bilby as tf2_for_bilby
# import Library.SWIG_TF2.SWIG_files.Swig_TaylorF2 as stf2
......@@ -233,7 +234,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, si
"""
# Initialization of the momentum: draw from a gaussian distribution with mean 0 and sigma=1
# p_mom_0 = pu.get_next_random_numbers(9)
# p_mom_0 = su.get_next_random_numbers(9)
p_mom_0 = randGen.normal(loc=0, scale=1.0, size=9)
# Main HMC routine with numerical gradients
......@@ -245,7 +246,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, si
# Compute waveforms
# h_wave_1, h_wave_2, h_wave_3 = tf2.TaylorF2_Waveform_Threedetectors(q_pos, n_freq)
# 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 = pu.q_pos_to_dictionary(q_pos, start_time)
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)
......@@ -254,7 +255,7 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, si
logL_final = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
H = computeHamiltonian(p_mom,logL_final)
# a = pu.get_next_random_numbers(1)
# a = su.get_next_random_numbers(1)
# a = randGen.normal(loc=0, scale=1.0) # draw random number
a = randGen.uniform(low=0, high=1) # draw random number
......@@ -360,7 +361,7 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
DlogLFit_traj[0] = dlogL_pos
if conf.debug:
parameters = pu.q_pos_to_dictionary(q_pos, start_time)
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)
logL = bilby_utils.loglikelihood(h_wave_1, h_wave_2, h_wave_3, interferometers)
......@@ -473,7 +474,7 @@ def NumericalGradientLeapfrog_ThreeDetectors_Trajectory(q_pos_0, p_mom_0, dlogL_
# dlogL_pos = tf2.dlogL_ThreeDetectors(signal_1, signal_2, signal_3, shf_L, shf_V, q_pos, n_freq)
# dlogL_pos = stf2.dlogL_ThreeDetectors(9, signal_1, signal_2, signal_3, shf_L, shf_V, q_pos, n_freq)
# 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 = pu.q_pos_to_dictionary(q_pos, start_time)
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 = tf2_for_bilby.dlogL_ThreeDetectors(signal_1, signal_2, signal_3, q_pos, n_freq, interferometers)
......@@ -519,7 +520,7 @@ def AnalyticalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, s
n_freq = int(n_time/2)
# Initialization of the momentum: draw from a gaussian distribution with mean 0 and sigma=1
# p_mom_0 = pu.get_next_random_numbers(9)
# p_mom_0 = su.get_next_random_numbers(9)
p_mom_0 = randGen.normal(loc=0, scale=1.0, size=9)
# Main HMC routine with analytical gradients
......@@ -551,7 +552,7 @@ def AnalyticalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, s
# 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)
H = computeHamiltonian(p_mom,logL_final)
# a = pu.get_next_random_numbers(1)
# a = su.get_next_random_numbers(1)
a = randGen.normal(loc=0, scale=1.0) # draw random number
# Hamiltonian at the start of the trajectory
......
......@@ -4,15 +4,6 @@ import sys
sys.path.append('../')
import os
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.BNS_Utility_Functions as buf
from decimal import Decimal
import time
import bilby.core.utils as bc_utils
import bilby.gw.utils as gw_utils
import Library.config as conf
import Library.set_injection_parameters as set_inj
def read_dat_file(fileName,columnSep=None):
"""
......@@ -150,335 +141,8 @@ def config_parser_to_dict(config):
return dictionary
# def as_dict(config):
# d = dict(config._sections)
# for k in d:
# d[k] = dict(config._defaults, **d[k])
# d[k].pop('__name__', None)
# return d
def binaryFile_to_sigpar(binary_file_name):
"""
Reads a binary.in file and returns the sigpar array and n
"""
params = np.loadtxt(binary_file_name)
(RA, Dec, inc, psi, phic, dist, m1, m2, GMST, f_low, Rmin, m1min, m1max, m2min, m2max, Dmin, Dmax)=params
sigpar, n = buf.CBC_Parameter_Allocation(m1, m2, Dec, RA, inc, psi, phic, dist, GMST, Rmin, m1min+m2min, f_low)
return sigpar, n
def binaryFile_to_boundary(binary_file_name):
"""
Reads a binary.in file and returns the boundary conditions for MCMC
"""
params = np.loadtxt(binary_file_name)
(RA, Dec, inc, psi, phic, dist, m1, m2, GMST, f_low, Rmin, m1min, m1max, m2min, m2max, Dmin, Dmax)=params
boundary=np.array([[m1min,m1max],[m2min,m2max],[Dmin * MPC,Dmax * MPC]])
return boundary
def sigpar_to_injection_parameters(sigpar, start_time):
"""
From the hmc sigpar array of parameters, returns the injection_parameters needed by bilby to run.
"""
# bilby.gw.utils.time_delay_geocentric(H.vertex, np.array([0,0,0]), injection_parameters['ra'], injection_parameters['dec'], start_time)
# time_shift_H1 = 0.012139641777406887
time_shift_H1 = gw_utils.time_delay_geocentric(conf.H1_vertex, np.array([0,0,0]), sigpar[7], sigpar[6], start_time)
F_plus_H1 = conf.H1.antenna_response(sigpar[7], sigpar[6], start_time + sigpar[8] - time_shift_H1, sigpar[2], 'plus')
F_cross_H1 = conf.H1.antenna_response(sigpar[7], sigpar[6], start_time + sigpar[8] - time_shift_H1, sigpar[2], 'cross')
phase_shift_H1 = np.angle((0.5*(1+np.cos(sigpar[0])**2)*F_plus_H1) - 1j * np.cos(sigpar[0])*F_cross_H1)
# if phase_shift_H1 < 0: phase_shift_H1 += 2 * np.pi
injection_parameters = dict(
mass_1=sigpar[9],
mass_2=sigpar[10],
chi_1=0.0,
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)),
phase=0.5 * (sigpar[1]),
# geocent_time=geocent_time,
geocent_time=start_time + sigpar[8] - time_shift_H1,#geocent_time,
# geocent_time=start_time + sigpar[8],#geocent_time,
ra=sigpar[7],
dec=sigpar[6],
lambda_1=0.,
lambda_2=0.
)
# print("\n2 * phase = {} * pi".format(2 * injection_parameters['phase']/np.pi))
# import IPython; IPython.embed();
return injection_parameters
def injection_parameters_to_sigpar(injection_parameters, minimum_frequency):
# time_shift_H1 = gw_utils.time_delay_geocentric(conf.H1_vertex, np.array([0,0,0]), sigpar[7], sigpar[6], start_time)
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 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']
GMST = bc_utils.gps_time_to_gmst(injection_parameters['geocent_time']) # is GMST !=0 possible in Yann's code??...
Rmin = 6.
m1min = 1.
m2min = 1.
f_low = minimum_frequency
sigpar, n = buf.CBC_Parameter_Allocation(m1, m2, Dec, RA, inc, psi, phic, dist, GMST, Rmin, m1min+m2min, f_low)
return sigpar, n
def injection_parameters_to_binaryFile(injection_parameters, minimum_frequency, binaryFileName):
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 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']
GMST = bc_utils.gps_time_to_gmst(injection_parameters['geocent_time']) # is GMST !=0 possible in Yann's code??...
Rmin = 6.
m1min = 1.
m2min = 1.
f_low = minimum_frequency
m1max = 2.6
m2max = 2.6
Dmin = 1e-6
Dmax = 70
params = (RA, Dec, inc, psi, phic, dist, m1, m2, GMST, f_low, Rmin, m1min, m1max, m2min, m2max, Dmin, Dmax)
np.savetxt(binaryFileName, params)
return 0
def injection_parameters_to_q_pos(injection_parameters, start_time):
q_pos = np.zeros(13)
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']
q_pos[3] = injection_parameters['luminosity_distance'] * MPC
q_pos[4] = injection_parameters['chirp_mass']
q_pos[5] = injection_parameters['reduced_mass']
q_pos[6] = injection_parameters['dec']
q_pos[7] = injection_parameters['ra']
q_pos[8] = injection_parameters['geocent_time']-start_time
q_pos[9] = injection_parameters['mass_1']
q_pos[10] = injection_parameters['mass_2']
q_pos[11] = injection_parameters['total_mass']
q_pos[12] = injection_parameters['symmetric_mass_ratio']
return q_pos
def q_pos_to_dictionary_original(q_pos, start_time):
"""
True version of this conversion function. Using in trajectories the version bellow so as to run effectively the hmc over a subset of parameters, ie dimensions.
The point of this is to target which parameter(s) mess-up the acceptance rate.
"""
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
dictionary['chirp_mass'] = q_pos[4]
dictionary['reduced_mass'] = q_pos[5]
dictionary['dec'] = q_pos[6]
dictionary['ra'] = q_pos[7]
dictionary['geocent_time'] = q_pos[8] + start_time
dictionary['mass_1'] = q_pos[9]
dictionary['mass_2'] = q_pos[10]
dictionary['total_mass'] = q_pos[11]
dictionary['symmetric_mass_ratio'] = q_pos[12]
dictionary['chi_1'] = 0 # bilby raises an error of this key is absent
dictionary['chi_2'] = 0 # bilby raises an error of this key is absent
return dictionary
def q_pos_to_dictionary(q_pos, start_time, dictionary=None):
if dictionary is 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:
dictionary['psi'] = q_pos[2]
if 3 in conf.param_indices:
dictionary['luminosity_distance'] = q_pos[3] / MPC
if 4 in conf.param_indices:
dictionary['chirp_mass'] = q_pos[4]
dictionary['mass_1'] = q_pos[9]
dictionary['mass_2'] = q_pos[10]
dictionary['total_mass'] = q_pos[11]
dictionary['symmetric_mass_ratio'] = q_pos[12]
if 5 in conf.param_indices:
dictionary['reduced_mass'] = q_pos[5]
dictionary['mass_1'] = q_pos[9]
dictionary['mass_2'] = q_pos[10]
dictionary['total_mass'] = q_pos[11]
dictionary['symmetric_mass_ratio'] = q_pos[12]
if 6 in conf.param_indices:
dictionary['dec'] = q_pos[6]
if 7 in conf.param_indices:
dictionary['ra'] = q_pos[7]
if 8 in conf.param_indices:
dictionary['geocent_time'] = q_pos[8] + start_time
dictionary['chi_1'] = 0 # bilby raises an error of this key is absent
dictionary['chi_2'] = 0 # bilby raises an error of this key is absent
return dictionary
def sigpar_to_ifo_parameters(sigpar):
"""
"""
duration = sigpar[16] # = Tobs = 1.1 * tchirp for us, cf buf.CBC_Parameter_Allocation()
sampling_frequency = 2.0 * sigpar[15]
# I don't use the + 2 here because it's kind of already done using Tobs = 1.1 * tchirp. Note sigpar[8] = tchirp = sigpar[16] / 1.1
# In bilby's example, the start_time is computed using the formula below where 'geocent_time' is the time of coalescence, ie our sigpar[8], but expressed in gps time when in the hmc sigpar[8] is counted from start_time
# start_time = injection_parameters['geocent_time'] + 2 - duration
# So we derive the start_time from our input GMST angle = sipar[18]. This will be the reference start_time through the whole algorithm, the geocent_time however will vary when sigpar[8] varies.
# However sigpar[18] is the GMST angle of the earth when the waveform arrives in Handford, what bilby wants is the start_time = time when the GW arrives at the center of the earth. Hence I should input my function `gmst_to_gps_time()` not with sigpar[18] but the GMST when the GW arrived at the earth center.
# bilby.gw.utils.time_delay_geocentric(H.vertex, np.array([0,0,0]), injection_parameters['ra'], injection_parameters['dec'], start_time)
H_time_delay_from_geocenter = 0.012139641777406887
omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
gmst_earth_center = sigpar[18] - H_time_delay_from_geocenter * omega_earth
start_time = gmst_to_gps_time(sigpar[18])
# start_time = geocent_time - sigpar[8]
# start_time = 1541967481.3621483#geocent_time - sigpar[8]
return duration, sampling_frequency, start_time
def gmst_to_gps_time(gmst):
"""
Converts a gmst angle to the corresponding gps_time closest to that of now.
"""
gps_time_now = time.time()
gmst_now = bc_utils.gps_time_to_gmst(gps_time_now)
gmst_diff = gmst_now - gmst
omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
number_of_seconds = gmst_diff / omega_earth
gps_time = gps_time_now - number_of_seconds
return gps_time
def sigpar_to_freq_and_indices(sigpar,n_freq):
"""
When computing waveforms and/or their derivatives, Fourier Transform purposes imposes that the array of frequencies has to start from zero and have a power of 2 length=n_freq. Hence the end frequency is higher than f_high.
We then set to zero all the waveform values corresponding to frequencies f not respecting : f_low <= f <= f_high; where f_low is an input of the binary.in file and f_high calculated from it.
This function returns:
return freqs, i_low, i_high
Where:
- freqs [array of floats]:
array of all the frequencies from 0 to n_freq*(delta_f-1) included
- i_low [int]:
s.t. freqs[i_low] is the 1st frequency where the waveform is non-zero
- i_high [int]:
s.t. freqs[i_high] is the last frequency where the waveform is non-zero
"""
delta_f = sigpar[17]
f_low = sigpar[13]
f_high = sigpar[14]
freqs = np.arange(0,n_freq*delta_f,delta_f)
i_low=int(f_low/delta_f)+1
i_high=int(f_high/delta_f)
return freqs, i_low, i_high
def sigpar_to_freq_and_indices_2(sigpar,n_freq):
"""
When computing waveforms and/or their derivatives, Fourier Transform purposes imposes that the array of frequencies has to start from zero and have a power of 2 length=n_freq. Hence the end frequency is higher than f_high.
We then set to zero all the waveform values corresponding to frequencies f not respecting : f_low <= f <= f_high; where f_low is an input of the binary.in file and f_high calculated from it.
This function returns:
return freqs, i_low, i_high
Where:
- freqs [array of floats]:
array of all the frequencies from 0 to n_freq*(delta_f-1) included
- i_low [int]:
s.t. freqs[i_low] is the 1st frequency where the waveform is non-zero
- i_high [int]:
s.t. freqs[i_high] is the last frequency where the waveform is non-zero
"""
duration, sampling_frequency, start_time = sigpar_to_ifo_parameters(sigpar)
import Library.set_injection_parameters as set_inj
duration = set_inj.injection_parameters['meta']['duration']
sampling_frequency = set_inj.injection_parameters['meta']['sampling_frequency']
freqs = bc_utils.create_frequency_series(sampling_frequency, duration)
f_low = sigpar[13]
f_high = sigpar[14]
delta_f = freqs[1] - freqs[0]
i_low=int(f_low/delta_f)+1
i_high=int(f_high/delta_f)
return freqs, i_low, i_high
def sigpar_to_n_freq_as_computed_in_bilby(sigpar):
"""
This comes from the computation in bilby.core.utils.create_frequency_series()
"""
duration, sampling_frequency, start_time = sigpar_to_ifo_parameters(sigpar)
number_of_samples = duration * sampling_frequency
number_of_samples = int(np.round(number_of_samples))
# prepare for FFT
number_of_frequencies = (number_of_samples - 1) // 2
if number_of_frequencies % 2 == 1:
number_of_frequencies += 2
else:
# no Nyquist frequency when N=odd
number_of_frequencies += 1
return number_of_frequencies
def scalarProduct(h1,h2,step=1):
"""
......@@ -507,7 +171,6 @@ def scalarProduct(h1,h2,step=1):
return arg
def overlap(h1,h2):
"""
return scalarProduct(h1,h2,1)/sqrt(scalarProduct(h1,h1,1)*scalarProduct(h2,h2,1))
......@@ -534,22 +197,3 @@ def overlap(h1,h2):
# h1_normalised = h1/h1_h1_norm
# h2_normalised = h2/h2_h2_norm
# return scalarProduct(h1_normalised,h2_normalised,1)
def overlap_LVC_RiemannSum(CA1, CA2, shf, Tobs, n, flow, fhigh):
d1=buf.LVC_RiemannSum(CA1, CA1, shf, Tobs, n, flow, fhigh)
d2=buf.LVC_RiemannSum(CA2, CA2, shf, Tobs, n, flow, fhigh)
num=buf.LVC_RiemannSum(CA1, CA2, shf, Tobs, n, flow, fhigh)
return num/np.sqrt(d1*d2)
# random_numbers = np.loadtxt('../.input_data/gsl_ran_numbers.dat')
import Library.config as conf
def get_next_random_numbers(n):
# global random_numbers
next_random_numbers = conf.random_numbers[0:n]
conf.random_numbers = conf.random_numbers[n:]
if n==1: next_random_numbers = next_random_numbers[0]
return next_random_numbers
......@@ -11,6 +11,7 @@ from Headers.PN_Coefficients import * # Headers.Constants already imported in PN
import Library.BNS_Utility_Functions as buf
import Library.LVC_Detector_Functions as lvdf
import Library.Python_Utility as pu
import Library.sigpar_utils as su
def ZipMerge(h_real,h_img):
"""
......@@ -144,7 +145,7 @@ def TaylorF2_Waveform_Threedetectors(sigpar, n_freq):
# Replacing the loop over the range of frequency with np.array manipulations
# Retrieve the frequencies f s.t. : f_low <= f <= f_high, for which we are not going to set the waveform values to zero
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices(sigpar, n_freq)
freq = all_freqs[i_low_num:i_high_num+1]
AmpWave = freq**c76
......@@ -268,7 +269,7 @@ def TaylorF2_Waveform_Threedetectors_complex(sigpar, n_freq):
# Replacing the loop over the range of frequency with np.array manipulations
# Retrieve the frequencies f s.t. : f_low <= f <= f_high, for which we are not going to set the waveform values to zero
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices_2(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices_2(sigpar, n_freq)
freq = all_freqs[i_low_num:i_high_num+1]
AmpWave = freq**c76
......
......@@ -11,6 +11,7 @@ from Headers.PN_Coefficients import * # Headers.Constants already imported in PN
import Library.BNS_Utility_Functions as buf
import Library.LVC_Detector_Functions as lvdf
import Library.Python_Utility as pu
import Library.sigpar_utils as su
import Library.TaylorF2 as tf2
import Library.bilby_utils as bilby_utils
import Library.config as conf
......@@ -196,7 +197,7 @@ def TaylorF2_h_dh_Generation_ThreeDetectors(sigpar, n_freq):
offset_lntc = conf.param_offsets[8]
if offset_lntc==0:
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices_2(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices_2(sigpar, n_freq)
two_pi_f_tc = -1j * 2*np.pi * all_freqs * sigpar[8]
dh_Wave_1[8] = two_pi_f_tc * h_wave_1
dh_Wave_2[8] = two_pi_f_tc * h_wave_2
......
......@@ -2,6 +2,7 @@ import bilby
import numpy as np
import Library.TaylorF2 as tf2
import Library.Python_Utility as pu
import Library.sigpar_utils as su
# Stuff coming from the binary_neutron_star_example.py
def WaveForm_ThreeDetectors(sigpar, n_freq):
......@@ -76,7 +77,7 @@ def WaveForm_ThreeDetectors(sigpar, n_freq):
# Those signal do not spand on the same frequency band as for us, bilby does not go up until a power of 2.
# Hence we are going to expand and padd with zeros those signals at the end.
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices(sigpar, n_freq)
freq = all_freqs[i_low_num:i_high_num+1]
n1=0
......@@ -169,7 +170,7 @@ def WaveForm_ThreeDetectors_with_noise(sigpar, n_freq):
# Those signal do not spand on the same frequency band as for us, bilby does not go up until a power of 2.
# Hence we are going to expand and padd with zeros those signals at the end.
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices(sigpar, n_freq)
freq = all_freqs[i_low_num:i_high_num+1]
n1=0
......
......@@ -6,12 +6,13 @@ import numpy as np
import bilby.gw.utils as gwutils
import Library.Python_Utility as pu
import Library.sigpar_utils as su
import bilby
def log_likelihood_from_gravitational_wave_transient(sigpar, start_time, interferometers):
# duration, sampling_frequency, start_time = pu.sigpar_to_ifo_parameters(sigpar)
injection_parameters = pu.sigpar_to_injection_parameters(sigpar, start_time)
# duration, sampling_frequency, start_time = su.sigpar_to_ifo_parameters(sigpar)
injection_parameters = su.sigpar_to_injection_parameters(sigpar, start_time)
# for interferometer in interferometers:
# interferometer.duration = duration
# interferometer.start_time = start_time
......
......@@ -3,6 +3,7 @@ import bilby
import sys
import Library.bilby_utils as bilby_utils
import Library.Python_Utility as pu
import Library.sigpar_utils as su
import Library.BNS_Utility_Functions as buf
import Library.bilby_detector as bilby_det
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
......@@ -605,8 +606,8 @@ def dlogL_ThreeDetectors2(signal1, signal2, signal3, parameters, start_time, min
def WaveForm_ThreeDetectors_sigpar(sigpar, start_time, interferometers):
# duration, sampling_frequency, start_time = pu.sigpar_to_ifo_parameters(sigpar)
injection_parameters = pu.sigpar_to_injection_parameters(sigpar, start_time)
# duration, sampling_frequency, start_time = su.sigpar_to_ifo_parameters(sigpar)
injection_parameters = su.sigpar_to_injection_parameters(sigpar, start_time)
# for interferometer in interferometers:
# interferometer.duration = duration
# interferometer.start_time = start_time
......@@ -1209,7 +1210,7 @@ def WaveForm_ThreeDetectors_padded(sigpar, n_freq, start_time, interferometers):
return signal_H, signal_L, signal_V
def padd_with_zeros(signal_H, signal_L, signal_V, sigpar, n_freq):
all_freqs, i_low_num, i_high_num = pu.sigpar_to_freq_and_indices(sigpar, n_freq)
all_freqs, i_low_num, i_high_num = su.sigpar_to_freq_and_indices(sigpar, n_freq)