{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "featlearn.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# feature learning \u2014 Python demo\n\nNumerical companion to the entry [feature learning](https://dictionaryofml.org/terms/featlearn.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 'featlearn' (feature learning). The GeoSphere Austria weather station Krems (station id 3805) records eight daily measurements: minimum, maximum and mean air temperature, precipitation, sunshine duration, relative humidity, air pressure and wind speed. This script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to featlearn_weather.csv. A value of -1 for precipitation marks a trace of rain too small to record and is set to 0.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/featlearn.py`](https://dictionaryofml.org/terms/featlearn.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(), \"featlearn.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Feature learning on a year of weather at Krems: forty raw measurements\nof the previous five days become two learned features, and a linear\nmodel on the two features predicts the next day's maximum temperature.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'featlearn' (feature learning).\nThe GeoSphere Austria weather station Krems (station id 3805) records\neight daily measurements: minimum, maximum and mean air temperature,\nprecipitation, sunshine duration, relative humidity, air pressure and\nwind speed.  This script downloads the records for 2024 from the\nGeoSphere data hub (dataset klima-v2-1d) and writes them to\nfeatlearn_weather.csv.  A value of -1 for precipitation marks a trace\nof rain too small to record and is set to 0.\n\nEach day from January 6 to December 31 is a data point.  Its label is\nthe maximum temperature of that day; its raw features are the eight\nmeasurements of each of the five previous days, a list of forty\nnumbers.  A day cannot be drawn as a point with forty coordinates, so\nthe demo learns a feature transformation that delivers two new\nfeatures: principal component analysis (PCA) on the forty raw features,\neach scaled to zero sample mean and unit sample variance because\ntemperatures, precipitation and pressure carry different units.  The\ndemo checks the entry's claims: the PCA transformation has the minimum\nlinear reconstruction error, equal to the sum of the dropped eigenvalues\nof the sample covariance matrix; the two learned features are the\ncoordinates of a scatterplot in which warm and cold days separate,\nmost of each class on its own side of the first feature's zero; and\nlinear regression on the two learned features predicts the label on\nheld-out days far better than the sample mean of the training labels.\n\nDeterministic: the data are a fixed archive year, PCA is an eigenvalue\ndecomposition, and the random linear transformations used for\ncomparison are drawn with a fixed seed.  Self-contained: numpy +\nmatplotlib only (stdlib urllib for the download).\n\nBlocks\n------\n[B-fetch] Download the eight daily measurements at Krems for 2024 and\n          write them to featlearn_weather.csv; check the count and one\n          pinned value against the archive.\n[B-dataset] Build the 361 data points: forty raw features (five previous\n          days times eight measurements) and the label (maximum\n          temperature of the day); scale every raw feature to zero\n          sample mean and unit sample variance.\n[B-pca] Learn the feature transformation: the two eigenvectors of the\n          sample covariance matrix with the largest eigenvalues; check\n          the reconstruction-error identity, that no random linear\n          transformation reconstructs better, and that warm and cold\n          days separate along the first learned feature.\n[B-linreg] Fit linear regression to the two learned features on the\n          days of January to August and validate on September to\n          December; compare with the forty raw features and with the\n          sample mean of the training labels.\n[B-plot] Preview: the scatterplot of the days in the two learned\n          features, warm and cold days marked differently, and the\n          largest ten eigenvalues of the sample covariance matrix (the\n          first is 16.4 of a total of 40, i.e., 41 percent).\n\nOutputs\n-------\nfeatlearn_weather.csv   : date and the eight measurements, 366 days of 2024\nfeatlearn_warm.csv      : z1, z2 of the days whose label is above the median\nfeatlearn_cold.csv      : z1, z2 of the days whose label is at or below it\nfeatlearn_eigvals.csv   : index, eigenvalue of the sample covariance matrix\nfeatlearn_valerr.csv    : model, validation error (mean squared error)\nfeatlearn.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-fetch]** Download the eight daily measurements at Krems for 2024 and write them to featlearn_weather.csv; check the count and one pinned value against the archive."
  },
  {
   "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 / \"featlearn_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\")\ncheck(\"[B-fetch] 366 days with eight measurements downloaded for 2024\",\n      M.shape == (366, 8))\ncheck(\"[B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4 deg)\",\n      stamps[31] == \"2024-02-01\" and np.allclose(M[31, :2], [-3.8, 10.4]))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-dataset]** Build the 361 data points: forty raw features (five previous days times eight measurements) and the label (maximum temperature of the day); scale every raw feature to zero sample mean and unit sample variance."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "LAG = 5\ndays = np.arange(LAG, len(M))                       # Jan 6 .. Dec 31\nX_raw = np.stack([M[t - LAG:t].ravel() for t in days])   # 361 x 40\ny = M[days, 1]                                      # tlmax of the day\nm, d = X_raw.shape\nmu, sigma = X_raw.mean(axis=0), X_raw.std(axis=0)\nX = (X_raw - mu) / sigma                            # zero mean, unit variance\ncheck(\"[B-dataset] 361 data points with forty raw features each\",\n      (m, d) == (361, 40))\ncheck(\"[B-dataset] every scaled feature has unit sample variance\",\n      np.allclose(X.var(axis=0), 1.0))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-pca]** Learn the feature transformation: the two eigenvectors of the sample covariance matrix with the largest eigenvalues; check the reconstruction-error identity, that no random linear transformation reconstructs better, and that warm and cold days separate along the first learned feature."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "Q = X.T @ X / m                                     # sample covariance matrix\neigvals, eigvecs = np.linalg.eigh(Q)\norder = np.argsort(eigvals)[::-1]\neigvals, eigvecs = eigvals[order], eigvecs[:, order]\nW = eigvecs[:, :2].T                                # 2 x 40\nZ = X @ W.T                                         # z = W x, 361 x 2\n\n\ndef reconstruction_error(W):\n    \"\"\"Minimum over R of sum_r ||x^(r) - R W x^(r)||^2 (least squares).\"\"\"\n    Zc = X @ W.T\n    R = np.linalg.lstsq(Zc, X, rcond=None)[0].T     # d x 2\n    return float(((X - Zc @ R.T) ** 2).sum())\n\n\nerr_pca = reconstruction_error(W)\nrng = np.random.default_rng(0)\nerr_random = [reconstruction_error(rng.standard_normal((2, d)))\n              for _ in range(20)]\nprint(f\"  eigenvalues 1-3: {eigvals[0]:.2f}, {eigvals[1]:.2f}, \"\n      f\"{eigvals[2]:.2f} of total {eigvals.sum():.0f}\")\nprint(f\"  reconstruction error PCA {err_pca:.0f}, best of 20 random \"\n      f\"transformations {min(err_random):.0f}\")\ncheck(\"[B-pca] the reconstruction error equals m times the sum of the \"\n      \"dropped eigenvalues\", np.isclose(err_pca, m * eigvals[2:].sum()))\ncheck(\"[B-pca] no random linear transformation reconstructs better\",\n      err_pca < min(err_random))\nwarm = y > np.median(y)\nfrac_warm = float((Z[warm, 0] > 0).mean())      # warm days right of z1 = 0\nfrac_cold = float((Z[~warm, 0] <= 0).mean())    # cold days left of it\nprint(f\"  first learned feature positive for {100 * frac_warm:.0f}% of the \"\n      f\"warm days, nonpositive for {100 * frac_cold:.0f}% of the cold days\")\ncheck(\"[B-pca] warm and cold days separate along the first learned feature \"\n      \"(at least 85% of each class on its side of z1 = 0)\",\n      min(frac_warm, frac_cold) >= 0.85)\nnp.savetxt(OUT_DIR / \"featlearn_warm.csv\", Z[warm], delimiter=\",\",\n           header=\"z1,z2\", comments=\"\", fmt=\"%.3f\")\nnp.savetxt(OUT_DIR / \"featlearn_cold.csv\", Z[~warm], delimiter=\",\",\n           header=\"z1,z2\", comments=\"\", fmt=\"%.3f\")\nnp.savetxt(OUT_DIR / \"featlearn_eigvals.csv\",\n           np.stack([np.arange(1, d + 1), eigvals], 1), delimiter=\",\",\n           header=\"index,eigenvalue\", comments=\"\", fmt=[\"%d\", \"%.4f\"])"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-linreg]** Fit linear regression to the two learned features on the days of January to August and validate on September to December; compare with the forty raw features and with the sample mean of the training labels."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "train = days < 244                                  # Jan 6 .. Aug 31\nval = ~train                                        # Sep 1 .. Dec 31\n\n\ndef linreg_valerr(F):\n    \"\"\"Least squares with intercept on the training days; validation MSE.\"\"\"\n    A = np.hstack([F, np.ones((m, 1))])\n    w = np.linalg.lstsq(A[train], y[train], rcond=None)[0]\n    return float(((A[val] @ w - y[val]) ** 2).mean())\n\n\nvalerr = {\"two learned features\": linreg_valerr(Z),\n          \"forty raw features\": linreg_valerr(X),\n          \"sample mean of training labels\":\n              float(((y[train].mean() - y[val]) ** 2).mean())}\nwith open(OUT_DIR / \"featlearn_valerr.csv\", \"w\") as f:\n    f.write(\"model,valerr\\n\")\n    for name, e in valerr.items():\n        f.write(f\"{name},{e:.2f}\\n\")\n        print(f\"  validation error {name}: {e:.2f}\")\ncheck(\"[B-linreg] two learned features beat the sample mean of the labels\",\n      valerr[\"two learned features\"] < valerr[\"sample mean of training labels\"] / 2)\ncheck(\"[B-linreg] forty raw features are not more than twice as good\",\n      valerr[\"forty raw features\"] > valerr[\"two learned features\"] / 2)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Preview: the scatterplot of the days in the two learned features, warm and cold days marked differently, and the largest ten eigenvalues of the sample covariance matrix (the first is 16.4 of a total of 40, i.e., 41 percent)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))\nax1.scatter(Z[~warm, 0], Z[~warm, 1], marker=\"s\", facecolors=\"none\",\n            edgecolors=\"tab:blue\", s=22, label=\"cold day (label at or below median)\")\nax1.scatter(Z[warm, 0], Z[warm, 1], marker=\"o\", color=\"tab:red\", s=18,\n            label=\"warm day (label above median)\")\nax1.set_xlabel(\"learned feature $z_1$\")\nax1.set_ylabel(\"learned feature $z_2$\")\nax1.set_title(\"361 days of 2024 at Krems in the two learned features\")\nax1.legend(frameon=False, fontsize=8)\nax2.bar(np.arange(1, 11), eigvals[:10], color=\"0.5\", edgecolor=\"black\")\nax2.set_xlabel(\"index of the eigenvalue\")\nax2.set_ylabel(\"eigenvalue of the sample covariance matrix\")\nax2.set_title(\"Largest ten eigenvalues: the first carries 41% of the total\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"featlearn.png\", dpi=110)\n\nprint()\nfailed = [n for n, ok in report if not ok]\nprint(f\"{len(report) - len(failed)}/{len(report)} checks passed\"\n      + (f\"; FAILED: {failed}\" if failed else \"\"))"
  }
 ]
}