Skip to content
Snippets Groups Projects
Commit de2b3831 authored by Jade Powell's avatar Jade Powell
Browse files

changes to SN and SG examples

parent 8a9a0dec
No related branches found
No related tags found
1 merge request!52Burst
Pipeline #
...@@ -43,13 +43,16 @@ priors = dict() ...@@ -43,13 +43,16 @@ priors = dict()
# so for this example we will set almost all of the priors to be equall to their injected values. This implies the # 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 # 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. # not sample any parameter that has a delta-function prior.
for key in ['hrss', 'psi', 'ra', 'dec', 'geocent_time']: for key in ['psi', 'ra', 'dec', 'geocent_time']:
priors[key] = injection_parameters[key] priors[key] = injection_parameters[key]
# The above list does *not* include frequency and Q, which means those are the parameters # The above list does *not* include frequency and Q, which means those are the parameters
# that will be included in the sampler. If we do nothing, then the default priors get used. # that will be included in the sampler. If we do nothing, then the default priors get used.
priors['Q'] = tupak.prior.create_default_prior(name='Q') #priors['Q'] = tupak.prior.create_default_prior(name='Q')
priors['frequency'] = tupak.prior.create_default_prior(name='frequency') #priors['frequency'] = tupak.prior.create_default_prior(name='frequency')
priors['Q'] = tupak.prior.Uniform(2, 50, 'Q')
priors['frequency'] = tupak.prior.Uniform(30, 1000, 'frequency')
priors['hrss'] = tupak.prior.Uniform(1e-23, 1e-21, 'hrss')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator # Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator) likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
......
...@@ -23,7 +23,7 @@ np.random.seed(170801) ...@@ -23,7 +23,7 @@ np.random.seed(170801)
# We are going to inject a supernova waveform. We first establish a dictionary of parameters that # 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. # 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 = 60.0, ra = 1.375, 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) dec = -1.2108, geocent_time = 1126259642.413, psi= 2.659)
# Create the waveform_generator using a supernova source function # Create the waveform_generator using a supernova source function
...@@ -44,8 +44,8 @@ realPCs = np.loadtxt('SupernovaRealPCs.txt') ...@@ -44,8 +44,8 @@ realPCs = np.loadtxt('SupernovaRealPCs.txt')
imagPCs = np.loadtxt('SupernovaImagPCs.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. # 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, coeff1 = 0.1, coeff2 = 0.1, simulation_parameters = dict(realPCs=realPCs, imagPCs=imagPCs, pc_coeff1 = 0.1, pc_coeff2 = 0.1,
coeff3 = 0.1, coeff4 = 0.1, coeff5 = 0.1, luminosity_distance = 60.0, 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) ra = 1.375, dec = -1.2108, geocent_time = 1126259642.413, psi=2.659)
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration, waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
...@@ -62,14 +62,13 @@ priors = dict() ...@@ -62,14 +62,13 @@ priors = dict()
for key in ['psi', 'geocent_time']: for key in ['psi', 'geocent_time']:
priors[key] = injection_parameters[key] priors[key] = injection_parameters[key]
# The above list does *not* include frequency and Q, which means those are the parameters # don't use default for luminosity distance because we want kpc not Mpc
# that will be included in the sampler. If we do nothing, then the default priors get used. priors['luminosity_distance'] = tupak.prior.Uniform(2, 20, 'luminosity_distance')
priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance') priors['pc_coeff1'] = tupak.prior.Uniform(-100, 100, 'pc_coeff1')
priors['coeff1'] = tupak.prior.create_default_prior(name='coeff1') priors['pc_coeff2'] = tupak.prior.Uniform(-100, 100, 'pc_coeff2')
priors['coeff2'] = tupak.prior.create_default_prior(name='coeff2') priors['pc_coeff3'] = tupak.prior.Uniform(-100, 100, 'pc_coeff3')
priors['coeff3'] = tupak.prior.create_default_prior(name='coeff3') priors['pc_coeff4'] = tupak.prior.Uniform(-100, 100, 'pc_coeff4')
priors['coeff4'] = tupak.prior.create_default_prior(name='coeff4') priors['pc_coeff5'] = tupak.prior.Uniform(-100, 100, 'pc_coeff5')
priors['coeff5'] = tupak.prior.create_default_prior(name='coeff5')
priors['ra'] = tupak.prior.create_default_prior(name='ra') priors['ra'] = tupak.prior.create_default_prior(name='ra')
priors['dec'] = tupak.prior.create_default_prior(name='dec') priors['dec'] = tupak.prior.create_default_prior(name='dec')
......
...@@ -438,15 +438,7 @@ def create_default_prior(name): ...@@ -438,15 +438,7 @@ def create_default_prior(name):
'iota': Sine(name=name), 'iota': Sine(name=name),
'cos_iota': Uniform(name=name, minimum=-1, maximum=1), 'cos_iota': Uniform(name=name, minimum=-1, maximum=1),
'psi': Uniform(name=name, minimum=0, maximum=2 * np.pi), 'psi': Uniform(name=name, minimum=0, maximum=2 * np.pi),
'phase': Uniform(name=name, minimum=0, maximum=2 * np.pi), 'phase': Uniform(name=name, minimum=0, maximum=2 * np.pi)
'hrss': Uniform(name=name, minimum=1e-23, maximum=1e-21),
'Q': Uniform(name=name, minimum=2.0, maximum=50.0),
'frequency': Uniform(name=name, minimum=30.0, maximum=2000.0),
'coeff1': Uniform(name=name, minimum=-1.0, maximum=1.0),
'coeff2': Uniform(name=name, minimum=-1.0, maximum=1.0),
'coeff3': Uniform(name=name, minimum=-1.0, maximum=1.0),
'coeff4': Uniform(name=name, minimum=-1.0, maximum=1.0),
'coeff5': Uniform(name=name, minimum=-1.0, maximum=1.0)
} }
if name in default_priors.keys(): if name in default_priors.keys():
prior = default_priors[name] prior = default_priors[name]
......
...@@ -59,7 +59,7 @@ def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi ...@@ -59,7 +59,7 @@ def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi
fm = frequency_array - frequency fm = frequency_array - frequency
fp = frequency_array + frequency fp = frequency_array + frequency
h_plus = (hrss / np.sqrt(temp * (1+np.exp(-Q**2)))) * ((np.sqrt(np.pi)*tau)/2.0) * (np.exp(-fm**2 * np.pi**2 * tau**2) + np.exp(-fp**2 * pi**2 * tau**2)) h_plus = (hrss / np.sqrt(temp * (1+np.exp(-Q**2)))) * ((np.sqrt(np.pi)*tau)/2.0) * (np.exp(-fm**2 * np.pi**2 * tau**2) + np.exp(-fp**2 * np.pi**2 * tau**2))
h_cross = -1j*(hrss / np.sqrt(temp * (1-np.exp(-Q**2)))) * ((np.sqrt(np.pi)*tau)/2.0) * (np.exp(-fm**2 * np.pi**2 * tau**2) - np.exp(-fp**2 * np.pi**2 * tau**2)) h_cross = -1j*(hrss / np.sqrt(temp * (1-np.exp(-Q**2)))) * ((np.sqrt(np.pi)*tau)/2.0) * (np.exp(-fm**2 * np.pi**2 * tau**2) - np.exp(-fp**2 * np.pi**2 * tau**2))
...@@ -71,13 +71,13 @@ def supernova(frequency_array, realPCs, imagPCs, file_path, luminosity_distance, ...@@ -71,13 +71,13 @@ def supernova(frequency_array, realPCs, imagPCs, file_path, luminosity_distance,
realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(file_path, usecols = (0,1,2,3), unpack=True) realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(file_path, usecols = (0,1,2,3), unpack=True)
# waveform in file at 10kpc # waveform in file at 10kpc
scaling = 1e-2 * (10.0 / luminosity_distance) scaling = 1e-3 * (10.0 / luminosity_distance)
h_plus = scaling * (realhplus + 1.0j*imaghplus) h_plus = scaling * (realhplus + 1.0j*imaghplus)
h_cross = scaling * (realhcross + 1.0j*imaghcross) h_cross = scaling * (realhcross + 1.0j*imaghcross)
return {'plus': h_plus, 'cross': h_cross} return {'plus': h_plus, 'cross': h_cross}
def supernova_pca_model(frequency_array, realPCs, imagPCs, coeff1, coeff2, coeff3, coeff4, coeff5, luminosity_distance, ra, dec, geocent_time, psi): def supernova_pca_model(frequency_array, realPCs, imagPCs, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5, luminosity_distance, ra, dec, geocent_time, psi):
""" Supernova signal model """ """ Supernova signal model """
pc1 = realPCs[:,0] + 1.0j*imagPCs[:,0] pc1 = realPCs[:,0] + 1.0j*imagPCs[:,0]
...@@ -87,10 +87,10 @@ def supernova_pca_model(frequency_array, realPCs, imagPCs, coeff1, coeff2, coeff ...@@ -87,10 +87,10 @@ def supernova_pca_model(frequency_array, realPCs, imagPCs, coeff1, coeff2, coeff
pc5 = realPCs[:,4] + 1.0j*imagPCs[:,5] pc5 = realPCs[:,4] + 1.0j*imagPCs[:,5]
# file at 10kpc # file at 10kpc
scaling = 1e-22 * (10.0 / luminosity_distance) scaling = 1e-23 * (10.0 / luminosity_distance)
h_plus = scaling * (coeff1*pc1 + coeff2*pc2 + coeff3*pc3 + coeff4*pc4 + coeff5*pc5) h_plus = scaling * (pc_coeff1*pc1 + pc_coeff2*pc2 + pc_coeff3*pc3 + pc_coeff4*pc4 + pc_coeff5*pc5)
h_cross = scaling * (coeff1*pc1 + coeff2*pc2 + coeff3*pc3 + coeff4*pc4 + coeff5*pc5) h_cross = scaling * (pc_coeff1*pc1 + pc_coeff2*pc2 + pc_coeff3*pc3 + pc_coeff4*pc4 + pc_coeff5*pc5)
return {'plus': h_plus, 'cross': h_cross} return {'plus': h_plus, 'cross': h_cross}
......
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