010_cosine_elbow_method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the elbow method to find the optimal number of clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFV5JREFUeJzt3W9sXGeVx/HfyWQKdmHrRLW06jRpSoVSmhpiZFFDJCRaRGAhwRuhLRHhxa60fbP8cYiMmt0uSaQsQQpiwwuEFJVFWiUCNyFrpSxLulLZF1ttKpzaWWNCpAKtm0kRZpspqB3IdHz2xXhSxxnHY/veuffO8/1IlWJn5urYTX5+cu5zn2PuLgBAdqxKugAAwNIQ3ACQMQQ3AGQMwQ0AGUNwA0DGENwAkDEENwBkDMENABlDcANAxqyO46K33367b9iwIY5LA0BbOnfu3O/cvbuZ18YS3Bs2bNDo6GgclwaAtmRmLzb7WlolAJAxBDcAZAzBDQAZQ3ADQMYQ3ACQMbHsKgGArBsZK+rwmYu6XCrrjq4ODW3dqIHeQtJlSSK4AeAGI2NF7TlxXtWZ2oSwYqmsPSfOX/v9pAPd4hhd1tfX5+zjBpAlI2NFHXhyUlderyz5vV0dee3fvmlFAW5m59y9r6nXEtwAQjcyVtTQyfOqVJefhybJJRWWuQpfSnBzcxJA8A6fubii0JZqoS3V2ip7T01oZKy48sIWQHADCN7lUjnS65UrVR0+czHSa85FcAMIXldnPvJrFiP+YTAXwQ0geDHc6lPOLPqLziK4AQTv1fLSd5IsphrHT4NZBDeA4N3R1RH5NQsxXLOO4AYQvKGtG9WRz0V2vY58TkNbN0Z2vfl4chJA8Op7rgeHx1d8rTWdee3btrKHcRZDcAOAauH9D/82odeuVpf1/uU+eLMctEoAYNY//WWPlrMXZFf/ej3z6IMtO7OEFTcAzKoH7+7hcTWzJyRnpp0PrNPBgZ54C5uH4AaAOerhvffUhMqVG9smnflV+uqOdyd6xCvBDQDz1EM56eNbF0JwA0ADA72F1AT1fNycBICMIbgBIGNolQBIrTTPfUwSwQ0glUbGitft7KgPKJAUfHgzugxAKsxfXV957U96vTKz4OujmPOYJksZXcaKG0DiGq2uF1MqVzQ0O3m9XcK7WdycBJC4w2cuNnzYZTGVGdf+05MxVJRuTQW3me02s0kz+5mZfc/M3hp3YQDCsZKZj6VyRXc/+u96bGQiworSbdHgNrOCpC9I6nP3+yXlJH067sIAhGOlMx9d0rGzU9r0lR/HOl09LZptlayW1GFmqyV1SrocX0kAQhPVHonXrla199RE24f3osHt7kVJX5c0JellSa+6+1NxFwYgHFHOfCxXqjp85mJk10ujZlolayR9UtLdku6QdKuZ7WrwukfMbNTMRqenp6OvFEDbinrm40p65lnQTKvkw5J+7e7T7l6RdErSB+a/yN2Punufu/d1d3dHXSeANhb1zMc4hv+mSTPBPSWp38w6zcwkPSTpQrxlAQjJQG9Bh3ZEM4wg7kG9adBMj/tZSSclPSdpYvY9R2OuC0BgBnoL6upY2e6Sro68Du3oafsHcpraVeLu+9z9Xne/390/6+5/irswAOHZv33Tst+7q3+9xvd9pO1DW+LJSQApMtBb0K7+9Ut+35GHN7d87mOSCG4AqXJwoEdHHt6sW3LNzVvf1b8+iFX2XBwyBSB15o8Nm38IVd2We9YGtdKuI7gBpF7ah/e2GsENIBPSPLy31ehxA0DGsOIGsCLMhWw9ghvAsjEXMhkEN4AFLbaabjS5pn46H8EdH4IbwHXqYT1/7mOxVNbu4XENDo8rZ6bqTQ7RbvfT+ZJGcAOBm7uqvq0jrz/86Q1VZxqHcv2zNwttqf1P50sawQ0EbH6PuhTBQIMQTudLGtsBgYAtd7r6QkwK4nS+pBHcQMCi7kVHNDoSiyC4gYDdtsLzrxtp93mPaUBwAwGz5g7gWxJ2lMSP4AYCduX16Kar17GjJH4EN4DI5FYZO0pagOAGEJm3v2U1O0pagOAGApaLuMn9agT7wLE4ghsI2M4H1kV6PfrbrUFwAwE7ONCjXf3rr628c2ba1b9eL3zt40tejfPEZOvwyDsQuIMDPQ3nNu58YJ2OnZ1a8H1b7lmrF/6vzDncCSC4ATRUD/PvPfuSqu7KmWnnA+uCHM6bNuaLnPK1HH19fT46Ohr5dQGgXZnZOXfva+a19LgBIGMIbgDIGHrcQIoxiBeNENxASjGIFwuhVQKk1EKDeAeHx/XYyERCVSENWHEDKTIyVtSBJycXPbXv2NkpHTs7pQLtkyCx4gZS4rGRCQ0Ojy/pqNViqazB4XG96x//QyNjxRirQ5oQ3EAKjIwVdfwmTykuplyZ0dCJ84R3IJpqlZhZl6THJd2v2li5v3H3/4mzMCAE9V0jxQimxlRmXF96YlwSNy/bXbM97m9K+rG7f8rMbpHUGWNNQBDm7xqJwoxLQyfPSyK829mirRIz+zNJH5T0HUly96vuXoq7MKDdNdo1EoVK1RnY2+aa6XG/Q9K0pO+a2ZiZPW5mt85/kZk9YmajZjY6PT0deaFAu4miPbIQBva2t2aCe7Wk90r6trv3SnpN0qPzX+TuR929z937uru7Iy4TaD9RT5+Z67aOfGzXRvKaCe5Lki65+7OzH59ULcgBrEA1hpM562L8mYAUWDS43f03kl4ys/poi4ck/TzWqoAAFGIc81Vawl5wZE+z+7g/L+m4mf2vpM2SvhpfSUAYhrZuVEc+F8u1mf3Y3poKbncfn+1fv9vdB9z9StyFAe1uoLegQzt6VOjqkKm2Aj/y8OYVX5fZj+2Ps0qABA30Fm7Ybz04PL7s63F2SRgIbiBl1nTml3ReiVRbZR/a0UNgB4KzSoCU2bdtk/K567eF5HOmXf3rtabzzW1+9VcUujoI7cCw4gZSph7AjSbfMGEdEsENpFKj3jdQR6sEADKGFTcwDwN6kXYENzAHA3qRBQQ3MMdCA3r3n568bhX+oXu79ZNfTKtYKitnpqo7e6jRMgQ3MMdCx6GWyhWVyrW91cVSWcfmjBmrHxbF6hytws1JYI6VnvFRrlQZYoDYEdzAHB+6d+VnyTPEAHEjuIFZj41MXNcCWS5O5kPcCG5Atd0kUYS2KZpVO3AzBDcgRdaXdkk/OFfUyFgxkusBjRDcgKLtS3ODEnEjuAFF35fmBiXiRHADUuQTY7hBiTgR3ICifWCG0WGIG09OAhHisXe0AsENzNpyz1o988tXlv3+NZ15PfPogxFWBDRGqwSYdfxv37/s966y2sgxoBUIbmCOwk1uKq7pzKszf+Nfma6OvL7xV5tpj6BlaJUAcwxt3XjdedwSE9SRPgQ3MMfNBvUCaUFwA/MwqBdpR48bADKGFTdSj+G9wPVYcSPV6sN7i6WyXLXxYIPD49p84ClO4EOwCG6kWqPhvVJtBuTeUxOEN4JEqwSJa9QKkWqhXbzJKXvlSlV7njgvieG8CAvBjUTVWyH1VXW9FdKsqruGThDeCEvTrRIzy5nZmJn9MM6CEJaFWiFLUZlxDQ6Pa8vXnqZ1giAspcf9RUkX4ioEYYpy4ECxVNbu4XE9NjIR2TWBNGoquM3sTkkfl/R4vOUgNFEPHHBJx85OsfJGW2t2xX1E0pclzcRYCwI0tHWjLIbrHnhyMoarAumwaHCb2Sck/dbdzy3yukfMbNTMRqenpyMrEO1toLcgj+G6V16vxHBVIB2aWXFvkbTdzF6Q9H1JD5rZsfkvcvej7t7n7n3d3d0Rl4l2lrM41txA+1o0uN19r7vf6e4bJH1a0tPuviv2yhCMqke/5u7qyEd+TSAteHISibvZ8ILl2r+daTRoX0sKbnf/L3f/RFzFIExDWzeqI5+L7Hq7+tfzMA7aGk9OInH1kN3zxPkVtU3WdOa1b9smQhttj+BGKtTDdvfw+JJ2mezqX6+DAz3xFAWkFD1upMZAb0Gf6V/f9L5uQhuhYsWNVDk40KO+u9ZeOy3wto68zGr7snNmqrqrwDAFBI7gRuow8xG4OVolAJAxBDcAZAzBDQAZQ48bS8bUdSBZBDeaNjJW1IEnJ687ea8+vGD0xVfYmge0CMENSbVQ3n96UqVyLZRvvSWnfG6VXi1X1NWZ1x8rVZUrjY9jd0nHz06p7661rLyBFqDHDT02MqHB4fFroS1Jr12tqlSuyFXbQ71QaNe5avMjAcSP4A7cyFhRx85ORXKtYoTzIwEsjOAOXNSr5M0HnmLeIxAzgjtwUU5Zl6RSuaLB4XECHIgRwR24qKes15XKFe09NUF4AzEguAM3tHVjbNcuV6rcsARiQHAHLu7te1G3YgAQ3FA8Mx/r4mrFACEjuBH5zMe6jnwu1lYMECqCGxroLejQjh7lrNnZMzfacs9aHXl4swpdHTLVVvGHdvTwJCUQAx55h6Q3e92Dw+NLet/8Ab0ENRA/Vty4ZqC3oDWd+aZfb5LGvvIRwhpoMYIb19m3bZPyueZaJtx4BJJBcOM6A70FHf7Ue65beXfkV90Q5tx4BJJDjxs3aDSsl+EJQHoQ3GgKk9eB9KBVAgAZw4o7heptiWKprJyZqu4q0J4AMIvgTpmRsaL2nppQuVKVJFXdJdWGFAydOC+p1rag5wyEi+BOmcNnLl4L7fkqM67B4XGdGJ3Sc1OvXntdsVTW3lMTkngABggBPe6UaWb81zO/fOWGcOcIVSAcBHfKLP+0EI5QBUKxaKvEzNZJ+ldJfy5pRtJRd/9m3IWFZO7NyJXgSUYgDM30uN+QtMfdnzOzt0s6Z2b/6e4/j7m2IMy/Gblc+ZzxJCMQiEVbJe7+srs/N/vrP0i6IIk7YBE58OTkikNbkt6Y8QiqAZAFS+pxm9kGSb2Sno2jmNCMjBV15fVKJNdyF8N5gUA0Hdxm9jZJP5A06O6/b/D7j5jZqJmNTk9PR1lj24p6Fwg7S4AwNBXcZpZXLbSPu/upRq9x96Pu3ufufd3d3VHW2Lbi2AXCzhKg/S0a3GZmkr4j6YK7fyP+ksIRxy4QdpYA7a+ZFfcWSZ+V9KCZjc/+9xcx1xWEqIf0ckY2EIZFtwO6+39rZc+FYAH1x9Oj2MPNIVRAODirJGFzz7keGSvqwJOTS9ppsqt/vQ4O9MRVHoAU4pH3FBnoLWjsKx/RkYc3q9BEr5rQBsLEijuFFluFd3XktX/7JtoiQKAI7pRjZBiA+WiVAEDGsOJeABNmAKQVwd3A/BP7mDADIE1olTTQaHwY54AASAuCu4GFHoYplsra8rWnOYEPQKJolTSQM7s2XX2+Yqms3cPjGhwev/a0oiT64QBahuCeo35DcqHQrqv/brFU1tCJ86q6qz7HoP45iX44gHjQKplVvyG51DNDKjNvhvbcz+0/PRlhdQDwJlbcqoX2nifOL7rSXopSOZrJNgAwX/Ar7vpKO8rQBoA4BR/cjbb+RWFNZz7yawKARHCv+BzshezbtimW6wJA8MGds3hmRLCjBEBcgg9uetsAsib44G5mYMFSdXXQ3wYQn+CDO+qBvflVpv3b6W8DiE/w+7ijGthrEo+7A2iJ4INbenPKzPzjXJvF7EcArURwz1FfKTc7aT1npp0PrCO0AbQUwT3P3NV3vX1SPy2wQCsEQAoQ3AtgSC+AtAp+VwkAZA3BDQAZk+pWychYUftPT147InVNZ177tm2ihQEgaKkN7pGxor40PK6ZOZ+78npFQyeZLgMgbKltlew/PXldaNdVqs60dQBBS21w32yCTLFUZtI6gGCZx3A6Xl9fn4+Oji7pPSNjxaYffJGkVSa9ZfUqlSu1dTn9bwBZZmbn3L2vmdemosc9MlbU0MnzqlSb/yEy47oW2hL9bwDhaKpVYmYfNbOLZva8mT0adRGHz1xcUmgvhP43gBAsGtxmlpP0LUkfk3SfpJ1mdl+URVyOcHxYlNcCgDRqZsX9PknPu/uv3P2qpO9L+mSURdwR4TCDVWbcuATQ1poJ7oKkl+Z8fGn2c9cxs0fMbNTMRqenp5dUxNDWjUt6/c1U3bX31AThDaBtNRPcjabp3tCQdvej7t7n7n3d3d1LKiLqm4nlSpVeN4C21UxwX5K0bs7Hd0q6HHUhnflot5TT6wbQrppJy59KeqeZ3W1mt0j6tKTTURfy1R3vjvR6UfbNASBNFg1ud39D0ucknZF0QdIT7j4ZdSFRtktM0fbNASBNmnoAx91/JOlHMdeiNZ35pp+cXIhJ+kz/eh7CAdC2UnVWyb5tm5b1vvrd00JXh/754c3MgATQ1lLxyHvdQG9Boy++ouNnp27ctrIA5kACCE2qgluSDg70qO+utTp85qIul8rq6szrj5XqdeeSSFJHPqdDO3oIbADBSV1wS40H9danrl8ulXUHq2wAAUtlcDfC1HUAqEnVzUkAwOIIbgDIGIIbADKG4AaAjCG4ASBjYhkWbGbTkl6M+LK3S/pdxNfMktC/fonvQehfv9Te34O73L2pM7FjCe44mNlosxOQ21HoX7/E9yD0r1/ie1BHqwQAMobgBoCMyVJwH026gISF/vVLfA9C//olvgeSMtTjBgDUZGnFDQBQBoLbzD5qZhfN7HkzezTpelrNzNaZ2U/M7IKZTZrZF5OuKQlmljOzMTP7YdK1JMHMuszspJn9YvbPwvuTrqmVzGz37J//n5nZ98zsrUnXlKRUB7eZ5SR9S9LHJN0naaeZ3ZdsVS33hqQ97v4uSf2S/i7A74EkfVG1maeh+qakH7v7vZLeo4C+F2ZWkPQFSX3ufr+knGpDy4OV6uCW9D5Jz7v7r9z9qqTvS/pkwjW1lLu/7O7Pzf76D6r9hQ3qfFszu1PSxyU9nnQtSTCzP5P0QUnfkSR3v+rupWSrarnVkjrMbLWkTkmXE64nUWkP7oKkl+Z8fEmBhdZcZrZBUq+kZ5OtpOWOSPqypJnFXtim3iFpWtJ3Z9tFj5vZrUkX1SruXpT0dUlTkl6W9Kq7P5VsVclKe3Bbg88FuQ3GzN4m6QeSBt3990nX0ypm9glJv3X3c0nXkqDVkt4r6dvu3ivpNUnB3O8xszWq/Uv7bkl3SLrVzHYlW1Wy0h7clyStm/PxnQrwn0hmllcttI+7+6mk62mxLZK2m9kLqrXKHjSzY8mW1HKXJF1y9/q/tE6qFuSh+LCkX7v7tLtXJJ2S9IGEa0pU2oP7p5LeaWZ3m9ktqt2QOJ1wTS1lZqZab/OCu38j6Xpazd33uvud7r5Btf//T7t7UKstd/+NpJfMbOPspx6S9PMES2q1KUn9ZtY5+/fhIQV0c7aRVM+cdPc3zOxzks6odif5X9x9MuGyWm2LpM9KmjCz8dnP/b27/yjBmtB6n5d0fHYB8ytJf51wPS3j7s+a2UlJz6m2y2pMgT9ByZOTAJAxaW+VAADmIbgBIGMIbgDIGIIbADKG4AaAjCG4ASBjCG4AyBiCGwAy5v8BqhkFesF+nq8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cluster import KMeans\n",
    "%matplotlib inline\n",
    "\n",
    "# number of clusters\n",
    "true_n_cluster = 10\n",
    "test_n_cluster = 20\n",
    "\n",
    "n_obs = 100\n",
    "\n",
    "X = None\n",
    "\n",
    "for i in range(true_n_cluster):\n",
    "    x = np.random.normal(loc=i, size=n_obs, scale=0.1)\n",
    "    y = np.random.normal(loc=i, size=n_obs, scale=0.1)\n",
    "    data = np.column_stack([x, y])\n",
    "    if X is None:\n",
    "        X = data\n",
    "    else:\n",
    "        X = np.row_stack([X, data])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(X[:, 0], X[:, 1], 'o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## computing within-cluster sum of squares (WCSS) measure for cluster numbers from 1 to cluster_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.65221814e+04, 4.02639090e+03, 1.81632322e+03, 1.01655314e+03,\n",
       "       5.14422635e+02, 4.13735851e+02, 3.14827936e+02, 2.15370088e+02,\n",
       "       1.16462173e+02, 1.95783668e+01, 1.87294460e+01, 1.79874591e+01,\n",
       "       1.72627946e+01, 1.64961841e+01, 1.58364681e+01, 1.50798745e+01,\n",
       "       1.46037751e+01, 1.39168622e+01, 1.32227080e+01, 1.27832795e+01])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "within_cluster_sum_squares = np.zeros(test_n_cluster)\n",
    "\n",
    "for i in range(test_n_cluster):\n",
    "    n = i+1\n",
    "    kmeans = KMeans(n_clusters = n, init = 'k-means++', max_iter=500, n_init=20, random_state = 0)\n",
    "    kmeans.fit(X)\n",
    "    within_cluster_sum_squares[i] = kmeans.inertia_\n",
    "    \n",
    "within_cluster_sum_squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## making a plot to visualize the Elbow method\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcXFWd9/HPt/ds3UlI6G4TMAGjgqgRwuIGCg4GRHAdYZyBURwcB2R05Bm3cRkdn8EFdJxx8EGNgAqIoIKKQgYVHZUlYQ17WNMkJMGskK3T/Xv+uKdI0XR3LV3Vle76vl+vetWtU/fcc26nUr865557jiICMzOzSmiodQXMzGz8cFAxM7OKcVAxM7OKcVAxM7OKcVAxM7OKcVAxM7OKcVCxcUXSZyV9fxTKmSMpJDWl17+V9L5qlzsaKnkuki6Q9G+VOJaNDQ4qNqZIeirv0S9pa97rd1e4rAsk7RhQ5u2VLKNceUHtlgHpM1KdHynyOKMShK1+OKjYmBIRk3MP4DHgzXlpP6hCkV/KLzMiXl6FMkZikqQD8l7/FfBwrSpj5qBi41GLpIskbZZ0l6QFuTckPU/SFZLWSnpY0pkVLHdfSTdJ2ijpSknT88o9PtVlQ+pe2i+lv0fSz/L2Wy7psrzXKyTNH6bM7wGn5L0+Gbgof4ehzlnSQuATwLsGaYU9X9If0t/wWkkzCp1Leu8Vkm5J+X4ItBX3p7PxwkHFxqPjgUuBqcBVwH8BSGoAfgbcDswCjgI+JOmNFSr3ZOC9wPOAncDXU7kvBC4BPgTMBK4GfiapBbgeeK2kBkndQDPw6pRvH2AycMcwZX4fOFFSY/pynwLcmHtzuHOOiF8B/xf44SCtsL8C3gPsCbQAZxU6l3Q+PyULdNOBHwFvL+kvaGOeg4qNR/8bEVdHRB/ZF1zuy/JgYGZEfC4idkTEQ8C3gBOHOdZZ6Rd57nHhMPt+LyKWRcTTwKeAv5TUCLwL+EVELI6IXuArwATgVakOm4H5wBHANcDjkl6cXv8+IvqHKbMHuA94A1mL5aIB75dzzgDfjYj7I2IrcFmqH8OdC3AYWVD8WkT0RsTlwM0FyrFxpqnWFTCrgifytrcAbWmU1vOB50nakPd+I/D7YY71lYj4lyLLXZG3/SjZF+wMspbLo7k3IqJf0gqylgNkrZXXAS9I2xvIAsor0+tCLgL+luyL/XBgXt575ZwzPPdvODltD3cufcDj8exZah/F6oqDitWTFcDDETGv4J7l2Stve2+gF3gSWAm8NPeGJKV9H09J1wNvBuaSdUdtAN5NFlT+q4hyr0j7LY2IRyXln1+hcy51mvLhziWAWZKUF1j2Bh4ssQwbw9z9ZfXkJmCTpI9KmpCuQxwg6eAKHf+vJe0vaSLwOeDy1AV3GfAmSUdJagY+AmwH/pjyXQ+8HpgQET1krYiFwB7ArYUKTd1tRwKD3VtS6JxXA3PStZdiDHcufyK7lnSmpCZJbwMOKfK4Nk44qFjdSF/wbya7PvAwWSvi20DHMNn+ecB9Kk8Os+/3gAvIuo7agDNTufcBfw38ZyrzzWRDoXek9+8HniJ1SUXEJuAh4A+pzsWc25KIeE6LoIhz/lF6/vPAe16GKGfIc0nn8zayrrj1ZNdfflxM/W38kBfpMjOzSnFLxczMKsZBxczMKsZBxczMKsZBxczMKqbu7lOZMWNGzJkzp9bVMDMbU5YuXfpkRMwstF/dBZU5c+awZMmSWlfDzGxMkVTU7Aju/jIzs4pxUDEzs4pxUDEzs4pxUDEzs4pxUDEzs4pxUDEzs4pxUDEzs4pxUCnST299nO/f4EXszMyG46BSpF/cuYrv/clBxcxsOA4qReruaGPVxq21roaZ2W7NQaVIne1tbNq2ky07dta6KmZmuy0HlSJ1d7QB8MTGbTWuiZnZ7stBpUhd7SmobHJQMTMbioNKkbpSS2W1g4qZ2ZAcVIqUCyqr3P1lZjYkB5UiTWxpYkpbE6sdVMzMhuSgUoJsWLGDipnZUBxUStDZ3uZrKmZmw3BQKUF3R5tHf5mZDcNBpQRd7W2s3bydnX39ta6KmdluyUGlBJ0dbfQHrH1qe62rYma2W3JQKUG3hxWbmQ3LQaUEnemueg8rNjMbXNWCiqRFktZIWpaX9llJj0u6LT2OzXvv45KWS7pP0hvz0hemtOWSPpaXPlfSjZIekPRDSS3VOpec7o4JgFsqZmZDqWZL5QJg4SDpX42I+elxNYCk/YETgZekPP8tqVFSI/AN4Bhgf+CktC/AF9Ox5gHrgVOreC4ATJvYTEtTg4cVm5kNoWpBJSJ+B6wrcvcTgEsjYntEPAwsBw5Jj+UR8VBE7AAuBU6QJOBI4PKU/0LgLRU9gUFIorO91cOKzcyGUItrKmdIuiN1j01LabOAFXn79KS0odL3ADZExM4B6YOSdJqkJZKWrF27dkSV726f4O4vM7MhjHZQOQ/YF5gPrALOSekaZN8oI31QEXF+RCyIiAUzZ84srcYDdHb4rnozs6GMalCJiNUR0RcR/cC3yLq3IGtp7JW362xg5TDpTwJTJTUNSK+63PxfEUPGMDOzujWqQUVSd97LtwK5kWFXASdKapU0F5gH3ATcDMxLI71ayC7mXxXZN/pvgHek/KcAV47GOXS2t7FjZz8btvSORnFmZmNKU+FdyiPpEuB1wAxJPcBngNdJmk/WVfUI8H6AiLhL0mXA3cBO4PSI6EvHOQO4BmgEFkXEXamIjwKXSvo34FbgO9U6l3z5K0BOm1T1UcxmZmNK1YJKRJw0SPKQX/wR8QXgC4OkXw1cPUj6Q+zqPhs1XXlr1e/X3T7axZuZ7dZ8R32JngkqvlhvZvYcDiol2nNKK5LvqjczG4yDSomaGxuYMbnV83+ZmQ3CQaUMXe1erMvMbDAOKmXo6mjjCbdUzMyew0GlDG6pmJkNzkGlDF0dbWzc2svWHX21roqZ2W7FQaUM+TdAmpnZLg4qZci/AdLMzHZxUCnDrhsgt9a4JmZmuxcHlTI80/21cXuNa2JmtntxUCnDpNYmprQ18cRGt1TMzPI5qJTJw4rNzJ7LQaVMXR1tPLHJ3V9mZvkcVMrU1d7m7i8zswEcVMrU1dHG2s3b2dnXX+uqmJntNhxUytTV0UZ/wNqn3AVmZpbjoFKmXcOKfbHezCzHQaVMuRsgV3sEmJnZMxxUypRrqXgFSDOzXRxUyjR9UgstjQ2+V8XMLE/VgoqkRZLWSFqWl/ZlSfdKukPSTyRNTelzJG2VdFt6fDMvz0GS7pS0XNLXJSmlT5e0WNID6Xlatc5liPOjs6PV11TMzPJUs6VyAbBwQNpi4ICIeBlwP/DxvPcejIj56fH3eennAacB89Ijd8yPAddFxDzguvR6VGX3qjiomJnlVC2oRMTvgHUD0q6NiJ3p5Q3A7OGOIakbaI+IP0VEABcBb0lvnwBcmLYvzEsfNV0dE3yh3swsTy2vqbwX+GXe67mSbpV0vaTXprRZQE/ePj0pDaAzIlYBpOc9hypI0mmSlkhasnbt2oqdQFd7K6s2biOLd2ZmVpOgIumTwE7gBylpFbB3RLwC+CfgYkntgAbJXvI3eEScHxELImLBzJkzy632c3S2t7F9Zz8bt/ZW7JhmZmPZqAcVSacAxwHvTl1aRMT2iPhz2l4KPAi8kKxlkt9FNhtYmbZXp+6xXDfZmtE5g126OyYAHlZsZpYzqkFF0kLgo8DxEbElL32mpMa0vQ/ZBfmHUrfWZkmHpVFfJwNXpmxXAaek7VPy0kdNV0cr4LXqzcxymqp1YEmXAK8DZkjqAT5DNtqrFVicRgbfkEZ6HQ58TtJOoA/4+4jIXeT/ANlIsglk12By12HOBi6TdCrwGPDOap3LULpSS2W1WypmZkAVg0pEnDRI8neG2PcK4Ioh3lsCHDBI+p+Bo0ZSx5Hac0orkru/zMxyfEf9CDQ3NrDHpFYPKzYzSxxURqi7o80tFTOzxEFlhDrb29xSMTNLHFRGqLujzaO/zMwSB5UR6upoY8OWXrb19tW6KmZmNeegMkKdXgHSzOwZDioj1N3hxbrMzHIcVEYo11LxxXozMweVEetyS8XM7BkOKiM0ubWJKa1NbqmYmeGgUhFdHV4B0swMHFQqoqujjVVuqZiZOahUQmd7m2cqNjPDQaUiujvaWLN5Gzv7+mtdFTOzmnJQqYDO9jb6A558aketq2JmVlMOKhWQuwHSc4CZWb1zUKmAXVO1bK1xTczMastBpQJyN0B6WLGZ1TsHlQqYPrGFlsYGDys2s7rnoFIBDQ1iz/ZWDys2s7rnoFIhXqzLzKzKQUXSIklrJC3LS5suabGkB9LztJQuSV+XtFzSHZIOzMtzStr/AUmn5KUfJOnOlOfrklTN8xlOZ7unajEzq3ZL5QJg4YC0jwHXRcQ84Lr0GuAYYF56nAacB1kQAj4DHAocAnwmF4jSPqfl5RtY1qjpas9aKhFRqyqYmdVcVYNKRPwOWDcg+QTgwrR9IfCWvPSLInMDMFVSN/BGYHFErIuI9cBiYGF6rz0i/hTZN/lFeccadV0dbWzr7Wfj1t5aVcHMrOZqcU2lMyJWAaTnPVP6LGBF3n49KW249J5B0p9D0mmSlkhasnbt2oqcxEBdvgHSzGy3ulA/2PWQKCP9uYkR50fEgohYMHPmzBFUcWjdvlfFzKwmQWV16roiPa9J6T3AXnn7zQZWFkifPUh6Tey6q95BxczqVy2CylVAbgTXKcCVeeknp1FghwEbU/fYNcDRkqalC/RHA9ek9zZLOiyN+jo571ijbs8p7v4yM2uq5sElXQK8DpghqYdsFNfZwGWSTgUeA96Zdr8aOBZYDmwB3gMQEeskfR64Oe33uYjIXfz/ANkIswnAL9OjJlqaGpgxudUtFTOra1UNKhFx0hBvHTXIvgGcPsRxFgGLBklfAhwwkjpWUldHq1sqZlbXdqcL9WNeV/sEt1TMrK45qFSQWypmVu8cVCqoq72NDVt62dbbV+uqmJnVhINKBXV1TAA8rNjM6tewQUXSwZK68l6fLOnKNHnj9OpXb2zpavewYjOrb4VaKv8P2AEg6XCy4cAXARuB86tbtbEnN1XLagcVM6tThYYUN+bdE/Iu4PyIuAK4QtJt1a3a2JMLKqvc/WVmdapQS6VRUi7wHAX8Ou+9qt7jMhZNbm1iSmuTr6mYWd0qFBguAa6X9CSwFfg9gKQXkHWB2QCdHV6sy8zq17BBJSK+IOk6oBu4NnatQNUAfLDalRuLcot1mZnVo2GDiqSJwNKI6E2vX0Q2P9ejEfHjUajfmNPV0cYflj9Z62qYmdVEoWsqvwLmwDNdXn8C9gFOl/Tv1a3a2NTV3saazdvp6/eywmZWfwoFlWkR8UDaPgW4JCI+SLae/HFVrdkY1dXRRl9/8ORT22tdFTOzUVcoqOT/3D6SbH14ImIH0F+tSo1luRsgPazYzOpRodFfd0j6CvA48ALgWgBJU6tdsbGqK39Z4b0K7GxmNs4Uaqn8HfAk2XWVoyNiS0rfH/hKFes1ZvmuejOrZ4VaKpOBn0XEXQPSN5FdxLcBpk9soblR7v4ys7pUqKXyn8CMQdJnAf9R+eqMfQ0NorO9zS0VM6tLhYLKSyPi+oGJEXEN8LLqVGns62pvY9XGrbWuhpnZqCsUVJrLfK+udXa0sXqThxSbWf0pFFQekHTswERJxwAPVadKY193ezb/165ZbczM6kOhC/UfAn4h6S+BpSltAfBKyrz5MU318sO8pH2ATwNTyUabrU3pn4iIq1OejwOnAn3Aman7DUkLya7tNALfjoizy6lTpXV1tLG1t49NW3fSMdENOjOrH4WCypvIvsxfDLwopV0PvD8iyroSHRH3AfMBJDWS3QPzE+A9wFcj4llDlSXtD5wIvAR4HvA/kl6Y3v4G8BdAD3CzpKsi4u5y6lVJz9yrsmmbg4qZ1ZVC3V+zgS8CXyJroewAVgMTK1T+UcCDEfHoMPucAFwaEdsj4mFgOXBIeiyPiIfSHf6Xpn1rbtdd9b5Yb2b1ZdigEhFnRcSrgE7gE8A64L3AMkmVaBGcSLZmS84Zku6QtEjStJQ2C1iRt09PShsq/TkknSZpiaQla9euHWyXiups9w2QZlafCrVUciYA7UBHeqwEbhxJwZJagOOBH6Wk84B9ybrGVgHn5HYdJHsMk/7cxIjzI2JBRCyYOXPmSKpdlE7P/2VmdarQeirnk13L2EwWRP4InBsR6ytQ9jHALRGxGiD3nMr9FvDz9LKHZ8+iNZssqDFMek21NDUwY3KLWypmVncKtVT2BlqBJ8guqPcAGypU9knkdX1J6s57763AsrR9FXCipFZJc4F5wE3AzcA8SXNTq+fEtO9uocvLCptZHSq0nPBCSSJrrbwK+AhwgKR1wJ8i4jPlFJpWlPwL4P15yV+SNJ+sC+uR3HsRcZeky4C7gZ3A6RHRl45zBnAN2ZDiRYPMUVYzXe1t9Kz3hXozqy+FhhST1qVfJmkDsDE9jiMbfVVWUEmzHe8xIO1vhtn/C8AXBkm/Gri6nDpUW2d7G0sfrUQvoZnZ2FHomsqZZC2UVwO9wB/IlhReBNxZ9dqNYd0dbazf0su23j7amhtrXR0zs1FRqKUyB7gc+HBErKp+dcaP/GHFz99jUo1rY2Y2OgpdU/mn0arIeNPdMQHIVoB0UDGzelHsfSpWoq6OViCbqsXMrF44qFRJrvvLw4rNrJ44qFTJlLZmJrc2+a56M6srDipV1Nne6rvqzayuOKhUUXfHBF9TMbO64qBSRZ3tnqrFzOqLg0oVdXe0sWbzdvr6vaywmdUHB5Uq6uxoo68/ePKp7bWuipnZqHBQqaIuDys2szrjoFJF3Xlr1ZuZ1QMHlSryDZBmVm8cVKpoj0ktNDfKLRUzqxsOKlXU0CD2nOJhxWZWPxxUqszLCptZPXFQqbKujjZP1WJmdcNBpcq62ttYtXEb2arMZmbjm4NKlXV3tLG1t49N23bWuipmZlXnoFJlHlZsZvWkZkFF0iOS7pR0m6QlKW26pMWSHkjP01K6JH1d0nJJd0g6MO84p6T9H5B0Sq3OZyhdvgHSzOpIrVsqr4+I+RGxIL3+GHBdRMwDrkuvAY4B5qXHacB5kAUh4DPAocAhwGdygWh3kZuqZbVbKmZWB2odVAY6AbgwbV8IvCUv/aLI3ABMldQNvBFYHBHrImI9sBhYONqVHk6u+8srQJpZPahlUAngWklLJZ2W0jojYhVAet4zpc8CVuTl7UlpQ6U/i6TTJC2RtGTt2rUVPo3htTQ1MGNyi7u/zKwuNNWw7FdHxEpJewKLJd07zL4aJC2GSX92QsT5wPkACxYsGPWxvdliXVtHu1gzs1FXs5ZKRKxMz2uAn5BdE1mdurVIz2vS7j3AXnnZZwMrh0nfrXS1t/HEJq+pYmbjX02CiqRJkqbktoGjgWXAVUBuBNcpwJVp+yrg5DQK7DBgY+oeuwY4WtK0dIH+6JS2W/Fd9WZWL2rV/dUJ/ERSrg4XR8SvJN0MXCbpVOAx4J1p/6uBY4HlwBbgPQARsU7S54Gb036fi4h1o3caxelqb2Pd0zvY1ttHW3NjratjZlY1NQkqEfEQ8PJB0v8MHDVIegCnD3GsRcCiStexknL3qqzZtJ2995hY49qYmVXP7jakeFzKBZVVvlhvZuOcg8ooeGatel9XMbNxzkFlFORaKr5Yb2bjnYPKKJjS1syklkbfVW9m456DyijxsGIzqwcOKqOkq6ONh9Y+7cW6zGxcc1AZJUfv38W9T2zmV8ueqHVVzMyqxkFllLz70L3Zv7udz/38bp7e7lUgzWx8clAZJU2NDXz+LQewauM2vv7rB2pdHTOzqnBQGUUHPX8af7lgNt/5/cM8sHpzratjZlZxDiqj7KMLX8yk1iY+deUyX7Q3s3HHQWWU7TG5lX9e+CJueGgdV92+283Sb2Y2Ig4qNXDiwXvz8tkd/Nsv7mHTtt5aV8fMrGIcVGqgsUF8/i0H8ORT2/nq4vtrXR0zs4pxUKmRl82eyrsP3ZsL//gId6/cVOvqmJlVhINKDf2fo1/MtIktfOrKZfT3+6K9mY19Dio11DGxmY8d82KWPrqey2/pqXV1zMxGzEGlxt5+4GwWPH8aZ//yXjZs2VHr6piZjYiDSo01pIv2G7f28uVr7qt1dczMRsRBZTewX3c7p7xyDhff9Bi3r9hQ6+qYmZXNQWU38eG/mMfMya186spl9PmivZmNUaMeVCTtJek3ku6RdJekf0zpn5X0uKTb0uPYvDwfl7Rc0n2S3piXvjClLZf0sdE+l0qa0tbMJ9+0H3f0bOSSmx6rdXXMzMpSi5bKTuAjEbEfcBhwuqT903tfjYj56XE1QHrvROAlwELgvyU1SmoEvgEcA+wPnJR3nDHp+Jc/j1fuswdfvuY+/vzU9lpXx8ysZKMeVCJiVUTckrY3A/cAs4bJcgJwaURsj4iHgeXAIemxPCIeiogdwKVp3zFLEp874SU8vX0nZ//y3lpXx8ysZDW9piJpDvAK4MaUdIakOyQtkjQtpc0CVuRl60lpQ6UPVs5pkpZIWrJ27doKnkHlzeucwqmvncuPlvaw5JF1ta6OmVlJahZUJE0GrgA+FBGbgPOAfYH5wCrgnNyug2SPYdKfmxhxfkQsiIgFM2fOHHHdq+3MI+fR3dHGv/x0GTv7+mtdHTOzotUkqEhqJgsoP4iIHwNExOqI6IuIfuBbZN1bkLVA9srLPhtYOUz6mDeptYlPH7c/9z6xmYv+9Gitq2NmVrRajP4S8B3gnog4Ny+9O2+3twLL0vZVwImSWiXNBeYBNwE3A/MkzZXUQnYx/6rROIfRsPCALg5/4UzOXXw/azZtq3V1zMyKUouWyquBvwGOHDB8+EuS7pR0B/B64MMAEXEXcBlwN/Ar4PTUotkJnAFcQ3ax/7K077ggiX89/iXs2NnPF66+p9bVMTMriuptSdsFCxbEkiVLal2Nop177X18/dfLufjvDuVV+86odXXMrE5JWhoRCwrt5zvqd3P/8PoXsNf0CXz6yrvYsdMX7c1s9+agsptra27ks29+CcvXPMU3r3+QemtZmtnY4qAyBhy1XyfHHNDFuYvv570X3MwjTz5d6yqZmQ3KQWWM+PpJr+Bf3rQfNz+ynqO/+jvOufY+tu7oq3W1zMyexUFljGhubOB9r92HX3/kCI59aRf/+evlvOHc67nmrifcJWZmuw0HlTFmz/Y2vnbiK7j0tMOY3NrE+7+3lL/97s08tPapWlfNzMxBZaw6bJ89+PmZr+FTx+3PLY+uZ+HXfs+XfnUvW3bsrHXVzKyOOaiMYc2NDZz6mrlcd9YRHPeybv77tw/yhnOu55d3rnKXmJnVhIPKOLDnlDbOfdd8fvT3r6R9QjMf+MEtnLzoJh50l5iZjTIHlXHk4DnT+fkHX8Nn37w/t63YwMKv/Y6zf3kvT293l5iZjQ4HlXGmqbGBv331XH79kddxwvxZfPP6BznqnOv5+R0r3SVmZlXnoDJOzZzSylfe+XKu+MArmT6phTMuvpW3n/dHLluygqfccjGzKvGEknWgrz+4+KbHWPS/D/Pwk08zobmRY17axTsOms1hc/egoWGw9c7MzHYpdkJJB5U6EhHc8th6Ll/aw89vX8Xm7TuZNXUCbz9oNu84cDZ77zGx1lU0s92Ug8oQ6jmo5Nu6o49r736Cy5f28L/LnyQCDp07nXccNJtjX9rNpNamWlfRzHYjDipDcFB5rpUbtvKTWx/n8qU9PPzk00xsaeSYA7p5x0GzOXTudHePmZmDylAcVIaW3z32s9tX8dT2new1fQJvP3A2bz9wNntNd/eYWb1yUBmCg0pxBusee2HnZGZMbmXaxBamTmx+5nn6pJZnpU2b1EJ7WxOSWzhm40WxQcUd5zaoCS2NnDB/FifMn8XKDVv58S093LZiIxu27ODeJzaxYUsv67fsoH+I3ySNDWLqhOZngs7UiS1MeyYQZdtTJzan7dzrFlqaPMrdbCxzULGCnjd1AmccOe856f39weZtO1m/Zceux9O9ea972ZDSVqzbwh09WdpwyyJPamnMAs2kgQGohakTmp9Jzz2mTmpmSqtbRWa7CwcVK1tDg+iY2EzHxGbmMKmoPBHB1t4+1m/pZf3TO9i4tXdXAHo6LxCltBXrtrB+Sy+btvUyVE9tU4PyWj3NQ7SKUtqkrJtu6gS3isyqYcwHFUkLgf8AGoFvR8TZNa6SDUMSE1uamNjSxKypE4rO19cfzwSgDXktolw3XH4wKrZVNLm16VnXhgYPSs1uFZmVYEwHFUmNwDeAvwB6gJslXRURd9e2ZlZpjQ1i+qQWpk9qKTrPwFZRLgBtSEFoYFDKtYo2bu0dth65a0W54NMx4bnXiKZOfPY+E5obHYysLozpoAIcAiyPiIcAJF0KnAA4qFhVWkUbtmYBaGMKRis3bOPulZvYsLWXLTv6hjxmS1MDk1ubKBRWCsed5+4wWJ6hDlPo+CpQw8L5C7w/wsA66LkOev6DlzPSuD7S8yvjn7ei5S865eCqz5wx1oPKLGBF3use4NCBO0k6DTgNYO+99x6dmtmYVU6rCGBbbx8bt/bmtYiyLrkNKUAVWoKg0Oj+wd4ePM/gByp4/ILlD79DOfUfcfnFJaXjj+z2iZHXv1D+An/fAvkL78CoXEcc60FlsLD8nD9tRJwPnA/ZfSrVrpTVp7bmRtqaG+lsb6t1VcxqZqwPf+kB9sp7PRtYWaO6mJnVvbEeVG4G5kmaK6kFOBG4qsZ1MjOrW2O6+ysidko6A7iGbEjxooi4q8bVMjOrW2M6qABExNXA1bWuh5mZjf3uLzMz2404qJiZWcU4qJiZWcU4qJiZWcXU3SJdktYCj5aZfQbw5AiKd37nd37nH6v5nx8RMwvuFRF+FPkAlji/8zu/89dj/mIf7v4yM7OKcVAxM7OKcVApzfnO7/zO7/x1mr8odXeh3szMqsctFTMzqxgHFTMzqxgHlSJIWiRpjaRlZebfS9JvJN0j6S5J/1hi/jZJN0m6PeX/1zLq0CjpVkk/LzVvyv+IpDsl3SZpSRn5p0q6XNK96e/wyhLyviiVm3tmG9bdAAAJ0ElEQVRskvShEsv/cPrbLZN0iaSSVtKS9I8p713FlD3YZ0bSdEmLJT2QnqeVmP+dqfx+SQvKKP/L6e9/h6SfSJpaYv7Pp7y3SbpW0vNKyZ/33lmSQtKMEsv/rKTH8z4Hx5ZavqQPSrov/R2/VGL5P8wr+xFJt5WYf76kG3L/hyQdUmL+l0v6U/p/+DNJ7cPkH/Q7p5TPYNlGY9zyWH8AhwMHAsvKzN8NHJi2pwD3A/uXkF/A5LTdDNwIHFZiHf4JuBj4eZnn8AgwYwR/wwuB96XtFmBqmcdpBJ4guxGr2DyzgIeBCen1ZcDflpD/AGAZMJFsZu//AeaV+pkBvgR8LG1/DPhiifn3A14E/BZYUEb5RwNNafuLZZTfnrd9JvDNUvKn9L3Ilqp4dLjP0xDlfxY4q8h/s8Hyvz7927Wm13uWWv+8988BPl1i+dcCx6TtY4Hflpj/ZuCItP1e4PPD5B/0O6eUz2C5D7dUihARvwPWjSD/qoi4JW1vBu4h+6IrNn9ExFPpZXN6FD3CQtJs4E3At4uudAWlX1SHA98BiIgdEbGhzMMdBTwYEaXOitAETJDURBYcSlkhdD/ghojYEhE7geuBtw6XYYjPzAlkwZX0/JZS8kfEPRFxXzEVHiL/tan+ADeQrZRaSv5NeS8nMcxncJj/M18F/nm4vAXyF2WI/B8Azo6I7WmfNeWUL0nAXwKXlJg/gFzrooNhPoND5H8R8Lu0vRh4+zD5h/rOKfozWC4HlVEmaQ7wCrLWRin5GlNzew2wOCJKyf81sv/I/aWUOUAA10paKum0EvPuA6wFvpu64L4taVKZ9TiRYf4zDyYiHge+AjwGrAI2RsS1JRxiGXC4pD0kTST7lblXgTyD6YyIValOq4A9yzhGpbwX+GWpmSR9QdIK4N3Ap0vMezzweETcXmq5ec5IXXCLyui6eSHwWkk3Srpe0sFl1uG1wOqIeKDEfB8Cvpz+fl8BPl5i/mXA8Wn7nRT5GRzwnVP1z6CDyiiSNBm4AvjQgF99BUVEX0TMJ/t1eYikA4os8zhgTUQsLbnCz/bqiDgQOAY4XdLhJeRtImvKnxcRrwCeJmt6l0TZktHHAz8qMd80sl9oc4HnAZMk/XWx+SPiHrLuosXAr4DbgZ3DZtqNSfokWf1/UGreiPhkROyV8p5RQpkTgU9SYiAa4DxgX2A+2Y+Dc0rM3wRMAw4D/g9wWWp1lOokSvxhk3wA+HD6+32Y1HIvwXvJ/u8tJevS2lEow0i+c8rloDJKJDWT/eP+ICJ+XO5xUrfRb4GFRWZ5NXC8pEeAS4EjJX2/jHJXpuc1wE+AIS8yDqIH6MlrXV1OFmRKdQxwS0SsLjHfG4CHI2JtRPQCPwZeVcoBIuI7EXFgRBxO1i1R6q9UgNWSugHS85DdL9Ui6RTgOODdkTrWy3Qxw3S/DGJfsqB+e/oszgZukdRV7AEiYnX6cdUPfIvSPoOQfQ5/nLqTbyJruQ85WGAwqfv0bcAPSywb4BSyzx5kP4xKqn9E3BsRR0fEQWRB7cECdR3sO6fqn0EHlVGQfg19B7gnIs4tI//M3EgdSRPIviTvLSZvRHw8ImZHxByyrqNfR0TRv9JTmZMkTcltk13wLXokXEQ8AayQ9KKUdBRwdyl1SMr9hfgYcJikienf4iiyPuaiSdozPe9N9qVSTj2uIvtiIT1fWcYxyiZpIfBR4PiI2FJG/nl5L4+nyM8gQETcGRF7RsSc9FnsIbuQ/EQJ5XfnvXwrJXwGk58CR6ZjvZBswEips/a+Abg3InpKzAfZNZQj0vaRlPjDJO8z2AD8C/DNYfYd6jun+p/BSl/5H48Psi+QVUAv2X+GU0vM/xqyaxJ3ALelx7El5H8ZcGvKv4xhRp0UOM7rKGP0F9k1kdvT4y7gk2UcYz6wJJ3DT4FpJeafCPwZ6Cjz3P+V7EtwGfA90gigEvL/niwQ3g4cVc5nBtgDuI7sy+Q6YHqJ+d+atrcDq4FrSsy/HFiR9xkcbvTWYPmvSH+/O4CfAbPK/T9DgdGEQ5T/PeDOVP5VQHeJ+VuA76dzuAU4stT6AxcAf1/mv/9rgKXpM3QjcFCJ+f+RbBTX/cDZpBlRhsg/6HdOKZ/Bch+epsXMzCrG3V9mZlYxDipmZlYxDipmZlYxDipmZlYxDipmZlYxDio27qQZcM/Je32WpM9W6NgXSHpHJY5VoJx3phlmf1PNekmaI+mvSq+h2eAcVGw82g68TcNMrV4LkhpL2P1U4B8i4vXVqk8yBygpqJR4HlZnHFRsPNpJth73hwe+MfAXvaSn0vPr0iSDl0m6X9LZkt6tbB2bOyXtm3eYN0j6fdrvuJS/Udl6JTenCQ/fn3fc30i6mOzGvYH1OSkdf5mkL6a0T5PdvPZNSV8eJM8/pzy3Szp7kPcfyQVUSQsk/TZtH6Fd64HcmmZJOJtsksXblK05U9R5pFkWfpHqsEzSu4r5h7Hxr6nWFTCrkm8Ad2iYhZgG8XKyae7XAQ8B346IQ5QtcPRBsllmIft1fwTZfFa/kfQC4GSy2Y8PltQK/EFSbibkQ4ADIuLh/MKULXL1ReAgYD3ZLNBviYjPSTqSbO2QJQPyHEM2XfmhEbFF0vQSzu8s4PSI+EOaaHAb2cSeZ0VELjieVsx5SHo7sDIi3pTydZRQDxvH3FKxcSmyGVkvIltMqlg3R7YOxXayyfpyX6Z3kgWSnMsioj+yqc8fAl5MNh/aycqWJ7iRbDqM3FxZNw0MKMnBZAs1rY1snZMfkK07M5w3AN+NNHdXRJSy5sgfgHMlnUm2SNpgMy0Xex53krXYvijptRGxsYR62DjmoGLj2dfIrk3kr92yk/S5T5PuteS9tz1vuz/vdT/PbtUPnNsoyFbn/GBEzE+PubFrzZanh6hfOdOua5DyB3rmHIFnlk2OiLOB9wETgBskvXiI4xc8j4i4n6yFdSfw76nLzsxBxcav9Cv+MrLAkvMI2ZchZGusNJdx6HdKakjXWfYB7iNbIvcDabpxJL1QhRciuxE4QtKMdPH7JLJVJYdzLfBeZeuTMET31yPsOsdnpqeXtG9kswV/kWxyzxcDm8nW5sgp6jxS192WiPg+2YJT5SxlYOOQr6nYeHcOz15M6lvAlZJuIpuldahWxHDuI/vy7ySbsXabpG+TdZHdklpAaymwVGtErJL0ceA3ZC2EqyNi2KnII+JXkuYDSyTtAK4GPjFgt38FviPpEzx7hdEPSXo90Ec24/IvyVphOyXdTjYD738UeR4vJVvFsJ9sJt0PDFdvqx+epdjMzCrG3V9mZlYxDipmZlYxDipmZlYxDipmZlYxDipmZlYxDipmZlYxDipmZlYx/x9vmjCo0Wn/cgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(within_cluster_sum_squares)\n",
    "\n",
    "xticks = np.arange(test_n_cluster)\n",
    "xticklabels = xticks + 1\n",
    "ax.set(xticks=xticks, xticklabels=xticklabels, title='The Elbow Method', \n",
    "       xlabel='Number of clusters', ylabel='WCSS');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.        , -0.95775705, -0.90309048, -0.97522825, -0.91109173,\n",
       "       -0.99999775, -1.        , -0.99999978, -0.99999708, -0.99331233,\n",
       "       -0.99999999, -1.        , -1.        , -0.99999999, -1.        ,\n",
       "       -0.99999994, -1.        , -1.        , -0.99999995, -1.        ])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Numbers ase\n",
    "\n",
    "# wcssScaled= [x*clusterN/max(wcss) for x in wcss ]\n",
    "within_cluster_sum_squares_scaled = test_n_cluster*within_cluster_sum_squares/within_cluster_sum_squares.max()\n",
    "wcssScaled = within_cluster_sum_squares_scaled\n",
    "\n",
    "cosines = -1 * np.ones(test_n_cluster)\n",
    "\n",
    "\n",
    "for i in range(test_n_cluster-1):\n",
    "   # check if the point is below a segment midpoint connecting its neighbors\n",
    "   if (wcssScaled[i] < (wcssScaled[i+1]+wcssScaled[i-1])/2 ):\n",
    "       cosines[i]= (-1+(wcssScaled[i-1]-wcssScaled[i])*(wcssScaled[i+1]-wcssScaled[i]))/ \\\n",
    "       ((1+(wcssScaled[i-1]-wcssScaled[i])**2)*(1+ (wcssScaled[i+1]-wcssScaled[i])**2))**.5\n",
    "cosines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.00000000e+01, 4.87392168e+00, 2.19864821e+00, 1.23053139e+00,\n",
       "       6.22705469e-01, 5.00824730e-01, 3.81097300e-01, 2.60704181e-01,\n",
       "       1.40976752e-01, 2.36994939e-02, 2.26718804e-02, 2.17737098e-02,\n",
       "       2.08965079e-02, 1.99685304e-02, 1.91699482e-02, 1.82540963e-02,\n",
       "       1.76777809e-02, 1.68462770e-02, 1.60060076e-02, 1.54740821e-02])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "within_cluster_sum_squares_scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.00000000e+01, 4.87392168e+00, 2.19864821e+00, 1.23053139e+00,\n",
       "       6.22705469e-01, 5.00824730e-01, 3.81097300e-01, 2.60704181e-01,\n",
       "       1.40976752e-01, 2.36994939e-02, 2.26718804e-02, 2.17737098e-02,\n",
       "       2.08965079e-02, 1.99685304e-02, 1.91699482e-02, 1.82540963e-02,\n",
       "       1.76777809e-02, 1.68462770e-02, 1.60060076e-02])"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "within_cluster_sum_squares_scaled[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "within_cluster_sum_squares_scaled[1:] < within_cluster_sum_squares_scaled[:-1] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}