Mobility Analysis and Clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Network Analytics Project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Project: CIS6930 Network Analytics & Machine Learning\n",
      "Group Members: ['Pankaj Chand', 'Keith Levengood', 'Nicholas Kroeger', 'Varun Puri']\n",
      "\n",
      "Library Versions\n",
      "### Python version: 3.7.6 (default, Jan  8 2020, 13:42:34) \n",
      "[Clang 4.0.1 (tags/RELEASE_401/final)]\n",
      "### Pandas Version: 1.2.2\n",
      "### MatPlotLib version:3.3.4\n",
      "------------\n"
     ]
    }
   ],
   "source": [
    "group_members=['Pankaj Chand', 'Keith Levengood', 'Nicholas Kroeger', 'Varun Puri']\n",
    "print('Project: CIS6930 Network Analytics & Machine Learning')\n",
    "print('Group Members:', group_members)\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from os import listdir\n",
    "from datetime import datetime\n",
    "import pickle\n",
    "\n",
    "print('\\nLibrary Versions')\n",
    "print('### Python version: ' + __import__('sys').version)\n",
    "print('### Pandas Version: ' + pd.__version__)\n",
    "print('### MatPlotLib version:' + matplotlib.__version__)\n",
    "print('------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "# Helper function to get the prefix of each APNAME \n",
    "def convertAPNameToPrefix(df):\n",
    "    prefixes = []\n",
    "    for i in range(len(df)):\n",
    "        tempPrefix = [\"\".join(x) for _, x in itertools.groupby(df.iloc[i], key=str.isdigit)][0]\n",
    "        prefixes.append(tempPrefix)\n",
    "    return prefixes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadPreprocessedData = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if loadPreprocessedData == False:\n",
    "    #load in csv, append APNAME to end of dataframe, save it out as .pkl\n",
    "    data = pd.read_csv('DHCP_April_2012_ANON_MAC_devType.csv')\n",
    "    print(\"loaded first dataframe\")\n",
    "    \n",
    "    print(\"Appending prefix name to table... will take a few minutes\")\n",
    "    \n",
    "    prefixes_col   = convertAPNameToPrefix(data['APNAME'])\n",
    "    data['prefix'] = prefixes_col\n",
    "    \n",
    "    print(\"Saving out new dataframe to .pkl file\")\n",
    "    data.to_pickle(\"./DHCP_April_2012_ANON_MAC_devType_wAPNAME.pkl\")\n",
    "else:\n",
    "    #load in preprecossed data\n",
    "    data = pd.read_pickle(\"./DHCP_April_2012_ANON_MAC_devType_wAPNAME.pkl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Look at the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(           userIP  userMAC                APNAME              APMAC  \\\n",
       " 0   10.131.172.58    26257  fli113-win-lap3502-1  44:e4:d9:00:ed:e0   \n",
       " 1        10.0.1.3    22699  hrb112-win-lap1231-1  00:1a:a2:09:93:90   \n",
       " 2  10.132.190.205    30264  soc1402-temp-lap3502  58:35:d9:d5:cd:c0   \n",
       " 3  10.130.146.252    72815  soc1402-temp-lap3502  58:35:d9:d5:cd:c0   \n",
       " 4   10.130.248.68    69249               unknown  00:00:00:00:00:00   \n",
       " \n",
       "     startTime     endTime deviceType   prefix  \n",
       " 0  1333238682  1333238682          f      fli  \n",
       " 1  1333238682  1333238687    unknown      hrb  \n",
       " 2  1333238685  1333238687          f      soc  \n",
       " 3  1333238688  1333238688          f      soc  \n",
       " 4  1333238690  1333238690          f  unknown  ,\n",
       " (22821939, 8))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head(), data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2012-04-01 00:04:41\n",
      "2012-05-01 00:05:20\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEFCAYAAAAWrxseAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbX0lEQVR4nO3de7Ac5Xnn8e8PZC5BRgLLewojlcWutWQJTih0Csg6dh0ZLwjHZZH4suAEhBdbu2Vw2PW6CnmrXJIvVHAqicuUbVxaS7HIxTJLUovKgBUFOOV11YJBjgzCmHDMJUaWIUGysHwjeJ/9Y94jRqPpufZ09/T8PlVTZ+btt7vfZ/rydL/d00cRgZmZWTvHlN0AMzOrLicJMzPL5CRhZmaZnCTMzCyTk4SZmWVaUHYD8rZkyZJYvnx52c3o209+8hNOOumkspsxcpMSJzjWuqprrLt27frniHh1a3ntksTy5ct58MEHy25G32ZnZ5mZmSm7GSM3KXGCY62rusYq6el25e5uMjOzTE4SZmaWyUnCzMwyOUmYmVkmJwkzM8vkJGFmZpmcJMzMLJOThJmZZXKSMDOzTE4SZmaWyUnCzMwyOUmYmVkmJwkbO8vX31F2E8wmRtckIWmLpOck7WkqO1XSTkmPp7+npHJJuknSnKSHJJ3bNM7aVP9xSWubyldKejiNc5MkdZqH2VMnvKfsJphNjF7OJL4ErG4pWw/cHRErgLvTZ4BLgBXptQ64GRo7fGADcD5wHrChaad/M/D+pvFWd5mHmZkVpGuSiIivA/tbitcAW9P7rcClTeW3RMN9wGJJpwEXAzsjYn9EHAB2AqvTsJMj4r6ICOCWlmm1m4eZmRVk0H86NBUR+9L7HwJT6f3pwPeb6j2TyjqVP9OmvNM8jiJpHY0zF6amppidne0znPIdOnRoLNvdr1ziPPNjMAbf1aQsU3CsdTb0f6aLiJAUeTRm0HlExCZgE8D09HSM43+Nqut/u2qVS5wb18DlB3NpzyhNyjIFx1png97d9GzqKiL9fS6V7wWWNdVbmso6lS9tU95pHmZmVpBBk8R2YP4OpbXA7U3lV6a7nC4ADqYuox3ARZJOSResLwJ2pGEvSLog3dV0Zcu02s3DzMwK0rW7SdKXgRlgiaRnaNyldCNwq6SrgaeBd6fqdwJvBeaAnwLvBYiI/ZI+ATyQ6n08IuYvhn+Axh1UJwJ3pRcd5mFmZgXpmiQi4vKMQRe2qRvANRnT2QJsaVP+IHB2m/Ln283DzMyK419cm5lZJicJMzPL5CRhZmaZnCTMzCyTk4SZmWVykjAzs0xOEmZmlslJwszMMjlJmJlZJicJMzPL5CRh9bZxUdktMBtrThJ58c7IzGrIScJGbvn6O8pugpkNyEnCRu6pE95TdhPMbEBOEmZmlslJwszMMjlJWLX4BgCzSnGSMDOzTE4SZmaWyUmiBnyLqZmNipNEDfgWUzMbFScJsxL5LNCqzknCrEQ+C7Sqc5IwM7NMThJmZpbJScLMzDI5Sdh48i+zzQrhJGFmZpmcJCrKt0aaWRU4SVSUb400sypwkjAzs0xDJQlJ/03SI5L2SPqypBMknSHpfklzkr4i6bhU9/j0eS4NX940nY+k8sckXdxUvjqVzUlaP0xbzcysfwMnCUmnA38ATEfE2cCxwGXAp4BPR8TrgAPA1WmUq4EDqfzTqR6Szkrj/RqwGvi8pGMlHQt8DrgEOAu4PNU1M7OCDNvdtAA4UdIC4FeAfcCbgdvS8K3Apen9mvSZNPxCSUrl2yLiFxHxJDAHnJdecxHxRES8CGxLdc3MrCCKiMFHlq4DbgB+BvwtcB1wXzpbQNIy4K6IOFvSHmB1RDyThn0POB/YmMb5i1S+GbgrzWJ1RLwvlV8BnB8R17ZpxzpgHcDU1NTKbdu2DRzTwPbthtPOGXj0Q4cOsXDhwsGmN+S8R66pfUfF2aFuxzrQW8wlfjddY4XqL7se9RRrTdQ11lWrVu2KiOmjBkTEQC/gFOAe4NXAK4D/Dfw+jaP/+TrLgD3p/R5gadOw7wFLgM8Cv99Uvhl4Z3p9san8CuCz3dq1cuXKKMWGk4ca/d577x18ekPOe+Sa2ndUnB3qdqzTa8wlfjddY42o/rLrUU+x1kRdYwUejDb71GG6m94CPBkR/xQR/wL8DfAGYHHqfgJYCuxN7/empEEavgh4vrm8ZZyschvAOPzuYhzaaMXx+lANwySJfwQukPQr6drChcB3gHtpnAUArAVuT++3p8+k4fek7LUduCzd/XQGsAL4JvAAsCLdLXUcjYvb24dor5mNEf9WqBoWdK/SXkTcL+k24FvAS8DfA5uAO4Btkj6ZyjanUTYDfy5pDthPY6dPRDwi6VYaCeYl4JqI+CWApGuBHTTunNoSEY8M2t5J19jgDpbdjI68UzCrnoGTBEBEbAA2tBQ/QePOpNa6PwfelTGdG2hcAG8tvxO4c5g2mpnZ4PyLazMzy+QkYWZmmZwkrBgl/P8H3x1jNjwniTHnHWE2Xwg3G56TRJX1cPTtHaGZjZKThJmZZXKSMDOzTE4SZmaWyUnCbFyUcIeYmZOEmZllcpIwM7NMThJmZpbJScLMqsvXYUrnJGFmZpmcJCaJj8oq5fAjVbxcrMKcJMzMLJOTRA78kD0zqysnCTMzy+QkUTCfddg8P8HXxoGTRA763thHcaHSFz/NbAScJArmo0cbOz4AmWhOEmZjwN2UVhYnCaseH7kexWegw3OiHYyThBk4MU0AJ9rBOEmY1YmTXWf+fvrmJGFmZpmcJMwmlPvorRdOEmYTyn301gsniapzH6qZlchJwsaXE6jZyDlJ2FgZRT+6++bNsg2VJCQtlnSbpO9KelTSb0o6VdJOSY+nv6ekupJ0k6Q5SQ9JOrdpOmtT/cclrW0qXynp4TTOTZI0THvN2nHfvFm2Yc8kPgN8LSJ+FfgN4FFgPXB3RKwA7k6fAS4BVqTXOuBmAEmnAhuA84HzgA3ziSXVeX/TeKuHbK+Z1ZG7Hkdm4CQhaRHwJmAzQES8GBE/AtYAW1O1rcCl6f0a4JZouA9YLOk04GJgZ0Tsj4gDwE5gdRp2ckTcFxEB3NI0LTMzK4Aa+98BRpTOATYB36FxFrELuA7YGxGLUx0BByJisaSvAjdGxDfSsLuB64EZ4ISI+GQq/yjwM2A21X9LKn8jcH1EvK1NW9bRODthampq5bZt2waKaWD7djf+nnbOwHUPHTrEwoULj67Xbbr91uuljXlrmvdRcbarl+q28/Deg7z+mCe71mudb6/ty1PHWJvn28u8RxHLvt251esa6yDyjiWn5TySWCtg1apVuyJiurV8wRDTXACcC3wwIu6X9Ble7loCICJC0mBZqA8RsYlGwmJ6ejpmZmZGPcsjbVzT+Hv5wYHrzs7OckS75+t1m26/9XppY96a5n1UnO3qpbrtzLR2K+QR84i+m46xNs+3l3mPIpaNa3Kr1zXWQeQdS07LeSSxVtgw1ySeAZ6JiPvT59toJI1nU1cR6e9zafheYFnT+EtTWafypW3KzcysIAMniYj4IfB9SWemogtpdD1tB+bvUFoL3J7ebweuTHc5XQAcjIh9wA7gIkmnpAvWFwE70rAXJF2Quq2ubJqWmZkVYJjuJoAPAn8p6TjgCeC9NBLPrZKuBp4G3p3q3gm8FZgDfprqEhH7JX0CeCDV+3hE7E/vPwB8CTgRuCu9zKzmlq+/g6dOKLsVBkMmiYjYDRx1oYPGWUVr3QCuyZjOFmBLm/IHgbOHaaNZrWxcBBtLuK5kE8u/uM6T79U2K563u5FykjAzs0xOEmWZxKOfSYzZBuJHpVSHk4SZmWVykjAzs0xOEhXkR1ebWVU4SZiZWSYnCTMzy+QkYVYFBd/55S5N65WThNWfb709im8xtV45SVSQN2AzqwonCTMzy+QkYVYTvs5go+AkYWZmmZwkzMwsk5OEWR31ekeX7/yyLpwkzMwsk5PEGPOFSjMbNScJs5rw72tsFJwkzMwsk5OEmZllcpIws8niO7r64iSRtwJXQPdB20Tw7bylcpIokO9GKo6/a7N8OElYJXinblZNThIFcveQjRsnb3OSMDOzTE4Sk8YX96wPh89+vd5MLCcJq4S8u+LctWeWDycJG6mx6tP20bLZUZwkzMws09BJQtKxkv5e0lfT5zMk3S9pTtJXJB2Xyo9Pn+fS8OVN0/hIKn9M0sVN5atT2Zyk9cO21czM+pPHmcR1wKNNnz8FfDoiXgccAK5O5VcDB1L5p1M9JJ0FXAb8GrAa+HxKPMcCnwMuAc4CLk91bYz42oDZeBsqSUhaCvw28MX0WcCbgdtSla3Apen9mvSZNPzCVH8NsC0ifhERTwJzwHnpNRcRT0TEi8C2VNfMDOj9mtdYXRurGEXE4CNLtwF/CLwS+DBwFXBfOltA0jLgrog4W9IeYHVEPJOGfQ84H9iYxvmLVL4ZuCvNYnVEvC+VXwGcHxHXtmnHOmAdwNTU1Mpt27YNHNNA9u0+8vNp5/Rd79ChQyxcuDC36XWs26le3lraeOiVr3s5zg71ho65At/NEcu023y7zbuXNo4i5h7rdY21Xz3G8vDeg7z+mCdzq9eL3GOtiFWrVu2KiOnW8gWDTlDS24DnImKXpJkh2ja0iNgEbAKYnp6OmZmCm7Ox5QTn8oN915udneVwu3OYXse6nerlraWNszO303b55B1zBb6bI5Zpt/l2m3cvbRxFzD3W6xprv3qM5ar1d/DUCRtyq9eL3GOtuGG6m94AvF3SUzS6gt4MfAZYLGk++SwF9qb3e4FlAGn4IuD55vKWcbLKzV42AbetuqskW6/XvHxtbHADJ4mI+EhELI2I5TQuPN8TEb8H3Au8M1VbC9ye3m9Pn0nD74lGX9d24LJ099MZwArgm8ADwIp0t9RxaR7bB22vmZn1bxS/k7ge+JCkOeBVwOZUvhl4VSr/ELAeICIeAW4FvgN8DbgmIn4ZES8B1wI7aNw9dWuqa9a/CTjjMBuFga9JNIuIWWA2vX+Cxp1JrXV+DrwrY/wbgBvalN8J3JlHG23MbFwEGwu8dlJh7iqxMvkX12aTzmdZ1oGTxDjwRmxmJXGSKJN3/mZWcU4SZmaWyUmiTnxmYjZZCtjmnSTMzMZRQQeFThKj4CN6M6sJJwmzuvLBysCq/iiUItvnJGHFan26p1kFVf0HjEW2z0nCbNxM4hnCJMZcEU4SZmbtODEBThJmZgOr+rWLPDhJ2ESbhI3cRqfq1y7y4CRhZmaZnCTMzCyTk4RNtEnoLrD+9N0FWfML3E4SZmZNfOBwJCcJM6uPmh/Vl8FJwszMMjlJmE0Y3/aLzzj64CRhxfMG2p6/F6sgJwkzM8vkJGFmNoBJ6bZzkjCbML7F0/rhJGHjz335VoJJSbZOEmZmlslJwszMMjlJ1E2dul7qFIvZmHKSKFsBO8JJuQvDas4HDaVwkrDJ4Z1MX3xwUU1FLxcnCTOzYZV5ADLieQ+cJCQtk3SvpO9IekTSdan8VEk7JT2e/p6SyiXpJklzkh6SdG7TtNam+o9LWttUvlLSw2mcmyRpmGBHobCsPsSKMCm36plNgqK352HOJF4C/ntEnAVcAFwj6SxgPXB3RKwA7k6fAS4BVqTXOuBmaCQVYANwPnAesGE+saQ6728ab/UQ7S2WuzZslApYv47aGY1inlWfprfjwZNEROyLiG+l9z8GHgVOB9YAW1O1rcCl6f0a4JZouA9YLOk04GJgZ0Tsj4gDwE5gdRp2ckTcFxEB3NI0LTMzK4Aa+98hJyItB74OnA38Y0QsTuUCDkTEYklfBW6MiG+kYXcD1wMzwAkR8clU/lHgZ8Bsqv+WVP5G4PqIeFub+a+jcXbC1NTUym3btg0dU68e3nuQ1x/zZPuBp51z5Od9uzPrHTp0iIULF3at1+v0OtZrHT5KLfM+dPxrWPiLH+QXS971cvxujlim7bRrY7v5DxNzr9Mbph49xNriiO0m75h7rZe1rLusE4djLWu76ieWPqxatWpXREy3li8YdsKSFgJ/DfzXiHih+bJBRISk4bNQFxGxCdgEMD09HTMzM6Oe5WEznU5HLz945OeNazLrzc7OcrjdHer1Or2O9VqHj1LLvGfP/Bgzj23IL5a86z0GbMzn+zlimbbTro3t5j9MzL1Or9060Ws9eoi1xRHbTa/zbld3mHpZ20GX7eVwrGVtV/3EkoOh7m6S9AoaCeIvI+JvUvGzqauI9Pe5VL4XWNY0+tJU1ql8aZtyM8ub+94twzB3NwnYDDwaEX/aNGg7MH+H0lrg9qbyK9NdThcAByNiH7ADuEjSKemC9UXAjjTsBUkXpHld2TQtM6upcbhjcJIM0930BuAK4GFJu1PZ/wBuBG6VdDXwNPDuNOxO4K3AHPBT4L0AEbFf0ieAB1K9j0fE/vT+A8CXgBOBu9LLzMwKMnCSSBegs363cGGb+gFckzGtLcCWNuUP0rgYXg0bF+XWX21mNg6GvnA9CZavv8M/SDOzieTHcvTACcIGNmn93oPE2zLOxG9vFVtnnCTMRq1iG71ZP5wkzEbk4b2Tdf3KT42tJyeJSeQj20Jk/hK/pia+m6iDvhNohbZRJwlrr0IraZX4aLkHVVh3qtCGJk+d8J7KtalXThJmrTpszBN/tDymO7qhjMNTZUe4XJwkzIowiTvXDLU9G8tjGTdPoyLrjJNEF16hS5pelU1SrCPgs7HxWn+cJLoYyxV6zFbCsdTlO67twYVNHCeJcVLUzj+n+UzyjnIsDy4mUcu6PsnrbBYniYLUeuXzmYtZbTlJ2BGOSGZD7vx9NF0RvS5HJ/vsdXaCv0MniX4NuBKMyw5zXNo5sBpuxJOk1mfk7VRgfXWSsOFUYCUuTNVirVp7ClCVg5hJSlZOEtZZhx3RURtKnXZadYqlDP7+Ohuj78dJogLG9aikkKO6Km9MVW6bjVQ/6/64bt/znCQ6KGrhVuUUuki5f7d12mFXPZY27Zu0J972Y9y3byeJDoZduON+BDFKg3y3/j6ra9gn3nrZVpeTxAhV/gii6kesLSr/fdrAartsx2wba8dJwkpXu6PIHnYMtYvZClHGeuMkYVaC2h45dzHWybECZwVlrDdOEpavATaksd5hVmDHMU5yWdZV/85zbl/ZidVJYhBVX0mtksre2Avn7aQvWetH2QdRThLN6rRS1ymWmhjpxu7lfbSqfyct7Ss7GWRxkqiKfbsbf7ut2GWs+FXf2Oaldo7iiL20s4C8vvtxWYZ9yG2Z1PC7yZOTxKQbcAMpdKfZZxu7HpEN8ETPqh7lTfJ/GOx1mdSim6/E5eIk0Wp+YYzRxlKGibgA2au6xNGPXv8XcwW+m8om+DHhJNFOWSt2P/PtpW4FNlBrMS7LZBK7Ncuef0U5SQyqjv+EJK8jwlHEPC7fY9WvIYzL95iXqh/wjUHPxYKyG2ADKHuHXWV1TN5FKzNBVX35bVwEZ36s97p5z7sElT+TkLRa0mOS5iStL7s9R6j6Cl2munw3ZcaR97zL2mlVfRnDeLSxJJVOEpKOBT4HXAKcBVwu6axyW1VDGxf1d5Q3aRtU1XfWZZrEmPNW8e+m0kkCOA+Yi4gnIuJFYBuwpuQ2TbaKr9Bmli9FRNltyCTpncDqiHhf+nwFcH5EXNtSbx2wLn08E3is0IbmYwnwz2U3ogCTEic41rqqa6yvjYhXtxbW4sJ1RGwCNpXdjmFIejAipstux6hNSpzgWOtqkmKF6nc37QWWNX1emsrMzKwAVU8SDwArJJ0h6TjgMmB7yW0yM5sYle5uioiXJF0L7ACOBbZExCMlN2tUxrq7rA+TEic41rqapFirfeHazMzKVfXuJjMzK5GThJmZZXKSyIGkLZKek7QnY/gaSQ9J2i3pQUm/lcpfK+lbqfwRSf+laZyvSfp2Kv9C+vU5kjZK2pvG2S3prcVEebhducfaNO725ulKOlXSTkmPp7+njC6ytrEUGWvtlquk2fRInfmY/lUqP17SV9Kjdu6XtLyQICk8zqsk/VNT+fuKiTJnEeHXkC/gTcC5wJ6M4Qt5+frPrwPfTe+PA45vqvMU8Jr0+eT0V8BfA5elzxuBD9cp1lT2u8BfNU8X+CNgfXq/HvhUjWOt3XIFZoHpNtP6APCF9P4y4Cs1jfMq4LNlLdO8Xj6TyEFEfB3Y32H4oUhrDXASEKn8xYj4RSo/nqYzu4h4Ib1dQGMFrcQdBqOIVdJC4EPAJ1smtwbYmt5vBS4dtv39KDjWUo0i1g6al+ttwIWSNFDD+1RwnLUwMYGWTdLvSPoucAfwn5rKl0l6CPg+jSPlHzQN2wE8B/yYxsY079p0Sryl6C6YXgwQ6yeAPwF+2jKpqYjYl97/EJgabcv7l2OsUL/lCvBnqavlo02J4PRUl4h4CTgIvKqQIHqQY5wA70jL9DZJzT8MHh9ln8rU5QUsJ+MUtqXem4C/a1P+GuCbNHaMzeUn0Ohu+g/p8xSN34wcA9xA47cjYxsrcA6wvd10gR+1jHegxrHWarmmz6env68E/ha4Mn3eAyxtGu97wJIaxvkqXu6i+s/APUUv0zxePpMoWDROd/+1pCUt5T+gsfG8saX858DtpKffRsSzEfHLiPh/wP+k8aTcSuox1t8EpiU9BXwD+LeSZlPVZyWdBpD+PldQ0/s2bKw1XK5ExN7098c0rsHMx3T4cTuSFgCLgOcLaXwfho0zIp6Pl7uovgisLKjpuXKSKICk182fgko6l0af5vOSlko6MZWfAvwW8JikhU07xwXAbwPfTZ9Pa5r079BYWSuj31gj4uaIeE1ELE9l/xARM2ly24G16f1aGsmyMvKMtW7LVdKC+Z2rpFcAb+PlmJqX6ztpHGFX4ppbnnG2LNO3A48WF0l+Kv1YjnEh6cvADLBE0jPABuAVABHxBeAdwJWS/gX4GfAfIyIk/TvgTyQFjbuY/jgiHpY0BWyXNH+B7F7gC2l2fyTpHBoX1J6icRpbmLxj7TK7G4FbJV0NPA28exQxZSk41lotV0knATvSjvNY4O9onCEBbAb+XNIcjYvIl9U0zj+Q9HbgpRTnVQWFmSs/lsPMzDK5u8nMzDI5SZiZWSYnCTMzy+QkYWZmmZwkzMzGmLo8tLCl7msl3Z1+BT4raWm3cZwkzMzG25eA1T3W/WPgloj4deDjwB92G8FJwsxsjEWbhxZK+jdq/LuBXZL+j6RfTYPOAu5J7+8lPcmhEycJM7P62QR8MCJWAh8GPp/Kv03jUfXQ+GX/KyV1fLiif3FtZlYjajyO/t8D/6vpgbTHp78fBj4r6Srg6zSeo/XLTtNzkjAzq5djaDxB+ZzWAenhhL8Lh5PJOyLiR90mZmZmNRGNf1j2pKR3AajhN9L7JZLm9/sfAbZ0m56ThJnZGEsPLfy/wJmSnkkPxPw94GpJ3wYe4eUL1DM0nl77DzT+h8kNXafvB/yZmVkWn0mYmVkmJwkzM8vkJGFmZpmcJMzMLJOThJmZZXKSMDOzTE4SZmaW6f8DrYmZWOGcW+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(pd.to_datetime(data['startTime'].min(), unit='s'))\n",
    "print(pd.to_datetime(data['endTime'].max(), unit='s'))\n",
    "data['startTime'].hist(bins=1000)\n",
    "data['endTime'].hist(bins=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load in prefix and building information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "prefix_lat_lon = pd.read_csv('prefix_lat_lon_name_category.csv')\n",
    "legalPrefixes = prefix_lat_lon.prefix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    prefix        lat        lon                             name  category\n",
      "0      aaf  29.650275 -82.345720         Academic Advising Center     admin\n",
      "1      adv  29.648792 -82.359303  Martin H. Levin Advocacy Center  academic\n",
      "2      aer  29.643272 -82.348300       Mec. Aerospace Engineering  academic\n",
      "3     aerv  29.643272 -82.348300       Mec. Aerospace Engineering  academic\n",
      "4      alf  29.642948 -82.348636   Alfred A. Ring Tennis Pavilion    sports\n",
      "..     ...        ...        ...                              ...       ...\n",
      "189   weim  29.647761 -82.347342                      Wiemer Hall  academic\n",
      "190  weimG  29.647761 -82.347342                      Wiemer Hall  academic\n",
      "191  weimg  29.647761 -82.347342                      Wiemer Hall  academic\n",
      "192     wm  29.648052 -82.346709                  Williamson Hall  academic\n",
      "193    yon  29.649968 -82.347888                         Yon Hall  academic\n",
      "\n",
      "[194 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "print(prefix_lat_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     aaf\n",
       "1     adv\n",
       "2     aer\n",
       "3    aerv\n",
       "4     alf\n",
       "Name: prefix, dtype: object"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "legalPrefixes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clean the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prefixes found from our preprocessing function\n",
    "#prefixesFound = np.unique(data.prefix)\n",
    "prefixesFound = data['prefix'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prefixes in the base table that don't match in the prefix table \n",
    "#TODO: more preprocessing to keep some of these?\n",
    "unknownPrefixes = list(np.setdiff1d(prefixesFound, np.asarray(legalPrefixes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['AP',\n",
       " 'APc',\n",
       " 'APe',\n",
       " 'APf',\n",
       " 'apopka',\n",
       " 'apopka-win-lap',\n",
       " 'b',\n",
       " 'cba-ftlauderdale-win-lap',\n",
       " 'ctx',\n",
       " 'dairy-win-lap',\n",
       " 'ecob',\n",
       " 'escv',\n",
       " 'flrec',\n",
       " 'hotidf-win-lap',\n",
       " 'irrec',\n",
       " 'mowv',\n",
       " 'ohl',\n",
       " 'osbs-win-lap',\n",
       " 'osbs-xxx-win-lap',\n",
       " 'osbstr',\n",
       " 'pep',\n",
       " 'presAttic-wout-lap',\n",
       " 'presStairs-win-lap',\n",
       " 'rl',\n",
       " 'rubygillhs-win-lap',\n",
       " 'treeo-win-lap',\n",
       " 'uic-win-lap',\n",
       " 'unknown',\n",
       " 'wrt']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unknownPrefixes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(18867555, 8)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#keep only valid prefixes (based on the prefix_lat_lon_name_category.csv table)\n",
    "#removes about 400k rows?\n",
    "data = data[~data['prefix'].isin(unknownPrefixes)]\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9179058, 8)\n"
     ]
    },
    {
     "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>userIP</th>\n",
       "      <th>userMAC</th>\n",
       "      <th>APNAME</th>\n",
       "      <th>APMAC</th>\n",
       "      <th>startTime</th>\n",
       "      <th>endTime</th>\n",
       "      <th>deviceType</th>\n",
       "      <th>prefix</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>10.131.163.198</td>\n",
       "      <td>12544</td>\n",
       "      <td>gri101-wout-lap1231-1</td>\n",
       "      <td>00:1a:e2:10:ba:90</td>\n",
       "      <td>1333238683</td>\n",
       "      <td>1333238723</td>\n",
       "      <td>f</td>\n",
       "      <td>gri</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>10.132.190.177</td>\n",
       "      <td>33782</td>\n",
       "      <td>soc-temp-lap3502-5</td>\n",
       "      <td>88:f0:77:f2:67:a0</td>\n",
       "      <td>1333238693</td>\n",
       "      <td>1333238733</td>\n",
       "      <td>f</td>\n",
       "      <td>soc-temp-lap</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>10.132.251.191</td>\n",
       "      <td>85756</td>\n",
       "      <td>soc-temp-lap3502-2</td>\n",
       "      <td>88:f0:77:f2:5f:a0</td>\n",
       "      <td>1333238692</td>\n",
       "      <td>1333238736</td>\n",
       "      <td>f</td>\n",
       "      <td>soc-temp-lap</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>10.132.190.116</td>\n",
       "      <td>59170</td>\n",
       "      <td>soc-temp-lap3502-4</td>\n",
       "      <td>88:f0:77:f2:60:e0</td>\n",
       "      <td>1333238694</td>\n",
       "      <td>1333238738</td>\n",
       "      <td>f</td>\n",
       "      <td>soc-temp-lap</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>10.132.255.218</td>\n",
       "      <td>23629</td>\n",
       "      <td>soc-temp-lap3502-1</td>\n",
       "      <td>88:f0:77:f2:6b:80</td>\n",
       "      <td>1333238691</td>\n",
       "      <td>1333238741</td>\n",
       "      <td>unknown</td>\n",
       "      <td>soc-temp-lap</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            userIP  userMAC                 APNAME              APMAC  \\\n",
       "32  10.131.163.198    12544  gri101-wout-lap1231-1  00:1a:e2:10:ba:90   \n",
       "41  10.132.190.177    33782     soc-temp-lap3502-5  88:f0:77:f2:67:a0   \n",
       "42  10.132.251.191    85756     soc-temp-lap3502-2  88:f0:77:f2:5f:a0   \n",
       "45  10.132.190.116    59170     soc-temp-lap3502-4  88:f0:77:f2:60:e0   \n",
       "46  10.132.255.218    23629     soc-temp-lap3502-1  88:f0:77:f2:6b:80   \n",
       "\n",
       "     startTime     endTime deviceType        prefix  \n",
       "32  1333238683  1333238723          f           gri  \n",
       "41  1333238693  1333238733          f  soc-temp-lap  \n",
       "42  1333238692  1333238736          f  soc-temp-lap  \n",
       "45  1333238694  1333238738          f  soc-temp-lap  \n",
       "46  1333238691  1333238741    unknown  soc-temp-lap  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Keep only sessions longer than 10 seconds. Removes ~ 7 million receords. \n",
    "data = data[data['endTime']-data['startTime'] > 30]\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3722227, 8)\n"
     ]
    }
   ],
   "source": [
    "# Remove all records with users who have had fewer than 500 sessions \n",
    "data = data.groupby('userMAC').filter(lambda x : len(x)>500)\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get information about each user (useful for prediction later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "numUsers = data['userMAC'].unique().shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prefix</th>\n",
       "      <th>name</th>\n",
       "      <th>category</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>aaf</td>\n",
       "      <td>Academic Advising Center</td>\n",
       "      <td>admin</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>adv</td>\n",
       "      <td>Martin H. Levin Advocacy Center</td>\n",
       "      <td>academic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>aer</td>\n",
       "      <td>Mec. Aerospace Engineering</td>\n",
       "      <td>academic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>aerv</td>\n",
       "      <td>Mec. Aerospace Engineering</td>\n",
       "      <td>academic</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>alf</td>\n",
       "      <td>Alfred A. Ring Tennis Pavilion</td>\n",
       "      <td>sports</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  prefix                             name  category\n",
       "0    aaf         Academic Advising Center     admin\n",
       "1    adv  Martin H. Levin Advocacy Center  academic\n",
       "2    aer       Mec. Aerospace Engineering  academic\n",
       "3   aerv       Mec. Aerospace Engineering  academic\n",
       "4    alf   Alfred A. Ring Tennis Pavilion    sports"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "location_info = prefix_lat_lon[['prefix','name','category']]\n",
    "location_info.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Association Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(721,)\n"
     ]
    }
   ],
   "source": [
    "data_start = data['startTime'].min()\n",
    "data_end = data['endTime'].max()\n",
    "step = 3600\n",
    "time_slots = np.arange(data_start, data_end, step)\n",
    "#time_slots[-1] = data_end\n",
    "print(time_slots.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1333238682 1335830719\n",
      "[1333238682 1333242282 1333245882 1333249482 1333253082 1333256682\n",
      " 1333260282 1333263882 1333267482 1333271082 1333274682 1333278282\n",
      " 1333281882 1333285482 1333289082 1333292682 1333296282 1333299882\n",
      " 1333303482 1333307082 1333310682 1333314282 1333317882 1333321482\n",
      " 1333325082 1333328682 1333332282 1333335882 1333339482 1333343082\n",
      " 1333346682 1333350282 1333353882 1333357482 1333361082 1333364682\n",
      " 1333368282 1333371882 1333375482 1333379082 1333382682 1333386282\n",
      " 1333389882 1333393482 1333397082 1333400682 1333404282 1333407882\n",
      " 1333411482 1333415082 1333418682 1333422282 1333425882 1333429482\n",
      " 1333433082 1333436682 1333440282 1333443882 1333447482 1333451082\n",
      " 1333454682 1333458282 1333461882 1333465482 1333469082 1333472682\n",
      " 1333476282 1333479882 1333483482 1333487082 1333490682 1333494282\n",
      " 1333497882 1333501482 1333505082 1333508682 1333512282 1333515882\n",
      " 1333519482 1333523082 1333526682 1333530282 1333533882 1333537482\n",
      " 1333541082 1333544682 1333548282 1333551882 1333555482 1333559082\n",
      " 1333562682 1333566282 1333569882 1333573482 1333577082 1333580682\n",
      " 1333584282 1333587882 1333591482 1333595082 1333598682 1333602282\n",
      " 1333605882 1333609482 1333613082 1333616682 1333620282 1333623882\n",
      " 1333627482 1333631082 1333634682 1333638282 1333641882 1333645482\n",
      " 1333649082 1333652682 1333656282 1333659882 1333663482 1333667082\n",
      " 1333670682 1333674282 1333677882 1333681482 1333685082 1333688682\n",
      " 1333692282 1333695882 1333699482 1333703082 1333706682 1333710282\n",
      " 1333713882 1333717482 1333721082 1333724682 1333728282 1333731882\n",
      " 1333735482 1333739082 1333742682 1333746282 1333749882 1333753482\n",
      " 1333757082 1333760682 1333764282 1333767882 1333771482 1333775082\n",
      " 1333778682 1333782282 1333785882 1333789482 1333793082 1333796682\n",
      " 1333800282 1333803882 1333807482 1333811082 1333814682 1333818282\n",
      " 1333821882 1333825482 1333829082 1333832682 1333836282 1333839882\n",
      " 1333843482 1333847082 1333850682 1333854282 1333857882 1333861482\n",
      " 1333865082 1333868682 1333872282 1333875882 1333879482 1333883082\n",
      " 1333886682 1333890282 1333893882 1333897482 1333901082 1333904682\n",
      " 1333908282 1333911882 1333915482 1333919082 1333922682 1333926282\n",
      " 1333929882 1333933482 1333937082 1333940682 1333944282 1333947882\n",
      " 1333951482 1333955082 1333958682 1333962282 1333965882 1333969482\n",
      " 1333973082 1333976682 1333980282 1333983882 1333987482 1333991082\n",
      " 1333994682 1333998282 1334001882 1334005482 1334009082 1334012682\n",
      " 1334016282 1334019882 1334023482 1334027082 1334030682 1334034282\n",
      " 1334037882 1334041482 1334045082 1334048682 1334052282 1334055882\n",
      " 1334059482 1334063082 1334066682 1334070282 1334073882 1334077482\n",
      " 1334081082 1334084682 1334088282 1334091882 1334095482 1334099082\n",
      " 1334102682 1334106282 1334109882 1334113482 1334117082 1334120682\n",
      " 1334124282 1334127882 1334131482 1334135082 1334138682 1334142282\n",
      " 1334145882 1334149482 1334153082 1334156682 1334160282 1334163882\n",
      " 1334167482 1334171082 1334174682 1334178282 1334181882 1334185482\n",
      " 1334189082 1334192682 1334196282 1334199882 1334203482 1334207082\n",
      " 1334210682 1334214282 1334217882 1334221482 1334225082 1334228682\n",
      " 1334232282 1334235882 1334239482 1334243082 1334246682 1334250282\n",
      " 1334253882 1334257482 1334261082 1334264682 1334268282 1334271882\n",
      " 1334275482 1334279082 1334282682 1334286282 1334289882 1334293482\n",
      " 1334297082 1334300682 1334304282 1334307882 1334311482 1334315082\n",
      " 1334318682 1334322282 1334325882 1334329482 1334333082 1334336682\n",
      " 1334340282 1334343882 1334347482 1334351082 1334354682 1334358282\n",
      " 1334361882 1334365482 1334369082 1334372682 1334376282 1334379882\n",
      " 1334383482 1334387082 1334390682 1334394282 1334397882 1334401482\n",
      " 1334405082 1334408682 1334412282 1334415882 1334419482 1334423082\n",
      " 1334426682 1334430282 1334433882 1334437482 1334441082 1334444682\n",
      " 1334448282 1334451882 1334455482 1334459082 1334462682 1334466282\n",
      " 1334469882 1334473482 1334477082 1334480682 1334484282 1334487882\n",
      " 1334491482 1334495082 1334498682 1334502282 1334505882 1334509482\n",
      " 1334513082 1334516682 1334520282 1334523882 1334527482 1334531082\n",
      " 1334534682 1334538282 1334541882 1334545482 1334549082 1334552682\n",
      " 1334556282 1334559882 1334563482 1334567082 1334570682 1334574282\n",
      " 1334577882 1334581482 1334585082 1334588682 1334592282 1334595882\n",
      " 1334599482 1334603082 1334606682 1334610282 1334613882 1334617482\n",
      " 1334621082 1334624682 1334628282 1334631882 1334635482 1334639082\n",
      " 1334642682 1334646282 1334649882 1334653482 1334657082 1334660682\n",
      " 1334664282 1334667882 1334671482 1334675082 1334678682 1334682282\n",
      " 1334685882 1334689482 1334693082 1334696682 1334700282 1334703882\n",
      " 1334707482 1334711082 1334714682 1334718282 1334721882 1334725482\n",
      " 1334729082 1334732682 1334736282 1334739882 1334743482 1334747082\n",
      " 1334750682 1334754282 1334757882 1334761482 1334765082 1334768682\n",
      " 1334772282 1334775882 1334779482 1334783082 1334786682 1334790282\n",
      " 1334793882 1334797482 1334801082 1334804682 1334808282 1334811882\n",
      " 1334815482 1334819082 1334822682 1334826282 1334829882 1334833482\n",
      " 1334837082 1334840682 1334844282 1334847882 1334851482 1334855082\n",
      " 1334858682 1334862282 1334865882 1334869482 1334873082 1334876682\n",
      " 1334880282 1334883882 1334887482 1334891082 1334894682 1334898282\n",
      " 1334901882 1334905482 1334909082 1334912682 1334916282 1334919882\n",
      " 1334923482 1334927082 1334930682 1334934282 1334937882 1334941482\n",
      " 1334945082 1334948682 1334952282 1334955882 1334959482 1334963082\n",
      " 1334966682 1334970282 1334973882 1334977482 1334981082 1334984682\n",
      " 1334988282 1334991882 1334995482 1334999082 1335002682 1335006282\n",
      " 1335009882 1335013482 1335017082 1335020682 1335024282 1335027882\n",
      " 1335031482 1335035082 1335038682 1335042282 1335045882 1335049482\n",
      " 1335053082 1335056682 1335060282 1335063882 1335067482 1335071082\n",
      " 1335074682 1335078282 1335081882 1335085482 1335089082 1335092682\n",
      " 1335096282 1335099882 1335103482 1335107082 1335110682 1335114282\n",
      " 1335117882 1335121482 1335125082 1335128682 1335132282 1335135882\n",
      " 1335139482 1335143082 1335146682 1335150282 1335153882 1335157482\n",
      " 1335161082 1335164682 1335168282 1335171882 1335175482 1335179082\n",
      " 1335182682 1335186282 1335189882 1335193482 1335197082 1335200682\n",
      " 1335204282 1335207882 1335211482 1335215082 1335218682 1335222282\n",
      " 1335225882 1335229482 1335233082 1335236682 1335240282 1335243882\n",
      " 1335247482 1335251082 1335254682 1335258282 1335261882 1335265482\n",
      " 1335269082 1335272682 1335276282 1335279882 1335283482 1335287082\n",
      " 1335290682 1335294282 1335297882 1335301482 1335305082 1335308682\n",
      " 1335312282 1335315882 1335319482 1335323082 1335326682 1335330282\n",
      " 1335333882 1335337482 1335341082 1335344682 1335348282 1335351882\n",
      " 1335355482 1335359082 1335362682 1335366282 1335369882 1335373482\n",
      " 1335377082 1335380682 1335384282 1335387882 1335391482 1335395082\n",
      " 1335398682 1335402282 1335405882 1335409482 1335413082 1335416682\n",
      " 1335420282 1335423882 1335427482 1335431082 1335434682 1335438282\n",
      " 1335441882 1335445482 1335449082 1335452682 1335456282 1335459882\n",
      " 1335463482 1335467082 1335470682 1335474282 1335477882 1335481482\n",
      " 1335485082 1335488682 1335492282 1335495882 1335499482 1335503082\n",
      " 1335506682 1335510282 1335513882 1335517482 1335521082 1335524682\n",
      " 1335528282 1335531882 1335535482 1335539082 1335542682 1335546282\n",
      " 1335549882 1335553482 1335557082 1335560682 1335564282 1335567882\n",
      " 1335571482 1335575082 1335578682 1335582282 1335585882 1335589482\n",
      " 1335593082 1335596682 1335600282 1335603882 1335607482 1335611082\n",
      " 1335614682 1335618282 1335621882 1335625482 1335629082 1335632682\n",
      " 1335636282 1335639882 1335643482 1335647082 1335650682 1335654282\n",
      " 1335657882 1335661482 1335665082 1335668682 1335672282 1335675882\n",
      " 1335679482 1335683082 1335686682 1335690282 1335693882 1335697482\n",
      " 1335701082 1335704682 1335708282 1335711882 1335715482 1335719082\n",
      " 1335722682 1335726282 1335729882 1335733482 1335737082 1335740682\n",
      " 1335744282 1335747882 1335751482 1335755082 1335758682 1335762282\n",
      " 1335765882 1335769482 1335773082 1335776682 1335780282 1335783882\n",
      " 1335787482 1335791082 1335794682 1335798282 1335801882 1335805482\n",
      " 1335809082 1335812682 1335816282 1335819882 1335823482 1335827082\n",
      " 1335830682]\n"
     ]
    }
   ],
   "source": [
    "print(data_start, data_end)\n",
    "print(time_slots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "137\n"
     ]
    }
   ],
   "source": [
    "numLocations = location_info['name'].unique().shape[0]\n",
    "print(numLocations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def encounter_length(sa, ea, sb, eb):\n",
    "    smax = sa if sa > sb else sb\n",
    "    emin = ea if ea < eb else eb\n",
    "    return emin - smax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ncls import NCLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create dict mapping prefix to building name\n",
    "name_prefix_codes = dict(zip(location_info.prefix, location_info.name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateAssociationMatrices():\n",
    "    ass_mats = [0]*numUsers\n",
    "    usermacs = [0]*numUsers\n",
    "    inds = np.arange(0, len(time_slots))\n",
    "    #print(\"indx\", inds)\n",
    "    #print(time_slots)\n",
    "    int_tree = NCLS(time_slots.astype(np.int64), (time_slots + step).astype(np.int64), inds.astype(np.int64))\n",
    "\n",
    "    counter = 0\n",
    "    for user in data['userMAC'].unique():\n",
    "        if (counter % 5 == 0):\n",
    "            print(\"User \", counter, '/', numUsers)\n",
    "        #ass_mat = np.zeroes((time_slots.shape[0], numLocations))\n",
    "        ass_mat = pd.DataFrame(0, index=time_slots, columns=location_info['name'].unique())\n",
    "        userSessions = data[data['userMAC'] == user]\n",
    "        \n",
    "        for session in userSessions.iterrows():\n",
    "            start = session[1]['startTime']\n",
    "            end =  session[1]['endTime']\n",
    "            overlaps = int_tree.find_overlap(start, end)\n",
    "            for o in overlaps:\n",
    "                enc_len = encounter_length(start, end, o[0], o[1])\n",
    "                prefix = session[1]['prefix']\n",
    "                #sessionBuildingName = location_info[location_info['prefix'] == prefix] ['name']\n",
    "                sessionBuildingName = name_prefix_codes[prefix]\n",
    "                ass_mat.loc[o[0]][sessionBuildingName] += enc_len\n",
    "        \n",
    "        ass_mats[counter] = ass_mat\n",
    "        usermacs[counter] = user\n",
    "        counter+=1\n",
    "\n",
    "    return ass_mats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateBuildingVectors():\n",
    "    bldg_vecs = [0]*numUsers\n",
    "    usermacs  = [0]*numUsers\n",
    "    inds = np.arange(0, len(time_slots))\n",
    "    #print(\"indx\", inds)\n",
    "    #print(time_slots)\n",
    "    int_tree = NCLS(time_slots.astype(np.int64), (time_slots + step).astype(np.int64), inds.astype(np.int64))\n",
    "    np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)\n",
    "    counter = 0\n",
    "    for user in data['userMAC'].unique():\n",
    "        if (counter % 5 == 0):\n",
    "            print(\"User \", counter, '/', numUsers)\n",
    "        if counter > 2500:\n",
    "            break\n",
    "        #ass_mat = np.zeroes((time_slots.shape[0], numLocations))\n",
    "        ass_mat = pd.DataFrame(0, index=time_slots, columns=location_info['name'].unique())\n",
    "        userSessions = data[data['userMAC'] == user]\n",
    "        \n",
    "        for session in userSessions.iterrows():\n",
    "            start = session[1]['startTime']\n",
    "            end =  session[1]['endTime']\n",
    "            overlaps = int_tree.find_overlap(start, end)\n",
    "            for o in overlaps:\n",
    "                enc_len = encounter_length(start, end, o[0], o[1])\n",
    "                prefix = session[1]['prefix']\n",
    "                #sessionBuildingName = location_info[location_info['prefix'] == prefix] ['name']\n",
    "                sessionBuildingName = name_prefix_codes[prefix]\n",
    "                ass_mat.loc[o[0]][sessionBuildingName] += enc_len\n",
    "        ass_mat = np.asarray(ass_mat)\n",
    "        \n",
    "        tempBldgVec = np.zeros((len(time_slots), 1))   \n",
    "        for i in range(len(time_slots)): #loop all the time slots\n",
    "            if np.sum(ass_mat[i]) == 0: #havent been in any buildings for this time slot\n",
    "                if counter == 0: #first iter\n",
    "                    tempBldgVec[i] = np.random.randint(numLocations) #randomly pick a location \n",
    "                else: #still not in a building\n",
    "                    tempBldgVec[i] = bldg_vecs[counter-1][i] #take prev building\n",
    "            else:\n",
    "                tempBldgVec[i] = np.argmax(ass_mat[i]) #get building index\n",
    "        bldg_vecs[counter] = tempBldgVec\n",
    "        usermacs[counter]  = user\n",
    "        counter+=1\n",
    "\n",
    "    return bldg_vecs, usermacs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mats = generateAssociationMatrices()\n",
    "#Save association matrices (takes a long time to make them)\n",
    "#mats.to_pickle(\"./AssociationMatrices.pkl\")\n",
    "\n",
    "if loadPreprocessedData == False:\n",
    "    bldgHistoryPerPerson, usermacs = generateBuildingVectors()\n",
    "    with open('BldgHistoryPerPerson.npy', 'wb') as f:\n",
    "        np.save(f, bldgHistoryPerPerson)\n",
    "    with open('userMacInfo_BldgHistory.npy', 'wb') as f:\n",
    "        np.save(f, usermacs)\n",
    "else:\n",
    "    bldgHistoryPerPerson = np.load('BldgHistoryPerPerson.npy', allow_pickle=True)\n",
    "    usermacs = np.load('userMacInfo_BldgHistory.npy', allow_pickle=True)\n",
    "\n",
    "\n",
    "#DESCRIPTION: Each entry in the list is a user. Each user has an array over time of \n",
    "#             their most used building in that time slot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Normalize per row\n",
    "# def normalizeAssMat(mats):\n",
    "#     for i in range(0, len(mats)):\n",
    "#         mats[i] = mats[i].div(mats[i].sum(axis=1), axis=0)\n",
    "#         mats[i].fillna(value=0, inplace=True)\n",
    "#     return mats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalizedMats = normalizeAssMat(mats[:5])\n",
    "# print(normalizedMats[1].head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from sklearn.decomposition import PCA \n",
    "def ClusterWithAssMats():\n",
    "    #Clustering with ass_mats\n",
    "\n",
    "    numTimeSlots = len(time_slots)\n",
    "    numLocs = mats[0].shape[1]\n",
    "    X = [0]*100 #np.zeros((numUsers*numTimeSlots, numLocs))\n",
    "    for i in range(100):\n",
    "        X[i] = mats[i]\n",
    "    X = np.asarray(X)\n",
    "    print(X.shape)\n",
    "    X = X.reshape(100*numTimeSlots, numLocs)\n",
    "    print(X.shape)\n",
    "    \n",
    "    kmeans = KMeans(n_clusters = 5, random_state = 0).fit(X)\n",
    "    clusterCenters = kmeans.cluster_centers_\n",
    "    labels = kmeans.labels_ \n",
    "    print(\"labels\", labels)\n",
    "    \n",
    "    MC_Cluster = labels.reshape(100, 2881).T\n",
    "    unique, counts = np.unique(labels, return_counts=True)\n",
    "    print(unique, counts)\n",
    "    \n",
    "    #Try PCA on data and cluster it\n",
    "    pca = PCA(n_components = 3)\n",
    "    X_pca = pca.fit_transform(X)\n",
    "    print(pca.explained_variance_ratio_.cumsum())\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(projection='3d')\n",
    "    ax.scatter(X_pca[:,0], X_pca[:,1], X_pca[:,2])\n",
    "    plt.show()\n",
    "\n",
    "    kmeans_pca = KMeans(n_clusters = 5, random_state=5).fit(X_pca)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2500, 721)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Cluster the building time-series vector for each person\n",
    "#Temporarily use part of this array\n",
    "X = np.stack(bldgHistoryPerPerson[:2500],axis=1).squeeze().T\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Custom distance metric (based on building location)\n",
    "if loadPreprocessedData == False:\n",
    "    BldgsCoords = [0]*numLocations\n",
    "    unique_bldgs = prefix_lat_lon['name'].unique()\n",
    "\n",
    "    for i in range(numLocations):\n",
    "        BldgsCoords[i] = prefix_lat_lon[prefix_lat_lon.name == unique_bldgs[i]][['lat','lon']].iloc[0]\n",
    "\n",
    "    with open('BldgsCoords.npy', 'wb') as f:\n",
    "        np.save(f, BldgsCoords)\n",
    "    with open('BldgDict.npy', 'wb') as f:\n",
    "        np.save(f, bldg_dict)\n",
    "else:\n",
    "    BldgsCoords = np.load('BldgsCoords.npy', allow_pickle=True)\n",
    "    bldg_dict = np.load('BldgDict.npy', allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 29.650275, -82.34572 ],\n",
       "       [ 29.648792, -82.359303],\n",
       "       [ 29.643272, -82.3483  ],\n",
       "       [ 29.642948, -82.348636],\n",
       "       [ 29.651573, -82.341939],\n",
       "       [ 29.631292, -82.351768],\n",
       "       [ 29.631292, -82.351768],\n",
       "       [ 29.646905, -82.343942],\n",
       "       [ 29.647656, -82.341748],\n",
       "       [ 29.648887, -82.342909],\n",
       "       [ 29.643915, -82.344403],\n",
       "       [ 29.639509, -82.359952],\n",
       "       [ 29.639518, -82.359949],\n",
       "       [ 29.647961, -82.351393],\n",
       "       [ 29.649216, -82.359055],\n",
       "       [ 29.641892, -82.347876],\n",
       "       [ 29.646552, -82.342104],\n",
       "       [ 29.648906, -82.345784],\n",
       "       [ 29.648906, -82.345784],\n",
       "       [ 29.645354, -82.347903],\n",
       "       [ 29.645256, -82.34771 ],\n",
       "       [ 29.645442, -82.348053],\n",
       "       [ 29.64903 , -82.341675],\n",
       "       [ 29.64903 , -82.341675],\n",
       "       [ 29.652971, -82.342913],\n",
       "       [ 29.643928, -82.343066],\n",
       "       [ 29.64272 , -82.347873],\n",
       "       [ 29.650947, -82.344073],\n",
       "       [ 29.64635 , -82.346664],\n",
       "       [ 29.639071, -82.368803],\n",
       "       [ 29.65013 , -82.341602],\n",
       "       [ 29.643859, -82.344859],\n",
       "       [ 29.64821 , -82.344387],\n",
       "       [ 29.633672, -82.361716],\n",
       "       [ 29.648189, -82.348498],\n",
       "       [ 29.649895, -82.344979],\n",
       "       [ 29.647025, -82.341302],\n",
       "       [ 29.641637, -82.348364],\n",
       "       [ 29.64353 , -82.343826],\n",
       "       [ 29.643476, -82.366184],\n",
       "       [ 29.634543, -82.366474],\n",
       "       [ 29.648127, -82.340636],\n",
       "       [ 29.648118, -82.339764],\n",
       "       [ 29.647768, -82.340343],\n",
       "       [ 29.648372, -82.340171],\n",
       "       [ 29.638503, -82.357321],\n",
       "       [ 29.638525, -82.35733 ],\n",
       "       [ 29.63847 , -82.361135],\n",
       "       [ 29.649463, -82.347294],\n",
       "       [ 29.651575, -82.34385 ],\n",
       "       [ 29.649981, -82.332703],\n",
       "       [ 29.650131, -82.343723],\n",
       "       [ 29.65047 , -82.347275],\n",
       "       [ 29.648006, -82.350075],\n",
       "       [ 29.650745, -82.340119],\n",
       "       [ 29.649287, -82.342056],\n",
       "       [ 29.637031, -82.370247],\n",
       "       [ 29.651103, -82.3486  ],\n",
       "       [ 29.650423, -82.340748],\n",
       "       [ 29.649817, -82.359215],\n",
       "       [ 29.651757, -82.334063],\n",
       "       [ 29.648251, -82.345586],\n",
       "       [ 29.644505, -82.352442],\n",
       "       [ 29.652301, -82.342142],\n",
       "       [ 29.652293, -82.341916],\n",
       "       [ 29.63781 , -82.362908],\n",
       "       [ 29.645697, -82.343192],\n",
       "       [ 29.649513, -82.346505],\n",
       "       [ 29.644393, -82.342214],\n",
       "       [ 29.643124, -82.347321],\n",
       "       [ 29.65112 , -82.34288 ],\n",
       "       [ 29.648846, -82.340697],\n",
       "       [ 29.64668 , -82.349201],\n",
       "       [ 29.647327, -82.349466],\n",
       "       [ 29.651649, -82.341016],\n",
       "       [ 29.646982, -82.345531],\n",
       "       [ 29.646754, -82.34492 ],\n",
       "       [ 29.646682, -82.345338],\n",
       "       [ 29.64646 , -82.344389],\n",
       "       [ 29.636058, -82.370845],\n",
       "       [ 29.648039, -82.346425],\n",
       "       [ 29.639963, -82.362885],\n",
       "       [ 29.641979, -82.344079],\n",
       "       [ 29.632906, -82.361672],\n",
       "       [ 29.647957, -82.343904],\n",
       "       [ 29.648244, -82.342908],\n",
       "       [ 29.644484, -82.3441  ],\n",
       "       [ 29.642764, -82.346123],\n",
       "       [ 29.642287, -82.347076],\n",
       "       [ 29.643913, -82.349779],\n",
       "       [ 29.646251, -82.337966],\n",
       "       [ 29.646947, -82.338164],\n",
       "       [ 29.646645, -82.337581],\n",
       "       [ 29.647381, -82.347996],\n",
       "       [ 29.645392, -82.345778],\n",
       "       [ 29.639754, -82.358956],\n",
       "       [ 29.639824, -82.358554],\n",
       "       [ 29.635312, -82.369309],\n",
       "       [ 29.642948, -82.369274],\n",
       "       [ 29.644585, -82.34878 ],\n",
       "       [ 29.639079, -82.358344],\n",
       "       [ 29.635726, -82.370236],\n",
       "       [ 29.648184, -82.36085 ],\n",
       "       [ 29.639481, -82.359603],\n",
       "       [ 29.641504, -82.347885],\n",
       "       [ 29.643418, -82.346266],\n",
       "       [ 29.649468, -82.345665],\n",
       "       [ 29.650076, -82.34634 ],\n",
       "       [ 29.646386, -82.34779 ],\n",
       "       [ 29.650148, -82.346799],\n",
       "       [ 29.647253, -82.348948],\n",
       "       [ 29.633165, -82.361168],\n",
       "       [ 29.633176, -82.360914],\n",
       "       [ 29.647983, -82.34811 ],\n",
       "       [ 29.646965, -82.343196],\n",
       "       [ 29.644495, -82.345433],\n",
       "       [ 29.649259, -82.344494],\n",
       "       [ 29.647647, -82.342621],\n",
       "       [ 29.650019, -82.344212],\n",
       "       [ 29.650878, -82.341791],\n",
       "       [ 29.649461, -82.351125],\n",
       "       [ 29.649961, -82.349511],\n",
       "       [ 29.650832, -82.341177],\n",
       "       [ 29.650976, -82.340971],\n",
       "       [ 29.638378, -82.368414],\n",
       "       [ 29.639601, -82.360412],\n",
       "       [ 29.649397, -82.339815],\n",
       "       [ 29.649333, -82.343903],\n",
       "       [ 29.633778, -82.366198],\n",
       "       [ 29.650447, -82.346867],\n",
       "       [ 29.648549, -82.351042],\n",
       "       [ 29.649455, -82.341582],\n",
       "       [ 29.648418, -82.348606],\n",
       "       [ 29.648399, -82.349265],\n",
       "       [ 29.647761, -82.347342],\n",
       "       [ 29.648052, -82.346709],\n",
       "       [ 29.649968, -82.347888]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BldgsCoords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyclustering.cluster.kmeans import kmeans\n",
    "from pyclustering.cluster.center_initializer import kmeans_plusplus_initializer\n",
    "from pyclustering.cluster.elbow import elbow\n",
    "from pyclustering.utils import read_sample\n",
    "from pyclustering.samples.definitions import SIMPLE_SAMPLES\n",
    "\n",
    "from pyclustering.cluster.kmedians import kmedians\n",
    "\n",
    "from haversine import haversine, Unit\n",
    "import math\n",
    "import pyclustering.utils.metric\n",
    "from pyclustering.utils.metric import distance_metric, type_metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def building_dist_metric(a, b):\n",
    "    coords_a = [BldgsCoords[int(x)] for x in a]\n",
    "    coords_b = [BldgsCoords[int(x)] for x in a]\n",
    "    \n",
    "    totalDist = 0\n",
    "    \n",
    "    for i in range(0,len(coords_a)):\n",
    "        totalDist += haversine(coords_a[i], coords_b[i])\n",
    "        \n",
    "    return math.sqrt(totalDist)\n",
    "bldg_dist_func = lambda point1, point2: building_dist_metric(point1, point2)\n",
    "bldg_dist_metric = distance_metric(type_metric.USER_DEFINED, func=bldg_dist_func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "object of type 'numpy.int64' has no len()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m--------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-43-bcc3a4fab639>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m# Run cluster analysis and obtain results.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mkmedians_instance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprocess\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      9\u001b[0m \u001b[0mclusters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkmedians_instance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_clusters\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     10\u001b[0m \u001b[0mmedoids\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkmedians_instance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_medoids\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~/opt/anaconda3/lib/python3.7/site-packages/pyclustering/cluster/kmedians.py\u001b[0m in \u001b[0;36mprocess\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    112\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    113\u001b[0m             \u001b[0;31m# Check for dimension\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 114\u001b[0;31m             \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__pointer_data\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__medians\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\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    115\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0mNameError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Dimension of the input data and dimension of the initial medians must be equal.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    116\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: object of type 'numpy.int64' has no len()"
     ]
    }
   ],
   "source": [
    "# Initialize initial medoids using K-Means++ algorithm\n",
    "initial_medoids = kmeans_plusplus_initializer(X, 4).initialize(return_index=True)\n",
    " \n",
    "# Create instance of K-Medoids (PAM) algorithm.\n",
    "kmedians_instance = kmedians(X, initial_medoids, metric = bldg_dist_metric)\n",
    "\n",
    "# Run cluster analysis and obtain results.\n",
    "kmedians_instance.process()\n",
    "clusters = kmedians_instance.get_clusters()\n",
    "medoids = kmedians_instance.get_medoids()\n",
    " \n",
    "# Print allocated clusters.\n",
    "print(\"medoids:\", medoids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished elbow\n"
     ]
    }
   ],
   "source": [
    "# create instance of Elbow method using K value from 1 to 10.\n",
    "kmin, kmax = 1, 75\n",
    "elbow_instance = elbow(X, kmin, kmax)\n",
    " \n",
    "# process input data and obtain results of analysis\n",
    "elbow_instance.process()\n",
    "amount_clusters = elbow_instance.get_amount()  # most probable amount of clusters\n",
    "wce = elbow_instance.get_wce()  # total within-cluster errors for each K\n",
    "print(\"finished elbow\")\n",
    "manhattan_metric = distance_metric(type_metric.MANHATTAN)\n",
    "\n",
    "# perform cluster analysis using K-Means algorithm\n",
    "\n",
    "centers = kmeans_plusplus_initializer(X, amount_clusters,\n",
    "                      amount_candidates=kmeans_plusplus_initializer.FARTHEST_CENTER_CANDIDATE,\n",
    "                      metric = bldg_dist_metric).initialize()\n",
    "kmeans_instance = kmeans(X, centers)\n",
    "kmeans_instance.process()\n",
    "\n",
    "# obtain clustering results and visualize them\n",
    "clusters = kmeans_instance.get_clusters()\n",
    "means    = kmeans_instance.get_centers()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/nkroeger/opt/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray\n",
      "  return array(a, dtype, copy=False, order=order, subok=True)\n"
     ]
    }
   ],
   "source": [
    "with open('kmeans_instance.npy', 'wb') as f:\n",
    "    np.save(f, kmeans_instance)\n",
    "with open('clusters.npy', 'wb') as f:\n",
    "    np.save(f, clusters)\n",
    "with open('means.npy', 'wb') as f:\n",
    "    np.save(f, means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "roundedMeans = [np.rint(x) for x in means]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEWCAYAAACEz/viAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4lElEQVR4nO3de7zUc/7A8de7u66oWEJh12WXRCdFqFPZUiTKbZXLhhBySSs6+7MsbWflkojcryHRDSV1KiV02dzKJRQWK9dSKafevz/e3zHfM2fOOTM1c2ZO834+HvOY+V7mO5858X3P5/b+iKrinHMut1XLdAGcc85lngcD55xzHgycc855MHDOOYcHA+ecc3gwcM45hwcDVwEROUdE5oW2VUR+n8kypUoqv4uIrBSRLqm4VqaJyJki8nKarj1bRM4r49j1IvJ4Oj7XVcyDgYvcyDaIyM+hx+hMlwt+C0YqIrfF7D8x2P9wgtcp8yaUbiLysIhsivn7npaia9cWkeEi8lnwb/iRiFwtIpLg+1sEf8cakX2q+oSq/jkV5XNVR42KT3E54gRVfSXThSjDx8CpInK1qhYH+84GPsxgmZJVqKrDtvbNIlIj9N3DxgO/A7oD7wN5wGPAnsBlW/t5Lvd4zcBtje4i8omIfCsi/xaRagAiUk1EhonIKhH5RkQeFZFGwbFHROSq4HWz4NfowGB7XxH5PnKdOL4G3gG6BufvDBwJTA6fJCLtROQ1EflRRN4SkY7B/puAo4HRcWo9XYJf0z+KyF2RX9TlfZfgeL/g2Hcict3W/iFF5HwRWRF8/8kisnvomIrIQBH5CPgozns7A38Geqvqu6parKqvA32BgZEmsKBWNFxE3hSRNSIyKfgbAswNnn8M/jZHlNE0eHHwd1orIjcG/2avBdd7RkRqBefuJCJTRWS1iPwQvN5jK/4uNUVknIhMiFzbpZcHA7c1TsJ+gR4GnAj8Ndh/TvDIB/YB6gORG+8coGPwugPwCXBMaPtVVd1Szmc+CpwVvD4dmARsjBwUkWbAC8A/gZ2BwcAEEWmqqtcBrwKXqGp9Vb0kdN3jgTZAS+BUgoBT3ncRkT8CY4B+wO5AY2BrbnidgOHB5+4GrAKeijmtF9AW+GOcSxwLvKGqn4d3quobwBdA59Dus7B/p92AYmBUsD/yb7Bj8LdZUEZxuwKtgXbAEGAsFnT2BA4CzgjOqwY8BDQH9gI2EP1vICEisgMwEfv3PVVVNyXzfrd1qmwwEJEHg19s7yZwbnMRmSkibwe/kpL+HzcHTAx+HUce55dz7ghV/V5VPwNuJ3ojOBO4VVU/UdWfgaHA6UF79BzgqODX/zFAIdA+eF+H4Hh5ngc6Br/Oz8KCQ1hf4EVVfVFVt6jqDGAR1nxSnn+p6o/BdykCWiXwXfoAU1V1rqpuBAqA8gIZwODQ3/bb0Gc8qKpLgusMBY4QkRah9w0P/tYb4lyzCfBVGZ/3VXA84rGg9rAuKO+pIlK9gjKHFarqGlV9D3gXeDn42/wEvAQcCqCq36nqBFVdr6prgZuwf99ENQSmYU2D56rq5iTe67ZBlQ0GwMNAtwTPvQV4VFVbAjdgv8ZcSb1UdcfQ475yzg3/El2F/TomeF4Vc6wGsKuqfgysw262RwNTgS9FZH8SCAbBzfAFYBjQWFXnx5zSHDglHNCAo7BfwuX5OvR6PVYDKPe7BMd++xsEN9jvKvicW0J/28hNusRnBEHnO6BZ6H0lfvXH+Jayv99uwfF411kF1KRksKjI/0KvN8TZrg8gInVF5N6gCW0N1gy1YxKBpx1WS/uXehbNSlVlg4GqzgW+D+8L2jGnichiEXlVRA4IDv0RmBW8LsKaNtzW2zP0ei/gy+D1l9hNOXysmOiNYw72q7qWqv432D4b2AlYmsDnPgpcBcQbfvg59us3HNDqqeq/guPJ3ljK+y5fEfobiEhdrKkoWSU+Q0TqBdf5b+ic8sr9CtBWRML/HohI26B8s0K7Y//NfsWCRapvuFcB+wNtVbUh0WaohEY3AS9jP9ZmisiuKS6bK0eVDQZlGAtcqqqtsTbju4P9bwEnB69PAhqIyNb8z+vM1UFH4Z7AIODpYP844AoR2VtE6gM3A0+HRsHMAS4h2mk5O9iel2BzwBysnfzOOMceB04Qka4iUl1E6ohIx1CT4P+wtv9ElfddngWOF5Gjgs7NG9i6/5fGAeeKSCsRqR18xhuqujKRNwejv2ZifSN/Cr53O+xvMUZVw53OfUXkj0HgugF4Nvibr8aauJL525SnAVZT+DHopP6/ZC+gqoXAk1hASKb24rbBdhMMgv9hjwTGi8hS4F6iVejBQAcR+Q/WJPFfwNsiS5oiJcfBP1/OuZOAxdiv+ReAB4L9D2LDGucCnwK/AJeG3jcHu1lEgsE8oG5ou1xqZqrq93GOfY7V+K7FbnCfA1cT/W/8DqBPMMJlVOz74yjzuwTt5gOxG9ZXwA9Yh21Sgpt5ATAhuM6+WOd4Mnpjtd1pwM9YIHiAkn93gu/yMNYsVodg2Kmqrsfa9ecHzWvtkv0eMW4HdsBqHa8H5Uqaqt6IdSK/Ehr55NJIqnKzXNDRNlVVDxKRhsAHqlpuG3EQNN5XVe9EdjlBRGYDj6vq/Zkui8te203NQFXXAJ+KyCkAYg4JXjeR6Bj2odivPuecc4EqGwxEZBywANhfRL4Qkf7YUL3+IvIW8B7RjuKOwAci8iE2GuSmDBTZOeeyVpVuJnLOOZcaVbZm4JxzLnWqZKK6Jk2aaIsWLTJdDOecq1IWL178rao2jXesSgaDFi1asGjRokwXwznnqhQRWVXWMW8mcs4558HAOeecBwPnnHN4MHDOOYcHA+ecc+RIMCgshKKikvuKimy/c865HAkGbdrAqadGA0JRkW23aZPZcjnnXLaokvMMkpWfD888AyefDAceCB99ZNv5+ZkumXPOZYecqBmA3fjbt4cFC+D00z0QOOdcWM4Eg6IimB+smvvoo6X7EJxzLpflRDCI9BE88YRtn3xyyT4E55zLdWkNBiKyp4gUicgyEXlPRAaVc24bESkWkT6pLsfChdZH0L077LMPrFtn2wsXpvqTnHOuakp3B3IxcJWqLhGRBsBiEZmhqsvCJ4lIdWAE8HI6CjFkSPR1y5bw9tvWZ+D9Bs45Z9JaM1DVr1R1SfB6LbAcaBbn1EuxRcG/SWd5AA4+2EYTbdiQ7k9yzrmqo9L6DILF6w8F3ojZ3ww4CRhTwfsvEJFFIrJo9erVW12Oli1hyxZYtqzic51zLldUSjAQkfrYL//Lg4Xrw24H/qaqW8q7hqqOVdU8Vc1r2jTu2gwJOfhge37nna2+hHPObXfSPulMRGpigeAJVX0uzil5wFMiAtAE6C4ixao6MR3l+f3voU4d6zdwzjln0hoMxO7wDwDLVfXWeOeo6t6h8x8GpqYrEABUrw5/+pPXDJxzLizdNYP2QD/gHRFZGuy7FtgLQFXvSfPnx9WyJbz4YiY+2TnnslNag4GqzgMkifPPSV9pog4+GB56CL75BnbZpTI+0TnnsltOzECO1bKlPXtTkXPOmZwMBpERRd6J7JxzJieDwS67wK67es3AOecicjIYgNUOvGbgnHMmZ4NBy5bw3nuweXOmS+Kcc5mXs8Hg4IPhl19gxYpMl8Q55zIvZ4OBjyhyzrmonA0GBx4I1ap5v4FzzkGOBoPCQnj9ddhvv2jNoKjI9jvnXC7KyWDQpo0te7nrrlYziCyL2aZNpkvmnHOZkfaspdkoP9+WvTz+eFi/Hk45BcaP95XPnHO5KydrBmA3/p497XXXrh4InHO5LWeDQVERvPIK1K4NEybYtnPO5aqcDAaRPoJnnoEzz7RRRaee6gHBOZe7cjIYLFxogSA/34LBhg1w0UW23znnclFOdiAPGRJ93aED7L47vPUWTJqUuTI551wm5WTNIKx6dTjjDHjpJfjuu0yXxjnnMiPngwFYU9Gvv9rwUuecy0UeDIBWrSw9xRNPZLokzjmXGR4MABGrHcybB6tWZbo0zjlX+TwYBNautecnn4zu83xFzrlc4cEg0LUr1KgBY8eCqucrcs7lloSDgYi0F5F6weu+InKriDRPX9EqV34+DBwIK1fa6KLIpDRPU+GcywXJ1AzGAOtF5BDgKuBj4NG0lCpD/vEPaNoUnn4aGjaELVtKHvdmI+fc9iqZYFCsqgqcCIxW1buABukpVmYsWWIB4KCD4JNPoFs3mDzZjnmzkXNue5ZMMFgrIkOBfsALIlINqJmeYlW+yM1+/Hhb42DQICguhl694NJLvdnIObd9SyYYnAZsBP6qql8DewD/TkupMiCcr0gEbr8d7rrLXo8ebbmLPBA457ZXCQeDIABMAGoHu74Fnk9HoTJhyJDSN/sDD4RatSxlxd13e1ZT59z2K5nRROcDzwL3BruaARPTUKasEGk2uvNO2LzZXnuaa+fc9iqZZqKBQHtgDYCqfgTsko5CZYNIs9F550H79jBjBjz1lKe5ds5tn5JJYb1RVTeJCAAiUgPQtJQqC4TTXF98saWr2LKl5H7nnNteJFMzmCMi1wI7iMixwHhgSnqKlV1697b5B3ffnemSOOdceiQTDK4BVgPvAAOAF4Fh6ShUtqldG/r3tzkHn3+e6dI451zqJRMMdgAeVNVTVLUP8GCwLyeoWjPRffdF9/mMZOfc9iKZYDCTkjf/HYBXUluc7NW1qw0zHT0aNm2CAQNsQlp4RrIHB+dcVZVMMKijqj9HNoLXdVNfpOyUnw/XXw8//GDrJo8bZxPSIjxdhXOuKksmGKwTkcMiGyLSGtiQ+iJlryFDYN994fXXbf2DXXaB44+HwYM9XYVzrmpLJhhcDowXkVdFZB7wNHBJWkqVpebOhZ9+gksugXr1bGby+vUwcqSnq3DOVW0JzzNQ1YUicgCwf7DrA1X9NT3Fyj6RZqDIr/+TT7Y+gzp14JdfbKZyfr4HBOdc1ZTsSmdtgJbAYcAZInJWeSeLyJ4iUiQiy0TkPREZFOecM0XkbRF5R0ReC9ZLyDrhRHYRInDaaVZLOPxwT1fhnKu6Eq4ZiMhjwL7AUmBzsFspf4GbYuAqVV0iIg2AxSIyQ1WXhc75FOigqj+IyHHAWKBtEt+hUsTOPF64EJ5/3oJDw4Zwzz3wxBO232sHzrmqJpl0FHnAH4MFbhKiql8BXwWv14rIcizB3bLQOa+F3vI6lho764WDw6BBNuR06VK46SbbV1hoI4vCgaGoyIKFp7RwzmWbZJqJ3gV+t7UfJCItgEOBN8o5rT/wUhnvv0BEFonIotWrV29tMdJi332t/+Cee2DdOtvXpk3JZiMfeuqcy2bJBIMmwDIRmS4ikyOPRN4oIvWxtRAuV9U1ZZyTjwWDv8U7rqpjVTVPVfOaNm2aRLErx5VXwvffw6NBo1l+vmU57dkTevTwoafOueyWTDPR9VvzASJSEwsET6jqc2Wc0xK4HzhOVb/bms/JtPnzYf/94bbbbHbyu+/acpk//wwvvmhNSR4InHPZKpmhpXOSvbhYvusHgOWqemsZ5+wFPAf0U9UPk/2MbHH44dZfsHatZTmdPNlyGUWGnt53H5x4ogcE51x2Smals3YislBEfhaRTSKyWUTiNvmEtAf6AZ1EZGnw6C4iF4rIhcE5fwcaA3cHxxdt3VfJrPx8eO45qFYNJk60CWkNG1qtoFs3qFvXh54657JXMs1Eo4HTsXUM8oCzgP3Ke4OqzgOkgnPOA85LohxZq0sXuOAC60g+5hi47joLEmvWWAfzjTf60FPnXHZKatKZqq4AqqvqZlV9COiWnmJVTUVF8OyzUFAAb70V3d+jBzRrBvPm+bBS51x2SiYYrBeRWsBSESkUkSuSfP92LZyu4oYb7DnSLFSjBpx/PkyfDp98kumSOudcacnczPsF518CrAP2BE5OR6Gqoth0Ffn5tr1woW2fd571I4wdm7kyOudcWSTRCcUiMkhV76hoX2XIy8vTRYuqVj9zYSFMmQIffGBLZ9au7TOSnXOVS0QWq2pevGPJ1AzOjrPvnK0qUQ5q0wbeeQdWr7acRrEzkgsLS4808pXTnHOVpcJgICJniMgUYO/wzGMRmQ18n/YSbify82HCBBt6evXVpWcke/oK51wmJTK09DUs2VwTYGRo/1rg7XQUanvVubOtpfzSS9C6NXTsGD0W6WM45RRbVnPuXE9f4ZyrPBXWDFR1larOBroArwYzkb/CsouWO4fAlRTpIzjiCFi8GI46CsJdNu+/b4nunnsOataEadO86cg5VzmS6TOYC9QRkWbAy9jooofTUajtUXjo6fz5Vit47TULCJs3Wz6jiy+GjRtt8trXX8OoUTZZzZuOnHPplkwwEFVdjw0nvVtVTwH+lJ5ibX/CQ09FYNasaEBo0sSGnNaoYaksZsyAxx+H4mKbvdyjBwwb5plPnXPpk1QwEJEjgDOBF4J91VNfpO3TkCGll8ycNQuOPBJ+/NEynk6fbimvAf7yF5g61Y5v2GBJ8C66yAOBcy49kgkGlwNDgedV9T0R2QfwtGvbYPZs+PBDy2H03XcWIMK6drV8RtWqWbAYM8YT3Tnn0iPhYKCqc1S1p6qOCLY/UdXL0le07Vu4D+Gf/yyZviJ8zmmnwQknwGef2cI5nvnUOZcOFQ4tFZHbVfXyYK5BqenKqtozLSXbzpWXviKyL3JOcTFMmgS//lr6HOecS4UK01GISGtVXSwiHeId35pFb7ZVVUxHsS02bYKmTaFPH3jggUyXxjlXVZWXjqLCmoGqLg6eK/2m70ytWjaiaMoUG4Za3bvtnXMplkgz0TvEaR6KUNWWKS2Ri+vEE2HcOFiwwOYmOOdcKiXSgXw8cAIwLXicGTxeAl5MX9Fc2HHH2azkSZOi+zy5nXMuVRJNR7EKOFZVh6jqO8Hjb8Cf019EB7aecqdONikt0s3jye2cc6mS7KSz9qGNI5N8v9tGvXrBihWwfLlt5+fbTOUePWz46fHHw9ChJUcaeU3BOZeIZG7m/YG7RWSliKwE7gb+mpZSubgis5MnTrTn77+Hm2+2GcpTp1o6i3/8w2Y2g9cUnHOJS2bS2WJVPQQ4BDhEVVup6pLIcRGJt/iNS6HHH7eZyJMmwcqV0L49zJtnq6adcoplPF2zxmYun3uu5zJyziUu6WYeVf1JVX+Kc2hQCsrjytGmjS2Z+eab0LatBQRVqx0884ylvK5f3zqaH37YhqT+/HPJa3izkXMunlS2+fvaBmmWnw+jR9vrX36xXEa33AJXXmn7unSxZiIRCxZffmlDUv/xDzvuzUbOubKkMhiUP5XZpcS551pG0zVrYPDgaCAAu9kPH279B6+/bontROD66yEvz5uNnHNl85pBFVNUBC+/DAUFpbOYxuY7uvBCazraYw9bWa1TJw8Ezrn4UhkM5qfwWi6OcKbTG24onek0ds0EsBFGv/wCu+4K48fD5MmVX27nXPZLJB3FleUdV9Vbg+dLUlUoF18imU7DwsFjxx2jTUUvvZS+GkJhofVJxM51WLjQgpVzLjslUjNoEDzygIuAZsHjQuCw9BXNxYr3yz8/v+ybbDh4HHqoLbO5cSM89lj0nFSPLvJZ0c5VTRWmsP7tRJG5QA9VXRtsNwBeUNVj0li+uHIthXWqTJ9us5V32slmMi9Zkp5O5VmzoFs3+POf4Y03vNPauWxRXgrrZPoMdgU2hbY3BftcFdG1K4waBd9+C+3aWWBIR/qKvfayhXheeMHXbXauqkgmGDwKvCki14vI9cAbwCNpKZVLm4svhg4d4P33LY3FtdfCnXfa5LVUNencc48916jh6zY7V1Ukk47iJuBc4Ifgca6q3pyugrn0KCqC996DgQNhhx1sHsJll0GLFqlpMioqik6MKy62UU++brNz2S/ZoaV1gTWqegfwhYjsnYYyuTQJjy4aPdqacerVg9at4bPP4MADt71J5803oUEDy5sE1mEdGfHknMteCQcDEfk/4G/A0GBXTeDxdBTKpUe8oanXXmtNRgcdBK++ajmNylPRgjp9+lifxJlnQvPmMH9++SOenHPZIZmawUlAT2AdgKp+iQ05dVVE7NDUSPqKKVNg5kz7RX/++fa6LBUNHZ092547drTawfz50cV4nHPZK5lgsEltHKoCiEi99BTJVZZwTWGXXWy0UXGxdfqC/dq/9daSo4v++19o0gQ6d4Z+/Ur3M8yebdc64AALBl99ZdlVnXPZLZlg8IyI3AvsKCLnA68A96WnWK4yxNYULrrI0l9PmwZff22jgQYPhurVrS+gY0cLACtW2HmPP15y6KiqBYMOHaxjOtJvMG9eZX8z51yyEgoGIiLA08CzwARgf+DvqnpnBe/bU0SKRGSZiLwnIqXWPBAzSkRWiMjbIuKzmjOkUye47z5bJKd9exg2DPbZB66+2lJiz50LZ5wB48ZZoADriI40GX36KXzxhQUNsH6Ihg2tqcg5l90qzE0EoKoqIi+q6sHAjCSuXwxcpapLghnLi0VkhqouC51zHPCH4NEWGBM8uww4+2x46CGYM8dqBI0bW16jxYstKHTrZk1Djz9ugeGoo6JNRZ9+ateIBIPq1eGIIzwYOFcVJNNMtEREkpqOpKpfRZbGDNJYLMfyGoWdCDyq5nWsGWq3ZD7HpU5kHsIll1jaitNOg1WrLGX2gw/CU0/Zjf+kk2xdhZkz4YEHrP9h9mxo2tSGqEa0b2/X+/HH+J9X0egk51zlSCYYtAUWiMjHQXPOOyLydqJvFpEWwKHYzOWwZsDnoe0vKB0wEJELRGSRiCxavXp1EsV2iQrPQ7jzTktVMXiwPUdSZj/3XPT8Sy+F9evho4+s1jB7ttUKJLSyRfv21pewYEH8z/TEds5lh2SCQVdgX6ATcAJwfPBcIRGpj/U1XK6qa5ItJICqjlXVPFXNa9q06dZcwlUgdh5CcbEtq1lcbNvhlNlgmVCPPtr6DVassPWZI01EEW3bWnNRWU1FkWv26mXX8tXYnMuMhPoMAFR1FYCI7ALUSfR9IlITCwRPqOpzcU75L7BnaHuPYJ+rZLETw+JNFMvPL3mjbt7cJqsNHmzbHTqUXL+gXj1o1ar8foM2bSzgzJsHV1zhgcC5TEhmBnJPEfkI+BSYA6wEXqrgPQI8ACyPLIITx2TgrGBUUTvgJ1X9KtFyucw6+2yoVs1WUGvSBP73v9LNPO3bWyrrX3+Nf40rrrDmJrDRTJ7HyLnKl0wz0Y1AO+BDVd0b6Ay8XsF72gP9gE4isjR4dBeRC0XkwuCcF4FPgBXYvIWLk/oGLqO6dIH+/e1106bW4Rxu5ikshEaNLEPqf/5j+8IdxFOmWAd027ZWy2jZ0hPbOZcJCTcTAb+q6nciUk1EqqlqkYjcXt4bVHUeIBWco8DAJMrhsszw4Tb3YPlyG3UUbuZp08byFYE1Fa1bF+0XAFtfQRXuvtuGtD7wAEyYUPZSns659EimZvBj0BE8F3hCRO4gyFPkctvbb0OdOhYIYtcvyM+HZ5+1pqQbb4TevaM1h+++s5nNvXvDYYdBz55Wgygu9sR2zlW2ZILBicAG4ApgGvAxCY4mctuv8HDUyPDT2Gae/HxLY/HDD7BmTXQ00i23wNq1tuJaYaF1PjdoYP0PzrnKlcziNutUdbOqFqvqI6o6SlW/S2fhXPaLlxY7dv2CoiJbO+HSS61J6G9/g3PPtSaiTp2sFtCmDdSqBccdB1OnwpYtmfk+zuUq0QTzC4vIWoKMpUAtbD2DdaraME1lK1NeXp4uWrSosj/WbYVwzSE/34JCr17WFCRiqS4mTIgGkyeegL59bfTR4YdnsuTObX9EZLGq5sU7lkzNoIGqNgxu/jsAvYG7U1RGt52KrTn06GFZUffd12oJl1xSsqP4uONskpo3FTlXuZJd9hKwEUCqOhGblexcmWLTZIN1Jv/0U/wO5513tpnIHgycq1wJDy0VkZNDm9WAPOCXlJfIbddim43y80unoOjZE6680hbFadEik6V1LnckUzM4IfToCqzFRhg5l7CKOpwLC23yGtiENPAsps5VhmRyE52bzoK43FBRvqNIFtO99rKmooMOKjlJLZ7CQntf7PrOkfxIzrmKJdNMNKq846p62bYXx+W6SE2hRw/LgrpokaXNLm82ciSARGoc4aYo51xikmkmqgMcBnwUPFphQ0wXBw/nUiI/HwYOtNFGv/4K++1X8fkPPADdu8PFF3sabOe2RjLBoCXQUVXvDNY+7gy0CiagPZKe4rlcVFQEDz8MF15ouYy6dIFfKhiq8PTTds6YMTBggAcC55KVTDDYCQhPMKsf7HMuZcJNPGPGwLHHwvvv20S1yPzI2A7lZ56BJ5+EmjVt+/bbPeupc8lKJhj8C/iPiDwsIo8AS4Cb01Msl6tiRxsNHQq1a8P06XDXXaWXxfzySzjvPKhRA158Ebp1g82bLVOqBwTnEpfMaKKHROQlbC1kgL+p6tfpKZbLVbGjfyIpLLp3t9xGNWpE02SrWiBYv94yn3bpYjObDzoI9t/fMqJ6c5FziUlmpbP2wFpVnQQ0AIaISPO0lcy5QOfOthoaWE6j66+3juK77oKXXoIddrD+BYC997YAsGAB/P730Wv4XAXnypdMM9EYYL2IHAJciaWwfjQtpXIupKjIRgsVFMBOO9noojFjrKZQsyZMnFiyBnD55VaDOP98+PHH0k1LzrnSkgkGxcGqZCcCd6nqXVgNwbm0iV0vYcIEWxSnUyc7fvHFVnMI69IFRo+29RO6dCk91LSwsHR/gtccXK5LJhisFZGhQF/gBRGphqWxdi5t4qWvGDrUUlwXFFjK63gdxQMGWArsxYst11HsUpzhBXgqo+bgAchlPVVN6AH8DmseOjrY3gs4K9H3p/LRunVrdblp1izVJk3sOd52+LzGjVXr11etUUN1xoySx195RbVuXdVu3eK/P1kjRpS8xogRqiNH2nOkPI0aqV5wQfnldi6dgEVaxn01mfUMvlbVW1X11WD7M1X9rc9ARBakMkg5F0+iK6udeiqMHw/332+dziedFP1lrmrH1q+3tRXatdv2UUextY0aNWDwYFi9Gl5+GV55BTZuhEcegUGDfJa0y0JlRYlkH8B/UnWtih5eM3DlCf9K37JF9Q9/UK1dW/W662z7sstUQbVmTdV99rHXN9+87Z8b+fV/9NF27Zo17drxHgUF2/55ziWLVNQMEokrKbyWc1stvKCOCAwbZr/KFyyw9ZdHBSkXhw+HpUuhcWO49lq4557oNbamPf/HH23RnldftWsOGGA1ErCO7okTbfTTLruUXtTHuYwrK0ok+wCWpOpaFT28ZuCS1a9f9Fd5jRqqt9wSPTZ+vKqIap06qh9/vHXt+Vu2qO63n2q1aqpDhtj7R46054IC1YYNrdZw7rn2WRMnep+Bq3xsS81ARGonGFdkW4KSc+k0diw0a2avr7kGrroqeqxPH3jwQas9HHnk1rXn33YbfPih9QeMGGEjngYPtucbboDTT7dQ1LKlPYuU7utwLpMSaSZaACAij1VwXr9tL45z6bFggd3sCwqsOSi2ieaccyylxf/+B3l5yXfs3nuvTYi7OcjWVVwMt9xiz5HjEydaFtbatWHuXPsMX3zHZYtEchPVEpG/AEfGrIMMgKo+Fzy/m+rCOZcKiay7XFQEs2dbOovp0+H556Pt/RVZutRqBTffDHXq2L7yVnR7+WULBs5lk0RqBhcCRwM7UnId5BOA49NWMudSpKLhqOFg8cIL1oTzl78k3sH7739D/frR/EgV6dABliyBtWuT/y7OpUuFwUBV56nqRcAQVT035vHXSiijc9skPLooItxEEw4WBx5oSfF++cWW24wnPJt45UpbWOe44+C++xIrzzHHWJrt117bqq/jXFokM7T0MRG5TESeDR6Xioino3BVXmywaNjQ2v/ffBO2bLF94aGm4QlmI0favpkzE09nccQRNiltzpzUfQfntlXC6xkAd2O5iO4OtvthmUzPS3WhnMuko4+2zt8337TlN994A556yjqAIZofqUcPW6O5Rg149tnEO53r1bNOau83cNkkmWDQRlUPCW3PEpG3Ul0g5zItP99u/F27Qv/+UK2a3fAXLLBaw6hRliBv82Y7//zzkx99dMwxNhx1/XqoWzflX8G5pCXTTLRZRPaNbIjIPsDm1BfJuczr1Akuu8xe77WXzQ247jr7Rf/oo7avXj2bVzB+fPKziTt0sFrFG2+kvuzObY1kgsHVQJGIzBaROcAs4KoK3uNclVRUZDf9ggL4+WeYMgX69rVjffrYSKApU+D2263zOZykLhHt21uNw/sNXLZIZg3kmSLyB2D/YNcHqroxclxEjlXVGakuoHOVLd68hF69bMhpQYF1Gt94Y/yhqok2FzVqBK1aeb+Byx5JJapT1Y2q+nbw2BhzeEQKy+VcxsTOSwALBKedZqklpk61JHfhmkCys4kLC22CW2RmNPhiNy6zkulArojnJnLbhdib+sKFNiN5W2oCsdq0gZtusvkMixbBpk3R2ohzmSCWyC4FFxJZoqqHpeRiFcjLy9NFixZVxkc5lzYTJ1rKi44d4d13fbEbl34islhV8+IdS+V6Bs65JPTqBc2bW06ks87yQOAyK5XBYGXsDhF5UES+EZG4SexEpJGITBGRt0TkPRE5N4XlcS6rFRXBmjU2quiuu7ZtsZtwiozw9b0PwiUq4WAgItVFpGeQkuLKyCNyXFVLZTQFHga6lXPZgcCyYDJbR2CkiNRKtEzOVVWREUsTJtj8hY0bS67TnKzYNZgj1080RYZzyXQgTwF+Ad4BtiTyBlWdKyItyjsFaCAiAtQHvgeKkyiTc1VSeMTSEUfY659+gvnzt665KNKp3bOnpdHesiW5FBnOJRMM9lDVlin+/NHAZOBLoAFwmqrGDTQicgFwAcBee+2V4mI4V7nCI5bq1LFO5HvvtYAA1rxTo4YtjhM5t6jIgkhZQ1jffNMmyP38M7RrVzoQFBZaTSG8v6JrutyRTJ/BSyLy5xR/fldgKbA70AoYLSIN452oqmNVNU9V85o2bZriYjiXWaedZkFh5EiYPNlWXBs8GFavhu+/hwEDrMM50uxTWAi33mrPqnYzv+YaqF4dDjjA0ly8+GLJz8iGpiTv28hiZS2OHPsATgLWARuANcBaYE0C72sBvFvGsReAo0Pbs4DDK7pm69att3FZaOeyz+TJqiKqdnuPPqpVU61ZU7V+fdVZs+zckSPt3H/9S7V//+i5t9yiumCBvQ6fHzFtmmrduqqXXabapEnp4+k2a5Z97syZJbcruxy5ClikZd2ryzpQ6kT4FGhJMDchifeVFwzGANcHr3cF/gs0qeiaHgzc9uryy+3/ylNOUZ0wQbV3b9uuWdOea9RQPfJI1Vq17CYaCQI1a1ogiGjbVrVZMwsWYddcE33PsGGV+90ibrnFAtl553kgqGypCgZzgWqJnh+8ZxzwFfAr8AXQH1tG88Lg+O7Ay1in9LtA30Su68HAbY8iv5ILCux55MjoduPGqoMGRQNA48aqJ5+sevTRtl1QUPJa48bZ/ilTovvmz7ebcLVqdqxBg8q/EX/6qZU9EpBiy+3SK1XB4OEgIAwFrow8En1/Kh8eDNz2Jra5JNIMNHJk9HjDhqqNGqlee23pYBH7C3vTJtU99lDt3Nm2f/5ZdffdLRBMmqR6yCH2nsaNKy8grFun2qqVar169t0OPNBrBpWtvGCQTAfyp8BMoBY28ifycM5to9jkeMXFttpacWigdSRZ3k032Uprgwfb8w03lE6jXbMmHHywLcf59tt23pdfWkf0++/bJLdvv4Vjj7XPTodwZ7GqLQK0dKl9p2OOgQ0bti79t0uTsqJENj+8ZuByzYgRJX9BjxhhNYMRI6L7Zs0quT1xotX999vPnk8+ueQv8datVatXV33//fjXiP3MeJ9RnnBtZ+RI/a1vI1JuUF29Orlrum1DipqJirDRPiUeib4/lQ8PBs4l5oQT9Lc+htgmmWefteaa1q1Vt2wp3VRV0XYiZs1S3XFHK0OtWtFRRDNn2r5p01LzPV1iygsGyUw6Gxx6XQfojc8Wdi6r3X67Nc18/rktzBOecNa7NwwcCKNHWyqM+fNLNlVFZjX37g39+sGTTyafWfWoo6zJCuDyy205UYDWre150SJba9plXsJ9Bqq6OPSYr6pXYvmEnHNZatUqa5svKIAxY0q3zd92G+y5J0yaBF26lL7RFxfbEp+jRlmbf7LpLS67zCbOnXIKPPhg9PMbNYL99rNg4LJDMonqdg49mohIN6BRGsvmnNsG4eU743UyA7z6KqxfD7vtBk89BRddFD02aRJ06wabN9v2qFHJdfQ+8wzccw+0bQtPP1368/PyPBhkk2RGEy0GFgWP17Chpf3TUSjn3LaLHaEUXqENosFi/Hh45x3YdVe7eV96qTUJnXSSJby74Qbo0cNGBJ1ySuIB4aabrInoySdtJFTs5+flwRdfwNdfp/67u61QVmdC5AG0AX4X2j4bSy43Cti5oven4+EdyM5tu9jRQl99VXJWc7VqqjffbMeWLbORRz17lj3yJ3y9GTPsGuecU/b5c+faOVOnpu47ufKxjfMM7gU2AYjIMcBw4BHgJ2BsyqOTc65SDBlSsg/gd7+DJUtgp52ix4cOtdcHHggXXghTp1otIZ5IIrzp061jevfd7fyyEuEdeqjVGLypKDskEgyqq+r3wevTgLGqOkFVC4Dfp69ozrnKtmKFZT4tKID77y/ZJNS4MeywA1x9dXRfOONofr71DfTqBR9+aKm0yxt9VL++BRkPBtkhoWAgIpEhqJ2x+QURyQxNdc5lsYo6nDt2tOeXXrJf//FSYC9ZAr/8Yq8HDap49FGkE9laoF0mJRIMxgFzRGQSlr76VQAR+T3WVOSc2w5U1OGcn2/LdFarBn36wHHHWTNS5Pxx46zWUL06DBsWfyhrrLw860D+8sv0fS+XGNEEQrKItAN2A15W1XXBvv2A+qq6JL1FLC0vL08Xed3SuYwYMADGBr2FtWrBY49BkyaW52jLFhgxwvobwjWNsmoICxbAkUfCxIlw4omV9hVylogsVtW8eMcSGlqqqq+r6vORQBDs+zATgcA5lzlFRfDcc3DttVCvnt38TzvNagmq8M9/RpfQjK1ZxHPIIVaTiPy285XQMieZeQbOuRwW/qV/000wZQo0bGijhjZtstnG111X8j35+eWvr1y3LvzpT9FgEBmR9PTT8Zf7jJTDg0PqeTBwziUkXp/CddfBTz9ZH8ETTySfirqwEJo1i3YiH3QQ7LsvnH66jV66/37YuDGa+joT6zbnjLImIGTzwyedOZd5qcpqWr++/rbqWb169vqAA+x5jz2iK7P5YjjbjhQtbuOcc7+paPRRIvLzo00+N95ozU1DhtjCOwUFNkx14kRo3x6WL4cOHZJLlud9EEkoK0pk88NrBs5tP379VbVdO/v137dv6dpGZLnP3Xaz9Reefz7xa6ei9rI9wWsGzrls9eqrNvO5oMBGKoXnLkB0uc9p0+z16adHf+1X9Ms/Ulvp08eGrlY01DWXeTBwzmVM7KznqVNh+PDoDX7hQnj+ebj3XmjZ0m7iGzfa3AawjuSTTrJRRxB/9NGSJbYmw+TJ9t7YCW7ebGQ8GDjnMqaifofYZHqRGc5Tp9rQ07fessV77r8fdtkFHn00eq4qXHABDB5sazJ062ZBoW/f6BBYH50UldAM5GzjM5Cdy12jR9uaCyJ2w69WzVZrW7XKjtesafuaN7eEeTVrWid09+6WMqNvX5ss17Kl1RJim40KCy04hPcVFVmAKm/ORFWwzTOQnXMuW1xyif3K12CxnUmTYN0663PYaSc7tmWLBYI//MH6Grp3t/eecYYl2WveHN5+20YpxfYfRCa+TZtmnxFbe9heRyh5MHDOVSlFRTZJraDAbux9+0b7HCZMgDlzLNX2FVfADz9YDSKsenULHk2bWnPTrFklj+fnwx13WIqNnj1LdzpHgkUkIGw3waKsYUbZ/PChpc7lptihoRdcYENPw0NHGzWy/fHOD2/ff78NZ23UqPRQ07PP1t9WfBs2LH45dtpJ9corSw9VDX/Gli3ZNZyVcoaWZvzGvjUPDwbO5abYpTpVbTuytGYyxzdsUG3aVLVt25JLcy5bZvMZqlcvO1j89JNqnTp2fNCg0uWcPl21Vi2bQZ0tgUDVg4FzzsV1/fV2F1y+PLrvmGNs35NPWq2jc+fSN/T+/e2cmjUtaIQnwm3apNqrV7Rm0adPpX2dCpUXDLzPwDmXsy66CGrXtj4CsDkJc+fCWWdZZ/NZZ9mkuPvuiw53feEFePBBaNsWZsywW37v3vDii5ZM7y9/sdFLtWpZh/bEiaX7JbKRBwPnXM7aZRfrgH7kEcuHNGwY7Lwz3HmnHR8wwPIlrVgRHVZ6110WAO64w3IlXX+9jV4aMAD69YNnn7Xzhg+Hf/3LAsTJJyef0bWyeTBwzuW0Ro1s4tpf/2rrO19zDSxebKN/DjoIjjrKZkBv2WKjkBYuhK5drWYANqpp8GD44gt46imb1zByJFx5pdUsfvc7S8udTAK/TPBg4JzLaccfbzfwKVNgt90sAISHil54odUMiorgnnusBvH3v5e8xr//bcNQwYLJlVfa6zp17PWSJdCpU+JlysTwVA8Gzrmclp9vK7cBtGplv+bD8wpWrrQV3W67zW76nTtbjqPwjbmoCF57zWoJY8aUvJFv2GBLhA4fXvL88m7sFc1lSIuyepaz+eGjiZxzqTZggP62yE7YrFmqO+wQHR10223lp8WOt123bnTUUqLzDmbNUt1xR9VrrrH3jxxZ+nh4SGwi8KGlzjlXtsgNuqAg/o36scfsbtm8eenjFc1tUFV97jl7/yGHJD7vYPr0aABq1Uq1ceNtX5fBg4FzzpUh0QVweveOX3NI1OGH2/tPPbXic7dsUf3Tn2zJz3331d/mNOywg01y29qJbOUFA+8zcM7ltESW7ywqspxH8foEElFUBB9/DLvvbtceNar882+7Dd57z7Kzrlhh8xyqVbP+hzvusPkRKV+gp6wokc0Prxk45yrLti6dGT5/9WpLbwGqDzxQ8pxws9K++6ruvLOlzIgYOdJqBkOGeM3AOecqXSI1h0Tf36SJ1SxE7Ff/Z5+VHin02mtWixg2zIamgp0zfLjNfh4xwq4XHm2UCmld3EZEHgSOB75R1YPKOKcjcDtQE/hWVTtUdF1f3MY5V5Xdf7+twrbzzjaZbcKEaLDp3t0CyMqVNiQVUrfgTnmL26Q7GBwD/Aw8Gi8YiMiOwGtAN1X9TER2UdVvKrquBwPnXFV3zjmWBqNaNbj8crjlFpv53KYNnHeeLcyT6pXVMrbSmarOBb4v55S/AM+p6mfB+RUGAuecq+qKiqzJZ9AgazK69VZLe3HddVYbeP75yl+XOdN9BvsBO4nIbBFZLCJnlXWiiFwgIotEZNHq1asrsYjOOZc6kT6CZ56B22+31drq1LG+gpdftpkF48enYbRQBTIdDGoArYEeQFegQET2i3eiqo5V1TxVzWvatGllltE551ImtkO6c2dLf3344bY9cGDlBwLIfDD4ApiuqutU9VtgLnBIhsvknHNpM2RI/Jv9J5/YPIaHHspMuutMB4NJwFEiUkNE6gJtgeUZLpNzzlWacLPRDTekZ9hoImqk8+IiMg7oCDQRkS+A/8OGkKKq96jqchGZBrwNbAHuV9V301km55zLJuXNY6jM5qK0Di1NFx9a6pxzycvY0FLnnHNVgwcD55xzHgycc855MHDOOYcHA+ecc1TR0UQishpYtZVvbwJ8m8LipENVKCNUjXJ6GVPDy5gamS5jc1WNm8KhSgaDbSEii8oaWpUtqkIZoWqU08uYGl7G1MjmMnozkXPOOQ8GzjnncjMYjM10ARJQFcoIVaOcXsbU8DKmRtaWMef6DJxzzpWWizUD55xzMTwYOOecy61gICLdROQDEVkhItdkujwAIvKgiHwjIu+G9u0sIjNE5KPgeacMl3FPESkSkWUi8p6IDMq2copIHRF5U0TeCsr4j2D/3iLyRvBv/rSI1MpUGUNlrS4i/xGRqdlYRhFZKSLviMhSEVkU7Muaf+tQOXcUkWdF5H0RWS4iR2RTOUVk/+BvGHmsEZHLs6mMYTkTDESkOnAXcBzwR+AMEfljZksFwMNAt5h91wAzVfUPwMxgO5OKgatU9Y9AO2Bg8LfLpnJuBDqp6iFAK6CbiLQDRgC3qervgR+A/pkr4m8GUXIRp2wsY76qtgqNic+mf+uIO4BpqnoAtkLicrKonKr6QfA3bIUt77seeD6byliCqubEAzgCW2Izsj0UGJrpcgVlaQG8G9r+ANgteL0b8EGmyxhT3knAsdlaTqAusARbOe9boEa8/wYyVLY9sBtAJ2AqIFlYxpVAk5h9WfVvDTQCPiUYBJOt5QyV68/A/GwuY87UDIBmwOeh7S+CfdloV1X9Knj9NbBrJgsTJiItgEOBN8iycgbNL0uBb4AZwMfAj6paHJySDf/mtwNDsJX9ABqTfWVU4GURWSwiFwT7surfGtgbWA08FDS53S8i9ci+ckacDowLXmdlGXMpGFRJaj8fsmL8r4jUByYAl6vqmvCxbCinqm5Wq5LvARwOHJDJ8sQSkeOBb1R1cabLUoGjVPUwrEl1oIgcEz6YDf/W2JK9hwFjVPVQYB0xzS1ZUk6CPqCewPjYY9lSRsitYPBfYM/Q9h7Bvmz0PxHZDSB4/ibD5UFEamKB4AlVfS7YnXXlBFDVH4EirMllRxGJrPWd6X/z9kBPEVkJPIU1Fd1BdpURVf1v8PwN1sZ9ONn3b/0F8IWqvhFsP4sFh2wrJ1hQXaKq/wu2s7GMORUMFgJ/CEZu1MKqbZMzXKayTAbODl6fjbXRZ4yICPAAsFxVbw0dyppyikhTEdkxeL0D1qexHAsKfYLTMlpGVR2qqnuoagvsv79ZqnomWVRGEaknIg0ir7G27nfJon9rAFX9GvhcRPYPdnUGlpFl5QycQbSJCLKzjLnTgRx01nQHPsTakq/LdHmCMo0DvgJ+xX7t9MfakWcCHwGvADtnuIxHYVXZt4GlwaN7NpUTaAn8Jyjju8Dfg/37AG8CK7Bqeu1M/5sH5eoITM22MgZleSt4vBf5/ySb/q1DZW0FLAr+zScCO2VbOYF6wHdAo9C+rCpj5OHpKJxzzuVUM5FzzrkyeDBwzjnnwcA555wHA+ecc3gwcM45hwcD51JGRFqEs886V5V4MHDOOefBwLl0EJF9ggRqbTJdFucSUaPiU5xzyQhSJDwFnKOqb2W6PM4lwoOBc6nVFMs1c7KqLst0YZxLlDcTOZdaPwGfYfmcnKsyvGbgXGptAk4CpovIz6r6ZKYL5FwiPBg4l2Kqui5YyGZGEBCyNVW6c7/xrKXOOee8z8A555wHA+ecc3gwcM45hwcD55xzeDBwzjmHBwPnnHN4MHDOOQf8P8gQTMrZ95orAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(amount_clusters)\n",
    "plt.plot(np.arange(len(wce)), wce, 'bx-')\n",
    "plt.xlabel('k')\n",
    "plt.ylabel('Sum_of_squared_distances')\n",
    "plt.title('Elbow Method For Optimal k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#KMedoids Sklearn (custom distance metric)\n",
    "from sklearn_extra.cluster import KMedoids\n",
    "kmedoids = KMedoids(n_clusters=2, random_state=0).fit(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov Chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MakeTransitionMatrix(history, NumComponents, normalize=True):\n",
    "#      Input: DataArray - array of length N of integers that belong to 1 of C classes\n",
    "#      Output: TransMat - is a right stochastic matrix (each row sums to 1)\n",
    "#                         Each entry of TransMat(i,j) represents the probability\n",
    "#                         of going from label i to label j\n",
    "    NumTimeSteps  = len(history)\n",
    "    \n",
    "    #Pre-allocate transition matrix\n",
    "    TransMat = np.ones((NumComponents+1, NumComponents+1)) #+1 because 0-indexed\n",
    "\n",
    "    for i in range(NumTimeSteps-1):\n",
    "        Current = history[i]\n",
    "        Next    = history[i+1]\n",
    "        TransMat[Current, Next] = TransMat[Current, Next] + 1\n",
    "        \n",
    "    #Normalize rows\n",
    "    if normalize:\n",
    "        TransMat = TransMat/np.sum(TransMat,1)[:,None] #broadcast vector for normalization by rows\n",
    "    return TransMat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FindNextState(TransMat, state, loc_dict, plotNextStateProb=False):\n",
    "    #assume state is 0-indexed\n",
    "    rand = np.random.random()\n",
    "    \n",
    "    #Make cdf of transitioning out of a given state\n",
    "    cdf = np.cumsum(TransMat[state])    \n",
    "\n",
    "    #find first state above random number on cdf\n",
    "    nextState = np.where(cdf > rand)[0][0]\n",
    "    if plotNextStateProb: #plot cdf and probability of transitioning\n",
    "        StateSize = len(TransMat)\n",
    "        plt.figure(dpi=100)\n",
    "        plt.plot(np.insert(cdf, 0, 0))\n",
    "        plt.title(\"Cumulative Distribution function\" + \" - Chose state \" + str(nextState) + ' ' + loc_dict[nextState])\n",
    "        plt.axhline(rand, c='red')\n",
    "        for i in range(StateSize):\n",
    "            plt.axvline(i, c='k')\n",
    "        plt.show()\n",
    "    \n",
    "    return nextState"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PredictFutureSteps(TransMat, MostRecentState, bldg_dict, NumFutureSteps, PlotCDF=False):\n",
    "    nextStates  = []\n",
    "    nextBldgs   = []\n",
    "    currState   = MostRecentState\n",
    "    for i in range(NumFutureSteps):\n",
    "        nextState = FindNextState(TransMat, currState, bldg_dict, PlotCDF)\n",
    "        #print(\"\\nNext state\", bldg_dict[nextState])\n",
    "        nextStates.append(nextState)\n",
    "        nextBldgs.append(bldg_dict[nextState])\n",
    "        currState = nextState\n",
    "\n",
    "    #print(\"Next states array:\", nextStates)\n",
    "    #print(\"Next buildings array:\", nextBldgs)\n",
    "    return nextStates, nextBldgs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Toy example for Markov Chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TransMat:\n",
      " [[0.375      0.375      0.125      0.125     ]\n",
      " [0.22222222 0.22222222 0.44444444 0.11111111]\n",
      " [0.2        0.2        0.2        0.4       ]\n",
      " [0.14285714 0.28571429 0.42857143 0.14285714]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([0, 0, 3, 0, 1],\n",
       " ['Not Connected', 'Not Connected', 'lib west', 'Not Connected', 'hub'])"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# demo of a Markov chain predictor\n",
    "\n",
    "#Toy data\n",
    "bldg_dict = {0: 'Not Connected', \n",
    "             1: 'hub',\n",
    "             2: 'reitz',\n",
    "             3: 'lib west'}\n",
    "\n",
    "personHistory = np.array([0,0,0,1,1,2,2,1,2,3,2,3,2,3,1,2,0,1,0])\n",
    "\n",
    "#Make transition matrix\n",
    "TransMat = MakeTransitionMatrix(personHistory,3)\n",
    "print(\"TransMat:\\n\", TransMat)\n",
    "\n",
    "#Predict future steps based on their most recent history\n",
    "NumFutureSteps  = 5\n",
    "MostRecentState = personHistory[-1]\n",
    "nextSteps, nextBldgs = PredictFutureSteps(TransMat, MostRecentState, bldg_dict, NumFutureSteps, PlotCDF=False)\n",
    "nextSteps, nextBldgs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MC prediction on each cluster\n",
    "Make Transition Matrix on entire cluster, leave 1 person out. \n",
    "Use remaining users in cluster to build transition matrix.\n",
    "User left-out user to predict 1 step ahead, compare to ground truth. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "bldgs     = location_info['name'].unique()\n",
    "bldg_dict = dict(zip(np.arange(len(bldgs)), bldgs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 'Academic Advising Center',\n",
       " 1: 'Martin H. Levin Advocacy Center',\n",
       " 2: 'Mec. Aerospace Engineering',\n",
       " 3: 'Alfred A. Ring Tennis Pavilion',\n",
       " 4: 'Anderson Hall',\n",
       " 5: 'Animal Science',\n",
       " 6: 'Animal Science Building C',\n",
       " 7: 'Aquatic Products Lab',\n",
       " 8: 'School of Architecture',\n",
       " 9: 'University Auditorium',\n",
       " 10: 'Bartram Hall',\n",
       " 11: 'BioContainment Facility',\n",
       " 12: 'Bio-Containment Facility',\n",
       " 13: 'Tolbert Hall',\n",
       " 14: 'Bruton-Geer Hall',\n",
       " 15: 'Black Hall',\n",
       " 16: 'Broward Hall',\n",
       " 17: 'Bryant Space Science Center',\n",
       " 18: 'Bryan Hall',\n",
       " 19: 'UF bookstore',\n",
       " 20: 'UF Bookstore',\n",
       " 21: 'UF Welcome Center',\n",
       " 22: 'Carleton Auditorium',\n",
       " 23: 'Carleton Auditorium - Attic',\n",
       " 24: 'Classroom Building 105',\n",
       " 25: 'Centrex Building',\n",
       " 26: 'Dept. of Chemical Engineering',\n",
       " 27: 'Chemistry Lab',\n",
       " 28: 'Constans Theatre',\n",
       " 29: 'Center for Outdoor Rec. & Educ.',\n",
       " 30: 'Criser Hall',\n",
       " 31: 'Carr Hall',\n",
       " 32: 'CISE',\n",
       " 33: 'Crops And Soils Field Lab',\n",
       " 34: 'Lacy R. Chilled Water Plt. CWP2',\n",
       " 35: 'Dauer Hall',\n",
       " 36: 'Fresh Food',\n",
       " 37: 'Engineering Design Building',\n",
       " 38: 'Env. Health and Safety EHS',\n",
       " 39: 'Elmore Hall',\n",
       " 40: 'Entomology-Nematology',\n",
       " 41: 'Fine Arts A',\n",
       " 42: 'Fine Arts B',\n",
       " 43: 'Fine Arts C',\n",
       " 44: 'Fine Arts D',\n",
       " 45: 'Forestry Resource and Conserv. FF1',\n",
       " 46: 'Forestry Resource and Conservation',\n",
       " 47: 'Fifield Hall',\n",
       " 48: 'Florida Gym',\n",
       " 49: 'Keene-Flint Hall',\n",
       " 50: 'Innovation Hub',\n",
       " 51: 'Griffin-Floyd Hall',\n",
       " 52: 'Florida Pool Offices/Lockers',\n",
       " 53: 'Gator Corner Dining Center',\n",
       " 54: 'Gerson Hall',\n",
       " 55: 'Grinter Hall',\n",
       " 56: 'Harn Museum Of Art',\n",
       " 57: 'Heritage Hall',\n",
       " 58: 'Hough Hall Business School',\n",
       " 59: 'Holland Law Center',\n",
       " 60: 'HR Building',\n",
       " 61: 'The Hub',\n",
       " 62: 'Hume Commons',\n",
       " 63: 'Institute Of Black Culture',\n",
       " 64: 'Hispanic/Latino Cultures Inst',\n",
       " 65: 'Mycology Lab',\n",
       " 66: 'Office Building',\n",
       " 67: 'Infirmary',\n",
       " 68: 'Jennings Hall Office',\n",
       " 69: 'Larsen Hall',\n",
       " 70: 'Library West',\n",
       " 71: 'Little Hall',\n",
       " 72: 'Materials Engineering',\n",
       " 73: 'MAE C',\n",
       " 74: 'Matherly Hall',\n",
       " 75: 'Mccarty Hall A',\n",
       " 76: 'Mccarty Hall B',\n",
       " 77: 'Mccarty Hall C',\n",
       " 78: 'Mccarty Hall D',\n",
       " 79: 'Mcguire Center',\n",
       " 80: 'Microkelvin Laboratory',\n",
       " 81: 'Microbiology/Cell Science',\n",
       " 82: 'Harrell Medical Education Building',\n",
       " 83: 'Mehrhof Hall',\n",
       " 84: 'Marston Science Library',\n",
       " 85: 'Music Building',\n",
       " 86: 'Dickinson Hall',\n",
       " 87: 'Nanoscale Research',\n",
       " 88: 'New Engineering Building',\n",
       " 89: 'Physics',\n",
       " 90: 'Norman Gym',\n",
       " 91: 'Norman Hall',\n",
       " 92: 'Norman Hall Addition',\n",
       " 93: 'Nuclear Sciences',\n",
       " 94: 'Newins-Ziegler Hall',\n",
       " 95: 'Ornamental Hort Laboratory 1',\n",
       " 96: 'Ornamental Hort Laboratory 2',\n",
       " 97: 'Phillips Center',\n",
       " 98: 'Physical Plant Central Stores',\n",
       " 99: 'Phelps Lab',\n",
       " 100: 'Agronomy Plant Science Lab',\n",
       " 101: 'Powell Hall',\n",
       " 102: 'Presidents Residence',\n",
       " 103: 'Plant Science Facility',\n",
       " 104: 'Particle Science & Technology',\n",
       " 105: 'Psychology',\n",
       " 106: 'Pugh Hall',\n",
       " 107: 'Racquet Club Dining Center',\n",
       " 108: 'Reitz Student Union',\n",
       " 109: 'Student Recreation and Fit. Center',\n",
       " 110: 'Rhines Hall',\n",
       " 111: 'Research Laboratory 1',\n",
       " 112: 'Research Laboratory 2',\n",
       " 113: 'Reed Lab',\n",
       " 114: 'Rinker Hall',\n",
       " 115: 'Rogers Frazier Hall',\n",
       " 116: 'Rolfs Hall',\n",
       " 117: 'Steinbrenner',\n",
       " 118: 'Sisler Hall Chemistry Research Bldg',\n",
       " 119: 'Smathers Library',\n",
       " 120: 'OConnell Center',\n",
       " 121: 'Ben Hill Griffin Stadium',\n",
       " 122: 'Stuzin Hall',\n",
       " 123: '# Stuzin Hall Walk Way',\n",
       " 124: 'SW Recreation Center',\n",
       " 125: 'Turfgrass Envirotron',\n",
       " 126: 'Tigert Hall',\n",
       " 127: 'Turlington Hall',\n",
       " 128: 'Urban Entomology',\n",
       " 129: 'Ustler Hall',\n",
       " 130: 'Van Fleet Hall',\n",
       " 131: 'Walker Hall',\n",
       " 132: 'Weil Hall',\n",
       " 133: 'Clock Tower',\n",
       " 134: 'Wiemer Hall',\n",
       " 135: 'Williamson Hall',\n",
       " 136: 'Yon Hall'}"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bldg_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "NumClusters = len(roundedMeans)\n",
    "nBuildings = 137"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster: 0\n",
      "Average MC(O=1) prediction accuracy for cluster: 0 is 1.42\n",
      "Cluster: 1\n",
      "Average MC(O=1) prediction accuracy for cluster: 1 is 2.69\n",
      "Cluster: 2\n",
      "Average MC(O=1) prediction accuracy for cluster: 2 is 3.28\n",
      "Cluster: 3\n",
      "Average MC(O=1) prediction accuracy for cluster: 3 is 1.1\n",
      "Cluster: 4\n",
      "Average MC(O=1) prediction accuracy for cluster: 4 is 1.06\n",
      "Cluster: 5\n",
      "Average MC(O=1) prediction accuracy for cluster: 5 is 3.41\n",
      "Cluster: 6\n",
      "Average MC(O=1) prediction accuracy for cluster: 6 is 1.16\n",
      "Cluster: 7\n",
      "Average MC(O=1) prediction accuracy for cluster: 7 is 0.72\n",
      "Cluster: 8\n",
      "Average MC(O=1) prediction accuracy for cluster: 8 is 0.75\n",
      "Cluster: 9\n",
      "Average MC(O=1) prediction accuracy for cluster: 9 is 1.89\n",
      "Cluster: 10\n",
      "Average MC(O=1) prediction accuracy for cluster: 10 is 1.49\n",
      "Cluster: 11\n",
      "Average MC(O=1) prediction accuracy for cluster: 11 is 2.71\n",
      "Cluster: 12\n",
      "Average MC(O=1) prediction accuracy for cluster: 12 is 1.33\n",
      "Cluster: 13\n",
      "Average MC(O=1) prediction accuracy for cluster: 13 is 1.97\n",
      "Cluster: 14\n",
      "Average MC(O=1) prediction accuracy for cluster: 14 is 1.65\n",
      "Cluster: 15\n",
      "Average MC(O=1) prediction accuracy for cluster: 15 is 8.53\n",
      "Cluster: 16\n",
      "Average MC(O=1) prediction accuracy for cluster: 16 is 0.8\n",
      "Cluster: 17\n",
      "Average MC(O=1) prediction accuracy for cluster: 17 is 1.84\n",
      "Cluster: 18\n",
      "Average MC(O=1) prediction accuracy for cluster: 18 is 1.74\n"
     ]
    }
   ],
   "source": [
    "#Get the average accuracy for each cluster\n",
    "AvgAccPerCluster = []\n",
    "for clusterIdx in range(NumClusters):\n",
    "    print(\"Cluster:\", clusterIdx)\n",
    "    #find all people within a cluster\n",
    "    AccPerPerson = []\n",
    "    UsersInCluster    = clusters[clusterIdx]\n",
    "    NumUsersInCluster = len(UsersInCluster)\n",
    "    for u in range(NumUsersInCluster):\n",
    "        TempUsersInCluster = UsersInCluster.copy()\n",
    "        UserToPredict      = TempUsersInCluster[u]\n",
    "        TempUsersInCluster.remove(UserToPredict) #remove person to use prediction\n",
    "        ClusterSubsetForTransMat = TempUsersInCluster\n",
    "\n",
    "        #Grab the data associated with a person's cluster assignment\n",
    "        ClusterHistory = X[ClusterSubsetForTransMat].astype(int)\n",
    "        NumUsersInSubCluster = len(ClusterSubsetForTransMat)\n",
    "        TransMat = np.zeros((nBuildings, nBuildings))\n",
    "        for i in range(NumUsersInSubCluster):\n",
    "            userHistory = ClusterHistory[i]\n",
    "            TransMat += MakeTransitionMatrix(userHistory, nBuildings-1, False)\n",
    "        TransMat = TransMat/np.sum(TransMat,1)[:,None] #normalize\n",
    "\n",
    "        #Predict each successive step for the left-out user\n",
    "        NumFutureSteps = 1\n",
    "        NumCorrect = 0\n",
    "        userHistory = X[UserToPredict].astype(int)\n",
    "        for t in range(len(userHistory)-1):  \n",
    "            MostRecentState = userHistory[i]\n",
    "            nextSteps, nextBldgs = PredictFutureSteps(TransMat, MostRecentState, bldg_dict, NumFutureSteps, PlotCDF=False)\n",
    "\n",
    "            if nextSteps[0] == userHistory[t+1]: #Does it match the ground truth?\n",
    "                NumCorrect += 1\n",
    "        UserAccuracy = NumCorrect/(len(userHistory)-1)*100\n",
    "        AccPerPerson.append(UserAccuracy)\n",
    "    AvgAcc = round(np.array(AccPerPerson).mean(),2)\n",
    "    print(\"Average MC(O=1) prediction accuracy for cluster: \" + str(clusterIdx) + \" is \" + str(AvgAcc))\n",
    "    AvgAccPerCluster.append(AvgAcc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Baseline: 0.73\n"
     ]
    }
   ],
   "source": [
    "BaselineAcc = round(1/137 *100, 3)\n",
    "print('Baseline:', BaselineAcc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAHPCAYAAACP7aS6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0WklEQVR4nO3deZhkVX3/8feHAUbFAdxFUVFxiUtE3Jcg+SlR4hL3DRPRxCjRJBpNBFeMScBdEUJMXDHGfUWCuMR9QUEiKLgAsqMgIBCEQeH7++Pc1qKo7umuvtO3a3i/nuc+PXXuufd863ZN17fOOfdUqgpJkiT1Z7OhA5AkSdrUmGBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzEyxJkqSemWBJWpQkeyapJDuMlH0pyZd6bGPfJH69hDYoyb8m+dzQcWxIkg8k+dDQcWjlmWBpVUnyV92b+JFDx7LaJDmluzZz2zlJvprkMUPHthRJrtMlUrsOHcskSV7bXd8PDh3LLBlJwCvJAyfsT5LTu/2fnrD/WklekOTIJBcmuSzJj5McmOT2Y3VvDfwF8C8TznODJK9L8qPuHOcnOSLJI3p+vi9N8qkkP++e077zVH0N8Lgkd+uzfa1+JlhabfYATgHunWTHgWNZjf4X+NNuez1wM+BjSZ4zUDx/1G1LcR3glcCuE/b9E3DtZcY0tSQBnkJ7DT4yybqhYplhlwFPnVD+IGB7YP34jiQ3BL4GvBE4B3gF8FzgE8CjgO+PHfK3wE+r6otj57kD8D3gb4AvAs+jJWE3Bg5N8rppn9QE/wTcCzhmoUpVdQxwFPDCHtvWDDDB0qrRfSq9P/B3wLm0ZGulY9gsybVWut0lOLOq/rPbXgs8ALgEeMF8ByTZPMmWGyOYqrq8qi7v8Xy/qarL+jrfFHalJQHPBDYHHjtgLAtKcp2hY5jHfwNPSLL5WPlTgaOBn0045t3A3YHHV9Ujq+otVfWOqvoH4HbAAXMVk2xB+9twlWG3rvwjwPWAXarqOVX19qp6PXBP4IPAi5I8qY8nCdy6qrYDnraIuh8CHpvkuj21rRlggqXVZA/gAuAw2h/K3yZYSbbouvrfNX5Qkq27oYDXj5StTfKqJCcmWd8NTbw2ydqxY6sbgtgjyQ9on64f1u17UZJvJDkvyaVJjk7y+AntXzvJAUl+keTibtjg5pOGDbryd3bDCuuT/CDJM6e9YFX1M+AE4Nbd+Xfo2n1RkucnOal7Tnfq9t8xyUe6a3lZkqOSPGrCc7pzkv/pnvcZSV7GhL8XmTAHqxvq2bcb3rksydlJPpbktmnzt87tqr5yZEhp3+7Yq8zBSvL9JFfppejKN0tyZpKPjJU9v7uml3XX+G1JrreES7oHcHzXM/J55knyu9/jO5Kc1f0ef5rk4NFENsm2Sd6UNrS7vruOh3S9NRPntHXlu3blu46Ufam7FvdI8pUkv6IbHkvyJ0kOG4nlpCQvT7JmQtz3SfLfSS5IckmSY5P8bbfvGV27d59w3EuSXJHk5ou4hu8HbgDsNnL8lsDjgf+aFBPwcOAdVfXR8f1Vtb6qXjRS9EDghrTfz6jHAXcB9q+qq0wxqKorgGcDvwT2XcRz2KCqOmUJ1T8HbMXINdGmb/wThjSkPYCPVdXlSd4P7JXkXlX1nar6dZKP0z4FPnus1+TRwFrgA9DeaIFP0f4Q/zstAbkrrZfn9l39Uf8PeCJwIPAL2vAQtGGITwHvA7YEngx8OMkjquqwkePf3R3/XuBbtKGQ0f10cd2k219dW+cCuwPvSLJ1Vb15sRdq5JxbALcAzhvb9QzgWrTnvx44P8mdga8DZwL703q+ngh8Isnjqurj3TlvShte2Xyk3l8Cly4injXAp4EH034fbwHW0d5Y7kJ7U9wLOBj4OPCx7tBj5znlB4F9k9y0SybnPJA2PPqBkbK3AXsC76L1eNyaNkR09yQPqKpfbyD2tbQ36Td0Re8H3jXedpKbAd8GtqVd3x8CN6clENcBLk/rqfgq8HvAO4Hv0pKCR9F6yH6xUCzzuAFwePec/xP4eVe+J/B/tOG1/6O9nv8R2Br4+5G4d6P9bs6m/V5+1sX3iO7xR4CDaP8Px4e99gC+VFVnLiLOU4Bv0oZaD+/Kdge26WL/m7H6cwn+exdxbmi93DUhxkd2Pw+ZdFBVXZjkk8DTk+xYVSfCb4cnF+Piqrra8OYiHU/7//MA2ute1wRV5eY2+Abcg/ZH8yHd4wCnA28eqfNHXZ1HjB17GHDSyOOnAVcADxyr9+zu+PuPlFVX904TYrr22OMtgOOAL4yU7dyd401jdd/Vle87UvZ24CzgBmN130/7ZH3t8RjG6p0CHEF7o74h8PvdsQUc0NXZoXt8IXCjseM/T0tk1o6UhZZ0/Xik7E3dOe49UnajLsYCdhgp/xLtjXfu8TO6Oi+YEH+6nzccvzYjdfZtf5Z++/j2Xd3njdU7CLh47prREq4CnjpW76GTyue5vo/r6u7YPV5He1N8/li993SvmXsu8Bxf1Z3rMQvU2XP8enblu3blu45d5wKevaHXaVf2b7TEeG33eA1wcvca2nZSPN2//4uWgG82Unb3ru09N3D95p7PPWnzpy4a+f18CPifkdfxp0eO+1h33LYLnX+k/nuBX0woPwb45QaOfUHX1iNHymqR28Tnv9Dreazej4D/XsxzdNs0NocItVrsQftE/kXo3mFb78WTR4Y6/of2yf+3cyi64Z/durpznkDrtfphkhvObd3xAH841vaXq+r48YCq6rc9Nl0729B6JXYeqfaw7ue/jh3+1tEHSUJ7Az+0ezga1xHduXdmw/6I1vN1Lm0y7xNobzgvHqv30aqaG4ojyfVpPRsfAtaNtH2Drv3bjQz//DHwrar69si1OJfWk7chj6P9jt46vqP7nS5JVf2YNrF/9He+htZbdOjI7+gJtKTyc2PX9mhar87473ySPYCjquvZqKqLacn76FD1ZrQe0EOr6qgFnuPjgO9V1ys4T52lWk9L3MfPN/o6Xdc976/SetPu2O26O61H781V9csF4jmE1jM4er32oCWaVxu+W8CHaDcrPCLtRoFHMGF4sLN19/PiRZ77BrSpBOPWLeIcc/u3HinbbZHbEYuMbz4X0JIxXUM4RKjBdW+YT6YlV7duuQgAR9LuvHkw8Nmq+k2SjwJPTbK2Wnf9Y2k9S6MJ1u1oQx/nMtmNxx7/dJ64HgG8DNiJNgQ5Z/QN6VbAlRPOceLY4xvRhpT+stsWE9ckR3YxFfAr4ITxN8zOeDw70nqrXt1t87V/Ju05TVom40eLiO+2wI+q6jeLqLtYHwT+JcnNqw1R7UqLdfx3vg3tDrRJFry2SbalJZYH5qp3r36ddov97btk70a0N+fxu9rG3ZalJSSLcWZNuKGgG/r9J1oCvfXY7m1G4oENx/052hDiHsAXuoTyKcAnu4RzUarq3CSfp01svw6tB+0j81S/qPu5jtZLuhiZUHYxG05g5u4K/e1zqarxuVwbS7jq3w5t4kywtBr8P2A7WpL15An79wA+2/37A7Shvt1pt3A/EfhhVX1vpP5mtKG8v5unvdPHHl9tblGSP6DNv/oK8Fe0N51f04bAJt2CviFzvcX/SRtimmS+eUijfrHIN4Tx5zTX/uuZ/5P4eFK4WnwQ2I/WS/Vm2u/8QuAzI3U2oyVX8915Ol+yPecJtCT6hUy+nX4P2tISfZrvzfZqk9M7k16n2wJfpiUprwBOoi2TsDNt/aUljVJU1RVJ/gt4VpK/os0ZuhntdbtU/wX8B3BT4PB5PghAm8MGbZ7kVxdx3vNodwqOOwHYKcktq+q0eY79/e7nb3usuzmHi3HhaG/hFK4H/GQZx2vGmGBpNdiD9ub43An7Hgs8Jslzuj9uX6ElO09K8jVacvbPY8ecBNyNNldq2k+Mj6O9UT20Ria2JnnGWL1TaW9it+aqfzzH1/A6l/apec0KfmIedXL389eLaP9UWo/QuDssop2TgPsk2aLmn1S+pN9JVf00ybdpv/MDaa+JT9RVJxyfBDwE+PqUb4J70Hp3XjVh37NpSfUrab/Hi2gT9hdy0iLqzA1zbTtWfqsNHDdqV9qQ2WOr6itzhWlLnozHA7+70WAhh9CSzEfSPsicy3TDYx+n3XhwX0aGeCc4FNiHNndyMQnWD4E9kmxTVReOlH+a1tv2Z7QevatIsjXwJ7QPZKMfJs5eRJvQPly9e5F1x9venHYzyqemOV6zyTlYGlSSa9PeMD9dVR8Z32h3262ju9Ooqq6kDTU8krbY5uZcdagI2vyPmwPPmtRekq0WEdoVtETgt70J3e30jx6rN/fG81dj5X89+qDabeIfpQ03Xe2NN8mNFhHT1KrqHNpE6Wcn2W4D7f83cN8k9x7bv5h1yT5KG6Z53oQ25oZ1ftX93HYxsXc+SHujfmZ3/km/8zXAyye0u3nX0zNRklsAuwAfmuc1+C5gxyT36V5/n6AtQnrPCeeae44fBe6WCavsj9SZS3p2Gdm3hvmHkCe5Yu7QkXNsydVfj9+lDRs/f/xajMQDQFUdS+tN/QvaB40PTDPkW1X/R7tjdF9aEjVfvW/SeiP/Ismjx/cn2TIjS7DQ7lAM7caYUR+h9UztPf676YY6D6b1Io0n0SsxB+tOtLt6v7GMc2jG2IOloT2KlkDN98nuW/xu0dG5N9UP0hKYVwHHVdUJY8e8lzaM9G9J/pA2j2YNbcLvE2l3ll1tgvKYw2hDjJ/phkxuTOthO5HfDTNQVUd388Ken+QG/G6Zhrmv9hjtrdmbNnn4yCT/QXszuD5tOOch3b83pufSVss+rmv/ZOAmwP1oSwfcrav3Wlry+pkkb+F3yzScyshzn8chtB6EN3YJ2ldp6/88hHYjwCer6tIkx9N6pH4MnA98v6oWmh/0Idrw5uu7+lfphamqLyd5G7BPkp1oQ8q/pvXEPYG25MZ8c4CeSnvDnu81+N/Ab2ivwSOBl9BuNvhykrllQLbr2nkgbR7R62gT8T+c5J20yfbXp73en0ObAP+DJN8C9utuQjifNkS+lL/L36D1hL0nyQG019ufMjZHqaquTLIXLdH537T15M6m/Z+4M+3/xKhDaNcaphsenGt3vuHwcX9G+519LMmhwBdor7vb0a7JdsDcWlhfow0TPoTf3bhCteVdHt8d+7XuOR5FS+SfSvt/9oaqGl3aY+o5WEn+lNbbOLfg6y5p68UBvLeqTh2pvhvtg8Wq/+5E9Wjo2xjdrtkb7U3tUuA6C9R5F3A53fIGtDeP02hvJi+d55gtgH+gDftcRnvzOoo2T2XrkXoFHDjPOZ4J/Lg7/gTabej7MrKMQFfvOrSetvNow4Af53fLC7x4rO6Nu7qndc/pbFqy8KxFXKtTGLm9fZ46O3Ttvmie/behzQE7u2v/DNqb7uPG6t2V1uN1aVfnZd31WHCZhq7s2rQhmpNHnuOHgduM1Llf9/tYz8gt7pOu78gxX+vq/scCz/9Z3Xl/RRvKO5Y2F2m7BY45Fjh1A9f1i7S7XDfvHt+yu47ndK+Pk7rf65Yjx1yfdjflGd3zPJ02xHSDkTq3ob3pXkZbl+qfaYnDpGUavj9PbPen9er8inaTwmv43ZImu47VfQAtkbmIdnfl9xhbAqOrd1NaUvmjJfxf3rNr82rLVyzmddy9bl5IW2Ps4u6a/Zi2ptltx+q+BfjJPOe/EW0ts5901/WC7ho/crHPZZHP90vMv6TD+HX/Fi3p6q19t9W/za3HIqlHXS/KMcDTqmoxyxtIq0a31MPZwD9W1Xx3nQ4myW1oc7F2r6ovDB3PQrq/Bd8Fdq6q/x02Gq0k52BJy9TNIxv3fNryDV+ZsE9a7fakDasvdnX1FVVVJwPvoA27r3Z7Ax8xubrmsQdLWqYkr6RNuP0ibVhl927796p69pCxSUuR5P/RJmS/GvhiVa3aL7uWVjsTLGmZuu94eyXtjem6tPlV7wX+ufpdcFPaqNK+uPv+tBtDnlaL++5BSROYYEmSJPXMOViSJEk9M8GSJEnq2UwvNNqtQHwzFv8t7JIkScu1DjirFphnNdMJFi25OmPoICRJ0jXO9rTFfSea9QTrYoDTTz+drbfeeuhYJEnSJu6iiy7iFre4BWxg9GzWEywAtt56axMsSZK0ajjJXZIkqWcmWJIkST0zwZIkSeqZCZYkSVLPTLAkSZJ6ZoIlSZLUMxMsSZKknplgSZIk9cwES5IkqWcmWJIkST0zwZIkSeqZCZYkSVLPTLAkSZJ6ZoIlSZLUMxMsSZKknm0+dACSJM2yHfY+rPdznrL/w3s/p1aWPViSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9GzTBSrImyauT/DTJpUlOSvLyJBkyLkmSpOXYfOD2XwzsBTwd+AFwT+BdwIXAAQPGJUmSNLWhE6z7A5+sqsO6x6ckeQpw7wFjkiRJWpah52B9A3hwktsDJLkb8EDg8EGjkiRJWoahe7D2B7YGfpjkCmAN8NKqet+kyknWAmtHitZt/BAlSZKWZugerCcCewBPBXamzcV6UZKnz1N/H9r8rLntjJUIUpIkaSmG7sF6HbB/VX2ge3xcklvREqn3TKi/H/DGkcfrMMmSJEmrzNAJ1nWAK8fKrmCenrWqWg+sn3vsag6SJGk1GjrBOhR4aZLTaMs03B34O+Cdg0YlSZK0DEMnWH8NvBr4V+DGwFnA24B/HDIoSZKk5Rg0waqqi4Hnd5skSdImYei7CCVJkjY5JliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ4NmmAlOSVJTdgOGjIuSZKk5dh84PbvBawZeXwX4HPAh4cJR5IkafkGTbCq6tzRx0n2Bk4CvjxMRJIkScu3auZgJdkSeBrwzqqqoeORJEma1tBDhKMeDWwLvHu+CknWAmtHitZt1IgkSZKmsGp6sIA/Bw6vqrMWqLMPcOHIdsZKBCZJkrQUqyLBSnIr4CHA2zdQdT9gm5Ft+40cmiRJ0pKtliHCZwDnAIctVKmq1gPr5x4n2chhSZIkLd3gPVhJNqMlWO+pqt8MHY8kSdJyDZ5g0YYGbwm8c+hAJEmS+jD4EGFVfRZwrE+SJG0yVkMPliRJ0ibFBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1LPBE6wkN0/yn0nOS3JpkuOS3HPouCRJkqa1+ZCNJ7ke8HXgi8DuwLnA7YALhoxLkiRpOQZNsIAXA6dX1TNGyn46VDCSJEl9GHqI8FHAUUk+nOScJMckedbAMUmSJC3L0AnWbYC9gJ8ADwUOBg5I8vRJlZOsTbL13AasW7lQJUmSFmfoIcLNgKOq6iXd42OS3AV4DvCeCfX3AV65UsFJkiRNY+gerLOB48fKTgBuOU/9/YBtRrbtN15okiRJ0xm6B+vrwB3Gym4PnDqpclWtB9bPPU6y8SKTJEma0tA9WG8C7pvkJUl2TPJU4C+BgwaOS5IkaWqDJlhV9R3gMcBTgO8DLweeX1XvGzIuSZKk5Rh6iJCq+jTw6aHjkCRJ6suSE6wktwb+ALgVcB3a6uvHAN+sqsv6DU+SJGn2LDrBSrIH8LfAPYGfA2cBlwLXB24LXJbkfcBrqmriJHVJkqRrgkUlWEmOAS4H3g08rqpOH9u/Frgf8GTayux/VVUf7jlWSZKkmbDYHqy9q+qI+XZ2yyd8CfhSkpcCOyw/NEmSpNm0qARroeRqQt3zgPOmjkiSJGnGTTPJ/abAfYCbdkU/A46sqp/1GZgkSdKsWsok962At9HmWRVwfrfr+m133g88u6p+1XuUkiRJM2QpC42+Bbg38HDgWlV1k6q6CXAt4I+7fW/pP0RJkqTZspQE63HAnlV1RFVdMVdYVVdU1WeBZwKP7ztASZKkWbOUBGsz2lIN87l8ieeTJEnaJC0lIfo08O9J7j6+oys7GDi0r8AkSZJm1VISrOfRVnA/Osl5SU7otvOAo4BzujqSJEnXaIu+i7CqLgB2T/J7wH256jIN36yqH26E+CRJkmbOktfBqqoTgBM2QiySJEmbhCUlWEm2BB5N+97B0R6sbwCfrKqFJsFLkiRdIyx6DlaSHWk9V+8B7t4du1n370OAH3R1JEmSrtGW0oN1MHAccPequmh0R5KtaUnWQcBD+wtPkiRp9iwlwXoAcO/x5Aqgqi5K8nLgyN4ikyRJmlFLWabhl8AOC+zfoasjSZJ0jbaUHqy3A4ckeTXwBdqaWAA3AR4MvAx4a7/hSZIkzZ6lrIP1iiSXAH8PvAGobldodxK+pqpe23+IkiRJs2VJyzRU1WuA1yS5NSPLNFTVT3uPTJIkaUYteaFRgC6hMqmSJEmaYCmT3BeU5BZJ3tnX+SRJkmZVbwkWcH3g6T2eT5IkaSYteogwyaM2UOU2y4xFkiRpk7CUOVifoN05mAXq1AL7JEmSrhGWMkR4NvDYqtps0gbsvJFilCRJmilLSbCOBu6xwP4N9W5JkiRdIyxliPB1wFYL7D8R+MPlhSNJkjT7lrKS+1c3sP8S4MvLjkiSJGnGLWuZhiRPSbJQr5YkSdI1znLXwXob7cuep5Jk3yQ1tv1wmTFJkiQNaqqvyhnRx6T2HwAPGXn8mx7OKUmSNJjlJlh9+E1V/WzoICRJkvqy3CHC3YEzl3mO2yU5K8nJSd6X5JbzVUyyNsnWcxuwbpltS5Ik9W6qBCvJNknuAJwLXGsZ7R8J7Ak8DNgLuDXw1STzJU77ABeObGcso21JkqSNYkkJVpK/SHI8cD5wfLedn+T4JH++1Mar6vCq+nBVHVtVRwB/DGwLPHGeQ/YDthnZtl9qm5IkSRvbUr7s+e+BfYEDgCOAn3e7bgL8EfCWJNerqtdPG0xV/TLJj4Ed59m/Hlg/EtO0TUmSJG00S5nk/jzgGVX1obHyE4AvJfkebbX3qROsJNcFbgu8d9pzSJIkDW0pQ4Q3Bo5bYP9xwA2X0niS1yd5UJIdktwf+DhwBfD+pZxHkiRpNVlKgvUdYO8kV+v1SrIGeHFXZym2pyVTPwI+BJwH3Leqzl3ieSRJklaNpQ4RHgH8LMlXuOocrF2Ay2lzsRatqp68lPqSJEmzYNE9WFV1LHB74OXAxcBtuu1i4GXAHavq+xsjSEmSpFmypJXcq+pi4OBukyRJ0gSL6sFKstVSTrrU+pIkSZuSxQ4Rnphk7yTbzVchzW5JDgf+pp/wJEmSZs9ihwh3Bf4F2Ldb7+oo4CzgMuB6wJ2A+wG/oa22/rbeI5UkSZoRi0qwqupHwOO6L2J+AvAHwP2BawO/AI4BngUcXlVXbKRYJUmSZsJSJ7mfBryh2yRJkjTBkr7sWZIkSRtmgiVJktQzEyxJkqSemWBJkiT1zARLkiSpZ0tOsJKckuQV3ZINkiRJGjNND9abgccCJyf5XJInJ1nbb1iSJEmza8kJVlW9uap2Au4NnAC8FTg7yYFJdu45PkmSpJkz9RysqvpuVf0NcDPgVcBfAN9J8r9JnpkkfQUpSZI0S5a0kvuoJFsAjwGeAewGfAt4B7A97XsLHwI8tYcYJUmSZsqSE6xuGPAZwFOAK4FDgBdU1Q9H6nwc+E5fQUqSJM2SaXqwvgN8DtgL+ERV/XpCnZ8CH1hOYJIkSbNqmgTrNlV16kIVquoSWi+XJEnSNc40k9xvnOQ+44VJ7pPknj3EJEmSNNOmSbAOAm4xofzm3T5JkqRrtGkSrDsB351Qfky3T5Ik6RptmgRrPXCTCeXbAb9ZXjiSJEmzb5oE67PAfkm2mStIsi1t7avP9RSXJEnSzJrmLsIXAV8BTk1yTFe2E/Bz4E97ikuSJGlmLTnBqqozk/w+sAdwN+BS4F3A++dZE0uSJOkaZaqvyunWufr3nmORJEnaJCznuwjvBNwS2HK0vKo+tdygJEmSZtk030V4G+DjwF2BAtLtqu7nmn5CkyRJmk3T3EX4Ftp3Dd4Y+BVwZ2AX4Chg194ikyRJmlHTJFj3A15RVb8ArgSurKqvAfsAB0wbSJK9k1SSN097DkmSpNVgmgRrDXBx9+9fADfr/n0qcIdpgkhyL+DZwLHTHC9JkrSaTJNgfZ+2PAPAkcA/JHkA8Arg5KWeLMl1gfcBzwIumCIeSZKkVWWaBOufRo57BXBr4KvAHwN/M8X5DgIOq6rPb6hikrVJtp7bgHVTtCdJkrRRTbPQ6BEj/z4RuGOS6wMXVFXNf+TVJXkysDNwr0Uesg/wyqW0IUmStNKW1IOVZIskv0lyl9Hyqjp/iuTqFrQ7EveoqssWedh+wDYj2/ZLaVOSJGklLKkHq6p+neQ0+lnr6h60pR6+m8wtpcUaYJckzwPWVtUVY+2vB9bPPR45TpIkadWYZg7WPwP/0g0LLscXaIuV7jSyHUWb8L7TeHIlSZI0K6b5qpznATsCZyU5FbhkdGdV7byYk1TVxbQ7En8rySXAeVX1/clHSZIkrX7TJFif6DsISZKkTck0dxG+amME0p171411bkmSpJUyzRwsSZIkLWDJPVhJrgTmXZKhqvq4w1CSJGlmTTMH6zFjj7cA7g48HRcBlSRJmmoO1icnFH8kyQ+AJwHvWHZUkiRJM6zPOVjfAh7c4/kkSZJmUi8JVpJr077o+cw+zidJkjTLppnkfgFXneQeYB3wK+BpPcUlSZI0s6aZ5P4CrppgXQmcCxxZVRf0EpUkSdIMm2aS+7s3QhySJEmbjCXPwUryjCRPmFD+hCRP7ycsSZKk2TXNJPd9gF9MKD8HeMnywpEkSZp90yRYtwR+OqH81G6fJEnSNdo0CdY5wO9PKL8bcN7ywpEkSZp909xF+H7ggCQXA1/pyh4EvAX4QF+BSZIkzappEqyXAzsAXwB+05VtBhyCc7AkSZKmWqbhcuBJSV4G7ARcChxXVaf2HJskSdJMmqYHC4Cq+gnwkx5jkSRJ2iRMsw7WR5O8eEL5PyT5cD9hSZIkza5perB2AfadUH448MJlRSOtcjvsfVjv5zxl/4f3fk5J0rCmWabhusDlE8p/DWy9vHAkSZJm3zQJ1nHAkyaUPxk4fnnhSJIkzb5phghfDXwsyW2B/+nKHgw8BbjadxRKkiRd00yzTMOhSR5NW/Pq8bRlGo4FHlJVX+43PEmSpNkz1TINVXUYcLXZvknuUlXfX3ZU2qicqC1J0sY1zRysq0iyLslfJvk28L0eYpIkSZppUydYSXZJcghwNvAi2nys+/YVmCRJ0qxa0hBhkpsCewJ/TluS4UPAWuDRVeUdhJIkSSyhByvJocCPgN8Hng/crKr+eiPFJUmSNLOW0oO1O3AAcHD3PYSSJEmaYClzsB4IrAOOTnJkkuclueFGikuSJGlmLTrBqqpvVdWzgO2At9FWbj+rO8duSdZtnBAlSZJmy5LvIqyqS6rqnVX1QOCuwBuAvYFzknxqKedKsleSY5Nc1G3fTLL7UmOSJElaTZa1DlZV/aiq/gHYnvZVOUt1Bi05uwdwT9pSD59McuflxCVJkjSkqVZyH1dVVwCf6LalHHfoWNFLk+xFW0/rB33EJkmStNJ6SbD6kGQN7cuitwK+OU+dtbR1t+Y470uSJK06y/6qnOVKctck/wesB/4NeMwCi5buA1w4sp2xMlFKkiQt3uAJFm3x0p2A+wAHA+9Jcqd56u4HbDOybb8SAUqSJC3F4EOEVXU5cGL38Ogk9wL+Fnj2hLrraT1dACRZkRglSZKWYjX0YI3bjKvOs5IkSZopg/ZgJdkPOBw4jTZh/anArsBDBwxLkiRpWYYeIrwxcAhtdfgLgWOBh1bV5waNSpIkaRkGTbCq6s+HbF+SJGljWI1zsCRJkmaaCZYkSVLPTLAkSZJ6ZoIlSZLUMxMsSZKknplgSZIk9WzodbA0Yoe9D+v9nKfs//Dezylp5fh3QZpN9mBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzEyxJkqSemWBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzEyxJkqSemWBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzEyxJkqSeDZpgJdknyXeSXJzknCSfSHKHIWOSJElarqF7sB4EHATcF9gN2AL4bJKtBo1KkiRpGTYfsvGqetjo4yR7AucA9wC+MkRMkiRJyzVogjXBNt3P8yftTLIWWDtStG6jRyRJkrREQw8R/laSzYA3A1+vqu/PU20f4MKR7YyViU6SJGnxVk2CRZuLdRfgyQvU2Y/WyzW3bb8CcUmSJC3JqhgiTHIg8Ahgl6qat1eqqtYD60eOW4HoJEl92mHvw3o/5yn7P7z3c0rLMWiClZYhvRV4DLBrVf10yHgkSZL6MHQP1kHAU4E/AS5OctOu/MKqunS4sCRJkqY3dIK1V/fzS2PlzwDevaKRSNokOPwkaTUYeh0sJ1FJkqRNztA9WJIkaZWxJ3j5VtMyDZIkSZsEEyxJkqSemWBJkiT1zDlYkiRpEJvyXC97sCRJknpmgiVJktQzEyxJkqSeOQdLkrRJ2pTn92j1swdLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzl2mQpCm4BICkhdiDJUmS1DMTLEmSpJ45RLgIDgVIkqSlsAdLkiSpZyZYkiRJPTPBkiRJ6plzsCRJgPNNpT7ZgyVJktQzEyxJkqSemWBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pnLNEirkLfLS9JsswdLkiSpZ4P2YCXZBfh74B7AdsBjquoTQ8akfvXdE2MvjCRpFgzdg7UV8D3guQPHIUmS1JtBe7Cq6nDgcIAkQ4aiGeecJUnSajJ0D5YkSdImZ6buIkyyFlg7UrRuqFgkSZLmM2s9WPsAF45sZwwbjiRJ0tXNWoK1H7DNyLb9sOFIkiRd3UwNEVbVemD93GMnxkuSpNVo6HWwrgvsOFJ06yQ7AedX1WnDRCVJkrQ8Q/dg3RP44sjjN3Y/3wPsueLRSJIk9WDodbC+BDjOJ0mSNimzNsldkiRp1TPBkiRJ6tnQc7AkSdIi+bVgs2PTSLAuuQTWrNlop7/25Zf1f9JLLtnk29koba1UO/O0tam1syna1H5HM/3/dZ62bGe6dlayrU2tnSHOn6rauIFsREm2Bi68ENh66GAkSdIm7yLaSufANlV10Xz1nIMlSZLUs01jiPCss2DrjdeH9Xsv/0zv5zzh1Q/b5NvZGG2tVDvztWU707Wzkm3ZznTtrGRbtjNdOyvZ1qbWTq8uughudrMNVts0EqyttmrbRnLpltfq/6QT4t3U2tkoba1UO/O0ZTvTtbOSbdnOdO2sZFu2M107K9nWptZOr664YlHVHCKUJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnJliSJEk9M8GSJEnqmQmWJElSz0ywJEmSemaCJUmS1DMTLEmSpJ6ZYEmSJPXMBEuSJKlnqyLBSvLcJKckuSzJkUnuPXRMkiRJ0xo8wUryJOCNwKuAnYHvAUckufGggUmSJE1p8AQL+DvgP6rqXVV1PPAc4FfAM4cNS5IkaTqDJlhJtgTuAXx+rqyqruwe32+ouCRJkpZj84HbvyGwBvj5WPnPgTuOV06yFlg7UrQO4KKLLtpY8QFw5fpf9X7OSTFvau1sjLZWqp352rKd6dpZybZsZ7p2VrIt25munZVsa1NrZ4jzp6o2aiALNp7cDDgTuH9VfXOk/LXAg6rqPmP19wVeuaJBSpIkXd32VXXmfDuH7sH6BXAFcJOx8psAP5tQfz/ahPhR1wfO7z+0JVsHnAFsD1w8cCyzxms3Pa/ddLxu0/PaTc9rN73Vdu3WAWctVGHQBKuqLk9yNPBg4BMASTbrHh84of56YP1Y8cbtC1ykJHP/vLiqVkVMs8JrNz2v3XS8btPz2k3Paze9VXjtNhjD0D1Y0Hqk3pPkKODbwPOBrYB3DRmUJEnStAZPsKrqg0luBPwjcFPgf4GHVdX4xHdJkqSZMHiCBVBVBzJhSHDGrKctljo+hKkN89pNz2s3Ha/b9Lx20/PaTW/mrt2gdxFKkiRtilbDSu6SJEmbFBMsSZKknplgSZIk9cwES5IkqWcmWD1J8twkpyS5LMmRSe49dEyrXZK9khyb5KJu+2aS3YeOaxYkuXmS/0xyXpJLkxyX5J5DxzULkqxL8uYkp3bX7htJ7jV0XKtNkl2SHJrkrCSV5NEj+7ZI8prudXdJV+eQ7uvPrvEWunbd/nd35aPbZwYKd9VYxHW7bpIDk5zR/d89PslzBgp3g0ywepDkSbQFU18F7Ax8DzgiyY0HDWz1OwPYG7gHcE/gf4BPJrnzoFGtckmuB3wd+DWwO3An4IXABUPGNUPeDuwG/ClwV+CzwOeT3HzQqFafrWh/y547Yd91aH/rXt39fCxwB+BTKxbd6rbQtZvzGWC7ke0pKxDXareh6/ZG4GHA04DfA94MHJjkUSsS3RK5TEMPkhwJfKeqntc93gw4HXhrVe0/aHAzJsn5wN9X1TuGjmW1SrI/8ICq+oOhY5k1Sa5N+x6zP6mqw0bKjwYOr6qXDRbcKpakgMdU1ScWqHMv2rdx3KqqTlup2Fa7SdcuybuBbavq0QOFterNc92+D3ywql49UrZq/+/ag7VMSbak9cB8fq6sqq7sHt9vqLhmTZI1SZ5M+wTzzaHjWeUeBRyV5MNJzklyTJJnDR3UjNgcWANcNlZ+KfDAlQ9nk7INUMAvB45jVuza/f/9UZKDk9xg6IBmwDeAR3VTJJLkD4Hb03qhVx0TrOW7Ie0P9vhX+/yc9tU/WkCSuyb5P9rqvP9G+8Ry/MBhrXa3AfYCfgI8FDgYOCDJ0weNagZU1cW0BP7lSW7WJfZPo30Y2m7Y6GZXkmsBrwHev0q+iHe1+wzwZ8CDgRcDDwIOT7Jm0KhWv78GjqdNL7mcdh2fW1VfGTSqeayKr8rRNdqPgJ1on34fT/vi7weZZC1oM+CoqnpJ9/iYJHcBngO8Z7iwZsafAu8EzgSuAL4LvJ/WE60lSrIF8CEgtMRfG1BVHxh5eFySY4GTgF2BLwwS1Gz4a+C+tF78U4FdgIOSnFVVn1/wyAHYg7V8v6D9kb7JWPlNgJ+tfDizpaour6oTq+roqtqHNsHxb4eOa5U7m/YpbtQJwC0HiGXmVNVJVfUg4LrALarq3sAWwMnDRjZ7RpKrWwG72Xs1nao6mfZesuPQsaxW3fzJfwH+rqoOrapju+8x/iDwomGjm8wEa5mq6nLgaFpXL/DbSe4PxrlE09gMWDt0EKvc12l3bI26Pe0TnRapqi6pqrO7uzIfCnxy6JhmyUhydTvgIVV13sAhzawk2wM3oH140mRbdNuVY+VXsEpzGYcI+/FG2tDWUbS7aJ5Pm6z9riGDWu2S7AccDpwGrAOeSusif+iAYc2CNwHfSPIS2hvcvYG/7DZtQJKH0oazfkTrMXgd8EP8/3oVSa7LVXtUbp1kJ+B8WiLwEdoSDY8A1iSZm3N6fvfB8xprA9fufOCVwEdpoxy3BV4LnAgcsbKRri4LXbeqOi3Jl4HXJbmU9oHyQbS5bH+34sEuRlW59bABz6P9wtcDRwL3GTqm1b4B7wBO6a7ZObQ7L3cbOq5Z2GhvasfR7oY7AXjW0DHNygY8kTbfZT0tUTgQ2GbouFbbRvuwUxO2dwM7zLOvgF2Hjn3obQPX7tq0ROoc2kTtU4B/B24ydNxDbwtdt27/TWkfhM6k3fn7Q1pylaFjn7S5DpYkSVLPVuW4pSRJ0iwzwZIkSeqZCZYkSVLPTLAkSZJ6ZoIlSZLUMxMsSZKknplgSZIk9cwES9KKSlJJHj10HCspyZ5Jfjl0HJJWjgmWpN4kuWmStyY5Ocn6JKcnOTTJgzd89FTt7dolbNtujPN3bVzjEkJJy+d3EUrqRZIdaF9E/Uvg72lf5bMF7bslDwLuOFRsG5IkwJqq+s3QsUjaNNiDJakv/0r73rB7V9VHq+rHVfWDqnojcN9JB0zqgUqyU1e2Q/f4Vl0v2AVJLknygyR/3O3/YnfYBd0x7+6O2SzJPkl+muTSJN9L8vgJ7e6e5Gja9xI+cENPMMkO3XGPTfLFJL/qzn2/sXp7Jjmt2/9x4AYTzvUnSb6b5LKux++VSTbv9r0iyVlJbjBS/7CuTf9uSzPA/6iSli3J9YGHAQdV1SXj+6vql8s4/UHAWmAX4K7Ai4H/A04HHtfVuQOwHfC33eN9gD8DngPcGXgT8J9JHjR27v2BvYHfA45dQkz/DLwe2An4MfD+keToPrQvMj+w2/9F4GWjByf5A+AQ4C3AnYBnA3sCLx05/ynA27v6zwXuDzy9qq5cQpySBuIQoaQ+7AiE9u32fbsl8NGqOq57fPLcjiTnd/88Zy6JS7IWeAnwkKr65twxSR5IS2S+PHLuV1TV56aI6fVVdVjX3iuBH9CuwQ9pSd5nquq1Xd0fJ7k/LQGd80pg/6p6z0h8LwdeC7yqqq5I8jTgf5PsD/wN8BdVddoUsUoagAmWpD5kI577AODgJH8EfJ6WbC3U27QjcB3gc21q1W9tCRwzVveoKWMabf/s7ueNaQnW7wEfH6v/Ta6aYN0NeECSl46UrQGuleQ6VfWrqjo5yYuAtwEfrKr/mjJWSQMwwZLUh5/Q5l8tdSL73HDXaCa0xWiFqnp7kiOAhwN/BOyT5IVV9dZ5znnd7ufDgTPH9q0fe3y14cxF+vVoiN3PpUy5uC6tF+tjE/ZdNvLvXYArgB2SbO4kfGl2OAdL0rJV1fnAEcBzk2w1vn+BZRTO7X5uN1K204Tzn15V/1ZVjwXeADyr23V593PNSPXjaYnULavqxLHt9MU+p2U4AbjPWNn4JP/vAneYEN+Jc3OskjwJeCywK22Y9OUbOW5JPbIHS1JfnktbpuHbSV5BG0bbHNgN2Is2dDbuRNpk9X274bLbAy8crZDkzcDhtMnk1wP+kJbEAJxK60F6RJL/Bi6tqouTvB54U3fH3deAbYAHABeNzHvaWA4Avt4N732StkzFw8bq/CPw6SSnAR+h9eTdDbhLVb0syfbAwcCLq+prSZ7R1T+8qr61keOX1AN7sCT1oqpOBnam3TX3BuD7wOeAB9MSrEnH/Bp4Cm1o8VjaHYIvG6u2hnYn4QnAZ2iJ1l91x59JN2Ec+Dntzj1ovT2vpt1NOHfcw4GfLvuJbkCXAD2LNtn9e7RhzX8aq3ME8Ihu33eAbwEvAE7t1uR6N/BtuufT1T+YdifkdZG06qWqNlxLkiRJi2YPliRJUs9MsCRJknpmgiVJktQzEyxJkqSemWBJkiT1zARLkiSpZyZYkiRJPTPBkiRJ6pkJliRJUs9MsCRJknpmgiVJktQzEyxJkqSe/X/w3Fk4v5GQwwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.ticker import MaxNLocator\n",
    "ax = plt.figure(figsize=(7,5), dpi=100).gca()\n",
    "ax.bar(np.arange(NumClusters), np.array(AvgAccPerCluster))\n",
    "ax.axhline(y=BaselineAcc, color='r', linestyle='-')\n",
    "ax.set_title('Average Predictive Accuracy MC(O=1)')\n",
    "ax.set_xlabel('Cluster Index')\n",
    "ax.set_ylabel('Accuracy (0-100)')\n",
    "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusterLens = [len(c) for c in clusters]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAHPCAYAAADj6/a/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvQUlEQVR4nO3deZhkZX238fvLgKgwgBsDiIgRRUUMiiKgARQQFDSKJjGYCMZoRN64Jgoq4pYAYlAjhCguYOKSBeM27Bo04ohCREBERBz2RR1hWAeE3/vHOW1qiu6e6erqrtPT9+e6zlVd5znLr+qa5dvPOc9zUlVIkiSpe9YadQGSJEkan0FNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFN0pySZGmSE0ddxygkOTHJ0lHXIWn2GNQkdUKSxyb5eJIrktyVZHmSc5K8McmDZqmGByd5T5LdZuN8Ped9RJKPJrk0yZ1Jbkry/SRHJVl/NmuR1C1rj7oASUqyD/AfwArgs8DFwAOAZwNHA9sAr52FUh4MHN7+fPYsnI8kDwXOAzYAPg1cCjwMeApwEHA8cFu7+WvwF2xpXjGoSRqpJI8BvghcCTy3qq7vaT4uyVbAPiMpbkiSrFdVt0/Q/GpgC+BZVfXdvv02AO4ee19V98xclZK6yN/MJI3a24D1gVf3hTQAquryqvroRDu3lyprnPUHJqkkW/ase3qS05P8qr3E+Iskn27btgR+2W56eLtvJXlPz/5PSPKfSZa1l2fPS/KiCc67a5J/SnITcM0kn/+xwL3A98b57Mur6q6eY690j1qSs3vq7F8O7NluoyQfSXJ1khVJLk/y9iT+HyB1nD1qkkbthcAV/b1Jw5ZkY+AMmjB2JHAzsCWwX7vJL/m/S43/BXypXX9hu/82wDnAte3+twN/DHw5yUur6r/6TvlP7THfB6w3SWlXAguAPwdOmuLH+jvgk33r/gzYC7iprfvBwLeARwIfB64CdgaOADYF3jTFc0qaRQY1SSPTXtp7JPCVWTjdzsBDgOdV1Xk9698FUFW3J/lPmqB2YVX9a9/+H6UJOc+oqhUASf4J+A5wFE2467UM2L2q7l1FXZ8G3gycmOQQmnvjvg2cUlW3TLZjVZ3Z+z7JzsBzgU9X1Snt6rfQ9No9tap+1q77eJLrgL9N8g9VdfUqapQ0InZ7SxqlDdrXW2fhXDe3r/smWWcqO7Y3/D8X+HdgYZKHJ3k4zU3/pwOPS/LIvt1OWI2QRlXdCPw+8M80QfJ1wOeBm5IcliSrWeMmwH8CFwCv72n6I+B/gN+M1d3WfhZNT94uq3N8SaNhUJM0Ssvb14WzcK5vASfTjOr8VZKvJHlVknVXY9+tgADvp7mc2bu8t91m4759frG6hVXV9VV1EM2lyK2BN/B/l01fvar9k6xNEyIXAPuN9fi1HgfsPU7dZ01Qt6QO8dKnpJGpquXtJbgnT+cwE6xf0HeuAl6WZEea++L2orns+NYkO1bVbeMcY8zYL7UfoulBG8/lfe/vnLTqcbQ1XgZclmQx8DPgFdz/PrR+RwM7AXtUVf/AhbWAM4EPTrDvZVOtU9LsMahJGrWvA69NslNVLRlg/99AM7Kxqm7uWf/o8Tauqu/RjLB8Z5L9gc8BL6cJQxOFviva13uq6qwJthmqqroiyW9oetkmlOTlNAMC3lRV3xpnk58D689W3ZKGy0ufkkbtgzQjKD+ZZFF/Y/vEgjdOsv/P29ff3WuVZD3ggL7jPGSc+70uaF/HLn/e0b5u1LtRVd1Ec5P/XyW5X3BK8ohJ6ptUkme29fav34HmHrifTrLvk2kC5r9OMoXJvwM7JdlrnP03ai+bSuoo/4JKGqmq+nnbs/VvwE+S9D6ZYGeam+FPnOQQZ9CMxvxUkqNp5iT7C5r7sLbo2e4A4PVJ/osm3C2kmel/OXBKW8udSS4B/iTJZTQjNy+uqouBg2lGeF6U5ASaXrZFNJccN6cZEDCIPwde0dZ1Ps0Et09sP8NdwN9Psu9n2tdvJ/mzvrbvVtUVNJdFXwR8Pc0zUs+nmS5kW+BlNFOU/GrA2iXNMIOapJGrqq8meQrwt8Af0sxntoJmDrO3AidMsu89SV5CM2/Z+4EbgI/QXBL9TM+m3wJ2oLnMuQi4Bfg+8Iqq6r3x/y+BjwEfpgmL76UJa5ckeTrNYIQDaXq7bgJ+SHPT/6A+TtOTtzvNZ9+AJmSeARxRVT+cZN9H0ISuT4zT9iqa+enuSLIr8A6a0PtKmnB6WftZJp0CRNJopbl3VZIkSV3jPWqSJEkdZVCTJEnqKIOaJElSRxnUJEmSOsqgJkmS1FEGNUmSpI5yHjWgna18M+DWUdciSZLmhYXAdbWKedIMao3NgP4HGUuSJM2kzYFrJ9vAoNa4FeDqq69mgw02GHUtkiRpDbZ8+XIe9ahHwWpcyTOo9dhggw0MapIkqTMcTCBJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUketPeoCpPlmy0MWD/2YS4/cZ+jHlCSNnj1qkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR11EiDWpKDklyYZHm7LEny/J72s5NU3/LPfcfYIsniJHckuSnJ0UmcH06SJM15ow401wCHAD8DAhwAfCXJU6vqx+02JwDv7tnnjrEfkiwAFgM3ADsDmwKfBe4B3jHj1UuSJM2gkQa1qvpa36p3JjkI2BEYC2p3VNUNExziecCTgD2q6kbggiSHAUcleU9V3T0jhUuSJM2CztyjlmRBkpcD6wFLeppekeRXSS5OckSSB/e07QRc1Ia0MacDGwDbTHKudZNsMLYAC4f4USRJkoZi1Jc+SbItTTB7IHAb8JKquqRt/jxwJXAd8BTgKGBrYL+2fRPgRlZ2Y0/bRA4FDp928ZIkSTNo5EEN+CmwHbAh8DLgpCS7VtUlVfWJnu0uSnI98I0kj62qn0/jnEcAx/S8X0hzv5wkSVJnjDyotfeRXd6+PT/JM4A3An81zubntq9bAT+nGUSwQ982i9rXie5ro6pWACvG3ieZeuGSJEkzrDP3qPVYC1h3grbt2tfr29clwLZJNu7ZZk9gOXAJkiRJc9hIe9SSHAGcClxFc/lxf2A3YK8kj23fnwL8muYetQ8D366qC9tDnEETyP4lydto7kv7AHBc22smSZI0Z4360ufGNPOebQrcAlwI7FVVZyZ5FLAH8CaakaBXAyfTBDEAqureJPsCx9P0rt0OnMTK865JkiTNSaOeR+3Vk7RdDey6Gse4EnjBMOuSJEnqgi7eoyZJkiQMapIkSZ1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUQY1SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUQY1SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUQY1SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUQY1SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUSMNakkOSnJhkuXtsiTJ83vaH5jkuCS/TnJbkpOTLOo7xhZJFie5I8lNSY5OsvbsfxpJkqThGnWP2jXAIcD2wNOBbwJfSbJN2/5h4IXAHwG7ApsBXxrbOckCYDHwAGBn4ADgQOB9s1O+JEnSzBlpz1NVfa1v1TuTHATsmOQa4NXA/lX1TYAkrwJ+kmTHqvoe8DzgScAeVXUjcEGSw4Cjkrynqu6evU8jSZI0XKPuUfudJAuSvBxYD1hC08u2DnDW2DZVdSlwFbBTu2on4KI2pI05HdgA2AZJkqQ5bOT3ciXZliaYPRC4DXhJVV2SZDvg7qq6uW+XG4FN2p83ad/3t9OzzXjnXBdYt2fVwoGKlyRJmkFd6FH7KbAd8EzgeOCkJE+a4XMeCtzSs1wzw+eTJEmaspEHtaq6u6our6rzq+pQ4EfAG4EbgAck2ahvl0VtG+3ronHa6dlmPEcAG/Ysmw/+CSRJkmbGyIPaONaiuSx5PnAPsPtYQ5KtgS1oLpXSvm6bZOOe/fcElgOXTHSCqlpRVcvHFuDW4X4ESZKk6RvpPWpJjgBOpRkgsBDYH9gN2KuqbknyKeCYJMtowtfHgCXtiE+AM2gC2b8keRvNfWkfAI6rqhWz+mEkSZKGbNSDCTYGPgtsSnOv2IU0Ie3Mtv3NwH3AyTS9bKcDrx/buaruTbIvzb1tS4DbgZOAd8/WB5AkSZopo55H7dWraL8LOLhdJtrmSuAFQy5NkiRp5Lp4j5okSZIwqEmSJHXWqO9RkyRpyrY8ZPHQj7n0yH2GfkxpuuxRkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRa4+6AEmSumzLQxYP/ZhLj9xn6MfUmskeNUmSpI4yqEmSJHWUlz41MC8HSJI0s+xRkyRJ6iiDmiRJUkcZ1CRJkjpqSkEtyYIkuyTZaIbqkSRJUmtKQa2q7gXOAB4yM+VIkiRpzCCXPi8Gfm/YhUiSJGllgwS1dwEfSrJvkk2TbNC7DLtASZKk+WqQedROaV+/ClTP+rTvF0y3KEmSJA0W1J4z9CokSZJ0P1MOalX1rZkoRJIkSSsbaB61JH+Q5F+TfDfJI9t1f57k2cMtT5Ikaf6aclBL8lLgdOBO4GnAum3ThsA7hleaJEnS/DboqM/XVdVrgHt61p9DE9xWW5JDk/wgya1Jbkry5SRb921zdpLqW/65b5stkixOckd7nKOT+MB5SZI0pw0SZrYGvj3O+luAjaZ4rF2B44AftLX8PXBGkidV1e09250AvLvn/R1jPyRZACwGbgB2BjYFPksTIu3hkyRJc9YgQe0GYCtgad/6ZwNXTOVAVbV37/skBwI3Aduzchi8o6pumOAwzwOeBOxRVTcCFyQ5DDgqyXuq6u6p1CRJktQVg1z6PAH4aJJn0sybtlmSVwAfAo6fZj0btq/L+ta/Ismvklyc5IgkD+5p2wm4qA1pY04HNgC2Ge8kSdbtm6R34TTrliRJGrpBetSOpAl43wAeTNPztQL4UFV9bNBCkqwFfAQ4p6ou7mn6PHAlcB3wFOAomsuv+7XtmwC9IY2e95tMcLpDgcMHrVWSJGk2DDKPWgF/l+Romkug6wOXVNVt06zlOODJNJdQe8/3iZ63FyW5HvhGksdW1c8HPNcRwDE97xcC1wx4LEmSpBkxyPQcn06ysKrurqpLqur7VXVbkvWSfHqQIpIcC+wLPKeqVhWYzm1ft2pfbwAW9W2zqKftfqpqRVUtH1uAWwcoW5IkaUYNco/aAcCDxln/IOCVUzlQGscCLwGeW1W/WI3dtmtfr29flwDbJtm4Z5s9geXAJVOpR5IkqUtW+9Jne9N92mVhkrt6mhcAL6AZsTkVxwH7A38I3Jpk7J6yW6rqziSPbdtPAX5Nc4/ah4FvV9WF7bZn0ASyf0nyNpr70j4AHFdVK6ZYjyRJUmdM5R61m2lGeRZw2TjtxdRv0D+ofT27b/2rgBOBu4E9gDcB6wFXAyfTBLHmpFX3JtmXZsTpEuB24CRWnndNkiRpzplKUHsOTW/aN4GXsvIUGncDV1bVdVM5eVVlFe1X00yKu6rjXEnToydJkrTGWO2gVlXfAkjyGOCqdvSnJEmSZsgggwmeCDxr7E2Sg5NckOTzSR4yvNIkSZLmt0GC2tE0s/6TZFua+chOAR7DynOTSZIkaRoGeTLBY/i/aS9eCnytqt6R5Gk0gU2SJElDMEiP2t00j46CZkTmGe3Py2h72iRJkjR9g/SofQc4Jsk5wA7An7TrH4+PYZIkSRqaQXrU/h/wW+BlwEFVdW27/vnAacMqTJIkab4b5KHsV9E8l7N//ZuHUpEkSZKAAYJaki0ma2+DnCRJkqZpkHvUltI8LmoiCwYrRZIkSb0GCWpP7Xu/TrvuLcA7p12RJEmSgMHuUfvROKvPS3Id8LfAl6ZdlSRJkgYa9TmRnwLPGOLxJEmS5rVBBhP0T2obYFPgPcDPhlCTJEmSGOwetZu5/2CCAFcDL59uQZIkSWoMEtSe0/f+PuCXwOVV9dvplyRJkiQYbDDBt2aiEEmSJK1stYJakhet7gGr6quDlyNJkqQxq9uj9uXV3K5wwltJkqShWK2gVlXDnMZDkiRJq8EAJkmS1FGrHdSSPDfJJePMo0aSDZP8OMkuwy1PkiRp/ppKj9qbgBOqanl/Q1XdAnwcePOQ6pIkSZr3phLUfh84bZL2M4Dtp1eOJEmSxkwlqC0C7pmk/bfAI6ZXjiRJksZMJahdCzx5kvanANdPrxxJkiSNmUpQOwV4f5IH9jckeRDwXuDrwypMkiRpvpvKI6Q+AOwHXJbkWOCn7fonAAfTTHT7d8MtT5Ikaf5a7aBWVTcm2Rk4HjgCyFgTcDpwcFXdOPwSJUmS5qcpPZS9qq4EXpDkIcBWNGHtZ1X1m5koTpIkaT6bUlAb0wazHwy5FkmSJPXwEVKSJEkdZVCTJEnqKIOaJElSRxnUJEmSOmqgwQRJHgc8B9iYvrBXVe8bQl2SJEnz3pSDWpLX0Myl9ivgBpp51MYUYFCTJEkagkF61N4FvLOqjhp2MZIkSfo/g9yj9hDgP4ZdiCRJklY2SFD7D+B5wzh5kkOT/CDJrUluSvLlJFv3bfPAJMcl+XWS25KcnGRR3zZbJFmc5I72OEcnGej+O0mSpK4YJMxcDrw/yY7ARcA9vY1V9Y9TONauwHE0TzlYG/h74IwkT6qq29ttPgzsA/wRcAtwLPAl4FkASRYAi2nul9sZ2BT4bFvXOwb4fJIkSZ0wSFB7LXAbTcjata+tgNUOalW1d+/7JAcCNwHbA99OsiHwamD/qvpmu82rgJ8k2bGqvkfTu/ckYI/2ofAXJDkMOCrJe6rq7gE+oyRJ0shNOahV1WNmopDWhu3rsvZ1e2Ad4Kye81+a5CpgJ+B77etFbUgbczrNyNRtgB/2nyTJusC6PasWDusDSJIkDUtnJrxNshbwEeCcqrq4Xb0JcHdV3dy3+Y1t29g2N47TTs82/Q6luYw6tlwzcOGSJEkzZLV61JIcAxxWVbe3P0+oqt4yYC3HAU8Gnj3g/lNxBND7ORZiWJMkSR2zupc+n0pzCXLs54nUJG0TSnIssC+wS1X1BqYbgAck2aivV21R2za2zQ59h1zU03b/IqtWACt6zj9I2ZIkSTNqtYJaVT1nvJ+nK01C+hjwEmC3qvpF3ybn04ze3B04ud1na2ALYEm7zRLgnUk2rqqb2nV7AsuBS4ZVqyRJ0mwb9VxjxwH7A38I3Jpk7J6yW6rqzqq6JcmngGOSLKMJXx8DlrQjPgHOoAlk/5LkbTT3pX0AOK7tOZMkSZqTRh3UDmpfz+5b/yrgxPbnNwP30fSorUszovP1YxtW1b1J9qUZ5bkEuB04CXj3TBUtSZI0G0Ya1KpqlTeHVdVdwMHtMtE2VwIvGGJpkiRJI9eZ6TkkSZK0MoOaJElSR005qCU5IMk+Pe8/mOTmJN9N8ujhlidJkjR/DdKj9g7gToAkO9HcO/Y24Fc0D1CXJEnSEAwymOBRwOXtzy8GTq6qTyQ5h/uP3pQkSdKABulRuw14WPvz84Az25/vAh40jKIkSZI0WI/amcAnk/wQeDxwSrt+G2DpkOqSJEma9wbpUTsY+C7wCOClVfXrdv32wBeGVZgkSdJ8N6UetSRrA28Ajup7eDpVdfgwC5MkSZrvptSjVlW/pRnhOepHT0mSJK3xBrn0+Q1g12EXIkmSpJUN0jN2KnBkkm2B82kegv47VfXVYRQmSZI03w0S1P6pfX3LOG0FLBi8HEmSJI2ZclCrKp8PKkmSNAumNSggyQOr6q5hFSNJk9nykMVDP+bSI/dZ9UaSNCKDPJR9QZLDklwL3Jbk99r170/y6qFXKEmSNE8NchnzncCBNNN03N2z/mLgL4dQkyRJkhgsqL0SeG1VfQ64t2f9j4AnDKUqSZIkDRTUHglcPsGx1pleOZIkSRozSFC7BPiDcda/DPjh9MqRJEnSmEFGfb4POCnJI2mC3n5Jtqa5JLrvMIuTJEmaz6bco1ZVXwFeCOxB81SC9wFPBF5YVWcOtzxJkqT5a6B51Krqf4A9h1yLJEmSegwyj9qjkmze836HJB9J8trhliZJkjS/DTKY4PPAcwCSbAKcBewA/F2Sdw+xNkmSpHltkKD2ZOD77c9/DFxUVTsDr6CZCFeSJElDMEhQWwdY0f68B/DV9udLgU2HUZQkSZIGC2o/Bl6X5A9oBhSc1q7fDPj1sAqTJEma7wYJam8H/go4G/hCVf2oXf8i/u+SqCRJkqZpytNzVNXZSR4ObFBVv+lp+gRwx9AqkyRJmucGnUftXuA3feuWDqMgSZIkNVY7qCX5DVDjNN0CXAZ8yCcTSJIkDc9UetTeNMH6jYDtga8neVlVfW26RUmSJGkKQa2qTpqsPckFwKGAQU2SJGkIBhn1OZGvA08Y4vEkSZLmtWEGtXWBu4d4PEmSpHltmEHt1cAFQzyeJEnSvDaVUZ/HTNC0IfA04PHALsMoSpIkSVPrUXvqBMvDgTOBJ1fV+VM5eZJdknwtyXVJKsmL+9pPbNf3Lqf1bfPQJJ9LsjzJzUk+lWT9qdQhSZLURVMZ9fmcGTj/esCPgE8DX5pgm9OAV/W8X9HX/jmah8HvSfPA+M/QPCVh/6FWKkmSNMsGejLBsFTVqcCpAEkm2mxFVd0wXkOSJwJ7A8+oqvPadX8NnJLkb6rquuFXLUmSNDuGOZhgpuyW5KYkP01yfJKH9bTtBNw8FtJaZwH3Ac+c1SolSZKGbKQ9aqvhNJpLor8AHgv8PXBqkp3a541uAtzUu0NV/TbJsrZtXEnWpZlOZMzCYRcuSZI0XZ0OalX1xZ63FyW5EPg5sBvwjWkc+lDg8GnsL0mSNOPmwqXP36mqK4BfAVu1q24ANu7dJsnawEPbtokcQTOtyNiy+dCLlSRJmqaBetSSbAY8myYkrRT2quofh1DXROfdHHgYcH27agmwUZLte6YGeW5b07kTHaeqVtAzenSSgQySJEkjM+WgluRA4OM0j4v6NVA9zQWsdlBr5zvbqmfVY5JsByxrl8OBk2l6xx4LfBC4HDgdoKp+0s6rdkKS19FMz3Es8EVHfEqSpLlukB619wPvA46oqvumef6nA//d837s6QcnAQcBTwEOADYCrgPOAA5re8TGvIImnH2DZrTnycAbplmXJEnSyA0S1B5M02M13ZBGVZ0NTHbdca/VOMYynNxWkiStgQYZTPAp4I+GXYgkSZJWNkiP2qHA15PsDVwE3NPbWFVvGUZhkiRJ892gQW0v4Kft+/7BBJIkSRqCQYLaW4G/qKoTh1yLJEmSegxyj9oK4JxhFyJJkqSVDRLUPgr89bALkSRJ0soGufS5A/DcJPsCP+b+gwn2G0ZhkiRJ890gQe1m4EtDrkOSJEl9phzUqupVM1GIJEmSVjbIPWqSJEmaBYM8lP0XTDJfWlX93rQqkiRJEjDYPWof6Xu/DvBUYG/g6OkWJEmSpMYg96h9dLz1SQ4Gnj7tiiQNzZaHLB76MZceuc/QjylJGt8w71E7FXjpEI8nSZI0rw0zqL0MWDbE40mSJM1rgwwm+CErDyYIsAnwCOD1Q6pLkiRp3htkMMGX+97fB/wSOLuqLp12RZIkSQIGG0zw3pkoRJIkSStzwltJkqSOWu0etST3MclEt62qqkEup0qSJKnPVELVSyZp2wl4A/bQSZIkDc1qB7Wq+kr/uiRbA0cCLwQ+B7x7eKVJkiTNbwP1gCXZLMkJwEU0YW+7qjqgqq4canWSJEnz2JSCWpINkxwFXA5sA+xeVS+sqotnpDpJkqR5bCqDCd4GvB24AfjT8S6FSpIkaXimMpjgSOBOmt60A5IcMN5GVbXfMAqTJEma76YS1D7LqqfnkCRJ0pBMZdTngTNYhyRJkvo475kkSVJH+RQBdd6Whywe+jGXHrnP0I8pSdKw2aMmSZLUUQY1SZKkjjKoSZIkdZT3qK2BvKdLkqQ1gz1qkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRIw1qSXZJ8rUk1yWpJC/ua0+S9yW5PsmdSc5K8ri+bR6a5HNJlie5Ocmnkqw/qx9EkiRpBoy6R2094EfAwRO0vw14A/A64JnA7cDpSR7Ys83ngG2APYF9gV2AT8xUwZIkSbNlpPOoVdWpwKkASVZqS7PiTcAHquor7bpXAjcCLwa+mOSJwN7AM6rqvHabvwZOSfI3VXXd7HwSSZKk4Rt1j9pkHgNsApw1tqKqbgHOBXZqV+0E3DwW0lpnAffR9MCNK8m6STYYW4CFwy5ekiRpuroc1DZpX2/sW39jT9smwE29jVX1W2BZzzbjORS4pWe5ZrrFSpIkDVuXg9pMOgLYsGfZfLTlSJIk3V+Xn/V5Q/u6CLi+Z/0i4IKebTbu3SnJ2sBDe/a/n6paAazo2Wf61UqSJA1Zl3vUfkETtnYfW9HeT/ZMYEm7agmwUZLte/Z7Ls3nOneW6pQkSZoRI+1Ra+c726pn1WOSbAcsq6qrknwEeFeSn9EEt/cD1wFfBqiqnyQ5DTghyeuAdYBjgS864lOSJM11o770+XTgv3veH9O+ngQcCHyQZq61TwAbAd8B9q6qu3r2eQVNOPsGzWjPk2nmXpMkSZrTRj2P2tnAhDeIVVUB726XibZZBuw/9OIkSZJGrMv3qEmSJM1rBjVJkqSOMqhJkiR1lEFNkiSpo0Y96lOS5rUtD1k89GMuPXKfoR9T0mjYoyZJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6qi1R12AJGnmbXnI4qEfc+mR+wz9mJJWZo+aJElSRxnUJEmSOsqgJkmS1FEGNUmSpI4yqEmSJHWUQU2SJKmjDGqSJEkdZVCTJEnqKCe8lSRJc9qwJ3Tu0mTO9qhJkiR1lEFNkiSpowxqkiRJHWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1lEFNkiSpozod1JK8J0n1LZf2tD8wyXFJfp3ktiQnJ1k0ypolSZKGpdNBrfVjYNOe5dk9bR8GXgj8EbArsBnwpdkuUJIkaSbMhYey/7aqbuhfmWRD4NXA/lX1zXbdq4CfJNmxqr43y3VKkiQN1VzoUXtckuuSXJHkc0m2aNdvD6wDnDW2YVVdClwF7DTZAZOsm2SDsQVYOFPFS5IkDarrQe1c4EBgb+Ag4DHA/yRZCGwC3F1VN/ftc2PbNplDgVt6lmuGV7IkSdJwdPrSZ1Wd2vP2wiTnAlcCfwzcOY1DHwEc0/N+IYY1SZKGZstDFg/9mEuP3Gfox+y6rveoraTtPbsM2Aq4AXhAko36NlvUtk12nBVVtXxsAW6dgXIlSZKmpdM9av2SrA88FvgX4HzgHmB34OS2fWtgC2DJqGqciL9ZSJKkqep0UEvyIeBrNJc7NwPeC9wLfKGqbknyKeCYJMuA5cDHgCWO+JQkSWuCTgc1YHPgC8DDgF8C3wF2rKpftu1vBu6j6VFbFzgdeP0I6pQkSRq6Tge1qnr5KtrvAg5uF0mSpDXKnBpMIEmSNJ8Y1CRJkjrKoCZJktRRBjVJkqSOMqhJkiR1VKdHfUqS5hYn95aGyx41SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUU54K2nanORUkmaGPWqSJEkdZVCTJEnqKIOaJElSRxnUJEmSOsqgJkmS1FEGNUmSpI4yqEmSJHWUQU2SJKmjDGqSJEkdZVCTJEnqKIOaJElSRxnUJEmSOsqgJkmS1FFrj7oASeqaLQ9ZPPRjLj1yn6EfU9Kazx41SZKkjjKoSZIkdZSXPiVJmke8tD+3GNSklv94SZK6xkufkiRJHWVQkyRJ6iiDmiRJUkd5j5okSR3gfbIajz1qkiRJHbXGBLUkBydZmuSuJOcm2WHUNUmSJE3HGhHUkvwJcAzwXuBpwI+A05NsPNLCJEmSpmGNCGrAW4ATquozVXUJ8DrgDuAvRluWJEnS4OZ8UEvyAGB74KyxdVV1X/t+p1HVJUmSNF1rwqjPhwMLgBv71t8IPGG8HZKsC6zbs2ohwPLly2eiPgDuW3HH0I85Ub2zdS7P0+3zzOa5PM9g55nNc3mewc4zm+fyPIOdZybONZN5YKrHT1XNYCkzL8lmwLXAzlW1pGf9B4Fdq+qZ4+zzHuDwWStSkiTp/javqmsn22BN6FH7FXAvsKhv/SLghgn2OYJm8EGvhwLLhlvaQBYC1wCbA7eOuJa5xu9uMH5vg/O7G5zf3WD83gbXte9uIXDdqjaa80Gtqu5Ocj6wO/BlgCRrte+PnWCfFcCKvtUz28+5mpKM/XhrVXWiprnC724wfm+D87sbnN/dYPzeBtfB7261apjzQa11DHBSkvOA7wNvAtYDPjPKoiRJkqZjjQhqVfVvSR4BvA/YBLgA2Luq+gcYSJIkzRlrRFADqKpjmeBS5xyzgmbi3v5Ls1o1v7vB+L0Nzu9ucH53g/F7G9yc/O7m/KhPSZKkNdWcn/BWkiRpTWVQkyRJ6iiDmiRJUkcZ1CRJkjrKoNYhSQ5OsjTJXUnOTbLDqGvquiQHJbkwyfJ2WZLk+aOua65I8sgk/5rk10nuTHJRkqePuq6uS7IwyUeSXNl+b99N8oxR19U1SXZJ8rUk1yWpJC/uaVsnyVHtn7nb220+2z4WcN6b7Ltr209s1/cup42o3E5Zje9u/STHJrmm/ft7SZLXjajcVTKodUSSP6GZuPe9wNOAHwGnJ9l4pIV13zXAIcD2wNOBbwJfSbLNSKuaA5I8BDgHuAd4PvAk4K3Ab0ZZ1xzxSWBP4M+BbYEzgLOSPHKkVXXPejT/lh08TtuDaf6te3/7uh+wNfDVWauu2yb77sacBmzas/zpLNQ1F6zquzsG2Bv4M+CJwEeAY5O8aFaqmyKn5+iIJOcCP6iq/9e+Xwu4GvhYVR050uLmmCTLgL+tqk+NupYuS3Ik8Kyq+oNR1zKXJHkQzXMC/7CqFvesPx84tareNbLiOixJAS+pqi9Pss0zaJ4u8+iqumq2auu68b67JCcCG1XVi0dU1pwwwXd3MfBvVfX+nnWd/ftrj1oHJHkATY/QWWPrquq+9v1Oo6prrkmyIMnLaX6bWjLqeuaAFwHnJfmPJDcl+WGS14y6qDlgbWABcFff+juBZ89+OWuUDYECbh5xHXPFbu3f3Z8mOT7Jw0Zd0BzxXeBF7a0fSfIc4PE0PeOdY1DrhofT/MPf/8irG2keiaVJJNk2yW00s03/M81vT5eMuKy54PeAg4CfAXsBxwP/mOSAkVbVcVV1K80vAocl2az9BeHPaH6p2nS01c1dSR4IHAV8oSMPzO6604BXArsDbwd2BU5NsmCkVc0Nfw1cQnPrzN003+XBVfXtkVY1gTXmEVKa134KbEfz2/jLgJOS7GpYW6W1gPOq6h3t+x8meTLwOuCk0ZU1J/w58GngWuBe4H+BL9D0jGuKkqwD/DsQml8etApV9cWetxcluRD4ObAb8I2RFDV3/DWwI81VhSuBXYDjklxXVWdNuucI2KPWDb+i+cd+Ud/6RcANs1/O3FJVd1fV5VV1flUdSnMT6RtHXdcccD3Nb5W9fgJsMYJa5pSq+nlV7QqsDzyqqnYA1gGuGG1lc09PSHs0sKe9aYOpqito/i/ZatS1dFl7j+nfA2+pqq9V1YXts8L/Dfib0VY3PoNaB1TV3cD5NF3YwO8GE+yO91oNYi1g3VEXMQecQzPKrtfjaX7D1Gqoqtur6vp2BO1ewFdGXdNc0hPSHgfsUVW/HnFJc1aSzYGH0fwCpomt0y739a2/l45mIi99dscxNJfszqMZ9fQmmpviPzPKorouyRHAqcBVwEJgf5qu/71GWNZc8WHgu0neQfOf5Q7Aa9tFk0iyF81lup/S9GAcDVyKf19XkmR9Vu7heUyS7YBlNIHiP2mm5tgXWJBk7J7cZe0vsPPWKr67ZcDhwMk0V10eC3wQuBw4fXYr7Z7JvruquirJt4Cjk9xJ84vprjT3+71l1otdHVXl0pEF+H80f2hWAOcCzxx1TV1fgE8BS9vv7CaakbJ7jrquubLQ/Ad5Ec0Ixp8Arxl1TXNhAf6Y5n6gFTSB41hgw1HX1bWF5pemGmc5EdhygrYCdht17aNeVvHdPYgmkN1EczP8UuATwKJR192FZbLvrm3fhOaXqmtpRmtfShPSMurax1ucR02SJKmjOnk9VpIkSQY1SZKkzjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkzTlJKsmLR13HbEpyYJKbR12HpNllUJPUKUk2SfKxJFckWZHk6iRfS7L7qvce6Hy7tcFvo5k4fnuOeRcsJQ2Hz/qU1BlJtqR5WPzNwN/SPN5qHZpntx4HPGFUta1KkgALquq3o65F0prDHjVJXfJPNM/k26GqTq6qy6rqx1V1DLDjeDuM1yOWZLt23Zbt+0e3vXK/SXJ7kh8neUHb/t/tbr9p9zmx3WetJIcm+UWSO5P8KMnLxjnv85OcT/Pcz2ev6gMm2bLdb78k/53kjvbYO/Vtd2CSq9r2/wIeNs6x/jDJ/ya5q+2BPDzJ2m3bu5Ncl+RhPdsvbs/pv/3SHOFfVkmdkOShwN7AcVV1e397Vd08jcMfB6wL7AJsC7wduA24Gnhpu83WwKbAG9v3hwKvBF4HbAN8GPjXJLv2HftI4BDgicCFU6jp74APAdsBlwFf6AlZzwQ+RfOw9+1owuS7endO8gfAZ4GPAk8C/go4EHhnz/GXAp9stz8Y2Bk4oKrum0KdkkbIS5+SumIrIMClM3DsLYCTq+qi9v0VYw1JlrU/3jQWBpOsC7wD2KOqloztk+TZNIHoWz3HfndVnTlATR+qqsXt+Q4HfkzzHVxKExZPq6oPtttelmRnmiA75nDgyKo6qae+w4APAu+tqnuT/BlwQZIjgTcAf1lVVw1Qq6QRMahJ6orM4LH/ETg+yfOAs2hC22S9X1sBDwbObG49+50HAD/s2/a8AWvqPf/17evGNEHticB/9W2/hJWD2u8Dz0ryzp51C4AHJnlwVd1RVVck+Rvg48C/VdXnB6xV0ogY1CR1xc9o7k+b6oCBsct4vYlqnd4NquqTSU4H9gGeBxya5K1V9bEJjrl++7oPcG1f24q+9/e7TLua7uktsX2dyu0o69P0qn1pnLa7en7eBbgX2DLJ2g52kOYW71GT1AlVtQw4HTg4yXr97ZNMn/HL9nXTnnXbjXP8q6vqn6tqP+AfgNe0TXe3rwt6Nr+EJpBtUVWX9y1Xr+5nmoafAM/sW9c/mOJ/ga3Hqe/ysXvQkvwJsB+wG83l38NmuG5JQ2aPmqQuOZhmeo7vJ3k3zeXBtYE9gYNoLgn2u5xmUMB72suAjwfe2rtBko8Ap9LctP8Q4Dk0YQjgSpoerX2TnALcWVW3JvkQ8OF2hOR3gA2BZwHLe+4Lmyn/CJzTXrb8Cs30JHv3bfM+4OtJrgL+k6Zn8feBJ1fVu5JsDhwPvL2qvpPkVe32p1bV92a4fklDYo+apM6oqiuAp9GMcvwH4GLgTGB3mqA23j73AH9Kc8n0QpoRne/q22wBzcjPnwCn0QS217f7X0t7Yz5wI81IS2h6n95PM/pzbL99gF9M+4OuQhukXkMzqOBHNJdrP9C3zenAvm3bD4DvAW8GrmzndDsR+D7t52m3P55m5Or6SJoTUlWr3kqSJEmzzh41SZKkjjKoSZIkdZRBTZIkqaMMapIkSR1lUJMkSeoog5okSVJHGdQkSZI6yqAmSZLUUQY1SZKkjjKoSZIkdZRBTZIkqaMMapIkSR31/wEwVa3cfmTjRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.ticker import MaxNLocator\n",
    "ax = plt.figure(figsize=(7,5), dpi=100).gca()\n",
    "ax.bar(np.arange(NumClusters), np.array(clusterLens))\n",
    "ax.set_title('Cluster Size')\n",
    "ax.set_xlabel('Cluster Index')\n",
    "ax.set_ylabel('Num Users in Cluster')\n",
    "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]\n",
      " [0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]\n",
      " [0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]\n",
      " ...\n",
      " [0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]\n",
      " [0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]\n",
      " [0.00729927 0.00729927 0.00729927 ... 0.00729927 0.00729927 0.00729927]]\n"
     ]
    }
   ],
   "source": [
    "#QUICK CODE FOR DEMO\n",
    "cluster1History = roundedMeans[clusterIdx].astype(int)\n",
    "TransMat = MakeTransitionMatrix(cluster1History, nBuildings-1)\n",
    "print(TransMat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([99, 61, 112, 18, 58, 8, 83, 89, 124, 16],\n",
       " ['Phelps Lab',\n",
       "  'The Hub',\n",
       "  'Research Laboratory 2',\n",
       "  'Bryan Hall',\n",
       "  'Hough Hall Business School',\n",
       "  'School of Architecture',\n",
       "  'Mehrhof Hall',\n",
       "  'Physics',\n",
       "  'SW Recreation Center',\n",
       "  'Broward Hall'])"
      ]
     },
     "execution_count": 260,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Predict future steps based on their most recent history\n",
    "NumFutureSteps  = 10\n",
    "MostRecentState = cluster1History[-1]\n",
    "bldgs = location_info['name'].unique()\n",
    "bldg_dict = dict(zip(np.arange(len(bldgs)), bldgs))\n",
    "nextSteps, nextBldgs = PredictFutureSteps(TransMat, MostRecentState, bldg_dict, NumFutureSteps, PlotCDF=False)\n",
    "nextSteps, nextBldgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[99, 61, 112, 18, 58, 8, 83, 89, 124, 16]"
      ]
     },
     "execution_count": 261,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open('cluster1History.npy', 'wb') as f:\n",
    "    np.save(f, cluster1History)\n",
    "with open('nextSteps.npy', 'wb') as f:\n",
    "    np.save(f, nextSteps)\n",
    "    \n",
    "cluster1History\n",
    "nextSteps"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "metadata": {
   "interpreter": {
    "hash": "d262b3cf53436df4d43148531b8b8c13029e50250a57d118b8471bcdaf78eff0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}