Skip to content
Snippets Groups Projects
Commit b5d0877f authored by Colm Talbot's avatar Colm Talbot
Browse files

Merge branch 'master' into check_signal_duration

parents 746a2e72 39111e38
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 981 additions and 177 deletions
.gitignore deleted 100644 → 0
+ 0
− 123
View file @ 746a2e72
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
.version
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
.pytest_cache/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
db.sqlite3
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
# Idea
/.idea
*.idea
*.idea*
#.txt files
#*.txt
# Images
*.jpg
*.png
# Output directories
**/outdir
\ No newline at end of file
.gitlab-ci.yml
+ 22
− 8
View file @ b5d0877f
......@@ -21,25 +21,39 @@ exitcode-jessie:
- pip install -r requirements.txt
- pip install coverage
- pip install coverage-badge
- python setup.py install
- coverage erase
- coverage run -a test/prior_tests.py
- coverage run -a test/tests.py
- coverage run -a test/waveform_generator_tests.py
- coverage run -a test/noise_realisation_tests.py
- coverage html --include=tupak/*
- coverage-badge -o coverage.svg
- 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 --omit=* -a test/example_tests.py
- coverage run --omit=* -a test/noise_realisation_tests.py
- coverage run --omit=* -a test/tests.py
# Make the documentation
- pip install -r docs/requirements.txt
- cd docs
- make clean
- make html
artifacts:
paths:
- htmlcov/
- coverage.svg
- coverage_badge.svg
- docs/_build/html/
pages:
stage: deploy
dependencies:
- exitcode-jessie
script:
- mkdir public/
- mv htmlcov/ public/
- mv coverage.svg public/
- mv /builds/Monash/tupak/coverage_badge.svg public/
- mv docs/_build/html/* public/
artifacts:
paths:
- public
......
README.md
+ 30
− 27
View file @ b5d0877f
[![pipeline status](https://git.ligo.org/Monash/tupak/badges/master/pipeline.svg)](https://git.ligo.org/Monash/tupak/commits/master)
[![coverage report](https://monash.docs.ligo.org/tupak/coverage.svg)](
https://monash.docs.ligo.org/tupak/)
[![coverage report](https://monash.docs.ligo.org/tupak/coverage_badge.svg)](
https://monash.docs.ligo.org/tupak/htmlcov/)
# Tupak
Fulfilling all your GW dreams.
Fulfilling all your gravitational wave dreams.
Documentation can be found
[here](https://monash.docs.ligo.org/tupak/index.html). This and the code in
general are a work in progress. We encourage you to submit issues via are
[issue tracker](https://git.ligo.org/Monash/tupak/issues) and contribute to the
development via a merge request (for help in creating a merge request, see
[this page](https://docs.gitlab.com/ee/gitlab-basics/add-merge-request.html) or
contact us directly.
## Example
......@@ -13,6 +21,7 @@ To get started with `tupak`, we have a few examples and tutorials:
1. [Examples of injecting and recovering data](https://git.ligo.org/Monash/tupak/tree/master/examples/injection_examples)
* [A basic tutorial](https://git.ligo.org/Monash/tupak/blob/master/examples/injection_examples/basic_tutorial.py)
* [Create your own source model](https://git.ligo.org/Monash/tupak/blob/master/examples/injection_examples/create_your_own_source_model.py)
* [Create your own time-domain source model](https://git.ligo.org/Monash/tupak/blob/master/examples/injection_examples/create_your_own_time_domain_source_model.py)
* [How to specify the prior](https://git.ligo.org/Monash/tupak/blob/master/examples/injection_examples/how_to_specify_the_prior.py)
* [Using a partially marginalized likelihood](https://git.ligo.org/Monash/tupak/blob/master/examples/injection_examples/marginalized_likelihood.py)
......@@ -26,8 +35,10 @@ To get started with `tupak`, we have a few examples and tutorials:
## Installation
In the following, we assume you have installed
[https](pip://pypa.io.en/stable/installing/pip/) and [git](https://git-scm.com/).
`tupak` is developed to work with both Python 2.7+ and Python 3+. In the
following, we assume you have a working python installation, [python
pip](https://packaging.python.org/tutorials/installing-packages/#use-pip-for-installing),
and [git](https://git-scm.com/).
### Install tupak
Clone the repository, install the requirements, and then install `tupak`.
......@@ -38,12 +49,18 @@ $ pip install -r requirements.txt
$ python setup.py install
```
Once you have run these steps, you have `tupak` installed.
Once you have run these steps, you have `tupak` installed. You can now try to
run the examples. However, for many gravitational-wave related examples you'll
also need to install `lalsuite`.
### Install lalsuite
The simple way: `pip install lalsuite`, or,
from source:
We recommend you install `lalsuite` the simple way:
```bash
$ pip install lalsuite.
```
If this doesn't work, or you prefer to, you can also install from source.
Head to
[https://git.ligo.org/lscsoft/lalsuite](https://git.ligo.org/lscsoft/lalsuite)
to check you have an account and SSH keys set up. Then,
......@@ -60,32 +77,18 @@ $ make; make install
Warning: in the configure line here, we have disabled everything except
lalsimulation. If you need other modules, see `./configure --help`.
### Install pymultinest
If you want to use the `pymultinest` sampler, you first need the
MultiNest library to be installed to work properly. The full instructions can
be found [here](https://johannesbuchner.github.io/PyMultiNest/install.html). We
have also written [a shortened tl;dr here](./TLDR_MULTINEST.md).
## Tests and coverage
To locally test the code
```bash
$ python tests.py
```
To locally generate a coverage report
```bash
$ pip install coverage
$ coverage run tests.py
$ coverage html
```
This will generate a directory `htmlcov`, to see detailed coverage navigate
from your browser to the file `tupak/htmlcov/index.html`.
## Tests and coverage
The coverage report for master can be seen here:
[https://monash.docs.ligo.org/tupak/](https://monash.docs.ligo.org/tupak/).
We have a variety of tests which can be found in the `tests` directory. You
can find a [current coverage report for master
here.](https://monash.docs.ligo.org/tupak/htmlcov/).
docs/Makefile 0 → 100644
+ 20
− 0
View file @ b5d0877f
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = python -msphinx
SPHINXPROJ = tupak
SOURCEDIR = .
BUILDDIR = _build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
docs/basics-of-parameter-estimation.txt 0 → 100644
+ 112
− 0
View file @ b5d0877f
==============================
Basics of parameter estimation
==============================
In this example, we'll go into some of the basics of parameter estimation and
how they are implemented in `tupak`.
Firstly, consider a situation where you have discrete data :math:`\{y_0,
y_1,\ldots, y_n\}` taken at a set of times :math:`\{t_0, t_1, \ldots, y_n\}`.
Further, we know that this data is generated by a process which can be modelled
by a linear function of the form :math:`y(t) = m t + c`. We will refer to
this model as :math:`H`. Given a set of data,
how can we figure out the coefficients :math:`m` and :math:`c`?
This is a well studied problem known as `linear regression
<https://en.wikipedia.org/wiki/Linear_regression>`_ and there already exists
many ways to estimate the coefficients (you may already be familiar with a
least squares estimator for example).
Here, we will describe a Bayesian approach using `nested sampling
<https://en.wikipedia.org/wiki/Nested_sampling_algorithm>`_ which might feel
like overkill for this problem. However, it is a useful way to introduce some
of the basic features of `tupak` before seeing them in more complicated
settings.
The maths
---------
Given the data, the posterior distribution for the model parameters is given
by
.. math::
P(m, c| \{y_i, t_i\}, H) \propto P(\{y_i, t_i\}| m, c, H) \times P(m, c| H)
where the first term on the right-hand-side is the *likelihood* while the
second is the *prior*. In the model :math:`H`, the likelihood of the data point
:math:`y_i, t_i`, given a value for the coefficients we will define to be
Gaussian distributed as such:
.. math::
P(y_i, t_i| m, c, H) = \frac{1}{\sqrt{2\pi\sigma^2}}
\mathrm{exp}\left(\frac{-(y_i - (t_i x + c))^2}{2\sigma^2}\right)
Next, we assume that all data points are independent. As such,
.. math::
P(\{y_i, t_i\}| m, c, H) = \prod_{i=1}^n P(y_i, t_i| m, c, H)
When solving problems on a computer, it is often convienient to work with
the log-likelihood. Indeed, a `tupak` *likelihood* must have a `log_likelihood()`
method. For the normal distribution, the log-likelihood for *n* data points is
.. math::
\log(P(\{y_i, t_i\}| m, c, H)) = -\frac{1}{2}\left[
\sum_{i=1}^n \left(\frac{(y_i - (t_i x + c))}{\sigma}\right)^2
+ n\log\left(2\pi\sigma^2\right)\right]
Finally, we need to specify a *prior*. In this case we will use uncorrelated
uniform priors
.. math::
P(m, c| H) = P(m| H) \times P(c| H) = \textrm{Unif}(0, 5) \times \textrm{Unif}(-2, 2)
the choice of prior in general should be guided by physical knowledge about the
system and not the data in question.
The key point to take away here is that the **likelihood** and **prior** are
the inputs to figuring out the **posterior**. There are many ways to go about
this, we will now show how to do so in `tupak`. In this case, we explicitly
show how to write the `GaussianLikelihood` so that one can see how the
maths above gets implemented. For the prior, this is done implicitly by the
naming of the priors.
The code
--------
In the following example, also available under
`examples/other_examples/linear_regression.py
<https://git.ligo.org/Monash/tupak/tree/master/examples/other_examples/linear_regression.py>`_
we will step through the process of generating some simulated data, writing
a likelihood and prior, and running nested sampling using `tupak`.
.. literalinclude:: /../examples/other_examples/linear_regression.py
:language: python
:linenos:
Running the script above will make a few images. Firstly, the plot of the data:
.. image:: images/linear-regression_data.png
The dashed red line here shows the simulated signal.
Secondly, because we used the `plot=True` argument in `run_sampler` we generate
a corner plot
.. image:: images/linear-regression_corner.png
The solid lines indicate the simulated values which are recovered quite
easily. Note, you can also make a corner plot with `result.plot_corner()`.
Final thoughts
--------------
While this example is somewhat trivial, hopefully you can see how this script
can be easily modified to perform parameter estimation for almost any
time-domain data where you can model the background noise as Gaussian and write
the signal model as a python function (i.e., replacing `model`).
docs/compact-binary-coalescence-parameter-estimation.txt 0 → 100644
+ 27
− 0
View file @ b5d0877f
===============================================
Compact binary coalescence parameter estimation
===============================================
In this example, we demonstrate how to generate simulated data for a binary
black hold coalescence observed by the two LGIO interferometers at Hanford,
Livingston and the Virgo detector.
.. literalinclude:: /../examples/injection_examples/basic_tutorial.py
:language: python
:linenos:
Running this script will generate data then perform parameter estimation for
the luminosity distance, masses and inclination angle :math:`\iota`. In doing
all of this, it prints information about the matched-filter SNRs in each
detector (saved to the log-file). Moreover, it generates a plot for each
detector showing the data, amplitude spectral density (ASD) and the signal;
here is an example for the Hanford detector:
.. image:: images/H1_frequency_domain_data.png
Finally, after running the parameter estimation. It generates a corner plot:
.. image:: images/basic_tutorial_corner.png
The solid lines indicate the injection parameters.
docs/conf.py 0 → 100644
+ 173
− 0
View file @ b5d0877f
# -*- coding: utf-8 -*-
#
# tupak documentation build configuration file, created by
# sphinx-quickstart on Fri May 25 12:08:01 2018.
#
# This file is execfile()d with the current directory set to its
# containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
import os
import sys
sys.path.insert(0, os.path.abspath('../tupak/'))
# -- General configuration ------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#
# needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.mathjax', 'numpydoc',
'nbsphinx', 'sphinx.ext.autosummary']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix(es) of source filenames.
# You can specify multiple suffix as a list of string:
#
source_suffix = ['.rst', '.md']
source_suffix = '.txt'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = u'tupak'
copyright = u'2018, Paul Lasky'
author = u'Paul Lasky'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = u'0.2'
# The full version, including alpha/beta/rc tags.
release = u'0.2'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#
# This is also used if you do content translation via gettext catalogs.
# Usually you set "language" from the command line for these cases.
language = None
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This patterns also effect to html_static_path and html_extra_path
exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', 'requirements.txt']
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# If true, `todo` and `todoList` produce output, else they produce nothing.
todo_include_todos = False
# -- Options for HTML output ----------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
html_theme = 'sphinx_rtd_theme'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#
# html_theme_options = {}
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# Custom sidebar templates, must be a dictionary that maps document names
# to template names.
#
# This is required for the alabaster theme
# refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars
html_sidebars = {
'**': [
'about.html',
'navigation.html',
'relations.html', # needs 'show_related': True theme option to display
'searchbox.html',
'donate.html',
]
}
# -- Options for HTMLHelp output ------------------------------------------
# Output file base name for HTML help builder.
htmlhelp_basename = 'tupakdoc'
# -- Options for LaTeX output ---------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#
# 'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#
# 'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#
# 'preamble': '',
# Latex figure (float) alignment
#
# 'figure_align': 'htbp',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'tupak.tex', u'tupak Documentation',
u'Paul Lasky', 'manual'),
]
# -- Options for manual page output ---------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
(master_doc, 'tupak', u'tupak Documentation',
[author], 1)
]
# -- Options for Texinfo output -------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
(master_doc, 'tupak', u'tupak Documentation',
author, 'tupak', 'One line description of project.',
'Miscellaneous'),
]
numpydoc_show_class_members = False
docs/images/H1_frequency_domain_data.png 0 → 100644
+ 0
− 0
View file @ b5d0877f
docs/images/H1_frequency_domain_data.png

35.9 KiB

docs/images/basic_tutorial_corner.png 0 → 100644
+ 0
− 0
View file @ b5d0877f
docs/images/basic_tutorial_corner.png

395 KiB

docs/images/linear-regression_corner.png 0 → 100644
+ 0
− 0
View file @ b5d0877f
docs/images/linear-regression_corner.png

110 KiB

docs/images/linear-regression_data.png 0 → 100644
+ 0
− 0
View file @ b5d0877f
docs/images/linear-regression_data.png

22.2 KiB

docs/index.txt 0 → 100644
+ 19
− 0
View file @ b5d0877f
Welcome to tupak's documentation!
=================================
.. automodule:: tupak
:members:
.. toctree::
:maxdepth: 3
:caption: Contents:
basics-of-parameter-estimation
compact-binary-coalescence-parameter-estimation
prior
likelihood
samplers
tupak-output
writing-documentation
docs/likelihood.txt 0 → 100644
+ 188
− 0
View file @ b5d0877f
.. _likelihood:
==========
Likelihood
==========
`tupak` likelihood objects are used in calculating the likelihood of the data
for some specific set of parameters. In mathematical notation, the likelihood
can be generically written as :math:`\mathcal{L}(d| \theta)`. How this is
coded up will depend on the problem, but `tupak` expects all likelihood
objects to have a `parameters` attribute (a dictionary of key-value pairs) and
a `log_likelihood()` method. In thie page, we'll discuss how to write your own
Likelihood, and the standard likelihoods in :code:`tupak`.
The simplest likelihood
-----------------------
To start with let's consider perhaps the simplest likelihood we could write
down, namely a Gaussian likelihood for a set of data :math:`\vec{x}=[x_1, x_2,
\ldots, x_N]`. The likelihood for a single data point, given the mean
:math:`\mu` and standard-deviation :math:`\sigma` is given by
.. math::
\mathcal{L}(x_i| \mu, \sigma) =
\frac{1}{\sqrt{2\pi\sigma^2}}\mathrm{exp}\left(
\frac{-(x_i - \mu)^2}{2\sigma^2}\right)
Then, the likelihood for all :math:`N` data points is
.. math::
\mathcal{L}(\vec{x}| \mu, \sigma) = \prod_{i=1}^N
\mathcal{L}(x_i| \mu, \sigma)
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):
def __init__(self, data):
"""
A very simple Gaussian likelihood
Parameters
----------
data: array_like
The data to analyse
"""
self.data = data
self.N = len(data)
self.parameters = {'mu': None, 'sigma': None}
def log_likelihood(self):
mu = self.parameters['mu']
sigma = self.parameters['sigma']
res = self.data - mu
return -0.5 * (np.sum((res / sigma)**2)
+ self.N*np.log(2*np.pi*sigma**2))
This demonstrates the two required features of a :code:`tupak`
:code:`Likelihood` object:
#. It has a `parameters` attribute (a dictionary with
keys for all the parameters, in this case, initialised to :code:`None`)
#. It has a :code:`log_likelihood` method which, when called returns the log
likelihood for all the data.
You can find an example that uses this likelihood `here <https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/gaussian_example.py>`_.
.. tip::
Note that the example above subclasses the :code:`tupak.Likelihood` base
class, this simply provides a few in built functions. We recommend you also
do this when writing your own likelihood.
General likelihood for fitting a function :math:`y(x)` to some data with known noise
------------------------------------------------------------------------------------
The previous example was rather simplistic, Let's now consider that we have some
dependent data :math:`\vec{y}=y_1, y_2, \ldots y_N$` measured at
:math:`\vec{x}=x_1, x_2, \ldots, x_N`. We believe that the data is generated
by additive Gaussian noise with a known variance :math:`sigma^2` and a function
:math:`y(x; \theta)` where :math:`\theta` are some unknown parameters; that is
.. math::
y_i = y(x_i; \theta) + n_i
where :math:`n_i` is drawn from a normal distribution with zero mean and
standard deviation :math:`sigma`. As such, :math:`y_i - y(x_i; \theta)`
itself will have a likelihood
.. math::
\mathcal{L}(y_i; x_i, \theta) =
\frac{1}{\sqrt{2\pi\sigma^{2}}}
\mathrm{exp}\left(\frac{-(y_i - y(x_i; \theta))^2}{2\sigma^2}\right)
As with the previous case, the likelihood for all the data is the produce over
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):
def __init__(self, x, y, sigma, function):
"""
A general Gaussian likelihood - the parameters are inferred from the
arguments of function
Parameters
----------
x, y: array_like
The data to analyse
sigma: float
The standard deviation of the noise
function:
The python function to fit to the data. Note, this must take the
dependent variable as its first argument. The other arguments are
will require a prior and will be sampled over (unless a fixed
value is given).
"""
self.x = x
self.y = y
self.sigma = sigma
self.N = len(x)
self.function = function
# These lines of code infer the parameters from the provided function
parameters = inspect.getargspec(function).args
parameters.pop(0)
self.parameters = dict.fromkeys(parameters)
def log_likelihood(self):
res = self.y - self.function(self.x, **self.parameters)
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
parameters are inferred from the arguments to the :code:`function` argument,
for example if, when instatiating the likelihood you passed in a the following
function::
def f(x, a, b):
return x**2 + b
Then you would also need to provide priors for :code:`a` and :code:`b`. For
this likelihood, the first argument to :code:`function` is always assumed to
be the dependent variable.
.. note::
Here we have explicitly defined the :code:`noise_log_likelihood` method
as the case when there is no signal (i.e., :math:`y(x; \theta)=0`).
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>`_.
Likelihood for transient gravitational waves
--------------------------------------------
In the examples above, we show how to write your own likelihood. However, for
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
We also provide a simpler likelihood
.. autoclass:: tupak.likelihood.BasicGravitationalWaveTransient
Empty likelihood for subclassing
--------------------------------
We provide an empty parent class which can be subclassed for alternative use
cases
.. autoclass:: tupak.likelihood.Likelihood
docs/prior.txt 0 → 100644
+ 34
− 0
View file @ b5d0877f
.. _priors:
======
Priors
======
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
particular, the :code:`parameters` attribute of the :ref:`likelihood`. Each key
can be either
- fixed number, in which case the value is held fixed at this value. In effect,
this is a Delta-function prior,
- or a :code:`tupak.prior.Prior` instance.
If the later, it will be sampled during the parameter estimation. Here is a
simple example that sets a uniform prior for :code:`a`, and a fixed value for
:code:`b`::
priors = {}
priors['a'] = tupak.prior.Uniform(minimum=0, maximum=10, name='a', latex_label='a')
priors['b'] = 5
Notice, that the :code:`latex_label` is optional, but if given will be used
when generating plots.
We have provided a number of standard priors. Here is a complete list
.. automodule:: tupak.prior
:members:
docs/requirements.txt 0 → 100644
+ 6
− 0
View file @ b5d0877f
html5lib==1.0b8 # Hack to ensure documentation generated
sphinx
numpydoc
nbsphinx
autodoc
sphinx_rtd_theme
docs/samplers.txt 0 → 100644
+ 11
− 0
View file @ b5d0877f
.. _run_sampler:
========
Sampling
========
Given a :ref:`likelihood` and :ref:`priors`, we can run parameter estimation using the
`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
docs/tupak-output.txt 0 → 100644
+ 137
− 0
View file @ b5d0877f
============
Tupak output
============
In this document, we will describe what :code::code:`tupak` outputs, where it is stored,
and how you can access it.
When you call :code:`run_sampler`, there are two arguments :code:`outdir` and :code:`label` which
are used in generating all the file names for saved data. In the rest of these
documents, we'll assume the defaults where used (which are :code:`outdir` and
:code:`label`).
The raw data
------------
First off, the primary data dump of :code:`tupak` goes into :code:`outdir/label_result.h5`.
This is a binary file, so isn't human readable, but you can use it using a
command-line tool :code:`ddls` [#f1]_ (you have this installed already if you have installed
the :code:`tupak` requirements). To have a look at what is inside run
.. code-block:: console
$ ddls outdir/label_result.h5
/fixed_parameter_keys list
/fixed_parameter_keys/i0 'dec' (3) [ascii]
/fixed_parameter_keys/i1 'psi' (3) [ascii]
/fixed_parameter_keys/i2 'a_2' (3) [ascii]
/fixed_parameter_keys/i3 'a_1' (3) [ascii]
/fixed_parameter_keys/i4 'geocent_time' (12) [ascii]
/fixed_parameter_keys/i5 'phi_jl' (6) [ascii]
/fixed_parameter_keys/i6 'ra' (2) [ascii]
/fixed_parameter_keys/i7 'phase' (5) [ascii]
/fixed_parameter_keys/i8 'phi_12' (6) [ascii]
/fixed_parameter_keys/i9 'tilt_2' (6) [ascii]
/fixed_parameter_keys/i10 'tilt_1' (6) [ascii]
/injection_parameters dict
/injection_parameters/a_1 0.4 [float64]
/injection_parameters/a_2 0.3 [float64]
/injection_parameters/dec -1.2108 [float64]
/injection_parameters/geocent_time 1126259642.413 [float64]
/injection_parameters/iota 0.4 [float64]
/injection_parameters/luminosity_distance 4000.0 [float64]
/injection_parameters/mass_1 36.0 [float64]
/injection_parameters/mass_2 29.0 [float64]
/injection_parameters/phase 1.3 [float64]
/injection_parameters/phi_12 1.7 [float64]
/injection_parameters/phi_jl 0.3 [float64]
/injection_parameters/psi 2.659 [float64]
/injection_parameters/ra 1.375 [float64]
/injection_parameters/reference_frequency 50.0 [float64]
/injection_parameters/tilt_1 0.5 [float64]
/injection_parameters/tilt_2 1.0 [float64]
/injection_parameters/waveform_approximant 'IMRPhenomPv2' (12) [ascii]
/kwargs dict
/kwargs/bound 'multi' (5) [ascii]
/kwargs/dlogz 0.1 [float64]
/kwargs/nlive 1000 [int64]
/kwargs/sample 'rwalk' (5) [ascii]
/kwargs/update_interval 600 [int64]
/kwargs/verbose True [bool]
/kwargs/walks 20 [int64]
/label 'basic_tutorial' (14) [ascii]
/log_bayes_factor 29.570224130853056 [float64]
/logz -12042.875089354413 [float64]
/logzerr 0.040985337772488764 [float64]
/noise_logz -12072.445313485267 [float64]
/outdir 'outdir' (6) [ascii]
/parameter_labels list
/parameter_labels/i0 '$d_L$' (5) [ascii]
/parameter_labels/i1 '$m_2$' (5) [ascii]
/parameter_labels/i2 '$m_1$' (5) [ascii]
/parameter_labels/i3 '$\\iota$' (7) [ascii]
/posterior DataFrame (4, 15073)
/search_parameter_keys list
/search_parameter_keys/i0 'luminosity_distance' (19) [ascii]
/search_parameter_keys/i1 'mass_2' (6) [ascii]
/search_parameter_keys/i2 'mass_1' (6) [ascii]
/search_parameter_keys/i3 'iota' (4) [ascii]
This shows the different data that is stored in the :code:`h5` file. You can think of
the file like a python dictionary - its a bag with lots of different kind of
data which can be accessed via a :code:`key` (a string). We use `deepdish
<http://deepdish.io/>`_ to handle the saving of :code:`h5` files in :code:`tupak`. In python,
can load any :code:`h5` file and access its contents like a dictionary::
>>> import deepdish
>>> output = deepdish.io.load('outdir/label_result.h5')
>>> print(output['logz'])
-146.23
Here we printed the stored data for the :code:`'logz'` key, but you could equally get
the :code:`posterior` or any other attribute. For a full list, print the
`output.keys()`.
Reading in a result file
------------------------
Rather than reading in the raw :code:`h5` file, may find it more convienient to
instead load a :code:`*result.h5` as a :code:`tupak` :code:`result` object (i.e., like the output
of :code:`run_sampler`. To do this::
>>> import tupak
>>> result = tupak.result.read_in_result(outdir=outdir, label=label)
>>> result = tupak.result.read_in_result(filename='outdir/label_result.h5)
Note, these two lines are equivalent, but show two different ways to read in
the data. Using this method, :code:`result` is now a :code:`tupak.result.Result` instance
and has all the methods and attributes. So, for example, you could call
`result.plot_corner()` to generate a corner plot.
Accessing samples directly
--------------------------
To get the samples for a particular parameter, use::
>>> result.posterior.[key]
where :code:`key` is a string of the parameter name you are interested in. The
`posterior` attribute is a :code:`pandas.DataFrame`, so if you want to return just
a :code:`numpy` array, use::
>>> result.posterior.[key].values
Saving to ASCII text
--------------------
You may wish to get hold of the samples in a simple text file, this can be
done via::
result.save_posterior_samples()
which will generate a file :code:`outdir/label_posterior.txt`.
.. rubric:: Footnotes
.. [#f1] :code:`ddls` is an acronym for deep-dish ls
docs/writing-documentation.txt 0 → 100644
+ 100
− 0
View file @ b5d0877f
=====================
Writing documentation
=====================
This is a short tutorial on how to contribute to the documentation of ``tupak``.
Writing the basics
------------------
First, open your terminal and head to the ``docs`` directory in your clone of
`tupak`. Once here, you'll notice there are a number of `*txt` files. Each one
is a page of the documentation.
Let's say you want to write documentation about the your new feature (or update
the existing documentation). Simply open a new file ``a-new-feature.txt`` in
your text editor. Then add a title and some text::
====================
A new tupak feature!
====================
Here we'll put a description of the new feature
Next, in order to get your new page known by the rest of the documentation,
open ``index.rst`` and, under ``toctree`` add the name of your file (without
the suffix). For the example above::
.. toctree::
:maxdepth: 3
:caption: Contents:
likelihood
samplers
writing-documentation
a-new-feature
Checking the results
--------------------
You can check what this will look like by (whilst in the ``docs`` directory)
running the command::
$ make html
This will create a directory ``./_build/html`` with the documentation. To
see the result, open the file ``./_build/html/index.html`` in your browser.
Pushing your changes to master
------------------------------
To contribute your documentation changes, you should create a branch and add in
all of the new/changed files::
$ git checkout -b adding-my-new-documentation
$ git add index.txt
$ git add a-new-feature.txt
$ git commit -m "Adding my documentation for the feature"
$ git push origin adding-my-new-documentation
Then, on the web interface create an merge request.
Using reStructured text
-----------------------
The help files are written in reStructured text format. To familiarise yourself
with these features visit http://www.sphinx-doc.org/en/master/usage/restructuredtext/basics.html.
A useful feature is the ability to `format code examples <http://www.sphinx-doc.org/en/stable/markup/code.html>`_. This is done through the use of ``::`` and indendentation. For
example, to make this code block::
import tupak
You would write::
to make this code block::
import tupak
reStructured text is very powerful, but can be quite particular. For example,
**all code blocks must be indented by 3 spaces**.
Sphinx autodoc
--------------
Most of the documentation for ``tupak`` should be written in the `docstrings
<https://www.python.org/dev/peps/pep-0257/>`_ of the functions themselves. We
can add these into the online documentation using `autodoc
<http://www.sphinx-doc.org/en/master/ext/autodoc.html>`_. To add the documentation
for a class, add a line such as::
.. autoclass:: tupak.likelihood.GravitationalWaveTransient
:members:
into your documentation file. Similarly, to document a function::
.. autofunction:: tupak.sampler.run_sampler
examples/injection_examples/basic_tutorial.py
+ 36
− 19
View file @ b5d0877f
......@@ -9,45 +9,62 @@ 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 = 'basic_tutorial'
tupak.utils.setup_logger(outdir=outdir, label=label)
np.random.seed(170809)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# We are going to inject a binary black hole waveform. We first establish a dictionary of parameters that
# includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2),
# spins of both black holes (a, tilt, phi), etc.
injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3,
luminosity_distance=4000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
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.source.lal_binary_black_hole,
parameters=injection_parameters)
hf_signal = waveform_generator.frequency_domain_strain()
# Set up interferometers.
# 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']]
sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1']]
# Set up prior
# Set up prior, which is a dictionary
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', 'iota', 'ra', 'dec', 'geocent_time']:
# 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 ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'dec', 'geocent_time', 'phase']:
priors[key] = injection_parameters[key]
# 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,
'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)
# Initialise Likelihood
likelihood = tupak.likelihood.Likelihood(interferometers=IFOs, waveform_generator=waveform_generator)
# 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,
injection_parameters=injection_parameters, outdir=outdir, label=label)
# Run 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()
result.plot_walks()
result.plot_distributions()
print(result)
examples/injection_examples/basic_tutorial_dist_time_phase_marg.py 0 → 100644
+ 66
− 0
View file @ b5d0877f
#!/bin/python
"""
Tutorial to demonstrate running parameter estimation on a reduced parameter space for an injected signal.
This example estimates the masses using a uniform prior in both component masses and distance using a uniform in
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
# 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_dist_time_phase_marg'
tupak.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# We are going to inject a binary black hole waveform. We first establish a dictionary of parameters that
# includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2),
# spins of both black holes (a, tilt, phi), etc.
injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3,
luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
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.WaveformGenerator(time_duration=time_duration,
sampling_frequency=sampling_frequency,
frequency_domain_source_model=tupak.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(
name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, time_duration=time_duration,
sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1']]
# 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 ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'dec', 'geocent_time', 'phase']:
priors[key] = injection_parameters[key]
# 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')
# 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)
# 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,
injection_parameters=injection_parameters, outdir=outdir, label=label)
# make some plots of the outputs
result.plot_corner()
print(result)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment