{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "autoencoder.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# autoencoder \u2014 Python demo\n\nNumerical companion to the entry [autoencoder](https://dictionaryofml.org/terms/autoencoder.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nNumerical companion to the glossary entry 'autoencoder'. An autoencoder learns an encoder and a decoder together, judged by how well the decoder rebuilds a data point from the code the encoder produced. Nothing in that criterion needs a label: what the decoder has to reproduce is the data point itself.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/autoencoder.py`](https://dictionaryofml.org/terms/autoencoder.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(), \"autoencoder.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Autoencoders on a year of days at Krems and on a curved toy dataset:\na linear autoencoder finds the subspace PCA finds, a smaller code costs\nreconstruction error, and a nonlinear autoencoder follows a curve that\nno linear one can.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'autoencoder'.  An autoencoder\nlearns an encoder and a decoder together, judged by how well the decoder\nrebuilds a data point from the code the encoder produced.  Nothing in\nthat criterion needs a label: what the decoder has to reproduce is the\ndata point itself.\n\nTwo datasets.  The GeoSphere Austria weather station Krems (station id\n3805) records eight measurements per day; this script downloads the\nrecords for 2024 from the GeoSphere data hub (dataset klima-v2-1d),\nwrites them to autoencoder_weather.csv, and scales each measurement to\nzero sample mean and unit sample variance because temperatures,\nprecipitation and pressure carry different units.  The curved dataset is\nsynthetic on purpose: points along a parabola in the plane, where the\nstructure a nonlinear autoencoder can follow and a linear one cannot is\nknown in advance.\n\nThe demo checks the entry's claims: (1) a linear autoencoder trained by\ngradient descent on the weather data reaches the reconstruction error\nPCA attains in closed form, and reconstructs into the same subspace, the\nlargest principal angle between the decoder's range and the principal\nsubspace staying below a hundredth of a degree; (2) the code\nsize controls what can be rebuilt, the error falling as the code grows\nand matching the sum of the eigenvalues the code drops; (3) on the\ncurved dataset a nonlinear autoencoder with a code of one number beats\nthe best linear one by more than a factor of ten, because a line cannot\nfollow a parabola; (4) the hand-written gradient used for the training\nagrees with a finite-difference gradient, so the fits rest on a\ngradient that was checked rather than assumed.\n\nDeterministic: the weather data are a fixed archive year, PCA is an\neigenvalue decomposition, and every initialization is drawn from a fixed\nseed.  Self-contained: numpy + matplotlib only (stdlib urllib for the\ndownload).\n\nBlocks\n------\n[B-data]      Download the 366 days with eight measurements each and\n              scale them; build the curved dataset.\n[B-linear]    Train a linear autoencoder (encoder and decoder both\n              matrices) by gradient descent with code size two: check\n              its reconstruction error reaches the PCA minimum, that the\n              decoder reconstructs into the principal subspace, and that\n              its reconstruction map agrees with the PCA projector on\n              the data. The comparison is made on the decoder because\n              the encoder's row space is not determined: the data carry\n              a direction of almost no variance.\n[B-code]      The code size is what forces the choice: the reconstruction\n              error of the best linear autoencoder equals the sum of the\n              dropped eigenvalues, and falls as the code grows.\n[B-gradcheck] The hand-written gradient of the nonlinear autoencoder\n              agrees with a finite-difference gradient to eight digits.\n[B-nonlinear] On the curved dataset, a nonlinear autoencoder with one\n              code number beats the best linear one by more than a\n              factor of ten.\n[B-plot]      Preview: the curved dataset with both reconstructions, and\n              the reconstruction error against the code size.\n\nOutputs\n-------\nautoencoder_weather.csv : date and the eight measurements\nautoencoder_curve.csv   : x1, x2 -- the curved dataset\nautoencoder_linear.csv  : x1, x2 -- its linear (PCA) reconstruction\nautoencoder_nonlinear.csv : x1, x2 -- its nonlinear reconstruction\nautoencoder_codesize.csv  : code, error -- error against code size\nautoencoder.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 = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** Download the 366 days with eight measurements each and scale them; build the curved dataset."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "PARAMS = [\"tlmin\", \"tlmax\", \"tl_mittel\", \"rr\", \"so_h\", \"rf_mittel\",\n          \"p_mittel\", \"vv_mittel\"]\nURL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       f\"klima-v2-1d?parameters={','.join(PARAMS)}&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\"]\nstamps = [t[:10] for t in payload[\"timestamps\"]]\nM = np.stack([np.array(params[p][\"data\"], dtype=float) for p in PARAMS], 1)\nM[:, 3] = np.maximum(M[:, 3], 0.0)             # -1 marks a trace of rain\nwith open(OUT_DIR / \"autoencoder_weather.csv\", \"w\") as f:\n    f.write(\"date,\" + \",\".join(PARAMS) + \"\\n\")\n    for day, row in zip(stamps, M):\n        f.write(day + \",\" + \",\".join(f\"{v:g}\" for v in row) + \"\\n\")\nX = (M - M.mean(axis=0)) / M.std(axis=0)\nm, d = X.shape\ncheck(\"[B-data] 366 days with eight measurements each\", (m, d) == (366, 8))\n\nrng = np.random.default_rng(0)\nt = np.linspace(-1.5, 1.5, 300)\nC = np.stack([t, t ** 2], 1) + 0.03 * rng.standard_normal((300, 2))\nC = C - C.mean(axis=0)                          # centered, for a fair line\ncheck(\"[B-data] the curved dataset lies along a parabola\", len(C) == 300)\n\n\ndef eig_spectrum(A):\n    vals, vecs = np.linalg.eigh(A.T @ A / len(A))\n    order = np.argsort(vals)[::-1]\n    return vals[order], vecs[:, order]\n\n\nlam, U = eig_spectrum(X)\n\n\ndef recon_error(A, W, R):\n    \"\"\"sum_r || a^(r) - R W a^(r) ||^2 for an encoder W and a decoder R.\"\"\"\n    return float(((A - A @ W.T @ R.T) ** 2).sum())"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-linear]** Train a linear autoencoder (encoder and decoder both matrices) by gradient descent with code size two: check its reconstruction error reaches the PCA minimum, that the decoder reconstructs into the principal subspace, and that its reconstruction map agrees with the PCA projector on the data. The comparison is made on the decoder because the encoder's row space is not determined: the data carry a direction of almost no variance."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def train_linear(A, code, steps=4000, lr=0.02, seed=1):\n    g = np.random.default_rng(seed)\n    n, p = A.shape\n    W = 0.1 * g.standard_normal((code, p))\n    R = 0.1 * g.standard_normal((p, code))\n    for _ in range(steps):\n        Z = A @ W.T\n        E = A - Z @ R.T                         # residual, n x p\n        gR = -2.0 * E.T @ Z / n\n        gW = -2.0 * (E @ R).T @ A / n\n        R -= lr * gR\n        W -= lr * gW\n    return W, R\n\n\nCODE = 2\nW_ae, R_ae = train_linear(X, CODE)\nerr_ae = recon_error(X, W_ae, R_ae)\nerr_pca = m * lam[CODE:].sum()\nprint(f\"  linear autoencoder: reconstruction error {err_ae:.1f} against the \"\n      f\"PCA minimum {err_pca:.1f}\")\ncheck(\"[B-linear] the trained linear autoencoder reaches the PCA minimum \"\n      \"(within one percent)\", err_ae < 1.01 * err_pca)\n\n\ndef principal_angles(B1, B2):\n    \"\"\"Angles in degrees between the subspaces spanned by the columns.\"\"\"\n    Q1 = np.linalg.qr(B1)[0]\n    Q2 = np.linalg.qr(B2)[0]\n    s = np.clip(np.linalg.svd(Q1.T @ Q2, compute_uv=False), -1.0, 1.0)\n    return np.degrees(np.arccos(s))\n\n\n# The subspace the autoencoder RECONSTRUCTS into is the decoder's column\n# space, and that is what PCA's principal subspace is compared with. The\n# encoder's row space is not determined here: the measurements carry a\n# direction of almost no variance (see below), and an encoder may read\n# along it without changing any reconstruction.\nang = principal_angles(R_ae, U[:, :CODE])\nang_enc = principal_angles(W_ae.T, U[:, :CODE])\nprint(f\"  largest principal angle to the PCA subspace: decoder \"\n      f\"{ang.max():.4f} degrees, encoder {ang_enc.max():.2f} degrees\")\ncheck(\"[B-linear] the decoder spans the subspace PCA spans (largest \"\n      \"principal angle below a hundredth of a degree)\", ang.max() < 0.01)\ncheck(\"[B-linear] the reconstruction map agrees with the PCA projector on \"\n      \"the data\",\n      np.allclose(X @ (R_ae @ W_ae).T, X @ (U[:, :CODE] @ U[:, :CODE].T),\n                  atol=1e-3))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-code]** The code size is what forces the choice: the reconstruction error of the best linear autoencoder equals the sum of the dropped eigenvalues, and falls as the code grows."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "curve = []\nfor k in range(1, 6):\n    Wk, Rk = train_linear(X, k, seed=2 + k)\n    curve.append((k, recon_error(X, Wk, Rk) / m, m * lam[k:].sum() / m))\nprint(\"  error per data point by code size: \"\n      + \", \".join(f\"{k}: {e:.2f} (PCA {p:.2f})\" for k, e, p in curve))\ncheck(\"[B-code] the error falls as the code grows\",\n      all(curve[i][1] > curve[i + 1][1] for i in range(len(curve) - 1)))\ncheck(\"[B-code] each equals the sum of the dropped eigenvalues (within two \"\n      \"percent)\", all(e < 1.02 * p + 1e-9 for _, e, p in curve))\nwith open(OUT_DIR / \"autoencoder_codesize.csv\", \"w\") as f:\n    f.write(\"code,error,pca\\n\")\n    for k, e, p in curve:\n        f.write(f\"{k},{e:.4f},{p:.4f}\\n\")\n\n\n# ---- the nonlinear autoencoder: one hidden layer on each side\ndef init_nonlinear(p, hidden, seed=3):\n    g = np.random.default_rng(seed)\n    s = 0.8\n    return {\"A1\": s * g.standard_normal((p, hidden)), \"b1\": np.zeros(hidden),\n            \"a2\": s * g.standard_normal((hidden, 1)), \"c2\": np.zeros(1),\n            \"a3\": s * g.standard_normal((1, hidden)), \"b3\": np.zeros(hidden),\n            \"A4\": s * g.standard_normal((hidden, p)), \"b4\": np.zeros(p)}\n\n\ndef forward(P, A):\n    H1 = np.tanh(A @ P[\"A1\"] + P[\"b1\"])         # encoder hidden\n    Z = H1 @ P[\"a2\"] + P[\"c2\"]                  # the code, one number\n    H2 = np.tanh(Z @ P[\"a3\"] + P[\"b3\"])         # decoder hidden\n    Xh = H2 @ P[\"A4\"] + P[\"b4\"]                 # reconstruction\n    return H1, Z, H2, Xh\n\n\ndef loss_and_grad(P, A):\n    n = len(A)\n    H1, Z, H2, Xh = forward(P, A)\n    E = A - Xh\n    loss = float((E ** 2).sum() / n)\n    dXh = -2.0 * E / n\n    g = {\"A4\": H2.T @ dXh, \"b4\": dXh.sum(0)}\n    dH2 = dXh @ P[\"A4\"].T\n    dpre2 = dH2 * (1.0 - H2 ** 2)\n    g[\"a3\"] = Z.T @ dpre2\n    g[\"b3\"] = dpre2.sum(0)\n    dZ = dpre2 @ P[\"a3\"].T\n    g[\"a2\"] = H1.T @ dZ\n    g[\"c2\"] = dZ.sum(0)\n    dH1 = dZ @ P[\"a2\"].T\n    dpre1 = dH1 * (1.0 - H1 ** 2)\n    g[\"A1\"] = A.T @ dpre1\n    g[\"b1\"] = dpre1.sum(0)\n    return loss, g"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-gradcheck]** The hand-written gradient of the nonlinear autoencoder agrees with a finite-difference gradient to eight digits."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "P0 = init_nonlinear(2, 12)\n_, g0 = loss_and_grad(P0, C)\nworst = 0.0\nprobe = np.random.default_rng(4)\nfor name in P0:\n    flat = P0[name].ravel()\n    for _ in range(3):\n        i = int(probe.integers(flat.size))\n        eps, keep = 1e-6, flat[i]\n        flat[i] = keep + eps\n        lp = loss_and_grad(P0, C)[0]\n        flat[i] = keep - eps\n        lm = loss_and_grad(P0, C)[0]\n        flat[i] = keep\n        num = (lp - lm) / (2 * eps)\n        worst = max(worst, abs(num - g0[name].ravel()[i]))\nprint(f\"  gradient check: largest difference to finite differences {worst:.2e}\")\ncheck(\"[B-gradcheck] the hand-written gradient agrees with finite differences\",\n      worst < 1e-8)\n\n\ndef train_nonlinear(A, steps=20000, lr=0.05, seed=3, hidden=12):\n    P = init_nonlinear(A.shape[1], hidden, seed)\n    vel = {k: np.zeros_like(v) for k, v in P.items()}\n    for _ in range(steps):\n        _, g = loss_and_grad(P, A)\n        for k in P:\n            vel[k] = 0.9 * vel[k] - lr * g[k]\n            P[k] = P[k] + vel[k]\n    return P"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-nonlinear]** On the curved dataset, a nonlinear autoencoder with one code number beats the best linear one by more than a factor of ten."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "lamC, UC = eig_spectrum(C)\nw_lin = UC[:, :1].T\nlin_recon = C @ w_lin.T @ w_lin\nerr_lin = float(((C - lin_recon) ** 2).sum() / len(C))\nP = train_nonlinear(C)\nnonlin_recon = forward(P, C)[3]\nerr_non = float(((C - nonlin_recon) ** 2).sum() / len(C))\nprint(f\"  curved dataset, code of one number: linear error {err_lin:.4f}, \"\n      f\"nonlinear error {err_non:.4f} (factor {err_lin / err_non:.1f})\")\ncheck(\"[B-nonlinear] the nonlinear autoencoder beats the best linear one by \"\n      \"more than a factor of ten\", err_lin > 10.0 * err_non)\n\nnp.savetxt(OUT_DIR / \"autoencoder_curve.csv\", C, delimiter=\",\",\n           header=\"x1,x2\", comments=\"\", fmt=\"%.4f\")\nnp.savetxt(OUT_DIR / \"autoencoder_linear.csv\", lin_recon, delimiter=\",\",\n           header=\"x1,x2\", comments=\"\", fmt=\"%.4f\")\norder = np.argsort(nonlin_recon[:, 0])\nnp.savetxt(OUT_DIR / \"autoencoder_nonlinear.csv\", nonlin_recon[order],\n           delimiter=\",\", header=\"x1,x2\", comments=\"\", fmt=\"%.4f\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Preview: the curved dataset with both reconstructions, and the reconstruction error against the code size."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2))\nax1.scatter(C[:, 0], C[:, 1], s=9, color=\"0.6\", label=\"data point\")\nax1.plot(lin_recon[np.argsort(lin_recon[:, 0]), 0],\n         lin_recon[np.argsort(lin_recon[:, 0]), 1], \"k--\", lw=2,\n         label=\"linear autoencoder (the PCA line)\")\nax1.plot(nonlin_recon[order, 0], nonlin_recon[order, 1], \"k-\", lw=2,\n         label=\"nonlinear autoencoder\")\nax1.set_aspect(\"equal\")\nax1.set_xlabel(\"first feature $x_1$\")\nax1.set_ylabel(\"second feature $x_2$\")\nax1.set_title(\"One code number: a line cannot follow a curve\")\nax1.legend(frameon=False, fontsize=8)\nks = [k for k, _, _ in curve]\nax2.plot(ks, [e for _, e, _ in curve], \"ko-\", label=\"trained linear autoencoder\")\nax2.plot(ks, [p for _, _, p in curve], \"s--\", color=\"0.5\",\n         label=\"PCA (sum of the dropped eigenvalues)\")\nax2.set_xlabel(\"code size\")\nax2.set_ylabel(\"reconstruction error per data point\")\nax2.set_title(\"The code size decides what can be rebuilt\")\nax2.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"autoencoder.png\", dpi=110)\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)"
  }
 ]
}