Skip to content
Snippets Groups Projects
Commit 7e05d8a6 authored by Moritz Huebner's avatar Moritz Huebner
Browse files

Merge branch 'master' of https://git.ligo.org/Monash/tupak

parents 753d7445 47b9d7d4
No related branches found
No related tags found
No related merge requests found
Showing
with 197 additions and 110 deletions
......@@ -23,15 +23,15 @@ exitcode-jessie:
- pip install coverage-badge
- python setup.py install
- coverage erase
- coverage run --include=tupak/* -a test/detector_tests.py
- coverage run --include=tupak/* -a test/prior_tests.py
- coverage run --include=tupak/* -a test/sampler_tests.py
- coverage run --include=tupak/* -a test/waveform_generator_tests.py
- coverage html
- coverage-badge -o coverage_badge.svg -f
- coverage run --source=tupak/ -a test/detector_tests.py
- coverage run --source=tupak/ -a test/prior_tests.py
- coverage run --source=tupak/ -a test/sampler_tests.py
- coverage run --source=tupak/ -a test/waveform_generator_tests.py
- coverage run --omit=* -a test/example_tests.py
- coverage run --omit=* -a test/noise_realisation_tests.py
- coverage run --omit=* -a test/tests.py
- coverage html
- coverage-badge -o coverage_badge.svg -f
# Make the documentation
- pip install -r docs/requirements.txt
......
.. _hyperparameters:
================
Hyper-parameters
================
This short guide will illustrate how to estimate hyper-parameters from
posterior samples using :code:`tupak` using a simplistic problem.
Setting up the problem
----------------------
We are given three data sets labelled by a, b, and c. Each data set consists of
:math:`N` observations of a variable :math:`y` taken at a dependent variable
:math:`x`. Notationally, we can write these three data sets as :math:`(y_i^a,
x_i^a),(y_i^b, x_i^b),(y_i^c, x_i^c)` where :math:`i \in [0, N]` labels the
indexes of each data set.
Plotting the data, we see that all three look like they could be modelled by
a linear function with some slope and some intercept:
.. image:: images/hyper_parameter_data.png
Given any individual data set, you could write down a linear model :math:`y(x)
= c_0 + c_1 x` and infer the intercept and gradient. For example, given the
a-data set, you could calculate :math:`P(c_0|y_i^a, x_i^a)` by fitting the
model to the data. Here is a figure demonstrating the posteriors on the
intercept, for the three data sets given above
.. image:: images/hyper_parameter_combined_posteriors.png
While the data looks noise, the overlap in the :math:`c_0` posteriors might
(especially given context about how the data was produced and the physical
setting) make you believe all data sets share a common intercept. How would
you go about estimating this? You could just take the mean of the means of the
posterior and that would give you a pretty good estimate. However, we'll now
discuss how to do this with hyper parameters.
Understanding the population using hyperparameters
--------------------------------------------------
We first need to define our hyperparameterized model. In this case, it is
as simple as
.. math::
c_0 \sim \mathcal{N}(\mu_{c_0}, \sigma_{c_0})
That is, we model the population :math:`c_0` (from which each of the data sets
was drawn) as coming from a normal distribution with some mean and some
standard deviation.
To do - write in details of the likelihood and its derivation
For the samples in the figure above, the posterior on these hyperparamters
is given by
.. image:: images/hyper_parameter_corner.png
docs/images/hyper_parameter_combined_posteriors.png

11.6 KiB

docs/images/hyper_parameter_corner.png

121 KiB

docs/images/hyper_parameter_data.png

30.1 KiB

......@@ -15,5 +15,6 @@ Welcome to tupak's documentation!
samplers
tupak-output
writing-documentation
hyperparameters
......@@ -36,7 +36,7 @@ Then, the likelihood for all :math:`N` data points is
In practise, we implement the log-likelihood to avoid numerical overflow
errors. To code this up in :code:`tupak`, we would write a class like this::
class GaussianLikelihood(tupak.Likelihood):
class SimpleGaussianLikelihood(tupak.Likelihood):
def __init__(self, data):
"""
A very simple Gaussian likelihood
......@@ -105,7 +105,7 @@ the likelihood for each data point.
In :code:`tupak`, we can code this up as a likelihood in the following way::
class GaussianLikelihood(tupak.Likelihood):
class GaussianLikelihoodKnownNoise(tupak.Likelihood):
def __init__(self, x, y, sigma, function):
"""
A general Gaussian likelihood - the parameters are inferred from the
......@@ -139,10 +139,6 @@ In :code:`tupak`, we can code this up as a likelihood in the following way::
return -0.5 * (np.sum((res / self.sigma)**2)
+ self.N*np.log(2*np.pi*self.sigma**2))
def noise_log_likelihood(self):
return -0.5 * (np.sum((self.y / self.sigma)**2)
+ self.N*np.log(2*np.pi*self.sigma**2))
This likelihood can be given any python function, the data (in the form of
:code:`x` and :code:`y`) and the standard deviation of the noise. The
......@@ -164,6 +160,20 @@ be the dependent variable.
You can see an example of this likelihood in the `linear regression example
<https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/linear_regression.py>`_.
General likelihood for fitting a function :math:`y(x)` to some data with unknown noise
--------------------------------------------------------------------------------------
In the last example, we considered only cases with known noise (e.g., a
prespecified standard deviation. We now present a general function which can
handle unknown noise (in which case you need to specify a prior for
:math:`\sigma`, or known noise (in which case you pass the known noise in when
instatiating the likelihood
.. literalinclude:: ../examples/other_examples/linear_regression_unknown_noise.py
:lines: 52-94
An example using this likelihood can be found `on this page <https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/linear_regression_unknown_noise.py>`_.
Likelihood for transient gravitational waves
--------------------------------------------
......@@ -173,11 +183,11 @@ routine gravitational wave data analysis of transient events, we have in built
likelihoods. The default likelihood we use in the examples is
`GravitationalWaveTransient`:
.. autoclass:: tupak.likelihood.GravitationalWaveTransient
.. autoclass:: tupak.GravitationalWaveTransient
We also provide a simpler likelihood
.. autoclass:: tupak.likelihood.BasicGravitationalWaveTransient
.. autoclass:: tupak.gw.likelihood.BasicGravitationalWaveTransient
Empty likelihood for subclassing
--------------------------------
......@@ -185,4 +195,4 @@ Empty likelihood for subclassing
We provide an empty parent class which can be subclassed for alternative use
cases
.. autoclass:: tupak.likelihood.Likelihood
.. autoclass:: tupak.Likelihood
......@@ -4,7 +4,7 @@
Priors
======
The priors object passed to :ref:`run_sampler <run-sampler>` is just a regular
The priors object passed to :ref:`run_sampler <run_sampler>` is just a regular
`python dictionary <https://docs.python.org/2/tutorial/datastructures.html#dictionaries>`_.
The keys of the priors objects should reference the model parameters (in
......@@ -28,7 +28,7 @@ when generating plots.
We have provided a number of standard priors. Here is a complete list
.. automodule:: tupak.prior
.. automodule:: tupak.core.prior
:members:
......@@ -8,4 +8,4 @@ Given a :ref:`likelihood` and :ref:`priors`, we can run parameter estimation usi
`run_sampler` function. This can be accessed via :code:`tupak.run_sampler` or
:code:`tupak.sampler.run_sampler`. Here is the detailed API:
.. autofunction:: tupak.sampler.run_sampler
.. autofunction:: tupak.run_sampler
......@@ -6,17 +6,20 @@ This example estimates the masses using a uniform prior in both component masses
comoving volume prior on luminosity distance between luminosity distances of 100Mpc and 5Gpc, the cosmology is WMAP7.
"""
from __future__ import division, print_function
import tupak
import numpy as np
import tupak
# 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 = 'basic_tutorial'
tupak.utils.setup_logger(outdir=outdir, label=label)
tupak.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
......@@ -31,13 +34,13 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
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(
IFOs = [tupak.gw.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']]
......@@ -53,13 +56,15 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, 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['geocent_time'] = tupak.prior.Uniform(injection_parameters['geocent_time'] - 1,
injection_parameters['geocent_time'] + 1,
priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
priors['geocent_time'] = tupak.core.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.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=False, phase_marginalization=False, distance_marginalization=False, prior=priors)
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
time_marginalization=False, phase_marginalization=False,
distance_marginalization=False, prior=priors)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
......
......@@ -10,13 +10,15 @@ 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
import tupak.gw.prior
time_duration = 4.
sampling_frequency = 2048.
# Specify the output directory and the name of the simulation.
outdir = 'outdir'
label = 'basic_tutorial_dist_time_phase_marg'
tupak.utils.setup_logger(outdir=outdir, label=label)
tupak.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
......@@ -31,13 +33,13 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
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(
IFOs = [tupak.gw.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']]
......@@ -52,10 +54,12 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, 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['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=True, phase_marginalization=True, distance_marginalization=True, prior=priors)
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
time_marginalization=True, phase_marginalization=True,
distance_marginalization=True, prior=priors)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
......
......@@ -10,13 +10,15 @@ 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
import tupak.gw.prior
time_duration = 4.
sampling_frequency = 2048.
# Specify the output directory and the name of the simulation.
outdir = 'outdir'
label = 'basic_tutorial_time_phase_marg'
tupak.utils.setup_logger(outdir=outdir, label=label)
tupak.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
......@@ -31,13 +33,13 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
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(
IFOs = [tupak.gw.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']]
......@@ -52,10 +54,12 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, 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['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=True, phase_marginalization=True, distance_marginalization=False, prior=priors)
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
time_marginalization=True, phase_marginalization=True,
distance_marginalization=False, prior=priors)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
......
......@@ -9,7 +9,8 @@ from __future__ import division, print_function
import tupak
import numpy as np
tupak.utils.setup_logger(log_level="info")
tupak.core.utils.setup_logger(log_level="info")
time_duration = 4.
sampling_frequency = 2048.
......@@ -22,16 +23,16 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.waveform_generator.WaveformGenerator(
waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(
sampling_frequency=sampling_frequency, time_duration=time_duration,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
parameter_conversion=tupak.conversion.convert_to_lal_binary_black_hole_parameters,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
parameter_conversion=tupak.gw.conversion.convert_to_lal_binary_black_hole_parameters,
non_standard_sampling_parameter_keys=['chirp_mass', 'mass_ratio', 'cos_iota'],
parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers.
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
IFOs = [tupak.gw.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']]
......@@ -40,14 +41,14 @@ priors = dict()
# These parameters will not be sampled
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'psi', 'ra', 'dec', 'geocent_time']:
priors[key] = injection_parameters[key]
priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance')
priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
# Initialise GravitationalWaveTransient
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
# Run sampler
result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, label='DifferentParameters',
outdir=outdir, conversion_function=tupak.conversion.generate_all_bbh_parameters)
result = tupak.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, label='DifferentParameters',
outdir=outdir, conversion_function=tupak.gw.conversion.generate_all_bbh_parameters)
result.plot_corner()
print(result)
......@@ -7,7 +7,7 @@ import tupak
import numpy as np
# First set up logging and some output directories and labels
tupak.utils.setup_logger()
tupak.core.utils.setup_logger()
outdir = 'outdir'
label = 'create_your_own_source_model'
sampling_frequency = 4096
......@@ -26,14 +26,14 @@ def sine_gaussian(f, A, f0, tau, phi0, geocent_time, ra, dec, psi):
# We now define some parameters that we will inject and then a waveform generator
injection_parameters = dict(A=1e-23, f0=100, tau=1, phi0=0, geocent_time=0,
ra=0, dec=0, psi=0)
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=sine_gaussian,
parameters=injection_parameters)
waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=sine_gaussian,
parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers.
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
IFOs = [tupak.gw.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)
......@@ -42,12 +42,12 @@ IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
# Here we define the priors for the search. We use the injection parameters
# except for the amplitude, f0, and geocent_time
prior = injection_parameters.copy()
prior['A'] = tupak.prior.PowerLaw(alpha=-1, minimum=1e-25, maximum=1e-21, name='A')
prior['f0'] = tupak.prior.Uniform(90, 110, 'f')
prior['A'] = tupak.core.prior.PowerLaw(alpha=-1, minimum=1e-25, maximum=1e-21, name='A')
prior['f0'] = tupak.core.prior.Uniform(90, 110, 'f')
likelihood = tupak.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
result = tupak.sampler.run_sampler(
result = tupak.core.sampler.run_sampler(
likelihood, prior, sampler='dynesty', outdir=outdir, label=label,
resume=False, sample='unif', injection_parameters=injection_parameters)
result.plot_corner()
......
......@@ -8,7 +8,7 @@ and then recovered.
import tupak
import numpy as np
tupak.utils.setup_logger()
tupak.core.utils.setup_logger()
# define the time-domain model
......@@ -33,9 +33,9 @@ outdir='outdir'
label='time_domain_source_model'
# call the waveform_generator to create our waveform model.
waveform = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration, sampling_frequency=sampling_frequency,
time_domain_source_model=time_domain_damped_sinusoid,
parameters=injection_parameters)
waveform = tupak.gw.waveform_generator.WaveformGenerator(time_duration=time_duration, sampling_frequency=sampling_frequency,
time_domain_source_model=time_domain_damped_sinusoid,
parameters=injection_parameters)
hf_signal = waveform.frequency_domain_strain()
#note we could plot the time domain signal with the following code
......@@ -48,7 +48,7 @@ hf_signal = waveform.frequency_domain_strain()
# inject the signal into three interferometers
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
IFOs = [tupak.gw.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)
......@@ -57,19 +57,19 @@ IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
# create the priors
prior = injection_parameters.copy()
prior['amplitude'] = tupak.prior.Uniform(1e-23, 1e-21, r'$h_0$')
prior['damping_time'] = tupak.prior.Uniform(0, 1, r'damping time')
prior['frequency'] = tupak.prior.Uniform(0, 200, r'frequency')
prior['phase'] = tupak.prior.Uniform(-np.pi/2, np.pi/2, r'$\phi$')
prior['amplitude'] = tupak.core.prior.Uniform(1e-23, 1e-21, r'$h_0$')
prior['damping_time'] = tupak.core.prior.Uniform(0, 1, r'damping time')
prior['frequency'] = tupak.core.prior.Uniform(0, 200, r'frequency')
prior['phase'] = tupak.core.prior.Uniform(-np.pi / 2, np.pi / 2, r'$\phi$')
# define likelihood
likelihood = tupak.likelihood.GravitationalWaveTransient(IFOs, waveform)
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(IFOs, waveform)
# launch sampler
result = tupak.sampler.run_sampler(likelihood, prior, sampler='dynesty', npoints=1000,
injection_parameters=injection_parameters,
outdir=outdir, label=label)
result = tupak.core.sampler.run_sampler(likelihood, prior, sampler='dynesty', npoints=1000,
injection_parameters=injection_parameters,
outdir=outdir, label=label)
result.plot_corner()
print(result)
......@@ -6,7 +6,9 @@ from __future__ import division, print_function
import tupak
import numpy as np
tupak.utils.setup_logger()
import tupak.gw.prior
tupak.core.utils.setup_logger()
time_duration = 4.
sampling_frequency = 2048.
......@@ -19,14 +21,14 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
parameters=injection_parameters)
waveform_generator = tupak.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers.
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
IFOs = [tupak.gw.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']]
......@@ -36,32 +38,32 @@ priors = dict()
for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra', 'dec', 'geocent_time', 'psi']:
priors[key] = injection_parameters[key]
# We can assign a default prior distribution, note this only works for certain parameters.
priors['mass_1'] = tupak.prior.create_default_prior(name='mass_1')
priors['mass_1'] = tupak.gw.prior.create_default_prior(name='mass_1')
# We can make uniform distributions.
priors['mass_2'] = tupak.prior.Uniform(name='mass_2', minimum=20, maximum=40)
priors['mass_2'] = tupak.core.prior.Uniform(name='mass_2', minimum=20, maximum=40)
# We can load a prior distribution from a file, e.g., a uniform in comoving volume distribution.
# If no path is given it will look in it's directory of known distributions.
# Note: that this file is used for the default prior distribution for distance.
# Also note: this special case is coded in as tupak.prior.UniformComovingVolume.
priors['luminosity_distance'] = tupak.prior.FromFile('comoving.txt', name='luminosity_distance',
minimum=1e3, maximum=5e3)
priors['luminosity_distance'] = tupak.core.prior.FromFile('comoving.txt', name='luminosity_distance',
minimum=1e3, maximum=5e3)
# We can make a power-law distribution, p(x) ~ x^{alpha}
# Note: alpha=0 is a uniform distribution, alpha=-1 is uniform-in-log
priors['a_1'] = tupak.prior.PowerLaw(name='a_1', alpha=-1, minimum=1e-2, maximum=1)
priors['a_1'] = tupak.core.prior.PowerLaw(name='a_1', alpha=-1, minimum=1e-2, maximum=1)
# We can define a prior from an array as follows.
# Note: this doesn't have to be properly normalised.
a_2 = np.linspace(0, 1, 1001)
p_a_2 = a_2 ** 4
priors['a_2'] = tupak.prior.Interped(name='a_2', xx=a_2, yy=p_a_2, minimum=0, maximum=0.5)
priors['a_2'] = tupak.core.prior.Interped(name='a_2', xx=a_2, yy=p_a_2, minimum=0, maximum=0.5)
# Additionally, we have Gaussian, TruncatedGaussian, Sine and Cosine.
# Finally, if you don't specify any necessary parameters it will be filled in from the default when the sampler starts.
# Enjoy.
# Initialise GravitationalWaveTransient
likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
# Run sampler
result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, outdir=outdir, label='specify_prior')
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, outdir=outdir, label='specify_prior')
result.plot_corner()
print(result)
......@@ -7,7 +7,7 @@ from __future__ import division, print_function
import tupak
import numpy as np
tupak.utils.setup_logger()
tupak.core.utils.setup_logger()
time_duration = 4.
sampling_frequency = 2048.
......@@ -20,13 +20,13 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = tupak.waveform_generator.WaveformGenerator(
waveform_generator = tupak.WaveformGenerator(
time_duration=time_duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole, parameters=injection_parameters)
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers.
IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
IFOs = [tupak.gw.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']]
......@@ -38,12 +38,12 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota
# Initialise GravitationalWaveTransient
# Note that we now need to pass the: priors and flags for each thing that's being marginalised.
likelihood = tupak.likelihood.GravitationalWaveTransient(
likelihood = tupak.GravitationalWaveTransient(
interferometers=IFOs, waveform_generator=waveform_generator, prior=priors,
distance_marginalization=True, phase_marginalization=True)
# Run sampler
result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, outdir=outdir, label='BasicTutorial')
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
injection_parameters=injection_parameters, outdir=outdir, label='BasicTutorial')
result.plot_corner()
print(result)
......@@ -7,15 +7,16 @@ This example estimates all 15 parameters of the binary black hole system using
commonly used prior distributions. This will take a few hours to run.
"""
from __future__ import division, print_function
import tupak
# Define some convienence labels and the trigger time of the event
import tupak.gw.likelihood
outdir = 'outdir'
label = 'GW150914'
time_of_event = tupak.utils.get_event_time(label)
time_of_event = tupak.core.utils.get_event_time(label)
# This sets up logging output to understand what tupak is doing
tupak.utils.setup_logger(outdir=outdir, label=label)
tupak.core.utils.setup_logger(outdir=outdir, label=label)
# Here we import the detector data. This step downloads data from the
# LIGO/Virgo open data archives. The data is saved to an `Interferometer`
......@@ -24,7 +25,7 @@ tupak.utils.setup_logger(outdir=outdir, label=label)
# makes sense, for each detector a plot is created in the `outdir` called
# H1_frequency_domain_data.png and LI_frequency_domain_data.png. The two
# objects are then placed into a list.
interferometers = tupak.detector.get_event_data(label)
interferometers = tupak.gw.detector.get_event_data(label)
# We now define the prior. You'll notice we only do this for the two masses,
# the merger time, and the distance; in each case choosing a prior which
......@@ -33,12 +34,12 @@ interferometers = tupak.detector.get_event_data(label)
# using the syntax below, or choose a fixed value by just providing a float
# value as the prior.
prior = dict()
prior['mass_1'] = tupak.prior.Uniform(30, 50, 'mass_1')
prior['mass_2'] = tupak.prior.Uniform(20, 40, 'mass_2')
prior['geocent_time'] = tupak.prior.Uniform(
prior['mass_1'] = tupak.core.prior.Uniform(30, 50, 'mass_1')
prior['mass_2'] = tupak.core.prior.Uniform(20, 40, 'mass_2')
prior['geocent_time'] = tupak.core.prior.Uniform(
time_of_event - 0.1, time_of_event + 0.1, name='geocent_time')
#prior['geocent_time'] = time_of_event
prior['luminosity_distance'] = tupak.prior.PowerLaw(
prior['luminosity_distance'] = tupak.core.prior.PowerLaw(
alpha=2, minimum=100, maximum=1000)
# In this step we define a `waveform_generator`. This is out object which
......@@ -47,13 +48,13 @@ prior['luminosity_distance'] = tupak.prior.PowerLaw(
# the waveform approximant and reference frequency.
waveform_generator = tupak.WaveformGenerator(time_duration=interferometers[0].duration,
sampling_frequency=interferometers[0].sampling_frequency,
frequency_domain_source_model=tupak.source.lal_binary_black_hole,
frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
parameters={'waveform_approximant': 'IMRPhenomPv2',
'reference_frequency': 50})
# In this step, we define the likelihood. Here we use the standard likelihood
# function, passing it the data and the waveform generator.
likelihood = tupak.GravitationalWaveTransient(interferometers, waveform_generator)
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
# Finally, we run the sampler. This function takes the likelihood and prio
# along with some options for how to do the sampling and how to save the data
......
......@@ -6,9 +6,9 @@ stimation on GW150914 using open data.
"""
import tupak
t0 = tupak.utils.get_event_time("GW150914")
prior = dict(geocent_time=tupak.prior.Uniform(t0-0.1, t0+0.1, name='geocent_time'))
interferometers = tupak.detector.get_event_data("GW150914")
likelihood = tupak.likelihood.get_binary_black_hole_likelihood(interferometers)
t0 = tupak.gw.utils.get_event_time("GW150914")
prior = dict(geocent_time=tupak.core.prior.Uniform(t0 - 0.1, t0 + 0.1, name='geocent_time'))
interferometers = tupak.gw.detector.get_event_data("GW150914")
likelihood = tupak.gw.likelihood.get_binary_black_hole_likelihood(interferometers)
result = tupak.run_sampler(likelihood, prior, label='GW150914')
result.plot_corner()
......@@ -8,7 +8,7 @@ import tupak
import numpy as np
# A few simple setup steps
tupak.utils.setup_logger()
tupak.core.utils.setup_logger()
label = 'gaussian_example'
outdir = 'outdir'
......@@ -22,7 +22,7 @@ outdir = 'outdir'
data = np.random.normal(3, 4, 100)
class GaussianLikelihood(tupak.Likelihood):
class SimpleGaussianLikelihood(tupak.Likelihood):
def __init__(self, data):
"""
A very simple Gaussian likelihood
......@@ -44,9 +44,9 @@ class GaussianLikelihood(tupak.Likelihood):
+ self.N*np.log(2*np.pi*sigma**2))
likelihood = GaussianLikelihood(data)
priors = dict(mu=tupak.prior.Uniform(0, 5, 'mu'),
sigma=tupak.prior.Uniform(0, 10, 'sigma'))
likelihood = SimpleGaussianLikelihood(data)
priors = dict(mu=tupak.core.prior.Uniform(0, 5, 'mu'),
sigma=tupak.core.prior.Uniform(0, 10, 'sigma'))
# And run sampler
result = tupak.run_sampler(
......
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