Skip to content
Snippets Groups Projects
Commit 47cec64a authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Add time prior and increase number of live points

parent 9c82e28c
No related branches found
No related tags found
1 merge request!52Burst
Pipeline #
#!/bin/python
"""
Tutorial to demonstrate running parameter estimation/model selection on an NR supernova injected signal.
Signal model is made by applying PCA to a set of supernova waveforms. The first few PCs are then linearly
combined with a scale factor. (See https://arxiv.org/pdf/1202.3256.pdf)
Tutorial to demonstrate running parameter estimation/model selection on an NR
supernova injected signal. Signal model is made by applying PCA to a set of
supernova waveforms. The first few PCs are then linearly combined with a scale
factor. (See https://arxiv.org/pdf/1202.3256.pdf)
"""
from __future__ import division, print_function
......@@ -21,49 +23,51 @@ tupak.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(170801)
# We are going to inject a supernova waveform. We first establish a dictionary of parameters that
# includes all of the different waveform parameters. It will read in a signal to inject from a txt file.
injection_parameters = dict(file_path = 'MuellerL15_example_inj.txt', luminosity_distance = 10.0, ra = 1.375,
dec = -1.2108, geocent_time = 1126259642.413, psi= 2.659)
# We are going to inject a supernova waveform. We first establish a dictionary
# of parameters that includes all of the different waveform parameters. It will
# read in a signal to inject from a txt file.
injection_parameters = dict(file_path='MuellerL15_example_inj.txt',
luminosity_distance=10.0, ra=1.375,
dec=-1.2108, geocent_time=1126259642.413,
psi=2.659)
# Create the waveform_generator using a supernova source function
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.supernova,
parameters=injection_parameters)
waveform_generator = tupak.waveform_generator.WaveformGenerator(
time_duration=time_duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.supernova,
parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers. In this case we'll use three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1),
# and Virgo (V1)). These default to their design sensitivity
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)). These default to
# their design sensitivity
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, time_duration=time_duration,
sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']]
name, injection_polarizations=hf_signal,
injection_parameters=injection_parameters, time_duration=time_duration,
sampling_frequency=sampling_frequency, outdir=outdir)
for name in ['H1', 'L1', 'V1']]
# read in from a file the PCs used to create the signal model.
realPCs = np.loadtxt('SupernovaRealPCs.txt')
imagPCs = np.loadtxt('SupernovaImagPCs.txt')
# now we have to do the waveform_generator again because the signal model is not the same as the injection in this case.
simulation_parameters = dict(realPCs=realPCs, imagPCs=imagPCs, pc_coeff1 = 0.1, pc_coeff2 = 0.1,
pc_coeff3 = 0.1, pc_coeff4 = 0.1, pc_coeff5 = 0.1, luminosity_distance = 10.0,
ra = 1.375, dec = -1.2108, geocent_time = 1126259642.413, psi=2.659)
# Now we make another waveform_generator because the signal model is
# not the same as the injection in this case.
simulation_parameters = dict(
realPCs=realPCs, imagPCs=imagPCs)
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.supernova_pca_model,
parameters=simulation_parameters)
search_waveform_generator = tupak.waveform_generator.WaveformGenerator(
time_duration=time_duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.supernova_pca_model,
parameters=simulation_parameters)
# Set up prior, which is a dictionary
# Set up prior
priors = dict()
# By default we will sample all terms in the signal models. However, this will take a long time for the calculation,
# so for this example we will set almost all of the priors to be equall to their injected values. This implies the
# prior is a delta function at the true, injected value. In reality, the sampler implementation is smart enough to
# not sample any parameter that has a delta-function prior.
for key in ['psi', 'geocent_time']:
priors[key] = injection_parameters[key]
# don't use default for luminosity distance because we want kpc not Mpc
priors['luminosity_distance'] = tupak.prior.Uniform(2, 20, 'luminosity_distance')
priors['luminosity_distance'] = tupak.prior.Uniform(
2, 20, 'luminosity_distance')
priors['pc_coeff1'] = tupak.prior.Uniform(-1, 1, 'pc_coeff1')
priors['pc_coeff2'] = tupak.prior.Uniform(-1, 1, 'pc_coeff2')
priors['pc_coeff3'] = tupak.prior.Uniform(-1, 1, 'pc_coeff3')
......@@ -71,13 +75,20 @@ priors['pc_coeff4'] = tupak.prior.Uniform(-1, 1, 'pc_coeff4')
priors['pc_coeff5'] = tupak.prior.Uniform(-1, 1, 'pc_coeff5')
priors['ra'] = tupak.prior.create_default_prior(name='ra')
priors['dec'] = tupak.prior.create_default_prior(name='dec')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=100,
outdir=outdir, label=label)
priors['geocent_time'] = tupak.prior.Uniform(
injection_parameters['geocent_time'] - 1,
injection_parameters['geocent_time'] + 1,
'geocent_time')
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveoform generator
likelihood = tupak.GravitationalWaveTransient(
interferometers=IFOs, waveform_generator=search_waveform_generator)
# Run sampler.
result = tupak.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
outdir=outdir, label=label)
# make some plots of the outputs
result.plot_corner()
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment