{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "kmeans.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# k-means \u2014 Python demo\n\nNumerical companion to the entry [k-means](https://dictionaryofml.org/terms/kmeans.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 'kmeans' ($k$-means). The GeoSphere Austria weather station Krems (station id 3805) records the minimum and maximum air temperature of each day; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to kmeans_weather.csv. Each day is a data point with feature vector (tmin, tmax), and $k$-means with $k = 2$ partitions the 366 days into a cold-season and a warm-season cluster.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/kmeans.py`](https://dictionaryofml.org/terms/kmeans.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(), \"kmeans.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Lloyd's algorithm on a year of days at Krems: the clustering error\nnever increases, and the iteration stops at a fixed point.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'kmeans' ($k$-means).  The\nGeoSphere Austria weather station Krems (station id 3805) records the\nminimum and maximum air temperature of each day; this script downloads\nthe records for 2024 from the GeoSphere data hub (dataset klima-v2-1d)\nand writes them to kmeans_weather.csv.  Each day is a data point with\nfeature vector (tmin, tmax), and $k$-means with $k = 2$ partitions the\n366 days into a cold-season and a warm-season cluster.\n\nThe demo checks the entry's central claims about Lloyd's algorithm:\neach iteration -- assign every data point to its nearest cluster\ncentroid, then recompute each centroid as the mean of its assigned\ndata points -- never increases the clustering error; after finitely\nmany iterations nothing changes anymore, i.e., the iteration reaches a\nfixed point, where each centroid coincides with the mean of the data\npoints assigned to it.\n\nDeterministic: the initial centroids are the coldest and the warmest\nday of the year (no randomness).  Self-contained: numpy + matplotlib\nonly (stdlib urllib for the download).\n\nBlocks\n------\n[B-fetch] Download the 366 daily temperature pairs at Krems for 2024\n          and write them to kmeans_weather.csv; check the count and\n          one pinned value against the archive.\n[B-lloyd] Run Lloyd's algorithm with k = 2 from the coldest/warmest\n          day; check that the clustering error never increases and\n          that assignments stop changing after finitely many\n          iterations.\n[B-fixedpoint] Verify the fixed-point property of the result: each\n          final cluster centroid equals the mean of the data points\n          assigned to it, so one more iteration changes nothing.\n[B-image] Image compression and image segmentation on a subsampled\n          photo of the Oetscher massif (assets/oetscher.jpg): k-means\n          on the pixel colors with k = 4 replaces each pixel's color by\n          the nearest palette color (compression factor ~12), and with\n          k = 2 partitions the pixels into a sky-and-mountain region\n          and a vegetation region.\n\nOutputs\n-------\nkmeans_weather.csv             : date, tmin, tmax for the 366 days of 2024\nkmeans_cluster1.csv            : tmin, tmax of the days in the cold cluster\nkmeans_cluster2.csv            : tmin, tmax of the days in the warm cluster\nkmeans_centroids.csv           : tmin, tmax of the two final cluster centroids\nkmeans_error.csv               : iter, error -- clustering error per iteration\nkmeans_oetscher_original.png   : the subsampled photo\nkmeans_oetscher_compressed.png : the photo quantized to 4 palette colors\nkmeans_oetscher_mask.png       : sky-and-mountain/vegetation mask (k = 2)\nkmeans.png                     : preview (checking only) -- the two clusters\n                                 with their centroids, and the monotone\n                                 clustering error\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 366 daily temperature pairs at Krems for 2024 and write them to kmeans_weather.csv; check the count and one pinned value against the archive."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "URL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       \"klima-v2-1d?parameters=tlmin,tlmax&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\"]]\nX = np.stack([np.array(params[\"tlmin\"][\"data\"], dtype=float),\n              np.array(params[\"tlmax\"][\"data\"], dtype=float)], 1)\nwith open(OUT_DIR / \"kmeans_weather.csv\", \"w\") as f:\n    f.write(\"date,tmin,tmax\\n\")\n    for day, (lo, hi) in zip(stamps, X):\n        f.write(f\"{day},{lo},{hi}\\n\")\ncheck(\"[B-fetch] 366 daily temperature pairs downloaded for 2024\",\n      len(X) == 366)\ncheck(\"[B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4)\",\n      stamps[31] == \"2024-02-01\" and np.allclose(X[31], [-3.8, 10.4]))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-lloyd]** Run Lloyd's algorithm with k = 2 from the coldest/warmest day; check that the clustering error never increases and that assignments stop changing after finitely many iterations."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def assign(X, centroids):\n    \"\"\"Index of the nearest cluster centroid for every data point.\"\"\"\n    d = ((X[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)\n    return d.argmin(axis=1)\n\n\ndef clustering_error(X, centroids, labels):\n    return float(((X - centroids[labels]) ** 2).sum())\n\n\nmean_temp = X.mean(axis=1)\ncentroids = np.stack([X[mean_temp.argmin()], X[mean_temp.argmax()]])\nlabels = assign(X, centroids)\nerrors = [clustering_error(X, centroids, labels)]\niterations = 0\nwhile True:\n    iterations += 1\n    centroids = np.stack([X[labels == c].mean(axis=0) for c in (0, 1)])\n    new_labels = assign(X, centroids)\n    errors.append(clustering_error(X, centroids, new_labels))\n    if np.array_equal(new_labels, labels):\n        break\n    labels = new_labels\nerrors = np.array(errors)\nprint(f\"  fixed point reached after {iterations} iterations, \"\n      f\"clustering error {errors[-1]:.0f}\")\ncheck(\"[B-lloyd] the clustering error never increases\",\n      bool(np.all(np.diff(errors) <= 1e-9)))\ncheck(\"[B-lloyd] assignments stop changing after finitely many iterations\",\n      iterations < 50)\ncheck(\"[B-lloyd] both clusters are nonempty\",\n      0 < int(labels.sum()) < len(X))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fixedpoint]** Verify the fixed-point property of the result: each final cluster centroid equals the mean of the data points assigned to it, so one more iteration changes nothing."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "recomputed = np.stack([X[labels == c].mean(axis=0) for c in (0, 1)])\ncheck(\"[B-fixedpoint] each centroid equals the mean of its assigned \"\n      \"data points\", np.allclose(recomputed, centroids))\ncheck(\"[B-fixedpoint] one more iteration changes no assignment\",\n      np.array_equal(assign(X, recomputed), labels))\n\nheader = \"tmin,tmax\"\nnp.savetxt(OUT_DIR / \"kmeans_cluster1.csv\", X[labels == 0],\n           delimiter=\",\", header=header, comments=\"\", fmt=\"%.1f\")\nnp.savetxt(OUT_DIR / \"kmeans_cluster2.csv\", X[labels == 1],\n           delimiter=\",\", header=header, comments=\"\", fmt=\"%.1f\")\nnp.savetxt(OUT_DIR / \"kmeans_centroids.csv\", centroids,\n           delimiter=\",\", header=header, comments=\"\", fmt=\"%.2f\")\nnp.savetxt(OUT_DIR / \"kmeans_error.csv\",\n           np.stack([np.arange(len(errors)), errors], 1),\n           delimiter=\",\", header=\"iter,error\", comments=\"\", fmt=\"%.1f\")\n\nfig, axes = plt.subplots(1, 2, figsize=(9, 3.4))\nax = axes[0]\nax.plot(X[labels == 0, 0], X[labels == 0, 1], \"o\", color=\"tab:blue\",\n        ms=3, label=\"cold-season cluster\")\nax.plot(X[labels == 1, 0], X[labels == 1, 1], \"^\", mfc=\"none\",\n        mec=\"tab:red\", ms=4, label=\"warm-season cluster\")\nax.plot(centroids[:, 0], centroids[:, 1], \"kx\", ms=10, mew=2,\n        label=\"cluster centroids\")\nax.set_xlabel(\"minimum temperature of the day (\u00b0C)\")\nax.set_ylabel(\"maximum temperature of the day (\u00b0C)\")\nax.set_title(\"k-means with k = 2 on the 366 days of 2024 (Krems)\")\nax.legend(frameon=False, fontsize=8)\nax = axes[1]\nax.plot(np.arange(len(errors)), errors, \"k.-\")\nax.set_xlabel(\"Lloyd iteration\")\nax.set_ylabel(\"clustering error\")\nax.set_title(\"monotone descent to a fixed point\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"kmeans.png\", dpi=150)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-image]** Image compression and image segmentation on a subsampled photo of the Oetscher massif (assets/oetscher.jpg): k-means on the pixel colors with k = 4 replaces each pixel's color by the nearest palette color (compression factor ~12), and with k = 2 partitions the pixels into a sky-and-mountain region and a vegetation region."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# A photo of the Oetscher massif, subsampled by keeping every 36th pixel\n# in each direction: the blue sky and gray summit, the dark forest and\n# the bright meadow give a small image whose pixel colors k-means can\n# quantize (compression) and partition (segmentation).\nfrom matplotlib.image import imread\n\nSUBSAMPLE = 36\nphoto = imread(OUT_DIR.parent / \"assets\" / \"oetscher.jpg\") / 255.0\nsub = photo[::SUBSAMPLE, ::SUBSAMPLE]\ncolors = sub.reshape(-1, 3)                    # one RGB vector per pixel\nluminance = colors.mean(axis=1)\ncheck(\"[B-image] subsampling keeps every 36th pixel\",\n      sub.shape == (97, 129, 3))\nprint(f\"  photo subsampled to {sub.shape[0]} x {sub.shape[1]} pixels\")\n\n\ndef lloyd_colors(k):\n    \"\"\"Lloyd's algorithm on the pixel colors; deterministic init at the\n    colors of the pixels whose luminances sit at k evenly spaced\n    quantiles.\"\"\"\n    qs = np.quantile(luminance, np.linspace(0.0, 1.0, k))\n    cents = np.array([colors[np.abs(luminance - q).argmin()] for q in qs])\n    labs = assign(colors, cents)\n    while True:\n        cents = np.stack([colors[labs == c].mean(axis=0) for c in range(k)])\n        new = assign(colors, cents)\n        if np.array_equal(new, labs):\n            return cents, labs\n        labs = new\n\n\nK_PALETTE = 4\npalette, plabels = lloyd_colors(K_PALETTE)\ncompressed = palette[plabels].reshape(sub.shape)\nnpix = colors.shape[0]\nbits_orig = 24 * npix                          # 8-bit RGB per pixel\nbits_comp = 2 * npix + K_PALETTE * 24          # 2-bit index + palette\nfactor = bits_orig / bits_comp\ncheck(\"[B-image] 4-color palette compresses by a factor of about 12\",\n      11.5 < factor < 12.0)\n\nmask_cents, mlabels = lloyd_colors(2)\n# the bluer of the two clusters collects the sky and the gray summit,\n# the other the vegetation\nblueness = [float((colors[mlabels == c][:, 2]\n                   - colors[mlabels == c][:, 0]).mean()) for c in (0, 1)]\nsky = int(np.argmax(blueness))\nmask = (mlabels == sky).reshape(sub.shape[:2])\ncheck(\"[B-image] the sky-and-mountain cluster is markedly bluer\",\n      blueness[sky] > 0.3 > 0.1 > blueness[1 - sky])\ncheck(\"[B-image] both regions are present\",\n      0.1 < mask.mean() < 0.9)\n\nUPSCALE = 6                                    # keep the pixels crisp\n\n\ndef save_pixels(img01, path):\n    from PIL import Image\n    rgb8 = (np.clip(img01, 0.0, 1.0) * 255.0 + 0.5).astype(np.uint8)\n    im = Image.fromarray(rgb8)\n    im = im.resize((im.width * UPSCALE, im.height * UPSCALE), Image.NEAREST)\n    im.save(path)\n\n\nsave_pixels(sub, OUT_DIR / \"kmeans_oetscher_original.png\")\nsave_pixels(compressed, OUT_DIR / \"kmeans_oetscher_compressed.png\")\nsave_pixels(np.repeat(mask[:, :, None], 3, axis=2).astype(float),\n            OUT_DIR / \"kmeans_oetscher_mask.png\")\nprint(f\"  wrote kmeans_oetscher_original/compressed/mask.png \"\n      f\"(compression factor {factor:.1f})\")\n\nfailed = [name for name, ok in report if not ok]\nprint(f\"{len(report) - len(failed)}/{len(report)} checks passed\"\n      + (f\", FAILED: {failed}\" if failed else \"\"))"
  }
 ]
}