Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Chris vitc updates #116

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Chris vitc updates #116

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bd2add1
All changes made before switching to updates from Chris
Jun 8, 2021
1e30cd7
latest version of vitamin for LLO runs
Jul 6, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 53 additions & 15 deletions vitamin_c/gen_benchmark_pe.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
import time
import h5py
from scipy.ndimage.interpolation import shift
from scipy import signal
import argparse
from gwpy.timeseries import TimeSeries

Expand Down Expand Up @@ -122,8 +123,11 @@ def gen_template(duration,
if sampling_frequency>4096:
print('EXITING: bilby doesn\'t seem to generate noise above 2048Hz so lower the sampling frequency')
exit(0)

real_data_psd = True
# compute the number of time domain samples
if real_data_psd:
# Always using a buffer of an additional 1s for real data psd
duration = 2
Nt = int(sampling_frequency*duration)

# define the start time of the timeseries
Expand Down Expand Up @@ -164,13 +168,20 @@ def gen_template(duration,
if len(psd_files) > 0:
type_psd = psd_files[0].split('/')[-1].split('_')[-1].split('.')[0]
for int_idx,ifo in enumerate(ifos):
if type_psd == 'psd':
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(psd_file=psd_files[0])
elif type_psd == 'asd':
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(asd_file=psd_files[0])
if real_data_psd:
#psd_array = np.array(h5py.File(psd_files[0], 'r')['%s_psd_array' % ifo.name])
psd_array = np.load('psd_estimates/O1_psd_%s_2sTS.npy' % ifo.name)
psd_frequencies = np.linspace(0, len(psd_array), num=len(psd_array))
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(
frequency_array=psd_frequencies, psd_array=psd_array)
else:
print('Could not determine whether psd or asd ...')
exit()
if type_psd == 'psd':
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(psd_file=psd_files[0])
elif type_psd == 'asd':
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(asd_file=psd_files[0])
else:
print('Could not determine whether psd or asd ...')
exit()

ifos.set_strain_data_from_power_spectral_densities(
sampling_frequency=sampling_frequency, duration=duration,
Expand All @@ -184,7 +195,29 @@ def gen_template(duration,
whitened_signal_td_all = []
whitened_h_td_all = []
# iterate over ifos
for i in range(len(pars['det'])):

"""
def real_data_whitening(time_series, estimated_psd, duration, sampling_frequency):
tukey_win = signal.windows.tukey(len(time_series), alpha=0.5)
time_series = tukey_win*time_series
# Transform time series to frequency domain
hf = np.fft.rfft(np.array(time_series))
# Whiten
real_noise = hf / np.sqrt(estimated_psd) * norm
# irfft
noise_sample = np.fft.irfft(real_noise, n=Nt)
# Apply butterworth high-pass filter to supress low frequency noise at 10Hz
fband = [40, 950]
b, a = signal.butter(4, [fband[0]/sampling_frequency, fband[1]/sampling_frequency], 'band')
normalization = np.sqrt((fband[1]-fband[0])/(sampling_frequency/2))
noise_sample_before_butterworth = np.array(noise_sample)
noise_sample = signal.filtfilt(b, a, np.array(noise_sample)) / normalization
real_noise_data = real_noise_data[int( ((duration*0.5)-0.5)*sampling_frequency ):int( ((duration*0.5)+0.5)*sampling_frequency )]
return real_noise_data
"""

for i_idx,i in enumerate(range(len(pars['det']))):

# get frequency domain noise-free signal at detector
signal_fd = ifos[i].get_detector_response(freq_signal, injection_parameters)

Expand All @@ -203,11 +236,16 @@ def gen_template(duration,

# inverse FFT noisy signal back to time domain and normalise
whitened_h_td = np.sqrt(2.0*Nt)*np.fft.irfft(whitened_h_fd)


if real_data_psd:
whitened_h_td = whitened_h_td[int( ((duration*0.5)-0.5)*sampling_frequency ):int( ((duration*0.5)+0.5)*sampling_frequency )]
whitened_signa_td = whitened_signal_td[int( ((duration*0.5)-0.5)*sampling_frequency ):int( ((duration*0.5)+0.5)*sampling_frequency )]

whitened_h_td_all.append([whitened_h_td])
whitened_signal_td_all.append([whitened_signal_td])

print('... Whitened signals')
exit()
return np.squeeze(np.array(whitened_signal_td_all),axis=1),np.squeeze(np.array(whitened_h_td_all),axis=1),injection_parameters,ifos,waveform_generator

def run(sampling_frequency=256.0,
Expand Down Expand Up @@ -313,12 +351,12 @@ def run(sampling_frequency=256.0,
priors['a_2'] = fixed_vals['a_2']

if np.any([r=='tilt_1' for r in rand_pars]):
priors['tilt_1'] = bilby.gw.prior.Uniform(name='tilt_1', minimum=bounds['tilt_1_min'], maximum=bounds['tilt_1_max'])
priors['tilt_1'] = bilby.core.prior.Sine(name='tilt_1')
else:
priors['tilt_1'] = fixed_vals['tilt_1']

if np.any([r=='tilt_2' for r in rand_pars]):
priors['tilt_2'] = bilby.gw.prior.Uniform(name='tilt_2', minimum=bounds['tilt_2_min'], maximum=bounds['tilt_2_max'])
priors['tilt_2'] = bilby.core.prior.Sine(name='tilt_2')
else:
priors['tilt_2'] = fixed_vals['tilt_2']

Expand All @@ -338,7 +376,7 @@ def run(sampling_frequency=256.0,
priors['ra'] = fixed_vals['ra']

if np.any([r=='dec' for r in rand_pars]):
pass
priors['dec'] = bilby.core.prior.Cosine(name='dec')
else:
priors['dec'] = fixed_vals['dec']

Expand All @@ -348,7 +386,7 @@ def run(sampling_frequency=256.0,
priors['psi'] = fixed_vals['psi']

if np.any([r=='theta_jn' for r in rand_pars]):
pass
priors['theta_jn'] = bilby.core.prior.Sine(name='theta_jn')
else:
priors['theta_jn'] = fixed_vals['theta_jn']

Expand Down Expand Up @@ -432,7 +470,7 @@ def run(sampling_frequency=256.0,

# Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveform generator
phase_marginalization=True
phase_marginalization=False
likelihood = bilby.gw.GravitationalWaveTransient(
interferometers=ifos, waveform_generator=waveform_generator, phase_marginalization=phase_marginalization,
priors=priors)
Expand Down Expand Up @@ -468,7 +506,7 @@ def run(sampling_frequency=256.0,
# Run sampler dynesty 1 sampler

result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, nact=50, npool=8, dlogz=27.0,
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, nact=50, npool=8, dlogz=0.1,
injection_parameters=injection_parameters, outdir=out_dir+'_'+samplers[-1], label=label,
save='hdf5', plot=True)
run_endt = time.time()
Expand Down
26 changes: 13 additions & 13 deletions vitamin_c/make_params_files.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@
'__definition__luminosity_distance': 'luminosity distance range'}

# arbitrary value for normalization of timeseries (Don't change this)
y_normscale = 16.638832624721797
y_normscale = 36.0

##########################
# Main tunable variables
Expand All @@ -94,7 +94,7 @@
#inf_pars = ["mass_1","mass_2","luminosity_distance","geocent_time","phase","theta_jn","psi","ra","dec","a_1","a_2","tilt_1","tilt_2","phi_12","phi_jl"]
batch_size = int(512)
weight_init = 'xavier'
n_modes=64; n_modes_q=1
n_modes=32; n_modes_q=1
initial_training_rate=1e-4
batch_norm=False

Expand Down Expand Up @@ -131,25 +131,25 @@
#############################
# optional tunable variables
#############################
run_label = 'vitamin_c_run19'#'demo_%ddet_%dpar_%dHz_hour_angle_with_late_kl_start' % (len(det),len(rand_pars),ndata)
gpu_num = 5
run_label = 'vitamin_c_run53'#'vitamin_c_run31_fast_testing'#'demo_%ddet_%dpar_%dHz_hour_angle_with_late_kl_start' % (len(det),len(rand_pars),ndata)
gpu_num = 2

# 1024 Hz label
#bilby_results_label = 'weichangfeng_theta_jn_issue'
# 256 Hz label
bilby_results_label = '1024Hz_full_15par'

r = 10 # 251
r = 2 # 251
pe_test_num = 256
tot_dataset_size = int(1e7)

tset_split = int(1e3)
save_interval = int(1000)
plot_interval = int(1000)
num_iterations=int(25000)+1 # 12500 for Green et al.
num_iterations=int(100000)+1 # 12500 for Green et al.
ref_geocent_time=1126259642.5
load_chunk_size = int(2e4)
samplers=['vitamin','dynesty','ptemcee','cpnest']
samplers=['vitamin','dynesty','ptemcee']#,'cpnest','emcee']
val_dataset_size = int(1e3)

# Directory variables
Expand All @@ -162,8 +162,8 @@
#pe_dir='/home/hunter.gabbard/CBC/public_VItamin/provided_models/vitamin_b/vitamin_b/test_sets/1024_khz_spins_included_15par/test'

# default training/testing directories
train_set_dir='./training_sets_realnoise_%ddet_%dpar_%dHz/tset_tot-%d_split-%d' % (len(det),len(rand_pars),ndata,tot_dataset_size,tset_split)
val_set_dir='./validation_sets_realnoise_%ddet_%dpar_%dHz/tset_tot-%d_split-%d' % (len(det),len(rand_pars),ndata,val_dataset_size,tset_split)
train_set_dir='./training_sets_TiltSineDist_%ddet_%dpar_%dHz/tset_tot-%d_split-%d' % (len(det),len(rand_pars),ndata,tot_dataset_size,tset_split)
val_set_dir='./validation_sets_TiltSineDist_%ddet_%dpar_%dHz/tset_tot-%d_split-%d' % (len(det),len(rand_pars),ndata,val_dataset_size,tset_split)
test_set_dir = './test_sets/%s/test_waveforms' % bilby_results_label
pe_dir='./test_sets/%s/test' % bilby_results_label

Expand Down Expand Up @@ -203,7 +203,7 @@ def get_params():
__definition__use_real_det_noise='If True, use real detector noise around reference time',
real_noise_time_range = [1126051217, 1137254417],
__definition__real_noise_time_range = 'time range to produce real noise samples over',
use_real_events=['GW150914'], # e.g. 'GW150914'
use_real_events=False, # e.g. 'GW150914'
__definition__use_real_events='list containing real ligo events to be analyzed',
convert_to_hour_angle = False,
__definition__convert_to_hour_angle='If True, convert RA to hour angle during trianing and testing',
Expand Down Expand Up @@ -374,11 +374,11 @@ def get_params():
__definition__figure_sampler_names='matplotlib figure sampler labels (e.g. [ptemcee,CPNest,emcee])',
y_normscale = y_normscale,
__definition__y_normscale='arbitrary normalization factor on all time series waveforms (helps convergence in training)',
gauss_pars=['luminosity_distance','geocent_time','theta_jn','a_1','a_2','tilt_1','tilt_2','ra','dec','phase','psi','phi_12','phi_jl'],
gauss_pars=['luminosity_distance','geocent_time','theta_jn','a_1','a_2','tilt_1','tilt_2','ra','dec'],
__definition__gauss_pars='parameters that require a truncated gaussian distribution',
vonmise_pars=[],
vonmise_pars=['phase','phi_12','phi_jl','psi'],
__definition__vonmises_pars='parameters that get wrapped on the 1D parameter',
sky_pars=[],
sky_pars=['ra','dec'],
__definition__sky_pars='sky parameters',
det=det,
__definition__det='LIGO detectors to perform analysis on (default is 3detector H1,L1,V1)',
Expand Down
41 changes: 21 additions & 20 deletions vitamin_c/params_files/params.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,7 @@
1137254417
],
"__definition__real_noise_time_range": "time range to produce real noise samples over",
"use_real_events": [
"GW150914"
],
"use_real_events": false,
"__definition__use_real_events": "list containing real ligo events to be analyzed",
"convert_to_hour_angle": false,
"__definition__convert_to_hour_angle": "If True, convert RA to hour angle during trianing and testing",
Expand Down Expand Up @@ -95,19 +93,19 @@
"__definition__load_plot_data": "Use plotting data which has already been generated",
"doPE": true,
"__definition__doPE": "if True then do bilby PE when generating new testing samples",
"gpu_num": 5,
"gpu_num": 2,
"__definition__gpu_num": "gpu number run is running on",
"ndata": 1024,
"__definition__ndata": "sampling frequency",
"run_label": "vitamin_c_run19",
"run_label": "vitamin_c_run53",
"__definition__run_label": "label of run",
"bilby_results_label": "1024Hz_full_15par",
"__definition__bilby_result_label": "label given to bilby results directory",
"tot_dataset_size": 10000000,
"__definition__tot_dataset_size": "total number of training samples available to use",
"tset_split": 1000,
"__definition__tset_split": "number of training samples in each training data file",
"plot_dir": "/home/hunter.gabbard/public_html/CBC/chris_dec2020_vitamin/vitamin_c_run19",
"plot_dir": "/home/hunter.gabbard/public_html/CBC/chris_dec2020_vitamin/vitamin_c_run53",
"__definition__plot_dir": "output directory to save results plots",
"hyperparam_optim_stop": 500000,
"__definition__hyperparam_optim_stop": "stopping iteration of hyperparameter optimizer per call",
Expand All @@ -121,7 +119,7 @@
"__definition__weight_init": "[xavier,VarianceScaling,Orthogonal]. Network model weight initialization (default is xavier)",
"n_samples": 8000,
"__definition__n_samples": "number of posterior samples to save per reconstruction upon inference (default 3000)",
"num_iterations": 25001,
"num_iterations": 100001,
"__definition__num_iterations": "total number of iterations before ending training of model",
"initial_training_rate": 0.0001,
"__definition__initial_training_rate": "initial training rate for ADAM optimiser inference model (inverse reconstruction)",
Expand All @@ -131,7 +129,7 @@
"__definition__batch_norm": "if true, do batch normalization in all layers of neural network",
"report_interval": 500,
"__definition__report_interval": "interval at which to save objective function values and optionally print info during training",
"n_modes": 64,
"n_modes": 32,
"__definition__n_modes": "number of modes in Gaussian mixture model (ideal 7, but may go higher/lower)",
"n_filters_r1": [
32,
Expand Down Expand Up @@ -270,7 +268,7 @@
"__definition__n_weights_q": "number fully-connected neurons in layers of encoders and decoders in the q model",
"duration": 1.0,
"__definition__duration": "length of training/validation/test sample time series in seconds",
"r": 10,
"r": 2,
"__definition__r": "number of GW timeseries to use for testing.",
"rand_pars": [
"mass_1",
Expand Down Expand Up @@ -334,9 +332,9 @@
"dec"
],
"__definition__inf_pars": "parameters to infer",
"train_set_dir": "./training_sets_realnoise_3det_15par_1024Hz/tset_tot-10000000_split-1000",
"train_set_dir": "./training_sets_TiltSineDist_3det_15par_1024Hz/tset_tot-10000000_split-1000",
"__definition__train_set_dir": "location of training set directory",
"val_set_dir": "./validation_sets_realnoise_3det_15par_1024Hz/tset_tot-1000_split-1000",
"val_set_dir": "./validation_sets_TiltSineDist_3det_15par_1024Hz/tset_tot-1000_split-1000",
"__definition__val_set_dir": "location of validation set directory",
"test_set_dir": "./test_sets/1024Hz_full_15par/test_waveforms",
"__definition__test_set_dir": "location of test set directory waveforms",
Expand All @@ -347,8 +345,7 @@
"samplers": [
"vitamin",
"dynesty",
"ptemcee",
"cpnest"
"ptemcee"
],
"__definition__samplers": "Samplers to use when comparing ML results (vitamin is ML approach) dynesty,ptemcee,cpnest,emcee",
"figure_sampler_names": [
Expand All @@ -359,7 +356,7 @@
"Emcee"
],
"__definition__figure_sampler_names": "matplotlib figure sampler labels (e.g. [ptemcee,CPNest,emcee])",
"y_normscale": 16.638832624721797,
"y_normscale": 36.0,
"__definition__y_normscale": "arbitrary normalization factor on all time series waveforms (helps convergence in training)",
"gauss_pars": [
"luminosity_distance",
Expand All @@ -370,16 +367,20 @@
"tilt_1",
"tilt_2",
"ra",
"dec",
"dec"
],
"__definition__gauss_pars": "parameters that require a truncated gaussian distribution",
"vonmise_pars": [
"phase",
"psi",
"phi_12",
"phi_jl"
"phi_jl",
"psi"
],
"__definition__gauss_pars": "parameters that require a truncated gaussian distribution",
"vonmise_pars": [],
"__definition__vonmises_pars": "parameters that get wrapped on the 1D parameter",
"sky_pars": [],
"sky_pars": [
"ra",
"dec"
],
"__definition__sky_pars": "sky parameters",
"det": [
"H1",
Expand Down
Loading