.gitignore

-# 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

@@ -23,17 +23,20 @@ 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/tests.py
+    - coverage run --include=tupak/* -a test/sampler_tests.py
     - coverage run --include=tupak/* -a test/waveform_generator_tests.py
-    - coverage run --include=tupak/* -a test/noise_realisation_tests.py
-    - coverage run --include=tupak/* -a test/example_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:

README.md

@@ -6,6 +6,14 @@ https://monash.docs.ligo.org/tupak/htmlcov/)
 
 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
 
 To get started with `tupak`, we have a few examples and tutorials:

docs/basics-of-parameter-estimation.txt

@@ -87,7 +87,6 @@ a likelihood and prior, and running nested sampling using `tupak`.
 
 .. literalinclude:: /../examples/other_examples/linear_regression.py
    :language: python
-   :emphasize-lines: 12,15-18
    :linenos:
 
 Running the script above will make a few images. Firstly, the plot of the data:

docs/compact-binary-coalescence-parameter-estimation.txt

+===============================================
+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

@@ -168,5 +168,6 @@ texinfo_documents = [
      'Miscellaneous'),
 ]
 
+numpydoc_show_class_members = False
 
 

docs/index.txt

@@ -9,15 +9,11 @@ Welcome to tupak's documentation!
    :caption: Contents:
 
    basics-of-parameter-estimation
-   tupak-output
+   compact-binary-coalescence-parameter-estimation
+   prior
    likelihood
    samplers
+   tupak-output
    writing-documentation
 
 
-Indices and tables
-==================
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`

docs/likelihood.txt

+.. _likelihood:
+
 ==========
 Likelihood
 ==========
@@ -7,24 +9,180 @@ 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.
+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
 
-The default likelihood we use in the examples is `GravitationalWaveTransient`:
+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
-    :members:
 
 We also provide a simpler likelihood
 
 .. autoclass:: tupak.likelihood.BasicGravitationalWaveTransient
-    :members:
+
+Empty likelihood for subclassing
+--------------------------------
 
 We provide an empty parent class which can be subclassed for alternative use
 cases
 
 .. autoclass:: tupak.likelihood.Likelihood
-    :members:
-
-
-
-

docs/prior.txt

+.. _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/samplers.txt

-==========
-Samplers
-==========
+.. _run_sampler:
 
-Given a `likelihood` and `prior`, we can run parameter estimation using `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

examples/injection_examples/basic_tutorial.py

@@ -19,27 +19,27 @@ label = 'basic_tutorial'
 tupak.utils.setup_logger(outdir=outdir, label=label)
 
 # Set up a random seed for result reproducibility.  This is optional!
-np.random.seed(170801)
+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.  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, which is a dictionary
 priors = dict()
@@ -47,19 +47,23 @@ priors = dict()
 # 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', 'phase', 'psi', 'ra', 'dec', 'geocent_time']:
+
+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)
+likelihood = tupak.likelihood.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.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
-                                   injection_parameters=injection_parameters, outdir=outdir, label=label)
+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()

examples/injection_examples/basic_tutorial_dist_time_phase_marg.py

+#!/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)

examples/injection_examples/basic_tutorial_time_phase_marg.py

+#!/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_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=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,
+                           injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+print(result)

examples/open_data_examples/GW150914.py

@@ -37,6 +37,7 @@ 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(
     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(
     alpha=2, minimum=100, maximum=1000)
 
@@ -44,20 +45,23 @@ prior['luminosity_distance'] = tupak.prior.PowerLaw(
 # creates the frequency-domain strain. In this instance, we are using the
 # `lal_binary_black_hole model` source model. We also pass other parameters:
 # the waveform approximant and reference frequency.
-waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=interferometers[0].duration,
-                                                                sampling_frequency=interferometers[
-                                                                    0].sampling_frequency,
-                                                                frequency_domain_source_model=tupak.source.lal_binary_black_hole,
-                                                                parameters={'waveform_approximant': 'IMRPhenomPv2',
-                                                                            'reference_frequency': 50})
+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,
+                                             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.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
+<<<<<<< HEAD
+likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers, waveform_generator, time_marginalization=False, distance_marginalization=True, prior=prior)
+=======
+likelihood = tupak.GravitationalWaveTransient(interferometers, waveform_generator)
+>>>>>>> e8b704afba65886f49cb448402a95ae90c468a24
 
 # 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
-result = tupak.sampler.run_sampler(likelihood, prior, sampler='dynesty',
-                                   outdir=outdir, label=label)
+result = tupak.run_sampler(likelihood, prior, sampler='dynesty',
+                           outdir=outdir, label=label)
 result.plot_corner()
 print(result)

examples/open_data_examples/GW150914_minimal.py

@@ -10,5 +10,5 @@ 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)
-result = tupak.sampler.run_sampler(likelihood, prior, label='GW150914')
+result = tupak.run_sampler(likelihood, prior, label='GW150914')
 result.plot_corner()

examples/other_examples/gaussian_example.py

+#!/bin/python
+"""
+An example of how to use tupak to perform paramater estimation for
+non-gravitational wave data consisting of a Gaussian with a mean and variance
+"""
+from __future__ import division
+import tupak
+import numpy as np
+
+# A few simple setup steps
+tupak.utils.setup_logger()
+label = 'gaussian_example'
+outdir = 'outdir'
+
+# Here is minimum requirement for a Likelihood class to run with tupak. In this
+# case, we setup a GaussianLikelihood, which needs to have a log_likelihood
+# method. Note, in this case we will NOT make use of the `tupak`
+# waveform_generator to make the signal.
+
+# Making simulated data: in this case, we consider just a Gaussian
+
+data = np.random.normal(3, 4, 100)
+
+
+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))
+
+
+likelihood = GaussianLikelihood(data)
+priors = dict(mu=tupak.prior.Uniform(0, 5, 'mu'),
+              sigma=tupak.prior.Uniform(0, 10, 'sigma'))
+
+# And run sampler
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
+    walks=10, outdir=outdir, label=label)
+result.plot_corner()
+print(result)

examples/other_examples/linear_regression.py

 #!/bin/python
 """
 An example of how to use tupak to perform paramater estimation for
-non-gravitational wave data
+non-gravitational wave data. In this case, fitting a linear function to
+data with background Gaussian noise
+
 """
 from __future__ import division
 import tupak
 import numpy as np
 import matplotlib.pyplot as plt
+import inspect
 
 # A few simple setup steps
 tupak.utils.setup_logger()
-label = 'linear-regression'
+label = 'linear_regression'
 outdir = 'outdir'
 
-# Here is minimum requirement for a Likelihood class to run linear regression
-# with tupak. In this case, we setup a GaussianLikelihood, which needs to have
-# a log_likelihood method. Note, in this case we make use of the `tupak`
-# waveform_generator to make the signal (more on this later) But, one could
-# make this work without the waveform generator.
-
-# Making simulated data
-
 
-# First, we define our signal model, in this case a simple linear function
+# First, we define our "signal model", in this case a simple linear function
 def model(time, m, c):
     return time * m + c
 
@@ -55,9 +50,11 @@ fig.savefig('{}/{}_data.png'.format(outdir, label))
 # our model.
 
 
-class GaussianLikelihood(tupak.likelihood.Likelihood):
-    def __init__(self, x, y, sigma, waveform_generator):
+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
         ----------
@@ -65,19 +62,25 @@ class GaussianLikelihood(tupak.likelihood.Likelihood):
             The data to analyse
         sigma: float
             The standard deviation of the noise
-        waveform_generator: `tupak.waveform_generator.WaveformGenerator`
-            An object which can generate the 'waveform', which in this case is
-            any arbitrary function
+        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.waveform_generator = waveform_generator
-        self.parameters = waveform_generator.parameters
+        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.waveform_generator.time_domain_strain()
+        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))
 
@@ -86,18 +89,9 @@ class GaussianLikelihood(tupak.likelihood.Likelihood):
                        + self.N*np.log(2*np.pi*self.sigma**2))
 
 
-# Here, we make a `tupak` waveform_generator. In this case, of course, the
-# name doesn't make so much sense. But essentially this is an objects that
-# can generate a signal. We give it information on how to make the time series
-# and the model() we wrote earlier.
-
-waveform_generator = tupak.waveform_generator.WaveformGenerator(
-    time_duration=time_duration, sampling_frequency=sampling_frequency,
-    time_domain_source_model=model)
-
-# Now lets instantiate a version of out GravitationalWaveTransient, giving it
-# the time, data and waveform_generator
-likelihood = GaussianLikelihood(time, data, sigma, waveform_generator)
+# Now lets instantiate a version of our GaussianLikelihood, giving it
+# the time, data and signal model
+likelihood = GaussianLikelihood(time, data, sigma, model)
 
 # From hereon, the syntax is exactly equivalent to other tupak examples
 # We make a prior
@@ -106,8 +100,9 @@ priors['m'] = tupak.prior.Uniform(0, 5, 'm')
 priors['c'] = tupak.prior.Uniform(-2, 2, 'c')
 
 # And run sampler
-result = tupak.sampler.run_sampler(
+result = tupak.run_sampler(
     likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
     walks=10, injection_parameters=injection_parameters, outdir=outdir,
-    label=label, plot=True)
+    label=label)
+result.plot_corner()
 print(result)