Merge branch 'burst' into 'master'

Burst

See merge request Monash/tupak!52

examples/other_examples/sine_gaussian_example.py

+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation on a sine gaussian injected signal.
+
+"""
+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 = 4.
+sampling_frequency = 2048.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'sine_gaussian'
+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 sine gaussian waveform.  We first establish a dictionary of parameters that
+# includes all of the different waveform parameters
+injection_parameters = dict(hrss = 1e-22, Q = 5.0, frequency = 200.0, ra = 1.375, dec = -1.2108, 
+                             geocent_time = 1126259642.413, psi= 2.659)
+
+# Create the waveform_generator using a sine Gaussian source function
+waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
+                                                                sampling_frequency=sampling_frequency,
+                                                                frequency_domain_source_model=tupak.source.sinegaussian,
+                                                                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']]
+
+# 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', 'ra', 'dec', '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['Q'] = tupak.prior.create_default_prior(name='Q')
+#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
+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)
+
+
+
+
+
+
+
+
+
+
+

examples/supernova_example/supernova_example.py

+#!/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=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)
+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 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)
+
+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
+priors = dict()
+for key in ['psi', 'geocent_time']:
+    priors[key] = injection_parameters[key]
+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')
+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')
+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()
+print(result)
+
+
+
+
+
+
+
+
+
+
+

tupak/source.py

 from __future__ import division, print_function
 
 import logging
+import numpy as np
 
 try:
     import lalsimulation as lalsim
@@ -59,6 +60,67 @@ def lal_binary_black_hole(
     return {'plus': h_plus, 'cross': h_cross}
 
 
+def sinegaussian(
+        frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi):
+
+    tau  = Q / (np.sqrt(2.0)*np.pi*frequency)
+    temp = Q / (4.0*np.sqrt(np.pi)*frequency)
+    t = geocent_time
+    fm = 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 * 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)))
+
+    return{'plus': h_plus, 'cross': h_cross}
+
+
+def supernova(
+        frequency_array, realPCs, imagPCs, file_path, luminosity_distance, ra,
+        dec, geocent_time, psi):
+    """ A supernova NR simulation for injections """
+
+    realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(
+        file_path, usecols=(0, 1, 2, 3), unpack=True)
+
+    # waveform in file at 10kpc
+    scaling = 1e-3 * (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, realPCs, imagPCs, pc_coeff1, pc_coeff2, pc_coeff3,
+        pc_coeff4, pc_coeff5, luminosity_distance, ra, dec, geocent_time, psi):
+    """ Supernova signal model """
+
+    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-23 * (10.0 / luminosity_distance)
+
+    h_plus = scaling * (pc_coeff1*pc1 + pc_coeff2*pc2 + pc_coeff3*pc3
+                        + pc_coeff4*pc4 + pc_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}
+
+
+
 #class BinaryNeutronStarMergerNumericalRelativity:
 #    """Loads in NR simulations of BNS merger
 #    takes parameters mean_mass, mass_ratio and equation_of_state, directory_path