{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Clustering\n", "\n", "A common task or goal within data analysis is to learn some kind of structure from data. One way to do so is to apply clustering analysis to data. This typically means trying to learn 'groups' or clusters in the data. \n", "\n", "This example is a minimal example of clustering adapted from the `sklearn` tutorials, to introduce the key points and show an introductory example of a clustering analysis. \n", "\n", "As with many of the other topics in data analysis, there are many resources and tutorials available on the clustering analyses. A good place to start is the extensive coverage in the `sklearn` documentation. If you are interested in clustering analyses, once you have explored the basic concepts here, we recommend you go and explore some of these other resources. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Clustering is the process of trying to find structure (clusters) in data.\n", "
\n", "\n", "
\n", "Clustering\n", "article from wikipedia. The sklearn \n", "user guide \n", "has a detailed introduction to and tutorial on\n", "clustering. \n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Imports\n", "%matplotlib inline\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sklearn import cluster, datasets\n", "from sklearn.cluster import KMeans\n", "from scipy.cluster.vq import whiten" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load an Example Dataset\n", "\n", "Scikit-learn has example datasets that can be loaded and used for example.\n", "\n", "Here, we'll use the iris dataset. This dataset contains data about different species of plants. It includes information for several features across several species. Our task will be to attempt to cluster the data, to see if we can learn a meaningful groupings from the data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Load the iris data\n", "iris = datasets.load_iris()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `iris`, as loaded by `sklearn` is an object. The data is stored in `iris.data`. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There are 150 samples of data, with 4 features and 3 labels.\n" ] } ], "source": [ "# Let's check how much data there is\n", "[n_samples, n_features] = np.shape(iris.data)\n", "n_labels = len(set(iris.target))\n", "print(\"There are {} samples of data, with {} features and {} labels.\".format(\\\n", " n_samples, n_features, n_labels))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sepal length (cm)\n", "sepal width (cm)\n", "petal length (cm)\n", "petal width (cm)\n" ] } ], "source": [ "# Check out the available features\n", "print('\\n'.join(iris.feature_names))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "setosa\n", "versicolor\n", "virginica\n" ] } ], "source": [ "# Check out the species ('clusters')\n", "print('\\n'.join(iris.target_names))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Let's set up some indexes, so we know what data we're using\n", "sl_ind = 0 # Sepal Length\n", "sw_ind = 1 # Septal Width\n", "pl_ind = 2 # Petal Length\n", "pw_ind = 3 # Petal Width" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEcCAYAAADDfRPAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7wcdX3/8dc7hwDhIjESEEIgSPmhSIzYSMBQDVK8oD9NEYQIFmgrrZaK8hMrlSpa8NKoFa8UFAHBaAWMaKlIlYvcogkXoyIKEsiFS7iEa4AQPr8/ZjaZ3TN7zszZ3dndc97Px2Mf2f3Od2Y+Mzm737l8P/NVRGBmZlYzrtsBmJlZb3HDYGZmddwwmJlZHTcMZmZWxw2DmZnVccNgZmZ13DDkkHSqpEhfp5aY79zMfNPaFMuyzDKfl/S0pHslXSfp45K2b3H5c9PtPVXSxHbEXHC9czLblX09Jul6Sce2uOzaNk1rU4ynFqi/YTtGus4qSTqmtp9yppXa9l4j6Y2Z+P+zYdo5mWmnNUy7NDNtz7Ss9h1cVmC92d+OYzLlTb9nkqZl5jm3hc1um026HYCVImAz4MXp6zXACZIOj4grRrjMucDR6ftzgTWtBtmirYH9gP0kvSIiPjiCZcwBPp6+vwpY1pbIRp9jgNel70/tXhgdcQPwPMnB7+yGaa/JvG827RHgtjbG02vfsyH5jKENJG0OEBHHRITS17IOrGpXkoZhOskfF8ALgR9IelkH1leVqyNCwOYkP1Y1J7TrzMvGloh4DPhN+nFPSS8EkLQtsEem6j6SNkmnvRR4UVp+faTZvxExLf1OT6sk+B7ghqGEhktFfyHpIkmPkh5ZNLuUJOk9khZLeljSM5JWSrpC0tFNVtVURDwbEb+JiGOBH6bFW7LxCBlJh6fLXy7pqXSdf5J0ZvbSU3rJIxvDXdn4JW0l6TxJSyU9JGmdpDWSrpF0eM7+qc27rOx2pdv2TEScByytLRKYmVn+vpJ+IOn+NJZV6T6flqmzLLsvgCszcc1J63xd0k2SVqfLeVzSLyW9V5JGEntZkt6V7sdH0/+fP0g6XdIWDfVqsV8l6WBJv5K0VtKdkj7cGG9a51Yllxxvl/Tuhr/LObVLF2w8Wxj2Mpik49MYn5J0i6Q3D7N9kyU9my7zJw3T3pxZ3yfTshmSLkm/G8+kf283SfpPSeNL7t6a62qrJDkLhY1nBIuB1cAWwN5pWfbs4fpMvLmXkiTtJ+mGdF8vk3RiXhDDfc9y6r9T0q/T/+fbRvI70bKI8KvhRXJaHenr1Ez5uZnyBzPvl+VMn5aWHZYpa3xdVCCWZY3LzEzbLzPtMWBcWn7mEOv8PbBp7WBoiNc0kstVQ9U5uiGeuv0xzHbNydS/qmHabzLTDk3L3gk81ySOh4A9cvZX42tOWufpIep8vEmMpxbYpg3LGabel4dY/6+ACTnLfJTk0khj/aMa4s3bRyuz+yD9v236/5qz7ffl1HsWeMkw23lJWncdMDlTfn5a/jzJWfAWJD/SzWLaaoTf4yMzyzg9LftM+vmLwA/S9yek087J1H9dzndwWaZsT+DJnFhXZd4fU/B7lv3/yNvXAexf5W+gzxhG7jGSH+YJwMFD1Htt+u8TJKewmwG7kPzQ/aTZTAX9PvN+azaeBn8HmAVsC4wHtge+lU7boxZvJJdvzsssY9eovxT2OHA4yR/uFiSXel4DPJXWH8n1/6YkbZoeHb08LQrgV+lR9NeBAeAm4KUk+/EAkh+oScD8dJumAZ/ILPaAzDZdlZb9DbA7yT7bFHgFsCKddkInzxok7Qscn348l6Tx3QI4KS2bCbw3Z9YXAJ8muXR4fKb83Zn3p5PsI0jOmrYB5gE7ZhcUEcvS//urM2W1fZS37ZOAQ4CJwIVp2XiSv42hnJP+uwnJ3zuSJpBcb4fkgOAu4GUkf6sAHyb5O5sM7J9u83PDrKeZ6zLv90//nZ2Zdl1DWa3OOpIGeij/SvL/BvA1kn3zepJ9XqfA9yxre+B96fI+myl/N1WqshXqlxfFzhjelTNfdvq0tOxENh4dnQ+cALwB2KZgLMsal5mZ9iLqjyomp+W7A98G7iH54Ww8+vjnoWLOTBPJH+mN5B+xrm1hH8/JiavxdUZa96ACdddmlp39/5uTs+7DSW5KPwysz1nW9jkxnlpgm+qOupvUOb3AtvwkZ5n3AQNp2VaZ8t+nZVtktuWhWt102nV5+yPdB7nxNmz7RZnyt2bKzxxmfwyw8Qj62rQsewZ9VFq2HcmPcQBLgI+RNCR/1obv8op0uU+l+21t+nlHYN/0/UqShqgW1y+bfAeXZcruz9TfJlP+7Uz5MQW/Z9My0xZnyvfK+5uo4uUzhpG7uWC9rwHfJ/lRfTfJKezlwP2SPtJiDC/NvH8MeEjSNsC1wFHAVJIju0YTCi7/n4Gvkpx9vICkocjavFS0xTwJLAKOAz6Qlm1XYL7NJW05XCVJRwDfJbm+/kLy77MV3T8jUWRbXpRTdmdErE/fP5kpr/0fTGLjtqzM1IXkAKEVt2fe5607VxrD+enH10jaBXhX+vlR4OK03gPAP5Bcnn0VyRnf94A/SvqFpBe0EHvtXsEE4O/SmJdFxCqSRmgtSSNxZGae6xhe7f/o8Yh4NFO+Iq9yCSPa1+3mhmHk1hapFBFPR8Q7Sb64+5NcxlhEcinkU5KmtBDDyZn3/x0Rz5NcXqn9+PwM2CGSU9n3NwtxiOUfkXk/F9gsXdZDI4y3matj46n1VhGxb0ScHelhE/BApu7ZmbrZyx/jIqL2RSq6Tf9Ecj1fJJeoqpDdliObbMs+OfOtq73J7Jesh0kOPgB2kJT9bk9tEstQ+yl33SXmqflW+q9IfvxrN60XRMSG71BEfJPkstp04B3Al9JJ+wP/WHKdWdkf+dqlz+vTdWYvGX2wyTzNPJj+u3V6MFazU5P6VezrtnHD0GGS3iHpeGAKcCvJ2cOttck0/0NqtrzxkvZSkgjzlrT4SeCT6fvs9dingSclvZzkRzBP9kd+RsP19eyy1gDjJf0r+Ue0LfdKGsL1JP3KAY5Oe/RsnfZ8mS1pPsmZWE12m17R8COZ3abHkrB1LBt7prSFpDflvLYDfpypdloa/+aSpqa9db5D/dFrIRHxFMklP0iu138o3UdHUN9vP2vDfpL0yrLrLBjX7Ww8av8QyQERbLz/gKRtJX2O5J7dA8CP2NjjDmDnTN2yf2PZH/mdc8qubVxHJt6hXJl5/ylJ20g6gOReTJ6hvme9p8rrVv3yotg9hmk58w2aDpySKWt8rSLTA6VJLMuGmD9IjhQPytR/IcmXq7HeH5ps06E5dZel0z6aM201yY/0oGvTjfMPs11zMvWvKlB/Hvn3A2qvczN1Z+bVSacdmTPtKWB5zv9dNsZTC8Q41P9TAHPTel8bpt4xOcu8qsm6ljXs07xeSdmeMq/L1P9QTt2rhtr2hvJzh9sn6Tx/27COpQ3Tdyqy38r+jaX1NyHp+JFd3ozM9IMbpt01xHcwu6+b9UrK9q7K/j8O9T2blrdPG8qH/Y608+Uzhs77GUkvoTtI/kDXA/eSXueOzOl0AQE8Q3Ij8nqSnicvi0zWc0Q8QnK6fi3JD94qkobuM02WeTHJDdF70tiyPgt8iuTm3FqSXiyvJ7k+XKmIWEByWeFikht/z5F8CRencX4+U3cxyaWzO6k/NSciLiS5bHAXyRnVYpL9dWfHN2JjDO8juQd0Ncm+XEdybfpKkl45/zPC5V4FvJ0kD+RZkoOBo6m/TJY9cv0qSdfme6Gjly2+R/318nMapj8C/AfJZZ0HSf4OHyf5G39XRCwc6Yoj4jngl5mix9mYJ0O6jucbPhdZ7u9IOkUsItnXy0kOpL7WZJahvmc9R2nLZGZ9Lk0EOwC4MpLr50h6E7CQ5BLOvcBOkdyLMmvKz0oyGz02I+nxtk7S/SR5GrUbo88B/+BGwYrwpSSz0eMZkkSqu0l6wU0guXRxAfDqiLi0i7FZH/GlJDMzq+MzBjMzq9P39xi23XbbmDZtWrfDMDPrK0uWLHkwIibnTev7hmHatGksXry422GYmfUVSXc3m+ZLSWZmVscNg5mZ1XHDYGZmddwwmJlZHTcMZmZWxw2DmZnVqay7qqSpJKM5vZjkaYZnRcQZDXXmkDyH/a606JKI+CRmZjkW3ryS+Zffzqo1a9lx4gROeuMezN272NhXefMCI17eaFJlHsNzwP+LiJskbQ0skXRF+vjarF9ExFsrjMvM+tDCm1dy8iVLWbsueYr1yjVrOfmS5Inaw/2Y58170vdvBcG69VF6eaNNZZeSIuLeiLgpff84cBvJqGZmZqXNv/z2DT/sNWvXrWf+5bc3mWPoedc9HxsahbLLG226co9B0jSSoRQX5UzeT9Ktkv4nHZIyb/7jJC2WtHj16tUdjNTMetWqNfljXDUrL1tnJHVHi8obBklbkYxm9IGIeKxh8k3ALhExA/gyyQAjg0TEWRExMyJmTp6c+6gPMxvldpw4oVR52TojqTtaVNowpCNMXQxcGBGXNE6PiMci4on0/WUkg89vW2WMZtYfTnrjHkwYP1BXNmH8wIabyGXnHT9OjB/QiJY32lTZK0nAN4HbIuILTeq8GLg/IkLSPiQN10N5dc1sbKvdEB5JL6Jm8450eaNNZQP1SNof+AXJQNy14QX/BdgZICLOlHQ88F6SHkxrgRMjYsjBuWfOnBl+uqqZWTmSlkTEzLxplZ0xRMS1gIap8xXgK9VEZGZmefp+PAYz665WksyKOmXhUhYsWs76CAYk5s2aymlzp7d1HbaRGwYzG7FWksyKOmXhUi648Z4Nn9dHbPjsxqEz/KwkMxuxVpLMilqwaHmpcmudGwYzG7FWksyKWt+kg0yzcmudGwYzG7FWksyKGlB+n5Vm5dY6NwxmNmKtJJkVNW/W1FLl1jrffDazEWslyayo2g1m90qqTmUJbp3iBDczs/J6IsHNzEanogPeFC0rc7ZRNIei3fW6qYoYfcZgZiPWmMcAycPosgPeAMnD6SIZ82CoehPGD/DpQ6YX+qHLW3fe/O2u103tjHGoMwbffDazESs64M269VHXKDSrVyYHomgORbvrdVNVMbphMLMR68QgNkWXWTSHot31uqmqGN0wmNmIdWIQm6LLLJpD0e563VRVjG4YzGzEig54M35AyT2FYeqVyYEomkPR7nrdVFWM7pVkZiNWZsCbomVFb6IWzaFod71uqipG90oyMxuDnMdgZn2tTN/9fshFyNNLcbthMLOeVmbMhyrGh+iEXovbN5/NrKeV6bvfD7kIeXotbjcMZtbTyvTd74dchDy9FrcbBjPraWX67vdDLkKeXovbDYOZ9bQyfff7IRchT6/F7ZvPZtbTyvTd74dchDy9FrfzGMzMxiDnMZhZT2plLIe5e0/p23EWei2eRj5jMLOuyB3LocS4De/48ylcvGRl342z0CvxeDwGM+s5uWM5lBi3YcGi5X05zkKvxZPHDYOZdUWrffTXN7na0evjLPRaPHncMJhZV7TaR39Ayi3v9XEWei2ePG4YzKwrcsdyKDFuw7xZU/tynIVeiyePeyWZWVe0OpbD3L2nMHOXSX03zkKvxZPHvZLMzMagnshjkDQVOB94MfA8cFZEnNFQR8AZwMHAU8AxEXFTVTGaWaJZP/tW8w5Gu1byE3opt6GyMwZJOwA7RMRNkrYGlgBzI+J3mToHA/9E0jDMAs6IiFlDLddnDGbt1ayffV7eQJm8g27lDVSllfyEbuQ29EQeQ0TcWzv6j4jHgduAxi1+O3B+JG4EJqYNiplVpFk/+7y8gTJ5B73UT78TWslP6LXchq70SpI0DdgbWNQwaQqwPPN5BYMbDyQdJ2mxpMWrV6/uVJhmY1Kz/vTN8gZaXe5o0Up+Qq/lNpRqGCTtKukASQdLerWkzcuuUNJWwMXAByLiscbJObMM+muMiLMiYmZEzJw8eXLZEMxsCM360zfLG2h1uaNFK/kJvZbbMGzDIGmapM9Kuge4A/gZ8GOSo/01kq6QdJikIssaT9IoXBgRl+RUWQFMzXzeCVhVYDvMrE2a9bPPyxsok3fQS/30O6GV/IRey20Y8sdc0hnArcBLgI8CewLbAJuS9C46GLgW+Dfg15JePcSyBHwTuC0ivtCk2qXAXyuxL/BoRNxbbpPMrBVz957Cpw+ZzpSJExAwZeIEPn3IdE6bO31Q+fxDZzD/sBn1ZYfNYP6hMwbNP5pvPEPz/VZku1uZtxOG7JUkaT7w2Yh4cNgFJT2KtoiIi5pM3x/4BbCUpLsqwL8AOwNExJlp4/EV4E0k3VWPjYghuxy5V5KZWXkjzmOIiJOKriQiLhtm+rXk30PI1gngH4uu08yqldfXfvHdD7Ng0XLWRzAgMW/WVE6bO73QvL12FtEPMVbBj8Qws0Ia+9qvXLOW//f9W1mf6a66PoILbrwHoK5xyJv35EuWAvTMD28/xFiVwr2SJL1Q0hmSfi3pPkkPZF+dDNLMui+vr/365/MvRS9YtLzuc6/108/TDzFWpcwZw/nAy4HzgPvJ6UZqZqNXmT71jTkPvdZPP08/xFiVMg3DHOB1fnaR2di048QJrCz4I9mY89Bs3l7KbeiHGKtSJsHtzpL1zWwUyetrPzAuvz/JvFlT6z73Wj/9PP0QY1XK/NCfAHxa0gxJA8PWNrNRJa+v/ecPm8FR++684QxhQOKofXce1Cup1/rp5+mHGKtS+OmqkqYA3wP2y5seEV1pLJzHYGZWXrvGY1hAkvX8fnzz2WzUOPLsG7juzoc3fJ692yR2nbxVodwEaH/f/1MWLh207iIjtZWNZ7SMndAJZc4YngL2iYjfdDakcnzGYDZyjY3CUPIuEbV7HIFTFi7dkAeRNY6Nj0sYah1F4+m3sRM6oV3jMfwOeEF7QjKzXlC0UYDBuQnQ/r7/eeuA+kZhqHUUjWc0jZ3QCWUahlOAL0j6S0nbS5qUfXUqQDPrDXnjMbS773+ZMR/y1lE0ntE0dkInlGkYLgP2AX5K8ijs1enrwfRfMxvF8sZjaPc4AmXGfMhbR9F4RtPYCZ1QpmE4IPN6feZV+2xmfWb2bsVP9htzE6D9ff/z1gGDf6iaraNoPKNp7IROKNwrKSKu7mQgZla9C9+zX0u9kmo3W9vVQ6e2jpH2SioaTytxt3ube1GZXknHA2si4oKG8qOAF0TE1zoQ37DcK8nMrLx29Ur6AJDXZWAZ8MERxGVmZj2oTILbTsDdOeUr0mlm1iVVJVyN9sQuS5RpGO4DXklyhpD1KpKeSWbWBVUNMOOBbMaOMpeSvgN8SdJBksanrzcAXwQu7Ex4ZjacqhKuxkJilyXKnDF8HNgVuByo/XWMA74P/Gub4zKzgqpKuBoLiV2WKHzGEBHrImIe8H+AdwFHAntExBERsa5TAZrZ0KpKuBoLiV2WKD3wTkTcERHfj4j/iog7OhGUmRVXVcLVWEjsssSQDYOkUyRtWWRBkmZL+r/tCcvMiqpqgBkPZDN2DJngJulc4P8CFwOXAosj4r502ubAnsD+wFHAi4CjI+LaDsdcxwluZmbljXignog4RtJ04Hjg28ALJAWwDtgUEHATcBZwbkQ829bIzayQKganaXXdVem1ePrRsL2SImIp8PeS3gu8AtgFmECSu3BLRDiHwayLiuYXdCIPoddyG3otnn5VplfS8xFxS0T8MCK+GxH/60bBrPuqGJym1XVXpdfi6VeleyWZWW+pYnCaVtddlV6Lp1+5YTDrc1UMTtPquqvSa/H0KzcMZn2uisFpWl13VXotnn5V5pEYZtaDqhicptV1V6XX4ulXhQfq6VXOYzAzK2/EeQw5CzocOBDYjobLUBHxtmHmPQd4K/BAROyVM30O8EPgrrTokoj4ZJn4zPpVq33vZ51+Bfc/vjGNaPutN+Xkg/cctEwodjR9ysKluUN75sVZdJnWP8oM7TmfZBS3K4FVQN2MEXHsMPO/FngCOH+IhuFDEfHWQgGlfMZg/a6x7z0k18WLPm6isVFoZvw4gWDd+o1f3bz1nLJwKRfceM+g+WfvNomb7nm0Ls6iy7Te064zhr8G5kXERSMJIiKukTRtJPOajWZD9b0v8uNapFEAWPf84IPAvPUsWJQ3gi9cd+fDI16m9ZcyvZLGAbd0KpDUfpJulfQ/kl7erJKk4yQtlrR49erVHQ7JrLO63fe+cT3r23Df0XkD/a1Mw3AWycPyOuUmYJeImAF8GVjYrGJEnBURMyNi5uTJkzsYklnndbvvfeN6BqS2L9P6y3CP3f5S7QVsA5wg6TpJX89OS6e3JCIei4gn0veXAeMlbdvqcs16Xat977ffetNC9caPE+MH6n/089Yzb9bU3Pln7zZpUJxFl2n9ZbgzhumZ18tJLiU9C7y0Ydr0VgOR9GIpOVSRtE8a20OtLtes17U6zsGijx40qHHYfutN+eLhr6xb5vzDZjD/0BnDrue0udM5at+dN5w5DEgcte/OXPie/QbFWXSZ1l8qy2OQtACYA2wL3E8yhvR4gIg4U9LxwHuB54C1wIkRcf1wy3WvJDOz8trSKynNQzghIh5vKN8S+HJE/M1Q86fjRQ81/SvAV4rGY9ZrqhoHoFmOQZF4Ft/98KB5Z+4yqefzEDzGQrXK5DGsB3aIiAcayrcF7ouIrjxew2cM1gtazUUoqlmOwVH77lzXOOTFMzBOrM/pXtpY3mt5CFXt27FmqDOGYXslSZok6UUko7W9MP1ce00myWa+v70hm/WXqsYBaJZj0FieF09eo5BX3mvjF3iMheoVOcp/kCTLOYDf5UwPkvsFZmNWVbkIzXIMGstbXW8v5SF0O89jLCrSMBxAcrbwc+AdQDb98Vng7ohY1YHYzPrGjhMnsDLnh6rd/fkHpNzGoTH3oFk8RfVSHkJV+9Y2GvZSUkRcHRFXAbsCC9PPtdcNbhTMqhsHoFmOQWN5XjwD4/IT1xrLey0PwWMsVG/IM4b0wXdZu6hJVmREXNOuoMz6TVXjANRuMA/XK6lZPP3YK8ljLFRvyF5Jkp4nuYdQaw1qlRs/ExH1TXpF3CvJzKy8VvIYsg8imgV8DjgduCEt2w/4F+DDrQZp1gva3V/+yLNvqHsq6ezdJrHr5K0GHbVD/llAXs5C3hE+DD6izivzUbYVUSaPYQnwkYi4oqH8IODfI2LvDsQ3LJ8xWLu0u798Y6NQ1u7bbckfH3hyUPk4QbaH6fgBQdQ/AtvjJNhwWspjyNgTWJFTvpLk2Ulmfa3d/eVbaRSA3EYB6hsFSH78G8dFWPd81DUK4L7/VlyZhuG3wMclbegjlr7/WDrNrK+Nhf7yo2lbrHPKPMbivcCPgZWSfp2WTQfWA29pd2BmVRsL/eVH07ZY5xQ+Y4iIX5HkMnyEZFCdm9P3u6bTzPpau/vLz95tUkvx7L7dlrnljekI4weU3FPIlnmcBGtBmUtJRMRT6ehpJ0bEByPi7IjIvxBq1mdaHReh0YXv2W9Q4zB7t0m5Yx3klV1x4pzc8i+8s2GchUNnMP+wGR4nwdpmuDyGQ4AfRcS69H1TEXFJu4Mrwr2SzMzKayWP4SLgxcAD6ftmAuhKgptZr8jLgYBiuQRl8idaybXwuAZWRGUjuHWKzxisF+TlQBTNJSiTP9FKroXHNbCsVsdj2Kz9IZmNLnk5EEVzCcrkT7SSa+FxDayoIt1VH5V0A8ljt68EboyI5zoblll/KZMf0Fi3TP5EK7kWYyFPw9qjSK+kfyLJbv574BrgEUk/kfTPkl4tqVTPJrPRqEx+QGPdZvPmlZep2855bWwpMh7D2RFxVETsBLwMOAlYA3wAuBF4WNIPOxumWW/Ly4EomktQJn+ilVwLj2tgRZXJfCYibgduB86UtAPwPuD9JOM+m41ZzcYMyCtrvNFbZryBVsYm8LgGVlSZp6tuC8whGerzAOAlwBKSy0tXRcTlHYpxSO6VZGZWXit5DEg6g6Qh2J2kIbgaOAG4LiKeamegZtAffe1byVkw63XDnjGko7jdTTJIz2URcVcVgRXlM4bRpR/62ufmLOSMidBrcZtltToew2uBbwKHAL+VdLek8yQdK2nXdgZq1g997XNzFnLGROi1uM2KKtIr6dqIOC0iDgQmAkcDd6X/1hqKczsbpo0V/dDXvpWcBbN+UPbpqs9GxFXAp4CPA18maSze3f7QbCzqh772reQsmPWDQg2DpE0kzZZ0iqSfkeQx/Bw4jOThekd3MEYbQ/qhr31uzkLOmAi9FrdZUUV6Jf0UeA2wBUkG9JXAPwI/j4i7OxuejTX90Ne+lZwFs35QpFfSAtLnJEXEHZVEVYJ7JZmZlddSr6SImJc+FqOlRkHSOZIekPSbJtMl6UuS7pD0a0mvamV9ZmY2MqUeidGic4GvAOc3mf5mkiS63YFZwNfTf80KO2XhUhYsWs76CAYk5s2aymlzp4+4HrR/YBzwJSfrbZU1DBFxjaRpQ1R5O3B+JNe2bpQ0UdIOEXFvJQFa3ztl4VIuuPGeDZ/XR2z4nP3RL1oPBiezrVyzlpMvWQpQemCclWvWctL3b60bvKfM8syq0kuPzJ4CLM98XpGWmRWyYNHyQuVF60H7B8YpOniPWTf1UsOgnLLcO+OSjpO0WNLi1atXdzgs6xfrm3SkaCwvWg86MzBOq3XNOq2XGoYVwNTM552AVXkVI+KsiJgZETMnT55cSXDW+waUd2wxuLxoPejMwDit1jXrtCEbBkmPS3qsyKsNsVwK/HXaO2lf4FHfX7Ay5s2aWqi8aD1o/8A4RQfvMeum4W4+H9+uFaX5EHOAbSWtIHmkxniAiDgTuAw4GLgDeAo4tl3rtrGhduN4uN5GRetBZwbGGenyzKpSeKCeXuUENzOz8lp97LaZmY0hhRsGSZtK+oSkP0h6WtL67KuTQZqZWXXKnDH8G8lTVD8PPA+cBHwVeAh4X/tDMzOzbijTMLwT+IeI+E9gPfDDiHg/yU3kgzoRnJmZVa9Mw7A98Lv0/RMkA/QA/AR4QzuDMjOz7inTMNwD7Ji+vwN4Y/p+P8Bpm2Zmo0SZhuEHwIHp+zOAT0i6i+Spqd9oc1xmZtYlhZ+uGhEnZ9Oj2DMAABBBSURBVN5fJGk5MBv4Q0T8uBPBmZlZ9Qo3DJJeC1wfEc8BRMQiYFE6HvRrI+KaTgVpZmbVKXMp6UpgUk75Nuk0MzMbBco0DCL/MdgvAp5sTzhmZtZtw15KknRp+jaACyQ9k5k8AOwFXN+B2MzMrAuK3GN4KP1XwCPUd019FrgWOLvNcZmZWZcM2zBExLEAkpYBn4sIXzYyMxvFCt9jiIhPRMSTkmZKOlzSlgCStpRUuHeTmZn1tjLdVbcnGWXt1ST3G3YH/gR8AXgaOKETAZqZWbXK9Er6D+A+kl5IT2XKv4+flWRmNmqUuQR0IHBgRDyi+kHT7wR2bmtUZmbWNWXOGCaQ9EJqNJnkUpKZmY0CZRqGa4BjMp9D0gDwz8DP2hmUmZl1T5lLSR8Grpb0amAzkpHcXk7ySIzZHYjNzMy6oEx31d8BrwBuAH4KbE5y43nviLizM+GZmVnVSuUfRMS9wMc6FIuZmfWAYc8YJG0h6auSVkp6QNJ3JG1bRXBmZla9ImcMnyC56XwhSe+jecDXgcM6F9botvDmlcy//HZWrVnLjhMncNIb92Du3lO6HZaZGVCsYTgE+NuI+C6ApAuA6yQNRMT6jkY3Ci28eSUnX7KUteuSXbdyzVpOvmQpgBsHM+sJRW4+TwV+UfsQEb8EngN27FRQo9n8y2/f0CjUrF23nvmX396liMzM6hVpGAYYnNj2HCVvXFti1Zq1pcrNzKpW5MddDB6gZ3PgbEkbnpkUEW9rd3Cj0Y4TJ7AypxHYceKELkRjZjZYkTOG84BVJAP21F4XAMsbyqyAk964BxPGD9SVTRg/wElv3KNLEZmZ1Ss8UI+1R+0Gs3slmVmv8n2CLpi79xQ3BGbWs8o8RK9lkt4k6XZJd0j6SM70YyStlnRL+vq7KuPrpoU3r2T2Z37Orh/5b2Z/5ucsvHllt0MyszGqsjOG9EmsXwUOAlYAv5J0afoMpqzvRcTxVcXVC5zbYGa9pMozhn2AOyLiTxHxLPBd4O0Vrr9nObfBzHpJlQ3DFJKeTDUr0rJG75D0a0kXSZqatyBJx0laLGnx6tWrOxFrpZzbYGa9pMqGQTll0fD5R8C0iHgF8L8kXWUHzxRxVkTMjIiZkydPbnOY1WuWw+DcBjPrhiobhhUkj9eo2YkkP2KDiHgoImqJdGcDf15RbF3l3AYz6yVVNgy/AnaXtKukTYEjgEuzFSTtkPn4NuC2CuPrmrl7T+HTh0xnysQJCJgycQKfPmS6bzybWVdU1ispIp6TdDxwOcnzl86JiN9K+iSwOCIuBd4v6W0kz2J6mPoxpkc15zaYWa9QRONl/v4yc+bMWLx4caXrLDOewpFn38B1dz684fPs3SZx2MydB80PxbKhPZaDmbWDpCURMTN3mhuGchpzDiC5H5B36aexUagR9Xfdxw8IAtY9v7E0b5ll1m1mNpShGoZKM59HgzI5B3mNAgzuirVufdQ1Cs2W6XwHM6uCG4aSqsw5aFym8x3MrApuGEqqMuegcZnOdzCzKrhhKKlMzsHs3SblLqMx02/8gBg/rr40b5nOdzCzKrhhKKlMzsGF79lvUOMwe7dJ/Mfhr6ybf/6hM5h/2Ixhl+l8BzOrgnslmZmNQe6VZGZmhXkEtxE4ZeFSFixazvoIBiTmzZrKXaufGJTIduF79iuckObENTPrFb6UVNIpC5dywY33FKq7+3ZbsuKRp4dNSHPimplVzZeS2mjBouXDV0r98YEnCyWkOXHNzHqJG4aS1rfhDMuJa2bWy9wwlDSgvPGGynHimpn1MjcMJc2blTvaaK7dt9uyUEKaE9fMrJe4YSjptLnTOWrfnTecOQxIHLXvzrmJbFecOKdQQpoT18ysl7hXkpnZGDRUr6Qxm8dQNG8gL2dh0Z8e4o8PPLmhzu7bbcldq5/kuUwbu4ngjk+/hZd+9DKeXr9xwuYDYpstxnP/489uKNt+6005+eA9ne9gZj1hTJ4xFM0bKJOz0G7OdzCzTnIeQ4OieQNlchbazfkOZtYtY7JhKJo30I6chVY438HMumFMNgxF8wbakbPQCuc7mFk3jMmGoWjeQJmchXZzvoOZdcuYbBiK5g00y1nYfbst6+rtvt2WbNJwcrGJYNln3sLmA/UTNh8Q22+9aV3Z9ltvyhcbBu9xvoOZdcuY7JVkZjbWOY+hoFZyBPLyHU6bO71puZlZr3LDkGrMEVi5Zi0nX7IUYNjGoTHfYX0EF9x4z6BEuFo54MbBzHrWmLzHkKeVHIFm+Q7ZRqFIfTOzXuCGIdVKjkDZfIdu50eYmQ3FDUOqlRyBsvkO3c6PMDMbihuGVCs5As3yHRq7tQ5X38ysF7hhSLWSI9As3+GKE+fklvvGs5n1MucxmJmNQT3zdFVJb5J0u6Q7JH0kZ/pmkr6XTl8kaVqV8ZmZWYUNg6QB4KvAm4E9gXmS9myo9rfAIxHxZ8B/AJ+tKj4zM0tUecawD3BHRPwpIp4Fvgu8vaHO24Hz0vcXAQdK7sJjZlalKhuGKUA2s2tFWpZbJyKeAx4FXtS4IEnHSVosafHq1as7FK6Z2dhUZcOQd+TfeOe7SB0i4qyImBkRMydPntyW4MzMLFFlw7ACyHbg3wlY1ayOpE2AbYCHK4nOzMyAah+i9ytgd0m7AiuBI4B3NdS5FDgauAE4FPh5DNOfdsmSJQ9KuruFuLYFHmxh/l4ymrYFRtf2jKZtgdG1PWN1W3ZpNqGyhiEinpN0PHA5MACcExG/lfRJYHFEXAp8E/i2pDtIzhSOKLDclq4lSVrcrC9vvxlN2wKja3tG07bA6Noeb8tglT52OyIuAy5rKPtY5v3TwGFVxmRmZvX8SAwzM6vjhgHO6nYAbTSatgVG1/aMpm2B0bU93pYGff+sJDMzay+fMZiZWR03DGZmVmfMNgySzpH0gKTfdDuWVkmaKulKSbdJ+q2kE7od00hJ2lzSLyXdmm7LJ7odU6skDUi6WdKPux1LqyQtk7RU0i2S+v5595ImSrpI0u/T789+3Y5pJCTtkf6f1F6PSfrAiJc3Vu8xSHot8ARwfkTs1e14WiFpB2CHiLhJ0tbAEmBuRPyuy6GVlj40ccuIeELSeOBa4ISIuLHLoY2YpBOBmcALIuKt3Y6nFZKWATMjYlQkhEk6D/hFRHxD0qbAFhGxpttxtSJ9kvVKYFZEjCj5d8yeMUTENYySx21ExL0RcVP6/nHgNgY/oLAvROKJ9OP49NW3Ry+SdgLeAnyj27FYPUkvAF5LklhLRDzb741C6kDgzpE2CjCGG4bRKh3caG9gUXcjGbn00sstwAPAFRHRt9sCfBH4MPB8twNpkwB+KmmJpOO6HUyLXgKsBr6VXur7hqT8gdr7yxHAglYW4IZhFJG0FXAx8IGIeKzb8YxURKyPiFeSPGhxH0l9ealP0luBByJiSbdjaaPZEfEqkgG3/jG9JNuvNgFeBXw9IvYGngQGjSzZT9LLYW8Dvt/KctwwjBLp9fiLgQsj4pJux9MO6Wn9VcCbuhzKSM0G3pZel/8u8HpJF3Q3pNZExKr03weAH5AMwNWvVgArMmekF5E0FP3szcBNEXF/KwtxwzAKpDdsvwncFhFf6HY8rZA0WdLE9P0E4C+B33c3qpGJiJMjYqeImEZyev/ziDiqy2GNmKQt084NpJdc3gD0ba++iLgPWC5pj7ToQKDvOmw0mEeLl5Gg4ofo9RJJC4A5wLaSVgAfj4hvdjeqEZsNvBtYml6bB/iX9KGF/WYH4Ly0Z8U44L8iou+7eY4S2wM/SEfb3QT4TkT8pLshteyfgAvTSzB/Ao7tcjwjJmkL4CDg71te1ljtrmpmZvl8KcnMzOq4YTAzszpuGMzMrI4bBjMzq+OGwczM6rhhsFFP0jGSnhi+Zm9qJX5Jr5P0h7T7b0dImi5p5Sh5nIThhsEqIulcSZG+1kn6k6TPlfkxSZfRkZyG9HHSH+rEsrscx3zg9IhY38Zl1omIpcCNwImdWodVyw2DVel/SRLYXgKcArwP+FxXIxrFJL0GeCktPjenoG8B75U0ZpNmRxM3DFalZyLivohYHhHfAS4E5tYmStpT0n9LejwdRGmBpBen004FjgbekjnzmJNO+4yk2yWtTY+4/13S5u0MfKjY0unnSvqxpBPSyyqPSPpWmo1aq7OlpPMlPSHpfkknp/Ocm06/CtgFmF/bxoYYDpT0G0lPKhmYaddhwn4X8L8R8VTDct4iaVG6vx6S9KPa/kr338fS7Xlc0nJJhysZ0Oa7aex/lPSGhnX9FJhE8jQB63NuGKyb1pKMt1AbbOgakmfv7EPyjKStgEsljSM5s/gvNp517ABcny7nSeBvgJeRnIUcAXy0XUEWiK3mL4C90umHA38FZEfT+zzwurT89cCMdJ6aQ0ge7PbJzDbWbAacTLKd+wETgTOHCf0vgLpR1iS9CfghcAXw58ABwNXU/xZ8APglyQPl/gs4D/gOcBnwynRfXJBtfCPiWeCWdPus30WEX351/AWcC/w483kf4EHge+nnTwI/a5jnhSTP/98nbxlDrOsfgDsyn48BnhhmnmXAh5pMKxrbcmCTTJ2zSY7YIWlIngWOyEzfEngEOHeoONL4A9gjU3ZkurxxQ2zTGuDYhrLrgO8Osx8WZD5vla77S5myaWnZzIZ5LwG+3e2/Nb9af/l6oFXpTWnvmk1IzhR+SPIQM0iOXl/bpPfNbiRHsLkkHUpylPtnJD9kA+mrXYrG9ruIeC4zbRUwK1NvfKYuEfGkio85/kxE3N6w7PEkZw7NRiKcADzdULY3SSM2lF9nYnxC0lPA0sz02iOdt2uYb226TutzbhisStcAxwHrgFURsS4zbRzw30Bej5ymz5aXtC/JWAefAD5IcpT8Ntp7U7tobOsapgUbL9EoUzYSzzV8ri1nqMvBD5Kc2ZSVtx3rGj7nrXsSyRmH9Tk3DFalpyLijibTbgLeCdzd0GBkPcvgM4HZwMqI+LdagaRdWo60fGzDuYPkx3Uf4C7Y8JjkvYA7M/XytnGkbgb2zCk7kOQyV7vtRXI5yfqcbz5br/gqsA3wPUmzJL1E0l9KOqs2OAzJ0ehekvaQtG06at0fgCmSjkzneS/JYCUjsaOkVza8ti0Y25Ai4gngHOCzae+iPYFvkHwHs2cRy4C/kDQlXXcrLgf2byg7HThM0mlpT6uXS/pgtvfUSCgZa3wKSe8k63NuGKwnRDJk5GzgeeAnwG9JfpCfSV+QHOXeRtLTZjXJ+MM/Ikni+iLJtfGDgI+NMIwPkhxRZ19HFIytiA8BvwAuBa5M411M/X2AjwFTSc4iVo9wO2ouAP6PpJfXCiIZvOmvSIaAvJmkR9IBJNvWinnATyPi7haXYz3AA/WYdYmkzYC7gfkR8fkOreMzwOSI+NtOLD9dx2bAH4F5EXFdp9Zj1fEZg1lFJO0t6V2S/kzS3iT5AVsD3+vgaj8F/EkdfFYSSVLe6W4URg+fMZhVJG0Mzgb2IOlldAtJzsKSrgZm1sANg5mZ1fGlJDMzq+OGwczM6rhhMDOzOm4YzMysjhsGMzOr8/8B0ksJLM4NI0AAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Let's start looking at some data - plotting petal length vs. petal width\n", "plt.scatter(iris.data[:, pl_ind], iris.data[:, pw_ind])\n", "\n", "# Add title and labels\n", "plt.title('Iris Data: Petal Length vs. Width', fontsize=16, fontweight='bold')\n", "plt.xlabel('Petal Length (cm)', fontsize=14);\n", "plt.ylabel('Petal Width (cm)', fontsize=14);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just from plotting the data, we can see that there seems to be some kind of structure in the data. \n", "\n", "In this case, we do know that there are different species in our dataset, which will be useful information for comparing to our clustering analysis. \n", "\n", "Note that we are not going to use these labels in the clustering analysis itself. Clustering, as we will apply it here, is an unsupervised method, meaning we are not going to use any labels to try and learn the structure of the data. \n", "\n", "To see the structure that is present in the data, let's replot the data, color coding by species." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEcCAYAAADDfRPAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgU1bn48e87CzKsiqAgKGBUUARFR0XxKsa4RL2iESLGDdFr1KgYf66J0VEx8YZcExKNXk0UF1QWFRe8alBxxWWG1Q1EHAPDNqAiOzPM+/ujqqG6p3q6erq7epn38zz90H3qVNVbzXSfrlPnrSOqijHGGBNRlO0AjDHG5BZrGIwxxkSxhsEYY0wUaxiMMcZEsYbBGGNMFGsYjDHGRLGGwYeIVIiIuo+KJNYb71mvV5piqfZss0FENovIchF5T0RuE5HdU9z+Ge7xVojIzumIOeB+h3iOy/v4QUTeF5GLUtx25Jh6pSnGigD1tx9Hc/cZJhEZGXmffJYldey5RkRO8sT/vzHLHvYsGxOz7AXPsgPcsshnsDrAfr3fHSM95XE/ZyLSy7PO+BQOO21Ksh2ASYoAOwFd3cdRwGgROVtV/9XMbZ4BXOg+Hw98n2qQKWoPHAkcKSIDVPXXzdjGEOA29/kMoDotkRWekcCx7vOK7IWRETOBBpwfv4Njlh3leR5v2XfA52mMJ9c+Z02yM4Y0EJHWAKo6UlXFfVRnYFe9cRqG/jh/XAC7AM+JyP4Z2F9Y3lJVAVrjfFlFjE7XmZdpWVT1B+AT9+UBIrILgIh0Bvp4qh4uIiXusr7Arm75++pm/6pqL/cz3SuU4HOANQxJiOkq+g8RmSIia3F/WcTrShKR/xKRShH5VkS2iEiNiPxLRC6Ms6u4VHWrqn6iqhcBz7vFbdnxCxkROdvd/hIR2ejuc7GIPODtenK7PLwxfO2NX0TaicijIjJfRNaISJ2IfC8ib4vI2T7vT2Td6mSPyz22Lar6KDA/skmg3LP9QSLynIisdGNZ5r7nvTx1qr3vBfCmJ64hbp37RWSWiNS621knIh+JyOUiIs2JPVki8gv3fVzr/v8sFJG7RKRNTL1I7DNE5BQR+VhENonIVyJyQ2y8bp254nQ5LhCR82P+LodEui7YcbaQsBtMRK50Y9woInNE5KcJjq+LiGx1t/lKzLKfevZ3h1t2kIg86342trh/b7NE5H9FpDTJtzfivcgucc5CYccZQSVQC7QBBrpl3rOH9z3x+nYliciRIjLTfa+rReRavyASfc586v9cROa5/8+fN+d7ImWqao+YB85ptbqPCk/5eE/5as/zap/lvdyy4Z6y2MeUALFUx27Ts+xIz7IfgCK3/IEm9vkF0CryY6iJRy+c7qqm6lwYE0/U+5HguIZ46s+IWfaJZ9kwt+znQH2cONYAfXzer9jHELfO5ibq3BYnxooAx7R9Ownq/a2J/X8MlPlscy1O10hs/fNi4vV7j2q874H7fxv3/9Xn2Ff41NsK7J3gOJ9169YBXTzlj7nlDThnwW1wvqTjxdSumZ/jcz3buMstu9t9/RfgOff5aHfZw576x/p8Bqs9ZQcAG3xiXeZ5PjLg58z7/+H3XitwdJjfgXbG0Hw/4HwxlwGnNFHvGPff9TinsDsBPXG+6F6Jt1JAX3iet2fHafCTwBFAZ6AU2B14xF3WJxKvOt03j3q20Vuju8LWAWfj/OG2wenqOQrY6NZvTv9/XCLSyv111M8tUuBj91f0/UAxMAvoi/M+HofzBdUJGOseUy/gds9mj/Mc0wy3bBSwL8571goYACx1l43O5FmDiAwCrnRfjsdpfNsA17tl5cDlPqt2AP6A03V4paf8fM/zu3DeI3DOmjoC5wB7eDekqtXu//1bnrLIe+R37J2AnwE7AxPcslKcv42mPOz+W4Lz946IlOH0t4Pzg+BrYH+cv1WAG3D+zroAR7vHXJ9gP/G853l+tPvvYM+y92LKInXqcBropvwO5/8N4O84782Pcd7zKAE+Z167A1e42/tvT/n5hCnMVihfHgQ7Y/iFz3re5b3csmvZ8evoMWA0cCLQMWAs1bHb9CzblehfFV3c8n2Bx4F/43xxxv76uLGpmD3LBOeP9AP8f7FuSuE9HuITV+xjnFv3hAB1N3m27f3/G+Kz77NxLkp/C2zz2dbuPjFWBDimqF/dcercFeBYXvHZ5gqg2C1r5yn/wi1r4zmWNZG67rL3/N4P9z3wjTfm2Kd4yk/zlD+Q4P0oZscv6HfdMu8Z9Hlu2W44X8YKVAG34jQk+6Ths7zU3e5G933b5L7eAxjkPq/BaYgicX0U5zNY7Slb6anf0VP+uKd8ZMDPWS/PskpP+YF+fxNhPOyMoflmB6z3d2Ayzpfq+TinsK8CK0XkphRj6Ot5/gOwRkQ6Au8C5wF74vyyi1UWcPs3AvfhnH10wGkovFonFW0wG4APgUuBa9yy3QKs11pE2iaqJCIjgKdx+td3wf86W9D3pzmCHMuuPmVfqeo29/kGT3nk/6ATO46lxlMXnB8IqVjgee63b19uDI+5L48SkZ7AL9zXa4Fn3HqrgMtwumcPwTnjmwh8KSLviEiHFGKPXCsoAy5xY65W1WU4jdAmnEbiXM8675FY5P9onaqu9ZQv9auchGa91+lmDUPzbQpSSVU3q+rPcT64R+N0Y3yI0xXyexHpnkIMN3ueT1PVBpzulciXz+tAN3VOZa+OF2IT2x/heX4GsJO7rTXNjDeet3THqXU7VR2kqg+p+7MJWOWp+5Cnrrf7o0hVIx+koMd0FU5/vuB0UYXBeyznxjmWw33Wq4s88bwvXt/i/PgA6CYi3s/2nnFiaep98t13EutEPOL+Kzhf/pGL1k+p6vbPkKr+E6dbrT9wFvBXd9HRwK+S3KeX90s+0vX5vrtPb5fRr+OsE89q99/27o+xiB5x6ofxXqeNNQwZJiJniciVQHdgLs7Zw9zIYuL/IcXbXqmIHChOIsypbvEG4A73ubc/djOwQUT64XwJ+vF+yR8U07/u3db3QKmI/A7/X7Qpj0pqwvs448oBLnRH9LR3R74MFpGxOGdiEd5jGhDzJek9ph+csOUidoxMSQsROdnnsRvwkqfaGDf+1iKypzta50mif70Goqobcbr8wOmvv859j0YQPW7fa/v7JCIHJ7vPgHEtYMev9utwfhDBjusPiEhnEfkTzjW7VcCL7BhxB7CXp26yf2PeL/m9fMrejd2HJ96mvOl5/nsR6Sgix+Fci/HT1Ocs94TZb5UvD4JdY+jls16j5cAtnrLYxzI8I1DixFLdxPqK80vxBE/9XXA+XLH1FsY5pmE+davdZb/1WVaL8yXdqG86dv0ExzXEU39GgPrn4H89IPIY76lb7lfHXXauz7KNwBKf/ztvjBUBYmzq/0mBM9x6f09Qb6TPNmfE2Vd1zHvqNyrJO1LmWE/963zqzmjq2GPKxyd6T9x1Lo7Zx/yY5T2CvG/J/o259UtwBn54t3eQZ/kpMcu+buIz6H2v441K8o6u8v4/NvU56+X3nsaUJ/yMpPNhZwyZ9zrOKKFFOH+g24DluP3c6jmdDkCBLTgXIt/HGXmyv3qynlX1O5zT9XdxvvCW4TR0d8fZ5jM4F0T/7cbm9d/A73Euzm3CGcXyY5z+4VCp6lM43QrP4Fz4q8f5EFa6cf6Pp24lTtfZV0SfmqOqE3C6Db7GOaOqxHm/vsr4QeyI4Qqca0Bv4byXdTh902/ijMr5v2ZudwYwFCcPZCvOj4ELie4m8/5yvQ9naPNyyGi3xUSi+8sfjln+HfBnnG6d1Th/h+tw/sZ/oapTm7tjVa0HPvIUrWNHngzuPhpiXgfZ7mc4gyI+xHmvl+D8kPp7nFWa+pzlHHFbJmNMnnMTwY4D3lSn/xwRORmYitOFsxzooc61KGPisnslGVM4dsIZ8VYnIitx8jQiF0brgcusUTBBWFeSMYVjC04i1Tc4o+DKcLoungAOU9UXshibySPWlWSMMSaKnTEYY4yJkvfXGDp37qy9evXKdhjGGJNXqqqqVqtqF79led8w9OrVi8rKymyHYYwxeUVEvom3zLqSjDHGRLGGwRhjTBRrGIwxxkTJ+2sMfurq6li6dCmbN2/OdigFoXXr1vTo0YPS0ubOsGiMyScF2TAsXbqU9u3b06tXL3L9Joa5TlVZs2YNS5cupXfv3tkOxxgTgoLsStq8eTO77rqrNQppICLsuuuudvZlTAtSkA0DYI1CGtl7aUzLElpXkojsiTPNX1ec29w+qKrjYuoMwZmg42u36FlVvQNjjPGYtnga42aNY8WGFXRt25XRh4zm1L1PTbxiE+sDKW2zkIR5jaEe+H+qOktE2gNVIvIv977mXu+o6mkhxpV148eP58QTT2SPPfbIdijG5Lxpi6dR8X4Fm7c53ZvLNyyn4v0KgEBf5H7r3/LuLYgIdQ11zdpmoQmtK0lVl6vqLPf5OuBznOkuW7zx48ezbNmybIdhTF4YN2vc9i/1iM3bNjNu1rg4ayRev17rtzcKzdlmocnKNQYR6YUzx+6HPouPFJG5IvJ/7lzFfutfKiKVIlJZW1ubcjxTZ9cw+O436H3TNAbf/QZTZ9ekvM0NGzZw6qmnctBBB3HggQcyceJEqqqqOPbYYzn00EM56aSTWL58OVOmTKGyspJzzz2Xgw8+mE2bNvH6668zcOBA+vfvz6hRo9iyZQsAN910EwcccAADBgzguuuuA+DFF1/kiCOOYODAgfzkJz9h5cqVKcduTC5bsWFFUuXNrZds3UISesMgIu1wprm7RlV/iFk8C+ipqgcBf8OZeaoRVX1QVctVtbxLF997QAU2dXYNNz87n5rvN6FAzfebuPnZ+Sk3Dq+88gp77LEHc+fO5ZNPPuHkk0/mqquuYsqUKVRVVTFq1Ch++9vfMmzYMMrLy5kwYQJz5sxBRBg5ciQTJ05k/vz51NfXc//99/Ptt9/y3HPP8emnnzJv3jxuueUWAI4++mg++OADZs+ezYgRI/jjH/+YUtzG5LqubbsmVd7cesnWLSShNgzu1IPPABNU9dnY5ar6g6qud5+/DJSKSOdMxjT21QVsqouegnVT3TbGvrogpe3279+f6dOnc+ONN/LOO++wZMkSPvnkE0444QQOPvhgxowZw9KlSxutt2DBAnr37s1+++0HwIUXXsjbb79Nhw4daN26NZdccgnPPvssbdq0AZycjZNOOon+/fszduxYPv3005TiNibXjT5kNK2LW0eVtS5uvf0CcnPWL5ESSouiEziT2WahCXNUkgD/BD5X1Xvi1OkKrFRVFZHDcRquNX5102XZ95uSKg9qv/32o6qqipdffpmbb76ZE044gX79+jFz5swm14s3cVJJSQkfffQRr7/+Ok8//TT33nsvb7zxBldddRXXXnstp59+OjNmzKCioiKluI3JdZGLwc0dQRRv/VS2WWjCHJU0GDgfmC8ic9yy3wB7AajqA8Aw4HIRqQc2ASM0w1PM7bFzGTU+jcAeO5eltN1ly5bRqVMnzjvvPNq1a8eDDz5IbW0tM2fO5Mgjj6Suro6FCxfSr18/2rdvz7p16wDo27cv1dXVLFq0iH322YfHH3+cY489lvXr17Nx40ZOOeUUBg0axD777APA2rVr6d7duYb/6KOPphSzMfni1L1PTelLO976LbUhiBVaw6Cq7wJNZkqp6r3AveFE5Lj+pD7c/Oz8qO6kstJirj+pT0rbnT9/Ptdffz1FRUWUlpZy//33U1JSwtVXX83atWupr6/nmmuuoV+/fowcOZLLLruMsrIyZs6cySOPPMLw4cOpr6/nsMMO47LLLuPbb79l6NChbN68GVXlz3/+MwAVFRUMHz6c7t27M2jQIL7++usEkRmTPqnmEwQx5oMxTF44mQZtoEiKGL7fcG4ZdEta92Gi5f2cz+Xl5Ro7Uc/nn3/O/vvvH3gbU2fXMPbVBSz7fhN77FzG9Sf14YyBNpLWK9n31BS+2HwAcPrlK46qSFvjMOaDMUxcMLFR+dl9zrbGIUUiUqWq5X7LCvImesk6Y2B3awiMSVJT+QTpahgmL5wct9wahswp2HslGWMyK9V8giAatCGpcpMe1jAYY5ol1XyCIIrE/ysqXrlJD3t3jTHNkmo+QRDD9xueVLlJD7vGYIxpllTzCYKIXEewUUnhslFJJhB7T40pLE2NSrKupDxx6623Mn369KTXmzFjBqed1qLuYm6MSZF1JeUQVUVVKSpq3F7fcUc48xXV19dTUmJ/FiaYoBPeBC0L2g2VTGJd0LphJOulKqwY7RsAYN4keP0OWLsUOvaA42+FAT9v9uZuvPFGevbsyRVXXAE42cnt27enoaGBSZMmsWXLFs4880xuv/12qqur+elPf8pxxx3HzJkzmTp1KrfddhuVlZWICKNGjeLXv/41I0eO5LTTTmPYsGF8/PHHjB49mg0bNrDTTjvx+uuvU1payuWXX05lZSUlJSXcc889HHfccVFxffvtt4waNYrFixfTpk0bHnzwQQYMGEBFRQXLli2jurqazp078+STT6b0dpqWIeiEN79773eoKvVa32S9oBPjJDNRT9C6qU7+E4YwY7SupHmT4MWrYe0SQJ1/X7zaKW+mESNGMHHijmzNSZMm0aVLF7788ks++ugj5syZQ1VVFW+//Tbg3FH1ggsuYPbs2axevZqamho++eQT5s+fz0UXXRS17a1bt3L22Wczbtw45s6dy/Tp0ykrK+O+++4DnFtxPPXUU1x44YVs3hydfHTbbbcxcOBA5s2bx+9//3suuOCC7cuqqqp4/vnnrVEwgQWd8KauoW57o9BUvaAT4yQzUU/QuqlO/hOGMGO0huH1O6Au5iZ6dZuc8mYaOHAgq1atYtmyZcydO5dddtmFefPm8dprrzFw4EAOOeQQvvjiC7788ksAevbsyaBBgwDYe++9Wbx4MVdddRWvvPIKHTp0iNr2ggUL6NatG4cddhgAHTp0oKSkhHfffZfzzz8fcG7E17NnTxYuXBi1rrfOj3/8Y9asWcPatWsBOP300ykrS+3GgaZlycQkNkG2mUxiXdC6YSTrpSrMGK1hWNt4ToQmywMaNmwYU6ZMYeLEiYwYMQJV5eabb2bOnDnMmTOHRYsWcfHFFwPQtm3b7evtsssuzJ07lyFDhnDfffdxySWXRG1XVXHuYE6j8kT86kS25Y3BmCAyMYlNkG0mk1gXtG4YyXqpCjNGaxg69kiuPKARI0bw9NNPM2XKFIYNG8ZJJ53Eww8/zPr16wGoqalh1apVjdZbvXo1DQ0NnHXWWdx5553MmjUrannfvn1ZtmwZH3/8MQDr1q2jvr6eY445hgkTJgCwcOFC/v3vf9OnT/QdYr11ZsyYQefOnRudkRgTVNAJb0qLSimRkoT1gibHJZNYF7RuGMl6qQozRrv4fPytzjUFb3dSaZlTnoJ+/fqxbt06unfvTrdu3ejWrRuff/45Rx55JADt2rXjiSeeoLi4OGq9mpoaLrroIhoanHvB/OEPf4ha3qpVKyZOnMhVV13Fpk2bKCsrY/r06VxxxRVcdtll9O/fn5KSEsaPH89OO+0UtW5FRQUXXXQRAwYMoE2bNjZ/g0lJMhPeBC0LchE1mcS6oHXDSNZLVZgxWoIbpH1UUiGyBDdjCovddjuRAT+3hsCYHFVIeQh+cjFuaxiMMTmrkPIQ/ORq3Hbx2RiTswopD8FPrsZtDYMxJmcVUh6Cn1yN2xoGY0zOKqQ8BD+5Grc1DMaYnFVIeQh+cjVuaxhCsmzZMoYNG5b0epdccgmfffZZk3UeeOABHnvsseaGZkzOOnXvU6k4qoJubbshCN3adqPiqArfPIQg9XJNrsZteQxZli+3uc6n99QYk5jlMSSQ7nHE8W67/cgjj/DJJ58wfvx4pk2bxubNm9mwYQPTp0/nyiuv5K233qJ37940NDQwatQohg0bxpAhQ/jTn/5EeXk57dq1Y/To0bz00kuUlZXx/PPPs/vuu1NRUUG7du247rrrWLRoEZdddhm1tbUUFxczefJkdt99d4YOHcp3331HXV0dY8aMYejQoel6+4xplrDmcsilPIFciqUpLb5hyMQ44hEjRnDNNddsbxgmTZrEAw88wCOPPLK9zsyZM5k3bx6dOnViypQpVFdXM3/+fFatWsX+++/PqFGjGm13w4YNDBo0iLvuuosbbriBhx56iFtuiZ779txzz+Wmm27izDPPZPPmzTQ0NNCqVSuee+45OnTowOrVqxk0aBCnn3667834jAmD3+cu6LwNfvVSnY8hDLkUSyIt/hpDJsYR+912e6+99oqqc8IJJ9CpUyfAuR328OHDKSoqomvXro0m2Ilo1arV9mk6Dz30UKqrq6OWr1u3jpqaGs4880wAWrduTZs2bVBVfvOb3zBgwAB+8pOfUFNTw8qVK5t9fMakyu9zF3TeBr96qc7HEIZciiWRFn/GkKlxxJHbbq9YsYIRI0Y0Wu69zXXQ6zylpaXbf+UXFxdTXx/94Yi3nQkTJlBbW0tVVRWlpaX06tWr0SQ+xoQprLkccilPIJdiSaTFnzFkahxx7G23m3L00UfzzDPP0NDQwMqVK5kxY0az9tmhQwd69OjB1KlTAdiyZQsbN25k7dq17LbbbpSWlvLmm2/yzTffNGv7xqRLWHM55FKeQC7FkkiLbxgyNY449rbbTTnrrLPo0aMHBx54IL/85S854ogj6NixY7P2+/jjj/PXv/6VAQMGcNRRR7FixQrOPfdcKisrKS8vZ8KECfTt27dZ2zYmXfw+d0HnbfCrl+p8DGHIpVgSseGq5MZIgfXr19OuXTvWrFnD4YcfznvvvUfXrrnzS8KGq5p0s1FJ2Y2lqeGq1jDkiCFDhvD999+zdetWbrjhBkaOHJntkKLk43tqjIkvJ/IYRGRP4DGgK9AAPKiq42LqCDAOOAXYCIxU1Vmx2ypEzb2uYEy6xftVm8ov/FwbjpluqZwJ5NJZRERoZwwi0g3opqqzRKQ9UAWcoaqfeeqcAlyF0zAcAYxT1SOa2m68M4a+ffvaOP00UVW++OILO2NoAWLH2oPTDz50n6E8v+j5qPLSotKofAJwrgl48w4i6+fCbR4yJd57FuSYU1k3VU2dMYR28VlVl0d+/avqOuBzoHtMtaHAY+r4ANjZbVCS0rp1a9asWRN4GKiJT1VZs2YNrVu3TlzZ5L14Y+0nL5zc7LyDXB2rny6p5Cfkam5DUl1JItIb6AWUAbXAfFVNekC8iPQCBgIfxizqDizxvF7qli2PWf9S4FKgUeIYQI8ePVi6dCm1tbXJhmZ8tG7dmh49emQ7DBOCeGPqG7QhI9stBKnkJ+RqbkPChsH9Er8cOAfnS9rbP7NVRN4BHgSeUU381yMi7YBngGtU9YfYxT6rNPrZr6oPuvukvLy80fLS0lJ69+6dKBRjTIyubbuyfMPyRuVFUpRS45CLY/XTJd57FuSYU1k3k5rsShKRccBcYG/gt8ABQEegFc5F5FOAd4E7gXkicliC7ZXiNAoTVPVZnypLgT09r3sAywIdiTEmZfHG2g/fb3iz8w5ydax+uqSSn5CruQ2Jzhi2Aj9S1dU+y1YBb7iP290Lxz2Bj/025I44+ifwuareE2d/LwBXisjTOBef16pq4+bUGJMRkQuefqNkBu420EYl+WjqPcvkupkU5qiko4F3gPk4w1UBfgPsBaCqD7iNx73AyTjDVS9S1UqfzW3nNyrJGGNM03Iij0FV38X/GoK3jgK/CiciY0wy/Mbbz141m8kLJ9OgDRRJEcP3G84tg24JtG62fxXHyocYwxL4jEFEdgEqgOOA3Yi5PqGqu6U7uCDsjMGYzPMbb18sxWzTbY3qnt3n7KjGIZtj9YPKhxjTLS23xBCRF4F+wKPASmJGC6nq/6YYZ7NYw2BM5p045UTf0TN+iqSIuRfMTbhut7bdeG3Ya2mLMRX5EGO6pasraQhwbEu5RYUxZodkxtXHDmvN1bH6XvkQY5iSyXz+Ksn6xpgCkcy4+iKJ/prIh3kI8iHGMCXzRT8a+IOIHCQixZkKyBiTe/zG2xfH+RoYvt/whOvmwlh9r3yIMUzJdCUtwrkVxiyg0Q3qVNUaC2MKVLzx9kFGJeXqWH2vfIgxTMlcfH4b2AV4AP+Lz8+kPboA7OKzMcYkL10Xn8uBw1X1k/SEZYzJtjEvjWTy6koacPqVh3cuh877BMpNgPSP/R/zwZhG+/bLuPbbRzKxFNr8CemWzBlDJXC1qr6f2ZCSY2cMxjTPmJdGMnF1JXi7hVWjX7ticxMg/WP/x3wwhokLJjYqFwT1dFD47SOZWPJ1/oR0S9d8DLcA94jIT0RkdxHp5H2kJ1RjTFgmxzYK4NsoAExeOLlRWbrnEvDbBxDVKMTbRzKxFOL8CemWTFfSy+6/rxF9fUHc13bx2Zg8ksxNtP1uuZ3usf/J3NY7dh/JxFKI8yekWzINw3EZi8IYE7oigjcOsbkJkP65BJKZ8yF2H8nEUojzJ6Rb4K4kVX2rqUcmgzTGpN/wzuXONQWvONccY3MTIP1j//32Ac41hkT7SCaWQpw/Id0CnzGIyJXA96r6REz5eUAHVf17uoMzxmTOLaeNhxRGJaV77H9kH80ZlZRMLIU4f0K6JTMqaRFwcezZgTvPwiOqum8G4kvIRiUZY0zy0jUqqQfwjU/5UneZMSYLpi2exolTTmTAowM4ccqJTFs8LTM7mjcJ/nwgVOzs/DtvUmb2Y7IumYvPK4CDgeqY8kMAv6k/jTEZFjuufvmG5VS8XwGQ3u6NeZPgxauhbpPzeu0S5zXAgJ+nbz8mJyRzxvAk8FcROUFESt3HicBfgAmZCc8Y05TQxtW/fseORiGibpNTbgpOMmcMtwG9gVeByLRNRcBk4HdpjssYE0Bo4+rXLk2u3OS1ZIar1qnqOcB+wC+Ac4E+qjpCVesyFaAxJr7Q5hHoGOcyYrxyk9eSnnhHVRep6mRVnaSqizIRlDEmmNDG1R9/K5SWRZeVljnlpuA02TCIyC0i0jbIhkRksIj8Z3rCMsYEcerep1JxVAXd2nZDELq17ZaZG7oN+Dn851+h456AOP/+51/tws5ZGGkAABtKSURBVHOBajKPQUTGA/8JPAO8AFSq6gp3WWvgAOBo4DxgV+BCVX03wzFHsTwGY4xJXrPzGFR1JDAE5yZ5jwM1IlIvIpuADUAlcAHwD2D/sBsFYwzB8wsykYeQQ7kNoeVztAAJRyWp6nzglyJyOTAA6IkzxedqYI6qWg6DMdkSNL8gE3kIOZTbEFo+RwuRzKikBlWdo6rPq+rTqjrdGgVjsixofkEm8hByKLehpcyTEJakRyUZY3JI0PyCTOQh5FBuQ0uZJyEs1jAYk8+C5hdkIg8hh3IbQsvnaCGsYTAmnwXNL8hEHkIO5Ta0lHkSwpLMLTGMMbkmcpH39TucLpyOPZwv5tiLv0HrZWLfIWgp8ySEJfB8DLnK8hiMMSZ5TeUxJHXGICJnA8cDuxHTDaWqpydY92HgNGCVqh7os3wI8DzwtVv0rKrarRtNyzBvUvN/ef+pL6z3zEPcrhuceEfj7UGwfbx0LVSNB90GUgyHjoTT7mHa4mmNfpGD/UovRMnM4DYWuAZ4E1iGk/S2napelGD9Y4D1wGNNNAzXqeppgQJy2RmDyXux+QDg9NUHueVEbKMQT1EpiMC2rU3v46VrofKfjVafduApVGxeFDUktERKEBHqGnbcQ7N1cevM3JLDpF26zhguAM5R1SnNCUJV3xaRXs1Z15iC1lQ+QKKGIUijANDgcwNkv31UjfddfdzaOWwuif66qNf6mJ+HO3IHrGHIb8mMSioC5mQqENeRIjJXRP5PRPrFqyQil4pIpYhU1tbWZjgkYzIsm/kAsfvQbb7VVhQXB96k5Q7kv2QahgdxbpaXKbOAnqp6EPA3YGq8iqr6oKqWq2p5ly5dMhiSMSHIZj5A7D7EvwHous2/wfCta7kDeS/Rbbf/GnkAHYHRIvKeiNzvXeYuT4mq/qCq693nLwOlItI51e0ak/NSyQdo1y3YPopKobhV4n0cOtJ39dEdD26UJ1AiJZQWlUaVWe5AYUh0xtDf8+iH05W0Fegbs6x/qoGISFcREff54W5sa1LdrjE5L5W5Dq77onHj0K4b/Oyh6O2d8XcYel/ifZx2D5RfvOPMQYqh/GJOHfZUo3kfxhw9hjsH35n5uSBM6ELLYxCRp3Bu4d0ZWIkzh3QpgKo+ICJXApcD9cAm4FpVfT/Rdm1UkjHGJC8to5LcPITRqrouprwt8DdVHdXU+u580U0tvxe4N2g8xuSUVPIQkhEnxyBQPP/+oPG6ew3KiczlpvjlT9hZSWYlk8ewDeimqqtiyjsDK1Q1K7fXsDMGk3Wp5CEkI06OAeUXRzcOfvEUFUODzwXk2PJMxJ2C2HkWwHIl0qXZM7i5K3cSkV0BAXZxX0ceXXCymVemN2Rj8khY8xLEyTFoVO4Xj1+j4FeepfkU4rF5FrIjyK/81ThpLAp85rNcca4XGNMyhZWHECfHoFF5qvvNwnwK8dg8C9kRpGE4Duds4Q3gLOBbz7KtwDequiwDsRmTHzr2cKa19CtPJyn2bxxicw/ixRNUFuZTiKdr264s39A4u9tyJTIrYVeSqr6lqjOA3sBU93XkMdMaBdPihTUvQZwcg0blfvEUxclcji3P0nwK8dg8C9nR5BmDe+M7r55uqkEjqvp2uoIyJq+ENS9B5AJzolFJ8eLJw1FJNs9CdjQ5KklEGnCuIURag0jl2NeoavCbqaSRjUoyxpjkpTIqqQvO3AuR0UcLcO6yuo/7uAD4AmhyLgZjjDH5o8muJFXdfksKEbkTJ8HtX54qi0VkFfBHYFpmQjQmJOlOUnv0dPj6rR2vex8Lu+7j3xXkl7gGwbp+oPmT8hjjI5kEt03AIar6eUz5AUCVqpb5r5lZ1pVk0iLdSWqxjUJTOveF1V8Eqxs7Mqm4FahGz7cQdFIe06KllODm8Slwm4hsbwDc57e6y4zJX+lOUgvaKEDwRgEaD1fdtrXxJDwNddGNAuRc4prJbcncxuJy4CWgRkTmuWX9gW2ADREw+S2bk+WEpZCOxWRU4IZBVT8Wkd44k/X0xRmZNAF4UlU3ZCg+Y8IRVpJaNhXSsZiMSurGd6q6EWcmN2MKy/G3+l9jaG6yV+9jc+8aQw4lrpnclmgGt5+JSKnnedxHOOEakyGpTJbj58IXnMbBq/exvpPgcOWH/uV+ZWc+EB3j0PucSXiaMymPMXEESXDrqqqr3OfxqCW4GWNM/mj2RD2qWuT33BgTwy8HAoLnEgTNoUgl1yKsyYRM3kt4jUFEdlLVLWEEY0xeis2BWLsEpl4R3c+/dolTBxp/Gfut71c3aL2gMQZd17Q4Qc4C1orImyLyOxE5WkSyMlObMTnLd2KcJHIJguZQpJJrEdZkQqYgBGkYrgJqgF8CbwPficgrInKjiBwmItbFZFq2ZPID/OoGzaFIJdeiJeRpmLQJMh/DQ6p6nqr2APYHrge+B64BPgC+FZHnMxumMTksmfwAv7rx1o8tD1ovmTqW22B8JPVrX1UXqOoDqjoCOAT4PU6i22mZCM6YvOA7MU6pk2PgFS+XIOhEP6lMCBTWZEKmIAS+XiAinYEhOFN9HgfsDVQBfwdmZCA2Y/JDvIlx/Mr8LvQGnegnlQmBwppMyBSEhHdXFZFxOA3BvjgNwVs4DcF7biZ0VlkegzHGJK/ZeQyuq4BvgOuAl1X163QGZ0yUfBhrn2rOgjE5LkjDcAxOF9LPgLEiUotzxjADmGENhUmbfBhr7xfj87+Kvl9RLsZtTBKCjEp6V1XHqOrxwM7AhcDX7r+fisg3IjI+s2GaFiEfxtr7xeg3J0KuxW1MEpIdlbRVVWfgjEa6DfgbTmNxfvpDMy1OPoy1TzVnwZg8EGhUkpvtfAQ7RiQdCeyEc+1hCvBmpgI0LUg+zIkQL8Z4dY3JQwnPGETkNZyEtndwsp9rgF8Be6vq3qp6sao+kdkwTYuQD2Pt/WIsbuXkLXjlWtzGJCHIGcMa4NfAm6q6KMPxmJYsH8bap5qzYEweSJjHkOssj8EYY5KXah5DuoJ4GOfWGatU9UCf5QKMA04BNgIjVXVWWPGZAvHStVA13pn+Uorh0JFw2j3Nr5fu+Q/AzixMzgvzFtrjgXuBx+Is/ylOdvW+OBe673f/NSaYl66Fyn/ueK3bdrz2fukHrZfu+Q+SmaPBmCwK7ZbZqvo28G0TVYYCj6njA2BnEekWTnSmIFSND1YetF665z9IZo4GY7Iol+ZS6A54xwEudcsaEZFLRaRSRCpra2tDCc7kAd0WrDxovUzMf5BqXWNCkEsNg/iU+V4ZV9UHVbVcVcu7dOmS4bBM3pDiYOVB62Vi/oNU6xoTgiYbBhFZJyI/BHmkIZalwJ6e1z2AZWnYrmkpDh0ZrDxovXTPf5DMHA3GZFGii89XhhKF4wXgShF5Guei81pVXR7i/k2+i1w4TjTaKGi9TMx/0NztGROi0PIYROQpnLu0dgZW4txrqRRAVR9wh6veC5yMM1z1IlVNmKBgeQzGGJO8nMhjUNVzEixXnFttGGOMyaLAF59FpJWI3C4iC0Vks4hs8z4yGaQxxpjwJDMq6U6cORj+B2gArgfuw7mX0hXpD80YY0w2JNMw/By4TFX/F9gGPK+qV+NcKzghE8EZY4wJXzINw+7AZ+7z9TgT9AC8ApyYzqCMMcZkTzINw7+BPdzni4CT3OdHApt81zDGGJN3kmkYngOOd5+PA24Xka9xbo73jzTHZYwxJksCD1dV1Zs9z6eIyBJgMLBQVV/KRHDGGGPCF7hhEJFjgPdVtR5AVT8EPhSREhE5xr17qjHGmDyXTFfSm0Ann/KO7jJjjDEFIJmGQfC/2+muwIb0hGOMMSbbEnYlicgL7lMFnhCRLZ7FxcCBwPsZiM0YY0wWBLnGsMb9V4DviB6auhV4F3gozXEZY4zJkoQNg6peBCAi1cCfVNW6jYwxpoAFvsagqrer6gYRKReRs0WkLYCItBWR0O7SaowxJrOSGa66O85kOofhXG/YF1gM3ANsBkZnIkBjjDHhSmZU0p+BFTijkDZ6yidj90oyxpiCkUwX0PHA8ar6nTPZ2nZfAXulNSpjjDFZk8wZQxnOKKRYXXC6kowxxhSAZBqGt4GRntcqIsXAjcDr6QzKGGNM9iTTlXQD8JaIHAbshDOTWz+cW2IMzkBsxhhjsiCZ4aqfAQOAmcBrQGucC88DVfWrzIRnjDEmbEnlH6jqcuDWDMVijDEmByQ8YxCRNiJyn4jUiMgqEXlSRDqHEZwxxpjwBTljuB3novMEnNFH5wD3A8MzF1Zhmzq7hrGvLmDZ95vYY+cyrj+pD2cM7J7tsIwxBgjWMPwMuFhVnwYQkSeA90SkWFW3ZTS6AjR1dg03PzufTXXOW1fz/SZufnY+gDUOxpicEOTi857AO5EXqvoRUA/skamgCtnYVxdsbxQiNtVtY+yrC7IUkTHGRAvSMBTTOLGtniQvXBvHsu83JVVujDFhC/LlLjSeoKc18JCIbL9nkqqenu7gCtEeO5dR49MI7LFzWRaiMcaYxoKcMTwKLMOZsCfyeAJYElNmArj+pD6UlRZHlZWVFnP9SX2yFJExxkQLPFGPSY/IBWYblWSMyVV2nSALzhjY3RoCY0zOSuYmeikTkZNFZIGILBKRm3yWjxSRWhGZ4z4uCTO+bJo6u4bBd79B75umMfjuN5g6uybbIRljWqjQzhjcO7HeB5wALAU+FpEX3HsweU1U1SvDiisXWG6DMSaXhHnGcDiwSFUXq+pW4GlgaIj7z1mW22CMySVhNgzdcUYyRSx1y2KdJSLzRGSKiOzptyERuVREKkWksra2NhOxhspyG4wxuSTMhkF8yjTm9YtAL1UdAEzHGSrbeCXVB1W1XFXLu3TpkuYwwxcvh8FyG4wx2RBmw7AU5/YaET1w8iO2U9U1qhpJpHsIODSk2LLKchuMMbkkzIbhY2BfEektIq2AEcAL3goi0s3z8nTg8xDjy5ozBnbnDz/rT/edyxCg+85l/OFn/e3CszEmK0IblaSq9SJyJfAqzv2XHlbVT0XkDqBSVV8ArhaR03HuxfQt0XNMFzTLbTDG5ApRje3mzy/l5eVaWVkZ6j6Dzqdw7kMzee+rb7e/HvyjTkz4ryN914dg2dA2l4MxJh1EpEpVy32XWcOQnNicA3CuB8R2/cQ2ChH77taWpd9tjlq/tFhAoa5hx/+F3zaD7tsYYxJpqmEINfO5EATNOfBrFAC+XLWh0fp12zSqUYi3Tct3MMaEwRqGJIWZcxC7Tct3MMaEwRqGJIWZcxC7Tct3MMaEwRqGJAXNORj8o06+6++7W9tG65cWC6VF0fl/ftu0fAdjTBisYUhS0JyDCf91ZKPGYfCPOvGva4c0Wn/ssIMYO/yghNu0fAdjTBhsVJIxxrRANirJGGNMYDaDWzPcMnU+T324hG2qFItwzhF78nXt+kbJbMPL9wqcjGaJa8aYXGFdSUm6Zep8nvjg34HqCtG3j42XjGaJa8aYsFlXUho99eGSxJVcsU1uvGQ0S1wzxuQSaxiStC3FMyy/ZDRLXDPG5BJrGJJULH7zDQXnl4xmiWvGmFxiDUOSzjnCd7ZRX7FNSLxkNEtcM8bkEmsYkjTmjP6cN2iv7WcOxSKcN2gv32S2P599cKBkNEtcM8bkEhuVZIwxLVBTo5JabB5D0LwBv5yFDxev4ctVG7bX2Xe3tnxdu4F6TxtbIlBSJGzetqOwdbHwxV2ncMRd/2Lluq3by3dv34qbTznAJuoxxuSEFnnGEDRvIJmchXSziXqMMZlkeQwxguYNJJOzkG42UY8xJltaZMMQNG8g1ZyFVNlEPcaYbGiRDUPQvIFUcxZSZRP1GGOyoUU2DEHzBpLJWUg3m6jHGJMtLbJhCJo3EC9nYd/d2kbV23e3tpTEnFyUiDMKyat1sVB996ns3r5VVPnu7VvxlwA5D5bvYIwJQ4sclWSMMS2d5TEElEqOgF++A9CobMwZ/TN5CMYYkzJrGFyxOQI132/i5mfnAyRsHGLzHbapNsp/8JZZ42CMyWUt8hqDn1RyBJLJd8hmboQxxgRhDYMrlRyBZPIdsp0bYYwxiVjD4EolRyCZfIds50YYY0wi1jC4UskRSCbfIZu5EcYYE4RdfHZFLjA3Z1RS5GKyjUoyxhQCy2MwxpgWKGfurioiJ4vIAhFZJCI3+SzfSUQmuss/FJFeYcZnjDEmxIZBRIqB+4CfAgcA54jIATHVLga+U9V9gD8D/x1WfMYYYxxhnjEcDixS1cWquhV4GhgaU2co8Kj7fApwvIgN4zHGmDCF2TB0B7zZXUvdMt86qloPrAV2jd2QiFwqIpUiUllbW5uhcI0xpmUKs2Hw++Ufe+U7SB1U9UFVLVfV8i5duqQlOGOMMY4wG4algHcQfw9gWbw6IlICdAS+DSU6Y4wxQLh5DB8D+4pIb6AGGAH8IqbOC8CFwExgGPCGJhhPW1VVtVpEvkkhrs7A6hTWzyV2LLmpkI4FCut4WvKx9Iy3ILSGQVXrReRK4FWgGHhYVT8VkTuASlV9Afgn8LiILMI5UxgRYLsp9SWJSGW8sbz5xo4lNxXSsUBhHY8di79QM59V9WXg5ZiyWz3PNwPDw4zJGGNMNLtXkjHGmCjWMMCD2Q4gjexYclMhHQsU1vHYsfjI+3slGWOMSS87YzDGGBPFGgZjjDFRWmzDICIPi8gqEfkk27GkSkT2FJE3ReRzEflUREZnO6bmEpHWIvKRiMx1j+X2bMeUKhEpFpHZIvJStmNJhYhUi8h8EZkjInl9r3sR2VlEpojIF+7n5shsx9QcItLH/f+IPH4QkWtS3m5LvcYgIscA64HHVPXAbMeTChHpBnRT1Vki0h6oAs5Q1c+yHFrS3JsmtlXV9SJSCrwLjFbVD7IcWrOJyLVAOdBBVU/LdjzNJSLVQLmq5n1CmIg8Cryjqv8QkVZAG1X9PttxpcK9g3UNcISqppL023LPGFT1bQrkdhuqulxVZ7nP1wGf0/gGhXlBHevdl6XuI29/vYhID+BU4B/ZjsU4RKQDcAxOQi2qujXfGwXX8cBXqTYK0IIbhkLlTm40EPgwu5E0n9v1MgdYBfxLVfP2WIC/ADcADdkOJA0UeE1EqkTk0mwHk4K9gVrgEbeL7x8i0jbbQaXBCOCpdGzIGoYCIiLtgGeAa1T1h2zH01yquk1VD8a50eLhIpKXXX0ichqwSlWrsh1LmgxW1UNwJtv6ldsdm49KgEOA+1V1ILABaDSjZD5xu8NOByanY3vWMBQItz/+GWCCqj6b7XjSwT29nwGcnOVQmmswcLrbN/808GMReSK7ITWfqi5z/10FPIcz+VY+Wgos9ZyJTsFpKPLZT4FZqroyHRuzhqEAuBds/wl8rqr3ZDueVIhIFxHZ2X1eBvwE+CK7UTWPqt6sqj1UtRfOaf4bqnpelsNqFhFp6w5swO12ORHIyxF9qroCWCIifdyi44G8G6gR4xzS1I0EId9EL5eIyFPAEKCziCwFblPVf2Y3qmYbDJwPzHf75gF+4960MN90Ax51R1gUAZNUNa+HeRaI3YHn3Jl2S4AnVfWV7IaUkquACW4XzGLgoizH02wi0gY4Afhl2rbZUoerGmOM8WddScYYY6JYw2CMMSaKNQzGGGOiWMNgjDEmijUMxhhjoljDYAqeiIwUkfWJa+amVOIXkWNFZKE7/DcjRKS/iNQUyG0lDNYwmJCIyHgRUfdRJyKLReRPyXyZuNvISE6De0vp6zKx7SzHMRa4S1W3pXGbUVR1PvABcG2m9mHCZQ2DCdN0nAS2vYFbgCuAP2U1ogImIkcBfUnT/XMSeAS4XERabNJsIbGGwYRpi6quUNUlqvokMAE4I7JQRA4QkWkiss6dROkpEenqLqsALgRO9Zx5DHGX3S0iC0Rkk/uL+48i0jqdgTcVm7t8vIi8JCKj3W6V70TkETcrNVKnrYg8JiLrRWSliNzsrjPeXT4D6AmMjRxjTAzHi8gnIrJBnImZeicI+xfAdFXdGLOdU0XkQ/f9WiMiL0beL/f9u9U9nnUiskREzhZnYpun3di/FJETY/b1GtAJ524CJs9Zw2CyaRPOfAuRyYbexrn/zuE490hqB7wgIkU4ZxaT2HHW0Q14393OBmAUsD/OWcgI4LfpCjJAbBH/ARzoLj8bOBPwzqb3P8CxbvmPgYPcdSJ+hnODtzs8xxixE3AzznEeCewMPJAg9P8AomZaE5GTgeeBfwGHAscBbxH9XXAN8BHOjeUmAY8CTwIvAwe778UT3sZXVbcCc9zjM/lOVe1hj4w/gPHAS57XhwOrgYnu6zuA12PW2QVnDoDD/bbRxL4uAxZ5Xo8E1idYpxq4Ls6yoLEtAUo8dR7C+cUOTkOyFRjhWd4W+A4Y31QcbvwK9PGUnetur6iJY/oeuCim7D3g6QTvw1Oe1+3cff/VU9bLLSuPWfdZ4PFs/63ZI/WH9QeaMJ3sjq4pwTlTeB7nZmbg/Ho9Js7omx/h/IL1JSLDcH7l7oPzRVbsPtIlaGyfqWq9Z9ky4AhPvVJPXVR1gwSfc3yLqi6I2XYpzplDvJkIy4DNMWUDcRqxpszzxLheRDYC8z3LI7d23i1mvU3uPk2es4bBhOlt4FKgDlimqnWeZUXANMBvRE7ce8yLyCCcuQ5uB36N8yv5dNJ7UTtobHUxy5QdXTTiKWuO+pjXke001R28GufMJll+x1EX89pv351wzjhMnrOGwYRpo6ouirNsFvBz4JuYBsNrK43PBAYDNap6Z6RARHqmHGnysSWyCOfL9XDga9h+u+QDga889fyOsblmAwf4lB2P082VbgfidCeZPGcXn02uuA/oCEwUkSNEZG8R+YmIPBiZIAbn1+iBItJHRDq7s9YtBLqLyLnuOpfjTFrSHHuIyMExj84BY2uSqq4HHgb+2x1ddADwD5zPoPcsohr4DxHp7u47Fa8CR8eU3QUMF5Ex7kirfiLya+/oqeYQZ67x7jijk0yes4bB5AR1po0cDDQArwCf4nwhb3Ef4PzK/RxnpE0tzhzEL+Ikcf0Fp2/8BODWZobxa5xf1N7HiICxBXEd8A7wAvCmG28l0dcBbgX2xDmLqG3mcUQ8AewnIv0iBepM3nQmzlSQs3FGJB2Hc2ypOAd4TVW/SXE7JgfYRD3GZImI7AR8A4xV1f/J0D7uBrqo6sWZ2L67j52AL4FzVPW9TO3HhMfOGIwJiYgMFJFfiMg+IjIQJz+gPTAxg7v9PbBYMnivJJykvLusUSgcdsZgTEjcxuAhoA/OKKM5ODkLVVkNzJgY1jAYY4yJYl1JxhhjoljDYIwxJoo1DMYYY6JYw2CMMSaKNQzGGGOi/H/uD3Kaj1VJyAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the data color coded by species\n", "for ind in range(n_labels):\n", " plt.scatter(iris.data[:, pl_ind][iris.target==ind],\n", " iris.data[:, pw_ind][iris.target==ind],\n", " label=iris.target_names[ind])\n", "\n", "# Add title, labels and legend\n", "plt.title('Iris Data: Petal Length vs. Width', fontsize=16, fontweight='bold')\n", "plt.xlabel('Petal Length (cm)', fontsize=14);\n", "plt.ylabel('Petal Width (cm)', fontsize=14);\n", "plt.legend(scatterpoints=1, loc='upper left');\n", "\n", "# Note that the data, colored by label, can also be plotted like this:\n", "# plt.scatter(iris.data[:, pl_ind], iris.data[:, pw_ind], c=iris.target)\n", "# It is, however, more difficult to add labelled legend when plotted this way" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this data, we know we have 3 different 'groups' of data, which are the different species. \n", "\n", "As we can see in the plots above, these different species seem to be fairly distinct in terms of their feature values.\n", "\n", "The question then is whether we can learn a clustering approach, based on the feature data and without using the labels, that can learn a meaningful grouping of this data. \n", "\n", "If this approach works, then we might be able to try to use it on other data, for which we have feature data, but might not be sure about the groupings present in the data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Apply K-Means Clustering\n", "\n", "Clustering is the process of trying to learn groups algorithmically. \n", "\n", "For this example, we are going to use the K-means clustering algorithm. \n", "\n", "K-means attempts to group the data into `k` clusters, and does so by labeling each data point to be in the cluster with the nearest mean. To learn the center means, after a random initialization, an iterative procedure assigns each point to a cluster, then updates the cluster centers, and repeats, until a final solution is reached. \n", "\n", "
\n", "K-means is a clustering algorithm that attempts to learn k clusters by grouping datapoints to the nearest cluster center mean. \n", "
\n", "\n", "
\n", "For more information on K-means, see the article on \n", "wikipedia\n", "or on the sklearn\n", "user guide. \n", "
" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Pull out the data of interest - Petal Length & Petal Width\n", "d1 = np.array(iris.data[:, pl_ind])\n", "d2 = np.array(iris.data[:, pw_ind])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Whitening Data\n", "\n", "In this example, we are using two features (or two dimensions) of data. \n", "\n", "One thing to keep in mind for clustering analyses, is that if different dimensions use different units (or have very different variances), then these differences can greatly impact the clustering. \n", "\n", "This is because K-means is isotropic, which means that it treats different in each direction as equally important. Because of this, if the units or variance of different features are very different, this is equivalent to weighting certain features / dimensions as more or less important.\n", "\n", "To correct for this it is common, and sometimes necessary to 'whiten' the data. 'Whitening' data means normalizing each dimension by it's respective standard deviation. By transforming the data to be on the same scale (have the same variance), we can ensure that the clustering algorithm treats each dimension with the same importance. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Check out the whiten function\n", "whiten?" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Whiten Data\n", "d1w = whiten(d1)\n", "d2w = whiten(d2)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Combine data into an array to use with sklearn\n", "data = np.vstack([d1w, d2w]).T" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# Initialize K-means object, and set it to fit 3 clusters\n", "km = KMeans(n_clusters=3, random_state=13)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,\n", " n_clusters=3, n_init=10, n_jobs=None, precompute_distances='auto',\n", " random_state=13, tol=0.0001, verbose=0)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Fit the data with K-means\n", "km.fit(data)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEcCAYAAADDfRPAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZfbA8e+ZPgkhlIQWSkCkF0WkIyCiiIrY6yq6rm1tu6uu67qWXd2fu669l7WtunbXXlAERQEFFASkN+klkJ7JlPP7406SmcmkTDIpwPt5nnmYufe9954JyZy5bxVVxTAMwzDK2Jo6AMMwDKN5MYnBMAzDiGISg2EYhhHFJAbDMAwjikkMhmEYRhSTGAzDMIwoJjEkgYjcLiIaftyewHHPRxyXnaRYNkScMyQiJSKyTUS+EZHbRKR9Pc8/Lfx+bxeRVsmIuZbXHR/xviIfeSLyrYhcFFO+Tv8nzUUV77VYRJaLyF9FJCVcrnfE/k9jznFrxL6XYvY9FLFvSi3icYjIdBH5VER2ikipiGwVkdkicp2IpIbLTY847/Qk/kjixTS97HexIa+TDDF/6+ObOp6aOJo6AKNBCeAGOoQfo4BrReQsVZ1Rx3NOAy4MP38e2FffIOspDRgJjBSRQar6uyaOpyF5gL7AX4DjRGSMqq4Ukd1ABjBCRGyqGgqXHxVx7OiYc5W9VmBudRcVkQzgPayfc6SO4cdRwCzgx8TeTr1NB8aFn9/eyNc+oJk7hiYgIh4AVZ2uqhJ+bGiAS3XHSgwDsT7EAVoD74hI3wa4XmOZraqC9UE5PWL7tcm682pOwu/VDowA9oQ3DwPOCT//NvxvS6z/a0TEFi5fJltEOoX3tQAGh7cvV9W9NYTwBhVJYRkwEfACrYApQF2/ZDRr4bske1PH0RRMYmhAMbePY0XkTRHJBX6Osz874rjfiMgCEckREZ+IbBGRGSJyYRWXqpKqlqrqUlW9CHg3vDkVuC3iemeFz/+LiBSFr7lORJ6IrHoSEaXibgFgfWT8ItJCRF4QkZ9EZI+I+EVkn4h8JSJnxfn5lB27IdH3FX5vPlV9Afip7JTA0LqcKyau/4XjCpZ9mIa328LVJypWlZ2IiFdE/k9EVopIgYgUhn92b4nIiOqukwhVDanqfOCViM3Dwv9+E7Gt7E5gAJAO7AIWxOwbjpVooCKpxCUixwPjwy+LgONUdaaqlqhqrqp+rKrHAktrOE/Z//WsmraLSIaIPBr+ORaJVV24UkT+G646yw7/Lo6Lcx6N2NZaRO4JH1sSPs9sEZkWE0Nk9dflInKviGwFSoEu4TKdReRxEVkvVjXaXhH5WESOivNeTxSRJeFrrhCR86v72TRHpiqp8bwDtA0/r/IbmoicATwVs7lT+JELvFCPGP4BnBx+PiWi2mECcExM2e7AZcD4cBVNaS3O3wK4IGZbOjAWGCsinvAHebIl+wvOs1g/JxtwFnB/ePt4rKoTgOdVVUXkX8CVMcd3Dz++BOYlObZ47zU2MTxGRTXSt8B6rIQ5Cuvb/+gqjo3nxIjnr6jqlniFVDVQw3kS8QLWnUikNKAX8DI1JCEAEWmH9d4Pidjsxqr2OkpEblTVe+IceicVf6dl5+oNzMGqrivTCpgMHCsi56rqa+GyE4H/UZF4ewP/AbbVFHNzYu4YGk8e1u24l8q/9JHKvoEUYP1SuYFuwJnAJ/WMYUXE8zQq/gBewfoWmQE4gfbAc+F9vcviDVdpRH6wd4+pCsvH+iDNBlKwqnpGYX3TBEhq/b+IuMJ3Uf3DmxT4Pgmn/gjYHn5+bsT2sudKRdVc2f/XPCAT626sD1ay+DkJsQDldyvDqag+Apgf/ncB4As/HxP+t+zD/xsqPvxHx5SBmhND94jny2sdcP2U/Uzfxvpi0RIYBPwB2KyqG8K/i7PLDoj4PZTwpr9iJYUgcBrW311n4Kvw/jtFpCzJR2qB9TNuAfQEdgIPYv1t5GJ9ifIAh2L9PdmAR0TEVXZeKpLCzeH4z8Jq49tvmDuGxnOLqpZ9e6zuD2x9+N9U4BZgIdYHzGeqmlvPGKr6IrANuBXr1rwDVnKI1LuW5y/C+gN6DauRNA2reifueSL+iBM1LrLKIMLDqrqxjucsp6oBEfkPcAMwVER6ApuAU8NFZka0Ca3HqrYpaxReGn78u5Z3WTWq4r0uAF4Nx+sTkYVYSbiriHQmOjGsCz8/XETSqGh72Kmqa5IRY5Ktx2orGYn1N7Acq7rwgYiG9ZqcFP7XDrwVZ78L6/f91ZjtL6pq2ba1IuKl4m46HesuMFYGMEREfqKiem838I9wvK+LyNVEJ+RmzdwxNJ4falnuMazb/RDwK+AB4FNgh4jcVM8Y+kQ8zwP2iEg61m3y+Vj1qbFJAaxvW7XxR+BRrLuPlkQnBbC+aSVbIdY350uB65J43mcjnp8LHI/VcA/w74h9v8e6S0kHrsGqBvwW2CIixyYxHrDuClYAfwcmqKo/Yl/kN//Tsb7p+4CFqrodKzk4gF9jJWyooX0hbH3E86R2WBCRqr6Y/gZYiVVtdwPW3esCYJ2IHFbL07erRZm2cbbF/p22oeIOoKZztabiM3VbTBLbXItzNBsmMTSe4toUCjfqnYn1CzkGuBjrg88N/F1EsuoRw58inn8Y0b5Q9kf0BdAx/E3+mqpCrOb8Z0c8nwa4w+faU0X5upodUXXQQlVHqOrTmsQ55FV1BRXtA+dQUYWzD6u9qKzcGlUdhtUGNAkrOW3D+hb5cJJiKXuvHlXtq6p/VtWCmGKRiaGsym6hqpZVMc2J2Rd7TFU+iHh+XhXVL0jNvXfK7p4ivxz0iFdQVeerah+s5HY8cBNW1Wo3rHay8qLVXG9n+N8Cwr+HMdVNNlV9NM5xsX+nOVjVUQCrY88Tca4PsdoOy5JBR7F6hpXpXE2szY5JDM2MiJwmIlcBWcBirLuHxWW7SfAXTEScIjJARJ4HTghvLsSqgwWIbDQsAQpFpD9wdRWnjPyQHywikXcFkefaBzhF5C/E/2ZW715JCegpIpNjHqNqPqz8rqEPVj01WA2wJWUFROSGcI8rD/A1VjXa1vDurhHlGnqAU+SHfNc42+bE7INa3DGo6sdU1MunAJ+KNdjQIyLpInK8iMwg3E22GmVVfANFpJuIOIG/xSsoIneJyElYH8gzgdep6LARGf+eiGNi7yTKEloL4BkR6SIiKSLSX0QupeJvqlqqWoz1hQngUBH5p4i0E5E0ERksIr8v26+qhcB34bIZwB/D5c6k8jiS5k1VzaOeD6zBNRp+3B6x/fmI7dlxjqu0H6tOVat4bAW8NcSyoZrjFesb0KSI8q2xvl3FlltVxXs6PU7ZDeF9f46zbxfWH7Vav25RsUYdX8P7Gh9RflaC/yfxHj/W4hxpWEk08rgjYsp8Xs01/lfF//X4Wly7/DwJ/B6ujLn+yRH7+sXsK8H6Jl2b82Zi3bVW9/M8LFx2esS26RHnuC1iux+rPSryZzsrouyaaq7zQES56+PsnxXe1x6r+qzKmCPOEzfmiP19sZJQVefaEFF2IlZCiy2zO5H//6Z+mDuG5ucLrF5Ca7Bug4NYVROvAuPU+gZTW4pVz7wd69vhbUBfjRj1rNbgpuOxvlEWYSWf24G7qzjnW8BdWI2xwZh9/8Cq/96CdUs+GzgaqzfHfkdV84E3IzYtUdWFMcVewOrFtBnrw9aP9X93H1YbUWOKrRqKvCP4mei7vchqpmqp6i6sas1fYyXC3VjvczvWXdLvgdU1nOZurG6/ZeMDvqbqb9GPYN0plJUtwRpYdxtWm0OZR4EnsP4+NCbmHVhddP+J1S7jw/p7Wo319xXZu6taqvozcBjwOFayKcX6nV6O1d50eUTZL7CqUZeGy63Bav/6gP2IhLOcYRiGYQCmjcEwDMOIYRKDYRiGEcUkBsMwDCOKSQyGYRhGlP1+SoyMjAzNzs5u6jAMwzD2KwsXLtytqpnx9u33iSE7O5sFCxbUXNAwDMMoJyJVzitmqpIMwzCMKCYxGIZhGFFMYjAMwzCimMRgGIZhRDGJwTAMw4hiEoNhGPstfzDIlrw8iv3+mgvHsbuoiF2FhVHbcktK2Jafz8E8j1yjdVcVkS7Ai1hLR4aAp1T1wZgy44F3qVg16m1V/SuGYRgxXlz8A/fOnUMgFCKkcGb/AdwydjxOe80Lrm3Yt5frPvmQFXt2A9C9VWtuH3c0jy/4jnmbf0FEaJvi5Z/HTGZUl641nO3A02izq4ZXfuqoqovC684uBKap6vKIMuOB61X1xNqed+jQoWrGMRjGweWj1Su5YcYnFAcq1obyOByc3X8gt447utpjfYEAY557mr0lxYQiPv9sIthECIQqVuT0Ohx8cO4FdG/VOt6p9msislBVh8bb12hVSaq6TVUXhZ/nY80PX59lKg3DOEg9/N28qKQAUBII8Oqyn/DFbI/1+bq1lAT8UUkBIKQalRTAqqp68cfaLtd+4GiSNgYRyQYOx1oVKtZIEVksIh+Hl5iMd/ylIrJARBbs2rWrASM1DKM52lEQu+S1RVXJLy2Nu6/M1oI8SoOxa0zFF1BlQ+7emgseYBo9MYhIC6xVwK5T1byY3YuAbqo6GGsh9f/FO4eqPqWqQ1V1aGZm3Kk+DMM4gA1s3z7u9lSnizZeb7XHDmrXoVbtEAAeu4MRWV0Sjm9/16iJIbwA+FvAy6r6dux+Vc1T1YLw84+wFpPPaMwYDcNo/m4cNRavw4FEbPM6HNw8dhw2kSqPAxiW1Zn+me1w2yv63rjsdtLdHryOim0Om400t5uzBwxKdvjNXmM2PgvW+rg5qnpdFWU6ADtUVUVkGNZ6u920miBN47NhHJx+3rWT++Z9w5Id28lKS+ea4SMZn929Vsf6AgGeXPg9b/68lFBIOblPXy4/YhjvrvyZ5xcvorC0lIndD+Ha4aPITE1t4HfSNKprfG7MxDAGawHwn7C6qwLcDHQFUNUnROQq4AoggLWY/O9V9ds4pytnEoNhGEbiqksMjTaOQVXnANXe46nqI8AjjRORYRj1tWTHdu6YPZMlO7bTwuXigsGHc/WwkThsya2lvunzT3nr52UEVXHYbFw0+HD+NHZ8Uq9hVNjv12MwDKNprNubw7lvv05ReNRxrs/H04sWsDU/n3smTU7ada7/7GPeXlE+3IlAKMTTPyzEabdz/aixSbuOUcFMiWEYRp08tfD7SmMGSgIB3l+1gl1FhVUclZhQKMQ7EUkh0tOLTBVyQzGJwTCMOlm2ayfBOG2UbruDjfv2JeUahf5SqmoF9ccMRjOSxyQGwzDqpG9GJvY4XUN9wQDd0lsl5RqpTleVDZPOJLdjGBXMT9YwjDq57IgjccUMFPM4HEw5tFfSunjabDZO6tUn7r6LDz8iKdcwKjOJwTCMOjmkTVtePvVMBrZrjwAtXC6mDx7CPyYel9TrPDD5BE7p07d84JpdhIsGD+GPo49K6nWMCo02jqGhmHEMhtH0VBWpYcRxMoRCIWymCikpmsU4BsMwDjxb8/N4cP5c5mzaSBuvl0uPOJKjunbj8QXf8dHqVbgdDs4bOJgz+w3g+cU/8ObypQRVmdanL78ZciQfrFrB84sXUVBayjE9enLVkSPISEmp8nqRSWHVnt3cP+8bFm/fTpf0dK4aNoKxXbMrHZPnK+Gx7+dHxfOrQYdhj0kwIVXeWPYTz5XF0/0Qrho2stp4msKMtWt4bMF8dhQWMDyrC9cOH0l2kqcFN3cMhmHUyc7CAia//AL5Pl957ySP3Y7b4aQ44C+fwdTrcOB1OCkK+CkJd2912+2kOJ2UBALl02c7bTbaeFP49PzptHS7q732z7t2csabr1Ls95f3WvI6HPx94rGc3LtvebmSgJ/jX36RbQX5UfFMyO7BI1NOijrnLTNn8M6K5XHiuZCWbk/9flhJ8vyPi7jn26/LY7SL4HU6ef/sX9GtVWIN/s1iPQbDMA4sTy9aQGFpaVSX1ZJgkFxfSdS01sWBADklxeVJAcAXDLK3pCRqTQV/KESur4TXli6p8dr//PZriiKSQtl17vxqVtQ6C++tXMGuwsJK8czcsI41OXvKt23Lz+fNn5dVEc9PNf8wGoEvEODeuXOiYgyqUuT388j385J6LZMYDMOok/mbf0n6WIKSQIC5mzfVWG7xju1xt+eX+sgpLi5/PX/LZooCldeDtonw4/Zt5a+X7txRqYdVWTzf1iKexrApNzfu9pAq32/ZnNRrmcRgGEaddG6ZXv3kZ3VgF6FrLcZAZKbE7w5rQ0hzucpfd01Pj/uBL0CHFmnlrzukpREMVa5Wr208jSEjJaXKRNwprWVSr2USg2EYdXLpEUfidkT3X3GIrdJ6CGVrKceyIThitrvsdi4YfHiN175q2IiotRPAGkNxRv8BUTGd2W8gDon+mLOL0MabwqguXcu3DchsR7dWrSpN/ue027mwFvE0htZeL8d0PwR3TKLzOhxcceSwpF7LJAbDMOrksA4d+dekybTxePE6nLjsdsZld+ffJ51CVlpLPA4HLrudwe078NIpp9M3IxO33Y7bbqdH69a8eMppHJnVGZfdjsfhoH1qC5444WR6tG5T47VP6tWH348YTarTRYrTidtuZ1rvvvw5ZsbVjmlpPD/ttKh4BrXvwH9POzMqWYkIL047nSM7ZUXF82Qt42ks90yazLGH9MRlt5PicNLS7eb2cUfH7Y1VH6ZXkmEY9RIMhdiSn0dLt5tWHmtZTVVla34+Loc9qtpnR0EBQQ1FVX3sKSqiyO+nc8uWCY+F8AUCbC8ooG1KCi0iqpBiVRVPPPWJp7Hk+UrYW1xCp7S0Wi9TGqtZLNTTUExiMIwD36bcfTz7w0JW7tnN4PYdmH7YkKg2gjIFpaW8unQJMzeso31qKhcOHsJhHTo2QcSJUVU+W7eGN5YtJaAhTu3TjxMO7V1prEUymcRgGMZ+a/GO7Zz39uuUBoMEQiGcNjseh523zzyXQ9q0LS+X7/Mx9dWX2FFYQEkggABuh4M7xh3NGf0HNt0bqIU/fv4pH65aWd6DyutwMqZrV5444eQGu2sx4xgMw9hv3TJzBkV+P4Fwjxx/KEhBaSl3fj0rqtyLS35ge0F++XgJxepuesdXX1ISp8tqc/Hz7l28v2pFVLfa4oCfbzZt4vutW5okJpMYDMNotkqDQX7evavSdgXmbf4latuMtWvwRQxkK2MTYdmunQ0VYr19s2kjwTjdUIsCfmZvXN8EEZnEYBhGM+aw2apcdyE1prG5tdcbt1wwFKJVM5nSIp50jyduA7Lbbqe1J/57amgmMRiG0WzZRDi1b79Kffc9DgfnDzwsatv0wUMqjW0oG6AW2RbR3Bx3yKFxt4sIU3vHX4uioZnEYBhGs3bL2AmM6NwVt91BmsuF225nUo9D+O2Rw6PKjcvuzlVHjsBtt5PmssY3ZLdqzTNTT2miyGunpdvNs1NPpZXbQwuXq/zx2JSptEtt0SQxmV5JhmHsFzbs28vGffvo2bYtWdVMAZFbUsKSHdtp7fXSP7Ndsx2LEMsfDPLD9m0EQyGGdOxUaVR5spn1GAzDaJZ+yc3lfyuXk1tSwoTuPRjVuSt7iot5Z8VytuXnMSyrC8f0OISSQIBvftnEmpw97Coq5MRevRGET9au5sft28hOb8W0Pv1I93hI93gY2y27/BqqyoJtW5ixdg1eh5OT+/Rt0tHMVcXjtNsZltW5yeKKZO4YDMNoEh+tXsn1Mz4hGArhD4VIcTrpn9mOZTt3EFTFFwyS4nSSlZbG7qJiSgJ+igMBUhxOWrhduOx2coqLKfL7y6e7eO30s+ndNqP8GqrKjZ9/wkerV1MS8GO32XCIjVvHTeDsAYMa/T2rqjVmYfWq8njs4XjOaeR4zDgGwzCalSK/nxtmfEJJIFA+Y2iR38+CrVsoCgTKu50W+f2syclhX0lx+ToERQE/uwoL2ZKXR5Hf6vtfEgiQ5/Nx/WcfR13n282b+Hj1aooD1toNgVCIkmCAO2bPZG/E9NyN5dvNm/ho9aqoeHzBAH9toniqYhKDYRiNbv6WX+JO9xCv/kLjbI+3DazlPnNLSspfR44mjuSw2fhq04YEIk6Oj1ZXHc/sjY0fT1VMYjAMo9HFTm+dLApRs6Y67fa4U36DVDk+oiE5bVXH47I3n4/j5hOJYRgHjeFZXeJ+QMb7yIy3noNNpNKHl12EoR07kRaxXvSpffrFXagnpCHGdetel9Dr5ZS+/ZtVPFUxicEwjEbnstt56sRppDqdpDqdeOwO3HYHp/btRxuPh1SnNV7B63Awpks3uqW3Kt+W4nTSp20Gh3fsVL4WQ6rTRfvUFvzr2OOjrjO4Q0euGDoMd3iNhRSnE6/DwaNTplYaOd0YBrfvwJUx8XgcDh6ZclKTxFMV0yvJMIwmU1Bayufr1pJf6mNM1250b9Wa0mCQmevXsbOwgCEdOzGgXXtCqszZtJH1+3Lo3TaT4eFunYu2b2XZzp1ktWzJuG7dq6yi2pKXx6yN6/E4HEzqcQgtm3iKjOYQj5l22zCMhIRUmbt5Exv27aN32wyO6NgJgCU7d7B05w66tExndJeu2G021ubsYf6WzbTxpjAhuztuh4Nt+fl8tWkDHoeDo7N7RFXvHKhUtfzn0zmtJWO6dktoPYV9JcXMXL+OoCrjs7vXuKBQfTWLAW4i0gV4EegAhICnVPXBmDICPAhMAYqA6aq6qLFiNAwDcoqLOPut19mWn0dQFZsIPVu3xetwsGTndpSKdZMHt2/P5+vWIQI2seG02zi5d19eXboEm9iwCajC4ydMjRp0dqDxBQL8+r13+GH7NhTFLjZaez28fvrZcRcUivXR6pX84bNPsNsEFG7VELeOO7rRxzaUabQ7BhHpCHRU1UUikgYsBKap6vKIMlOAq7ESw3DgQVUdHveEYeaOwTCS68oP3+OL9WvLxxeAlQhUlcjJoa3vwkIopuOoULkraYrTyXeXXEGK09kwQTexB+Z9y5MLv4ua9tsuwrCszrx86pnVHru7qIijnnuakmAgarvbbueT86bTrVWrBom5WQxwU9VtZd/+VTUf+BnIiil2MvCiWuYBrcIJxTCMRhAMhfg8JikABGOSAli3/bFJAeKPL7AhzNrQNGsLNIbXl/9UaS2IoCoLtm4h3+er9thP164mXg/WkCofrl6ZzDBrrUl6JYlINnA4MD9mVxYQufrGZionD0TkUhFZICILdu2qvIiHYRh1E1Il1AC1CIpSGmcRnQNFIFh5oZ3yfXEW4YnkDwbj/syDIcUXCMQ5ouEllBhEpLuITBCRKSJypIgk3JQuIi2At4DrVDUvdnecQyr9xFT1KVUdqqpDMzMzEw3BMIwqOO12juyUVekPUYj/x1nbeUv9oRBHdetWv+Casck9D407YO6QNm2rXECozITsHnG3ux12Jh3SMynxJarGxCAi2SLyDxHZBKwBvgA+wPq2v09EZojIGSJSm3M5sZLCy6r6dpwim4EuEa87A1tr8T4Mw0iSv088lnSPp3zRmxSHk4yUFDq2SCPFYbUReBwOWrhc9G6bWd5u4LTZcNvtjMjqTIrDiWDVs3scDm4ZO5423pSmeksN7ncjRsf9+dw7aXKNx3Zr1YrfHjkcj8OBDUEAr8PBWf0HMaBd+waOPL5qG59F5EFgOvAZ8B7wHdYHdTHQBhgAjAXOAQLARar6fRXnEuAFIEdVr6uizAnAVVQ0Pj+kqsOqewOm8dkwki/PV8LbPy9ndc4eBrZrz9TefbGL8OHqlSzatpXsVq05rW9/0txuZq5fy5xNG8lISeX0fv3p2CKNeZt/4dO1q0lxujilTz8Obdt8V1BLFl8gEPXzObVvv4SS4bKdO3h35c8EVTnh0N4MCXcRbih1HscgIvcA/1DV3bW4yBQgRVXfrGL/GOBr4Ccob8e6GegKoKpPhJPHI8BkrO6qF6lqtZ/6JjEYRuPZVVTIqj27yUprSXar1oA1WOujNavo3LIlx/Xoia2KvvvBUIjFO7YT1BCHte8Yd53jpqaqLNu1k4LSUga173DA9qKCeoxjUNUbansRVf2ohv1zqKFKUq0s9dvaXtMwjMYRUuX2WV/w+vKluO12/KEQQzp0wm4Tvt60sbyc227n9dPPZmD7DlHHL9y2hcs+eBdfIIgANpvw8PEnMrZrduO+kWqs37eXi959i91FRdhECIZC3DbuaM7sP7CpQ2t0Zq4kwzBq9NKSH3nr52WUBoPkl5ZSEggwf8svUUkBwBcMcuabr0Vty/f5mP6/t8gpLqbQX0qBv5Q8n4/LP3iXXUWFjfk2qhRS5VfvvMEvubkU+f0UlJZSHAhw++yZLNmxvanDa3S1Tgwi0lpEHhSRJSKyXUR2Rj4aMkjDMJrWsz8uLF8op0ywimpoXzDA7IgxC5+sXR13bENIlfdXrkhmmHW2aNtWcktKKsVZGgzy8k+LmySmppTIlBgvAv2xGpB3EH8ci2EYB6C8GgZpxdqSX9ETPbekBH+cMQy+YJCcZrJq2b6SYiTOKLOQarO5q2lMiSSG8cA4M3eRYRx8xnTpxkerV1Ya/VyVKYf2Kn8+KjzZXuxo6hSHkzFdm8fYhiEdO8VNXl6Hg0ndD2mCiJpWIm0MaxMsbxjGAeKGUWNJc7tx2ayeRGXjE+J9IEzr3ZdWnopBXf0y2zH5kEOjevh4HU6Gd+5cPn12U2vjTeGqYSPKx26ANRaha3orTunbrwkjaxq1nkRPRMYBtwDXA0tVtVmMbzfdVQ2jcewsLOD5Hxfx/dYt9Gjdhl8ffgSpThd/njmDRdu20sLt4sqhwzl/0GGVjg2p8vHqVby+/CeCIeW0vv2Z2rtPQtNSN4Y5mzbynyU/sq+kmON79uKs/gPxHqBdVpOyHoOIZAGvASPj7VfVJumUbBKDYdTfur05rM7Zw1FduuENryS2r6SYtTk59GrblrRqFpIJqbKtIJ80l5uWSVp3YUdBAVvy8xiQ2Q5X+Ft8ScDPnqJiMlJScDuqrgVPJJ7dRUWoKpmpia99UNt4mqtkrcfwXyAduAbT+GwYB4TtBfkc8+KzFEX0OBqRlUVxIMjiiG6aY7t247mpp1YavPbpmtX8ZdbnFJSWEhdV9TcAACAASURBVFLl6Owe/OOY4+q8MM++kmKmvvoSm/OsxmsBzuo/kHSPhxcW/2DN2STCFUOHceXQ4ZUajD9bu5q/fPk5+eF4JmR355/HTK4Uz/p9e7nukw9Zuccau9u9VWvuP24KfTJqnnstpMo9337Ni4t/AKz1py+vIp79VSJ3DEXAMFVd2rAhJcbcMRhG3fV79MFK6wBUZWqvPjww+YTy14t3bOect16jJCKpuOx2hmd15oVpp9cpntHPPsW2gvxK2x0iBCI+q7wOB38eO55zBw4u37Zkx3bOrkU8vkCAMc89zd6S4qhZTVu63Xw1/Tc13mU8NP9bnlz4fVT3Xa/Dwc1jxnFenGq05ipZ6zEsB1omJyTDMJraj9u31TopAJXWBnh64feVpoUuDQb5bstmtuTFTpxcs/V798ZNCkBUUgAoDgR4bEH0rP1PLapdPJ+vW0tJwF9pqmt/MMgHq6ofV6GqPPND5TEdVjzfVXvs/iSRxHALcJ+IHCMi7UWkTeSjoQI0DKNhLNqe2MTFsQPaNubui1uf7LLbq/yAr87qnBqnZIuyu6go6vWmfbWLZ2tBXty1IYoDgajxF/H4QyEKS0vj7ttTXBR3+/4okcTwETAMa6bVrcCu8GN3+F/DMPYjU3r2qrlQBK8junfO8KzOcdcgKA0G6VWH2VSHJdh1tW9Me8CwWsYzqF2HuBP4pTidHNa++gUjXXY7nVumx93Xp+2BszZMIolhQsTj6IhH2WvDMPYjHVqk0bN17W/2bxozNur1JUOGkuJ0YotocPU6HPz68CNoWU0vpqq08ng5Ort73H2xH/geh4M/jRkXte03Q44kxemqFM/Fh0XHMyyrM/0z2+G2V/S9cdntdEtvxYTu8RfNiXTrURPwxPRC8jgc3Dx2XBVH7H9q3fjcXJnGZ8Oon0vee4cvN6xDsWZHvWP80RSU+nlg/rcUlpaS7vHwp9FHcUacWUY35+Vy/7xv+eaXjbT2eLl0yJFM69O3Xr1z7vpqFi8vXUxpIEi7FqncPfFYXHYH98/7hg379tI7I5M/jBjN4A6Vv93XNh5fIMCTC7/nzZ+XElLl5N59uXLocFLDXXVrMm/zLzww71vW78upNp7mLFnjGK4C9qnqSzHbzwdaqupj9Y60DkxiMAzDSFyyxjFcB/w6zvYNwHNAkyQGwzjYbVu3g4ev/jeLPl+C0+Vg4nljuexfF+BtUf1aw4mavWE9d82Zxbq9e8nwpvDbYcM5f+BhB0zffaNCIomhM7AxzvbN4X2GYTSy/L0FXDX8T+TvLUBDStAf5LMXZrPup008OOfOpH1oz/1lE1d89F75GIGdRYXcPecrikr9XDa02tV3jf1QIo3P24F4ozeGYPVMMgyjkX363Jf4in1oqKJK2O/zs37JRlZ+vyZp17l37jdRA8fA6t756IL5BEK1nXPV2F8kkhheAR4SkUki4gw/jgUeAF5umPAMw6jO6h/W4yuK069ehI3LNyftOmv37om73R8MsrekeaypYCRPIonhNuAb4FOgKPz4GPgW+EvyQzMMoyY9D8vG7Y3Tk0aVrn2TV8PbvVX8bq1Om51WdeiaajRvtU4MqupX1XOAXsC5wHlAb1U9W1X9DRWgYRhVO+6iCbi8LsRW0ZbgdDvo1r8LfYb1TNp1/jBqdKW++16Hg8uHDos7WMzYvyU8GbqqrlHVN1T1dVVNXiWmYRgJa9kmjYfm/p3Djx6IzW7D5XEy8dyx/OOzvyS1t9DoLt14+PgT6d6qNQBtvSlcP3IMV5iG5wNSteMYROQW4H5VrXHRUxEZDbRR1feTGF+NzDgGw7CoaqN0HW2s6xgNqz7jGHoCm0TkLeA9YIGqbg+f1AP0A8YA5wNtgQuTFrVhGLWSv7eA//7f23z1xjzcXhcnXnEsU688DntMFU8oFOKTf8/knYc/oiivmJFTh3LeLafTul38uX+qEpkUdm7axX/++gYLP1tCemZLzrh+KhPOHt1kiWNrfh4Pzp/LnE0baeP1csmQoUzt1ccksgTVOPJZRAYCVwFnYk27rYAfcGGto7EIeAp4XlXjTzvYgMwdg3Ew8xX7uHTw9ezctJtAqdWd1J3iZvgJQ/jLa7+PKvvA5U/y+Utf4yvyAeBw2mnVLp1nlt5HanriK5jt2baX3wz6PUW5RQQDVpdVT6qb0/9wEhfeflY931nidhYWMPnlF8j3+cpngvU6HFwyZCi/GzG60eNp7uq1HoOq/qSql2HdERwBnAZcBBwHtFPVoar6VFMkBcM42H3532/I2ba3PCkA+Ip8zP9gIRt/ruiuuvOX3cx4cXZ5UgAI+IPk5xTw0TNf1Onab9z7HsX5JeVJAaCk0Mfr/3yXwtwaa5+T7plFCyksLY2aHrw4EOCphd+T5/NVc6QRK5FeSSFV/VFV31XVV1X1c1U1A9sMown9OGspJYWVP/TEJqz8rqJvyJpF63G4Ktcc+4pL+XFm3RZlXDJrWVRCKuNwO1m/9Jc6nbM+5m35BX+cwXYuu51Ve8xHVSIS7pVkGEbz0bFHe5zuyh/4YhMyOlesQdA2qw2hYOUPTbvDRsdD2tfp2h26tyde1X3A5ycjq/HX7uqc1pJ4LQmlwRAdWrRo9Hj2ZyYxGMZ+bMolx2B3RDcy2+w20jNactiE/uXbeh3Rg06HdKhU1uFyMO2q4+t07TOun4orZnCd0+Wg/+g+dMhuV6dz1sdlRxyJO2ashdNm54iOnapcXMeIzyQGw9iPZXZuy/99/GfaZ2fi8rpwuh30Gd6Te2fdgS1icRsR4e7P/sLAo/ridDtwp7ho26kNd7xzI517darTtfsOP5Qbn7+K9IyWuFPcON1Ohk4+jNve/EOy3l5CBnfoyL2TJtPG48XrcOCy2xmXnc1jJ0xtknj2Z2ahHsM4AKgqOzftxuVx0rp9q2rL5u7Ooyi/mA7Z7ZLSjTMYDLJz425atE4lrXXTV9kEQyG25OeR7vaQ7jHTdVQlWesxGIbRQLat28FbD3zA+p820fvInpx67RQysmq3bvL2Tbu4ethN7NtpLWTfY3A3nlh0D3Pe+Y5Pn5tJKBjimF+NY9yZI1n69Qree+wT8vcWMvbU4daUGp7Kcy0tnrWMR697lp0bd5PVqyPXPnoJmV0zeffhj/hpzgq69M7i1OtOICOrDR8+NYP5Hy6ibac2TLv6ePoOPzSpP5tE2W02uqZXnxyN6iV0xyAiZwETgXbEVEOparX3ayLyLHAisFNVB8TZPx54F1gf3vS2qv61ppjMHYOxv1vx3WpumPhX/KV+gv4gDpcDt9fFQ3P/Ttc+WdUem7s7l9PbXRJ3nyfVXd5jyZPqpn23TLZv2FXeZdWd4qZL7048+M2dUcnhk+dmcu+vH690Pm+al0BpAL/Pj81uw+FykNY6lYK9hfiKSxERXF4XVz18MZMvMsvAN3f1GscQcZJ7gJeAbGAfsCfmUZPngck1lPlaVQ8LP2pMCoZxIHjwiqcoKSwh6A8CECgNUJRXxBO/f6HGY68ecXOV+yK7sZYU+ti4fHPUOAZfkY/Nq7Yy85U5Ucc9fNW/456vOL8Yv8+aLzMUDFFaXMqebXvxFVtDmFQVX5GPR699Dl+xGTewP0ukKukC4BxVfbMuF1LVr0Qkuy7HGsaByl/qZ+3iygsjqsLi2ctqPH7b+p31un5JoY8578xn8sXWN/yigmJKixMYqxqnwsFmE9b+uIF+I3vXKzaj6STSK8kG/NhQgYSNFJHFIvKxiPSvqpCIXCoiC0Rkwa5duxo4JMNoOHaHHWecgWcAKWk1N5za7PXrWCg2iZoryeVx1ut8AMFAiBbNoBHaqLtEfquewposr6EsArqp6mDgYeB/VRUMT8ExVFWHZmZmNmBIhtGwbDYbx1wwrtIHstvr4qTLj63x+BMvq7lMOQk/Irg8Tk6MuI7D4aBLFe0akWs+ANidtkqJyWa3kdWzQ41tI0bzVm1iEJGHyh5AOnCtiHwjIo9H7gvvrxdVzVPVgvDzjwCniGTU97yG0dxdcd90Bo/vj8vrIjU9BafHyciTj+TcP59W47FXP/xrMrJaV9p+2NEDSGvTgpSWXlJaekltmcJ1T1xK+66ZeFt4SG2Zgtvr4vL7ptP7yOgFfe776q+kZ6RFbcvo0pZJF4zD6XGSmm4d239kH869+RRc4W2eVDede3Xkb+/fVL8fiNHkalqP4cvankhVJ9R4MauN4YMqeiV1AHaoqorIMOBNrDuIartNmV5JxoFiy5ptbF2znW79OtOua2J3wqt/XM8Dlz5Ji1ap3Pa/P5CSkkLAH2D53FWEgiH6j+6N0+VEVVnx3RqK8oroO6IXKWneKs+56Isl/Dx3FYdPHFjeXrB7aw7rf9pE+26Z5XcF+XsLWPn9WtIz0uh5eHczxfV+orpeSY02wE1E/guMBzKAHVhrSDsBVPUJEbkKuAIIAMXA71X125rOaxKD0VysW7KRL1+dQygY4qgzRtF76CFJv0ZpqZ//3vU23763gNbtW3LxnefQa2j8JTzjxbPkq2W8ePsbFOYWMf6sUZxx/VRyd+Ux4z+z2bU5h0Fj+zLq5CMrTZ3R1Arzipj5yhw2Lt/MoUO6M+7MUXhS3E0d1n4tKYkhPA7hWlXNj9meCjysqhfXO9I6MInBaA5e+ftbvHLX2/hLA6gqLo+Tk387md/841dJu0ZRQTHndL6MorziqO2X3H0eZ904rcZ4uvbtzOqF66LKpbVpQaA0QDAQpLTEj7eFh869OnLfV39rNh+8W9du55qRN+MrLqWk0Icn1U2L1qk8Mv9u2nasXI1m1E5SxjFgrc4W777Ti9WV1TAOSlvXbuflO9/CV1xKKBhCQ4qvqJR3H/mEtYs3JO06913yRKWkAPDvP71MaUlFF9Nt63bEjSc2KQDk5xRQXFBCaYk1PqG4oISNy7fw9gMfJC3u+rr/sifJyykoH5dRUuhj7/Z9PHl9zeM8jLqpMTGISBsRaYvVn6F1+HXZIxNrNPOOhg7UMJqreR8sJN6Nt780wLfvfp+068z/aFHc7aow67Vvyl/Pfb9+d9ClJaV8/tLX9TpHsgQDQZbMXo6GNGZ7iLnvmZqChlKbAW67sYaxKLA8zn7Fai8wjIOSw+nAZqvc4GqzSZVjFOoi3jXKOCOmtHA4HfVuAE5m3PUiVPlemls7yIGkNlVJE7DmRxLgdODoiMcYoKuq3tVgERpGMzf6lGFx7xhsdhtHnTEyadeZeO7YuNttdhvjIq4z5tT48dSWO8XNlEuPqfsJkshutzPypKHYndFJwOl2cPS5Y5ooqgNfbdZ8nq2qs4DuwP/Cr8sec1V1a4NHaRjNWNuOrfn9M5fj8jjxpLpxp7hxeZxcft+FdDqkQ9Kuc+VDF9Ghe8wCOAI3vXRN1NoLbTrEjydekurWL4tW7dJJSfPg8jhxp7g54thBnHjppKTFXV/XPvEbOvZojzcco7eFh+wBXbnk7oYcb3twq2kcw1G1PZGqfpWUiBJkeiUZzcW+XbnMfW8BoWCIEScNbbAeM3Pemc+Xr35DRqfWnPeX02nZJi1uuXjx5Gzfy0t/fZO8nAKm/GYiQyYOwl/qZ94Hi9izNYf+o3pz6JAeDRJ3fYRCIRbOWMKWVdvoPrArg8b1M+Ml6qnO3VVFJITVhlD2P1BWOPY1qtokFX4mMRjJoqqs/H4Nqxetp0N2JkMmDcJur9+v9bb1O3j4t89QsLeQU66bwoSzxlBSVMIrf3+bHRt2MXraMI46fSSBQIC37/+QNT+sZ8CYPuXTVHz2wiwWzVhCt36dOePGk3G5nGxasYUls5eTntmS4ScMweV2svOX3Sz4dDGeFBfDTzyC1JYp5O3JZ/6HiwgGQww/YUjUnEiGUZ/EELlSyHDgX8BdwNzwtpHAzcCNqvphcsJNjEkMRjKUlpRy8wl/Z+V3a9CQYnPYSc9I4/6v/0ZGp7otbP/E9S/w1n3R3T5btE6lMLcoqpdNq3bp5O8tKJ92G8Cd4sLudFCUW1S+zeawMWzKEBbNWIIANocdp8vO0eeO4cOnv8Bms2GzCarKtGum8Pb9H2Jz2ECVUDDEbx+6mCmXNI+2A6PpJWuA20LgJlWdEbN9EvBPVT283pHWgUkMRjI8f9trvHHPu+X9+cFq1B08vj//nHFrwufLzcnn9IxGGvMpxJ3+OpbL4+SZpffTsUf7Bg/JaP6SNcCtH7A5zvYtQJ+6BGYYzcWnz86MSgpgLUbz01fLKcqvPKisJs/c8J9khVazWvZACgVDzHq9xllmDCOhxLAMuE1Eykc/h5/fGt5nGPutQGmgyn3BQLDKfVXxJbLYTSMJBUP4S5pfXEbzk0hiuAJrTMMWEZklIrOw7iCODu8zjP3WmNOG43BWbmju1q8LaXVYdGb6385KRlhJ5fQ4GXXysKYOw9gP1DoxqOr3WGMZbsJaVOeH8PPu4X2Gsd+66G/nkNG5LZ4W1qpp7vDaCDe+cFWdztfpkI4MOWZgpe3xVlxzuOL0fJL4ZTv36oQ3HKPT7cCd4mLw+P54Ut2IgN1hw+V1MeLEI3CnuLDZBLEJ7hQ3J146iZ6Hd6/T+zEOLo027XZDMY3PRrL4in3Mfn0uP89fRdahHTn2gvG0bBt/jEBtffbiLJ675b/4inyMOGko1z15KZuWb+GZm15i1+Ychk4axEV3nUPurjyevP5FNi7fTK+hPbn0nvOxO+w888eXWTZ3BVk9O3LZPRfQoXs75n2wkIUzFtOmQ2uOnT6ezM5tWTxrGXPe+Q5vCzfHnH8U3fp1YfWidXz56hyCgRDjzhxFvxG9kvSTMg4E9emueirwvqr6w8+rpKpv1y/MujGJwWgugsEgK+avIRQM0XfEoTicDlSV1YvWUZRXTO9hPfGmVr2O8/YNO9m6Zjtd+2aRkdW2ynKFuYWsXLCO1u1akj2ga0IDvRKJxziwVZcYapop602gA7Az/LwqCpgZrYyD1tJvVnD7Kf/E7wuUVwNdfu+FvHznW+zdmYvdJgQCIa56+GImX3R01LG+Yh93nf0AC2csxul2UlriZ9xZo7j+mSsqTRT3+r/e44VbX8PpdhAMBOnYoz1//+jmahNJmc2rt3HzlLvYtyMXWzXxGIapSjKMeirMLeScrpdTnF8SvUNAsAaclXGnuLj/q79FTTvx4JVP89nzX0Z1l3WnuDj7plM4/5bTy7ct+Gwxt596D74iX/k2m91G94FdeWLRPdXGGAqFOL/7b9m9eU90PF4X9331V3odkfzV5ozmrV7jGESkeSzjZBjN1Ndvf1dpvQAAFGK/ePlL/Lz32Kflr0OhUKWkAJQv9BPp7Qc/iEoKYHVB3bxyK7+s3FJtjMu/XUnBvoLK8fj8vP/4Z9Ueaxx8ajPpeq6IzAVmAl8C81S16k7fhnGQyc8pqHYcRKRQSNm7Y1/562AgiL+KY2MH1uXuzItbzu60k59TUO1183IK4rZFxMZjGFC77qpXY41uvgz4CtgrIp+IyB9F5EgRSWQshGEccA4/ekCl9QKq4k5xR40lcLqcdB/YtVI5ERg0tm/UtpFTh+LyOCuVDQVDHHJYdrXX7T+qd9zk5YmJxzCgdusxPK2q56tqZ6AvcAOwD7gOmAfkiMi7DRumYTRfPQ/vzphTR+BJrah19aS66dK7E+6Uim3uFBdZh3bgmPOjF9y59vFL8aS6sTusP0eH04E3zcvl910YVW7a1VNo06F1eXIQsRLNFfdPx+2tvsY3PaMl591yWqV4OvXswMTzzII3RrQ6Nz6LSEfgSuAaoIWZdts4mIVCIb5+az4f//sLQoEgx06fwIRzRrP4y2W8++gn5OcUcNQZI5l88dF4Uip/iG9evY237nufdT9tou/wnpx63Ym065JRqVxhbiEfPDmDeR8spG3H1pxy7Qn0H9W71nEu+nxJreIxDnzJml01AxiPNS3GBKAHsBCremmWqn5a9dENxySGA08oFGL35j2kpqeQmp7a1OHEpars2ZqDy+OKGgS3d2cuGgrRpkPDLNJjGMlSn3EMiMiDWIngUKxEMBu4FvhGVYuqO9YwEvX12/N56MqnKc4vJhQKMfyEI7jhud+Skuat+eBGsnzuSu6+4GH2bMkhFFL6jujFRXeezeO/e54NSzcBQlbPDtz8yrV0H9itqcM1jITVeMcQXsVtI9YiPR+p6vrGCKy2zB3DgWPFd6u5fsLtUTOTOt0OBo/vz/99fEsTRlZh95Y9XNTnWkoKo8cSoNadTqQWrVJ5acNjpLZMaewwDaNG9V2P4Sjg38CpwDIR2SgiL4jIRSJiZuQykub1e96jNGZaaL8vwJLZy9m5aVcTRRXto2e+IBiITgChYKhSUgAI+APMfs2sf2Dsf2rTK2mOqt6pqhOBVsCFwPrwv2WJ4vmGDdM4GGxdu514N7AOl5PdW3IaP6A4tq7Zjt/nr7kgUFLoY8fG5pHQDCMRCY1BUNVSVZ0F/B24DXgYK1n8KvmhGQebQeP6xV0TIVDqp2vfzk0QUWWDjuoX1S21Ot4WHnoP69nAERlG8tUqMYiIQ0RGi8gtIvIF1jiGmcAZWJPrXVjtCQyjFs68fiqeFp6odQg8KW5O/8NJtGjVPHonHX3eWFplpuNwVfTbcKe4aNm2RdTgM6fbSaeeHRh+wpCmCNMw6qU2jc+fAaOAFKwR0F+GHzNVdWODR1gD0/h8YNm+YSfP3/oaP3zxE+mZaZx5/clMPG9sQlNLN7S8Pfm8dOebfP3WPNxeNydeNokTLj2Gt+7/gE+fn0UoGGLieWM550+n4G3RfHpTGUakeo1jEJH/Ep4nSVXXNEB89WISg2EYRuLq1StJVc8JT4tRr6QgIs+KyE4RWVrFfhGRh0RkjYgsERFzD24kJC8nn8uH3MAk2xlMsp3BKW2nM/eDhXHL3vubxznOeRaTbGcw2X02T9/0n7jlVi5YyzWjbuY451mc0nY6z9/2KsFAsFbxlJaU8tjvnuPk9AuY7DqLG465g/XLNvHf/3ub09tdzHHOs7hi6I0s/WZFnd+zYTSERluPQUSOAgqAF1V1QJz9U7Am7JsCDAceVNXhNZ3X3DEYZU5vdzG5u/MrbX980T/peVhFz+q7f/UQX7z8daVy59x8ChffeW75682rtnLFETdGjVlwp7gYf9Zorv/3lTXGc/OUu1g8a1n5lNoiYHc6sNttUWM13CkuHvj6TrMes9Go6juOISlU9Suguj6HJ2MlDVXVeUCr8HxMhlGjuR8sjJsUAB69+t/lz0OhEDNfmRO33Bv/ej/q9Wv/+F/cdRJmvjKHvTtzq41n04otLJm9POp4VQiUBqKSAkBpsZ+X/lbdAomG0bia05TZWcAvEa83h7cZRo2Wzvm5yn2/rNxW/rwor6jSYjVlYqelXv3DekLBygPXXB4nW1Zvq7Q90qafN9d6Km5VZf1PTd6PwzDKNafEEK/bSdy/YBG5VEQWiMiCXbvMACIDBozpW+W+zr0qbjxTWqZU2cMpsgsqQM/DsqO6zpbx+/xk9exQbTxd+mQR9NeuLUIEsgdUXpPBMJpKc0oMm4EuEa87A1vjFVTVp1R1qKoOzczMbJTgjOZt5IlHkNamRdx9Vz50cflzm83GuLNGxi132u9OjHp95o3TKi2M4/a6OOqMkbRu36raeLr17cyAMX0qHe9w2nF5XVHbXF4X5//ldAyjuag2MYhIvojk1eaRhFjeAy4I904aAeSqavX364YR4dkVD5I9oOK7RWp6Cre/cyO9hvSIKvfnV37HxPPHIjbrzsHusHHqdSdwyf+dF1Wua58s/vn5bfQaeggiQkpLL9OumcIfnrmiVvHc8b8bmfzro/GkuhGbMHBsXx757m7OuWkaLVqngkCPwd2468ObOTQmRsNoStX2ShKRWo9oVtUXqr2QNR5iPJAB7MCaUsMZPvYJse7vHwEmA0XARapaY3cj0yvJiCcUCmGz1XxDXNtyqlqvQXbxjq/vOQ2jPuq8HkNNH/aJUNVzativwG+TdT3j4FabD/tEytX3Azze8SYpGM1Vc2pjMAzDMJqBWicGEXGJyB0iskpESkQkGPloyCANwzCMxpPIHcPfsGZRvRcIATcAjwJ7gJqHgRqGYRj7hUQSw5nA5ar6JBAE3lXVa7AakSc1RHCGYRhG40skMbQHloefF2At0APwCXBsMoMyDMMwmk4iiWET0Cn8fA1wXPj5SKA4mUEZhmEYTSeRxPAOMDH8/EHgDhFZDzwPPJPkuAzDMIwmUu04hkiq+qeI52+KyC/AaGCVqn7QEMEZhmEYja/WiSG8nsK3qhoAUNX5wPzwetBHhafVNgzDMPZziVQlfQm0ibM9PbzPMAzDOAAkkhiE+NNgtwUKkxOOYRiG0dRqrEoSkffCTxV4SUR8EbvtwADg2waIzTAMw2gCtWlj2BP+V4C9RHdNLQXmAE8nOS7DMAyjidSYGFT1IgAR2QD8S1VNtZFhGMYBrNZtDKp6h6oWishQETlLRFIBRCRVRGrdu8kwDMNo3hLprtoea5W1I7HaGw4F1gH3ASXAtQ0RoGEYhtG4EumVdD+wHasXUlHE9jcwcyUZhmEcMBKpApoITFTVvTErT60FuiY1KsMwDKPJJHLH4MXqhRQrE6sqyTAMwzgAJJIYvgKmR7xWEbEDfwS+SGZQhmEYRtNJpCrpRmC2iBwJuLFWcuuPNSXG6AaIzTAMw2gCiXRXXQ4MAuYCnwEerIbnw1V1bcOEZxiGYTS2hMYfqOo24NYGisUwDMNoBmq8YxCRFBF5VES2iMhOEXlFRDIaI7gDlQbWoEWvoiWfohqvPd8wDKPp1OaO4Q6sRueXsXofnQM8DpzRcGEdmFRDaN6foPhja4PYARe0eRFx9m7S2AzDMMrUJjGcCvxaVV8FEJGXgG9ExK6qwQaN7kBT8iEUf0J5714FKET3XgGZXxAzPsQwDKNJ1KbxuQvwddkLVf0OCACdGiqoA5UWCdcW+wAAEzJJREFU/ZfoyWnLduyBwKpGj8cwDCOe2iQGO5UHtgVIsOHaAKpsT7ARf+ygYRhG46vNh7tQeYEeD/C0iJTPmaSqU5Md3AHHexLkr6LyQHEnOPo2RUSGYRiV1CYxvBBn20vJDuRgIClnoyUfWtVGWgS4ABvS6j7MzOWGYTQXtV6ox6g/ETe0eQV8X6K+OWBvh3hPQ+wdmjo0wzCMcuZraiMTcYBnEuKZFLVd1QelP4J4wDkQkUSmsTIMw0ieRv30EZHJIrJSRNaIyE1x9k8XkV0i8mP4cUljxtdUQsUfojuHo/uuQPdeiO4aj/pXNnVYhmEcpBrtjiE8E+ujwCRgM/C9iLwXnoMp0muqelVjxdXUNLAWcv9EVIO0FqE5F0K7rxFxNllshmEcnBrzjmEYsEZV16k1D8SrwMmNeP1mSYvewOr9G8sHpd80djiGYRiNmhiygF8iXm8Ob4t1mogsEZE3RaRLvBOJyKUiskBEFuzatashYm08oT3ETwxAKLdRQzEMw4DGTQzx5nvQmNfvA9mqOgj4nPhdZVHVp1R1qKoOzczMTHKYjUvc44GUyjvUD65hjR2OYRhGoyaGzVjTa5TpDGyNLKCqe1S1bCDd08ARjRRb0/EcC85eWCunhokXUi5A7B2bLCzDMA5ejdld9XvgUBHpDmwBzgbOjSwgIh3Daz4ATAV+bsT4moSIE9q8hBa9DSUfgKQiqeeAa1xTh2YYxkGq0RKDqgZE5CrgU6z5l55V1WUi8ldggaq+B1wjIlOxKt1ziF5julnRUAFoHtjaY3W4ii8UzIXSueAahM1uzTuoGoLQdpA0xJaGiAtSzgLPeBA3Ymtdw7X3gZaEr21mZDUMI7lENbaaf/8ydOhQXbBgQaNdT7UYzb0FSj4FbFa1T9qfsaVUnioqtOskCEaOR0iHtFug8G4IFQIhcE+AlLMh71YI7rS2OQ+zpsmwt4++digH3fcHKP3OuratLZJ+N+Ie0YDv2DCMA5GILFTVoXH3mcSQmNDeq8E3C4ieU1BaP424h1eU2/Nr8H8de3gcTiAIhCK22cHeGcn4tHwEtKqie06GwBqiezF5kYz3EEe3Or4jwzAORtUlBjPvQgI0lAP/396dB7lZ33ccf38erXa1h836IGBswhFSN8YpGIgpuFwxITAwHC0JBpIJKYUGCuUoQwLNkEJDh4SE0jSEFDCYcBgIx3CE4UpIoaFAbCCAYwjG4TCGYBeMWe+p1bd//J71SlrtIWlX8mq/rxkNfn7P9f3JWN/n+B1dj5ObFAA6sU0/zS0aUVIA6CE3KQD0QmY9dP+2vyj9MqTfYmDT1h6s/ZYRnss554bniaEYve/DYD2Re9eM8skMMu/2L/auhYLjJ6Wh941RPrdzbiLzxFCMuh0YeHUPkID6/DuywV9Ij4hlIDm3fzk5N/RtGCAFyb0LlDvnXGk8MRRBaoTmM8jpcxC/gFbzabkbt5w5yFEicr/2RtAkwruGPiloOBDV7dJ/7sRMaDwi79x1EE1GTV8uvjLOOTcIH3a7SFHLKVjdLKztvyCzDpKfQ5PORnWfzNvudDIkoO1HhPcIEaSOQi1nYm3/Ad1PQTQVNf8dVr8fbLoKOh8Nw243HY+avjrg3Jp8KVY3F9pvAtsEDQtRyxkomlSZyjvnJgRvleSccxPQUK2S/I6hSJlMGj48BXqeIgz1lIJJ34KO+yD9XP+GyYNRy3HYxsug948QTYPm01HTCQM6pVn6LWzjxaEjnJKQOhpNOh9FzRWtm3POgd8xFC2z/khIvzLCrSNyX1Y3Qss/ErWcvLnEMhuwdYeEXtSbt62H5Fw0dan3bHbOjQnvxzBKMum1RSQFGNiCqQM2/QSz/r4I1n5XGN4iZ9tuSK+E9EtlROucc6XxxFCMnt+Vfwzriu8OYukV5Mzetpkg/Xr553POuSJ5YihGcrfyj6F60OT+5brPAKmB25lBYufyz+ecc0XyxFCEqG47SHy6mD3ylhuh+RtI/e/81XQsqIHceYzqITkbkn9RRrTOOVcaTwzFmnZPXk/jBmi+EBJ5dxPJA1Hr1ZDYMSxH02DSP6HmU3I2UzQFTfs51O9D+OtogMaj0JTr/cWzc64qvFVSGTKZDFGUm1szvb1EidzhMMxsRD/yI93OOefK5f0Y8lhmA9Z2dZhTQY1xT+MTB0y4k8mk4aNvQtfDQDo8RppyJWw4f3OLoQwJaDkPNl0Pti4uAxK7gKZC+tlwTgBtC00nwqYr+kpAM6D1TuhaDJ0PYUPEY5bB2u+AjpvCfA4NC1HL6Sgxbcy+K+fcxDPh7hjMOrD1R0Dve4ShKgBSkDqYqPWKnG0z674YOqdVRDIrnkZIHUTUemVuPB99GzruBzr694mmoekPoqilQnE652qB92PIYu33Q+96+n+EATqh81EsvXpzSab7+QomBfLi6YDOX2JZzVWtdy103Et/Uoj3yWwIdxHOOTdKJlxioOdpcn9cY0pAz4v9y52PVSykgpSAnqwObj0vDzIXRGdcJ+ecGx0TLzEkPgnUF1ghiGb0LyY/VamIBiGItu1fTGxL4bkg6uI6Oefc6JhwiUFNx4Wr8RyJ0Jy0/nP9RQ1HUziBVEJfPPP7i+o+C4lZDGgvoCRq+kpFo3PO1baJlxgSM9CU6yCaCTQQOpPthqbejLKmzoyiCKbdndtLmQhSJzHwa5tVRARbDyxquXDYeCShKUvimeLqgRRE26DWn6C6HYs4v3PODW3CtUrqY2aQWQukhm3umUmvht4PILnH5n4Lmc7lkH4RUscQ1bWGso7HoP0OaD6VKBVe9me6VsKmxdB4OFHjQaGs5yPoWAx184iaDio6Hst8AJl2SMz0fg/OuZJ4P4YCJEFi5uZlM4OuR8Jop/SixmMgdRhSgqhuZ6jLHbcoSu0J7Ll5OdP1NLRdCZl3oe0yMtF3wlAXH18Mva9B7woyiohSBxAlt4LkuUPGM2Ts0VSIppZcd+ecG8qETQz57KNvQucjQHtY7lkGnQ9B638Oe1Weab8TNl7YX5B+ET74m9yNej+GDaeQabmQqOWk0Q3eOedG0YR7x1CI9awISSBOCqGwA7qfhJ7lwx9g48UjP1nb98lkCrUucs65LYMnBghTapIeWG6dWNdTQ+6aybQBXUWcLA2Zt4qJzjnnKsoTA4SWRwU7j9WjaKthdi6hSWvUWvw+zjlXIZ4YAFKHFi5XBKnDh9w1iuohsdPIzxVtS+SJwTm3BfPEACiajKZcC9oK1BJ/JqHWH6PE9OEPMG0paEreQbeDaIe8skkw1cc1cs5t2bxVUkz18+ETT0HPC2AZqJ+HNLLHRFE0FbZ5hkxXvH/9AqL6MHFPpnsFdD8BdXOIUgeMZRWcc25UeGLIIiVzh8UoUtSwLzTsm1tWvyvU71puaM45VzEVfZQk6VBJr0paJelbBdY3SLo9Xv+MpB0rGZ9zzrkKJgaF6ciuAg4D5gDHS5qTt9nJwIdmtgvw78D3KhWfc865oJJ3DPOBVWa22sy6gduAo/K2OQq4Mf7zncBC+WBAzjlXUZVMDDOBt7OW18RlBbcxszTwETBgRDlJp0paJmnZunXrxihc55ybmCqZGApd+ecP7TqSbTCza8xsLzPba+utCwxj7ZxzrmSVTAxrgO2zlmcBawfbRlIdsBXwQUWic845B1S2uepvgU9L2gl4B1gEnJC3zX3A14D/BY4FfmXDTBixfPny9ZLeLCOu6cD6MvbfktRSXaC26lNLdYHaqs9ErcsOg62oWGIws7SkM4CHgQRwvZmtkHQJsMzM7gMWAzdJWkW4U1g0guOW9SxJ0rLBJqsYb2qpLlBb9amlukBt1cfrMlBFO7iZ2YPAg3llF2X9uRP4UiVjcs45l8vHSnLOOZfDEwNcU+0ARlEt1QVqqz61VBeorfp4XfJomHe7zjnnJhi/Y3DOOZfDE4NzzrkcEzYxSLpe0vuSXq52LOWStL2kxyWtlLRC0lnVjqlUklKSnpX0u7guF1c7pnJJSkh6XtID1Y6lXJLekPSSpBckLat2POWS1CrpTkmvxP9+9ql2TKWQNDv+O+n7bJR0dsnHm6jvGCTtD7QBPzOzudWOpxySZgAzzOw5SZOA5cDRZvb7KodWtHjQxGYza5OUBP4HOMvMnq5yaCWTdC6wFzDZzI6odjzlkPQGsJeZ1USHMEk3Ak+a2XUKM3M1mdmGasdVjngk63eAvc2spM6/E/aOwcyeoEaG2zCzd83sufjPHwMrGThA4bhgQVu8mIw/4/bqRdIs4HDgumrH4nJJmgzsT+hYi5l1j/ekEFsIvF5qUoAJnBhqVTy50TzgmepGUrr40csLwPvAo2Y2busCXAmcD2SqHcgoMeARScslnVrtYMq0M7AOuCF+1HedpOZqBzUKFgFLyzmAJ4YaIqkFuAs428w2VjueUplZr5ntThhocb6kcfmoT9IRwPtmtrzasYyiBWa2B2HCrX+IH8mOV3XAHsDVZjYP2AQMmFlyPIkfhx0J/Lyc43hiqBHx8/i7gFvM7O5qxzMa4tv6XwOHVjmUUi0Ajoyfy98GfF7SzdUNqTxmtjb+7/vAPYQJuMarNcCarDvSOwmJYjw7DHjOzP5UzkE8MdSA+IXtYmClmV1R7XjKIWlrSa3xnxuBg4FXqhtVaczsAjObZWY7Em7vf2VmX6lyWCWT1Bw3biB+5HIIMG5b9ZnZe8DbkmbHRQuBcddgI8/xlPkYCSo8iN6WRNJS4EBguqQ1wHfMbHF1oyrZAuCrwEvxs3mAC+NBC8ebGcCNccuKCLjDzMZ9M88asQ1wTzzbbh1wq5k9VN2QynYmcEv8CGY18PUqx1MySU3AF4C/L/tYE7W5qnPOucL8UZJzzrkcnhicc87l8MTgnHMuhycG55xzOTwxOOecy+GJwdU8SSdJaht+yy1TOfFLOkDSH+Lmv2NC0mclvVMjw0k4PDG4CpG0RJLFnx5JqyX9oJgfk/gYY9KnIR5O+ryxOHaV47gcuNTMekfxmDnM7CXgaeDcsTqHqyxPDK6SHiN0YNsZ+DZwOvCDqkZUwyTtC/w5ZY6bM0I3AKdJmrCdZmuJJwZXSV1m9p6ZvW1mtwK3AEf3rZQ0R9IvJH0cT6K0VNK28bp/Ab4GHJ5153FgvO4ySa9K6oivuL8vKTWagQ8VW7x+iaQHJJ0VP1b5UNINcW/Uvm2aJf1MUpukP0m6IN5nSbz+18AOwOV9dcyLYaGklyVtUpiYaadhwj4BeMzM2vOOc7ikZ+Lv6/8k3d/3fcXf30VxfT6W9Lak4xQmtLktjv01SYfknesRYCphNAE3znlicNXUQZhvoW+yoScIY+/MJ4yR1ALcJyki3FncQf9dxwzgqfg4m4C/BT5DuAtZBPzzaAU5gtj67AfMjdcfBxwDZM+m90PggLj888Bu8T59/powsNslWXXs0wBcQKjnPkAr8NNhQt8PyJllTdKhwL3Ao8CewEHAf5P7W3A28CxhQLk7gBuBW4EHgd3j7+Lm7ORrZt3AC3H93HhnZv7xz5h/gCXAA1nL84H1wO3x8iXAL/P2mUIY/39+oWMMca5vAKuylk8C2obZ5w3gvEHWjTS2t4G6rG2uJVyxQ0gk3cCirPXNwIfAkqHiiOM3YHZW2Ynx8aIh6rQB+Hpe2W+A24b5HpZmLbfE5/5RVtmOcdleefveDdxU7f/X/FP+x58Huko6NG5dU0e4U7iXMIgZhKvX/QdpffMpwhVsQZKOJVzl7kL4IUvEn9Ey0th+b2bprHVrgb2ztktmbYuZbdLI5xzvMrNX846dJNw5DDYTYSPQmVc2j5DEhvJiVoxtktqBl7LW9w3p/Im8/Tric7pxzhODq6QngFOBHmCtmfVkrYuAXwCFWuQMOra8pL8kzHVwMXAO4Sr5SEb3pfZIY+vJW2f0P6JRVlkp0nnLfccZ6nHwesKdTbEK1aMnb7nQuacS7jjcOOeJwVVSu5mtGmTdc8CXgTfzEka2bgbeCSwA3jGzf+0rkLRD2ZEWH9twVhF+XOcDf4TNwyTPBV7P2q5QHUv1PDCnQNlCwmOu0TaX8DjJjXP+8tltKa4CtgJul7S3pJ0lHSzpmr7JYQhXo3MlzZY0PZ617g/ATEknxvucRpispBTbSdo97zN9hLENyczagOuB78Wti+YA1xH+DWbfRbwB7CdpZnzucjwM/FVe2aXAlyR9N25ptaukc7JbT5VCYa7xmYTWSW6c88TgtggWpoxcAGSAh4AVhB/krvgD4Sp3JaGlzTrC/MP3EzpxXUl4Nv4F4KISwziHcEWd/Vk0wthG4jzgSeA+4PE43mXkvge4CNiecBexrsR69LkZ+DNJu/YVWJi86RjCFJDPE1okHUSoWzmOBx4xszfLPI7bAvhEPc5ViaQG4E3gcjP74Rid4zJgazM7eSyOH5+jAXgNON7MfjNW53GV43cMzlWIpHmSTpC0i6R5hP4Bk4Dbx/C0/was1hiOlUTolHepJ4Xa4XcMzlVInAyuBWYTWhm9QOizsLyqgTmXxxODc865HP4oyTnnXA5PDM4553J4YnDOOZfDE4Nzzrkcnhicc87l+H+NXfT8BwoLAAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Let's check out the clusters that KMeans found\n", "plt.scatter(d1, d2, c=km.labels_);\n", "\n", "# Add title, labels and legend\n", "plt.title('Iris Data: PL vs. PW Clustered', fontsize=16, fontweight='bold')\n", "plt.xlabel('Petal Length (cm)', fontsize=14);\n", "plt.ylabel('Petal Width (cm)', fontsize=14);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plot above, each data point is labeled with it's cluster assignment, that we learned from the data. \n", "\n", "In this case, since we do already know the species label of the data, we can see that it seems like this clustering analysis is doing pretty well! There are some discrepancies, but overall a K-means clustering approach is able to reconstruct a grouping of the datapoints, using only information from a couple of the features. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other Clustering Approaches\n", "\n", "Clustering is a general task, and there are many different algorithms that can be used to attempt to solve it. \n", "\n", "For example, below are printed some of the different clustering algorithms and approaches that are available in sklearn." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Clustering approaches in sklearn:\n", " AffinityPropagation\n", " AgglomerativeClustering\n", " Birch\n", " DBSCAN\n", " FeatureAgglomeration\n", " KMeans\n", " MeanShift\n", " MiniBatchKMeans\n", " OPTICS\n", " SpectralBiclustering\n", " SpectralClustering\n", " SpectralCoclustering\n" ] } ], "source": [ "print('Clustering approaches in sklearn:')\n", "for name in dir(cluster):\n", " if name[0].isupper():\n", " print(' ', name)" ] } ], "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }