2_kmeans-with-text-data_blank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# k-means with text data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": {},
   "source": [
    "## Import necessary packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code block will check if you have the correct version of GraphLab Create. Any version later than 1.8.5 will do. To upgrade, read [this page](https://turi.com/download/upgrade-graphlab-create.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO] graphlab.cython.cy_server: GraphLab Create v2.1 started. Logging: /tmp/graphlab_server_1476898937.log\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This non-commercial license of GraphLab Create for academic use is assigned to rojour@gmail.com and will expire on August 06, 2017.\n"
     ]
    }
   ],
   "source": [
    "import graphlab\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",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "'''Check GraphLab Create version'''\n",
    "from distutils.version import StrictVersion\n",
    "assert (StrictVersion(graphlab.version) >= StrictVersion('1.8.5')), 'GraphLab Create must be version 1.8.5 or later.'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data, extract features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 582,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "wiki = graphlab.SFrame('people_wiki.gl/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 583,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "wiki['tf_idf'] = graphlab.text_analytics.tf_idf(wiki['text'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 586,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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 x[column_name].dtype() == dict, \\\n",
    "        'The chosen column must be dict type, representing sparse data.'\n",
    "        \n",
    "    # Create triples of (row_id, feature_id, count).\n",
    "    # 1. Add a row number.\n",
    "    x = x.add_row_number()\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 words into integers using a OneHotEncoder feature transformation.\n",
    "    f = graphlab.feature_engineering.OneHotEncoder(features=['feature'])\n",
    "    # 1. Fit the transformer using the above data.\n",
    "    f.fit(x)\n",
    "    # 2. The transform takes 'feature' column and adds a new column 'feature_encoding'.\n",
    "    x = f.transform(x)\n",
    "    # 3. Get the feature mapping.\n",
    "    mapping = f['feature_encoding']\n",
    "    # 4. Get the feature id to use for each key.\n",
    "    x['feature_id'] = x['encoded_features'].dict_keys().apply(lambda x: x[0])\n",
    "\n",
    "    # Create numpy arrays that contain the data for the sparse matrix.\n",
    "    i = np.array(x['id'])\n",
    "    j = np.array(x['feature_id'])\n",
    "    v = np.array(x['value'])\n",
    "    width = x['id'].max() + 1\n",
    "    height = x['feature_id'].max() + 1\n",
    "\n",
    "    # Create a sparse matrix.\n",
    "    mat = csr_matrix((v, (i, j)), shape=(width, height))\n",
    "\n",
    "    return mat, mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 587,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The conversion will take about a minute or two.\n",
    "tf_idf, map_index_to_word = sframe_to_scipy(wiki, 'tf_idf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 588,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<59071x547979 sparse matrix of type '<type 'numpy.float64'>'\n",
       "\twith 10379283 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 588,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf_idf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": "markdown",
   "metadata": {},
   "source": [
    "## Normalize all vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 589,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import normalize\n",
    "tf_idf = normalize(tf_idf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implement k-means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 590,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": {},
   "source": [
    "### Assigning clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 591,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "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"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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",
   "execution_count": 592,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Students should write code here\n",
    "centroids = tf_idf[:3,:]\n",
    "dist_to_centroids = pairwise_distances(tf_idf, centroids, metric='euclidean')\n",
    "dist = dist_to_centroids[430][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 593,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pass\n"
     ]
    }
   ],
   "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')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 594,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "closest_cluster = np.arange(len(wiki['tf_idf']))\n",
    "distances = dist_to_centroids\n",
    "# Students should write code here\n",
    "for i in closest_cluster:\n",
    "    closest_cluster[i] = np.argmin(distances[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 595,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pass\n"
     ]
    }
   ],
   "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')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "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",
   "execution_count": 596,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Students should write code here\n",
    "centroids = tf_idf[:3,:]\n",
    "distance = pairwise_distances(tf_idf, centroids, metric='euclidean')\n",
    "cluster_assignment = map(lambda x: np.argmin(distance[x]), range(len(distance)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 597,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pass\n"
     ]
    }
   ],
   "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.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to fill in the blanks in this function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 598,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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')\n",
    "    # Compute cluster assignments for each data point:\n",
    "    # Fill in the blank (RHS only)\n",
    "    cluster_assignment = map(lambda x: np.argmin(distances_from_centroids[x]), range(len(distances_from_centroids)))\n",
    "    \n",
    "    return np.array(cluster_assignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 599,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pass\n"
     ]
    }
   ],
   "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.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Revising clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's turn to Step 2, where we compute the new centroids given the cluster assignments. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "source": [
    "To develop intuition about filtering, let's look at a toy example consisting of 3 data points and 2 clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 600,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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.]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assign these data points to the closest centroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 601,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 1]\n"
     ]
    }
   ],
   "source": [
    "cluster_assignment = assign_clusters(data, centroids)\n",
    "print np.bincount(cluster_assignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 602,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False,  True, False], dtype=bool)"
      ]
     },
     "execution_count": 602,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster_assignment==1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likewise for cluster 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 603,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True, False,  True], dtype=bool)"
      ]
     },
     "execution_count": 603,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster_assignment==0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 604,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.]])"
      ]
     },
     "execution_count": 604,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[cluster_assignment==1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 605,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  0.],\n",
       "       [ 2.,  2.,  0.]])"
      ]
     },
     "execution_count": 605,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[cluster_assignment==0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 606,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.5,  2. ,  0. ])"
      ]
     },
     "execution_count": 606,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[cluster_assignment==0, :].mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to complete this function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 607,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def revise_centroids(data, k, cluster_assignment):\n",
    "    new_centroids = []\n",
    "    for i in xrange(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, :]\n",
    "        # Compute the mean of the data points. Fill in the blank (RHS only)\n",
    "        centroid = member_data_points.mean(axis=0)\n",
    "        # Convert numpy.matrix type to numpy.ndarray type\n",
    "        centroid = centroid.A1\n",
    "        new_centroids.append(centroid)\n",
    "        \n",
    "    new_centroids = np.array(new_centroids)\n",
    "\n",
    "    return new_centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 608,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pass\n"
     ]
    }
   ],
   "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')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assessing convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 609,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_heterogeneity(data, k, centroids, cluster_assignment):\n",
    "    \n",
    "    heterogeneity = 0.0\n",
    "    for i in xrange(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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 610,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.25"
      ]
     },
     "execution_count": 610,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_heterogeneity(data, 2, centroids, cluster_assignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combining into a single function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 611,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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 xrange(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",
    "          np.array_equal(prev_cluster_assignment, cluster_assignment):\n",
    "            break\n",
    "        \n",
    "        # Print number of new assignments \n",
    "        if prev_cluster_assignment is not None:\n",
    "            num_changed = np.sum(np.not_equal(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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting convergence metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the above function to plot the convergence metric across iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 612,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 647,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "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",
      "[19595 10427 29049]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAEKCAYAAAAo150IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXm8VVX5/98fwTlnSMkENMsE01TEIYeLQ1pGDn0dKTUN\nmxwqLLVCcSJzyq/m9Evrm0hWJDjklNNVEscAE5xSwAEcQBREBYP7/P541vHuu+85995z77n33OF5\nv177tc9e+9lrPWuffc6z11rPepbMjCAIgiAIOg8rVVuBIAiCIAgaEsY5CIIgCDoZYZyDIAiCoJMR\nxjkIgiAIOhlhnIMgCIKgkxHGOQiCIAg6GWGcg26DpAGS6iT9vtq6tARJp0h6VtKHklZI+nqF8h2T\n7sPulcivMyHpmFS3o6qtS08h3e/7q61HTyOMcyLzx35rEzJHJ5mfVKCcLmFAuiCWto+RNEfSrCrp\nUxRJ3wQuAN4DLgbOAp6rUPaN7kF70sH3t0Pr1hOQVCuprgmRLn3P0//tv0uc+356MZ4tafMKlbej\npHGSZkhaKOkDSc9Luk7S51qaT+9KKNPD6LIPaQ9gLrAlsCiX3hm/s/1wvfY3s/nVVqaNdOT9nQg8\nArzegWV2d5ozvlsCH3SQLh2GpNOAsfhL8d5mNq9CWe8K1ACPAvcAS4HPA98EjpS0n5k92FwmYZzL\nR1W+PiiBmS0HXqi2Hi2kH0A3MMwdipm9h/c2BB2EmXWV31SLkXQ+8DNgKrCfmS2oYPaXm9nFRcqs\nAe4HfgXs0lwm0a1dISRtJukPkl6VtCztr5DUJyNzNDALf0stjJ3VpW6V3TNyq0j6qaTpkt6X9K6k\ne4qNIaYuqRWSVpU0VtIsSR9JOikjs7WkmyS9JWmppBcknStpzSL59ZZ0RsrnQ0nPSPqepD2SrmcU\nuWZbSRMkvZHyfzHpsmZO7uM8JG2f6rQ41W+ipAEl7m1L828wZFA4BvoDAzP3u6BDTfp8eYlyd0nn\nry52voh8s/e5MDQCDEvHBX1a1C0saU9Jt2XKmJ260AY1c11hSKbRWG2p71bSkPS9vJLKelPSFEk/\nTOebvL+5vL6RntV35d18UyWNLKLLx+Plko5Nch9ImpjONxpzzn7vkj4jaZK8O3FJesa2LnFPDpM0\nLT3nr0m6WNJqKnOMVdJASX+UNE/+258j6TI1/O2vLuk9SU+XyKO3pAX550DSWukZKvgmLEj1a1Sn\nVO4sSetKulL+H7RcTfgypO9vd//Y4Pv7fVYmfz/U8H/n1+kZeV/SPyXtkGQ+LelGSfPTdzFJUr8S\negyTdEeq34eSZko6VVKvUrq3FklX4ob5n8CeFTbMmNlHJdJrgXeAz7Qkn2g5VwBJOwN3AasAtwKz\n8W6M7wFflrSDmb0LTAcuBX6UPt+cyWZOymtVvCvkS8ATwDXAmsABwH2SDjGz7HWF7qhJwOCkx2Lg\ntZTfHsAd+IvYX4B5uGH4edJtNzNblsnveuBwvKvnf4F18De9f1Kk60vSQcCNwLJUnzeA7YDTgBpJ\nu6cWbZahwKn4W+TVwLbAgcBWkrbKPtytzL/Au8AY4MdJ999Q33NRa2YPSfoP3tU0qsiP6rh03e9K\n5J+9Dy29z9OSTt/GjdqZSad3W1DGKOBC/PudlMrYBNgTf1aeaSaLFnc/S/oi/p0vAW7Bu5H7AFsD\nRwFX0Mz9zeR1ITAKf8b/jHeR7gNcI+nzZjYqp6Ph3++uwG34M/1e7nwxNsW7EmcA1+F/ggcC90va\nMttLIel4/NlbgH+/y4CDgBaPCaZ8Po/fp3Xx7+QFYAhwAvBVSTuZ2QIz+1DSJGCEpK3NLD8Guh+w\nPnBVJu8NgMnAFsADwN+BDYBvAPtI2svMHsvkYcCq+O9q9aQPwMImqjCG+mdxDPXf3/Rmql74Dv6C\n/9dNTPofAdwt6Uv49zYH+D9gG/w/bF3Si2mmnifg/zXz8d/4O/h3/ytgB+B/mtGlRSRD/0fgSOBu\n4CAzW1qJvFtY/s7Aevj30zxmFpsv/jEAqMON0pkltknACuAnmetWBl7Gf+Sfy+X5jZTnZUXK+X0J\nPX6VyvhpLn0D3Oi/CayaSX8g5fcEsFbumpWAl4DlwK65c39I5YzOpO2T8von0DuT/lng/SR/Rk6n\nRcB/gA1z+Y8qcq/2SPmvAP4nJ//HlH5oG/Ivem/TfZtV4n7/NOVzeC59TdwITmvBs1PWfc58byvK\neD63SfnPAjYqUn7fzPGZqczdM2lHp7SjiuRd+F6y3+0lSX6rIvLrlXF/90153wSsnEnvRf3vafuc\n7nW44f98kfwa1SPzva8ATsnJn53Sf5ZJWxc39m8DG2fSVweeSvL3t/B7qU3yI3LpZyWdrsukfTml\nXVAknxtTPp8vknZoTnazdH+eKvI9rMBfaFYp49lq8llMOt9f5Jo63NBk/49+ktIXAr/KXXNL0m/b\nTNog4CPgYeATOfnLk/zBLa1LE/o/hxv+OmACmf+3Er+1Ujag2HZAiXx2TufHpjI/xF9yt26R3m2p\ndHfacj/w5rasQTg4XXdKiXyfAN4sUk4j44y/tS4Eni6R1w9T+V/NpD2Q0r5SRH73VNbEIuf64Y4K\nL2bS/i/l9eUi8lfS2Dj/hCKGNlOXN4HHM2kFI9Dojy9z7sI25N8a49wXbzX9I5d+XMrrhBY8O2Xd\n5+z3VsbzeVW6F0e0QLYSxvniJL95C8pr6v7eir9U9C1ybjA5Y0W9cf51ifyaMs4vFpEvnJuQy6MO\nGFtE/tBSz2gR2f5JdmqRc6ul5/N9kiHAX6JeB17Nya6J91A8mUnbIN23v5co+8J0HwblvocVwJYt\nfa5a8iwWux/U/+/skkvfOMkvAlbLnRuRrjk6k3ZZShtSpNy10rm/llOfEvqvSPtngZWakS88Yy3d\nSjW0RuXkXiDzYtLcFt3ajbndzIqO0cjHjP+QSx6Kd/FsJenMIpetDvSRtL6ZNdW9BN59tS7wcom8\nPosbpc/jXahZ/lVEfpuk20P5E2b2eurSHSRpTTN7H++yBO8azDMF76bPMjTtd5U0OHdOwH+Trnmm\nFkl7Le3XrUD+LcbM5ku6GfiGpAFm9nI6dRxuVG9oQTbl3ufWMCTt72nl9eUyATgZeEzSjcB9wGQr\nf3xuKN5K/YHUyBdylbTPf4dG8ee5OYp1xRZ7rgrfV6nnvKVsk/bFvvelkh4Hvor/rmeaWZ2kPwMn\nSaoxH4MEf8Ffg4bP2g64MV+zxH/Blmn/eRoOZyw1s2fLqENbyXfPF7zo/2ONu4zfwH+3n8qkFf4/\nh0vaPycvvLXZpt944jX8BWhLfPjl5FKCZvZHvCevTZg7hV0saTW8h+AMYIqko8xsQnPXh3FuO+vj\nD9G3mpAx/O24OeO8ftpvTb2hLJVXw0Szt4rIrp32b5bI6w38oVkbf8NfC/ivmS0uIlss/0LdTyyR\nf0HfPMXyL4wbZx1AWpt/ufw/4BB87G2MpC2BnYAbzX0FmqPc+9wa1gGWt8I4tgoze1RSYcz8OOD7\nAJIewnuJWmo818e/00aOhIWicMOUp9S9bIpGz5WZrUgvBdnnaq20L+YpX+w5L0VLvvesHLgBPhmf\nVlOb0gotyj9n5Ar/BbunrRT5/4Jy9G8zZrYkd1yX7ndTv/GVM2mF3/gvmyim2PNRLouAvfF7fqKk\nOjP7cQXybZb0kjJV0jfwntTfSbrXzN5p6rowzm1nMf4Hs5+ZtbVVU3ig/2JmR7Yxr0J+AjYscX7D\njBx4C2dlSWsXMdCfLJG/4eNk/2mrslXIHwAzu0/SS3h31hjqHcGubWEW5d7n1vAu8BlJfVppoOtw\nHYv95tcpdoGZTQa+kt78d8adq74H3JkcuZp72QSv80dm9qlmJXPFlylfDgXnsr5FzhV7zktR+D5b\n/L2b2b8kPY/31PwAb9HvhXcbv5G5tnDNr83s52Xo1J73rT0o/MbXtIaOqRXHzN5KL5z3AydLopiB\nlrQN/qy3lOlmdksLyl8h6QG8x2U7vDeqJDGVqu08jv/p7dxC+RVpX2yKwLP4H8cQFekDbAWFbr5i\nU7A2wj1TZ2W6Wp9K+52K5LULjX/45da9XCqV/wqK3+8s1wL9U9faN4HZZvZAC/Mv9z63hifS/sut\nvL7wlr5xkXPbNXWhmS01swfM7GTc92ADfDZBgabu7+PAhpL6l6lve/IU/lyVes5bSuF73y1/Is26\nGIoPjTyfOz0eb00PBw7D/4fzwydP4L+3YjpWmhXgc6k6oKw8hd/4jh1RmJm9iXuLP4sb6EuKiH0R\n7+lp6XZAGSoUfn+lZph8TBjntnMz8CrwU0mNHjD5vMmhmaR38B/dJnlZM1uBT+/YHDhfUqPvR9LQ\n1JJpCf/EvXu/nqY2ZBmLdy9lx1ZuxH8ooyV93PUk6bP49Jk8f8DHcX6lImHpJK2dpuS0lkrlvxAf\n91+lCZk/4D+Ya/AWVTnhVcu9z63h/+HPzXmSGrRCJfWSVKwVmOVf6frDk+EoXPtZ4CQahzzdqcT9\n2ijts+OJTd3fy/Fn6jpJjVro8jnKA5rRvdLcik/nOl7Sxy8rktbAu/Fb1Po0s1eBB4EvSjoid/o0\nvBV+ozWe6jcevyffxLu0P8SnImXzfhP4G7BHmmrUCFUudnqhB6TRf1IHcCXeq3N5sTnQkj6Zpqtl\n0+bI51636oUvDQEWDPSP8gbazP5oZr3K2I7N6Vf0ZVfS3vh0vXeAx4rJZOnwbu00H7RYi+RdM1s/\nIzcIOBd/o1oHny/3B+DSZMRIX85l+JvOJ/HxvJl4V9CduXLn4N6VWQyf63ZrTnYk7im8aSr3N2Z2\nTbH6mNlHkg7BHbSmSPoH7qDRGxiIe8JOwR1DMLP3JT0B7C7penyaUB1wffqxnwFsD5yC/9lPxn88\nn07pn6PeA7hJzMwkHQvcCdwr6a94iMsa/I38CdzrsyD/jyRzCPCUpNvwe38Y3hW0f9K1ID9f0gh8\nrOxpSXek+qyZ7l0N7gH+g+Z0LaF/pfK/H793d6X7+RHwUOq2LZT1VqrvwXhL4v/K0LOs+9wazOwp\nST9L+TwrD8oxD38T3yulX9bE9a8nx64jgH9Jugv/zRyU9M7PJT0Vn0f+EO4F/BH+W9wNeJKGczVL\n3l8zu1PSr3Bj9R9Jd+POOX1x55wd8XmnL9Ny2tTCM7N3JJ2Cz9Wenpy0luL34hnc36OpWNNZvo/P\nRb5e0sHUz3PeB59ed1qR8mdLegT/PfXGh7GK9ap8H3cm+19J38b/0Jfg/2M74d9fJcZj78e//4mS\n7sTvxVNm9vcK5N0kZjZD0on4S9wLkm7H/3PXwxspu+Hj0dm488L/u5ttfTZRbqGL+wHcQGNmrV4z\nIcdNkj7CHV9fwZ2Cv4Dbgo+A44o4yxVVskO3pOAK/A91aGbbLiPTD3dsmIrPFa4BzknX/SojNwgP\nIDAi5TscfyuuAw7MlTsbN6BDc9s66fyAlP/0tD875VmYJ/ld6l3sf1KkXp/G/xz/g78Jv413n11K\nZh5nkt0cn4v4Nv6A5ae9rIT/MKfgY43vAy/ib9cjyEwFwB+u5c3c863xt/D5+A/vhXQ/1ygi2xt/\nQZid6vFMqnthytjJRa75PN7SfDnl/xb+B34embnfme9+dJE8Cvf/ujbkXzQP3JhfjRuFj8hNCcvI\nDU91vK2Vz3Y597nZ761EGXun5/htvPU3C3+R2DIj02gqVUpfFfdUnZeunYa/eDX6XnDj8ge8dbEI\nf9ufjhub/HzUZu8vPt/5tvTdLcV7mx7AA/Ks35zumfOlplIVfXbS+RXAfUXSD0v34IOk+0X477gO\nmFTGdzIwfQfzUt3m4EE1+jRxzfeTXsvx+Oql5FZP9/xf+JDXYrybfDy5+bX4b/alVjxTvfD4CrPx\naYUNpgcVu39NPb9N3O+mfv874gFN5qZ7OA////sF8OmM3Drpnj1YRv3qyM0Jz5z7JN6gWwFcXO69\nK5Hnd3A7NCc9Wx+k7+xqypjm1mZFWqF44QvaswmZ4ykyxxLvdp3bggftFeCWIg/u9S249k0az5O9\nDv9T6dXR96uzbNS/HO1XbV3asY6jUx2LBhWIrftveC9EHXB+tXWJrej389X0/exbbV3ae6vWmHNz\n3VKF8c786kKLaGac3LzLexGt6/LYGQ9ROD6XPg53gtm1FXl2KSQ18jyVtAU+nWkxmbCM3QlJq+Oe\nyPPwFl7QjZG0Xn6MPI2Jn4d3mTbrfRtUhS/hQV/urrYi7U01HcLGy4OyL5A0XlLWGWECHg7zCnlQ\n+bXkMZZH4F1PDZDTS9KG8oD7n8XHMPIMlwdnXyrpEUl5L7tCoIsZufSZ+AtFk4sLdBN+Kek5SddK\nOj+Nn07H54b+xDowFm1HIOlLkn6BB/fYCDjPzFo63hh0XfYE5kq6IT3nhS78HYAbzOyR6qoXFMPM\nfmFmQ5qX7PpUY57zItzAPoi3xLbFxxWmSNrWPEj8W5J2wd9eC6u01AFjrMhSXPii9YXg+e/hcZJr\nczK34o45s/H5hycAkyR908z+lGQKDmn5yeELc+e7M3fiY7xfw50y3scdXi4xs7uqqVg7sTc+xv4W\n7s/QohWogi7P0/g8093x3rI6fFzwV8Bvq6hXEACg1I9fXSWkbfH5bmPN7Ez5Umv34wPp5+PGcU/g\ndOCXZnZh7vpP4a2ejfApPwcA3zCzfIjL7DUr4eH7PmlmA1Pa6biH+OrWcGWkXnioyNFmdl5FKh0E\nQRAEJegUEcLMbJqkF6iPpXwqPl2gv9VHqnpIUm/gHEnXWSY6kZnNw8cKAe6QR2G5iMbxp7Nl1kma\ngM8n3tB8XmGhxbweDUPyFVrMJSMiSar+W04QBEHQoZhZuwRv6axBSLbCpwTkwx0+jjuLbd7M9U+2\nQKYYhbHl/CILhbHmJtfLrbZ3X7W2M888s+o6RP2j/lH3qH9Hb+1JpzDOkobgk+0Lq8S8gccRzkcU\nKoSym9tEXsInrr/UTJm9gMOBV8xbzQCP4I5oI3Li38LnlT7cdE2CIAiCoO1UI0LYONxwTsMdwrbD\nJ9m/Sr2H9dV41KB7JF2IG8ZhuNPXRDObm/I6E+9yfhg36BvhE8CH4JGQCmUejjs43YEb9n742shf\nxA00AGa2XNJo3Et8HnAvPu/xGHxd31ZHpAmCIAiCllKNMeeZuEE8GQ899wYeVWmMpXFkM3tM0m64\nF+2leJD4OfiKQdk4qFNTPofhkWPewKNy7Wpm2bVaZ+OG+2LcmL+Pd33va2b3ZpUzs2sk1eEvAqfg\nAU1+aCXCdwZQU1NTbRWqStS/ptoqVI2eXHeI+rcnncJbuzsgyeJeBkEQ9BxSTO4e5RAWBEEQBD2W\nMM5BEARB0MkI4xwEQRAEnYwwzkEQBEHQyQjjHARBEASdjDDOQRAEQdDJCONcQQ4/vHmZIAiCIGiO\nMM4V5JVXqq1BEARB0B0I41xB5s+vtgZBEARBdyCMcwUJ4xwEQRBUgjDOFWTRIli2rNpaBEEQBF2d\nMM4VZsGCamsQBEEQdHXCOFeY6NoOgiAI2koY5woTxjkIgiBoKx1unCXtIamuyLYwJzdI0kRJcyUt\nkTRD0ihJvTIy/SXdLGmOpA8kzZdUK+krRcqVpNMlzZb0oaTpkg4uoeNISc9KWirpOUnfbWn9wjgH\nQRAEbaV3lco14ETgyUza8sIHSf2AWuA14CTgbWAv4AKgD3B6Ev0EMB/4RZJdGxgJ3C7pYDO7OZP/\nucBPgJ8DU4HDgQmS9jezuzJljwSuBs4D7kvlXpnW7bymuYq99VbLbkAQBEEQlEJm1rEFSnsA9wP7\nmNn9JWSOB64CtjCzFzPpNwK7m9nGTeTfC5gNTDOzA1JaX+BVYKyZnZ2RvRfoY2ZfzFw7D7jdzI7N\nyF0HDAf6mdmKEuUaGD//OZx3XkvuRBAEQdCVSY02tUfe1Rpzbq4yK6f9olz6IprRORnPRWRa4sB+\nKc/xOfEbgC9IGpCOd8Zb5nm5ccAGwK7N6B3d2kEQBEGbqaZD2HhJyyUtkDRe0iaZcxOABcAVkgZK\nWkvSQcAI4KJ8Rmk8uZekDSWdAXwWuDwjMghYZmYv5S6dib8oDErHg9N+RjNyJQnjHARBELSVaow5\nL8IN7IPAYmBbfMx4iqRtzWyBmb0laRfgFmBWuq4OGGNmFxfJ8wJgVPr8HnC4mdVmzq8PvFvkuoWZ\n89n9O83IlSSMcxAEQdBWOtw4m9l0YHomabKkycDjuJPYmZL6AJOAJcDBuHHcExgtaZmZXZjL9jfA\njcBGwFHAjZK+YWZ3tG9tGhMOYUEQBEFbqZa3dgPMbJqkF4ChKelUoD/Q38wWp7SHJPUGzpF0nZkt\nzFw/D3fkArhD0gN467xgnN8B1i1SdKElvDAjB7Ae8GYTciWJlnMQBEHQVjqFcS7CVsBLGcNc4HHc\nsWvz9LkUTwInZ45nAqtK2szMZmXSB+PTup7JyCmlZ41zYaz5GZpkDO++C6NHw1571VBTU9O0eBAE\nQdBlqK2tpba2tkPK6vCpVEWVkIYAjwLnmNlZkv4AHAQMMLNFGbmxpFa1mc0tkZeAR4B1zGzLlNYX\nnwd9rpmdk5G9F+hrZtuk4954C/w2MzsuI3ctcAA+lSrrBZ4t19zOw7x50K9f6+5FEARB0DVoz6lU\nHd5yljQOeAmYhjuEbQechs9DLnhYXw0cCdwj6UI8CMkw3OlrYsEwSzoT73J+GHgDH3P+DjAEOKJQ\nppnNl3QJcLqkJdQHIanB5y8X5JZLGo17ic8D7sWDkBwDnFDKMOeZPz+McxAEQdB6qtGtPRM3jCcD\na+BG9W+4J/ZCADN7TNJuwBnApXjkrznAGOCSTF5TUz6HAeukvJ4CdjWzR3Pl/hz35D4JN+LPA4eY\n2Z1ZITO7RlId/iJwCvAK8MOWRAcrEE5hQRAEQVvoFN3a3YFst/af/gRHHNHMBUEQBEGXpjtGCOvW\nhMd2EARB0BbCOLcDYZyDIAiCthDGuR0I4xwEQRC0hTDO7UA4hAVBEARtIYxzOxAt5yAIgqAthHFu\nB8I4B0EQBG0hjHM7EMY5CIIgaAsxz7lCZOc5A/z3v9C7s0YuD4IgCNpMzHPugixYUG0NgiAIgq5K\nWcZZ0vGS1mwvZboT0bUdBEEQtJZyW85XAfMkXSFp6/ZQqLsQxjkIgiBoLeUa588AVwIHA9MkPSLp\naEmrVV61rk0Y5yAIgqC1lGWczWyOmZ0ObIKvLPUB8HtgrqTfSNqyHXTskoRxDoIgCFpLqxzCzGy5\nmU0ws72ALYCn8aUYZ0h6UNL+lVSyKxJRwoIgCILW0mpvbUlrSfoBcBOwOzAd+AW+RvStks6ujIpd\nk2g5B0EQBK2lbOMsaYik3wHzgItxo7yzmW1vZueb2ZeAMcAPS1y/h6S6ItvCnNwgSRMlzZW0RNIM\nSaMk9crpcq2kFyS9L+llSTdIGlik3DlFylwh6etFZEdKelbSUknPSfpuufcpjHMQBEHQWsoKkyFp\nKrANMBs4G7jOzBYWEb0HOKuJrAw4EXgyk7Y8U04/oBZ4De8ufxvYC7gA6AOcnkQPAwYBlwIzgE8B\nZwBPStrGzObmyrwLf3HI8nyujiOBq4HzgPtSuVemyebXNFGnBoRxDoIgCFpLuTGsXsO7ru+ypkOL\nTQU2bSav58zs8RLnhgMbALuY2YsprVbS5sBR1BvnX5tZg3AfkqbgLw8jaWyIFzRRJqlVfi7wRzM7\nIyU/KGlj4BxJ15rZimbqBYRxDoIgCFpPud3aFwGTixlmSZ+QtDuAmX1kZi83kU9z4c5WTvtFufRF\nZHTOG+aU9gowH9i4mTKKsTPeMh+fSx+Hvyzs2tKMwiEsCIIgaC3lGucH8G7kYmyRzreU8ZKWS1og\nabykTTLnJgALgCskDUzOZwcBI/AXhJKk6VyfBJ4pcnp4GptemuZoH5A7PzjtZ+TSZ+IvFKXq3oiF\nC2FFi9rYQRAEQdCQco1zUy3eVYGWmKNFuIH9DjAMH7veG5giqQ+Amb0F7IIbw1npmgl4N/bFJZXz\nbumrgbfw+ddZbsXHub8MHAl8CEySdGRGZv20fyd37cLc+aKst179ZzN4++2mpIMgCIKgOM2OOSfP\n580ySUMkfSIntjpwLPBKc/mZ2XTcw7vAZEmTgcdx43lmMtKTgCV4NLKFwJ7AaEnLzOzCEtlfAewE\nfNXMGnSJm9nJuXrdDDwKjAX+1JzeLaFvX3gnY9bnz4dPfrISOQdBEAQ9iZY4hB0NnIl7OxtwOQ1b\n0JaOl1Ni+lRzmNk0SS8AQ1PSqUB/oL+ZLU5pD0nqjTtmNfISl3Q+3ho/yszua0GZdZImAOdL2tDM\n3qS+xbwe8GZGvNBiLuaZ/jEffjgmc1TD/Pk1zakRBEEQdBFqa2upra3tkLJaYpz/D5/WJOB+3ADn\nx3OXAS+UmFbVGrYCXsoY5gKP485im6fPAEj6BfBT4AQza0sruDC2PJiGxrkw1lxsHPtjtt9+DK++\nWn8cTmFBEATdh5qaGmpqaj4+PuuspmYMt41mjXPyun4ZQNIwYKqZvVdJJSQNwR3K/pqS3gB2lrRO\nrnt6p7Sfm7n2JOAc4HQzu6qMMnvh8cFfSa1mgEdwR7QR+ItIgW/hc60fbirPfBd2TKcKgiAIWkNZ\n85zN7MG2FihpHPASMA1YDGwHnAa8ineZgzt1HQncI+lC3DAOA0YBEwvBRSQdDvwGuBOfB71jpqjF\nZvZsRu5rwB24Ye+H9wB8ETfQhfotlzQa9xKfB9yLByE5Bm+VL6cJ+vZteBzGOQiCIGgNLXEImwUc\nZGZPSZqNjzGXwszsM81kORM3iCcDa+Ct5L8BYwrd4mb2mKTd8GhflwJrA3PwoCKXZPLaN+33S1uW\nB3EnMvCgJBvh4UbXB97Ho5Pta2b35ipwjaQ6/EXgFNzJ7YctiQ4WxjkIgiCoBC1pOT+It3ALn5sy\nzs1iZucD57dA7nG8tduUzLeBb7cgr8fw6Vot1fF3wO9aKl8gjHMQBEFQCVoy5vztzOdj2lWbLk7e\nOIdDWBAEQdAaWr1kZNCYcAgLgiAIKkFrlozcNi3luCCF39wupY+VlB/37VFEt3YQBEFQCcoyzpJ2\nxacbfR7tB8OSAAAgAElEQVSPqpW9vg74XuVU63r06dPw+O23oa6uOroEQRAEXZdyW87nA3fjQTp+\nkjs3FZ8W1WNZZRVYZ53647o6XwAjCIIgCMqhXOO8HXBVWjIy77W9AOjb+JKeRTiFBUEQBG2lXOO8\nFJ+bXIx+NF5/uccRTmFBEARBWynXOP8T+FEKfVmg0II+joYhL3sk4RQWBEEQtJWywncCo/H40k/h\nUb0MOFrSJcD2wA6VVa/rEcY5CIIgaCtltZzN7Clgd3zFpl/gKzidkE7vYWbPV1a9rkcY5yAIgqCt\nlNtyxsymAntJWg2PU/2umX1Qcc26KOEQFgRBELSVso1zATNbCsyroC7dgnAIC4IgCNpK2cZZ0mbA\noUB/YLXcaTOz4yqhWFclurWDIAiCtlKWcZZ0IPBXfKz6LWBZTqRNK1Z1B8I4B0EQBG2l3KlU5wC1\nQD8z+5SZbZrbNmsuA0l7SKorsi3MyQ1KMbznSloiaYakUdlpXJKGSLpW0guS3pf0sqQbJA0sUq4k\nnS5ptqQPJU2XdHAJHUdKelbSUknPSfpuS29QGOcgCIKgrZTbrb0ZMMrM2mpyDDgReDKTtrzwQVI/\n/CXgNeAk4G1gL+ACoA9wehI9DBgEXArMAD4FnAE8KWkbM5ubyf9cPOToz/FQo4cDEyTtb2Z3Zcoe\nCVwNnAfcl8q9UhJmdk1zFcsb5wULPIznSrH+VxAEQdBC5JE4Wygs/Qu4yMxubHWB0h54sJJ9zKxo\n0BJJxwNXAVuY2YuZ9BuB3c1s43Tcx8wW5K7tD8wGzjGzMSmtL/AqMNbMzs7I3gv0MbMvpuNeuJPb\n7WZ2bEbuOmA43mOwooTOVriX66wDixfXn1uwADbYoPl7EwRBEHQdUqNN7ZF3ue25nwE/T05hbaG5\nyqyc9vlwoIvI6Jw3zCntFWA+sHEmeb+U5/ic+A3AFyQNSMc74y3zvNw4YANg12b0BqJrOwiCIGgb\n5RrnMbiRejaNAT+U2x4sI6/xaT3oBZLGS9okc24CvpDGFZIGSlpL0kHACOCipjKVtCXwSeCZTPIg\nYJmZvZQTn4m/KAxKx4PTfkYzck0SxjkIgiBoC+WOOa8A2hoFbBFuYB8EFgPb4tHGpkja1swWmNlb\nknYBbgFmpevqgDFmdnGpjFO39NW4J/nvM6fWB94tcsnCzPns/p1m5JokjHMQBEHQFsoyzmZW09YC\nzWw6MD2TNFnSZOBx3EnsTEl9gEnAEuBg3DjuCYyWtMzMLiyR/RXATsBXzaxqK2RFlLAgCIKgLbQ6\nQlglMbNpkl4AhqakU/EgJ/3NrOBa9ZCk3sA5kq4zs/zUq/OB7wBHmdl9uSLeAdYtUnShJbwwIwew\nHh4/vJRck0SUsCAIgqAttCZC2MbAKHwBjA2A4WY2Q9KPgEfM7LEK6LUV8FLGMBd4HHfs2jx9Luj0\nC+CnwAlm9qci+c0EVpW0mZnNyqQPxqd1PZORU0rPGufCWHN2HLsRY8aMAWDqVICatIVxDoIg6A7U\n1tZSW1vbIWWVO5VqMDAZH3t+BNgf2MHMpkr6DbChmR1ZthLSEOBRfPrTWZL+ABwEDMh2T0saS2pV\nF+YwSzoJn+d8upn9ukT+ffE50+ea2TmZ9HuBvma2TTrujU+lui0bhlTStcAB+FSq5RQhO5Vq3Dg4\n6qj6c4cfDje2evJZEARB0Blpz6lU5bacLwaeBfYFlgIfZc5NAYoaxyySxgEvAdNwh7DtgNPweciX\nJ7GrgSOBeyRdiAchGYa32CdmDPPhwG+AO4FaSTtmilpsZs8CmNn8tOb06ZKWUB+EpAafv0ySWy5p\nNO4lPg+4Fw9CcgzeKi9qmPOEQ1gQBEHQFso1zrsCR5jZkmwYzcSbwEYtyGMmbhhPBtYA3gD+hnti\nLwQws8ck7YZH+7oUWBuYg0/luiST175pv1/asjyIO5EV+DnwHh5xbCPc6/wQM7sze5GZXSOpDn8R\nOAV4BfhhS6KDFQiHsCAIgqAtlNutvRgYYWa3JeP8X2BI6tY+GPidmfXIWFjZbu1XX4X+/evPbbQR\nvP56lRQLgiAI2oXOFCHsceDbJc4dCjzcNnW6B8Xia5fxDhQEQRD0cFqzKtVwSf8AvoV7Ou8t6Y+4\nA9d5FdavS7LaavCJT9QfL18O7xYLgRIEQRAERSjLOJvZg8CBwKZ4BC4B5wO7AQdWaBpVtyDfen7z\nzeJyQRAEQZCn7IUMzex2M/ss8DncQWxLM9ss71jV08mOOQM88UR19AiCIAi6Hq1eZdjMXjSzKWbW\n1ljb3ZLddmt4/MAD1dEjCIIg6HqU6619VBOn6/BFLaaZ2WttVayrkfXWBrjvPth77/rzm24Ks2YV\nuTAIgiDokrSnt3a5xrkOdwKDhmsyZ9PqgL8A3zazbJCSbk3eOH/wAay3HnyUuQNz5sCAAY2vDYIg\nCLoenWkq1ZeAl4HfAnsAn0/7K/FgHfvj0b4OwgOG9FjWWAN23LFh2oPlrHYdBEEQ9FjKNc6nAH82\ns5PNbLKZvZD2JwI3Aseb2UV4mM/DK61sV6OmpuFxjDsHQRAELaFc4/xlIL8cY4H78TjUAA8BG7dW\nqe5C3jh30GImQRAEQRenXOO8DNi+xLntqV8IYyXg/dYq1V3YeWdYZZX64zlzfAuCIAiCpijXOE8A\nzpI0StIASaun/Sn4GPNfktwX8YUlejSrrw477dQwLcadgyAIguYo1zj/BLgJuACYBSxJ+1/jK0uN\nSnIz8HWXezwx7hwEQRCUS1lTqT6+SPocsBO+9OLrwOM9PRhJfipVgdpaGDas/njAgOjaDoIg6A50\nmnnOQWlKGeelS2HddWHZsvq02bNh4MCO0y0IgiCoPJ1pnjOS1pB0gqQJku5L+x9IWr2F1+8hqa7I\ntjAnN0jSRElzJS2RNCONdffKyY2VdLekBSmfolHMJM0pUuYKSV8vIjtS0rOSlkp6TtJ3y7lHWVZb\nrfG4c3htB0EQBE1RlnGWtBEwFbgMGAKskfa/BaZK2rCFWRlwAt41Xtg+DnYpqR9QCwwETgK+BkzC\nx7rPzeV1ArAacBv1kcpKlXlXrsydgQYuWpJGAlfjzm/7An8FrmyLgY5x5yAIgqAcyg3feT1usA42\ns4cz6bvgjmJ3m9kxzeSxBz4neh8zu7+EzPHAVcAWZvZiJv1GYHczazSHWtJngP8Ax5jZ9UXOzwYm\nm1nJ+OCpVT4PuN3Mjs2kXwcMB/qZ2YoS1xbt1gb30M4a6P79fdxZ7dIZEgRBEHQEnalb+yvA6VnD\nDGBmU4Bf4uE7W0JzlVk57Rfl0hfRhpW0WsDOQB9gfC59HLABvkRm2ey4I6y6av3xK6+EU1gQBEFQ\nmnIN3SfwlmUxXkvnW8p4ScvTWPF4SZtkzk0AFgBXSBooaS1JBwEjgIvK1DnLcEnvp7HkRyQdkDs/\nOO1n5NJn4i8Ug1pT6GqreUCSLDHuHARBEJSiXOP8PPCtEue+CTzXgjwW4Qb2O8Aw4Gx8vHmKpD4A\nZvYWsAtuDGelayYAvzazi8vUucCtwIl4CNIjgQ+BSZKOzMisn/bv5K5dmDtfNjHuHARBELSU3mXK\nXwRcnxy//oTPcd4IX+Rib0ob7o8xs+nA9EzSZEmTgcdx43lmMtKT8CAnB+PGcU9gtKRlZnZhmXpj\nZidnjyXdDDwKjE11aVeGDYMxY+qPa2vBLMadgyAIgsaUZZzN7AZJa+Ct3Wszp94EvmdmrTJyZjZN\n0gvA0JR0KtAf6G9mi1PaQ5J6A+dIus7MFhbLq4wy6yRNAM6XtKGZvUl9i3k9vE4FCi3mJssck7G+\nNTU11GSay0OHevf20qV+/OqrPt95s83aUosgCIKgo6itraW2g8Yky205Y2b/T9K1wBa40VoIPG9m\ndRXUayvgpYxhLvA47iy2efpcaQpjy4NpaJwLY83PNHVx1jjnKYw7Z7uzH3ggjHMQBEFXId/oOuus\ns9qtrBaPOUtaRdJUSV82szoze9bMHk77NhlmSUNwY/9oSnoD+IykdXKihXAec9tSXiqzF94d/0pq\nNQM8gjuijciJfwt4G3iYNhBLSAZBEAQtocUtZzP7SNKmwPK2FChpHPASMA1YDGwHnAa8ClyexK7G\nnbbukXQhbhiH4QtrTDSzuZn8dgf6Av1S0g6S3k8635RkDscDmdyBG/Z+wA/x1bMOz9RxuaTRuJf4\nPOBefI3qY4ATzKxNdR82DM48s/44xp2DIAiCYpQbhOSvwCwzO63VBUqn4QZxAB5h7A3caI7JtGCR\nNBQ4A9gWWBuYgztuXWJmyzJyDwC7FyvLzHolmR2B8/Du6vXxtaafBC4ws3uL6DgSfxEYALySyrym\nmXqVDEJSYNkyj7NdGHcGePFF+MxnmrwsCIIg6IR0moUvJO0G3IBPa7oZ99ZukIGZzaqkgl2Flhhn\ngL32gvszcdGuvhq+2+rAoEEQBEG16EzGOTu2XPTCQmu1p9FS43zuuTB6dP3xoEHw9NOwUnvGPQuC\nIAgqTmcyzkc3J2Nmf2yTRl2UlhrnF1+ELbaAusxrzl/+Aoce2o7KBUEQBBWn0xjnoDQtNc4A3/oW\n3HBD/fHgwfDvf0frOQiCoCvRmRa+AEDSSpK2Smszr1lppbo7v/xlQ0M8cyZMnFg9fYIgCILORdnG\nWdIPcQ/rf+NLP26R0m+WdFJl1euebLEFHH54w7Szz27Y1R0EQRD0XMoyzmmK0f/intqH0nDpx8nA\nNyqnWvfml79sOL/56adh0qTq6RMEQRB0HsptOf8EuNjMjscXpsjyHKkVHTTPllvCYYc1TIvWcxAE\nQQDlG+dNgbtLnHsfWLdt6vQsRo9u2Hr+97/hlluqp08QBEHQOSjXOC8ABpY4twUViHndkxg0CA45\npGHa2Wd7SM8gCIKg51Kucf47cIak7FpKltZf/jE+Fh2UQTYgCcD06XDrrdXRJQiCIOgclBuEpA++\nMtMmwGN4TOspwOeBt4BdzGxRO+jZ6SlnnnOeQw6Bv/2t/njbbeFf/4oFMYIgCDoznWaes5ktAIYA\nv8LXVX4JX9nqt8DOPdUwt5Uzzmh4PG0a/P3v1dElCIIgqD4RIaxCtKXlDPA//wM33VR/3KcP3HEH\n7LBDBZQLgiAIKk6naTlLmiVpmxLntpLUI1ekqgT51vOCBVBTA3fdVRV1giAIgipSrkPYQGDVEudW\nw9c/bpIU8rOuyLYwJzdI0kRJcyUtkTRD0ihJvXJyYyXdLWlByueoEuVK0umSZkv6UNJ0SQeXkB0p\n6VlJSyU9J6ndF3Xcemv48Y8bpn3wAQwfDuPGtXfpQRAEQWeiNbG1S/XdDgHeLSOPE4CdMtvehZOS\n+gG1+MvAScDX8KAnFwDn5vI6AX8xuK0J3UjXnQFcBuwHPAJMkLRfVihFQbsaX7N6X+CvwJUdYaAv\nvhhOP71h2vLlcNRRcMEFMcUqCIKgp9DsmLOkH+PTpAA2BuYDH+XEVgfWB/5sZiOayW8PPCb3PmZ2\nfwmZ44GrgC3M7MVM+o3A7ma2cZFrPgP8BzjGzK7PnesLvAqMNbOzM+n3An3M7IvpuBcwD7jdzI7N\nyF0HDAf6mdmKEjq3acw5y+WXw8knNzbGP/qRG/BYvSoIgqD6VHvMeRZwX9oEPJk5Lmw34QZ8ZAvL\nba4yK6d93vt7Ea1r7e+X8hyfS78B+IKkQnf8zkCfInLjgA2AXVtRdtmceCL8+c+wyioN0y+9FA4+\nGN56qyO0CIIgCKpF7+YEzOwW4BbwtwTgbDObXYGyx6cW7bt4SNDTzOzVdG4C3gV9haSfAW/j3d4j\ngDGtKGsQsMzMXsqlz8RfFAYBLwODU/qMJuQebEX5ZXPoodC3Lxx4ICxeXJ9+yy3w8MPeuj7ssJgL\nHQRB0B0pd57ztwuGWdInJA2QtHJz1+VYBFwEfAcYBpyNG94pKcgJZvYWsAtuDGelayYAvzazi8ss\nD7zLvdh4+MLM+ez+nWbkOoRhw+Chh6Bfv4bpCxbAEUd4K/qNNzpSoyAIgqAjaM16zl+TNBU3mLOA\nL6T0ayUd2dz1ZjbdzH5mZreb2WQzKzhobQScmPLqgzuALQEOBmpwh67Rkn5ars5dmW22gSlTYPvt\nG5+7+WaPz3399eEsFgRB0J0od57zgXgX9wLgVBqOHc8Gjm6NEmY2DXgBGJqSTgX6A182s5vN7CEz\nGwNcCJwjqdwW7DsUXzGrkM/CjBzAes3IdSgDB8Kjj8LYsY3Hod95B44+2qdcvZNv7wdBEARdkmbH\nnHOcCfzBzL4jqTc+tanADOAHFdJrK+AlM1ucS38cd+zaPH1uKTOBVSVtZmbZQCmD8elXz2TklNLf\nzMgNSvtnaIIxY8Z8/LmmpoaampoyVGya3r19mtUBB8Cxx8JjjzU8f/vtsNtucOedsMkmFSs2CIIg\nSNTW1lJbW9shZZW78MVSYLiZ3ZOmHf0XGGJmUyXtDvzDzFYrWwlpCPAocI6ZnSXpD8BBwIBsvG5J\nY0mtajObm8ujualUrwHnmtk5mfR7gb5mtk067o1PpbrNzI7LyF0LHIBPpVpeog4Vm0rVHCtWuOf2\nL38JS5c2PPfpT3tUscGDi18bBEEQVIb2nEpVbst5MT7VqBgD8TnQTSJpHL5gxrSU33bAafg85MuT\n2NXAkcA9ki7EvbWHAaOAiVnDnF4K+gIFt6kdJL0PYGY3pf18SZcAp0taAkwFDsfHsocX8jKz5ZJG\n417i84B7gb2AY4ATShnmjqZXLxg1yruyjzqqYSv6tddg11192cnddquejkEQBEHrKbflPB53ANsd\neA9vOW+Pd/dOBqab2fHN5HEabhgHAGsAbwB3AGPM7M2M3FB8OtW2wNrAHOBPwCVmtiwj90DSpxFm\n1isjJ+B0fC72RsDzwFlmNqmIjiPxF4EBwCupzGuaqVeHtZyzfPghjBgBk3K1WHVVuPFGOOigDlcp\nCIKgR9CeLedyjfNAfKzXcIN6FPA3YGtgHbyLe17FtewCVMs4g3dzn3giXHVVw/SVVoIrroDvfa8q\nagVBEHRrqh0h7GPMbA7eDf13YB9gBd5qfRTYsaca5mrTq5cb4XNzUcfr6uD734dLLqmOXkEQBEHr\niPWcK0Q1W85Zfv97OP54b00XWHllePpp2GKL6ukVBEHQ3ahqt7akM5oUaIhlvaF7Ep3FOINPqzrk\nEB+PLrDffnDHHRHuMwiCoFJU2zjXFUk2ii9eYVknrJ5EZzLO4AtnHHFEw7Rbb3UP7yAIgqDtVNs4\n541tb+BDYEd8SlIDSi2p2N3pbMbZzGNzP5hZpmOzzWDmTFit7JnoQRAEQZ6qOoSZ2YrsBhTm+q7I\nn+uphrkzIsFllzVc+3nWLF8POgiCIOjctGZt5KCLsPXW8INcQNWxY+HVV4vLB0EQBJ2DMM7dnLPO\ngg02qD/+4AM45ZTq6RMEQRA0Txjnbs7663trOctf/wodFLs9CIIgaAUtcQjbLJfUCw99eQC+ilMD\ncqs+9Rg6m0NYlhUrYOhQmJpx3/vCF/y4d7nR1YMgCAKg+t7adfjUqQbJRdKAhvGsexKd2TgDTJkC\nX/pSw7TLLvOwn0EQBEH5VNs4H11Ohmb2xzZp1EXp7MYZfAWrcePqj9ddF158seGYdBAEQdAyOs3C\nF0FpuoJxfv11D+H53nv1aRdcAD/9afV0CoIg6Kp0moUvgq5Nv35w2mkN06691gOWBEEQBJ2HaDlX\niK7QcgaYPx823hj++9/6tNpa2GOPqqkUBEHQJelWLWdJe0iqK7ItzMkNkjRR0lxJSyTNkDQqH05U\n0qqSLpQ0T9IHkqZI2q1IuXOKlLlC0teLyI6U9KykpZKek/Tdyt+J6tC3Lxx4YMO0a6+tji5BEARB\ncarVrW3ACcBOmW3vwklJ/YBaYCBwEvA1YBJwAZBbtZjfA8cBvwT2B14H7pa0dZEy78qVuTPwYFZI\n0kjgamACsC/wV+DK7mSgR45sePy3v8E771RHlyAIgqAxHd6tLWkP4H5gHzO7v4TM8cBVwBZm9mIm\n/UZgdzPbOB1vA0wDjjGz61NaL3z+9XNmdmDm2tnAZDM7qgndegHzgNvN7NhM+nXAcKBfqfjhXaVb\nG6CuDjbfHGbPrk+LaVVBEATl0a26tRPNVWbltF+US19EQ52/DnyEt26Bj1fF+jOwr6SVKY+dgT7A\n+Fz6OGADYNcy8+uUrLQSHHdcw7Tf/S4cw4IgCDoL1fTWHi9puaQFksZL2iRzbgKwALhC0kBJa0k6\nCBgBXJSRGwTMNrOlubxnAqsAm+fSh0t6P40lPyLpgNz5wWk/o0h+SuV1C445puGKVU8/DY8/XjV1\ngiAIggzVMM6LcAP7HWAYcDY+3jxFUh8AM3sL2AU3hrPSNROAX5tZdtHD9YFio6ULM+cL3AqcCHwZ\nOBJfk3qSpCNz+VEkz2L5dWk23hj2379hWjiGBUEQdA463Dib2XQz+5mZ3W5mk83sMmA/YCPceJKM\n9CRgCXAwUIM7go2W9LNWlnuymd1gZg+b2UT8heBJYGwzl3Zb8o5hN97YMEBJEARBUB06xbIHZjZN\n0gvA0JR0KtAf6G9mi1PaQ5J6A2dLutbMFuIt3P5Fsiy0cBcWOVcos07SBOB8SRua2ZvUt5jXA94s\nJz+AMWPGfPy5pqaGmpqapsSrzle+Ap/6FMyb58fvvw9//nNjox0EQRBAbW0ttR20pF+nMM5F2Ap4\nKWOYCzyOO4ttnj7PBA6UtFpu3Hkw7ij2IuVRGFseTEPjXBhrfqapi7PGuSvQuzcceyycm5mc9rvf\nhXEOgiAoRr7RddZZZ7VbWZ0ifKekIcAWwKMp6Q3gM5LWyYnulPZz0/423PHrkExevYBDgbvN7L+U\nIMkdDrySWs0Aj+COaCNy4t8C3gYeLqNaXYJjj214/MQT8NRT1dElCIIgcDq85SxpHPASPj95MbAd\ncBrwKnB5Ersad9q6R9KFuGEcBowCJprZXPDxa0l/AS6VtAowG/gBHrzkiEyZh+OBTO7ADXs/4IfA\nF3EDTcpvuaTRuJf4POBeYC/gGOAEM1te4dtRdTbdFPbZB+65pz7t2mvh8stLXxMEQRC0L9UIQnIa\nbhAHAGvgreQ7gDGZFiyShgJnANsCawNzgD8Bl5jZsozcqsB5uDFfF3gK+JmZTc7I7JhkBuPjx+/j\nzmAXmNm9RXQcib8IDABeSWVe00y9ukwQkjwTJsChh9Yfr7uuj0Ovvnr1dAqCIOjsxJKRXYCubJw/\n+gg+/WlfFKPAN77hK1gNGVI9vYIgCDoz3TFCWNCJWGUVOProhmk33QQ77AA77QQ33ADLlhW/NgiC\nIKg80XKuEF255Qzwn//AoEGwvMSoet++MGIEbLQRrLmmb2us4fvVV4eVV3bv78LWq1f9XvJoZCut\nVP9ZarxByz+XSmvN5zylzlXymqbyqiSVLKc1eZV7TUfd4466/0H3Jrq1uwBd3TgDTJwIP/0pzJpV\nbU2CoOtT6RfA1pTT3vSkl5wpUxoP84Vx7gJ0B+MMsGIF3HUX/Pa3vg+CIAh87YEddmiYFmPOQYfR\nq5fH3L7zTnj+eTj5ZFh77WprFQRB0LOIlnOF6C4t52IsWQI33wwvvggffOBhPrP7Dz7wsepiW12d\nb2aNP5vVb9Dyz6XSWvM5T6lzlbymox6TSpbTmrzKvaYr3uOg59DRLecwzhWiOxvnIOjJVPrFpJIv\nFNX8y+lpf3e9ezdcZhfa1zh31tjaQRAEnYKu5qQVdA9izDkIgiAIOhlhnIMgCIKgkxHGOQiCIAg6\nGWGcgyAIgqCTEcY5CIIgCDoZYZyDIAiCoJPR4cZZ0h6S6opsC3NygyRNlDRX0hJJMySNktQrJ7eq\npAslzZP0gaQpknYrUq4knS5ptqQPJU2XdHAJHUdKelbSUknPSfpuZe9CEARBEJSmWvOcDTgReDKT\n9vF6SJL6AbXAa8BJwNvAXsAFQB/g9Mx1vwe+ApwCzAZOAO6WtJOZ/Tsjdy7wE+DnwFTgcGCCpP3N\n7OMo0pJGAlcD5wH3pXKvTJPNr2lzzYMgCIKgGTo8QpikPYD7gX3M7P4SMscDVwFbmNmLmfQbgd3N\nbON0vA0wDTjGzK5Pab2AmcBzZnZgSusLvAqMNbOzM/ndC/Qxsy9mrp0H3G5mx2bkrgOGA/3MbEUJ\nnSNCWBAEQQ+iOy580VxlVk77Rbn0RTTU+evAR8BfCwnJeP4Z2FdSIZ/9Up7jc/ndAHxB0oB0vDPe\nMs/LjQM2AHZtRu8gCIIgaDPVdAgbL2m5pAWSxkvaJHNuArAAuELSQElrSToIGAFclJEbBMw2s6W5\nvGcCqwCbZ+SWmdlLReSUzgMMTvsZzcgFGWpra6utQlWJ+tdWW4Wq0ZPrDlH/9qQaxnkRbmC/AwwD\nzgb2BqZI6gNgZm8Bu+DGcFa6ZgLwazO7OJPX+sA7RcpYmDlf2L/bQjmK5JmXCzL09B9o1L+22ipU\njZ5cd4j6tycd7hBmZtOB6ZmkyZImA4/jTmJnJiM9CVgCHIwbxz2B0ZI+MrMLOljtIAiCIOgwOsWq\nVGY2TdILwNCUdCrQH+hvZotT2kOSegNnS7rWzBbiLdz+RbIstHALLd53gHVbKAewHvBmE3JBEARB\n0H6YWafY8HHdO9PnO4F/FZH5OrACGJqORwNLgdVycmOAD4GV0/G30nWb5eSOSekD0vFuQB2wZ05u\nj5S+RxP6W2yxxRZbbD1ray+b2ClazpKGAFtQ73X9BrCzpHXMLOuxvVPaz03724CzgENwj+rCdKhD\ngbvN7L9J7i58HvUI4JxMft8EZpjZy+n4EdwRbQQ+3avAt/C51g+XqkN7udMHQRAEPY8ON86SxgEv\n4fOTFwPbAafh85AvT2JXA0cC90i6EDeMw4BRwEQzmwtgZtMl/QW4VNIqeBCSHwADgSMKZZrZfEmX\nAKdLWkJ9EJIafP5yQW65pNG4l/g84F48CMkxwAlm9nGglCAIgiBoL6oRhOQ03DAOANbAW8l3AGPM\n7GpLadQAAAsJSURBVM2M3FDgDGBbYG1gDvAn4BIzW5aRWxWP5nUkPq78FPAzM5ucK1d4ZLGRwEbA\n88BZZjapiI4j8ReBAcArqcyIDhYEQRB0DNUea+7KG/Bp4G/4NK1FwE3AJtXWqx3quTHeqzEFeB8f\nf+9fRG5d4FpgPu5pfw+wVbX1b2Pd/wefOfAK8AHwHDAW+ER3r3uq15fxMLav4/4drwJ/AbbsCfUv\ncj/uSs//2d29/tT72uS3hd297rn6fRV4EHgv/c8/DtS0d/1jVapWIml14AHgc/iY9DeBzwL3p3Pd\nic1xI7UQeAh3hCjG3/E/8x/iU+BWBh6Q9KmOULKdGIX7K5yGR5q7Evg+8I+cXHesO/hMhSfxeu2D\n34fBwCO5wEHdtf4fI+kIYGuKP//dtf6Gr1ewU2bbOyfTXetOWvToZuAJ4ED8f3AC3utboH3qX+23\nkq66AScD/wU2zaQNTGk/qrZ+7Vjv43AP9/659ANS+u6ZtLVxf4FLq613G+q7QZG0gvd/TXeuexP3\n5HN4C+rHPaX++PTK14HDyLWcu2v98ZbzCnKzV3Iy3bLuqR4D8N6yE6tR/2g5t57hwKNmNruQYGZz\ncI/uA6qlVBUZDswzs4cKCeZz1G+jC98PM3u7SPITeDjXjdNxt6x7ExTm+xdmQ3yd7l//XwP/v717\njbGrrMI4/n9asRVQQIEqRGhJCIhESEEhKUJBtDaSoKKCCGIakYuioFK5WWJBoljjh2rggwGl1ICI\nEUIJUiDDRUFAaAvlIqXQSkttS2sEsRfa5Yf1Hub0zJnpdDoz3bP7/JITzr69Z69zyqx9efe75kXE\nzW2W1fn339xTKHWOvXEi0lN/owGL38m57z5M1zG4IZ/X3h7H4O7p+9hH0o5tlg1V48nLfc+U6drH\nLmmYpB0k7U/+sVpKFpiB/Pde2/glHUXetvpmN6vU/ffvqQ5CnWMfR/Yx+bKkBZLWS3pB0rlN6wxY\n/E7OfdfTuN67DfK+VMHmxjmvxXciaW/y2frZEfFkmb09xP43YC35lMPBwCciYmVZVtv4S2W7a4Gf\nRVP52hZ1jX+zdRCob+wAe5G3cK4mO4F+kuxr8ktJ55V1Biz+SgxCYjYUSNoJuI0sUzppM6vXzWnk\nvbT9gO8D90gaFxGLt+1uDbgfACPJP87blehFHYRtsmODZxiwM/DViLitzOuQNIZ8LHd6t1v204db\n36ym/VFRd0dSddfT99FYPmRJGkn2yhwNTIiIpU2Lax07QEQ8HxGPlXuux5N/tC4qi2sZf7l8ewk5\nTPBISbtIaozRP6JMD6Om8bdTrhY110Goc+yN/ib3tMy/GxglaRQDGL+Tc9/Np7P+c7OD6LwXuT3p\n6ftYHBFvDvL+9JtScOVWcjS7iRHR+vvWNvZ2IofUXUBnvfS6xr8fMAK4kfwju5q8XBnAheX9wdQ3\n/t6oc+zze7nOgMTv5Nx3twNHShrdmFHejyMvfW5vbgf2lvTxxgxJ7yF7Mw7Z76OMLPc7shPYiRHx\nWJvVahl7d8oZw4Fkgob6xv8kea/1WPL3b7xEjuU/nvwO6hp/F011EB4ps+oce2P0yAkt8ycCr0SO\naDlg8Q/68J11UXrhzSGrX/2wzJ4K7AQcMsSPGLuQdFJ5ezxwFjmG+QpgRUQ8UJLYQ+SoaZPJUdMu\nJs8sDokyHvpQI+kaMt4rgVkti1+JiCV1jR1A0h/JsejnkWPhHwCcD+wJHBERC+ocfzuSNgJXRsSU\nMl3L+Huog/AGcFhErKpr7A2S7iUHnrkMWEgWVZoEfC0iZgxo/Nv6Qe+h/Co/yC1sOnxnl2Et6/Ai\nB17Y0OZ1X9M6jWHsVpL/A9/NEB/Gjyym0i7uDcCUOsde4rqQfK57VYnrWXKUtNZBaGoZfzffyQZy\nXP5ax08m4jnk5fy1wCLgGmBU3WNvim1nsuNXY/jaOcDJgxG/z5zNzMwqxveczczMKsbJ2czMrGKc\nnM3MzCrGydnMzKxinJzNzMwqxsnZzMysYpyczczMKsbJ2ayCJJ0uaVHT9HxJZ/dy25cl3dA0fYik\ny5uKNgw6SSdKuqDN/GMkbZR09LbYL7OqcnI2q6axwOPwdqnKA4C/93Lb1pGFDiXL+723zbqD5bNA\nl+RMxnQkOUSomRVOzmbVdBidyXgsOWTk3D62Jbom7K1SKnVttYh4IyIejYg3+qM9s7pwcjarmFIj\n+FA6k/NHgWciYl0f2joDuK5MLiiXkDdI2qcsHy7pYknPSlojaYmkaZJGNLWxb9nuHEk/lbQEWFPq\nGe8u6VpJz0v6r6TFkmZK2qtp++uBM8jqPRvLa2FZNr7dZW1JF0h6TtJaSUslTZf07pZ1NkqaKuk8\nSQsl/UdSh6SDWtabIOkvkv4t6fXS7mVb+l2aDaZ+Ofo1s60n6SVg3zIZwJ2l6g1AlGpIAYyJiMW9\nbHYWWVHrUuAkoFEl59Xy35nAZ4CfAA8DHyrr7wt8saWtS8giGGcCw8lCAKPIogiXAMuBDwDfAx6S\ndGA5oJgK7AEcTpbSU9mmEecmZ/WSriKLLkwH7iBr415JVgc6pmWfTgOeB74NvBOYBvypfPZGSWPI\n0n2/B34ErAP2J2s1m1WWk7NZdUwkE8wZwKeAU8lE9iBZlrSjrLe0tw1GxEpJL5bJuRGxsLGs1KD9\nEnB6RMwss++TtBqYIekjETGvqbllEfH5lo/4B/CdpjaHAX8FFpd4bouIlyStANZF+3rYNG2/G/Bd\n4PqIaLQ7W9LKsk8nRMQdTZusB06IiA1le5GJ+GNkzeGxwA7AuU2Xzjt62gezKvBlbbOKiIjnSjL8\nINAREU8Bb5Jl626JiHnl9VY/feQE8gz21nJ5e7ik4cBs8qCgtQd12+Lx5XL3HEmvA2+RiTnITmxb\n6kgymc5smX9Tabv1zHl2IzEXT5V936dMzyET+M2STpK0Rx/2yWzQOTmbVYCkYSU5vgMYBzxcEuXR\n5KXo5WW6P+0JjCAPANY3vf5FJtf3taz/ass0ks4DfkXWsP0ceX/8CDJBjuzDPjV6lG/yWSUBv0bX\nHuerWqYbl8tHlu1eJA9CBNwALJP0sB/dsqrzZW2zariXzrPCAGYANzZNryfvOx8bEQ/002e+BvwP\nOIpMXq1aL5+36/F9MnBPRExuzJA0eiv2aVXZl/cDzza1OZw8WGhNxpsVEfcD90vagTzwuQK4Q9Lo\niNji9swGg8+czarhG2SHqWnAgvL+cGAF2ZnrcPKstLfPOjdrnE2+q2X+XeQZ5q4R8USb17JetL0j\neeDQbBJdE/naNp/fziNkp61TWuafQnZC6+hFG21FxPqI6ACuBnYCxvS1LbOB5jNnswqIiBcAJE0B\nZkXEk5IOAHYHrouI5VvR/DPk2ei3JP2WTKZzI+J+STcBf5D0C+BRYCOZtCYCkyNiwWbavguYLOni\nsv1xwBe62YczyyhnjwNrIuLpsuzts/aIWC3p58BFkt4E7iR7a18BPBgRs7YkcElnkbcG7gT+SfYa\nv4i8VfB0D5uabVNOzmYVUS67HkfewwX4NPBEHxLzJo8nRcQ8SZeTZ+dfJ6+YjQEWR8RXyn3jSeTj\nUGuBl4E/k/eem9tsZyqwC3A+eRbeQfY0X9iyza/Je9E/BnYFFtH5ONMmbUfEpZKWA2cD55CX339T\n9q/bOLvZ17nk93gVeY99Fdn7/dSIWNtmW7NKUES/DhxkZmZmW8n3nM3MzCrGydnMzKxinJzNzMwq\nxsnZzMysYpyczczMKsbJ2czMrGKcnM3MzCrGydnMzKxinJzNzMwq5v8CxQI05RMFzwAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x120899150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "print np.bincount(cluster_assignment)\n",
    "\n",
    "plot_heterogeneity(heterogeneity, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**. (True/False) The clustering objective (heterogeneity) is non-increasing for this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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",
   "execution_count": 648,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[19595 10427 29049]\n"
     ]
    }
   ],
   "source": [
    "print np.bincount(cluster_assignment)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beware of local maxima"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 650,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seed=000000, heterogeneity=57457.52442\n",
      "[18047  3824  5671  6983  1492  1730  3882  3449  7139  6854] 1492 18047\n",
      "seed=020000, heterogeneity=57533.20100\n",
      "[ 3142   768  3566  2277 15779  7278  6146  7964  6666  5485] 768 15779\n",
      "seed=040000, heterogeneity=57512.69257\n",
      "[ 5551  6623   186  2999  8487  3893  6807  2921  3472 18132] 186 18132\n",
      "seed=060000, heterogeneity=57466.97925\n",
      "[ 3014  3089  6681  3856  8080  7222  3424   424  5381 17900] 424 17900\n",
      "seed=080000, heterogeneity=57494.92990\n",
      "[17582  1785  7215  3314  6285   809  5930  6791  5536  3824] 809 17582\n",
      "seed=100000, heterogeneity=57484.42210\n",
      "[ 6618  1337  6191  2890 16969  4983  5242  3892  5562  5387] 1337 16969\n",
      "seed=120000, heterogeneity=57554.62410\n",
      "[ 6118  5841  4964  8423  4302  3183 16481  1608  5524  2627] 1608 16481\n",
      "252.686903\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "heterogeneity = {}\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",
    "    print('seed={0:06d}, heterogeneity={1:.5f}'.format(seed, heterogeneity[seed]))\n",
    "    sys.stdout.flush()\n",
    "    print np.bincount(cluster_assignment), np.min(np.bincount(cluster_assignment)), np.max(np.bincount(cluster_assignment))\n",
    "end = time.time()\n",
    "print(end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": "markdown",
   "metadata": {},
   "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",
   "execution_count": 617,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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 xrange(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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 618,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seed=000000, heterogeneity=57468.63808\n",
      "[19618  3857  4173  5219  1743  6900  5530  1348  4384  6299]\n",
      "seed=020000, heterogeneity=57486.94263\n",
      "[ 7132  1726  7419 20581  1941  5121  3777  3741   776  6857]\n",
      "seed=040000, heterogeneity=57454.35926\n",
      "[17602  3415  3535  1736  6445  2552  7106  7155   599  8926]\n",
      "seed=060000, heterogeneity=57530.43659\n",
      "[ 2655  9080  8047  2751 16942  1163  2810  6788  3389  5446]\n",
      "seed=080000, heterogeneity=57454.51852\n",
      "[ 6222  2024  6461  6765  5538 18352  3936  1563  1751  6459]\n",
      "seed=100000, heterogeneity=57471.56674\n",
      "[ 6459  5477  5746 17792  3436  1346  1743  3792  8075  5205]\n",
      "seed=120000, heterogeneity=57523.28839\n",
      "[ 3791  5186  1349 16549  7119  2625  8441  1131  6223  6657]\n",
      "309.646239042\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "heterogeneity_smart = {}\n",
    "start = time.time()\n",
    "for seed in [0, 20000, 40000, 60000, 80000, 100000, 120000]:\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",
    "    print np.bincount(cluster_assignment)\n",
    "\n",
    "end = time.time()\n",
    "print(end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 619,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAFSCAYAAADRt4TFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGC1JREFUeJzt3XuwZVV9J/Dvj8ZHiY80CU6NlqBRkvKRGR8JUwQjOCbq\noMZEaowxBJgKCDqlVc6UikOEDgoSJ8kYRscICkp8myrUUUEFuej4KCYaNGAQAz0xvqKkUQENIL3m\nj71vezh9bve9/bhn9b2fT9Wu22fvtfde+3fO3efbe69zbrXWAgDQs/3m3QEAgJ0RWACA7gksAED3\nBBYAoHsCCwDQPYEFAOje/vPuAMtXVT6DDsDctNZqXvt2hWUf01ozTU1nnHHG3PvQ26QmaqIuarKn\np3kTWACA7gksAED3BBb2eUcdddS8u9AdNdmemsymLttTkz5VD/elWJ6qap4vAOahqtIMugUAWJrA\nAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDu\nCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA\n6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNY\nAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9\ngQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0b/95dwD2hgMP\nTG6+ed69GLRUKm3e3WAP2bgx2bJl3r2A9adacyLdV1RV83wtT1XSTam66gy7y9PJelVVaa3VvPbv\nlhAA0D2BBQDonsACAHRvWYGlqjZV1daqEnAAgFW33ADSxgl2W9XcxmwBnXN+YCmumOwhVbW5qk5f\n4TonVNXWvdUnAFgrdjmwVNXTq+qWqjp3B21OGG8lHV5V762qH1bVd6rq1IltfLGqbq2qq6rq8TO2\n8Zyq+lxV3VZVN1fV+6rqIVNtfqeqLq+q7459+mJVHTdjW1ur6syqenFV3Tj2Z6GqHjXV7mlV9Zmq\n+v64veuq6g93tVY74MoVACzDLgWWMQx8MMnZrbWX7KDp4pvx25J8OclvJbk4ydlVdU6SP07y2iTP\nTXJAkouratuX2VXVKUn+Ksk1SY5J8oIkj0myUFUHTOzn4eN2j03y7CQfSnJ+Vb1gRp+OTXJ0kpck\nOSHJwUk+sDg+p6oeNh7bDWO/npXkT8f+AQBzsOJvuq2qlyd5dZJTWmsXLnO1i1prZ43rX5nkOUle\nmuTQ1trXx/kbknwgyeFJPj0GknOSvLW1dtLE/q9Kcn2SP0hybpK01s6eWF5JrkzyoCQvTHLeVF/u\nTPLM1tpdE+3fl+SwJJ9P8vgk90jyotbareM6CzPqsGHy4fhzv6n5rbW2dWKd/SbaJsmGGdvKYt8A\ngMFKr7C8PsmmJMdMhpWq2q+qNixOU+u0JJduezC8Gf99kusXw8rougxv5ou3ew5Pcr8k75ra9jfH\ntk+a2P8jqurdVfWNDIHkziQnJvnFGcfwialA8Lfjfg8eH189rv/eqjqmqg6a3kBVHTmxnzuT3JHk\nkCSnT82/bGrVy6eWv2Xc9922VVVPCgCwzUoCSyV5XoY3+Munlt2Qu7/hTo8fmf6rLncsMS9J7j3+\nfOC4z+k3+Tsy3Bb62SQZr8RcluSXkrw8yROT/HKSC5Lca8ZxTP8VkNsn99tauyHJ08Z9X5TkO+MY\nmskQ8dfjPianb2e4mjM57+Spfb1gavkfZQh0T5iY9ytJvjCj30mSTZs2bZsWFhaWata9TZuWnl+1\n/bTS9rA37YnXqPaz29OPhYWFu73nzNuy/pZQVZ2R4erBY5N8IslXkxzdWrttXP7o3D0cbG6t3VxV\nx2cIDoe21m6c2N4VSTa01iavkhySZHOSE1trF1TV05JckuS4JF+Z0a1bWmtfq6pfT/KxJE9srX1u\nYntvT3Jsa23DxLytSV7TWjt9Yt7ifk9orV00ddz3SHJEhltg/zbJQ1trM//sWVVtTnJha+3MWcuX\nWOf4JBdM9nEn7dfE3xIa/x7FXt5HR3/vpavOsLs8nXvXapwf2DXz/ltCKx3Dcm2So5J8MsklVfX0\n1tqPWmvX7vGeJZ9NckuGsPOOHbS7z/jzJ4szqmpjkt/c3Q601u7MMMD3dRnG1zws21+hAQD2shUP\num2tXTeO4VhI8vExtNy6k9VWrLV2S1W9LMkbquqBGa62/CDJg5McmeSK1tp78tNg88aq2pTkvklO\nS/K9JPdf6X6r6uQM42M+muQfkxyU5NQMY2eu2c3DAgB2wUrGsGy7Rtda+1qGN/WDk3ysqu67C/ue\ndc3vbvNaa+dluFLyCxnGk3wkyRkZPl1z9djmpgwfl96Q5P1JzkpyfpJ3LrH9ne33Sxmu2pyd4VbT\nuRnG6DyltXb7jHV3dDwAwB6wrDEs9MEYlpXso6NxBl11ht3l6dy7jGHp17zHsPhqfgCgewILq87/\nnoClOD+wFIEFAOiewAIAdE9gAQC6t+LvYYF9RS9f893ST1/YfRs3zrsHsD4JLKxJfY3ba76kB2A3\nuSUEAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gA\ngO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2B\nBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADd\nE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA\n0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiew\nAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7\nAgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEA\nuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQW\nAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRP\nYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA\n9wQWAKB7AgsA0D2BBQDonsACAHRPYAEAuiewAADdE1gAgO4JLABA9wQWAKB7AgsA0D2BBQDo3v7z\n7gAw24EHJjffPO9e9KGlUmnz7sa6tnFjsmXLvHvBelatOQnsK6qqeb7Wj6rE0z1SjLnzFFBVaa3V\nvPbvlhAA0D2BBQDonsACAHSvu8BSVZuqamtVddc3AGA+egwFbZwgVXMb3wWwz1tL59AeAwsAwN3s\nE4Glqp5eVbdU1bk7aHPCeCvp8Kp6b1X9sKq+U1WnTmzji1V1a1VdVVWPn7GN51TV56rqtqq6uare\nV1UPmWrzO1V1eVV9d+zTF6vquBnb2lpVZ1bVi6vqxrE/C1X1qKl2T6uqz1TV98ftXVdVf7jr1QKA\ntaf7wDKGgQ8mObu19pIdNF28jfS2JF9O8ltJLk5ydlWdk+SPk7w2yXOTHJDk4qra9sV5VXVKkr9K\nck2SY5K8IMljkixU1QET+3n4uN1jkzw7yYeSnF9VL5jRp2OTHJ3kJUlOSHJwkg8sjs+pqoeNx3bD\n2K9nJfnTsX8AwKjrb7qtqpcneXWSU1prFy5ztYtaa2eN61+Z5DlJXprk0Nba18f5G5J8IMnhST49\nBpJzkry1tXbSxP6vSnJ9kj9Icm6StNbOnlheSa5M8qAkL0xy3lRf7kzyzNbaXRPt35fksCSfT/L4\nJPdI8qLW2q3jOgvLPE4AWDd6vsLy+iSbkhwzGVaqar+q2rA4Ta3Tkly67cEQFP4+yfWLYWV0XZJK\nsni75/Ak90vyrqltf3Ns+6SJ/T+iqt5dVd/IEEjuTHJikl+ccQyfWAwro78d93vw+Pjqcf33VtUx\nVXXQzoqyadOmbdPCwsLOmq8JVcO0adPs5Zs2/bTN5LSvt4fe7Au/N9rfvf3uWFhYuNt7zrx199X8\nVXVGkjOS3JRkc5KjWms/nli+Ockh48OW5D+11i6qquOTXJDhSsqNE+2vSLKhtTYZOg4Zt31ia+2C\nqnp+kncs0aWW5FOttSePV2KuTXJrkrOT3JjkjiQvGvuxLUBV1dYkr2mtnT5jvye01i4a5x2Z5BVJ\njkxy7yRXJXlFa+1TM2qz7r6af/wq6Hl3Yy58FfoExZg7T8G+aU+eQ+f91fy93hJqSZ6S5BNJLq2q\no1trt43LnpnkXhNtN++B/f3z+PO4JF+ZsfyW8efhGa7KPLG19rnFhVV1j13dcWvtyiRXjts4IsMt\nsA9X1UNba/7UGACk38CSDFcyjkryySSXVNXTW2s/aq1duxf29dkMoeTQ1tpSV1qS5D7jz58szqiq\njUl+c3c70Fq7M8MA39dlGF/zsCQCCwCk78CS1tp14y2ThSQfH0PLrTtZbVf2c0tVvSzJG6rqgUku\nSfKDJA/OcKvmitbae/LTYPPGqtqU5L5JTkvyvST3X+l+q+rkDONjPprkH5MclOTUDGNnrtnNwwKA\nNaPXQbfbbri11r6W4U394CQfq6r77s72lprXWjsvw5WSX0hyUZKPZBhLsyHD4Ni01m7K8HHpDUne\nn+SsJOcneecS29/Zfr+U4arN2Uk+luGTSDckeUpr7fblHRoArH3dDbplaQbdri8GOU5QjLnzFOyb\n1tKg216vsECSrNuwArAnrKVzqMACAHRPYAEAuiewAADdE1gAgO51/T0ssN75m0KDFrWYt40b590D\n1juBBTq1hgb37wFt5pcaAeuHW0IAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDu\nCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA\n6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNY\nAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9\ngQUA6J7AAgB0T2ABALonsAAA3RNYAIDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA\n3RNYAIDuCSwAQPcEFgCgewILANA9gYV93sLCwry70B012Z6azKYu21OTPgks7POcXLanJttTk9nU\nZXtq0ieBBQDonsACAHSvWmvz7gPLVFWeLADmprVW89q3wAIAdM8tIQCgewILANA9gWUPqqojq2rr\njGnLRJsLl2iztaq+soNtnzq2+dRO+vC8sd3Xl1h+UlX9XVX9S1VdV1Un7/oRL88861JVD6qqC6rq\n2+Mx31hVZ81ot6p1mVdNqurAqvrzqrqhqn401uN/VtXPzWi7T9dkiTZ3VdW/mWpXVfXKqtpcVT+u\nqqur6jlL9HGf//1ZTl2q6tDxdXFtVd1SVd+qqg9O126i/bp4rUyts+bPtSupS63CuXb/lRaEnWpJ\nXpzkryfm/WTi32cmedPUOg9L8u4kH5y1war6+SSnJfmnHe24qh6Q5H8k+fYSy09K8hdJzkpyeZKn\nJPlfVZXW2pt3tO09YNXrUlWHJPlMkhvHff9TkocmecRUu3nVZR6vlf+d4fhfleS6JI9K8uokT0jy\nqxPbWSs1uSDJeVPzrp96/Jok/yXJf0vyxSTPS/L+qnpGa+3SxUZr7PdnZ3V5apKjxnZfSPKAJK9I\n8vmqOqK19jeLDdfZayXJujvX7rQuq3auba2Z9tCU5MgkdyX59ytc71Xjeo9cYvml4wvsiiSf2sF2\nzktySZILk3x9atmG8UV0wdT8tyb5bpINa60u4/LPJ9lvB/uYS13mUZMkhybZmuTEqfknj9s8dC3V\nZDzWM3ey7kFJ/iXJ6VPzL0ty9bxfJ3Osy4Ez5t0/yZYkb5t3XeZRk6n26+Jcu9y6rNa51i2hPW9X\nPvL1+0m+0Fr7u+02VvX8JI9L8sod7rTqiCTPT/Kfl2hyeJKfS/LOqfl/meRnkzxxhX1eqVWty3il\n4alJzm2tbd3BPuZZl9V+rdxz/PmDqfmLjxfPB2umJsvw9CT3yPbH+o4kvzT+zzFZY78/O9Na2zJj\n3g8z/M/6wROz19NrZdjpOjvX7nRnq3iuFVj2jndW1U+q6qaqemdVPWSphuOL/xFJ3jZj2c8k+bMk\nL2utfX8H29g/yZuTvK61duMSzR49/rxmav61GV7gj1pq+3vQatbliAyXRm+vqo+P90u3VNXbq+rA\niXbzrsuq1aS1dm2SK5O8qqqeUFUHVNVhGf539dHW2lfHpmuiJqMXjs/9bVV1eVVNnxQfleT21toN\nU/Onj3XeNUlWty6ztrkxyWOSTI5zmHddVrUm6+1cO9pZXVbtXGsMy571gyR/kuFN4YcZ/rd7WpLP\nVtXjWms3zVjnuCR3JHnPjGV/kuSrrbWLdrLfUzP87/mcHbRZfOHcPDV/y9TyvWEedXlQhl+Ct2ZI\n8Gdn+KU8J8kjkxw2tptXXeb1WnlGhnr834l5H07yHycer5Wa/GWGY/tWkkOSvCzJJ6vq11triwOS\nD0wyK+BNH+ta+v1ZTl1mecP4888n5q2n10qy/s61y6nL6p1r99a9NNO2+3OPS3Jnkj+asexe45P1\n/hnLfi3DvfVHTsybNS7hEUl+lOQ3JubNuq/6ygz3KO85NX9DhvuUp62xurxyPK6Lp+Y/d6zD03qr\ny96uyTj/XUm+keTEDJdgT8owcPDDPb5WdrUmS2zrvkn+X5IrJ+a9Ocm3ZrR9+Hisv9dbTVajLjPa\nLB7/8UvMn3tdVuG1sq7OtSuoy6qda90S2svaMJr++vw0ZU56dobR92+fsewvMiTWb1XVA8ZL/vsn\n2TA+XhyPcG6G0dZXTbS7Z4ZPaz6gqu49tltMtRun9rOYare7Z703rUJd/nn8ednU+h/P8L+Bx46P\nu6nL3q5JVT0jwydgjm2tvaW19n9aa+dnuH99dFU9a9zeWqjJrG3dmuQjSX5lYvbNSX5mRvPpY+2m\nJsmq1GWbqjolw6c6TmutTW+zm7qsQk3W27l21rZm1WXVzrUCy3wdn+SmDKPNpz0yySkZnuSbMzyZ\nR2QYuLRlXLbY7uipdr+bYWDclgyX55Kf3idcvI+4aPG+4ZLf6zEHe6Iu1y5zX/tKXfZETR6T4V7z\nF6bWv2piO8naqMlyXZvkXuPAwUmPzlCrr0y02xdqkuyZuiRJqur3k7wxyX9vrc26DbKv1GVP1GS9\nnWuXa/XOtat5aWo9Tkl+OcNn4M+Ymv/ADJfp/myJ9Z40Y/qbJF/KcAvgQWO7w2a0uyTDx8d+LcnP\nj+32z/DRsbdO7ectSb6XZP81VpcNGe67fmhq/d/NxMf+eqrLKtTk+PHYnzy1/lNz99sf+3xNltjW\n/ZP8Q5IrJuYdlOT2JK+aantZki9NPO6mJqtRl3H+b4/betMO1u2mLqvwWllX59oV1GXVzrWrVrT1\nMGUYcLQpw2W2Jyf5r+MTsTlT32uQ4Yuq7kry2BVsf+a4hBnttruvOs4/eXzhvjrD5/XPHB+fshbr\nkmEw2V0ZvpfkN5K8KMP/hC6bd13mUZMk98swfuUbGa66HJXkhRnGsGxOcp+1UpNx3TdluI9+ZIaw\n9uUMY31+dartazOMTXjp2PZN47H+h3m/TuZVlwxvxj/OMDj78CT/bmJ67NQ2181rZca6a/Zcu8Lf\noVU51+61wq3HKcMI8qszXC68PUMSfVOSfzWj7dWZ+GKqZW7/iuxgYNxEuwuT/MMSy07K8A2nP07y\n1SQnr+W6JPm98Zfsx0m+meT1mXhjnldd5lWTDJevz09yQ4Y36RsyjIH512upJkmemeTTGf5Hd/t4\n0r44yRNmtK0M33K7eTzWq5P89hLb3ad/f5ZblyRnZHgDmjXdOO+6zOu1MmPdNXuuXWldsgrn2ho3\nAADQLYNuAYDuCSwAQPcEFgCgewILANA9gQUA6J7AAgB0T2ABALonsAAA3RNYAIDu/X+xDlaiil4A\nzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x120a028d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.boxplot([heterogeneity.values(), 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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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",
   "execution_count": 620,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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 xrange(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_list)\n",
    "        \n",
    "        # Run k-means\n",
    "        # YOUR CODE HERE\n",
    "        centroids, cluster_assignment = kmeans(data, k, initial_centroids, maxiter=400,\n",
    "                                       record_heterogeneity=heterogeneity, verbose=True)\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to choose K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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 one hour to finish. We highly suggest that you use the arrays that we have computed for you.**\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",
   "execution_count": 621,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#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",
    "#start = time.time()\n",
    "#centroids = {}\n",
    "#cluster_assignment = {}\n",
    "#heterogeneity_values = []\n",
    "#k_list = [2, 10, 25, 50, 100]\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),\n",
    "#                                                               seed_list=seed_list,\n",
    "#                                                               verbose=True)\n",
    "#    score = compute_heterogeneity(tf_idf, k, centroids[k], cluster_assignment[k])\n",
    "#    heterogeneity_values.append(score)\n",
    "\n",
    "#plot_k_vs_heterogeneity(k_list, heterogeneity_values)\n",
    "\n",
    "#end = time.time()\n",
    "#print(end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the pre-computed NumPy arrays, first download kmeans-arrays.npz as mentioned in the reading for this assignment and load them with the following code.  Make sure the downloaded file is in the same directory as this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 622,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "10\n",
      "25\n",
      "50\n",
      "100\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAEKCAYAAADO/lZMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXmYFNXVh9/DqqKICCoSAXEHARHEFRhF1MTg+mncEDfE\nNUZNVFwQcUuiMTER4wZREKPirhhRBBTBnUXBBUUWBWWXfZ/z/XGqneqa6p7p2bp75rzPU0/Tt07d\nOnWHmV/de889V1QVx3Ecx3GyQ61sO+A4juM4NRkXYsdxHMfJIi7EjuM4jpNFXIgdx3EcJ4u4EDuO\n4zhOFnEhdhzHcZws4kLsOI5TCYjIeBEpzLYfTu7jQuzkPSLSUkQKReSVFOfvDs5PFpGmVe1fWUj8\nEReRxmlsCkXks4q4T3nqcFKiQFLbisjA4OfWLUs+OTlInWw74DiViYg8CFwCTAR+o6qrsuxSadHg\nqC73qYn0BraJlHl7O8VwIXaqJSJSG3gCOAt4EzhZVddl1yunJqGqP8QUS5U74uQ8PjTtVDtEpB7w\nAibCzwO9SiPCIjIkGDbslOL8/cH5glDZaSLyrogsEpF1IvKDiIwSkaMr6HEyRkTqicifRGSqiKwR\nkZ9F5K3ocGgwJN3N/imFoWNoxO5IEXldRJYEzzhDRK4PXnbCdn2C688VkRNFZKKIrBKRySGbBiJy\np4jMFJH1Qbs9JyLtUzxLDxGZFDzHwuBn1FhE5ojIdzH224nIHSLyZeDrEhF5Ma7+RB2BT/eLyPzA\np2kicmoKfzKpP2nYX0TGAQOCr4mph8LEc4jIbBFZLCJ1U9z7axFZJiL14847+Yv3iJ1qhYg0AF4G\njgQeBy7U0idUHwGcD5wNfBqptxZwOjBPVccHZZcD/wK+BZ4GVgG7AocDxwBjyvc0mRP8kX4r8OFj\n4GGgAXAi8LaInKaqLwXmA7HnbRH8O9Fbmxqq7wrgfmAx8BKwHDgCuBs4CPi/iAsKnAH0AF4B3iX4\nOyMiWwHvAAcCHwAjgd2wdv21iByrqu+F7v0b7Ge5HvvZLAV+HTxfHWBj5Nl3BCYA+wDjgNeAHYFT\ngZ4i0kNVP4z4WhcbMWkEPIcNJZ8BPCMix6nqmHLWH/6/95/gsxv2f3NO8H158PkYMAg4KWib8LN1\nBfYCBqvqBpzqhar64UdeH0BLLCjmXWASsAX4RxnqEeAHYAEgkXPHBve4O1T2KfA9UD+mrh3K+Uzj\nguf4C3BrzDEw8OezyHV3B9f9KVK+IzAbWBj2N3GfFD60wcRuIrBt5Ny/gvucEirrE/i0CegaU1/C\n50ci5UcG5TNDZbWAecAGoEPkZzQ6sP8uUs9/A59Oj5S3Bn4GpkXKZwf2zwN1QuVHBfW/Xs76i7Vt\n8LPbAnSLaZ9mQXu/EXPu8eC6A7L9++ZHxR9Zd8APP8p7hIR4S/A5phx13RPU0zNSPiwobxcq+xSY\nBdSrhGcaF3qmVEeSEAcitQz4PEWdlwfX/SZyn1RC/M/AvnPMue2Cc8+GyhJC/EyK+r4D1gE7xZx7\nNajviOB796Cup2Jsu0SFGHvR2Ay8VsLPtU2oLCHELWLsZwOLy1l/RkIcnH8Re5H5VaStVwOfVObv\nkR/ZO3xo2qlOTAeaAEeJyDWqel8Z6ngSuBYbnn4LQES2xoYLp6vq5yHbZ7Ae6HQReRoYD7yvFRcU\npkATVV0edzJm2dE+2BDrXBG5NeaSvTCx3hd4vRT37xL40EtEjo/eHhPVfWOu+zRaICLbAa2wF4dF\nMdeMB34DdADeCz4VG8KO8jEmimEOwnrRDVI8+37B577AF6Hyn1V1Xoz9D8AhFVB/pjyCTSOcD9we\nlJ2JDZk/Vo56nRzGhdipTswFTsP+qN8rIoWq+o9MKlDVaSIyAzhZRC5R1fXYH8ZtgeER27+KyFJs\nedRNwM3ABhF5Drg2heBkSiZRtok1x+2DIw7F5oxLW59gz5WK6PIcsOHvKA3TnAP4KbhXwm674HNx\n1FBVVUSWxPgKNv+abo1u9NlXpLDbTHIwa1nrz5Q3sCH5sBBfiL30PFXOup0cxaOmnWqFqs7E5hx/\nBO4TkT+UoZoRmPCeGHw/GxsK/W/M/Yao6kHATljg0ujA/pky3Le8rAw+n1HV2mmO29PWklyfAtuk\nqWvPyDWp1skmfNs5xb12jtitwoS5WAIWERFsqDiu/r+U8OzDKRuVXT9gLxnAEKCliBwlIm2x3vhz\nqroy/dVOvuJC7FQ7VPVrksX4qgyrGBF8niOW2eoY4B1VnZ/mnktV9QVVPQmLOu4WDMdWJV9iAtY5\nEKvSsAV+EbcoH2FieHB5HVNLpDIb2Fvis5sVYAKeiNieFnweEmN7EBbtHObj4Po4+4qgourfEnzW\nTmMzBHvxuzA4FB+Wrta4EDvVkkjP+O+ZiLGqfo8tUzkGC3Cqi80dJxFdlxuUbQ1sjw1tbgmVtxaR\nfaJrbysSVd0CPATsCfw5WHIV9a9LsIwowbLgc7eYKh/EBOFfItIspq6dRCRujjgVw4CtKBpyTdRT\nABwPfKuqE4PiCdg87aki0i5kWwtb4pOEqi7Elh91D5ZcFSPu51VaKrD+ZdjLTVx7J+61AJvDPwU4\nF2uXCRk77eQPVR0dRlE0ZPRYFrFrgyVlmI9FDE7HgmhqR+zi6toCtI/YCdAfeytfh715n5LCx75Y\n72I98BXQr6rbyY+M/k8loqZfiTm3d/B/qBC4KoM6+wbXbADWAg1jbJYH/5+expYZ3Q/MJGb5FLZm\ntJCYCN0U908sX2qcxiZu+VJ9bP3yluD/8CPAn7EXiS+D8p1C9pcE9XyCCeRNwG8j5zdhPe2ng7oe\nBt7GltpcF7JNRE2fm8LfrbCe5RZsSdRd2Lz7+uB3/PCI/fHBPVYFz3F38Hv7CSbS30bsGwfnt2AB\nYw8Cfw38ngOsjdjPJrIEKtr+5aw/ro79guvnB9feBFwec//fUvS37Lo4H/2oPkfV39CEeAtwGRaV\nmTgODNk0AxYBk7HF8gXBH4kthNZxBraF2LBNl8ixVcTuTkyArw58+HdQ33ERu75B+aDAblDw3cU4\nRw9MiLcAL6c4nxDjLcDvS1lnI0yAk5boRGz6YUkuvgtsF2Jrmc+KsZ2N9ZIzEeLNlCzE02LKawGX\nYmuqfwbWYElHXsDmr2uFbGtjAjcbe+nYAgyN1HcwNuc9HxPNBUHdN5G8zKZPcH2sEAc2DYA7sBeW\n9Vgw1khCy8Ii9j0w0V4TtO9jwA7YnO2UGPutgRswoVwV2H2NTTecGPMzmZWu/ctZf6o6emOCnvj/\nVexlIPgZLsJeRHbJ9u+YH5V7SPBDrzJEpDswFlunOTaFzcWYUO6jqt+Gyv+Lrb9rHiorBO5Q1QHF\na/rFpimWeOEuVR0UKh+DLQ85IPheG/sjM0pVLwjZDQF6Ac3Uhv8cx8kSIrIH8A32gnRGtv2pDESk\nBfaC96qqnpxtf5zKJVtzxCUFkiQCMaJLC1ZQNp+PC+ocESl/EmgnIi2D74di61CjdsOxKM0jynBv\nx3HKQJADukGkrB7wNyyA6aXYC6sHV2N/Jx/KtiNO5ZPNYK0RIrI5SJo+QkTCwQsjgSXAYBFpFSRa\nPxkbVrs3pq5Lg2Tta0TkbRGJCmYbYIOqzoqUz8D+s7cJvrcNPqeXYOc4TuWzL/CjiDwrIn8RkUew\nZBm9sLXi2VgiVmmISEMRuUFEHgKuAD5S1dHZ9supfLKR0GMFJqbvYPMrHbG5pkki0lFVl6jqIhE5\nDEv4nthhpRAYqKp/i9Q3HEu+vgCbK/wTMFZEjlbVdwObxthcWZRlofPhz2gmo6id4ziVz/dYHuhu\nWNatWtjfg4HAX7Wq59Uqnx2wALa1WDBcv+y641QVVS7EqjqV0O4uwAQRmYCtWbwSuFVEmmA5V1dj\nIfzLsETst4jIBlW9J1Rfn1BdE0XkFaxHezsWbOU4Th6ilpns/Gz7UVWo6lx8SWmNJCdSXKrqFBGZ\niUU7A1yPbc3WQouyybwrInWA20VkiKouS1HXahEZRfIv8HIsCjZKooe7LGQH9ma6MI1dMUSkur2d\nO47jOBFUNZO0s6UiV9++9seWFURTun2EBV1F0+qVxAygvoi0jpS3xYI+vgjZCUVzxQkSc8Npk7ln\nOwQ+n45bb7016z7ky+Ft5e3l7ZUbR2WRE0IsIp2xnWMSO638BOwhIttHTBPp5VKmGhSRhthi+PAG\n3W9gazLPjpifg+2oMzf4/j4WJBa1641tSj4Rx3Ecx6lAqnxoWkSGY3u4TsGCtQ7EFsh/j202Dhay\nfxbwlojcg4ngkVhmrRc0yPkrItdiveNx2FByq8BmZ2zrMABUdbGI3Af0F5HVWKKQM7BEIb1CdptF\n5BYsWnsBlqGoB3AecIWqRrdecxzHcZxykY054hmYCF6FbaH2E5bDdaAG876q+qGIdAUGAP/Atkab\ng0VLhveY/RrbJ/ZULL/vSmwv0/NVNbon6o1YJpzfA7sE156mqv8LG6nqw0GSkGuBP2Jbkl2uqg9X\nwLM7AQUFBdl2IW/wtsoMb6/M8PbKPlWeWau6IiLqbek4jlN9ERG0BgVrOY7jOE6NwIXYcRzHcbKI\nC7HjOI7jZBEXYsdxHMfJIi7EjuM4jpNFXIgdx3EcJ4u4EDuO4zhOFnEhdhzHcZws4kLsOI7jOFnE\nhdhxHMdxsogLseM4juNkERfiKsBTUDuO4zipcCGuZObMgeOOg+efz7YnjuM4Ti7iuy9VENHdlwoL\nYfBg6N8f1qyBnXeGL76Axo2z6KTjOI5TZipr9yUX4goiKsQzZ8L++8OmTUU2ffrA449XvW+O4zhO\n+fFtEPOMvfeGG29MLnviCRg9Ojv+OI7jOLmJ94griGiPGGDjRujUCaZPLypr0cK+b7ddFTvoOI7j\nlAvvEech9erBkCFQK9TK8+YV7yk7juM4NRcX4kqmSxe4+urksgcegPfey44/juM4Tm5R5UIsIt1F\npDDmWBaxayMiL4jIfBFZLSLTReRaEakdsasvIveIyAIRWSsik0Ska8x9RUT6i8hsEVknIlNF5JQU\nPvYVkS9FZL2IfCUi/crzzIMGQevWyWUXXgjr15enVsdxHKc6kK0esQJXAIeEjqMTJ0WkGTAeaAX8\nHvgt8CLwV+COSF1DgQuBm4HjgR+B0SLSPmJ3BzAA+CdwHPA+MFJEjgsbiUhf4CFgJHAs8CzwYHnE\neJtt4LHHkstmzjSBdhzHcWo2VR6sJSLdgbFAT1Udm8LmYuDfwD6q+m2o/L9AN1VtHnzvAEwBzlPV\nYUFZbWAG8JWqnhSUNQW+B+5S1UGh+sYATVT1gNC1C4BRqnpByG4I0AtopqpbUvhcLFgrSr9+8Mgj\nRd9r14aPPoIDD0x7meM4jpMDVLdgrZIepG7wuSJSvoJkn08ANmK9VgACoXwaOFZEEvUcF9Q5IlLf\nk0A7EWkZfD8UaBJjNxzYETiiBL/T8te/QvPmRd+3bLEh6vBaY8dxHKdmkc1grREisllElojICBHZ\nLXRuJLAEGCwirURkOxE5GTgbuDdk1waYrarR2dYZQD1gz5DdBlWdFWMnwXmAtsHn9BLsysT228O/\n/51cNnUq3HNPeWp1HMdx8plsCPEKTEwvAo4EBmHzw5NEpAmAqi4CDsOE77vgmpHAX1T1b6G6GgPL\nY+6xLHQ+8flzKe2IqTNqV2Z69YIzz0wuu+02+Oqr8tbsOI7j5CNVLsSqOlVVr1PVUao6QVUTwVO7\nAFcCBIL8IrAaOAUowIKtbhGR66ra54rm/vuhSZOi7xs32hD1ltjZZ8dxHKc6UyfbDgCo6hQRmQl0\nCYquB1oALVR1ZVD2rojUAQaJyGOqugzrubaIqTLRc030ZJcDjUppB7ADsDCNXSwDBw785d8FBQUU\nFBTE2jVtCv/8J5x1VlHZpEnw4INw5ZXp7uA4juNUFePHj2f8+PGVfp+cSXEpIjOAear6axH5H7CT\nqnaK2JyA9ZQPVdWPROQW4CagUXieWEQGYmLeUFU3iUhv4HFgL1X9LmR3HjAEaK2qc4P1x+8AR4cj\nuoNI73HAkar6Tgr/S4yaDqMKJ5wAr71WVNaggaW/bNWq1NU4juM4VUR1i5pOQkQ6A/sAHwRFPwF7\niMj2EdNDgs/5weerWFDWaaG6agOnA6NVNRGP/AawGQv2CnMOMF1V5wbf38eCxKJ2vYGlwMTMniw1\nIha41bBhUdmaNXDxxSbSjuM4Ts2gyoemRWQ4MAtb/7sSOBC4AVvn+6/A7CHgLOAtEbkHE8EjgWuB\nF1R1Pth8s4g8A/xDROoBs4HLsEQgv4REqepiEbkP6C8iq4HJwBnY3HOvkN3moJc9WEQWAGOAHsB5\nwBWqurki2+JXv7KI6X6hVCFvvWW7NJ13XkXeyXEcx8lVspHQ4wZMBFsC22C939eBgaq6MGTXBcuE\n1RFoCMwBngLuU9UNIbv6wJ2YcDcCpgHXqeqEyH0F6A/0xQLDvgZuU9UXY3zsi4l+S2BecM+HS3iu\njIamE6hCjx4wblxRWaNG8MUX0KxZxtU5juM4lURlDU3nzBxxvlNWIQaYNQvatYN164rKTjkFnn++\ngpxzHMdxyk21niOu6eyxB9wRyaD9wgvw3HPZ8cdxHMepOrxHXEGUp0cMtob4sMMs93SCnXaCL7+E\nxuVOI+I4juOUF+8RV3Nq14YhQ6Bu3aKyRYvgmmuy55PjOI5T+bgQ5xD77w833ZRc9sQT8MYb2fHH\ncRzHqXwyGpoOticcoaprKs+l/KS8Q9MJNm6ETp0ssUeCFi3s+3bblbt6x3Ecp4zkytD0v4EFIjJY\nRNpXtDMO1KtnQ9S1Qj+ZefOgf//s+eQ4juNUHpkK8R7Ag9hGDFNE5H0R6SMiW1W8azWXLl3g6quT\nywYPhgkT4u0dx3Gc/KVMUdPB5gsnA5dg2al+BoYBj6jqlxXpYL5QUUPTCdauhfbtbY1xgr33tv2L\nt966wm7jOI7jlJJcGZoGLBWkqo5U1R5YjujPgd8D00XkHRE5viKdrIlssw089lhy2cyZMGhQdvxx\nHMdxKocyR02LyHYichnwPNANmIrthFQHeEVEXDLKSUFBch5qsNzUkydnxR3HcRynEsh4aDrYKakf\nli+6DvAs8KCqfhiyuQX4g6ruWIG+5jQVPTSdYMUKaNsW5s8vKjvgAEv8EV5z7DiO41QuOTE0LSKT\ngQ+xnZAGAc1VtU9YhAPeAnaoGBdrNttvDw9HtpuYOtV6xo7jOE7+k+k64lewJUxvpOv+BVsSNgvt\n81vtqawecYKzz4annir6Xq+eCfJ++1XaLR3HcZwQObH7koh0Ayar6uqYc9sCB6rquxXoX95Q2UK8\neDG0aQNLlhSVHXqoLWmqXbvSbus4juME5MTQNDAOaJPi3D7BeacSaNoU/vWv5LL334cHH8yOP47j\nOE7FkKkQp3sTqA9sKYcvTgn87nfQq1dyWf/+MGdOVtxxHMdxKoASh6ZFpBXQOvg6BrgC+CpitjVw\nAdBBVfesWBfzg8oemk4wf74NUa9cWVTWsyeMHg1S4QMmjuM4ToKszRGLyK3ArUDYMOyIBt83A5er\n6qMV7WQ+UFVCDPDoo3DxxcllQ4fC+edXye0dx3FqJNkU4pZAK0xsxwKXA19EzDYAM1V1WUU7mC9U\npRCrQo8eMC40I9+oEXzxBTRrViUuOI7j1DiyFqylqnNV9R1VHY+tH34y+B4+PiitCItIdxEpjDmW\nhWz+k8KmUES+iNQXZ7MlujuUGP1FZLaIrBORqSJySgof+4rIlyKyXkS+EpF+cXbZQsR6xeGc0z//\nDJdfbiLtOI7j5A8ZBWsForuqAu6r2FzzIaHj6ND5QZFzhwBnBudejqlvaMT2UGBmxOYOYADwT+A4\n4H1gpIgcFzYSkb7AQ8BI4FiCzGG5JsZ77AF33JFc9uKL8Pzz2fHHcRzHKRulGZr+DjhZVaeJyGyS\n54qjqKruUUJ93bEh7p6qOrbUjlrazIHA/uEdnkSkELhDVQekubYp8D1wl6oOCpWPAZqo6gHB99rA\nAmCUql4QshsC9MKSlMRGhlfl0HSCLVvgsMMs3WWCnXaCL7+Exo2r1BXHcZxqTzbXEb8DrAz9O91R\n2mQeZXmQ3sCnZdxm8TigLjAiUv4k0C6YBwfrSTeJsRsO7AgcUYZ7Vxq1a1uQVjjn9KJFxfcydhzH\ncXKX0swRn6+qs4N/nxd8T3lkcO8RIrJZRJaIyAgR2S2VoYgcDuwJPJ7C5NJgPneNiLwtIlHBbANs\nUNVZkfIZ2EtBIklJ2+Bzegl2OUPbtnDzzcllw4bBG29kxx/HcRwnM8q8DWI5WAHcC1xE0eYRRwOT\nRKRJimvOBTYCT8ecGw5cBvQA+gKNgbFBOs4EjYGfY65dFjof/lxegl1OccMN0K5dclm/frCqImbz\nHcdxnEolYyEWkY4i8kLQk90sIgcG5XdFA5/iUNWpqnqdqo5S1Qmqmgie2gW4MuZ+9YHTgFfjIrOD\n3Z9GqupEVX0K6IrN896e6bPlK/XqwZAhUCv005w3z7JuOY7jOLlNnUyMgyHfMcB3wFNY5HOCQuAS\nIONBUVWdIiIzgS4xp08EtgeeKGVdq0VkFBAeJl8ONIoxT/Rwl4XswLZwXJjGLpaBAwf+8u+CggIK\nCgpK43KFcNBBcM01cO+9RWWDB1tazK5dq8wNx3GcasP48eMZP358pd8n092X3gOWAicBtbHh4s6q\nOjlYk/sPVW1RJkdEZgDzVPXXkfJRQGdg11QRyzF1DQbOV9Vtgu+9sfnlvVT1u5DdecAQoLWqzhWR\nrljQ2dHhiO4g0nsccKSqvpPinlUeNR1l7Vro0AG+/baobK+9bLvEbbbJnl+O4zjVgVzZfelA4N+B\n4kRVZwnQtCxOiEhnbPemDyLlOwHHACMyEOGGwG+BD0PFb2ApOM+OmJ8DTA/tm/w+9hxRu97YC8jE\n0viQLbbZxhJ9hPnmG+jUCT78MP4ax3EcJ7tkNDQNrAdS9a2aYYFYaRGR4cAsYAq2LOpA4AZsnW9k\noz/OwV4WhqWo61osmnocNpTcCrgW2JmiBCCo6mIRuQ/oLyKrgcnAGUABtj44Ybc5WK88WEQWYMPw\nPYDzgCtUdXNJz5dtCgosUOvhh4vKvvrK1htfdx0MHAj162fLO8dxHCdKpkPTr2BzrUcGRZuATsEc\n75vAElU9q4Q6bsBEsCUm6j8BrwMDVXVhxHYqQCLhRkxdvwWux3rT22PC/h6W4OPTiK0A/bHI6l2A\nr4HbVPXFmHr7YoLeEpgH3KeqD0ftItdkfWg6wcqVcPjhMD26CAtb7vTEE9ZLdhzHcUpP1jZ9iDjR\nARuenQM8B9yC9WI7AJ2Ag1T164p2Mh/IJSEGE+M//rH4UDVYIpCbbrKjXr2q981xHCcfyQkhDhw5\nELgH6IYFbBUCE4BrVHVKRTuYL+SaECcYPRouvND2MY5ywAHw+OMW4OU4juOkJ2eE+JcLRbYiSJSh\nqmsr1Ks8JFeFGGxnpquvNtGNUrcuDBgA11+fnCrTcRzHSSbnhNhJJpeFOMGrr8LFF8NPPxU/16mT\nzR23bVv8nOM4jpNDQiwirYHTgRbAVpHTqqoXVpBveUU+CDHAsmVw5ZXw1FPFz9WrB4MG2dxy7dpV\n75vjOE4ukxNCLCInYfvz1gIWARsiJqqqrSvOvfwhX4Q4wQsvwCWXwOLFxc8dcogNY++zT5W75TiO\nk7PkihB/DvwInK2qMX/Cay75JsRgInzZZfDcc8XPbbUV3HUXXHVVcg5rx3GcmkquCPEa4GRVfbOi\nHcl38lGIAVTh2WdNkJfFZNLu2hX+8x/YY4+q981xHCeXyJUUl18BO1a0E072ELGNIWbMgBNOKH5+\nwgRo3942kCgsrHr/HMdxqjuZCvF1wI1BwJZTjdhlF3jpJRg2DBpF9qlauxauuAJ69oQ5c7LinuM4\nTrUl06HpCcAeWK/4G4pvC6iq2r3i3Msf8nVoOo7586FvX/jf/4qf23Zb+Nvf7LxU+ACN4zhO7pIr\nQ9NbsBzNk4DFwffw4YOX1YDmzWHUKBgyBLbbLvnc6tW2qcRxx8H332fHP8dxnOqEJ/SoIKpTjzjM\nvHmWInPMmOLnGjaE+++HPn28d+w4TvUnV3rETg2jRQt4803497+hQYPkcytXwvnnW5DXggXZ8c9x\nHCffyViIRaS5iNwnIp+IyGwR2T8o/4OIHFzxLjrZRsSSf3z2GXSPiQB47TXYf38YMcKWQzmO4zil\nJyMhFpG2wOdAb2ABluYysZFeS+CqCvXOySlat4axY204euutk88tXw7nnAOnngoLF8Zf7ziO4xQn\n0x7x34Avgd2BU4DwWPkk4JAK8svJUWrVgt//HqZNg8MOK37+xRdt44iRI6veN8dxnHwkUyE+Aviz\nqq4GooOQC4FdKsQrJ+fZay949124916oXz/53NKlcPrpcMYZsGRJdvxzHMfJFzIV4nTLk5oA68rh\ni5Nn1K4N114LU6dCly7Fzz/zjPWOX3qp6n1zHMfJFzIV4o+A81OcOx2YWD53nHxk331h4kTbJKJu\n3eRzixbBySdD7942j+w4juMkk6kQ3w70EpE3sYAtBY4WkSeAk4E7S6pARLqLSGHMsSxk858UNoUi\n8kWkvvoico+ILBCRtSIySUS6xtxXRKR/EOm9TkSmisgpKXzsKyJfish6EflKRPpl1kw1jzp1oH9/\n+PRT6Nix+Pknn7Te8ahRVe+b4zhOLpOREKvqO8BJWLDWUCxY689AV+AkVf2wtFUBV2DBXYnj6ND5\nQZFzhwBnBudejtQ1FLgQuBk4HtumcbSItI/Y3QEMAP4JHAe8D4wUkePCRiLSF3gIGAkci+2//KCL\ncelo1w4+/BBuu83EOcyPP8JvfwsXXAArVmTHP8dxnFyjzJm1RGRPYCdgqap+ncF13YGxQE9VHZvB\ndbcAA4H9VfXLoKwDMAU4T1WHBWW1gRnAV6p6UlDWFPgeuEtVB4XqHAM0UdUDQtcuAEap6gUhuyFA\nL6CZqm5J4V+1zKxVHqZMsaxbn39e/NyvfmUpNI85pur9chzHKQs5l1lLVb9V1UmZiHCIsjxIb+DT\nhAgHnAC5WHm7AAAgAElEQVRsxHqtCb+2AE8Dx4pIYsbyOKAuMCJS55NAOxFpGXw/FAs6i9oNxza6\nOKIMftdYOnaETz6Bm26ywK4wP/wAxx5riUJWrcqOf47jOLlApgk9zk1znCMivUTkV6WsboSIbBaR\nJSIyQkR2S3Pfw4E9gccjp9oAs1V1faR8BpZoZM+Q3QZVnRVjJ8F5gLbB5/QS7JxSUq8e3HEHvP8+\n7Ldf8fMPP2z7HY8bV/W+OY7j5AJ1SjZJ4nGK1g+He7XhskIReQY4X1U3xtSxArgXeAdYCXQEbgIm\niUhHVY1beXou1vN9OlLeGIiLxV0WOp/4/LmUdsTUGbVzMuSgg2DyZLj1Vlt7XBhaCDdnDhx1lO15\n/Oc/F89p7TiOU53JdGj6cGAu8ADQHdg3+HwQmIcFS92ARVAPjKtAVaeq6nWqOkpVJ6hqInhqF+DK\nqL2I1AdOA15V1ej+x04esdVW8Je/wHvvWUKQKA88AAccYOcdx3FqCpn2iP8IPK2qN4bKZgITRGQV\ncLGqniwi2wNnAzfGVRJFVaeIyEwgJi0EJwLbA0/EnFuO5buOkui5LgvZNSqlHcAOWKawVHaxDBw4\n8Jd/FxQUUFBQkM68xnLooZYE5KabLG91OMbt22+hWze4+mob0o7mtHYcx6kqxo8fz/jx4yv9PhlF\nTQdie5Kqvh1z7mjgBVVtKCI9gddUtX6xSlLXPQOYp6q/jpSPAjoDu0YjloNI6puARuF5YhEZCFwP\nNFTVTSLSGxtW30tVvwvZnQcMAVqr6txg/fE7wNHhiO4g0nsccGSwhCvOf4+aLgPvvmtbKX73XfFz\n++wDTzwBB/ueXo7j5AC5EjW9AeiU4lwnbB43Ue+a0lYqIp2BfYAPIuU7AccAI1IsG3oVC8o6LXRN\nbSzL12hV3RQUvwFsxnrpYc4Bpqvq3OD7+8CSGLvewFI8c1iF062bbSBx2WXFz339tW0s0b8/bNhQ\n9b45juNUBZkOTY8EbhORLcBzwCJsLfFp2Jzw0MDuACB2WZOIDAdmYet/VwIHYvPK3wP/ipifg4n6\nsLi6VHVqEBj2DxGpB8wGLgNaUZQABFVdLCL3Af1FZDUwGTgDKMDWByfsNge97MEisgAYA/QAzgOu\nUNXN6RrHKRvbbguDB8Mpp1iyj3nzis4VFloA16uvWu+4U6rXQMdxnDwl06HprYFHCYlciKeAvqq6\nXkSOB1ap6rsxddyAiWBLYBvgJ+B1YKCqLozYTgVIJNxI4VN9LLXmWdg88DTgOlWdELEToD/QFwsM\n+xq4TVVfjKmzL3Bt4OM84D5VfTiVD8E1PjRdAaxcCX/8Izz6aPFztWvDjTfCzTfbsijHcZyqpLKG\npsuUWUtE9sbSTu6CpZT8qIyJPaoNLsQVy+jRcOGFMH9+8XMdOljvuEOHqvfLcZyaS04JsVMcF+KK\n5+ef4Q9/MNGNUqcODBgAN9xQfMcnx3GcyiBnhFhEtgEuwNYPN8aW9IwD/qOqNXY/YhfiyuOVV6Bf\nP/jpp+LnOnWCRx6BAw+ser8cx6lZ5IQQi8guwHhgbyyxx0/Y8HRLbM61IDrPW1NwIa5cli6F3/8e\nnnoq/nyXLtC3L5xxhgV/OY7jVDS5IsTDsK0BT1HViaHyw4DnsSVD51W0k/mAC3HV8PzzcOmlsHhx\n/Pltt4UzzzRR7twZpMJ/ZRzHqankihAvBq5X1aEx5y4E/qyqTSvQv7zBhbjqWLzY1h0/91x6uw4d\n4KKL4JxzoFFcXjXHcZwMyJWEHtti+/XG8UNw3nEqlaZN4dlnbe746KNT202bBldeCc2awbnnwoQJ\nyek0HcdxcoFMe8RTgRmqGs08lUjUsb+qdqxA//IG7xFnj+++g6FD7fjxx/S2++xjveQ+fUzQHcdx\nSkuuDE2fg2W5Gosl8PgRC9Y6Azga6K2qKcJpqjcuxNln82Z4/XVLBvL668lbLUapWxdOOsnmknv0\ngFqZjg05jlPjyAkhDhy5GBiEpbZMsBAYoKox+ZBqBi7EucUPP8Djj8Njj8HcueltW7Wy5CHnnw/N\nm1eFd47j5CM5I8SBM7WwTRoS64i/VtU0/Y/qjwtxblJYCGPGWC/55Zdh06bUtrVqwfHHWy/517+2\npCGO4zgJsi7EwaYKHwA3qOqbFe1IvuNCnPssWgTDhpkoz5yZ3nbXXa2HfOGFsPvuVeOf4zi5TdaF\nOHBiOXBqeK9ex3Ahzh9ULYL60UdtCdT69altRSwyu29fOPFE32zCcWoyuSLEzwLfqeoNFe1IvuNC\nnJ8sXw5PPmmi/Pnn6W2bNLFo6759LfracZyaRa4IcVfgSWxf4pewqOmkClT1u4p0MF9wIc5vVOHj\nj02Q//tfWLMmvX3XribI//d/sPXWVeOj4zjZJVeEOByQFXuhqtYur1P5iAtx9WHVKnj6aRPljz9O\nb9uokWXuuugi35bRcao7uSLEfUqyUdWYTeuqPy7E1ZNp02wJ1JNP2raM6TjooKKNJ7bbrmr8cxyn\n6sgJIXZS40JcvVm3zgK7HnsM3n03vW2DBrbxxEUX2a5QvvGE41QPckqIg3XEbYAdgU9UtYQZteqP\nC3HN4euvTZCfeCL1LlAJ2rWzXvI558AOO1SNf47jVA45I8QicjlwK9AEmyc+SFUni8hLwFhV/WdF\nO5kPuBDXPDZutCQhjz0Gb72VfkOJ+vUtsKtvX+jWzXvJjpOP5MTuSyLSF7gfi5g+HQg7NAE4tRR1\ndBeRwphjWYztISLyPxFZLiKrRWSaiJwesYmra4uItI/YiYj0F5HZIrJORKaKyCmpnlNEvhSR9SLy\nlYj0K7l1nJpGvXpw2mkwejTMmgU332yJQOLYsAFGjICCAth3X7jnHksw4jiOk2mw1pfAK6p6vYjU\nBjYBnYMe8fHAEFXdpYQ6umObRlwJfBI6tVlVJ4fsjgdeoGi51EZsOHylqg4L2RUCQ4FHIrf6TFXX\nh+zuBK4BbgQmYxtVXAwcr6pvhOz6Ag8BdwJvAz2Am4DLVPXhNM/lPWKHzZvhf/+ziOtRo0reeOLE\nE20uuWdP33jCcXKdnBiaFpH1wG9UdWyMEBcAb6jqViXUkRDinqkydInItsAs4ElVvbaE+gqBO1R1\nQBqbpsD3wF2qOihUPgZooqoHBN9rY/stj1LVC0J2Q4BeQDNV3ZLiHi7EThLz5xdtPDFnTnrbli2L\nNp741a+qwjvHcTIlJ4amgSVAqxTn9gHml7Kekh7kdGwO+r5S1lcSxwF1gRGR8ieBdiLSMvh+aHDf\nqN1wLDDtiAryx6kBNG8ON91kw9ZvvmnD2HXrxtvOnQsDBpgg9+oFr7xivWvHcao/mQrxa8AAEWkd\nKlMRaQJcjc0dl5YRIrJZRJaIyAgR2S107nBsV6f2IvKZiGwSkXkiMiCI2I5yaTCfu0ZE3haRqGC2\nATao6qxI+QzspaBN8L1t8Dm9BDvHKTW1atnQ87PPWi/53ntTp8gsLITXXrMh6xYtTMi/q5G56hyn\n5pCpEN8MbMCEagwWNf1P4EtgC7ZPcUmsAO4FLgKODK45GpgUCDrArkADrGc6FJunfRy4BbgnUt9w\n4LLApi+2NeNYEekWsmkMxKVjWBY6H/5cXoKd45SJpk3h2mvhyy9tPXLv3rBVismcH3+Eu+6CPfYo\nEvING6rWX8dxKp+yLF/aDvgDcCywE7AUeAP4u6quLJMTIh2Bj7A53FtFZDQmzteo6v0huweBC4Cm\nqroqRV3bYi8Kc1W1e1D2MNBLVXeN2O4BfAP0VtURItIfuAPYWlU3huwS8+G3qOqdKe7rc8ROmVi+\nHJ56ygK8pk1Lb9ukCZx7rgV47bdf1fjnOI6RK3PEqOoqVb1dVY9Q1b1V9VBVva2sIhzUOQWYCXQJ\nipYGn2Mipm9ic70ph4hVdTUwCjgoVLwcaBRjnujhLgvZAURTL0TtHKfC2GEHuPxymDIFPvrI1hpv\nu2287ZIlcN990KaNbTzxxBOwdm3V+us4TsVSJxNjEfkOOFlVi723i8j+2NKm1sWvzJgZFVBHtL76\nItI6sjtUW2x4/YuQnQTlC0N2CeH/gjQMHDjwl38XFBRQUFBQLqedmoWI5as+6CAT22eesV7yhx/G\n27/3nh1XXQVnn20CfsABVeuz41Rnxo8fz/jx4yv9PmXZfekQVf0o5lxn4MOy7L4UXPsBcLuq3iYi\nbYHPgT+q6n0hu4eBc7Ch6dh+gIg0DK79TlWPDMqaAj9gy5xuD9mOCerqEHyvgy1felVVLwzZPQac\niC1fio1l9aFpp7L4/HMT5OHDS954onNnG7Y+80xo2LBq/HOcmkKurCMuBA5W1WKbw4nIJcCdqrpj\nCXUMx9YITwFWAgcCNwCrgU6quiywG4otY7oNS8DRE7gWGJQQUxG5FtgTGIf1YFsFNnsDR6nqpNB9\n7wauwpJzJBJ69MXmjv8XsusHDAbuxobGe2BJQK5Q1YfSPJcLsVOprFsHL7xgovzOO+ltGzSA3/3O\neskHH+wpNR2nIsiaEIvI1djSJIDmwGIsy1WYrbF51KdV9ewS6rsBE8GWwDbAT8DrwEBVXRiyqwMM\nAPoAOwNzgAdU9YGQzW+B67E1zNtjwv4e1vP9NHJfAfpj4rsL8DVwm6q+GONjX0zQWwLzgPvSZdUK\nrnEhdqqMmTMtUcjjj5e88cT++xdtPNHY4/4dp8xkU4hPBE4KvvbBRDP6q78Bmz99LNWQcXXHhdjJ\nBhs3wquvWi/5zTdL3nji1FNNlLt3916y42RKrgxN/wcbGp5d0Y7kOy7ETraZMweGDrVjfgk57vba\ny+aS+/SBnXeuEvccJ+/JCSFOutDW6+4ILFDVTRXqVR7iQuzkCps3wxtv2ND1a6/Bltjs6EadOnDC\nCdZL7tkTamccauk4NYecEeJgXnYQ0CEoSuxH/Bi2H/FTFexjXuBC7OQiCxYUbTwxu4RxrBYtijae\n2G239LaOUxPJiYQeInIS8DK2+cP1JG/eMBubQ3YcJ0fYdVe48Ub49lsYM8YiqevVi7edNw9uvRVa\ntYLjj4eXXoJNNX6sy3Eqn0zniKcAn6rqRUFU80aKtkE8EXhQVZtXkq85jfeInXxhyRIYNswCvL76\nKr3tLrtYD/nCCy3ntePUZHJiaDrYj7iXqr4Vsx9xN+DNkvYjrq64EDv5hipMmmSC/Oyztk45HUcd\nZXPJJ59sEdiOU9PIiaFpbJ1ukxTnWlF8WZPjODmKCBx+uM0hL1gAgwenT5E5dqxl7GreHK65BqZP\nT79cynGc0pFpj3gE0A7oBqzCesSdsDXEE4CpqnpxJfiZ83iP2KkOqMKnn1pw11NPwarYPc6K2Hln\nOOIIO7p2hQ4dLBLbcaojuTI03QrbrlCxxB7nAs8B7bHMVp1VdUFFO5kPuBA71Y3Vq23I+tFH4YMP\nSndNgwZw6KFF4nzwwal3knKcfCMnhDhw5FdY/ufofsQDVPX7inYwX3Ahdqoz06dbL3nYMNs/ubTU\nrg0dOxYJ8+GHWwCY4+QjOSPETjwuxE5NYP1623hi2DCYMKFseyHvuWeRMB9xBOy9t6fbdPKDbOaa\nHpBBfRreZrAm4ULs1DQ2bYJp04r2RX7vPVi4sOTrojRpkizMHTumXuvsONkkm0JcGFOsJCfz+KW8\nLPsRVwdciJ2ajirMmpUszF9/nXk9W29tc8sJYT70UN9b2ckNsinEUWGtA6wDDsb29U1CVdNktq2+\nuBA7TnEWL4aJE4uE+dNPLRd2JtSqBe3bJ/eam9fItEFOtsmZOeJoIo+KdihfcSF2nJJZuxY++qhI\nmCdNKnmJVBytWiUL8377mWA7TmXiQpzjuBA7TuZs2QKff14kzBMmWHKRTNlhB4vITghz586e/cup\neFyIcxwXYscpP6owd26RKL/3HnzxReb11K8PBx1UlGjk0ENNrB2nPLgQ5zguxI5TOSxdakPYiV7z\nxx9nviuUCOy/f/JwdosWleOvU33JZrBW60hRbeBr4ERgRtReVb+rMO/yCBdix6ka1q2DTz4pEuaJ\nE2HFiszr2W23ZGFu29YSkDhOKrK9fClqJDFlAJS0fElEugPjYk79rKqNI7aHALcChwB1gVnAnar6\nbMimPnAHcDbQCJgKXK+qEyJ1CXADcDGwC/YyMUhVX4jxsS9wDbA7MAf4u6o+XMJzuRA7ThYoLIQZ\nM5KXTc2bl3k9228Phx1WJMwHHWRLqRwnQTaFuE8mFarqEyXU1x0YC1wJfBI6tTk81C0ixwMvAE8C\nI7G9j9sAK1V1WMhuBPBr4I/AbOCK4PshqvpZyO5OTFxvxJZdnYGJ8vGq+kbIri/wEHAn8DbQA7gJ\nuCydGLsQO07uMG9e8rKpzz/PfKeounUt6CucnnPHHSvHXyc/yJk54nLfsEiIe6rq2BQ222K93ydV\n9do0dXUApgDnJcQ5mMOeAXylqicFZU2B74G7VHVQ6PoxQBNVPSB07QJglKpeELIbAvQCmqVaJ+1C\n7Di5y88/w/vvFwnzhx/Chg2Z17PffsnD2bvv7uk5axLVTYjHAUenEeILgEeBFqo6P01dt2C91Uaq\nuj5UPhC4HmioqptEpDfwOLC3qs4K2Z0HDAFaq+pcETkCeAc4RlXfDtkVYL3jo1T1nRS+uBA7Tp6w\nYQNMnpw8nL1sWeb1NGuWLMzt2/s2kNWZyhLibC6BHyEim0VkiYiMEJHdQucOB5YB7UXkMxHZJCLz\nRGSAiIR9bgPMDotwwAygHrBnyG5DWIRDdhKcB2gbfE4vwc5xnDymfn1b0vSnP8HLL1sGsC++gEce\ngXPPtZ5uafjxRxg5Eq66Cjp1siVSxxwDgwbB2LGwZk3lPodTPcjGu9sK4F6s57kS6Ij1aieJSEdV\nXQLsCjQARgCDsDndo4FbsH2PE8PVjYG4TdmWhc4nPn8upR0xdUbtHMepRtSqZcPO++0Hffta2fz5\nyfPM06ZZYFg6Vq+Gt96yAywK+8ADi9YzH3447LRT5T6Lk39UuRCr6lQssjnBBBGZAHyEBXDdivXU\n6wP9VfX+wO5dEWkCXC4iA1W1DInxHMdxSkfz5nD66XYArFwJH3xQJMwffGBLqdKxZYute/74Y/j7\n361s772Th7P33NPnmWs6OTGboapTRGQm0CUoWhp8jomYvgn0w4aIP8R6rnHL8hM910RPdjm2tKk0\ndgA7AAvT2DmOU8No2NCGnY85xr5v2gRTpiTPMy9eXHI9M2faMXSofd9pp2RhPuAAi9h2ag45IcQx\nFEsUksbuJBHZKjJP3BZb7vRtyK6+iLSOJBxpi62H/iJkJ0F5WIgTc8Npk+0NHDjwl38XFBRQUFBQ\nysdwHCffqFsXunSx45prbHnUN98kC/M335Rcz6JF8MILdgBssw0cckiRMB9yCGy3XeU+ixPP+PHj\nGT9+fKXfp8qjpmOdEOkMfADcrqq3iUhb4HPgj6p6X8juYeAcoKmqrhWRA7D54z6qOjywqR1cOzOy\nfOkH4A5VvT1U35igrg7B9zrY8qVXVfXCkN1jWCaxZqoau4mbR007jhNl4cLkeebJk224OhNq1bJe\ncrjX3KxZ5fjrpKc6LV8ajq0RnoIFax2IZbxaDXRS1WWB3VDgdOA2TGx7YkFagyJi+l/gGOA6LKHH\nZcBvgENVdVrI7m7gKiwwLJHQoy/QS1X/F7LrBwwG7saGxntgSUCuUNWH0jyXC7HjOGlZs8bWMCeE\n+f33LcArU1q3Thbmfff1eeaqoDoJ8Q2YCLYEtgF+Al4HBqrqwpBdHWAA0AfYGUs1+YCqPhCprz6W\nBessbB54GnBdihSX/THxTaS4vE1VX4zxsS8m+i2BecB9nuLScZyKZvNm+Oyz5G0gf/op83p23DF5\nG8hOnaBevYr3t6ZTbYS4uuJC7DhOeVGF2bOT55m//DLzerbayuauE8J86KHQKC5c1ckIF+Icx4XY\ncZzKYMmS5HnmTz8t2zaQ7dolD2fvtlvJ1znJuBDnOC7EjuNUBWvX2rrkhDBPmmRrnDOlRYuiRCNH\nHAFt2lhgmJMaF+Icx4XYcZxssGULTJ+ePM88P2WG/tQ0apQ8z9y5sw1xO0W4EOc4LsSO4+QCqrYN\nZHieeXo0e34pqFfP9mROCPNhh0HjGp7k14U4x3EhdhwnV1m2LHkbyI8+go0bM6+nbdvkeeaWLWvW\nsikX4hzHhdhxnHxh/XoL+koI88SJsDxu+5wSaN48WZjbtbONLqorLsQ5jgux4zj5SmGhLZMKD2fP\nmZN5PdttZ0PYCWHu0sVSdlYXXIhzHBdix3GqEz/8kCzMn31m88+ZUKeOJRdJCPPhh0PTppXjb1Xg\nQpzjuBA7jlOdWbEieZ75ww9tiDtT9tkneTh7jz3yZ57ZhTjHcSF2HKcmsXGjbWIR7jUvXVrydVF2\n3rlIlLt2hQ4drCedi7gQ5zguxI7j1GRU4euvk4V51qzM62nQwFJyJsT54INh220r3t+y4EKc47gQ\nO47jJPPjj8npOadMscCwTKhdGzp2TJ5n3mWXyvG3JFyIcxwXYsdxnPSsWlV8G8i1azOvZ889k+eZ\n9967auaZXYhzHBdix3GczNi0CaZNS07PuWhR5vU0aZIszB07Vs42kC7EOY4LseM4TvlQtXnlhCi/\n9x7MnJl5PVtvbXPL4W0gGzYsv38uxDmOC7HjOE7Fs2hR8jzz5MmweXNmddSqBe3bJ/eamzfP3BcX\n4hzHhdhxHKfyWbPGcmWHt4FcvTrzelq1Shbm/fYreRtIF+Icx4XYcRyn6tm8GT7/PHme+ccfM69n\nhx2KbwNZv36yjQtxjuNC7DiOk31ULU92eD3zF19kXk/9+kXbQHbtavPMjRtXEyEWke7AuJhTP6tq\n48CmJTA7xkaBHVR1Zai+uFVpCnRU1c9CdgLcAFwM7AJ8DQxS1RdifOwLXAPsDswB/q6qD5fwXC7E\njuM4OcjSpTaEnRDmjz+2iO1MEAHVyhHibCUSU+BK4JNQWdz0+53Aq5GyVTF2Q4FHImXRWLs7MHG9\nEZgMnAGMFJHjVfWNhFEgwg8F934b6AE8GAxJpBVjx3EcJ/fYcUfo1csOgHXr4JNPkreBXLEifR2V\n2c/KVo94LNBTVcemsEn0iC9S1aEl1FcI3KGqA9LYNAW+B+5S1UGh8jFAE1U9IPheG1gAjFLVC0J2\nQ4BeQDNV3ZLiHt4jdhzHyUMKC2HGjOR55u+/j7OsnB5xCTFilUZV77VxHFAXGBEpfxJoFwg/wKFA\nkxi74cCOwBGV6aTjOI5T9dSqBe3awaWXwogRMG8ezJ1r/770UjtXmZm7siXEACNEZLOILBGRESKy\nW4zN3SKySUR+FpGXRWT/FHVdKiLrRWSNiLwtIlHBbANsUNVoCvIZ2EtBm+B72+Bzegl2TjkZP358\ntl3IG7ytMsPbKzO8veJp0QLOOgsefND2Yi7LzlKlJRtCvAK4F7gIOBIYBBwNTBKRJoHNBmyeth9Q\nAFwLtAMmisjekfqGA5dhc7l9gcbAWBHpFrJpDPwc48uy0Pnw5/IS7Jxy4r/8pcfbKjO8vTLD26t0\n7LBD5dVd5cFaqjoVmBoqmiAiE4CPsACuW1X1J0xcE0wUkdFYz/QmoE+ovj4Ru1ewHu3tQPfKeQrH\ncRzHqRiyOTT9C6o6BYty7pLG5gfgvXQ2gd1qYBRwUKh4OdAoxjzRw10WsgOIvvtE7RzHcRynYlDV\nnDiw3u7/SrAZBXxZiroGA2tD33sDW4DWEbvzgvKWwfeuQCFwVMSue1DePc091Q8//PDDj+p9VIb+\nZWsdcRIi0hnYB3g2jU0LLGq5WAKOiF1D4LfAh6HiN7B1ymdjQ9YJzgGmq+rc4Pv7wJLALry0qjew\nFJiY6r6VEdLuOI7jVH+qXIhFZDgwC5gCrAQOxDJefQ/8K7C5F+uBfoANB+8b2GwG7grVdS2wJ5ap\nayHQCgvs2hk4M2GnqotF5D6gv4ispiihRwG2Pjhht1lEbgEGi8gCYAwWBHYecIWqZrjnh+M4juOk\nJxsJPW7ARLAlsA3wE/A6MFBVFwY25wOXYCK7LdYbfRtLSflNqK7fAtdjventMWF/D0vw8WnkvgL0\nxyKrEykub1PVF2N87IsJektgHnCfZ9VyHMdxKoVszw3n8wH8CngOWxq1Ange2C3bfmX7AP4PeBF7\niVkLfIWNZGwbsWsEPAYsBlYDbwH7Z9v/bB/YVEoh9uLp7RXfRr8B3sFS3q7AVl0UeFvFttXhwGhs\n1HAl8Clwfk3/vwU0x0ZhJwFrgt+5FjF2pWoboD5wD5adcW1Qb9fS+JITUdP5iIhsjQ2J743NIZ8D\n7IWtYd46m77lANdi0wg3YFnNHgQuBd6M2L0GHANcDpyCZT8bJyK7Vp2ruYWInAm0xwJDonh7ASLS\nD3gJ+Bg4CXvxG4mNsCXwtgJEpB0mHHWw3A0nYy8tQ4J2TFAT22tP7P/OMuBd4n/noPRtMxS4ELgZ\nOB74ERgtIu1L9CTbbyX5egBXAZuA3UNlrYKyP2Tbvyy3zY4xZYnI9YLg+4nB924hm4bYNMQ/sv0M\nWWq3HYJf3t8R6RF7e/3yzC2x3saVaWy8rYqe+y5gPbB1pHwSMNHb65fnvTBogxaR8lK1DdAh+J09\nN1RWGxsNfKmk+3uPuOz0Aj5Q1dmJAlWdg0VWn5gtp3IBVY1LBvcxlia0efC9F7BAVd8NXbcS222r\nprbfX4DPVPWZmHPeXkbiD2a6mA1vqyLqAhtVdV2kfAVFeSROwNsrFaX9v3QCsJHQyh+1DYKeBo4V\nkbrpbuJCXHbaUjwnNdh6aM9JXZwCbOgnsUV3uvZrISLbxJyrtgT50c/Bhr/i8PYyDsd6GWeKyLdB\nLvpvRCScic/bqojHsVjVf4pIMxHZPghGPQq4L7Bpg7dXKkr7f6kNMFtV18fY1cOGwVPiQlx2GlM8\nJzXYfEMlZiXNP0SkOXAb8JZaFjVI335Qg9oweFt+CLhHVb9NYebtZeyKxWX8FRt27YnFHjwgIlcG\nNsloGbgAAARISURBVN5WAao6A8vpfzIwH2uXfwGXqOrIwMzbKzWlbZuS7NLuU5ATCT2c6ouINABe\nxoZtLijBvKZyPbAVoTXyTkpqYUsaz1XVl4Oy8SKyO7Y88V9Z8ywHEZE9sdUcnwMXY/PFJwIPi8h6\nVf1vNv1zDBfisrOc+DfFVG9GNQ4R2QqLOGyFBTssCJ1O136J89WeYPvPG7G5z62CNktkaasvIttj\nS3S8vYyl2DDfmEj5m9hc3M54W4W5G3sJPkGLEhKNC3a6ux/4L95e6Sht2ywHWqSxS7tPgQ9Nl50Z\nFO1fHKYNRfOgNRYRqYO9iR8I/FpVo22Srv3mqeraSnYxV2iNrT98EvtlXo790irwp+Df++PtlWBG\nKW28rYz9sQDAaFbAj4AdRWQnvL3SUdq2mQHsHrxIh2mLvQilmnICXIjLwyvAISLSKlEQ/PtwbCi2\nxhJkMXsKC9A6UVU/jjF7BWguIl1D1zXEohRrUvtNwebwjsTaK3EIttd2AfZL7O1lJDLhHRsp/zXw\ng1p2Pm+rIn4C2gcvxmEOwYapl+HtlY7Sts2rWFDWaSG72sDpwGhV3ZT2Ltlev5WvB5Y8YCYwDQtd\nPwHbZ/kbYJts+5fltvk3wTpY4ODI0TywEWyp11xs3eyxwHhs043m2X6GbB8UX0fs7VXUFm9jWY76\nYcFaj2JLmnp7WxVrq1ODtnkj+BvVE3ggKLunprdX0D6nhv5mXRJ875Zp22DD/EuxaaajsKyLa4EO\nJfqR7YbI5wNLcTmS5BSXxVKk1bQDmB38oscdA0J2idRxS7DUcW9SzdPqZdCGW7Bc6OEyby9rh22x\noKwfsV7dVOB33lYp2+tYbDe5hcHfqcnBS4zU9PYKxDfu79TYTNsGm2K6l6IUl+9TyhSXVb7pg+M4\njuM4RfgcseM4juNkERdix3Ecx8kiLsSO4ziOk0VciB3HcRwni7gQO47jOE4WcSF2HMdxnCziQuw4\njuM4WcSF2HGcUiEifUSkUERaR8oPEpFlIvKpiKTd7s1xnOK4EDuOkwlJGYBE5DDgLeAr4EhVTbvL\njOM4xXEhdhynTIhIdyyH8TTgGFVdmWWXHCcvcSF2HCdjRKQn8DrwIXCcqq7OskuOk7e4EDvO/7d3\nhyheRmEUh88LIu7D7UwwDWKx2a0Ws1gGi6DFRQhW9+Ie1BH0GP4LGL4BefnkefINp/24cOFy1FUu\n38N9TXLV9sfuHDg3IQaOmCQ3Sb4ledL2dnkPnJ4QA0d9TvI4yavtIfA/eLA9ADiVJnmZy9+2r2fm\ne9u3y5vg1IQYOKpJXiR5lOTNzNy2fbe8CU5LiIHD2nZmnid5mORmZn62/bC9C85IiIF7aftnZp7l\nEuP3M/Or7aflWXA6HmsB99b2d5LrJF+SfJyZp8uT4HSm7d2nAIB/wo0YABYJMQAsEmIAWCTEALBI\niAFgkRADwCIhBoBFQgwAi4QYABb9BRqA3nQxM4QqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1254ca250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "filename = 'kmeans-arrays.npz'\n",
    "\n",
    "heterogeneity_values = []\n",
    "k_list = [2, 10, 25, 50, 100]\n",
    "\n",
    "if os.path.exists(filename):\n",
    "    arrays = np.load(filename)\n",
    "    centroids = {}\n",
    "    cluster_assignment = {}\n",
    "    for k in k_list:\n",
    "        print k\n",
    "        sys.stdout.flush()\n",
    "        '''To save memory space, do not load the arrays from the file right away. We use\n",
    "           a technique known as lazy evaluation, where some expressions are not evaluated\n",
    "           until later. Any expression appearing inside a lambda function doesn't get\n",
    "           evaluated until the function is called.\n",
    "           Lazy evaluation is extremely important in memory-constrained setting, such as\n",
    "           an Amazon EC2 t2.micro instance.'''\n",
    "        centroids[k] = lambda k=k: arrays['centroids_{0:d}'.format(k)]\n",
    "        cluster_assignment[k] = lambda k=k: arrays['cluster_assignment_{0:d}'.format(k)]\n",
    "        score = compute_heterogeneity(tf_idf, k, centroids[k](), cluster_assignment[k]())\n",
    "        heterogeneity_values.append(score)\n",
    "    \n",
    "    plot_k_vs_heterogeneity(k_list, heterogeneity_values)\n",
    "\n",
    "else:\n",
    "    print('File not found. Skipping.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "source": [
    "## Visualize clusters of documents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
   "execution_count": 623,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def visualize_document_clusters(wiki, tf_idf, centroids, cluster_assignment, k, map_index_to_word, 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",
    "    \n",
    "    print('==========================================================')\n",
    "\n",
    "    # Visualize each cluster c\n",
    "    for c in xrange(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 xrange(5): # Print each word along with the TF-IDF weight\n",
    "            print('{0:s}:{1:.3f}'.format(map_index_to_word['category'][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",
    "            nearest_neighbors = distances.argsort()\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 xrange(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('==========================================================')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us first look at the 2 cluster case (K=2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 624,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========================================================\n",
      "Cluster 0     she:0.025 her:0.017 music:0.012 he:0.011 university:0.011 \n",
      "\n",
      "* Anita Kunz                                         0.97401\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",
      "* Janet Jackson                                      0.97472\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",
      "* Madonna (entertainer)                              0.97475\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",
      "* %C3%81ine Hyland                                   0.97536\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",
      "* Jane Fonda                                         0.97621\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",
      "* Christine Robertson                                0.97643\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",
      "* Pat Studdy-Clift                                   0.97643\n",
      "  pat studdyclift is an australian author specialising in historical fiction and nonfictionb\n",
      "  orn in 1925 she lived in gunnedah until she was sent to a boarding\n",
      "\n",
      "* Alexandra Potter                                   0.97646\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",
      "Cluster 1     league:0.040 season:0.036 team:0.029 football:0.029 played:0.028 \n",
      "\n",
      "* Todd Williams                                      0.95468\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",
      "* Gord Sherven                                       0.95622\n",
      "  gordon r sherven born august 21 1963 in gravelbourg saskatchewan and raised in mankota sas\n",
      "  katchewan is a retired canadian professional ice hockey forward who played\n",
      "\n",
      "* Justin Knoedler                                    0.95639\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",
      "* Chris Day                                          0.95648\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",
      "* Tony Smith (footballer, born 1957)                 0.95653\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",
      "* Ashley Prescott                                    0.95761\n",
      "  ashley prescott born 11 september 1972 is a former australian rules footballer he played w\n",
      "  ith the richmond and fremantle football clubs in the afl between\n",
      "\n",
      "* Leslie Lea                                         0.95802\n",
      "  leslie lea born 5 october 1942 in manchester is an english former professional footballer \n",
      "  he played as a midfielderlea began his professional career with blackpool\n",
      "\n",
      "* Tommy Anderson (footballer)                        0.95818\n",
      "  thomas cowan tommy anderson born 24 september 1934 in haddington is a scottish former prof\n",
      "  essional footballer he played as a forward and was noted for\n",
      "==========================================================\n"
     ]
    }
   ],
   "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both clusters have mixed content, although cluster 1 is much purer than cluster 0:\n",
    "* Cluster 0: artists, songwriters, professors, politicians, writers, etc.\n",
    "* Cluster 1: baseball players, hockey players, soccer (association football) players, etc.\n",
    "\n",
    "Top words of cluster 1 are all related to sports, whereas top words of cluster 0 show no clear pattern.\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",
   "execution_count": 625,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========================================================\n",
      "Cluster 0     film:0.020 art:0.014 he:0.011 book:0.010 television:0.010 \n",
      "\n",
      "* Wilson McLean                                      0.97479\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",
      "* Anton Hecht                                        0.97748\n",
      "  anton hecht is an english artist born in london in 2007 he asked musicians from around the\n",
      "   durham area to contribute to a soundtrack for\n",
      "\n",
      "* David Salle                                        0.97800\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",
      "* Vipin Sharma                                       0.97805\n",
      "  vipin sharma is an indian actor born in new delhi he is a graduate of national school of d\n",
      "  rama new delhi india and the canadian\n",
      "\n",
      "* Paul Swadel                                        0.97823\n",
      "  paul swadel is a new zealand film director and producerhe has directed and produced many s\n",
      "  uccessful short films which have screened in competition at cannes\n",
      "\n",
      "* Allan Stratton                                     0.97834\n",
      "  allan stratton born 1951 is a canadian playwright and novelistborn in stratford ontario st\n",
      "  ratton began his professional arts career while he was still in high\n",
      "\n",
      "* Bill Bennett (director)                            0.97848\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",
      "* Rafal Zielinski                                    0.97850\n",
      "  rafal zielinski born 1957 montreal is an independent filmmaker he is best known for direct\n",
      "  ing films such as fun sundance film festival special jury award\n",
      "==========================================================\n",
      "Cluster 1     league:0.052 rugby:0.044 club:0.042 cup:0.042 season:0.041 \n",
      "\n",
      "* Chris Day                                          0.93220\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",
      "* Gary Hooper                                        0.93481\n",
      "  gary hooper born 26 january 1988 is an english professional footballer who plays as a forw\n",
      "  ard for norwich cityhooper started his career at nonleague grays\n",
      "\n",
      "* Tony Smith (footballer, born 1957)                 0.93504\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",
      "* Jason Roberts (footballer)                         0.93527\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",
      "* Paul Robinson (footballer, born 1979)              0.93587\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.93732\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",
      "* Neil Grayson                                       0.93748\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",
      "* Sol Campbell                                       0.93759\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 2     championships:0.040 tour:0.037 championship:0.032 world:0.029 won:0.029 \n",
      "\n",
      "* Alessandra Aguilar                                 0.94505\n",
      "  alessandra aguilar born 1 july 1978 in lugo is a spanish longdistance runner who specialis\n",
      "  es in marathon running she represented her country in the event\n",
      "\n",
      "* Heather Samuel                                     0.94529\n",
      "  heather barbara samuel born 6 july 1970 is a retired sprinter from antigua and barbuda who\n",
      "   specialized in the 100 and 200 metres in 1990\n",
      "\n",
      "* Viola Kibiwot                                      0.94617\n",
      "  viola jelagat kibiwot born december 22 1983 in keiyo district is a runner from kenya who s\n",
      "  pecialises in the 1500 metres kibiwot won her first\n",
      "\n",
      "* Ayelech Worku                                      0.94636\n",
      "  ayelech worku born june 12 1979 is an ethiopian longdistance runner most known for winning\n",
      "   two world championships bronze medals on the 5000 metres she\n",
      "\n",
      "* Morhad Amdouni                                     0.94763\n",
      "  morhad amdouni born 21 january 1988 in portovecchio is a french middle and longdistance ru\n",
      "  nner he was european junior champion in track and cross country\n",
      "\n",
      "* Krisztina Papp                                     0.94776\n",
      "  krisztina papp born 17 december 1982 in eger is a hungarian long distance runner she is th\n",
      "  e national indoor record holder over 5000 mpapp began\n",
      "\n",
      "* Petra Lammert                                      0.94869\n",
      "  petra lammert born 3 march 1984 in freudenstadt badenwrttemberg is a former german shot pu\n",
      "  tter and current bobsledder she was the 2009 european indoor champion\n",
      "\n",
      "* Hasan Mahboob                                      0.94880\n",
      "  hasan mahboob ali born silas kirui on 31 december 1981 in kapsabet is a bahraini longdista\n",
      "  nce runner he became naturalized in bahrain and switched from\n",
      "==========================================================\n",
      "Cluster 3     baseball:0.110 league:0.103 major:0.052 games:0.047 season:0.045 \n",
      "\n",
      "* Steve Springer                                     0.89300\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",
      "* Dave Ford                                          0.89547\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",
      "* Todd Williams                                      0.89820\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.90035\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.90643\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",
      "* James Baldwin (baseball)                           0.90648\n",
      "  james j baldwin jr born july 15 1971 is a former major league baseball pitcher he batted a\n",
      "  nd threw righthanded in his 11season career he\n",
      "\n",
      "* Joe Strong                                         0.90655\n",
      "  joseph benjamin strong born september 9 1962 in fairfield california is a former major lea\n",
      "  gue baseball pitcher who played for the florida marlins from 2000\n",
      "\n",
      "* Javier L%C3%B3pez (baseball)                       0.90691\n",
      "  javier alfonso lpez born july 11 1977 is a puerto rican professional baseball pitcher for \n",
      "  the san francisco giants of major league baseball he is\n",
      "==========================================================\n",
      "Cluster 4     research:0.038 university:0.035 professor:0.032 science:0.023 institute:0.019 \n",
      "\n",
      "* Lawrence W. Green                                  0.95957\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.96057\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.96100\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.96323\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.96374\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",
      "* Stephen Park Turner                                0.96405\n",
      "  stephen turner is a researcher in social practice social and political theory and the phil\n",
      "  osophy of the social sciences he is graduate research professor in\n",
      "\n",
      "* Robert Bates (political scientist)                 0.96489\n",
      "  robert hinrichs bates born 1942 is an american political scientist he is eaton professor o\n",
      "  f the science of government in the departments of government and\n",
      "\n",
      "* Georg von Krogh                                    0.96505\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",
      "Cluster 5     football:0.076 coach:0.060 basketball:0.056 season:0.044 played:0.037 \n",
      "\n",
      "* Todd Curley                                        0.92731\n",
      "  todd curley born 14 january 1973 is a former australian rules footballer who played for co\n",
      "  llingwood and the western bulldogs in the australian football league\n",
      "\n",
      "* Ashley Prescott                                    0.92992\n",
      "  ashley prescott born 11 september 1972 is a former australian rules footballer he played w\n",
      "  ith the richmond and fremantle football clubs in the afl between\n",
      "\n",
      "* Pete Richardson                                    0.93204\n",
      "  pete richardson born october 17 1946 in youngstown ohio is a former american football defe\n",
      "  nsive back in the national football league and former college head\n",
      "\n",
      "* Nathan Brown (Australian footballer born 1976)     0.93561\n",
      "  nathan daniel brown born 14 august 1976 is an australian rules footballer who played for t\n",
      "  he melbourne demons in the australian football leaguehe was drafted\n",
      "\n",
      "* Earl Spalding                                      0.93654\n",
      "  earl spalding born 11 march 1965 in south perth is a former australian rules footballer wh\n",
      "  o played for melbourne and carlton in the victorian football\n",
      "\n",
      "* Bud Grant                                          0.93766\n",
      "  harry peter bud grant jr born may 20 1927 is a former american football and canadian footb\n",
      "  all head coach grant served as the head coach\n",
      "\n",
      "* Tyrone Wheatley                                    0.93885\n",
      "  tyrone anthony wheatley born january 19 1972 is the running backs coach of michigan and a \n",
      "  former professional american football player who played 10 seasons\n",
      "\n",
      "* Nick Salter                                        0.93916\n",
      "  nick salter born 30 july 1987 is an australian rules footballer who played for port adelai\n",
      "  de football club in the australian football league aflhe was\n",
      "==========================================================\n",
      "Cluster 6     she:0.138 her:0.089 actress:0.014 film:0.013 miss:0.012 \n",
      "\n",
      "* Lauren Royal                                       0.93445\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",
      "* Barbara Hershey                                    0.93496\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.93559\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",
      "* Jane Fonda                                         0.93759\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",
      "* Janine Shepherd                                    0.93833\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",
      "* Ellina Graypel                                     0.93847\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.93858\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",
      "* Melissa Hart (actress)                             0.93913\n",
      "  melissa hart is an american actress singer and teacher she made her broadway debut in 1966\n",
      "   as an ensemble member in jerry bocks the apple\n",
      "==========================================================\n",
      "Cluster 7     music:0.057 album:0.040 band:0.035 orchestra:0.023 released:0.022 \n",
      "\n",
      "* Brenton Broadstock                                 0.95722\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",
      "* Prince (musician)                                  0.96057\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",
      "* Will.i.am                                          0.96066\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",
      "* Tom Bancroft                                       0.96117\n",
      "  tom bancroft born 1967 london is a british jazz drummer and composer he began drumming age\n",
      "  d seven and started off playing jazz with his father\n",
      "\n",
      "* Julian Knowles                                     0.96152\n",
      "  julian knowles is an australian composer and performer specialising in new and emerging te\n",
      "  chnologies his creative work spans the fields of composition for theatre dance\n",
      "\n",
      "* Dan Siegel (musician)                              0.96223\n",
      "  dan siegel born in seattle washington is a pianist composer and record producer his earlie\n",
      "  r music has been described as new age while his more\n",
      "\n",
      "* Tony Mills (musician)                              0.96238\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",
      "* Don Robertson (composer)                           0.96249\n",
      "  don robertson born 1942 is an american composerdon robertson was born in 1942 in denver co\n",
      "  lorado and began studying music with conductor and pianist antonia\n",
      "==========================================================\n",
      "Cluster 8     hockey:0.216 nhl:0.134 ice:0.065 season:0.053 league:0.047 \n",
      "\n",
      "* Gord Sherven                                       0.83598\n",
      "  gordon r sherven born august 21 1963 in gravelbourg saskatchewan and raised in mankota sas\n",
      "  katchewan is a retired canadian professional ice hockey forward who played\n",
      "\n",
      "* Eric Brewer                                        0.83765\n",
      "  eric peter brewer born april 17 1979 is a canadian professional ice hockey defenceman for \n",
      "  the anaheim ducks of the national hockey league nhl he\n",
      "\n",
      "* Stephen Johns (ice hockey)                         0.84580\n",
      "  stephen johns born april 18 1992 is an american professional ice hockey defenceman he is c\n",
      "  urrently playing with the rockford icehogs of the american hockey\n",
      "\n",
      "* Mike Stevens (ice hockey, born 1965)               0.85320\n",
      "  mike stevens born december 30 1965 in kitchener ontario is a retired professional ice hock\n",
      "  ey player who played 23 games in the national hockey league\n",
      "\n",
      "* Tanner Glass                                       0.85484\n",
      "  tanner glass born november 29 1983 is a canadian professional ice hockey winger who plays \n",
      "  for the new york rangers of the national hockey league\n",
      "\n",
      "* Todd Strueby                                       0.86053\n",
      "  todd kenneth strueby born june 15 1963 in lanigan saskatchewan and raised in humboldt sask\n",
      "  atchewan is a retired canadian professional ice hockey centre who played\n",
      "\n",
      "* Steven King (ice hockey)                           0.86129\n",
      "  steven andrew king born july 22 1969 in east greenwich rhode island is a former ice hockey\n",
      "   forward who played professionally from 1991 to 2000\n",
      "\n",
      "* Don Jackson (ice hockey)                           0.86661\n",
      "  donald clinton jackson born september 2 1956 in minneapolis minnesota and bloomington minn\n",
      "  esota is an ice hockey coach and a retired professional ice hockey player\n",
      "==========================================================\n",
      "Cluster 9     party:0.028 election:0.025 minister:0.025 served:0.021 law:0.019 \n",
      "\n",
      "* Doug Lewis                                         0.96516\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",
      "* David Anderson (British Columbia politician)       0.96530\n",
      "  david a anderson pc oc born august 16 1937 in victoria british columbia is a former canadi\n",
      "  an cabinet minister educated at victoria college in victoria\n",
      "\n",
      "* Lucienne Robillard                                 0.96679\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",
      "* Bob Menendez                                       0.96686\n",
      "  robert bob menendez born january 1 1954 is the senior united states senator from new jerse\n",
      "  y he is a member of the democratic party first\n",
      "\n",
      "* Mal Sandon                                         0.96706\n",
      "  malcolm john mal sandon born 16 september 1945 is an australian politician he was an austr\n",
      "  alian labor party member of the victorian legislative council from\n",
      "\n",
      "* Roger Price (Australian politician)                0.96717\n",
      "  leo roger spurway price born 26 november 1945 is a former australian politician he was ele\n",
      "  cted as a member of the australian house of representatives\n",
      "\n",
      "* Maureen Lyster                                     0.96734\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",
      "* Don Bell                                           0.96739\n",
      "  donald h bell born march 10 1942 in new westminster british columbia is a canadian politic\n",
      "  ian he is currently serving as a councillor for the\n",
      "==========================================================\n"
     ]
    }
   ],
   "source": [
    "k = 10\n",
    "visualize_document_clusters(wiki, tf_idf, centroids[k](), cluster_assignment[k](), k, map_index_to_word)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clusters 0, 1, and 5 appear to be still mixed, but others are quite consistent in content.\n",
    "* Cluster 0: artists, actors, film directors, playwrights\n",
    "* Cluster 1: soccer (association football) players, rugby players\n",
    "* Cluster 2: track and field athletes\n",
    "* Cluster 3: baseball players\n",
    "* Cluster 4: professors, researchers, scholars\n",
    "* Cluster 5: Austrailian rules football players, American football players\n",
    "* Cluster 6: female figures from various fields\n",
    "* Cluster 7: composers, songwriters, singers, music producers\n",
    "* Cluster 8: ice hockey players\n",
    "* Cluster 9: politicians\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",
   "execution_count": 626,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([17602,  3415,  3535,  1736,  6445,  2552,  7106,  7155,   599,  8926])"
      ]
     },
     "execution_count": 626,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.bincount(cluster_assignment[10]())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**. Which of the 10 clusters above contains the greatest number of articles?\n",
    "\n",
    "1. Cluster 0: artists, actors, film directors, playwrights\n",
    "2. Cluster 4: professors, researchers, scholars\n",
    "3. Cluster 5: Austrailian rules football players, American football players\n",
    "4. Cluster 7: composers, songwriters, singers, music producers\n",
    "5. Cluster 9: politicians"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**. Which of the 10 clusters contains the least number of articles?\n",
    "\n",
    "1. Cluster 1: soccer (association football) players, rugby players\n",
    "2. Cluster 3: baseball players\n",
    "3. Cluster 6: female figures from various fields\n",
    "4. Cluster 7: composers, songwriters, singers, music producers\n",
    "5. Cluster 8: ice hockey players"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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",
   "execution_count": 627,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========================================================\n",
      "Cluster 0     law:0.077 district:0.048 court:0.046 republican:0.038 senate:0.038 \n",
      "==========================================================\n",
      "Cluster 1     research:0.054 professor:0.033 science:0.032 university:0.031 physics:0.029 \n",
      "==========================================================\n",
      "Cluster 2     hockey:0.216 nhl:0.134 ice:0.065 season:0.052 league:0.047 \n",
      "==========================================================\n",
      "Cluster 3     party:0.065 election:0.042 elected:0.031 parliament:0.027 member:0.023 \n",
      "==========================================================\n",
      "Cluster 4     board:0.025 president:0.023 chairman:0.022 business:0.022 executive:0.020 \n",
      "==========================================================\n",
      "Cluster 5     minister:0.160 prime:0.056 cabinet:0.044 party:0.043 election:0.042 \n",
      "==========================================================\n",
      "Cluster 6     university:0.044 professor:0.037 studies:0.035 history:0.034 philosophy:0.031 \n",
      "==========================================================\n",
      "Cluster 7     election:0.066 manitoba:0.058 liberal:0.051 party:0.045 riding:0.043 \n",
      "==========================================================\n",
      "Cluster 8     racing:0.095 formula:0.056 championship:0.054 race:0.052 poker:0.051 \n",
      "==========================================================\n",
      "Cluster 9     economics:0.146 economic:0.096 economist:0.053 policy:0.048 research:0.043 \n",
      "==========================================================\n",
      "Cluster 10     championships:0.075 olympics:0.050 marathon:0.048 metres:0.048 she:0.048 \n",
      "==========================================================\n",
      "Cluster 11     she:0.144 her:0.092 miss:0.016 actress:0.015 television:0.012 \n",
      "==========================================================\n",
      "Cluster 12     he:0.011 radio:0.009 show:0.009 that:0.009 his:0.009 \n",
      "==========================================================\n",
      "Cluster 13     baseball:0.109 league:0.104 major:0.052 games:0.047 season:0.045 \n",
      "==========================================================\n",
      "Cluster 14     art:0.144 museum:0.076 gallery:0.056 artist:0.033 arts:0.031 \n",
      "==========================================================\n",
      "Cluster 15     football:0.125 afl:0.060 nfl:0.051 season:0.049 played:0.045 \n",
      "==========================================================\n",
      "Cluster 16     music:0.097 jazz:0.061 piano:0.033 composer:0.029 orchestra:0.028 \n",
      "==========================================================\n",
      "Cluster 17     league:0.052 rugby:0.044 club:0.043 cup:0.042 season:0.042 \n",
      "==========================================================\n",
      "Cluster 18     poetry:0.055 novel:0.045 book:0.042 published:0.039 fiction:0.035 \n",
      "==========================================================\n",
      "Cluster 19     film:0.095 theatre:0.038 films:0.035 directed:0.029 television:0.028 \n",
      "==========================================================\n",
      "Cluster 20     album:0.064 band:0.049 music:0.037 released:0.033 song:0.025 \n",
      "==========================================================\n",
      "Cluster 21     bishop:0.075 air:0.066 force:0.048 church:0.047 command:0.045 \n",
      "==========================================================\n",
      "Cluster 22     orchestra:0.146 opera:0.116 symphony:0.106 conductor:0.077 music:0.064 \n",
      "==========================================================\n",
      "Cluster 23     basketball:0.120 coach:0.105 nba:0.065 head:0.042 season:0.040 \n",
      "==========================================================\n",
      "Cluster 24     tour:0.256 pga:0.213 golf:0.142 open:0.073 golfer:0.062 \n",
      "==========================================================\n"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Looking at the representative examples and top words, we classify each cluster as follows. Notice the bolded items, which indicate the appearance of a new theme.\n",
    "* Cluster 0: **lawyers, judges, legal scholars**\n",
    "* Cluster 1: **professors, researchers, scholars (natural and health sciences)**\n",
    "* Cluster 2: ice hockey players\n",
    "* Cluster 3: politicans\n",
    "* Cluster 4: **government officials**\n",
    "* Cluster 5: politicans\n",
    "* Cluster 6: **professors, researchers, scholars (social sciences and humanities)**\n",
    "* Cluster 7: Canadian politicians\n",
    "* Cluster 8: **car racers**\n",
    "* Cluster 9: **economists**\n",
    "* Cluster 10: track and field athletes\n",
    "* Cluster 11: females from various fields\n",
    "* Cluster 12: (mixed; no clear theme)\n",
    "* Cluster 13: baseball players\n",
    "* Cluster 14: **painters, sculptors, artists**\n",
    "* Cluster 15: Austrailian rules football players, American football players\n",
    "* Cluster 16: **musicians, composers**\n",
    "* Cluster 17: soccer (association football) players, rugby players\n",
    "* Cluster 18: **poets**\n",
    "* Cluster 19: **film directors, playwrights**\n",
    "* Cluster 20: **songwriters, singers, music producers**\n",
    "* Cluster 21: **generals of U.S. Air Force**\n",
    "* Cluster 22: **music directors, conductors**\n",
    "* Cluster 23: **basketball players**\n",
    "* Cluster 24: **golf players**\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",
   "execution_count": 628,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========================================================\n",
      "Cluster 0     brazilian:0.137 brazil:0.082 de:0.056 rio:0.053 paulo:0.050 \n",
      "==========================================================\n",
      "Cluster 1     bishop:0.170 diocese:0.085 archbishop:0.083 church:0.072 ordained:0.058 \n",
      "==========================================================\n",
      "Cluster 2     zealand:0.247 new:0.069 auckland:0.056 wellington:0.031 zealands:0.029 \n",
      "==========================================================\n",
      "Cluster 3     comics:0.181 comic:0.121 strip:0.042 graphic:0.036 book:0.034 \n",
      "==========================================================\n",
      "Cluster 4     puerto:0.309 rico:0.220 rican:0.066 juan:0.041 ricos:0.031 \n",
      "==========================================================\n",
      "Cluster 5     bbc:0.192 radio:0.127 presenter:0.054 show:0.046 news:0.042 \n",
      "==========================================================\n",
      "Cluster 6     senate:0.059 district:0.053 county:0.051 committee:0.049 state:0.044 \n",
      "==========================================================\n",
      "Cluster 7     labor:0.105 australian:0.099 liberal:0.071 election:0.067 seat:0.061 \n",
      "==========================================================\n",
      "Cluster 8     economics:0.065 university:0.048 research:0.045 professor:0.043 economic:0.043 \n",
      "==========================================================\n",
      "Cluster 9     foreign:0.086 ambassador:0.076 affairs:0.061 nations:0.053 united:0.040 \n",
      "==========================================================\n",
      "Cluster 10     she:0.188 her:0.052 women:0.026 womens:0.020 council:0.019 \n",
      "==========================================================\n",
      "Cluster 11     rowing:0.246 sculls:0.097 rower:0.081 olympics:0.073 championships:0.068 \n",
      "==========================================================\n",
      "Cluster 12     fashion:0.086 photography:0.085 photographer:0.057 photographs:0.038 art:0.025 \n",
      "==========================================================\n",
      "Cluster 13     republican:0.098 governor:0.051 district:0.044 election:0.043 senate:0.043 \n",
      "==========================================================\n",
      "Cluster 14     orchestra:0.227 symphony:0.177 philharmonic:0.084 music:0.080 conductor:0.057 \n",
      "==========================================================\n",
      "Cluster 15     air:0.375 force:0.242 command:0.106 commander:0.094 base:0.080 \n",
      "==========================================================\n",
      "Cluster 16     baseball:0.098 league:0.097 era:0.083 pitcher:0.083 pitched:0.075 \n",
      "==========================================================\n",
      "Cluster 17     church:0.114 theology:0.072 theological:0.066 seminary:0.047 christian:0.037 \n",
      "==========================================================\n",
      "Cluster 18     song:0.071 songs:0.043 music:0.041 album:0.030 singer:0.025 \n",
      "==========================================================\n",
      "Cluster 19     basketball:0.165 nba:0.113 points:0.067 season:0.044 rebounds:0.044 \n",
      "==========================================================\n",
      "Cluster 20     art:0.209 museum:0.186 gallery:0.082 arts:0.046 contemporary:0.044 \n",
      "==========================================================\n",
      "Cluster 21     poetry:0.213 poems:0.083 poet:0.069 poets:0.044 literary:0.040 \n",
      "==========================================================\n",
      "Cluster 22     guitar:0.215 guitarist:0.045 music:0.045 guitars:0.037 classical:0.028 \n",
      "==========================================================\n",
      "Cluster 23     novel:0.127 published:0.045 novels:0.044 book:0.039 fiction:0.030 \n",
      "==========================================================\n",
      "Cluster 24     jazz:0.205 music:0.048 band:0.034 pianist:0.025 recorded:0.023 \n",
      "==========================================================\n",
      "Cluster 25     polish:0.211 poland:0.097 warsaw:0.091 sejm:0.039 she:0.023 \n",
      "==========================================================\n",
      "Cluster 26     trinidad:0.259 tobago:0.178 calypso:0.058 caribbean:0.033 soca:0.027 \n",
      "==========================================================\n",
      "Cluster 27     tour:0.261 pga:0.220 golf:0.140 open:0.073 golfer:0.063 \n",
      "==========================================================\n",
      "Cluster 28     afl:0.177 football:0.128 australian:0.092 adelaide:0.064 season:0.062 \n",
      "==========================================================\n",
      "Cluster 29     skating:0.263 skater:0.107 speed:0.095 she:0.066 ice:0.060 \n",
      "==========================================================\n",
      "Cluster 30     party:0.073 election:0.035 elected:0.029 candidate:0.022 parliament:0.021 \n",
      "==========================================================\n",
      "Cluster 31     rugby:0.198 cup:0.049 against:0.046 played:0.045 wales:0.040 \n",
      "==========================================================\n",
      "Cluster 32     book:0.039 books:0.029 published:0.026 editor:0.021 author:0.017 \n",
      "==========================================================\n",
      "Cluster 33     piano:0.150 music:0.071 orchestra:0.056 competition:0.053 pianist:0.051 \n",
      "==========================================================\n",
      "Cluster 34     wrestling:0.299 wwe:0.163 wrestler:0.092 championship:0.079 tag:0.078 \n",
      "==========================================================\n",
      "Cluster 35     opera:0.269 she:0.067 la:0.041 sang:0.040 operatic:0.036 \n",
      "==========================================================\n",
      "Cluster 36     radio:0.080 show:0.069 host:0.038 sports:0.030 television:0.028 \n",
      "==========================================================\n",
      "Cluster 37     music:0.131 composition:0.038 composer:0.037 orchestra:0.026 ensemble:0.023 \n",
      "==========================================================\n",
      "Cluster 38     drummer:0.099 band:0.092 album:0.040 drums:0.039 rock:0.034 \n",
      "==========================================================\n",
      "Cluster 39     moore:0.306 moores:0.034 her:0.021 she:0.020 sports:0.012 \n",
      "==========================================================\n",
      "Cluster 40     computer:0.086 engineering:0.072 research:0.045 science:0.044 technology:0.042 \n",
      "==========================================================\n",
      "Cluster 41     minister:0.164 prime:0.068 cabinet:0.043 party:0.039 government:0.038 \n",
      "==========================================================\n",
      "Cluster 42     research:0.062 professor:0.035 university:0.034 science:0.031 psychology:0.030 \n",
      "==========================================================\n",
      "Cluster 43     news:0.127 anchor:0.062 reporter:0.059 she:0.045 correspondent:0.045 \n",
      "==========================================================\n",
      "Cluster 44     league:0.088 town:0.060 season:0.060 club:0.059 football:0.055 \n",
      "==========================================================\n",
      "Cluster 45     football:0.046 cup:0.044 club:0.042 team:0.041 league:0.033 \n",
      "==========================================================\n",
      "Cluster 46     football:0.108 vfl:0.099 australian:0.068 melbourne:0.067 goals:0.064 \n",
      "==========================================================\n",
      "Cluster 47     design:0.166 architecture:0.119 architectural:0.058 architects:0.038 architect:0.037 \n",
      "==========================================================\n",
      "Cluster 48     philosophy:0.227 philosophical:0.045 university:0.044 professor:0.041 philosopher:0.041 \n",
      "==========================================================\n",
      "Cluster 49     physics:0.121 mathematics:0.072 mathematical:0.060 theory:0.053 professor:0.043 \n",
      "==========================================================\n",
      "Cluster 50     baron:0.070 lord:0.060 lords:0.054 chairman:0.035 british:0.034 \n",
      "==========================================================\n",
      "Cluster 51     chef:0.143 food:0.136 restaurant:0.095 wine:0.086 cooking:0.064 \n",
      "==========================================================\n",
      "Cluster 52     fiction:0.138 stories:0.069 short:0.054 fantasy:0.048 writers:0.043 \n",
      "==========================================================\n",
      "Cluster 53     poker:0.477 wsop:0.121 event:0.091 limit:0.078 winnings:0.072 \n",
      "==========================================================\n",
      "Cluster 54     canadian:0.122 canada:0.068 toronto:0.053 ontario:0.049 curling:0.028 \n",
      "==========================================================\n",
      "Cluster 55     sri:0.282 lanka:0.183 lankan:0.094 colombo:0.046 ceylon:0.027 \n",
      "==========================================================\n",
      "Cluster 56     conductor:0.207 orchestra:0.136 conducting:0.087 music:0.080 symphony:0.073 \n",
      "==========================================================\n",
      "Cluster 57     prison:0.035 police:0.027 sentenced:0.026 court:0.025 convicted:0.023 \n",
      "==========================================================\n",
      "Cluster 58     blues:0.234 band:0.047 music:0.039 album:0.037 guitar:0.035 \n",
      "==========================================================\n",
      "Cluster 59     dj:0.093 hop:0.052 hip:0.051 music:0.048 album:0.037 \n",
      "==========================================================\n",
      "Cluster 60     de:0.127 la:0.059 el:0.035 mexico:0.026 y:0.025 \n",
      "==========================================================\n",
      "Cluster 61     jewish:0.193 rabbi:0.132 israel:0.052 hebrew:0.038 jews:0.032 \n",
      "==========================================================\n",
      "Cluster 62     ballet:0.362 dance:0.109 dancer:0.084 she:0.057 danced:0.044 \n",
      "==========================================================\n",
      "Cluster 63     hockey:0.220 nhl:0.138 ice:0.067 season:0.053 league:0.048 \n",
      "==========================================================\n",
      "Cluster 64     law:0.148 court:0.093 judge:0.071 district:0.051 justice:0.043 \n",
      "==========================================================\n",
      "Cluster 65     coach:0.205 head:0.086 basketball:0.059 coaching:0.052 football:0.046 \n",
      "==========================================================\n",
      "Cluster 66     armenian:0.278 armenia:0.168 yerevan:0.100 sargsyan:0.055 genocide:0.031 \n",
      "==========================================================\n",
      "Cluster 67     album:0.088 released:0.044 music:0.040 records:0.033 albums:0.027 \n",
      "==========================================================\n",
      "Cluster 68     she:0.158 her:0.152 music:0.020 album:0.016 singer:0.013 \n",
      "==========================================================\n",
      "Cluster 69     theatre:0.194 directed:0.034 production:0.031 play:0.029 actor:0.027 \n",
      "==========================================================\n",
      "Cluster 70     health:0.099 medical:0.089 medicine:0.086 research:0.039 clinical:0.039 \n",
      "==========================================================\n",
      "Cluster 71     european:0.145 parliament:0.115 party:0.053 member:0.049 committee:0.048 \n",
      "==========================================================\n",
      "Cluster 72     marathon:0.459 half:0.087 she:0.082 hours:0.063 championships:0.062 \n",
      "==========================================================\n",
      "Cluster 73     she:0.147 her:0.105 actress:0.098 film:0.063 role:0.054 \n",
      "==========================================================\n",
      "Cluster 74     she:0.101 her:0.065 women:0.012 show:0.010 television:0.009 \n",
      "==========================================================\n",
      "Cluster 75     lds:0.196 church:0.177 churchs:0.099 latterday:0.074 byu:0.073 \n",
      "==========================================================\n",
      "Cluster 76     quebec:0.242 qubcois:0.064 universit:0.061 minister:0.059 parti:0.051 \n",
      "==========================================================\n",
      "Cluster 77     film:0.233 festival:0.085 films:0.048 documentary:0.048 feature:0.045 \n",
      "==========================================================\n",
      "Cluster 78     hong:0.288 kong:0.268 chinese:0.068 china:0.037 wong:0.035 \n",
      "==========================================================\n",
      "Cluster 79     soccer:0.296 league:0.072 indoor:0.065 team:0.053 season:0.052 \n",
      "==========================================================\n",
      "Cluster 80     he:0.011 that:0.009 his:0.009 world:0.008 it:0.007 \n",
      "==========================================================\n",
      "Cluster 81     ireland:0.092 northern:0.072 election:0.072 irish:0.066 gael:0.054 \n",
      "==========================================================\n",
      "Cluster 82     comedy:0.048 series:0.047 actor:0.043 television:0.038 role:0.037 \n",
      "==========================================================\n",
      "Cluster 83     racing:0.128 formula:0.080 race:0.066 car:0.061 driver:0.055 \n",
      "==========================================================\n",
      "Cluster 84     election:0.096 manitoba:0.086 liberal:0.071 party:0.067 conservative:0.060 \n",
      "==========================================================\n",
      "Cluster 85     business:0.038 company:0.031 chairman:0.027 ceo:0.025 management:0.023 \n",
      "==========================================================\n",
      "Cluster 86     chess:0.414 grandmaster:0.085 olympiad:0.066 championship:0.064 fide:0.059 \n",
      "==========================================================\n",
      "Cluster 87     tennis:0.077 doubles:0.068 boxing:0.057 title:0.048 open:0.047 \n",
      "==========================================================\n",
      "Cluster 88     president:0.038 served:0.028 board:0.028 university:0.026 education:0.022 \n",
      "==========================================================\n",
      "Cluster 89     campaign:0.061 presidential:0.054 political:0.047 republican:0.037 bush:0.037 \n",
      "==========================================================\n",
      "Cluster 90     football:0.120 nfl:0.106 yards:0.081 bowl:0.052 quarterback:0.041 \n",
      "==========================================================\n",
      "Cluster 91     baseball:0.117 league:0.108 runs:0.061 major:0.052 batted:0.044 \n",
      "==========================================================\n",
      "Cluster 92     album:0.115 her:0.073 billboard:0.066 chart:0.064 singles:0.064 \n",
      "==========================================================\n",
      "Cluster 93     film:0.087 films:0.050 directed:0.029 television:0.024 actor:0.022 \n",
      "==========================================================\n",
      "Cluster 94     championships:0.106 metres:0.086 she:0.059 m:0.059 athletics:0.054 \n",
      "==========================================================\n",
      "Cluster 95     art:0.109 gallery:0.040 artist:0.036 paintings:0.032 painting:0.032 \n",
      "==========================================================\n",
      "Cluster 96     band:0.120 album:0.040 bands:0.035 bass:0.031 rock:0.030 \n",
      "==========================================================\n",
      "Cluster 97     miss:0.361 pageant:0.209 usa:0.127 she:0.110 teen:0.063 \n",
      "==========================================================\n",
      "Cluster 98     freestyle:0.155 swimming:0.120 m:0.119 swimmer:0.090 heat:0.075 \n",
      "==========================================================\n",
      "Cluster 99     army:0.081 commander:0.080 command:0.076 military:0.076 staff:0.058 \n",
      "==========================================================\n"
     ]
    }
   ],
   "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 ;)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class of soccer (association football) players has been broken into two clusters (44 and 45). Same goes for Austrialian rules football players (clusters 26 and 48). The class of baseball players have been also broken into two clusters (16 and 91).\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
   },
   "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?\n",
    "\n",
    "Hint: Use `cluster_assignment[100]()`, with the extra pair of parentheses for delayed loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 645,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29\n"
     ]
    }
   ],
   "source": [
    "count = 0\n",
    "for i in range(100):\n",
    "    if np.sum(cluster_assignment[100]()==i) < 236:\n",
    "        count = count + 1\n",
    "print count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 480,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = [[0., 0., 1.], [0.5, 0., 0.], [0.5, 1., 0.]]\n",
    "B = [[0.3333333333, 0.3333333333, 0.3333333333], [0.3333333333, 0.3333333333, 0.3333333333], [0.3333333333, 0.3333333333, 0.3333333333]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.array(A)\n",
    "b = np.array(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 482,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.33333333,  0.33333333,  0.33333333])"
      ]
     },
     "execution_count": 482,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = np.array([0.3333333333, .3333333333, .3333333333])\n",
    "c = 1 * c\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 487,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta = .85"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 484,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.33333333  0.16666667  0.5       ]\n",
      "[ 0.5         0.16666667  0.33333333]\n",
      "[ 0.33333333  0.25        0.41666667]\n",
      "[ 0.41666667  0.16666667  0.41666667]\n",
      "[ 0.41666667  0.20833333  0.375     ]\n",
      "[ 0.375       0.20833333  0.41666667]\n",
      "[ 0.41666667  0.1875      0.39583333]\n",
      "[ 0.39583333  0.20833333  0.39583333]\n",
      "[ 0.39583333  0.19791667  0.40625   ]\n",
      "[ 0.40625     0.19791667  0.39583333]\n",
      "[ 0.39583333  0.203125    0.40104167]\n",
      "[ 0.40104167  0.19791667  0.40104167]\n",
      "[ 0.40104167  0.20052083  0.3984375 ]\n",
      "[ 0.3984375   0.20052083  0.40104167]\n",
      "[ 0.40104167  0.19921875  0.39973958]\n",
      "[ 0.39973958  0.20052083  0.39973958]\n",
      "[ 0.39973958  0.19986979  0.40039062]\n",
      "[ 0.40039062  0.19986979  0.39973958]\n",
      "[ 0.39973958  0.20019531  0.4000651 ]\n",
      "[ 0.4000651   0.19986979  0.4000651 ]\n",
      "[ 0.4000651   0.20003255  0.39990234]\n",
      "[ 0.39990234  0.20003255  0.4000651 ]\n",
      "[ 0.4000651   0.19995117  0.39998372]\n",
      "[ 0.39998372  0.20003255  0.39998372]\n",
      "[ 0.39998372  0.19999186  0.40002441]\n",
      "[ 0.40002441  0.19999186  0.39998372]\n",
      "[ 0.39998372  0.20001221  0.40000407]\n",
      "[ 0.40000407  0.19999186  0.40000407]\n",
      "[ 0.40000407  0.20000203  0.3999939 ]\n",
      "[ 0.3999939   0.20000203  0.40000407]\n",
      "[ 0.40000407  0.19999695  0.39999898]\n",
      "[ 0.39999898  0.20000203  0.39999898]\n",
      "[ 0.39999898  0.19999949  0.40000153]\n",
      "[ 0.40000153  0.19999949  0.39999898]\n",
      "[ 0.39999898  0.20000076  0.40000025]\n",
      "[ 0.40000025  0.19999949  0.40000025]\n",
      "[ 0.40000025  0.20000013  0.39999962]\n",
      "[ 0.39999962  0.20000013  0.40000025]\n",
      "[ 0.40000025  0.19999981  0.39999994]\n",
      "[ 0.39999994  0.20000013  0.39999994]\n",
      "[ 0.39999994  0.19999997  0.4000001 ]\n",
      "[ 0.4000001   0.19999997  0.39999994]\n",
      "[ 0.39999994  0.20000005  0.40000002]\n",
      "[ 0.40000002  0.19999997  0.40000002]\n",
      "[ 0.40000002  0.20000001  0.39999998]\n",
      "[ 0.39999998  0.20000001  0.40000002]\n",
      "[ 0.40000002  0.19999999  0.4       ]\n",
      "[ 0.4         0.20000001  0.4       ]\n",
      "[ 0.4         0.2         0.40000001]\n",
      "[ 0.40000001  0.2         0.4       ]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n",
      "[ 0.4  0.2  0.4]\n"
     ]
    }
   ],
   "source": [
    "d = (beta * a) + ((1-beta)/ 3)\n",
    "for i in range (100):\n",
    "    m = np.dot(d , c)\n",
    "    print m\n",
    "    c = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 488,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.05 ,  0.05 ,  0.9  ],\n",
       "       [ 0.475,  0.05 ,  0.05 ],\n",
       "       [ 0.475,  0.9  ,  0.05 ]])"
      ]
     },
     "execution_count": 488,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = (beta * a) + ((1-beta)/ 3)\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.36125   ,  0.12041667,  0.24083333])"
      ]
     },
     "execution_count": 267,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.20470833,  0.15353125,  0.25588542])"
      ]
     },
     "execution_count": 269,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.2175026 ,  0.08700104,  0.2175026 ])"
      ]
     },
     "execution_count": 271,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.18487721,  0.09243861,  0.16638949])"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.14143107,  0.07857282,  0.15714563])"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.13357379,  0.0601082 ,  0.1268951 ])"
      ]
     },
     "execution_count": 277,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.10786083,  0.05676886,  0.10786083])"
      ]
     },
     "execution_count": 279,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.09168171,  0.04584085,  0.09409438])"
      ]
     },
     "execution_count": 281,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.07998023,  0.03896473,  0.07792945])"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.06624003,  0.0339916 ,  0.06711161])"
      ]
     },
     "execution_count": 285,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.dot(a, b)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:gl-env]",
   "language": "python",
   "name": "conda-env-gl-env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}