.gitlab-ci.yml

@@ -138,6 +138,12 @@ plotting:
 
     - pytest test/gw/plot_test.py
 
+authors:
+  stage: test
+  image: bilbydev/bilby-test-suite-python37
+  script:
+    - python test/check_author_list.py
+
 pages:
   stage: deploy
   dependencies:

AUTHORS.md

+# Authors
+
+This file lists all the authors in first-name alphabetical order who have
+contributed (either by code contribution or indirectly). If your name is not
+listed here, please contact anyone on this list and raise your concern.
+
+Abhirup Ghosh
+Aditya Vijaykumar
+Andrew Kim
+Andrew Miller
+Antoni Ramos-Buades
+Avi Vajpeyi
+Bruce Edelman
+Carl-Johan Haster
+Cecilio Garcia-Quiros
+Charlie Hoy
+Christopher Berry
+Christos Karathanasis
+Colm Talbot
+Daniel Williams
+David Keitel
+Duncan Macleod
+Eric Thrane
+Ethan Payne
+Francisco Javier Hernandez
+Gregory Ashton
+Hector Estelles
+Ignacio Magaña Hernandez
+Isobel Marguarethe Romero-Shaw
+Jade Powell
+James A Clark
+John Veitch
+Katerina Chatziioannou
+Kaylee de Soto
+Khun Sang Phukon
+Kshipraa Athar
+Liting Xiao
+Maite Mateu-Lucena
+Marc Arene
+Marcus Edward Lower
+Margaret Millhouse
+Marta Colleoni
+Matthew Carney
+Matthew David Pitkin
+Michael Puerrer
+Michael Williams
+Monica Rizzo
+Moritz Huebner
+Nicola De Lillo
+Nikhil Sarin
+Nirban Bose
+Paul Easter
+Paul Lasky
+Philip Relton
+Rhys Green
+Roberto Cotesta
+Rory Smith
+S. H. Oh
+Sacha Husa
+Scott Coughlin
+Serguei Ossokine
+Shanika Galaudage
+Sharan Banagiri
+Shichao Wu
+Simon Stevenson
+Soichiro Morisaki
+Sumeet Kulkarni
+Sylvia Biscoveanu
+Tathagata Ghosh
+Virginia d'Emilio
+Vivien Raymond

CHANGELOG.md

 # All notable changes will be documented in this file
 
+## [1.0.3] 2020-10-23
+
+Version 1.0.3 release of bilby
+
+### Added
+- SlabSpikePrior and examples (!857)
+- Authors file (!885)
+- CDF function to conditional priors (!882)
+- Waveform plot in visualising_the_results.ipynb (!817)
+- Addition of dnest4 sampler (!849, !883)
+- Loaded modules added to meta-data (!881)
+
+### Changes
+- Constraint to Uniform priors in ROQ tutorial (!884)
+- Fix to CDF and PDF for SymmetricLogUniform prior (!876)
+- Fix bug in evidence combination (!880)
+- Typo fixes (!878, !887, !879)
+- Minor bug fixes (!888)
+
 ## [1.0.2] 2020-09-14
 
 Version 1.0.2 release of bilby
@@ -8,7 +27,7 @@ Version 1.0.2 release of bilby
 - Template for the docker files (!783)
 - New delta_phase parameter (!850)
 - Normalization factor to time-domain waveform plot (!867)
-- JSON encoding for int and float types (!866) 
+- JSON encoding for int and float types (!866)
 - Various minor formatting additions (!870)
 
 ### Changes

bilby/core/prior/__init__.py

@@ -4,3 +4,4 @@ from .conditional import *
 from .dict import *
 from .interpolated import *
 from .joint import *
+from .slabspike import *

bilby/core/prior/analytical.py

@@ -337,7 +337,8 @@ class SymmetricLogUniform(Prior):
         -------
         float: Prior probability of val
         """
-        return (np.nan_to_num(0.5 / np.abs(val) / np.log(self.maximum / self.minimum)) *
+        val = np.abs(val)
+        return (np.nan_to_num(0.5 / val / np.log(self.maximum / self.minimum)) *
                 self.is_in_prior_range(val))
 
     def ln_prob(self, val):
@@ -354,6 +355,18 @@ class SymmetricLogUniform(Prior):
         """
         return np.nan_to_num(- np.log(2 * np.abs(val)) - np.log(np.log(self.maximum / self.minimum)))
 
+    def cdf(self, val):
+        val = np.atleast_1d(val)
+        norm = 0.5 / np.log(self.maximum / self.minimum)
+        cdf = np.zeros((len(val)))
+        lower_indices = np.where(np.logical_and(-self.maximum <= val, val <= -self.minimum))[0]
+        upper_indices = np.where(np.logical_and(self.minimum <= val, val <= self.maximum))[0]
+        cdf[lower_indices] = -norm * np.log(-val[lower_indices] / self.maximum)
+        cdf[np.where(np.logical_and(-self.minimum < val, val < self.minimum))] = 0.5
+        cdf[upper_indices] = 0.5 + norm * np.log(val[upper_indices] / self.minimum)
+        cdf[np.where(self.maximum < val)] = 1
+        return cdf
+
 
 class Cosine(Prior):
 
@@ -605,6 +618,7 @@ class TruncatedGaussian(Prior):
             / self.sigma / self.normalisation * self.is_in_prior_range(val)
 
     def cdf(self, val):
+        val = np.atleast_1d(val)
         _cdf = (erf((val - self.mu) / 2 ** 0.5 / self.sigma) - erf(
             (self.minimum - self.mu) / 2 ** 0.5 / self.sigma)) / 2 / self.normalisation
         _cdf[val > self.maximum] = 1

bilby/core/prior/conditional.py

@@ -110,9 +110,41 @@ def conditional_prior_factory(prior_class):
             return super(ConditionalPrior, self).prob(val)
 
         def ln_prob(self, val, **required_variables):
+            """Return the natural log prior probability of val.
+
+            Parameters
+            ----------
+            val: Union[float, int, array_like]
+                See superclass
+            required_variables:
+                Any required variables that this prior depends on
+
+
+            Returns
+            -------
+            float: Natural log prior probability of val
+            """
             self.update_conditions(**required_variables)
             return super(ConditionalPrior, self).ln_prob(val)
 
+        def cdf(self, val, **required_variables):
+            """Return the cdf of val.
+
+            Parameters
+            ----------
+            val: Union[float, int, array_like]
+                See superclass
+            required_variables:
+                Any required variables that this prior depends on
+
+
+            Returns
+            -------
+            float: CDF of val
+            """
+            self.update_conditions(**required_variables)
+            return super(ConditionalPrior, self).cdf(val)
+
         def update_conditions(self, **required_variables):
             """
             This method updates the conditional parameters (depending on the parent class

bilby/core/prior/slabspike.py

+import numpy as np
+
+from bilby.core.prior.base import Prior
+from bilby.core.utils import logger
+
+
+class SlabSpikePrior(Prior):
+
+    def __init__(self, slab, spike_location=None, spike_height=0):
+        """'Slab-and-spike' prior, see e.g. https://arxiv.org/abs/1812.07259
+        This prior is composed of a `slab`, i.e. any common prior distribution,
+        and a Dirac spike at a fixed location. This can effectively be used
+        to emulate sampling in the number of dimensions (similar to reversible-
+        jump MCMC).
+
+        `SymmetricLogUniform` and `FermiDirac` are currently not supported.
+
+        Parameters
+        ----------
+        slab: Prior
+            Any instance of a bilby prior class. All general prior attributes
+            from the slab are copied into the SlabSpikePrior.
+            Note that this hasn't been tested for conditional priors.
+        spike_location: float, optional
+            Location of the Dirac spike. Must be between minimum and maximum
+            of the slab. Defaults to the minimum of the slab
+        spike_height: float, optional
+            Relative weight of the spike compared to the slab. Must be
+            between 0 and 1. Defaults to 0, i.e. the prior is just the slab.
+
+        """
+        self.slab = slab
+        super().__init__(name=self.slab.name, latex_label=self.slab.latex_label, unit=self.slab.unit,
+                         minimum=self.slab.minimum, maximum=self.slab.maximum,
+                         check_range_nonzero=self.slab.check_range_nonzero, boundary=self.slab.boundary)
+        self.spike_location = spike_location
+        self.spike_height = spike_height
+        try:
+            self.inverse_cdf_below_spike = self._find_inverse_cdf_fraction_before_spike()
+        except Exception as e:
+            logger.warning("Disregard the following warning when running tests:\n {}".format(e))
+
+    @property
+    def spike_location(self):
+        return self._spike_loc
+
+    @spike_location.setter
+    def spike_location(self, spike_loc):
+        if spike_loc is None:
+            spike_loc = self.minimum
+        if not self.minimum <= spike_loc <= self.maximum:
+            raise ValueError("Spike location {} not within prior domain ".format(spike_loc))
+        self._spike_loc = spike_loc
+
+    @property
+    def spike_height(self):
+        return self._spike_height
+
+    @spike_height.setter
+    def spike_height(self, spike_height):
+        if 0 <= spike_height <= 1:
+            self._spike_height = spike_height
+        else:
+            raise ValueError("Spike height must be between 0 and 1, but is {}".format(spike_height))
+
+    @property
+    def slab_fraction(self):
+        """ Relative prior weight of the slab. """
+        return 1 - self.spike_height
+
+    def _find_inverse_cdf_fraction_before_spike(self):
+        return float(self.slab.cdf(self.spike_location)) * self.slab_fraction
+
+    def rescale(self, val):
+        """
+        'Rescale' a sample from the unit line element to the prior.
+
+        Parameters
+        ----------
+        val: Union[float, int, array_like]
+            A random number between 0 and 1
+
+        Returns
+        -------
+        array_like: Associated prior value with input value.
+        """
+        val = np.atleast_1d(val)
+
+        lower_indices = np.where(val < self.inverse_cdf_below_spike)[0]
+        intermediate_indices = np.where(np.logical_and(
+            self.inverse_cdf_below_spike <= val,
+            val <= self.inverse_cdf_below_spike + self.spike_height))[0]
+        higher_indices = np.where(val > self.inverse_cdf_below_spike + self.spike_height)[0]
+
+        res = np.zeros(len(val))
+        res[lower_indices] = self._contracted_rescale(val[lower_indices])
+        res[intermediate_indices] = self.spike_location
+        res[higher_indices] = self._contracted_rescale(val[higher_indices] - self.spike_height)
+        return res
+
+    def _contracted_rescale(self, val):
+        """
+        Contracted version of the rescale function that implements the `rescale` function
+        on the pure slab part of the prior.
+
+        Parameters
+        ----------
+        val: Union[float, int, array_like]
+            A random number between 0 and self.slab_fraction
+
+        Returns
+        -------
+        array_like: Associated prior value with input value.
+        """
+        return self.slab.rescale(val / self.slab_fraction)
+
+    def prob(self, val):
+        """Return the prior probability of val.
+        Returns np.inf for the spike location
+
+        Parameters
+        ----------
+        val: Union[float, int, array_like]
+
+        Returns
+        -------
+        array_like: Prior probability of val
+        """
+        res = self.slab.prob(val) * self.slab_fraction
+        res = np.atleast_1d(res)
+        res[np.where(val == self.spike_location)] = np.inf
+        return res
+
+    def ln_prob(self, val):
+        """Return the Log prior probability of val.
+        Returns np.inf for the spike location
+
+        Parameters
+        ----------
+        val: Union[float, int, array_like]
+
+        Returns
+        -------
+        array_like: Prior probability of val
+        """
+        res = self.slab.ln_prob(val) + np.log(self.slab_fraction)
+        res = np.atleast_1d(res)
+        res[np.where(val == self.spike_location)] = np.inf
+        return res
+
+    def cdf(self, val):
+        """ Return the CDF of the prior.
+        This calls to the slab CDF and adds a discrete step
+        at the spike location.
+
+        Parameters
+        ----------
+        val: Union[float, int, array_like]
+
+        Returns
+        -------
+        array_like: CDF value of val
+
+        """
+        res = self.slab.cdf(val) * self.slab_fraction
+        res = np.atleast_1d(res)
+        indices_above_spike = np.where(val > self.spike_location)[0]
+        res[indices_above_spike] += self.spike_height
+        return res

bilby/core/result.py

@@ -1644,24 +1644,26 @@ class ResultList(list):
             The result object with the combined evidences.
         """
         self.check_nested_samples()
-        if result.use_ratio:
-            log_bayes_factors = np.array([res.log_bayes_factor for res in self])
-            result.log_bayes_factor = logsumexp(log_bayes_factors, b=1. / len(self))
-            result.log_evidence = result.log_bayes_factor + result.log_noise_evidence
-            result_weights = np.exp(log_bayes_factors - np.max(log_bayes_factors))
-        else:
-            log_evidences = np.array([res.log_evidence for res in self])
-            result.log_evidence = logsumexp(log_evidences, b=1. / len(self))
-            result_weights = np.exp(log_evidences - np.max(log_evidences))
+
+        # Combine evidences
+        log_evidences = np.array([res.log_evidence for res in self])
+        result.log_evidence = logsumexp(log_evidences, b=1. / len(self))
+        result.log_bayes_factor = result.log_evidence - result.log_noise_evidence
+
+        # Propogate uncertainty in combined evidence
         log_errs = [res.log_evidence_err for res in self if np.isfinite(res.log_evidence_err)]
         if len(log_errs) > 0:
             result.log_evidence_err = 0.5 * logsumexp(2 * np.array(log_errs), b=1. / len(self))
         else:
             result.log_evidence_err = np.nan
+
+        # Combined posteriors with a weighting
+        result_weights = np.exp(log_evidences - np.max(log_evidences))
         posteriors = list()
         for res, frac in zip(self, result_weights):
             selected_samples = (np.random.uniform(size=len(res.posterior)) < frac)
             posteriors.append(res.posterior[selected_samples])
+
         # remove original nested_samples
         result.nested_samples = None
         result.sampler_kwargs = None

bilby/core/sampler/__init__.py

@@ -4,7 +4,7 @@ import datetime
 from collections import OrderedDict
 
 import bilby
-from ..utils import command_line_args, logger
+from ..utils import command_line_args, logger, loaded_modules_dict
 from ..prior import PriorDict, DeltaFunction
 
 from .base_sampler import Sampler, SamplingMarginalisedParameterError
@@ -21,13 +21,14 @@ from .pymc3 import Pymc3
 from .pymultinest import Pymultinest
 from .ultranest import Ultranest
 from .fake_sampler import FakeSampler
+from .dnest4 import DNest4
 from . import proposal
 
 IMPLEMENTED_SAMPLERS = {
-    'cpnest': Cpnest, 'dynamic_dynesty': DynamicDynesty, 'dynesty': Dynesty,
-    'emcee': Emcee, 'kombine': Kombine, 'nestle': Nestle, 'ptemcee': Ptemcee,
-    'ptmcmcsampler': PTMCMCSampler, 'pymc3': Pymc3, 'pymultinest': Pymultinest,
-    'pypolychord': PyPolyChord, 'ultranest': Ultranest,
+    'cpnest': Cpnest, 'dnest4': DNest4, 'dynamic_dynesty': DynamicDynesty,
+    'dynesty': Dynesty, 'emcee': Emcee, 'kombine': Kombine, 'nestle': Nestle,
+    'ptemcee': Ptemcee, 'ptmcmcsampler': PTMCMCSampler, 'pymc3': Pymc3,
+    'pymultinest': Pymultinest, 'pypolychord': PyPolyChord, 'ultranest': Ultranest,
     'fake_sampler': FakeSampler}
 
 if command_line_args.sampler_help:
@@ -141,6 +142,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
     if meta_data is None:
         meta_data = dict()
     meta_data['likelihood'] = likelihood.meta_data
+    meta_data["loaded_modules"] = loaded_modules_dict()
 
     if command_line_args.bilby_zero_likelihood_mode:
         from bilby.core.likelihood import ZeroLikelihood

bilby/core/sampler/base_sampler.py

 from __future__ import absolute_import
 import datetime
+import distutils.dir_util
 import numpy as np
+import os
+import tempfile
 
 from pandas import DataFrame
 
-from ..utils import logger, command_line_args, Counter
+from ..utils import logger, check_directory_exists_and_if_not_mkdir, command_line_args, Counter
 from ..prior import Prior, PriorDict, DeltaFunction, Constraint
 from ..result import Result, read_in_result
 
@@ -541,7 +544,8 @@ class Sampler(object):
 
 
 class NestedSampler(Sampler):
-    npoints_equiv_kwargs = ['nlive', 'nlives', 'n_live_points', 'npoints', 'npoint', 'Nlive', 'num_live_points']
+    npoints_equiv_kwargs = ['nlive', 'nlives', 'n_live_points', 'npoints',
+                            'npoint', 'Nlive', 'num_live_points', 'num_particles']
     walks_equiv_kwargs = ['walks', 'steps', 'nmcmc']
 
     def reorder_loglikelihoods(self, unsorted_loglikelihoods, unsorted_samples,
@@ -601,6 +605,27 @@ class NestedSampler(Sampler):
         else:
             return np.nan_to_num(-np.inf)
 
+    def _setup_run_directory(self):
+        """
+        If using a temporary directory, the output directory is moved to the
+        temporary directory.
+        Used for Dnest4, Pymultinest, and Ultranest.
+        """
+        if self.use_temporary_directory:
+            temporary_outputfiles_basename = tempfile.TemporaryDirectory().name
+            self.temporary_outputfiles_basename = temporary_outputfiles_basename
+
+            if os.path.exists(self.outputfiles_basename):
+                distutils.dir_util.copy_tree(self.outputfiles_basename, self.temporary_outputfiles_basename)
+            check_directory_exists_and_if_not_mkdir(temporary_outputfiles_basename)
+
+            self.kwargs["outputfiles_basename"] = self.temporary_outputfiles_basename
+            logger.info("Using temporary file {}".format(temporary_outputfiles_basename))
+        else:
+            check_directory_exists_and_if_not_mkdir(self.outputfiles_basename)
+            self.kwargs["outputfiles_basename"] = self.outputfiles_basename
+            logger.info("Using output file {}".format(self.outputfiles_basename))
+
 
 class MCMCSampler(Sampler):
     nwalkers_equiv_kwargs = ['nwalker', 'nwalkers', 'draws', 'Niter']

bilby/core/sampler/cpnest.py

 from __future__ import absolute_import
 
+import array
 import copy
 
 import numpy as np
@@ -88,8 +89,8 @@ class Cpnest(NestedSampler):
                 prior_samples = self.priors.sample()
                 self._update_bounds()
                 point = LivePoint(
-                    self.names, [prior_samples[name]
-                                 for name in self.names])
+                    self.names, array.array(
+                        'f', [prior_samples[name] for name in self.names]))
                 return point
 
         self._resolve_proposal_functions()

bilby/core/sampler/dnest4.py

+import os
+import shutil
+import distutils.dir_util
+import signal
+import time
+import datetime
+import sys
+
+import numpy as np
+import pandas as pd
+
+from ..utils import check_directory_exists_and_if_not_mkdir, logger
+from .base_sampler import NestedSampler
+
+
+class _DNest4Model(object):
+
+    def __init__(self, log_likelihood_func, from_prior_func, widths, centers, highs, lows):
+        """Initialize the DNest4 model.
+        Args:
+            log_likelihood_func: function
+                The loglikelihood function to use during the Nested Sampling run.
+            from_prior_func: function
+                The function to use when randomly selecting parameter vectors from the prior space.
+            widths: array_like
+                The approximate widths of the prior distrbutions.
+            centers: array_like
+                The approximate center points of the prior distributions.
+        """
+        self._log_likelihood = log_likelihood_func
+        self._from_prior = from_prior_func
+        self._widths = widths
+        self._centers = centers
+        self._highs = highs
+        self._lows = lows
+        self._n_dim = len(widths)
+        return
+
+    def log_likelihood(self, coords):
+        """The model's log_likelihood function"""
+        return self._log_likelihood(coords)
+
+    def from_prior(self):
+        """The model's function to select random points from the prior space."""
+        return self._from_prior()
+
+    def perturb(self, coords):
+        """The perturb function to perform Monte Carlo trial moves."""
+        idx = np.random.randint(self._n_dim)
+
+        coords[idx] += (self._widths[idx] * (np.random.uniform(size=1) - 0.5))
+        cw = self._widths[idx]
+        cc = self._centers[idx]
+
+        coords[idx] = self.wrap(coords[idx], (cc - 0.5 * cw), cc + 0.5 * cw)
+
+        return 0.0
+
+    @staticmethod
+    def wrap(x, minimum, maximum):
+        if maximum <= minimum:
+            raise ValueError("maximum {} <= minimum {}, when trying to wrap coordinates".format(maximum, minimum))
+        return (x - minimum) % (maximum - minimum) + minimum
+
+
+class DNest4(NestedSampler):
+
+    """
+    Bilby wrapper of DNest4
+
+    Parameters
+    ----------
+    TBD
+
+    Other Parameters
+    ----------------
+    num_particles: int
+        The number of points to use in the Nested Sampling active population.
+    max_num_levels: int
+        The max number of diffusive likelihood levels that DNest4 should initialize
+        during the Diffusive Nested Sampling run.
+    backend: str
+        The python DNest4 backend for storing the output.
+        Options are: 'memory' and 'csv'. If 'memory' the
+        DNest4 outputs are stored in memory during the run. If 'csv' the
+        DNest4 outputs are written out to files with a CSV format during
+        the run.
+        CSV backend may not be functional right now (October 2020)
+    num_steps: int
+        The number of MCMC iterations to run
+    new_level_interval: int
+        The number of moves to run before creating a new diffusive likelihood level
+    lam: float
+        Set the backtracking scale length
+    beta: float
+        Set the strength of effect to force the histogram to equal bin counts
+    seed: int
+        Set the seed for the C++ random number generator
+    verbose: Bool
+        If True, prints information during run
+    """
+
+    default_kwargs = dict(max_num_levels=20, num_steps=500,
+                          new_level_interval=10000, num_per_step=10000,
+                          thread_steps=1, num_particles=1000, lam=10.0,
+                          beta=100, seed=None, verbose=True, outputfiles_basename=None,
+                          backend='memory')
+
+    def __init__(self, likelihood, priors, outdir="outdir", label="label", use_ratio=False, plot=False,
+                 exit_code=77, skip_import_verification=False, temporary_directory=True, **kwargs):
+        super(DNest4, self).__init__(
+            likelihood=likelihood, priors=priors, outdir=outdir, label=label,
+            use_ratio=use_ratio, plot=plot, skip_import_verification=skip_import_verification,
+            exit_code=exit_code, **kwargs)
+
+        self.num_particles = self.kwargs["num_particles"]
+        self.max_num_levels = self.kwargs["max_num_levels"]
+        self._verbose = self.kwargs["verbose"]
+        self._backend = self.kwargs["backend"]
+        self.use_temporary_directory = temporary_directory
+
+        self.start_time = np.nan
+        self.sampler = None
+        self._information = np.nan
+        self._last_live_sample_info = np.nan
+        self._outputfiles_basename = None
+        self._temporary_outputfiles_basename = None
+
+        signal.signal(signal.SIGTERM, self.write_current_state_and_exit)
+        signal.signal(signal.SIGINT, self.write_current_state_and_exit)
+        signal.signal(signal.SIGALRM, self.write_current_state_and_exit)
+
+        # Get the estimates of the prior distributions' widths and centers.
+        widths = []
+        centers = []
+        highs = []
+        lows = []
+
+        samples = self.priors.sample(size=10000)
+
+        for key in self.search_parameter_keys:
+            pts = samples[key]
+            low = pts.min()
+            high = pts.max()
+            width = high - low
+            center = (high + low) / 2.0
+            widths.append(width)
+            centers.append(center)
+
+            highs.append(high)
+            lows.append(low)
+
+        self._widths = np.array(widths)
+        self._centers = np.array(centers)
+        self._highs = np.array(highs)
+        self._lows = np.array(lows)
+
+        self._dnest4_model = _DNest4Model(self.log_likelihood, self.get_random_draw_from_prior, self._widths,
+                                          self._centers, self._highs, self._lows)
+
+    def _set_backend(self):
+        import dnest4
+        if self._backend == 'csv':
+            return dnest4.backends.CSVBackend("{}/dnest4{}/".format(self.outdir, self.label), sep=" ")
+        else:
+            return dnest4.backends.MemoryBackend()
+
+    def _set_dnest4_kwargs(self):
+        dnest4_keys = ["num_steps", "new_level_interval", "lam", "beta", "seed"]
+        self.dnest4_kwargs = {key: self.kwargs[key] for key in dnest4_keys}
+
+    def run_sampler(self):
+        import dnest4
+
+        self._set_dnest4_kwargs()
+        backend = self._set_backend()
+
+        self._verify_kwargs_against_default_kwargs()
+        self._setup_run_directory()
+        self._check_and_load_sampling_time_file()
+        self.start_time = time.time()
+
+        self.sampler = dnest4.DNest4Sampler(self._dnest4_model, backend=backend)
+        out = self.sampler.sample(self.max_num_levels,
+                                  num_particles=self.num_particles,
+                                  **self.dnest4_kwargs)
+
+        for i, sample in enumerate(out):
+            if self._verbose and ((i + 1) % 100 == 0):
+                stats = self.sampler.postprocess()
+                logger.info("Iteration: {0} log(Z): {1}".format(i + 1, stats['log_Z']))
+
+        self._calculate_and_save_sampling_time()
+        self._clean_up_run_directory()
+
+        stats = self.sampler.postprocess(resample=1)
+        self.result.log_evidence = stats['log_Z']
+        self._information = stats['H']
+        self.result.log_evidence_err = np.sqrt(self._information / self.num_particles)
+
+        if self._backend == 'memory':
+            self._last_live_sample_info = pd.DataFrame(self.sampler.backend.sample_info[-1])
+            self.result.log_likelihood_evaluations = self._last_live_sample_info['log_likelihood']
+            self.result.samples = np.array(self.sampler.backend.posterior_samples)
+        else:
+            sample_info_path = './' + self.kwargs["outputfiles_basename"] + '/sample_info.txt'
+            sample_info = np.genfromtxt(sample_info_path, comments='#', names=True)
+            self.result.log_likelihood_evaluations = sample_info['log_likelihood']
+            self.result.samples = np.array(self.sampler.backend.posterior_samples)
+
+        self.result.sampler_output = out
+        self.result.outputfiles_basename = self.outputfiles_basename
+        self.result.sampling_time = datetime.timedelta(seconds=self.total_sampling_time)
+
+        self.calc_likelihood_count()
+
+        return self.result
+
+    def _translate_kwargs(self, kwargs):
+        if 'num_steps' not in kwargs:
+            for equiv in self.walks_equiv_kwargs:
+                if equiv in kwargs:
+                    kwargs['num_steps'] = kwargs.pop(equiv)
+
+    def _verify_kwargs_against_default_kwargs(self):
+        self.outputfiles_basename = self.kwargs.pop("outputfiles_basename", None)
+        super(DNest4, self)._verify_kwargs_against_default_kwargs()
+
+    def _check_and_load_sampling_time_file(self):
+        self.time_file_path = self.kwargs["outputfiles_basename"] + '/sampling_time.dat'
+        if os.path.exists(self.time_file_path):
+            with open(self.time_file_path, 'r') as time_file:
+                self.total_sampling_time = float(time_file.readline())
+        else:
+            self.total_sampling_time = 0
+
+    def _calculate_and_save_sampling_time(self):
+        current_time = time.time()
+        new_sampling_time = current_time - self.start_time
+        self.total_sampling_time += new_sampling_time
+
+        with open(self.time_file_path, 'w') as time_file:
+            time_file.write(str(self.total_sampling_time))
+
+        self.start_time = current_time
+
+    def _clean_up_run_directory(self):
+        if self.use_temporary_directory:
+            self._move_temporary_directory_to_proper_path()
+            self.kwargs["outputfiles_basename"] = self.outputfiles_basename
+
+    @property
+    def outputfiles_basename(self):
+        return self._outputfiles_basename
+
+    @outputfiles_basename.setter
+    def outputfiles_basename(self, outputfiles_basename):
+        if outputfiles_basename is None:
+            outputfiles_basename = "{}/dnest4{}/".format(self.outdir, self.label)
+        if not outputfiles_basename.endswith("/"):
+            outputfiles_basename += "/"
+        check_directory_exists_and_if_not_mkdir(self.outdir)
+        self._outputfiles_basename = outputfiles_basename
+
+    @property
+    def temporary_outputfiles_basename(self):
+        return self._temporary_outputfiles_basename
+
+    @temporary_outputfiles_basename.setter
+    def temporary_outputfiles_basename(self, temporary_outputfiles_basename):
+        if not temporary_outputfiles_basename.endswith("/"):
+            temporary_outputfiles_basename = "{}/".format(
+                temporary_outputfiles_basename
+            )
+        self._temporary_outputfiles_basename = temporary_outputfiles_basename
+        if os.path.exists(self.outputfiles_basename):
+            shutil.copytree(
+                self.outputfiles_basename, self.temporary_outputfiles_basename
+            )
+
+    def write_current_state_and_exit(self, signum=None, frame=None):
+        """ Write current state and exit on exit_code """
+        logger.info(
+            "Run interrupted by signal {}: checkpoint and exit on {}".format(
+                signum, self.exit_code
+            )
+        )
+        self._calculate_and_save_sampling_time()
+        if self.use_temporary_directory:
+            self._move_temporary_directory_to_proper_path()
+        sys.exit(self.exit_code)
+
+    def _move_temporary_directory_to_proper_path(self):
+        """
+        Move the temporary back to the proper path
+
+        Anything in the proper path at this point is removed including links
+        """
+        self._copy_temporary_directory_contents_to_proper_path()
+        shutil.rmtree(self.temporary_outputfiles_basename)
+
+    def _copy_temporary_directory_contents_to_proper_path(self):
+        """
+        Copy the temporary back to the proper path.
+        Do not delete the temporary directory.
+        """
+        logger.info(
+            "Overwriting {} with {}".format(
+                self.outputfiles_basename, self.temporary_outputfiles_basename
+            )
+        )
+        if self.outputfiles_basename.endswith('/'):
+            outputfiles_basename_stripped = self.outputfiles_basename[:-1]
+        else:
+            outputfiles_basename_stripped = self.outputfiles_basename
+        distutils.dir_util.copy_tree(self.temporary_outputfiles_basename, outputfiles_basename_stripped)

bilby/core/sampler/pymultinest.py

 import importlib
 import os
-import tempfile
 import shutil
 import distutils.dir_util
 import signal
 import time
 import datetime
+import sys
 
 import numpy as np
 
@@ -175,7 +175,7 @@ class Pymultinest(NestedSampler):
         self._calculate_and_save_sampling_time()
         if self.use_temporary_directory:
             self._move_temporary_directory_to_proper_path()
-        os._exit(self.exit_code)
+        sys.exit(self.exit_code)
 
     def _copy_temporary_directory_contents_to_proper_path(self):
         """
@@ -239,26 +239,6 @@ class Pymultinest(NestedSampler):
         self.result.sampling_time = datetime.timedelta(seconds=self.total_sampling_time)
         return self.result
 
-    def _setup_run_directory(self):
-        """
-        If using a temporary directory, the output directory is moved to the
-        temporary directory.
-        """
-        if self.use_temporary_directory:
-            temporary_outputfiles_basename = tempfile.TemporaryDirectory().name
-            self.temporary_outputfiles_basename = temporary_outputfiles_basename
-
-            if os.path.exists(self.outputfiles_basename):
-                distutils.dir_util.copy_tree(self.outputfiles_basename, self.temporary_outputfiles_basename)
-            check_directory_exists_and_if_not_mkdir(temporary_outputfiles_basename)
-
-            self.kwargs["outputfiles_basename"] = self.temporary_outputfiles_basename
-            logger.info("Using temporary file {}".format(temporary_outputfiles_basename))
-        else:
-            check_directory_exists_and_if_not_mkdir(self.outputfiles_basename)
-            self.kwargs["outputfiles_basename"] = self.outputfiles_basename
-            logger.info("Using output file {}".format(self.outputfiles_basename))
-
     def _check_and_load_sampling_time_file(self):
         self.time_file_path = self.kwargs["outputfiles_basename"] + '/sampling_time.dat'
         if os.path.exists(self.time_file_path):

bilby/core/sampler/ultranest.py

@@ -6,7 +6,6 @@ import inspect
 import os
 import shutil
 import signal
-import tempfile
 import time
 
 import numpy as np
@@ -287,6 +286,7 @@ class Ultranest(NestedSampler):
         stepsampler = self.kwargs.pop("step_sampler", None)
 
         self._setup_run_directory()
+        self.kwargs["log_dir"] = self.kwargs.pop("outputfiles_basename")
         self._check_and_load_sampling_time_file()
 
         # use reactive nested sampler when no live points are given
@@ -326,30 +326,6 @@ class Ultranest(NestedSampler):
 
         return self.result
 
-    def _setup_run_directory(self):
-        """
-        If using a temporary directory, the output directory is moved to the
-        temporary directory and symlinked back.
-        """
-        if self.use_temporary_directory:
-            temporary_outputfiles_basename = tempfile.TemporaryDirectory().name
-            self.temporary_outputfiles_basename = temporary_outputfiles_basename
-
-            if os.path.exists(self.outputfiles_basename):
-                distutils.dir_util.copy_tree(
-                    self.outputfiles_basename, self.temporary_outputfiles_basename
-                )
-            check_directory_exists_and_if_not_mkdir(temporary_outputfiles_basename)
-
-            self.kwargs["log_dir"] = self.temporary_outputfiles_basename
-            logger.info(
-                "Using temporary file {}".format(temporary_outputfiles_basename)
-            )
-        else:
-            check_directory_exists_and_if_not_mkdir(self.outputfiles_basename)
-            self.kwargs["log_dir"] = self.outputfiles_basename
-            logger.info("Using output file {}".format(self.outputfiles_basename))
-
     def _clean_up_run_directory(self):
         if self.use_temporary_directory:
             self._move_temporary_directory_to_proper_path()

bilby/core/utils.py

@@ -4,6 +4,7 @@ from distutils.spawn import find_executable
 import logging
 import os
 import shutil
+import sys
 from math import fmod
 import argparse
 import inspect
@@ -1269,6 +1270,15 @@ def get_function_path(func):
         return func
 
 
+def loaded_modules_dict():
+    module_names = sys.modules.keys()
+    vdict = {}
+    for key in module_names:
+        if "." not in key:
+            vdict[key] = str(getattr(sys.modules[key], "__version__", "N/A"))
+    return vdict
+
+
 class IllegalDurationAndSamplingFrequencyException(Exception):
     pass
 

bilby/gw/likelihood.py

@@ -1371,7 +1371,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
 
 def get_binary_black_hole_likelihood(interferometers):
-    """ A rapper to quickly set up a likelihood for BBH parameter estimation
+    """ A wrapper to quickly set up a likelihood for BBH parameter estimation
 
     Parameters
     ----------

docs/basics-of-parameter-estimation.txt

@@ -6,7 +6,7 @@ In this example, we'll go into some of the basics of parameter estimation and
 how they are implemented in :code:`bilby`.
 
 Firstly, consider a situation where you have discrete data :math:`\{y_0,
-y_1,\ldots, y_n\}` taken at a set of times :math:`\{t_0, t_1, \ldots, y_n\}`.
+y_1,\ldots, y_n\}` taken at a set of times :math:`\{t_0, t_1, \ldots, t_n\}`.
 Further, we know that this data is generated by a process which can be modelled
 by a linear function of the form :math:`y(t) = m t + c`. We will refer to
 this model as :math:`H`. Given a set of data,

docs/prior.txt

@@ -45,7 +45,7 @@ which provides extra functionality. For example, to sample from the prior:
 
 .. code:: python
 
-   >>> priors = bilby.core.priors.PriorDict()
+   >>> priors = bilby.core.prior.PriorDict()
    >>> priors['a'] = bilby.prior.Uniform(minimum=0, maximum=10, name='a')
    >>> priors['b'] = bilby.prior.Uniform(minimum=0, maximum=10, name='b')
    >>> priors.sample()
@@ -89,7 +89,7 @@ matrix and standard deviations, e.g.:
    >>> names = ['a', 'b']  # set the parameter names
    >>> mu = [0., 5.]  # the means of the parameters
    >>> cov = [[2., 0.7], [0.7, 3.]]  # the covariance matrix
-   >>> mvg = bilby.core.priors.MultivariateGaussianDist(names, mus=mu, covs=cov)
+   >>> mvg = bilby.core.prior.MultivariateGaussianDist(names, mus=mu, covs=cov)
 
 It is also possible to define a mixture model of multiple multivariate Gaussian modes of
 different weights if required, e.g.:
@@ -100,7 +100,7 @@ different weights if required, e.g.:
    >>> mu = [[0., 5.], [2., 7.]]  # the means of the parameters
    >>> cov = [[[2., 0.7], [0.7, 3.]], [[1., -0.9], [-0.9, 5.]]]  # the covariance matrix
    >>> weights = [0.3, 0.7]  # weights of each mode
-   >>> mvg = bilby.core.priors.MultivariateGaussianDist(names, mus=mu, covs=cov, nmodes=2, weights=weights)
+   >>> mvg = bilby.core.prior.MultivariateGaussianDist(names, mus=mu, covs=cov, nmodes=2, weights=weights)
 
 The distribution can also have hard bounds on each parameter by supplying them.
 

docs/samplers.txt

@@ -89,6 +89,11 @@ UltraNest
 
 .. autoclass:: bilby.core.sampler.ultranest.Ultranest
 
+Dnest4
+=========
+
+.. autoclass:: bilby.core.sampler.dnest4.Dnest4
+
 -------------
 MCMC samplers
 -------------

examples/core_examples/slabspike_example.py

+#!/usr/bin/env python
+"""
+An example of how to use slab-and-spike priors in bilby.
+In this example we look at a simple example with the sum
+of two Gaussian distributions, and we try to fit with
+up to three Gaussians.
+
+"""
+
+import bilby
+import numpy as np
+import matplotlib.pyplot as plt
+
+outdir = 'outdir'
+label = 'slabspike'
+bilby.utils.check_directory_exists_and_if_not_mkdir(outdir)
+
+
+# Here we define our model. We want to inject two Gaussians and recover with up to three.
+def gaussian(xs, amplitude, mu, sigma):
+    return amplitude / np.sqrt(2 * np.pi * sigma**2) * np.exp(-0.5 * (xs - mu)**2 / sigma**2)
+
+
+def triple_gaussian(xs, amplitude_0, amplitude_1, amplitude_2, mu_0, mu_1, mu_2, sigma_0, sigma_1, sigma_2, **kwargs):
+    return \
+        gaussian(xs, amplitude_0, mu_0, sigma_0) + \
+        gaussian(xs, amplitude_1, mu_1, sigma_1) + \
+        gaussian(xs, amplitude_2, mu_2, sigma_2)
+
+
+# Let's create our data set. We create 200 points on a grid.
+
+xs = np.linspace(-5, 5, 200)
+dx = xs[1] - xs[0]
+
+# Note for our injection parameters we set the amplitude of the second component to 0.
+injection_params = dict(amplitude_0=-3, mu_0=-4, sigma_0=4,
+                        amplitude_1=0, mu_1=0, sigma_1=1,
+                        amplitude_2=4, mu_2=3, sigma_2=3)
+
+# We calculate the injected curve and add some Gaussian noise on the data points
+sigma = 0.02
+p = bilby.core.prior.Gaussian(mu=0, sigma=sigma)
+ys = triple_gaussian(xs=xs, **injection_params) + p.sample(len(xs))
+
+plt.errorbar(xs, ys, yerr=sigma, fmt=".k", capsize=0, label='Injected data')
+plt.plot(xs, triple_gaussian(xs=xs, **injection_params), label='True signal')
+plt.legend()
+plt.savefig(f'{outdir}/{label}_injected_data')
+plt.clf()
+
+
+# Now we want to set up our priors.
+priors = bilby.core.prior.PriorDict()
+# For the slab-and-spike prior, we first need to define the 'slab' part, which is just a regular bilby prior.
+amplitude_slab_0 = bilby.core.prior.Uniform(minimum=-10, maximum=10, name='amplitude_0', latex_label='$A_0$')
+amplitude_slab_1 = bilby.core.prior.Uniform(minimum=-10, maximum=10, name='amplitude_1', latex_label='$A_1$')
+amplitude_slab_2 = bilby.core.prior.Uniform(minimum=-10, maximum=10, name='amplitude_2', latex_label='$A_2$')
+# We do the following to create the slab-and-spike prior. The spike height is somewhat arbitrary and can
+# be corrected in post-processing.
+priors['amplitude_0'] = bilby.core.prior.SlabSpikePrior(slab=amplitude_slab_0, spike_location=0, spike_height=0.1)
+priors['amplitude_1'] = bilby.core.prior.SlabSpikePrior(slab=amplitude_slab_1, spike_location=0, spike_height=0.1)
+priors['amplitude_2'] = bilby.core.prior.SlabSpikePrior(slab=amplitude_slab_2, spike_location=0, spike_height=0.1)
+# Our problem has a degeneracy in the ordering. In general, this problem is somewhat difficult to resolve properly.
+# See e.g. https://github.com/GregoryAshton/kookaburra/blob/master/src/priors.py#L72 for an implementation.
+# We resolve this by not letting the priors overlap in this case.
+priors['mu_0'] = bilby.core.prior.Uniform(minimum=-5, maximum=-2, name='mu_0', latex_label='$\mu_0$')
+priors['mu_1'] = bilby.core.prior.Uniform(minimum=-2, maximum=2, name='mu_1', latex_label='$\mu_1$')
+priors['mu_2'] = bilby.core.prior.Uniform(minimum=2, maximum=5, name='mu_2', latex_label='$\mu_2$')
+priors['sigma_0'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_0', latex_label='$\sigma_0$')
+priors['sigma_1'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_1', latex_label='$\sigma_1$')
+priors['sigma_2'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_2', latex_label='$\sigma_2$')
+
+# Setting up the likelihood and running the samplers works the same as elsewhere.
+likelihood = bilby.core.likelihood.GaussianLikelihood(x=xs, y=ys, func=triple_gaussian, sigma=sigma)
+result = bilby.run_sampler(likelihood=likelihood, priors=priors, outdir=outdir, label=label,
+                           sampler='dynesty', nlive=400)
+
+result.plot_corner(truths=injection_params)
+
+
+# Let's also plot the maximum likelihood fit along with the data.
+max_like_params = result.posterior.iloc[-1]
+plt.errorbar(xs, ys, yerr=sigma, fmt=".k", capsize=0, label='Injected data')
+plt.plot(xs, triple_gaussian(xs=xs, **injection_params), label='True signal')
+plt.plot(xs, triple_gaussian(xs=xs, **max_like_params), label='Max likelihood fit')
+plt.legend()
+plt.savefig(f'{outdir}/{label}_max_likelihood_recovery')
+plt.clf()
+
+# Finally, we can check what fraction of amplitude samples are exactly on the spike.
+spike_samples_0 = len(np.where(result.posterior['amplitude_0'] == 0.0)[0]) / len(result.posterior)
+spike_samples_1 = len(np.where(result.posterior['amplitude_1'] == 0.0)[0]) / len(result.posterior)
+spike_samples_2 = len(np.where(result.posterior['amplitude_2'] == 0.0)[0]) / len(result.posterior)
+print(f"{spike_samples_0 * 100:.2f}% of amplitude_0 samples are exactly 0.0")
+print(f"{spike_samples_1 * 100:.2f}% of amplitude_1 samples are exactly 0.0")
+print(f"{spike_samples_2 * 100:.2f}% of amplitude_2 samples are exactly 0.0")