Commit 0ac0ce10 authored by Marc Arene's avatar Marc Arene
Browse files

- Correcting: a = randGen.normal() to a = randGen.uniform() !!

- 1D-HMC with distance
parent 2ea675f4
......@@ -22,6 +22,7 @@ import Library.config as conf
import bilby
import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
import Library.set_injection_parameters as set_inj
if __name__=='__main__':
......@@ -60,18 +61,25 @@ if __name__=='__main__':
print("The corresponding random number file: {} doesn't exist.".format(gsl_ran_numbers_filename))
# Create a subfolder to store output files
output_dir = "../.output_data/GW170817" + "_bilbyoriented" + "/"
# sys.argv[6] = sub directory to create
if sys.argv[6] == 'None':
sub_dir = ''
else:
sub_dir = sys.argv[6] + '/'
output_dir = "../.output_data/"+ sys.argv[1] + "/bilbyoriented/" + sub_dir
sub_output_path = "{}{:.0e}_{}_{}_{}_{}_{}/".format(output_dir, n_traj_HMC_tot, n_traj_fit, n_fit_1, n_fit_2, lengthTnum, epsilon0/1000)
# sub_output_path += '_bilbyoriented/'
try:
subprocess.check_call(["mkdir", output_dir])
print("Creating directory: {}".format(output_dir))
except subprocess.CalledProcessError:
try:
subprocess.check_call(["mkdir", sub_output_path])
print("Creating directory: {}".format(sub_output_path))
except subprocess.CalledProcessError:
pass
pu.mkdirs(sub_output_path)
# try:
# subprocess.check_call(["mkdir", output_dir])
# print("Creating directory: {}".format(output_dir))
# except subprocess.CalledProcessError:
# try:
# subprocess.check_call(["mkdir", sub_output_path])
# print("Creating directory: {}".format(sub_output_path))
# except subprocess.CalledProcessError:
# pass
# sub_output_path += 'DhReal_FIMbilby/'
# Create output files: chain file, info file, file containing the coefficients of global cubic fit, file containing the positions and gradients of points used for fits
......@@ -86,37 +94,38 @@ if __name__=='__main__':
#################################################
randGen = np.random.RandomState(seed=2)
injection_parameters = dict(
mass_1=1.47,
mass_2=1.27,
chi_1=0.0,
chi_2=0.0,
luminosity_distance=40,
# iota=146*np.pi/180,
iota=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.
)
injection_parameters = bilby.gw.conversion.generate_mass_parameters(injection_parameters)
injection_parameters['reduced_mass'] = (injection_parameters['mass_1'] * injection_parameters['mass_2']) / (injection_parameters['mass_1'] + injection_parameters['mass_2'])
injection_parameters = set_inj.injection_parameters
# injection_parameters = dict(
# mass_1=1.47,
# mass_2=1.27,
# chi_1=0.0,
# chi_2=0.0,
# luminosity_distance=40,
# # iota=146*np.pi/180,
# iota=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.
# )
#
# injection_parameters = bilby.gw.conversion.generate_mass_parameters(injection_parameters)
# injection_parameters['reduced_mass'] = (injection_parameters['mass_1'] * injection_parameters['mass_2']) / (injection_parameters['mass_1'] + injection_parameters['mass_2'])
#
#
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into. For the
# TaylorF2 waveform, we cut the signal close to the isco frequency
minimum_frequency=30.0
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
tc_3p5PN = bilby_wv.ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
tc_3p5PN = set_inj.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
f_high = bilby_wv.Compute_fmax(R_isco, m1, m2)
f_high = set_inj.Compute_fmax(R_isco, m1, m2)
# import IPython; IPython.embed();
duration = tc_3p5PN + 2
sampling_frequency = 4096
......@@ -279,7 +288,7 @@ if __name__=='__main__':
FisherMatrix = bilby_wv.FisherMatrix_ThreeDetectors(injection_parameters, start_time, minimum_frequency, interferometers)
n_freq = len(interferometers.frequency_array)
FIM_from_hmc = True
FIM_from_hmc = False
if FIM_from_hmc:
###################
# tf2_for_bilby FIM
......@@ -289,13 +298,20 @@ if __name__=='__main__':
# tf2_for_bilby 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
ixgrid = np.ix_(conf.param_indices, conf.param_indices)
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
# if optionScale == 1:
# InverseFisherMatrix = np.linalg.inv(FisherMatrix)
# scale = np.sqrt(InverseFisherMatrix.diagonal())
# else:
# scale = np.sqrt(1 / FisherMatrix.diagonal())
if optionScale == 1:
InverseFisherMatrix = np.linalg.inv(FisherMatrix)
scale = np.sqrt(InverseFisherMatrix.diagonal())
else:
scale = np.sqrt(1 / FisherMatrix.diagonal())
# re-adjust scales if FIM is incorrect
if scale[0] >1.0: scale[0] = 1.0 # assume 50% error in cosinc
......
......@@ -12,6 +12,7 @@ import Library.TaylorF2_for_bilby as tf2_for_bilby
import Library.SWIG_TF2.SWIG_files.Swig_TaylorF2 as stf2
import Library.bilby_waveform as bilby_wv
import Library.bilby_utils as bilby_utils
import Library.config as conf
def FindClosestPoints(q_pos_fit, PointFit_Sort, n_pt_fit, n_param, n_1_fit, param_index):
"""
......@@ -248,11 +249,19 @@ def NumericalGradientLeapfrog_ThreeDetectors(q_pos_0, logL_0, Snr_0, dlogL_0, si
# 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)
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 = randGen.normal(loc=0, scale=1.0) # draw random number
# a = randGen.normal(loc=0, scale=1.0) # draw random number
a = randGen.uniform(low=0, high=1) # draw random number
# Hamiltonian at the start of the trajectory
H_0 = computeHamiltonian(p_mom_0, logL_0)
# print("logL_0 = {}".format(logL_0))
# print("H_0 = {}".format(H_0))
# print("logL_final = {}".format(logL_final))
# print("H = {}".format(H))
# import IPython; IPython.embed();
# sys_exit = [False]
# import IPython; IPython.embed();
# if sys_exit[0]: sys.exit()
......@@ -621,12 +630,18 @@ def FindClosestPointParameter(Parameter, SortedVectorParameter, n_pt_fit):
return IndexClosest
def computeHamiltonian(p_mom, logL):
def computeHamiltonian_original(p_mom, logL):
"""
Routine that computes the Hamiltonian given the momenta and the log-likelihood
"""
return 0.5 * (p_mom ** 2).sum() - logL
def computeHamiltonian(p_mom, logL):
"""
Routine that computes the Hamiltonian given the momenta and the log-likelihood
"""
return 0.5 * (p_mom[conf.param_indices] ** 2).sum() - logL
def BoundaryCheck_On_Position(q_pos, p_mom, boundary, epsilon, scale, cos_inc, sin_theta):
"""
......
......@@ -2,6 +2,7 @@ import numpy as np
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import os
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.BNS_Utility_Functions as buf
......@@ -10,6 +11,8 @@ 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):
"""
......@@ -49,6 +52,15 @@ def write_dat_file(dataArray,fileName,columnSep='\t',floatPrecision=12):
outputFile.write("{:.{}e}\n".format(dataArray[i][-1],floatPrecision,columnSep))
outputFile.closed
def mkdirs(path):
"""
Helper function. Make the given directory, creating intermediate
dirs if necessary, and don't complain about it already existing.
"""
if os.access(path,os.W_OK) and os.path.isdir(path): return
else: os.makedirs(path)
def binaryFile_to_sigpar(binary_file_name):
"""
Reads a binary.in file and returns the sigpar array and n
......@@ -178,7 +190,11 @@ def injection_parameters_to_q_pos(injection_parameters, start_time):
return q_pos
def q_pos_to_dictionnary(q_pos, start_time):
def q_pos_to_dictionnary_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.
"""
dictionnary = dict()
dictionnary['iota'] = q_pos[0]
dictionnary['phase'] = q_pos[1] * 0.5
......@@ -198,6 +214,26 @@ def q_pos_to_dictionnary(q_pos, start_time):
return dictionnary
def q_pos_to_dictionnary(q_pos, start_time):
dictionnary = set_inj.injection_parameters
# dictionnary['iota'] = q_pos[0]
# dictionnary['phase'] = q_pos[1] * 0.5
# dictionnary['psi'] = q_pos[2]
dictionnary['luminosity_distance'] = q_pos[3] / MPC
# dictionnary['chirp_mass'] = q_pos[4]
# dictionnary['reduced_mass'] = q_pos[5]
# dictionnary['dec'] = q_pos[6]
# dictionnary['ra'] = q_pos[7]
# dictionnary['geocent_time'] = q_pos[8] + start_time
# dictionnary['mass_1'] = q_pos[9]
# dictionnary['mass_2'] = q_pos[10]
# dictionnary['total_mass'] = q_pos[11]
# dictionnary['symmetric_mass_ratio'] = q_pos[12]
# dictionnary['chi_1'] = 0 # bilby raises an error of this key is absent
# dictionnary['chi_2'] = 0 # bilby raises an error of this key is absent
return dictionnary
def sigpar_to_ifo_parameters(sigpar):
"""
......
......@@ -6,7 +6,7 @@ import Library.Python_Utility as pu
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
# Stuff coming from the binary_neutron_star_example.py
import Library.config as conf
def WaveForm_ThreeDetectors(parameters, start_time, minimum_frequency, interferometers):
......@@ -73,169 +73,173 @@ def h_dh_Generation_ThreeDetectors(parameters, start_time, minimum_frequency, in
# dh_Wave_2 = np.empty((number_of_parameters,waveform_length))
# dh_Wave_3 = np.empty((number_of_parameters,waveform_length))
# Deriving w.r.t. np.cos(inc) where inc = parameters[0]
parameters_minus = parameters.copy()
parameters_minus['iota'] = np.arccos(np.cos(parameters['iota'])-offset)
if 0 in conf.param_indices:
# Deriving w.r.t. np.cos(inc) where inc = parameters[0]
parameters_minus = parameters.copy()
parameters_minus['iota'] = np.arccos(np.cos(parameters['iota'])-offset)
parameters_plus = parameters.copy()
parameters_plus['iota'] = np.arccos(np.cos(parameters['iota'])+offset)
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[0] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[0] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[0] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. phi_c = parameters['phase'] * 2
# h(f) = A * np.exp()
# dh/dphic = (h(phic+offset)-h(phic-offset))/(2*offset)
# = (h(phi_c)*np.exp(i*offset)-h(phi_c)*np.exp(-i*offset))/2*offset # if -1j convention used
# = i*sin(offset)*h(phi_ref)/offset
# = i * h # since offset is really small
dh_Wave_1[1] = 1j * h_wave_1
dh_Wave_2[1] = 1j * h_wave_2
dh_Wave_3[1] = 1j * h_wave_3
# parameters_minus = parameters.copy()
# parameters_minus['phase'] = parameters['phase']-offset/2
#
# parameters_plus = parameters.copy()
# parameters_plus['phase'] = parameters['phase']+offset/2
#
# h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
# h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
#
# dh_Wave_1[1] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
# dh_Wave_2[1] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
# dh_Wave_3[1] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
parameters_plus = parameters.copy()
parameters_plus['iota'] = np.arccos(np.cos(parameters['iota'])+offset)
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[0] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[0] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[0] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
if 1 in conf.param_indices:
# Deriving w.r.t. phi_c = parameters['phase'] * 2
# h(f) = A * np.exp()
# dh/dphic = (h(phic+offset)-h(phic-offset))/(2*offset)
# = (h(phi_c)*np.exp(i*offset)-h(phi_c)*np.exp(-i*offset))/2*offset # if -1j convention used
# = i*sin(offset)*h(phi_ref)/offset
# = i * h # since offset is really small
dh_Wave_1[1] = 1j * h_wave_1
dh_Wave_2[1] = 1j * h_wave_2
dh_Wave_3[1] = 1j * h_wave_3
# parameters_minus = parameters.copy()
# parameters_minus['phase'] = parameters['phase']-offset/2
#
# parameters_plus = parameters.copy()
# parameters_plus['phase'] = parameters['phase']+offset/2
#
# h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
# h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
#
# dh_Wave_1[1] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
# dh_Wave_2[1] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
# dh_Wave_3[1] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
if 2 in conf.param_indices:
# Deriving w.r.t. psi = parameters['psi']
parameters_minus = parameters.copy()
parameters_minus['psi'] = parameters['psi']-offset
parameters_plus = parameters.copy()
parameters_plus['psi'] = parameters['psi']+offset
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[2] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[2] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[2] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. psi = parameters['psi']
parameters_minus = parameters.copy()
parameters_minus['psi'] = parameters['psi']-offset
if 3 in conf.param_indices:
# Deriving w.r.t. ln(D) where D = parameters['luminosity_distance']
# Here since h ~ constant/D, dh/dD = -h/D. Hence dh/dlnD = D*dh/dD = -h
dh_Wave_1[3] = -h_wave_1
dh_Wave_2[3] = -h_wave_2
dh_Wave_3[3] = -h_wave_3
parameters_plus = parameters.copy()
parameters_plus['psi'] = parameters['psi']+offset
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[2] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[2] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[2] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
if 4 in conf.param_indices:
# Deriving w.r.t. ln(Mc)
parameters_minus = buf.parameters_new_dlnMc_old(parameters,-offset)
parameters_plus = buf.parameters_new_dlnMc_old(parameters,offset)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
# Deriving w.r.t. ln(D) where D = parameters['luminosity_distance']
# Here since h ~ constant/D, dh/dD = -h/D. Hence dh/dlnD = D*dh/dD = -h
dh_Wave_1[3] = -h_wave_1
dh_Wave_2[3] = -h_wave_2
dh_Wave_3[3] = -h_wave_3
# If parameters['symmetric_mass_ratio'] = 0.25: a backward step will give parameters_minus['symmetric_mass_ratio'] > 0.25 which is unphysical, hence we only use forward difference
if parameters_minus['symmetric_mass_ratio'] > 0.25:
dh_Wave_1[4] = (h_wave_1_plus-h_wave_1)/offset
dh_Wave_2[4] = (h_wave_2_plus-h_wave_2)/offset
dh_Wave_3[4] = (h_wave_3_plus-h_wave_3)/offset
else:
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
dh_Wave_1[4] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[4] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[4] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. ln(Mc)
parameters_minus = buf.parameters_new_dlnMc_old(parameters,-offset)
parameters_plus = buf.parameters_new_dlnMc_old(parameters,offset)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
if 5 in conf.param_indices:
# Deriving w.r.t. ln(mu)
parameters_minus = buf.parameters_new_dlnmu_old(parameters,-offset)
parameters_plus = buf.parameters_new_dlnmu_old(parameters,offset)
# If parameters['symmetric_mass_ratio'] = 0.25: a backward step will give parameters_minus['symmetric_mass_ratio'] > 0.25 which is unphysical, hence we only use forward difference
if parameters_minus['symmetric_mass_ratio'] > 0.25:
dh_Wave_1[4] = (h_wave_1_plus-h_wave_1)/offset
dh_Wave_2[4] = (h_wave_2_plus-h_wave_2)/offset
dh_Wave_3[4] = (h_wave_3_plus-h_wave_3)/offset
else:
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
dh_Wave_1[4] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[4] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[4] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# If parameters['symmetric_mass_ratio'] = 0.25: a forward step will give parameters_plus['symmetric_mass_ratio'] > 0.25 which is unphysical, hence we only use backward difference
if parameters_plus['symmetric_mass_ratio'] > 0.25:
dh_Wave_1[5] = (h_wave_1-h_wave_1_minus)/offset
dh_Wave_2[5] = (h_wave_2-h_wave_2_minus)/offset
dh_Wave_3[5] = (h_wave_3-h_wave_3_minus)/offset
else:
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[5] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[5] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[5] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. ln(mu)
parameters_minus = buf.parameters_new_dlnmu_old(parameters,-offset)
parameters_plus = buf.parameters_new_dlnmu_old(parameters,offset)
if 6 in conf.param_indices:
# Deriving w.r.t. np.sin(dec) where dec = parameters['dec']
parameters_minus = parameters.copy()
parameters_minus['dec'] = np.arcsin(np.sin(parameters['dec'])-offset)
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
parameters_plus = parameters.copy()
parameters_plus['dec'] = np.arcsin(np.sin(parameters['dec'])+offset)
# If parameters['symmetric_mass_ratio'] = 0.25: a forward step will give parameters_plus['symmetric_mass_ratio'] > 0.25 which is unphysical, hence we only use backward difference
if parameters_plus['symmetric_mass_ratio'] > 0.25:
dh_Wave_1[5] = (h_wave_1-h_wave_1_minus)/offset
dh_Wave_2[5] = (h_wave_2-h_wave_2_minus)/offset
dh_Wave_3[5] = (h_wave_3-h_wave_3_minus)/offset
else:
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[5] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[5] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[5] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. np.sin(dec) where dec = parameters['dec']
parameters_minus = parameters.copy()
parameters_minus['dec'] = np.arcsin(np.sin(parameters['dec'])-offset)
dh_Wave_1[6] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[6] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[6] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
parameters_plus = parameters.copy()
parameters_plus['dec'] = np.arcsin(np.sin(parameters['dec'])+offset)
if 7 in conf.param_indices:
# Deriving w.r.t. ra = parameters['ra']
parameters_minus = parameters.copy()
parameters_minus['ra'] = parameters['ra']-offset
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
parameters_plus = parameters.copy()
parameters_plus['ra'] = parameters['ra']+offset
dh_Wave_1[6] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[6] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[6] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. ra = parameters['ra']
parameters_minus = parameters.copy()
parameters_minus['ra'] = parameters['ra']-offset
parameters_plus = parameters.copy()
parameters_plus['ra'] = parameters['ra']+offset
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
dh_Wave_1[7] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[7] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[7] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
# Deriving w.r.t. ln(tc) where tc = parameters['meta']['tc_3p5PN']
# dh/dtc = (h(tc+offset)-h(tc-offset))/(2*offset)
# = (h(f_ref)*np.exp(-i*2*pi*f*offset)-h(f_ref)*np.exp(-i*2*pi*f*-offset))/2*offset # if -1j convention used
# = -i*sin(2*pi*f*offset)*h(f_ref)/offset
# = -i*sin(2*pi*f*offset)*h(f_ref) * 2*pi*f/(2*pi*f*offset)
# = -i * 2*pi*f * h(f_ref) # since 2*pi*f*offset is very small
# dh/dlntc = tc*dh/dtc
# dh/dlntc = -i * 2*pi*f*tc * h(f_ref)
two_pi_f_tc = -1j * 2*np.pi*interferometers.frequency_array*(parameters['geocent_time']-start_time)
dh_Wave_1[8] = two_pi_f_tc * h_wave_1
dh_Wave_2[8] = two_pi_f_tc * h_wave_2
dh_Wave_3[8] = two_pi_f_tc * h_wave_3
h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
# parameters_minus = parameters.copy()
# parameters_minus['geocent_time'] = np.exp(np.log(parameters['geocent_time'] - start_time) - offset) + start_time
#
# parameters_plus = parameters.copy()
# parameters_plus['geocent_time'] = np.exp(np.log(parameters['geocent_time'] - start_time) + offset) + start_time
#
# h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
# h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
#
# dh_Wave_1[8] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
# dh_Wave_2[8] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
# dh_Wave_3[8] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
dh_Wave_1[7] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
dh_Wave_2[7] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
dh_Wave_3[7] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
if 8 in conf.param_indices:
# Deriving w.r.t. ln(tc) where tc = parameters['meta']['tc_3p5PN']
# dh/dtc = (h(tc+offset)-h(tc-offset))/(2*offset)
# = (h(f_ref)*np.exp(-i*2*pi*f*offset)-h(f_ref)*np.exp(-i*2*pi*f*-offset))/2*offset # if -1j convention used
# = -i*sin(2*pi*f*offset)*h(f_ref)/offset
# = -i*sin(2*pi*f*offset)*h(f_ref) * 2*pi*f/(2*pi*f*offset)
# = -i * 2*pi*f * h(f_ref) # since 2*pi*f*offset is very small
# dh/dlntc = tc*dh/dtc
# dh/dlntc = -i * 2*pi*f*tc * h(f_ref)
two_pi_f_tc = -1j * 2*np.pi*interferometers.frequency_array*(parameters['geocent_time']-start_time)
dh_Wave_1[8] = two_pi_f_tc * h_wave_1
dh_Wave_2[8] = two_pi_f_tc * h_wave_2
dh_Wave_3[8] = two_pi_f_tc * h_wave_3
# parameters_minus = parameters.copy()
# parameters_minus['geocent_time'] = np.exp(np.log(parameters['geocent_time'] - start_time) - offset) + start_time
#
# parameters_plus = parameters.copy()
# parameters_plus['geocent_time'] = np.exp(np.log(parameters['geocent_time'] - start_time) + offset) + start_time
#
# h_wave_1_minus, h_wave_2_minus, h_wave_3_minus = WaveForm_ThreeDetectors(parameters_minus, start_time, minimum_frequency, interferometers)
# h_wave_1_plus, h_wave_2_plus, h_wave_3_plus = WaveForm_ThreeDetectors(parameters_plus, start_time, minimum_frequency, interferometers)
#
# dh_Wave_1[8] = (h_wave_1_plus-h_wave_1_minus)/(2*offset)
# dh_Wave_2[8] = (h_wave_2_plus-h_wave_2_minus)/(2*offset)
# dh_Wave_3[8] = (h_wave_3_plus-h_wave_3_minus)/(2*offset)
return h_wave_1, h_wave_2, h_wave_3, dh_Wave_1, dh_Wave_2, dh_Wave_3
......@@ -1064,45 +1068,3 @@ def WaveForm_ThreeDetectors_with_noise(sigpar, n_freq):
############################################################################################################
##################### END of sigpar version of the functions ###############################################
############################################################################################################
def ComputeChirpTime3p5PN(f_low, m1, m2):
"""
Computes chirp times at 3.5 PN.
Inputs:
-------
f_low: [float] low frequency cutoff of the detector
sigpar: [array_like] parameters of the system
Outputs:
--------
"""
Mt = m1 + m2
Mt_sec = Mt * GEOM # Total mass in s
etaM = m1 * m2 / Mt**2
etaM2 = etaM*etaM
Tau_0 = tf1 * Mt_sec / etaM
Tau_2 = tf2 + tf3 * etaM
Tau_3 = tf4 * PI
Tau_4 = tf5 + tf6 * etaM + tf7 * etaM2
Tau_5 = PI*(tf9 + tf8*etaM)
Tau_6 = tf11 + tf12*PI2 + tf10*(EULER + np.log(4.0)) + ( tf13 + tf14*PI2 )*etaM + tf15*etaM2 + tf16*etaM2*etaM
Tau_6_log = tf10
Tau_7 = PI*(tf19 + tf18*etaM + tf17*etaM2 )
v_0 = pow(PI*Mt_sec*f_low,1.0/3.0)
v_0_2 = v_0*v_0
v_0_3 = v_0_2*v_0
v_0_4 = v_0_3*v_0
v_0_5 = v_0_4*v_0