{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "dataaug.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# data augmentation \u2014 Python demo\n\nNumerical companion to the entry [data augmentation](https://dictionaryofml.org/terms/dataaug.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 'dataaug'. The GeoSphere Austria weather station Krems (station id 3805) records the maximum air temperature of each day. This script downloads the 2024 records from the GeoSphere data hub (dataset klima-v2-1d) and writes them to dataaug_maxtemp.csv.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/dataaug.py`](https://dictionaryofml.org/terms/dataaug.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(), \"dataaug.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Three measured days, a degree-five polynomial that interpolates them,\nand what happens when each day is replaced by copies that lie within the\naccuracy of the thermometer.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'dataaug'.  The GeoSphere\nAustria weather station Krems (station id 3805) records the maximum air\ntemperature of each day.  This script downloads the 2024 records from\nthe GeoSphere data hub (dataset klima-v2-1d) and writes them to\ndataaug_maxtemp.csv.\n\nThe learning task pairs consecutive days: the feature of a data point is\nthe maximum temperature measured today, its label is the maximum\ntemperature measured tomorrow.  Three days far apart in the year form\nthe training set, so the three features are far apart as well.\n\nThree data points leave a polynomial of degree five underdetermined: it\nhas six coefficients, so infinitely many such polynomials pass through\nall three and every one of them has zero training error.  ERM delivers any of them and the training error cannot choose.\nThe curve drawn is one of them, built by adding to the smoothest\ninterpolant a degree-five polynomial that vanishes at the three measured\nfeatures, scaled so that it reaches 15 degrees Celsius.  That curve\nleaves the range the three labels occupy by far, which is what\noverfitting looks like here.\n\nA thermometer reports the air temperature only to within its accuracy,\nassumed here to be 0.5 degrees Celsius.  Both the feature and the label\nof a data point are readings of that instrument, so moving each of them\nanywhere inside that tolerance gives a day that is just as consistent\nwith what was measured.  Replacing each of the three days by copies\ndrawn that way is the augmented training set.  The copies carry many\ndistinct feature values, so the feature matrix of the augmented set has\nfull column rank and ERM on it has a unique\nsolution.  A curve that is steep where the copies lie pays for it, since\na copy displaced in the feature direction then misses its label, so the\nsteep interpolants are the ones the augmented criterion rejects.\n\nThe demo checks the claims the entry makes.  (1) The three training days\nlie in different seasons and their features span a wide range.  (2) Both\nthe smoothest interpolant and the drawn one have zero training error on\nthe three days, so the training error does not separate them.  (3) The\ndrawn one leaves the range of the three labels by more than ten degrees\nbetween them.  (4) Every copy lies within the assumed accuracy of the\nday it came from, and the augmented feature matrix has full column rank,\nso the augmented problem has a unique solution.  (5) That solution stays\nwithin one degree of the range of the three labels, an excursion smaller\nby more than a factor of ten.  (6) On the day pairs of 2024 that were\nheld out and whose feature lies between the smallest and the largest\ntraining feature, it has a far smaller average squared error than the\ndrawn interpolant, and it lands close to the straight line through the\nthree days.\n\nDeterministic: the data are a fixed archive year and the copies are\ndrawn with a fixed seed.  Self-contained: numpy + matplotlib only\n(stdlib urllib for the download).\n\nBlocks\n------\n[B-data]      Download the maximum temperature of every day of 2024 and\n              pair each day with the next one.\n[B-three]     Pick the three training days, one per season, and report\n              their features and labels.\n[B-overfit]   Two of the infinitely many ERM solutions: both with zero\n              training error, differing widely between the three days.\n[B-augment]   Replace each day by copies inside the accuracy of the\n              thermometer and fit the same degree-five model to them.\n[B-compare]   The fits on the held-out day pairs whose feature lies\n              between the smallest and the largest training feature,\n              and the straight line through the three days.\n[B-fig]       The preview figure: the picture on the left, the held-out\n              errors on the right.\n\nOutputs\n-------\ndataaug_maxtemp.csv : date and maximum temperature, 366 days\ndataaug_points.csv  : the three training days, feature and label\ndataaug_copies.csv  : the copies drawn inside the sensor accuracy\ndataaug_curves.csv  : the two fitted curves on a grid of features\ndataaug_zoom.csv    : copies of the middle day, relative to that day\ndataaug_valerr.csv  : held-out error of the three fits\ndataaug.png         : preview figure\n\"\"\"\n\nimport json\nimport urllib.request\nfrom pathlib import Path\n\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt                        # noqa: E402\nimport numpy as np                                     # noqa: E402\n\nOUT_DIR = Path(__file__).parent\nDEGREE = 5\nTOL = 0.5                 # assumed accuracy of the thermometer, in Celsius\nNR_COPIES = 200           # copies drawn per training day\nSWING = 15.0              # amplitude, in Celsius, of the drawn ERM solution\nRNG = np.random.default_rng(20240101)\n\nFAILED = []\n\n\ndef check(name, ok):\n    FAILED.append(name) if not ok else None\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef design(x, lo, hi):\n    \"\"\"Powers of the feature, rescaled to [-1, 1] so the fit is stable.\"\"\"\n    u = 2.0 * (x - lo) / (hi - lo) - 1.0\n    return np.stack([u ** k for k in range(DEGREE + 1)], axis=1)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** Download the maximum temperature of every day of 2024 and pair each day with the next one."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "URL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       \"klima-v2-1d?parameters=tlmax&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=180) as resp:\n    payload = json.load(resp)\ntmax = np.array(\n    payload[\"features\"][0][\"properties\"][\"parameters\"][\"tlmax\"][\"data\"],\n    dtype=float)\ndates = [t[:10] for t in payload[\"timestamps\"]]\nwith open(OUT_DIR / \"dataaug_maxtemp.csv\", \"w\") as f:\n    f.write(\"date,tlmax\\n\")\n    for day, v in zip(dates, tmax):\n        f.write(f\"{day},{v:g}\\n\")\n# today's maximum predicts tomorrow's maximum\nFEAT, LAB = tmax[:-1], tmax[1:]\nprint(f\"[B-data] {len(FEAT)} day pairs of 2024 at Krems an der Donau\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-three]** Pick the three training days, one per season, and report their features and labels."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "IDX = [20, 105, 200]                      # late January, mid April, late July\nxt, yt = FEAT[IDX], LAB[IDX]\nfor i, d in zip(IDX, [dates[i] for i in IDX]):\n    print(f\"[B-three] {d}: today {FEAT[i]:.1f} C, tomorrow {LAB[i]:.1f} C\")\ncheck(\"[B-three] the three training days lie in different seasons\",\n      len({dates[i][5:7] for i in IDX}) == 3)\ncheck(\"[B-three] their features span more than twenty degrees\",\n      xt.max() - xt.min() > 20.0)\n\nLO, HI = FEAT.min(), FEAT.max()\nGRID = np.linspace(xt.min(), xt.max(), 400)\nc_lin = np.polyfit(xt, yt, 1)\nLINE = np.polyval(c_lin, GRID)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-overfit]** Two of the infinitely many ERM solutions: both with zero training error, differing widely between the three days."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "A = design(xt, LO, HI)\nc_smooth = np.linalg.lstsq(A, yt, rcond=None)[0]   # smallest coefficients\n# every polynomial that vanishes at the three features can be added to it\nnull = np.linalg.svd(A)[2][len(IDX):]              # directions spanning it\nwiggle = design(GRID, LO, HI) @ null[0]\nc_int = c_smooth + (SWING / np.max(np.abs(wiggle))) * null[0]\nfit_smooth = design(GRID, LO, HI) @ c_smooth\nfit_int = design(GRID, LO, HI) @ c_int\ntrain_smooth = np.mean((A @ c_smooth - yt) ** 2)\ntrain_int = np.mean((A @ c_int - yt) ** 2)\nexc_int = max(fit_int.max() - yt.max(), yt.min() - fit_int.min())\nprint(f\"[B-overfit] {null.shape[0]} independent directions leave the three \"\n      f\"labels untouched, so ERM has infinitely many solutions\")\nprint(f\"[B-overfit] training error of the smoothest {train_smooth:.2e} and \"\n      f\"of the drawn one {train_int:.2e}\")\ngap = np.max(np.abs(fit_int - fit_smooth))\nMIDX = 12.0                       # a mild spring day, between the measured ones\nmid = design(np.array([MIDX]), LO, HI)\nspan = np.max(np.abs(mid @ null.T))\nprint(f\"[B-overfit] between the measured days the two solutions differ by \"\n      f\"up to {gap:.1f} C\")\nprint(f\"[B-overfit] at a today-maximum of {MIDX:.0f} C the smoothest \"\n      f\"predicts {float(mid @ c_smooth):.1f} C and the drawn one \"\n      f\"{float(mid @ c_int):.1f} C\")\ncheck(\"[B-overfit] six coefficients against three data points, so ERM \"\n      \"has infinitely many solutions\", null.shape[0] == DEGREE + 1 - len(IDX))\ncheck(\"[B-overfit] the training error does not separate them: both are \"\n      \"zero\", train_smooth < 1e-12 and train_int < 1e-12)\ncheck(\"[B-overfit] yet they differ by more than ten degrees between the \"\n      \"measured days\", gap > 10.0)\ncheck(\"[B-overfit] at a feature between the measured ones the solutions \"\n      \"take every real value, since a direction that leaves the three \"\n      \"labels untouched is nonzero there\", span > 1e-6)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-augment]** Replace each day by copies inside the accuracy of the thermometer and fit the same degree-five model to them."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "xa = np.repeat(xt, NR_COPIES) + RNG.uniform(-TOL, TOL, NR_COPIES * len(IDX))\nya = np.repeat(yt, NR_COPIES) + RNG.uniform(-TOL, TOL, NR_COPIES * len(IDX))\ncheck(\"[B-augment] every copy lies within the assumed accuracy of the \"\n      \"day it came from\",\n      np.all(np.abs(xa - np.repeat(xt, NR_COPIES)) <= TOL)\n      and np.all(np.abs(ya - np.repeat(yt, NR_COPIES)) <= TOL))\nAa = design(xa, LO, HI)\nc_aug = np.linalg.lstsq(Aa, ya, rcond=None)[0]\nfit_aug = design(GRID, LO, HI) @ c_aug\nexc_aug = max(fit_aug.max() - yt.max(), yt.min() - fit_aug.min())\nprint(f\"[B-augment] the augmented feature matrix is {Aa.shape[0]} by \"\n      f\"{Aa.shape[1]} of rank {np.linalg.matrix_rank(Aa)}\")\ndev_aug = np.max(np.abs(fit_aug - LINE))\ndev_int = np.max(np.abs(fit_int - LINE))\nprint(f\"[B-augment] the augmented fit leaves the range of the three \"\n      f\"labels by {exc_aug:.2f} C\")\nprint(f\"[B-augment] it stays within {dev_aug:.1f} C of the straight line \"\n      f\"through the three days, against {dev_int:.1f} C for the drawn \"\n      f\"solution\")\ncheck(\"[B-augment] the augmented feature matrix has full column rank, so \"\n      \"the augmented problem has a unique solution\",\n      np.linalg.matrix_rank(Aa) == DEGREE + 1)\ncheck(\"[B-augment] its excursion beyond the range of the three labels \"\n      \"stays below one degree\", exc_aug < 1.0)\ncheck(\"[B-augment] it departs from the straight line through the three \"\n      \"days by less than the drawn solution does\",\n      dev_aug < dev_int)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-compare]** The fits on the held-out day pairs whose feature lies between the smallest and the largest training feature, and the straight line through the three days."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# held out, and inside the feature range the three training days cover\ninside = (FEAT >= xt.min()) & (FEAT <= xt.max())\ninside[IDX] = False\nrest = np.flatnonzero(inside)\n\n\ndef val(coef, poly1d=False):\n    pred = (np.polyval(coef, FEAT[rest]) if poly1d\n            else design(FEAT[rest], LO, HI) @ coef)\n    return float(np.mean((pred - LAB[rest]) ** 2))\n\n\nv_int, v_aug, v_lin = val(c_int), val(c_aug), val(c_lin, True)\nv_smooth = val(c_smooth)\nprint(f\"[B-compare] held-out error on {len(rest)} day pairs inside the \"\n      f\"training range: drawn solution {v_int:.1f}, smoothest solution \"\n      f\"{v_smooth:.1f}, augmented fit {v_aug:.1f}, straight line \"\n      f\"{v_lin:.1f}\")\ncheck(\"[B-compare] the two ERM solutions differ widely on the held-out \"\n      \"days although both have zero training error\",\n      abs(v_int - v_smooth) > 20.0)\ncheck(\"[B-compare] the augmented fit beats the drawn solution\",\n      v_aug < v_int / 2.0)\ncheck(\"[B-compare] and lands close to the straight line through the \"\n      \"three days\", abs(v_aug - v_lin) < 0.3 * v_lin)\nwith open(OUT_DIR / \"dataaug_valerr.csv\", \"w\") as f:\n    f.write(\"fit,valerr\\n\")\n    f.write(f\"drawn,{v_int:.4f}\\n\")\n    f.write(f\"smoothest,{v_smooth:.4f}\\n\")\n    f.write(f\"augmented,{v_aug:.4f}\\n\")\n    f.write(f\"line,{v_lin:.4f}\\n\")\n\n# ---- the figure data\nwith open(OUT_DIR / \"dataaug_points.csv\", \"w\") as f:\n    f.write(\"x,y\\n\")\n    for a, b in zip(xt, yt):\n        f.write(f\"{a:.4f},{b:.4f}\\n\")\nwith open(OUT_DIR / \"dataaug_copies.csv\", \"w\") as f:\n    f.write(\"x,y\\n\")\n    for a, b in zip(xa[::25], ya[::25]):       # a readable subset\n        f.write(f\"{a:.4f},{b:.4f}\\n\")\none = slice(NR_COPIES, 2 * NR_COPIES)          # the middle measured day\nwith open(OUT_DIR / \"dataaug_zoom.csv\", \"w\") as f:\n    f.write(\"dx,dy\\n\")\n    for a, b in zip(xa[one][:28], ya[one][:28]):\n        f.write(f\"{a - xt[1]:.4f},{b - yt[1]:.4f}\\n\")\nwith open(OUT_DIR / \"dataaug_curves.csv\", \"w\") as f:\n    f.write(\"x,interpolant,augmented\\n\")\n    for a, b, c in zip(GRID, fit_int, fit_aug):\n        f.write(f\"{a:.4f},{b:.4f},{c:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fig]** The preview figure: the picture on the left, the held-out errors on the right."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (axP, axV) = plt.subplots(1, 2, figsize=(9.6, 3.8))\naxP.plot(GRID, fit_int, \"k-\", lw=1.6,\n         label=\"degree-5 fit on the three days\")\naxP.plot(GRID, fit_aug, \"k--\", lw=1.6,\n         label=\"degree-5 fit on the augmented set\")\naxP.plot(xa[::25], ya[::25], \"s\", mfc=\"none\", mec=\"0.45\", ms=4,\n         label=\"copy inside the sensor accuracy\")\naxP.plot(xt, yt, \"ko\", ms=7, label=\"measured day\")\naxP.set_xlabel(\"maximum temperature today in $^\\\\circ$C\")\naxP.set_ylabel(\"maximum temperature tomorrow in $^\\\\circ$C\")\naxP.set_title(\"Three measured days and two fitted hypotheses\",\n              fontsize=10)\naxP.set_ylim(yt.min() - 14, yt.max() + 14)\naxP.legend(frameon=False, fontsize=7, loc=\"lower right\")\nnames = [\"drawn ERM\\nsolution\", \"smoothest ERM\\nsolution\",\n         \"fit on the\\naugmented set\", \"straight line\\non three days\"]\naxV.bar(names, [v_int, v_smooth, v_aug, v_lin], color=\"0.55\",\n        edgecolor=\"black\")\naxV.set_xlabel(\"hypothesis compared\")\naxV.set_ylabel(\"average squared error in $(^\\\\circ$C$)^2$\")\naxV.set_title(f\"Error on the {len(rest)} day pairs held out\", fontsize=10)\naxV.set_yscale(\"log\")\naxV.tick_params(axis=\"x\", labelsize=7)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"dataaug.png\", dpi=110)\ncheck(\"[B-fig] the preview figure was written\",\n      (OUT_DIR / \"dataaug.png\").exists())\n\nprint()\nif FAILED:\n    print(f\"{len(FAILED)} check(s) FAILED: \" + \"; \".join(FAILED))\n    raise SystemExit(1)\nprint(\"all checks passed\")"
  }
 ]
}