{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "sobolevspace.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Sobolev space \u2014 Python demo\n\nNumerical companion to the entry [Sobolev space](https://dictionaryofml.org/terms/sobolevspace.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nRunning example: forecast the maximum temperature of a day at the GeoSphere station Krems (station id 3805) from its minimum temperature. A deep ReLU network learns a non-linear hypothesis map from the 366 days of 2024; a depth-2 decision tree learns a piecewise constant one on the same data. Every claim of the entry is verified on these two maps: the integration-by-parts identity of the weak derivative, the equivalence \"bounded weak gradient <=> Lipschitz\", the robustness reading of that bound, the spectral-norm bound for the network, the absence of any such bound for the tree, and the two regularizers built from Sobolev seminorms. Self-contained (numpy/matplotlib only), fixed seeds.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/sobolevspace.py`](https://dictionaryofml.org/terms/sobolevspace.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"sobolevspace.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nsobolevspace.py \u2014 numerical companion to the glossary entry 'Sobolev space'.\n\nRunning example: forecast the maximum temperature of a day at the\nGeoSphere station Krems (station id 3805) from its minimum temperature.\nA deep ReLU network learns a non-linear hypothesis map from the 366 days\nof 2024; a depth-2 decision tree learns a piecewise constant one on the\nsame data.  Every claim of the entry is verified on these two maps: the\nintegration-by-parts identity of the weak derivative, the equivalence\n\"bounded weak gradient <=> Lipschitz\", the robustness reading of that\nbound, the spectral-norm bound for the network, the absence of any such\nbound for the tree, and the two regularizers built from Sobolev\nseminorms.  Self-contained (numpy/matplotlib only), fixed seeds.\n\nBlocks\n------\n[B-forecast]  Download the 366 days of 2024 at Krems (tlmin, tlmax) and\n              train a ReLU network with three hidden layers of width 8\n              by gradient descent on the squared error loss.  The learned map is\n              non-linear: its slope varies across the feature range.\n[B-weakderiv] The network is piecewise linear, so its classical\n              derivative is undefined at the kinks.  Its slope g, a\n              step function, satisfies the defining identity\n              int f(x) phi'(x) dx = - int g(x) phi(x) dx for infinitely\n              differentiable test functions phi vanishing at both\n              ends: g is the weak\n              derivative of f.\n[B-lipschitz] The largest value of |g| equals the largest difference\n              quotient |f(x) - f(x')| / |x - x'| to three decimals: a\n              bounded weak derivative IS a Lipschitz constant L.  Hence\n              a perturbation delta of the minimum temperature moves the\n              forecast by at most L |delta|; 20000 random perturbations\n              of up to 0.5 degC confirm the bound.\n[B-locallinear] Near almost every minimum temperature x the forecast\n              agrees with the linear map x + delta -> f(x) + g(x) delta:\n              for 20000 random pairs with |delta| <= 0.1 degC the\n              first-order error is exactly zero whenever no kink lies\n              between x and x + delta, and never exceeds 2 L |delta|.\n[B-spectral]  The product of the largest singular values of the weight\n              matrices bounds L (valid but loose).\n[B-jump]      The depth-2 tree is piecewise constant with jumps: a\n              perturbation of 2e-9 degC across its root threshold moves\n              the forecast by the jump height.  No constant L exists.\n[B-nonparam]  Linear regression with a Sobolev penalty (nonparametric\n              regression): the squared error loss over piecewise linear maps on\n              200 knots with the penalty lambda * int f'^2, solved as one\n              linear system.  Raising lambda lowers the H^1 energy of the\n              fit and raises its training error; the fit's slope stays\n              below a bound L of its own.\n[B-tvh1]      Sharpening a transition of width eps: the squared H^1\n              seminorm grows like 4 / (3 eps) while the total variation\n              stays at the height of the limiting jump \u2014 why the L^1\n              gradient penalty tolerates jumps and the H^1 penalty does\n              not.\n[B-graph]     Graph counterpart: for node values sampled from a sine\n              the edge-difference energy (the Laplacian\n              quadratic form used by GTVMin) is small, while node values\n              with one jump make it large.\n\nOutputs\n-------\nsobolevspace_weather.csv : date, tmin, tmax of the 366 downloaded days\nsobolevspace_scatter.csv : tmin, tmax of the 366 days\nsobolevspace_net.csv     : tmin, pred, slope at the breakpoints of the network\nsobolevspace_tree.csv    : tmin, pred -- the tree's step function\nsobolevspace_spline.csv  : tmin, pred -- the Sobolev-penalized linear regression fit\nsobolevspace.png         : preview (checking only)\n\"\"\"\n\nimport json\nimport urllib.request\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nOUT_DIR = Path(__file__).parent\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-forecast]** Download the 366 days of 2024 at Krems (tlmin, tlmax) and train a ReLU network with three hidden layers of width 8 by gradient descent on the squared error loss. The learned map is non-linear: its slope varies across the feature range."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-forecast] download the Krems data and train the network\")\nURL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       \"klima-v2-1d?parameters=tlmin,tlmax&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=120) as resp:\n    payload = json.load(resp)\nparams = payload[\"features\"][0][\"properties\"][\"parameters\"]\nrecords = [(stamp[:10], lo, hi)\n           for stamp, lo, hi in zip(payload[\"timestamps\"],\n                                    params[\"tlmin\"][\"data\"],\n                                    params[\"tlmax\"][\"data\"])\n           if lo is not None and hi is not None]\nwith open(OUT_DIR / \"sobolevspace_weather.csv\", \"w\") as f:\n    f.write(\"date,tmin,tmax\\n\")\n    for day, lo, hi in records:\n        f.write(f\"{day},{lo},{hi}\\n\")\ncheck(f\"[B-forecast] {len(records)} days downloaded for 2024\", len(records) == 366)\n\ntmin = np.array([lo for _, lo, _ in records], dtype=float)\ntmax = np.array([hi for _, _, hi in records], dtype=float)\nnp.savetxt(OUT_DIR / \"sobolevspace_scatter.csv\", np.stack([tmin, tmax], 1),\n           delimiter=\",\", header=\"tmin,tmax\", comments=\"\", fmt=\"%.1f\")\n\nmu_x, sd_x = tmin.mean(), tmin.std()\nmu_y, sd_y = tmax.mean(), tmax.std()\nZ = ((tmin - mu_x) / sd_x)[:, None]           # standardized feature\nY = ((tmax - mu_y) / sd_y)[:, None]           # standardized label\n\nrng = np.random.default_rng(20260918)\nWIDTHS = [1, 8, 8, 8, 1]\nWs = [rng.normal(size=(a, b)) * np.sqrt(2.0 / a)\n      for a, b in zip(WIDTHS[:-1], WIDTHS[1:])]\nbs = [np.zeros(b) for b in WIDTHS[1:]]\n\n\ndef forward(Zin):\n    \"\"\"Return the output and the ReLU masks of the hidden layers.\"\"\"\n    a, masks = Zin, []\n    for W, b in zip(Ws[:-1], bs[:-1]):\n        pre = a @ W + b\n        masks.append(pre > 0)\n        a = np.maximum(pre, 0.0)\n    return a @ Ws[-1] + bs[-1], masks\n\n\n# gradient descent on the average squared error loss over all 366 days\nLR = 0.1\nfor step in range(5000):\n    acts, masks = [Z], []\n    a = Z\n    for W, b in zip(Ws[:-1], bs[:-1]):\n        pre = a @ W + b\n        masks.append(pre > 0)\n        a = np.maximum(pre, 0.0)\n        acts.append(a)\n    out = a @ Ws[-1] + bs[-1]\n    delta = 2.0 * (out - Y) / len(Y)                     # d loss / d out\n    grads_w, grads_b = [], []\n    for layer in range(len(Ws) - 1, -1, -1):\n        grads_w.insert(0, acts[layer].T @ delta)\n        grads_b.insert(0, delta.sum(0))\n        if layer > 0:\n            delta = (delta @ Ws[layer].T) * masks[layer - 1]\n    for k in range(len(Ws)):\n        Ws[k] -= LR * grads_w[k]\n        bs[k] -= LR * grads_b[k]\n\n\ndef f_net(x):\n    \"\"\"Forecast of tmax (degC) for minimum temperatures x (degC).\"\"\"\n    out, _ = forward(((np.asarray(x, dtype=float) - mu_x) / sd_x)[:, None])\n    return out[:, 0] * sd_y + mu_y\n\n\ndef slope_net(x):\n    \"\"\"Weak derivative of f_net at x (degC per degC), read off the masks.\"\"\"\n    _, masks = forward(((np.asarray(x, dtype=float) - mu_x) / sd_x)[:, None])\n    J = np.repeat(Ws[0], len(x), axis=0)                 # (n, width)\n    for k, mask in enumerate(masks):\n        J = (J * mask) @ Ws[k + 1] if k + 1 < len(Ws) - 1 else (J * mask) @ Ws[-1]\n    return J[:, 0] * sd_y / sd_x\n\n\ngrid = np.linspace(tmin.min() - 1.0, tmin.max() + 1.0, 40001)\npred, slope = f_net(grid), slope_net(grid)\nerr_net = float(np.mean((tmax - f_net(tmin)) ** 2))\ncheck(f\"[B-forecast] the learned map is non-linear: its slope ranges from \"\n      f\"{slope.min():.2f} to {slope.max():.2f} degC per degC\",\n      slope.max() - slope.min() > 0.1)\nprint(f\"  average squared error loss of the network on the 366 days: {err_net:.2f}\")\n# the map is piecewise linear: its breakpoints describe it exactly\nkinks = np.flatnonzero(np.diff(slope) != 0.0)\nkeep = np.unique(np.concatenate([[0], kinks, kinks + 1, [len(grid) - 1]]))\nnp.savetxt(OUT_DIR / \"sobolevspace_net.csv\",\n           np.stack([grid[keep], pred[keep], slope[keep]], 1),\n           delimiter=\",\", header=\"tmin,pred,slope\", comments=\"\", fmt=\"%.4f\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-weakderiv]** The network is piecewise linear, so its classical derivative is undefined at the kinks. Its slope g, a step function, satisfies the defining identity int f(x) phi'(x) dx = - int g(x) phi(x) dx for infinitely differentiable test functions phi vanishing at both ends: g is the weak derivative of f."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-weakderiv] integration-by-parts identity for the network\")\na_, b_ = grid[0], grid[-1]\nfor p in (1, 2, 3):\n    # infinitely differentiable test function vanishing at both ends\n    phi = np.sin(p * np.pi * (grid - a_) / (b_ - a_))\n    dphi = p * np.pi / (b_ - a_) * np.cos(p * np.pi * (grid - a_) / (b_ - a_))\n    lhs = np.trapezoid(pred * dphi, grid)\n    rhs = -np.trapezoid(slope * phi, grid)\n    check(f\"[B-weakderiv] int f phi' = -int g phi for phi_{p} \"\n          f\"({lhs:.4f} vs {rhs:.4f})\", abs(lhs - rhs) < 1e-2 * max(1.0, abs(lhs)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-lipschitz]** The largest value of |g| equals the largest difference quotient |f(x) - f(x')| / |x - x'| to three decimals: a bounded weak derivative IS a Lipschitz constant L. Hence a perturbation delta of the minimum temperature moves the forecast by at most L |delta|; 20000 random perturbations of up to 0.5 degC confirm the bound."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-lipschitz] the bound on the weak derivative is a robustness guarantee\")\nL_emp = float(np.abs(slope).max())\ni = rng.integers(0, len(grid), 30000)\nj = rng.integers(0, len(grid), 30000)\nsep = np.abs(grid[i] - grid[j]); ok = sep > 1e-9\nquot = float(np.max(np.abs(pred[i][ok] - pred[j][ok]) / sep[ok]))\ncheck(f\"[B-lipschitz] sup |weak derivative| ({L_emp:.3f}) equals the largest \"\n      f\"difference quotient ({quot:.3f})\", abs(L_emp - quot) < 2e-3)\nx0 = rng.uniform(tmin.min(), tmin.max(), 20000)\ndlt = rng.uniform(-0.5, 0.5, 20000)\nchange = np.abs(f_net(x0 + dlt) - f_net(x0))\ncheck(f\"[B-lipschitz] 20000 perturbations of up to 0.5 degC: every forecast \"\n      f\"change <= L |delta| (largest {change.max():.3f} <= {0.5 * L_emp:.3f})\",\n      np.all(change <= L_emp * np.abs(dlt) + 1e-9))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-locallinear]** Near almost every minimum temperature x the forecast agrees with the linear map x + delta -> f(x) + g(x) delta: for 20000 random pairs with |delta| <= 0.1 degC the first-order error is exactly zero whenever no kink lies between x and x + delta, and never exceeds 2 L |delta|."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-locallinear] first-order approximation of the forecast\")\nx1 = rng.uniform(tmin.min(), tmin.max(), 20000)\nd1 = rng.uniform(-0.1, 0.1, 20000)\nerr1 = np.abs(f_net(x1 + d1) - f_net(x1) - slope_net(x1) * d1)\nk_lo, k_hi = grid[kinks], grid[kinks + 1]             # grid cells holding a kink\nlo_, hi_ = np.minimum(x1, x1 + d1), np.maximum(x1, x1 + d1)\ncrossed = np.array([np.any((k_hi > lo) & (k_lo < hi)) for lo, hi in zip(lo_, hi_)])\ncheck(f\"[B-locallinear] the approximation is exact when no kink lies between x \"\n      f\"and x + delta ({100 * (1 - crossed.mean()):.1f}% of the pairs)\",\n      np.all(err1[~crossed] < 1e-9))\ncheck(f\"[B-locallinear] across a kink the error stays below 2 L |delta| \"\n      f\"(largest ratio {np.max(err1 / (2 * L_emp * np.abs(d1))):.3f})\",\n      np.all(err1 <= 2 * L_emp * np.abs(d1) + 1e-9))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-spectral]** The product of the largest singular values of the weight matrices bounds L (valid but loose)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-spectral] spectral-norm bound on the weak derivative\")\nL_spec = float(np.prod([np.linalg.norm(W, 2) for W in Ws]) * sd_y / sd_x)\ncheck(f\"[B-spectral] product of largest singular values bounds L: \"\n      f\"{L_emp:.3f} <= {L_spec:.3f}\", L_emp <= L_spec + 1e-9)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-jump]** The depth-2 tree is piecewise constant with jumps: a perturbation of 2e-9 degC across its root threshold moves the forecast by the jump height. No constant L exists."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-jump] the tree's forecast jumps\")\nMINLEAF = 5\n\n\ndef best_split(x, y):\n    order = np.argsort(x)\n    xs, ys = x[order], y[order]\n    best, best_sse = None, np.inf\n    for k in range(MINLEAF, len(xs) - MINLEAF + 1):\n        if xs[k - 1] == xs[k]:\n            continue\n        left, right = ys[:k], ys[k:]\n        sse = ((left - left.mean()) ** 2).sum() + ((right - right.mean()) ** 2).sum()\n        if sse < best_sse:\n            best_sse, best = sse, (xs[k - 1] + xs[k]) / 2\n    return best\n\n\nroot_t = best_split(tmin, tmax)\ncuts = [root_t]\nfor side in (tmin <= root_t, tmin > root_t):\n    t = best_split(tmin[side], tmax[side])\n    if t is not None:\n        cuts.append(t)\ncuts = sorted(cuts)\nedges = [-np.inf] + cuts + [np.inf]\nmeans = np.array([tmax[(tmin > lo) & (tmin <= hi)].mean()\n                  for lo, hi in zip(edges[:-1], edges[1:])])\nf_tree = lambda xq: means[np.searchsorted(np.array(cuts), np.asarray(xq, dtype=float))]\neps = 1e-9\njump = float(abs(f_tree([root_t + eps]) - f_tree([root_t - eps]))[0])\ncheck(f\"[B-jump] across the root threshold {root_t:.2f} degC a perturbation of \"\n      f\"{2 * eps:.0e} degC moves the forecast by {jump:.2f} degC \u2014 no constant L\",\n      jump > 1.0)\nbounds = [tmin.min() - 1.0] + cuts + [tmin.max() + 1.0]\nrows = []\nfor lo, hi, mv in zip(bounds[:-1], bounds[1:], means):\n    rows += [(lo, mv), (hi, mv)]\nnp.savetxt(OUT_DIR / \"sobolevspace_tree.csv\", np.array(rows),\n           delimiter=\",\", header=\"tmin,pred\", comments=\"\", fmt=\"%.4f\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-nonparam]** Linear regression with a Sobolev penalty (nonparametric regression): the squared error loss over piecewise linear maps on 200 knots with the penalty lambda * int f'^2, solved as one linear system. Raising lambda lowers the H^1 energy of the fit and raises its training error; the fit's slope stays below a bound L of its own."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-nonparam] squared error loss over piecewise linear maps with the H^1 penalty\")\nknots = np.linspace(tmin.min() - 1.0, tmin.max() + 1.0, 200)\nh_k = knots[1] - knots[0]\nidx = np.clip(np.searchsorted(knots, tmin) - 1, 0, len(knots) - 2)\nwgt = (tmin - knots[idx]) / h_k\nS_mat = np.zeros((len(tmin), len(knots)))          # piecewise linear map at the days\nS_mat[np.arange(len(tmin)), idx] = 1.0 - wgt\nS_mat[np.arange(len(tmin)), idx + 1] = wgt\nD_mat = (np.eye(len(knots))[1:] - np.eye(len(knots))[:-1]) / h_k   # slopes between knots\n\n\ndef spline_fit(lam):\n    \"\"\"Knot values minimizing sum (y - f(x))^2 + lam * int f'(x)^2 dx.\"\"\"\n    A = S_mat.T @ S_mat + lam * h_k * D_mat.T @ D_mat\n    return np.linalg.solve(A, S_mat.T @ tmax)\n\n\nh1_energy = lambda v: float(h_k * np.sum((D_mat @ v) ** 2))\nfits = {lam: spline_fit(lam) for lam in (1.0, 10.0, 100.0)}\nenergies = [h1_energy(v) for v in fits.values()]\nerrors = [float(np.mean((tmax - S_mat @ v) ** 2)) for v in fits.values()]\ncheck(\"[B-nonparam] a larger penalty weight lowers the H^1 energy of the fit \"\n      f\"({', '.join(f'{e:.1f}' for e in energies)})\",\n      all(b < a for a, b in zip(energies, energies[1:])))\ncheck(\"[B-nonparam] and raises its training error \"\n      f\"({', '.join(f'{e:.2f}' for e in errors)})\",\n      all(b > a for a, b in zip(errors, errors[1:])))\nv_sp = fits[10.0]\nL_sp = float(np.abs(D_mat @ v_sp).max())\nprint(f\"  penalty weight 10: training error {errors[1]:.2f}, largest slope {L_sp:.2f} degC per degC\")\nnp.savetxt(OUT_DIR / \"sobolevspace_spline.csv\", np.stack([knots, v_sp], 1),\n           delimiter=\",\", header=\"tmin,pred\", comments=\"\", fmt=\"%.4f\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-tvh1]** Sharpening a transition of width eps: the squared H^1 seminorm grows like 4 / (3 eps) while the total variation stays at the height of the limiting jump \u2014 why the L^1 gradient penalty tolerates jumps and the H^1 penalty does not."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-tvh1] squared H^1 seminorm vs total variation of a sharpening transition\")\nxs = np.linspace(-3.0, 3.0, 600001)\nh1, tv = [], []\nfor width in (0.2, 0.05, 0.0125):\n    df = (1.0 / width) / np.cosh(xs / width) ** 2     # derivative of tanh(x/width)\n    h1.append(np.trapezoid(df ** 2, xs))              # squared H^1 seminorm\n    tv.append(np.trapezoid(np.abs(df), xs))           # total variation\ncheck(\"[B-tvh1] H^1 energy grows without bound as the transition sharpens \"\n      f\"({', '.join(f'{v:.1f}' for v in h1)})\",\n      all(b > 3.5 * a for a, b in zip(h1, h1[1:])))\ncheck(\"[B-tvh1] total variation stays at the height of the jump \"\n      f\"({', '.join(f'{v:.4f}' for v in tv)})\",\n      all(abs(v - 2.0) < 1e-3 for v in tv))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-graph]** Graph counterpart: for node values sampled from a sine the edge-difference energy (the Laplacian quadratic form used by GTVMin) is small, while node values with one jump make it large."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-graph] edge-difference energy on a chain graph\")\nn = 60\nnodes = np.linspace(0.0, 1.0, n)\nchain = [(k, k + 1) for k in range(n - 1)]          # chain of nodes, unit weights\nsmooth = np.sin(2 * np.pi * nodes)                  # samples of a sine\njumpy = np.where(nodes < 0.5, -1.0, 1.0)            # one jump of height 2\nenergy = lambda s: sum((s[a] - s[b]) ** 2 for a, b in chain)\ncheck(f\"[B-graph] edge-difference energy: sine node values {energy(smooth):.4f} \"\n      f\"<< node values with one jump {energy(jumpy):.4f}\",\n      energy(smooth) < 0.1 * energy(jumpy))\n\n# ---- preview figure (checking only)\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.8))\nax1.plot(tmin, tmax, \"o\", ms=2.5, color=\"0.6\", label=\"days of 2024\")\nax1.plot(grid, pred, \"-\", color=\"black\", lw=2, label=\"ReLU network\")\ntr = np.array(rows)\nax1.plot(tr[:, 0], tr[:, 1], \"--\", color=\"black\", lw=1.2, label=\"depth-2 tree\")\nax1.plot(knots, v_sp, \":\", color=\"black\", lw=1.6, label=\"Sobolev-penalized linear regression\")\nax1.set_xlabel(\"minimum temperature of the day (degC)\")\nax1.set_ylabel(\"maximum temperature of the day (degC)\")\nax1.set_title(\"(a) two hypothesis maps for the Krems forecast\")\nax1.legend(frameon=False, loc=\"upper left\")\nax2.plot(grid, slope, \"-\", color=\"black\", lw=2, label=\"weak derivative of the network\")\nax2.axhline(L_emp, ls=\":\", color=\"0.4\", label=f\"L = {L_emp:.2f}\")\nfor c in cuts:\n    ax2.annotate(\"\", xy=(c, 2.6), xytext=(c, 0.0),\n                 arrowprops=dict(arrowstyle=\"->\", color=\"black\", ls=\"--\"))\nax2.plot([], [], \"--\", color=\"black\", label=\"tree: unbounded at its thresholds\")\nax2.set_xlabel(\"minimum temperature of the day (degC)\")\nax2.set_ylabel(\"slope (degC per degC)\")\nax2.set_ylim(-1.0, 2.8)\nax2.set_title(\"(b) slope of the network stays below L\")\nax2.legend(frameon=False, loc=\"lower right\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"sobolevspace.png\", dpi=120)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}