{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "unsupervisedlearning.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# unsupervised learning \u2014 Python demo\n\nNumerical companion to the entry [unsupervised learning](https://dictionaryofml.org/terms/unsupervisedlearning.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 'unsupervisedlearning'. 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 unsupervisedlearning_weather.csv. Each day is a data point with eight features and no label: nothing in the record says what should be predicted for a day. 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/unsupervisedlearning.py`](https://dictionaryofml.org/terms/unsupervisedlearning.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(), \"unsupervisedlearning.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Unsupervised learning on a year of days at Krems: the same days\nwithout any label, grouped into clusters, compressed to two features,\nand fitted with a density -- and no criterion that says which of these\nis right.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'unsupervisedlearning'.  The\nGeoSphere Austria 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\nunsupervisedlearning_weather.csv.  Each day is a data point with eight\nfeatures and no label: nothing in the record says what should be\npredicted for a day.  A value of -1 for precipitation marks a trace of\nrain too small to record and is set to 0.\n\nThe demo checks the entry's claims about the three tasks the entry\nnames and about their evaluation: (1) clustering the days by k-means\nwith k = 2 separates the cold from the warm half of the year, although\nno month, season or other label ever entered the computation; (2)\nprincipal component analysis compresses the eight features to two with\na reconstruction error equal to the sum of the dropped eigenvalues of\nthe sample covariance matrix; (3) a Gaussian density fitted to the days\nassigns a higher average log-density to held-out days than a density\nthat ignores the correlation between the measurements; and (4) none of\nthese has a direct measure of success: the smallest clustering error\nfound keeps falling as clusters are added, so it cannot say how many\nclusters the data have, while in supervised learning the validation\nerror does say when a hypothesis is worse.\n\nDeterministic: the data are a fixed archive year, the initial centroids\nof the two-cluster run are the coldest and the warmest day, the\nrestarts of the last block draw from a fixed seed, and PCA is an\neigenvalue decomposition.  Self-contained: numpy + matplotlib only (stdlib urllib\nfor the download).\n\nBlocks\n------\n[B-data]     Download the 366 days with eight measurements each, scale\n             every feature to zero sample mean and unit sample\n             variance, and check that the data carry no label.\n[B-cluster]  k-means with k = 2 on the eight scaled features: check the\n             clustering error never increases, that the iteration\n             reaches a fixed point, and that the two clusters agree\n             with the cold and warm half of the year on 90 percent of the days -- a comparison\n             made only after the clustering, never during it.\n[B-dimred]   Principal component analysis to two features: check the\n             reconstruction error equals the sum of the dropped\n             eigenvalues and that the two components carry more than\n             half of the total variance.\n[B-density]  A Gaussian fitted to the two temperature features: check\n             its average log-density on held-out days exceeds that of a\n             Gaussian with the same means but no correlation.\n[B-nocriterion] No direct measure of success: the smallest clustering\n             error found decreases with every added cluster (ten restarts\n             per number of clusters, the best kept), so it cannot choose\n             the number of clusters, whereas the validation error of\n             supervised learning does grow when a hypothesis is worse.\n[B-plot]     Preview: the days in the two learned features marked by\n             cluster, and the clustering error against the number of\n             clusters.\n\nOutputs\n-------\nunsupervisedlearning_weather.csv : date and the eight measurements\nunsupervisedlearning_cluster1.csv, _cluster2.csv : z1, z2 per cluster\nunsupervisedlearning_centroids.csv : z1, z2 of the two centroids\nunsupervisedlearning_error.csv   : nrcluster, error\nunsupervisedlearning.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, scale every feature to zero sample mean and unit sample variance, and check that the data carry no label."
  },
  {
   "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 / \"unsupervisedlearning_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))\ncheck(\"[B-data] 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]))\ncheck(\"[B-data] the data carry features only, no label column\",\n      M.shape[1] == len(PARAMS))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-cluster]** k-means with k = 2 on the eight scaled features: check the clustering error never increases, that the iteration reaches a fixed point, and that the two clusters agree with the cold and warm half of the year on 90 percent of the days -- a comparison made only after the clustering, never during it."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def assign(X, centroids):\n    dist = ((X[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)\n    return dist.argmin(axis=1)\n\n\ndef clustering_error(X, centroids, labels):\n    return float(((X - centroids[labels]) ** 2).sum())\n\n\ndef lloyd(X, centroids):\n    labels = assign(X, centroids)\n    errors = [clustering_error(X, centroids, labels)]\n    for _ in range(100):\n        centroids = np.stack([X[labels == c].mean(axis=0)\n                              if (labels == c).any() else centroids[c]\n                              for c in range(len(centroids))])\n        new = assign(X, centroids)\n        errors.append(clustering_error(X, centroids, new))\n        if np.array_equal(new, labels):\n            break\n        labels = new\n    return centroids, labels, np.array(errors)\n\n\nmean_temp = X[:, 2]\nstart = np.stack([X[mean_temp.argmin()], X[mean_temp.argmax()]])\ncentroids, labels, errors = lloyd(X, start)\ncheck(\"[B-cluster] the clustering error never increases\",\n      bool(np.all(np.diff(errors) <= 1e-9)))\ncheck(\"[B-cluster] the iteration reaches a fixed point\",\n      np.array_equal(assign(X, centroids), labels))\nmonth = np.array([int(s[5:7]) for s in stamps])\ncold_half = (month <= 3) | (month >= 11)       # used only for checking\nwarm_half = (month >= 5) & (month <= 9)\nknown = cold_half | warm_half\nagree = max(np.mean(labels[known] == cold_half[known].astype(int)),\n            np.mean(labels[known] == warm_half[known].astype(int)))\nprint(f\"  k-means: {len(errors) - 1} iterations, clustering error \"\n      f\"{errors[-1]:.0f}; the two clusters agree with the cold and warm \"\n      f\"half of the year on {100 * agree:.1f}% of those days\")\ncheck(\"[B-cluster] the two clusters agree with the cold and warm half of \"\n      \"the year on more than 85 percent of the days\", agree > 0.85)\ncheck(\"[B-cluster] no month or season entered the clustering\",\n      start.shape == (2, d))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-dimred]** Principal component analysis to two features: check the reconstruction error equals the sum of the dropped eigenvalues and that the two components carry more than half of the total variance."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "Q = X.T @ X / m\neigvals, eigvecs = np.linalg.eigh(Q)\norder = np.argsort(eigvals)[::-1]\neigvals, eigvecs = eigvals[order], eigvecs[:, order]\nW = eigvecs[:, :2].T\nZ = X @ W.T\nrecon = float(((X - Z @ W) ** 2).sum())\nprint(f\"  PCA: eigenvalues 1-3 {eigvals[0]:.2f}, {eigvals[1]:.2f}, \"\n      f\"{eigvals[2]:.2f} of total {eigvals.sum():.0f}; reconstruction \"\n      f\"error {recon:.0f}\")\ncheck(\"[B-dimred] the reconstruction error equals m times the sum of the \"\n      \"dropped eigenvalues\", np.isclose(recon, m * eigvals[2:].sum()))\ncheck(\"[B-dimred] the two components carry more than half of the total \"\n      \"variance\", eigvals[:2].sum() > 0.5 * eigvals.sum())"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-density]** A Gaussian fitted to the two temperature features: check its average log-density on held-out days exceeds that of a Gaussian with the same means but no correlation."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "temps = M[:, :2]\nfit_days, held_days = temps[:244], temps[244:]\nmu = fit_days.mean(axis=0)\nC = np.cov(fit_days, rowvar=False)\nC_diag = np.diag(np.diag(C))                   # same means, no correlation\n\n\ndef mean_log_density(U, mean, cov):\n    diff = U - mean\n    inv = np.linalg.inv(cov)\n    quad = np.einsum(\"ij,jk,ik->i\", diff, inv, diff)\n    return float(np.mean(-0.5 * quad - 0.5 * np.log(np.linalg.det(cov))\n                         - np.log(2.0 * np.pi)))\n\n\nll_full = mean_log_density(held_days, mu, C)\nll_diag = mean_log_density(held_days, mu, C_diag)\nprint(f\"  density: average log-density on held-out days {ll_full:.3f} with \"\n      f\"the fitted covariance against {ll_diag:.3f} without correlation\")\ncheck(\"[B-density] the fitted Gaussian beats the one that ignores the \"\n      \"correlation between the measurements\", ll_full > ll_diag)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-nocriterion]** No direct measure of success: the smallest clustering error found decreases with every added cluster (ten restarts per number of clusters, the best kept), so it cannot choose the number of clusters, whereas the validation error of supervised learning does grow when a hypothesis is worse."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# The k-means iteration finds a local minimum, so the error of a single run\n# need not fall when a cluster is added; ten restarts per number of\n# clusters, the best kept, approximate the smallest error attainable.\nrng = np.random.default_rng(0)\ncurve = []\nfor k in range(1, 7):\n    best = np.inf\n    for _ in range(10):\n        init = X[rng.choice(m, size=k, replace=False)]\n        _, _, err_k = lloyd(X, init)\n        best = min(best, err_k[-1])\n    curve.append((k, best))\nerrs = np.array([e for _, e in curve])\nprint(\"  clustering error by number of clusters: \"\n      + \", \".join(f\"k={k}: {e:.0f}\" for k, e in curve))\ncheck(\"[B-nocriterion] the smallest clustering error found decreases with \"\n      \"every added cluster, so it cannot choose their number\",\n      bool(np.all(np.diff(errs) < 0)))\n\n# ---- CSVs for the entry's figure\nnp.savetxt(OUT_DIR / \"unsupervisedlearning_cluster1.csv\", Z[labels == 0],\n           delimiter=\",\", header=\"z1,z2\", comments=\"\", fmt=\"%.3f\")\nnp.savetxt(OUT_DIR / \"unsupervisedlearning_cluster2.csv\", Z[labels == 1],\n           delimiter=\",\", header=\"z1,z2\", comments=\"\", fmt=\"%.3f\")\nnp.savetxt(OUT_DIR / \"unsupervisedlearning_centroids.csv\",\n           np.stack([centroids[c] @ W.T for c in (0, 1)]), delimiter=\",\",\n           header=\"z1,z2\", comments=\"\", fmt=\"%.3f\")\nwith open(OUT_DIR / \"unsupervisedlearning_error.csv\", \"w\") as f:\n    f.write(\"nrcluster,error\\n\")\n    for k, e in curve:\n        f.write(f\"{k},{e:.1f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Preview: the days in the two learned features marked by cluster, and the clustering error against the number of clusters."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2))\nax1.scatter(Z[labels == 0, 0], Z[labels == 0, 1], marker=\"s\",\n            facecolors=\"none\", edgecolors=\"tab:blue\", s=20,\n            label=\"first cluster\")\nax1.scatter(Z[labels == 1, 0], Z[labels == 1, 1], marker=\"o\", color=\"tab:red\",\n            s=16, label=\"second cluster\")\ncent = np.stack([centroids[c] @ W.T for c in (0, 1)])\nax1.scatter(cent[:, 0], cent[:, 1], marker=\"X\", color=\"black\", s=90,\n            label=\"cluster centroid\")\nax1.set_xlabel(\"first learned feature $z_1$\")\nax1.set_ylabel(\"second learned feature $z_2$\")\nax1.set_title(\"366 unlabeled days, grouped without any label\")\nax1.legend(frameon=False, fontsize=8)\nax2.plot([k for k, _ in curve], errs, \"ko-\")\nax2.set_xlabel(\"number of clusters\")\nax2.set_ylabel(\"clustering error\")\nax2.set_title(\"The error keeps falling: it cannot choose the number\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"unsupervisedlearning.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)"
  }
 ]
}