.gitlab-ci.yml

@@ -23,13 +23,15 @@ exitcode-jessie:
     - pip install coverage-badge
     - python setup.py install
     - coverage erase
+    - coverage run --include=tupak/* -a test/detector_tests.py
     - coverage run --include=tupak/* -a test/prior_tests.py
-    - coverage run --include=tupak/* -a test/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

docs/likelihood.txt

@@ -9,9 +9,169 @@ 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 default likelihood we use in the examples is `GravitationalWaveTransient`:
+The simplest likelihood
+-----------------------
+
+To start with let's consider perhaps the simplest likelihood we could write
+down, namely a Gaussian likelihood for a set of data :math:`\vec{x}=[x_1, x_2,
+\ldots, x_N]`. The likelihood for a single data point, given the mean
+:math:`\mu` and standard-deviation :math:`\sigma` is given by
+
+.. math::
+
+   \mathcal{L}(x_i| \mu, \sigma) =
+   \frac{1}{\sqrt{2\pi\sigma^2}}\mathrm{exp}\left(
+   \frac{-(x_i - \mu)^2}{2\sigma^2}\right)
+
+Then, the likelihood for all :math:`N` data points is
+
+.. math::
+
+   \mathcal{L}(\vec{x}| \mu, \sigma) = \prod_{i=1}^N
+   \mathcal{L}(x_i| \mu, \sigma)
+
+In practise, we implement the log-likelihood to avoid numerical overflow
+errors. To code this up in :code:`tupak`, we would write a class like this::
+
+   class GaussianLikelihood(tupak.Likelihood):
+       def __init__(self, data):
+           """
+           A very simple Gaussian likelihood
+
+           Parameters
+           ----------
+           data: array_like
+               The data to analyse
+           """
+           self.data = data
+           self.N = len(data)
+           self.parameters = {'mu': None, 'sigma': None}
+
+       def log_likelihood(self):
+           mu = self.parameters['mu']
+           sigma = self.parameters['sigma']
+           res = self.data - mu
+           return -0.5 * (np.sum((res / sigma)**2)
+                          + self.N*np.log(2*np.pi*sigma**2))
+
+
+This demonstrates the two required features of a :code:`tupak`
+:code:`Likelihood` object:
+
+#. It has a `parameters` attribute (a dictionary with
+   keys for all the parameters, in this case, initialised to :code:`None`)
+
+#. It has a :code:`log_likelihood` method which, when called returns the log
+   likelihood for all the data.
+
+You can find an example that uses this likelihood `here <https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/gaussian_example.py>`_.
+
+.. tip::
+
+   Note that the example above subclasses the :code:`tupak.Likelihood` base
+   class, this simply provides a few in built functions. We recommend you also
+   do this when writing your own likelihood.
+
+
+General likelihood for fitting a function :math:`y(x)` to some data with known noise
+------------------------------------------------------------------------------------
+
+The previous example was rather simplistic, Let's now consider that we have some
+dependent data :math:`\vec{y}=y_1, y_2, \ldots y_N$` measured at
+:math:`\vec{x}=x_1, x_2, \ldots, x_N`. We believe that the data is generated
+by additive Gaussian noise with a known variance :math:`sigma^2` and a function
+:math:`y(x; \theta)` where :math:`\theta`  are some unknown parameters; that is
+
+.. math::
+
+   y_i = y(x_i; \theta) + n_i
+
+where :math:`n_i` is drawn from a normal distribution with zero mean and
+standard deviation :math:`sigma`. As such, :math:`y_i - y(x_i; \theta)`
+itself will have a likelihood
+
+.. math::
+
+   \mathcal{L}(y_i; x_i, \theta) =
+   \frac{1}{\sqrt{2\pi\sigma^{2}}}
+   \mathrm{exp}\left(\frac{-(y_i - y(x_i; \theta))^2}{2\sigma^2}\right)
+
+
+As with the previous case, the likelihood for all the data is the produce over
+the likelihood for each data point.
+
+In :code:`tupak`, we can code this up as a likelihood in the following way::
+
+   class GaussianLikelihood(tupak.Likelihood):
+       def __init__(self, x, y, sigma, function):
+           """
+           A general Gaussian likelihood - the parameters are inferred from the
+           arguments of function
+
+           Parameters
+           ----------
+           x, y: array_like
+               The data to analyse
+           sigma: float
+               The standard deviation of the noise
+           function:
+               The python function to fit to the data. Note, this must take the
+               dependent variable as its first argument. The other arguments are
+               will require a prior and will be sampled over (unless a fixed
+               value is given).
+           """
+           self.x = x
+           self.y = y
+           self.sigma = sigma
+           self.N = len(x)
+           self.function = function
+
+           # These lines of code infer the parameters from the provided function
+           parameters = inspect.getargspec(function).args
+           parameters.pop(0)
+           self.parameters = dict.fromkeys(parameters)
+
+       def log_likelihood(self):
+           res = self.y - self.function(self.x, **self.parameters)
+           return -0.5 * (np.sum((res / self.sigma)**2)
+                          + self.N*np.log(2*np.pi*self.sigma**2))
+
+       def noise_log_likelihood(self):
+           return -0.5 * (np.sum((self.y / self.sigma)**2)
+                          + self.N*np.log(2*np.pi*self.sigma**2))
+
+
+This likelihood can be given any python function, the data (in the form of
+:code:`x` and :code:`y`) and the standard deviation of the noise. The
+parameters are inferred from the arguments to the :code:`function` argument,
+for example if, when instatiating the likelihood you passed in a the following
+function::
+
+   def f(x, a, b):
+       return x**2 + b
+
+Then you would also need to provide priors for :code:`a` and :code:`b`. For
+this likelihood, the first argument to :code:`function` is always assumed to
+be the dependent variable.
+
+.. note::
+    Here we have explicitly defined the :code:`noise_log_likelihood` method
+    as the case when there is no signal (i.e., :math:`y(x; \theta)=0`).
+
+You can see an example of this likelihood in the `linear regression example
+<https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/linear_regression.py>`_.
+
+
+Likelihood for transient gravitational waves
+--------------------------------------------
+
+In the examples above, we show how to write your own likelihood. However, for
+routine gravitational wave data analysis of transient events, we have in built
+likelihoods.  The default likelihood we use in the examples is
+`GravitationalWaveTransient`:
 
 .. autoclass:: tupak.likelihood.GravitationalWaveTransient
 
@@ -19,6 +179,9 @@ We also provide a simpler likelihood
 
 .. autoclass:: tupak.likelihood.BasicGravitationalWaveTransient
 
+Empty likelihood for subclassing
+--------------------------------
+
 We provide an empty parent class which can be subclassed for alternative use
 cases
 

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/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
 
@@ -56,8 +51,10 @@ fig.savefig('{}/{}_data.png'.format(outdir, label))
 
 
 class GaussianLikelihood(tupak.Likelihood):
-    def __init__(self, x, y, sigma, waveform_generator):
+    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):
             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):
                        + 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.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
@@ -109,5 +103,6 @@ priors['c'] = tupak.prior.Uniform(-2, 2, 'c')
 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)

examples/other_examples/linear_regression_unknown_noise.py

+#!/bin/python
+"""
+An example of how to use tupak to perform paramater estimation for
+non-gravitational wave data. In this case, fitting a linear function to
+data with background Gaussian noise with unknown variance.
+
+"""
+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_unknown_noise'
+outdir = 'outdir'
+
+
+# First, we define our "signal model", in this case a simple linear function
+def model(time, m, c):
+    return time * m + c
+
+
+# New we define the injection parameters which we make simulated data with
+injection_parameters = dict(m=0.5, c=0.2)
+
+# For this example, we'll inject standard Gaussian noise
+sigma = 1
+
+# These lines of code generate the fake data. Note the ** just unpacks the
+# contents of the injection_paramsters when calling the model function.
+sampling_frequency = 10
+time_duration = 10
+time = np.arange(0, time_duration, 1/sampling_frequency)
+N = len(time)
+data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)
+
+# We quickly plot the data to check it looks sensible
+fig, ax = plt.subplots()
+ax.plot(time, data, 'o', label='data')
+ax.plot(time, model(time, **injection_parameters), '--r', label='signal')
+ax.set_xlabel('time')
+ax.set_ylabel('y')
+ax.legend()
+fig.savefig('{}/{}_data.png'.format(outdir, label))
+
+
+# Parameter estimation: we now define a Gaussian Likelihood class relevant for
+# our model.
+
+class GaussianLikelihood(tupak.Likelihood):
+    def __init__(self, x, y, function, sigma=None):
+        """
+        A general Gaussian likelihood for known or unknown noise - the model
+        parameters are inferred from the arguments of function
+
+        Parameters
+        ----------
+        x, y: array_like
+            The data to analyse
+        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).
+        sigma: None, float, array_like
+            If None, the standard deviation of the noise is unknown and will be
+            estimated (note: this requires a prior to be given for sigma). If
+            not None, this defined the standard-deviation of the data points.
+            This can either be a single float, or an array with length equal
+            to that for `x` and `y`.
+        """
+        self.x = x
+        self.y = y
+        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)
+        self.function_keys = self.parameters.keys()
+        self.parameters['sigma'] = None
+
+    def log_likelihood(self):
+        model_parameters = {k: self.parameters[k] for k in self.function_keys}
+        res = self.y - self.function(self.x, **model_parameters)
+        sigma = self.parameters['sigma']
+        return -0.5 * (np.sum((res / sigma)**2)
+                       + self.N*np.log(2*np.pi*sigma**2))
+
+    def noise_log_likelihood(self):
+        return np.nan
+        sigma = self.parameters['sigma']
+        return -0.5 * (np.sum((self.y / sigma)**2)
+                       + self.N*np.log(2*np.pi*sigma**2))
+
+
+# Now lets instantiate a version of our GaussianLikelihood, giving it
+# the time, data and signal model
+likelihood = GaussianLikelihood(time, data, model)
+
+# From hereon, the syntax is exactly equivalent to other tupak examples
+# We make a prior
+priors = {}
+priors['m'] = tupak.prior.Uniform(0, 5, 'm')
+priors['c'] = tupak.prior.Uniform(-2, 2, 'c')
+priors['sigma'] = tupak.prior.Uniform(0, 10, 'sigma')
+
+# And run sampler
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
+    walks=10, injection_parameters=injection_parameters, outdir=outdir,
+    label=label)
+result.plot_corner()
+print(result)

examples/other_examples/sine_gaussian_example.py

+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation on a sine gaussian injected signal.
+
+"""
+from __future__ import division, print_function
+import tupak
+import numpy as np
+
+# Set the duration and sampling frequency of the data segment that we're going to inject the signal into
+time_duration = 4.
+sampling_frequency = 2048.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'sine_gaussian'
+tupak.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility.  This is optional!
+np.random.seed(170801)
+
+# We are going to inject a sine gaussian waveform.  We first establish a dictionary of parameters that
+# includes all of the different waveform parameters
+injection_parameters = dict(hrss = 1e-22, Q = 5.0, frequency = 200.0, ra = 1.375, dec = -1.2108, 
+                             geocent_time = 1126259642.413, psi= 2.659)
+
+# Create the waveform_generator using a sine Gaussian source function
+waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
+                                                                sampling_frequency=sampling_frequency,
+                                                                frequency_domain_source_model=tupak.source.sinegaussian,
+                                                                parameters=injection_parameters)
+hf_signal = waveform_generator.frequency_domain_strain()
+
+# Set up interferometers.  In this case we'll use three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1),
+# and Virgo (V1)).  These default to their design sensitivity
+IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
+    name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, time_duration=time_duration,
+    sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']]
+
+# Set up prior, which is a dictionary
+priors = dict()
+# By default we will sample all terms in the signal models.  However, this will take a long time for the calculation,
+# so for this example we will set almost all of the priors to be equall to their injected values.  This implies the
+# prior is a delta function at the true, injected value.  In reality, the sampler implementation is smart enough to
+# not sample any parameter that has a delta-function prior.
+for key in ['psi', 'ra', 'dec', 'geocent_time']:
+    priors[key] = injection_parameters[key]
+
+# The above list does *not* include frequency and Q, which means those are the parameters
+# that will be included in the sampler.  If we do nothing, then the default priors get used.
+#priors['Q'] = tupak.prior.create_default_prior(name='Q')
+#priors['frequency'] = tupak.prior.create_default_prior(name='frequency')
+priors['Q'] = tupak.prior.Uniform(2, 50, 'Q')
+priors['frequency'] = tupak.prior.Uniform(30, 1000, 'frequency')
+priors['hrss'] = tupak.prior.Uniform(1e-23, 1e-21, 'hrss')
+
+# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
+likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
+
+# Run sampler.  In this case we're going to use the `dynesty` sampler
+result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
+                                   injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+print(result)
+
+
+
+
+
+
+
+
+
+
+

examples/supernova_example/supernova_example.py

+#!/bin/python
+"""
+
+Tutorial to demonstrate running parameter estimation/model selection on an NR
+supernova injected signal.  Signal model is made by applying PCA to a set of
+supernova waveforms. The first few PCs are then linearly combined with a scale
+factor. (See https://arxiv.org/pdf/1202.3256.pdf)
+
+"""
+from __future__ import division, print_function
+import tupak
+import numpy as np
+
+# Set the duration and sampling frequency of the data segment that we're going to inject the signal into
+time_duration = 3.
+sampling_frequency = 4096.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'supernova'
+tupak.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility.  This is optional!
+np.random.seed(170801)
+
+# We are going to inject a supernova waveform.  We first establish a dictionary
+# of parameters that includes all of the different waveform parameters. It will
+# read in a signal to inject from a txt file.
+
+injection_parameters = dict(file_path='MuellerL15_example_inj.txt',
+                            luminosity_distance=10.0, ra=1.375,
+                            dec=-1.2108, geocent_time=1126259642.413,
+                            psi=2.659)
+
+# Create the waveform_generator using a supernova source function
+waveform_generator = tupak.waveform_generator.WaveformGenerator(
+    time_duration=time_duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.source.supernova,
+    parameters=injection_parameters)
+hf_signal = waveform_generator.frequency_domain_strain()
+
+# Set up interferometers.  In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).  These default to
+# their design sensitivity
+IFOs = [tupak.detector.get_interferometer_with_fake_noise_and_injection(
+    name, injection_polarizations=hf_signal,
+    injection_parameters=injection_parameters, time_duration=time_duration,
+    sampling_frequency=sampling_frequency, outdir=outdir)
+    for name in ['H1', 'L1', 'V1']]
+
+# read in from a file the PCs used to create the signal model.
+realPCs = np.loadtxt('SupernovaRealPCs.txt')
+imagPCs = np.loadtxt('SupernovaImagPCs.txt')
+
+# Now we make another waveform_generator because the signal model is
+# not the same as the injection in this case.
+simulation_parameters = dict(
+    realPCs=realPCs, imagPCs=imagPCs)
+
+search_waveform_generator = tupak.waveform_generator.WaveformGenerator(
+    time_duration=time_duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.source.supernova_pca_model,
+    parameters=simulation_parameters)
+
+# Set up prior
+priors = dict()
+for key in ['psi', 'geocent_time']:
+    priors[key] = injection_parameters[key]
+priors['luminosity_distance'] = tupak.prior.Uniform(
+    2, 20, 'luminosity_distance')
+priors['pc_coeff1'] = tupak.prior.Uniform(-1, 1, 'pc_coeff1')
+priors['pc_coeff2'] = tupak.prior.Uniform(-1, 1, 'pc_coeff2')
+priors['pc_coeff3'] = tupak.prior.Uniform(-1, 1, 'pc_coeff3')
+priors['pc_coeff4'] = tupak.prior.Uniform(-1, 1, 'pc_coeff4')
+priors['pc_coeff5'] = tupak.prior.Uniform(-1, 1, 'pc_coeff5')
+priors['ra'] = tupak.prior.create_default_prior(name='ra')
+priors['dec'] = tupak.prior.create_default_prior(name='dec')
+priors['geocent_time'] = tupak.prior.Uniform(
+    injection_parameters['geocent_time'] - 1,
+    injection_parameters['geocent_time'] + 1,
+    'geocent_time')
+
+# Initialise the likelihood by passing in the interferometer data (IFOs) and
+# the waveoform generator
+likelihood = tupak.GravitationalWaveTransient(
+    interferometers=IFOs, waveform_generator=search_waveform_generator)
+
+# Run sampler.
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
+    outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+print(result)
+
+
+
+
+
+
+
+
+
+
+

test/detector_tests.py

@@ -381,7 +381,7 @@ class TestDetector(unittest.TestCase):
             self.assertEqual(self.ifo.time_delay_from_geocenter(1, 2, 3), 1)
 
     def test_vertex_position_geocentric(self):
-        with mock.patch('tupak.utils.vertex_position_geocentric') as m:
+        with mock.patch('tupak.utils.get_vertex_position_geocentric') as m:
             m.return_value = 1
             self.assertEqual(self.ifo.vertex_position_geocentric(), 1)
 
@@ -389,3 +389,7 @@ class TestDetector(unittest.TestCase):
         self.ifo.data = np.array([2])
         self.ifo.power_spectral_density.power_spectral_density_interpolated = MagicMock(return_value=np.array([1]))
         self.assertTrue(np.array_equal(self.ifo.whitened_data, np.array([2])))
+
+
+if __name__ == '__main__':
+    unittest.main()

test/example_tests.py

@@ -6,7 +6,10 @@ import logging
 # Required to run the tests
 from past.builtins import execfile
 from context import tupak
+
+# Imported to ensure the examples run
 import numpy as np
+import inspect
 
 tupak.utils.command_line_args.test = True
 

test/prior_tests.py

@@ -189,7 +189,6 @@ class TestPriorClasses(unittest.TestCase):
                 domain = np.linspace(-1e2, 1e2, 1000)
             else:
                 domain = np.linspace(prior.minimum, prior.maximum, 1000)
-            print(prior.minimum, prior.maximum)
             self.assertAlmostEqual(np.trapz(prior.prob(domain), domain), 1, 3)
 
 

test/sampler_tests.py

@@ -126,5 +126,9 @@ class TestSampler(unittest.TestCase):
 
     def test_base_run_sampler(self):
         sampler_copy = copy.copy(self.sampler)
-        self.sampler.run_sampler()
-        self.assertDictEqual(sampler_copy.__dict__, self.sampler.__dict__)
\ No newline at end of file
+        self.sampler._run_external_sampler()
+        self.assertDictEqual(sampler_copy.__dict__, self.sampler.__dict__)
+
+
+if __name__ == '__main__':
+    unittest.main()

test/tests.py

@@ -8,14 +8,14 @@ from past.builtins import execfile
 
 
 class Test(unittest.TestCase):
-    outdir = 'TestDir'
+    outdir = './outdir'
     dir_path = os.path.dirname(os.path.realpath(__file__))
     dir_path = os.path.abspath(os.path.join(dir_path, os.path.pardir))
 
     # Run a script to produce standard data
     msd = {}  # A dictionary of variables saved in make_standard_data.py
     execfile(dir_path + '/test/make_standard_data.py', msd)
-
+    '''
     @classmethod
     def setUpClass(self):
         if os.path.isdir(self.outdir):
@@ -33,7 +33,7 @@ class Test(unittest.TestCase):
             except OSError:
                 logging.warning(
                     "{} not removed prior to tests".format(self.outdir))
-
+    '''
     def test_make_standard_data(self):
         " Load in the saved standard data and compare with new data "
 
@@ -57,6 +57,7 @@ class Test(unittest.TestCase):
         priors['luminosity_distance'] = tupak.prior.Uniform(
             name='luminosity_distance', minimum=dL - 10, maximum=dL + 10)
 
+
         result = tupak.sampler.run_sampler(
             likelihood, priors, sampler='nestle', verbose=False, npoints=100)
         self.assertAlmostEqual(np.mean(result.samples), dL,
@@ -65,5 +66,3 @@ class Test(unittest.TestCase):
 
 if __name__ == '__main__':
     unittest.main()
-
-

test/waveform_generator_tests.py

@@ -8,6 +8,7 @@ def gaussian_frequency_domain_strain(frequency_array, amplitude, mu, sigma, ra,
           'cross': amplitude * np.exp(-(mu - frequency_array) ** 2 / sigma ** 2 / 2)}
     return ht
 
+
 def gaussian_frequency_domain_strain_2(frequency_array, a, m, s, ra, dec, geocent_time, psi):
     ht = {'plus': a * np.exp(-(m - frequency_array) ** 2 / s ** 2 / 2),
           'cross': a * np.exp(-(m - frequency_array) ** 2 / s ** 2 / 2)}
@@ -69,20 +70,6 @@ class TestParameterSetter(unittest.TestCase):
         for key in self.simulation_parameters:
             self.assertEqual(self.waveform_generator.parameters[key], self.simulation_parameters[key])
 
-    def test_parameter_setter_with_unexpected_keys(self):
-        self.waveform_generator.parameters['foo'] = 1337
-
-        def parameter_setter(wg_params, sim_params):
-            wg_params = sim_params
-        self.assertRaises(KeyError, parameter_setter(self.waveform_generator.parameters, self.simulation_parameters))
-
-    def test_parameter_setter_raises_type_error(self):
-        a = 4
-
-        def parameter_setter(wg_params, sim_params):
-            wg_params = sim_params
-        self.assertRaises(TypeError, parameter_setter(self.waveform_generator.parameters, a))
-
     def test_parameter_setter_none_handling(self):
         self.waveform_generator.parameters = None
         self.assertItemsEqual(self.waveform_generator.parameters.keys(), self.simulation_parameters.keys())
@@ -109,4 +96,4 @@ class TestSourceModelSetter(unittest.TestCase):
 
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    unittest.main()

tupak/detector.py

@@ -323,8 +323,8 @@ class Interferometer(object):
         return self.power_spectral_density.power_spectral_density_interpolated(self.frequency_array)
 
     def set_data(self, sampling_frequency, duration, epoch=0,
-                 from_power_spectral_density=False,
-                 frequency_domain_strain=None):
+                 from_power_spectral_density=False, frame_file=None,
+                 frequency_domain_strain=None, channel_name=None, overwrite_psd=True, **kwargs):
         """
         Set the interferometer frequency-domain stain and accompanying PSD values.
 
@@ -341,6 +341,15 @@ class Interferometer(object):
         from_power_spectral_density: bool
             If frequency_domain_strain not given, use IFO's PSD object to
             generate noise
+        frame_file: str
+            File from which to load data.
+        channel_name: str
+            Channel to read from frame.
+        overwrite_psd: bool
+            Whether to overwrite the psd in the interferometer with one calculated
+            from the loaded data, default=True.
+        kwargs: dict
+            Additional arguments for loading data.
         """
 
         self.epoch = epoch
@@ -356,6 +365,14 @@ class Interferometer(object):
             frequency_domain_strain, frequencies = \
                 self.power_spectral_density.get_noise_realisation(
                     sampling_frequency, duration)
+        elif frame_file is not None:
+            logging.info('Reading data from frame, {}.'.format(self.name))
+            strain = tupak.utils.read_frame_file(
+                frame_file, t1=epoch, t2=epoch+duration, channel=channel_name, resample=sampling_frequency)
+            frequency_domain_strain, frequencies = tupak.utils.process_strain_data(strain, **kwargs)
+            if overwrite_psd:
+                self.power_spectral_density = PowerSpectralDensity(
+                    frame_file=frame_file, channel_name=channel_name, epoch=epoch, **kwargs)
         else:
             raise ValueError("No method to set data provided.")
 
@@ -407,16 +424,33 @@ class Interferometer(object):
 
 class PowerSpectralDensity:
 
-    def __init__(self, asd_file=None, psd_file='aLIGO_ZERO_DET_high_P_psd.txt'):
+    def __init__(self, asd_file=None, psd_file='aLIGO_ZERO_DET_high_P_psd.txt', frame_file=None, epoch=0,
+                 psd_duration=1024, psd_offset=16, channel_name=None, filter_freq=1024, alpha=0.25, fft_length=4):
         """
         Instantiate a new PowerSpectralDensity object.
 
         Only one of the asd_file or psd_file needs to be specified.
         If multiple are given, the first will be used.
-        FIXME: Allow reading a frame and then FFT to get PSD, use gwpy?
 
-        :param asd_file: amplitude spectral density, format 'f h_f'
-        :param psd_file: power spectral density, format 'f h_f'
+        Parameters:
+        asd_file: str
+            File containing amplitude spectral density, format 'f h_f'
+        psd_file: str
+            File containing power spectral density, format 'f h_f'
+        frame_file: str
+            Frame file to read data from.
+        epoch: float
+            Beginning of segment to analyse.
+        psd_duration: float
+            Duration of data to generate PSD from.
+        psd_offset: float
+            Offset of data from beginning of analysed segment.
+        channel_name: str
+            Name of channel to use to generate PSD.
+        alpha: float
+            Parameter for Tukey window.
+        fft_length: float
+            Number of seconds in a single fft.
         """
 
         self.frequencies = []
@@ -435,6 +469,24 @@ class PowerSpectralDensity:
                 logging.warning("The minimum of the provided curve is {:.2e}.".format(
                     min(self.amplitude_spectral_density)))
                 logging.warning("You may have intended to provide this as a power spectral density.")
+        elif frame_file is not None:
+            strain = tupak.utils.read_frame_file(frame_file, t1=epoch - psd_duration - psd_offset,
+                                                 t2=epoch - psd_duration, channel=channel_name)
+            sampling_frequency = int(strain.sample_rate.value)
+
+            # Low pass filter
+            bp = filter_design.lowpass(filter_freq, strain.sample_rate)
+            strain = strain.filter(bp, filtfilt=True)
+            strain = strain.crop(*strain.span.contract(1))
+
+            # Create and save PSDs
+            NFFT = int(sampling_frequency * fft_length)
+            window = signal.windows.tukey(NFFT, alpha=alpha)
+            psd = strain.psd(fftlength=fft_length, window=window)
+            self.frequencies = psd.frequencies
+            self.power_spectral_density = psd.value
+            self.amplitude_spectral_density = self.power_spectral_density**0.5
+            self.interpolate_power_spectral_density()
         else:
             self.power_spectral_density_file = psd_file
             self.import_power_spectral_density()

tupak/likelihood.py

 from __future__ import division, print_function
+
 import numpy as np
 
 try:
@@ -7,7 +8,6 @@ except ImportError:
     from scipy.misc import logsumexp
 from scipy.special import i0e
 from scipy.interpolate import interp1d
-from scipy.integrate import quad as gaussian_quadrature
 from scipy.interpolate import interp2d
 import tupak
 import logging
@@ -60,6 +60,7 @@ class GravitationalWaveTransient(Likelihood):
         some model parameters
 
     """
+
     def __init__(self, interferometers, waveform_generator, time_marginalization=False, distance_marginalization=False,
                  phase_marginalization=False, prior=None):
 
@@ -77,7 +78,7 @@ class GravitationalWaveTransient(Likelihood):
             self.delta_distance = 0
             self.distance_prior_array = np.array([])
             self.setup_distance_marginalization()
-            prior['luminosity_distance'] = 1 # this means the prior is a delta function fixed at the RHS value
+            prior['luminosity_distance'] = 1  # this means the prior is a delta function fixed at the RHS value
 
         if self.phase_marginalization:
             self.bessel_function_interped = None
@@ -111,8 +112,7 @@ class GravitationalWaveTransient(Likelihood):
 
         matched_filter_snr_squared = 0
         optimal_snr_squared = 0
-        matched_filter_snr_squared_tc_array = None
-
+        matched_filter_snr_squared_tc_array = np.zeros(self.interferometers[0].data[0:-1].shape, dtype=np.complex128)
         for interferometer in self.interferometers:
             signal_ifo = interferometer.get_detector_response(waveform_polarizations,
                                                               self.waveform_generator.parameters)
@@ -121,69 +121,52 @@ class GravitationalWaveTransient(Likelihood):
 
             optimal_snr_squared += tupak.utils.optimal_snr_squared(
                 signal_ifo, interferometer, self.waveform_generator.time_duration)
-
             if self.time_marginalization:
-
                 interferometer.time_marginalization = self.time_marginalization
-
-                if matched_filter_snr_squared_tc_array is None:
-
-                    matched_filter_snr_squared_tc_array = 4. * (1./interferometer.duration) \
-                                                            * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
-                                                            / interferometer.power_spectral_density_array[0:-1]) * len(interferometer.data[0:-1])
-                else:
-
-                    matched_filter_snr_squared_tc_array += 4. * (1./interferometer.duration) \
-                                                        * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
-                                                        / interferometer.power_spectral_density_array[0:-1] ) * len(interferometer.data[0:-1])
+                matched_filter_snr_squared_tc_array += 4. * (1. / interferometer.duration) * np.fft.ifft(
+                    signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1]
+                    / interferometer.power_spectral_density_array[0:-1]) * len(interferometer.data[0:-1])
 
         if self.time_marginalization:
 
-            delta_tc = 1./self.waveform_generator.sampling_frequency
-            tc_log_norm = np.log(interferometer.duration*delta_tc)
+            delta_tc = 1. / self.waveform_generator.sampling_frequency
+            tc_log_norm = np.log(self.interferometers[0].duration * delta_tc)
 
             if self.distance_marginalization:
 
-                rho_opt_ref = optimal_snr_squared.real * \
-                             self.waveform_generator.parameters['luminosity_distance'] ** 2 \
-                            / self.ref_dist**2.
-                rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \
-                            self.waveform_generator.parameters['luminosity_distance'] \
-                            / self.ref_dist
-
+                rho_mf_ref_tc_array, rho_opt_ref = self.__setup_rho(matched_filter_snr_squared_tc_array,
+                                                                    optimal_snr_squared)
 
                 if self.phase_marginalization:
 
-                    phase_margd_rho_mf_tc_array = self.bessel_function_interped(abs(rho_mf_ref_tc_array))
+                    phase_marged_rho_mf_tc_array = self.bessel_function_interped(abs(rho_mf_ref_tc_array))
 
-                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(phase_margd_rho_mf_tc_array, rho_opt_ref)
+                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(phase_marged_rho_mf_tc_array,
+                                                                                      rho_opt_ref)
                     log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - tc_log_norm
 
                 else:
 
-                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array.real, rho_opt_ref)
+                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array.real,
+                                                                                      rho_opt_ref)
 
-                    log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) 
+                    log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc)
 
             elif self.phase_marginalization:
-
-                log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
-                                                    - optimal_snr_squared/2. - tc_log_norm
+                log_l = logsumexp(self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
+                        - optimal_snr_squared / 2. - tc_log_norm
 
             else:
 
-                log_l = logsumexp(matched_filter_snr_squared_tc_array.real, axis=0, b=delta_tc) - optimal_snr_squared / 2. - tc_log_norm
+                log_l = logsumexp(matched_filter_snr_squared_tc_array.real, axis=0,
+                                  b=delta_tc) - optimal_snr_squared / 2. - tc_log_norm
 
         elif self.distance_marginalization:
 
-            rho_opt_ref = optimal_snr_squared.real * \
-                         self.waveform_generator.parameters['luminosity_distance'] ** 2 \
-                        / self.ref_dist**2.
-            rho_mf_ref = matched_filter_snr_squared.real * \
-                        self.waveform_generator.parameters['luminosity_distance'] \
-                        / self.ref_dist
+            rho_mf_ref, rho_opt_ref = self.__setup_rho(matched_filter_snr_squared.real,
+                                                       optimal_snr_squared)
 
-            log_l=self.interp_dist_margd_loglikelihood(rho_mf_ref, rho_opt_ref)[0]
+            log_l = self.interp_dist_margd_loglikelihood(rho_mf_ref, rho_opt_ref)[0]
 
         elif self.phase_marginalization:
             matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared))
@@ -194,11 +177,20 @@ class GravitationalWaveTransient(Likelihood):
 
         return log_l.real
 
+    def __setup_rho(self, matched_filter_snr_squared, optimal_snr_squared):
+        rho_opt_ref = optimal_snr_squared.real * \
+                      self.waveform_generator.parameters['luminosity_distance'] ** 2 \
+                      / self.ref_dist ** 2.
+        rho_mf_ref = matched_filter_snr_squared * \
+                     self.waveform_generator.parameters['luminosity_distance'] \
+                     / self.ref_dist
+        return rho_mf_ref, rho_opt_ref
+
     def log_likelihood(self):
         return self.log_likelihood_ratio() + self.noise_log_likelihood()
 
     def setup_distance_marginalization(self):
-        if 'fluminosity_distance' not in self.prior.keys():
+        if 'luminosity_distance' not in self.prior.keys():
             logging.info('No prior provided for distance, using default prior.')
             self.prior['luminosity_distance'] = tupak.prior.create_default_prior('luminosity_distance')
         self.distance_array = np.linspace(self.prior['luminosity_distance'].minimum,
@@ -208,36 +200,35 @@ class GravitationalWaveTransient(Likelihood):
                                               for distance in self.distance_array])
 
         ### Make the lookup table ###
-        self.ref_dist = 1000 # 1000 Mpc
-        self.dist_margd_loglikelihood_array = np.zeros((400,800))
-
-        self.rho_opt_ref_array = np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[0]) # optimal filter snr at fiducial distance of ref_dist Mpc
-        self.rho_mf_ref_array = np.hstack( (-np.logspace(2,-3,self.dist_margd_loglikelihood_array.shape[1]/2), \
-                                np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[1]/2)) ) # matched filter snr at fiducial distance of ref_dist Mpc
+        self.ref_dist = 1000  # 1000 Mpc
+        self.dist_margd_loglikelihood_array = np.zeros((400, 800))
 
+        self.rho_opt_ref_array = np.logspace(-3, 4, self.dist_margd_loglikelihood_array.shape[
+            0])  # optimal filter snr at fiducial distance of ref_dist Mpc
+        self.rho_mf_ref_array = np.hstack((-np.logspace(2, -3, self.dist_margd_loglikelihood_array.shape[1] / 2), \
+                                           np.logspace(-3, 4, self.dist_margd_loglikelihood_array.shape[
+                                               1] / 2)))  # matched filter snr at fiducial distance of ref_dist Mpc
 
         for ii, rho_opt_ref in enumerate(self.rho_opt_ref_array):
 
             for jj, rho_mf_ref in enumerate(self.rho_mf_ref_array):
-
-                optimal_snr_squared_array = rho_opt_ref * self.ref_dist**2. \
+                optimal_snr_squared_array = rho_opt_ref * self.ref_dist ** 2. \
                                             / self.distance_array ** 2
 
                 matched_filter_snr_squared_array = rho_mf_ref * self.ref_dist \
-                                                / self.distance_array
+                                                   / self.distance_array
 
                 self.dist_margd_loglikelihood_array[ii][jj] = \
-                                                    logsumexp(matched_filter_snr_squared_array - \
-                                                     optimal_snr_squared_array / 2, \
-                                                     b= self.distance_prior_array * self.delta_distance)
+                    logsumexp(matched_filter_snr_squared_array - \
+                              optimal_snr_squared_array / 2, \
+                              b=self.distance_prior_array * self.delta_distance)
 
-        log_norm = logsumexp(0. / self.distance_array  - 0./self.distance_array**2.,\
-                            b=self.distance_prior_array* self.delta_distance)
+        log_norm = logsumexp(0. / self.distance_array - 0. / self.distance_array ** 2., \
+                             b=self.distance_prior_array * self.delta_distance)
         self.dist_margd_loglikelihood_array -= log_norm
 
-
-        self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array, self.dist_margd_loglikelihood_array)
-
+        self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array,
+                                                        self.dist_margd_loglikelihood_array)
 
     def setup_phase_marginalization(self):
         if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
@@ -249,7 +240,6 @@ class GravitationalWaveTransient(Likelihood):
                                                  bounds_error=False, fill_value=-np.inf)
 
 
-
 class BasicGravitationalWaveTransient(Likelihood):
     """ A basic gravitational wave transient likelihood
 
@@ -272,6 +262,7 @@ class BasicGravitationalWaveTransient(Likelihood):
         some model parameters
 
     """
+
     def __init__(self, interferometers, waveform_generator):
         Likelihood.__init__(self, waveform_generator.parameters)
         self.interferometers = interferometers

tupak/prior.py

@@ -430,7 +430,7 @@ class FromFile(Interped):
             logging.warning(r"x\tp(x)")
 
     def __repr__(self, subclass_keys=list(), subclass_names=list()):
-        return Prior._subclass_repr_helper(self, subclass_args=['xx', 'yy', 'id'])
+        return Prior._subclass_repr_helper(self, subclass_args=['id'])
 
 
 class UniformComovingVolume(FromFile):