From 19e9a671e9ce10b9d8c1b19b97252a35ab4fcd01 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 11:34:08 -0500
Subject: [PATCH 01/28] Automatically reduce local example test runtime.

---
 docs/test_examples.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

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 = []

From 540e7d538a0ce25a7ba78c3a4f94981ec28f9454 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 12:15:46 -0500
Subject: [PATCH 02/28] Consistent function naming.

---
 gptools-stan/docs/poisson_regression/poisson_regression.md  | 2 +-
 .../poisson_regression_fourier_non_centered.stan            | 2 +-
 .../poisson_regression_graph_non_centered.stan              | 2 +-
 gptools-stan/docs/trees/trees.stan                          | 2 +-
 gptools-stan/docs/tube/tube.stan                            | 2 +-
 gptools-stan/gptools/stan/gptools_fft1.stan                 | 2 +-
 gptools-stan/gptools/stan/gptools_fft2.stan                 | 2 +-
 gptools-stan/gptools/stan/gptools_graph.stan                | 6 +++---
 gptools-stan/gptools/stan/profile/fourier_non_centered.stan | 2 +-
 gptools-stan/gptools/stan/profile/graph_non_centered.stan   | 2 +-
 gptools-stan/tests/test_stan_functions.py                   | 4 ++--
 11 files changed, 14 insertions(+), 14 deletions(-)

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_non_centered.stan b/gptools-stan/docs/poisson_regression/poisson_regression_fourier_non_centered.stan
index 6e73cc1..518b90c 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), gp_evaluate_rfft_scale(cov_rfft, n));
 }
 
 model {
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..9accbf8 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
@@ -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.stan b/gptools-stan/docs/trees/trees.stan
index cb463dd..f5c4793 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -29,7 +29,7 @@ transformed parameters {
         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]');
     // 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),
+    matrix[num_rows_padded, num_cols_padded] f = gp_transform_inv_rfft2(z, mu + zeros_matrix(num_rows_padded, num_cols_padded),
                                              gp_evaluate_rfft2_scale(rfft2_cov, num_cols_padded));
 }
 
diff --git a/gptools-stan/docs/tube/tube.stan b/gptools-stan/docs/tube/tube.stan
index 98593f5..824bc92 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -28,7 +28,7 @@ parameters {
 }
 
 transformed parameters {
-    vector[num_stations] f = gp_graph_exp_quad_cov_transform(
+    vector[num_stations] f = gp_transform_inv_graph_exp_quad_cov(
         z, zeros_vector(num_stations), station_locations, sigma, length_scale, edge_index, degrees, epsilon);
     vector[num_stations] log_mean = mu + f
         + include_zone_effect * one_hot_zones * zone_effect
diff --git a/gptools-stan/gptools/stan/gptools_fft1.stan b/gptools-stan/gptools/stan/gptools_fft1.stan
index c04868a..d427815 100644
--- a/gptools-stan/gptools/stan/gptools_fft1.stan
+++ b/gptools-stan/gptools/stan/gptools_fft1.stan
@@ -119,6 +119,6 @@ with structure expected by the fast Fourier transform. The input vector :math:`z
 
 :returns: Realization of the Gaussian process with :math:`n` elements.
 */
-vector gp_transform_irfft(vector z, vector loc, vector rfft_scale) {
+vector gp_transform_inv_rfft(vector z, vector loc, vector rfft_scale) {
     return get_real(inv_rfft(rfft_scale .* gp_pack_rfft(z), size(z))) + loc;
 }
diff --git a/gptools-stan/gptools/stan/gptools_fft2.stan b/gptools-stan/gptools/stan/gptools_fft2.stan
index 1257637..ebe0165 100644
--- a/gptools-stan/gptools/stan/gptools_fft2.stan
+++ b/gptools-stan/gptools/stan/gptools_fft2.stan
@@ -90,7 +90,7 @@ matrix gp_transform_rfft2(matrix y, matrix loc, matrix rfft2_scale) {
 /**
 Transform white noise in the Fourier domain to a Gaussian process realization.
 */
-matrix gp_transform_irfft2(matrix z, matrix loc, matrix rfft2_scale) {
+matrix gp_transform_inv_rfft2(matrix z, matrix loc, matrix rfft2_scale) {
     complex_matrix[rows(z), cols(z) %/% 2 + 1] y = gp_pack_rfft2(z) .* rfft2_scale;
     return inv_rfft2(y, cols(z)) + loc;
 }
diff --git a/gptools-stan/gptools/stan/gptools_graph.stan b/gptools-stan/gptools/stan/gptools_graph.stan
index bf020c4..08d2a2a 100644
--- a/gptools-stan/gptools/stan/gptools_graph.stan
+++ b/gptools-stan/gptools/stan/gptools_graph.stan
@@ -131,9 +131,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)) {
diff --git a/gptools-stan/gptools/stan/profile/fourier_non_centered.stan b/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
index 0de764e..168befb 100644
--- a/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
+++ b/gptools-stan/gptools/stan/profile/fourier_non_centered.stan
@@ -17,7 +17,7 @@ transformed parameters {
     {
         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);
+        eta = gp_transform_inv_rfft(eta_, zeros_vector(n), rfft_scale);
     }
 }
 
diff --git a/gptools-stan/gptools/stan/profile/graph_non_centered.stan b/gptools-stan/gptools/stan/profile/graph_non_centered.stan
index 108ab50..a7413dc 100644
--- a/gptools-stan/gptools/stan/profile/graph_non_centered.stan
+++ b/gptools-stan/gptools/stan/profile/graph_non_centered.stan
@@ -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/tests/test_stan_functions.py b/gptools-stan/tests/test_stan_functions.py
index 0856c39..f3cfab5 100644
--- a/gptools-stan/tests/test_stan_functions.py
+++ b/gptools-stan/tests/test_stan_functions.py
@@ -173,7 +173,7 @@ 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},
@@ -259,7 +259,7 @@ 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},

From 78215db33cae98ac6a9be12bd833614320432b57 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 12:36:27 -0500
Subject: [PATCH 03/28] Ignore `*.pdf` globally.

---
 .gitignore         | 1 +
 figures/.gitignore | 1 -
 2 files changed, 1 insertion(+), 1 deletion(-)
 delete mode 100644 figures/.gitignore

diff --git a/.gitignore b/.gitignore
index ce0076b..9e32ab0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -141,3 +141,4 @@ workspace/
 *.o
 *.ipynb
 .cmdstanpy_cache
+*.pdf
diff --git a/figures/.gitignore b/figures/.gitignore
deleted file mode 100644
index a136337..0000000
--- a/figures/.gitignore
+++ /dev/null
@@ -1 +0,0 @@
-*.pdf

From 4d325f8c925a30b27f29ff4c3a61c1c08d58364f Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 12:36:58 -0500
Subject: [PATCH 04/28] Exclude figures from docs.

---
 conf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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

From 8d2b7cb1ab7b9d042c9f6466f6b07e882e7afe7f Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 12:37:18 -0500
Subject: [PATCH 05/28] Make sphinx build nitpicky.

---
 dodo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dodo.py b/dodo.py
index 8dfcdab..e324c74 100644
--- a/dodo.py
+++ b/dodo.py
@@ -57,7 +57,7 @@
 
 # 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="html", actions=["sphinx-build -n . docs/_build"])
     manager(name="tests", actions=["sphinx-build -b doctest . docs/_build"])
 
 # Compile example notebooks to create html reports.

From 38ea15c19293fc19038ae7d4d204a9d8a689be84 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 12:38:08 -0500
Subject: [PATCH 06/28] Update docs and example code.

---
 figures/profile.md                           |  2 ++
 gptools-stan/README.rst                      |  1 +
 gptools-stan/docs/trees/trees.md             |  8 +++++---
 gptools-stan/docs/tube/tube.md               | 10 ++++++----
 gptools-stan/docs/tube/tube.stan             |  2 +-
 gptools-stan/gptools/stan/gptools_graph.stan | 14 ++++++++++++--
 gptools-stan/gptools/stan/gptools_util.stan  |  9 ++++++---
 7 files changed, 33 insertions(+), 13 deletions(-)

diff --git a/figures/profile.md b/figures/profile.md
index 13195c1..9c259e8 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
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/trees/trees.md b/gptools-stan/docs/trees/trees.md
index 15aad64..e19f29c 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
@@ -89,8 +91,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 +169,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$")
diff --git a/gptools-stan/docs/tube/tube.md b/gptools-stan/docs/tube/tube.md
index f78d479..7e39a8c 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
@@ -93,7 +95,7 @@ data = {
 ```{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())
 ```
@@ -208,7 +210,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 +227,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")
diff --git a/gptools-stan/docs/tube/tube.stan b/gptools-stan/docs/tube/tube.stan
index 824bc92..5598daa 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -29,7 +29,7 @@ parameters {
 
 transformed parameters {
     vector[num_stations] f = gp_transform_inv_graph_exp_quad_cov(
-        z, zeros_vector(num_stations), station_locations, sigma, length_scale, edge_index, degrees, epsilon);
+        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;
diff --git a/gptools-stan/gptools/stan/gptools_graph.stan b/gptools-stan/gptools/stan/gptools_graph.stan
index 08d2a2a..2641c5b 100644
--- a/gptools-stan/gptools/stan/gptools_graph.stan
+++ b/gptools-stan/gptools/stan/gptools_graph.stan
@@ -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.
@@ -144,3 +143,14 @@ vector gp_transform_inv_graph_exp_quad_cov(vector z, vector mu, array [] vector
     }
     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,
+                                               in_degrees(size(z), edges), 0.0);
+}
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);

From eaf08e6fa301e86297577e733508b335ba7134a0 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 14:07:33 -0500
Subject: [PATCH 07/28] Simplify `rfft` interface for Stan.

---
 .../poisson_regression_fourier_centered.stan  |  2 +-
 ...isson_regression_fourier_non_centered.stan |  2 +-
 gptools-stan/gptools/stan/gptools_fft1.stan   | 41 +++++++------------
 .../stan/profile/fourier_centered.stan        |  5 +--
 .../stan/profile/fourier_non_centered.stan    |  5 +--
 gptools-stan/tests/test_stan_functions.py     | 13 +++---
 6 files changed, 28 insertions(+), 40 deletions(-)

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 518b90c..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_inv_rfft(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/gptools/stan/gptools_fft1.stan b/gptools-stan/gptools/stan/gptools_fft1.stan
index d427815..f702bfd 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 covariancer 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: Real fast Fourier transform of the covariance kernel.
 
 :returns: Realization of the Gaussian process with :math:`n` elements.
 */
-vector gp_transform_inv_rfft(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/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 168befb..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_inv_rfft(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/tests/test_stan_functions.py b/gptools-stan/tests/test_stan_functions.py
index f3cfab5..97d7337 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,
@@ -175,8 +176,8 @@ def assert_stan_python_allclose(
     add_configuration({
         "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),

From d630ad1472b7980aa9cb0e8c00b08fdfde211c30 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 14:12:54 -0500
Subject: [PATCH 08/28] Update `cmdstanpy` for better handling of errors in
 included files.

---
 dev_requirements.txt               | 2 +-
 gptools-stan/setup.py              | 2 ++
 gptools-stan/test_requirements.in  | 2 +-
 gptools-stan/test_requirements.txt | 2 +-
 4 files changed, 5 insertions(+), 3 deletions(-)

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/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

From c116b10722f739ca04315393cbb16c2438ddf158 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 14:44:44 -0500
Subject: [PATCH 09/28] Fix compilation test after upgrading `cmdstanpy`.

---
 gptools-stan/tests/test_stan.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/gptools-stan/tests/test_stan.py b/gptools-stan/tests/test_stan.py
index e41a0fe..e682025 100644
--- a/gptools-stan/tests/test_stan.py
+++ b/gptools-stan/tests/test_stan.py
@@ -10,7 +10,7 @@ def test_include(filename: str) -> None:
     assert include.is_file()
 
 
-def test_needs_compilation(tmp_path: pathlib.Path, caplog: pytest.LogCaptureFixture):
+def test_needs_compilation(tmp_path: pathlib.Path):
     # Create files.
     include_file = tmp_path / "include.stan"
     with include_file.open("w") as fp:
@@ -40,7 +40,5 @@ def test_needs_compilation(tmp_path: pathlib.Path, caplog: pytest.LogCaptureFixt
 
     # Remove the include file and check we get an error.
     include_file.unlink()
-    with caplog.at_level("DEBUG"):
+    with pytest.raises(ValueError):
         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)

From 6169cf958da677acfc4830d831eb692000a8c058 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 14:50:24 -0500
Subject: [PATCH 10/28] Expand doc tests.

---
 dodo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dodo.py b/dodo.py
index e324c74..b2770e7 100644
--- a/dodo.py
+++ b/dodo.py
@@ -58,7 +58,7 @@
 # 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 -n . docs/_build"])
-    manager(name="tests", actions=["sphinx-build -b doctest . 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"):

From 649534c7ce3e4e61e5962d73b531e19006002556 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 14:51:55 -0500
Subject: [PATCH 11/28] Simplify `rfft2` interface for Stan.

---
 gptools-stan/docs/trees/trees.stan          |  4 +--
 gptools-stan/gptools/stan/gptools_fft2.stan | 34 +++++++++------------
 gptools-stan/tests/test_stan_functions.py   | 13 ++++----
 3 files changed, 23 insertions(+), 28 deletions(-)

diff --git a/gptools-stan/docs/trees/trees.stan b/gptools-stan/docs/trees/trees.stan
index f5c4793..546ce71 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -29,8 +29,8 @@ transformed parameters {
         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]');
     // Transform from white-noise to a Gaussian process realization.
-    matrix[num_rows_padded, num_cols_padded] f = gp_transform_inv_rfft2(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, mu + zeros_matrix(num_rows_padded, num_cols_padded), rfft2_cov);
 }
 
 model {
diff --git a/gptools-stan/gptools/stan/gptools_fft2.stan b/gptools-stan/gptools/stan/gptools_fft2.stan
index ebe0165..5f09f77 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_inv_rfft2(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_inv_rfft2(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
@@ -141,7 +135,7 @@ Evaluate the log probability of a two-dimensional Gaussian process with zero mea
 
 :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/tests/test_stan_functions.py b/gptools-stan/tests/test_stan_functions.py
index 97d7337..9f2afa2 100644
--- a/gptools-stan/tests/test_stan_functions.py
+++ b/gptools-stan/tests/test_stan_functions.py
@@ -246,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,
@@ -262,8 +263,8 @@ def assert_stan_python_allclose(
     add_configuration({
         "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)],
@@ -273,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()),

From 40acf3861ff9e9f83cc14f35fdfe3d6b19aec922 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 15:33:37 -0500
Subject: [PATCH 12/28] Minor FFT documentation fixes.

---
 gptools-stan/gptools/stan/gptools_fft1.stan | 4 ++--
 gptools-stan/gptools/stan/gptools_fft2.stan | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/gptools-stan/gptools/stan/gptools_fft1.stan b/gptools-stan/gptools/stan/gptools_fft1.stan
index f702bfd..19f0837 100644
--- a/gptools-stan/gptools/stan/gptools_fft1.stan
+++ b/gptools-stan/gptools/stan/gptools_fft1.stan
@@ -59,7 +59,7 @@ space.
 
 :param y: Random variable whose likelihood to evaluate.
 :param loc: Mean of the Gaussian process.
-:param cov_rfft: Real fast Fourier transform of the covariancer kernel.
+:param cov_rfft: Real fast Fourier transform of the covariance kernel.
 
 :returns: Log probability of the Gaussian process.
 */
@@ -102,7 +102,7 @@ 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: Real fast Fourier transform of the covariance kernel.
+:param cov_rfft: Real fast Fourier transform of the covariance kernel.
 
 :returns: Realization of the Gaussian process with :math:`n` elements.
 */
diff --git a/gptools-stan/gptools/stan/gptools_fft2.stan b/gptools-stan/gptools/stan/gptools_fft2.stan
index 5f09f77..10fc43c 100644
--- a/gptools-stan/gptools/stan/gptools_fft2.stan
+++ b/gptools-stan/gptools/stan/gptools_fft2.stan
@@ -131,7 +131,7 @@ 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.
 */

From 15f9b3806c0ab7d176494bd6c04b11c2f56c40de Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 15:34:35 -0500
Subject: [PATCH 13/28] Rename `in_degrees` to `out_degrees`.

We evaluate the number of connections from successors to predecessors.
---
 .../poisson_regression_graph_centered.stan                  | 2 +-
 .../poisson_regression_graph_non_centered.stan              | 2 +-
 gptools-stan/docs/tube/tube.stan                            | 2 +-
 gptools-stan/gptools/stan/gptools_graph.stan                | 6 +++---
 gptools-stan/gptools/stan/profile/graph_centered.stan       | 2 +-
 gptools-stan/gptools/stan/profile/graph_non_centered.stan   | 2 +-
 gptools-stan/tests/test_stan_functions.py                   | 2 +-
 7 files changed, 9 insertions(+), 9 deletions(-)

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 9accbf8..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 {
diff --git a/gptools-stan/docs/tube/tube.stan b/gptools-stan/docs/tube/tube.stan
index 5598daa..3a8f89f 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -16,7 +16,7 @@ data {
 }
 
 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 {
diff --git a/gptools-stan/gptools/stan/gptools_graph.stan b/gptools-stan/gptools/stan/gptools_graph.stan
index 2641c5b..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])) {
@@ -152,5 +152,5 @@ 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,
-                                               in_degrees(size(z), edges), 0.0);
+                                               out_degrees(size(z), edges), 0.0);
 }
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 a7413dc..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 {
diff --git a/gptools-stan/tests/test_stan_functions.py b/gptools-stan/tests/test_stan_functions.py
index 9f2afa2..87501e4 100644
--- a/gptools-stan/tests/test_stan_functions.py
+++ b/gptools-stan/tests/test_stan_functions.py
@@ -487,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},

From 48817e4c34bd504839dde7696e15dd26efc2f89b Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 16:02:03 -0500
Subject: [PATCH 14/28] Add isotropic `rfft` evaluation for Matern kernels.

---
 gptools-stan/docs/trees/trees.stan             | 10 +++++-----
 gptools-stan/gptools/stan/gptools_kernels.stan |  5 +++++
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/gptools-stan/docs/trees/trees.stan b/gptools-stan/docs/trees/trees.stan
index 546ce71..2fd0a8c 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -20,17 +20,17 @@ 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, length_scale, kappa;
 }
 
 transformed parameters {
     // 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_inv_rfft2(
-        z, mu + zeros_matrix(num_rows_padded, num_cols_padded), rfft2_cov);
+        z, rep_matrix(mu, num_rows_padded, num_cols_padded), rfft2_cov);
 }
 
 model {
@@ -38,7 +38,7 @@ model {
     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);
             }
         }
     }
@@ -55,5 +55,5 @@ model {
     // 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);
+    kappa ~ cauchy(0, 1);
 }
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);
+}

From d37072b07447356014ac20916fb762a50c4d4255 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 16:14:03 -0500
Subject: [PATCH 15/28] Use keyword arguments for forwards compatibility.

---
 gptools-stan/gptools/stan/__init__.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

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__":

From 704928a74671657448b1ceb34e482f763e383b4e Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 22 Nov 2022 16:15:20 -0500
Subject: [PATCH 16/28] Remove compile test as it is now in `cmdstanpy` (see
 stan-dev/cmdstanpy#638).

---
 gptools-stan/tests/test_stan.py | 44 ---------------------------------
 1 file changed, 44 deletions(-)
 delete mode 100644 gptools-stan/tests/test_stan.py

diff --git a/gptools-stan/tests/test_stan.py b/gptools-stan/tests/test_stan.py
deleted file mode 100644
index e682025..0000000
--- a/gptools-stan/tests/test_stan.py
+++ /dev/null
@@ -1,44 +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):
-    # 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 pytest.raises(ValueError):
-        stan.compile_model(stan_file=stan_file)

From 364c640cee3b7775f1028bc54117770f499e5ba3 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Mon, 5 Dec 2022 16:46:49 -0500
Subject: [PATCH 17/28] Use log-uniform distribution for tree length scale.

---
 gptools-stan/docs/trees/trees.md   | 3 +++
 gptools-stan/docs/trees/trees.stan | 7 +++++--
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/gptools-stan/docs/trees/trees.md b/gptools-stan/docs/trees/trees.md
index e19f29c..a219711 100644
--- a/gptools-stan/docs/trees/trees.md
+++ b/gptools-stan/docs/trees/trees.md
@@ -41,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
@@ -63,6 +64,8 @@ data = {
     "num_cols_padded": padded_cols,
     "frequency": np.where(train_mask, frequency, -1).astype(int),
     "epsilon": 0,
+    "length_scale_lower": 2,
+    "length_scale_upper": 12.5,
 }
 padded_rows, padded_cols
 ```
diff --git a/gptools-stan/docs/trees/trees.stan b/gptools-stan/docs/trees/trees.stan
index 2fd0a8c..9341918 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -11,7 +11,8 @@ data {
     // Number of trees in each quadrant. Masked quadrants are indicated by a negative value.
     array [num_rows, num_cols] int frequency;
     // Nugget variance.
-    real<lower=0> epsilon;
+    real<lower=0> epsilon, length_scale_lower;
+    real<lower=length_scale_lower> length_scale_upper;
 }
 
 parameters {
@@ -20,10 +21,12 @@ 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, kappa;
+    real<lower=0> sigma, kappa;
+    real<lower=log(length_scale_lower), upper=log(length_scale_upper)> log_length_scale;
 }
 
 transformed parameters {
+    real 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,

From 76954580f31cab2a3d9dc47a6dfc52ab181017cd Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Mon, 5 Dec 2022 16:47:26 -0500
Subject: [PATCH 18/28] Use log-uniform distribution for tube length scale.

---
 gptools-stan/docs/tube/tube.md   |  4 ++++
 gptools-stan/docs/tube/tube.stan | 11 ++++++-----
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/gptools-stan/docs/tube/tube.md b/gptools-stan/docs/tube/tube.md
index 7e39a8c..1260c10 100644
--- a/gptools-stan/docs/tube/tube.md
+++ b/gptools-stan/docs/tube/tube.md
@@ -23,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")
@@ -76,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],
@@ -89,6 +91,8 @@ data = {
     "epsilon": 0,
     "include_zone_effect": 1,
     "include_degree_effect": 1,
+    "length_scale_lower": 2 * distances.min(),
+    "length_scale_upper": distances.max() / 2,
 }
 ```
 
diff --git a/gptools-stan/docs/tube/tube.stan b/gptools-stan/docs/tube/tube.stan
index 3a8f89f..80b4488 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -9,10 +9,10 @@ data {
     array [2, num_edges] int edge_index;
     matrix[num_stations, num_zones] one_hot_zones;
     matrix[num_stations, num_degrees] one_hot_degrees;
-    real<lower=0> epsilon;
+    real<lower=0> epsilon, length_scale_lower;
+    real<lower=length_scale_lower> length_scale_upper;
 
-    int include_zone_effect;
-    int include_degree_effect;
+    int include_zone_effect, include_degree_effect;
 }
 
 transformed data {
@@ -22,12 +22,14 @@ transformed data {
 parameters {
     vector[num_stations] z;
     real mu;
-    real<lower=0> sigma, length_scale, kappa;
+    real<lower=0> sigma, kappa;
+    real<lower=log(length_scale_lower), upper=log(length_scale_upper)> log_length_scale;
     vector[num_zones] zone_effect;
     vector[num_degrees] degree_effect;
 }
 
 transformed parameters {
+    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
@@ -39,7 +41,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);

From 4357d8785a7b4742022de6319f6983677dc03907 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Mon, 5 Dec 2022 16:48:07 -0500
Subject: [PATCH 19/28] Add decision tree figure for practitioners.

---
 figures/decision_tree.md | 94 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 94 insertions(+)
 create mode 100644 figures/decision_tree.md

diff --git a/figures/decision_tree.md b/figures/decision_tree.md
new file mode 100644
index 0000000..9022aff
--- /dev/null
+++ b/figures/decision_tree.md
@@ -0,0 +1,94 @@
+---
+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")
+```

From c5ad26b993b4049ec92dc131d69a12888c58da65 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 13 Dec 2022 13:15:32 -0500
Subject: [PATCH 20/28] Ignore .png files.

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 9e32ab0..d143762 100644
--- a/.gitignore
+++ b/.gitignore
@@ -142,3 +142,4 @@ workspace/
 *.ipynb
 .cmdstanpy_cache
 *.pdf
+*.png

From a306619fff2472d0b075931851f9b750a25e96a0 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 13 Dec 2022 13:17:17 -0500
Subject: [PATCH 21/28] Save figures as .png for presentations.

---
 figures/decision_tree.md         | 1 +
 figures/profile.md               | 1 +
 gptools-stan/docs/trees/trees.md | 1 +
 gptools-stan/docs/tube/tube.md   | 1 +
 4 files changed, 4 insertions(+)

diff --git a/figures/decision_tree.md b/figures/decision_tree.md
index 9022aff..247be9a 100644
--- a/figures/decision_tree.md
+++ b/figures/decision_tree.md
@@ -91,4 +91,5 @@ 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/profile.md b/figures/profile.md
index 9c259e8..6757ce8 100644
--- a/figures/profile.md
+++ b/figures/profile.md
@@ -267,4 +267,5 @@ legend = fig.legend(*zip(*handles_labels), fontsize="small", loc="upper center",
                     ncol=5, bbox_to_anchor=bbox_to_anchor)
 
 fig.savefig("scaling.pdf", bbox_inches="tight")
+fig.savefig("scaling.png", bbox_inches="tight")
 ```
diff --git a/gptools-stan/docs/trees/trees.md b/gptools-stan/docs/trees/trees.md
index a219711..e1aa5c3 100644
--- a/gptools-stan/docs/trees/trees.md
+++ b/gptools-stan/docs/trees/trees.md
@@ -195,4 +195,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/tube/tube.md b/gptools-stan/docs/tube/tube.md
index 1260c10..b9b8e73 100644
--- a/gptools-stan/docs/tube/tube.md
+++ b/gptools-stan/docs/tube/tube.md
@@ -205,6 +205,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.

From 3e665f80ccc0661e1f646ec6f19d985a6f41421a Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 13 Dec 2022 13:17:46 -0500
Subject: [PATCH 22/28] Move color bar inside panel.

---
 figures/profile.md | 38 ++++++++++++++++++++++++--------------
 1 file changed, 24 insertions(+), 14 deletions(-)

diff --git a/figures/profile.md b/figures/profile.md
index 6757ce8..9227205 100644
--- a/figures/profile.md
+++ b/figures/profile.md
@@ -127,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])
 
@@ -163,7 +162,6 @@ for i, ax in [(0, ax1), (-1, ax2)]:
     ax.set_xlabel("size $n$")
     ax.set_ylabel("duration (seconds)")
 
-
 ax = ax3
 ax.set_xlabel(r"noise scale $\kappa$")
 ax.set_ylabel(r"runtime (seconds)")
@@ -179,7 +177,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$")
@@ -206,23 +205,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)
@@ -258,7 +255,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()
@@ -266,6 +263,19 @@ 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")
 ```

From 7f6683f8769dbfd12b73a5c145630816cf22c2e6 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 13 Dec 2022 13:20:34 -0500
Subject: [PATCH 23/28] Log uniform prior for tree example length scale.

---
 gptools-stan/docs/trees/trees.md   |  6 ++----
 gptools-stan/docs/trees/trees.stan | 27 +++++++++------------------
 2 files changed, 11 insertions(+), 22 deletions(-)

diff --git a/gptools-stan/docs/trees/trees.md b/gptools-stan/docs/trees/trees.md
index e1aa5c3..7627589 100644
--- a/gptools-stan/docs/trees/trees.md
+++ b/gptools-stan/docs/trees/trees.md
@@ -64,8 +64,6 @@ data = {
     "num_cols_padded": padded_cols,
     "frequency": np.where(train_mask, frequency, -1).astype(int),
     "epsilon": 0,
-    "length_scale_lower": 2,
-    "length_scale_upper": 12.5,
 }
 padded_rows, padded_cols
 ```
@@ -73,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())
 ```
 
diff --git a/gptools-stan/docs/trees/trees.stan b/gptools-stan/docs/trees/trees.stan
index 9341918..d795db1 100644
--- a/gptools-stan/docs/trees/trees.stan
+++ b/gptools-stan/docs/trees/trees.stan
@@ -11,8 +11,7 @@ data {
     // Number of trees in each quadrant. Masked quadrants are indicated by a negative value.
     array [num_rows, num_cols] int frequency;
     // Nugget variance.
-    real<lower=0> epsilon, length_scale_lower;
-    real<lower=length_scale_lower> length_scale_upper;
+    real<lower=0> epsilon;
 }
 
 parameters {
@@ -22,11 +21,11 @@ parameters {
     real mu;
     // Kernel parameters and averdispersion parameter for the negative binomial distribution.
     real<lower=0> sigma, kappa;
-    real<lower=log(length_scale_lower), upper=log(length_scale_upper)> log_length_scale;
+    real<lower=log(2), upper=log(28)> log_length_scale;
 }
 
 transformed parameters {
-    real length_scale = exp(log_length_scale);
+    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,
@@ -37,6 +36,12 @@ transformed parameters {
 }
 
 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) {
@@ -45,18 +50,4 @@ model {
             }
         }
     }
-    // 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.
-    kappa ~ cauchy(0, 1);
 }

From 4754161cc0f51aead1775f66f7c97058a304782f Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Tue, 13 Dec 2022 13:21:02 -0500
Subject: [PATCH 24/28] Align priors and use log uniform length scale prior for
 tube example.

---
 gptools-stan/docs/tube/tube.md              | 5 ++---
 gptools-stan/docs/tube/tube.stan            | 5 ++---
 gptools-stan/docs/tube/tube_without_gp.stan | 6 +++---
 3 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/gptools-stan/docs/tube/tube.md b/gptools-stan/docs/tube/tube.md
index b9b8e73..af58fd5 100644
--- a/gptools-stan/docs/tube/tube.md
+++ b/gptools-stan/docs/tube/tube.md
@@ -91,9 +91,8 @@ data = {
     "epsilon": 0,
     "include_zone_effect": 1,
     "include_degree_effect": 1,
-    "length_scale_lower": 2 * distances.min(),
-    "length_scale_upper": distances.max() / 2,
 }
+distances.min(), distances.max()
 ```
 
 ```{code-cell} ipython3
@@ -246,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 80b4488..21980cb 100644
--- a/gptools-stan/docs/tube/tube.stan
+++ b/gptools-stan/docs/tube/tube.stan
@@ -9,8 +9,7 @@ data {
     array [2, num_edges] int edge_index;
     matrix[num_stations, num_zones] one_hot_zones;
     matrix[num_stations, num_degrees] one_hot_degrees;
-    real<lower=0> epsilon, length_scale_lower;
-    real<lower=length_scale_lower> length_scale_upper;
+    real<lower=0> epsilon;
 
     int include_zone_effect, include_degree_effect;
 }
@@ -23,7 +22,7 @@ parameters {
     vector[num_stations] z;
     real mu;
     real<lower=0> sigma, kappa;
-    real<lower=log(length_scale_lower), upper=log(length_scale_upper)> log_length_scale;
+    real<lower=log(0.32), upper=log(31)> log_length_scale;
     vector[num_zones] zone_effect;
     vector[num_degrees] degree_effect;
 }
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);

From 1f5c5d3dd4c7c32fc1d077fe1d22ce85ebf81369 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Wed, 14 Dec 2022 21:48:54 -0500
Subject: [PATCH 25/28] Print log-log slope for fit to last three finite
 runtimes.

---
 figures/profile.md | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/figures/profile.md b/figures/profile.md
index 9227205..f3ebd9e 100644
--- a/figures/profile.md
+++ b/figures/profile.md
@@ -147,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(
@@ -155,12 +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$")
@@ -279,3 +287,11 @@ 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],
+])
+```

From 5e0cfe09919a585654fa3aafedfecf0ca6408a4c Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Wed, 14 Dec 2022 21:49:31 -0500
Subject: [PATCH 26/28] Remove whitespace.

---
 figures/decision_tree.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/figures/decision_tree.md b/figures/decision_tree.md
index 247be9a..6b88d5b 100644
--- a/figures/decision_tree.md
+++ b/figures/decision_tree.md
@@ -74,7 +74,7 @@ def connect(xy1, xy2, frac=0.25, color="k", label=None):
             }
         }
         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")

From 7de6bcc235224a2fe3d11e6b8cbe2e7ee3fcd619 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Wed, 14 Dec 2022 21:49:55 -0500
Subject: [PATCH 27/28] Add notebook to illustrate different kernels.

---
 figures/kernels.md | 87 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 87 insertions(+)
 create mode 100644 figures/kernels.md

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")
+```

From 6aed7648f069c51adccddd7dccb106c3b02925c4 Mon Sep 17 00:00:00 2001
From: Till Hoffmann <tillahoffmann@gmail.com>
Date: Wed, 14 Dec 2022 22:04:07 -0500
Subject: [PATCH 28/28] Remove duplicate notebook tests.

---
 .github/workflows/main.yml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

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