bilby/gw/likelihood.py

@@ -15,14 +15,14 @@ except ImportError:
     from scipy.misc import logsumexp
 from scipy.special import i0e
 
-from ..core.likelihood import Likelihood, MarginalizedLikelihoodReconstructionError
+from ..core.likelihood import Likelihood
 from ..core.utils import BilbyJsonEncoder, decode_bilby_json
 from ..core.utils import (
     logger, UnsortedInterp2d, create_frequency_series, create_time_series,
     speed_of_light, radius_of_earth)
 from ..core.prior import Interped, Prior, Uniform
 from .detector import InterferometerList
-from .prior import BBHPriorDict, CBCPriorDict
+from .prior import BBHPriorDict, CBCPriorDict, Cosmological
 from .source import lal_binary_black_hole
 from .utils import noise_weighted_inner_product, build_roq_weights, blockwise_dot_product
 from .waveform_generator import WaveformGenerator
@@ -143,6 +143,9 @@ class GravitationalWaveTransient(Likelihood):
                  for distance in self._distance_array])
             self._setup_distance_marginalization(
                 distance_marginalization_lookup_table)
+            for key in ['redshift', 'comoving_distance']:
+                if key in priors:
+                    del priors[key]
             priors['luminosity_distance'] = float(self._ref_dist)
             self._marginalized_parameters.append('luminosity_distance')
 
@@ -214,6 +217,18 @@ class GravitationalWaveTransient(Likelihood):
                 self.priors[key] = Uniform(
                     self.interferometers.start_time,
                     self.interferometers.start_time + self.interferometers.duration)
+            elif key == 'luminosity_distance':
+                for key in ['redshift', 'comoving_distance']:
+                    if key in self.priors:
+                        if not isinstance(self.priors[key], Cosmological):
+                            raise TypeError(
+                                "To marginalize over {}, the prior must be specified as a "
+                                "subclass of bilby.gw.prior.Cosmological.".format(key)
+                            )
+                        self.priors['luminosity_distance'] = self.priors[key].get_corresponding_prior(
+                            'luminosity_distance'
+                        )
+                        del self.priors[key]
             else:
                 self.priors[key] = BBHPriorDict()[key]
 
@@ -392,17 +407,13 @@ class GravitationalWaveTransient(Likelihood):
             np.exp(time_log_like - max(time_log_like)) * time_prior_array)
 
         keep = (time_post > max(time_post) / 1000)
+        if sum(keep) < 3:
+            keep[1:-1] = keep[1:-1] | keep[2:] | keep[:-2]
         time_post = time_post[keep]
         times = times[keep]
 
-        if len(times) > 1:
-            new_time = Interped(times, time_post).sample()
-            return new_time
-        else:
-            raise MarginalizedLikelihoodReconstructionError(
-                "Time posterior reconstruction failed, at least two samples "
-                "are required."
-            )
+        new_time = Interped(times, time_post).sample()
+        return new_time
 
     def generate_distance_sample_from_marginalized_likelihood(
             self, signal_polarizations=None):

bilby/gw/prior.py

@@ -2,10 +2,13 @@ import os
 import copy
 
 import numpy as np
-from scipy.interpolate import InterpolatedUnivariateSpline
+from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
+from scipy.integrate import cumtrapz
+from scipy.stats import norm
 
 from ..core.prior import (PriorDict, Uniform, Prior, DeltaFunction, Gaussian,
-                          Interped, Constraint, conditional_prior_factory)
+                          Interped, Constraint, conditional_prior_factory,
+                          BaseJointPriorDist, JointPrior, JointPriorDistError)
 from ..core.utils import infer_args_from_method, logger
 from .conversion import (
     convert_to_lal_binary_black_hole_parameters,
@@ -133,31 +136,7 @@ class Cosmological(Interped):
 
     @minimum.setter
     def minimum(self, minimum):
-        cosmology = get_cosmology(self.cosmology)
-        self._minimum[self.name] = minimum
-        if self.name == 'redshift':
-            self._minimum['luminosity_distance'] =\
-                cosmology.luminosity_distance(minimum).value
-            self._minimum['comoving_distance'] =\
-                cosmology.comoving_distance(minimum).value
-        elif self.name == 'luminosity_distance':
-            if minimum == 0:
-                self._minimum['redshift'] = 0
-            else:
-                self._minimum['redshift'] = cosmo.z_at_value(
-                    cosmology.luminosity_distance, minimum * self.unit)
-            self._minimum['comoving_distance'] = self._minimum['redshift']
-        elif self.name == 'comoving_distance':
-            if minimum == 0:
-                self._minimum['redshift'] = 0
-            else:
-                self._minimum['redshift'] = cosmo.z_at_value(
-                    cosmology.comoving_distance, minimum * self.unit)
-            self._minimum['luminosity_distance'] = self._minimum['redshift']
-        try:
-            self._update_instance()
-        except (AttributeError, KeyError):
-            pass
+        self._set_limit(value=minimum, limit_dict=self._minimum)
 
     @property
     def maximum(self):
@@ -165,21 +144,44 @@ class Cosmological(Interped):
 
     @maximum.setter
     def maximum(self, maximum):
+        self._set_limit(value=maximum, limit_dict=self._maximum)
+
+    def _set_limit(self, value, limit_dict):
+        """
+        Set either the limits for redshift luminosity and comoving distances
+
+        Parameters
+        ----------
+        value: float
+            Limit value in current class' parameter
+        limit_dict: dict
+            The limit dictionary to modify in place
+        """
         cosmology = get_cosmology(self.cosmology)
-        self._maximum[self.name] = maximum
+        limit_dict[self.name] = value
         if self.name == 'redshift':
-            self._maximum['luminosity_distance'] = \
-                cosmology.luminosity_distance(maximum).value
-            self._maximum['comoving_distance'] = \
-                cosmology.comoving_distance(maximum).value
+            limit_dict['luminosity_distance'] = \
+                cosmology.luminosity_distance(value).value
+            limit_dict['comoving_distance'] = \
+                cosmology.comoving_distance(value).value
         elif self.name == 'luminosity_distance':
-            self._maximum['redshift'] = cosmo.z_at_value(
-                cosmology.luminosity_distance, maximum * self.unit)
-            self._maximum['comoving_distance'] = self._maximum['redshift']
+            if value == 0:
+                limit_dict['redshift'] = 0
+            else:
+                limit_dict['redshift'] = cosmo.z_at_value(
+                    cosmology.luminosity_distance, value * self.unit)
+            limit_dict['comoving_distance'] = (
+                cosmology.comoving_distance(limit_dict['redshift']).value
+            )
         elif self.name == 'comoving_distance':
-            self._maximum['redshift'] = cosmo.z_at_value(
-                cosmology.comoving_distance, maximum * self.unit)
-            self._maximum['luminosity_distance'] = self._maximum['redshift']
+            if value == 0:
+                limit_dict['redshift'] = 0
+            else:
+                limit_dict['redshift'] = cosmo.z_at_value(
+                    cosmology.comoving_distance, value * self.unit)
+            limit_dict['luminosity_distance'] = (
+                cosmology.luminosity_distance(limit_dict['redshift']).value
+            )
         try:
             self._update_instance()
         except (AttributeError, KeyError):
@@ -761,3 +763,362 @@ def secondary_mass_condition_function(reference_params, mass_1):
 ConditionalCosmological = conditional_prior_factory(Cosmological)
 ConditionalUniformComovingVolume = conditional_prior_factory(UniformComovingVolume)
 ConditionalUniformSourceFrame = conditional_prior_factory(UniformSourceFrame)
+
+
+class HealPixMapPriorDist(BaseJointPriorDist):
+    """
+    Class defining prior according to given HealPix Map, defaults to 2D in ra and dec but can be set to include
+    Distance as well. This only works with skymaps that include the 2D joint probability in ra/dec and that use the
+    normal LALInference type skymaps where each pixel has a DISTMU, DISTSIGMA, and DISTNORM defining the conditional
+    distance distribution along a given line of sight.
+
+    Parameters
+    ----------
+
+    hp_file : file path to .fits file
+        .fits file that containes the 2D or 3D Healpix Map
+    names : list (optional)
+        list of names of parameters included in the JointPriorDist, defaults to ['ra', 'dec']
+    bounds : dict or list (optional)
+        dictionary or list with given prior bounds. defaults to normal bounds on ra, dev and 0, inf for distance
+        if this is for a 3D map
+
+    Returns
+    -------
+
+    PriorDist : `bilby.gw.prior.HealPixMapPriorDist`
+        A JointPriorDist object to store the joint prior distribution according to passed healpix map
+    """
+    def __init__(self, hp_file, names=None, bounds=None, distance=False):
+        self.hp = self._check_imports()
+        self.hp_file = hp_file
+        if names is None:
+            names = ["ra", "dec"]
+        if bounds is None:
+            bounds = [[0, 2 * np.pi], [-np.pi / 2.0, np.pi / 2.0]]
+        elif isinstance(bounds, dict):
+            bs = [[] for _ in bounds.keys()]
+            for i, key in enumerate(bounds.keys()):
+                bs[i] = (bounds[key][0], bounds[key][1])
+            bounds = bs
+        if distance:
+            if len(names) == 2:
+                names.append("distance")
+            if len(bounds) == 2:
+                bounds.append([0, np.inf])
+            self.distance = True
+            self.prob, self.distmu, self.distsigma, self.distnorm = self.hp.read_map(
+                hp_file, verbose=False, field=range(4)
+            )
+        else:
+            self.distance = False
+            self.prob = self.hp.read_map(hp_file, verbose=False)
+
+        super(HealPixMapPriorDist, self).__init__(names=names, bounds=bounds)
+        self.distname = "hpmap"
+        self.npix = len(self.prob)
+        self.nside = self.hp.npix2nside(self.npix)
+        self.pixel_area = self.hp.nside2pixarea(self.nside)
+        self.pixel_length = self.pixel_area ** (1 / 2.0)
+        self.pix_xx = np.arange(self.npix)
+        self._all_interped = interp1d(x=self.pix_xx, y=self.prob, bounds_error=False, fill_value=0)
+        self.inverse_cdf = None
+        self.distance_pdf = None
+        self.distance_dist = None
+        self.distance_icdf = None
+        self._build_attributes()
+        name = self.names[-1]
+        if self.bounds[name][1] != np.inf and self.bounds[name][0] != -np.inf:
+            self.rs = np.linspace(self.bounds[name][0], self.bounds[name][1], 1000)
+        else:
+            self.rs = np.linspace(0, 5000, 1000)
+
+    def _build_attributes(self):
+        """
+        Method that builds the inverse cdf of the P(pixel) distribution for rescaling
+        """
+        yy = self._all_interped(self.pix_xx)
+        yy /= np.trapz(yy, self.pix_xx)
+        YY = cumtrapz(yy, self.pix_xx, initial=0)
+        YY[-1] = 1
+        self.inverse_cdf = interp1d(x=YY, y=self.pix_xx, bounds_error=True)
+
+    @staticmethod
+    def _check_imports():
+        """
+        Static method to check that healpy is installed on the machine running bibly
+        """
+        try:
+            import healpy
+        except Exception:
+            raise ImportError("Must have healpy installed on this machine to use HealPixMapPrior")
+        return healpy
+
+    def _rescale(self, samp, **kwargs):
+        """
+        Overwrites the _rescale method of BaseJoint Prior to rescale a single value from the unitcube onto
+        two values (ra, dec) or 3 (ra, dec, dist) if distance is included
+
+        Parameters
+        ----------
+        samp : float, int
+            must take in single value for pixel on unitcube to recale onto ra, dec (distance), for the map Prior
+        kwargs : dict
+            kwargs are all passed to _rescale() method
+
+        Returns
+        -------
+        rescaled_sample : array_like
+            sample to rescale onto the prior
+        """
+        if self.distance:
+            dist_samp = samp[:, -1]
+            samp = samp[:, 0]
+        else:
+            samp = samp[:, 0]
+        pix_rescale = self.inverse_cdf(samp)
+        sample = np.empty((len(pix_rescale), 2))
+        dist_samples = np.empty((len(pix_rescale)))
+        for i, val in enumerate(pix_rescale):
+            theta, ra = self.hp.pix2ang(self.nside, int(round(val)))
+            dec = 0.5 * np.pi - theta
+            sample[i, :] = self.draw_from_pixel(ra, dec, int(round(val)))
+            if self.distance:
+                self.update_distance(int(round(val)))
+                dist_samples[i] = self.distance_icdf(dist_samp[i])
+        if self.distance:
+            sample = np.row_stack([sample[:, 0], sample[:, 1], dist_samples])
+        return sample.reshape((-1, self.num_vars))
+
+    def update_distance(self, pix_idx):
+        """
+        Method to update the conditional distance distributions at given pixel used for distance handling in the
+        JointPrior Parameters. This function updates the current distance pdf, inverse_cdf, and sampler according to
+        given pixel or line of sight.
+
+        Parameters
+        ----------
+        pix_idx : int
+            pixel index value to create the distribtuion for
+
+        Returns
+        -------
+        None : None
+            just updates these functions at new pixel values
+        """
+        self.distance_pdf = lambda r: self.distnorm[pix_idx] * norm(
+            loc=self.distmu[pix_idx], scale=self.distsigma[pix_idx]
+        ).pdf(r)
+        pdfs = self.rs ** 2 * norm(loc=self.distmu[pix_idx], scale=self.distsigma[pix_idx]).pdf(self.rs)
+        cdfs = np.cumsum(pdfs) / np.sum(pdfs)
+
+        def sample_distance(n):
+            gaussian = norm(loc=self.distmu[pix_idx], scale=self.distsigma[pix_idx]).rvs(size=100 * n)
+            probs = self._check_norm(gaussian[gaussian > 0] ** 2)
+            ds = np.random.choice(gaussian[gaussian > 0], p=probs, size=n, replace=True)
+            return ds
+
+        self.distance_dist = sample_distance
+        self.distance_icdf = interp1d(cdfs, self.rs)
+
+    @staticmethod
+    def _check_norm(array):
+        """
+        static method to check if array is properlly normalized and if not to normalize it.
+
+        Parameters
+        ----------
+        array : array_like
+            input array we want to renormalize if not already normalized
+
+        Returns
+        -------
+        normed_array : array_like
+            returns input array normalized
+        """
+        norm = np.linalg.norm(array, ord=1)
+        if norm == 0:
+            norm = np.finfo(array.dtype).eps
+        return array / norm
+
+    def _sample(self, size, **kwargs):
+        """
+        Overwrites the _sample method of BaseJoint Prior. Picks a pixel value according to their probabilities, then
+        uniformly samples ra, and decs that are contained in chosen pixel. If the PriorDist includes distance it then
+        updates the distance distributions and will sample according to the conditional distance distribution along a
+        given line of sight
+
+        Parameters
+        ----------
+        size : int
+            number of samples we want to draw
+        kwargs : dict
+            kwargs are all passed to be used
+
+        Returns
+        -------
+        sample : array_like
+            sample of ra, and dec (and distance if 3D=True)
+        """
+        pixel_choices = np.arange(self.npix)
+        pixel_probs = self._check_norm(self.prob)
+        sample_pix = np.random.choice(pixel_choices, size=size, p=pixel_probs, replace=True)
+        sample = np.empty((size, self.num_vars))
+        for samp in range(size):
+            theta, ra = self.hp.pix2ang(self.nside, sample_pix[samp])
+            dec = 0.5 * np.pi - theta
+            if self.distance:
+                self.update_distance(sample_pix[samp])
+                dist = self.draw_distance(sample_pix[samp])
+                ra_dec = self.draw_from_pixel(ra, dec, sample_pix[samp])
+                sample[samp, :] = [ra_dec[0], ra_dec[1], dist]
+            else:
+                sample[samp, :] = self.draw_from_pixel(ra, dec, sample_pix[samp])
+        return sample.reshape((-1, self.num_vars))
+
+    def draw_distance(self, pix):
+        """
+        Method to recursively draw a distance value from the given set distance distribution and check that it is in
+        the bounds
+
+        Parameters
+        ----------
+
+        pix : int
+            integer for pixel to draw a distance from
+
+        Returns
+        -------
+        dist : float
+            sample drawn from the distance distribution at set pixel index
+        """
+        if self.distmu[pix] == np.inf or self.distmu[pix] <= 0:
+            return 0
+        dist = self.distance_dist(1)
+        name = self.names[-1]
+        if (dist > self.bounds[name][1]) | (dist < self.bounds[name][0]):
+            self.draw_distance(pix)
+        else:
+            return dist
+
+    def draw_from_pixel(self, ra, dec, pix):
+        """
+        Recursive function to uniformly draw ra, and dec values that are located in the given pixel
+
+        Parameters
+        ----------
+        ra : float, int
+            value drawn for rightascension
+        dec : float, int
+            value drawn for declination
+        pix : int
+            pixel index for given pixel we want to get ra, and dec from
+
+        Returns
+        -------
+        ra_dec : tuple
+            this returns a tuple of ra, and dec sampled uniformly that are in the pixel given
+        """
+        if not self.check_in_pixel(ra, dec, pix):
+            self.draw_from_pixel(ra, dec, pix)
+        return np.array(
+            [
+                np.random.uniform(ra - self.pixel_length, ra + self.pixel_length),
+                np.random.uniform(dec - self.pixel_length, dec + self.pixel_length),
+            ]
+        )
+
+    def check_in_pixel(self, ra, dec, pix):
+        """
+        Method that checks if given rightacension and declination values are within the given pixel index and the bounds
+
+        Parameters
+        ----------
+        ra : float, int
+            rightascension value to check
+        dec : float, int
+            declination value to check
+        pix : int
+            index for pixel we want to check in
+
+        Returns
+        -------
+        bool :
+            returns True if values inside pixel, False if not
+        """
+        for val, name in zip([ra, dec], self.names):
+            if (val < self.bounds[name][0]) or (val > self.bounds[name][1]):
+                return False
+        phi, theta = ra, 0.5 * np.pi - dec
+        pixel = self.hp.ang2pix(self.nside, theta, phi)
+        return pix == pixel
+
+    def _ln_prob(self, samp, lnprob, outbounds):
+        """
+        Overwrites the _lnprob method of BaseJoint Prior
+
+        Parameters
+        ----------
+        samp : array_like
+            samples of ra, dec to evaluate the lnprob at
+        lnprob : array_like
+            array of correct length we want to populate with lnprob values
+        outbounds : boolean array
+            boolean array that flags samples that are out of the given bounds
+
+        Returns
+        -------
+        lnprob : array_like
+            lnprob values at each sample
+        """
+        for i in range(samp.shape[0]):
+            if not outbounds[i]:
+                if self.distance:
+                    phi, dec, dist = samp[0]
+                else:
+                    phi, dec = samp[0]
+                theta = 0.5 * np.pi - dec
+                pixel = self.hp.ang2pix(self.nside, theta, phi)
+                lnprob[i] = np.log(self.prob[pixel] / self.pixel_area)
+                if self.distance:
+                    self.update_distance(pixel)
+                    lnprob[i] += np.log(self.distance_pdf(dist) * dist ** 2)
+        lnprob[outbounds] = -np.inf
+        return lnprob
+
+    def __eq__(self, other):
+        skip_keys = ["_all_interped", "inverse_cdf", "distance_pdf", "distance_dist", "distance_icdf"]
+        if self.__class__ != other.__class__:
+            return False
+        if sorted(self.__dict__.keys()) != sorted(other.__dict__.keys()):
+            return False
+        for key in self.__dict__:
+            if key in skip_keys:
+                continue
+            if key == "hp_file":
+                if self.__dict__[key] != other.__dict__[key]:
+                    return False
+            elif isinstance(self.__dict__[key], (np.ndarray, list)):
+                thisarr = np.asarray(self.__dict__[key])
+                otherarr = np.asarray(other.__dict__[key])
+                if thisarr.dtype == np.float and otherarr.dtype == np.float:
+                    fin1 = np.isfinite(np.asarray(self.__dict__[key]))
+                    fin2 = np.isfinite(np.asarray(other.__dict__[key]))
+                    if not np.array_equal(fin1, fin2):
+                        return False
+                    if not np.allclose(thisarr[fin1], otherarr[fin2], atol=1e-15):
+                        return False
+                else:
+                    if not np.array_equal(thisarr, otherarr):
+                        return False
+            else:
+                if not self.__dict__[key] == other.__dict__[key]:
+                    return False
+        return True
+
+
+class HealPixPrior(JointPrior):
+    def __init__(self, dist, name=None, latex_label=None, unit=None):
+        if not isinstance(dist, HealPixMapPriorDist):
+            raise JointPriorDistError("dist object must be instance of HealPixMapPriorDist")
+        super(HealPixPrior, self).__init__(dist=dist, name=name, latex_label=latex_label, unit=unit)

bilby/gw/prior_files/aligned_spins_waveform_tides_off.prior

+chirp_mass = Uniform(name='chirp_mass', minimum=3.56, maximum=3.68)
+mass_ratio = Uniform(name='mass_ratio', minimum=0.05, maximum=1)
+mass_1 = Constraint(name='mass_1', minimum=1.00, maximum=22.0) #optional
+mass_2 = Constraint(name='mass_2', minimum=1.00, maximum=2.95) #optional
+chi_1 = bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.5))
+chi_2 = bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.05))
+luminosity_distance = bilby.core.prior.PowerLaw(alpha=2, name='luminosity_distance', minimum=1, maximum=750, unit='Mpc')
+ra = Uniform(minimum=0, maximum=2 * np.pi, name='ra',  boundary='periodic')
+dec = Cosine(name='dec')
+cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')

bilby/gw/prior_files/aligned_spins_waveform_tides_on.prior

+chirp_mass = Uniform(name='chirp_mass', minimum=3.56, maximum=3.68)
+mass_ratio = Uniform(name='mass_ratio', minimum=0.05, maximum=1)
+mass_1 = Constraint(name='mass_1', minimum=1.00, maximum=22.0) #optional
+mass_2 = Constraint(name='mass_2', minimum=1.00, maximum=2.95) #optional
+chi_1 = bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.5))
+chi_2 = bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.05))
+luminosity_distance = bilby.core.prior.PowerLaw(alpha=2, name='luminosity_distance', minimum=1, maximum=750, unit='Mpc')
+ra = Uniform(minimum=0, maximum=2 * np.pi, name='ra',  boundary='periodic')
+dec = Cosine(name='dec')
+cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')
+lambda_1 = Uniform(name='lambda_2', minimum=0, maximum=5000)
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=5000)

bilby/gw/prior_files/precessing_spins_waveform_tides_off.prior

+chirp_mass = Uniform(name='chirp_mass', minimum=3.56, maximum=3.68)
+mass_ratio = Uniform(name='mass_ratio', minimum=0.05, maximum=1)
+mass_1 = Constraint(name='mass_1', minimum=1.00, maximum=22.0)
+mass_2 = Constraint(name='mass_2', minimum=1.00, maximum=22.0)
+a_1 = Uniform(name='a_1', minimum=0, maximum=0.89)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.05)
+tilt_1 = Sine(name='tilt_1')
+tilt_2 = Sine(name='tilt_2')
+phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic')
+phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic')
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1, maximum=750, unit='Mpc')
+ra = Uniform(minimum=0, maximum=2 * np.pi, name='ra',  boundary='periodic')
+dec = Cosine(name='dec')
+cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')

bilby/gw/prior_files/precessing_spins_waveform_tides_on.prior

+chirp_mass = Uniform(name='chirp_mass', minimum=3.56, maximum=3.68)
+mass_ratio = Uniform(name='mass_ratio', minimum=0.05, maximum=1)
+mass_1 = Constraint(name='mass_1', minimum=1.00, maximum=22.0)
+mass_2 = Constraint(name='mass_2', minimum=1.00, maximum=22.0)
+a_1 = Uniform(name='a_1', minimum=0, maximum=0.89)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.05)
+tilt_1 = Sine(name='tilt_1')
+tilt_2 = Sine(name='tilt_2')
+phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic')
+phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic')
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1, maximum=750, unit='Mpc')
+ra = Uniform(minimum=0, maximum=2 * np.pi, name='ra',  boundary='periodic')
+dec = Cosine(name='dec')
+cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')
+lambda_1 = Uniform(name='lambda_2', minimum=0, maximum=5000)
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=5000)

bilby/gw/result.py

@@ -9,7 +9,10 @@ from matplotlib import rcParams
 import numpy as np
 
 from ..core.result import Result as CoreResult
-from ..core.utils import infft, logger, check_directory_exists_and_if_not_mkdir
+from ..core.utils import (
+    infft, logger, check_directory_exists_and_if_not_mkdir,
+    latex_plot_format, safe_save_figure
+)
 from .utils import plot_spline_pos, spline_angle_xform, asd_from_freq_series
 from .detector import get_empty_interferometer, Interferometer
 
@@ -133,6 +136,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
             logger.info("No injection for detector {}".format(detector))
             return None
 
+    @latex_plot_format
     def plot_calibration_posterior(self, level=.9, format="png"):
         """ Plots the calibration amplitude and phase uncertainty.
         Adapted from the LALInference version in bayespputils
@@ -148,12 +152,6 @@ class CompactBinaryCoalescenceResult(CoreResult):
         """
         if format not in ["png", "pdf"]:
             raise ValueError("Format should be one of png or pdf")
-        _old_tex = rcParams["text.usetex"]
-        _old_serif = rcParams["font.serif"]
-        _old_family = rcParams["font.family"]
-        rcParams["text.usetex"] = True
-        rcParams["font.serif"] = "Computer Modern Roman"
-        rcParams["font.family"] = "Serif"
 
         fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(15, 15), dpi=500)
         posterior = self.posterior
@@ -214,17 +212,11 @@ class CompactBinaryCoalescenceResult(CoreResult):
 
         filename = os.path.join(outdir, self.label + '_calibration.' + format)
         fig.tight_layout()
-        try:
-            plt.savefig(filename, format=format, dpi=600, bbox_inches='tight')
-        except RuntimeError:
-            logger.debug(
-                "Failed to save waveform with tex labels turning off tex."
-            )
-            rcParams["text.usetex"] = False
-            plt.savefig(filename, format=format, dpi=600, bbox_inches='tight')
-        rcParams["text.usetex"] = _old_tex
-        rcParams["font.serif"] = _old_serif
-        rcParams["font.family"] = _old_family
+        safe_save_figure(
+            fig=fig, filename=filename,
+            format=format, dpi=600, bbox_inches='tight'
+        )
+        logger.debug("Calibration figure saved to {}".format(filename))
         plt.close()
 
     def plot_waveform_posterior(
@@ -268,6 +260,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
                 save=True, format=format, start_time=start_time,
                 end_time=end_time)
 
+    @latex_plot_format
     def plot_interferometer_waveform_posterior(
             self, interferometer, level=0.9, n_samples=None, save=True,
             format='png', start_time=None, end_time=None):
@@ -328,13 +321,6 @@ class CompactBinaryCoalescenceResult(CoreResult):
                     "HTML plotting requested, but plotly cannot be imported, "
                     "falling back to png format for waveform plot.")
                 format = "png"
-        else:
-            _old_tex = rcParams["text.usetex"]
-            _old_serif = rcParams["font.serif"]
-            _old_family = rcParams["font.family"]
-            rcParams["text.usetex"] = True
-            rcParams["font.serif"] = "Computer Modern Roman"
-            rcParams["font.family"] = "Serif"
 
         if isinstance(interferometer, str):
             interferometer = get_empty_interferometer(interferometer)
@@ -383,7 +369,6 @@ class CompactBinaryCoalescenceResult(CoreResult):
             len(frequency_idxs))
         )
         plot_times = interferometer.time_array[time_idxs]
-        # if format == "html":
         plot_times -= interferometer.strain_data.start_time
         start_time -= interferometer.strain_data.start_time
         end_time -= interferometer.strain_data.start_time
@@ -663,7 +648,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
             fig.update_xaxes(title_text=f_domain_x_label, type="log", row=1)
             fig.update_yaxes(title_text=f_domain_y_label, type="log", row=1)
             fig.update_xaxes(title_text=t_domain_x_label, type="linear", row=2)
-            fig.update_yaxes(title_text=t_domain_x_label, type="linear", row=2)
+            fig.update_yaxes(title_text=t_domain_y_label, type="linear", row=2)
         else:
             axs[0].set_xlim(interferometer.minimum_frequency,
                             interferometer.maximum_frequency)
@@ -684,19 +669,12 @@ class CompactBinaryCoalescenceResult(CoreResult):
                 plot(fig, filename=filename, include_mathjax='cdn', auto_open=False)
             else:
                 plt.tight_layout()
-                try:
-                    plt.savefig(filename, format=format, dpi=600)
-                except RuntimeError:
-                    logger.debug(
-                        "Failed to save waveform with tex labels turning off tex."
-                    )
-                    rcParams["text.usetex"] = False
-                    plt.savefig(filename, format=format, dpi=600)
+                safe_save_figure(
+                    fig=fig, filename=filename,
+                    format=format, dpi=600
+                )
                 plt.close()
-                rcParams["text.usetex"] = _old_tex
-                rcParams["font.serif"] = _old_serif
-                rcParams["font.family"] = _old_family
-            logger.debug("Figure saved to {}".format(filename))
+            logger.debug("Waveform figure saved to {}".format(filename))
         else:
             return fig
 
@@ -714,14 +692,14 @@ class CompactBinaryCoalescenceResult(CoreResult):
         Parameters
         ----------
         maxpts: int
-            Number of samples to use, if None all samples are used
+            Maximum number of samples to use, if None all samples are used
         trials: int
             Number of trials at each clustering number
         jobs: int
             Number of multiple threads
         enable_multiresolution: bool
             Generate a multiresolution HEALPix map (default: True)
-        objid: st
+        objid: str
             Event ID to store in FITS header
         instruments: str
             Name of detectors
@@ -766,16 +744,16 @@ class CompactBinaryCoalescenceResult(CoreResult):
         if load_pickle is False:
             try:
                 pts = data[['ra', 'dec', 'luminosity_distance']].values
-                cls = kde.Clustered2Plus1DSkyKDE
+                confidence_levels = kde.Clustered2Plus1DSkyKDE
                 distance = True
             except KeyError:
                 logger.warning("The results file does not contain luminosity_distance")
                 pts = data[['ra', 'dec']].values
-                cls = kde.Clustered2DSkyKDE
+                confidence_levels = kde.Clustered2DSkyKDE
                 distance = False
 
             logger.info('Initialising skymap class')
-            skypost = cls(pts, trials=trials, multiprocess=jobs)
+            skypost = confidence_levels(pts, trials=trials, jobs=jobs)
             logger.info('Pickling skymap to {}'.format(default_obj_filename))
             with open(default_obj_filename, 'wb') as out:
                 pickle.dump(skypost, out)
@@ -788,7 +766,8 @@ class CompactBinaryCoalescenceResult(CoreResult):
             logger.info('Reading from pickle {}'.format(obj_filename))
             with open(obj_filename, 'rb') as file:
                 skypost = pickle.load(file)
-            skypost.multiprocess = jobs
+            skypost.jobs = jobs
+            distance = isinstance(skypost, kde.Clustered2Plus1DSkyKDE)
 
         logger.info('Making skymap')
         hpmap = skypost.as_healpix()
@@ -844,12 +823,12 @@ class CompactBinaryCoalescenceResult(CoreResult):
             cb.set_label(r'prob. per deg$^2$')
 
         if contour is not None:
-            cls = 100 * postprocess.find_greedy_credible_levels(skymap)
-            cs = ax.contour_hpx(
-                (cls, 'ICRS'), nested=metadata['nest'],
+            confidence_levels = 100 * postprocess.find_greedy_credible_levels(skymap)
+            contours = ax.contour_hpx(
+                (confidence_levels, 'ICRS'), nested=metadata['nest'],
                 colors='k', linewidths=0.5, levels=contour)
             fmt = r'%g\%%' if rcParams['text.usetex'] else '%g%%'
-            plt.clabel(cs, fmt=fmt, fontsize=6, inline=True)
+            plt.clabel(contours, fmt=fmt, fontsize=6, inline=True)
 
         # Add continents.
         if geo:
@@ -875,7 +854,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
                 text.append('event ID: {}'.format(objid))
             if contour:
                 pp = np.round(contour).astype(int)
-                ii = np.round(np.searchsorted(np.sort(cls), contour) *
+                ii = np.round(np.searchsorted(np.sort(confidence_levels), contour) *
                               deg2perpix).astype(int)
                 for i, p in zip(ii, pp):
                     text.append(
@@ -884,15 +863,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
 
         filename = os.path.join(self.outdir, "{}_skymap.png".format(self.label))
         logger.info("Generating 2D projected skymap to {}".format(filename))
-        plt.savefig(filename, dpi=500)
-
-
-class CompactBinaryCoalesenceResult(CompactBinaryCoalescenceResult):
-
-    def __init__(self, **kwargs):
-        logger.warning('CompactBinaryCoalesenceResult is deprecated use '
-                       'CompactBinaryCoalescenceResult')
-        super(CompactBinaryCoalesenceResult, self).__init__(**kwargs)
+        safe_save_figure(fig=plt.gcf(), filename=filename, dpi=dpi)
 
 
 CBCResult = CompactBinaryCoalescenceResult

bilby/gw/source.py

@@ -16,8 +16,8 @@ try:
     import lal
     import lalsimulation as lalsim
 except ImportError:
-    logger.warning("You do not have lalsuite installed currently. You will"
-                   " not be able to use some of the prebuilt functions.")
+    logger.debug("You do not have lalsuite installed currently. You will"
+                 " not be able to use some of the prebuilt functions.")
 
 
 def lal_binary_black_hole(
@@ -301,6 +301,9 @@ def _base_lal_cbc_fd_waveform(
     pn_tidal_order = waveform_kwargs['pn_tidal_order']
     pn_phase_order = waveform_kwargs['pn_phase_order']
     pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
+    waveform_dictionary = waveform_kwargs.get(
+        'lal_waveform_dictionary', lal.CreateDict()
+    )
 
     approximant = lalsim_GetApproximantFromString(waveform_approximant)
 
@@ -327,7 +330,6 @@ def _base_lal_cbc_fd_waveform(
     longitude_ascending_nodes = 0.0
     mean_per_ano = 0.0
 
-    waveform_dictionary = lal.CreateDict()
     lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(
         waveform_dictionary, int(pn_spin_order))
     lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(
@@ -383,14 +385,26 @@ def _base_lal_cbc_fd_waveform(
     h_cross = np.zeros_like(frequency_array, dtype=np.complex)
 
     if len(hplus.data.data) > len(frequency_array):
-        raise ValueError("Waveform longer than frequency array")
-
-    h_plus[:len(hplus.data.data)] = hplus.data.data
-    h_cross[:len(hcross.data.data)] = hcross.data.data
+        logger.debug("LALsim waveform longer than bilby's `frequency_array`" +
+                     "({} vs {}), ".format(len(hplus.data.data), len(frequency_array)) +
+                     "probably because padded with zeros up to the next power of two length." +
+                     " Truncating lalsim array.")
+        h_plus = hplus.data.data[:len(h_plus)]
+        h_cross = hcross.data.data[:len(h_cross)]
+    else:
+        h_plus[:len(hplus.data.data)] = hplus.data.data
+        h_cross[:len(hcross.data.data)] = hcross.data.data
 
     h_plus *= frequency_bounds
     h_cross *= frequency_bounds
 
+    if wf_func == lalsim_SimInspiralFD:
+        dt = 1. / delta_frequency + (hplus.epoch.gpsSeconds + hplus.epoch.gpsNanoSeconds * 1e-9)
+        h_plus *= np.exp(
+            -1j * 2 * np.pi * dt * frequency_array)
+        h_cross *= np.exp(
+            -1j * 2 * np.pi * dt * frequency_array)
+
     return dict(plus=h_plus, cross=h_cross)
 
 

bilby/gw/utils.py

@@ -15,15 +15,15 @@ from ..core.utils import (ra_dec_to_theta_phi,
 try:
     from gwpy.timeseries import TimeSeries
 except ImportError:
-    logger.warning("You do not have gwpy installed currently. You will "
-                   " not be able to use some of the prebuilt functions.")
+    logger.debug("You do not have gwpy installed currently. You will "
+                 " not be able to use some of the prebuilt functions.")
 
 try:
     import lal
     import lalsimulation as lalsim
 except ImportError:
-    logger.warning("You do not have lalsuite installed currently. You will"
-                   " not be able to use some of the prebuilt functions.")
+    logger.debug("You do not have lalsuite installed currently. You will"
+                 " not be able to use some of the prebuilt functions.")
 
 
 def asd_from_freq_series(freq_data, df):

docs/bilby-output.txt

@@ -95,3 +95,11 @@ done via::
    result.save_posterior_samples()
 
 which will generate a file :code:`outdir/label_posterior.txt`.
+
+
+Visualising the results  
+--------------------
+Bilby also provides some useful built-in plotting tools. Some examples on how
+to visualise results using these tools (and  how to extend them) are shown in
+one of the tutorials at `visualising_the_results.ipynb  <https://git.ligo.org/lscsoft/bilby/-/blob/master/examples/tutorials/visualising_the_results.ipynb>`_.
+

docs/gw_prior.txt

+.. gw_prior:
+
+===================================
+Transient Graviatiaonal wave priors
+===================================
+
+A Cosmological GW prior, :code:`Cosmological`:
+
+.. autoclass:: bilby.gw.prior.Cosmological
+   :members:
+
+Uniform in Comoving Volume GW Prior (inherited from Cosmological) :code:`UniformComovingVolume`:
+
+.. autoclass:: bilby.gw.prior.UniformComovingVolume
+   :members:
+
+Uniform in Source Frame GW Prior :code:`UniformSourceFrame`:
+
+.. autoclass:: bilby.gw.prior.UniformSourceFrame
+   :members:
+
+Aligned Spine GW Prior :code:`AlignedSpin`:
+
+.. autoclass:: bilby.gw.prior.AlignedSpin
+   :members:
+
+HealPixMap JointPriorDist (See JointPriors in bilby.core.prior.joint) :code:`HealPixMapPriorDist`:
+
+.. autoclass:: bilby.gw.prior.HealPixMapPriorDist
+   :members:
+
+

docs/index.txt

@@ -20,6 +20,7 @@ Welcome to bilby's documentation!
    compact-binary-coalescence-parameter-estimation
    transient-gw-data
    gw_likelihood
+   gw_prior
    conversion
    gw_references
    writing-documentation

docs/prior.txt

@@ -178,8 +178,9 @@ First thing: define a function which generates z=x-y from x and y.
         -------
         dict: Dictionary with constraint parameter 'z' added.
         """
-        parameters['z'] = parameters['x'] - parameters['y']
-        return parameters
+        converted_parameters = parameters.copy()
+        converted_parameters['z'] = parameters['x'] - parameters['y']
+        return converted_parameters
 
 Create our prior:
 
@@ -216,3 +217,64 @@ We provide default conversion functions for the BBH and BNS PriorDicts.
 For BBHs this generates all the BBH mass parameters so constraints can be placed on any mass parameters.
 
 For BNSs it also generates the tidal deformability parameters.
+
+-----------------
+Prior Examples 
+-----------------
+
+Here we show some examples of prior files for different waveform families.
+To constrain a certain parameter to a fixed value, you just need:
+
+.. code:: python
+
+  parameter_name = <value>
+
+------
+
+To constrain the prior to a certain range , you can use:
+
+.. code:: python
+
+ parameter_name = Constraint(name='parameter_name', minimum=<value>, maximum=<value>)
+
+------
+
+Note that to activate the tidal effect you need to specify in your configuration
+file:
+
+.. code:: python
+
+ frequency_domain_source_model=lal_binary_neutron_star 
+
+------
+
+
+Aligned spins waveform with tides off
+==============
+
+.. literalinclude:: /../bilby/gw/prior_files/aligned_spins_waveform_tides_off.prior
+
+Aligned spins waveform with tides on
+==============
+
+.. literalinclude:: /../bilby/gw/prior_files/aligned_spins_waveform_tides_on.prior
+
+Precessing spins waveform with tides off
+==============
+
+.. literalinclude:: /../bilby/gw/prior_files/precessing_spins_waveform_tides_off.prior  
+
+Precessing spins waveform with tides on
+==============
+
+.. literalinclude:: /../bilby/gw/prior_files/precessing_spins_waveform_tides_on.prior  
+
+
+-----------------
+Priors using a Jupyter notebook 
+-----------------
+
+Bilby saves as output the prior volume sampled. You might also find useful to
+produce priors directly from a Jupyter notebook. You can have a look at one of
+the Bilby tutorials to check how you define and plot priors in a Jupyter notebook:
+`making_priors.ipynb <https://git.ligo.org/lscsoft/bilby/-/blob/master/examples/tutorials/making_priors.ipynb>`_.

examples/core_examples/15d_gaussian.py

-from scipy.stats import multivariate_normal
-
 import bilby
 import numpy as np
 
+from bilby.core.likelihood import AnalyticalMultidimensionalCovariantGaussian, \
+    AnalyticalMultidimensionalBimodalCovariantGaussian
+
 logger = bilby.core.utils.logger
 
 cov = [
@@ -58,68 +59,6 @@ cov = [
 dim = 15
 mean = np.zeros(dim)
 
-
-class AnalyticalMultidimensionalCovariantGaussian(bilby.Likelihood):
-    """
-        A multivariate Gaussian likelihood
-        with known analytic solution.
-
-        Parameters
-        ----------
-        mean: array_like
-            Array with the mean value of distribution
-        cov: array_like
-            The ndim*ndim covariance matrix
-        """
-
-    def __init__(self, mean, cov):
-        super(AnalyticalMultidimensionalCovariantGaussian, self).__init__(parameters=dict())
-        self.cov = np.array(cov)
-        self.mean = np.array(mean)
-        self.sigma = np.sqrt(np.diag(self.cov))
-        self.pdf = multivariate_normal(mean=self.mean, cov=self.cov)
-
-    @property
-    def dim(self):
-        return len(self.cov[0])
-
-    def log_likelihood(self):
-        x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)])
-        return self.pdf.logpdf(x)
-
-
-class AnalyticalMultidimensionalBimodalCovariantGaussian(bilby.Likelihood):
-    """
-        A multivariate Gaussian likelihood
-        with known analytic solution.
-
-        Parameters
-        ----------
-        mean_1: array_like
-            Array with the mean value of the first mode
-        mean_2: array_like
-            Array with the mean value of the second mode
-        cov: array_like
-        """
-
-    def __init__(self, mean_1, mean_2, cov):
-        super(AnalyticalMultidimensionalBimodalCovariantGaussian, self).__init__(parameters=dict())
-        self.cov = np.array(cov)
-        self.mean_1 = np.array(mean_1)
-        self.mean_2 = np.array(mean_2)
-        self.sigma = np.sqrt(np.diag(self.cov))
-        self.pdf_1 = multivariate_normal(mean=self.mean_1, cov=self.cov)
-        self.pdf_2 = multivariate_normal(mean=self.mean_2, cov=self.cov)
-
-    @property
-    def dim(self):
-        return len(self.cov[0])
-
-    def log_likelihood(self):
-        x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)])
-        return -np.log(2) + np.logaddexp(self.pdf_1.logpdf(x), self.pdf_2.logpdf(x))
-
-
 label = "multidim_gaussian_unimodal"
 outdir = "outdir"
 

examples/core_examples/logo/sample_logo.py

@@ -8,14 +8,14 @@ from skimage import io
 class Likelihood(bilby.Likelihood):
     def __init__(self, interp):
         self.interp = interp
-        self.parameters = dict(x=None, y=None)
+        super().__init__(parameters=dict(x=None, y=None))
 
     def log_likelihood(self):
         return -1 / (self.interp(self.parameters['x'], self.parameters['y'])[0])
 
 
 for letter in ['B', 'I', 'L', 'Y']:
-    img = 1 - io.imread('{}.png'.format(letter), as_grey=True)[::-1, :]
+    img = 1 - io.imread('{}.png'.format(letter), as_gray=True)[::-1, :]
     x = np.arange(img.shape[0])
     y = np.arange(img.shape[1])
     interp = si.interpolate.interp2d(x, y, img.T)

examples/tutorials/compare_samplers.ipynb

 %% Cell type:markdown id: tags:
 
 # Compare samplers
 
 In this notebook, we'll compare the different samplers implemented in `bilby`. As of this version, we don't compare the outputs, only how to run them and the timings for their default setup.
 
 ## Setup
 
 %% Cell type:code id: tags:
 
 ``` python
 import numpy as np
 import pylab as plt
 
 %load_ext autoreload
 %autoreload 2
 
 import bilby
 
 bilby.utils.setup_logger()
 
 time_duration = 1.         # set the signal duration (seconds)
 sampling_frequency = 4096. # set the data sampling frequency (Hz)
 
 injection_parameters = dict(
-mass_1=36.,                          # source frame (non-redshifted) primary mass (solar masses)
-mass_2=29.,                          # source frame (non-redshifted) secondary mass (solar masses)
+mass_1=36.,                          # detector frame (redshifted) primary mass (solar masses)
+mass_2=29.,                          # detector frame (redshifted) secondary mass (solar masses)
 a_1=0,                               # primary dimensionless spin magnitude
 a_2=0,                               # secondary dimensionless spin magnitude
 tilt_1=0,                            # polar angle between primary spin and the orbital angular momentum (radians)
 tilt_2=0,                            # polar angle between secondary spin and the orbital angular momentum
 phi_12=0,                            # azimuthal angle between primary and secondary spin (radians)
 phi_jl=0,                            # azimuthal angle between total angular momentum and orbital angular momentum (radians)
 luminosity_distance=100.,            # luminosity distance to source (Mpc)
 iota=0.4,                            # inclination angle between line of sight and orbital angular momentum (radians)
 phase=1.3,                           # phase (radians)
 waveform_approximant='IMRPhenomPv2', # waveform approximant name
 reference_frequency=50.,             # gravitational waveform reference frequency (Hz)
 ra=1.375,                            # source right ascension (radians)
 dec=-1.2108,                         # source declination (radians)
 geocent_time=1126259642.413,         # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds)
 psi=2.659                            # gravitational wave polarisation angle
 )
 
 
 # initialise the waveform generator
 waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
     sampling_frequency=sampling_frequency,
     duration=time_duration,
     frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
     parameters=injection_parameters)
 
 # generate a frequency-domain waveform
 hf_signal = waveform_generator.frequency_domain_strain()
 
 # initialise a single interferometer representing LIGO Hanford
 H1 = bilby.gw.detector.get_empty_interferometer('H1')
 # set the strain data at the interferometer
 H1.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, duration=time_duration)
 # inject the gravitational wave signal into the interferometer model
 H1.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters)
 
 IFOs = [H1]
 
 # compute the likelihood on each of the signal parameters
 likelihood = bilby.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Prior
 
 For this test, we will simply search of the sky position, setting the other parameters to their simulated values.
 
 %% Cell type:code id: tags:
 
 ``` python
 # set the priors on each of the injection parameters to be a delta function at their given value
 priors = bilby.gw.prior.BBHPriorDict()
 for key in injection_parameters.keys():
     priors[key] = injection_parameters[key]
 
 # now reset the priors on the sky position co-ordinates in order to conduct a sky position search
 priors['ra'] = bilby.prior.Uniform(0, 2*np.pi, 'ra')
 priors['dec'] = bilby.prior.Cosine(name='dec', minimum=-np.pi/2, maximum=np.pi/2)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## PyMultinest
 
 %% Cell type:code id: tags:
 
 ``` python
 %%time
 result = bilby.core.sampler.run_sampler(
     likelihood, priors=priors, sampler='pymultinest', label='pymultinest',
     npoints=200, verbose=False, resume=False)
 fig = result.plot_corner(save=False)
 # show the corner plot
 plt.show()
 print(result)
 ```
 
 %% Cell type:markdown id: tags:
 
 # dynesty
 
 %% Cell type:code id: tags:
 
 ``` python
 %%time
 result = bilby.core.sampler.run_sampler(
     likelihood, priors=priors, sampler='dynesty', label='dynesty',
     bound='multi', sample='rwalk', npoints=200, walks=1, verbose=False,
     update_interval=100)
 fig = result.plot_corner(save=False)
 # show the corner plot
 plt.show()
 print(result)
 ```
 
 %% Cell type:markdown id: tags:
 
 # Dynamic Nested Sampling (Dynesty)
 
 See [the dynesty docs](http://dynesty.readthedocs.io/en/latest/dynamic.html#). Essentially, this methods improves the posterior estimation over that of standard nested sampling.
 
 %% Cell type:code id: tags:
 
 ``` python
 %%time
 result = bilby.core.sampler.run_sampler(
     likelihood, priors=priors, sampler='dynesty', label='dynesty_dynamic',
     bound='multi', nlive=250, sample='unif', verbose=True,
     update_interval=100, dynamic=True)
 fig = result.plot_corner(save=False)
 # show the corner plot
 plt.show()
 print(result)
 ```
 
 %% Cell type:markdown id: tags:
 
 # ptemcee
 
 %% Cell type:code id: tags:
 
 ``` python
 %%time
 result = bilby.core.sampler.run_sampler(
     likelihood, priors=priors, sampler='ptemcee', label='ptemcee',
     nwalkers=100, nsteps=200, nburn=100, ntemps=2,
     tqdm='tqdm_notebook')
 fig = result.plot_corner(save=False)
 # show the corner plot
 plt.show()
 print(result)
 ```

setup.cfg

 [flake8]
 exclude = .git,docs,build,dist,test,*__init__.py
 max-line-length = 120
-ignore = E129 W504 W605
+ignore = E129 W503 W504 W605 E203
 
 [tool:pytest]
 addopts =
@@ -9,6 +9,8 @@ addopts =
     --ignore test/gw_example_test.py
     --ignore test/example_test.py
     --ignore test/sample_from_the_prior_test.py
+    --ignore test/gw_plot_test.py
+    --ignore test/sampler_test.py
 
 [metadata]
 license_file = LICENSE.md

setup.py

@@ -64,7 +64,7 @@ def readfile(filename):
     return filecontents
 
 
-VERSION = '0.6.4'
+VERSION = '0.6.7'
 version_file = write_version_file(VERSION)
 long_description = get_long_description()
 

test/gw_plot_test.py

+import os
+import shutil
+import unittest
+
+import pandas as pd
+
+import bilby
+
+
+class TestCBCResult(unittest.TestCase):
+
+    def setUp(self):
+        bilby.utils.command_line_args.bilby_test_mode = False
+        priors = bilby.gw.prior.BBHPriorDict()
+        priors['geocent_time'] = 2
+        injection_parameters = priors.sample()
+        self.meta_data = dict(
+            likelihood=dict(
+                phase_marginalization=True,
+                distance_marginalization=False,
+                time_marginalization=True,
+                frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+                waveform_arguments=dict(
+                    reference_frequency=20.0,
+                    waveform_approximant='IMRPhenomPv2'),
+                interferometers=dict(
+                    H1=dict(optimal_SNR=1, parameters=injection_parameters),
+                    L1=dict(optimal_SNR=1, parameters=injection_parameters)),
+                sampling_frequency=4096,
+                duration=4,
+                start_time=0,
+                waveform_generator_class=bilby.gw.waveform_generator.WaveformGenerator,
+                parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+            )
+        )
+        self.result = bilby.gw.result.CBCResult(
+            label='label', outdir='outdir', sampler='nestle',
+            search_parameter_keys=list(priors.keys()), fixed_parameter_keys=list(),
+            priors=priors, sampler_kwargs=dict(test='test', func=lambda x: x),
+            injection_parameters=injection_parameters,
+            meta_data=self.meta_data,
+            posterior=pd.DataFrame(priors.sample(100))
+        )
+        if not os.path.isdir(self.result.outdir):
+            os.mkdir(self.result.outdir)
+        pass
+
+    def tearDown(self):
+        bilby.utils.command_line_args.bilby_test_mode = True
+        try:
+            shutil.rmtree(self.result.outdir)
+        except OSError:
+            pass
+        del self.result
+        pass
+
+    def test_calibration_plot(self):
+        calibration_prior = bilby.gw.prior.CalibrationPriorDict.constant_uncertainty_spline(
+            amplitude_sigma=0.1,
+            phase_sigma=0.1,
+            minimum_frequency=20,
+            maximum_frequency=2048,
+            label="recalib_H1_",
+            n_nodes=5,
+        )
+        calibration_filename = f"{self.result.outdir}/{self.result.label}_calibration.png"
+        for key in calibration_prior:
+            self.result.posterior[key] = calibration_prior[key].sample(100)
+        self.result.plot_calibration_posterior()
+        self.assertTrue(os.path.exists(calibration_filename))
+
+    def test_calibration_plot_returns_none_with_no_calibration_parameters(self):
+        self.assertIsNone(self.result.plot_calibration_posterior())
+        calibration_filename = f"{self.result.outdir}/{self.result.label}_calibration.png"
+        self.assertFalse(os.path.exists(calibration_filename))
+
+    def test_calibration_pdf_plot(self):
+        calibration_prior = bilby.gw.prior.CalibrationPriorDict.constant_uncertainty_spline(
+            amplitude_sigma=0.1,
+            phase_sigma=0.1,
+            minimum_frequency=20,
+            maximum_frequency=2048,
+            label="recalib_H1_",
+            n_nodes=5,
+        )
+        calibration_filename = f"{self.result.outdir}/{self.result.label}_calibration.pdf"
+        for key in calibration_prior:
+            self.result.posterior[key] = calibration_prior[key].sample(100)
+        self.result.plot_calibration_posterior(format="pdf")
+        self.assertTrue(os.path.exists(calibration_filename))
+
+    def test_calibration_invalid_format_raises_error(self):
+        with self.assertRaises(ValueError):
+            self.result.plot_calibration_posterior(format="bilby")
+
+    def test_waveform_plotting_png(self):
+        self.result.plot_waveform_posterior(n_samples=200)
+        for ifo in self.result.interferometers:
+            self.assertTrue(os.path.exists(
+                f"{self.result.outdir}/{self.result.label}_{ifo}_waveform.png")
+            )
+
+    def test_plot_skymap_meta_data(self):
+        from ligo.skymap import io
+        expected_keys = {
+            "HISTORY", "build_date", "creator", "distmean", "diststd",
+            "gps_creation_time", "gps_time", "nest", "objid", "origin",
+            "vcs_revision", "vcs_version", "instruments"
+        }
+        self.result.plot_skymap(
+            maxpts=50, geo=False, objid="test", instruments="H1L1"
+        )
+        fits_filename = f"{self.result.outdir}/{self.result.label}_skymap.fits"
+        skymap_filename = f"{self.result.outdir}/{self.result.label}_skymap.png"
+        pickle_filename = f"{self.result.outdir}/{self.result.label}_skypost.obj"
+        hpmap, meta = io.read_sky_map(fits_filename)
+        self.assertEqual(expected_keys, set(meta.keys()))
+        self.assertTrue(os.path.exists(skymap_filename))
+        self.assertTrue(os.path.exists(pickle_filename))
+        self.result.plot_skymap(
+            maxpts=50, geo=False, objid="test", instruments="H1L1",
+            load_pickle=True, colorbar=True
+        )
+
+
+if __name__ == '__main__':
+    unittest.main()

test/gw_result_test.py

-from __future__ import absolute_import, division
+import os
+import shutil
+import unittest
+
+import pandas as pd
 
 import bilby
-import unittest
-import shutil
 
 
 class TestCBCResult(unittest.TestCase):
 
     def setUp(self):
         bilby.utils.command_line_args.bilby_test_mode = False
-        priors = bilby.prior.PriorDict(dict(
-            x=bilby.prior.Uniform(0, 1, 'x', latex_label='$x$', unit='s'),
-            y=bilby.prior.Uniform(0, 1, 'y', latex_label='$y$', unit='m'),
-            c=1,
-            d=2))
+        priors = bilby.gw.prior.BBHPriorDict()
+        priors['geocent_time'] = 2
+        injection_parameters = priors.sample()
         self.meta_data = dict(
             likelihood=dict(
-                phase_marginalization=True, distance_marginalization=False,
+                phase_marginalization=True,
+                distance_marginalization=False,
                 time_marginalization=True,
                 frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
                 waveform_arguments=dict(
                     reference_frequency=20.0,
                     waveform_approximant='IMRPhenomPv2'),
                 interferometers=dict(
-                    H1=dict(optimal_SNR=1, parameters=dict(x=0.1, y=0.3)),
-                    L1=dict(optimal_SNR=1, parameters=dict(x=0.1, y=0.3)))))
+                    H1=dict(optimal_SNR=1, parameters=injection_parameters),
+                    L1=dict(optimal_SNR=1, parameters=injection_parameters)),
+                sampling_frequency=4096,
+                duration=4,
+                start_time=0,
+                waveform_generator_class=bilby.gw.waveform_generator.WaveformGenerator,
+                parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+            )
+        )
         self.result = bilby.gw.result.CBCResult(
             label='label', outdir='outdir', sampler='nestle',
-            search_parameter_keys=['x', 'y'], fixed_parameter_keys=['c', 'd'],
+            search_parameter_keys=list(priors.keys()), fixed_parameter_keys=list(),
             priors=priors, sampler_kwargs=dict(test='test', func=lambda x: x),
-            injection_parameters=dict(x=0.5, y=0.5),
-            meta_data=self.meta_data)
+            injection_parameters=injection_parameters,
+            meta_data=self.meta_data,
+            posterior=pd.DataFrame(priors.sample(100))
+        )
+        if not os.path.isdir(self.result.outdir):
+            os.mkdir(self.result.outdir)
         pass
 
     def tearDown(self):
@@ -82,6 +94,36 @@ class TestCBCResult(unittest.TestCase):
         with self.assertRaises(AttributeError):
             self.result.reference_frequency
 
+    def test_sampling_frequency(self):
+        self.assertEqual(
+            self.result.sampling_frequency,
+            self.meta_data['likelihood']['sampling_frequency'])
+
+    def test_sampling_frequency_unset(self):
+        self.result.meta_data['likelihood'].pop('sampling_frequency')
+        with self.assertRaises(AttributeError):
+            self.result.sampling_frequency
+
+    def test_duration(self):
+        self.assertEqual(
+            self.result.duration,
+            self.meta_data['likelihood']['duration'])
+
+    def test_duration_unset(self):
+        self.result.meta_data['likelihood'].pop('duration')
+        with self.assertRaises(AttributeError):
+            self.result.duration
+
+    def test_start_time(self):
+        self.assertEqual(
+            self.result.start_time,
+            self.meta_data['likelihood']['start_time'])
+
+    def test_start_time_unset(self):
+        self.result.meta_data['likelihood'].pop('start_time')
+        with self.assertRaises(AttributeError):
+            self.result.start_time
+
     def test_waveform_approximant(self):
         self.assertEqual(
             self.result.waveform_approximant,
@@ -107,6 +149,26 @@ class TestCBCResult(unittest.TestCase):
         with self.assertRaises(AttributeError):
             self.result.frequency_domain_source_model
 
+    def test_parameter_conversion(self):
+        self.assertEqual(
+            self.result.parameter_conversion,
+            self.meta_data['likelihood']['parameter_conversion'])
+
+    def test_parameter_conversion_unset(self):
+        self.result.meta_data['likelihood'].pop('parameter_conversion')
+        with self.assertRaises(AttributeError):
+            self.result.parameter_conversion
+
+    def test_waveform_generator_class(self):
+        self.assertEqual(
+            self.result.waveform_generator_class,
+            self.meta_data['likelihood']['waveform_generator_class'])
+
+    def test_waveform_generator_class_unset(self):
+        self.result.meta_data['likelihood'].pop('waveform_generator_class')
+        with self.assertRaises(AttributeError):
+            self.result.waveform_generator_class
+
     def test_interferometer_names(self):
         self.assertEqual(
             self.result.interferometers,