B-faithful writes when the script runsDictionary of Applied Machine Learning · explanation
Numerical companion to the entry explanation: it recomputes what the entry states and prints one line per check
An 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 48 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. They are binned to a 40x40 grid of 1.2 km cells and committed as explanation_radar.csv, so this script needs no network.
Run it with python3 pythondemos/explanation.py, from the repository root — it writes its data files under pythondemos/. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download explanation.py
One cell per block of the script: the code, and what that code printed when it last ran here
"""
explanation.py — numerical companion to the glossary entry 'explanation'.
Purpose
-------
An 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 48 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. They are binned to a 40x40 grid of 1.2 km cells and committed as
explanation_radar.csv, so this script needs no network.
The prediction task is the one a reader can check against the sky: given
the precipitation over the box at hour t, will it rain at Krems at
hour t+2? A small convolutional network is trained on the first days and
tested on the later ones. Its explanation is a class activation map
obtained by Grad-CAM: the gradient of the predicted score with respect to
the convolutional feature maps, averaged per map to weights, combined and
rectified. The map says which parts of the radar image the network used.
Self-contained (numpy + matplotlib only), fixed seed.
Blocks
------
[B-data] The committed radar frames, read back and checked against the
storm they record: 83 mm falls on Krems 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".
[B-train] A convolutional network (one 3x3 layer, four filters, ReLU,
global average pooling, linear head) trained by gradient
descent on the first 70% of the hours. Its accuracy on the
held-out later hours beats the majority-class baseline.
[B-cam] Grad-CAM for one held-out hour: gradients of the score with
respect to the feature maps, 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.
[B-faithful] The map is tested, not asserted. Setting the precipitation
to zero in the cells the map 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.
Outputs
-------
explanation_frame.csv : x, y, rr (mm/h), cam (0-1) for the explained hour
explanation.png : preview (checking only) — the radar image, the
class activation map, and the perturbation test
"""
import os
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
rng = np.random.default_rng(20240914)
N = 40 # grid is N x N cells, 1.2 km each
KR, KC = 19, 19 # cell holding Krems an der Donau
LEAD = 2 # forecast lead time in hours
WET = 0.1 # mm/h counting as rain
The committed radar frames, read back and checked against the storm they record: 83 mm falls on Krems 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".
# The committed radar file is an INPUT, unlike the CSVs the other demos
# generate, so it is read from beside this script rather than from the
# working directory: the demo then runs from anywhere, including the
# throwaway directory the site export uses (which carries an empty
# pythondemos/ so that OUTPUTS do not overwrite the committed ones).
RADAR = os.path.join(os.path.dirname(os.path.abspath(__file__)),
"explanation_radar.csv")
ncol = len(open(RADAR).readline().split(",")) - 2
if ncol != N * N:
raise SystemExit(f"{RADAR} holds {ncol} cells per hour, this demo expects "
f"{N} x {N} = {N * N}: the committed grid and the script "
f"disagree")
raw = np.loadtxt(RADAR, delimiter=",", skiprows=1,
usecols=range(2, 2 + N * N))
stamps = [ln.split(",")[1] for ln in
open(RADAR).read().splitlines()[1:]]
image = raw.reshape(-1, N, N) / 10.0 # back to mm/h
krems = image[:, KR, KC]
daily = [krems[d * 24:(d + 1) * 24].sum() for d in range(10)]
check("[B-data] 14 September brings over 80 mm to Krems (storm Boris)",
max(daily) > 80.0)
check("[B-data] the wettest hour in the box exceeds 20 mm/h",
image.max() > 20.0)
X = image[:-LEAD] # input at hour t
y = (krems[LEAD:] > WET).astype(float) # rain at Krems at t+2
check(f"[B-data] {len(y)} samples, {int(y.sum())} of them with rain ahead",
len(y) == 238 and 0.35 < y.mean() < 0.65)
ntr = int(0.7 * len(y))
Xtr, ytr, Xte, yte = X[:ntr], y[:ntr], X[ntr:], y[ntr:]
check("[B-data] rain falls in both halves of the split, so neither is trivial",
0.3 < ytr.mean() < 0.7 and 0.3 < yte.mean() < 0.7)
scale = Xtr.std()
Xtr, Xte = Xtr / scale, Xte / scale
print(f" train {len(ytr)} hours ({stamps[0]} ..), test {len(yte)} hours "
f"(.. {stamps[len(y) - 1]})")
[ok] [B-data] 14 September brings over 80 mm to Krems (storm Boris) [ok] [B-data] the wettest hour in the box exceeds 20 mm/h [ok] [B-data] 238 samples, 110 of them with rain ahead [ok] [B-data] rain falls in both halves of the split, so neither is trivial train 166 hours (2024-09-08T00:00 ..), test 72 hours (.. 2024-09-17T21:00) loss 1.176 -> 0.351, train acc 0.86, test acc 0.81, majority baseline 0.53
A convolutional network (one 3x3 layer, four filters, ReLU, global average pooling, linear head) trained by gradient descent on the first 70% of the hours. Its accuracy on the held-out later hours beats the majority-class baseline.
K, F = 4, 3 # filters, kernel size
M = N - F + 1 # feature map side
def to_grid(m):
"""Place an M x M feature map back on the N x N input grid.
The convolution is valid, so feature-map cell (i, j) is computed from
input cells i..i+F-1, and the map covers the input minus a one-cell
border. Padding it back into the centre keeps map and field aligned.
Upsampling by repetition and cropping to [:N, :N] does NOT: it keeps
the upper-left corner of a 2M x 2M array and shifts the explanation
away from the weather it explains.
"""
out = np.zeros((N, N))
pad = (N - M) // 2
out[pad:pad + M, pad:pad + M] = m
return out
def forward(imgs, W, b, v, c):
"""Conv -> ReLU -> global average pool -> linear. Returns (score, A)."""
n = len(imgs)
A = np.zeros((n, K, M, M))
for i in range(F):
for j in range(F):
patch = imgs[:, i:i + M, j:j + M] # (n, M, M)
A += W[None, :, i, j, None, None] * patch[:, None]
A = np.maximum(A + b[None, :, None, None], 0.0) # (n, K, M, M)
g = A.mean(axis=(2, 3)) # (n, K)
return g @ v + c, A
def backward(imgs, A, err, W, v):
"""Gradients of the mean logistic loss w.r.t. W, b, v, c."""
n = len(imgs)
g = A.mean(axis=(2, 3))
gv, gc = g.T @ err / n, err.mean()
dA = (err[:, None] @ v[None, :])[:, :, None, None] / (M * M)
dA = np.repeat(np.repeat(dA, M, axis=2), M, axis=3) * (A > 0)
gW = np.zeros_like(W)
for i in range(F):
for j in range(F):
patch = imgs[:, i:i + M, j:j + M]
gW[:, i, j] = np.einsum("nkxy,nxy->k", dA, patch) / n
return gW, dA.sum(axis=(2, 3)).mean(axis=0), gv, gc
W = rng.normal(0, 0.5, (K, F, F))
b = np.zeros(K)
v = rng.normal(0, 0.5, K)
c = 0.0
lrate, losses = 0.5, []
for step in range(600):
s, A = forward(Xtr, W, b, v, c)
p = 1.0 / (1.0 + np.exp(-s))
losses.append(-np.mean(ytr * np.log(p + 1e-9) + (1 - ytr) * np.log(1 - p + 1e-9)))
gW, gb, gv, gc = backward(Xtr, A, p - ytr, W, v)
W -= lrate * gW
b -= lrate * gb
v -= lrate * gv
c -= lrate * gc
acc = lambda Xs, ys: float((( forward(Xs, W, b, v, c)[0] > 0) == (ys > 0.5)).mean())
base = max(yte.mean(), 1 - yte.mean())
print(f" loss {losses[0]:.3f} -> {losses[-1]:.3f}, "
f"train acc {acc(Xtr, ytr):.2f}, test acc {acc(Xte, yte):.2f}, "
f"majority baseline {base:.2f}")
check("[B-train] the loss decreases", losses[-1] < 0.6 * losses[0])
check("[B-train] test accuracy beats the majority-class baseline",
acc(Xte, yte) > base + 0.05)
[ok] [B-train] the loss decreases [ok] [B-train] test accuracy beats the majority-class baseline explained hour 2024-09-15T08:00 UTC: predicted rain at Krems in 2 h with score 1.00; observed rain
Grad-CAM for one held-out hour: gradients of the score with respect to the feature maps, 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.
# The explained hour: the held-out hour with the highest predicted score,
# i.e. the one the network is most sure will bring rain to Krems.
ste, _ = forward(Xte, W, b, v, c)
pick = int(np.argmax(ste))
img = Xte[pick]
score, A = forward(img[None], W, b, v, c)
prob = float(1.0 / (1.0 + np.exp(-score[0])))
# Grad-CAM: d(score)/d(A) pooled per feature map gives one weight per map;
# the map is the rectified weighted sum of the feature maps.
dA = np.repeat(np.repeat((v[:, None, None] / (M * M))[None], M, axis=2),
M, axis=3) * (A > 0)
alpha = dA[0].mean(axis=(1, 2)) # one weight per map
cam = np.maximum((alpha[:, None, None] * A[0]).sum(axis=0), 0.0)
cam = cam / cam.max() if cam.max() > 0 else cam
cam_full = to_grid(cam)
hour = stamps[ntr + pick]
truth = "rain" if yte[pick] > 0.5 else "no rain"
print(f" explained hour {hour} UTC: predicted rain at Krems in {LEAD} h "
f"with score {prob:.2f}; observed {truth}")
check("[B-cam] the map is nonzero, bounded to [0,1] and aligned to the grid",
cam_full.max() == 1.0 and cam_full.min() >= 0.0
and cam_full[0, 0] == 0.0)
check("[B-cam] the explained hour is one the network predicts as rain",
prob > 0.5)
with open("pythondemos/explanation_frame.csv", "w") as fh:
fh.write("x,y,rr,cam\n")
for r in range(N):
for cix in range(N):
fh.write(f"{cix},{r},{img[r, cix] * scale:.2f},"
f"{cam_full[r, cix]:.3f}\n")
print(" wrote pythondemos/explanation_frame.csv")
[ok] [B-cam] the map is nonzero, bounded to [0,1] and aligned to the grid [ok] [B-cam] the explained hour is one the network predicts as rain wrote pythondemos/explanation_frame.csv zeroing the 160 highest-scoring cells drops the score by 3.32, the 160 lowest-scoring by 0.64
The map is tested, not asserted. Setting the precipitation to zero in the cells the map 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.
def zero_cells(image, mask_order, k):
out = image.copy()
out[np.unravel_index(mask_order[:k], (N, N))] = 0.0
return out
order_hi = np.argsort(cam_full.ravel())[::-1] # highest-scoring first
order_lo = np.argsort(cam_full.ravel()) # lowest-scoring first
KCELLS = 160 # a tenth of the box
drop_hi = score[0] - forward(zero_cells(img, order_hi, KCELLS)[None],
W, b, v, c)[0][0]
drop_lo = score[0] - forward(zero_cells(img, order_lo, KCELLS)[None],
W, b, v, c)[0][0]
print(f" zeroing the {KCELLS} highest-scoring cells drops the score by "
f"{drop_hi:.2f}, the {KCELLS} lowest-scoring by {drop_lo:.2f}")
check("[B-faithful] the highlighted cells move the score further",
drop_hi > 2 * drop_lo)
wins = 0
for t in range(len(yte)):
it = Xte[t]
s_t, A_t = forward(it[None], W, b, v, c)
dA_t = np.repeat(np.repeat((v[:, None, None] / (M * M))[None], M, axis=2),
M, axis=3) * (A_t > 0)
cm = np.maximum((dA_t[0].mean(axis=(1, 2))[:, None, None] * A_t[0]).sum(0), 0)
cm = to_grid(cm)
oh, ol = np.argsort(cm.ravel())[::-1], np.argsort(cm.ravel())
dh = s_t[0] - forward(zero_cells(it, oh, KCELLS)[None], W, b, v, c)[0][0]
dl = s_t[0] - forward(zero_cells(it, ol, KCELLS)[None], W, b, v, c)[0][0]
wins += dh > dl
print(f" the map-guided cells move the score further on {wins} of "
f"{len(yte)} held-out hours")
check("[B-faithful] the map wins on most held-out hours", wins > 0.7 * len(yte))
# ------------------------------------------------------------------ figure
fig, ax = plt.subplots(1, 3, figsize=(11.4, 3.6))
ext = [0, N, 0, N]
im0 = ax[0].imshow(img * scale, origin="lower", extent=ext, cmap="YlGnBu",
vmin=0, vmax=max(1.0, (img * scale).max()))
ax[0].set_title(f"radar: precipitation, {hour} UTC")
fig.colorbar(im0, ax=ax[0], label="mm/h", fraction=0.046)
im1 = ax[1].imshow(cam_full, origin="lower", extent=ext, cmap="YlOrRd",
vmin=0, vmax=1)
ax[1].contour(np.arange(N) + 0.5, np.arange(N) + 0.5, img * scale,
levels=[0.5], colors="k", linewidths=1.0, linestyles="dashed")
ax[1].set_title(f"explanation: Grad-CAM, p(rain in {LEAD} h) = {prob:.2f}")
fig.colorbar(im1, ax=ax[1], label="relevance", fraction=0.046)
XL, YL = "west - east (cells of 1.2 km)", "south - north (cells of 1.2 km)"
ax[0].set_xlabel(XL)
ax[0].set_ylabel(YL)
ax[1].set_xlabel(XL)
ax[1].set_ylabel(YL)
# unrolled rather than looped: check_demo_plots reads the calls statically
# and cannot tell that a loop variable is one of the axes above
ax[0].plot(KC + 0.5, KR + 0.5, marker="*", ms=13, color="crimson",
markeredgecolor="k")
ax[0].annotate("Krems", (KC + 0.5, KR + 0.5), textcoords="offset points",
xytext=(7, 5), fontsize=9)
ax[1].plot(KC + 0.5, KR + 0.5, marker="*", ms=13, color="crimson",
markeredgecolor="k")
ax[1].annotate("Krems", (KC + 0.5, KR + 0.5), textcoords="offset points",
xytext=(7, 5), fontsize=9)
ax[2].bar([0, 1], [drop_hi, drop_lo], color=["0.25", "0.7"],
edgecolor="k", width=0.6)
ax[2].set_xticks([0, 1])
ax[2].set_xticklabels([f"{KCELLS} highest-scoring", f"{KCELLS} lowest-scoring"],
fontsize=9)
ax[2].set_ylabel("drop of the predicted score")
ax[2].set_xlabel("cells set to zero precipitation")
ax[2].set_title("faithfulness: which cells matter")
fig.tight_layout()
fig.savefig("pythondemos/explanation.png", dpi=110)
print(" wrote pythondemos/explanation.png")
# ----------------------------------------------------------------- summary
bad = [n for n, ok in report if not ok]
print(f"\n{len(report) - len(bad)}/{len(report)} checks passed"
+ (f"; FAILED: {bad}" if bad else ""))
[ok] [B-faithful] the highlighted cells move the score further the map-guided cells move the score further on 58 of 72 held-out hours [ok] [B-faithful] the map wins on most held-out hours wrote pythondemos/explanation.png 10/10 checks passed
B-faithful writes when the script runs