Commit 6cfb77d8 authored by Marc Arene's avatar Marc Arene
Browse files

Comments main.py until line 518.

parent 8430493b
......@@ -185,11 +185,15 @@ if __name__ == '__main__':
files_in_event_dir = os.listdir(hdf5_file_path_prefix)
# INITIALIZE THE THREE INTERFEROMETERS
# We use a home-made encapsulation of the standard Interferometer Bilby object
# because we have added some optimized ifo methods for the HMC.
# interferometers = bilby.gw.detector.InterferometerList(ext_analysis_dict['ifos'])
interferometers = gwdn.InterferometerList(ext_analysis_dict['ifos'])
# The minimun frequency for each interferometer needs to be set if we want the strain derived from the psd to start being non zero at that frequency and not at the default one which is 20Hz
# For each ifo we will extract the time domain strain segment, attach it
# to the ifo and then set its PSD.
for i, ifo in enumerate(interferometers):
# The minimun frequency for each interferometer needs to be set if we want the strain derived from the psd to start being non zero at that frequency and not at the default one which is 20Hz
ifo.minimum_frequency = ext_analysis_dict['minimum_frequency_ifos'][ifo.name]
ifo.maximum_frequency = ext_analysis_dict['maximum_frequency_ifo']
......@@ -203,6 +207,7 @@ if __name__ == '__main__':
hdf5_file_name = list(filter(lambda file: ifo_prefix in file, files_in_event_dir))[0]
hdf5_file_path = f'{hdf5_file_path_prefix}/{hdf5_file_name}'
# Extract the strain and other information from the hdf5 file
hdf5_strain, hdf5_start_time, hdf5_duration, sampling_frequency, hdf5_number_of_samples = psd_utils.get_strain_from_hdf5(hdf5_file_path)
if ext_analysis_dict['sampling_frequency'] > sampling_frequency:
cut.logger.error('The sampling frequency required for this analysis: {}Hz, is higher than the sampling frequency of the input strain file: {}Hz'.format(ext_analysis_dict['sampling_frequency'], sampling_frequency))
......@@ -214,45 +219,55 @@ if __name__ == '__main__':
hdf5_strain = hdf5_strain[::downsampling_factor]
sampling_frequency = ext_analysis_dict['sampling_frequency']
# Find the index of the first sample of the segment which will be used for the analysis.
idx_segment_start_float = (start_time - hdf5_start_time) * sampling_frequency
if idx_segment_start_float == int(idx_segment_start_float):
idx_segment_start = int(idx_segment_start_float)
# If the `start_time` does not fall exactly on a sample (which is generally the case),
# we shift its value to be that of the sample.
else:
idx_segment_start = int(idx_segment_start_float) + 1
start_time_offset = (idx_segment_start - idx_segment_start_float) / sampling_frequency
start_time +=start_time_offset
cut.logger.info(f'Due to numerical precison when selecting the segment, the start_time has been offset by +{start_time_offset:.4e} to equal {start_time}')
# Normally everything should be fine for the end of the segment if the
# sampling_frequency and duration are both integer values
idx_segment_end_float = (start_time + duration - hdf5_start_time) * sampling_frequency
idx_segment_end = int(idx_segment_end_float)
if idx_segment_end_float != idx_segment_end:
duration_real = (idx_segment_end_float - idx_segment_start) / sampling_frequency
cut.logger.warning(f'Due to numerical precision when selecting the segment, the last sample chosen corresponds to a duration of {duration_real} and not {duration} but we keep the latter value as we need an integer...')
idx_segment_end = int((start_time + duration - hdf5_start_time) * sampling_frequency)
# By convention due to fft routines used, `idx_segment_end` is not included
# in the selected time segment; meaning the segment we select is said to have
# lasted a duration = `duration` but its corresponding time_domain_series will give:
# `tds[-1] - tds[0] = duration - 1/sampling_frequency`
# `tds[-1] - tds[0] = duration - 1/sampling_frequency`.
desired_hdf5_strain = hdf5_strain[idx_segment_start:idx_segment_end]
# desired_hdf5_strain = hdf5_strain[idx_segment_start:idx_segment_end+1]
# Create the bilby strain object, set it from the time domain strain
# we selected from the hdf5 file and attach it to the ifo.
strain = bilby.gw.detector.strain_data.InterferometerStrainData(minimum_frequency=ifo.minimum_frequency)
strain.set_from_time_domain_strain(time_domain_strain=desired_hdf5_strain, sampling_frequency=sampling_frequency, duration=duration)
# strain.low_pass_filter(filter_freq=ifo.minimum_frequency)
ifo.strain_data = strain
ifo.strain_data.start_time = start_time
# import IPython; IPython.embed(); sys.exit()
# import matplotlib.pyplot as plt
# plt.plot(desired_hdf5_strain)
# plt.savefig(outdir + 'desired_hdf5_strain_4096.png')
# SET THE PSDs
# If a single file containing the, per say 3 PSDs, the following
# order is assumed column-wise: freq - H1 - L1 - V1.
if set_psd_from_single_file:
cut.logger.info(f'Setting PSD of {ifo.name} from {psds_file_path}')
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(psd_array=psds_event[:, i+1], frequency_array=psds_event[:, 0])
# If a path for each individual PSD has been given in the config.ini file,
# ie `psd_dict` is not an empty dict, we retrieve the path of the file.
elif psd_dict != {}:
psd_file_path = psd_dict[f'{ifo.name}']
freq_and_psd_array_ifo = np.loadtxt(psd_file_path)
cut.logger.info(f'Setting PSD of {ifo.name} from {psd_file_path}')
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(psd_array=freq_and_psd_array_ifo[:, 1], frequency_array=freq_and_psd_array_ifo[:, 0])
# Otherwise we compute the PSD from the time-domain strain file using the Welch method.
else:
cut.logger.info(f'Computing the PSD of {ifo.name} from the strain data in {hdf5_file_path}')
# 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.
......@@ -267,23 +282,30 @@ if __name__ == '__main__':
NFFT = int(4 * sampling_frequency)
psd_utils.set_interferometer_psd_from_hdf5_file(interferometer=ifo, hdf5_file_path=hdf5_file_path, signal_start_time=start_time, signal_duration=duration, psd_duration=psd_duration, psd_offset=psd_offset, filter_freq=filter_freq, NFFT=NFFT)
# I used to use a copy of `ifo.power_spectral_density_array` because
# at some point in Bilby the latter was hidding a function which was
# recomputing the PSD each time... it's been fixed by now i think.
# [TO BE MODIFIED]: Hence the two following could/should be removed
# from the entire code.
ifo.psd_array = ifo.power_spectral_density_array
ifo.inverse_psd_array = 1 / ifo.psd_array
# WAVEFORM GENERATION
# Set some waveform args which will be used when setting the
# waveform_generator
waveform_arguments = {}
waveform_arguments['minimum_frequency'] = min(ext_analysis_dict['minimum_frequency_ifos'].values())
waveform_arguments['waveform_approximant'] = ext_analysis_dict['approximant']
waveform_arguments['reference_frequency'] = ext_analysis_dict['reference_frequency']
waveform_arguments['maximum_frequency'] = ext_analysis_dict['maximum_frequency_injected_waveform']
# It's essentially a renaming of the variable here.
starting_point_parameters = trigger_parameters.copy()
del hdf5_strain, desired_hdf5_strain
# WARNING: this `inj_file` branch is certainly broken
# WARNING: this `inj_file` branch is certainly broken.
elif opts.inj_file is not None:
# SET THE INJECTION PARAMETERS AND THE ANALYSIS PARAMETERS
injection_parameters, ext_analysis_dict = set_inj.get_inj_parameters_and_analysis_dict(opts.inj_file, **config_dict['analysis'])
......@@ -351,9 +373,14 @@ if __name__ == '__main__':
elif opts.trigger_time is not None:
pass
# Set the
starting_point_parameters['geocent_duration'] = starting_point_parameters['geocent_time'] - start_time
# Note: Calling the `ifo.strain_data.frequency_domain_strain` for the first time sets the `ifo.strain_data.window_factor` to a value potentially != 1.
# Similar to ifo.psd_array, I used a copy of the `.frequency_domain_strain`
# to avoid some potentially unecessary back-end call that bilby was doing.
# [TO BE MODIFIED]: It should be checked if this is still necessary and if not
# remove `ifo.fd_strain` from the entire code.
for ifo in interferometers:
ifo.fd_strain = ifo.strain_data.frequency_domain_strain
......@@ -362,7 +389,10 @@ if __name__ == '__main__':
# SET THE LIKELIHOOD OBJECT, two different cases whether ROQ basis is used or not
# SET THE LIKELIHOOD AND PRIOR OBJECTS
# Two different cases whether ROQ basis is used or not.
# The ROQ case is probably broken by now but could be fixed quickly I reckon,
# and it treats some tricky steps.
if config_dict['analysis']['roq']:
if ext_analysis_dict['approximant'] != 'IMRPhenomPv2':
raise ValueError("Need to set approximant to IMRPhenomPv2 to work with the ROQ basis, approximant here is {}".format(ext_analysis_dict['approximant']))
......@@ -470,6 +500,7 @@ if __name__ == '__main__':
# write the weights to file so they can be loaded multiple times
likelihood.save_weights(weights_file_path_no_format, format=weight_format)
else:
# Set the frequency_source_model and priors according to BBH or BNS case.
prior_file = config_dict['parameters']['prior_file']
if config_dict['analysis'].get('frequency_domain_source_model', 'None') == 'lal_binary_black_hole':
frequency_domain_source_model = bilby.gw.source.lal_binary_black_hole
......@@ -479,7 +510,7 @@ if __name__ == '__main__':
frequency_domain_source_model = bilby.gw.source.lal_binary_neutron_star
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters
priors = bilby.gw.prior.BNSPriorDict(filename=prior_file)
# make waveform generator
# Set the waveform generator, needed for the likelihood object
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
duration=ext_analysis_dict['duration'], sampling_frequency=ext_analysis_dict['sampling_frequency'],
frequency_domain_source_model=frequency_domain_source_model,
......
......@@ -28,6 +28,9 @@ psd = 1
psds_single_file = ../__input_data/GWTC1_PSDs/GWTC1_GW170817_PSDs.dat
#psd_dict = {H1: /path/to/H1_PSD.dat, L1: /path/to/H1_PSD.dat, V1: /path/to/H1_PSD.dat}
phase_marginalization = True
# `lal_binary_black_hole` or `lal_binary_neutron_star`. If not set, `lal_binary_neutron_star` is selected by default.
# These are the waveform models defined in bilby.gw.source
frequency_domain_source_model = lal_binary_neutron_star
[parameters]
prior_file=../gw/prior_files/GW170817_LALInf_IMRPhenomD_NRTidal.prior
......
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