saidmrad/ML0101EN-DBSCAN-Clustering-py-v1.ipynb

## ML0101EN-DBSCAN-Clustering-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",
    "\n",
    "# <center>Density-Based Clustering</center>\n"
   ]
  },
  {
   "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",
    "In this section, the main focus will be manipulating the data and properties of DBSCAN and observing the resulting clustering.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the following libraries:\n",
    "\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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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 matplotlib==3.1 -y\n",
    "# Notice: you maight have to refresh your page and re-run the notebook after installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "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",
    "\n",
    "The function below will generate the data points and requires these inputs:\n",
    "\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>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = createDataPoints([[4,3], [2,-1], [-1,4]] , 1500, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling\n",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 0, ..., 0, 1, 2])"
      ]
     },
     "execution_count": 4,
     "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",
    "\n",
    "Lets Replace all elements with 'True' in core_samples_mask that are in the cluster, 'False' if the points are outliers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True, ...,  True,  True,  True])"
      ]
     },
     "execution_count": 5,
     "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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 6,
     "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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{-1, 0, 1, 2}"
      ]
     },
     "execution_count": 7,
     "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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create colors for the clusters.\n",
    "colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABbaklEQVR4nO39eZRk6XmfBz7vd++NNdfK2reu6g1AV4MEGjVAg6SsJtmUwBZNkDJJgzyWZVk+faQhZVs9PEcay0eyPZ4xRc40TYu0QJiGaY45pGjQECETXACKFCgBLfaCBljV6KW6q6pr66rcM2O9y/fOH9+NyMisyK0yK9fvOScrMyNuxL0RFfm7332X3yuqisfj8Xj2Pma7D8Dj8Xg8W4MXfI/H49kneMH3eDyefYIXfI/H49kneMH3eDyefUK43QewEgcPHtQzZ85s92F4PB7PruHll1+eUNVD/e7b0YJ/5swZXnrppe0+DI/H49k1iMjV5e7zIR2Px+PZJ3jB93g8nn2CF3yPx+PZJ3jB93g8nn3Cjk7aelan3W5z8eJFJiYmOHjwIOfOnaNYLG73YXk8nh2IF/xdzLvvvsvzzz/PzMwMIoKqMjIywnPPPcfp06e3+/A8Hs8Ow4d0dintdpvnn38eay1nzpzhgQce4MyZM1href7554njeLsP0ePx7DC84O9SLl68yMzMDAcOHFh0+4EDB5iZmeHChQvbdGQej2en4gV/lzIxMYGI9L1PRJicnNziI/J4PDsdL/i7lIMHD7Lc8BpVZWxsbIuPyOPx7HS84O9Szp07x8jICFNTU4tun5qaYmRkhMcff3ybjszj8exUvODvUorFIs899xzGGK5cucLVq1e5cuUKxhiee+45CoXCPT931o6ZeOVNbnzpJSZeeZOs7RPAHs9eQHbyTNvz58/rXjFPy9ox0xev0J6cozg2xOi5MwTFexflDnEcc+HCBSYnJxkbG+Pxxx/fkNjX3r3NxZ//HPF0DURAlcLoAOf+7o8wcPrIho/X4/HcX0TkZVU93+8+X4e/BdxPES0UCjzxxBObcpxZO+biz38OzZSBM0e7t7en5rn485/j/D9+lqAQbcq+PB7P1uNDOveZpSI68MARBs4cRTPl4s9/jixOtvsQu0xfvEI8XaN4YHDR7cUDg8TTNaYvXN6mI/N4PJvBpgi+iHxWRO6ISN/ibxF5SkRmReTV/OsfbsZ+dwO7SUTbk3PuCqQfIu5+j8eza9mskM6vAr8I/NoK2/ypqv7AJu1v17CbRLQ4NgTL5XRU3f0ej2fXsikrfFX9CjC16ob7kN0koqPnzlAYHaA9Nb/o9vbUPIXRAUYfP7tNR+bxeDaDrYzhf1xEviEivyci55bbSESeFZGXROSl8fHxLTy8+8NuEtGgWODc3/0RJBBqV96jdvU2tSvvIYFw7u/+iE/Yejy7nE0ryxSRM8D/qap3dfyIyBBgVbUmIs8Av6Cqj6z2nHulLHO3lTpmccL0hcsLJaSPn/Vi7/HsElYqy9wSwe+z7RXgvKpOrLTdXhF88CLq8Xi2hm2vwxeRo8BtVVUR+SgulLRr3L02o2kqKEQcfOLRZZ83GqygKOl8c1Mbszwej6fDpgi+iPwG8BRwUESuA/8IiABU9dPAjwB/W0RSoAl8Sndyi28P9ysc0/u8SaPFzMUroDBy7gxRtbSmfdyv7l2Px7M38dYKK5C1Y176+/8TmumiOvr21DwSyD13nvY+bzRc4eaXXwYL1lpsK2b0g2dRC6VDQ3z0//O3++7jXk5E6zlB+JOJx7M72faQzm6l0zTVazMArmmqduU9pi9cvitMsxpZO+byb3+FqT+/TGG4Svpmi9bEHEExovHeFBontGfrBFGIqnLsez/MmR/6C3c9x3otENZzgthtSWaPx7M2vOCvwGY3TXWEdPzFN5j6xttgLZpZbGYBxQQBNsuwcUo0WMEmKS/99C9z4PEHGXr4RPd5VjsRTb7yFqYQdlfnQw8fX/MJwvvpeDx7Fy/4feiEM+beukFrcpbKyYOYIFi80TqbpjpC2pqaZ+a1K2TNtrvDLoTUrM0gEAgMWZwQlApkWcbX/s7/wAf+zg9TOXqA0XNnVjwRJY0Wf/5zv0lhsIpapT09T1pvIoHh0JOPLdq235XK/biq8Xg8OwMv+EvoDWeoVeYv3aR+7Q7HnvoQheEB4N6apqYvXqE1OcfkK2+RNWMkMKDqmnA7eRRVBIOJQjRJsVmG1ltMf/Md3vj0v6B0YIjC6ADHv+8jfbt3bZYxc/EKox98iMLoALf/zQXimTrt6XlsO6Y9U+Pod32w+zoA1CrjL3yrezXQvDW1a6wgPB7P+vCC30O/cEZxbJBbf/INbvzhS4ydfxQjphvPXk9ooz05RzxTI56tAa6rNW22gcXCrYCmmbvVKqCEAxUKIwNUTh5i7tJN3vjl/5OgXKA5MUNxdJDWnRnSZpvm7WmwysDZI9z4/Rdp3pp04m0tWSth7s3rpPNNTv/wdxEWC8SzNSZeeoP2zDylF4dBlSxOscs5eO4wKwiPx7M+vOD30C+cURge4NQPPMnUq29z+GMf4NCTj91T01RxbIisnQupgBghKEVkrQRQJ+4KYalAOFAmayeotYSlAqilfn2cqa9fQgIhqbUYOHuU6QtXsHGCWkVEyOKUoBhSu/wejRvjmCgEEWySoqqkrZja1fe49oWvceSpb+PWl19BU0v5yCjlowcwQUBzYoY7X71I+cgBSoeGu8e/0lWNr+jxeHYHXvB7WC42boKA0sFhhh45ec/x69FzZyiNDWPCABAn6KouNGPExe5Ti6o6Ec8yxAg2Tmi3YhdiUiUoFygMVSkfHqV1ZwbEcOhjjxINlNHMcvvfXGD8334Lm2YEpQLtqXl3QggMIoKmlvqtCd75374MCNWTB5l48Q2CcoHDTz5G+eAIQw8dpz0zR1pv3lWls5HqH4/Hs714we/hfjlbdlbAR7/727j1L1/BtuPFkZxO4jYMyJptsmabcKhC0mgRVcsURwZoT88TlIpgldb4DGmjhYgQFCOyVowYwRQLRENlF7NPMtqTcwtXFVbRPG+Q1pqIMRQPDCJhgIQBNkm5+UevMPrBs2TtlIf+2vdROjzC1CuXAGX0gw9SvznB7BvXuqt4wFf0eDy7CC/4PfQ6Wy5ttLpXZ8vOCrh+bZypC5dpTc4uDdsvkGYLP87UAUjiGslsHQJD2owhX603bk7SnqnRnpyjcWuCsFJCgoCgECKBIWu1u2GiLpl1341B45TWe9Mk9RbGCDbNMGFI1k4hy3jrs79HYcS9F/FsjfrP/RYA1VOHKYxUGXjgCCe+/6O+osfj2UV4we+hYw988ec/R+3Ke6uGM/rRG88OB8u88+t/hFqlfn3cibDm4ZueuH1nP8tiFWyG4k4ImsLky29Art9Zo008U8dEAWG5SJakmCjEphmkts9BZpAfQlZrYgODKtggI200KQwPMPP6NeKpOcpHDtC4OUnaaAGQzDWonDhI/dod5i7doHxo1B1imnWvPMJKyZWE+ooej2dH4QV/CQOnj3D+Hz97T86WS+PZzdvTzL7+LtUzR2hNztIanwW1SBCgllx0dWWx74eyeOUeCCKCTTJSbXfvNoUIm8V3P7/m+zbuRKOZe4Rmlni2jolCmu9NEeQnj6zVdglgIGu1kSggCENm37iOMYbi2CB3vvYaWSvu7sLGKWc+9dT6XpfHswpWM+rJDIltEZkS1WgEI8HqD/QAXvD70s/ZcjWWlnTGszXuvPAa7al5mhOz2HaMphbUlV12RXczsBZVJ97WWjDGhW9KEUG1SFZr3f0YI4gxLnFsLYiL75cPjyJGUKvYVkwy18BmlsAY5/WTpDRuTVI9NgaixPUmt/7kG5hCSGHE1fdnzRgEbvz+i5x4+iM+ju/ZFFppjWv110jtwsIiNAVOVR+jFA6s8EhPh62ceLXrydoxE6+8yY0vvcTEK2+StRc+eL3Dym2acedrryGhIShFaJKiVkE6DVYs/r5RbM+TKd1YvW0lXfG9+8W4Fb17jCCBwYQBYtzGXdFPEjRzdfy2FUNqiSfmqF8bJ56qMXD6MDZOsO2EeLZOPFMDA8ee+hDpXGNHDWn37F6sZlyrv4aqUgoGul+qyrX6a1jNVn8Sj1/hr5XVyg97Szpb4zNkrZhosEL93TvYdtKNmd83lnvurE8Mv4sghRDJ1IWFFMKBinsydU+axSmg7oTlHoJEAWmzjQkDalfvMPbhRwirJdJmm7BcpHR4BBMErsvXx/E9m0A9mSG1MaVg8Uo+MkVaWY16MsNgYWybjm734Ff4a2BpuGbggSMMnDmKZsrFn/8cWZwsKulM5hu0p+eZevWSq2WH+yv2a6Gzyjey8KUWsfltmcWUCtg0Ja21MMUIJAAFCYKFPIDJ6/mzjPKRUaKBMvFcncqxMYYePE7l2NiC75DvzPVsEontE5ZcdH97i45kd+MFfw30hmt6KR4YJJ6uMX3hcreks3b1NlOvvk3j+rhb2a+0wN5KOjkDEcKBMtFQBQkDNMvQzKJBQFgqktbdH1bx4BAEgCrRQBmJAiQKCaslglKBsFpm5NwD3Way3TCk3bN2rGbMx5NMtW4wH0/et5DJWvcTmdKKzxOZ4v04vD2HD+n0sJxFQCdcs7T0sHRopGsoFhQLvP//+kn++N//b2hNzGz/ir4f6mLzWctV7ki+Eo+GKmia0ZqYQYwweOYYxbFBDp1/P+N/9i1QKIwOUBiqunr9KETTlNLBEWw74ZH/5Blufenley5l9ewstio5up79VKMRAgmpp7MYBCMhBVMi1YTQFKhGI5t2XHsZL/g5K8Xoi2NDJPMNbn755UWlh0GpQPXkoW7YIplvMHD6CLNvXNuul7EqmmXdqw5NMiQIyNoJxQMuJGVTi00zjn73hzAmoHblPdozNYJCiAlDosFKHr8vYsKQsFLkxNMf4cTTH/FD2vcAS5OjHRLb5lr9NR4a+simlEGudz9x1sxLMqe7/ShCwFDhEGcHvt2XZq6RTQnpiMhnReSOiFxY5n4Rkf9BRC6JyDdF5InN2O9msVqMvnRkhKlvvkPt+h2SepO0FaPW0pqeY+7SDYYePQm4k8b4i6+7ssudypIQk2YZaaNN48Y4aTMmKLpJW/HkPCYMOPzxxzCBoTAyQFJr0Lg+TlprUj15CFMIuqv4Tinrie87z8EnHvViv0vpJEeXhkgiUyS1MfVkZsv30zk5GAk5VD7NaOEYQ9FhBsIRQgkpBOVNOab9wGat8H8V+EXg15a5//uBR/KvjwH/NP++6dyLc+NKQz+mL17mX/34/5PGzQlXhpmpi4WHeVdrtcybn/09TnzvE7z63/5vtG5N3Y+XdX+xthuBSuotgnKRuXduAlA6PMLBJx7l6NMfJplr0Lg5QeXYGAc/+n4v7HuQ1ZKj7awJ8eSGG5/Wk4SdTyZppfMEEqGaEQUlCvla1VforI9NEXxV/YqInFlhk08Cv6ZuYvoLIjIiIsdU9dZm7L/DvTo3LueSabOMya9fonFrMq9ZNxBYZ4+QWWdMpvDN/+43+Mb/4/9L2ojvfvLdQP7SNU3RVIhnatTfHad9Z4a0FVM+eoDs9/4MU4gwQUDr1jS1y+9RPX6Q8pFRb428h1gpOZrZhNvNtxcJ/Eqx/dTGTLau00prlMIBxkonCU1h1f24+93Kv5XWeHf+AvV0hkDc4sKkAUPRoe5z+QqdtbNVMfwTQG9g+3p+212CLyLPAs8CnD59es072Mgs1uVcMhu3JmnensYExhla2gxMgCYpncJ6CQ1ZvYm1iqbpmo93R5FaV5pJXq4pzpZ5/r1JbJpRvzFB+dAwpcOjHPmOcxSGB2hOzPDCf/pPANBMKY4O5p2+GUef+hCtOzOoKoMPHuP0DzxJYch3Qu4GqtEIoSmQ2PaicEtsmzSzeYaiwxSCBbFeLuY+277Dm7MvkOnCMJ13axd4dPhJhouHl91PYtvdJGwnlCMSEBARihP4TFPmknFGi8cAX6GzHrZK8Pv1evatY1HVzwCfATh//vyaa106YZnKyUM0bk0uqqRpXB9f0bmxU1LZnJhBk6zbQDT31nVQCCpFbK2JqnabZVFFE1fSaFVRa9FkB8fuV8P2tP9apXF7CrX51UyaEc/UyVoxt5KUI9/1QSZefIO5N69TGB7AFEJsZjFRQNaMufrb/9o1dAWG4nCVC//v3+I7fuk/4/DHz23rS/SsjpGAU9XHuFZ/jVZW696uaikHg4vEHvo3PqU25s3ZF1BViqba3Taxbd6cfYEPH/xE98pg6X46txsJmI8nSW1MJRikndbINCWQkEBCUo1pJnMUwoqv0FkHWyX414FTPb+fBG5u5g7ak3Mk9Vb/SppTh1fs+AyKBc7+2FN89Sd/gXSugYogqmTthKBaRESIBspumEivqFsn9KK4mvu9gIW0ljeLBeJcNzPI2gkmCpl/+ybJXN2FbdR56ZiWIZmtuzN4loExFAbKbrpiZtHE8tWf/AWe+cp/T2Ggsp2vzrMGSuEADw19hHoyQztrkuXiOpdMolikT61Hb1hlsnWdTJNFYg/u5NC2dSZb1zlSeXDRfjor/d6cQGJbroDAtiiYMs1srrsf1YwsGOieHDxrY6sE/wvAT4nIb+KStbObHb8PB8vMXLxCOFjumngBpM02MxevEA0tLzRZO+byb/0JRz7+uOs0zVf4rclZbv3LV9FCgASB6zA10l0NSylCwpC01tjMl7JzyBSrLkxlk5R4ro6mGY04JayUSJsxhUKIGEHCAJIUm1oIcMNbKkU0STGFgGS2zrtf+CoP/8TT2/yiPGvBSEBkirzXfJvUxiS2TSObIWm3FsXPO/SGVZrJHFYzkqyNMQYhQMmwarE2ox7PMh+unvhVVeaTSQwGi6WdNeh0EAqQpC3qyQyFoOxFf41siuCLyG8ATwEHReQ68I+ACEBVPw18EXgGuAQ0gL+xGftddAyI83lZEj3q3L4Sy1XplA6PMP3Nd7CpJZ5e3EmKOFdNtZn7ZS2+9ruRnlCPs9dRNMuIa40FO4Z26uL3WX71k1nntJllYJW03kZFqL97Z5tehGe9LK2TL5oKSdYiswvxc8EsirmDS7LOJuOkNkYl7+sgRTAYDCkJ77Uu0chmCIzLqfVL/FrNmGrfwIhBENK0BVgUm38ShUY2wxszX2UgGuNE9VFGS8e98K/CZlXp/Pgq9yvwk5uxr+VI5huMnDtD/dod59iYi29QLjBy7gzJ3PKr8JVm2R740MNkcUrz1iRz79zEiMFmGbYVkzWc93x3UtVeE/u7UPc+ibgpWYmS1lqotfmErZ4tM4sUI6zNaE/PEQ2UqZ4+vH2H7lkTHb/52XicZjrPQOgG3IgYhgqHmIvHSWyLWjxNaAoolrHiSerJDOVwkGv11ygFVWoSYdWSkeAWC7b78UhtQi2Z4mDpNIGJFiV+gXz/d2hlNYYLh5lu3SIlWTwCgghLhtU2tWSKt2ZfYrB5gLHiSarRiPfJX4Y902lbHBsiqpQ4/vRHaN2ZWeTc2Lg23tfEq1OzP/fWDVqTs1ROHlww/sqJKiU+8Hc+wbUvfI35y7eQQLCtlCy1zlIYUGNchU5nNRyIq9ffa5iAqFLsTtOySZ201kTVLm7oEhBxFg5BsYCNUyQKOf2D37Fth77XuB+DQHqtDtpZnVZWI8laDBVcCCc0BUaLx6il05TNEI1sBlXLneYVoriAqmI1oxqNcqB4gjuty/SrzbAktGzCe81LjBVPUwzK1NJp3p27QD2bQTDEWZOWnSe2LaKgRNs2MAgZLsSYkXRzCaKQ0mauPU4tmWaoMEZoit4nvw97RvA7lTbJbIPKsYUmjOVMvHpr9tUq85duUr92h2NPfYjC8ED3sRIY3vnNP+bmH75EOt9Abe4jL6BWXNI2CgkqJbJ6032+7R4U+zCgMFRm6OETqLXMvHbVzdctGLCCGnVXOuJyHJq5UJcWFGOED/ztH/QJ201iM7xuOqZl88kkIAxEI9xpXgGEUjCAycM1ijIX5yEcMe5Lhcn4WteQT1UJsoDIlGhlNQSDklGQMom6JKulU7Lc8QkXMs0Yb10hMiWsZtSSKQIJCaVAKRjoJm3bWR3Qrtg7tGuxkBIjagiCCBHBEHR98jfLCmKvsGcEfz3zaPvV7BfHBrn1J9/gxh++xNj5RzFiCAfLKBDfmSOZrRONuBOKph0/GifsGqdkvXYKxuTCZ/NxhLv0BBAaCAzEKQgUDwxhopC03iIoFdzQ9FKBrBkjYUDWbKEK5cMjSD50ffjh41QfOMKRv/DB7X41e4LN8LpppTUuz3+DuXh8YU6yKoJhrORsQqKghEmD7qo9ti2KQYXYNqlnM27hQ+YmrKHEpDSyOUCYTydRtaQaY8T0fPx7h0K4xyvQtk0MgqhQDCsoSiubX1jBA5aVSp4FEeP2hZJpRjmodMtFq9GIH4uYs2cEH9Y+j7ZfkrYwPMCpH3iSqVff5vDHPsChJx/DJilv/NN/QdpwY/7ColtB9J1m0jtntt/QkSBvakp2il/yChg3CCVLMrQdu9eVZNSu3qY9MUdYLYFVhh48xoEPPcztf3OBoBjRvD2NZpbC8ICzjp5xfRHlQyPeJnmT2OggEKsZ79YuUEsm85p2V1ffzhokts1cfIfR4nEXs48OMZeMk+Qx91ZaQ0QItUBL5wHBkqF5jN6hpFlMFLiqnUzTRff1JwNCRIS2rVMKBslsQiUappXOE+vKNgwZmUsKS5DX6i+IeT2Z6VYaddjPYxH3lODD2ubRrpSkLR0cZuiRkxx84lFufOmlfDs3DSqZa7jJUIG52yCt3yq+p4STTN1VQe9tO5BgsAyZJW3FkKTdy3YJ3UclbbcxxYjSoVGO/sVvJxwouzkBFgYfPMb8pZtkrZjWnRkQKB0a8TbJm8hGB4HUkxlaWd01FJqFP//ARCTaItW4u5oPTYHB6CCTmQvfWDKyLKZtG7nI5uMwWSzlKTFZlvQI/eoo6q4myLrdudZmrqlxxdU9gEXVUk9miIISYbdEVJlsXycypfvq/Lmb2HOCvxaWs1IAFk1p6mxXPjbmPoxJhimESCDoWlwUrC75a9CFQSQ7SfN7ykk1ySgcHCJ7b9rdFxhMISQarJK12tg4oTBcZeD0YbIkpRAEHH7yMe688BrxTJ3SoWEGzh4jKEU8+jef4bgfYr6p3MsgkN4EbyOdJ7N3f3gDCREMmU1pp3WsphgJmIsnMBIiApnGJMSk3cob8nj9YmHvd9tqKBkZeQg0qxOZEk07T6YphjAXfbOM+LtkriAkts1M+xaVcATFhan6OXLuV9O1fSn4nQRve2p+0RSrToJ36JETTLzyJs1bU2RxSliOKB0cpnZtnKzRWpjvuiZ61H0niTw4oTe52OfHZpPUOX52B5wrIkI6X8cmKZoptau3iQbLNL96kaGHjhNUihTHhgjLRY59zxMce/rD3knzPrEWD5peliZ4E9umldXuWtkKQigRsW1Rz6YxNsDajFRTikEFq5bEtsk0o/eD3E/YF0TZsL6Rb06gM83IshohBQITEVEmoZHH8futloQA0124xFkLyxQnqh9gNr697N5Wuxq6H5VQ282+FPyVErxnf+wpvv4Pf7VbvdO8Pc38n32LqFoCtSCChKbrpbMqfUM9xiV0txPBhaaybKExLTBEA2WS+Z6eBVWydox0EtEGJHAll4e/83EaNyexSYaIoXryMPOXbtCenKV6/OCKDqWee2M5r5teD5oO/RK8RVOhldaIbZNAAoyEWE1JbUJs2xRNlaHCIZSMxLapxdN5lQx5Jf16Prfr/4wrFkOY7ykDazDGEGmFts7Tf9WUr+VVsFhK4QDFoLoolq9qiW0rv3JxMx9WMl3bqqlfW82+FHzon+AdeuQEX/+Hv4pmTvzvfO01V0sehbQnZpEopDQ2RDhYJplr0LwxcW+r9u0We+g2UCECQQBkFA8OQ2YJq2VnF9EN9QiaWmefgCUsl9xtScb8pRuMfvAhRt6/YJW0FodSz72zmgdNh34JXhHDaPEYE63rxLZFZpO8e9Wt8o0YQhMRmioks6TEa4ihby6uhFPIsK4U03bioMuhWFKMFBAVl6MQIbYtAglppjWa6Wx3Xq4lw0jQ7fS9a/9bNPVrO9i3gg93J3gnXnmz67h588svo6pdXx4bJxQDQ9aMSestsjh1K/V+FTm7AauoprktsnWDzQfKtCfnCCtFbJJg0wxS7b5GzTKCYoGBs0excUrj1qQzVSss/vAXDwxSu/Leig6lno1hJFg1/twvwatqyTSlaCpYjZGgmgufEtsmIoa5eJyRwhEa6RzbF4dcut/VjkNJtQ2IKxlNMqYlJJJC3mugXZEOpUAlGOZG/fW+4r3RSqidzL4W/KV0qnda4zNkrXiRCZuEIQbImjGF4SoqQuP6HVSM85PfLQgQBUSVMqMffoh4fJbS4VGmvvE27clZslaMjV2VRFgpkdZbmDDAJinhYJnhR08RlovEuTuoiBCW+1wa58PdPdtHIM62QNW6xCswG98h0zQ3xTOUggoFUyBRt5I3EpHYFjPxbdq2Y1a2m3CVPiJlBsJRmvlJazg6jMUSSEAUlBDMsuK90UqonYwX/B46VTlp4+7/8LAUkdQtmqSE1TISBbRuBVhRNNPd01yV19SHQxUe/onv5d0vfI2Z1666IS8CEgR5p4u7hA6rJQoHh0km56gcP+jGOjZjgnKBcKCEiUJKh0f67Ef72ll4toZWWuNO8zLN1MW9VSG2dYQgb4a2CJZmOk8zqxFJidi2aWV1DIYWjZ5Qzk4rK1sNRTWjmTlrhkwTFKUcDt61ZT/xvpdKqN3Cpgwx3yt0qnfskuaorBkTVIqMPn4Gco+YeKbuSjRFkChw8e3eMGO4g2N8Cq2bk7z1v/wBaaNFe3zWVeCkFhunZM02ai22HRMNVcjm6hRGqqhaGtfHSWoNqicPUT48yuGPP0Yyu9iYbjk7C8/W0I1BAweKxwklIrFNLNatfgkoBlVnPpaLusn7UjS3MFgct99NYu9eR6YJ0+1btLM6iY2Zjm/STGbvqioKTMh8PMlU6wbz8WTuBbRQCdVLpxKqHA7e9Zjdgl/h99Cp3vnzn/tnTH8zpXlnhqAQEpQLHH7yMRrvTVE+doAj3/E4jVuTzF0ymCik+d5U118na8bYJMVEgbML7vytdE4GO+FvR0CtpfneFIc+9gEGHjjimlzilCxPSLsKHktQiCiODjF45giP/M3vJ6k1KY4OUj4yyujjZ2m+N7UmOwvP1rE0Bl2JRmim8xhnYEChZ4VqsaApjWw+F8Od8AHdCO74LRlGIdEYUFKbMNm+SSEtUwmGEBMQEHC7cXnRGMZOJc5ylVCHS2e4PP/qrq3eEd3BoYjz58/rSy+9tOX7zeKEG19+mbd+5YvYNKMwOtD11skabcJqmWiows0vv0w8X6d5cwoJDeUjB2jcmCCZbyBhiOaxcEzenXu/RyAG+QXbWhLJYUDp0AgH/y+PUr9ym8LIAI1bk6BQPjpKWm/Rnqkx9qGHGXnsARrXxvnAf/rDfZOwWZysamfh2TqmWje43bxMKRggtTFT7Ru0s0Z3dWtwV5/OFmEX5Z/uCUPoRnPkVs3Omd+NajQMRgcph4srcURkkVVzpxKqHA5yef7Vu0o6ex+zE6p3RORlVT3f7z6/wu9DUIg4/cyTnHj6I3cJWXdF++5tJAqovzuObSeYUkTz5iSmGFEdPULWThAjxPN1tJ26K4D7jaqbLyvqqmtWwkA632D29WtEFReX1yxDVRBjkDCgfGiEkccecJbRKyRh12Jn4dk6OjFoxTKXjONELshr2+0qRmR7jbzbHUE0IDIFRAyV0Bmq1ZIJVDMCE1EwpbsqcXoTup0Zu7u5escL/gr0E7KB00f4wE/+EC//l/8z9XfvUDo4DCjxXINoqAootp0QDZQYPXcWBCZeeYt4eo60mQ8JsZbcV3b1g1hPvswqYbWEhCHtyTkn/s22e7zp2V8YEhYjgnKRIAw7np+u1FQhnqkRlAoc/vhjoNC4NUn9+jjN8RmytvO49+wM+nWDdmLQnVGDkSmS2RirLmTTMUXY/eGbu+ln65BqjCFw1sumCChJ5nyD1FqsWgIJMRIwVDgE9E/m7oXqnc0acfgJ4BeAAPgVVf2ZJfc/BfwOcDm/6f9Q1f9mM/Z9L3QGn3RX7ufOrFnEsnbMt37pnwOG4ugg0VCVpNYgizOa4zMUhquIMQw+eIygVODkX/koj/3nf5Vv/re/zvS3rmLbCRIYNxA9XtmQx5QiZ8LZTpxgr2GoimZKUBBMaCgfPYAAzdvT2CxDkwwTufm8phBhopDCgQGSuQbVM0cYOHs0F/QSQ48cJ601ufnll51pnMCN33+RyRff4Nzf/RHfRbsDWKkb9FT1MS7N/pmLTyuIhKAxguSCv/fEHvpbPXQas1DXhwAQa5NObb4hIDSFPNF7k0jKJMUWVl2DVuek2kjnSWyboqkgcne9y26o3tmw4ItIAPwS8H3AdeBFEfmCqr62ZNM/VdUf2Oj+Nkrv4JOlSca1iFjHWtlEBpukzL11nWSuTtZ2f1hxZgkrRcJKkdLBYSZfepOhh09QOjJKda5O1kpojc9gwsCtqJcjT/KKCBqYNTtsps02IlA+PIoJDOXDo1QfOMLsxSvEtSaVw6OYUkQ8XaMwWCGZa5DWW2iS8dH/3k2hvPjzn6N+5TYTr7yJWqUwUuXwk49RGB7wXbQ7hLV0g54a+CBXa98kIKSeTmNMlVQTMk26YR13AtjLcfyFvxtLSjObwxBQCEoYAjLNMMY4z/+sSaYJmXHzdOeSCQ6XznCndcXN6EVppvO0szojhaPdQe7L+RjtRDZjhf9R4JKqvgMgIr8JfBJYKvjbTr/BJ7C6FUDvFcHcWzdQq5higeadaRf7RnKvGUCEpNFm7o3rjDx2hsa1cSZffhPT4ypZHHVe8V0CV+1jO97z4OxoUzcEHKuYkgvBJNO1u46vl+rxMQYfPk5QKnZ9cbJ2wui3P4imlsoJV0tfGBsknpwnbbaJp2t88O99qnvCO/+Pn+Xy5/4VjTvTDJ45SunwSHf0o++i3RmsrRv0AOVwkEbi6tEF6VoQO6Oyzmp/L2JgmSsZRYmzNsWggtIisxmJuqlaIoaBcJRSMEicNXlz9gWGCoe673NULDITv8dU+yaD0ZhrPOzjY7RT2QzBPwFc6/n9OvCxPtt9XES+AdwEflpVL27CvtdFv8EnsLKILb0iaE3OMn/pJqPf5gRUczdJha71sREhjZPcEz4fmKJKYXiA409/hJnXrtKeqREbIWu2MVHozMyMK4WUKCAousamtOEsiQ98+GHEGCZefKN/KMgIhUPDfNt/8RNMfeMdBFnkBNqcmOHOVy8y8v4HKB0aBiA8NkZ7ap7S2BBjH36ku21QiCgfGqF6/OCicZFdfBfttrOWeLKRgBOV9/GNyS+hmiHi+kaMBljSPSz2cNeQ5a6Zc0CnGxeUQ6UHaKcN5pN4wXohX7m76VnJolxbaAqMlU5ST6YZKhxiuHBoV7lobobg93M1WvpJegV4QFVrIvIM8M+BR+56FCAizwLPApw+fXoTDm+B5Qaf5Du+S8T6XRFUTh6kfu0OEy+9QTBYxk4mzrtDrRuzFgbuKsFa0mYbYwwHnniYubeud+2Yw2qJ0oEh0noTDUOi4Qq2nZKpuvm4xYjCgUFG3n+aaKjCjd9/kWS2QengMCPvP830xSuLSi+lEFI9dYhj/863M/r4We786YW7TmrlgyMMPXSc9swcab25as1878wAm2a0xmdIGy3CSglNM99Fu8UsTc4GsnI4rRNPzjSlGo3QTOdcly1C2zZQtXs6lr+Yxa/RSISoUjQDCIZGNktGSmYzRIRGOkNoCrn9BLkl9AKCITRFKuHQjq/KWcpmCP514FTP7ydxq/guqjrX8/MXReR/FJGDqjqx9MlU9TPAZ8DV4W/C8XVZ6+CTDv2uCEwQcOypD3H1d74K7RQpRgQiJKpE1SJhtZwnZgNsnFE6NsTBJx6levxgt0EpnqmR1puUD4/SujODdObeKoiB6qlDiDFUT7qKgWioQlprUs9nxw4+dIzmrWlULeUjoww9fILhR07y+E//GLNvXFv2pBYNVDjzY3+R8pHRVWvmO13Htau3mfnWVbKWSwzaOMWUorwiybMV9E3OSgT5wI+VfPET2yKUIqEU8/ocJZ/htidj98sPX3GreyMBiqUajRKYkKn2TTLrBD0Q14HcGdxeCUe6t/djNyRpl7IZgv8i8IiInAVuAJ8CfqJ3AxE5CtxWVRWRj+ICbJObsO91sdrgk6VWAMtdERSGBzjyHY9RuzZB48YdgnKRwUqJxs0J0kYbrMUYQ/nwcHf13GvH3Hxvird//ctEg1Wmvv4WSaNF1k5oj88SVkv5gHBD6fAI8WydpNakMFRFAoNNU6JKidM/8HEUy0P/wfd1u16DQkRrYnbFk1r5yGj/5ql2zMTX32LqlUuAMvbEozzynzzDV/7af+esJUoFUCUarjDy/tO8/kv/3Cdut4CVkrOdT+ZKvviRKSEiDBUOMRePk9gWNh9YTj5F6l5863cqis1FH5a+LiPuCseSUTAlTgy8j2u1ixiJqCdTmLw0EyC1MZnG7kpqiQbspiTtUjYs+KqaishPAX+AK8v8rKpeFJG/ld//aeBHgL8tIinQBD6l29Diu9Lgk9XCGksxYci3/b1/n6u//afc+dprxFPzBMUIE0UMv+8kH/g7P8yJJeP9euv6Rx8/y8Wf/xzVU4eZ+dZVyCxiDEHZif3hJx8D4PaffpPK0QOEA+7KoePlM/2tK3zPb/0jhh4+sei41ntSA5enePW//jXufO0156mjSlCMGHrkBINnj1E6OEzabBOWi90Erk/cbg39krOKRTWjndU5XH6QSjREZtO+vvidmnxVZbR4jHo6Qy2ecmt9dSv9jLjfrnch0q23d6MRIxSb5ytceAsJCaXEWOkkVjNCU6QUDFAwJeaScVJ170WGM1x7dPhJ7rSurDpsZrewKXX4qvpF4ItLbvt0z8+/CPziZuxro/QbfLJaWGM58Tz+9Ec4/vRHmPz6W0y+/BagjH3kUcY+/MiqK9+lK/72TA1NMm78wUtIaIhn6sy9fRPbTjnxl88TDpRp3ZnpCm/abBPP1e963vWe1LJ2zJ//3D9j4pW3CAfLXavjrBkz/edXMMWQM//ev9Ot0uniE7dbwtLkbGpj5pJxrLph37cb71AtjKzo5XKgcJxbjTdpoxiJMCZEsUQqWM0QirmX/G7H1dWLGNx60l0Vad5oVQhKRKYM4sJZjXSO1LbRfGD7aPEYSdZyIxY15vTg4wwXDjNYGFt12MxuYV922q7VCmCt4nn4Y49x+GOPbcpxPPjj39M9Gc29dYPxF75FYTg3weqpmKldvb2s4K50UlvadGbjlNqV26C6yNc+KBeQ0JDUmrTuzNxdrePtj7eEXqvejlWCqhJKAdRNv1LVvpOYGukcV+ZedScNNYgBqymnB84xH08w2b7uVsIS3NPg8Z1IqjGoq8rJUNq2nhfouEq6ZjZPqjGttEY1HKaRzdPqqasvBBUn7FJkMHKf+bUMm9kt7EvBXw/ruSLYDHpPAsWxISb+7PX+G64iuP1OJv2azuK5OvFc/9r+sFQka8U0by8WfG9/vHX0WvWqZi4MIYWuUBdMCZG7h3k0klkuTv+rbrcogLEBlWCYWjrNg8MfodIY5t3aRTJNMLkUOIOx3Vy5s1B5FEoJJSMjRbRTjummfRkJiG2T0cJRZpM7rq6+MOaGue/ikM1qeMFfA9tlDnYv8fjlWK7pbOb1azSuTRANle9+kMDwoycJqwVvf7xN9A4tb6a1rlVCx/elt8W/4+WS2pg3Z/8tSdYiCsoEErpYvaY0slkqMkQzned49X3MtseZS8YJJMIYgyrUs6ntermbhiHAEFINx6hnM+5za0pgpVtdk9oYi2WseJJaOs1QtPvq6teLF/wdzHrj8SuxXNPZ0MPHmXr1EmkzxhTbi2L4CAx/4AE+8v/6m8y9ed3bH28TnaHl482r3KpfohQOdFf2vUSmSCut8c7cK9STaaxaNGu4JrygSiAhqcYkWdxtzBorObET3NCQxO6hBK4xqFhCEzkfHZxVSS9WU0QM0S6tq18vXvB3OJsVUlquxNSEAaMfPEsy36D+7h3iqflulc6hJx/jgz/9YxQGKr4aZ5sxEnCo/ABzyUTe3b14Zd+ZxOSGcyRuMLlNcvMvSzurU8pH/ClZd5VbjUaoBEM07TzW2m6Vym5HMCRZk6KpdG+LTIE0W/z6jIQ99+++uvr14gV/F7AZIaWVSkyjaonHf/pHAdZdbeS5vyztsD1ReR83Gm90ywRVFc3DEpOt6yRZm2JQppXVsJLmA8xNXtWT5nbJpW4NeTkcpGVroBBKgTbORXK307GOUDKstag4S3LNp1+JSDcHspvr6teLF/x9wmr5gI6430u1kef+sJz98Ynq+8lsQj2ZYap9AyOG6fgWcdaglTUYLR4lMCFYSGljNcOqJcmaREGJM0Mf6saom+k85WCQZjZPK51n9zdhuQqdTsVRPZkFUQKNqKVTqFVSEkIpMhiN0raNPZ2kXYoX/H3CZuYDPPeflTpsb9Rf5+zgh3iv+TahKXZDEQZDK6szH08yEB5gVu+ANYgookolGuHR4SephEM9z9fCSIioITIlrG2hrDynYWej3U7bIPe5N2IIJKQajWLVut4DgSPlhykG5T2dpF2KF/z7QLvd5uLFi0xMTHDw4EHOnTtHsbj98cGtLjH13Dur2R9Ptq7fdX8UlAhNgThrMZvcdjeKktmM0IQ8NPQElWhB7K1mJLZNI50l1ZjIlJ0vvO5ewXe9thEqSmQKuegXSfPhL+U8j9HKahSD8p5P0i7FC/4m8+677/L8888zMzPjbJNVGRkZ4bnnntt09897wc+f3R2sZn/cSu/unRAMg+EYt9PLWGtQcaENl+AVXp/5KudGn6ISDfWEi9rEtulGIJJiJAJts7NCO7KmSV2GMB/MnmEkpGJGXONVzlLXy90wknCzuXtOl+eeabfbPP/881hrOXPmDA888ABnzpzBWsvzzz9PHO+NCgjP/ae3w7Yfy9koWCyhRAQmJJQi5WCQSjhMKahiNePy/NdJbdwTLhpkpHAYkQBrM2LbyOXVIKw1zLGM5fgmYAgJpbCM2DvfT0NIORh0A8oxBFLgUOkBCuHi3pKO66Wqza9s5piPJ7G6f4a6e8HfRC5evMjMzAwHDhxYdPuBAweYmZnhwoUL23Rknt1Gb4dtL52KkrHSyb73x1kTEIyEFIMygUTdSVeqSiOe4e3ZV6jH06hmKJZSOEg5GKCQbx+ZItVwNB8WshbuX1WPxZnClc2gE38KlIMRIilTMtXuCW2ocJhKOEIoEZVgyNkkmBJGgm7PQRSUSG3MZPs6zXSeuWSc6/Vv8fbcy32vmPYiXvA3kYmJibsaOzqICJOTW+4I7dmldDpsRYRWVut+iQinqo91K0uW3h+YkGJQWfQ5tJrRsvMktkXTznOn+Q7z6QQz8R2m27fIbMpQ4TCBiQhMgKprwioHgwwEY9zPFfxaaGXOCkQQAlNgKBojMgVAcmfMhPl4wr1+iWjZBuPNd2mkc933omgqtNI6U203quNA8TjlYNCZq+VeRPthpe9j+JvIwYMHWc71WVUZG9tfCSLPxuh02C7n1Njv/nI4yJszLzCfTOaTNZ2BmLVZ7hXvBng7V8gYYwPmknFGi8ecfXIyzWB0kHo6g2BopfMYDBZl9bh+Z5Tg5qIo7ayej4x2E6kU8q7hFpmmRFKiFFQZKh4ikJBGNo+1KYcqpwkpYMmI89dcjUbzeQCOxXOA9/bfqBf8TeTcuXOMjIwwNTW1KKwzNTXFyMgIjz/++DYenWc3sppTY7/7zwx9iItTf0Jsm6gqmXX2AaKGQAIiKZJJQqYZBZNX62QtRAJK4SCnB93ndD6e5PLcK0RSwhhX8tkZINKP9Q9FX+vwFWd6JgSotaTGhXkiKZJoK1/ZB4wUjmCMk7SiKTOV3GSi+S5hj3eOG/PYP7CxH5K4PqSziRSLRZ577jmMMVy5coWrV69y5coVjDE899xzFAqF7T5Ezz6gEg5xbvQpBqIDBIQEJnQGamIoBtWF7wpJ1ibTZFG4yEjQ/SoEVQpBuXvlapcVaFmX1Belk3ReT7hISbRFatu00vl8EIw64zdCktwWQtUyF48DSiAFSsEApWAAkYB6OrOsDbS3VvCsm9OnT/OzP/uzXLhwgcnJScbGxnj88ce92Hu2lEo0xAdGv4vx5lXea1wikIhGOrdglSwBxaBCKRwksylHKw9zqPzAogakxLa64xGn2zdx4rxULDthnLUNRO+usEUwKj2P6vdYV4GjZF33y5S266NVSzObp2BKRFJyVx55/0BsW11b6N55tJVgkHoyTTOZo9Jjo+CtFTwbolAo8MQTT2z3YXj2Ob2Ga1YzAtsg03RhBKCEFEwFE5q7xB4WSkM7VUG3G1dItIXB5L750E+oO52ukv+mWCIpMRSN0bJ1WplLIGuPBcLdA1jEDXnBWUAFJiQyJbI0ITJFRA2WjMHoII1sltTGXSM0qykWNzcgChbKW0UM1XAYxe7YkYX3u2lzUwRfRD4B/AJupu2vqOrPLLlf8vufARrAf6Sqr2zGvj0ez/L0+ukXTYV6OkuqbYSAajiMMWZZsestDY1MkeHiIabbt1B1RmT9V/yudl7E5INVXMLXiFAIylQLo9TiKQJTYKZ9C0NIRoLVrGel775ULYGEWIGicaEok48vLARFUhujWDfURadcNVJWI9UEEIaiQ3fF6wMTcbzyPoyYHTeycCuaNjcs+CISAL8EfB9wHXhRRL6gqq/1bPb9wCP518eAf5p/93g895neap521iTTmEAKq/rIdE8WtYvUkklnuyxuHm5ARKZpLvlOoF0y1yVYQy3mCd58BS9hXgefUI6GOFg6xXw82bU8AOcFBHQTtJEpoSgFKSBisJrmk7mUVNtkpKSaUA4HOTf0FJlNSGybQELuNC/flUDuhG4GCwd2hMD3srRps8PU1BTPP/88P/uzP7spYeHNWOF/FLikqu8AiMhvAp8EegX/k8Cvqcv8vCAiIyJyTFVvbcL+PR7PKnSqeQZX33QRVjPaWZN21qIzFLyV1VzNuihCQJiPBYxtk0zTXPrdKhssgqEUVGhnC86UjXSWUlChlra7oR3JrxoCCl1LCIOgZDTTOXdlIeJONjalHA7xwMC39RXwYlDhWv21HRu6WUqnabNX7ME1bV65coULFy5sSph4MwT/BHCt5/fr3L1677fNCeAuwReRZ4FngR3hPePx7Fca6dzCXNy869aSUQ5GqKeThFIgNEWCPHaeatx1qYxMhYIpYkyASMDR8kOLrihcCaRQMGWspqjqomHqRiJULWnu4a9YDMaViIoBMRSCyrKr9dV6GHYaW9W0uRmC3+8o+5lerLaNu1H1M8BnAM6fP7/7JzF4PLsQqxlX5l7FakbBLHjSZJqS2lZXOENZcFqNpITSRDAUwzKRKXZX1Uu9f8rhICJCatsEJsrj8246l+Sel0CeL3DzaAVDJRoiMiUKpkTbNlZsllqth2EnsVVNm5sh+NeBUz2/nwRu3sM2Ho9nh9CZsrV0Rezm4rYJpUhsG7TSmhNs3ASp0egYAIfLDy6bI+g4dRoCMjLSLMlr6l2SNswbqiJTIgqKCG7mbKYpcdakGo50RzzulWaprWra3IzGqxeBR0TkrIgUgE8BX1iyzReA/1AcTwKzPn7v8excEtvq65ZpNaNtG8S2QTkcBhESG+e2DkNEQYkHh57gYPkkg4Wxu8S+d7DLQOEAR0pnKQVVQimi+fCSjBRDQGBCBqOxbqgjkND59fdYR++VZqmtatrc8ApfVVMR+SngD3BlmZ9V1Ysi8rfy+z8NfBFXknkJV5b5Nza6X4/Hc/9wq+sCiba6tfuK8+VxZZElRgqHAWhk86hmnBh4H4PR3SLfy9LBLlFQ4mDpNLFt0UrnGSkeQ1WZTyYYKIwC0Mzmu8cArs5+LzZLbUXT5qbU4avqF3Gi3nvbp3t+VuAnN2NfHo/n/uNq8ItUgmHX2KQxmU27pZkjhaPdsEo1HKaV1bpx9/l4sjt0fWlIp99gF2f1UEGxDBUOOjMzO9+toR+KDjGXjLtjICHVhEhKO7biZiPc76ZN32nr8XjuordhqyJDJFlMW+tEpshY8SRRcHcopZ7M8F7z7buGrvcmbVcb7NKppult+ApNgdHiMZrJHIrl1MAHd2Qt/W7AC77H4+nL0tLGxLaYbN3o2hUotmtP7AaLXCMy5buGrl+rv8ZDQx/BSHCXmPdu1wnR9J5seuvoC2Glb8WPZ+14wfd4PMvSW9poNWMumSCxbQRhLhnHapZ/WYwYDhQrix6/1Gt+OTFf2hS12+rodwte8D0ez5rotVpwk6MUQ0AoBQpBmUY22x2mstTDprd8crnBLc10nkZrdlHsf7fU0e8WvOB7PJ41UwoHOFx+kGZWI5Qod9wsEdsWJrdcSLIWheDulX4vvWLeSmtcnn91xdi/Z3PwA1A8Hs+6yHSh7t7NjDXdgeFWM7Ke2bCrlU/21uV3BpXstzmzW4kXfI/Hsy76VdqIGIYKhwAh0/iuoevLxd47dflLrwAi4+yP68nMfXgF+xcf0vF4POtiuUobRRktHuNI5SyZTdeUaO1Xl7/4/r1hnbBT8ILv8XjWxWqVNuuJu6+lLt+zeXjB93g862azyibXUpfv2Ty84Hs8nntiM8om11qX79kcvOB7PJ5txTdZbR1e8D0ez7bjm6y2Bl+W6fF4PPsEL/gej8ezT/CC7/F4PPsEL/gej8ezT/CC7/F4PPuEDVXpiMgB4J8BZ4ArwI+p6nSf7a4A80AGpKp6fiP79Xg8Hs/62egK/+8Df6SqjwB/lP++HN+tqh/yYu/xeDzbw0YF/5PA/5r//L8CP7TB5/N4PB7PfWKjgn9EVW8B5N8PL7OdAn8oIi+LyLMrPaGIPCsiL4nIS+Pj4xs8PI/H4/F0WDWGLyJfBo72uesfrGM/36mqN0XkMPAlEXldVb/Sb0NV/QzwGYDz58/rOvbh8Xg8nhVYVfBV9enl7hOR2yJyTFVvicgx4M4yz3Ez/35HRD4PfBToK/gej8fjuT9sNKTzBeCv5z//deB3lm4gIlURGez8DPwl4MIG9+vxeDyedbJRwf8Z4PtE5C3g+/LfEZHjIvLFfJsjwL8WkW8Afwb8rqr+/gb36/F4PJ51sqE6fFWdBL63z+03gWfyn98Bvn0j+/F4PB7PxvGdth6Px7NP8ILv8Xg8+wQv+B6Px7NP8ILv8Xg8+wQ/4tDj2SBqU6jfQpM6ElWhegwx/k/Ls/Pwn0qPZwNoaxq99seQNt3vAGEZTn03Uhrd1mPzeJbiQzoezz2iNnVirwqlAwtfqui1P3Yrf49nB+EF3+O5V+q33Mo+qiy+Paq42+u3tue4PJ5l8ILv8dwjmtQ3dL/Hs9X4GL7Hw/KJ15USshJV0YUngLgGNgZTALVu+006Do9nM/CfJM++Z7nEqx5+Au68snxCtnrM/d6cguZt6MTsbQomQoMSsgnH4RPAns3Ch3Q8+5plE682gzf/d7AJmNCt4E0INkOv/Uvs7BV05hKMPAL1G5C28idUEANRFa58Edu5/V6PYxMSwGpTdP4adup1dP6aTybvY/wK37O/6SReSwcW364W0gbMXnYC3ksWQ3PSiXpSdyv6wZOQJdAcBxFIatCahDf+GfrQD66+Ql/uOKIKtKbc/YOn1v3y/FWDpxe/wvfsa5ZNrGZtJ/g2cwIZliEoQTznxNxm7qSg6r6aE9CecfH7sJJvX4QsWdMK/X4kgH3ZqGcpfoXv2dcsSrz2YhMnlGFx4bYsdiJvExezT+bdqj5tgKYgARSHFz9PYRDSBjr+KhoNLJuIXfY4eu5fN/fpqsGze/GC79nfdBKvSWNxPb3NnID3pl01hTR2IRtTWlj1Z00X2gl7RLlzcojnoH4TmtMQFt2qOqqiZ5/BVA65p7Wpuz1rQeMOlA8uhJGShttP9di6X9parhrWk1T27H684Ht2JZtVvigmhFPf7UIframFO8IiDJ2BtN6Nf5O2QNSFaoJ8XyJuBV2/5Vb+8Zxb9dsYwgGYuwxZ6kJE0YA7iTTH4Rufxj78SSfkN/4034dA7QbUbkL1KAQFCMvIqe++t9d2P64aPLsaL/ieXUffRGRQRMceAzHrPgFIaRQe+sFFJxAtH4LLXwR7IA/jxBDXYa7pVv9xHYLYJW4ViAZd8jZtu22tBTubH50CZkHUycNClz6fr95PQGnE5QYGTrhVvWZw/DuRwZP3Xoe/3NXLBq4aPLubDQm+iPwo8F8BHwA+qqovLbPdJ4BfAALgV1T1ZzayX8/+RW2KvvtHTnCNcUlSE8DsOzD9FgycQG0LJECPPYmMvm9NgikmhMFT3RCHANpZ+aexC7FkrYVa+3jWfW9NuRW+CAychOYdJ9zgcgAoSOjuz5pginm4RvKwDy7k03jPiXyHdgZpY8XmL1dueR3mr7nHDJ5EBk8tNIYte/Vy71cNnt3NRv/HLwB/Ffjl5TYQkQD4JdyQ8+vAiyLyBVV9bYP79uxDdPoNmH7T1cS7W1yHazTofq3fdIJqU7j0O+joo3D6e9dVgtgrsBx+wu2jPQ+zX8iTuJJX6ORfzSkX5gkb0Jp22yP5l+Ylnk13ArAtJ/5q85NA4k4khWEo9IRY4nm4+VVs+TDc+Er38Zo13Mls7HGYuQRz77rcArjbh87A2e/vvt7O1YvOX4P56267wVMumezZd2x0iPm3AERWTP18FLiUDzNHRH4T+CTgBd+zLtSmcOsF90tYdt/TtlsZJzUXJokGoJgLZ9KAxjj61m+jp59GBk+621eI/dvGuAvlJHV3UgnLLmxTPuj2I8bF8jWD3gh5lkKtN0naGz3PTwzdu6RnkySP8ceQiLtaMYW82SuDy7/rXlNYcfkAm7qTxPRb7lgKwxCU8zJRhflr7gro4R9aeF3xPNz5+kIuYv5d1Nfi70u24pruBHCt5/frwMeW21hEngWeBTh9+vT9PTLP7qJ+y4lgbyiiI8I2caIa5GWU3QRq7Bba7/4hGg0uhFk6//YIn22Mw8VfXeiuhTw8Y2D2JVd+afIVPgawsGJadDk6VwC68BVPQVpwz2uMK+9U3ImnfBCm33C/h2X3euM599pbE3mYqPPUCvXb3ZLLu2rxOyQNd/tDP+hDO/uIVRuvROTLInKhz9cn17iPfsv/Zf9KVPUzqnpeVc8fOnRojbvw7Ac0qUNQcWKcJe5GCfI7c/ENCi5h2pzMPW3yFbMUYOZtFwYRA8WRRU1INm3BlS86sS8MLjRbKa7mPoshy9x39O7u2/W/mrtvsrEL+aQtFyZqz7r9JDX3WoIo367T9JXmVg4CJnJfKDTH0fac29ZbOHt6WPXUrqpPb3Af14He7o6TwM0NPqdnHyJRFRWBobMuvNGJi3di6UHRrcLjeRcXR9zPNnNJ1iyPdU+94cRx8BSUx1yH7OQFt61aFwrqDa20pyFpApl7rnta1fdjuefRvLlrHDRxJ7Usdq/VBAsnGwWM5FctnTfJuDuyvILJ1+J7etiKa7kXgUdE5CxwA/gU8BNbsF/PXqNTZqgWRt/nVr5xzdW4d0onG+N0QybG5CvyzFXadIhTd3KYet2t5iVwol+/vXBVAO45rELWwIVvtoq8jDOquBPM7Nvupk7tv+bHJiys8CFf+eNEvzHpErVBacU9+Vr8/cVGyzJ/GPgnwCHgd0XkVVX9yyJyHFd++YyqpiLyU8Af4MoyP6uqFzd85J59x6Iyw/aME7vGHQhLMPaY87OZvZyv+BOwec17N17eweQrdeseg0LT5JUzQb6KDiGtLZRhbjl2IdSPuIYvtU7MBSgMuPi+5mGmrO2EXwJ3//xVtDXhTmxi3IkR67YN8lyBr8Xfd2y0SufzwOf73H4TeKbn9y8CX9zIvjweWNIkNXvVhXAGTuTVM7lvjE2dU6UxboVOtuRZ0rxSRoFgIQZuU3db2gaT5ieF7ULdFYwJIBpaiOULrmInbbljlgLdNJnN6/rDysJ7kjTynMC4Oyl0CErw6I/6hO0+w/9ve3YkK06aypukNKnDfHUhph0UcvEOnRiq5W6x7+6B/MkW4t5R3pUqZvPC9Bsizxm0J4DQncCSzB1jUHKreZu64w2KLu4vAYsOPiy5Cp/qMXci6EzkEgN3XkE30snr2XX4/2nPjkJt6pqrbr3gwhVBxSVq+9SN3zViUG3eyJTkFSxrqKQJSwthG9W8HDNw4t/OFpqa1nUGWBpC2gxS0JBuLiFrOK+e4gF3jGnL9QsUhlzCOqm5n+PaQilrcWjxU3rHzH2HF3zPjkFb065paPpNd4MJwdRcVY7aRXXjzmEyc2GcuWvOqlgzQCCZoycAvjJJI78aSN0KGe2Jl0se+hGWv1Lo+0rW+crX+rRL6v7TmjtuEy30DXQSuFmepLb596DQ/yl9lc6+wgu+Z0fQbRCKezpcwa3W5y67qpz2jAvzRAM95mkG5t522xYGAZvbFhfy1W2nQWo5sjzOT76d6QnzSB7e2SmDQvq8jo6/j4QukdvpT+gIvLXuvmig7zP6Kp39hRd8z5axoqVxp0HILAnDBJG7PTck0/YcvPfiQudoe8aFMmzmVujlQ65uPmvj+gqXrMylU9rYK+JLhDQouquFzjSrnU7H5qE4Bmned5AlLmRTqEJ4Jg/5eMfM/Y4XfM+WsNps1W6DkOkTelCF1owTttr1BbuBtOVcMpO6q2bJMldmmTZdmCOIgNB52i9CnPAv8sPphG7It9e88arswkVbWoe/XvIQVNaA4QddiWqPTbS2pp0/UOP2In8g75i5//D/2577zpr8XDoJ2MLAgnVCEPV44uRliDOXnOWAZm5YiO3Uzue2A1m2UK2StXumVuXCbkJ3W9pksdgveOx0jc5s3ry1o8U+R0I3NOXBH8CEC81W2pp2A1ZM/h5pAlg48Re8cdo+xAu+5/6zltmqnS7atOUmTc1dcWGH9ozbtjDokrfJvLP5nXp9ofs0bTmRt9bZF6u6+nwkbzySPMRj3YlDeqpdCFg2IasJJMmmvx33jXge3vgNbPV4Xot/3LlkIlAag855IGnAjT9FvXHavsP/b3vuO2vxczG9XbRp08XiW1OQFmHkQSdYALPjefgmyVffeQVNp7Qy7RFoU3BVKoVBKB9eGFpSPQZzV3OHzV0k6CtRGHRXPnMz7mQYVnEVR4EL8UDuE1TLO3Ob6Px1ZPjMNh60Z6vxgu+576x1turSUYM0xp0wl3Ox71gClw44S4U0AVp0SzA7lTWdhKzN3PmgNQNBA0qjTuyrR91Vwfzl+/WStxYJ3VVNXHNhMLVuhZ+13RXR3GUYPAPzVxdOjFkb3v0y+si/50M7+4iNerx6PKvTO1u1lz6VImJCZPAU5sD73eqzd7hOp7bcRFAYWbALlnDBJrkj+u4BeWI2cwLXnnUr/qwNtd4RDbscTXN3TO2xi87y2QDinD5n3lpw3BSTdyVH6LU/dtVTnn2BX+F77jv3PFt16RDuTm152nSVNBZchUpHsBZNpF18W2fs4dw1mH57B9XWbxJZCwgWylolWLCaSJru9XZOkB2TuOGBBU983227L/CC79kSloZr+o0XvOsxS08Uqk60k5qzFbBzS/KtuqQbtXeubG6i1vG1JwT2mOh3Xl9UXhD7aMDN2e1tpzWRi/HPX4XSQXT26pr/Tzy7G/8/69kyOqZn62nlvyuuP3QKrv/rnhV+ZzW/pLSy+wS5mVinAkjjnu33IJpAeIDugHS1TvQ7YRwTgOTlru05l99ImzA/cFdvhGfv4WP4nh1Pb1yfyhHXdBWWnRmY5GGKZdPCtmdiVArdebR7DcnfCwPtmitnzdouwV0+4IRe8xNdexLimYX+hqTm3s+ekY8+rr832YuffM8eRqJq7v4YgW0vuFt2LBMIcBeuueWxZvms2Lb7vTt4fBdYJqwHCdycXlOAymE49CE48/256+eM2yaZc/H6LMkrmCTvjRDX96DWz7rd43jB9+wuqsecd07WyoU/T1SKcWIXlvJu2s4JIJ8CteijvsfEHuNeY2vandza066c9fK/cKMdo6rzzy+OLrwvQTF/vyp5R3OaT8VyrNY74dmdeMH37CrEhHD2GdyYwqTH4EwXEpVoHqvHbdepSqG3ZHOvEDrh1szZOktuNjd72cXn5y7nsfzcSrpT5iqB88vvLXvtWCnjXTT3KhsSfBH5URG5KCJWRM6vsN0VEflzEXlVRF7ayD49+xe1qRvM3ZqEE98JQdmtTsMKIJCl7kstTgiL+SM7zVhL/OT3ApKv7jXLh8C08/BVb//BNKDOaiGsuCTu4ClnQ5Et6Uz2Lpp7mo1W6VwA/irwy2vY9rtVdWKD+/PsU5a6bXZX9DYfANL5OW04S2CM85YxuXmaba/w7LsZzVfmK0zZymK32h844d6btOVi+8Wz7gognnfb2RiCyLto7mE2OsT8WwAie+0y2bOT6Lpt2nxUXxYvTHjKWlAYzkMT1q14w0puIWBBVzBHuy+jCLcYzdwKX0zeVdyP/L64BuUjUL+Zm82lLsGrGRx4DKke9XX4e5yt+p9V4A9FRIFfVtXPLLehiDwLPAtw+vTpLTo8z46mfsvZIrSne7xgYrfajwagcnDBKK2RX0QWhvKyw051Ti+dhqxdYHu8IvlrUBZyF32xbqpXWofKITj3N5Cs6Zut9iGr/i+LyJeBo33u+geq+jtr3M93qupNETkMfElEXlfVr/TbMD8ZfAbg/Pnzu3z55dkMtJ2XEwalhdGHqkDDiXo0BNUBV5bYnnUhi86qN6y4ShVSFgRxaVnmblzp9wxswbLswPag7O6XEI5/F3Lkia64++vy/ceqgq+qT290J6p6M/9+R0Q+D3wU6Cv4Hs9dZM3cDCzquVEXOklbE27akwldsnH2ilvZqwUVV5qwqCt3CRLsTm8dCXBin+czsoTF4au8F0ECOPD+RWLv2Z/c97JMEamKyGDnZ+Av4ZK9Hs/aCMp5t2xeUaLqqkk6oYzOyl+B+Xdh9GE4/GGoHAWTd9bKCh/1HT+3VlhYm+VumJ2msm44B9eD0C09Nfn7UoIDH0Ae/Cte7D0bi+GLyA8D/wQ4BPyuiLyqqn9ZRI4Dv6KqzwBHgM/nid0Q+P+p6u9v8Lg9+wgpDqHVo65jNG3mCccEopKL6XdqzMEldsceQwpD6PQbYIfdiSKuAWnPHNve1f5yyc4dgsktE2zK4mPNV++dkYxB2Z3YCsMw9riryBk8hQye9GLvATZepfN54PN9br8JPJP//A7w7RvZj2efUz3m4vPRoBO35ngexig73R44npdn5olbCdxjxLgTQFhyg86tpRv77tSub3nidr35gs7A9d6QU8/Vimbu96jqpoSJwAN/CZNPslKbQv0W1idoPXi3TM8uYJFNchrngz3UhWuGziwkcgFaU0hURUyIHv8OuPQ7uZ9Op/HKuFj/doVxTJj72OSi3b3i6EdnmEnqrmg6LC3BNPnqP4th8DgyeNKVsk6/AbdeyPMfFVTEu2Huc7zge3YFi2yS23Nw+yVXgdMr9ku6RGX0UfTAoy6cowKN93pCIJ1KnV575WUsljvJzw1X8oQLVxkmWojBL+vLn28rkbNNQN3jTOge23EAtTEkeQ9CeQydv+7en+k388MPwdTcEHi17sTpB5jvS7yXjmfX0LVJPngOeegHXRdtawqaEzB/zdXpH3j/4u1PfY8b8F0cdlU+Jp8EVRjKY+MdsTe5mObxcnBCG5RdqKRyxCVBN/Qnk8fge62K+4r9kpNP5RAUBlzIKojcfaawYKsgkXt9h59wx//m/77QZRwU3L6y2HXVhiXvhrmP8YLv2ZVIadSJ/uEP5XF7nCDf+Tr69hfQ1vTi7c78ZTfIOxrIO3NDFoQ1r+IRWfiOWRipaNsu8VscdvmC8pGFebHrO2qQgou3C4vDNItYcjXRnMjn+A64SV9CbvmcuecMinDgA/nK3+Ye9/W8WW3GDTJPatCcdD5EeDfM/YoXfM/uZuoNN9B84BSUx/oO8RATYobPIO//lLsCKI068Y4qzmAtyE3WbJrH+ntKPUVcOEjEhWM6JmXhQI8521rI4/FRPnmrk4dYafvOcPbiMAw94PIVmuVJXHGx++KoO+l13C1tPqi8NZFvE/Vcuai7EkK9G+Y+xQfxPLuX+i0XnigdWHx7VIHmBDr+KhoNLFSnlEbh4R9aGJeoGdx51XnHp828Q1dcNVDlkBPHIJ8P25mWJREkU8B6mrXyRLFN3Ekia7sQy4qCb90+TAiDJ91J5vh3uK/56877vn7T3dfbY9C1Sg4h0AWvHXDbZXlpqnfD3Jd4wffsWpYNS6RNJ9bJPBSG757V2jNXV0ff55Kc89edEEYlt3rPmnlCWLqrYiRwiV+bQXkU0jYk06sfqAnzx4uzJNZKLr7NFR6U+/yXDrh8Q3sGspYb8zh8xlXhvP2FvAHNupW9KfTsL6/Hb8/mfQs233/kjNJ8wnZf4v/XPbsWiap3r5HVunF9AKUxJ5YASaNvdYqYEBk+A3ndegc79boT+OYdJ/xB5ATeBE7w27N5F+9qdEYtkodi8jLMwrALu7ijoO9qP6pAdcHGqjcMIyZEDz/hErRZa+ExQQmOf6crx7Rp7hyaQCDdqwGp9rPG8uwHvOB7di/VY06Mk4YTR3AlmGnLVaNEAwvbRhVX0VO/5YZ/rIJEVTRr5KLZMWzLm5xMPlDF9hP8jnGP4OLw5DnYfIWdJe721gTLl4FCd5JV/b38mE8vCsOoTeHOKy6uj3UhoqDg9t+ehpGHIWm6OL8puIRvJ2Tlwzn7Fi/4nl3Looas1pS7sT3rvg+d7eufo0l9bbU11WNuhW974vQSAHmHb9pafB/k9fKSnxjC/KQzCMmsW90Pn3WJ5SyG8T9nxa5b6akSQu5uFFsufwHuvTj8YZh63W2jqQsJhWU/3GSf4//nPbuaRQ1ZSd2VH069nhuJ9dl+jdUpYkL02JPw1j93JxHNG5swTsQlcKWPprg4pNIR8MFjrhpo8AzMvAHVEwsnoLAEBx6ByddyQbZLHi95OOeYq+YpDDjB7rk6WbWsUowrR83fF2+r4AEv+J49gJiwm4hVm6Jz7y4O88C9zWotH85j9026NfvFodywTV24RATCKl2bYtQJezSAPPiDaOM21K7ffbXRafwKxTWGdSp30jaEBVe+GRYXchAsvjrpm7/ofU9ye4neBLXH4+vwPXsK11373U6IW1MLXyLrCmeoTeHGn8LASSfIYcmJv01dFcypp6B63AlycdAlYaOq23bwATj9tGv6Wu6KQowzhDOdMFEh/wrzMI3kpZs9D+l9rt78RS9+CLlnBfwK37PnWBrmuadwRm+MvDjsQkWdxKhNXFVPezrvdl2SNDUBMnjSPU+/xDK430vDUBpxYZvOnN76zYV4fSeG30fE++YvwMfpPSviPxWePclGwxmLYuRiFoVWaE1B1kZ6HTzFdCt6egV3JWHmzDPuKkIVyvnzB0WYueROJFniHrOMiG/Kic2zr/CfDI8np+Md303+rhAll6i6ZsFdaTvtdzIYfhDGHnOJ11VE3MfpPevBC75nX9Er6ouEtzWdr9abnQ1dDTzGVdt0WGrBvEbBXW47v0r3bCX+U+XZNywV9Y7lgp74CwuhlUV17QZqN1hopGLTYuT9TjzGi7znPrPRmbY/B/y7QAy8DfwNVZ3ps90ngF/AWQb+iqr+zEb26/GsF7WpE/ulop404MoXAQPlg4sfVD4AWOew2WvCtlGxX+bE4ydRee43Gy3L/BLwuKp+G/Am8H9fuoGIBMAvAd8PPAb8uIg8tsH9ejzro1N101spA+73pL4QyllKXlNvDrwfGTy1KSv7RSeeztcSS2eP536wIcFX1T9U7XrEvgCc7LPZR4FLqvqOqsbAbwKf3Mh+PZ71smJnqkR32yQsunsTveNXOvH4SVSe+8xmNl79x8Dv9bn9BHCt5/fr+W19EZFnReQlEXlpfHx8Ew/Ps59ZUbSjsmua2oImptUsEfwkKs/9ZFXBF5Evi8iFPl+f7NnmH+CGc/56v6foc9uy9W6q+hlVPa+q5w8dOrSW1+DxrM6KnakVOPtXNtyduxZWu1rwk6g895NVP8mq+vRK94vIXwd+APhe1aWWfoBb0ff60Z4Ebq7nID2ejbJqZ2ppFN2K8siVOm+9JYLnPrPRKp1PAH8P+Iuq2lhmsxeBR0TkLHAD+BTwExvZr8dzL6xW874VTUzeEsGznWz00/WLQBH4kogAvKCqf0tEjuPKL59R1VREfgr4A1xZ5mdV9eIG9+vx3BM7oTPVN1t5tosNfcJU9eFlbr8JPNPz+xeBL25kXx7PXmInnHg8+w9vj+zxeDz7BC/4Ho/Hs0/wgu/xeDz7BC/4Ho/Hs0+Q/qXzOwMRGQeu3sNDDwITm3w4exX/Xq0N/z6tDf8+rY37+T49oKp9u1Z3tODfKyLykqqe3+7j2A3492pt+Pdpbfj3aW1s1/vkQzoej8ezT/CC7/F4PPuEvSr4n9nuA9hF+Pdqbfj3aW3492ltbMv7tCdj+B6Px+O5m726wvd4PB7PErzgezwezz5hzwq+iPyciLwuIt8Ukc+LyMh2H9NORER+VEQuiogVEV9OtwQR+YSIvCEil0Tk72/38exUROSzInJHRC5s97HsZETklIj8sYh8K/+7+8+2cv97VvBZw4B1DwAXgL8KfGW7D2SnISIB8EvA9wOPAT8uIo9t71HtWH4V+MR2H8QuIAX+b6r6AeBJ4Ce38jO1ZwV/jQPW9z2q+i1VfWO7j2OH8lHgkqq+o6ox8JvAJ1d5zL5EVb8CTK264T5HVW+p6iv5z/PAt1hhxvdms2cFfwnLDVj3eFbiBHCt5/frbOEfp2dvIyJngA8D/3ar9rmrR+yIyJeBo33u+geq+jv5NisNWN8XrOV98vSl33wSX8fs2TAiMgD8NvCfq+rcVu13Vwv+JgxY3xes9j55luU6cKrn95PAzW06Fs8eQUQinNj/uqr+H1u57z0b0ukZsP6DKwxY93hW4kXgERE5KyIF4FPAF7b5mDy7GHHDv/9n4Fuq+vxW73/PCj5uwPogbsD6qyLy6e0+oJ2IiPywiFwHPg78roj8wXYf004hT/r/FPAHuOTab6nqxe09qp2JiPwG8DXgfSJyXUT+5nYf0w7lO4G/BnxPrkuvisgzqz1os/DWCh6Px7NP2MsrfI/H4/H04AXf4/F49gle8D0ej2ef4AXf4/F49gle8D0ej2ef4AXf4/F49gle8D0ej2ef8P8HtwoqDmq2OHQAAAAASUVORK5CYII=\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",
    "\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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABVjklEQVR4nO29eZBl2V3f+fmdu7wt38ulMitr6+qq7q7eKbWkshbTgAAhS8IgWwaHYGxrgAmFsRX2OJgAT0DMKBw4BtvBxGADlhWMbDFmMYQFCEtiERZIjSVQqymVunrvrurqrCUrs3J7+13OmT/OfS9fZr3cKrNyPR9FRuZ7775773vq+p3f+S3fnxhjcDgcDsf+R+30DTgcDodje3AG3+FwOA4IzuA7HA7HAcEZfIfD4TggOIPvcDgcBwR/p29gNUZHR82pU6d2+jYcDodjz/D1r3992hgz1u+1XW3wT506xdNPP73Tt+FwOBx7BhF5faXXXEjH4XA4DgjO4DscDscBwRl8h8PhOCA4g+9wOBwHBGfwHQ6H44Cwq6t0HGszc+FVLn/6KepXJimdHOfUB59k5Oz9O31bDodjF+I8/D3MzIVXefbnf5totkrpxBjRbJVnf/63mbnw6k7fmsPh2IU4g7+HufzppwiHBsgNlxGlyA2XCYcGuPzpp3b61hwOxy7EGfw9TP3KJOFgaclz4WCJ+pXJHbojh8Oxm3EGfw9TOjlONF9f8lw0X6d0cnyH7sjhcOxmnMHfw5z64JNEczXas1WM1rRnq0RzNU598MmdvjWHw7ELcQZ/DzNy9n4e/4kfJBwuU5+YIhwu8/hP/KCr0nE4HH1xZZl7nJGz92+5gXelng7H/sQZ/G1irxjRTqlnODSwpNTT7Rwcjr2PC+lsA3upXt6Vejoc+5ctMfgi8kkRuSkiz67w+rtEZF5Ezmc//8dWXHevsJeMqCv1dDj2L1vl4f8n4L1rHPNlY8wT2c+/2KLr7gn2khF1pZ4Ox/5lSwy+MeZLwMxWnGs/speMqCv1dDj2L9sZw3+niHxDRD4vIo+tdJCIfEREnhaRp6emprbx9u4ee8mIulJPh2P/IsaYrTmRyCngvxljHu/zWgXQxpiaiLwf+AVjzJm1znnu3DmzX2ba7pUqHYfDsbcRka8bY871e21byjKNMQs9f39ORH5ZREaNMdPbcf3dwN2ol3c4HI6NsC0hHRE5IiKS/f227Lq3tuPaDofD4bBsiYcvIr8BvAsYFZEJ4P8EAgBjzMeBHwB+XEQSoAl8yGxVLGkbuFvhmN7zerkAbcBEsQv5OByOu8KWxfDvBrshht/beRoOlojm60RztU0nMnvPm7Yibnz5AgDj3/Yt+Pncuq7h8gIOh2M5Ox7D38v0Nk0B3d+XP/3Upoxr73mv/+l5wkqJuN7i2h89TW6ojMoFPPfLn+HJj/+zvu+/EwmEjSwQbjFxOPYfzuCvQf3KJKUTY0ue20zTVMeQvvwfPw8ieLmQ1q15RISk1gQgaUWYJGXhpQmOv/stnP6B77jtPBtdiDayQDg9HYdjf+IM/hqUTo4TzVa7BhXuvGmqY0ib03M0b85h4sS+oA0oAU+BNsQLDfxSHgk9zv/s/8fggyduM7SrLUT9vPONLBB3a1fjcDh2Fmfw+9BrMCUMaF6z1aO9MfwHf+x9Gz7v5U8/hU5Spv/ihUUjn2j7ojagU1QhRERIG20GH76H9kyVr/6TX6R839EloZWVFiIvF3S9cwk8rv7BX/Lar30B8RRHvuNNS+5npZ3KVu9qHA7H7sCpZS5jubKlnwvAgI7iTXee1q9MUp+YImm28XIBXi4EYfEHEBHEU/iFHGkrIWm0ad6cu01lc6XuXW2wieB2zI0//Qb1iSmi+RqNa7eY+OOnaU4uKmBE83UkDHjmY5/iyz/6r3nmY59i5sKre0oKwuFwrB9n8JfRT9myfN9R8uMjfNsnf5K3fOzDdxzWKJ0cpzVl4/UAylOI79mf0AdPEZQL+OUiKhfQnp4jHCxRODyEKEXajpl/8Q2+9OF/xeVPP8U93/dOdBRz5fe/wo0vXUAVcjRv3CIcLDH99ItEM1V7YU+RRgnNq9Nc+q9fon59mvZslepr12lem75Ntnn48VN7RgrC4XCsHxfSWcbdDGec+uCTvP67T+HlQ+JGC6UUojx0miIIxXvGSKoNTDOi/MAxotkaxhjiWpNXf+O/E83XyI8NIZ4imq3yyq/+ERg48u1nu+Gm2uVJ/GKextUpJBdgtCFZaOCFPiofkNRb3Pjv5xl58/00rt8iabSJFxoMPXKSwvgIALPPXubxn/jBJXmAB3/sfa6ix+HY4ziDv4ytTNIuZ+Ts/TzxM3+Pr/6TX8TM1Ug7PRC+h5cLINWUTx2heOoIuXKR2W++RuvWAsr30O0IgOaNGcqnxskNl5n66vMgMPrWBwGbXB1+7BSzFy9hUoNJYuKFBiZNUTnQ9RQx4BVCZr7xGsr3yI9WSFsRN7/yHIff+Sj5sSHqVybXLQXhKnocjr2DM/jLOPXBJ3n2538b2HyStkPHA751/mUWXp4gmq2hwgBjDCZJbeJWhGi+TuPqLW498zIqF6KKeSRNiWtNkloL8RTK9+isE0k7Im22uf6n54nn6wSDJSoPnaB86ig6Sll4eQK0PbduRWBAFULmX7kKqSGoFGlOzREU80jgMfEHf4lfzJMfG+pO41qrE9hV9Dgcewdn8JfRkQdebzhjLToecBon1C7dYOG16xBrVOAhIhilQAzJsiSpbkXWSHdQNu6fAnPPXiJtRlSvTKKbbWqv30QFCuX7zF28zJHvfBPDj99L/Y1J0lYE6WI3tY5TRClMkhLNVlG+R1RtINrg5UOCSoni8VG+/jOfJJqvo1sxrbkq7ZtzAOQOD5G/+Dq3nnmJt/7sj7qKHodjD+EMfh82o2y5PJ7dmpwhHBpg5huvYlIDcQpKMFpbTz1NIUvirorOjLaATlLmX7piDbkSTJqSJgnaT/FyATPfeA0QcoeHaV6/hUnjxfMkKUZpMAaUAt9DkgSTpKStiLQd07g6TePqNO25GpX7j5FUG6StCAN41SYMl1l45Sov/PvPdENgaTtm/oUrxPN1VC5g+IkH7uj7czgcdw9n8LeQ5fHs+Zfe4NoffI384SFa0/MkrbY1tKnBKGU9fANsRM+oc3x3kTCI8uwCkqSkIjbZqzVePiQ/UiGar5E22tjaT5MZe0GUstWgojBYz784PkzaiqhN3MTP52yCudpEAg8xkFQbdidgDDe/+jzf/qmf4us//UkWXr1KWC4ioU+0UKd5bZqZC6+6sI5jS6lGM0y3rtBOa+S8AUbzJymHIzt9W3sGZ/C3kN54duPGDPPPvY6EPq2ZBZJGGx3F4AkkBlKNWYdjvyKdRUKD8ckWD4OJE9qzVUht7F4CDxX4i+9RYhcNJYivsoXDIIGP+J6Ve8iHgGBS2xSm0wRSMNruDGaevYSfD/FLBdsEdmKU1s050igmHCwx+pYzqDBwcXzHllKNZpioX8STkFCViHWbifpFTvCYM/rrxBn8DbBW+WFvPHv+hSt4uZDSiVHmLr6Oyvk2Jq97TrhVQqVRsvRUvljDru0CoG/bQdhdBjnF4EMnWHjlGmkrwivkMMagWzFBMU/ajmhNLyBGoROb9MUT0mYbHSX4pTwzF14lbccc/xvnELXY1mG0dnF8x5Yy3bqCJyGBygEQSA60fd4Z/PXhGq/WyfIO3N6u1w69HarxfJ32Qo2Fl66i2zFJrbV4MhHwFZIPusnYLSXRLK4mmace+Ejoo8IAb6CABL5N3FabBAMFcocqDJw8bGv2CyGjb3+YwvFDRHNVvHxg4/0iqCBARPBCn9FzD3H500+5zlzHttBOa/gSLnnOl5B2WtuhO9p7OIO/Tvp14IZDA1z+9FPdYzpyB3MvvUH1jZtUX5ogaWSGvpN0DTxUPoBEY1rx4vMbRei/WHSe0oCvIPCs7U9tjD8cGmDowRMMP34KPx+SLNQJBgr4A3kqD57g0NseIm1FVF++yqE3PUAwUKB07zjhYIn8+BDhYIncoQr50UEq9x+jfmVyTw1pd+xdct4AiYmWPJeYiJw3sEN3tPdwIZ11sp7yw5Gz93PP972T8z/7n4kX6osJ1l7iFB2nm78hpRAlGJMuDQ2JLF4z0SA6M/gGPEVSa1K9dINwpIyEPn4xT+HIMPU3prj6R0+jfEX53qMc/a434+VD5l96g+HHTpHLBrV4+ZC0GeEVwq4Xv9WlrI6dZ7uSoxu5zmj+JBP1i6CtZ5+YiNREHM2f2fL72q84g9/DajH60slx5l96g+a1W90mp8KxQww+eM+Sc8w+e5nisUPMPXvp7t5sqjFabs8DLN8x9D40mrQdkdSbRLNVJLRSC14xT/HEGHG9Balm5M0PUDxi/9ENPXqK2YuXbQfvs5dI2xFGGwbuP7qkIc0Nad8/bFdy9E6uo/CZjyYBQzkY40TJJWw3wpaEdETkkyJyU0SeXeF1EZF/KyKviMgFEXnLVlx3K1krRj/8+Cmm/+J5ovkafrlANF9j+i+eZ/jxU0vPc/4Vbj3z8vpq6zfLRsdTajBxtiMw2IRulFJ95SoLL1/FJCkGw/wLV7pvGXzgOAOnxqmcOUH59FHCwQEGTh9h8MF7nHzCPqU3OSoiBCqHJyHTrStrv/kuXaezOIgoDuXuYTA8gmELdsoHjK3y8P8T8IvAr67w+vuAM9nP24F/n/3eNawmEQDwzX/1G7RnazSuz2AwBMU8haOHuP4nz3QnUs1ceJWpv3ye9tTc1lXg3E20AZ1gxCNeqGOMwS/liXsSsNF8nUNPnOEtH/vwDt6oYztppzVCVVry3N1Ijm7kOq5CZ2vYEoNvjPmSiJxa5ZAPAL9q7MT0r4rIkIgcNcZc34rr93Knyo0rxehnzr/C9S/+FbPfvIQRsdo3IqRxQnNyhld//U/QBipnjvHiJz5L6+bc6sbekyVSB7uFtB3jFUJ0K6I5Ncfl330KPxeSHx/irT/2PqeIeYDIeQPEum2NakYnObqVsf3VrtNL55rGGAKVo+gPEnoFV6FzB2xXDP848EbP44nsudsMvoh8BPgIwMmTJzd0kc0oN66kklmfmKL++iTie+goyZKimrTehFQTDg1w88+f5fJvfdHWoWf17yuyC429STR4VosnjWLC4TICxI0WrRev8NV/+ou0Z2sMP3aKyv3HlnyvgFsI9hkrJUcHvcMbirlP1i9xrfEikW4QqiLHig8xXjq95nV6k7DdUA4KEUhNwkJ8kwqHEVGuQmeDbJfB7xfQ7mv5jDGfAD4BcO7cuQ1Zx80oN66kkhlXmzaEUymSTs1b4bHMaJvUkB8fpnltGhEhrjcXRxbuJYxBlG2oKj9wnPFvfZzpp1+kcfkGaZJSf2OKsFIibbYJK0UK4yO0pud56sf+DdFcndyhCsNvuo/5l97gz/7B/4UXBrSm5tCppnLfMR756N/qO4jdsTsphyOc4LElnvzR/JkNhVUm65d4rfoMnvj4kifWbV6rPgPQNforXaf3XJ1rDgQjVKMplPgYo6jFMxSDiqvQ2SDbZfAngN5ylhPAta2+SP3KJBJ4zPzpq0vkguMrjTXf2ympfP4Xf5fGtWmKx0Z55KN/i5kLrxHXsmajfEjajrpLlfgK3U4A29lq2slWf6Rtw2TduvPPvc7881fsTkXAC3x0lNj6+pkq00+/yOi5h5h99hLN6zPkDg9Rn5hi7oUrKN9DAo9otkZQKSBK0bh+i6d/6hMAzujvIcrhyG1G/Gr9uXXH3K81XsQTv7s4qGxxuNZ4cYmX3+86vXTi/KIEwjGayQKaFE3iKnTugO1qvPoM8A+yap13APN3I34vYcDkl79J2ozwK0XSZsTkl7+JhMGa75258Cpv/P5XGHnT/dz3oe9i5E3388bvf4XKmWN4uYC0HaMCr0e1UvBydsB5OFwmrjWXzKbdsxhjNfSTFOIUYwziKZJ6C92Omb34OhOf/0urnqk17el5+7ZUk7YiopkqRhuCYgHl++goJqgUef4Xf3dnP5fjjqhGM1xaOM8Ls0/RTGo0k+qS11dqfIp0A0+W/rvzJCDSaztfS1HMtK8y3bpCNb5FomO00fg9cX/H+tkSD19EfgN4FzAqIhPA/wkEAMaYjwOfA94PvAI0gB/Ziusup9t4KpndlWXPr8JK4SC/mGP48dPULlste/E9wA4a9/Ih4fCAvVCq8StFknprd4R1PGUF1DaJbkVIIYduRmiVgDE0JmfxQh8JAtvw63vd400m/9y+tQCeIL5P6Z7DNK5Nb/peHNvL8jr5RKUsxFMAFPzyio1P1WgGbTQNPY8vAQoPTUqiI3wVMlm/RCOdXzPxW41miHST1EQYI7R1A9tCLsS6xYXpL3CocILjpYedp79OtqpK54fWeN0A/3grrrUaaTvmyLedZf7FN4jm64SDJY686X7Sdrzme1fspJ2Y4q0/+6Nc/vRTvPbrX6B47BCDD99Le2aeyacuUn31OmiNESFtRnYHsFbi9m4jIIGHMXqpWNsdnUswrch6/sbuasSzU7ckTcETdJKCMXbYCljRNgwkBpNGNCdnKR4b3fTHcmwvy2P2pWAQyLx3rfrG3DuLRNEfZD66Sdu0sGN7xP5Pe7xWfYZKMEbBL6+a+J1uXSHvDZBTRaaaV4A0i6Yaq/9HwlTzdeai6wyHxzhWesgZ/jXYV522nUqbo+96ovtce7ZK4eihdb+33yzbThdp48YMs+df4cYXz9O6NY9JNV7BSg0oJXburMEGyjoediZBv60YrE7PViDYyiKBcKhE+fQRABYu3SCuNcmFPu1bC1b62WD1exJtv5MwwMvZBO5b/+WPbc39OO4qvWWX1fgWA/4YQU/gt+hXiLTHw8NPdo+9Wn+u66l3FolKWEbrlPnkZve9Cp8WVdCGW+02w+YoA+HwksTv0uvPUPZHyflFqzdIjtQkGAwigjYJguAxwFx7koX4JnmvzEAw4nTyV2BfGfyNzqPtrS2X0Mbj+7135sKrPPfLn+HK736ZuNoEbdBp5rX4irBcRAU+SaNtNXS0seGMwLPlnM2o7/V3PYJV1ZQU5fuUTx/BLxUACIcGSGoNdDtGfI+gmCepNTAieKWwm+gNKkUq9x9zCdst5m5o3SwP4SjmWYgnERkn5xWBxZj9ZP0SV2oX0GgCyZNqzUR6kTiNGAiy+xDwCVDik5gITYL1CgRNwlxs03i+yjEfvU4tnqGVVil4QxT9CsIc8/ENhuQoYEsy7TkgNml2CQ+DITJNcqZInLaJPaeTvxL7yuBvRMTrtpr9eSt2pqOY+sRU970Az/78bzP9zMuQGvxCjmiu1p06ZVopUVrDGIPn++TGhqynj6ByPu3pBXvB3jzC7ivFv53QJxwoUDhyiObNGcRgJ2tlevnxfJ3KmXs48uTjvP67fw7GYPIBXj7H0CP3gjFECw2OfPtZwp5dk2PzbIXWTTWa4Wr9BarxFCBUgjEMZkkIZyAYYS66Ti2eIVSFJfX4V+oXSEyKISXW8zTTefJemUg3aaYLhKpApJso8bOqGo1CFkMyCKlJuRVNoPDwxQdjX28mc/gqYCAYYT6apBbPIHhoWn0/S6xbeHj4KiQ1sb1/14Xbl31l8GH9Il59k7T3QThcXiIj8MzHPkU4NEB7eg7JBShjbAw7TrpaNkbbCVZpEqG1xivlSWtNTM2OE1QFOxw8rTdJWgkku7B8M1uQvFIBtG0oS6OY6qVr6DjBCwOql693vzO/kGP8rz9GYXyEwYfusfNw44T65RukzcgajzBYdYfluDM2KzNgK2+eoZlW8STAGMNcdJ3UpAyHx7shnNArMBgcoZpME+n6knr8RCek2u5cNRqDpppM4xGS8wJi3SbWbZuwNWl2VC+ZoBOgETwV0tI18qqMEkUzWWAod4RKME4tmSIx7b6fxS4kKaEqoU3a1cvvLRd1YxEXObB6+PUrk4SDS2uKl8sdLz1O0FFMNF+3owF76cTqsYJkyVyN4snDlE4eBkA3I9pT86hcSOWBo0iwy9bZTFu/E3MfuHccow3RzAK6GSHKwyvm0FFC2mwz/MQD3PO9b89GIcLQIydtwlZrBk4fAYFopsqhtzzgBNbuApsdBDLdukKkW/gS2sYoFeBJCMZQT2aWHKuURyUYI+cN0E5rtjwyuoXJjLzOBMxMZs5TYhIdkZoYbQyxad1m6m9H44mHwifWTQTV1b3XJgFUVyhN4SFLzJbOvpMGiY4o+jax3CsFMVG/SKzbS3ZD1Wjp5zwo7DLLs32slqTtd1zx+Chzz72OKFuSeVsljln6u355Mqts6cye1XaYeDtGAoXZopzqpvGyubZpCkqRGykTlIu0Z6rYhjLIjQwgSpEaCAYHKB4ZWZIvyY8NMfSYlVEunRjj0BNnnMTCXWS9GjQr0U5rGDSKxTp5JR6CR2KsZ96ROmilNnwZ64hIN5iPpkl0i9RoQGOQFVpPBENsE6yoNaKYhnoyj0KhSVHaJ/QKNOIFFuIpFCqL1dtFxp4vm8WMAuOR0ibVCQvRFEV/CKVkw93BB4EDa/DXm+DtHFd58ATzL09g4tTOpoVsILjpH5Pv1MD7nWYtY6tXak28Uh4JfUy0g6EdEZtYxo5bNFqj09ROrrq10M1ThCNlQEhqTZJWRDRf4+X/+HkA7vm+dzL77GXqVyapnDnB2Z/6kDPy28BGB4EsD2lI5iVb42l7KLRJCb08Oa9EPZ7t6t94EqLEp5nMoTGkOiY2CfRIEy+Zv5MlZLVOe3ygtWuDDSka68HHpkVoikS6QTkYpZnOE+gCbWrYgt/Fayt8AuUTEJKYtq3VTw0n82cphyMb6g4+CBxYg79agne5MmTHsN38HxdpTc0RVEp4+ZD2XHXt8ketrf5ONgAcrVE530o07BQCKvQxRltdoDgBTxGWClYJNEqyclIhqjZRqg1iB5OjDWk7ZuHlCWqXri8RULv4//xXJ6C2DaxHg6ZDvwRvO613DauHjeFrYgJVQJuEUjDMoIyTmIiZ1gSehGhMj2e/sg69QZOalN4Y/XoxpHiSxxOPFNv7kZcKvoQo3yNJWqQ9lT4KDxGhrZt44iMofOUzGB6hkdoO8M3shvZj7F/MRodobCPnzp0zTz/99LZes7d6J2m1mf3Ga7RvLXDs3W+hMT3P/MXLhOUSKh8w9/wVokxaAFi55l6JLW/sNIB5ysoX7MRXrwBjW5El9O09aNtPMHByHBUG1F6fJG22SePESiyoTDPCGFQuoHzfMfKHKoy86X50FJM02oRDA0t2Si52vzu4tHDeGjy1aPBi3c4SqTDfniQxEZ4EeBJQ8CvdBiuAm43LNPUC2/kfa0FV0EZjSElJCKQARhOZVne3sHzXEEg+Cx/BodxJEN3tFZioXyTVmkg3bCJZhJOls0s0fZbTu1D27qL2gn6PiHzdGHOu32sHNmm7Ep3qnbQdM/0XL4CBcKTMrWdeYe78qww+ei86TVl4aYKk1kTyQaah02fcYAdtFo092HDPdht7AQIPUQrJeTZBGwaIwMB9RykeOUQ4XKY9PY+XC/DyIcFAIftcNl+hQp+Be8fJj5S7ncw3v/r8msPdHTvHSgle0BwvPUwlN8pY4RSHC6fRJDSSOdppI3uvTYRu93+sTb1A29SITJPUxLT0Am3TQJN04/jLiU2L1MRgFDPtCarxLS4tnAdgJDxBM5kj1i0CyVH0hpmJJlZN3G7X1K/txhn8ZXSqcuZfuIKXC/EKIX4+JI1iciMVFl6cwCQpAyfHCQaLdiBKltzc1RhsYtZTBKUCZz78Hv7+3O/zLT/5Q5x4z18jKBeovnqN5s05ooU6Xi6geGyUoFwgNzKAVwwp33+cwuFh0nbc9eaBdVU7OXaGnDfQrXgBa8Rn2lepxjO8Mv+XpFp3jVqobFNdM1no/taZLMJOs2jk7W5T4fU9JqZFbCIG/LFuRc5sdJ3B3BEOF04znD9KMaisabw3Wwm1W9ntZmrbKZ0cJ5qvE8/XUXlbxdAxcCNP3G+7cQW8fADaIEqhCgG3/aMQ7Lfb+7S3yX84y//f2ujpNJh2gkG6BUanPvgk1UvXad5asBVEStBRQtJoY6KYo+9+K7nRQQbuPYooIZqvkTTbFI+PEs3VGH37I13D36FftZNj+6lGM0Rpk5n2BLdaE8w2bzDdukIjmccYaKV1avEtpptv2Np6ExPrNvVkjqnm69SSGVJiPALUrkn39csNqGWv241pxyuvxlMbNt7LF0rYWCXUbmW3/L+4a+hU5ahcYKtxREhbEYfe/AAqDMiPDRFUSkQLDUQpyvcdpXHtFibR+KU8aZxg4hQV+ug4sdOkMFvjJPXuZLthpGX/8a9Duyett5j++os887FPUb8ySePaLYJSAT8IiKoNkmqDNE5oTM8TXpkiPzbE8KP3Ur9+i3iuTjg0QOXMCU598EmADclZOLaH3hj0UHCUhfgmtXQGD5+CV7HJzjRCo/El6QqZxSZCIbTSarfGXpOsWHy5U+gVHwkKxVzrOp4XAoYobdFIFpbkJtYy3huthNorOIO/jE71znO//Blu/Mkz5EYqjL39EVTWNTr+rifwcwG54TLX//Q8aTOiHS4QDpYoHB5i4ZVrxO0Y0zG8HYOsVGb0tygealY411qn9xQm1cw/e5n5MycYfOA4N7/6PCrwGP/OJxCBa//9r0jqLdKmrc4JKyUe/vHvXzEJu145C8f2saT+XEGQ5vF0Ew/baAW2hNKQYvBIdEw7rUO3PNKn4z3YuPluo1+lkGSfCRJidJriSYgxmpn2G9TjWQbDwyjldY33SpU4q1VC7eXqHVelswr9BncDS6p4Jr/8TVpT8xSOjBCU8rRnq7RnqiS1FgjoJLGzbkUwSXLXZtpK6GOUIElqSy2NsTvd5fmtrP7eL+Y59MQDHH3XE1z/0/NE8zXCQevxpC2r+ukVQo6+6wnas9XbJCccu5sXZp+yk6JEiNJmFrKxIYqSP4QBG9ohzeLhknn0u8Ee3InEbGf6kM4eSfe5vDcABiLdBDGEqsBQeISh8Agz0cSGKnH2QvXOalU6zsNfhZV0eToebXylweBjp+A523wUlIuMvfNRvDBg6qvPMfr2R0ijmOtfeAbdsnLBOrFTpJTngafQjf4aIbcReBCvUv+cpLbypmibuqK5Gsrz0MtnARir36zjlPkX32DokZMMPXKSyazHQAIfP5vwNfLmBwCXhN2LdOrPjdYsxDexJtA2WzWShcwgdpQrDavV1m8ngm9nKWz4fuwCoQjQWYevfeyhMkVNBHKqhC8hIoortQsU/CHywfq7cPd6565L2t4BI2fv59QHn0TCgLnzrxIOlDj2PW+ldGKM2fOvElaKnPtXH2HwwXvIVUqc/tB3ceZH3kt4qIJXyOEXcoBBt9qrx/Y9Bdmc2OWVMLehO5r1A3hhgJcLKZ08THioYs/jKQh9CHxU4KHyAeIpbn7lOQCGHz9NfmwIEQERxt7xKMUj9j9gl4Tde4zmT5KaiFo8A0bhqxBBCCiQtUdlmjQ7MbBhZQxLu3g3iidLtXYMtv4+MT2KmlhFzVjHLMQ3mW5dYa59g3baWDOZu9erd7ZqxOF7gV8APOBXjDE/t+z1dwG/B1zKnvq0MeZfbMW1d4JOc9b8i2+QGymDQPWVa1QePAEGbn71efLjI0s6TmcuvMrkV57ryiurwCea7cgsZ7rzvVILgYefCzDaRlDj+fqa/zZ1OyZpRegkpXhkBJOkHPuet5JGMZN/doGk1iAYLOGXi8QzVcKhMioXcOuZVxh8+B7e8W8/CtiQlZcLMFoz/8pV5p67TPnUUZ752KdcF+0uY60Y9AtzXwYRQskTSJ62rll1V8AntKMHWecuc9djxdp6h0sbUmJjw1Z5r9JV1IzSJtrEGCDPANqkVKMpCv7QkuTucjarY7TTbNrgi4gH/BLwPcAE8DUR+Ywx5rllh37ZGPM3N3u9raBfbH4jRqzTnKXbMUYJzeszRLNVZi9eJhweID9SIZqt8uzP/3a34/Typ59i+LFTzF28jAoDkmqDaDYbCG1sSMZ2tBrwFEoErQ1eIYdptjGe4OULpNXmyjemDcr38MKAY+9+C0e/+y1drZt7/863M/mn5+0c3qEBcm9+gPbUvNXMEVnSGdsJWd06/zK1y5MMP3aayv3HbvtMjp1lLV38cjjCaP7ebmdtNZoiVAVSnVrPlyYeAf2TPXubTgzfdCuNUhrpHB4+w7mjNJJ5fJUjNZEVbBOPRKfMR9cxpFxaON9dPHsXVcGzyW1/b1bvbIWH/zbgFWPMawAi8pvAB4DlBn9XcNvgkzswYt35t56i9uo1dBRbnRljuhOx0nbc7TgdOXu/FRi7/xhhpcjc81doxIkN1xg7fxYBk2hMmiJAOFIhqTVJGy28MMCvVCgdPcStb7yCafcXXfPLBR78X763m1xevqgVj4zcphDaScj2fvZO7uKZj32K8r1Hbhvs3vlMjp1lPfHkTnnhQvsWiWln4mOLiU1NmqlU7i+Db2dheSRkZdEA2aevRjMYDL4KKPnjVnQtbZDqNp4KKfnD3cVzJD7RTeyGqmQT3yJokxKZ+qo6RruRrTD4x4E3eh5PAG/vc9w7ReQbwDXgfzPGXNyCa2+YvoNPWN2ILd8RSBh0m43SKEZErMH3FKIElDD/whWOfPvZbrKzI7NcGB+hMD5iK2PmasTVJsoXjMYGxDAo38dozdAj95IfH7J1/mmKlw+p3HeM+ReuLJMoFIKRMqf+zrfzlo99eMVF7Z7veydv/P5XgPXVzK842N0lcHcF7bS2phJkORxhJD7BrdaEncomAb5XoJkuZInN3VhyeScIHsGSSiOdJaUF1Rmhjie+bSYTn2JP+GaufQMl3qKUQrZ4Xmu8SCkYXrqoYpu6Tlee2IHPuTm2ImnbL+24/L+hZ4B7jTFvAv4d8LsrnkzkIyLytIg8PTU1tQW3t5T1Dj7p0DGe0Wy1azyb16apvnadtNm2Q0CU/Y9J+T5BpYQoRTxfX5LsPPXBJ6308GwVozXNm3OElRLKV/gDRXKHKrZyB2Hw4ZMUjx7ixPveRvn0UeJaE6MNaTMiNzpI4fhoFv4Br5Rn4L4jjL/jUR79R98PLF3UevVtZp+9zOM/8YOEw2XqE1OEw+VVdzadruNeXAJ397DebtBGNn6w6A9S8MuI2H/2dgD40gTkXqBXx78XjaYSjGVJ285ORme9BnYHE6g8Oa/IQDCCpxSxbtNKGtSTWdq6QWLirpaQLyGRbuzpJO1ytsLDnwDu6Xl8AuvFdzHGLPT8/TkR+WURGTXGTC8/mTHmE8AnwNbhb8H9LWG9g086rLQj0FFMXG10DXd4z2Hat6yqoPgeKrd0vN9yOebC4SEKxw6x8NIEzRszNh+QasLhso3DD2TDwgdL+IUcg4+cpHntlp0le98xhh47hW5GlO87elseYjXPfL0jIGGx67g5NUfz2jStqXkk8HjiZ/7++r9wx5bQLzm73m7Qdlqj5I9Qi60DFactFD4p0Z4M5WhsqbHCy+6/U4IpRGmjW6Wjl1X7iCjaaZ1AFRnIjzCaP8nV+gvMR9dt+aYEKBTVaArCMZR4hKpIYqI9m6RdzlYY/K8BZ0TkNHAV+BDww70HiMgRYNIYY0Tkbdidxa0tuPaGWe/gkw4rGs+JKd7xbz/K13/6kyy8etUqSyqheWMGP/AZfuIBHv1H3983Ng6LO4fKgyfQcYIowRhDOFjsSjmAXYwqZ44z99xlTJySHxukcOwQXuCv6J1vdFFbiZGz93PP972T8z/7n9FRQn5skNKJMd74/a8w+OAJF8ffJlZMzpYe40RpbV38TmVJORzLBNESlCgwAZ74RGYXzlheBxpNKCVi08j8+YS2bqDwGc4dYa59M6tAsrF7MWLn75qkm5DNtYqM5E9gjO1XMBhEPGrxDKVgkGPFh5iJJvaNxMKmDb4xJhGRjwJ/iI1Cf9IYc1FE/mH2+seBHwB+XEQSoAl8yOxQi+9qg0/6sZrxHDl7P2/9lz/Kc7/8Gab/4nmU73P6B75jVRmCfveRNFpdjZr2bI2hx06RHxuiPVuleuk6GBh+7DT1iSmaN+dpz9V54mf+3orX2OiiBnYB6nwOgMPveISHf/z7mX32Mke/4023JXpd4nb7WC05e7ryxJoJw85OwJOQwXCc1CSkJgJCFB5JGqPZi0bfEBkbWhEUgeTRxkooN5MaOb9AaApEaZOUOGu8GiDvlbrfWScPIkqocJhGMk+s24hIt3u2GAyua9jMXsBJK6xBbwJ0OwZ8LE8QN27MdLV7OqxH6mAjpaczF17t7lTCchEDxNUGlQeOgxIOnb3fykNkGK2pT0zxbZ/8yS373I6V6ZVJAIjSJvV4jtg0OVy4b11aLr0hIVBEukmUNmkndYwYUpMsGR24FylkdfapThARiv4gqbHCb0o8hnJHusNgOgnXlQbE7NWkLDhphU2x0R3BVlyv99xf/tF/TXh4aMkx66mU2Uis/vKnn6I1NUdYLuEVbIJKRGhNzuEVQqL5+qbDQ447p7fZJ0qbNvRgIFCF22rvl7M89n+89Gi3tvxq/QUayTwePjkpkpg2Cb1J4N3VhbsW7bROqEpo2mAMraROaiI8FVDyRphpXSXWLUr+CC/M/g/IZOIibXtb9kPIZi2cwV8HGzGeW81WxeM79PP861cmSVsRQU/1kpcLaM/XKRwdsc1ZOPnjnaI3OVuP5zrN2RT9QeuZrqDlMlm/xJXaBTSaQPKkWjORLi4OD4d/HYD59hSxaSDikaNM29RhSf363kCT0tJVwKAI8L0ApT0SHVHjFqEqUPSGaaSzkAqVYBxPeWAMxug9WVe/UZyWzi5neTlne7ZKNFfrNldthH4lph0ZBS8fonsGsqftGD8XMvLEAxsq5XRsPeVwhBOlx6z+i2niq5ByOEbOKwL9ywSr0QxX6hcw2FJE22k6i9ZmyaSn46WHUQKhKlLwyoRejoAcexdbex9KjnIwyuHiKfJ+iVAVOJQ/QUrUVbpspQsEKkfeLxN6BR4efnJdOZG9jPPwdzlbGVJaqcQ0acfkx4ZsDD9rx+nE8Duxf2fgd5aOVAJwW8x5eZlgNbLjC1tpHY8AEYUvAamxYQ+vJx9TDke69fypia2S5C4bdnIniPJoJgvkvKKN42f5j87AdgzdHoa9XFe/UZzB3wNslcFdqcQ0nphaUm0EcPQ73rSuaiPH9rJW7X2nhDPKhp0YkxKldfBKeOIT6xY57+iSc5bDQzTiBSLdsMlcGjvx0bYIm3eI0yadkJQnHcll+511Ermdhqq9XFe/UZzBP0CsVWL65Mf/2Q7enWM9rDaJCRZLOENVIDLYOnRjm62MClEoRvMnl5yz6A1yo/Eqnvgkae/8hL2VtF1ESNEok9jdkJe3Ole6TcGrMB/fAISKP06s2/s6SbscZ/APEHdSn+/YWVaTP+59/Wr9OXLeANXoFgPBCEV/kMS08XWOlNhqwEvIydLZ22LUjXSecjBKM11YVqWz94y9Tw5NgiYh0YZ6PMtweJy2rjEfTQKGojdM4OUATaCK+zpJuxxn8A8Q211i6tgca8kf93u9ndZQ4lMKBruNRJGGgqrwwODb+hq2dlrL4tqaQEISk3TlC3Yeu8tQ+Nls3ZWlIASve4zCZ8AfIVAFrjdfohKMcSh3TzcEdrz08IEx8r04g3/AcAnYvcNa8sf9Xi/6WdkhEOkGsW6jRHGs+NCKBi7nDXCrNYESn9BTpGltA879ylr61gBvtpnLVt0o8cAoEmJWujmVqWIqfJQIpWDI9hmIT6QbFKWy50YSbjXO4N8FNjtgxeGAteWP+71e8Mu00wbNZI5YW+Oos/mtAOOl091jO+GgWjxDM10goEDo53uqdNaK4cuqR2xF565VtfdJu3o/kommLT13KHmMMSRE+BJQCcYJvQIL8RSeBEtURQ9SVc5yXB3+FrNSrfvMhVd3+tYce4y15I9Xel0Eiv4wgRcQegVynpXLuFK/QDWaARbDRbFuU/KHyakiCW2itIkAealkVSwrl2guLgwqM8NBNjN3q1AUvApK7PlVZvxV5qcKQsGrMJq7h4JfIfDy+BIyGB5lIBwGOlU58RKJ44NUlbMcZ/C3mJW06C9/+qmdvjXHHqMziDzWbUxWZZKaqFtls9Lrvgppp3XrHYvfLUHUZrHpqjccJCIMhuMEKsRXIQVvEKUUfteA9zf6Bno057UNu0jH319PLf9ax2iMSSn4FUr+IDmvRKBCDhdPUQnGssVACFSeodwRhnJHGM2f7OrcG2MIVIHUJISq2P2OWkmVKG3ywuxTXFo4310EDwLO4G8xGx2w4nCsRG+HbaTrBCrXVXBc7fWBYITYtKwBzrCLQcyNxsv81dQfcL3+ErVohii1OjKhV2AwOIJBE3g5RGAgHOVo6QFCitnUKMVSI700dh+bXn39tZMAah0R5aau0k6bCIpINwhUnkDlKPqDNpxkrNxEI15gvn2DdtqwuQOjiXSdUjDIfeW3UAwqRLqONimIIKKWJMIPitF3MfwtZqu1bxwHm94SzI28Pt18nURH+BIS6zaRbuARkJIQSxuDtnNciahwmNAroJTHaP4kpytPLCkH9TyPgEFaaZWUlJUStQJZU5PKGp1Wj/93xMtWH6BuaOt693wlf5gobdJMFkhNijYpMS1SIgr+EEW/0q3E6V0cO1xaON8dZQj95wDvZ5yHv8VspfaNw3EnlMMRTg6cRQRi0yYlxld5jKQEXoFA5fBVHk3c9ZCXh4vK4QinK09wvPQoSlQ2QqQTt7/dQNuySehU1azl4edVOds1rI1Bo03EUHiU1CQsxDdtmMYr4KsQMBQ8O59WRAhUDk9CrtZf4NLC+SWhm3Za21cjCzeKM/hbTKfW3YmNOXaS8dJpHh76Ng4XTuGJR04Vbivh9FUOX4XEpnlbuKjDdOsKBW8oM8yyYh28Tan63eNWw6dAziugxOtKHqyETdYqBA9tNDPRVdpJE200xmgEW8kT6aVyEKlOmWtfJ9btJaEbwVvXHOD9igvp3AVcrbtjN7BccK0ez5KaJNOW0YSqQCkYXnXYRzutUfQr+CogacdEOl0SrunMj01JWEzXmu7znV1BqPKUg0MsRFOkpKQmsfF07FByD5+YNou7B1uTY+vqA7SJaabziPEQgUhbKeNyOEYjtlOqeqknM914PyyGbrRJ0SbaNyMLN4rz8B2OfU6nmidURXSmL5OaJKtgiW7T1umlU/qZ84qM5I7jiY9HgE9uSWfr4iLQESxbVN1UKHwVUgqGGcmdIJAchtSGX8jjiYeRlJwqUFQj+BLiS2Bj7V4BIwmhV8SXEE8pcl6RglfBl4CcVyTnlVAiS6qVEmPLTXuxoRy9aiJ8v7MlHr6IvBf4BexM218xxvzcstcle/39QAP4n40xz2zFtR0Ox+r0Cq5ZjRkrEVwKBtccj9irzhmoPIEU7bCULNBiCLIwTyfKrwgkRERlHrzdEXTi5kp5jBfv43TlCV6c/Qpz0XUwQkpCYmIi08wknQVPAkKVB4yN1RuyMFCKMQpNajuJlXAyf5ZGOt/VHBoOjyGy1J/thG7WSoTvJCtpJ20Vmzb4IuIBvwR8DzABfE1EPmOMea7nsPcBZ7KftwP/PvvtcDi2gTs1csvVOQdzY915uKlOiE0T8MipAZrpgjXCJsbvmha7EBS8ym3KlAZDqlPbNWuMHTSePR9IHkRT9AchoStpHKoCRX+QWjyDJiFQub7iZ53Gsr0UullLO2kr2IqQztuAV4wxrxljIuA3gQ8sO+YDwK8ay1eBIRE5uvxEDodjd5Pziozn78OgiXULwSPnDeCrwIZc8PDFGntPPEJK1miL7hM+0QQqj8J66zbtm0Mg6yEQavEMBa9CpFs00yqtpE4tniHw8jw89G0rTqhaq4dhN7K8Ga5TbdQ7oWyzbEVI5zjwRs/jCW733vsdcxy4vvxkIvIR4CMAJ0+uHFt0OBx3n8n6Ja7UL6CNIVA5Ep3SSqtUgjHmzA0Ctai946sArVN8lWMsf++q9fBg8wPzTJP3yjTTBRuCMQbwQQwVf5xaMmWbyFB4kkOprLLHrN3YtZtDN/1YSztpK9gKD79fDdby/zfWc4x90phPGGPOGWPOjY2N9TvE4XBsA9Vohiu1Cxhjq1yM0TSTOVJtk7QKRaIjjDFokyIIA94hQpVfl1c9mj+JErH6P9iYv0HjZbsFq4GTo5XW8FXIYG6c0fxJDuVPkPfLW+r57gbW0k7aCrbC4E8A9/Q8PgFcu4NjHA7HLsImeTW+ChGxEsUqkxoGvaS5S0RR8IfIBXkeGHzbmgPBO8lJj5BYZ+MIjUF1qnuMohpPdzVw7HumaKe23n4/NkutpZ20FWxFSOdrwBkROQ1cBT4E/PCyYz4DfFREfhMb7pk3xtwWznE4HLuHdlojkDzapHjSUahUmV6NAFcpB2MIgiFdd1VJb3JyKHeEMCnSTOfwxdbMexKQmphKMJZp4DRITYISvzuYfD82S601vnIr2LTBN8YkIvJR4A+xZZmfNMZcFJF/mL3+ceBz2JLMV7BlmT+y2es6HI67S84bINWaRjpLmpVERrpJrNuUg1FCVVozTt+P5YNbSsEgoZdf0gD2wuxT3Xh20R9kIb65pBRzt1fc3Cl3O++wJXX4xpjPYY1673Mf7/nbAP94K67lcDi2h9H8SSbSixQZpp3WiXWL1ESU/VFKwSBwZ+Jj60lO5rwBO4BccoRegQqH1yzFdKyNk1ZwOBx96Q0x2A7Xo90h6b30Guv1NA71GvMOy0M0vQ1fftbIVQwqu760crfjDL7D4ViR5SGGSwvnVzTW620cWm7M+zVFbUc8+yDiDL7D4Vg3vcZa65R6MkusWwzljnKt/uKqQ9c7rNeY77U6+r2AM/gOh2PddIz1tfqLzMXX8SXHYHgUJR4z7QmGgqNLir1XKp90xnxncAbf4XBsiHI4QtgqWOVLtRja8SVHPZkl5xe7z623fPJui4Y5LE4e2eFwbJh+k6NK/gixbm24cagT+18+rOSgzJndTpzBdzgcG6afDICnPIZyRzcsWLYdomEOiwvpOByODbNSpc2dlE1uh2iYw+I8fIfDsWG2Un54O0TDHBbn4Tscjjtiqypt1lOX79ganIfvcDh2lL04rGSv4jx8h8Ox47i6/O3BefgOh8NxQHAG3+FwOA4IzuA7HA7HAcEZfIfD4TggOIPvcDgcBwRn8B0Oh+OAsKmyTBEZAf4LcAq4DPxdY8xsn+MuA1UgBRJjzLnNXNfhcDgcG2ezHv4/B/7EGHMG+JPs8Up8pzHmCWfsHQ6HY2fYrMH/APCp7O9PAX9rk+dzOBwOx11iswZ/3BhzHSD7fXiF4wzwRyLydRH5yGonFJGPiMjTIvL01NTUJm/P4XA4HB3WjOGLyBeAI31e+ukNXOdbjTHXROQw8Mci8oIx5kv9DjTGfAL4BMC5c+fMBq7hcDgcjlVY0+AbY9690msiMikiR40x10XkKHBzhXNcy37fFJHfAd4G9DX4DofD4bg7bDak8xngw9nfHwZ+b/kBIlISkXLnb+A9wLObvK7D4XA4NshmDf7PAd8jIi8D35M9RkSOicjnsmPGgadE5BvAXwKfNcb8wSav63A4HI4Nsqk6fGPMLeC7+zx/DXh/9vdrwJs2cx2Hw+FwbB7XaetwOBwHBGfwHQ6H44DgDL7D4XAcEJzBdzgcjgOCM/gOh8NxQHBDzB2OTaKrV2H6ArRnITcMo2dR5eM7fVsOx204D9/h2AS6ehUmvghxA8Ih+3vii/Z5h2OX4Qy+w7EZpi+AV4CgCCL2t1ewzzscuwxn8B2OzdCeBb+w9Dm/YJ93OHYZzuA7HJshNwxJc+lzSdM+73DsMlzS1uFg5cTrmgnZ0bM2hh/VIJqHuA5Kwcn3bOl9OBxbgfPwHQeeFROvk3+1ZkJWlY/DyGPQnIS4BkEJikdh5uKGE7cuAey42zgP3+HoTbzC4u9rT0E4CK1rNkzjF+zja3+ODsuLXnhUhcH77fuiKjQm7e9XPo1+4IPr99BXuo/pC7AJL9/tGhwdnIfvcKyUeG3egvo1SGPw8vb3wutw6+JSL3z2RdCJNfILl+xxfsmGeTbiod+FBLDbNTh6cQbf4Vgp8SoACrzAllx6AaRta9zr1+DWs/a3+FC/aj17CexxpBAObKxE824kgF3ZqKMHZ/AdjtGzkDat92uM/Z02IXcIMKBj+7yO7U+a/fby2XMptOashy9e9lwChcPWQ69dRV/6PPqFX7e/V/KuV7qP0bN3/tlc2aijBxfDdxx4VPk4+sR3Lo1zH32HfVy/CfHCYgxfPPAUqMC+WQLwQvDz4OdsGCccgGAImjezcE8ESQwmhvnXYfoC+uR7UONvBpbF2CUAk0I0172PTcXbc8N24ejkA8CVjR5gnMF37Em2OhGpysdvS4xqgNYXITxmjX3ShOa0DeGkMSjfevJoyI/A4bfClT+G5gyYySzu37aLQ+0KhBX73vY8vPxb6LmXYegMzFy0YZZwyF4jbcKJ79yaxGqnbBQWP0PatAua48CxqZCOiPygiFwUES0i51Y57r0i8qKIvCIi/3wz13Q4+iYiL30W/eJvrR022QCqfBxOfGdWfTNnfxfGIW3Z2H3tDUhaUDpmjfnMRfu6TiBuQXvO/q1TUJ6t0Y/nbSzdALUJu0CkqfXq51+1P42bcO3PN33/K36GrVpMHHuOzXr4zwIfBP7DSgeIiAf8EnbI+QTwNRH5jDHmuU1e23FQufqUDbWYxHqtXsl63nENSsdh5nmYfBo99CAcf3LLPH89+Vdw9StgtL2u1taItgfAG7deejxvPWgM4Nljk4Z9zUTW29eRjdG3ZkBCYAqaaRYeyqqBZl9EV6+u2vylq1ftwrBw2d5s+d6+n7ff7sVxMNnsEPPnAURktcPeBrySDTNHRH4T+ADgDL5jw+jqVZh7CfziYtK08Zotg0yaUL1sDadfhJrdCegNerQrhouuPQV+CInJqmmM9cyrV2yVjgZMOzuLWKOOtn+nTevZd94nHqQRkEDasAleL8sLKAEp2Vg/2N2MVwC8xcWseBySGsRVUHlbUTT/KkTz6NPfe9vndbX4DtieGP5x4I2exxPA21c6WEQ+AnwE4OTJk3f3zhx7j+kLtpvVYA2oBJkX3bRhE69gDacxEC3YncALv4Ye/ZbFapdVDJ+e/Cu48kfWew9KNtnasosGrVuQRPa6XYOO3Wmkmizqn2Gyn87f6eJD8e1zxmQLhLHnFgFRoEIYfMDeY6es0qRLF7PqZXtdLw+6asNGIrZ3YFmjVjcE1skTdGrxXWjnwLFmDF9EviAiz/b5+cA6r9HP/Td9nrMvGPMJY8w5Y8y5sbGxdV7CcWBoz9qwjcnKI00WOtFta0hV5sPEdRtr7xjWLM7Ppc+u2ISkq1dtTN1ItqgkVjIhTZfWrZvMuC75L3vF/6RvxyT2HCZZfF9St/cc1ex9tefsgtQpq+yt8VdB9vlTWwpqtN0xdBa56sTS67lafEfGmh6+Mebdm7zGBHBPz+MTwLVNntNxUOmUGVZOWyOYNDPjnIOgsFgj354BlK2ICSvWwNYmbBI1bUNxHMKyPWfHI56+YF8PBhZ3D2BF0XRkbbNuWwOLx+0e/WYx9p7TFsy9bMs8/aL9jEnTevOQefY50A27ExEve79koaLW0tO2Z+0C14urxT+QbEdI52vAGRE5DVwFPgT88DZc17Ef6ZQZegWrX5M0bTgEH5rXrYdrMsNnlE2YJg1o3Mg8cw/qN2yVTf4QlE6AzgxfdcLG1eNJ60n7A1ly9eZiuMUoIAWSLfgw/RaJ1P5SOWuQdWw/Y1yzXr/ybcindBLmXrDH6uw9JrHfy/JGK1eL78jYbFnm3xaRCeCdwGdF5A+z54+JyOcAjDEJ8FHgD4Hngd8yxlzc3G07Diq3lRmazNjlB2HszdaIKZXFxyNrIE1W9ojJ4u1tuyi0ZmDqPEw/i/7qv4T51xYvlCZ2IWnetF592rJhF+K7/Amz/IAK7W7Dy9uwjV+03rzW2UbAszuX3DB2p6EhNwLFMUiaS8tT70YHr2NPIsZsxVb07nDu3Dnz9NNP7/RtOHYx+tLnl3qv09+0XnxrFtBZHXy07F0e1rBqwFivWGVVM0ZDrmK95miBbpXNloRsNoD4mRZPlks4/MRSJc6wDMeehMmvWe+/k89Im1C5r2v4O01ctrb/KRueCgfh2JPdTl/H/kJEvm6M6dsX5TptHXub5fFpv2AraUghHLFJ19vIqmZQ4BVtGSTGSiQYbX+Cog0F6WxnsKbBzzzzbly/s3nW/Q9fC5PYhSyqWqN/4y/tLsPLQ27U7jZufj1L5vqLjV2ibC6jk5wF27dgYtsg1gmDzVxEFw+7Kp0DhhNPc+xKdHWdgmPLFSaL46BbgLIaOKafwc2e6xhKFfRo44hdNEa/Jav46RjytejU3Hs22duJ+W8GHdtQVLIAUZagbc/Dwst2MYgbNtavPNt0FQ5AULaLQAe/ANXXXZWOA3AevmOXoatXrUc695KtvikdX71uvN+IwSTCNjtFizIG/TARaAXhqD02qdrnxYfGVBY3TzOvP13vJ8iul4WLNkXP+00bCBYfR3OLgm5e3uYa/EL2OXoWwM7fTjHTgfPwHbuIboNQ7ao1toaswSi9zSPt7AC4+mc2fl17Y3HEoF+wYY3hh7LSy14PXbH4n73Y9zamsjJOscZeFDSuw5F3wsCxrNFqvXTCRXch5t/t3M3Q7cXqnWjBduum7SxR3ZOcrZxyg9YdgDP4jt1Et6s0sSEWL7DhkcbkEo/0NvG09qw10oP3w/CZrEpH2UVAejexkpU1eqAyj1dlw0o6lTy5QRh52M6lnfqrxcamdbGRheFO6FMKarKkdNy0n60wCgMnlgilcexbXZWOA3AhHcc2sqaeSycB6xds/FoCa8Q6jUcdj7S3czSq2vJKrWH2eRh5FBu/n7P2t3jEesAmM5YiNu6dtLJKmMpik5Z41mOefcmGhXQbioehvZDd4A5U66xJdk86tt/HCvr5/fT+XcL24OE8fMe2sK7Zqp0EbOFwVk7ZmS6Vwvwr3clRVCfsotCZIdsx1mlsH+uITGzHLgamx/OWMKvRT+0uoD2Txee1XRSShk2Mpg3b6eqFWfJ1PcZ+hxYDCaFwCHX6ff2NvRNOc2Q4g+/YHtaj59JpEFI+DJy09jOaB1IbYikdswtFcxKmvmHfG9UzyYFMcx4/K0/MRNV0DGEJvEqWwI2s7Q4Gs3JMsP8MeuLjJl5s0Grc7FlAdhvZImSsnES/aiY3xNzRiwvpOLaHdei5LBk1qGfh0COLc2I7NeUmtaJhurGomNmaAbT9O27YY1RoPfSwkunOxGCG4NAjqNPvQ7/4W7aOXXeSrB16DHvatj+7luxelW/VNRtTcPE/oskGrldO2d1MZ6GFnga1C04j/wDiDL5je1innsvyYR36hV+3YmYdGpP2sfYzkbQ61mB71tNPGoDY1yJtdwB+3lbtVB5YXGBMYhua0pm79pG3B2UTzSK23j5u2AXAy1sBtjSC4UdsRWenUzdugAjahXYOHC6k49ge7lTPZXljVdK0kYywYhO0CsCz4Rkd2Xi7BNghI1mHbNK0hr8jOQx2BGG0H+rQtd3xNCYX5+xi7KAWr2D/rl9bzHd0jhHfhXYOIM7gO7aFO56tunyhUL4Ns/gla+RMii2rzIy7nxm5zjSqzt9aWw+4eAR9+U9g9gV2Z1z+DohmbCK7U4nU6RpWvk3oJnUrDBfVoDWdJW+HXLftAcSFdBzbxp3MVl0S12/P2hrz+nXbGGVUZrOzISdgyy1NR9HSYzE+r0FyNswxfSGrz18uqraH0W3bNawCK44GttIpLNpdzexLWUlqaHMarSnbpOa6bQ8UzuA7dj23x/X/i/VYk4VsKHjM4gjB3hJMrGH38tnuQGU19s3FBWLfkAm3GZPNzo2splBhzHYtF8bsIR3vP43tHN6BE7bU1ZVsHghcSMexB0mtd1oYsz9LjLxnwxhgyzK9PIs19CYr2fRuP+WuolewzWNd/0xVDvLjtjS1E+YaOgOnv9f+XTq22NsQt6yBr12H2Rft3F9XsnkgcB6+Y++RG4b51633Hs9nWjcBVjrBs5Os4qZNysYNurF6Ubb+XswuD2V0p53b8EvaWON4ZTX8vRC8YdRf+8nuK7p61Q42j2rWu49bmUicZ8+tcravwc/bc4Ar2dzHOIPv2HuMns1i+vM9owe1NXgq895F2wXBpFhRtEw1M21xxxr120420GQ1JLAJ7LQFrTnwQvQ3fwUGjltZiZmLNqYf1+13lDRsOaofWI/fz9tkb/OmNfhORXNf40I6jj2HKh+Hk+8BUmu0VC5Lwho7hzZp2OQtOvPqi9b4+aGNa3cqevYEK92nWEOfH7EJ2yS2x/p5G6q5+Q146bdg4Q373OB9thfBJPanctqWturEfned0lenormv2ZSHLyI/CHwMeAR4mzGm7zxCEbkMVMmmP680fsvhWC9q/M3ouZehNpEZraxiJ25Y+QS/AImyofC0befTpp3kbub1b9jo91b97DQGUHboSafXwKhs3GF78ZhWDLMpDD8IQw9kPRC1TDYaW5ufZJ5+pzfi6Dt26DM57jabDek8C3wQ+A/rOPY7jTHTm7ye4wBzmwjY0JlMC75mvVbxbSXOwBk77ap+IzN+spjIBJbKGG/AiK82TGUnSJs2KV0oWqXQbu9BDya1TVfVK3DoMRveSVtZ1/OA1Shq3MhmCKystunYH2zK4BtjngeQDQ2IcDg2TlcEzCssVpRMfm2x81YbG6DUsU1GBuNQu0ZXNdP0aMmLb423jtlQPN/00aPfScS3eYxWNrzltqlc2b9Lo+1Q97gBngcnv8ca+fasHXZ+8rudkT8gbFfS1gB/JCIG+A/GmE+sdKCIfAT4CMDJkye36fYcu57pC7ZbtnVtcbRfXLeNRIceWTxu5nkrJXDoEevBxvUs8Wmw3rzJBpp0Jl/tlhDNHaB820XbEYtL+yxIKkvKop0H71jb4IvIF4AjfV76aWPM763zOt9qjLkmIoeBPxaRF4wxX+p3YLYYfALg3Llzu2kD7dhJalehOZ11iuZtPL41B7ny0uNKx+083Lhhj6WayS2oTGmzxaLR75RrBj2J3N28APTkHSS0Bj+qYheumNsXMM969EZg9E2o0+/b7ht27DLWNPjGmHdv9iLGmGvZ75si8jvA24C+Bt/h6EvSBJSV/YVs/KHYkMb0N61RL45bIzj8kPVmTTbEXOWtoe8Ye7LyTd05Zy47/24s1+w0YSm6Iw4lgKBs6/MlyATjsiQufqbfT/YZNZTG4fiTO3L3jt3FXQ/piEgJUMaYavb3e4B/cbev69hneFkViY6zMsJM915r2zmbRFYnpzAGp7/XavC0ZyFqQfW1zAhmIR3irNs2C+t0UlDiZZ7/dhl+tb5riW/LL6P5rBGrCUkN22gW2MR0fiSrtY9tnkMFdvBL+V44/qQL4ziAzZdl/m3g3wFjwGdF5Lwx5m+IyDHgV4wx7wfGgd/JErs+8OvGmD/Y5H07DhrlEzZEEy9YbzxtQVCx5YReaJ/zChAOLhq33DDUnrfGsDPqMMm6VnWCXQCSLPTT3OakbGfK1jqO83yIs2HqcdPuXADwsh2Mb3864a2H/ydn4B192WyVzu8Av9Pn+WvA+7O/XwPetJnrOByMnoXWFyE8ZsM3U+etR16+Z1ESwBgrvdz7nsmnwStmMgw1rEfdMZTZHNz2Ng5B6YRgwkGrWLnWsSbrGE7jnqqiTr9kCqm2C5ufs3o5QdEZe8eKOGkFx57gNpnksGw9/I6xh9u6RFX5OHr4Iduc5QXYhKbKhpJ3uJMGrDskKNtKGpMpWYqfhZB6+wM6SVnP7lyMzp6O7esqZ183WQOZKLtbiWqQv71pyg0wd/TiDL5jz9Ark7xkOLdfyMI8fbpEj33rYv1+9HQWyklBFbMu3KxO3wttPqB3aErHwOYq9vxxna4M8YYXicAqfMY3bRzeaBtK0i2r079kULrY3UrStCGrwhjIrI3Vq8BeX7LdSdrKdi5pNumLxe/n6lO2Yiko2eqljhrmegbPOPYlTkvHsSdZ7wStJceFA9bTz4/ahqOuxIKX6chnA0TEzzzsog37BAOQH7a5gPyw1aCRcIMyy4ktofRykBuxT3meDe3oziLTvWu6i4pXsqEaP2d3CF7eLha9C5OXg5HHbRx/4ovoyb+yi1ztqtXCN0D1sn2Pm3J1oHEevmPP0vH4u2GLq3+G7hO26B43ehYufdaGP9IspIKXzYHNDGx3VGIxC58k9li/AAP3QOmwPencZTtAJKlt7KaDki2TNKN2GEvaySssJwvvRLN2V3LyPbazuDmVDV9vLH6GwQcgP7j41mtP2UUibdpdgMm0hqpX7Bxgp4Z5YHEevmNPsyS0s8YQD1U+bgeCHHrEVv0cegwqJ60xFz/ztLX17AeO278LR6xXHVWtbnzxiP1pTWWe/np8JpXtBoxdNGpXbDewMfQ39mkWusk0600MxcP23ofOZDuDso3Jj70ZykcX3+oXbPmmTmxIqKOGabC6+M1pp4Z5gHEevmNvM33BhimCon3c+b3CEI/bxiVWr8K1P4eFy9Yb9nKQH7IGf/AMTD1tY/vBgA2/zFy0x3VKRG/Tr+lDR7pZxzaUVHvDGmaTaeD0ywcYbSuJ8sPdMIw6/T546O8u3vulz2cDXnpImvY+61etfHJSp7t7Ec9q6Jz87rXv2bEvcQbfsbdpz1rPvhe/ANWJdc1qVeXjS4xoL/rS5224xKTQmLQevtY2zJI/DIXDEF+h2wHbFy+z6cYmVZOmTRCPPJwNYl+h8cqkkBu0O4KVhpKMnrW7m85n7iSujz0Jl/6bDUuJbydcmcTmDnJDLmF7gHEhHcfeJje8OLyjQ2PKGsh1hHlWpZ3FzxcuZd5/3lbq6NiGWdIokzRYCbVo7MkStGnTdr8mTRY7f/vg52H4kcUKoT5hGFU+DiOP2fDQ1DP298hjqPE3w9CD2Elfxlb5jD5hh6CUT2zsO3DsK5yH79jb9PNymzeszvs6wzwrkhu26psSLGr4kFptnriWDV3pkQZXmaAZYvMByrP3lEbY3IBvO2VzoR3IorwVlBWyaqGFy/a6wUDfoSS6etWGmErHYPB++9lnLqKLh612TqccdbWyVceBwnn4jj1N3/LM/CEojC498E5mtY6ezWrvs+odHWddsuVFxc5ei23IauNjGDwFpSMwdJ9NDpdP2rGCA8cXyzklyGrnu2I+dCUXvMLaA1d68xci9ncn3r/OslXHwcJ5+I49z22J2E4ys+PZwx3NalXl4+ihB209e6c0s3yf3UEExayJSxa7XtH2sRfYmv/MwOpLn7ee//IdR9oGxm3Vjk6yoSzaev5D9y+OIYwb/XcnK+UvsoVt+fficDgP37H/GD1rwxdxI5vhms1qHT278XMNP2jj9Tq151Ke3T0MZLHwTiVNcTxLknq3e9Pt2UyTvwe/E2pZsAtR+Z5s0HoAuSGbJO49tt/upF/+wg0hd6yCM/iOfcdWhTO6MfLi0UwWoWYrdcb/GuqhvwtH355p75cAY0suDz0Kh9+y9ForGeaB4zb85Gd6/V5gH3fUP3uP7WfEt3JhcxwIXEjHsS/ZknBGJ0aeL1pjDtaoNm7YvzsJ49KxpYnR5QZ3pfLJTgK1N/zUXoD5V+0i0NHTWSHZepugXG7YjTB0rIoz+A7HSqwjRr4eg7vacRqWLgbKtyGj3JDdnaxhxF2c3rERnMF3OFYiN7xm8ne9Bnel45yX7thOnMF3ODJu044vHrExfLirtezOS3dsF87gOw4UKw0E6YqweYXF7tyZi7aTtXHDed+OfcFmZ9r+G+D7gAh4FfgRY8xcn+PeC/wCtpf8V4wxP7eZ6zocd0Jfo54NBFlRhK1xw4qW3Y17cZOoHNvMZssy/xh43BhzFngJ+N+XHyAiHvBLwPuAR4EfEpFHN3ldh2PjrNKZumKt/F3Qjt+IpLPDsZVsyuAbY/7IGNORCvwq0E+Z6W3AK8aY14wxEfCbwAc2c12H445YzahvZxPTaguPw3EX2crGqx8FPt/n+ePAGz2PJ7Ln+iIiHxGRp0Xk6ampqS28PceBZzWjvp1NTNu4m3A4elnT4IvIF0Tk2T4/H+g55qexouC/1u8UfZ5bURLKGPMJY8w5Y8y5sbGx9XwGh2N9rGLUt1VszEkiOHaINZO2xph3r/a6iHwY+JvAdxtj+hnyCeCenscngGsbuUmHYytYq+Z928oj1+q8dTjuEput0nkv8FPAdxhjGisc9jXgjIicBq4CHwJ+eDPXdTjulN1Q8+6arRw7xWbr8H8RyAF/LCIAXzXG/EMROYYtv3y/MSYRkY8Cf4gty/ykMebiJq/rcOxpdsPC4zh4bMrgG2MeWOH5a8D7ex5/DvjcZq7lcDgcjs3h5JEdDofjgOAMvsPhcBwQnMF3OByOA4Iz+A6Hw3FAkP6l87sDEZkCXr+Dt44C01t8O/sV912tD/c9rQ/3Pa2Pu/k93WuM6du1uqsN/p0iIk8bY87t9H3sBdx3tT7c97Q+3Pe0Pnbqe3IhHYfD4TggOIPvcDgcB4T9avA/sdM3sIdw39X6cN/T+nDf0/rYke9pX8bwHQ6Hw3E7+9XDdzgcDscynMF3OByOA8K+Nfgi8m9E5AURuSAivyMiQzt9T7sREflBEbkoIlpEXDndMkTkvSLyooi8IiL/fKfvZ7ciIp8UkZsi8uxO38tuRkTuEZEvisjz2b+7f7qd19+3Bp91DFh3APAs8EHgSzt9I7sNEfGAXwLeBzwK/JCIPLqzd7Vr+U/Ae3f6JvYACfATxphHgHcA/3g7/5vatwZ/nQPWDzzGmOeNMS/u9H3sUt4GvGKMec0YEwG/CXxgjfccSIwxXwJmdvo+djvGmOvGmGeyv6vA86wy43ur2bcGfxkrDVh3OFbjOPBGz+MJtvEfp2N/IyKngDcDf7Fd19zsxKsdRUS+ABzp89JPG2N+LztmtQHrB4L1fE+Ovkif51wds2PTiMgA8F+B/9UYs7Bd193TBn8LBqwfCNb6nhwrMgHc0/P4BHBth+7FsU8QkQBr7H/NGPPp7bz2vg3p9AxY//5VBqw7HKvxNeCMiJwWkRD4EPCZHb4nxx5G7PDv/xd43hjzf2/39fetwccOWC9jB6yfF5GP7/QN7UZE5G+LyATwTuCzIvKHO31Pu4Us6f9R4A+xybXfMsZc3Nm72p2IyG8AXwEeEpEJEfmxnb6nXcq3An8f+K7MLp0Xkfev9aatwkkrOBwOxwFhP3v4DofD4ejBGXyHw+E4IDiD73A4HAcEZ/AdDofjgOAMvsPhcBwQnMF3OByOA4Iz+A6Hw3FA+P8BsxrPFinQVNsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# write your code here\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"
   ]
  },
  {
   "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",
    "-->\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<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",
    "&lt;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",
    "\n",
    "1.  Loading data\n",
    "\n",
    "-   Overview data\n",
    "-   Data cleaning\n",
    "-   Data selection\n",
    "-   Clusteing\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### About the dataset\n",
    "\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",
    "\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",
    "</table>\n",
    "\n",
    "</body>\n",
    "</html>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-Download data\n",
    "\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)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-10-20 01:32:07--  https://cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud/IBMDeveloperSkillsNetwork-ML0101EN-Coursera/labs/Data_files/weather-stations20140101-20141231.csv\n",
      "Resolving cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud (cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud)... 67.228.254.196\n",
      "Connecting to cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud (cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud)|67.228.254.196|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 129821 (127K) [text/csv]\n",
      "Saving to: ‘weather-stations20140101-20141231.csv’\n",
      "\n",
      "weather-stations201 100%[===================>] 126.78K  --.-KB/s    in 0.08s   \n",
      "\n",
      "2020-10-20 01:32:07 (1.64 MB/s) - ‘weather-stations20140101-20141231.csv’ saved [129821/129821]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O weather-stations20140101-20141231.csv https://cf-courses-data.s3.us.cloud-object-storage.appdomain.cloud/IBMDeveloperSkillsNetwork-ML0101EN-Coursera/labs/Data_files/weather-stations20140101-20141231.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2- Load the dataset\n",
    "\n",
    "We will import the .csv then we creates the columns for year, month and day.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 12,
     "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",
    "\n",
    "Lets remove rows that dont have any value in the **Tm** field.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 13,
     "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",
    "\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. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting dedent\n",
      "  Downloading https://files.pythonhosted.org/packages/a8/0a/1cd1e8aa840995e0dc761ab4b801da54fb4ea41795d161b80597737aa607/dedent-0.5.tar.gz\n",
      "Building wheels for collected packages: dedent\n",
      "  Building wheel for dedent (setup.py) ... \u001b[?25ldone\n",
      "\u001b[?25h  Stored in directory: /home/jupyterlab/.cache/pip/wheels/2f/c7/0d/7a4a0c56ea2992c33d6316524582a9fdd51de51e3e9c080019\n",
      "Successfully built dedent\n",
      "Installing collected packages: dedent\n",
      "Successfully installed dedent-0.5\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install dedent\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: basemap in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (1.2.0)\n",
      "Requirement already satisfied: pyproj>=1.9.3 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from basemap) (1.9.6)\n",
      "Requirement already satisfied: six in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from basemap) (1.15.0)\n",
      "Requirement already satisfied: numpy>=1.2.1 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from basemap) (1.19.1)\n",
      "Requirement already satisfied: pyshp>=1.2.0 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from basemap) (2.1.2)\n",
      "Requirement already satisfied: matplotlib>=1.0.0 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from basemap) (3.3.2)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (2.8.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (1.2.0)\n",
      "Requirement already satisfied: certifi>=2020.06.20 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (2020.6.20)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (7.2.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib>=1.0.0->basemap) (2.4.7)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install basemap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: matplotlib in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (3.3.2)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (1.2.0)\n",
      "Requirement already satisfied: certifi>=2020.06.20 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (2020.6.20)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (2.8.1)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (7.2.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (2.4.7)\n",
      "Requirement already satisfied: numpy>=1.15 in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from matplotlib) (1.19.1)\n",
      "Requirement already satisfied: six in /home/jupyterlab/conda/envs/python/lib/python3.6/site-packages (from cycler>=0.10->matplotlib) (1.15.0)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'mpl_toolkits.cartopy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-39-e3a85e1fa6c1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcartopy\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mCartopy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mdedent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'mpl_toolkits.cartopy'"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.cartopy import Cartopy\n",
    "from matplotlib.cbook import dedent\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 = Cartopy(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. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"None of [Index(['xm', 'ym'], dtype='object')] are in the [columns]\"",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-6300e51d4d85>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpreprocessing\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mutils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_random_state\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpdf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'xm'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'ym'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnan_to_num\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2906\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mis_iterator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2907\u001b[0m                 \u001b[0mkey\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2908\u001b[0;31m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_listlike_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2909\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2910\u001b[0m         \u001b[0;31m# take() does not accept boolean indexers\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36m_get_listlike_indexer\u001b[0;34m(self, key, axis, raise_missing)\u001b[0m\n\u001b[1;32m   1252\u001b[0m             \u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_indexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reindex_non_unique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1253\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1254\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_validate_read_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mraise_missing\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1255\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1256\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36m_validate_read_indexer\u001b[0;34m(self, key, indexer, axis, raise_missing)\u001b[0m\n\u001b[1;32m   1296\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmissing\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1297\u001b[0m                 \u001b[0maxis_name\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_axis_name\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1298\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"None of [{key}] are in the [{axis_name}]\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1299\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1300\u001b[0m             \u001b[0;31m# We (temporarily) allow for some missing keys with .loc, except in\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{-1, 0, 1, 2}"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6- Visualization of clusters based on location\n",
    "\n",
    "Now, we can visualize the clusters using basemap:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "cannot import name 'dedent'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-19-814b2feba22e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'figure.figsize'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m14\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0m__version__\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0m_matplotlib_version\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mdedent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m \u001b[0;31m# check to make sure matplotlib is not too old.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[0m_matplotlib_version\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLooseVersion\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_matplotlib_version\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: cannot import name 'dedent'"
     ]
    }
   ],
   "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",
    "\n",
    "In this section we re-run DBSCAN, but this time on a 5-dimensional dataset:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"['ym', 'xm'] not in index\"",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-20-544e7ecd4e3f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpreprocessing\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mutils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_random_state\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpdf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'xm'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'ym'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'Tx'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'Tm'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'Tn'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnan_to_num\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mClus_dataSet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mStandardScaler\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mClus_dataSet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2906\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mis_iterator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2907\u001b[0m                 \u001b[0mkey\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2908\u001b[0;31m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_listlike_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2909\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2910\u001b[0m         \u001b[0;31m# take() does not accept boolean indexers\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36m_get_listlike_indexer\u001b[0;34m(self, key, axis, raise_missing)\u001b[0m\n\u001b[1;32m   1252\u001b[0m             \u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_indexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reindex_non_unique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1253\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1254\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_validate_read_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mraise_missing\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1255\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mkeyarr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1256\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36m_validate_read_indexer\u001b[0;34m(self, key, indexer, axis, raise_missing)\u001b[0m\n\u001b[1;32m   1302\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mraise_missing\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1303\u001b[0m                 \u001b[0mnot_found\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1304\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"{not_found} not in index\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1305\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1306\u001b[0m             \u001b[0;31m# we skip the warning on Categorical\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: \"['ym', 'xm'] 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": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "cannot import name 'dedent'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-814b2feba22e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mpylab\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mrcParams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'figure.figsize'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m14\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/conda/envs/python/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0m__version__\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0m_matplotlib_version\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mdedent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m \u001b[0;31m# check to make sure matplotlib is not too old.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[0m_matplotlib_version\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLooseVersion\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_matplotlib_version\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: cannot import name 'dedent'"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Date (YYYY-MM-DD) | Version | Changed By | Change Description    |\n",
    "| ----------------- | ------- | ---------- | --------------------- |\n",
    "| 2020-08-04        | 0       | Nayef      | Upload file to Gitlab |\n",
    "|                   |         |            |                       |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<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/).\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "language": "python",
   "name": "conda-env-python-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.11"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}