Customer Segmentation for Market Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Customer Segments\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project we will analyze a dataset containing annual spending amounts for internal structure, to understand the variation in the different types of customers that a wholesale distributor interacts with.\n",
    "\n",
    "The dataset can be downloaded from : https://archive.ics.uci.edu/ml/datasets/Wholesale+customers\n",
    "\n",
    "It contains the folliwing attributes:\n",
    "1.\tFRESH: annual spending (m.u.) on fresh products (Continuous)\n",
    "2.\tMILK: annual spending (m.u.) on milk products (Continuous)\n",
    "3.\tGROCERY: annual spending (m.u.)on grocery products (Continuous)\n",
    "4.\tFROZEN: annual spending (m.u.)on frozen products (Continuous) \n",
    "5.\tDETERGENTS_PAPER: annual spending (m.u.) on detergents and paper products (Continuous) \n",
    "6.\tDELICATESSEN: annual spending (m.u.)on and delicatessen products (Continuous)\n",
    "7.\tCHANNEL: customer™ Channel - Horeca (Hotel/Restaurant/Cafe) or Retail channel (Nominal) \n",
    "8.\tREGION: customers™ Region -  Lisnon, Oporto or Other (Nominal)\n",
    "\n",
    "We would not be using the 2 columns 'Channel' and 'Region' as they represent classes. Instead we would use the other 6 attributes for customer clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import libraries: NumPy, pandas, matplotlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "# Tell iPython to include plots inline in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# read .csv from provided dataset\n",
    "csv_filename=\"Wholesale customers data.csv\"\n",
    "\n",
    "# df=pd.read_csv(csv_filename,index_col=0)\n",
    "df=pd.read_csv(csv_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Channel</th>\n",
       "      <th>Region</th>\n",
       "      <th>Fresh</th>\n",
       "      <th>Milk</th>\n",
       "      <th>Grocery</th>\n",
       "      <th>Frozen</th>\n",
       "      <th>Detergents_Paper</th>\n",
       "      <th>Delicassen</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>12669</td>\n",
       "      <td>9656</td>\n",
       "      <td>7561</td>\n",
       "      <td>214</td>\n",
       "      <td>2674</td>\n",
       "      <td>1338</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>7057</td>\n",
       "      <td>9810</td>\n",
       "      <td>9568</td>\n",
       "      <td>1762</td>\n",
       "      <td>3293</td>\n",
       "      <td>1776</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>6353</td>\n",
       "      <td>8808</td>\n",
       "      <td>7684</td>\n",
       "      <td>2405</td>\n",
       "      <td>3516</td>\n",
       "      <td>7844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>13265</td>\n",
       "      <td>1196</td>\n",
       "      <td>4221</td>\n",
       "      <td>6404</td>\n",
       "      <td>507</td>\n",
       "      <td>1788</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>22615</td>\n",
       "      <td>5410</td>\n",
       "      <td>7198</td>\n",
       "      <td>3915</td>\n",
       "      <td>1777</td>\n",
       "      <td>5185</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Channel  Region  Fresh  Milk  Grocery  Frozen  Detergents_Paper  Delicassen\n",
       "0        2       3  12669  9656     7561     214              2674        1338\n",
       "1        2       3   7057  9810     9568    1762              3293        1776\n",
       "2        2       3   6353  8808     7684    2405              3516        7844\n",
       "3        1       3  13265  1196     4221    6404               507        1788\n",
       "4        2       3  22615  5410     7198    3915              1777        5185"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Fresh', 'Milk', 'Grocery', 'Frozen', 'Detergents_Paper', 'Delicassen'], dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features = df.columns[2:]\n",
    "features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Fresh  Milk  Grocery  Frozen  Detergents_Paper  Delicassen\n",
      "0  12669  9656     7561     214              2674        1338\n",
      "1   7057  9810     9568    1762              3293        1776\n",
      "2   6353  8808     7684    2405              3516        7844\n",
      "3  13265  1196     4221    6404               507        1788\n",
      "4  22615  5410     7198    3915              1777        5185\n"
     ]
    }
   ],
   "source": [
    "data = df[features]\n",
    "print(data.head(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first PCA dimension is the dimension in the data with highest variance. Intuitively, it corresponds to the 'longest' vector one can find in the 6-dimensional feature space that captures the data, that is, the eigenvector with the largest eigenvalue.\n",
    "\n",
    "The first component will carry a high load of the 'Fresh' feature, as this feature seems to vary more than any of the other features (according to the README.md-file, 'Fresh' has the highest variance). Moreover, this feature seems to vary independently of the others, that is, a high or low value of 'Fresh' is not very informative for the values of the other features. A pratical interpretation of these observations could be that some of the supplied customers focus on fresh items, whereas other focus on non-fresh items.\n",
    "\n",
    "ICA, as opposed to PCA, finds the subcomponents that are statistically independent. ICA also finds 'Fresh' as one of the first components. The other components, however, may differ, as they need not be orthogonal in the feature space (in contrast to PCA). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.97653685  0.12118407  0.06154039  0.15236462 -0.00705417  0.06810471]\n",
      " [-0.11061386  0.51580216  0.76460638 -0.01872345  0.36535076  0.05707921]\n",
      " [-0.17855726  0.50988675 -0.27578088  0.71420037 -0.20440987  0.28321747]\n",
      " [ 0.04187648  0.64564047 -0.37546049 -0.64629232 -0.14938013  0.02039579]\n",
      " [-0.015986   -0.20323566  0.1602915  -0.22018612 -0.20793016  0.91707659]\n",
      " [-0.01576316  0.03349187  0.41093894 -0.01328898 -0.87128428 -0.26541687]]\n",
      "[ 0.45961362  0.40517227  0.07003008  0.04402344  0.01502212  0.00613848]\n"
     ]
    }
   ],
   "source": [
    "# Apply PCA with the same number of dimensions as variables in the dataset\n",
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=6) # 6 components for 6 variables\n",
    "pca.fit(data)\n",
    "\n",
    "# Print the components and the amount of variance in the data contained in each dimension\n",
    "print(pca.components_)\n",
    "print(pca.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The explained variance is high for the first two dimensions (45.96 % and 40.52 %, respectively), but drops significantly beginning with the third dimension (7.00 % for the third, 4.40 % for the fourth dimension). Thus, the first two components explain already 86.5 % of the variation in the data.\n",
    "\n",
    "How many dimension to choose for the analysis really depends on the goal of the analysis. Even though PCA reduces the feature space (with all advantages that brings, such as faster computations) and makes interpreting the data easier for us by projecting them down to a lower dimension, it necessarily comes with a loss of information that may or may not be desired.\n",
    "\n",
    "It the case at hand, assuming interpretation is the goal (creating customer segments) and given the sharp drop of the explained variance after the second component, we would choose the first two dimensions for analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGHCAYAAAAdnkAlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XecFPX9x/HX5xBEitgFRWk2LKhgQ43GEmvQRIyKEjuK\nqOjdGZMolqgk+amAqCiWKBKVaMQCNuwlihrBRsQGqLFhA4IUKff5/fGd9fb29u52h92bK+/n4zGP\n2536mdnZ2c995/v9jrk7IiIiIoVWknQAIiIi0jQpyRAREZGiUJIhIiIiRaEkQ0RERIpCSYaIiIgU\nhZIMERERKQolGSIiIlIUSjJERESkKJRkiIiISFEoyagnZjbOzObEXPZjM7ut0DHluO3YcRdLksej\nsTCzS82sIuk48mFmO5nZS2b2g5mtNLNeSceUTWM4/8xsMzN7wszmR8fysKRjkuapWSUZZnaCmVXU\nMKw0s12KuHmPhjgqVmHZVbUqcRdLksejwTCzNczsEjPbK8tkJxynRsHMVgPuA9YGzgV+C3ySYDx9\no2O7ZpbJjeH8Gw9sA1xAOJavZ5vJzLpkXAdXmNknZna/mW2fZf7VzazUzF6JEpglZva+mV1nZpvX\nsI0ro3VPKOQOSnxmdrCZXVIf21qtPjbSwDhwEfBxlmkf1W8oOduSRvSDUQ90PII2wCWEc/qFjGmX\nA3+p94ji6wFsCpzi7rcnHQywO3AxcDvwv4xpDfr8M7PWwG7A5e5+Q46L3Q08CrQAegJDgIPMbDd3\nfzta77rAFGBH4GHgLuAHwvE4BhgEtM6y7mOAOUA/M2vr7ovi7psUzCGEz/hPxd5Qc0wyAB539+lJ\nB5Erd1+edAwNgZm1dvelTfV4mFkLoCSP/bOaJrh7BbCsIIHVjw2jvwsSjaJSbce2oZ9/G0R/8zmW\n09397tQbM3sZmAScEQ0AdwDbA/3d/cH0hc3sImB45krNbB9gY2Bf4AngCODvecQlxVHj+V1w7t5s\nBuAEYCXQu475Lo3m2ydj/M3Aj8B20fu9Cf/RHAX8GfiSkNk/BHTOWPZ2YHbGuPOAl4BvgcWEIs3+\nWeL5GLgtYz8qCP9tjQS+jrZ7P7BuluUPJvyn+wPhv7KHga2zzPcrYAawBHg7el8t7izLTQZm1TBt\nKvBa2vuTgKeBucBS4D/A4Br2eRJwAPDvKKahNRyPtYGro5gXEi6ujwK9MtaZ+rx+A1wI/Dda71NA\njywx7Bqt5/vo2L2ViiFtni0JxfzfRev6N9Avh3OxSxRLGXAOoRRtOdALaAlcFp0P86NtvwD8PMvy\nK6O/qeHitHO4ImObLQileB9Fx34O4YehVQ7xbhedC7Oi/fwS+BuwTsZ87YBronUvjT7nJ4Adaln3\n7Vn25Zlo2nOp1xnLjAPm1HA8B6Xt42vATlmW3xK4l/DdWQy8B1wRTbskSzwrgU2znX/RuG7AP6Pz\nYBHhvD9kVc6/Go7VjsBjhHN8YbTsrmnTs8Ve4/c3/bhljG8TjX887btQAdyYS5xp67kVeCd6/Uhq\nfXksPxB4NTqm3wPPA/tnzDOEcN1aCnwOXA90yJjnOcL1Ybvo9SLgQ6LrbfTZvJJ2LuyXsfyl0f6n\nzpsFhOv2NcDqcb5nVF7j9oj2cQnh+/XbLMehQ7StT6N1fgicD1i+3wGyf99Wpk0/hnDt+V+0n2+T\ncd3L6zOMu2BjHKhMMvYB1s0Y1kmbbzVgGjAbaBuNOzD6MP6YNl/qovEW8Abhx2J4dKLOTD/5yJ5k\nfApcR/hP4RzChWklcHDGfHPInmRMA56MvmRXEn6kJmQs+9tonY9E850XncjfEV00o/kOAFYAb0ax\nXAbMA97JjLuGC8FKoE/G+E2jOEvTxr1K+HEaGsXzWDTPGVn2+QPCF3l49KXZq4bj0SeadzhwKuEC\n/inhotQxy+c1LfriDSVcDH4ApmZs/xeEL+hsQrH5acAoYEraPNukHaPzos/x2ehYHF7HMUtdEGYQ\nLhi/i+LpTDgfPwOuirZbDrwbxdMrWr5NNK2CkOQcGw3bRtMvIe3CEY0bF83/D2AwlRebiTl8d8oI\nF+cLgVMIye2iLMftLsLF8kpCQnke8CAwoJZ170q4vVMRHeNjiS7y0fHMlmRU+T6lHc9pwPvRdssJ\nScQnQIu0eXsRLp5fR9s9lXBr6c1o+nbRfqwEzk47tmvUcP5tAHxFSAj/RPj+TCd8nw5Pmy/n86+G\n47Q1IbH4DPhjdM58FB3vnaN5to3WWwHcGcV9WA7nYWaS0Su1juj98Oh47JHH9bYV4Tv4x7TrxDJg\ngxyXTyVML0bn31mEUpA/p81zaTTP44TryWjCdfCVjM/82ei4fQz8NZr3nSieo4AvgGHR5/3fKO62\nWWJ5i3A+n0Eo2akAxsX5nkXn0cxo25dH6/x3dN70TJtvjWi7XxOuy4Oida4ERub7HSB836ZEyw+I\nzpFj0657FdH0wdEwGvhHrp97tc8x7oKNcaDyxznbsDhj3m0IF/WbCFnkZ9GJW5I2T+qi8SnQJm38\nkdH4s9LGZUsysmXAbwNPZjkZsyUZj2fMNyL60rSP3reNviw3Zsy3PuHHcWzauDeifWyXNm4/6vhP\nKJqvPdEPS8b430VfmM417XM07jHgwyz7vJKM/1pqOB4ts8yzaRTThVk+rxlUvQCdHW1r6+h9CSG5\nmJU6ljXs91PRcVstY/y/gPfqOGapC8I8qpcGWJZ1rkkoPbglbdy6pJVeZMxfJcmg8kdjbMZ8V0b7\nvncd8Wb73I4m44cn2p9rY3w3U5/NERnj800yvgbWTBvfL4rxkLRxzxMSgo1riaectNKLOs6/UdG8\nfdPGtY3On1lZ9rHW86+WmB6IzukuaeM6EhKmZ7Mci7La1pcx77DofNoginM6ackyMDF6v2Zd60xb\nd/9ome7R+3aEf8Dq/K+YUEdnBfDPWuZZj3CNfjRj/JBouydknEcrgaPSxm0R7ftyqv6nn/qhPT7j\n+1QB3J+xreuj9aaS+5y/Z1Re43bP2Kcq19Los/lf6jimjf8z4Xq/cYzvwHVk/BOSdi7Py/f7W9vQ\nrFqXRJyQMe6fMRxcZSb3/xBOrEGErG4dwkmbrcLXHe6+OG3Z+wg/CIfUGoj7j6nXZrYWodj/RaB3\njvtxc8a4FwmJSpfo/QGEBOkfZrZuaoiWfZVQooOZdSTcax3n7j+kxfc04T/o2gNxX0hIFI7KmHQU\n8Iq7f5Y2b/o+rxnF8wLQ3czaZyw/x92fymH7P90jN7MSM1uHcDF7n+zH8jZ3X5n2/kXCD3v36P2O\nQFfgmmjfqjGztQnH759Ah4zj+wSwuZl1qit24D53/z5jf9zdV0TbsWhbrQhFmLmcG9kcQvjcR2WM\nH0HY90NrWzjjc1s92s9Xo2XTY5oP7JrjvhfDP9w9vaJmlc/WzNYDfgb8zd0/L9A2DybcEpyaGuGh\ncuPNQFcz2zpj/rrOv2rMrITw4/eAu//U6sbdvyJU2tzTzNqtwj78CfiGUCLzDOH2z/nu/lA0PdXK\nJuv3oQbHAq+7++wo1h8IJarH5bDsrwnH5LJa5tmfcGvxmozxt0RxZp7TP7j7vak37v4B4Xyd6e7p\nrW9ejf5mfh4OjMkYd10UZ+pafyj5fc/edfeX02L6lnDdSt/2kYRzZEHGdeZpQql7ZuuyWr8DdZgP\ntDWzA3OYNyfNteLnvz23ip9XEe5P7Qxc4O7v1zBftlYpHxF+qGpkZr8kFD/vAKyeNinXmuv/zXg/\nL/q7dvR3M8LJ9WyWZZ3KimGppCTbfrxP+NGtyz3A4VFt9FfMrDvhNsbQ9JnMbA/CBW03QpF/ejwd\nqHoRm5PDdjEzIzR7PINwcWyRts5vsyxS13HrES37n1o2mzq2lwNXZJnuhP8Kv6wj/I+zjTSzEwhF\nxFsRLqQps+tYX01S/+VU+Yzdfa6ZzafyHMgqSnQuJZRebJA2KfW5pZxPKC7+r5lNI9RpGe/uOX2W\nBVDls3X3+eH0+OmzTV1oa/ts89WFUMqZaWba9PRkva7zL5v1Cd+XD2rYTgmwSdo283UzIWGuIPzQ\n/MerVnBN/Wi1p3prm2rMrAPhh/c6M+uRNull4Agz28zda2vN1z2Kpbb9SZ2zVY6Juy83s9lUP6c/\no7oFVD9n/pdxzqTLjHlWFGfX6H3qFnGu37NPs2xjXsa2Nyfcwvsmy7yp60y6ur4DtbmBUGfoUTP7\ngvAP073uPiWHZbNqrklGrnoQPmAIH3LBmNnPCBVEnyP8OH5JKLY7mXCfLBcrs4wzKmsOlxBOwoGE\nCniZVuQecZ0mE4r5jiJccFNF6ff9FFhIPJ4iXDhKCV+GZYTs/lyq99uyJMdtX0j4j+dWQtHi94Qv\n+ugs64Tsxw3yq3GdWu/VhJKubHJpEl1tH81sIOF2wP2EYtavCTFfQG7/jdTGYy73T0JieCXh/vAP\nhGMwhbRj7O7/NLMXCP+JHkC4L/x7M/t1zAtVTfG2qGF8IT7bYmuIMX7o7s/UMv296O92hMrqdTmK\n8I9TOeEcSOeE0oyiN5/MUNNxX5XPo6bzM9fvWS7bLiHUvfu/GmLKTDxj74+7f2NmOxDqIB4cDSeZ\n2R3uflJdy2ejJKMG0X/H4wiZ7ijgQjO7zzOabkWydUKzGeFiXJMjCD8wB6aKxqPtnhI76CD95J5F\nOLG+qeMCkip+zbYfW+a0UffFZvYw8BszKydcZF6MinNT+hGK/fulF1Wb2X65bKMW/Qn37U9LHxnd\ngsqW/dclddy2JRQdZ5MqUVhex7GNoz/hXv6R6SPNLLPoOJ+E4RPCxWpzQulUap0bAGtRS8dX0XHc\nF7jI3Yenjd8s2/zuPhcYC4yNbk+8QUgE4yQZ8wilU5lqLXmpRepz27aO+fI9ttm+Jz3Tpq+qbwi3\nAGvaTgXVS0gKaTKhsulAcksyjiVUrMyWSAyOpteWZMwinK9bE+qpZZM6rluSViJoZi0J58yTOcSZ\nr82p+nluRogzVVIX+3tWi1mEunLZSqTjqvH8jn6PHokGzOxG4DQzuzx16ysfzbFORq7KCf+5DSK0\nLngZuDG635/p+PT7oWb2G6AToai4JisJH/RPiZ6ZdQUOX9XA00whFG1eEPWoWEX0A5C6r/smcEJ6\nvQgz+wXhS56re4CNCLX1tyfUrk6XyrB/Ou+iYtUT89hGNivJyNKjz2DjmOubTrhonBvFV427f0Mo\nhTo9qtNSRerYxlTtPxEz2xXomzE6VQ9orRzW+SjhGJ2bMb6ccB4+kkM8mdeLUtIuVlF9mCo9ZEb3\nmL+g6u3AfMwCtoruQae2sz2h2V/eonheAE42s01qmTXVYVSux3aX6DNKxdiW0PpnjrvXWa+pLlFd\nsCcItyQ3TdvOhoSSzxfT61MVmru/QmjBcaqZVbtGmVkrM7sqet2ZUE/gHne/P3MglNJtZmY717LJ\nBwnn1sXRP3zZPEUo/R2aMf5UQh2Sh/PYxVwYcGbGuKGEOB+P3q/K96wm9wJ9zeyAagGZdYj618nX\nomj5Kt/XGn7f3on+xvoON8eSDAMOMbOeWaa97O5zommXAbe7+6MAZnYi4Yf4RsKtgHTfA/8ys9sJ\ntb3PIRRh3VpLHI8Q7rlPMbO7CZ0RDSE0Z8zlmQ01ffF+Gu/uC83sDEIXw9PN7B+E/4g2Jdyi+BeV\nX9A/Er6UL1l4LsO6hCZjMwi1wnPxKKEY/WrCrZj7M6Y/QbgoPGxmNxHu755KuJVT7Yc6Dw8DF0Vx\nv0wo0j2O8AOVN3f36LhNAt6MPtcvCfUjtnb3VCXhMwmVqt4xs1sI/yVvSEgGNia3uiw17c8RZvYg\n4TzpDpxOqEfw02fh7kvN7F3gaDP7kHAezogqLWfu09tmdgfhP5K1CS0sdgWOJ9SYf76W47EwugVy\nvpm1IvRFcADhPnT6edge+MzM7qPylsovgJ0I53oct0XLPmFmfyMc39MJ52W2Lr9zMZTwuU03s5sJ\nCWU3Qu371Gc2jbBvf46+N8uBSe6e7RbeXwk/9I+b2bWEz+FEQmnLETFjzGYYobLjS2Z2AyH5O41Q\nOnh+AbdTk+MJ/7hMjEotnyb8WG1OqLvWkdCiLFWxc3IN63mUEPtxhCab1bj7LDMbTtjnF83sfkIf\nRTsDn7v7he7+rZn9hZCIPE74vm5FuP38GqEZcqF1M7OHCEnF7tE+3Onu70Rxx/6e1eIq4DDCdXMc\n4dxsS/idOILwPfy+poVrkDq/rzOzKYSWJvcAt0aJxjOEOixdCb8Db7h7vPo+hWyq0tAHKvvJqGk4\nnvDf2quE4rf2GcunmpodGb3fO3p/FKHyX12dcc3KGHci4V7nYsIPyPFk799gNqE2fOZ+9M6YLxXP\nXhnj96KyU6lFhATob8COGfOlOuNaTMheD88Wdx3H+O9RDFk73SEkN29EccwiZPgnktFcMNrnh2pY\nR+bxaEWoK/BZdPyfB3YhfFGeznJ8MptJdkl9/hnj+xIuJvMJJUJvUL0/j67RMfqc0Jzu0+jz/1Ud\nxym1zdIapv8+2s9UJ20H13AO7Uq4oC6J1pfqjOsSYEXGvCWEi3aqk56PCRVXqzUBzhJPJyo7Hfse\nmED4wV9JuI0CoYLqXwklQaljNh04LYf1Z/1somkDCMn3EsLFcf/MY1Hb8UyPMW1cz7T9WUSomHlJ\nxjwXRJ/n8vTzM/P8SzsP7qFqZ1wH5bKPNZ1/NRyn7Qnf5VRnXE8Cu+RzbsWdN5p/dUIJ1itRDEsI\n17BRQNdonreou9n7M4TrZUkd850Qnf+LCZW4nwH2zZjnDML1cymh1Ow6MpraEiq/v5Vl/VmvM9Ex\nGZ32/pJoXKozrvlUdsaV2clWTt+zWrb9LGnXrWhcG8JvzPvRMZ9LSJTPpbL/i5y/A1GM1xBaE60g\n+s0hJC2PRZ/NEkICPoYc+zbJNli0YonBzPYmnBBHeigGFBGRJsbCw8QuBtb3jCbnUrsGUyfDzM40\nszkWnur3Sm3366zyaaorrfLpgYtrml9ERETqX4NIMszsaEJnJZcQ7mO/RairUFvluQWEe4CpIW5t\ncxERESmCBpFkEO7x3eTu4939PUITp8WEPiNq4u7+jbt/HQ1xmioWgu43iYiIZJF4khG1ae5DqKkM\nhOyB0Dwps8leunZm9rGZfWpmD1r1rnuLzt2fd/cWqo8hItJ0ufufomu96mPkKfEkg/BAmBZU75Gy\ntmaN7xNKOQ4jNCEqAV42s42KFaSIiIjkp1H2k+GhY5ifnhVgZlMJXVWfTqjXUU3Umc+BhOZES4sf\npYiISJPRmtBUe4q7f5frQg0hyfiW0IZ3w4zxGxLa8NbJ3VeY2RuELl5rciDF6ZxFRESkuTiO8OTf\nnCSeZHh4Yt40YD9Cj22p54bsB1ybyzosPAZ5O2rvsvVjgDvvvJOePbN19inFUFpayqhRmU89lmLS\nMa9/Oub1T8e8fs2cOZOBAwdCDU+OrkniSUZkJDAuSjZeI7Q2aUN4QBlmNh74zN0viN5fRLhd8hHh\n2QLnE7rKrq0b76UAPXv2pHfv3sXZC6mmQ4cOOt71TMe8/umY1z8d88TkVd2gQSQZ7n5v1CfGZYTb\nJG8Snk6aapbamaqPJV8buJlQMXQeoavhvlHzVxEREWkAGkSSAeDuNwA31DBt34z3ZcR/4JKIiIjU\ng4bQhFVERESaICUZUlQDBgxIOoRmR8e8/umY1z8d88ah2TyF1cx6A9OmTZumykIiIiJ5mD59On36\n9AHo4+7Tc11OJRkiIiJSFEoyREREpCiUZIiIiEhRKMkQERGRolCSISIiIkWhJENERESKQkmGiIiI\nFIWSDBERESkKJRkiIiJSFEoyREREpCiUZIiIiEhRKMkQERGRolCSISIiIkXR7JKMX/5yMEOHXsLC\nhQuTDkVERKRJa3ZJxpdf3siYMX3p27e/Eg0REZEianZJBhgVFQcxc2Ypw4aNSDoYERGRJqsZJhlB\nRcVBTJr0UtJhiIiINFnNNskAY/78Nsyf70kHIiIi0iQ14yTDmT9/ER07GkcdBZMnw/LlScckIiLS\ndDTbJKOk5HFOPnlPLr8c3nsPDjsMNtoIzj4bXnsNXAUcIiIiq6QZJhlOSclj9Ow5imuuKed3v4O3\n34Y334QTToCJE2HXXWGrreCKK+Djj5OOV0REpHFqdklGp05DOOusV5k6dSLt27f/afz228PVV8N/\n/wtTpsAuu8Bf/gLdusHPfgY33wzz5iUYuIiISCPT7JKMhx++kdGjL62SYKRr0QIOOAD+/neYOzf8\nbdMGzjgDOnaEI4+Ehx6CZcvqOXAREZFGptklGflo1w4GDgwlG//9L/z5z/DRR/CrX4X6G2eeCa+8\novobIiIi2SjJyNFGG0F5eai78fbbcPLJ8OCD0LcvbLklXHYZzJ6ddJQiIiINh5KMGLbbDq68Ej79\nFJ58EnbfHa66Cnr0gD33hJtuUv0NERERJRmroEUL2H9/GDcOvvoK7roL2reHIUNC/Y0jjoAHHoAf\nf0w6UhERkfqnJKNA2raFY4+Fxx6Dzz+Hv/41NH894gjo1ClUHH35ZdXfEBGR5kNJRhF07AilpTB9\nOsyYAaedBg8/DHvsAZtvDpdeGiqQioiINGVKMopsm21CqcYnn8DTT8Nee8HIkSHZ2H13uPFG+O67\npKMUEREpPCUZ9aSkBPbdF267LdTfmDAB1l47dGPeqVNoFjtxoupviIhI06EkIwFt2sAxx8Ajj4T6\nG1ddBZ99Fjr66tgRTj8d/vUv1d8QEZHGTUlGwjbcEM45B15/Hd59N1QQffzx0JV5jx5w8cXw4YdJ\nRykiIpI/JRkNSM+eoVfROXPguefC7ZXRo2GLLWC33WDMGPj226SjFBERyY2SjAaopAT23htuvTXU\n37jnHlh/fTj33FB/4/DD4b77YOnSpCMVERGpmZKMBm6NNeCoo2DyZPjii9Ay5csv4Te/CfU3Bg2C\nF16AioqkIxUREalKSUYjsv76oTXKa6/BzJlw1lmhW/O994bu3WHYMHj//aSjFBERCZRkNFJbbQVX\nXBEeyvbCC+Hx9NdfH8bvsgtcdx18803SUYqISHOmJKORKykJLVFuvjnU3/jnP0O9jbKy8OTYfv3g\n3nthyZKkIxURkeZGSUYT0rp16GvjoYdCvY1rrgmlGUcfHepvnHJKaLWi+hsiIlIflGQ0UeutB2ee\nCa+8EuppnHMOPPss7LMPdOsGF1wQ6nWIiIgUi5KMZmCLLeCyy2DWrNCT6MEHh2embL017LRT6Ivj\n66+TjlJERJoaJRnNiFl4EuzYsaH+xsSJsMkm8Lvfhfobhx4K//gHLF6cdKQiItIUKMloplZfHY44\nAh54INTfuO46mDcPBgwI9TdOOgmeeUb1N0REJD4lGcK664Znprz8Mnz0EZSXw4svwn77QZcu8Ic/\nwH/+k3SUIiLS2CjJkCp69IBLLgkPZXv55dAE9pZbYNttoXdvGDUq3GrJletRsiIizZaSDMnKDPr2\nhRtuCLdTHnggtEr5wx9g441D5dG7785ef2PhwoUMHXoJ3brtzyab/Ipu3fZn6NBLWLhwYf3viIiI\nJEZJhtSpVSv41a9CRdGvvgqJx8KFcNxx4VH1J54ITz8NK1eGBKNv3/6MGdOXjz9+ks8/f4iPP36S\nMWP60rdvfyUaIiLNSINJMszsTDObY2ZLzOwVM9s5x+WOMbMKM7u/2DEKrL02nH56aAo7axacf364\nrbL//qH+xs9/fjUzZ5ZRUXEQYNFSRkXFQcycWcqwYSOSDF9EROpRg0gyzOxoYARwCbAj8BYwxczW\nq2O5rsBVwAtFDlGy6N4dLroodPb1yiuhtOPNN1+iouLArPNXVBzEpEkv1XOUIiKSlAaRZAClwE3u\nPt7d3wMGA4uBk2tawMxKgDuBi4E59RKlZGUGu+4K113ndOzYlsoSjGpzsnx5G1UGFRFpJhJPMsys\nJdAHeDo1zsOv0FNA31oWvQSY6+63FzdCyZWZ0arVIqCmJMJp2XIRZjUlISIi0pQknmQA6wEtgLkZ\n4+cCHbMtYGZ7AicBpxY3NMlXv357UFIyJeu0kpLHOeywPes5IhERSUpDSDLyYmbtgPHAIHefl3Q8\nUtXw4efRs+dISkoeo7JEwykpeYyePUdxxRXlSYYnIiL1aLWkAwC+BVYCG2aM3xDI1u1TD6ALMNkq\ny91LAMxsGbClu9dYR6O0tJQOHTpUGTdgwAAGDBgQL3qpon379kydOpFhw0YwadJI5s9vw/z5iznu\nuD0YM2Yi7du3TzpEERGpxYQJE5gwYUKVcQsWLIi1LmsIlfDM7BXgVXc/J3pvwKfAte5+Vca8rYDN\nMlYxHGgHDAU+dPcVWbbRG5g2bdo0evfuXYS9kGyWLnW6djV+/evw5FcREWl8pk+fTp8+fQD6uPv0\nXJdrKLdLRgKDzOx4M9sKGAu0AcYBmNl4M/szgLsvc/d30wdgPrDQ3WdmSzAkOa1bG2edBePGwbff\nJh2NiIjUpwaRZLj7vcB5wGXAG0Av4EB3/yaapTM1VAKVhm/w4NDMdezYpCMREZH61BDqZADg7jcA\nN9Qwbd86lj2pKEFJQay3HpxwAlx/PZx3HrRunXREIiJSHxpESYY0faWl8PXXkFGXSEREmjAlGVIv\nttgiPDZ+5EhoAHWNRUSkHijJkHpTVgYzZsATTyQdiYiI1AclGVJv9toL+vQJpRkiItL0KcmQemMW\nSjOeeALeeSfpaEREpNhiJRlm1sLM+pvZsGj4tZm1KHRw0vT85jfQuTOMGpV0JCIiUmx5Jxlmthnw\nLuH5IUdEw53Af8ysR2HDk6amZUsYOhTuugu+ytZpvIiINBlxSjKuBWYDm7h7b3fvDWwKzImmidRq\n0CBo1QrGjEk6EhERKaY4ScbewPnu/n1qhLt/B/whmiZSq7XWglNPDc8yWbw46WhERKRY4iQZPwLZ\nHqXZDli2auFIczF0KMybB3fckXQkIiJSLHGSjIeBm81sV6u0G+GhZpMKG540Vd26Qf/+oQJoRUXS\n0YiISDHESTKGArOAqcDSaHgJ+Ag4p3ChSVNXVgYffggPP5x0JCIiUgx5PyDN3ecDh5vZ5sBW0eiZ\n7v5RQSMI0KTJAAAgAElEQVSTJm+33WD33UPnXIcdlnQ0IiJSaLGfwuruHwIfFjAWaYbKy8Ntk2nT\nQm+gIiLSdOSUZJjZSOAid18Uva6Ru5cVJDJpFg4/HLp3hxEj4O67k45GREQKKdeSjB2BlmmvRQqi\nRQs499zwKPj/+z/YZJOkIxIRkULJKclw932yvRYphJNOgosvhmuvhauuSjoaEREplDjdit9mZtX6\nyTCztmZ2W2HCkuakXTs4/XS4+WZYuDDpaEREpFDiNGE9AVgjy/g1gONXLRxprs46K/T++be/JR2J\niIgUSs5JhpmtaWYdAAPaR+9Tw9rAIcDXxQpUmrbOneGYY2D0aFixIuloRESkEPIpyZgPfA848AEw\nL234FrgN0COvJLayMvj4Y3jggaQjERGRQsinn4x9CKUYzwD9CQlHyjLgE3f/ooCxSTOz446wzz6h\nOeuRR4JZ0hGJiMiqyDnJcPfnAcysG/Bfd9cTJ6Tgysvhl7+EqVNDb6AiItJ4xelW/BMAM2sDbAq0\nypj+dmFCk+bo4INhyy1DaYaSDBGRxi3vJMPM1gduBw6uYZYWqxSRNGslJaFuxuDBMGsW9OiRdEQi\nIhJXnCas1wBrAbsCS4CDCM1aPwT0mCtZZb/9Lay7bmhpIiIijVecJGNfoMzdXwcqCBU+7wTOB/5Y\nyOCkeVpjDRgyBG67DebNSzoaERGJK06S0ZbK/jDmAetHr98BehciKJEhQ0J/GTffnHQkIiISV5wk\n431gy+j1W8DpZrYxMBj4slCBSfO24YYwcGB4nsmyZUlHIyIiccRJMkYDnaLXfyJUAP0UGApcUKC4\nRCgthS++gHvvTToSERGJI+8kw93vdPdx0etpQBdgZ2ATd7+nsOFJc7bNNnDQQaE5q3vS0YiISL7y\nSjLMrKWZzTKznqlx7r7Y3ae7+7eFD0+au/JyePNNeO65pCMREZF85ZVkuPtyoHWRYhGpZr/9oFev\nUJohIiKNS5w6GWOA35tZ3h15ieTLLHTO9cgjMHNm0tGIiEg+4iQZOwNHAJ+a2RQzuz99KHB8IgwY\nAJ06wTXXJB2JiIjkI06SMR+YCEwBvgAWZAwiBdWqFZx1FowfD998k3Q0IiKSqzgPSDupGIGI1Gbw\nYBg+HG68ES6+OOloREQkF3FKMkTq3TrrwIknwpgxsHRp0tGIiEgulGRIo1FaGm6X3HVX0pGIiEgu\nlGRIo7HZZnD44TBypDrnEhFpDJRkSKNSVgbvvguPP550JCIiUpdVSjLMTB1zSb3ac0/YeedQmiEi\nIg1b3kmGmZWY2UVm9jnwg5l1j8ZfbmanFDxCkTSpzrmeegreeivpaEREpDZxSjKGAScC5wPpD+Ge\nAZxagJhEanXkkbDppjBqVNKRiIhIbeIkGccDp7n7XcDKtPFvAVsVJCqRWqy2GgwdCnffHR4FLyIi\nDVOcJGNj4KMa1tVy1cIRyc2pp0Lr1nD99UlHIiIiNYmTZLwL/CzL+COBN1YtHJHcdOgAgwbB2LGw\naFHS0YiISDZxkozLgOvN7PfR8keY2S3AhdE0kXoxdCj8738wblzSkYiISDZ5Jxnu/hDQD9gfWERI\nLHoC/dz9ycKGJ1KzLl1CJdBrroGVK+ueX0RE6lesfjLc/UV3/4W7b+Dubdx9T3d/otDBidSlrAw+\n+ggmT046EhERyRSnn4ydzWzXLON3NbOdChOWSG522SV00KXOuUREGp44JRljgE2yjN84mhaLmZ1p\nZnPMbImZvWJmO9cy76/N7N9mNs/MfjCzN8xsYNxtS+NWVgYvvgj//nfSkYiISLo4ScbWwPQs49+I\npuXNzI4GRgCXADsS+tyYYmbr1bDId8AVwG7AdsDtwO1m9os425fG7bDDoEcPGDEi6UhERCRdnCTj\nR2DDLOM7AStixlEK3OTu4939PWAwsBg4OdvM7v6Cuz/k7u+7+xx3vxZ4G9gz5valEWvRIjwG/r77\n4JNPko5GRERS4iQZTwB/MbMOqRFmthbwZyDv1iVm1hLoAzydGufuDjwF9M1xHfsBWwDP57t9aRpO\nPBHWXBOuvTbpSEREJCVOknEeoU7GJ2b2rJk9C8wBOgLlMda3HtACmJsxfm60zqzMbE0zW2hmy4DJ\nwNnu/kyM7UsT0LYtDB4Mt9wS+s4QEZHkrZbvAu7+uZn1Ao4DtgeWEOpETHD35QWOrzYLo+23A/YD\nRpnZbHd/obaFSktL6dChQ5VxAwYMYMCAAUULVOrHWWfB1VfDrbeGyqAiIpK/CRMmMGHChCrjFixY\nEGtdFu5MJCe6XbIY6O/uk9LGjwM6uPuvc1zPLUBndz+4hum9gWnTpk2jd+/eqx64NEgnnADPPQez\nZoUHqYmIyKqbPn06ffr0Aejj7tkaf2QV6zJsZpsD+wAbkHHLxd3z6lrc3Zeb2TRCacSkaP0Wvc/n\nDnsJsHo+25amp6wMxo+HiRPh6KOTjkZEpHnLO8kws0HAjcC3wFdAelGIE+/5JSOBcVGy8RqhtUkb\nYFy0zfHAZ+5+QfT+D8DrwCxCYnEoMJDQKkWase23h/32C81ZjzoKzJKOSESk+YpTkjEMuNDd/69Q\nQbj7vVGfGJcRmse+CRzo7t9Es3SmavPYtoSOvzoT6oS8Bxzn7vcVKiZpvMrL4ZBD4KWXQm+gIiKS\njDhJxtrAPwsdiLvfANxQw7R9M95fBFxU6BikaTjwQOjZM5RmKMkQEUlOnCas/wQOKHQgIoVSUhLq\nZjz0UHh4moiIJCNOScZHwOVmthvwDlCl2WrU+6ZIogYOhAsuCI+Bv/76pKMREWme4iQZpwE/AHtH\nQzonvxYhIkXRujWceSZceSVcdhmss07SEYmIND953y5x9261DN2LEaRIHEOGwMqVcNNNSUciItI8\nxamTIdIorL8+HH88XHcdLFuWdDQiIs1P3M64OgOHAZsCrdKnubs6dJYGo7Q0PM/kH/8ICYeIiNSf\nOJ1xpXrmnA1sBcwAugIG5NzVqEh96Nkz9JkxYgT89rfqnEtEpD7FuV3yF+Bqd98OWAr0JzyV9XmK\n0H+GyKoqL4e334Zn9IxeEZF6FSfJ6AmMj16vANZw9x+Ai4HfFyowkULZZ5/Q3fiIEUlHIiLSvMRJ\nMhZRWQ/jS6BH2rT1VjkikQIzC6UZjz0G776bdDQiIs1HnCTjFSDVWfOjwAgzuxC4LZom0uAcfTRs\ntBGMGpV0JCIizUecJKMMeDV6fQnwNHA08DFwSmHCEimsVq3g7LPh73+HuXOTjkZEpHmI0xnXbHd/\nO3q9yN0Hu3svd+/v7p8UPkSRwjj9dGjRAm68MelIRESaB3XGJc3G2mvDySfDmDGwZEnS0YiINH05\nJRlm9r2ZrRe9nhe9zzoUN1yRVXPuufDdd3DnnUlHIiLS9OXaGVcpsDB6fW6RYhEpuh494Fe/gpEj\n4ZRTwmPhRUSkOHJKMtz9DgAzW43wpNUp7q7qc9IolZfDnnuGJq2HHpp0NCIiTVde/8e5+wpgLNC6\nOOGIFN/uu8Ouu4bSDBERKZ44hcWvATsWOhCR+mIGZWWhm/E330w6GhGRpivOU1hvIHTA1RmYRugB\n9Cep5q0iDdkRR0CXLqE0Y/z4uucXEZH8xSnJ+AfQDbgWeAl4E3gj7a9Ig7faanDOOTBhAnz+edLR\niIg0TXGSjG5Zhu5pf0UahVNOgTZt4Lrrko5ERKRpitPj5ye1DcUIUqQY1lwTTjsNbroJfvgh6WhE\nRJqe2L0EmNnWZnaQmR2WPhQyOJFiO/tsWLgQbr896UhERJqevCt+mll34AFgO0KfGRZN8uhvi8KE\nJlJ8m24KRx0F11wDQ4aEZ5uIiEhhxCnJGA3MATYAFgPbAHsBrwM/L1hkIvWkrAxmz4aHHko6EhGR\npiVOktEXuNjdvwUqgAp3/xfwR0KLE5FGZaedYK+91DmXiEihxUkyWlD5HJNvgY2i158AWxYiKJH6\nVlYGL70Er76adCQiIk1HnCRjBrB99PpV4Hwz2wO4GJhdqMBE6lO/frD55jBiRNKRiIg0HXGSjCvS\nlruY0D/Gi8AhwNACxSVSr0pKoLQUJk6Ejz9OOhoRkaYhTj8ZU9z9/uj1R+6+FbAesIG7P1PoAEXq\nywknwFprwejRSUciItI05J1kmNlAM2ubPs7dv3d3r2kZkcagTRs44wy49VZYsCDpaEREGr84t0tG\nAXPN7G4zO8TM1LOANBlnngnLlsEttyQdiYhI4xcnyegEHEPofOte4EszG2Nmuxc0MpEEdOoExx4b\nbpksX550NCIijVucOhkr3P1hdz+O0CFXKdAVeNbMZhU4PpF6V1oKn30G992XdCQiIo1b7GeXALj7\nYmAK8BjwISHZEGnUevWCX/wiNGdVTSMRkfhiJRlm1sbMjjOzR4HPgXMJzzPZppDBiSSlvBymTYMX\nX0w6EhGRxitO65J/AF8TKoDOBn7u7pu5+0Xu/l6hAxRJwgEHwDbbqHMuEZFVEackYyVwFNDJ3c9y\n96kFjkkkcWahq/HJk+GDD5KORkSkcYpT8fM4d3/U3VcWIyCRhuLYY2H99cNj4EVEJH+rVPFTpClr\n3RrOOgvGjYPvvks6GhGRxkdJhkgtBg8OLUzGjk06EhGRxkdJhkgt1l8/PNPk+uvhxx+TjkZEpHFR\nkiFSh9JS+OormDAh6UhERBqXnJIMM1sz16HYAYvUty23hF/+EkaOVOdcIiL5yLUkYz4wL8dBpMkp\nL4d33oGnnko6EhGRxmO1HOfbJ+11V+CvwDgg1UdGX+AE4I+FCkykIdl7b9hxx9A51y9+kXQ0IiKN\nQ05Jhrs/n3ptZhcDZe6efod6kpm9A5wG3FHYEEWSZxZKMwYOhBkzYNttk45IRKThi1Pxsy/wepbx\nrwO7rFo4Ig3XUUfBxhvDqFFJRyIi0jjESTL+CwzKMv7UaJpIk9SyJQwdCnfeGVqbiIhI7eIkGaXA\n2Wb2jpndGg1vA2dH00SarNNOC8nGDTckHYmISMMX59kljwJbAJOBdaJhMrBFNC0WMzvTzOaY2RIz\ne8XMdq5l3lPN7AUz+z4anqxtfpFCWWstOOWUkGQsXpx0NCIiDVuszrjc/b/ufoG7HxENF7p77Fsl\nZnY0MAK4BNgReAuYYmbr1bDI3sDdwM+B3Qi3aZ4ws05xYxDJ1TnnwLx5MH580pGIiDRssZIMM/uZ\nmd1pZi+b2cbRuN+a2Z4x4ygFbnL38e7+HjAYWAycnG1md/+tu49197fd/QNCfZASYL+Y2xfJWffu\n8OtfhwqgFRVJRyMi0nDlnWSYWX9gCrAE6A2sHk3qAFwQY30tgT7A06lx7u7AU4SWLLloC7QEvs93\n+yJxlJfDBx/AI48kHYmISMMVpyRjGDDY3QcBy9PGv0RIOvK1HtACmJsxfi7QMcd1/B/wOSExESm6\nvn3DMHJk0pGIiDRccZKMLYEXsoxfAKy1auHkz8z+ABwF/Mrdl9X39qX5KiuD556D6dOTjkREpGHK\ntVvxdF8BmwEfZ4zfE5gdY33fAiuBDTPGbxhtq0Zmdh5wPrCfu/8nl42VlpbSoUOHKuMGDBjAgAED\ncg5YBEK9jG7dQmnGnXcmHY2ISGFMmDCBCRmPnV6wYEGsdZnn+VhJM/sjMJBQKfNJ4BCgCzAKuNzd\nr8s7CLNXgFfd/ZzovQGfAte6+1U1LHM+4VkpB7j7v3PYRm9g2rRp0+jdO85dHZHqRo+G886DOXOg\nc+ekoxERKY7p06fTp08fgD7unnP5bZzbJX8lNB99GmhHuHVyK6F1SN4JRmQkMMjMjjezrYCxQBvC\nQ9gws/Fm9ufUzGb2e+AyQqLzqZltGA1tY25fJJaTT4a2beHaa5OORESk4YnTGZe7+3BCJ1zbEvqp\nWN/dL4obhLvfC5xHSBzeAHoBB7r7N9EsnalaCXQwoTXJfcAXaUN53BhE4mjfPvQCevPNsHBh0tGI\niDQssfrJAHD3Ze7+rru/5u4/rGog7n6Du3d19zXcva+7v542bV93PzntfTd3b5FluGxV4xDJ19Ch\nsGgR3HZb0pGIiDQscfrJaGtml0cdcX1kZrPTh2IEKdKQde4MRx8N11wDK1cmHY2ISMMRp3XJrYRu\nvf8OfAnkV3NUpAkqK4O77oIHHoAjj0w6GhGRhiFOknEwcKi7v1ToYEQaq9694ec/D81ZlWSIiARx\n6mTMQ913i1RTVgZTp4ZBRETiJRkXAZeZWZtCByPSmB16KGyxBYwYkXQkIiINQ5zbJeVAD2CumX1M\n1eeX4O7q6UqapZISKC2FM8+E2bPD01pFRJqzOEnGgwWPQqSJOP54GDYs9AQ6enTS0YiIJCvvJMPd\n/1SMQESagjZtYMiQUAH0T3+Cter9kYEiIg1H7M64RCS7IUNg+fLQC6iISHOWU5JhZt+b2XrR63nR\n+6xDccMVafg6doSBA8PzTJYtSzoaEZHk5Hq7pBRIPZnh3CLFItJklJaGbsb/+U847rikoxERSUZO\nSYa735HttYhkt+22cOCBoTnrsceCWdIRiYjUv1Wqk2Fmrc1szfShUIGJNHbl5fDGG/D880lHIiKS\njLgPSLvezL4GFhF6AE0fRATYf3/Ybjt1ziUizVeckowrgX2BM4AfgVOBS4AvgOMLF5pI42YWuhp/\n+GF4//2koxERqX9xkox+wBB3nwisAF509yuACwBVcRNJM2BAaG0yalTSkYiI1L84ScY6wOzo9f+i\n9wD/AvYqRFAiTcXqq8NZZ8Edd8A33yQdjYhI/YqTZMwGukWv3wOOil73A+YXIiiRpmTw4HDrZOzY\npCMREalfcZKM24Hto9d/Bc40s6XAKOCqQgUm0lSsuy6ceCJcfz0sXZp0NCIi9SfvJMPdR7n7tdHr\np4CtgGOBHd1dj4QSyeLcc8PtkrvvTjoSEZH6s8rPLnH3T9z9fnd/uxABiTRFW2wB/fqFB6e5Jx2N\niEj9yKnHTzMbmusKU6UcIlJVeTnsvTdMmQIHHZR0NCIixZfPs0ty4YCSDJEsfvYz6NMnlGYoyRCR\n5iDXZ5d0q3suEamNWSjNOPZYePtt6NUr6YhERIprVZ9dYmZ69JNIro48Ejp3VudcItI8xEoyzOwU\nM5sBLAWWmtkMMzu1sKGJND0tW8I558Bdd8GXXyYdjYhIccV5QNplwGhgMvCbaJgMjIqmiUgtBg0K\nPYGOGZN0JCIixRWnJOMMYJC7/9HdJ0XDH4HTgCGFDU+k6enQAU49FW68ERYtSjoaEZHiiZNktARe\nzzJ+Grm3VhFp1s45B+bPD880ERFpquIkGX8nlGZkOg24a9XCEWkeunaF/v1DBdCKiqSjEREpjrgl\nD6eY2QHAK9H7XYFNgfFmNjI1k7uXrWJ8Ik1WeTnsthtMngyHH550NCIihRcnydgWmB697hH9/TYa\ntk2bT50ni9Ri111h991D51xKMkSkKco7yXD3fYoRiEhzVF4ebpu8/jrstFPS0YiIFFacJqzr1zJt\nu1ULR6R5Ofxw6N49lGaIiDQ1cSp+vmNmh2aONLPzgNdWPSSR5qNFi/AY+HvvhU8/TToaEZHCipNk\njAQmmtmNZraGmW1sZk8D5wPHFjY8kabvpJOgfXu4Vo8WFJEmJu8kw92vBPoCPwPejoYfgV7u/kBh\nwxNp+tq1g9NPh1tugf/9L+loREQKJ+4D0j4CZgBdgTWBe9z9q0IFJdLcnH02LF4Mf/tb0pGIiBRO\nnIqfexBKLzYHehE65rrOzO4xs7ULHJ9Is7DxxjBgAIweDStWJB2NiEhhxCnJeAa4B9jN3We6+63A\njoTOuN4pZHAizUlZGXzyCdx/f9KRiIgURpwk4wB3/4O7L0+NcPdZwB7ATQWLTKSZ2WEH2HdfGDEC\nXF3ZiUgTEKfi5/M1jK9w98tXPSSR5qusDF57DV5+OelIRERWXc5Jhpk9amYd0t7/wczWSnu/rpm9\nW+gARZqTgw+GrbYKpRkiIo1dPiUZBwKrp72/AFgn7f1qwJaFCEqkuSopgdJSePBBmDUr6WhERFZN\nPkmG1fFeRArgt7+F9daDa65JOhIRkVUTt58MESmSNdaAIUPgtttg3rykoxERiS+fJMOp/vh21YEX\nKYIhQ2DlSrhJ7bVEpBHL51HvBowzsx+j962BsWa2KHq/evbFRCRfG2wQbptcd11ocdKqVdIRiYjk\nL5+SjDuAr4EF0XAn8EXa+6+B8YUOUKS5Ki2FL76Ae+5JOhIRkXhyLslw95OKGYiIVLX11qFJ64gR\nMHAgmKpai0gjo4qfIg1YWRm89RY8+2zSkYiI5K/BJBlmdqaZzTGzJWb2ipntXMu8W5vZfdH8FWY2\ntD5jFakv++0HvXqpcy4RaZwaRJJhZkcDI4BLCA9bewuYYmbr1bBIG2AW8Hvgy3oJUiQBZlBeDo8+\nCjNnJh2NiEh+GkSSAZQCN7n7eHd/DxgMLAZOzjazu7/u7r9393uBZfUYp0i9O+YY6NQJRo1KOhIR\nkfwknmSYWUugD/B0apy7O/AU0DepuEQailat4OyzYfx4+PrrpKMREcld4kkGsB7QApibMX4u0LH+\nwxFpeE4/HVq0gBtvTDoSEZHcNYQkQ0TqsM46cNJJMGYMLFmSdDQiIrnJp8fPYvkWWAlsmDF+Q+Cr\nQm+stLSUDh06VBk3YMAABgwYUOhNiRTUuefCDTfAXXfBqacmHY2INFUTJkxgwoQJVcYtWLAg1ros\nVH9Ilpm9Arzq7udE7w34FLjW3a+qY9k5wCh3v7aO+XoD06ZNm0bv3r0LFLlI/fr1r+H992HGjPBY\neBGR+jB9+nT69OkD0Mfdp+e6XEO5TI0EBpnZ8Wa2FTCW0Ex1HICZjTezP6dmNrOWZra9me0AtAI2\njt73SCB2kXpTXh6asj7+eNKRiIjUrUEkGVFT1POAy4A3gF7Age7+TTRLZ6pWAt0omm9aNP48YDpw\nS33FLJKEPfaAnXeGkSOTjkREpG4NoU4GAO5+A3BDDdP2zXj/CQ0kQRKpT6nOuY45Bt58E3bYIemI\nRERqph9qkUamf3/YdFN1ziUiDZ+SDJFGZrXV4JxzYMKE8Ch4EZGGSkmGSCN06qnQujVcf33SkYiI\n1ExJhkgjtOaaMGgQjB0LixYlHY2ISHZKMkQaqaFD4X//g9tvTzoSEZHslGSINFJdusCRR8I118DK\nlUlHIyJSnZIMkUasrAxmzYJJk5KORESkOiUZIo3YLrvAnnuqcy4RaZiUZIg0cuXl8K9/wWuvhfcN\n4XlEIiLQgHr8FJF4+vWD7t0XctxxV7NixUssX96Wli0X0a/fHgwffh7t27dPOkQRaaaUZIg0cosX\nL2TJkv58+WUZcClggDNmzBSeeaY/U6dOVKIhIonQ7RKRRu7CC69m7twy4CBCggFgVFQcxMyZpQwb\nNiLB6ESkOVOSIdLITZ78EhUVB2adVlFxEJMmvVTPEYmIBLpdItKIuTvLl7elsgQjk/HFF2044wyn\nVy9j221h221h7bXrM0oRaa6UZIg0YmZGy5aLACd7ouG0bLmIf/3LuPVWWLEijN1445BsbLcdPyUe\nW28Na6xRj8GLSJOnJEOkkevXbw/GjJlCRcVB1aaVlDzOKafsyejRsGwZfPABzJgB77wT/k6cCFdf\nnZoXevSomnhstx1stll48quISL506RBp5IYPP49nnunPzJkeJRqhdUlJyeP07DmKK66YCECrVpXJ\nwzHHVC7/ww/w7ruViceMGXDTTTB3Lj8t17NnZfKR+rvJJmA13aUREUFJhkij1759e6ZOnciwYSOY\nNGkky5e3oWXLxRx22B5ccUXdzVfbtQs9h+6yS9Xx33xTmXSkSj8eeggWLgzT11yzaolH6vV66xVp\nR0Wk0bHm0jugmfUGpk2bNo3evXsnHY5I0bg7VqQiBnf49NOqiceMGTBzZrgdA9CxY/XkY5ttoG3b\nooQkIvVg+vTp9OnTB6CPu0/PdTmVZIg0McVKMMK6w9Nfu3SBQw+tHL9iBXz4YdXE45FHYPTokJgA\ndO9eNfHYbjvYYgto2bJo4YpIwpRkiMgqW221UG+jZ0/4zW8qxy9eHEo50pOP22+HL74I01u2hC23\nrF7ZtEuXUBFVRBo3JRkiUjRt2kCfPmFI9/338J//VCYe77wDjz0G8+eH6e3ahVssmbddNtyw/vdB\nROJTkiEi9W6ddeBnPwtDinso4Uhv5fLGG3DXXbB0aZhn/fWz1/dYc81k9kNEaqckQ0QaBLPQSdjG\nG8NBaV1+rFwJs2dXLfV48kkYMwYqKsI8XbpUb+Wy1Vaw+urJ7IuIBEoyRKRBa9ECNt88DEccUTl+\n6VJ4772q9T3uvju0fkktt8UW1et7dOsWphVSMVv0iDRmSjJEpFFq3Rp22CEM6RYsCPU90pOP0aPh\nu+/C9DXWyF7fo1On/DoXW7hwIRdeeDWTJ7/E8uVtadlyEf367cHw4efV2TeJSHOhJENEmpQOHWD3\n3cOQ4h56ME1PPGbMgHvvDS1gINQTyda52FprVd/GwoUL6du3PzNnllFRcSmpXlbHjJnCM8/0Z+rU\nujtBE2kOlGSISJNnFjoJ69gR9t+/cnxFBXz8cdXE4/nn4eabKx8m17lz9cTj1luvjhKM9OfFGBUV\nBzFzpjNs2AhGj760HvdQpGFSj58iIhmWLYP336/es+mcOak59geepKYn33btegBz5jxZb/GKFJt6\n/BQRKZBWrULJxXbbVR3/ww8wY4Zz8MFtmT+/pgocxrfftmHCBKdPH2OzzdSxmDRfSjJERHLUrh3s\ntpux1lqLmD/fqakkY8mSRRx7rP20zA47wI47Qu/e4e/WW6s7dWkelF+LiOSpX789KCmZknVaScnj\nnHnmnnz7bejP4+KLQ72OJ56Ak08OCUe7drDTTnDaaXDjjfDqq5UVUEWaEpVkiIjkafjw83jmmf7M\nnOlR5c/QuqSk5HF69hzFFVdMpH37UMk0vaLpDz/AW2/B9OmhN9PXXgvPclmxItxS6dmzaonHDjtk\nb4nOEkEAAA5mSURBVN0i0lgoyRARyVP79u2ZOnUiw4aNYNKkkSxf3oaWLRdz2GF7RAlG9uar7drB\nHnuEIeXHHyu7UE8lHxMnwpIlYXr37pVJR+qvnuEijYVal4iIrKJC9/i5YgV88EFl0pH6u2BBmL7R\nRlWTjt69YdNN8+tMTCQfal0iIpKQQncpvtpqoXLo1lvDwIFhnHtoQpuedNx8c+hkDGDttauXeGy+\neeG7UBfJh5IMEZFGwCzcOuneHfr3rxz/5ZdVSzzuuw+uvjpMa9sWtt++auKxzTahia5IfVCSISLS\niHXqBIceGoaU77+HN9+sTD6efhpuuCGUhrRsGXotTS/16NUrJCQihaYkQ0SkiVlnHdh33zCk/PAD\nvP12ZYnH9OkwfjwsXx5atmy5ZdUSjx13DLdgRFaFkgwRkWagXbvqD45btiw8sTb9dsuDD1b22dG1\na/V6Hp06JRK+NFJKMkREmqlWrSpLLVJWrgwtW9IrmI4YAfPnh+kdO1Zv2dK1q1q2SHZKMkRE5Cct\nWoROwXr2hGOPDePc4ZNPqpZ43HZbqHQKocOwVLKSSj623FItW0RJhoiI1MEslFZ07QpHHFE5/quv\nqpZ4PPAAjBwZprVpEyqUppd4bLMNrL56YWIqdN8kUhxKMkREJJaOHeHgg8OQMm9eaNmSSj6eew7G\njoWKitCyZZttqpZ4bL99qC+Si4ULF3LhhVczefJLLF/elpYtF9Gv3x4MH35ejb2sSrKUZIiISMGs\nvTbss08YUhYvDi1b0m+33HVXqHhqBltsUb2C6TrrVF3vwoUL6du3PzNnllFRcSmp58WMGTOFZ57p\nz9SpNXfnLslRkiEiIkXVpg3stlsYUpYtg3ffrXq7ZdIkWLQoTO/SpWrS8cADV0cJxkFpazYqKg5i\n5kxn2LARjB59aX3uluRAzy4REZEGYeXK/2/v/oOsKus4jr8/y68FUuwXiGk6agJGkfgjTaEs1H6M\naz/EtBzTpilKqqksrbFIB/ujzBrFsmkak6k0++k6hPgDtTFRTMnSlsSS0EREF1ZCEIFvfzznxt3r\n7nrv7p57lns/r5kzc89znvOc7z0wu999zvOcBx59dFfSUUpAOjsBZgG3kHowKgUTJ57I0qW3MG4c\njBsHo0d7xstg8tolZma2Wxs2LM1KmTQJzjgjlUXAmjXB9Olj6ezsLWsQa9eOYcqUoJSEDB+eko29\n9uL/iUfl1tuxUnlra12+dkNzkmFmZkOWBPvvL/bcczOdnbuSiO6CffbZzHXXia4uXrJt3Ljr86pV\n3Y9t2tT7tUeOrC4p6evYYM2m2V05yTAzsyHv5JOP5corl1SMyUhaWm7i1FOPY8aM2tvdsQOee67v\nxKSy/Mknu5eVxpH0pLW1f70o5duIEbV/r8FSmtHz618v7tf5TjLMzGzIu+SS81i69EN0dESWaKTZ\nJS0tNzFlyveYP/83/Wp32LA0I2Yg67Rs354Slb4Sk8qyJ57oXr5lS+/tjx5de2JSXr7nnunxUa26\nz+hpA46ouY0hk2RIOhc4D9gbeBD4bETc10f92cDFwAHAI8AFEdG/VMtyc+2113JG6eGq1YXvef35\nnudvjz32YNmy33Dhhd+lvf0yuro6GTfuVbS1Hcv8+cVOXx0+PE25rZx2W4tt27r3qLxcwrJhA6xe\n3b38hRd6b3/s2Np7Uq66qnxGT9VjPbsZErNLJH0YuAb4JLAc+AIwGzgkIp7pof7bgDuB84FFwEez\nz4dFxN97uYZnlxSgra2N9vb2osNoKr7n9ed7Xn++5y/1wgvVP/bp7di2beUtls/oeQDYfWeXfAH4\nUUQsBJA0B3gf8HHg2z3U/xywOCKyF9jyDUknAHOBz9QhXjMzsyFl1CgYPz5t/bV1a0o+Nm4MZs4c\ny/r1A5sH3DKgsweBpBGk9Oi2Ulmk7pVbgWN6Oe2Y7Hi5JX3UNzMzs5fR2ppeFz95shg7djMwsKcd\nhScZwGuAYcC6ivJ1pPEZPdm7xvpmZmZWg5NPPpaWliUDamOoPC6ph1aAjo6OouNoKl1dXTzwQP8G\nDFn/+J7Xn+95/fme52/27FksWvRlHntsFRH/H9Va0yvKhkKS8QywA5hQUT4BeKqXc56qsT6kWSic\neeaZtUdoA5K9itbqyPe8/nzP68/3vF7uLd85ALi72jMLTzIi4kVJ9wPvAtoBJCnbv7yX05b1cPyE\nrLw3S0izUFYDWwcWtZmZWVNpJSUYNT0/GSpTWE8DfgrMYdcU1lOByRGxXtJC4ImI+FpW/xjgDuCr\npCmsZwAXANN7m8JqZmZm9VV4TwZARFwv6TWkl2tNAP4CnBQR67Mq+wLby+ovk/QR4JJsWwWc4gTD\nzMxs6BgSPRlmZmbWeIbCFFYzMzNrQE4yzMzMLBdNkWRIOlfSY5K2SLpH0pFFx9TIJM2Q1C7pP5J2\nSmorOqZGJumrkpZLek7SOkm/k3RI0XE1MklzJD0oqSvb7pb00jXILTeSLsh+vlz28rWtPyTNy+5x\n+VbT2MeGTzKyxde+C8wDDiOt8LokG2hq+RhLGrz7GQb6TlqrxgzgCuCtpBWNRgA3SxpdaFSN7XHS\noozTScsiLAVukDSl0KiaRPaH4idJP88tXw+RJmTsnW3H1XJyww/8lHQPcG9EfD7bF+kHxOUR0dPi\nazaIJO0E3h8RXi6xTrIE+mlgZkTcVXQ8zULSs8B5EXF10bE0MkmvAO4HPg18HVgREV8sNqrGJGke\naeZmv5cub+iejH4uvma2u9uL1IPUWXQgzUBSi6TTgTH0/UJAGxxXAjdGxNKiA2kSb8geff9T0s8k\n7VfLyUPiPRk56mvxtUn1D8csX1lP3feBu/zemHxJmkpKKlqBTcAHImJlsVE1tiyZewtwRNGxNIl7\ngLOBfwATgW8Cf5Q0NSI2V9NAoycZZs3mB8ChwLFFB9IEVgLTgHGkNxQvlDTTiUY+JO1LSqBnRcSL\nRcfTDCKi/BXiD0laDvwbOA2o6rFgoycZ/Vl8zWy3JGkB8F5gRkSsLTqeRhcR24F/ZbsrJB0FfJ40\nVsAG3+HAa4EHsh47SD3VMyXNBUZFow8yLFhEdEl6BDi42nMaekxGlu2WFl8Dui2+VvUqcmZDXZZg\nnAIcHxFrio6nSbUAo4oOooHdCryJ9LhkWrb9GfgZMM0JRv6yQbcHAVX/EdPoPRkAlwE/zVZ6LS2+\nNoa0IJvlQNJYUqZb+mvjQEnTgM6IeLy4yBqTpB+QFglsAzZLKvXcdUWEVxzOgaRvAYuBNcAepBWe\n3w6cWGRcjSwbA9BtnJGkzcCzEdFRTFSNTdJ3gBtJj0heB1xEWkfs2mrbaPgko4rF12zwHQHcTprh\nEKT3lABcA3y8qKAa2BzSfb6jovwcYGHdo2kO40n/nycCXcBfgRM946Hu3HuRr32BXwCvBtYDdwFH\nR8Sz1TbQ8O/JMDMzs2I09JgMMzMzK46TDDMzM8uFkwwzMzPLhZMMMzMzy4WTDDMzM8uFkwwzMzPL\nhZMMMzMzy4WTDDMzM8uFkwwzMzPLhZMMsyYjaYKkKyT9U9JWSf+W1C7pnUXHNlRIulrSb4uOw2x3\n1/Brl5jZLpL2J61A3Al8CXgIGAG8G1gAHFpcdGbWaNyTYdZcfgjsAI6MiN9HxKMR0RER3wOOBpC0\nn6QbJG2S1CXpl5LGlxqQNE/SCknnZL0gmyQtkNQi6SuS1kpaJ+lr5ReWtFPSHEl/kPR81pPyoYo6\nUyXdlh1/RtKPslV9S8evlvQ7SV+S9GRWZ4GkYWV1Rkq6VNITkv4raZmkt5cd/5ikDZJOlPT3LP7F\npdVrJc0DPgacksW8Q9LMQf1XMGsSTjLMmoSkVwInAQt6WgI+Ip6TJKAd2AuYAcwCDgSuq6h+EKn3\n4yTgdOATwCJgH2AmcD4wX9KRFeddDPwKeDPwc+A6SZOy+MYAS4BngcOBU7PrX1HRxvFZTO8AzgLO\nzraSK4G3AqcBb8qut1jSQWV1xpB6cj6afc/XA5dmxy4FrgduIq3cPJHU+2NmNfLjErPmcTAg4B99\n1JkFvBE4ICKeBJB0FvCwpMMj4v6snoBzIuJ5YKWk24FDIuI92fFVks4nJQT3lbV/fURcnX3+hqQT\ngM8Cc0m/8EcBZ2VJUIekucCNks6PiPXZeZ3A3EhLSD8iaRHwLuAnkl5PSjj2i4insvqXSXoPcA5w\nYVY2HPhURKzOvuMC4OsAEbFZ0hZgZNk1zawfnGSYNQ9VUWcy8HgpwQCIiA5JG4EpQCnJWJ0lGCXr\ngO0Vba0DxleU3VOxvwyYVnbtByt6Wf5E6nGdBJR+4T+cJRgla4Gp2eepwDBS8lH+fUcCz5TtP19K\nMMraqIzVzAbISYZZ81gFBOmX+Q0DbOvFiv3opSyPR7J9XecVpGRnOrCzot5/X6aNapIwM6uBx2SY\nNYmI2EAa83CupNGVxyWNAzqA/SS9rqz8UNIYjYcHIYyje9jvyD53ANMqYjuONFC1r0c85VaQejIm\nRMS/Krana4hzW9aOmQ2Akwyz5nIu6ZfnckkflHSwpMmSPgfcHRG3kqa1/lzSYZKOAq4Bbo+IFYNw\n/dnZrJQ3SLoIOJI0dRbSQNCtwDWS3ijpeOByYGG1YyMiYhXwC2ChpA9IOkDSUZIuyMZlVGs18GZJ\nh0h6tST3+pr1g5MMsyYSEY+RHiXcTppF8TfgZtIAzTlZtTZgA3BnduxR0gySmi/XQ9m8rK0HgTOB\n0yNiZRbbFtJslVcBy0kzPG4hDQytxdnAQtL3Wwn8FjgCWFNDGz8m9Z78GXgaeFuNMZgZoO7jp8zM\n8iFpJ/D+iGgvOhYzqw/3ZJiZmVkunGSYWb2429SsyfhxiZmZmeXCPRlmZmaWCycZZmZmlgsnGWZm\nZpYLJxlmZmaWCycZZmZmlgsnGWZmZpYLJxlmZmaWCycZZmZmlgsnGWZmZpaL/wF/sp/G3bSPgAAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9dfaba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(list(pca.explained_variance_ratio_),'-o')\n",
    "plt.title('Explained variance ratio as function of PCA components')\n",
    "plt.ylabel('Explained variance ratio')\n",
    "plt.xlabel('Component')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The first dimension seems to basically represent only the 'fresh'-feature, as this feature has a strong negative projection on the first dimension. The other features have rather weak (mostly negative) projections on the first dimension. That is, the first dimension basically tells us whether the 'fresh'-feature value is high or low, mixed with a little bit of information from the other features.\n",
    "\n",
    "The second dimension is mainly represented by the features 'Grocery', 'Milk' and 'Detergents', in the order of decreasing importance, and has rather low correlation with the other features.\n",
    "\n",
    "There are two main uses of this information. The first use is feature interpretation and hypothesis formation. We could form initial conjectures about the customer segments contained in the data. One conjecture could be that the bulk of customers can be split into customers ordering mainly 'fresh' items and customers mainly ordering 'Grocery', 'Milk' and 'Detergents' from the wholesale distributor. The second use is that, given knowledge of the PCA components, new features can be engineered for further analysis of the problem. These features could be generated by applying an exact PCA-transformation or by using some heuristic based on the feature combinations recovered in PCA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ICA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ -2.74750657e-07   2.48856352e-06   1.16059751e-05  -1.48757046e-06\n",
      "   -2.80619738e-05  -5.69826563e-06]\n",
      " [  3.97798492e-06  -8.91221623e-07  -7.37147347e-07  -6.67184699e-07\n",
      "    2.31765570e-06  -9.68872424e-07]\n",
      " [  3.88102558e-07   2.41511953e-07   5.81460867e-07   5.27946533e-07\n",
      "   -5.00502570e-07  -1.81079160e-05]\n",
      " [  8.64027159e-07   1.38829776e-07  -7.70771017e-07  -1.11465618e-05\n",
      "    5.39460191e-07   5.94210415e-06]\n",
      " [  2.00851727e-07  -1.79251002e-06   7.14897994e-06   3.15277260e-07\n",
      "   -2.47481579e-06  -1.77717934e-06]\n",
      " [  1.61978460e-07   9.81519765e-06  -5.87798423e-06  -3.57537381e-07\n",
      "    3.54446077e-06  -5.97842711e-06]]\n",
      "Independent components scaled with mean\n",
      "[[-0.00329709  0.01442438  0.09228233 -0.00456972 -0.08086039 -0.00868912]\n",
      " [ 0.047737   -0.00516576 -0.00586126 -0.00204955  0.00667831 -0.0014774 ]\n",
      " [ 0.00465735  0.00139987  0.00462336  0.00162182 -0.00144219 -0.02761223]\n",
      " [ 0.01036858  0.00080469 -0.00612861 -0.03424148  0.00155445  0.00906094]\n",
      " [ 0.00241028 -0.01038986  0.05684352  0.00096851 -0.00713116 -0.00270997]\n",
      " [ 0.00194379  0.0568915  -0.04673748 -0.00109833  0.01021334 -0.00911633]]\n"
     ]
    }
   ],
   "source": [
    "# Fit an ICA model to the data\n",
    "# Note: Adjust the data to have center at the origin first!\n",
    "def center_data(data, rescale = 0):\n",
    "    centeredData = data.copy()\n",
    "    for col in centeredData.columns:\n",
    "        centeredData[col] = (centeredData[col] - np.mean(centeredData[col]))/ (1 + rescale * np.std(centeredData[col]))\n",
    "    return centeredData\n",
    "\n",
    "from sklearn.decomposition import FastICA\n",
    "#data_centered = center_data(data)\n",
    "\n",
    "ica = FastICA(n_components=6, whiten=True)\n",
    "ica.fit(center_data(data,0))\n",
    "\n",
    "# Print the independent components\n",
    "print(ica.components_)\n",
    "\n",
    "# Print the independent components (rescaled again)\n",
    "print('Independent components scaled with mean')\n",
    "print(np.multiply(ica.components_,list(np.mean(data))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The first vector [-0.04771087  0.00496636  0.00492989  0.00208307 -0.0059068   0.00159593] again represents mainly the 'fresh'-feature, with a coefficient of -0.0477. The other features have a rather weak projection on the first dimension.\n",
    "\n",
    "The second vector [ 0.00182027  0.0571306  -0.04596392 -0.00113553  0.00928388 -0.00925863] corresponds mainly to the features 'Milk' and 'Grocery', but in different directions. This indicates that, other things equal, high 'Milk'-spending is associated with low 'Grocery'-spending and vice versa.\n",
    "\n",
    "The third vector [ 0.00360762 -0.01311074 -0.09638513  0.00448148  0.08132511  0.00872532] has as strong association with the 'Grocery'- and 'Detergents_Paper'-features, again in opposite directions. This indicates a negative association between these features across the wholesalers customers. \n",
    "\n",
    "The main charactistic of the fourth vector [ 0.00463807  0.00127625  0.00476776  0.00160491 -0.00146026 -0.02758939] is that this vector has a relatively strong negative association with 'delicatessen' (and only rather weak associations with the other features). Even though the coefficient are very low, the vector permits the interpretation that 'delicatessen' are negativly related to the 'fresh'- and 'grocery'-features.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering\n",
    "\n",
    "In this section we will choose either K Means clustering or Gaussian Mixed Models clustering, which implements expectation-maximization. Then we will sample elements from the clusters to understand their significance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choosing a Cluster Type\n",
    "\n",
    "#### K Means Clustering or Gaussian Mixture Models?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before discussing the advantages of K Means vs Gaussian Mixture models, it is helpful to observe that both methods are actually very similar. The main difference is that Gaussian Mixture models make a probabilistic assignment of points to classes depending on some distance metric, whereas K Means makes a deterministic assignment depending on some metric. Now, when the variance of the Gaussian mixtures is very small, this method becomes very similar to K Means, since the assignment probabilities to a specific cluster converge to 0 or 1 for any point in the domain. Because of the probabilistic assignment, Gaussian Mixtures (in contrast to K Means) are often characterized as soft clustering algorithms.\n",
    "\n",
    "An advantage of Gaussian Mixture models is that, if there is some a priori uncertainty about the assignment of a point to a cluster, this uncertainty is inherently reflected in the probabilistic model (soft assignment) and assignment probabilities can be computed for any data point after the model is trained. On the other hand, if a priori the clusters assignments are expected to be deterministic, K Means has advantages. An example would be a data generating process that actually is a mixture of Gaussians. Applying a Gaussian mixture model is more natural given this data generating process. When it comes to processing speed, the EM algorithm with Gaussian mixtures is generally slightly slower than Lloyd's algorithm for K Means, since computing the normal probability (EM) is generally slower than computing the L2-norm (K Means). A disadvantage of both methods is that they can get stuck in local minima (this can be considered as the cost of solving NP-hard problems (global min for k-means) approximately).\n",
    "\n",
    "Since there is no strong indication that the data are generated from a mixture of normals (this assesment may be different given more information about the nature of the spending data) and the goal is to \"hard\"-cluster them (and not assign probabilities), I decided the use the general-purpose k-means algorithm.\n",
    "\n",
    "A decision on the number of clusters will be made by visualizing the final clustering and deciding whether k equals the number of data centers found by visual inspection. Note that many other approaches for this task could be utilized, such as silhoutte analysis (see for example http://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is some starter code to help you visualize some cluster data. The visualization is based on [this demo](http://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_digits.html) from the sklearn documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import clustering modules\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.mixture import GMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[   650.02212207   1585.51909007]\n",
      " [ -4426.80497937   4042.45150884]\n",
      " [ -4841.9987068    2578.762176  ]\n",
      " [   990.34643689  -6279.80599663]\n",
      " [ 10657.99873116  -2159.72581518]\n",
      " [ -2765.96159271   -959.87072713]\n",
      " [  -715.55089221  -2013.00226567]\n",
      " [ -4474.58366697   1429.49697204]\n",
      " [ -6712.09539718  -2205.90915598]\n",
      " [ -4823.63435407  13480.55920489]]\n"
     ]
    }
   ],
   "source": [
    "# First we reduce the data to two dimensions using PCA to capture variation\n",
    "pca = PCA(n_components=2)\n",
    "reduced_data = pca.fit_transform(data)\n",
    "print(reduced_data[:10])  # print upto 10 elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,\n",
      "    n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',\n",
      "    random_state=None, tol=0.0001, verbose=0)\n"
     ]
    }
   ],
   "source": [
    "# Implement your clustering algorithm here, and fit it to the reduced data for visualization\n",
    "# The visualizer below assumes your clustering object is named 'clusters'\n",
    "\n",
    "# TRIED OUT 2,3,4,5,6 CLUSTERS AND CONCLUDED THAT 3 CLUSTERS ARE A SENSIBLE CHOICE BASED ON VISUAL INSPECTION, SINCE \n",
    "# WE OBTAIN ONE CENTRAL CLUSTER AND TWO CLUSTERS THAT SPREAD FAR OUT IN TWO DIRECTIONS.\n",
    "kmeans = KMeans(n_clusters=3)\n",
    "clusters = kmeans.fit(reduced_data)\n",
    "print(clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the decision boundary by building a mesh grid to populate a graph.\n",
    "x_min, x_max = reduced_data[:, 0].min() - 1, reduced_data[:, 0].max() + 1\n",
    "y_min, y_max = reduced_data[:, 1].min() - 1, reduced_data[:, 1].max() + 1\n",
    "hx = (x_max-x_min)/1000.\n",
    "hy = (y_max-y_min)/1000.\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, hx), np.arange(y_min, y_max, hy))\n",
    "\n",
    "# Obtain labels for each point in mesh. Use last trained model.\n",
    "Z = clusters.predict(np.c_[xx.ravel(), yy.ravel()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** K MEANS CENTROIDS ***\n",
      "[[ 24220.71188261  -4364.45560022]\n",
      " [ -1497.13461172  24998.27760147]\n",
      " [ -4106.90273941  -3168.41202086]]\n",
      "*** CENTROIDS TRANSFERED TO ORIGINAL SPACE ***\n",
      "[[ 36135.48457282   6480.2347537    6104.73878363   6844.02904962\n",
      "    1116.07889738   2925.29532792]\n",
      " [  7773.13472743  18509.00259201  26972.9856871    2375.76735765\n",
      "   12025.19401164   2849.79055267]\n",
      " [  8340.22615312   3664.30095668   5275.94880254   2505.50876436\n",
      "    1752.88223845   1064.32056813]]\n"
     ]
    }
   ],
   "source": [
    "# Find the centroids for KMeans or the cluster means for GMM \n",
    "\n",
    "centroids = kmeans.cluster_centers_\n",
    "print('*** K MEANS CENTROIDS ***')\n",
    "print(centroids)\n",
    "\n",
    "# TRANSFORM DATA BACK TO ORIGINAL SPACE FOR ANSWERING 7\n",
    "print('*** CENTROIDS TRANSFERED TO ORIGINAL SPACE ***')\n",
    "print(pca.inverse_transform(centroids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAFzCAYAAAD/m0kvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8FPW5P/DPA0kIEUIg3LEQENFotd6KnnoDD9Zq1bZU\na4sexVq19lcstrVeqqKtvVjrUWuPVvFCrcXWemhVjtZ6AcXWSgEvaEHkEkHkEgIJJCHk9v398Xxn\nMzs7uzub7GZndz/v14tX2Nm5fGd2dubZ53sZMcaAiIiIKKz6ZLsARERERIkwWCEiIqJQY7BCRERE\nocZghYiIiEKNwQoRERGFGoMVIiIiCjUGK0RERBRqDFaIiIgo1BisEBERUahlNFgRkRoReTiT28gU\nEZkpIp0iMjbbZeltInKy3ffp2S6LHxGZJyJ70rzOrJ+rIjLOHvcLs1mOQuQ650/KdlmyQUQmi8g+\nEflEtsuSbWG4FiQiIjeLSGcPll8sIi+ns0wpbPtyEflQRIpTXbZbwYqITBCR+0VknYjsFZEGEXlN\nRK4UkVLXrBkby19E+ovInAxeXAwyWP4wEJGvich34rwd5n3PxGcT5v2l3tGtc0BETheROekuTHeJ\nyHUi8oUUF7sVwO+NMZtc61lsAzjnX52ILBWRi0VEfLY7RUQWiMgWG/hsE5GnReRLccp5sF1vs4iU\np1jeTAr7taCn17/unufpuOfOA1AC4PJUF0w5WBGRzwNYCeAcAE8D+DaAawF8COAXAO5KdZ3dVAZg\nDoApGVr/owD6G2M2Zmj9YTADQLxgJeZiRES+zgBwU7YL4XI9gMDBiogcAWAagN943jIANgE4H8AF\nAH4EoC+AhwD8xLOOWwC8DOAQu57LofeD/QA8KSJf9dn0BQC22P+fE7S8lDU9vucaY/YB+C2A76a6\nbFEqM4tIFYDHAWwAcIoxZrvr7ftE5EYAn0+1EN2UkZupiJQZY5qNPuGxNRPbIMoFznch2+VwiEg/\nAK0mfE9fzfXA/mIAHxpjlvq812CMedx5ISIPAHgfwLdF5EZjTIeInAPgRgBPADjfGNPhWv4OETkV\ngF/afwaA+QDGQwOible92ExPib0ZUmak6zx/AsAPRGSKMWZx4KWMMYH/AbgPQAeAYwPOvwHAw67X\nNwPo9JlvJoBOAGNd044B8DyAWgDNANYDeMi+N87O32H/Ov9uci1/EIAnAdQB2AvgXwDO8mz3Irvc\nSQDuBbANQF2CMtVAs0nHA3jDrncdgP/y2afDAbxiy74JwA+hF4WodSY4dqcAWAKgEcAuAH8BcLBn\nnpvt+g6Aptd2AaiHfulLk6x/kc8xXG/fO9lOP8eWe5Pd1xcBHOCzrmMB/NVuuwnAYgCfCbCPtQB+\n6Xotdh1tAMpd06+x08rs60cA7AYw2h6XPQC2A7gdgHi2UQbgDgAbAbQAWA3ge8nOVTttEDRT6Cz7\nAYAf+GzjqwCW2TI1AHgHwJWu9wcD+KWdvsfO8yyAwz3rcc7rCz3Tk57LCY7xEAC/s9vcZY/d4d7t\n2PNnD4AJtmy7ASxwvX+u3cdm+7n9DsBon+0dBL0YbbfzrgZwq2ee0fYc3WqP67sALvbMc7It43nQ\nKoqPALQDONJO/47Ptj/jLJPkmIyx500j9Dv/3wA+a5c9yTXfCXZfPrTl3GjnLXXN8whiv0cdrve/\nD+DvAHbY47EMwJd9ynQq9Pu+y34OqwH8xDNPCYBb7HnolOc26E3amcfvuvhwkuNRA3tt9blGvOMz\n/Qm7jZH29Sp7TuwX5Jx0HdsOAEfbc6vN73xKsHwngF9BA553AewDcLbrOjLbTt9rz7PfAKjwWc8N\n0OtbE4CXoJmhbt+37PTTodd+53qwFMDXPPMEumba4/Qvux8fALgMmt3oCHicLgOw1p57/7TrWwTg\nZdc8xdCs2TJbnkYArwKY4nNt8r3nAjgM+l1YZ8u6BZqBGxKnXDsA3Bn08zbGpJZZAXAm9Ib2RorL\nOeLVtUVNF5Fh0EBlO4CfQQ9gFQCnwWctgG9CT8AF9h+gNwOIyKEAXoNe4H4GPRm+AuAvIjLdGPOU\nZ/v32m3dAk1bxiurAXAggD9BP4h5AL4O4BERWWaMWWW3Pxp6QnRA06XNAL4BzdQk/VUoItOgN4x1\n0BOzP4ArAbwmIkeZrqopZ11PQIO5awEcZbe1DcB1CTZzK/RmPAb6xRboSRophl2+AxoEDIIGDY8B\n+A9XWU+xZV2GruDpYgAvi8gJxphlCcrwd2ig6DgcQLnd5vEAnrPTTwCwwkT/yi+CniP/BPA9aBr7\nu9Av5v2u+Z6B3vgeBPA2gNMA3C4io40x34tXMBHpD/3CjoKeZ5ugN8OfARhptwX7q3E+gBeggQwA\nVNt5f2VfTwBwNvS82QBgBDRNvlhEDjHGbE1QjlTPZfeyAmAhNPC/F/qL+AvQNKzfue0c0yXQY9ps\n1zMTGly8AT3HRkDPmc+IyJHGmN12vsPtsvugn8GH0ED6TOhNASIy3K6nwx6fHdCL+0MiMtAY4xwz\nx412fbdDb9aroefN+QDu9sx7PvQGkeiYlEKrK/a3y28B8F/QHwfeY3Iu9Lt3LzRQnAxgFvQ7c56d\n5zfQ4Gua3b731+eVtjyP2fJ/FcATInKmMeY5W6ZDoOfpW679nQg9h5xyi53nM9Bjuxp6g7gKek1y\nro0XQK9NbwB4wE5bl+B4jAYwFsCKePP4OAD6+dWLyERogPqgMaYphXWcD2CdMWa5iLwHvcF9DfrD\nIqj/hH4Xfg09j2rs9AcAXAg9Z++GZm5mAThCRI43NvMjIj+G/hhbCL3WHAXgb4jNAgW6b9l1zoQe\n/3cB/BR67zoSet153M4T6JopIp9E133wJluum+3rpETkEuj5+RqAO6HXoacB7IQGuo5y6H3sceix\nGwjgEgB/FZHJxph3kOSeCw22x6PrR8ih0GvcIXDdL1xWQK/xwaUQyQ6EHtQFKSzjjVB9I0JohqMD\nNkKFXlA7AByZYN2V8GRTXO+9COBNAEWe6a8BWO3Zbic0qvX+Wo4qk2t/OuCKgAEMhX7RfuGa9ivo\nr8DDXNMqoF+oqHXG2bc3oRfRQa5ph9l1PuI5np0AHvAs/78Atgf4fJ6BzaZ4pju/at8F0Nc1fZYt\n/yGuae8D+D/P8v2gF8i/Jtn+96AB3H729behQdfrAH5qpwn0y+XOwDxiy3G9Z33LASx1vf6C3Y9r\nPfM9YY/l+ATn6g3QG98Ez7I/tWUeY1/fCWBXkv0s9pk21p43P3RNi8msBD2X42x3ul3ft32+Hx2e\n7TjH1JsFKYJefN5C9C/4M+y657imvQK9OI9JUKYHoYFXhWf6fPs59/Ocgx+4t2vfu9SWdZKnnNvh\nkyHwLPsdu+x017RSAGvsdHdmpZ/P8tfYc2d/17R7EOeXrncd0DYf7wB4wadMgxOU+wJo9uE/PNMv\ns8se55q2B0myKa55T7HH+Qyf9xYBeA96ra2EBiV32/n/bOc5y76+Msj2XJ9VLYBbXNMeg/4gCbqO\nTns8DvJMPwE+2TXozbQTwFft66HQ7NRTnvluhScbheD3rXJoJuXv3nPWs1ygayaAP0N/nIxxTTvI\n7nfCzAq6vrfL4Lp2QIOQTkRnVgSx15dy6D1ormtaonuu33flPHt8jvd57zcAGoN+3saYlBrYOq21\n09plNI566AE8W0RSbVczGMBU6K/YQSJS6fyDRs0Hisgo1yIG+oEkzXhY/zbG/COysDE7oCffBNc8\npwF43Riz0jVfPYDfByj/SACfggYlDa7lV0J/vZ/hWcQgOpMA6K/bShEZEGiP4nvYRNc/L4F+LhNs\nWY+A/qp73HOcB0JTqslajS+BfqmcX5An2mlL7P8BDdIq7DQvv/12fw6nQ28s93jmuwPauPz0BGU7\nx66vwbNvL9kyO/tWD2A/ETkt3oqMMW3O/0Wkj4gMgWYt3of+mvPVjXPZ6zRoYPWgZ/r/IH79s7eR\n5TEAhgO41xgTacNljHkW+uv+87asQ6Gf2UPGmM0JyjQdGiT39dmfQYg9HvPc27WegGYfzndN+xz0\nYvpYgm0D+plvMcY4vwxhjGlBVxYCrumR9g8iUmbL+Tr03DkyyXb81lEBrRJcguj9rLd/v+TXy8Y6\nB1rdssZz3BZBP8upQcrjoxJ6DdkV5/1qaGBRa7f//6Cf3yX2/e7cF86AVk8+7pr2OIBPiUh1CutZ\nbIx53zPtHOjxfMlznN6EZo6d4+S0o/FeG3rSQeRUAAMA/NznnAUQ/JopIn2gVZN/dn+f7P4+H6As\nzvf2N8aYdtf030IDqgij2u12xV53SqCBTtzrk2cd7vO8n92nN6Dnpt86dgHoL9G9hxNKJRDYbf8O\nTGGZbjHGvCIiT0JTX1eJyGJoHfP8eCeBy0ToAfoxNEqOWT30Q9zimlaTQvH8egftgl6EHOMA/MNn\nvrUB1j/O/l3j894qAJ8Vkf7GmL0JyuRceAYjumonVZs8r93rBfRLB2jPKT+dIjLIHXR5rIDetE+E\nBmInQj/zbQBmiUiJnWagmQS3FmNMnU/5vJ/DxyY2Pb3K9X48B0IDpVqf95xzCNAqgnMBPCsiH0Nv\nuk8YYyIXFHsDmg3gCmiqtK9rPTsSlKE757LbOOiNucUzPd552G6M+chnHQb+5+NqdKVynSDxvTjr\ndqp3K6DZAL+ui+7j6qiJmcmYBhF5BtpeYY6dfD6AzcaYRfG2b42D//57b3oQHXPkx9Dsgfu8MtDA\nKikRORNa1XAE9Nezwz1Oxh+hN/+5AH4uIi9B0+xPun5EHQjgYCQ/H7srXpC0AVqtDNh2W/YHmqM7\n94UL7HrbROQAO209NNN4PrqqDMuh1XCOVmOMO6iq8Vn3gdBzzK+qxH2cnPGzos4FY8wOEYkXuCXj\n7Evc7wACXjOh2b7+3vJZ7yPxDy2g63vr3b92EVnvnVlELoJWbR+M6GqwmHn92ADnZmg2xX0uxvuu\niOv9QAIHK8aYPfZi/Mmgy/itJs70vt4JxpiviMhk6IXiNGhd2HdF5DiTuIeCky36JeJHoN4TYK/v\nXP464kzPZo+ATJUp2XqdY/09aHsQP3GDJfvFeQPASfaiNRLaTqQW+oU5FprWXe0TmMQrW7r0gQZQ\nt8H/OK4BAGNMrf21dBr0AnI6gItF5LfGmIvtvD+ENmB7EHoh3gm9Wd2NxMMHdOdc7olM96Rw9ucx\n6C88P+94Xsf7bj4K4BwROQ5aXXkWtO1CWthfti9Cb3w/g94gmqDtVX6LAMM+iMiJ0PYqi6GB6hZo\nCv/r0PYZACKZnZNEZCo0U/U56EX/JRH5rA1Y+kCHjLgK/uej94dFUHV2fYPjvN+UJABcbf8eFmRj\nIjIQ2oapH7SKz81AA9Ab7Ou7oVUtjsXQaiuH37nRB/pjZwb8j5NfsJdM4PtWAEGvmYEzDj0lIhdA\nq4EXQLubb4etZkd0pjqRPwE4zi7/NnQf+kCvW37flcEAmk0KvbdSbWC7EMClInKs6V4j212ARszG\nNsqzqvxmNtqVbimAG0Xka9BqlK9CA5d4J5ATCbYZY7IySh+0YeFEn+kH+kzzWxbQukmvgwHs8GRV\neiJwVBuH03BvTw+O9RJow9RpAGqNMWsAwDa6OwmaWXmmm+v+EMB/ish+nuxKtev9eNYBGBDglzps\nCvX/7D+IyH0ALhORHxtj1gP4MrSO+DL3crZaINHFs6fn8ocApohIqSe7EuQ8dK9DoOfjYs97B6Hr\nGDplTfRjphZaXdA3Dd/Nv0KzUudDrxH9kbwKCNDyHuoz/WDP68Ogx+m/jDGR6lvRxu9e8b5H06E3\n1NPcqXjb8DF2JXquLQLwfRG5DppNmwptELwO2nss6fmYoDx+nGBjfArLdG3ImA9E5H0AXxCR7yT5\nIQnod6EftLGm9wfIQQBuFZHP2Kr226C9zhxBMh7roA1v/5HkRuictwfClaGx1ZnewC3ofWsd9Lvy\nScTPSAS6ZopILfTc8fuues9VP8739kC4vre2WcV4aBs0x5ehjZ2jxroRkR951ul7Xtnr2CkAbjTG\n/MQ13e8e6BiPrgx3IKkOCvcLaNr+QduqP4qIHCAiVyZY3vkwI20ZRGQ/aMtt93oqfJZ1olAnlep8\nKaLmNcbUQj+cy237D28ZhyYoX7o8D+A/RHtHONsdAo32EzLaM+QtABeJa1RH2zL8s7A3xDRpQsB0\ndhzLoZ/p9+3nGCXgsV4C/RUxG9FVPa9Be2mMgn97lSCehQbk3/ZMvwqa2XguZokuT0A/w8963xCR\nQSLS1/5/iM+yTlsl51ztgOdXnoicC/2VHlcazuXnoXXPl7qWEWi7g6A3tGXQX1rfFNcQ2SJyOjTo\nW2jLugOaFfu6xBmy3RjTCW38/WXRXk6p7o97XR3Qdg7nQbuQrjTGvBtg0WcBjBaRL7u2WwbXMbKc\nzJ33Gjkbsceuya7HOwprB7p6WTnbqoJnwDabQvd6G3rOOOfQEwD2FxFvOSEipXYf3OXxu4bGMMZ8\nDM3KHBNk/jjmQBusPuR8LzzlO1VEnLZ2F0Ab9c81xixw/4O2JWuCbYtkjFltjHnZ9e/NAGV5Anq8\nYwbpE5G+tooF0KxZO7TTgNtVPusMdN+CVgHvAXCd6JhAfgJdM+135XkAXxSR/V3vV0PvA8ksg+3B\nI9HtPi9G7LkRk6UWkWMR24vH956L+N+VqxD/OnMU/JtKxJVSZsUYs15EZgD4A4BVIvIoNAVbAm0k\neS40nRTP36DtKx4WkdvR1WVrOwD3Be4iEfkWtDX0Omh96KXoGp8CxpgWEfk3gPNE5ANoav1dY8x7\n0IvxEgArRWQuNModAT34YxDdOC4T1Te/gH4pXxSRe6BfwG9Ao93BSH6juBq6n/8UkYegY4V8Gxrh\n35LGci4H8BURuQPal7/RGLMw6MLGGCMi37BlfU9EHgGwGXqMp0I/r2Qjab4OvWhMQnSD2VehqXOD\n7gcrz0B/qf5ERMajq+vyWdA+/hsSLHs7tLvxQhGZBz1W+0G7V0+H/qraCQ3ch0B//X5kp38bwJvG\ndmWH3tBvFH3eyD+gv9rPR4IupS6pnMtef4FmHe4QkQOhv6LPRtfFJmnAYqvqroFmM18Vkceh1XVX\n2rK4GyReacu6QnTwsA3QX1BnGGOccl4LHQHzDbs//4Y2tjwa+usslR8Tj9ptTkFXt/Fk5kI/n9+J\nyDHo6rrsbde0Gvr53GFvFruhv0D9goDl0OvIPSLyPLSnxh+hPyy+C+B5EZkP/dy+Ba3+ONy1/E2i\nQ5j/H/QaMQJ67m9EVwD/O2g33ftsddHfodUQ1dDr7mfR1f14OYBpInIVgI8BbDD+A745ngLwxQTv\nJ2SMeUJEDoNWGxxpz5EPoY13Pwf9XGeINgafgjiNWI0xrfb4nSsiV3oa9wcty6sicj+Aa2317N+g\nVW+ToI1vr4T2aN0hIr+08y2EXsOOtOX1ZjsD3bdsU4mroOfYv+xnvgvaYaK/MebiFK+Zc2x5XhOR\ne6FV49+G3nPd54/fcWgXkRugDeYXicgfod/FixF73VkIYLqI/AV6Dk6Atil7D9pg2Fln3HuuiLwK\nHeitxO7PZ6HXQr/HMhwN/c7/JdE++O1Uyv+gDYl+g64BYOrRdXMpds0XGcjNNe0I6AV7L/RidiVi\nu4AdAU3pboBGc1vsjh3pWdex0IvxXru8e1C4KmjgtBldAyg9BeBLrnmc7R7ls49RZXLtz1M+8y4C\n8JJn2uHQX8XN0C/u1dATrQPAsADHeKo9ps6gcH9GbDe9OXZ9Q5KVPc42yqAXwTo7v3dQuOme+cfB\n0+XVta9/QtdAYOuhv3qnBDyf3oAGLMe4po2229rgM/8j0JE1vdPnQBuJevfxl9Bfj86gcFf5LOt3\nrpZBU/Hv23NsG/RmPBu2SzeAL0EzNFtc5/T/ABjuWk8JNID9yH6er0DH7HjZfd4kOL7xzuUvBji2\nzqBw9bDBFTTQ6QRwbrJj6nr/HEQPCvdbAKN85qtG1wB2TdBgZI5nnqHQ7v01dn82Q28IX3fN43sO\n+mxvJfRmFFOWBMvsD/0+7bGf6R3QnhzerssHQX/dNtj57oOm+L3dvvtAb8Bb7XnsHhRupj3nmqEX\n/wvh6QoLvYEvQNfgi5vsZ3aAp9x9oYPMvWPXtwN6/fshtMrSmW8S9JrUaMuabFC4I+AZksFOXwTg\n7RSOq7MfW6Dtn7ZCr9uft+9fZbcT97pgj08HgDOTbKsDwN0J3r/EHptGe+6/BR12YIRnvhvQ9b18\n0Z6/3bpvueb9PKIH9HwdwFc88wS6ZkLb7Dn3uA+gP9qjzp8kx+lydA0K9wa0QXzUdcfOd40tgzNw\n4emwg7x55vO950Kv1873fqfdlxF2nhs96/g5fK7ryf6JXZh6gYjcBT3ZBhgeeMoSEfkitDrmBGPM\n69kuT0+IyAroqNOnZrssuUxEXoT2nPNWbRCljc281EDH0UqpQXy3nrpMyYmn/7hov/MLACxhoEK9\nxec87AOtp9+N1EYtDR1bjXME4vcsouCuh1YJ+7Y3IkqTi6FjP3nHyEqKmZUMEZE3odVAq6B1/F+H\nNhY9xRjz9ywWjQqIbRfSH5qK7gdtd3EcgOuMMb/IZtm6yzbOPQbaHmQItLqEDx0lymOpdl2m4P4P\nWs9/KbQh43Low9oYqFBvehl6U/88tNfVWujw+/dltVQ94zzldzX0AXEMVIjyHDMrREREFGpss0JE\nREShxmCFiIiIQo3BCuUkEZknIokGdXPmGycinSLCLplpIiI1IvJ0L23n4Qyte6Y9L8YGmPciO2+g\nJ9Cmi93mrwLMF3hfiHIVg5U8JyITROR+EVknIntFpEFEXhORK73dWtO83VEiMkdcjxxIM4PoJ9dS\n7+mthm6Z3I7xrl9ErhB9+mxvl6WnUt0XopzDYCWPicjnoSN8ngPgaegIutdCR9T9BeIMe50mo6Ej\nLR6RofV/A8Ee6EXk51HoEOgbXdO+hein/OaKfNoXIl/supyn7APTHocODX2KMWa76+37RORGaHfW\njBUhpZlF+psUniZt9LkhKT87JFtEpMwkfyJtqKX6GYWZHZgxL7o858K+5NO5Q9nBzEr+ugb64L1L\nPIEKAH0opTHmHvc0EblARJaJSLOI1InI4+4nftp5FovIOyJSLSKLRKRJRD4Skatd85wMfX6EATDP\n1qd3OO1GXOs4SkReFZEmAO5Hi39LRN4VkRYR2Swiv5aup6U688S0WRF9GvI8EakXkV32IWExD54T\nkREi8oiIbLLb+FhE/pKszl9EDrPLOVVqW0TkIfE8eVlEbrb7XC0i80VkJ1wPYxSRg0TkSXuM94rI\nv0TkrETbtss57W++a4/ROnv8nxeRMXaeG+1+Ndt9qvCs42wRWWiPa4uIrBWRG+zItu75En5GPmW7\nSETaROQ21zQRkdn2s9wrIltF5DfeMtl5b7DlbhKRl0TkkGTHwy63XESe9ExbaY/TJ13TzrPTDrKv\no9p52HPpUABT7PROEXnZs7l+IvLfIrJdRBpFZIHoyNSJyneWT1mm22necq8SfQigdx1fsPvUYo/l\naZ73U9oX+z25S0Q22nV+ICI/EJFAPzBE5HQReUVEdotWKy8Vka+53k/H93uiiPyv/Y7ttefG4yIy\n0DXPqSKyRPS7vkdEVotI3HOUchszK/nrTOiDCd8IMrOI/BDAj6BP1J4LYBj0YV2viMiRxpjddlYD\nHTX0OehDy/4ArWb6uYi8Y4x5Hjpq7012ffej60b9D9c6hkKfPPoHaBp7my3HzXbZvwG4F/owuW8B\nOEZEjjddT2KNqaeHVnV9BvrAudXQhwz+1me+BdAHlv0KWiU2HPogu7HQhwTGcyr0yaUPQx/Sdij0\nQWGHIPpx6s72/gRgDYDrYDNNoqOvvgZ9eNrPoA/7+wqAv4jIdGPMUwm277gA+gTWX0E/i2sA/Mne\nkE6GPihsIvTz+yW0yswxE/oAvzugD1o7Bfo5DbTrce+D72fkJSKXQY/5rcaYOa63HoA+mO5hAHdD\nj90sAEe4P0sR+TH0YXwLoefVUdDPvzjAsVgC4KuusgyGfh4dAE6EPqEW0AfCbTfGvO/aP/d58R0A\nv4Yem1uhn5d7f8W+vxPAzdCHS15lp30N8b1mt3OSqywnQttbneAq91DouX63Z/kToU/5vteW7UoA\nT4rIWGPMrlT3RUT6Qx+QOgr6MNpN0O/Mz6AjbX83wb5ARGYCeMjuy0+hDwk8Evo0cyfQ6tH3W0SK\n0fX5/wr6XRsDvaZVANhjg9lnoA8ovBH64MSJdl8oH6X65EP+C/8/6I2nE/oo9CDzj4U+ufYaz/RD\noOnla13TFkFvBDNc04qhj6J/wjXtaFuGC32256zjG57pQ6FP4X3WM/1bdv6LXNMegX1KtH39Bbu9\n77qmCfQJx5Gn5AIY5J0vhePaz2faeXb9x7umzbHb+J3P/C8CeBNAkWf6awBWJ9n+OLverYh+yu5P\n7PQVAPq4pv8e+nTU4iT7cB/0xuaez/czsu9tAPC0/f+Vdr7rPPOcYMt0nmf6qXb6Vz2f+VOe+W61\n8yV7YvCX7fYPsq/PtPv8ZwDzXfO9BeBJ1+uYJ+ZC23e97LONi2xZ/uqZfgf0+zEwSRlXAnjc9XoZ\n9CbeAWCSnfYl+/qTrvk67b5UuaYdZqd/q5v7cgP0uVATPNN/avdlTIL9KIc+gfrvAEoSzNej7zeA\nT9l9/FKCbXzHLjM4yHeX/3L/H6uB8lO5/bsn4Pxfht7Y/yQilc4/6OPLPwAw1TN/ozFmvvPCGNMG\nrfaZkEIZ9wGY55k2DRr4eBv+zoXuS6I2NqdDA67fuMplANyD6PYze6EX5Sl+1RGJGGP2Of8XkX72\nGL1h1+/t1mrgeViX/dU/FZpxGeQ51n8DcKCIjApQlCeMMY2u10727HfGmE7P9BLor1K/fRhgt/0a\ngDLENlj2+4zc+3M19LO62hjzM8/b50B/db/k2c83oRkd55w6FfqZ3+NZPmjj7yXQ43+SfX0i9Fx8\nwf4ftorhk3BVxXWDgWaKvNvuCw0ik5XRKctA6M34AQB1znT7t94Y865n2ReMMTWRQhizEjbYSH0X\nAOjnsgQOvtywAAAgAElEQVRAg+dzeQmaaT8pwbKnAhgA4Ocm+SMOevL9brB/P2czQX7q7d8vBa2+\notzGYCU/OVU2AxPO1WUi9FxYC6DW9W879AY23DP/Rz7r2AVgcApl3GyMafdMcy76a9wTbTC0Holv\nCuMAbDGxjVjfd7+wF9lroMHNNlv3frWIjEhWYBEZLCJ3i8hWaNBTa8tloBkbL+84MBOhN9YfI/o4\n10KrFoDYY+1nk+e1c3H3fi7O9MjnIiKHiMifRaQeep7UAvidfdu7D36fkWMKgNugN67/9nn/QGjK\nfjtiz6n90LWfTjuhte6FjTE7oOdUQkbbY32A6Jv+EvtvjGhD8xOgx70nwQoQe9yd8iU775cAGCUi\nE6DVFJ3QB0tGghhbRr/nhnm36Ww3le+a24EAPofY8+8F6Hmc6Pw7wP59L8B2uv39tsHZHdDqyx0i\n8lfbzqXctdgfocdrLvR7/LiInMvAJX+xzUoeMsbsEZGPob8mg+gDvYB+Dv5jlzR6XsfrhZPKhSJr\nPQOMMXeLDmr2RWhd+48AXCciU40xbydY9E/QJxb/AsDb0OPSB8Dz8A/8vfvozPNLu4yftXGmu8U7\n/gk/F5theBX6q/QG6A2iBVpl93PE7kOiz+hdaDDyXyLygPvXv9UH2k5hBvzPi9oE607VawBOER03\n6Gho4PcudD9PhFZnNkKzOj3R3fP+NXRlfw4AsMIYs1dElgCYJSL7Qdt9XJ/GbcbTBxqY3BZnHWt8\npnVHj77fxpirRWQetHr3s9C2K9eKyHHGmI+NMS0AThKRqdCMzOegVbIvichnbVaV8giDlfy1EMCl\nInKsSd7Idh30wlVjjAlyswyiOxeLD+3fgwDUOBNtg7vx0ItsomVPkdguwr5jsRhjNgC4E8CdInIA\nNPj4HrRBaAxbZXQKgBuNMe6eDROT7JPbevu3zRjj7WnSG6ZAf5F/wbie/m33P1U7oFUKf4feII43\nxmx1vb8OwH8C+Ie76smH85kfiOjPfCiCZw+WQBsOfxV6M37dGGNE5DVogFBty5HsnMzIDc4Ys0lE\nNtqyTEBXhudVaAbhXFvuV9O52TjT10HbOy3qxjqd68Qn0XUupyKl77cx5j1oFuenInIctIH+N6EN\ndJ15FkHbyHxfRK6DtnWaCn3aOOURVgPlr18AaAbwoIjEpHZF5AARudK+XADNqMzxzmfnHeI3PYkm\n+zeVdiEvQtudXOmZ/g1oO5yFCZZ9FloffoUzQbQ77iy4Ltwi0l9E+nmW3QCtM/dOd3N+4Xq/M1ch\n4E3OGFMLYDGAy0VkpPd9e4POpA7ozSayDyJSAm3gmDJjzMfQdgj9Abxg2+Q4noD+GLrJu5yI9HV1\nVX0RQDv0c3K7KoWiOO1WrgHwjjFmj2v6f0KzLUGqgJqQ2vmaiiXQYPfTrrK8Bc34XAvNRCxP4/bi\n7csTAP5DRD7rfUO0S3PfBOv8G/R7cp3PdyiIQN9vERnoU473oNeofnYev0D2beh50J2yUcgxs5Kn\njDHrRWQGtNfBKhF5FJoaL4HWm58L7VHjzHsD9BfMeAB/gV6UJkCrSu4H4NcuIZF10DT8N0WkEXrx\n/Kcx5sN4CxhjdojIzwDcJCJ/hXZFPhgagCyF9m6J5xnor/yf2334N7TLp7fdziRoJuAJO0+7nW84\nurpe+pVtj4i8CuAH9ga/GZqerkJqKfn/B71ZrRSRudBfqCOgXZ/HQKsD0sldtn9A2zs8Kl3PnLkA\nPcgoGGPW2RvfYgB/E5FTjDF7jDGvisj90NT9EdAbXRv0+J8DvWEtsJ/5L+18C6FB55HoalcRtAxb\n7brdDXVfhVZ3GAQLVpZDz9cfQqvjtrsyEPE+46Cf/RIA50NvuK/ZcneKyD+gVZGLErQP6o54+3I7\ngLMBLLTVLMuhbYgOh34PqqDds2PY78BV0HYi/xKR+dDz6VPQEXQvTlSgFL7fpwD4tYg4Xf+LoBnP\ndgDO2DQ3ichJAP4PmrEZYdezEfb4Up7Jdnck/svsP2gd+W+gwcNeaADxKvSLXeyZ94vQrr677b/3\noOM+THTNswjA2z7beQTAOs+0M6FdKPchuvuw7zpcy11ht90C7RJ9D4DyANurgPZA2AW94D4CvQi7\ntz0EWv/9nt3HndCb+PQAx3IU9GJZZ5d7HHqR7IBWDznzzbHThsRZT5Ut22a7jxsBPAXgi0m2P86u\n9yrP9JPt9Ome6U6X1qNc046DBnWN0MabP4VmRzoAnJTsc7bvrUdsV+NP23NrEVzdowFcAr0RNdr3\n37LbHOFZ/gZoA+FG6C/warudhwKe53+0+3COa1qRXV8zPF1t4d/ddzj0Blpv33s53nH0HPeTApSv\n2s670jP9ejt9js8yHQDujnP8H+rOvtj3yqDVJe9DrwnboMHUbAB9A+zL5+38jdDv2usAvhLk3Any\n/YZ+P+ZCA5UmaND6IoAprnmmQDPCm+w+bII2FD8gyPnCf7n3T+wHT0RERBRKbLNCREREocZghYiI\niEKNwQoRERGFGoMVIiIiCrWc67psn2NxGnRQoZbsloaIiCinlEJ7XD1vjKnLclkCy7lgBRqoJBpv\ng4iIiBI7H8D8pHOFRC4GKzUA8MDnDsOkIQOyXBSKZ0fzPsz/98eYcchoDC3jgJJERGGwZmcjLvvr\nSsD1yINckIvBSgsATBoyAEcML082L2XRtKph2S4CERH5y6lmFGxgS0RERKHGYIWIiIhCjcEKERER\nhRqDFSIiIgo1BitEREQUagxWiIiIKNQYrBAREVGoMVghIiKiUGOwQkRERKHGYIWIiIhCjcEKERER\nhRqDFSIiIgo1BitEREQUagxWiIiIKNQYrBAREVGoMVghIiKiUGOwQkRERKHGYIWIiIhCjcEKERER\nhRqDFSIiogLw7IUP4u+fvzHbxeiWomwXgIiIiDLn2QsfzHYReozBChERUR7KhyDFwWogIiKiPJNP\ngQrAzAoREVHeyLcgxcHMChERUR7I10AFYGaFiIgop+VzkOJgsEJERJSDCiFIcbAaiIiIKMcUUqAC\nMLNCRESUMwotSHEwWCEiIgq5Qg1SHKwGIiIiCrFCD1QAZlaIiIhCiUFKFwYrREREIcIgJRargYiI\niEKCgYo/BitEREQhwEAlPlYDERERZQkDlGCYWSEiIsoCBirBMbNCRETUixikpI6ZFSIiol7CQKV7\nmFkhIiLKMAYpPcNghYiIKEMYpKQHgxUiIqI0Y5CSXmyzQkRElEYMVNKPmRUiIqI0YJCSOQxWiIiI\neoBBSuaxGoiIiKibGKj0DmZWiIiIUsQgpXcxs0JERJQCBiq9j5kVIiKiABikZA+DFSIiogQYpGQf\nq4GIiIjiYKASDsysEBEReTBICRcGK0RERBaDlHBisEJERAWPQUq4sc0KEREVNAYq4cfMChERFSQG\nKbmDmRUiIio4DFRyCzMrRERUMBik5CYGK0RElPcYpOQ2VgMREVFeY6CS+5hZISKivMQgJX8wWCEi\norzCICX/sBqIiIjyBgOV/MTMChER5TwGKfmNwQoREeUsBimFgdVARESUkxioFA5mVoiIKKcwSCk8\nDFaIiCgnMEgpXKwGIiKi0GOgUtiYWSEiotBikEIAgxUiIgohBinkxmogIiIKFQYq5MXMChERhQKD\nFIqHwQoREWUVgxRKhtVARESUNQxUKAgGK0RElBUMVCgoVgMREVGvYYBC3cHMChER9QoGKtRdzKwQ\nEVFGMUihnmJmhYiIMoaBCqUDMytERJR2DFIonRisEBFR2jBIoUxgsEJERD3GIIUyiW1WiIioRxio\nUKYxs0JERN3CIIV6C4MVIiJKCYMU6m2sBiIiosAYqFA2MLNCRERJMUihbGJmhYiibGlswTWLV2FL\nY0u2i0IhwUCFso2ZFSKKcteyDbj/rY0AgNumVGe5NJRNDFIoLBisEFGU2ceMj/pLhYdBCoUNgxUi\nijJqQCkzKgWMgQqFEYMVIiJikEKhxmCFiKiAMUihXMBghYioADFIoVzCrstERAWGgQrlGmZWiIgK\nBIMUylXMrBBlGAdZozBgoEK5jJkVogzjIGuUTQxSKB8wWCHKMA6yRtnAIIXyCYMVogzjIGvU2xio\nUL5hsEJElCcYpFC+YrBCRJTjGKRQvmNvICKiHMZAhQoBMytERDmIQQoVEgYrREQ5hEEKFSJWAxER\n5QgGKlSomFkhIgo5BilU6JhZyRIOwd5zPIaU75698EEGKkRgsJI1zhDsdy3bkO2i5CweQ8pnDFKI\nurAaKEs4BHvP8RhSPmKQQhSLwUqWcAj2nuMxpHzDQIXIH4MVIqIsY5BClBiDFSKiLGGQQhQMgxUi\nol7GIIUoNewNRETUixioEKWOwQoRUS9hoELUPawGIiLKIAYoRD3HzAoRUYYwUCFKD2ZWiIjSjEEK\nUXoxs0JElEYMVIjSj5kVIqI0YJBClDkMVoiIeoBBClHmMVghIuoGBilEvYdtVoiIUsRAhah3MbNC\nRBQQgxSi7GCwQkSUBIMUouxiNRARUQIMVIiyj5kVIiIfDFKIwoOZFSIiDwYqROHCzAoRkcUghSic\nGKwQUcFjkEIUbqwGIqKCxkCFKPyYWSGigsQghSh3MFghooLCIIUo97AaiIgKBgMVotzEzAoR5T0G\nKUS5jcEKEeUtBilE+YHVQESUlxioEOUPZlaIKK8wSCHKPwxWiCgvMEghyl+sBiKinMdAhSi/FXxm\nZUtjC+5atgGzjxmPUQNKs10cIkoBgxSiwlDwwcpdyzbg/rc2AgBum1Kd5dIQURAMUogKS8EHK7OP\nGR/1l4jCjYEKUeEp+GBl1IBSZlSIcgCDFKLCVfDBChGFG4MUImJvICIKLQYqRAQwWCGikGKgQkQO\nVgMRUWgwQCEiP8ysEFEoMFAhoniYWSGirGKQQkTJMLNCRFnDQIWIgmBmJc/w8QGUCxikEFEqGKzk\nGT4+gMKMQQoRdQeDlTzDxwdQGDFIIaKeYLCSZ/j4AAobBipE1FNsYEtEGcNAhYjSgZkVIkorBihE\nlG7MrBBR2jBQIaJMYGaFiHqMQQoRZRIzK0TUIwxUiCjTmFkhom5hkEJEvYXBChGlhEEKEfU2VgMR\nUWAMVIgoG5hZIaKkGKQQUTYxWCGiuBikEFEYsBqIiHwxUCGisGBmhYiiMEghorBhsBJSWxpbcNey\nDZh9zHiMGlCa7eJQAWCQQkRhxWqgkLpr2Qbc/9ZG3LVsQ7aLQgWAgQoRhRkzKyE1+5jxUX+JMoFB\nChHlAgYrvSiVqp1RA0px25TqXioZFRoGKUSUS1gN1ItYtUNhwECFiHINMyu9iFU7lE0MUogoVzFY\n6UWs2qFsYJBCRLmO1UBEeYyBChHlA2ZW0oBjolDYMEghonzCYCUNnIazAFjNQ1nFIIWI8hGDFR+p\nZkrYcJbCgIEKEeUrBis+Us2UsOEsZRODFCLKdwxWfDBTQrmAQQoRFQoGKz6YKaGwY6BCRIWEwQpR\nDmGQQkSFiMEKUQ5gkEJEhYyDwhGFHAMVIip0zKwQhRSDFCIixWCFKGQYpBARRWM1EFGIMFAhIorF\nYKWHtjS24JrFq7ClsSXbRaEcx0CFiMgfq4F6yD3a7d133onOhh1oXbwg5fWUTJmOPoOGouWpB9Jd\nRAoxBihERMkxWOkhZ5Tbn/7idvSbMj0yPZWApWTKdPT/wmWR1wxYCgMDFSKiYBis9NCoAaW4+847\nowIVJ/AIErB4AxVnPQxY8heDFCKi1LDNShp0NuyImdb/C5ehxBXA+PEGKonWR/mhNwOVXbVb8dvb\n52BX7dZe2yYRUSYws5IGTgbFG3gkyrDEC1T2PvVAt9q8ULhlI5vy9Lz78PzjDwMALrr6ll7fPhFR\nujBYSZNUAhYGKoUjm1U+Z8+8IuovEVGuYjVQGrUuXoC9Pm1N3FVCmQpUCrELdZj3+dkLH8x625TB\nw0bioqtvweBhI7NaDiKinmJmJc0SZVj6nTwdfSqGxiyTjoyKuwv1bVOqe7SuXBHWfc52kEJElG8Y\nrGRAvIAlU4EK0NWF2vlbCMK4zwxUiIjSj8FKhsQLWNzS2UZl1IDSUGUXekOY9plBChFR5rDNSjcE\nbSvRungBOuv9uyF31ndvpFsKlzC0TSEiynfMrHSD01aiqa0d+xUXYfYx4zFqQGnMfCVT/NuoAFol\nVDJles4GLFsaW3DXsg1x9z3fMUAhIuo9DFa6wWkj0dzWEbeBZ7xeP26pjHQbNmFt3NobGKgQEfUu\nBivd4LSV2NLYgrLivjENPOOOTFu/IybTkqsBSxgbt2YagxQiouwo2GAlHdUYfg08k42j4vd+LgYs\nYWrcmmkMUoiIsqtgG9g61Rh3LduQtnUGGfAtyMBxFB4MVIiIsq9gMyvprsZIZWTaoEPzF3oj1mxi\nkEJEFB4Fm1lxV2MEHbI9Xpfl7gyhHyTDkonsDyXGrshEROFTsJkVRyq9WuLN22dQ90amjZdh2VtW\ngesXr8KM6tEACqsRazYxSCEiCqeCD1ZSqQ6KN++G3/8KH320E9Mu+AaA1Eam9QYs+xYvwPU/uLpg\nuwVnA4MUIqJwK/hgJZVeLfHmvWvZBtz/4KV4AcDxAztT7tXTungB9uxrxzutJdj/lccKsltwNjBI\nISLKDQUfrKSDE1Ts/8pjaO1mQ9ibfnYb7n9rIy4/Yixum1LNjEqGMVAhIsodBRWsZKp3TTrGHGE2\npXcwSCEiyj0FFayEeYj4QhpkLRsYpBAR5a6CClaYvShMDFSIiHJbQQUrzF4UFgYpRET5oaCClaA4\ncmxuY5BCRJRf8nIEW/dIs/H+nwhHjs1dDFSIiPJPXmZW3A1pAUT+39TWjsfe+xhNbe349amHxV0+\naNsWZmDCg0EKEVH+ystgxS/YmH3MePz09bUAAIEkXD5o25Yw9y4qFAxSiIjyX14GK95gw/n/9f8x\nEWXFfRNmTFLJlrB3UXYxUCEiKgx5GazEEyRjkkq2hL2LsoNBChFRYSmoYCUIZkvCi0EKEVFhysve\nQD3hZEv8qoD8ehMF7WFEPcNAhYiocDGzkgK/KiI2ss0sBilERFSQwUp3uxzH62XknUbpwUCFiIiA\nAg1WupsN8WtQy0a26ccghYiI3AoiWPFmUpgNCScGKURE5KcgghV3JmX2MeO7PeosR6zNDAYpRESU\nSEH0Bpp9zHhcfsTYSKCSynN/3L19gizL3kGpYaBCRETJ5FVmxZ35ABBT9XPXsg2YUT0aQPAqIG9W\nJtmyYegdlCsZIAYqREQURF4FK/EeYHjblOpuBxHuAMVpTOtkT/yCAff8PQkaerJsGAKmRBikEBFR\nKvIqWAnStTjVRrV+vX0SBQPu+a9ZvKrbQUNPAo6wNiBmkEJERN2RV8FKvAcY+r3XE0GDgZ4EDT1Z\nNmzdqRmkEBFRT+RtA9tMNnRNNCR/d+ZL97JhwkCFiIh6Kq8yK25hb7eR7xikEBFRuuRssLKibgsa\nzK647x87rhibmytw7LhiAOHvIRP28gXFIIWIiNItZ4MVt7rmNixYVYfp1ZWoLNPgpLKsGJcePRIA\n8Mq2DzF3+VYsXFOPzc27I9ODOHnEuIyU2cubCcq14IVBChERZUpeBCsLVtVh4Zp6AIgbiEyvroz6\nG9Qr2z7sVpm8QU6y4MPboDaXqrEYqBARUSblRbAyvboSLe2daGnvRF1zWyS74ubOtLj5ZWXS4ZVt\nH0at2wmoEmV2zqguw5qmbRg1YFxoux97MVAhIqJMy4tgpbKsGKVFfbBwTT3W7mzBTSd/ApVlxYEC\nkSBZme5yrzuVzI6TzXGClzVNwbbXW1VWAIMUIiLqPXkRrAAaBLy7vRk19a1YsKoOlx49MqPVQ0HL\n5PyNl9lJp+5WWcXjF/wwSCEiyl27d9VluwjdkjfBSmVZMWZNHoV7lm7B1KpBAPwDEW+2JZNBRG8E\nKJmU7uCHiIiy69Vn/pTtInRLXg0Kt6imATX1rVhU0wCgK1hwVwE52ZYFq3IzuiQiIuquk846N9tF\n6Ja8yawAwRraOlmWqVWDMHf51rQ3rM1nTVe/kO0iEBFRD5QPTn+Th96QV5kVJ+h4cf1uzF9ZG3ee\nS48eiUU1DTEZlrrmNsxdvhV1zW29Ut5sKqR9JSKi3JZXwUpQdc1taGnvxLQJ5VHtWXqjiigsQQKr\nw4iIKFfkXbAy47BhmDahHACwtm6vb2CwYFUdXly/G6VFuvvOPNOrK3HmpIqYnkHpDDBSDRIyFdzE\n21ciIqKwyas2K0DsmCs19a0AEBmYbXp1ZVQvISd4aGnvRGlRH0ytGhQzNks6x2JJtat0psaByfWe\nSkREVDjyLlgBohvRLqppiApKAL3pOzdqZ96W9k4sXFMfGavFmc89TzqyEKkGCYm2nanRd4mIiMIk\n76qBAA0IpldXRgIV57VftYcTPMw4bBjOnFSBWZNHReZzqmAAxHSB7s19ibdttjsBdtVuxW9vn4Nd\ntVuzXZSkcqmsRERhknfBytq6vfjOc+vxyJvbom7kiW76boP7F0Xm661goLvtUjLV7iRIecJy4316\n3n14/vGH8fS8+7JajiByqaxERGGSd9VA9yzdgpr6VnR0mpgbubvaBNDMhFNV1NLeiRfX7wbQveqf\nnlTJdLddSqbanQQpj3PjBYCLrr4l7WUI6uyZV0T9DbNcKisRUZjkXbDiDLk/a/IoTKzsH5le19yG\nH72yKdIeBUBUG5VpE8qjgptUg4+eNITN5POJuiNIeTJ9491VuxVPz7sPZ8+8AoOHxT+eg4eNzGqw\nlIpcKisRUZjkXbAysbI/7j59QqQqwwk2FqyqQ019K6oqSqJuwu5GuH7D8gPBgo+eBBxh65kTpDyZ\nvvGGJXNDRETZl3dtVhze9iZTqwahqqIEsyaPinqA4cTK/r5tWVJtDxKkTUxYBoRzpLM86W7DcvbM\nK3Da176eMHMTlnYzRESUWXkbrEyvrsS0CeWR5wQ5Dzm8Z+mWrAULYeu9k87ypLvxqJO5SVQFxAar\nRESFIa+qgbztTJzB4UqL+mB6dWWkfcqCVXVJqzm81UDedXenQW0utk0JKhuNR9lglYioMOR8ZsVd\nleHNFLircirLijFr8ihUVZRgatWgpFUg3mog77qTZSX81u9XVZTNqqGg3bmDCJIJSbdsbJOIiHpf\nzmdWnKDh3e3NmDV5FICuTIG3oai7KmjikNKYrsqAjtPi9CZyT/dmIZJlJRI10HVnZTI1nH66NV39\nQraLQEREBSrnMyvTqytRVVGCmvpWLKppSJgpcM8LwLcBrTNOyz1Lt0RN92YhkmUlEjXQdQKU+Str\nfZ/+nO/YMJaIiFKR88FKZVkxbjr5E4F67rirgk6fONg32HDev/DwYXGrZ4JU3SQKZpxABkDk6c+F\n9GwfNowlIqJU5Hw1EJDaOCVOVdCimoaoQeMczjgtc5dvjdvAtqdVN05565rbIo1/85l3gDc2jCUi\nolTkRbASj1+PnenVlahv6cBbW5uwfPMerNja5NvDx9smxR2geN/r7lD73iArLE9RTnc5nExKS3Mz\nSsvKcPbMKzjQGxERBZbXwYpfBqSyrBgf7d6Hj3a34c5/bsGe1s7I+878Le2dkYyHO8hx/jpPcZ6/\nsjayLb/GuukobxB+zzzyBhqpBCDpbvTrZFD2NTdxVFoiIkpZXgcr8XrsXHj4MNz5zy34xpHD8cGu\nlpj5Wto7fYMc9417waq6SIDifa5QusubjDu4AOAbaKQSgKR7PBini/Gu2q3oV7Yfq3+IiCgleRWs\neLMH8dqyrNjahD2tnfhgV0vU+962JFOrBkU9X8htenUlWto1KzPjsGFpqS7p7jOC/IILb6CRSgCS\nqWcV9eaD/II+CJGIiMIvZ3sDNbS0Y+7yrVhbtzfSM2f+ytpIl+BEkj33x7lZL6ppwMI19Rj8hatw\nwBkzY+aZdexozDp2dMJA5YAzZuLQC65JWJ6eDgzn7nkUb+C5MLSF6U3scURElD9yNrPy8oZ6LNnY\nGBlCHwBa2k3UXzcnmAE0ExL0ScqXXHsrvnjx/4tMW/fsvEBtRIDYQOW9x27z3U462ogkCkiCrD/f\nAhr2OCIiyh85G6ycMr4Cg0qLMLVqEBbVNEQ1eC0tkpj53W1MSov6BAoKTrrshqiMihN4zP3xz5O2\nEfEGKs56/AIWp0rJeehid4KFRAFJkCqgXBlJN6jerHIiIqLMytlgZVBpUaR9iWPGYcPijlvibmMS\ntOFoy85tMdMOveAa3NjWCfz8F5haNQjPrd0VMwLtyFMv9K368Vufk9EAugaISxQsxMuAJApIgrRB\nCdtDFomIiBw522blz6t2xDy80H1T9nuIoBPMBLXu2Xm+mZDJF1+HG6/9Ae5ZuiVmBNoDzpiJyRdf\nF7PMe4/dhnXPzouZ7s5oBOlR5Mz/o1c2JX1IYirS+VDDdOvp8PxHjhmEg4cP6NayBw8fgCPHDOrW\nskRElB45m1lZsrERg0qLfDMC8ao0ulPV4QQY3kzJ5Iuvw/R3d2DBw7+ObDteY9p4gYq73EHbikyv\nroy001mwqg5TqwZFHrzoNyJvUGFus+I0lgVSH5+lenARDh4+MPJ69fbGwMtqoFIRef3m5oaUtk1E\nROmRs8HKiWMHRG6s3sBjatUgvLu9GVOron8RO1VB9S0duOeNjwN3OY4XsNxxxx2YVFkKLP1jSoGK\nNzBIpY2I8ywkZ/kfvbIp8uDFu0+fEHg9XmFus9LdxrKaURnoeq2BR5CAxRuoOOthwEJE1Ptythro\nlPEVWLCqLqoqxOkC/NzaXZHn/7i7BVeWFaO0qA9e27gHL67fHVOV4l2P+714VUKXX/8TnH3/qyll\nVNxPXe5Ol2V3lY3z4MVZk0eltA6vRN25m65+IaWqmHQ/VdlpLJvqeCl72zpiph05piJplZA3UEm0\nPiIiyrycDVZe3lAfaavi8Gv/4W7TAuhNedqEcuxfXhypSvHyLuOIF7BUjhgdM80JVNyBj/P/qVWD\nIiSUwnAAABEFSURBVE9d9ttOKpwHL7qrgLozbkuyNiupjFsSljFOVm9vxJub62OmJwpY4gUqb26u\nT6kKiYiI0idnq4Gcrst+o7a6xz9xqoKcac5gbt6xUoCu6hnvMm7xqoTc3BmV+Str8eL63ZHnDXmf\n5Oz0Xkpnm5FMVOmkUhUTpjFOnADDG4B4q4R21W5F0fZ1OPLI82LWwUCFiCi7cjZYcbouu7nbf9zz\nxseRIGHWsV2ZD3dQ4F0+0U0+Krh5dh4mnDET/YeMiCnX3p3bfKt+Glo68H7jXkweXRY1noq791K6\nAoxMdENOZdySsI1xEiRgKdq+Dpdd0LuBCh8JQEQUTM4GKw53NsQZHC5eZqKuuS3SIBVIPniaO0Bx\nZ0juuuV630AFAPoPGYEDzpgZCVic7tJvbW3CR7vbUN/SgT2tnVi7swU3nfwJ36c695Rfo90w9/bp\nDYkCloOGDURZFjIqPenlRERUSHK2zQrQFXwsXFOPmxdvimq0evrEwThzUgVmHDYs0oZj/spa1NS3\nYv/y4kh2w6ulvRPzV9ZGjeHyo1c2YZ8dUO6cS76d9Fk/h15wTWTEWidwuOq40aiqKMFVx2mDWG97\nmcqy4kgbm6CNfhPN450/XjucVKS74WxvbyteG5aykr4x03qj6ufsmVfgtK99PRTVZYWsN89rIuqe\nnM2s/HnVDvQr6oOa+lYMLOmDPa2dqKooAaCNVt/d3oxZk0dhwao6NLS0Y8nGRpw4dgDOnFSB+paO\nyND77ioi75D87jFNJg4pxYM/+yEuufbWmLLs3bktJtPiBDROhsVpCAsAVYNLY9rLONuPVxWUrB2K\nN2vkHAdn/nRkbtKZCfCrAnFPS3VbQatU4mVY3HqrjUrYqssKFTNcROHX9+abb852GVJyyy23jAJw\n+caGVkwY3A+HDOuPi48YjpK+gvMOHYq1O1vQ2tGJjQ1tWFO3F0s3N8EYg937OjGpshRXfHoUVmxp\nxPpd+9DW0Ylj9x+IsmL9ZT2oX1+sqm3Gp0aURcZgmTxmADqMwS3XX4OzZs2JKc97j92GZXfPRsPu\nPRhz5IlR7w0//AS0723Crg/eippeVtwXR48eENmuo6qiHzqMwfTqypTeA4Dfv1OLpZubUFVRgm8e\nMxLVw8qi5o+3TUADnd+/U4uqin6+77cdfyEAYNykanR0dODsmVeg/37am2ZX7VY8ce/tGFQ5FE8/\nci/GTapG//0GRKY7r72euPd2PP/4w+jo6MARx0+NmXb2zCtitpWI3/ri2dHUigMq90Nx39jEYnNr\nB/5eszPp9jIp2bGj9PI7r4nyVf2O7Xh5we8B4IGbb755S7bLE1TOBisnjh2ACz81HCeMG4QhZcU4\nevQA/O+qOvx1bQMG9euLT40ow6B+fTGpshT/dbgGM185dCjKivti4pBSvL2tCRsb2tBhDI4erReo\n/11Vh2UfN+PwEWWoHlaG379Ti+phZbjoiu/gyIuujSmLu3vyrAeexqbaepx22mlR88QLWIDYICFR\nQJHoPaArmPnmMdr9ONn8br9/pxYL19RHHQs3J1jpv98AHHH81KgLuhMkLH3pWaxa9jo6OjowblI1\nbpt1EVa88re4wYPfDcKZdtJZ5+BvTzyaUsPTVG44Bw8fgLGDy3zfK+7bB+2dndjR1Or7fm9IJfCi\nnvM7r4nyFYOVXuIEK7OPG40x5f0i0+ua2/Dm1ia0dXRi0+42dBqDd7drluWMAwejelgZFqyqQ1VF\nP1SWFePAIaVYU7cX5x06FAZ6w55aNQglfSVqfJZPn3Mpvjz75phyvPfYbVj65Fw8tGIb/vBuHT7a\n3YaP3luG40cWBc6wJAsSvBJlQPyCk2QZE0eyrI0TrPgZN6ka7y79O2o3b8LYSYfgkut/gqfn3YcV\nr/wt8trvJuB3g3CmPf3IvSnfrIPecOKNo+I2qrw0qwELf+kTUaYwWOklTrBy2sQKDOnf1eTmoRXb\n8OL63fjUiDJMqizFtsZ2NLZ2Yve+Try9rQmrapuxqGYPdu5tw2fGluN/V9Vh6eYmlPQVrNqxFwvX\n1KOkr+DSo0eirLgvqir64dPnXIof3ho7CJwTqPzolU1YsaUZja3a+HbymAEYuet9rKjZjqNOOCVq\nmeGHn4CG3Xvw4bvLIgFE9bAyNLV1oL3TYOKQUpQV900YYDjBTVNbB97Z1uQ7j3t5J+BKFgwly8Ik\nClb67zcAR598Kjo6OnDJ9T/B4GEjIzdb53WqMnWzjheoNLd2xFQJJQpYMl1Nk8u/9FmFRRRuuRqs\n5GwD24aW9shosItqGiK9dfrZpypvaWzD8LIiNLd1oKa+FTuatFdMTf0+ANHPDxpsgx5349PKsmJ8\n+tMHxWzXqfpZsKoONfWt6CtAhwFK7X1+/spavPjUTwAgpjHuxyiPaSjrDBRXWtQHlx49MqqLtNP4\n1ztYXUt7Z6CGuIka1aazK7O3oWhPG45mouFpspFp/d6P9yyhfGmQmYlxXvLl2BBRuORssHLfsq3Y\n3tQe6a0zbUI5pk0oBwBX4CLY3mxQVVGCLx00BA++uR0XHzEcdc1tuGfplsjzg5zeMk5A8NzaXQCA\nz2y+Ca9v2o1Lvz0bQPTItM5DEd/Z2oTtzR0Y3L8IL67fjWkTynHmpAr0+ecf8N5jAyO9gv7yyP9g\n11N3xgQQ3r9O2Z2/AGICGPfIt17epzjHG2AuzA8uTLcgQ+gHHekWCNcIvT2RicAiX44NEYVLzgYr\n25vaUdwH+NJBQ/DBrpaodiZOwOBkXaZWDcI9S7dgT2snHn2nFhOHlKKmvhVVFSVRQ/M7XZ6d7r+L\nN+xG+6tXYU9rJ7Z+vBntr83HB7taIuvd3dKO7c0d6F8EfPXQofjz+zvxmf0HYlFNA+Ys2ojZe/XZ\nOKVDRqDPC7+OGgBu/spaADponDP0/tzl0eM8zF2+NSYgqWtuiyzrnhZvYLx4GZRMjHLbW1LJCKTy\nrJ9EAcu/Fj2P9uEHYPCwkXnT5TgTgUW+HBsiCpecDVZKi4CWduDP7++MjF/i92wgJ4hxxmOpqW9F\nR6fBtAnlOH3i4Mg8U6sG4e2tTRheVoT9y/vhw/oWbNqtVUff+973AABFArQbRAKaAcUCANjbDvxp\nlTaynbtiG7Y0tgMA7vrnx5hYOS8mYPCO53Lp0SMjwdIJYweiqqIEBhLJfDij4Dr74iy7dmcLZk0e\nFckSOeVyxphZVNMQGWPGr1opE6PZ9sYQ8kEzAt15KGG8gOWyC87DA4/9EcijYfEZWBBRrsjZYOWK\nY0biyX/X4RPlJZHn7MR71o7TPuXCw4fh0Xd0FNtPjdwPz63dFfWQwU2727BpdxumTSjHLVPH4pE3\nt6Gmfh8GlvTFv3e04NDh/TFsv2Lsa+/ExCGl+Mz+A/HIW9tRVdEPLe1Gh9Pf2xEp46gBxZi7fCta\n2jujtjO1ahBabDWPE1g5ZRQY1NS34hPlJaiqKMHUqkFR++VUP63esRc19a2RQKWqogQXHj4Md/5z\nS9T0/WxABXQFKU55nPWlI3DRhpV3YO3K5di8/gMAmWuzEDQj0N+nwXCQAd/iBSyTTzgZq3a1p1JU\nIiJKg5wdbv8P7+7AuIpSLNnY6DuE/PTqSpw5qQLTqyuxqKYBNfWtePSdWsyaPArTJpSjoaUd723f\nCwBoaTdoae/E6AEau72/Q6cPKi3Cpt1t2NPagWkTyvGd40ajtKgPlmxsRGlRHxw9ZiBumToWg0qL\nIg1793YYfKK8GNMmlKO8tAgL19RjX3sn9i8vxr+3N2Phmno8t3YXSov6RAaeAxApY7+iPjhzUkVk\ndN5FNQ0x+1Za1AdXHTcaZ06qwKzJo3DmpArcdPInsGJrE/a0dmLUgCLsazcYNaAITW3aZmfGYcMi\n2ZuWdhMJhOavrI08psBP09UvBPo8np53H1556g/YvP4DjJ10CM6eeUXGhjF3MgLJMjdvbm7A6u17\nXK+Dj0zrHZp/9fY9DFSIiLIkZzMrWxrb0bdPSyQgAWKrN9zZCKeKZFFNA0qL+kQyC1UVJSgtEry4\nfjdOHDsAe1qbsWl313qc5T41cr/I83sAzYS4syYnjh2A/cuLUVVRiq8fORyVZcWRhrAt7Z34yFYp\nfaJcgxNv41Yns3L6xMGYWNk/biNav4axEyv7R9bx1tYm7GxuR3O7wYBiwYljB0QCKWddLe2dqKlv\nxXNrd0UCMzf3cSyNc/y91T1nz7wCLc2aGTr3W9/H4GEj8dvb52S9Z8ibmzXY29vWkfIQ+s78/Yv7\nRtYTNnxyMxEVgpwNVsqKBF8/YjiOHjMwMi1eD5fKsmLcdPInop7H41TDzDhsGABEgoo9rZ0o7gMc\nZYMTp02I023YCYKcaianF9L7O/bio91tOMIu5553bd1evL5pD5raDA4a2j+qDYrDyawsqmmIBB9+\nEjWMXVTTEAmKivsAjW0Gm3a3oqa+NdI2xmnce+LYAfj39mZ83NiO/cuLI8fB+4yhGXHK4W03MnjY\nSFw+5/aoeRJV1/TmTbYngUZvPCOoJ9hVmIgKQc4GK83tBiu2NkUFK4lu5N5uvO4HGAKI9Mh546PG\nSK+ho8cMjBtEuLe1YFUdNu1ui+pd5LaopgFNbQb7lxdHbQ+IHUPF2zvJaefibTTs177Eac8CAKdP\nHBzpCeX0EHKWf3H9blRVlOBj2xD44KH9I+tzGiM71US/vX2Ob0ARpN1IogacYb3J5lqmgl2FiagQ\n5OwItieOHYAZhw2LGnU1lefh+Ckr7ovDR5RhTd1ezJo8CkPKiuMORe/elve5PF7O+xWlOhbL29ua\nMHmMLuuMSrumbm/U8s4y7Z0Gf13bgA5jIiPtxhuRtqy4L47dfyCO3X9g5HlJzl+n7M56zzt0KDqN\nwYTB/aKOo3tfnlu7C8+98rrvsPdBR1n1jmjqvD7prHNQXNKvV4eUDzK6ahify5Oo3Lk82i0R9b5c\nHcE2Z4MV77OB0mVIWTFOP3AwhtigIUgAlGwe5/2JQ0qxYksjNja0oamtA8fuPxBVFf3w9rYm7VLt\nCkLcyzjBkvdJyt3hrHdIWXEksEkUhLUc9aUeBRTem7/zurikHy66+pZevcl6y+IXBCR6snSQIeQz\nMdx8GAMoSg8+noB6W64GKzlbDZSLKsuKcdDQ/pHxW5xp3vY03mXc1VfpGG026DgrlWXFPa6i8VZT\npFptkc5qGe+2/aqi/KquUqmy4qiwlIqwVocShQ2DlV7m17g20bD4mdCbQ+339LlB6byYe7cdNAhI\nJVjI91Fhc61NT9gxECUKRowx2S5DSkTkKADL//u0KhwwJF7HWkokWWYlquvynMW9X0CXDatX4v6b\nv4/Lb/4lxh98WFbLQoh0Rz/ta18PTQBFRMFtWLUSPzz/DAA42hizItvlCYqZlQKULJPjzrzE67rc\nW1595klsXPNvvPrMkwxWQoCZACLKBgYrFCNMDznM5s0xaJVHIVWNhKlKiogKR84Ot0+Z42Re0v2Q\nw+4IOrR+JjjtZZ6ed19a5iMiou5hZoUojkw0wCUiotQxs0IUR9CsjvNspKfn3Zf2hzYSERGDFaK0\nYFUQEVHmsBqI4mq6+oVsFyFnsCqIiChzcjFYKQWAj3bvy3Y58t7eVSuzXYScctKZ56B+Ry3qd9Rm\nuygZs3tXHV595k846axzUT44+73FiCg1mzesdf6bUwOV5eKgcDMA/D7b5SAiIsph5xtj5me7EEHl\nYrBSCeA0ADUAWrJbGiIiopxSCqAKwPPGmLoslyWwnAtWiIiIqLCwNxARERGFGoMVIiIiCjUGK0RE\nRBRqDFaIiIgo1BisEBERUagxWCEiIqJQY7BCREREofb/Ac1YUGe5mGp7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa697fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Put the result into a color plot\n",
    "Z = Z.reshape(xx.shape)\n",
    "plt.figure(1)\n",
    "plt.clf()\n",
    "plt.imshow(Z, interpolation='nearest',\n",
    "           extent=(xx.min(), xx.max(), yy.min(), yy.max()),\n",
    "           cmap=plt.cm.Paired,\n",
    "           aspect='auto', origin='lower')\n",
    "\n",
    "plt.plot(reduced_data[:, 0], reduced_data[:, 1], 'k.', markersize=2)\n",
    "plt.scatter(centroids[:, 0], centroids[:, 1],\n",
    "            marker='x', s=169, linewidths=3,\n",
    "            color='w', zorder=10)\n",
    "plt.title('Clustering on the wholesale grocery dataset (PCA-reduced data)\\n'\n",
    "          'Centroids are marked with white cross')\n",
    "plt.xlim(x_min, x_max)\n",
    "plt.ylim(y_min, y_max)\n",
    "plt.xticks(())\n",
    "plt.yticks(())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The first cluster contains customers that have vastly (around 3 times) higher spendings in the 'Fresh'-category compared to the average, indicating that those customers specialize in selling fresh products. Also, customers in this cluster tend to place many orders in the 'Frozen'- and 'Delicatessen'-Category, but relatively few in the 'Detergents and Paper'-category.\n",
    "\n",
    "Customers in the second cluster tend to spend the most overall, with particularly high spendings in the categories 'Milk', 'Grocery' and 'Detergents and Paper' and relatively low spendings in the 'Fresh' and 'Frozen' categories. Overall, this indicates that customers in this segment sell products that are more durable (i.e. not fresh).\n",
    "\n",
    "The last cluster reflects small customers that have below-average annual spendings for each of the items. Appart from the low total spending, it is apparent that the spending distribution across categories is not pathological, that is, there is no category for which spendings are particularly low or high (given that spendings are low overall).\n",
    "\n",
    "Regarding the question targeted at distinguishing the clusters visually: I generally had no problems distinguishing the clusters. Besides that, one observation is that the PCA reduction does not result in clusters that are well separated from each other. Reducing the data to three or four dimensions only (instead of two) may result in clusters that have more separation, but adds the complexity of having to visually represent the data using an (hyper-)cube instead of a plane. Of course, one could try to improve cluster representation using a 3-component PCA and a cube.\n",
    "\n",
    "*** CENTROIDS TRANSFERED TO ORIGINAL SPACE ***\n",
    "\n",
    "[\n",
    " [ 35908.28;     6409.09;   6027.84;    6808.70;    1088.15;   2904.19]  (first cluster)\n",
    "   \n",
    " [  7896.20;  18663.60;  27183.75;   2394.58;    12120.22;   2875.42]  (second cluster)\n",
    "   \n",
    " [  8276.38;   3689.87;   5320.73;   2495.45;   1776.40;   1063.97]]  (third cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "## Elbow Method\n",
    "\n",
    "### Using the elbow method to find the optimal number of clusters \n",
    "\n",
    "One of the main challenges in unsupervised learning is that we do not know the definitive answer. We don't have the ground truth class labels in our dataset that allow us to apply the techniques in order to evaluate the performance of a supervised model. Thus, in order to quantify the quality of clustering, we need to use intrinsic metrics—such as the within-cluster SSE (distortion) to compare the performance of different k-means clusterings. Conveniently, we don't need to compute the within-cluster SSE explicitly as it is already accessible via the inertia_ attribute after fitting a KMeans model.\n",
    "\n",
    "Based on the within-cluster SSE, we can use a graphical tool, the so-called elbow method, to estimate the optimal number of clusters k for a given task. Intuitively,\n",
    "we can say that, if k increases, the distortion will decrease. This is because the\n",
    "samples will be closer to the centroids they are assigned to. The idea behind the elbow method is to identify the value of k where the distortion begins to increase most rapidly, which will become more clear if we plot distortion for different\n",
    "values of k:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = df[features]\n",
    "y = df['Region']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4VeW5/vHvkxAQMKBVCwooOILFoglHjCLWoUKtWC1O\ncWgd+0O0KmLVKooTaFWgWnFo60SPpnWoip6DE9raIsWaoFbBoSpqURHFYkQqSJ7fH+/KyU7IsLOz\nd9Ye7s91rYtk7bX2fhYouXnXu57X3B0RERERaVtR3AWIiIiI5AoFJxEREZEkKTiJiIiIJEnBSURE\nRCRJCk4iIiIiSVJwEhEREUmSgpOIiIhIkhScRERERJKk4CQiIiKSJAUnERERkSQVXHAys73NbI6Z\nLTOzOjM7pJ3ndzOzO8zsZTNbZ2Z/bOaYvmZ2t5m9bmbrzWxG+q5ARERE4lJwwQnoCbwITABSWaiv\nGPgSuB54soVjugEfA1dEnyUiIiJ5oEvcBXQ2d38MeAzAzKzp62bWFZgGHA1sAvwDuMDd/xyd/yVw\nenTsSKB3M5/xLjAxOubkjFyIiIiIdLpCHHFqyyxgBHAksAtwHzDXzLaLtSoRERGJnYJTAjMbAJwA\nHOHuz7n7O+4+A5gPnBhrcSIiIhK7grtV14ZdCHOY3mhyG68r8Ek8JYmIiEi2UHBqbGPga6AMqGvy\n2hedX46IiIhkEwWnxhYRRpz6uPv8uIsRERGR7JIVc5xS6a1kZl3NbKqZLTWz/5jZ22Z2QhLn9TSz\nYWa2a7Rr2+j7Ae7+JnAPMNvMDjOzgWa2u5ldYGbfS3iPIdH53wB6R+cPa/I59Z+xMbBF9P2Q5H9X\nREREJNuYeyqtjNJchNkYYE+gGvgjcJi7z2njnIeBLYCLgLeALYEid1/Qxnn7AM+wYQ+nu9z9JDMr\nBiYDPwL6EeY2/Q2Y4u6vRu/xDrB14tsC7u7FCZ9T18xnvOvu27ZWn4iIiGSvrAhOiaLAcWhrwSkK\nWvcA27r7vzutOBERESloWXGrLgVjgReA883sX9HSJtea2UZxFyYiIiL5K1cnh28L7A38BzgU2By4\nmTDnSJ26RUREJCNyNTgVEdoFHOPuXwCY2TnAfWY2wd2/anqCmW0GjAaWEgKXiIiI5L+NgIHA4+7+\naUffLFeD04fAsvrQFFlCmKTdnzBZvKnRwN2dUJuIiIhkn2MJ86M7JFeD03zgcDPrES26C7ATYRTq\nXy2csxRg8813Zv/99+L004+nZ8+ema+0k02cOJGZM2fGXUZG5fs15vv1ga4xX+ga80O+X+OSJUs4\n7rjjIMoBHZUVwcnMegLbE0aMIOqtBKx09/fN7CpgK3f/cfT6PYSWAXeY2aWEtgTXALc1d5su8h+A\nTz6ZzX33reCVV65gwYIHKC0tzdBVxaN3796UlZXFXUZG5fs15vv1ga4xX+ga80MhXGMkLdN0suWp\nuuGErt3VhN5H04Ea4LLo9b7AgPqD3X018F1gE+DvwO+Ah4Gz2v4oo65uDEuWTGTy5OnpuwIRERHJ\ne1kx4uTuf6aVEOfuJzaz7w3CvKWU1NWNYc6cGVx/farvICIiIoUmW0acYmCsW9eDbGsAKiIiItmr\ngIOTU1KyGjNr+9AcUllZGXcJGZfv15jv1we6xnyha8wPhXCN6ZR1S65kipmVAdVhGlUZRUVzOeOM\nhVx//aUxVyYiIiKZUlNTQ3l5OUC5u9d09P2yYo5T53KKiuYyZMhMrrzygbiLERERkRxScLfqiosn\nsNVWC3nuufxrRSAiIiKZVXDBafr0m/nXvy7lT39SaBIREZH2KbjgNHIkHHggTJoEX7XUKlNERESk\nGQUXnMxgxgx4+2248ca4qxEREZFcUnDBCeBb34Lx4+Hyy2HFirirERERkVxRkMEJ4LLLoKgILr44\n7kpEREQkVxRscNp8c5gyBX7zG3j55birERERkVxQsMEJYMIE2H57OOccKJA+oCIiItIBBR2cunYN\nE8XnzYNHHom7GhEREcl2BR2cAA46SO0JREREJDkFH5zUnkBERESSVfDBCdSeQERERJKj4BRRewIR\nERFpi4JTRO0JREREpC0KTgnUnkBERERao+CUQO0JREREpDUKTk2oPYGIiIi0RMGpCbUnEBERkZYo\nODVD7QlERESkOQpOLahvT3DJJXFXIiIiItlCwakF9e0Jfv1rtScQERGRQMGpFWpPICIiIokUnFqh\n9gQiIiKSSMGpDWpPICIiIvUUnNqg9gQiIiJST8EpCWpPICIiIqDglDS1JxAREREFpySpPYGIiIgo\nOLWD2hOIiIgUNgWndlB7AhERkcKWFcHJzPY2szlmtszM6szskHacu5eZrTOzmkzWWE/tCURERApX\nVgQnoCfwIjABSPommJn1Bu4CnspQXc18ptoTiIiIFKqsCE7u/pi7X+LuDwPWjlNvAe4G/paZypqn\n9gQiIiKFKSuCUyrM7ERgEHBZHJ+v9gQiIiKFJyeDk5ntAEwDjnX3ujhqUHsCERGRwpNzwcnMigi3\n56a4+1v1u+OoRe0JRERECkuXuAtIQSkwHNjVzGZF+4oAM7O1wIHu/qeWTp44cSK9e/dutK+yspLK\nysp2F1LfnuDgg0N7gkOSfhZQRERE0q2qqoqqqqpG+1atWpXWzzDPsqESM6sDDnX3OS28bsCQJrtP\nB/YFxgFL3X1NM+eVAdXV1dWUlZWlrV53GDMmPGX36qshTImIiEh2qKmpoby8HKDc3TvcuigrbtWZ\nWU8zG2Zmu0a7to2+HxC9fpWZ3QXgweLEDfgY+I+7L2kuNGW29ob2BL/6VWd+soiIiHS2rAhOhFtv\ni4BqQh+n6UANDU/M9QUGxFNa29SeQEREpDBkRXBy9z+7e5G7FzfZTopeP9Hd92vl/MvcPX3331Kg\n9gQiIiL5LyuCUz5QewIREZH8p+CURmpPICIikt8UnNKovj3BvHmhPYGIiIjkFwWnNDvoIDjwQJg0\nCdaujbsaERERSScFpzRTewIREZH8peCUAWpPICIikp8UnDJE7QlERETyj4JThqg9gYiISP5RcMog\ntScQERHJLwpOGaT2BCIiIvlFwSnD1J5AREQkfyg4ZZjaE4iIiOQPBadOoPYEIiIi+UHBqZOoPYGI\niEjuU3DqJIntCf7xj7irERERkVQoOHWi+vYEEyeqPYGIiEguUnDqRGpPICIiktsUnDqZ2hOIiIjk\nLgWnTqb2BCIiIrlLwSkGak8gIiKSmxScYqL2BCIiIrlHwSkmak8gIiKSexScYqT2BCIiIrlFwSlG\nak8gIiKSWxScYqb2BCIiIrlDwSlmak8gIiKSOxScsoDaE4iIiOQGBacsofYEIiIi2U/BKUuoPYGI\niEj2U3DKImpPICIikt0UnLKI2hOIiIhkNwWnLKP2BCIiItlLwSnLqD2BiIhI9lJwykJqTyAiIpKd\nsiI4mdneZjbHzJaZWZ2ZHdLG8YeZ2RNm9rGZrTKz58zswM6qtzOoPYGIiEj2yYrgBPQEXgQmAMk8\nTzYKeAL4HlAGPAM8YmbDMlZhJ1N7AhERkeyTFcHJ3R9z90vc/WHAkjh+ortf5+7V7v6Wu18EvAmM\nzXixnUjtCURERLJLVgSnjjIzA0qBlXHXkk5qTyAiIpJd8iI4AT8j3O67N+5C0k3tCURERLJHzgcn\nMzsGuBg4wt0/ibuedEtsT3DjjXFXIyIiUti6xF1AR5jZ0cCvgcPd/Zlkzpk4cSK9e/dutK+yspLK\nysoMVJgeie0Jjj8ettgi7opERESyT1VVFVVVVY32rVq1Kq2fYZ5ls47NrA441N3ntHFcJfBb4Ch3\nfzSJ9y0DqqurqykrK0tPsZ3ok09ghx3g6KPh5pvjrkZERCQ31NTUUF5eDlDu7jUdfb+suFVnZj3N\nbJiZ7Rrt2jb6fkD0+lVmdlfC8ccAdwGTgL+bWZ9o69X51XcOtScQERGJX1YEJ2A4sAioJvRxmg7U\nAJdFr/cFBiQcfypQDMwCPkjYftlJ9cZC7QlERETilRVznNz9z7QS4tz9xCbf75vxorJQfXuCgw8O\n7QkOabW/uoiIiKRbtow4SZLUnkBERCQ+Ck45Ru0JRERE4qPglIMS2xOsWBF3NSIiIoVDwSlHXXZZ\nGH265JK4KxERESkcCk45Su0JREREOp+CUw5rrj1BtjU0FRERyScKTjmsvj3BvHm1HHLIFAYNOoAB\nAw5l0KADOPPMKdTW1sZdooiISF5RcMpxe+9dS8+e43j00QqWLn2SZcseZunSJ5k1q4KKinEKTyIi\nImmk4JTjJk++jjVrzgHGABbtNerqxrBkyUQmT54eY3UiIiL5RcEpxz3yyHzq6kY3+1pd3RjmzJnf\nyRWJiIjkLwWnHOburFvXk4aRpqaMdet6aMK4iIhImig45TAzo6RkNWFd5OY4JSWrMWspWImIiEh7\nKDjluLFj96Ko6PFmXysqeoxDDhnZyRWJiIjkLwWnHDd16rkMGTKDoqK5NIw8OTCX7befyZVXToqx\nOhERkfyi4JTjSktLWbDgAc44YyEDBx5Iv34/YOutD6S0dCHduz9ASUlp3CWKiIjkjS5xFyAdV1pa\nyvXXX8r114cJ42bGokVQUQFnnQW33hp3hSIiIvlBI055pn4i+G67waxZYS27u+6KuSgREZE8oeCU\nx04+GU46CcaPh5dfjrsaERGR3KfglOduvBEGD4Zx42DVqrirERERyW0KTnmue3e4/35YsQJOOAHU\nC1NERCR1Ck4FYLvtYPZseOghmK6l60RERFKm4FQgDjkELrggbM8+G3c1IiIiuUnBqYBccQWMGgVH\nHQUffhh3NSIiIrlHwamAdOkCVVVgBkcfDV9/HXdFIiIiuUXBqcD06QP33gvz58OFF8ZdjYiISG5R\ncCpAI0fCtdeG7cEH465GREQkdyg4Faizz4bDDw8tCt58M+5qREREcoOCU4Eyg9tug759Q3PML7+M\nuyIREZHsp+BUwHr1ggcegLfegtNOU3NMERGRtig4FbihQ8NCwLNnw29+E3c1IiIi2U3BSTj22DDi\n9NOfwgsvxF2NiIhI9lJwEgBmzoRhw8KE8ZUr465GREQkO6UUnMysj5n9zsw+MLOvzWx94pbuIiXz\nunWD++6D2lo4/nioq4u7IhERkezTJcXz7gS2Bq4APgQ0rTgPbLMN3H03HHQQTJsGkyfHXZGIiEh2\nSTU4jQT2dvcX01GEme0N/AwoB7YEDnX3OW2c8x1gOvAt4D1gqrvflY56CtmYMXDJJWEbMQK++924\nKxIREckeqc5xeh+wNNbRE3gRmEASo1dmNhB4FJgHDAOuB35rZvoxnwYXXxwC0zHHwPvvx12NiIhI\n9kg1OJ0NXB0FmA5z98fc/RJ3f5jkAtlpwNvufp67v+7us4D7gYnpqKfQFReHW3bdu8ORR8LatXFX\nJCIikh1SDU5/AL4DvGVmtWa2MnFLX3kt2gN4qsm+x4GKTvjsgrD55nD//VBdDeeeG3c1IiIi2SHV\nOU5np7WK9usLLG+ybznQy8y6uftXMdSUd3bfHX75Szj9dNhzTzj66LgrEhERiVdKwUmTsAvHaafB\n/Plwyinw7W/DzjvHXZGIiEh8Uh1xwsyKgUOBIdGuV4E57t4ZfZw+Avo02dcH+Lyt0aaJEyfSu3fv\nRvsqKyuprKxMb4V5wiwsyfLSS6E55vPPw8Ybx12ViIjIhqqqqqiqqmq0b9WqVWn9DPMUVnY1s+2B\n/wX6Aa9Hu3ciPG33fXd/K+WCzOpoox2BmV0NfM/dhyXsuwfYxN0PauGcMqC6urqasrKyVMsrWK+/\nDsOHw8EHwz33hEAlIiKS7WpqaigvLwcod/eajr5fqpPDbwDeAga4e5m7lxEaYr4TvdYuZtbTzIaZ\n2a7Rrm2j7wdEr19lZom3B2+JjvmFme1kZhOAw4EZKV6PtGGnneD22+H3v4cbb4y7GhERkXikGpz2\nAc5z9/97gs7dPwUuiF5rr+HAIqCa0MdpOlADXBa93hcYkPBZS4HvAwcQ+j9NBE5296ZP2kkaHXEE\nnH02TJoECxbEXY2IiEjnS3WO01dAaTP7Nwba3fXH3f9MKyHO3U9sZt+zhE7j0omuuQb+/vfQ36mm\nBrbYIu6KREREOk+qI06PAr82sxHWYA/CLbRWl0qR3FZSAn/4Q2iKecwxsF5LOouISAFJNTidSZjj\ntAD4T7TNB/4JnJWe0iRb9esHVVXw9NNw6aVxVyMiItJ5Uu3j9G/gB2a2AzA42r3E3f+Ztsokq+23\nH1x5JVx4IeyxB3z/+3FXJCIiknkp93ECcPc3gTfTVIvkmPPPD5PEjz8+zHcaODDuikRERDIr6eBk\nZjOAi919dfR1i9z9nA5XJlmvqAjuugvKy0NzzL/+FTbaKO6qREREMqc9I067ASUJX4uw6abwwANQ\nUQFnnQW33hp3RSIiIpmTdHBy932b+1pkt91g1qywnt2ee8KPfxx3RSIiIpmR0lN1Zna7mW3Qxynq\nAH57x8uSXHPyyXDSSTB+PLz8ctzViIiIZEaq7Qh+DHRvZn934EeplyO57MYbYfBgGDcO0rymooiI\nSFZoV3Ays15m1hswoDT6vn7bFDgI+DgThUr2694d7r8fVqyAE06AFNaPFhERyWrtHXH6N7CSsJ7c\nG8BnCdsnwO3ArHQWKLllu+1g9mx46CGYPj3uakRERNKrvX2c9iWMNj0NjCOEqHprgXfd/YM01SY5\n6pBD4IILwrb77jBqVNwViYiIpEe7gpO7/9nMugB3AS+4+/uZKUty3RVXwMKFcNRRoTnmllvGXZGI\niEjHtXtyuLt/DRwOFKe/HMkXXbqE9ezM4Oij4euv465IRESk41J9qu5pYJ90FiL5p08fuPdemD8/\nrGknIiKS61Jdq24ucLWZ7QJUA6sTX3T3OR0tTPLDyJFw7bVwzjmhu/hhh8VdkYiISOpSDU43Rb82\ntyado9t4kuDss+G550KLgqFDYYcd4q5IREQkNSndqnP3olY2hSZpxAxuuw369g3NMb/8Mu6KRERE\nUpPqHCeRdunVKywG/NZbMGGCmmOKiEhuSjk4mdk+ZvaImf0z2uaY2d7pLE7yy9ChcOutcNdd8Nvf\nxl2NiIhI+6W6yO9xwFPAl8AN0bYGmGdmx6SvPMk3xx0Hp50GZ5wB1dVxVyMiItI+qY44XQSc5+5H\nufsN0XYUcAFwcfrKk3w0cyYMGwaHHw4rV7Z9vIiISLZINThtCzzSzP45wKDUy5FC0K0b3HcffP45\nHH881NXFXZGIiEhyUg1O7wP7N7P/gOg1kVZtsw3cfTfMnQtXXRV3NSIiIslJtY/TdOAGM9sVeC7a\ntxdwAnBWGuqSAjBmDFxyCVx8MYwYAQccEHdFIiIirUspOLn7zWb2ETAJODLavQQ4yt0fTldxkv8u\nvhgWLIDKSli0CPr3j7siERGRlqXcjsDdH3T3ke6+WbSNVGiS9iouDrfsuneHI46AtWvjrkhERKRl\nqbYjeNvMNmtm/yZm9nbHy5JCsvnmcP/9oT3Bz34WdzUiIiItS3XEaSDNr0fXDeiXcjVSsHbfHX75\nS7jhBvj97+OuRkREpHntmuNkZockfDvazFYlfF9MeNJuaRrqkgJ02mkwfz6cckro8zRkSNwViYiI\nNNbeyeEPRb86cFeT19YRQtOkDtYkBcoMfv1reOmlsBjw88/DxhvHXZWIiEiDdt2qc/cidy8C3gO+\nWf99tHVz953c/dHMlCqFoGfPsBjw++/DqadqMWAREckuKc1xcvdB7v5J4j4z2yQ9JUmh22knuP32\nMNfpxhvjrkZERKRBqk/VnW9mRyV8fx+w0syWmdmwtFUnBeuII+Dss2HSpNDnSUREJBuk+lTdeKKl\nVczsu4SlVsYAc4FrU3lDMzvdzN4xszVm9jcz+682jj/bzF4zsy/N7D0zm2Fm3VL5bMlO11wTnrY7\n8khYsaJhv+v+nYiIxCTVJVf60rAm3cHAve7+hJktBRa2982i0avpwE+A54GJwONmtmPTW4LR8ccA\nVxGWeFkA7EiYrF4HnNvez5fsVFICf/gDlJXBkUfWMnTodTz66HzWretJSclqxo7di6lTz6W0tDTu\nUkVEpECkOuL0GTAg+noM8FT0tdF8f6e2TARudffZ7v4aYUTrS+CkFo6vAP7q7n9w9/fc/SmgCtg9\nhc+WLNavH9x2Wy1/+tM4bryxgqVLn2TZsodZuvRJZs2qoKJiHLW1tXGXKSIiBSLV4PRH4B4zexLY\njHCLDmA34J/teSMzKwHKgXn1+zzci3mKEJCa8xxQXn87z8y2BQ4C/qc9ny254YknrsPsHEJGt2iv\nUVc3hiVLJjJ58vQYqxMRkUKSanCaCNwILAa+6+5fRPu3BG5q53ttThilWt5k/3LCLcENuHsVMAX4\nq5mtBd4EnnH3X7TzsyUHPPLIfNxHN/taXd0Y5syZ38kViYhIoUppjpO7rwOua2b/zA5XlAQz+w5w\nIeGW3vPA9sANZvahu1/ZGTVI53B31q3rScNIU1PGunU9cHfMWjpGREQkPZIOTtFyK3PdfV2TpVc2\n4O5z2lHDJ8B6oE+T/X2Aj1o453JgtrvfEX3/qpltDNwKtBqcJk6cSO/evRvtq6yspLKysh0lS2cx\nM0pKVhOa1TcXjJwuXVYrNImICFVVVVRVVTXat2rVqhaOTk17RpweItw6+5iGpVea47RjgngUxKoJ\n69zNAbDwU3B/4IYWTutBeIIuUV39ud7K8+ozZ86krKws2fIkC4wduxezZj1OXd2YZl59jNWrR7Jw\nIYwY0emliYhIFmluIKSmpoby8vK0fUbSc5yiZVU+Tvi6pS2Vp+pmAKea2Y/MbDBwCyEc3QlgZrPN\nbFrC8Y8Ap5nZUWY2MOoldTkwp7XQJLlp6tRzGTJkBkVFcwm5HMApKprLttvOpF+/SVRUwPjxsHJl\nnJWKiEi+a/fkcDMrMrOTzOxRM3vFzP5hZg9HoSel+yXufi+h/9LlwCLg28Bod69ve9ifxhPFryD0\nfboCeBX4DeHJvvGpfL5kt9LSUhYseIAzzljIwIEH0q/fDxg48EDOOGMhL774ANXVpdxwA1RVweDB\ncNddWuNOREQyw9ozQBMFo0cIj/6/BLxGmHgyBNiFMOJzaAbq7DAzKwOqq6urdasux7U0Efyjj8IS\nLffcA3vvDTfdBEOHxlCgiIhkjYRbdeXuXtPR92vviNMJwChgf3ffzd0r3f1odx9GWHZlPzP7UUeL\nEmlNSwObffvC3XfDvHnw8cew225w3nnwxRfNHi4iItJu7Q1OlcA0d3+m6Qvu/jRwNXBsOgoTSdV+\n+8FLL8Fll8GvfgVDhsAf/6jbdyIi0nHtDU7fBh5r5fW5wLDUyxFJj27d4MILYfFi2HVXGDcODj4Y\n3n477spERCSXtTc4fYMNO3wnWg5smno5Iuk1aBA88gg89BC88gp861twxRXw1VdxVyYiIrmovcGp\nGPi6ldfXk2I3cpFM+sEPwujT2WfD5ZfDLrvAU0+1fZ6IiEii9oYcA+40s5b+vd6tg/WIZEzPnnDV\nVXD88TBhAnz3u3DUUTBjBmy1VdzViYhILmjviNNdhM7hq1rYPgZmp7NAkXTbeWd45hn43e/Cr4MH\nw/XXw9etjaWKiIjQzhEndz8xU4WIdCYzOO64MGH8ootg4kS48064+WbYY4+4qxMRkWzV7s7hIvlk\nk01g1ix4/nkoLoaKCvjJT+DTT+OuTEREspGCkwgwfDgsXBhC1L33htt3d9wBdU2XkhYRkYKm4CQS\nKS4Ok8Zffx3GjIGTToJRo+Dll+OuTEREsoWCk0gTffo0TBxfuRLKysIaeLW1cVcmIiJxU3ASacF3\nvgMvvghTp4ZJ40OGwP33a+kWEZFCpuAk0oquXeH882HJkjAP6ogj4KCD4J//jLsyERGJg4KTSBK2\n2SYs2zJnTghRQ4eGRYT/85+4KxMRkc6k4CTSDmPHhqVbJk0Kt/B22QWeeCLuqkREpLMoOIm0U48e\nITS9/DIMGACjR8ORR8KyZXFXJiIimabgJJKiwYNh3jy4+2549tnw/cyZWrpFRCSfKTiJdIAZHHMM\nvPYanHACnHsulJfD/PlxVyYiIpmg4CSSBptsAr/6VVi6pVs3GDkSTj4ZPvkk7spERCSdFJxE0qi8\nHBYsCH2f/vhH2Gkn+O1vtXSLiEi+UHASSbPiYhg/PizdMnYsnHpqGIF66aW4KxMRkY5ScBLJkG9+\nE+68E/78Z1i1KoxGnXOOlm4REcllCk4iGTZqVFi6Zdo0uPXW8PTdvfdq6RYRkVyk4CTSCUpK4Lzz\nQtfxESPgqKNgzBh4880Nj3UlKhGRrKXgJNKJtt46TBp/9NEQmoYOhSlT4OOPaznzzCkMGnQAAwYc\nyqBBB3DmmVOo1X09EZGsouAkEoPvfx9eeSWMQk2bVkv//uO48cYKli59kmXLHmbp0ieZNauCiopx\nCk8iIllEwUkkJj16wBVXQGXldaxbdw7uYwCLXjXq6sawZMlEJk+eHmeZIiKSQMFJJGZ/+ct8YHSz\nr9XVjWHOHLUhFxHJFgpOIjFyd9at60nDSFNTxpo1PTRhXEQkSyg4icTIzCgpWQ20FIyc5ctXM26c\nsXBhZ1YmIiLNUXASidnYsXtRVPR4s68VFT3G/vuP5NVXYY89YN994bHH1ANKRCQuCk4iMZs69VyG\nDJlBUdFcGkaenKKiuQwZMpMHH5zE4sXwwAOwejV873uw225QVQVffx1n5SIihUfBSSRmpaWlLFjw\nAGecsZCBAw+kX78fMHDggZxxxkIWLHiA0tJSiovhhz+EhQvh6aehb1845hjYcUe46SZYsybuqxAR\nKQyWLZNOzex04FygL/AS8FN3/3srx/cGpgGHAd8AlgJnu/tjLRxfBlRXV1dTVlaW5upF0sfdMWtp\nsniDRYvgmmvC8i2bbQZnnQUTJsCmm3ZCkSIiOaKmpoby8nKAcnev6ej7ZcWIk5kdBUwHpgC7EYLT\n42a2eQvHlwBPAVsDPwR2BE4FlnVKwSIZlExogobbdW+8AYcfHnpCbb01nHsuLNP/CSIiGZEVwQmY\nCNzq7rMnrG1oAAAdu0lEQVTd/TVgPPAlcFILx58MbAIc6u5/c/f33P0v7v6PTqpXJGtst124Xffu\nu3DmmfDb38KgQXDyyfDaa3FXJyKSX2IPTtHoUTkwr36fh/uHTwEVLZw2FlgA3GRmH5nZP8zs52YW\n+/WIxKVPH5g6Fd57D6ZNg7lzYeedG+ZGiYhIx2VD0NgcKAaWN9m/nDDfqTnbAkcQ6v8ecDkwCbgo\nQzWK5IxevcLtunfegd/8BrUyEBFJo2wITqkoIgSrn7j7Ine/D5hKuMUnIkC3buF2nVoZiIikT5e4\nCwA+AdYDfZrs7wN81MI5HwJrvfEjgUuAvmbWxd1b/JEwceJEevfu3WhfZWUllZWV7S5cJBfUtzI4\n7DB45hn4xS9CK4OLLgojUyeeCN27x12liEjHVVVVUVVV1WjfqlWr0voZWdGOwMz+Bix097Oi7w14\nD7jB3a9t5vipQKW7b5uw7yzgZ+7ev4XPUDsCkciiRSFA3XefWhmISH7Ly3YEwAzgVDP7kZkNBm4B\negB3ApjZbDOblnD8zcA3zOwGM9vBzL4P/By4sZPrFslJu+0Gv/99aGUwbpxaGYiIJCsrgpO730to\nfnk5sAj4NjDa3VdEh/QnYaK4u/8LGA0MJ/R8+iUwE/hFJ5YtkvO22w5uvjm0MvjpT9XKQESkLVlx\nq64z6FadSNs+/xxuvRVmzoSPPoJDD4Xzz4cRI+KuTEQkNfl6q05EskCvXvCzn4VWBr/+tVoZiIg0\npeAkIhvo1g1OOSW0Mrj/frUyEBGpp+AkIi0qLg6TxxcuhHnzQnfyY46BHXcMy7ysWRN3hSIinUvB\nSUTaZAb77QePPw7V1bD77mEy+TbbhGVePvss7gpFRDqHgpOItEtZmVoZiEjhUnASkZTUtzJYujSM\nPv3mN8m3MiiUp3lFJP8oOIlIh/TtC9Omwfvvh9t2c+fCzjuHZV4WLmw4rra2ljPPnMKgQQcwYMCh\nDBp0AGeeOYXa2tr4ihcRaScFJxFJi6atDF55paGVwYMP1lJRMY5ZsypYuvRJli17mKVLn2TWrAoq\nKsYpPIlIzlBwEpG0qm9lsGRJaGXwxRfwwx9ex6uvnkNd3RjAoiONuroxLFkykcmTp8dZsohI0hSc\nRCQj6lsZPP889Okzn7BK0obq6sYwZ878zi1ORCRFCk4ikmFOly49aRhpaspYubIHf/mLqy+UiGS9\nLnEXICL5zcwoKVkNOM2HJ6e2djWjRhklJaHdwZ57NmxbbdXJBYuItEIjTiKScWPH7kVR0ePNvlZU\n9Binnz6SmpqwuPB228GDD8IRR0C/fqHFwbHHwqxZsGiRlnsRkXhpxElEMm7q1HN5+ulxLFniCRPE\nnaKixxgyZCbTpj1AaWlYC+/008M5H3wACxbAc8/B/Plw332wbh307AkjRjSMSO2xB2y6aZxXJyKF\nxAqlEZ2ZlQHV1dXVlJWVxV2OSMGpra1l8uTpzJkzn3XrelBS8iWHHLIXV145idLS0jbPX7MmLPfy\n3HMN24oV4bWdd258e2/HHcMyMSIiNTU1lJeXA5S7e01H30/BSUQ6nbtjHUw27vDWW42D1CuvhP2b\nbdY4SA0fDj16pKl4Eckp6Q5OulUnIp2uo6EpvAdsv33YfvSjsG/VqtCtvD5ITZsGtbXQpUu4DZgY\npvr373AJIlKAFJxEJG/07g0HHhg2gPXr4dVXG4LUo4/C9deH1wYMaBykhg2DkpL4aheR3KDgJCJ5\nq7gYvv3tsI0fH/YtXx4mnc+fH8LUz34Ga9eGW3m779540vlmm6X+2em4HSki2UfBSUQKSp8+cOih\nYQP46iuoqWkYlbr99nCLD2Dw4MajUjvtBEWtNHGpra3loouu45FH5rNuXU9KSlYzduxeTJ16blIT\n4EUk+2lyuIhIAndYurTxpPOXX4a6utD2oKKiIUjtvntojwAhNFVUjGPJknOoqxtNQ8uFxxkyZAYL\nFjyg8CQSA00OFxHJILPQdLO+8SaECebPP98QpK69NkxELy4Oc6P23BPefPO6KDSNSXy3aCFjZ/Lk\n6Vx//aVxXJKIpJGCk4hIG0pLYf/9wwZh9GnJkoYg9cQT8MYb84FLmz0/LGQ84/8mpotI7lJwEhFp\np6Ii+Na3wnbqqWEi+FZb9eSjj1peyPiTT3rw4IPOnnsaffp0arkikkYKTiIiHWRmbLRR6wsZr1mz\nmh/+MLy23Xbh9l79fKmhQ8NtPxHJflrkV0QkDZJZyPi99+APf4CDD4bXX4czz4Rdd4VNNoEDDoBL\nLoHHHoN//7uTixeRpOmpOhGRNGh4qm5iswsZN/dU3Zo18MILDXOlFizQ+nsi6aan6kREslBpaSkL\nFjwQLWQ8o8lCxs23IujeHfbeO2yw4fp7CxbAbbeF/d/4RuNWCP/1Xw2tEESk82jESUQkA9LVOfzz\nzxvW31uwIGyff964FUL9fKltttGolEhTGnESEckB6VpupVcv+O53wwahFcLixQ1B6okn4MYbw2tb\nbtl40nlZGXTrlpYyRCSi4CQikkOKisJTeEOHwk9+EvZ98knDaNRzz8HFF4f5U127wvDhDUGqoiKE\nKxFJnYKTiEiO23xzGDs2bADr1sFLLzUEqfvug+nTw2uDBjWeK7XLLtAlxZ8EWshYCpGCk4hInikp\nCSNNw4fDT38a9i1b1hCk6sPUunVhgvnuuzcEqT32CBPRW6KFjKXQaXK4iEgBWrMGamoaL2b88cfh\ntcGDG086Hzw43CLUQsaSi/J2criZnQ6cC/QFXgJ+6u5/T+K8o4F7gIfc/YeZrVJEJD907w577RU2\nCC0P3nmncZC6884wGX2TTUKA+uyz61i8+BzctZCxFK6s6BxuZkcB04EpwG6E4PS4mW3exnkDgWuB\nZzNcoohIXjODbbeF446Dm26CF18MHcyfegrOOScEq4UL5+M+utnzw0LG8zu5apHOly0jThOBW919\nNoCZjQe+D5wEXNPcCWZWBPw3cAkwCujdOaWKiBSG0lLYf/+wuTv9+/fkgw9aXsh42bIeTJjgDB9u\nlJeH7uclJZ1askjGxR6czKwEKAem1e9zdzezp4CKVk6dAix39zvMbFSGyxQRKWhmRteurS9k3LXr\nav70J+OWW8II1UYbhSad5eUNm8KU5LpsuFW3OVAMLG+yfzlhvtMGzGwkcCJwSmZLExGRem0tZHzy\nySNZvBhWrYJnn4WrroIddoBnnoFTTgkLGvfqFZ7cO/10uP12ePll+PrrTr4QkQ6IfcSpvcxsY2A2\ncKq7fxZ3PSIihWLq1HN5+ulxLFnizS5kfOWVDwDhFl/iGnwAtbVh3lR1ddiefhpuvrnxyNTw4Y1H\nplLtLyWSSbG3I4hu1X0JjHP3OQn77wR6u/thTY4fBtQA62kYL64fOVsP7OTu7zTzOWVA9ahRo+jd\nu/F0qMrKSiorK9NzQSIieay2tjZayHh+k4WMJ7W7FUFimHrhhfDr6683hKldd93wNp/ClLSmqqqK\nqqqqRvtWrVrFs88+C2lqRxB7cAIws78BC939rOh7A94DbnD3a5sc2xXYvslbTAU2Bs4E3nT3DQZ+\n1cdJRCS9MtE5vD5M1QepxDDVvXvzc6YUpqQ1+drHaQZwp5lVA88TnrLrAdwJYGazgX+5+4XuvhZY\nnHiymf2bMKd8SadWLSJSwDKx3EpLt/kWLWoIUvPmhZYJzYWp4cNhyBCFKcmcrPhPy93vjXo2XQ70\nAV4ERrv7iuiQ/oCmD4qIFKDSUhg1Kmz1kglTiXOmUg1TWo9PmsqKW3WdQbfqRETyW9Mw9cIL8MYb\nDWGq6ZyplsKU1uPLL/l6q05ERKRD2hqZeuEFePJJmDWr5TDVv38te+9dvx7fpdQ/OThr1uM8/fQ4\nrccnCk4iIpK/WgtT9RPQE8NUcfF1rF9/DqD1+KR5Ck4iIlJQmgtTn38ewtRhh83ns88ubfa8urox\n3HHHDLbZBnbcMTT3HDQIunbtnLolOyg4iYhIwevVC0aNcnr06Mlnn7W8Ht+aNT2YPNlZsyYcU1wM\nAweGEFUfpup/3Xrr8LrkFwUnERERQnuFkpLW1+Pr3381b79tfPBBmHj+5pthe+ONcMvvlltg7dpw\ndNeusN12jcNU/ddbbQV6WC83KTiJiIhExo7di1mzHo+WlGmsqOgxDjlkJGbQr1/Y9t238THr18N7\n720Yqh58EN55B+rqwnE9ejQOUolfb765QlU2U3ASERGJJLseX0uKi8O8p0GDYPToxq+tXRvCU32Y\nqv/1d7+D999vOK537w3DVP3Xm2yS9ktWr6p2UnASERGJlJaWsmDBA9F6fDOarMfXsVYEXbvCTjuF\nrakvv4S33towVM2bB8uXNxy3xRbNh6rtt4eePZOvRb2qUqcGmCIiIi3IhtGYzz9vfNsv8evPPms4\nrl+/5iepb7stdOvWcFxtbS0VFfW9qkbTMKr2OEOGzMi7XlVqgCkiItJJ4g5NEJ74q2/Q2dSnn24Y\npv7+d7j7bli9OhxTVATbbNMQpl5++ToWLz4Hd/WqSoWCk4iISI7abDOoqAhbInf46KMNQ9Uzz8Cr\nr84HLm32/erqxnDPPTM47LDQTqF/f/WpakrBSUREJM+YwZZbhm2ffRr2uzv9+/fkgw9a7lX1ySc9\n2Hff0JLBLLRO2HrrMGpV/2vi1716dcYVZQ8FJxERkQJhZnTt2nqvqq23Xs2TTxrvvgvvvhvaK9T/\nunBheALw668bzujde8Mwlfhr377hdmG+UHASEREpIG31qjr00JHsuGOYD9Wc9evDbcDEUFX/9Z//\nHL6urW04vmtXGDCg5VGrAQNgo40ydLEZoKfqRERECkjDU3UTm+1VlY6n6v797+aDVf2vH37Y+Pg+\nfVoetdpmm9C/qr3z9OtbLtx//1w+/PDvoKfqREREpL0y2auq3iabhG3YsOZf/+qrcMsvMUzVB6xF\ni8L39UvXAGy8cctzrLbeOszDSlwXsHHLhUOA4R2+pnoacRIRESlg2dCrqqm6utD4s7lgVf/1v//d\ncHxxcXgCsD5ULVkyherqiqjlQg2gPk4iIiKSBtkWmiBMJq9/KnDEiOaP+fzzDYPVe++FDuyLFs3H\n/dKM1KbgJCIiIjmnVy8YOjRsidydAQN6smxZZgJhHj0gKCIiIoXOzCgpqW+5kH4KTiIiIpJXxo7d\ni6KixzPy3gpOIiIiklemTj2XIUNmUFQ0l3SPPCk4iYiISF6pb7lwxhkL2XLLCWl9bwUnERERyTul\npaVcf/2lPProzWl9XwUnERERkSQpOImIiIgkScFJREREJEkKTiIiIiJJUnASERERSZKCk4iIiEiS\nFJxEREREkqTgJCIiIpIkBScRERGRJCk4iYiIiCQpa4KTmZ1uZu+Y2Roz+5uZ/Vcrx55iZs+a2cpo\ne7K14wtJVVVV3CVkXL5fY75fH+ga84WuMT8UwjWmU1YEJzM7CpgOTAF2A14CHjezzVs4ZR/gHuA7\nwB7A+8ATZrZl5qvNboXwP0C+X2O+Xx/oGvOFrjE/FMI1plNWBCdgInCru89299eA8cCXwEnNHezu\nx7v7Le7+sru/AZxCuJb9O61iERERKTixByczKwHKgXn1+9zdgaeAiiTfpidQAqxMe4EiIiIikdiD\nE7A5UAwsb7J/OdA3yff4BbCMELZEREREMqJL3AV0lJldABwJ7OPua1s5dCOAJUuWdEpdcVm1ahU1\nNTVxl5FR+X6N+X59oGvMF7rG/JDv15jwc3+jdLyfhbti8Ylu1X0JjHP3OQn77wR6u/thrZx7LnAh\nsL+7L2rjc44B7k5L0SIiIpJrjnX3ezr6JrGPOLn7OjOrJkzsngNgZhZ9f0NL55nZecDPgQPbCk2R\nx4FjgaXAfzpYtoiIiOSGjYCBhBzQYbGPOAGY2ZHAnYSn6Z4nPGV3ODDY3VeY2WzgX+5+YXT8+cBl\nQCXwXMJbfeHuqzuzdhERESkcsY84Abj7vVHPpsuBPsCLwGh3XxEd0h/4OuGU8YSn6O5v8laXRe8h\nIiIiknZZMeIkIiIikguyoR2BiIiISE7I++BkZnub2RwzW2ZmdWZ2SNw1pZOZ/dzMnjezz81suZk9\naGY7xl1XOpnZeDN7ycxWRdtzZjYm7royycwuiP57nRF3LeliZlOia0rcFsddV7qZ2VZm9jsz+8TM\nvoz+2y2Lu650idYUbfrnWGdmv4q7tnQxsyIzu8LM3o7+DP9pZpPjrivdzGxjM/ulmS2NrvOvZjY8\n7rpSlczPezO73Mw+iK73STPbvr2fk/fBidBV/EVgApCP9yX3Bn4FjAAOIMz9esLMusdaVXq9D5wP\nlBG6zD8NPGxmQ2KtKkOiBat/QlizMd+8QpjH2DfaRsZbTnqZ2SbAfOArYDQwBJgEfBZnXWk2nIY/\nv77Adwl/t94bZ1FpdgHw/wg/NwYD5wHnmdkZsVaVfrcRnmA/FhgKPAk8lcPrvrb68z56sOwMwt+v\nuwOrCevidm3PhxTUHCczqwMOTewXlW+iSfYfA6Pc/a9x15MpZvYpcK673xF3LelkZhsD1cBpwMXA\nInc/J96q0sPMpgA/cPe8GX1pysyuBircfZ+4a+ksZvZL4CB3z5uRbjN7BPjI3U9N2Hc/8KW7/yi+\nytLHzDYCaoGx7v5Ywv4XgP9190tiKy4Nmvt5b2YfANe6+8zo+16EVUp+7O5JB/9CGHEqNJsQknZe\nrtsXDaEfDfQAFsRdTwbMAh5x96fjLiRDdoiG0d8ys/82swFxF5RmY4EXzOze6NZ5jZmdEndRmRI1\nMD6WMHKRT54D9jezHQDMbBiwF/C/sVaVXl0Iy5191WT/GvJsJBjAzAYRRkgT18X9HFhI8uviAlnS\njkDSI2oc+kvgr+6eV3NHzGwoISjV/yvpMHd/Ld6q0isKhLsSboXko78BJwCvA1sClwLPmtnQPOq/\nti1htHA6MJVwO+AGM/vK3X8Xa2WZcRjQG7gr7kLS7GqgF/Cama0nDDJc5O6/j7es9HH3L8xsAXCx\nmb1GGHk5hhAi3oy1uMzoSxhU6Mi6uICCU765CdiZ8C+jfPMaMIzwl/ThwGwzG5Uv4cnM+hNC7wHu\nvi7uejLB3RO79r5iZs8D7xLWmsyXW65FwPPufnH0/UtR6B8P5GNwOgmY6+4fxV1Imh1FCBFHA4sJ\n/6C53sw+yLMAfBxwO7CM0CuxBriHMJdUWqBbdXnCzG4EDgK+4+4fxl1Purn71+7+trsvcveLCBOn\nz4q7rjQqB7YAasxsnZmtA/YBzjKztdFoYl5x91XAG0C7n2rJYh8CTVcSXwJsHUMtGWVmWxMeSPlN\n3LVkwDXAVe5+n7u/6u53AzMJy3zlDXd/x933JUyqHuDuewBdgbfjrSwjPgKM8HBKoj7Ra0lTcMoD\nUWj6AbCvu78Xdz2dpAjoFncRafQUsAvhX7bDou0F4L+BYZ6HT3FEE+G3I4SNfDEf2KnJvp0II2v5\n5iTCbY58mvdTrwcbPpVVR57+zHT3Ne6+3Mw2JTwN+lDcNaWbu79DCEj71++LJoePoPHSbW3K+1t1\nZtaT8C/a+n+xbxtN9Fvp7u/HV1l6mNlNhDX7DgFWm1l9ml7l7nmxmLGZTQPmAu8BpYTJqPsAB8ZZ\nVzpFc3wazUszs9XAp+7edAQjJ5nZtcAjhBDRj7BE0tdAVZx1pdlMYL6Z/ZzweP4I4BTg1FbPyjHR\nCOgJwJ3uXhdzOZnwCHCRmb0PvEpohTIR+G2sVaWZmR1I+Nn4OrADYaRtMWHt2JyTxM/7XwKTzeyf\nwFLgCuBfwMPt+Zy8D06EibbPEP714IRJmxAmM54UV1FpNJ5wXX9qsv9EYHanV5MZ3yT8eW0JrAJe\nBg7M4yfP6uXbKFN/wvyJzYAVwF+BPdz901irSiN3f8HMDiNMLr4YeAc4K58mFUcOAAaQP3PTmjqD\n8EN1FuHvnw+Am6N9+aQ3cBXhHzIrCeu/Tnb39bFWlbpWf967+zVm1gO4lfAE+l+A77n72vZ8SEH1\ncRIRERHpiLy8XysiIiKSCQpOIiIiIklScBIRERFJkoKTiIiISJIUnERERESSpOAkIiIikiQFJxER\nEZEkKTiJiIiIJEnBSURERCRJCk4iknFmto2Z1ZnZt+OupZ6Z7WRmC8xsjZnVpHB+1l2TiGSegpNI\nATCzO6Mf8uc12f8DM+usRVqzbX2ny4AvCIub7t/GsS1J2zWZ2Y/N7LN0vZ+IZIaCk0hhcGANcL6Z\n9W7mtc5gbR/Szjc0K+nA6dsBf3X3f7l7qoElnddkpOHPwsyKzCztv9ciEig4iRSOp4CPgAtbOsDM\nppjZoib7zjKzdxK+v8PMHjSzn5vZR2b2mZlNNrNiM7vGzD41s/fN7IRmPmKImc2Pbo/9w8xGNfms\noWb2v2ZWG733bDPbLOH1Z8zsV2Y208xWAI+1cB1mZpdEdfzHzBaZ2eiE1+uAMmCKma03s0taeZ/z\nzOzN6H2WmtnPWzj2hKYjRk1H9Mzs22b2tJl9bmarzOzvZlZmZvsAtwO9o5HB/6vJzLqa2XVm9i8z\n+yK6vbhPwnv+OPozGGtmrwL/AQaY2XfMbGF0zmdm9hczG9Bc7SKSPAUnkcKxnhCafmpmW7VyXHOj\nHk337QdsCewNTAQuBx4FVgK7A7cAtzbzOdcA1wK7AguAR8xsU4BoJGweUE0INaOBbwL3NnmPHwFf\nAXsC41u4hrOjus4BdgEeB+aY2XbR632BxcB10XVc18L7XA2cR7itNwQ4BljewrFO2793dwPvA+WE\na7waWAfMj2r+HOjTpKZZwAjgyOha7gPmJlwLQI+ozpOBbwGfAQ8CzwBDgT2AX7dQn4i0Q5e4CxCR\nzuPuD5vZi4QgcGoH3upTdz8z+vpNMzsf6O7uVwOY2VXABcBIGgefX7n7Q9ExpwFjCD/srwPOAGrc\n/eL6g83sFOA9M9ve3f9Z/3nufkEb9U0Crnb3+6LvLzCzfQnh5Kfu/rGZfQ184e4fN/cGZrYxcCYw\nwd3/O9r9DvBcG5/dmq2Ba9z9zej7txI+bxXg7r4iYd8A4ARggLt/FO2eYWbfA04EJkf7ugCnufsr\n0XmbAr2A/3H3pdExr3egbhGJKDiJFJ7zgXlm1tIoSzJebfL9cuAf9d+4e52ZfUoYMUr0t4Rj1pvZ\nC4SRHIBhwH5mVtvkHCfMR6oPTtWtFWZmpcBWbBhw5gPteQJuCNAVeLod57RlBnCbmf2IcOv0Pnd/\nu5XjdwGKgTeazFvqCnyS8P3a+tAE4O6fmdldwBNm9mT0WfcmhC8RSZFu1YkUGHf/C+HW1dXNvFzH\nhhOem5uAva7p27awrz1/x2wMzCGEm2EJ2w7AswnHrW7He3bEmnYe3+bvnbtfBuxMuK25H7DYzH7Q\nyntuDHxNuK2X+HsyBDirtVrd/STCLbr5wFHA62a2ezuuR0SaoeAkUph+DowFKprsX0GY/5NotzR+\n7h71X5hZMWGuz+JoVw1hfs677v52ky3pEOPutcAHwF5NXtor4bOS8SZhonWyrQpWAKVm1j1h3wa/\nd+7+T3e/3t1HA38k3HIDWEsYXUq0KNrXp5nfk2ZvMTb5rJfc/RfuvhdhlPCYJK9FRFqg4CRSgKLb\nOncT5vAk+hOwRfQk2bZmdjphHlK6nG5mh5rZTsBNwCbAHdFrs4BvAL83s+HR5482s9tTeLz+WkLr\nhSPNbEczu5owUnN9sm/g7l8BvwCuMbPjo3pGmNlJLZyyEPgSuCo69hjgx/UvmtlG0ROB+5jZ1ma2\nF/BfNIS5pcDGZrafmW1mZt2juVD3ALPN7DAzG2hmu5vZBdE8p2ZFx00zsz2izzqQMHLXnuAoIs1Q\ncBIpXJcQ/g74vyet3P01YEK0vQgMJ4SQtiTzJJ4TJoxfEL33nsBYd18ZffaHhFGhIsKtxJcJc4I+\nc3dv4T1bckN07nXR+xwYfdZbCce0+V7ufjkwnTCZfjHwe2CL5t4j6gV1HPC96DOPAqYkHLse2Ay4\nizBR+/fA/wCXRucvIDyN+AfgY+Bn0XknALOja3mNMEo1HHivldK/BAYD90efdQthYv6v27pmEWmd\nNfx9JCIiIiKt0YiTiIiISJIUnERERESSpOAkIiIikiQFJxEREZEkKTiJiIiIJEnBSURERCRJCk4i\nIiIiSVJwEhEREUmSgpOIiIhIkhScRERERJKk4CQiIiKSJAUnERERkST9f+kDujqj9jIwAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa561128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "distortions = []\n",
    "for i in range(1, 11):\n",
    "    km = KMeans(n_clusters=i, \n",
    "                init='k-means++', \n",
    "                n_init=10, \n",
    "                max_iter=300, \n",
    "                random_state=0)\n",
    "    km.fit(X)\n",
    "    distortions .append(km.inertia_)\n",
    "plt.plot(range(1,11), distortions , marker='o')\n",
    "plt.xlabel('Number of clusters')\n",
    "plt.ylabel('Distortion')\n",
    "plt.tight_layout()\n",
    "#plt.savefig('./figures/elbow.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantifying the quality of clustering via silhouette plots\n",
    "\n",
    "Another intrinsic metric to evaluate the quality of a clustering is silhouette analysis, which can also be applied to clustering algorithms other than k-means that we will discuss later in this chapter. Silhouette analysis can be used as a graphical tool to plot a measure of how tightly grouped the samples in the clusters are. To calculate the silhouette coefficient of a single sample in our dataset, we can apply the following three steps:\n",
    "1. Calculate the cluster cohesion a(i) as the average distance between a sample x(i) and all other points in the same cluster.\n",
    "2. Calculate the cluster separation b(i) from the next closest cluster as the average distance between the sample x(i) and all samples in the nearest cluster.\n",
    "3. Calculate the silhouette s(i) as the difference between cluster cohesion and separation divided by the greater of the two as shown : \n",
    "\n",
    "s(i) = b(i) - a(i) / max(b(i),a(i))\n",
    "\n",
    "The silhouette coefficient is bounded in the range -1 to 1. Based on the preceding formula, we can see that the silhouette coefficient is 0 if the cluster separation and cohesion are equal (b(i)=a(i)). Furthermore, we get close to an ideal silhouette coefficient of 1 if (b(i)>>a(i)) since b(i) quantifies how dissimilar a sample is to other clusters, and a(i) tells us how similar it is to the other samples in its own cluster, respectively.\n",
    "The silhouette coefficient is available as silhouette_samples from scikit-learn's metric module, and optionally the silhouette_scores can be imported. This calculates the average silhouette coefficient across all samples, which is equivalent to numpy.mean(silhouette_samples(…)). By executing the following code, we will now create a plot of the silhouette coefficients for a k-means clustering with k=3:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmUHGd57/Hv40UZ2WCxCTSJD4sslpkRmwYwwRizXDaH\nyGGPgLATDM69jkhCCJCjEeAQgrETwCa+zoUTY5iAcwiYgG3Wiy9gSJDMoplhEcI24DEo2AgcNAjb\n7/2jauzWeEZT1dPdVd31/ZzTp6ffqq56ul/P6Of3fbs6UkpIkiRpeYdVXYAkSVK/MDhJkiQVZHCS\nJEkqyOAkSZJUkMFJkiSpIIOTJElSQQYnSZKkggxOkiRJBR1RdQFlRcRdgScDVwFz1VYjSZJqagi4\nN3BZSulnnTpo3wUnstD0gaqLkCRJfeH5wAc7dbB+DE5XAVx44YWMjIxUXIqWsnXrVs4+++yqy9Ah\n2Ef1Zx/Vm/1TbzMzM7zgBS+APDd0Sj8GpzmAkZERNm3aVHUtWsKaNWvsn5qzj+rPPqo3+6dvdHRZ\nj4vDJUmSCjI4SZIkFWRwkiSVNzsLExPZvdQgBid1xZYtW6ouQcuwj+qv1n00Owvbtzc6ONW6f9Q1\nkVKquoZSImITsGPHjh0uyquz2ai6Aknd9E3gKcClwIPytuH++vdEg23nzp2Mj48DjKeUdnbquI44\nSZJWztCkhjA4SZLKW7vjtntDkxrE4CRJKm44GZTUaP14AUxJUq8ZliTA4CRJWsiQJC3JqTpJ0m2K\nhqahIRgdze6lBnHESZKUKTPSNDoKU1Pdq0WqKYOTJDWZ03JSKU7VSVJTGZqk0gxOktREhiapLU7V\nSVKTGJikFTE4SdIgMyhJHWVwkqR+ZCCSKmFwkqR+YmCSKmVwkqSiDC23mZ6GZz8bLroou6aT1BCV\nB6eIOBV4FXDvvGkKeFNK6dLKipLUTAaj4ubmsvA0N1d1JVJPVR6cgB8Cfwl8DwjgxcDHIuIhKaWZ\nKguTNKAMSJLaVHlwSil9YkHTGyPiVcAjAYOTpJUzKEnqkMqDU6uIOAx4DnAUcEXF5UjqJ4YjST1Q\ni+AUERvJgtIQ8Evg6Smlb1dblaRaMBBJqpFaBCfg28CDgTXAs4ALIuIxhidpgBiAJA2AWgSnlNJN\nwJ784ZUR8QjgdLJP2y1q69atrFmz5qC2LVu2sGXLlq7VKalNhiZJXTQ5Ocnk5ORBbfv27evKuSKl\n+v1Bi4jPAlenlF66yLZNwI4dO3awadOm3henYmaj6gpUFUNSM8zOwnnnwStfCcPDVVcj3c7OnTsZ\nHx8HGE8p7ezUcSsfcYqIvwEuAa4B7gg8HzgJeFKVdUk6BMORhodhYqLqKqSeqzw4AXcH/hkYBvYB\n3wSelFL6XKVVSTIgSdIClQenlNLLq65BagyDkCStSOXBSVJBhh5JqtxhVRcgqQBDkyTVgiNOUt0Y\nkiSptgxOUqcYeCRp4BmcpDIMR1Jm/37YswfWr4fVq6uuRuoZ1zhJixlOi98kZWZmYOPG7F5qEIOT\ntJABSZK0BKfqJIOSJKkgg5MGn8FIktQhBif1P4ORJKlHDE7qL4YkSVKFDE6qhgFIktSHDE7qDoOR\nJGkAGZxUzGyU29/gJA22kRHYtSu7AKbUIAanflA2tEhSt61eDWNjVVch9ZzBqSqGIUmS+o7BqSoD\nPpUV01VXIEmqShqtuoLuMThJkqQVGeSgtJDfVSdJklSQwUmSJLWtSaNNYHCSJEltalpoAoOTJKkN\n6/bOsu2cCdbtna26FKmnDE6SpNKG984yce52hg1OahiDkyRJKq2J03RgcJIkSSrM4CRJklSQwUmS\nJKkgg5MkSVJBBidJklRKUxeGg8FJktSGuVVDTB03ytyqoapLkXrKL/mVJJU2s2GUjRdPVV2GKtDk\n0SZwxEmSJKkwg5MkSSqk6aNNYHCSJEkqzOAkSZJUkMFJkiSpIIOTJElaluubMgYnSZKkggxOkqTS\nRnZPs2vzGCO7p6suReopg5MkqbShA3OMfX+aoQNzVZci9ZTBSZIkqSCDkyRJOiQXht/G4CRJklSQ\nwUmSJKkgg5MkSVJBBidJkqSCDE6SpNJm1w4z8eptzK4drroUdZkLww92RNUFSJL6z3Vrh9l+2kTV\nZUg954iTJElSQQYnSZKkggxOkiRpUa5vuj2DkyRJUkEGJ0mSdDuONi3O4CRJkg5iaFqawUmSVNrQ\n3H5Gd08xNLe/6lLUYYamQzM4SZJKG9kzw9QpGxnZM1N1KVJPGZwkSZIKMjhJkiTAaboiDE6SJEkF\nGZwkSZIKMjhJkiQVZHCSJEmubyrI4CRJklTQEVUXIEnqPzPrRxj72C72HLu+6lKknjI4SZJKmxta\nzfSGsarLUIc4TVecU3WSJEkFGZwkSZIKMjhJkiQVZHCSJKnBXN9UjsFJkqSGMjSVZ3CSJEkqyOAk\nSSpt3d5Ztp0zwbq9s1WXojY52tQeg5MkqbThvbNMnLudYYNTXzI0tc/gJEmSVJDBSZKkBnG0aWUM\nTpIkSQUZnCRJkgoyOEmSJBVkcJIkSSrI4CRJKm1u1RBTx40yt2qo6lJUggvDV+6IqguQJPWfmQ2j\nbLx4quoypJ5zxEmSJKkgg5MkSVJBBidJkqSCDE6SJDWAC8M7w+AkSZJUkMFJkiSpIIOTJEkDzmm6\nzjE4SZJKG9k9za7NY4zsnq66FKmnDE6SpNKGDswx9v1phg7MVV2KluFoU2cZnCRJkgoyOEmSJBVk\ncJIkaUA5Tdd5BidJkgaQoak7DE6SJEkFGZwkSZIKMjhJkkqbXTvMxKu3Mbt2uOpSpJ46ouoCJEn9\n57q1w2w/baLqMrQI1zZ1lyNOkiRJBRmcJEmSCjI4SZI0IJym6z6DkyRJUkEGJ0mSBoCjTb1hcJIk\nqc8ZmnrH4CRJKm1obj+ju6cYmttfdSmNZ2jqLYOTJKm0kT0zTJ2ykZE9M1WXIvWUwUmSpD7laFPv\nlQ5OEXF4RDwmIu7UjYIkSZLqqnRwSindDHwKuHPny5EkSUU42lSNdqfqdgHrO1mIJElS3bUbnN4I\nnBkRT4uI4Yg4pvXWyQIlSdLBHG2qzhFtPu+T+f3FQGppj/zx4SspSpIkqY7aDU6P62gVkiRJfaCt\n4JRS+kKnC5Ek9Y+Z9SOMfWwXe451uauape3rOEXEiRFxYUR8OSJ+J2/7o4h4dOfKkyTV0dzQaqY3\njDE3tLrqUhrH9U3Vais4RcQzgcuA/cAm4LfyTWuA13emNEmSpHpZyafqTk0pvQL4TUv7l8iClCRJ\n0sBpNzjdH7h8kfZ9gFcUlySpC5ymq167wek6YMMi7Y8G9rRfjiRJWoyhqR7aDU7nA/8QEceTXbfp\ntyPi+cCZwHs6VZwkSVKdtHsdp78lC12fBY4im7b7NXBmSuldHapNkiThaFOdtDXilDJnAHcBNgKP\nBNamlP66k8VJkupp3d5Ztp0zwbq9s1WXIvVUu5cjeG9E3DGldCClNJ1S+o+U0o0RcXREvLfTRUqS\n6mV47ywT525n2ODUdY421Uu7a5xeBCx21bPVwAvbL0eSJKm+Sq1xiohjyL7IN4A7RsRcy+bDgZOB\nn3auPEmSpPoouzj852SfokvAdxfZnoBtKy1KkqSmc4qunsoGp8eRjTZ9DngmcH3LtgPA1SmlaztU\nmyRJjWRoqq9SwSml9AWAiLgPcE1KKXWlKkmSGsSg1D/aXRw+Apww/yAiTouIr0fEByPizp0pTZKk\nwWdo6i/tBqe3A8cARMQDgbOATwL3yX+WJA2wuVVDTB03ytyqoapL6WuGpv7T7pXD7wNM5z8/E/h4\nSun1EbGJLEBJkgbYzIZRNl48VXUZUs+1O+J0gOyrVgD+B/Cp/OfryUeiioqIv4qI/4iIX0TETyLi\n3yLifm3WJUmS1DXtBqcvAmdFxF8DjwA+kbffD/hRyWOdCLwLOJ4shB0JfCoiFrvApiRJA8Fpuv7U\n7lTdnwDnAs8CXpVS+nHe/lTg0jIHSimd3Po4Il5MdhHNcbKAJknSQDE09a+2glNK6RrgaYu0b11x\nRXAnsgtpXr/cjpIk9RtDU39rKzhFxD0PtT0PVu0cN4C/B76YUppebn9JkqReaneq7iqyUaGlHN7m\ncc8FRmm5RpQkSYPC0ab+125weuiCx0fmba8B3tDOASPi3WRfEnxiSml2uf23bt3KmjVrDmrbsmUL\nW7Zsaef0kiR1hWGp+yYnJ5mcnDyobd++fV05V3TyW1Mi4veAv0gpPbbk894NnAKclFLas8y+m4Ad\nO3bsYNOmTW3Xqu4KJ1qlgTaye5qLXvNsnn3WRcxsMBkcisGpGjt37mR8fBxgPKW0s1PHbXfEaSnf\nAR5e5gkRcS6wBdgM/HdE3CPftC+lNNfh+iRJHTB0YI6x708zdMA/00sxMA2mdheHL7zIZQDDwATw\nvZKHO5VsvdT/XdD+EuCCNsqTJKlShqbB1e6I08+5/eLwAH4I/GGZA6WU2r0IpyRJtWNoGmztBqfH\nLXh8C7AX2J1SumllJUmS1J8MTYOv3QtgfqHThUiS1M8MTc1QODhFxOai+6aULm6vHEmSpPoqM+L0\n0YL7Jdq/AKYkSX3H0abmKBycXMQtSZo3u3aYiVdvY3btcNWlSD1VKgxFxOMjYnqRyxEQEWsiYioi\nTuxceZKkOrpu7TDbT5vgOoOTo00NU3YU6U+B81NKv1i4IaW0DziP7GtXJEkaeIam5ikbnB4MXHqI\n7Z8CxtsvR5Ikqb7KXo7gHsBvDrH9JmBt++VIklR/jjQ1V9kRpx8DGw+x/UHAbPvlSJJUb4amZisb\nnD4JvDkihhZuiIjVwHbg3ztRmCRJUt2Unap7C/AM4LsR8W7gO3n7A4DTyK7fdEbnypMkqT4cbVKp\n4JRS+klEPAp4D/BWsi/2heyil5cBp6WUftLZEiVJdTM0t5/1P9rDnmPXMze0uupyesLQJCg/VUdK\n6eqU0snA3YDjgUcCd0spnZxS+kGnC5Qk1c/InhmmTtnIyJ6ZqkuReqqtL/kFSCndAPxnB2uRJKmW\nHG3SvLaDkyRJg87ApIX8/jlJkqSCDE6SJC3C0SYtxqk6SZJyhiUtxxEnSZIwNKkYg5MkqfEMTSrK\nqTpJUmkz60cY+9gu9hy7vupSVsTApLIMTpKk0uaGVjO9YazqMqSeMzhJkhrHkSa1y+AkSWoEw5I6\nweAkSRpIBiV1g5+qkyQNHEOTusXgJEkaKIYmdZPBSZIkqSCDkySptHV7Z9l2zgTr9s5WXQqQjTLN\n36RuMjhJkkob3jvLxLnbGa44OBmW1GsGJ0lSXzIwqQoGJ0mSpIIMTpIkSQUZnCRJfcdpOlXF4CRJ\nklSQwUmS1FccbVKVDE6SpNLmVg0xddwoc6uGenpeQ5Oq5pf8SpJKm9kwysaLp7p6DkOS6sjgJEmq\nBYOS+oHBSZJUKQOT+onBSZLUM4Yk9TsXh0uSesLQpEFgcJIkSSrI4CRJ6jpHmzQoXOMkSeoaA5MG\njSNOkqTSRnZPs2vzGCO7p5fcx9CkQeSIkySptKEDc4x9f5qhA3MHtRuWNOgMTpKkFTMwqSkMTpKk\n0naub7k3NKlBXOMkSSrF0SU1mcFJklSYoUlNZ3CSJEkqyOAkSSrE0SbJxeGSpGUsGpiGh2Hbtuxe\nahCDkyRpSUuOMg0Pw8REL0uRasGpOknSopyak27PESdJ0q0MS9KhOeIkSQIMTVIRjjhJUsMYkKT2\nOeIkSZJUkMFJkhrE0SZpZQxOkqTy9u+HqansXmoQ1zhJ0gDq+sjSzAyMj8OOHbBpU5dPJtWHI06S\nNGCcjpO6x+AkSQPE0CR1l8FJkiSpIIOTJA0IR5uk7nNxuCT1McOS1FsGJ0nqM4YlqToGJ0nqEwYm\nqXoGJ0nqA7ULTSMjsGsXrF9fdSVSTxmcJKnmaheaAFavhrGxqquQes5P1UlSjdUyNEkNZnCSpJoy\nNEn141SdJNWAIUnqD444SVLFDE1S/zA4SVKFDE1Sf3GqTpIqYGCS+pPBSZK6aGAD0uwsnHcevPKV\nMDxcdTVSzzhVJ0kqb3YWtm/P7qUGccRJkjpsYEeZJBmcJKkTDEtSMzhVJ0krZGiSmsMRJ0lqg2FJ\naiZHnCSpJEOT1FwGJ0kqwdAkNZvBSZJU3tAQjI5m91KDuMZJkhbhyNIyRkdhaqrqKqSeMzhJUgsD\nk6RDMThJaiQDkqR2GJwkNYqBSdJKGJwkNYKBSVInGJwkDQzDkaRuMzhJ6lsGJUm9ZnCS1BcMSZLq\nwAtgSqq1NGpoqqXpaRgby+6lBnHESVItGZZqbm4uC01zc1VXIvWUI06SasfQJKmuHHGS1HMGI0n9\nyhEnSZKkghxxktQ1jixJGjSOOEmSJBVkcJLUFY42SRpEBidJHWdoaoDhYdi2LbuXGsQ1TpJux+Cj\nZQ0Pw8RE1VVIPWdwkmRQkqSCDE7SgDIMSVLnucZJGkCGJknqDoOTJElSQQYnacA42iRJ3eMaJ2kA\nGJYkqTcMTlKfMBypVvbvhz17YP16WL266mqknjE4SRUxCKmvzczA+Djs2AGbNlVdjdQzBiepBwxJ\nkjQYDE7qjrHtVVdQK1F1ARoIKW2rugSp8QxOklRjhiWpXrwcgSRJUkGOOElSDTiyJPUHg5MkdZmh\nSBocBidJ6hIDkzR4DE6S1KZGB6OREdi1K7sAptQgBidJUnmrV8PYWNVVSD1ncJLUeI0eOZJUisFJ\nUl8w3EiqA6/jJEmSVJAjTpJ6ztEjSf3K4CSpZwxMkvpdLabqIuLEiLg4In4cEbdExOaqa5LUWYYm\nSYOgFsEJOBr4OvBqIFVci6QOMzQNoNlZmJjI7qUGqcVUXUrpUuBSgIiIisuR1CEGpgE2Owvbt8Pm\nzTA8XHU1Us/UIjhJ6l+GI0lNYnCSVJphSVJT1WWNk6Q+YWiS1GR9O+K0detW1qxZc1Dbli1b2LJl\nS0UVSYPNwCSpriYnJ5mcnDyobd++fV05V98Gp7PPPptNmzZVXYY0sAxKkvrFYgMnO3fuZHx8vOPn\nqkVwioijgQ3A/Cfq1kfEg4HrU0o/rK4yabAZjiSpnFoEJ+BhwOfJruGUgHfk7f8MvLSqoqRBY1BS\nxwwNwehodi81SC2CU0rpC7hQXQIMN+oTo6MwNVV1FVLP1SI4SU1lSJKk/mJwknKGGEnScgxOagyD\nkSRppVxXJEmSVJAjTqo9R4okSXVhcFItGZYkSXVkcFLXGYIkSYPCNU7qKkOTNKCmp2FsLLuXGsQR\nJ3VcxHYDkzTo5uay0DQ3V3UlUk8ZnAZUxPZKzz85OXm7L1xUvdhH9Wcf1Zv900xO1Q2AiO23u1Vt\ncnKy6hK0DPuo/uyjerN/mskRpwrVIeBIkqTiDE4dZhiSJGlwOVXXQYYmSZIGWz+OOA0BzMzMdOXg\n4+PndeW4TbNv3z527txZdRk6BPuo/mrdR/N/g7v0t7gf1Lp/1JoThjp53EgpdfJ4XRcRzwM+UHUd\nkiSpLzw/pfTBTh2sH4PTXYEnA1cBXkBEkiQtZgi4N3BZSulnnTpo3wUnSZKkqrg4XJIkqSCDkyRJ\nUkEGJ0mSpIL6IjhFxJ0j4gMRsS8iboiIf4qIo5fZ/50R8e2I+FVEXB0R/xARx/Sy7kEWEadFxA8i\nYn9EfCUiHr7M/s+OiJl8/29ExFN7VWtTlemjiHh5RFweEdfnt08v16daubK/Ry3P+8OIuCUiPtLt\nGpusjb9zayLinIi4NiLm8n+DntKrepuojT7605ZscE1EnBURv1XmnH0RnIAPAiPAE4DfAx4DHOqC\nS78NDAOvAcaAFwFPAf6pu2U2Q0Q8F3gHsA14KPAN4LKIuNsS+z+KrA/PBx4CfAz4aESM9qbi5inb\nR8BJZH30WOCRwA+BT0XEcPerbaY2+mj+efcG3g5c3uUSG62Nv3NHAp8B7gk8A7gf8Argxz0puIHa\n6KPnAW/N938A8FLgucAZpc5b90/VRcQDgGlgPKV0Zd72ZOATwLEppesKHudZwPuBo1NKt3Sr3iaI\niK8AX00pnZ4/DrJ/aN+ZUvq7Rfb/F+ColNLmlrYrgCtTSq/uUdmNUraPFnn+YcANwGkppQu7WmxD\ntdNHeb9cDvwfsv+BXJNSekaPSm6UNv7OnQr8GfCAlNLNPS22odroo3eR9c8TW9rOBB6RUnpM0fP2\nw4jT7wI3zIem3GeABBxf4jh3An5haFqZ/P+qxoHPzrelLH1/hqyvFvO7+fZWlx1if61Am3200NHA\nkcD1HS9QK+mjbcBPUkrv626FzdZm//w+cAVwbkRcFxHfioi/ysOuOqzNPvoyMD4/nRcR64GTyQZi\nCuuHr1xZB/y0tSGldHNEXJ9vW1Y+bPdGDj29p2LuBhwO/GRB+0+A+y/xnHVL7F+o/1RaO3200NvI\nphgWBl51Ruk+iohHAy8BHtzd0kR7v0PrgccDFwJPBTYA7yH7d/bN3Smz0Ur3UUppMs8DX8xHpw4H\n/jGl9LYyJ64sCUfEW/PFjUvdbo6I+3XgPHckS5O7AL+FV1pGRLwOeA7wBymlA1XXI4iIOwAXAK9I\nKd1QdT1a1GFk/2j/cUrpypTSRWRrZ06ttizNi4jHAq8n65OHkq1Fe1pEvLHMcaoccToTWG64eQ9w\nHXD31saIOBy4S75tSfkfm8uAnwPPcN65I/4LuBm4x4L2e7B0f1xXcn+tTDt9BEBE/DnwWuAJKaWp\n7pQnyvfRccC9gI/n/6cM+f/4RsQB4P4ppR90qdYmaud3aBY4kA5eODwDrIuII1JKN3W+zEZrp4/e\nBFzQMtU9leeE84C3FD1xZSNOKaWfpZS+u8ztJrI54ztFxENbnv4EIICvLnX8fKTpU8B+YLP/59wZ\nKaXfADvI+gC4dUHeE8jmjxdzRev+uSfm7eqwNvuIiHgt8AbgyQvWFKrD2uijGeCBZJ9KfXB+uxj4\nXP7zD7tccqO0+Tv0JbLpuVb3B2YNTZ3XZh8dBSxc53xLy3MLn7z2N+CTwNeAhwMnAN8B3t+y/bfJ\n/rA8LH98R+ArwNeB+5Al0PnbYVW/nn6/kU3j/Ap4IdlHOs8DfgaszbdfAPxNy/6/C/ya7PIQ9wcm\nyL6gebTq1zKotzb66C/zPnn6gt+Xo6t+LYN6K9tHizz/fcBHqn4dg3pr43foWLLZjXcC9yW7dM51\nwOuqfi2Demujj7blffRcsi//fSLwPeCDZc7bD4vDAZ4HvJtsoeotwL8Cp7dsP5LsmhlH5Y83kYUs\ngN35fZB9Eu8+wDVdrnegpZQ+nC+wexPZP65fJxul2JvvcixwU8v+V+TXzzgjv30POCWlNN3bypuj\nbB+RzfkfSfa71Wp7fgx1WBt9pB5q4+/cj/JL5ZxNdj2hH+c/L3v5D7Wnjd+hN5NliDcDvwPsJRu5\nLbXGqfbXcZIkSaoLry8hSZJUkMFJkiSpIIOTJElSQQYnSZKkggxOkiRJBRmcJEmSCjI4SZIkFWRw\nkiRJKsjgJEmSVJDBSRoQEXFLRGzOf75X/vhB+eOT8sfHVFvlYIiIR0XENyPiQER8JG87obUtf89v\nLvqeR8TnI+Ks7lYuaaX65bvqpEbLv4/pzcDJZN/JdAPZ9zK9KaV0Rb7burx93sLvU6rF9ytFxOeB\nK1NKr2lpOwn4PHCnlNIvKiuuuLOAncCTgf/O296xoO1XwHCJ1/N04DedLDIi3gesSSk9o5PHlZrM\n4CT1h4+Q/b7+EfADsvD0BOCu8zuklH664DnRs+pWbv5LuPul5uOA96SUZpdpW9gnS0op/bxTxUnq\nHqfqpJqLiDXAo4G/TCldnlL6YUrpaymlt6WU/r1lv1un6g7hYRHxnxHx3xHxpYi474JzvSoidkfE\nryNiJiJe0LLtoOm/+drytse0tG2MiE9GxC8j4rqIuCAi7pJvex9wEnB6/rybI+JewOfyp9+Qt703\n3z8i4q8iYk9E/CoiroyIZy7zfq2KiLdFxDURMRcR342Il7RsPykivppvuzYi3hoRh7VsX/Kc8+8B\ncBfgfXmtL1qk7YWLTY/m03mfz9//6yPikrx/bzdVl7+OMyPiRxFxY0RckY/MzW9/UUTcEBFPiojp\n/P2+JCLukW/fBrwIOKXlvb61nyS1x+Ak1d+N+e0PImLVCo4TwFuArcA4cBPw3ls3Rjwd+Hvg7cAY\n8L/JgsBJLcc45HRfHgI+C+wANpFNW90duCjf5XTgCuB8sqnFYeAaYD4M3TdvOz1//HrgBcAfA6PA\n2cD7I+LEQ5TxfuC5wJ8ADwBeSfb+ERG/A3wC+CrwIOBU4GXAG1uef6hzXpPX/Uvgf+W1fniRtg8t\nfL8i4iHAZ4BdwCOBE4CPA4cv8TrOAY4HngM8kOw9vCQijmvZ5yjgz4DnAycC9wTOzLedmdd2KdkI\n5TDw5SXfNUnFpJS8efNW8xvZ+pf/Ils380XgDOCBC/a5Bdic/3yv/PGD8scnATcDj23Z/6l526r8\n8RfJpppaj/kh4OOLHTNvW5O3PSZ//AbgkgXHODbfZ0P++PPAWQv2ma/vmJa2VWSB5/gF+54PXLjE\n+3Tf/FyPW2L7GcD0grZXAfvKnJNsLdkLF+xzUNvC1wR8ALj8EH186/tCFoB+A6xbsM+ngbfkP78o\nP/69F7yWa1sevw/4SNX//XrzNkg31zhJfSCl9G8R8QmyUYVHkoWe10bEy1JKF5Q41Ldafp5fi3N3\n4EfACHDegv2/RDaKUtSDgcdHxC8XtCeyNUC7SxxrA9mIyqcjonXt05HAlUs85yFkI2mXL7H9AWQj\nXq2+BNwhIo4FjmnjnEU9hGwEqIiNZCNR311QxyqyAD3vVymlq1oez5L1p6QuMThJfSKldIBsGuyz\nwBkRcT6wHSgTnFo/tTU/jVR0yv6W/H5hoGh1B+Bi4LXcfqH3LOXcIb8/Gbh2wbZfL/Gc/SXP0Ylz\nFlWmtjuQBcBN3Pa+z7ux5eeFn8LrpwX2Ul8yOEn9awY4pcPHO4FsjdC8E4Dp/Oe9+f0w8I3854dy\n8LqnncAzgKtTSgv/wZ93gNuv6zmQ37e2T5OFlXullL5Y8DV8iywInsRtC85bzeT1tXo08MuU0o8i\n4udtnLNCF73HAAACM0lEQVSob5J9EnJ7gX2vJHsv7pFS+tIKzrnYey1pBQxOUs3ln0i7iGwh9zfJ\nFiE/HPgL4KNlDrVM29uBD0XE18kWMW8mW1v1BICU0lxEfAV4XURcRbbg+M0LjncO8HLgXyLi74Dr\nydYdPRd4WUopAVcBx+efprsxpfQz4GqyAPb7EfFJYH9K6caIOBM4OyIOJ1uDtYYszO1LKb1/wblJ\nKV0dERcA742I08kC3r2Au6eULgLOJftE37uAd5NN3U2QXYOJds65jNb3963ANyPiHOAfyUaLHgt8\nOKV0/YLX8b2I+CBwQUT8OVmQujvweOAbKaVLCp7/KuBJEXE/4Gf5a7ip5GuQ1MJP1Un1dyPwFeBP\ngS+QjapsJ1uP9D9b9lvugpeLfSLu1raU0sfIPs32Z2Sf/HoF8OKU0v9r2f+lZP/D9TWyi0C+4aCD\nZdcwOoHsb8tlZEHvLOCGPDRB9mmvm8lGlH4aEfdMKV0LbAP+FrgOeFd+vL8mC2evy/e/hGwa7QeL\nvJZ5pwL/ShbiZsg+HXhUfrxr8+c/nOwCoueSLfw+o+U1FDnnId/LxdpSSt8DnkT2ab6vkq2t2kw2\nJbfY819MNg17JvBtsmt5PYzsk31FnQ98h6y/fgo8qsRzJS0ibvtbJkmSpENxxEmSJKkgg5MkSVJB\nBidJkqSCDE6SJEkFGZwkSZIKMjhJkiQVZHCSJEkqyOAkSZJUkMFJkiSpIIOTJElSQQYnSZKkggxO\nkiRJBf1/tQmRvwACf4cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9e09160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import cm\n",
    "from sklearn.metrics import silhouette_samples\n",
    "\n",
    "km = KMeans(n_clusters=3, \n",
    "            init='k-means++', \n",
    "            n_init=10, \n",
    "            max_iter=300,\n",
    "            tol=1e-04,\n",
    "            random_state=0)\n",
    "y_km = km.fit_predict(X)\n",
    "\n",
    "cluster_labels = np.unique(y_km)\n",
    "n_clusters = cluster_labels.shape[0]\n",
    "silhouette_vals = silhouette_samples(X, y_km, metric='euclidean')\n",
    "y_ax_lower, y_ax_upper = 0, 0\n",
    "yticks = []\n",
    "for i, c in enumerate(cluster_labels):\n",
    "    c_silhouette_vals = silhouette_vals[y_km == c]\n",
    "    c_silhouette_vals.sort()\n",
    "    y_ax_upper += len(c_silhouette_vals)\n",
    "    color = cm.jet(i / n_clusters)\n",
    "    plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, \n",
    "            edgecolor='none', color=color)\n",
    "\n",
    "    yticks.append((y_ax_lower + y_ax_upper) / 2)\n",
    "    y_ax_lower += len(c_silhouette_vals)\n",
    "    \n",
    "silhouette_avg = np.mean(silhouette_vals)\n",
    "plt.axvline(silhouette_avg, color=\"red\", linestyle=\"--\") \n",
    "\n",
    "plt.yticks(yticks, cluster_labels + 1)\n",
    "plt.ylabel('Cluster')\n",
    "plt.xlabel('Silhouette coefficient')\n",
    "\n",
    "plt.tight_layout()\n",
    "# plt.savefig('./figures/silhouette.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus our clustering with 3 centroids is good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 1, 2], dtype=int64)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also the regions are 3 which validates our assumption."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply different unsupervised clustering techniques using Scikit-Learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying agglomerative clustering via scikit-learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster labels: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cluster import AgglomerativeClustering\n",
    "\n",
    "ac = AgglomerativeClustering(n_clusters=3, affinity='euclidean', linkage='complete')\n",
    "labels = ac.fit_predict(X)\n",
    "print('Cluster labels: %s' % labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\hp\\appdata\\local\\programs\\python\\python35\\lib\\site-packages\\sklearn\\cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X = df[features]\n",
    "y = df['Region']\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y ,test_size=0.25, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### K Means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(330,)\n",
      "[1 1 2 2 1 2 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 1 1 1 2 1 1 1 1 0\n",
      " 1 0 1 1 1 2 1 1 1 1 1 1 2 1 2 1 1 2 1 1 1 2 1 1 2 2 1 1 1 0 1 1 1 1 1 1 1\n",
      " 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 0 1 1 1 0 1 1 1 2 1 1 0 1 1 1 1 0 1 2 1 1 2\n",
      " 1 2 1 1 1 1 2 0 2 2 1 1 2 1 1 0 1 1 1 2 1 1 2 0 1 1 1 2 2 1 1 1 1 1 1 2 1\n",
      " 2 1 1 2 1 1 1 1 1 1 2 2 1 2 1 1 1 2 1 1 1 1 1 1 2 0 2 0 2 1 1 2 1 0 1 2 2\n",
      " 1 1 1 1 2 1 2 1 2 1 1 2 1 1 2 1 1 1 2 1 1 0 1 1 1 0 0 2 2 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 2 1 2 2 1 0 1 2 1 1 1 1 1 2 0 1 2 1 1 1 0 1 1 2 2 1 2 1\n",
      " 1 1 1 0 1 1 1 1 1 1 1 2 1 1 2 2 1 1 2 2 1 2 2 1 1 2 0 1 2 1 1 1 1 2 1 1 1\n",
      " 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 2 1 2 1 1 1 2 1]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import cluster\n",
    "clf = cluster.KMeans(init='k-means++', n_clusters=3, random_state=5)\n",
    "clf.fit(X_train)\n",
    "print (clf.labels_.shape)\n",
    "print (clf.labels_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Predict clusters on testing data\n",
    "y_pred = clf.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addjusted rand score:-0.023\n",
      "Homogeneity score:0.028 \n",
      "Completeness score: 0.022 \n",
      "Confusion matrix\n",
      "[[ 0  0  0  0]\n",
      " [ 0  9  1  0]\n",
      " [ 1  5  2  0]\n",
      " [ 3 66 23  0]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "print (\"Addjusted rand score:{:.2}\".format(metrics.adjusted_rand_score(y_test, y_pred)))\n",
    "print (\"Homogeneity score:{:.2} \".format(metrics.homogeneity_score(y_test, y_pred)) )\n",
    "print (\"Completeness score: {:.2} \".format(metrics.completeness_score(y_test, y_pred)))\n",
    "print (\"Confusion matrix\")\n",
    "print (metrics.confusion_matrix(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Affinity Propogation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(33,)\n"
     ]
    }
   ],
   "source": [
    "# Affinity propagation\n",
    "aff = cluster.AffinityPropagation()\n",
    "aff.fit(X_train)\n",
    "print (aff.cluster_centers_indices_.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_pred = aff.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addjusted rand score:-0.001\n",
      "Homogeneity score:0.26 \n",
      "Completeness score: 0.052 \n",
      "Confusion matrix\n",
      "[[ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  1  0  1  2  0  0  0  2  0  0  0  0  0  0  1  0  1  1  1]\n",
      " [ 1  0  0  0  1  0  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  1  2]\n",
      " [ 3  0  3  0  4  9  5 10  1  1  2 13  4  4  9  1  2  4  5  3  2  2  5]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "print (\"Addjusted rand score:{:.2}\".format(metrics.adjusted_rand_score(y_test, y_pred)))\n",
    "print (\"Homogeneity score:{:.2} \".format(metrics.homogeneity_score(y_test, y_pred)) )\n",
    "print (\"Completeness score: {:.2} \".format(metrics.completeness_score(y_test, y_pred)))\n",
    "print (\"Confusion matrix\")\n",
    "print (metrics.confusion_matrix(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MeanShift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  8.11076522e+03   3.49355217e+03   4.81676522e+03   2.34956087e+03\n",
      "    1.67056522e+03   1.07935217e+03]\n",
      " [  1.82850000e+04   1.19360000e+04   3.25175000e+04   2.39500000e+02\n",
      "    1.37160000e+04   3.14850000e+03]\n",
      " [  8.56500000e+03   4.98000000e+03   6.72980000e+04   1.31000000e+02\n",
      "    3.81020000e+04   1.21500000e+03]\n",
      " [  3.27170000e+04   1.67840000e+04   1.36260000e+04   6.08690000e+04\n",
      "    1.27200000e+03   5.60900000e+03]\n",
      " [  2.63730000e+04   3.64230000e+04   2.20190000e+04   5.15400000e+03\n",
      "    4.33700000e+03   1.65230000e+04]\n",
      " [  2.29250000e+04   7.34980000e+04   3.21140000e+04   9.87000000e+02\n",
      "    2.00700000e+04   9.03000000e+02]\n",
      " [  7.62370000e+04   3.47300000e+03   7.10200000e+03   1.65380000e+04\n",
      "    7.78000000e+02   9.18000000e+02]\n",
      " [  1.82910000e+04   1.26600000e+03   2.10420000e+04   5.37300000e+03\n",
      "    4.17300000e+03   1.44720000e+04]\n",
      " [  1.61170000e+04   4.61970000e+04   9.27800000e+04   1.02600000e+03\n",
      "    4.08270000e+04   2.94400000e+03]\n",
      " [  3.59420000e+04   3.83690000e+04   5.95980000e+04   3.25400000e+03\n",
      "    2.67010000e+04   2.01700000e+03]\n",
      " [  8.50000000e+01   2.09590000e+04   4.58280000e+04   3.60000000e+01\n",
      "    2.42310000e+04   1.42300000e+03]\n",
      " [  4.44660000e+04   5.42590000e+04   5.55710000e+04   7.78200000e+03\n",
      "    2.41710000e+04   6.46500000e+03]\n",
      " [  3.68470000e+04   4.39500000e+04   2.01700000e+04   3.65340000e+04\n",
      "    2.39000000e+02   4.79430000e+04]]\n"
     ]
    }
   ],
   "source": [
    "ms = cluster.MeanShift()\n",
    "ms.fit(X_train)\n",
    "print (ms.cluster_centers_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_pred = ms.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addjusted rand score:0.035\n",
      "Homogeneity score:0.1 \n",
      "Completeness score: 0.1 \n",
      "Confusion matrix\n",
      "[[ 0  0  0  0  0  0  0]\n",
      " [10  0  0  0  0  0  0]\n",
      " [ 5  3  0  0  0  0  0]\n",
      " [80  3  0  0  4  4  1]\n",
      " [ 0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "print (\"Addjusted rand score:{:.2}\".format(metrics.adjusted_rand_score(y_test, y_pred)))\n",
    "print (\"Homogeneity score:{:.2} \".format(metrics.homogeneity_score(y_test, y_pred)) )\n",
    "print (\"Completeness score: {:.2} \".format(metrics.completeness_score(y_test, y_pred)))\n",
    "print (\"Confusion matrix\")\n",
    "print (metrics.confusion_matrix(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixture of Guassian Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adjusted rand score for covariance=spherical:-0.016\n",
      "Adjusted rand score for covariance=tied:-0.08\n",
      "Adjusted rand score for covariance=diag:0.029\n",
      "Adjusted rand score for covariance=full:0.055\n"
     ]
    }
   ],
   "source": [
    "from sklearn import mixture\n",
    "\n",
    "# Define a heldout dataset to estimate covariance type\n",
    "X_train_heldout, X_test_heldout, y_train_heldout, y_test_heldout = train_test_split(\n",
    "        X_train, y_train,test_size=0.25, random_state=42)\n",
    "for covariance_type in ['spherical','tied','diag','full']:\n",
    "    gm=mixture.GMM(n_components=3, covariance_type=covariance_type, random_state=42, n_init=5)\n",
    "    gm.fit(X_train_heldout)\n",
    "    y_pred=gm.predict(X_test_heldout)\n",
    "    print (\"Adjusted rand score for covariance={}:{:.2}\".format(covariance_type, \n",
    "                                                               metrics.adjusted_rand_score(y_test_heldout, y_pred)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GMM(covariance_type='tied', init_params='wmc', min_covar=0.001,\n",
       "  n_components=3, n_init=1, n_iter=100, params='wmc', random_state=42,\n",
       "  tol=0.001, verbose=0)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gm = mixture.GMM(n_components=3, covariance_type='tied', random_state=42)\n",
    "gm.fit(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addjusted rand score:-0.044\n",
      "Homogeneity score:0.015 \n",
      "Completeness score: 0.025 \n",
      "Confusion matrix\n",
      "[[ 0  0  0  0]\n",
      " [ 0  9  1  0]\n",
      " [ 0  8  0  0]\n",
      " [ 1 83  8  0]]\n"
     ]
    }
   ],
   "source": [
    "# Print train clustering and confusion matrix\n",
    "y_pred = gm.predict(X_test)\n",
    "print (\"Addjusted rand score:{:.2}\".format(metrics.adjusted_rand_score(y_test, y_pred)))\n",
    "print (\"Homogeneity score:{:.2} \".format(metrics.homogeneity_score(y_test, y_pred)) )\n",
    "print (\"Completeness score: {:.2} \".format(metrics.completeness_score(y_test, y_pred)))\n",
    "\n",
    "print (\"Confusion matrix\")\n",
    "print (metrics.confusion_matrix(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl=plt\n",
    "from sklearn import decomposition\n",
    "# In this case the seeding of the centers is deterministic,\n",
    "# hence we run the kmeans algorithm only once with n_init=1\n",
    "\n",
    "pca = decomposition.PCA(n_components=2).fit(X_train)\n",
    "reduced_X_train = pca.transform(X_train)\n",
    "\n",
    "# Step size of the mesh. Decrease to increase the quality of the VQ.\n",
    "h = .01     # point in the mesh [x_min, m_max]x[y_min, y_max].\n",
    "\n",
    "# Plot the decision boundary. For that, we will asign a color to each\n",
    "x_min, x_max = reduced_X_train[:, 0].min() + 1, reduced_X_train[:, 0].max() - 1\n",
    "y_min, y_max = reduced_X_train[:, 1].min() + 1, reduced_X_train[:, 1].max() - 1\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "\n",
    "gm.fit(reduced_X_train)\n",
    "#print np.c_[xx.ravel(),yy.ravel()]\n",
    "Z = gm.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "\n",
    "# Put the result into a color plot\n",
    "Z = Z.reshape(xx.shape)\n",
    "pl.figure(1)\n",
    "pl.clf()\n",
    "pl.imshow(Z, interpolation='nearest',\n",
    "          extent=(xx.min(), xx.max(), yy.min(), yy.max()),\n",
    "          cmap=pl.cm.Paired,\n",
    "          aspect='auto', origin='lower')\n",
    "#print reduced_X_train.shape\n",
    "\n",
    "pl.plot(reduced_X_train[:, 0], reduced_X_train[:, 1], 'k.', markersize=2)\n",
    "# Plot the centroids as a white X\n",
    "centroids = gm.means_\n",
    "\n",
    "pl.scatter(centroids[:, 0], centroids[:, 1],\n",
    "           marker='.', s=169, linewidths=3,\n",
    "           color='w', zorder=10)\n",
    "\n",
    "pl.title('Mixture of gaussian models on the seeds dataset (PCA-reduced data)\\n'\n",
    "         'Means are marked with white dots')\n",
    "pl.xlim(x_min, x_max)\n",
    "pl.ylim(y_min, y_max)\n",
    "pl.xticks(())\n",
    "pl.yticks(())\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot of the data fitted by Guassian models looks like :\n",
    "    \n",
    "<img src=\"images/gmplot.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "PCA in combination with K-means-clustering was the most insightful solution technique. Even though only six features are present in the data, PCA identified the the two features with the highest variance and enabled me two cluster the data in a lower dimensional space. We used PCA over ICA for dimensionality reduction, since the PCA eigenvalues provide a convenient ordering of the most important variance directions (eigenvectors) of the data that can later be used in clustering. ICA does not provide a similar rank order for the most important components. After preprocessing the data with PCA, K-means uncovered (clustered) interpretable customer groups that are helpful to the client. We'hard-clustered' the data with K-means (instead of GMM) because there is no clear indication about the underlying data generating process and hence I prefered a method that does not assume Gausssian distributions. \n",
    "Due to the clustering, we can now tell the client that, among the large customers, one group focusses on the 'fresh' category and the other on more durable products. The small customers, however, have a more balanced assortment (and we also concluded that small customers actually exist.). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How would this help the company?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the client plans to change, retire or introduce new (delivery) strategies in the future, she could use the segmentation k-means found to try the strategies on subsets of the segments. That is, the client could choose a small number of customers in each segment and introduce the change for those customers only. It could then, after some time, see how the customer segments respond to the change (for example by asking them for feedback). Based on this, the client may conclude to not proceed with the new strategy, proceed only for one or two customer segments, or implement it for the whole customer base. The ability to receive feedback on a new strategy from a small subset of customers before making a final decision is advantageous for the client."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How would we use that data to help predict future customer needs?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the client implemented a strategy for a subset of each of the customer segements and asked them to provide feedback on a, say, 1 (very poor) to 5 (great) scale, the client could then run a linear regression of the feedback provided on a set of cluster membership features (dummy variables) to identify the needs of each of the customer segments (note that this linear regression basically results in segment averages for the response). This will help the client to determine whether the new strategy suits the customers needs (by segment) or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "___"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}