Interactive_DBSCAN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a6885e33",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Interactive Density-based Clustering with DBSCAN \n",
    "\n",
    "### Michael J. Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60d9f326",
   "metadata": {},
   "source": [
    "#### DBSCAN Clustering\n",
    "\n",
    "The DBSCAN (Density-Based Spatial Clustering of Applications with Noise) approach finds clusters of samples (high density) is well suited to working with data with an arbitrary (even non-convex) shape and with noise.\n",
    "\n",
    "For a complete lecture with linked Python workflows check out:\n",
    "\n",
    "* [Density-based Clustering Lecture](https://www.youtube.com/watch?v=3GaLe8HaDMc&list=PLG19vXLQHvSC2ZKFIkgVpI9fCjkN38kwf&index=15)\n",
    "\n",
    "Here's a complete Python workflow and more details DBSCAN:\n",
    "\n",
    "* [DBSCAN workflow](https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/SubsurfaceDataAnalytics_advanced_clustering.ipynb)\n",
    "\n",
    "The two DBSCAN hyperparameters include:\n",
    "\n",
    "* **eps** - maximum distance between to samples in features space for a sample to be included in the neighbourhood of the other sample. Too small a value will result in many clusters and outliers, and too large a value will result in one cluster with all the data together.\n",
    "\n",
    "* **min_samples** - minimum number of samples in a neighbourhood to spawn a new cluster. Use a larger value for noisy data to prevent identifying clusters due to a few outliers.\n",
    "\n",
    "The method proceeds through the dataset and assigns the samples as:\n",
    "\n",
    "* **core** point if there are $\\ge$ **min_samples** within **eps** distance from the sample\n",
    "\n",
    "* **border** point if there are $lt$ **min_samples** within **eps** distance from the sample, but the sample is within **eps** distance of a core point.\n",
    "\n",
    "* **outlier** point if there are $lt$ **min_samples** within **eps** distance from the sample and the sample is not within **eps** distance of a core point\n",
    "\n",
    "Once the points are assigned these labels, all connected core points and their associate border points are assigned to an unique cluster.  Outlier points are left as outliers without an assigned cluster.  To understand the cluster assigments we should explain the following forms of connection.\n",
    "\n",
    "**directly density reachable** - point X is directly density reachable from A, if A is a core point and X belongs to the neighborhood (distance $le$ eps) from A.\n",
    "\n",
    "**density-reachable** - point Y is density reachable from A if Y belongs to a neighborhood of a core point that can reached from A. This would require a chain of core points each belonging the previous core points and the last core point including point Y.\n",
    "\n",
    "**density-connected** - points A and B are density-connected if there is a point Z that is density-reachable from both points A and B.\n",
    "\n",
    "**density-based cluster** - a nonempty set where all points are density-connected to eachother. \n",
    "\n",
    "The approach iterates as follows:\n",
    "\n",
    "1. all points are labled as unvisited\n",
    "\n",
    "2. randomly visit an unvisited sample\n",
    "\n",
    "3. check if a core point ($ge$ min_sample within eps distance), if so label as core otherwise label as outlier\n",
    "\n",
    "4. now visit all points within eps distance of the core point, determine if core, otherwise label as border point\n",
    "\n",
    "5. recusive operation where all points within eps distance of new core points are checked\n",
    "\n",
    "6. once this is exhausted then randomly visit an unvisited point\n",
    "\n",
    "This apporach may be thought of as identify and grow/merge clusters guided by local point density.\n",
    "\n",
    "\n",
    "After some careful interations of these parameters we get the following result.\n",
    "\n",
    "#### Load and Configure the Required Libraries\n",
    "\n",
    "The following code loads the required libraries and sets a plotting default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ecc3ed57",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "supress_warnings = False\n",
    "import os                                               # to set current working directory \n",
    "import sys                                              # supress output to screen for interactive variogram modeling\n",
    "import numpy as np                                      # arrays and matrix math\n",
    "import pandas as pd                                     # DataFrames\n",
    "import matplotlib.pyplot as plt                         # plotting\n",
    "from sklearn.model_selection import train_test_split    # train and test split\n",
    "from sklearn import tree                                # tree program from scikit learn (package for machine learning)\n",
    "from sklearn import metrics                             # measures to check our models\n",
    "import scipy.spatial as spatial                         #search for neighbours\n",
    "from matplotlib.patches import Rectangle                # build a custom legend\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "import math                                             # sqrt operator\n",
    "from ipywidgets import interactive                      # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "cmap = plt.cm.inferno                                   # default color bar, no bias and friendly for color vision defeciency\n",
    "plt.rc('axes', axisbelow=True)                          # grid behind plotting elements\n",
    "if supress_warnings == True:\n",
    "    import warnings                                     # supress any warnings for this demonstration\n",
    "    warnings.filterwarnings('ignore')                  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaacbb36",
   "metadata": {},
   "source": [
    "#### Declare Functions\n",
    "\n",
    "The following functions for clean code. \n",
    "\n",
    "* Just a improved grid for the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "64952c7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_grid():\n",
    "    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36bf017",
   "metadata": {},
   "source": [
    "#### Recursive Method to Perform DBSCAN \n",
    "\n",
    "I wrote this recursive method to perform DBSCAN. Then I used it in the dashboard below.\n",
    "\n",
    "* Recusursion can be tricky to code so I left my debugging outputs below to demonstrate how this was made and checked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "088d9559",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "at point:5, step: 0\n",
      "[]\n",
      "step: 1\n",
      "at point:37, step: 1\n",
      "[43, 30, 42]\n",
      "at point:43, step: 2\n",
      "[37, 30, 14, 42, 62]\n",
      "at point:30, step: 3\n",
      "[37, 43, 76, 74]\n",
      "at point:76, step: 4\n",
      "[30, 74]\n",
      "step: 5\n",
      "at point:74, step: 5\n",
      "[30, 76, 6]\n",
      "at point:6, step: 6\n",
      "[94, 74]\n",
      "step: 7\n",
      "step: 7\n",
      "step: 7\n",
      "at point:14, step: 7\n",
      "[67, 43, 42, 62]\n",
      "at point:67, step: 8\n",
      "[17, 89, 14]\n",
      "at point:17, step: 9\n",
      "[79, 97, 67]\n",
      "at point:79, step: 10\n",
      "[66, 17]\n",
      "step: 11\n",
      "at point:97, step: 11\n",
      "[17, 38, 15, 89, 41]\n",
      "at point:38, step: 12\n",
      "[91, 97, 15, 41]\n",
      "at point:91, step: 13\n",
      "[29, 13, 20, 23, 38]\n",
      "at point:29, step: 14\n",
      "[77, 13, 20, 91]\n",
      "at point:77, step: 15\n",
      "[29, 13]\n",
      "step: 16\n",
      "at point:13, step: 16\n",
      "[77, 29, 20, 91]\n",
      "at point:20, step: 17\n",
      "[29, 13, 91, 23]\n",
      "at point:23, step: 18\n",
      "[20, 91]\n",
      "step: 19\n",
      "step: 19\n",
      "step: 19\n",
      "step: 19\n",
      "step: 19\n",
      "at point:15, step: 19\n",
      "[97, 38, 41, 86]\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n",
      "step: 20\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAF9CAYAAAAOfLFQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABZp0lEQVR4nO3deVxV1f74/9diOoCAoiKoYBo4hwOmacaNUsoKse79ebNJK6tPRZQZmQ037/10y65hwzWab2Vl+bU+3USuWlqiWKKW8zymOKDiBAIyrt8fnMMFBGU40z7n/Xw8zgP2Pvus8z7rDO+99l57LaW1RgghhBCuz8PRAQghhBDCPiTpCyGEEG5Ckr4QQgjhJiTpCyGEEG5Ckr4QQgjhJiTpCyGEEG5Ckr4QokmUUluVUnGOjkMI0XSS9IVwY0qp35VSxUqpc0qpY0qpT5RSARd7jNa6r9Y6swnlj7RKsEKIFpOkL4QYrbUOAGKAwcALDo5HCGEjkvSFEABorQ8Di4ArlFKJ5sP4Z5RSmUqp3pbtarbelVJ/VUrNU0p9ppQqMD/mSvN9nwNdgAXmIwlTlFK+SqkvlFInzWWvVUqFOuL1CuGOJOkLIQBQSkUANwMFwFfAJCAEWEhV4vZp4KGJwFygDZAOvA2gtb4HOIj5SILWegYwAWgNRADtgIeBYtu8IiFEXZL0hRDfKaXOACuB5cA24D9a6yVa6zIgFfADrm7g8Su11gu11hXA50D/izxXGVXJPkprXaG1/k1rnW+tFyKEuDhJ+kKIW7XWbbTWl2mtHwU6AQcsd2qtK4EcoHMDj8+t8X8R4KuU8mpg28+B74G5SqkjSqkZSinvlr8EIURjSNIXQtR1BLjMsqCUUlQdjj/cjLJqTeOptS7TWv9Na92HqiMHCcD4FsQqhGgCSfpCiLrmAbcopUaYW+FPASXAL80o6xhwuWVBKXWdUipaKeUJ5FN1uL/CCjELIRpBkr4Qohat9U7gbmAWkAeMpqozXmkzipsOvGDuqZ8ChAHfUJXwt1PVh+ALqwQuhLgkpbW+9FZCCCGEMDxp6QshhBBuwmZJXyn1sVLquFJqS411bZVSS5RSu81/g2vc96xSao9SaqdS6sYa6wcppTab7/unuVOREEIIIZrIli39T4FRddZNBX7UWncHfjQvo5TqA4wD+pof8465ow/Au8BDQHfzrW6ZQgghhGgEmyV9rfUK4FSd1WOA2eb/ZwO31lg/V2tdorXeD+wBhiilOgJBWutVuqrzwWc1HiOEEEKIJrD3Of1QrfVRAPPfDub1naka/MPikHldZ/P/ddcLIYQQookaGjXL3uo7T68vsv7CApR6iKrTABfw9fUd1Lmz7CsIIYQwlr179+ZprUOsVZ69k/4xpVRHrfVR86H74+b1h6ga8csinKpRwQ6Z/6+7/gJa6w+AD+q7LyoqSu/Zs6elsddSXFw1R4ifn59Tl2nLcgHmz5/PmDFjrFaekerASPUKxqoDW5VrlHq1VblG+swaqV5tWa5S6sClt2o8ex/eT6dqli3Mf+fXWD9OKWVSSnWjqsPeGvMpgAKl1FBzr/3xNR4jhBBCiCawWUtfKfUVEAe0V0odAqYBrwLzlFITqZpycyyA1nqrUmoeVbN7lQNJ5hm7AB6h6koAP6rm+l5kq5iFEEIIV2azpK+1vqOBu0Y0sP3LwMv1rP8VuMKKoQkhhDCAjIwMVq1aRVxcHImJiY4OxyXIiHxCCCGcTkZGBikpKYSFhZGcnEx6erqjQ3IJztJ7XziJ9PR0MjMzZc9aCOFQWVlZJCUlMWXKFAAyMzPlN8kKpKXv4tLT03nmmWfIyMho1LbJyclN2rNes2YNkydPlr1wIRopPT1dvjONEBsbS1paGjNmzCAtLY24uDhHh+QSJOnbkKO/3JYkHh4eTkpKyiXjyMzMrN6zTkpKIjMz85Llz549Ww6/CdFIzdmxdlcJCQmkpqaSm5vLrFmzpJVvJXJ430YsX+6kpCSSk5MB7P6hrZnELcsXiyEuLq461rS0NGbNmnXJ8h9//HE5/CZEIzX1O+nuEhISGDt2rKPDcCmS9G3EGb7cTU3ilvgyMzMbtWcdFxfHo48+2ujyhXB3Tf1OCmFtkvQvoiWd2pzhy22JeenSpaSmpjbqNSQmJjb6tSYmJrJ69Wo5/CYazd07ijZ1x1oIa5Ok34CWHp53li93YmIi8fHxNit/yJAhVh/WVLgmZzjl5QyasmMthLVJ0m+ANQ7Py5e7fu7e2nNXznDKSwh35zZJ3zIZQmMNGzaMlJQUoOrwfGpqaq0yzp8/b9X4bFWmLcu1aErdWgbcsLT2SkpKSEhIqL7fSHXgTPXaGI6ug0t9p5pbblMZoV5tVa6RPrNGqldblmttLpn0lVKjgdGW5bCwsCaXYUlEWVlZpKam1kpMovnqDriRlZUldesm5DslhOO5ZNLXWi8AFliWo6KiHmzOdIdjx4695OUitpii0hZlOku5I0eOvKCDY32Pd4ZYHVmm0cptbJmN+U41p9ymMFK92qpcidV45VqLSyZ9V+GK576dpYOjEEK4I0n6TsqVezpLB0chhHAMSfpOSno6CyGEsDZJ+k7KGQb3cSRXPLUhhBCOJhPuOKnExERmzZrllqPdZWRkyKQkQghhA9LSd2Lueu5b5tEWQgjbkKQvnE5sbGytQVzc7dSGEVhOvwwbNkyutxfCQOTwvnA6CQkJbntqwwhqzgmfkpJCRkaGo0NyCunp6UyePFlORwmnJi194ZTc9dSGEdS9siQrK8vt5zyvO7w0uM4ltsK1SNIXQjRJ3StLUlNTHRyR40k/FGEUkvSFEE1Sc1RFGUO/ivRDEUYhSV8I0WSW0y/WnrHOqBISEjCZTDK8tHB6kvSFEMIKpB+KMAK3SfpGmEPbqPM8u+sc2kaqVzBWHdiybo1Qr7Yq10ifWSPVqy3LtTaXTPpKqdHAaMtyWFiYA6MRQgghnINLJn2t9QJggWU5KirqQSPNnWykWG1VrsRqrHIlVmOVK7Ear1xrkcF5hBBCCDchSV8IIYRwE5L0hRBCiBZIT0/nmWeeMcSQ1JL0hRBWJWPQC3dimYsiPDyclJQUp//cS9IXQlhNzcl4kpOTnf4HUIiWqjkXRVJSEpmZmY4O6aJcsve+EMIx6k7GI2PQC1dXdy4KZx+CWZK+EM106NAh3n57Jtu3ryM6+iomTZpCRESEo8NyKKP9AArRUpad2qVLl5Kamur0O7mS9IVohpycHOLihjBuXAGTJ1ewZMk6hg6dQ3b2RrdO/DUn45Ex6B0rPT2dzMxM4uLi5H2wscTEROLj4x0dRqNI0heiGd58cwbjxhUwc2YFAKNGVQD5vPnmDGbOdO/WrYxB73iWvhVJSUnVR17kPREgHfmEaJbNm1cTH19Ra118fAWbN69xUERC/JfROpcJ+5GWvhDNEB19FUuWrDO38KssWeJJdPQQB0YlRBXpWyEaIklfiGZ4+OHHGT78E6CY+PhKlizxZO7cILKzpzg6NCGkb4VokCR9IZpIa83Ro0f5+OO5fPrpB/zjH4eIiRlOdrb03hfOQ/pWiPq4TdI3whzaljIzMjLIysoiNjaWhIQEq5VrK+42h3ZOTg6nT59m+PDhFBcXExcXR0BAAGDdujDSZ9Yo5YIx6tVW5cpvgfHKtTaX7MinlBqtlPrAcissLHR0SI22aNEiUlJSqod0NMJYzu6kpKSEjRs3MmjQIDw8PDh//rzTT6UphBAWLtnS11ovABZYlqOioh40ytzJ2dnZtUY0W7VqFWPHjrVK2UapA1uVaY1y169fT2RkJOHh4Zw8eRI/P7/qVr61OWsd2KtMW5VrpFhtVa7E6lzl2nNMBZds6RtZbGwsaWlpzJgxg7S0NOLi4hwdkjA7fPgweXl59OvXD6g6lCmtfCFES9h7vgqXbOkbWUJCAiaTSXrdOpnS0lLWrl3LsGHD8PKq+toUFRXRqlUrB0cmhDAye89XIUnfCUmvW+ezfv16OnXqRGhoaPW6oqIiaekLIVrE3mMqSNIX4hJyc3PJzc3l5ptvrrW+qKiIwMBAB0UlhHAF9h5TQZK+EBdRXl7OmjVrGDx4MN7e3rXuKy4urtXyF0KI5rDn0V3pyCfERWzcuJGQkBA6dep0wX1yeF8IYTSS9IVowIkTJzh48CAxMTH13l9UVIS/v7+doxJCiOaTpC9EPSoqKli9ejWDBg3CZDJdcH9ZWRmVlZX4+Pg4IDohhGgeSfpC1GPLli20bt26wbH0LYf2lVJ2jkwIIZpPkr4QdZw+fZq9e/dy5ZVXNpjUCwsL5dC+EMJwJOkLUUNlZSXZ2dkMGDDgop305Hy+EMKIHJL0lVJPKqW2KqW2KKW+Ukr5KqXaKqWWKKV2m/8G19j+WaXUHqXUTqXUjY6IWbiH7du34+fnR7du3S66nbT0hRBGZPekr5TqDDwOXKm1vgLwBMYBU4EftdbdgR/Nyyil+pjv7wuMAt5RSnnaO27h+s6ePcvOnTsZPHjwJc/VS0tfCGFEjjq87wX4KaW8AH/gCDAGmG2+fzZwq/n/McBcrXWJ1no/sAcYYt9whavTWrN69Wqio6MbNZ5+YWGhXKMvhDAcpbW2/5Mq9QTwMlAM/KC1vkspdUZr3abGNqe11sFKqbeBbK31F+b1/wIWaa2/qVPmQ8BD9T1fSEjIoA8//NBGr0YYmdYapRRnz56lsLCQjh07NqpH/sGDB+nYseMFo/QJ0Rxr1qxh69at9O3blyFDpE3jzOz9Xt16662/aa2vtFqBWmu73oBg4CcgBPAGvgPuBs7U2e60+W8acHeN9f8C/tSU54yMjNTWVlRUpIuKipy+TFuWq7XW3333nVXLs1cdBAcGauCCW3Bg4CXLqqys1HPnztUFBQWGqVet5TOrtXPW6/z583WXLl30P/7xD92lSxc9f/58q5RbH/ktaFm5Db1XLS33YoBftRVzsCPG3h8J7NdanwBQSn0LXA0cU0p11FofVUp1BI6btz8E1LxYOpyq0wFCNNvpggLqO8alCgou+dji4mJ8fHzw9JSuJaLl7D21qmg+V3ivHJH0DwJDlVL+VB3eHwH8ChQCE4BXzX/nm7dPB75USr0OdAK6A2vsHbRwHm2DgjhdT3IODgzkVH5+vY/RWlNaWsrJkyc5fvx4vdtYHD16FC8vLyorK9FaU1lZSWVlJdG9enGmsLBJzyvEpdh7alXRfK7wXtk96WutVyulvgHWAeXAeuADIACYp5SaSNWOwVjz9luVUvOAbebtk7TWFfaOWziPi7XSg4ODad++PcHBwQQHB9O6dWsCAwPx9/fHZDJRUlJCUVHRRct//fXXKSsro7CwkIKCAs6cOcOpU6c4U1jY7KMDQjTE3lOriuZzhffKIVPraq2nAdPqrC6hqtVf3/YvU9XxT4iLevrpp/Hy8sLLywsfHx9MJhMmkwlfX198fX2rO+l9/PHHDZYxZMgQzp8/T3FxMSUlJZSWllJWVsbatWsbfMyWLVu44oorrP56hHuw59SqomWM/l45JOkLcSlaa4qLi8nPz+fs2bPk5+eTn59PwSVa1X369MHDw6M6uSulqm+N5e/v3+Rr8EeMGMHw4cN54YUXGpyVTwghHE2SvrCqnJwc3nxzBps3ryY6+iomTZrS4KQ1UDWbXUFBQXVSP3nyJAUFBRQXF+Pt7U1gYCBBQUEEBQURHh5OUFDQRZ//zjvvrD4HX/N8PFTtAGit+frrr2llMqFKSi54vDcwfvz46h0HDw8PPDw86Nq160Wf98MPP2Tv3r0kJiYyYMAAXnjhBYYOHXrJ+rKH5vSBEEK4Jkn6wmpycnIYOrQ/48blM3lyBUuWrGPo0DmsWrWBDh06VCf2mrfi4mICAgKqE3toaChRUVGEhIQ0a9rakydPXpCwlVKcP3+ehx9+GF9fX06dOsW8b79l4sSJzJs3j9jY2AbLKy8vZ/Hixezdu5dVq1Y1uF1xcTFKKd5++20OHTrEHXfcQVRUFC+88ALXXnttk1+HNbXkSgUhhGuRpC+s5s03ZzBuXD4zZ1b1sxw1qgKt80lO/h/uuWdidWIPCgqiQ4cOBAUFERAQgIfHfweGLC4uBrhowg8ODKw3YQUHBjY4St7ChQvJysrirbfe4sSJE3z66afcddddF034x44dY/78+Xh7ezN+/Hj+9pe/1Pu8Jg8PXn75ZZ566inKysooLy9n5syZ5OXl8cADD9CxY0deeOEF4uPj7ToVb1FREV9//bXdnk8I4fwk6Qur2bx5NZMn176w4oYbKkhNzeNPf/qT1RLeqfx8HnzwQQYOHEhZWRk7duwgOjqaRx99tN7tT548yWOPPcZ7773H3r17CQ4OZsuWLXz22WcNPkd2djbZ2dlER0dz3XXX4eHhccGhcMsOiq+vL4sXL+all17i1KlTPPXUU0BV0p0+fToFBQU8+eSTBAQE8MILL5CQkGDT5L9hwwY+/PBD5s6dy7Bhw2z2PEII45GpdYXVREdfxZIltQes+eEHTy6/vK/Vnys0NJQTJ05QXl5OaGgox44da3Dbxx57jHHjxnHkyBF69erF008/zezZs/H19b1g26KiIr766it+++03br31VkaMGFHrSER9lFLcdNNN/Pzzz7zzzjt8+eWXPP3005hMJjw9PTl79izTpk0jOTmZF198kYEDB/LNN99U9zWwhoKCAj788EOGDBnCmDFj6NChA1988QVjxoyx2nMI55Cens7kyZNJT093dCjCgCTpC6uZNGkKc+cG8dRTnixeDE895cncuQEMH34dP/74I6dPn7bac1mSvqenJ+3atWsw6c+bN48NGzZw7bXX4u/vT1paGg8//DCDBw++YNvdu3fz8ccf4+XlxcSJEy/Zea8upRTXX389y5YtY/bs2cyfP58nnngCDw8PTCYTeXl5TJkyhalTp/Laa68RHR3Nl19+SXl5eXOqAK01a9eu5cEHHyQiIoJFixaRkpLC66+/TuvWrTl27BghISGXLEMYR3p6OsnJyYSFhZGcnCyJXzSZJH1hNREREWRnbwQe4fXXhwCPsHr1ZsaPH0/Xrl3JzMxk9erV1YfFW8LSuvfx8aFNmzb1Jv3c3FySk5N57bXXOHDgACUlJRw+fJgXXnih1naVlZUsXryYhQsXMnToUMaOHVvvUYCmuOaaa1i8eDHffPMNP/74Iw8//DBlZWUEBARw7NgxkpOT+etf/8p7771H7969+eSTTygrK2tU2WfOnCEtLY2BAwdyxx130K1bN7766itGjRpFbm4u3t7etG7dmi+//JIff/wRfx8fFFxwC/TzIzs7u9k7HcL+ag4Dm5SURGZmpqNDEgYjSV9YVUREBDNnzuKHH1Yzc+YsIiIiUEoRFRXFLbfcgo+PDwsXLmTbtm1UVDR/YEVL0jeZTPj7+1+Q9LXW/M///A8PPvggBw4coGfPnjz//PPMnj27VifBvLw8Pv74Yw4fPsxdd93FlVdabzIrqBroZ/78+SxatIi1a9cyceJECswjBx45coQHHniA1157jS+//JLu3bvz3nvvUVLPpYRaa37++WcmTJhA165dycrK4oUXXmD69On4+/tz4sQJOnfuzL59+7j//vv59ttvuf766xkxYgQFxcVorSkqKqKgoIBvv/22apTB/Hy01vzwww/ky6V7hhAXF0daWhozZswgLS2NuLg4R4ckDMZtOvJZo3VZ0/nz561anq3KtGW5Fk2p2169ehEeHs7mzZurO+B17ty5umNbY2O1HL729fXF29ub3NzcWnF8/vnnHDhwgAkTJlBcXMxrr71GUlISUVFR1dtt2LCB1atX0717d+Li4vDw8GjSa2lKvfbs2ZPZs2ezfft2UlNTefPNN0lKSqJ79+7s3buXu+66i7Zt2/LBBx/w+KOPUlbPYXcfDw/+/uqr3H777Rw8eJCjR4/SrVs3WrVqxXvvvcf+/fsZP34833//PatXr6Zjx47Ex8dX70RY4u3SpQtbt25l4MCBDBgwgP3797N48WJiYmLo3Llzo19TU+vAGcoFY/wWNFRufHw8qampZGVlkZqaSnx8vM0+s81hzbq1ZqwZGRlkZWURGxvLyJEjrVZuTbauW2txyZa+Umq0UuoDy62wnklShOMEBAQwbNgwBg0axI4dO1ixYkWTz/d36NCBY8eO4e/vj6enZ61JdHJycnj++eeZPn06hw4dIi8vj/z8fB577DGg6ss5f/58fv31V0aOHMn1119/yc561tK7d2/+9a9/kZmZye+//87dd9/N4cOHCQ0N5cCBA9x2222UaX3hnL9AaWUlXl5enDp1is6dO7N9+3bGjx/Pv//9b5588kl27drF5MmTqxP+DTfcUO/r6tatGzk5OZSWlqKU4vLLL2f48OFs2rSJjRs3WrWDobC+hIQE/vGPf5CQkODoUAwhIyODlJQUwsPDSUlJYdGiRY4OyaFcsqWvtV4ALLAsR0VFPdjQ9dstZYtyjRRrS8q97LLL6NKlC3v37mX16tV06tSJ7t274+fnd8kyfX19KS4uxmQyUVZWVj3ynp+fH0lJSTzxxBPs2bOHnj17ctddd/HTTz/RqlUrjh07xn/+8x/atWvHxIkTmzzcbn2a8/qvuOIKZs+ezYEDB/jHP/7BnXfeyb333ssNN9xw0ceZTCbefPNNTpw4wcSJE9m0aRPh4eEA5Ofn8+9//5vw8HASEhIa3JHx8/MjIiKCI0eO0Lt3bwA6d+7MLbfcwqpVq/jll18YPnx4k+rG2T5b9i7TaOW6U6yrVq2qNR1udnY2t912m6HqwJpcsqUvjKPm+X6TycTSpUvZsWPHJc/3K6Xo0KFD9Rj9lnP87733HgUFBfTo0YOQkBD+9re/MW3aNCIjI1m+fDnz589n0KBBjBs3zioJv6Uuu+wy3nnnHbZu3QpUDSN8MUuXLuXvf/87+/bt48UXX6yV8L/88ks6dep00YRv0bNnT3bt2lWr977JZOLaa6+lU6dOfP/99+Tm5rbw1QnheHX7QVxsQC534JItfWE8Pj4+DBgwgM6dO7N582YyMjIYOHBgdUfA+oSGhlJaWlqd9H/77TdefPFF5syZw65duzh37hy+vr7ceeedfPXVV1RUVDBu3DhCQ0Pt/OourVOnTrzxxhs8++yzF43vm2++uSChNzXhA7Rr1w4/Pz8OHz5cveMAVTtTffv2pV27dqxatYoePXrQp08fu44kKIQ11Z0ONz4+3sEROZYkfeFULOf78/PzWbduHbt27SImJoa2bdtesG1oaChFRUWUlJQQGhrKiy++yHPPPceuXbvo3r0799xzD1999RVffPEFXbt25frrrycwMNABr6rxOnTocNH7FyxYQK9evYiMjMTLy6tZCd+iZ8+e7Ny5s1bStwgLC+PGG2/k559/5sSJEwwbNgyTydTk1yOEM6g5Ha61O3IajRzeF04pNDSUUaNG0a1bN5YvX052dvYFX9bQ0FAKCgooLS0lNDSUsLAwwsLCqse6f/rpp9m+fTs33ngjN9xwA15ezr+Pe+7cOby48Jp6RdXcAtdccw3Hjx8nPT2dtWvX8sUXXzQr4UPV5ZUFBQWcOXOm3vv9/f0ZMWIEQUFBLF68mJMnT7botQkhHE+SvnBaSikiIyNJSEjA19eXhQsXsnXr1urz/aGhoeTl5VFWVkZ4eDhTp07l9OnTbNq0iREjRtCpUyfuu+8+evXq5eBX0njZ2dkMveYaPv74Y+bPn0+PHj3QWqO15lR+Pu3atSM2NparrrqKFStWUFZWRkRERL3X9l+Kh4cHUVFR7Ny586LbxMTEEBMTw/Lly9m9e7eM4ieEgTl/08dJpKenk5mZSVxcXPVhImEf3t7eDBgwgMjISDZs2FB9vn/WG29QWE+y8/f25pv587nxxhvtdimetaxcuZJrrrmGM2fOcO7cuXo7HeXn57NgwQK6d+/Oddddx86dO1m4cCFdunShd+/eBAQENPr5oqKiyMjIYMCAARc9fB8REUHr1q1ZuXIlJ06cYMiQIYY4ciKEqM1Yv4gOIuNdO4fAwEBiY2MZOnQoW7dupbCkpN7r2YvKyrjpppsMl/ABsrKyiI6OxmQy8csvv3DNNdfUur/uOfzAwECuvPLK6tEOv//+e1atWsXZs2cb9Xy+vr6Eh4ezd+/eS24bFBTEDTfcgFKK77//XkbxE8KAjPer6AAy3rVzsZzvdzVlZWWsWbOGNm3aEBISUt3qt7Ak/I4dO15wDt/X15f+/fszevRogoKC+Omnn1ixYgV5eXmXfN6ePXuye/fuRg3K4+XlxdChQ+nZsydLly7l0KFDzXuxQgiHkONzjRAXF0dycjIAaWlpzJo1y8ERCVe8hGzDhg1069aNU6dOERwcTHl5OZGRkUDthD969OgGj2L4+PjQt29fevXqxd69e/n5558JCAigb9++BAUF1VtvwcHBtGrVikOHDtGlS5dLxmkZW6Ft27ZkZmaSl5fHVVddZcgjK0K4G0n6jVD3Os/ExES3v+xDWF9WVhbXXHMNJ0+epKSkhGuuuQalVKMTfk2enp706NGDqKgoDhw4wG+//QZUteojIyMvSP49evRg586djUr6Fm3btmXEiBH8+uuvLF26lHPnzrFy5Urp9yKEE5Nd80ZKTEzk9ddflx8zYTMrV65k8ODBKKVYu3Yt11xzTbMSfk0eHh5069aNm2++ufq6/IULF7Jv375ah/PDw8MpKiri1KlTTSrfx8eHYcOGcfToUf7973/ToUMH6fcihBOTpC8MKzgwsMHr2Y1Ga83KlSsJDQ2lXbt21TsALUn4NSml6Ny5M9dddx2DBg1i//79LFiwgF27dlFRUYGHhwfdu3dn165dzSp748aN9O3bl8mTJ0u/FyGcmBzeF4Z1yoV6j+/atQs/Pz8KCgoIDg7m2LFjbNy4kU6dOrU44deklKoexCgvL49t27axdetWevToQZcuXfj444+ZN29ekw/RW/q9eHh4SL8XIZyY2yR9I8yhbcS5ycF559C2dbnWKrNzhw6cPncOgHHjxlWvf2byZI7m5TVr4J361I23VatWDB48mPz8fHbu3MkPP/zA+vXr6dmzJ8nJyZSUlFxy+lZLmS2d5/1SsVqTEX4LbFWu/BYYr1xrc8mkr5QaDYy2LIeFhTkwGiEu7vS5c9Q3xp06f77RLXytNWVlZZSXl9f6a7mVl5dTWFhIRUXFBdta/s/LyyMyMhKTycSjjz5KVlZWk+ZsT0hIkDnehXByLpn0tdYLgAWW5aioqAeNNHeykWK1VbkSa5WDBw/WSt41b6WlpbUSt5eXF97e3vXefHx88PX1xcvLi4CAgHrv9/b25oknnuCRRx7h3XffZdasWY1+bfJ+GatcidV45VqLSyZ9IYyiocluLAoKCvD29sZkMhEQEICPjw9eXl7VSbrm7VJjF1gOvTb0o3Trrbfi4eFR69JUIYRrkaQvhIMUFRUxevToi27TtWtX2rdvb6eIak9BKoRwPXLJnhAOUFZWxp///Ge6du160e1+/vlnlixZQk5OjsxuJ4RoMUn6QthZZWUl9913H0opEhISCPD1bXC8gdGjR9OrVy927NjBggUL2LFjB2VlZQ6NXwhhXJL0hbAjrTWTJk3i4MGD3H///Zw5c4Ypzz1HdHQ0p06dQmuN1prvvvuOU/n5eHh4EBERQXx8PFdffTV5eXmkp6ezfv16CgsLHf1yhBAGI+f0hbCjl156iRUrVvD3v/+dgwcPUllZyaeffsrKlSsJDg6+6GPbt2/PNddcQ2FhITt37mTx4sWEhobSq1cvu573F6Kp0tPTyczMlHkZnIC09IWwk7S0ND7//HNeffVV9u/fT2BgINOnT2fJkiV07Nix0eW0atWKmJgYEhMTad++vZz3F04tPT2d5ORkwsLCZF4GJyAtfSHs4KuvvmL69Ol8/PHHbN++nfDwcB599FGWLFnC5Zdf3qwyvb296dWrFz169ODQoUPs2LGjekS9yy+/HG9vbyu/CiGaLjMzk6SkJKZMmVK9LK19x5GkL4SNLVq0iCeeeILPP/+c7du3ExkZyX333cd3331Hv379Wly+h4cHXbp0oUuXLuTl5bFjxw62bNnC5ZdfTo8ePWjVqpUVXoUQzWOZlwGQeRmcgCR9IWzol19+Yfz48XzyySfs2LGDXr16MX78eGbPns3w4cOt/nz1nfcPCwujV69e+Pv7W/35hLgUS6teBn1yDpL0hbCRzZs3c+utt/Lee++xd+9eevfuzf33388bb7zBzTffbNPntpz3j46OZu/evaxcuRIfHx+6d+9OZGTkJUfvE8KaZNAn5yFJXwgb2LdvH6NGjSI1NZUjR47Qq1cvkpKSmDp1Knfeeafd4qh53n/v3r3s3r2bbdu2VZ/3X7RokfSqFsKNuE3SN8J0mkad8tFdp9NsqMzc3FxGjhzJs88+S0FBAd26dWPq1KncfvvtTJw4sdH1Ze3PbPv27Wnfvj2FhYXs2bOHhQsXkp2dzcCBAxs9lW5dRnq/LIzwW2CrcuW3wPHl5uTk8OWXXzJlyhSUUmitmTFjBnfeeScRERE2ia0ml0z6MrWucJQzZ84wZswYJkyYgLe3Nx06dGDGjBkMHTqU5557ztHhAdCuXTvatm3LwoULueKKK0hJScHDw6PJU+kKIZruyy+/5G9/+xtr167l8ssvZ9++fSxcuBCAZ555xubP75JJX6bWNX65Roy1qKiIsWPHcv311xMREUGrVq345JNPiIiIIC0tDQ+Ppg2LYcs62LBhA127diU1NRWTyVTdq7q5z2nE98udy5VYHVfutGnTWLduHQsWVKcoRo8ezbRp0+zS18Ylk74Q9lZWVsbYsWO5/PLLGTRoEAD/+c9/qKio4NNPP21ywrel7du3c/jwYSZOnEh4eLj0qhbCjpRSREZG1lpnz861kvSFaKHKykruvfdePDw8uOWWWygsLGTTpk3s27ePH374wakGyTlw4AC7du0iPj4ek8kkvaqFsDOtNXv37q21bu/evWitpaUvnFPNcbTdndaaJ554gpycHJ588kmOHz/OsWPHWLZsGcuXL3eqa+OPHDnCli1buOGGG5wqLiHcySuvvMKCBQsYPXo0kZGR7N27lwULFvDKK6/w/PPP2/z5JemLJrGMo52UlERycjJ33303Y8aMcXRYDjN9+nRWrlzJSy+9VD2BzuzZs1m5ciVt2rRxdHjVjh8/zrp167j66qsJCgpydDhCuK177rkHgOeee6669/4rr7xSvd7WJOmLJqk7jvYvv/zi4Igc591332Xu3LnMmjWL3bt307p1a5599lmysrKaNIGOrZ0+fZqVK1cyZMgQ2rZt6+hwhHBrXbp0qdWiV0rZpYVvIUlfNEndcbTvvvtuB0fkGF9++SUzZ87kk08+qZ5A55FHHmHp0qXNnkDHFgoKCli+fDmDBw+W6XeFEJL0RdPUHUfbHadyXbRoEZMmTeKzzz6rnkDn3nvvZf78+VaZQMdaiouLWbZsGVdccQURERFWH5RGCGE8kvRFk9Xs8T1//nwHR2NfP//8c60JdHr37s0999zD559/bpMJdJqrtLSUZcuWERkZSVRUlKPDEVZSsxOtXHUhmkOSvhCX0DYoiNMFBbXWjR49miA/PwLbtuWtt97ipptuclB0FyovL2f58uWEhYXRp08fR4cjrKRuJ1pAEr9oMkn6QlzC6YIC6juJoYqL+etTT3HHHXfYPaaGVFZWsnLlSgICAhg4cKDMpudC6naizczMlKQvmkySvhAtMGfOHLZt20b//v3p378//fr1o3Xr1g6JRWtNdnY2SimuuuoqSfgupm4n2lmzZjk4ImFEDkn6Sqk2wEfAFYAG7gd2Av8P6Ar8DvxZa33avP2zwESgAnhca/293YMWbqe0tJS8vLyLbvPcc89x+PBhtmzZwpw5c9i8eTMhISHVOwGWW7du3Ww6FK/WmnXr1lFUVMR1113nVMP+Cuuo24lWWvmiORzV0n8LWKy1/v+UUj6AP/Ac8KPW+lWl1FRgKvCMUqoPMA7oC3QCliqlemitKxwUu3BRxcXFnDhxghMnTnD8+HHOnTt3yevaT506RXl5OX369GHo0KEEBwfj6enJqVOn2L17N59++ikbN27k1KlT9OvXr9aOQHR0NK1atbJK7Fu3buX48eOMGDECT09Pq5QpnI8Mmyxayu5JXykVBPwBuBdAa10KlCqlxgBx5s1mA5nAM8AYYK7WugTYr5TaAwwBVtk1cOFStNacO3euOsGfOHGC0tJS2rdvT0hICIMHD6Zt27aXbDHPnTuX/v3707t37+qEf/r0afLz8wkODuaWW27hrrvuIjAwkKKiIg4dOsSaNWv48MMP2bZtG+Hh4bVODfTv358uXbo06bXs3r2b/fv3M3LkSHx8fFpSLUIIF6fsfZ21UmoA8AGwDegP/AY8ARzWWrepsd1prXWwUuptIFtr/YV5/b+ARVrrb+qU+xDwUH3PGRISMujDDz+0wasRRqG1prS0lPPnz1ffoGoaTF9fX0wmEz4+PvWeB7/njjsoqOca9wBfXyZPmcLvv//O/v37+f333zl+/Dg9e/asbslfdtllBAQEUFRUxKlTpzh37hze3t4EBQXRtm1btNacPHmSXbt2sXHjRjZu3MjZs2fp2rVr9a1bt25ERERgMpkuiOHcuXOcPHmSTp06OdXEPkII67j11lt/01pfaa3yHJH0rwSygeFa69VKqbeAfCC5gaSfBqyqk/QXaq3/r7HPGRUVpffs2WPV12EZ6MSaczLbokxblgtV1+lbc+x9a8VaWVnJyZMnq1vyx44dw2Qy0bFjR0JCQujQoQOtWrVqUWe3+mI9f/48W7durU7glpufn191i75Hjx60b98epRSnT5/m7NmzlJWVERAQQHBwMK1ateLcuXMcPHiQzZs3s3HjRnbs2IEuK6O0svKCOIIDAjhV55LCxsbbUkb7zFr78wrGqgN3/C0werlKKasmfUec0z8EHNJarzYvf0PV+ftjSqmOWuujSqmOwPEa20fUeHw4cMRu0QpDKCsrIy8vrzrJnz59moCAADp06EBkZCT9+/fH19fXJj92Nfn6+jJo0CAGDRpUvU5rTU5OTvUOwKJFi9i4cSOHDh2id+/e1UcF2rZti7+/PwUFBZw5cwYPDw8GDx7MyJEjCQ4OZvTo0fVfOnjunE1fkxDCddg96Wutc5VSOUqpnlrrncAIqg71bwMmAK+a/1qGeksHvlRKvU5VR77uwBp7xy2cy/nz52t1uisoKCA4OJiQkBD69OlDSEhIrcPdjhyCVilFly5d6NKlC6NHj65eX1hYyObNm9m0aRMbN27k//7v/9i0aRNt2rSpPirQvXt32rRpQ0XFxfutfvfdd4wZM0Yu0xNCXJSjeu8nA3PMPff3AfcBHsA8pdRE4CAwFkBrvVUpNY+qnYJyIKmpPffz8vJIT0+XXq9OLCcnhzffnMHGjdn07TuIlJTniYioOsCjtaaoqKi6w93x48c5f/58dae7QYMG0bZtW8P1Wm/VqhVDhw5l6NCh1esqKyv5/fff+fDDD/Hy8uK7775j48aNHDt27KJlTZs2jddee43XXnuNq6++2tahi2aQIXSFM3BI0tdabwDqO0cxooHtXwZebu7zBQQEyLCVTiwnJ4ehQ/szblw+KSkVLFmynquu+n/MnTsfLy8vTpw4QWVlJSEhIYSEhFS3fl2xVevh4cHll1/O0KFDa50fPXv2LG3atGnwcc899xwFBQXccccdxMTEMH36dHr16mWHiEVjyBC6wlm4xYh8YWFh/PnPf5ZhK53Um2/OYNy4fGbOrDqAM2pUBVoX8P77s/jf/32V6OhoAgICXDLJN9alRvkzmUwcPXqUt956i127dhEbG8sf//hH/vrXv9KxY0c7RSkaIkPoCmfhFsN25ebmkpaWRlxcnKNDEfXYvHk18fG1z9jccEMFJ04cJDIyksDAQLdO+BbBgYEouODm6+nJ5MmTad26NefPn8dkMvHJJ5/QunVrrrjiCv7yl7+Qn5/v0NjdXVxcHGlpacyYMUN+i4RDuUXSP3funAxb6cSio69iyZLa5+OXLPEkOnqIgyJyTqfy89FaX3ArLi9n9uzZvP/++8yYMYPw8HDy8vLo0qULX3zxBTk5OfTo0YN3332X0tJSR78Mt5SYmMisWbPIzc2V3yLhUG5xeL99+/byJXNikyZNYejQOUA+8fEVLFniydy5QWRnT3F0aIYRGxvLqlWr+Pbbb3n22Wfp2rUrkyZN4uDBgwwbNoyxY8eSlpZGWloar7zyCn/+85+d/uhJRkYGq1atcpmObzKErnAGbtHSF84tIiKC7OyNwCOkpl5JefkDZGdvrO69LxpHKcWf/vQntm7dym233cbEiRPJysrisssuY//+/YwdO5bXXnuN1NRUhgwZwrJlyxwdcoMyMjJ49dVX6dChA8nJyaSnpzs6JCFcgiR94RQiIiKYOXMWCxas4NVX35CE3wLe3t488sgj7N69m+7du3PXXXexb98+QkNDycnJYfLkyTz22GM88MAD3HzzzWzevNnRIVfTWnPkyBF+/fVXxowZw7XXXktSUhKZmZmODk20UHp6OpMnT5YdOAeTpC+EiwoICGDatGls27aNkpIS7rnnHgoLC/H39+f06dO89tpr3HjjjYwcOZJ7772XgwcPOizW8vJydu/ezX/+8x82btxI7969+eyzz/jxxx+ZP3++dHwzOMsli2FhYXLkxsEk6Qvh4kJDQ0lLS2PZsmX8+uuvPPbYY/j5+VFeXo7Wmo8++oiIiAgGDBjAlClTOH36tN1iKyoqYsOGDaSnp3P06FGGDBnCqFGjuP3225k+fTqnT5/m9ttvrzWssTCempcsypEbx3KLjnxCCOjevTtz5sxh06ZNTJkyhZMnT/KXv/yFM2fO0KFDB+bMmcN3331Hjx49eOaZZ3jsscfw9fW1SSx5eXns3LmT3NxcunXrxg033EBAQECtbRISEhg7dix5eXmsWLGCa6+9lnbt2tkkHmFbcXFx1YMSpaWlMWvWLAdH5L7cJulbe+x1y9Sszl6mLcu1sGbdGqkOjFSv8N94+/Xrx6JFi1i4cCEvvvgiISEhPP300xw6dIgBAwYwZswYPvzwQ/75z3/y4osvMm7cuAaHOG5KHVRWVnL48GH27NnD+fPniYqKIjo6unqOhJqvt2a5rVq1on///ixbtoy4uDhatWrVnJdfzQi/BbYq11Gf2fj4eFJTU8nKyiI1NZX4+PhLvg9GqldblmttLpn0lVKjgeqZTcLCwhwYjRDORynFLbfcwo033sgXX3zBQw89xPDhw7n//vvZt28ft9xyC/feey8zZszgn//8Jy+99BLx8fHNusyvtLSU/fv3s3fvXlq1akXPnj3p2LFjk8rq1KkTRUVF/Pzzz8TFxeHj49PkOIRjJSQkkJCQ4Ogw3J5LJn2t9QJggWU5KirqQVtNqWqLco0Uq63KlVjtV+4jjzzChAkTeOutt7j77rsZN24co0ePZteuXTzyyCMopXjmmWf45z//yYwZM+o9v15frPn5+ezcuZMDBw7QuXNnrrvuOtq2bdvsWKOjoykrK2Pt2rVcd911zZ5gyejvl7OWaatyjRSrLcu1FunIJ4TA39+fZ599lp07d+Lj48Odd97J6dOnad26NSdOnOCll17ij3/8IwkJCdxxxx3s27ev3nK01uTm5pKZmcnSpUsxmUzccsstDBs2rMkJvz4DBw7EZDKRnZ2N1rrF5QnhbiTpCyGqtW/fnjfeeIPffvuNXbt28eCDD+Lp6YlSirKyMt5//30WfPstkZGRKKXw9/fH398fpRRtAgJYuHAh69atIyIigjFjxtCvXz+rtnyUUlx99dUUFhayadMmq5UrhLuQpC+EuEC3bt2YM2cOCxcuJD09neeff5727dtTUFBAYWkpGi64nS0spLi4mIEDBxIZGdnsw++X4unpybXXXsvBgwfZs2ePTZ5DCFd1yXP6SilfIAGIBToBxcAW4D9a6622DU8I4UgxMTEsWbKEH374gWeeeeaSl/C9//77PPzww3h7exMTE8PAgQOJiYkhJiaGyy67zGrj/ZtMJuLi4li6dCn+/v506tTJKuUKYQ/p6elkZmY6ZF6JiyZ9pdRfqeoFnwmsBo4DvkAP4FXzDsFTWms5ziaEC7vhhhsYOXIkX375JdnZ2Q1ud9NNN3HnnXcSGBhIUVERe/fu5ZNPPiE5Obn6KEDNnYHu3bs3+4hAYGAgsbGxrFixgri4OKv0GRDC1iyjEyYlJVWPXWDPxH+plv5arfVfG7jvdaVUB6CLdUMSQjgjDw8P7r77bu65554Gt2nXrh1aaw4fPkxRUREmk4nrrruOP/7xj7Rp04bz58+zd+9e/v3vf/OXv/yFY8eO0b9//+qjAQMHDqRPnz6Njql9+/YMHjyYFStWEB8f3+Jr+IXtObKV6wxqjk5oWXaapK+1/s8l7j9OVetfCCFIS0tjy5YtdOnShZiYGPr160fr1q3RWnPw4EGKi4vx9PRk2LBh3HTTTbRr147S0lIOHDjATz/9RGpqKvv376dnz54MGDCAIUOGMHDgQPr164e/v3+9zxkREUFhYSGZmZnEx8fLNfxOzNGtXGfg6NEJG3WdvlLqSuB54DLzYxSgtdb9bBibEMIJBQcGogoK6l2/Zs0aysrK2LFjB+vWrWPdunUsXLiQ9evX0759ewYNGkR0dDTt2rWjsrKS/fv3U1paSkVFBQMHDiQuLo6QkBDKy8s5dOgQv/32Gx999BHbtm0jMjKy1umBgQMH0rp1awB69epFYWEhWVlZxMXF2awToWgZR7dynYHl9WZmZjJr1iznOqdfwxzgaWAzUGm7cIQQzu5Ufn71/5ahVGteluft7U10dDTR0dFMmDABqBqCd8+ePaxfv55169bx9ddfs379ery9vRkyZAj9+vUjKCiI8vJytm/fTllZGV5eXvTq1YshQ4bQrl07PD09yc3NZcuWLXz99dds3LiRjh071toR0FqzZs0ahg4dypkzZ2jTpo3VOg+KlnN0K9dZJCYmOmxnp7FJ/4TWWuZCFEI0i4eHBz169KBHjx7cfvvtQNVAPjk5OdU7AosXLyYzMxNfX19uvfVWQkJCqofrPXr0KEVFRVRUVBAREUGvXr147LHHMJlMnDlzhu3btzNjxgw2bdpE4dmznC8vvyCG4MBAPvn8c3u/dFGDvVu57t5/oD6NTfrTlFIfAT8CJZaVWutvbRKVEMLlKaXo0qULXbp0YcyYMUDVj/Sjjz5KZWUlGzduZN++ffz666/k5uYyYMAAYmJiiIqKolWrVhQWFnL06FFKSkpo164diYmJjB8/nvvuu4/6xuqr75SEsD97tXLt3X8gIyODrKwsRo4c6dQ7GI1N+vcBvQBv/nt4XwOS9IUQVmP5sVy6dCnjx4+ndevWPPvss1x22WVs2LCB9evXs3btWtatW8e+ffvo06cPMTEx9O3bl7Zt28rQvKKaPfsPpKenk5KSYogOio1N+v211tE2jcSO5JCPEM4rMTGR+Ph4oOq0wJIlS/Dx8SEuLo64uLjq7YqKiti0aRPr1q2rPkWwfft2B0UtnI09+w8YqYNiY5N+tlKqj9Z6m02jsSFLh6OMjIxae2QlJSXNmu7RSHNdG2nedyPVgZHqFYxVB5ZyfX19GTZsGMuXL0drTbdu3aq3UUrRv39/+vfvz3333QdUTePbpk2bBss9dOiQIerVVuUa6TPb0ljj4+NJTU0lKyuLcePGsXTpUkpKShg5cqSVIvyvYcOGkZKSAlTtYKSmplr9c2YtjR17/xpgg1Jqp1Jqk1Jqs1LKaUfhU0qNVkp9YLkVFhZW35eVlVW9R5aUlERWVpYDIxVCXIq/vz+xsbFs27aNQ4cOXXRbL6+Lt2P+9re/ceTIEWuGJ5xYQkICsbGxzJ07l/DwcFJSUli0aJFNnmf69On8/vvvpKamNqshaS+NbemPsmkUVqa1XgAssCxHRUU9aLmkaOTIkRcc8mnJLGAyf7TUgZFitVW5to7Vz8+PkSNH8tNPP9GqVasGx9rftWsXPkqh6jm3H2AyMWnSJP74xz+yYsWK6mv8rR2rtRnx/XKmMletWlXr0Ht2dja33Xab1WO97bbbbFKutV1q7P0ArfU5rfWBS21j/dBsw9EDIwghmqdNmzb84Q9/YMWKFcTGxhISEnLBNuvWrSPhttsYOXIkwcHBfPTRR3Tr1o2cnBzuv/9+zp07x7XXXsutt97K4sWLMZlMDnglwp7qnttPTU11cESOdanD+/OVUjOVUn9QSlUPaq2UulwpNVEp9T0GOwoAVYn/9ddfl4QvhMG0b9+eq6++mqysLE6fPn3B/evWrWPAgAEUFxeze/duBg0axLvvvouPjw8ZGRl4eXlx1VVXERISwj333ENlpYw15uoSExOZNWsWubm5zJo1y6kPvdvDRZO+1noEVdfm/w+wVSl1Vil1EvgCCAMmaK2/sX2YQghRJSwsjMGDB7N8+XLya4wOCFVJv3v37vj5+bFhwwZiYmLw8vJi7ty57N69m61bt1JUVMTtt9/OiRMnePLJJ+UyPzcgDb3/umRHPq31Qq31XVrrrlrr1lrrdlrrq7XWL2utc+0RpBBC1BQREUF0dDTLli2jqKgIqBrhb926dQQFBREcHMy6deuIiYkBqjoDZmRkMGfOHM6dO8exY8d46qmnyMzM5LXXXnPkSxHCrhrbe18IIZxKZGQkPXv25KeffuL8+fMcOHAAPz8/CgsLCQoK4uTJk0RGRlZv365dO6ZNm0ZqaioBAQHs37+fGTNm8M477/DZZ5858JUIYT+N7b0vhBBOp1evXpSWlrJs2TIKCwuJiYnh1KlTaK0ZOHAgHh612zUdOnRg0aJFjBgxgg8//JAdO3bw/vvvM378eDp06MCoUYbroiREk1yq9/5C4FGt9e/2CUcIIZomOjqaPwwbxlnzYf6MjIzq+9oGBdWaFRCgX79+fP3114wdO5bPPvuMHTt28K9//Yt77rmHhQsXMnjwYLvGL4Q9Xerw/qfAD0qp55VS3naIRwghmkQpxdmiIjRccDvdwCQ7cXFxvPPOOzzwwAP07t2bffv28fbbb5OYmMju3bvtF7wQdnap3vvzgIFAEPCrUipFKTXZcrNLhEIIYQNjx45l6tSpPP7443Tv3p0TJ04wbdo0brzxRnJzpY+ycE2N6chXBhQCJiCwzk0IIZxaTk5Og/clJydz22238b//+7+Eh4ejlGLChAnccsstFMhUvMIFXeqc/ijgdSAdiNFaF9klKiGEsJJ///vf+Pr60rdvX8rLyy+4/5VXXuHee+/lo48+4u6778bHx4chQ4bwpz/9iYyMDHx8fOwWa05ODm++OYPNm1fTu3cMjz32FN27d7fb8wvXd6mW/vPAWK31VEn4QggjevTRR7nyyivZs2cPO3fu5JtvvuH333+vvl8pxUcffURFRQVLly7Fw8ODP/zhD7Rq1Yr777/fbqP25eTkMHRof+BdJk9ei5fXR8TFDbnokQohmupS5/RjtdZb7RWMEEI0R3BgIAouuAUHBuLl5UVMTAwTJkwgMjKSVq1asWDBAt5//31WrFhBUVER3t7efP3112zatIk9e/ZQUFDA+PHj+f3333nmmWdsHr/Wmtdee5lx4/KZObOCUaNg5swKxo0r4M03Z9j8+YX7cJvr9I0wh7aR5uWuyZnm0LZnuUaqVzBWHTS13MPHjjV4X8169PX1JS4ujj/84Q9s376drVu38uuvv9KxY0eio6OZN28eI0eO5IknnuDIkSM899xzPP3007Rv357HH3+8xbFqrSkqKiI/P5+CgoJaf3/5ZSl//3tFre3j4yuYNm0Jv/32G+3ataNdu3a1pg/u3KEDp89dON9ZcEAAh48fb1GszeGuvwW2LNfaXDLpK6VGA6Mty2FhYQ6MRgjhbDw8POjbty99+/blzJkzrF+/nmXLluHp6ckrr7zCM888wyuvvMLtf/wj50pKmDp1KlOnTq1+fENJFaoSe0lJCadOnSIvL48zZ85w9uxZzp07R1FREZWVlXh4eKC1pqKigrKyMkpKSjCZ2vD994pRo/47F8DixZCfX8Y//vEPTCYTvr6+5Ofnc+LECY4dO8bpc+eob+YAVc+OgBDgoklfa70AWGBZjoqKetDe80enp6eTmZlJXFxckyd5MNJc17Yq15GxNvW9M1K92qpcI8fq5+dHx44dqaysZMeOHaxfv56HH36YxYsXc66kpMGkOmLECIqKitBa4+XlhY+PDyaTCX9/f1q3bk15eTnnz5+npKSE8vLy6r4Blm38/Pyqb/7+/gwefA2ffbYVpUq54YZKfvjBkzlzTEyb9iRhYWHVowuWl5dTUlLC+fPneeihhxr9Oht7X0u4+2fLluVai0smfUdLT08nOTmZpKSk6nmcZXYnY5D3zn15eHjQp08f+vTpQ35+PvPmzbvomPx33HEH5eXleHl5YTKZqm8+Pj54enpSUVFBRUUF5eXllJeXVy9XVFRQWVl5wd/KykoefTSFX37JZOnSg4SGXsakSTfi7e3NmTNn8PT0rL75+Pjg7+9vx9oRrkKSvg1kZmaSlJTElClTqpclcRiDvHcCICgoiAceeIAHH3ywwW1MJhNBQUF4enri5eVVfbMkZm9vb7TW1QndkuDLy8spKyujrKyM0tJSKisrKSsr4/z585SXlxMdPYji4j4UFRWxfv16iouLKS4upqioqPp/y7IQTSVJ3wbi4uKqW4lpaWnMmjXLwRGJxpL3TjTW7NmzL0jElpuHh0etw/eWQ/gXW/bz8yMoKIjQ0FD8/Pzw9PSsPlXQ0ONCQkIcXQ3CYCTp24ClZZiZmcmsWbOkpWgg8t6JxpoxYwZ+fn7VCT44OLg6GdfsYd9clp7wFztHHBwYiKpn5MDgQBkwVdRPkr6NJCYmSsIwKHnvhMXFkurw4cOBxiVnW6k5g+DGjRsB6N+/v93jEMbRmLH3hRDCLZ3Kz0drfcGt7nS9zqBjx44cPXrU0WEIJydJXwghXED79u05d+6cYQaJEY4hSV8IIVyAh4cHHTp0kGmBxUVJ0hdCCBfRsWNHSfrioiTpCyGEiwgLCyM3Nxet6xtHUAgHJn2llKdSar1SKsO83FYptUQptdv8N7jGts8qpfYopXYqpW50VMxCCOHMAgMD8fDw4OzZs44ORTgpR7b0nwC211ieCvyote4O/GheRinVBxgH9AVGAe8opTztHKsQQhiC9OIXF+OQpK+UCgduAT6qsXoMMNv8/2zg1hrr52qtS7TW+4E9wBA7hSqEEIYiSV9cjKNa+m8CU4DKGutCtdZHAcx/O5jXdwZyamx3yLxOCCFEHaGhoZw8eZLy8nJHhyKckN1H5FNKJQDHtda/KaXiGvOQetZd0EtFKfUQUO88kyEhIcyfP78pYYomkLq1DalX23CHes3Ly+Prr7+2+0x87lC3RueIYXiHA4lKqZsBXyBIKfUFcEwp1VFrfVQp1RE4bt7+EBBR4/HhwJG6hWqtPwA+qO8Jo6Ki9JgxY6z5Gmwy9KathvNsTrmNnVN+/vz5WLNunakOHFGmhbXrFYxVB7Yq1yj12tJyt2zZQmlpKTExMVYr81Lc+bfAluVam90P72utn9Vah2utu1LVQe8nrfXdQDowwbzZBMCyy5gOjFNKmZRS3YDuwBo7h+1WLHPKh4WFkZycTHp6uqNDEkI0gZzXFw1xpgl3XgXmKaUmAgeBsQBa661KqXnANqAcSNJaVzguTNcnc8oLYWxt27bl/PnzFBUV2f0Qv3BuDk36WutMINP8/0lgRAPbvQy8bLfA3JzMKS+EsSmlCA0NJTc3l8svv9zR4Qgn4kwtfXEJjT3P3lIXm1PeXjEIIVrGcohfkr6oSZK+QVjOsyclJVW3wm2d+OuWX18MQgjn1LFjRzZs2IDWGqXquwhKuCNJ+gbR1PPsNVvk8fHxNovh2muvtUrZQgjr8vf3x9fXl1OnTtGuXTtHhyOchEy4YxBxcXGkpaUxY8YM0tLSiIuLa3Dbur3vMzIy7B6DEMLxZNY9x0tPT2fy5MlOcxWUJH2DSExMZNasWeTm5l5wnr2umi3ypKQksrKy7B6DEMLxwsLC3OLSPWdLrBbOePmzHN43kPrOs9enbu/71NRUu8cgrEs6UIrm6NChAz///DNlZWV4e3s7OhybsHd/p6ZwxsufJem7oLq97611Tl84hjP/qAnn5uXlRbt27Th27Bjh4eGODscmnDGxWjjj5c+S9F1UzRa5ZXhIYUzO/KMmnJ/l0j1XTfrOmFgtLnb5s6NI0hfCyTnzj5pwfhUVFUyf/jfOnTtO794xPPbYU3Tv3t3RYVmNMybWmpztlKjbJH1rt3bPnz9v1fJsVaYty7WwZt0aqQ7sVa/x8fGkpqaSlZVFamoq8fHxzapzI9WBLevWCL8F1ir30KFD3HhjLLffns8NN1SyZMk6rr32/7F8+RqbtPwd9VsQHx9ffRrzUjE48/tlDy6Z9JVSo4HRluWwsDAHRiNEyyUkJJCQkODoMITBvP32TMaNK2DmzEoARo2qAAp4++2ZvPrqG44NTjiESyZ9rfUCYIFlOSoq6kFbTXdoi3KNFKutypVYjVWuxOqc5W7fvo7Jk2vPTxYfX8Hrr683TN06Y706olxrkev0bcxZrx91BlI3QthWdPRVLFniWWvdkiWeREcPcVBEwtEk6duQMw7M4CykboSwvUmTpjB3bhBPPeXJ4sXw1FOezJ0byKRJUxwdmnAQSfo2VHdkvMzMTEeH5DSkboSwvYiICLKzNwKP8PrrQygvf4DMzDVEREQ4OjThIC55Tt9ZyKVWDZO6sR4ZrU9cTEREBDNnVn2/ZMwOIUnfhpz9+lFHkrqxDhmtTwjRFJL0bcxZBmbIyMhg1apVTtUadJa6MTIZrU8I0RRyTt8NZGRkkJKSIp3mXJA1pjuWqyiEcB/S0ncDWVlZ0hp0US09TSKnB4RwL5L03UBsbCwpKSmAdJpzRS05TSKnB4RwL5L03YBl+NZVq1ZJpzlRS82rKN544w3uv/9+B0ckhLAlSfoGYI1LshISEhg7dqyVIxNGl5iYyOrVq3n11VdJTEzkiy++4KqrrpIdQyFclNt25DNK5yUZuU7YWnFxMVOnTuXTTz+VgZKEcHFu2dI3UuclOecqbE0GShLCfbhN0q85EtXSpUtrJdKlS5dWz8XcWPaam3zYsGG1OuGlpqY2eVQte837bg1GmuvaSPUKDccbHx9PamoqWVlZpKamEh8f3+jnNtL7ZWGvenXGco30mTVSvdqyXGtzyaSvlBoNjLYsh4WF1bq/bm/21NRUu8bXFJZOeJYfZJlTXdhCQkKCfLYMIiMjg6ysLGJjY+U9E03mkklfa70AWGBZjoqKerDmHMdjx47FZDJZZQhYe8wfPXbsWKt0wjPS/NESq7HKlVjtU256ejopKSkkJSWRkpKCyWRq1u+XkerASLHaslxrccmk3xgyBKwQwmikj49oKbdN+kIIYTTS6VK0lCR9IYQwCJmdUrSUJH0hhDAQOTUpWsJtB+cRQggh3I0kfSGEEMJNSNIXQggh3IQkfSGEEMJNSNIXQggh3IQkfSGEYRlltkwhnIUkfSGEIcm000I0nVynL4QwJBmSVoimk6TvIOnp6WRmZhIXF+dWP1Tu+rqF9cmQtEI0ndskfWeaQzsjI6N6pqzk5GRKSkpISEgw7DzPja3bhl53TUaqA2ep18YyUh00ptz4+HhSU1Orp52Oj49vVJ0ZoV5tVa6RPrNGqldblmttLpn0lVKjgdGW5bCwMAdGc6GsrKxahyWzsrLcYl5sd33djuAuc64nJCS49OsTwtpcMulrrRcACyzLUVFRDzrT3MkjR4684LBkzXKcKVZrlnup192cMpvKHeb7vtSc6+5QB/Yu02jlGj1Wa5wmNFIdWJNLJn1n564zZbnr67Y36eAmmspIfW0sV21YThMCTh+zM5Gk7yDuOlOWu75ue5IObqIpjJZEZae2ZSTpC+Fi5IiKaAqjJVHZqW0ZSfpCuCA5oiIay2hJVHZqW0aSvhBC1MNI57lbwohJVHZqm0+SvhBC1GG089wtJUnUfUjSF0KIOox2nluIxpKkL4QQdRjtPLcQjSVJXwgh6jDieW4hGkOSvhBC1EPOcwtX5OHoAIQQQghhH3ZP+kqpCKXUMqXUdqXUVqXUE+b1bZVSS5RSu81/g2s85lml1B6l1E6l1I32jlkIIYRwBY5o6ZcDT2mtewNDgSSlVB9gKvCj1ro78KN5GfN944C+wCjgHaWUpwPiFkIIIQzN7klfa31Ua73O/H8BsB3oDIwBZps3mw3cav5/DDBXa12itd4P7AGG2DVoIVogPT2dyZMnk56e7uhQhBBuzqEd+ZRSXYGBwGogVGt9FKp2DJRSHcybdQayazzskHld3bIeAh6q73lCQkKYP3++FSMXNUndNmzNmjXMnj2bxx9/nEcffZTVq1czZEjj9lmlXm1D6tV2pG6dn8OSvlIqAPg/YJLWOl8p1eCm9azTF6zQ+gPgg/oKiIqK0mPGjGluqPUqLi4GrDt3si3KtGW5UPUlt2bdGqkOGlPm8uXLefzxx6sHecnNzW1UfVm7XkE+s2CcerVVufJbYLxyrc0hSV8p5U1Vwp+jtf7WvPqYUqqjuZXfEThuXn8IiKjx8HDgiP2iFaL5ZJAXIYzF1edccETvfQX8C9iutX69xl3pwATz/xOA+TXWj1NKmZRS3YDuwBp7xStESyQmJjJr1ixyc3NlkBchnJxlzoWwsDCSk5Ndsh+OI1r6w4F7gM1KqQ3mdc8BrwLzlFITgYPAWACt9Val1DxgG1U9/5O01hV2j1qIZnLnQV5cvdUkXIs7zLlg96SvtV5J/efpAUY08JiXgZdtFpQQwurcbaY6YXzucDpOhuEVAmmR2oI7tJqEa3GHORck6Qu3Jy1S23CHVlNNsuPoGlz9dJwkfeH2pEVqG+7QarLIyMggJSVFdhyF05OkL9yeu7VI7cnVW00WWVlZsuMoDEGSvnB77tQiFbYRGxtLSkoKIDuOwrlJ0hcC92mRCttISEjAZDLJjqNwepL0hRDCCmTHURiBI6bWFUIIIYQDSNIXQggh3ITbHN63zIBkLefPn7dqebYq05blWlizbo1UB0aqVzBWHdiybo1Qr7Yq10ifWSPVqy3LtTaXTPpKqdHAaMtyWFiYA6MRQgjbycjIICsri9jYWBISEhwdjnByLpn0tdYLgAWW5aioqAdtNcexLco1Uqy2KldiNVa5Eqtjyk1PT68eFCglJQWTydSozoRGqgMjxWrLcq3FJZO+EEK4AxlN0hicaYhm6cgnhBAGFRcXR1paGjNmzCAtLY24uDhHhyTqsMztERYWRnJyMunp6Q6NR1r6QghhUDKapPNztqMxkvSFEMLAZFAg5+Zsc3tI0hdCCCFsxNmOxkjSF0IIIWzImY7GSEc+IYQQwk1I0hdCCCHchCR9IYQQwk1I0hdCCCHchCR9IYQQwk1I0hdCCCHchCR9IYQQwk24zXX6RphD26jzPLvrHNpGqlcwVh3Ysm6NUK+2KtdIn1kj1asty7U2l0z6SqnRwGjLclhYmAOjEUIIIZyDSyZ9rfUCYIFlOSoq6kEjzZ1spFhtVa7EaqxyJVZjlSuxGq9ca5Fz+kII4WbS09OZPHmyw6d5FfYnSV8IIdxIRkaGU83vLuzLJQ/v21p6ejpLly4lNjaWsWPHOjocIYRotKysLKea313Yl7T0myg9PZ3k5GTCw8NJSUmRvWQhhKHExsaSlpbGjBkzSEtLIy4uztEhCTuSln4TZWZmyl6yEMKwEhISMJlMTjO/u7AvSfpNFBcXR3JyMgBpaWnMmjXLwREJIUTTONP87sK+JOk3keWLsnTpUlJTU+WLI4QQwjAk6TdDYmIi8fHxjg5DCCGEaBLpyCeEEEK4CUn6wqpk0A/rkzoVQliLJH1hNZbLGWXQD+uROhVCWJOc0xdWI5czWp/UqRDCmtwm6RthOk2jTvloqdthw4aRkpICVF3OmJqa2uR6N1Id2KNerVGnFkaqA5la11jvl4VMrev8XDLpy9S6jpGQkABUDfOZmppavSyaT+pUCGFNLpn03XFq3fT0dDIzM4mLi2vw8K896mDs2LFWmY/A1d+vppRrrTqtW6412ev9asznvKllWouRypVYjVeutUhHPhcgnb2ErTjTlQPyORdG4Uzfm7ok6buAmp29kpKSyMzMdHRIwgU4W5KVz7kwAmf73tTlkof33Y3MByBswdmuHJDPuTACZ/ve1CVJ3wVYPlAya5awJmdLsvI5F0bgbN+buiTpuwiZNUtYmzMmWfmcC2fnjN+bmiTpCyEaJElWiKZz5u+NdOQTQggh3IQkfSGEEMJNSNIXQggh3IQkfSGEEMJNSNIXQggh3IRhkr5SapRSaqdSao9Saqqj4xFCCCGMxhBJXynlCaQBNwF9gDuUUn0cG5UQQghhLIZI+sAQYI/Wep/WuhSYC4xxcExCCCGEoRhlcJ7OQE6N5UPAVTU3UEo9BDzUwONLlFJbbBBXa+CsAcq0ZbntgTwrl2mkOjBSvYKx6sAW5RqpXm1VrpE+s0aqV1uV29OqpWmtnf4GjAU+qrF8DzCrCY//1UZxfWCEMm1crtXr1kh1YKR6NWAd2CJWw9SrAd8vt/4tsOH7ZdV6Ncrh/UNARI3lcOCIg2IRQgghDMkoSX8t0F0p1U0p5QOMA5xrkmIhhBDCyRninL7Wulwp9RjwPeAJfKy13urgsAAWGKRMW5ZrC0aqAyPVKxirDoxUt0aqA6lX45VrNcp8zsClKaV+1Vpf6eg4XJHUrW1IvdqG1KvtSN3ahrXr1SiH94UQQgjRQu6S9D9wdAAuTOrWNqRebUPq1Xakbm3DqvXqFof3hRBCCOE+LX0hhBDC7UnSF0IIIdyEyyd9mZ2v+ZRSEUqpZUqp7UqprUqpJ8zr2yqlliildpv/Btd4zLPmut6plLrRcdE7P6WUp1JqvVIqw7ws9WoFSqk2SqlvlFI7zJ/dYVK3LaeUetL8O7BFKfWVUspX6rV5lFIfK6WO1xwevjl1qZQapJTabL7vn0opdanndumkL7PztVg58JTWujcwFEgy199U4EetdXfgR/My5vvGAX2BUcA75vdA1O8JYHuNZalX63gLWKy17gX0p6qOpW5bQCnVGXgcuFJrfQVV46WMQ+q1uT6lql5qak5dvkvVnDPdzbe6ZV7ApZM+Mjtfi2itj2qt15n/L6Dqx7MzVXU427zZbOBW8/9jgLla6xKt9X5gD1XvgahDKRUO3AJ8VGO11GsLKaWCgD8A/wLQWpdqrc8gdWsNXoCfUsoL8KdqKHSp12bQWq8ATtVZ3aS6VEp1BIK01qt0VY/8z2o8pkGunvTrm52vs4NiMTSlVFdgILAaCNVaH4WqHQOgg3kzqe/GexOYAlTWWCf12nKXAyeAT8ynTj5SSrVC6rZFtNaHgVTgIHAUOKu1/gGpV2tqal12Nv9fd/1FuXrSr+/8hlyj2ERKqQDg/4BJWuv8i21azzqp7zqUUgnAca31b419SD3rpF7r5wXEAO9qrQcChZgPkzZA6rYRzOeXxwDdgE5AK6XU3Rd7SD3rpF6bp6G6bFYdu3rSl9n5Wkgp5U1Vwp+jtf7WvPqY+dAS5r/HzeulvhtnOJColPqdqlNO1yulvkDq1RoOAYe01qvNy99QtRMgddsyI4H9WusTWusy4FvgaqRerampdXnI/H/d9Rfl6klfZudrAXNP0H8B27XWr9e4Kx2YYP5/AjC/xvpxSimTUqobVR1L1tgrXqPQWj+rtQ7XWnel6jP5k9b6bqReW0xrnQvkKKV6mleNALYhddtSB4GhSil/8+/CCKr6+Ei9Wk+T6tJ8CqBAKTXU/J6Mr/GYhmmtXfoG3AzsAvYCzzs6HiPdgGuoOly0Cdhgvt0MtKOqd+lu89+2NR7zvLmudwI3Ofo1OPsNiAMyzP9LvVqnTgcAv5o/t98BwVK3VqnXvwE7gC3A54BJ6rXZdfkVVX0jyqhqsU9sTl0CV5rfj73A25hH2b3YTYbhFUIIIdyEqx/eF0IIIYSZJH0hhBDCTUjSF0IIIdyEJH0hhBDCTUjSF0IIIdyEJH0hhBDCTUjSF0JUU1XTKe9XSrU1Lwebly+rZ1s/pdTypsyeppR6TCl1nzVjFkI0nlynL4SoRSk1BYjSWj+klHof+F1rPb2e7ZIAL631W00o2x/4WVeNiy+EsDNp6Qsh6nqDqiFXJ1E1KuPMBra7C/Own0qpOHOrf55SapdS6lWl1F1KqTVKqc1KqUgArXUR8LtSSqZZFcIBJOkLIWrRVROqPE1V8p+ktS6tu415LovLtda/11jdH3gCiAbuAXporYcAHwHJNbb7FYi1TfRCiIuRpC+EqM9NVI0NfkUD97cHztRZt1ZrfVRrXULVWOA/mNdvBrrW2O44VdOzCiHsTJK+EKIWpdQAIB4YCjxpme6zjmLAt866khr/V9ZYrqRqnnsLX/PjhRB2JklfCFHNPEXnu1Qd1j8IvAak1t1Oa30a8FRK1U38jdGDqpnBhBB2JklfCFHTg8BBrfUS8/I7QC+l1LX1bPsDVR39mmo4sLSZ8QkhWkAu2RNCNItSaiAwWWt9jy0fI4SwHmnpCyGaRWu9HljWlMF5qOoA+BcbhSSEuARp6QshhBBuQlr6QgghhJuQpC+EEEK4CUn6QgghhJuQpC+EEEK4CUn6QgghhJv4/wEcz/T1+dMteAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ndata = 100; seed= 13013\n",
    "np.random.seed(seed=seed)\n",
    "index = np.arange(0,ndata,1); loc = np.random.rand(ndata,2)*1000; lab = np.zeros(ndata)\n",
    "visited = np.zeros(ndata)\n",
    "KDtree = spatial.cKDTree(loc)\n",
    "plt.scatter(loc[:,0],loc[:,1],marker = 'o',s=10,c='white',edgecolor='black',cmap=plt.cm.inferno)\n",
    "\n",
    "i = 7; radius = 100; minpts = 4\n",
    "\n",
    "lab = np.zeros(ndata)\n",
    "\n",
    "step = 0; maxstep = 20\n",
    "\n",
    "def search(i,loc,lab,radius,minpts,direct,iprev,step,maxstep):\n",
    "    if visited[i] == 0 and step < maxstep:\n",
    "        print('at point:' + str(i) + ', step: ' + str(step))\n",
    "        visited[i] = 1\n",
    "        l = KDtree.query_ball_point(loc[i], radius)\n",
    "        #print('at: ' + str(i) + ', found ' + str(len(l)) + ' points')\n",
    "        l.remove(i)\n",
    "        print(l)\n",
    "        step = step + 1\n",
    "        if len(l) + 1 >= minpts:\n",
    "            lab[i] = 2\n",
    "            plt.scatter(loc[i,0],loc[i,1],s=40,marker='s',c='red',edgecolor='black',cmap=plt.cm.inferno,zorder=100)\n",
    "            for il in l:\n",
    "                plt.plot([loc[i,0],loc[il,0]],[loc[i,1],loc[il,1]],color='grey',alpha=0.7,lw=1,zorder=15)\n",
    "            \n",
    "            if direct == 0:\n",
    "                plt.plot([loc[i,0],loc[iprev,0]],[loc[i,1],loc[iprev,1]],color='black',lw=4.0,zorder=10)\n",
    "                plt.plot([loc[i,0],loc[iprev,0]],[loc[i,1],loc[iprev,1]],color='white',lw=2.0,zorder=13)\n",
    "            #print('anchor: ' + str(i))\n",
    "            #print(l)\n",
    "            for ii in l:\n",
    "                if step < maxstep:\n",
    "                    step = search(ii,loc,lab,radius,minpts,0,i,step,maxstep)\n",
    "        \n",
    "        else:\n",
    "            if direct == 0:\n",
    "                lab[i] = 1\n",
    "                plt.scatter(loc[i,0],loc[i,1],s=30,c='yellow',edgecolor='black',cmap=plt.cm.inferno,zorder=100)\n",
    "            elif direct == 1:\n",
    "                plt.scatter(loc[i,0],loc[i,1],marker='x',s=30,c='black',cmap=plt.cm.inferno)\n",
    "        print('step: ' + str(step))\n",
    "    return step\n",
    "            \n",
    "while np.sum(visited == 0) > 0 and step < maxstep:\n",
    "    #print('unvisited remaining: ' + str(np.sum(visited == 0)))\n",
    "    i = np.random.choice(index[visited == 0],replace=False)  \n",
    "    step = search(i,loc,lab,radius,minpts,1,-1,step,maxstep)\n",
    "        \n",
    "plt.xlim([0,1000]); plt.ylim([0,1000]);plt.xlabel('X (m)');plt.ylabel('Y (m)');plt.title('Points');add_grid()\n",
    "plt.subplots_adjust(left=0.0,bottom=0.0,right=1.0,top=1.1); plt.show() # set plot size  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e02a56a",
   "metadata": {},
   "source": [
    "### Interactive BDSCAN Dashboard\n",
    "\n",
    "The following code includes:\n",
    "\n",
    "* the dashboard with widgets linked to run DBSCAN over a set number of steps for visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ce02c52c",
   "metadata": {},
   "outputs": [],
   "source": [
    "l = widgets.Text(value='                                                              Interactive DBSCAN Demo, Prof. Michael Pyrcz, The University of Texas at Austin',\n",
    "                 layout=Layout(width='890px', height='30px'))\n",
    "\n",
    "maxstep = widgets.IntSlider(min=1, max = 101, value=10, step = 1, description = '$maxstep$',orientation='horizontal', style = {'description_width': 'initial'}, continuous_update=False)\n",
    "radius = widgets.FloatSlider(min=10, max = 500, value=110, step = 10, description = r'$r$',orientation='horizontal', style = {'description_width': 'initial'}, continuous_update=False)\n",
    "minpts = widgets.IntSlider(min=2, max = 20, value=4, step = 1, description = r'$min_{pts}$',orientation='horizontal', style = {'description_width': 'initial'}, continuous_update=False)\n",
    "\n",
    "ui = widgets.HBox([maxstep,radius,minpts],)\n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def run_plot(maxstep,radius,minpts):\n",
    "    ndata = 100; seed= 13013\n",
    "    np.random.seed(seed=seed)\n",
    "    index = np.arange(0,ndata,1); loc = np.random.rand(ndata,2)*1000; lab = np.zeros(ndata)\n",
    "    visited = np.zeros(ndata)\n",
    "    KDtree = spatial.cKDTree(loc)\n",
    "    plt.scatter(loc[:,0],loc[:,1],marker = 'o',s=10,c='white',edgecolor='black',cmap=plt.cm.inferno)\n",
    "    \n",
    "    step = 0\n",
    "    def search(i,loc,lab,radius,minpts,direct,iprev,step,maxstep):\n",
    "        if visited[i] == 0 and step < maxstep:\n",
    "#             print('at point:' + str(i) + ', step: ' + str(step))\n",
    "            visited[i] = 1\n",
    "            l = KDtree.query_ball_point(loc[i], radius)\n",
    "            #print('at: ' + str(i) + ', found ' + str(len(l)) + ' points')\n",
    "            l.remove(i)\n",
    "#             print(l)\n",
    "            step = step + 1            \n",
    "            lab[i] = 2\n",
    "            if len(l) + 1 >= minpts:\n",
    "                plt.scatter(loc[i,0],loc[i,1],s=40,marker='s',c='red',edgecolor='black',cmap=plt.cm.inferno,zorder=100)\n",
    "                \n",
    "                if step > maxstep - 1:\n",
    "                    for il in l:\n",
    "                        plt.plot([loc[i,0],loc[il,0]],[loc[i,1],loc[il,1]],color='black',alpha=0.9,lw=1,zorder=15)\n",
    "                    plt.gca().add_patch(plt.Circle(( loc[i,0] , loc[i,1] ), radius,color='red',fill=False))   \n",
    "                \n",
    "                if direct == 0:\n",
    "                    plt.plot([loc[i,0],loc[iprev,0]],[loc[i,1],loc[iprev,1]],color='black',lw=4.0,zorder=10)\n",
    "                    plt.plot([loc[i,0],loc[iprev,0]],[loc[i,1],loc[iprev,1]],color='white',lw=2.0,zorder=13)\n",
    "                #print('anchor: ' + str(i))\n",
    "                #print(l)\n",
    "                for ii in l:\n",
    "                    if step < maxstep:\n",
    "                        step = search(ii,loc,lab,radius,minpts,0,i,step,maxstep)\n",
    "            \n",
    "            else:\n",
    "                if direct == 0:\n",
    "                    lab[i] = 1\n",
    "                    plt.scatter(loc[i,0],loc[i,1],s=30,c='yellow',edgecolor='black',cmap=plt.cm.inferno,zorder=100)\n",
    "                    plt.plot([loc[i,0],loc[iprev,0]],[loc[i,1],loc[iprev,1]],color='grey',alpha=0.7,lw=1,zorder=15)\n",
    "                elif direct == 1:\n",
    "                    plt.scatter(loc[i,0],loc[i,1],marker='x',s=30,c='black',cmap=plt.cm.inferno)\n",
    "#             print('step: ' + str(step))\n",
    "        return step\n",
    "                \n",
    "    while np.sum(visited == 0) > 0 and step < maxstep:\n",
    "        #print('unvisited remaining: ' + str(np.sum(visited == 0)))\n",
    "        i = np.random.choice(index[visited == 0],replace=False)  \n",
    "        step = search(i,loc,lab,radius,minpts,1,-1,step,maxstep)\n",
    "            \n",
    "    plt.xlim([0,1000]); plt.ylim([0,1000]);plt.xlabel('$X_1$');plt.ylabel('$X_2$');plt.title('Interactive Density-based Clustering with BDSCAN');add_grid()\n",
    "    plt.subplots_adjust(left=0.0,bottom=0.0,right=2.0,top=2.9); plt.show() # set plot size  \n",
    "        \n",
    "# connect the function to make the samples and plot to the widgets    \n",
    "interactive_plot = widgets.interactive_output(run_plot, {'maxstep':maxstep,'radius':radius,'minpts':minpts})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70ea3f08",
   "metadata": {},
   "source": [
    "### Interactive BDSCAN Demonstation \n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "Set the radius and minimum number of points and step through the recursive BDSCAN method\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "* **maxstep** - recursive steps, **radius** - local neighbourhood, **$min_{pts}$** - minimum number of points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "54e3b1af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f6df612647b43c6b065a83077e2511c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                                              Interactive DBSCAN De…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "abff09c9463f4c6681b0a6ed4d56a7a2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                           # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1792a3",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration density-based clustering with BDSCAN. I have many other demonstrations and even basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "#### The Author:\n",
    "\n",
    "### Michael J. Pyrcz, Professor, The University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Professor, The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)  \n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce851a64",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}