diff --git a/examples/supernova_example/supernova_example.py b/examples/supernova_example/supernova_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..d2b0a8751fddbc4cc886f91917f4514e13d437b0
--- /dev/null
+++ b/examples/supernova_example/supernova_example.py
@@ -0,0 +1,96 @@
+#!/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)
+
+"""
+from __future__ import division, print_function
+import tupak
+import numpy as np
+
+# Set the duration and sampling frequency of the data segment that we're going to inject the signal into
+time_duration = 3.
+sampling_frequency = 4096.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'supernova'
+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 = 60.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)
+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
+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']]
+
+# 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, coeff1 = 0.1, coeff2 = 0.1,
+ coeff3 = 0.1, coeff4 = 0.1, coeff5 = 0.1, luminosity_distance = 60.0,
+ ra = 1.375, dec = -1.2108, geocent_time = 1126259642.413, psi=2.659)
+
+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
+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]
+
+# 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.
+priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance')
+priors['coeff1'] = tupak.prior.create_default_prior(name='coeff1')
+priors['coeff2'] = tupak.prior.create_default_prior(name='coeff2')
+priors['coeff3'] = tupak.prior.create_default_prior(name='coeff3')
+priors['coeff4'] = tupak.prior.create_default_prior(name='coeff4')
+priors['coeff5'] = tupak.prior.create_default_prior(name='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=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+#result.plot_corner()
+print(result)
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tupak/prior.py b/tupak/prior.py
index 00e5b417e39cb7a6969dd96d1d42188a5c149bab..9324442ecd8be9d4e4d7f2024b5ac0e3ba5e5282 100644
--- a/tupak/prior.py
+++ b/tupak/prior.py
@@ -441,7 +441,12 @@ def create_default_prior(name):
'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)
+ '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():
prior = default_priors[name]
diff --git a/tupak/source.py b/tupak/source.py
index 5c1a50f0496d50038d93fe4c9000d9ac3caefe66..5cdf3e43efd672de7a5f29e2b2aa52cf235e4842 100644
--- a/tupak/source.py
+++ b/tupak/source.py
@@ -53,7 +53,6 @@ def lal_binary_black_hole(
def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi):
- pi = 3.14159
tau = Q / (np.sqrt(2.0)*np.pi*frequency)
temp = Q / (4.0*np.sqrt(np.pi)*frequency)
t = geocent_time
@@ -62,35 +61,30 @@ def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi
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_cross = -1j*(hrss / np.sqrt(temp * (1-np.exp(-Q**2)))) * ((np.sqrt(pi)*tau)/2.0) * (np.exp(-fm**2 * pi**2 * tau**2) - np.exp(-fp**2 * 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))
return{'plus': h_plus, 'cross': h_cross}
-def supernova(frequency_array, file_path, luminosity_distance, ra, dec, geocent_time, psi):
+def supernova(frequency_array, realPCs, imagPCs, file_path, luminosity_distance, ra, dec, geocent_time, psi):
""" A supernova NR simulation for injections """
-# realhplus, imaghplus = np.loadtxt(file_path , usecols = (0,1), unpack=True)
- realhplus, imaghplus, realhcross, imaghcross = np.loadtxt('MuellerL15_example_inj.txt', 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
- scaling = 10.0 / luminosity_distance
+ scaling = 1e-2 * (10.0 / luminosity_distance)
h_plus = scaling * (realhplus + 1.0j*imaghplus)
h_cross = scaling * (realhcross + 1.0j*imaghcross)
return {'plus': h_plus, 'cross': h_cross}
-def supernova_pca_model(frequency_array, coeff1, coeff2, coeff3, coeff4, coeff5, luminosity_distance, ra, dec, geocent_time, psi):
+def supernova_pca_model(frequency_array, realPCs, imagPCs, coeff1, coeff2, coeff3, coeff4, coeff5, luminosity_distance, ra, dec, geocent_time, psi):
""" Supernova signal model """
- # this is slow reading in the file every time
- realpc1, realpc2, realpc3, realpc4, realpc5 = np.loadtxt('SupernovaRealPCs.txt', usecols = (0,1,2,3,4), unpack=True)
- imagpc1, imagpc2, imagpc3, imagpc4, imagpc5 = np.loadtxt('SupernovaImagPCs.txt', usecols = (0,1,2,3,4), unpack=True)
-
- pc1 = realpc1 + 1.0j*imagpc1
- pc2 = realpc2 + 1.0j*imagpc2
- pc3 = realpc3 + 1.0j*imagpc3
- pc4 = realpc4 + 1.0j*imagpc4
- pc5 = realpc5 + 1.0j*imagpc5
+ pc1 = realPCs[:,0] + 1.0j*imagPCs[:,0]
+ pc2 = realPCs[:,1] + 1.0j*imagPCs[:,1]
+ pc3 = realPCs[:,2] + 1.0j*imagPCs[:,2]
+ pc4 = realPCs[:,3] + 1.0j*imagPCs[:,3]
+ pc5 = realPCs[:,4] + 1.0j*imagPCs[:,5]
# file at 10kpc
scaling = 1e-22 * (10.0 / luminosity_distance)