{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "gmm.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Gaussian mixture model \u2014 Python demo\n\nNumerical companion to the entry [Gaussian mixture model](https://dictionaryofml.org/terms/gmm.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 'gmm' (Gaussian mixture model). A photograph of the Oetscher massif (assets/oetscher.jpg, the same photograph the 'kmeans' entry uses) is cut into square patches. Each patch is a data point whose feature vector holds two numbers, how green and how blue the patch is on average, so the covariance matrix of a component is a genuine matrix rather than a single number and each component shows as an ellipse in the plane.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/gmm.py`](https://dictionaryofml.org/terms/gmm.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(), \"gmm.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"A three-component GMM on the greenness and blueness of image patches:\nmeadow, foliage and sky, with graded membership where they meet.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'gmm' (Gaussian mixture\nmodel).  A photograph of the Oetscher massif (assets/oetscher.jpg, the\nsame photograph the 'kmeans' entry uses) is cut into square patches.\nEach patch is a data point whose feature vector holds two numbers, how\ngreen and how blue the patch is on average, so the covariance matrix of\na component is a genuine matrix rather than a single number and each\ncomponent shows as an ellipse in the plane.\n\nThe demo checks the entry's claims: the cluster probabilities sum to\none; the three components are the sunlit meadow, the dark foliage and\nthe sky with the mountain; each covariance matrix is symmetric positive semi-definite;\nthe posterior distribution grades the membership of a patch in each\ncluster instead of assigning it to exactly one; and the EM algorithm\nreduces to k-means when the cluster probabilities are equal and every\ncovariance matrix is a shrinking multiple of the identity matrix.\n\nDeterministic: the components are initialized by splitting the patches\ninto three equal parts ordered by blueness (no randomness).\nSelf-contained: numpy + matplotlib only.\n\nBlocks\n------\n[B-patches]    Cut the photograph into 64x64 patches and measure how\n               green and how blue each one is; check the patch count.\n[B-em]         Fit the three-component GMM by the EM algorithm: E-step\n               (posterior probability of each cluster index per patch),\n               M-step (re-weighted cluster probabilities, means and\n               covariance matrices).  Check that the cluster\n               probabilities sum to one and that the fit stops changing.\n[B-components] The three components are the sunlit meadow, the dark\n               foliage, and the sky with the mountain.  Check each against the average color\n               of its patches and against where they sit in the\n               photograph, and check that each covariance matrix is\n               symmetric positive semi-definite.\n[B-soft]       Soft clustering: the posterior distribution grades\n               membership.  Check that most patches are graded\n               decisively and that the rest are shared between clusters.\n[B-images]     The photograph with the patch grid drawn on it, the\n               photograph at patch resolution, and one copy per\n               cluster whose brightness is that cluster's posterior\n               probability.  Because the posterior probabilities of a\n               patch sum to one, the three copies add back up to the\n               photograph.\n[B-kmeans]     With equal cluster probabilities and covariance matrices\n               that shrink to a multiple of the identity matrix, the\n               posterior concentrates on the nearest mean and the EM\n               update becomes the k-means update, so the algorithm\n               reduces to k-means.  Check both.\n[B-plot]       Write the patch scatter and the three component ellipses\n               for the entry's figure, plus the preview.\n\nOutputs\n-------\ngmm_patches.csv  : green, blue -- every patch\ngmm_points.csv   : x1, x2 -- the patches graded decisively\ngmm_between.csv  : x1, x2 -- the patches shared between clusters\ngmm_ellipse1.csv : x1, x2 -- contour of the meadow component\ngmm_ellipse2.csv : x1, x2 -- contour of the foliage component\ngmm_ellipse3.csv : x1, x2 -- contour of the sky component\ngmm_oetscher_raster.png   : the photograph with the patch grid drawn on it\ngmm_oetscher_original.png : the photograph at patch resolution\ngmm_oetscher_vegetation.png : brightness = posterior of the vegetation cluster\ngmm_oetscher_mountain.png   : brightness = posterior of the mountain cluster\ngmm_oetscher_sky.png        : brightness = posterior of the sky cluster\ngmm.png          : preview (checking only) -- the patches with the\n                   three component ellipses drawn over them\n\"\"\"\n\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom matplotlib.image import imread\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-patches]** Cut the photograph into 64x64 patches and measure how green and how blue each one is; check the patch count."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "PATCH = 128\nphoto = imread(OUT_DIR.parent / \"assets\" / \"oetscher.jpg\") / 255.0\nrows = photo.shape[0] // PATCH\ncols = photo.shape[1] // PATCH\ntiles = photo[:rows * PATCH, :cols * PATCH].reshape(\n    rows, PATCH, cols, PATCH, 3).mean(axis=(1, 3))\nrgb = tiles.reshape(-1, 3)\ngreen = rgb[:, 1]                   # average greenness, on a 0 to 1 scale\nblue = rgb[:, 2]                    # average blueness, on the same scale\nX = np.column_stack([green, blue])\nwith open(OUT_DIR / \"gmm_patches.csv\", \"w\") as f:\n    f.write(\"green,blue\\n\")\n    for a, b in X:\n        f.write(f\"{a:.4f},{b:.4f}\\n\")\nprint(f\"  photograph cut into {rows} x {cols} = {len(X)} patches \"\n      f\"of {PATCH} x {PATCH} pixels\")\ncheck(\"[B-patches] every patch has a two-number feature vector\",\n      X.shape == (rows * cols, 2))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-em]** Fit the three-component GMM by the EM algorithm: E-step (posterior probability of each cluster index per patch), M-step (re-weighted cluster probabilities, means and covariance matrices). Check that the cluster probabilities sum to one and that the fit stops changing."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "NRCLUSTER = 3\n\n\ndef normal_pdf(A, mean, cov):\n    \"\"\"Density of the multivariate normal distribution at each row of A.\"\"\"\n    d = A.shape[1]\n    diff = A - mean\n    quad = np.einsum(\"ij,jk,ik->i\", diff, np.linalg.inv(cov), diff)\n    return np.exp(-0.5 * quad) / np.sqrt(((2.0 * np.pi) ** d)\n                                         * np.linalg.det(cov))\n\n\ndef fit_quality(A, p, means, covs):\n    \"\"\"Negative log of the likelihood of the whole dataset.\"\"\"\n    mix = sum(p[c] * normal_pdf(A, means[c], covs[c]) for c in range(len(p)))\n    return float(-np.log(mix).sum())\n\n\ngroups = np.array_split(np.argsort(X[:, 1]), NRCLUSTER)   # ordered by blueness\nmeans = np.stack([X[g].mean(axis=0) for g in groups])\ncovs = np.stack([np.cov(X[g].T) for g in groups])\np = np.full(NRCLUSTER, 1.0 / NRCLUSTER)\n\ntrace = [fit_quality(X, p, means, covs)]\nfor _ in range(300):\n    # E-step: posterior probability of each cluster index per patch\n    joint = np.stack([p[c] * normal_pdf(X, means[c], covs[c])\n                      for c in range(NRCLUSTER)])\n    posterior = joint / joint.sum(axis=0)\n    # M-step: re-weighted cluster probabilities, means, covariance matrices\n    weight = posterior.sum(axis=1)\n    p = weight / len(X)\n    means = (posterior @ X) / weight[:, None]\n    covs = np.stack([\n        (posterior[c][:, None] * (X - means[c])).T @ (X - means[c]) / weight[c]\n        for c in range(NRCLUSTER)])\n    trace.append(fit_quality(X, p, means, covs))\n\norder = np.argsort(means[:, 1])                  # greenest first, sky last\np, means, covs = p[order], means[order], covs[order]\nposterior = posterior[order]\nlabel = posterior.argmax(axis=0)\n\ncheck(\"[B-em] the cluster probabilities sum to one\",\n      abs(p.sum() - 1.0) < 1e-12)\ncheck(\"[B-em] the fit stops changing (last update tiny)\",\n      abs(trace[-1] - trace[-2]) < 1e-9)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-components]** The three components are the sunlit meadow, the dark foliage, and the sky with the mountain. Check each against the average color of its patches and against where they sit in the photograph, and check that each covariance matrix is symmetric positive semi-definite."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "NAMES = (\"vegetation\", \"mountain and haze\", \"bright sky\")\nwhere = np.argwhere(np.ones((rows, cols), dtype=bool))    # (row, col) per patch\nfor c in range(NRCLUSTER):\n    sel = label == c\n    print(f\"  {NAMES[c]:<14} probability {p[c]:.2f}, \"\n          f\"greenness {means[c, 0]:+.3f}, blueness {means[c, 1]:+.3f}, \"\n          f\"average color RGB {np.round(rgb[sel].mean(axis=0), 2)}, \"\n          f\"{(where[sel][:, 0] < rows / 2).mean():.0%} of its patches in the \"\n          f\"upper half\")\ncheck(\"[B-components] every average lies on the 0 to 1 scale\",\n      float(X.min()) >= 0.0 and float(X.max()) <= 1.0)\ncheck(\"[B-components] the sky is the bluest component\",\n      means[2, 1] == means[:, 1].max())\ncheck(\"[B-components] the sky sits in the upper half of the photograph\",\n      float((where[label == 2][:, 0] < rows / 2).mean()) > 0.8)\ncheck(\"[B-components] the vegetation is the least blue and sits in the \"\n      \"lower half\", means[0, 1] == means[:, 1].min()\n      and float((where[label == 0][:, 0] >= rows / 2).mean()) > 0.8)\ncheck(\"[B-components] each covariance matrix is symmetric\",\n      all(np.allclose(c, c.T) for c in covs))\ncheck(\"[B-components] each covariance matrix is positive semi-definite\",\n      all(np.linalg.eigvalsh(c).min() > 0 for c in covs))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-soft]** Soft clustering: the posterior distribution grades membership. Check that most patches are graded decisively and that the rest are shared between clusters."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "top = posterior.max(axis=0)\nshared = top < 0.8\nprint(f\"  {int((~shared).sum())} of {len(X)} patches are graded above 0.8 \"\n      f\"for one cluster; {int(shared.sum())} are shared\")\ncheck(\"[B-soft] the posterior distribution of each patch sums to one\",\n      np.allclose(posterior.sum(axis=0), 1.0))\ncheck(\"[B-soft] most patches are graded decisively\", (~shared).mean() > 0.8)\ncheck(\"[B-soft] a shared patch is not assigned to exactly one cluster\",\n      bool(np.all(posterior[:, shared].max(axis=0) < 0.8)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-images]** The photograph with the patch grid drawn on it, the photograph at patch resolution, and one copy per cluster whose brightness is that cluster's posterior probability. Because the posterior probabilities of a patch sum to one, the three copies add back up to the photograph."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def save_image(arr, name, zoom=6):\n    \"\"\"Write an RGB array as a PNG, enlarged so the patches stay visible.\"\"\"\n    img = np.clip(arr, 0.0, 1.0).repeat(zoom, axis=0).repeat(zoom, axis=1)\n    plt.imsave(OUT_DIR / name, img)\n\n\n# the photograph at full detail with the patch grid drawn on it, so\n# the reader can see what one data point covers\nSHRINK = 4\nview = photo[:rows * PATCH, :cols * PATCH:, :][::SHRINK, ::SHRINK].copy()\ncell = PATCH // SHRINK\nview[::cell, :, :] = 1.0                       # horizontal rules\nview[:, ::cell, :] = 1.0                       # vertical rules\nview[-1, :, :] = 1.0\nview[:, -1, :] = 1.0\nsave_image(view, \"gmm_oetscher_raster.png\", zoom=1)\ncheck(\"[B-images] the raster has one cell per patch\",\n      view.shape[:2] == (rows * cell, cols * cell))\n\nsave_image(tiles, \"gmm_oetscher_original.png\")\nCOPIES = (\"gmm_oetscher_vegetation.png\", \"gmm_oetscher_mountain.png\",\n          \"gmm_oetscher_sky.png\")\ndimmed = [tiles * posterior[c].reshape(rows, cols, 1)\n          for c in range(NRCLUSTER)]\nfor fname, arr in zip(COPIES, dimmed):\n    save_image(arr, fname)\ncheck(\"[B-images] one copy of the photograph per cluster\",\n      all((OUT_DIR / f).exists() for f in COPIES))\ncheck(\"[B-images] the three copies add back up to the photograph\",\n      np.allclose(sum(dimmed), tiles))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-kmeans]** With equal cluster probabilities and covariance matrices that shrink to a multiple of the identity matrix, the posterior concentrates on the nearest mean and the EM update becomes the k-means update, so the algorithm reduces to k-means. Check both."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def kmeans(A, cents):\n    \"\"\"Lloyd's algorithm from the given starting means.\"\"\"\n    while True:\n        lab = ((A[:, None, :] - cents) ** 2).sum(axis=2).argmin(axis=1)\n        new = np.stack([A[lab == c].mean(axis=0) for c in range(len(cents))])\n        if np.allclose(new, cents):\n            return lab, new\n        cents = new\n\n\ndef spherical_em(A, cents, var):\n    \"\"\"EM with equal cluster probabilities and covariance matrices var*I.\"\"\"\n    post = None\n    for _ in range(300):\n        d2 = ((A[:, None, :] - cents) ** 2).sum(axis=2)\n        logp = -d2 / (2.0 * var)\n        post = np.exp(logp - logp.max(axis=1, keepdims=True))\n        post /= post.sum(axis=1, keepdims=True)\n        new = (post.T @ A) / post.sum(axis=0)[:, None]\n        if np.allclose(new, cents):\n            break\n        cents = new\n    return post.argmax(axis=1), cents, post\n\n\nstart = np.stack([X[g].mean(axis=0) for g in groups])\nlab_km, cent_km = kmeans(X, start.copy())\nfor var in (1e-2, 1e-5):\n    lab_em, cent_em, post = spherical_em(X, start.copy(), var)\n    agree = float((lab_em == lab_km).mean())\n    hard = float((post.max(axis=1) > 0.99).mean())\n    print(f\"  variance {var:g}: {hard:.0%} of posteriors above 0.99, \"\n          f\"assignment agrees with k-means on {agree:.1%} of patches\")\n    if var == 1e-5:\n        check(\"[B-kmeans] the posterior becomes a hard assignment\", hard > 0.99)\n        check(\"[B-kmeans] the assignment agrees with k-means\", agree > 0.99)\n        check(\"[B-kmeans] the means agree with the k-means centroids\",\n              np.abs(cent_em - cent_km).max() < 0.01)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Write the patch scatter and the three component ellipses for the entry's figure, plus the preview."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def ellipse(mean, cov, sigma=2.0, n=200):\n    \"\"\"Contour at `sigma` standard deviations of one component.\"\"\"\n    vals, vecs = np.linalg.eigh(cov)\n    t = np.linspace(0.0, 2.0 * np.pi, n)\n    circle = np.stack([np.cos(t), np.sin(t)])\n    return (mean[:, None] + sigma * vecs @ (np.sqrt(vals)[:, None]\n                                            * circle)).T\n\n\nfor name, sel in ((\"gmm_points.csv\", ~shared), (\"gmm_between.csv\", shared)):\n    with open(OUT_DIR / name, \"w\") as f:\n        f.write(\"x1,x2\\n\")\n        for a, b in X[sel]:\n            f.write(f\"{a:.4f},{b:.4f}\\n\")\nfor c, name in enumerate((\"gmm_ellipse1.csv\", \"gmm_ellipse2.csv\",\n                          \"gmm_ellipse3.csv\")):\n    with open(OUT_DIR / name, \"w\") as f:\n        f.write(\"x1,x2\\n\")\n        for a, b in ellipse(means[c], covs[c]):\n            f.write(f\"{a:.4f},{b:.4f}\\n\")\n\nfig, ax = plt.subplots(figsize=(6.2, 5.0))\nax.plot(X[~shared, 0], X[~shared, 1], \"o\", color=\"0.6\", markersize=2.2,\n        linestyle=\"none\", label=\"patch graded to one cluster\")\nax.plot(X[shared, 0], X[shared, 1], \"^\", color=\"black\", markersize=4,\n        markerfacecolor=\"none\", linestyle=\"none\", label=\"patch shared\")\nfor c, style in enumerate((\"-\", \"--\", \":\")):\n    e = ellipse(means[c], covs[c])\n    ax.plot(e[:, 0], e[:, 1], color=\"black\", linewidth=1.7, linestyle=style,\n            label=NAMES[c])\nax.set_xlabel(\"average greenness of the patch\")\nax.set_ylabel(\"average blueness of the patch\")\nax.set_title(\"Patches of the Oetscher photograph and a three-component GMM\")\nax.legend(frameon=False, loc=\"upper right\", fontsize=\"small\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"gmm.png\", dpi=150)\nplt.close(fig)\ncheck(\"[B-plot] the three ellipses and the patch scatter were written\",\n      (OUT_DIR / \"gmm_ellipse3.csv\").exists()\n      and (OUT_DIR / \"gmm_points.csv\").exists())\n\npassed = sum(1 for _, ok in report if ok)\nprint(f\"\\n{passed}/{len(report)} checks pass\")"
  }
 ]
}