{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Choosing K from the spectrum\n",
    "\n",
    "The singular value spectrum and the two criteria that read it, which do not always agree.\n",
    "\n",
    "From chapter 7, [Count Vectors, PPMI, and SVD](https://nlp.jcrlabz.com/book/countvectors/), of the course notes.\n",
    "\n",
    "Source: `book/code/worked_examples/svd_rank.py`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What you need\n",
    "\n",
    "`pip install numpy`. The next cell does it. Colab usually has numpy already, but nothing here assumes that."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "%pip install -q numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The script\n",
    "\n",
    "Everything the example defines, in one cell. Run it and the definitions\n",
    "are live in the notebook."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "\"\"\"Choosing K from the singular value spectrum, as Chapter 7 does it.\n",
    "\n",
    "The singular values tell you how many dimensions your data actually has. Each\n",
    "one carries a share of the total, and the running total says how much you keep\n",
    "if you cut at K.\n",
    "\n",
    "    contribution   r_i = sigma_i / sum(sigma)\n",
    "    cumulative     R_K = sum(sigma_1..K) / sum(sigma)\n",
    "\n",
    "    python3 svd_rank.py                 # the book's 8x8 example\n",
    "    python3 svd_rank.py --squared       # the variance convention\n",
    "    python3 svd_rank.py --threshold 0.9 # cut wherever you like\n",
    "\n",
    "The demonstration matrix has two hidden groups, animals and vehicles, sharing\n",
    "no context. Watch the spectrum discover that without being told.\n",
    "\n",
    "Install:  pip install numpy\n",
    "\"\"\"\n",
    "\n",
    "import argparse\n",
    "\n",
    "try:\n",
    "    import numpy as np\n",
    "except ImportError:\n",
    "    raise SystemExit(\"this one needs numpy:  pip install numpy\")\n",
    "\n",
    "\n",
    "WORDS = [\"dog\", \"cat\", \"car\", \"truck\", \"pet\", \"vet\", \"drive\", \"fuel\"]\n",
    "CONTEXTS = [\"pet\", \"feed\", \"fur\", \"vet\", \"drive\", \"fuel\", \"road\", \"tyre\"]\n",
    "\n",
    "# Rows 1-2 and 5-6 are animals, rows 3-4 and 7-8 are vehicles. The two blocks\n",
    "# share no nonzero column, so the true structure is two dimensional.\n",
    "M = np.array([\n",
    "    [4, 3, 3, 2, 0, 0, 0, 0],   # dog\n",
    "    [4, 3, 4, 2, 0, 0, 0, 0],   # cat\n",
    "    [0, 0, 0, 0, 4, 3, 3, 2],   # car\n",
    "    [0, 0, 0, 0, 4, 3, 4, 2],   # truck\n",
    "    [3, 2, 2, 3, 0, 0, 0, 0],   # pet\n",
    "    [2, 1, 1, 4, 0, 0, 0, 0],   # vet\n",
    "    [0, 0, 0, 0, 3, 2, 3, 2],   # drive\n",
    "    [0, 0, 0, 0, 2, 3, 2, 3],   # fuel\n",
    "], dtype=float)\n",
    "\n",
    "\n",
    "def spectrum(matrix, squared=False):\n",
    "    \"\"\"Return (singular values, share of each, cumulative share).\"\"\"\n",
    "    sigma = np.linalg.svd(matrix, compute_uv=False)\n",
    "    weight = sigma ** 2 if squared else sigma\n",
    "    share = weight / weight.sum()\n",
    "    return sigma, share, np.cumsum(share)\n",
    "\n",
    "\n",
    "def choose_k(cumulative, threshold):\n",
    "    \"\"\"Smallest K whose cumulative share reaches the threshold.\"\"\"\n",
    "    return int(np.searchsorted(cumulative, threshold) + 1)\n",
    "\n",
    "\n",
    "def elbow(share):\n",
    "    \"\"\"Index after the largest drop between consecutive shares.\n",
    "\n",
    "    A blunt instrument, and it agrees with the eye on clean spectra. On real\n",
    "    corpora the curve is smooth and you should use a threshold instead.\n",
    "    \"\"\"\n",
    "    drops = share[:-1] - share[1:]\n",
    "    return int(np.argmax(drops) + 1)\n",
    "\n",
    "\n",
    "def main():\n",
    "    ap = argparse.ArgumentParser(description=__doc__)\n",
    "    ap.add_argument(\"--squared\", action=\"store_true\",\n",
    "                    help=\"use sigma^2, the explained-variance convention\")\n",
    "    ap.add_argument(\"--threshold\", type=float, default=0.8)\n",
    "    a = ap.parse_args()\n",
    "\n",
    "    sigma, share, cum = spectrum(M, a.squared)\n",
    "    convention = \"sigma^2 / sum(sigma^2)\" if a.squared else \"sigma / sum(sigma)\"\n",
    "    print(f\"\\n{M.shape[0]}x{M.shape[1]} co-occurrence matrix\")\n",
    "    print(f\"contribution measured as {convention}\\n\")\n",
    "\n",
    "    k_elbow = elbow(share)\n",
    "    k_thresh = choose_k(cum, a.threshold)\n",
    "\n",
    "    print(f\"  {'i':>2}{'sigma':>10}{'share':>9}{'cumulative':>12}\")\n",
    "    print(\"  \" + \"-\" * 33)\n",
    "    for i, (s, r, c) in enumerate(zip(sigma, share, cum), 1):\n",
    "        mark = \"\"\n",
    "        if i == k_elbow:\n",
    "            mark += \"  <-- elbow\"\n",
    "        if i == k_thresh:\n",
    "            mark += f\"  <-- reaches {a.threshold:.0%}\"\n",
    "        print(f\"  {i:>2}{s:>10.3f}{r:>8.1%}{c:>11.1%}{mark}\")\n",
    "\n",
    "    print(f\"\\n  sum of sigma = {sigma.sum():.3f}\")\n",
    "    print(f\"  elbow says          K = {k_elbow}\")\n",
    "    print(f\"  {a.threshold:.0%} threshold says  K = {k_thresh}\")\n",
    "\n",
    "    print(\"\\n  This matrix has two hidden groups, animals and vehicles, and\")\n",
    "    print(\"  they share no context. Nobody told the decomposition that. It\")\n",
    "    print(\"  read the number of themes off the data.\")\n",
    "    if not a.squared:\n",
    "        _, sq_share, sq_cum = spectrum(M, squared=True)\n",
    "        print(f\"\\n  For comparison, the squared convention gives\"\n",
    "              f\" {sq_share[0]:.1%} for the first\")\n",
    "        print(f\"  dimension and {sq_cum[1]:.1%} for the first two. Squaring\"\n",
    "              f\" exaggerates the\")\n",
    "        print(\"  lead of the top values. Both conventions are used; say which.\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running it\n",
    "\n",
    "The script is a command line tool, so each run below sets the arguments\n",
    "and calls `main()`. The heading shows the equivalent shell command."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "def run(*args):\n",
    "    sys.argv = [\"svd_rank.py\", *args]\n",
    "    main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `python3 svd_rank.py`\n",
    "\n",
    "The elbow says K = 2, an 80 per cent threshold says K = 3. The truth is 2: the matrix has exactly two hidden groups."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `python3 svd_rank.py --squared`\n",
    "\n",
    "Switch to the explained variance convention and the threshold changes its mind. Thresholds are conventions, not findings."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "run(\"--squared\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Read the chapter this comes from: [Count Vectors, PPMI, and SVD](https://nlp.jcrlabz.com/book/countvectors/)."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "svd_rank.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}