Port the legacy GBM script and add its regression test

CHANGELOG.md

@@ -49,3 +49,6 @@ entries here only refer to their errata ID.
 - `asarray_like` to convert problem parameters to the backend of the state.
 - Theory pages for the Euler–Maruyama scheme, geometric Brownian motion and
   the linear system.
+- Port of the legacy script `num_exp_GBM.py` (Euler–Maruyama part) as
+  `examples/strong_error_gbm.py`, and a regression test that reproduces the
+  legacy errors from stored reference data.

docs/development/architecture.md

@@ -616,6 +616,8 @@ the remaining items follow in phase 2 after the lecture.
 | 2026-11-02 | Review of PR #12, P1: supplied and sampled noise pass the same compatibility check (noise dimension, grid starting at `t0`, namespace, dtype and device of an array `x0`); explicit sampling options must agree with an array `x0` |
 | 2026-11-02 | Review of PR #12, P2: noise-structure flags describe the dependence of $g$ on $x$ only and certify no time regularity; `expected_strong_order` states the order under the assumptions of the theory page, and strong order 1 of Euler–Maruyama for additive noise requires $g$ Lipschitz continuous in $t$ |
 | 2026-11-02 | Text files in the repository contain no control characters other than tab, line feed and carriage return; `tests/unit/test_repository_text.py` enforces this, because Sphinx does not warn when LaTeX such as `\big` was turned into a control character |
+| 2026-11-02 | First legacy port: `num_exp_GBM.py` (Euler–Maruyama part). Regression data come from the unmodified legacy function at reduced sample size; the test rebuilds the legacy Brownian paths (generator `default_rng(seed)`, increments of shape `(n_paths, m, n)`, one stream for all step sizes) and wraps them as `BrownianPath`; tolerance `rtol = 1e-12` for platform-dependent rounding of `exp` |
+| 2026-11-02 | Ported examples use the library's coupling (one finest path, coarsened) rather than the legacy practice of drawing new paths per step size; examples run as smoke tests with small parameters |
 | 2026-11-03 | The `gpu` nox session fails unless PyTorch can use a CUDA device (otherwise the GPU tests would only be skipped), reuses its environment because the CUDA wheels take several GB, shows the installation output, and accepts `NUMSDE_TORCH_INDEX` for a PyTorch build that matches the local driver |
 | 2026-11-03 | Literature check of the Euler–Maruyama theory page by the maintainer: Kloeden–Platen Theorem 10.2.2 gives the final-time error in $L^1$, the fundamental theorem of Milstein–Tretyakov (2021, Theorem 1.1.1) the final-time error in $L^2$, Mao (2008, Chapter 2, Theorem 7.3) the error uniformly on the grid for $p = 2$, Kloeden–Platen Section 14.1 the weak order; the divergence for superlinear coefficients follows Hutzenthaler–Jentzen–Kloeden (2011). Kanagawa (1988) for $p \ge 2$ stays unchecked and outside `references.bib` |
 

examples/strong_error_gbm.py

@@ -0,0 +1,166 @@
+r"""Strong error of the Euler–Maruyama scheme for geometric Brownian motion.
+
+Port of the book project's script ``num_exp_GBM.py`` (Euler–Maruyama part;
+the Milstein part follows with milestone M2). The experiment approximates
+
+.. math::
+
+    \mathrm{d}X(t) = \mu X(t)\,\mathrm{d}t + \sigma X(t)\,\mathrm{d}W(t),
+    \qquad X(0) = 1, \quad t \in [0, 1],
+
+with :math:`\mu = 1/4`, :math:`\sigma = 1`, and estimates the strong error
+uniformly on the grid,
+
+.. math::
+
+    e_h = \Big(\mathbb{E}\Big[\max_n \big|X_h^n - X(t_n)\big|^2\Big]\Big)^{1/2},
+
+for the step counts :math:`N = 10 \cdot 2^k`, against the exact solution on
+the same Brownian path.
+
+Unlike the legacy script, which draws new Brownian paths for every step size,
+this version samples one path on the finest grid and coarsens it, so all step
+sizes see the same realization (see the theory page on Brownian motion). The
+estimated errors therefore differ from the legacy ones within the Monte Carlo
+error; the regression test ``tests/regression/test_legacy_gbm.py`` checks
+that the library reproduces the legacy numbers exactly when it is given the
+legacy paths.
+
+Run with the parameters of the book (50 000 paths, about 3 GB of memory)::
+
+    python examples/strong_error_gbm.py
+
+or smaller, with a log-log plot (needs ``numsde[plot]``)::
+
+    python examples/strong_error_gbm.py --n-paths 5000 --plot
+"""
+
+from __future__ import annotations
+
+import argparse
+from dataclasses import dataclass
+from statistics import NormalDist
+
+import numpy as np
+
+import numsde as ns
+
+
+@dataclass(frozen=True)
+class ErrorTable:
+    """Estimated strong errors with 95 % confidence intervals."""
+
+    n_steps: np.ndarray
+    step_sizes: np.ndarray
+    errors: np.ndarray
+    lower: np.ndarray
+    upper: np.ndarray
+
+    @property
+    def eoc(self) -> np.ndarray:
+        """Experimental orders of convergence between consecutive step sizes."""
+        return np.diff(np.log(self.errors)) / np.diff(np.log(self.step_sizes))
+
+    def format(self) -> str:
+        """Return the table as text, in the layout of the legacy output."""
+        lines = [
+            f"{'N':>6} {'h':>10} {'error':>10} {'EOC':>6}   95% conf. interval",
+            "-" * 62,
+        ]
+        eoc = [float("nan"), *self.eoc]
+        for row in zip(
+            self.n_steps,
+            self.step_sizes,
+            self.errors,
+            eoc,
+            self.lower,
+            self.upper,
+            strict=True,
+        ):
+            n, h, e, order, lo, up = row
+            order_text = "" if np.isnan(order) else f"{order:6.2f}"
+            lines.append(
+                f"{n:6d} {h:10.6f} {e:10.6f} {order_text:>6}   [{lo:.6f}, {up:.6f}]"
+            )
+        return "\n".join(lines)
+
+
+def strong_errors(
+    n_paths: int = 50_000, k_max: int = 8, seed: int = 20250718
+) -> ErrorTable:
+    """Estimate the strong error of Euler–Maruyama for ``N = 10 * 2**k``.
+
+    Parameters
+    ----------
+    n_paths
+        Number of Monte Carlo paths.
+    k_max
+        Number of step sizes; ``k = 0, ..., k_max - 1``.
+    seed
+        Seed of the Brownian path.
+
+    Returns
+    -------
+    ErrorTable
+        Errors and confidence intervals, finest step size last.
+    """
+    problem = ns.benchmarks.GeometricBrownianMotion(mu=0.25, sigma=1.0, x0=1.0)
+    solver = ns.Solver(problem, ns.schemes.EulerMaruyama())
+    n_list = np.array([10 * 2**k for k in range(k_max)])
+    finest = ns.WienerProcess(1).sample_path(
+        ns.uniform_grid(1.0, int(n_list[-1])), n_paths=n_paths, seed=seed
+    )
+    z = NormalDist().inv_cdf(0.975)
+
+    errors, lower, upper = [], [], []
+    for n in n_list:
+        path = finest.coarsen(int(n_list[-1] // n))
+        x = solver.solve(path).x
+        exact = problem.exact_solution(path.times[:, None, None], path.values)
+        # Pathwise maximum over the grid of the squared error.
+        y = np.max(np.sum((x - exact) ** 2, axis=-1), axis=0)
+        mean, half = y.mean(), z * y.std(ddof=1) / np.sqrt(n_paths)
+        errors.append(np.sqrt(mean))
+        lower.append(np.sqrt(max(mean - half, 0.0)))
+        upper.append(np.sqrt(mean + half))
+    return ErrorTable(n_list, 1.0 / n_list, *map(np.array, (errors, lower, upper)))
+
+
+def plot(table: ErrorTable) -> None:
+    """Show the errors in a log-log plot with an order line 1/2."""
+    import matplotlib.pyplot as plt
+
+    _, ax = plt.subplots()
+    ax.errorbar(
+        table.n_steps,
+        table.errors,
+        yerr=[table.errors - table.lower, table.upper - table.errors],
+        ls="none",
+        marker="o",
+        label="Euler–Maruyama",
+    )
+    reference = table.errors[0] * (table.n_steps / table.n_steps[0]) ** -0.5
+    ax.loglog(table.n_steps, 0.9 * reference, ":", color="gray", label="order 1/2")
+    ax.set_xlabel("number of time steps $N$")
+    ax.set_ylabel("root-mean-square error, max over the grid")
+    ax.legend()
+    plt.show()
+
+
+def main(argv: list[str] | None = None) -> ErrorTable:
+    """Run the experiment from the command line and print the table."""
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("--n-paths", type=int, default=50_000)
+    parser.add_argument("--k-max", type=int, default=8)
+    parser.add_argument("--seed", type=int, default=20250718)
+    parser.add_argument("--plot", action="store_true")
+    args = parser.parse_args(argv)
+    table = strong_errors(args.n_paths, args.k_max, args.seed)
+    print(table.format())
+    if args.plot:
+        plot(table)
+    return table
+
+
+if __name__ == "__main__":
+    main()

tests/regression/data/README.md

@@ -0,0 +1,47 @@
+# Reference data from legacy scripts
+
+Reference numbers from seeded NumPy scripts of the book project are stored
+here once, together with seed, parameters and the SHA-256 checksums of the
+script files that produced them. The scripts themselves stay with the book
+project and are not part of this repository.
+
+## `legacy_gbm_em_l2error.npz`
+
+| Field | Content |
+| --- | --- |
+| Source | `num_exp_GBM.py`, function `compute_L2_errors`, with `StochODEs.py` and `SODESolvers.py` |
+| Change to the script | none; only the arguments differ from `main()` |
+| Solvers | `EulerMaruyama` only (Milstein follows in milestone M2) |
+| Problem | geometric Brownian motion, `mu = 1/4`, `sigma = 1`, `X0 = 1`, `t_span = [0, 1]` |
+| Seed | 20250718, as in `main()` |
+| Sample size | 2000 paths instead of 50 000, to keep the test fast |
+| Step counts | `N = 10 * 2**k` for `k = 0, ..., 5` instead of `k = 0, ..., 7` |
+| Error | `L2_error()`: (E[max_i \|X(t_i) − X_h(t_i)\|²])^(1/2), the library's default strong error |
+| NumPy | 2.5.4 |
+| Generated | 2026-11-02 |
+
+The test `tests/regression/test_legacy_gbm.py` rebuilds the legacy Brownian
+paths and checks that the library reproduces the errors.
+
+To regenerate the data, run the following with the legacy scripts on the
+path (they need SciPy and matplotlib):
+
+```python
+import numpy as np
+import matplotlib
+
+matplotlib.use("Agg")
+from num_exp_GBM import compute_L2_errors
+from SODESolvers import EulerMaruyama
+from StochODEs import GeometricBrownianMotion
+
+GBM = GeometricBrownianMotion(
+    t_span=[0.0, 1.0],
+    num_samples=2000,
+    seed=20250718,
+    param={"mu": 1 / 4, "sigma": 1.0},
+    X0=1,
+)
+N_list = np.array([10 * 2**k for k in range(6)])
+errors = compute_L2_errors(N_list, GBM, [EulerMaruyama])[0]
+```

tests/regression/test_legacy_gbm.py

@@ -0,0 +1,78 @@
+"""Regression test against the legacy script ``num_exp_GBM.py`` of the book project.
+
+The reference numbers in ``data/legacy_gbm_em_l2error.npz`` were produced by
+the unmodified legacy function ``compute_L2_errors`` with the Euler–Maruyama
+solver only, at a reduced sample size (see ``data/README.md``). The test feeds
+the legacy Brownian paths into the library and requires the same errors.
+
+Legacy sampling convention (``StochODEs.Wiener.simulate``): one generator
+``numpy.random.default_rng(seed)`` serves all step sizes in the order of
+``n_list``; for each step size it draws the increments with shape
+``(n_paths, m, n)`` and sums them cumulatively along time. These are other
+random numbers than :meth:`numsde.WienerProcess.sample_path` draws for the
+same seed, so the paths are rebuilt here and wrapped as
+:class:`numsde.BrownianPath`.
+
+The legacy error is the default strong error of the library, uniformly on the
+grid with p = 2. On the machine that generated the data the results agree
+bitwise; the tolerance covers platform-dependent rounding of ``exp``.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+import numsde as ns
+
+DATA = Path(__file__).parent / "data" / "legacy_gbm_em_l2error.npz"
+
+
+@pytest.fixture(scope="module")
+def reference():
+    with np.load(DATA, allow_pickle=False) as data:
+        return {key: data[key] for key in data.files}
+
+
+def legacy_brownian_paths(seed, n_list, n_paths, t0, t_end):
+    """Yield the Brownian paths of the legacy script, one per step size."""
+    rng = np.random.default_rng(int(seed))
+    for n in n_list:
+        t = np.linspace(t0, t_end, int(n) + 1)
+        dW = rng.normal(loc=0.0, scale=np.sqrt(np.diff(t)), size=(int(n_paths), 1, n))
+        W = np.zeros((int(n_paths), 1, int(n) + 1))
+        W[:, :, 1:] = np.cumsum(dW, axis=2)
+        # Legacy layout (n_paths, m, n + 1) -> library layout (n + 1, n_paths, m).
+        yield ns.BrownianPath(t, np.moveaxis(W, 2, 0), seed=0)
+
+
+def strong_error_on_grid(problem, path, x):
+    """(E[max_n |X_h^n - X(t_n)|^2])^(1/2), computed as in the legacy script."""
+    exact = problem.exact_solution(path.times[:, None, None], path.values)
+    return np.sqrt(np.mean(np.max(np.sum((exact - x) ** 2, axis=-1), axis=0)))
+
+
+def test_reference_data_documents_its_origin(reference):
+    assert str(reference["script"]).startswith("num_exp_GBM.py")
+    assert int(reference["seed"]) == 20250718
+    assert len(str(reference["sha256_num_exp_GBM"])) == 64
+    assert reference["errors"].shape == reference["n_list"].shape
+
+
+def test_euler_maruyama_reproduces_the_legacy_errors(reference):
+    problem = ns.benchmarks.GeometricBrownianMotion(
+        mu=float(reference["mu"]),
+        sigma=float(reference["sigma"]),
+        x0=float(reference["x0"]),
+        t0=float(reference["t0"]),
+    )
+    solver = ns.Solver(problem, ns.schemes.EulerMaruyama())
+    paths = legacy_brownian_paths(
+        reference["seed"],
+        reference["n_list"],
+        reference["n_paths"],
+        float(reference["t0"]),
+        float(reference["t_end"]),
+    )
+    errors = [strong_error_on_grid(problem, p, solver.solve(p).x) for p in paths]
+    np.testing.assert_allclose(errors, reference["errors"], rtol=1e-12, atol=0)

tests/unit/test_examples.py

@@ -0,0 +1,32 @@
+"""Smoke tests: the ported legacy scripts in ``examples/`` run."""
+
+import importlib.util
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+EXAMPLES = Path(__file__).parents[2] / "examples"
+
+
+def load(name):
+    path = EXAMPLES / f"{name}.py"
+    if not path.exists():  # pragma: no cover - only for the installed-wheel job
+        pytest.skip(f"{path} is not available")
+    spec = importlib.util.spec_from_file_location(f"example_{name}", path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+def test_strong_error_gbm_runs_and_prints_a_table(capsys):
+    example = load("strong_error_gbm")
+    table = example.main(["--n-paths", "400", "--k-max", "3", "--seed", "1"])
+    np.testing.assert_array_equal(table.n_steps, [10, 20, 40])
+    assert np.all(table.lower <= table.errors)
+    assert np.all(table.errors <= table.upper)
+    assert np.all(np.diff(table.errors) < 0)
+    assert table.eoc.shape == (2,)
+    output = capsys.readouterr().out
+    assert "95% conf. interval" in output
+    assert len(output.splitlines()) == 5  # header, rule and three rows