ML0101EN-Clus-DBSCN-weather-py-v1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://www.bigdatauniversity.com\"><img src = \"https://ibm.box.com/shared/static/cw2c7r3o20w9zn8gkecaeyjhgw3xdgbj.png\" width = 400, align = \"center\"></a>\n",
    "# <center>Density-Based Clustering</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the traditional clustering techniques, such as k-means, hierarchical and fuzzy clustering, can be used to group data without supervision. \n",
    "\n",
    "However, when applied to tasks with arbitrary shape clusters, or clusters within cluster, the traditional techniques might be unable to achieve good results. That is, elements in the same cluster might not share enough similarity or the performance may be poor.\n",
    "Additionally, Density-based Clustering locates regions of high density that are separated from one another by regions of low density. Density, in this context, is defined as the number of points within a specified radius.\n",
    "\n",
    "\n",
    "\n",
    "In this section, the main focus will be manipulating the data and properties of DBSCAN and observing the resulting clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the following libraries:\n",
    "<ul>\n",
    "    <li> <b>numpy as np</b> </li>\n",
    "    <li> <b>DBSCAN</b> from <b>sklearn.cluster</b> </li>\n",
    "    <li> <b>make_blobs</b> from <b>sklearn.datasets.samples_generator</b> </li>\n",
    "    <li> <b>StandardScaler</b> from <b>sklearn.preprocessing</b> </li>\n",
    "    <li> <b>matplotlib.pyplot as plt</b> </li>\n",
    "</ul> <br>\n",
    "Remember <b> %matplotlib inline </b> to display plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notice: For visualization of map, you need basemap package.\n",
    "# if you dont have basemap install on your machine, you can use the following line to install it\n",
    "# !conda install -c conda-forge  basemap==1.1.0  matplotlib==2.2.2  -y\n",
    "# Notice: you maight have to refresh your page and re-run the notebook after installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\damaro\\AppData\\Roaming\\Python\\Python38\\site-packages\\sklearn\\utils\\deprecation.py:144: FutureWarning: The sklearn.datasets.samples_generator module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.datasets. Anything that cannot be imported from sklearn.datasets is now part of the private API.\n",
      "  warnings.warn(message, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np \n",
    "from sklearn.cluster import DBSCAN \n",
    "from sklearn.datasets.samples_generator import make_blobs \n",
    "from sklearn.preprocessing import StandardScaler \n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data generation\n",
    "The function below will generate the data points and requires these inputs:\n",
    "<ul>\n",
    "    <li> <b>centroidLocation</b>: Coordinates of the centroids that will generate the random data. </li>\n",
    "    <ul> <li> Example: input: [[4,3], [2,-1], [-1,4]] </li> </ul>\n",
    "    <li> <b>numSamples</b>: The number of data points we want generated, split over the number of centroids (# of centroids defined in centroidLocation) </li>\n",
    "    <ul> <li> Example: 1500 </li> </ul>\n",
    "    <li> <b>clusterDeviation</b>: The standard deviation between the clusters. The larger the number, the further the spacing. </li>\n",
    "    <ul> <li> Example: 0.5 </li> </ul>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createDataPoints(centroidLocation, numSamples, clusterDeviation):\n",
    "    # Create random data and store in feature matrix X and response vector y.\n",
    "    X, y = make_blobs(n_samples=numSamples, centers=centroidLocation, \n",
    "                                cluster_std=clusterDeviation)\n",
    "    \n",
    "    # Standardize features by removing the mean and scaling to unit variance\n",
    "    X = StandardScaler().fit_transform(X)\n",
    "    return X, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use <b>createDataPoints</b> with the <b>3 inputs</b> and store the output into variables <b>X</b> and <b>y</b>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = createDataPoints([[4,3], [2,-1], [-1,4]] , 1500, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling\n",
    "DBSCAN stands for Density-Based Spatial Clustering of Applications with Noise. This technique is one of the most common clustering algorithms  which works based on density of object.\n",
    "The whole idea is that if a particular point belongs to a cluster, it should be near to lots of other points in that cluster.\n",
    "\n",
    "It works based on two parameters: Epsilon and Minimum Points  \n",
    "__Epsilon__ determine a specified radius that if includes enough number of points within, we call it dense area  \n",
    "__minimumSamples__ determine the minimum number of data points we want in a neighborhood to define a cluster.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 1, ..., 1, 1, 1], dtype=int64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epsilon = 0.3\n",
    "minimumSamples = 7\n",
    "db = DBSCAN(eps=epsilon, min_samples=minimumSamples).fit(X)\n",
    "labels = db.labels_\n",
    "labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distinguish outliers\n",
    "Lets Replace all elements with 'True' in core_samples_mask that are in the cluster, 'False' if the points are outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True, ...,  True,  True,  True])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Firts, create an array of booleans using the labels from db.\n",
    "core_samples_mask = np.zeros_like(db.labels_, dtype=bool)\n",
    "core_samples_mask[db.core_sample_indices_] = True\n",
    "core_samples_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of clusters in labels, ignoring noise if present.\n",
    "n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)\n",
    "n_clusters_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{-1, 0, 1, 2}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Remove repetition in labels by turning it into a set.\n",
    "unique_labels = set(labels)\n",
    "unique_labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.61960784, 0.00392157, 0.25882353, 1.        ],\n",
       "       [0.99346405, 0.74771242, 0.43529412, 1.        ],\n",
       "       [0.74771242, 0.89803922, 0.62745098, 1.        ],\n",
       "       [0.36862745, 0.30980392, 0.63529412, 1.        ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create colors for the clusters.\n",
    "colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))\n",
    "colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABggklEQVR4nO39eXBk2X3fiX5+5y65YCmgUPvW1eyFZHdLXLrERXqSSJPU0BxZlD2STTsc0njs6CeHFS/CjFGYHk08hfXHC+vJjxp7JI1MaWRLfnqWPBxR5AybkkhatiSTHPbCbnZX9d5d+4YdyPUu5/f+ODcTCSCxFVAACjifDjQSeW/ee/JW5u+e81u+P1FVPB6Px7P3MTs9AI/H4/FsD97gezwezz7BG3yPx+PZJ3iD7/F4PPsEb/A9Ho9nnxDu9ABW49ChQ3r27NmdHobH4/HcMzzzzDMTqnq437ZdbfDPnj3L008/vdPD8Hg8nnsGEbm00jbv0vF4PJ59gjf4Ho/Hs0/wBt/j8Xj2Cd7gezwezz5hVwdt9xrtdpvz588zMTHBoUOHePTRRymVSjs9LI/Hs0/wBn+buHz5Mp/97GeZmZlBRFBVRkZG+PSnP82ZM2d2engej2cf4F0620C73eazn/0s1lrOnj3Lfffdx9mzZ7HW8tnPfpYkSXZ6iB6PZx/gDf42cP78eWZmZjh48OCi5w8ePMjMzAwvvvjiDo3M4/HsJ7zB3wYmJiYQkb7bRITJycltHpHH49mPeIO/DRw6dIiVGs2oKmNjY9s8Io/Hsx/xBn8bePTRRxkZGWFqamrR81NTU4yMjPDYY4/t0Mg8Hs9+whv8baBUKvHpT38aYwwXL17k0qVLXLx4EWMMn/70p4njeKeH6PF49gGym3vanjt3TveSeFqSJLz44otMTk4yNjbGY4895o29x+PZUkTkGVU912+bz8PfRuI45r3vfe9OD8Pj8exTvEvH4/F49gl+hr/HUJtB/Qaa1pFoAAaOI8b/M3s8ni0y+CLy28CPArdVdVnKiYh8CPgi8Fbx1B+q6i9uxbk9C2hrGr3yZ5A13d8AYQVOfxgpj+7o2Dwez86zVS6dfwt8fI19/kJV3138eGO/xajNnLFXhfLBhR9V9MqfuZm/x+PZ12yJwVfVPwem1tzRc/eo33Az+6i6+Pmo6p6v39iZcXk8nl3Ddjp3PygizwPXgf9eVc/320lEngCeALyK5AbQtL769vYccGVV3773/3s8e5vt+jY/C9ynqjUR+QTwR8BD/XZU1c8BnwOXh79N47ur9DOkwJYaV4kGWPFi5W24/TQq7vj9fPve/+/x7H22xeCr6lzP4ydF5NdF5JCqTmzH+XeSvoZUjPO1Fyb6To1r742EoARhGZIaYCFPIIjdwZuTMHwW4sGFF6cNN64Hfswdq9f/32cfP9P3eO59tuVbLCLHgFuqqiLyPlzsYFdLRG7GvdF9bXsObj8NQXXBkKqFyRcBgbFHC+NvoTGBvva/o2c+igydWvNcS28kgDPyjZuwKEBroHJosbEH59tvTS349rPmYmO/dJ+h0+t67x6PZ/eyVWmZ/x74EHBIRK4CvwBEAKr6G8BPAP9QRDKgCXxKd7Gmw2bcG4tem9ahcQvCKlQPAwI2cQcU3Gw8iGDuojPSeRsu/ylaPrTquZZl5LgnYfK8O/DIg+54Job2LLQmYOiUu7ksPdZS379aSGs9KwR1N70NXD+Px7M72RKDr6p/e43tvwr86lac627T15jCutwby16rFiSA9rSbKcfDYNuQtSAcgNY0tCbdPp3smqDSTaXsnGvpakNttnxGntTcbzHup3JoYVvzljPi8fCyMUs04MYOblxzby1eIdgcjno5CI9nL7DvHbPLjWl+5+6NTmpk57UmdLN8CUAsmABMFdKGuwnY1M3qTQB50/nhg3jRuTQaXL7ayFvQO+dWC60ZZ7DFQNp0O9oEJCxuOkUYJWuDpmAthDFaOexuYGEZZl4FidxqBiBPC+fbBeyBtyHN8SJeUHYjyds+m8fjuYfY199S2xiHt550RtmEaFhxBliCFV+zmntjuXuk+J/m7idrQTQImrlVgATFTSByxlWbzv3TeXl7Dm4+tXy10bgNtWtQPerGO3fR3UTS+oJLJh52x7epuwHkbagZyBsL5x6+z73/0x+Gg++AqVfACGRZMW6FymGo34JXfh81kXP11G+67QPH0KDks3k8nnuEfSueZhvjcP7fOsOZ1qA9A81xsOoMmtq+r+u4QNa1LSsMsG07g57OuaCqhG4mb1PnMrEpGOOMfdpYeH3eXFxMpRaSOVwsIHXjnLu4sFKwmbuZ5KlzFyU1SBpgW+5mY1sQDED1CJQPQdYAmxduKNzKZei+boUuCDQnoHYFZi8X57wBIm410p6B0oiv5vV47hH25QxfbQYXn3QGLB5a2JCnkM05w9qccIaxQ1IDm3ULmPq6MQaOu9lu2nAukuYEUMzgjYVoqDC8CZSGnZ+9dhM0AYwz2s0JSOYhjFzQtcNS/7oCc5ecC8cm7rkgLoLC6tw2WeYMs5SdS8kqmCKYK+JuJmohS9zNRYzL5qlfc/sEkRsTuO1TL7vHJnTjNSHMX4XyiLt2PpvHs0Xk7YTp8xdpT85RGhtm9NGzBCXfO2Kz7EuDT/1G142ziCAqfPAH3Ey5VahF5G2Xy145BOPPLcvasVnLpVq2pqEytuByyZPC1a4LNxYjzsDa1PnC86bzyasCxQy+dAAGT8K1/+L2E+PGrCz41zupPq1pd0OJKs5Pr7jjdPbJUyBx+6KQWpB5Z9jzxK1qghgkdsduTLibSuc8Wav4u+rGAj1xgtzFBJrjgKCjt5Aeg+8rdz13Qu3yLc7/yudJpmtuYqJKPDrIo//4Jxg8c3Snh3dPsy+/fZrWXXByRQyc/ghigoVc+hUKl/TQu+D1PywCqZ2XxzB4whnN8qhzvbRnCuNbFFylTZh4vjDWQ86thDjDm7edIU/rRWrnuDOu1cPu+FnTrThM4RbS3J1fAhaMfYd8yZtPoD0JyYxzIyHuNRPfgZM/DJe/tjCrdy8o4gxBsbqwxcogd49tCpmFPIPr38AeegwTln3lrueOyNsJ53/l82iuDJ491n2+PTXP+V/5POd+6QmCeLXvrmc19qUPX6IBNyM2YWGEe7AZRAOu+GnoNFIadj73foVLaQ1e+X3nFimNuJ9oAJJZmHjR+cjrNztndQY6GnC/NSty75PCwIrbFpbcrLo56baLuEBr3nLB07ThXD7RkFsJmLBw4ag71soCC71vsji/LSp0q5Dn8Ob/CeUxF/CtHnPumfKYO2SeLOTx2yKoixQ3m9Qds3YNXvw32Ov/F/rGl1x8wit3ejbA9PmLJNM1SgeHFj1fOjhEMl1j+sW3VnilZz3syxm+87VXnVFr3FyoVrWZm3Hf/193XQ+ripK1Z5whLhW+flWXHoksGGDVYlmaQjjobg5BqTDORc68iYvZeeHHxxYB2CK7R0zxXOL2DQfcc1mziBnUnXFdNrtfjWKlkbXc7H3mVXfO+k13zvY8hPFCfKAz6zcxC5Vjxc1FjHMLZU3n/7807p6PB93KqOMeCstQu4Ze/yYcOOtdPJ5ltCfn3PelHyJuu+eO2ZffNjEhnP6wczlUjxT+6NTNsM9+AlNdKFpaVZQsay+uXs0azugJzmha42bSEriAqW3jFlVaGFLrbLQUxjRLFwws2cJkvVuUrNAa77wJ50YBMGXQVuFmWfdVcMdNZooCsc64itl7ewqS0LmkTAgaQVR218rE7vU2cysdCQo3Dy4u0b3J4bKIRt/ekz5ad8Hv+cvexeNZRmlsuOfzvgRVt91zx+xLgw84I/PAj60dVOzNvOnVmk8bLhCbNpwxy1JnJDvVtSLObZS33QfYBGBKRapm7wc6Bw2LgGjPrHkRPc9ljWK/YhZkInezwgIdH75hme9+4Z33HE8XbhKdNEwt/PRYty1tuFTNdH7BfZMV7p2w7G4GNlu4aQQld9Oy+UIQvD0HjSLoHJRcYDse9uJsnmWMPnqWeHSQ9tT8IrdOe2qeeHSQ0cfu38HR3fvsSx9+BzEhMnQac/Adzl+/gj48Q6chqzkZhNaU+xGBUx92Br1+C1q3C8Nb5NV3DD9SLFGlm/feDXh2jLYt3DgrGumlFBk96MLxtPf1qx2nx9gvurl0jtkx+MVx8rbLwx864wq1qoddQDou3FM2LW5W1hl4zd3NI4gX4iPJ3EI6qQld8Rn45iyeZQSlmEf/8U8ggVC7eJPapVvULt5EAuHRf/wTPmC7Sfy0agWWq1EaZ1zHHkUGjqKVw/DaH7pZfKvJcuNpnI8/iJ2fvJMNs2y/TY0SyJa47rdak8668U++WPjjYzj2QzD9CkxdcIFfYwqDLguFYfFwsSJoOzdO3oYohOH7l4m4eXE2Ty+DZ45y7peeYPrFtxby8B+73xv7LcAb/D6sJqDG/GU48m6Yv+J80hIWQdg2i42tBZUFf7bEzjBmnUranRYLDYtUy/bqu2kxM09mXeFVNAjTr8KJH3Txj1vPuhWPCQuXTeSyhxQ38x867Yz8zCswcLKvYudq1cue/UkQRxx678M7PYw9hzf4BYuaiaQ1Z5jLY52NRaVtAlkTnb/iKkxt6gK3eT+jWfjKpfDdR9WF9MsNZdPcLXIXSN4Q4vz2aRNe+ncwWGQ7yYxT+exk76SNIraQw7HvQ0bfjrYmCy2hJXGQsLLQAczjWYHeyttwqIIgpPMNX4W7QbzBp4/7pj3r3DHhAKALevXgjPaV/wjVE4ubj1CkTi71kecJxFVnXLMWu8PYw4LPfgOUht1r6tfcTdBEzoAHpYXK26hKt5tX9ShIsDgrqlO9DBBWkNMf9gFbz6r0Vt6m9RYz5y+6tg+PniWqln0V7gbY10Fb6OO+KR9cmNnPvun0azqSBmGl0KYJ3ba8yE4BXKC0z4xZc5elkkw5F8+uYmVV0L4kc5C03czcZkXNgSlkKqQQgCs7bR0Tu+YvySx26mW3arr/E8ipH4KjjyOnfgh54Md8SqZnVXorb6unDlO/cptwqEI4UGb25cvkaUbj5hQv/PIfkCfp2gfc52xVx6vfBn4UuK2qj/XZLsC/BD4BNID/VlWf3Ypzb5qlGvbgMlDCsqtoNeFC45A8XdDfqV/HZcl05BL6Ie4mkdXdPgawK6VebjcCohsbSt4uKm4L/39z3Bn1brWyFOmYON9/Mg+3noN4cJG0gllFYM3r7+x+tlPYrFN5O3j2GI0bk+StBFOKqF26Rd5MSOYaxAcG0Czn+tee4fQnPnBXxrFX2Kpv0r/FdbT63RW2/1XgoeLn/cD/UvzecfpW0opxGSm3n3cFTVlPdkCWuDRM7czoV7OYnRuC9DzeLfTk4G/kNQBYp0VkE+em0hzywp2Vx+7ttucXJKA70hRZa9W8e6+/s/vpda+oVdrT85go4OF/8AlOfOS9W274eytvs0aLPEmZe+M6Wnze0rk6mucE5ZhX/9cnOfHRx302zypsiUtHVf8cmFpll08Cv6uObwEjIrIrInUrZoiEFZdvXjnsgorVI841k9YWZA/W4wNXnK+bQjphV8zuN0FQonuz6BRSCS4wbWI3848KX3+euFjI/CWXxtmRZO6Td9/fteb1d3YTve6VeHSQ2VevUL9ym9mXr/DUz/1rvv1z/5ra5Vtbes7eyltTimnenAJrCeIQMYagXAKFdLZO3kq91s4abJcP/yRwpefvq8VzO09vJW0vacNVhA6edHIB9Zu4mW2A89ev1x+f91SvbuR1u5GiSUtYdfn32IUYhubupighNG+754PYucPCyoLMgtr+q6qOa603iweWFWepzdD5K9ipl9H5K/5GsI103CvRcJXb37yAqhKPDFI+MoKJQlrjM5z/lc9vqS+9t/IWFLXqtAKtIsZg4nChcl24J7R28nbCxLOvcu2rTzPx7Kvk7WTtF20R2+Uc7Wfl+k51ReQJ4AmAM2fO3M0xufOtmkHyV9xA3/iiy9xZpEbZuVeu5dbp7NujUX/PYiGpFxXFgMbuZmjiQmgtLaqJC/2g0oEFIayOzELW7LuqWlWkrrPdu3x2lI57pTU+Q95KiEd6FGRFMFHYVbRcLYd+IzGATuXt+V/5PPUr4wTVEradoDkEA2VsO0XCgPLhEWya73qtnZ3W+t8ug38V6I3UnQKu99tRVT8HfA7g3Llz22IdV9PV0dY0tGYLvfve4RSVpcuPtni/RX5yYXWdm3sA7ak5yGpF8/UyHHjErWRm3iyanMMyESybufhIn7z7VUXqAILSisVwXo9ne+i4V7JGa/lGVcJKCZtkq86y78TgdSpv3/r8f+a13/4KjRtT5E03+QrKMaWDw7TGZxExNG5OMfHsq7syN383aP1vl0vnS8BPieMDwKyq7ioBlX66Omoz9PLXnPAXsqCN02Xp7F6WPO4x/h2FyWhgjeYr9xid4PXsWzB32f3dcbO0Z1ymTjGzB4ET37/IMHdcNNqeKzJ7aouP3ynOQtbl8vHcPTruFZu61aqqktaaNG9NobmlNDbEaoqWSw3e4H1HGTx7jDzJePqf/iaXn/zWii6OII64/7/5IYbuO4YAWa1BMlundXOa6fNv0Z6YJRwocel/+3Ne+ldf4OnP/OaWxxM2y27Q+t+qtMx/D3wIOCQiV4FfACIAVf0N4ElcSubruLTMv7cV573r1G+49MtO4FU7hVVL5qLSUazs6MQvnfmbhZl+lrhjLW28cs9SpGK2xiE+AGOPuiCtspClUz0GKMSDyOjCUn9ZwZvNoX7RtZIMitlZWIaDb4eZ15zxL430lWfwejx3n4575YVf/gMmn21Rv3LbFZOHhvKhEa5+5duMvffhrqLlUtdNnqTdFMsOyWyNyWdepT0xSzrfpHxweNUZvwImDqkcGyNvp+RpSjpdo3xkhJFH7sMErrZkN3bI2g1a/1ti8FX1b6+xXYF/tBXn2mpWy/vWtF5IAYvLRslaPYVWHcqFQFoRkDURTtCsN5hoC/eGcWme+V4yTVoY4MApakZVl9I6d9EZcJuCbUH50KKq2hX1iqIhFwc4/G533ScvwO3nilaPt5xM8/D97kbQg9fj2R4Gzxzl8f/X36d+5TbTL15EAkNQKSGqZI2ExrUJJr/zGqWDw7z8619c5LpJ5upIsFDsl7USrn/tWfI0gyAgGqoyePbYisZ6+vxFyC2nf/SDtG7PkDXbJHN1Zi9cxiYpMxcudY1+6eAQtYs314wnbJTN1CDsBq3/fe30XCvvW6IBNCz+Mbta70uyQiQvXPNRYdyWiqjR8/deCNz20lNfEA4sXJ/2tHNhBTEEVRh9J3Lig4t97EsL3tS6lNc8cdcwKMHEd905ygfdzD6tuZXR3FuuqYoYr8ezDfSbqZdGhnjb3/4rtG7P0JqYYe7VawSliLnXr/HcL/47GjcmGXvXA4tm8zOvXGb6+Tc58I7TZLUm17/2LPXrEwRRiE0ypr/7JpUjI5QODjH/xnXe+vx/pnJ4pGtYOzNkEwRUj4/Rmpzj1p9/l9bUPCLC1Hdep3F9giMfeIT4wOCyWXPeTpj4zmtMPfs6oIy992HG3vPgug32WvGHznVq3Jgkma1TOjBI5fjB7k1hN2j971uDv5oiZicIyMBxGDjhCq3ypMjSgQWXjRQFSGkxsxc3m91vaJGiaRO4+eLym2LwvuW9BnqzcrKWM+K9ekUXv+KUOQcKg9Ephut0zapdd6sJr8dzV+ln5DozdRO47Jip598gqJYIKyUkcO4220qZefky1VOHu26W4QdPMv3dN5l95Qr1y7fJ04wgCpEwIIgCgkrM7W9d4ND3vYOJZ1+lcXuagROHuob1+MceB1VsljP3+jWu/enTaJphkxRUaU5YJAy49Y3znPyRc4tmzbXLt3jun/0ut795AZtmqCpBKeLwBx7hPb/wU2tmyKwVcH3nP/pxXvq1P6J+dZyZC5ew7RRTihh5530MnD7cvSl0Mo5qF28uu2lsh+tp/35L+kkqgDMirSm3feC480nPXXFNQJbOzE0MUcll8Ky7gfheovN+i3aHt58rmqJ3snSKAO6V/4Q9/L2YHjdMNytH7WK9ot5D128Uuf2Fzz6suJl97RocOIscfIeXXriLrGTkemfqvSmaal0GT2tiFptbskab1u0ZqsedNpUJAkYePUtrYpb2xCwEAXmaEUYBA6cOE5Ri2lPz3Pyz51CrDJ09RvnwCK3xGebfukXjd/+UIIq49MVvUL90g6yZuHqQ4mOYN9rMv3WDaKJK5egoww+cYPSx+8nbCS/88h8w8exrBINlAqvO6FvLxNOv8MK/+A+871/8zKoGt1fioZfSwSHm3rzOs7/wb4hHh6hfHScarBIcjsmabepXx6kcPcgLv/wHPPB3P0o61+D+T33YFYt11D63Uet/335T+uZ9d2SQ27PoxAW4+W13U+ikHuaFqyEou3x0mxXtDTupinvU4IcD7qamWmQbUaxkDC5uERZGOXfXCnUVtxK4AGxadw1Ujp5bOGan4K054a5jx9h39Iqqh2Bmxrlx4h7fphiIBpAim8pz91jJyHVm6nOvXcdExkksTM1Rv+r6Lau1tCbnEGB62P27lg+PYMKAqFrmyI8/xo2vPUs0VGX6u29iKnHXrZK3U/JWm+qJMUwp4vrXniFvuZV1Wmti04xktkbWaC+W9+4suq0lna8z8cyrjD3+MK/+1pcJByrMv3kDm6QkEzNomndn12otM+cvLnMfLXXzrBZwTWZq5K2EaLBK3lyoTwgrJZKZGu2ZOaaee4PaxVuUDx3Y9tz7XvatwV+W9501F2SQs5bLOgnKTjkzrdPNvtEUNHKKkO3pBQMVDbm/95zRD5yhDytOAlnA/a/o22tiCEL3N63CrW8hPljILhRfktb0oqN2C95e/8PF/QRM6Fw3QcndMNrziw2+99lvGysZuc5M3aYJrYkGtcu3yFsJYoRosEI638A2E2yWM/X8WyQzdQShcmKMeHSQaKBC3kyoHh/j2IfezcRTL9OemidvJbQmZggrZQ6+50Emvv1yt5oXnJZOe2oeCQNMOcY2i89NJ3GuExDNldaNKV761T/CRO7zq4AmGUGlRFBeMObJfIPxp14BIwyeOrKiMS6NDWPznMaNSbJmm7BSchXGQUDeTgnKJbJme9n1UmDy2dcQMcSjgwze5465U1lE+9bgL5JUCMtF2T9FB6u46FIVwswbbn8TOqVLLYKUeZNuNyg1kM4543dP+/BDug3MwenjaAao08eJhtzjZH6hTy84qQUtZk1SzPq7vXwL+lTCSnkUPf0RuPSn7t/CxE5krePCGTjmVgpeQ39HWC2rJKqWefDvf5yXfu2LRIMVRITowAA2zWiPz6DWNbW37YS5164RxCH1q+OUxoaZfv5NbJ5TvzqOCQMqJ8ZoT86RtRPCaglTirj9Fy+i1lI+PIKqktWdsdfcfc6000lulSY+yWyNeHQITTLSmkvMKJcilMDNSxRsKwEjaJbTvO0mJWmzxQu//Ae87//zD7vGOBqqMv3iW+TNxN0wVAkqMSPvOEN57ACmyFZaer3yZhsUgmpIWCktXNu7lEW0Fvv2W7NIUqF2zc3ig5Iz7JXDrmG5dkTC8h7j3nFbBMWPOtmArkxyR2vnXqQ32GrcTL2TfaO5M8BB5GbrtetuuxjndumugHIIgoW6g7ThVkpjy1SzAZChU2jFBeaWdcMqjTgN/ea4l0veAdbKKgniiKhc4vD738Gtv3iBvJmQztW6OjfRyCDZfBMJDBIE2CQja7QojQwjBlrj06T1NrOvXMGUI4I4Yuy9D9G8OUU63yCttYgGK9SvT5A12mTzTTTLISg+p2t1bMst2XzT9WrI3SSmNTlPWI4IBytktSZ5OwURJp59DRMGSFFXEI8Mcu2vvo8zn/gAeTvh5V//ImPvepCZly4tuJhmG0w+/wY//P/9H3j1N7/sgtCV2N0UKs6HL4FBFYJyTPnIyOLxbVPufS/7+pvTkVTQ6990rpzKmMsMSWqFwc8WDBdS9IBloQArrLjKWZsVvWrFzVBb03T1Znac9ejv9+5TyD9UDrqnKgfhzI+4954XfvtLX12YXQnFFy9bOESWgmm4rJ2gDA//5KKA7aIzr9UNKyzD0GlfVLWN9KZhHv/Y41z7yrf7ZpXMvnIFta7YMBysOF/+nEBgkMBgs7zwkyvpfAOskqcp7ekammaEA2WCUkRWb6JpjoYhU8+9TunQAdK5BlmzzfT5i4gRbJoTVEtkdSdxImGA5msI5ynYJEWMQeIQTTJElbyVuBtB92Ov5PUW0eERglKEzXLak3O8+ptf5uRHH18Uy6ieOtStAQgrJfJGm7yVdLNvBk4dZubCJdqTs5hSRPXkIZrXJjnywUe62UoL49ue3Pte9rXBh8LgHDiLzl9e8BXHg26mn9QW/Nd5UsziLV1XxuAJZwSz4maguKyBsALpbjH464gpSLhwg5LQzeqrx9yMvj0DxiAmQOduwZU/c7NvVSCFdEkZfDjofPqH3g3Dp2DssRWNfff0q2gZebaXfmmY0XCVs3/zhzFRuCirZOrFt5h4+hVMHJDO1p3bJHOTIZtbtPgR01GzFGeokwyskrUS5+awLniqWYbEZTfTDgNskoOmSBA4F5FZSIfWNIPAdGfuK2K1m1tAIEgpcr5/AVMpYRvt7tjak3NEw1U0t9g0o3V7hukX31oUy+jUAHSv16VbtCfnOPTehzn3S08w/eJbNG9N056epzQySOngMG/+/76O5ou/h9uZe9+L/0bBYn9+VF3I+Z54ga50AMVsVqIFETATudfOvglJ7oKPNin271XIXILpFHBZdoWYWnc1Im710olDiOn28FUTweRLPb0AltJx6WRw4B2ARQ6/e91GW0zoZ/I7zGq55te/+syiAGPeTrj2lW9jSiEmCp1RznXBh62Kqqs+197nWj0ToTQnK2TJNbcuzdIY5t/qrCbsgr9e1QVd4xBrFc1AohDVdE3XjloLuVuJ2OZCdk/X2ANYiyok0/O4z7Ey+dzr3Piz73D8w+9ZV4VsEEfL/PF5O6E9Ncerv/Uk82/doDQ6hBjZ1tz7XvZ9T1twxkZOf9h9yFpT7idrwvAZGDxV2LIAgooL6JrQ/d2adH+PPuwCvyYqgo/B4qAmhoVLXeSsB6Ez/FG1Z7+dpCfQqta1L7SZ6wMg4UK3rk4DlO7Nquf1IoXPP/KCZruU1bTYNyLuNX3+Itl8k+MfejdZo00613Da9IFxM/k4hLDI8Mrt2v526/bL5hrOTVJrFh+z4vXqbgqmHBOWY0wUYszqQdtFxwZMOSIcLPesFHrQ3t8KxmCT3BVrfftl8jSlOT6z6CVrzdJrl2/x9Gd+k4v/4T8TDlVAIK01uO9v/jDnfumJHWm67mf4Bf3cClo5DK/8AeSNImMld4Ze86LzVeoCvh1phWPfB81J91zeZqHxSUdUDbo5ZNGQy0ppTrDrUjlN5Iz3/FVAXC79bD8lv94VTOH7Fyn63hovaLbLWEsaYCPiXp194wODDJx2zcU7WSpZrYkK0M4WPtlrhZIWbV95x3Sm5mb2yQYb3xhDWCm5LJ+1bhKmmPRgyOotvvM//jYH3/U2Zs5fYvjBk0RD1TUrZPuulu53N4kbX32Gkx99fGPj3yK8we9hqVtBADv2qDPg2OLLYN0MduTdbiVw4KwLTAYXiwBu7vze3YBvJ1Cjxee4EFiLh1zMIJl1Am27gmJJHhQxC9NyfnwxEFUW9unnhhJxq51OpyubeUGzXcR6tNhLY8NollO7Ok7zxiSgVI6PUT02tizA2JuyKcZgwoCwHGPz3E1pUutm+7b43ogsPO5ncHWFx33227CxB8ht4a5Zz76Lta+SWoO5164joaF5e5rD3/8IYTlm7PGHqRxdnm4Mq1fm7kQ6Zoc9Z/BXU7+8E2TgKDp8ppA0Tpwxi4pc8bzlyvvBBX3BBXyDCIIB58/XHr8mpuj3mkH5EIy93blNps5v7k1vCUXFbKdYqjwKR7/PyRIDVI7A7MWF+IV2gtedFcuAW/EEMa4Voi+O2k2sxwBFQ1UmvvMajRuTiIIKmMBQPjzCkQ880pUpmD5/keaNKbJWm+kLF13AVbWQP3ZNxRXFhCEdE9M10iv5wreDOz11Zknm6mCV+pVx5t64zpEPPMLEt19ZsWJ2N0gh92NPGfy11C/viIHjRQ/X3BmzPAFqLDNqvUHf4bNOHyYNQZOFYq7OjN6mrql3c9y1TtwyVgkUrwcp6gqyOiQGxh6B2tWF93XwEZg8D3SCXYXfP4gWUjYrR8AEvjhql7GWAWremOK1f/PHTlQsDgsvpJAnKY0bk8xfvMVLv/oFbv7n76KqmChk+oU3aU/OEw1VsElGOltfZFRtlriP5NJ0xN3GWu4mBdtO0SIWkc7Vuf2N89z/kz+M5tq3YnY3SCH3Y898I9ejfnknBkhMiB55L7z6vxVGraDIL+8ec2kueeUIEEBjHEwRtG2O4wJCkSv0yttsbYbOJow9tnh/FTdemzn9+YNvhxvfgvakkzo++LBzWWWJCzwH5SK2YeDIe5HBEz6lcheylgFqz9aoXbpFEIUMP3yarN4qqlvnyNsp499+ifFvv4SJIypHRkhmalSOHSQ6MFikOQr1q7edOyQUyIpzWYqstV3MWjN/Vay17n5phCB2ufo3/uw5zvz4D9C4Ms70i28x+ujZbv1CNFQlHKrsqBRyP/bOt3I96pfrENta6hLSymG4/aybtWshAWxiZ+BuP4sOnXJZPkuCvmgO174BwayrzM1TnHEPnKHvpHWqYXGF652wyZl9p32jFtk3laOueOytJ931NHEhItdyLQo7Xat83vw9w1pVs/HIYLfqVEQIyhH1K7fRPCdvp67RSRwRliPaU3PY3NIan2HogRM05+qota56db6xuY/iTtKbV7HkeTHiRNcKeWYJXXvH1u0ZEGHmpUu88e++uiggTmAQIJ2rLwuS71QXrq1qcfhx4F/iIpS/par/fMn2DwFfBDqpHn+oqr+4Fefu0Ff9csn2tTJG+rqEbOYM/cCx5S9YciPpBH2xGfrGlwr9nWShn23eovtt0OK4JnQdsEzsZkJSbOuy0nqzkxK5FVIORRaRMe5hnhRaQSxbLTH1Mow+7PPm7zE67QlX0mJvTcwSlCIoKlHn37rpqmMRNElRY5z3rhyRJykmCNAsJ6u30DQnrbWKKth1pkruRjrDDqQncOue1/bCd9I1cRc0t7QmZomqZa5++duUDx1YFhC3ecah97+T5vVJBs4c4cxf+yDx0M4lM2za4ItIAPwa8DHgKvCUiHxJVS8s2fUvVPVHN3u+FcexVP2yz/bVWNElNH8NmrcX67L3vq7fjaSz2gjLdD9Fqiwz3DYpZvjGpX1KgvvCUBh9Kap8U7oiZuBy3d1Bi9WGFFXB66nuNXQDrYvfSXEDEif9jDiJ4o5kdOfGZZN1r5Y8u4NeqYT7P/VhbDstcuqFsccfonJ0lMrRUQbvO0rt0i1qF29hswwxBtupZFXFthJSQDOLRq5tZ/P2NOlMzbUptPbend33kq9sSSQMivpEJW8k3PovL3Lg7aeIhqoMPbA4SUEC4ebXn6d+eZzyoQPMvXyFyade2RFZ5A5bMcN/H/C6qr4JICK/D3wSWGrw7y5Lq2U7rFdOdyWXUDwEjRvLddkL+t1IuqsNm7mMnrzVp7aqU5nauSEUGT3dNorFzSVrFa6fntepOv95OFh02wohBtpzuNn+0pP1foB7v5E9rqDODSZ33YMIq+68068s7mBlM7R+02vR3yMszb1Paw3m3rjO8NtOuKycb7/cneV/z8/9LWqXb1O/fBvNLXmSFkoiQlgtOdGwIrvYNhPn6qk3QUxXFZvCz72a0bxnECAOISl0ooy4qt1iBaNAe2KWmdwSDZZpT8512yvaPOf2ty6AsiFZ5M30zF0PW2HwTwJXev6+Cry/z34fFJHngevAf6+qW5qLuKYI1xr+5RVdQvEgG9Vl7642OlW50SH3epv0fCtgwZVSBD9N5J7LTTFbL9w5NgOJ6BrnqDDGQehcTY2brro3qrgbBPQE6IosGozrFauFTASwzPhL6Nw64YA79dwluquMDsk8TF1Aj7xn7Wu6xSmyno2xNPfe5jnXv/YM5Er96jgnPvY4JggWGaG3/99/FHLrtOy/+6ZzX3QqVaPQ6eNknQYi7jxBOSqe67So3APGPjIEUeRuXkHgbnYdscCOl7WjzyNCe3qecKDcba/Yuj1D3mhj85z25CyNSony4ZFV8/DXKozbCrbi29fPjbv0X/xZ4D5VrYnIJ4A/Ah7qezCRJ4AnAM6cObOxgWxChGtFl5CYjeuyd1YbNi/8+IWUgsQguXPjaOpuJCaAjoZNEHfz+xdJMpMX7paSc6UEJSgNudUDFAY/gihwy2oTLOj/QBEnSHBGPQDtTassCqnCAYicMiWlUReoTmuLVzx56txUEqzp1rkrKbKeDbE09751e8Z1ZBodJJmpddsPRsNVpp5/g1d/68tEI4OURocYeuAEYbXEzb94ARFI6y00s5hSSFC4NcSYYhWgSBwiSeYM/x5AgsClnyKoOleOphkYg4nDIoALmhc5+gr1y7fBGEoHhwhKEfWr40hgqF+ZoHFtkqAcc+SDj/TNw19PYdxWBHq3wuBfBXq/+adws/guqjrX8/hJEfl1ETmkqhNLD6aqnwM+B3Du3LkNTxXuOJi4mktog7rsi1YbpVFnHPM2iEI04gxy+ZCbPUgIyYybPdiiL65SxAsiunIOpRFnn8PYFUVlLZceCYXej7iOVO05dxwJF3zumi+sHrImi1cYQeESqjrN+s57Ko06mYjCYANFN6r7IWusGgS/Wymy+43NLu+X5t4v7ciUNdskszVuf/MC7ck58iSlcnSU6fNvIaFrUD5wcgwJQ9JaA82UgVOHqF2+RfPGFBK5j7TmObbVdnnqewRtp2AMmRGMFOqcSnflrKruOxsIJnCxDptkqFVu/cULrpm7CCMPnep22MqabW5/8wIHHjq1LA9/uypzt+Jb9xTwkIjcD1wDPgX8nd4dROQYcEtVVUTeh3MeT27BubeMrdZlX7TaaM+52fL4d90svnJoIQCcNlyVansGGrfpavWgRUl65Iy1CYuy9BSOPO7kintuPoA7V+0xuPWMWxFo5oy2Zu4m0ero9vTM7MG5fMQ4Tf+O2yoqQ3jC3WiWVhhnjdWD4FuUIruf2Yrl/dLc+3BJRyYTh9z+5gVUlXCgzNDZY1SPj2ECw+TzrzP6jvuwaY62M0qjQxz5wCMk8w0mv+Oqr6ORIZLpOfJGe+1c9nuNYvaOBVuKCjXQIi2zlXZjFWE5Jiv+jg8Oo0laBHbFdf+aqhENlQkHyoSVEq3bM1ibL8vD367K3E0bfFXNRORngT/B5Qj+tqqeF5GfKbb/BvATwD8UkQxoAp9SXakKZOfYal327mqjSHvW0YfRy1932jw2c0a8fBA58xG0fgNe/6Jz9WjuXDMd/37eBlspjL6rapUD71h+8xk6jQydRo++d+FGc+XPXDVvZ0YinfdSPO7cRJIGtGacUc9arpevFjeHysZ6ym5Fiux+ZquW90tz78tHRggqMcl0jaDiUoXzVkJQigkqCx2ZqqcO0xyfYeC+Ixx45AzTL76FiSOSmbprA6hKUCkVH528v7FfKaf9XkOBdpGWGgSYUrSo8Y9LvrMEcezUO6slRAzJXB0RoTU+TWtSMEHgGpgLnPrE+5f9+21XZe6WrKtV9UngySXP/UbP418FfnUrznW3uev55SILXwahe1eX0YfRgw87l0yWgK3hArZt97s15fzn1UNrppj23mis4G4kWRM6ssYSQlhyfn617uZj61C/viAId+ajrjL4+jecDHRYcbP7dQTBN5siu9/ZquV9v9z7gVOHmXvjOgOnDlO/fJu01iSoxBz5gOvIlMzWuP2tC7TGZ7FpzsDxMSpHRjn58e/DRCFzr10ja7WZeu4NmremVxYyu9cNfS8CoJggICxHgGDCgGSm5uSlFXcDNUL1xCHmXr0KQFCKXKPzOCKtuVaPo9/7Nkbeed+yU6xVGLdVlbnekbpNLPi1gcEed0bh15YHfgw5/VfQK/+x8KN3cuVlYSauFpqTaOXw+l1Lo29HRx92Eg+dSmEpCqxs6rKDwLlbBk+7G0OzqLIV47bnRbvC4x9ARt++9opnsymy+5ytXN4Pnjna7cTUiQUMP3yKuVevMv6tl7j6x9/m4LsewATBQiqhhWiw0nXx9DZAKY0Nc+2rT5O3EufyWEuHZg8gUYgJApellFnCwTKDZ49Rv3wbm2aYKKR8eIRwoEx7ag4TBa6lYjvDxBHRYIVosOKqctFFInS98ZnVCuO2qjLXG/ztYh1+bRk6jR553OnQqy20dgpfuzFuyRcNIs3xdfvAxYRw5iPOldSeXpwFJLKQSTR8vysws2mxmqi6GoThs67Xb9qAqVdg9O3rO+cmUmT3O1u9vO/XienQex9m9NGzzL12lXS2QengUDeLx5SiRU23e1cWww+eYPalywRRiIgUufk9Y12pKcm9fGMQg7UWkhQNDHEUcuxD76J8eITZ85e49ZcvkDfaaJajmaV8aARTCmnenEbTjGS27q5R4c5p3pxaMT6z9ObcaSe5Vfhv3jaxyK+t1gVxO8FQ1QW/dl4UWkVVd3PotEzsplqu7gPvm/teHoUHfxydfsUJodnMVc+2Z1zA+MDbnCqmKx4Aipm9AnMXnZHfYMDV96m9c7Zreb/U5VO/MUlWb1Eq0gcXNd0uVhZzQPnwCK3JOYyIC9j20vFvK11j1vmwSjed8d6y/Npy/So6CafNLOeVX/8/OPz+dzD84Em+5+f+Fpf/j29gkxxVZe7ly0QHqpz9yR/GttNuw/Os2Wb4gZNrxmfupk6+//ZtE12/dtZy0smLqldzOPrehf06RlHE5dx3e6g0IQhX9IGvlfsuY4+io293Rnj2EkxdgMGThXRCtuCGAXeD6bQqTGpQcrPKjQRcvd7OnbGW7s1Wzvh6XT7j37rA1a88xcF3P7DY2EN3ZdGenMOUIlQteau94oxeShFiBJukhS4+mHKMDLhVgW3ulqY/G8AIQTmmdOgAeSth8juvM3jmGFPPv8EP/s5nmHv1Ks0bU7zxe18nPjhI+eDCSqw9NU9YLYGwo41RvMFfB1tSMTpw3AVdZ151mTad6tU8LZJULzhf+8BxN7Nvz7htQbSwH+qyZ/r4wNeb+94VeBs4jjZvuxuQ7fnyqTq3j+nJ9+7Z7gOu20M/3/tWL+87dFw+zsVzrevi6dC7sph49lUaVycQdWJ74UBE1uk/C3T62Wrq2huaKMRmLs3Y2BATRgRDA7Ta6aK89l2PCBIGxKODrhymHJPMN2iOTxMPVZl79WrXUI9+z/0r3qxnX7myLemXK+EN/hpsVcWomBA9+A7nBzcCWTHDNyEceLDb9FuGTjuf+1uJc6dkhStIQjhwFjn9V/rfbDaY+77Iz541i3gB7gZjRorq4MK4mNgHXHeAfr73u3q+dawsBEECIR4ZIL/ZJu/ozHRRtJ0iUeAklUuxc4nkFptbDn3P27A2hzeMC2JmK6R17ga62XRCUCkap/esfEQhmasTDw8sMtSr3axbE7Pbkn65Et7gr8KWV4xK4IxupwLWxIVWj4Gs2XWXSHkU3v430fmrMH8FEBg6hRTa+33Hege57x0/u85fgSv/0d1UKodcbGHuLaebA26sQeQDrvuAtVYW6XyDkUfPMvvyZYgCtN6ZKDjNGRMG2FYCxhCUSoSDJafEef8xbJbzwN/9KK1b0zz3i//OuYJ2rbEXMOLcV0UwlljIkwy1FjEGFSUeHuhrqFe6WW9XfGYl/Ld3Nba4YlSiAVRM1x/eb3v3sQmRA2ddTvx6j73G9r7PmxA5cD9a+hvuJtaRa6gecZk8Bx9BBo75gOs+YrWVRWlsmKha5ugPfg+137sNKBhBxDglSQIomqUE1ZjD73snI4/chwkCps9f5PIf/iWzL19xBVu7SYqh03C98yUybnUjRlwGs7XQTMhbiZONTnNMYAjL8YYM9XbGZ/rhv8GrsOUVo3czP32Tx/ZZNZ71MProWcKhCrf+8gVKYwewaYamuTP2hWa+hAF5O0ECYeDkIUwQkLUSJp55ldLoIEmt0TWkO07nC5zbwugX9S9FtpFtpwRxRNZsu2pbq2iWI4GhfPgA0y9d4sP/7L/dkKHezvjMUvy3eRW2umL0buanb8WxfVaNZy2CUszJv/o+bnz9O5jAEJQiNArJak2kyM0HN1EuHx5l/KmXOfR97+DanzxNMj2P1ZzWrWm3Q2gg26FZfrcuQJw7KnM3rK5/PTLEw1VQJRoZQm9NE48MMvTgCYwxlA+PEB8YIG+0XQvDDbLd8ZkO3uCvxl2Ykd/NmbSfpXvuBkurQkXh0Lm3I3HA7W+cxyYZaSlCAkPeSgmrJTTLCUoRzdszXPnytwAoHz5Aa2JuIY0zW2E6tR1FWt0aAcHEEVYTl1GUZGCVeKhKWK2Q1Ztk801MFFAaHWTk4dNUj491D1O7dOuuZ9ZsJd4SrMLdmpHfzZm0n6V7tpJ+qp15kmGTlNG33U/8X1W5+sdPYbPcNc/Mc0wYcPj/9j2Eg2Vmzl8ClAPvOMPVL3+LdX0wpbD4W230O0YeFm461mKTFEGc8Y9CbDslbyWuR6+qqx8AGjenaE3MLDL425FZs5V4g78Gftbs2a+spNrZnJjh9jfOUzl6kPLhAxz9gce4/vVnSWZrSGAIh6vMvnyZoOxSGePRQSQInA1fT4OUrW6C3mlKvlLrxcyiRhBrUavdVoZBJaZ66jDxgQFqb9100tDPvk7p4DCVYwdJZxvbklmzlXirtQ78rNmzH1lJtbNyaISh+48x9+Y1Zl+9Uqhs1pHAMPzgyW6TlmS6RjJTY+DMEWw7oXT4AHOvX+93qruL4gKyWsg75IVUSZYttIcQQeIIAfJ6CxHD8EOnCCsl8rbr36tpRmtilqt/8hTxUJXDH3iE9/zCT21LsHWr8Abf4/H0ZSXVzmS2xuzLV6icGHNNPcZniYYqlA4OkTcT8lbqNPOrMUMPHgecGyhvtDGhKSpvtwEpfpzDBjXuibBaxpRj0pn6QicrcIqYcUh5dIhkrkHtzRuUjx+kdXMKEwSMPnY/yUydgfuOUB4bIRooUzl2cMXT70a8wfd49jCbaZPYT7WzK6GsMPy2E1SPjzH3xjXS2Tomjjj4vW8jb6eElRLlIyM0rk5w8N0P8Oq/+WNaE7PYdAX9/A692jxRUOTqK/RL5umX5dMb8O3+tu6hEcQIQaVE3nQ3nzyThewcI1QOHWDksfuZuXCJvJUQHxjAtlIqx0Zd/n2Sdd/3dmjfbDXe4Hs8e5TNtknsVxXauj1DMlMnPjDQlU8Oq2VM7IKdYgzDbzvRPYbNc8a//RJHPvAI0UCZyeded0VXK83yrTptKQsigoogEjg9q6X0lShYkuJjDBIIEroMHFMtkdVaiAEJDBIF3bx6cgtGqL11AwkNJgwIyjHhQBkxhqzZXiQbvR3aN1uN2YqDiMjHReQVEXldRD7TZ7uIyL8qtn9XRN67Fef1eDz9WRpwHbzvKINnj6G5cv5XPk+epGseo1MVKoFQu3iT2qVbzF+8iRhZJJ9cPjxCUI7Jk8w1Si9oT82jmSVvJkw99zpZrUlQKRFEYY+7pQ/a8yC3aLbCqqAIsAJIFLjjdSyacdIIQSWmdGiE4bedJBooE8QRNk3J26m7BqquYfuREYJShImjngWCupTT+QbJTA2Rxe/7XsvQgS2Y4YtIAPwa8DHgKvCUiHxJVS/07PZXgYeKn/cD/0vx2+Px3AW2qk3i0qrQ5vgM17/yFPGBwe4+Jgw48sFHuPbHT1G7dJPmrWmCUsTgfUc5+bHHeeFf/AHhYIV4dIgD1RKzr10jz7KVZ/mFxdXMQlj0gejn0hHBVCJILeVjB5HQkNfbtGdqrhuXKhIElIYHCAdijn/kPVSOjjLx9Kuue1WeU782Qd5ou/z7ovF4aXCI1nyTg9/7AO/8R5/k/P/0eUwUM/zQia6x3yrtm8243O6ErXDpvA94XVXfBBCR3wc+CfQa/E8Cv1s0Lv+WiIyIyHFVvbEF5/d4PEvYyjaJvVWheTth8qlXlol/JTN1pBQRDVZdUFadS6Y9X8e2U8LDI+5YpZjB+44y/8Z1snqLoBRhc+vSNZdq61glrMbF8ZzapgkDCAMMgrWWeHiQtNbAhAHHP/Qubn/rJbJ6C40UAcpHRzj4ngcZettxGlfGOfIDj5HX28Sjg1z942+T11pIaMhqTdSqc98UfadP/7UPcOwHv5fB+45y/lc+T+PK+Lq0b9ZrxDfrcrsTtsLgnwSu9Px9leWz9377nASWGXwReQJ4AuDMmTNbMDyPZ/+x1W0SO/QT/9IsZ+r8Wxx6z4MM3teTrz8+w+U//AvSdkL9yjhBNSaII8JKjBgnsxyNDDqpBbvE2BfB26zhDHJYKSNJSnxwCE1zSkdGsO2U4QdPkicppdEhbn/zgsufL6p+B88eA4X5N68Tjw4Sjw5y5q99kPFvXuDGf3oOE4SunWMpQoqCqyAOOfR9bydvp91m4xvRvlmvEV+pxqG389XdSPfcCoPfbxqx9JO2nn3ck6qfAz4HcO7cuV2hr+Tx3GvcTRnefm4eFRYZ+2S2xuQzrzJ/+RbZXIN2niMKQaVUyBQMObeLtQguOEsgXSNnopC01iQ6UCUarBINlIvesBS9mJXS6BBDZ49j4oD7/psf5IVf+gPi0UE0z5l77brrtwtktSbJVI1z//wfEA8NdLWAJAzQPCetZQTlmOGHTpI3E9Jak+qxg4uu0Xq0bzZixLfK5bZRtsLgXwV6NYJPAUurK9azj8fj2SLutgxvrwG89tWnFzUGsVnO7W9ewOY5muWYOCKuDGCzHLVK+fAoBMLAQIX4QNUZ8ty1D+yoVAZVV/AUVsocevwhGtcnAWjenEJzi1rLwKnDmDjodpIqHzrA4H1uFj30thO0bs+QNdskM3Ue+Lsf6c6wgyjk0Lm3E1RiWhOzzL16FQmEvOnkFGw7vaNrtBEjvpUut42wFQb/KeAhEbkfuAZ8Cvg7S/b5EvCzhX///cCs9997PHeX7ZLhXeo+ao3PON34wKC5ZeDUIdrT8y7N0ubYLMMQ8D0/97e49tWnmX/9OsnMfJHWKd00yIFTh7GZxSY5ww+fIpmuMfzwKY59+N0cePAklWIW3q+TlAmCruZN7eLNRQVSpbFhxAjV42NUj48x8sh9CzeH6Rrf808+dUc+9I0Y8bvlcluLTRt8Vc1E5GeBP8G12/5tVT0vIj9TbP8N4EngE8DrQAP4e5s9r8fjWZvtkOFd6j7KGi0AsnoLCQzx6BDx6BBZvUUyU2P4wZOUDg1TPXmI7/+Nf8x/+QefZfzplxARgmqZIAoRE2CzlCPf/xgP/tTHSOcaq96wNuLCWrpv5+bQnpqnPDbM2HseuqPrsBEjvlOdr0R3cRPhc+fO6dNPP73Tw/B4PGvQG6xsTc0ze+EiQbWE5pZykaEDkMzUOPz+d2LbKe/8f/x1Dr33YWqXb/HcP/td5wZKM1SVoBR1tWrWO9veSNbL3ciQydsJT3/mN9FclxlxCWRZIPZuZemIyDOqeq7vNm/wPR7PVpAnKdMvvkXzxhRv/N7XCQ9UmP7OG6iqEyFrJmBg7PGHCaJwkQHMk5TJ77zG5DOvAcrY4w8z9p6HNux+6h1De7ZGPDLoArB9UiM7+26lu2ujRvxujMEbfI/Hs610DF/9yjgzL13CtlNMKWLkkfsYOHX4ruaa70R+ey93w4hvBG/wPR7PttOdbd+coj1TozQ6ROXo6F01gBt1q+xFVjP4XjzN4/HcFXaib+tO5bffK2yJeJrH4/HsBnYqv/1ewRt8j8ezZ9ip/PZ7BW/wPR7PnqE3v72Xu53ffq/gDb7H49kz9NPwr128iQSyJZIS9zo+aOvxePYU2yUpcS/iDb7H49lz7ESG0L2Ad+l4PB7PPsEbfI/H49kneIPv8Xg8+wRv8D0ej2ef4A2+x+Px7BO8wfd4PJ59gjf4Ho/Hs0/YVB6+iBwE/gA4C1wE/qaqTvfZ7yIwD+RAtpJ0p8fj8XjuHpud4X8G+LqqPgR8vfh7JT6squ/2xt7j8Xh2hs0a/E8Cv1M8/h3gxzd5PI/H4/HcJTZr8I+q6g2A4veRFfZT4E9F5BkReWK1A4rIEyLytIg8PT4+vsnheTwej6fDmj58EfkacKzPpp/fwHl+QFWvi8gR4Ksi8rKq/nm/HVX1c8DnwLU43MA5PB6Px7MKaxp8Vf3oSttE5JaIHFfVGyJyHLi9wjGuF79vi8gXgPcBfQ2+x+PxeO4Om3XpfAn46eLxTwNfXLqDiAyIyFDnMfAjwIubPK/H4/F4NshmDf4/Bz4mIq8BHyv+RkROiMiTxT5Hgb8UkeeBbwNfVtU/3uR5PR6Px7NBNpWHr6qTwEf6PH8d+ETx+E3gXZs5j8fj8Xg2j6+09Xg8nn2CN/gej8ezT/AG3+PxePYJvqetx3OHWM2ppzOktkVkygxEIxgJdnpYHs+KeIPv8dwBrazGlfoFMpt0nwtNzOmBRyiHgzs4Mo9nZbxLx+PZAFZzZpPbvDb7bdp5k5KpUg4GKQeDqCpX6hewmu/0MD2evvgZvsezTjqz+lY2Tz2bISCindUYjg8TmpjIlGjlNerpDEPx2E4P1+NZhp/hezzrwGrOlfoFVJVAIgKJCE2Moswl46ja7r6pbe/gSD2elfEG3+NZB/V0hswmRKZEIAsL40BCrOYkttV9LjKlnRiix7Mm3qXj8ayDdt4gzduoWowECIZcs67xt5qR2jahiRmIRrb8/D4jyLMVeIPv8axBK6txq/kmjXyGwEaAYtWiasklwaJkmhJJmdMDjywzxEuNdSUcopnNr9t4+4wgz1bhDb7Hswod330oMZGUsZqTarubiWM0wohwvPIgY5VTywz3UmOd25RmPk8lGCIwEbC68baac7n2IoltYxCMhMSmTKYpV+oXeGD4cT/T96wb78P3eFah47uPgwpD8RiJNsk1RzCoKqEJGY6PMpPeWvba3kBvORikZKq08waqSts2KAXVVdM5reZcr73CePMKtWSSWjLNfDLBdPsGgpDZhHo6s2j/2eQ2V2svc7X2ErPt2z5F1LMIP8P3eFagY0DbWQPBgFoiKWOMwVqLlZyBcJRKONg3HbNzsygHbuae2BZWc0ITk9gm88kU5aBKFJRp541Fr29lNa7UzjPRvIzVBAixklGSgW5mUCUc7mYEtbIaF+efZzYZR3FGXggYjg9z/9C7vOvHA3iD7/H0pTfnvmXnSZMWlhwQl5YZQGaTrlsGlqdjpj2ZO+ACu5acdtbAakZuM5rZHCJCKajSzpsMsbAySPImxgQYG2IkQNXStnXKZohcUzJ1WUNWc67UzjOfThJISCBlAHLNqKWTXK6d58ED57zrx+MNvsfTi9Wc+WSSK7UXEQzVcITEtlBVRC2JbRFJCUuOkYDYOOOqqqS2xXjjEvVsFtSCuJUAQefYllZeK86kZNom0wSDIcmb3DRvMBiNkNo2mU0QMQgBSkJmU4w4N5LF3TgEYSAaoZ7OdI+7NGU0swmtvMZ8MomRwGf57HO8wfd4Cjqz+mY2TyObIZCIlq0zEI5Qz2YQBNWcZj5HaGIGw4MANLMa9WyatNFiLplAyRERQonJNceSUw2GmE8mULQ4myIIoFgsAUHXlz8aH++OKbFNAGfkFRRo5y1CE3Gi+hAAs8k47axBTkYgUXHcBbK8zZXaiwQm7j53N7N8fArp7mVTQVsR+UkROS8iVkTOrbLfx0XkFRF5XUQ+s5lzejx3g44bxdocay2qAIKqpZ7NcCA+QmyqgCHXnNxmzGcT3G5eYj6doBocoJZOISJEpkxAhKKExMwkN7nZeItEmziT7Yy+Yosf5yrKbKubrgnq3D0IgiEyZUKJMRgEGI2PUQkP8MbcM0y1r9HOGyR5k3o6TWrbqCqqSq4pjWyezKbbovvTymq8MfcMV+svcav5FlfrL/HG3DO0straL/bcdTY7w38R+BvAv15pBxEJgF/D9by9CjwlIl9S1QubPLfHsyVYzRlvXmK+PUlWpFxmto3VzJlbCZhqXSuqaS2xKWMkoBIeILcJbdvoullCiZ2xRcltThyEaJ4TSAyUyLRdzPI7M32Dm7+3aeZzaKZkeUIYxIWPvkxiG6iCKogIqFIOB7lavwAIZTPAvEygmmMRmvk8gUTFOHIUpZnPk7SajJSOEZnyXdH9WZqV1CG1bZ9CukvYbE/bl6D4EK7M+4DXi962iMjvA58EvMH37DgdN049maaeTyFqAIOIYAt9nFTbRKYMqgQmIjZVLBnNbJZQSuQ2pZHOFTPqDEtWzNxdha6SY9WSk6LYJSNY+Du3GYrSyGbIsxxQTBEAsFhAiSghJmC8eZFMM0bjY8ynk8SmAhYyTVFyMs3d+yj+s5qT25Tx1iUOl84SBU7+YSt1f5ZmJXXwonK7h+3w4Z8ErvT8fRV4/0o7i8gTwBMAZ86cubsj8+xremekoYmgmM2nmrggLQCd2XrmsmnMAJacVl530go0yUgRbRfe+GzROToGPqGx5nhysuJ8pvs6C0RSQciKG0oOWGzm3EETzSuoKIEExKZCTJl23uyuIkITY9UFmJGATFNm0pscCk4DEJiQ+WRymb/9TvzwS7OSlm/3onI7zZoGX0S+Bhzrs+nnVfWL6zhHv+m/9nnObVD9HPA5gHPnzq24n8ezEfoZsN4ZaWpboEJKUsyJwRCSkwMZIobYlLozcC1MakZH7mClcNhGPsId3772PKOk2sQQYMndjUAhX/KyTA0ghBI5rR+BzGYYQncjUuvel0KSNZlrjxObClfmXRWwMYGrJg5KHCmf5Xbr4oalHKIiY2nl7V5UbqdZ0+Cr6kc3eY6rwOmev08B1zd5TI9n3aykRTMULrgXnNxxRGLzwh2jKCmGACUkFJdv38zmsFhXiNV1x0i32GlrWDD4gqBY8kWuIGHpjaRzk0i1jSEkkhhEMWIQcS6VzvsSDPVslnmmCCV2781CKi1Uh3l19lsMx4ddcBdLmrdopvO8OfcsD498gLAn26eXgWiE0MSktr3IuN9NUTnPxtgOl85TwEMicj9wDfgU8He24bwez6qBxKn2NToL0NiUEVk8SzeYnjm308FfcLv0Gvi7txBd7vNf6XwLgWBLRmaVnLyIIWiPsRdiUwIMuc3ISTEaEAQBVi21dJJME1pZndwkNLL57hgaecorM9/kgeHH+870jQScHnjEFazlC1k5ndWBD9juPJsy+CLy14H/GTgMfFlEnlPV/0pETgC/paqfUNVMRH4W+BNcCcpvq+r5TY/c41kHqwUSc02wakltu8ixZ5Eht9hiBlwl0QZ5njoNnSKAunvR4sbkXE5uleACwB1ZZ7dysUVMQMlyd7PLNAVwgeMiS6kcDDpjrS6wvFrGTTkc5IHhxwv3mZvp+zz83cNms3S+AHyhz/PXgU/0/P0k8ORmzuXx3AmrBxKFsdIp5tMJptrXUXX58C7z3X01OrNj0YDQRIVB3M3GfikLEQHBEJgAVMg0YaH4y73jtAg8C6YbuAWhndcph0MAxEGlK9q2UsaNkcBn4+xSfKWtZ0+xNDjbq3XT8UcvNC5RBqIRKuEwzbxGblMkd8VWLo0yQwvdexBadn7H3tdmEQyhiQgkLAy7FqsVBaTILnJuISkKzgCMGCcJbduEEhGbMm3boJ03oU92j2d34w2+Z8/QLzgbSIggtLIa9WyGrEi5VJTIuN60qZ0nMiViUykE0UIa+RxGXSpjbCq08nvX2HdiDqltYSUiV1cPYIhwyaCJ0/4pCImQwtCjFDN9ZTg+jIghtym3m28uinn4hiz3Bl4P37MnWBqc7fyAK6CaSydo5nOkeYvUtshsm8ymXJx/nqBTHavWqVnaJgKEJsLgDJzeU26c/gimKOBy5KTYYqYvhX8/kJggiAkIiU2VQCIiU2Y4OuJknfOmq+Q10aLrfLekGjxbizf4nj1Bb5PxXiJTIrVtrHUtCMFl3xgCcpsy2b7KfDpBPZtmLh1HFRLbINU2bdsk1YRU29xbfvv+5KTkmrGQ1qlQBKA7FbmqiiluDEPRKIEJEXGB6lZeI9OESjDkKnt7iExpWUMWz+7Du3Q8e4LVgrOpbbvceVWMhJgeV0Rq21yef4GR0nHq2TStvN6TCul+3/umvhe7aJYPFHo7DkFdSqUoUR4zEh/lYOkkIobIlGjnDW63LpLkjW4sJArKRV2Cr6bd7XiD79kTrFblqaqFLo66RiLd7BRHrpkrUFJBBAKNyUlWPN5exmXphIgGGAm5f/g9iwqtkmaj0NZfuGmaLGA4Ogz4atrdjnfpePYEvVWeHVQtc+kkia1jSbFkrlOUTbqZOADGhLTzQhMH44w/Af1VQfY+GSkqSmRimtlCsNpqzmT7WnF9TFGj4OIfM8lNQokYiEa6TWSmWteYTya9X38X4Wf4nj3B0irP3KbU0hlSbREQIQSLUg9TTQgIiExMbi1ZUYQFdFcDQNFxar8ZLMVqTiObZzYZ76Zc1lNXjDUSH2MuGV+UDWWxHCydJMmbfWUsfAbP7sAbfM89xWoqjp0qz06LwnJQJbCG0JTIbZl6PoNiCQidH1sAFVSs6xilCR1dnP6SBvsHS0Yzn2eqfY1mPsfpgUe6cZLQxIyWjhdN2TOMhFh1cQCvh7+78Qbfc8+wkgha7+zRSICRgMDEGAm7Lp7AhAwwQjOfK2QFFNGAOKwwEh6nkU1Ty1qLZvOyzz2eAQGD4SiZplypX+BI5Wx3m4ihFFS7f7fyGrkmXg9/l7O/P9GeXUk/H/BKefb98r87M1FTNPTutPqzmhOZEqVgkEBChuMxBsIRlJyBcLTbo9YhLLQf3J90rl8n5RJlWZwEFtQwA4n6HWbRfp6dxc/wPbuKlWbxB+MT6549djJ24uJ3I5/rhl+t5mSSgkIzny9UHV1eejOvu+IjYow4MbUU11927yFElFyAts9NLaREIAGJbXVn8rlmq6phrmXQfQbPzuMNvmfXsJqU8Y3Gq5hVZpC9xmZRxo6AuGawbqVADlo0D8+VwASUzIBbAfRIH+cqezY10+nll1FR6GOkDSFRUCLHtUUsuiwSmdKqapixVrwe/i7HG3zPrmE1KeN2UeW5Er0GppOx88bcM+SaEgWVrluoZAYIJKCRzyLFLL6V192NAHCuHLB71NiDUDUj5KQ07RwBETlpsa2jp5+T5E1XomWdpn5gQlqZm9UPRCN9ffFeD3/34w2+Z9ewWrVsYGJym6579lgOBzlaeRu5zYq+rhnNfI5QSqR527UMVIv0VJ5K0fTbdg3g3iMgRCUvrnWPpMKirCTt3gQSbZBkDSKp8Ob8dxgID1AKB1ZMs/R6+LsbH7T17BpWq5YVhBPVhxARWnmt+yMiK84eS0GVKChRCYeK7dI5GJbc/SVSVN5SaMhky46zl8hJaRYN2AEy2utKQbXkBBLQto2u621pQVUn2D7TvgnASOkoQ/GYN/a7CD/D9+wa1uqJOlo+wWj5xLpnj73HC2ShoUmSN6AQUOu8VtX2uDb2Nht3V0khqayuSYwqmS5ugrKelFnPzrOpGb6I/KSInBcRKyLnVtnvooi8ICLPicjTmzmnZ+/S8QGvNovvdFM6WD6x5uyx93iWHKuWdl4HEarBMGLEdbBS9vzMfivIrbtGeTGz7wTKN5Iy69lZNjvDfxH4G8C/Xse+H1bViU2ez7PH2WofcO/x6vEMNxqvFWqQLvUwNBAQUcum9lHOfUceufN7YwQSYLHdVdhqwXZfcLW72GxP25cARPanyJTn7rDVPVE7xxuKx6hEQ1yef8E1+pCAKCgzn0wWgcvO53gv5t0v0Glt6Mz9et5rEdoV1+sWEUJZCJSv3jfYF1ztJrYraKvAn4rIMyLyxGo7isgTIvK0iDw9Pj6+TcPz7BeGojHK4RChiYmDqjN+SreV91439oaw2/bRGf2g+HFXYOne7ke716Zkql1XWWfVtVqw3W33BVe7hTVn+CLyNeBYn00/r6pfXOd5fkBVr4vIEeCrIvKyqv55vx1V9XPA5wDOnTu3t799nm2nX664iIAIomaPi6YVchHd9+mKzJwiaG9aag4IIZErXEMYjY8zFB9iuHSIoWhx7GStYLsvuNo9rGnwVfWjmz2Jql4vft8WkS8A7wP6GnyP526zNE4QSMil2neZbt8s3B17xeibnjx7LZqbmCIl1RTG3nSNvsG5bMpmAHDBWJEAVct8NgUitGyNkqkuyrzxBVf3Dnc9LVNEBgCjqvPF4x8BfvFun9fjWY2lcYIj2f3dG0BHIvnex+JcMoIRQ6AlgsA1ezEY8sJV01lGa9HQPDJlFKWRzRJJQCmoktnEdQsrMm+WSh37gqt7g82mZf51EbkKfBD4soj8SfH8CRF5stjtKPCXIvI88G3gy6r6x5s5r8ez1QxEIwxHRygHg4QSs/5uV7s7YSEkBhRUCCRC1PQUmAmG0AnGFXUKmSZFS8isWBEsvD+r2arNyjeSMuvZGTabpfMF4At9nr8OfKJ4/Cbwrs2cx+O52wxEI8RhGZEx5tMJ0ry5zlfu7jBTXtQXCEJCnZAStmhaHknssnSUbhYOarFk3e5f0tu7VhbMhc+8uTfx0gqePc96eqx2/NBREBMQsje+GkJsKpSDIUZKx4ikgioYgm5KphbGXtHieUNmUyjcO4GERZP3oCs3DT7z5l7FSyt49jQbKfnv+KHf4jnajQZxUCbJmoUm/mI23ut2ccqne30nOHx3VgkBIaGJUFVCEzNSOsp461KRiWOx2kYQYlPGqiU0EdaaIsCrhCYitS1CiRmODyNifObNPY43+J49y2r6+iv1WDUScKhymonWFVQVFdvHHi8N6q5UsOXy2wMMcVAls23XQ0sTDMZlwUuEqJDQ2KJ3vTCmjgvGiGvWPpPcIjZlSjJIS+uoOh+9kYDDlRO0szqK5eTAoxgjNNN5JttX3axfE7I88Zk39zje4Hv2LHda8j8UjTEcHWKyfbXwZfcryDKARQiITAmrFqvZEk0el+5YCgYITYyIMBgeZCa5RSgRuWbEpkKmaSHotlkWbjyCcefEMBwfJrMJVnMMIeVokAEZYS4Zx2ruFDSzOcrh0KKVz4H4CEeqZ33mzR7CG3zPnuVOS/6NBByp3MdscotUO26XTka7IRBnxBvZHII4Y4qFJfn7hhCDIdUWWMGYAGMCRkvHEAwzyS3atkG+SmOXjVAxw0W2jYA4F1U1GELE0MxaWM0JJSY2ZUQMo6XjJLZFK6sxWjrJiYGH+q54vA7O3sEbfM+eZTMl/4rLUCkFVRTtuTmom/WrMBwfYja5vUqhlpKRITbDogyaEQThvqF3EZqYG/VXuVw7D+puDispdq63GCyxDcrBABblzMBj1LIp2tatHHJNAen64ln0/iwH4sN+5r4P8Abfs2fZTMm/03+3xEXVaWTKzmWjlsy2GaucZC6ZYPU8fCdWoKibVZuA1LYYb15iLp0gsylGwm6f3UgqpNovHVTWYfTdOCyKEUMtneT+4ffQzOZdNbEJudV4a5lYmg/C7i+8wffsWTZT8h9IhBCQa9YVGwskAjICExNLlVo2vWqmjiUv5AoMpdAFbVt2jqn2DSIpQREwDU1MO8/JNEEIoJBB6KCFts1aOGOujMTHyDShmc0vcseUTNXLH+xzvMH37GnutOS/FFQZCA/QzhuLUjqNBAyEB5hJbxISk7J6gZaiiEIzmyc0cXHTaBGaiFxdb9lqMExH2CwgJl+xI5ULFC9Firz6UGKGo8OEJibLk2UxCi9/4PEG37PnuZPA40A0QikcIDIVFJeBY8TJELggq1AKqjTtXPGKfrn0HQ1KgwgEEjojrIK1roFISotUEyJTIs8TlGyFY0EnK8h09ew7ypbO6FtyGtkMoYmB/jEKH4Td3+yFckKPZ8vpar4b5zsXcb+NMRwsnUREKAdVAqJFejNLURQRgxSzaCMGRDHGICJEUobCFSMYjISu1y5h8eOqXwMiQokJJeZA6SgHSycITNTNGgpMQCUYQoGZ5CaBRN4v71mGn+F7PCuwkgukns4wndwgDquUwwFaWb1ogL54Zt7pLGVMuNBREArVSfe3EcNgNEZOglVLJRiikc6Rk3V994IQmBBVS2BCrM05XL2v8PMLzWweEbqzfauWsdIJ76rxLMMbfI9nFfq5QDrZP5lNGYmPMctt2nnd5dtTtA6RGCMBgQnJNXVZP9YSSszB+BSNfIbENun47kvBAKdHH+VW802aeQ1rczqdpkJx8giuatfl8id5k9CUKAeDDEQjpHmLXHMCCcjJuysKj6cXb/A9ng3Sm/2T2YRqdMA1CLSGgXCE0JQIJSQKyqBQy6apBgdo2zqCc+lUw2FUlYOlkwxEIwvuF4GIEtear5Br6lYDheGPTNll9UhMORygkc8WL3HSDR1aec2Lm3n64g2+Z19jNS9cNi0iU1531spSd09aajHVvkY5GFq8o7jg6dHq/V13UL8MmaUibwfiw8wktwp3TuRm9hJQCQ8QBTFj5VPMpRO+raBnQ3iD79m3bERJsx+97h6r+ZoGeKUMmb4ib8EgATFz2TglM0BgQkKJiQI3vs44fV69ZyN4g+/Zl9yJkuZqbKbIayWRt0o0hBg4WDpJZMrLVgX9gsqVcIhmNk+jNbuhFYtnf7Apgy8ivwz8NSAB3gD+nqrO9Nnv48C/BALgt1T1n2/mvB7PZrlTJc1e+rmD7qSwaXWRN+e7P1g+0Xdr76qhldV4a/65O16xePY+m53hfxX4p6qaicgvAf8U+Ce9O4hLF/g14GPAVeApEfmSql7Y5Lk9njvmTpU0O6zmDtpoYdNmRN46bPWKxbM32VThlar+qap2JP6+BZzqs9v7gNdV9U1VTYDfBz65mfN6PJtlM0Z2qXHt/KgqV+oX+rZQXI1ekbdeNhKA7axYlo57tabjnv3HVlba/nfAV/o8fxK40vP31eK5vojIEyLytIg8PT4+voXD83gW2IyR3Wrj2vH/iwitvNb9EZF1B2A3u2Lx7A/WdOmIyNeAY302/byqfrHY5+eBDPi9fofo89yKTTxV9XPA5wDOnTt3d5p9evY9mwmy3g3jullhs61wC3n2PmsafFX96GrbReSngR8FPqKq/Qz0VeB0z9+ngOsbGaTHcze4UyPba1xVLYltdcXVVPWOjetmhM02o/3v2T9sNkvn47gg7Q+r6kpNOZ8CHhKR+4FrwKeAv7OZ83o8W8WdKmmGJqaZ1Whms12fvSUv5BSiuzHUVdnMisWzf9isD/9XgSHgqyLynIj8BoCInBCRJwGKoO7PAn8CvAT8B1U9v8nzejw7hpGAk9W3U8+myDRxTktxxnUwPMi1+ssbDtxuBZ0Vy6mBd3K08jZODbyTB4Yf9ymZni6bmuGr6oMrPH8d+ETP308CT27mXB7PbiLXjMFoFEPQFS2LgjKCWXce/93A6917VsNX2no8d4AL3Moi0bLF231WjGf34RugeDx3gM+K8dyLeIPv8dwBW1Es5fFsN97gezx3wFYUS3k824334Xs8d8hmi6U8nu3GG3yPZxP4rBjPvYR36Xg8Hs8+wRt8j8fj2Sd4g+/xeDz7BG/wPR6PZ58g/QUudwciMg5cuoOXHgImtng4W4kf3+bY7eOD3T9GP77Ns1vHeJ+qHu63YVcb/DtFRJ5W1XM7PY6V8OPbHLt9fLD7x+jHt3nuhTEuxbt0PB6PZ5/gDb7H4/HsE/aqwf/cTg9gDfz4NsduHx/s/jH68W2ee2GMi9iTPnyPx+PxLGevzvA9Ho/HswRv8D0ej2efcM8bfBH5ZRF5WUS+KyJfEJGRFfb7uIi8IiKvi8hntnmMPyki50XEisiKaVwiclFEXij6Az+9C8e3I9dQRA6KyFdF5LXi9+gK+23r9VvreojjXxXbvysi773bY7qDMX5IRGaLa/aciPw/t3Fsvy0it0XkxRW274brt9YYd+z63RGqek//AD8ChMXjXwJ+qc8+AfAG8DYgBp4HHtnGMb4TeDvwn4Bzq+x3ETi0A9dwzfHt5DUE/t/AZ4rHn+n3b7zd12891wPX1/kruDbnHwD+r23+d13PGD8E/J/b/Zkrzv1DwHuBF1fYvqPXb51j3LHrdyc/9/wMX1X/VFWz4s9vAaf67PY+4HVVfVNVE+D3gU9u4xhfUtVXtut8G2Wd49vJa/hJ4HeKx78D/Pg2nXc11nM9Pgn8rjq+BYyIyPFdNsYdQ1X/HJhaZZedvn7rGeM9xT1v8Jfw3+FmBEs5CVzp+ftq8dxuQ4E/FZFnROSJnR7MEnbyGh5V1RsAxe8jK+y3nddvPddjpz936z3/B0XkeRH5iog8uj1DWxc7ff3Wy269fsu4JxqgiMjXgGN9Nv28qn6x2OfngQz4vX6H6PPcluajrmeM6+AHVPW6iBwBvioiLxczjN0wvrt6DVcb3wYOc9euXx/Wcz3u+uduDdZz/mdx2is1EfkE8EfAQ3d7YOtkp6/fetjN128Z94TBV9WPrrZdRH4a+FHgI1o41pZwFTjd8/cp4PrWjXDtMa7zGNeL37dF5Au4JfmWGKwtGN9dvYarjU9EbonIcVW9USzpb69wjLt2/fqwnutx1z93a7Dm+VV1rufxkyLy6yJySFV3gyjYTl+/Ndnl128Z97xLR0Q+DvwT4MdUtbHCbk8BD4nI/SISA58CvrRdY1wPIjIgIkOdx7hgdN/MgB1iJ6/hl4CfLh7/NLBsRbID12891+NLwE8V2SYfAGY7rqltYs0xisgxEZHi8ftwNmFyG8e4Gjt9/dZkl1+/5ex01HizP8DrOD/fc8XPbxTPnwCe7NnvE8CruKyFn9/mMf513GylDdwC/mTpGHGZFM8XP+e3c4zrGd9OXkNgDPg68Frx++BuuH79rgfwM8DPFI8F+LVi+wuskqG1g2P82eJ6PY9Levj+bRzbvwduAGnx+fv7u/D6rTXGHbt+d/LjpRU8Ho9nn3DPu3Q8Ho/Hsz68wfd4PJ59gjf4Ho/Hs0/wBt/j8Xj2Cd7gezwezz7BG3yPx+PZJ3iD7/F4PPuE/z8RnMy0jqNMpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the points with colors\n",
    "for k, col in zip(unique_labels, colors):\n",
    "    if k == -1:\n",
    "        # Black used for noise.\n",
    "        col = 'k'\n",
    "\n",
    "    class_member_mask = (labels == k)\n",
    "\n",
    "    # Plot the datapoints that are clustered\n",
    "    xy = X[class_member_mask & core_samples_mask]\n",
    "    plt.scatter(xy[:, 0], xy[:, 1],s=50, c=col, marker=u'o', alpha=0.5)\n",
    "\n",
    "    # Plot the outliers\n",
    "    xy = X[class_member_mask & ~core_samples_mask]\n",
    "    plt.scatter(xy[:, 0], xy[:, 1],s=50, c=col, marker=u'o', alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice\n",
    "To better underestand differences between partitional and density-based clusteitng, try to cluster the above dataset into 3 clusters using k-Means.  \n",
    "Notice: do not generate data again, use the same dataset as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n",
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2-D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABfWUlEQVR4nO29eZBl2V3f+fnd7W358mVmZVbW1tXV6lXdotRqikaNJJCMBEIGhGXLFh4PGnCEAmx5wgRjYELEmCFwGIYgPLYFaDocGoONzQSgpQEJIWGppRZqpFZ3dfVW3aruqq7KWrJyffn2u5wzf5z7Xr7MfJn5snKtyvNRpPIt99173u3K3znnt3x/orXGYrFYLLc+zm4PwGKxWCw7gzX4FovFsk+wBt9isVj2CdbgWywWyz7BGnyLxWLZJ3i7PYC1GB0d1SdOnNjtYVgsFstNw7e//e1prfVYr/f2tME/ceIETz311G4Pw2KxWG4aROT11d6zLh2LxWLZJ1iDb7FYLPsEa/AtFotln2ANvsVisewTrMG3WCyWfcKeztK51Zg98yoXPvUEtYuTFI6Pc+IDb2fk5J27PSyLxbJPsCv8HWL2zKs8/9t/TDhXoXBsjHCuwvO//cfMnnl1t4dmsVj2Cdbg7xAXPvUEwdAAmeEi4jhkhosEQwNc+NQTuz00i8WyT7AGf4eoXZwkKBWWvBaUCtQuTu7SiCwWy37DGvwdonB8nLBcW/JaWK5ROD6+SyOyWCz7DWvwd4gTH3g74XyV1lwFrRStuQrhfJUTH3j7bg/NYrHsE6zB3yFGTt7Jm37hgwTDRWoTUwTDRd70Cx+0WToWi2XHsGmZO8jIyTutgbdYLLuGXeFbLBbLPsEafIvFYtknWINvsVgs+4Qt8eGLyCeBHwWua63f1OP9dwKfBc6nL31Ka/1rW3Fty1KsfIPFYlmNrVrh/2fgvesc8zWt9YPpjzX224CVb7BYLGuxJQZfa/1VYHYrzmW5cax8g8ViWYud9OE/IiLPisjnReSB1Q4SkY+IyFMi8tTU1NQODu/mx8o3WCyWtdipPPyngdu11lUReR/wGeDuXgdqrR8FHgU4deqU3qHx3RIUjo8TzlXIDBc7r7XlG/rx7Vv/v8VyayNab41NFZETwJ/3Ctr2OPYCcEprPb3WcadOndJPPfXUloxvN+llSIEtN65tH34wNEBQKhCWa4TzVW77sUe49GffWPF6d6Xvap+11cAWy82FiHxba32q13s7ssIXkUPApNZai8jDGFfSzE5ce7fpNqTtQOq3f+WToKH4hsNLgqsbNa69JpLbfuwRXvr4Z6hfmSZ/ZJQ3fvQnmHv+Qse3D3R+X/jUE53rdfv/VzvGYrHc3GxVWuZ/B94JjIrIBPCvAR9Aa/0J4B8APyciMdAAPqS3amuxTWzGvdH92cprV8kfHV1iSKeefAkERr/7HgCSVkT55Ut89cO/yfH3v62va/WcSD72SRAYefOdHPr+k4TlGpf+7BuE5RojJ9+w5PPLffu1i5MUjo2teYzFYrm52RKDr7X+yXXe/zjw8a241k7Qy5j2uwJf/tmpv32J2pVpFl65hIoVQalAc76Cn8sAUL82y9STL+JmfLTWfV+r14p8auol0DD60D2d1wCqr08Slms9fftt2v5/FUbMv3SRsFzDDXwOPHTXjd1Ei8Wy57DiaT3YjHtj+WfdXIb61VmSeovBO4+QNEPC2QV0aYCrXzlN+eVLiOsQDBfJHhjsea1eu41eK/KkGbJ82xSUCvhDBcL5aud52z9/zz/9kc5xJz7wdr79K59k4dxl/GIeN/AIKzVqE9PMnnnVunUslluAfW/w+zWm/bo3en3WcR2SMDLGWIMT+DSn5gmKObRSaK1pXJ5m6I23rbjWarsNJ5fprNrr12Ypn71I/co0TuAz9+J5mtfLROUaTsZn+MG7OPruh3jp45+hcv4q4giFE4c6+fltFc/ckVGa1+dRrQivVODQW+7Gzfid49r3yc34KA06jGw2j8VyE7GvDf75P3mc07/+X9BRQnasRNwMef63/xgvn1nXBbIay1Mj42oDBahGyMIrExSOjZI/PEI4VyUoDdC4Noe4DpmDQzSvl+H+pddabbehwohwvkpjap75Fy6QNEOSOCGsNrj42Dco3HYQJ+PTnJiicuEak19+hoG7jhAMDSCOEM5WKL9yaYn7SIcRR3/oFOI4nUkknK8SN0Jmnn6F4h2HcXyPq48/C8D4O77rhgPOFotl59m34mmzZ17l9K//V0SEzGiJpBFRfvF1kihGaW64O1V3Z6va1WmaMwuQKEr33sbA7eOoKCGcq5I/Nsbhdz7I7T/xNrKjJdxcpvO57mt1F1M1Jme5+pXTXH38Wa4/+RK3/dgjNK7MEFUbhOUaA8fGyJQGENehenGS5rVZ8odHcDyXVrnGzLe/g1aaoDSAl83QuDKzpBK33YaxHVdIGiFOxieuNZh7/jyTX3+ei4/9DdFCHcdzWXh5wlbzWiw3Eft2hX/hU0+gwpjcWAlEcHMBAI0r03jZgDf9wgeXuHrG3/FdXPjUE7zwf//pmm6MdmerC596gmufPUNufAQdRTiui5PxSVohUa3Zcfs0Z8o0rs8TzlVAhCSMGLzn2FIjnAZTr3/jRdxsgJfxQaSTW1/QGtWKUHFCXG2glYJYEVWbRAt1kmaEeA46UcSVOhwYxMn6ROXaEvfRiQ+8ned/+48pv3wJN+ODgGpGOL5HuFBHhTGkyVW1y9M0rs4Slmv4g/kluyGLZbNUwlmmmxdpJVUy7gCj2eMUg5HdHtZNz741+LWLk2THSiStCDdrjL2T9WlOlTn0zrcs6U61XtbO+T95fEXu+0O/+uGOP3/updeZeeplwrkqCGhHmHvhApWL1yi/8DriObi5DCpJaFyZoXT/7ZRfucTjP/VvyQwXac1VQSncXAa0SeMcvOco5bOXqF68hkoU2dES4WwFpTQ6jCH1sTeuz6GaEYiAQDhXIQkj3FwWVMKFP/0q2bGhTmD2Tb/wQb764d9Ea01maIADb7mLib/8Jo7nohOFmw2IGy2iagNBkGuzNK7N4A3mVwR3beWu5UaohLNM1F7AlYDAKRCpFhO1FzjGA9bob5J969IpHB+ncGzM+L6bIWhNtFBDfHeF62YtUbLzf/I4T/3So4TlKrlDI4TlKn/7Lz/O//iH/yczz5zjwqefYO7Ma3gDeZyMj/guru+TOzzC7DOvohOFn8/hZgOCYh4nGzD1Ny8y8+1XaE7OUT57CZ0oqhcmico13FzA4D1HWXjlMhqNP1hANFTOX0MphUqMsW+jmpF5oLX5UZrmVJnaxWtorRHPJX90tKOqOXLyTo6//20c/oE3c/idD5IbH8EJAkiUySY6MEirXINEocHEDRbqNK/O8rWf+a2OMqdV7rTcKNPNi7gS4DsZRATfyeBKwHTz4m4P7aZn367w2+6LoQdOUJuYojFVxgk8HvyV/3nFKnStrJ1rXzmNP5gnKA0AIK5LUg+5/jcvcNvffSsXPvVVSLQx9o4g2iV7cAjVCHF9D28wz8CJQ5RfvoTjeeA6NK/NoaMYLRDXm6gkQSUJSZxw+J0PcvUrp82uREP+8AGGT76BiT//BkmtiWjQngsqMYZ/SZ6mGNeOShDPQxDGv8/o2E1+/Xm+8MO/SO7QAQbuOoJuRTSny9QmpmjNLaCimMzYEKI0juOgfUAE1WjhFXI4vkvl1Ss8/lO/wcCJcaL5GrkjB2zlrmXDtJIqgbNUBNCTgFZS3aUR3TrsW4Pf7Wt3Mz6H3/WWVV0Oa4mSTT7xHLlDi9vM+pUZkjgmmm5SfvkS/kCOJEoIZxfIHhwmN1ZCKU355UuoOKF5fZ7c2BBuxieqN4nrLXScEFYa6DhGHAflxTieS/3SFOXvTFC7OkNSa5LUWwy84TCl+44z9MAJqhevoxfqaK1wggDVCpd8DzcXoJIElOA4DkkrpDlTZvqbL1O/PgdKg+sSlqsEw0XiagPHcxk4Pk5zdgHQHHj4XmqXp1FKdWIJjucS1VskzRARiOarNK6Xac1VCAYL5NP7EzdbXPvsGevisaxJxh0gUi18yXRei3VIxh3YxVHdGuxbgw8s8dOvRXs3ACsLl6595bQxkKUBky0zX0WjcbMBSTMkboRkRkuQKPKHRggXalRfnzT+74yHboVUXruKV8oTlWsm4CrG/w6gdUJcaRh3jOcwe/oc4XQZxHy+enGS2sQUA3ceITNcREUJuhmitWJ5FZZSCh0r0BqlNXGtycXHvmEmFs/DywboKKYxW6FxdZbMgUFu+9FHyB8aoTE5y8zT55j+5ssUjh+kNbtAXG2aSSROiGsN/ME8/mCBqFwnd3CIsFylfPYi+UMj1K/NMvm15wiKhU3pB1lufUazx5movQDKrOxjHZLokMPZngK7lg2wb334G6G9GwiGi9QmpgiGix1DdezH3kr19UlmTp9j7vnzqFaIDmOyo4O42YDcwSGak7Nkx0o0Zxeonr+GThRKK5J6y0gqoGlenQVHEMcxK+02mk5mDLFi4dwVECGcr9GaraDCGNWKWHjpdQq3HzSuJ50GbpehmxEoBY6DjhOcwEcnCSjQUUxUa9KaqxoXUhjRmilz8TNPMPvCeXLjIxz94VMceMtdvOOT/4rhN92BE3jE9SYAjudRuO0gqhnhlwoMvfE4WmmaU/NopZg9fQ6AAw/dZZuzWNakGIxwrPAAvpMhVDV8J8Oxgg3YbgX7eoW/HutlmcyeeZXJr5zBK+ZpXJmBRAHgBB71K7O05qq4uQyO51K89zYmv3waFaYZM56L4zqoZoRuhKCMn18lyZpj0lFMOFsBAWIxBtd1cByHa48/h+M66ChesbpfPAGI5+LlMiStEJLuycUEddsTjFaaOIqZ/NpzZA+UcDM+hePjjJy8k+/+9Z/hxd99jGt//TSZkcFUhmEKHcYMvOEwAMNvuoP65WlqE1MkrYhD7zhJbnzxj9aKs1lWoxiMWAO/DdgV/ir0k2Xy4u8+xsKrlyFRBMMDSOCaQGYrdceEMSQJEnjoVkRufBh3IAeeC1GCaoToRHVW9KoVLV3dr0V75R8raMWoRgjpynxVYy+A65gxeS7ZAyVwZfF8XcYeQGtNvFCnNVPmtf/211z49NeoX5vtZPO8/RM/z9/501/lwEN3EVeb6ESRO3IAcR2uPv4sUaXOW//DR3nHJ3+R4+9/Wyf9tU2/1csWi2VrsCv8VehHQG36b19CXNe4VqIYETEGUwAxK+nc+AgqUTQn58mODZmgbqIWDauwuoHeKK4D0Ro7BI1JpxQIyxVzeccxbp1exKrzUMUJXjZDtFDj+d/+Y277sUeYe/5CRwJ6/B3fhV/MUz570RR0DRbIHRnt3Ku14iAWi2VnsAY/Zbn7Zvb0uY6GfLeuDCId107SaNGYLqMTZYy9yOIJU6W0sbfez9STLxK3QvLHx1Bx3HH9LB63Raxl7LsRgUQT1xpoR/qadNx8lsxwkcaVGXJHDnD61/8rh3/gzUskoDOlAjpW+KUCg/ce6wSeYWlWVPse3/NPf8QGbC190V15Cw6CoElsFe4GsQaf3pW0lQtX8QpZ/GI+1asPTC490sku8QZycH0ecQQVJ9C2t9qsiFGKi5/9OlopxPdoXps1QdndRmv84SJeNjBB1T4mnaTWJIljVDkCNCqMl0hAN1IJ6OKdR0gaIZNfe47x73/zknP0mxVlsXTTXXmrtcNCNAloSv4hW4W7Qbaq49UngR8FrvfqaSsiAvx74H1AHfhftNZPb8W1t4Je7puh+08w98IFvGywRFfm4CP34wQ+Z3/vMaJKA50WRS0JfgLECeQCtFaEC3V0nOANZHDzGZNmuZukrp3swSFaMwtoUaxn9VUUUfnOZfLHRqm9OIW4Lle/cprSfcfNAa6QtCIEzI4BCGcXePpXf9/m3Vs2RXflbS26lqZqRsyFV/CcAAeXK7WXuTd4ZLeHuufZqhX+f8Z0tPqDVd7/EeDu9Od7gd9Lf+8JelXSlu46SlxvUnv9+hJdGYDJrz/HwndMoxBvMEdcay1+sMs9ohOFlzX6N2G5BlqIq00TtI17uF9cWTlxbAPeQBZEmDtz3qRpdruiVkWTtEJql64jrkPhtmGaMwvMf+YJVDNEAg/HMyJrQanA0G1jzD5zjvzhAxvKu7f6O3ufnRY26668jXVIksSEuo5Go7XGFZ+58AqVcNau8tdhS/wLWuuvArNrHPJ+4A+04UlgSEQOb8W1t4K2LHA3YbnGgQfvXqIr05xZMNo4L7yOihOiSp3mlbnUaGIMpwZcQbI+ju8xcOIQIoK4Dm7gm93CaqvpHTD2CCTNiKhSAzSS9cF3+/qguC5eIUfh6EFUGNO4Zv6Ta6WN776Y4+Bb38jhdz5I+ewl4kaLya+e4dpXz5C0onXz7q3+zt6n7V6JVAutHWaaE7ww92VenvsGlXAtE3DjZNwBYm2qxgWhpWtoNIKDRhOpBqI9q7XTBzvlUD4KXOp6PpG+tifo1rBfrn/ffq/8nQmuPv4sWiWmUjVRxI3WorHUGKPvOaaQKYxJmiHlly8R11s4votW6WeUWms424frgAg6Tl1QWpvJKFl/POK5FG47SOmeY7iBRzCYN5MXIIFHbqyEl88w9+LrRv7h9UmCoQG8wTxJIzSB62Zrzbz7tUTqLHuDtntF6YRqNIUguARUIzMRbIfRH80eJ9EhkWoZNdj0fw6u2VBLOhFYrZ112amgbS+fQc/lrIh8BPgIwPHjx7dzTB3WyyB50y98kCf/14+jwxi/VCCph2al3FKI1uC5xohqcHwPFYedCSCJEyOVECmcYoDnCGEjXHtA20W7sKrznI47hnWMvg5jWrMLtKbn8QfzOEGQ9uiNcHMBpfuOUz77OvUrM4gIhePj+ANZpKvXwNyzr3H0vQ+veo1+Wktal8/u0navlKNJHPFwxEXEIdFRR9Fyq90qxWCEYzzAdPMiWmI8ApP2jCDikJFsJ2PnZmA3tf53yuBPALd1PT8GXOl1oNb6UeBRgFOnTu2Aj8OwXgZJ9fVJolqDqFIHzIpXRNDKCJVplQY+NWaVL0agTDVDBo4fJAlj6penCEoDZEYHTSC3h/zBttKrqEuDbvU3jmihhk4SwkodEkX1whW8gTyH/85byB8awc34HHrnW6hdnKR0721M/e1LAB35iHCmsmbXsLVE6mD9vgSW7actbBbrEBezw1M6wZOgb0XLGzF43ZW3tahMI57HEQ/BIdYhDg6j2Z1ZIG6G3db63ymXzmPAT4nhrUBZa311h669KWbPvMo3/sV/pHF1BsQYctImIx3XSCqW5uazeLkAv5Ale2AQr5BFJ4raxBRxvYmXz5I7NEwwXMTN+rv91TaMjhJwnMWgdKRQrYjZZ19l/pVLVF67SnNylplnzjF7+lUG7zmGmw0IF+oIwpF3P7SmYV7LtQbW5bMXaLtXBAdFQqJjNAl5r9SXomV3DKDb4PXrChrNHkfpCKUVzaRKPSmjdMJY9g1MNy9ydu4Jzi+c3rZ4wmbZba3/rUrL/O/AO4FREZkA/jWY6V9r/Qngc5iUzHOYtMyf3orr7gQXPvUElVev4BSySCtCa0G1Fl0judsPohoh4VyV7PgQQTFP5fxVklZsVv1Ko9KUTgTiusno0fEu+fE3g+uYCc4VxPeMLo8jhHNVFl6ZIH/kAE7gM/rwvUx+7TnC589z6B0njdGfr3Lfz/34ilMud9F0V/Aud6314/KxbC9t98rl2llmW5dNHAiXKrP4bpbD+bUVLbsNHoAvGcKoybnyN8l5A/2t+EXwnADRDm66yl+Ir5N1B/Z8h6zd1vrfEoOvtf7Jdd7XwD/fimvtNLWLkySNFsFADhX4JI0WOlZoFOK65IaK+LcXiKOYhZcuosIYt5Ajmq+aScERI3gWxkZvx3UI52umy9bNhkollxXoJAYEnfMpHB1FJ4riHYfJDBfJUOTQ97+Z2dPnmPrmWY6//209q2p7uWgu/dk3lrhoZs+82snlr7x2laQVUbr7WOccVo9n5ykGIxzlPsKkTqiaaBQavSL7rJfrZrnBayV16vE8GkXJGV/XWE83L5J1Byj6BzqvzTQniJJm5zVfMqDYlnjCZtltrf99X2m7XhCwcHwcN5chaUV42QBxhKjeBAEvHzBw1xFaU2Vqr14xssRxkrp60hV8ojspm3G1SXWhAa7gBB5aL/af3VX61fNpH6O0cW05kNRappLYdQlKi3/I+UMj5H7oFLWJKR761Q/3PN16ekXLJ4S4abJ9AAbvPGL1eHaR6eZFsl6RojMKGMNdjWY5O/8Eo9nj5N0Ss+HECl+14BLrEF8ytJI6c60rxCrCc3wi1SRwc2sa614TRiupoXSC03LJeyUCN7dtq+bNBlx3W+t/Xxv8foKAJz7wdq5++Rnmnn2VKE6MsVfa5NUXclz762fIHR1FK2WqaGtN4lozzclvO7sxmTqi0VEEIYjvoeM9YOzhBsdgvpQIxK2Iw+96E2G5Rv3qDNPfPEtUqePmMhxMWyj2ottF05icZf6li7Tmq0iqV7R8Qhi6x8T965enO1LNVo9n+1lvpd5K6lTCKURctNZEqsXF5hny7jBZf9F1gwKtFYkOCaMm9XieWBu9JVd8FqLrDHIQ38lSCWc4v3B6hWHtXiFXwlkq0fX0HEKkWp1ziDgrVs2VcJbLtbNUoilAGPTHOFK4t2+DvV7AtX2fqtEssQpxxacYHFgyKXRnHLW/2+Hs3bdcls6epB9FzJGTd/LIf/wX/O3/9gkmH38WnSicjI83mDcL98AjnKsgjkP2wCDl2YpZ3Tuy1JAqjW4sionpaIczdG4UgeDAYGclj2syj8D073V8l6TRwhHTF3fhOxO4uQAnG5DUW1z/m+c5/yePc8c/+IEVp25n5agw4vo3XsTNBp22ic//9h8TlmsdAbs2pbuO4mUD3vHJX9yRr7/fWc3Ida/UG/ECjnhoNL6TwXcyKK1pJTXy/mDnXJ4EhLrGscIDnCt/E43CEx/X8fCdDImOqcdlfCeklVTJqpU++dHscc4vPM1cfIWmMt3l2n9o9WQeB7NzLgYHlqyaK+Es5xeeppFUcMVHa818eJWWqnNH8S19Gdxe8Yf2bgRgovYCSmka8QIgRLRwIo9mstRFtZta//va4K8XBGy7e2ZPn6N5fY780VEKx8ZQYUxrrkJrpmI0ZJoh2fEh6hPT6FQLv29d+72Am+5GlseRHcHxPfNTcgnnqgTFHHJgkHihYfq4DOZBa64+fobqhWs4gcnLdgOP/OED6CThpY9/pqfBb0sml89ewsn4oCFpRYy99X7cjE/19UnCcm3VNE3L9rOakWuv1FHGB40WEh2Bo5luXiRREUqWyoe0fdXFYIScN0DJGSdUDSrhFEonCA6hahDrFnlvGN8xbp9GvECompyLv8mR/L3EKqKl6ugl/2CNX1IR00xq3B68eYlRnW5eJFRNBCFSzfR6QjOu9O3rXyvg2r5PTTWHIx6ueCidEKkGgTvMldrLBM3cruTeLxnvjl9xD7E877t+bZbZ0+dIWhFP/Oy/o3FlmuIdh6leuk7j2hxRuU5UbTBwfNz0jw0jnIxPOF8jKBZAaSOhDKv7xbdS/36ztMeizWpduyyRWPaHCiSNEBVG5A4dIBgu0ppZQFUbgMIbHCBaqJE/OkZ2pEjltSuoSMjfNYo/kANAK0X9ynTPy7cL3r764d8ErQmGBhh5y13kD42glcIfKnTup9XQ3x1WM3LtlXrbuMc6QpOQJDG+k0VpTaLrXKu/SsbNEzh5XMfprLrbrpmMm4dgrGPUAyeH5wTkvGLHVeSIZ66pGpyvPGP89eKidMLiH9PiH5Ui4ULlNBdrZ2i7bkLVINERUdJE0mIxrRXNpMZ861pP99Fy1gq4tu9TrI0rB+jUCCiVMB9dZSRzbNeziPa1we9uyhE3W0x+7TkADr3jpGkWvlDDzQbULlzDzWXwilniWtPIBhwYpDVfI2m08AdytMpVwukFI7UQJYtSC8uNe/u5I7u+C3CCAEl3KLgOrueC75G0QsR1EW0MrZsN0EnCoXecJKo2mHvhAq3ZBdDaGPsDZtvuZn1UK6I1Xe4Y/KhSJ39kdNUxjJy8k+Pvf1vPgqsDD97d8eVbDf3dYS0j1zZWlWiGhdZ1TFmPpp4sYLaLDpFqonVCgzI5t8Tl2lmkJoSqQTOpkHOHyHuDxE5ErFt4TkCsQurxApFqdKp5Ex0TODlq0RwiLq54xLrVa8iApqEW8FQGRDOV1BDtGlVYBEfa5UeC4FCP5yn4w+sa47xb4mLjDAqFL1kybgHHEQ5n72a6eZFItfAkINExrnipyyqgFs/hSaanK8ga/B2kW1Lh2mfPEBQLHHjoLnLjI0x/8yx+Mc/0N88inku0UEcphQBKKWoXrpEZLeGX8hArahcm8Yp5Bu86SuXVK2v3ld1u2k1NEg2BB1G8ZDUPgICT883k5Ai6GZFIDA44gU/hyCg4QumNx8kdOoAjxt0yePcxTv7Shzj7e49x/k+/impFxJU6mdESubFhahNTRJU6WimiSp1ooc6bP/ZP1hzuWt2wrIb+7rJeVkk7TbLulE0hlooxjSEcAsmQ6IhEx3hOlkTHlMNFLXvlJJTDq8y3riAiDPgHKHjDNOIKC9EUAIGTI1ItI5CGR4xxH3k6k4qnrV7PEhPi6QCNIqGJ1kZwDaWXfFZrh3LrOqEyVfSeZLnMWe4Lvq9zrko4y2w4Qc4bIlR1ItUiSUKOZ092jPZE7QUCJ08tnk13H2bCrEbTlIKlWpE7mXu/5Lo7fsU9RtugtP357QYlfqlA0miZrBHXQSVJqnrpEjdauIHH8R97hPmXLpI0Q5N/nwvIjQ2hwojG1Vm0gKqttgrZRpRGPAdvZICB28aoXLqOqjZRUfcWOO29q7WRgHB0x2efOzhE5sAg2bEhvvvXfrpn/nxtYhpJlSSSKKb6+jUyI4OMfPfdNC5P07g2S/7IKG/+2D/p6b/vxnbD2rusl1XSdmVk3DyJjoloopT5d6a7VjyhqtFCp2tqj9nkMlqSNINf4emAajRDM6kSODny3iD1uEyoGigd4eCTYLJx2r56ve6KSpMQdR1nrp4QI6m8l0q7FqkkTIO5EOo6s63LS+SW2z76rJ+hQAmASLWoJ+UV90kRd7J0Cn4Jz/FxZKki7U7m3nez7w1+m+X+/NJ9x7n21WdBKxzfR9wcSb2JBB4SJXiFHLlx84/h+jdexAk84nqTpBniF3JkT76B2dOprG9bOllps4L2TSaKSpRZfTtpe8R2YdMmEd9FI0TzVebmKiCC+C4SuKbCN9X9UY0QJxsgvgtK4fgexTuOkB0d5PA7H6Q1V1mSsQRw/k8e51v/6hOEc1Vj7CshWhnVTdUyk973/NbPrmvkl2NX8nuHXmmYdww+uMrRDrOtyyQ6IkwaXUYZlI7STaXqWok7JIQkhLjax8EjIiKkgat9EhWTiAm8aq2JdSudOoycg4OHWnK+tTHHGUPvikvGyRMmLWKaOPi029QpEhzt4jo+iTZFZN0ul34qZFfLvmlnOu1W7n031uCnLHcruBmfwTuPEjdaxAt1cB38Yh6dKBJXTC9YIDc+wsFH7ufa156nduk6tYkpxHNovjKRSiJ76HYRlmOUI/3BfMd/n8QJuhmiEoUKw80bfEk1b7r1+bVeFEhbli6qwhAnMBkyqhVRm7jeSRltZyy1s5WuPX6a2WdfRSsjrxCXG2bycATtOiRhRP7oKJf+7BuU7jlmDfhNyEbEvSrhbBoMDY3kSMcI6/T/jaknNbrmtcWkgIQoXWGb4xUJThoHiFUrdbu0z2d+BE1G8kS6gSsuWqt0Alirn/NiPUxL1Y1bqPP6YqAtJiRWUef1+da1zhk2UyG727n33eyBBqt7g7ZbIRguUpuYIhgu8t3/5me4+8Pv5cDD9yGOceu4uYDs6CCqETL/yqVULdOn+IZDjLz5TpO22Yoh3Ui2/5lLWoiVNFok9RaFY6MUjo1CYgq2tLAyLXIjCEuCxOL7xrgvZ3mgWJnXRAQcMSmn5Sr1a7OE5RoS+J2mJJXz1xDPIWmFxtg77WtqXM8lGCyQ1FpW0OwmZiPiXm3//VBw2Bj9NFArOLj4qbHvChr1UEnvXqlrEmIdUY3nSIiIaaGJF33tmAItjQmKemRwxF3H2LevrVEoRC+aPLMbWS2rQlOOprhQNokc3Zr87eKyRIfrKnSa/P/TXK6ZCvGjhfu5Y/BBm4e/F1jNrfD6Z75G/vCIMWitiKQZMvTACRpXZvCyAYXj4+RaUUdL5uJnv46T8WlNl8F1cFzHNDVHd9oJjjx0N3G1SePaHM35KnoLNPLF9031LqDjqO8JRCXGnaPTzl3ZsRKzp89Ruvc2vHwGLy1OiysNvIEsca2FVrFJ5Wyrh2Iko8NyzQqa7WHWkwbYiLhX+1hJFxa+ZFIVzRhXfJRqpo1KvCUr+bVY23grIh3hiINolyaVvs4paaMUVzzMFnjdj3R2Fxdrz7IQXacYHGAkOEY9Kfe9St9tKeRe2BX+OoycvJPiicP4gwXChTpuNuDgI/cz9t33UnzDYd7xyV/koV/9MDqMOloyfqmAP2DcP2hMwxDfw8kEDN59jOzBYXSUMHj3MR7+d//MCJE50ntF3i+OdP7wUD2KqJYj3b81bsYnGMgxcPs44rokrYg3/cIHTWvC9vcazJlG5Y6zmAmkTVwiKBXQcdzJsrHFUXuPfqSJu9sJtlnNdbH02EX3iOCQqKSzMjd56d0r/fVY/RhNTKJDYhr0ZbnTzyhiIt1E6RhhvZae7a2yQpGwEE0x27zMZPM1RrPHuW/47X2t0ndbCrkXdoXfByMP3rUiT7w1V1li1LqDvqX7jhM++SJOPoOO4tTFA9mxIYLBAsMP3sXbP/Hznc9+9ad+A3y5MbkFJ62S1YITmDaKfTVWcQQvm8HJBagoIRgq4GYCcgeHyR05QOme2xg5eeeS7zX6Pfdy9UvPGDdWxjcre1fj57OoOAZxyB8dtcVRe5S1pAHaxivvljhff4ZYhwiCJ6YI6vDgygBjd8qm7+Q7uwCdeu8BPHxc8Ui0h8K4OlcWqLQNfL8BLIcb9X8mnTGsxdL8n0SbIHI9mefF2ccpZcYBtW7F7G5LIffiljP429ECb6088V7H5A4OUbr/duJGy+TvxwnFE4fwClka1+doffkZvvgTv9IpLHIzPtFCrePuWcFq1bmeixt4qChCi+AEPnGXXs+aJJokjHEyPiqKaS3UcTyThlp59QonPvCOFd9r6L7biapNrv/tS+hWhFvMkR0dJG6GJPUWw/ffzuDdx2zbwT3KegaoEs4y2Xwt9XWbIGxEC19lO8d3u4QElyiJqKpZNAoRl1g3V1w3VM2ugC6s7jdf7flyNhPsWnnu9fL5NQqtExISGnqBsNmk5I3jyNoumt2WQu6FaN3vrLrznDp1Sj/11FN9H9+tftltmLeiBd75P3mclz7+GepXpskfGeWNH/2JFamHvSabs7/3GDNPnyMJIxzPoX59Hh1GqCjBzQbgCHG1QbRQQ5y0N+56mEQG0OAN5nAzgQmwui4qjFjRu3Y1hLSDlTL6N4UM/kAeJ/A49sPfw/f93s+v+r0A21v2JuP8wmljgJxFA9R+fsfgg5xfOM1s8zIanfq7IUyaxDpMc8k9wrhB1isQODnqyTxGumCcREfMt64QEyFIlwHdS1oiN4pJCDWTlsYnAwIHssdxxO3cv+V0+/C70zGPFbbXhy8i39Zan+r13i21wu9H/fJGmD3zKpf+7BuMvPlODn3/ScJyrWfq4fKg7+yZV7n+5EsmM6cV0pyroKMEN+MT11uI56ITlRr/DFqbvAZUwmp/LxJ4i64fR3AzAf5QASfj07g8Y4x990ahl4RDmpff1u0X3zVB21iRPzwCjsP1J19aYegf+Jd/f8X3tdw8rFc120qqaVWsaTofq4hYN0m0IkkiHFy0aGIV0kwqeE4WT3yayUJ6BUn/171i3ovGvtcktNyt1FsXRXARcVEkNOIFSsH4kh3S8oB4W29ot9Mx22xVi8P3Av8ecIH/pLX+jWXvvxP4LHA+felTWutf24prd7MVLfB6rWZvZCJp7zZ0ktCYLuN4Lkm1hfguUbWBG/h42QAVJyTNkPzRUVQYkRkqUp+cozk1D4BfyKISRdJsISI4nos7kCNphqaBitYExQK1cp3ufOMl+A6O66FaptBKR4nx87ddSBrEddCxojlVJjc+TNxo2Ybhtxjr5YNn3AFcMT1itdY0kkonu8YjMJXmWlCSoNEkOsKXTCdw2y6IWr8CdrfpNT7zmqnAddJXks577YpcjSLUdcChldQ7LppVM3IKD6xRtLbzbNrgi4gL/A7wHmAC+JaIPKa1fnHZoV/TWv/oZq+3FsurZWFjcrqrNUTppcu+3kTSniScwO8YanHENEBRmMrEchXxPbxsgJsNUFGMX8qTXLiGOGL04QdyaKVMa0WtUUrheS7+QI7M2BDZA4MkzZDcwSGiSh3Vnd7ppL4fpfGKGVTgoaIYJP0H7wgkZtLQiUJ8h6jexK3U8Qdy27Jbsuwua2mxj2aPU4lmqIYzJDpJDZ7xHSZEoCXVkkpM2qVuZ+IIsWqlQVnYm6v6/jCTlTHugotHQERjyRHt3y3VoBlXODx496oB8b0ii9xmK9IyHwbOaa1f01qHwB8B79+C826YEx94O+F8ldZcBa0UrbkK4Xy143Nej+6VvDgOmeEiwdAA4XyVsFxbcux6E0nt4qRJVUwUA7cfxPHc1L2CaQYOqCghWqiTPThE/ugocb1F43oZfzCPX8wjjtNJ7QxKBcQzCpZuxiczNoTjOgy98bhRu3QcBm47iGR8I5XgOYgjeIN5cuMjZA8OkTRbqFaExlnM7vHMpKBaEXE9RDUjsqNDFI6NEZQK1K/NcvUrp7n42a8z++w5Zk5/50b/81h2iXbxz9m5Jzi/cJrJ2vklz9tpmcVghDuKbyFwc50yKQd3ScWrSlMVlVYkOqIWz9OIKyTK+LdvJQQH1/FwCTqvObi03T2ahEo0w4XKaaOho5bG35RKmAuvrJkGu9NshcE/Clzqej6RvracR0TkWRH5vIis3vduE/Sqlt2IC6JtpLsJSoWOLvtGJpLC8fFOEZLjeQycOIQ/mEcyPo7v0W5/6OUy4DjUr0xz5Acf4tD3nyR3cAjViogqDZpzFcR18EsFxr73jQw9cAKvkCN7YJCDj9xPbnwk7RQVkDs0gpcL8HIZHN/Dy2fIHiilRV9C4fZDZEYGcV0HN5fBCTzcbAbHdXEH8/gDWQZOjBMMFcgeGqF87jJTT75I0gjxBvOE5RrVC0ZqwXJzsDz3vh4t8FrlaWpRuacRKgYjDPjDHCncx8Hcnbji4uB1PNwuPiDEmPM5aVGTiCIgb9Qob3oW4xBKG9FEQ6e0PH2uiXVMI15AKU05ukaYLO4GumWR+83DXz45b/XksBU+/F65hMun+qeB27XWVRF5H/AZoKdykIh8BPgIwPHja5ct92IzIlyruYRuRJe9nc6YPzrK3PPnSVohJIrckQOoWhM3lwE0Evg4nsvA7eN4A1mmnnzRNAVvyxpHMY3r83jFHPf980c49p5TSzKRWnMVsmNDIOAX88y99DpohZ8NyI4O0ZyeN5NWIUNzqozWGncgi+t7+MU89SszOLkMI9/1BobeeJzc+AittO3g/IsXEBHcrE/SikDD8AMn+nbrbEeKrGVjLHc1hKqOK56RG5ZSz1z87uYknpMxksfpAsW0L0wQcRjJHE6Nl9NZDSfKJdlU2uReYLFWQKNS7f1oUSaFbjkHRStJi8BEmA8nGcveTqxDWkkdz8kw3byY1jMMEji5VfPwd6IydysM/gRwW9fzY8CV7gO01gtdjz8nIr8rIqNa6xWtkLTWjwKPgknL3ILx9c1W6rJ3S/7GtSbhfJV4LCRTGuDAD93VUdpszVUI0gnm8l9+EzcT0JicxfFcHNdBvAyZ0RKZ4QFe/S9f4th7Tq2QEv7uf2PqAS586gkaV2eYe/F1s7ofLRLXG2itCRfqeMU8rbmKabKuNUlkilAyY6WOsW9/99rEFMUTh81uZqFOUCpw4C13kR0b6isI3k+DeMv2szz3vt2RqbuaNlEJ5fD1jp8575aYTSbSgqosQogvAYP+QQI3x/XGebSG6ebrxNpIDSckRLp1yzh12rEJgERF0KUM1F2pq1PXDjiI1jSSBaaaFyh4wx0hOJcApRMq4RQ5b4iCX+p5zX4K4zbLVhj8bwF3i8gdwGXgQ8A/7j5ARA4Bk1prLSIPY/ZGM1tw7S1lq3XZe6VpPv/bf4wT+GilVhRwvfaHXyIYKZpCpkaIAF4xj2pFBIMFGlNlLnzqCR761Q/3HFP7te6VNQjNmQWT7RMnRhO/LZTWjExQ2RHmX7rYMfjt+ESvHc/yCuPV2K4U2f3Eero3/bC8+MeTYEkufiupsxBN4nWtKmeTiY5ujCcBCSF5dxjfyXbEw0LVwMHFxSNmF3o+7BCRbuFLlkAXiGikE4FeYfQFhYiHrzO44tOIy2TdQSLqRgVUXLSGRjLP7cXv6nmtnajM3bTB11rHIvJR4AuYtMxPaq1fEJGfTd//BPAPgJ8TkRhoAB/Se7Tiazt12debUI68+yFmnj6HiKCTBL9UQBwHxzfaNtmxUl+r6+7vMHvmVf7H3/9VyGc6qZ5O4JlYhYbMaInm1Lwp+uoxCa1XYbwaW5Eiu5/Zqu398tz7wMnTTGrknSG01lSjWUAz4I8YP3O6qqwn5U464fKJp+ANE4YNRNIssFX+kterYN37pLLJuoWDR8YZME4eHeI5WZpJpZOC2o5zBG4OjSJSIZGEaTPzCMEhcHL4bmbV/347UZm7JREWrfXntNb3aK3v1Fr/m/S1T6TGHq31x7XWD2it36y1fqvW+m+24rq3Gvf93I9Tuu82Dv+dtyBZn+Zshca1ORrX5lg4dxmvkN2wKNnIyTs59IMPgQhxtYkTeHiFHOI4uBkfN+OTtELCuQUu/tk3UGHUcbtsJgjeDlp3Y0XV+merhLeKwQjHCg/gOxlCVSPvD/KG4kMU/BKhqqFRlPxDBG6u8xmlEqabFzuBQ4A7Bh/siIb5rk/WLaK1ceOsxs1t7BfRaWaSce2A4KXKoG6n1MwVz/S4FZPFo9N2j75kCZwcrrgEbp4Bf/XJ+kYlmDfCLVVpu9dZz6/dNrBnf+8xs3qKE5PC6bmoKGb+ufPc8zMbFyW7/5/9eGelHlUbNCfnSFoh/uAQldeu4noet/349+FlM4Tz1c5Y16qyXY9+9Icsq7OV2/u1cu/bcgttwqRBObq25s6iXaDVW4jsxoXN9ipOaiZFTHg68AbIeQMU/CFaSY1ItVBaGWloHQMOrgRA2Lc7B3amUcqtkEN107Bann93s5CRk3eSHR+hdNdRhu6/ncLhAwQDOfxinmB4gLnnL2z4uu2JZPjBu1CtiOz4CAO3H6I1vYC4DuPvfDOFw6MkrYjyy5f48j/8NR7/qX9L+ZVLSyamjaRjbjZFdr+zEZnizbB8VWlcPLLo4umxs8i7JepJuQ+Z4eVsQv57F9Ht0KxO0KK4q/Qw9w2/nXuHH+GOwYcY8EdICNFoBvyx1Mg7lILxVCk0wpOArFvsGO/V0i+LwciS3dRWF2nZFf4O0q9fu3ZxkqQZEgwXyYwMmhe1plWu3bAPfOTknbz9Ez+/ZOWuE8Xow/dSODxK/dosU0++iJvxac1XyR8eofzi6wSDBfKHzD+6jQZcbZ/aG2c93ZutYvmqUhEz6I+TcfOdY5bvLOpJ2RQjiSbR7c5RbVSai59W5XbSFyV9vd3cfE+G8FZBp+IKCp0onp/9H2ScAsPZw4xmj3Pv8CNL4hyBkyVw8uT9wc4ZugPlu9kYxRr8HaSd9aLCiPmXLhKWa7iBz4GH7lpx3PyLr6OaEW7OVPklrQgvE6zrA18v973bCD/9q79POFcBoHz2Im4mAAERMd29mhHlsxfJHxqxAdcdZif7oHa7fJa7eGDlzqKVVPGdTDpBrFR3dXARabc61ES6hZO2PRQRRBxC1bhpfPztcUq6Q4lVRKzmkKZLM6l0DHX36r1d7NZrst6J9MvVsC6dHeTEB95O5fxVrj7+LHGjhRt4hJUatYnpJe6SEx94O9mxIcJKjaTRIm60CBdqZMeH1qzubccIwrlKX66YbimKcL5qtq3NkMIx495xsj5RGni1AdedZ7u3973oL3DoEKlm6sNfjqTSCzGJjlDa6PEoNEoUrgQUvCECJ8fN5eIxzWBc8Y1UNA5NtdAzkL48UO47mSWSyK2kiifBks/sVGMUu8Lvg62qGB05eSe5I6M0r8+jWhFeqcCht9yNm/GXuEtGTt7Jd/+bn+HF332M6b99CYDDP/Bm7vu5H1/zuhvNfe9OE0VMtsHBR+4H4Po3XiRphZ1qXhtw3R/0s7MQBNdxiZO2quTSlXra4RgHB0cctFYIQsEd6rhBzs4/gVZ05bbvHUw6qUYQAidnYinaTE7dMgtt33wvQ71WoHw3G6NYg78OW10xqsOIoz90yvSFbb+m1Ap3SdvnvhFuJPe97eJpZ9U4gU9QKjD0wAnmXriAPzRAMFzcVAGa5eZiLWMFRjuz5B9iTl9J3T9Gak0wWjsxIS4uCOS9EjlvcEWjkKxTpMGVPWfsYXECMw3YXZROcPBwxDXZOGKiE+2K5Y0a6p2Kz/TCGvx12OqK0c1KOG/XuZcXhQ3efYyTv/Qha+QtK2ivUAe8A8yFE3Tar+EiAq72EHHIeyWGMocA0FrTSqpUwlku185SDid7+v+X6tTsBIJPliRtdN4dTG4Hn118FKZHQKSbRjYHoUDxhgz1TsZnlmMN/jpsdcXoduanb/bcNqvG0g+j2eOcrzxDo5OaaQy+UZfUxnjqZEkrRZNi6jBRe4FaVKZ3M3N22NiTti5MUMSdyuB2cFajaakqWXeQJIlRhKk8ssnN107MSHDshgz1eruo7cIa/HXY6hX5Vuv17NS5LZY2xWCEwMkRSgNHIlw8wCHSDRQJLj4a3ZEHdhyXZlwhVA0SHRMm5jjT+3Z36Mgfo7rMe7urlXFNOThpXFmjJKQUHF4ifBapFvWkvPOD3wTW4K/DdqzIt1uvxxp4y/ajGMkcJVQNKuEUkWqljUIUGTdPxi0SqjqVeJqiPwoiJDoB7axSobuTODh4JIQYITTSEG1b+0d3BWcl1Rr0CVWdAosGf6cya7YSa/DXwa6aLfudXqqd3Zr5BGPMNC+BmEBmMRgj4+bRephQ1ci4eRxxacV1WrrWsz34IqbPlqQuoq1CcPElS6hrmHKwxSpm1cO9JFpIiHHwUFqhtaKV1Jecc6cya7YSa/D7wK6aLfuV1apCR4JjHc38wMmRcQuESR3fyVCJpmnEAb6To+CXFnWBBNDtStveODh4EhAuk5XYDJKeU+lolSPMRCXa7QileU4WlxjPyaTZOQ0i1aQeLZDzijuaWbOV2MIri8WyKqupdtaT8pLiosDJgwhKKxw8ItWiEk2Td0tdukAaryv3vBeKhFA3Ur+60xEu2wyCQ6xbaa2Jm1YA+0ved3BxRNAoMk4BTYIjnqkV1goH8AiYD68w07qE1mpJMdXNgjX4FotlVVarCq1GS908GTdHyR83LRCJ8Z0Mg/4Y9aTcqd6VtFipH9q+dNWzmrc/2v112zuK9qTj4uM5AR6ZTpFVpJsorfCdnBFDUyZrJyFK64Q1vpshcPOUgvFNjWs3sS4di8WyKr2qQhtxhWZSIeMWOm6eufAKQ/7hTt49LObeF4MRRqJjVMN5onW7Y6308JvcfLXidQdvTcO7WECl0nRKhYupEVA6WVIH4OBR8IfIeYNk3DytxEiNHMjexnzrGjimEr1dQLZT2jdbjTX4FsstzGbbJPaqCq3Hc+S8oSXiX55kqMVzZLxFlc12ULMSzjIbTlDKHES3NM2kimK1lf5iALW9+l4s7FqKMfZrh4DTdiU4kpo6ESLV7Mq3b3/eCLq154+CN8x8dJVItYhUK/XjxxT8YeDmzNCBLXLpiMh7ReRlETknIr/c430Rkf+Qvn9GRB7aiutaLJbV6VZt7A64trXX+6GXEFjGHSDvDS45ruCNEKlmT9G16eZFEqWoRXNoElxxELw19PS72wau1ybR+OUd/PS30ymcAvAlwHeyuI5L0R/lQOYorrhpEqZDximYgC2QqLjTbDzWEcPBkU7swtyLsY5s9M2YoQNbsMIXERf4HeA9wATwLRF5TGv9YtdhPwLcnf58L/B76W+LxbJNbJUM7/Kq0LaEcrebx3VchjKHO7LJ3XIB5xeeoZlUccXDlyw4Qpz2g3Vw07TI5UbdpEW6eF16+zp9pz0JmB68SicEboah4DAazUzzEkonpoBKAgI3R84dBFHcO/wIei4hcAqICPOtazhKCFXTFIOJ09Wd6h0Ug5HOLsd0rtI3bYYObI1L52HgnNb6NQAR+SPg/UC3wX8/8Adp4/InRWRIRA5rra9uwfUtFksPtrJNYjeriX+NBMd6Vp62A7Wm3ysEkiVSTZROyLoDNJOFtN612zXTrWkj6cSQ4OLhOr5J3VQNEyTWCTmvZIy1WqzgzblFBoNRAjeXNiAxq/N2XEIrRT0uo7U2fSCQjgJmd7Px7dS+2azLbaNshcE/Clzqej7BytV7r2OOAisMvoh8BPgIwPHjW9e812LZb2yXDG8vA1hyDzIbTvTs4uQ5gdHP13HqCzeGW2Gem2ya7uBr2yWjO64Xt5MiiVGvdFzyXomCP4zvZBjNHudK7WXmo6u4EnQCveVwkoI3guNIZ0Xe0QKKuycnjSsBg/4YIs4SHaD2d+7XEPdrxHej89VW+PB7dTFYHkXp5xjzotaPaq1Paa1PjY2N9TrEYrH0QX/NTG6M5c1Z6kl5Rb5+ohTnyt+kFs0TJk2aSYVGvIDWmsHgIDnXyBToLleNg4eLh0eAg0fWG+Bg7g34bgbfyXZSJBMd4zu5zvcpBiMEbo6RzDEO5k8wnDXuJQ2Eqr4kZ76tBeSmDU1EwHdyuOJTjWY3dY82EjdZrcZheUOVrWQrVvgTwG1dz48BV27gGIvFsoXspAzvcvdRK6lTj+eJVWh83yjQdPL0Ex1yfOBNzIYThEmTSDdSs28Kt4zMsk+sW4g4FL0xavEcrmPiABk3T8EvLVk9d4/B6Pnk0VoTqlqP72y0gESEVlKnES8Q6xBFvKmCqo3ETbbL5bYWW2HwvwXcLSJ3AJeBDwH/eNkxjwEfTf373wuUrf/eYtl+dkqGd7n7qBEvAHRWrp7jEyYNYt0k6xTJuAXGC3eQ90s04xpx1DIqOuLgio+Di+cGDPgjJhCsq4xkj67p496IC6v72Pbk0G40vpn7tREjvhudrzZt8LXWsYh8FPgC4AKf1Fq/ICI/m77/CeBzwPuAc0Ad+OnNXtdisewdlgdyQ9VMDbjx0TvikvU8Eh0xnDlCqEyv5GIwwl2lhzm/8DSNpIIrPlprFBG+m+VI4d6+DfBGOkltV9epjRjx3eh8JSZxZm9y6tQp/dRTT+32MCwWSx90BysbcZXAyROqugmyikuiY1zxOoHWdrvD9mcv185SiaYAYdAf25Cx7zWG9bJetiNDpjsQ223EV3MTbccYROTbWutTPd+zBt9isWw1bcOXKEU9nk9f1Z2MmZtReKxfdjrVcjlrGXwrrWCxWLac7oCx0rHJb3cC8v7gjhjA3TS6u9W+sB+swbdYLNvCbhm+3chvv1mw8sgWi+WWYjfy228WrMG3WCy3FKtp+N+M6pZbjTX4FovllmKxw9YiN6u65VZjDb7FYrml2E5JiZsda/AtFsstRS8N/1s5DXQj2Cwdi8Vyy7GXUyN3E7vCt1gsln2CNfgWi8WyT7AG32KxWPYJ1uBbLBbLPsEafIvFYtknWINvsVgs+wRr8C0Wi2WfYA2+xWKx7BM2VXglIiPA/wecAC4A/1BrPdfjuAtABUiAeDVxfovFYrFsH5td4f8y8Nda67uBv06fr8a7tNYPWmNvsVgsu8NmDf77gd9PH/8+8BObPJ/FYrFYtonNGvxxrfVVgPT3wVWO08Bfici3ReQja51QRD4iIk+JyFNTU1ObHJ7FYrFY2qzrwxeRLwGHerz1sQ1c521a6ysichD4ooic1Vp/tdeBWutHgUfBNDHfwDUsFovFsgbrGnyt9btXe09EJkXksNb6qogcBq6vco4r6e/rIvJp4GGgp8G3WCwWy/awWZfOY8CH08cfBj67/AARKYhIsf0Y+CHg+U1e12KxWCwbZLMG/zeA94jId4D3pM8RkSMi8rn0mHHgCRF5Fvgm8Bda67/c5HUtFovFskE2lYevtZ4BfrDH61eA96WPXwPevJnrWCwWi2Xz2Epbi8Vi2SdYg2+xWCz7BGvwLRaLZZ9gDb7FYrHsEzYVtLVY9jOqchmmz0BrDjLDMHoSp3h0t4dlsayKXeFbLDeAqlyGiS9DVIdgyPye+LJ53WLZo9gVvsWyATqr+unnQDwYOAoi4OfNAdNnwK7yLXsUu8K3WPpkyapea0BD5XVoLZgDvJxx71gsexRr8C2Wfpk+A27OrOb9PCDgeNBIJaTihvHlWyx7FGvwLZZ+ac2ZVTxAfhx0BEobQx/VIWnA6MndHaPFsgbWh2+x9IGqXIbGDJQvQFA0Bn/wDqhOgMKs+A+/dduydGxGkGUrsAbfYlmHju/eHzQr+VYZatfAyYCfhePvwRl/y9qf7zbW+UNQv9a38e5c380tzQg69i5r9C0bwrp0LJb1aPvuCwchNwZJE7QCFOTGYfaFVdMxV6Rv1q7Da38G9an+0zkvP2E+V37V/KjYjGf6zHZ8W8stjF3hWyzr0ZozxhkgrkFmxARrk6aZBKL66umY3YFegGgB3CyEZciPrZnOqSqX4crXYerpdDdRBBWZzKCB46Dmeh+/cMG8ULwdjr7d7gIsHazBt1hWoeOKqUyAXDM593HDGGwdLwZw10rH7J4swHweoD5pHns5yB1c8fnOzqB+HZzA7CiiMsiQmWxqV+DAG5cef/4voDEFThYEsxsIy6g7/q41+hbAunQslp4sccUUbzcZOOVXQQvETeNWyR00B6+VjpkZXjTyAEmcGvsWhBUTC7j+NFSuLHXrtHcGKoagZAy4BqJqmhlUW5oRNH3GvOfmwAvADcDLQli1rh9LB7vCt1i6WFJJ63gwcAwygzB0t8nIiRvm9dwhk63TTscsvQF1/vMwfw6aM6YwKzcKw/dDMmNOnkQm4IsGcc1nwZxPRUsDsZ2dgQPhglnhq9ikgmYGYfDepav21pw5f3vXAaYSWDWhNWezfCzAJlf4IvJBEXlBRJSInFrjuPeKyMsick5Efnkz17RYtoslq3owK+mF82YlHhRh+L7FdMzGJMw8b4zwyAMw+wIsvG7cP3FoXm+W4erXIXfYPJ993riCxAcVYvI5NYgDjrs0EJsZhsZ0GiCOjfF2AjNR+ANw5G2dMavzn0+vW1scOyx+DndHdX/aY1Jn/5v5bfWF9gybdek8D3wA+OpqB4iIC/wO8CPA/cBPisj9m7yuxbL1TJ+BJDH+8XDB+MyVNi6YsGLen3/VGNeB41C606y4579jjHV9EhwfvIwx4MTG3z/9jDnOK6Qr8Bhj6D1jwJOmOX/5VZh+zhjI0ZMmddMJjPEXQCcmcJsZwikeXTpBZQ6YiaZxHarXjCunVYG4DuVzJstHxYu6P9uU5WNF5fY2mzL4WuuXtNYvr3PYw8A5rfVrWusQ+CPg/Zu5rsWylXRWyVefhPJ3IKyBXzIGMlyA+rQx6uGCycUXoHpxMT1y4YIx5Ekzzd6JUx991Rjc+vU0h34A44gX86OTNL0T49IJK9Cah7N/uPiZuA6tWbMLGDgGfgHmXzHjvfzEop+/NQOZITPBqNCcRyvjVorqZrcw8xxUr5rrbZfuT3dW0jZPLpaNsxM+/KPApa7nE8D3rnawiHwE+AjA8ePHt3dkln3PkqImBJQyKZOOZ/zwKjYGNFMyz4OCMWRxCHMvgZeH5pwxaCo19KTCajjG569imHoO4qpxs/QiaZnfbs7k6L/0X4yR9wcgf9gY/uol8AaMDz+qw/wrULobmlNmvJ5vDGx74olD813ETXcIChZeMxOPuNuj+7M8KwmsqNweYl2DLyJfAg71eOtjWuvP9nEN6fGaXu1grfWjwKMAp06dWvU4i2VL6F6Rur7JwNGxWX27WXOMjk1KZuO6CYzGTbPqRkPUSv31c5g/p6jr5Dr11WuI5tcZiAaVQLJgJhS0MdDhAuACiZlwWrPGeFYumr+i5rdB0l2DlzVuo6BgJoS4amoGgkEI5815dWKCz/mDcPitqwZzbzjImxk2127XF4AVldtDrGvwtdbv3uQ1JoDbup4fA65s8pwWy4ZY1YB1r0iDdOWslDG2jmsMV1SH2mXjUpl7JXWXYN7XURpwzZi0yCWIWUnriL5oH6fbOZiOSa+Mq2bZ1J48xE9dQUnX8Z7ZBSQhFMZTF5EyK30RYMhMFioyufoDx2DmLEw9Zb6vX4A4guaXUe0g9I1IOYyeNDsmMCv7uGGymA6/tb97YNlWdiIP/1vA3SJyh4gEwIeAx3bguhYLsE4gsTtPPncwNaRuGkwNjRHPHYSoZl5r+95JFicGlZgKWhTmT8oxAdlgYHUXzpq0N7bp+RzHuI6Ajv+/bew7x0dmLDoxuf3+gPlMa8GIvjWmjdvIK6QuoiZc+pIx8n7BjLMxaYLWV55Y3PVEVRPEXrgI5z61bvDVKR6FY+8ynw3nzW+r+bNn2JQPX0T+HvAfgTHgL0TktNb6h0XkCPCftNbv01rHIvJR4AuYvekntdYvbHrkFku/LJc36JYz6F6RQmrAW4CbZtwMGGOeP2qMabRgXtep0dUak3XTRi2eJ0nTLm8YDUnNyCpIvHguHfY4tisQ3NbnF9cEc3HMKl9SOYigaHz7KoGoYr6v45rr1K9DOGsmRcc1E4OXNxNFWO1rpe8Uj9quX3uUTRl8rfWngU/3eP0K8L6u558DPreZa1ksN8wagUSneBR17F1Gg2b+FeMCAWPU0WmwU0NmwBjKYDh16TTXv25S25rxq1a6q1iL9sSSuoMa08bIeznzPAlTw501q/ZWuiPRpG6nNF6QRBDkzX0Iy2kmUirVEAwsZtxYg35TYqUVLLc+y+UNYEkg0SkeNcZ85H6TjZMdNQFclZjJQmuTejn9AjRnTfVqX2iQzObH72TSCaiP67V3GOh0R1I12jp+3mT3eHnz3XWcZu9I6sZyjLEXoHjCfFbFgGt2NW0pCZtxc1NjpRUstxS9grMdt01UM6vzqGpWsMffs/jB9i7Ay5nAphSMAUSlOfYRNCdT49gOqq6FYNZTap3j+vlSYVqk5fcZAG6PLb1+OJ+mj86aSczNmO8hAl4xzTyqArGZHIKBNIibBoG1Y/SE2umgmWEr1XCTYlf4lluG1YKzgJE/qFxMhcvqxvBNfmsxCNneBeQOmtVstGBW9o63WCDVblzel18+1ctZ1xXTD7rr2v2u0dqrfW1W80kjVdF0zA/K6AF5efO+n/rpxTNyEiJmx+MX0tTOLt2g/CFbTXuTYg2+5dZhrSrPuVfSQqaiMWpRxbQrPP8F89n8ISNBUH4tVcRMtW6CQciMGn/2RgOwOmZLVvgAxOnuYiNZP13jVZHx3w/dBdlhs2qvXzZCb46fqmzmzDFaTKaPuKbhy8DRJRk31K/ZatqbFOvSsdw6rFXlWZkAJM1pd4xGjYqh/B3U5DMm7zx/GGqT0JxedKO0yqbyVYX058rppl9j3z7veuff5ORRv2Z+oGvnEZl0zbgG3iB4aT5/XIPhu3v26VWXH7fVtDcp1uBbbh3WqvKsTKSphk6aTw+mMMoxeeeFI+AkJiArbuovDyFKWHTj9GOUbwS97PdqbPTa6fdrG/eeNQHtqt60GlgXTFaOF6zulxcP5s6aCdPLGQXR7ZJqsGwp1qVjuXUYPWl8zFE9NWBdPmdxUvngtHmJSounMsPGXeOlapdJBG67OtVN0zTTqlcnm77e/Wezk39CG51o2r7/9c7Z9RPXjLsrGOzpolGVyybwHacTYxwaYbnmzNKGLJY9iV3hW24ZOjn17ewRUoN0/s9TPXk/TVWsG1+1nzc+6mjB7ATiBsZtkgZqHc9MBGG6Mu4EbjHndjNmQtkzOCxKV7VdNht1A6WTRFjp7aKZPgPZA6YLV+N62vIxB0GpsxuwGTx7F2vwLbcU7SrPTsZOaz6VJZA0HTFKF+y+CcY2Jo0Lo/xqVz56i45scVw3k4VOUzSNqA2QQFLvPYhdQ5vK4Hhh8+epXwe0kWHuNtjtOImISdOEdIKYB5apj25Uh8ey7ViDb7mp6Hv12M7Y0e1KUTF68VoZQbLmrGlwkjsE+TETmK1eMnnnKkp9+PGi39vJdhn9PUY7AI2GuJK+uIlYg47M5weOrzTY66lhriVjYQ3+rmN9+Jabhg11U2rNpWmGuS6j7QHKpBnmDkDpLigcNJNB4SCMvDGtRs0CSZc0gZMa/537rhtCKTo590sCwD2UySXVCGrHJFZb8zlB75TL1eIkbf99+753YzN49gx2hW/Zc6y6it/I6rG9Es0dNEFFlfaZdTwTYGxXz9bOmRWqlwO3YAK42YOptELcVW2asLdW990r+NVy89Ng8/JxuxkYfwha06YD1pIaAw2StlVsXE8nwEWDvSJOkhlemrpp9fD3NNbgW/YUa/mAN9RNqS2nkCRp0akyP25g/NNJZFwXQaovoyKovGwMe2sm9fdLWnF7IxLH202/KaLLjL04pqK2NQ3j3wPOadOMXbf19zPmnvj5Rf2hZQZ7TTVMq4e/p7EuHcveYq1q2XVE0Lrp6LJHCyb3vjAOI/eaNoBOgFGKbCtEhiabp6NZ45i+tuHs1hj7Tt7/VtKvj76tn58+9odMk/XGNFz+G7OKd7NGPz83ZibEsAy1q2njl+tLXTbrYPXw9zZ2hW/ZW6y1ij/6AxtaPTrFo6jcASjdaSaP+XMmGOt6EM5B5oDJOQ/nF+UFRAAf9PLuVZtgS/R0etGv0U+PcTJABM4g0ILKa5AbTxulzJp7YQ409yGJYPZ58Ieh+SmUlzPxj3XSLK0e/t7FrvAte4s1VvE3tHrsPl/cSMXQYrOSjyt0XB6OZ1wZKjbZO5umVyvnrUSb77Dudbp08lXLFExBKpGsjZFvziyblNRizEIlxv3TKptdQX3KCqXdxNgVvmVvsY4PeMOrx+7zuVlj8FQIOMZv317NRjWMUdyq1MsdSOnpt1du904gaUL1snFjkWbZrPt9VdoeMW/cPYUjPQPltuBq77OpFb6IfFBEXhARJSKn1jjugog8JyKnReSpzVzTcmuz1T7gJefz8yCpbHFmMFXOxKRhujkjK7Bncy83w7LvpBIW9fr7cTeJmSgcbzGjaVmgfEMps5ZdY7Mr/OeBDwD/Tx/HvktrPb3J61n2AVvtA+4+n6pchrN/aFw3wQAM3Wm03qeeZzHAeSsafTBiap4J2opv3Dm6n++brgvbYmm9AuW24OqmYLM9bV8CENluf6XFsjU4xaOo0e9amSsezqcBXd/k3vdt9B2z8lW9GovvFbomsnbhlV8wrrKkq1hLnK6spPZnXCCNeagQcqO9A+UbSZm17Bo7FbTVwF+JyLdF5CNrHSgiHxGRp0TkqampqR0anmVf0ataFGVW/G2Jgr5IjWJf/Wa3i+60y7VIdYAcz2TrxM20gVbO+OTb1bFOYN7PHTRSFEHWvJc7YIx9fqy3i20DKbOW3WPdFb6IfAk41OOtj2mtP9vndd6mtb4iIgeBL4rIWa31V3sdqLV+FHgU4NSpU7fq3tqyi/SsFh26BxYupga/T8TFNFXZTcXMtphbt2umq6FKezJSMTiuMeZJM53sorTQbNYUYwVFE7QORlIV0ASCMTjydpzxt6w9DFtwdVOwrsHXWr97sxfRWl9Jf18XkU8DDwM9Db7FshMsjxOoymV44f9l0Y3RDmyutuaQXa7AbY+tnWHTPU69KAftZkxGjutC4Rg0ri7uYsQzRh4x2TeZEcgfMQqiaCjdbbJzLn4RNfkUFI+tmnmzruSCZU+w7WmZIlIAHK11JX38Q8Cvbfd1LZaN4BSPovLjaXOUtO2fhjV7yLr5tLHKFhZp9UUqgaDi1VMzdZIGZDOm0MzNgSjIjqTVxTEmDTVOM3AC4+5J0l4BGqMe2j5P3FhX6tgWXO19NpuW+fdEZAJ4BPgLEflC+voREflcetg48ISIPAt8E/gLrfVfbua6Fsu2UDwGg3dApmQMas/8dIfOn03S3KGEnuV+emWuLV1jWYE2PWrdwDQsyR2AsAq4poJWxWlxlTZpmiKmAKtVTncHnuklIL6ZAJKmbVZ+C7DZLJ1PA5/u8foV4H3p49eAN2/mOhbLjjB6EppfNnIDzTl6G/y0SQoOxlg26O0/X4seCpZrkmbKLD+Hl4Ow1wo/HUNSB3zzuHAYWgum1sBN/fqOa4bqeiZ47Xb1AGh/hXZlcjuoazNvbmqstILFktIp0iocTAuyCmZFu4KONWSl7ny34V9OW0RNpSJtQZ8jW27spas6uDP63tfyS8aX3y42y41C8Q4jHJckZpXvDaSNYnJmRR83TcDVL6WVybHJ2gGbeXOTY6UVLLc8Gyn577RInH4OcMELYOH1ZQHaVI5Bq7QzVirR0EEDXmrz2w3Q02boWhkD6mbSzloZ0x+2eX0D30hSOePuMXXvGNLxuT44Ak5XQVS4kEo+p20a3QJkh8A/Zu6PiqF4u/l8uGDcPLlDJrjbbnZiM29uWqzBt9zS3HCP1cETpnFKuLygyl3MfGkbcy9nHqt4Ma3TDRYX+Zlh8zxpmraBs8+b3UNcM6toP2sMa18yB9DbHdS9o0gnHzdvxlNIv6eKoToBI/dDdgwWXjPjzh00E5Lrr8ixXzFZ2sybmxpr8C23Njda8n/kbUYTvj6ZemlcFlsIauMa0cqkM4qXVrB6EDVBh6CaqYxBwUgZJFHaVcs3sQIVGq15tHnPz0G0PNun39iAa3z0Ydmcy00DrW5g5IzbzcZrl814/Lxx7cudRkit8jqMfldPY24zb24trMG33NrcYMm/UzyKKhyGqAGkTVLajULcjPHzJxHMV9OMl7SFIonxzbv+oiBb0zGB0bYswbF3mYtcfgLmX1nFn99t4F2jVNlpUL6cxOTQjz5gXD3h/GLvACcNyMYN4/Mv3b34scygcdWE8zh3/EgfN9Nys2ODtpZbm82U/IcLxu5mRyF/yGS1tMqLRtvLQP6omQxUSMcl0/HpJ2mla7xClsApHsW57x/BHT9qzquaaVPxblI3jZs1Kp9rMf8SXH8aZs8Cbm/V0aF7FjN0NnovLLcEdoVvubXZTMl/3MAUOaVG0hszK34/Z3YAmWEjx+C4oLpEyNCmcCsoYDpMOUuF2ujyjU8/Z1bhfsHINdcnU/G2Nv6ipPNaKGXcSQBhGVW53Lua2Mof7Guswbfc0myq5N/NmiCviozbRaeSBG66Es8fgsaX0gKmZWmYOjK7AZ0Yf3owZLpFnf1DlJszOfK5VKJKpS4XnfbZ7WQBJUAEiWNSJ31Mhs0KxBRhiYahu83k0CNGYeUPLNbgW255bjjwWDxmJAeihcXGH5kDUDhoVsuzL6SL+l6FVHpRs8YJTEC2ftXIFDRnjE++MWk+64gx6K0ZI27m+KYdYWfHEJtVf+EoNANznKmYMobeDYwv3k8F0LReNUZhg7D7G2vwLZbVaFfeBkeWukBGTy5m/7iB2QHobpcOmFU3oB0j1VCfXNpA3RsxOwbVzsv3zWQgYnYHTmBW6lqlvWjrpun48fdCbgQufjHtv1s3O4GwnMohLJjf1i9v6YEN2losq7Bmu8XWnJkEsgcwKZvL9W5S2QIRyI8vbaDeFj4TD1Cm0AkxQeDMsPnxcqnxj82KPxg0x089BfmDcN//BKU7TfxABPxBc47yq2YHMXpyx++XZe9jV/gWyxqs6gLJDJvV9eBxs+qPm2mQt11ZmwZZ3VSqoN1AXTDFV80piGPj33c8k+Z527uMm6h2HZrTJu0TFoXc3IzZEUyfwbnjR1CZQWPYdZIGexvmfJkh65e39MQafIvlRmhn/7g5k+5Yuwy1yTTjJm9W5Plxk9UTls1rScMEavNj5rj6NbOS9/Od4KnKHzT5+Y0p487x8osyDl7aiL3tn2/XGIikuvYs5uFbLD2wBt9iuQFWZLyMvNEY+cIRY4Db+AMQ+jj3/eOlMgX5MTj+g70rW+/7R6jhe+Dcn6aB34wx9o5jdHfa/vn2LqM75dPm1VvWwBp8i+UGWZHnfv7zaxrg9TJkVujWHH2X8dkrle4aSqZzVds/b9sKWjaINfiWfc1GlDTXZRMGuKfIWzIDx3/IuH56jM/m1Vs2ijX4ln3LDStprsKmDPBqIm/1a2vq3Ni8estG2JTBF5HfAn4MCIFXgZ/WWs/3OO69wL/H5K/9J631b2zmuhbLlnCjSpprcMMG+AZF3pazpTsWyy3HZvPwvwi8SWt9EngF+N+XHyAiLvA7wI8A9wM/KSL3b/K6FsvmaefSd7MBI6sql1HnP486+9/M78rlGx/LZkTeusbDxJfNTqV7x7KZcVluKTZl8LXWf6V1pxXQk8CxHoc9DJzTWr+mtQ6BPwLev5nrWixbwiaM7JYb19GTxt8f1U1qZbu71EYKqLp3LCK26bhlBVtZafszwOd7vH4UuNT1fCJ9rSci8hEReUpEnpqamtrC4Vksy9iMkd1i47pmVW+/bHLHYrn1WdeHLyJfAg71eOtjWuvPpsd8DNNg8w97naLHa6uKe2utHwUeBTh16tQ6IuAWy42zqSDrFvncl49nUwFYm5dvWYd1Db7W+t1rvS8iHwZ+FPhBrXUvAz0B3Nb1/BhwZSODtFi2ixs2snvRuNq8fMs6bMqlk2bf/BLw41rr+iqHfQu4W0TuEJEA+BDw2Gaua7HsOm13UH0K5r4D158xwmX5XpvhnWFL3EKWW5rN5uF/HMgAXxRTTv6k1vpnReQIJv3yfVrrWEQ+CnwBk5b5Sa31C5u8rsWyqzjFo6iRB+DiX6WVsAOmEnb2BVT+4K4ZWZuXb1mLTRl8rfVdq7x+BXhf1/PPAZ/bzLUslj1H/RqU7lrq1onqm8rjt1i2E6uHb7HcKDYrxnKTYQ2+xXKjbEGxlMWyk1iDb7HcKFtRLGWx7CDW4FssN4jNirHcbFi1TItlE9isGMvNhF3hWywWyz7BGnyLxWLZJ1iDb7FYLPsEa/AtFotln2ANvsVisewTpLfA5d5ARKaA12/go6PA9BYPZyux49sce318sPfHaMe3efbqGG/XWo/1emNPG/wbRUSe0lqf2u1xrIYd3+bY6+ODvT9GO77NczOMcTnWpWOxWCz7BGvwLRaLZZ9wqxr8R3d7AOtgx7c59vr4YO+P0Y5v89wMY1zCLenDt1gsFstKbtUVvsVisViWYQ2+xWKx7BNueoMvIr8lImdF5IyIfFpEhlY57r0i8rKInBORX97hMX5QRF4QESUiq6ZxicgFEXlORE6LyFN7cHy7cg9FZEREvigi30l/9+wwstP3b737IYb/kL5/RkQe2u4x3cAY3yki5fSenRaR/2MHx/ZJEbkuIs+v8v5euH/rjXHX7t8NobW+qX+AHwK89PFvAr/Z4xgXeBV4AxAAzwL37+AY3wjcC3wFOLXGcReA0V24h+uObzfvIfB/Ab+cPv7lXv+Nd/r+9XM/MH2dPw8I8Fbgb3f4v2s/Y3wn8Oc7/W8uvfb3Aw8Bz6/y/q7evz7HuGv370Z+bvoVvtb6r7TWcfr0SeBYj8MeBs5prV/TWofAHwHv38ExvqS1fnmnrrdR+hzfbt7D9wO/nz7+feAndui6a9HP/Xg/8Afa8CQwJCKH99gYdw2t9VeB2TUO2e37188YbypueoO/jJ/BrAiWcxS41PV8In1tr6GBvxKRb4vIR3Z7MMvYzXs4rrW+CpD+PrjKcTt5//q5H7v9767f6z8iIs+KyOdF5IGdGVpf7Pb965e9ev9WcFN0vBKRLwGHerz1Ma31Z9NjPgbEwB/2OkWP17Y0H7WfMfbB27TWV0TkIPBFETmbrjD2wvi29R6uNb4NnGbb7l8P+rkf2/7vbh36uf7TGO2Vqoi8D/gMcPd2D6xPdvv+9cNevn8ruCkMvtb63Wu9LyIfBn4U+EGdOtaWMQHc1vX8GHBl60a4/hj7PMeV9Pd1Efk0Zku+JQZrC8a3rfdwrfGJyKSIHNZaX0239NdXOce23b8e9HM/tv3f3Tqse32t9ULX48+JyO+KyKjWei+Igu32/VuXPX7/VnDTu3RE5L3ALwE/rrWur3LYt4C7ReQOEQmADwGP7dQY+0FECiJSbD/GBKN7ZgbsErt5Dx8DPpw+/jCwYkeyC/evn/vxGPBTabbJW4Fy2zW1Q6w7RhE5JCKSPn4YYxNmdnCMa7Hb929d9vj9W8luR403+wOcw/j5Tqc/n0hfPwJ8ruu49wGvYLIWPrbDY/x7mNVKC5gEvrB8jJhMimfTnxd2coz9jG837yFwAPhr4Dvp75G9cP963Q/gZ4GfTR8L8Dvp+8+xRobWLo7xo+n9ehaT9PB9Ozi2/w5cBaL0398/3YP3b70x7tr9u5EfK61gsVgs+4Sb3qVjsVgslv6wBt9isVj2CdbgWywWyz7BGnyLxWLZJ1iDb7FYLPsEa/AtFotln2ANvsVisewT/n8wzTsISJNbNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# write your code here\n",
    "\n",
    "from sklearn.cluster import KMeans \n",
    "k = 3\n",
    "k_means3 = KMeans(init = \"k-means++\", n_clusters = k, n_init = 12)\n",
    "k_means3.fit(X)\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "for k, col in zip(range(k), colors):\n",
    "    my_members = (k_means3.labels_ == k)\n",
    "    plt.scatter(X[my_members, 0], X[my_members, 1],  c=col, marker=u'o', alpha=0.5)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Double-click __here__ for the solution.\n",
    "\n",
    "<!-- Your answer is below:\n",
    "\n",
    "\n",
    "from sklearn.cluster import KMeans \n",
    "k = 3\n",
    "k_means3 = KMeans(init = \"k-means++\", n_clusters = k, n_init = 12)\n",
    "k_means3.fit(X)\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "for k, col in zip(range(k), colors):\n",
    "    my_members = (k_means3.labels_ == k)\n",
    "    plt.scatter(X[my_members, 0], X[my_members, 1],  c=col, marker=u'o', alpha=0.5)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "<h1 align=center> Weather Station Clustering using DBSCAN & scikit-learn </h1>\n",
    "<hr>\n",
    "\n",
    "DBSCAN is specially very good for tasks like class identification on a spatial context. The wonderful attribute of DBSCAN algorithm is that it can find out any arbitrary shape cluster without getting affected by noise. For example, this following example cluster the location of weather stations in Canada.\n",
    "<Click 1>\n",
    "DBSCAN can be used here, for instance, to find the group of stations which show the same weather condition. As you can see, it not only finds different arbitrary shaped clusters, can find the denser part of data-centered samples by ignoring less-dense areas or noises.\n",
    "\n",
    "let's start playing with the data. We will be working according to the following workflow: </font>\n",
    "1. Loading data\n",
    "- Overview data\n",
    "- Data cleaning\n",
    "- Data selection\n",
    "- Clusteing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### About the dataset\n",
    "\n",
    "\t\t\n",
    "<h4 align = \"center\">\n",
    "Environment Canada    \n",
    "Monthly Values for July - 2015\t\n",
    "</h4>\n",
    "<html>\n",
    "<head>\n",
    "<style>\n",
    "table {\n",
    "    font-family: arial, sans-serif;\n",
    "    border-collapse: collapse;\n",
    "    width: 100%;\n",
    "}\n",
    "\n",
    "td, th {\n",
    "    border: 1px solid #dddddd;\n",
    "    text-align: left;\n",
    "    padding: 8px;\n",
    "}\n",
    "\n",
    "tr:nth-child(even) {\n",
    "    background-color: #dddddd;\n",
    "}\n",
    "</style>\n",
    "</head>\n",
    "<body>\n",
    "\n",
    "<table>\n",
    "  <tr>\n",
    "    <th>Name in the table</th>\n",
    "    <th>Meaning</th>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"green\"><strong>Stn_Name</font></td>\n",
    "    <td><font color = \"green\"><strong>Station Name</font</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"green\"><strong>Lat</font></td>\n",
    "    <td><font color = \"green\"><strong>Latitude (North+, degrees)</font></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"green\"><strong>Long</font></td>\n",
    "    <td><font color = \"green\"><strong>Longitude (West - , degrees)</font></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Prov</td>\n",
    "    <td>Province</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Tm</td>\n",
    "    <td>Mean Temperature (°C)</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwTm</td>\n",
    "    <td>Days without Valid Mean Temperature</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>D</td>\n",
    "    <td>Mean Temperature difference from Normal (1981-2010) (°C)</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"black\">Tx</font></td>\n",
    "    <td><font color = \"black\">Highest Monthly Maximum Temperature (°C)</font></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwTx</td>\n",
    "    <td>Days without Valid Maximum Temperature</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"black\">Tn</font></td>\n",
    "    <td><font color = \"black\">Lowest Monthly Minimum Temperature (°C)</font></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwTn</td>\n",
    "    <td>Days without Valid Minimum Temperature</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>S</td>\n",
    "    <td>Snowfall (cm)</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwS</td>\n",
    "    <td>Days without Valid Snowfall</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>S%N</td>\n",
    "    <td>Percent of Normal (1981-2010) Snowfall</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td><font color = \"green\"><strong>P</font></td>\n",
    "    <td><font color = \"green\"><strong>Total Precipitation (mm)</font></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwP</td>\n",
    "    <td>Days without Valid Precipitation</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>P%N</td>\n",
    "    <td>Percent of Normal (1981-2010) Precipitation</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>S_G</td>\n",
    "    <td>Snow on the ground at the end of the month (cm)</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Pd</td>\n",
    "    <td>Number of days with Precipitation 1.0 mm or more</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>BS</td>\n",
    "    <td>Bright Sunshine (hours)</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>DwBS</td>\n",
    "    <td>Days without Valid Bright Sunshine</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>BS%</td>\n",
    "    <td>Percent of Normal (1981-2010) Bright Sunshine</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>HDD</td>\n",
    "    <td>Degree Days below 18 °C</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>CDD</td>\n",
    "    <td>Degree Days above 18 °C</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Stn_No</td>\n",
    "    <td>Climate station identifier (first 3 digits indicate   drainage basin, last 4 characters are for sorting alphabetically).</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>NA</td>\n",
    "    <td>Not Available</td>\n",
    "  </tr>\n",
    "\n",
    "\n",
    "</table>\n",
    "\n",
    "</body>\n",
    "</html>\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-Download data\n",
    "To download the data, we will use **`!wget`**. To download the data, we will use `!wget` to download it from IBM Object Storage.  \n",
    "__Did you know?__ When it comes to Machine Learning, you will likely be working with large datasets. As a business, where can you host your data? IBM is offering a unique opportunity for businesses, with 10 Tb of IBM Cloud Object Storage: [Sign up now for free](http://cocl.us/ML0101EN-IBM-Offer-CC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'wget' n'est pas reconnu en tant que commande interne\n",
      "ou externe, un programme ex‚cutable ou un fichier de commandes.\n"
     ]
    }
   ],
   "source": [
    "!wget -O weather-stations20140101-20141231.csv https://s3-api.us-geo.objectstorage.softlayer.net/cf-courses-data/CognitiveClass/ML0101ENv3/labs/weather-stations20140101-20141231.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2- Load the dataset\n",
    "We will import the .csv then we creates the columns for year, month and day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Stn_Name</th>\n",
       "      <th>Lat</th>\n",
       "      <th>Long</th>\n",
       "      <th>Prov</th>\n",
       "      <th>Tm</th>\n",
       "      <th>DwTm</th>\n",
       "      <th>D</th>\n",
       "      <th>Tx</th>\n",
       "      <th>DwTx</th>\n",
       "      <th>Tn</th>\n",
       "      <th>...</th>\n",
       "      <th>DwP</th>\n",
       "      <th>P%N</th>\n",
       "      <th>S_G</th>\n",
       "      <th>Pd</th>\n",
       "      <th>BS</th>\n",
       "      <th>DwBS</th>\n",
       "      <th>BS%</th>\n",
       "      <th>HDD</th>\n",
       "      <th>CDD</th>\n",
       "      <th>Stn_No</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMAINUS</td>\n",
       "      <td>48.935</td>\n",
       "      <td>-123.742</td>\n",
       "      <td>BC</td>\n",
       "      <td>8.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>13.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1011500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>COWICHAN LAKE FORESTRY</td>\n",
       "      <td>48.824</td>\n",
       "      <td>-124.133</td>\n",
       "      <td>BC</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>104.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>307.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LAKE COWICHAN</td>\n",
       "      <td>48.829</td>\n",
       "      <td>-124.052</td>\n",
       "      <td>BC</td>\n",
       "      <td>6.8</td>\n",
       "      <td>13.0</td>\n",
       "      <td>2.8</td>\n",
       "      <td>16.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>-2.5</td>\n",
       "      <td>...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>168.1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012055</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DISCOVERY ISLAND</td>\n",
       "      <td>48.425</td>\n",
       "      <td>-123.226</td>\n",
       "      <td>BC</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>12.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1012475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DUNCAN KELVIN CREEK</td>\n",
       "      <td>48.735</td>\n",
       "      <td>-123.728</td>\n",
       "      <td>BC</td>\n",
       "      <td>7.7</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.4</td>\n",
       "      <td>14.5</td>\n",
       "      <td>2.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>267.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012573</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 25 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Stn_Name     Lat     Long Prov   Tm  DwTm    D    Tx  DwTx  \\\n",
       "0               CHEMAINUS  48.935 -123.742   BC  8.2   0.0  NaN  13.5   0.0   \n",
       "1  COWICHAN LAKE FORESTRY  48.824 -124.133   BC  7.0   0.0  3.0  15.0   0.0   \n",
       "2           LAKE COWICHAN  48.829 -124.052   BC  6.8  13.0  2.8  16.0   9.0   \n",
       "3        DISCOVERY ISLAND  48.425 -123.226   BC  NaN   NaN  NaN  12.5   0.0   \n",
       "4     DUNCAN KELVIN CREEK  48.735 -123.728   BC  7.7   2.0  3.4  14.5   2.0   \n",
       "\n",
       "    Tn  ...  DwP    P%N  S_G    Pd  BS  DwBS  BS%    HDD  CDD   Stn_No  \n",
       "0  1.0  ...  0.0    NaN  0.0  12.0 NaN   NaN  NaN  273.3  0.0  1011500  \n",
       "1 -3.0  ...  0.0  104.0  0.0  12.0 NaN   NaN  NaN  307.0  0.0  1012040  \n",
       "2 -2.5  ...  9.0    NaN  NaN  11.0 NaN   NaN  NaN  168.1  0.0  1012055  \n",
       "3  NaN  ...  NaN    NaN  NaN   NaN NaN   NaN  NaN    NaN  NaN  1012475  \n",
       "4 -1.0  ...  2.0    NaN  NaN  11.0 NaN   NaN  NaN  267.7  0.0  1012573  \n",
       "\n",
       "[5 rows x 25 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import csv\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "filename='weather-stations20140101-20141231.csv'\n",
    "\n",
    "#Read csv\n",
    "pdf = pd.read_csv(filename)\n",
    "pdf.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-Cleaning\n",
    "Lets remove rows that dont have any value in the __Tm__ field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Stn_Name</th>\n",
       "      <th>Lat</th>\n",
       "      <th>Long</th>\n",
       "      <th>Prov</th>\n",
       "      <th>Tm</th>\n",
       "      <th>DwTm</th>\n",
       "      <th>D</th>\n",
       "      <th>Tx</th>\n",
       "      <th>DwTx</th>\n",
       "      <th>Tn</th>\n",
       "      <th>...</th>\n",
       "      <th>DwP</th>\n",
       "      <th>P%N</th>\n",
       "      <th>S_G</th>\n",
       "      <th>Pd</th>\n",
       "      <th>BS</th>\n",
       "      <th>DwBS</th>\n",
       "      <th>BS%</th>\n",
       "      <th>HDD</th>\n",
       "      <th>CDD</th>\n",
       "      <th>Stn_No</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CHEMAINUS</td>\n",
       "      <td>48.935</td>\n",
       "      <td>-123.742</td>\n",
       "      <td>BC</td>\n",
       "      <td>8.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>13.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1011500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>COWICHAN LAKE FORESTRY</td>\n",
       "      <td>48.824</td>\n",
       "      <td>-124.133</td>\n",
       "      <td>BC</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>104.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>307.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LAKE COWICHAN</td>\n",
       "      <td>48.829</td>\n",
       "      <td>-124.052</td>\n",
       "      <td>BC</td>\n",
       "      <td>6.8</td>\n",
       "      <td>13.0</td>\n",
       "      <td>2.8</td>\n",
       "      <td>16.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>-2.5</td>\n",
       "      <td>...</td>\n",
       "      <td>9.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>168.1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012055</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DUNCAN KELVIN CREEK</td>\n",
       "      <td>48.735</td>\n",
       "      <td>-123.728</td>\n",
       "      <td>BC</td>\n",
       "      <td>7.7</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.4</td>\n",
       "      <td>14.5</td>\n",
       "      <td>2.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>11.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>267.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ESQUIMALT HARBOUR</td>\n",
       "      <td>48.432</td>\n",
       "      <td>-123.439</td>\n",
       "      <td>BC</td>\n",
       "      <td>8.8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>13.1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.9</td>\n",
       "      <td>...</td>\n",
       "      <td>8.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>12.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>258.6</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1012710</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 25 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Stn_Name     Lat     Long Prov   Tm  DwTm    D    Tx  DwTx  \\\n",
       "0               CHEMAINUS  48.935 -123.742   BC  8.2   0.0  NaN  13.5   0.0   \n",
       "1  COWICHAN LAKE FORESTRY  48.824 -124.133   BC  7.0   0.0  3.0  15.0   0.0   \n",
       "2           LAKE COWICHAN  48.829 -124.052   BC  6.8  13.0  2.8  16.0   9.0   \n",
       "3     DUNCAN KELVIN CREEK  48.735 -123.728   BC  7.7   2.0  3.4  14.5   2.0   \n",
       "4       ESQUIMALT HARBOUR  48.432 -123.439   BC  8.8   0.0  NaN  13.1   0.0   \n",
       "\n",
       "    Tn  ...  DwP    P%N  S_G    Pd  BS  DwBS  BS%    HDD  CDD   Stn_No  \n",
       "0  1.0  ...  0.0    NaN  0.0  12.0 NaN   NaN  NaN  273.3  0.0  1011500  \n",
       "1 -3.0  ...  0.0  104.0  0.0  12.0 NaN   NaN  NaN  307.0  0.0  1012040  \n",
       "2 -2.5  ...  9.0    NaN  NaN  11.0 NaN   NaN  NaN  168.1  0.0  1012055  \n",
       "3 -1.0  ...  2.0    NaN  NaN  11.0 NaN   NaN  NaN  267.7  0.0  1012573  \n",
       "4  1.9  ...  8.0    NaN  NaN  12.0 NaN   NaN  NaN  258.6  0.0  1012710  \n",
       "\n",
       "[5 rows x 25 columns]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf = pdf[pd.notnull(pdf[\"Tm\"])]\n",
    "pdf = pdf.reset_index(drop=True)\n",
    "pdf.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-Visualization\n",
    "Visualization of stations on map using basemap package. The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python. Basemap does not do any plotting on it’s own, but provides the facilities to transform coordinates to a map projections. \n",
    "\n",
    "Please notice that the size of each data points represents the average of maximum temperature for each station in a year. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'mpl_toolkits.basemap'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-6b98a7110c83>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'matplotlib'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'inline'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'figure.figsize'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m14\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'mpl_toolkits.basemap'"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize'] = (14,10)\n",
    "\n",
    "llon=-140\n",
    "ulon=-50\n",
    "llat=40\n",
    "ulat=65\n",
    "\n",
    "pdf = pdf[(pdf['Long'] > llon) & (pdf['Long'] < ulon) & (pdf['Lat'] > llat) &(pdf['Lat'] < ulat)]\n",
    "\n",
    "my_map = Basemap(projection='merc',\n",
    "            resolution = 'l', area_thresh = 1000.0,\n",
    "            llcrnrlon=llon, llcrnrlat=llat, #min longitude (llcrnrlon) and latitude (llcrnrlat)\n",
    "            urcrnrlon=ulon, urcrnrlat=ulat) #max longitude (urcrnrlon) and latitude (urcrnrlat)\n",
    "\n",
    "my_map.drawcoastlines()\n",
    "my_map.drawcountries()\n",
    "# my_map.drawmapboundary()\n",
    "my_map.fillcontinents(color = 'white', alpha = 0.3)\n",
    "my_map.shadedrelief()\n",
    "\n",
    "# To collect data based on stations        \n",
    "\n",
    "xs,ys = my_map(np.asarray(pdf.Long), np.asarray(pdf.Lat))\n",
    "pdf['xm']= xs.tolist()\n",
    "pdf['ym'] =ys.tolist()\n",
    "\n",
    "#Visualization1\n",
    "for index,row in pdf.iterrows():\n",
    "#   x,y = my_map(row.Long, row.Lat)\n",
    "   my_map.plot(row.xm, row.ym,markerfacecolor =([1,0,0]),  marker='o', markersize= 5, alpha = 0.75)\n",
    "#plt.text(x,y,stn)\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5- Clustering of stations based on their location i.e. Lat & Lon\n",
    "\n",
    "__DBSCAN__ form sklearn library can runs DBSCAN clustering from vector array or distance matrix. In our case, we pass it the Numpy array Clus_dataSet to find core samples of high density and expands clusters from them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"None of [Index(['xm', 'ym'], dtype='object')] are in the [columns]\"",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-19-6300e51d4d85>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpreprocessing\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0msklearn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mutils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheck_random_state\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'xm'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'ym'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnan_to_num\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m   2906\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mis_iterator\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2907\u001b[0m                 \u001b[0mkey\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2908\u001b[1;33m             \u001b[0mindexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_listlike_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2909\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2910\u001b[0m         \u001b[1;31m# take() does not accept boolean indexers\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\indexing.py\u001b[0m in \u001b[0;36m_get_listlike_indexer\u001b[1;34m(self, key, axis, raise_missing)\u001b[0m\n\u001b[0;32m   1252\u001b[0m             \u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnew_indexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_reindex_non_unique\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1253\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1254\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_validate_read_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mraise_missing\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1255\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1256\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\indexing.py\u001b[0m in \u001b[0;36m_validate_read_indexer\u001b[1;34m(self, key, indexer, axis, raise_missing)\u001b[0m\n\u001b[0;32m   1296\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mmissing\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1297\u001b[0m                 \u001b[0maxis_name\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_axis_name\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1298\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mf\"None of [{key}] are in the [{axis_name}]\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1299\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1300\u001b[0m             \u001b[1;31m# We (temporarily) allow for some missing keys with .loc, except in\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: \"None of [Index(['xm', 'ym'], dtype='object')] are in the [columns]\""
     ]
    }
   ],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "import sklearn.utils\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "sklearn.utils.check_random_state(1000)\n",
    "Clus_dataSet = pdf[['xm','ym']]\n",
    "Clus_dataSet = np.nan_to_num(Clus_dataSet)\n",
    "Clus_dataSet = StandardScaler().fit_transform(Clus_dataSet)\n",
    "\n",
    "# Compute DBSCAN\n",
    "db = DBSCAN(eps=0.15, min_samples=10).fit(Clus_dataSet)\n",
    "core_samples_mask = np.zeros_like(db.labels_, dtype=bool)\n",
    "core_samples_mask[db.core_sample_indices_] = True\n",
    "labels = db.labels_\n",
    "pdf[\"Clus_Db\"]=labels\n",
    "\n",
    "realClusterNum=len(set(labels)) - (1 if -1 in labels else 0)\n",
    "clusterNum = len(set(labels)) \n",
    "\n",
    "\n",
    "# A sample of clusters\n",
    "pdf[[\"Stn_Name\",\"Tx\",\"Tm\",\"Clus_Db\"]].head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see for outliers, the cluster label is -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{-1, 0, 1, 2}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6- Visualization of clusters based on location\n",
    "Now, we can visualize the clusters using basemap:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'mpl_toolkits.basemap'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-22-814b2feba22e>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'matplotlib'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'inline'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'figure.figsize'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m14\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'mpl_toolkits.basemap'"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize'] = (14,10)\n",
    "\n",
    "my_map = Basemap(projection='merc',\n",
    "            resolution = 'l', area_thresh = 1000.0,\n",
    "            llcrnrlon=llon, llcrnrlat=llat, #min longitude (llcrnrlon) and latitude (llcrnrlat)\n",
    "            urcrnrlon=ulon, urcrnrlat=ulat) #max longitude (urcrnrlon) and latitude (urcrnrlat)\n",
    "\n",
    "my_map.drawcoastlines()\n",
    "my_map.drawcountries()\n",
    "#my_map.drawmapboundary()\n",
    "my_map.fillcontinents(color = 'white', alpha = 0.3)\n",
    "my_map.shadedrelief()\n",
    "\n",
    "# To create a color map\n",
    "colors = plt.get_cmap('jet')(np.linspace(0.0, 1.0, clusterNum))\n",
    "\n",
    "\n",
    "\n",
    "#Visualization1\n",
    "for clust_number in set(labels):\n",
    "    c=(([0.4,0.4,0.4]) if clust_number == -1 else colors[np.int(clust_number)])\n",
    "    clust_set = pdf[pdf.Clus_Db == clust_number]                    \n",
    "    my_map.scatter(clust_set.xm, clust_set.ym, color =c,  marker='o', s= 20, alpha = 0.85)\n",
    "    if clust_number != -1:\n",
    "        cenx=np.mean(clust_set.xm) \n",
    "        ceny=np.mean(clust_set.ym) \n",
    "        plt.text(cenx,ceny,str(clust_number), fontsize=25, color='red',)\n",
    "        print (\"Cluster \"+str(clust_number)+', Avg Temp: '+ str(np.mean(clust_set.Tm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7- Clustering of stations based on their location, mean, max, and min Temperature\n",
    "In this section we re-run DBSCAN, but this time on a 5-dimensional dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"['xm', 'ym'] not in index\"",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-23-544e7ecd4e3f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpreprocessing\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0msklearn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mutils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheck_random_state\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'xm'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'ym'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'Tx'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'Tm'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'Tn'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnan_to_num\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m   2906\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mis_iterator\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2907\u001b[0m                 \u001b[0mkey\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2908\u001b[1;33m             \u001b[0mindexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_listlike_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2909\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2910\u001b[0m         \u001b[1;31m# take() does not accept boolean indexers\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\indexing.py\u001b[0m in \u001b[0;36m_get_listlike_indexer\u001b[1;34m(self, key, axis, raise_missing)\u001b[0m\n\u001b[0;32m   1252\u001b[0m             \u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnew_indexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_reindex_non_unique\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1253\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1254\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_validate_read_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mraise_missing\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1255\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mkeyarr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1256\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\pandas\\core\\indexing.py\u001b[0m in \u001b[0;36m_validate_read_indexer\u001b[1;34m(self, key, indexer, axis, raise_missing)\u001b[0m\n\u001b[0;32m   1302\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1303\u001b[0m                 \u001b[0mnot_found\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mset\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mset\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1304\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mf\"{not_found} not in index\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1305\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1306\u001b[0m             \u001b[1;31m# we skip the warning on Categorical\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: \"['xm', 'ym'] not in index\""
     ]
    }
   ],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "import sklearn.utils\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "sklearn.utils.check_random_state(1000)\n",
    "Clus_dataSet = pdf[['xm','ym','Tx','Tm','Tn']]\n",
    "Clus_dataSet = np.nan_to_num(Clus_dataSet)\n",
    "Clus_dataSet = StandardScaler().fit_transform(Clus_dataSet)\n",
    "\n",
    "# Compute DBSCAN\n",
    "db = DBSCAN(eps=0.3, min_samples=10).fit(Clus_dataSet)\n",
    "core_samples_mask = np.zeros_like(db.labels_, dtype=bool)\n",
    "core_samples_mask[db.core_sample_indices_] = True\n",
    "labels = db.labels_\n",
    "pdf[\"Clus_Db\"]=labels\n",
    "\n",
    "realClusterNum=len(set(labels)) - (1 if -1 in labels else 0)\n",
    "clusterNum = len(set(labels)) \n",
    "\n",
    "\n",
    "# A sample of clusters\n",
    "pdf[[\"Stn_Name\",\"Tx\",\"Tm\",\"Clus_Db\"]].head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8- Visualization of clusters based on location and Temperture\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'mpl_toolkits.basemap'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-24-814b2feba22e>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'matplotlib'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'inline'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'figure.figsize'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m14\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'mpl_toolkits.basemap'"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize'] = (14,10)\n",
    "\n",
    "my_map = Basemap(projection='merc',\n",
    "            resolution = 'l', area_thresh = 1000.0,\n",
    "            llcrnrlon=llon, llcrnrlat=llat, #min longitude (llcrnrlon) and latitude (llcrnrlat)\n",
    "            urcrnrlon=ulon, urcrnrlat=ulat) #max longitude (urcrnrlon) and latitude (urcrnrlat)\n",
    "\n",
    "my_map.drawcoastlines()\n",
    "my_map.drawcountries()\n",
    "#my_map.drawmapboundary()\n",
    "my_map.fillcontinents(color = 'white', alpha = 0.3)\n",
    "my_map.shadedrelief()\n",
    "\n",
    "# To create a color map\n",
    "colors = plt.get_cmap('jet')(np.linspace(0.0, 1.0, clusterNum))\n",
    "\n",
    "\n",
    "\n",
    "#Visualization1\n",
    "for clust_number in set(labels):\n",
    "    c=(([0.4,0.4,0.4]) if clust_number == -1 else colors[np.int(clust_number)])\n",
    "    clust_set = pdf[pdf.Clus_Db == clust_number]                    \n",
    "    my_map.scatter(clust_set.xm, clust_set.ym, color =c,  marker='o', s= 20, alpha = 0.85)\n",
    "    if clust_number != -1:\n",
    "        cenx=np.mean(clust_set.xm) \n",
    "        ceny=np.mean(clust_set.ym) \n",
    "        plt.text(cenx,ceny,str(clust_number), fontsize=25, color='red',)\n",
    "        print (\"Cluster \"+str(clust_number)+', Avg Temp: '+ str(np.mean(clust_set.Tm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Want to learn more?\n",
    "\n",
    "IBM SPSS Modeler is a comprehensive analytics platform that has many machine learning algorithms. It has been designed to bring predictive intelligence to decisions made by individuals, by groups, by systems – by your enterprise as a whole. A free trial is available through this course, available here: [SPSS Modeler](http://cocl.us/ML0101EN-SPSSModeler).\n",
    "\n",
    "Also, you can use Watson Studio to run these notebooks faster with bigger datasets. Watson Studio is IBM's leading cloud solution for data scientists, built by data scientists. With Jupyter notebooks, RStudio, Apache Spark and popular libraries pre-packaged in the cloud, Watson Studio enables data scientists to collaborate on their projects without having to install anything. Join the fast-growing community of Watson Studio users today with a free account at [Watson Studio](https://cocl.us/ML0101EN_DSX)\n",
    "\n",
    "### Thanks for completing this lesson!\n",
    "\n",
    "Notebook created by: <a href = \"https://ca.linkedin.com/in/saeedaghabozorgi\">Saeed Aghabozorgi</a>\n",
    "\n",
    "<hr>\n",
    "Copyright &copy; 2018 [Cognitive Class](https://cocl.us/DX0108EN_CC). This notebook and its source code are released under the terms of the [MIT License](https://bigdatauniversity.com/mit-license/).​"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.8.5"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}