Compare revisions

MoritzThomasHuebner · MoritzThomasHuebner · MoritzThomasHuebner · MoritzThomasHuebner · MoritzThomasHuebner · MoritzThomasHuebner
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,11 @@

 ## Unreleased

+### Changes
+- Renamed "prior" to "prior" in bilby.gw.likelihood.GravtitationalWaveTransient
+  for consistency with bilby.core. **WARNING**: This will break scripts which
+  use marginalization.
+
 ## [0.3.3] 2018-11-08

 Changes currently on master, but not under a tag.
@@ -14,6 +19,7 @@ Changes currently on master, but not under a tag.
 - Added method to result to get injection recovery credible levels
 - Added function to generate a pp-plot from many results to core/result.py
 - Fixed a bug which caused `Interferometer.detector_tensor` not to update when `latitude`, `longitude`, `xarm_azimuth`, `yarm_azimuth`, `xarm_tilt`, `yarm_tilt` were updated.
+- Extracted time and frequency series behaviour from `WaveformGenerator` and `InterferometerStrainData` and moved it to `series.gw.CoupledTimeAndFrequencySeries`

 ### Changes
 - Switch the ordering the key-word arguments in `result.read_in_result` to put

--- a/bilby/core/likelihood.py
+++ b/bilby/core/likelihood.py
 from __future__ import division, print_function
-
 import copy
+
 import numpy as np
 from scipy.special import gammaln
+
 from .utils import infer_parameters_from_function


@@ -47,6 +48,20 @@ class Likelihood(object):
        """
        return self.log_likelihood() - self.noise_log_likelihood()

+    @property
+    def meta_data(self):
+        try:
+            return self._meta_data
+        except AttributeError:
+            return None
+
+    @meta_data.setter
+    def meta_data(self, meta_data):
+        if isinstance(meta_data, dict):
+            self._meta_data = meta_data
+        else:
+            raise ValueError("The meta_data must be an instance of dict")
+

 class Analytical1DLikelihood(Likelihood):
    """

--- a/bilby/core/prior.py
+++ b/bilby/core/prior.py
 from __future__ import division

-import numpy as np
-from scipy.interpolate import interp1d
-from scipy.integrate import cumtrapz
-from scipy.special import erf, erfinv
-import scipy.stats
 import os
 from collections import OrderedDict
 from future.utils import iteritems

-from .utils import logger, infer_args_from_method
-from . import utils
+import numpy as np
+import scipy.stats
+from scipy.integrate import cumtrapz
+from scipy.interpolate import interp1d
+from scipy.special import erf, erfinv
+
+# Keep import bilby statement, it is necessary for some eval() statements
 import bilby  # noqa
+from .utils import logger, infer_args_from_method, check_directory_exists_and_if_not_mkdir


 class PriorDict(OrderedDict):
@@ -50,7 +51,7 @@ class PriorDict(OrderedDict):
            Output file naming scheme
        """

-        utils.check_directory_exists_and_if_not_mkdir(outdir)
+        check_directory_exists_and_if_not_mkdir(outdir)
        prior_file = os.path.join(outdir, "{}.prior".format(label))
        logger.debug("Writing priors to {}".format(prior_file))
        with open(prior_file, "w") as outfile:
@@ -121,8 +122,7 @@ class PriorDict(OrderedDict):
        likelihood: bilby.likelihood.GravitationalWaveTransient instance
            Used to infer the set of parameters to fill the prior with
        default_priors_file: str, optional
-            If given, a file containing the default priors; otherwise defaults
-            to the bilby defaults for a binary black hole.
+            If given, a file containing the default priors.


        Returns
@@ -257,8 +257,7 @@ def create_default_prior(name, default_priors_file=None):
    name: str
        Parameter name
    default_priors_file: str, optional
-        If given, a file containing the default priors; otherwise defaults to
-        the bilby defaults for a binary black hole.
+        If given, a file containing the default priors.

    Return
    ------
@@ -1219,17 +1218,25 @@ class StudentT(Prior):


 class Beta(Prior):
-    def __init__(self, alpha, beta, name=None, latex_label=None, unit=None):
+    def __init__(self, alpha, beta, minimum=0, maximum=1, name=None,
+                 latex_label=None, unit=None):
        """Beta distribution

        https://en.wikipedia.org/wiki/Beta_distribution

+        This wraps around
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.beta.html
+
        Parameters
        ----------
        alpha: float
            first shape parameter
        beta: float
            second shape parameter
+        minimum: float
+            See superclass
+        maximum: float
+            See superclass
        name: str
            See superclass
        latex_label: str
@@ -1238,7 +1245,7 @@ class Beta(Prior):
            See superclass

        """
-        Prior.__init__(self, minimum=0., maximum=1., name=name,
+        Prior.__init__(self, minimum=minimum, maximum=maximum, name=name,
                       latex_label=latex_label, unit=unit)

        if alpha <= 0. or beta <= 0.:
@@ -1246,6 +1253,8 @@ class Beta(Prior):

        self.alpha = alpha
        self.beta = beta
+        self._loc = minimum
+        self._scale = maximum - minimum

    def rescale(self, val):
        """
@@ -1256,7 +1265,8 @@ class Beta(Prior):
        Prior.test_valid_for_rescaling(val)

        # use scipy distribution percentage point function (ppf)
-        return scipy.stats.beta.ppf(val, self.alpha, self.beta)
+        return scipy.stats.beta.ppf(
+            val, self.alpha, self.beta, loc=self._loc, scale=self._scale)

    def prob(self, val):
        """Return the prior probability of val.
@@ -1270,7 +1280,8 @@ class Beta(Prior):
        float: Prior probability of val
        """

-        spdf = scipy.stats.beta.pdf(val, self.alpha, self.beta)
+        spdf = scipy.stats.beta.pdf(
+            val, self.alpha, self.beta, loc=self._loc, scale=self._scale)
        if np.all(np.isfinite(spdf)):
            return spdf

@@ -1283,7 +1294,8 @@ class Beta(Prior):
            return 0.

    def ln_prob(self, val):
-        spdf = scipy.stats.beta.logpdf(val, self.alpha, self.beta)
+        spdf = scipy.stats.beta.logpdf(
+            val, self.alpha, self.beta, loc=self._loc, scale=self._scale)
        if np.all(np.isfinite(spdf)):
            return spdf


--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -2,13 +2,15 @@ from __future__ import division

 import os
 from distutils.version import LooseVersion
+from collections import OrderedDict, namedtuple
+
 import numpy as np
 import deepdish
 import pandas as pd
 import corner
+import scipy.stats
 import matplotlib
 import matplotlib.pyplot as plt
-from collections import OrderedDict, namedtuple

 from . import utils
 from .utils import (logger, infer_parameters_from_function,
@@ -74,7 +76,7 @@ class Result(object):
                 log_bayes_factor=np.nan, log_likelihood_evaluations=None,
                 sampling_time=None, nburn=None, walkers=None,
                 max_autocorrelation_time=None, parameter_labels=None,
-                 parameter_labels_with_unit=None):
+                 parameter_labels_with_unit=None, version=None):
        """ A class to store the results of the sampling run

        Parameters
@@ -110,6 +112,9 @@ class Result(object):
            The estimated maximum autocorrelation time for MCMC samplers
        parameter_labels, parameter_labels_with_unit: list
            Lists of the latex-formatted parameter labels
+        version: str,
+            Version information for software used to generate the result. Note,
+            this information is generated when the result object is initialized

        Note:
            All sampling output parameters, e.g. the samples themselves are
@@ -139,6 +144,7 @@ class Result(object):
        self.log_bayes_factor = log_bayes_factor
        self.log_likelihood_evaluations = log_likelihood_evaluations
        self.sampling_time = sampling_time
+        self.version = version
        self.max_autocorrelation_time = max_autocorrelation_time

    def __str__(self):
@@ -245,7 +251,19 @@ class Result(object):
    def posterior(self, posterior):
        self._posterior = posterior

+    @property
+    def version(self):
+        return self._version
+
+    @version.setter
+    def version(self, version):
+        if version is None:
+            self._version = 'bilby={}'.format(utils.get_version_information())
+        else:
+            self._version = version
+
    def _get_save_data_dictionary(self):
+        # This list defines all the parameters saved in the result object
        save_attrs = [
            'label', 'outdir', 'sampler', 'log_evidence', 'log_evidence_err',
            'log_noise_evidence', 'log_bayes_factor', 'priors', 'posterior',
@@ -253,7 +271,7 @@ class Result(object):
            'fixed_parameter_keys', 'sampling_time', 'sampler_kwargs',
            'log_likelihood_evaluations', 'samples', 'nested_samples',
            'walkers', 'nburn', 'parameter_labels',
-            'parameter_labels_with_unit']
+            'parameter_labels_with_unit', 'version']
        dictionary = OrderedDict()
        for attr in save_attrs:
            try:
@@ -514,36 +532,45 @@ class Result(object):
        """
        if isinstance(parameters, dict):
            plot_parameter_keys = list(parameters.keys())
-            truths = list(parameters.values())
+            truths = parameters
        elif parameters is None:
-            plot_parameter_keys = self.injection_parameters.keys()
-            truths = [self.injection_parameters.get(key, None) for key
-                      in plot_parameter_keys]
+            plot_parameter_keys = self.posterior.keys()
+            if self.injection_parameters is None:
+                truths = dict()
+            else:
+                truths = self.injection_parameters
        else:
            plot_parameter_keys = list(parameters)
-            truths = [self.injection_parameters.get(key, None) for key
-                      in plot_parameter_keys]
+            if self.injection_parameters is None:
+                truths = dict()
+            else:
+                truths = self.injection_parameters

        if file_base_name is None:
            file_base_name = '{}/{}_1d/'.format(self.outdir, self.label)
            check_directory_exists_and_if_not_mkdir(file_base_name)

        if priors is True:
-            priors = getattr(self, 'priors', False)
-        elif isinstance(priors, dict) or priors in [False, None]:
+            priors = getattr(self, 'priors', dict())
+        elif isinstance(priors, dict):
            pass
+        elif priors in [False, None]:
+            priors = dict()
        else:
            raise ValueError('Input priors={} not understood'.format(priors))

        for i, key in enumerate(plot_parameter_keys):
            if not isinstance(self.posterior[key].values[0], float):
                continue
+            prior = priors.get(key, None)
+            truth = truths.get(key, None)
            for cumulative in [False, True]:
-                self.plot_single_density(
-                    key, prior=priors[i], cumulative=cumulative, title=titles,
-                    truth=truths[i], save=True, file_base_name=file_base_name,
+                fig = self.plot_single_density(
+                    key, prior=prior, cumulative=cumulative, title=titles,
+                    truth=truth, save=True, file_base_name=file_base_name,
                    bins=bins, label_fontsize=label_fontsize, dpi=dpi,
                    title_fontsize=title_fontsize, quantiles=quantiles)
+                plt.close(fig)

    def plot_corner(self, parameters=None, priors=None, titles=True, save=True,
                    filename=None, dpi=300, **kwargs):
@@ -911,6 +938,44 @@ class Result(object):
                    return np.all(A == B)
        return False

+    @property
+    def kde(self):
+        """ Kernel density estimate built from the stored posterior
+
+        Uses `scipy.stats.gaussian_kde` to generate the kernel density
+        """
+        try:
+            return self._kde
+        except AttributeError:
+            self._kde = scipy.stats.gaussian_kde(
+                self.posterior[self.search_parameter_keys].values.T)
+            return self._kde
+
+    def posterior_probability(self, sample):
+        """ Calculate the posterior probabily for a new sample
+
+        This queries a Kernel Density Estimate of the posterior to calculate
+        the posterior probability density for the new sample.
+
+        Parameters
+        ----------
+        sample: dict, or list of dictionaries
+            A dictionary containing all the keys from
+            self.search_parameter_keys and corresponding values at which to
+            calculate the posterior probability
+
+        Returns
+        -------
+        p: array-like,
+            The posterior probability of the sample
+
+        """
+        if isinstance(sample, dict):
+            sample = [sample]
+        ordered_sample = [[s[key] for key in self.search_parameter_keys]
+                          for s in sample]
+        return self.kde(ordered_sample)
+

 def plot_multiple(results, filename=None, labels=None, colours=None,
                  save=True, evidences=False, **kwargs):

--- a/bilby/core/sampler/__init__.py
+++ b/bilby/core/sampler/__init__.py
@@ -109,6 +109,11 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',

    priors.fill_priors(likelihood, default_priors_file=default_priors_file)

+    # Generate the meta-data if not given and append the likelihood meta_data
+    if meta_data is None:
+        meta_data = dict()
+    meta_data['likelihood'] = likelihood.meta_data
+
    if isinstance(sampler, Sampler):
        pass
    elif isinstance(sampler, str):

--- a/bilby/core/sampler/base_sampler.py
+++ b/bilby/core/sampler/base_sampler.py
 from __future__ import absolute_import
 import datetime
 import numpy as np
+
 from pandas import DataFrame
+
 from ..utils import logger, command_line_args
 from ..prior import Prior, PriorDict
 from ..result import Result, read_in_result
@@ -143,9 +145,9 @@ class Sampler(object):
        external_sampler_name = self.__class__.__name__.lower()
        try:
            self.external_sampler = __import__(external_sampler_name)
-        except ImportError:
-            raise ImportError(
-                "Sampler {} not installed on this system".format(external_sampler_name))
+        except (ImportError, SystemExit, ModuleNotFoundError):
+            raise SamplerNotInstalledError(
+                "Sampler {} is not installed on this system".format(external_sampler_name))

    def _verify_kwargs_against_default_kwargs(self):
        """
@@ -473,3 +475,15 @@ class MCMCSampler(Sampler):
        except emcee.autocorr.AutocorrError as e:
            self.result.max_autocorrelation_time = None
            logger.info("Unable to calculate autocorr time: {}".format(e))
+
+
+class Error(Exception):
+    """ Base class for all exceptions raised by this module """
+
+
+class SamplerError(Error):
+    """ Base class for Error related to samplers in this module """
+
+
+class SamplerNotInstalledError(SamplerError):
+    """ Base class for Error raised by not installed samplers """
--- a/bilby/core/sampler/cpnest.py
+++ b/bilby/core/sampler/cpnest.py
 from __future__ import absolute_import
+
 import numpy as np
 from pandas import DataFrame
-from ..utils import logger, check_directory_exists_and_if_not_mkdir
+
 from .base_sampler import NestedSampler
+from ..utils import logger, check_directory_exists_and_if_not_mkdir


 class Cpnest(NestedSampler):

--- a/bilby/core/sampler/dynesty.py
+++ b/bilby/core/sampler/dynesty.py
@@ -2,9 +2,11 @@ from __future__ import absolute_import

 import os
 import sys
+
 import numpy as np
 from pandas import DataFrame
 from deepdish.io import load, save
+
 from ..utils import logger, check_directory_exists_and_if_not_mkdir
 from .base_sampler import Sampler, NestedSampler

@@ -391,6 +393,7 @@ class Dynesty(NestedSampler):
        self.sampler.saved_scale = [self.sampler.saved_scale[-1]]

    def generate_trace_plots(self, dynesty_results):
+        check_directory_exists_and_if_not_mkdir(self.outdir)
        filename = '{}/{}_trace.png'.format(self.outdir, self.label)
        logger.debug("Writing trace plot to {}".format(filename))
        from dynesty import plotting as dyplot

--- a/bilby/core/sampler/emcee.py
+++ b/bilby/core/sampler/emcee.py
 from __future__ import absolute_import, print_function
+
 import numpy as np
 from pandas import DataFrame
+
 from ..utils import logger, get_progress_bar
 from .base_sampler import MCMCSampler


--- a/bilby/core/sampler/nestle.py
+++ b/bilby/core/sampler/nestle.py
@@ -2,6 +2,7 @@ from __future__ import absolute_import

 import numpy as np
 from pandas import DataFrame
+
 from .base_sampler import NestedSampler



--- a/bilby/core/sampler/ptemcee.py
+++ b/bilby/core/sampler/ptemcee.py
 from __future__ import absolute_import, division, print_function

 import numpy as np
+
 from ..utils import get_progress_bar, logger
 from . import Emcee

@@ -57,7 +58,6 @@ class Ptemcee(Emcee):
    def run_sampler(self):
        import ptemcee
        tqdm = get_progress_bar()
-
        sampler = ptemcee.Sampler(dim=self.ndim, logl=self.log_likelihood,
                                  logp=self.log_prior, **self.sampler_init_kwargs)
        self.pos0 = [[self.get_random_draw_from_prior()

--- a/bilby/core/sampler/pymc3.py
+++ b/bilby/core/sampler/pymc3.py
 from __future__ import absolute_import, print_function

 from collections import OrderedDict
+
 import numpy as np

 from ..utils import derivatives, logger, infer_args_from_method
-from ..prior import Prior
+from ..prior import Prior, DeltaFunction, Sine, Cosine, PowerLaw
 from ..result import Result
 from .base_sampler import Sampler, MCMCSampler
+from ..likelihood import GaussianLikelihood, PoissonLikelihood, ExponentialLikelihood, \
+    StudentTLikelihood
+from ...gw.likelihood import BasicGravitationalWaveTransient, GravitationalWaveTransient


 class Pymc3(MCMCSampler):
@@ -32,11 +36,13 @@ class Pymc3(MCMCSampler):
        chains.
    step: str, dict
        Provide a step method name, or dictionary of step method names keyed to
-        particular variable names (these are case insensitive). If no method is
-        provided for any particular variable then PyMC3 will automatically
-        decide upon a default, with the first option being the NUTS sampler.
-        The currently allowed methods are 'NUTS', 'HamiltonianMC',
-        'Metropolis', 'BinaryMetropolis', 'BinaryGibbsMetropolis', 'Slice', and
+        particular variable names (these are case insensitive). If passing a
+        dictionary of methods, the values keyed on particular variables can be
+        lists of methods to form compound steps. If no method is provided for
+        any particular variable then PyMC3 will automatically decide upon a
+        default, with the first option being the NUTS sampler. The currently
+        allowed methods are 'NUTS', 'HamiltonianMC', 'Metropolis',
+        'BinaryMetropolis', 'BinaryGibbsMetropolis', 'Slice', and
        'CategoricalGibbsMetropolis'. Note: you cannot provide a PyMC3 step
        method function itself here as it is outside of the model context
        manager.
@@ -56,14 +62,28 @@ class Pymc3(MCMCSampler):
        live_plot_kwargs=None, compute_convergence_checks=True, use_mmap=False)

    def __init__(self, likelihood, priors, outdir='outdir', label='label',
-                 use_ratio=False, plot=False, draws=1000,
+                 use_ratio=False, plot=False,
                 skip_import_verification=False, **kwargs):
        Sampler.__init__(self, likelihood, priors, outdir=outdir, label=label,
                         use_ratio=use_ratio, plot=plot,
                         skip_import_verification=skip_import_verification, **kwargs)
-        self.draws = draws
+        self.draws = self.__kwargs['draws']
        self.chains = self.__kwargs['chains']

+    @staticmethod
+    def _import_external_sampler():
+        import pymc3
+        from pymc3.sampling import STEP_METHODS
+        from pymc3.theanof import floatX
+        return pymc3, STEP_METHODS, floatX
+
+    @staticmethod
+    def _import_theano():
+        import theano  # noqa
+        import theano.tensor as tt
+        from theano.compile.ops import as_op  # noqa
+        return theano, tt, as_op
+
    def _verify_parameters(self):
        """
        Change `_verify_parameters()` to just pass, i.e., don't try and
@@ -217,8 +237,6 @@ class Pymc3(MCMCSampler):
        Map the bilby delta function prior to a single value for PyMC3.
        """

-        from ..prior import DeltaFunction
-
        # check prior is a DeltaFunction
        if isinstance(self.priors[key], DeltaFunction):
            return self.priors[key].peak
@@ -230,17 +248,10 @@ class Pymc3(MCMCSampler):
        Map the bilby Sine prior to a PyMC3 style function
        """

-        from ..prior import Sine
-
        # check prior is a Sine
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
        if isinstance(self.priors[key], Sine):
-            import pymc3
-
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")

            class Pymc3Sine(pymc3.Continuous):
                def __init__(self, lower=0., upper=np.pi):
@@ -277,18 +288,10 @@ class Pymc3(MCMCSampler):
        Map the bilby Cosine prior to a PyMC3 style function
        """

-        from ..prior import Cosine
-
        # check prior is a Cosine
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
        if isinstance(self.priors[key], Cosine):
-            import pymc3
-
-            # import theano
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")

            class Pymc3Cosine(pymc3.Continuous):
                def __init__(self, lower=-np.pi / 2., upper=np.pi / 2.):
@@ -324,23 +327,15 @@ class Pymc3(MCMCSampler):
        Map the bilby PowerLaw prior to a PyMC3 style function
        """

-        from ..prior import PowerLaw
-
        # check prior is a PowerLaw
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
        if isinstance(self.priors[key], PowerLaw):
-            import pymc3

            # check power law is set
            if not hasattr(self.priors[key], 'alpha'):
                raise AttributeError("No 'alpha' attribute set for PowerLaw prior")

-            # import theano
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")
-
            if self.priors[key].alpha < -1.:
                # use Pareto distribution
                palpha = -(1. + self.priors[key].alpha)
@@ -385,11 +380,8 @@ class Pymc3(MCMCSampler):
            raise ValueError("Prior for '{}' is not a Power Law".format(key))

    def run_sampler(self):
-        import pymc3
        # set the step method
-
-        from pymc3.sampling import STEP_METHODS
-
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
        step_methods = {m.__name__.lower(): m.__name__ for m in STEP_METHODS}
        if 'step' in self.__kwargs:
            self.step_method = self.__kwargs.pop('step')
@@ -403,13 +395,26 @@ class Pymc3(MCMCSampler):
                    if key not in self.__search_parameter_keys:
                        raise ValueError("Setting a step method for an unknown parameter '{}'".format(key))
                    else:
-                        if self.step_method[key].lower() not in step_methods:
-                            raise ValueError("Using invalid step method '{}'".format(self.step_method[key]))
+                        # check if using a compound step (a list of step
+                        # methods for a particular parameter)
+                        if isinstance(self.step_method[key], list):
+                            sms = self.step_method[key]
+                        else:
+                            sms = [self.step_method[key]]
+                        for sm in sms:
+                            if sm.lower() not in step_methods:
+                                raise ValueError("Using invalid step method '{}'".format(self.step_method[key]))
            else:
-                self.step_method = self.step_method.lower()
+                # check if using a compound step (a list of step
+                # methods for a particular parameter)
+                if isinstance(self.step_method, list):
+                    sms = self.step_method
+                else:
+                    sms = [self.step_method]

-                if self.step_method not in step_methods:
-                    raise ValueError("Using invalid step method '{}'".format(self.step_method))
+                for i in range(len(sms)):
+                    if sms[i].lower() not in step_methods:
+                        raise ValueError("Using invalid step method '{}'".format(sms[i]))
        else:
            self.step_method = None

@@ -419,18 +424,6 @@ class Pymc3(MCMCSampler):
        # set the prior
        self.set_prior()

-        # set the step method
-        if isinstance(self.step_method, (dict, OrderedDict)):
-            # create list of step methods (any not given will default to NUTS)
-            self.kwargs['step'] = []
-            with self.pymc3_model:
-                for key in self.step_method:
-                    curmethod = self.step_method[key].lower()
-                    self.kwargs['step'].append(pymc3.__dict__[step_methods[curmethod]]([self.pymc3_priors[key]]))
-        else:
-            with self.pymc3_model:
-                self.kwargs['step'] = None if self.step_method is None else pymc3.__dict__[step_methods[self.step_method]]()
-
        # if a custom log_likelihood function requires a `sampler` argument
        # then use that log_likelihood function, with the assumption that it
        # takes in a Pymc3 Sampler, with a pymc3_model attribute, and defines
@@ -442,6 +435,107 @@ class Pymc3(MCMCSampler):
            # set the likelihood function from predefined functions
            self.set_likelihood()

+        # get the step method keyword arguments
+        step_kwargs = self.kwargs.pop('step_kwargs')
+        nuts_kwargs = self.kwargs.pop('nuts_kwargs')
+        methodslist = []
+
+        # set the step method
+        if isinstance(self.step_method, (dict, OrderedDict)):
+            # create list of step methods (any not given will default to NUTS)
+            self.kwargs['step'] = []
+            with self.pymc3_model:
+                for key in self.step_method:
+                    # check for a compound step list
+                    if isinstance(self.step_method[key], list):
+                        for sms in self.step_method[key]:
+                            curmethod = sms.lower()
+                            methodslist.append(curmethod)
+                            args = {}
+                            if curmethod == 'nuts':
+                                if nuts_kwargs is not None:
+                                    args = nuts_kwargs
+                                elif step_kwargs is not None:
+                                    args = step_kwargs.pop('nuts', {})
+                                    # add values into nuts_kwargs
+                                    nuts_kwargs = args
+                                else:
+                                    args = {}
+                            else:
+                                if step_kwargs is not None:
+                                    args = step_kwargs.get(curmethod, {})
+                                else:
+                                    args = {}
+                            self.kwargs['step'].append(pymc3.__dict__[step_methods[curmethod]](vars=[self.pymc3_priors[key]], **args))
+                    else:
+                        curmethod = self.step_method[key].lower()
+                        methodslist.append(curmethod)
+                        args = {}
+                        if curmethod == 'nuts':
+                            if nuts_kwargs is not None:
+                                args = nuts_kwargs
+                            elif step_kwargs is not None:
+                                args = step_kwargs.pop('nuts', {})
+                                # add values into nuts_kwargs
+                                nuts_kwargs = args
+                            else:
+                                args = {}
+                        else:
+                            if step_kwargs is not None:
+                                args = step_kwargs.get(curmethod, {})
+                            else:
+                                args = {}
+                        self.kwargs['step'].append(pymc3.__dict__[step_methods[curmethod]](vars=[self.pymc3_priors[key]], **args))
+        else:
+            with self.pymc3_model:
+                # check for a compound step list
+                if isinstance(self.step_method, list):
+                    compound = []
+                    for sms in self.step_method:
+                        curmethod = sms.lower()
+                        methodslist.append(curmethod)
+                        args = {}
+                        if curmethod == 'nuts':
+                            if nuts_kwargs is not None:
+                                args = nuts_kwargs
+                            elif step_kwargs is not None:
+                                args = step_kwargs.pop('nuts', {})
+                                # add values into nuts_kwargs
+                                nuts_kwargs = args
+                            else:
+                                args = {}
+                        else:
+                            args = step_kwargs.get(curmethod, {})
+                        compound.append(pymc3.__dict__[step_methods[curmethod]](**args))
+                        self.kwargs['step'] = compound
+                else:
+                    self.kwargs['step'] = None
+                    if self.step_method is not None:
+                        curmethod = self.step_method.lower()
+                        methodslist.append(curmethod)
+                        args = {}
+                        if curmethod == 'nuts':
+                            if nuts_kwargs is not None:
+                                args = nuts_kwargs
+                            elif step_kwargs is not None:
+                                args = step_kwargs.pop('nuts', {})
+                                # add values into nuts_kwargs
+                                nuts_kwargs = args
+                            else:
+                                args = {}
+                        else:
+                            args = step_kwargs.get(curmethod, {})
+                        self.kwargs['step'] = pymc3.__dict__[step_methods[curmethod]](**args)
+                    else:
+                        # re-add step_kwargs if no step methods are set
+                        self.kwargs['step_kwargs'] = step_kwargs
+
+        # check whether only NUTS step method has been assigned
+        if np.all([sm.lower() == 'nuts' for sm in methodslist]):
+            # in this case we can let PyMC3 autoinitialise NUTS, so remove the step methods and re-add nuts_kwargs
+            self.kwargs['step'] = None
+            self.kwargs['nuts_kwargs'] = nuts_kwargs
+
        with self.pymc3_model:
            # perform the sampling
            trace = pymc3.sample(**self.kwargs)
@@ -470,8 +564,7 @@ class Pymc3(MCMCSampler):
        self.setup_prior_mapping()

        self.pymc3_priors = OrderedDict()
-
-        import pymc3
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()

        # set the parameter prior distributions (in the model context manager)
        with self.pymc3_model:
@@ -534,18 +627,10 @@ class Pymc3(MCMCSampler):
        Convert any bilby likelihoods to PyMC3 distributions.
        """

-        try:
-            import theano  # noqa
-            import theano.tensor as tt
-            from theano.compile.ops import as_op  # noqa
-        except ImportError:
-            raise ImportError("Could not import theano")
-
-        from ..likelihood import GaussianLikelihood, PoissonLikelihood, ExponentialLikelihood, \
-            StudentTLikelihood
-        from ...gw.likelihood import BasicGravitationalWaveTransient, GravitationalWaveTransient
-
        # create theano Op for the log likelihood if not using a predefined model
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
+
        class LogLike(tt.Op):

            itypes = [tt.dvector]
@@ -605,8 +690,6 @@ class Pymc3(MCMCSampler):

                outputs[0][0] = grads

-        import pymc3
-
        with self.pymc3_model:
            #  check if it is a predefined likelhood function
            if isinstance(self.likelihood, GaussianLikelihood):

--- a/bilby/core/sampler/pymultinest.py
+++ b/bilby/core/sampler/pymultinest.py
 from __future__ import absolute_import

 from ..utils import check_directory_exists_and_if_not_mkdir
-from .base_sampler import NestedSampler
 from ..utils import logger
+from .base_sampler import NestedSampler


 class Pymultinest(NestedSampler):
@@ -56,7 +56,7 @@ class Pymultinest(NestedSampler):
        https://github.com/JohannesBuchner/PyMultiNest/issues/115
        """
        if not self.kwargs['outputfiles_basename']:
-            self.kwargs['outputfiles_basename'] =\
+            self.kwargs['outputfiles_basename'] = \
                '{}/pm_{}/'.format(self.outdir, self.label)
        if self.kwargs['outputfiles_basename'].endswith('/') is False:
            self.kwargs['outputfiles_basename'] = '{}/'.format(

--- a/bilby/core/utils.py
+++ b/bilby/core/utils.py
 from __future__ import division
+
 import logging
 import os
-import numpy as np
 from math import fmod
 import argparse
 import traceback
 import inspect
 import subprocess

+import numpy as np
+
 logger = logging.getLogger('bilby')

 # Constants

--- a/bilby/gw/calibration.py
+++ b/bilby/gw/calibration.py
@@ -49,6 +49,9 @@ class Recalibrate(object):
        self.params.update({key[len(self.prefix):]: params[key] for key in params
                            if self.prefix in key})

+    def __eq__(self, other):
+        return self.__dict__ == other.__dict__
+

 class CubicSpline(Recalibrate):


--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
 from __future__ import division
+
 import numpy as np
 from pandas import DataFrame

@@ -181,6 +182,16 @@ def convert_to_lal_binary_black_hole_parameters(parameters):
                converted_parameters['mass_2'] /\
                converted_parameters['mass_ratio']

+    for idx in ['1', '2']:
+        key = 'chi_{}'.format(idx)
+        if key in original_keys:
+            converted_parameters['a_{}'.format(idx)] = abs(
+                converted_parameters[key])
+            converted_parameters['cos_tilt_{}'.format(idx)] = \
+                np.sign(converted_parameters[key])
+            converted_parameters['phi_jl'] = 0.0
+            converted_parameters['phi_12'] = 0.0
+
    for angle in ['tilt_1', 'tilt_2', 'iota']:
        cos_angle = str('cos_' + angle)
        if cos_angle in converted_parameters.keys():
@@ -262,6 +273,18 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
            * converted_parameters['mass_1']**5\
            / converted_parameters['mass_2']**5

+    for idx in ['1', '2']:
+        mag = 'a_{}'.format(idx)
+        if mag in original_keys:
+            tilt = 'tilt_{}'.format(idx)
+            if tilt in original_keys:
+                converted_parameters['chi_{}'.format(idx)] = (
+                    converted_parameters[mag] *
+                    np.cos(converted_parameters[tilt]))
+            else:
+                converted_parameters['chi_{}'.format(idx)] = (
+                    converted_parameters[mag])
+
    added_keys = [key for key in converted_parameters.keys()
                  if key not in original_keys]


--- a/bilby/gw/detector.py
+++ b/bilby/gw/detector.py
 from __future__ import division, print_function, absolute_import

 import os
+import sys

 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.signal.windows import tukey
 from scipy.interpolate import interp1d
+import deepdish as dd

 from . import utils as gwutils
 from ..core import utils
@@ -38,12 +40,12 @@ class InterferometerList(list):

        list.__init__(self)
        if type(interferometers) == str:
-            raise ValueError("Input must not be a string")
+            raise TypeError("Input must not be a string")
        for ifo in interferometers:
            if type(ifo) == str:
                ifo = get_empty_interferometer(ifo)
            if type(ifo) not in [Interferometer, TriangularInterferometer]:
-                raise ValueError("Input list of interferometers are not all Interferometer objects")
+                raise TypeError("Input list of interferometers are not all Interferometer objects")
            else:
                self.append(ifo)
        self._check_interferometers()
@@ -108,7 +110,7 @@ class InterferometerList(list):
        """
        if injection_polarizations is None:
            if waveform_generator is not None:
-                injection_polarizations =\
+                injection_polarizations = \
                    waveform_generator.frequency_domain_strain(parameters)
            else:
                raise ValueError(
@@ -207,6 +209,54 @@ class InterferometerList(list):
        super(InterferometerList, self).insert(index, interferometer)
        self._check_interferometers()

+    @property
+    def meta_data(self):
+        """ Dictionary of the per-interferometer meta_data """
+        return {interferometer.name: interferometer.meta_data
+                for interferometer in self}
+
+    @staticmethod
+    def _hdf5_filename_from_outdir_label(outdir, label):
+        return os.path.join(outdir, label + '.h5')
+
+    def to_hdf5(self, outdir='outdir', label='ifo_list'):
+        """ Saves the object to a hdf5 file
+
+        Parameters
+        ----------
+        outdir: str, optional
+            Output directory name of the file
+        label: str, optional
+            Output file name, is 'ifo_list' if not given otherwise. A list of
+            the included interferometers will be appended.
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        label = label + '_' + ''.join(ifo.name for ifo in self)
+        utils.check_directory_exists_and_if_not_mkdir(outdir)
+        dd.io.save(self._hdf5_filename_from_outdir_label(outdir, label), self)
+
+    @classmethod
+    def from_hdf5(cls, filename=None):
+        """ Loads in an InterferometerList object from an hdf5 file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        res = dd.io.load(filename)
+        if res.__class__ == list:
+            res = cls(res)
+        if res.__class__ != cls:
+            raise TypeError('The loaded object is not a InterferometerList')
+        return res
+

 class InterferometerStrainData(object):
    """ Strain data for an interferometer """
@@ -242,6 +292,21 @@ class InterferometerStrainData(object):
        self._time_domain_strain = None
        self._time_array = None

+    def __eq__(self, other):
+        if self.minimum_frequency == other.minimum_frequency \
+                and self.maximum_frequency == other.maximum_frequency \
+                and self.roll_off == other.roll_off \
+                and self.window_factor == other.window_factor \
+                and self.sampling_frequency == other.sampling_frequency \
+                and self.duration == other.duration \
+                and self.start_time == other.start_time \
+                and np.array_equal(self.time_array, other.time_array) \
+                and np.array_equal(self.frequency_array, other.frequency_array) \
+                and np.array_equal(self.frequency_domain_strain, other.frequency_domain_strain) \
+                and np.array_equal(self.time_domain_strain, other.time_domain_strain):
+            return True
+        return False
+
    def time_within_data(self, time):
        """ Check if time is within the data span

@@ -807,6 +872,23 @@ class Interferometer(object):
        self._strain_data = InterferometerStrainData(
            minimum_frequency=minimum_frequency,
            maximum_frequency=maximum_frequency)
+        self.meta_data = dict()
+
+    def __eq__(self, other):
+        if self.name == other.name and \
+                self.length == other.length and \
+                self.latitude == other.latitude and \
+                self.longitude == other.longitude and \
+                self.elevation == other.elevation and \
+                self.xarm_azimuth == other.xarm_azimuth and \
+                self.xarm_tilt == other.xarm_tilt and \
+                self.yarm_azimuth == other.yarm_azimuth and \
+                self.yarm_tilt == other.yarm_tilt and \
+                self.power_spectral_density.__eq__(other.power_spectral_density) and \
+                self.calibration_model == other.calibration_model and \
+                self.strain_data == other.strain_data:
+            return True
+        return False

    def __repr__(self):
        return self.__class__.__name__ + '(name=\'{}\', power_spectral_density={}, minimum_frequency={}, ' \
@@ -1242,7 +1324,7 @@ class Interferometer(object):

        if injection_polarizations is None:
            if waveform_generator is not None:
-                injection_polarizations =\
+                injection_polarizations = \
                    waveform_generator.frequency_domain_strain(parameters)
            else:
                raise ValueError(
@@ -1271,12 +1353,16 @@ class Interferometer(object):
                sampling_frequency=self.strain_data.sampling_frequency,
                duration=self.strain_data.duration,
                start_time=self.strain_data.start_time)
-        opt_snr = np.sqrt(self.optimal_snr_squared(signal=signal_ifo).real)
-        mf_snr = np.sqrt(self.matched_filter_snr_squared(signal=signal_ifo).real)
+
+        self.meta_data['optimal_SNR'] = (
+            np.sqrt(self.optimal_snr_squared(signal=signal_ifo)).real)
+        self.meta_data['matched_filter_SNR'] = (
+            np.sqrt(self.matched_filter_snr_squared(signal=signal_ifo)))
+        self.meta_data['parameters'] = parameters

        logger.info("Injected signal in {}:".format(self.name))
-        logger.info("  optimal SNR = {:.2f}".format(opt_snr))
-        logger.info("  matched filter SNR = {:.2f}".format(mf_snr))
+        logger.info("  optimal SNR = {:.2f}".format(self.meta_data['optimal_SNR']))
+        logger.info("  matched filter SNR = {:.2f}".format(self.meta_data['matched_filter_SNR']))
        for key in parameters:
            logger.info('  {} = {}'.format(key, parameters[key]))

@@ -1569,6 +1655,48 @@ class Interferometer(object):
            fig.savefig(
                '{}/{}_{}_time_domain_data.png'.format(outdir, self.name, label))

+    @staticmethod
+    def _hdf5_filename_from_outdir_label(outdir, label):
+        return os.path.join(outdir, label + '.h5')
+
+    def to_hdf5(self, outdir='outdir', label=None):
+        """ Save the object to a hdf5 file
+
+        Attributes
+        ----------
+        outdir: str, optional
+            Output directory name of the file, defaults to 'outdir'.
+        label: str, optional
+            Output file name, is self.name if not given otherwise.
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of Interferometer is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        if label is None:
+            label = self.name
+        utils.check_directory_exists_and_if_not_mkdir('outdir')
+        filename = self._hdf5_filename_from_outdir_label(outdir, label)
+        dd.io.save(filename, self)
+
+    @classmethod
+    def from_hdf5(cls, filename=None):
+        """ Loads in an Interferometer object from an hdf5 file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of Interferometer is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+
+        res = dd.io.load(filename)
+        if res.__class__ != cls:
+            raise TypeError('The loaded object is not an Interferometer')
+        return res
+

 class TriangularInterferometer(InterferometerList):

@@ -1637,6 +1765,15 @@ class PowerSpectralDensity(object):
        self.psd_file = psd_file
        self.asd_file = asd_file

+    def __eq__(self, other):
+        if self.psd_file == other.psd_file \
+                and self.asd_file == other.asd_file \
+                and np.array_equal(self.frequency_array, other.frequency_array) \
+                and np.array_equal(self.psd_array, other.psd_array) \
+                and np.array_equal(self.asd_array, other.asd_array):
+            return True
+        return False
+
    def __repr__(self):
        if self.asd_file is not None or self.psd_file is not None:
            return self.__class__.__name__ + '(psd_file=\'{}\', asd_file=\'{}\')' \
@@ -1758,12 +1895,12 @@ class PowerSpectralDensity(object):

    @property
    def asd_file(self):
-        return self.__asd_file
+        return self._asd_file

    @asd_file.setter
    def asd_file(self, asd_file):
        asd_file = self.__validate_file_name(file=asd_file)
-        self.__asd_file = asd_file
+        self._asd_file = asd_file
        if asd_file is not None:
            self.__import_amplitude_spectral_density()
            self.__check_file_was_asd_file()
@@ -1777,12 +1914,12 @@ class PowerSpectralDensity(object):

    @property
    def psd_file(self):
-        return self.__psd_file
+        return self._psd_file

    @psd_file.setter
    def psd_file(self, psd_file):
        psd_file = self.__validate_file_name(file=psd_file)
-        self.__psd_file = psd_file
+        self._psd_file = psd_file
        if psd_file is not None:
            self.__import_power_spectral_density()
            self.__check_file_was_psd_file()
@@ -2144,16 +2281,16 @@ def load_data_from_cache_file(
            frame_duration = float(frame_duration)
            if frame_name[:4] == 'file':
                frame_name = frame_name[16:]
-            if not data_set & (frame_start < segment_start) &\
+            if not data_set & (frame_start < segment_start) & \
                    (segment_start < frame_start + frame_duration):
                ifo.set_strain_data_from_frame_file(
                    frame_name, 4096, segment_duration,
                    start_time=segment_start,
                    channel=channel_name, buffer_time=0)
                data_set = True
-            if not psd_set & (frame_start < psd_start) &\
+            if not psd_set & (frame_start < psd_start) & \
                    (psd_start + psd_duration < frame_start + frame_duration):
-                ifo.power_spectral_density =\
+                ifo.power_spectral_density = \
                    PowerSpectralDensity.from_frame_file(
                        frame_name, psd_start_time=psd_start,
                        psd_duration=psd_duration,

--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -44,7 +44,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
    phase_marginalization: bool, optional
        If true, marginalize over phase in the likelihood.
        This is done analytically using a Bessel function.
-    prior: dict, optional
+    priors: dict, optional
        If given, used in the distance and phase marginalization.

    Returns
@@ -56,7 +56,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
    """

    def __init__(self, interferometers, waveform_generator, time_marginalization=False, distance_marginalization=False,
-                 phase_marginalization=False, prior=None):
+                 phase_marginalization=False, priors=None):

        self.waveform_generator = waveform_generator
        likelihood.Likelihood.__init__(self, dict())
@@ -64,32 +64,33 @@ class GravitationalWaveTransient(likelihood.Likelihood):
        self.time_marginalization = time_marginalization
        self.distance_marginalization = distance_marginalization
        self.phase_marginalization = phase_marginalization
-        self.prior = prior
+        self.priors = priors
        self._check_set_duration_and_sampling_frequency_of_waveform_generator()
+        self.meta_data = self.interferometers.meta_data

        if self.time_marginalization:
            self._check_prior_is_set(key='geocent_time')
            self._setup_time_marginalization()
-            prior['geocent_time'] = float(self.interferometers.start_time)
+            priors['geocent_time'] = float(self.interferometers.start_time)

        if self.phase_marginalization:
            self._check_prior_is_set(key='phase')
            self._bessel_function_interped = None
            self._setup_phase_marginalization()
-            prior['phase'] = float(0)
+            priors['phase'] = float(0)

        if self.distance_marginalization:
            self._check_prior_is_set(key='luminosity_distance')
-            self._distance_array = np.linspace(self.prior['luminosity_distance'].minimum,
-                                               self.prior['luminosity_distance'].maximum, int(1e4))
+            self._distance_array = np.linspace(self.priors['luminosity_distance'].minimum,
+                                               self.priors['luminosity_distance'].maximum, int(1e4))
            self._setup_distance_marginalization()
-            prior['luminosity_distance'] = float(self._ref_dist)
+            priors['luminosity_distance'] = float(self._ref_dist)

    def __repr__(self):
        return self.__class__.__name__ + '(interferometers={},\n\twaveform_generator={},\n\ttime_marginalization={}, ' \
-                                         'distance_marginalization={}, phase_marginalization={}, prior={})'\
+                                         'distance_marginalization={}, phase_marginalization={}, priors={})'\
            .format(self.interferometers, self.waveform_generator, self.time_marginalization,
-                    self.distance_marginalization, self.phase_marginalization, self.prior)
+                    self.distance_marginalization, self.phase_marginalization, self.priors)

    def _check_set_duration_and_sampling_frequency_of_waveform_generator(self):
        """ Check the waveform_generator has the same duration and
@@ -113,28 +114,28 @@ class GravitationalWaveTransient(likelihood.Likelihood):
            setattr(self.waveform_generator, attr, ifo_attr)

    def _check_prior_is_set(self, key):
-        if key not in self.prior or not isinstance(
-                self.prior[key], Prior):
+        if key not in self.priors or not isinstance(
+                self.priors[key], Prior):
            logger.warning(
                'Prior not provided for {}, using the BBH default.'.format(key))
            if key == 'geocent_time':
-                self.prior[key] = Uniform(
+                self.priors[key] = Uniform(
                    self.interferometers.start_time,
                    self.interferometers.start_time + self.interferometers.duration)
            else:
-                self.prior[key] = BBHPriorDict()[key]
+                self.priors[key] = BBHPriorDict()[key]

    @property
-    def prior(self):
+    def priors(self):
        return self.__prior

-    @prior.setter
-    def prior(self, prior):
-        if prior is not None:
-            self.__prior = prior.copy()
+    @priors.setter
+    def priors(self, priors):
+        if priors is not None:
+            self.__prior = priors.copy()
        elif any([self.time_marginalization, self.phase_marginalization,
                  self.distance_marginalization]):
-            raise ValueError("You can't use a marginalized likelihood without specifying a prior")
+            raise ValueError("You can't use a marginalized likelihood without specifying a priors")
        else:
            self.__prior = None

@@ -229,7 +230,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):

    @property
    def _ref_dist(self):
-        """ Smallest distance contained in prior """
+        """ Smallest distance contained in priors """
        return self._distance_array[0]

    @property
@@ -253,7 +254,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):

    def _create_lookup_table(self):
        """ Make the lookup table """
-        self.distance_prior_array = np.array([self.prior['luminosity_distance'].prob(distance)
+        self.distance_prior_array = np.array([self.priors['luminosity_distance'].prob(distance)
                                              for distance in self._distance_array])
        logger.info('Building lookup table for distance marginalisation.')

@@ -284,9 +285,9 @@ class GravitationalWaveTransient(likelihood.Likelihood):
            self.interferometers.start_time + np.linspace(
                0, self.interferometers.duration,
                int(self.interferometers.duration / 2 *
-                    self.waveform_generator.sampling_frequency) + 1)[1:]
+                    self.waveform_generator.sampling_frequency + 1))[1:]
        self.time_prior_array =\
-            self.prior['geocent_time'].prob(times) * delta_tc
+            self.priors['geocent_time'].prob(times) * delta_tc


 class BasicGravitationalWaveTransient(likelihood.Likelihood):

--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
 import os
+
 import numpy as np
 from scipy.interpolate import UnivariateSpline
+
 from ..core.prior import (PriorDict, Uniform, FromFile, Prior, DeltaFunction,
                          Gaussian, Interped)
 from ..core.utils import logger

--- a/bilby/gw/series.py
+++ b/bilby/gw/series.py
@@ -20,6 +20,12 @@ class CoupledTimeAndFrequencySeries(object):
        self.start_time = start_time
        self._frequency_array_updated = False
        self._time_array_updated = False
+        self._frequency_array = None
+        self._time_array = None
+
+    def __repr__(self):
+        return self.__class__.__name__ + '(duration={}, sampling_frequency={}, start_time={})'\
+            .format(self.duration, self.sampling_frequency, self.start_time)

    @property
    def frequency_array(self):
@@ -31,7 +37,7 @@ class CoupledTimeAndFrequencySeries(object):
        """
        if self._frequency_array_updated is False:
            if self.sampling_frequency and self.duration:
-                self.frequency_array = utils.create_frequency_series(
+                self._frequency_array = utils.create_frequency_series(
                    sampling_frequency=self.sampling_frequency,
                    duration=self.duration)
            else:
No results found