diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 52a7426..30bf0d0 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -77,7 +77,5 @@ jobs:
         run: pip install -r doc_requirements.txt
       - name: Install cmdstanpy
         run: python -m cmdstanpy.install_cmdstan --version ${{ env.cmdstanVersion }}
-      - name: Test the example notebooks
-        run: pytest docs/test_examples.py -v
-      - name: Build the documentation
+      - name: Build the documentation, including example notebooks
         run: doit docs
diff --git a/.gitignore b/.gitignore
index ce0076b..d143762 100644
--- a/.gitignore
+++ b/.gitignore
@@ -141,3 +141,5 @@ workspace/
 *.o
 *.ipynb
 .cmdstanpy_cache
+*.pdf
+*.png
diff --git a/conf.py b/conf.py
index 65437c3..b68b701 100644
--- a/conf.py
+++ b/conf.py
@@ -20,7 +20,7 @@
 ]
 html_theme = "sphinx_rtd_theme"
 html_sidebars = {}
-exclude_patterns = ["docs/_build", "docs/jupyter_execute", ".pytest_cache", "playground"]
+exclude_patterns = ["docs/_build", "docs/jupyter_execute", ".pytest_cache", "playground", "figures"]
 
 # Configure autodoc to avoid excessively long fully-qualified names.
 add_module_names = False
diff --git a/dev_requirements.txt b/dev_requirements.txt
index 9dffe9e..6de130f 100644
--- a/dev_requirements.txt
+++ b/dev_requirements.txt
@@ -133,7 +133,7 @@ cloudpickle==2.2.0
     #   -r gptools-torch/test_requirements.txt
     #   -r gptools-util/test_requirements.txt
     #   doit
-cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@develop
+cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@1612c04
     # via
     #   -r doc_requirements.txt
     #   -r gptools-stan/test_requirements.txt
diff --git a/docs/test_examples.py b/docs/test_examples.py
index 2a88296..ba45fb8 100644
--- a/docs/test_examples.py
+++ b/docs/test_examples.py
@@ -6,6 +6,7 @@
 import pathlib
 import pytest
 from typing import Any
+from unittest import mock
 
 
 def run_myst_notebook(path: str) -> Any:
@@ -19,8 +20,9 @@ def run_myst_notebook(path: str) -> Any:
         content = fp.read()
     reader: NbReader = create_nb_reader(path, md_config, nb_config, content)
     notebook = reader.read(content)
-    preprocessor = ExecutePreprocessor(timeout=timeout)
-    return preprocessor.preprocess(notebook, {"metadata": {"path": pathlib.Path(path).parent}})
+    with mock.patch.dict(os.environ, CI="true"):
+        preprocessor = ExecutePreprocessor(timeout=timeout)
+        return preprocessor.preprocess(notebook, {"metadata": {"path": pathlib.Path(path).parent}})
 
 
 notebooks = []
diff --git a/dodo.py b/dodo.py
index 8dfcdab..b2770e7 100644
--- a/dodo.py
+++ b/dodo.py
@@ -57,8 +57,8 @@
 
 # Build documentation at the root level (we don't have namespace-package-level documentation).
 with di.defaults(basename="docs"):
-    manager(name="html", actions=["sphinx-build . docs/_build"])
-    manager(name="tests", actions=["sphinx-build -b doctest . docs/_build"])
+    manager(name="html", actions=["sphinx-build -n . docs/_build"])
+    manager(name="tests", actions=["sphinx-build -b doctest . docs/_build", "pytest docs"])
 
 # Compile example notebooks to create html reports.
 for path in pathlib.Path.cwd().glob("gptools-*/**/*.ipynb"):
diff --git a/figures/.gitignore b/figures/.gitignore
deleted file mode 100644
index a136337..0000000
--- a/figures/.gitignore
+++ /dev/null
@@ -1 +0,0 @@
-*.pdf
diff --git a/figures/decision_tree.md b/figures/decision_tree.md
new file mode 100644
index 0000000..6b88d5b
--- /dev/null
+++ b/figures/decision_tree.md
@@ -0,0 +1,95 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+```{code-cell} ipython3
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+import numpy as np
+
+
+mpl.style.use("../jss.mplstyle")
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots()
+ax.set_aspect(1)
+
+offset = 2
+xy = {
+    (0, 0): ("how strong are\nthe data?", {"question": True}),
+    (-1, -1): ("non-centered\nparameterization", {"rotate": True}),
+    (0, -1): ("centered\nparameterization", {"rotate": True}),
+    (1, -1): ("other method?", {"rotate": True}),
+
+    (offset + 0.5, 0): ("is the GP on a\nregular grid?", {"question": True}),
+    (offset, -1): ("Fourier\nmethods", {"rotate": True}),
+    (offset + 1, -1): ("are the data\n\"dense\"?", {"question": True}),
+    (offset + 0.5, -2): ("Fourier inducing\npoints", {"rotate": True}),
+    (offset + 1.5, -2): ("structured\ndependencies", {"rotate": True}),
+}
+ax.scatter(*np.transpose(list(xy)))
+for (x, y), (text, kwargs) in xy.items():
+    kwargs = kwargs or {}
+    rotate = kwargs.pop("rotate", False)
+    question = kwargs.pop("question", False)
+    kwargs = {
+        "fontsize": "small",
+        "ha": "center",
+        "va": "top" if rotate else "center",
+        "rotation": 90 if rotate else 0,
+        "bbox": {
+            "edgecolor": "none" if question else "k",
+            "boxstyle": "round,pad=.5",
+            "facecolor": "k" if question else "w",
+        },
+        "color": "w" if question else "k",
+    } | kwargs
+    ax.text(x, y, text, **kwargs)
+
+
+def connect(xy1, xy2, frac=0.25, color="k", label=None):
+    x1, y1 = xy1
+    x2, y2 = xy2
+    xy = [xy1, (x1, y1 - frac), (x2, y1 - frac), xy2]
+    ax.plot(*np.transpose(xy), color=color, zorder=0, linewidth=1)
+    if label:
+        kwargs = {
+            "fontsize": "small",
+            "rotation": 90,
+            "ha": "center",
+            "va": "center",
+            "bbox": {
+                "edgecolor": "none",
+                "facecolor": "w",
+            }
+        }
+        ax.text(x2, (y1 - frac + y2) / 2, label, **kwargs)
+
+connect((0, 0), (-1, -1), label="weak")
+connect((0, 0), (0, -1), label="strong")
+connect((0, 0), (1, -1), label="very\nstrong")
+
+connect((offset + 0.5, 0), (offset, -1), label="yes")
+connect((offset + 0.5, 0), (offset + 1, -1), label="no")
+connect((offset + 1, -1), (offset + 0.5, -2), label="yes")
+connect((offset + 1, -1), (offset + 1.5, -2), label="no")
+
+ax.set_xlim(-1.25, 3.75)
+ax.set_ylim(-3, 0.25)
+ax.set_axis_off()
+ax.text(-1, 0, "(a)", ha="center", va="center")
+ax.text(1.8, 0, "(b)", ha="center", va="center")
+fig.tight_layout()
+fig.savefig("decision_tree.pdf", bbox_inches="tight")
+fig.savefig("decision_tree.png", bbox_inches="tight")
+```
diff --git a/figures/kernels.md b/figures/kernels.md
new file mode 100644
index 0000000..df1a8cd
--- /dev/null
+++ b/figures/kernels.md
@@ -0,0 +1,87 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+```{code-cell} ipython3
+from gptools.util.kernels import ExpQuadKernel, MaternKernel
+from gptools.util.fft import transform_irfft, evaluate_rfft_scale
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+import numpy as np
+
+
+mpl.style.use("../jss.mplstyle")
+```
+
+```{code-cell} ipython3
+np.random.seed(9)  # Seed picked for good legend positioning. Works for any though.
+fig, axes = plt.subplots(2, 2)
+length_scale = 0.2
+kernels = {
+    "squared exp.": lambda period: ExpQuadKernel(1, length_scale, period),
+    "Matern ³⁄₂": lambda period: MaternKernel(1.5, 1, length_scale, period),
+}
+
+x = np.linspace(0, 1, 101, endpoint=False)
+z = np.random.normal(0, 1, x.size)
+
+for ax, (key, kernel) in zip(axes[1], kernels.items()):
+    value = kernel(None).evaluate(0, x[:, None])
+    line, = axes[0, 0].plot(x, value, ls="--")
+    rfft = kernel(1).evaluate_rfft([x.size])
+    value = np.fft.irfft(rfft, x.size)
+    axes[0, 1].plot(rfft, label=key)
+    axes[0, 0].plot(x, value, color=line.get_color())
+    
+    for maxf, ls in [(x.size // 2 + 1, "-"), (5, "--"), (3, ":")]:
+        rfft_scale = evaluate_rfft_scale(rfft=rfft, size=x.size)
+        rfft_scale[maxf:] = 0
+        f = transform_irfft(z, np.zeros_like(z), rfft_scale=rfft_scale)
+        ax.plot(x, f, ls=ls, color=line.get_color(), label=fr"$\xi_\max={maxf}$")
+        
+    ax.set_xlabel("position $x$")
+    ax.set_ylabel(f"{key} GP $f$")
+    
+ax = axes[0, 0]
+ax.set_ylabel("kernel $k(0,x)$")
+ax.set_xlabel("position $x$")
+ax.legend([
+    mpl.lines.Line2D([], [], ls="--", color="gray"),
+    mpl.lines.Line2D([], [], ls="-", color="gray"),
+], ["standard", "periodic"], fontsize="small")
+ax.text(0.05, 0.05, "(a)", transform=ax.transAxes)
+ax.yaxis.set_ticks([0, 0.5, 1])
+    
+ax = axes[0, 1]
+ax.set_yscale("log")
+ax.set_ylim(1e-5, x.size)
+ax.set_xlabel(r"frequency $\xi$")
+ax.set_ylabel(r"Fourier kernel $\tilde k=\phi\left(k\right)$")
+ax.legend(fontsize="small", loc="center right")
+ax.text(0.95, 0.95, "(b)", transform=ax.transAxes, ha="right", va="top")
+
+ax = axes[1, 0]
+ax.legend(fontsize="small", loc="lower center")
+ax.text(0.95, 0.95, "(c)", transform=ax.transAxes, ha="right", va="top")
+
+ax = axes[1, 1]
+ax.legend(fontsize="small", loc="lower center")
+ax.sharey(axes[1, 0])
+ax.text(0.95, 0.95, "(d)", transform=ax.transAxes, ha="right", va="top")
+
+for ax in [axes[0, 0], *axes[1]]:
+    ax.xaxis.set_ticks([0, 0.5, 1])
+
+fig.tight_layout()
+fig.savefig("kernels.pdf", bbox_inches="tight")
+fig.savefig("kernels.png", bbox_inches="tight")
+```
diff --git a/figures/profile.md b/figures/profile.md
index 13195c1..f3ebd9e 100644
--- a/figures/profile.md
+++ b/figures/profile.md
@@ -11,6 +11,8 @@ kernelspec:
   name: python3
 ---
 
+# Profiling of different methods and parameterizations
+
 ```{code-cell} ipython3
 import cmdstanpy
 import itertools as it
@@ -125,16 +127,15 @@ ls_by_parameterization = {
     "non_centered": "--",
 }
 
-fig = plt.figure(layout="constrained")
+fig = plt.figure()
 # First ratio is manually fiddled to match the heights of ax1 and ax3
 gs = fig.add_gridspec(2, 2)
 ax1 = fig.add_subplot(gs[0, 0])
 ax2 = fig.add_subplot(gs[0, 1], sharex=ax1, sharey=ax1)
 ax3 = fig.add_subplot(gs[1, 0], sharey=ax1)
 
-gs4 = mpl.gridspec.GridSpecFromSubplotSpec(2, 1, gs[1, 1], hspace=0,
-                                           wspace=0,
-                                           height_ratios=[2, 1])
+gs4 = mpl.gridspec.GridSpecFromSubplotSpec(2, 1, gs[1, 1], height_ratios=[2, 1],
+                                           hspace=0.1)
 ax4t = fig.add_subplot(gs4[0])
 ax4b = fig.add_subplot(gs4[1])
 
@@ -146,6 +147,7 @@ ax.set_xlabel(r"size $n$")
 ax.set_ylabel(r"runtime (seconds)")
 
 for i, ax in [(0, ax1), (-1, ax2)]:
+    print(f"kappa = {LOG10_NOISE_SCALES[i]}")
     for key, value in durations.items():
         method, parameterization = key.split("_", 1)
         line, = ax.plot(
@@ -154,13 +156,19 @@ for i, ax in [(0, ax1), (-1, ax2)]:
             ls=ls_by_parameterization[parameterization],
         )
         line.set_markeredgecolor("w")
+        
+        f = np.isfinite(value[i])
+        fit = np.polynomial.Polynomial.fit(np.log(SIZES[f])[-3:], np.log(value[i][f])[-3:], 1)
+        fit = fit.convert()
+        print(f"{key}: n ** {fit.coef[1]:.3f}")
+        
     ax.set_xscale("log")
     ax.set_yscale("log")
     ax.text(0.05, 0.95, fr"({'ab'[i]}) $\kappa=10^{{{LOG10_NOISE_SCALES[i]:.0f}}}$",
             transform=ax.transAxes, va="top")
     ax.set_xlabel("size $n$")
     ax.set_ylabel("duration (seconds)")
-
+    print()
 
 ax = ax3
 ax.set_xlabel(r"noise scale $\kappa$")
@@ -177,7 +185,8 @@ for parameterization in ["non_centered", "centered"]:
         ax.plot(NOISE_SCALES, y, color=mappable.to_rgba(size),
                 ls=ls_by_parameterization[parameterization])
 ax.text(0.05, 0.95, "(c)", transform=ax.transAxes, va="top")
-fig.colorbar(mappable, ax=ax, label="size $n$")
+ax.set_ylim(top=150)
+# We'll add the colorbar later below after laying out the figure.
 
 ax = ax4b
 ax.set_xlabel(r"noise scale $\kappa$")
@@ -204,23 +213,21 @@ for ax in [ax4b, ax4t]:
     ax.set_xscale("log")
     ax.ticklabel_format(axis="y", scilimits=(0, 0), useMathText=True)
 ax4b.legend(loc="lower right")
-ax4b.set_ylim(-10e3, -7e3)
+ax4b.set_ylim(-10e3, -6.5e3)
 ax4t.set_ylim(-1300, 600)
 ax4b.yaxis.get_offset_text().set_visible(False)
-ax4b.set_yticks([-10e3, -8e3])
+ax4b.set_yticks([-9e3, -7e3])
 
 # Visibility adjustment must happen after plotting.
-ax4t.set_ylabel(r"log p.p.d. difference $\Delta$", y=0.2)
+ax4t.set_ylabel(r"log p.d. difference $\Delta$", y=0.2)
 ax4t.set_xticklabels([])
 ax4t.spines["bottom"].set_visible(False)
 plt.setp(ax.xaxis.get_majorticklines(), visible=False)
 plt.setp(ax.xaxis.get_minorticklines(), visible=False)
 ax4t.text(0.05, 0.95, "(d)", transform=ax.transAxes, va="top")
 
-# Disable automatic layout.
-fig.get_layout_engine().set(rect=[0, 0, 1, 0.93])
-fig.draw_without_rendering()
-fig.set_layout_engine(None)
+# Lay out the figure before adding extra elements.
+gs.tight_layout(fig, rect=[0, 0, 1, 0.93], h_pad=0)
 
 # Add the broken axis markers.
 angle = np.deg2rad(30)
@@ -256,7 +263,7 @@ handles_labels = [
 for method, color in color_by_method.items():
     handles_labels.append((
         mpl.lines.Line2D([], [], color=color),
-        method,
+        "Fourier" if method == "fourier" else method,
     ))
 bbox1 = ax1.get_position()
 bbox2 = ax2.get_position()
@@ -264,5 +271,27 @@ bbox_to_anchor = [bbox1.xmin, 0, bbox2.xmax - bbox1.xmin, 1]
 legend = fig.legend(*zip(*handles_labels), fontsize="small", loc="upper center",
                     ncol=5, bbox_to_anchor=bbox_to_anchor)
 
+# Finally add the colorbar.
+rect = ax3.get_position() 
+rect = (
+    rect.xmin + 0.525 * rect.width, 
+    rect.ymin + rect.height * 0.89, 
+    rect.width * 0.45,
+    rect.height * 0.06,
+)
+cax = fig.add_axes(rect)  # rect=(left, bottom, width, height)
+cb = fig.colorbar(mappable, cax=cax, orientation="horizontal")
+cax.tick_params(labelsize="small")
+cax.set_ylabel("size $n$", fontsize="small", rotation=0, ha="right", va="center")
+
 fig.savefig("scaling.pdf", bbox_inches="tight")
+fig.savefig("scaling.png", bbox_inches="tight")
+```
+
+```{code-cell} ipython3
+# Show the typical runtimes for the "compromise" noise scale for the standard approach.
+np.transpose([
+    SIZES,
+    durations["standard_centered"][LOG10_NOISE_SCALES.size // 2],
+])
 ```
diff --git a/gptools-stan/README.rst b/gptools-stan/README.rst
index 3651fe9..a67cbc2 100644
--- a/gptools-stan/README.rst
+++ b/gptools-stan/README.rst
@@ -6,6 +6,7 @@
 
     docs/poisson_regression/poisson_regression
     docs/trees/trees
+    docs/tube/tube
 
 The interface definitions below provide a comprehensive overview of the functionality offered by the Stan library. Please see the example :doc:`docs/poisson_regression/poisson_regression` for an illustration of how to use the library.
 
diff --git a/gptools-stan/docs/poisson_regression/poisson_regression.md b/gptools-stan/docs/poisson_regression/poisson_regression.md
index d654e1e..edd9acf 100644
--- a/gptools-stan/docs/poisson_regression/poisson_regression.md
+++ b/gptools-stan/docs/poisson_regression/poisson_regression.md
@@ -218,7 +218,7 @@ The non-centered parameterization of the graph Gaussian process is even faster,
 
 ## Gaussian process using fast Fourier transforms
 
-If, as in this example, observations are regularly spaced, we can evaluate the likelihood using the fast Fourier transform. Similarly, we can construct a non-centered parameterization by transforming Fourier coefficients with {stan:func}`gp_transform_irfft`.
+If, as in this example, observations are regularly spaced, we can evaluate the likelihood using the fast Fourier transform. Similarly, we can construct a non-centered parameterization by transforming Fourier coefficients with {stan:func}`gp_transform_inv_rfft`.
 
 ```{literalinclude} poisson_regression_fourier_centered.stan
    :language: stan
diff --git a/gptools-stan/docs/poisson_regression/poisson_regression_fourier_centered.stan b/gptools-stan/docs/poisson_regression/poisson_regression_fourier_centered.stan
index 6b6892f..908cb55 100644
--- a/gptools-stan/docs/poisson_regression/poisson_regression_fourier_centered.stan
+++ b/gptools-stan/docs/poisson_regression/poisson_regression_fourier_centered.stan
@@ -22,6 +22,6 @@ transformed parameters {
 
 model {
     // Fourier Gaussian process and observation model.
-    eta ~ gp_rfft(zeros_vector(n), gp_evaluate_rfft_scale(cov_rfft, n));
+    eta ~ gp_rfft(zeros_vector(n), cov_rfft);
     y ~ poisson_log(eta);
 }
diff --git a/gptools-stan/docs/poisson_regression/poisson_regression_fourier_non_centered.stan b/gptools-stan/docs/poisson_regression/poisson_regression_fourier_non_centered.stan
index 6e73cc1..0c79642 100644
--- a/gptools-stan/docs/poisson_regression/poisson_regression_fourier_non_centered.stan
+++ b/gptools-stan/docs/poisson_regression/poisson_regression_fourier_non_centered.stan
@@ -18,7 +18,7 @@ transformed parameters {
     vector[n] eta;
     vector[n %/% 2 + 1] cov_rfft = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n, 10)
         + epsilon;
-    eta = gp_transform_irfft(z, zeros_vector(n), gp_evaluate_rfft_scale(cov_rfft, n));
+    eta = gp_transform_inv_rfft(z, zeros_vector(n), cov_rfft);
 }
 
 model {
diff --git a/gptools-stan/docs/poisson_regression/poisson_regression_graph_centered.stan b/gptools-stan/docs/poisson_regression/poisson_regression_graph_centered.stan
index fcb1689..9f6a02a 100644
--- a/gptools-stan/docs/poisson_regression/poisson_regression_graph_centered.stan
+++ b/gptools-stan/docs/poisson_regression/poisson_regression_graph_centered.stan
@@ -12,7 +12,7 @@ data {
 }
 
 transformed data {
-    array [n] int degrees = in_degrees(n, edge_index);
+    array [n] int degrees = out_degrees(n, edge_index);
 }
 
 parameters {
diff --git a/gptools-stan/docs/poisson_regression/poisson_regression_graph_non_centered.stan b/gptools-stan/docs/poisson_regression/poisson_regression_graph_non_centered.stan
index 56d2315..e595f36 100644
--- a/gptools-stan/docs/poisson_regression/poisson_regression_graph_non_centered.stan
+++ b/gptools-stan/docs/poisson_regression/poisson_regression_graph_non_centered.stan
@@ -12,7 +12,7 @@ data {
 }
 
 transformed data {
-    array [n] int degrees = in_degrees(n, edge_index);
+    array [n] int degrees = out_degrees(n, edge_index);
 }
 
 parameters {
@@ -21,7 +21,7 @@ parameters {
 
 transformed parameters {
     // Transform white noise to a sample from the graph Gaussian process.
-    vector[n] eta = gp_graph_exp_quad_cov_transform(
+    vector[n] eta = gp_transform_inv_graph_exp_quad_cov(
         z, zeros_vector(n), X, sigma, length_scale, edge_index, degrees, epsilon);
 }
 
diff --git a/gptools-stan/docs/trees/trees.md b/gptools-stan/docs/trees/trees.md
index 15aad64..7627589 100644
--- a/gptools-stan/docs/trees/trees.md
+++ b/gptools-stan/docs/trees/trees.md
@@ -11,6 +11,8 @@ kernelspec:
   name: python3
 ---
 
+# Density of T. panamensis on a 50 ha plot in Panama
+
 ```{code-cell} ipython3
 from gptools import util
 from gptools.stan import compile_model
@@ -39,6 +41,7 @@ ax.set_xlabel("easting (km)")
 ax.set_ylabel("northing (km)")
 fig.colorbar(im, ax=ax, location="top", label="tree density")
 fig.tight_layout()
+frequency.shape
 ```
 
 ```{code-cell} ipython3
@@ -68,8 +71,8 @@ padded_rows, padded_cols
 ```{code-cell} ipython3
 # Compile and fit the model.
 model = compile_model(stan_file="trees.stan")
-niter = 3 if "CI" in os.environ else 500
-fit = model.sample(data, chains=1, iter_warmup=niter, iter_sampling=2 * niter, seed=seed)
+niter = 3 if "CI" in os.environ else 1000
+fit = model.sample(data, chains=1, iter_warmup=niter, iter_sampling=niter, seed=seed)
 print(fit.diagnose())
 ```
 
@@ -89,8 +92,8 @@ def evaluate_error(actual, prediction, error, num_bs=1000):
         return (np.square(bs_actual - bs_prediction) / np.maximum(bs_actual, 1)).mean(axis=-1)
     else:
         raise ValueError(error)
-        
-        
+
+
 def filter_estimate(frequency, train_mask, scale):
     smoothed_mask = ndimage.gaussian_filter(train_mask.astype(float), scale)
     smoothed_masked_frequency = ndimage.gaussian_filter(np.where(train_mask, frequency, 0), scale)
@@ -167,7 +170,7 @@ cax.xaxis.set_ticks_position("top")
 cax.xaxis.set_label_position("top")
 
 ax3 = fig.add_subplot(gs2[0])
-ax3.scatter(fit.length_scale * delta * 1e3, fit.sigma, marker=".", 
+ax3.scatter(fit.length_scale * delta * 1e3, fit.sigma, marker=".",
             alpha=0.25)
 ax3.set_xlabel(r"correlation length $\ell$ (m)")
 ax3.set_ylabel(r"marginal scale $\sigma$")
@@ -190,4 +193,5 @@ ax3.text(0.05, 0.95, "(b)", va="top", transform=ax3.transAxes)
 ax4.text(0.05, 0.95, "(d)", va="top", transform=ax4.transAxes)
 
 fig.savefig("trees.pdf", bbox_inches="tight")
+fig.savefig("trees.png", bbox_inches="tight")
 ```
diff --git a/gptools-stan/docs/trees/trees.stan b/gptools-stan/docs/trees/trees.stan
index cb463dd..d795db1 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -20,40 +20,34 @@ parameters {
     // Mean log rate for the trees.
     real mu;
     // Kernel parameters and averdispersion parameter for the negative binomial distribution.
-    real<lower=0> sigma, length_scale, phi;
+    real<lower=0> sigma, kappa;
+    real<lower=log(2), upper=log(28)> log_length_scale;
 }
 
 transformed parameters {
+    real<lower=0> length_scale = exp(log_length_scale);
     // Evaluate the RFFT of the Matern 3/2 kernel on the padded grid.
     matrix[num_rows_padded, num_cols_padded %/% 2 + 1] rfft2_cov = epsilon +
         gp_periodic_matern_cov_rfft2(1.5, num_rows_padded, num_cols_padded, sigma,
-        ones_vector(2) * length_scale, [num_rows_padded, num_cols_padded]');
+        length_scale, [num_rows_padded, num_cols_padded]');
     // Transform from white-noise to a Gaussian process realization.
-    matrix[num_rows_padded, num_cols_padded] f = gp_transform_irfft2(z, mu + zeros_matrix(num_rows_padded, num_cols_padded),
-                                             gp_evaluate_rfft2_scale(rfft2_cov, num_cols_padded));
+    matrix[num_rows_padded, num_cols_padded] f = gp_transform_inv_rfft2(
+        z, rep_matrix(mu, num_rows_padded, num_cols_padded), rfft2_cov);
 }
 
 model {
+    // Implies that eta ~ gp_rfft(mu, rfft2_cov) is a realization of the Gaussian process.
+    z ~ std_normal();
+    // Weakish priors on all other parameters.
+    mu ~ student_t(2, 0, 1);
+    sigma ~ student_t(2, 0, 1);
+    kappa ~ student_t(2, 0, 1);
     // Observation model with masking for negative values.
     for (i in 1:num_rows) {
         for (j in 1:num_cols) {
             if (frequency[i, j] >= 0) {
-                frequency[i, j] ~ neg_binomial_2(exp(f[i, j]), 1 / phi);
+                frequency[i, j] ~ neg_binomial_2(exp(f[i, j]), 1 / kappa);
             }
         }
     }
-    // Implies that eta ~ gp_rfft(mu, rfft2_cov) is a realization of the Gaussian process.
-    z ~ std_normal();
-    // We bound the correlation length below because the likelihood is flat and place an upper bound
-    // of 200m correlation which is broad enough to avoid the posterior samples hitting the
-    // boundary.
-    length_scale ~ uniform(1, 10);
-    // Prior on the overall density based on Gelman's default prior paper (although we here use
-    // count rather than binary regression).
-    mu ~ cauchy(0, 2.5);
-    // We want to allow a reasonable amount of variation but restrict excessively variable
-    // functions. So we use slightly lighter tail than the Cauchy to constrain the sampler.
-    sigma ~ student_t(2, 0, 1);
-    // The overdispersion parameter again has a weak prior.
-    phi ~ cauchy(0, 1);
 }
diff --git a/gptools-stan/docs/tube/tube.md b/gptools-stan/docs/tube/tube.md
index f78d479..af58fd5 100644
--- a/gptools-stan/docs/tube/tube.md
+++ b/gptools-stan/docs/tube/tube.md
@@ -11,6 +11,8 @@ kernelspec:
   name: python3
 ---
 
+# Passengers on the London Underground network
+
 ```{code-cell} ipython3
 from gptools.stan import compile_model
 from gptools.util import encode_one_hot
@@ -21,6 +23,7 @@ from matplotlib import pyplot as plt
 import networkx as nx
 import numpy as np
 import os
+from scipy.spatial.distance import pdist
 from scipy import stats
 
 mpl.style.use("../../../jss.mplstyle")
@@ -74,6 +77,7 @@ one_hot_degrees = encode_one_hot(degrees.clip(max=max_degree) - 1)
 np.random.seed(seed)
 train_mask = np.random.uniform(0, 1, y.size) < train_frac
 
+distances = pdist(X)
 data = {
     "num_stations": graph.number_of_nodes(),
     "num_edges": edge_index.shape[1],
@@ -88,12 +92,13 @@ data = {
     "include_zone_effect": 1,
     "include_degree_effect": 1,
 }
+distances.min(), distances.max()
 ```
 
 ```{code-cell} ipython3
 niter = 3 if "CI" in os.environ else 1000
 model_with_gp = compile_model(stan_file="tube.stan")
-fit_with_gp = model_with_gp.sample(data, chains=1, iter_warmup=niter, iter_sampling=niter, 
+fit_with_gp = model_with_gp.sample(data, chains=1, iter_warmup=niter, iter_sampling=niter,
                                    seed=seed, adapt_delta=0.9)
 print(fit_with_gp.diagnose())
 ```
@@ -199,6 +204,7 @@ ax4.text(0.95, 0.95, "(d)", transform=ax4.transAxes, va="top", ha="right")
 
 fig.tight_layout()
 fig.savefig("tube.pdf", bbox_inches="tight")
+fig.savefig("tube.png", bbox_inches="tight")
 ```
 
 On the one hand, the three northern stations of the [Hainault Loop](https://en.wikipedia.org/wiki/Hainault_Loop) ([Roding Valley](https://en.wikipedia.org/wiki/Roding_Valley_tube_station), [Chigwell](https://en.wikipedia.org/wiki/Chigwell_tube_station), and [Grange Hill](https://en.wikipedia.org/wiki/Grange_Hill_tube_station)) are underused because they are serviced by only three trains an hour whereas nearby stations (such as [Hainault](https://en.wikipedia.org/wiki/Hainault_tube_station), [Woodford](https://en.wikipedia.org/wiki/Woodford_tube_station), and [Buckhurst Hill](https://en.wikipedia.org/wiki/Buckhurst_Hill_tube_station)) are serviced by twelve trains an hour. On the other hand, [Canary Wharf](https://en.wikipedia.org/wiki/Canary_Wharf_tube_station) at the heart of the financial district has much higher use than would be expected for a station that only serves a single line in zone 2.
@@ -208,7 +214,7 @@ Let's compare the predictive performance on the held out data with and without t
 ```{code-cell} ipython3
 model_without_gp = compile_model(stan_file="tube_without_gp.stan")
 niter = 3 if "CI" in os.environ else 1000
-fit_without_gp = model_without_gp.sample(data, chains=1, iter_warmup=niter, 
+fit_without_gp = model_without_gp.sample(data, chains=1, iter_warmup=niter,
                                          iter_sampling=niter, seed=seed, adapt_delta=0.9)
 print(fit_without_gp.diagnose())
 ```
@@ -225,13 +231,13 @@ bs_samples = []
 for label, fit in [("with GP", fit_with_gp), ("without GP", fit_without_gp)]:
     log_mean = fit.stan_variable("log_mean")[..., ~train_mask]
     kappa = fit.stan_variable("kappa")
-    
+
     # Evaluate the score and bootstrapped error.
     log_score = stats.norm(log_mean, kappa[:, None]).logpdf(test).mean(axis=0)
     x = np.random.dirichlet(np.ones(log_score.shape[0]), 1000) @ log_score
     print(f"{label}: {log_score.mean():.3f} +- {np.std(x):.3f}")
     bs_samples.append(x)
-    
+
     ax.errorbar(test, log_mean.mean(axis=0), log_mean.std(axis=0), ls="none", marker=".", label=label)
 
 ax.legend(fontsize="small")
@@ -239,5 +245,5 @@ ax.set_aspect("equal")
 fig.tight_layout()
 
 delta = np.diff(bs_samples, axis=0).squeeze()
-delta.mean(), delta.std()
+print(f"difference: {delta.mean():.3f} +- {delta.std():.3f}")
 ```
diff --git a/gptools-stan/docs/tube/tube.stan b/gptools-stan/docs/tube/tube.stan
index 98593f5..21980cb 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -11,25 +11,26 @@ data {
     matrix[num_stations, num_degrees] one_hot_degrees;
     real<lower=0> epsilon;
 
-    int include_zone_effect;
-    int include_degree_effect;
+    int include_zone_effect, include_degree_effect;
 }
 
 transformed data {
-    array [num_stations] int degrees = in_degrees(num_stations, edge_index);
+    array [num_stations] int degrees = out_degrees(num_stations, edge_index);
 }
 
 parameters {
     vector[num_stations] z;
     real mu;
-    real<lower=0> sigma, length_scale, kappa;
+    real<lower=0> sigma, kappa;
+    real<lower=log(0.32), upper=log(31)> log_length_scale;
     vector[num_zones] zone_effect;
     vector[num_degrees] degree_effect;
 }
 
 transformed parameters {
-    vector[num_stations] f = gp_graph_exp_quad_cov_transform(
-        z, zeros_vector(num_stations), station_locations, sigma, length_scale, edge_index, degrees, epsilon);
+    real length_scale = exp(log_length_scale);
+    vector[num_stations] f = gp_transform_inv_graph_exp_quad_cov(
+        z, zeros_vector(num_stations), station_locations, sigma, length_scale, edge_index);
     vector[num_stations] log_mean = mu + f
         + include_zone_effect * one_hot_zones * zone_effect
         + include_degree_effect * one_hot_degrees * degree_effect;
@@ -39,7 +40,6 @@ transformed parameters {
 model {
     z ~ std_normal();
     sigma ~ student_t(2, 0, 1);
-    length_scale ~ inv_gamma(2.0, 2.5);
     zone_effect ~ student_t(2, 0, 1);
     degree_effect ~ student_t(2, 0, 1);
     kappa ~ student_t(2, 0, 1);
diff --git a/gptools-stan/docs/tube/tube_without_gp.stan b/gptools-stan/docs/tube/tube_without_gp.stan
index 2aa9ef9..219f215 100644
--- a/gptools-stan/docs/tube/tube_without_gp.stan
+++ b/gptools-stan/docs/tube/tube_without_gp.stan
@@ -23,9 +23,9 @@ transformed parameters {
 }
 
 model {
-    zone_effect ~ cauchy(0, 1);
-    degree_effect ~ cauchy(0, 1);
-    kappa ~ cauchy(0, 1);
+    zone_effect ~ student_t(2, 0, 1);
+    degree_effect ~ student_t(2, 0, 1);
+    kappa ~ student_t(2, 0, 1);
     for (i in 1:num_stations) {
         if (passengers[i] >= 0) {
             log(passengers[i]) ~ normal(log_mean[i], kappa);
diff --git a/gptools-stan/gptools/stan/__init__.py b/gptools-stan/gptools/stan/__init__.py
index 5efa8db..4357b24 100644
--- a/gptools-stan/gptools/stan/__init__.py
+++ b/gptools-stan/gptools/stan/__init__.py
@@ -18,12 +18,14 @@ def compile_model(
         exe_file: OptionalPath = None, compile: Union[bool, str] = True,
         stanc_options: Optional[dict[str, Any]] = None,
         cpp_options: Optional[dict[str, Any]] = None,
-        user_header: OptionalPath = None) -> cmdstanpy.CmdStanModel:
+        user_header: OptionalPath = None, **kwargs) -> cmdstanpy.CmdStanModel:
     # Add gptools includes by default.
     stanc_options = stanc_options or {}
     stanc_options.setdefault("include-paths", []).append(get_include())
-    return cmdstanpy.CmdStanModel(model_name, stan_file, exe_file, compile, stanc_options,
-                                  cpp_options, user_header)
+    return cmdstanpy.CmdStanModel(
+        model_name=model_name, stan_file=stan_file, exe_file=exe_file, compile=compile,
+        stanc_options=stanc_options, cpp_options=cpp_options, user_header=user_header, **kwargs,
+    )
 
 
 if __name__ == "__main__":
diff --git a/gptools-stan/gptools/stan/gptools_fft1.stan b/gptools-stan/gptools/stan/gptools_fft1.stan
index c04868a..19f0837 100644
--- a/gptools-stan/gptools/stan/gptools_fft1.stan
+++ b/gptools-stan/gptools/stan/gptools_fft1.stan
@@ -2,9 +2,9 @@
 // I.e., x[1:3] includes x[1] to x[3]. More generally, x[i:j] comprises j - i + 1 elements. It could
 // at least have been exclusive on the right...
 
-vector gp_evaluate_rfft_scale(vector rfft_, int n) {
+vector gp_evaluate_rfft_scale(vector cov_rfft, int n) {
     int nrfft = n %/% 2 + 1;
-    vector[nrfft] result = n * rfft_ / 2;
+    vector[nrfft] result = n * cov_rfft / 2;
     // Check positive-definiteness.
     real minval = min(result);
     if (minval < 0) {
@@ -20,18 +20,6 @@ vector gp_evaluate_rfft_scale(vector rfft_, int n) {
     return sqrt(result);
 }
 
-/**
-Evaluate the scale of Fourier coefficients.
-
-The Fourier coefficients of a zero-mean Gaussian process with even covariance function are
-uncorrelated Gaussian random variables with zero mean. This function evaluates their scale.
-
-:param cov: Covariance between the origin and the rest of the domain.
-*/
-vector gp_evaluate_rfft_scale(vector cov) {
-    return gp_evaluate_rfft_scale(get_real(rfft(cov)), size(cov));
-}
-
 
 /*
 Unpack the complex Fourier coefficients of a real Fourier transform with `n` elements to a vector of
@@ -50,15 +38,16 @@ vector gp_unpack_rfft(complex_vector x, int n) {
 /**
 Transform a Gaussian process realization to white noise in the Fourier domain.
 */
-vector gp_transform_rfft(vector y, vector loc, vector rfft_scale) {
-    return gp_unpack_rfft(rfft(y - loc) ./ rfft_scale, size(y));
+vector gp_transform_rfft(vector y, vector loc, vector cov_rfft) {
+    return gp_unpack_rfft(rfft(y - loc) ./ gp_evaluate_rfft_scale(cov_rfft, size(y)), size(y));
 }
 
 
 /**
 Evaluate the log absolute determinant of the Jacobian associated with :stan:func:`gp_transform_rfft`.
 */
-real gp_rfft_log_abs_det_jacobian(vector rfft_scale, int n) {
+real gp_rfft_log_abs_det_jacobian(vector cov_rfft, int n) {
+    vector[n %/% 2 + 1] rfft_scale = gp_evaluate_rfft_scale(cov_rfft, n);
     return - sum(log(rfft_scale[1:n %/% 2 + 1])) -sum(log(rfft_scale[2:(n + 1) %/% 2]))
         - log(2) * ((n - 1) %/% 2) + n * log(n) / 2;
 }
@@ -70,16 +59,15 @@ space.
 
 :param y: Random variable whose likelihood to evaluate.
 :param loc: Mean of the Gaussian process.
-:param cov: Covariance between the origin and the rest of the domain (see
-    :stan:func:`gp_evaluate_rfft_scale(vector)` for details).
+:param cov_rfft: Real fast Fourier transform of the covariance kernel.
 
 :returns: Log probability of the Gaussian process.
 */
-real gp_rfft_lpdf(vector y, vector loc, vector rfft_scale) {
+real gp_rfft_lpdf(vector y, vector loc, vector cov_rfft) {
     int n = size(y);
     int nrfft = n %/% 2 + 1;
-    vector[n] z = gp_transform_rfft(y, loc, rfft_scale);
-    return std_normal_lpdf(z) + gp_rfft_log_abs_det_jacobian(rfft_scale, n);
+    vector[n] z = gp_transform_rfft(y, loc, cov_rfft);
+    return std_normal_lpdf(z) + gp_rfft_log_abs_det_jacobian(cov_rfft, n);
 }
 
 
@@ -114,11 +102,12 @@ with structure expected by the fast Fourier transform. The input vector :math:`z
 
 :param z: Fourier-domain white noise comprising :math:`n` elements.
 :param loc: Mean of the Gaussian process.
-:param cov: Covariance between the origin and the rest of the domain (see
-    :stan:func:`gp_evaluate_rfft_scale(vector)` for details).
+:param cov_rfft: Real fast Fourier transform of the covariance kernel.
 
 :returns: Realization of the Gaussian process with :math:`n` elements.
 */
-vector gp_transform_irfft(vector z, vector loc, vector rfft_scale) {
-    return get_real(inv_rfft(rfft_scale .* gp_pack_rfft(z), size(z))) + loc;
+vector gp_transform_inv_rfft(vector z, vector loc, vector cov_rfft) {
+    int n = size(z);
+    vector[n %/% 2 + 1] rfft_scale = gp_evaluate_rfft_scale(cov_rfft, n);
+    return get_real(inv_rfft(rfft_scale .* gp_pack_rfft(z), n)) + loc;
 }
diff --git a/gptools-stan/gptools/stan/gptools_fft2.stan b/gptools-stan/gptools/stan/gptools_fft2.stan
index 1257637..10fc43c 100644
--- a/gptools-stan/gptools/stan/gptools_fft2.stan
+++ b/gptools-stan/gptools/stan/gptools_fft2.stan
@@ -2,12 +2,12 @@
 // I.e., x[1:3] includes x[1] to x[3]. More generally, x[i:j] comprises j - i + 1 elements. It could
 // at least have been exclusive on the right...
 
-matrix gp_evaluate_rfft2_scale(matrix rfft2_, int width) {
-    int height = rows(rfft2_);
+matrix gp_evaluate_rfft2_scale(matrix cov_rfft2, int width) {
+    int height = rows(cov_rfft2);
     int n = width * height;
     int fftwidth = width %/% 2 + 1;
     int fftheight = height %/% 2 + 1;
-    matrix[height, fftwidth] fftscale = n * rfft2_ / 2;
+    matrix[height, fftwidth] fftscale = n * cov_rfft2 / 2;
     // Check positive-definiteness.
     real minval = min(fftscale);
     if (minval < 0) {
@@ -31,13 +31,6 @@ matrix gp_evaluate_rfft2_scale(matrix rfft2_, int width) {
     return sqrt(fftscale);
 }
 
-/**
-Evaluate the scale of Fourier coefficients.
-*/
-matrix gp_evaluate_rfft2_scale(matrix cov) {
-    return gp_evaluate_rfft2_scale(get_real(rfft2(cov)), cols(cov));
-}
-
 
 matrix gp_unpack_rfft2(complex_matrix z, int m) {
     int height = rows(z);
@@ -82,16 +75,17 @@ complex_matrix gp_pack_rfft2(matrix z) {
 /**
 Transform a Gaussian process realization to white noise in the Fourier domain.
 */
-matrix gp_transform_rfft2(matrix y, matrix loc, matrix rfft2_scale) {
-    return gp_unpack_rfft2(rfft2(y - loc) ./ rfft2_scale, cols(y));
+matrix gp_transform_rfft2(matrix y, matrix loc, matrix cov_rfft2) {
+    return gp_unpack_rfft2(rfft2(y - loc) ./ gp_evaluate_rfft2_scale(cov_rfft2, cols(y)), cols(y));
 }
 
 
 /**
 Transform white noise in the Fourier domain to a Gaussian process realization.
 */
-matrix gp_transform_irfft2(matrix z, matrix loc, matrix rfft2_scale) {
-    complex_matrix[rows(z), cols(z) %/% 2 + 1] y = gp_pack_rfft2(z) .* rfft2_scale;
+matrix gp_transform_inv_rfft2(matrix z, matrix loc, matrix cov_rfft2) {
+    complex_matrix[rows(z), cols(z) %/% 2 + 1] y = gp_pack_rfft2(z)
+        .* gp_evaluate_rfft2_scale(cov_rfft2, cols(z));
     return inv_rfft2(y, cols(z)) + loc;
 }
 
@@ -99,12 +93,12 @@ matrix gp_transform_irfft2(matrix z, matrix loc, matrix rfft2_scale) {
 /**
 Evaluate the log absolute determinant of the Jacobian associated with :stan:func:`gp_transform_rfft`.
 */
-real gp_rfft2_log_abs_det_jacobian(int width, matrix rfft2_scale) {
-    int height = rows(rfft2_scale);
+real gp_rfft2_log_abs_det_jacobian(matrix cov_rfft2, int width) {
+    int height = rows(cov_rfft2);
     int n = width * height;
     int fftwidth = width %/% 2 + 1;
     int fftheight = height %/% 2 + 1;
-    matrix[height, fftwidth] log_rfft2_scale = log(rfft2_scale);
+    matrix[height, fftwidth] log_rfft2_scale = log(gp_evaluate_rfft2_scale(cov_rfft2, width));
     real ladj = 0;
 
     // For the real part, we always use the full height of the non-redundant part. For the imaginary
@@ -137,11 +131,11 @@ Evaluate the log probability of a two-dimensional Gaussian process with zero mea
 
 :param y: Random variable whose likelihood to evaluate.
 :param loc: Mean of the Gaussian process.
-:param cov: First row of the covariance matrix.
+:param cov_rfft2: Real fast Fourier transform of the covariance kernel.
 
 :returns: Log probability of the Gaussian process.
 */
-real gp_rfft2_lpdf(matrix y, matrix loc, matrix rfft2_scale) {
-    return std_normal_lpdf(gp_transform_rfft2(y, loc, rfft2_scale))
-        + gp_rfft2_log_abs_det_jacobian(cols(y), rfft2_scale);
+real gp_rfft2_lpdf(matrix y, matrix loc, matrix cov_rfft2) {
+    return std_normal_lpdf(gp_transform_rfft2(y, loc, cov_rfft2))
+        + gp_rfft2_log_abs_det_jacobian(cov_rfft2, cols(y));
 }
diff --git a/gptools-stan/gptools/stan/gptools_graph.stan b/gptools-stan/gptools/stan/gptools_graph.stan
index bf020c4..140cc83 100644
--- a/gptools-stan/gptools/stan/gptools_graph.stan
+++ b/gptools-stan/gptools/stan/gptools_graph.stan
@@ -1,5 +1,5 @@
 /**
-Evaluate the in-degree of nodes in the directed acyclic graph induced by :code:`edges`. The node
+Evaluate the out-degree of nodes in the directed acyclic graph induced by :code:`edges`. The node
 labels of successors must be ordered and predecessors must have index less than successors.
 
 :param n: Number of nodes.
@@ -9,7 +9,7 @@ labels of successors must be ordered and predecessors must have index less than
 
 :returns: In-degree of each node.
 */
-array [] int in_degrees(int n, array [,] int edge_index) {
+array [] int out_degrees(int n, array [,] int edge_index) {
     array [n] int count = rep_array(0, n);
     int previous = 0;
     for (i in 1:size(edge_index[2])) {
@@ -114,8 +114,7 @@ real gp_graph_exp_quad_cov_lpdf(vector y, vector mu, array [] vector x, real sig
 
 
 /**
-Transform white noise to a sample from a graph Gaussian process with zero mean. A mean can be
-added after the transformation.
+Transform white noise to a sample from a graph Gaussian process with zero mean.
 
 :param z: White noise for each node.
 :param mu: Mean for each node.
@@ -131,9 +130,9 @@ added after the transformation.
 
 :returns: Sample from the Graph gaussian process.
 */
-vector gp_graph_exp_quad_cov_transform(vector z, vector mu, array [] vector x, real sigma,
-                                       real length_scale, array [,] int edges, array [] int degrees,
-                                       real epsilon) {
+vector gp_transform_inv_graph_exp_quad_cov(vector z, vector mu, array [] vector x, real sigma,
+                                           real length_scale, array [,] int edges,
+                                           array [] int degrees, real epsilon) {
     vector[size(z)] y;
     int offset_ = 1;
     for (i in 1:size(x)) {
@@ -144,3 +143,14 @@ vector gp_graph_exp_quad_cov_transform(vector z, vector mu, array [] vector x, r
     }
     return y + mu;
 }
+
+/**
+Transform white noise to a sample from a graph Gaussian process with zero mean. See
+:stan:func:`gp_transform_inv_graph_exp_quad_cov(vector, vector, array [] vector, real, real,
+array [,] int, array [] int, real)` for details.
+*/
+vector gp_transform_inv_graph_exp_quad_cov(vector z, vector mu, array [] vector x, real sigma,
+                                           real length_scale, array [,] int edges) {
+    return gp_transform_inv_graph_exp_quad_cov(z, mu, x, sigma, length_scale, edges,
+                                               out_degrees(size(z), edges), 0.0);
+}
diff --git a/gptools-stan/gptools/stan/gptools_kernels.stan b/gptools-stan/gptools/stan/gptools_kernels.stan
index ce65860..198220a 100644
--- a/gptools-stan/gptools/stan/gptools_kernels.stan
+++ b/gptools-stan/gptools/stan/gptools_kernels.stan
@@ -151,3 +151,8 @@ matrix gp_periodic_matern_cov_rfft2(real dof, int m, int n, real sigma, vector l
         * tgamma(dof + ndim / 2) / tgamma(dof)
         * result .^ -(dof + ndim / 2) * prod(to_array_1d(length_scale ./ period));
 }
+
+matrix gp_periodic_matern_cov_rfft2(real dof, int m, int n, real sigma, real length_scale,
+                                    vector period) {
+    return gp_periodic_matern_cov_rfft2(dof, m, n, sigma, [length_scale, length_scale]', period);
+}
diff --git a/gptools-stan/gptools/stan/gptools_util.stan b/gptools-stan/gptools/stan/gptools_util.stan
index f675b1e..08b36df 100644
--- a/gptools-stan/gptools/stan/gptools_util.stan
+++ b/gptools-stan/gptools/stan/gptools_util.stan
@@ -36,7 +36,8 @@ int is_close(real actual, real desired, real rtol, real atol) {
 }
 
 /**
-Assert that two values are close. See :stan:func:`is_close` for description of parameters.
+Assert that two values are close. See :stan:func:`is_close(real, real, real, real)` for a
+description of the parameters.
 */
 void assert_close(real actual, real desired, real rtol, real atol) {
     if (!is_close(actual, desired, rtol, atol)) {
@@ -45,7 +46,8 @@ void assert_close(real actual, real desired, real rtol, real atol) {
 }
 
 /**
-Assert that two values are close. See :stan:func:`is_close` for description of parameters.
+Assert that two values are close. See :stan:func:`is_close(real, real, real, real)` for a
+description of the parameters.
 */
 void assert_close(real actual, real desired) {
     assert_close(actual, desired, 1e-6, 0);
@@ -70,7 +72,8 @@ int is_finite(complex x) {
 // Vectors -----------------------------------------------------------------------------------------
 
 /**
-Assert that two vectors are close. See :stan:func:`is_close` for description of parameters.
+Assert that two vectors are close. See :stan:func:`is_close(real, real, real, real)` for a
+description of the parameters.
 */
 void assert_close(vector actual, vector desired, real rtol, real atol) {
     int n = size(desired);
diff --git a/gptools-stan/gptools/stan/profile/fourier_centered.stan b/gptools-stan/gptools/stan/profile/fourier_centered.stan
index d723f57..d80dfcc 100644
--- a/gptools-stan/gptools/stan/profile/fourier_centered.stan
+++ b/gptools-stan/gptools/stan/profile/fourier_centered.stan
@@ -13,8 +13,7 @@ parameters {
 }
 
 model {
-    vector[n %/% 2 + 1] rfft_cov = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n, 10);
-    vector[n %/% 2 + 1] rfft_scale = gp_evaluate_rfft_scale(rfft_cov, n);
-    eta ~ gp_rfft(zeros_vector(n), rfft_scale);
+    vector[n %/% 2 + 1] cov_rfft = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n, 10);
+    eta ~ gp_rfft(zeros_vector(n), cov_rfft);
     y[observed_idx] ~ normal(eta[observed_idx], noise_scale);
 }
diff --git a/gptools-stan/gptools/stan/profile/fourier_non_centered.stan b/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
index 0de764e..524696e 100644
--- a/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
+++ b/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
@@ -15,9 +15,8 @@ parameters {
 transformed parameters {
     vector[n] eta;
     {
-        vector[n %/% 2 + 1] rfft_cov = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n, 10);
-        vector[n %/% 2 + 1] rfft_scale = gp_evaluate_rfft_scale(rfft_cov, n);
-        eta = gp_transform_irfft(eta_, zeros_vector(n), rfft_scale);
+        vector[n %/% 2 + 1] cov_rfft = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n, 10);
+        eta = gp_transform_inv_rfft(eta_, zeros_vector(n), cov_rfft);
     }
 }
 
diff --git a/gptools-stan/gptools/stan/profile/graph_centered.stan b/gptools-stan/gptools/stan/profile/graph_centered.stan
index 2edf266..4a46e87 100644
--- a/gptools-stan/gptools/stan/profile/graph_centered.stan
+++ b/gptools-stan/gptools/stan/profile/graph_centered.stan
@@ -7,7 +7,7 @@ functions {
 #include data.stan
 
 transformed data {
-    array [n] int degrees = in_degrees(n, edge_index);
+    array [n] int degrees = out_degrees(n, edge_index);
 }
 
 parameters {
diff --git a/gptools-stan/gptools/stan/profile/graph_non_centered.stan b/gptools-stan/gptools/stan/profile/graph_non_centered.stan
index 108ab50..c600397 100644
--- a/gptools-stan/gptools/stan/profile/graph_non_centered.stan
+++ b/gptools-stan/gptools/stan/profile/graph_non_centered.stan
@@ -7,7 +7,7 @@ functions {
 #include data.stan
 
 transformed data {
-    array [n] int degrees = in_degrees(n, edge_index);
+    array [n] int degrees = out_degrees(n, edge_index);
 }
 
 parameters {
@@ -15,7 +15,7 @@ parameters {
 }
 
 transformed parameters {
-   vector[n] eta = gp_graph_exp_quad_cov_transform(
+   vector[n] eta = gp_transform_inv_graph_exp_quad_cov(
     eta_, zeros_vector(n), X, sigma, length_scale, edge_index, degrees, epsilon);
 }
 
diff --git a/gptools-stan/setup.py b/gptools-stan/setup.py
index 860f71b..540609b 100644
--- a/gptools-stan/setup.py
+++ b/gptools-stan/setup.py
@@ -12,6 +12,8 @@
     packages=find_namespace_packages(),
     version="0.1.0",
     install_requires=[
+        # Required because of a bug in how complex numbers are handled (see
+        # https://github.com/stan-dev/cmdstanpy/pull/612).
         "cmdstanpy>=1.0.7",
         "gp-tools-util",
         "numpy",
diff --git a/gptools-stan/test_requirements.in b/gptools-stan/test_requirements.in
index 03aa22d..14719e3 100644
--- a/gptools-stan/test_requirements.in
+++ b/gptools-stan/test_requirements.in
@@ -1,4 +1,4 @@
 -e file:gptools-stan[tests]
 -e file:gptools-util
 -r ../shared_requirements.in
-cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@develop
+cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@1612c04
diff --git a/gptools-stan/test_requirements.txt b/gptools-stan/test_requirements.txt
index 7388648..98804a6 100644
--- a/gptools-stan/test_requirements.txt
+++ b/gptools-stan/test_requirements.txt
@@ -45,7 +45,7 @@ charset-normalizer==2.1.1
     # via requests
 cloudpickle==2.2.0
     # via doit
-cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@develop
+cmdstanpy @ git+https://github.com/stan-dev/cmdstanpy@1612c04
     # via
     #   -r gptools-stan/test_requirements.in
     #   gp-tools-stan
diff --git a/gptools-stan/tests/test_stan.py b/gptools-stan/tests/test_stan.py
deleted file mode 100644
index e41a0fe..0000000
--- a/gptools-stan/tests/test_stan.py
+++ /dev/null
@@ -1,46 +0,0 @@
-from gptools import stan
-import pathlib
-import pytest
-
-
-@pytest.mark.parametrize("filename", ["gptools_graph.stan", "gptools_fft.stan",
-                                      "gptools_kernels.stan"])
-def test_include(filename: str) -> None:
-    include = pathlib.Path(stan.get_include()) / filename
-    assert include.is_file()
-
-
-def test_needs_compilation(tmp_path: pathlib.Path, caplog: pytest.LogCaptureFixture):
-    # Create files.
-    include_file = tmp_path / "include.stan"
-    with include_file.open("w") as fp:
-        fp.write("real x;")
-    stan_file = tmp_path / "model.stan"
-    with stan_file.open("w") as fp:
-        fp.write("""parameters {
-            #include include.stan
-        }
-        model { x ~ normal(0, 1); }
-        """)
-
-    # Check we compile as usual first.
-    model = stan.compile_model(stan_file=stan_file)
-    exe_file = pathlib.Path(model.exe_file)
-    assert exe_file.exists()
-    last_compiled = exe_file.stat().st_mtime
-
-    # Check the model is not recompiled.
-    stan.compile_model(stan_file=stan_file)
-    assert exe_file.stat().st_mtime == last_compiled
-
-    # Touch the included file and check we compile anew.
-    include_file.touch()
-    stan.compile_model(stan_file=stan_file)
-    assert exe_file.stat().st_mtime > last_compiled
-
-    # Remove the include file and check we get an error.
-    include_file.unlink()
-    with caplog.at_level("DEBUG"):
-        stan.compile_model(stan_file=stan_file)
-    assert any("--info" in record.message and "returned non-zero exit status 1" for record in
-               caplog.records)
diff --git a/gptools-stan/tests/test_stan_functions.py b/gptools-stan/tests/test_stan_functions.py
index 0856c39..87501e4 100644
--- a/gptools-stan/tests/test_stan_functions.py
+++ b/gptools-stan/tests/test_stan_functions.py
@@ -159,13 +159,14 @@ def assert_stan_python_allclose(
         + kernels.DiagonalKernel(.1, 1)
     cov = kernel.evaluate(np.arange(n)[:, None])
     lincov = cov[0]
+    cov_rfft = np.fft.rfft(lincov).real
     rfft_scale = fft.evaluate_rfft_scale(cov=lincov)
     z = fft.transform_rfft(y, loc, rfft_scale=rfft_scale)
     add_configuration({
         "stan_function": "gp_transform_rfft",
         "arg_types": {"n_": "int", "y": "vector[n_]", "loc": "vector[n_]",
-                      "rfft_scale": "vector[n_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "y": y, "loc": loc, "rfft_scale": rfft_scale},
+                      "cov_rfft": "vector[n_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "y": y, "loc": loc, "cov_rfft": cov_rfft},
         "result_type": "vector[n_]",
         "includes": ["gptools_util.stan", "gptools_fft1.stan"],
         "desired": z,
@@ -173,10 +174,10 @@ def assert_stan_python_allclose(
 
     # ... and back again.
     add_configuration({
-        "stan_function": "gp_transform_irfft",
+        "stan_function": "gp_transform_inv_rfft",
         "arg_types": {"n_": "int", "z": "vector[n_]", "loc": "vector[n_]",
-                      "rfft_scale": "vector[n_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "z": z, "loc": loc, "rfft_scale": rfft_scale},
+                      "cov_rfft": "vector[n_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "z": z, "loc": loc, "cov_rfft": cov_rfft},
         "result_type": "vector[n_]",
         "includes": ["gptools_util.stan", "gptools_fft1.stan"],
         "desired": [y, fft.transform_irfft(z, loc, rfft_scale=rfft_scale)],
@@ -186,8 +187,8 @@ def assert_stan_python_allclose(
     add_configuration({
         "stan_function": "gp_rfft_lpdf",
         "arg_types": {"n_": "int", "y": "vector[n_]", "loc": "vector[n_]",
-                      "rfft_scale": "vector[n_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "y": y, "loc": loc, "rfft_scale": rfft_scale},
+                      "cov_rfft": "vector[n_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "y": y, "loc": loc, "cov_rfft": cov_rfft},
         "result_type": "real",
         "includes": ["gptools_util.stan", "gptools_fft1.stan"],
         "desired": [fft.evaluate_log_prob_rfft(y, loc, rfft_scale=rfft_scale),
@@ -245,13 +246,14 @@ def assert_stan_python_allclose(
     xs = coordgrid(np.arange(n), np.arange(m))
     cov = kernel.evaluate(xs)
     lincov = cov[0].reshape((n, m))
+    cov_rfft2 = np.fft.rfft2(lincov).real
     rfft2_scale = fft.evaluate_rfft2_scale(lincov)
     z = fft.transform_rfft2(y, loc, rfft2_scale=rfft2_scale)
     add_configuration({
         "stan_function": "gp_transform_rfft2",
         "arg_types": {"n_": "int", "m_": "int", "y": "matrix[n_, m_]", "loc": "matrix[n_, m_]",
-                      "rfft2_scale": "matrix[n_, m_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "m_": m, "y": y, "loc": loc, "rfft2_scale": rfft2_scale},
+                      "cov_rfft2": "matrix[n_, m_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "m_": m, "y": y, "loc": loc, "cov_rfft2": cov_rfft2},
         "result_type": "matrix[n_, m_]",
         "includes": ["gptools_util.stan", "gptools_fft1.stan", "gptools_fft2.stan"],
         "desired": z,
@@ -259,10 +261,10 @@ def assert_stan_python_allclose(
 
     # ... and back again.
     add_configuration({
-        "stan_function": "gp_transform_irfft2",
+        "stan_function": "gp_transform_inv_rfft2",
         "arg_types": {"n_": "int", "m_": "int", "z": "matrix[n_, m_]", "loc": "matrix[n_, m_]",
-                      "rfft2_scale": "matrix[n_, m_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "m_": m, "z": z, "loc": loc, "rfft2_scale": rfft2_scale},
+                      "cov_rfft2": "matrix[n_, m_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "m_": m, "z": z, "loc": loc, "cov_rfft2": cov_rfft2},
         "result_type": "matrix[n_, m_]",
         "includes": ["gptools_util.stan", "gptools_fft1.stan", "gptools_fft2.stan"],
         "desired": [y, fft.transform_irfft2(z, loc, rfft2_scale=rfft2_scale)],
@@ -272,8 +274,8 @@ def assert_stan_python_allclose(
     add_configuration({
         "stan_function": "gp_rfft2_lpdf",
         "arg_types": {"n_": "int", "m_": "int", "y": "matrix[n_, m_]", "loc": "matrix[n_, m_]",
-                      "rfft2_scale": "matrix[n_, m_ %/% 2 + 1]"},
-        "arg_values": {"n_": n, "m_": m, "y": y, "loc": loc, "rfft2_scale": rfft2_scale},
+                      "cov_rfft2": "matrix[n_, m_ %/% 2 + 1]"},
+        "arg_values": {"n_": n, "m_": m, "y": y, "loc": loc, "cov_rfft2": cov_rfft2},
         "result_type": "real",
         "includes": ["gptools_util.stan", "gptools_fft1.stan", "gptools_fft2.stan"],
         "desired": [stats.multivariate_normal(loc.ravel(), cov).logpdf(y.ravel()),
@@ -485,7 +487,7 @@ def assert_stan_python_allclose(
 ]:
     edges = np.asarray(edges)
     add_configuration({
-        "stan_function": "in_degrees",
+        "stan_function": "out_degrees",
         "arg_types": {"num_nodes": "int", "num_edges_": "int",
                       "edges": "array [2, num_edges_] int"},
         "arg_values": {"num_nodes": num_nodes, "num_edges_": edges.shape[1], "edges": edges},