Merge pull request #58 from tillahoffmann/interface

Update docs and clean up interface.

.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

.gitignore

@@ -141,3 +141,5 @@ workspace/
 *.o
 *.ipynb
 .cmdstanpy_cache
+*.pdf
+*.png

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

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

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

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

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

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

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,13 +263,35 @@ 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()
 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],
+])
 ```

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.