From afb2d1295fadc7960de058edc12bac0a3b39421f Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Mon, 4 Apr 2022 19:23:29 +0000
Subject: [PATCH] Add missing GW docstrings

---
 bilby/gw/conversion.py         | 134 +++++++++++++++-----
 bilby/gw/cosmology.py          |  36 ++++--
 bilby/gw/prior.py              | 218 ++++++++++++++++++++++++++++-----
 bilby/gw/result.py             |   4 +
 bilby/gw/source.py             | 103 +++++++++++++++-
 bilby/gw/utils.py              |  76 +++++++++++-
 bilby/gw/waveform_generator.py |  68 +++++-----
 7 files changed, 529 insertions(+), 110 deletions(-)

diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index bdfd51c13..e811612f6 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -1,3 +1,8 @@
+"""
+A collection of functions to convert between parameters describing
+gravitational-wave sources.
+"""
+
 import os
 import sys
 import multiprocessing
@@ -24,14 +29,12 @@ def redshift_to_comoving_distance(redshift, cosmology=None):
     return cosmology.comoving_distance(redshift).value
 
 
-@np.vectorize
 def luminosity_distance_to_redshift(distance, cosmology=None):
     from astropy import units
     cosmology = get_cosmology(cosmology)
     return z_at_value(cosmology.luminosity_distance, distance * units.Mpc)
 
 
-@np.vectorize
 def comoving_distance_to_redshift(distance, cosmology=None):
     from astropy import units
     cosmology = get_cosmology(cosmology)
@@ -50,33 +53,44 @@ def luminosity_distance_to_comoving_distance(distance, cosmology=None):
     return redshift_to_comoving_distance(redshift, cosmology)
 
 
-def bilby_to_lalsimulation_spins(
-        theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
-        reference_frequency, phase):
-    if (a_1 == 0 or tilt_1 in [0, np.pi]) and (a_2 == 0 or tilt_2 in [0, np.pi]):
-        spin_1x = 0
-        spin_1y = 0
-        spin_1z = a_1 * np.cos(tilt_1)
-        spin_2x = 0
-        spin_2y = 0
-        spin_2z = a_2 * np.cos(tilt_2)
-        iota = theta_jn
-    else:
-        iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
-            transform_precessing_spins(
-                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
-                mass_2, reference_frequency, phase)
-    return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
+_cosmology_docstring = """
+Convert from {input} to {output}
+
+Parameters
+----------
+{input}: float
+    The {input} to convert.
+cosmology: astropy.cosmology.Cosmology
+    The cosmology to use for the transformation.
+    See :code:`bilby.gw.cosmology.get_cosmology` for details of how to
+    specify this.
+
+Returns
+-------
+float
+    The {output} corresponding to the provided {input}.
+"""
+
+for _func in [
+    comoving_distance_to_luminosity_distance,
+    comoving_distance_to_redshift,
+    luminosity_distance_to_comoving_distance,
+    luminosity_distance_to_redshift,
+    redshift_to_comoving_distance,
+    redshift_to_luminosity_distance,
+]:
+    input, output = _func.__name__.split("_to_")
+    _func.__doc__ = _cosmology_docstring.format(input=input, output=output)
 
 
-@np.vectorize
-def transform_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1,
-                               a_2, mass_1, mass_2, reference_frequency, phase):
+def bilby_to_lalsimulation_spins(
+    theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
+    reference_frequency, phase
+):
     """
-    Vectorized version of
-    lalsimulation.SimInspiralTransformPrecessingNewInitialConditions
+    Convert from Bilby spin parameters to lalsimulation ones.
 
-    All parameters are defined at the reference frequency
+    All parameters are defined at the reference frequency and in SI units.
 
     Parameters
     ==========
@@ -95,9 +109,9 @@ def transform_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1,
     a_2: float
         Secondary dimensionless spin magnitude
     mass_1: float
-        Primary mass _in SI units_
+        Primary mass in SI units
     mass_2: float
-        Secondary mass _in SI units_
+        Secondary mass in SI units
     reference_frequency: float
     phase: float
         Orbital phase
@@ -109,13 +123,40 @@ def transform_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1,
     spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z: float
         Cartesian spin components
     """
+    if (a_1 == 0 or tilt_1 in [0, np.pi]) and (a_2 == 0 or tilt_2 in [0, np.pi]):
+        spin_1x = 0
+        spin_1y = 0
+        spin_1z = a_1 * np.cos(tilt_1)
+        spin_2x = 0
+        spin_2y = 0
+        spin_2z = a_2 * np.cos(tilt_2)
+        iota = theta_jn
+    else:
+        from numbers import Number
+        args = (
+            theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+            mass_2, reference_frequency, phase
+        )
+        float_inputs = all([isinstance(arg, Number) for arg in args])
+        if float_inputs:
+            func = lalsim_SimInspiralTransformPrecessingNewInitialConditions
+        else:
+            func = transform_precessing_spins
+        iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = func(*args)
+    return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
 
-    iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = (
-        lalsim_SimInspiralTransformPrecessingNewInitialConditions(
-            theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
-            reference_frequency, phase))
 
-    return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
+@np.vectorize
+def transform_precessing_spins(*args):
+    """
+    Vectorized wrapper for
+    lalsimulation.SimInspiralTransformPrecessingNewInitialConditions
+
+    For detailed documentation see
+    :code:`bilby.gw.conversion.bilby_to_lalsimulation_spins`.
+    This will be removed from the public API in a future release.
+    """
+    return lalsim_SimInspiralTransformPrecessingNewInitialConditions(*args)
 
 
 def convert_to_lal_binary_black_hole_parameters(parameters):
@@ -517,8 +558,10 @@ def chirp_mass_and_mass_ratio_to_total_mass(chirp_mass, mass_ratio):
 
     Returns
     =======
-    total_mass: float
-        Total mass of the binary
+    mass_1: float
+        Mass of the heavier object
+    mass_2: float
+        Mass of the lighter object
     """
 
     with np.errstate(invalid="ignore"):
@@ -904,6 +947,29 @@ def generate_all_bns_parameters(sample, likelihood=None, priors=None, npool=1):
 
 
 def generate_specific_parameters(sample, parameters):
+    """
+    Generate a specific subset of parameters that can be generated.
+
+    Parameters
+    ----------
+    sample: dict
+        The input sample to be converted.
+    parameters: list
+        The list of parameters to return.
+
+    Returns
+    -------
+    output_sample: dict
+        The converted parameters
+
+    Notes
+    -----
+    This is _not_ an optimized function. Under the hood, it generates all
+    possible parameters and then downselects.
+
+    If the passed :code:`parameters` do not include the input parameters,
+    those will not be returned.
+    """
     updated_sample = generate_all_bns_parameters(sample=sample.copy())
     output_sample = sample.__class__()
     for key in parameters:
diff --git a/bilby/gw/cosmology.py b/bilby/gw/cosmology.py
index c05c7ada3..c7a22f004 100644
--- a/bilby/gw/cosmology.py
+++ b/bilby/gw/cosmology.py
@@ -1,3 +1,7 @@
+"""
+Wrappers to :code:`Astropy` functionality for specifying the cosmology to use.
+"""
+
 DEFAULT_COSMOLOGY = None
 COSMOLOGY = [None, str(None)]
 
@@ -11,6 +15,25 @@ def _set_default_cosmology():
 
 
 def get_cosmology(cosmology=None):
+    """
+    Get an instance of a astropy.cosmology.FLRW subclass.
+
+    To avoid repeatedly instantiating the same class, test if it is the same
+    as the last used cosmology.
+
+    Parameters
+    ==========
+    cosmology: astropy.cosmology.FLRW, str, dict
+        Description of cosmology, one of:
+            None - Use DEFAULT_COSMOLOGY
+            Instance of astropy.cosmology.FLRW subclass
+            String with name of known Astropy cosmology, e.g., "Planck13"
+
+    Returns
+    =======
+    cosmo: astropy.cosmology.FLRW
+        Cosmology instance
+    """
     from astropy import cosmology as cosmo
     _set_default_cosmology()
     if cosmology is None:
@@ -22,10 +45,8 @@ def get_cosmology(cosmology=None):
 
 def set_cosmology(cosmology=None):
     """
-    Get an instance of a astropy.cosmology.FLRW subclass.
-
-    To avoid repeatedly instantiating the same class, test if it is the same
-    as the last used cosmology.
+    Set an instance of a astropy.cosmology.FLRW subclass as the default
+    cosmology.
 
     Parameters
     ==========
@@ -36,11 +57,6 @@ def set_cosmology(cosmology=None):
             String with name of known Astropy cosmology, e.g., "Planck13"
             Dictionary with arguments required to instantiate the cosmology
             class.
-
-    Returns
-    =======
-    cosmo: astropy.cosmology.FLRW
-        Cosmology instance
     """
     from astropy import cosmology as cosmo
     _set_default_cosmology()
@@ -74,4 +90,4 @@ def z_at_value(func, fval, **kwargs):
     for detailed documentation.
     """
     from astropy.cosmology import z_at_value
-    return float(z_at_value(func=func, fval=fval, **kwargs))
+    return z_at_value(func=func, fval=fval, **kwargs).value
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 824bd4d4e..54f803401 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -92,6 +92,9 @@ def convert_to_flat_in_component_mass_prior(result, fraction=0.25):
 
 
 class Cosmological(Interped):
+    """
+    Base prior class for cosmologically aware distance parameters.
+    """
 
     @property
     def _default_args_dict(self):
@@ -105,6 +108,29 @@ class Cosmological(Interped):
 
     def __init__(self, minimum, maximum, cosmology=None, name=None,
                  latex_label=None, unit=None, boundary=None):
+        """
+
+        Parameters
+        ----------
+        minimum: float
+            The minimum value for the distribution.
+        maximum: float
+            The maximum value for the distribution.
+        cosmology: Union[None, str, dict, astropy.cosmology.FLRW]
+            The cosmology to use, see `bilby.gw.cosmology.set_cosmology`
+            for more details.
+            If :code:`None`, the default cosmology will be used.
+        name: str
+            The name of the parameter, should be one of
+            :code:`[luminosity_distance, comoving_distance, redshift]`.
+        latex_label: str
+            Latex string to use in plotting, if :code:`None`, this will
+            be set automatically using the name.
+        unit: Union[str, astropy.units.Unit]
+            The unit of the maximum and minimum values, :code:`default=Mpc`.
+        boundary: str
+            The boundary condition to apply to the prior when sampling.
+        """
         from astropy import units
         self.cosmology = get_cosmology(cosmology)
         if name not in self._default_args_dict:
@@ -210,6 +236,22 @@ class Cosmological(Interped):
             pass
 
     def get_corresponding_prior(self, name=None, unit=None):
+        """
+        Utility method to convert between equivalent redshift, comoving
+        distance, and luminosity distance priors.
+
+        Parameters
+        ----------
+        name: str
+            The name of the new prior, this should be one of
+            :code:`[luminosity_distance, comoving_distance, redshift]`.
+        unit: Union[str, astropy.units.Unit]
+            The unit of the output prior.
+
+        Returns
+        -------
+        The corresponding prior.
+        """
         subclass_args = infer_args_from_method(self.__init__)
         args_dict = {key: getattr(self, key) for key in subclass_args}
         self._convert_to(new=name, args_dict=args_dict)
@@ -267,6 +309,16 @@ class Cosmological(Interped):
 
 
 class UniformComovingVolume(Cosmological):
+    r"""
+    Prior that is uniform in the comoving volume.
+
+    Note that this does not account for time dilation in the source frame
+    use `bilby.gw.prior.UniformSourceFrame` to include that effect.
+
+    .. math::
+
+        p(z) \propto \frac{d V_{c}}{dz}
+    """
 
     def _get_redshift_arrays(self):
         zs = np.linspace(self._minimum['redshift'] * 0.99,
@@ -280,8 +332,11 @@ class UniformSourceFrame(Cosmological):
     Prior for redshift which is uniform in comoving volume and source frame
     time.
 
-    For redshift this is p(z) \propto dVc/dz 1 / (1 + z), where the extra 1+z
-    is due to doppler shifting of the source frame time.
+    .. math::
+
+        p(z) \propto \frac{1}{1 + z}\frac{dV_{c}}{dz}
+
+    The extra :math:`1+z` is due to doppler shifting of the source frame time.
     """
 
     def _get_redshift_arrays(self):
@@ -292,15 +347,25 @@ class UniformSourceFrame(Cosmological):
 
 
 class UniformInComponentsChirpMass(PowerLaw):
+    r"""
+    Prior distribution for chirp mass which is uniform in component masses.
+
+    This is useful when chirp mass and mass ratio are sampled while the
+    prior is uniform in component masses.
+
+    .. math::
+
+        p({\cal M}) \propto {\cal M}
+
+    Notes
+    -----
+    This prior is intended to be used in conjunction with the corresponding
+    :code:`bilby.gw.prior.UniformInComponentsMassRatio`.
+    """
 
     def __init__(self, minimum, maximum, name='chirp_mass',
                  latex_label=r'$\mathcal{M}$', unit=None, boundary=None):
         """
-        Prior distribution for chirp mass which is uniform in component masses.
-
-        This is useful when chirp mass and mass ratio are sampled while the
-        prior is uniform in component masses.
-
         Parameters
         ==========
         minimum : float
@@ -333,15 +398,25 @@ class WrappedInterp1d(interp1d):
 
 
 class UniformInComponentsMassRatio(Prior):
+    r"""
+    Prior distribution for chirp mass which is uniform in component masses.
+
+    This is useful when chirp mass and mass ratio are sampled while the
+    prior is uniform in component masses.
+
+    .. math::
+
+        p(q) \propto \frac{(1 + q)^{2/5}}{q^{6/5}}
+
+    Notes
+    -----
+    This prior is intended to be used in conjunction with the corresponding
+    :code:`bilby.gw.prior.UniformInComponentsChirpMass`.
+    """
 
     def __init__(self, minimum, maximum, name='mass_ratio', latex_label='$q$',
                  unit=None, boundary=None):
         """
-        Prior distribution for mass ratio which is uniform in component masses.
-
-        This is useful when chirp mass and mass ratio are sampled while the
-        prior is uniform in component masses.
-
         Parameters
         ==========
         minimum : float
@@ -388,20 +463,26 @@ class UniformInComponentsMassRatio(Prior):
 
 
 class AlignedSpin(Interped):
+    r"""
+    Prior distribution for the aligned (z) component of the spin.
 
-    def __init__(self, a_prior=Uniform(0, 1), z_prior=Uniform(-1, 1),
-                 name=None, latex_label=None, unit=None, boundary=None,
-                 minimum=np.nan, maximum=np.nan):
-        """
-        Prior distribution for the aligned (z) component of the spin.
+    This takes prior distributions for the magnitude and cosine of the tilt
+    and forms a compound prior using a numerical convolution integral.
 
-        This takes prior distributions for the magnitude and cosine of the tilt
-        and forms a compound prior.
+    .. math::
 
-        This is useful when using aligned-spin only waveform approximants.
+        \pi(\chi) = \int_{0}^{1} da \int_{-1}^{1} d\cos\theta
+        \pi(a) \pi(\cos\theta) \delta(\chi - a \cos\theta)
 
-        This is an extension of e.g., (A7) of https://arxiv.org/abs/1805.10457.
+    This is useful when using aligned-spin only waveform approximants.
 
+    This is an extension of e.g., (A7) of https://arxiv.org/abs/1805.10457.
+    """
+
+    def __init__(self, a_prior=Uniform(0, 1), z_prior=Uniform(-1, 1),
+                 name=None, latex_label=None, unit=None, boundary=None,
+                 minimum=np.nan, maximum=np.nan):
+        """
         Parameters
         ==========
         a_prior: Prior
@@ -437,8 +518,9 @@ class ConditionalChiUniformSpinMagnitude(ConditionalLogUniform):
     if the distribution of spin orientations is isotropic.
 
     .. math::
-        p(a) &= \frac{1}{a_{\max}}
-        p(\chi) &= - \frac{1}{2 a_{\max}} \ln(|\chi|)
+
+        p(a) &= \frac{1}{a_{\max}} \\
+        p(\chi) &= - \frac{\ln(|\chi|)}{2 a_{\max}}  \\
         p(a | \chi) &\propto \frac{1}{a}
     """
 
@@ -471,9 +553,10 @@ class ConditionalChiInPlane(ConditionalBasePrior):
     if the distribution of spin orientations is isotropic.
 
     .. math::
-        p(a) &= \frac{1}{a_{\max}}
-        p(\chi_\perp) = 2 N \chi_\perp / (\chi ** 2 + \chi_\perp ** 2)
-        N^{-1} &= 2 \ln(a_\max / |\chi|)
+
+        p(a) &= \frac{1}{a_{\max}} \\
+        p(\chi_\perp) &= \frac{2 N \chi_\perp}{\chi^2 + \chi_\perp^2} \\
+        N^{-1} &= 2 \ln\left(\frac{a_\max}{|\chi|}\right)
     """
 
     def __init__(self, minimum, maximum, name, latex_label=None, unit=None, boundary=None):
@@ -635,6 +718,18 @@ class CBCPriorDict(ConditionalPriorDict):
 
     @property
     def minimum_component_mass(self):
+        """
+        The minimum component mass allowed for the prior dictionary.
+
+        This property requires either:
+        * a prior for :code:`mass_2`
+        * priors for :code:`chirp_mass` and :code:`mass_ratio`
+
+        Returns
+        -------
+        mass_2: float
+            The minimum allowed component mass.
+        """
         if "mass_2" in self:
             return self["mass_2"].minimum
         if "chirp_mass" in self and "mass_ratio" in self:
@@ -870,13 +965,15 @@ class BNSPriorDict(CBCPriorDict):
         Default parameter conversion function for BNS signals.
 
         This generates:
-        - the parameters passed to source.lal_binary_neutron_star
-        - all mass parameters
-        - all tidal parameters
+
+        * the parameters passed to source.lal_binary_neutron_star
+        * all mass parameters
+        * all tidal parameters
 
         It does not generate:
-        - component spins
-        - source-frame parameters
+
+        * component spins
+        * source-frame parameters
 
         Parameters
         ==========
@@ -966,9 +1063,15 @@ Prior._default_latex_labels = {
 
 
 class CalibrationPriorDict(PriorDict):
+    """
+    Prior dictionary class for calibration parameters.
+    This has methods for simplifying the creation of priors for the large
+    numbers of parameters used with the spline model.
+    """
 
     def __init__(self, dictionary=None, filename=None):
-        """ Initialises a Prior set for Binary Black holes
+        """
+        Initialises a Prior dictionary for calibration parameters
 
         Parameters
         ==========
@@ -1140,6 +1243,35 @@ class CalibrationPriorDict(PriorDict):
 
 
 def secondary_mass_condition_function(reference_params, mass_1):
+    """
+    Condition function to use for a prior that is conditional on the value
+    of the primary mass, for example, a prior on the secondary mass that is
+    bounded by the primary mass.
+
+    .. code-block:: python
+
+        import bilby
+        priors = bilby.gw.prior.CBCPriorDict()
+        priors["mass_1"] = bilby.core.prior.Uniform(5, 50)
+        priors["mass_2"] = bilby.core.prior.ConditionalUniform(
+            condition_func=bilby.gw.prior.secondary_mass_condition_function,
+            minimum=5,
+            maximum=50,
+        )
+
+    Parameters
+    ----------
+    reference_params: dict
+        Reference parameter dictionary, this should contain a "minimum"
+        attribute.
+    mass_1: float
+        The primary mass to use as the upper limit for the prior.
+
+    Returns
+    -------
+    dict:
+        Updated prior limits given the provided primary mass.
+    """
     return dict(minimum=reference_params['minimum'], maximum=mass_1)
 
 
@@ -1502,7 +1634,29 @@ class HealPixMapPriorDist(BaseJointPriorDist):
 
 
 class HealPixPrior(JointPrior):
+    """
+    A prior distribution that follows a user-provided HealPix map for one
+    parameter.
+
+    See :code:`bilby.gw.prior.HealPixMapPriorDist` for more details of how to
+    instantiate the prior.
+    """
     def __init__(self, dist, name=None, latex_label=None, unit=None):
+        """
+
+        Parameters
+        ----------
+        dist: bilby.gw.prior.HealPixMapPriorDist
+            The base joint probability.
+        name: str
+            The name of the parameter, it should be contained in the map.
+            One of ["ra", "dec", "luminosity_distance"].
+        latex_label: str
+            Latex label used for plotting, will be read from default values if
+            not provided.
+        unit: str
+            The unit of the parameter.
+        """
         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)
diff --git a/bilby/gw/result.py b/bilby/gw/result.py
index 5c95fe0f1..6dee3d6bd 100644
--- a/bilby/gw/result.py
+++ b/bilby/gw/result.py
@@ -14,6 +14,10 @@ from .detector import get_empty_interferometer, Interferometer
 
 
 class CompactBinaryCoalescenceResult(CoreResult):
+    """
+    Result class with additional methods and attributes specific to analyses
+    of compact binaries.
+    """
     def __init__(self, **kwargs):
         super(CompactBinaryCoalescenceResult, self).__init__(**kwargs)
 
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index f2aea7409..4e044a0b3 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -828,6 +828,44 @@ def _base_waveform_frequency_sequence(
 
 
 def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
+    r"""
+    A frequency-domain sine-Gaussian burst source model.
+
+    .. math::
+
+        \tau &= \frac{Q}{\sqrt{2}\pi f_{0}} \\
+        \alpha &= \frac{Q}{4\sqrt{\pi} f_{0}} \\
+        h_{+} &=
+            \frac{h_{\rm rss}\sqrt{\pi}\tau}{2\sqrt{\alpha (1 + e^{-Q^2})}}
+            \left(
+                e^{-\pi^2 \tau^2 (f + f_{0})^2}
+                + e^{-\pi^2 \tau^2 (f - f_{0})^2}
+            \right) \\
+        h_{\times} &=
+            \frac{i h_{\rm rss}\sqrt{\pi}\tau}{2\sqrt{\alpha (1 - e^{-Q^2})}}
+            \left(
+                e^{-\pi^2 \tau^2 (f + f_{0})^2}
+                - e^{-\pi^2 \tau^2 (f - f_{0})^2}
+            \right)
+
+    Parameters
+    ----------
+    frequency_array: array-like
+        The frequencies at which to compute the model.
+    hrss: float
+        The amplitude of the signal.
+    Q: float
+        The quality factor of the burst, determines the decay time.
+    frequency: float
+        The peak frequency of the burst.
+    kwargs: dict
+        UNUSED
+
+    Returns
+    -------
+    dict:
+        Dictionary containing the plus and cross components of the strain.
+    """
     tau = Q / (np.sqrt(2.0) * np.pi * frequency)
     temp = Q / (4.0 * np.sqrt(np.pi) * frequency)
     fm = frequency_array - frequency
@@ -848,7 +886,33 @@ def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
 
 def supernova(
         frequency_array, realPCs, imagPCs, file_path, luminosity_distance, **kwargs):
-    """ A supernova NR simulation for injections """
+    """
+    A source model that reads a simulation from a text file.
+
+    This was originally intended for use with supernova simulations, but can
+    be applied to any source class.
+
+    Parameters
+    ----------
+    frequency_array: array-like
+        Unused
+    realPCs: UNUSED
+    imagPCs: UNUSED
+    file_path: str
+        Path to the file containing the NR simulation. The format of this file
+        should be readable by :code:`numpy.loadtxt` and have four columns
+        containing the real and imaginary components of the plus and cross
+        polarizations.
+    luminosity_distance: float
+        The distance to the source in kpc, this scales the amplitude of the
+        signal. The simulation is assumed to be at 10kpc.
+    kwargs: UNUSED
+
+    Returns
+    -------
+    dict:
+        A dictionary containing the plus and cross components of the signal.
+    """
 
     realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(
         file_path, usecols=(0, 1, 2, 3), unpack=True)
@@ -864,7 +928,42 @@ def supernova(
 def supernova_pca_model(
         frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5,
         luminosity_distance, **kwargs):
-    """ Supernova signal model """
+    r"""
+    Signal model based on a five-component principal component decomposition
+    of a model.
+
+    While this was initially intended for modelling supernova signal, it is
+    applicable to any situation using such a principal component decomposition.
+
+    .. math::
+
+        h_{A} = \frac{10^{-22}}{d_{L}} \sum_{i=1}^{5} c_{i} h_{i}
+
+    Parameters
+    ----------
+    frequency_array: UNUSED
+    pc_coeff1: float
+        The first principal component coefficient.
+    pc_coeff2: float
+        The second principal component coefficient.
+    pc_coeff3: float
+        The third principal component coefficient.
+    pc_coeff4: float
+        The fourth principal component coefficient.
+    pc_coeff5: float
+        The fifth principal component coefficient.
+    luminosity_distance: float
+        The distance to the source, the amplitude is scaled such that the
+        amplitude at 10 kpc is 1e-23.
+    kwargs: dict
+        Dictionary containing numpy arrays with the real and imaginary
+        components of the principal component decomposition.
+
+    Returns
+    -------
+    dict:
+        The plus and cross polarizations of the signal
+    """
 
     realPCs = kwargs['realPCs']
     imagPCs = kwargs['imagPCs']
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 1d87e971a..20b483706 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -251,6 +251,9 @@ def matched_filter_snr(signal, frequency_domain_strain, power_spectral_density,
 
 def optimal_snr_squared(signal, power_spectral_density, duration):
     """
+    Compute the square of the optimal matched filter SNR for the provided
+    signal.
+
 
     Parameters
     ==========
@@ -271,6 +274,34 @@ def optimal_snr_squared(signal, power_spectral_density, duration):
 
 def overlap(signal_a, signal_b, power_spectral_density=None, delta_frequency=None,
             lower_cut_off=None, upper_cut_off=None, norm_a=None, norm_b=None):
+    r"""
+    Compute the overlap between two signals
+
+    .. math::
+
+        {\cal O} = \frac{4 \Delta f}{N_{a} N_{b}} \sum_{i} \frac{h^{*}_{a,i} h_{b,i}}{S_{i}}
+
+    Parameters
+    ----------
+    signal_a: array-like
+    signal_b: array-like
+    power_spectral_density : array-like
+    delta_frequency: float
+        Frequency spacing of the signals
+    lower_cut_off: float
+        Minimum frequency for the integral
+    upper_cut_off: float
+        Maximum frequency for the integral
+    norm_a: float
+        Normalizing factor for signal_a
+    norm_b: float
+        Normalizing factor for signal_b
+
+    Returns
+    -------
+    float
+        The overlap
+    """
     low_index = int(lower_cut_off / delta_frequency)
     up_index = int(upper_cut_off / delta_frequency)
     integrand = np.conj(signal_a) * signal_b
@@ -536,6 +567,34 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=0
 
 
 def get_gracedb(gracedb, outdir, duration, calibration, detectors, query_types=None, server=None):
+    """
+    Download information about a trigger from GraceDb and create cache files
+    for finding gravitational-wave strain data.
+
+    Parameters
+    ----------
+    gracedb: str
+        The GraceDb ID of the trigger.
+    outdir: str
+        Directory to write output.
+    duration: int
+        Duration of data to find about the trigger time (units: :code:`s`).
+    calibration: int
+        Calibration label of the data, should be one of :code:`1, 2`.
+    detectors: list
+        The detectors to look for data for.
+    query_types: str
+        The LDR query type
+    server: str
+        The LIGO datafind server to look for data on.
+
+    Returns
+    -------
+    candidate: dict
+        Information downloaded from GraceDb about the trigger.
+    cache_files: list
+        List of cache filenames, one per interferometer.
+    """
     candidate = gracedb_to_json(gracedb, outdir=outdir)
     trigger_time = candidate['gpstime']
     gps_start_time = trigger_time - duration
@@ -608,7 +667,7 @@ def gw_data_find(observatory, gps_start_time, duration, calibration,
         The start time in gps to look for data
     duration: int
         The duration (integer) in s
-    calibrartion: int {1, 2}
+    calibration: int {1, 2}
         Use C01 or C02 calibration
     outdir: string
         A path to the directory where output is stored
@@ -966,6 +1025,21 @@ def plot_spline_pos(log_freqs, samples, nfreqs=100, level=0.9, color='k', label=
 
 
 def greenwich_mean_sidereal_time(time):
+    """
+    Compute the greenwich mean sidereal time from a GPS time.
+
+    This is just a wrapper around :code:`lal.GreenwichMeanSiderealTime` .
+
+    Parameters
+    ----------
+    time: float
+        The GPS time to convert.
+
+    Returns
+    -------
+    float
+        The sidereal time.
+    """
     from lal import GreenwichMeanSiderealTime
     time = float(time)
     return GreenwichMeanSiderealTime(time)
diff --git a/bilby/gw/waveform_generator.py b/bilby/gw/waveform_generator.py
index b068156fe..f4cd5a6fe 100644
--- a/bilby/gw/waveform_generator.py
+++ b/bilby/gw/waveform_generator.py
@@ -7,6 +7,11 @@ from .conversion import convert_to_lal_binary_black_hole_parameters
 
 
 class WaveformGenerator(object):
+    """
+    The base waveform generator class.
+
+    Waveform generators provide a unified method to call disparate source models.
+    """
 
     duration = PropertyAccessor('_times_and_frequencies', 'duration')
     sampling_frequency = PropertyAccessor('_times_and_frequencies', 'sampling_frequency')
@@ -18,37 +23,38 @@ class WaveformGenerator(object):
                  time_domain_source_model=None, parameters=None,
                  parameter_conversion=None,
                  waveform_arguments=None):
-        """ A waveform generator
-
-    Parameters
-    ==========
-    sampling_frequency: float, optional
-        The sampling frequency
-    duration: float, optional
-        Time duration of data
-    start_time: float, optional
-        Starting time of the time array
-    frequency_domain_source_model: func, optional
-        A python function taking some arguments and returning the frequency
-        domain strain. Note the first argument must be the frequencies at
-        which to compute the strain
-    time_domain_source_model: func, optional
-        A python function taking some arguments and returning the time
-        domain strain. Note the first argument must be the times at
-        which to compute the strain
-    parameters: dict, optional
-        Initial values for the parameters
-    parameter_conversion: func, optional
-        Function to convert from sampled parameters to parameters of the
-        waveform generator. Default value is the identity, i.e. it leaves
-        the parameters unaffected.
-    waveform_arguments: dict, optional
-        A dictionary of fixed keyword arguments to pass to either
-        `frequency_domain_source_model` or `time_domain_source_model`.
-
-        Note: the arguments of frequency_domain_source_model (except the first,
-        which is the frequencies at which to compute the strain) will be added to
-        the WaveformGenerator object and initialised to `None`.
+        """
+        The base waveform generator class.
+
+        Parameters
+        ==========
+        sampling_frequency: float, optional
+            The sampling frequency
+        duration: float, optional
+            Time duration of data
+        start_time: float, optional
+            Starting time of the time array
+        frequency_domain_source_model: func, optional
+            A python function taking some arguments and returning the frequency
+            domain strain. Note the first argument must be the frequencies at
+            which to compute the strain
+        time_domain_source_model: func, optional
+            A python function taking some arguments and returning the time
+            domain strain. Note the first argument must be the times at
+            which to compute the strain
+        parameters: dict, optional
+            Initial values for the parameters
+        parameter_conversion: func, optional
+            Function to convert from sampled parameters to parameters of the
+            waveform generator. Default value is the identity, i.e. it leaves
+            the parameters unaffected.
+        waveform_arguments: dict, optional
+            A dictionary of fixed keyword arguments to pass to either
+            `frequency_domain_source_model` or `time_domain_source_model`.
+
+            Note: the arguments of frequency_domain_source_model (except the first,
+            which is the frequencies at which to compute the strain) will be added to
+            the WaveformGenerator object and initialised to `None`.
 
         """
         self._times_and_frequencies = CoupledTimeAndFrequencySeries(duration=duration,
-- 
GitLab