Commit 7c5dbb56 authored by Marc Arene's avatar Marc Arene
Browse files

not much

parent cb3069e4
......@@ -67,6 +67,7 @@ if __name__=='__main__':
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# Parse the inputs relative to the HMC algorithm
n_traj_hmc_tot = config_dict['hmc']['n_traj_hmc_tot']
n_traj_fit = config_dict['hmc']['n_traj_fit']
n_phase1_for_this_run = config_dict['hmc']['n_phase1_for_this_run']
......@@ -75,6 +76,7 @@ if __name__=='__main__':
length_num_traj = config_dict['hmc']['length_num_traj']
epsilon0 = config_dict['hmc']['epsilon0']
# Parse the inputs relative to the analysis
ifos_possible = ['H1', 'L1', 'V1']
conf.ifos = config_dict['analysis']['ifos'].split(',')
if set.intersection(set(ifos_possible), set(conf.ifos)) != set(conf.ifos):
......@@ -348,30 +350,31 @@ if __name__=='__main__':
logL = likelihood.log_likelihood_ratio()
print("")
print("likelihood.log_likelihood_ratio() = {}".format(logL))
breakpoint()
likelihood_gradient = lg.GWTransientLikelihoodGradient(likelihood, search_parameters_keys)
dlogL_dict = likelihood_gradient.calculate_log_likelihood_gradient()
pu.print_dict({"dlogL_ROQ": dlogL_dict})
dlogL = likelihood_gradient.convert_dlogL_dict_to_np_array(dlogL_dict)
else:
# Here we add constraints on chirp mass and mass ratio to the prior, these are
# determined by the domain of validity of the ROQ basis.
priors = bilby.gw.prior.BNSPriorDict()
# for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'luminosity_distance']:
priors[key] = injection_parameters[key]
for key in ['mass_1', 'mass_2']:
priors[key].minimum = max(priors[key].minimum, priors['comp-min'])
priors['chirp_mass'] = bilby.core.prior.Constraint(name='chirp_mass', minimum=float(params['chirpmassmin']), maximum=float(params['chirpmassmax']))
priors['mass_ratio'] = bilby.core.prior.Constraint(0.125, 1, name='mass_ratio')
priors['geocent_time'] = bilby.core.prior.Uniform(injection_parameters['geocent_time'] - 0.1, injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
interferometers=interferometers, waveform_generator=search_waveform_generator,
priors=priors, phase_marginalization=config_dict['analysis']['phase_marginalization'])
search_parameters_keys = [conf.fparam_names[i] for i in conf.param_indices]
likelihood_gradient = lg.GWTransientLikelihoodGradient(likelihood, search_parameters_keys)
# # Here we add constraints on chirp mass and mass ratio to the prior, these are
# # determined by the domain of validity of the ROQ basis.
# priors = bilby.gw.prior.BNSPriorDict()
# # for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
# for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra', 'dec', 'luminosity_distance']:
# priors[key] = injection_parameters[key]
# for key in ['mass_1', 'mass_2']:
# priors[key].minimum = max(priors[key].minimum, priors['comp-min'])
# priors['chirp_mass'] = bilby.core.prior.Constraint(name='chirp_mass', minimum=float(params['chirpmassmin']), maximum=float(params['chirpmassmax']))
# priors['mass_ratio'] = bilby.core.prior.Constraint(0.125, 1, name='mass_ratio')
# priors['geocent_time'] = bilby.core.prior.Uniform(injection_parameters['geocent_time'] - 0.1, injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
#
# likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
# interferometers=interferometers, waveform_generator=search_waveform_generator,
# priors=priors, phase_marginalization=config_dict['analysis']['phase_marginalization'])
#
# search_parameters_keys = [conf.fparam_names[i] for i in conf.param_indices]
# likelihood_gradient = lg.GWTransientLikelihoodGradient(likelihood, search_parameters_keys)
likelihood_gradient = None
......@@ -896,7 +899,6 @@ if __name__=='__main__':
# has_NaN = np.isnan(oluts_dict[key][0]).any() or np.isnan(oluts_dict[key][1]).any() or has_NaN
# if has_NaN:
# print(f'At traj {traj_index}, oluts_dict has NaNs !!')
# breakpoint()
fileOutChain.closed
......
......@@ -146,19 +146,33 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
rescaled_params = roqu.rescale_params(params_best_basis, scale_factor)
# minimum_frequency = max(conf.minimum_frequency, rescaled_params['flow'])
# maximum_frequency_ifo = min(int(conf.maximum_frequency), int(rescaled_params['fhigh']))
# sampling_frequency = 2 * maximum_frequency_ifo
# tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
# # int() function because we need to make sure that sampling_freq * duration = integer to the 14 decimal precision
# duration = min(int(tc_3p5PN + 2 + 1), int(rescaled_params['seglen']))
# maximum_frequency_injected_waveform = rescaled_params['fhigh']
# maximum_frequency_search_waveform = maximum_frequency_ifo
# reference_frequency = conf.reference_frequency * scale_factor
minimum_frequency = rescaled_params['flow']
maximum_frequency_injected_waveform = rescaled_params['fhigh']
# maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh'])
maximum_frequency_ifo = min(conf.maximum_frequency, rescaled_params['fhigh'])
maximum_frequency_search_waveform = maximum_frequency_injected_waveform
# maximum_frequency_search_waveform = maximum_frequency_ifo
# maximum_frequency_ifo = rescaled_params['fhigh']
sampling_frequency = 2 * max(conf.maximum_frequency, rescaled_params['fhigh'])
duration = rescaled_params['seglen']
reference_frequency = conf.reference_frequency * scale_factor
injection_parameters['roq']['params'] = params_best_basis
injection_parameters['roq']['rescaled_params'] = rescaled_params
injection_parameters['roq']['scale_factor'] = scale_factor
......
# How to set input parameters when using an ROQ basis
- An ROQ basis has certain ranges of validity for different parameters which must be taken into account when setting the values of the minimum, maximum frequency, chirp mass etc.
- Once `roq_params` is set we can derive the other parameters as follow:
-
[analysis]
# Interferometers to run the analysis on.
ifos = H1,L1,V1
approx = IMRPhenomD
minimum_frequency = 30.0
maximum_frequency = 2048.0
reference_frequency = 20.0
roq = False
roq_b_matrix_directory = /Users/marcarene/roq/ROQ_data/IMRPhenomPv2/
# 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.
psd = 3
# Indices of the parameters to run the hmc on. Other are kept fixed.
param_indices = 0,1,2,3,4,5,6,7,8
# Indices of the parameters where numerical gradients are going to be used instead of analytical formula if using the `FullHybridGradient_LeapFrog()` function
param_indices_fullhybrid = 0,2,3
# -log10() of the offset used for central differencing on each parameter. Default is '707077777' meaning conf.param_offsets = [1e-07, 0, 1e-07, 0, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07]. 0 Means analytical formula is used if available.
param_offsets = 7,0,7,0,7,7,7,7,7
# Method to use to compute the gradient of the log-likelihood. Use `dlogL1` to use dlogL = <s|dh> - <h|dh>, `dlogL2` for dlogL = logL(p+offset) - logL(p-offset).
dlogL = dlogL1
phase_marginalization = False
[hmc]
# Total number of trajectories of HMC
n_traj_hmc_tot = 1500
# Number of trajectories using numericalgradients in phase 1.
n_traj_fit = 1500
# Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
n_phase1_for_this_run = 800
# Number of points used to sub-select points in Ordered Look-Up Tables for local fit bimodal parameters like (cos(inc), psi, logD)
n_fit_1 = 2000
# Number of points sub-selected from the above n1 ones to perform a linera fit
n_fit_2 = 200
# Length of numerical trajectories during phase1
length_num_traj = 200
# = epsilon * 1000 where espilon is the stepsize of trajectories used for phase1.
epsilon0 = 5
[priors]
# Priors
# For all parameters known to lalinference, the min and max default prior can be overwritten with
# parname-min=MIN
# parname-max=MAX
# The starting point for the MCMC chain(s) can be specified with
# parname=START
# Parameters can be fixed to some value with
# fix-parname=FIXVALUE
# currently known parameters, together with default [min-max] are (radians for angle, Mpc for distance, Msun for masses)
# time Waveform time [trigtime-0.1-trigtime+0.1]
# chirpmass Chirpmass [1,15.3]
# eta Symmetric massratio (needs --use-eta) [0,0.0312]
# q Asymmetric massratio (a.k.a. q=m2/m1 with m1>m2) [0.033,1]
# phase Coalescence phase [0,2Pi]
# costheta_jn Cosine of angle between J and line of sight [-1,1]
# logdistance Log Distance [log(1),log(2000)]
distance-max=70
distance-min=1e-6
# rightascension Rightascension [0,2Pi]
# declination Declination [-Pi/2,Pi/2]
# polarisation Polarisation angle [0,Pi]
# Spin Parameters:
# a_spin1 Spin1 magnitude [-1,1] for aligned, [0,1] for precessing
# a_spin2 Spin2 magnitude [-1,1] for aligned, [0,1] for precessing
# tilt_spin1 Angle between spin1 and orbital angular momentum [0,Pi]
# tilt_spin2 Angle between spin2 and orbital angular momentum [0, Pi]
# phi_12 Difference between spins' azimuthal angles [0,2Pi]
# phi_jl Difference between total and orbital angular momentum azimuthal angles [0,2Pi]
# Tidal parameters (requires tidal= or tidalT=):
# lambda1 lambda1 [0,3000]
# lambda2 lambda2 [0,3000]
# lambdaT lambdaT [0,3000]
# dLambdaT dLambdaT [-500,500]
# Equation of State parameters (requires 4PolyEOS=):
# logp1 logp1 [33.6,35.4]
# gamma1 gamma1 [2.0,4.5]
# gamma2 gamma2 [1.1,4.5]
# gamma3 gamma3 [1.1,4.5]
# Settings for allowed component masses in Solar Masses, with default values
comp-max=2.6
comp-min=1.0
# Allowed total masses in Msun (note, used together with component masses, mc,q,eta priors may lead to inconsistencies. Be careful!)
# mtotal-max=35
# mtotal-min=2
# Setting time prior [seconds]
# dt=0.1
......@@ -10,7 +10,7 @@ roq_b_matrix_directory = /Users/marcarene/roq/ROQ_data/IMRPhenomPv2/
# 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.
psd = 3
# Indices of the parameters to run the hmc on. Other are kept fixed.
param_indices = 1,2
param_indices = 0,1,2,3,4,5,6,7,8
# Indices of the parameters where numerical gradients are going to be used instead of analytical formula if using the `FullHybridGradient_LeapFrog()` function
param_indices_fullhybrid = 0,2,3
# -log10() of the offset used for central differencing on each parameter. Default is '707077777' meaning conf.param_offsets = [1e-07, 0, 1e-07, 0, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07]. 0 Means analytical formula is used if available.
......@@ -22,11 +22,11 @@ phase_marginalization = False
[hmc]
# Total number of trajectories of HMC
n_traj_hmc_tot = 100
n_traj_hmc_tot = 10
# Number of trajectories using numericalgradients in phase 1.
n_traj_fit = 100
n_traj_fit = 10
# Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
n_phase1_for_this_run = 1500
n_phase1_for_this_run = 800
# Number of points used to sub-select points in Ordered Look-Up Tables for local fit bimodal parameters like (cos(inc), psi, logD)
n_fit_1 = 2000
# Number of points sub-selected from the above n1 ones to perform a linera fit
......
......@@ -26,3 +26,5 @@ lambda_1=0
lambda_2=0
tilt_1=0
tilt_2=0
a_1=0
a_2=0
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment