bilby/core/sampler/ptmcmc.py

@@ -41,6 +41,8 @@ class PTMCMCSampler(MCMCSampler):
 
     """
 
+    sampler_name = "ptmcmcsampler"
+    abbreviation = "ptmcmc_temp"
     default_kwargs = {
         "p0": None,
         "Niter": 2 * 10**4 + 1,

bilby/core/sampler/pymc.py

@@ -52,6 +52,7 @@ class Pymc(MCMCSampler):
 
     """
 
+    sampler_name = "pymc"
     default_kwargs = dict(
         draws=500,
         step=None,

bilby/core/sampler/pymultinest.py

@@ -34,6 +34,8 @@ class Pymultinest(_TemporaryFileSamplerMixin, NestedSampler):
 
     """
 
+    sampler_name = "pymultinest"
+    abbreviation = "pm"
     default_kwargs = dict(
         importance_nested_sampling=False,
         resume=True,

bilby/core/sampler/ultranest.py

@@ -38,6 +38,8 @@ class Ultranest(_TemporaryFileSamplerMixin, NestedSampler):
         stepping behaviour is used.
     """
 
+    sampler_name = "ultranest"
+    abbreviation = "ultra"
     default_kwargs = dict(
         resume=True,
         show_status=True,

bilby/core/sampler/zeus.py

@@ -38,6 +38,7 @@ class Zeus(Emcee):
 
     """
 
+    sampler_name = "zeus"
     default_kwargs = dict(
         nwalkers=500,
         args=[],

bilby/core/utils/calculus.py

 import math
-from numbers import Number
+
 import numpy as np
-from scipy.interpolate import interp2d
+from scipy.interpolate import RectBivariateSpline
 from scipy.special import logsumexp
 
 from .log import logger
@@ -189,79 +189,34 @@ def logtrapzexp(lnf, dx):
     return C + logsumexp([logsumexp(lnfdx1), logsumexp(lnfdx2)])
 
 
-class UnsortedInterp2d(interp2d):
-    def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
-        """Modified version of the interp2d call method.
-
-        This avoids the outer product that is done when two numpy
-        arrays are passed.
-
-        Parameters
-        ==========
-        x: See superclass
-        y: See superclass
-        dx: See superclass
-        dy: See superclass
-        assume_sorted: bool, optional
-            This is just a place holder to prevent a warning.
-            Overwriting this will not do anything
+class BoundedRectBivariateSpline(RectBivariateSpline):
 
-        Returns
-        =======
-        array_like: See superclass
+    def __init__(self, x, y, z, bbox=[None] * 4, kx=3, ky=3, s=0, fill_value=None):
+        self.x_min, self.x_max, self.y_min, self.y_max = bbox
+        if self.x_min is None:
+            self.x_min = min(x)
+        if self.x_max is None:
+            self.x_max = max(x)
+        if self.y_min is None:
+            self.y_min = min(y)
+        if self.y_max is None:
+            self.y_max = max(y)
+        self.fill_value = fill_value
+        super().__init__(x=x, y=y, z=z, bbox=bbox, kx=kx, ky=ky, s=s)
 
-        """
-        from scipy.interpolate.dfitpack import bispeu
-
-        x, y = self._sanitize_inputs(x, y)
+    def __call__(self, x, y, dx=0, dy=0, grid=False):
+        result = super().__call__(x=x, y=y, dx=dx, dy=dy, grid=grid)
         out_of_bounds_x = (x < self.x_min) | (x > self.x_max)
         out_of_bounds_y = (y < self.y_min) | (y > self.y_max)
         bad = out_of_bounds_x | out_of_bounds_y
-        if isinstance(x, Number) and isinstance(y, Number):
+        result[bad] = self.fill_value
+        if result.size == 1:
             if bad:
-                output = self.fill_value
-                ier = 0
+                return self.fill_value
             else:
-                output, ier = bispeu(*self.tck, x, y)
-                output = float(output)
+                return result.item()
         else:
-            output = np.empty_like(x)
-            output[bad] = self.fill_value
-            if np.any(~bad):
-                output[~bad], ier = bispeu(*self.tck, x[~bad], y[~bad])
-            else:
-                ier = 0
-        if ier == 10:
-            raise ValueError("Invalid input data")
-        elif ier:
-            raise TypeError("An error occurred")
-        return output
-
-    @staticmethod
-    def _sanitize_inputs(x, y):
-        if isinstance(x, np.ndarray) and x.size == 1:
-            x = float(x)
-        if isinstance(y, np.ndarray) and y.size == 1:
-            y = float(y)
-        if isinstance(x, np.ndarray) and isinstance(y, np.ndarray):
-            original_shapes = (x.shape, y.shape)
-            if x.shape != y.shape:
-                while x.ndim > y.ndim:
-                    y = np.expand_dims(y, -1)
-                while y.ndim > x.ndim:
-                    x = np.expand_dims(x, -1)
-            try:
-                x = x * np.ones(y.shape)
-                y = y * np.ones(x.shape)
-            except ValueError:
-                raise ValueError(
-                    f"UnsortedInterp2d received incompatibly shaped arrays: {original_shapes}"
-                )
-        elif isinstance(x, np.ndarray) and not isinstance(y, np.ndarray):
-            y = y * np.ones_like(x)
-        elif not isinstance(x, np.ndarray) and isinstance(y, np.ndarray):
-            x = x * np.ones_like(y)
-        return x, y
+            return result
 
 
 def round_up_to_power_of_two(x):

bilby/core/utils/io.py

@@ -154,13 +154,16 @@ def decode_bilby_json(dct):
         try:
             cls = getattr(import_module(dct["__module__"]), dct["__name__"])
         except AttributeError:
-            logger.debug(
+            logger.warning(
                 "Unknown prior class for parameter {}, defaulting to base Prior object".format(
                     dct["kwargs"]["name"]
                 )
             )
             from ..prior import Prior
 
+            for key in list(dct["kwargs"].keys()):
+                if key not in ["name", "latex_label", "unit", "minimum", "maximum", "boundary"]:
+                    dct["kwargs"].pop(key)
             cls = Prior
         obj = cls(**dct["kwargs"])
         return obj
@@ -264,9 +267,9 @@ def encode_for_hdf5(key, item):
 
     if isinstance(item, np.int_):
         item = int(item)
-    elif isinstance(item, np.float_):
+    elif isinstance(item, np.float64):
         item = float(item)
-    elif isinstance(item, np.complex_):
+    elif isinstance(item, np.complex128):
         item = complex(item)
     if isinstance(item, np.ndarray):
         # Numpy's wide unicode strings are not supported by hdf5

bilby/gw/conversion.py

@@ -2261,8 +2261,8 @@ def compute_snrs(sample, likelihood, npool=1):
                     snr_update = new_sample[ii].snrs_as_sample
                     for key, val in snr_update.items():
                         snr_updates[f"{ifo.name}_{key}"].append(val)
-            for k, v in snr_updates.items():
-                sample[k] = v
+                for k, v in snr_updates.items():
+                    sample[k] = v
     else:
         logger.debug('Not computing SNRs.')
 

bilby/gw/detector/interferometer.py

@@ -616,19 +616,95 @@ class Interferometer(object):
             power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask],
             duration=self.strain_data.duration)
 
+    def whiten_frequency_series(self, frequency_series : np.array) -> np.array:
+        """Whitens a frequency series with the noise properties of the detector
+
+        .. math::
+            \\tilde{a}_w(f) = \\tilde{a}(f) \\sqrt{\\frac{4}{T S_n(f)}}
+
+        Such that
+
+        .. math::
+            Var(n) = \\frac{1}{N} \\sum_k=0^N n_W(f_k)n_W^*(f_k) = 2
+
+        Where the factor of two is due to the independent real and imaginary
+        components.
+
+        Parameters
+        ==========
+        frequency_series : np.array
+            The frequency series, whitened by the ASD
+        """
+        return frequency_series / (self.amplitude_spectral_density_array * np.sqrt(self.duration / 4))
+
+    def get_whitened_time_series_from_whitened_frequency_series(
+        self,
+        whitened_frequency_series : np.array
+    ) -> np.array:
+        """Gets the whitened time series from a whitened frequency series.
+
+        This ifft's and also applies a windowing factor,
+        since when f_min and f_max are set bilby applies a mask to the series.
+
+        Per 6.2a-b in https://arxiv.org/pdf/gr-qc/0509116 since our window
+        is just a band pass,
+        this coefficient is :math:`w/W` where
+
+        .. math::
+
+            W = \\frac{1}{N} \\sum_{k=0}^N w^2[j]
+
+        Since our window :math:`w` is simply 1 or 0, depending on the mask, we get
+
+        .. math::
+
+            W = \\frac{1}{N} \\sum_{k=0}^N \\Theta(f_{max} - f_k)\\Theta(f_k - f_{min})}
+
+        and accordingly the termwise window factor is
+
+        .. math::
+            w = \\sqrt{N W} = \\sqrt{\\sum_{k=0}^N \\Theta(f_{max} - f_k)\\Theta(f_k - f_{min})}}
+
+        """
+        frequency_window_factor = (
+            np.sum(self.frequency_mask)
+            / len(self.frequency_mask)
+        )
+
+        whitened_time_series = (
+            np.fft.irfft(whitened_frequency_series)
+            * np.sqrt(np.sum(self.frequency_mask)) / frequency_window_factor
+        )
+
+        return whitened_time_series
+
     @property
     def whitened_frequency_domain_strain(self):
-        """ Calculates the whitened data by dividing the frequency domain data by
-        ((amplitude spectral density) * (duration / 4) ** 0.5). The resulting
-        data will have unit variance.
+        r"""Whitens the frequency domain data by dividing through by ASD,
+        with appropriate normalization.
+
+        See `whiten_frequency_series()` for details.
 
         Returns
         =======
         array_like: The whitened data
         """
-        return self.strain_data.frequency_domain_strain / (
-            self.amplitude_spectral_density_array * np.sqrt(self.duration / 4)
-        )
+        return self.whiten_frequency_series(self.strain_data.frequency_domain_strain)
+
+    @property
+    def whitened_time_domain_strain(self) -> np.array:
+        """Calculates the whitened time domain strain
+        by iffting the whitened frequency domain strain,
+        with the appropriate normalization.
+
+        See `get_whitened_time_series_from_whitened_frequency_series()` for details
+
+        Returns
+        =======
+        array_like
+            The whitened data in the time domain
+        """
+        return self.get_whitened_time_series_from_whitened_frequency_series(self.whitened_frequency_domain_strain)
 
     def save_data(self, outdir, label=None):
         """ Creates save files for interferometer data in plain text format.

bilby/gw/likelihood/base.py

@@ -7,7 +7,7 @@ import numpy as np
 from scipy.special import logsumexp
 
 from ...core.likelihood import Likelihood
-from ...core.utils import logger, UnsortedInterp2d, create_time_series
+from ...core.utils import logger, BoundedRectBivariateSpline, create_time_series
 from ...core.prior import Interped, Prior, Uniform, DeltaFunction
 from ..detector import InterferometerList, get_empty_interferometer, calibration
 from ..prior import BBHPriorDict, Cosmological
@@ -752,7 +752,7 @@ class GravitationalWaveTransient(Likelihood):
             d_inner_h_ref = np.real(d_inner_h_ref)
 
         return self._interp_dist_margd_loglikelihood(
-            d_inner_h_ref, h_inner_h_ref)
+            d_inner_h_ref, h_inner_h_ref, grid=False)
 
     def phase_marginalized_likelihood(self, d_inner_h, h_inner_h):
         d_inner_h = ln_i0(abs(d_inner_h))
@@ -891,9 +891,9 @@ class GravitationalWaveTransient(Likelihood):
                 self._create_lookup_table()
         else:
             self._create_lookup_table()
-        self._interp_dist_margd_loglikelihood = UnsortedInterp2d(
+        self._interp_dist_margd_loglikelihood = BoundedRectBivariateSpline(
             self._d_inner_h_ref_array, self._optimal_snr_squared_ref_array,
-            self._dist_margd_loglikelihood_array, kind='cubic', fill_value=-np.inf)
+            self._dist_margd_loglikelihood_array.T, fill_value=-np.inf)
 
     @property
     def cached_lookup_table_filename(self):

bilby/gw/likelihood/roq.py

@@ -136,13 +136,28 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
                           linear_matrix['duration_s'][()])],
                         dtype=[('flow', float), ('fhigh', float), ('seglen', float)]
                     )
-                elif is_hdf5_quadratic:
-                    self.roq_params = np.array(
-                        [(quadratic_matrix['minimum_frequency_hz'][()],
-                          quadratic_matrix['maximum_frequency_hz'][()],
-                          quadratic_matrix['duration_s'][()])],
-                        dtype=[('flow', float), ('fhigh', float), ('seglen', float)]
-                    )
+                if is_hdf5_quadratic:
+                    if self.roq_params is None:
+                        self.roq_params = np.array(
+                            [(quadratic_matrix['minimum_frequency_hz'][()],
+                              quadratic_matrix['maximum_frequency_hz'][()],
+                              quadratic_matrix['duration_s'][()])],
+                            dtype=[('flow', float), ('fhigh', float), ('seglen', float)]
+                        )
+                    else:
+                        self.roq_params['flow'] = max(
+                            self.roq_params['flow'], quadratic_matrix['minimum_frequency_hz'][()]
+                        )
+                        self.roq_params['fhigh'] = min(
+                            self.roq_params['fhigh'], quadratic_matrix['maximum_frequency_hz'][()]
+                        )
+                        self.roq_params['seglen'] = min(
+                            self.roq_params['seglen'], quadratic_matrix['duration_s'][()]
+                        )
+            if self.roq_params is not None:
+                for ifo in self.interferometers:
+                    self.perform_roq_params_check(ifo)
+
             self.weights = dict()
             self._set_weights(linear_matrix=linear_matrix, quadratic_matrix=quadratic_matrix)
             if is_hdf5_linear:
@@ -158,9 +173,10 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
         for basis_type in ['linear', 'quadratic']:
             number_of_bases = getattr(self, f'number_of_bases_{basis_type}')
             if number_of_bases > 1:
-                self._verify_prior_ranges_and_frequency_nodes(basis_type)
+                self._verify_numbers_of_prior_ranges_and_frequency_nodes(basis_type)
             else:
                 self._check_frequency_nodes_exist_for_single_basis(basis_type)
+            self._verify_prior_ranges(basis_type)
 
         self._set_unique_frequency_nodes_and_inverse()
         # need to fill waveform_arguments here if single basis is used, as they will never be updated.
@@ -171,7 +187,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             self._waveform_generator.waveform_arguments['linear_indices'] = linear_indices
             self._waveform_generator.waveform_arguments['quadratic_indices'] = quadratic_indices
 
-    def _verify_prior_ranges_and_frequency_nodes(self, basis_type):
+    def _verify_numbers_of_prior_ranges_and_frequency_nodes(self, basis_type):
         """
         Check if self.weights contains lists of prior ranges and frequency nodes, and their sizes are equal to the
         number of bases.
@@ -205,6 +221,35 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
                 raise ValueError(
                     f'The number of arrays of frequency nodes does not match the number of {basis_type} bases')
 
+    def _verify_prior_ranges(self, basis_type):
+        """Check if the union of prior ranges is within the ROQ basis bounds.
+
+        Parameters
+        ==========
+        basis_type: str
+
+        """
+        key = f'prior_range_{basis_type}'
+        if key not in self.weights:
+            return
+        prior_ranges = self.weights[key]
+        for param_name, prior_ranges_of_this_param in prior_ranges.items():
+            prior_minimum = self.priors[param_name].minimum
+            basis_minimum = np.min(prior_ranges_of_this_param[:, 0])
+            if prior_minimum < basis_minimum:
+                raise BilbyROQParamsRangeError(
+                    f"Prior minimum of {param_name} {prior_minimum} less "
+                    f"than ROQ basis bound {basis_minimum}"
+                )
+
+            prior_maximum = self.priors[param_name].maximum
+            basis_maximum = np.max(prior_ranges_of_this_param[:, 1])
+            if prior_maximum > basis_maximum:
+                raise BilbyROQParamsRangeError(
+                    f"Prior maximum of {param_name} {prior_maximum} greater "
+                    f"than ROQ basis bound {basis_maximum}"
+                )
+
     def _check_frequency_nodes_exist_for_single_basis(self, basis_type):
         """
         For a single-basis case, frequency nodes should be contained in self._waveform_generator.waveform_arguments or
@@ -701,6 +746,8 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
                         roq_scale_factor = 1.
                     prior_ranges[param_name] = matrix[key][param_name][()] * roq_scale_factor
                 selected_idxs, selected_prior_ranges = self._select_prior_ranges(prior_ranges)
+                if len(selected_idxs) == 0:
+                    raise BilbyROQParamsRangeError(f"There are no {basis_type} ROQ bases within the prior range.")
                 self.weights[key] = selected_prior_ranges
                 idxs_in_prior_range[basis_type] = selected_idxs
             else:
@@ -725,7 +772,6 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             ifo_idxs = {}
             for ifo in self.interferometers:
                 if self.roq_params is not None:
-                    self.perform_roq_params_check(ifo)
                     # Get scaled ROQ quantities
                     roq_scaled_minimum_frequency = self.roq_params['flow'] * self.roq_scale_factor
                     roq_scaled_maximum_frequency = self.roq_params['fhigh'] * self.roq_scale_factor
@@ -1110,11 +1156,11 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             f_high: float
                 The maximum frequency which must be considered
             """
-            from scipy.integrate import simps
+            from scipy.integrate import simpson
             integrand1 = np.power(freq, -7. / 3) / psd
-            integral1 = simps(integrand1, freq)
+            integral1 = simpson(y=integrand1, x=freq)
             integrand3 = np.power(freq, 2. / 3.) / (psd * integral1)
-            f_3_bar = simps(integrand3, freq)
+            f_3_bar = simpson(y=integrand3, x=freq)
 
             f_high = scaling * f_3_bar**(1 / 3)
 

bilby/gw/prior.py

@@ -1362,6 +1362,7 @@ class HealPixMapPriorDist(BaseJointPriorDist):
         else:
             self.distance = False
             self.prob = self.hp.read_map(hp_file)
+        self.prob = self._check_norm(self.prob)
 
         super(HealPixMapPriorDist, self).__init__(names=names, bounds=bounds)
         self.distname = "hpmap"
@@ -1436,7 +1437,7 @@ class HealPixMapPriorDist(BaseJointPriorDist):
                 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])
+            sample = np.vstack([sample[:, 0], sample[:, 1], dist_samples])
         return sample.reshape((-1, self.num_vars))
 
     def update_distance(self, pix_idx):
@@ -1458,7 +1459,7 @@ class HealPixMapPriorDist(BaseJointPriorDist):
         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)
+        pdfs = self.rs ** 2 * self.distance_pdf(self.rs)
         cdfs = np.cumsum(pdfs) / np.sum(pdfs)
         self.distance_icdf = interp1d(cdfs, self.rs)
 
@@ -1501,9 +1502,7 @@ class HealPixMapPriorDist(BaseJointPriorDist):
         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 = random.rng.choice(pixel_choices, size=size, p=pixel_probs, replace=True)
+        sample_pix = random.rng.choice(self.npix, size=size, p=self.prob, replace=True)
         sample = np.empty((size, self.num_vars))
         for samp in range(size):
             theta, ra = self.hp.pix2ang(self.nside, sample_pix[samp])

bilby/gw/result.py

@@ -377,10 +377,6 @@ class CompactBinaryCoalescenceResult(CoreResult):
         logger.debug("Downsampling frequency mask to {} values".format(
             len(frequency_idxs))
         )
-        frequency_window_factor = (
-            np.sum(interferometer.frequency_mask)
-            / len(interferometer.frequency_mask)
-        )
         plot_times = interferometer.time_array[time_idxs]
         plot_times -= interferometer.strain_data.start_time
         start_time -= interferometer.strain_data.start_time
@@ -451,11 +447,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
                 fig.add_trace(
                     go.Scatter(
                         x=plot_times,
-                        y=np.fft.irfft(
-                            interferometer.whitened_frequency_domain_strain
-                            * np.sqrt(np.sum(interferometer.frequency_mask))
-                            / frequency_window_factor
-                        )[time_idxs],
+                        y=interferometer.whitened_time_domain_strain[time_idxs],
                         fill=None,
                         mode='lines', line_color=DATA_COLOR,
                         opacity=0.5,
@@ -478,11 +470,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
                     interferometer.amplitude_spectral_density_array[frequency_idxs],
                     color=DATA_COLOR, label='ASD')
                 axs[1].plot(
-                    plot_times, np.fft.irfft(
-                        interferometer.whitened_frequency_domain_strain
-                        * np.sqrt(np.sum(interferometer.frequency_mask))
-                        / frequency_window_factor
-                    )[time_idxs],
+                    plot_times, interferometer.whitened_time_domain_strain[time_idxs],
                     color=DATA_COLOR, alpha=0.3)
             logger.debug('Plotted interferometer data.')
 
@@ -493,10 +481,10 @@ class CompactBinaryCoalescenceResult(CoreResult):
             wf_pols = waveform_generator.frequency_domain_strain(params)
             fd_waveform = interferometer.get_detector_response(wf_pols, params)
             fd_waveforms.append(fd_waveform[frequency_idxs])
-            td_waveform = infft(
-                fd_waveform * np.sqrt(2. / interferometer.sampling_frequency) /
-                interferometer.amplitude_spectral_density_array,
-                self.sampling_frequency)[time_idxs]
+            whitened_fd_waveform = interferometer.whiten_frequency_series(fd_waveform)
+            td_waveform = interferometer.get_whitened_time_series_from_whitened_frequency_series(
+                whitened_fd_waveform
+            )[time_idxs]
             td_waveforms.append(td_waveform)
         fd_waveforms = asd_from_freq_series(
             fd_waveforms,

containers/env-template.yml

@@ -3,6 +3,7 @@
 channels:
   - conda-forge
   - defaults
+  - pytorch
 dependencies:
   - pip
   - setuptools
@@ -21,7 +22,7 @@ dependencies:
   - black
   - pytest-cov
   - pytest-requires
-  - arviz
+  - conda-forge::arviz
   - parameterized
   - scikit-image
   - celerite
@@ -31,6 +32,7 @@ dependencies:
   - cython
   - pytables
   - pytorch
+  - cpuonly
   - python-lalsimulation
   - bilby.cython
   - pyseobnr
@@ -67,6 +69,7 @@ dependencies:
   - nbconvert
   - twine
   - glasflow
+  - myst-parser
   - pip:
     - autodoc
     - ipykernel

docs/citing-bilby.txt

@@ -58,7 +58,7 @@ If you use any of the accelerated likelihoods like `ROQGravitationalWaveTransien
 - `ROQGravitationalWaveTransient`
  .. code:: bibtex
 
-     @article{roq_paper,
+     @article{roq_paper_1,
         author = {Smith, Rory and Field, Scott E. and Blackburn, Kent and Haster, Carl-Johan and P\"urrer, Michael and Raymond, Vivien and Schmidt, Patricia},
         title = "{Fast and accurate inference on gravitational waves from precessing compact binaries}",
         eprint = "1604.08253",
@@ -73,6 +73,21 @@ If you use any of the accelerated likelihoods like `ROQGravitationalWaveTransien
         year = "2016"
     }
 
+    @article{roq_paper_2,
+        author = "Morisaki, Soichiro and Smith, Rory and Tsukada, Leo and Sachdev, Surabhi and Stevenson, Simon and Talbot, Colm and Zimmerman, Aaron",
+        title = "{Rapid localization and inference on compact binary coalescences with the Advanced LIGO-Virgo-KAGRA gravitational-wave detector network}",
+        eprint = "2307.13380",
+        archivePrefix = "arXiv",
+        primaryClass = "gr-qc",
+        doi = "10.1103/PhysRevD.108.123040",
+        journal = "Phys. Rev. D",
+        volume = "108",
+        number = "12",
+        pages = "123040",
+        year = "2023"
+    }
+
+
 - `MBGravitationalWaveTransient`
  .. code:: bibtex
 
@@ -93,8 +108,18 @@ If you use any of the accelerated likelihoods like `ROQGravitationalWaveTransien
 
 - `RelativeBinningGravitationalWaveTransient`
  .. code:: bibtex
+    
+    @article{relbin_bilby,
+        author = "Krishna, Kruthi and Vijaykumar, Aditya and Ganguly, Apratim and Talbot, Colm and Biscoveanu, Sylvia and George, Richard N. and Williams, Natalie and Zimmerman, Aaron",
+        title = "{Accelerated parameter estimation in Bilby with relative binning}",
+        eprint = "2312.06009",
+        archivePrefix = "arXiv",
+        primaryClass = "gr-qc",
+        month = "12",
+        year = "2023"
+    }
 
-    @article{relbin_paper,
+    @article{relbin_cornish,
         author = "Cornish, Neil J.",
         title = "{Fast Fisher Matrices and Lazy Likelihoods}",
         eprint = "1007.4820",
@@ -102,10 +127,64 @@ If you use any of the accelerated likelihoods like `ROQGravitationalWaveTransien
         primaryClass = "gr-qc",
         month = "7",
         year = "2010"
-    }   
+    }
+
+    @article{relbin_zackay,
+        author = "Zackay, Barak and Dai, Liang and Venumadhav, Tejaswi",
+        title = "{Relative Binning and Fast Likelihood Evaluation for Gravitational Wave Parameter Estimation}",
+        eprint = "1806.08792",
+        archivePrefix = "arXiv",
+        primaryClass = "astro-ph.IM",
+        month = "6",
+        year = "2018"
+    } 
+
+If you use the :code:`bilby_mcmc` sampler, please additionally cite the following paper
 
-Please also ensure that you appropriately cite the following:
+ .. code:: bibtex
+
+    @article{bilby_mcmc_paper,
+        author = "Ashton, Gregory and Talbot, Colm",
+        title = "{B\,ilby-MCMC: an MCMC sampler for gravitational-wave inference}",
+        eprint = "2106.08730",
+        archivePrefix = "arXiv",
+        primaryClass = "gr-qc",
+        doi = "10.1093/mnras/stab2236",
+        journal = "Mon. Not. Roy. Astron. Soc.",
+        volume = "507",
+        number = "2",
+        pages = "2037--2051",
+        year = "2021"
+    }
 
-- The sampler(s) used for your analysis
-- Any additional package(s) that you use to generate waveforms e.g. `lalsuite` (https://doi.org/10.7935/GT1W-FZ16)
-as well as the original scientific papers for each waveform.
+Additionally, :code:`bilby` builds on a number of open-source packages. If you
+make use of this functionality in your publications, we recommend you cite them
+as requested in their associated documentation.
+
+**Samplers**
+* `cpnest <https://github.com/johnveitch/cpnest>`__
+* `dnest4 <https://github.com/eggplantbren/DNest4>`__
+* `dynesty <https://github.com/joshspeagle/dynesty>`__
+* `emcee <https://github.com/dfm/emcee>`__
+* `kombine <https://github.com/bfarr/kombine>`__
+* `nestle <https://github.com/kbarbary/nestle>`__
+* `nessai <https://github.com/mj-will/nessai>`__
+* `PyMC3 <https://github.com/pymc-devs/pymc3>`__
+* `pymultinest <https://github.com/JohannesBuchner/PyMultiNest>`__
+* `pypolychord <https://github.com/PolyChord/PolyChordLite>`__
+* `ptemcee <https://github.com/willvousden/ptemcee>`__
+* `ptmcmcsampler <https://github.com/jellis18/PTMCMCSampler>`__
+* `ultranest <https://github.com/JohannesBuchner/UltraNest>`__
+* `zeus <https://github.com/minaskar/zeus>`_
+
+
+**Gravitational-wave tools**
+
+* `gwpy <https://github.com/gwpy/gwpy>`__
+* `lalsuite <https://git.ligo.org/lscsoft/lalsuite>`__
+* `astropy <https://github.com/astropy/astropy>`__
+
+**Plotting**
+
+* `corner <https://github.com/dfm/corner.py>`__ for generating corner plot
+* `matplotlib <https://github.com/matplotlib/matplotlib>`__ for general plotting routines

docs/conf.py

@@ -66,6 +66,7 @@ extensions = [
     'sphinx.ext.autosectionlabel',
     'sphinx_tabs.tabs',
     "sphinx.ext.linkcode",
+    'myst_parser'
 ]
 autosummary_generate = True
 

docs/contribution-guide.md

+```{include} ../CONTRIBUTING.md
+```

docs/index.txt

@@ -29,6 +29,7 @@ Welcome to bilby's documentation!
    writing-documentation
    hyperparameters
    containers
+   contribution-guide
    faq
 
 

docs/plugins.txt

@@ -61,6 +61,19 @@ here, we encourage you to open a
 - This could be your sampler
 
 
+Bilby-native samplers
+---------------------
+
+Some samplers are implemented directly in :code:`bilby` and these are avertised
+under two possible names:
+
+- :code:`bilby.<sampler name>`: always available, indicates the sampler is implemented in bilby,
+- :code:`<sampler name>`: only refers to the native bilby implementation if an external plugin does not already provide this sampler.
+
+This allows for an external plugin to provide a sampler without introducing
+namespace conflicts.
+
+
 --------------------------------
 Information for bilby developers
 --------------------------------

docs/requirements.txt

@@ -5,3 +5,4 @@ numpydoc
 nbsphinx
 autodoc
 sphinx_rtd_theme
+myst_parser