4 Clustering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Clustering is the task of gathering samples into groups of similar\n",
    "samples according to some predefined similarity or dissimilarity\n",
    "measure (such as the Euclidean distance).\n",
    "In this section we will explore a basic clustering task on some synthetic and real datasets.\n",
    "\n",
    "Here are some common applications of clustering algorithms:\n",
    "\n",
    "- Compression, in a data reduction sens\n",
    "- Can be used as a preprocessing step for recommender systems\n",
    "- Similarly:\n",
    "   - grouping related web news (e.g. Google News) and web search results\n",
    "   - grouping related stock quotes for investment portfolio management\n",
    "   - building customer profiles for market analysis\n",
    "- Building a code book of prototype samples for unsupervised feature extraction\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A closer look to some Clustering Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### K-Means\n",
    "\n",
    "To cluster well, we need to identify distinct groups such that the instances within a group are similar to each other but different from instances in other groups. \n",
    "\n",
    "One such algorithm is `k-means` clustering. \n",
    "\n",
    "With this algorithm, we specify the **number of desired clusters** `k`, and the algorithm will assign each instance to exactly one of these `k` clusters. \n",
    "\n",
    "It optimizes the grouping by minimizing the **within-cluster** variation (also known as *inertia*) such that the sum of the **within-cluster** variations across all k clusters is as small as possible.\n",
    "\n",
    "To speed up this clustering process, k-means randomly assigns each observation to one of the k clusters and then begins to reassign these observations to minimize the Euclidean distance between each observation and its cluster’s center point, or centroid. As a result, different runs of k-means - each with a randomized start —will result in slightly different clustering assignments of the observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Hierarchical Clustering\n",
    "\n",
    "An alternative clustering approach — one that does not require us to precommit to a particular number of clusters — is known as **hierarchical clustering**. \n",
    "\n",
    "One version of hierarchical clustering called **agglomerative clustering** uses a tree-based clustering method, and builds what is called a dendrogram. \n",
    "\n",
    "A dendrogram can be depicted graphically as an upside-down tree, where the leaves are at the bottom and the tree trunk is at the top.\n",
    "\n",
    "The leaves at the very bottom are individual instances in the dataset. \n",
    "\n",
    "Hierarchical clustering then joins the leaves together — as we move vertically up the upside-down tree — based on how similar they are to each other. \n",
    "\n",
    "Once the hierarchical clustering algorithm has finished running, we can view the dendrogram and determine where we want to cut the tree — the lower we cut, the more individual branches we are left with (i.e., more clusters)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### DBSCAN\n",
    "\n",
    "An even more powerful clustering algorithm (based on the **density of points**) is known as `DBSCAN` (*density-based spatial clustering of applications with noise*).\n",
    "\n",
    "Given all the instances we have in space, `DBSCAN` will group together those that are **packed closely** together, where close together is defined as a minimum number of instances that must exist within a certain distance. We specify both the minimum number of instances required and the distance.\n",
    "\n",
    "If an instance is within this specified distance of multiple clusters, it will be grouped with the cluster to which it is most densely located. \n",
    "\n",
    "Any instance that is not within this specified distance of another cluster is labeled an **outlier**.\n",
    "\n",
    "Unlike `k-means`, we do not need to prespecify the number of clusters. \n",
    "\n",
    "We can also have arbitrarily shaped clusters. \n",
    "\n",
    "**DBSCAN** is much less prone to the distortion typically caused by outliers in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Clustering Hands-on"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's start of with a very simple and obvious example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "X, y = make_blobs(random_state=42)\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAeFklEQVR4nO3df5AcZ53f8fdH6zWsDXids8BobZ18iU8OYJDwnjClIsHmh4xNsGJfwNSFI1wqOiigMMUpJ8GFM6mkrELHcVyZw1HAyV0VAZOzLRRskA12wuEqc5aQhCwsHcJnsHYdLO5YAaclXq2/+WNn5NnZ7pmemd7p3unPq2prZ7p7up/F4tvPfJ/v87QiAjMzG3zLim6AmZn1hwO+mVlFOOCbmVWEA76ZWUU44JuZVcQZRTeglfPOOy9WrVpVdDPMzJaMvXv3/iQiliftK3XAX7VqFXv27Cm6GWZmS4akH6btc0rHzKwiHPDNzCrCAd/MrCIc8M3MKsIB38ysIkpdpWPVsHPfBNt3H2FyapoVoyNs3rCajWvHim6W2cBxwLdC7dw3wdY7DzI9MwvAxNQ0W+88COCgb5Yzp3SsUNt3Hzkd7OumZ2bZvvtIQS0yG1zu4VuhJqemO9ru9I9Z9xzwrVArRkeYSAjuK0ZHFmzLkv7xDcEsnVM6VqjNG1YzMjw0b9vI8BCbN6xecGy79E/9hjAxNU3w7A1h576JeZ/ZuW+C9dvu56Itd7N+2/0L9psNKgd8K9TGtWPcfN2ljI2OIGBsdISbr7s0sVfeLv2TZTwg603BbBA5pWOF27h2LFPapV36J8t4QKubglM/Nujcw7clo136Jynv37y900Fis0GSOeBLuk3SU5Ieadj2jyTdJ+n7td/npnz2KklHJB2VtCWPhttgapVf37h2jOsvG2NIAmBI4vrLnv12sHnDaoaHNO98w0OaNx6Q5aZgNqg66eH/d+Cqpm1bgG9ExMXAN2rv55E0BHwaeBPwEuDtkl7SVWttoLXLr+/cN8EdeyeYjQBgNoI79k7Mz79H00mb3qd9S7jikuUeyLWBlzngR8Q3gb9v2nwt8Oe1138ObEz46DrgaEQ8FhFPA1+sfc5snnaDrln2zzwzP8LPPBPzBm2TBomvv2yMO/ZOeCDXBl6vg7YviognASLiSUkvTDhmDHii4f0x4FVpJ5S0CdgEsHLlyh6bZ0tJu/x6r/vrmgeJ12+7v6OBXNf621LVjyodJWxr/uL97I6IHcAOgPHx8dTjbOlrDpznjAwzNT2z4Lh6fr1dlU4nk7gadTKQ67V/bCnrtUrnx5JeDFD7/VTCMceACxveXwBM9nhdW+KS8vX/8PQpljV1DxoHXdtV6XQyiatRJwO5XvvHlrJeA/4u4J211+8EvpxwzMPAxZIuknQmcEPtc1ZhSYFzZjZ4psWga7tJWp1M4mrUyY3CZZ22lGVO6Uj6AvBa4DxJx4A/BLYBX5L0b4EfAf+qduwK4LMRcXVEnJL0PmA3MATcFhGH8v0zrGwa0zXnjAwjwdTJmdM576wBsj7o2hjUWwXwrJO4mj8DZMrLd5s2MisDRZQ3TT4+Ph579uwpuhmWQWOAHz1rmF/88tSCipm6keEhnju8jJ+eXJivTyLgb7ddk2Nru9ecw4e5vyfLNwmzfpC0NyLGk/Z5aQXrWXMQbBfIp2dmec4ZyxgZHlqQ1kmSZ++51wqbTr4NmJWNA771LCkf386J6Rk++bY1bN99hImpaURy6VaWQdes8qqw6SZtZFYGXkvHetbNgOWK0RE2rh3jwS1XMjY6khjsh6RcUyWusLGqc8C3nnWacmnutafdMGYjcu1Ju8LGqs4B33qWVNY4PCRGR4YRMDoyzLlnDaeWSp4zMpx4XsG85Q16fXBJrwun+cEpttQ5h28962Qgsz5o+sHb97NidIQrLlnOPzx9KvG8UTvnxrVjmfPvrQZlN29YnVhhk2WMwDNsbRC4LNP6JqmkMW2wttHj265h/bb7E+vfx0ZHeHDLlannby6Z7LZKJ8v1zcrAZZnWszwWDEsaNG0X7OtpnbyeZtVthY3z/zYInMO3tvJ6Dmw3wbGe1in6aVZ+cIoNAgd8ayut53zTrkMdDWKmBcek5VQbTU5NZ1rvJm3wN217J7pdmM2sTBzwra20HvLU9ExHvf60oPlbl69krEVPuV6z325hNKXcOdK2d6LbhdnMysQ5fGsrbcGwZq0eGgLtq3nSBl3rveh2+feplCUd0rZ3yjNsbalzwLe2ksoZ07TLl7cKmlnLO9MGkL2SpVlrDvjWVlIgPvn0qcRF0noNru160a3q4XupszerAgd8y6Q5ELdLv9SPyXtVyVall/V6+E4mgHnFS6sSB3zrSqf5+LxmprYrvcySZ/esWasqB3zrWqvgmtYT/9CXDpxeVqGbXnVann70rOyll1kmaJkNop7LMiWtlrS/4ednkm5sOua1kk40HPPRXq9r5dZqBcxeJm9t3rCa4aGFdZa/+OWpzOfyrFmrqp4DfkQciYg1EbEGuAw4CdyVcOhf1Y+LiP/Y63Wt3LIM3nazFv3GtWOcfebCL6b1Z9/20jZX89igy3vi1euAH0TED3M+ry0xSZOsknTTqz4xnVxXPzk1nWkJY8+atarKO+DfAHwhZd+rJR2Q9FVJL007gaRNkvZI2nP8+PGcm2f90jwzdShlums3veq0z5wzMpxpzR/PmrWqym15ZElnApPASyPix037XgA8ExG/kHQ18KmIuLjdOb088uDIsnRxL+dqxUsYW5W0Wh45zx7+m4DvNAd7gIj4WUT8ovb6HmBY0nk5XttKLs9edeO5svBgrNmcPMsy305KOkfS+cCPIyIkrWPuRvN3OV7bloA816KpnyvtwSSNPBhrNieXgC/pLOANwO82bHs3QETcCvwm8B5Jp4Bp4IYo86O2bMlo13v3YKzZs3IJ+BFxEviVpm23Nry+Bbglj2uZNWq1kueQ5MFYswZeD9+WtLSJWMPLxCfe+goHe7MGXlrBlrR6QP/Y/zp0evXO0ZFhbnrLSx3szZo44NuS5weTmGXjlI6ZWUU44JuZVYQDvplZRTjgm5lVhAO+mVlFOOCbmVWEA76ZWUU44JuZVYQDvplZRTjgm5lVhAO+mVlFOOCbmVWEA76ZWUXkEvAlPS7poKT9khY8dVxz/lTSUUnflfTKPK5rZmbZ5bk88hUR8ZOUfW8CLq79vAr4TO23mZn1Sb9SOtcCfxFzHgJGJb24T9c2MzPyC/gB3Ctpr6RNCfvHgCca3h+rbTMzsz7JK6WzPiImJb0QuE/S4Yj4ZsP+hQ8dnbtJLFC7YWwCWLlyZU7NMzOzXHr4ETFZ+/0UcBewrumQY8CFDe8vACZTzrUjIsYjYnz58uV5NM/MzMgh4Es6W9Lz66+BNwKPNB22C/jtWrXO5cCJiHiy12ubmXVr574J1m+7n4u23M36bfezc99E0U1adHmkdF4E3CWpfr7/ERFfk/RugIi4FbgHuBo4CpwE3pXDdc3MurJz3wRb7zzI9MwsABNT02y98yAAG9cO7vBizwE/Ih4DXpGw/daG1wG8t9drmZnlYfvuI6eDfd30zCzbdx9xwDczGySTU9Mdbd+5b4Ltu48wOTXNitERNm9YvSRvDF5awcwqZ8XoSObt9fTPxNQ0wbPpn6WY83fAN7PK2bxhNSPDQ/O2jQwPsXnD6gXHtkr/LDVO6ZhZ5dTTMVnSNJ2mf8rMAd/MKmnj2rFMefgVoyNMJAT3tLRQmTmlY2bWQifpn7JzD9/MrIVO0j9l54BvZktCkaWRWdM/ZeeAb2alV9WZsXlzDt/MSm+QSiOL5IBvZqU3SKWRRXLAN7PS62RmrKVzwDez0huk0sgiedDWzEpvsUojB2VRtKwc8M2spbIExbxLI6tY+eOAb2apBjkoLtaa+GW5QSZxDt/MUg1yOeRiVP6UfSnlPJ5pe6GkByQ9KumQpA8kHPNaSSck7a/9fLTX65rZ4ltK5ZCdPqN2MSp/yn6DzCOlcwr4UER8p/Yw872S7ouI7zUd91cR8eYcrmdmfVL0SpFZ0yNZUk/N57rikuXcsXdiXoDutfKn7DfIPJ5p+yTwZO31zyU9CowBzQHfzJaYzRtWzwuk0L9yyE7GD9r1rLfe+V2mZ545vW9iapo79k5w/WVjPHD4eNf59uabyOhZw/z05MyC48oyXyDXQVtJq4C1wLcTdr9a0gFgEvi9iDiUco5NwCaAlStX5tk8M+tQkStFdjKomtaDnpiaZvP/PMDMM7Fg3/TMLA8cPs6DW648va2eFsrytybdkIaXieEhMTP77PXKNF8gt4Av6XnAHcCNEfGzpt3fAX41In4h6WpgJ3Bx0nkiYgewA2B8fHzhfyUz66uiVorsJD2SlnoakhKDfdK5Oq1ISrohzTwTjI4Mc/ZzzujqBrnYFT65BHxJw8wF+89HxJ3N+xtvABFxj6Q/k3ReRPwkj+ub2dKSJbB1Mn6QlnpqDsitztVpmWbaDenE9Az7//CNLa+bpB8lsHlU6Qj4HPBoRPxxyjHn145D0rradf+u12ubWXmlVc1kLV3sZDmFjWvHuPm6SxkbHUHA2OjI6fdpVLtGXdLNBdIDe95VPv2o8Mmjh78eeAdwUNL+2rYPAysBIuJW4DeB90g6BUwDN0SE0zVmA6pVbzVrT7rT8YO01FNaDv+3Ll85r4JHQFJQSgvgeQ9o96PCJ48qnW8xd7NsdcwtwC29XsvMloZWQb2TwNbr+EH9szftOsTU9Fz1zLlnDXPNy1/MA4ePc9GWu1kxOsLJp08lBvvmbwFJ584r596PElgvrWBmuWsV1PMIbPUxgImpaYYkZiMYSwm4zTeNpG8faYLW+fM8B7T7UQLrpRXMLHet8ttJuXkxF3izzJBtHAMAmK1lh7MuY5D07SNNqzGAvKWNQ5SuSsfMrFGr3mpjKmRianpe7jytMqWxqgdB2ghglsXPsubEF6N+vl110mKXwLqHb2a5a9db3bh2jAe3XMnY6MiC3HlzZUpzVU+7co92AT1L6mgxetdlWFjNPXwzWxRZeqtZBnA7ScFA+4Ce9O2j0ZA0b/ZtXhZrOeZOOOCbWWGyDOB2UpbYmIZJS5/Ug+uNt+9PPMfsIlWMl2FhNad0zKwwWSZXZa3eaUzDtEufbFw7ljogm2WgttOlmFv9Hf1cWM0B38wKk6UyJemm0GhkeIg/edsaHtxy5bwB4XazVrt9MHq3ufgyPIjdKR0zK1S7XH/zBKdzRoaRYOrkTOpkpyzpk24nTnWbiy9y5dE6B3wzK71OyxVbjQ30uiJlL7n4olYerXNKx8wGTlr65IpLlvdcGjl61nBH28vEAd/MBk7a2MADh4/3vCJlWhFPq+KebgZ5F4MDvpkNnLS0TR6lkSemFz7CsNX2Mky4qnPAN7OB0irA5lEa2ek5+rHOfVYO+GY2UFoF2DxKIzs9RxkmXNW5SsfMBkqrAJtHaWSn5+jHOvdZOeCb2UBpF2DzKI3s5Bz9WOc+q1xSOpKuknRE0lFJWxL2S9Kf1vZ/V9Ir87iumVmzTlIu/aie6cc691n13MOXNAR8GngDcAx4WNKuiPhew2FvAi6u/bwK+Eztt5lZrrKmXFo9dzfvYFz0hKu6PFI664CjEfEYgKQvAtcCjQH/WuAvag8uf0jSqKQXR8STOVzfzGyeLAG2DMsV91seKZ0x4ImG98dq2zo9BgBJmyTtkbTn+PHjOTTPzGyhMlXP9EseAV8J25rnnGU5Zm5jxI6IGI+I8eXLl/fcODOzJGVYrrjf8gj4x4ALG95fAEx2cYyZWd+UYbnifssj4D8MXCzpIklnAjcAu5qO2QX8dq1a53LghPP3ZlakMlXP9EvPg7YRcUrS+4DdwBBwW0QckvTu2v5bgXuAq4GjwEngXb1e18ysV2WpnumXXCZeRcQ9zAX1xm23NrwO4L15XMvMzLrjtXTMzCrCAd/MrCIc8M3MKsIB38ysIhzwzcwqwgHfzKwiHPDNzCrCAd/MrCIc8M3MKsIB38ysIhzwzcwqwgHfzKwiHPDNzCrCAd/MrCIc8M3MKsIB38ysInp6AIqk7cC/AJ4GfgC8KyKmEo57HPg5MAuciojxXq5rZmad67WHfx/wsoh4OfA3wNYWx14REWsc7M3MitFTwI+IeyPiVO3tQ8AFvTfJzMwWQ545/N8BvpqyL4B7Je2VtKnVSSRtkrRH0p7jx4/n2Dwzs2prm8OX9HXg/IRdH4mIL9eO+QhwCvh8ymnWR8SkpBcC90k6HBHfTDowInYAOwDGx8cjw99gZmYZtA34EfH6VvslvRN4M/C6iEgM0BExWfv9lKS7gHVAYsA3M7PF0VNKR9JVwO8Db4mIkynHnC3p+fXXwBuBR3q5rpmZda7XHP4twPOZS9Psl3QrgKQVku6pHfMi4FuSDgB/DdwdEV/r8bpmZtahnurwI+KfpGyfBK6uvX4MeEUv1zEzs955pq2ZWUU44JuZVYQDvplZRTjgm5lVhAO+mVlFOOCbmVWEA76ZWUU44JuZVURPE6+qbue+CbbvPsLk1DQrRkfYvGE1G9eOFd0sM7NEDvhd2rlvgq13HmR6ZhaAialptt55EMBB38xKySmdLm3ffeR0sK+bnpll++4jBbXIzKw1B/wuTU5Nd7TdzKxoDvhdWjE60tF2M7OiOeB3afOG1YwMD83bNjI8xOYNqwtqkZlZax607VJ9YNZVOma2VDjg92Dj2jEHeDNbMpzSMTOriF6faXuTpIna4w33S7o65birJB2RdFTSll6uWQU7902wftv9XLTlbtZvu5+d+yaKbpKZDYA8UjqfjIg/StspaQj4NPAG4BjwsKRdEfG9HK49cDyhy8wWSz9SOuuAoxHxWEQ8DXwRuLYP112SPKHLzBZLHgH/fZK+K+k2Secm7B8Dnmh4f6y2LZGkTZL2SNpz/PjxHJq3tHhCl5ktlrYBX9LXJT2S8HMt8BngHwNrgCeBTySdImFbpF0vInZExHhEjC9fvjzjn1EuveTgPaHLzBZL2xx+RLw+y4kk/VfgKwm7jgEXNry/AJjM1LolqNcc/OYNq+d9Hjyhy8zy0WuVzosb3v5L4JGEwx4GLpZ0kaQzgRuAXb1ct8zScvAf+tKBTD3+jWvHuPm6SxkbHUHA2OgIN193qQdszaxnvVbpfFzSGuZSNI8DvwsgaQXw2Yi4OiJOSXofsBsYAm6LiEM9Xre00nLtszGXxcrS4/eELjNbDD0F/Ih4R8r2SeDqhvf3APf0cq2lYsXoCBNtBljrVTedBHU/bMXMeuWZtjlLWlQtSSdVN/VxgYmpaYJnvyV4QpaZdcJr6XQgSy+7eVG1ZdLpdE6jTqpuWtXmu5dvZlk54JMtkHdSfdOYg2/+XN2qX8ke8F2bb2Z5qHxKJ2u6JK2XfdOuQy1r7jeuHeOVK89ZcN0Hf/D3/MHOg5na6Np8M8tD5QJ+86Som3YdyrSUQVpvemp6pu3N4qHHfpr42S98+4nE7c38sBUzy0OlAn5Sb35qeibx2OYAn7U3nXSzSMrh17dnGXjNozbfK3CaWaVy+ElpmTTNAT5pBmyaxptFu8Daria/eXzhk29b0/FArVfgNDOoWA8/6yBnUrokqZd97lnDiZ8PYP22+/mDnQdPB9Y00zOz3Hj7/sRed17lmF6B08xgAHv4rSpuskyKAk6nS5LO9eCWK+ddK63XPzE1zecf+lH6KnEJxzf3uvMqx3SVj5nBgPXw2/WIs0yKGpL44O37WfOxe9n8lwda9q4be/1Jsgb7uuZed16B2lU+ZgYDFvDbpS6aA3TSus2zEQRz1Tczs/NDdlIaZOPaMR7ccmXiubrRGMzzCtSu8jEzGLCAn6VHXA/Qj2+7hk++bc3pnPyQsoXstGukBeFObwSN58krUHsFTjODAcvhp+Xo04Jx44zYi7bcnfkaSdLWsb/+sjG+cuDJ1PLPRsPLNC+YNy/T0MuiaV6B08wGKuD38vCQLAO6rc7VKjj/p42XzhsATltf53nPPaPlMg1mZr1QpEwKKoPx8fHYs2dPR5/pdhnhpIqb4WXiec89g6mTM7kuSXzRlrsTB3QF/O22a3o+v5lVl6S9ETGetG+gevjQfY84z/RJO52mnszM8tBTwJd0O1DPcYwCUxGxJuG4x4GfA7PAqbS7T9H6lT7xc2vNrAi9PvHqbfXXkj4BnGhx+BUR8ZNerjco+vltwsysLpeUjiQBbwWubHfsoMk6ZtBu1q6Z2WLLK4f/GuDHEfH9lP0B3CspgP8SETtyum6uOh3wzbooWbvjOrmun21rZt1qG/AlfR04P2HXRyLiy7XXbwe+0OI06yNiUtILgfskHY6Ib6ZcbxOwCWDlypXtmpdJ3k+0qsu61k27GcBZr+tVL82sF21n2kbE6yPiZQk/XwaQdAZwHXB7i3NM1n4/BdwFrGtx7I6IGI+I8eXLl3f69yzQ6xOtWq0omXWtm1bHdXJdr3ppZr3IY2mF1wOHI+JY0k5JZ0t6fv018EbgkRyum0nWINnNQmWt1rppfODIspRlG1aMjnR0Xa96aWa9yCPg30BTOkfSCkn31N6+CPiWpAPAXwN3R8TXcrhuJlmDZDcLlaWtdXPFJcvnfatImlVbL8Ps5Lpe9dLMetFzwI+IfxMRtzZtm4yIq2uvH4uIV9R+XhoR/7nXa3YiLRguk+aldbpZqCxtUbIHDh9v+WQsAddfNlfz38l1veqlmfVi4GbaNkt7NOFsxLwBz25r45Mma33w9v0tPxPAA4ePn/581uu6ft/MejHwAb8eDD/0pQMLUivdPD0qiywLsTUv2Zy1DV5Mzcy6NVDr4afZuHaMZ1IWiasH3ryeHwvZnqzlvLuZ9VslAj60H/DMs+Sx3ZO1nHc3syJUJuC3G/DMu+Qx7claftqUmRVl4HP4de0GPBdzyWLn3c2sDCoT8KF14PWSxWY26CoV8FtxyaOZDToH/AZOvZjZIKvMoK2ZWdU54JuZVYQDvplZRTjgm5lVhAO+mVlFKFLWmCkDSceBHy7yZc4DfrLI1+hF2dsH5W9j2dsH5W9j2dsH5W9jv9r3qxGR+LjAUgf8fpC0JyLGi25HmrK3D8rfxrK3D8rfxrK3D8rfxjK0zykdM7OKcMA3M6sIB3zYUXQD2ih7+6D8bSx7+6D8bSx7+6D8bSy8fZXP4ZuZVYV7+GZmFeGAb2ZWEQ74gKQ1kh6StF/SHknrim5TM0nvl3RE0iFJHy+6PWkk/Z6kkHRe0W1pJGm7pMOSvivpLkmjRbcJQNJVtf+uRyVtKbo9zSRdKOkBSY/W/u19oOg2JZE0JGmfpK8U3ZYkkkYl/WXt3+Cjkl5dRDsc8Od8HPhYRKwBPlp7XxqSrgCuBV4eES8F/qjgJiWSdCHwBuBHRbclwX3AyyLi5cDfAFsLbg+ShoBPA28CXgK8XdJLim3VAqeAD0XEPwUuB95bwjYCfAB4tOhGtPAp4GsRcQnwCgpqqwP+nABeUHt9DjBZYFuSvAfYFhH/DyAiniq4PWk+Cfx75v73LJWIuDciTtXePgRcUGR7atYBRyPisYh4Gvgiczf20oiIJyPiO7XXP2cuUJXqoRGSLgCuAT5bdFuSSHoB8M+AzwFExNMRMVVEWxzw59wIbJf0BHO958J7f01+HXiNpG9L+j+SfqPoBjWT9BZgIiIOFN2WDH4H+GrRjWAucD7R8P4YJQumjSStAtYC3y62JQv8CXMdjWeKbkiKXwOOA/+tlnb6rKSzi2hIZZ54JenrwPkJuz4CvA74YETcIemtzN2JX1+i9p0BnMvcV+rfAL4k6deizzW1bdr4YeCN/WxPs1bti4gv1475CHNpis/3s20plLCtdN+OACQ9D7gDuDEiflZ0e+okvRl4KiL2Snpt0e1JcQbwSuD9EfFtSZ8CtgD/od8NcR0+IOkEMBoRIUnAiYh4QbvP9YukrzGX0vnftfc/AC6PiOOFNqxG0qXAN4CTtU0XMJcWWxcR/7ewhjWR9E7g3cDrIuJku+P70J5XAzdFxIba+60AEXFzoQ1rImkY+AqwOyL+uOj2NJJ0M/AO5m7iz2UuNXtnRPzrQhvWQNL5wEMRsar2/jXAloi4pt9tcUpnziTwz2uvrwS+X2Bbkuxkrl1I+nXgTEq0KmBEHIyIF0bEqto/6mPAK0sW7K8Cfh94SxmCfc3DwMWSLpJ0JnADsKvgNs1T6wB9Dni0bMEeICK2RsQFtX93NwD3lynYA9T+f/CEpNW1Ta8DvldEWyqT0mnj3wGfknQG8EtgU8HtaXYbcJukR4CngXf2O50zAG4BngPcNxfDeCgi3l1kgyLilKT3AbuBIeC2iDhUZJsSrGeuB31Q0v7atg9HxD0Ftmkpej/w+dqN/THgXUU0wikdM7OKcErHzKwiHPDNzCrCAd/MrCIc8M3MKsIB38ysIhzwzcwqwgHfzKwi/j+jNYy73QHy4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are clearly three separate groups of points in the data, and we would like to recover them using clustering.\n",
    "Even if the groups are obvious in the data, it is hard to find them when the data lives in a high-dimensional space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now we will use one of the simplest clustering algorithms, K-means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "kmeans = KMeans(n_clusters=3, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** what would you expect the output to look like?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can get the cluster labels either by calling fit and then accessing the \n",
    "``labels_`` attribute of the K means estimator, or by calling ``fit_predict``.\n",
    "Either way, the result contains the ID of the cluster that each point is assigned to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "labels = kmeans.fit_predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 2, 0, 1, 0, 2, 0, 0, 2, 2, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2, 1,\n",
       "       1, 2, 2, 2, 0, 1, 1, 1, 1, 0, 0, 1, 2, 2, 2, 2, 0, 0, 1, 2, 0, 2,\n",
       "       2, 0, 1, 1, 1, 0, 0, 0, 2, 1, 1, 1, 2, 2, 0, 2, 1, 0, 1, 0, 1, 1,\n",
       "       0, 1, 0, 0, 0, 1, 1, 2, 0, 1, 0, 1, 0, 0, 2, 0, 2, 1, 2, 2, 2, 0,\n",
       "       2, 0, 0, 0, 2, 0, 2, 2, 2, 0, 1, 2], dtype=int32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all(labels == kmeans.labels_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's visualize the assignments that have been found"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5hTZfbA8e9JncYMQ+8CCipFUUcUsQso2HvZta0rYlvXXV3brmvZZl31Z8Wy6q69d4oioiLIoIA06SC9Mz2ZJOf3R8IwmSRTmECGyfk8zzyTvPfm3pMRT27e+77nFVXFGGNM8+dIdQDGGGN2D0v4xhiTJizhG2NMmrCEb4wxacISvjHGpAlXqgOoTZs2bbR79+6pDsMYY/YY06dP36iqbeNta9IJv3v37hQWFqY6DGOM2WOIyPJE26xLxxhj0oQlfGOMSROW8I0xJk1YwjfGmDRhCd+knKqiwY1oqCTVoRjTrDXpUTqm+VP/NHTbrRBcByjqOQJpeR/iaJXq0IxpduwK36SMBn5Bt1wJwV8AP1AJ/m/Rzb/Bqrgak3yW8E3KaNkroJU1WgMQXAqBOYlfF9qCBpagMa81xtTGEr5JneBSIF7SdkJwVUyrhkoIbbkWXX8UuuksdP3hhMreid6n/GNCG08mtG4goc0j0cqfd03sxuyBLOGb1HEXABmx7VoJrj6xzVv/AL6vAD9oGWgxFN2N+qYAECp9Ht12BwQWgm4F/1fo5vPRwKLYY6naNwSTdizhm5SRrPPBkUP02IEMyDgRcXWN2leDG8A/mXBff3UVaOloVP1Q8n9AefVXgVagxY/taFElVPoquuEIdF0/QuuPJFT2bnLfmDFNlCV8kzLiyEVavw+ZZ4KjNTi7QYsbkbz7YncObQBxxz9QcE3cLqDIC6FyVtUzLXsdiu+D0CZAIbQeiu4mVP5xo9+PMU2dDcs0KSXOdkje3+ve0dUDCMXbAJ7DwNEWNBD/tc4uOx6XPkb0twDCz0segcxT6he0MXsou8I3TYb6ZxLadD6htf0IrR9MqOR5VMNJXiQTsm8AvNVe4QTJQnKuQhw5kHkasfcEMpCca8LH12Dkyj6O4Opkvx1jmpx6J3wReUFE1ovI7GptrURkvIgsjPzOT/Dak0TkZxFZJCK3JiNw07xo5QJ08yVQ+SPgD3fhlDyGFt8f3q7l4BsPKCCRn0zIfxFxdgwfxD04sn07D7T4K+I9AgARJzjaxw9AWhHaMorQ1htR33e75D0ak2oNucJ/ETipRtutwBeq2gv4IvI8iog4gSeA4UAf4EIRiR2CYdKaljwO+Gq0lkPZK2ioBC1+FCpnE75pq5GfCih9PPx6/wwouq3GMRR8Y6MPmXMTsd8CJDyqxzcBKj5Bt4wiVPxI0t6bMU1FvRO+qk4CNtdoPh14KfL4JeCMOC8dCCxS1SWq6gdej7zOmB0C84jbRy8uCK6E8neJ/UAIgG8Sqn60dHSc7ZXg/w4Nrq1qcWSdDnn/AudegBscHQAP0aN/yqH0uajXbafBjYS2/I7Q2r7hny03oMEE3UTGNDGN7cNvr6prACK/28XZpzPwS7XnKyNtcYnISBEpFJHCDRs2NDI805SpBlDfV2jZ2+DoRLibpuZOAXB2AhLckEWBUKQ8Q5xyDOKO1OnZwZE5Akfb8Tg6zAHvUcR+UBD+oPFNrhGvH918Hvg+JzxhrBJ849HN59qYfrNH2B03beP8Xxzv/8zIBtXRqlqgqgVt28ZdltE0AxpYgW44Dt36e7ToXqicTuw/FS9knoE4csE7hNhBZQLu/ohkgPvQONuJTOLaO3Egjrz4r0MicwSq8U2A0BaiP3wC4TbfhMTnMKaJaGzCXyciHQEiv9fH2WclUH0WTRfAhkSkOd16XfjGrJYSHibpJ7ZLJwiZ5wMgLW4GRxsgM7ItA6QFkveP8PacK0Eyif4nnQnZvw2P4ElAMs8ifsJ3gPeY6KbA4vAM35g3UxbeZkwT19iE/yFwaeTxpcAHcfaZBvQSkR4i4gEuiLzONFPh4ZXnElrbh9C6wwhtu5dQ+Udo5bzw9uAqCCwl/rj66gJQfC8A4myLtBmD5P45/CHQ4g9I2y8Q1z6R7R2R1u9BxsnhMfmu3kjePUjO72o9g7j2gdx7CH+A5IBkg+Qj+S8g4o3e2dWTHR841Q+SVfu3CGOaCKlvGVoReQ04FmgDrAP+CrwPvAl0A1YA56rqZhHpBDynqiMirx0BPAI4gRdUtR4zbaCgoEALCwsb9IZM6mhwE1o5H7ZeDVTU2OoK/7j7Q4vbYfNFxE6AikeQ9vMRidczmDwaKoHKQiADPAWIxF71q/rRjSdG7gls79ZxgbMD0mYskmgmsDG7kYhMV9WCuNuact1xS/h7Bg2uRbfeCJU/EU6EtV25eyDjDPBPDJc1qFMmjg4zkxNnqCQ8/LPiY0DC9wdyrglP6qrvMYIb0aJ7wPdFuMF7ApJ7J+Jsk5QYjWms2hK+lVYwjaIaQjdfHB46SbAer/BDxQeQ/yxsvQo0SGxBtO28kHl2kuIMoJsvgMCyHecrfRH1T4FWb9b7G4Q42yD5j9W9ozFNkJVWMI1TWQihjdQv2Ve9CPEchrQZC1mXEe7pi8N9CJL7p8bHCOD7MvKhVP3DxRcupeyfkpxzGNPEWcI3jRNnclLtBNwFiAji7IC4OhGe+FSTB8k8KTzkMgm08qcEI2z8kRm8xjR/lvBN47j7J6hSWb2LZPsVvBckB8m7a8em0DZib/AC+NGKcTGtqr6qgmoNIc4uxB9h4wVnwnmAO84bXE1o218IbRhCaNNFaMWXDY7BmFSzPnzTKOLqgWYMhYov2DHqxg3O9tD6QwR/uAZ95ZxwiePAGnTTuSgankzl7JX44P5paMXnSMYQ1DcpfLM0uBIkA828CGnxh5jRNBoqC+/jbI848nZsyBgBxQ+AVrBj3p8jPHY/Y0it71GDa9CNp4OWAEEIrkC3zkFb/BFH9iUN/IsZkzqW8E2jSd4DqPs1KHslnFAzhiE5V1dNeJKca1ANoptOi9w0jZQhqPiE8IieRCPFfGjpc+BojW65jqpvAloGZf9DtQTJuyfcpIqWPg4lz4I4QSvRjJORvHsR8YRjaf0auvVPEIisc+vui+Q9SHh6SGJa8nRkglj1+xTlUPIwmnV+7Hh9Y5ooS/imXlQrwPcNqA+8RyCOHZWwRZxI9q8h+9eJD+D/JlJzvnrNmXrc6A2uQ0ueILbbpwLK30Vb3Iw4WqDlb0PJc+H27Z8fFZ+hkoXk/TUcp2sfpM27aGgrINHfAGrjn0L8Wj4CgSXg3r9+xzEmxSzhmzqpbyq69WqqMqkG0BZ34Mi+oP4HqVwQ/rBosEoILo2/SdwQWgeOFlA6mtiJXBVQ/jaae1vUVbw4WjYsBGeH+DFoZaTcgzF7Bkv4plYaKkO3XhU7wqX474S0BHwTAR9knIZknZ+4e8TVPXyDNNEyhImENoPnqMiatTVu1mowUmUzsl9cfjS0BXEmWPikHiR7ZLjeftQHigc8RyBOK/Bn9hw2SsfUzjeR+AVPfVDyMFR+D5UzofgBdPPFaKKE7j0OHPlEj7l3RH5qu+5wQvYlRC9tSPhma/ZliCMr/Nx9YILXK5S9Vcvx6ybewZB7e6TWThbgAe+RSMuHG3VcY3Y3S/imdlpO4lIJ1ZN7RfhmaIIywSIupNUb4D2WcNJ3gucIaP0uZF0USaRxONsiniOQ1v8FdwHgBUdHyLkZyblxx35Zlyd+DxXvJ95WT46s85F2U5DWbyPtJuHIf7rWKpzGNEXWpWNq5x0MRfUc965lqO9bJGNY3M3ibIvkPxVeTJzIGrMA7j5ozjXoxjMhtJVw14kbxIXk3Rcue+A+AGn9avg0oW1QOQsCc1FXn/AkLncvNGblqu2Sc10j4oFIdU5j9kSW8E2txNkBzbkOSp5gR816D+ERNjVH2bjDpYnrOqbEllIQRyto8yla/h74vwfXXkjWBUiNSVGhkmeg5PHwDVtC4av9Vs8jzk6oq3u4VELUMM+MpNXjMWZPZ9UyTb2ofyZa/k64i8c7FIr+Arqlxl4ZSNsx4eRbOQ8tfgQCs8DREcm5Dsk4vnEx+CahW64n+uapA1y9cLT5CA0sQjf9CvCHRwSJB1x9kVb/ibqZrFoB5R+j/sng7IJknoe4ujQqNmOaCquWaRpNPAcinh03RtW1F7rlKtDwmHZwIS0f2pHsN11AeOy8QmhTeCnD3D/jyDpvp2PQ0peJHXoZgsByNLAovJhJu6+gYly4Zr3nQHAfGlUJU0PF6KazIzXtw11HWvoS5D+NeAftdGzG7Aks4ZudIu59oe2XEJgfLkDm7ltV5kCLH6Yq2VepgKJ7CVV8Dq4eSNavEFe3hp00VPMbxQ4a3IK4CBdbyzwt8X6lz0YmgG3v6w8vRq7bboa2kxCxcQym+bJ/3WanhW+W7h+++q9e06byJ+KXS/CFFz4p+x+66VTUN7VhJ8w4gZjhmUD4w+TO8KpVdakYQ9wbu6FiCC5vWDzG7GEanfBFZF8RmVHtp0hEfl9jn2NFZFu1fe5s7HlNE+bsUMcOlaDl6LZbaMg9JMm6JFyULd68gOAv4dWs6jxIonLLocgi6MY0X41O+Kr6s6oOUNUBwCFAGfBenF2/3r6fqt7T2POapktyrgXqUcc+tClcGqG+x3XkQP6LxJ8I5o8sXRi+KauJyjhk/prYMsnhG79S5weVMXu2ZHfpnAAsVlX7bpzGJGNoeKFyySWc+BMtH9jwq2pxtCDhClkaIrTp1+i6g9B1AwhtvhwNrol+fdY5kHEi4dr8WSDZ4OyE5P9fg+IwZk+U1GGZIvIC8IOqPl6j/VjgHWAlsBq4SVXnJDjGSGAkQLdu3Q5Zvtw+O/ZUqgEIbUDLP4OSR4keYeMCz2E4Wv2nwccNbTo/XM4hagbw9g+B6uWWneBoi7T9PKbGjwaWhY/haAeew+xmrWk2ahuWmbSEL+H/o1YDfVV1XY1tuUBIVUtEZATwqKrWsvJFmI3Dbx5UQ2jRbVD+6Y4JU86uSP6LiLN1w48X+CW8ILmWRRY0qaWmvmQjef9CMk5szFswZo+xu8bhDyd8dR/TKauqRdUefyoiT4pIG1XdmMTzmyZKxIHk3YdmXw+BOeDsCK7+UePjG3Q8V9fwkFDfBLToAQj9knhn9UFgxU5GbkzzksyEfyHwWrwNItIBWKeqKiIDCd872JTEc5s9gLi6QJJmtIp4UPcACK2uY0cPuPdNyjmN2dMlJeGLSBYwFLiqWtsoAFV9GjgHuFpEAoQ7ci/QplzTwewZgr8QvnZItHKWE5xdwXPkbgzKmKYrKQlfVcuA1jXanq72+HGgHoOkjWkAZ8/at2ecieTebjdkjYmw/xPMHkucrSHzLOJet7T4M46W/7Ca9cZUY7V0zB5Ncu9CnV2h9AXQInB2g9y7cXgPS3VoxjQ5lvDNHk3EieSMhJyRqQ7FmCbPunSMMSZNWMI3xpg0YQnfGGPShCV8Y4xJE5bwjTEmTVjCN8aYNGEJ3xhj0oQlfGOMSROW8I0xJk1YwjfGmDRhCd8YY9KEJXxjjEkTlvCNMSZNJCXhi8gyEflJRGaISMyq4xL2mIgsEpFZInJwMs5rjDGm/pJZHvm4WhYlHw70ivwcBjwV+W2MMWY32V1dOqcDL2vYFKCliHTcTec2xhhD8hK+AuNEZLqIxFuJojPwS7XnKyNtxhhjdpNkdekMVtXVItIOGC8i81V1UrXtEuc1Gu9AkQ+MkQDdunVLUnjGGGOScoWvqqsjv9cD7wEDa+yyEuha7XkXYHWCY41W1QJVLWjbtm0ywjPGmBiqytLZK1j603JU415/NjuNvsIXkWzAoarFkcfDgHtq7PYhcJ2IvE74Zu02VV3T2HMbY8zOmP/9Qu4+5yFKtpQAkJOfw1/f/iP7DeyV4sh2rWR06bQH3hOR7cd7VVXHiMgoAFV9GvgUGAEsAsqAy5NwXmOMabDSbaXcMuxeyorKq9oqSn3cMvReXl3xFNl52SmMbtdqdMJX1SXAgXHan672WIFrG3suY4xprK/emkIoGIppD4VCfPXmd4y4ckjMtkU/LmXe1IW07pTPwOEH4XInc0T77rNnRm2MMTtpy7qt+Mv9Me2+cj9b1m2LagsGgtxz7kNMHz8TFBwuJ5k5Gfx70j102rvD7go5aay0gjEmrfQ/an88mZ6Ydm+mh35H7RfV9uGTY5g+fia+Mj++cj/lxeVsWbeVe89/eHeFm1SW8I0xaaX/UfvTb/B+eLO8VW3eLC/9Bu/HAUf3idr3k9Gf4yuL/jagIWX53JVsXLVpt8SbTNalY4xJKyLCvR/dymfPT2DMCxMAOOny4xj+2xOIDD6pUukLxD2GQ4RKf/xtTZklfGNM2nG5XZw6ahinjhpW637HXTiYNx/4kEpfZVR7fvs8OnRvtytD3CWsS8cYYxI47+bT6dyrAxk5GQB4Mtxk5mRw2ys3xHwb2BPYFb4xpskLhUJ8/Mx4Pvi/zygrLmfQaQVcfOe55LdvuUvPm9UikycL7+Pb975n1qS5tO/elmGXHLvLz7urSFOeUlxQUKCFhTHl9Y0xaebBK55k4huT8ZX5AHC5neS1zeX5Of9u1hOldoaITFfVgnjbrEvHGNOkrVu+gQmvflOV7AEClUFKtpTy2fMTUhjZnscSvjGmSVswfQlub2zvs6/cz4yJc1IQ0Z7LEr4xpklr17V13FIITreTzvvsebNdU8kSvjGmSetdsDed9u6A0+2Mane5XZx2zYkpimrPZAnfGNOkiQj/GvcXDjymL26PC0+mh3bd2nDvh7fQeZ+dXyl15cI13H3Og5zV+jIu7XU9Hz41ttnXxbdROsaYWq1cuIbZX8+jZbs8Ck48MKWVIos2F1NR6qNtl9aNGge//peNjDzgj5QVl6OhcA70Znk5ZeQQRj18WZKiTY3aRunYOHxjTFyqyiOjnuHz/07C4XQgDgfeTDcPTriLvfp0rfsAu0BuqxbktmrR6OO89cCH+Mp8VckewFfm48OnxnHRHWeT23rnzhEMBpn01hQmvPo13kwPJ11xAocMPaDJTNKyhG+MiWviG5OZ8Oo3+Ct2lBWoKCnnztPv58UFjzWZJLYz5kyeT6AyGNPuyXCzYv4q+g3eL86rahcKhbjz9PuY9dVcKkrDQ0infvoDp4waxlUPXNLomJPB+vCNMXF99PTYqsS1nSpsXruF5XNXpiiqxFSVJbOWs2D6YoLB2GReXZfenRBH7AdWpa+Sdt3a7NT5p4+byaxJ86L+ZhWlPj58YgxrlqzbqWMmW6MTvoh0FZEvRWSeiMwRkRvi7HOsiGwTkRmRnzsbe15jzK5VfaJTdQ6HA1+cBUR2ha/fncr1g27n4p7X8sioZ9iwMn5J4kUzlvLrHtdww5F3cNPxd3F+p5H8OOGnqH0q/ZVsXLWJSn8l5/3pdDwZ0TXxPRluDjqhP+267lzCn/rpD1SUVMS0i8PBD5/P2qljJlsyunQCwB9V9QcRaQFMF5Hxqjq3xn5fq+opSTifMWY3OP7CI1k+Z2VMcne6newzoPsuP//r973HK/e+Q0Xkg2fMC18y6Z0pjJ75EG06taraz1fu4+YT7qZkS2lVW3lxBXeedh8vLfo/Ksp8PHj5k8z57mccIrg8Ls696TTufPuPPHr1aDav2YoIHH3eEfzuid82KEZVZc7kn9m8ZgsigtPtJFijq8jhFLLzshrxl0ieZKxpuwZYE3lcLCLzgM5AzYRvjNmDnDJqGBNe+5YV81dRUVKBy+PE6XJyy8vX43Q56z5AI5QVl/O/e96O+rAJBoKUF5fz1oMfcnW1kTSTPygkGIjtwgmFQrz01zcZ858JhALhiVshwmUZ3rj/Ay69+zz+t+RJijYVk5mTASJMfP1bJn84jVYdWnLKVcPoecBeCWPcsHITfxpyN5tWb0Ecgr+iMu6wThHh8FPjDprZ7ZJ601ZEugMHAVPjbB4kIjOB1cBNqhp3TrSIjARGAnTr1i2Z4RljGsCb6eXRb//Gt+99T+HYGbTu3IqTfnP8bqkDv2LeyvBEq/Lo9oA/yIwJs6Patq7fFnNVDeCvqOSLVyZVJfuobeV+3rjvA8676XTy2uTiK/dx45F/5pefV1NR6sPhdDDupYnc8NRIhl58TNwY7zn3IVYvXhc1C9jlceLwuPBkuAFwupzc++GtZFRbXSuRQGWAz56fwLiXJuJwCsOvGMLQS47G6Uzeh2vSEr6I5ADvAL9X1aIam38A9lLVEhEZAbwP9Ip3HFUdDYyG8Dj8ZMVnjGk4l9vFMecdwTHnHZHU465atIb3/+8zfvl5NQcc04dTrhoaNdyyVcd8AglWlGq/V9uo5wcc0yfuDVhvlpfKisqY9u2KNhUTCoVwOByM+c+XrJi/uuq+RSgYwlfm57FrnuWosw+PSdgbVm5i8cxlMSUfAv4gew/ozm/+fhFur4v+R+1fr3kLqsrtJ/+DuZMXVMWwZOZypnw0jb++c3PSRkQlZZSOiLgJJ/tXVPXdmttVtUhVSyKPPwXcIrJzd0aMMU1eMBBk/H+/4rbhf+PO0+/ju48Kq7o7Zk6cw1UDbuajp8YxfdxMXrn3ba7oc2PUGrHturah35H74/ZEJ0tvlpfzbj4tqm3vA7sz+IyBZGRXX6PWQ4/+3Qhp7NX9dvnt83A4wilw0pvfxb1J7XA6mD91YUx7eUkFTmf89FlR6mPg8IM46Pj+9Z6k9uOE2cybsjAqhopSH9PHz2L+94vqdYz6SMYoHQGeB+apatyl3EWkQ2Q/RGRg5Lx73grAxpg6hUIh7jjlnzx2zbMUjp3Jdx8V8o+LHuGJG15AVXnwiifxlfmq+t39FZUUby7mpb++EXWcO9/6AwcPPQC3101GTgY5+dn8/pmR9Dty/5hz3vLy9fzuiSvpd+R+7DdwH668/2Ie/urumG8D1V14+1lVjxNVHNCQktkiM6a9c68OZGRnxLS7vS6OOvuwhOdMZObEOXFH+FT6A/w0KXm3Q5PRpTMYuBj4SURmRNpuB7oBqOrTwDnA1SISINwrd4E25ZoOxpidVjh2JnMm/xwzHv2z5ydw3AVHsmn1lpjXBAMhpnzyQ1Rbdl42f/voNrZtLKJoUzEde7ZPeMXscDgYeskxDL0kur/9H5/ewZX9/xDT9dKpVwfWLFnHf/78GoNOP5QF0xfHPW5Oq2x6H9Izpt3pdPKnF6/l7nMeIuAPEAwEycjy0rJ9HufdfHr8P0wt8tvn4cn04K8xIsrjddOyXV6Dj5dIMkbpfAPU2sGkqo8Djzf2XMaYpu/7z36Me7UKMG/qwoRX01lxrqQB8trkktcmt+p5MBjkm3emMu7lr9i8Zgs5+dkcdHw/hl9xQszSg93268z/lj7Bg795ijnfzsed4aZ1p1asXbqe9x79FIBX/xHTC11lxG+HJOw/P/Skg3j6h/v56OlxrF26noOHHsCwS44hMyf++6jNcRcM5vnbX41pdzgcO/WNIRErrWCMSaq8Ni1weVwxN12dLgdtOuVTMOxACsfOiCpt4HA68GZ6mPDaNxxz7qCEwz5DoRB/PfMBZkz4CV/ZjqvhWV/N4Y37P+Dfk+6NGUrZtksb7hv3FwAKx83k7rMfTDiprDqHQ+q8Wdqld6eoIaI7K69NLv/89HbuOfdhKkorUCCnZTZ3vXvzTn2AJGLVMo0xSbV22Xp+2/fGmAlb2XlZvL5qNP4KP7ed9HeWz11JKBii0rdjJE1Gtpf9B/Xmn5/dETUcMVAZYPXidSyZtYwHLnsiqr5Pdb0L9uaJ7/+VMLaHrnyKMfVcFjEj28vfP7mdA47uU6/961JeWsG4lyZSOGYG7bq14dSrT6R73+gidKFQiCUzlyMOoecBe+3U6ByrlmmM2W06dG/Hrf/7Hfdf9ng4YSm4M9xV49Ezsrw88f2/+OHzWdw+4h9Rr60o9TFvykK++7CQI88Md2WM/+9XPHnDf8ITr0oroJZr1MUzllJeWkFmnBuqEC6fIA6JqpIZjzfLQ9/B+9H/qNgbxDujdFsp1xTcyuY1W6goC4/zH/vil9z2vxsYfMbAqv0cDgf7HNQjKeeMxxK+MSbpjjzzMAYOP4g5k3/G7XGx/6DeMROItq7fhifDTXlJ9KSpipIKvn3/e4488zBmTZrLo1ePjuq+qY2I4HInnqg09OJjGPufL2s9Xvf+3Tj1qmGMuPKEpI1/f/vfH7Nh5aaqbzPbx/k/9NunOPyUQ3b5zOXtrFqmMWaX8GR4OOj4/vQ7cv+4s0UzW2TGTagOp4MW+dkAvPnAB/VO9k6Xg8NPLcDtcVNeUs6bD37AtYfdyi3D7mXyB9NQVfYb2IsLbzsr4THcXhdX/uvXnHbNiUld6OWbd6ZGdV1tF/AHWDbnl6Sdpy52hW+MSYlDhh2I0xV7zen2uBh+xQkArF+xsd7H26tPV2585ip85T6uH3QHa5esq7qPMPe7nznjuuFc8c9f8as7zmbpT8v56s3vYo4hDgetO+XXep5KfyUTX5/Mt+9PJa9tHqdcNZReB8cO3awuK0HxtGAwRFZu8m7K1sWu8I0xKeHxuvnnmD+T26YFWS0yycrNxJPp4ZpHL6dH//BImwHH96u1iwbA5XZy0R1n8/SPD5DbugVfvPIN65atj7ppXFHq451HPmHTmvAcgAtvOwtvVnR5ZKfLSdd9O7H3gd0Tnsvvq+TGo/7CY9c+y7fvT2PM819w49F/4dPnPq81xjOvHx41ExjC32T22r8zHXu0r/W1yWQJ3xiTMvseug9vrBrNX9+9mVtevp43V49mxG+HVG0/7+bTyczNjOrjdrocOF0OMnMyyMj2cuk953PZPedXdQ9N/WR6zMItEP7mMPe7BUC4HMOt//0dua1bkJGTgdvrrhodVJvxL01k+dyVVccPhRRfmZ8nb/gP5SXlCV93zHlHcPJVQ/FkuMnKzSQzJ4OOPdtx17s31/+PlQTWpWOMSSmX28XBJ/SPu61Np1Y88+ODvPK3d5g+fiatOrTkvJtPZ8BxfdmybhvturWJWcikdad8HE5HzOxaRfFkuHng8ieY+Ma3BAMhDhl6AGf/4RS67d8lqsZ+IpPe/i7uh4nT7WTO5IzIF8YAABxXSURBVAUUDDsw7utEhFEPXso5fziV+VMX0qpDS/Y/vPduXybSEr4xpklr26U1v396ZEx7dl523P1PHTWMcS9NjLrZKw4ht1ULXrjjVVbMW1U1Kaxw3EwWTF/MSwvrVwggp2X8c/rL/eGa+nVo06lV1XDTVLAuHWNMs9Kj/17c9Pw1VfcFvFleuvTuxOV/u5A1i9dFzQAOBUNUlPr44n+T6nXsU0YNi3tVHggE2bwmtkbQdsvnreSZm1/mwd88wTfvTa1zzd1dxRK+MaZZWblgNWNfnIjfV0koGOLocw7jqen3UVZcHjfRVpT6WDJreb2O3bJdHo44I4tQ+PiZ8XFf8/n/vuLaglt479FPGPviRO679HFuHfa3uKt07WqW8I0xzcaW9du4ftDtTB83k4A/QEWpj6/e/I67znqAbvt1jjsfICPbW+tShtVVlPrwZnribistKotpKy8p55FRo/GV+wlGVt6qKKlg/vcL+fL1bxvwzpLDEr4xptn45Jlx+Mv9URU5/RWVzJo0j7y2uXTap0PUoioOp4OMbC8n/Proeh1/n4O6x233Zno4Ns6qYLO/mR93Fm1Fqc8SvjHGNMbPhYvjFlZzuZ2smLeSByfcxfEXHYkn04PT5eTQkwbw+NR/JSzNXJPb4+am56/Bm+WpmjSWke2lc++OnDJqWMz+ngTfBgAyc+pe5zbZbJSOMabZ2HtAD6aPm0mlL7o0czAQpOu+nchpmc1NL1zLTS9cu9PnOOrsw9mrb1c+GT2eTau3MHD4QRx7wWA8XnfMvv0G74fH66asxmrsGdleRlw5dKdj2FlWHtkY02xsWrOF3+x/A2VFOxKs2+umz+G9ePDLu6P29ZX7mPDqN8yaNJdOe3dg+G9PqNdY/Ib6uXAxt514L4FACFUlWBnk3JtP4/J7Lkj6uaD28shJSfgichLwKOAEnlPVf9XYLpHtI4Ay4DJV/SHmQDVYwjfGNNTS2St47Jpnqyp1Drn4aEY9dGnUQiJFm4u5buBtbFm3lYpSH26vC5fbxb/G/pk+g/ZNekx+XyWFY2dQuq2MAcf1o22X1kk/x3a7NOGLiBNYAAwFVgLTgAtVdW61fUYA1xNO+IcBj6pqnbMPLOEbY3ZWKBRCJP6qVU/94UU+enIslTVW5eq4d3teWvB/u30GbDLVlvCTcdN2ILBIVZeoqh94Hai5iu/pwMsaNgVoKSIdk3BuY4yJy+FwJEzcX78zJSbZA2xatZmNqzbv6tBSJhkJvzNQvaDzykhbQ/cBQERGikihiBRu2LAhCeEZY0y0mvV3tlMN19tprpKR8ON9hNbsJ6rPPuFG1dGqWqCqBW3btm10cMYYU9MpVw2NKY/scDrY99B9yGuTm6Kodr1kJPyVQPWVeLsAq3diH2OM2S3O/N0IDj3pIDyZHjKyvWS2yKBdtzbc/urvUx3aLpWMcfjTgF4i0gNYBVwAXFRjnw+B60TkdcI3bbep6poknNsYYxrM6XLy17dvYtmcX/h52iLadm3DgOP64nA077mojU74qhoQkeuAsYSHZb6gqnNEZFRk+9PAp4RH6CwiPCzz8sae1xhjGqt7365079u17h2biaTMtFXVTwkn9eptT1d7rMDOT20zxhjTaM37+4sxxpgqlvCNMSZNWMI3xpg0YQnfGGPShCV8Y4xJE5bwjTEmTVjCN8aYNGEJ3xhj0oQlfGOMSROW8I0xJk1YwjfGmDRhCd8YY9KEJXxjjEkTlvCNMSZNWMI3xpg0YQnfGGPSRKMWQBGRB4BTAT+wGLhcVbfG2W8ZUAwEgYCqFjTmvMYYYxqusVf444F+qnoAsAC4rZZ9j1PVAZbsjTEmNRqV8FV1nKoGIk+nAF0aH5IxxphdIZl9+L8BPkuwTYFxIjJdREbWdhARGSkihSJSuGHDhiSGZ4wx6a3OPnwR+RzoEGfTHar6QWSfO4AA8EqCwwxW1dUi0g4YLyLzVXVSvB1VdTQwGqCgoEDr8R6MMcbUQ50JX1WH1LZdRC4FTgFOUNW4CVpVV0d+rxeR94CBQNyEb4wxZtdoVJeOiJwE3AKcpqplCfbJFpEW2x8Dw4DZjTmvMcaYhmtsH/7jQAvC3TQzRORpABHpJCKfRvZpD3wjIjOB74FPVHVMI89rjDGmgRo1Dl9V90nQvhoYEXm8BDiwMecxxhjTeDbT1hhj0oQlfGOMSROW8I0xJk1YwjfGmDRhCd8YY9KEJXxjjEkTlvCNMSZNWMJvpCJfBetKSkhQVcIYY5qMRk28Smeby8v4w7jP+O6XFThEaJudzQNDTuKwLl1THZoxxsRlV/g7QVW55L23mfzLCipDIXzBICuLivjNh++yfGvMgl/GGNMkWMLfCT+tX8eyrVsJhEJR7YFQiP/O+jFFURljTO0s4e+E1cXFOBwS014ZCrF065YURGSMMXWzhL8T+rVrR2UwGNOe4XJxeGfrwzfGNE2W8HdCl9w8Tu61L5muHfe8XQ4HuR4v5/frn8LIjDEmMRuls5PuG3Ii/dq15+VZP1Lqr2RIz735/WFHkOvNSHVoxhgTlyX8neR0OLhswMFcNuDgVIdijDH1Ygm/CZq2eiWvzJrJVl8Fw/fpzRn77o/XZf+pjDGN06gsIiJ3AVcCGyJNt6vqp3H2Owl4FHACz6nqvxpz3ubs2enTeGTqZCoCARSYtmoVr82exRtnn29J3xjTKMm4aftvVR0Q+YmX7J3AE8BwoA9woYj0ScJ5m50t5eU8POVbyiPJHqA8UMnCTZv4aMH8lMZmjNnz7Y5ROgOBRaq6RFX9wOvA6bvhvHucwtWrcDudMe3lgUrGLl6UgoiMMc1JMhL+dSIyS0ReEJH8ONs7A79Ue74y0haXiIwUkUIRKdywYUOi3Zq0kCori7axtaK8Qa/L9XqJV4PNIUJ+ho3+McY0Tp2dwiLyOdAhzqY7gKeAewGN/H4I+E3NQ8R5bcLSkqo6GhgNUFBQsMeVoJywdAm3fTGOYr+PkCqDunTj3ycOp2VGZp2vLejUmWyPm9JKf1S7x+nkV/0P3FUhG2PSRJ0JX1WH1OdAIvIs8HGcTSuB6tNPuwCr6xXdHmb+xg1c99lHVAQCVW2Tf1nOhe+8yYAOHSnx+xjRa1+G9dwHpyP2y5XT4eDlM87h0vffocTvRwQqgyFuG3w0B3bouDvfijGmGWrsKJ2Oqrom8vRMYHac3aYBvUSkB7AKuAC4qDHnbape+HF6TMmFylCInzdtZMGmjSjw5bKlvN6xEy+cdlbcpN+7dRu+/c1Ipq9ZRYnfzyEdO5Pr9e6md2CMac4aO87vfhEZQLiLZhlwFYCIdCI8/HKEqgZE5DpgLOFhmS+o6pxGnrdJWr5tK8EEC6Fsby2rrGT6mtV8vnQxJ+7dK+6+DhEO7dQl9hiq+IIBvE4XIvF6yowxJrFGJXxVvThB+2pgRLXnnwIxQzabm0FdujFz3Vr8cQqrVVdWWcmYRQsTJvyaVJUXZvzAE99Pocjvo3VmFjcfcSTn9OmXjLCNMWnCiqfthGCNOvjbXXLgAFp4vLjquPp2itDC46n3+V6Y8QMPf/cNW30VhFTZUFbKXyd+wcc2Nt8Y0wCW8IGtFeVMWLqE6WtWEaplbdrXZs/isOeeotfj/2bQ80/zztzoWxatMrP4+MKLObdvfzrk5NC7VWsy4s2OVeXM/eo390xVeWLaFMqr3QgGKA8E+PeUyfU6hjHGgNXS4Znp3/PIlMm4nU5UlZYZmbx85jn0aBk9peD12bP426QvqxLvutJS/jLxC9aUllDqCw+jPLn3vvRr156/Hz+06nUz1q7h3DdfpXonTxD485ef8/75v4o70ao6fzBIkc8Xd9vq4qKGv2FjTNpKq4Svqrw+exaPT5vChrIyOua0YF1JMf7IurQQ7l+/7P13mHjpFVE3Rh+ZOjnmKrsiEODf332LRKYavDTrR648uIAbDx9ctc+6khLi9ejP37iBTxf+zOl1XOl7nE5aZ2axoaw0ZlvNDyVjjKlNWnXpPPdDIX/7eiJrSkoIhEL8UrQNf43+eAU2lZfx0/p1VW0hVdaXxibc7fuHUEIoFYEAo6cXsmjzpqrtz/wwLeHrHp7ybdyVs6oTEf50xJFRi61AeHWtWwYfXetra6oIVLKhrLTWbitjTPOVNgk/EArxeJy+8HgcIhT7fVHPO+a0qNd5ghri8yWLq56XV1Ym3HdVURE3jf+s1uOtKS6mpNLP8H160zU3D6/TxX5t2vL0yadzTPce9YrJFwhw6+djGfDMExz1n2cZ9PwzjFm0oF6vNcY0H2nTpbO1oqLO4ZLbBUIhBrSPntn6p8FHcdsX46Jm0cYTDIV4Y84sOuS0oHfr1qwrKUm4bwgYs2ghH/08n6F7702Gyx21/Y3ZP3HXV18A4W8EAtx4+GBGHnJovd5HVeyfj2Hc4kVV739DWSl/HPcZbbKyKeiUsKyRMaaZaZZX+OtKSthQowsmz+vF5aj9BimE+8z/fNSxuJ1O3p8/lxvGfMLfv55Iv7btuH/IieyV1xKXw0G33DzccWbKKrB82zbu+GIcZ77xClt9FbWerzIU4tYvxlLw7FN8tvDnqva1JcXc9dUX+IJBfMEg/sjvR6ZMZsmWzfX7QwCby8sYu3hR1T2K7coDAZ6YNrXexzHG7Pma1RX+/I0b+P2YT1i+bSsK7NOqNY+ddDI981vhdjoZVXAoT06bWmu3TobLxYcL5vHcD4WsLy2lLFCJS4RXfprJQ0NP4stLr6ja9605P3HnxPAVeExCDdbddVS1bySeP44fw/5t29G9ZT7jFi8iXt25ylCQTxcu4LqBh9fr2OtKSvA4nXG/3azYtrXeMRpj9nzN5gq/2OfjgnfeYMHmTVVXxPM2rOe8t1/HF0mo1xQcxo2HDyarRtdJdUU+H9+vWsWybVspC4T73wMaviF7yxfjohLnuX37M+nyKxnUpVtS3kMgFOLtyNh+RUlUVLQhN133aplPIM5EMYcIB3e0gmzGpJNmk/A/XvhzzIgXJXzDctyS8OIhIsJvDy7g+yuvpmtuLp56dPHU9NP6tVHP22ZlM7jbXnjrGE9fH4FQiC2RGvpDeu4Tdx+308mJ+9SvJANAltvNqEMGRo3yESDT5eLaQ+v3LcEY0zw0m4S/urgobleNLxhkTXFxVFuW281HF17MlQcXsHd+K1xx+uLjCanG/XZw2r774YhTTsEpErfMQqLzZbndHNe9JwCdW+Ryy+Cj8TpduB0OnCJ4nS6uOvhQ9m3dpl7xbnf9wMO559gh7JPfipYZGZzQc2/ePe9XdLdx/MaklWbThz+gfUey3G7KagyD9DidHNA+dv2WXG8GfzziSP54xJGc8urLzN1Y++paArTJymK/Nm1jtrXNyubZU8/k+s8+ojIYQlGy3R6eGHEqc9av4/FpU9lUXlY1qzYQjF+Lp0fL/KqED3DZgIM5tnsPPlu0gJAqw3r2olfr1nX9KWJjF+HsPn05u0/fBr/WGNN8NJuEf2z3HvRomc+iSB8+gNfppG/bdhzWObbUcHXXHHoYN48fE/UNwQmoSKQrRMh2u3n+1DMTliU+oms3pv72amavX4fT4aBv23Y4RDikU2cuGXAwpX4/K7Zt5ce1a/jnN19RWuODSYAD2rePqZHfvWU+Vxcc1uC/hzHG1NRsEr7T4eCNcy7g6cLveX/+XESEc/v05cqDD62zdvyIXvuyuriYR6ZORgj3pZ+0T29uGDiImevX0iojk0Fdu9XZ9eNyOBiQYGWqbI+H/du2Y2G1Wbgx70GaTQ+bMaYJEm3C0+wLCgq0sLBwt53PFwiwYts22mZn1WsN2p1R5Kvg8OefiZnAlely8dIZ59hEKGNMo4jIdFUtiLetUZeUIvKGiMyI/CwTkRkJ9lsmIj9F9tt9GbyBvC4XvVq33mXJHsL3Du4bciJepwuv04nL4SDD5eKifgdasjfG7FKNXfHq/O2PReQhYFstux+nqhsbc76mxh8M8vmSRczbuIHuLfMZsU9vMt3xx/iPX7yIh6Z8w6qiInrmt+LBYSeyobSMikCA43v0pHcDR94YY0xDJaVLR8Kd5CuA41V1YZzty4CChib83dmlM3XlLzw05VsWb95Ej5b5/GHQkRzRNfGEqi3l5Zz15qtsLCultLKSLJebLLebd867iK55eVH7fvjzPG6tUYcnw+XimVNO56hu3Sny+RizaAFFPh9H7dU94bDLikAlz/5QyHvz5yII5/bpx+UDDsYbb5EVY0xaqq1LJ1kJ/2jg4YQnEVkKbCE8F+oZVR1dn+MmI+GX+P28PXc2369eSY+W+VzU/0A6t8iN2ufrFcu46uMPYhLy48NP5fgePWseEoBbPh/L+/PnUlltFqtDhMM6d+GVs86ralNVBr8wmrWlsUXU9m/TljuPPo4rPnoPNFw2welwcPq++/OP44dG3WwOqXLOW68xb8MGfJGyDRkuFwe078BrZ51ni5obY4DaE36dl4Yi8jkQO5Ad7lDVDyKPLwReq+Uwg1V1tYi0A8aLyHxVnZTgfCOBkQDdujWuZMHGsjJOe/2/bKuooDwQwO1w8uKMH3npjLOj+sv/PmlizE3UikCAv389MWHCH7NoQVSyh3BS/n7VSnyBQNVVtz8YZH2cxUsAlmzZzFWffBA1d6AyFOKjn+dzQo+eUbNtv16+jAWbNlYl++0xzl6/jqmrVnJ4l671/KsYY9JVnTdtVXWIqvaL8/MBgIi4gLOAN2o5xurI7/XAe8DAWvYdraoFqlrQtm3sJKeGeHTqZDaVlVWNr68MBSkPVHLz+DFU/2azOEH1yaVbt5DoG1DiK2pBgS+XLeH+b7/mtdmzyHbHX7A8LyMjbl2cskAlb82dE9X249rVMZPKIDyy6Me1qxPEYowxOySj83cIMF9VV8bbKCLZgENViyOPhwH3JOG8dfp8yaKYq3AIlx5eX1pK+5wcAFpnZcVd0apVRmbCxH5q7/14a+7sqGJqThGO6NqVi997i3kbN1BWWUmG00VQQzhFCFZL7pkuF2fu24dXZs+Me/xgjbg75LQg0+WKKR+R4XLVe3EWY0x6S8ZMnwuo0Z0jIp1E5NPI0/bANyIyE/ge+ERVxyThvHWquaDIdv5gkP/7/jtWFYUXAb/20MNilhDMdLkYVZDwiwg3H3EUPVrmk+1245LwTNz22TkMaN+RORvWV12NVwQDVIZCUck+2+3hjqOO5feHHxH3G0SW282ZNda6PbnXvjELngvgdjg5ce/6F1MzxqSvZj3x6vkfC3n4u2/jFlVzOxx4nC7eOvcC9m3dhmemT+PJaVMJRK7Grzz4UK4feHitN0NDqny9fBnzN21gr7x8TujRk9PfeIX5ddTlyXS5mHblNWS53UxYuoTrPvuIkCr+YJAst5uju3Xn8RGnxhRkm7dhPdeP+ZhVRcWA0i2vJY8PP3Wn6usYY5qnXT5KZ1dpbMIPhEL8YeynjI907cTrLx/UpWvVqJrZ69fxdOFUVpeUMKhLVy4bcDBts7IbdM5TX/svczasr3WfHI+HB4eexLDIlfm6khI+WjCfrRUVHL1Xdw7t1LnWD5pVxUUI0KnGaCNjjGnUKJ09mcvh4LHhp7B482aG/u8/cfcpXL0KgInLlnLNpx/iDwYJqTJ3/Xpemz2Ljy68OGYYZ20u6Nuff3zzVZ2LpVf/6Gmfk8NvD4773yeuhsRjjDHbpUW1rp75+TF99Ntluz2oKrdPCE+M2v4twB8KUuTz8fB33zboXOf3O4AjunYj0+XGmeAqPRAKMbjrXg17E8YY00jN+gp/u3DlzP68MWdW1NqzGS4XvzrgQNaXlrKlvDzmdSFVvl6xrEHncjkcjD7lDGauW8u0VSv5ctlSfly7mspgCLfTAQgPDx1Ojif+UE1jjNlV0iLhA9x25NGsKSli0vJleJwu/MEAQ3rsze8GDsIXDBJKcCsjz5vR4HOJCAM6dGRAh45cecihzFy3lknLl5LpcnNK733pYMMojTEpkDYJ3+ty8cwpZ7CyaBtLt25h7/xWVTc93U4nx/fowZfLlkaNq890ubjioEMafe4D23fgwDirbhljzO6UFn341XXJzeOobt1jRrjcN+QkDu7QiQyXixYeD16nk/P79uf8vv1TFKkxxiRX2lzh1yXX6+XVs89j6dYtrC4uYt/WbWmTlZXqsIwxJmks4dfQo2U+PVrmpzoMY4xJurTr0jHGmHRlCd8YY9KEJXxjjEkTlvCNMSZNWMI3xpg00aSrZYrIBmD5Lj5NG6BBi6vvZk09Pmj6MTb1+KDpx9jU44OmH+Puim8vVY27XGCTTvi7g4gUJiol2hQ09fig6cfY1OODph9jU48Pmn6MTSE+69Ixxpg0YQnfGGPShCV8GJ3qAOrQ1OODph9jU48Pmn6MTT0+aPoxpjy+tO/DN8aYdGFX+MYYkyYs4RtjTJqwhA+IyAARmSIiM0SkUEQGpjqmmkTkehH5WUTmiMj9qY4nERG5SURURNqkOpbqROQBEZkvIrNE5D0RaZnqmABE5KTIf9dFInJrquOpSUS6isiXIjIv8m/vhlTHFI+IOEXkRxH5ONWxxCMiLUXk7ci/wXkiMigVcVjCD7sfuFtVBwB3Rp43GSJyHHA6cICq9gUeTHFIcYlIV2AosCLVscQxHuinqgcAC4DbUhwPIuIEngCGA32AC0WkT2qjihEA/qiq+wOHA9c2wRgBbgDmpTqIWjwKjFHV/YADSVGslvDDFNi+BFYesDqFscRzNfAvVfUBqOr6FMeTyL+BPxH+ezYpqjpOVQORp1OALqmMJ2IgsEhVl6iqH3id8Ad7k6Gqa1T1h8jjYsKJqnNqo4omIl2Ak4HnUh1LPCKSCxwNPA+gqn5V3ZqKWCzhh/0eeEBEfiF89Zzyq78aegNHichUEflKRA5NdUA1ichpwCpVnZnqWOrhN8BnqQ6CcOL8pdrzlTSxZFqdiHQHDgKmpjaSGI8QvtAIpTqQBHoCG4D/RLqdnhOR7FQEkjYrXonI50C8lcTvAE4AblTVd0TkPMKfxEOaUHwuIJ/wV+pDgTdFpKfu5jG1dcR4OzBsd8ZTU23xqeoHkX3uINxN8crujC0BidPW5L4dAYhIDvAO8HtVLUp1PNuJyCnAelWdLiLHpjqeBFzAwcD1qjpVRB4FbgX+srsDsXH4gIhsA1qqqoqIANtUNbeu1+0uIjKGcJfOxMjzxcDhqrohpYFFiEh/4AugLNLUhXC32EBVXZuywGoQkUuBUcAJqlpW1/67IZ5BwF2qemLk+W0AqvrPlAZWg4i4gY+Bsar6cKrjqU5E/glcTPhDPINw1+y7qvrrlAZWjYh0AKaoavfI86OAW1X15N0di3XphK0Gjok8Ph5YmMJY4nmfcFyISG/AQxOqCqiqP6lqO1XtHvlHvRI4uIkl+5OAW4DTmkKyj5gG9BKRHiLiAS4APkxxTFEiF0DPA/OaWrIHUNXbVLVL5N/dBcCEppTsASL/H/wiIvtGmk4A5qYilrTp0qnDlcCjIuICKoCRKY6npheAF0RkNuAHLt3d3TnNwOOAFxgfzmFMUdVRqQxIVQMich0wFnACL6jqnFTGFMdgwlfQP4nIjEjb7ar6aQpj2hNdD7wS+WBfAlyeiiCsS8cYY9KEdekYY0yasIRvjDFpwhK+McakCUv4xhiTJizhG2NMmrCEb4wxacISvjHGpIn/BzikH04tfWbiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1], c=labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here, we are probably satisfied with the clustering. But in general we might want to have a more quantitative evaluation. How about we compare our cluster labels with the ground truth we got when generating the blobs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "[[ 0  0 34]\n",
      " [33  0  0]\n",
      " [ 0 33  0]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import confusion_matrix, accuracy_score\n",
    "print(accuracy_score(y, labels))\n",
    "print(confusion_matrix(y, labels))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y == labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Even though we recovered the partitioning of the data into clusters perfectly, the cluster IDs we assigned were arbitrary,\n",
    "and we can not hope to recover them. Therefore, we must use a different scoring metric, such as ``adjusted_rand_score``, which is invariant to permutations of the labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.metrics import adjusted_rand_score\n",
    "adjusted_rand_score(y, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Clustering comes with assumptions**: A clustering algorithm finds clusters by making assumptions with samples should be grouped together. \n",
    "\n",
    "Each algorithm makes different assumptions and the quality and interpretability of your results will depend on whether the assumptions are satisfied for your goal. \n",
    "\n",
    "For K-means clustering, the model is that all clusters have equal, spherical variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**In general, there is no guarantee that structure found by a clustering algorithm has anything to do with what you were interested in**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can easily create a dataset that has non-isotropic clusters, on which kmeans will fail:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3gUVdvA4d/Z2Z6e0HsHFZUuKjbADooFRVFsiPW1vZ+9vfbeu2IHBMUGqCh2UURpSi+KFCmBkL51Zs73x4RAsrshZTf13NeVSzMzO3NmgWfPnnnOc4SUEkVRFKXhstV1AxRFUZSaUYFcURSlgVOBXFEUpYFTgVxRFKWBU4FcURSlgbPXxUWbNWsmO3XqVBeXVhRFabAWLly4U0rZvPz2OgnknTp1YsGCBXVxaUVRlAZLCLEh2nY1tKIoitLAqUCuKIrSwKlAriiK0sCpQK4oitLAqUCuKIrSwNVJ1kpTIM1CZOHjEPgMMMF9PCLlRoQts66bpihKI6MCeQJIaSJ3nQv6eiBkbfR/igz9Bs2+QAhnnbZPUZTGRQXyBJCBWaD/Beh7bdXByEbuPAVp5oDWErzjsXlH1VUzFUVpJNQYeQ1JIxsZ/BmpW3n6Ut8M+XdQNojvFgTjb5D5oK+Bgpswd45EmsW12mZFURoX1SOvJikNZMHd4P8UhBNkGOnsByITCFb+RPoaZOGDiLQHEtZWRVEaN9UjryZZ/Db4ZwJBkIVAAEILIDgHqMqqSxL8M1ArNSmKUl2qR15dvncAf7mNIUBU42QhrOBfndcqitLUqR55dcnCGDtsgLvcNido3Yj5djv6I4T6o1AUpXpU9Kgu52FEffu0duC9AHCBSLH+6xwEmVPBc1b5g0EkI1L/l/DmKorSeKmhlWoSKTciQ/NA+oEwVlB2IlLvR7gOQSZfCvpa0FohtLbWi9LuxUy+AgqfhNDvIFzgGQP2TnV4J4qiNHQqkFeTsHeAZp8jfW9DaBHYuyKSLkLYu1r7bang7B/5wuA8CHyJNS5uQvHTyMAnkDUNIcoPySiKouybCuQ1ILQWiJQbK328lH4ovAcI7LXRD/p6pG86Ium8+DdSUZRGT42R16bQn4AWZUcAAl9U+FIpA0hjC1KGEtI0RVEaLtUjr022JMCMsS816mYpDWThY+CbAggQGjL5amxJFyNl0KrnYstEaC0S1mxFUeo3Fchrk/0Aa+an9FN20pAH4R0b9SWy6BnwvUfpcIwECp/BDK2G0Fcl23SkcxAi/SlrbF5RlCalwQytSDMPs+h1zLwbMYvfQpoFdd2kKhNCIDIngq0FiCTrByckj0e4hkQcL6URY+KRH4Ifgyy2fghC6Fdk3nW1cBeKotQ3DaJHLvX1yJyzQAaxxpO/RBa9BFnTEfb2dd28KhH2LtD8e2s6v8yzJgNpzaIfLP0l91wZYQj9jjS2I7SW8WquoigNQIPokcuCu0EWsCfbIwAyH1l4f102q9qE0BCuQxDu42MHcbB67LYK9kccbwczp+YNVBSlQan3gVxKE0K/EVmIyoTg3LpoUq0RQkDKbUTPdIlGgr1LIpukKEo9VO8DuVVIKsYIUCNcaUdKExleaw0nSYnNc2LJ1P59FdTyQPL/qUlFitIE1fsxciEE0n0yBGZhTYXfzQnuxrW6jgz9hsy7vuQhpglaK8h4EZHyX2TgS2tMvTR90Qki3Upp1FoiksYjXEdWfH4pIbwYQvPBlg7uExG29D37zSJk0fMQmAkI8JyCSLoKYUtK2D0rilJzoi7qYA8YMEAuWLCg0sdLsxC560Iw/irZIMGxHyLj9UYTZKSRjdx5bElq4m4CRAaixY9g7kAWPAjBH61vIp7TECn/RQhPxectPZ8TmXdNyXBUEKuol0BkvIZwDrTy1XNGlV1nFCfYuyGyPlLVGRWlHhBCLJRSDii/vd73yAGELQWypkP4D9D/Bkc3sB9ojSE3EtL/Mcjyy8NJkH7krotBX2Htdw1DpN66z8wUafyLzLsZwous82gdwdjCngfGfuv0uVdDi1+sDwhjE3uCONb/G/9A6GdwHRGvW1UUJc4aRCCHkgd/zj7WT2Okb6Xs0NFufgj/tufX4BfInT8hm32BLcZsTimDVrqmmUPpUIzxd4wLhyC8FPRlIH1RThaA8DIVyBWlHlPfl+sLLbmSB0prUYsdQ5GhJdEPCcwpmSgUoxxA1Ou3A7yR24UbtIaVq68oTY0K5PWFvqGKLwgh86620jPLMzZG711H5QTHgeA+waqPXiY7xgbCA+5jq9g2RVFqkwrk9YWobK74XmQR6KsjN9syKvFiNwgvIuN5a4KS8CCypoHjIMBh/TgORmROQwiXdV4pkfpfyPCa6B8giqLUiQYzRt7YCc/pyMB3RNZV2ecrgZIg6/8Qip4Hc/s+XmMH59GI9HvLpB8KeydE1geldWz2LsAlw6uQeVeDsQOEAJEM6U8jnBEP0BVFqWWqR15fOI8AzyishZsdVGo2p0gBe08Aa6WigvvA3AIY+3ihDqHZyB3DMAuftYpzlTlvEpgFSNNaYFrKAHLX+daQDX5r2MbMRuaORxqqJICi1DXVI68nhBCItHuQ3nMhNBeJ15oA5P8IzGLAXjIhKISVA65B2nNI/8fI4oklOfZVnBMgC6F4ItLcjkh7AADT/xkU3Gtlq2AgXUPBdRRQPjUSkAbSPwORfFGN7l1RlJpRgbyeEY6e4Oi555Fj8hUAmKYJgY+tXHr7/gjPCGTRy+B7l6oPx+wtAP4PMMOrIelCyL+NMkvRBb8F/S+Q0VIjg2Bm1+DaiqLEgwrkDYAMr4bcy0p65AL4HGlLAd/bWLM098HWGswdRO1V76b/Cfk3RTmmZFJQ1KEeL8I5qFL3oChK4sRtjFwIoQkhFgshZsXrnApIGULuGmeNfUtfSR2WAsi/nkoPpZh54L2YfRfeihHohQscBwN7lwNwg6MX7KO+i6IoiRfPHvm1wEpArTUWT8EfKDttfm+VXYjZD/5pVHkMfTepQ/oLEJhdch4dPKMQ3rGI6qRNKooSV3EJ5EKIdsDJwAPADfE4p1LC3GVVQqwpWVzNF3og+TJsWioknWX9KIpSr8SrR/40cBOQEusAIcQEYAJAhw4d4nTZJsA5iCpNtY+pgvHxCMJKQdTaIZIuA605Zv7NIMMI9whwHa2qISpKPVLjf41CiBFAtpRyYUXHSSlflVIOkFIOaN68eU0vW23SyEYaO+vs+lUl7J1L8stjLaKRiKENac0aNQuQoeXIXZeC/xMIzELm34DMu566KH+sKEp08ehWHQ6cIoT4B5gKDBVCTIrDeeNKhldh7jgJuWMocsfRmDtPQ+r/1HWzKkWk3gupDxH9C5SNsg8h48jcAv6JWOmNJYFb+iA4G7lzBGbw18RcV1GUKqlxIJdS3iqlbCel7ASMAb6VUp5X45bFkbUwxXlgrMN6QBgCfQVy1zlIWdkHhnVHCIHNOxKR9SnY2oLwWlPkRQqk3k98hl6qQoKxFnLHYe48T9VdUZQ61jTyyCOWiQNr+CAAga/Bc1JdtKrKhKM7NP8W9JUgg+A4ACGcmOYWKHqVPROD3FiphkESHuT135C51yIyn0vsdRRFiSmugVxK+T3wfTzPGQ/S2FJuCbXdO4Jgbq39BtWAEAIc+5f+LvXNIEPg6Asyv6Ts7EngOhp2XQ7mmsQ3KjQHaeQgtKzEX0tRlAhNokcuHAchhTeyRrdwlpRtbZhkcC4y90qsIllhK4hrHRCe0xE2L2bGa5BzVO00Rl8HKpArSp1oGjlkrmOsNSvLZH64wH4AOMqWYZVmEdI3FbPgIaR/Rr0dQ5fSQOb/H1ZdlJJhI+kH/R+rEiIg7OkxXx9fArTWtXQtRVHKaxo9cmGHzCnI4lfBPwOEDTxnIJIuKbOAs9TXI3POtoZc8Fu9+KKnIWs6wpZZdzcQjb6upEJheUEIfAbJVyCEB+kYVHbNz6i0kh+TquWbl3AMQNituQFSmirHXFFqWZMI5ADCloRIuR5Sro95jMy/qaQw1e4NPjBCyMInSsu81hvCHXvGp9iTjigyJyJzxoEebX1PzcqASboabClQ+KCVP14VWi9IuRmz+D0ofgvM9UiRCkkXg3c8IvQj6KtA6wDu40tXG1IUJX6aTCDfF9PYYpWIjaCDfzbSdSQYW631LR39yvTk64Kwd0Rq7UtSKveanCM84DkTs+h1K1tHeBBJ45COF8HcAPom8E8FMwecR0LSeAR+ZHiVVVOlqozVsOv0sttkARS9DEVvIoVufSAKDxQ+DJnvI+ztanTviqKUpQL5bkUTK9qJzL/Fyg4RDrAfCJkT67x3KTJeRO4aW1JHxbR66K4Twfce6H+zu664zF8OnlOxpd0Lzv7gHQWA6f8Cck5FEioZTtrXykLRxJrhGbB+du+WPpABZMGtiMx3q3EdRVFiUYOZu4V+qGCnLAmWYSsghf+wVuWpY8LeEdH8e0T6s4jUuxHNZiJcQ0D/hzKLQ+AH/8dIfWPpFjP0J+TfbA0lSR/VC+JVZUJoATLq2L6iKNWlAvluIma9rygC4PswYU2pCiHsCNcRVsqhvRMy+BPgi3KgZgVRswAz70bYNZqywV5RlIZKBfISImkckTVLBLEXY6jceLKUZu0WmNJaEnXETAaQaNYiFYHPqXZt8kqL9t7ZwDkQIdwJvraiNC0qkO/mPg28Z2ItbJwMeMHeA7RO0Y/X2lZ4OhlajLlzFHL7fsjsvpj592HWworzwnMm0f9YTSh6qWTZtmjrb8abpOyHhQDs4B1XC9dWlKZF1EU50gEDBsgFCxbU+nUrQxrbIbzc6tna9wd9FTLndCLHkN2I9CcR7uGR59D/Qu48jahDF/bukHI3Nlfi1ro0d55ipfxFXhwrX7wS63wmhABciIwXrbH8vUgZBgzVW1eUCgghFkopB5Tfrnrk5QitJcI9FOE4wEoxtKURveZ3AFn8ZtRzyOLXiDn+rK+F3PMwd4yqgzK6dva9bmciSSCALLhnzxazCDPvv8jtfZDb+1jfYsLL6q6JitIAqUC+L2a+lXIYdV9u9O2haJNvyjFWIHPGJCaDw3U00ReisFnfCGIuUlFLjE2lpQ9k7ngIfIk13GOWlBc+H2k0rGJmilKXVCDfF3tXor9NGti7Ic0oMyErnQEThMBXNWhcdCLpIrBlArvz3AXggdQ7EJlvg/vUKpzNSdx78cIN2JHhlRBeScQi0jKM9E2J7zUVpRFTgXwfhHBCyt1YGS17BzQTQj8hsw/D9H1c9kXuYZU7ufSD8W+cWrqHsGUgms2E5MvBcTC4jkNkvoHNeybClmyVKqCuJjO5wXOuVY/F2GilRUYIgV4L5XcVpZFQgbwSbN5TEJnvgOt49oyX754kFICCu5H6utLjhfdcKhcoXdaU/wpIfSPS9z4y8EWVhmGELQ1b8lXYsj7AlvEcwtl/rxtKjz1chCip1eK12ucaDsQ6thrcwxEpN1j/b+8BMloGjcv6AFIUpVJUIK8k4TwY4TmlZFigvDDSN33PsbZUSP1fJc4aRPo+xsx/COmbjDQLSvdIKTELHkbuPBlZ8AAy/zZk9mHIyoy/7+tehB2SLicyb94Nac9A8vUl6ZUmBL8iYuijJqRuXZ+ShaVdR1jXLWWz6sN4x8TvmorSyKlaK1Uhi4g+kcYAmY8ML0MWvw76Rqu3uU8mBGdap/Z7oPBJyJyMcPSC0Fzwv0dpquDutY9zL4MWP5cGw+oSSZciRRIUv2QV0LJ3QaTcCs7DkTmjwFhPtUra7kvw+7LtSH8aWfSiVR9G+qyVjlLvr39lgxWlHlOBvCqcg6NXCBResDVH5pyL1Xs1QV9exZP7rdGa/BsRzWYife9HX56OEIQWguuQqrd/7yYLgUgaC0ljy2yX4RUlk4YSEMShzJi41DciC5+H4DdAEeCG8J+QfwUy402E1iIxbVCURkYNrVSB0FpB0gTKDkl4rPFc3wdYueO7a4RXc9FjfT3SyCmpRhi1FcR1qKM8YwvIBOaae04Ddi/iMQqCnwKFWF85/IAP9L+Ref9JXBsUpZFRPfIqsqVcg3QOQvrfB1mMcI9AOvrAzpPjdAUJwobwjECG5mMFt70Z4IyY2BU/jt4kbOanaA1JVyH9nyGL3yh5WBxjqCq8AmlsA1tzMLNBpCJsSYlpl6I0cCqQV4NwDUa4Bu/ZYBYi99kD97BnuCJWrRMbOA5A2DKQ7pPA/wmEF5WUmbVbP6kPIkT5h5TxI7RWSOeh1hh9XLkg9XrYMRQpbCVBvMKGIH2fgu/1kiXtTKT7ZETavXVeB15R6hsVyONA2FKQrqMg+B3Rx5Y1cB6GSH8IwsuRhY+D/lfJPhOrF+4EkYRIe8I6p7BDxkQI/YgMfAe2dITnjNK1MRMq/UXIPoTIbwM1EYT8m6z/rVR5HxsUP0+ZbweBz5HoiPQn4tguRWn4VCCPE5H2MDJnDBhro+w1QGgIWzq4Dke4DgesFEMZ/Nqa3am1gaQJZYYPhLCB62iE6+jauYkSNpsbmTXFmj4vA4CwHvJq7cHYZP1OgMSVwnWD1tqqS1NGEAKzMH2HIDyj63y5PUWpL1QgjxNhS4GM55A7TyHyYaQbtC7I4M9g72o9NAVk8QvW2pbYQQjwTUVmvgn2XhBeCOYucPRBaC1r+3YQjgOg+dySoR0/OPojbEnWtHp9FdIshsJ7E3PxrGmQe0WMnRIK7rWya1LvLlkP1I0oyYaRoYXIomdAX2e918nXIhL5TEFR6gFVxjbOzNyrIfgje6of7k63MyhNEnKdDN7TIPcqIoYvRDpIJ1Y6HkAYvOMQKTfVix6oNAuReddAaAHWPcU7TTEVhFnBg9DdbGDLsj7shNuqc+4YAHlXUbbypBuR8VLptyBFachilbFVgTzOpAxbk4J871k9WRkGyj/YE6B1izEME42A5GuxJV8Z59ZWnbnrUgjNI6EpkNXiwSqLkBe5S+uKrfkXtd0gRYk7VY+8lgjhwJZ8ObYWP0DGW0QGcQAJxt9VOKuEomeRRjbSLMYsehVz52mYOeOQga9qbSk5ae6qp0EcrG82eYQCULBLo8xbUqX3WlEaHjVGnkjh+RXsNLB6kZXNDJFI3zQIzrZKAJRkc8j8P6xqgqk316ytlWHmgrCDrI+BHEJBwWXDepC92UlmS51rHtnMwGMKQWTUddMUJaFUjzyBhK0FFdbytmWUVBqsDAmhX0HfTJmUPOkH3yRr8kyiaR2I/dkvqOu/TnpIkL3ZiR62kb3ZyX3jO7FuWTokX1qn7VKURFOBPJHcw6w6LLGYBZD6IDiPIPpycnvzgOknag9eOCC8uPRXKWXJMEyU8eIaEMIBKXcS/cNJWotWa/vF9ZqV5ffZmPRUS/Twnr/SwYDg/Zf7I7wX10mbFKW2qECeQEK4EZmTKzjCj3Adgy3zdXAOIXYNc6eVZ25ujLFfWhkcgAwtQe48DrljGDL7cMyc860FpePE5j2VmH9tZGFJnnn8hEOCHVvshMuN5vw5z8OKBV78xTa2/OPk+Vvb8OHLzcu1R/DPane9yPZRlERSY+QJJhz7I7XeYERZUNiWBcKDlMGS+t9GyQ5HSc/Wb23zjIDQnyVVCaNdJA0cA5DGNmTuhSVT+kuEFyB3nQ/NZlsTjOLB3gP0lZHbtbbWf40oy9+V4cBKLdSRUgCS8rFWSpj0RAumPtcSwxAIIRl6eh7/fWoTNhs8eUNHtm7Y91R9zb6vbzqK0vCpQF4LRNotyF3jKZ/fTPKNCCEwc28oyT3fqxaLuRaR9SnC3gkpQ8iig9kT6MudP/MdMHchc68pG8TBeo2xEel7zypbG4/7SbkZmXt5lPu5FcytUPgEsR/iuiHzQ9YtXs765SZ60QzmzdrCkrkpOD0mJ56bw/n/t51Zb2cy6alWpZUYJYKvP8jAV2Tjrokb2LqhcgtIh/z188GsosSTCuS1QDgHQebryMLHrBmHWltE8nUI9zCk8W9JEC9XcVD6kTlnId3DwHshsYI4gAyvg8I7rAUiojKh8EGk6xCEvVvN78d1WMn9PFkyg7KTNYPSNQQpTaT+F/inE1kczE7QfgfXD36VLeu2YegG4WAIRCrSFAQDNj56rTlLf03mn1WuKOV0Bb98kUbQL8hobpK7Y9+97VadVU1zpfFTgbyWCOdARNb7kTv0jVbBrGj1x2WeVQHRP7Pik+dfiTVUUVE+eRiZezU0+yIuY8bW/bwXuV3YEGn3IJP/gwwvh8DX1iIb9s6IpEt47tJv2bhiM+HQ7m8fokyz9ZCNlQsreECMYPbMxykqfI/YVSQtLq+L0f93apXvTVEamhoHciFEe+AdoBVWKb9XpZTP1PS8TYa9SwWLSIDVE4/dG7dUchELYz0EPrPG3BNMaM0Q2lHgPqp0m5SS76fdv1cQj/nqCvdOvGUK4WDF53B6nFz+5AX0G1bx4taK0hjEo0euA/+VUi4SQqQAC4UQc6SUK+Jw7kZPaC2R7hMh8Dn76mHWnETm3w7OfgitTYKvFZ0ernltln0FcYCn595H975danwtRWkIapzGIKXcKqVcVPL/hcBKoG1Nz9uUiLSHwNaqlq7mR+4YgVn8tjXlvhYJIejcuxbqqQPTHv20Vq6jKPVBXPPIhRCdgL5AxNx0IcQEIcQCIcSCHTt2xPOyDZ4QdrCl1eIVi6yHn9lHI4O/1OJ1oVvfTrVynZ8/ms9Zrcdz7ZA7WPT1n7VyTSklv8z4ndtOeoDrj7yTmS99SSiY6G9ZihLHQC6ESAY+BK6TUhaU3y+lfFVKOUBKOaB58+aRJ2jq3MeTyPlZ4ZBg3VIP2zY5SrZIIIDMuxYpdQK+IL99sZiFc/4gHEpc8OnevysuT2TqoOaw4XA5Sn+3aTV7L/SwQe72fFb8spq7Rj3Cj9Pn1eh8lfHK/73NQ2Of4ffZS1g2dxWv3Pgu/3fM3XEZTlKUisQlcgghHFhBfLKU8qN4nLPJER4q/dCyir79KJ2zDjyAG8/oyqVH9eLaEd3Ylb378YjOH1+/x+iWl/DAOU9x7+gnGN1yPH98vzwhbRk29ogyARusoJ3ZKoNb3v0P/YYfRLe+neN6zaAvxEs3vJ3QKpHb/slmxktfESje8+A66AuyfulGfv74t4RdV1EgDoFcWLlsrwMrpZRP1rxJTY80i6Dw8UocWfW0wdVLPDx9Y3t8hRq+Io1QwMbqJV5uO8d6EKiHdd7933QCxUF8BX58BX6K833cecrDFBeUn1xUcykZyTz5471069sZu9OO3aFx4BH78fTc+znyzEN55Ks7GX7+kZhGfD/UcrflEigO7PvAalr640rs9sh/ToHiIL9+vihh11UUiE+P/HDgfGCoEGJJyc9JcThv0xFebJWHjcoGSZcjWizGL+7kr3VnUqyPx5rmvodpwvSXs7jw0F6c228/XryzDQW5Gh+/1pxQoOwHgDQF61e6WfJzEoV5JisXuqNeOVE9yc69O/DSwkeZtuVVbp96PUbY4NrD7+Des55gw4pNfDtlbtyv6XQ7cUYZ0omX1KxkhC3yn5Pm0MhqlZ6w6yoKxCH9UEo5l+p0FRVkeDnSNwX09bFrfDsPQSRfz8RbJ/PJs19gd2rooX847oJLuPLe1WjGz0CIey/uyKKfUgj6rdmOn72bxfw5qaQ3213PJNKr97Tmn1VeDD2y92voBsX58e+R7+3XWQt59srXCPqse5/70S5+n70EfZ955lXj8roYdc2JaFri6q70P+7gkiGjsqUJ7HaNEy4ZlrDrKgqo6od1xvRNQ+acA/4PIbyA6DnkHkTyVXzy/Bd8+vxsQoEQvgI/oUCYOe8s5PVHDsPWahmTX72JhT+klgZxsGZIZv/rZNM6F9FnfAr+WubF0GN8BgvBgOMPjsOdRmcYBq/8953SIA4gTUkwzsMfml1jxOXHcsE9ZwPgL/KTuz0v7uPldoedx76+i+bts3Anu/GmevCmeLj5nf/QrnvruF5LUcpTU/TrgDSLoOAByhad2h1YnCVT9kOQfCXCOYgPHr+coK/s7M+gP8Rnr8xhwqPnM+PFeYSCkQHZNATFBRX9EUcP4nannRMvGUr7nombDpC7PT/insCqemh3aJi6iWnWfJzc0A2kKQkUB3li/Ev8OmMBCEFGqzRuePVy+h9b9sOqOL+YH96fR+72fHoP6cVBR+1f6ZIGnQ/syOR/XmLtor8J+kL0HNQNZ7kHu4qSCCqQ14XdY+LROoX2Hojky8A5EGHLBKAgpzDqaUKBMOFgmFAgRDxHt44++3CueiaxizEkpyfF7BW36dYKT7Kb5T+vjsu1Pn1+NmsX/c2q+WtLywNkb9jJ3ac9xvPzH6LTAe2RUrJi3mpuPfEBTEMS8odwJbnYb3B3Hph1Kw5n5QKyEIIe/bvGpd2KUllqaKUuiCRiFrjS2iPcx5cGcYBeg7pHPbRV5xa4PC76H9cnxnXAm15RAapIdqedEZcdm/DFGNxeFwcesV9EHXKX18XY28/g8W//R69Dot93VRmGwerf1kVM7Q8Hw0x/YgYfPTOLM1tczHVD7sRfGCDoCyKlJFAUYMUva/j8ta/j0g5FSRQVyOuCo4+1LFoED8J7bsTWyx4fhzvJhc1mRT0hwOV18p/nx5O7PY/u/brgcDvKdsoFjL3tDAKFVRtzNg2T/QbHJ4BW5JspP7F07krKd8qPu/BojjzzUOwOO0//dF/pPdeEQGB3Rn75NA2TP35YwZu3T6UgJ/piGEFfkC/f/K7GbVCURFKBvA4IYUNkvAG25lbvXCQDLmtM3HVIxPHd+3Xh+fkPcdTZh9O2R2sGjxzAE9/dQ87WXZzX+UqmPPAhNiGw2Ww0b5/FkNMP4aWFj1ZrerhpmNx16iN8M/mnhM1IlFLy6o3vEvJHtm/lr2tK/1+za7SMUU9cs9vQouRtRzPkjEOiFtrSHBrbN+wgEGWsXlEaEpHI2W6xDBgwQC5YsKDWr1vfSGlA6HeQ+WXGxCtjx+YcLux5TcQKOHaXnTOuG8GCL5ewfukGTKN6f77uJBc9B3bjka/ujPtyaQFfkFPTxkWd9OPyOJlVvGed0++n/czjl7xYJrvF5XVx/EXH8PU73+PbxzeOI886lNsmX8tTE17h+2k/l55H2ATS3Pd74/+jglAAACAASURBVPK6uPTR8zj1yhMqe3uKkjBCiIVSygHlt6uHnXVICA1cg6v12p+m/0rEuASgB3U+eHxGjWdGBoqDLPt5FXeNeoQTLh7GYacMqDCgF+wqZP6sRRiGySEn9SWjZexJMC6PE2+Kh6K84oh9zdpllfn96LMPxzRNXr91Ctkbd9KsbSYX3Hs2Q88ZwreTf4p4vU2zccyYwznm3CH0HXogq39bx8W9riN70w4M3cSd5EJzaPgLA8gKF+KwPsz2P6wnJ186vMLjqqswt4jXbp7ED+//gqGbtOjYjIHH92Hk5cfRrkfdlBlWGibVI2+g3n/sU964fUrUyTyJ4E3xcO+Mm3B5XORuzyc1K4UuB3XAk+zhhw/m8eiFz6NpNqSUmIbJlc9cxMmXHhvzfNMe/YR3751eJgXR5XVxw2uXM/ScIVFfI6Us8xB29e/ruGPEQwQDIYLFQcySHrbdqaHZ7Vz97MW8cO0bZeqfaHYbUlLhB53doXH0OUM48eKhJQ9k4//g19ANLj3ov2z5axtGuOzCIU6Pk+tfvYzhY4+M+3WVhi1Wj1wF8gZq85otXNL7esxaCuS7CWF9ERBC4HDZOee203nvoY8jhnicHiev/fkEbbpGr7MupWTqI58w7ZFPCPqCJKV5ueiBcyoM/gD/rtvKO/d8wNIfV5DVJpOzbjyVP75fxuevfR0xDu5wOTANo0ofdsImuG/GLRxyUr9Kv6Y6fpnxOw+f/yz+GENDLq+TD7ZNxJPsSWg7lIZFDa00Mu16tKFb386s+f2vKr0uKd1LcV71p97v/tyXUhIKhJl83/SoNUYM3eC7qT8z9vYzop5HCME5t5zGWTeeQqAogCfFgy3Kefa29e/tXDngZgLFQUzDZMemHB4Z9xxpzVOjPsw0dKNKQ0yaQ+OsG09JeBAHWP/nxjLfFCLaYtf44/sVDB7Rv8z2LX9tY82Cv2jevhn7H9oDIQS5fj9rd+XQOjmFFklJODQNmxCYUvLJqhW8++cS/LrOiO49uKhPf5Kcias5o9QNFcgbsMsfv4Bbjr+PUKDy2SnDxh7BzBe/itsUdT1sIERksDQNk3AlsmY0TSMpLalS15p03/TSIL5b0Bdk56adUY+3aTZsdht6JZaGA+gztHfpVP5Ea9u9FW6vC39R7Ie1TveeSUiGYfD4xS/y4wfz0BwaUkLzdll0ePA4pv27BpsQBHQdCXjsdsYd3Jccn4/P1q7Gr1v3/0JeHjPXrGbGmPNw2dU//cZEpR82YL2H9OL4i4bi8joRNoHdoZWk5cV+KLlu8T+4k1xxbUe0DwWn28Fhpw6M63WW/rQyag9bc9ijLlaR2Sodb3L0yo7lCSE4aMh+CS2stbfDRg3Cm+pBxMiT1+waBx21P2C9v6/dPInvp/5MKBDGXxggUBRg89otfHvzhwQNA39JEAfw6zpvLlnEx6tWlAZxgKCh829BAZ+vXRPlikpDpgJ5AyaE4JoXxvPEd/dwzq2ncd5do3lj5TM8Nfe+mMF85785BPzxzZvObJ2Oy+PEZhMIIXB5XZx4ybC4T1Vv0aFZ1O1SSo44czBOtwN3kgtvioe05qncP+tWrn35skrlm0spWbPo77i2tyJOl4Nn5z1I/+EHlU7kstmE1f5UD/fNuAW7w044FObWEx/go6c+Qy/3UFQaEuf6QmxFkd98QoaBEeUD1qeH+WnjP4m4JaUOqe9XjUDPgd3oObBb6e+tu7SkZadmbFm3veyBAnZsyqlU/nRVFO4q4rbJ1/Ht1J9weVyMvPw49j+0Z1yvATDmltNY9du6MpkuTreDwSP6c/Pb/2HsHWey9McVpDVPZeAJfXA4HaRmpSCEjcqsvtTpgPZxb3NFWrRvxkOz7yAcCrP9nx0s/Wkl3lQvh5zcD7fX+tb06fOzWTZ3ZcyhMCkEtrBZ6bWlHDYbbVJS4nQHSn2hAnkjJISIPuVcRh8GqQzNoUWkye1m6Cb3j3kKd5ILPayzefUW7p91K99OmcvkBz6kON9Hy07NueaF8fQbdlC1rg8w4LiDuerZi3jl/97BCFvZKIeNGsR/J14BQLvurSNKxma1zqDvsAP5ffbifZ5/8Mj++zwmERxOB+16tImaOz77jW/LTIYqz0i2o6dX/uGl3WZjTO/q/xko9ZNKP2yEPn1hNs//5/W4nlNz2DDClev32R0a6S3T2Ll5V8S+WydfGzNPvLL0sE72xp2kZqWQnL7vB6UBX4DTm11MuIKHwg63g49z3sTlie/zg5q6qNe1bF6zJeo+h8fB9gm9yO1cuYfFAnjr1DM4omOn+DVQqVWx0g/VGHkjNPOlL+N+zsoGcbAyWaIFcYCnL3ulRhkzOVtzWf7LapxuR6WCOIDb6+aOqddb4/hRxstdXhejbxhZ74I4wLHjjsTpiSyhm5TmZeLSp5h82xWV/kec7nZz5ecz6Pbck4yZPo1VO3fEt7FKnVGBvBFyeutfQNot4AuSvTF6umBFDN3gkQue4/wuV3H3qEcZ1+0/3D/mScKhyqVeHnbKQF5e8jhnXHcyBx99AFltMwBIa5bCuP+N5oJ7ayftsKpOv24EnXt3wFOSfePyWuUNHplzF226tKRX8+Zc0ncAnr3SCaPlwdiEoDAUojgcxpSS37ZsZvQHU9lSWFBLd6IkkhpaaYR+nD6P+856sq6bEZWwCT7c8QYpGdHK+Mb21t1Tmf74TIJ7zSB1epyccsXxXPb4uGq1pfyU//rKMAx++3wxy39ZTfN2WRxzzuGkZu55YCml5JNVK5m4eAF5gQBHduyIYUpmrVkNwnrA6QuHI7JYHDaNi/r045YhqhRAQ6Gm6DcxNx17D4u/WVZmW8f927H17+1VmkBUXZrdFnVqfKvOLXju1wdJb55WpfOd3uwiCndFPsB1J7mZUfBOgwjI8bZiRzYrd+6gY1o6/Vu3iXgPAnqYXH+AZdnb+e+cLygKRT40Pax9ByadNrq2mqzUkJqi38Q8Oudulny3jI+f/Ry7w845t51Gtz6dyc3OZ2zHK2LOuhRCxGXWZ6z6JtvWZ3NW60u56L4xnHPr6ZU+n6/QH3X77tV8mlIgD+o6l878hAVb/8VWct8d09KZdNpoMjx7arO47Q5apzgIGjrhKBOpHDaN3s2j13tXGhY1Rt6I9TmmN/d8fBN3vn8D3fp0BiCjRRqnXn1C1BVzoPrpiVUhTckbt7/HZ6/NqfRrYq1a1K1f533WaGlsnpk/j9+3bCag6/jCYXzhMOt25XD7t9Hfz07pGRzeoQOucrNWnZrGuIP71kaTlQRrWv8CFADGPzQ2ou53XXj5hrfZvqFymRNXPXMxnmR36YxVm2bDneTimhfGJ7KJ9dL7K5YSNMrm9IdNk6/X/0XIiJ7r/8KJIxl74MF47HYEMLB1Wz4YPYY2Kam10GIl0VQgb4I0u0ZaVtVm99mddpq1j2/wDxQHGdf9am454b59LkvXrU9nXlnyOCeNH0bPgV05/qJjeGnhozEXpm7MYgVrKSVhw4j6rWprUSE/bdxgFdVyOPgrbxc7iiMX9lAaJvWwswkwDINfZy5k3ozfSclK4YSLh/LT9HlWHfFaePBZGanNUnhl8WM0a1vxh0U4FGbtwr9xeV10Oahjkxob3+3a2Z/x+drVEVkoSQ4Hfl3HbrNxSo9e3HnkMaS4XOimyRFvvkZ2cVGZNZGcNhv/O3oYAAPbtKVrZt1/S1MqprJWmig9rHPbiQ+w8rd1BIoCaHYbdoedy5+8gKkPf0Lu9rx6E8w7H9iBV/94Iub+uR/P57GLXgDANE0yWqRx38xb6bhfu9pqYr2wraiQU6ZOojgUwq/rOG02QqaJgNJA7dQ09m/egg9Hn8NPGzdw6cyPCZuRDzztQmAvGTs/uVsPHjn2hNIHqEr9o2Z2NlHfT/uFlfPXEiipe23oJkF/iJf/+zZPz72PsXecQa9DupPaLHGFlDSHRkar2Gt47rZl3TbWL9sYdd/mtVt5+Lxn8RX48RX4CRQF2bY+mxuH/g89XLl6441BYTDIs/PnURwKETIMWiYlc1DLVjg1rUxvO2QYrN65k6XZ29lUkB81iAPoUhLQdQK6zmfr1jBj9crauRElrlQgb+S+m/pz1JVoNLvG+qUbOfe2M3hu3oNccM/ZMTNZaurQEf258c0rScmqeBKQ5rBTsLMw6r4vJn6Nrpcr4yoh6A+x6OulcWtrfSal5PxPpvPRyuX4dR1DSnb6ivkze3vUcXPDNJiy9A+8WuX+XAO6zuSlf8S72UotUHnkjZwn1iIS0poZuVtymheHy44ein/v9vcv/2DRN0sJVrC0GUA4GKJbv85R9+3amhe1+qI0Jfk7m8Y084Vbt7BuVw6hvXrXhpQIKbELgV5umDRkmry/Yhnvr1hW/lQx5fqj5+sr9ZvqkTdyJ004trS29d6cbge9h/Qq/X3wyAGQoMclQV8QX4EfYx/rZ+ohnQVfLom6b+CJfXEnR96HoRscdOT+cWlnfbduV07UjBTdNBFCxGVsu0+r1hHbQobB+rxcCoKxl6VT6pYK5I1cv2EHcvr1J+NwO3Anu/GmeEjJSOL+z24rs6yZN8Valcab6sGb6sHpdljBQau9B19SwgPnPB21h33kmYPp0KsdLu+ebxHuJBcnTziWlh2b11ob61KXjMyoWToeu53L+w+iZ1b0FZSq4qI+ZReenvznEga89iIj33uXQRNf5rrZnxHQ68fDcWUPlbXSROzYnMOSb5eRlOZlwAl9cLoiS6MChAIhFn+7DD2k02dob1b//hd3nfJwmWJVwibivsrQ3sY/fB5n33RqxPagP8jnr33Nd1N/wZPsZsTlxzHktEFNJgVRSsnIqZNYm7Oz9OGlTQjSXC6+u2A8J015my2F0Z8xVIZNCMYccCD3Dz0WgG/W/8U1X8wqs+6nS7NzQtduPHXCyTW7GaVaVPqhUm2/z17MKze+w+Y1W8lslcEZN5zMazdNirliUE2decPIalc0bOwKggHu+eE7Plu7Gt00Oax9B+47ejgd09Pp9tyTmDX89+y221lx5bUAnPH+FBZv2xpxjFPTmH/JZbjsdlyavcl8kNYHqmiWUm0DT+jLwBPK1uT46q3v+fuPDQm53px3f6B1lxaMvOJ4FSTKSXW5eeK4E3n82BMAyrw/HdLS+CcvL+ZrBeAqqVse0KM/1A7qOlsLC2idksq2oijLBWJ9Mxj6zpvkBwOkOl1cNWgwF/fpp/6s6pAaI1eqZcRlx+GKtoBFHP4t5+8o4LWbJzHlwY9qfrJGSggREThvPfyoiMJYu2lCcFTHzswYcx6LJlzJATGqHkpgyJuvcfPXsxnQpm3UB6hh0yQ34MeUkrxggCfnzeWNJYtqfE9K9cUlkAshThBCrBZCrBNC3BKPcyr120mXDuOwUwbg8jhxeV14Ujxktc3kupcmxOX8geIgUx/5ZJ81WMrbuWUX30z+iV9nLaz06kGNxbFdu3Fy955RP0s1m437hg6nW2YWbruD/x01FLfdjj1K5UgJfLBiOclOBx67HW0fPW2/rvPC77/WSuVMJboaD60IITTgBeBYYDPwuxBihpRyRU3PrdRfmqZx25Tr+Gf5Jlb+uobM1hkMOO5gvnjj2/hdxJTkbsurdFbKO/e8z9SHP8Hu0BBCoDk0Hv7yDnr07xq/NtVzOX5f1CxSp6axIjubtimpzFi9ipvnzEYIq9RBLO8tW1r6WmkYSGJnqBYEg4QMo3ToRqld8XjXBwHrpJR/AwghpgKnAiqQNwGdDmhPpwPal/7erU+nuGW1SCCjZeVWElry3TI+eGwG4WC4zKIZt5/0IFP/fbW0/G1j1zEtPerkIN0wKQgFOXnKO6ys4qLLsaot7q2Z14szxrCOknjxGFppC2za6/fNJdvKEEJMEEIsEEIs2LFDrd7dWPUc2I3M1hlVfp3dUTYIuLwuTrvmRJxuZ4xXlPXZa18T8EXOHA0Fwiz7eVWV29NQjTu4L44oATVg6Nw4Z3aVg3hleOx2bjrsCIQQ+MNh3l++lNu//Yo3Fi8kL6BmitaGePTIow2gRXTHpJSvAq+ClX4Yh+sq9Ywe1nn84hcp2FH1KfNturdCDxlsWbeNpFQvo28cWaWl4PwxloJDELXWTGPVJSOTV0eO4qY5X5Lj91WqN10TTpuNCf0GsCk/n6d//YVpy/6kMBTCp4dxa3aemT+PD0aPoUccJispscWjR74ZaL/X7+2ALXE4r9LAvH33+8z9aD7hatRryWiRRu8hvXC4HAQDIT6f+A0Lv6p8AaejzjoMd5S6MkbY4MAj9qtyexqyw9t3ZO5Fl9K7ecuEXytkmjz/+3ye+20ez//+K9t9xfhKZn4GDJ3CUJAR773LdbM/Y0MFqZFKzcQjkP8OdBdCdBZCOIExwIw4nFdpYGa+/GWZGaBVsXTuKr6dMpdwMIwe0snesJPbT36Iaw6/nUXf7Lu64TFjDqfnwG64k92AtRScy+PkPy+Mx5vi2cerGx8hBJsL8mvlWoaUmBBzMpJumsxas4pTpr7LpvzaaVNTU+NALqXUgauBL4GVwPtSyuU1Pa/S8PgLq19UydTNiMqLUkpWzlvDXac+zIyXvqzw9XaHnUe+upOb3rqao8ccTpeDO+FOcvH6rZN58bo3Kc5vesua9WhWf4YzTMAXDvPigvl13ZRGKS555FLKz6WUPaSUXaWUD8TjnErD03Ngt4ScN+gL8epN7/Laze8yutV4RmVcwEPnPcPOf3PKHKfZNQ4fNZAt67axYfkm8ncWkrs9n5kvf8W1h9+BoSd2vLi+uf6Qw3DXo3RAQ0p++3dzXTejUVIzO5W4uerZi3EnudHsZf9aaVrN/5oFi4O8/9gM8rLzKc738f20X7hy4C0UF/jKHLf4m6VsWvVvmRREPaSTvXEn82Y2rfo+fVu34ZpBh8Zjsm3ctE1NresmNEoqkCtx03NAV15a9CjHX3QMPQZ05cjRhzLm5lFMeGwcbbu1wumOXnGxOkzDxFfg56u3viuzfd3if6KuQeovCrB20fq4Xb+haJWcjMcRv/e9JjQhCOsGH69ckfBsmqam/nzvUhqFdt1bc/0rl0dsP+GSoUx/ciY/vP9LSa9dY+2iv8uMi7u8ToK+yj8sDfqCLP9lNadds6ekaqtOzXG6HfiLygYKd5KL1l0Sn8VR3wxs0w6jgtmbtcmQkvlbNrM0eztv/bGIaWeejdtePz5kGjrVI1dqhTfFw7i7z+L15U/zwm8Pc+mj50UcU52Ml4Vf/UneDisTInd7Hq4kJ2a5WaVCCJxuJ0eNHly9xjdgbVNTGXvgwXj2Cphuu512KXU3xOHTw6zdlcP0FSonIl5Uj1ypE1Me/ChyfdBqTBPzFwV49cZ3sTvtzHnnB0zDxCy3pJzm0Hj067vwJDe9NESA2484msHt2jPpzz8oDoc4pUcvzjrgQEZOfZd1u3bVSZsCus7MNas476A+dXL9xkYFcqVWrZy/li/f/I7lc+Mzbd7QDb6f9gs2TcRcONru0Fi78G+6HtwpLtdsaIQQDO/SjeFdymYVTT79LMbP+JjVOTsxTBOjlqsXpjhjLAyuVJkK5Eqtee/hj5l8/3RCgXCVimoJIUBIZIyh3vA+St0GioOsnL+GEy4eWpXmNnrNvUl8OuY8/s7dxU6fj/M+/gC9lsbTXZqd7llZzNu0kUFt26FFKaerVJ5695RasXPLLibd+wFBXyhmEHe4HNiddpweq1CW3WWnc+8OzCh6l6/0Dzjk5H7YqpHK6PI4adejTY3a35h1ychkUNt29G+duPdIQ5DkcJJkd2ADdNNg0p9LmDDrU456a6Kavl9DKpArtWLRnD+x2aP/ddMcGq27tOTq5y7mvU0vM/6hsYy5eRT3fHQTLy95DLfH+gp+5TMX4XBW/Uuk3Wnn+AuPqVH7m4KHhh0Xc4WhmnJoNu444ijapKZgYmWwFIfDFIdDbCsuYsKsT6K+TjdN1u3KYeHWf5nz1zrW7cqJelxTp4ZWlFrhTnJhE5GB3KbZOO2ak7jssT2LLZ92zUkAhENhfvxgHsvnrSElM4lZL31FqIJVf4QmkMae3r5Ns9Hl4I6Mvf0MPp/4DXaHxpFnDqZFh8otVNHUdErP4McLL+W6Lz/j182bqvPsOaagYXD3D99GzR83peSfvFzW5uyk+15VEj9fu5rbv/2aolAQQ0o0IbDbbAxo05ZXR4yqN/nx9YGoi+WZBgwYIBcsaFqz7Jq6gC/IWa0vjSg36/Q4eX7+Q3Tu3aHM9uICH9ccejvZm3YSKApgs4mItMK9OdwOhowaxO+zl+Ar8OEomXyUlOqlKK8YQzexaTaEgGtevFT10GOY89c6rv3ys5iLMydS/9ZteP/MMQgh+HP7NsZ8OC1qO5yaxpn7HcD9Q4+t9TbWNSHEQinlgPLb1dCKUivcXhf3zbgZb6rH+knx4HQ7uOLJCyKCOMDk+z9k69/bCBRZhbgqCuJOj4PDThlIh/3bEQ6GMU1J0Bci6Auxa1seoUAYQzcIB8OEAmGevfI1crNVFb5oHp83t06COMAf27fx62ZrjZo3Fi8kGKMdIcPgw5XL1Rqhe1FDK0qtOfioA/hg20QWzvmTkD9E32EHkpqVEvXY76f9TDhYiYAi4JIHx3Lq1ScwtuOVlZpUZNNszJ+1UGWxRFFbpW+j0U2TO76bwzunjebfwoIKh3ZChoFZMtyiqECu1DKn28mhIyO+GUYov/RbrGOOHH0op19rTdEPRlnqTamazhmZrNiRHbHdpWm0TUllfV5uXMfOy1ufl8dRb02kT8tWODUtZk2Wvq3bqJTFvah3QqmXTrhkGC5PlPU6hVWu1pXkpE23Vlz17MWluwae2KdS6YmmYTJ4ZH9ytuby0TOfMeXBj1i3uOkV1CovL+BnYOu22MsFyN1rch7ZsTNJtVAbxZSSRdu2xgziSQ4n9x09LOHtaEhUj1ypl868YSRLvl3Gyl/XYOgmmsOGHtRBCIRNYIRNOh/YkaRUb+lrxj98Hou/XkpxoZ9wlAqIYKU6XvPSpSz9aRWPnP8sUkr0sMGUBz9k+PlHcu2LE6wJSE3M3I0buGzWp4DElLK09G3blFT+c8ihvPPHYtbm5BAy675qod1mY9qKZVx3yKGku5tm2YXyVI9cqZecLgePfHUnj359NxMeO59ufTojbNY0/HDAWg7u15kLmPbonvzjFu2b8caqZzj6rMMiaqLv1nfogRx5xmAeGfccQX+IUCCMaZgEfSG+mfQTiyuxrFxjE9R1rvp8Bn49jF/XMaVEAh6Hg7uPGkqrpGT+ycutF0EcID8Y4L2lf3D6tCkxH4g2NSqQK/WWEIL9DunOyROGs/q3dREPP4P+EDNe/KrMtpSMZIafdyQuT2QdDyEgo2UaC+f8GTXQB4qDfD3px/jeRAPw+5Z/o457+8JhPli5jBU7s2MGzM5p6YltXAxh0yTbV8wX69bWyfXrGxXIlXpPD+kx0w/9RX7WLvqbP75fjr/YSlU86Kj9sbsiRw2dHhcnXTq8wqGTpjisUlEan5SS9qnpuKIsGed1OLhi4CF8OPqcOske8YXDLNm2pdavWx+pQK7Ue55kD+17RdYBETaBaZjccPTd3DXqEUa3HM8Xb3yD3WHnwc9vJyUzGW+qB0+KG4fbwdg7zqD34b3oN/xATD2yOJQ7ycWw846sjVuqVwa1bRc1mHsdDk7f7wCGd+lKstNVJljbhMBtt3NS9570bd2GB4ceiz3KzN1EcmkandIzavWa9ZUK5EqDcP2rl+NOcmF3WD3D3cvGBYqDBIoC+Ar8BH1BXrjmDVYv+IueA7oybcur3DblOq5/5XImr3+Rc245DbA+GG586yocLgdOt8PKgvE4Oe6Co+k7tHed3WNdcdntPHviCNx2Oy5NwwZ47A6Gde7KsV264dQ0PjzrHAa1bYfdZsNus9G/dRs+HH0u3pJp8mfu35tTevaq1fVBbUIwqtd+tXjF+ktN0VcajK3rt/Pp81+wYcVmmrXN5PtpvxAoLps7LmyC48Ydxf+9cVXUc0gpef/xGUy+bzqmKdHDOt37deGaF8bTvV+X2riNemtHcTGz1q6mMBjkiA4d6dOqdcRQkz8cRkJpAN/bLr+Pwa+/UoulcDX2b96C6wcfzpAOHWvlmnUt1hR9lX6oNBitO7fk8icuBGDezAX8OP3XiGOkKWNOv5dS8tzVE5n1ypwypXTXL93AN5N/bPKBvHlSEhf16VfhMRUVqioKhXDYbLUWyIOGweJtW5kw6xOeOf4k8oNB3lyyiMJQkGM7d+OKgYfQzOvd94kaATW0ojRIBxzWM+qKQC6vi0NHDoz6mom3TGLWy3Mi6qEHfSFmvTIHPaxS2WqiXWoayc4ok7gSLKDr/HfObO7+/htW7tzB5oICJi1dwsj33qUg2DRm+6pArjRIqVkpnHfnmbi9e9IMXR4nbbq25NhxkQ8sc7fn8fFzX8TM0DB0E39JgS6lemxC8PDw43FHyXBJtKJQCP9eKZJh0yQ/EOC9ZX/WelvqggrkSoN1zq2nc88nN3HYqQM58Ij9uOThsTw778GoOeRrFvyFs4JFKdKyUkhOT0pkc5uEYzp14ZOzx5JRD2ZcBgydXzZtqOtm1Ao1Rq40aP2GH0S/4Qft87isNpnoUVIOdxtyxuAmmUOeCD2ymvHCSSO5eMZHdVYSF0ATgg51NGGptqkeudIkdO3TifTmqTH3Z2/cUYutafwGt2vPW6eeQYfUtDprg1PTuODgvnV2/dqkeuRKo5S7PY/v3vuZ/JxCOu7fjrkfzScvRjaLEJCU1jSyG2rToLbtOLxDRzbWwTh1c6+XR4afQLfMrFq/dl1QgVxpdBZ9s5S7T30E0zQJxaiCuDenx8XJlw6vhZY1PQc0b4HH7sCv7/vPIV5aJScz96IJ2JrQUJkaWlEaFUM3uP/sJwn4gvsM4q4kFw6Xg/PvOpPeQ9QMwUQ4ted+4KniuAAADexJREFUpDidEbVYEhViPXY7l/cf1KSCOKhArjQyq35bh6Hvu9yqy+vk1CuOZ/KGlzj7plG10LKmKcnp5JMxYzm+a3c8djspTifn9j6I508cSYbbHffrtU1J5aiOneN+3vpODa0ojYpNs1V6Ud7h444io0XdPYxrKlolp/D8SSMjtndIT+e0qZPQ41gm5K/cXYyc+i5PHHsCS7O3Y7dpjOjRky4ZmXG7Rn2keuRKo9JjQJeoeeR7czjtdO/Xhc69O9RSq5RoDmjegtdPPZ2WSfHL35dYk4Ou+HwmL/7+G8/9No+Tp7zLm0sWxe0a9ZEK5Eqjomka93x8I54Ut1Ut0anhcNlJSvei2TUcTjtDzhjM/bNureumKsCh7TpwWPuOOOK8kLIpJSYSQ0qChs6DP33PpvzoWUuNgap+qDRKxQU+5n40n4KcIvoccwDd+3XBX+TH7rTjcCZ+AWGlcp7/bR4vLvitViYO9cjMYvZ5Fyb8OomUkOqHQojHgJFACPgLuEhKmVeTcypKPCSlejn+wmPKbPMk1/20caWsd/5cUmuzP9fn5bI8ezsHtGhZK9erTTX9PjMH6C2lPAhYA6jvq4qiVFpRKBRzX7eMTJ4/cUTcKioapmTB1sa5NFyNArmU8isp5e6P01+BdjVvkqIoTcWgNu2i5pT3ymrGV+dfxEndezJ99Dl4HQ60kqXknNUcTzeRzF63Jub+Hb5iXl+0gEd//omfNv6DWQfDztUVtzFyIcRMYJqUclKM/ROACQAdOnTov2FD06hKpihKbOt25XD6+1MI6jph00QTAqem8e5po+nXes86rTuKi5m0dAnLs7PpmpnJxEULqG7k+u6CS0hyONlSWECn9Iz/b+/eg6MqzziOf5/dPbubkEAwARISJIEit3AnkeKIQoSiMuAdEAtWq4WqU6ftVFtmnHZaZ9pqbTvVYovSOq2jQkVAoVpoUSwMyFW5JdzKJYR7gATIZXfz9o9dKWl2A0lOcnbh+fyV3ezZffJO8tvNe97zvLT3+Vhz6CCPvf8edcZQEwqRbFkMzezK6xPvxnK77flhbRBrjvyyQS4iK4DMKN+abYxZHHnMbGA4cI+5gncGPdmplPpSWWUF8zZv5PNjR+mdnsGjQ4eT18imyhU11RTMnUOgmTsRZbZLoby6Cp/bTW0oxEMDBrOweDunq+v3o0/yeHhu1Ggm51++u2ZbaXaQX8ETzwBmAkXGmAtXcowGuVKqJaa++w4byg4TsmFGwed2Y4DaUMMrgodndWX+/VNb/Bp2iRXkLZojF5HxwDPAxCsNcaWUaqmXxt1BZkoqPhumPWpCoaghDuC2eX17a2lplS8DqcByEdkiIq/aUJNSSjUqKzWVj2c8ys9G34anlcLW53Yzuf+AVnluu7VoHbkx5it2FaKUUrEE6+pYtruExSU78bk9TMkfyM3Xd+fefvmsPniQRbt22v6a7X0+JvZOjK6Y2jRLKRXX6ozhsfcXsb6slAuBcGviTw7s58H8gcwedStPjxjJ0j0lzT75GYsgCdMONzEmgJRS16xVB/bXC3GAqmCAv27dwoEzZ+jWoQNZKam2v24idUzUIFdKxbWV+/fVC/EvuURYU3oQEeG3t0+gneXF7wlPMrSzLPpldGLrzKeY2r/pywctl4unR4xsce1tRadWlFJxLc3nx3K5GkyduERIjVy+P6hLJp88/CiLindSVllBQXYORXk98bhczCoo5K3tTds39Ffj7qAwO3EuVNcgV0rFtXv69mfu5g0NglxEKMrrefH2dUnJPDJkWIPjs1JS8UYu/rkSHpeLFfv2UJTXgyQrMTpl6tSKUiqudU9L48Wx40n2WKR4vaRYXtL8fv486d4rClq3y8XU/IFX3KMlWFfHkl3FDJv7ew6eTYxmrtqPXCmVEKoCAdaXHcZyuSjIzsHjchEIhViyq5ilu0pI9XmZmj+IETndGhxbGwrxrQ8W8cmB/U16zb4ZnVj64PSLtzcdKWPe5o0cOVfJqO65zBg0hDR/27VHbpV+5Eop1VaSLItR3XMv3g6EQkx7bwHbjx+nKhg+Gbpi315mDS/kycKv1jvW63Yzc1ghG8oOcz7KidNY9pSf4vj5c3Rul8K7O7bx3Mf/pDoYxAA7Thzn7W1fsPTB6VyXlGzHj9hsOrWilEpIH+3dzY4T/wtxgKpgkFfWr+PE+fMNHj8oM5OmTkAE6+o4X1tLTTDIT1atpCoS4hC+tL+8qpq5m5yfXdAgV0olpOX79kZdluhxuVh3+FCD+/0ei+eLxjbpNQywpKSYvafLo74JBOpC/Os/+5r0nK1Bg1wplZA6+v0xrrwUUr2+qMdM6t2XwV2ideWObWHxDjr4/QTroq96uS7J+S0ENciVUglpcv5AvFG6H3rdbkZ2uz7mcc8XjSPZsvBE3gTcIvg9HtwxLscPmTqyU9uT37nLxWO+lOSx+OaQBuceAThTXcVrmzbw3Y+W8actG6moqY76ODtokCulElLfjE78+JYx+D0eUiPLEjOSk/nL3fc1uqtP34xOLHtwOg/0H8CAzl24u08/lkx5iJ5RLsn3ud3cFWmcNefOSfTr3IWkyOv53G6eLLiRoh49Gxy373Q5o9+Yx0trV7OoZCcvrvk3o9+Yx6GzZ+0bgEvo8kOlVEI7V1vLhrLDJFsWw7K6NruHePHJE0x59x0CoRBVwSDJlkVuhzTeuW8KSZbFrlMnERHcCKeqLtC3U2fa+6JP4UxbOJ+1pYfqbUfnEuHW7nm8NvHuZtUHuvxQKXWVSvF6uTU3r8XP0yejE58+/Bgf7C6htOIsQzKzGJ3bgy+OHeXby96nsrYGgDS/nzl3TooZ4nXGsO5waYM9ReuM4dOD+1tcZzQa5EopFZHq8zH1kj06z1RXMX3R3+qtPb8QCPDQwgWsfuRxUiK9Xi4lhFfORGsJEG1O3w46R66UUjF8sKsk6r6gIVPHh3t2RT1GRJjQq3eD0Pa63dzVp1+r1KlBrpRSMZy8cJ7qYLDB/bXBICcvxN6m+LlbxtAnPYNkyyLZskjyWOR37sKzN41qlTp1akUppWIoyM4h2bIaXHjk9XgoyM6OeVx7n4/3Jk9jy9Ej7D1dTq/0DAZ27oK00o5DGuRKKRXDyJzrGdwli81Hy6iKfDJP8ni4MTuHoZldGz1WRBiS1ZXBmVl8uHc3z69ayenqaoryevD4sAJb+7NokCulVAwiwrxJ9zB/+1YW7NiGAA/0H8AD/QfU+3S9p/wUv/tsLVuOHiE3LY0nCkZc3JjiN+vW8NqmjRd7whw6e5bFJcV8OG0GHfx+e+rUdeRKKdV8O0+e4P4Fb1EdDFIXydMvrxbt4Pdz7Nx5Qqb+phjhi4lG8EThiCa9Vqx15HqyUymlWuAXq1dxIRC4GOIAIWM4HwhQVlnZIMQh3DlxlY1ryjXIlVKqBbYcPdLkY1widE1tb1sNGuRKKdUC6c04ael1u/nG4KG21aBBrpRSLTBzeCFJnsb3DhUEy+WinWXR3ufjF0VfY2AT2+k2RletKKVUC9zXtz/HzlUyZ8NniEiDNeeWy8UN6Rn8YcIkKmpq6Nnxuka7MzaHrlpRSikbVAUClFZUcLiygp+uWklpRQVguPn6XF4YO56ONmxAod0PlVKqFSVZFr3S0+mVns4t3XMpr6rC7/HQLkpjLbtpkCullM1EhPRk+67cvBw92amUUglOg1wppRKcBrlSSiU4DXKllEpwGuRKKZXgHFlHLiIngAMxvp0BnGzDchKJjk10Oi7R6bhEl8jj0t0Y0+n/73QkyBsjIhuiLXhXOjax6LhEp+MS3dU4Ljq1opRSCU6DXCmlElw8BvkfnS4gjunYRKfjEp2OS3RX3bjE3Ry5UkqpponHT+RKKaWaQINcKaUSXNwGuYg8JSIlIrJdRH7pdD3xRES+LyJGRDKcriUeiMgLIlIsIl+IyHsikuZ0TU4SkfGRv509IvKs0/XECxHpJiIrRWRnJFe+43RNdonLIBeR0cAkYKAxpj/wosMlxQ0R6QaMBQ46XUscWQ7kG2MGAruAHzpcj2NExA28AtwO9AOmikg/Z6uKG0Hge8aYvsAI4ImrZWziMsiBWcDPjTE1AMaY4w7XE09+DfwA0LPUEcaYfxhjgpGba4EcJ+txWCGwxxizzxhTC7xN+EPRNc8Yc8QYsynydSWwE8h2tip7xGuQ3wDcLCLrROQTESlwuqB4ICITgcPGmM+driWOPQL83ekiHJQNHLrkdilXSVjZSURygSHAOmcrsYdjOwSJyAog2jbSswnX1ZHwvz8FwHwR6WGugbWSlxmXHwHj2rai+NDYuBhjFkceM5vwv89vtmVtcUai3HfV/900hYikAO8CTxtjKpyuxw6OBbkx5rZY3xORWcDCSHB/JiJ1hBvdnGir+pwSa1xEZACQB3wuIhCePtgkIoXGmKNtWKIjGvt9ARCRGcAEoOhaeMNvRCnQ7ZLbOUCZQ7XEHRGxCIf4m8aYhU7XY5d4nVpZBIwBEJEbAC+J263MFsaYrcaYzsaYXGNMLuE/2KHXQohfjoiMB54BJhpjLjhdj8PWA71EJE9EvMAUYInDNcUFCX8Ceh3YaYx5yel67BSvQT4P6CEi2wifrJlxjX/KUo17GUgFlovIFhF51emCnBI56fsk8BHhk3nzjTHbna0qbtwEfB0YE/k92SIidzhdlB30En2llEpw8fqJXCml1BXSIFdKqQSnQa6UUglOg1wppRKcBrlSSiU4DXKllEpwGuRKKZXg/guvhM5nkGDafgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "X, y = make_blobs(random_state=170, n_samples=600)\n",
    "rng = np.random.RandomState(74)\n",
    "\n",
    "transformation = rng.normal(size=(2, 2))\n",
    "X = np.dot(X, transformation)\n",
    "\n",
    "y_pred = KMeans(n_clusters=3).fit_predict(X)\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1], c=y_pred)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Some Notable Clustering Routines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The following are two well-known clustering algorithms. \n",
    "\n",
    "- `sklearn.cluster.KMeans`: <br/>\n",
    "    The simplest, yet effective clustering algorithm. Needs to be provided with the\n",
    "    number of clusters in advance, and assumes that the data is normalized as input\n",
    "    (but use a PCA model as preprocessor).\n",
    "- `sklearn.cluster.MeanShift`: <br/>\n",
    "    Can find better looking clusters than KMeans but is not scalable to high number of samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- `sklearn.cluster.DBSCAN`: <br/>\n",
    "    Can detect irregularly shaped clusters based on density, i.e. sparse regions in\n",
    "    the input space are likely to become inter-cluster boundaries. Can also detect\n",
    "    outliers (samples that are not part of a cluster).\n",
    "- `sklearn.cluster.AffinityPropagation`: <br/>\n",
    "    Clustering algorithm based on message passing between data points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- `sklearn.cluster.SpectralClustering`: <br/>\n",
    "    KMeans applied to a projection of the normalized graph Laplacian: finds\n",
    "    normalized graph cuts if the affinity matrix is interpreted as an adjacency matrix of a graph.\n",
    "- `sklearn.cluster.Ward`: <br/>\n",
    "    Ward implements hierarchical clustering based on the Ward algorithm,\n",
    "    a variance-minimizing approach. At each step, it minimizes the sum of\n",
    "    squared differences within all clusters (inertia criterion).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Of these, Ward, SpectralClustering, DBSCAN and Affinity propagation can also work with precomputed similarity matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"images/cluster_comparison.png\" width=\"80%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For more see**: [Clustering Algorithm Comparison](https://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise: digits clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Perform K-means clustering on the digits data, searching for ten clusters.\n",
    "Visualize the cluster centers as images (i.e. reshape each to 8x8 and use\n",
    "``plt.imshow``)  Do the clusters seem to be correlated with particular digits? What is the ``adjusted_rand_score``?\n",
    "\n",
    "Visualize the projected digits as in the last notebook, but this time use the\n",
    "cluster labels as the color.  What do you notice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "digits = load_digits()\n",
    "# ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/08B_digits_clustering.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7 (DL-TORCH)",
   "language": "python",
   "name": "dl-torch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}