diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index b4a090f..41e71af 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -16,12 +16,12 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies (Linux)
diff --git a/flexknot/__init__.py b/flexknot/__init__.py
index ad1a4a7..288973b 100644
--- a/flexknot/__init__.py
+++ b/flexknot/__init__.py
@@ -1,5 +1,5 @@
 """
-# Flex-Knot
+# Flex-Knot.
 
 This repo contains flex-knots and associated likelihoods used in my
 toy_sine project, and likelihoods associated with them.
@@ -10,3 +10,9 @@
 from flexknot.core import AdaptiveKnot, FlexKnot
 from flexknot.likelihoods import FlexKnotLikelihood
 from flexknot.priors import AdaptiveKnotPrior, FlexKnotPrior
+
+__all__ = [
+    "AdaptiveKnot", "FlexKnot",
+    "FlexKnotLikelihood",
+    "AdaptiveKnotPrior", "FlexKnotPrior",
+]
diff --git a/flexknot/core.py b/flexknot/core.py
index 9611cd6..d85e3e5 100644
--- a/flexknot/core.py
+++ b/flexknot/core.py
@@ -1,16 +1,14 @@
 """
-Linear INterpolation Functions.
+Flex-knot.
 
-theta refers to the full set of parameters for an adaptive linear interpolation model,
-[n, y0, x1, y1, x2, y2, ..., x_n, y_n, y_n+1],
-where n is the greatest allowed value of ceil(n).
-
-The reason for the interleaving of x and y is it avoids the need to know n.
+x and y are interleaved so that N does not need to be provided for
+a non-adaptive flex-knot.
 """
+
 import numpy as np
 from scipy.integrate import quad
 
-from flexknot.helper_functions import (
+from flexknot.utils import (
     get_theta_n,
     get_x_nodes_from_theta,
     get_y_nodes_from_theta,
@@ -24,10 +22,14 @@ class FlexKnot:
     x_min: float
     x_max: float > x_min
 
-    Returns:
-    flexknot(x, theta)
+    Returns
+    -------
+    flexknot(x, theta): callable
+
+    theta has format
+    [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)]
+    for N nodes.
 
-    theta in format [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] for N nodes.
     """
 
     def __init__(self, x_min, x_max):
@@ -36,14 +38,26 @@ def __init__(self, x_min, x_max):
 
     def __call__(self, x, theta):
         """
-        Flex-knot with end nodes at x_min and x_max
+        Flex-knot with end nodes at x_min and x_max.
+
+        For N nodes:
+        theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)].
+
+        y0 and y_(N-1) are the y values at x_min and x_max respecively.
+
+        If theta only contains a single element, the flex-knot is constant.
+        If theta is empty, the flex-knot is constant at -1 (cosmology!)
+
+        Parameters
+        ----------
+        x : float or array-like
 
-        theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] for N nodes.
+        theta : array-like
 
-        y0 and y_(N-1) are the y values corresponding to x_min and x_max respecively.
+        Returns
+        -------
+        float or array-like
 
-        If theta only contains a single element, the flex-knot is constant at that value.
-        If theta is empty, the flex-knot if constant at -1 (cosmology!)
         """
         if 0 == len(theta):
             return np.full_like(x, -1)
@@ -63,25 +77,31 @@ def __call__(self, x, theta):
 
     def area(self, theta0, theta1):
         """
-        Calculate the area between the flex-knot with parameters
-        theta_0 and theta_1.
+        Area between two flex-knots.
+
+        Parameters
+        ----------
+        theta0 : array-like
+
+        theta1 : array-like
         """
         return quad(lambda x: np.abs(self(x, theta0)-self(x, theta1)),
-                self.x_min, self.x_max)[0] / (self.x_max - self.x_min)
+                    self.x_min, self.x_max)[0] / (self.x_max - self.x_min)
 
 
 class AdaptiveKnot(FlexKnot):
     """
-    Adaptive flex-knot which allows the number of parameters being used to vary.
+    Adaptive flex-knot which allows the number of knots to vary.
 
     x_min: float
     x_max: float > x_min
 
-    Returns:
-    adaptive_flexknot(x, theta)
+    Returns
+    -------
+    adaptive_flexknot(x, theta): callable
 
-    The first element of theta is N; floor(N)-2 is number of interior nodes used in
-    the linear interpolation model.
+    The first element of theta is N; floor(N)-2 is number of interior nodes
+    used by the flexknot.
 
     theta = [N, y0, x1, y1, x2, y2, ..., x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
     where Nmax is the greatest allowed value of floor(N).
@@ -92,14 +112,29 @@ class AdaptiveKnot(FlexKnot):
 
     def __call__(self, x, theta):
         """
-        The first element of theta is N; floor(N)-2 is number of interior nodes used in
-        the linear interpolation model. This is then used to select the
+        Adaptive flex-knot with end nodes at x_min and x_max.
+
+        The first element of theta is N; floor(N)-2 is number of interior nodes
+        used in the flexknot. This is then used to select the
         appropriate other elements of params to pass to flexknot()
 
-        theta = [N, y0, x1, y1, x2, y2, ..., x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
+        For a maximum of Nmax nodes:
+        theta = [N, y0, x1, y1, x2, y2, ...,
+                 x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
         where Nmax is the greatest allowed value of floor(N).
 
         if floor(N) = 1, the flex-knot is constant at theta[-1] = y_(Nmax-1).
         if floor(N) = 0, the flex-knot is constant at -1 (cosmology!)
+
+        Parameters
+        ----------
+        x : float or array-like
+
+        theta : array-like
+
+        Returns
+        -------
+        float or array-like
+
         """
         return super().__call__(x, get_theta_n(theta))
diff --git a/flexknot/helper_functions.py b/flexknot/helper_functions.py
deleted file mode 100644
index 0ee5c98..0000000
--- a/flexknot/helper_functions.py
+++ /dev/null
@@ -1,93 +0,0 @@
-"""
-Helpful functions for dealing with the interleaved arrangement of theta.
-"""
-
-import numpy as np
-
-
-def validate_theta(theta, adaptive):
-    """
-    Check that theta contains an odd/even number of elements for an adaptive/
-    non-adaptive flex-knot, and in the adaptive case that floor(theta[0])-2 isn't
-    greater than the provided number of internal nodes.
-    """
-    if adaptive:
-        if len(theta) % 2 != 1:
-            raise ValueError(
-                "theta must contain an odd number of elements for an adaptive flex-knot."
-            )
-        if np.floor(theta[0]) > (len(theta) + 1) / 2:
-            raise ValueError("n = ceil(theta[0]) exceeds the number of internal nodes.")
-
-    if not adaptive:
-        if len(theta) > 1 and len(theta) % 2 != 0:
-            raise ValueError(
-                "theta must contain an even number of elements for a non-adaptive flex-knot."
-            )
-
-
-def get_theta_n(theta):
-    """
-    Extracts the first n parameters from
-
-    theta = [N, y0, x1, y1, x2, y2, ..., x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)]
-
-    where Nmax is the maximum value of floor(N).
-
-    returns theta_n = [y0, x1, y1, x2, y2, ..., x_floor(N)-2, y_floor(N)-2, y_(Nmax-1)]
-    """
-    validate_theta(theta, adaptive=True)
-
-    n = np.floor(theta[0]).astype(int) - 2
-    # w = -1 case, return empty
-    if -2 == n:
-        return np.array([])
-    theta_n = np.concatenate(
-        (
-            theta[1 : 2 * n + 2],  # y0 and internal x and y
-            theta[-1:],  # y end node
-        )
-    )
-    return theta_n
-
-
-def create_theta(x_nodes, y_nodes):
-    """
-    Takes x_nodes = [x1, ... x_(N-2)] and y_nodes = [y0, y1, ..., y_(N-2), y_(N-1)] to
-    return theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)], where
-    N is the number of nodes.
-    """
-    if len(y_nodes) == 1:
-        return y_nodes
-    if len(x_nodes) + 2 != len(y_nodes):
-        raise ValueError("y_nodes must have exactly two more elements than x_nodes")
-    n = len(x_nodes)
-    theta = np.zeros(len(x_nodes) + len(y_nodes))
-    theta[1 : 2 * n + 1 : 2] = x_nodes
-    theta[0 : 2 * n + 2 : 2] = y_nodes[:-1]
-    theta[-1] = y_nodes[-1]
-    return theta
-
-
-def get_x_nodes_from_theta(theta, adaptive):
-    """
-    Takes theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] to return
-    x_nodes = [x1, ... x_(N-2)], where N is the number of nodes.
-    """
-    validate_theta(theta, adaptive)
-    if adaptive:
-        theta = get_theta_n(theta)
-    n = len(theta) // 2 - 1
-    return theta[1 : 2 * n + 1 : 2]
-
-
-def get_y_nodes_from_theta(theta, adaptive):
-    """
-    Takes theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] to return
-    y_nodes = [y0, y1, ..., y_(N-2), y_(N-1)], where N is the number of nodes.
-    """
-    validate_theta(theta, adaptive)
-    if adaptive:
-        theta = get_theta_n(theta)
-    n = len(theta) // 2 - 1
-    return np.concatenate((theta[0 : 2 * n + 2 : 2], theta[-1:]))
diff --git a/flexknot/likelihoods.py b/flexknot/likelihoods.py
index 0a5ea30..3c09b47 100644
--- a/flexknot/likelihoods.py
+++ b/flexknot/likelihoods.py
@@ -1,26 +1,24 @@
-"""
-Likelihoods using flex-knots.
-"""
+"""Likelihoods using flex-knots."""
 
 import numpy as np
 from scipy.special import erf, logsumexp
-from flexknot.helper_functions import (
-    get_x_nodes_from_theta,
-    get_theta_n,
-)
+from flexknot.utils import get_x_nodes_from_theta
+
 from flexknot.core import AdaptiveKnot, FlexKnot
 
 
 class FlexKnotLikelihood:
     """
-    Likelihood for a flex-knot, relative to data described by xs, ys, and sigma.
+    Likelihood for a flex-knot, with data described by xs, ys, and sigma.
 
-    sigma is either sigma_y, [sigma_x, sigma_y], [sigma_ys] or [[sigma_xs], [sigma_ys]].
+    sigma is either sigma_y, [sigma_x, sigma_y], [sigma_ys],
+    or [[sigma_xs], [sigma_ys]].
 
-    (obviously the middle two are degenerate when len(sigma) = 2, in which case
+    (Obviously the middle two are degenerate when len(sigma) = 2, in which case
     [sigma_x, sigma_y] is assumed.)
 
-    Returns likelihood(theta) -> log(L), [] where [] is the (lack of) derived parameters.
+    Returns likelihood(theta) -> log(L), [] where [] is the (lack of)
+    derived parameters.
     """
 
     def __init__(self, x_min, x_max, xs, ys, sigma, adaptive):
@@ -30,28 +28,56 @@ def __init__(self, x_min, x_max, xs, ys, sigma, adaptive):
 
     def __call__(self, theta):
         """
-        Likelihood relative to a flex-knot with parameters theta.
+        Likelihood of the data being described by flex-knot(theta).
 
-        If self.adaptive = True, the first element of theta is N; floor(N) is the number of
-        nodes used to calculate the likelihood.
+        If self.adaptive = True, the first element of theta is N; floor(N) is
+        the number of nodes in the flex-knot.
 
-        theta = [N, y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] for N nodes.
+        For N nodes,
+        theta = [N, y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)].
 
         Otherwise, theta is the same but without N.
 
         theta = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)].
+
+        Parameters
+        ----------
+        theta : array-like
+
+        Returns
+        -------
+        tuple(float, [])
+
         """
         return self._likelihood_function(theta)
 
 
 def create_likelihood_function(x_min, x_max, xs, ys, sigma, adaptive):
     """
-    Creates a likelihood function for a flex-knot, relative to data descrived by xs, ys, and sigma.
+    Create a likelihood function for a flex-knot, for data xs, ys, sigma.
 
-    sigma is either sigma_y, [sigma_x, sigma_y], [sigma_ys] or [[sigma_xs], [sigma_ys]].
+    sigma is either sigma_y, [sigma_x, sigma_y], [sigma_ys],
+    or [[sigma_xs], [sigma_ys]].
 
-    (obviously the middle two are degenerate when len(sigma) = 2, in which case
+    (Obviously the middle two are degenerate when len(sigma) = 2, in which case
     [sigma_x, sigma_y] is assumed.)
+
+    Returns likelihood(theta) -> log(L), [] where [] is the (lack of)
+    derived parameters.
+
+    Parameters
+    ----------
+    x_min : float
+    x_max : float > x_min
+    xs : array-like
+    ys : array-like
+    sigma : float or array-like
+    adaptive : bool
+
+    Returns
+    -------
+    likelihood(theta) -> log(L), []
+
     """
     LOG_2_SQRT_2πλ = np.log(2) + 0.5 * np.log(2 * np.pi * (x_max - x_min))
 
@@ -79,7 +105,9 @@ def xy_errors_likelihood(theta):
             x_nodes = np.concatenate(
                 ([x_min], get_x_nodes_from_theta(theta, adaptive), [x_max])
             )
-            # use flex-knots to get y nodes, as this is simplest way of dealing with N=0 or 1
+            # use flex-knots to get y nodes, as this is
+            # the simplest way to deal with N=0 or 1
+
             y_nodes = flexknot(x_nodes, theta)
 
             ms = (y_nodes[1:] - y_nodes[:-1]) / (x_nodes[1:] - x_nodes[:-1])
@@ -92,8 +120,9 @@ def xy_errors_likelihood(theta):
             beta = (xs * var_y + (delta * ms).T * var_x).T / q
             gamma = (np.outer(xs, ms) - delta) ** 2 / 2 / q
 
-            t_minus = (np.sqrt(q / 2).T / (sigma_x * sigma_y)).T * (x_nodes[:-1] - beta)
-            t_plus = (np.sqrt(q / 2).T / (sigma_x * sigma_y)).T * (x_nodes[1:] - beta)
+            t = (np.sqrt(q / 2).T / (sigma_x * sigma_y)).T
+            t_minus = t * (x_nodes[:-1] - beta)
+            t_plus = t * (x_nodes[1:] - beta)
 
             logL = -len(xs) * LOG_2_SQRT_2πλ
             logL += np.sum(
diff --git a/flexknot/priors.py b/flexknot/priors.py
index 0fc6226..f67c195 100644
--- a/flexknot/priors.py
+++ b/flexknot/priors.py
@@ -8,18 +8,15 @@
 """
 import numpy as np
 from pypolychord.priors import UniformPrior, SortedUniformPrior
-from flexknot.helper_functions import (
+from flexknot.utils import (
     create_theta,
-    get_theta_n,
     get_x_nodes_from_theta,
     get_y_nodes_from_theta,
 )
 
 
 class FlexKnotPrior(UniformPrior):
-    """
-    Interleaved uniform and sorted uniform priors appropriate for a flex-knot.
-    """
+    """Interleaved uniform and sorted uniform priors for a flex-knot."""
 
     def __init__(self, x_min, x_max, y_min, y_max):
         self._x_prior = SortedUniformPrior(x_min, x_max)
@@ -29,10 +26,22 @@ def __call__(self, hypercube):
         """
         Prior for flex-knot.
 
-        hypercube = [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)] for N nodes.
+        For N nodes:
+        hypercube -> [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)].
+
+        Parameters
+        ----------
+        hypercube : array-like of Uniform(0, 1).
+
+        Returns
+        -------
+        theta : array-like of SortedUniform(x_min, x_max)
+        and Uniform(y_min, y_max).
+
         """
         if len(hypercube) > 2:
-            _x_prior = self._x_prior(get_x_nodes_from_theta(hypercube, adaptive=False))
+            _x_prior = self._x_prior(get_x_nodes_from_theta(hypercube,
+                                                            adaptive=False))
         else:
             _x_prior = np.array([])
         return create_theta(
@@ -45,7 +54,8 @@ class AdaptiveKnotPrior(FlexKnotPrior):
     """
     Interleaved uniform and sorted uniform priors appropriate for a flex-knot.
 
-    N_max: int is the maximum number of nodes to use with an interactive flex-knot.
+    N_max: int
+    The maximum number of nodes to use with an adaptive flex-knot.
     """
 
     def __init__(self, x_min, x_max, y_min, y_max, N_min, N_max):
@@ -60,7 +70,7 @@ def __init__(self, x_min, x_max, y_min, y_max, N_min, N_max):
             lambda hypercube_x: np.concatenate(
                 (
                     self.__used_x_prior(hypercube_x[: self.__n_x_nodes]),
-                    self.__unused_x_prior(hypercube_x[self.__n_x_nodes :]),
+                    self.__unused_x_prior(hypercube_x[self.__n_x_nodes:]),
                 )
             )
             if self.__n_x_nodes > 0
@@ -71,11 +81,14 @@ def __call__(self, hypercube):
         """
         Prior for adaptive flex-knot.
 
-        hypercube = [N, y0, x1, y1, x2, y2, ..., x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
-        where Nmax is the greatest allowed value of floor(N).
+        hypercube = [N, y0, x1, y1, x2, y2, ...,
+                     x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
+        where Nmax is the greatest allowed value of floor(N), i.e.
+        the maximum number of nodes.
+
         """
         prior = np.empty(hypercube.shape)
-        prior[0] = self._N_prior(hypercube[0:1])
+        prior[[0]] = self._N_prior(hypercube[0:1])
         self.__n_x_nodes = int(prior[0])
         prior[1:] = super().__call__(hypercube[1:])
         return prior
diff --git a/flexknot/utils.py b/flexknot/utils.py
new file mode 100644
index 0000000..81a703d
--- /dev/null
+++ b/flexknot/utils.py
@@ -0,0 +1,143 @@
+"""Utilities for dealing with the interleaved arrangement of theta."""
+
+import numpy as np
+
+
+def validate_theta(theta, adaptive):
+    """Check theta for mistakes.
+
+    theta must contain an odd/even number of elements for an adaptive/
+    non-adaptive flex-knot,
+
+    In the adaptive case, floor(theta[0])-2 must not be greater than the
+    provided number of internal nodes.
+
+    Parameters
+    ----------
+    theta : array-like
+
+    adaptive : bool
+
+    """
+    if adaptive:
+        if len(theta) % 2 != 1:
+            raise ValueError(
+                "theta must contain an odd number of elements "
+                "for an adaptive flex-knot."
+            )
+        if np.floor(theta[0]) > (len(theta) + 1) / 2:
+            raise ValueError("n = ceil(theta[0]) exceeds the "
+                             "number of internal nodes.")
+
+    if not adaptive:
+        if len(theta) > 1 and len(theta) % 2 != 0:
+            raise ValueError("theta must contain an even number of elements "
+                             "for a non-adaptive flex-knot.")
+
+
+def get_theta_n(theta):
+    """Extract the first N parameters from theta.
+
+    Parameters
+    ----------
+    theta : array-like
+    [N, y0, x1, y1, x2, y2, ..., x_(Nmax-2), y_(Nmax-2), y_(Nmax-1)],
+    where Nmax is the maximum value of floor(N).
+
+    Returns
+    -------
+    theta_n : array-like
+    [y0, x1, y1, x2, y2, ..., x_floor(N)-2, y_floor(N)-2, y_(Nmax-1)]
+
+    """
+    validate_theta(theta, adaptive=True)
+
+    n = np.floor(theta[0]).astype(int) - 2
+    # w = -1 case, return empty
+    if -2 == n:
+        return np.array([])
+    theta_n = np.concatenate(
+        (
+            theta[1:2*n+2],  # y0 and internal x and y
+            theta[-1:],  # y end node
+        )
+    )
+    return theta_n
+
+
+def create_theta(x_nodes, y_nodes):
+    """Interleave x and y nodes to create theta.
+
+    Parameters
+    ----------
+    x_nodes : array-like
+    [x1, ... x_(N-2)]
+
+    y_nodes : array-like
+    [y0, y1, ..., y_(N-2), y_(N-1)]
+
+    Returns
+    -------
+    theta : array-like
+    [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)]
+
+    """
+    if len(y_nodes) == 1:
+        return y_nodes
+    if len(x_nodes) + 2 != len(y_nodes):
+        raise ValueError("y_nodes must have exactly two "
+                         "more elements than x_nodes")
+    n = len(x_nodes)
+    theta = np.zeros(len(x_nodes) + len(y_nodes))
+    theta[1:2*n+1:2] = x_nodes
+    theta[0:2*n+2:2] = y_nodes[:-1]
+    theta[-1] = y_nodes[-1]
+    return theta
+
+
+def get_x_nodes_from_theta(theta, adaptive):
+    """Get the x nodes from theta.
+
+    Parameters
+    ----------
+    theta : array-like
+    [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)]
+
+    adaptive : bool
+    Whether theta is for an adaptive or non-adaptive flex-knot.
+
+    Returns
+    -------
+    x_nodes : array-like
+    [x1, ... x_(N-2)]
+
+    """
+    validate_theta(theta, adaptive)
+    if adaptive:
+        theta = get_theta_n(theta)
+    n = len(theta) // 2 - 1
+    return theta[1:2*n+1:2]
+
+
+def get_y_nodes_from_theta(theta, adaptive):
+    """Get the y nodes from theta.
+
+    Parameters
+    ----------
+    theta : array-like
+    [y0, x1, y1, x2, y2, ..., x_(N-2), y_(N-2), y_(N-1)]
+
+    adaptive : bool
+    Whether theta is for an adaptive or non-adaptive flex-knot.
+
+    Returns
+    -------
+    y_nodes : array-like
+    [y0, y1, ..., y_(N-2), y_(N-1)]
+
+    """
+    validate_theta(theta, adaptive)
+    if adaptive:
+        theta = get_theta_n(theta)
+    n = len(theta) // 2 - 1
+    return np.concatenate((theta[0:2*n+2:2], theta[-1:]))
diff --git a/pyproject.toml b/pyproject.toml
index 4c0990a..fe9afd0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,5 +10,5 @@ dependencies = ["numpy", "scipy", "pypolychord"]
 repository = "https://github.com/Ormorod/flexknot"
 
 [project.optional-dependencies]
-dev = ["black"]
+dev = ["pytest", "flake8", "pydocstyle"]
  
diff --git a/tests/test_area.py b/tests/test_area.py
index b6c22bf..f898a64 100644
--- a/tests/test_area.py
+++ b/tests/test_area.py
@@ -1,6 +1,6 @@
 import numpy as np
 from flexknot import FlexKnot, AdaptiveKnot
-from flexknot.helper_functions import get_theta_n
+from flexknot.utils import get_theta_n
 
 
 def test_area():
diff --git a/tests/test_flexknot.py b/tests/test_flexknot.py
index a363d94..95986b0 100644
--- a/tests/test_flexknot.py
+++ b/tests/test_flexknot.py
@@ -18,8 +18,8 @@ def test_flexknot():
     x_nodes = np.sort(rng.uniform(x_min, x_max, n))
     y_nodes = rng.uniform(-10, 10, n + 2)
     theta = np.zeros(len(x_nodes) + len(y_nodes))
-    theta[1 : 2 * n + 1 : 2] = x_nodes
-    theta[0 : 2 * n + 2 : 2] = y_nodes[:-1]
+    theta[1:2*n+1:2] = x_nodes
+    theta[0:2*n+2:2] = y_nodes[:-1]
     theta[-1] = y_nodes[-1]
     xs = np.linspace(x_min, x_max, 100)
     assert np.all(
@@ -30,7 +30,7 @@ def test_flexknot():
 
 def test_adaptive_flexknot():
     """
-    Test that AdaptiveKnot returns the same results as FlexKnot with appropriate arguments.
+    Test that AdaptiveKnot returns the same results as FlexKnot.
     """
     x_min = 0
     x_max = 6
@@ -38,7 +38,8 @@ def test_adaptive_flexknot():
     theta_n = np.array([0, 1, 1, 2, 2, 3, 3, 6])
     xs = np.linspace(x_min, x_max, 100)
     assert np.all(
-        FlexKnot(x_min, x_max)(xs, theta_n) == AdaptiveKnot(x_min, x_max)(xs, theta)
+        FlexKnot(x_min, x_max)(xs, theta_n)
+        == AdaptiveKnot(x_min, x_max)(xs, theta)
     )
 
 
diff --git a/tests/test_helper_functions.py b/tests/test_helper_functions.py
index 5164364..5d7beec 100644
--- a/tests/test_helper_functions.py
+++ b/tests/test_helper_functions.py
@@ -1,9 +1,9 @@
 """
-Test each of the functions from helper_functions.py
+Test utils.py
 """
 
 import numpy as np
-from flexknot.helper_functions import (
+from flexknot.utils import (
     create_theta,
     get_theta_n,
     get_x_nodes_from_theta,
@@ -17,7 +17,8 @@
 
 def test_get_theta_n():
     """
-    Check that get_theta_n() extracts the correct elements from [3, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6].
+    Check that get_theta_n() extracts the correct elements from
+    [3, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6].
     """
 
     theta = np.array([5, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6])
diff --git a/tests/test_likelihood.py b/tests/test_likelihood.py
index f37b8b8..6a18260 100644
--- a/tests/test_likelihood.py
+++ b/tests/test_likelihood.py
@@ -4,17 +4,18 @@
 import numpy as np
 from scipy.special import erf
 from flexknot import FlexKnotLikelihood
-from flexknot.helper_functions import create_theta
+from flexknot.utils import create_theta
 
 
 def test_likelihood():
     """
     Test the likelihood function against a trivial case with just sigma_y.
 
-    There are two datapoints at (0, 0) and (1, 1), and the interpolation function
-    is a single straight line between the same two points. Sigma is set to 1.
+    There are two datapoints at (0, 0) and (1, 1), and the flex-knot
+    is a single straight line between the same two points. sigma is set to 1.
 
     The value for the likelihood in this case should be -ln(2π).
+
     """
     x_min, x_max = 0, 1
     x_nodes = np.array([])
@@ -25,19 +26,22 @@ def test_likelihood():
     y_data = np.array([0, 1])
     sigma = 1
 
-    l = FlexKnotLikelihood(x_min, x_max, x_data, y_data, sigma, adaptive=False)
-    assert l(theta)[0] == -np.log(2 * np.pi)
+    logl = FlexKnotLikelihood(x_min, x_max, x_data, y_data, sigma,
+                              adaptive=False)
+    assert logl(theta)[0] == -np.log(2 * np.pi)
 
 
 def test_likelihood_sigma_x():
     """
-    Test the likelihood function against a trivial case with sigma_x and sigma_y.
+    Test the likelihood function against a trivial case
+    with sigma_x and sigma_y.
 
-    There are two datapoints at (0, 0) and (1, 1), and the interpolation function
-    is a single straight line between the same two points. Sigma is set to 1 for
-    both x and y.
+    There are two datapoints at (0, 0) and (1, 1), and the flex-knot
+    is a single straight line between the same two points. Sigma is set to 1
+    for both x and y.
 
-    The value for the likelihood in this case should be ln[(1/16π)(erf(1)-erf(0))(erf(0)-erf(-1))].
+    The value for the likelihood in this case should be
+    ln[(1/16π)(erf(1)-erf(0))(erf(0)-erf(-1))].
     """
 
     x_min, x_max = 0, 1
@@ -49,7 +53,9 @@ def test_likelihood_sigma_x():
     y_data = np.array([0, 1])
     sigma = np.array([1, 1])
 
-    l = FlexKnotLikelihood(x_min, x_max, x_data, y_data, sigma, adaptive=False)
+    logl = FlexKnotLikelihood(x_min, x_max, x_data, y_data, sigma,
+                              adaptive=False)
     assert np.isclose(
-        l(theta)[0], np.log((erf(1) - erf(0)) * (erf(0) - erf(-1)) / (16 * np.pi))
+        logl(theta)[0], np.log((erf(1) - erf(0)) * (erf(0) - erf(-1))
+                               / (16 * np.pi))
     )
diff --git a/tests/test_prior.py b/tests/test_prior.py
index 0c29513..a3689b0 100644
--- a/tests/test_prior.py
+++ b/tests/test_prior.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 from flexknot import AdaptiveKnotPrior, FlexKnotPrior
-from flexknot.helper_functions import get_x_nodes_from_theta
+from flexknot.utils import get_x_nodes_from_theta
 
 rng = np.random.default_rng()
 x_min = 0
@@ -31,6 +31,8 @@ def test_adaptiveknotprior_x_nodes_are_sorted():
     """
 
     hypercube = rng.random(2 * N_max - 1)
-    prior = AdaptiveKnotPrior(x_min, x_max, y_min, y_max, N_min, N_max)(hypercube)
+    prior = AdaptiveKnotPrior(
+        x_min, x_max, y_min, y_max, N_min, N_max
+    )(hypercube)
 
     assert np.all(np.diff(get_x_nodes_from_theta(prior, adaptive=True)) >= 0)