{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "explanation.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# explanation \u2014 Python demo\n\nNumerical companion to the entry [explanation](https://dictionaryofml.org/terms/explanation.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nAn explanation accompanies a prediction, and this demo produces both from real data. The data are hourly precipitation analyses of GeoSphere Austria (INCA, 1 km) over a 120 km box around Krems an der Donau for 8-17 September 2024, the days of storm Boris: 83 mm fell on Krems on 14 September alone (79 mm averaged over the 2 km cell holding the town). They are binned to a 60x60 grid of 2 km cells and committed as explanation_radar.csv, so this script needs no network. (The box was fetched once from the open data hub's grid endpoint, inca-v1-1h-1km, parameters=RR, output_format=netcdf, and binned 2x2 with Krems in the central cell.)\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/explanation.py`](https://dictionaryofml.org/terms/explanation.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(), \"explanation.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nexplanation.py \u2014 numerical companion to the glossary entry 'explanation'.\n\nPurpose\n-------\nAn explanation accompanies a prediction, and this demo produces both from\nreal data.  The data are hourly precipitation analyses of GeoSphere\nAustria (INCA, 1 km) over a 120 km box around Krems an der Donau for\n8-17 September 2024, the days of storm Boris: 83 mm fell on Krems on\n14 September alone (79 mm averaged over the 2 km cell holding the town).  They are binned to a 60x60 grid of 2 km cells and\ncommitted as explanation_radar.csv, so this script needs no network.  (The\nbox was fetched once from the open data hub's grid endpoint,\ninca-v1-1h-1km, parameters=RR, output_format=netcdf, and binned 2x2 with\nKrems in the central cell.)\n\nThe prediction task is the one a reader can check against the sky: given\nthe precipitation over the box at hour t, will it rain at Krems at\nhour t+2?  A small convolutional network is trained on the first days and\ntested on the later ones.  Its explanation is a class activation map\nobtained by Grad-CAM: how strongly the predicted score responds to each\nconvolutional activation map, averaged per map to one weight, then the\nmaps combined with these weights and rectified.  The CAM says which parts of the radar image the network used,\nand the figure shows it by modulating the intensity of the radar image's\npixels: each cell's precipitation is scaled by its CAM score, so the cells\nthe network used stay visible and the rest fade.\n\nSelf-contained (numpy + matplotlib only), fixed seed.\n\nBlocks\n------\n[B-data]  The committed radar frames, read back and checked against the\n          storm they record: 79 mm falls on the Krems cell on 14 September, the\n          wettest hour of the box exceeds 20 mm/h, and the 2-hour-ahead\n          labels are close to balanced. The ten days are chosen so that\n          rain falls in both halves of the chronological split -- with\n          the storm alone, the held-out days are dry and the task there\n          is trivially solved by answering \"no rain\".\n[B-train] A convolutional network (one 3x3 layer, four channels, ReLU,\n          global average pooling, linear head) fitted to the first 70% of\n          the hours by lowering the logistic loss step by step.  Its accuracy on the\n          held-out later hours beats the majority-class baseline.\n[B-cam]   Grad-CAM for one held-out hour: the response of the score to\n          each activation map, pooled to one weight per map, combined and\n          rectified.  The map is written together with the\n          radar image to explanation_frame.csv, which the entry's figure\n          plots.\n[B-text]  The same explanation as one sentence, read off the CAM: the\n          radius around Krems that holds nine in ten of the cells the CAM\n          scores highest (the tenth of the box the faithfulness test\n          perturbs; the last tenth of them sit at the box's corner) and\n          the compass direction of their centre of mass, checked against\n          the precipitation actually recorded there.  The\n          sentence is printed and drawn under the figure's panels.\n[B-faithful] The CAM is tested, not asserted.  Setting the precipitation\n          to zero in the cells the CAM scores highest moves the predicted\n          score much further than zeroing as many cells it scores lowest,\n          and does so for the majority of held-out hours.\n\nOutputs\n-------\nexplanation_frame.csv : x, y, rr (mm/h), cam (0-1), rrcam (mm/h, the\n                        precipitation scaled by the CAM) for the explained hour\nexplanation.png       : preview (checking only) \u2014 the radar image, the\n                        radar image with its intensity modulated by the CAM,\n                        the perturbation test, and the sentence under them\n\"\"\"\n\nimport os\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\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}\")\n\n\nrng = np.random.default_rng(20240914)\nN = 60                              # grid is N x N cells, 2 km each\nKR, KC = 30, 30                     # cell holding Krems an der Donau\nLEAD = 2                            # forecast lead time in hours\nWET = 0.1                           # mm/h counting as rain"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** The committed radar frames, read back and checked against the storm they record: 79 mm falls on the Krems cell on 14 September, the wettest hour of the box exceeds 20 mm/h, and the 2-hour-ahead labels are close to balanced. The ten days are chosen so that rain falls in both halves of the chronological split -- with the storm alone, the held-out days are dry and the task there is trivially solved by answering \"no rain\"."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# The committed radar file is an INPUT, unlike the CSVs the other demos\n# generate, so it is read from beside this script rather than from the\n# working directory: the demo then runs from anywhere, including the\n# throwaway directory the site export uses (which carries an empty\n# pythondemos/ so that OUTPUTS do not overwrite the committed ones).\nRADAR = os.path.join(os.path.dirname(os.path.abspath(__file__)),\n                     \"explanation_radar.csv\")\nncol = len(open(RADAR).readline().split(\",\")) - 2\nif ncol != N * N:\n    raise SystemExit(f\"{RADAR} holds {ncol} cells per hour, this demo expects \"\n                     f\"{N} x {N} = {N * N}: the committed grid and the script \"\n                     f\"disagree\")\nraw = np.loadtxt(RADAR, delimiter=\",\", skiprows=1,\n                 usecols=range(2, 2 + N * N))\nstamps = [ln.split(\",\")[1] for ln in\n          open(RADAR).read().splitlines()[1:]]\nimage = raw.reshape(-1, N, N) / 10.0                 # back to mm/h\nkrems = image[:, KR, KC]\n\ndaily = [krems[d * 24:(d + 1) * 24].sum() for d in range(10)]\ncheck(\"[B-data] 14 September brings over 70 mm to the Krems cell (storm Boris)\",\n      max(daily) > 70.0)\ncheck(\"[B-data] the wettest hour in the box exceeds 20 mm/h\",\n      image.max() > 20.0)\n\nX = image[:-LEAD]                                     # input at hour t\ny = (krems[LEAD:] > WET).astype(float)                # rain at Krems at t+2\ncheck(f\"[B-data] {len(y)} samples, {int(y.sum())} of them with rain ahead\",\n      len(y) == 238 and 0.35 < y.mean() < 0.65)\n\nntr = int(0.7 * len(y))\nXtr, ytr, Xte, yte = X[:ntr], y[:ntr], X[ntr:], y[ntr:]\ncheck(\"[B-data] rain falls in both halves of the split, so neither is trivial\",\n      0.3 < ytr.mean() < 0.7 and 0.3 < yte.mean() < 0.7)\nscale = Xtr.std()\nXtr, Xte = Xtr / scale, Xte / scale\nprint(f\"  train {len(ytr)} hours ({stamps[0]} ..), test {len(yte)} hours \"\n      f\"(.. {stamps[len(y) - 1]})\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-train]** A convolutional network (one 3x3 layer, four channels, ReLU, global average pooling, linear head) fitted to the first 70% of the hours by lowering the logistic loss step by step. Its accuracy on the held-out later hours beats the majority-class baseline."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "K, F = 4, 3                                           # channels, window size\nM = N - F + 1                                         # activation map side\n\n\ndef to_grid(m):\n    \"\"\"Place an M x M activation map back on the N x N input grid.\n\n    The convolution is valid, so activation-map cell (i, j) is computed from\n    input cells i..i+F-1, and the map covers the input minus a one-cell\n    border.  Padding it back into the centre keeps map and field aligned.\n    Upsampling by repetition and cropping to [:N, :N] does NOT: it keeps\n    the upper-left corner of a 2M x 2M array and shifts the explanation\n    away from the weather it explains.\n    \"\"\"\n    out = np.zeros((N, N))\n    pad = (N - M) // 2\n    out[pad:pad + M, pad:pad + M] = m\n    return out\n\n\ndef forward(imgs, W, b, v, c):\n    \"\"\"Conv -> ReLU -> global average pool -> linear.  Returns (score, A).\"\"\"\n    n = len(imgs)\n    A = np.zeros((n, K, M, M))\n    for i in range(F):\n        for j in range(F):\n            patch = imgs[:, i:i + M, j:j + M]         # (n, M, M)\n            A += W[None, :, i, j, None, None] * patch[:, None]\n    A = np.maximum(A + b[None, :, None, None], 0.0)   # (n, K, M, M)\n    g = A.mean(axis=(2, 3))                           # (n, K)\n    return g @ v + c, A\n\n\ndef backward(imgs, A, err, W, v):\n    \"\"\"How the mean logistic loss changes with W, b, v, c (for the update).\"\"\"\n    n = len(imgs)\n    g = A.mean(axis=(2, 3))\n    gv, gc = g.T @ err / n, err.mean()\n    dA = (err[:, None] @ v[None, :])[:, :, None, None] / (M * M)\n    dA = np.repeat(np.repeat(dA, M, axis=2), M, axis=3) * (A > 0)\n    gW = np.zeros_like(W)\n    for i in range(F):\n        for j in range(F):\n            patch = imgs[:, i:i + M, j:j + M]\n            gW[:, i, j] = np.einsum(\"nkxy,nxy->k\", dA, patch) / n\n    return gW, dA.sum(axis=(2, 3)).mean(axis=0), gv, gc\n\n\nW = rng.normal(0, 0.5, (K, F, F))\nb = np.zeros(K)\nv = rng.normal(0, 0.5, K)\nc = 0.0\nlrate, losses = 0.5, []\nfor step in range(600):\n    s, A = forward(Xtr, W, b, v, c)\n    p = 1.0 / (1.0 + np.exp(-s))\n    losses.append(-np.mean(ytr * np.log(p + 1e-9) + (1 - ytr) * np.log(1 - p + 1e-9)))\n    gW, gb, gv, gc = backward(Xtr, A, p - ytr, W, v)\n    W -= lrate * gW\n    b -= lrate * gb\n    v -= lrate * gv\n    c -= lrate * gc\n\nacc = lambda Xs, ys: float((( forward(Xs, W, b, v, c)[0] > 0) == (ys > 0.5)).mean())\nbase = max(yte.mean(), 1 - yte.mean())\nprint(f\"  loss {losses[0]:.3f} -> {losses[-1]:.3f}, \"\n      f\"train acc {acc(Xtr, ytr):.2f}, test acc {acc(Xte, yte):.2f}, \"\n      f\"majority baseline {base:.2f}\")\ncheck(\"[B-train] the loss decreases\", losses[-1] < 0.6 * losses[0])\ncheck(\"[B-train] test accuracy beats the majority-class baseline\",\n      acc(Xte, yte) > base + 0.05)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-cam]** Grad-CAM for one held-out hour: the response of the score to each activation map, pooled to one weight per map, combined and rectified. The map is written together with the radar image to explanation_frame.csv, which the entry's figure plots."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# The explained hour: the held-out hour with the highest predicted score,\n# i.e. the one the network is most sure will bring rain to Krems.\nste, _ = forward(Xte, W, b, v, c)\npick = int(np.argmax(ste))\nimg = Xte[pick]\nscore, A = forward(img[None], W, b, v, c)\nprob = float(1.0 / (1.0 + np.exp(-score[0])))\n\n# Grad-CAM: d(score)/d(A) pooled per activation map gives one weight per map;\n# the map is the rectified weighted sum of the activation maps.\ndA = np.repeat(np.repeat((v[:, None, None] / (M * M))[None], M, axis=2),\n               M, axis=3) * (A > 0)\nalpha = dA[0].mean(axis=(1, 2))                       # one weight per map\ncam = np.maximum((alpha[:, None, None] * A[0]).sum(axis=0), 0.0)\ncam = cam / cam.max() if cam.max() > 0 else cam\ncam_full = to_grid(cam)\n\nhour = stamps[ntr + pick]\ntruth = \"rain\" if yte[pick] > 0.5 else \"no rain\"\nprint(f\"  explained hour {hour} UTC: predicted rain at Krems in {LEAD} h \"\n      f\"with score {prob:.2f}; observed {truth}\")\ncheck(\"[B-cam] the CAM is nonzero, bounded to [0,1] and aligned to the grid\",\n      cam_full.max() == 1.0 and cam_full.min() >= 0.0\n      and cam_full[0, 0] == 0.0)\ncheck(\"[B-cam] the explained hour is one the network predicts as rain\",\n      prob > 0.5)\n\nwith open(OUT_DIR / \"explanation_frame.csv\", \"w\") as fh:\n    fh.write(\"x,y,rr,cam,rrcam\\n\")\n    for r in range(N):\n        for cix in range(N):\n            fh.write(f\"{cix},{r},{img[r, cix] * scale:.2f},\"\n                     f\"{cam_full[r, cix]:.3f},\"\n                     f\"{img[r, cix] * scale * cam_full[r, cix]:.2f}\\n\")\nprint(\"  wrote explanation_frame.csv\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-text]** The same explanation as one sentence, read off the CAM: the radius around Krems that holds nine in ten of the cells the CAM scores highest (the tenth of the box the faithfulness test perturbs; the last tenth of them sit at the box's corner) and the compass direction of their centre of mass, checked against the precipitation actually recorded there. The sentence is printed and drawn under the figure's panels."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# A textual explanation generated from the CAM: where, in distance and\n# direction from Krems, the cells it scores highest lie.  Read off the CAM\n# rather than written freely, so it is faithful by construction.\nCELL_KM = 2.0\nKCELLS = 360                                          # a tenth of the box\nrows, cols = np.indices((N, N))\ndist_km = np.hypot(rows - KR, cols - KC) * CELL_KM\ntop = np.zeros((N, N), bool)\ntop[np.unravel_index(np.argsort(cam_full.ravel())[::-1][:KCELLS], (N, N))] = True\n# the radius that holds nine in ten of those cells, rounded up to 5 km: a\n# few high-scoring cells sit at the box's far corner and would otherwise\n# stretch it to the whole box\nradius_km = float(np.ceil(np.percentile(dist_km[top], 90) / 5.0) * 5.0)\nwts = cam_full[top]\ndn = float(((rows - KR)[top] * wts).sum() / wts.sum())      # north is +row\nde = float(((cols - KC)[top] * wts).sum() / wts.sum())      # east is +col\nbearing = (np.degrees(np.arctan2(de, dn)) + 360.0) % 360.0  # 0 = N, 90 = E\nSECTORS = [\"north\", \"north-north-east\", \"north-east\", \"east-north-east\",\n           \"east\", \"east-south-east\", \"south-east\", \"south-south-east\",\n           \"south\", \"south-south-west\", \"south-west\", \"west-south-west\",\n           \"west\", \"west-north-west\", \"north-west\", \"north-north-west\"]\ndirection = SECTORS[int(((bearing + 11.25) % 360.0) // 22.5)]\nwet_share = float((img[top] * scale > WET).mean())\nsentence = (f\"Rain is likely at Krems in {LEAD} hours since precipitation has \"\n            f\"been recorded within a radius of {radius_km:.0f} km, mainly to \"\n            f\"the {direction} of the town.\")\nprint(f\"  textual explanation: {sentence}\")\nprint(f\"  ({int(top.sum())} highest-scoring cells, bearing {bearing:.0f} deg, \"\n      f\"{100 * wet_share:.0f}% of them wet)\")\ncheck(\"[B-text] the radius holds nine in ten of the highest-scoring cells\",\n      (dist_km[top] <= radius_km).mean() >= 0.9)\ncheck(\"[B-text] precipitation is recorded in most of the cells named\",\n      wet_share > 0.5)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-faithful]** The CAM is tested, not asserted. Setting the precipitation to zero in the cells the CAM scores highest moves the predicted score much further than zeroing as many cells it scores lowest, and does so for the majority of held-out hours."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def zero_cells(image, mask_order, k):\n    out = image.copy()\n    out[np.unravel_index(mask_order[:k], (N, N))] = 0.0\n    return out\n\n\norder_hi = np.argsort(cam_full.ravel())[::-1]         # highest-scoring first\norder_lo = np.argsort(cam_full.ravel())               # lowest-scoring first\ndrop_hi = score[0] - forward(zero_cells(img, order_hi, KCELLS)[None],\n                             W, b, v, c)[0][0]\ndrop_lo = score[0] - forward(zero_cells(img, order_lo, KCELLS)[None],\n                             W, b, v, c)[0][0]\nprint(f\"  zeroing the {KCELLS} highest-scoring cells drops the score by \"\n      f\"{drop_hi:.2f}, the {KCELLS} lowest-scoring by {drop_lo:.2f}\")\ncheck(\"[B-faithful] the highlighted cells move the score further\",\n      drop_hi > 2 * drop_lo)\n\nwins = 0\nfor t in range(len(yte)):\n    it = Xte[t]\n    s_t, A_t = forward(it[None], W, b, v, c)\n    dA_t = np.repeat(np.repeat((v[:, None, None] / (M * M))[None], M, axis=2),\n                     M, axis=3) * (A_t > 0)\n    cm = np.maximum((dA_t[0].mean(axis=(1, 2))[:, None, None] * A_t[0]).sum(0), 0)\n    cm = to_grid(cm)\n    oh, ol = np.argsort(cm.ravel())[::-1], np.argsort(cm.ravel())\n    dh = s_t[0] - forward(zero_cells(it, oh, KCELLS)[None], W, b, v, c)[0][0]\n    dl = s_t[0] - forward(zero_cells(it, ol, KCELLS)[None], W, b, v, c)[0][0]\n    wins += dh > dl\nprint(f\"  zeroing the cells the CAM scores highest lowers the score more than \"\n      f\"zeroing the cells it scores lowest on {wins} of {len(yte)} held-out hours\")\ncheck(\"[B-faithful] the CAM wins on most held-out hours\", wins > 0.7 * len(yte))\n\n# ------------------------------------------------------------------ figure\nfig, ax = plt.subplots(1, 3, figsize=(11.4, 3.6))\next = [0, N, 0, N]\nim0 = ax[0].imshow(img * scale, origin=\"lower\", extent=ext, cmap=\"YlGnBu\",\n                   vmin=0, vmax=max(1.0, (img * scale).max()))\nax[0].set_title(f\"radar: precipitation, {hour} UTC\")\nfig.colorbar(im0, ax=ax[0], label=\"mm/h\", fraction=0.046)\n\n# the CAM shown by modulating the intensity of the radar image's pixels:\n# precipitation scaled by the CAM score, on the same colour scale as ax[0]\nim1 = ax[1].imshow(img * scale * cam_full, origin=\"lower\", extent=ext,\n                   cmap=\"YlGnBu\", vmin=0, vmax=max(1.0, (img * scale).max()))\nax[1].set_title(f\"radar image x CAM, p(rain in {LEAD} h) = {prob:.2f}\")\nfig.colorbar(im1, ax=ax[1], label=\"mm/h, scaled by the CAM\", fraction=0.046)\n\nXL, YL = \"west - east (cells of 2 km)\", \"south - north (cells of 2 km)\"\nax[0].set_xlabel(XL)\nax[0].set_ylabel(YL)\nax[1].set_xlabel(XL)\nax[1].set_ylabel(YL)\n# unrolled rather than looped: check_demo_plots reads the calls statically\n# and cannot tell that a loop variable is one of the axes above\nax[0].plot(KC + 0.5, KR + 0.5, marker=\"*\", ms=13, color=\"crimson\",\n           markeredgecolor=\"k\")\nax[0].annotate(\"Krems\", (KC + 0.5, KR + 0.5), textcoords=\"offset points\",\n               xytext=(7, 5), fontsize=9)\nax[1].plot(KC + 0.5, KR + 0.5, marker=\"*\", ms=13, color=\"crimson\",\n           markeredgecolor=\"k\")\nax[1].annotate(\"Krems\", (KC + 0.5, KR + 0.5), textcoords=\"offset points\",\n               xytext=(7, 5), fontsize=9)\n\nax[2].bar([0, 1], [drop_hi, drop_lo], color=[\"0.25\", \"0.7\"],\n          edgecolor=\"k\", width=0.6)\nax[2].set_xticks([0, 1])\nax[2].set_xticklabels([f\"{KCELLS} highest-scoring\", f\"{KCELLS} lowest-scoring\"],\n                      fontsize=9)\nax[2].set_ylabel(\"drop of the predicted score\")\nax[2].set_xlabel(\"cells set to zero precipitation\")\nax[2].set_title(\"faithfulness: which cells matter\")\nfig.text(0.5, 0.015, \"textual explanation: \" + sentence, ha=\"center\",\n         va=\"bottom\", fontsize=9, style=\"italic\")\nfig.tight_layout(rect=(0, 0.06, 1, 1))\nfig.savefig(OUT_DIR / \"explanation.png\", dpi=110)\nprint(\"  wrote explanation.png\")\n\n# ----------------------------------------------------------------- summary\nbad = [n for n, ok in report if not ok]\nprint(f\"\\n{len(report) - len(bad)}/{len(report)} checks passed\"\n      + (f\"; FAILED: {bad}\" if bad else \"\"))"
  }
 ]
}