Week_22-Clustering_text_data_with_k_means.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.0"
    },
    "colab": {
      "name": "Clustering text data with k-means_CLU03-NB01.ipynb",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qhKCDDoyJoqS"
      },
      "source": [
        "# k-means with text data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WwqnnAD1Joqm"
      },
      "source": [
        "In this assignment you will\n",
        "* Cluster Wikipedia documents using k-means\n",
        "* Explore the role of random initialization on the quality of the clustering\n",
        "* Explore how results differ after changing the number of clusters\n",
        "* Evaluate clustering, both quantitatively and qualitatively\n",
        "\n",
        "When properly executed, clustering uncovers valuable insights from a set of unlabeled documents."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mrkeu3YOJoqr"
      },
      "source": [
        "**Note to Amazon EC2 users**: To conserve memory, make sure to stop all the other notebooks before running this notebook."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "12MJPyhqJoqt"
      },
      "source": [
        "## Import necessary packages"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "205HpPr_JtOS",
        "outputId": "7b92ed55-07dd-43d5-892d-072742954651"
      },
      "source": [
        "pip install turicreate"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting turicreate\n",
            "  Downloading turicreate-6.4.1-cp37-cp37m-manylinux1_x86_64.whl (92.0 MB)\n",
            "\u001b[K     |████████████████████████████████| 92.0 MB 5.3 kB/s \n",
            "\u001b[?25hRequirement already satisfied: decorator>=4.0.9 in /usr/local/lib/python3.7/dist-packages (from turicreate) (4.4.2)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from turicreate) (1.19.5)\n",
            "Requirement already satisfied: pandas>=0.23.2 in /usr/local/lib/python3.7/dist-packages (from turicreate) (1.1.5)\n",
            "Collecting coremltools==3.3\n",
            "  Downloading coremltools-3.3-cp37-none-manylinux1_x86_64.whl (3.5 MB)\n",
            "\u001b[K     |████████████████████████████████| 3.5 MB 28.0 MB/s \n",
            "\u001b[?25hCollecting tensorflow<2.1.0,>=2.0.0\n",
            "  Downloading tensorflow-2.0.4-cp37-cp37m-manylinux2010_x86_64.whl (86.4 MB)\n",
            "\u001b[K     |████████████████████████████████| 86.4 MB 28 kB/s \n",
            "\u001b[?25hRequirement already satisfied: pillow>=5.2.0 in /usr/local/lib/python3.7/dist-packages (from turicreate) (7.1.2)\n",
            "Requirement already satisfied: six>=1.10.0 in /usr/local/lib/python3.7/dist-packages (from turicreate) (1.15.0)\n",
            "Collecting prettytable==0.7.2\n",
            "  Downloading prettytable-0.7.2.zip (28 kB)\n",
            "Requirement already satisfied: requests>=2.9.1 in /usr/local/lib/python3.7/dist-packages (from turicreate) (2.23.0)\n",
            "Collecting resampy==0.2.1\n",
            "  Downloading resampy-0.2.1.tar.gz (322 kB)\n",
            "\u001b[K     |████████████████████████████████| 322 kB 45.6 MB/s \n",
            "\u001b[?25hRequirement already satisfied: scipy>=1.1.0 in /usr/local/lib/python3.7/dist-packages (from turicreate) (1.4.1)\n",
            "Collecting numba<0.51.0\n",
            "  Downloading numba-0.50.1-cp37-cp37m-manylinux2014_x86_64.whl (3.6 MB)\n",
            "\u001b[K     |████████████████████████████████| 3.6 MB 51.3 MB/s \n",
            "\u001b[?25hRequirement already satisfied: protobuf>=3.1.0 in /usr/local/lib/python3.7/dist-packages (from coremltools==3.3->turicreate) (3.17.3)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from numba<0.51.0->turicreate) (57.4.0)\n",
            "Collecting llvmlite<0.34,>=0.33.0.dev0\n",
            "  Downloading llvmlite-0.33.0-cp37-cp37m-manylinux1_x86_64.whl (18.3 MB)\n",
            "\u001b[K     |████████████████████████████████| 18.3 MB 1.3 MB/s \n",
            "\u001b[?25hRequirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=0.23.2->turicreate) (2.8.2)\n",
            "Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.7/dist-packages (from pandas>=0.23.2->turicreate) (2018.9)\n",
            "Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.7/dist-packages (from requests>=2.9.1->turicreate) (2021.5.30)\n",
            "Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /usr/local/lib/python3.7/dist-packages (from requests>=2.9.1->turicreate) (1.24.3)\n",
            "Requirement already satisfied: idna<3,>=2.5 in /usr/local/lib/python3.7/dist-packages (from requests>=2.9.1->turicreate) (2.10)\n",
            "Requirement already satisfied: chardet<4,>=3.0.2 in /usr/local/lib/python3.7/dist-packages (from requests>=2.9.1->turicreate) (3.0.4)\n",
            "Collecting numpy\n",
            "  Downloading numpy-1.18.5-cp37-cp37m-manylinux1_x86_64.whl (20.1 MB)\n",
            "\u001b[K     |████████████████████████████████| 20.1 MB 1.5 MB/s \n",
            "\u001b[?25hCollecting gast==0.2.2\n",
            "  Downloading gast-0.2.2.tar.gz (10 kB)\n",
            "Collecting tensorboard<2.1.0,>=2.0.0\n",
            "  Downloading tensorboard-2.0.2-py3-none-any.whl (3.8 MB)\n",
            "\u001b[K     |████████████████████████████████| 3.8 MB 42.5 MB/s \n",
            "\u001b[?25hCollecting h5py<=2.10.0\n",
            "  Downloading h5py-2.10.0-cp37-cp37m-manylinux1_x86_64.whl (2.9 MB)\n",
            "\u001b[K     |████████████████████████████████| 2.9 MB 51.2 MB/s \n",
            "\u001b[?25hRequirement already satisfied: termcolor>=1.1.0 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (1.1.0)\n",
            "Requirement already satisfied: grpcio>=1.8.6 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (1.39.0)\n",
            "Requirement already satisfied: google-pasta>=0.1.6 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (0.2.0)\n",
            "Requirement already satisfied: astor>=0.6.0 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (0.8.1)\n",
            "Collecting tensorflow-estimator<2.1.0,>=2.0.0\n",
            "  Downloading tensorflow_estimator-2.0.1-py2.py3-none-any.whl (449 kB)\n",
            "\u001b[K     |████████████████████████████████| 449 kB 58.7 MB/s \n",
            "\u001b[?25hRequirement already satisfied: absl-py>=0.7.0 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (0.12.0)\n",
            "Collecting keras-applications>=1.0.8\n",
            "  Downloading Keras_Applications-1.0.8-py3-none-any.whl (50 kB)\n",
            "\u001b[K     |████████████████████████████████| 50 kB 9.4 MB/s \n",
            "\u001b[?25hRequirement already satisfied: opt-einsum>=2.3.2 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (3.3.0)\n",
            "Requirement already satisfied: wheel>=0.26 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (0.37.0)\n",
            "Requirement already satisfied: keras-preprocessing>=1.0.5 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (1.1.2)\n",
            "Requirement already satisfied: wrapt>=1.11.1 in /usr/local/lib/python3.7/dist-packages (from tensorflow<2.1.0,>=2.0.0->turicreate) (1.12.1)\n",
            "Requirement already satisfied: markdown>=2.6.8 in /usr/local/lib/python3.7/dist-packages (from tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (3.3.4)\n",
            "Requirement already satisfied: werkzeug>=0.11.15 in /usr/local/lib/python3.7/dist-packages (from tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (1.0.1)\n",
            "Requirement already satisfied: google-auth-oauthlib<0.5,>=0.4.1 in /usr/local/lib/python3.7/dist-packages (from tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (0.4.5)\n",
            "Requirement already satisfied: google-auth<2,>=1.6.3 in /usr/local/lib/python3.7/dist-packages (from tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (1.34.0)\n",
            "Requirement already satisfied: rsa<5,>=3.1.4 in /usr/local/lib/python3.7/dist-packages (from google-auth<2,>=1.6.3->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (4.7.2)\n",
            "Requirement already satisfied: pyasn1-modules>=0.2.1 in /usr/local/lib/python3.7/dist-packages (from google-auth<2,>=1.6.3->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (0.2.8)\n",
            "Requirement already satisfied: cachetools<5.0,>=2.0.0 in /usr/local/lib/python3.7/dist-packages (from google-auth<2,>=1.6.3->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (4.2.2)\n",
            "Requirement already satisfied: requests-oauthlib>=0.7.0 in /usr/local/lib/python3.7/dist-packages (from google-auth-oauthlib<0.5,>=0.4.1->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (1.3.0)\n",
            "Requirement already satisfied: importlib-metadata in /usr/local/lib/python3.7/dist-packages (from markdown>=2.6.8->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (4.6.4)\n",
            "Requirement already satisfied: pyasn1<0.5.0,>=0.4.6 in /usr/local/lib/python3.7/dist-packages (from pyasn1-modules>=0.2.1->google-auth<2,>=1.6.3->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (0.4.8)\n",
            "Requirement already satisfied: oauthlib>=3.0.0 in /usr/local/lib/python3.7/dist-packages (from requests-oauthlib>=0.7.0->google-auth-oauthlib<0.5,>=0.4.1->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (3.1.1)\n",
            "Requirement already satisfied: typing-extensions>=3.6.4 in /usr/local/lib/python3.7/dist-packages (from importlib-metadata->markdown>=2.6.8->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (3.7.4.3)\n",
            "Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.7/dist-packages (from importlib-metadata->markdown>=2.6.8->tensorboard<2.1.0,>=2.0.0->tensorflow<2.1.0,>=2.0.0->turicreate) (3.5.0)\n",
            "Building wheels for collected packages: prettytable, resampy, gast\n",
            "  Building wheel for prettytable (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for prettytable: filename=prettytable-0.7.2-py3-none-any.whl size=13714 sha256=38013f76b418cc63fe340ea13d1cf8d8411946a6932e593d37a3fff2df84699b\n",
            "  Stored in directory: /root/.cache/pip/wheels/b2/7f/f6/f180315b584f00445045ff1699b550fa895d09471337ce21c6\n",
            "  Building wheel for resampy (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for resampy: filename=resampy-0.2.1-py3-none-any.whl size=320860 sha256=b2506a172c29e62f130cc8eeb56c07f9ab45e1696c03b64911a15ce9e307b35d\n",
            "  Stored in directory: /root/.cache/pip/wheels/71/74/53/d5ceb7c5ee7a168c7d106041863e71ac3273f4a4677743a284\n",
            "  Building wheel for gast (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for gast: filename=gast-0.2.2-py3-none-any.whl size=7554 sha256=c641a65e000a8452658e2725cbbf1cad23c2b0287680a372e2f54e6dc86de054\n",
            "  Stored in directory: /root/.cache/pip/wheels/21/7f/02/420f32a803f7d0967b48dd823da3f558c5166991bfd204eef3\n",
            "Successfully built prettytable resampy gast\n",
            "Installing collected packages: numpy, llvmlite, h5py, tensorflow-estimator, tensorboard, numba, keras-applications, gast, tensorflow, resampy, prettytable, coremltools, turicreate\n",
            "  Attempting uninstall: numpy\n",
            "    Found existing installation: numpy 1.19.5\n",
            "    Uninstalling numpy-1.19.5:\n",
            "      Successfully uninstalled numpy-1.19.5\n",
            "  Attempting uninstall: llvmlite\n",
            "    Found existing installation: llvmlite 0.34.0\n",
            "    Uninstalling llvmlite-0.34.0:\n",
            "      Successfully uninstalled llvmlite-0.34.0\n",
            "  Attempting uninstall: h5py\n",
            "    Found existing installation: h5py 3.1.0\n",
            "    Uninstalling h5py-3.1.0:\n",
            "      Successfully uninstalled h5py-3.1.0\n",
            "  Attempting uninstall: tensorflow-estimator\n",
            "    Found existing installation: tensorflow-estimator 2.6.0\n",
            "    Uninstalling tensorflow-estimator-2.6.0:\n",
            "      Successfully uninstalled tensorflow-estimator-2.6.0\n",
            "  Attempting uninstall: tensorboard\n",
            "    Found existing installation: tensorboard 2.6.0\n",
            "    Uninstalling tensorboard-2.6.0:\n",
            "      Successfully uninstalled tensorboard-2.6.0\n",
            "  Attempting uninstall: numba\n",
            "    Found existing installation: numba 0.51.2\n",
            "    Uninstalling numba-0.51.2:\n",
            "      Successfully uninstalled numba-0.51.2\n",
            "  Attempting uninstall: gast\n",
            "    Found existing installation: gast 0.4.0\n",
            "    Uninstalling gast-0.4.0:\n",
            "      Successfully uninstalled gast-0.4.0\n",
            "  Attempting uninstall: tensorflow\n",
            "    Found existing installation: tensorflow 2.6.0\n",
            "    Uninstalling tensorflow-2.6.0:\n",
            "      Successfully uninstalled tensorflow-2.6.0\n",
            "  Attempting uninstall: resampy\n",
            "    Found existing installation: resampy 0.2.2\n",
            "    Uninstalling resampy-0.2.2:\n",
            "      Successfully uninstalled resampy-0.2.2\n",
            "  Attempting uninstall: prettytable\n",
            "    Found existing installation: prettytable 2.1.0\n",
            "    Uninstalling prettytable-2.1.0:\n",
            "      Successfully uninstalled prettytable-2.1.0\n",
            "\u001b[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.\n",
            "tensorflow-probability 0.13.0 requires gast>=0.3.2, but you have gast 0.2.2 which is incompatible.\n",
            "librosa 0.8.1 requires resampy>=0.2.2, but you have resampy 0.2.1 which is incompatible.\n",
            "datascience 0.10.6 requires folium==0.2.1, but you have folium 0.8.3 which is incompatible.\n",
            "albumentations 0.1.12 requires imgaug<0.2.7,>=0.2.5, but you have imgaug 0.2.9 which is incompatible.\u001b[0m\n",
            "Successfully installed coremltools-3.3 gast-0.2.2 h5py-2.10.0 keras-applications-1.0.8 llvmlite-0.33.0 numba-0.50.1 numpy-1.18.5 prettytable-0.7.2 resampy-0.2.1 tensorboard-2.0.2 tensorflow-2.0.4 tensorflow-estimator-2.0.1 turicreate-6.4.1\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "application/vnd.colab-display-data+json": {
              "pip_warning": {
                "packages": [
                  "numpy"
                ]
              }
            }
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LD4XIo80Joqw"
      },
      "source": [
        "from __future__ import print_function # to conform python 2.x print to python 3.x\n",
        "import turicreate\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import sys\n",
        "import os\n",
        "from scipy.sparse import csr_matrix\n",
        "from sklearn.preprocessing import OneHotEncoder, LabelEncoder\n",
        "\n",
        "%matplotlib inline"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gjoCRORPJoq2"
      },
      "source": [
        "## Load data, extract features"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nMTTc2tVJoq5"
      },
      "source": [
        "To work with text data, we must first convert the documents into numerical features. As in the first assignment, let's extract TF-IDF features for each article."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IE0KD-XGKY39",
        "outputId": "20bf05df-d85b-4b98-bbd7-28900d014f03"
      },
      "source": [
        "from google.colab import drive\n",
        "drive.mount ('/content/gdrive')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/gdrive\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k3Rht7YOJoq8"
      },
      "source": [
        "# wiki = turicreate.SFrame('people_wiki.sframe/')\n",
        "\n",
        "wiki = turicreate.SFrame('/content/gdrive/My Drive/Turicreate/Week 22/people_wiki.sframe/')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2ly8D-ulJoq_"
      },
      "source": [
        "wiki['tf_idf'] = turicreate.text_analytics.tf_idf(wiki['text'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DOYLjb_7JorF"
      },
      "source": [
        "For the remainder of the assignment, we will use sparse matrices. Sparse matrices are matrices that have a small number of nonzero entries. A good data structure for sparse matrices would only store the nonzero entries to save space and speed up computation. SciPy provides a highly-optimized library for sparse matrices. Many matrix operations available for NumPy arrays are also available for SciPy sparse matrices.\n",
        "\n",
        "We first convert the TF-IDF column (in dictionary format) into the SciPy sparse matrix format. We included plenty of comments for the curious; if you'd like, you may skip the next block and treat the function as a black box."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OwYZ1dQHJorH"
      },
      "source": [
        "def sframe_to_scipy(x, column_name):\n",
        "    '''\n",
        "    Convert a dictionary column of an SFrame into a sparse matrix format where\n",
        "    each (row_id, column_id, value) triple corresponds to the value of\n",
        "    x[row_id][column_id], where column_id is a key in the dictionary.\n",
        "       \n",
        "    Example\n",
        "    >>> sparse_matrix, map_key_to_index = sframe_to_scipy(sframe, column_name)\n",
        "    '''\n",
        "    assert type(x[column_name][0]) == dict, \\\n",
        "        'The chosen column must be dict type, representing sparse data.'\n",
        "    \n",
        "    # 1. Add a row number (id)\n",
        "    x = x.add_row_number()\n",
        "\n",
        "    # 2. Stack will transform x to have a row for each unique (row, key) pair.\n",
        "    x = x.stack(column_name, ['feature', 'value'])\n",
        "\n",
        "    # Map feature words to integers \n",
        "    unique_words = sorted(x['feature'].unique())\n",
        "    mapping = {word:i for i, word in enumerate(unique_words)}\n",
        "    x['feature_id'] = x['feature'].apply(lambda x: mapping[x])\n",
        "\n",
        "    # Create numpy arrays that contain the data for the sparse matrix.\n",
        "    row_id = np.array(x['id'])\n",
        "    col_id = np.array(x['feature_id'])\n",
        "    data = np.array(x['value'])\n",
        "    \n",
        "    width = x['id'].max() + 1\n",
        "    height = x['feature_id'].max() + 1\n",
        "    \n",
        "    # Create a sparse matrix.\n",
        "    mat = csr_matrix((data, (row_id, col_id)), shape=(width, height))\n",
        "    return mat, mapping"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jyHNQcEvJorK",
        "outputId": "a2365be9-90db-41ea-8b3d-5b0e8de597f4"
      },
      "source": [
        "%%time\n",
        "# The conversion will take about a minute or two.\n",
        "tf_idf, map_index_to_word = sframe_to_scipy(wiki, 'tf_idf')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 1min 52s, sys: 4.83 s, total: 1min 57s\n",
            "Wall time: 1min 25s\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Owhcxs4fJorM"
      },
      "source": [
        "The above matrix contains a TF-IDF score for each of the 59071 pages in the data set and each of the 547979 unique words."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "85mW2tqnJorO",
        "outputId": "75af423c-f408-44d0-c7aa-6340cf7e288d"
      },
      "source": [
        "tf_idf.shape"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(59071, 547979)"
            ]
          },
          "metadata": {},
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NH93C_-eJorQ"
      },
      "source": [
        "## Normalize all vectors"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R9nE02lRJorR"
      },
      "source": [
        "As discussed in the previous assignment, Euclidean distance can be a poor metric of similarity between documents, as it unfairly penalizes long articles. For a reasonable assessment of similarity, we should disregard the length information and use length-agnostic metrics, such as cosine distance.\n",
        "\n",
        "The k-means algorithm does not directly work with cosine distance, so we take an alternative route to remove length information: we normalize all vectors to be unit length. It turns out that Euclidean distance closely mimics cosine distance when all vectors are unit length. In particular, the squared Euclidean distance between any two vectors of length one is directly proportional to their cosine distance.\n",
        "\n",
        "We can prove this as follows. Let $\\mathbf{x}$ and $\\mathbf{y}$ be normalized vectors, i.e. unit vectors, so that $\\|\\mathbf{x}\\|=\\|\\mathbf{y}\\|=1$. Write the squared Euclidean distance as the dot product of $(\\mathbf{x} - \\mathbf{y})$ to itself:\n",
        "\\begin{align*}\n",
        "\\|\\mathbf{x} - \\mathbf{y}\\|^2 &= (\\mathbf{x} - \\mathbf{y})^T(\\mathbf{x} - \\mathbf{y})\\\\\n",
        "                              &= (\\mathbf{x}^T \\mathbf{x}) - 2(\\mathbf{x}^T \\mathbf{y}) + (\\mathbf{y}^T \\mathbf{y})\\\\\n",
        "                              &= \\|\\mathbf{x}\\|^2 - 2(\\mathbf{x}^T \\mathbf{y}) + \\|\\mathbf{y}\\|^2\\\\\n",
        "                              &= 2 - 2(\\mathbf{x}^T \\mathbf{y})\\\\\n",
        "                              &= 2(1 - (\\mathbf{x}^T \\mathbf{y}))\\\\\n",
        "                              &= 2\\left(1 - \\frac{\\mathbf{x}^T \\mathbf{y}}{\\|\\mathbf{x}\\|\\|\\mathbf{y}\\|}\\right)\\\\\n",
        "                              &= 2\\left[\\text{cosine distance}\\right]\n",
        "\\end{align*}\n",
        "\n",
        "This tells us that two **unit vectors** that are close in Euclidean distance are also close in cosine distance. Thus, the k-means algorithm (which naturally uses Euclidean distances) on normalized vectors will produce the same results as clustering using cosine distance as a distance metric.\n",
        "\n",
        "We import the [`normalize()` function](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.normalize.html) from scikit-learn to normalize all vectors to unit length."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mSsrEdCzJorU"
      },
      "source": [
        "from sklearn.preprocessing import normalize\n",
        "tf_idf = normalize(tf_idf)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iAHwO2rXJorX"
      },
      "source": [
        "## Implement k-means"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FBq08-pNJorZ"
      },
      "source": [
        "Let us implement the k-means algorithm. First, we choose an initial set of centroids. A common practice is to choose randomly from the data points.\n",
        "\n",
        "**Note:** We specify a seed here, so that everyone gets the same answer. In practice, we highly recommend to use different seeds every time (for instance, by using the current timestamp)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Aj7iXNMkJora"
      },
      "source": [
        "def get_initial_centroids(data, k, seed=None):\n",
        "    '''Randomly choose k data points as initial centroids'''\n",
        "    if seed is not None: # useful for obtaining consistent results\n",
        "        np.random.seed(seed)\n",
        "    n = data.shape[0] # number of data points\n",
        "        \n",
        "    # Pick K indices from range [0, N).\n",
        "    rand_indices = np.random.randint(0, n, k)\n",
        "    \n",
        "    # Keep centroids as dense format, as many entries will be nonzero due to averaging.\n",
        "    # As long as at least one document in a cluster contains a word,\n",
        "    # it will carry a nonzero weight in the TF-IDF vector of the centroid.\n",
        "    centroids = data[rand_indices,:].toarray()\n",
        "    \n",
        "    return centroids"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s-LWYyS5Jorc"
      },
      "source": [
        "After initialization, the k-means algorithm iterates between the following two steps:\n",
        "1. Assign each data point to the closest centroid.\n",
        "$$\n",
        "z_i \\gets \\mathrm{argmin}_j \\|\\mu_j - \\mathbf{x}_i\\|^2\n",
        "$$\n",
        "2. Revise centroids as the mean of the assigned data points.\n",
        "$$\n",
        "\\mu_j \\gets \\frac{1}{n_j}\\sum_{i:z_i=j} \\mathbf{x}_i\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hr2pyE4_Jord"
      },
      "source": [
        "In pseudocode, we iteratively do the following:\n",
        "```\n",
        "cluster_assignment = assign_clusters(data, centroids)\n",
        "centroids = revise_centroids(data, k, cluster_assignment)\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kJRvgdpsJore"
      },
      "source": [
        "### Assigning clusters"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2S6KqgB3Jorf"
      },
      "source": [
        "How do we implement Step 1 of the main k-means loop above? First import `pairwise_distances` function from scikit-learn, which calculates Euclidean distances between rows of given arrays. See [this documentation](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.pairwise_distances.html) for more information.\n",
        "\n",
        "For the sake of demonstration, let's look at documents 100 through 102 as query documents and compute the distances between each of these documents and every other document in the corpus. In the k-means algorithm, we will have to compute pairwise distances between the set of centroids and the set of documents."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TPk-6PYUJorg",
        "outputId": "75faaec7-5295-4d76-89e5-71dfddf4b332"
      },
      "source": [
        "from sklearn.metrics import pairwise_distances\n",
        "\n",
        "# Get the TF-IDF vectors for documents 100 through 102.\n",
        "queries = tf_idf[100:102,:]\n",
        "\n",
        "# Compute pairwise distances from every data point to each query vector.\n",
        "dist = pairwise_distances(tf_idf, queries, metric='euclidean')\n",
        "\n",
        "print(dist)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[1.41000789 1.36894636]\n",
            " [1.40935215 1.41023886]\n",
            " [1.39855967 1.40890299]\n",
            " ...\n",
            " [1.41108296 1.39123646]\n",
            " [1.41022804 1.31468652]\n",
            " [1.39899784 1.41072448]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3tDau1x5MDJ9",
        "outputId": "86d6cb0b-cd7e-43a0-9de8-8ef4dff1edc5"
      },
      "source": [
        "print (queries.shape)\n",
        "print (tf_idf.shape)\n",
        "print (dist.shape)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(2, 547979)\n",
            "(59071, 547979)\n",
            "(59071, 2)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RsxYQtZpJorh"
      },
      "source": [
        "More formally, `dist[i,j]` is assigned the distance between the `i`th row of `X` (i.e., `X[i,:]`) and the `j`th row of `Y` (i.e., `Y[j,:]`)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bzg3X7EmJori"
      },
      "source": [
        "**Checkpoint:** For a moment, suppose that we initialize three centroids with the first 3 rows of `tf_idf`. Write code to compute distances from each of the centroids to all data points in `tf_idf`. Then find the distance between row 430 of `tf_idf` and the second centroid and save it to `dist`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wvdqdo5nJorj",
        "outputId": "07579425-b846-4ded-d679-828fdd7230e1"
      },
      "source": [
        "# Students should write code here\n",
        "\n",
        "first_3_centroids = tf_idf[:3,:]\n",
        "distances = pairwise_distances(tf_idf, first_3_centroids, metric='euclidean')\n",
        "dist = distances[430, 1]\n",
        "\n",
        "print (dist)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "1.4071310658540348\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "aG19bOmJJork",
        "outputId": "10be563b-d930-43dd-999d-0f4f32d17776"
      },
      "source": [
        "'''Test cell'''\n",
        "if np.allclose(dist, pairwise_distances(tf_idf[430,:], tf_idf[1,:])):\n",
        "    print('Pass')\n",
        "else:\n",
        "    print('Check your code again')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Pass\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JLKEind0Jorl"
      },
      "source": [
        "**Checkpoint:** Next, given the pairwise distances, we take the minimum of the distances for each data point. Fittingly, NumPy provides an `argmin` function. See [this documentation](http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.argmin.html) for details.\n",
        "\n",
        "Read the documentation and write code to produce a 1D array whose i-th entry indicates the centroid that is the closest to the i-th data point. Use the list of distances from the previous checkpoint and save them as `distances`. The value 0 indicates closeness to the first centroid, 1 indicates closeness to the second centroid, and so forth. Save this array as `closest_cluster`.\n",
        "\n",
        "**Hint:** the resulting array should be as long as the number of data points."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "O21RHuRVJorm",
        "outputId": "d8094a52-9cdf-47cf-a591-31c1d1841ee6"
      },
      "source": [
        "# Students should write code here\n",
        "\n",
        "distances = distances.copy()\n",
        "closest_cluster = np.argmin(distances, axis=1)\n",
        "\n",
        "print (closest_cluster)\n",
        "print (closest_cluster.shape)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[0 1 2 ... 0 0 0]\n",
            "(59071,)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Q-aHEB5uJoro",
        "outputId": "624116d6-d69c-4ae3-d394-ff4acebb4701"
      },
      "source": [
        "'''Test cell'''\n",
        "reference = [list(row).index(min(row)) for row in distances]\n",
        "if np.allclose(closest_cluster, reference):\n",
        "    print('Pass')\n",
        "else:\n",
        "    print('Check your code again')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Pass\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g8aZXzKoJorp"
      },
      "source": [
        "**Checkpoint:** Let's put these steps together.  First, initialize three centroids with the first 3 rows of `tf_idf`. Then, compute distances from each of the centroids to all data points in `tf_idf`. Finally, use these distance calculations to compute cluster assignments and assign them to `cluster_assignment`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r-25hNu0Jorr"
      },
      "source": [
        "# Students should write code here\n",
        "\n",
        "first_3_centroids = tf_idf[:3,:]\n",
        "distances = pairwise_distances(tf_idf, first_3_centroids, metric='euclidean')\n",
        "cluster_assignment = np.argmin(distances, axis=1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CH0IPDU_Jorr",
        "outputId": "d3427c1b-6742-49be-894c-fa6eadfe1108"
      },
      "source": [
        "if len(cluster_assignment)==59071 and \\\n",
        "   np.array_equal(np.bincount(cluster_assignment), np.array([23061, 10086, 25924])):\n",
        "    print('Pass') # count number of data points for each cluster\n",
        "else:\n",
        "    print('Check your code again.')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Pass\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pamKcA-wJors"
      },
      "source": [
        "Now we are ready to fill in the blanks in this function:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mAiwKbl-Jort"
      },
      "source": [
        "def assign_clusters(data, centroids):\n",
        "    \n",
        "    # Compute distances between each data point and the set of centroids:\n",
        "    # Fill in the blank (RHS only)\n",
        "    distances_from_centroids = pairwise_distances(data, centroids, metric='euclidean')   # YOUR CODE HERE\n",
        "    \n",
        "    # Compute cluster assignments for each data point:\n",
        "    # Fill in the blank (RHS only)\n",
        "    cluster_assignment = np.argmin(distances_from_centroids, axis=1)   # YOUR CODE HERE\n",
        "    \n",
        "    return cluster_assignment"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s2xplv0BJoru"
      },
      "source": [
        "**Checkpoint**. For the last time, let us check if Step 1 was implemented correctly. With rows 0, 2, 4, and 6 of `tf_idf` as an initial set of centroids, we assign cluster labels to rows 0, 10, 20, ..., and 90 of `tf_idf`. The resulting cluster labels should be `[0, 1, 1, 0, 0, 2, 0, 2, 2, 1]`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "clGXQkLFJorv",
        "outputId": "507cd096-9ab9-4e5b-ab11-77c7028db1db"
      },
      "source": [
        "if np.allclose(assign_clusters(tf_idf[0:100:10], tf_idf[0:8:2]), np.array([0, 1, 1, 0, 0, 2, 0, 2, 2, 1])):\n",
        "    print('Pass')\n",
        "else:\n",
        "    print('Check your code again.')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Pass\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6w4sbVtWJorw"
      },
      "source": [
        "### Revising clusters"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SRS8-AQHJorx"
      },
      "source": [
        "Let's turn to Step 2, where we compute the new centroids given the cluster assignments. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VQJJh0WHJorx"
      },
      "source": [
        "SciPy and NumPy arrays allow for filtering via Boolean masks. For instance, we filter all data points that are assigned to cluster 0 by writing\n",
        "```\n",
        "data[cluster_assignment==0,:]\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y4kZVsuFJory"
      },
      "source": [
        "To develop intuition about filtering, let's look at a toy example consisting of 3 data points and 2 clusters."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lu50VLk1Jory"
      },
      "source": [
        "data = np.array([[1., 2., 0.],\n",
        "                 [0., 0., 0.],\n",
        "                 [2., 2., 0.]])\n",
        "centroids = np.array([[0.5, 0.5, 0.],\n",
        "                      [0., -0.5, 0.]])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d7azESN0Jorz"
      },
      "source": [
        "Let's assign these data points to the closest centroid."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4tWMrzZxJor0",
        "outputId": "87f84d21-7ccc-45f5-d2b4-a9f43534cca5"
      },
      "source": [
        "cluster_assignment = assign_clusters(data, centroids)\n",
        "print(cluster_assignment)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[0 1 0]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9B9I6RvaJor1"
      },
      "source": [
        "The expression `cluster_assignment==1` gives a list of Booleans that says whether each data point is assigned to cluster 1 or not:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EnidkQf6Jor2",
        "outputId": "eefaa037-abe3-4e44-9968-21141f6e4afe"
      },
      "source": [
        "cluster_assignment==1"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([False,  True, False])"
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m1J_uSxOJor2"
      },
      "source": [
        "Likewise for cluster 0:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UzpzopcOJor3",
        "outputId": "138a0237-2edd-4da9-c8c1-a5bc3f1eb5da"
      },
      "source": [
        "cluster_assignment==0"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([ True, False,  True])"
            ]
          },
          "metadata": {},
          "execution_count": 24
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "el1qcF1oJor4"
      },
      "source": [
        "In lieu of indices, we can put in the list of Booleans to pick and choose rows. Only the rows that correspond to a `True` entry will be retained.\n",
        "\n",
        "First, let's look at the data points (i.e., their values) assigned to cluster 1:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2BwOSGjmJor4",
        "outputId": "75a3d630-55b5-4a8c-c994-234c3e9aea9b"
      },
      "source": [
        "data[cluster_assignment==1]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[0., 0., 0.]])"
            ]
          },
          "metadata": {},
          "execution_count": 25
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kz1v-9KwJor5"
      },
      "source": [
        "This makes sense since [0 0 0] is closer to [0 -0.5 0] than to [0.5 0.5 0].\n",
        "\n",
        "Now let's look at the data points assigned to cluster 0:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "w8w0NCOSJor6",
        "outputId": "b68e3f0e-474b-4894-f235-d5de256513c9"
      },
      "source": [
        "data[cluster_assignment==0]"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[1., 2., 0.],\n",
              "       [2., 2., 0.]])"
            ]
          },
          "metadata": {},
          "execution_count": 26
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LkXjtCEmJor6"
      },
      "source": [
        "Again, this makes sense since these values are each closer to [0.5 0.5 0] than to [0 -0.5 0].\n",
        "\n",
        "Given all the data points in a cluster, it only remains to compute the mean. Use [np.mean()](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.mean.html). By default, the function averages all elements in a 2D array. To compute row-wise or column-wise means, add the `axis` argument. See the linked documentation for details. \n",
        "\n",
        "Use this function to average the data points in cluster 0:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "H4cVJe5SJor7",
        "outputId": "37dab491-a32b-45bf-a7f9-e94fe211030d"
      },
      "source": [
        "data[cluster_assignment==0].mean(axis=0)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([1.5, 2. , 0. ])"
            ]
          },
          "metadata": {},
          "execution_count": 27
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x9brzkopJor7"
      },
      "source": [
        "We are now ready to complete this function:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kJjZlOuVJor8"
      },
      "source": [
        "def revise_centroids(data, k, cluster_assignment):\n",
        "    new_centroids = []\n",
        "    for i in range(k):\n",
        "        # Select all data points that belong to cluster i. Fill in the blank (RHS only)\n",
        "        member_data_points = data[cluster_assignment==i]   # YOUR CODE HERE\n",
        "        # Compute the mean of the data points. Fill in the blank (RHS only)\n",
        "        centroid = member_data_points.mean(axis=0)   # YOUR CODE HERE\n",
        "        \n",
        "        # Convert numpy.matrix type to numpy.ndarray type\n",
        "        centroid = centroid.A1\n",
        "        new_centroids.append(centroid)\n",
        "    new_centroids = np.array(new_centroids)\n",
        "    \n",
        "    return new_centroids"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tyk0sACqJor9"
      },
      "source": [
        "**Checkpoint**. Let's check our Step 2 implementation. Letting rows 0, 10, ..., 90 of `tf_idf` as the data points and the cluster labels `[0, 1, 1, 0, 0, 2, 0, 2, 2, 1]`, we compute the next set of centroids. Each centroid is given by the average of all member data points in corresponding cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KU2KwGxPJor9",
        "outputId": "b451f0c5-1877-46c1-900a-c2b37dc30451"
      },
      "source": [
        "result = revise_centroids(tf_idf[0:100:10], 3, np.array([0, 1, 1, 0, 0, 2, 0, 2, 2, 1]))\n",
        "if np.allclose(result[0], np.mean(tf_idf[[0,30,40,60]].toarray(), axis=0)) and \\\n",
        "   np.allclose(result[1], np.mean(tf_idf[[10,20,90]].toarray(), axis=0))   and \\\n",
        "   np.allclose(result[2], np.mean(tf_idf[[50,70,80]].toarray(), axis=0)):\n",
        "    print('Pass')\n",
        "else:\n",
        "    print('Check your code')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Pass\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pQeLfkxAJor_"
      },
      "source": [
        "### Assessing convergence"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UnK2mquSJor_"
      },
      "source": [
        "How can we tell if the k-means algorithm is converging? We can look at the cluster assignments and see if they stabilize over time. In fact, we'll be running the algorithm until the cluster assignments stop changing at all. To be extra safe, and to assess the clustering performance, we'll be looking at an additional criteria: the sum of all squared distances between data points and centroids. This is defined as\n",
        "$$\n",
        "J(\\mathcal{Z},\\mu) = \\sum_{j=1}^k \\sum_{i:z_i = j} \\|\\mathbf{x}_i - \\mu_j\\|^2.\n",
        "$$\n",
        "The smaller the distances, the more homogeneous the clusters are. In other words, we'd like to have \"tight\" clusters."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tloFkaAAJosA"
      },
      "source": [
        "def compute_heterogeneity(data, k, centroids, cluster_assignment):\n",
        "    \n",
        "    heterogeneity = 0.0\n",
        "    for i in range(k):\n",
        "        \n",
        "        # Select all data points that belong to cluster i. Fill in the blank (RHS only)\n",
        "        member_data_points = data[cluster_assignment==i, :]\n",
        "        \n",
        "        if member_data_points.shape[0] > 0: # check if i-th cluster is non-empty\n",
        "            # Compute distances from centroid to data points (RHS only)\n",
        "            distances = pairwise_distances(member_data_points, [centroids[i]], metric='euclidean')\n",
        "            squared_distances = distances**2\n",
        "            heterogeneity += np.sum(squared_distances)\n",
        "        \n",
        "    return heterogeneity"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4BnBkqR9JosB"
      },
      "source": [
        "Let's compute the cluster heterogeneity for the 2-cluster example we've been considering based on our current cluster assignments and centroids."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wGlLGjGjJosB",
        "outputId": "ed461ff0-a37b-4e59-e0aa-32c5b00795b2"
      },
      "source": [
        "compute_heterogeneity(data, 2, centroids, cluster_assignment)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "7.25"
            ]
          },
          "metadata": {},
          "execution_count": 31
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B0RzW85RJosB"
      },
      "source": [
        "### Combining into a single function"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9vRnFPm8JosC"
      },
      "source": [
        "Once the two k-means steps have been implemented, as well as our heterogeneity metric we wish to monitor, it is only a matter of putting these functions together to write a k-means algorithm that\n",
        "\n",
        "* Repeatedly performs Steps 1 and 2\n",
        "* Tracks convergence metrics\n",
        "* Stops if either no assignment changed or we reach a certain number of iterations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sBZK_lEsJosC"
      },
      "source": [
        "# Fill in the blanks\n",
        "def kmeans(data, k, initial_centroids, maxiter, record_heterogeneity=None, verbose=False):\n",
        "    '''This function runs k-means on given data and initial set of centroids.\n",
        "       maxiter: maximum number of iterations to run.\n",
        "       record_heterogeneity: (optional) a list, to store the history of heterogeneity as function of iterations\n",
        "                             if None, do not store the history.\n",
        "       verbose: if True, print how many data points changed their cluster labels in each iteration'''\n",
        "    centroids = initial_centroids[:]\n",
        "    prev_cluster_assignment = None\n",
        "    \n",
        "    for itr in range(maxiter):        \n",
        "        if verbose:\n",
        "            print(itr)\n",
        "        \n",
        "        # 1. Make cluster assignments using nearest centroids\n",
        "        # YOUR CODE HERE\n",
        "        cluster_assignment = assign_clusters(data, centroids)\n",
        "            \n",
        "        # 2. Compute a new centroid for each of the k clusters, averaging all data points assigned to that cluster.\n",
        "        # YOUR CODE HERE\n",
        "        centroids = revise_centroids(data, k, cluster_assignment)\n",
        "            \n",
        "        # Check for convergence: if none of the assignments changed, stop\n",
        "        if prev_cluster_assignment is not None and \\\n",
        "          (prev_cluster_assignment==cluster_assignment).all():\n",
        "            break\n",
        "        \n",
        "        # Print number of new assignments \n",
        "        if prev_cluster_assignment is not None:\n",
        "            num_changed = np.sum(prev_cluster_assignment!=cluster_assignment)\n",
        "            if verbose:\n",
        "                print('    {0:5d} elements changed their cluster assignment.'.format(num_changed))   \n",
        "        \n",
        "        # Record heterogeneity convergence metric\n",
        "        if record_heterogeneity is not None:\n",
        "            # YOUR CODE HERE\n",
        "            score = compute_heterogeneity(data, k, centroids, cluster_assignment)\n",
        "            record_heterogeneity.append(score)\n",
        "        \n",
        "        prev_cluster_assignment = cluster_assignment[:]\n",
        "        \n",
        "    return centroids, cluster_assignment"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YE189-aQJosD"
      },
      "source": [
        "## Plotting convergence metric"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5IeBkCHtJosE"
      },
      "source": [
        "We can use the above function to plot the convergence metric across iterations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OHYhH5czJosF"
      },
      "source": [
        "def plot_heterogeneity(heterogeneity, k):\n",
        "    plt.figure(figsize=(7,4))\n",
        "    plt.plot(heterogeneity, linewidth=4)\n",
        "    plt.xlabel('# Iterations')\n",
        "    plt.ylabel('Heterogeneity')\n",
        "    plt.title('Heterogeneity of clustering over time, K={0:d}'.format(k))\n",
        "    plt.rcParams.update({'font.size': 16})\n",
        "    plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sfy9FehDJosG"
      },
      "source": [
        "Let's consider running k-means with K=3 clusters for a maximum of 400 iterations, recording cluster heterogeneity at every step.  Then, let's plot the heterogeneity over iterations using the plotting function above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "R3Hf4g67JosH",
        "outputId": "b12e4ae7-f010-46aa-c5f6-ce1afb3e4d84"
      },
      "source": [
        "k = 3\n",
        "heterogeneity = []\n",
        "initial_centroids = get_initial_centroids(tf_idf, k, seed=0)\n",
        "centroids, cluster_assignment = kmeans(tf_idf, k, initial_centroids, maxiter=400,\n",
        "                                       record_heterogeneity=heterogeneity, verbose=True)\n",
        "plot_heterogeneity(heterogeneity, k)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0\n",
            "1\n",
            "    19157 elements changed their cluster assignment.\n",
            "2\n",
            "     7739 elements changed their cluster assignment.\n",
            "3\n",
            "     5119 elements changed their cluster assignment.\n",
            "4\n",
            "     3370 elements changed their cluster assignment.\n",
            "5\n",
            "     2811 elements changed their cluster assignment.\n",
            "6\n",
            "     3233 elements changed their cluster assignment.\n",
            "7\n",
            "     3815 elements changed their cluster assignment.\n",
            "8\n",
            "     3172 elements changed their cluster assignment.\n",
            "9\n",
            "     1149 elements changed their cluster assignment.\n",
            "10\n",
            "      498 elements changed their cluster assignment.\n",
            "11\n",
            "      265 elements changed their cluster assignment.\n",
            "12\n",
            "      149 elements changed their cluster assignment.\n",
            "13\n",
            "      100 elements changed their cluster assignment.\n",
            "14\n",
            "       76 elements changed their cluster assignment.\n",
            "15\n",
            "       67 elements changed their cluster assignment.\n",
            "16\n",
            "       51 elements changed their cluster assignment.\n",
            "17\n",
            "       47 elements changed their cluster assignment.\n",
            "18\n",
            "       40 elements changed their cluster assignment.\n",
            "19\n",
            "       34 elements changed their cluster assignment.\n",
            "20\n",
            "       35 elements changed their cluster assignment.\n",
            "21\n",
            "       39 elements changed their cluster assignment.\n",
            "22\n",
            "       24 elements changed their cluster assignment.\n",
            "23\n",
            "       16 elements changed their cluster assignment.\n",
            "24\n",
            "       12 elements changed their cluster assignment.\n",
            "25\n",
            "       14 elements changed their cluster assignment.\n",
            "26\n",
            "       17 elements changed their cluster assignment.\n",
            "27\n",
            "       15 elements changed their cluster assignment.\n",
            "28\n",
            "       14 elements changed their cluster assignment.\n",
            "29\n",
            "       16 elements changed their cluster assignment.\n",
            "30\n",
            "       21 elements changed their cluster assignment.\n",
            "31\n",
            "       22 elements changed their cluster assignment.\n",
            "32\n",
            "       33 elements changed their cluster assignment.\n",
            "33\n",
            "       35 elements changed their cluster assignment.\n",
            "34\n",
            "       39 elements changed their cluster assignment.\n",
            "35\n",
            "       36 elements changed their cluster assignment.\n",
            "36\n",
            "       36 elements changed their cluster assignment.\n",
            "37\n",
            "       25 elements changed their cluster assignment.\n",
            "38\n",
            "       27 elements changed their cluster assignment.\n",
            "39\n",
            "       25 elements changed their cluster assignment.\n",
            "40\n",
            "       28 elements changed their cluster assignment.\n",
            "41\n",
            "       35 elements changed their cluster assignment.\n",
            "42\n",
            "       31 elements changed their cluster assignment.\n",
            "43\n",
            "       25 elements changed their cluster assignment.\n",
            "44\n",
            "       18 elements changed their cluster assignment.\n",
            "45\n",
            "       15 elements changed their cluster assignment.\n",
            "46\n",
            "       10 elements changed their cluster assignment.\n",
            "47\n",
            "        8 elements changed their cluster assignment.\n",
            "48\n",
            "        8 elements changed their cluster assignment.\n",
            "49\n",
            "        8 elements changed their cluster assignment.\n",
            "50\n",
            "        7 elements changed their cluster assignment.\n",
            "51\n",
            "        8 elements changed their cluster assignment.\n",
            "52\n",
            "        3 elements changed their cluster assignment.\n",
            "53\n",
            "        3 elements changed their cluster assignment.\n",
            "54\n",
            "        4 elements changed their cluster assignment.\n",
            "55\n",
            "        2 elements changed their cluster assignment.\n",
            "56\n",
            "        3 elements changed their cluster assignment.\n",
            "57\n",
            "        3 elements changed their cluster assignment.\n",
            "58\n",
            "        1 elements changed their cluster assignment.\n",
            "59\n",
            "        1 elements changed their cluster assignment.\n",
            "60\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAELCAYAAADqYO7XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5wddb3/8dd7dxNSSANCTQNpBgkBIggiTaQICl75of4AQxMLKHBVLNefgIJebIgdkKZeehEsSJMiLZBcOkgPhBAIhPS+2c/vj/lucvZsO3vOnp3snvfzwSEz32mfmVM++/3Od2YUEZiZmVl+6vIOwMzMrNY5GZuZmeXMydjMzCxnTsZmZmY5czI2MzPLmZOxmZlZzpyMzSok6UhJt/XQtj4o6QVJiyQd1oXljpF0XzVjK5WkpyXtnXcc1Zbeoy3yjsN6BydjK5mk6ZL2Kyor+Ude0pmS/lSd6PITEf8TEfs3j0sKSVtWaXPfA34VEetGxJ+rtI12dce+RcR2EXF3N4W0VpB0t6QTCsvSe/RyDrG0+J5K+rSkuZL26sI6NpB0v6Q5kuZJelDSB6sTsYGTsfUikhryjmEtMBZ4Ou8gytFX3j9J9XnHUCpJk4FfAwdHxD1dWHQRcBwwEhgBnAv8pa+8h2sjJ2PrVpI2lXS9pLclvSLpK6n8QODbwKdS893jqXyYpIslzZI0U9LZzT92qdZ9v6TzJM0Bzkzz/yGt/1VJ35FUl+avl/RTSe+kbZ+canINJW7rPkk/SbWIVyQdVLBfnS6bhu9Nizye9vNTkp6S9LGCdfVLMe7YzjH8nKQXJb0r6WZJm6byl4AtyH4UF0lap41lR0u6IR2fOZJ+1cY84wqPSypbXbOTtKWkeyTNT3Fe3d6+pfJDJD2WalAPSJpQsN7pkr4h6QlgsaSGwppbai25Jr2nC5U1YU8qWH4nSY+maddKulrS2e0ct7r0eXhV0uy0zmFp2i2STi6a/3FJ/5GGt5V0ezrmz0k6omC+yyT9VtLfJS0G9ilazznAh4BfpePyq1S+uhUhreM3KY5F6XO9saSfp8/bvws/D2rne9QVkj4P/BQ4ICIe6MqyEbEsIp6LiCZAwCqypLxeV+OwEkWEX36V9AKmA/sVlR0D3JeG64BpwHeB/mSJ42WyHwOAM4E/FS1/I3ABMBjYEHgY+HzBuhuBLwMNwEDgD8BNwBBgHPA8cHya/wvAM8Aosh+OO4AAGkrc1krgc0A98EXgDUAlLntfwT4FsGXB+OnA1QXjhwJPtnOM9wXeAXYC1gF+Cdzb0XtQMK0eeBw4L8U5ANijjfdpXOFxSWV3Ayek4SuB/0rv5+p1tLNvOwKzgV3T9ienGNcpiPcxYDQwsHgf0mdiGfDRtPwPgYfStP7Aq8ApQD/gP4AVwNnt7P9xwItkn7t1gRuAP6ZpnwXuL5h3PDAvHePBwAzgWLLP2Y7pPRif5r0MmA98sPmYtLHt1cevrWOV1vEOsHM6pv8EXklx1QNnA3eV8j0q8Xt6PfAWsEMb0+d18Ppm0bxPpGMewEV5/wb15VfuAfjVe17pS76o6Mu7hDU/8rsCrxUt8y3g0jR8JgXJGNgIWN78I53KPlPwo3RM4frSj9aK5h/JVPZ54O40/E9Sgkzj+6UfkYYSt/ViwbRBadmNS1y2o2S8KbAQGJrGrwNOb+cYXwz8qGB8XbI/EsYVvAftJePdgLcpSLIF01bHSOfJ+A/AhcCoNtZTvG+/Bb5fNM9zwF4F8R7XxueoMBnfUTBtPLA0De8JzCT9QZTK7qP9ZHwn8KWC8W3SsWsg++NtMTA2TTsHuCQNfwr4V9G6LgDOSMOXAX/o5Lux+vi1dazSOi4qmPZl4NmC8e2BeaV8j0r8ni4g+6O1rpRlOlnfALLP++RK1+VX+y83U1tXHRYRw5tfwJcKpo0FNk3NlfMkzSNrmt6onXWNJavxzCqY/wKymmezGQXDG6T5Xy0oexXYLA1vWjR/4XAp23qzeSAilqTBdUtctl0R8QZwP/BJScOBg4D/aWf2TQv3LyIWAXMK9rEjo4FXI6KxlLg6cDpZ0+TDqdn4uA7mHQt8teg9H022H81mtL3oam8WDC8BBqQm9E2BmZEyQgnranHs0nADsFFELAT+Bnw6TfsMa96DscCuRftwJNkfYqXuQyneKhhe2sb4ugXxdOV71JYvAlsDv5ekCmImsibrK4FvStqhknVZ+3wy3rrTDOCViNiqnenFjwibQVbj3KCDBFK4zDtkNZ2xZM3RAGPIak8As8iaqJuN7uK22lPJss0uB04g+849GBEz25nvDbL9A0DSYGB91uxjZ3GOkdTQSZyL07+DyGpQUJB4IuJNsuZ6JO0B3CHp3oh4sZ1tnhMR53SwvXIfDTcL2EySChLyaOClduZvcezIPhuNrEl6VwJnpHPfA4C7UvkM4J6I+EgHsXS2D935+LvOvkeleAv4MHAP8Buy5Axkl1x1sNwPIuIH7UzrR9Zk/ngFcVk7XDO27vQwsDB12BmorEPV+yS9P01/Cxin1OEqImYBtwE/lTQ0dcB5j9q5BCMiVgHXAOdIGiJpLPCfQPPlUtcAp0jaLNVAv1GwbJe2VbTdri77FtmPVqE/k50HPoWsGbg9VwLHSpqorIPWD4ApETG9szjJjv8s4L8lDZY0QG1cjhIRb5Ml96PSe3Qc8J7m6ZL+j6TmP2rmkiWapnb27SLgC5J2VWawpIMlDSkh3s48SNZx6GRlHb8OBXbpYP4rgdMkbS5pXbJjd3XBHyZ/J0vW30vlzfv0V2BrSUcr61zXT9L7Jb23C7G29Z6Xq8PvkaS9JXWa/FOLzIeBAyWdV1C+bgevH6RtfEDSHpL6pxi+QVYzn9JN+2hFnIyt26RkeQgwkaxzyjvA74FhaZZr079zJP1vGv4sWSeVZ8h++K8DNulgM18mq9m9THb+8ArgkjTtIrKk+QTwKNmPbyPZD3o52yrUlWXPBC5PTYxHAETEUrJONZuTdSxqU0TcAfy/NO8ssiT56fbmL1p2FfAxYEvgNeB1svOhbfkc8HWyJvDtgMLetu8HpqQa1M3AKbHmetkW+xYRU9O6fkV2XF4kOz9dsYhYQdZp63iy/glHkSXO5e0scgnwR+Bess/fMrLPS/P6lpMd+/3IPjfN5QuB/cmO8xtkzebnknXuKtX5wOGpZ/QvurBcKyV8j0bT8v3qaF2vkXUKPFzSD7sQxjpkl0TNIfvD7aNkl0e90YV1WBc09xQ163OUXZr0u4gY2+nMPUDSd4GtI+KovGPprSRNIXtPL807lrxI+j1wbUTcmncs1n18ztj6DEkDya4BvY2sSe0MskuScidpPbIa3tF5x9KbpFMBz5HVDo8EJgD/yDWonEXECZ3PZb2Nm6mtLxFwFllz6aPAs2TXauZK0ufIOuXcEhH3dja/tbANWYehecBXgcPTOXyzPsXN1GZmZjlzzdjMzCxnNXfOeIMNNohx48blHYaZmdWgadOmvRMRI4vLay4Zjxs3jqlTp+YdhpmZ1SBJr7ZV7mZqMzOznDkZm5mZ5czJ2MzMLGdOxmZmZjlzMjYzM8tZzfWmrtT8pSuZOXcp85auYP6SlYwaMYjtRw3rfEEzM7N2OBl30bVTZ3D2355dPX7M7uOcjM3MrCJupu6i4YP6txifv3RlTpGYmVlf4WTcRcMH9msxPnfJipwiMTOzvsLJuItGDG6ZjOctcc3YzMwq42TcRcMGupnazMy6l5NxFw0f5GZqMzPrXk7GXVR8znj+0pU0NfmZ0GZmVj4n4y5qqK9jyDprrgiLgIXLGnOMyMzMejsn4zIMc1O1mZl1IyfjMowoutZ4njtxmZlZBZyMy+BOXGZm1p2cjMswrLgTl681NjOzCjgZl6FVM7VrxmZmVgEn4zK0bqZ2zdjMzMpX1WQsabqkJyU9JmlqKpso6aHmMkm7pPJDJT1RUL5HwXomS3ohvSYXlO+c1v+ipF9IUjX3p1mrZmp34DIzswr0xCMU94mIdwrGfwScFRG3SPpoGt8buBO4OSJC0gTgGmBbSesBZwCTgACmSbo5IuYCvwU+B0wB/g4cCNxS7R1yM7WZmXWnPJqpAxiahocBbwBExKKIaL6V1eA0H8ABwO0R8W5KwLcDB0raBBgaEQ+l5f4AHNYTO+BmajMz607VrhkHcJukAC6IiAuBU4FbJf2E7I+B3ZtnlvQJ4IfAhsDBqXgzYEbBOl9PZZul4eLyViSdCJwIMGbMmIp3qjgZ+zpjMzOrRLVrxntExE7AQcBJkvYEvgicFhGjgdOAi5tnjogbI2Jbshru97sriIi4MCImRcSkkSNHVry+4UXN1PPdTG1mZhWoajKOiJnp39nAjcAuwGTghjTLtamseLl7gS0kbQDMBEYXTB6Vymam4eLyqit+WISbqc3MrBJVS8aSBksa0jwM7A88RXaOeK80277AC2meLZt7Q0vaCVgHmAPcCuwvaYSkEWk9t0bELGCBpA+k5T4L3FSt/SlU3Jt6wbKVrPKTm8zMrEzVPGe8EXBjyq8NwBUR8Q9Ji4DzJTUAy0jncoFPAp+VtBJYCnwqdcx6V9L3gUfSfN+LiHfT8JeAy4CBZL2oq96TGtKTmwY0rH5aU/bkppWtmq/NzMxKUbVkHBEvAzu0UX4fsHMb5ecC57azrkuAS9oonwq8r+JgyzB8UL8Wj06cu8TJ2MzMyuM7cJVp+EBfa2xmZt3DybhMvrzJzMy6i5NxmYqbpF0zNjOzcjkZl6n48qZ5vrzJzMzK5GRcphHFzdROxmZmViYn4zINczO1mZl1EyfjMrVqpnYHLjMzK5OTcZlGDHYztZmZdQ8n4zIN83XGZmbWTZyMy+TrjM3MrLs4GZdpRKsOXE7GZmZWHifjMg0d0PK23n5yk5mZlcvJuEzNT25qFgEL3FRtZmZlcDKuQKumaidjMzMrg5NxBYo7cc11j2ozMyuDk3EFhhXd+GO+O3GZmVkZnIwrUNxM7ZqxmZmVw8m4Aq2uNXbN2MzMyuBkXAHfn9rMzLqDk3EFhvvJTWZm1g2cjCvgZmozM+sOTsYV8P2pzcysOzgZV8DN1GZm1h2cjCvQqgOXm6nNzKwMTsYVcM3YzMy6g5NxBYrvwLVgWSONq5pyisbMzHorJ+MK1NepjUcpNuYUjZmZ9VZOxhVyU7WZmVXKybhCI1o9ucmduMzMrGucjCs0rKhmPH+pa8ZmZtY1TsYV8uVNZmZWKSfjCrmZ2szMKuVkXKFWzdTuwGVmZl3kZFwhP0bRzMwq5WRcoRGD3UxtZmaVcTKu0PCBvs7YzMwq42RcoWFFHbjmu5nazMy6yMm4QiOKOnDNdc3YzMy6yMm4Qr7O2MzMKuVkXKGhA/shrRlf6Cc3mZlZF5WUjCXdIOlgSU7eRbInN/m8sZmZla/U5Pob4P8CL0j6b0nbVDGmXmf4IF9rbGZm5SspGUfEHRFxJLATMB24Q9IDko6V1K/jpfs+nzc2M7NKlNzsLGl94BjgBOBR4Hyy5Hx7VSLrRfxMYzMzq0RDKTNJuhHYBvgj8LGImJUmXS1parWC6y1aNVO7ZmxmZl1QUjIGLoqIvxcWSFonIpZHxKQqxNWr+P7UZmZWiVKbqc9uo+zB7gykN3MztZmZVaLDmrGkjYHNgIGSdgSar6gdCgyqcmy9hpupzcysEp01Ux9A1mlrFPCzgvKFwLerFFOvU5yMfUtMMzPrig6TcURcDlwu6ZMRcX0PxdTrFDdT+6YfZmbWFZ01Ux8VEX8Cxkn6z+LpEfGzNharOb7O2MzMKtFZM/Xg9O+61Q6kN/OTm8zMrBKdNVNfkP49q2fC6Z2KzxnPd83YzMy6oNQHRWwt6U5JT6XxCZK+U93Qeo8hA4qe3LS8kZV+cpOZmZWo1OuMLwK+BawEiIgngE9XK6jepr5ODBvoJzeZmVl5Sk3GgyLi4aKyxu4OpjdzJy4zMytXqcn4HUnvAQJA0uHArI4XAUnTJT0p6bHme1hLmijpoeYySbuk8iMlPZHmf0DSDgXrOVDSc5JelPTNgvLNJU1J5VdL6t86ip4xrNXlTe7EZWZmpSk1GZ8EXABsK2kmcCrwxRKX3SciJhbcw/pHwFkRMRH4bhoHeAXYKyK2B74PXAggqR74NXAQMB74jKTxaZlzgfMiYktgLnB8iTF1uxHFN/5Y7JqxmZmVptTnGb8cEfsBI4FtI2KPiJhe5jaD7HaaAMOAN9I2HoiIuan8IbK7fgHsAryYYlgBXAUcKknAvsB1ab7LgcPKjKlifliEmZmVq9RHKK4DfBIYBzQodR2OiO91smgAt0kK4IKIuJCsVn2rpJ+Q/TGwexvLHQ/ckoY3A2YUTHsd2BVYH5gXEY0F5Zu1E/+JwIkAY8aM6STk8vhhEWZmVq5SH6F4EzAfmAYs78L694iImZI2BG6X9G/gcOC0iLhe0hHAxcB+zQtI2ocsGe/Rhe10KP0RcCHApEmTorvWW8gPizAzs3KVmoxHRcSBXV15RMxM/86WdCNZk/Nk4JQ0y7XA75vnlzQhjR8UEXNS8UxgdGEsqWwOMFxSQ6odN5fnonUztWvGZmZWmlI7cD0gafuurFjSYElDmoeB/YGnyM4R75Vm2xd4Ic0zBrgBODoini9Y1SPAVqnndH+y65tvjogA7iKraUOW5G/qSozdqXUztWvGZmZWmlJrxnsAx0h6hayZWkBExIQOltkIuDGdX24AroiIf0haBJwvqQFYRjqXS9azen3gN2mZxoiYFBGNkk4GbgXqgUsi4um0zDeAqySdDTxK1uSdCzdTm5lZuUpNxgd1dcUR8TKwQxvl9wE7t1F+AnBCO+v6O/D3draxS1djq4ZWNWM3U5uZWYlKvbTpVbLztvum4SWlLlsrfAcuMzMrV6kPijiDrEn4W6moH/CnagXVGxU/RtHJ2MzMSlVq7fYTwMeBxQAR8QYwpFpB9UZDBjS0eHLTIj+5yczMSlRqMl6Rei8335t6cPVC6p3q6tSqdvz63KU5RWNmZr1Jqcn4GkkXkF3X+zngDrLHKlqBbTdu2Vjw2Iy57cxpZma2RqkduH5Cdg/o64FtgO9GxC+rGVhvtOOY4S3GH3ttXk6RmJlZb1LqpU1ExO3A7VWMpdebOHpEi/FHZzgZm5lZ50rtTb1Q0oKi1wxJN0raotpB9hYTR7esGT/zxgKWrVyVUzRmZtZblHrO+OfA18meijQK+BpwBdnjDC+pTmi9z8gh6zBqxMDV441NwdNvzM8xIjMz6w1KTcYfj4gLImJhRCxIT0E6ICKuBkZ0tnAt2XFMUVO1zxubmVknSk3GSyQdIakuvY4gu680pMudLLNjUVO1zxubmVlnSk3GRwJHA7PT62jgKEkDgZOrFFuvNNE9qs3MrItK6k2dHsjwsXYm39d94fR+2206lP71daxId9+aOW8psxcsY8OhA3KOzMzM1lal9qYelXpOz06v6yWNqnZwvdE6DfWM33RoizI3VZuZWUdKbaa+FLgZ2DS9/pLKrA3Flzi5E5eZmXWk1GQ8MiIujYjG9LoMGFnFuHq1Vnfi8m0xzcysA6Um4zmSjpJUn15HAXOqGVhvtmPRnbieeH0+q5rc6dzMzNpWajI+DjgCeBOYBRwOHFutoHq70esNZP3Ba57gtGTFKp5/a2GOEZmZ2dqs02QsqR74QUR8PCJGRsSGEXFYRLzWA/H1SpJaNVX7vLGZmbWn02QcEauAsZL6dzavrdG6E5fPG5uZWdtKfWrTy8D9km4GFjcXRsTPqhJVH1B8W8zHfHmTmZm1o9Rk/FJ61QFDqhdO3zFh1DAkiNRv68W3F7Fg2UqGDuiXb2BmZrbWKfUOXGcBSBoUEUuqG1LfMGRAP7bacF2ef2sRkCXlJ2bMZ4+tNsg5MjMzW9uUegeu3SQ9A/w7je8g6TdVjawPKL7EyeeNzcysLV15nvEBpGuLI+JxYM9qBdVXFD80wrfFNDOztpSajImIGUVFq7o5lj6n9Z245hHhm3+YmVlLpSbjGZJ2B0JSP0lfA56tYlx9wlYbDmFw//rV4+8uXsFr7/qUu5mZtVRqMv4CcBKwGTATmAh8qVpB9RX1dWLCqNa1YzMzs0KlJuNtIuLIiNgo3YHrKOC91Qysr/CduMzMrDOlJuNfllhmRVrdics1YzMzK9LhdcaSdgN2B0ZK+s+CSUOB+raXskLFPaqfeWM+y1auYkA/Hz4zM8t0VjPuD6xLlrSHFLwWkD25yTqx4ZABjBoxcPX4ylXB028syDEiMzNb23RYM46Ie4B7JF0WEa/6DlzlmTh6OK/PXbp6/MGX3mHnsSM6WMLMzGpJqeeMN/UduMq3y+brtRi/8uEZNK5qyikaMzNb2/gOXD3g4ztsyoB+aw71zHlLuePZt3KMyMzM1ia+A1cPGD6oP5/YcbMWZZfePz2fYMzMbK3jO3D1kMm7j2sxPuWVd3nGHbnMzIzK7sB1UrWC6ou23Xgou22xfouyyx+Ynk8wZma2VikpGUfEO8V34IqIOdUOrq855oPjWoz/+bGZvLt4RT7BmJnZWqOzm378Emj3MUMR8ZVuj6gP2++9G7HZ8IHMnJdd5rS8sYmrHnmNL+29Zc6RmZlZnjqrGU8FpqXXxwuGm1/WBfV1YvLuY1uU/fHBV32Zk5lZjevsph+XNw9LOrVw3MrzqUljOO/2F1i6MuuMPmv+Mm575i0+uv0mOUdmZmZ5KfnSJjporrbSDRvUj0/s1PIyp8t8mZOZWU3rSjK2bjJ5t3Etxh+e/i5PzZyfTzBmZpa7DpOxpIWSFkhaAExoHm4u76EY+5xtNh7C7u/xZU5mZpbpMBlHxJCIGJpeDQXDQyJiaE8F2RcdU3QTkJsef4M5i5bnE4yZmeXKzdQ5+fB7N2rxaMUVjU1cMeW1HCMyM7O8OBnnpL5Orc4d/+KfL3DXv2fnE5CZmeXGyThHR0wazeD+9avHV64KPv+nadz/4js5RmVmZj3NyThHwwb145xPbN+ibEVjEydcPpVHpr+bU1RmZtbTnIxzdtiOm3H2Ye9rUbZ05SqOvfQRHp8xL6eozMysJzkZrwWO+sBYvnPwe1uULVreyGcvediPWTQzqwFOxmuJEz60BV8/YJsWZfOXruToi6fwwlsLc4rKzMx6gpPxWuSkfbbk5H1aPsFpzuIVfPrCh9ypy8ysD3MyXst8df+tOX6PzVuUzVm8gqMvnsIv73yBpibfItzMrK+pajKWNF3Sk5IekzQ1lU2U9FBzmaRdUvm2kh6UtFzS14rWc6Ck5yS9KOmbBeWbS5qSyq+W1L+a+9MTJPGdg9/LkbuOaVHeFPDT25/n2Mse4d3FK3KKzszMqqEnasb7RMTEiJiUxn8EnBURE4HvpnGAd4GvAD8pXFhSPfBr4CBgPPAZSePT5HOB8yJiS2AucHxV96SHSOL7h76PUz68FVLLafc8/zYH/+JfTHt1bj7BmZlZt8ujmTqA5vtaDwPeAIiI2RHxCLCyaP5dgBcj4uWIWAFcBRwqScC+wHVpvsuBw6odfE+pqxOnfWRrLj92F9Yb3LLCP2v+Mj51wYNcfN8rRLjZ2syst6t2Mg7gNknTJJ2Yyk4FfixpBlkt+FudrGMzYEbB+OupbH1gXkQ0FpW3IunE1CQ+9e233y5zV/Kx59Yj+dtX9mDnsSNalDc2Bd//6zOc9ZdnnJDNzHq5aifjPSJiJ7Im5pMk7Ql8ETgtIkYDpwEXVzkGIuLCiJgUEZNGjhxZ7c11u02GDeSqEz/A5z60eatplz0wnXP/8ZwTsplZL1bVZBwRM9O/s4EbyZqcJwM3pFmuTWUdmQmMLhgflcrmAMMlNRSV90n96uv4r4PHc8HROzNkQEOLab+75yV+ceeLOUVmZmaVqloyljRY0pDmYWB/4Cmyc8R7pdn2BV7oZFWPAFulntP9gU8DN0dWFbwLODzNNxm4qXv3Yu1zwHYbc+0XdmP4oH4tys+743kuuOelnKIyM7NKNHQ+S9k2Am7M+lnRAFwREf+QtAg4P9VolwEnAkjaGJhK1rmrSdKpwPiIWCDpZOBWoB64JCKeTtv4BnCVpLOBR+mBJu+1wbYbD+WPx+3K/73oIRYub1xd/sNb/s06DXUc88HWzdlmZrb2Uq2da5w0aVJMnTo17zC6xbRX3+Xoix9myYpVLcp/+B/b85ldxrSzlJmZ5UXStIJLfVfzHbh6sZ3HrsfFk9/POg0t38Zv3/gkf360z54+NzPrc5yMe7nd3rM+F312Ev3r17yVEXD69U/w+twlOUZmZmalcjLuA/bceiS/OXInGurW3K5rRWMTv77LPazNzHoDJ+M+Yr/xG3HGx7drUXbt1Nd5bY5rx2Zmazsn4z7kM+8fzdj1B60eb2wKfvHPzq4cMzOzvDkZ9yEN9XWc8uGtWpTd8L+v88o7i3OKyMzMSuFk3MccOnEzthg5ePV4U8D5dzyfY0RmZtYZJ+M+pr5OnLrf1i3Kbnr8DV6cvTCniMzMrDNOxn3QIdtvwtYbrbt6PALOu8Pnjs3M1lZOxn1QXZ04rah2/LcnZvHvNxfkFJGZmXXEybiPOmC7jRm/ydAWZefd7nPHZmZrIyfjPqquTpz2kZa141uffounZs7PKSIzM2uPk3Eftt97N2TCqGEtyn7untVmZmsdJ+M+TGpdO77j2dk8NmNeThGZmVlbnIz7uL23HsmOY4a3KLv8gen5BGNmZm1yMu7jJLW6K9ftz7zFspWr2lnCzMx6mpNxDfjQViMZOWSd1eOLljdy93Nv5xiRmZkVcjKuAfV14qPv27hF2V+feCOnaMzMrJiTcY04ZIdNW4zf+exslqxozCkaMzMr5GRcI3YeM4JNhg1YPb505SrufHZ2jhGZmVkzJ+MaUVcnDt5+kxZlbqo2M1s7OBnXkOKm6ruee5uFy1bmFI2ZmTVzMq4hO4waxuj1Bq4eX9HYxB3PvpVjRGZmBk7GNUUSh0xoWTv+y+OzcorGzMyaORnXmEMmtDxv/K8X3mb+EjdVm5nlycm4xozfZChbbDB49fjKVcGtT7+ZY0RmZuZkXGOypuqWteO/uFe1mVmunIxr0MeKelU/8NIc5ixanlM0ZmbmZFyDttpoCNtsNGT1+IphTDEAAAl9SURBVKqm4Jan3FRtZpYXJ+MaVdxU7RuAmJnlx8m4RhXfAGTKK+8ye8GynKIxM6ttTsY1avMNBrPdpkNXj0fA3570NcdmZnlwMq5hxR25fnv3S1zzyAxWrmrKKSIzs9rkZFzDih8cMXvhck6//gn2/vHd/OmhV1neuCqnyMzMaouTcQ0bvd4g9tx6ZKvymfOW8p0/P8VeP7qbS+9/heffWsjrc5cwd/EKljeuIiJyiNbMrO9Srf2wTpo0KaZOnZp3GGuN+UtW8sNbnuW6aa/T2FTaZ6GhTgzqX886/eppqBMN9aKhro6GOlGfxuslkKgT1KV/JSFAAiHSLKuHIZsHWD1f83DztDXDzdGoxXjr6dn6aaO8WHvL0NEy7a6r7SkdrKosHe1Pl9dV1va7vlS7S3R4nNs5nu0s01FU7S9Txr6U+Hkq3lKX11VyRJXpzs/T2qqc9xlg202GcOSuYyvfvjQtIiYVlzdUvGbr1YYN6sd/f3ICJ++7Jb+75yWueeR1VnRyzrixKViwrBGWNfZQlGZm+frI+I26JRm3x83UBsCoEYM4+7Dtuff0fTj2g+NYp8EfDTOznuKasbWw8bABnPGx7fjS3lty2QOv8MBLc1i4rJElyxtZvGIVi5c3ltycbWZmpXEytjaNHLIOXz9g2zanrWhsYvHyRlasaqKxKWhc/W/Q2NRE46oggKYIIoIIaIrstptBkP4jAoJI/7K6Y1is/l82neZ5V5elfwvnZ830NSWFZYWlrbWcr+3lWy3T7rp65o+V7txMdHh0um/77S3T8XvT9tR2l+nwPWtnXeXsS0fT2t3Prm+/x/70rYH+Q5Xs4egRg7otjrY4GVuX9W+oo39D/7zDMDPrM3xi0MzMLGdOxmZmZjlzMjYzM8uZk7GZmVnOnIzNzMxy5mRsZmaWs5q7N7Wkt4FXu2l1GwDvdNO6eiPvv/e/lvcffAy8/13f/7ER0eoJPTWXjLuTpKlt3fC7Vnj/vf+1vP/gY+D97779dzO1mZlZzpyMzczMcuZkXJkL8w4gZ97/2lbr+w8+Bt7/buJzxmZmZjlzzdjMzCxnTsZmZmY5czIug6QDJT0n6UVJ38w7np4g6RJJsyU9VVC2nqTbJb2Q/h2RZ4zVJGm0pLskPSPpaUmnpPKaOAaSBkh6WNLjaf/PSuWbS5qSvgtXS+rTz9aUVC/pUUl/TeM1s/+Spkt6UtJjkqamspr4/DeTNFzSdZL+LelZSbt11zFwMu4iSfXAr4GDgPHAZySNzzeqHnEZcGBR2TeBOyNiK+DONN5XNQJfjYjxwAeAk9L7XivHYDmwb0TsAEwEDpT0AeBc4LyI2BKYCxyfY4w94RTg2YLxWtv/fSJiYsG1tbXy+W92PvCPiNgW2IHss9Atx8DJuOt2AV6MiJcjYgVwFXBozjFVXUTcC7xbVHwocHkavhw4rEeD6kERMSsi/jcNLyT7Em5GjRyDyCxKo/3SK4B9getSeZ/dfwBJo4CDgd+ncVFD+9+Omvj8A0gaBuwJXAwQESsiYh7ddAycjLtuM2BGwfjrqawWbRQRs9Lwm8BGeQbTUySNA3YEplBDxyA10T4GzAZuB14C5kVEY5qlr38Xfg6cDjSl8fWprf0P4DZJ0ySdmMpq5vMPbA68DVyaTlX8XtJguukYOBlbt4jsGrk+f52cpHWB64FTI2JB4bS+fgwiYlVETARGkbUQbZtzSD1G0iHA7IiYlncsOdojInYiO0V3kqQ9Cyf29c8/0ADsBPw2InYEFlPUJF3JMXAy7rqZwOiC8VGprBa9JWkTgPTv7JzjqSpJ/cgS8f9ExA2puKaOAUBqmrsL2A0YLqkhTerL34UPAh+XNJ3s1NS+ZOcPa2X/iYiZ6d/ZwI1kf5DV0uf/deD1iJiSxq8jS87dcgycjLvuEWCr1IuyP/Bp4OacY8rLzcDkNDwZuCnHWKoqnR+8GHg2In5WMKkmjoGkkZKGp+GBwEfIzpvfBRyeZuuz+x8R34qIURExjuw7/8+IOJIa2X9JgyUNaR4G9geeokY+/wAR8SYwQ9I2qejDwDN00zHwHbjKIOmjZOeP6oFLIuKcnEOqOklXAnuTPTLsLeAM4M/ANcAYssdSHhERxZ28+gRJewD/Ap5kzTnDb5OdN+7zx0DSBLLOKfVkf8RfExHfk7QFWU1xPeBR4KiIWJ5fpNUnaW/gaxFxSK3sf9rPG9NoA3BFRJwjaX1q4PPfTNJEsg58/YGXgWNJ3wcqPAZOxmZmZjlzM7WZmVnOnIzNzMxy5mRsZmaWMydjMzOznDkZm5mZ5czJ2KwPkPRDSftIOkzSt9qZ50xJX0vDx0jatBu3v7ek3QvGvyDps921frO+zsnYrG/YFXgI2Au4t4T5jwG6lIwL7jTVlr2B1ck4In4XEX/oyvrNapmvMzbrxST9GDiA7Cb2LwHvAV4BrouI7xXNeyawCJhO9kjMmcBSsttajgd+BqwLvAMcExGzJN0NPAbsAVwJPA98h+ymB3OAI4GBZH8IrCK7kf6Xye5OtCgifpJulPA7YFCK8biImJvWPQXYBxgOHB8R/5K0HXBp2kYd8MmIeKGbDpnZWsk1Y7NeLCK+TvYM3cuA9wNPRMSE4kRctMx1wFTgyPTgh0bgl8DhEbEzcAlQeFe5/hExKSJ+CtwHfCDdKP8q4PSImE6WbM9Lz7r9V9Em/wB8IyImkN3B7IyCaQ0RsQtwakH5F4DzU2yTyO4JbNanddTsZGa9w07A42RPUXq2k3nbsg3wPuD27Bbc1AOzCqZfXTA8Crg63RC/P1ktvF3pGbDDI+KeVHQ5cG3BLM0P3JgGjEvDDwL/lZ4ffINrxVYLnIzNeqnU/HsZWYJ8h6wZWOmZw7tFxNJSVwU8HRG7tTN9ccHwL4GfRcTN6R7NZ5YReqHm+zivIv0eRcQVkqYABwN/l/T5iPhnhdsxW6u5mdqsl4qIx1JT7vNk53z/CRyQmoo7S8QLgSFp+DlgpKTdIHtUZDpv25ZhrHlM4OSC8sL1FcY4H5gr6UOp6GjgnuL5CqWHErwcEb8gewLOhE72xazXczI268UkjQTmRkQTsG1EPFPiopcBv0u16HqyxwCeK+lxsg5bu7ez3JnAtZKmkdXGm/0F+ISkxwoSb7PJwI8lPQFMBNo9n50cATyVYnsf2Tlnsz7NvanNzMxy5pqxmZlZzpyMzczMcuZkbGZmljMnYzMzs5w5GZuZmeXMydjMzCxnTsZmZmY5+/+wGVtdYm0l8AAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 504x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "z69fzxAlNYLc",
        "outputId": "670e8201-1eac-461a-ae32-7e60db20587e"
      },
      "source": [
        "print(heterogeneity)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[58324.51355736017, 58184.88702561689, 58138.8351095433, 58115.5443794838, 58105.144360097955, 58097.525184447324, 58087.152190572655, 58067.523328335345, 58038.163658158126, 58033.277654160556, 58032.80538698504, 58032.63119637119, 58032.554427394636, 58032.51498266415, 58032.48655252265, 58032.4628390735, 58032.44795796612, 58032.43449158974, 58032.42602967142, 58032.41911185921, 58032.41058106258, 58032.4034507918, 58032.39943856925, 58032.39710978097, 58032.39550452751, 58032.39362039754, 58032.391293004, 58032.38939093164, 58032.38753814069, 58032.3854990268, 58032.3821917849, 58032.378699442255, 58032.37264035171, 58032.365997387846, 58032.359590041546, 58032.3537596153, 58032.34721772186, 58032.342702778304, 58032.338476927755, 58032.33478185108, 58032.330043465045, 58032.322585288915, 58032.31617340006, 58032.31114087295, 58032.30835655161, 58032.30617328345, 58032.304782515435, 58032.30379008189, 58032.30282286614, 58032.301868322225, 58032.30089023452, 58032.29970841954, 58032.299265305046, 58032.298842522, 58032.29839474874, 58032.298170224196, 58032.29785759554, 58032.297509798605, 58032.29722383869, 58032.29707355844]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bPL9Wdb8JosI"
      },
      "source": [
        "**Quiz Question**. (True/False) The clustering objective (heterogeneity) is non-increasing for this example."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AfqDUfXGJosI"
      },
      "source": [
        "**Quiz Question**. Let's step back from this particular example. If the clustering objective (heterogeneity) would ever increase when running k-means, that would indicate: (choose one)\n",
        "\n",
        "1. k-means algorithm got stuck in a bad local minimum\n",
        "2. There is a bug in the k-means code\n",
        "3. All data points consist of exact duplicates\n",
        "4. Nothing is wrong. The objective should generally go down sooner or later."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9FL8vullJosJ"
      },
      "source": [
        "**Quiz Question**. Which of the cluster contains the greatest number of data points in the end? Hint: Use [`np.bincount()`](http://docs.scipy.org/doc/numpy-1.11.0/reference/generated/numpy.bincount.html) to count occurrences of each cluster label.\n",
        " 1. Cluster #0\n",
        " 2. Cluster #1\n",
        " 3. Cluster #2"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "W7kGbLvwJosJ",
        "outputId": "68c404ce-7299-474e-d2d7-af4c6d6e63e4"
      },
      "source": [
        "print (np.bincount(cluster_assignment))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[19595 10427 29049]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "syHuVtq6JosK"
      },
      "source": [
        "## Beware of local maxima"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M8KsG_nvJosK"
      },
      "source": [
        "One weakness of k-means is that it tends to get stuck in a local minimum. To see this, let us run k-means multiple times, with different initial centroids created using different random seeds.\n",
        "\n",
        "**Note:** Again, in practice, you should set different seeds for every run. We give you a list of seeds for this assignment so that everyone gets the same answer.\n",
        "\n",
        "This may take several minutes to run."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": false,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "z2bzIudKJosK",
        "outputId": "ee62cf48-4703-4d36-f800-fac6c24b5420"
      },
      "source": [
        "k = 10\n",
        "heterogeneity = {}\n",
        "cluster_assignment_dict = {}\n",
        "import time\n",
        "start = time.time()\n",
        "for seed in [0, 20000, 40000, 60000, 80000, 100000, 120000]:\n",
        "    initial_centroids = get_initial_centroids(tf_idf, k, seed)\n",
        "    centroids, cluster_assignment = kmeans(tf_idf, k, initial_centroids, maxiter=400,\n",
        "                                           record_heterogeneity=None, verbose=False)\n",
        "    # To save time, compute heterogeneity only once in the end\n",
        "    heterogeneity[seed] = compute_heterogeneity(tf_idf, k, centroids, cluster_assignment)\n",
        "\n",
        "    # This is the line we added for the next quiz question\n",
        "    cluster_assignment_dict[seed] = np.bincount(cluster_assignment)\n",
        "    \n",
        "#    print('seed={0:06d}, heterogeneity={1:.5f}'.format(seed, heterogeneity[seed]))\n",
        "    # And this is the modified print statement\n",
        "    print('seed={0:06d}, heterogeneity={1:.5f}, cluster_distribution={2}'.format(seed, heterogeneity[seed], \n",
        "                                           cluster_assignment_dict[seed]))\n",
        "    sys.stdout.flush()\n",
        "end = time.time()\n",
        "print(end-start)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "seed=000000, heterogeneity=57457.52442, cluster_distribution=[18047  3824  5671  6983  1492  1730  3882  3449  7139  6854]\n",
            "seed=020000, heterogeneity=57533.20100, cluster_distribution=[ 3142   768  3566  2277 15779  7278  6146  7964  6666  5485]\n",
            "seed=040000, heterogeneity=57512.69257, cluster_distribution=[ 5551  6623   186  2999  8487  3893  6807  2921  3472 18132]\n",
            "seed=060000, heterogeneity=57466.97925, cluster_distribution=[ 3014  3089  6681  3856  8080  7222  3424   424  5381 17900]\n",
            "seed=080000, heterogeneity=57494.92990, cluster_distribution=[17582  1785  7215  3314  6285   809  5930  6791  5536  3824]\n",
            "seed=100000, heterogeneity=57484.42210, cluster_distribution=[ 6618  1337  6191  2890 16969  4983  5242  3892  5562  5387]\n",
            "seed=120000, heterogeneity=57554.62410, cluster_distribution=[ 6118  5841  4964  8423  4302  3183 16481  1608  5524  2627]\n",
            "142.62409734725952\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MmWL4JuKJosL"
      },
      "source": [
        "Notice the variation in heterogeneity for different initializations. This indicates that k-means sometimes gets stuck at a bad local minimum."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R3iRSRDHJosL"
      },
      "source": [
        "**Quiz Question**. Another way to capture the effect of changing initialization is to look at the distribution of cluster assignments. Add a line to the code above to compute the size (# of member data points) of clusters for each run of k-means. Look at the size of the largest cluster (most # of member data points) across multiple runs, with seeds 0, 20000, ..., 120000. How much does this measure vary across the runs? What is the minimum and maximum values this quantity takes?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tLv3U3Z0ORxz",
        "outputId": "763d7042-63fe-4c59-9ebc-a172b4e9d59b"
      },
      "source": [
        "for k, v in cluster_assignment_dict.items():\n",
        "  print (k, np.max(v))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0 18047\n",
            "20000 15779\n",
            "40000 18132\n",
            "60000 17900\n",
            "80000 17582\n",
            "100000 16969\n",
            "120000 16481\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KSWMuZIyJosM"
      },
      "source": [
        "One effective way to counter this tendency is to use **k-means++** to provide a smart initialization. This method tries to spread out the initial set of centroids so that they are not too close together. It is known to improve the quality of local optima and lower average runtime."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Vcfp2jouJosN"
      },
      "source": [
        "def smart_initialize(data, k, seed=None):\n",
        "    '''Use k-means++ to initialize a good set of centroids'''\n",
        "    if seed is not None: # useful for obtaining consistent results\n",
        "        np.random.seed(seed)\n",
        "    centroids = np.zeros((k, data.shape[1]))\n",
        "    \n",
        "    # Randomly choose the first centroid.\n",
        "    # Since we have no prior knowledge, choose uniformly at random\n",
        "    idx = np.random.randint(data.shape[0])\n",
        "    centroids[0] = data[idx,:].toarray()\n",
        "    # Compute distances from the first centroid chosen to all the other data points\n",
        "    squared_distances = pairwise_distances(data, centroids[0:1], metric='euclidean').flatten()**2\n",
        "    \n",
        "    for i in range(1, k):\n",
        "        # Choose the next centroid randomly, so that the probability for each data point to be chosen\n",
        "        # is directly proportional to its squared distance from the nearest centroid.\n",
        "        # Roughtly speaking, a new centroid should be as far as from ohter centroids as possible.\n",
        "        idx = np.random.choice(data.shape[0], 1, p=squared_distances/sum(squared_distances))\n",
        "        centroids[i] = data[idx,:].toarray()\n",
        "        # Now compute distances from the centroids to all data points\n",
        "        squared_distances = np.min(pairwise_distances(data, centroids[0:i+1], metric='euclidean')**2,axis=1)\n",
        "    \n",
        "    return centroids"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9xQs38qGJosO"
      },
      "source": [
        "Let's now rerun k-means with 10 clusters using the same set of seeds, but always using k-means++ to initialize the algorithm.\n",
        "\n",
        "This may take several minutes to run."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-2TDv4YpJosO",
        "outputId": "56985067-9e66-4bf0-f6ae-2d7e648f671c"
      },
      "source": [
        "%%time \n",
        "\n",
        "k = 10\n",
        "heterogeneity_smart = {}\n",
        "seeds = [0, 20000, 40000, 60000, 80000, 100000, 120000]\n",
        "for seed in seeds:\n",
        "    initial_centroids = smart_initialize(tf_idf, k, seed)\n",
        "    centroids, cluster_assignment = kmeans(tf_idf, k, initial_centroids, maxiter=400,\n",
        "                                           record_heterogeneity=None, verbose=False)\n",
        "    # To save time, compute heterogeneity only once in the end\n",
        "    heterogeneity_smart[seed] = compute_heterogeneity(tf_idf, k, centroids, cluster_assignment)\n",
        "    print('seed={0:06d}, heterogeneity={1:.5f}'.format(seed, heterogeneity_smart[seed]))\n",
        "    sys.stdout.flush()\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "seed=000000, heterogeneity=57468.63808\n",
            "seed=020000, heterogeneity=57486.94263\n",
            "seed=040000, heterogeneity=57454.35926\n",
            "seed=060000, heterogeneity=57530.43659\n",
            "seed=080000, heterogeneity=57454.51852\n",
            "seed=100000, heterogeneity=57471.56674\n",
            "seed=120000, heterogeneity=57523.28839\n",
            "CPU times: user 2min 55s, sys: 2.19 s, total: 2min 57s\n",
            "Wall time: 2min 58s\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GkaKTqTTJosP"
      },
      "source": [
        "Let's compare the set of cluster heterogeneities we got from our 7 restarts of k-means using random initialization compared to the 7 restarts of k-means using k-means++ as a smart initialization.\n",
        "\n",
        "The following code produces a [box plot](http://matplotlib.org/api/pyplot_api.html) for each of these methods, indicating the spread of values produced by each method."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 356
        },
        "id": "3vOCHPabJosQ",
        "outputId": "7948e3c7-396a-4b70-8de4-731009fc3e04"
      },
      "source": [
        "plt.figure(figsize=(8,5))\n",
        "plt.boxplot([list(heterogeneity.values()), list(heterogeneity_smart.values())], vert=False)\n",
        "plt.yticks([1, 2], ['k-means', 'k-means++'])\n",
        "plt.rcParams.update({'font.size': 16})\n",
        "plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAFTCAYAAAD4N0wZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYJklEQVR4nO3dfbAlZX0n8O8PEFij4CC4upYjEqJZ0UowIxtiQNTsgopGXbfU+E4SkN2sqbJ8w1UZ3yCWFTUbkhWWIAi6qahZIQpRzDBYSCwdRAQUdsmOxpf4goy8RAODPvtH95Xj4dyBuTNzz3NnPp+qUz23++k+z/lNn3O+t/vpvtVaCwBAr3abdwcAALZEWAEAuiasAABdE1YAgK4JKwBA1/aYdwdYmv33378deOCB8+4GAGw3V1xxxY2ttQOm5wsrK9SBBx6YDRs2zLsbALDdVNXXZ813GggA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga8IKANA1YQUA6JqwAgB0TVgBALomrAAAXRNWAICuCSsAQNeEFQCga3vMuwPsuvbbb79s2rRp3t1YVDt5n9Rbbpl3N+jIqlWrctNNN827G7DLEVaYm02bNqW1Nu9uLG7tvn33j2VXVfPuAuySnAYCALomrAAAXRNWAICuCSsAQNfuVVipqrVV1arKgNydlIGDwErks2vX4MgKANA1YWU7qaqvVdXarVznZVXl2lgA2IIlh5WqOqaqbquq06pq5nYWvoyr6jeq6q+q6taq+m5VnTSxjSur6p+r6gtV9WsztvGcqvpcVf2oqn5YVR+uqtVTbZ5fVeuq6vtjn66sqpfO2FarqrdX1SurauPYn0ur6pCpdkdX1eVVdfO4veur6s1LrRUAsHRLCitV9ZIkFyT5o9baH7TWfnoPq5yT5Ookz07ysSSnVNU7k7wryTuTPC/JLyT5WFXtOfE8r0jy0SRfSfLcJCckeUySS6vq/hPbPyjJR5K8MMmzkvxNkjPH9ae9KMnTk/xhkpcnWZ3k/IXxOFV10PjaNo79emaSd4/9AwCW2VYPmK2q1yZ5R5ITW2tn3svVzm2tvW1cf32G0PKqJI9srW0c5++W5Pwkh2cII/fLEGTe31o7buL5P5/k+iS/m+S9SdJaO2Vi+W5J1id5SJITk7xvqi+bkxzbWts8tk+SDyc5LMnlSR6XZM/x9S3ca33dVA0qye4zXuduU4OQW2vtJxPr7Z5kcjTYbuP86f+Hn7QZt06tquOTHJ8kq1evnl68zQxUg3vmfQLLb2vDynuS/F6S57bWzl+YOeNLePrL9qKFf7TW7qyqG5LsuxBURteN04eN08OT7JPkg1Nf5t8Y2x6ZMaxU1S8lees478G564jR7TNew8ULQWV09ThdnSGsfClDoPnLqjoryWdaa9+b2sYTk1wyY9tvGh8LLk1y1MTPfzeuO23z1M9PyhC4fk5r7YwkZyTJmjVrtvtYl+W+tbwPfVYif4KhLz5Hdg1bG1ZekOSaJJ+emv8PSR4+8fPLk5w98fP0X6u7Y5F5SbL3OH3QOJ1+rp/b5ngE5uIkP0ry+rEvd2Q4qnLcjPWm/wrZQqDZO0laazdU1dFJXpfk3CR7jUdzXtdau3Rse0WSx09t54IkH88YJka3TrU5Icnk6atjk5w8Y1vXz+g3AOyStjasPCXJp5JcVFVPa63dNs5/RpK9JtptvNuaW+8H4/RlSa6dsXwhCByeISgd0Vq7bGHhttwTprV2SZJLqmqvJE/IcNTmE1V1YGvtxtbarUk2TK5TVXck+XZrbcPdt/iz7f5cCKmqx4zzF10HAHZ1W/uFfm2G0xrrMgSWp7bWbmutXb3l1Zbk8gyB5ODW2jlbaHffcfqzUylVtSrJb29rB1prtydZNx69OT/JI5LcuK3bBQDuva0++tBa+2pVHZVhzMYnq+qY8UjDdtVau6WqXpPkz6rqgAzjXm5O8tAM4z7Wt9Y+lCHU3DK2OznDVTtvzBAq9t3a5x2vIDoyyYUZxsfsn+SkJN/OcAoMAFhGS7p0eTyd8cQMp18+VVX7bNde3fU8p2e4dPhRGcaPXJhkbYaQ9aWxzfczXF20e4bLl09NcmaS85b4tFdlCDynZjjldVqG01pPbq39eInbBACWqIxsX5nWrFnTNmxY2UNdqqrvKyvW7pusvXnevaAj3e+zsMJV1RWttTXT891uHwDomrACAHRNWAEAuiasAABdW/KN02B76PlW2e3kfbruH8tv1apV8+4C7JKEFeZmJVxV0dbOuwcAOA0EAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0LU95t0B2Fntt99+2bRp07y7sWK0k/dJveWWeXeDJVq1alVuuummeXeDnZSwAjvIpk2b0lqbdzdWjrX7qtcKVlXz7gI7MaeBAICuCSsAQNeEFQCga8IKS+YcNcCubbm+B7oLK1W1tqpaVRn8CwD0F1YAACYJKwBA11ZEWKmqY6rqtqo6rapm9rmqXjaePvqNqvqrqrq1qr5bVSdNbOPKqvrnqvpCVf3ajG08p6o+V1U/qqofVtWHq2r1VJvnV9W6qvr+2Kcrq+qlM7bVqurtVfXKqto49ufSqjpkqt3RVXV5Vd08bu/6qnrztlUMAHYe3YeVqnpJkguS/FFr7Q9aaz+9h1XOSXJ1kmcn+ViSU6rqnUneleSdSZ6X5BeSfKyq9px4nlck+WiSryR5bpITkjwmyaVVdf+J7R+U5CNJXpjkWUn+JsmZ4/rTXpTk6Un+MMnLk6xOcv7CeJyqOmh8bRvHfj0zybvH/gEA6fwOtlX12iTvSHJia+3Me7naua21t43rr88QWl6V5JGttY3j/N2SnJ/k8Axh5H4Zgsz7W2vHTTz/55Ncn+R3k7w3SVprp0ws3y3J+iQPSXJikvdN9WVzkmNba5vH9kny4SSHJbk8yeOS7Dm+voX7jK/bQj2OT3J8kqxevXqxZsvKFUHAAp8H7Cg9h5X3JPm9JM9trZ2/MLOqdk8y+Y74Sfv5e3RftPCP1tqdVXVDkn0XgsrounH6sHF6eJJ9knxw6iqkb4xtj8wYVqrql5K8dZz34Nx1dOr2Ga/h4oWgMrp6nK7OEFa+lCHQ/GVVnZXkM621783YzsLrOSPJGUmyZs2aLu5L7vboi/PBza7G58GuZ5e9dHnCC5Jck+TTU/P/IcMX/MJjerzI9F+Ou2OReUmy9zh90Dj99NS2Nyd5bJIHJsl4BObiJL+S5PVJjkjy+CRnJdlrxmuY/qteC4Fm7yRprd2Q5OgM/w/nJvnOOGbmiTO2BQC7pJ6PrDwlyaeSXFRVT2ut3TbOf0Z+PhhsvNuaW+8H4/RlSa6dsfzWcXp4kocnOaK1dtnCwm25J0xr7ZIkl1TVXkmekOGozSeq6sDW2o1L3S4A7Cx6DivXJjkqwxiOi6rqqa2121prV295tSW5PEMgObi1ds4W2t13nP7s1E5VrUry29vagdba7UnWjUdvzk/yiCTCCgC7vJ7DSlprX62qo5JckuSTVXVMa+3We1htKc9zS1W9JsmfVdUBGca93JzkoUmemGR9a+1DGULNLWO7kzNctfPGDKFi36193vEKoiOTXJhhfMz+SU5K8u0Mp8AAYJfX85iVJElr7foMgeHhST5VVfvsoOc5PcOlw4/KMH7kwiRrMwS6L41tvp/h6qLdM1y+fGqSM5Oct8SnvSpD4Dk1wymv0zKc1npya+3HS9wmAOxUyujtlWnNmjVtw4YN8+4GW1BVro7YGmv3TdbePO9esET2d7aHqrqitbZmen73R1YAgF2bsAIAdE1YAQC61vXVQLDSuYvtvddO3ke9VrBVq1bNuwvsxIQV2EEMNtx6be28ewD0yGkgAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDomrACAHRNWAEAuiasAABdE1YAgK5Va23efWAJqur7Sb4+737cS/snuXHenVih1G7bqN/Sqd22Ub+leXhr7YDpmcIKO1xVbWitrZl3P1Yitds26rd0ardt1G/7choIAOiasAIAdE1YYTmcMe8OrGBqt23Ub+nUbtuo33ZkzAoA0DVHVgCArgkrAEDXhJVdXFUdVVVtxuOHE23OXqRNq6rrtrDt149tLltk+UOr6qyq+k5V3V5VG6vq1Bntfr+qrhvbXF9Vr9g+r37bzat+VfXAqvqTqvp/VfXjsXanVdXd709Q9ayqurKq/qWqvl5Vb6yq3bdfFZZue9dvC+1+dardblV1UlV9bazLVVX1HxfpY5f73zxqV1WPHPe7L1fVbVX1T1V1QVX9yiJ97LJ2yfz2val1nj+2+eYiy7ut33LbY94doBuvTPKFiZ/vnPj325K8b6r9gUn+V5ILZm2sqg5K8sYk31tk+YFJPptk4/jc3x23efBUu99PcnqSU5N8OslTkvx5VVVr7X/c04taRstWv6qqcb1HJnlzkq8meXSStyZZU1WHt3EwWlUdneSjSf4iyauSHJrklCT3T/K6rXmBO9j2rN/ZGfaZSf9n6ue3JXl1kv+W5Iokz0/y4ao6trV24UKjFbL/LWft/kOSJyU5J8kXkzwgyWuTfK6qfrO1dsVCwxVSu2T5970kSVU9IMl7k3xnkeUrpX7Lo7XmsQs/khyVpCX5ra1c703jeocssvyTGd5o65NcNmP53yb5fJL7bOE59sjwZX3O1PyzMtwZctF1d+b6ZQgpLcnxU/NfMc5/1MS8K5NcOtXuzUnuSPLgna1+47y338O6D0pye5K3TM3/uyRfXin735xqt3/GCzMm5u2bZFOSD6yU2s2rflPtzxjf52cn+ebUsu7rt9wPp4FYqpckuaK1du30gqr6nSSPS3LSrBWr6heTHJ3kT1trm7fwHIcnOSDJeVPzz03ywCS/uYR+92LJ9Uuy5zi9ZWr+wuHr3cbtPCzJr2Z2/e6T5Klb3+1uLFq/e+HoDDWcrst5SR5bVY8Yf95Z978l1661dmMbvzUn5t2c4ejBQydm76y1S7Zt30uSVNUTkrwoyX9ZpMnOXL8lEVZY8MGq+klV/aCqPlRVqxdrOL7RDs5wKHh62aok70ny2tbaTYts4gnj9MdVdfF4PnZTVX2gqh440e6QcXrN1PoLHxKPvqcXtYyWs37XJvlMkjdV1Zqqul9VHZbhiMlFrbWvju1m1q+1tjHJj7IT1m904rhP/aiq1lXVEVPLD8lwZOWGqfnT+9VK2f+Ws3aztrlfksdkOB25YKXULlnm+lXVfTIcVXlXa216H1ywkuq3LIxZ4eYkf5zk0gy/qR+a5A1J/r6qDm2tzRpz8pIkmzOct532rgy/ZZ29hef8N+P0rAy/KZya4QPg1CSPrqrDWms/TbLf2G7T1PoLX+L7Zf6WvX6ttVZVT8tQu8lz7Z9I8p8mfl6sfgvzdsb6nZfk40m+neThSV6TZF1V/fvW2vqxzX5Jfjh9hCB336963//mUbtZ/jRJZRh/saD32iXzq9/rkuyV4fNuMSuhfstr3uehPPp7ZDgFcWdmnH9NsneGN9Bfz1h2RIaxEI+ZmLc+dx9z8YYM53cvmJr/vHH+U6fa7T3Vbo9x/pvmXat51G+c/6EMH4onJDlynH4nQ2DZbWzzO2OdfnnG+t9M8hfzrtX2rN8i27p/hr9OftnEvDOSfGdG24PHer14pe5/O7p2M9qcNNbiuKn5K652y7TvHZzkx0mOmZh3du4+ZmVF1m9HPpwG4m5aa1/M8Nv942csfmaGKwBmHQY9PcNVJ9+sqgeMo933SLL7+PNeY7sfjNOLp9b/1Dg9dJwu/Faxaqrdwm8Vi50mmasdXb+qenqSF2T4Uj29tfaZ1trpSV6c5GlJnjFub7H6Lczb2eo3a1u3Zghwk9valOQB41VVk6b3qxW3/y1D7X5mvIz2lCRvbK2dNbV4xdUuWZb6/fck6zJcPbXwHt8zw0V+D6iqfzW2W5H125GEFbZk1t9ieGmG0egXzlj2bzNckbJp4vGEJL8+/vvEsd09DUz76VS7Q6aWL5yv/co9bGfedlT9HjtOvzC1/ucntpMsUr/xsvH7Zuer373d1rUZDsP/4lSb6f1qJe9/O6p2SZKqenGSP0/yx621d8xYZyXXLtlx9Xt0hl8oJt/jL8hwanxT7jo1tNLrt90JK9xNVa1J8qjc9eW3MP9fZ7iS4kNt9lU8T5rxuCrDILEnJfnI2O5zGU5ZHD21/jHjdOFL+O8zfDi8cKrdizL8ZvHZrXldy2UZ6rdwX4bDptb/d+P0W0nSWvvHcf1Z9duc5KKteV3LZRvqN2tb+yQ5dmpbf5vh9c+qyzVtGICcrMD9bxlql6p6dpL3JzmztfbqRVZfcbVLlqV+z8/d3+OfzFCrJyU5bWy3Iuu3Ixlgu4urqg9muDHbFzNc+npohvPQ38pwyHLSC5PsnkUOg7YZg/BquBvkHpPLWmt3VtXrk5xdVe9L8tcZzuW+I8MYjXVju81V9aYMN0L6VoYbIz05yXFJ/mtr7Y4lvejtaB71y1CvdyT5QFW9Lcl1SX45yclJvpHkf0+0fUOSj1fV6RkGBR6a4WZzf9Jam3kzquW0PetXVa/O8EVzSe4a5PjqJA/OxId+a+17VfXuJCdV1a3jcz8vw771zIl2Xe9/86hdVR2ZYT+6KsP799cnNnN7a+3KpP/aJXPb9z43Y92XZajd+ol23ddv2c170IzHfB8Z3pxfzjAyfnOGL7szkjxkRturkly9ldtfn0UG6GUYY3FNhstI/ynDVQX3m9HuhAznkW9P8n+T/Od5123e9UvysAzjWzYm+Zdx+j+TPHRG2+eMz317kn/McInz7vOu3fauX4axOp/N8Bvp5gxjoy5IctiMtrtnCG1fH+vy5STPXWS7Xe5/86hdkrUZTmvMenxtpdRunvvejHXPztQA25VQv+V+1FgQAIAuGbMCAHRNWAEAuiasAABdE1YAgK4JKwBA14QVAKBrwgoA0DVhBQDo2v8HpJqmbBPJHPYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 576x360 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AAdidDwJJosQ"
      },
      "source": [
        "A few things to notice from the box plot:\n",
        "* On average, k-means++ produces a better clustering than Random initialization.\n",
        "* Variation in clustering quality is smaller for k-means++."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bvt8QBfKJosR"
      },
      "source": [
        "**In general, you should run k-means at least a few times with different initializations and then return the run resulting in the lowest heterogeneity.** Let us write a function that runs k-means multiple times and picks the best run that minimizes heterogeneity. The function accepts an optional list of seed values to be used for the multiple runs; if no such list is provided, the current UTC time is used as seed values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IzGszq63JosR"
      },
      "source": [
        "def kmeans_multiple_runs(data, k, maxiter, num_runs, seed_list=None, verbose=False):\n",
        "    heterogeneity = {}\n",
        "    \n",
        "    min_heterogeneity_achieved = float('inf')\n",
        "    best_seed = None\n",
        "    final_centroids = None\n",
        "    final_cluster_assignment = None\n",
        "    \n",
        "    for i in range(num_runs):\n",
        "        \n",
        "        # Use UTC time if no seeds are provided \n",
        "        if seed_list is not None: \n",
        "            seed = seed_list[i]\n",
        "            np.random.seed(seed)\n",
        "        else: \n",
        "            seed = int(time.time())\n",
        "            np.random.seed(seed)\n",
        "        \n",
        "        # Use k-means++ initialization\n",
        "        # YOUR CODE HERE\n",
        "        initial_centroids = smart_initialize(data, k, seed=seed)\n",
        "        \n",
        "        # Run k-means\n",
        "        # YOUR CODE HERE\n",
        "        centroids, cluster_assignment = kmeans(data, k, initial_centroids, maxiter=400, record_heterogeneity=None, verbose=False)\n",
        "        \n",
        "        # To save time, compute heterogeneity only once in the end\n",
        "        # YOUR CODE HERE\n",
        "        heterogeneity[seed] = compute_heterogeneity(data, k, centroids, cluster_assignment)\n",
        "        \n",
        "        if verbose:\n",
        "            print('seed={0:06d}, heterogeneity={1:.5f}'.format(seed, heterogeneity[seed]))\n",
        "            sys.stdout.flush()\n",
        "        \n",
        "        # if current measurement of heterogeneity is lower than previously seen,\n",
        "        # update the minimum record of heterogeneity.\n",
        "        if heterogeneity[seed] < min_heterogeneity_achieved:\n",
        "            min_heterogeneity_achieved = heterogeneity[seed]\n",
        "            best_seed = seed\n",
        "            final_centroids = centroids\n",
        "            final_cluster_assignment = cluster_assignment\n",
        "    \n",
        "    # Return the centroids and cluster assignments that minimize heterogeneity.\n",
        "    return final_centroids, final_cluster_assignment"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H_dVA9oMJosS"
      },
      "source": [
        "## How to choose K"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OKaNH-iEJosT"
      },
      "source": [
        "Since we are measuring the tightness of the clusters, a higher value of K reduces the possible heterogeneity metric by definition.  For example, if we have N data points and set K=N clusters, then we could have 0 cluster heterogeneity by setting the N centroids equal to the values of the N data points. (Note: Not all runs for larger K will result in lower heterogeneity than a single run with smaller K due to local optima.)  Let's explore this general trend for ourselves by performing the following analysis."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iOll_iHQJosT"
      },
      "source": [
        "Use the `kmeans_multiple_runs` function to run k-means with five different values of K.  For each K, use k-means++ and multiple runs to pick the best solution.  In what follows, we consider K=2,10,25,50,100 and 7 restarts for each setting.\n",
        "\n",
        "**IMPORTANT: The code block below will take about 10 minutes to finish**\n",
        "\n",
        "In order to speed up the computations, we run them with only one random seed, but for better performance, one should use more seeds and compare the results. If you don't mind running the code for approximately one hour, feel free to uncomment the following line of code below:\n",
        "\n",
        "`seed_list = [0]#, 20000, 40000, 60000, 80000, 100000, 120000]`\n",
        "\n",
        "Side note: In practice, a good implementation of k-means would utilize parallelism to run multiple runs of k-means at once. For an example, see [scikit-learn's KMeans](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": false,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DIvP2YWKJosU",
        "outputId": "ef0fa0a4-18c5-46e6-dec0-bf2aa4b6e7bc"
      },
      "source": [
        "%%time\n",
        "import numpy as np \n",
        "\n",
        "def plot_k_vs_heterogeneity(k_values, heterogeneity_values):\n",
        "    plt.figure(figsize=(7,4))\n",
        "    plt.plot(k_values, heterogeneity_values, linewidth=4)\n",
        "    plt.xlabel('K')\n",
        "    plt.ylabel('Heterogeneity')\n",
        "    plt.title('K vs. Heterogeneity')\n",
        "    plt.rcParams.update({'font.size': 16})\n",
        "    plt.tight_layout()\n",
        "\n",
        "centroids = {}\n",
        "cluster_assignment = {}\n",
        "heterogeneity_values = []\n",
        "k_list = [2, 10, 25, 50, 100]\n",
        "seed_list = [0]\n",
        "# Uncomment the following line to run the plot with all the seeds (it may take about an hour to finish).\n",
        "#seed_list = [0, 20000, 40000, 60000, 80000, 100000, 120000]\n",
        "\n",
        "for k in k_list:\n",
        "    heterogeneity = []\n",
        "    centroids[k], cluster_assignment[k] = kmeans_multiple_runs(tf_idf, k, maxiter=400,\n",
        "                                                               num_runs=len(seed_list),                                                               seed_list=seed_list,\n",
        "                                                               verbose=True)\n",
        "    score = compute_heterogeneity(tf_idf, k, centroids[k], cluster_assignment[k])\n",
        "    heterogeneity_values.append(score)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "seed=000000, heterogeneity=58224.59913\n",
            "seed=000000, heterogeneity=57468.63808\n",
            "seed=000000, heterogeneity=56913.24052\n",
            "seed=000000, heterogeneity=56399.72145\n",
            "seed=000000, heterogeneity=55649.66538\n",
            "CPU times: user 7min 6s, sys: 5.25 s, total: 7min 12s\n",
            "Wall time: 7min 11s\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 284
        },
        "id": "A7gGekYKJosV",
        "outputId": "1200102b-5c03-448b-e184-7c17bb95110f"
      },
      "source": [
        "plot_k_vs_heterogeneity(k_list, heterogeneity_values)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAELCAYAAADqYO7XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3gc5bXH8e+RZLl3yca94IJtukV1oxNKKAGbEBJaAoGEUFLuDaEGwg0hCSQECHAJkBsIxUBopppiGzAYuQA24N6rZLlXlXP/mJFZrVZlrZV2pf19nmcfWe+8M3t2bOtoZt7zvubuiIiISPJkJDsAERGRdKdkLCIikmRKxiIiIkmmZCwiIpJkSsYiIiJJpmQsIiKSZErGIiJJYmZ9zczN7PFkxyLJpWQsaSPiB9+rVWz/fbh9lpl1bej49oaZvR/GnFNNHzez2Yl4n7ocQ2rPzJaY2ZJkxyENJyvZAYgkm5kZcD9wJTAVONXdNyY3KkkTK4EhwKZkByLJpWQsac3MsoDHgQuAicBZ7r4tqUFJ2nD3YuDrZMchyafb1JK2zKw58DxBIv4PcHptErGZPRre+j2siu1/DbcfF9E21symmFmBme00s5Vm9pqZnZiozxMvM2tuZr8ys8/MbLuZbTKziWY2JqqfA2PK/xzxejyq37HhZ1offsYvzezX4S88kf0uDve/2MzONLOPzGyrmc2K6HOQmT0fnq9dZjbfzO4wszYxPkczM7vZzBaH7/uVmV1pZseE73NrjH0OMbPxZrY2PP7C8DFFm6h+e45hZnlm9raZbQnP1X/MrG8V57a2x6/wzLj8e6AP0CfqfN8anmM3s/ureN+jw+0Px9ouqUtXxpKWwh+KLwLHA/8HXOrupbXc/UngEoIk/mnUcTOB84AVwPth20+B+4CFwDPAFqA7MAI4GXi7bp8mfmbWAngLGEXwGR4GWgNnAu+Y2Th3fyHs/lvgYoIE8duIw0Qmz6uAe4FC4CVgAzAS+D1wOPCdGGGcB5wAvAxMIfx5FP4y8DqQCTxLcCv3WOA3wMlmNsrdd0Qc55/A+cDcMIb2wJ3hMWN99rOBp4HdBP8G1gCHAr8GjjWz0e6+O2q3w4D/At4DHgIOAc4CDjCz/d19Zx2PX24jwTm+Nvz+LxHb3gcmAfOB75nZLyLfN/TD8OsjVRxfUpW766VXWryAvoADHxA8G3bgb4DFeZwMggSxBsiM2nZyeNy7ItpmhP1bxThW5zp+pvfL3w+4tYqXA7Oj9vt92P7rqPZcYAlQALSMfp8qYhgGFAMfA+0j2o3glxAHzo1ovzhsKwWOiTpWJsEvLWXAmKhjPR7ud0tE+4lh21QgO6J9MLA93HZrRHsOsBlYBHSPeu9fhf1/GdF2TNjmwHlR/f8vbP9uHY5f/m/y8ai+S4AlVZzv/wr3+V5UexuCX/Q+T/b/Nb3ifyU9AL30aqhXxA++8te7dTjWn8JjnBzVXv4D+qCIthnAYqB5PXym96M+U1Wv2RH7ZBBcuc6p4phXhfucHv0+VfS/N+x/eIxt7cLE+lxEW3kyfi5G/zHhtpdibOsO7AIWRbSVJ+hTYvT/e4xk/POwbVyM/hnAOiA/oq08GU+K0b9825/rcPy9ScZdCK66J0a1/zA81jX1+f9Ir/p56Ta1pKPZBFcwx4a3+v68F8d4AvgFwa3qNwHMrBVwNkHi+yyi7zMEt01nm9nTBIltqrtv3/uPUEmuuxfG2hCjJGkw0AFYFut5KjAw/LofELMMLMoRhMnbzE6NsX1HeKxo+THaDgq/Tore4O6rzGw+MMzM2rr7loj+U2Mc6yPgihixAowws6Ex9imuItbpMdpWhF87JOD4tebu68zsReBcM+vn7ovDTT8k+GXlX3U5viSHkrGko6XAWILnf38yM3f3u+M5gLvPMrM5wNlm1ipMrGcQ3Cp8Iqr7XUARQWK4MXztMrPngF+4+9q6fZy4dQq/Hhi+qtI6juMZcFOcx4r1udtVsw2CRwPDwn5bgLZAsccuRVtXRawAV1cdakybY7SVhF8zE3D8eD1M8G/4EuBmMxsCHAU84+5F9fzeUg80mlrSkrt/TTAoaA3wZzP7+V4c5kmC5Htm+P0FBFeI/456L3f3/3X34QS3GMcSDJ66gGCAUkMrTyzPuLtV8/pttUepeDwHWldzrH4x9os1iUh5bFVNutI1qt8WoJmZdYjRt0s1xx9S3Wev4r1ro76PX+4dgmfrF5lZBhq41egpGUvaipGQr4vzEE8SJJTvm1lngsFbk919eTXvWeDuz7n7GQSjkUebWfu9+wR77SuCJJYXjv6ujVLYM1o82jSCK+PDExBb+Qjt0dEbzKwbMIjgmfGWsLn8ccCRMY51VIy2adX0T4REHb+UilfcFbi7EyTe3sCpwA8IxiW8U8f3lSRRMpa0FpWQ744nIbv7MoLymZOAnwLNqHyLurxUJ7qtJUEJTglhogvb9zWz/cysWZwfpdbcvQR4ENgX+H2sBGtmR4TPwMuV3/rsFeOQDxB8hr+ZWfcYx+oa3katjQ8IRiKfYWYjo7bdAWQTDJIr91T49SYzy454z4HARTGO/xiwleBzD44Ra3szO6SWscaSqOMXATlhCVp171VMUGrVBXgsTNLSCOmZsaQ9d//azI4leIZ8d/gM+S817Rd6kuAq7kaCwTPPxejzkpltBD4heF7dAjgF6Afc5+5bI/q+Q1DP249gRG19uRnIIyi3OcPMphAkgJ7AcIJBXt0IyoMA3gXOBZ43s9eBncBn7v6Ku39hZj8jKGOaZ2YTwtg7AgMIaplvIrgir5a7l5nZpcAbwEQzK68zPobganM6wTP48v5vhn3GAZ+Z2SsEz5O/S/D3eRrBaO7y/uvM7AKCQXVfmNlrBHW7rQnO+TEEdcvRA79qJYHHf5fg7+f18O9mN8Fdl8kR77U2/LzfCT/jY3sTs6SIZA/n1kuvhnrxTRnJq1Vs3w9YHfa5tpbH7EiQmBx4voo+VxJMbLEk7FtAcAX4A6JqnMM+DvSt5fu/H/bPqaZPpTrjsD2L4Ir+Y4JnnTsIrkpfBC4EsqL6/oHgl4liYpfjHAWMD8/hboK7DR8TJP7eEf0uDve/uJqYDwZeANaHx1oA/A/QJkbfbIJ66qUEvxB9DfwEOCd8n+ti7DOUoCxqeXj8QoIStDuB/SL6HUNUeVSMf0+P1+H4MY9BMBbhYWAVwd2TqmI4I9z2WrL/f+lVt5eFf6EiIk2Kmd1OcMfiNHd/Ldnx1Aczuwm4DTjHv5kxTRohJWMRadTMbB93XxPVNpjgsYAB3TyxNd0pIXymPz/8to8HYwGkkdIzYxFp7G60YMGNDwhuB/cjuH3bHLisqSXicGDbGIJxB92Bq5SIGz8lYxFp7F4nWBP42wTP8LcBHwJ3u/uEZAZWT04AbiGY1OROgmk/pZHTbWoREZEk05VxguXk5Hjfvn2THYaIiKSg6dOnF7p7bnS7knGC9e3bl/z8WPPfi4hIujOzpbHaNQOXiIhIkikZi4iIJJmSsYiISJIpGYuIiCSZknEKWrVxB//6OOYzfhERaYI0mjqFbN9dwoOTFvHw5IXsLC5jWPd2HNq7Y7LDEhGReqYr4xTykydncO8789lZHKz4dtsrX6JJWUREmj4l4xRy+aj+Fb6ftXwjL3+2KknRiIhIQ1EyTiFHD8jhxKFdK7Td+frX7NhdmqSIRESkISgZp5jfnDqEZpm25/vVm3by8ORFSYxIRETqm5JxiumX05qLj+5boe3BSQtZvWlHcgISEZF6p2Scgq46biCdWmfv+X5HcSl/fGNuEiMSEZH6pGScgtq3bMbPTxxUoe2FmSuZtXxjkiISEZH6pGScor57WC8Gd21boe22V+ao1ElEpAlSMk5RWZkZ3Hj6kAptM5Zt5JXPVycpIhERqS9Kxils1MBcThjSpULbna99xc5ilTqJiDQlSsYp7jenDiEr45tSp1WbdvK/KnUSEWlSGjQZm9kxZuYxXhuj+g0zsxfMbJWZbTOzOWb2SzPLiuqXYWbXm9kSM9tpZp+Z2TlVvPdlZva1me0ys7lmdkUV/c4ys5nh8Zaa2Y1mlpm4sxCf/rltuCiq1OmB9xeydvPO5AQkIiIJl6wr46uBoyJeJ5RvMLPuwPtAf+Ba4NvAi8BdwB1Rx7kduBW4DzgF+BgYb2anRnYys8uAh4DngW8B44EHzOzKqH4nh30+DY/3V+BG4H/q9nHr5urjBtKxVbM93+8oLuUulTqJiDQZ1pCjc83sGOA94ER3n1hFn8sJEudgd58X0f40MMbdu4XfdwGWA3e6+y0R/d4Bct39wPD7LGAV8Lq7XxTR71HgDKCbuxeHbTOBze4+JqLfzQQJube7r6npM+bl5Xl+fn5tTkdc/jV1CTe9NKdC28tXjeDAnh0S/l4iIlI/zGy6u+dFt6fiM+Py2S42R7VvpGK8J4d9n4jq9wRwgJn1C78/CsiN0e9fQGdgJICZ9QIOrqJfM4Ir5aQ5//DeDOzSpkKbVnUSEWkakpWMnzSzUjNbb2b/NrPeEdvGA4XAfWbWz8zamdnZwA+AP0f0GwbsAhZEHbv88nFoRD+A2XvTz90XA9sj+iVFVmYGN51eMYT8pRuY8IVKnUREGruGTsabCBLqj4DjCJ75ngBMDW874+5rCa5mhwCLwn2eB/7g7ndFHKsTsNErXxoWRWyP/LphL/uVt3WK0Q4Et9bNLN/M8gsKCqrqVmejB+Vy3H4VS51+/9rXKnUSEWnkGjQZu/tMd/+lu7/i7pPc/S8EA6q6EgzqwsxygReAbcC5wLHA74Abzey/GzLe2nL3h909z93zcnNz6/W9okudVm7cwT8+WFyv7ykiIvUr6c+M3X0GMA84LGz6L6AvcLK7P+/u77v7zcAfgdvNLCfstwHoYGYWdcjyK9iiiH4AHfeyX3lbUYz2BjegSxt+cFSfCm33v7eAdSp1EhFptJKejCOU324+AFjg7tG3i6cRDKQaEH4/B2gO7BvVr/zB6pcR/eCbZ8Jx9TOzvkCriH5Jd83xA+kQUeq0fXcpf3xTpU4iIo1V0pOxmeUBgwmSLcAaYICZRV+hHhF+XRl+fQMoBi6I6vd9YHY48ApgKsGAsFj9ioAPAdx9GfBZFf2Kgddr/6nqV4dW2Vx3QsVVnZ6bsYIvVmxKUkQiIlIXWTV3SRwzexJYDMwgKFU6BLieIMHeG3Z7kCAhvmVmfwTWA8cAvwT+4+7LAdx9nZndDVxvZlvCY55HMDDsjPL3dPdiM7uJYJKPlcDEsM+lwM/cfXdEiL8BXjWzh4CnwvhuBP5amxrjhnTBEb154uOlzF+3FQB3uP3VL3nmx0dS+c69iIiksoa+Mp5NkCgfA94kmGHrBeAIdy8EcPePgVFAAcEMWC8D3wFuo/JV6w0Eg7uuCY83Ahjn7q9GdnL3B4ErgXFhv/OBq9z9/qh+rxEMGjsy7Hcdwexbv677R0+sYFWniqVO05YU8frslPqdQUREaqFBZ+BKB/U1A1dVLnlsGu/N/aacqmfHlkz8+RhaNEvadNoiIlKFxjQDl8ThhtOGkhlR6rRiww4e/VClTiIijYmScSM3oEsbfnBkVKnTuwtYt0WlTiIijYWScRNw7QkDad/ym1KnbbtL+fOb86rZQ0REUomScRMQlDoNrND27PTlzF6pUicRkcZAybiJuODIPuyb23rP9+WlThqgJyKS+pSMm4hmMUqdPllcxJtzVOokIpLqlIybkGMHd2HMoIoLVdzx2lfsKtGqTiIiqUzJuIm58bQhFUqdlhft4LEPlyQvIBERqZGScRMzsGtbvn9E7wpt9727gIItu5IUkYiI1ETJuAm69oRBtGvxzbTjW3eVcPfbWtVJRCRVxZWMzax9fQUiidOxdTbXRq3q9PSny5mzSqVOIiKpKN4r41Vm9g8zO6xeopGE+cFRfegfVep02ysqdRIRSUXxJuM/AicCH5vZTDO73Mza1ENcUkfNMjO48bQhFdqCUqe1SYpIRESqElcydvdbgb7A2cAq4AFgpZn93cwOTnh0UifHDu7CqIE5Fdr+R6VOIiIpJ+4BXO5e5u4vu/tpwL7AvQRrFE83s0/M7GIza57oQCV+ZsZNpw8lotKJZUXbeVylTiIiKaWuo6k3A0XAVsCA9sA/gAVmNrKOx5YEGNS1LRccUXFVp7+p1ElEJKXsVTI2sxFm9n/ASuC3wLvAQe6+HzAUWAQ8lLAopU6uO3EQbSuVOmlVJxGRVBFvadPPzOwLYDJwKPAroIe7X+nuXwC4+1zgFmC/RAcre6dT62yuOb7iqk7PfLqML1dtTlJEIiISKd4r4z8BXwLHufv+7n6/u2+J0W8+cFudo5OEufCovvTL+abUqUyrOomIpIx4k3Fvdz/P3SdV18ndV7r7b+sQlyRYdlYGN5xasdRp6qL1vP2lSp1ERJIt3mQ81cwOirXBzPY3s0UJiEnqyfFDujByQMVSJ63qJCKSfPEm475AVWVLLYA+VWyTFGBm3Hj6kAqlTkvXb+f/PlqavKBERGSvRlNX9ZAxD9hYh1ikAey3TzvOP7ziqk73vjOf9VtV6iQikiw1JmMzu87MlpnZMoJE/Er59xGvAuB+4I36Dljq7ucnDqJt829Knbao1ElEJKlqc2W8CHgnfBmQH/F9+et54DrgsvoJUxKpc5vmXB1V6vTUtGV8vUalTiIiyZBVUwd3fwl4CYJnjsBt7r64nuOSenbR0X158pOlLFm/Hfim1OmJHx5R/vcsIiINJN6FIi5RIm4asrMy+E1UqdOHC9bz+uw1SYpIRCR91XhlbGY3A4+4+6rwz9Vxd789MaFJfTtxaFeO3rczHy1cv6ftqn/P4CfHDODq4weSnVXXqctFRKQ2rKYZmMysDDjS3aeFf66Ou3tmwqJrhPLy8jw/Pz/ZYdTaV6s3c9q9UyiL+mcwtFs77jnvYAbv0zY5gYmINEFmNt3d86Lba7z0cfcMd58W8efqXmmdiBujId3a8auTK08j/uXqzXz7bx/w0KSFlEZnahERSSjdhxSuPGZfHrvkMLq0rTify+7SMn7/+td89+GpLAsHeomISOLFnYwtcIaZ/cnMHjOzPmH7GDPrnvgQpSEcO7gLb103mm8fVPmv8NMlG/jWXyfz1LRlWlhCRKQexLuEYkfgI+BFgpriC4HO4ebLgF8nNDppUB1aZfO38w/h3vMPoX3LZhW2bd9dyvUvfMGlj3/Kus07kxShiEjTFO+V8R+BXsAIgiQcWZA6ETg+QXFJEp1xUHfeum40YwblVtr23twCTvrLZCZ8vjoJkYmINE3xJuMzgRvcfSqV56heRpCopQno2q4Fj19yGHecvT8tm1Ucl7dxezE//fcMrnl6Jpu2FycpQhGRpiPeZNwGWFnFthZUvFKWRs7MuOCIPrx+zSiG9+lYaftLs1Zx0l8mMXleQRKiExFpOuJNxnOBk6rYNgb4om7hSCrqm9OaZ398FL8+ZT+yMyv+k1m7eRcXPjqNG1/8gu27S5IUoYhI4xZvMn4AuNbMbgDK1+HrYGaXAFcRrNwkTVBmhnHFmH156aoR7BdjIpAnPl7GqX+dwvSlRUmITkSkcatxBq5KO5jdCfyS4Ja0ETw7LgPucvcbEh5hI9PYZuDaG7tKSvnrxPk8OGlhpZm7MgyuGLMv154wSNNpiohEqWoGrriTcXiwPgS3q3OB9cDb7r6ozlE2AemQjMtNX1rEz5/9jKUxJgQZ0q0dd487iCHd2iUhMhGR1JTQZCxVS6dkDLB9dwn/89pXPPHxskrbsjMzuO7EQVw+uj+ZGRrbJyKS6CvjfQieGbeI3ubuk/cqwiYi3ZJxuUnzCviv5z5j7eZdlbbl9enIn8cdRJ/OrZMQmYhI6khIMjazHsC/CEZOwzelTB7+Was2pWkyBti0vZibX57NS7NWVdrWKjuTG04bwvcO742ZrpJFJD0lKhm/DBwF3ElQxlTpMsjdJ9UhzkYvnZNxuVc/X8WNL85mY4wJQcYMyuWucw+ka7tKN1VERJq8RCXjDcDV7v6vRAbXlCgZB9Zt3sl/P/85782tPCFI+5bN+N1Z+8dclEJEpCnb6/WMo+wA1iUmJGnKurRrwaMXH8bvv3MArbIrPrnYtKOYnz01k589NZON23cnKUIRkdQRbzL+X+AH9RGIND1mxvmH9+aNa0ZzWN/K02m+8tkqTrpnMu/P1e93IpLe4r1NfTnBMomLgdeBStMtufujCYuuEdJt6thKy5x/fLCIP705j92lZZW2f++I3txw6hBaN89KQnQiIg0jUc+MK/8UrUijqZWMqzV3zRaue2YWX67eXGlb706tuHvcQeT17ZSEyERE6l+iknGfmvq4+9I4Y2tSlIxrtrukjHvfmc8D7y+oNJ2mGfx49L5cd+JAmmel9e91ItIEaQauBqJkXHszlm3gF89+xuLCbZW27bdPW+4edzBDu2s6TRFpOhI1mrr8YAea2VVmdks4GxdmNsDMKi/nI1KFQ3t3ZMLVI7nwqMo3XL5es4Uz7/+AB95fQGn05bOISBMTVzI2s+ZmNh6YCdwL3AyUF4veBaT9qk0Sn1bZWdx25v7864eHs0/URCDFpc5db8xl3ENTWRLj6llEpKmI98r4DuAEgvKmrnwzHSYEo6tPTlBckmZGDczlzWtHc/YhPSptm750A6f8dQr/+ngpeqwiIk1RvMn4fOBGd/83lcuaFgN9ExGUpKf2rZpxz3kH88AFh9KxVbMK23YUl3LTi7O56LFPWbNpZ5IiFBGpH/Em487AV9Ucq3lNBzCzY8zMY7w2RvR5vIo+bmZfRx2vhZn90cxWm9kOM5tqZqNjvG+GmV1vZkvMbKeZfWZm51QR42Vm9rWZ7TKzuWZ2RU2fSxLn1AO68eZ1ozl+vy6Vtk2eV8BJ90zipVkrdZUsIk1GvMl4McFCEbEcDsyN41hXh8cqf50Qse32qG1HEVyVA7wcdZx/AJcRPL8+HVgNvGlmB0f1ux24FbgPOAX4GBhvZqdGdjKzy4CHgOeBbwHjgQfM7Mo4PpvUUZe2LXjkojz+cM4BtI6aTnPzzhKueXoWVz01kw3bNJ2miDR+8dYZXw/8BriCIFltB4YDHYDngFvd/W81HOMY4D3gRHefGMd73wTcBuzv7nPCtoOAWcCl7v5Y2JYFzAHmuvsZYVsXYDlwp7vfEnHMd4Bcdz8wYt9VwOvuflFEv0eBM4Bu7l55KaIIKm1KvOVF2/nF+M+YtrjShG/ktm3OH845gOP265qEyERE4pOo0qa7gAkEaxpvCNs+ACYCb9SUiOvoQmB6eSIOnQEUA8+UN7h7CfA0cLKZld82PxnIBp6IOuYTwAFm1i/8/iggN0a/fxHcoh+ZgM8hcerVqRVPX3YkN5w6hOzMiv9kC7bs4tLH87n+hc/ZuqskSRGKiNRNXMnY3Uvd/bvAGODPwCMEJU7HufsFcb73k2ZWambrzezfZta7qo5mNgIYAPwzatMwYLG7b49qn0OQfAdE9NsFLIjRD2BoRD+A2TX0kwaWkWFcNro/r149kmExJgJ5atpyTvnr5JhXzyIiqW6vZuV39ynAlL18z00EiXwSsBk4hODW91QzO8TdYy3hcyHBFfBTUe2d+OYKPVJRxPbyrxu98j35WP2IcczofhWEC2hcDtC7d5W/U0gCDOralv/8ZAT3vTuf+99fWGFCkOVFOzjv4alcPqo/1504iBbNNJ2miDQOezUDV124+0x3/6W7v+Luk9z9LwQDpboSDOqqwMxaAOOAV929sIHDrRV3f9jd89w9Lzc3N9nhNHnZWRn8/KTBPHfFUfTPaV1hmzs8NHkRZ973IXNWbUpShCIi8Yl3Bq6y8NZyrFdJeMv5bTM7KZ7juvsMYB5wWIzNZxAMEIu+RQ3BFWzlhXK/uYItiujXwcysFv2IcczofpICDundkQlXj+Lio/tW2jZ37RbOuv9D7n9vASUxlmwUEUkl8V4Z304wKrkAeBz4A0GSLABWEAx0ygVeN7PT9yKeWEO7LwIKgddibJsD9DOzVlHtQ4HdfPOMeA5BDfS+MfoBfBnRD755dlxVP0kRLbMzufWMYTz5oyPo1r7ydJp/fHMuYx+aysKCrUmKUESkZvEm452EM225+w/d/TfufinQD1hCkJQPBd4ieA5cK2aWBwwGpkW1dyUYCf3vKkqKXgGaAWMj9skCzgPecvddYfMbBM+coweZfR+Y7e6Lw++nEiT+WP2KgA9r+5mkYY0YkMMb147mO4dWnk5z5rKNnHD3JC58dBqvfr6KncWlSYhQRKRq8Q7gugL4mbtXmI/Q3XeY2T3Afe5+h5k9QuzbypjZkwQJfQawkWAA1/XASoKR2ZEuADKrOpa7zzSzZ4C/mFmz8LhXEvxycEFEv3VmdjdwvZltCd/7POA4gtvg5f2Kw3rmB8xsJUHJ1nHApeHn1gwTKax9y2bcPe5gThrald/8ZzZFEROCuAezd02eV0D7ls046+DujM3rxf492icxYhGRQLzJOJfgSjSWbIJaXAiuLqOfz5abTTCb1s+AVsAa4AXglhgDtC4iuHKdUU1MlxAsYPE7gmfLnwHfirHPDcBW4BpgH4LZwsa5+6uRndz9QTNz4BfAr4BlwFXu/kA1MUgK+db+3RjepxPXv/A5E7+qPDh/045i/jl1Kf+cupQh3doxLq8nZx7cg06ts5MQrYhI/DNwTSFIeCe5++qI9u4Et6aL3H20mV0I3OTuAxMdcKrTDFypw915adYqHvlgEbNXbq62b7NM48ShXRk7vBejBuaQldnghQYikgaqmoEr3mR8KPAO0IJgbud1QBeCmau2E0z+McvMbgM8curJdKFknJrmrNrE+PwVvDRrJRu2VzujKV3bNec7h/Zk7PCe9M9t00ARikg6SEgyDg/UmeAW7hFAN4KFGT4G7nb39QmItVFTMk5tu0pKeeerdYzPX86keQWU1fDP/7C+HRk7vBenHtiNNs33ao4cEZE9EpaMpXpKxo3Hmk07eWHmCsbnr2Bx4bZq+7bKzuTUA7oxLq8Xh/XtSOWSdRGRmiU0GZtZJ4Jb052A9cDH7q4JMVAybozcnfylGxifv5xXP1/N9t3Vlz717W3/MvgAAB0uSURBVNyKsXm9+M6hPejWvmUDRSkiTUEib1P/juA2dfOI5l3An9z9pjpF2QQoGTdu23aV8NoXqxmfv4JpS6r//TLDYNTAXMbl9eKEoV1onqW5sEWkeokawHUtcDfwD4JlBtcQlAl9n6AW9zp3j64VTitKxk3H4sJtPDd9Oc9PX8mazTur7duhVTPOOrgH5w7vqdplEalSopLx18Dr7n5djG33AKe4+351irSRUzJuekrLnCnzCxifv4K3v1zL7hrmuh7arR1j83py1sE96KjaZRGJkKhkvBM43d0nxth2AsHKSi0q75k+lIybtg3bdvPSrJWMn76COauqr13OzszgxKFdOTevJ6MH5pKZoUFfIumuqmQcb63GemB/gmkiow0Lt4s0WR1bZ3PxiH5cPKLfntrlF2etZGOM2uXdpWVM+GI1E75YzT7tWvCdQ3swNq8X/aKWfRQRiffK+D6CKSp/Ajzl7iXhwgxjgQeBf7p7pTWJ04mujNPPrpJSJn65jvHTlzO5trXLeb047YButFbtskhaSdRt6rYESxmOAEoJVjLqRLCYwwfAqe6e1mvVKRmntzWbdvL8jBWMz1/OkvXbq+3bKjuT0w7oxrjDepHXR7XLIukgkaVNBpwGjCJIxEXAJIKBXWk/g4iSscA3tcvPfrqcCV/UXLvcL6c15w7vyTmH9mSf9mk97EKkSatzMjazbIJpL3/t7m8lOL4mQ8lYom3bVcKEL1YzPn85ny7ZUG3fDIPRg3IZO1y1yyJNUaJuU28AznH3dxMZXFOiZCzVWVSwleemr+D5GStYu3lXtX3La5fH5vVkWHfVLos0BYlKxs8Ci9z914kMrilRMpbaKC1zJs8vYHz+ct7+ci3FpdX/PxzWvR1jhwfrLqt2WaTxSlQyHkUw89Z44EWCFZsqHMDdF9Ut1MZNyVjiVRTWLj+bv4KvVteidnlYV8YO78ko1S6LNDqJSsaRUw/F3NHd0/ohl5Kx1MXslZt4bvoK/jNzJZt2VL/u8j7tWnDO8B6MHd6LvqpdFmkUEpWML6qpj7v/M87YmhQlY0mEncWlTPxqLePzVzB5fgE1/Tc9vG8nxub15FTVLoukNK1n3ECUjCXRVm/awQszVvJs/nKW1lC73Do7k9MODNZdHq7aZZGUk+j1jDOAoUBnIN/dq1+ZPY0oGUt9cXc+XbKBZ/OXM+Hz1eworr52uX9Oa87NC2qXu7ZT7bJIKkjkpB8/BW4BcgieGx/m7jPM7EXgXS2hqGQs9W/rrhJe+3w1z+YvJ39pzbXLYwYF6y4fP6Qr2VkZDRSliERL1DPjy4C/A48CbwHPAnlhMv4FcIa7j0lQzI2SkrE0tIXltcvTV7BuS/W1yx1bNeOsQ4JBX0O7t2ugCEWkXKKS8VfAy+7+32aWCRTzTTI+DfiHu++TsKgbISVjSZaS0jKmzC/k2fzlTPyq5trl/Xu0Y+zwXpx5cHc6tFLtskhDSNQSiv2AN6vYtg3oEG9gIpIYWZkZHLtfF47drwtF23bz4sxg0NfXa7bE7D975WZmr5zDHRO+4sRhXRmX14uRA3JUuyySBPEm40KgbxXbBgMr6xSNiCREp9bZXDqyH5eM6MucVZt5Nn85L81aFbN2eXdpGRM+X82Ez1fTrX0Lzjm0J+cO76naZZEGFO9t6geBbwHHAUsJblMPB5YTLKE4wd1/UQ9xNhq6TS2pamdxKW9/uZbx01cwpTa1y/06MS6vF6cesA+tslW7LJIIiXpmnAN8CPQCPgFGAx8B+wHrgKPdfVNCIm6klIylMVi1cQcvzFjBs/krWFZUc+3y6Qd2Z2xeT9Uui9RRIkub2gLXAicDXYD1wBvAPe5e/cS6aUDJWBqTsjJn2pIixuev4LUvale7PDavF985tIdql0X2gmbgaiBKxtJYbd1VwoTPV/Fs/gqm16J2+ZjBXRg7vKdql0XikKjb1IuAs939sxjb9icoe+pfp0gbOSVjaQoWFmxlfH6w7nJBDbXLnVpnc+bB3RmX14sh3VS7LFKdRK7adKS7T4uxLQ/4RKs2KRlL01FSWsbk+QU8++kK3vm6drXLZx3cgzGDchnQpY2eL4tESVSdMVSxdCKQB2zci+OJSIrKyszguP26ctx+XVm/dRcvzlrF+Bprlzfzuwlf0b19C0YPymX0oFxGDMihfctmDRy9SONR45WxmV0HXBd+2wMoAHZHdWsJdAKedvcLEh1kY6IrY2nq3J3ZKzczfvpyXpy5ks07S2rcJ8PgkN4dGT0wl9GDcjiwZwdNLiJpaa9vU5vZmcBZ4bcXAa8RJORIu4AvgUfcvfo6iSZOyVjSyc7iUt76ci3j85fzwYLCGmuXy3Vo1YyRA3IYPSiXMYNyNTJb0kainhk/Btzm7osTGVxTomQs6Wrlxh28MXsNk+cV8Mni9ewsLqv1voO7tmXM4FxGD8wlr29HWjRL66En0oQlvLTJzNoQrGe8yt0rz7GXppSMRYIr5k+XFDF5XgGT5hUwb+3WWu/bolkGR/XvvOd5c/+c1hoIJk1GIif9OB24DTgobCpfz/gRgvWM/13naBsxJWORylZv2sGUeYVMml/AB/MLY86RXZUeHVruuZ199IDOtGuhgWDSeCXqNvVZwPPAOwTrGd/FN0so3gCMdveTExRzo6RkLFK90jLnsxUbmTyvgMnzCpi1fCNltfwxlJlhHNq7A6MH5jJmcC77d29PhgaCSSOSqGQ8E5ju7j8ysyyCUdXlyfhM4AF375GwqBshJWOR+GzaXswHCwqD5Dy/gNWbdtZ6306tsxk5IIcxg3IZNSiHLm01EExSW6KS8U7g2+7+tpllEqzaVJ6MRwNvuXta/29QMhbZe+7O/HVb9zxr/mRxEbtLaj8QbEi3doweFCTnvD6dNE2npJxETfqxGcipYltfKpc8iYjUmpkxqGtbBnVty49G9WfH7lI+WbyeyfMKmTy/gAXrqh8I9tXqzXy1ejMPTVpEq+zMPQPBxgzK1frMktLivTJ+EjiAYOnELXyznvGXwBRglrtfXg9xNhq6MhapPys37tjzrPmDBYVsqcWEI+V6d2rF6EE5jB6Yy9EDcmjTXGs0S8NL1G3qvsA0gikxXwMuBJ4DDgTaE9yyXpWAeBstJWORhlFSWsas5cFAsEnzC/l8xcZaTzqSlWEM79Nxz1Xz0G7tNBBMGkQiS5t6Ar+l8nrGN7v78gTE2qgpGYskx4Ztu5lSPhBsXgHralhtKlJOm2xGhVN1jhqYS06b5vUYqaQzrWfcQJSMRZLP3Zm7dguT5gYjtD9dvIHdpbUfCLZ/j3bhPNq5HNq7owaCScLUZW7qm+N4H3f32+MNrilRMhZJPdt3l/DJoiImhVfNiwq31Xrf1tmZHF0+j/bAXHp3blWPkUpTV5dkHOvXSQdiPWBxrWesZCyS6pYXbWfy/AImzS3go4Xr2bqr9gPB+nZuxZhwqs4j+3emtQaCSRzqkoyjk2sWsAM4ApgR3d/dS+sQZ6OnZCzSuBSXljFz2cY9tc1frNxU632bZRp5fTrtWeRiSLe2mkdbqpXIAVwVJvtIUHxNhpKxSOO2fusuPlhQGD5vLqRwa+0HguW2bc6ogeGMYANz6dQ6ux4jlcZIybiBKBmLNB1lZc5XazYHk47MKyB/aRHFpbX7mWkGB/Rov2ce7YN7daBZpgaCpTsl4waiZCzSdG3bVcLUheuZPD8YCLZk/fZa79u2eRZHDwiXhhyYS69OGgiWjpSMG4iSsUj6WLp+W/isuZCPFhayfXfth8z0z20dXDWHA8FaZqf12Ne0UZcBXP2jmjKBucCZwJzo/u6+qA5xNnpKxiLpaXdJGdOXbthz1Txn1eZa75udlcHhfTsF03UOymVwVw0Ea6rqWtoU3clitAGg0iYlYxGBgi27mBIm5inzC1m/bXet9+3arvmeSUdGDsihowaCNRl1ScYXxfNG7v7PGo53DPBejE2b3L1DVN8jgVuBI4FmwCLgDnd/OqJPC+B24PtAB2AW8N/uPjnqWBnAfwM/BvYhuLq/zd2fjxHjZcAvgH7AEuAed3+wus9VTslYRKKVlTlzVm0OapvnFTBj6QZKymo/EOygnh3CebRzOKhnB7I0EKzRSpnpMCOS8dXApxGbStw9P6LfacB/gH8DzwK7gaHAZnd/PKLfk8BpwK8IkvVPgVOAo9x9VkS/O4BfAjcA04HvApcBp7v7axH9LgMeAn4PTASOB34D/NTd/17T51MyFpGabNlZzEcL1wfzaM8vYHnRjlrv265FFiMH5uy5cu7eoWU9RiqJlorJ+ER3n1hFn7bAQuDf7n5tNcc6iOBK+FJ3fyxsyyJ4lj3X3c8I27oAy4E73f2WiP3fAXLd/cCIfVcBr7v7RRH9HgXOALq5e3F1n0/JWETi4e4sWb+dSXPXMXl+IVMXrmdHce0Hgg3s0iYYoT0olyP6daJFs7R+UpjyqkrGqTqP21ggF/hzDf3OIBjZ/Ux5g7uXmNnTwK/NrLm77yJYYSobeCJq/yeAR82sn7svBo4K3ze637+AS4CRxL7FLiKyV8yMfjmt6ZfTj4tH9GNXSSnTl2xgUjgj2NdrtlS7//x1W5m/biv/+GAxzbMyOLxfJ8aES0MO6NJGA8EaiWQm4yfNLAfYCLwJ/Nrdl4XbRgJFwAFm9howBFgNPAL8LmLKzWHAYnePLvabQ5B8B4R/HgbsAhbE6AfB7e/FYT+A2dX0UzIWkXrTPCtYmOLoATlcf+oQ1m7eGd7OLuSD+QVs2F71zbldJWVMmV/IlPmF/G7CV3Rr32LPpCMj9s2hfatmDfhJJB7JSMabCK54JwGbgUMInslONbND3H0d0B1oRfC8+HaCZ7wnADcRDNK6LjxWJ2BDjPcoithe/nWjV74nH6sfMY4Z3a8CM7scuBygd+/esbqIiOyVru1aMDavF2PzelFa5sxeuWnP6lMzl2+ktJqBYKs37eSZ/OU8k7+cDIODe3VgzKAujB6Uw4E9O5CZoavmVNHgydjdZwIzI5ommdlkYBrBoK4bgQygBXCDu98d9nvfzDoDPzWzW9299rO51zN3fxh4GIJnxkkOR0SaqMwM46BeHTioVweuPn4gm3YUM3VhYZicC1m5seqBYGUOM5ZtZMayjdwzcR4dWjVjxIBgHu3RA3PZp32LBvwkEi0lnhm7+wwzmwccFjatD7++HdX1LeAKgtvJHxFcwfaJccjyK9jyK9oNQAczs6ir41j9ADoS3Bavqp+ISNK1b9mMb+3fjW/t3w13Z2HBtj2rT32yeD07i2OtgBvYuL2YCZ+vZsLnwY+6wV3bMnpQDmMGdSGvb0cNBGtgKZGMI5Qnykoze0Up/xc2BzjbzFpFPTceSlAKtSCiX3NgXyo+Nx4afv0y6n2HUTEZR/cTEUkpZsaALm0Y0KUNl47sx87iUj5dUrQnOc9bu7Xa/eeu3cLctVv43ymLadEsgyP7d95TPrVvbmsNBKtnDV7aFDMIszzgE4IJPW42s/2BL4BfufufIvo9BFxIUI601cwOIVhT+eLyyUbC8qQvgAXu/u2wrQuwIjz+byOONxHo6u4HhN83IyhtetXdL4no9whwNkFpU7XT6Ki0SURS0epNO5gyr5BJ8wv4YH4hm3ZUW6VZQY8OLfdMOnL0gBzatdBAsL2VMqVN4SQdiwmS6EaCAVzXAyuBewHcfbaZPQ7cFs6cNYNgANePgNvdfWvYb6aZPQP8JUyki4ErCWbOuqD8Pd19nZndDVxvZlvC450HHEdQHlXer9jMbgIeMLOVBJN+HAdcCvyspkQsIpKqurVvybjDejHusGAg2GcrNgajtOcVMGv5RqqbEGzlxh08NW0ZT01bRmaGcWjvDnuumg/o0Z4MDQSrs2RM+nE9cD7Bs95WwBrgdeAWd18d0S8buBm4COhKMC3l/e7+16jjtQTuAL5HMNL6M4LpMN+P6pdJkPQvo+J0mM/FiPHHBNNh9gGWEUyH+UBtPp+ujEWksdm4fTcfLli/55b2ms07a71vp9bZjByQEy4NmUOXdhoIVp2UmYGrqVMyFpHGzN2Zv25rxECwInaXVD0QLNqQbu2CgWADcxnetyPNszQQLJKScQNRMhaRpmTH7lI+WbyeyfMKmTy/gAXrqh8IFqlVdiZH9e+8Z7rOvp1bpf1AMCXjBqJkLCJN2cqNO/Y8a/5gQSFbdpbUet9enVruqWs+ekAObZqnWkFP/VMybiBKxiKSLkpKy5i1PBgINml+IZ+v2EhtU0pWhnFon4575tEe2q1dWgwEUzJuIErGIpKuNmzbzZQFhXuunNdt2VXrfXPaBAPBxgzOZdTAXHLaNK/HSJNHybiBKBmLiAQDwb5es2XPms2fLt7A7tLaDwQb1r1dWNucy6G9O5KdlVGP0TYcJeMGomQsIlLZ9t0lfLKoaM8iF4sKt9V639bZmRy1bw5jwuk6e3duVY+R1i8l4waiZCwiUrPlRduZPL+ASXML+Gjherbuqv1AsL6dW4V1zbkctW9nWjeigWBKxg1EyVhEJD7FpWXMXLZxT23zFytrvyhfs0wjr0+nsHwqh6Hd2qV0+ZSScQNRMhYRqZv1W3fxwYJCJs0tYPL8Qgq31n4gWG7b5owaGCwNOXJADp1TbCCYknEDUTIWEUmcsjLnqzWbg0lH5hWQv7SI4tLa5S0zOKBH+z3zaB/SuwPNMpM7EEzJuIEoGYuI1J9tu0qYunB98Lx5XgFL12+veadQ2+ZZHD2g857nzb06NfxAMCXjBqJkLCLScJau3xY+ay7ko4WFbN9dWut9++e03lM+dUT/TrTKrv+BYErGDUTJWEQkOXaXlDF96QYmzw/Kp+as2lzrfbMzMzi8XydGDwpWoBrctW29DARTMm4gSsYiIqmhYMsupoSJefL8Qoq21X5J+q7tmu951jxyQA4dW2cnJCYl4waiZCwiknrKypw5qzbvedY8Y+kGSspqPxDswJ4dGDMwh0tG9KtTYlYybiBKxiIiqW/LzmI+Wrh+z3Sdy4t21LhPVoYx8+YTadui2V6/b1XJuPFMWyIiIpIgbVs04+Rh+3DysH1wd5as386kueuYPL+QqQvXs6O48kCwQ/t0rFMiro6SsYiIpDUzo19Oa/rl9OPiEf3YVVJK/pINe2YE+3rNFgDGDMqttxiUjEVERCI0z8pkxIAcRgzI4fpTh7B2804mzyvgsL6d6u09lYxFRESq0bVdC8bm9arX92gaC0SKiIg0YkrGIiIiSaZkLCIikmRKxiIiIkmmZCwiIpJkmoErwcysAFhay+45QGE9htPY6fzUTOeoZjpHNdM5qlmizlEfd69UsKxknERmlh9rWjQJ6PzUTOeoZjpHNdM5qll9nyPdphYREUkyJWMREZEkUzJOroeTHUCK0/mpmc5RzXSOaqZzVLN6PUd6ZiwiIpJkujIWERFJMiVjERGRJFMybkBm1svMnjOzTWa22cxeMLPeyY4rGczsXDN73syWmtkOM5trZr83s7ZR/Tqa2SNmVmhm28xsopkdkKy4k83M3jAzN7PfRbWn9Xkys1PNbLKZbQ3/b+Wb2XER29P9/Iwws7fMbJ2ZbTGzGWZ2aVSfFmb2RzNbHf6fnGpmo5MVc30xs55m9rfw820P/z/1jdGvVufDzDLM7HozW2JmO83sMzM7J964lIwbiJm1At4F9gMuAn4ADATeM7PWyYwtSX4JlAK/Ab4F/B24EnjbzDIAzMyAV8LtPwPOAZoRnLOeyQg6mczsfOCgGO1pfZ7M7MfAS8B04GxgLDAeaBVuT/fzcyAwkeAzXwZ8B/gU+IeZXRnR9R/h9puB04HVwJtmdnDDRlzvBgDjgA3AlGr61fZ83A7cCtwHnAJ8DIw3s1Pjisrd9WqAF3ANQfIZENHWDygBfp7s+JJwPnJjtF0IOHBc+P2Z4ffHRvRpDxQB9yb7MzTw+eoIrAHOD8/J7yK2pe15AvoCO4Brq+mTtucn/Kz/A+wG2kS1TwWmhn8+KDxHl0RszwLmAi8n+zMk+HxkRPz5R+Hn7hvVp1bnA+gC7AJ+G7X/O8Dn8cSlK+OGcwbwsbsvKG9w98XAhwQ/LNKKuxfEaP40/Noj/HoGsMrd34vYbxPBVU66nbM/ALPd/akY29L5PF0KlAEPVtMnnc8PQDZQTPBLS6RNfHN39IywzzPlG929BHgaONnMmjdAnA3C3ctq0a225+NkgvP7RNT+TwAHmFm/2salZNxwhgGzY7TPAYY2cCypakz49avwa3XnrLeZtWmQqJLMzEYS3DX4aRVd0vk8jQS+Br5rZgvNrMTMFphZ5LlK5/MD8Hj49V4z625mHczsMuB44J5w2zBgsbtvj9p3DkGyGdAgkaaO2p6PYQRXxgti9IM4frYrGTecTgTPKKIVEdyCTGtm1gO4DZjo7vlhc3XnDNLgvJlZNvAQ8Cd3n1tFt3Q+T90Jxl78EbgTOAl4G7jPzK4J+6Tz+cHdZwPHENwFWElwLu4HrnD3p8NuNZ2jTvUcZqqp7fnoBGz08N50Nf1qlBVXeCL1ILwyeYng+fklSQ4n1fwX0BK4I9mBpKgMoC1wsbu/ELa9G46Ovd7M7k1WYKnCzAYCzxNcrV1BcLv6TOBBM9vp7k8mMz4JKBk3nA3E/g28qt/A0oKZtSR4dtcfGOPuKyI2V3fOyrc3WWHZ2w0Eg0yaRz23a25mHYAtpPd5Wk9wZfx2VPtbBKOnu5He5weCAVzFwOnuXhy2vWNmnYG/mtlTBOegT4x9y89RUYxtTVltz8cGoIOZWdTVcdznTbepG84cgucL0YYCXzZwLCnBzJoBzwF5wKnu/kVUl+rO2TJ331rPISZbf6AFwWCQDREvCErDNgAHkN7naU4N28tI7/MDwb+RzyIScblpQGeCEcFzgH5hCWakoQQjsaOfiTZ1tT0fc4DmwL4x+kEcP9uVjBvOy8CRZta/vCG8lTYi3JZWwlriJ4HjgLPc/eMY3V4GepjZmIj92gHfJj3O2Szg2BgvCBL0sQQ/FNL5PP0n/HpyVPu3gBXuvob0Pj8QlMQdHI4/iHQEsJPg6u0VgjrkseUbzSwLOA94y913NVCsqaK25+MNgrsOF0Tt/32C6ofFtX7HZNd8pcsLaE3wg/MLguc1ZwCfAYuIqv9LhxfBJB8O/A44MurVM+yTAXwELAe+S/AD932CHx69kv0ZknjuouuM0/Y8AUYwmc56guehJwH/G56ji9P9/ISf/9zwfLwZ/uw5iWCCCgfujuj3NMHdlh8RjLR+jiBZH5rsz1BP5+TciJ9DV4bfj4n3fBAMHNwJ/JxgoNzfCe7InB5XTMk+Ken0AnoTDKTYTPCs70Wiis3T5QUsCf8TxHrdGtGvE/Bo+INzO0Ex/UHJjj/J565CMk738wS0IxgdvJbgFuLnwPd0fip8/lPCX0AKwp89s4CfAJkRfVoCdxNcSe8EPgGOSXbs9XQ+qvrZ83685wPIBG4ElhKUOX0OnBtvTFpCUUREJMn0zFhERCTJlIxFRESSTMlYREQkyZSMRUREkkzJWEREJMmUjEVERJJMyVhEEsrMLjYzN7MBUe2HmVmRmc00s5xkxSeSipSMRaTemdnRwERgPnCcuxcmOSSRlKJkLCL1KpwT+k2CqWBPdPemvkqSSNyUjEWk3pjZicDrwKfAye6+OckhiaQkJWMRqS+nEax+Mxk4zd23JTkekZSlZCwi9eUvwArgTHffkexgRFKZkrGI1JcJBIuuX5/sQERSXVayAxCRJus6guXnbjGzHe7+h2QHJJKqlIxFpL44cDnQArjTzHa5+1+SHJNISlIyFpF64+5lZnYRkA3cY2Y73f3BZMclkmqUjEWkXrl7qZl9jyAhPxBeIT+W7LhEUokGcIlIvXP3EmAc8AbwSJicRSRk7p7sGERERNKaroxFRESSTMlYREQkyZSMRUREkkzJWEREJMmUjEVERJJMyVhERCTJlIxFRESSTMlYREQkyf4f2A9Y5SbldckAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 504x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U_H4s9hDJosX"
      },
      "source": [
        "In the above plot we show that heterogeneity goes down as we increase the number of clusters. Does this mean we should always favor a higher K? **Not at all!** As we will see in the following section, setting K too high may end up separating data points that are actually pretty alike. At the extreme, we can set individual data points to be their own clusters (K=N) and achieve zero heterogeneity, but separating each data point into its own cluster is hardly a desirable outcome. In the following section, we will learn how to detect a K set \"too large\"."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JXtnCPXLJosY"
      },
      "source": [
        "## Visualize clusters of documents"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ED0-clbyJosZ"
      },
      "source": [
        "Let's start visualizing some clustering results to see if we think the clustering makes sense.  We can use such visualizations to help us assess whether we have set K too large or too small for a given application.  Following the theme of this course, we will judge whether the clustering makes sense in the context of document analysis.\n",
        "\n",
        "What are we looking for in a good clustering of documents?\n",
        "* Documents in the same cluster should be similar.\n",
        "* Documents from different clusters should be less similar.\n",
        "\n",
        "So a bad clustering exhibits either of two symptoms:\n",
        "* Documents in a cluster have mixed content.\n",
        "* Documents with similar content are divided up and put into different clusters.\n",
        "\n",
        "To help visualize the clustering, we do the following:\n",
        "* Fetch nearest neighbors of each centroid from the set of documents assigned to that cluster. We will consider these documents as being representative of the cluster.\n",
        "* Print titles and first sentences of those nearest neighbors.\n",
        "* Print top 5 words that have highest tf-idf weights in each centroid."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FJ6ARWIsJosa"
      },
      "source": [
        "def visualize_document_clusters(wiki, tf_idf, centroids, cluster_assignment, k, map_word_to_index, display_content=True):\n",
        "    '''wiki: original dataframe\n",
        "       tf_idf: data matrix, sparse matrix format\n",
        "       map_index_to_word: SFrame specifying the mapping betweeen words and column indices\n",
        "       display_content: if True, display 8 nearest neighbors of each centroid'''\n",
        "    map_index_to_word =  {v:k for k,v in map_word_to_index.items()}\n",
        "    print('==========================================================')\n",
        "    # Visualize each cluster c\n",
        "    for c in range(k):\n",
        "        # Cluster heading\n",
        "        print('Cluster {0:d}    '.format(c)),\n",
        "        # Print top 5 words with largest TF-IDF weights in the cluster\n",
        "        idx = centroids[c].argsort()[::-1]\n",
        "        for i in range(5): # Print each word along with the TF-IDF weight\n",
        "            print('{0:s}:{1:.3f}'.format(map_index_to_word[idx[i]], centroids[c][idx[i]])),\n",
        "        print('')\n",
        "        \n",
        "        if display_content:\n",
        "            # Compute distances from the centroid to all data points in the cluster,\n",
        "            # and compute nearest neighbors of the centroids within the cluster.\n",
        "            distances = pairwise_distances(tf_idf, centroids[c].reshape(1, -1), metric='euclidean').flatten()\n",
        "            distances[cluster_assignment!=c] = float('inf') # remove non-members from consideration\n",
        "            \n",
        "            nearest_neighbors = distances.argsort()\n",
        "            \n",
        "            # For 8 nearest neighbors, print the title as well as first 180 characters of text.\n",
        "            # Wrap the text at 80-character mark.\n",
        "            for i in range(8):\n",
        "                text = ' '.join(wiki[nearest_neighbors[i]]['text'].split(None, 25)[0:25])\n",
        "                print('\\n* {0:50s} {1:.5f}\\n  {2:s}\\n  {3:s}'.format(wiki[nearest_neighbors[i]]['name'],\n",
        "                    distances[nearest_neighbors[i]], text[:90], text[90:180] if len(text) > 90 else ''))\n",
        "        print('==========================================================')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J0GcSkbxJosb"
      },
      "source": [
        "Let us first look at the 2 cluster case (K=2)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qtjB_e3WJosc",
        "outputId": "12ad518f-9b57-4113-9ac4-19b99c5ceff8"
      },
      "source": [
        "'''Notice the extra pairs of parentheses for centroids and cluster_assignment.\n",
        "   The centroid and cluster_assignment are still inside the npz file,\n",
        "   and we need to explicitly indicate when to load them into memory.'''\n",
        "visualize_document_clusters(wiki, tf_idf, centroids[2], cluster_assignment[2], 2, map_index_to_word)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "==========================================================\n",
            "Cluster 0    \n",
            "she:0.021\n",
            "university:0.015\n",
            "her:0.013\n",
            "he:0.012\n",
            "served:0.010\n",
            "\n",
            "\n",
            "* Kayee Griffin                                      0.97358\n",
            "  kayee frances griffin born 6 february 1950 is an australian politician and former australi\n",
            "  an labor party member of the new south wales legislative council serving\n",
            "\n",
            "* %C3%81ine Hyland                                   0.97370\n",
            "  ine hyland ne donlon is emeritus professor of education and former vicepresident of univer\n",
            "  sity college cork ireland she was born in 1942 in athboy co\n",
            "\n",
            "* Christine Robertson                                0.97373\n",
            "  christine mary robertson born 5 october 1948 is an australian politician and former austra\n",
            "  lian labor party member of the new south wales legislative council serving\n",
            "\n",
            "* Anita Kunz                                         0.97471\n",
            "  anita e kunz oc born 1956 is a canadianborn artist and illustratorkunz has lived in london\n",
            "   new york and toronto contributing to magazines and working\n",
            "\n",
            "* Barry Sullivan (lawyer)                            0.97488\n",
            "  barry sullivan is a chicago lawyer and as of july 1 2009 the cooney conway chair in advoca\n",
            "  cy at loyola university chicago school of law\n",
            "\n",
            "* Margaret Catley-Carlson                            0.97534\n",
            "  margaret catleycarlson oc born 6 october 1942 is a canadian civil servant she was chair an\n",
            "  d is now a patron of the global water partnership\n",
            "\n",
            "* Vanessa Gilmore                                    0.97579\n",
            "  vanessa diane gilmore born october 1956 is a judge on the united states district court for\n",
            "   the southern district of texas she was appointed to\n",
            "\n",
            "* James A. Joseph                                    0.97624\n",
            "  james a joseph born 1935 is an american former diplomatjoseph is professor of the practice\n",
            "   of public policy studies at duke university and founder of\n",
            "==========================================================\n",
            "Cluster 1    \n",
            "she:0.023\n",
            "music:0.017\n",
            "her:0.017\n",
            "league:0.016\n",
            "season:0.016\n",
            "\n",
            "\n",
            "* Patricia Scott                                     0.97143\n",
            "  patricia scott pat born july 14 1929 is a former pitcher who played in the allamerican gir\n",
            "  ls professional baseball league for parts of four seasons\n",
            "\n",
            "* Madonna (entertainer)                              0.97181\n",
            "  madonna louise ciccone tkoni born august 16 1958 is an american singer songwriter actress \n",
            "  and businesswoman she achieved popularity by pushing the boundaries of lyrical\n",
            "\n",
            "* Janet Jackson                                      0.97257\n",
            "  janet damita jo jackson born may 16 1966 is an american singer songwriter and actress know\n",
            "  n for a series of sonically innovative socially conscious and\n",
            "\n",
            "* Natashia Williams                                  0.97343\n",
            "  natashia williamsblach born august 2 1978 is an american actress and former wonderbra camp\n",
            "  aign model who is perhaps best known for her role as shane\n",
            "\n",
            "* Todd Williams                                      0.97384\n",
            "  todd michael williams born february 13 1971 in syracuse new york is a former major league \n",
            "  baseball relief pitcher he attended east syracuseminoa high school\n",
            "\n",
            "* Marilyn Jenkins                                    0.97430\n",
            "  marilyn a jenkins jenks born september 18 1934 is a former catcher who played in the allam\n",
            "  erican girls professional baseball league listed at 5 ft\n",
            "\n",
            "* Kayla Bashore Smedley                              0.97496\n",
            "  kayla bashore born february 20 1983 in daegu south korea is an american field hockey defen\n",
            "  der and midfielder now living in san diego california she\n",
            "\n",
            "* Cher                                               0.97510\n",
            "  cher r born cherilyn sarkisian may 20 1946 is an american singer actress and television ho\n",
            "  st described as embodying female autonomy in a maledominated industry\n",
            "==========================================================\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jeydihm8Josd"
      },
      "source": [
        "Both clusters have mixed content, although cluster 1 is much purer than cluster 0:\n",
        "* Cluster 0: academia, law\n",
        "* Cluster 1: female figures, baseball players\n",
        "\n",
        "Roughly speaking, the entire dataset was divided into athletes and non-athletes. It would be better if we sub-divided non-atheletes into more categories. So let us use more clusters. How about `K=10`?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "99wFQtwVJose",
        "outputId": "1c56c141-3c0f-44c4-bce3-daf50bc65645"
      },
      "source": [
        "k = 10\n",
        "visualize_document_clusters(wiki, tf_idf, centroids[k], cluster_assignment[k], k, map_index_to_word)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "==========================================================\n",
            "Cluster 0    \n",
            "he:0.012\n",
            "art:0.011\n",
            "his:0.009\n",
            "book:0.008\n",
            "that:0.008\n",
            "\n",
            "\n",
            "* Wilson McLean                                      0.97661\n",
            "  wilson mclean born 1937 is a scottish illustrator and artist he has illustrated primarily \n",
            "  in the field of advertising but has also provided cover art\n",
            "\n",
            "* Tang Xiyang                                        0.97988\n",
            "  tang xiyang born january 30 1930 in miluo hunan province is a chinese environmentalist he \n",
            "  was awarded the 2007 ramon magsaysay award for peace and\n",
            "\n",
            "* David Salle                                        0.98168\n",
            "  david salle born 1952 is an american painter printmaker and stage designer who helped defi\n",
            "  ne postmodern sensibility salle was born in norman oklahoma he earned\n",
            "\n",
            "* Alberto Blanco (poet)                              0.98172\n",
            "  alberto blanco is considered one of mexicos most important poets born in mexico city on fe\n",
            "  bruary 18 1951 he spent his childhood and adolescence in\n",
            "\n",
            "* John Donald (jewellery designer)                   0.98290\n",
            "  john donald is a british jeweller designer whose work is strongly identified in the 1960s \n",
            "  and 1970s in london princess margaret and the queen mother\n",
            "\n",
            "* David Elliott (curator)                            0.98298\n",
            "  david stuart elliott born 29 april 1949 is a britishborn art gallery and museum curator an\n",
            "  d writer about modern and contemporary arthe was educated at\n",
            "\n",
            "* Chris Hunt                                         0.98307\n",
            "  chris hunt is a british journalist magazine editor and author he has worked in journalism \n",
            "  for over twenty years most often writing about football or\n",
            "\n",
            "* Kcho                                               0.98342\n",
            "  kchosometimes spelled kcho born alexis leiva machado on the isla de pinos 1970 is a contem\n",
            "  porary cuban artist kcho has had art showings around the\n",
            "==========================================================\n",
            "Cluster 1    \n",
            "film:0.088\n",
            "theatre:0.037\n",
            "films:0.032\n",
            "television:0.028\n",
            "actor:0.027\n",
            "\n",
            "\n",
            "* Shona Auerbach                                     0.93531\n",
            "  shona auerbach is a british film director and cinematographerauerbach began her career as \n",
            "  a stills photographer she studied film at manchester university and cinematography at\n",
            "\n",
            "* Singeetam Srinivasa Rao                            0.93748\n",
            "  singeetam srinivasa rao born 21 september 1931 is an indian film director producer screenw\n",
            "  riter composer singer lyricist and actor known for his works in telugu\n",
            "\n",
            "* Justin Edgar                                       0.93801\n",
            "  justin edgar is a british film directorborn in handsworth birmingham on 18 august 1971 edg\n",
            "  ar graduated from portsmouth university in 1996 with a first class\n",
            "\n",
            "* Laura Neri                                         0.94151\n",
            "  laura neri greek is a director of greek and italian origins born in brussels belgium she g\n",
            "  raduated from the usc school of cinematic arts in\n",
            "\n",
            "* Bill Bennett (director)                            0.94260\n",
            "  bill bennett born 1953 is an australian film director producer and screenwriterhe dropped \n",
            "  out of medicine at queensland university in 1972 and joined the australian\n",
            "\n",
            "* Robert Braiden                                     0.94344\n",
            "  robert braiden is an australian film director and writer born in sydney he grew up in moor\n",
            "  ebank liverpool new south wales and now currently lives\n",
            "\n",
            "* Nitzan Gilady                                      0.94369\n",
            "  nitzan gilady also known as nitzan giladi hebrew is an israeli film director who has writt\n",
            "  en produced and directed the documentary films in satmar custody\n",
            "\n",
            "* Robb Moss                                          0.94484\n",
            "  robb moss is an independent documentary filmmaker and professor at harvard university nota\n",
            "  ble work includes such films as the same river twice secrecy film and\n",
            "==========================================================\n",
            "Cluster 2    \n",
            "league:0.061\n",
            "baseball:0.048\n",
            "season:0.046\n",
            "coach:0.042\n",
            "games:0.034\n",
            "\n",
            "\n",
            "* Todd Williams                                      0.92759\n",
            "  todd michael williams born february 13 1971 in syracuse new york is a former major league \n",
            "  baseball relief pitcher he attended east syracuseminoa high school\n",
            "\n",
            "* Justin Knoedler                                    0.93295\n",
            "  justin joseph knoedler born july 17 1980 in springfield illinois is a former major league \n",
            "  baseball catcherknoedler was originally drafted by the st louis cardinals\n",
            "\n",
            "* Kevin Nicholson (baseball)                         0.93579\n",
            "  kevin ronald nicholson born march 29 1976 is a canadian baseball shortstop he played part \n",
            "  of the 2000 season for the san diego padres of\n",
            "\n",
            "* Dave Ford                                          0.93642\n",
            "  david alan ford born december 29 1956 is a former major league baseball pitcher for the ba\n",
            "  ltimore orioles born in cleveland ohio ford attended lincolnwest\n",
            "\n",
            "* Steve Springer                                     0.93649\n",
            "  steven michael springer born february 11 1961 is an american former professional baseball \n",
            "  player who appeared in major league baseball as a third baseman and\n",
            "\n",
            "* Chris Young (pitcher)                              0.93772\n",
            "  christopher ryan chris young born may 25 1979 is an american professional baseball rightha\n",
            "  nded pitcher who is a free agent he made his major league\n",
            "\n",
            "* Eric Fox                                           0.93875\n",
            "  eric hollis fox born august 15 1963 in lemoore california is an american professional base\n",
            "  ball coach the 5 ft 10 in 178 m 180 lb\n",
            "\n",
            "* Ted Silva                                          0.93948\n",
            "  theodore a silva born august 4 1974 in inglewood california has held numerous roles in ama\n",
            "  teur and professional baseball he has played in the minor\n",
            "==========================================================\n",
            "Cluster 3    \n",
            "party:0.046\n",
            "election:0.042\n",
            "minister:0.039\n",
            "elected:0.028\n",
            "member:0.020\n",
            "\n",
            "\n",
            "* Stephen Harper                                     0.95128\n",
            "  stephen joseph harper pc mp born april 30 1959 is a canadian politician who is the 22nd an\n",
            "  d current prime minister of canada and the\n",
            "\n",
            "* Lucienne Robillard                                 0.95307\n",
            "  lucienne robillard pc born june 16 1945 is a canadian politician and a member of the liber\n",
            "  al party of canada she sat in the house\n",
            "\n",
            "* Marcelle Mersereau                                 0.95379\n",
            "  marcelle mersereau born february 14 1942 in pointeverte new brunswick is a canadian politi\n",
            "  cian a civil servant for most of her career she also served\n",
            "\n",
            "* Maureen Lyster                                     0.95450\n",
            "  maureen anne lyster born 10 september 1943 is an australian politician she was an australi\n",
            "  an labor party member of the victorian legislative assembly from 1985\n",
            "\n",
            "* Bruce Flegg                                        0.95554\n",
            "  dr bruce stephen flegg born 10 march 1954 in sydney is an australian former politician he \n",
            "  was a member of the queensland legislative assembly from\n",
            "\n",
            "* Doug Lewis                                         0.95583\n",
            "  douglas grinslade doug lewis pc qc born april 17 1938 is a former canadian politician a ch\n",
            "  artered accountant and lawyer by training lewis entered the\n",
            "\n",
            "* Paul Martin                                        0.95583\n",
            "  paul edgar philippe martin pc cc born august 28 1938 also known as paul martin jr is a can\n",
            "  adian politician who was the 21st prime\n",
            "\n",
            "* Gordon Gibson                                      0.95612\n",
            "  gordon gibson obc born 1937 is a political columnist author and former politician in briti\n",
            "  sh columbia bc canada he is the son of the late\n",
            "==========================================================\n",
            "Cluster 4    \n",
            "music:0.095\n",
            "orchestra:0.087\n",
            "symphony:0.057\n",
            "opera:0.050\n",
            "conductor:0.041\n",
            "\n",
            "\n",
            "* Heiichiro Ohyama                                   0.89162\n",
            "  heiichiro ohyama yama heiichir born 1947 in kyoto japan is a japanese conductor and violis\n",
            "  the has a longestablished reputation as a remarkable conductor and one\n",
            "\n",
            "* Brenton Broadstock                                 0.90284\n",
            "  brenton broadstock ao born 1952 is an australian composerbroadstock was born in melbourne \n",
            "  he studied history politics and music at monash university and later composition\n",
            "\n",
            "* Toshiyuki Shimada                                  0.90301\n",
            "  toshiyuki shimada is a japanese american orchestral conductor he is music director of both\n",
            "   the eastern connecticut symphony orchestra in new london ct and the\n",
            "\n",
            "* David Porcelijn                                    0.90357\n",
            "  david porcelijn born 7 january 1947 in achtkarspelen is a dutch composer and conductordavi\n",
            "  d porcelijn studied flute composition and conducting at the royal conservatoire of\n",
            "\n",
            "* Hugh Wolff                                         0.90731\n",
            "  hugh wolff born 21 october 1953 in paris is an american conductorhe was born in paris whil\n",
            "  e his father was serving in the u s\n",
            "\n",
            "* Daniel Meyer (conductor)                           0.90849\n",
            "  daniel meyer was born in cleveland ohio and has been conductor and musical director of sev\n",
            "  eral prominent american orchestrashe is a graduate of denison university\n",
            "\n",
            "* Paul Hostetter                                     0.90967\n",
            "  paul hostetter is the ethel foley distinguished chair in orchestral activities for the sch\n",
            "  wob school of music at columbus state university the conductor and artistic\n",
            "\n",
            "* Peter Ruzicka                                      0.91097\n",
            "  peter ruzicka born july 3 1948 is a german composer and conductor of classical musicpeter \n",
            "  ruzicka was born in dsseldorf on july 3 1948 he\n",
            "==========================================================\n",
            "Cluster 5    \n",
            "she:0.140\n",
            "her:0.088\n",
            "miss:0.012\n",
            "actress:0.011\n",
            "womens:0.011\n",
            "\n",
            "\n",
            "* Lauren Royal                                       0.93427\n",
            "  lauren royal born march 3 circa 1965 is a book writer from california royal has written bo\n",
            "  th historic and novelistic booksa selfproclaimed angels baseball fan\n",
            "\n",
            "* Janine Shepherd                                    0.93681\n",
            "  janine lee shepherd am born 1962 is an australian pilot and former crosscountry skier shep\n",
            "  herds career as an athlete ended when she suffered major injuries\n",
            "\n",
            "* Barbara Hershey                                    0.93708\n",
            "  barbara hershey born barbara lynn herzstein february 5 1948 once known as barbara seagull \n",
            "  is an american actress in a career spanning nearly 50 years\n",
            "\n",
            "* Janet Jackson                                      0.93735\n",
            "  janet damita jo jackson born may 16 1966 is an american singer songwriter and actress know\n",
            "  n for a series of sonically innovative socially conscious and\n",
            "\n",
            "* Ellina Graypel                                     0.93837\n",
            "  ellina graypel born july 19 1972 is an awardwinning russian singersongwriter she was born \n",
            "  near the volga river in the heart of russia she spent\n",
            "\n",
            "* Alexandra Potter                                   0.93867\n",
            "  alexandra potter born 1970 is a british author of romantic comediesborn in bradford yorksh\n",
            "  ire england and educated at liverpool university gaining an honors degree in\n",
            "\n",
            "* Dorothy E. Smith                                   0.93904\n",
            "  dorothy edith smithborn july 6 1926 is a canadian sociologist with research interests besi\n",
            "  des in sociology in many disciplines including womens studies psychology and educational\n",
            "\n",
            "* Jane Fonda                                         0.93914\n",
            "  jane fonda born lady jayne seymour fonda december 21 1937 is an american actress writer po\n",
            "  litical activist former fashion model and fitness guru she is\n",
            "==========================================================\n",
            "Cluster 6    \n",
            "album:0.053\n",
            "band:0.044\n",
            "music:0.042\n",
            "released:0.028\n",
            "jazz:0.023\n",
            "\n",
            "\n",
            "* Will.i.am                                          0.95336\n",
            "  william adams born march 15 1975 known by his stage name william pronounced will i am is a\n",
            "  n american rapper songwriter entrepreneur actor dj record\n",
            "\n",
            "* Tony Mills (musician)                              0.95359\n",
            "  tony mills born 7 july 1962 in solihull england is an english rock singer best known for h\n",
            "  is work with shy and tnthailing from birmingham\n",
            "\n",
            "* Keith Urban                                        0.95379\n",
            "  keith lionel urban born 26 october 1967 is a new zealand born australian country music sin\n",
            "  ger songwriter guitarist entrepreneur and music competition judge in 1991\n",
            "\n",
            "* Prince (musician)                                  0.95420\n",
            "  prince rogers nelson born june 7 1958 known by his mononym prince is an american singerson\n",
            "  gwriter multiinstrumentalist and actor he has produced ten platinum albums\n",
            "\n",
            "* Steve Overland                                     0.95503\n",
            "  steve overland is a british singermusician who was the lead vocalist and songwriter for th\n",
            "  e bands wildlife fm the ladder shadowman and his own group\n",
            "\n",
            "* Jesse Johnson (musician)                           0.95606\n",
            "  jesse woods johnson born june 1 1960 in rock island illinois is a musician best known as t\n",
            "  he guitarist in the original lineup of the\n",
            "\n",
            "* Mark Cross (musician)                              0.95609\n",
            "  mark cross born 2 august 1965 london is a hard rock and heavy metal drummer he was born to\n",
            "   an english father and german mother\n",
            "\n",
            "* Stewart Levine                                     0.95888\n",
            "  stewart levine is an american record producer he has worked with such artists as the crusa\n",
            "  ders minnie riperton lionel richie simply red hugh masekela dr\n",
            "==========================================================\n",
            "Cluster 7    \n",
            "law:0.133\n",
            "court:0.081\n",
            "judge:0.060\n",
            "district:0.042\n",
            "justice:0.040\n",
            "\n",
            "\n",
            "* Barry Sullivan (lawyer)                            0.89228\n",
            "  barry sullivan is a chicago lawyer and as of july 1 2009 the cooney conway chair in advoca\n",
            "  cy at loyola university chicago school of law\n",
            "\n",
            "* William G. Young                                   0.89433\n",
            "  william glover young born 1940 is a united states federal judge for the district of massac\n",
            "  husetts young was born in huntington new york he attended\n",
            "\n",
            "* Bernard Bell (attorney)                            0.89617\n",
            "  bernard bell is the associate dean for academic affairs and faculty professor of law and h\n",
            "  erbert hannoch scholar at rutgers school of lawnewark bell received\n",
            "\n",
            "* George B. Daniels                                  0.89796\n",
            "  george benjamin daniels born 1953 is a united states federal judge for the united states d\n",
            "  istrict court for the southern district of new yorkdaniels was\n",
            "\n",
            "* Robinson O. Everett                                0.90432\n",
            "  robinson o everett march 18 1928 june 12 2009 was an american lawyer judge and a professor\n",
            "   of law at duke universityeverett was born in\n",
            "\n",
            "* James G. Carr                                      0.90595\n",
            "  james g carr born july 7 1940 is a federal district judge for the united states district c\n",
            "  ourt for the northern district of ohiocarr was\n",
            "\n",
            "* John C. Eastman                                    0.90764\n",
            "  john c eastman born april 21 1960 is a conservative american law professor and constitutio\n",
            "  nal law scholar he is the henry salvatori professor of law\n",
            "\n",
            "* Jean Constance Hamilton                            0.90830\n",
            "  jean constance hamilton born 1945 is a senior united states district judge of the united s\n",
            "  tates district court for the eastern district of missouriborn in\n",
            "==========================================================\n",
            "Cluster 8    \n",
            "football:0.050\n",
            "league:0.044\n",
            "club:0.043\n",
            "season:0.042\n",
            "played:0.037\n",
            "\n",
            "\n",
            "* Chris Day                                          0.93686\n",
            "  christopher nicholas chris day born 28 july 1975 is an english professional footballer who\n",
            "   plays as a goalkeeper for stevenageday started his career at tottenham\n",
            "\n",
            "* Jason Roberts (footballer)                         0.93775\n",
            "  jason andre davis roberts mbe born 25 january 1978 is a former professional footballer and\n",
            "   now a football punditborn in park royal london roberts was\n",
            "\n",
            "* Tony Smith (footballer, born 1957)                 0.93839\n",
            "  anthony tony smith born 20 february 1957 is a former footballer who played as a central de\n",
            "  fender in the football league in the 1970s and\n",
            "\n",
            "* Neil Grayson                                       0.94178\n",
            "  neil grayson born 1 november 1964 in york is an english footballer who last played as a st\n",
            "  riker for sutton towngraysons first club was local\n",
            "\n",
            "* Richard Ambrose                                    0.94220\n",
            "  richard ambrose born 10 june 1972 is a former australian rules footballer who played with \n",
            "  the sydney swans in the australian football league afl he\n",
            "\n",
            "* Paul Robinson (footballer, born 1979)              0.94245\n",
            "  paul william robinson born 15 october 1979 is an english professional footballer who plays\n",
            "   for blackburn rovers as a goalkeeper he is a former england\n",
            "\n",
            "* Alex Lawless                                       0.94269\n",
            "  alexander graham alex lawless born 26 march 1985 is a welsh professional footballer who pl\n",
            "  ays for luton town as a midfielderlawless began his career with\n",
            "\n",
            "* Sol Campbell                                       0.94275\n",
            "  sulzeer jeremiah sol campbell born 18 september 1974 is a former england international foo\n",
            "  tballer a central defender he had a 19year career playing in the\n",
            "==========================================================\n",
            "Cluster 9    \n",
            "research:0.038\n",
            "university:0.035\n",
            "professor:0.030\n",
            "science:0.023\n",
            "institute:0.019\n",
            "\n",
            "\n",
            "* Lawrence W. Green                                  0.95862\n",
            "  lawrence w green is best known by health education researchers as the originator of the pr\n",
            "  ecede model and codeveloper of the precedeproceed model which has\n",
            "\n",
            "* Timothy Luke                                       0.96028\n",
            "  timothy w luke is university distinguished professor of political science in the college o\n",
            "  f liberal arts and human sciences as well as program chair of\n",
            "\n",
            "* Ren%C3%A9e Fox                                     0.96119\n",
            "  rene c fox a summa cum laude graduate of smith college in 1949 earned her phd in sociology\n",
            "   in 1954 from radcliffe college harvard university\n",
            "\n",
            "* Francis Gavin                                      0.96213\n",
            "  francis j gavin is first frank stanton chair in nuclear security policy studies and profes\n",
            "  sor of political science at mit before joining mit he was\n",
            "\n",
            "* Catherine Hakim                                    0.96358\n",
            "  catherine hakim born 30 may 1948 is a british sociologist who specialises in womens employ\n",
            "  ment and womens issues she is currently a professorial research fellow\n",
            "\n",
            "* Daniel Berg (educator)                             0.96361\n",
            "  daniel berg is a scientist educator and was the fifteenth president of rensselaer polytech\n",
            "  nic institutehe was born on june 1 1929 in new york city\n",
            "\n",
            "* Georg von Krogh                                    0.96375\n",
            "  georg von krogh was born in oslo norway he is a professor at eth zurich and holds the chai\n",
            "  r of strategic management and innovation he\n",
            "\n",
            "* Martin Apple                                       0.96381\n",
            "  martin a apple is president of the council of scientific society presidents cssp an organi\n",
            "  zation of presidents of some sixty scientific federations and societies whose\n",
            "==========================================================\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BbaNecYIJosf"
      },
      "source": [
        "Clusters 0, 1, and 5 appear to be still mixed, but others are quite consistent in content.\n",
        "* Cluster 0: artists, book, him/his\n",
        "* Cluster 1: film, theatre, films, tv, actor \n",
        "* Cluster 2: baseball players\n",
        "* Cluster 3: elections, ministers\n",
        "* Cluster 4: music, orchestra, symphony \n",
        "* Cluster 5: female figures from various fields\n",
        "* Cluster 6: composers, songwriters, singers, music producers\n",
        "* Cluster 7: law, courts, justice \n",
        "* Cluster 8: football \n",
        "* Cluster 9: academia\n",
        "\n",
        "Clusters are now more pure, but some are qualitatively \"bigger\" than others. For instance, the category of scholars is more general than the category of baseball players. Increasing the number of clusters may split larger clusters. Another way to look at the size of the clusters is to count the number of articles in each cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ZmhOFxp8Josg",
        "outputId": "d731494f-03fa-43c5-cbc0-6d70315ad0e0"
      },
      "source": [
        "np.bincount(cluster_assignment[10])"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([19618,  3857,  4173,  5219,  1743,  6900,  5530,  1348,  4384,\n",
              "        6299])"
            ]
          },
          "metadata": {},
          "execution_count": 53
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h_mWs6QfJosg"
      },
      "source": [
        "**Quiz Question**. Which of the 10 clusters above contains the greatest number of articles?\n",
        "\n",
        "1. * Cluster 0: artists, book, him/his\n",
        "2. * Cluster 4: music, orchestra, symphony \n",
        "3. * Cluster 5: female figures from various fields\n",
        "4. * Cluster 7: law, courts, justice \n",
        "5. * Cluster 9: academia"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "eaFF1-tvTnz_",
        "outputId": "cfc84db3-aebb-49a6-d9e2-84dab749922e"
      },
      "source": [
        "np.argmax(np.bincount(cluster_assignment[10]))\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0"
            ]
          },
          "metadata": {},
          "execution_count": 63
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ax2ZHsU7Josh"
      },
      "source": [
        "**Quiz Question**. Which of the 10 clusters contains the least number of articles?\n",
        "\n",
        "1. * Cluster 1: film, theatre, films, tv, actor \n",
        "2. * Cluster 3: elections, ministers\n",
        "3. * Cluster 6: composers, songwriters, singers, music producers\n",
        "4. * Cluster 7: law, courts, justice \n",
        "5. * Cluster 8: football "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "t9rfEvS5VQQu",
        "outputId": "8d8b25fb-7a6a-46bf-92ef-51b555fc30ea"
      },
      "source": [
        "np.argmin(np.bincount(cluster_assignment[10]))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "7"
            ]
          },
          "metadata": {},
          "execution_count": 62
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZhjnPO2_Josh"
      },
      "source": [
        "There appears to be at least some connection between the topical consistency of a cluster and the number of its member data points."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t0eEv9gxJosi"
      },
      "source": [
        "Let us visualize the case for K=25. For the sake of brevity, we do not print the content of documents. It turns out that the top words with highest TF-IDF weights in each cluster are representative of the cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": false,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "aGyd-20fJosi",
        "outputId": "f2f103c2-6edc-4cb2-ed07-2d10bc28407a"
      },
      "source": [
        "visualize_document_clusters(wiki, tf_idf, centroids[25], cluster_assignment[25], 25,\n",
        "                            map_index_to_word, display_content=False) # turn off text for brevity"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "==========================================================\n",
            "Cluster 0    \n",
            "poetry:0.053\n",
            "novel:0.043\n",
            "book:0.042\n",
            "published:0.039\n",
            "fiction:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 1    \n",
            "film:0.100\n",
            "theatre:0.039\n",
            "films:0.036\n",
            "directed:0.029\n",
            "actor:0.028\n",
            "\n",
            "==========================================================\n",
            "Cluster 2    \n",
            "law:0.143\n",
            "court:0.087\n",
            "judge:0.066\n",
            "district:0.045\n",
            "justice:0.042\n",
            "\n",
            "==========================================================\n",
            "Cluster 3    \n",
            "republican:0.061\n",
            "senate:0.050\n",
            "district:0.044\n",
            "state:0.039\n",
            "democratic:0.037\n",
            "\n",
            "==========================================================\n",
            "Cluster 4    \n",
            "music:0.114\n",
            "piano:0.047\n",
            "orchestra:0.039\n",
            "composition:0.038\n",
            "composer:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 5    \n",
            "album:0.116\n",
            "released:0.058\n",
            "her:0.056\n",
            "single:0.046\n",
            "music:0.040\n",
            "\n",
            "==========================================================\n",
            "Cluster 6    \n",
            "music:0.055\n",
            "jazz:0.038\n",
            "album:0.028\n",
            "song:0.020\n",
            "records:0.019\n",
            "\n",
            "==========================================================\n",
            "Cluster 7    \n",
            "board:0.028\n",
            "business:0.027\n",
            "economics:0.026\n",
            "chairman:0.025\n",
            "president:0.025\n",
            "\n",
            "==========================================================\n",
            "Cluster 8    \n",
            "he:0.011\n",
            "his:0.009\n",
            "that:0.009\n",
            "world:0.007\n",
            "book:0.007\n",
            "\n",
            "==========================================================\n",
            "Cluster 9    \n",
            "research:0.050\n",
            "university:0.039\n",
            "professor:0.038\n",
            "science:0.030\n",
            "institute:0.021\n",
            "\n",
            "==========================================================\n",
            "Cluster 10    \n",
            "foreign:0.075\n",
            "ambassador:0.063\n",
            "affairs:0.057\n",
            "security:0.044\n",
            "nations:0.042\n",
            "\n",
            "==========================================================\n",
            "Cluster 11    \n",
            "baseball:0.110\n",
            "league:0.103\n",
            "major:0.052\n",
            "games:0.047\n",
            "season:0.045\n",
            "\n",
            "==========================================================\n",
            "Cluster 12    \n",
            "art:0.146\n",
            "museum:0.078\n",
            "gallery:0.057\n",
            "artist:0.033\n",
            "arts:0.032\n",
            "\n",
            "==========================================================\n",
            "Cluster 13    \n",
            "air:0.028\n",
            "military:0.027\n",
            "police:0.024\n",
            "force:0.023\n",
            "commander:0.022\n",
            "\n",
            "==========================================================\n",
            "Cluster 14    \n",
            "party:0.064\n",
            "minister:0.063\n",
            "election:0.054\n",
            "parliament:0.031\n",
            "elected:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 15    \n",
            "radio:0.072\n",
            "show:0.052\n",
            "news:0.051\n",
            "bbc:0.033\n",
            "television:0.030\n",
            "\n",
            "==========================================================\n",
            "Cluster 16    \n",
            "church:0.120\n",
            "bishop:0.091\n",
            "diocese:0.044\n",
            "lds:0.044\n",
            "archbishop:0.043\n",
            "\n",
            "==========================================================\n",
            "Cluster 17    \n",
            "opera:0.212\n",
            "ballet:0.088\n",
            "she:0.061\n",
            "la:0.035\n",
            "her:0.033\n",
            "\n",
            "==========================================================\n",
            "Cluster 18    \n",
            "orchestra:0.203\n",
            "symphony:0.146\n",
            "conductor:0.107\n",
            "philharmonic:0.077\n",
            "music:0.076\n",
            "\n",
            "==========================================================\n",
            "Cluster 19    \n",
            "she:0.146\n",
            "her:0.092\n",
            "miss:0.017\n",
            "actress:0.015\n",
            "women:0.011\n",
            "\n",
            "==========================================================\n",
            "Cluster 20    \n",
            "racing:0.127\n",
            "formula:0.078\n",
            "race:0.066\n",
            "car:0.060\n",
            "driver:0.054\n",
            "\n",
            "==========================================================\n",
            "Cluster 21    \n",
            "championships:0.057\n",
            "tour:0.055\n",
            "pga:0.041\n",
            "olympics:0.035\n",
            "metres:0.035\n",
            "\n",
            "==========================================================\n",
            "Cluster 22    \n",
            "league:0.052\n",
            "rugby:0.049\n",
            "club:0.046\n",
            "cup:0.045\n",
            "season:0.041\n",
            "\n",
            "==========================================================\n",
            "Cluster 23    \n",
            "band:0.104\n",
            "album:0.049\n",
            "rock:0.031\n",
            "guitar:0.031\n",
            "bands:0.030\n",
            "\n",
            "==========================================================\n",
            "Cluster 24    \n",
            "football:0.057\n",
            "coach:0.053\n",
            "season:0.047\n",
            "basketball:0.042\n",
            "played:0.039\n",
            "\n",
            "==========================================================\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "By80g2QJJosj"
      },
      "source": [
        "Looking at the representative examples and top words, we classify each cluster as follows.\n",
        "\n",
        "* Cluster 0: Literature\n",
        "* Cluster 1: Film and theater\n",
        "* Cluster 2: Law\n",
        "* Cluster 3: Politics\n",
        "* Cluster 4: Classical music\n",
        "* Cluster 5: Popular music\n",
        "* Cluster 6: Jazz music\n",
        "* Cluster 7: Business and economics\n",
        "* Cluster 8: (mixed; no clear theme)\n",
        "* Cluster 9: Academia and research\n",
        "* Cluster 10: International affairs\n",
        "* Cluster 11: Baseball\n",
        "* Cluster 12: Art\n",
        "* Cluster 13: Military\n",
        "* Cluster 14: Politics\n",
        "* Cluster 15: Radio and TV\n",
        "* Cluster 16: Catholic church\n",
        "* Cluster 17: Opera and ballet\n",
        "* Cluster 18: Orchestra music\n",
        "* Cluster 19: Females from various fields\n",
        "* Cluster 20: Car racing\n",
        "* Cluster 21: General sports\n",
        "* Cluster 22: Rugby\n",
        "* Cluster 23: Rock music\n",
        "* Cluster 24: Team sports\n",
        "\n",
        "Indeed, increasing K achieved the desired effect of breaking up large clusters.  Depending on the application, this may or may not be preferable to the K=10 analysis.\n",
        "\n",
        "Let's take it to the extreme and set K=100. We have a suspicion that this value is too large. Let us look at the top words from each cluster:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8nUCmTIZJosk",
        "outputId": "52ed38d2-a334-41c8-e807-160c0cef7b1a"
      },
      "source": [
        "k=100\n",
        "visualize_document_clusters(wiki, tf_idf, centroids[k], cluster_assignment[k], k,\n",
        "                            map_index_to_word, display_content=False)\n",
        "# turn off text for brevity -- turn it on if you are curious ;)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "==========================================================\n",
            "Cluster 0    \n",
            "psychology:0.195\n",
            "psychological:0.066\n",
            "research:0.057\n",
            "psychologist:0.045\n",
            "cognitive:0.041\n",
            "\n",
            "==========================================================\n",
            "Cluster 1    \n",
            "film:0.213\n",
            "festival:0.060\n",
            "films:0.054\n",
            "directed:0.039\n",
            "feature:0.037\n",
            "\n",
            "==========================================================\n",
            "Cluster 2    \n",
            "law:0.146\n",
            "court:0.097\n",
            "judge:0.074\n",
            "district:0.051\n",
            "justice:0.045\n",
            "\n",
            "==========================================================\n",
            "Cluster 3    \n",
            "mayor:0.137\n",
            "city:0.049\n",
            "council:0.040\n",
            "elected:0.032\n",
            "election:0.030\n",
            "\n",
            "==========================================================\n",
            "Cluster 4    \n",
            "music:0.147\n",
            "composition:0.048\n",
            "composer:0.047\n",
            "orchestra:0.026\n",
            "composers:0.025\n",
            "\n",
            "==========================================================\n",
            "Cluster 5    \n",
            "album:0.108\n",
            "her:0.078\n",
            "billboard:0.075\n",
            "chart:0.072\n",
            "singles:0.067\n",
            "\n",
            "==========================================================\n",
            "Cluster 6    \n",
            "music:0.052\n",
            "songs:0.025\n",
            "records:0.023\n",
            "song:0.022\n",
            "album:0.022\n",
            "\n",
            "==========================================================\n",
            "Cluster 7    \n",
            "chairman:0.057\n",
            "board:0.048\n",
            "president:0.035\n",
            "executive:0.033\n",
            "ceo:0.025\n",
            "\n",
            "==========================================================\n",
            "Cluster 8    \n",
            "german:0.120\n",
            "germany:0.042\n",
            "der:0.030\n",
            "berlin:0.025\n",
            "die:0.017\n",
            "\n",
            "==========================================================\n",
            "Cluster 9    \n",
            "india:0.092\n",
            "indian:0.084\n",
            "sabha:0.038\n",
            "lok:0.033\n",
            "singh:0.028\n",
            "\n",
            "==========================================================\n",
            "Cluster 10    \n",
            "czech:0.207\n",
            "prague:0.124\n",
            "republic:0.046\n",
            "czechoslovakia:0.032\n",
            "vclav:0.021\n",
            "\n",
            "==========================================================\n",
            "Cluster 11    \n",
            "soccer:0.294\n",
            "league:0.071\n",
            "indoor:0.063\n",
            "team:0.055\n",
            "season:0.052\n",
            "\n",
            "==========================================================\n",
            "Cluster 12    \n",
            "novel:0.096\n",
            "fiction:0.079\n",
            "published:0.044\n",
            "stories:0.043\n",
            "short:0.039\n",
            "\n",
            "==========================================================\n",
            "Cluster 13    \n",
            "prison:0.036\n",
            "police:0.034\n",
            "sentenced:0.026\n",
            "court:0.025\n",
            "convicted:0.024\n",
            "\n",
            "==========================================================\n",
            "Cluster 14    \n",
            "labor:0.104\n",
            "australian:0.097\n",
            "liberal:0.073\n",
            "election:0.067\n",
            "minister:0.064\n",
            "\n",
            "==========================================================\n",
            "Cluster 15    \n",
            "radio:0.112\n",
            "show:0.069\n",
            "host:0.041\n",
            "station:0.035\n",
            "sports:0.026\n",
            "\n",
            "==========================================================\n",
            "Cluster 16    \n",
            "bishop:0.147\n",
            "church:0.084\n",
            "diocese:0.075\n",
            "archbishop:0.073\n",
            "ordained:0.055\n",
            "\n",
            "==========================================================\n",
            "Cluster 17    \n",
            "de:0.118\n",
            "la:0.050\n",
            "french:0.026\n",
            "el:0.021\n",
            "paris:0.018\n",
            "\n",
            "==========================================================\n",
            "Cluster 18    \n",
            "clarinet:0.087\n",
            "bass:0.086\n",
            "saxophone:0.079\n",
            "flute:0.077\n",
            "music:0.062\n",
            "\n",
            "==========================================================\n",
            "Cluster 19    \n",
            "book:0.045\n",
            "books:0.032\n",
            "published:0.027\n",
            "editor:0.025\n",
            "magazine:0.021\n",
            "\n",
            "==========================================================\n",
            "Cluster 20    \n",
            "racing:0.135\n",
            "nascar:0.104\n",
            "car:0.094\n",
            "race:0.077\n",
            "series:0.075\n",
            "\n",
            "==========================================================\n",
            "Cluster 21    \n",
            "tour:0.259\n",
            "pga:0.216\n",
            "golf:0.139\n",
            "open:0.073\n",
            "golfer:0.062\n",
            "\n",
            "==========================================================\n",
            "Cluster 22    \n",
            "league:0.089\n",
            "town:0.063\n",
            "season:0.061\n",
            "club:0.059\n",
            "football:0.054\n",
            "\n",
            "==========================================================\n",
            "Cluster 23    \n",
            "album:0.119\n",
            "released:0.058\n",
            "music:0.034\n",
            "records:0.027\n",
            "single:0.026\n",
            "\n",
            "==========================================================\n",
            "Cluster 24    \n",
            "football:0.054\n",
            "cup:0.048\n",
            "club:0.045\n",
            "team:0.039\n",
            "league:0.036\n",
            "\n",
            "==========================================================\n",
            "Cluster 25    \n",
            "league:0.096\n",
            "era:0.093\n",
            "baseball:0.089\n",
            "innings:0.086\n",
            "pitcher:0.085\n",
            "\n",
            "==========================================================\n",
            "Cluster 26    \n",
            "air:0.373\n",
            "force:0.241\n",
            "command:0.105\n",
            "commander:0.094\n",
            "base:0.080\n",
            "\n",
            "==========================================================\n",
            "Cluster 27    \n",
            "physics:0.173\n",
            "quantum:0.059\n",
            "theoretical:0.045\n",
            "research:0.043\n",
            "theory:0.039\n",
            "\n",
            "==========================================================\n",
            "Cluster 28    \n",
            "sierra:0.276\n",
            "leone:0.220\n",
            "koroma:0.061\n",
            "freetown:0.056\n",
            "leonean:0.046\n",
            "\n",
            "==========================================================\n",
            "Cluster 29    \n",
            "russian:0.172\n",
            "soviet:0.067\n",
            "russia:0.056\n",
            "moscow:0.056\n",
            "vladimir:0.022\n",
            "\n",
            "==========================================================\n",
            "Cluster 30    \n",
            "harris:0.398\n",
            "university:0.014\n",
            "alabama:0.013\n",
            "state:0.013\n",
            "he:0.012\n",
            "\n",
            "==========================================================\n",
            "Cluster 31    \n",
            "theatre:0.194\n",
            "directed:0.033\n",
            "production:0.031\n",
            "she:0.029\n",
            "play:0.029\n",
            "\n",
            "==========================================================\n",
            "Cluster 32    \n",
            "linguistics:0.164\n",
            "language:0.147\n",
            "linguistic:0.064\n",
            "languages:0.044\n",
            "research:0.036\n",
            "\n",
            "==========================================================\n",
            "Cluster 33    \n",
            "economics:0.150\n",
            "economic:0.098\n",
            "economist:0.053\n",
            "policy:0.047\n",
            "research:0.045\n",
            "\n",
            "==========================================================\n",
            "Cluster 34    \n",
            "news:0.131\n",
            "anchor:0.068\n",
            "reporter:0.057\n",
            "she:0.044\n",
            "correspondent:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 35    \n",
            "rights:0.172\n",
            "human:0.129\n",
            "law:0.056\n",
            "she:0.036\n",
            "civil:0.029\n",
            "\n",
            "==========================================================\n",
            "Cluster 36    \n",
            "foreign:0.062\n",
            "ambassador:0.056\n",
            "affairs:0.051\n",
            "security:0.044\n",
            "secretary:0.043\n",
            "\n",
            "==========================================================\n",
            "Cluster 37    \n",
            "mathematics:0.151\n",
            "mathematical:0.113\n",
            "theory:0.055\n",
            "professor:0.047\n",
            "mathematician:0.047\n",
            "\n",
            "==========================================================\n",
            "Cluster 38    \n",
            "mexico:0.132\n",
            "mexican:0.118\n",
            "de:0.040\n",
            "pri:0.031\n",
            "cartel:0.026\n",
            "\n",
            "==========================================================\n",
            "Cluster 39    \n",
            "film:0.141\n",
            "documentary:0.087\n",
            "films:0.056\n",
            "festival:0.041\n",
            "cinema:0.033\n",
            "\n",
            "==========================================================\n",
            "Cluster 40    \n",
            "hong:0.283\n",
            "kong:0.267\n",
            "chinese:0.067\n",
            "china:0.038\n",
            "wong:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 41    \n",
            "actor:0.049\n",
            "role:0.047\n",
            "film:0.043\n",
            "series:0.042\n",
            "appeared:0.037\n",
            "\n",
            "==========================================================\n",
            "Cluster 42    \n",
            "football:0.120\n",
            "afl:0.113\n",
            "australian:0.083\n",
            "season:0.058\n",
            "club:0.058\n",
            "\n",
            "==========================================================\n",
            "Cluster 43    \n",
            "band:0.116\n",
            "album:0.046\n",
            "bands:0.035\n",
            "guitar:0.033\n",
            "rock:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 44    \n",
            "puerto:0.308\n",
            "rico:0.217\n",
            "rican:0.068\n",
            "juan:0.042\n",
            "ricos:0.032\n",
            "\n",
            "==========================================================\n",
            "Cluster 45    \n",
            "giants:0.257\n",
            "baseball:0.104\n",
            "league:0.085\n",
            "francisco:0.065\n",
            "san:0.064\n",
            "\n",
            "==========================================================\n",
            "Cluster 46    \n",
            "racing:0.104\n",
            "jockey:0.078\n",
            "race:0.069\n",
            "stakes:0.064\n",
            "horse:0.053\n",
            "\n",
            "==========================================================\n",
            "Cluster 47    \n",
            "archaeology:0.281\n",
            "archaeological:0.114\n",
            "ancient:0.072\n",
            "archaeologist:0.059\n",
            "excavations:0.055\n",
            "\n",
            "==========================================================\n",
            "Cluster 48    \n",
            "art:0.082\n",
            "artist:0.032\n",
            "gallery:0.031\n",
            "painting:0.028\n",
            "paintings:0.028\n",
            "\n",
            "==========================================================\n",
            "Cluster 49    \n",
            "she:0.136\n",
            "her:0.120\n",
            "actress:0.024\n",
            "film:0.018\n",
            "television:0.013\n",
            "\n",
            "==========================================================\n",
            "Cluster 50    \n",
            "comics:0.191\n",
            "comic:0.122\n",
            "strip:0.039\n",
            "graphic:0.036\n",
            "book:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 51    \n",
            "comedy:0.165\n",
            "show:0.068\n",
            "comedian:0.060\n",
            "standup:0.050\n",
            "series:0.030\n",
            "\n",
            "==========================================================\n",
            "Cluster 52    \n",
            "formula:0.167\n",
            "racing:0.129\n",
            "car:0.079\n",
            "driver:0.078\n",
            "championship:0.074\n",
            "\n",
            "==========================================================\n",
            "Cluster 53    \n",
            "church:0.186\n",
            "lds:0.185\n",
            "churchs:0.094\n",
            "latterday:0.070\n",
            "byu:0.068\n",
            "\n",
            "==========================================================\n",
            "Cluster 54    \n",
            "design:0.169\n",
            "architecture:0.121\n",
            "architectural:0.058\n",
            "architects:0.038\n",
            "architect:0.038\n",
            "\n",
            "==========================================================\n",
            "Cluster 55    \n",
            "university:0.048\n",
            "philosophy:0.042\n",
            "professor:0.041\n",
            "studies:0.038\n",
            "history:0.037\n",
            "\n",
            "==========================================================\n",
            "Cluster 56    \n",
            "food:0.260\n",
            "cooking:0.049\n",
            "she:0.038\n",
            "cookbook:0.031\n",
            "culinary:0.028\n",
            "\n",
            "==========================================================\n",
            "Cluster 57    \n",
            "oklahoma:0.212\n",
            "oregon:0.169\n",
            "portland:0.040\n",
            "district:0.032\n",
            "law:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 58    \n",
            "piano:0.093\n",
            "music:0.071\n",
            "orchestra:0.065\n",
            "chamber:0.045\n",
            "symphony:0.040\n",
            "\n",
            "==========================================================\n",
            "Cluster 59    \n",
            "iraqi:0.160\n",
            "iraq:0.150\n",
            "baghdad:0.060\n",
            "saddam:0.044\n",
            "hussein:0.035\n",
            "\n",
            "==========================================================\n",
            "Cluster 60    \n",
            "business:0.034\n",
            "company:0.023\n",
            "technology:0.023\n",
            "management:0.023\n",
            "global:0.019\n",
            "\n",
            "==========================================================\n",
            "Cluster 61    \n",
            "league:0.120\n",
            "baseball:0.108\n",
            "major:0.058\n",
            "minor:0.057\n",
            "season:0.042\n",
            "\n",
            "==========================================================\n",
            "Cluster 62    \n",
            "freestyle:0.154\n",
            "swimming:0.124\n",
            "m:0.117\n",
            "swimmer:0.090\n",
            "heat:0.074\n",
            "\n",
            "==========================================================\n",
            "Cluster 63    \n",
            "song:0.162\n",
            "eurovision:0.112\n",
            "contest:0.073\n",
            "she:0.057\n",
            "her:0.038\n",
            "\n",
            "==========================================================\n",
            "Cluster 64    \n",
            "bbc:0.235\n",
            "radio:0.119\n",
            "news:0.053\n",
            "presenter:0.053\n",
            "she:0.051\n",
            "\n",
            "==========================================================\n",
            "Cluster 65    \n",
            "army:0.078\n",
            "command:0.078\n",
            "commander:0.078\n",
            "military:0.074\n",
            "staff:0.059\n",
            "\n",
            "==========================================================\n",
            "Cluster 66    \n",
            "jazz:0.214\n",
            "music:0.047\n",
            "band:0.035\n",
            "pianist:0.026\n",
            "trio:0.024\n",
            "\n",
            "==========================================================\n",
            "Cluster 67    \n",
            "chef:0.195\n",
            "restaurant:0.130\n",
            "wine:0.102\n",
            "cooking:0.060\n",
            "food:0.059\n",
            "\n",
            "==========================================================\n",
            "Cluster 68    \n",
            "turkish:0.176\n",
            "turkey:0.104\n",
            "istanbul:0.072\n",
            "ankara:0.029\n",
            "she:0.026\n",
            "\n",
            "==========================================================\n",
            "Cluster 69    \n",
            "thai:0.158\n",
            "cpc:0.063\n",
            "china:0.060\n",
            "party:0.057\n",
            "thailand:0.055\n",
            "\n",
            "==========================================================\n",
            "Cluster 70    \n",
            "miss:0.358\n",
            "pageant:0.206\n",
            "usa:0.122\n",
            "she:0.111\n",
            "her:0.063\n",
            "\n",
            "==========================================================\n",
            "Cluster 71    \n",
            "republican:0.070\n",
            "senate:0.054\n",
            "district:0.048\n",
            "state:0.041\n",
            "house:0.039\n",
            "\n",
            "==========================================================\n",
            "Cluster 72    \n",
            "wrestling:0.069\n",
            "championship:0.046\n",
            "tennis:0.045\n",
            "doubles:0.039\n",
            "champion:0.036\n",
            "\n",
            "==========================================================\n",
            "Cluster 73    \n",
            "poetry:0.214\n",
            "poems:0.081\n",
            "poet:0.069\n",
            "poets:0.043\n",
            "literary:0.041\n",
            "\n",
            "==========================================================\n",
            "Cluster 74    \n",
            "baseball:0.235\n",
            "league:0.077\n",
            "she:0.055\n",
            "aagpbl:0.045\n",
            "allamerican:0.039\n",
            "\n",
            "==========================================================\n",
            "Cluster 75    \n",
            "poker:0.477\n",
            "wsop:0.121\n",
            "event:0.091\n",
            "limit:0.078\n",
            "winnings:0.072\n",
            "\n",
            "==========================================================\n",
            "Cluster 76    \n",
            "rugby:0.195\n",
            "cup:0.049\n",
            "against:0.046\n",
            "played:0.044\n",
            "wales:0.039\n",
            "\n",
            "==========================================================\n",
            "Cluster 77    \n",
            "he:0.010\n",
            "that:0.009\n",
            "his:0.009\n",
            "it:0.007\n",
            "has:0.006\n",
            "\n",
            "==========================================================\n",
            "Cluster 78    \n",
            "blues:0.121\n",
            "drummer:0.066\n",
            "band:0.054\n",
            "rock:0.041\n",
            "drum:0.037\n",
            "\n",
            "==========================================================\n",
            "Cluster 79    \n",
            "marathon:0.064\n",
            "olympics:0.059\n",
            "championships:0.057\n",
            "olympic:0.055\n",
            "she:0.041\n",
            "\n",
            "==========================================================\n",
            "Cluster 80    \n",
            "hockey:0.218\n",
            "nhl:0.136\n",
            "ice:0.066\n",
            "season:0.053\n",
            "league:0.048\n",
            "\n",
            "==========================================================\n",
            "Cluster 81    \n",
            "party:0.066\n",
            "minister:0.061\n",
            "election:0.044\n",
            "parliament:0.036\n",
            "elected:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 82    \n",
            "computer:0.096\n",
            "engineering:0.074\n",
            "research:0.047\n",
            "science:0.046\n",
            "systems:0.038\n",
            "\n",
            "==========================================================\n",
            "Cluster 83    \n",
            "election:0.083\n",
            "manitoba:0.071\n",
            "minister:0.069\n",
            "liberal:0.068\n",
            "canadian:0.055\n",
            "\n",
            "==========================================================\n",
            "Cluster 84    \n",
            "orchestra:0.221\n",
            "symphony:0.156\n",
            "conductor:0.132\n",
            "music:0.081\n",
            "philharmonic:0.080\n",
            "\n",
            "==========================================================\n",
            "Cluster 85    \n",
            "sri:0.288\n",
            "lanka:0.187\n",
            "lankan:0.098\n",
            "colombo:0.048\n",
            "ceylon:0.029\n",
            "\n",
            "==========================================================\n",
            "Cluster 86    \n",
            "basketball:0.157\n",
            "nba:0.086\n",
            "coach:0.077\n",
            "points:0.049\n",
            "season:0.042\n",
            "\n",
            "==========================================================\n",
            "Cluster 87    \n",
            "cricket:0.194\n",
            "firstclass:0.114\n",
            "cricketer:0.073\n",
            "batsman:0.069\n",
            "wickets:0.061\n",
            "\n",
            "==========================================================\n",
            "Cluster 88    \n",
            "runs:0.116\n",
            "league:0.100\n",
            "baseball:0.089\n",
            "batted:0.067\n",
            "home:0.064\n",
            "\n",
            "==========================================================\n",
            "Cluster 89    \n",
            "she:0.178\n",
            "her:0.053\n",
            "women:0.021\n",
            "member:0.017\n",
            "university:0.017\n",
            "\n",
            "==========================================================\n",
            "Cluster 90    \n",
            "research:0.061\n",
            "medical:0.047\n",
            "medicine:0.046\n",
            "professor:0.033\n",
            "chemistry:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 91    \n",
            "columbia:0.096\n",
            "vancouver:0.091\n",
            "bc:0.075\n",
            "british:0.073\n",
            "canadian:0.071\n",
            "\n",
            "==========================================================\n",
            "Cluster 92    \n",
            "metres:0.179\n",
            "championships:0.145\n",
            "athletics:0.096\n",
            "she:0.079\n",
            "m:0.070\n",
            "\n",
            "==========================================================\n",
            "Cluster 93    \n",
            "health:0.228\n",
            "medical:0.072\n",
            "medicine:0.071\n",
            "care:0.057\n",
            "research:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 94    \n",
            "football:0.118\n",
            "nfl:0.089\n",
            "yards:0.065\n",
            "coach:0.058\n",
            "bowl:0.046\n",
            "\n",
            "==========================================================\n",
            "Cluster 95    \n",
            "jewish:0.219\n",
            "rabbi:0.172\n",
            "israel:0.041\n",
            "yeshiva:0.037\n",
            "hebrew:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 96    \n",
            "film:0.045\n",
            "producer:0.038\n",
            "television:0.037\n",
            "series:0.033\n",
            "directed:0.031\n",
            "\n",
            "==========================================================\n",
            "Cluster 97    \n",
            "chess:0.416\n",
            "grandmaster:0.085\n",
            "olympiad:0.066\n",
            "championship:0.064\n",
            "fide:0.059\n",
            "\n",
            "==========================================================\n",
            "Cluster 98    \n",
            "opera:0.227\n",
            "ballet:0.088\n",
            "she:0.065\n",
            "la:0.036\n",
            "her:0.034\n",
            "\n",
            "==========================================================\n",
            "Cluster 99    \n",
            "art:0.208\n",
            "museum:0.149\n",
            "gallery:0.085\n",
            "arts:0.043\n",
            "contemporary:0.041\n",
            "\n",
            "==========================================================\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zVyMzdedJosl"
      },
      "source": [
        "The class of team sports has been broken into several clusters, soccer (association football) (11, 22, 24), rugby (76), hockey (80), basketball (86), cricket (87), and American football (85).\n",
        "\n",
        "The class of baseball has been broken into San Francisco Giants (45), baseball (61, 74), and baseball stats (88).\n",
        "\n",
        "The class of car racing has been broken into Nascar (20) and Formula 1 (52).\n",
        "\n",
        "**A high value of K encourages pure clusters, but we cannot keep increasing K. For large enough K, related documents end up going to different clusters.**\n",
        "\n",
        "That said, the result for K=100 is not entirely bad. After all, it gives us separate clusters for such categories as Brazil, wrestling, computer science and the Mormon Church. If we set K somewhere between 25 and 100, we should be able to avoid breaking up clusters while discovering new ones.\n",
        "\n",
        "Also, we should ask ourselves how much **granularity** we want in our clustering. If we wanted a rough sketch of Wikipedia, we don't want too detailed clusters. On the other hand, having many clusters can be valuable when we are zooming into a certain part of Wikipedia.\n",
        "\n",
        "**There is no golden rule for choosing K. It all depends on the particular application and domain we are in.**\n",
        "\n",
        "Another heuristic people use that does not rely on so much visualization, which can be hard in many applications (including here!) is as follows.  Track heterogeneity versus K and look for the \"elbow\" of the curve where the heterogeneity decrease rapidly before this value of K, but then only gradually for larger values of K.  This naturally trades off between trying to minimize heterogeneity, but reduce model complexity.  In the heterogeneity versus K plot made above, we did not yet really see a flattening out of the heterogeneity, which might indicate that indeed K=100 is \"reasonable\" and we only see real overfitting for larger values of K (which are even harder to visualize using the methods we attempted above.)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "qmF_4rY6Josl"
      },
      "source": [
        "**Quiz Question**. Another sign of too large K is having lots of small clusters. Look at the distribution of cluster sizes (by number of member data points). How many of the 100 clusters have fewer than 236 articles, i.e. 0.4% of the dataset?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Gz4GEKbpJosm",
        "outputId": "bad2438d-c5d9-4eea-805a-078669efe89e"
      },
      "source": [
        "cluster_assignment_bincount = np.bincount(cluster_assignment[100])\n",
        "len (cluster_assignment_bincount[cluster_assignment_bincount <= 236])"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "33"
            ]
          },
          "metadata": {},
          "execution_count": 70
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1yH2rXgFJosn"
      },
      "source": [
        "### Takeaway\n",
        "\n",
        "Keep in mind though that tiny clusters aren't necessarily bad. A tiny cluster of documents that really look like each others is definitely preferable to a medium-sized cluster of documents with mixed content. However, having too few articles in a cluster may cause overfitting by reading too much into a limited pool of training data."
      ]
    }
  ]
}