Commit 09a919ce authored by Marc Arene's avatar Marc Arene
Browse files

Using ROQ in configuration where:

dlogL_roq.py gives likelihood.log_likelihood_ratio()=674.13
parent 993c10ed
......@@ -83,7 +83,6 @@ if __name__=='__main__':
inj_file_name = opts.inj_file.split('/')[-1]
inj_name = inj_file_name.split('.')[0]
binary_file_name = '../Codes/' + inj_name
SkipPhase1 = opts.SkipPhase1
verbose = opts.verbose
......
......@@ -9,7 +9,7 @@ from Headers.PN_Coefficients import * # Headers.Constants already imported in PN
import Library.bilby_waveform as bilby_wv
import Library.param_utils as paru
import Library.config as conf
# import Library.bilby_detector as bilby_detector
import Library.python_utils as pu
import Codes.set_injection_parameters as set_inj
import Codes.set_psds as set_psds
......@@ -17,7 +17,7 @@ import Codes.logL_snr as logL_snr
def main(injection_parameters, start_time, minimum_frequency, interferometers, do_prints=False):
import IPython; IPython.embed()
dlogL, logL = bilby_wv.dlogL_ThreeDetectors(injection_parameters, start_time, minimum_frequency, interferometers)
prefix = '\ndlogL with bilby '
......@@ -68,7 +68,10 @@ if __name__=='__main__':
# SET THE INJECTION PARAMETERS
injection_parameters = set_inj.ini_file_to_dict()
pu.print_dict(injection_parameters)
minimum_frequency = injection_parameters['meta']['minimum_frequency']
maximum_frequency = 2048
# maximum_frequency = injection_parameters['meta']['maximum_frequency']
duration = injection_parameters['meta']['duration']
sampling_frequency = injection_parameters['meta']['sampling_frequency']
f_high = injection_parameters['meta']['f_high']
......@@ -83,6 +86,7 @@ if __name__=='__main__':
for ifo in interferometers:
# ifo.__class__ = bilby_detector.Interferometer
ifo.minimum_frequency = minimum_frequency
ifo.maximum_frequency = maximum_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
......
......@@ -7,6 +7,7 @@ import bilby
import Library.bilby_utils as bilby_utils
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.bilby_waveform as bilby_wv
import Library.python_utils as pu
import Codes.set_injection_parameters as set_inj
import Codes.set_psds as set_psds
......@@ -46,7 +47,9 @@ if __name__=='__main__':
# SET THE INJECTION PARAMETERS
injection_parameters = set_inj.ini_file_to_dict()
pu.print_dict(injection_parameters)
minimum_frequency = injection_parameters['meta']['minimum_frequency']
maximum_frequency = injection_parameters['meta']['maximum_frequency']
duration = injection_parameters['meta']['duration']
sampling_frequency = injection_parameters['meta']['sampling_frequency']
f_high = injection_parameters['meta']['f_high']
......@@ -60,6 +63,7 @@ if __name__=='__main__':
interferometers = bilby.gw.detector.InterferometerList(ifos_chosen)
for ifo in interferometers:
ifo.minimum_frequency = minimum_frequency
# ifo.maximum_frequency = maximum_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
......
......@@ -10,6 +10,33 @@ from Headers.PN_Coefficients import * # Headers.Constants already imported in PN
import Library.python_utils as pu
import Library.param_utils as paru
import Library.config as conf
import Library.roq_utils as roqu
LAL_PI = 3.141592653589793238462643383279502884
LAL_MSUN_SI = 1.988546954961461467461011951140572744e30
LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
def Compute_fISCO(Rmin, m1, m2):
"""
Compute the frequency at the Innermost Stable Circular Orbit
"""
vISCO = np.sqrt(1.0 / Rmin)
Mt = m1 + m2
fmax = vISCO**3 / (PI * GEOM * Mt) # fmax signal
return fmax
def Compute_LAL_fISCO(m1_SI, m2_SI):
m1 = m1_SI / LAL_MSUN_SI
m2 = m2_SI / LAL_MSUN_SI
m_sec = (m1 + m2) * LAL_MTSUN_SI
piM = LAL_PI * m_sec
vISCO = 1 / np.sqrt(6)
fISCO = vISCO**3 / piM
return fISCO
def ComputeChirpTime3p5PN(f_low, m1, m2):
......@@ -48,11 +75,20 @@ def ComputeChirpTime3p5PN(f_low, m1, m2):
Tau_chirp = Tau_0*pow(v_0,-8.0)*(1.0 + Tau_2*v_0_2 + Tau_3*v_0_3 + Tau_4 *v_0_4 + Tau_5 *v_0_5 + (Tau_6 + Tau_6_log*np.log(v_0))*v_0_6 + Tau_7 *v_0_7 )
return Tau_chirp
def Compute_fmax(Rmin, m1, m2):
vmax = np.sqrt(1.0 / Rmin)
Mt = m1 + m2
fmax = vmax**3 / (PI * GEOM * Mt) # fmax signal
return fmax
def compute_time_from_flow_to_fhigh_3p5PN(f_low, f_high, mass_1, mass_2):
"""
Computes the time it takes for the binary to go from f_low to f_high, using 3.5 PN approximation.
"""
# Compute the time it takes from flow to get to merger
tc_flow = ComputeChirpTime3p5PN(f_low, mass_1, mass_2)
# Compute the time it takes from fhigh to get to merger
tc_fhigh = ComputeChirpTime3p5PN(f_high, mass_1, mass_2)
# Substracting the two gives the time it take for the binary to go from f_low to f_high
return tc_flow - tc_fhigh
def ini_file_to_dict(file_path='../examples/GW170817.ini'):
......@@ -71,6 +107,7 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
# to inject the signal into. For the
# TaylorF2 waveform, we cut the signal close to the isco frequency
minimum_frequency=30.0
# minimum_frequency=30.0
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
......@@ -78,11 +115,58 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
R_isco = 6. # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
m1_min = 1
m2_min = 1
approximant = 'TaylorF2'
# approximant = 'TaylorF2'
# approximant = 'SpinTaylorF2'
# approximant = 'IMRPhenomD'
# approximant = 'IMRPhenomPv2'
approximant = 'IMRPhenomPv2_ROQ'
# approximant = 'IMRPhenomPv2_NRTidal'
f_ISCO = Compute_fISCO(R_isco, m1, m2)
f_ISCO_run = Compute_fISCO(R_isco, m1_min, m2_min)
if approximant == 'TaylorF2':
f_high = Compute_fmax(R_isco, m1, m2)
f_high_run = Compute_fmax(R_isco, m1_min, m2_min)
elif approximant = 'IMRPhenomD':
f_high = f_ISCO
f_high_run = f_ISCO_run
# maximum_frequency = min(4096 * 4, f_high_run)
maximum_frequency = 0
sampling_frequency = 2 * f_high_run
duration = tc_3p5PN + 2
reference_frequency = 0
elif approximant == 'IMRPhenomPv2_ROQ':
approximant = 'IMRPhenomPv2'
# cf:
# - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
# - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
roq_directory = '/Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/'
params_best_basis, scale_factor = roqu.get_best_segment_params_from_chirp_mass(injection_parameters['chirp_mass'], roq_directory=roq_directory)
injection_parameters['roq'] = {}
injection_parameters['roq']['directory'] = roq_directory + str(int(params_best_basis['seglen'])) + 's/'
injection_parameters['roq']['params_path'] = injection_parameters['roq']['directory'] + 'params.dat'
rescaled_params = roqu.rescale_params(params_best_basis, scale_factor)
minimum_frequency = rescaled_params['flow']
maximum_frequency = rescaled_params['fhigh']
# maximum_frequency = min(2048, rescaled_params['fhigh'])
# maximum_frequency = min(4 * 4096, f_high_run)
sampling_frequency = 2 * max(2048, rescaled_params['fhigh'])
duration = rescaled_params['seglen']
reference_frequency = 20 * scale_factor
injection_parameters['roq']['params'] = params_best_basis
injection_parameters['roq']['rescaled_params'] = rescaled_params
injection_parameters['roq']['scale_factor'] = scale_factor
tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
f_CUT = 0.2
M_sec = (m1 + m2) * LAL_MTSUN_SI
M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
f_high = f_CUT / M_sec
f_high_run = f_CUT / M_min_sec
else:
# cf:
# - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
# - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
......@@ -92,25 +176,37 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
f_high = f_CUT / M_sec
f_high_run = f_CUT / M_min_sec
maximum_frequency = min(2048, f_high_run)
# maximum_frequency = min(4 * 4096, f_high_run)
sampling_frequency = 2 * maximum_frequency
duration = 64
reference_frequency = 20
# import IPython; IPython.embed();
duration = tc_3p5PN + 2
# maximum_frequency = min(4096 * 4, f_high_run)
# maximum_frequency = max(4096, f_high)
# sampling_frequency = 4096
# sampling_frequency = 4397
sampling_frequency = 2 * f_high_run
start_time = injection_parameters['geocent_time'] - tc_3p5PN
# breakpoint()
meta_params = dict(
minimum_frequency = minimum_frequency,
maximum_frequency = maximum_frequency,
sampling_frequency = sampling_frequency,
duration = duration,
tc_3p5PN = tc_3p5PN,
start_time = start_time,
f_high = f_high,
f_high_run = f_high_run,
duration = duration,
sampling_frequency = sampling_frequency,
start_time = start_time
approximant = approximant,
reference_frequency = reference_frequency
)
injection_parameters['meta'] = meta_params
conf.minimum_frequency = minimum_frequency
conf.approximant = approximant
conf.maximum_frequency = maximum_frequency
conf.reference_frequency = reference_frequency
# These two lines have been added to be able to plot the injected value on phi_c on the corner plots. Otherwise the `phi_c` and `tc_3p5PN` keys are not used in the code.
injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
......@@ -123,7 +219,17 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
if __name__=='__main__':
from optparse import OptionParser
usage = """%prog [options]
Plots the three waveforms in the time domain"""
parser=OptionParser(usage)
parser.add_option("--inj_file", default='../examples/GW170817.ini', action="store", type="string", help="""Injection file path leading to the `.ini` file specifying the parameters of the injection.""")
injection_parameters = ini_file_to_dict()
(opts,args)=parser.parse_args()
injection_parameters = ini_file_to_dict(opts.inj_file)
pu.print_dict(injection_parameters, indent=3, align_keys=True)
breakpoint()
......@@ -13,16 +13,10 @@ import Library.config as conf
def get_source_frame_polarizations(parameters, minimum_frequency, interferometers):
# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
# approximant = 'TaylorF2'
# approximant = 'SpinTaylorF2'
# approximant = 'IMRPhenomPv2'
approximant = 'IMRPhenomD'
if approximant == 'TaylorF2':
maximum_frequency = 0
if 'meta' in parameters.keys():
waveform_arguments = dict(waveform_approximant=parameters['meta']['approximant'], reference_frequency=parameters['meta']['reference_frequency'], minimum_frequency=minimum_frequency, maximum_frequency=parameters['meta']['maximum_frequency'])
else:
maximum_frequency = 4096
# waveform_arguments = dict(waveform_approximant=approximant, reference_frequency=50, minimum_frequency=minimum_frequency)
waveform_arguments = dict(waveform_approximant=approximant, reference_frequency=0, minimum_frequency=minimum_frequency, maximum_frequency=maximum_frequency)
waveform_arguments = dict(waveform_approximant=conf.approximant, reference_frequency=conf.reference_frequency, minimum_frequency=conf.minimum_frequency, maximum_frequency=conf.maximum_frequency)
# Create the waveform_generator using a LAL Binary Neutron Star source function
waveform_generator = bilby.gw.WaveformGenerator(
......@@ -32,7 +26,7 @@ def get_source_frame_polarizations(parameters, minimum_frequency, interferometer
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
waveform_arguments=waveform_arguments
)
pu.print_dict({"waveform_arguments": waveform_arguments})
# [Marc]: h+ and hx in source frame
source_frame_polarizations = waveform_generator.frequency_domain_strain(parameters)
......@@ -768,6 +762,7 @@ if __name__=='__main__':
# SET THE INJECTION PARAMETERS
injection_parameters = set_inj.ini_file_to_dict()
pu.print_dict(injection_parameters)
minimum_frequency = injection_parameters['meta']['minimum_frequency']
duration = injection_parameters['meta']['duration']
sampling_frequency = injection_parameters['meta']['sampling_frequency']
......@@ -788,7 +783,7 @@ if __name__=='__main__':
ifo.strain_data.duration = duration
ifo.strain_data.sampling_frequency = sampling_frequency
ifo.strain_data.start_time = start_time
# import IPython; IPython.embed()
template_ifos = WaveForm_ThreeDetectors(injection_parameters, minimum_frequency, interferometers)
......@@ -809,7 +804,7 @@ if __name__=='__main__':
f_low = minimum_frequency
indices = np.where((f_low<interferometers.frequency_array) & (interferometers.frequency_array<f_high))[0]
# indices = np.where((f_low<interferometers.frequency_array) & (interferometers.frequency_array<f_high))[0]
h_wave_td = []
......@@ -824,23 +819,34 @@ if __name__=='__main__':
number_of_samples = len(h_wave_td[ifo_index])
duration_td = number_of_samples / sampling_frequency
time = np.linspace(start=0, stop=duration_td, num=number_of_samples)
time = np.linspace(start=0, stop=duration_td, num=number_of_samples)
idx_start = np.where(time>56.915)[0][0]
idx_stop = np.where(time<56.935)[0][-1]
fISCO = set_inj.Compute_fISCO(6, injection_parameters['mass_1'], injection_parameters['mass_2'])
time_until_fISCO = set_inj.compute_time_from_flow_to_fhigh_3p5PN(minimum_frequency, fISCO, injection_parameters['mass_1'], injection_parameters['mass_2'])
if plot_on_same_figure:
plt.figure(figsize=(15,5))
plt.figure(figsize=(15,10))
for ifo_index in range(len(interferometers)):
plt.plot(time, h_wave_td[ifo_index], linewidth=1, color=ifo_colors[ifo_index], label=interferometers[ifo_index].name)
# plt.plot(time[idx_start:idx_stop], h_wave_td[ifo_index][idx_start:idx_stop], linewidth=1, color=ifo_colors[ifo_index], label=interferometers[ifo_index].name)
plt.axvline(x=time_until_fISCO, linewidth=1, color='blue', label='time to reach fISCO')
plt.axvline(x=injection_parameters['meta']['tc_3p5PN'], color='black', linewidth=1, label='tc')
text = 'f_low = {:.2f}Hz\nf_ISCO = {:.2f}Hz\nmax_freq = {:.2f}Hz\nsamp_freq= {:.2f}Hz\ntc_3p5PN = {:.4f}s'.format(minimum_frequency, fISCO, injection_parameters['meta']['maximum_frequency'], sampling_frequency, injection_parameters['meta']['tc_3p5PN'])
plt.text(duration-7, h_wave_td[0].max()*2/3, text, fontsize=10)
# import IPython; IPython.embed()
plt.xlabel("Time (sec)")
plt.ylabel("h")
plt.legend()
plt.title(" TaylorF2 GW170817, fmin={}Hz".format(f_low))
plt.title(" GW170817 - {}".format(injection_parameters['meta']['approximant']))
else:
for ifo_index in range(len(interferometers)):
plt.figure(figsize=(15,5))
plt.plot(time, h_wave_td[ifo_index], linewidth=1, color=ifo_colors[ifo_index])
plt.xlabel("Time (sec)")
plt.ylabel("h")
plt.title(title_prefix + " TaylorF2 GW170817, fmin={}Hz, in {}".format(f_low, interferometers[ifo_index].name))
plt.title(title_prefix + "GW170817 - {}, fmin={}Hz".format(injection_parameters['meta']['approximant'], f_low))
plt.show()
......
......@@ -14,7 +14,7 @@ random_numbers = []
# 7 = ra
# 8 = ln(tc)
param_names = ["iota", "phi_c", "psi", "D", "Mc", "mu", "dec", "ra", "tc"]
fparam_names = ["cos(iota)", "phi_c", "psi", "ln(D)", "ln(Mc)", "ln(mu)", "sin(dec)", "ra", "ln(tc)"]
fparam_names = ["cos(theta_jn)", "phi_c", "psi", "ln(D)", "ln(Mc)", "ln(mu)", "sin(dec)", "ra", "ln(tc)"]
param_dict_keys = ["iota", "phase", "psi", "luminosity_distance", "chirp_mass", "reduced_mass", "dec", "ra", "geocent_time"]
parameter_keys=['theta_jn', 'phi_c', 'psi', 'luminosity_distance', 'chirp_mass', 'reduced_mass', 'dec', 'ra', 'tc_3p5PN']
......@@ -35,6 +35,10 @@ param_indices = [0,1,2,3,4,5,6,7,8]
# param_indices = [7]
dlogL = 'dlogL1'
minimum_frequency = 30
approximant = 'TaylorF2'
maximum_frequency = 0
reference_frequency = 0
debug = False
......
import sys
# # cf https:#docs.python.org/3/tutorial/modules.html
sys.path.append('../')
import numpy as np
from bilby.gw.likelihood import GravitationalWaveTransient
from bilby.core.likelihood import Likelihood
import Library.param_utils as paru
class LikelihoodGradient(object):
"""
"""
def __init__(self, likelihood, search_parameters_keys, offset=1e-7):
self.likelihood = likelihood
self.search_parameters_keys = search_parameters_keys
self.offset = offset
self.log_likelihood_gradient = dict()
for key in search_parameters_keys:
self.log_likelihood_gradient[key] = None
def calculate_log_likelihood_gradient(self):
self.likelihood.parameters_copy = self.likelihood.parameters.copy()
for key in self.search_parameters_keys:
self.likelihood.parameters[key] -= self.offset
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[key] += self.offset
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
return self.log_likelihood_gradient
class GWTransientLikelihoodGradient(LikelihoodGradient):
"""
"""
def __init__(self, likelihood, search_parameters_keys, offset=1e-7):
LikelihoodGradient.__init__(self, likelihood, search_parameters_keys, offset)
self.log_likelihood_at_point = self.likelihood.log_likelihood_ratio()
self.start_time = likelihood.interferometers[0].strain_data.start_time
def calculate_log_likelihood_gradient(self):
self.likelihood.parameters_copy = self.likelihood.parameters.copy()
for key in self.search_parameters_keys:
if key=='cos(theta_jn)':
lalsim_key = 'theta_jn'
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
# If we are at the boundary where iota = +/- pi, we only do a forward difference in cos(iota)
if np.cos(self.likelihood.parameters_copy[lalsim_key]) - self.offset < -1:
self.likelihood.parameters[lalsim_key] = np.arccos(np.cos(self.likelihood.parameters_copy[lalsim_key]) + self.offset)
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - self.log_likelihood_at_point) / self.offset
# If we are at the boundary where iota = 0, we only do a backward difference in cos(iota)
elif np.cos(self.likelihood.parameters_copy[lalsim_key]) + self.offset > 1:
self.likelihood.parameters[lalsim_key] = np.arccos(np.cos(self.likelihood.parameters_copy[lalsim_key]) - self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (self.log_likelihood_at_point - log_likelihood_minus) / self.offset
# Otherwise we keep the central difference
else:
self.likelihood.parameters[lalsim_key] = np.arccos(np.cos(self.likelihood.parameters_copy[lalsim_key]) - self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.likelihood.parameters[lalsim_key] = np.arccos(np.cos(self.likelihood.parameters_copy[lalsim_key]) + self.offset)
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='phi_c':
lalsim_key = 'phase'
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[lalsim_key] = self.likelihood.parameters_copy[lalsim_key] - self.offset / 2
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.likelihood.parameters[lalsim_key] = self.likelihood.parameters_copy[lalsim_key] + self.offset / 2
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='ln(D)':
lalsim_key = 'luminosity_distance'
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[lalsim_key] = np.exp(np.log(self.likelihood.parameters_copy[lalsim_key]) - self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.likelihood.parameters[lalsim_key] = np.exp(np.log(self.likelihood.parameters_copy[lalsim_key]) + self.offset)
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='ln(Mc)':
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters = paru.parameters_new_dlnMc(self.likelihood.parameters_copy, self.offset)
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
parameters_minus = paru.parameters_new_dlnMc(self.likelihood.parameters_copy, -self.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:
self.log_likelihood_gradient[key] = (log_likelihood_plus - self.log_likelihood_at_point) / self.offset
else:
self.likelihood.parameters = parameters_minus
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='ln(mu)':
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters = paru.parameters_new_dlnmu(self.likelihood.parameters_copy, -self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
parameters_plus = paru.parameters_new_dlnMc(self.likelihood.parameters_copy, self.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:
self.log_likelihood_gradient[key] = (self.log_likelihood_at_point - log_likelihood_minus) / self.offset
else:
self.likelihood.parameters = parameters_plus
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='sin(dec)':
lalsim_key = 'dec'
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[lalsim_key] = np.arcsin(np.sin(self.likelihood.parameters_copy[lalsim_key]) - self.offset)
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.likelihood.parameters[lalsim_key] = np.arcsin(np.sin(self.likelihood.parameters_copy[lalsim_key]) + self.offset)
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
elif key=='ln(tc)':
lalsim_key = 'geocent_time'
# Reset the parameters to their original value
self.likelihood.parameters = self.likelihood.parameters_copy.copy()
self.likelihood.parameters[lalsim_key] = np.exp(np.log(self.likelihood.parameters_copy[lalsim_key] - self.start_time) - self.offset) + self.start_time
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
# Reset the parameters to their original value
self.likelihood.parameters[lalsim_key] = self.likelihood.parameters_copy[lalsim_key]
self.likelihood.parameters[lalsim_key] = np.exp(np.log(self.likelihood.parameters_copy[lalsim_key] - self.start_time) + self.offset) + self.start_time
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
# The other parameters handled by this else case are: psi, ra
else:
self.likelihood.parameters[key] = self.likelihood.parameters_copy[key] - self.offset
log_likelihood_minus = self.likelihood.log_likelihood_ratio()
self.likelihood.parameters[key] = self.likelihood.parameters_copy[key] + self.offset
log_likelihood_plus = self.likelihood.log_likelihood_ratio()
self.log_likelihood_gradient[key] = (log_likelihood_plus - log_likelihood_minus) / (2 * self.offset)
return self.log_likelihood_gradient
import sys
sys.path.append('../')
import Library.python_utils as pu
import numpy as np
def get_bases_params(roq_directory = '/Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/'):
segment_lenghts = ['4s', '8s', '16s', '32s', '64s', '128s']
params_bases = {}
for segment in segment_lenghts:
# params_bases[segment] = np.genfromtxt(roq_directory + segment + '/' + "params.dat", names=True)
params_bases[segment] = {}
params_of_segment = np.genfromtxt(roq_directory + segment + '/' + "params.dat", names=True)
for key in params_of_segment.dtype.names:
params_bases[segment][key] = params_of_segment[key] * 1
return params_bases
def find_best_segment_from_chirp_mass(chirp_mass, params_bases):
"""
Return the best ROQ segment adapted to the chirp_mass which has triggered together with the scale factor needed to apply to that segment if the chirp mass is out of the segment range.
Since segments overlap in chirp mass ranges, the best segment is taken to be the one where the chirp mass triggered is closest to its average chirp mass range value.
"""
min_distance = 999 #randomly big number
best_segment = ''
segments_scale_factors = {}
for segment in params_bases.keys():
chirp_mass_range_average = (params_bases[segment]['chirpmassmin'] + params_bases[segment]['chirpmassmax']) / 2
distance_to_average_chirp_mass = abs(chirp_mass - chirp_mass_range_average)
if distance_to_average_chirp_mass < min_distance:
min_distance = distance_to_average_chirp_mass
best_segment = segment
segments_scale_factors[segment] = chirp_mass_range_average / chirp_mass
# If the chirp mass triggered is in the range of chirp masses of the best segment then we leave the scale factor to one
if params_bases[best_segment]['chirpmassmin'] < chirp_mass and chirp_mass < params_bases[best_segment]['chirpmassmax']:
scale_factor = 1
else:
scale_factor = round(segments_scale_factors[best_segment], 1)
return params_bases[best_segment], scale_factor
def get_best_segment_params_from_chirp_mass(chirp_mass, roq_directory = '/Users/marcarene/Documents/APC/Code/ROQ_data/IMRPhenomPv2/'):
params_bases = get_bases_params(roq_directory)
params_best_basis, scale_factor = find_best_segment_from_chirp_mass(chirp_mass, params_bases)