{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a296acfa",
   "metadata": {},
   "source": [
    "# Compatibility with HDSBCAN\n",
    "\n",
    "A fast multicore version of HDSBCAN is not much use if it can't reproduce the clustering we would have gotten with the classic ``hdbscan`` library. The bad news is that, because it uses a different stack for compiling things (numba instead of cython), combined with floating point precision issues, and slightly different (and vastly more efficient) approches to condensed tree construction an exactly identical clustering cannot be guaranteed. On the other hand we are executing the same underlying algorithm, so aside from potentially permuting the label numbers of cluster and modulo the rare cases where the tree get juggled by precision issues, the results will be essentially the same. And we can measure that! Let's load some libraries and get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "505ad214",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sklearn.datasets\n",
    "import sklearn.metrics\n",
    "\n",
    "import hdbscan\n",
    "import fast_hdbscan\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set(rc={\"figure.figsize\":(8, 8)})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f52275ab",
   "metadata": {},
   "source": [
    "Since the assigned label numbers from ``fast_hdbscan`` and ``hdbscan`` can be different, we need a smarter way to determine if we have produced the same clustering than just checking if the output vectors are identical. Conveniently there are metrics for comparing clusterings. The two relevant ones here are the [Adjusted Rand Score](https://scikit-learn.org/stable/modules/clustering.html#adjusted-rand-score) and the [Adjusted Mutual Information Score](https://scikit-learn.org/stable/modules/clustering.html#mutual-info-score). You can read more about their details in the scikit-learn documentation linked. For our purposes it suffices to note that they measure how similar two clusterings are, with a perfect score of 1.0 if the clusterings are identical, and a score of 0.0 if they are essentially unrelated to each other (as good as random). We can run ``fast_hdbscan`` and ``hdbscan`` on the same randomly generated dataset and compare the ARI and AMI scores. In fact let's make a function to do that a whole bunch of times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1b812e6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_hdbscan_models(n_runs=100):\n",
    "    aris = np.zeros(n_runs)\n",
    "    amis = np.zeros(n_runs)\n",
    "    for i in range(n_runs):\n",
    "        data, _ = sklearn.datasets.make_blobs(n_samples=20000, centers=25, cluster_std=0.5)\n",
    "        classic_hdbscan_labels = hdbscan.HDBSCAN(approx_min_span_tree=False).fit_predict(data.astype(np.float64))\n",
    "        fast_hdbscan_labels = fast_hdbscan.HDBSCAN().fit_predict(data.astype(np.float64))\n",
    "        aris[i] = sklearn.metrics.adjusted_rand_score(classic_hdbscan_labels, fast_hdbscan_labels)\n",
    "        amis[i] = sklearn.metrics.adjusted_mutual_info_score(classic_hdbscan_labels, fast_hdbscan_labels)\n",
    "    results = pd.DataFrame({\"Adjusted Rand Index\": aris, \"Adjusted Mutual Information\": amis})\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bfdedd1",
   "metadata": {},
   "source": [
    "Now we can simply plot the distributions of those scores over the 100 different randomly generated datasets. The first thing to note is that the vast majority of the time the scores are 1.0 for a perfect match, or so close to 1.0 that it hardly matters. There are a few outlying cases where the scores are notably lower however. What happened there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "266fca8f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7f77fb4dfb20>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAMQCAYAAAC3+YP9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABUKklEQVR4nO3deXxU9aH///csmSwkIQkk7GEJGECC7LuAyqIQquJy9SrFVvkp+EWrjYper0ulF6tYLVjbK1qLpW5VpAVFQZRV9n1fQgghLEnIMtlnO78/IrmmAfwM22D7ej4ePpwc5pz5nElyMq85y9gsy7IEAAAAAAbsoR4AAAAAgB8PAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDFnqAdwvvz+gAoLy0M9DAAAAOBHLTExxuh+7IEAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGAs6IDIzs7WM888oxtvvFGdO3dWenq68byffvqprr/+eqWlpSk9PV0LFy4M9uEBAAAAhFDQAbF//34tW7ZMrVu3VkpKivF8X3zxhaZMmaLhw4dr1qxZ6tevnx555BGtXLky2CEAAAAACBGbZVlWMDMEAgHZ7TXdMWXKFO3YsUMLFiz4wfluuOEGXXHFFfrd735XO+3ee+9VaWmpPvrooyCH/X/8/oAKC8vPeX4AAAAAUmJijNH9gt4DcSoegpGTk6ODBw/WO9wpPT1d27ZtU2FhYdDLBAAAAHDpXZKTqA8ePChJateuXZ3pKSkpsiyr9t8BAAAAXN6cl+JBSkpKJEmxsbF1pjds2LDOv58rp5OLSUnSko1HNG/5QV3Xs6XsdpsWr8/RTwa1VV5xpdbsOK47h1+h7ZkntetQoe65oaOWbTmq7OOlum9MZ/1jZZYK3VWaMOZKvffVPlV7/bovvbP+8uVeOR12/efwDnp7wW7Fx4Trpqvb6Z3Pd6tlYrSu7dFC7365V51ax+uq9o310TcH1KdTkpokRGnBqkMa2r2FHA6bvtmUq1H9W+tkSZXW7jqhW4emKDO3RDuyCvWfw6/Qmp0nlH28VD8dmapF6w/rpLta40amat6Kg6r2+nXnsA76aMkBOZ123TIkRX9dvFfxMRG6vk+y3l+yT62SYtS3cxN99PUBdWwdry5tE7Tg20PqkZqoxg0jtHj9EQ1MayqH3aaV247pul4tVVrh1aa9+bq+X7KOFpRrT3aRRvdvo/25Jco5Uaob+rXWpn35KnRX6fq+rbVy2zFVe/0a3qullm05qjCnXQPTmmnF1qOKiwlXl3YJWr3juJo1aqB2zWO1J7tIzRo3UFx0uCqqfIoMd8hms/3g9/HYyXK9+uFWNWzg0i/v7KYIV/1fU38gIIfdroKSKhW5q9ShVZyyjrlVWuFV15RGWr8nT4GApT6dkrRy2zFFhDvVrX1jLVp3WI3jItU1pZG+XHtYyU1j1L5FQy1ad1gdW8crKT5SX204oh5XJCrC5dDSLUc1KK2Zqr1+rd5xXMN6t1R+UaW27C9Q+sA22p9Tor05RRo7OEUrth3V0fxy/cewDvpi7WG5yzz6j2Ht9emyg/L5Ld12TYre+2q/Il0OjRnYRn9dvE+JcZG6tkdLzflyr9o2i1XPjkn6cMl+dWnXSCnNY/Xp8oPq16WpYhu49MWabF3Xq5Wqqn1avvWoxgxso6MF5dq4N1+3X9teWw8UaO/hYo0bmapvNuUqt6BM42/oqE+XHVRJuUf33NBRf120T76ApZ9en6p3Pt+tSJdTt12TolnzdykxLlI39GutWfN3qm2zWF3dtbn+vHC30lIaqXPrBH2wZL8GpDVTQmy4Fqw6pBF9kuXx+rV0c67GDknR4ROl2rA3T+NGpqpPpyYX41cc/2YCFSUqWzBdsgKKHHyPKpa+LZsrUpH971DF17Nkj01UePd0VXwzS86kdnJdMVAVy2crLLmrnM07qnL1h3J1HCR7VENVbVqgiKuuly0iRtV7liuy140KlBbIe2izIvrcIt/RvfId36/IfrfLs3+1AkXHFDHgdnm2L1GgqlSRfW9T1aYFkgKK6H2LqrculC28gcI7DVHV1oVyxCYprE13VW1bJGdiGzmatFP1zqUKa95R9thEefatUlhyV8kVKe+BdQpL6S1ZAXlztsuV0leWp0K+4wfkatdLgdIC+QtzFdaupwKlBQqUFSqsRSf5C3NleavkbJIif/FxSZYccc0UqHTLZnfKFh6lQFW5bM4w2ZwuWd4qyemSzWaX5fPI5nRJkizLMtoWn+LN3ipbRLScTczP+wyGFQjIZrfLCvilgL9m7H6vZFmyOV0KVJXJZrPLFh4lX8Fh2cMbyB7TSL4TB2WPTpAtKla+nJ1yNGopW0S0vIe314zVbpf30FaFtepS+1yHte0hq8Itb+5uua7or0DhUfnyMhXeaYh8x/bJX3RU4V2uk/fgBgUqShSeNlzVu5ZKfq9cV15X830PC5er49WqWvM32WMayXXFQFWu+UiOxDYKa91NlevmKqxlZzkat1bVpvlypfSWLTJW1Vu/lKvzUMlXrerdyxTRbbS8hzbLc3C9IgfeLd+hTfIe3aOooT+XZ8dX8hfmKuraCapaN1eByhJFXXOfKlfOkfx+RQ65RxVL/yRbWLiiBv6nype8KXtsY0WPerT2+4yL65IExCn//At76vSLYH6R/5ndblN8fIPzGte/ihXbjqmk3KNvNufKYbeppNyjFduO6UhemXz+gFZsO6qdB2sOF1u5/bg27cuvvb07u6hmGTuO6eBRtyRp+bb/u710y1Hl5JUpJ69MsdHhtbf9smpvHy+qVPbxUh3JK1PrZrE6dLxUf1+ZJYfdpoKSKi34Nrt2LJ+vydaO78ayYHW21u86UXN7TbZWbj0qSYqPDdfKbcckSdFRLq357j5Op10b9tSMvcob0IY9+dqwJ1/HTlZo0758bdqXr25XJGrLvnxt2JOn1s1itT+nWHtzimW325RXWKGjhRU6ml+mymq/iss92nagQP6ApbJqn77ddkyWJZVV+fXV+sOSam4vXH1IklRS7tVn32ZJko4XV2r+iizZbNLogW21YGWWnA670ge21bzlmUqIDdfoge30l4W71adzU/33vX1rv18b95xQ5pESpQ9qq+MnK1RQUqk+nZvq681HdbywQscLK1RQ6lVa+4b66xd7tGJLru4Z3UnzV2Zp58GTemBsV/15wU6VV/n0nyNT9eHiffIHLN16XQd9vGS/JCl3cDvNW16zh2/UgDb6/NuadRjWu5W+Wp8ju00a0LW5Vm49qoVrD+vKto20ZX++Vm0/psT4KB3MLdHew0WqrPar0F2lE8WVyjxSIp8/oIpqvzbty5NlSf6A9O32mu+VK9yp5ZtzJUkREWFatqXm+xke7tS6776Hstm0eV+BJKm82q8dWYXakVWoo0UV2p1dpN3ZRbqybSPtzSnWoROlSoqPUs6JUpVUeFVV7ZO73KPPVmfrYG6JApa0cG2Otu6v+Zn4Yl2O1u48LklavP6INu+veZwlm49qb06xJOmbzUd16FipJGnZ1mM6drJCx05WKLqBS4XuahW68yWbTcVlHq3Yekwn3dUqKffo601HlBQfpZJyj5ZtyVWVxy93uUcrtx9T5pFiBSxpza48jRxQd28rcC7Kjm1RSUF2zRf7lylQVPO7ZO1bqoA7TwF3nlwNomWVFcpbViinzS+roliePctlK8uTVVkiz44lcsY2klXplnfPUgWqKxWoLJV3099VfSxTsgJy2m2qzNoiSQpz2lW1d03NbZdDVTtXSJJcTpuqt31Tc9vmU9WWJZIkR+lxVW39SpIU3XWoqrYtlewONUjtq6rd36o6LEIRyZ1VmblJnh1fyRmToOqj++XPWqdAdYV8Rcel3B2qPnZAgQq37Pn7VL53razqCsX2Tlfplq9keasUN/g/5F75iRTwKeHan8q99K+SZanx9f+fiha9JbsrQo2G3aP8z/4gZ0yC4gbdpoLP/yhXkzaK6XatTn7xtiLbpikypbuKlvxFMV2vkbNhYxV/+6ni+t8smytcpVu/VsKQO+WvcKviwEbFD7lTnrxsFX02Q7LZ1eT2J1Wxd60aXNFHUR16SpKqjx/Usfd+JUd0nJr/9NeyqitUun2ZGlzRR1U5u+Te9KXiBoxVZfZOle9apYQRP1fF7tWqzN6hxNEPqnjN3+UtyFHiTyarcMm78peXKPHGh1Sw8E1Zfp8Sb7hf+Qt+L5vDqfihd6r0y7dlc0UobsDNKl36vuxRsYpJG6KytfPlaBCniHbdVL59qcISmsnZMEmVWVsV3ixFlt8nT162IrM3qvr4wZrnumC/KvZvlOXzyFF6VGXbl0tWQI7KApVvXixJclYVqnLTlzW3KwtU9d10R+kxVe+s+Xmwu4/Js+9bad+3suXtl/fAenkPrFVE6yvlzd4hX852hcU3lTcvW1bREVl+n/zlxar2Vsmbd0iyAgrsXy7Poe01P9/7l8ubvfW7n/Vl8h3b+93vwHL58w997/eh5u9MYO+y2t+HKE+BIhI7XIhfP/yASxIQ39/T0Lhx49rpbnfNi9N/3jMRjEDAkttdcX4D/Bcxul+y5q86pGt6tJDNZtPXm45o9Hfv+q/ZdVyj+rZWu2ax2plVqOv7tFJsVJgOnyjV9X1aSlZAJ0uqdH3vVior96ja69eIXi11LL9MYU6HhvdsqZwTpYqLDtd1PVoo+5hbLROjNfiqZjp01K1OreOV1q6Rck6Uqk+nJDVr1EAniyt1TfcWkk1avD5Hg9KaKq+oUqt3HtegtGaKigjTzqxC9evcRJWVXmWfKFXfTkk6ml+mk+4q9emUpL3ZRar2+tWnY5K27s9XmMOu3qmJ2rA7T/Ex4erRoZHW7jymlonRuqp9I63ffUKpyXHq2KqhtuzLV6c28WreOFr7c4p1RcuGstmkvMIKpTSPVUSYQzuzCtWmSbTyiyp0JL9cLRrV7DEoKq1WUly4oiKcqqzyKbFhuFxOu/wBS41iXLLbbLLbpQbhNb9CToe99nhAp8OmwpJKSVJZpVebdte8aN6874Te+NsWHS0o09ghKXrhnfXyBywdPlaiVduPyee3dM8NHRUZ7lC75rGKjgzTe1/uVpd9jfS3bw7IH7D0ydf7tedwsSRp7fajqqz2S5LyC8sV+C7IfR6fbDZJlhQR5pBNNaEd1yBMkhQV4VSj2AhJUsPocCU2DJckNW4YoWaNIrVlv9SicQM1bdRAB3NLlJwUI38goHXuKrVrFiuHTdp1qEjtW8TqZElNNKYmxynnRKmOF1boqnYJ2ptdpNIKj3p0aKRt+/MVsCz17NBYm/fmKcLlUO+Oidq4J0+N4yLUp2OituzNU9tmserXqYl2HyxUl5RG6tI2QQePFmtgWjM1ahih+cWVGtilqfwBS99sOqJBac3UukmMNuzN06C0JmoQ7tCew0W6umtTebw+Hc0v1+BuzVRSVi13uUfX9WiugqIK+f0BDevZQkfzyxQZ7tDwXi2Um1+mxg0jdH2fVsorrFDbZrHq2zlJRe4qpbVrpPYtG6qiyqsBXZqqYYNwLVybrWE9W8nrD2jZllyN6pes3PxG2rg3X8N6tlBRERd2wPmzGqUqrE13WYGA7Gmj5SwulC08Uva0dDmLC2SPTZK9yzA5ik/KmdROzvZ95CgtqdkD0SRFvuoqua4YIHtkrAJbv5Arbbj8J4+oetc3srfvr7DwWHmzt0jJPeSoKJM/P0tqeZXseTkKuPNktbxKtuydsqorZDVPk23vOkmSv3EHyb5UtrBw+Ru2kGSTLSJavqgkSZItsqF8EXE1t6MayhcWXbNCkQ0VCKt5wy/gipblD0iSvLYwWd9tQT3egKxAzfTqyipZPo8kqbKoSAr4JEnlBSekQM22z314n+T3KVBZpuJ9W6SAX76SfBXv2SBZAXmOH1TJ7hhJlioPbZenslKyAirbs1q2yFhZPo9Kti2VVVkqy1OhglWfyn88U5Ilr18Ka9Oj9vtRsPxj+Y/tVemO5Yq6+qfyHFgje2ySApWlClSW6mTmXlWt/1S+3F0q2bxElt8rq7xIJ1fOlb8wVwr4VbxpiXy5uyVJRVu+kfd4Zs3tbavkK6l5E6RoxxoFKmpeHxXvqXmBb/k8Ks2peTPI8lSpPK/mDZtAVbkqi4trvi+eKlWX12x7fNVVsvzf/U0ISLLs3922SWGRktzy2sKl8AaSzyOvM0a28ChZVWXyhcdLYRGSt1q+Bk0kp0sK+OVv2EqyOyWHU4GkVMmxUvbIWFnNukh7VssR30JW8y7SgQ1yNrtCapEmZe+QM/kq2eObS3mH5WjdQ7Lb5N/xtZzt+sjRpqe8BzfI2XW0wuNayXdsr+xpo+RyRstfeET2tFEK8wVkVZTU3C4rlRXwy5Y2Ws6SQtnCImTvOkrO4nzZYxqpIqKJKtn+nhfTN+WDvgrT95lehSknJ0fDhg3T66+/ruHDh9dO//TTT/Xkk0/q22+/VUJCwjmNgaswIWBZsn+3F+v7t79/2FBphUfRkTUvoEsrvYqNcikQsFRS7lF8TLg83pp3khvHRaqiyit3hVdNE6LkLveo0uNTk/goFZVWy+cPKDEuUscLK+Sw25QYF6mDR91qEOlUUlykdmYVKjEuUg2jXVq57ZjaNW8ol9Ouz9dmq03TGH2w5IAkaXivllq/J0/FZR7dOKitFnx7SP6ApQFXNtG3O0/I6bCpS9tG2nKg5t3z6/u00qb9Bbr56nbKPuFWXlGVbhnSTsWl1TpWWKGruzbXoeNuucu96pmaqMMnSmVZUuumMTp8olSuMIeaJkTp2MlyRUWEqWEDl3ILyhUX7VKDiDDl5pepUcMIRbicOnayXI0bRtaEkLta8bE1gVFa4VXDBjXPW0W1T9GRYQoELFV5fIqKqHluTx1a9c/fCwCXp+8fzmP5fbI5nDVHBwT8NbcDfskKyOYIq3tYTUWxbI4w2cIbKODOk80VJVtEtPyFObJHxUvhDRTIPyh7w6ZSWIT8Jw7I0ShZsjtqbjdJkXxe+fMy5WjeSValW/78LDlbd6t5N7kwV862PeU/cUCWO0/ODgPkP7JDVlWZnB0GyLev5hLwzpS+8myeL1t4A4WlXq3qDXNlj20iZ9ueql4/V44mKXI06yjPxnlyJneVLTZJns0LFNahv+QIk3fHYoVdOUyBolx5965QeO9b5M1cK3/Odjk79Jc/d7cczTsprNNQ+XK2y7v5H3I07yR/Xqbk88jepIPsUQ1rXtSGR8uqKpX/8BY5mqXK0SpN3p1fy9XjJ5LfK9/B9XL1Git/Xqb8R3crvN8d8h3aqMDJHLn63SHfvpWyyovk6n2rvNsWyvJ75eo5Vt5tn0vOcIV1HSnf7qWyNUiQs1WavLu+kT2hhRxNOsi7f5UcSSmyxybKm7lOzuYdZYuIke/ITjmad5RsNvmP7pGzZRdZfq8CBdlytOgkq6pcAfcJOZp0kFVZIquiWI7GbRSoKJY8VbLHNa25HfDLHt2o5vtud8oWES3LUyk5wmp+TnweyeH87nCxasnhks1mk+X3yuao+ftw6lAtXN5Mr8J0SQJCqrmMa8eOHfXqq6/WTuMyrvh34vX59fIHW3Q0v0w9U5PkdNg0oEszFZVWq6CkSolxkXKXV+svi/bJbrMpfUBrLVx7WG2bxqh54wbq27mJPll+UAeOlOiaHi00bkRqqFcJAH7Uqjf9Q/4jOxTe7w45kv7v0EPLCqjik2cVKMyRrUGC7LFJ8h/bK1ff2+Vs013Vaz6UP2ebwvv/p1xXXqfS2Q9K1eVyNOuo8D63yp7QSraw8BCuGXBuTAMi6EOYKisrtWzZMklSbm6uysrK9MUXX0iS+vTpo4SEBD311FOaN2+edu3aVTvfQw89pEceeUTJyckaMGCAlixZolWrVumtt94KdgjAj1KY06Gn7u6pnVmFeuXDLd9Ntemb784XuH9MZ207eFLdOzRWXLRL/1h1SC0ToxWwas5B2Zp5svYwpfyiytCsBAD8i7C81fJsmCtJqt66UGEpvWWPbixHUjsFio8prMsweXd8pbCOg1W9+j1JlvxH98iz4VPJ71HEiIdkC4tQ1cp3ZU9opcCxPXI07SBHk/ayAgFZAZ9sdqf8BYcku1P2qDhVr/1QtphEhff4SUjXHThfQQfEyZMn9fDDD9eZdurrd999V3379lUgEJDf769znxtuuEFVVVX64x//qLffflutW7fWq6++qkGDBp3H8IEfn2aNohQdGSaPz68WiTXHGtps0tbMAm09cFKSdFX7RpJqrsg0tHtzHcgtUZumMUof0EZbDxTo6q7NQzZ+APhXYAsLl7N9P/lzdsjmilTVV29Idqdc3dPl2ThPtpjGihhynzxbP5fzisGSr1rOdr3lz/nuBN+KElVveEdWVanszToq8qb/ljyV8pcVqnLer2T5qhXe9z9UvWK2ZLPJ2XGIfHtrTkp3Jl8lR+PWoVx94Lyc1yFMlwMOYcKPUbXXr0DAUmS4UwePuuWw21RZ7dOMT7apaaMo3f+TK/XVhiOKjgzT8F4tVeXxKy46XHY75xQAwIXm2fmVqlfNkRxhcqb0rTnHwuaQvWl7BY7tlRwuhfe+Wb6sTXK06ir/yWzJsmRVlSlwfK/CuoyQd/8qqbpczisGyrdvlSTJ2WmofLuXSpJc/e+UZ93HsjVIUIOxz8rmigrhGgOnd0nOgbgcEBD4V7J4fY7eX7JfHZPjFNvApXW7a65O9N/je4V6aADwL8uyLPkPb5UtppHsDRLk2b6o5lCm0nxVf/tXOVP6yZe5RrIs2Zt3VuBozSHazpS+ks0uW4N4ebcvkgI+hV01quYqUd5qhfe/Q77sLZIjTGFte9acbGx3cjIxLlsX7RwIABdWbkG5lm7KVe9OSdpzuObzOPYfKVHH1vGSpIpqXyiHBwD/8mw2m5ytu9V+Hd7r5trbYZ2vlc1mV5UzTN6DG2ouSeqtVODkYcnmkO/AtzXzXH2PZLMprMOA2isPSVJY+37/9zh8yBn+RRAQQIj95Ys92nekRBv35emnI1NV6fFpaLcWSk2O14Y9eUprd26XOAYAnD+brWZvQcTgn8sWESPPls8kV5Qa/PR1+XO2y3dog2zh0XK26SF75Ll/rhXwY0JAACGW3DRG+46UqHnjBvrD33fK6wuo5xVJatjApet6tgz18AAAp9Qe9W3Js+UzebcskC2xrRrc+F+y2XlJhX8f/LQDIfafw67Q8F6tJEn/NWuNpJpzewAAl4eqNR/If2SnXP3vlLO6XPJV1xzCJMkqypVsjhCPELi0OIkauIxkHXMrv7hSvTom8SnOAHAZsLxVKnvnAUmSo9VVtZdxdXYcIltYhJytusrRorMU8NU59wH4MTI9iZrLAACXkbbNYtWnUxPiAQAuE7awCIV1GipbTKLCOg+VLbZJzYf32GyySvOl8Aaq+PQ5lb3zgHzZm0M9XOCS4BAmIMR8/oD8AUvhYfV3ge/OLlJMVJhaJkaHYGQAAEmKuPoeWQGfvDsWy3XVKDladVHFexmSLFk+jwIF2ZIkX+4uOVt3D+1ggUuAgABCqKzSq+ffWafSCq8e/Y9uuqJVXO2/rdl5XG/O3yWnw6ap9/VVUjwfOgQAoeLds1zVaz6UJEXe+LTsLWo+DyIspa+cLa6UvyBbrq7Xh3iUwKVBQAAhdKKoQifd1ZKkzNwSXdEqTlUen1ZsPaaSco8kyR+w5PP/qE9VAoAfPXtMYs2hS85weTZ8qkDuTrl6jVVY6tWhHhpwyXESNRBi/1iVpWMF5TpeWKnGcRGKiXJp6eZcucLsunv4FUqMi1RqcnyohwkA//YCJcelsAiVfzhF8lbJ0eJKRY1+LNTDAi4YTqIGfiR+MrCtmidGK/tEqTbuzZc/UHMJ16hwp/p2bko8AMBlwt6wqexRcYq4ZoKc7fspvN9/hHpIQEhwCBNwGeiVmqj1u08oMtypK1vHq3/nJmrWqIHCnDQ+AFxunM06yuaKkj2+RaiHAoQEr06Ay0CzRg30yO3dlHWsVH/8xy7l5JerYXR4qIcFADiNigUvqnLBb1S98i+hHgoQEgQEcDn6UZ+ZBAD/2qzK0pr/V7lDPBIgNDiJGriMZB8vVUFJlXpc0Vg2PkwOAC5L/sIc+XN2yHnFQNkjY0M9HOCCMT2JmoAAAAAAwFWYAAAAAFx4BARwGflkWaZen7tdhe6qUA8FAADgtLiMK3CZOJJfps9WZ0uSmjWK0i1DUkI8IgAAgPrYAwFcJhLjItW6aYzCXQ51aZsQ6uEAAACcFidRA5cZy7K4AhMAALjkOIka+JEiHgDg8uPL3qKKBb+R9+A6WZ5K+Y7vlxUIhHpYQEhwDgQAAMAPqF73NwWKchVw58kT8ZkCBdkK63ytIgb9NNRDAy459kAAAAD8AGdKX8kRprD2/WWVF0lS7f+BfzecAwEAABAEf0G2fDnbFZY6SPaouFAPB7hgTM+B4BAmAACAIDgat5ajcetQDwMIGQ5hAgAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAuI5v35Wvh2mx5ff5QDwUAAOC0uAoTcJnIK67U63O3y5Lk81saM6BNqIcEAABQD3sggMtEpMuhqIiapm8cGxHi0QAAAJweHyQHXEZKyqpVWuFVy6ToUA8FAAD8m+GD5IAfoYbR4WoYHR7qYQAAAJwRhzABAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQACXkdU7j+vT5QdV5fGFeigAAACn5Qz1AADUOFFYoVnzd0mS7HabbhzUNsQjAgAAqI89EMBlIjoqTLENXLJJat64QaiHAwAAcFo2y7KsUA/ifPj9ARUWlod6GMAFUVHlVXmVT4lxkaEeCgAA+DeTmBhjdD8OYQIuI1ERYYqKCAv1MAAAAM6IQ5gAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAAAYIyAAAAAAGCMgAAAAABgjIAAAAIJQ9e1fVTr7QXn3rgj1UICQICAAAACC4N29TKoul3f/t6EeChASBAQAAEAQwvveJntSilzd0kM9FCAkbJZlWaEexPnw+wMqLCwP9TAAAACAH7XExBij+7EHAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAsaADIisrS/fee6+6deum/v37a+rUqaqqqvrB+Twej15++WUNGjRIXbt21a233qrVq1ef06ABAAAAhEZQAeF2uzV+/HiVl5drxowZeuKJJzR//nw9/fTTPzjv//zP/+i9997ThAkT9Pvf/16tWrXShAkTtHPnznMePAAAAIBLyxnMnT/44AO53W7NmzdPCQkJkiSHw6GMjAxNnDhRKSkpp53vxIkT+uijj/Tkk09q3LhxkqRBgwbpJz/5iV5//XX94Q9/OM/VAAAAAHApBLUHYvny5erfv39tPEjSyJEj5XK5tGzZsjPOt2fPHvn9fg0aNKh2ms1m06BBg7Ry5Up5PJ5zGDoAAACASy2ogMjMzKy3l8Hlcik5OVmZmZlnnO9UIISFhdWb1+Px6MiRI8EMAwAAAECIBHUIk9vtVmxsbL3psbGxKikpOeN8bdq0kSRt27ZNLVu2rJ2+ZcsWSTrrvCacTi4mBQAAAFwKQQXEmViWJZvNdsZ/79Chg/r06aPp06eradOmatu2rebOnav169dLkuz2cw8Au92m+PgG5zw/AAAAAHNBBURsbKzcbne96aWlpWc8gfqUF198UQ8//LDuvPNOSVKLFi00adIkzZw5U40bNw5mGHUEApbc7opznh8AAACAjN+UDyogUlJS6p3r4PF4dPjwYd1yyy1nnbdFixb6+OOPdeTIEVVVValt27Z65513lJiYqBYtWgQzjHp8vsB5zQ8AAADATFDHDg0ePFhr1qxRUVFR7bTFixfL4/FoyJAhRsto2bKl2rdvL6/Xq48//li33XZbcCMGAAAAEDI2y7Is0zu73W6lp6fXHn508uRJvfjiixo0aJCmT59ee7+nnnpK8+bN065du2qnzZkzR9HR0WrWrJlyc3P1zjvvyG636/3331dUVNQ5r4DfH1BhYfk5zw8AAABASkyMMbpf0OdAzJ49W1OnTtXkyZMVERGh9PR0ZWRk1LlfIBCQ3++vM83j8ej111/X8ePHFRcXpxEjRujhhx8+r3gAAAAAcGkFtQficsQeCAAAAOD8me6B4AMUAAAAABgjIAAAAAAYIyAAAACCULn0bZW+da88O5eEeihASBAQAAAAQfAdXCsF/PJlbQj1UICQICAAAACCEDFwnBwtu8jVa2yohwKEBFdhAgAAAMBVmAAAAABceAQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAEAQLMuSVV0e6mEAIUNAAAAABKFq0QyVzX5Q1ZsXhHooQEgQEAAAAEHwHd0tSfJ/93/g340z1AMAAAD4MYm4ZoJ8BzfIddWoUA8FCAmbZVlWqAdxPvz+gAoLOQ4RAAAAOB+JiTFG9+MQJgAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggI4DJSWe1Tobsq1MMAAAA4Iz4HArhMVFR59fRba1VS5tHEm7qoV8ekUA8JAACgHvZAAJeJ0gqviss8siQdyS8L9XAAAABOiz0QwGWiSUKU7h3dSccLKzSyT3KohwMAAHBafBI1AAAAAD6JGgAAAMCFR0AAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMEZAAAAAADBGQAAAAAAwRkAAAAAAMBZ0QGRlZenee+9Vt27d1L9/f02dOlVVVVU/OF9FRYWmT5+uYcOG6aqrrtKIESM0c+ZMeTyecxo4AAAAgEvPGcyd3W63xo8fr+bNm2vGjBkqLCzUtGnTVFxcrOnTp5913ueee05fffWVHnnkEXXo0EHbtm3TjBkzVFJSoqeffvq8VgIAAADApRFUQHzwwQdyu92aN2+eEhISJEkOh0MZGRmaOHGiUlJSTjufz+fTF198ofvuu0/jxo2TJPXr109Hjx7V559/TkAAAAAAPxJBHcK0fPly9e/fvzYeJGnkyJFyuVxatmzZGeezLEt+v18xMTF1psfGxsqyrCCHDAAAACBUggqIzMzMensZXC6XkpOTlZmZecb5wsLCNHbsWP3lL3/R1q1bVV5erjVr1uijjz7SXXfddW4jBwAAAHDJBX0ORGxsbL3psbGxKikpOeu8zz33nJ599lndfvvttdPGjRun//f//l8wQzgtp5OLSQEAAACXQlABcSaWZclms531PtOnT9fSpUv1wgsvqG3bttq5c6dmzJih2NhYPfTQQ+f82Ha7TfHxDc55fgAAAADmggqI2NhYud3uetNLS0vPeAK1JO3bt09/+tOf9MYbb+i6666TJPXu3Vs2m00vvfSS7rrrLjVq1CjIodcIBCy53RXnNC8AAACAGqZvygcVECkpKfXOdfB4PDp8+LBuueWWM8534MABSVKnTp3qTO/UqZN8Pp9yc3PPOSAkyecLnPO8AAAAAMwFdfLA4MGDtWbNGhUVFdVOW7x4sTwej4YMGXLG+Vq0aCFJ2rlzZ53pO3bskCS1bNkymGEAAAAACJGg9kDccccdmjNnjiZNmqRJkybp5MmTevHFFzVmzJg6hzA99dRTmjdvnnbt2iVJ6tKli7p27apnn31WBQUFatu2rbZv36433nhDo0aNqnNZWAAAAACXr6DPgZg9e7amTp2qyZMnKyIiQunp6crIyKhzv0AgIL/fX/u1w+HQH//4R/3ud7/TrFmzVFBQoGbNmunuu+/WAw88cGHWBAAAAMBFZ7N+5J/k5vcHVFhYHuphAAAAAD9qiYkxP3wnBXkOBAAAAIB/bwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjBAQAAAAAYwQEAAAAAGMEBAAAAABjzmBnyMrK0tSpU7Vx40ZFRkZq9OjRysjIUERExBnnOXLkiK677rrT/ltYWJh27NgR7DAAAAAAhEBQAeF2uzV+/Hg1b95cM2bMUGFhoaZNm6bi4mJNnz79jPMlJSXpww8/rDPNsixNmDBBffv2PbeRAwAAALjkggqIDz74QG63W/PmzVNCQoIkyeFwKCMjQxMnTlRKSspp53O5XOrWrVudaWvXrlVpaanS09PPbeQAAAAALrmgzoFYvny5+vfvXxsPkjRy5Ei5XC4tW7YsqAdesGCBoqOjde211wY1HwAAAIDQCSogMjMz6+1lcLlcSk5OVmZmpvFyvF6vFi1apOHDhys8PDyYIQAAAAAIoaDPgYiNja03PTY2ViUlJcbLWb58uYqLiy/Y4UtOJxeTAgAAAC6FoK/CdDqWZclmsxnff/78+WrcuLH69+9/3o9tt9sUH9/gvJcDAAAA4IcFFRCxsbFyu931ppeWlp7xBOp/Vl5erqVLl+rWW2+Vw+EI5uFPKxCw5HZXnPdyAAAAgH9npm/KBxUQKSkp9c518Hg8Onz4sG655RajZSxevFiVlZUaM2ZMMA99Vj5f4IItCwAAAMCZBXXywODBg7VmzRoVFRXVTlu8eLE8Ho+GDBlitIwFCxYoOTlZV111VXAjBQAAABByQQXEHXfcoZiYGE2aNEkrVqzQvHnz9MILL2jMmDF1DmF66qmn1Llz53rzFxYWavXq1Ro9evT5jxwAAADAJRdUQMTGxmr27NmKiorS5MmT9eKLLyo9PV1Tp06tc79AICC/319v/oULF8rn813Qw5eAfyVLNh7RnEV7VVbpDfVQAAAATstmWZYV6kGcD78/oMLC8lAPAzhvx06W679mrZUkpQ9oo7GD24V4RAAA4N9JYmKM0f34AAXgMhEXHa7GDSPksNuU0rz+560AAABcDtgDAVxGvL6APD6/GkSEhXooAADg34zpHogL8kFyAC6MMKddYXyyOgAAuIzxSgUAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgDECAgAAAIAxAgIAAACAMQICAAAAgLGgAyIrK0v33nuvunXrpv79+2vq1Kmqqqoymre4uFjPPfecBg0apLS0NI0cOVIffPBB0IMGAAAAEBrOYO7sdrs1fvx4NW/eXDNmzFBhYaGmTZum4uJiTZ8+/azzlpeXa9y4cQoPD9dTTz2lRo0aKTs7W16v97xWAAAAAMClE1RAfPDBB3K73Zo3b54SEhIkSQ6HQxkZGZo4caJSUlLOOO///u//qqqqSn/7298UEREhSerbt+95DB0AAADApRbUIUzLly9X//79a+NBkkaOHCmXy6Vly5addd5PPvlEt956a208AAAAAPjxCSogMjMz6+1lcLlcSk5OVmZm5hnny8nJUUFBgWJjY3X//ferS5cu6tu3r55//nnj8ycAAAAAhF7Q50DExsbWmx4bG6uSkpIzzldQUCBJeumll3T99ddr1qxZOnDggH7729/K6/Vq6tSpQQ67LqeTi0kBAAAAl0JQAXEmlmXJZrOd8d8DgYAkKSUlRdOmTZMk9e/fXz6fTy+99JIefvhhJSYmntNj2+02xcc3OKd5AQAAAAQnqICIjY2V2+2uN720tPSsJ1DHxcVJkvr161dner9+/RQIBJSZmXnOAREIWHK7K85pXgAAAAA1TN+UDyogUlJS6p3r4PF4dPjwYd1yyy1nnK9Vq1YKCwurN92yLEmS3X5+hyD5fIHzmh8AAACAmaBeuQ8ePFhr1qxRUVFR7bTFixfL4/FoyJAhZ5zP5XJp4MCBWr16dZ3pq1evltPpVPv27YMcNgAAAIBQCCog7rjjDsXExGjSpElasWKF5s2bpxdeeEFjxoypcwjTU089pc6dO9eZ98EHH9TevXv1+OOPa+XKlfrzn/+smTNn6q677qpzWVgAAAAAly+bdeo4IkNZWVmaOnWqNm7cqIiICKWnpysjI6PO5ztMmTJFn376qfbu3Vtn3lWrVumVV17Rvn37FBcXp5tuukkPP/zwaQ9vMuX3B1RYWH7O8wMAAACQEhNjjO4XdEBcbggIAAAA4PyZBgQfoAAAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMAYAQEAAADAGAEBAAAAwBgBAQAAAMCYM9gZsrKyNHXqVG3cuFGRkZEaPXq0MjIyFBERcdb5xo0bp3Xr1tWb/vnnnyslJSXYYQAAAAAIgaACwu12a/z48WrevLlmzJihwsJCTZs2TcXFxZo+ffoPzt+jRw898cQTdaa1bNkyuBEDAAAACJmgAuKDDz6Q2+3WvHnzlJCQIElyOBzKyMjQxIkTf3BPQmxsrLp163bOgwUAAAAQWkGdA7F8+XL179+/Nh4kaeTIkXK5XFq2bNkFHxwAAACAy0tQAZGZmVlvL4PL5VJycrIyMzN/cP5169apW7duSktL0913363169cHN1oAAAAAIRX0ORCxsbH1psfGxqqkpOSs8/bu3Vs33nij2rRpo7y8PL399tv62c9+pr/85S/q3r17cKP+J04nF5MCAAAALoWgr8J0OpZlyWaznfU+Dz30UJ2vhw4dqvT0dL3xxhuaNWvWOT+23W5TfHyDc54fAAAAgLmgAiI2NlZut7ve9NLS0qAvxRoVFaUhQ4boyy+/DGq+fxYIWHK7K85rGQAAAMC/O9M35YMKiJSUlHrnOng8Hh0+fFi33HJLMIuSVLPn4kLw+QIXZDkAAAAAzi6okwcGDx6sNWvWqKioqHba4sWL5fF4NGTIkKAeuKKiQsuWLVNaWlpQ8wEAAAAInaAC4o477lBMTIwmTZqkFStWaN68eXrhhRc0ZsyYOocwPfXUU+rcuXPt1xs2bNDEiRM1d+5crVmzRv/4xz901113KT8/Xw8++OCFWxvgX8CF2jMHAABwMQR9DsTs2bM1depUTZ48WREREUpPT1dGRkad+wUCAfn9/tqvExMT5fF49Nvf/lbFxcWKjIxU9+7d9fzzz6tr164XZk2AH7lqr18vztmk40UV+sWtXZWaHB/qIQEAANRjs37kb3f6/QEVFpaHehjAeTuSX6Zn3l4nSRrdv7VuGRLchQkAAADOR2JijNH9LshlXAGcv5aJ0Rrdv7WOn6zQtT1ahno4AAAAp8UeCAAAAADGeyD4CGcAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGCAgAAAAAxggIAAAAAMYICAAAAADGbJZlWaEexPmwLEuBwI96FQAAAICQczjM9i386AMCAAAAwKXDIUwAAAAAjBEQAAAAAIwREAAAAACMERAAAAAAjBEQAAAAAIwREAAAAACMERAAAAAAjBEQAAAAAIwREAAAAACMERAAAAAAjBEQAAAAAIwREAAAAACMERC4JG6++WalpqZq7dq1xvOkp6drypQptV9PmTJF6enpF3xsc+fO1fz58y/oMu+//36NGzfurPeZOXOmUlNTa//r27ev7rzzTi1btuyCjsXEPz/Xp3Mhn//CwkKlpqZq7ty5F2R5wL8itpv1ndpuXn311QoEAvX+fcKECUpNTdX9998f9ONfjHX6Zybfj7lz5yo1NVWFhYVBL3/VqlUaM2aMunTpol69ep3rMC+ZI0eOaObMmTpx4kS96ampqfriiy9CNDL8EGeoB4B/fZmZmdq1a5ckaf78+erbt+85LWfSpEmqqKi4kEOTJH366aeKiorSmDFjLviyf0hERIRmz54tScrLy9Obb76pBx54QH/961/Vo0ePSz4eAJcHtptnFhYWpqKiIq1du1b9+/evnV5YWKhvv/1WUVFR57TcUK7ThfLEE0+oU6dOevbZZxUeHh7q4fyg3Nxcvf766xo6dKiaNGlSOz0pKUkffvih2rRpE7rB4azYA4GLbv78+XI4HOrfv7++/PJLeTyec1pOcnKyOnbseIFHF1p2u13dunVTt27dNGLECL3xxhuyLEvz5s0L9dAAhBDbzTMLCwvT4MGDtWDBgjrTFy5cqKSkJF155ZUhGlloud1u5efna+TIkerVq5fS0tLOeVmWZZ3zz9yF4HK51K1bN8XFxYVsDDg7AgIX3YIFC9SvXz/97Gc/k9vt1vLly+vdZ9OmTRo7dqzS0tKUnp5+2sN4/nnX78yZM9W9e/d69+vevbtmzpxZ+/XGjRt11113qWfPnurevbvGjBmjTz/9VJI0btw4rVu3TkuXLq09lOj78y5dulS33Xabunbtqn79+unZZ5+t925eZmam7r77bqWlpWnYsGHn9eI/KSlJCQkJOnr0aO20vLw8Pfnkk7ruuuvUtWtXjRgxQr/97W/rbdxTU1M1a9YszZgxQwMGDFDfvn315JNP1huvyXNtYu3atUpNTdXKlSv1y1/+Ut27d9c111yjWbNm1bvvRx99pGuvvVZXXXWVxo8fr8OHD592mXPnztWYMWOUlpamq6++Wq+++qp8Pp8kqaysTNdcc40eeuihOvM8//zz6t27t44fP35O6wFcjthunl16eroWLVpUZzu4YMECjRo1Sjabrc59Tdb5bOt07bXX6le/+lWdeb/44gulpqbqyJEjtdOmT5+uMWPGqHv37rr66qv16KOPKi8vL6j1OpMf2r7PnTtXvXv3liT913/9l1JTU2sPZauurtaLL76oq6++Wl26dNGYMWPqHap16udk2bJl+slPfqK0tDQtWbKk9rnbsWNH7ff0pptu0o4dO1RdXa1nn31Wffr00eDBg/XnP/+5zjI3b96sBx54QIMGDVK3bt1044031vk+r127Vj/96U8lSbfeemvt8y6d/hCmQCCgP/7xj7r22mvVpUsXjRgxot5jnhrvnj17dOedd+qqq65Senq6VqxYcd7fA9TFIUy4qLZs2aKcnBxNnDhRAwcOVHx8vP7xj39o2LBhtffJz8/Xvffeq9TUVL322mtyu916/vnnVV5eft6PX1ZWpvvvv189e/bUb3/7W7lcLh04cEBut1uS9Oyzz+qxxx5TRESEnnjiCUlS06ZNJdX8gXjkkUc0duxYTZ48Wfn5+XrllVfkdrv16quvSqrZMP/85z9XZGSkXnrpJUnSa6+9pvLycrVt2zbo8ZaXl6ukpETJycm104qKihQXF6cnn3xSsbGxOnTokGbOnKn8/HxNmzatzvx//etf1bNnT7344ovKysrSyy+/rEaNGikjI0PSxXmun3vuOd144436/e9/r0WLFmn69OlKTU3V4MGDJUnffPON/vu//1tjx47VqFGjtGPHDj366KP1lvPOO+/o5Zdf1vjx4zVlyhRlZmbq1Vdfld/vV0ZGhqKjozVt2jTdc889mjdvnm666SatWLFC7733nl555ZXa7xvwY8d284dde+21evrpp7V8+XINGzZMubm52rx5s5599llt27Yt6HU+2zqZOnnypO6//34lJSWpsLBQ77zzjsaNG6fPPvtMTuf5v9w62/Z96NCheuutt3Tfffdp4sSJGjp0qBISEiRJGRkZWrZsmX7xi1+oQ4cO+vzzz5WRkSG/36+bbrqpdvl5eXn69a9/rYkTJ6pp06Zq1qyZDhw4IK/Xq6eeekr33HOPGjVqpOnTp2vy5Mnq0aOHGjdurFdffVVLlizRtGnT1LVr19rDb48ePaoePXrozjvvlMvl0qZNm/T0009Lkm666SZdeeWVeuaZZ/SrX/1K06ZNU7t27c66/i+99JJmz56t+++/X7169dKqVas0bdo0lZeX68EHH6y9n9fr1WOPPaaf/vSnmjRpkt5880099NBD+vrrrxUfH3/e3wfUICBwUc2fP18ul0sjRoyQ0+nUDTfcoE8++URlZWWKjo6WJM2ePVs2m01vvvmmYmNjJUmJiYm69957z/vxs7KyVFpaqkcffbT2nY3vHzPbvn17RUdHKyoqSt26daudblmWXnrpJY0aNUq//vWva6c3btxY999/vyZNmqQOHTpo7ty5ysvL08KFC2uP1UxNTdWoUaOM/xCeeoc9Pz9f06dPV3R0dO27MqeWd+oPmiT16NFDkZGRmjJlip555hlFRkbWGd8rr7wiSRo8eLC2b9+uL7/8sjYgLsZzPWLECE2ePFmS1K9fPy1dulRffvllbUD84Q9/UK9evWpj5+qrr1ZlZaX+93//t3YZZWVlmjFjhu67777auBg4cKAcDodeeukl3XvvvYqPj1e/fv00fvx4TZ06VampqXrqqac0atSoi3KSKBAqbDd/WEREhIYNG6YFCxbU/j8lJeWcD9c60zoF4/tv6Pj9fnXv3l2DBw/WmjVrNGjQoHNa5vedbfuekJBQe+hWcnJy7Trs2bNHixYt0jPPPKO77rpLUs02OC8vTzNmzKgTECUlJXrrrbfUtWvXOo/r9XqVkZFRu00PBAJ64IEH1K1bNz355JOSarb9X3zxhb744ovagBg9enTtMizLUu/evXXixAl98MEHuummmxQdHa327dtLkjp06HDWQ64KCws1Z84c/exnP9MvfvELSdKgQYNUXl6ut956S/fcc48aNGhQZ7xDhgypfT5GjBih5cuX68YbbwzuSccZcQgTLhq/36+FCxdq6NChiomJkSSNGTNG1dXVWrRoUe39tm7dqr59+9b+EZRqNgyn/lCej+TkZEVHR+u5557T559/bnxVi6ysLOXm5uqGG26Qz+er/a93796y2WzasWOHJGnbtm3q0KFDnRO92rVrpw4dOhg9TkVFha688kpdeeWVGjp0qBYuXKiXXnqpzvIsy9Kf//xnjRo1Sl27dtWVV16pjIwM+Xw+5eTk1FnewIED63zdvn37Oof2XIzn+vt/GO12u9q1a1f7mH6/Xzt37tTw4cPrzDNy5Mg6X2/evFkVFRW6/vrr6zzf/fr1U1VVlfbv319730cffVRNmjTR7bffLqlmDwjwr4LtprkxY8bom2++UXl5uRYsWBDyk5+XLVumO+64Qz179lTnzp1rX3AfOnTogiz/h7bvp7Nx40ZJ0qhRo+pMHz16tHJzc3Xs2LHaafHx8fXiQarZrvfr16/261PftwEDBtROczgcSk5OrjOekpISTZ06Vddcc03t37kPP/xQWVlZP7Cm9W3btk1er/e061FRUaHdu3fXGe/3g7d169YKCwurd6UnnB/2QOCiWbVqlU6ePKlrrrmmdtd3+/bt1bRpU82fP19jx46VVPPOe+vWrevN36hRo/MeQ8OGDfXOO+9oxowZevzxx+X3+9WrVy89/fTTte+snU5RUZEk1dkt+n2nNrp5eXmnHWfjxo1r9yycTUREhObMmSPLsnTo0CG98sorevzxxzV//nwlJSVJqnmn8Te/+Y3uu+++2hcM27dv169+9StVV1fXWd73X0xINScbfv8Y4YvxXJ96kfP9xzx1XG5hYaF8Pl/trvRTGjduXOfrU8/3zTfffNrH+P4fufDwcA0fPlx/+MMfNGbMGDVs2PCcxw5cbthu/vB285QBAwaoQYMGeuONN7Rv3z698cYbxvNeaNu2bdOkSZN03XXXacKECWrUqJFsNptuv/32etvpc/VD2/fTKSkpkdPprHfozqltcElJiZo1aybpzD87ERERcrlcdR5XOv22//vrOmXKFG3evFkPPvhg7R6e999/XwsXLjzrmM+0HlLNXrbTrUdxcfEZx3u6seH8ERC4aE6dpPXkk0/W7uY8JS8vT/n5+UpMTFRiYqJOnjxZb/7TTfu+8PBweb3eOtM8Ho8qKyvrTOvataveeustVVVVae3atfrNb36jBx98UF999dUZl33qyg/PPPPMad+ROfXiPikpSTt37qz37wUFBUZXj7Db7bW7bbt27ap27drptttu0+9//3s9//zzkmqOKb722mv1y1/+sna+zMzMH1z26Zzrc32uEhIS5HQ6672DWVBQUOfrUxHw+uuvn/a445YtW9be3rdvn95++2117txZc+bM0dixY2t3gwM/dmw34846/u9zOBy64YYb9Kc//Undu3dXq1atTns/03U+E5fLVW/+Uy9oT/nqq68UHR2t1157TXZ7zcEdubm5pqty0TRs2FA+n0/FxcV1nttT2+DvvwHzzyefn4/q6motW7ZMTzzxRJ3P9njvvffOaXmnxl5QUFDncq+n1oOrNV16HMKEi6KyslJfffWVhg0bpnfffbfOf6+99poCgYA+++wzSTV/qNauXavS0tLa+VeuXKmysrKzPkaTJk3k9XrrXNHn22+/lWVZp71/RESEhgwZojvvvFNHjhypfTfidO9MtGvXTk2bNlVOTo7S0tLq/XdqA5aWlqb9+/fX2UV98ODBOofcBKNLly4aPXq05s6dq/z8fElSVVVV7Ts+p5zrhx2d63N9rhwOhzp37qzFixfXmf7ll1/W+frUeR3Hjx8/7fN96t0zj8ejxx9/XF26dNGHH36oK664Qo8//ni9P+7AjxHbzeC3m7feequuueYa3XPPPee9zmd6l7pp06b13rRZtWpVna9Pbae//yL8Yn8onYmePXtKUr13/T///HO1aNGidu/DhebxeOT3++v87SorK9PXX39d536n/v2H9g6kpaUpLCzstOsRFRWlzp07X6CRwxR7IHBRfP3116qoqNC4ceNO+wFIb7/9tubPn6977rlH48eP13vvvacJEyZowoQJcrvdmjlz5g++ozB48GBFRUXp6aef1oQJE3T8+HG9++67dTZYS5cu1ccff6xhw4apefPmKigo0Jw5c9SjR4/aD9lp166d5s2bp6+//lqJiYlKSkpSkyZNNGXKFGVkZKiiokJDhw5VZGSkjh49qmXLlumRRx5R27ZtNXbsWP3hD3/QAw88oF/84heyLEu/+93v6h2iE4xJkybps88+0+zZs5WRkaEBAwbo3Xff1Zw5c9SmTRvNnz9f2dnZ57Tsc32uz8cDDzygSZMm6cknn6y9CtM/X789JiZGDz30kF5++WUdP35cffv2ld1uV05OTu2lBCMjIzVz5kxlZ2fr73//u1wul1566SXdfPPN+v3vf197Yh3wY8V2M/jtZqdOnX7w0CWTdT7bOo0cOVLPPfecXn/9dXXv3l1Lly7V9u3b68w7cOBAzZ49Wy+88IKGDx+uzZs36+9//3vQ63OhdezYUSNHjtSLL76oqqoqtW/fXgsXLtSKFSv0m9/85qI9bkxMjNLS0jRr1qzaPdFvvvmmoqOj6+yRbtOmjRwOhz755BM5HA45nc7TnkydkJCgcePG6U9/+pNcLpd69Oih1atX68MPP9TkyZPP+cMDce7YA4GLYv78+WrevPkZPz315ptv1o4dO5SVlaWkpCTNmjVLVVVVevjhhzVr1iw988wz9Y51lOruYo2Pj9eMGTNUWFioBx98UB9//LFefvnlOn8UkpOTZbfb9dprr+nnP/+5pk2bph49euh3v/td7X0mTJigHj166IknntCtt96qjz76SJJ0ww036M0331RWVpZ++ctfatKkSXrnnXfUokWL2j90ERER+tOf/lR7Kb3p06drwoQJ5/VBRu3atdPo0aP1/vvvq7S0VA8++KDGjBmjGTNm6NFHH5XL5aq9FF6wgnmuL5TrrrtOzz//vFavXq0HH3xQ3377be2VRL7v1Pdn7dq1mjx5sh5++GF99NFHte88bdq0SW+99ZaeeOKJ2svctmvXTo899pjefPNNbd269aKtA3ApsN28OB8AZ7LOZ1un2267TT//+c/1/vvv6+GHH1Z1dbUefvjhOvMOGTJEGRkZWrJkiSZOnKgNGzbUudJcKL388su688479fbbb2vixInatm2bXn755TpXYLoYXnnlFbVq1UpTpkzR1KlTNXLkyHqPmZCQoGeeeUbr16/X3XffrVtvvfWMy3vsscc0efJkzZs3Tw888IAWLVqkKVOmnPGcG1xcNutM+y2By8zkyZNVVFSkOXPmhHooAPCjwHYTwMXAHghc9kpKSvTVV19p3bp1Z71ONACgBttNABcTAYHL3vr16/XYY4+pe/fumjBhQqiHAwCXPbabAC4mDmECAAAAYIw9EAAAAACMERAAAAAAjBEQAAAAAIwREAAAAACMERAAAAAAjBEQAAAAAIwREAAAAACMERAAAAAAjBEQAAAAAIz9/3gf2K1bIv7mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.catplot(compare_hdbscan_models(), height=8, kind=\"swarm\", size=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170f6c3e",
   "metadata": {},
   "source": [
    "It turns out that it is not uncommon for multiple points to have the same mutual reachability distance to a given point. Given such ties they can end up in a slightly different order in the tree. Mix that with floating point precision issues that can also subtley rearrange the order of points in the tree, and the result can be somewhat different condensed trees. So yes, in rare cases it is possible to end up with different clusterings between ``hdbscan`` and ``fast_hdbscan``. It is worth noting, however, that qualitiatively these clustering are, in fact, very similar -- it is usually a matter of mergeing or splitting a single cluster that was right on the cusp, or a change in noise points with respect to a cluster or two.\n",
    "\n",
    "So, yes results can differ. For practical purposes they are almost always effectively identical; and in the rare cases where that's not the case they are qualitiatively very very similar. So, unless you need to exactly replicate ``hdbscan`` the ``fast_hdbscan`` library should be quite sufficient for your needs."
   ]
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}