.gitlab-ci.yml

@@ -79,6 +79,28 @@ python-3.7:
       - coverage_badge.svg
       - docs/_build/html/
 
+# Tests run at a fixed schedule rather than on push
+scheduled-python-3.7:
+  stage: test
+  image: bilbydev/bilby-test-suite-python37
+  only:
+    - schedules
+  before_script:
+    # Install the dependencies specified in the Pipfile
+    - pipenv install --three --python=/opt/conda/bin/python --system --deploy
+  script:
+    - python setup.py install
+
+    # Run pyflakes
+    - flake8 .
+
+    # Run tests
+    - pytest
+
+    # Run tests which are only done on schedule
+    - pytest test/example_test.py
+    - pytest test/gw_example_test.py
+
 pages:
   stage: deploy
   dependencies:

CHANGELOG.md

@@ -14,6 +14,7 @@
 ## [0.4.0] 2019-02-15
 
 ### Changed
+- Changed the logic around redundancy tests in the `PriorDict` classes
 - Fixed an accidental addition of astropy as a first-class dependency and added a check for missing dependencies to the C.I.
 - Fixed a bug in the "create-your-own-time-domain-model" example
 - Added citation guide to the readme
@@ -33,7 +34,7 @@
 - Improve the load_data_from_cache_file method
 
 ### Removed
--
+- Removed deprecated `PriorSet` classes. Use `PriorDict` instead.
 
 ## [0.3.5] 2019-01-25
 

bilby/core/prior.py

@@ -245,10 +245,28 @@ class PriorDict(OrderedDict):
         """
         return [self[key].rescale(sample) for key, sample in zip(keys, theta)]
 
-    def test_redundancy(self, key):
+    def test_redundancy(self, key, disable_logging=False):
         """Empty redundancy test, should be overwritten in subclasses"""
         return False
 
+    def test_has_redundant_keys(self):
+        """
+        Test whether there are redundant keys in self.
+
+        Return
+        ------
+        bool: Whether there are redundancies or not
+        """
+        redundant = False
+        for key in self:
+            temp = self.copy()
+            del temp[key]
+            if temp.test_redundancy(key, disable_logging=True):
+                logger.warning('{} is a redundant key in this {}.'
+                               .format(key, self.__class__.__name__))
+                redundant = True
+        return redundant
+
     def copy(self):
         """
         We have to overwrite the copy method as it fails due to the presence of

bilby/core/result.py

@@ -8,6 +8,7 @@ import numpy as np
 import deepdish
 import pandas as pd
 import corner
+import json
 import scipy.stats
 import matplotlib
 import matplotlib.pyplot as plt
@@ -19,7 +20,7 @@ from .utils import (logger, infer_parameters_from_function,
 from .prior import Prior, PriorDict, DeltaFunction
 
 
-def result_file_name(outdir, label):
+def result_file_name(outdir, label, extension='json'):
     """ Returns the standard filename used for a result file
 
     Parameters
@@ -28,17 +29,27 @@ def result_file_name(outdir, label):
         Name of the output directory
     label: str
         Naming scheme of the output file
+    extension: str, optional
+        Whether to save as `hdf5` or `json`
 
     Returns
     -------
     str: File name of the output file
     """
-    return '{}/{}_result.h5'.format(outdir, label)
+    if extension == 'hdf5':
+        return '{}/{}_result.h5'.format(outdir, label)
+    else:
+        return '{}/{}_result.json'.format(outdir, label)
 
 
 def read_in_result(filename=None, outdir=None, label=None):
-    """ Wrapper to bilby.core.result.Result.from_hdf5 """
-    return Result.from_hdf5(filename=filename, outdir=outdir, label=label)
+    """ Wrapper to bilby.core.result.Result.from_hdf5
+        or bilby.core.result.Result.from_json """
+    try:
+        result = Result.from_json(filename=filename, outdir=outdir, label=label)
+    except (IOError, ValueError):
+        result = Result.from_hdf5(filename=filename, outdir=outdir, label=label)
+    return result
 
 
 class Result(object):
@@ -155,7 +166,7 @@ class Result(object):
             if (outdir is None) and (label is None):
                 raise ValueError("No information given to load file")
             else:
-                filename = result_file_name(outdir, label)
+                filename = result_file_name(outdir, label, extension='hdf5')
         if os.path.isfile(filename):
             dictionary = deepdish.io.load(filename)
             # Some versions of deepdish/pytables return the dictionanary as
@@ -169,6 +180,50 @@ class Result(object):
         else:
             raise IOError("No result '{}' found".format(filename))
 
+    @classmethod
+    def from_json(cls, filename=None, outdir=None, label=None):
+        """ Read in a saved .json data file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+        outdir, label: str
+            If given, use the default naming convention for saved results file
+
+        Returns
+        -------
+        result: bilby.core.result.Result
+
+        Raises
+        -------
+        ValueError: If no filename is given and either outdir or label is None
+                    If no bilby.core.result.Result is found in the path
+
+        """
+        if filename is None:
+            if (outdir is None) and (label is None):
+                raise ValueError("No information given to load file")
+            else:
+                filename = result_file_name(outdir, label)
+        if os.path.isfile(filename):
+            dictionary = json.load(open(filename, 'r'))
+            for key in dictionary.keys():
+                # Convert some dictionaries back to DataFrames
+                if key in ['posterior', 'nested_samples']:
+                    dictionary[key] = pd.DataFrame.from_dict(dictionary[key])
+                # Convert the loaded priors to bilby prior type
+                if key == 'priors':
+                    for param in dictionary[key].keys():
+                        dictionary[key][param] = str(dictionary[key][param])
+                    dictionary[key] = PriorDict(dictionary[key])
+            try:
+                return cls(**dictionary)
+            except TypeError as e:
+                raise IOError("Unable to load dictionary, error={}".format(e))
+        else:
+            raise IOError("No result '{}' found".format(filename))
+
     def __str__(self):
         """Print a summary """
         if getattr(self, 'posterior', None) is not None:
@@ -303,9 +358,9 @@ class Result(object):
                 pass
         return dictionary
 
-    def save_to_file(self, overwrite=False, outdir=None):
+    def save_to_file(self, overwrite=False, outdir=None, extension='json'):
         """
-        Writes the Result to a deepdish h5 file
+        Writes the Result to a json or deepdish h5 file
 
         Parameters
         ----------
@@ -314,9 +369,11 @@ class Result(object):
             default=False
         outdir: str, optional
             Path to the outdir. Default is the one stored in the result object.
+        extension: str, optional
+            Whether to save as hdf5 instead of json
         """
         outdir = self._safe_outdir_creation(outdir, self.save_to_file)
-        file_name = result_file_name(outdir, self.label)
+        file_name = result_file_name(outdir, self.label, extension)
 
         if os.path.isfile(file_name):
             if overwrite:
@@ -341,8 +398,19 @@ class Result(object):
                 if hasattr(dictionary['sampler_kwargs'][key], '__call__'):
                     dictionary['sampler_kwargs'][key] = str(dictionary['sampler_kwargs'])
 
+        # Convert to json saveable format
+        if extension != 'hdf5':
+            for key in dictionary.keys():
+                if isinstance(dictionary[key], pd.core.frame.DataFrame):
+                    dictionary[key] = dictionary[key].to_dict()
+                elif isinstance(dictionary[key], np.ndarray):
+                    dictionary[key] = dictionary[key].tolist()
+
         try:
-            deepdish.io.save(file_name, dictionary)
+            if extension == 'hdf5':
+                deepdish.io.save(file_name, dictionary)
+            else:
+                json.dump(dictionary, open(file_name, 'w'), indent=2)
         except Exception as e:
             logger.error("\n\n Saving the data has failed with the "
                          "following message:\n {} \n\n".format(e))

bilby/core/sampler/__init__.py

@@ -85,6 +85,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
         overwritten.
     save: bool
         If true, save the priors and results to disk.
+        If hdf5, save as an hdf5 file instead of json.
     result_class: bilby.core.result.Result, or child of
         The result class to use. By default, `bilby.core.result.Result` is used,
         but objects which inherit from this class can be given providing
@@ -183,7 +184,10 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
 
     result.samples_to_posterior(likelihood=likelihood, priors=priors,
                                 conversion_function=conversion_function)
-    if save:
+    if save == 'hdf5':
+        result.save_to_file(extension='hdf5')
+        logger.info("Results saved to {}/".format(outdir))
+    elif save:
         result.save_to_file()
         logger.info("Results saved to {}/".format(outdir))
     if plot:

bilby/core/sampler/base_sampler.py

@@ -241,6 +241,10 @@ class Sampler(object):
             Likelihood can't be evaluated.
 
         """
+
+        if self.priors.test_has_redundant_keys():
+            raise IllegalSamplingSetError("Your sampling set contains redundant parameters.")
+
         self._check_if_priors_can_be_sampled()
         try:
             t1 = datetime.datetime.now()
@@ -502,3 +506,7 @@ class SamplerError(Error):
 
 class SamplerNotInstalledError(SamplerError):
     """ Base class for Error raised by not installed samplers """
+
+
+class IllegalSamplingSetError(Error):
+    """ Class for illegal sets of sampling parameters """

bilby/core/sampler/emcee.py

@@ -51,6 +51,11 @@ class Emcee(MCMCSampler):
                  use_ratio=False, plot=False, skip_import_verification=False,
                  pos0=None, nburn=None, burn_in_fraction=0.25, resume=True,
                  burn_in_act=3, **kwargs):
+        import emcee
+        if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'):
+            self.prerelease = True
+        else:
+            self.prerelease = False
         MCMCSampler.__init__(
             self, likelihood=likelihood, priors=priors, outdir=outdir,
             label=label, use_ratio=use_ratio, plot=plot,
@@ -82,7 +87,6 @@ class Emcee(MCMCSampler):
     @property
     def sampler_function_kwargs(self):
         import emcee
-
         keys = ['lnprob0', 'rstate0', 'blobs0', 'iterations', 'thin', 'storechain', 'mh_proposal']
 
         # updated function keywords for emcee > v2.2.1
@@ -93,7 +97,7 @@ class Emcee(MCMCSampler):
         function_kwargs = {key: self.kwargs[key] for key in keys if key in self.kwargs}
         function_kwargs['p0'] = self.pos0
 
-        if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'):
+        if self.prerelease:
             if function_kwargs['mh_proposal'] is not None:
                 logger.warning("The 'mh_proposal' option is no longer used "
                                "in emcee v{}, and will be ignored.".format(emcee.__version__))
@@ -109,8 +113,6 @@ class Emcee(MCMCSampler):
 
     @property
     def sampler_init_kwargs(self):
-        import emcee
-
         init_kwargs = {key: value
                        for key, value in self.kwargs.items()
                        if key not in self.sampler_function_kwargs}
@@ -122,7 +124,7 @@ class Emcee(MCMCSampler):
         updatekeys = {'dim': 'ndim',
                       'lnpostfn': 'log_prob_fn'}
 
-        if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'):
+        if self.prerelease:
             for key in updatekeys:
                 if key in init_kwargs:
                     init_kwargs[updatekeys[key]] = init_kwargs.pop(key)
@@ -193,8 +195,10 @@ class Emcee(MCMCSampler):
 
         for sample in tqdm(sampler.sample(**self.sampler_function_kwargs),
                            total=self.nsteps):
-            points = np.hstack([sample[0], np.array(sample[3])])
-            # import IPython; IPython.embed()
+            if self.prerelease:
+                points = np.hstack([sample.coords, sample.blobs])
+            else:
+                points = np.hstack([sample[0], np.array(sample[3])])
             with open(out_file, "a") as ff:
                 for ii, point in enumerate(points):
                     ff.write(template.format(ii, *point))

bilby/gw/detector.py

@@ -147,38 +147,8 @@ class InterferometerList(list):
         if utils.command_line_args.test:
             return
 
-        fig = plt.figure()
-        for ii, interferometer in enumerate(self):
-            ax = fig.add_subplot(len(self) // 2, 2, ii + 1)
-            ax.loglog(interferometer.frequency_array,
-                      gwutils.asd_from_freq_series(freq_data=interferometer.frequency_domain_strain,
-                                                   df=(interferometer.frequency_array[1] -
-                                                       interferometer.frequency_array[0])),
-                      color='C0', label=interferometer.name)
-            ax.loglog(interferometer.frequency_array,
-                      interferometer.amplitude_spectral_density_array,
-                      color='C1', lw=0.5, label=interferometer.name + ' ASD')
-            ax.grid('on')
-            ax.set_ylabel(r'strain [strain/$\sqrt{\rm Hz}$]')
-            ax.set_xlabel(r'frequency [Hz]')
-            ax.set_xlim(20, 2000)
-            ax.legend(loc='best')
-        if signal is not None:
-            ax.loglog(self.frequency_array,
-                      gwutils.asd_from_freq_series(freq_data=signal,
-                                                   df=(self.frequency_array[1] -
-                                                       self.frequency_array[0])
-                                                   ),
-                      color='C2',
-                      label='Signal')
-        fig.tight_layout()
-        if label is None:
-            fig.savefig(
-                '{}/frequency_domain_data.png'.format(outdir))
-        else:
-            fig.savefig(
-                '{}/{}_frequency_domain_data.png'.format(
-                    outdir, label))
+        for interferometer in self:
+            interferometer.plot_data(signal=signal, outdir=outdir, label=label)
 
     @property
     def number_of_interferometers(self):
@@ -468,7 +438,9 @@ class InterferometerStrainData(object):
         strain = strain.filter(bp, filtfilt=True)
         self._time_domain_strain = strain.value
 
-    def create_power_spectral_density(self, fft_length, name='unknown', outdir=None):
+    def create_power_spectral_density(
+            self, fft_length, overlap=0, name='unknown', outdir=None,
+            analysis_segment_start_time=None):
         """ Use the time domain strain to generate a power spectral density
 
         This create a Tukey-windowed power spectral density and writes it to a
@@ -478,12 +450,17 @@ class InterferometerStrainData(object):
         ----------
         fft_length: float
             Duration of the analysis segment.
+        overlap: float
+            Number of seconds of overlap between FFTs.
         name: str
             The name of the detector, used in storing the PSD. Defaults to
             "unknown".
         outdir: str
             The output directory to write the PSD file too. If not given,
             the PSD will not be written to file.
+        analysis_segment_start_time: float
+            The start time of the analysis segment, if given, this data will
+            be removed before creating the PSD.
 
         Returns
         -------
@@ -491,11 +468,28 @@ class InterferometerStrainData(object):
             The frequencies and power spectral density array
 
         """
-        strain = gwpy.timeseries.TimeSeries(self.time_domain_strain, sample_rate=self.sampling_frequency)
+
+        data = self.time_domain_strain
+
+        if analysis_segment_start_time is not None:
+            analysis_segment_end_time = analysis_segment_start_time + fft_length
+            inside = (analysis_segment_start_time > self.time_array[0] +
+                      analysis_segment_end_time < self.time_array[-1])
+            if inside:
+                logger.info("Removing analysis segment data from the PSD data")
+                idxs = (
+                    (self.time_array < analysis_segment_start_time) +
+                    (self.time_array > analysis_segment_end_time))
+                data = data[idxs]
+
+        # WARNING this line can cause issues if the data is non-contiguous
+        strain = gwpy.timeseries.TimeSeries(data=data, sample_rate=self.sampling_frequency)
         psd_alpha = 2 * self.roll_off / fft_length
-        logger.info("Creating PSD with non-overlapping tukey window, alpha={}, roll off={}".format(
-            psd_alpha, self.roll_off))
-        psd = strain.psd(fftlength=fft_length, overlap=0, window=('tukey', psd_alpha))
+        logger.info(
+            "Tukey window PSD data with alpha={}, roll off={}".format(
+                psd_alpha, self.roll_off))
+        psd = strain.psd(
+            fftlength=fft_length, overlap=overlap, window=('tukey', psd_alpha))
 
         if outdir:
             psd_file = '{}/{}_PSD_{}_{}.txt'.format(outdir, name, self.start_time, self.duration)
@@ -1831,7 +1825,8 @@ class PowerSpectralDensity(object):
     @staticmethod
     def from_frame_file(frame_file, psd_start_time, psd_duration,
                         fft_length=4, sampling_frequency=4096, roll_off=0.2,
-                        channel=None):
+                        overlap=0, channel=None, name=None, outdir=None,
+                        analysis_segment_start_time=None):
         """ Generate power spectral density from a frame file
 
         Parameters
@@ -1849,15 +1844,25 @@ class PowerSpectralDensity(object):
             This is twice the maximum frequency.
         roll_off: float, optional
             Rise time in seconds of tukey window.
+        overlap: float,
+            Number of seconds of overlap between FFTs.
         channel: str, optional
             Name of channel to use to generate PSD.
+        name, outdir: str, optional
+            Name (and outdir) of the detector for which a PSD is to be
+            generated.
+        analysis_segment_start_time: float, optional
+            The start time of the analysis segment, if given, this data will
+            be removed before creating the PSD.
 
         """
         strain = InterferometerStrainData(roll_off=roll_off)
         strain.set_from_frame_file(
             frame_file, start_time=psd_start_time, duration=psd_duration,
             channel=channel, sampling_frequency=sampling_frequency)
-        frequency_array, psd_array = strain.create_power_spectral_density(fft_length=fft_length)
+        frequency_array, psd_array = strain.create_power_spectral_density(
+            fft_length=fft_length, name=name, outdir=outdir, overlap=overlap,
+            analysis_segment_start_time=analysis_segment_start_time)
         return PowerSpectralDensity(frequency_array=frequency_array, psd_array=psd_array)
 
     @staticmethod
@@ -2286,23 +2291,30 @@ def get_event_data(
 
 
 def load_data_from_cache_file(
-        cache_file, start_time, segment_duration, psd_duration,
-        channel_name=None, sampling_frequency=4096):
+        cache_file, start_time, segment_duration, psd_duration, psd_start_time,
+        channel_name=None, sampling_frequency=4096, roll_off=0.2,
+        overlap=0, outdir=None):
     """ Helper routine to generate an interferometer from a cache file
 
     Parameters
     ----------
     cache_file: str
         Path to the location of the cache file
-    start_time: float
-        GPS start time of the segment
+    start_time, psd_start_time: float
+        GPS start time of the segment and data stretch used for the PSD
     segment_duration, psd_duration: float
         Segment duration and duration of data to use to generate the PSD (in
         seconds).
+    roll_off: float, optional
+        Rise time in seconds of tukey window.
+    overlap: float,
+        Number of seconds of overlap between FFTs.
     channel_name: str
         Channel name
     sampling_frequency: int
         Sampling frequency
+    outdir: str, optional
+        The output directory in which the data is saved
 
     Returns
     -------
@@ -2313,7 +2325,6 @@ def load_data_from_cache_file(
 
     data_set = False
     psd_set = False
-    psd_start_time = start_time - psd_duration
 
     with open(cache_file, 'r') as ff:
         for line in ff:
@@ -2337,10 +2348,17 @@ def load_data_from_cache_file(
             if not psd_set & psd_in_cache:
                 ifo.power_spectral_density = \
                     PowerSpectralDensity.from_frame_file(
-                        cache.path, psd_start_time=psd_start_time,
+                        cache.path,
+                        psd_start_time=psd_start_time,
                         psd_duration=psd_duration,
+                        fft_length=segment_duration,
+                        sampling_frequency=sampling_frequency,
+                        roll_off=roll_off,
+                        overlap=overlap,
                         channel=channel_name,
-                        sampling_frequency=sampling_frequency)
+                        name=cache.observatory,
+                        outdir=outdir,
+                        analysis_segment_start_time=start_time)
                 psd_set = True
     if data_set and psd_set:
         return ifo

bilby/gw/likelihood.py

@@ -412,10 +412,13 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
     """
     def __init__(self, interferometers, waveform_generator,
-                 linear_matrix, quadratic_matrix, priors):
+                 linear_matrix, quadratic_matrix, priors,
+                 distance_marginalization=False, phase_marginalization=False):
         GravitationalWaveTransient.__init__(
             self, interferometers=interferometers,
-            waveform_generator=waveform_generator, priors=priors)
+            waveform_generator=waveform_generator, priors=priors,
+            distance_marginalization=distance_marginalization,
+            phase_marginalization=phase_marginalization)
 
         if isinstance(linear_matrix, str):
             logger.info("Loading linear matrix from {}".format(linear_matrix))
@@ -434,7 +437,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
     def log_likelihood_ratio(self):
         optimal_snr_squared = 0.
-        matched_filter_snr_squared = 0.
+        d_inner_h = 0.
 
         indices, in_bounds = self._closest_time_indices(
             self.parameters['geocent_time'] - self.interferometers.start_time)
@@ -470,19 +473,27 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             if not in_bounds:
                 return np.nan_to_num(-np.inf)
 
-            matched_filter_snr_squared_array = np.einsum(
+            d_inner_h_tc_array = np.einsum(
                 'i,ji->j', np.conjugate(h_plus_linear + h_cross_linear),
                 self.weights[ifo.name + '_linear'][indices])
 
-            matched_filter_snr_squared += interp1d(
+            d_inner_h += interp1d(
                 self.time_samples[indices],
-                matched_filter_snr_squared_array, kind='quadratic')(ifo_time)
+                d_inner_h_tc_array, kind='quadratic')(ifo_time)
 
             optimal_snr_squared += \
                 np.vdot(np.abs(h_plus_quadratic + h_cross_quadratic)**2,
                         self.weights[ifo.name + '_quadratic'])
 
-        log_l = matched_filter_snr_squared - optimal_snr_squared / 2
+        if self.distance_marginalization:
+            rho_mf_ref, rho_opt_ref = self._setup_rho(d_inner_h, optimal_snr_squared)
+            if self.phase_marginalization:
+                rho_mf_ref = abs(rho_mf_ref)
+            log_l = self._interp_dist_margd_loglikelihood(rho_mf_ref.real, rho_opt_ref)[0]
+        else:
+            if self.phase_marginalization:
+                d_inner_h = self._bessel_function_interped(abs(d_inner_h))
+            log_l = d_inner_h - optimal_snr_squared / 2
 
         return log_l.real
 

bilby/gw/prior.py

@@ -207,63 +207,48 @@ class BBHPriorDict(PriorDict):
                 filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
         PriorDict.__init__(self, dictionary=dictionary, filename=filename)
 
-    def test_redundancy(self, key):
+    def test_redundancy(self, key, disable_logging=False):
         """
         Test whether adding the key would add be redundant.
+        Already existing keys return True.
 
         Parameters
         ----------
         key: str
             The key to test.
+        disable_logging: bool, optional
+            Disable logging in this function call. Default is False.
 
         Return
         ------
         redundant: bool
             Whether the key is redundant or not
         """
-        redundant = False
         if key in self:
             logger.debug('{} already in prior'.format(key))
-            return redundant
+            return True
+
         mass_parameters = {'mass_1', 'mass_2', 'chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio'}
-        spin_magnitude_parameters = {'a_1', 'a_2'}
         spin_tilt_1_parameters = {'tilt_1', 'cos_tilt_1'}
         spin_tilt_2_parameters = {'tilt_2', 'cos_tilt_2'}
         spin_azimuth_parameters = {'phi_1', 'phi_2', 'phi_12', 'phi_jl'}
         inclination_parameters = {'iota', 'cos_iota'}
         distance_parameters = {'luminosity_distance', 'comoving_distance', 'redshift'}
 
-        for parameter_set in [mass_parameters, spin_magnitude_parameters, spin_azimuth_parameters]:
-            if key in parameter_set:
-                if len(parameter_set.intersection(self)) > 2:
-                    redundant = True
-                    logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                        parameter_set.intersection(self)))
-                    break
-            elif len(parameter_set.intersection(self)) == 2:
-                redundant = True
-                break
-        for parameter_set in [inclination_parameters, distance_parameters, spin_tilt_1_parameters,
-                              spin_tilt_2_parameters]:
+        for independent_parameters, parameter_set in \
+                zip([2, 2, 1, 1, 1, 1],
+                    [mass_parameters, spin_azimuth_parameters,
+                     spin_tilt_1_parameters, spin_tilt_2_parameters,
+                     inclination_parameters, distance_parameters]):
             if key in parameter_set:
-                if len(parameter_set.intersection(self)) > 1:
-                    redundant = True
-                    logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                        parameter_set.intersection(self)))
-                    break
-                elif len(parameter_set.intersection(self)) == 1:
-                    redundant = True
-                    break
-
-        return redundant
-
-
-class BBHPriorSet(BBHPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BBHPriorDict INSTEAD"""
-        logger.warning("The name 'BBHPriorSet' is deprecated use 'BBHPriorDict' instead")
-        super(BBHPriorSet, self).__init__(dictionary, filename)
+                if len(parameter_set.intersection(self)) >= independent_parameters:
+                    logger.disabled = disable_logging
+                    logger.warning('{} already in prior. '
+                                   'This may lead to unexpected behaviour.'
+                                   .format(parameter_set.intersection(self)))
+                    logger.disabled = False
+                    return True
+        return False
 
 
 class BNSPriorDict(PriorDict):
@@ -286,34 +271,32 @@ class BNSPriorDict(PriorDict):
                 filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
         PriorDict.__init__(self, dictionary=dictionary, filename=filename)
 
-    def test_redundancy(self, key):
-        logger.info("Performing redundancy check using BBHPriorDict().test_redundancy")
-        bbh_redundancy = BBHPriorDict().test_redundancy(key)
+    def test_redundancy(self, key, disable_logging=False):
+        logger.disabled = disable_logging
+        logger.info("Performing redundancy check using BBHPriorDict(self).test_redundancy")
+        logger.disabled = False
+        bbh_redundancy = BBHPriorDict(self).test_redundancy(key)
+
         if bbh_redundancy:
             return True
         redundant = False
 
-        tidal_parameters =\
+        tidal_parameters = \
             {'lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda'}
 
         if key in tidal_parameters:
             if len(tidal_parameters.intersection(self)) > 2:
                 redundant = True
-                logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                    tidal_parameters.intersection(self)))
+                logger.disabled = disable_logging
+                logger.warning('{} already in prior. '
+                               'This may lead to unexpected behaviour.'
+                               .format(tidal_parameters.intersection(self)))
+                logger.disabled = False
             elif len(tidal_parameters.intersection(self)) == 2:
                 redundant = True
         return redundant
 
 
-class BNSPriorSet(BNSPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BNSPriorDict INSTEAD"""
-        super(BNSPriorSet, self).__init__(dictionary, filename)
-        logger.warning("The name 'BNSPriorSet' is deprecated use 'BNSPriorDict' instead")
-
-
 Prior._default_latex_labels = {
     'mass_1': '$m_1$',
     'mass_2': '$m_2$',
@@ -418,13 +401,13 @@ class CalibrationPriorDict(PriorDict):
         nodes = np.logspace(np.log10(minimum_frequency),
                             np.log10(maximum_frequency), n_nodes)
 
-        amplitude_mean_nodes =\
+        amplitude_mean_nodes = \
             UnivariateSpline(frequency_array, amplitude_median)(nodes)
-        amplitude_sigma_nodes =\
+        amplitude_sigma_nodes = \
             UnivariateSpline(frequency_array, amplitude_sigma)(nodes)
-        phase_mean_nodes =\
+        phase_mean_nodes = \
             UnivariateSpline(frequency_array, phase_median)(nodes)
-        phase_sigma_nodes =\
+        phase_sigma_nodes = \
             UnivariateSpline(frequency_array, phase_sigma)(nodes)
 
         prior = CalibrationPriorDict()
@@ -506,11 +489,3 @@ class CalibrationPriorDict(PriorDict):
                                         latex_label=latex_label)
 
         return prior
-
-
-class CalibrationPriorSet(CalibrationPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BNSPriorDict INSTEAD"""
-        super(CalibrationPriorSet, self).__init__(dictionary, filename)
-        logger.warning("The name 'CalibrationPriorSet' is deprecated use 'CalibrationPriorDict' instead")

bilby/gw/source.py

@@ -8,7 +8,9 @@ from .utils import (lalsim_SimInspiralTransformPrecessingNewInitialConditions,
                     lalsim_GetApproximantFromString,
                     lalsim_SimInspiralChooseFDWaveform,
                     lalsim_SimInspiralWaveformParamsInsertTidalLambda1,
-                    lalsim_SimInspiralWaveformParamsInsertTidalLambda2)
+                    lalsim_SimInspiralWaveformParamsInsertTidalLambda2,
+                    lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame,
+                    lalsim_SimIMRPhenomPFrequencySequence)
 
 try:
     import lal
@@ -410,25 +412,25 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
         spin_2z = a_2
     else:
         iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
-            lalsim.SimInspiralTransformPrecessingNewInitialConditions(
+            lalsim_SimInspiralTransformPrecessingNewInitialConditions(
                 iota, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
                 reference_frequency, phase)
 
     chi_1_l, chi_2_l, chi_p, theta_jn, alpha, phase_aligned, zeta =\
-        lalsim.SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+        lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
             mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
             spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version)
 
     waveform_polarizations = dict()
 
-    h_linear_plus, h_linear_cross = lalsim.SimIMRPhenomPFrequencySequence(
+    h_linear_plus, h_linear_cross = lalsim_SimIMRPhenomPFrequencySequence(
         frequency_nodes_linear, chi_1_l, chi_2_l, chi_p, theta_jn,
         mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency, version, None)
-    h_quadratic_plus, h_quadratic_cross = lalsim.SimIMRPhenomPFrequencySequence(
+        alpha, phase_aligned, reference_frequency, version)
+    h_quadratic_plus, h_quadratic_cross = lalsim_SimIMRPhenomPFrequencySequence(
         frequency_nodes_quadratic, chi_1_l, chi_2_l, chi_p, theta_jn,
         mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency, version, None)
+        alpha, phase_aligned, reference_frequency, version)
 
     waveform_polarizations['linear'] = dict(
         plus=(np.cos(2 * zeta) * h_linear_plus.data.data +

bilby/gw/utils.py

@@ -130,6 +130,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
     elif mode.lower() == 'breathing':
         return np.einsum('i,j->ij', m, m) + np.einsum('i,j->ij', n, n)
 
+    # Calculating omega here to avoid calculation when model in [plus, cross, breathing]
     omega = np.cross(m, n)
     if mode.lower() == 'longitudinal':
         return np.sqrt(2) * np.einsum('i,j->ij', omega, omega)
@@ -138,8 +139,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
     elif mode.lower() == 'y':
         return np.einsum('i,j->ij', n, omega) + np.einsum('i,j->ij', omega, n)
     else:
-        logger.warning("{} not a polarization mode!".format(mode))
-        return None
+        raise ValueError("{} not a polarization mode!".format(mode))
 
 
 def get_vertex_position_geocentric(latitude, longitude, elevation):
@@ -380,7 +380,7 @@ def get_open_strain_data(
     return strain
 
 
-def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1, **kwargs):
+def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=0, **kwargs):
     """ A function which accesses the open strain data
 
     This uses `gwpy` to download the open data and then saves a cached copy for
@@ -416,23 +416,22 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1
         except RuntimeError:
             logger.warning('Channel {} not found. Trying preset channel names'.format(channel))
 
-    while not loaded:
-        ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
-                              'DCH-CLEAN_STRAIN_C02']
-        virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
-        channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
-        for detector in channel_types.keys():
-            for channel_type in channel_types[detector]:
-                if loaded:
-                    break
-                channel = '{}:{}'.format(detector, channel_type)
-                try:
-                    strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
-                                             **kwargs)
-                    loaded = True
-                    logger.info('Successfully read strain data for channel {}.'.format(channel))
-                except RuntimeError:
-                    pass
+    ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
+                          'DCH-CLEAN_STRAIN_C02']
+    virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
+    channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
+    for detector in channel_types.keys():
+        for channel_type in channel_types[detector]:
+            if loaded:
+                break
+            channel = '{}:{}'.format(detector, channel_type)
+            try:
+                strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
+                                         **kwargs)
+                loaded = True
+                logger.info('Successfully read strain data for channel {}.'.format(channel))
+            except RuntimeError:
+                pass
 
     if loaded:
         return strain
@@ -770,6 +769,33 @@ def lalsim_SimInspiralChooseFDWaveform(
         waveform_dictionary, approximant)
 
 
+def lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
+        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version):
+    [mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
+     spin_1y, spin_1z, spin_2x, spin_2y, spin_2z] = convert_args_list_to_float(
+        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
+        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z)
+    return lalsim.SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
+        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version)
+
+
+def lalsim_SimIMRPhenomPFrequencySequence(
+        frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn,
+        mass_1, mass_2, luminosity_distance,
+        alpha, phase_aligned, reference_frequency, version):
+    [chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance,
+     alpha, phase_aligned, reference_frequency] = convert_args_list_to_float(
+        chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance,
+        alpha, phase_aligned, reference_frequency)
+
+    return lalsim.SimIMRPhenomPFrequencySequence(
+        frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2,
+        luminosity_distance, alpha, phase_aligned, reference_frequency, version,
+        None)
+
+
 def lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
         waveform_dictionary, lambda_1):
     try:

setup.cfg

@@ -6,6 +6,8 @@ ignore = E129 W504 W605
 [tool:pytest]
 addopts =
     --ignore test/other_test.py
+    --ignore test/gw_example_test.py
+    --ignore test/example_test.py
 
 [metadata]
 license_file = LICENSE.md

test/detector_test.py

@@ -39,7 +39,8 @@ class TestInterferometer(unittest.TestCase):
                                                     xarm_tilt=self.xarm_tilt, yarm_tilt=self.yarm_tilt)
         self.ifo.strain_data.set_from_frequency_domain_strain(
             np.linspace(0, 4096, 4097), sampling_frequency=4096, duration=2)
-        bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
+        self.outdir = 'outdir'
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir)
 
     def tearDown(self):
         del self.name
@@ -55,7 +56,7 @@ class TestInterferometer(unittest.TestCase):
         del self.xarm_tilt
         del self.yarm_tilt
         del self.ifo
-        rmtree('outdir')
+        rmtree(self.outdir)
 
     def test_name_setting(self):
         self.assertEqual(self.ifo.name, self.name)
@@ -798,7 +799,8 @@ class TestInterferometerList(unittest.TestCase):
         self.ifo2.strain_data.set_from_frequency_domain_strain(
             self.frequency_arrays, sampling_frequency=4096, duration=2)
         self.ifo_list = bilby.gw.detector.InterferometerList([self.ifo1, self.ifo2])
-        bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
+        self.outdir = 'outdir'
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir)
 
     def tearDown(self):
         del self.frequency_arrays
@@ -829,7 +831,7 @@ class TestInterferometerList(unittest.TestCase):
         del self.ifo1
         del self.ifo2
         del self.ifo_list
-        rmtree('outdir')
+        rmtree(self.outdir)
 
     def test_init_with_string(self):
         with self.assertRaises(TypeError):
@@ -896,6 +898,13 @@ class TestInterferometerList(unittest.TestCase):
         with self.assertRaises(ValueError):
             self.ifo_list.inject_signal(injection_polarizations=None, waveform_generator=None)
 
+    def test_meta_data(self):
+        ifos_list = [self.ifo1, self.ifo2]
+        ifos = bilby.gw.detector.InterferometerList(ifos_list)
+        self.assertTrue(isinstance(ifos.meta_data, dict))
+        meta_data = {ifo.name: ifo.meta_data for ifo in ifos_list}
+        self.assertEqual(ifos.meta_data, meta_data)
+
     @patch.object(bilby.gw.waveform_generator.WaveformGenerator, 'frequency_domain_strain')
     def test_inject_signal_pol_none_calls_frequency_domain_strain(self, m):
         waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
@@ -986,6 +995,15 @@ class TestInterferometerList(unittest.TestCase):
             with self.assertRaises(TypeError):
                 bilby.gw.detector.InterferometerList.from_hdf5(filename)
 
+    def test_plot_data(self):
+        ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+        ifos.set_strain_data_from_power_spectral_densities(2048, 4)
+        ifos.plot_data(outdir=self.outdir)
+
+        ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+        ifos.set_strain_data_from_power_spectral_densities(2048, 4)
+        ifos.plot_data(outdir=self.outdir)
+
 
 class TestPowerSpectralDensityWithoutFiles(unittest.TestCase):
 

test/gw_likelihood_test.py

@@ -509,7 +509,12 @@ class TestROQLikelihood(unittest.TestCase):
             interferometers=ifos, waveform_generator=roq_wfg,
             linear_matrix=linear_matrix_file,
             quadratic_matrix=quadratic_matrix_file, priors=self.priors)
-        pass
+
+        self.roq_phase_like = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+            interferometers=ifos, waveform_generator=roq_wfg,
+            linear_matrix=linear_matrix_file,
+            quadratic_matrix=quadratic_matrix_file,
+            phase_marginalization=True, priors=self.priors.copy())
 
     def tearDown(self):
         pass
@@ -527,6 +532,20 @@ class TestROQLikelihood(unittest.TestCase):
         self.assertEqual(
             self.roq_likelihood.log_likelihood_ratio(), np.nan_to_num(-np.inf))
 
+    def test_phase_marginalisation_roq(self):
+        """Test phase marginalised likelihood matches brute force version"""
+        like = []
+        self.roq_likelihood.parameters.update(self.test_parameters)
+        phases = np.linspace(0, 2 * np.pi, 1000)
+        for phase in phases:
+            self.roq_likelihood.parameters['phase'] = phase
+            like.append(np.exp(self.roq_likelihood.log_likelihood_ratio()))
+
+        marg_like = np.log(np.trapz(like, phases) / (2 * np.pi))
+        self.roq_phase_like.parameters = self.test_parameters.copy()
+        self.assertAlmostEqual(
+            marg_like, self.roq_phase_like.log_likelihood_ratio(), delta=0.5)
+
 
 class TestBBHLikelihoodSetUp(unittest.TestCase):
 

test/gw_prior_test.py

@@ -16,42 +16,174 @@ class TestBBHPriorDict(unittest.TestCase):
         self.base_directory =\
             '/'.join(os.path.dirname(
                 os.path.abspath(sys.argv[0])).split('/')[:-1])
-        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'prior_files/binary_black_holes.prior')
-        self.default_prior = bilby.gw.prior.BBHPriorDict(
-            filename=self.filename)
+        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                                     'prior_files/binary_black_holes.prior')
+        self.bbh_prior_dict = bilby.gw.prior.BBHPriorDict(filename=self.filename)
+        for key, value in self.bbh_prior_dict.items():
+            self.prior_dict[key] = value
 
     def tearDown(self):
         del self.prior_dict
         del self.filename
+        del self.bbh_prior_dict
+        del self.base_directory
 
     def test_create_default_prior(self):
         default = bilby.gw.prior.BBHPriorDict()
-        minima = all([self.default_prior[key].minimum == default[key].minimum
+        minima = all([self.bbh_prior_dict[key].minimum == default[key].minimum
                       for key in default.keys()])
-        maxima = all([self.default_prior[key].maximum == default[key].maximum
+        maxima = all([self.bbh_prior_dict[key].maximum == default[key].maximum
                       for key in default.keys()])
-        names = all([self.default_prior[key].name == default[key].name
+        names = all([self.bbh_prior_dict[key].name == default[key].name
                      for key in default.keys()])
 
         self.assertTrue(all([minima, maxima, names]))
 
     def test_create_from_dict(self):
-        bilby.gw.prior.BBHPriorDict(dictionary=self.prior_dict)
+        new_dict = bilby.gw.prior.BBHPriorDict(dictionary=self.prior_dict)
+        for key in self.bbh_prior_dict:
+            self.assertEqual(self.bbh_prior_dict[key], new_dict[key])
+
+    def test_redundant_priors_not_in_dict_before(self):
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_iota',
+                      'comoving_distance', 'redshift']:
+            self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_redundant_priors_already_in_dict(self):
+        for prior in ['mass_1', 'mass_2', 'tilt_1', 'tilt_2',
+                      'phi_1', 'phi_2', 'iota', 'luminosity_distance']:
+            self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_masses(self):
+        del self.bbh_prior_dict['mass_2']
+        for prior in ['mass_2', 'chirp_mass', 'total_mass', 'mass_ratio',  'symmetric_mass_ratio']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_magnitudes(self):
+        del self.bbh_prior_dict['a_2']
+        self.assertFalse(self.bbh_prior_dict.test_redundancy('a_2'))
+
+    def test_correct_not_redundant_priors_spin_tilt_1(self):
+        del self.bbh_prior_dict['tilt_1']
+        for prior in ['tilt_1', 'cos_tilt_1']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_tilt_2(self):
+        del self.bbh_prior_dict['tilt_2']
+        for prior in ['tilt_2', 'cos_tilt_2']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_azimuth(self):
+        del self.bbh_prior_dict['phi_12']
+        for prior in ['phi_1', 'phi_2', 'phi_12']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_inclination(self):
+        del self.bbh_prior_dict['iota']
+        for prior in ['iota', 'cos_iota']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_distance(self):
+        del self.bbh_prior_dict['luminosity_distance']
+        for prior in ['luminosity_distance', 'comoving_distance',
+                      'redshift']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_add_unrelated_prior(self):
+        self.assertFalse(self.bbh_prior_dict.test_redundancy('abc'))
+
+    def test_test_has_redundant_priors(self):
+        self.assertFalse(self.bbh_prior_dict.test_has_redundant_keys())
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_iota',
+                      'comoving_distance', 'redshift']:
+            self.bbh_prior_dict[prior] = 0
+            self.assertTrue(self.bbh_prior_dict.test_has_redundant_keys())
+            del self.bbh_prior_dict[prior]
+
+
+class TestBNSPriorDict(unittest.TestCase):
 
-    def test_create_from_filename(self):
-        bilby.gw.prior.BBHPriorDict(filename=self.filename)
+    def setUp(self):
+        self.prior_dict = dict()
+        self.base_directory =\
+            '/'.join(os.path.dirname(
+                os.path.abspath(sys.argv[0])).split('/')[:-1])
+        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                                     'prior_files/binary_neutron_stars.prior')
+        self.bns_prior_dict = bilby.gw.prior.BNSPriorDict(filename=self.filename)
+        for key, value in self.bns_prior_dict.items():
+            self.prior_dict[key] = value
 
-    def test_key_in_prior_not_redundant(self):
-        test = self.default_prior.test_redundancy('mass_1')
-        self.assertFalse(test)
+    def tearDown(self):
+        del self.prior_dict
+        del self.filename
+        del self.bns_prior_dict
+        del self.base_directory
 
-    def test_chirp_mass_redundant(self):
-        test = self.default_prior.test_redundancy('chirp_mass')
-        self.assertTrue(test)
+    def test_create_default_prior(self):
+        default = bilby.gw.prior.BNSPriorDict()
+        minima = all([self.bns_prior_dict[key].minimum == default[key].minimum
+                      for key in default.keys()])
+        maxima = all([self.bns_prior_dict[key].maximum == default[key].maximum
+                      for key in default.keys()])
+        names = all([self.bns_prior_dict[key].name == default[key].name
+                     for key in default.keys()])
 
-    def test_comoving_distance_redundant(self):
-        test = self.default_prior.test_redundancy('comoving_distance')
-        self.assertTrue(test)
+        self.assertTrue(all([minima, maxima, names]))
+
+    def test_create_from_dict(self):
+        new_dict = bilby.gw.prior.BNSPriorDict(dictionary=self.prior_dict)
+        self.assertDictEqual(self.bns_prior_dict, new_dict)
+
+    def test_redundant_priors_not_in_dict_before(self):
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio',
+                      'symmetric_mass_ratio', 'cos_iota', 'comoving_distance',
+                      'redshift', 'lambda_tilde', 'delta_lambda']:
+            self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_redundant_priors_already_in_dict(self):
+        for prior in ['mass_1', 'mass_2', 'chi_1', 'chi_2',
+                      'iota', 'luminosity_distance',
+                      'lambda_1', 'lambda_2']:
+            self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_masses(self):
+        del self.bns_prior_dict['mass_2']
+        for prior in ['mass_2', 'chirp_mass', 'total_mass', 'mass_ratio',  'symmetric_mass_ratio']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_magnitudes(self):
+        del self.bns_prior_dict['chi_2']
+        self.assertFalse(self.bns_prior_dict.test_redundancy('chi_2'))
+
+    def test_correct_not_redundant_priors_inclination(self):
+        del self.bns_prior_dict['iota']
+        for prior in ['iota', 'cos_iota']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_distance(self):
+        del self.bns_prior_dict['luminosity_distance']
+        for prior in ['luminosity_distance', 'comoving_distance',
+                      'redshift']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_tidal(self):
+        del self.bns_prior_dict['lambda_1']
+        for prior in['lambda_1', 'lambda_tilde', 'delta_lambda']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_add_unrelated_prior(self):
+        self.assertFalse(self.bns_prior_dict.test_redundancy('abc'))
+
+    def test_test_has_redundant_priors(self):
+        self.assertFalse(self.bns_prior_dict.test_has_redundant_keys())
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_iota', 'comoving_distance', 'redshift']:
+            self.bns_prior_dict[prior] = 0
+            self.assertTrue(self.bns_prior_dict.test_has_redundant_keys())
+            del self.bns_prior_dict[prior]
 
 
 class TestCalibrationPrior(unittest.TestCase):

test/gw_utils_test.py

+from __future__ import absolute_import, division
+
+import unittest
+import os
+from shutil import rmtree
+
+import numpy as np
+import gwpy
+import lal
+import lalsimulation as lalsim
+
+import bilby
+from bilby.gw import utils as gwutils
+
+
+class TestGWUtils(unittest.TestCase):
+
+    def setUp(self):
+        self.outdir = 'outdir'
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir)
+
+    def tearDown(self):
+        try:
+            rmtree(self.outdir)
+        except FileNotFoundError:
+            pass
+
+    def test_asd_from_freq_series(self):
+        freq_data = np.array([1, 2, 3])
+        df = 0.1
+        asd = gwutils.asd_from_freq_series(freq_data, df)
+        self.assertTrue(np.all(asd == freq_data * 2 * df ** 0.5))
+
+    def test_psd_from_freq_series(self):
+        freq_data = np.array([1, 2, 3])
+        df = 0.1
+        psd = gwutils.psd_from_freq_series(freq_data, df)
+        self.assertTrue(np.all(psd == (freq_data * 2 * df ** 0.5)**2))
+
+    def test_time_delay_from_geocenter(self):
+        det1 = np.array([0.1, 0.2, 0.3])
+        det2 = np.array([0.1, 0.2, 0.5])
+        ra = 0.5
+        dec = 0.2
+        time = 10
+        self.assertEqual(
+            gwutils.time_delay_geocentric(det1, det1, ra, dec, time), 0)
+        self.assertEqual(
+            gwutils.time_delay_geocentric(det1, det2, ra, dec, time),
+            1.3253791114055397e-10)
+
+    def test_get_polarization_tensor(self):
+        ra = 1
+        dec = 2.0
+        time = 10
+        psi = 0.1
+        for mode in ['plus', 'cross', 'breathing', 'longitudinal', 'x', 'y']:
+            p = gwutils.get_polarization_tensor(ra, dec, time, psi, mode)
+            self.assertEqual(p.shape, (3, 3))
+        with self.assertRaises(ValueError):
+            gwutils.get_polarization_tensor(ra, dec, time, psi, 'not-a-mode')
+
+    def test_inner_product(self):
+        aa = np.array([1, 2, 3])
+        bb = np.array([5, 6, 7])
+        frequency = np.array([0.2, 0.4, 0.6])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        ip = gwutils.inner_product(aa, bb, frequency, PSD)
+        self.assertEqual(ip, 0)
+
+    def test_noise_weighted_inner_product(self):
+        aa = np.array([1e-23, 2e-23, 3e-23])
+        bb = np.array([5e-23, 6e-23, 7e-23])
+        frequency = np.array([100, 101, 102])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        psd = PSD.power_spectral_density_interpolated(frequency)
+        duration = 4
+        nwip = gwutils.noise_weighted_inner_product(aa, bb, psd, duration)
+        self.assertEqual(nwip, 239.87768033598326)
+
+        self.assertEqual(
+            gwutils.optimal_snr_squared(aa, psd, duration),
+            gwutils.noise_weighted_inner_product(aa, aa, psd, duration))
+
+    def test_matched_filter_snr(self):
+        signal = np.array([1e-23, 2e-23, 3e-23])
+        frequency_domain_strain = np.array([5e-23, 6e-23, 7e-23])
+        frequency = np.array([100, 101, 102])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        psd = PSD.power_spectral_density_interpolated(frequency)
+        duration = 4
+
+        mfsnr = gwutils.matched_filter_snr(
+            signal, frequency_domain_strain, psd, duration)
+        self.assertEqual(mfsnr, 25.510869054168282)
+
+    def test_get_event_time(self):
+        events = ['GW150914', 'GW151012', 'GW151226', 'GW170104', 'GW170608',
+                  'GW170729', 'GW170809', 'GW170814', 'GW170817', 'GW170818',
+                  'GW170823']
+        for event in events:
+            self.assertTrue(isinstance(gwutils.get_event_time(event), float))
+
+        self.assertTrue(gwutils.get_event_time('GW010290') is None)
+
+    def test_read_frame_file(self):
+        start_time = 0
+        end_time = 10
+        channel = 'H1:GDS-CALIB_STRAIN'
+        N = 100
+        times = np.linspace(start_time, end_time, N)
+        data = np.random.normal(0, 1, N)
+        ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+        ts.channel = gwpy.detector.Channel(channel)
+        ts.name = channel
+        filename = os.path.join(self.outdir, 'test.gwf')
+        ts.write(filename, format='gwf')
+
+        # Check reading without time limits
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None, channel=channel)
+        self.assertEqual(strain.channel.name, channel)
+        self.assertTrue(np.all(strain.value==data))
+
+        # Check reading with time limits
+        start_cut = 2
+        end_cut = 8
+        strain = gwutils.read_frame_file(
+            filename, start_time=start_cut, end_time=end_cut, channel=channel)
+        idxs = (times > start_cut) & (times < end_cut)
+        # Dropping the last element - for some reason gwpy drops the last element when reading in data
+        self.assertTrue(np.all(strain.value==data[idxs][:-1]))
+
+        # Check reading with unknown channels
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None)
+        self.assertTrue(np.all(strain.value==data))
+
+        # Check reading with incorrect channel
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None, channel='WRONG')
+        self.assertTrue(np.all(strain.value==data))
+
+        ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+        ts.name = 'NOT-A-KNOWN-CHANNEL'
+        ts.write(filename, format='gwf')
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None)
+        self.assertEqual(strain, None)
+
+    def test_convert_args_list_to_float(self):
+        self.assertEqual(
+            gwutils.convert_args_list_to_float(1, '2', 3.0), [1.0, 2.0, 3.0])
+        with self.assertRaises(ValueError):
+            gwutils.convert_args_list_to_float(1, '2', 'ten')
+
+    def test_lalsim_SimInspiralTransformPrecessingNewInitialConditions(self):
+        a = gwutils.lalsim_SimInspiralTransformPrecessingNewInitialConditions( 
+            0.1, 0, 0.6, 0.5, 0.6, 0.1, 0.8, 30.6, 23.2, 50, 0)
+        self.assertTrue(len(a) == 7)
+
+    def test_get_approximant(self):
+        with self.assertRaises(ValueError):
+            gwutils.lalsim_GetApproximantFromString(10)
+        self.assertEqual(gwutils.lalsim_GetApproximantFromString("IMRPhenomPV2"), 69)
+
+    def test_lalsim_SimInspiralChooseFDWaveform(self):
+        a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+            35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+            45., 0.1, 10, 0.01, 10, 1000, 20, None, 69)
+        self.assertEqual(len(a), 2)
+        self.assertEqual(type(a[0]), lal.COMPLEX16FrequencySeries)
+        self.assertEqual(type(a[1]), lal.COMPLEX16FrequencySeries)
+
+        with self.assertRaises(ValueError):
+            a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+                35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+                45., 0.1, 10, 0.01, 10, 1000, 20, None, 'IMRPhenomPV2')
+
+    def test_lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(self):
+        version = lalsim.IMRPhenomPv2_V
+        a = gwutils.lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+            25.6, 12.2, 50, 0.2, 0, 0, 0, 0, 0, 0, 0, version)
+        self.assertEqual(len(a), 7)
+
+
+if __name__ == '__main__':
+    unittest.main()

test/likelihood_test.py

@@ -31,6 +31,18 @@ class TestLikelihoodBase(unittest.TestCase):
     def test_base_log_likelihood_ratio(self):
         self.assertTrue(np.isnan(self.likelihood.log_likelihood_ratio()))
 
+    def test_meta_data_unset(self):
+        self.assertEqual(self.likelihood.meta_data, None)
+
+    def test_meta_data_set_fail(self):
+        with self.assertRaises(ValueError):
+            self.likelihood.meta_data = 10
+
+    def test_meta_data(self):
+        meta_data = dict(x=1, y=2)
+        self.likelihood.meta_data = meta_data
+        self.assertEqual(self.likelihood.meta_data, meta_data)
+
 
 class TestAnalytical1DLikelihood(unittest.TestCase):
 

test/prior_files/binary_neutron_stars.prior

+# These are the default priors we use for BNS systems.
+# Note that you may wish to use more specific mass and distance parameters.
+# These commands are all known to bilby.gw.prior.
+# Lines beginning "#" are ignored.
+mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$')
+mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
+# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$')
+# total_mass =  Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
+# mass_ratio =  Uniform(name='mass_ratio', minimum=0.5, maximum=1)
+# symmetric_mass_ratio =  Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
+chi_1 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$')
+chi_2 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$')
+luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
+dec =  Cosine(name='dec')
+ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
+iota =  Sine(name='iota')
+# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
+phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000 )
+# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000)
+# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000)

test/prior_test.py

@@ -4,7 +4,7 @@ import unittest
 from mock import Mock
 import numpy as np
 import os
-import shutil
+import copy
 from collections import OrderedDict
 
 
@@ -130,6 +130,7 @@ class TestPriorClasses(unittest.TestCase):
             bilby.core.prior.Gaussian(name='test', unit='unit', mu=0, sigma=1),
             bilby.core.prior.Normal(name='test', unit='unit', mu=0, sigma=1),
             bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=0, minimum=0, maximum=1),
+            bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=-1, minimum=0.5, maximum=1),
             bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=2, minimum=1, maximum=1e2),
             bilby.core.prior.Uniform(name='test', unit='unit', minimum=0, maximum=1),
             bilby.core.prior.LogUniform(name='test', unit='unit', minimum=5e0, maximum=1e2),
@@ -205,6 +206,28 @@ class TestPriorClasses(unittest.TestCase):
                 outside_domain = np.linspace(prior.minimum - 1e4, prior.minimum - 1, 1000)
                 self.assertTrue(all(prior.prob(outside_domain) == 0))
 
+    def test_prob_and_ln_prob(self):
+        for prior in self.priors:
+            sample = prior.sample()
+            self.assertAlmostEqual(np.log(prior.prob(sample)), prior.ln_prob(sample), 12)
+
+    def test_log_normal_fail(self):
+        with self.assertRaises(ValueError):
+            bilby.core.prior.LogNormal(name='test', unit='unit', mu=0, sigma=-1)
+
+    def test_studentt_fail(self):
+        with self.assertRaises(ValueError):
+            bilby.core.prior.StudentT(name='test', unit='unit', df=3, mu=0, scale=-1)
+        with self.assertRaises(ValueError):
+            bilby.core.prior.StudentT(name='test', unit='unit', df=0, mu=0, scale=1)
+
+    def test_beta_fail(self):
+        with self.assertRaises(ValueError):
+            bilby.core.prior.Beta(name='test', unit='unit', alpha=-2.0, beta=2.0),
+
+        with self.assertRaises(ValueError):
+            bilby.core.prior.Beta(name='test', unit='unit', alpha=2.0, beta=-2.0),
+
     def test_probability_in_domain(self):
         """Test that the prior probability is non-negative in domain of validity and zero outside."""
         for prior in self.priors:
@@ -319,6 +342,15 @@ class TestPriorDict(unittest.TestCase):
         del self.default_prior_file
         del self.prior_set_from_file
 
+    def test_copy(self):
+        priors = bilby.core.prior.PriorDict(self.priors)
+        self.assertEqual(priors, priors.copy())
+
+    def test_prior_set(self):
+        priors_dict = bilby.core.prior.PriorDict(self.priors)
+        priors_set = bilby.core.prior.PriorSet(self.priors)
+        self.assertEqual(priors_dict, priors_set)
+
     def test_prior_set_is_ordered_dict(self):
         self.assertIsInstance(self.prior_set_from_dict, OrderedDict)