hw4_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSE 482: Big Data Analysis (Spring 2020) Homework 4 Part 2\n",
    "\n",
    "Due date: Sunday, April 5, 2020 (before midnight US Eastern Time)\n",
    "\n",
    "Submit your homework using the D2L system. Use the notebook below to write the solution of your homework. Make sure you submit the notebook along with its HTML version."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.** For this question, you need to write the code for a variation of k-means clustering algorithm known as bisecting k-means. The algorithm works as follows. Suppose you want to create 4 clusters. The algorithm will first apply standard k-means algorithm to bisect the entire data into 2 clusters, say, C1 and C2. It then computes the sum-of-squared error (SSE) of each cluster and choose the cluster with higher SSE to be further partitioned into 2 smaller clusters. For example, if cluster C1 has larger SSE than C2, the algorithm will apply standard k-means to all the data points in cluster C1 and divide the cluster into 2 smaller clusters, say, C1a and C1b. At this time, you have 3 clusters: C1a, C1b, and C2. Next, you will compare the SSE of the 3 clusters and choose the one with highest SSE. You will bisect the cluster with highest SSE into 2 smaller clusters, thereby creating the 4 clusters needed. You may review the example given in the attached slides, bisecting.pptx or bisecting.pdf.\n",
    "\n",
    "**(a)** Download the 2-dimensional dataset *2d_data.csv* from D2L and load it into a pandas DataFrame object named *data*. Assign the names of its two columns as 'x1' and 'x2', respectively. Display the first 5 rows of the dataframe and draw a 2-dimensional scatter plot of the data points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x20f50b68978>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df4xdxZXnv6fbD2iTDE2W1gYaGrM7yGgIi5u0wJGlETAzmF+B1iQ7wCYZbbSSlxVZxUzUK2eEAKNI8cqaJMOQgKzATlBYMDFMjwnOmkgwSoZdmLRpN44DXpEQsB/eYELahLiTtLvP/vHebd++79a9VfdnVd3zkVrufu92v3Ldut86dc6pU8TMEARBENynr+4GCIIgCMUggi4IguAJIuiCIAieIIIuCILgCSLogiAInrCirg8+44wzeNWqVXV9vCAIgpPs2bPnHWYeinuvNkFftWoVpqam6vp4QRAEJyGiN1TvictFEATBE0TQBUEQPEEEXRAEwRNE0AVBEDxBBF0QBMETastyEYQok9NtbN19AG/NzuGswQFMrF+N8dHhupslCM4ggi5YweR0G198ch/m5hcAAO3ZOXzxyX0AIKIuCJqIy0Wwgq27DyyJecDc/AK27j5QU4sEwT1E0AUreGt2zuh1QRB6EUEXrOCswQGj1wVB6EUEXbCCifWrMdDqX/baQKsfE+tX19QiQXAPCYoKVhAEPiXLRRCyI4Iu1E40XfGrN60RIReEDIigC7Ui6YpCnfi290F86EKtSLqiUBeBMdGenQPjhDExOd2uu2mZEUEXakXSFYW68NGYEEEXakXSFYW68NGYSBV0IjqFiP6FiGaIaD8RbY65hojoXiJ6jYheJqJLymmu4AqT022s2/Isztv0NNZteVa5jJV0RaEufDQmdCz03wG4kpkvBrAGwNVEtDZyzTUAzu9+bQBwf6GtFJzCxDc5PjqML//5RRgeHAABGB4cwJf//CKnA1OCG/hoTKRmuTAzA3i/+2Or+8WRy24E8HD32heIaJCIzmTmw4W2VnCCJN9knFCPjw6LgAuV4+PeB620RSLqB7AHwB8C+Dozvxi5ZBjAwdDPh7qvLRN0ItqAjgWPkZGRjE0WbMdH36TgJ74ZE1pBUWZeYOY1AM4GcCkRfSRyCcX9Wszf2cbMY8w8NjQ0ZN5awQl89E0KggsYZbkw8yyAfwJwdeStQwDOCf18NoC3crVMcBYffZOC4AKpLhciGgIwz8yzRDQA4E8B/PfIZTsBfI6IHgNwGYCj4j9vLj76JrPi205EwW50fOhnAvhW14/eB+BxZv4uEd0KAMz8AIBdAK4F8BqAYwA+W1J7rUIeVjW++SazIGUNhKrRyXJ5GcBozOsPhL5nALcV2zS7qeJhlQnDbUyzfVxExqhdSHGujJT9sPpu3TVBCMrI9rGp33wfoy4iW/8zkvaw6u6UVOFjnYkAH4sixVF0to9t/ebzGHUVEXRkE9+kh7WIB8/nXO6mCEHR2T629ZvPY9RVGi/oWcX3iguGepLvg4e1iAfP51zupghB0WUNbOs3l8Zo3hWzKzTeh775qf2x4rtx+15s3X0g1kc5Od3GE3vay3ZOEYBPfLST2XH79r2xnxV+8NJ8oRPrVy/zTwL+5HKfNTiAdowI2SgEeSky28e2fnNljNbl668j3tFoC31yuo1fHZtXvt+encPG7Xsxes8zy2b0OAucATz36hEA6ZaLzqrAhaJVWa0e2XiUDdv6zYUxCtTjqqor3tFoC133hv7q2PyyGT1t6ZtmuehmyNicy53H6pGNR9moq9+SLE2bx2hAUa4qE4u7rpTVRgu6yQ0N34y0pW/ag2ebLzQLeQesC0JgI1X3m6upiWHx7SPCAveUljJyVZn2g+pZbs/OYXK6XVrfNVrQVcKsQtcCB5IfPNt8oVmoelKyKf+6Sbi4OSoqvnFibuqqMu2HJG0pc0JstA89zieZRNgCz+M7tM0XmoUyMxyivvk7JvdZlX9tC1Vkbri4mowTXwDoJ8rs6zfthyRtKdN/32gLPc41csUFQ/juzGHMzi0PlppY4Fk+1zWLs6wMh7il7SMvvNlTi9l2K9EU0xVIVa4QF1eTKpFdZMbrW65T/l7SPTDth+D3NmpkvBUJccxypArGxsZ4amqqks/KslwP/87gyhaYgaNz806Kb1mY9Kvuteu2PKvtBiMg8QF1hag4A53JMcmKVPXT8OAAnt90Za1tqxtV3/QTYZE5dvyl/T+z9kMZ94mI9jDzWNx73lvoWS2ZwAJ3NShUBbqrFJM+NLFcbLYSTcjip67KFeLiajJu9Qic8KXHjb+0e5C1H6rO1fde0PMGdVwMCtmGSR+qlraE5UdguRZzSCKLOKv6qY8I5216urAVE+BeRlJUfOOyXKLjT+ceZOmHqidE74OieS0ZF4NCtmHSh6qSCp9aO2L9BpasZAkwq4JuC8ypgWPbinyVwfjoMJ7fdCVe33IdFhVu5fD4KyvIX3V2lvcWet6gjotBIdvQ7cOkkgpfGr+o3EbWSJZleRYrNKBpq06d8Zd0D7KKch3uWu8tdJUl85vfHdeySOJ+n9C5OT4X+SmCyek2Ru95JvZhihOstJIKvpI1DdbUCk16Lel119FJE1bdAwCZVzN1lBzw3kIPHorNT+1fVrdldm5eOzgKdG5OIEzBo9OencPEjpnUv9FEJqfbmNgxg/mFXqE5fWULd338Qu2AqG9Co7L48owhk5Vk01adSX7sNOt73ZZnM69m6hjP3lvoQOeGrjypd+5Kmi3DmzaCqounr2z1XDe/wNj81P7C2+w6W3cfiBVzAFh50grl7ro4fBKasvzXJpvVfNjYZkp4NfP8piuXZbAl3QtVCq2OKNcxnhsh6IDZbKm60arKjEkVG5tK0oBXuauaIDRJy/A8Oz9N3DauVEksmzSXyOR0uydAH8BA6j2qYzx773IJMFlmqm60r+SJxKt+N61OTlyAyMWcZ1OSijZFA2i3b9+LqTfe1Q4Im+wL8LmPdUkz8rbuPtCzQzlMWpCzjvHcGEE3ySQw9XENDvS6YlwhTyQ+6Xcn1q9W+tADdEoGB1arL+Kjmuj6iWIDwo+88CbGzv1QYf9n2Sh3gjQjT0cH0vzpVefwp7pciOgcInqOiF4hov1E9PmYay4noqNEtLf7dWc5zc2OyTJT5eM6fWULrb7eRdjs3LyzGS95IvFp6W9bP3lxbNwhTNJD42O+tGoZHlcREOiIepFZEbadS1onaS4RXV+3TUF7HQv9OIAvMPNLRPRBAHuI6PvM/JPIdT9k5uuLb2Jx6M6WKmv+un93Jp5++XCszzxq6biyrM0TiU/73XB/q2paJD00PuZLq5bh4SyqKEUKRlMyiXRIc4lMrF+tLK4VxqagfaqgM/NhAIe73/+aiF4BMAwgKui1U5SIxt3oVf9qILbqX5iwpePKsjZPCpvJ72bZPKMSOJMa9jaiMixu3743dnwVKRhNS1mMQ1cnxkeHe9Kdo7T6yaqgvVGWCxGtAjAK4MWYtz9GRDNE9D0iurCAthlR9PI8nOY0sX41/vdP300U84C3ZuecWtbmicSb/G6WzIp+is8xUL3uMuOjw/jU2pHYsgfR/syTDdOETKIkTHXiro9fmHxmQj3FapVoB0WJ6AMAngCwkZnfi7z9EoBzmfl9IroWwCSA82P+xgYAGwBgZGQkc6PjKHN5nhbtDnPW4IBTy9o8kXjT3zUNEKn8yqrXXedL4xdh7NwPJfZn3qBmEzKJkjDVibiNhWHmF9kqF6CWoBNRCx0xf4SZn4y+HxZ4Zt5FRN8gojOY+Z3IddsAbAM69dBztTxCmSKq+zcCS0d1821d1uaJxJcZxR9WuAeGLe3HIkjrzzRB0nEnuFY9sUiy6ETQX+dtejrWsLPJUEsVdCIiAA8CeIWZv6K45sMAfsHMTESXouPK+WWhLU2hTN9gUk716StbmD3We/BFlTWQ81JVADfpc+Leq7qWtAskCVKc9R6Xy+5KwL4MqooZ1YWOD30dgM8AuDKUlngtEd1KRLd2r/kkgB8T0QyAewHczBUfhVSmb1BVoOvTa0cwfedVy7YTA27txKsqNTDpc1TvAXCmH4skyUeuEo8+Imx+ar8ylz34Gz6mgppQVcyoLrw6gq5My8NXq6aqo8ySPgeIz1wpug0ukOUoNB2CvqzqfttMUTujTxtogQixK/QyacwRdGX6Bn31O1YVwM3yOTb5JqtC9yi0Lzw+YxQcDvrSpYB9WRQRM7J1x21jinMJ8ZhWhMuaMpf0OU2osqiL7lFoqvrnKga7O3alr4vB1tRkEfSGY+IXzON/TfocF3yTVaEruKYCzNy5f8d+f7znvab2dR5sXel45XIRzDHJS9bJ4U07vCHpc3yMUZiiyuy54oKhZUXKrrhgCE/safdcp/KtBwe6RN8fHGjh7hsuTL1/NlNHm23NeHFK0OsabC4OchN0fYppVkmaXzH8OUGf3r5975JACfETbFS827NzeGJPG5/46DCee/WIVk2YuGqOAHDqySuW3RMb/cJJVNnmyek27t65H7Nz8aUAdFY6ZWuJM4Je12Cr+nNtnjzSrBLdXXhxffrtF95cet8FISmT6ASrOgbt0RcP4m/+4uKePoqz8FWWe3iSdrEYWtY2mz5nk9NtTHxnBvOL8bGLYc2/UbaWOONDrysIUeXn2p4jnObr1vUrxvVpFBsCTLag6tcF5p7xodoDodpdG3YR2OoXTiLpwBAVWZ6zrbsPJIp5eB9K0t8oW0ucEfS6BluVn2tr5DwgbcOUbkBPt+9sFpIqSfLLxo2PcGG5QGh0As8uZsCo2kaAUqCzPGdFpNdWoSXOCHpdg63Kz3XBQooTiwDdbBXdvrNZSKokrl/D6IwPnd3LLmYbTaxfHXvuZ9LBIFmes6SxmHc8FznOnRH0ugZblZ/rooUURrfkgeohDJOUOpm1dKyrBP2qKhusOz6SJuPw57hUamF8dFhZCVUl0Fmes4n1q2NPKzOph16FljgTFDVJrysysFhluVEfilHpZMyMjw4nngSjCjC5mIVRFEUUffO1EqOqKqdKoLM8Z0GfhLNcTl/Zwl0fv1C7v6rQEmcEHdAbbGU89FUN8ibVqk4qjauqKeJiFkaR5BkfPk+GpgKt6kcAGL3nmaUTiqI5+i5Mdl4V5wKqKzYl5COtCFUcqnrUBOD1LdeV09Aa0F1hmqxEfX8u8q7KJ6fbmNgxg/mF5SOs1UfY+u97U0OzttF0zMfRmOJcgBuBRSGbtWnr7rwi0bWkk64DevvV9+cir/W8dfeBHjEHek8kyjNxVLHC9E7Qm/DQ+4LpQ+hDjCEN3Ydedd3dO/fjd8cXe4T+tIFW7A5HeS46JE1sgZ7kdVtJ2mIGXEy9cpkqs05czMIwRfehV103OzcfK/REkOcigaSJLchpz7tPRNIWM9CEh94W6tjZmpZ65zplVVucPTYvz0UCE+tXo9UfnxYa5LTntbAlbTEjLkSjfaDpWSdloOtWUl13SqtvKUsjzFmDA/JcJBD0iyqdNvCZ53HnStqiBdhcLKtufA+01YHuQ5+Ueud7nKEsxkeHldUqg/7N27dlT6oi6An4nLtbBFUFoJs2qeo+9EnXNam/iiRJtF3YJ+JdHnqR+J67m5ei8mpVfzuwlghYln9e1GcIQhy2GxCNykMvEnEpJFOWxRKdKKImh+9++iIExXZRshmXYw0i6AlITns6ZQx+nXrpvk6qRbj5fHEVyqRkTmraIhGdQ0TPEdErRLSfiD4fcw0R0b1E9BoRvUxEl5TT3GqRnPZ60BFrXyfVImri215XXwfbD3uxFR0L/TiALzDzS0T0QQB7iOj7zPyT0DXXADi/+3UZgPu7/zpNkktBrIfyUK2MAnyeVItw8/ngKpSU2GykCjozHwZwuPv9r4noFQDDAMKCfiOAh7kTYX2BiAaJ6Mzu7zpNnEvBlyWtrcRlGgSBUZ2zG11Gx82XZkz44Cr0YVKqA6OdokS0CsAogBcjbw0DOBj6+VD3tejvbyCiKSKaOnLkiFlLLcKHJa3NxO32/epNa/BzT3eHhklz8+m4InxwFbp+2EtdaAdFiegDAJ4AsJGZ34u+HfMrPfmQzLwNwDagk7Zo0E6rEOvBjCzuKZczDfKQljmk44pwIV86jSYUYisDLUEnohY6Yv4IMz8Zc8khAOeEfj4bwFv5m2cnPixpq0LcU+YkTWa6xoTrE6IPk1IdpAo6ERGABwG8wsxfUVy2E8DniOgxdIKhR33wn6sQ60EfCW4VS5OMCdcnpTrQsdDXAfgMgH1EFFSu+WsAIwDAzA8A2AXgWgCvATgG4LPFNzWdqjJPxHrQR9xTxSLGhJCETpbLPyPeRx6+hgHcVlSjslD10l6sBz2aZFFWgRgTQhLe7BSVpb2diEVZPGJMFItPe0q8EXRZ2tuJWJSCzfgWtPdG0GVpby9FWZQ+WVKCHfi2svfmCLqyN1NUeXam0IvU9hDKwLeVvTcWeplLe9+WZS7imyVVBbKiSce3lb03gg6UFywSMakf3yypshEjRA/fgvbeuFzKREdMxCVTLlLbwwypN6RHXN0gl0/D8spCL4u0ZZlYQ+XjmyVVNrKi0cenNFCx0DVIC7iKNVQ+vllSZSMrmmYiFroGqoAroD5IGhBrqAgksJcNWdE0ExF0TaLLsrgT76OINZQPcWVlRzZ0NRMR9IykHWQs1lB+JLsoHz75hgU9RNAzkuRO8f2YtKqQwJ4gmCGCnhFV5svw4ACe33RlDS3yD982fQhC2UiWS0Z8OLfRdqSPBcEMsdAzIkGn8pE+FgQzqHM2RfWMjY3x1NRULZ8tCILgKkS0h5nH4t4Tl4sgCIIniKALgiB4ggi6IAiCJ0hQVKgE2cIvuIwr41cEXSgd2cJfHa4Ij0u4NH5TXS5E9BARvU1EP1a8fzkRHSWivd2vO4tvpuAyUo2yGuSYvmyknWXg0vjV8aH/PYCrU675ITOv6X7dk79ZfiCHXnSQLfzV4JLw2ILOJKgap+3ZOeue61RBZ+YfAHi3grZ4hVhLJ5Da3NUgE6c5OpNg0ji17bkuKsvlY0Q0Q0TfI6ILVRcR0QYimiKiqSNHjhT00cnUZSWLtXSCPFv4ZZWjj0yc5uhMgnHjN4xNz3URgv4SgHOZ+WIAfwdgUnUhM29j5jFmHhsaGirgo5Op00oWa+kEWU8bklWOGVL7xhydSTA8flXY8lznznJh5vdC3+8iom8Q0RnM/E7ev52XOutp21Qp0IbMh7ja3Gntknro8aj6TWrfmKN7slPQv6oTymxZBeUWdCL6MIBfMDMT0aXoWP2/zN2yAqjTSrblCDBbU6502iWrnF7S+k0OtTDDdBK05blWoZO2+CiA/wNgNREdIqL/RES3EtGt3Us+CeDHRDQD4F4AN3NdFb8i1OlTtOVQY1t9+XmCUbZYQ3Vg6/10lWC1056dQx8R2rNz2Lr7gNKtZ8tzrSLVQmfmW1Levw/AfYW1qEDqnk1tsJZstXJ1g1E2W0N1YOv9dJHoameha4emrWJteK5VeF3LxfbZtApstXJNg1FNvX9RbL2fLpJ0LrCrqx7vt/4nzaZ1BAvvmNyHR188iAVm9BPhlsvOwZfGLyrt82y1ck2DUUIHW++ni6Stalxc9Xgv6CrqCBbeMbkP337hzaWfF5iXfi5L1G3NfLC1XbYj/VYcqky08Puu0dgTi1TpR2Ue8vxvv7hryU8Xpp8IP/3ytaV8puAHNqSe+kbUqAsz0Oq31r2XdGJRYy30OoJLcWKe9LpLiOCUh62pp64TXu20Z+fQT4QFZgw7PH4bK+hFbfwxEbJgwMS97jIiOOUiG6zKw7cYjddZLkkUsU3adGv6LZedY/S6K6gEZ+P2vVJ/pQCKXk1KfRx/aaygF5ESZ7rJ40vjF+HTa0eWLPJ+Inx67UipWS5VkCQsUn8lP6pVYx+Rcb9KfRy/aWxQtAjO2/Q04nqPALy+5bqqm1MbqgBzGJ1g8+R0G5uf2o9fHZsHAAwOtHD3DRd6tSTOQpHBuzqSAYRiSQqKNtZCLwLZ5NEhrbwokO4emJxuY2LHzJKYA8Ds3DwmvjPTeOsxWE3GxVpMN8DITlO/EUHPQVY/vG8+TJ3yommT3NbdBzC/0LvemV9kJ3fsFc346DAWFatpEzEWI8RvRNBzkMUP76sPc3x0GM9vuhJfu2lNpkkuSZTEeuxQhBhLzXS/aWzaYlGYpj35noKWdSdj0q49sR47FLHtv8qdprI3oXpE0CvGVR+mycOZJbd3Yv1qTOyY6XG7tPpIrMcuRYlxFbnXsjehHkTQK8amk4x0qeLhDP5Ok7NcdCZNVzbC+LoStX3VIWmLMZR50+JS0GyuGwFIqlsVqFITXZ3UXE/pjdMAAFY8u1LLxYA7JvfhkRfeXBqMJtaoroUFuFUtz1U3kUuoanPPzs33jD/brUTAzZVogGpFevKKPutXHSLoISan28vEPCDYxr519wHlw2PilnBl2Ryg+3C6IDS2kjQ5hkUjbpxN7JjB3Tv34+jcvDX97nLddpW7SHUYhk2GjaQthti6+0DsMjEgKcXQ57MedVLdfE3HrIo0yzUQjbhxNr/AmJ2bt6rfXT5tylSgbVp1iIUeQudGqpZYPrsldNxEpkEwseaXE2fRhukjUvqlo9jiBnBtJRqQdvBFGNtWHSLoIXRvZJxIu+wz1CHt4TSZ0HTcU00T/LgsnzCmNfN9MCTqIm1yDWPbqqMRLhfdrfY6NUmAeJFu+g68tF2M4XvwV4/vTXRPNdV9Mz46jOk7r8LXblqz5KrIWivfF0OiDqLuItU9GB4csErMgQYIuok4RG/k4EALrf7lN1Ml0i77DIsgaUKL3oNFhbGZ5Cf2JR6hQ1BG4fUt1ynrtwCdFMDTV7bQ6tMbo4I+46PDmFi/GmcNDmCBGVFJt7WPU10uRPQQgOsBvM3MH4l5nwD8LYBrARwD8B+Z+aWiG5oVU99u2LWgU861aa4BFUl+9nVbntVavgZWpc/xCFNUrrzwHgDdMZh2nYzlE0TdgozOBMqA1UfU6fjQ/x7AfQAeVrx/DYDzu1+XAbi/+68VZBWHuI0evzu+mHiNr9ubdR90lZ9dV4gDi8f3eIQJOul/OsHHtLHalLGsS5whGIi5zZvpUl0uzPwDAO8mXHIjgIe5wwsABonozKIamIfJ6Tb6FP4vnXKuacv+JrgGivBn6wjx4EBrSTiaFo9IivHEufI+8dFhbN19wKj8ctpYbcJYNsHVVWIRWS7DAA6Gfj7Ufe1w9EIi2gBgAwCMjIwU8NFqAiGKyw7IU841bDm6etNNKKImR1rWwECrH3ffcOHSzy7ups2KjmUcdQNmyRBKG6tNGMsm1HGIfBEUIehxJnBsJIeZtwHYBnRquRTw2UpUW6n7ibSClaobSuiUB3ju1SPKnGCfXANFPOjjo8OYeuNdPPriQSwwo4+Ak1f04bfzizhtoAUi4PbITlxXc5hNMZ0w065XCf4prT7MzS/2/L1grPrk5sojosHvtmfnlnzmAVkPka/SjVWEoB8CED62/mwAbxXwd3OhEpxFZq3OnFi/Grdv39sj2gzElgcIqMM1UKYVUMSDPjndxhN72kurpU6WC+FTa0fwxJ52o/22phNm2usqwY8jXJrY5a36YbKIaFjEw2QJhIafxT6iHg9B2Zu+ikhb3AngL6nDWgBHmbnH3VI1eU93GR8dVoq26vU6UhXLztkuwp+tEplHXzzYeL+t6Tg9baCVeL3JyukDp6xY5tbxIe3WNBYQfn7iiAZCk/azRJ9F1WawMt1YOmmLjwK4HMAZRHQIwF0AWgDAzA8A2IVOyuJr6KQtfrasxppQhMUxbLAFmIBaot9l150uwp+tGsB1DHjbMBmnk9Nt/Ob3x3teD1vaJtvWZyM7Un1wc5mubFSu2ejv6lj+On8LKNeNlSrozHxLyvsM4LbCWlQQRQhR3MMW9asF1OVrrCKYpXrQdV09KpHpj1mSBtc3BZNxqjpIO2xpuzBmy8TURajznJw1OKBlOOn8rbLdWF7XcslrccQ9bFdcMLTM7wvU62usK5hl4qtUWaGf+OiwVX1ZJGUc2acSjLCl7cKYLRPTFU+cnzsM4UQ8LY727BzO2/Q0zhocwGkDLczO9dbh6SfCIrMzWS5eE/ewjZ37IWtS6uoKZpm4epKsUJv6sijKym7QnbxtH7NlEJ5ATxto4ZRWH2aPqevDJ6U1h2F0xu3gylZs0bTgmvbsHFr9hFYfYT5U26LqE43kCDoPqGPLdpVHjLm2JX30nmdiH/68uwxdPL6wCrL0i+pYRRWtPgIIsS6vMKevbGHlSStKHatyBJ3nVBHMioqqanlZtKvHtS3pk9NtpSWXN67RpA1XJmRJDDC9F/OLjMGBFk49uSPWKlmfPTaP6TuvMvrbRSKCLqQSJ6qq5WXRrh6dh9UmCz4p5bKIyc6HTJSiyZIYYJINFHB0bh577+qItcrCrzvQLIIupKI69qyK5WXaw2qbBZ8kIrYFIW2aCPOQJTHA5BCLACJgzeZncHRuHqd1S2uHXTA2BJpF0DMS9zAAfi6Hk7Iryl5epj2sZefhm6Jqb7j4mA3YNhHmIUtiQNR9FZSgmD02j8GVLbz/2+PLVp9AZ4dz4GacnZtHq49w+spWYvC1akTQMxB78vp3ZpYFTdqzc9i4fS/u3rm/p4a6a5SdGplkKaY9rLYVlVK1N1x8LI6qrWXbJsI8ZI0tJLmvJqfb+MLjM4lZMPOLjJUnrajVZx5FBD0DsS4IxTE8s3Pzzlo+AWWmRiZZisCJvg42IUVrathWVCqLuMT1we3b92Lj9r2lHaZg20SYlzRxziL2qtzzMLb1lwh6BkxvoquWT0CZ2RUqS3HzU/vx2/nFpfcWmJcmkfDn2lhUyjRwqTpMASjPFaLKUlLVinGVPK4lncBp3UHQKCLoGcgSIbdtJjelrOwKVb/Epf7FTYw+pPKljY0yDALV2dMZz6S2ljyuJZ06/ldcMIR1W561ZuyJoGcgS4TctpncFkwnxzjxcz2VT6cPijYIooW50l53lTyupaixMLiyBeZO+mJcSQUbAstOCbotaVbRG51WD6JuF4DNqFwmJ6/oM880J/wAAAnzSURBVNq4ZMvYyIKOgVC0QWBb7KEsdP6fSWMnyViIO/y8bveqM4JuW5pV+EZPTrcxsWMmdlvw4EBLmeXisggVhcplAsCoyJJNY0OF6n6H+6CIk3J0sDH2UAZp/888Y8fGwLIztVxUO7NsOYV7crqNzU/tX/L9Jgl5cH2cVZb2e7ZQxWSk+xm2jw3ArN5IVRN9UwyKpP9nnrGj+t1+IvzNX1xcWl96UcvFxtkwTNTSPPXk5K5VFcN3Ic3xjsl9y47hK8sizltW1paxAZhXp6zi3rsee9Al6f+ZZ+yoXGULzLU9w0UcQVcJeY+UKxvTo+CSBozNx7BNTrdjz1Sts822jw1Afb+DetpxR5oJ2ZmcbiceFxeQZ+yMj3aO7euPSQ2q63lwRtCLONuyTEzPMkwbMDZZl2G27j6grDRX5+5Mm8cGkHy/854FqyteTcHEuMo7dsZHh7Fo0VGKzgh6MBvaeoit6dItbiCFscm6DJNWwa4ObB8bQPr9BrJZdWUfEu4iJsZVEWPHphWiMz50wG6fn2kaWPD/CAdSA2yzLsOo/p/BUV11YfPYAHpjLEWtcnyqyVIUpsZV3rFjU8aQMxa67WRZuo2PDmP6zqvwtZvWWG1dhon7fxKAT60dsbbNdRO4RILaIF/t3u84TK06FwLCVVO1xWzTCtEpC91m8mxBt926DOPDVvsqUeU5F3VAdlM2CJmQ12LOWszLhmdAKw+diK4G8LcA+gF8k5m3RN6/HMA/Ani9+9KTzHxP0t+UM0WFJpCU5zyxfnXuiVHOGY0na469C/2ZKw+diPoBfB3AnwE4BOBHRLSTmX8SufSHzHx97tYKgkckuUSKsOpkxRRP1r51PSah43K5FMBrzPwzACCixwDcCCAq6IIgRKjCJWLLct8HXI9J6ARFhwEcDP18qPtalI8R0QwRfY+IYo9nIaINRDRFRFNHjhzJ0FwhQHKP3cCFHHnhBDalIGZBR9DjKiRHHe8vATiXmS8G8HcAJuP+EDNvY+YxZh4bGhoya6mwhOQeZ6fqidCmDAghHdcnYB2XyyEA54R+PhvAW+ELmPm90Pe7iOgbRHQGM79TTDOFMK77+eqirqqM4hJxB9djEjqC/iMA5xPReQDaAG4G8B/CFxDRhwH8gpmZiC5Fx/L/ZdGNFTq47uerC5kIBR1cnoBTBZ2ZjxPR5wDsRidt8SFm3k9Et3bffwDAJwH8FyI6DmAOwM1cV13eBmBr7rHt5VhlIhR8R2tjETPvArAr8toDoe/vA3BfsU0TVNi01TggzZ1hg9jbOhEKQlHI1n8HsTHQluTOsCWI63rASxDScObEIsFuztv0dGzBKYLaMq7jRCEbVgqCkAcvTiwS7CbJnVGm79pUoF0OeAlCGuJyEQohyZ1R1mYNW1w5gmALIuhCIST59cvyXZueEiUIviMuF6EwVO6MsjZrSBqiICxHBF2ohDJ8101NQ5TArqBCXC6CszQxDVHiBkISIuiCs9iYj182EjcQkhCXi+A0TUtDlLiBkIRY6ILgEK7X6xbKRQRdEByiiXEDQR9xuQiCQ7her1soFxF0QXCMpsUNBH3E5SIIguAJIuiCIAieIIIuCILgCSLogiAIniBB0QxILQ1BEGxEBN2QtLMzBUEQ6kIE3ZCkWhoi6IJgD01cSYugGyK1NATBfpq6ktYKihLR1UR0gIheI6JNMe8TEd3bff9lIrqk+KbagdTSEAT7aWpVylRBJ6J+AF8HcA2APwJwCxH9UeSyawCc3/3aAOD+gttpDVJLQxDsp6kraR0L/VIArzHzz5j59wAeA3Bj5JobATzMHV4AMEhEZxbcVitoYg1uQXCNpq6kdXzowwAOhn4+BOAyjWuGARzO1TpLkVoagmA3E+tXL/OhA81YSesIOsW8xhmuARFtQMclg5GREY2PFgRBMKepVSl1BP0QgHNCP58N4K0M14CZtwHYBgBjY2M9gi8IglAUTVxJ6/jQfwTgfCI6j4hOAnAzgJ2Ra3YC+MtutstaAEeZ2Ut3iyAIgq2kWujMfJyIPgdgN4B+AA8x834iurX7/gMAdgG4FsBrAI4B+Gx5TRYEQRDi0NpYxMy70BHt8GsPhL5nALcV2zRBEATBBKm2KAiC4Aki6IIgCJ5AHW9JDR9MdATAG90fzwDwTi0NMUPaWSzSzuJwoY2AtLMIzmXmobg3ahP0ZY0gmmLmsbrbkYa0s1ikncXhQhsBaWfZiMtFEATBE0TQBUEQPMEWQd9WdwM0kXYWi7SzOFxoIyDtLBUrfOiCIAhCfmyx0AVBEISciKALgiB4Qu2Cnna8nQ0Q0UNE9DYR/bjutiRBROcQ0XNE9AoR7Seiz9fdpihEdAoR/QsRzXTbuLnuNiVBRP1ENE1E3627LSqI6OdEtI+I9hLRVN3tUUFEg0S0g4he7Y7Rj9XdpjBEtLrbh8HXe0S0se52mVCrD717vN3/BfBn6JTg/RGAW5j5J7U1KgYi+mMA76NzKtNH6m6Piu4pUWcy80tE9EEAewCM29SfREQATmXm94moBeCfAXy+e9KVdRDRXwEYA/AHzHx93e2Jg4h+DmCMmW3dCAMAIKJvAfghM3+zW7l1JTPP1t2uOLra1AZwGTO/kXa9LdRtoescb1c7zPwDAO/W3Y40mPkwM7/U/f7XAF5B5+Qoa+geU/h+98dW98vKyDwRnQ3gOgDfrLstrkNEfwDgjwE8CADM/HtbxbzLnwD4qUtiDtQv6Kqj64ScENEqAKMAXqy3Jb103Rh7AbwN4PvMbF0bu3wNwH8DsFh3Q1JgAM8Q0Z7uqWA28m8AHAHwP7ourG8S0al1NyqBmwE8WncjTKlb0LWOrhPMIKIPAHgCwEZmfq/u9kRh5gVmXoPOyVaXEpF1biwiuh7A28y8p+62aLCOmS8BcA2A27ouQttYAeASAPcz8yiA3wCwNWZ2EoAbAHyn7raYUregax1dJ+jT9Us/AeARZn6y7vYk0V1y/xOAq2tuShzrANzQ9U8/BuBKIvp2vU2Kh5nf6v77NoB/QMeVaRuHABwKrcZ2oCPwNnINgJeY+Rd1N8SUugVd53g7QZNuwPFBAK8w81fqbk8cRDRERIPd7wcA/CmAV+ttVS/M/EVmPpuZV6EzLp9l5k/X3KweiOjUbgAcXRfGVQCsy8Zi5v8H4CARre6+9CcArAnWR7gFDrpbAM0Ti8pCdbxdnW2Kg4geBXA5gDOI6BCAu5j5wXpbFcs6AJ8BsK/rowaAv+6eOGULZwL4VjeLoA/A48xsbUqgA/xrAP/QmcuxAsD/ZOb/VW+TlPxXAI90jbefwcKjKoloJTpZd/+57rZkQbb+C4IgeELdLhdBEAShIETQBUEQPEEEXRAEwRNE0AVBEDxBBF0QBMETRNAFQRA8QQRdEATBE/4/QDwBhX/c1AgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "name=['x1','x2']\n",
    "#matplotlib inline\n",
    "data = pd.read_csv('2d_data2.csv',encoding='utf-8',header=None,names=[name[i] for i in range(2)])\n",
    "data.head(5)\n",
    "plt.scatter(data.iloc[:,0],data.iloc[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)** Write the function bisect() that takes three input arguments: data to be clustered, k (number of clusters), and seed (for random number generator). The function should return the cluster labels as output. To ensure your code is correct, you should maintain a list of sum-square-error (SSE) of the clusters and display the SSE list at each iteration. As the cluster is bisected into 2 smaller clusters, the number of elements in the SSE list should increase by 1. The process continues until there are k elements in the SSE list.     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans \n",
    "import numpy as np\n",
    "import math\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "#Calculate list to add elements multiple times\n",
    "def fun1(list1,v,n):\n",
    "    i=0\n",
    "    while i<n:\n",
    "        list1.append(v)\n",
    "        i=i+1\n",
    "    return list1\n",
    "#Dictionary sorting to get tags\n",
    "def fun2(list2):\n",
    "    values = []\n",
    "    y = []\n",
    "    for i in range(len(list2)):\n",
    "        values=fun1(values,i,len(list2[i]))\n",
    "    keys=[n for a in list2 for n in a ]\n",
    "    dic=dict(zip(keys,values))\n",
    "    for i in sorted (dic) :\n",
    "            y.append(dic[i])\n",
    "    return y\n",
    "\n",
    "def bisect(data, k, seed = 1):\n",
    "    init_coor=[np.mean(data.iloc[:,0]),np.mean(data.iloc[:,1])]\n",
    "    SSE=[]\n",
    "    sum=0\n",
    "    for i in range(data.shape[0]):\n",
    "        ds=(data.iloc[i,0]-init_coor[0])**2+(data.iloc[i,1]-init_coor[1])**2\n",
    "        sum=sum+ds\n",
    "    SSE.append(sum)\n",
    "    print('Iteration 0  SSE =', SSE)\n",
    "    datall=data\n",
    "    datall_list=[]\n",
    "    datall_list.append([i for i in range(0,240)])\n",
    "    max_S_index=0\n",
    "    for numClusters in range(1,k):\n",
    "        SSE.pop(max_S_index)\n",
    "        datall_list.pop(max_S_index)\n",
    "        clusters = KMeans(2, max_iter=1,random_state=seed)\n",
    "        clusters.fit_predict(np.array(datall))\n",
    "        lable_pred=clusters.labels_\n",
    "        centroids=clusters.cluster_centers_\n",
    "        SSE1=0\n",
    "        SSE2=0\n",
    "        dataA=pd.DataFrame(columns=['x1','x2'])\n",
    "        dataB=pd.DataFrame(columns=['x1','x2'])\n",
    "        for j in range(len(lable_pred)):\n",
    "            if lable_pred[j]==False:\n",
    "                d=(datall.iloc[j,0]-centroids[0][0])**2+(datall.iloc[j,1]-centroids[0][1])**2\n",
    "                SSE1=SSE1+d\n",
    "                dataA=dataA.append(datall.iloc[j])\n",
    "            elif lable_pred[j]==True:\n",
    "                e=(datall.iloc[j,0]-centroids[1][0])**2+(datall.iloc[j,1]-centroids[1][1])**2\n",
    "                SSE2=SSE2+e\n",
    "                dataB=dataB.append(datall.iloc[j])\n",
    "        dataA_list= dataA.index.tolist() \n",
    "        dataB_list= dataB.index.tolist() \n",
    "        #print(dataB_list)\n",
    "        #print(dataA_list)\n",
    "        SSE.append(SSE1)\n",
    "        SSE.append(SSE2)\n",
    "        datall_list.append(dataA_list)\n",
    "        datall_list.append(dataB_list)\n",
    "        max_S_index=SSE.index(max(SSE))\n",
    "        data_tmp=pd.DataFrame(columns=['x1','x2'])\n",
    "        data_tmp=data_tmp.append(data.iloc[[i for i in datall_list[max_S_index]]])\n",
    "        datall=data_tmp\n",
    "        print('Iteration', numClusters, ' SSE =', SSE)\n",
    "    labels = fun2(datall_list)\n",
    "    data['labels'] = labels\n",
    "    return labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** Apply the bisecting k-means algorithm to generate 8 clusters from the data. Assign the clustering result as another column, named 'labels' of the dataframe. Draw a scatter plot of the data using the cluster labels as color of the data points in the scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0  SSE = [1365.6945979485834]\n",
      "Iteration 1  SSE = [206.10668865383343, 215.70502668666666]\n",
      "Iteration 2  SSE = [206.10668865383343, 43.43106116466667, 44.50749403133333]\n",
      "Iteration 3  SSE = [43.43106116466667, 44.50749403133333, 47.96253796372882, 47.72318725278689]\n",
      "Iteration 4  SSE = [43.43106116466667, 44.50749403133333, 47.72318725278689, 7.436857757666665, 6.254233716551726]\n",
      "Iteration 5  SSE = [43.43106116466667, 44.50749403133333, 7.436857757666665, 6.254233716551726, 7.481366151612905, 5.297583828333333]\n",
      "Iteration 6  SSE = [43.43106116466667, 7.436857757666665, 6.254233716551726, 7.481366151612905, 5.297583828333333, 4.395181067, 4.958932097666666]\n",
      "Iteration 7  SSE = [7.436857757666665, 6.254233716551726, 7.481366151612905, 5.297583828333333, 4.395181067, 4.958932097666666, 5.751973757241379, 5.042813212903225]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAD8CAYAAAA2Y2wxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3xUZfaHn/dOT68QeuggVUCqIIKIAvaGil0RdXetP+u6rrrFta19xd5QUAErqChFEAHpvfcECAnpyczc8v7+mBCSzEwyIZPMJMzz+VzN3PvOve8NkzPnnvec7xFSSiJEiBAhQv2hhHoCESJEiNDUiRjaCBEiRKhnIoY2QoQIEeqZiKGNECFChHomYmgjRIgQoZ6JGNoIESJEqGcihjZChAgRaoEQoqsQYm2FrUAIcU+174nk0UaIECHCySGEMAEZwCAp5T5/4yIebYQIESKcPKOBXdUZWQBzA03Gi5SUFJmenh6qy0eIEKERsWrVqmwpZWpdzjH27GiZc0wP7HrrXZsAZ4Vdb0kp3/IxdCLwWU3nC5mhTU9PZ+XKlaG6fIQIERoRQohqPcZAyDmms+LHtgGNNbXY4ZRSDqhhTlbgQuCRms4XMkMbIUKECA2JBAyMYJ7yfGC1lPJITQMjhjZChAinBBKJKgMLHQTI1QQQNoCIoY0QIcIpRLA8WiFEFDAGuD2Q8RFDGyFChFMCiUQPUjqrlLIESA50fMTQRqiWUrfK5kNZxNptdG6WjBAi1FOKEOGkMQhN3UDE0Ebwy8w1G/nHnIWYFIFuGLSMj2PqtRfTOjE+1FOLEKHWSEAPkaGNFCxE8MmGjMM8PWcBpapKkctNqaqxJyeXWz6eRaSaMEJjxUAGtAWbiKGN4JNpK9bi1iqv0BpScrSomA0ZNWazRIgQdkhAlTKgLdhEQgcRfJJVWIzh4wOnCEFeaWkIZhQhQt2QyEjoIEJ4cXaXDtgt3t/Dqq7Tu1WLEMwoQoQ6IkEPcAs2EUMbwSeX9etJi7hYbGZT+T6HxcydIwaTEGUP4cwiRDg5PJVhgW3BJhI6iOCFYUhmrdmIW9cBQZzdRtfmKUwePpDhndJDPb0IEU4SgU5o0hMjhjaCF//6YSEz12ykVNUAcGsaWw8fpVNqwPnZERoxBWoxH+z5ht+y12ESCqOancG16eNwmGyhnlqd8CyGhcbQRkIHESqRW1zKF6s3lBtZ8HxAS1WN95auCt3EIjQIqqFy75oX+OXICoq0EvLVIr7LXMyj619r9Gl9njxaEdAWbCKGNkIldmXnYDGZvPZrhsGaA5khmFGEhuS37HXkuQvQKoivqFJjf/FhNuTvDOHMgoMhRUBbsIkY2giVaBkfh6p7KxwpQtAhJTEEM4rQkOwsPIDTcHvt16TGnuKMEMwoeEQ82ggNTm5xKYu272FDxuFKj4QtE+IY0r4ttiperdVs4pah1eogR2gCtHI0w6ZYvfZbFDNp9sYdo5cIdJSAtmBT42KYEMIO/ArYysZ/KaV8osoYAbwMjANKgBullKuDPtsIQeGVBb/z7m9/YDGZMKQkLS6Wd6+7lBbxsQD894rxPPn9fOZs3IaUkuZxsTx5wWi6ptWpk0iERsBZzfrx4d5vcRsqsiy5X0Eh1hzNgKTTQjy7ulMfYYFACCTrwAWMklIWCSEswBIhxFwp5bIKY84HOpdtg4D/lf0/Qpgxf+su3l+6Epem4yorsd2bk8udn33N7CmTAHBYLTxzyVieumA0papGnN0WUe06RYgyO3ih7338d/s0thd6usf0ju/MPV2vxSS8Y/eNCYnALUNzDzUaWul5riwqe2kp26ouP14EfFQ2dpkQIkEI0UJKeSios41QZz5avqZSRgF4NAz2ZOeyLyePdskJ5futZjNWcyQD8FSjVVQznu97L07dhUBgM3mHEhojnoKF0ERLA7qqEMIkhFgLZAHzpJTLqwxpBRyo8Ppg2b6q55kshFgphFh59OjRk51zhDqQV+L0ud+kKBQ4fR+LcGpiN9majJE9TlgvhkkpdSllX6A1MFAI0bPKEF8z80q6k1K+JaUcIKUckJoaifeFgnO6d/Ja6AIQAro2j/ybRGi6SCnQpRLQFmxqdUYpZR6wEDivyqGDQJsKr1sDkaTLMOSGwafTLC4Ge1lIQBECu8XME+NHYTU38hiclEhZ2ugT6yPUHwYioC3YBJJ1kAqoUso8IYQDOAf4T5Vh3wB/EkJMx7MIlt9Y47PZWQVMe2chq37fRUJiNFdcP4zh5/QI9bSCRqzdxuwpk5i5eiOLduwhLS6Gawf2pUfL5qGe2kkjpUSWTIOiV0EWgIhHxtyDEj0x1FOLEEZ4FsNCs+YQyFVbAB8KIUx4PODPpZTfCSGmAEgp3wTm4Ent2oknveumeppvvXIsu5A7rv4fxYVOdN3gSGYezz0xm/17jnLtbSNrfT63S2X+3PX8tmArCYlRTLhiIF17eIWuG5wYm5UbhvTjhiH9Qj2VoCBLpkPRcyDLdHLlMSj8N4awokRdGtrJ1YDU9oC2G8wdEeb0UE+nSRPKxbBAsg7WA6f72P9mhZ8lcFdwp9bwzPxkKSXFLnT9hFCay6ky/f3FXHz1YKJjApcHdDlV7r35HTL25eB0qgghWDRvE7ffN5bxl51RH9MPClsPH2XW2k2UuFXGdOvE8E7pKEqYp3YVv3rCyJZTCkWvQC0MrdQzwLUURDTYz8bzAFc/SOlC5t4F7uUgLCBVpG0wIuE1hGjc4i3hjB7GebSnDGv/2IOmepefmi0m1q3cw+HMPNxOlUHDu9K+c/WP2vO+XcPBfTm4nCrgebx1OVWmvvADo87vjSMq/P6YPlm+hufnLcGt6xhSMmfDNoZ2bMsrV14QtsZWSglGtu+DRuAtd4zCF6H4fUABoUCBgMS3Edb+wZloFWThcx4jiwuky7PTtQxZ+Dwi7rF6ueapzvHKsGAghEgA3gF64nGWb5ZS/u5v/ClTgqvrBls3HGTLhgPomrcxBWjeIh5feflup8q/HvmC9179mY/eXMBfbniLqS/+UO31Fv+yudzIVsRkMbF5/cGTuof65FhxCc/NW4xT08pb2JSoKkt37WfRjj0hnp1/hBCg+AnHmNr43l8F6VoGxR/iqc0pBVkMsgiZOwUpvev+g0Lpl2XXq4irbH+E+sKQSkBbALwM/CCl7Ab0AbZUN/iUMLQb1uzj6rHP8fCdH/LoXR8xcezzbFi912vcmAtO96qAMpsVdN1Adeuobg1dN3C7NL6fuZJNa/f7vWZsvO/HTmlIomPCz5tduns/ZsX741Ciqvy0eUcIZlQLYv8PqBrWsSNiHwro7bL0C8BXHzQd3H/UcXI+riclSD85y/721yMu3c2ynA38dnQtxVrT7QfnEZWpu9aBECIOGAG8CyCldJdlZPmlyYcOCgtK+eufP8FZWtEzcfPwnR/x1EvX0G9QR4QQOEvdvPrMd14NCRWTgsms4HJWrqZyuzTm/7CeHn3boqk6a1bsprjISe/+6SSlxHLBFQNZsWRHJa9WCIiNjwqLBbGqeNK9vN15RQgcVkvDT6gWKI5xSGFFFr4I+gEwt0PE3I+wnx3YCWRVzzLAYyeJEAJp6Q9qVSMuwNqwwj1rcrfxz83veK6NRJcGf+58FaOaDwzK+bcU7OGbjEXkuQsZlNyTsS2GhkxAXCJQAy/BTRFCrKzw+i0p5VtlP3cAjgLvCyH6AKuAu6WUxf5O1qQNrZSSF5/8qoqR9aCpOo//5RP6D+nE31+4ml/nbaK40OVVZlFTSubOrYd45K6P0FQdKSWaZnDNLSO45tazmDR5JB+/uQCzxYSUkphYO/98dVKD6gYUu9z8tGUHRwqK6ds6jUHt2/i8/pmd0n2WnVjNJi7pG/5iIsJ+DsJ+zsm91zEe6V4CsqTyAamBNTgGx+uacU8gj00E6QbcgAWEDRH7t3q5ni+KtVL+seltL1nEV3fMoHtce1o46lbAMvfQb7y9a1aZQA1sK9zLnEO/8dLpDxBlbvi+c1JSm2KEbCmlv289M9AP+LOUcrkQ4mXgYeBxfydr0oZ29qe/s3zJdr/HdV2yYskOXvn3t9jtVp8GGSnRfbTFtNrMjBzbk7/++WMK8ir/gU5/fzE9+rblyhvO5LyL+7Fp7X5i4xyc1qcNio/H8/pi25Fsrnv/czTDwKlq2M1murdI5b3rLsNWpcOt3WLmf1dfxB2ffgV4vm80w+DPI4fQq1Vag805JNjOBetscK8oM7Zmzxb3NEKJCcolpFEM2lZQkhHmdISlC6TMQZZ8AupmsPRARF2LMDXc73pZzgaf+3VpMD9rJde2O7/a9xeqJWS7cmluTyLKXDlU5tRdvL1rNi7jxBOdy1A56srlh8NLubT1qLrfQK0JWjHCQeBgBSmCL/EYWr80aUM7/b3F6FrNPS1//Hotl103BLvDgrO08gKW2WLiwqsGMvvTZZ52xbqByaww/rIBIMHp0rzO53KqzJm5kj4D2hMXH8WQs7oF7Z5qw71ffE+B88Sjb4mqsjEziw+XrWHycO8Us4HprVnywO38unMvTlVlaId2pMZGN+SUQ4IQJkh4E9xLkM5fQIlHOC5BmNsH5fxG8XtQ+BIIM0gNaemKSHgTYUpDxD4QlGvUhl1FB/lk7/dszN9VyRCWz1fqOHX/sWJd6ryx4wt+ObICi2JGkzoTWg7npvYXogiPI7GjcD8m4e1UuA2VpdnrQmJoJbXyaP2fR8rDQogDQoiuUsptwGhgc3XvabKG1jAM8qt4mv6QUrJ7+xHsditul4ZheDxYs9lESrM4brxzNBddNYhf523C7dYYdGYX0js154/fdvj9fiwpCX5srzZk5hWQkVfgtd+lacxeu8mnoQWPROLY0zpXe+6Dufm8uuB3lu89QEpMNLcOG8B5PboEZd6hQggFbCMQthFBPa90LYLClwHnibCUugmZdxcieXpQrxUIOwr389C6V3Abbm8xkjKsipWBSVXlTE7w8d7vWZD1B6rUUHWPo/F95hKSrPFc0toTF48xR6FL305ORkkWWc5jNLMn1eleToYginr/GZgmhLACu6mhSKvJZh0oikKLVoG3Xsk6lMfLH93G6QM7oCgCk1lh6MhuvPDOzSiKQnJqHJdcM4SrbhxOeidPDu1pfdr4TBWzOyyMGOP/g9oQyAr/9X2sdhiGZGdWDiv3HeSSqdP4dsNWDhcUsTHzCI989SNvLwn+6nwokEYRRsEzGFnDMbLOwih8GVmHTABZ/B7eGQ2ax9hqDZ/m9/6eb3BVY2TtipVByT3pGd/J53EpJd9mLvbyhF2Gm1kHfyl/nR7dkma2RJ+OSKFWwv1rX0T14U3XJ5LA+oUFIg4upVxbJpDVW0p5sZQyt7rxTdajBZh871ieeWwmLlf1/6CKIjitT1vSWibyr9evxzAMhBA1LlpFx9i586HxvPHsHFS3xxO2O6x06taCs8f2Cuat1JpWCXG0Sohjd3blf3+b2Vzrxa3lew7wwMy5FLncuDUNvcoKYamq8caiZVw7sC9RYZ6hII0Sz6KTDxFrKTXksatA24dngQoofgfpXgpJ009uEdPI8b1fWEDm4dFfajh2FvpOSRQIhqb0ZkzaYAYknlbpXrNduRhSkmpLRJcGLt13bnGheuIJUgjBk72mcP+aF8lVCyuNk0hKdRfLcjYwPLXhysA97cbDV+ug0VBa4mLRTxvJ2J9Dx24tGXZ2N/72/EQ+/N8vZBw4Rpt2KaS1SmDJ/C2VKsCkhG2bDjJn9krGXtgPkylwR/+8i/rR5bSWzJ21ioL8UoaO7MaZo7pjCgMlrP9ePp7rPvgCVTdwqioOq4WuzVO4cbD3h/tIQRFFLhfpyYmYKizYHSkoYsqnX3mJhVfFJBT25uRyWotmQb+PYCBdvyELngA9A7Ago65AxD6E58mvDNf8suMVDYkLtG2eVKyTyUCwjQRtD1D1y16CufoQTX2QZIun2IcmsUkojG8xnF4JncqN7P6Swzyz+X0OOY8CglRbIg91v5FWjlQOlmZ5naNTbOUCkeb2ZMamDWX6gR+9xrp0lUOlfr6E6o360ZoNhCZjaDP253DPTe/gdqk4S1XsDisfvvELL39wG69+fHv5OCkl8+euZ/p7izm4PxspPUUE+3YdZerzP7Dujz088q8ranXtDp3TuOuh8cG+pTrTNS2V+ffeyo+btnOksIjerVowtEPbSuW02UXF3P35d2zIOIJZUbBZTPzjwnMZ3a0jADPXbEQ3ag42uHWd1JjwXDiT6iZk7p2ceITXoeQLpJGPSHj+xDj3Ou8UL88BUDeclKEV0TcjS78CI48TBtwBsY+ERNNgYtuxvLJ9Oi4fnW6f3vw2MeYo/tX7TyRb43lo7csUasXlYYaM0iweWfcK93W9jme3flCetiUQWBULkzt660p0im2Dw2SjVK+8ZmEzWegQ07D55BICrfoKOk0mRvvCk19RmF9anjXgLHWTdSif91//udI4IQSjx/Vh0uSR2KwWZAUj4nSqLF24lb27vL+tGysxNiuX9evJnWcN5sxO7SoZWSklt3w8i3UHD+HWdUpUldwSJ/fPnMO2Ix79gMMFRbh9tB+viNVkYljH8M1QkEVvAlW9OCc4f0Aax8r3CHMbwFdFn44s/gSjeDpSev8upLYPWfIpsvQbTxpXBYSShEj5FqJvA3NPsJ2DSHoHJerKOt/XyTCy2QCuSx+Hw2TDLE74WZrUKdVdZLvy+PvGN1mavQ5Val6xXE0a5LoL+E+fuxmU3IuWjlTOTOnLi6ffR5fYdl7XG5jcgxRbAuYKoRqLMJFmT6FfYsNn44R1h4Vwx1nqZuuGg16Cz5qms/hn31kXa1fuptRH3qwQgs3r/JfWNiW2HD7K/mP5aFU8VlXT+Xj5GgB6tWzuszQXwGY2YTWZGNW1A89fVn3OZUjRduJzCVBYy0IFZdgneGKnXn9oEowMKPw3Mv9EOpaUEqPgGWT2BGTBM8iCJ5BHz0RWKdsVSiJK7N0oKbNQEt9AWEOr3nZJ61F8NuTfdItL9zomkWS78tlZdBC37r224TLcHHXl0jm2LY/3uI23z3ich0+7ifTolj6vZRIm/t37LwxN6UOMOYp4SwzjWp7Js33uLk8FayikFMHUOqgVTSJ0UN0ihT/VqdRm8VisJlR3ZQ9FMQmSUmKDOr9w5WhhMSYfvx9dSjLyCjhWXMLL85eiGZXTdCwmhTHdOnHv6GHEO+zEORq+yqdWWHtD6V6gijcqVTCll78USgwkfYbMvx+0Hd7jKQXnL0h1B8LSGdy/QelnlIvDlNlymXsHNFtaOf4bZlgUM5rhO+5uEgot7MlYFDN6lRCDTVjZlL+LS5c8gEkonNWsHze1v4hos29tj58O/87bu2YjkWiGTq+ETlzT7vzQVIZBbUpwg0qT8Ghtdgu9B3jrpiomhcFndfX5nnMvPN3noldpsUfz4Kvpy5p0S5RCp4tNmUcocXt7LTazmTM7tuP9pav9Nmz8x0Xn0iYpIfyNLCCip4BXPNQBUdcilMpfqsLSGSXlG4936/tsoK4GQJZ+6UMHF8CoFzGaYDMstS9WxTtLRCIZkzaYLrHtKh23YEZHZ0vBHlyGmxLdybzDy3l43Ss+/1bW5+3gzZ1fUqI7KdVdqFJjfd4O/r35vXq9L/80kp5h4cz9T1xMcrM4LJYT31iGbjB/7nq+n+n9oU9pFsdTL11LUkoMVltlxz77SAHvvfoz7732s9f7mgKLtu9h6HNv8srC371EdCwmhYQoO1f068miHXtw695J5zazmV1HG3rFuHZIKZHuNcjSOYAJkfQZWIeAcIDSAmLvq17dy9we8OGRCgWUsswKWV3aYMPmiJ4M41oMo4U9BZviuU8FgU2x0D+xOzeveJItBXtItsbTzJZES3sqZySfhlmY0DnxmdCkTqYzm3V53qXuMw/87JVvq0mdzQV7OOqsNu20XvAshgUnj7a2NBlDm9o8nqdeusbrm1V167z5/A8cPZJfvq8gv4SP3pzP2y/9ROfuLWnb3ls8w+VU+eqzZZQUh7bCK9gUOl385fPvUH0YUIALenXjq9snEeewkxIT5XOMquskRtVf94G6IvWjyOxxyNybkAWPIbPHI0s+QCS+j9J8HUqzRSjRN1QbchKOS8Er11aAiALbmWVjJnhee08ALPUjRhNM7CYb/z39fm7veCmDknoytsUQ+id2549jm8hTC1GlxiFnNgVaMU/1uoMEa5yXAA2AbujsK/FuEXjU5Vs50CxM5FXJrW0ogiGTeDI0GUMLsGzRNt9qWwKWLtwKeIzsHRP/x+cf/saOLZksX7ydnVt995E0KQqHMhr+m7c+Wbh9t5cXW5HkmGgSoz1G9OahA7BXEZ8xKwo9WjandWJ8vc6zLsi8e0Hf60nVksV4BLXnIkumI40ipM/H/coIUxoi8e0y79UB2MDcGZE0DSHKHqdt55Z5yceNrQWwQ/y/EYrvL6lww2ayMrbFUP7WczKT2o3nj2ObvbxQVVeZefAXfs9e5/McqtRYfWwrRVrl1Li+iV0rZRscR5cGbaIavhloMCvDakuTMrRSSqS/4sIy4zL709/JzytGdVefgA/gcqukNo8L5hSDQk5RCav2ZXCkoCjg9xiGxK1puDSd6opw31u6kvu/nIOq65zZqR33jR6Gw2ImxmbFZjbTq1VzXrvqgiDcRf0gjWOgrsV7IcvpyRrIGog80h8jZyJGDSWwwjoQkforImUWImUOSsp3lRooCmGCmIc88VxLP3BMRKR8j+IYF/T7aggOlmZhVbzXx3UM1uRu9cqFrcia3K08tPZljAr6Bpe3Hk2UyY6pgpmxKVauSx+HPUSatAZKQFuwaRJZB8cZNuo0pr+/GLde2YhKSbmC1oolO7wyDfyRlBJLXHz4eCa6YfD3737h6/VbsJlMuDSdUd068uwlY7Gaff9TulSNf/+4iNlrN6HqBu2SEqrV2NUNyS9bd/HSL0v5v3OHc/3gflzerxfbjhwl1mZl1f5MHp79I8nRUVx9Rh96tw4zCUVZin//oYKhUFdD9miMmCdQYq7xezohFDB39HnMKP4ACl+g3Kirm5CmZoiY232OD3fS7Mm4fWQiKChEmewc8VdOjMcYH3YeY3XuVgYkeUq8k2zxvNb/YWbs/5E1edtIssZxWevRDEwOjQ6IlKAakYKFOpPesRlX3zwCq82MyaRgNitYbWZuu/tcmrVIACAxOTB9USFgwqUNq3ZfE28v+YNv12/FrekUuty4dZ0F23bzws9L/L7nvplzmL12Ey7N03BxT04uAjBVE590ahrTV64vfx1ltdCteSr3fTmXZ35cxK879/L1+i1c/+EXzKgwLixQWoISqJiQhKJ/IV2/1voyUjtYZmRdgFa2uaDoNaS2u9bnCwdSbAkMSu7llYlgUcyc1WwAVlG9joVqqOwtzqy0L9kWz52dr+TtMx7nP33uDpmRheOhg9Dk0dZ4RiFEGyHEAiHEFiHEJiHE3T7GjBRC5Ash1pZtDScTX4Vrbj2L16dN4fopZ3PDnaOZOuMuLpo4qPz4ZZOGYrPXLHwiJcz4cAkfvP4zuubpnrB1w0Gmvb2Qr2csJy/Xb9eKeuOj5WtwapU9DpemMWPVBp/pNZl5BSzZubcsXHACA8morh25op//D32J6q50zi9Xb2R/bl655oEhJc4yb7nYVU8NDE8CIQQi/j944qrHvfzqcifdyKKptb+Q62d8h2B0cP5U+/OFCQ90m8TYtCHYFAsCQbuoFjzd6w4ubDUCn51LK2A1WWhZx64M9U2oKsMCCR1owP1SytVCiFhglRBinpSyasnVYimlv+TDoONyqiz8aSNb1h+gTXoyYyacTlyC5zG/bftUn5kEAKcP7MDNfzqH9177GZNJwVnqLtefrUppiZtZny4jL7eE0hIXyxZtw+XSsFhNvPvKPJ54fiL9h/iWk6sPipy+DZpL1dAMA4upskHZn5uPtSzEUBHdkBwpLOLViRewIyuHtQe9FwN7tUyrtCo/b+tOnD6EZTTdYH3GYYZ0aHsyt1QvCNtgSPkWWfIZ6Ps9KV3On/Auwy3DOHwyVznJY+FLRmkWsw78wq6iDEamDuDi1mfTNvpEaGhwci8WHV3l9/1RJnu1Orah5nh6Vyio0aOVUh6SUq4u+7kQT1vdkHYXzMst5rbLX+ONZ+cwd/YqPnxjPjde9BJ7dhwJ6P0XXz2Yz39+kNvuPtevkT2Oy6ky79s1/L5oG06nipQSt0vD5VT5x8Ofo9agahVM+rZp4XN/l+YpXkYWoENKopeRBU/mQJ+y2Orj484mymopDyWYhCDKYuHxcZUbGzosvp8CNMMIy5xaYW6LEvcQSuLriPh/gpLiZ6QC1oFIqSOdP2Dk3o2R/xjSvbb6C9jH4NugmsA+to6zb3i2F+7jL6ue5afDy9lRtJ+fjyzn3jXPs6OCrOLo5gOxK/6r3VRDq3bBLPSEceig0jSFSAdOB5b7ODxECLFOCDFXCNEjCHPzywev/0LO0YLyHl8ul0ZxkYvnn5gd8DnsDiuL51fbfaIcKanUzfbEAdi0puF0ER49b2Qlo6gIgcNi5onxvtuCNIuNYUKvrpVStAQejYIbh/QHoEfL5sy+fRKX9etJz5bNuaxfT2ZNudarT1j3NP+PhJsyA/uCCxVC2BApX3oyAyqhgIiG6CnI3MnI/IfBNRdKv0Qeux6j6G3/5zS1hNiHARuewgar5+fYeyplJjQW/rfzC5yGG6OsGEHHwGm4mbprZvmYfond6J/UHZuPajIAp+Hmx8NLG2S+J4tR1jespi3YBJx1IISIAWYC90gpq/ZIWQ20k1IWCSHGAV8BXmKbQojJwGSAtm1P/lHztwVb0Hz0Atuz6wjFhU6iYwMrCz24NzugcdXlnTYk3dJS+WrKJN75bSWbMo/QpVkKtwwbQKdmyX7f8/SFY2iblMAny9dS5HLTv11LHh57Fq0STqSttUtO4KkLqu8ge1qLZlgUBbWK7oEA7H683XBCKEmI5OlI5wJk8VugHwHbIET0naBtRbpXAcfzQCXghKJXMOwXIdyLkSUfe/Jy7WMR0bchlDiU6GuR9pFlYQnpUeYyh08IJVCklOwoPODz2LaCfSU8BvsAACAASURBVOU/CyF4uPtNfLT3O2YdmF+pQgw8/cA25e/m8jZVzxIeeLIOQqN1EJChFZ4M7ZnANCnlrKrHKxpeKeUcIcQbQogUKWV2lXFvAW8BDBgw4KStl7kaUW2TOXAnvVXbZLIO51c7xma30Lt/OhtW7/Vq3CiEoMfpDfuHlRYXwxntWmNSFNolxpPsp3rrOCZFYcqIQUwZMcjvmE2ZR3h+3mLWZx4hJTqK24efwSV9e/Dt+q28tuh3sgqK6ZCS5OngW8XQ2ixmLulbrw8wQUXYz0bYK4dFjKJXOGFkK2KC/PuR2voTmgbF7yGdcyH5G09RgpIGthEgYhCmhk/CDwZCCOwmq8/Hfoe5cr6rIhSGpfTl24xfvQRnzMIUkkKEQDlesBAKajS0wrMi8i6wRUr5op8xacARKaUUQgzEE5Kot8DdeRedzpefLMVdoQOtyaTQZ0B77I7AFZOuufUsNq8/4BUWiI6x4SxViYmzc+k1Q7j8+mG8+ORX5Z0ZzBYTQsDjz12FxdJwqcgFpU6ufGc6WYVFlLhV7GYzry1axkc3XnHSnQ22H8lm0vufl2cTFLvcPD1nAfO27GTZngPl+zcfzsJiUrCZTZgVBUNKdENy+/CB4ZdLWw1S3Y4setlT1GBq6fFoRSyej2zVpySnp7NCpf0q6AeRJTPB3AqZ/6hnnNSQlp6IhFcQpvDsMlEd41qcybeZv+KuUBVmUyyMb3Gm19hOMW1oFdWMfcWH0Cro85qFiQkthzfIfE+W+ggLBEIgVmIYcB2wQQhxfIXgUaAtgJTyTeBy4A4hhIZHxn6irEfpq6tvGcGmdQfYtvEgEo8UYkJSDA/8/ZJanad3/3Qe+PvF/O/5uRQWlCKEYMyEvkx54Hys1sq/mgefvoxLrslk9bJdRMfaGTGmR4MXM7y+aBkZefnlOgVOTQMNHpr1A9/edf1JnfO1hb97ZROUqhoLtu/xGqvqBn1apTFp0OmUqirDO6XTIr7xSEpKdTvy2JVl3qkE4ygy7x6IvhlPjLVqVoK/VvUGFE/F8yBX4T3qOmTuTZD83cn1Fwsh16WP56grl9+z12NVLLgNlaEpfbi2nXeVmxCCf/S6i5e2TWNV7haklESbHfSI78COwv0k2+Ix+Si9DTWhzDoQoZICHDBggFy5cuVJv19KybZNGezefpjmLRI4fVAHz6PtSWAYBgV5JURF27DawjfeOOKFt8gq9M7ftZpMLLj31hrDCL4Y9d93yMwPXOAj1m7jj4fvrPV1wgEjdwq4FuCV/yriIfovUPQP72PV4sMLFlGIpI8RltA25zxZclz5ZJZm0crRjCRbzXoWK3I28cyW9zCkRJUadsVKm6g0/tPnL9hMwdPjFUKsklLWqYIoqXuqHPPeZQGN/Xzo1DpfryKNtgRXCEG3nq3p1rPuXUQVRSEhKbCKsVDir9OBRFZqqCilZH3GYTZkHKFlfCzDO6f7TP8Cz0JYbQxtm4Tw034IGHUdPg2pdCEcY5BFz+E319YnvjxeBfQsj75MIyTZFk9yAAYWPOIwL2//tJIIjdNws6/kEN9lLuayNqPra5onhZQCLUipW0KIvUAhnvprrSaj3GgN7anI5f16MnXxH7gqVIcpQtC7VRoJUZ5MC7emccdnX7N6fya6lFgUhRibjWk3X+lTceuOEYNZc+BQpfCB3WKmQ0oSe7KPVep+a7eY+cuoofV4h/WMkuan/bcEJR6UJDAyfRz3hQPP31iVIhKpQiP1ZmvL/uJDPmUT3YbKwqyVYWdoIeihg7OrLvj7o0lpHTR1bh02gH5tWuCwmLGZzURbLTSLjea5S0/063r/99Ws2pdJqarh1nSK3SpHi4q5f+YcDuUX8tjXPzHyxbe5+H+f8M26LZzRrhXPXXoeaXExmBUFu8XMlf178dnNV3H78IHE2qyYhKBZbDT/vHAMI7t0COFvoG6ImDvxbr5oA1MbZNbZYBQQuO/hxsvI4oCoiYDAKHwF49hkjMKXkHrom30WqsXsLDzgJWVYF8yK2W8XEoufXNtQEkrh70bh0RqGgdulYbNb6mWRweVUEYrwWgALN6xmM+9dfxkbMo6wIfMwLePjGN4pHXOFljxfrt7opYdgSMnmQ1lc9L+PKXa50aXkcEERT3z3M7uO5nDvOWdyTrdOFDpdZOQVsGTXPj5buY4ilxuzyQRCo0V8LG0SExr6loOKsI9Bxj4IRS8CmkegG7NHu5bjvzMFj8B3tKfNOCp+NQ0qnx2ibwX7ecjs80C6ADe4lyJLPoKkGZ4+Yw2MLnVe3/E584/84ekTJjXGpg1lcsdL69QcMct5jNd3zPDZttymWBnXYpjXfiklWwr2sDxnIw6TjZHNBpDm8J8DXh/UwoimCCEqLiK9VZaeehwJ/CSEkMDUKse8CGvLYhgG095exKxpv+MsdZPcLI4p95/HmaNOC8r5D+zN5sWnvmLrxoMIIeg/pBP3PX5RwApf1ZF3rIjdO47QvEUCrdoG78MkhKB36zS/KVWan84JhiEpcXuM7HFKVY0Plq3m5mEDiHfYmbr4D6atWItmGBjSoGJ18rqDh7nhwy/44rZrqi2QCHeU6GuRUVeCfgTp/h0K/wmy4heTAQiwXwwxf4acCWAcDeDMwqOrUPg0yCJOGGc3SBVZ8CQi+ZOg309NfLz3exZmrUSVGmqZfOhPh5eRbE3giraVi1R0qfPlgV/4LvNXSjQXfRI6c0vHi2nlqJyullmSxZRV/0KX3p81izAzJKU3o5pX7vQrpeTFbZ+wNHsdTsONWZiYceAn/tJ5Imc3b5iuwLXMo82uIe46TEqZKYRoBswTQmyVUvqVgQvr0MH7r//CFx/9RkmxC8OQHD2cz7OPz2L18l11PndxoZN7b3qHLesPYugSXTNYtXQn993yLrofYxUIUkreeG4Ok8b/l388OIM7Jr7B/01+n+LC2iyynDzjenbB6mPhy6QoXm3FASwmEzuysll74BDT/liLU9PKDK33uV2aztTFK+pj2g2KEBaEuTWo2zzVXl5IKJ0Brl/LymwDqTQ0wMgqa8pY9ZcnQV3Z4M0+pZR8l7nYq2OCy3Az6+AvXuNf2vYpM/b/yDF3AU7DxYpjm7h39QvkuCoX9Ty6/nWfRtYkFO7sfCX/1+16L295Ve6WciMLnt5hbkPllR3TKdZq7ngRLIJVgiulzCz7fxYwG6i2d1HYGlq3S+Xr6cu9iglcTpWP31xQ5/P/PGcdbrdW6cOv6wa5OUWsWX7yeqLfz1zJD1+tRnV79BdcLo0t6w/wwpNf1XnOgTBlxCDaJMYTZfXEyOxmT3eEM9q18vnxUXWD5rExfLdhK64aBHIMKdl8OPTxxqBhTsc7ZnscFYqeQ9gnIBJfAXMvEAlgGYDvlAI72M4G4ccoC1uD59YaSJx+RF4KtGLW5Z5oqJjtymXx0TWVjLJE4jLcfJ2xsHxfRmkWOW7fvcAMKSsVPFRkUdYqnwtnJqGwNndbILdTZ6QEzVAC2qpDCBFdpmSIECIaOBfYWN17wjZ0kJ/nP2ifeeBYnc9/cG+2T6EYXTPIPJgDnJz84exPf/c6r6rqrFiynZJiF1HR9dvCI9ZuY/aUSfy8dSdrDmTSJjGeC3p352BuPqsOZFbKLrCaFPq2bkGbpASqaQJUjiIEXZv5U8FqfAjHhciil/y0DAeMbMCNsI1E2Eae2F30Pyh6E09tDngW1JojHFcg9YNQMoNK3Rywgf3SermH6jAJhbZRLXw2TgR4fttHfDjoKRShsK/4MFbFXB5eOI4mdbYW7C1/XaiWYBImDOn7S7mNw3cJrlmYEHj7+gJRp1hxbQnSQldzYHbZF6cZ+FRK+UN1bwhbjzYxKQaTyff00jvVvcSxc/eW2B3enoliErTvfPIlpUX+QgRCUFrSMALZVrOJcT278tj5Z3P94H4kRjno1SqN5y49j+RoB3aLGavJxJmd0nltoqf/14Re3VBqEnY2m7h9eOUnJClljVKT4YpQ4jxtyH21FQcQcT6PKTF3IBJfBesIj6cb8ydE8myEEo2IfQCsZwB2EDGe/1v7I+IerMc78c8dnS73+yBcrDnJLPXEn1s4UlClt6SmCYV2FTRp20e3RPFzxjhLNL0TfC/4jW4+EKsPiUVDGpye2K2GuwgOwWrOKKXcLaXsU7b1kFL+s6Zrh61Ha7aYuObWs/h46oJKHqLNbuGmu04E8Z2lbvbuzCIhOZq0loG2MIGzzu3JR1MXoLoLy2OyFquZ9E7N6dn35IVi+g/uyIIfNngZn4TEKJJSQlsUMaZ7Z0Z17cih/EJi7TbiHScec09v09LzaOsnjti1eQqPjzubrmVyiSVulWd+XMTX67bg1jVOb92Sv08YTZfmjcvjFZbOyPhn4LhmQTkOiLnD5+O+lCrSvQ60zR5vWN0ARi4oMQhhRyS9h9R2grYLzB0R5oYTh69Kr4TOtLSnkuH0XtCTGNjKjF9LRyo94jqyMX8nagVv1ayYubj1CREem8nKbZ0u5a2ds3DLiroIVl7oe5/f8EjPhE5c2GoEX2cs8qQtC4EEHj3tFuxBrCCrCRmuojKh5PLrhhIX7+DTdxZxLKeIDp2bc9s9Y+nWy1MNNnPaUj58Yz4mk4Km6XTt0Yq/PTexvNNCddjsFl77eDLvvDyP3xZuwWwycc6EPtxwx6g6xdJuuHM0K37bgbPUjerWURSBxWrmnr9eFBb17yZF8dsq3N/sFOCrKZMqzf+OT79izYFDuHWPF7T6QCbXvDeDOX+6gWax4V9lVxHFMQFDOqHoBTDyPJ5ozBRE1I0+x8u8+8C1iHLD7PoZ6V4BqT8glCQAj3ENoYGtyKVtRvPWrlmVUrEUBG2jWpBqP+Gc/LXHLby580sWZq1Cx6CNozl3db6S/cWH+Xz/PGLNUYxqPpDzWwyjbVQaXx1cSI4rjzOSejCh1QhiLdX/3d3Y/kJGpPbnp8O/E2+J4YJWI4gxN6xeSKhEZRqt1sGKJdv5x0OfV/J2zWYTPfq25dmpNwZhhidPbk4RX322jPWr9tK6XTKXTRpKeqfwlY87zpRPv+LXHXsr6e8qQjCsYzvennRCsGf7kWyufOczLzEaq8nErcMGhHX1mJSaxzC6FoBIQkRdgTB3KDsmkTIfSuaA63sQ0Yioq8E2svxLRmr7kNkTqByDBbBBzO0oMX9q2BsKAEMavLRtGouz16KgIATEmKP4T5+/0NzunaqnGTqa1CjSSrh3zQscc5+Qn1YQjGo2kHu6XlNrx+HbjF95f8/XmIQJXRq0cKTwZM/bSbHV/CQaDK2DmC5psu8bgYkv/TbmuYjWAcCXH//mteikaTpb1h/g6JF8UpsHVq9dHyQmx3DTn6oX0g5HHh83iive+pRSVaVU1XBYzNgtFq8ODntycn120XXrOpsPB5JzGhqkdCOP3QDalrK0LjOyZBoy/t8ojvGABrmTPWlfZQtd0r0coq6G6BvBOIbU9oCwlBUkVMQF7nUNe0MBogiF+7pdx1WlY9lWsJckazy9Ezr7XYQyKybMmPjrhjcqGVnwZDIsyPqDoam9aeVoxuKja9ClwdCUPnSI8d/han3eDt7f83VZVoPn73Z/8WGe2DiV1/s/HLR7rR6BHqJ2443W0B7LKfK532RRKMgrCamhbay0Sojjp7/cxNfrt7Dl8FG6NU/l4j7dibFXzpTolJqE7mMBzGYy0bNlGHvupV+fiKsC5W3CCx5D2keDcz5o2zmRTYDn55L3kSWflBlYg+OGojIWsHSp7zuoE60czbyKD/yR7y5ke+E+n8d0DD7c/S2HXDnoho6BZNbBX7iw1UhubH+Bz/d8nbHQK5/XwOBQaTb7ig/RLtp3P7xgE4nR1pIBQzpx6MAx75Y2Etr46YB7KuNUNQwpy/Nr/RFjt3HtwL7VjumYmswZ6a1ZsfdAefNHgafbwsQB4SuoIku/9ZPKpYB7LdK10H8BA+6yktyy8Zg5UbYLCAsi6togzzh0uA0NUU08c19p5c7BLkPlm4yFjEjt59OzzXP7VogzCYUC1bfTFGzCugtuuHLlDWcSE+fAbDlRBWWzW7j9/vPCXrOgITmcX8jNH81kwL9fZ+Azb3D1uzPYm5Nb5/O+etUFTBzQhxibFYuiMKxjO2bcejUpMdFBmHU9ofhbeJGeQgMlmcB8D6NsrAUwgbk7IukjT8PGJkKKLYEka+2eClVD47ds392DByX3xOpDaEaXOp1iG6gdlPQk1QSyBZtGa5GSUmKZOuNOvvx4KauW7SK1eRyXTRpKnwHtQz21sEHVda5+bwZHCorKF7jWHsjk6nen8/PdtxBtO/m0GrvFzCPnncUj550VrOnWOyJqokffoKpXK2LA0huhxCFLplHJU/WHkoRInQ9oCOGvuqzxIoTg/7pdz4PrXg6glKXsPdUUH4xvOZwfD//OMXdBefWYTbFyU/sLcZjqt4inIuHcyiZsSUiK4da7z+XWu0/u/Tu2ZDLzk6UczsilzxntueTqwY1CADxQFm3fQ0Gps1IWgQScqs6cjdu4on/4PubXC9azwHEtlHwMwsTxEIBInIoQCpg7IOP/AwWPAqJM3ctXqMEKtuF4epaGnxxgsDgtvgNj0gbx0+FlAY03KSZGpFZt6e4h2uzg1X4P8v2hJSzL3kCCNZaLW42kl58Ch/pARhbDGp7fFmzmP4/NKtc72LntEHNnreKNz+4gpVkj7iJQgQO5+eV5rhUpVVX2HfNdrx4I+aVO9h/Lo1VCHEnRDZsHWReEEIi4B5HRk8C93KNdYBuGECc8e8VxPtI+ytONQTiQrj+g6GVOGFwLKHGIqJtCcg8NzSWtR7Eoa1WlhSyBINrkwGW4PYUHEhCS69PHV9sFN8rs4Io2Y7iizZgGmLlvQpTNemoaWl03eOVf3+FynfjwqG6dIt3JtLcXcvdjF4ZwdsGje1oqFpOpvJnjcaKsFnrUsmvunuxcvli1gUU79rDvWB52ixm3pnN+zy48fcEYrNW0gA83hKklOPw38hTCBlZPqbGw9EKa2yNL3gf9KNhGIKJvQ5gar1RkbWgblcb9Xa/j5e2fYSAxpEGqLZEnek7GbrKxLHs9ujQYlNyTZvakUE+3RiJZBw1I1qE8nD50B3Td4I/fdoRgRvXDoPZt6JCSxPYj2eWerUVRSI2JZnS3wKuW5mzcxqNf/4RL1cqjdUUuz+/vx007SHTYefi8kUGeff0h9Sxwfoc0ChC2YWAZUG3yvbCfDdb+yNJvPCLh7t+R9rGVPOGmzLDUvgxK7sXu4oPYTTbaOJqX/77GtfRuRx6ueBa6Ioa2wYiOsfvVnA2kfLexIITgwxsu55UFv/PN+i0YUnLeaZ25d/SZAXugTlXjsa/neVWBlR/XNKav2sCD556FooS+xLgmpHMBMu9ujqdsyZL3PeIwCS974rS+3qPtROZM9PQDoxQpojzhhOQvEErg+hr1xRFnDl8c+JnN+btp6UjlijZj6BrXLqjXMCsmusQG95yhIFTpXaekoY1LiKLvwA6sWbEbTT0Rw7TbLVw+ybsFR2Mm2matU3bA2gOZmGowoG5NRzN0rEp4f5ykdCHz76OSeIws9Qh8u34E+/m+35f3IMhCykX+ZAnoh5CFLyPi/17f066WjJIs7lnzPC7djY7B/pJDrM7dwsPdb2Jgcs+Qzi0cCVWMtsYlOCFEGyHEAiHEFiHEJiGE1xq/8PCKEGKnEGK9EML30mMY8dA/LqNbz9bYbGaiYmxYrGZGnteLed+tYcLgp7j87Gd477V5qDWIYTd1bBb/DfiOk56cgNUc3kYWAPdKfEvnlCJLfAuzS6MQtK14K6mq4Jwb5AnWng/2fEup7kIva30u8RQPvL7j8wbv6BDuSASGoQS0BZtA/jo04H4p5eoyVfFVQoh5UsrNFcacD3Qu2wYB/yv7f9gSG+fghXduJmN/DtlZBcTFR3HfLe9SWuJCSlDVUmZ/uozM/cf467NXhXq6IaNPqxZE26wUu73LTo9Xg/1tfPi1lfZNNeES4Tkm9cOgZ3jkDZUEqvVFROi/XDbk7/CZ55qvFpGvFpFgjQ3BrMKXUH311Gi6pZSHpJSry34uBLYAVWvsLgI+kh6WAQlCiIYpXq4jrdom02dAe+bPXYfq1io9WrhdGssXb+dwRt0rqRoriiKYes3FJDjsRFut2M1mFAEJDjtje3Tms5uvYnD7NqGeZmBY++HzIy8cYL8AI3cK8ugYZO5kZNaZGAX/AhEF1gE+3mcDR8N3TahKnMV/3ndDFgI0CsoWwwLZgk2tvpKFEOnA6cDyKodaAQcqvD5Ytq9SDw0hxGRgMkDbtg1Udhcg2zZnoqreOacWq4l9u4+S1ir0ix71RU5RCYoQJEb7rnDq3qIZv94/mSW79pJf6mRQehtaJjS+XGMhrJDwOjLv9rJgnQaYwD4enAvA9RvgOqHMVTIDaUpHxD9bthiW71kQE2Ywd0PE3BXCu/FwWetRTK2iNWtVzAxP6YetAQW1Gw3hnkcrhIgBZgL3SCkLqh728RavWyrrff4WePRoazHPeqdjlzQ2rdmPplU2tqqq07pd08yZ3JGVzQMz57InOxeJp4vCC5eNo11ygtdYq9nEqK4dG36SQUbYBkHqr+D8EWQBWIeBuT3ySH+gaspfKZS8i4i+BlJ/Btdi0A+CpTtY+oeFkPu5aUM47Mzhq4yFmIUJ1dDon3gad3a+MtRTC0vCOr1LeGoNZwLTpJSzfAw5CFR8fmwNZNZ9eg3HRRMHMXf2qkqG1mo106d/e1q1bXqGtsjl5tr3PqfQ6Sr/RtyUeYRr3pvB/HtuwWYJffyxvhBKHERdUf5aGrn4dXUMT6ttIcxgP9v3mBAihOCG9hdweZtzyCjNIsWaQJItIhHqCwkYRpiqdwnP1/a7wBYp5Yt+hn0DXF+WfTAYyJdS+m69GUTyc4tZu2I3B/Zm1/lcaS0Tee6tm+jaoxVCgNVmZsyFffnrs03TM5i7cRuqrlcyLx4dBJWft+4K1bRCg0gAxVelnABrWK/plhNtdtAltl3EyFaHBKQIbAsygbgtw4DrgA1CiOMaaI8CbQGklG8Cc4BxwE6gBKjXQnApJe++Oo+vPluG1WpGVXU6dWvBk/+9hrj4ky846Ny9Ja98NBldN1AUERaPhvVFRl4BpT5S11yaRmZ+1chQ00YIAfFPIXP/hKdFjQFYQNg9XW0jNBnCVutASrkE/337jo+RQIOtDMyfu55vZ6xAdeuobs+j/vZNGfznrzP556vX1fn8/tqcNyV6t0ojymqhpEraltVsplerk2+33lgRtuGQ/Bmy+B3Q9oDldETMrU1KY7Y25LkL+fzAT6zI2USMOYqLW5/NWan9Gr/zEURDK4QwASuBDCnlhOrGNspA3JcfL8Xp1S/MYN0feyjIK2lSZbT1xVld2tM2MYHd2cfKdRBsZhNdm6cwKL11iGcXGoTlNESCv+jYqUOhWsyfV/+HArUYTXo+G69u/4zdRQe5ucNFIZ5dXQh66tbdeNJda0zBaZSuW1GBL41QUEwKRYVOn8ciVMakKEy7+UpuHNKPtLgYWsbHcuuwAbx//eUBey2GITlSUOTlFUdo3HyXuZhCtaTcyAI4DTffZCwiv4HaztQbMsCtBoQQrYHxwDuBXLZRerQDhnXix6/XoFfpF+aIspLWyjs1qSaOHslnxZLtmM0mhpzV7ZTxiKNtVu4750zuO6f2Ckzfb9jKP+cupMStYkjJuJ5d+fuE0dibcLbCqcLavO2o0jt+b1HM7Co8QL+k7iGYVRCQIIOXdfAS8CAQUOldo/yrmHTbSH6bv5WSIieqqiOEwGozc/djF6IotXPSZ05bygev/YIoW/x67T/f8+DTlzJ8dI96mn3jZ/meAzz2TWVFr7mbPFkML1w+LoQzqz+k1JHF70HJR578W8sARNwjCHPgcpONhea2JDYjMKq4dro0mkBWQ8CGNkUIsbLC67fK6gAQQkwAsqSUq4QQIwM5WaM0tMmpcbz1xV18PX0Za//YQ1rLRC6bNJRO3WpX9bt3VxYfvP4Lbnflb+9nH59Fn/7tTxnPtrZMXbzCSzbRpenM27qT3JJSEqOaXg8tWfCkp1358U4L7iXInCsg+VuEuWnFtC9qPZIl2WsqdVUwodAmqjnp0Y18cTDwxbBsKeUAP8eGARcKIcYBdiBOCPGJlHKSv5M1SkMLkJAYzQ13jOaGO07+HAt+WO9VCQagKApLFmwhJ6uAObNW4XapDBrehZv/PKbJtLmpCwdz833utygKRwuLm5yhlXoOlM6icuWYBOlCFr+LiH8iVFOrFzrGtOb+rtfx2o4ZuA0VXRp0j2vPw91vDPXU6k4Qsg6klI8AjwCUebQPVGdkoREb2mCgqTrS8P7NSyn5evoyMg8cw+3yeG4Lf9zI6uW7eXfmn4mOtTf0VMOKfm1bcjCvoFLTRwBDStom1T5GHvbou0HYQFYt0dU8vcWaIMNS+zI4pTeZpUeJNjtIsjYBB+N4wUIIaJRZB8Fi+OgeWG0+es1rBpn7TxhZ8LS5KSl28eO3axpyimHJHSMG47CYUSpkJzgsZu48a3DTXAwztT4hNFMJBSwN18W1oTEJT7igSRjZMjztbGreAj+fXFhTDi2c4oa2W6/WjL3odGx2C0J4ChWsNjMDhnTA8NE91uVU2bL+gI8znToczi5g05ZMHhk6nFFdOpAcHUXX5in886Jzue3MM0I9vXpBmFqA7SygiuygsCGibg3JnCKcJIYIbAsyTdD9qB13PTie0ef3YfEvm7BYzcQpBu88Nh29TUuoUiFmsZpp2z41RDMNLVJKXpuxmM9/WoNJUVCEwGxWeP/hy+nSrnYddRsjIuFFZME/oXQ2oIGpHSL+KUQT9mibIiJcS3BPBbr1ak23Xq1xO91c3vwW1EInolkK2K2ICuliZrPCuEv9LUQ2bZau28PMn9fhVnXghLd/7/Oz+fbla7P6OwAAIABJREFUyY2iMWNdEMKGiH8KGfc3QEWIprXgd0oQYDFCfXBKhw6qsmX5jvKqKLlzHxQUIQ2JlBK7As+9dRPJqadma5DZ89dT6vKuACtxutm0u96F2sIGIcwRI9toCVC5K9QdFpo6Vrv1REM7XUfuOQhCgBD0Or8vnbs38hzCOuDLyIJH+crlOrUbWEZoREQ82tDT9YyORMVWKVKQErvDwvjJY0IzqTBh7JBu2K3e38uGIenZuVG0h4sQwaOAGcgWZCKGtgKKovCP7x4mLjkWR6wDe7QNq93C+MnnMHhC/1BPL6Scf+ZpdE1vjqMsHc6kCGxWM4/dOga71TtFLkKEsCPMhb9PKTr1bc/0jKn88cNaCnKK6Ht2D9LSm/6qek1YzCbeePQKfl29kyWrd5MYF8WFI3vSrkVSqKcWIULARLIOwgiL1cLQC5tmTmhdMJsURp3RhVFndAn1VCJEODkiMdoIESJEaJpEPNoIESKcMkRCBxFCxuLVu5jx0xoKip2cPaATV5x7OjEOW81vjHDKUqy6OVRcQFp0LDGWRvJZkdRLeW0gRAztKc5bM5cybe5KnGW5sHsycvh+yWY+enoSUXZriGcXegypsb/oV/YVLcCiRNM5/kJS7aeFelohw5CSf69awEfbVmMWCpo0mNSlL4/2H4WplqL7ISHi0UZoaHILSvj4uz9wV9Dkdas6R48V8e2iTVw19vQQzi70GFL///bOOzyqMvvjn3futPRKCST0XkQ6iCJFEFCEtfeKve5a+K3r2rGsLvaGylooNlBQUFBEQJDee0sggXTSp96Z9/fHhJAwk5CEmcwE7ud55nkyt547mTn3vec953v47cg/yLXtQJVWQHCwZDG9E+6ke9x1wTYvoORYS5m9dwsHi/Pp1ySZy9v3IMJg5KMda5ixZyN2l8pxPbNZezcTawrjwXOGBNXm2hCs0MEpb0FCiOlCiBwhxPZq1g8TQhQJITaXv572v5mhQ2ZqNge2pKE6G3811I4DWRgMitdym0Nl5ZaDQbAotDhcuqySkwWQuKSNTfkfYXMVBtW2QLItP4vhP0zj/W2rmJe6k5c2LmXkvI/JtZbx8Y61WF1Vv/tWl8rHO9YGydo64qfmjHWlNmP9z4Axp9hmhZTy3PLX86dvVuiRk57Hff2eYFKPf/D3of/mquaTWDF3TbDNOi3iosNw+xA+1wlB07jIIFgUWhwq/aOSkz2BDj1Zlg1BsKhheGzlAsqcDuxuz5OOVXWSZy3j9U3LKHL47kBd7LTzwLIfsKoh3hE5VB2tlHI5cMz/p/Y/+zensmLOao7s96/IiZSSyaOe58CWQzisDqwlNkoLynj15rdJ29F49Wm7tWtOk7hIL+Utg0HhqlG+wwZWm5M12w6xZe8RXO4A1CqGEEZdJKKan4hed2b2kyu0WzlYnO+1XJVuFqfvo0tc9cU7v2bs49GVCwJp3mkhZO1f/sZfMdrBQogtwFE8/XN2+NpICHEXcBdAq1at/HRqKC0s48mxUzi47TCKXofqUBl4SR+enPUIej8o/u9as4/8owW4XVUdi9OuMu/dn3n4g7tO+xzBQAjB25Ov4LGpP5CRXVgxmfHErSPp7KMa7peVu3jlf7+h0wmklISZDLzx2OU+tz0T6BgzngMlv+CStirLhVBICj8z5TINOu9Q0nHMip5n+l/ELb997RU+ALC7XCxJ30eBzUqcOUQVzhpx1sFGoLWUsrS8K+QPgE815PJ2vdMA+vXr57f7xht3f8T+Tak4K3WzXbtwE1+9+gM3PnXlaR+/IKsQ4UNv1e1yk5vuffc/XYrLbLw1cxm/rd2D2y05v3c7Hr1pOImx/n+cT0qMZuZLN3Mo8xilFjsdWzVhV2o2z374M6UWOyMGdGL0oM6kZxfy0vRfsVf6jC02Jw+++h0L3rkbg776H2hjQEpJlnUjaSVL0AkD7aMvJtHcjT6J97Ah730U9CAEAoWLWk5FEWemvkOEwciQ5m34MzMNVZ4YWJgVPdd2PJcBzVL4dsyNXPXLDJ/OVq9TyLeVhayjbbR5tFLK4kp/LxRCvC+ESJRS5p3usWuDw+5k1bx1qCe1DLdbHfz04WK/ONouAzvi9CEFaAo30W/Muad9/Mq43ZK7X/ya9MwCnOUj6GXr97N9fybfvnZbwARcjmsWfLlgHZ/M/Qu7Q0UC63em88PSbXRp29Rnx2BVdbNm2yHO790uIHY1BFJKVuW8QlrJr6jSjkCwr3g+PeNuolfC7bSLuphsy0b0ujCah/c9Y53scV4fcgnXLp5FZlmxJ/VUSs5r3pp7ewwCoEdCc8a17sL3B3fg9hHQTIkK4QadoZp1cCqEEM1FuVq2EGJA+TH9P8yrBtWh+uxkC2Ar89VQr+4kJMUx4YExmCNOJGYbTAYSkmK5+NZhfjnHcdbuOERWXnGFkwVwuSUlZXZ+X7vPr+c6mYJiC9PmrMJW7mTBo0O791AOu1NzcPn4nN1IistsXssbE7m27eVO1gZIJG5c0s62gi8odWZiVmJoHTWclhGDzngnC5AYFsH/RlzFPT0GMTK5A7d37cfT/S/CqJx4anm41/mEGwzoqNSgUzHweO+hmJQQzRoN5RitEGI2MAxIFEJkAM8ABgAp5YfAlcC9QggVsALXSlmXPpKnR3hUGCldWpC2veqklE4n6O/H0eZd/7mJLgM68v3bCykrKuP8ywdyxSOXEhbp30ek1Ix8nKr3JJPV7mTf4Vy/nutkNu3JwKAo5e1qqp5bVV2EmQxeAuAul5u+3VICalegOVy2HNVnl1tBRtlfdIm9vMFtChZ2l8rDK+azJGM/zvLJTh2CT3et55YufXiy7wgAWkXFsuDS23hry0pWZx+mWVgk9/cczEUpId5DLVRDB1LKGjOzpZTvAu/6zaJ68Ogn9/LERc/jdKioDhWj2UBYpJk7X73Rb+cQQnDhVYO58KrBfjumL1olxWHQ63Ce9JgeZjLQLjkhoOeOMJt8fg91QtA+JRGDQWF3WnZFFZnZpOfa0X1oFt+42/vohQmBDknVz1ygQ69rJOWlfuK1TctYmnGgwsmC56nF7lL5cs8mhrVox3lJbQBoHRXH1PNP2Wk7pBB+SJQRQpiB5XjaIuuB76SUz9S0T4iO8etGlwEd+WT7VOa/v4hDOzPoNrgTl9w1CofdyZTr3mD1TxvQG/SMvPEC7njper+MQqWULP58KbOmzMVSYmPQ+L5MevkGYhKjT+u4g85pQ2JsJEdzi1DLwwc6nSDcbGDUwM6nbXdN9O2WgtGgYDkpEmAwKFw+shcdUhJZtGo3i//aTbjZwMQR5zCoZ5uA2tQQtI0azfaCmbhkVUcrkaREXBAkq4LD7L1bKvJnT8aqOvl2/7YKR3sWYwdGlCcAGIA/hRA/SylXV7fDGeFoAZq2asKkV06MYK1lNu7q9SiFOUW4VDdgZ+HHS9i7/gBvrZxS0YSxvky57k2WfbOq4v0vn/7Osq9XMevwh0TGRtT7uIpOx8dPX8t/PlvCsg37cUvJoJ5tmHzrSMymwMYH9YqOdyZfwcOvzcVmVxECVNXFw9cNpWvbZgBcOrQ7lw7tHlA7GpoYYysGNHmEtblvIFAQCNy4ubD5C5iU07txNjZOVXDgqMYJNxr8EDooD42Wlr81lL9qPPIZ42hPZumsPykrspQ7WQ9Ou5O07ensWLWHHkO61PvYWYdyqjjZ41hLbUx74kv+Me2eeh8bIDYqjJcevBQpJVLSoK28O7Vuyk9v38XmPUewWB2c27klURHmBjt/Q1PqzGJN7lSOlq0GdCSYO9EhehxtokZg0NX/htlY6d2kBRtyj/hcF643MKFtIxbU8eNElxBCATYAHYD3pJQ1lomesY5238aDPrMO3C43qdsOn5aj/eXT36tdt/rH9fU+7skIITjNgXcFe9JymDZnJXsO5ZLSLJZJlw+mb1ffk1iKTlftujMJp9vCgvRJ5boFnhtynm0nTncZHaIbV+zRXzw/cDRX/zIDq+qs0qPQrOgZ1rJd6E92nYraO9pEIUTlH/O08joAz2GkdAHnCiFige+FED2klD71YOAMdrStuiVjCjditziqLNfpdbSsRdfWguxC1v2yGb1Rz8BL+hARfaLkMiKm+vJLcwiO/rbvz+T+l7+tyI3NLSjl769/z/P3jmVYv0b+wzkNDhb/iuq2UrntqRsnJc4jZFk3khR+9jXk7B7fjF8uu4NPd65jY+4RFCHoHNuUie26M7BZymmH3IJO7R1tnpTylOV/UspCIcQfePRgzj5HO+qmC/nyuW9xWJ0czzZTDApNkhM4d3jNMcbv31nIJ5NnoNPrEELgdkme/vZRBoz11P+PmzSSaU986fOfduWj4/1+LafL27OXYTu5oMOh8vL03xhybrtGX9VVXwrs+3yKxkhcFDnSTuloXW4H2wpmsL/4R9xSpU3kCHol3NHo47opkbE8O2BUsM3wOwK/ZR00AZzlTjYMuAh4taZ9GoFSb/2IjI3g7VVT6HFBF3Q6HYpBYfD4fkxd9jy6GgSKU7cf5tP/m4nD5sRWasdaYsNusfP8Vf+lrNgCQERMBA+9f6fX3f3cET0Yf8/ogF5Xfdh7yHf+bWGJlSsfn05eYanP9Wc6cab26IX3E4iUkgh9sxr3lVLy29FH2V7wBWVqNlZXPnuKfmBh+p243I4a99UIEv4rWEgClgohtgLrgF+llD/VtMMZO6IFSO7Ugql/PI/qVBFCoNRi5Pbbl8uraCYcR6cTrP5xAyNv8KT7jL97NAPH9WHBR79ScqyUUbdcSNeBodkdNi46HGtukc91OfklvPTpr0x99G8+11tsDpZvOEB2fgl9uyfTo32LQJraoLSLHs3mY5+guuxUfjxxo7Ii+zlG698m0dzV5755th3k2Xbgko5K+zmxqHkcKv2DdtGhd8PVwF9ZB1uBOqnin9GO9jh1UfBy2BxIH/J/UkoctqojlaYpidz2Yugr7d962QCmfrnUK3wA4JaweushnKrLK4Swbf9RHnhlDnaHEymBbyExNoLPnr+eJnGNu0gBwKCL4JKUT/jtyGMUOVMrrZE43WX8kfkkV7SZ6zMumW/fg8T7e6JKKzm2bZqjDVUaq9bBmcb5fxuIKdy7Gsilun2W9LrdbvZvTuXg1kM0YOVxnbjswh7cNmFgDVtIL9tdbjeP/vcHbPZyJ1tOXmEZt/x7pk/B8MZIpCEJvc73BKbdVUyhw3eniQh9c3Q+ximKMBFtOPMzNhorwdI60BztSZxzYTcuuHIQ5ggTQoBO0WEKM3Lbi9eS2LJqCez2lbu5Lvlu/jH0aR4e8i9uaHMve9YfCJLl1SOE4NbLBjJmSBeUk3JyhRD06twS40mj/h37s7DYfCevF5RY2bir8Qqee1P3X1bLiIEYlSgvYXCPxOKpGpJoBI0gdVg4K0IHdUEIwePT72f0zcNY/t1fGMOMjLrpQtr3alNlu+L8Ev45dgq20hP1qrYyO5Mvep5Z6R8SHhV6epyPXD+MrXszKSixYLU5CTPpMRsN/OsO78dc1eWqfoQuITO/2Pe6Rkj76HEU5qV5CXwbdZHEGtv63Ecn9IxN/pDlWc+SZ9sBCGKMrTi/+dOnnXUgpRu7qwiDEnlWqIU1GNI/WQf1QXO0PhBCcO7wHpw7vEe12yz9aqXPWK7L5WLFnNVcfOvwQJpYLVl5xSxdvx+3283Qvu1JaRZXsS4uOpyvX72Fpev3szcth9ZJcVw0qLPPtuI9O7ZA0ekq9BYqo+gEnVvXPCvfmOgcM5HDpcvIs+1ElVYUYUKgMCxpCkJU/9AXYWjG2JQPcLhKcOPCrJy+Duu+op/YkPc+TmlBoKNL7BX0SbgHnTg7U/D8Tqiqd2n4pjCnyKsYAjztbQpzgjPam/v7Vt6YsbQ8pir56LuV3PG3wdwyfkDFNkaDnosHd+HiwTVXxhn0Cs/fO5bJb/9YZbkQ0K97Kzq1bhKAKwgspc5MdhTMJs++izhje7rHXU+MsRU6oWd0y7fJsm4g27qZMCWBNlEjaz0yNSr+mRg8XLqcNblTq4ysdxfOAaBf4v1+OYc/ybGW8v2B7RwtK2FQ8xRGpXRCX0PqZCjQaDssNFYKcoqQbjfxzeNOvbEPeg7thjnSXCV0AGAw6ek51HdKUCDJPlbCGzOWemnJfvr9aob2aU/blnWXWBzWvyMzptzIy5/+xv6MPMJNBq4efS43XTrg1DuHGAX2A/ycfjeqdCBRybPtZH/xQgY3fZyOMeMRQpAU3i+ovcC25H/qFb5wSRu7C+fQO/5OFJ33k0dD4XS7WJK+nw25R2gVGUtKZAz3Lf8Bl3Rjd7n47sBW2kb/xbdjbiRMH8LhDs3RNgwHNqfy0g1vcfRANkJAyw5JPDnrYdr2bF3tPiUFpSh6pUrctfeIHnQd1JGdq/Zit3g0FUzhJnqP7EnXgQ1f1rp8wwEE3mlIqsvF7+v2cUc9HC1Ax1ZNmf7c9adrXtBZm/smTmmptEQiUVmV8wr59j0MbPIoNtcxjlrWodeF0TJ8UINr0Zaq2dWsceNwlxKmi29Qe45T6rRzxc8zSC8txKI6CVP02FxqFZ9VpjrZX5TP/3at576egdVsrjcBmuiqDWeVo53x4nd8/vTXVZal7Ujnrl6P0SQlgVueu6ZKbHX/plReu+09Du8+AlJyzrDuPPHZAyQkxSGE4KUFT/LzJ0tY/PkfCJ2OsXeMYPQtw4JUDy6R1X2LpMThVJm7ZCsL/tyBTgguG9aDCcPOQa+E9qOev8ixbatmjWR/0QKkdLG/5Gd0FTKJKmZdPG6cJIX3p3fCnUQaTq2RcTrEmzqRZfUWJVKEGZMSE9Bz18R7W/8irfhYhU6tr6aMADaXyg+pO0LW0Qq00EHAWTFnNbOmzKl2fW56Pu888Cm2MhsT7h9LQU4Rjw5/BkvxiVr4LUu38+iwZ5i+6010Oh16g57x917M+HsvbohLqJGhfTrw9uzlXsv1eoUL+3bg4f/MZcfBrIoutmmzlvPnplSmPjqx8QuF1AKDCMcufZfGurCzt/hHwF2lBKHMlQXAwZLFZJStYkLrGYTrEwNmY5/Ee1iU8UCV8IEizPRJDO5k2Ly0ndWKgZ9MTe3KQ4FgOdqzYzgDfP2feT472VbGbrHzxbPf4Ha7WfS/pV6ddV2qm2NZBWxeuiOQptaLZglRPHTdUEwGBb1eh17RYTLouXX8APKLLOxKza7SKtzmUNm4O4Nt+zKrHOdITiGzf9nIN4s3kX2spKEvI2B0ib0CHTXFDmvK+3Gjuq3sLPjG32ZVoYm5Gxcnv0vzsL4YdVHEGdtzQfOn6RwzMaDnPRVKLW/EYYqe6zv6tyu039HyaANLQXZhrbazlNiwFFvJ2HsUh4+EfbdLkp2W42/z/MJVo3pzXq+2/L5uH263ZFi/DrROimfanJVeTRUBnKqLLXuPcE4nj37BzIXr+fC7lZ4vmoB3vlrOYzePYMKwngDkF5aRlV9Mq+ZxjU4M/Jz4WyhyHCKt9Ld67e/GSbZto5+t8sbjbN8J+HnqwtUdzuG9bX9hqyZkcJwhSW24tmOvBrKqnmihA/9w9EAWGXszSenSgqS2J3I9e4/syeLP/6i2NflxTOFGwqLMdBvcmWXfrPIpHt6ht+8k9lCgZdNYbrqkf5Vl8TERmI16L60Do14hLjqMhX/u5JvFm9idmu31PXz9i9/p2zWZD75dyfKNBzDoFZyqi6tH9+aBay5oNGEHndBzYdLzJBX2Y03uf5EIJLVX2RLoztrS2ju7D2Rl5iG25mdiqabVjVGn8NyAUSihnN4VoPLa2hDCn0rdsFvtPDX+Ze4851Feuv5NJnX/O89d8TpOh+eLcdPTVxERE46ir/6SzREmrvvn5SiKwojrzyc6MQq94UTMyRhmpMf5XejYp13Ar8efjBrU2ac0pF7RsXzjAV79bAm7fDhZACQ899EiVmw8gMPposzqwOF08d1vm5mzZEvAbfc3nWIv46q2P9Aj7tpThBKqohNGuseFvoBQIDAremaPvo4vLrqGnvHNfW4TYzSTFNEIdHiDFDo4YxztJ5NnsGnJNhxWB2VFFhw2J2t/2cTnz3jias1aN2Halv9y6T2jadM9hYGX9OHKR8eT2MKTMhOdGMVtL17H1Y9dRn5mAd++Pp+OfdrReUAHohOiSEyO55rJE3l+3uRgXma9iIkM4+3Jl9M0PpIwkwGzSU+LJtE8dvNw1mw7hM1HWOE4DtXFjgNZ2E/Kz7XZVWYu3BBo0wOCWR9HcsSQastbzUo8zcP6osOAIoyEKYkMS3qBeFNVGcwiRxoHixeTY91WpVxZSjcOVxlSBqneMwAIIThQlE+mpWrcXuBxxC8OHI2uETzdCHftXv7mlKEDIcR04FIgR0rpVZMqPM+ObwHjAAtwq5Qy8MGsSkgp+Xn6Uq+YqsPqYMFHvzLp5RsAaJKcwANv38HRA1m8cPVUNi7Zhk4Imrdpyj9nPUy3QZ3Ys/4Aj494FtXpwml3Yg43ERkXwXtrX6l3cUNDYrU5KbHYSIyNrNLUsWeHFsx/805Sj+QjhKBNi3hmLtzgs8TWi2p+P0Wl3t0JGgsJ5i6e8tqTRi+KMNMr/na6xF6Ow1WGU1oIVxKqlOK6pcqyzH9zxLIagQJIIg0tGN3iLfaXLGBbwReobjtGXSR9Eu+mU8yEhr24ADBzzyZeXP87VteJ35gAejdpyTP9L6JXYmBT3/xFKIcOPsPTD6c6xgIdy193AR+cvll1Q3Wq2K3esVTwtB0/edu/D32aA1vScNqc2K0OstJy+L/RL1CQU8Trt72HtdSGs3yUZ7PYKcgu4n9PzQ74dZwONoeT5z78hdH3vs+Vj/2PcQ9+xJI1e6psI4SgXXIibVsmIIQgLjrslG1s4qLNNI2L9FougF6dWvrzEgKKy+2gyJGG3eUpj1aEgQuaP4MiTBUhBL0II87YjlaRQwEwKhFE6Jt46R3sKJjFEctqXNKOKi2o0kqR4xC/ZNzPlmP/w+kuQ6JidxeyNvctDhb/2rAX62eklLyxZUUVJwuee1Sh3dponGytwwbBCB1IKZcDx2rYZALwhfSwGogVQjTYJy+l5Lkr/1vtpEzPC6qWw679eRPWUpvXpJjqVPnpw0Vk7D3qdQyX6mLVfP91tw0Ez330C0vW7sGhurA7VQqKLTw3bRGb9/huHQ0wrF/HGluZGw0Kf79xBJNvG4nZqK/oyKvoBGFmAw9eN9TflxEQdhfO4euD4/jp8CS+Sb2MZZlPobptJEcMYWLrWfSIv5HWkSMxK/Hk2/cwJ/Vyfjp8O0WOwxXHsKkFZFk3UerMYk/R97hk1Ru7RKVYPeSzhHbzsY8b5DoDhdPtJt9m8bkuo9R3544T+7r4/uB27lo6h8dXLmBTrvfvq0EJVUdbC1oClcVJM8qXeSGEuEsIsV4IsT4313cfq7qydflOtizd7jObwBxu4v63bq+yLP/IMVyqd5qK066yaemOaqUBjT4UrkKFgmILf2466BVHtTtUPptffbv5iDAj/3n4MszGqhEkRSfQCYHJqOevLamkNIvjw6euYVi/jrRPTuDSod358sWbaFfPst6GJL10Jevz3sMpLajSgls6SC/7k5XZUwCP8HfPuJvIsm6gVD2KxIUblXz7Hn7OuAeny8LqnP/yXdrf+P3oZH44dC0WNb9ONlhU/3zX/UG2pYQHl8+j++yp9P76LV7esBRbNZkEx9ELQazJt+xnSmT1FWtOt4vrFs3iX6sXsTh9H98d2Mb1i2cxfee607qG+nK8MiwYwt/+SO/yNSTyPYHt6Ys+DaBfv35+uZxNS7b5TMFCwN8eHkeb7lVTcroO6lStM922fKfPkbEpzMgld13kD3Or4Ha7+f6thXz73/kU55fSuX977pl6K537ta/TcXILSjHoFS9BGYAjOdWPOFxuN1Nn/oGr0k1KJwRuKZESSsrsLFq1m+UbD/D5CzfwykOh1+H3VGwr+NzHKNPB4bLl2F3FmJRoDpUuxeWu2jcMJC63nb9y/kN62XJc0lGlP5gvBAoS7/9BjLF6HY2GpNRpZ/yCz8m3leGSkjLgs13r2ZqXyeyLr/fa9rv925i5dzP7i/J8/qDNip4n+gzzeS6Hy8UbW1awLT+roqpM4inffXXTMq5o35MYU8PnYosgdQbxx4g2A6jszZKBBns+iEmMxmj2nj02h5to0d47FaVD77Y0qW4kJvEaGesNCr1H9uSaJ/w/ofHx5Bl89u+vyD9agNPuZPufu3ls2DOk7ahb94KUZnG4qtGN7dWp+maKKzencjS3CKd6wjkcd7KV31vtTj6Zu7pONoUKFjXP53Ideuwuz02o2HHEZ9txVdrIsKxCPclR+z6ekXPib0URVYVoFGGib4hIHM49sJ0Shx1XpX+w3e1ic95RtuVnVSxLLy3kwrkf8dy639hXjZNNjozh5s59KHM6yLFU7aK8MjONvt+8zbQda3yW7hp0OtZkH/ZaHnBCOUZbC+YDNwsPg4AiKWXmqXbyF8OvG4LOhzCK0AkuuHKQz30Gja+9FF5S+2a8MP//MBj9K/1WVmxh/nu/YLNUHY07bA5mvvhdnY4VZjZwy2UDMJtOPKAIITCbDDX2CtuTlo21mnY1lXG7JZv3ZNTJplChWVgvfH3NhdARYfDciI9afN9EFExIWXM11HHiTR04N+EOLmj+LLHGtuhFGAmmLoxs8Rotwvuf+gANwOa8TK8JLQCEYFfBiWrHp1Yv5pjdWq2/MeoU8qxlzNq3mX+vWcQF33/A+9v/AjyTY3cunUOJs6pDr4yUEGlsWGW044Rs6EAIMRsYBiQKITKAZ8AzTSul/BBYiCe1az+e9K7b/G9m9cQ2ieH5eZN54eqpuJwuJBKj2chz3z9BRHS4z31GXHc+C6b9ViFvWBN2a+2rh+pCVmoOeqPeKyXN7Zbs25Ra5+PddtlAWjaJ4fMf13GsuIwa6ceCAAAgAElEQVTeXZK558ohtGwai6q60Ol0XhNfLZrEEGYy+CzPPZmm8Y2z622v+DtIL/sTp9vKcT0DRZjpm3g/ijBQ6EijwOG7z5tRiSLB1IkMyypONcw5PiJuHXkhrSMv9Ocl+I1OsYmYyyUOK6MD2kR5UhellPyZmVq9EhzgOD5KrTRYfWfLSgY1a8WegtxTDgjD9AYGNg1SlV2oluBKKWssh5GegGdQn416j+jJt1mfsHvtfoRO0Ll/exSl+rSlzv07MPGBMfzwzs+oDhWhE6g+4puKXmHA2Dq1b681TVsl+hS5EQJad02u8/GEEFx8XlcuPu9ElsXeQ7nc/swsdqZmo1d0jBrUmcduHkFEmGdib+SATrw9ezk2h1pjB1+zUc8tl4Wu2LfNVUi+bTdmJY54U6cqcfZoYzLjW33GlvzpZFs3E6FvRs/4m2kZ4XnaKbDvQ4fiI7IKieYu9E28j+z0zahuOxLfo1sdBlIiLgjEpfmV45oF9kpasgadjuTIWPo3PfGd0wkdLlk7ta7jWF0qV/8yk6bhETiq0UQw6BRijCa+uOjqoJXq+mO0KoRIAb4AmuO5e0+TUr5V0z5njNaBolfofl7nWm8/6ZUbGXnjUP6avx6DUU9EbATvPzIdp82J2y0xmg2ER4dz47+vCoi9UXGRXHTTUH6ftaJKSxxjmJHr/3XFaR8/t6CUe6Z8TVn5iNypuvh19R6O5BQy7d/XAmA2Gfj46Wt57qNf2J3qEZ3u2rYp0ZFhrNt+GEXxjILvv/p8zj83NMuON+d/wvaCGeiEASldRBiSGNXyTSL0J1rtRBlacn7zf/vcP9LQ0ufoTYeBOGMHYoytuazVl+wsmE2ObTuKMJBn240bJyBRhAmzEkf3uNAXR3dJyWtDxvHBttVsP5aFTghGJnfg5UFjK25OeTYLUQYjx+zVF6PoELh9fGaqdHO0zLfim1Gn8GS/4dzYqU9w2934Z0SrAo9KKTcKIaKADUKIX6WUO6vb4YxxtPWhbY9WtO3RquJ9p37tmPvmAjIPZnPuiB5MfGAsMYmBq99+6L1JRMVFMP/9xditdpI7JfHA23fUOevAF3N+24JDrToqcaou9qTlsPdQbkXPr1bN4/j0mesoLQ+jRIZ7YmfFZTYKi60kNYmuUtSwZO1eps1ZRVZeMa2T4rn/mgsYWEN3ikByuHQFOwpmVckIKHYcZunRyVzaanqtjpFo6kq0sRWF9gO4K41YdUJPp1iPPGGkoTkDmv69Yl2+bTc7C7+hTM2mZfggOsf8DaPiXdQRKpQ47Dy8Yj4rM9PQ6xQkkn/2Hc4tXfpW0Y+VUnLDr7MprMHJAj6d7MlUdsbhegMjkttzS+e+wRUhkv4pry2fg8os/7tECLELT0prtY5W1PTIGEj69esn168P7SKAhkJKiepUT3vCTUrJr6v3MOvnDRw8kl9Ff/Y4EWFG/nXHaEYO7OTjCDXz04odvPbZkioqYCajntcemRAUZ/tL+n1k2zZ7LVeEiQmtZxBlqF3lmt1VzKrsV8go+xOJJNbYliHNniTBXHMDy8bCHb9/x4qjqSdiq0CYYuD9YRMZ3vLETX1z3lGuXzy7WoWukxFUP0A0K3oubNEOt3RzRfuejG7V6bS0EIQQG6SUp9XQLTIhRfYY+/dTbwismfnoIaByysq08vTUk+1qAywHekgpq+3KelaPaEMFIYRfsho+nvsXM39ej60GgXPV5aZ9St27BEgpef/rP72kFu0OlXe/Xs7AnjdhsTlYsnYv2XkldGnXjMHntAloLC7fvsfncoEOh6uE2opzmZRohrd4CZfbgRsnBl2EH60MLnnWMi8nC2B1Oflw++oqjjazrLhOzrBtVDyHSwtRfYjnpETG8NHwy+tveKCo/cAy71SOXQgRCcwBHqnJyYLmaM8YSspsfLlgnc+iheOYDHr6dk2hTYu6N/mzO1QKSnyXYR7KLCD1SD53v/g1DqcLq91JuNlASvM4PvzX1YQHoKqu2JFebQGBWzqJNdU9/KLojCjU3Van20Jqya8U2A8Qb+pEm6iRGHS+K6kammN2Cwad4uVoAbJOyn/tmZCE01W7STAdgsSwcIyKwoHifJzuE842TNHzj3NDc3LQX6lbQggDHic7U0o591Tba462EtYyG/Pf+4Vl3/5FeFQYE+4fQ7terVmzYCMGo57zLx9IXLPYYJvpk/0ZedVWhwkBURFmJg7ryZ2X169xnsmoJ8JspMRHSlyz+Cie/mAhxaW2ikdJi81J6pF8Pp+/lnuvPr9e56yJEucRFGFCld7OP9LQsloJRH9T6sxkQfokVLcNVVrRizA25X/EJSmfEmFo2iA21ETrqDh8DVL1Qsd5zauGe5IjY5jQrjs/pu7ynW9bCTeSdTkZSE60utELHTFGM5P7DGNs6xAMu/ipGKFcsfBTYJeUcmpt9tEcbTkOm4OHBj/J0QPZOMpn6rev3I10S3SKDrfbzTsPfsrwa4fw4LuTiIwNrcfLJrGRqKr3I5wQMLx/R15+sO7ls3aHitPlIjLMhBCCWy8bwMdz/6oSPjAb9dx4ST9e+/x3r++ww+ni51W7AuJoY41tfaZbCfQV6luVKbSnsiHvfXJsWzErcfSIu5EO0Zec9uTM6pzXsLuKkOU5uqq0orrsrM19g+EtXj6tY/sDk6Lnn32G8+L6JRXda/VCR4TByIPnnOe1/SuDx9IrIYnPdm+g1GGjR2ISvRNbMG3HGooctir/4+N/Hy9M6BnfnLnjbgrprht+0podAtwEbBNCHJ8keFJKubC6HTRHW87vs1eSlZpT4WQBXOWjQ3el8tbfZ/3JpqXb+WjT68Q1DV4L6JNJbhZLt/bN2bb/aBWHazLouXFc3eYQikqtTPnkV1ZuPogE2raI56k7L+aGcf2w2lW++mUDVruTmMgw7rvGk/r12ue/+zxWoH5yEYZmtI4czqHSPyopaQkMOjNdY6+ssm2xI50F6XeWFxVIHO4SVuW8xOrc12kVcQF9E++rVytxKSVHLWsrnOwJ3OVFDqHBDZ17kxIVywfb/yLLUsp5zVtxf8/zaBERzZHSItbnHiHRHM6gZq1QdDpu6NybGzpXzR+PNpp4Yd2SGrvhbjmWiSrdGILYsfdU+Cnr4E/q+NXWHG05a3/e6FucxgdFOUXMmjLHSxks2Pznkct4+v2FrN+ZjqLoMBkUJt86ku7ta+9EpJQ88MocDmbkVYiC70/P454pX5OUEE1GThFut5vOrZvwwv2XkNzMU1HUPjmBPYdyqsw1GA0K4y7o5tdrrMyQZv8i2tiaPYVzcLotNA/vQ//EhwivlEMLsPXYZ+V6BVXH3G7pIK10KZmW9UxsPQuzvj7C7jp8dtANUgVSdQxt0ZahLU70upNS8tzaX5m5d3NFileU0cTs0dfRNto7hn+w+NgpW44rQqAL2K3VD0jqMhnmVzRHW05iy3gUvQ6Xj8fvk3G7JKvmrQs5RxsdYebNxy+nsMRKcZmNlk1j6jzrv/NgFulZBV6dF2x2ldSjJ2SJd6XlctcL3zDvzUkY9Aov3HcJd77wFQ6nis2hYjLqadsigVsuDVxFmU7o6RV/K73ib61xu1zbDqpvJ+5GlVb2FH1Pr4S6/T+FECSYOpNn924/L3HjcttRdMGp6T8VPx/ew1f7t+Jwuyomyiyqg0m/f8dvE+70evzvFt+McL2h2tQvvRDEGcM45+s3iTWaubPbAG7u0jfk2tuEcoeFs4JL7hqF3lj7+05EjG8dhVAgNiqMVs3j6pValZFTVKsYmyxX9Vq+0aMT0CopjnlvTmLyrRdx9xXn8epDl/HpM9dhNjXMpFRNRBtqLml2SQc5tm31Ora+GkeqCBM5tu31OmZD8MXujVhV744JR8tKOFDkrbd7aesuxBjN6E/qNiEAk07BLSHXVkaZ08GRsmJe3bSMKeuXBPAK6kkjVu8KKGXFFhZ+/Bv/e2o2q+atw6XWrQa7trTumszkLx4iIjYcc2TNOpkGk4GJD44LiB3BplOrJrjctQtkOZwqmbkn0gfNRgOd2jTFZDSQnV9MaTXthRqanvG3eMkXVkagJ9bYpl7H1leTxiXwjLhDlTKn79Q4RScoU51YnA5m7tnEwyvm89/NyymwW5k37hbGtu5MmKIn0mBkTKtOPD9gFENbtEWIqv7JqjqZsXczRfZTS0w2FI1d+DtgpO1I5+9D/43qULGV2TFHmmnRvhlvLH+B8Cj/5ylecPlABo/vy4HNaTx/1X/JOexby3T0zRcy9o4RFe+PZRWwYNpvHNx6iC4DOjB20kiiG6naVduWCQzs0Zo12w9VVJYJ4Tu0ZdArdGvXDPCMcP/75VLmL9uO2+1GryhMnfEH/3nkMgb0aNiqMZurkLSSJdhdxSSF96OJuQcXNHuOtbn/xeLK4+QhiyL0dDlpAq22dIweT5Zlg5dmraIz0sTs1cs0ZLi0TRf2FeV5K3kJQVJEFKN//JRjNgsW1YlRp/DpznV8cdE1vDPUW5d55r7NPiURjTqFg8XH6N2kek3kBkXKRi38HTBeufFtygrLKiapbKU20ncfZfbLp8wPrjd6g57O/Tvw+tJnCTtpZKtTdDz3wxM88tHdFY/XqdsPc1uXh/nqle/5c+4aPv3nLK5IvJ2rW0zixw8X1aiKFWhWb03jhie/4ILb3+LKx6az+K/dtdrvpQcv5Zbx/WkSF0FUhInRg7vQOim2iuaB0aDQqXUTenfxPJav2XaIH5fv8KSEqW6sdidWu5P/e/tHn6XAgSLTsp45qZezPu9dNh/7hF+PPMwfmf8iJXIIV7b9gb+1+YZWEcPRYUCgJ8bYhlEt36p1ue7JpERcQPvocSjChCJM6EU4Bl0EI1q8hi6EZ99v6tyH1lFxhOuPN6bUEaboef28S/hg22qyLSUV8ViH24VFdfKPP3/y+X1uH53gcxLM7lJpGRE4rZB6EaTQQciOaAtyiji8O8NrJOW0O1kyYwV3vHRDQM+f1LYZ3+Z8wgePfMa6RVtIaBHHPa/fTLfBVRXC3rjrIyzFJ0Q4jn8RC7KKmPbYl1iKrVzzxMSA2uqLJWv28NT7C3GX38HTswuZ8sli7E6V8UNrHmkZ9Ap3TBzMHRNPFDeUWu1M/2E1i//ag04nuPSC7tw8vn/FDeenFTuwVaNru3F3OoPPaetznT9xS5U/Mp+qMrpUpY2jljWklS6hXdRoog0tGd5iCi63HZd0nrYYjBCCQU0fo2vsVWRa1mNSokmOOD9kKsOqI9xgZN4lt/Bj6k7+OHKQpPAorut0Lu1jEnhm7a9VKr2Ok20tJctSQtJJzvPeHoNYkrHfa3QMklLVQfDLNk4QrMmwkHW0Op2o9s4iaujc6i+cDieTR73Agc1p2Mrs5GXk88So57n3jVu55M5RgKc77u61+6o9hs1iZ9ZLc7ni75eiNzTcR11UYuXfH/xc4WQr7HGovP/1n1x6Qfc6J5VHhpl46LoLeeg636LWNcV1XQ30uJZr2+6zZ5cqrazOeR29MNGqXJRb0ZlQ8F9GQIyxdcj0BjsVGaVFvLrxD5YfTSVcb+DmLn24s9vACvlCs973d9UtJSbFe12PhOZc1eEcZuzZWOUn63C7ufePuSy6bFIgLqPuSEALHVQlJjGadr1aezlVo9nA6FuHBfz8S2evrHCy4ClasFscfPDIZ5QVe8o+hU6g6Gt+PHQ5XRTm1qg34Xe+WrTRZw8xgIISS0Ae5cec15UwHxkGbrekX9eGUtOv/ubsdJeyIus5Nud/UqsjuaWLHOtWjpatLe/OcGZwzGZh/ILPWJC2myKHjUxLCW9tWck/Vv5Usc2NnfoQdpJDVYSgd2IL4s2+s23WZqf7/OjTSgpP2ZK8QdGyDrz558yHiUmMJizKjKJXMEea6dCnHddODvyj+LJv//JZwKAYFLb/6Yl16nQ6hl1zXo1pYULREZPYsBNjK7dU3wonzGTAVIc0ttpyQe/2nN+7HWaTAQEY9DpMBj3P3D2mwVK8mpi71xgXVaWN7QUzsLlq/uEfs+/lu9SJ/HrkH/yR+S++PngJB4p/8be5QWHGnk1YVEcVTVmbS2XR4T2klxQCcEuXvoxI7oBZ0ROuNxChN5ISGcs7Qy+r9rjOaooZdDWsCwZa1oEPWnZIYmba+6yat47sQ3l07t+eXsPq/thbHyJiwn3PtkuqTJI98M4dHNmXxYHNqV79v0zhJq78x3i/N3Y8FYmxEfgWEISrRvUOyOen0wleuG8cW/cdZdWWVCLDTIw+rwvNGjD7Qif0DEt6mSVHH/NZCebZxkC+bTctI3w3rXRJJ4szHsburuqM/8p5lQRTZ2JNgY81B5L1uRnYfSh0GXQKuwtzSYmKRa/T8d6FE9lfmMfW/CySIqIY2KxVjcUHE9t2573tnjY5lYk3h1f0IwsFgpV1ENKOFsBoNjLsmiENft5L7x7FqnnrvBo4miJMdB9yYkIsIjqct1dNYe+GAyz9aiXLvl5FbkY+UXERXDN5Ilc/7v825afi+jF92bAz3Us7tnlCFPdeFbjPUghBr04t6dWpfjP4p4OUkr1FP7D12Ge4pANFGCtpIFTezkWYvpp280CmZS1uH51v3VJlb/F8BjR52K92NzQdYxL5K/MQzpM0ZF1uNymRVbU7OsQm0iG2dtrFk7r1Z1H6XlKLj2FRnZgUPXoheHvohNARmQlQWKA2hLyjDRbnDO3GDU9dwZfPfYuh/FHbYDbw8s//8tn4sVPf9nTq2567X7sZl8tVY3PIQNOveysevHYo73y9AkUnUFU3HVsn8tojE0PnS+9ndhTMZPOx6bjKMw58NRcUKEQZkokzVq9Va3eV+OwhJnFhdxX6z+AgcUuXvszatxlnpVJzg06ha3wzusR55wdIKfkxbRfTdqyl0G5laMu2PHTOEJqHV31SCTcYmTfuFn5N38ua7HRaRERzebseJIaFjsqdp2AhOJ5Wa2VzCgpyiti6bCcRMeGcO7x7g2YPnC42u5N96bnERYWTHKI6uv7ALVW+OjAGpw9tWoEOvTDjxkWcsQPDW7xMuL76UVqZM4fvD13tJSquF2EMafYv2kSNqGbPxsOGnAz+76+fSSspAASjUzry8uAxRBu9KyKnblrOxzvXVejT6oWOaKOJRZdNokkDOlF/tLKJjk6W/fo/UKttl/7+z9M+X2Vq5TWEEGOAtwAF+ERK+cpJ64cB84DjszBzpZTP+8vIYBLXNIYLr6qfWHawMZsM9OwQIlU5AcTmKsTts2E4GHVRjGr5JiYlulZSiBGGpnSLvZ5dhV9V5OPqhZl4UyefOreNkb5Nk/l1wp0UOWyYdApmve85hCK7jY92rq0Sd1Wlm1Kng//tWscTfYY1kMX+I1gj2lM6WiGEArwHjAIygHVCiPk+WuuukFJeGgAbNTRqxKTEIKpJoIk2ppBgrn0beoA+iXfRPOxc9hbNwynLaBN5Ee2iLw5p7YL6EONjBFuZ3YU5GHWK1wSXw+1iZeahQJoWGEI8RjsA2C+lPAgghPgKmEANrXU1aubw7iN8/eoP7N+cSodz23LN5Im06tLwE0hnCoow0D3uerYXzKyI0XqWm+idcFe9jtkiYgAtIgIn8dgYaB4e5TM1SwApUaEjel97QlvroCWQXul9RvmykxkshNgihPhZCNHd14GEEHcJIdYLIdbn5ubWw9zGz+61+7i//2R+m7Gcg1sO8duM5dzff3KNFWZnGoUlVr8re/WKv51e8bdj1HnKQ6MMLRna/DmSwv0WZjvraB0VR6/EJAwnyW2aFT13dvOdHhfySFm7l5+pzYjW1zT1yZZsBFpLKUuFEOOAH4COXjt5+qJPA89kWB1tPSN496HpVQoh3C43tjI77z40nXdXB7/HVCDZlZrNcx/9TEZ2IVJCn67JPHP3GBJjT09vADypZT3jb6Rn/I24pSukBV0aE9OGXcHDK+azKusQik6HWdEzZeDF9Eqse+ufoCP91jOsztTG0WYAlWsok4GjlTeo3NNcSrlQCPG+ECJRSulbZ/AsZu/6A3VafqaQX1jGfS99g6VSUceGnencO+Vbvn71Vo+2hZ/QnKz/iDGZ+eyiqzlms1DssJMSWfeuHSFFkCbDavOJrQM6CiHaCiGMwLXA/MobCCGal7fgRQgxoPy43jLtGoRH+1Z1qm65P1FdbtZsO8Rva/aQV1ga8PNVZv6ybV7tcVxuSV5hKRt3p1ezl0aoEG8Op010/bp2hBShKpMopVSFEA8Ai/Ckd02XUu4QQtxTvv5D4ErgXiGECliBa2UwhVhDmAn3jWHOGz9hr9Rt1xRmZMJ9YwJ63oNH8nng5W+x2lVAorrc3DiuP8P6dWD6vNUcSM+jY+sm3D5hEB1bNTnl8erK4cxCHE7viRUpITOvYUV3NM5eRC27h/ibWuWrlPcrX3jSsg8r/f0u8K5/TTszufnZq8k7eoyls1diNBtw2JxceM153Pzs1QE7p5SSv782l/yiqgn9Mxau48sF61BdLqT09AtbtTmVd//vSnp29G/+7TmdWvD7+n1emrVSSrq0aebXc2lo+ERSfY/OOiKEmA5cCuRIKU/ZSqORPwc0PhS9wuPT72fmoQ+YsuBJZh76gMen339KucXTYVdqNsVl3r2bHE4XTtVVEbaSUmJzqEydsdTvNow5ryuxkWb0yomvnMmop2+3lICMoDU0TkYgEbJ2r1rwGVDrx9AzKwO7ERHXNIa4pg2Ti2ixOeqkcbAnLafO53C53ew8kIXqctOjQ1KVtjcAYWYDnz1/A9PmrOKP9fsxGhQmDu/JjZf0r/O5NDTqjZ8imlLK5UKINrXdXnO0ZwE92id5dVuoiaiImiuGTmb7/kwee2MedofqyQUU8MJ94xhybrsq28VFhzP5touYfNtFdTq+hobfqL2jTRRCVBZjmVaenlovtNDBWYDZZGDyrSMxGfUVaVRhJgOJsRFeIuBmo57rx/at9bGtNicP/WcOBcUWLDYHZTYHZVYHT77zEznHSvx6HRoap8XxGG1tXpAnpexX6VVvJwvaiPasYez53ejYugnfL91KQZGV83u3Y8SAjrw9azk/Lt+OQa/gVF1MGN6Tm+rwOL9sw36fnVFdbskvK3dx8/izu4xVI7QI6awDjTODDilNePzmkVWWPXHrSO69agiZeSW0aBJNZHjdGhYWl9m88mMBnKqLY8WNq9eWVT3G4dJlqNJOcsTgRtNsUaO2BKa8tjZojlaDqAhzneOyx+nXLcXnRFuYycCgcxqPozpUsowV2c8CAildbMr/iK6xV9M38d5gm6bhLyR+c7RCiNnAMDyx3AzgGSnlp9Vtr8VoNU6LdsmJjBncpUoDxjCTgV6dWjCge+NwtA5XKSuyn8Ul7bikDTdOXNLOrsJvybFuC7Z5Gv6k9jHaGpFSXielTJJSGqSUyTU5WdBGtBp+4J93jGJwr7bM+2MbTtXFuPO7cfF5Xf2qXxBIjlhWo0Pxkg53STsHixfRNKxnUOzS8D8hK/ytoXEqhBAM79+R4f29BNsaBbJK8+2T14VOq2wNP6A5Wg2N4NAyfJBPh6oXJtpGjQqCRRoBQUrwMXHbEGgxWo2zHpMSzeCmk1GEER0GQKAIM+2jx9IsrHewzdPwJyEs/K2hccbTPnoMzcJ6kVb6O6rbSkrE+SSYuwTbLA1/o4UONDSCS6QhiR5xNwTbDI1AIYEg9QzTHK2GhsZZggSpVYZpaGhoBA5J0CbDNEeroaFx9hDCPcPOKgpzi9i8dDuZB7ODbYqGhoa/0bIOgouUkvcf+R8LPv4No8mA06HSY0hnnpnzOOFRgW+cqKGhEWiCJyqjjWjL+fHDxfz86e84bU7Kiiw4rA62rdjF1Ls+PPXOGhqNDNXtZtHhvTy1ehFvbF5BRmlRsE0KPBJwu2v38jPaiLacuW/8hN1ir7LMaVdZ9cNarGU2wuqpbqWhEWrYXSo3LJ7NzoIcLKoTg07hox1reP/CiYxI7hBs8wKLNqINLqWFZdWus5XZq12nodHY+Hb/VnYcy8aiejoSO90ubC6VR/78Eaf7TNZ2KC/Brc3Lz9TK0Qohxggh9ggh9gsh/s/HeiGEeLt8/VYhRB+/Wxpgeo/s6VNtKq55LLFNooNgkYZGYPj+4A6sLtVruVtKtuZnBcGiBkKClO5avfzNKR2tEEIB3gPGAt2A64QQ3U7abCzQsfx1F/CBn+0MOLe/dD3hMeHoy3to6RQdpnATf//onjp1kNXQCHXMiu+IoVtKTLrAtb0PCdyydi8/U5sY7QBgv5TyIIAQ4itgArCz0jYTgC+kp3nUaiFErBAiSUqZ6XeLA0RS22Z8vG0qc99cwPY/d5PcKYmrHh1P256NQ7xaQ6O2XN+pN5vyjlaEDo4Tawyje3yzIFnVQISw1kFLIL3S+wxgYC22aQk0GkcLkNginrv+c1OwzdDQCCjjWnfmz8w05h7cjg5QhA69TscnI644s5/epAxIRkFtqI2j9fXJn3xbqM02CCHuwhNaoFWrVrU4tYaGhr8RQvDy4DFM6taftdnpxJnDGN6yPaZqQgpnFCE8os0AUiq9TwaO1mMbynujTwPo169fcK5YQ0MDgPYxCbSPSQi2GQ2IRLqCk1VRm6yDdUBHIURbIYQRuBaYf9I284Gby7MPBgFFjSk+q6GhcRZwXCYxFCfDpJSqEOIBYBGgANOllDuEEPeUr/8QWAiMA/YDFuA2v1uqoaGhcbqEskyilHIhHmdaedmHlf6WwP3+NU1DQ0PDf0hA+mm0KoQYA7yFZ/D5iZTylZq21yrDNDQ0zg5kufB3bV41UMvagiqcBdOMGhoaGh78NBlWm9qCKgTN0W7YsCFPCHGo/G0ikBcsW0IM7bPwoH0OJ9A+CzjtyqESChb9Jr9LrOXmZiHE+krvp5VnTUHtaguqEDRHK6VscvxvIcR6KWW/YNkSSmifhQftcziB9ln4BynlGD8dqlZ1A7zfMgYAAAJPSURBVJXRYrQaGhoadaNWdQOV0RythoaGRt2oTW1BFUJlMmzaqTc5a9A+Cw/a53AC7bMIIaqrLahpHyGDVPuroaGhcbaghQ40NDQ0AozmaDU0NDQCTFAd7ala5JwtCCFShBBLhRC7hBA7hBAPB9umYCOEUIQQm4QQPwXblmBSLqL/nRBid/n3Y3CwbdKoO0GL0ZaXse0FRuFJl1gHXCelrLa64kxFCJEEJEkpNwohooANwMSz8bM4jhDiH0A/IFpKeWmw7QkWQojPgRVSyk/KZ7jDpZSFwbZLo24Ec0RbUcYmpXQAx8vYzjqklJlSyo3lf5cAu/BUn5yVCCGSgUuAT4JtSzARQkQDQ4FPAaSUDs3JNk6C6Wira39zViOEaAP0BtYE15Kg8ibwBBAcTbvQoR2QC/yvPIzyiRAiIthGadSdYDraOpexnekIISKBOcAjUsriYNsTDIQQlwI5UsoNwbYlBNADfYAPpJS9gTLgrJ3LaMwE09HWuYztTEYIYcDjZGdKKecG254gMgS4TAiRhiecNEIIMSO4JgWNDCBDSnn86eY7PI5Xo5ERTEdb5zK2MxXhaT36KbBLSjk12PYEEynlP6WUyVLKNni+E79LKW8MsllBQUqZBaQLITqXLxpJDVJ8GqFLMNW76lzGdgYzBLgJ2CaE2Fy+7MnyzhYaZzcPAjPLByMH0dpENUq0ElwNDQ2NAKNVhmloaGgEGM3RamhoaAQYzdFqaGhoBBjN0WpoaGgEGM3RamhoaAQYzdFqaGhoBBjN0WpoaGgEmP8HFpvXoMocyMkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data['labels'] = bisect(data, 8, 1)\n",
    "plt.scatter(data['x1'],data['x2'],c=data['labels'])\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.** For this question, you will apply the distance-based anomaly detection algorithm to detect anomalies in a given dataset.  \n",
    "\n",
    "**(a)** Download the lymphography dataset (lympho.csv) from the D2L page. The original data was obtained from http://odds.cs.stonybrook.edu/lympho/.\n",
    "\n",
    "**(b)** Use pandas to load the file lympho.csv into a DataFrame object named *data*. Display the first 10 rows of the data frame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    0    1    2    3    4    5    6    7    8    9    10   11   12   13   14  \\\n",
      "0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0   \n",
      "1  1.0  1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0  2.0  2.0  1.0  1.0  2.0  1.0   \n",
      "2  3.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  3.0  1.0  1.0  2.0  2.0  8.0  1.0   \n",
      "3  3.0  1.0  1.0  1.0  2.0  2.0  2.0  1.0  3.0  1.0  1.0  4.0  2.0  5.0  3.0   \n",
      "4  3.0  2.0  2.0  2.0  2.0  2.0  2.0  1.0  2.0  2.0  2.0  4.0  2.0  4.0  3.0   \n",
      "5  3.0  1.0  1.0  1.0  2.0  2.0  2.0  1.0  3.0  1.0  1.0  2.0  1.0  5.0  3.0   \n",
      "6  3.0  2.0  1.0  1.0  2.0  2.0  1.0  2.0  1.0  3.0  3.0  2.0  3.0  4.0  2.0   \n",
      "7  3.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  2.0  4.0  3.0  5.0  1.0   \n",
      "8  2.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  1.0  3.0  3.0  3.0  3.0  6.0  3.0   \n",
      "9  2.0  2.0  1.0  1.0  1.0  1.0  1.0  2.0  1.0  2.0  3.0  2.0  3.0  8.0  2.0   \n",
      "\n",
      "    15   16   17  18  \n",
      "0  1.0  1.0  1.0   1  \n",
      "1  1.0  1.0  2.0   1  \n",
      "2  2.0  2.0  4.0   1  \n",
      "3  1.0  2.0  4.0   1  \n",
      "4  2.0  2.0  7.0   1  \n",
      "5  1.0  1.0  7.0   1  \n",
      "6  2.0  2.0  2.0   0  \n",
      "7  2.0  2.0  1.0   0  \n",
      "8  1.0  2.0  4.0   0  \n",
      "9  1.0  1.0  1.0   0  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('lympho.csv', header=None)\n",
    "print(data[:10]) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** Note that the last column in each row corresponds to the actual class (1: anomaly and 0: normal). Extract the last column into a pandas series object named *label* and count the number of normal and anomaly data points in the data. After extracting the label, drop the last column from *data*. DIsplay the size of the remaining dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of data = (148, 19)\n",
      "Class distribution:\n",
      "0  142\n",
      "1  6\n",
      "Name:  18 , dtype:  int64\n"
     ]
    }
   ],
   "source": [
    "...\n",
    "\n",
    "print('Size of data =', data.shape)\n",
    "row, col = data.shape\n",
    "label = data.loc[:,col-1]\n",
    "\n",
    "print('Class distribution:')\n",
    "print('0 ', len(label[label == 0]))\n",
    "print('1 ', len(label[label == 1]))\n",
    "print('Name: ', label.name, ', dtype: ' ,label.dtype)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(d)** Apply the nearest-neighbor approach to detect outliers in the data. Specifically, for each data point, find its top-10 nearest neighbors (excluding the data point itself). Use the distance to the 10-th nearest neighbor as outlier score. Store the outlier score and predicted class of each data point in a dataframe object named *result*, with two columns named 'k-distance' and 'predicted', respectively. Display the result dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     predict  distance\n",
      "0          0  4.123106\n",
      "1          0  3.605551\n",
      "2          1  4.582576\n",
      "3          1  4.358899\n",
      "4          0  3.872983\n",
      "..       ...       ...\n",
      "143        0  3.605551\n",
      "144        0  3.316625\n",
      "145        0  2.645751\n",
      "146        0  2.828427\n",
      "147        0  2.449490\n",
      "\n",
      "[148 rows x 2 columns]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# # Get the first n-1 column of data, which is the training sample train and the return type is dataframe\n",
    "train = data.loc[ : , 0 : col - 2]\n",
    "# # The value of K in k-nearest neighbor algorithm\n",
    "k = 10\n",
    "# # n_neighbors indicates the number of neighbors to be searched. The default value of algorithm is auto, and P is the parameter in Minkowski (P = 2 takes Euclidean distance)\n",
    "neighbor = NearestNeighbors(n_neighbors=k, algorithm='auto', p=2)\n",
    "# # neighbor.fit (x) takes x as the basic sample point to calculate the neighbor\n",
    "neighbor.fit(train, label)\n",
    "# # kneighbors([X,n_neighbors,return_distance]) Find K neighbors in each X\n",
    "# # By default, the row number array array of K adjacent points is returned\n",
    "# # return_distance=True, Returns the distance between fit (x) and the nearest K adjacent points of each point of the current x\n",
    "distanceArray, indexArray = neighbor.kneighbors(X=train, n_neighbors=k, return_distance=True)\n",
    "# # Distance of the k-th adjacent point\n",
    "distance = distanceArray[:, k-1]\n",
    "# # argsort(-array) Sort array from large to small, and return the sorted index\n",
    "indexSort = np.argsort(-distance)\n",
    "# # Number of exception Tags\n",
    "n = 5 \n",
    "# # Take the n-th distance as the threshold.\n",
    "threshold = distance[indexSort[n - 1]]\n",
    "# # Take the first n large distance sample points as the outliers\n",
    "# # (if the threshold value is greater than or equal to, the outliers will be labeled as 1)\n",
    "# # (if the threshold value is less than, the regular values will be labeled as 0)\n",
    "predict = []\n",
    "for i in range(len(distance)):\n",
    "    if distance[i] >= threshold:\n",
    "        predict.append(1)\n",
    "    else:\n",
    "        predict.append(0)\n",
    "result = pd.DataFrame({'predict': predict, 'distance': distance})\n",
    "print(result)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(e)** Show the confusion matrix by comparing the true class of the data points against the predicted class obtained in step (d). Display the accuracy of the method as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy_score:  0.9662162162162162\n",
      "confusion_matrix: \n",
      "[[140   2]\n",
      " [  3   3]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "print('accuracy_score: ', accuracy_score(label, predict))\n",
    "print('confusion_matrix: ')\n",
    "print(confusion_matrix(label, predict))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}