{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "featuremtx.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# feature matrix \u2014 Python demo\n\nNumerical companion to the entry [feature matrix](https://dictionaryofml.org/terms/featuremtx.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 'featuremtx'. The GeoSphere Austria weather station Krems (station id 3805) records eight measurements per day: 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 featuremtx_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/featuremtx.py`](https://dictionaryofml.org/terms/featuremtx.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(), \"featuremtx.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"The feature matrix of one weather year, summarized by a handful of\nnumbers: its shape, its rank, its singular values and its condition number.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'featuremtx'.  The GeoSphere\nAustria weather station Krems (station id 3805) records eight\nmeasurements per day: 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\nfeaturemtx_weather.csv.  A value of -1 for precipitation marks a trace\nof rain too small to record and is set to 0.\n\nCollecting the eight measurements of each day as a row gives the feature\nmatrix of the year.  Every measurement is scaled to zero sample mean and\nunit sample variance first, because temperatures, precipitation and\npressure carry different units and the singular values of an unscaled\nmatrix would report the choice of units rather than the data.\n\nThe demo checks the claims the entry makes.  (1) The shape is the pair\n(number of data points, number of features), and the whole year is 366\ntimes 8 numbers while the shape is two of them.  (2) The singular values\nare the square roots of the eigenvalues of X^T X, and X^T X divided by\nthe number of data points is the sample covariance matrix of the scaled\nmeasurements, so the singular values carry the spectrum that PCA reads.\n(3) The sum of the squared singular values is the squared Frobenius norm\nof the feature matrix, which for scaled measurements is the number of\ndata points times the number of features.  (4) The condition number is\nthe ratio of the largest to the smallest singular value.  (5) The\nsmallest singular value is far below the largest because the mean\ntemperature of a day is close to the midpoint of its minimum and its\nmaximum, a near-dependency among three columns; dropping the mean\ntemperature raises the smallest singular value by more than an order of\nmagnitude.\n\nDeterministic: the data are a fixed archive year and the decomposition\nis a singular value decomposition.  Self-contained: numpy + matplotlib\nonly (stdlib urllib for the download).\n\nBlocks\n------\n[B-data]      Download the 366 days with eight measurements each and\n              scale every measurement to zero sample mean and unit\n              sample variance.\n[B-shape]     The shape of the feature matrix, and how many numbers it\n              takes to state it.\n[B-spectrum]  The singular values, their link to the eigenvalues of\n              X^T X and to the sample covariance matrix, and the share\n              of the squared Frobenius norm that the leading ones carry.\n[B-cond]      The rank and the condition number, and the near-dependency\n              among the three temperature columns that makes the\n              smallest singular value small.\n[B-fig]       The preview figure: the spectrum as a stem plot and the\n              share of the squared Frobenius norm that k singular values\n              carry.\n\nOutputs\n-------\nfeaturemtx_weather.csv   : date and the eight measurements, 366 days\nfeaturemtx_spectrum.csv  : index, singular value, cumulative share\nfeaturemtx.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\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}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** Download the 366 days with eight measurements each and scale every measurement to zero sample mean and unit sample variance."
  },
  {
   "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=180) as resp:\n    payload = json.load(resp)\nparams = payload[\"features\"][0][\"properties\"][\"parameters\"]\nstamps = [t[:10] for t in payload[\"timestamps\"]]\nRAW = np.stack([np.array(params[p][\"data\"], dtype=float) for p in PARAMS], 1)\nRAW[:, 3] = np.maximum(RAW[:, 3], 0.0)       # -1 marks a trace of rain\nwith open(OUT_DIR / \"featuremtx_weather.csv\", \"w\") as f:\n    f.write(\"date,\" + \",\".join(PARAMS) + \"\\n\")\n    for day, row in zip(stamps, RAW):\n        f.write(day + \",\" + \",\".join(f\"{v:g}\" for v in row) + \"\\n\")\n\nX = (RAW - RAW.mean(0)) / RAW.std(0)         # zero mean, unit variance\nprint(\"[B-data] the feature matrix of one weather year\")\ncheck(\"[B-data] every measurement has zero sample mean and unit variance\",\n      np.allclose(X.mean(0), 0.0, atol=1e-12)\n      and np.allclose(X.std(0), 1.0, atol=1e-12))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-shape]** The shape of the feature matrix, and how many numbers it takes to state it."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "m, d = X.shape\nprint(f\"[B-shape] shape ({m}, {d}): {m} data points, {d} features\")\nprint(f\"[B-shape] {m * d} numbers in the matrix, 2 in its shape\")\ncheck(\"[B-shape] 366 days with eight measurements each\",\n      (m, d) == (366, 8))\ncheck(\"[B-shape] more data points than features\", m > d)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-spectrum]** The singular values, their link to the eigenvalues of X^T X and to the sample covariance matrix, and the share of the squared Frobenius norm that the leading ones carry."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "sv = np.linalg.svd(X, compute_uv=False)\neigs = np.linalg.eigvalsh(X.T @ X)[::-1]\ncum = np.cumsum(sv ** 2) / np.sum(sv ** 2)\nprint(\"[B-spectrum] singular values: \"\n      + \", \".join(f\"{s:.2f}\" for s in sv))\nprint(f\"[B-spectrum] two of them carry {100 * cum[1]:.0f}% of the \"\n      f\"squared Frobenius norm, four carry {100 * cum[3]:.0f}%\")\ncheck(\"[B-spectrum] the squared singular values are the eigenvalues \"\n      \"of X^T X\", np.allclose(sv ** 2, eigs, rtol=1e-9, atol=1e-9))\nsamplecov = (X.T @ X) / m\ncheck(\"[B-spectrum] X^T X divided by the number of data points is the \"\n      \"sample covariance matrix of the scaled measurements\",\n      np.allclose(samplecov, np.cov(X, rowvar=False, bias=True),\n                  atol=1e-12))\ncheck(\"[B-spectrum] the squared singular values sum to the squared \"\n      \"Frobenius norm, which is m times d here\",\n      abs(np.sum(sv ** 2) - m * d) < 1e-6)\nwith open(OUT_DIR / \"featuremtx_spectrum.csv\", \"w\") as f:\n    f.write(\"j,sigma,share\\n\")\n    for j, (s, c) in enumerate(zip(sv, cum), start=1):\n        f.write(f\"{j},{s:.6f},{c:.6f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-cond]** The rank and the condition number, and the near-dependency among the three temperature columns that makes the smallest singular value small."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rank = int(np.linalg.matrix_rank(X))\ncond = sv[0] / sv[-1]\nprint(f\"[B-cond] rank {rank}, condition number {cond:.1f}\")\ncheck(\"[B-cond] the rank is the number of features, so no measurement \"\n      \"is an exact combination of the others\", rank == d)\ncheck(\"[B-cond] the condition number is the ratio of the largest to \"\n      \"the smallest singular value\",\n      abs(cond - np.linalg.cond(X)) < 1e-6)\n# the mean temperature is close to the midpoint of the extremes\nmid = 0.5 * (RAW[:, 0] + RAW[:, 1])\ngap = np.abs(RAW[:, 2] - mid)\nprint(f\"[B-cond] mean temperature differs from the midpoint of the \"\n      f\"extremes by {gap.mean():.2f} degrees on average\")\nkeep = [i for i in range(d) if i != 2]            # drop mean temperature\nsv_keep = np.linalg.svd(X[:, keep], compute_uv=False)\ncond_keep = sv_keep[0] / sv_keep[-1]\nprint(f\"[B-cond] dropping the mean temperature: smallest singular value \"\n      f\"{sv[-1]:.2f} -> {sv_keep[-1]:.2f}, condition number \"\n      f\"{cond:.1f} -> {cond_keep:.1f}\")\ncheck(\"[B-cond] dropping the mean temperature raises the smallest \"\n      \"singular value by more than an order of magnitude\",\n      sv_keep[-1] > 10 * sv[-1])\ncheck(\"[B-cond] and lowers the condition number by more than an order \"\n      \"of magnitude\", cond_keep < cond / 10)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fig]** The preview figure: the spectrum as a stem plot and the share of the squared Frobenius norm that k singular values carry."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (axL, axR) = plt.subplots(1, 2, figsize=(9.4, 3.5))\naxL.stem(np.arange(1, d + 1), sv, linefmt=\"k-\", markerfmt=\"ko\",\n         basefmt=\" \")\naxL.set_xlabel(\"index $j$ of the singular value\")\naxL.set_ylabel(\"singular value $\\\\sigma_j$\")\naxL.set_title(\"Spectrum of the feature matrix\\n(366 days, 8 scaled \"\n              \"measurements)\", fontsize=10)\naxL.set_xticks(np.arange(1, d + 1))\naxR.plot(np.arange(1, d + 1), 100 * cum, \"k-o\")\naxR.axhline(90, color=\"gray\", linestyle=\"--\", label=\"90 percent\")\naxR.set_xlabel(\"number $k$ of leading singular values kept\")\naxR.set_ylabel(\"percent of the squared Frobenius norm\")\naxR.set_title(\"Share the leading singular values carry\", fontsize=10)\naxR.set_xticks(np.arange(1, d + 1))\naxR.set_ylim(0, 105)\naxR.legend(frameon=False)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"featuremtx.png\", dpi=110)\ncheck(\"[B-fig] the preview figure was written\",\n      (OUT_DIR / \"featuremtx.png\").exists())\n\nprint()\nif FAILED:\n    print(f\"{len(FAILED)} check(s) FAILED: \" + \"; \".join(FAILED))\n    raise SystemExit(1)\nprint(\"all checks passed\")"
  }
 ]
}