Merge remote-tracking branch 'origin/master' into Simplify_wg_redundant_code

.gitlab-ci.yml

@@ -42,6 +42,7 @@ python-3:
     - coverage --version
     - coverage erase
     - coverage run --debug=trace --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/conversion_tests.py
+    - coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/calibration_tests.py
     - coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/detector_tests.py
     - coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/utils_tests.py
     - coverage run --include=/opt/conda/lib/python3.6/site-packages/tupak* -a test/prior_tests.py

examples/injection_examples/calibration_example.py

+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation with calibration
+uncertainties included.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+import tupak
+
+# Set the duration and sampling frequency of the data segment
+# that we're going to create and inject the signal into.
+
+duration = 4.
+sampling_frequency = 2048.
+
+# Specify the output directory and the name of the simulation.
+outdir = 'outdir'
+label = 'calibration'
+tupak.core.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, ra=1.375, dec=-1.2108)
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+                          reference_frequency=50.)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+waveform_generator = tupak.gw.WaveformGenerator(
+    duration=duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
+    parameters=injection_parameters,waveform_arguments=waveform_arguments)
+
+# 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.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+for ifo in ifos:
+    injection_parameters.update({
+        'recalib_{}_amplitude_{}'.format(ifo.name, ii): 0.1 for ii in range(5)})
+    injection_parameters.update({
+        'recalib_{}_phase_{}'.format(ifo.name, ii): 0.01 for ii in range(5)})
+    ifo.calibration_model = tupak.gw.calibration.CubicSpline(
+        prefix='recalib_{}_'.format(ifo.name),
+        minimum_frequency=ifo.minimum_frequency,
+        maximum_frequency=ifo.maximum_frequency, n_points=5)
+ifos.set_strain_data_from_power_spectral_densities(
+    sampling_frequency=sampling_frequency, duration=duration)
+ifos.inject_signal(parameters=injection_parameters,
+                   waveform_generator=waveform_generator)
+
+# Set up prior, which is a dictionary
+# Here we fix the injected cbc parameters and most of the calibration parameters
+# to the injected values.
+# We allow a subset of the calibration parameters to vary.
+priors = injection_parameters.copy()
+for key in injection_parameters:
+    if 'recalib' in key:
+        priors[key] = injection_parameters[key]
+for name in ['recalib_H1_amplitude_0', 'recalib_H1_amplitude_1']:
+    priors[name] = tupak.prior.Gaussian(
+        mu=0, sigma=0.2, name=name, latex_label='H1 $A_{}$'.format(name[-1]))
+
+# Initialise the likelihood by passing in the interferometer data (IFOs) and
+# the waveform generator
+likelihood = tupak.gw.GravitationalWaveTransient(
+    interferometers=ifos, waveform_generator=waveform_generator)
+
+# Run sampler.  In this case we're going to use the `dynesty` sampler
+result = tupak.run_sampler(
+    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
+    injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# make some plots of the outputs
+result.plot_corner()
+

examples/injection_examples/eccentric_GW150914_inspiral.py

+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation on a reduced parameter space for an injected eccentric binary 
+black hole signal with masses & distnace similar to GW150914.
+
+This uses the same binary parameters that were used to make Figures 1, 2 & 5 in Lower et al. (2018) -> arXiv:1806.05350.
+
+For a more comprehensive look at what goes on in each step, refer to the "basic_tutorial.py" example.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+
+import tupak
+
+import matplotlib.pyplot as plt
+
+duration = 64.
+sampling_frequency = 2048.
+
+outdir = 'outdir'
+label = 'eccentric_GW140914'
+tupak.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility.
+np.random.seed(150914)
+
+injection_parameters = dict(mass_1=35., mass_2=30., eccentricity=0.1, luminosity_distance=440.,
+                            iota=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
+
+waveform_arguments = dict(waveform_approximant='EccentricFD', reference_frequency=10., minimum_frequency=10.)
+
+# Create the waveform_generator using the LAL eccentric black hole no spins source function
+waveform_generator = tupak.gw.WaveformGenerator(
+    duration=duration, sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=tupak.gw.source.lal_eccentric_binary_black_hole_no_spins,
+    parameters=injection_parameters, waveform_arguments=waveform_arguments)
+
+hf_signal = waveform_generator.frequency_domain_strain()
+
+# Setting up three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)) at their design sensitivities.
+# The maximum frequency is set just prior to the point at which the waveform model terminates. This is to avoid any biases 
+# introduced from using a sharply terminating waveform model.
+minimum_frequency = 10.0
+maximum_frequency = 133.0
+
+def get_interferometer(name, injection_polarizations, injection_parameters, duration, sampling_frequency,
+                       minimum_frequency, maximum_frequency, outdir):
+    """
+    Sets up the interferometers & injects the signal into them
+    """
+    start_time = injection_parameters['geocent_time'] + 2 - duration
+    
+    ifo = tupak.gw.detector.get_empty_interferometer(name)
+    if name == 'V1':
+        ifo.power_spectral_density.set_from_power_spectral_density_file('AdV_psd.txt')
+    else:
+        ifo.power_spectral_density.set_from_power_spectral_density_file('aLIGO_ZERO_DET_high_P_psd.txt')
+    
+    ifo.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, 
+            duration=duration, start_time=start_time)
+    
+    injection_polarizations = ifo.inject_signal(parameters=injection_parameters,
+                              injection_polarizations=injection_polarizations,
+                              waveform_generator=waveform_generator)
+
+    signal = ifo.get_detector_response(injection_polarizations, injection_parameters)
+
+    ifo.minimum_frequency = minimum_frequency
+    ifo.maximum_frequency = maximum_frequency
+    
+    ifo.plot_data(signal=signal, outdir=outdir, label=label)
+    ifo.save_data(outdir, label=label)
+    
+    return ifo
+
+# IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection(name, injection_polarizations=hf_signal, 
+#     injection_parameters=injection_parameters, duration=duration,
+#     sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']]
+
+name = ['H1', 'L1', 'V1']
+IFOs = []
+
+for i in range(0,3):
+    IFOs.append(get_interferometer(name[i], injection_polarizations=hf_signal, injection_parameters=injection_parameters,
+                                   duration=duration, sampling_frequency=sampling_frequency, 
+                                   minimum_frequency=minimum_frequency, maximum_frequency=maximum_frequency, 
+                                   outdir=outdir))
+
+# Now we set up the priors on each of the binary parameters.
+priors = dict()
+priors["mass_1"] = tupak.core.prior.Uniform(name='mass_1', minimum=5, maximum=60)
+priors["mass_2"] = tupak.core.prior.Uniform(name='mass_2', minimum=5, maximum=60)
+priors["eccentricity"] = tupak.core.prior.PowerLaw(name='eccentricity', latex_label='$e$', alpha=-1, minimum=1e-4, maximum=0.4)
+priors["luminosity_distance"] =  tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=2e3)
+priors["dec"] =  tupak.core.prior.Cosine(name='dec')
+priors["ra"] =  tupak.core.prior.Uniform(name='ra', minimum=0, maximum=2 * np.pi)
+priors["iota"] =  tupak.core.prior.Sine(name='iota')
+priors["psi"] =  tupak.core.prior.Uniform(name='psi', minimum=0, maximum=2 * np.pi)
+priors["phase"] =  tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+priors["geocent_time"] = tupak.core.prior.Uniform(1180002600.9, 1180002601.1, name='geocent_time')
+
+# Initialising the likelihood function.
+likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
+                                              time_marginalization=False, phase_marginalization=False,
+                                              distance_marginalization=False, prior=priors)
+
+# Now we run sampler (PyMultiNest in our case).
+result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='pymultinest', npoints=1000,
+                           injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# And finally we make some plots of the output posteriors.
+result.plot_corner()

examples/injection_examples/marginalized_likelihood.py

@@ -36,7 +36,7 @@ IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection(
 # Set up prior
 priors = tupak.gw.prior.BBHPriorSet()
 # These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra', 'dec', 'geocent_time']:
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra', 'dec', 'geocent_time']:
     priors[key] = injection_parameters[key]
 
 # Initialise GravitationalWaveTransient

examples/open_data_examples/GW150914.py

@@ -36,13 +36,15 @@ prior = tupak.gw.prior.BBHPriorSet(filename='GW150914.prior')
 # creates the frequency-domain strain. In this instance, we are using the
 # `lal_binary_black_hole model` source model. We also pass other parameters:
 # the waveform approximant and reference frequency.
-waveform_generator = tupak.gw.WaveformGenerator(frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
-                                             waveform_arguments={'waveform_approximant': 'IMRPhenomPv2',
-                                                                 'reference_frequency': 50})
+waveform_generator = tupak.WaveformGenerator(
+    frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole,
+    waveform_arguments={'waveform_approximant': 'IMRPhenomPv2',
+                        'reference_frequency': 50})
 
 # In this step, we define the likelihood. Here we use the standard likelihood
 # function, passing it the data and the waveform generator.
-likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
+likelihood = tupak.gw.likelihood.GravitationalWaveTransient(
+    interferometers, waveform_generator)
 
 # Finally, we run the sampler. This function takes the likelihood and prior
 # along with some options for how to do the sampling and how to save the data

test/calibration_tests.py

+from tupak.gw import calibration
+import unittest
+import numpy as np
+
+
+class TestBaseClass(unittest.TestCase):
+
+    def setUp(self):
+        self.model = calibration.Recalibrate()
+
+    def tearDown(self):
+        del self.model
+
+    def test_calibration_factor(self):
+        frequency_array = np.linspace(20, 1024, 1000)
+        cal_factor = self.model.get_calibration_factor(frequency_array)
+        assert np.alltrue(cal_factor.real == np.ones_like(frequency_array))
+
+
+class TestCubicSpline(unittest.TestCase):
+
+    def setUp(self):
+        self.model = calibration.CubicSpline(
+            prefix='recalib_', minimum_frequency=20, maximum_frequency=1024,
+            n_points=5)
+        self.parameters = {'recalib_{}_{}'.format(param, ii): 0.0
+                           for ii in range(5)
+                           for param in ['amplitude', 'phase']}
+
+    def tearDown(self):
+        del self.model
+        del self.parameters
+
+    def test_calibration_factor(self):
+        frequency_array = np.linspace(20, 1024, 1000)
+        cal_factor = self.model.get_calibration_factor(frequency_array,
+                                                       **self.parameters)
+        assert np.alltrue(cal_factor.real == np.ones_like(frequency_array))
+
+
+class TestCubicSplineRequiresFourNodes(unittest.TestCase):
+
+    def test_cannot_instantiate_with_too_few_nodes(self):
+        for ii in range(6):
+            if ii < 4:
+                with self.assertRaises(ValueError):
+                    calibration.CubicSpline('test', 1, 10, ii)
+            else:
+                calibration.CubicSpline('test', 1, 10, ii)
+
+
+if __name__ == '__main__':
+    unittest.main()

tupak/core/prior.py

@@ -26,7 +26,11 @@ class PriorSet(dict):
         dict.__init__(self)
         if type(dictionary) is dict:
             self.update(dictionary)
-        elif filename:
+        elif type(dictionary) is str:
+            logger.debug('Argument "dictionary" is a string.'
+                         + ' Assuming it is intended as a file name.')
+            self.read_in_file(dictionary)
+        elif type(filename) is str:
             self.read_in_file(filename)
 
     def write_to_file(self, outdir, label):
@@ -41,7 +45,7 @@ class PriorSet(dict):
         """
 
         utils.check_directory_exists_and_if_not_mkdir(outdir)
-        prior_file = os.path.join(outdir, "{}_prior.txt".format(label))
+        prior_file = os.path.join(outdir, "{}.prior".format(label))
         logger.debug("Writing priors to {}".format(prior_file))
         with open(prior_file, "w") as outfile:
             for key in self.keys():
@@ -147,7 +151,7 @@ class PriorSet(dict):
         """
         return self.sample_subset(keys=self.keys(), size=size)
 
-    def sample_subset(self, keys=list(), size=None):
+    def sample_subset(self, keys=iter([]), size=None):
         """Draw samples from the prior set for parameters which are not a DeltaFunction
 
         Parameters
@@ -377,7 +381,7 @@ class Prior(object):
         """
         return self._subclass_repr_helper()
 
-    def _subclass_repr_helper(self, subclass_args=list()):
+    def _subclass_repr_helper(self, subclass_args=iter([])):
         """Helps out subclass _repr__ methods by creating a common template
 
         Parameters

tupak/gw/__init__.py

@@ -6,6 +6,8 @@ import tupak.gw.prior
 import tupak.gw.source
 import tupak.gw.utils
 import tupak.gw.waveform_generator
+from . import calibration
+
 
 from tupak.gw.waveform_generator import WaveformGenerator
 from tupak.gw.likelihood import GravitationalWaveTransient

tupak/gw/calibration.py

+""" Functions for adding calibration factors to waveform templates.
+"""
+
+import numpy as np
+from scipy.interpolate import UnivariateSpline
+
+
+class Recalibrate(object):
+
+    name = 'none'
+
+    def __init__(self, prefix='recalib_'):
+        """
+        Base calibration object. This applies no transformation
+
+        Parameters
+        ----------
+        prefix: str
+            Prefix on parameters relating to the calibration.
+        """
+        self.params = dict()
+        self.prefix = prefix
+
+    def get_calibration_factor(self, frequency_array, **params):
+        """Apply calibration model
+
+        This method should be overwritten by subclasses
+
+        Parameters
+        ----------
+        frequency_array: array-like
+            The frequency values to calculate the calibration factor for.
+        params : dict
+            Dictionary of sampling parameters which includes
+            calibration parameters.
+
+        Returns
+        -------
+        calibration_factor : array-like
+            The factor to multiply the strain by.
+        """
+        self.set_calibration_parameters(**params)
+        return np.ones_like(frequency_array)
+
+    def set_calibration_parameters(self, **params):
+        self.params.update({key[len(self.prefix):]: params[key] for key in params
+                            if self.prefix in key})
+
+
+class CubicSpline(Recalibrate):
+
+    name = 'cubic_spline'
+
+    def __init__(self, prefix, minimum_frequency, maximum_frequency, n_points):
+        """
+        Cubic spline recalibration
+
+        see https://dcc.ligo.org/DocDB/0116/T1400682/001/calnote.pdf
+
+        This assumes the spline points follow
+        np.logspace(np.log(minimum_frequency), np.log(maximum_frequency), n_points)
+
+        Parameters
+        ----------
+        prefix: str
+            Prefix on parameters relating to the calibration.
+        minimum_frequency: float
+            minimum frequency of spline points
+        maximum_frequency: float
+            maximum frequency of spline points
+        n_points: int
+            number of spline points
+        """
+        Recalibrate.__init__(self, prefix=prefix)
+        if n_points < 4:
+            raise ValueError('Cubic spline calibration requires at least 4 spline nodes.')
+        self.n_points = n_points
+        self.spline_points = np.logspace(np.log10(minimum_frequency), np.log10(maximum_frequency), n_points)
+
+    def get_calibration_factor(self, frequency_array, **params):
+        """Apply calibration model
+
+        Parameters
+        ----------
+        frequency_array: array-like
+            The frequency values to calculate the calibration factor for.
+        prefix: str
+            Prefix for calibration parameter names
+        params : dict
+            Dictionary of sampling parameters which includes
+            calibration parameters.
+
+        Returns
+        -------
+        calibration_factor : array-like
+            The factor to multiply the strain by.
+        """
+        self.set_calibration_parameters(**params)
+        amplitude_parameters = [self.params['amplitude_{}'.format(ii)] for ii in range(self.n_points)]
+        amplitude_spline = UnivariateSpline(self.spline_points, amplitude_parameters)
+        delta_amplitude = amplitude_spline(frequency_array)
+
+        phase_parameters = [self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
+        phase_spline = UnivariateSpline(self.spline_points, phase_parameters)
+        delta_phase = phase_spline(frequency_array)
+
+        calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)
+
+        return calibration_factor
+

tupak/gw/detector.py

@@ -4,12 +4,13 @@ import os
 
 import matplotlib.pyplot as plt
 import numpy as np
-import scipy
+from scipy.signal.windows import tukey
 from scipy.interpolate import interp1d
 
 import tupak.gw.utils
 from tupak.core import utils
 from tupak.core.utils import logger
+from .calibration import Recalibrate
 
 try:
     import gwpy
@@ -85,7 +86,7 @@ class InterferometerList(list):
            `waveform_generator.frequency_domain_strain()`. If
            `waveform_generator` is also given, the injection_polarizations will
            be calculated directly and this argument can be ignored.
-        waveform_generator: tupak.gw.waveform_generator
+        waveform_generator: tupak.gw.waveform_generator.WaveformGenerator
             A WaveformGenerator instance using the source model to inject. If
             `injection_polarizations` is given, this will be ignored.
 
@@ -353,7 +354,7 @@ class InterferometerStrainData(object):
             self.roll_off = roll_off
         elif alpha is not None:
             self.roll_off = alpha * self.duration / 2
-        window = scipy.signal.windows.tukey(len(self._time_domain_strain), alpha=self.alpha)
+        window = tukey(len(self._time_domain_strain), alpha=self.alpha)
         self.window_factor = np.mean(window**2)
         return window
 
@@ -741,9 +742,9 @@ class InterferometerStrainData(object):
 class Interferometer(object):
     """Class for the Interferometer """
 
-    def __init__(self, name, power_spectral_density, minimum_frequency,
-                 maximum_frequency, length, latitude, longitude, elevation,
-                 xarm_azimuth, yarm_azimuth, xarm_tilt=0., yarm_tilt=0.):
+    def __init__(self, name, power_spectral_density, minimum_frequency, maximum_frequency,
+                 length, latitude, longitude, elevation, xarm_azimuth, yarm_azimuth,
+                 xarm_tilt=0., yarm_tilt=0., calibration_model=Recalibrate()):
         """
         Instantiate an Interferometer object.
 
@@ -774,6 +775,9 @@ class Interferometer(object):
             ellipsoid earth model in LIGO-T980044-08.
         yarm_tilt: float, optional
             Tilt of the y arm in radians above the horizontal.
+        calibration_model: Recalibration
+            Calibration model, this applies the calibration correction to the
+            template, the default model applies no correction.
         """
         self.__x_updated = False
         self.__y_updated = False
@@ -791,6 +795,7 @@ class Interferometer(object):
         self.yarm_tilt = yarm_tilt
         self.power_spectral_density = power_spectral_density
         self.time_marginalization = False
+        self.calibration_model = calibration_model
         self._strain_data = InterferometerStrainData(
             minimum_frequency=minimum_frequency,
             maximum_frequency=maximum_frequency)
@@ -1187,6 +1192,9 @@ class Interferometer(object):
         signal_ifo = signal_ifo * np.exp(
             -1j * 2 * np.pi * dt * self.frequency_array)
 
+        signal_ifo *= self.calibration_model.get_calibration_factor(
+            self.frequency_array, prefix='recalib_{}_'.format(self.name), **parameters)
+
         return signal_ifo
 
     def inject_signal(self,  parameters=None, injection_polarizations=None,

tupak/gw/prior.py

 import os
-from tupak.core.prior import PriorSet, FromFile, Prior
-
-from tupak.core.utils import logger
+from ..core.prior import PriorSet, FromFile, Prior, DeltaFunction, Gaussian
+from ..core.utils import logger
+import numpy as np
+from scipy.interpolate import UnivariateSpline
 
 
 class UniformComovingVolume(FromFile):
@@ -114,3 +115,168 @@ Prior._default_latex_labels = {
     'psi': '$\psi$',
     'phase': '$\phi$',
     'geocent_time': '$t_c$'}
+
+
+class CalibrationPriorSet(PriorSet):
+
+    def __init__(self, dictionary=None, filename=None):
+        """ Initialises a Prior set for Binary Black holes
+
+        Parameters
+        ----------
+        dictionary: dict, optional
+            See superclass
+        filename: str, optional
+            See superclass
+        """
+        if dictionary is None and filename is not None:
+            filename = os.path.join(os.path.dirname(__file__),
+                                    'prior_files', filename)
+        PriorSet.__init__(self, dictionary=dictionary, filename=filename)
+        self.source = None
+
+    def write_to_file(self, outdir, label):
+        """
+        Write the prior to file. This includes information about the source if
+        possible.
+
+        Parameters
+        ----------
+        outdir: str
+            Output directory.
+        label: str
+            Label for prior.
+        """
+        PriorSet.write_to_file(self, outdir=outdir, label=label)
+        if self.source is not None:
+            prior_file = os.path.join(outdir, "{}.prior".format(label))
+            with open(prior_file, "a") as outfile:
+                outfile.write("# prior source file is {}".format(self.source))
+
+    @staticmethod
+    def from_envelope_file(envelope_file, minimum_frequency,
+                           maximum_frequency, n_nodes, label):
+        """
+        Load in the calibration envelope.
+
+        This is a text file with columns:
+            frequency median-amplitude median-phase -1-sigma-amplitude
+            -1-sigma-phase +1-sigma-amplitude +1-sigma-phase
+
+        Parameters
+        ----------
+        envelope_file: str
+            Name of file to read in.
+        minimum_frequency: float
+            Minimum frequency for the spline.
+        maximum_frequency: float
+            Minimum frequency for the spline.
+        n_nodes: int
+            Number of nodes for the spline.
+        label: str
+            Label for the names of the parameters, e.g., recalib_H1_
+
+        Returns
+        -------
+        prior: PriorSet
+            Priors for the relevant parameters.
+            This includes the frequencies of the nodes which are _not_ sampled.
+        """
+        calibration_data = np.genfromtxt(envelope_file).T
+        frequency_array = calibration_data[0]
+        amplitude_median = 1 - calibration_data[1]
+        phase_median = calibration_data[2]
+        amplitude_sigma = (calibration_data[4] - calibration_data[2]) / 2
+        phase_sigma = (calibration_data[5] - calibration_data[3]) / 2
+
+        nodes = np.logspace(np.log10(minimum_frequency),
+                            np.log10(maximum_frequency), n_nodes)
+
+        amplitude_mean_nodes =\
+            UnivariateSpline(frequency_array, amplitude_median)(nodes)
+        amplitude_sigma_nodes =\
+            UnivariateSpline(frequency_array, amplitude_sigma)(nodes)
+        phase_mean_nodes =\
+            UnivariateSpline(frequency_array, phase_median)(nodes)
+        phase_sigma_nodes =\
+            UnivariateSpline(frequency_array, phase_sigma)(nodes)
+
+        prior = CalibrationPriorSet()
+        for ii in range(n_nodes):
+            name = "recalib_{}_amplitude_{}".format(label, ii)
+            latex_label = "$A^{}_{}$".format(label, ii)
+            prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
+                                   sigma=amplitude_sigma_nodes[ii],
+                                   name=name, latex_label=latex_label)
+        for ii in range(n_nodes):
+            name = "recalib_{}_phase_{}".format(label, ii)
+            latex_label = "$\\phi^{}_{}$".format(label, ii)
+            prior[name] = Gaussian(mu=phase_mean_nodes[ii],
+                                   sigma=phase_sigma_nodes[ii],
+                                   name=name, latex_label=latex_label)
+        for ii in range(n_nodes):
+            name = "recalib_{}_frequency_{}".format(label, ii)
+            latex_label = "$f^{}_{}$".format(label, ii)
+            prior[name] = DeltaFunction(peak=nodes[ii], name=name,
+                                        latex_label=latex_label)
+        prior.source = os.path.abspath(envelope_file)
+        return prior
+
+    @staticmethod
+    def constant_uncertainty_spline(
+            amplitude_sigma, phase_sigma, minimum_frequency, maximum_frequency,
+            n_nodes, label):
+        """
+        Make prior assuming constant in frequency calibration uncertainty.
+
+        This assumes Gaussian fluctuations about 0.
+
+        Parameters
+        ----------
+        amplitude_sigma: float
+            Uncertainty in the amplitude.
+        phase_sigma: float
+            Uncertainty in the phase.
+        minimum_frequency: float
+            Minimum frequency for the spline.
+        maximum_frequency: float
+            Minimum frequency for the spline.
+        n_nodes: int
+            Number of nodes for the spline.
+        label: str
+            Label for the names of the parameters, e.g., recalib_H1_
+
+        Returns
+        -------
+        prior: PriorSet
+            Priors for the relevant parameters.
+            This includes the frequencies of the nodes which are _not_ sampled.
+        """
+        nodes = np.logspace(np.log10(minimum_frequency),
+                            np.log10(maximum_frequency), n_nodes)
+
+        amplitude_mean_nodes = [0] * n_nodes
+        amplitude_sigma_nodes = [amplitude_sigma] * n_nodes
+        phase_mean_nodes = [0] * n_nodes
+        phase_sigma_nodes = [phase_sigma] * n_nodes
+
+        prior = CalibrationPriorSet()
+        for ii in range(n_nodes):
+            name = "recalib_{}_amplitude_{}".format(label, ii)
+            latex_label = "$A^{}_{}$".format(label, ii)
+            prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
+                                   sigma=amplitude_sigma_nodes[ii],
+                                   name=name, latex_label=latex_label)
+        for ii in range(n_nodes):
+            name = "recalib_{}_phase_{}".format(label, ii)
+            latex_label = "$\\phi^{}_{}$".format(label, ii)
+            prior[name] = Gaussian(mu=phase_mean_nodes[ii],
+                                   sigma=phase_sigma_nodes[ii],
+                                   name=name, latex_label=latex_label)
+        for ii in range(n_nodes):
+            name = "recalib_{}_frequency_{}".format(label, ii)
+            latex_label = "$f^{}_{}$".format(label, ii)
+            prior[name] = DeltaFunction(peak=nodes[ii], name=name,
+                                        latex_label=latex_label)
+
+        return prior

tupak/gw/prior_files/GW150914.prior

@@ -14,3 +14,35 @@ iota =  Sine(name='iota')
 psi =  Uniform(name='psi', minimum=0, maximum=2 * np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time')
+# These are the calibration parameters as described in
+# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041015
+# recalib_H1_frequency_0 = 20
+# recalib_H1_frequency_1 = 54
+# recalib_H1_frequency_2 = 143
+# recalib_H1_frequency_3 = 383
+# recalib_H1_frequency_4 = 1024
+# recalib_H1_amplitude_0 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_0), '$\\delta A_{H0}$'
+# recalib_H1_amplitude_1 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_1), '$\\delta A_{H1}$'
+# recalib_H1_amplitude_2 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_2), '$\\delta A_{H2}$'
+# recalib_H1_amplitude_3 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_3), '$\\delta A_{H3}$'
+# recalib_H1_amplitude_4 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_4), '$\\delta A_{H4}$'
+# recalib_H1_phase_0 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_0', '$\\delta \\phi_{H0}$')
+# recalib_H1_phase_1 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_1', '$\\delta \\phi_{H1}$')
+# recalib_H1_phase_2 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_2', '$\\delta \\phi_{H2}$')
+# recalib_H1_phase_3 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_3', '$\\delta \\phi_{H3}$')
+# recalib_H1_phase_4 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_4', '$\\delta \\phi_{H4}$')
+# recalib_L1_frequency_0 = 20
+# recalib_L1_frequency_1 = 54
+# recalib_L1_frequency_2 = 143
+# recalib_L1_frequency_3 = 383
+# recalib_L1_frequency_4 = 1024
+# recalib_L1_amplitude_0 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_0), '$\\delta A_{L0}$'
+# recalib_L1_amplitude_1 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_1), '$\\delta A_{L1}$'
+# recalib_L1_amplitude_2 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_2), '$\\delta A_{L2}$'
+# recalib_L1_amplitude_3 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_3), '$\\delta A_{L3}$'
+# recalib_L1_amplitude_4 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_4), '$\\delta A_{L4}$'
+# recalib_L1_phase_0 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_0', '$\\delta \\phi_{L0}$')
+# recalib_L1_phase_1 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_1', '$\\delta \\phi_{L1}$')
+# recalib_L1_phase_2 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_2', '$\\delta \\phi_{L2}$')
+# recalib_L1_phase_3 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_3', '$\\delta \\phi_{L3}$')
+# recalib_L1_phase_4 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_4', '$\\delta \\phi_{L4}$')

tupak/gw/source.py

@@ -108,6 +108,86 @@ def lal_binary_black_hole(
 
     return {'plus': h_plus, 'cross': h_cross}
 
+def lal_eccentric_binary_black_hole_no_spins(
+        frequency_array, mass_1, mass_2, eccentricity, luminosity_distance, iota, phase, ra, dec, 
+        geocent_time, psi, **kwargs):        
+    """ Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
+
+    Parameters
+    ----------
+    frequency_array: array_like
+        The frequencies at which we want to calculate the strain
+    mass_1: float
+        The mass of the heavier object in solar masses
+    mass_2: float
+        The mass of the lighter object in solar masses
+    eccentricity: float
+        The orbital eccentricity of the system
+    luminosity_distance: float
+        The luminosity distance in megaparsec
+    iota: float
+        Orbital inclination
+    phase: float
+        The phase at coalescence
+    ra: float
+        The right ascension of the binary
+    dec: float
+        The declination of the object
+    geocent_time: float
+        The time at coalescence
+    psi: float
+        Orbital polarisation
+    kwargs: dict
+        Optional keyword arguments
+
+    Returns
+    -------
+    dict: A dictionary with the plus and cross polarisation strain modes
+    """
+    
+    waveform_kwargs = dict(waveform_approximant='EccentricFD', reference_frequency=10.0,
+                           minimum_frequency=10.0)
+    waveform_kwargs.update(kwargs)
+    waveform_approximant = waveform_kwargs['waveform_approximant']
+    reference_frequency = waveform_kwargs['reference_frequency']
+    minimum_frequency = waveform_kwargs['minimum_frequency']
+
+    if mass_2 > mass_1:
+        return None
+
+    luminosity_distance = luminosity_distance * 1e6 * utils.parsec
+    mass_1 = mass_1 * utils.solar_mass
+    mass_2 = mass_2 * utils.solar_mass
+    
+    spin_1x = 0.0
+    spin_1y = 0.0
+    spin_1z = 0.0 
+    spin_2x = 0.0
+    spin_2y = 0.0
+    spin_2z = 0.0
+    
+    longitude_ascending_nodes = 0.0
+    mean_per_ano = 0.0
+
+    waveform_dictionary = None
+
+    approximant = lalsim.GetApproximantFromString(waveform_approximant)
+
+    maximum_frequency = frequency_array[-1]
+    delta_frequency = frequency_array[1] - frequency_array[0]
+
+    hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
+        mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, luminosity_distance, iota, phase,
+        longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
+        minimum_frequency, maximum_frequency, reference_frequency,
+        waveform_dictionary, approximant)
+
+    h_plus = hplus.data.data
+    h_cross = hcross.data.data
+
+    return {'plus': h_plus, 'cross': h_cross}
+
 
 def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi):
     tau = Q / (np.sqrt(2.0)*np.pi*frequency)

tupak/gw/utils.py

@@ -2,9 +2,9 @@ from __future__ import division
 import os
 
 import numpy as np
-from scipy import signal
 
-from tupak.core.utils import gps_time_to_gmst, ra_dec_to_theta_phi, speed_of_light, nfft, logger
+from ..core.utils import (gps_time_to_gmst, ra_dec_to_theta_phi, speed_of_light,
+                          nfft, logger)
 
 try:
     from gwpy.signal import filter_design