ljmartin/balanced_cut2.ipynb

## balanced_cut2.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.cluster.hierarchy import dendrogram as show_dendrogram\n",
    "from sklearn.metrics import homogeneity_score\n",
    "\n",
    "import networkx as nx\n",
    "import sknetwork as skn\n",
    "\n",
    "adjacency = skn.data.karate_club()\n",
    "paris = skn.hierarchy.Paris(engine='python')\n",
    "dendrogram = paris.fit_transform(adjacency)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Balanced cut:\n",
    "Straight cuts often have one large group with many small groups. \n",
    "\n",
    "As opposed to a straight cut, balanced cuts leave all clusters having around the same size. \n",
    "\n",
    "----------------\n",
    "\n",
    "Example with karate club:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def balanced_cut(dendrogram, max_cluster_size):\n",
    "    n_nodes = dendrogram.shape[0] + 1\n",
    "    labels = np.zeros(n_nodes, dtype=int)\n",
    "    cluster = {node: [node] for node in range(n_nodes)}\n",
    "    completed_clusters = list()\n",
    "    \n",
    "    for t in range(n_nodes - 1):\n",
    "        currentID = n_nodes+t\n",
    "        left = cluster[int(dendrogram[t][0])]\n",
    "        right = cluster[int(dendrogram[t][1])]\n",
    "        if len(left)+len(right) > max_cluster_size:\n",
    "            for clust in [left, right]:\n",
    "                if len(clust)<max_cluster_size:\n",
    "                    completed_clusters.append(clust)\n",
    "                    \n",
    "        cluster[currentID] = cluster.pop(int(dendrogram[t][0])) + cluster.pop(int(dendrogram[t][1]))\n",
    "    \n",
    "    for count, indices in enumerate(completed_clusters):\n",
    "        labels[indices]=count\n",
    "    return labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using max cluster size of 10\n",
      "Cluster ID 0 has 6 items\n",
      "Cluster ID 1 has 6 items\n",
      "Cluster ID 2 has 8 items\n",
      "Cluster ID 3 has 9 items\n",
      "Cluster ID 4 has 5 items\n"
     ]
    }
   ],
   "source": [
    "adjacency = skn.data.karate_club()\n",
    "paris = skn.hierarchy.Paris(engine='python')\n",
    "dendrogram = paris.fit_transform(adjacency)\n",
    "\n",
    "\n",
    "maxClusterSize = 10\n",
    "labels = balanced_cut(dendrogram, maxClusterSize)\n",
    "ids, counts = np.unique(labels, return_counts=True)\n",
    "print(f'Using max cluster size of {maxClusterSize}')\n",
    "for j,k in zip(ids, counts):\n",
    "    print(f'Cluster ID {j} has {k} items')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yN5//H8ddZSYgRNEaMIEasoCQpEbNWWi2htEpbSqlRsUvxM0trxGi/takatWLvvaJCEAQJiRiJlYQQWWfcvz9OnYok5IwEJ9fz8chDmtznPldUzvtc43NdMkmSJARBEAQhj5C/6QYIgiAIQm4SwScIgiDkKSL4BEEQhDxFBJ8gCIKQp4jgEwRBEPIUEXyCIAhCniKCTxAEQchTRPAJgiAIeYoIPkEQBCFPEcEnCIIg5Cki+ARBEIQ8RQSfIAiCkKeI4BMEQRDyFBF8giAIQp6ifNMNEARBEHLW0ycpHD8Qwc0b8SQ9U5PPXkXpsg54t3ChSNH8b7p5uU4mzuMTBEGwTjcj49m6/iLnz9xBJpOhTtMavqdSyZGA6m4ladepFlWqFX9zDc1lIvgEQRCs0NH91/lrYRBqtZZXvsrLwMZGQfsubvh0qIFMJsu1Nr4pYo5PEATByhw7EMFfC4NIS3tN6AFIkJaqZfPaC+zcFJor7XvTRPAJgiBYkdtRj1ix4BRpLwxrZsfz8Lsaej+HWvb2EMEnCIJgRXZuCkWj0Zn02LRULVvXXbBwi94+IvgEQRCsxLPENE6fvIVOZ/rSjfDLD4h7+MyCrXr7iOATBEGwEv8cu4G5a1MkCQ7vu2aZBr2lRPAJgiBYiTu3HpOWatzc3ss0Gh23ox5ZqEVvJxF8giAIViLpWZpF7pOcpLbIfd5WIvgEQRCshH0BW4vcJ7+9yiL3eVuJ4BMEQbASZco5YGOrMOseSpWcsuWLWKhFbycRfIIgCFbiA+/yry9Yfw0ZMpq2rGyZBr2lRPAJgiBYifz2Nng0dEYuN31pZ9UaxSn6nr0FW/X2EcEnCIJgRXx8a6BQmhZ8NrYKPvmsloVb9PYRwScIgmBFUtJiuRSxHpnMuN1bbGwV+H5Rm6o1SuRQy94eIvgEQRCsxKlTp/D29ubbvh/R64dG2NgqkL1m2FP27+kMnb6sS9v2NXKppW+WOIhWEATBCmzdupVvv/2W5cuX89FHHwFQrkJRtm+8xOnAKDQaNXLZf2UKKpUCkKhZ14mPO9akUlXHN9Ty3CfO4xMEQXjHzZ8/nwkTJrB161bc3d0zfL9/v0HYqypTxeV9EhNTsbe3walMYRq1cKGwQ7430OI3SwSfIAjCO0qSJMaMGcO6devYvXs3Li4umV7j4uLCli1bqFXL+heuZIcY6hQEQXgHpaWl0bt3b8LCwggMDMTRMfOhyvDwcNRqNTVr1szlFr69RPAJgiC8Y548eULHjh3Jnz8/Bw8eJH/+/Fleu3PnTnx8fJCZe2yDFRGrOgVBEN4hMTExNG7cmEqVKrFx48ZXhh7og69t27a51Lp3g5jjEwRBeEdcvnwZHx8f+vTpw48//vjaXlxiYiKlSpUiJiaGggUL5lIr335iqFMQBOEdcOzYMTp16sSMGTPo3r17th5z4MABPD09Rei9RASfIAjCW279+vX079+f1atX8+GHH2b7cbt27cLHxycHW/ZuEkOdgiAIb7HZs2czY8YMtm/fTp06dbL9OEmScHZ2Zu/evbi6uuZgC989oscnCILwFtLpdAwbNow9e/YQGBhIuXLljHp8aGgoCoWCqlWr5lAL310i+ARBEN4yKSkpfP3119y7d4/jx49TpIjxB8OKMoasiXIGQRCEt8ijR49o3bo1AHv27DEp9OC/4BMyEsEnCILwlrh16xZeXl7Uq1ePNWvWYGdnZ9J9EhISCA4OplmzZhZuoXUQwScIgvAWCAkJoWHDhvTu3ZtZs2Yhl5v+8rx//34aNWr02uL2vErM8QmCILxh+/fvp2vXrvz+++989tlnZt9PDHO+mihnEARBeIP++usvhg0bxoYNG/D29jb7fpIk4eTkxLFjx6hUqZIFWmh9RI9PEAThDZAkiWnTprFgwQIOHTpE9erVLXLf8+fPU7BgQRF6ryCCTxAEIZdptVoGDhxIYGAggYGBODk5WezeYpjz9cTiFkEQhFyUlJSEr68v169f5+jRoxYNPRDBlx0i+ARBEHLJw4cPad68OQ4ODmzfvp1ChQpZ9P7x8fFcvHiRxo0bW/S+1kYEnyAIQi6IiIjAy8uLFi1asHz5cmxsbCz+HHv37qVp06Ym1//lFSL4BEEQctjp06fx9vZmyJAhTJkyJce2EROHzmaPKGcQBEEwQuS1WPZuu8rtm49ITdFga6fEqUxhWrVzpVJVxwyhtmPHDnr06MGSJUto165djrVLp9NRsmRJgoKCKF++fI49jzUQqzoFQRCy4dTxKALWhPAoNok0tRZJ91+fIfrWY0LORFPIwY72n7vh1bQiMpmMRYsWMW7cOLZt24anp2eOtu/MmTM4OjqK0MsGEXyCIAivIEkSq5ac4ci+a6SlarO4BlJTNTy8n8if808RFnqfqPu7WbNmNUePHqVy5co53k5x6Gz2iTk+QRCEV1j759lXht7L0lK1HN0fzvmTyQQGBuZK6IEoYzCGCD5BEIQsXL5wlwO7wrIdev9RULJYXaKjUnOkXS978OABYWFheHl55crzvetE8AmCIGRh+8ZLJoSeXlqqlm0bLlm4RZnbs2cPzZs3z5ESCWsk5vgEQRAyEffwGeGXH5h1j1tRj7gX/YSSpS1XqK7V6gi//IAnCSloNDrs7W3YuWO/GOY0ggg+QRCETJw4HIG5tV46rY4j+6/T5ev3zW7P4/gkDu4OZ//OMLRaCZCQJJDLZCiSmvDwZinCQu9TpXrxHKsTtBYi+ARBEDJx/24iGrXOrHtotRL37z4xuy2H9oSzavEZQEKdSZsUchVhl2KZOekgFSoVw++nZuTLpzL7ea2VmOMTBEHIRGqy2jL3SdGY9fjtAZdYvfQMarU209B7TpL0zxUR9pCJw3eRYqH2WyMRfIIgCJkoUMjWIveJi7/HtWvX0OmM7z2eOXmTLX9fMGqBjVqt4+H9p8z++RBiY67MieATBEHIhEsVR2ztzJsNksl1RN0OoWXLlhQtWpQWLVowcuRI1q9fz40bN14ZTJIksXpJMGlpxq8qVat1RIbHEREea07zrZYIPkEQhEx4NHLG3NUtSqWKpSunEBUVxbVr1xg2bBgFChRg5cqVNGrUiPfee4/WrVszZswYNm/ezJ07dwxhGBb6gMRE0+sA09I07Np82bwfwEqJTaoFQRCy8NfCUxzac+3fVZTGkctlNGhSge8GZV1UHhMTQ3BwMGfOnOHMmTOcPn0auVxO/fr1KVPkY5Kf2pvTfJQqObOXdKRgIXFM0YvEqk5BEIQstPm0BscORqLVGr9ARamS83HHmq+8xsnJCScnJ8OpDZIkcefOHc6cOcOWvxJMavOLVCoFNyPjqVnHsqe8v+vEUKcgCEIWHEsUwG90U2xsFEY9zsZGQf9hjXEqU9iox8lkMsqWLUuHDh1QyM0vR5AkiaRnYnXny0TwCYIgvEJ1t1LICoQgoUb1mgBU2cixtVPyw6im1HEvY9bzKhTmvzzLZDKUKvEy/zIx1CkIgvAKK1eu5PDxTRw6eILN605z4uAtChYsjEajQ6vRoVDKUSrlKBRyWn7sSrPWVXAoks/s5y1Q0IbUVPNqACVJorCDmN97mQg+QRCELISGhjJ48GAOHjxI6TKOnLscgFujyvi0ase9mCckJ6nJl09F8VIFqVXXySK9tOeatKzMto2XUJtQzvCcra2SCpXes1ibrIVY1SkIgpCJp0+f4u7uzqhRo/j66695+vQp5cqV48qVK5QsWTLHnz/hcTJDewe8creWV7GxUdDhi9r4dKhh4Za9+8TgryAIwkskSeK7777D29ubr7/+GoC///6bZs2a5UroARR2yEetuqVRKEzbcFoCGreoZNlGWQkx1CkIQp5yNzqBo/sjuH/3CSkpGgoUtKVSVUcaNatIfnv9eXb/+9//uHr1KoGBgYbHLVq0iAkTJuRqW3v082SM30MSHicD2Q9AG1sF3w5oYLFt16yNGOoUBCFPOBd0m60bLnE76hFarQ7dC0XpNrYKJAk8GjrjXEVH9286cvLkSVxcXAAICQmhXbt23LhxA4XCuNIGc/08eTYX/7HFzrZgujZnxcZGwec96tGibdVcaN27SQSfIAhWTafVsWJhEIGHb7x2laRMLkOjTsWrVWG+H9jJ8PWBAwdSrFgxxo8fn8OtTW/NmjWMHDmSPbsPcWD7Hc6fjgYZGRa8yOQyVCo5xd6z58te7tSqKwrWX0UEnyAIVkuSJJb+fpJ/jkUZdcKBjY2CH0Y1pVZdJ5KTkylbtizBwcE4OzvnYGvTO3jwIJ9//jkHDhygVq1aADx9ksKRfddZtewwBewdkMkU2Nkpca1ZgjafVqdiZbGCMztE8AmCYLWOH4pgxfwgk+rhbO2U/PpHe7bv2MiqVavYtWtXDrQwcxcuXODDDz9k7dq1NGvWLN33JEmiaNGihIeH4+jomGttsiZiVacgCFZJkiS2rL1gchG4Tqfj4O4wFi9eTK9evSzcuqzdvn2bjz76iLlz52YIPYC7d++iVCpF6JlBBJ8gCFbpethDEh6lmPx4dZqOPVsvc/VKmGET6Zz26NEj2rZty+DBg/n8888zveby5cvUqCFq88whgk8QBKu0d/tV0tLM2/IrNTWNLh37YWNjY6FWZS0lJYX27dvTsmVLhgwZkuV1oaGhVK9ePcfbY81EHZ8gCFYp5nYC5q5g0Gh01H+/sWUa9Ao6nY6vvvqKEiVKMHPmzHTfS3yayqljUdy/+4SkJDUhp5JwrliHtDSt0adGCHoi+ARBsEqpKeb19gDkMiWFChaxQGtebejQody/f589e/Ygl+sH4qIi4ti56TJnT91GJueFValluBUmY8BX62jyYSVatauGY4kCOd5GayKCTxAEq2SXzxIvb1pCr1zA5ZKcqlWrolKZf0bey2bNmsW+ffs4duwYdnb6kxT2bLvChr/OodbokHQZu60atYRGreHg7nCO7L/OwJFNRO2eEUQ5gyAIVmmB/3FOHovKNDiySybXkiwL4tzF/dy5cwdXV1dq166d7qNo0aIm3//vv/9m+PDhBAYGUrZsWQB2bQ4lYE2I0XWHg0Y3FSetZ5MIPkEQrFLktVimjtlrVIC8rGAhW+Yu64RcIefZs2dcvHiRkJAQw8fFixcpXLhwhjCsVKnSa7c2O3z4MJ07d05XoH7l4j1mTT5oUptt7ZRMnfcJxRztTfpZ8xIRfIIgWK2R/bdwL/qJSY9VKmV82qU2n3xWK8trdDodUVFR6cIwJCSEBw8eUKNGjXRh6ObmRqFChQC4ePEiLVq0yFCgPmX0HsIvPzCxvXJatatGl6/fN+nxeYkIPkEQrNap41EsnhdoUg/KLp+S6fM7UKiw8SeYJyQkZOgdXrp0iRIlSlClShVOnjzJt99+y4ABAyhfvjxyuZz7d5/y0w/bUKtN76Hmy6fitxWfoVSJ1Z6vIoJPEASrtnLxaY7su2ZU+MnlEiMntsa1ZgmLtUOr1XL27Fl8fX2pUqUK+fLlIyQkhISEBNzc3HB1bocmqRSSZNr5e6AP6579G+DZqLzF2m2NxKpOQRCs2pff1sdGJWffzrDXhp9cISMtNRm/Mc0tGnoAGo2G4cOH07FjR/z9/ZHJ9AEXFxfHhQsX2PDnDdRmhB5ASrKGu3cSLNFcqyZ2bhEEwarJZDI6f12PIWOaU6N2KRRKGVqdOt01tnZKbO2UODgmYlPsHO4fuFi0Dc8L1B0dHZk1a5Yh9ACKFStGs2bNKFbUMkGb+DTVIvexZqLHJwhCnlCtVkmq1SrJHP8FnP0nGu9GbUhJVlOwkC0uVRxxb1gOt9o1Wb58ucWfe9iwYdy9e5e9e/caCtRfZmtnmZfjfPY5v73au04EnyAIeco/QYdp2bolPXs2TPf1o0ePolQq+eCDDyz6fP7+/uzZs4fjx48bCtQz41SmMNevPkCnM/25bO2UlChZ0PQb5BFiqFMQhDxDkiQOHz5M06ZNM3xv0aJF9O7dO90wpLnWrl3LrFmz2LVrF0WKvHrrsxY+VVEqzVuNKekk6jcoZ9Y98gIRfIIg5BnXrl1DqVRSoUKFdF9/9OgR27Zto1u3bhZ7rsOHDzNw4EB27NhBuXKvD6Ny5YtQwqmQyc8nV8ho2LQidvksv62atRHBJwhCnnHkyBGaNGmSoVe3atUq2rZty3vvvWeR57l06RKdO3fm77//xs3NLduP++SzmiAzrY5PqZTT+pNqJj02rxHBJwhCnvE8+F4kSRKLFi2y2Cnrd+7cwcfHh9mzZ9O8efNsPy4+Pp6xE/vwNDUMlY1xL802Ngq+7uuJU5nCxjY3TxLBJwhCniBJEkeOHMkwv3fmzBkSExPTbR1mqsePH9O2bVsGDhxI165ds/24ixcv4uHhQY0aNfh780SatqqCjW025vtk+tDr1tudRs0sW4JhzcTOLYIg5AkRERF4e3sTHR2dbqizT58+ODs7M3r0aLPun5qaSps2bahVqxZz5szJ9iKZDRs28P333+Pv759ujvFs0G22rrtI9K3HaDQ6dC+cMqFSyZGAGm6l+LSLGy5VLDNEm1eI4BMEIU9YunQp+/btY82aNYavJSYmUq5cOS5duoSTk+lH+uh0Orp27YpGo2Ht2rWvPZkB9FuYjRs3jpUrVxIQEEC9evUyvS7mdgIH94RzL/oJKclq8hewoUKlYjRrVRmHovlNbnNeJur4BEHIEzKb31u3bh3e3t5mhR7AiBEjiI6OZt++fdkKvcePH/Pll1/y7NkzTp8+TfHixbO81qlsYbr1cjerfUJ6Yo5PEIQ8IbP5vcWLF5u9qGX27Nns3LmTLVu2vLJA/bnLly/j4eGBi4sL+/bte2XoCTlDBJ8gCFYvKiqKlJQUqlatavhaaGgoN2/epG3btibfd926dcyYMYPdu3dn6yT2zZs306RJE0aPHs3cuXNRqUTN3ZsghjoFQbB6R44coXHjxukWnCxevJgePXqgVJr2MnjkyBEGDBjAvn37XlugrtPpmDhxIkuWLGHHjh14eHiY9JyCZYjgEwTB6r08v5eSksLKlSsJCgoy6X6hoaF07tyZNWvWULt27Vde++TJE7p37058fDynT5+mZMmSJj3nu06SJKQHN5Aex0BaMqjskBVxQla8okW3icsOEXyCIFi9I0eOMHToUMN/b968mTp16mTYuiw7nheoz5o1ixYtWrzy2rCwMNq3b0+zZs1Yv349NjZ57+QESZ2KLuw42tOb4Fk8yOSg04FcDpIO8hdB4d4Buas3MpVtrrRJlDMIgmDVbt++zfvvv8+DBw8MPYsWLVrQp08fOnfubNS9EhIS8Pb2plu3bowYMeKV127fvp2ePXsyZcoUevfubXL732W6B5FoNk4ErRrUKVlfqLIDhRKl7zjkJXK+EF8sbhEEwaq9PL8XERHBhQsX+PTTT426T2pqKh06dKBJkyYMHz48y+t0Oh2TJ0+mb9++bNmyJe+G3r1raNaNhZSnrw490H8/JRHN+rHo7obneNtE8AmCYNVent9bsmQJ3bt3x9Y2+8NqOp2Ob775hiJFijB79uws56SePn3KZ599xo4dOwgKCqJBgwZmt/9dJCXGoQmY+PrAe5k6FU3AJKSnsTnTsH+J4BMEwaq9WL+n0WhYvny50bV7I0eO5Pbt26xcuTLLAvXr16/zwQcfULRoUQ4fPmx2Ufy7TBu8DdRppj1Yk4Y2eKtlG/QSEXyCIFitmJgY4uLiqFmzJgA7duygQoUKVK9ePdv3mDNnDtu3b2fr1q3ky5cv02t2796Nl5cXAwcOZOHChUb1Jq2NpFGju7gfdBrTbqDToLt0AEmTatmGvUCs6hQEwWodOXIEb29v5HL9e/zFixcbNee2YcMGpk+fzokTJzItUJckiV9++YW5c+eyceNGGjVqZLG2v6t0104C5q+Z1IWfRFG9qdn3yYwIPkEQrNaL83vR0dGcOHGCv//+O1uPPXr0KP369WPv3r04Oztn+P6zZ8/o2bMnN27cICgoiDJlyli07e8qKTrU+Lm9l6lTkG5fghwKPjHUKQiC1Xpxfm/ZsmV06dIFe3v71z4uNDSUzz77jNWrV1OnTp0M34+MjKRBgwbkz5+fo0ePitB7gZT0xDL3SbbMfTIjgk8QBKt0//597t27h5ubGzqdjiVLlmRrUUt0dDQ+Pj7MnDmTDz/8MMP39+3bR4MGDfjuu+9YunRptjamzlOUFprfVOZcsb8Y6hQEwSodPXqURo0aoVAo2LdvH0WKFMnyzLvnEhISaNu2Ld9//326Q2FBP583c+ZMZs6cydq1azOc9JBXqdVqLl68yKlTpwgKCqK+NoKv3YphozCjXyWTIyvkaLlGvkQEnyAIVunw4cOG+b3sHD/0vEC9cePGjBw5Mt33kpKS6NWrF2FhYZw6deq1m1JbK0mSiIqKMoTcqVOnCAkJoXz58nh6evLBBx/gXfNLbIKXgNbEcgYAhRJ59WaWa/hLxJZlgiBYpZo1a7Js2TLKly9P5cqViYqKwsHBIdNrdTod3bp1IyUlhfXr16er1bt58ybt27enZs2aLFy4MMuSBmv06NEjTp8+zalTpwxhp1Qq8fT0NHzUq1ePQoUKpXucevVIpPvXTX/i4hWw+XKGma3Pmgg+QRCsTmxsLC4uLsTFxTF37lxCQkL4888/s7x+xIgRnDhxgv3796cLtkOHDtG1a1dGjBiBn59frp8ikJvS0tIICQkx9OROnTpFTEwM9erVw9PTEw8PDzw9PbO1kEd37R80e+aC2oRaPJUtylYDkFdpaMJPkT1iqFMQBKtz9OhRvLy8UCgULFq0iIULF2Z57dy5c9m6dSsnTpwwhJ4kScyZM4dp06axatWq157C8K6RJInIyMh0PbkLFy5QqVIlPDw8aNy4McOHD6d69epZ7lTzKrJKnsiuuSNFBIHGiCFPpQ2yCvWQVc7Zrd5E8OU4CUgF1P/+txKwA6z3naMgvGnP5/cCAwORJCnLwvKNGzfy66+/cvz4cYoVKwZAcnIyffv2JSQkhJMnT5p0dNHbJi4ujqCgIENvLigoiHz58hmGK319falXrx4FChSwyPPJZDKUrQei2TkLKeo8ZGcXFqUtMmc3lG1+yPGetRjqzDEa4B5wG0jjv6CTABVQFij57+eCIFhS7dq1WbBgAfPnz6dmzZoMGzYswzXHjh2jY8eO7N2711Crd/v2bTp06EClSpVYsmRJtmr+3japqamcP38+XW/uwYMH1K9f3zBc6eHhkSt7iUqShDZ4K7rTm7I+mkhlBwoV8vrtUdT/BJks56vsRPDliBjg2r+f67K45vn/3ApA3lwhJgg5IT4+nvLlyxMZGUmlSpUIDw+nePHi6a65fPkyzZo1Y+XKlbRs2RLQD49+/vnnDB48mGHDhr0T83mSJHHt2rV0qyxDQ0OpUqVKunk5V1dXk4YsLdZOnRbpxlm0wVuRHt/Vz/2pbJE5lEJRr51+eFOee+0TQ50WdwO4RdaB99zz799APxRaOScbJQh5xrFjx2jQoAEbNmygZcuWGUIvOjqatm3bMmPGDFq2bIkkSfzvf/9j4sSJ/PXXX7Rq1eoNtfz1Hj58mG7xyenTpylUqJAh4Lp06cL7779P/vz533RT05HJFchc3JG7uL/ppgAi+CzsLtkLvRfp0PcQ7dAPfwqCYI7n83uLFi1i6tSp6b6XkJCAj48Pffv2pXv37qSkpNC/f3+CgoIIDAzExSXnT//OruTkZM6dO5euNxcfH4+7uzuenp70798fDw8PSpYs+aab+s4RQ50WowOOA1oTHy8HvBDvRQTBPO+//z4DBgxg4sSJREZGGk5mSEtLw8fHh6pVq/Lbb78RExODr68v5cqVY9myZRZb2GEKnU5HWFhYut7c1atXqVatmqE35+npSZUqVQw/j2A68SprMQ8tcI97gNjsVhBM9fjxY8Oc17fffmsICZ1OR8+ePSlYsCBz584lMDCQzp07M2DAAH788cdcn8+7f/9+usUnp0+fpmjRooaA69atG3Xr1s1TxfK5SQSfxdzE9N4e6HuMt4HSiFIHQcgeKSkB3e1LkPIUZHLCQ8Np/IE769evJyQkxHDdqFGjiIyM5MCBAyxevJixY8eyfPlyfHx8cryNSUlJnD171hB0p06d4unTp3h4eODh4YGfnx8eHh44Oubc3pRCeiL4LEILPLPAfdL+/ci7pzcLwutIkoQUE4Y2eLO+RkyuBEn/ptM1Tc3a1g6cqONNaTv91+bNm8eWLVs4ePAgfn5+HDt2jBMnTlC5suUXlGm1Wq5evZpuXi48PJwaNWrg6enJJ598wuTJk6lcufI7sWrUWok5PotIBf7BuEUtmVEA9YB3r3ZIEHKDpEnTF0XfugDqNLI66VuHDLlSxQ07Z1qM/4uNAQEMGTKE4sWLs2LFCgoWLGiR9sTExKSblwsODqZ48eLpSgnq1KmDra14M/s2EcFnESL4BCGnSVo1mvXjkB5GZXsbrKQ0LbH2ZWg6bTPfffcdY8aMMXlxSGJiIsHBwel6c8nJyekWn7i7uxt2gBHeXiL4LEIHHLHAfWRAQyDnDmAUhHeVetccpOv/GLf3I/AsTcvdEnWp/s3/ZfsxWq2W0NDQdL25iIgI3Nzc0vXmKlasaNEhS0mSIC0Z1Mn6HU1s8uXKTiZ5jZjjswg5UBhIMOsukpQfmUyEniC8THryAOnaSf22V0ayt1FQ6WkYUloyMpuMqyQlSSI6OjrdKsvg4GCcnJwMPbk+ffrg5uaGjU3O/H5KifFoQ/agC9ml39ZLrgCdTn8uXc0PUdRpi8xB1OtZiujxWUwcWu0FTN0V6NmzFH766S9cXBrQo0ePN1pTJAhvG83RFejO7wCtxrQbqGxRNP4GhVsrnjx5wpkzZ9L15jQaTbqenLu7O0WKFLHsD5EJSZ2KZu9vSBGn9V/ILNjlSpDJkJWujvKjIcjsxDNmydoAACAASURBVGuDuUTwWUBMTAwjR45k1qyOODoWNukekqTgn38UzJzpz+HDh+nZsycDBw6kbFmxm4uQt0laNer5PfRDgGZ4mCqj9cqrREVFUadOnXRzc87Ozrm+ylJKfYZ63Rh4dC97p5UrlJDfAdXn05AVyPlQtmZi8NgMqamp/PLLL7i5uVG2bFkKFHDHlL9SSZIjk1WjQQMvNmzYwOnTp1Gr1dSuXZsvvviCoKAgyzdeEN4VTx6CBd6fv2cr8dfyZTx69Ijjx48za9YsunTpQvny5XM/9LQaNJsmQ3xM9kIP9L3dZ49QbxiHZOabgLxOBJ+Jdu7cSa1atTh+/Dj//PMPP//8M/nylQaqYsxfa0pKGtu2XQH+K16tUKEC/v7+3LhxA3d3dzp37kyjRo3YuHEjWq05RfKC8O6RUpPAEgs8FCpqV6uMSvXmjwLThR1Hir0JOiOHbnVaePIQ7bkdOdOwPCIPDnVK6AvOZegDyrh3etevX8fPz4/w8HBmz56dxc4P8cCVf58nq6BSAHIePy5F/fqtGTlyJL179870So1Gw+bNm5k1axb37t3jhx9+oGfPnhQqVMiotgvCu0SSJO7cucOtcyeoHRGAjVk7I0GaVkfZ8btQ5S9EqVKlKF68OMWLF8fR0RFHR0fD5y/+WahQoRzpDab9OQji75h+g3yFUH23OFeP8rEmeST4dEAs+pMTnpI+7BzRn4rw6hBJTExkypQpLFq0iBEjRjBo0KDXFKVKwKN/n/MR6Q+iLQw4A0UBGdevX8fb25vFixfz0UcfvbId//zzD/7+/uzfv59vvvmGgQMHUr58+Vc+RhDeZhqNhoiICK5cuZLu4+rVqxQoUADP2tVZ1dIBlZmdPrUko/OBZC5fvkx0dDTvvfcexYsXx8HBgXz58iGXy9FoNDx58oQHDx7w8OFDUlNTDde9GIqZBaWjo2O2glJ3PwLNurHZO5U8Kyo7lD5DkFesZ/o98rA8EHx3gev819PLjBz9sUA1gPQrpiRJYs2aNYwYMYJmzZrxyy+/mHBy8fPnltBXkGT8xTh16hQff/wxO3fuxN399WdW3bp1i3nz5rF06VJatGjB4MGDadCggZHtEoTc8+zZM8LCwjKEW2RkJKVLl6ZatWrpPuzs7Dhw4AAbN25kqocdDZzNWNAhkyOv1hRl6/6A/qSGyMhIwsLCuHr1KmFhYYbPtVotVatWxdXVlYoVK1KyZEmKFCmCnZ0dCQkJPHz40BCMz/98/nlaWlqWofj889opYZSMPoUsi11nsv0juTZG1XaQWffIq6w8+CLRb/yc3R1V5IAboP8FO3/+PAMHDuTZs2fMmzcPLy+vnGnmv7Zu3Urfvn05duxYts8Fe/r0KcuWLWPOnDk4OjoyZMgQfH19USpFiabwZsTFxWXovV25coUHDx5QqVKlDAFXpUoV7OzskCSJy5cvExAQQEBAADExMXz66ad07NiR+o4KVIfmYyc37eUqWaPjYpXOeH3y+Wt7ZLGxsYYgfDEYo6KiKF26NFWrVjUE4/PPS5YsiUwmIzk52RCEmQXkgwcP6FFeQzsX83dnkpVzQ9Ux+0X5wn+sOPhuow8+Y7cRk/P4sQujR//Mxo0bmTRpEt9++y0KUwv0jDR//nxmzZrFiRMnjNqtXavVsnXrVvz9/bl58yYDBw6kV69eODg45GBrhbxKkiRu375t6LW9GHBpaWnpgs3V1ZVq1apRoUKFDL9HkiQRHBzMxo0bCQgIIDk5GV9fX3x9ffHy8kKSJBYuXMikiRM47+dFEVsT5ttkch4rC9P4t+OUKlWKadOm4enpafRt1Gq1oZf4ciimpaUZQvDFUKxcuTJ2dnbp7qPZ8xu6y4eM/zle/rHK1kTVaYLZ93mdpLtxhC/ZSfz5CNIeJ6IqZE+RWuWp8q0PBcqVyPHnzwlWGnxpwElM2TtTkiQuXLjBokXnmDhxIkWLFrV4615n9OjRHDp0iAMHDpA/f36jH3/mzBn8/f3ZtWsX3bt3Z9CgQVSsWDEHWipYO7Vanen8W1hYGAULFszQe6tWrZqh95MVrVbLiRMnDD07Ozs7OnbsSMeOHalXr57hsQcOHMDPzw9HR0dq1qzJhf2b2fz1++S3MfJNqG1+VF/OQGtfjOXLlzN+/Hg8PDz4+eefcXV1NeevxyA+Pj7DsGlYWBiRkZGUKlUqXe+wZaFHlH0QbPbhY7KqjVD5DLZI+zPz4J/LhExeScyBsyCToUv5r+xCbqtCJpNR3KsGtX/qRqmmdXKsHTnBSoMvCv35eKZtGq3TgVzuwZvaLFqSJL7++muePHnCxo0bTe5t3rlzh99++43FixfTuHFjhgwZgpeXlzgORcggs/m3K1eucOPGjUzn31xdXY0aTUhLS+PQoUNs3LiRLVu24OTkZOjZVa9ePd2/yYiICIYNG0ZISAj/93//x9q1a7l27RoJCQkc+HMuVSO2gTqVrE5meE4ngdyuAMrPJiB3LG/4enJyMr/99hvTp0/nk08+Yfz48ZQpkzMHQGs0Gm7cuJEuFBUPI5jaoBD2xgb4i1R2KFsPQF45Z+b1r/xvC6dHLECbnPbaGkpFfltqj/4St1Fd35nXFisMPgk4ARi/p196pQDLvBs0RVpaGh999BGVK1fm999/N+sfVGJiIn/++SezZ8/GwcGBIUOG0KlTp7einknIXbGxsZnOvz18+JDKlStnCLfn82+mSEpKYs+ePQQEBLBjxw6qVauGr68vHTp0yHQE4unTp0yZMoXFixczdOhQvLy86NatG87Ozty+fZtdu3ZRrVo1dLG30B5ZhhR9BbVGjeqlX41ktRZbGxt2XblHfb+5lK2WeW/k8ePH/PrrryxYsICePXsyatSoXBnhkSSJlEV9UTyLNf0mtvao+ixFprD8XP7V+VsJGjYfbVL2V50q89vh9tOX1B7V1eLtyQlWGHxPgXOYdxo66FdfepvfHDM8efIEb29vvvjiC3788Uez76fT6dixYwezZs3i+vXrDBgwgO+++y5X9iQUcs+L828vf6jV6kyHJ8uXL2+ReeyEhAR27NjBxo0b2b9/P+7u7vj6+tK+ffssV0PrdDr+/PNPfvrpJ1q1asWUKVMM8+uenp7cvHmT3bt3U7p06fQ/59NYFvt9Tud65cmn0AEysCvAzutPuGZTlmRJycWLF1m/fv0r2xwTE8OkSZPYsGEDgwcPZtCgQdjb58xoT1JSEjNnziTu6Domt6mKjQmLddK0OnbfVxJaoKZhPrFKlSomTYu8LPZsODu9/dAmG19qochvS6sdUynZpLbZ7chpVhh88cAlzA8+gKYYW+BuadHR0TRs2JApU6bQrVs3i9333Llz+Pv7s337drp27cqgQYNy5ERqIeeo1WquX7+e6fxb4cKF0y0sye78mykePnzIli1bCAgI4Pjx4zRt2hRfX1/atWv32rPpTpw4waBBg7CxsWH27Nm4urrSq1cvwsPDKV26NImJiWzZsiXTYdXHjx9TtmxZYmNj09XUnjx5kp49exIcHEyNGjVYsmQJzZs3f+3Pce3aNcaOHcvRo0cZO3YsvXr1stioiE6nY8WKFYwZM4ZGjRoxdfJESgf+AY+i9buxZJMkk6FWFWC7fSMuXrthmEu8fv06jo6O6eYSn4di6dKls30G4aEuE4naeFQ/TmyCUh++T5u90016bG6ywuCLA0IxN/gkSeLWrYqUKVM211Z0ZiU0NJTmzZuzevVqWrRoYdF7x8TE8Pvvv7Nw4UIaNmzI4MGDadKkyTszVp8XPHv2LMPKyefzb2XLls10/q1wYdM2S8+u27dvs2nTJgICAjh//jytW7fG19cXHx+fbJ1ufuvWLUaOHMnx48f55Zdf+OKLL7h48SKdOnXCy8uLyMhIHB0dWblyZZZDrVu3bmXevHns27cv3dclSaJChQps3bqViIgIxo4dy/nz57Nd4hMcHMzo0aOJiIhg8uTJdO7c2eTDawEOHjzI0KFDyZcvH7NmzeKDDz7QtzMpAfWakfDsUfZOnZAr9EOcn0/NcESRVqvl5s2bGVabhoWF8eTJEypXrpyhBKNKlSrpToFJiUtgbZku6FJNnyZS2NngG/YnBcoWN/keucEKg+8JcB5zgy8pKZUqVXoRGxtLuXLlcHFxoWLFiri4uBg+r1ixYo4NibzsyJEjfPbZZ+zfvx83NzeL3z8pKYm//voLf39/8ufPz+DBg+nSpUuOnT8mZPTw4cN0hd2vmn+rVq1apkvlc9K1a9cMKzEjIiJo164dvr6+tGzZMtvtSEpK4tdff2XevHkMGDCAESNGYG9vz9KlSxk5ciTjxo1j8eLFeHt7M2fOnFe+6RwyZAjvvfceo0ePzvC9kSNHolAomDJlCq1ataJdu3b88MMPRv28Bw8e5Mcff0Sj0TB16lRatWpl1BvCK1euMGLECC5fvsy0adPo1KlThsdLKYlotk1HuheuDz8pkwV5MjkoVFDECVX70cgKGDcP+eTJk0xLMK5fv06xYsUMQVjrnhLbHVchzcSjnwC5jYpaI7vw/oQeJt8jN1hh8GnRL24xd6jzPaAWycnJREVFERERQUREBJGRkYY/b9y4QeHChdOF4Yt/lihRwqI9p7Vr1zJs2DACAwNz7LginU7H7t27mTVrFleuXKF///706dPntUNWQvbodLos59+0Wm2G2jdLzr8ZS5IkLl68aKixi42NpUOHDvj6+tKkSROjhgElSWLt2rWMGDGCBg0a8Ouvv+Ls7ExSUhL9+/fnn3/+Yfr06QwcOJDevXszatSo1/7u1K1bl//973+Z7lh07tw5OnbsaCjFaNKkCZcvXzaqNvZ5uwMCAvjpp59wcnJi6tSpr60BfPDgAePHj2f9+vWMGjWK/v37v2Z7Q5DibqM9ux3d1aMgk+nDTpJAp0VWyRNFvU+Ql8jephbZpdPpuHXrliEIUxYeofgV8w7TBnD29ab5hvHmNzAHWWHwgX6Lsju8brlz1uRAbeDVy7V1Oh0xMTHpwvDFgExKSsoQhs8/L1++vEm9qVmzZrF06VKOHz+e48XpFy5cYPbs2WzatInPP/8cPz8/qlatmqPPaS2yM//28oel3yiZQqfTERQUZOjZabVaOnbsiK+vLx988IFJQ37BwcEMGjSIpKQk5syZg7e3ftFYeHg4nTp1olatWvTq1YsvvviCqVOn0qPH63sL8fHxlC9fnri4uEwDWJIkXF1dWbFiBZ6engwePJhnz56xcOFCo9sP+rKE5cuXM2HCBDw8PJgyZUqGGsDk5GTmzJnDjBkz6NatG2PHjjX6DaOkSYXER/+dFm/vgEyVO736Q59PImrdYbPvU6pFXdrsm2F+g3KQlQZfMhCEqXV8+n07P8DchS0JCQncuHEj097inTt3KFWqVJbBmNVKS0mSGDx4MOfPn2fPnj2vfSdpCffu3eOPP/5g/vz5uLu7M3jwYJo3b27Ci7QEJKDfVScR/f8fOZAPKAMU400vJjJWYmJipvNvUVFRb2z+zVgajYZjx46xceNGNm3ahIODgyHsateubXIY37t3j9GjR7Nr1y4mT57MN998Y+i5rlu3jv79+zNp0iTKlClDjx49WL58+Ws3aX9u06ZNLFiwgN27d2d5zfjx40lISMDf35+EhARcXV3Zvn079eqZvrFzZjWATk5OrFmzhtGjR1O/fn2mTZv2Ti4UC+w3m7D528y+j3PHxjRf/3ZvpWalwQcQjn6DauO3LINa6E9OyDlqtZpbt25l6CU+/1ylUmUIw+d/Ojk50bVrV5RKJatWrTJr4t0YycnJrFq1Cn9/f5RKJYMHD+aLL77IZvjeA26gr6/MbBhaf0wTlEN/WsbbFYAvzr+9+BEbG0uVKlUynX8z5k2JpNMiRZ1HG7wVKe4WaNJAoUJWyBF53Y+QV2mITGm5+dbU1FT2799PQEAAW7dupXz58oYaO3N3M0lNTWX27NlMnz6dnj17MmbMGMMRWmlpaQwbNozt27ezfv16QkJCGD16NFu2bDFqG7EffviB0qVLM3LkyCyvuXLlCh9++CG3bt1CoVCwdOlSFi9ezPHjx83+nXleA/jbb79hb29PqVKl0vVm3wWSJBEZGcnRo0c5evQoidvO0SLeAVtMH1ZX2NlQe0w3ao/+0oIttTwrDj4JfVlDPMZtUl0FffH6myNJErGxsRl6ic9DMS4ujnLlyhEbG4uzszPdu3c3BGSFChVyfMGNJEns3bsXf39/QkJC6NevH3379s1i/kQCrpH9NyFy9JuE1yS3z0l+PueR2fE4L86/vfjh7Oxs1vybJEloQ3aj+2e9PuzUmZys/e9Ql7xOWxQNPje5aDkxMZFdu3YREBDA7t27cXNzM9TYOTs7m/wzPCdJElu2bGHYsGHUqFGDGTNmpOv53Lx5k86dO1OyZEmWLVvGH3/8weLFi9m9e7fRQ+hubm4sXrwYDw+PV15Xp04d5syZQ5MmTdDpdHh6evLDDz/QvXt3k37G565du8bIkSMJCgrC1dWV8+fPM2TIkBytATSXTqfjypUrhqA7evQoAE2aNKFx48Y0rOvO+WZj0KZk80T4TMhtVXx2YzX5S+b+Vo/GsOLgA/2L7nUgGn0PIvMX3pQUDXZ2tkB19Ita3m7JycncuHGDkJAQ/Pz8cHV1pWDBgkRERBAVFYWDg0OWvUVLzyOFhoYye/ZsNmzYwGeffYafnx/Vq1d/4YoI9POtxvS85eh73DXJiZ5fWlpalvNvRYoUyVD7llPzb5KkQ7v3d3ThJ7N3NpvSBlnxiig7jNHP/2TDo0eP2LZtGwEBARw6dIgGDRrg6+vLp59+SokSlttg+NKlS/j5+XH37l1mz55Ny5Yt031/x44d9OzZk2HDhhmKxE+cOMHOnTuNPubr4cOHVKpUibi4uNeWKEybNo2bN2/yxx9/APrjv3x9fbl69Wq2yi5eFhcXx8SJE1m1ahXDhw9n0KBB2NnZ5WgNoKk0Gg3nz5/n6NGjHDt2jGPHjuHg4EDjxo0NHxUqVEj37/ro19OIXH0ASWvCNJFMRtmPPPlw6xQL/hQ5w8qD77lUIIaMC14k0tLk+Pn5M2PGCvLnL5D5w99iUVFRNGrUiDlz5tCxY0fDgpvMFttERESQkpJiKMV4ORidnZ1NLl948OAB8+fP548//qB27doMGTKEli3dkckuYNpcq/m976dPnxrm316ch4uKiqJcuXKZzr/l5qn2moOL0IUeMu5AUoUKWakqKDv+X5anb9+7d4/NmzcTEBDAqVOnaNGiBb6+vnz88ccWXxAVFxfHuHHjWL9+PePGjaNv377pwkij0TBu3DhWrFjBmjVrcHd3p1u3bsTHx7Np0yaT5js3bNjAsmXL2LFjx2uvvXHjBh4eHsTExBiCqEePHhQvXpxffvkl28+ZmprKb7/9xrRp0+jcuTPjx4/PdIQjODiYUaNGERkZaZEaQGOkpqZy+vRpQ2/u5MmTlC1b1hBy3t7eGXa/edmj0Ci2efYzaruy5xT5bWl7YCaOntVM/RFyTR4Jvud0QAr6eSY5oALs+Pjjj+nQoQPffvvtG22dqc6dO0fr1q3ZtGnTa88MTEhIyDB8+vzP6OhoSpUqlWlphouLS7ZeNFNSUlizZg3+/v7MmdOLJk1qIpeb2lPKB3jyql6fJElZzr/FxcVZZP4tJ+huX0Kz+WfTTuFW2iBv1A1l3f8WgkRFRRkKyi9duoSPjw++vr60adMmR4be1Go1f/zxB5MnT6ZLly6MHz8+wwrGu3fv8sUXX6BUKlm9ejU2NjZ8+umnlCxZkhUrVpj8/6B///5UqFCBYcOGZev6Dz74gPHjx9OmTRtA/8agZs2anDhx4rVDrJIksX79en788Udq1qzJL7/8QrVqr39hP3DgAKNGjTK5BjA7nj17xsmTJw1Bd+bMGVxdXQ1B16hRI957z/gRrIg1BznRa4ZR25Yp8tviMaMvrn0/Mfr53oQ8FnyZ2717N6NGjeLs2bNvfDm5qfbu3ctXX33F4cOHTV6c8HzBTVa9RZVKlekKVBcXF0qXLp1urkuSUtDpAlEozPn7lAN1gUJZzr9duXIFSZIy7b2ZO/+Wk9SbJiNFnTP9BgWKEeHlR8C/YXf79m0+/fRTfH19ad68eY4G+549exg8eDBlypTB39+fGjVqZLjm8OHDdO3ald69ezNu3Dju3btHmzZtaN68Of7+/mb1gmrUqMGKFSuyvTpzzpw5nDt3juXLlxu+NnPmTPbv38/OnTuz/J0/efIkQ4cOJSUlhRkzZmRr27MXSZLExo0b+emnnyhdunS2agBf5dGjR5w4ccIQdBcvXqRu3bqGoGvYsKHFRiwi/z7I8V4z0KWqXznsKZPLkdup8JjVD9fvPrbIc+cGEXzoJ31dXV1ZunQpjRo1etPNMdnzOqOTJ09SsmTJ1z/ACM97VlnVLMbFxeHs7GwIw06dPGnUqDRKpenBp9PpOHw4nKFDFxIeHk6RIkUyXWBSvHjxd+oNi5QYj3ppP9CavjXUM7WWXgFhlPZoSceOHfHy8sr2llymCg8PZ+jQoVy9epWZM2fSrl27DH/vOp2OadOmMXfuXFasWEGrVq24cuUKbdu2pV+/fgwfPtys/1f379/H1dWV2NjYbL+piYmJoUaNGty9e9eww0xaWhpubm7MmDGDjz/+GE1KGncPniPlwSMe3L3P2u2b2HfjPD9OnUD37t3NCmqNRsOyZcuYMGECnp6emdYAZvWzHjt2zBB0EREReHp6GoLO09OTfPmyN9driseXo7g4fR031h5CppCjeZZi+J4ivy3oJMq196LWiM8pVqdSjrUjJ4jg+9ecOXM4efIkf//995tuilkmT57Mpk2bOHz4sEmT96ZKSkpKt8NN7dr5adbM/FqmmJgU7twpluvzbzlJe3Y72uMrzQo+CZBX9Ubl42e5hmUhISGBSZMmsXz5ckaOHMkPP/yQaY8yLi6Or776isePH7N27VrKlClDYGAgHTp0YPr06Xz11Vdmt2Xt2rWsWrWKrVu3GvW4Zs2aMWjQINq3b2/42p49e/jpu0HM/aQfEX/uBRmkpqSiUatR2dqgRI7Th+9Ta1gXSjR2M/vNVXJyMvPmzWP69Ol8+umnGc4BvHXrVroVl/fu3aNRo0aGoHv//fffyBaCaU+eEbn6II8uRpIa/wSbIgVxqO6MS9cW2BZ9N38nRfD9KyEhgfLlyxMaGmr0KrO3iSRJ9OnTh1u3brFt27Y3uLLsKvoSBnMVAkwvOH4baY78ie6scS/cmZE5VUPVZbIFWpQ5rVbL0qVLGTt2LB9//DFTpkzJciXoqVOn6Ny5M506dWLatGmoVCq2bt3Kt99+y4oVK2jbtq1F2tS3b1+qVq3K4MHGnTy+YMECDh06ZHhjK0kSF6au5vS4pchlMmTaLF4GZTKU9rY4urvSYstkVAXM72E9evSIX375hfnz5+Pp6YmDgwOnTp0iKSkp3YrLWrVqvbVD9e86EXwv6NevH46OjkyYMOFNN8UsGo2G9u3bU7x4cZYsWZLJO1UtkIZ+sY8CsMHyNXORwE0L3KcYYPlNud8kzcFF6EKy3nEku2QlKqHqmv2VicY4cuQIfn5+FChQgDlz5vD+++9nep0kScybN4/JkyezYMECOnToAMCiRYsYN24cW7ZseW2tnTGqVq3K2rVrqVMn88NlsxIbG0ulSpWIjo7G3t6eoGHzuTp/a7ZXLyrsbChUuTQfBc5DZW98+Ol0Oi5dupSuR6dQKLC3tycmJobevXszadKkdKclCDlHBN8LLl++TIsWLbh58+Y7fyrBs2fPaNq0KT4+Pv8G+YvbhcWhXykp+/frcsAJKI1+uzZLeAhcwbzNwuVAecD84uq3iebkWnT/rDP7PnfkxXjs3ZdatWpZrGcfFRXF8OHDCQoKYvr06Xz22WdZDvE9efKEXr16cf36ddavX4+LiwuSJBmGRffs2WPRrbtiYmKoVasWDx8+NGnOzcfHh+7du/N+UmFODpyLlGLcULPCzoaSTWvTaue0116rVqs5d+6cIeSOHz+Oo6Njuh7d800Drl27xpgxYzh27NhbUwNo7XJ3a4y3XPXq1alevTobNmx4000xm729PTt27GDVqlWsW/cncAoIAWLRh50OfSjpAA36QDyF/ixDU/c4fZGl9t18d4edsyJ3cjXsxmIqNXJO3knkq6++wsHBgYYNG+Ln58fq1au5du0axr6fTUxMZMyYMdSrVw83NzeuXr1K586dswy9CxcuUL9+fYoUKUJgYCAuLi5otVr69evH5s2bCQwMtPh+lYcPH6Zx48YmLzT5/PPPWbZkCYf7zzI69AC0KWncO3qBuPPXM3wvJSWFI0eOMGnSJFq2bEnRokXp3bs3N2/epHv37ly+fJmwsDAWLVpE9+7d0+2UU7lyZdauXcu2bdvYtGkT1apV4++//0ans8TvoZAZ0eN7yebNm/n1118JDAx8002xiFu3QnFwiKJAgXzZrKeTA/boywjMnV+4gX6409R/Yo7od2+xLpKkQ73oO/0BpKZSqFB9txiZXQGePn1KcHAwQUFBho/ExETc3d3x8PAwfGQ2P6fT6Vi9ejU//vgjTZs2Zdq0aekWXGTm+dl5/v7+dOvWDdAv3OjatStPnz4lICAgRxYi9e7dm1q1ahl9rh7oe6fjx4/n0Oy/6K+qg0JtWqjIFHIqdm1B3d/6ExgYaOjRnTt3jho1ahh6c15eXhQtatq2XebWAEqpSehuBP93wK2tPbLSrsjfs66RE3OI4HuJRqPBxcWFgIAAs3ZxfzukAKfR9+iMIUd/JJMb5vXa0tD3Ik052FIO1EcfwtYn4fBq5Gc2YmvKewuZHFkVr1eu6Lx37x6nT59OF4aFChVKF4Q6nY5Ro0ah0+mYM2dOpufavejFs/M2bNhgqN+Lj4/nk08+oVy5cixfvjzHpgkqVarEpk2bqFWrJLbk9AAAIABJREFUVrYfo9FoWLRoERMmTKBt27ZU33Kb9x6b95KnkUmMsgumRv06hqBr0KCBRefnTKkB1MXeRBe8DV3YCZDL/zvYVq7U/xo7lELh3gF55Q+QKfL2UKoIvkz88ssvXL16lWXLlr3pppjpMnDfxMda6pSKp8A5jJvrkwM1eBf2TTWWJEmsW7eOsSOGcHKgFwUUOmTG9ohVdqi6/oqs6Ku3n3r5ea9fv05QUBAHDx5k69atxMXFUbp0aT788EM8PT3x8PDIcr7wxbPzFixYYHiRv337Nm3atKF169bMmDHDIttzSanP0F4+jO7ifkhOAJ0WrdyGtYGhdJ+xBoXj63sukiSxY8cOhg8fjpOTEzNnzqROnToss2sFaeYdUi23t6XF9imUaVLXrPtkR3ZqACVJQvvPOnRnNmd9ivtzKjsoUBRVp4nICmR+9FleIIIvE7GxsVSuXJnw8HCjT2x+e6iBQMybrysCGLd6LnOJ6MNP95r2yNG/Na1JTh8L9Sbcv3+ffv36ceXKFZYtW4ZHlbKoV4/Un8iQ3V9DpQ3KT35E7lzb6OdPTk5m1qxZzJo1iz59+jB06FCioqLS9QqjoqKoXbt2up5hcHAwAwYMYNKkSfTp08cw7BYaGkrbtm0ZNGgQQ4cONbo9L5NSnqI98qe+xyKTZdjOTasDhY0NFCmNsmkP5GUy7hgDcP78eYYOHcrdu3eZPn06Pj4+yGQyJEliueJDs9upKmxP09U/Uaat6buwGOtVNYCaI8vQXdiX/e3v5AqwK4iq20xk9jl7mPXbSixuycR7771Hhw4dWLJkyZtuihnuWeAeCeiHS81VAGgAVEKrtSExMZm0NC36kJOj1cp4+DABnc4Z/QHA1hV6kiSxZs0a3NzcqFKlCmfPnsXT0xNZESdUX/4K9kVfv9hFZQs2+VH6jjM69CRJYsOGDVSvXp2zZ89y+vRpfv75Z4oVK0a9evX4/vvvWbZsGaGhody9e5eff/4ZJycnNmzYQN26denatSuVK1cmJiaGHTt2GHYUad68OVOnTrVM6D15gHrlMHRXj4E2LdMXcYUc/dFND2+g2TQZ7eVD6b4fHR1Njx49aNOmDZ06deLChQt89NFHhqCOj4//9ybmk9vm7qrvfPnyMWLECMOb8dq1azN8+HCenNpmXOgB6LSQ/BT1xglIr+odWjHR48vC2bNnad++PZGRkTm+FVTOOIN+mNEccqAyllxZuWXLFubOncz8+TOpXLkC+gU0dtSp04y5c+fRuHFjiz3X2+DevXt8//33hIeHs3z5ctzd3TNcI2k1SBFBaE9vQoq/o39HLkn6Xo8kQf5CKOp3QO7qne3jiJ47f/48fn5+PHr0iNmzZ9OsWbNsPe7Fs/OmTZvGtWvXDL3CwMBAkpOT8fLyol27dnh4eFCvXj2T57ik5KeoVw7VL8Yw5oVYaYPSZzBJJarz66+/8vvvv9OnTx/8/Py4c+cOFy9eNHycO3eOhIQEpqs9KYh581vK/La0O/0HDtXe3GKRmJgYJk2cwGinaEoUNHFfVpUdyo+HIS+f80O2bxsRfK/QsGFDhg8fbijKfbf8A2RyqKnRKmLJOroBAwawcOFCHj16lO7UgJ9//pno6Gh+//13iz3Xm/S8lzd48GB69erFuHHjsrVxtBR/Byk+Gik1CZnKDgo5IivhYvR2WQ8ePGDs2LFs3vz/7J13fE3nH8ffd+RmWFmCJGSYqVW199aW2kGMpCo1Su0Ro1RVf7VKlaIaSlRj+2mNInYVIbaEVERkkRiRxE1y1/P741YqMtzcG6E/9/16ndeVM57zPce953Oe5/mO/zJnzhw+/vjjXFlAhBAkHrnItW+28OBSFJr0TGRWFmhKK9h47zztAz5i0rSAHOdetWoVX3zxBd9//z1KpZLQ0FDOnj3L5cuX8fT0zDFEWqtWLYPi0dS7FyGizoKu8E5QaiGl5sJjlHGqgKurK1FRUcTExODi4oKNjQ0PHz4kKSkJjUZDiRIlGGrfkFqJMiQa43s6pau60Ov6+leeH1YXcxHVrvlItcYXjpVUrIWF9787YYcxmIWvAIKDg/nxxx85fPjwqzbFCEwXPp0ONJqKKBRFl4C2SpUqCCGIiorKsT4qKopmzZoRHx//L+1h/0NiYiKffPIJUVFR/PTTTzRo0KDYzq1SqVi+fDlff/01gwYNYtasWdjZ5XZiuBm0n7DP1qJKeYImPff3RGKtQCaTUWNEN96Z+xFSCzmzZ89m48aN7N+/n8qVK+c675UrV3LNF7799ts5xNDT0zOHYAjlY9SBw43OW/pEpWFtuJKrUleUSiURERFERkai0+mQSqXZVeZ9fHxwd3cnPTaJHdX80GYZdz55SWsaLxlJNf/ORh1flKi3f4G4c9m0RmQWWAxehqT0v9WXwTjMwlcAKpUKNzc3QkJC8iy98noTBqSa1IJSmcWECSvZv/8yb731Fl5eXtlB/l5eXoUuIpqUlIS7uzu9e/dmw4YNubY3aNCAefPm0aGD6Q4IrwIhBBs3bmTChAkMHz6czz77rNjq/gkh2Lt3LxMmTKBy5cosXrw4zwoAQgjOjP+evwL3ojEgXZfM2hLbWu7srZHFuWuX2Lt3L05OTgbZlJqamiu+UKlU5ogvbGnzEOtr+/Rzd0YSn5JBzYVHUavVVKtWjW7duuHt7U29evXy9DLd/14ACYfCIL/8nAUgL2mNT+JWo9KWFTWqlYMh08TpDIUN8vfHIvUsvpez1wGz8L2A2bNnc+/ePVauXPmqTSkUQsSi093EtBy3EjSahkRHJxAeHk5ERATh4eGEh4dz/fp1bG1tc4jh0+X5gqRP2bJlC5MnT2bixIl5BiEvWrQoO7vFv42EhARGjBjB7du3+emnn4o1BjQiIoIJEyYQHR3N4sWL6dw5/97I+ZlrubpkW6EqbGsl8LA0jLi1ndJ2pnkBJiQkcPbsWU6fPs2xY8dY36EMbnY2JrWZqZNwrUpP6nbsiY3Ni9tS3n3IrrpDUSY/QlqIOFWZtSXtd3yBy7v6eVqtVotarUalUqFWq3Mt+a03dtvz61e+rURhqq+OhTWy9sOQef1/za2/CLPwvYCntbyio6MNqkD+qnlaAX3t2tUcPDgbKytTvM/sgbw9CHU6HbGxsdlC+HSJiIhAoVDkEMKn4jh79mz27t3L5s2badasWa4279y5wzvvvENiYjwWFgJ94LsMsABez+FPIQQbNmxg0qRJfPLJJ8yYMaPY8rw+evSIL774go0bNzJjxgxGjRpV4Jza/bBI9rYeVyjRe4rM2pK3P/ejzhSfQh0nhCAhISGHo8nVq1e5fv06FSpU4MIndbEyqVgxZCHjsLQGNzVlDBYWxcNMGhx6iLVOhswA8VOhI9j6NuckydntAFhYWKBQKLCwsMi15Lfe2G3Pr+/9cA8WwpjkEM+gsEbe6VOkVZuY1s6/DLPwGUD//v1p3Lgx48a9/NpnxnLnzh1WrlzJmjVraNCgAaNHj+a999yQSEwJYK+DPpbPcIQQJCYm5hDC8PBwrl27pncnBz788ENq166dLYwVK1b8e94nk59/nk+/fi2xsJDxT9YYHfrcnxWBMhRNDlDTiY+PZ/jw4cTGxrJu3Trq1Sse7ziNRsPq1av54osv6NWrF3PmzDEo3vRo/y+5vfU4wsgckFZOtvgkbEWST5B6amoqV69ezSFyV65cQS6XU7t2bWrVqkXt2rXx8vLCwcGBlJQU6p76BpmJuWEztLD1vh3hqtKFEpaoC1dJ+m4ftTW2+ji/rOdERCZFKpdRsrorNed+iFOzmjnaedUlg1TrRsOjBNMasbBC7v0F0vL/rkKypmIWPgM4efIkgwcP5saNG0WSmaKoEEJw7Ngxli1bxtGjR/H19WXUqFHPJAfOQp+yrLAT+VL0vb1aFJXIxMXFUbNmTRwcHJgyZUqOXqJKlUlw8Cxat9bPoyoU+fXupIAlekE2bXjMFIQQrF+/nilTpjBy5EimT59ebL28Q4cOMW7cOBwdHfn222+pW9ewmL6sh6lsdu2HNtP4uTR5KWvabpqJU/t6REZG5hK45ORkqlevjoeHB05OTpQuXRqZTMbjx4+Jj4/PXpKSkihdujQKhYLLYxpha21i+iyFDfIPJiJ1K1yyBZ1Oh5ubG79t2obs1G3+Wvs7D2ISkAKlnBxwebcBNcf2fqVhCwWhubgP3YkNhYvhe55Sjlj4r3rlHqrFzes5fvSa0axZM0qVKsX+/fuLrKCmKSiVSjZu3MiyZctQq9V8+umnrFu3Lo+K65bok02fx/B8mVKgFPqUYUX3Yzhy5AiVK1emTp06jBgx4pktWjSas0gkGQbMR+oQIgM4h0TyDvrA+OIlLi6OYcOGkZCQwIEDBwpdF85YoqKimDRpEpcuXWLRokX07NmzUA+rmJ1/IDExeFuTlsHyPuNYoj6PnZ0dtra2WFlZIYRAJpMhlUq5ceMG6enpuLi4ZC+VKlWidGl9b+z+/ftIpVJ0Oh0PHz7kanIGzSpaYFD+9PzQqpGU9Sj0YVKplH79+rH99918+eWX1J7Uj/Hjx1OxYkX8J0wwwaDiQebVGt3xIOMbkFsird/tjRM9MAufQUgkEj799FOWL1/+SoUvOjqaFStW8NNPP9GsWTO++eYbOnTo8IIvbgmgIXAFUKKvlJC7k5+VpUYmkyGXO6MPWi/anu2RI0ewsbF5zrVfAJeRyw3PDiORgE6n5vHjowQEbKN8ebfsOcRq1aphZVVU9QRzIoTgp59+IiAggNGjRzNt2rRiqZmWlpbGV199RWBgIBMnTiQ4ONioa3wUnYD6SabJrzJltHLat29PpUqVcoibi4sLrq6ulClThvv373P8+PHs+dyEhAQsLS3JysqiZs2a+Pj40KZNG5o0aUKplNuods3TZ2sxAq1OcNeyLBUVJYx6mPn4+NC/f3/mzJmDRCJBrVb/a2pxSixtkNZoie76cePCQSQge6tNkdv1b8A81GkgGRkZVKpUiVOnTlGlSvGNhwshOHToEMuWLcsech05ciSenp5GtJaGvu5eEv/05gQg584dQY8eIzlz5nyRP9CFELi7u1OyZEkCAwOfqQJwH30i7cInDRYCoqOVrFt3JnvINDo6mooVK+byNK1Ro0aOYPnCEhsby7Bhw7h37x4//fSTwcOLpqDT6QgKCmL69Ol06tQpO42YEAJtRhbaLDUWpW2QGjDPdPnyZb5rM4TmKYULP8kL27fc6Hl1bY51iYmJHD16lF27dnH8+HHu37+f3QNs0KABXbp0oXXr1tSrVy+HqCQmJjLni9lMdY6ngpHZRzQSOROOJBNy7Q4BAQH4+fkVKoRECEFDdy9mNe+HuHGXhJu3UZSwoWLd6tT4pBuuXRobdI9fFUKVgTo4AFLu6lORGYpcgbzrlDcyawuYha9QBAQEoFarWbx48Us/V3p6OkFBQSxfvhypVMqYMWMYOHCgSQ/wf3hafFaH3mtSDkho164dQ4YMya6xVlTcunWLZs2akZqayv37959xOT+PPh+osUiB5jwduFCr1dy8eTOXl2lkZCROTk65wi5eFIsohGDt2rVMnTqVsWPHEhAQUCy9vJMnTzJu3DjkcjlLly6lUaNGPLxyi2tLthG96Qg6jRaJVIJOrcXWqxK1A/rj3qc18jw8eJ/WzlvUcQhi+wWE2rTKBE7Na1Lnl8mEhISwY8cOTp06RWqqPl60dOnSNG3alK5du9K6dWuqVauW52hEWloaCxcu5Pvvv+ejjz5ilk97rE4ZMVcls0BSrjLyvnP5448/+Prrr7l06RITJ05k2LBhL0yhdu/kVcKmBZJ46ioIgfQ5Hxt5KWtkCgveGteb2pP7IVO8nqV8hPIx6m2f68XPkJ6f3BLZe6ORVS24DNX/M2bhKwS3b9+mfv36xMTEULKkNXrnES36B68C0wu3ws2bN/n+++8JCgqidevWjB49mjZt2hTLOPyBAweYMGECly9fLlInnjVr1rB161YSEhK4fPlppokMIBTTqkdIgcpAwYVTtVot0dHRuWIRIyIiKFOmTJ6hF0qlkqFDh3L//n3WrVtXqBpwxhIbG0tAQAAnTpxg3rx5+sKuUQkc6TeHxzdi0WWpEdrc90teUh9MXe9zP2pO6INEIslROy84OJg7h86RNG0TMo3xP3eNFELkCezURaHT6XB2dqZVq1Z0796dli1bUr58+QKPV6lUrF69mrlz59KpUye+/PLL7Erkmj+D0YX9Zrj4ySyglIO+PJPlPy+DFy5cYN68eRw+fJhPP/2U0aNH51kQ9q/1+zk1cinaDMOC+B3erkLHfV+jKP161ocU6ix9aaLL+/XDIernpg9k+pdbiXMNZC0GvXFenM9jFr5CMnHiMD75pAtVqtiS0/lDAOXQu9wX7seh0+k4cOAAy5YtIzQ0FH9/fz755JPsh0JxIYSgfv36fPHFF3Tt2rXI2h04cCBSqRS5XP5MjcN44CamCR/owxveMerIp7GIz4rhtWvXuHTpEpmZmbi5udGxY8ccoRfly5fP8RKi1mmIeHSHsPs3SFNnoBVaLKQWONs40MipBi42ji98aVEqlSxcuJDvvvuOUaNGERAQQIkSJbgfFsnv7SehTleC7sU/U7mNFZ4D2uE48QN69+6No6MjLi4u/P7773i4e+AXXppSWca/QGkkgmu+VflgUB+aNGlicFJqIQRbt25l+vTpVKlShfnz5+c5XKwJ+w3dyY2A0NeVyw8LKyQOFZH3nInEKu/f2o0bN1iwYAE7d+7E39+fCRMmUKFCBQCitx7jxOD5BoneU6SWFjjUq8r7Rxe/tj0/AKFVo7t5Bt3lA4jsCuw2SCvVRfb2+29carL8MAufwWQBl9Fq0wBdATE8UqA0+lCAgn8gqamprFu3juXLl1OiRAlGjx5N//79sbZ+demQtm7dypIlSzh58mSR9DKFEDg7O9O6dWtatmzJqFGj/t5yG4g2uX19WEPR1EWLiYlh6NChPHr0iAULFqDT6XLFI6rVary8vKhZuxZVujQE1xJIpVI0eQi4hVROCbkVbSu8TQ27Srm2CyHYvHkzAQEBNGnShAULFmS/7KRFJ/LrO8NRPX5SqGvQyaX8rothjyKOBg0a0LhxYyIjIzl58iST6nfD+Y9EtBlGOJJIJLi+34iOu/9TqMOOHj3KlClT0Gq1LFiwgPbt2xe4v3ichPbiXnRXDv5TneJppQqdFomLF7IGPZBUqo1E8uJRidjYWL755huCgoLo27cvYwcP42z7zwolek+RWVtSc4I39b8cUuhjzbxemIXPIDLQ577UkJdHZG4k6EMJ6qMfAs3J9evXWb58Ob/88gsdOnRgzJgxNG/e/LVwK9ZqtdSoUYM1a9YUSYmg69ev895772Fvb8+KFSto0uRphojbFIXwRUbG4+29AAcHB+zt7Q36fN5rTwjB6tWr+eyzz5gwYQKTJ0/ON1F2cnIylyOuESaLRW0JUvmLh7flEhlNynnRsnyd7HVhYWGMHTsWpVLJ0qVLadmyZY5jDrw/lYSDYUYFm+vkEqzme/PTnm3cuHGD8ePHM3ToUKwkcnZ4DUaZ8MDwwrd/I7OxpMsf3+HwtmFDZFeuXGHq1KlERETw1Vdf0a9fv0INnwuNGhF7BaFMAa0GiVVJJOWrGt1jSU5O5rvvvuPmNztoryqPzMiBBosyJeh/b/tr3esz82LMwvdC1Ojnogr7lixB3xupD8jQarXs27ePZcuWcfHiRYYOHcqIESOyqyi/Tvz444/s3LmTvXv3mtzWihUrOHXqFNu3b+fBgwfP9GaLZqgzM9OK69etefjwIQ8ePMj38+m/Hz58iKWlZbYQ2tjYcPPmTYQQ9OjRgxo1auQpmPb29sjlcjQ6LUF/HeB+xmO0hbBdp9IgIlMon2nDvn37OHnyJF999RWDBw/ONXrwJD6Z7VX9jA421yC4bJ9BoyUj6d27N3/88Qe//PILJ06cQHn7HjPEO1ghQ2ZgyIrM2pLWG6fj1qPFC/eNjY1l1qxZ7N27l+nTpzNixIhiS9T9InQaLb849UKdkm50G/KS1rQInIRH3zZFZ5iZYscsfC8kCn0IgDG3SYpS6cKqVb/y/fffY29vz+jRo+nbt+9LizcrCrKysvD09GTv3r0mu+57e3tTu3ZtduzYwaVLl57ZkgmcwXTnliqAi8FHCCFIS0sjOTmZwMBAVqxYQZcuXWjSpAmPHz/OVzxTUlIoWbIknT7tS+1uzZEb8cavVWn4fsBnZD1MJyMjA09Pz1yhF9WqVePa3F+4+s0WdEaWzgHQyOArp0hi7+pTWrm6utK0aVP69u1LS696HOk8g6wHqXmWJHqKzMYSiURCu+1f4NKp4Oz9jx494uuvv2bNmjWMGDGCKVOmFLp6x8sm/sA5jvT5AnWa0qR2nJrXosuJpUVklZlXgTmAvUB0QALGiZ7++Lt3QwkLC2Pjxo00btz4tRjOfBGWlpaMHz+eefPmERwcbHQ7Op2Oo0ePUr9+/Txq0lmhd0x5ZIqpQMGehM8jkUh48OABQ4cO5cmTJ5w6dYq33nrrhcfpdDoepjxi3Z0Q1EbEHepPDt3HDuJJaCylS5dGCEFGRgbXrl3jyJEj2Sm9vhaNsdOYNpSm1enoU7c1XTZ+TMuWLXP1Kr0jg7iz6ySX528iJTwGqVyGTqtDKpMgBChK21BzYl+qDn4XS7vnMwL9Q2ZmJsuXL2f+/Pn06NGDy5cv4+Ji+ItIUSCEQAiBTqfLXrRabY6/dTodd69FodOYmNQZeBKbVARWm3mVmHt8BZIEXMeYAOun6HQSpNJ66B/y/x7S0tLw8PDgzJkzuYqOGsqlS5fo06cPrVq14p133mHkyJHP7fEQfUYZY3p9EvSil7vmXH7odDpWrVrFrFmzCAgIYPz48YUqehv+KIZ9sWdQGVEp/CkSHZS7piblft5DssnJyUyN98TGxHdSnULGg86VSavpWKAYaLVaFA8ysXqQiTRLg0YGmTYyHtnL0BUgJlqtlri4OG7cuEGpUqWoXLky1tbWLzzX8+uKYl+dTodEIkEqleZYngbRP11aZpWls7ICFiZmJbJysqX/3e0mtWHm1WLu8RVIMqaIHoBUKoAH/NuEr1SpUnzyyScsXLiQVatWGdXG4cOHadeuHWfOnGHo0KF57GH39/KIwoufHDA8e82tW7fw9/cnMzOTP/74I88irS/iwv2/TBI9AAu5nAbvN8E+y5Lk5GSSkpKIiori6NGjREVFce/ePaSSysYPMmQjkCBBoVAUKAb5rSto/aVLlwgKCsLS0lKf47J2bYPbLcy5CrOvISMpN4MOcOrTpWjSDU+RlxcWJV99EVozpmEWvgIxPpP9y2mneBkzZgzVq1fn888/z46BKgxHjhyhb9++BAUFUadOnTz2kKBPhn0JIdKQSAwRPwn6r2098vKYfR6dTseKFSuYPXs206ZNY9y4cUaXk0lTmzY3BKDMUPLJuE+JPHAOhUKBUqkkPT0dNzc3GjVqxLhx47CeE4J4WLgwhuextLLCx9+Pil2Krs5aWFgYU6ZMIS4ujm+++abQibJfNfZvV0YYEA9ZEBKpFMdGhX9pMvN68frU2DHz2lG2bFkGDhzI0qWFn8jXaDQcP34cR0dHqlWrVkBsogx4mwsXElCpNBT8lZSir8jQAEOSBERFRdG2bVuCg4M5efIkEydONKmGmlaYGmwPcgsLqlarikqlolWrVgQGBpKSksJff/3Fpk2bGDNmDM6dGyFM/GXqVBqcmtcy2V7Q95b79+9P165d6dOnD1evXqVXr17/KtEDsK9TmdKezia1IbNWUGtCnyKyyMyrwix8BVJUWdpfD3duY5g0aRI//vgjKSkphTru/PnzVKxYkejo6DwcW3ISHR1Dp06juXevMuCGPvBfgv7r+fSzPPrQkAboHWPyR6fT8d1339G4cWN69OjB8ePHqV69eqHszwtLmemxWzqNltbNW3L37l1++eUXevXqhY2NDVlZWWzbto3333+fMb+uQJggKhKZFLdeLbG0Na1sU3JyMmPHjqVhw4Z4eXkRGRnJiBEjiiVf6cuidoBPdoo3YyjhWhbHBqZ/l8y8WszCVyDlMDX/ZlaWhjNnbqHVmjZX+Kpwc3OjS5curFy5slDHHTlyhHbt2nHu3Dnq169f4L7jxo1jwoQJVKzoCbijTzzdAn1GluZAK8ALQ+rv3bx5kzZt2rBlyxb+/PNPxo8fX2SVsl1typo896awtOT9Zu2zw1muXLmSXQNuxYoVDBo0iAuJUTjVN/7hKrW0oOYEb6OPVyqVfPXVV3h5eaHVaomIiGDWrFkGpyh7nXH3boWlXUmMKQAos7ak/ry85qrN/NswC1+BOGDqLUpJyeSTTybj7u7O1KlTCQ8PLxrTipGAgACWLl1KRkb+MV/Pc/jwYdq2bUtYWFiBPb7du3dz/fp1Jk6c+Mzap/N4VvzT+ysYrVbLt99+S5MmTejduzfHjh2jWrVqBttbELdu3WLWrFlM8h6GRmV8bB2AnWVJrFRSVq1aRaNGjejcuTMlS5bk9OnTHD58mIEDB2JtbU2TZaOR2RR+pEBmbUnFD5rg+E7hr12j0fDjjz9StWpVLl26xKlTp1i+fDlOTk6Fbut1RWap4P0ji1GUKVEo8ZPZWFI7wAe37s1fonVmiguz8BWIBH3mf+Nuk0qlZfbstbRt25ZNmzYhhKBjx440aNCApUuXkpT074gHqlmzJo0bN34mwXTBqFQq/vzzTxo1akRkZGS+lQ0yMjIYM2YMy5cvNym7R2RkJK1bt2b79u2cOnWKsWPHmtzLy8jIYOPGjbRv357GjRtz/PhxkmMSUd4r3JDvs0h1Ei5sP46bmxuHDx9mzpw53L59my+//DJXfcWyDWvQ5pfPkFkbfl9k1pY4NqxOq6BphbJLCMGuXbuoU6cOGzduZOfOnWzZsoWqVauRJXC7AAAgAElEQVQWqp1/C6U8nekauhIbZwfkJQoeNpfIZcisFdT/cgj1ZvkVk4VmXjZm4Xshrujn6Ao7NCJBobDj889XkJaWRs+ePXF2dubmzZvMmzePsLAwqlWrxgcffMCWLVvIzDTNxfplM3XqVBYuXIjGgADg0NBQqlevTmxsLNWrV883S828efNo0KABHTt2NMomrVbL4sWLadasGX379uXYsWMmPayFEJw7d46RI0fi6urK+vXrqVSpEhYWFpQpU4a9e/fyQZXmeZYGeqGtGi2Pkx9StZT+O7Blyxbee++9AgW6UrdmdNo3D4V9KeSl8p+XkioskFkp8OjXhvcOLixUHslTp07RqlUrPvvsMxYuXMiRI0do1KhRoa7t30jpys70vr6ext+OolQVF+QlrJBZKUAqQWoh19fis7Gk6pD36Rb2AzXHGz90bOb1wxzAbhBZ6JNUqzBskkcKWKMvl6OPGLl69SqTJk3i1q1bLFiwgO7du/PkyRN27txJUFAQYWFheHt74+vrS4sWLV5Lj7nWrVszfPhwBgwYUOB+c+bMIS0tDXd3dy5evMiPP/6Ya5+bN2/SpEkTLl68aFS+0hs3bvDRRx9hYWHB2rVrjQ6yB7h//z4bN25k7dq1pKWl4evrC+hzljZs2JBBgwYRHh7O1q1befToEYMChmPfsio6iWE/HaHVoZBaMPStLpSxLPw8mU6t4c6uk5ybG8TDy7ewKV0SJBKEVodEKqH68K54jepOyUrlDG7zxo0bTJs2jbNnzzJnzhz8/PyKbC7034YQgvuh10k+ewNVSjpyawU2zo5U7NYUixLmmL3/R8zCZzBq4BpZWUnI5XJksryE6WkH2hF9RpHcD5L9+/czceJEHBwcWLx4cbbjR1xcHL/88gtBQUEolUp8fX3x9fWlSpXXp2Dk77//zpQpU7h06VKBwtymTRsCAgLYunUrjRo1YsSIETm2CyHo0qULbdu2ZfLkyYWyQavVsmTJEubNm8fs2bMZOXKkUUVztVotBw8eZO3atRw4cICuXbsyYMAArl27xjfffMNbb72Fp6cnf/zxB0+ePMHb2xtvb2+aNGmCVCol4tEddt85BYBG5O+4ZCGRUdLCmgFVOlBaYZPvfoawbds2tqxex8r5S9BmqVHYlqSUZ4VC9fASExP54osv2L59O5MmTWLMmDGvtAyWGTOvAvNQp8FYcPasmvbtpwHO6EXtWZd7OVAJaII+KDvvt+d3332XixcvMnDgQD744AP8/PyIjY3F1dWVKVOmcOXKFbZv305qairNmzenWbNmrFy5kocPHxbLVRbEu+++i1QqLbBqQ0ZGBufOnaNFixb5Orb897//JSYmhnHjxhXq/NevX6dFixbs2bOH0NBQPv3000KL3q1bt5g5cybu7u7MnDmTdu3ace3aNerWrYuvry+BgYEoFApu3bqFra0t69atIyYmJntI9en5vOwqMfKt7rxTpjIStQ51hgqVMhOh0SFDijpLhZOlLV3dmjHU6wOTRQ/g9OnT1G3dFId6VXFq8ha2NSoZLHppaWnMmjWLWrVqUaJECa5fv05AQIBZ9My8mQgzBtO9e3fx3Xff/f2XTgihFkJk/P2pK3R7qampYsaMGcLe3l7MmDFDpKam5tiuUqnEnj17RL9+/USZMmVEr169xM6dO0VWVpaJV2I8wcHBonnz5vluDwkJEU2bNhVKpVJYW1uLzMzMHNvT09NFpUqVxJEjRww+p0ajEQsWLBAODg7i+++/F1qttlA2P3nyRGzYsEG0bdtWODo6inHjxolLly6JtLQ0MXr0aGFjYyNKlCghKlasKKZOnSrOnTsndLr8/z91Op04cuSI8PX1FWXKlBHefbzF5pBd4lLyTXE+OVJce3hbdPXpKX799ddC2fkiWrRoIUJCQgp1TFZWlli2bJkoV66cGDRokIiOji5Sm8yY+TdiFj4DuXLliihfvrxQKpVF3vadO3eEr6+vKF++vPjhhx+EWq3OtU9KSooIDAwUrVq1Eo6OjmLUqFHi9OnTBT6gXwZqtVpUrlxZnDhxIs/t06dPFzNmzBCnTp0S77zzTq7t06ZNEwMGDDD4fOHh4aJx48aibdu24tatWwYfp9PpRGhoqBgxYoSws7MT7733nti6datQKpUiJCREtGjRQkilUlGqVCkxdOhQceHChRfey7i4OPHVV1+JypUri1q1aoklS5aIpKSkPPf9+uuvxZgxYwy290WoVCpRokQJ8fjxY4P21+l0YvPmzaJy5cqiU6dO4sKFC0Vmixkz/3bMwmcgAwYMEPPmzXup5zh79qxo1aqVqFWrlvj999/z3S86Olp8+eWXomrVqqJatWriyy+/LNY3+VWrVokuXbrkua1p06bi0KFDYtmyZWLo0KE5tkVERAhHR0eRkJDwwnOo1Woxb9484ejoKFauXGlwLy85OVksWbJE1K5dW3h6eoq5c+eK27dvixMnTogRI0aI0qVLC5lMJmrWrCl27tz5QrHLysoS27ZtE507dxZ2dnZi+PDhIjQ09IXHnTt3Tnh5eRlksyGEhYWJmjVrGrTvkSNHRMOGDUW9evXEwYMHi8wGM2b+XzALnwH89ddfwsHBweC3bVPQ6XRi586domrVquLdd98VV65cKXDf06dPi5EjRwoHBwfRunVrERgYKFJSUl6qjRkZGaJChQri0qVLOdanpqaKEiVKCKVSKQYPHix++OGHHLZ26NBBLFmy5IXtX716VTRs2FC0b9/eIEHXaDRi7969wtvbW5QpU0YMGjRIhISEiEOHDolRo0aJ8uXLi/LlywsbGxvRuXNnce3aNYNsGD9+vChbtqxo06aNCAoKEk+ePHnhcc/aZG9vL+Li4gw+piC+//574e/vX+A+ly9fFp07dxbu7u5i48aNhR4SNmPmTcEsfAbw8ccfi5kzZxbrObOyssS3334rypYtK4YOHSoSExNfuP/OnTtFz549RZkyZYSPj4/Ys2dPnsOmRcH8+fNzDVnu2bNHtGnTRgghRK1atcS5c+eyt23evFnUrl27QHvUarX4z3/+IxwcHMSqVate2Ku6efOmmDFjhnBxcRENGzYU33//vdi1a5cYMWKEKFeunKhTp47o2LGjsLOzEwMGDBDh4eEFtpeSkiJ++OEH0ahRI+Hs7CymT58u/vrrrxfdinzx9vYW69evN/r4Z/Hz8xOrV6/Oc9udO3fE4MGDhZOTk1iyZEmueVUzZszkxCx8L+DOnTvCzs5O3L9//5Wc/+HDh2LChAnCwcFBzJ0716A5xvv374sVK1aIJk2aiHLlyolx48aJ8+fPF+l84OPHj4WDg4OIiorKXjdp0iQxZ84c8eTJkxyOLampqcLV1TXfeUEh9HOoDRo0EB06dBC3b9/Od78nT56IoKAg0aZNG1G2bFkxevRosXLlSvHxxx8LR0dH0bBhQ/HFF1+IcePGCUdHRzFo0CBx/fr1fNvT6XTi6NGjws/PT5QpU0b07t1b7NmzR2g0GiPuSk5++OEH4evra3I7QghRrVo1cfny5RzrHj58KKZMmSLs7e3FtGnTXnpP34yZ/xfMwvcCxowZIyZOnPiqzRA3b94U3t7eomLFiiIoKMjgYazIyEgxc+ZM4e7uLmrVqiXmz59fZMNv06ZNE5988okQQoj7GSmis18vsePP/WL36RDRokPr7P0mTZokPvzwwzzbUKlUYu7cucLR0VGsXr06T3HW6XTizJkzYvjw4cLOzk68++67Yvr06eLDDz8U9vb2okmTJmLRokXi0qVLYvbs2cLR0VH4+fmJGzdu5Gv7U0eVKlWqiJo1a4rFixfn66hiLFFRUaJChQomv3A8ePBAlCpVKluMMzIyxKJFi0TZsmWFv79/kf1/mjHzpmAWvgK4d++esLOzE/Hx8a/alGxOnDghGjVqJOrXry+OHj1q8HFarVYcP35cfPzxx8LOzk507NhRBAUFibS0NKNtiU9MEM28O4qVV/4rFlwMFrOO/ygWXdos/hO6Qcw9GyR+jjwoDl3+UziWLSvu3r2b6/jLly+Ld955R3Tq1EnExMTk2p6UlCQWL14satWqJdzd3cWgQYNE7969hZ2dnWjRooX49ttvxZ07d8SjR4/E559/LhwcHMSHH34oIiMj87Q3KytLbN++PdtRZdiwYeLMmTMv1TPWw8NDXL161aQ29u3bJ9q2bSu0Wq0ICgoSbm5uomvXrgbNVZoxYyY35swtAOjQpyPTog9EtwCkTJ8+nUePHhW6JM/LRqfTsXnzZqZNm0a9evVYsGBBoXJUZmRk8NtvvxEUFMTJkyfp2rUrfn5+tG3b1uC0VYnKB2yOOsKTDCVShTzf/TSZKiyEjJH1e1FaoS8eq1armT9/PkuXLmXevHkMGTIkOxOMVqtl//79rF27loMHD1KvXj3kcjlhYWHUrVsXb29vevXqhbOzM48ePWLJkiWsWLGCbt26MX369Dwz3YSHh7NmzRp+/vlnatSogb+/P71796ZEiRcXszWVkRNG49GiNjUb1EWlU2Mts8S5hCNetpWQSw27159//jnXr1/nxo0bWFtbs2DBAlq2bPmSLTdj5v+XN1z4ngCxwL2//5bwNBdnZqY9rVr1Z/Pm3/Dw8HhF9hVMZmYmS5cuZeHChQwcOJBZs2bh4OBQqDaSkpIIDg5mw4YN3L17l4EDB+Ln50fNmjXzPSY2PYnNt46i1r04YTWABAlWMgWDq73LnRu3GDx4MOXKlWP16tVUrFgR0Ofu/Omnn/jpp58oUaIEtra2REZG0qBBA/r06UOPHj0oX748AA8fPmTJkiWsXLmS7t27M3369Fy5OlNTU9m8eTNr1qwhNjaWwYMHM3jw4GKrOBCTdo8/k64R8/guGq0WmcU/IqeQyhFAXfvKNHKqQRlF/gIcFhbGu+++i0KhYNmyZf/KyudmzLxuvKHCpwauAGnoe3u50Wp1aLU6FIqy6FOQvb5Vp5OSkpg9ezZbt25l2rRpjBo1yqgyP9euXWPDhg38/PPPlCtXDl9fX/r370+5cv8kP36UlcbaG/tQGSh6T5EAOqWa5T4z+OqLLxk8eDAZGRls27aN1atXc/nyZcqVK8fdu3dp1qwZ3t7e9OjRg7Jly2a38eDBAxYvXsyqVavo1asX06dPz/FSIoTgxIkTrF27ll27dtGuXTv8/f3p1KkTcnn+vdKiRAjBkYSLhN2PLDCHJ4AMKTKplD6ebahUMmfNu+joaGbMmMHhw4dJS0vj+vXr2S8JZsyYMY03UPgKW2lBgr4sUX1A8RLtMp2IiAgmT55MREQE8+fPp3fv3kb1DrRaLUePHiUoKIhdu3bRvHlz/Pz86NatGyFJF7n2KNqoQuRalZrGttUp/UjCqlWr2LJlC6VKlSItLY2WLVvSr18/unXrlqvXev/+fRYvXswPP/xA7969mT59Ou7u7tnbExISWL9+PWvXrkWhUODv78+gQYNeSQHVg3FhXHpwE/ULRO9Z5BIZPpXbUrGkE/fv32fu3Lls2LCBMWPG0K1bN3r27Mnt27dfntFmzLxhvGHCpwHOAZkYJnpPkQA26MXv9S/dcujQISZOnEjJkiVZvHixSfXVni2ddDniKqN3zkcqN/4epCU9YmnPKWi1Wlq1asWgQYPo2rUrdnZ2ufZNTk7mm2++4ccff6RPnz5MmzYNNzc3QF/sds+ePaxZs4aTJ0/Sp08f/P39adSo0SsbCgx/FMPe2NOodYaL3lMUUjlZ+2+xeP4i+vXrx6xZsyhXrhxBQUHs3buXTZs2vQSLzZh5M3nDhC8GuE1+w5sFIwUqoy9M+/qj1WpZv349M2fOpHXr1nz99dfZomEsIX+d4ezjv0BufFEPTaaKColy+rbvRunSpfPcJzk5mUWLFhEYGEi/fv2YOnUqlSpVAvSOKmvXrmXDhg3UqFGDIUOG4O3tXSyOKgUhhOCHiN08UqUZdbw6U8W9k5GM6/ZRjnnIUaNGUaVKFcaPH19Uppox88bzBpUlEugdWYwRPf4+7g6F6ym+OmQyGUOGDOHGjRtUq1aNd955h6lTp/L48WOj23wgVZokegByKwU1mtbNU/SSkpKYPHky1atXJz09nYsXL7JixQpsbW0JDAykadOmdOjQAQsLC/744w+OHTvGhx9++MpFDyBB+YB0tdLo4y2sFFR/t0Eur9TTp0/TpEkTU80zY8bMM7xBPb4HwDX0IQvGIgPqALZFYlFxEh8fz/Tp0/n999+ZPHkyPj4+aLVasrKy8lxUKlWudWm1SiDKFME8Z0wa1tEZKBQKFAoFSqWS/fv3c+zYMdq0acOAAQNwdnYmIiKCvXv3cuzYMRo3boyPjw8dO3bExsYGhUKBpaUlFhYWRhWiLWp23T5JREqMSa9FCqmcXh4t8ShVAQClUknZsmV58OABVlZWRWOoGTNm3iThuw4kFkE7rkDBLvFCCNRqdb4C8iKBMWSbMdsBFAoFGo0GIQT29vaUKVMGS0vLXMtTUXoqjhkZGdQf2RmHqi4m38HM8CQeh94mJSWFsLAwIiMjqVSpEpUrV0aj0RATE0N8fDxCCBwcHChZsiQ6nQ6VSpV9TU//rVKpkMvl2fY+a/uzS2HWG9PGaYs7pJNl0n2RSaS0c65Hg7LVAfjjjz+YMGECoaGhJt9zM2bM/EPx+Hi/FqiKpJWQkH189plvgeKjUqmQyWTZD8X8FkO3FyROhWn7qUu/EII9e/YwadIkypUrx/Dhw5HL5URHR3Pr1i1u3bpFeHg4sbGxODk54eHhwVtN38bexQmEABOcR2QSKS0bNWN36G1+/fVXfH192b17NxcuXGDNmjWcOXMGb29v/P39ady48QsdVYQQaDSaXGKYn0gWtP7ZdampqYVqo9N/PqJkWdNGAnRClyNM5PTp0zRu3NikNs2YMZObN0j4igaFQkG9evUoX748Li4uuLi4ULFiRRwdHbGyssoWHUMzoLxsMjMziYmJyRa0p+IWHR1NYmIiMTEx/Pnnn7i4uPD+++9Tt25devbsiaenJ25ubkjkUrZHHyf+yX29t6KJDpNZmVkM6Nqdtk1bsXPnTnbv3k3jxo2pVq0aQ4YMYdOmTYWas5NIJFhYWGBh8fLiLJ88eUJcXBzx8fHExcVlL8/+rUxNN1n4pBIpCuk/P8kzZ87Qo0cPU803Y8bMc7xBwmf63JQQgpIl7SldujR//fUXx44dIz4+nvj4eABcXFxwdnbO8fnsvytUqIBCUbSxgDqdjsTExByC9qzIJScnU6lSJTw8PPD09MTT05PGjRvj6emJh4cHdnZ2pKam8p///IfAwEDKly+f7TCSpVUTGL6HFFU6EpnUZNETOh2ZiY+p51WH3377je3bt9OwYUNWr15N165diz0MQQjBo0ePcghYXuKWmZmJq6srLi4uuLq64urqipeXFx07dsxefyrzJpGpcSbZI5NIsbUsmf336dOnmTdvnqmXacaMmed4I+b4hBBcu3aCKlUysbIypWcgBeryvHOLEIK0tLRsEUxISMjx+fTf9+7dw9bWNpcgPi+Sjo6OOUQgNTU1l6A9/bx9+zZlypTJFrJnPz09PXFxcTG493n79m2mTZvGiRMn6N69OyXaVaaUqz1yy6IRa61Kw7bpK/EoVZ4hQ4ZQoUIFtm3bxqZNm7C2tqZ///7079+/SNKK6XQ6kpKSChS0uLg4LCwsssXseXF7+m97e/sXinJh07jlhbVMwZhavZBKpMTHx/P222+TlJRkTlFmxkwR84p7fDr03YiX88O+desWGzZsICgoCEtLS0JDl5jYogIok2utRCKhdOnSlC5dGi8vr3yP1mq1JCUl5RLGkydPcuvWLWJjY0lKSiIrKwsrKyskEgkqlQqdToeDgwPOzs54enpmx695eXnh7u5eZO78SqUSZ2dnnjx5wrErZ+g7qEGRiZ46U4XioZrff96ZnXMToEGDBnz99decPn2a4OBgWrZsiYuLC/3796dfv355pulSq9XZ9y6/ocfExERsbW1ziVmHDh2y/3ZxcaFUqVJFcn2uJcpiI7PksZHCJ5PIaOBYHalE76F65swZg+Y4zZgxU3iKuccngIfo4+Ee809MnBQoC1QETHsQPX78mK1btxIUFERERAQ+Pj58+OGH1K9fH4kkFoimuAPYhRAkJyfnOxyZkJBAhQoVsntrFStWxNbWFktLS6RSKWlpaSQmJubqUVpZWb2w91iuXLkC81SmpaWxZcsWAgMDiYmJYfDgwQwZMoTzkjiiUhOK5J1EKiRUtXWhh3uL7Ad7Qfbs2LGDLVu2cOzYMRwdHXFzc8PGxobk5GTi4+N58OABTk5OBfbUnJ2diz0E4NKDKA7GnzMqc4uFVM4Ir66UtLAGICAggJIlSzJz5syiNtOMmTeeYhS+JOAv9HF0+T0YpIA18BZQMp99cqPRaAgJCWH9+vXs3buX9u3b4+fnR+fOnZ+bU9MCZ3kZKcuUSiW3b9/Oczjy1q1bWFpaZg8/Pj8kWalSpUI7ZwghePjwYa7h1OeHWJ+KxLOCWKFCBTIyMggLC+PPP/+kZcuWDB8+nM6dOyOXy0lXZ7AifBdaYWywvx6tRotCKqeFcx2aOHmRmpqaZ+8sh5OIUpktYs7OzmRlZXHnzh0iIiKoVasW/fr1w8/PD0dHR5NsexkIIdh95zTXU+68MEH1s8glMnp7tMKzdIXsda1bt+azzz6jY8eOL8NUM2beaIpJ+GKBWxje05KiDxTPnb/xWa5cucL69evZuHEjlSpVws/PDx8fnxeU5jEuSbVWW4+EhOR8he3Ro0e4u7vnmmPz8PDAw8ODMmVyD5EWB2q1mrt37xIfH09ERAS//vorJ06cQKVS4eTkhEQi4d69e2i12mxxrNGqHhU/eAepwjTPVFV6Jiv6TMe2dBni4uKQSqVUrFgxR+/s+d6ag4NDnsN7T548Yc+ePQQHB3P48GHatWuHj48PXbt2xcbGxiQ7ixKd0PF7bCjhKTEv7PlJkCCXSunh1oIqZf6Jj9RoNNja2hIXF4et7b8vWYIZM687xSB8d4EbFH54UYq+h5Wz53fv3j2Cg4NZv3499+/fx9fXF19f3wLn1nKjBq4CqejFL/ct0Gh0CCG4cOE2w4d/R0REJA4ODvn22ipUqPBaZBB5Hq1Wy4EDBwgMDOTQoUN0794df39/WrZsmUNg0tLSsnuKN9LjeFBOILM0MURAI5jXfiSbN2+mWbNm+ebmLCwpKSn897//JTg4mDNnztC5c2f69++fXbfuVSOE4K/Hcfx57xrJmY/RCh3ime+YhUSGAGrYVqRZuZo4WOV8Kbp48SIDBgwgPDy8mC03Y+bN4CULnwY4ifH5MUsAjcjMzOS3335j/fr1/PHHH3Tv3h0/Pz/atGlTqHg5lUqVI6YtLe0etWrZ0bZtdSQSCRqNFgsLOTqdIDQ0jjt3tDg6VsTDwwN3d/d/Vdqo6Oho1q5dy7p166hQoQL+/v74+Pi8sOeZnJzMpMWfU7l7IyysTBMRa5mCuxvPkZWVxZIlpjoW5U1SUhLbtm0jODiY8PBwevbsiY+PT6Gqyb9M7mc+5tKDKB5lpaPWqbGSW+JqU5baDh5YyfK+v6tWrSI0NJS1a9cWs7VmzLwZvGThiwOiMFb4tFr4+uvfWbz4R+rVq8eHH35Ir169KFky7/k/IQT37t3Ldzjy3r17uLq65tFrc6dyZTfs7EohkcjRO7u+fr23F5GZmcnOnTtZs2YNFy9eZODAgfj7+1OnTh2Djj9+/DgDBw6ki19vKvWoj8TCNOGwtyzF+yXqUq9ePW7fvl1kHpT5ERsby+bNm9m0aRNxcXH07dsXHx8fmjZt+q/yjvzoo49o0qQJw4cPf9WmmDHzf8lLFD4BnEbvSGIcGo2WiIiHlCnTJLssTXp6er7ekdHR0ZQqVSrPeDYPDw9cXV2LrRJ3cXLp0iUCAwMJDg6mXr16+Pv706NHD4N6qEIIwsPDCQgIICQkBJlMRg2vGvRaMRapwvh7JZfIaFm+Nk3KvUWfPn1o06YNo0aNMrq9whIZGcnmzZsJDg5GqVTi4+ODj48PdevWfe1F0MvLi02bNlG3bt1XbYoZM/+XvEThU6L3oDTNM1Cl0vLhh2uyRS49PT2XoD3rRJJfb/D/jZSUFIKDg1mzZg1JSUl89NFHfPTRRzkqk+fH3bt3CQkJISQkhP379/Po0SPs7OyYOXMmffr0oWzZshxJuEBo0nV0RtYbkEukfFqzJ9ZyS06cOMHQoUMJDw8v9nlQIQRXrlwhODiYTZs2YWVlhY+PD/3796datWrFaoshpKSkULFiRR49evR/+ZJmxszrwEv8ZakoiiAwCwspXbt+gIeHXtzKlSv32r+xvyyEEBw/fpw1a9bw66+/0rFjR+bOnUvHjh0LnM9KT0/n+PHjHDx4kJCQEOLj42nbti2urq7odDqmTJnCrFmzsh+0KpWKEz/vQ9LaBbkRDi5SJFQt44q13BKAFi1aYG1tzYEDB3jvvfeMu3gjkUgk1KlThzp16vCf//yHM2fOEBwcTOvWralQoUJ2oPzTEYVXTWhoKPXr1zeLnhkzL5GX2ONLAS5jWv27p7Tm3zjnVlQkJiayfv161qxZg0KhwN/fH19fX8qWLZvn/hqNhrNnzxISEsLBgwc5f/48DRs2pGPHjnTo0IG6desyd+5c1q5dS1BQEO3bt88+9tSpUwwdOhRPT09GL5rOxYyYQsWkAZSUW+Nf431s5P8Mtf70009s3bqVvXv3GncTihitVsuxY8cIDg5mx44deHl50b9/f7y9vSlXrtwrs+vLL78kPT2d+fPnvzIbzJj5f+clCl86cB7ThU8CtDHZmn8barWavXv3smbNGk6cOFFgqR4hBDdu3Mgevjx69Cju7u506NCBjh070qJFi+y0ZvHx8QwYMACFQsHPP/+c/ZBPS0tj+vTpbN++naVLl+Lt7Y1EImoxANwAABffSURBVOH7Axt5UFptkIenVqOllKUNvtU6YW+Z05ElMzMTNzc3Tpw48doNMapUKg4cOMCmTZvYvXs3jRo1wsfHh169ehV7HF2XLl34+OOP6dmzZ7Ge14yZN4mXKHw69KEMxift1WML1DPdnCJCCIFKp0EnBFYyiyIfdo2MjGTt2rWsX78eT09P/P396du3b665y3v37mULXUhICBKJJLtH1759e5ycnHK1vXfvXoYMGcLo0aOZNm1a9nzb7t27GTlyJJ06dWLhwoXY2ekTByQmJlKzZk1qtKtP31nD0ZKzXlz2PdHq53FTo+6hvZTMwq/yrijw2WefkZqaynfffWfSPXqZKJVKdu/ezaZNmzh06BBt2rShf//+dO3atchyoj5Lxr2HPI6MQ52WgdzGkhY93+fEtTCcnZ2L/FxmzJjR85LDGaLQZ20x9hQyoCZQUCaWl48Qgnjlfc4kRXDzsb4EkUQiQSt0lLe2p0m5t6hWxhXZC/JQ5odSqWTbtm0EBgZy48YNfH198ff3zxGU/3Se7qnQxcbG0rZtWzp06ECHDh2oWrVqviKsVquZMWMGwcHB/PLLL7Rs2RLQi+fYsWMJCwvjhx9+oF27djmu+YMPPuDUqVOsXr2a3r17czv9LqfvRRCXnkyWRoWlhQJruYKoo5eobetBt3e7UKdOHbZs2UKLFi1y2REfH0/t2rWJjo5+ZZlsCsPjx4+zA+VPnTqVI1De0tLS6HaFECQeucjVhZu5e+wS0r/nUXU6HZlpT3BpVpvaU3xw7dIY6WsQi2jGzP8bL1n4MoEzGO/ZaQE052VVbzCEROVD/nv7D55oMvMtOaOQypFIJHRwfoc6DpUNalcIQVhYGIGBgWzZsoWmTZvi7+/PBx98gEKhQKPRcO7cuWyHlPPnz9OgQYNsoTPUASImJgYfHx/s7e1Zv349jo6OCCFYt24dAQEB+Pv7M2vWLKytrXMcFxQUxPjx43n//ff5+eefc2wLCwtj2LBhhIWFAdC0aVMWLlxIixYt+O9//8vkyZO5ePFinj2kp3F1Y8eONeg+vS4kJydnB8pfu3aNHj160L9/f9q0aVMoR5Qn8cns7zSFJ7HJaNIz8t1PXsoaRekSvHtgAbZebkVxCWbMmPmbYkhZ9heQgHEpy7yA3EN2xcWt1ER2RB9HbaBzh1wio4mTFy0r5B8w/uDBAzZu3MiaNWtIS0tjyJAhDB48GBcXFyIjI7MdUo4dO0alSpWyhy9btmxZ6KG2Xbt2MWzYMCZPnsyECROQSqVERUUxbNgwHj9+TGBgIG+//Xau4+Lj46lZsyY2NjZERETk6p0dO3aMWbNmcezYMQCcnJy4ePFi9vCcr68vdnZ2eQ5p/vnnn/j5+REZGflapngzhLi4uOxA+djYWPr06ZMt6AVdU9rtu/zWaCSqlHSExoDvlESCRUkr3j+6BId6ptcoNGPGjJ5iePJUAewLeSop4MGrFL27yoeFEj0AjdByOimC8/cjc6zX6XSEhITg4+ND5cqVOX36NIsXL+bkyZN4enoya9Ys3Nzc6NChA2FhYfTr14/r169z6dIlFi1axHvvvVco0cvKymLcuHGMHTuWXbt2MWnSJHQ6HQsWLKBx48Z06dKF06dP5yl6QggGDx6MVqslODg4zyHJtLS07DnHtLQ00tPTqVDhn8oC3333HTt27ODIkSO5jm3atCl2dnbs27fP4Ot53XB1dWXixImcPXuWEydOUK5cOYYNG4aHhwcBAQFcuHCB598n1WlK9rUZj+phmmGiByAE6rQMfu8wmSfxyS/hSsyYeTMpBuGTALUA579PV9CwpfTvpTrwauOqfrtzqlCi9xSN0HIo/gIZGhWxsbHMmTOHypUrM2nSJBo2bMjq1aupUKECEyZM4K233mLr1q3Ur1+fkJAQ7ty5w9q1a+nfv7/RLvVRUVE0b96cmJgYLly4QJMmTTh//jyNGjUiJCSE0NBQJkyYkO/w3Lp16wgNDWXYsGG0bt06z33S09Ozhe/WrVt4enrmmF+0s7Pjxx9/ZMiQIaSlpeU4ViKRMGbMGJYuXWrU9b1uVK1alc8++4yrV6/y22+/IZPJ6NWrF15eXsyePZsbN24AcCNwL5nJjxG6wg/7q9OVXJm/qahNN2PmjaWYC9FmAvF/L88i0M/nVQLK86oLwycqH7Lxr4NGCR+ARAc3951j5+L1tGnTBicnJ8LDwzl//jz169fPnqdr0KBBkQYqb9myhU8//ZSZM2fy6aefkpGRweeff05QUBCLFi1i0KBBBXqhxsXF4eXlRfny5bl69Wq+DhyBgYGcOnWKNWvWsGPHDtavX8+uXbty7ffxxx8jk8n44YcfcqzPysrCzc2Nw4cP89Zbb5l20a8hQghCQ0MJDg5my5YtlHMqx6jbFZClZhndpryEFf2TdiC3Nt6pxowZM3qKWWGs0Fcx90Cf0kyNvgdogb7Q6+uRkSU0KQKNCUVYhRScm9dAu0jHzZs3cXd3Z+rUqbRq1erluMRnZDBhwgQOHjzIvn37qF+/PgcPHmT48OE0a9aMq1ev5hvsnm2zEAwaNAidTseOHTsK9FpMT0/PTjgdFRVF5cp5O/R888031KlTh/379/Puu+9mr7e0tGT48OEs/197dx9d07nnAfz77L3Pzsk5OXmRN5IQESRBUdVES5V66aXopcgwU5WaKhlSupg1tZZWV9dSt3NZbdGsxWhZt4PVYnW1xtDW5XbGpSVXSiISDiLeEvIm5yTHOftl/khlEhI5+4UrOb/Pn8fZz9ki8s3z7Of5/TZswGeffabjb/x4Y4whIyMDGRkZWLt2LQ58vA3X3tlucFDg4s5D6JP1aCvfENIZ/Z2mVhy0dFh/1M7dvtqif5oeNkcIjp46jrTuvU26q9YVFxdj5syZSE1NRV5eHiRJwty5c3H48GHk5uZiwoQJfo2zZcsWnDhxAu+99x6eeOKJB763+VKn0+ls8/1hYWHYsmULsrKycPr06RaHwRcsWID+/ftj9erVnbrZKs/ziKlUUC4Z+36SXB6U/dcxCj5CTNAxt9U9RKqqtnlsQQvRYoEt/OG24fnyyy8xYsQIZGdnY8eOHdi3bx8GDBiAiIgIFBQU+B16ZWVlWLJkCdLS0rBs2bJ2339v8LU14wOAsWPHYvLkyVi6dGmL17t164aJEycGRM85T3m1OePcqjVlHEICHVXCvYfRmd7DGqs5t9uNnJwcHDlyBAcPHkR4eDgmT56MsrIyfPPNN8jIyPD/HlUVmZmZUFUVu3bt8uuIgcvlQkJCAoD2gw8APvroIwwaNAh79+7FpEmTml7PyclBZmYm3nrrrceiaezDYqS9U4txLPTflRAz0IzvHhzjIDDjP4RVqG122DaisLAQ6enp8Pl8+Pnnn3H48GEMGTIEw4cPR15enqbQA4Dc3FycPHkSGzduRGKifwel7x5n8Pl8uHr1arvXhYSE4IsvvsCbb76JysrKptfT09MRGxuLvXv3arrnjsaWEA0mGP+essdHmXA3hBAKvlb0CDF+fpAxhuhg855dqaqKzz//HKNGjcLy5cuxbNkyvPjii9i9ezeOHDmCFStWwGLR1kLo8uXLWLZsGUaOHInXXnvN7+vuLnWWlpaiW7duEMX2A37kyJGYMWMGFi9e3OL1nJycx7p2pxmSpj8PzmA3eyEkGMlzxpt0R4QENgq+VgyL6QcLp39ZiWccnorsq7t2573q6urw6quvYt26dfj+++/hdDoxZswYzJs3D4cOHUJKSormMVVVxdSpU2GxWLB9+3ZNxbbv7uq8cOFCu8ucza1evRp5eXnYvXt302vTp09HUVERCgoKNN1/RxKW0h0RA3sZGkMMs6Pb6PsLDhBCtKPga0WPkBjDy5RDoswpMZWfn4+hQ4ciODgY69atw6xZs1BUVIRTp07hjTfe0F326+OPP27qTB4Zqa0I+N0Znz/P95qz2WzYunUrFi1ahIqKCgCAKIpYuHBhp5/1DfzXWRDs1vbf2Ao+OAgDls0M2AbMhJiNgq8VjDFM6J6u61mfhfEYGpUCh2gzdA+qqiI3Nxfjxo3D8uXLwfM8Xn/9dXz44YfYtWtXixJhWl26dAnvvPMOZsyYgYkTJ2q+Xm/wAY0ly+bMmYPs7Oymsl7z58/H119/3eL5X2fT4+VnETd2CHiNB9A50YKIAT2RumDyQ7ozQgIPBV8bkkPjMC7+KU3hZ2E8+oTFY3ScsSWp2tpaZGZmYtOmTVi1ahVWrVoFxhgKCwsNNyhVFAWTJk2Cw+HA5s2bdY1hJPgA4P3330dRURF27mwswxUbG4spU6Zgy5Ytuu6nI2Ach1E7VyLm2X4QbP6FH28VEZbaHeP3/wF8kPkbpQgJVBR8DzA4qjde7jkcIic88JmfwPjG53rRKZiSONzQktTx48cxZMgQ2Gw2JCYmYv369dixYwdyc3NN6WG3Zs0aFBcX47vvvoPNpm9WajT4rFYrtm3bhiVLluD69esAGje5bNy4EZJk/Azl44oPEjF+/x+Q+i+/h2C3QghpfemTtwWBt4roNesFTDq6AUERD/c8KCGB5hHX6uyYJEVGUc1lHCs/gxpvHTjGgwFQoMDCBDwdk4rBkcmwCfqe4QCNS5uffPIJVq9ejalTp2LPnj1YuHAhVqxYAatV/7jNOZ1OpKWlYcGCBYaeqdntdpSXl6Nr1664du0aQkNDdY2zcuVK5Ofn49tvvwVjDCNGjMDSpUvxyiuv6L63jsLnbsDFnYdw+o9fwXXxOhSvBM4iIDguEv1ypqFP1u8QFP74VjcipCOj4NPotrce9ZIHym/n9MJFO7jfdm/Kiozi2jKcrDwPl68BiqogiBeR5OiKIVF9ESa2XqezqqoKWVlZuHDhAoKDg8FxHDZv3txu6TAtFEVBSkoKvF4vnE6n7uLYsixDFEVcuXIFAwcOxM2b+tvleL1epKenY8mSJZg7dy6++uorbNy4sanPXyBRFQWsg/YnJKSjoVIQGoWKNoTes3HFp0j43xun8bdb5wGo8LYoeebGLU8tjt8sRnd7NEbHPYmuti5Nf3r06FHMmjUL8fHxuHbtGlatWoXs7GzTK5m8++67uHTpEgoKCgx1hHC73bDb7ZqPMrRGFEVs27atqdnu1KlT8fbbbyM/P7/VXoGdGYUeIY8O/W8zqF7yYGvJAZy4WQKv4rsn9BrJqgJZVXDJVY4/nfsBJTVlTY1hX3rpJaiqivDwcJw8eRKLFy82PfSKioqwZs0arFy5UteZv+buPt8zI/gAYNCgQcjJycG8efMgCAKys7Oxfv16w+MSQkhbaKnTAK8sYWvJflTfqYOioS4nDw5/+4/9+GnP95AkCZ9++ikyMzN1bYrxKRLOVJfifO1V1Mt3wDEGh8WGARFJSHJ0haqq6NmzJ0JCQlBYWGj4LFhJSQkmTZqE2bNnQ5ZlfPDBB4bGAwBJkjBs2DDMnz8f06ZNQ58+fXDu3DlERVGJLkKI+Wip04Afrp5AjdelKfQAQIaCtNkjEXyHwx/X/LvmA+QAUOetx1/LC3G66gLA2H0dJUpqr0DkBNw4WoJbVZU4duyYKQegm+/oHDNmjOHxAEAQBGzb1ti0d/z48Zg2bRo2bdqEFStWmDI+IYQ0R8Gn0x3ZhzPVpZB1Nqy1BluxaPW/6Qq9G/VV2OE8CK8sNYZuK7nrUyT4FAlBA2Lx3v5cRMRo/5zWNA+++fPnmzImAPTv3x/Lly9HVlYW1q5diylTpmDhkkU4776GWq8bPlmGzRKEOFskkkPjmjYUEUKIVrTUqdOJm8U4fP1XQ737Qi02ZPd7WdNM7JanFttKDrT6LLEtPDhEWB2Y2/dFQzVIAWDfvn3YsGED8vLycPLkScTFxRkarzlZlvHcc89h5oI5KOWrEZuWCI7nWvxyYeEECIzH09EpeDKqt6EjJISQwES/NuuUd6vEcMPaBtmLGw3+NylVVAU7nH/WFHpA49Jq9Z06HCg7rvUW71NXVwer1Yq6ujpDZdNaw3Ec3t70AW73DUJMvx5QOdw3o/YpEhrkO/hreSE2Fe1FuYavHyGEABR8url9HsNjMDDU+er9fr/z9jV4ZZ+uz5JVBUU1pWiQ7ui6/i6XywVFUdCrVy/TiyYfuHIcl+RbsFiDgHbGllQZDbIXX577ATcbaky9D0JI50bBp5OkyiaMomqaNR4tP6N5ttcSw69VTgPXNwafz+cz5ShDc6cqnSiovgifxq+rV5Gw3XnQ8OybEBI4KPh0MvqsDGic8QXx/jWPve1140ZDlaHPk1QZJ26WGBrD5XLB4/GYGnyqquKnG6fhU/T9MuFTZBRVl5p2P4SQzo2CT6fY4AjDY0iqjBirf13aa71uXW2S7uX2NRi63uVyweVymRp8l10V8Mhe3df7FAlHK86Ydj+EkM6Ngk+njJg0iAZnfQn2aIS2Ub/zXsaWOP+fAhWKziMYqqrC7alHbW2tqcH3y82zhpcq63z1uF5vbEZMCAkMdI5Pp16ObrBwgu5AEjkBw2LS/H6/v0ui7eEZp+kMXL3kQX6lEyduFsMtedD11acx59WncZarQ0RFEQZ2SUawYKxXXIUpOzMZKj216NasDiohhLSGgk8nxhhGdhuIH6/kad6QwYEhVLQjyeH/cYAI0WHKhpq2OkTcy6dI+O+yX3C25jIYWNNn3y2m3AAffrp+Cn+5/iue6JKE8fFDwXP6lmLN2JiiqgruKPp2vBJCAgstdRowOLI3BkYmw6Lh2RsHhmAhCLOTX9B0HMBusaJnSFc9t9nEwgnI8GOW2SB5sbXkAM7WXIasKm0GrqTKkFUFBVWX8KfzP8Ir6wswwYyNQoyZsuGIENL5UfAZNC7+KaTHpEJgPBgeHGQWTkB4UAheT5kAuyVY82dlxKQZ/OGuol94zwe+Q1Jk7HT+GdWeOr/LsUmqjIqGauy6+Bddzw/DReMNVxmYKeMQQjo/+hXZoMYlz0FICe+BXyrO4mzNZXCMQVJkqFDB/zYbjAhy4JnYfkgJ6w5B55Jgj5AYhIl2VHluay6MbWE8hkT2hcg/+J88v/I8bnlqIUNbgMmqgmv1lThTXYoBXZI0XZsenYLyhipDG3hEzoLu9mjd1xNCAgcFn0ligyMwOfEZjEt4Cudrr8IteSCrCqy8iHh7lCnHHxhjmJX8ArYU70OD5IXqZ/gJjEe8PRqj4gY98H2qquJYRZHuZ4k+RcKxijOag693WLyhotMC45Eek2p6JRlCSOdEwWcyKy9q/sGvRYglGHP7/g7/ef5HuH2edkPKwglIcnTFy4nD2w0Xo+fpAKD6jgvl9dWItfkf9BzjkB6diqMVhboOsTPGMLCLuZVkCCGdFz3j64DCRDv+OWUiRscNhsNiu++5HwfWOMuzRWFK4rOY1vM5v5ZXz1SXGt5hKasyimq0V1F5JrYfEuzREDTO/ATGY3rSSMNHKgghgYNmfB2UyFswNDoFT0X1xWVXBUpd5XD7POAZg0O0ISW8O7oEhWoaU0vB7LaoAFw6qsNwjMP0pOex59L/4HJdebtHRBgYBI7H73sOR0+Hsd2uhJDAQsHXwTHGkOiIRaIj1vBYWjfMtEVvc16B4zEj6Xn8WunE0YozcEue+2agjWXbVPQJS8CIrk8gyhpmwh0TQgIJBR9pYheCzBlHx1GNuxhjGBzVG4Mik3G1/hZ+rXQ2dmBXJFj5ICQ6YjCoSzKCTbpXQkjgoeAjTXqHxqOk9orBYwUCemmoSNMWxhgS7NFIoCMKhBCT0eYW0qRvWEK7h/DbI3IWJNEzN0LIY4yCjzThOR5PRvUBr/NMHZ2nI4R0BBR8pIWMmFRdnSAYGOyCFYMjez+EuyKEEPNQ8JEWbIIV/9R7LKy8xe9lTw4MwbyIf+wz1rT2SYQQ8rAwVVXN2cNOOpWaOy7sdB6CS2p44KF2CycgXAzBPySPRoiB3ZyEEPKoUPCRNqmqilJXOX6uKEKpq/y3M3SNJFVGL0cchsWkId4eRc/1CCEdBgUf8Uu95EGt1w2vLEHkLQgX7XSWjhDSIVHwEUIICSi0uYUQQkhAoeAjhBASUCj4CCGEBBQKPkIIIQGFgo8QQkhAoeAjhBASUCj4CCGEBBQKPkIIIQGFgo8QQkhAoeAjhBASUCj4CCGEBBQKPkIIIQGFgo8QQkhAoeAjhBASUP4P2HrHefQ0n1gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "graph = nx.from_scipy_sparse_matrix(adjacency)\n",
    "\n",
    "cmap = plt.cm.get_cmap('Spectral')\n",
    "colors = cmap(np.arange(labels.max()))\n",
    "nx.draw(graph, node_color=[cmap(i) for i in labels/labels.max()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing equivalence to `ward_cut_tree_balanced`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from scipy.cluster.hierarchy import ward, cut_tree\n",
    "from scipy.stats import gamma\n",
    "import numpy as np\n",
    "\n",
    "def ward_cut_tree_balanced(linkage_matrix_Z, max_cluster_size, verbose=False):\n",
    "    \"\"\"This function performs a balanced clustering by using the linkage matrix from a Ward histogram. \n",
    "       It builds upon the scipy and numpy libraries. \n",
    "       \n",
    "       The function looks recursively along the hierarchical tree, from the root (single cluster gathering \n",
    "       all the samples) to the leaves (i.e. the clusters with only one sample), retrieving the biggest \n",
    "       possible clusters containing a number of samples lower than a given maximum. In this way, if a \n",
    "       cluster at a specific tree level contains a number of samples higher than the given maximum, it is \n",
    "       ignored and its offspring (smaller) sub-clusters are taken into consideration. If the cluster contains \n",
    "       a number of samples lower than the given maximum, it is taken as result and its offspring sub-clusters \n",
    "       not further processed.\n",
    "       Input parameters:\n",
    "          \n",
    "          linkage_matrix_Z: linkage matrix resulting from calling the method scipy.cluster.hierarchy.ward()\n",
    "             I.e. it contains the hierarchical clustering encoded as a linkage matrix.\n",
    "          max_cluster_size: maximum number of data samples contained within the resulting clusters. Thus, all \n",
    "             resulting clusters will contain a number of data samples <= max_cluster_size.\n",
    "             Note that max_cluster_size must be >= 1.\n",
    "          verbose: activates (True) / deactivates (False) some output print commands, which can be useful to \n",
    "             test and understand the proposed tree cut method.\n",
    "       \n",
    "       Returns:\n",
    "          vec_cluster_id: one-dimensional numpy array of integers containing for each input sample its corresponding \n",
    "             cluster id. The cluster id is an integer which is higher for deeper tree levels.\n",
    "          vec_last_cluster_level: one-dimensional numpy array of arrays containing for each input sample its \n",
    "             corresponding cluster tree level, i.e. a sequence of 0s and 1s. Note that the cluster level is longer for \n",
    "             deeper tree levels, being [0] the root cluster, [0, 0] and [0, 1] its offspring, and so on. Also note that \n",
    "             in each cluster splitting, the label 0 denotes the bigger cluster, while the label 1 denotes the smallest.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        # Assert that the input max_cluster_size is >= 1\n",
    "        assert max_cluster_size >= 1\n",
    "        \n",
    "        # Perform a full cut tree of the linkage matrix, i.e. containing all tree levels\n",
    "        full_cut = cut_tree(linkage_matrix_Z)\n",
    "        if verbose:\n",
    "            print(\"Interim full cut tree (square matrix)\")\n",
    "            print(\"Shape = \" + str(full_cut.shape))\n",
    "            print(full_cut)\n",
    "            print('')\n",
    "                \n",
    "        # Initialize the vble containing the current cluster id (it will be higher for each newly \n",
    "        # found valid cluster, i.e. for each found cluster with <= max_cluster_size data samples)\n",
    "        last_cluster_id = 1\n",
    "        \n",
    "        # Initialize the resulting cluster id vector (containing for each row in input_data_x_sample \n",
    "        # its corresponding cluster id)\n",
    "        vec_cluster_id = np.zeros(full_cut.shape[1], dtype=int)\n",
    "        \n",
    "        # Initialize the resulting cluster level vector (containing for each data sample its \n",
    "        # corresponding cluster tree level, i.e. a string of '0's and '1's separated by '.')\n",
    "        vec_last_cluster_level = np.empty((full_cut.shape[1],), dtype=object)\n",
    "        for i in range(full_cut.shape[1]): vec_last_cluster_level[i] = np.array([0],int)\n",
    "\n",
    "        # Scan the full cut matrix from the last column (root tree level) to the first column (leaves tree level)\n",
    "        if verbose:\n",
    "            print(\"Note about columns: within the full cut tree, the column \" + str(full_cut.shape[1]-1) +\n",
    "                  \" represents the root, while 0 represent the leaves.\")\n",
    "            print(\"We now scan the full cut tree from the root (column \" + str(full_cut.shape[1]-1) + \") \"\n",
    "                  \"to the leaves (column 0).\")\n",
    "            print('')\n",
    "\n",
    "        for curr_column in range(full_cut.shape[1]-1,-1,-1):\n",
    "            \n",
    "            # Get a list of unique group ids and their count within the current tree level\n",
    "            values, counts = np.unique(full_cut[:,curr_column], return_counts=True)\n",
    "            \n",
    "            # Stop if all samples have been already selected (i.e. if all data samples have been already clustered)\n",
    "            if (values.size==1) and (values[0]==-1):\n",
    "                break\n",
    "                \n",
    "            # For each group id within the current tree level\n",
    "            for curr_elem_pos in range(values.size):\n",
    "            \n",
    "                # If it is a valid group id (i.e. not yet marked as processed with -1) ...\n",
    "                # Note: data samples which were alredy included in a valid cluster id (i.e. at a higher tree level) \n",
    "                # are marked with the group id -1 (see below)\n",
    "                if (values[curr_elem_pos] >= 0):\n",
    "                \n",
    "                    # Select the current group id\n",
    "                    selected_curr_value = values[curr_elem_pos]\n",
    "                    \n",
    "                    # Look for the vector positions (related to rows in input_data_x_sample) belonging to \n",
    "                    # the current group id\n",
    "                    selected_curr_elems = np.where(full_cut[:,curr_column]==selected_curr_value)\n",
    "                    \n",
    "                    # Major step #1: Populate the resulting vector of cluster levels for each data sample\n",
    "                    # If we are not at the root level (i.e. single cluster gathering all the samples) ...\n",
    "                    if curr_column < (full_cut.shape[1]-1):\n",
    "                        # Get the ancestor values and element positions\n",
    "                        selected_ancestor_value = full_cut[selected_curr_elems[0][0],curr_column+1]\n",
    "                        selected_ancestor_elems = np.where(full_cut[:,curr_column+1]==selected_ancestor_value)\n",
    "                        \n",
    "                        # Compute the values and counts of the offspring (i.e. curr_elem + brothers) and sort them\n",
    "                        # by their count (so that the biggest cluster gets the offspring_elem_label = 0, see below)\n",
    "                        offspring_values, offspring_counts = np.unique(full_cut[selected_ancestor_elems,curr_column], \n",
    "                                                                       return_counts=True)\n",
    "                        count_sort_ind = np.argsort(-offspring_counts)\n",
    "                        offspring_values = offspring_values[count_sort_ind]\n",
    "                        offspring_counts = offspring_counts[count_sort_ind]\n",
    "                        \n",
    "                        # If the number of descendants is > 1 (i.e. if the curr_elem has at least one brother)\n",
    "                        if (offspring_values.shape[0] > 1):\n",
    "                            # Select the position of the current value (i.e. 0 or 1) and append it to the cluster level\n",
    "                            offspring_elem_label = np.where(offspring_values==selected_curr_value)[0][0]\n",
    "                            for i in selected_curr_elems[0]:\n",
    "                                vec_last_cluster_level[i] = np.hstack((vec_last_cluster_level[i], offspring_elem_label))\n",
    "\n",
    "                    # Major step #2: Populate the resulting vector of cluster ids for each data sample, \n",
    "                    # and mark them as already clustered (-1)\n",
    "                    # If the number of elements is below max_cluster_size ...\n",
    "                    if (counts[curr_elem_pos] <= max_cluster_size):\n",
    "                    \n",
    "                        if verbose:\n",
    "                            print(\"Current column in full cut tree = \" + str(curr_column))\n",
    "                            print(\"list_group_ids:     \" + str(values))\n",
    "                            print(\"list_count_samples: \" + str(counts))\n",
    "                            print(\"selected_curr_value: \" + str(selected_curr_value) + \", count_samples = \" + \n",
    "                                  str(counts[curr_elem_pos]) + \", marked as result\")\n",
    "                            print('')\n",
    "                \n",
    "                        # Relate these vector positions to the current cluster id \n",
    "                        vec_cluster_id[selected_curr_elems] = last_cluster_id\n",
    "                    \n",
    "                        # Delete these vector positions at the lower tree levels for further processing \n",
    "                        # (i.e. mark these elements as already clustered)\n",
    "                        full_cut[selected_curr_elems,0:curr_column] = -1\n",
    "                    \n",
    "                        # Update the cluster id\n",
    "                        last_cluster_id = last_cluster_id + 1\n",
    "                    \n",
    "        # Return the resulting clustering array (containing for each row in input_data_x_sample its \n",
    "        # corresponding cluster id) and the clustering level\n",
    "        return vec_cluster_id, vec_last_cluster_level\n",
    "\n",
    "    except AssertionError:\n",
    "        print(\"Please use a max_cluster_size >= 1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "homogeneity_score(ward_cut_tree_balanced(dendrogram, 10)[0], balanced_cut(dendrogram, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Better example\n",
    "Larger graphs give more exaggerated differences between the biggest cluster and the long tail of smaller clusters that occurs with single straight cuts. \n",
    "\n",
    "Using the wikivitals graph here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = skn.data.load_wikilinks_dataset('wikivitals')\n",
    "adjacency = data.adjacency\n",
    "paris = skn.hierarchy.Paris(engine='python')\n",
    "dendrogram = paris.fit_transform(adjacency)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adjacency = skn.data.miserables()\n",
    "#paris = skn.hierarchy.Paris(engine='python')\n",
    "#dendrogram = paris.fit_transform(adjacency)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, a straight cut:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using max cluster size of 10\n",
      "Cluster ID 0 has 1234 items\n",
      "Cluster ID 1 has 1003 items\n",
      "Cluster ID 2 has 912 items\n",
      "Cluster ID 3 has 687 items\n",
      "Cluster ID 4 has 659 items\n",
      "Cluster ID 5 has 590 items\n",
      "Cluster ID 6 has 586 items\n",
      "Cluster ID 7 has 582 items\n",
      "Cluster ID 8 has 542 items\n",
      "Cluster ID 9 has 503 items\n",
      "Cluster ID 10 has 360 items\n",
      "Cluster ID 11 has 334 items\n",
      "Cluster ID 12 has 313 items\n",
      "Cluster ID 13 has 309 items\n",
      "Cluster ID 14 has 287 items\n",
      "Cluster ID 15 has 286 items\n",
      "Cluster ID 16 has 252 items\n",
      "Cluster ID 17 has 224 items\n",
      "Cluster ID 18 has 219 items\n",
      "Cluster ID 19 has 130 items\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 20 artists>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASaElEQVR4nO3dfYxc13nf8e+vpK34JY6laqXQJFHKBeGWMpJaIVQlbg0DCirGMkw1iAIaTUw0CggHchMHKWqqBuL8Q0BpWqNxUblgY8V0q1hmHbsi4qixwNYwCkRSV7JsiaJl0REjbcSQm6S13QZQQvnpH3MFjFcz3OXeeVnxfD/AYu6ce+6ch2eG89t778zdVBWSpHb9jXkXIEmaL4NAkhpnEEhS4wwCSWqcQSBJjds87wJWc+WVV9aOHTvmXYYkvaI88sgjf1ZVC2vpu+GDYMeOHSwuLs67DEl6RUnyx2vt66EhSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq3Ib/ZnEfOw5+4aL6n77z5ilVIkkbl3sEktS4VYMgyd1JziV5YqjtN5J8PcnXknw+yRuH1t2R5FSSp5LcNNT+I0ke79Z9LEkm/8+RJF2stewRfBLYs6LtAeCtVfVDwDeAOwCS7AL2Add229yVZFO3zceBA8DO7mflY0qS5mDVIKiqLwN/saLti1V1vrv7ILCtW94L3FtVL1TVM8Ap4PokW4A3VNUfVlUBnwJumdQ/QpK0fpM4R/BzwP3d8lbguaF1S13b1m55ZftISQ4kWUyyuLy8PIESJUnj9AqCJB8GzgP3vNQ0oltdoH2kqjpcVburavfCwpr+roIkaZ3W/fHRJPuBdwM3dod7YPCb/vahbtuA57v2bSPaJUlztq49giR7gA8B76mqvxxadQzYl+SyJNcwOCn8cFWdAb6T5Ibu00LvA+7rWbskaQJW3SNI8mngncCVSZaAjzD4lNBlwAPdp0AfrKr3V9WJJEeBJxkcMrq9ql7sHuoXGHwC6TUMzincjyRp7lYNgqp674jmT1yg/yHg0Ij2ReCtF1WdJGnq/GaxJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDVu1SBIcneSc0meGGq7IskDSZ7ubi8fWndHklNJnkpy01D7jyR5vFv3sSSZ/D9HknSx1rJH8Elgz4q2g8DxqtoJHO/uk2QXsA+4ttvmriSbum0+DhwAdnY/Kx9TkjQHm1frUFVfTrJjRfNe4J3d8hHgS8CHuvZ7q+oF4Jkkp4Drk5wG3lBVfwiQ5FPALcD9vf8FU7Lj4BcuepvTd948hUokabrWe47g6qo6A9DdXtW1bwWeG+q31LVt7ZZXto+U5ECSxSSLy8vL6yxRkrQWkz5ZPOq4f12gfaSqOlxVu6tq98LCwsSKkyS93HqD4GySLQDd7bmufQnYPtRvG/B8175tRLskac7WGwTHgP3d8n7gvqH2fUkuS3INg5PCD3eHj76T5Ibu00LvG9pGkjRHq54sTvJpBieGr0yyBHwEuBM4muQ24FngVoCqOpHkKPAkcB64vape7B7qFxh8Auk1DE4Sb9gTxZLUkrV8aui9Y1bdOKb/IeDQiPZF4K0XVZ0kaer8ZrEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDVu1W8Wa30u9u8Z+LcMJM2LewSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXG9giDJLyc5keSJJJ9O8n1JrkjyQJKnu9vLh/rfkeRUkqeS3NS/fElSX+sOgiRbgV8EdlfVW4FNwD7gIHC8qnYCx7v7JNnVrb8W2APclWRTv/IlSX31PTS0GXhNks3Aa4Hngb3AkW79EeCWbnkvcG9VvVBVzwCngOt7ji9J6mndQVBVfwL8a+BZ4Azwrar6InB1VZ3p+pwBruo22Qo8N/QQS12bJGmO+hwaupzBb/nXAG8CXpfkZy60yYi2GvPYB5IsJllcXl5eb4mSpDXoc2jox4Fnqmq5qv4a+BzwY8DZJFsAuttzXf8lYPvQ9tsYHEp6mao6XFW7q2r3wsJCjxIlSavpEwTPAjckeW2SADcCJ4FjwP6uz37gvm75GLAvyWVJrgF2Ag/3GF+SNAGb17thVT2U5LPAo8B54CvAYeD1wNEktzEIi1u7/ieSHAWe7PrfXlUv9qxfktTTuoMAoKo+AnxkRfMLDPYORvU/BBzqM6YkabL8ZrEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJalyvy1BrOnYc/MJF9T99581TqkRSC9wjkKTGGQSS1DiDQJIa5zmCS8zFnl8AzzFIrXOPQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDWuVxAkeWOSzyb5epKTSX40yRVJHkjydHd7+VD/O5KcSvJUkpv6ly9J6qvvHsFvAv+tqv4O8MPASeAgcLyqdgLHu/sk2QXsA64F9gB3JdnUc3xJUk/rDoIkbwDeAXwCoKr+qqr+D7AXONJ1OwLc0i3vBe6tqheq6hngFHD9eseXJE1Gnz2CNwPLwG8n+UqS30ryOuDqqjoD0N1e1fXfCjw3tP1S1/YySQ4kWUyyuLy83KNESdJq+gTBZuA64ONV9Tbg/9EdBhojI9pqVMeqOlxVu6tq98LCQo8SJUmr6XOtoSVgqaoe6u5/lkEQnE2yparOJNkCnBvqv31o+23A8z3G1xT0+VsI/h0F6ZVp3UFQVX+a5Lkkb6mqp4AbgSe7n/3And3tfd0mx4DfSfJR4E3ATuDhPsXr0mKQSPPR9+qj/wy4J8mrgT8C/imDw01Hk9wGPAvcClBVJ5IcZRAU54Hbq+rFnuNLknrqFQRV9Riwe8SqG8f0PwQc6jOmJGmy/GaxJDXOIJCkxhkEktQ4g0CSGuffLNYlwY+eSuvnHoEkNc4gkKTGGQSS1DiDQJIaZxBIUuP81JCad7GfOAI/daRLi3sEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxvYMgyaYkX0nye939K5I8kOTp7vbyob53JDmV5KkkN/UdW5LU3yT2CH4JODl0/yBwvKp2Ase7+yTZBewDrgX2AHcl2TSB8SVJPfQKgiTbgJuB3xpq3gsc6ZaPALcMtd9bVS9U1TPAKeD6PuNLkvrru0fwb4F/AXx3qO3qqjoD0N1e1bVvBZ4b6rfUtb1MkgNJFpMsLi8v9yxRknQh6w6CJO8GzlXVI2vdZERbjepYVYerandV7V5YWFhviZKkNejzN4vfDrwnybuA7wPekOQ/A2eTbKmqM0m2AOe6/kvA9qHttwHP9xhfkjQB694jqKo7qmpbVe1gcBL4v1fVzwDHgP1dt/3Afd3yMWBfksuSXAPsBB5ed+WSpInos0cwzp3A0SS3Ac8CtwJU1YkkR4EngfPA7VX14hTGlyRdhIkEQVV9CfhSt/znwI1j+h0CDk1iTEnSZExjj0Bqyo6DX7io/qfvvHlKlUjr4yUmJKlxBoEkNc4gkKTGGQSS1DhPFktz5IlmbQTuEUhS4wwCSWqch4akV6iLPawEHlrSaAaB1CjPT+glHhqSpMYZBJLUOINAkhrnOQJJF63P+QVPcm887hFIUuMMAklqnEEgSY0zCCSpcQaBJDXOTw1JekXxG9GT5x6BJDXOIJCkxhkEktQ4zxFIaobfah5t3UGQZDvwKeAHge8Ch6vqN5NcAXwG2AGcBn66qv53t80dwG3Ai8AvVtUf9KpekmboUj1R3efQ0HngV6rq7wI3ALcn2QUcBI5X1U7geHefbt0+4FpgD3BXkk19ipck9bfuIKiqM1X1aLf8HeAksBXYCxzpuh0BbumW9wL3VtULVfUMcAq4fr3jS5ImYyIni5PsAN4GPARcXVVnYBAWwFVdt63Ac0ObLXVtox7vQJLFJIvLy8uTKFGSNEbvIEjyeuB3gQ9W1bcv1HVEW43qWFWHq2p3Ve1eWFjoW6Ik6QJ6BUGSVzEIgXuq6nNd89kkW7r1W4BzXfsSsH1o823A833GlyT1t+4gSBLgE8DJqvro0KpjwP5ueT9w31D7viSXJbkG2Ak8vN7xJUmT0ed7BG8HfhZ4PMljXdu/BO4Ejia5DXgWuBWgqk4kOQo8yeATR7dX1Ys9xpckTcC6g6Cq/iejj/sD3Dhmm0PAofWOKUmaPC8xIUmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhrn3yOQpBnYyJewdo9AkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjZh4ESfYkeSrJqSQHZz2+JOl7zTQIkmwC/j3wE8Au4L1Jds2yBknS95r1HsH1wKmq+qOq+ivgXmDvjGuQJA1JVc1usOSngD1V9fPd/Z8F/n5VfWBFvwPAge7uW4CnJlzKlcCfTfgxJ2Gj1gUbt7aNWhds3No2al2wcWvbqHXB+Nr+VlUtrOUBNk+2nlVlRNvLkqiqDgOHp1ZEslhVu6f1+Ou1UeuCjVvbRq0LNm5tG7Uu2Li1bdS6YDK1zfrQ0BKwfej+NuD5GdcgSRoy6yD4X8DOJNckeTWwDzg24xokSUNmemioqs4n+QDwB8Am4O6qOjHLGjpTO+zU00atCzZubRu1Lti4tW3UumDj1rZR64IJ1DbTk8WSpI3HbxZLUuMMAklq3CUdBKtdziIDH+vWfy3JdTOoaXuS/5HkZJITSX5pRJ93JvlWkse6n1+ddl1DY59O8ng37uKI9fOYs7cMzcVjSb6d5IMr+sxszpLcneRckieG2q5I8kCSp7vby8dsO7VLrIyp6zeSfL17rj6f5I1jtr3g8z6l2n4tyZ8MPWfvGrPtrOfsM0M1nU7y2JhtpzZn494npvY6q6pL8ofByehvAm8GXg18Fdi1os+7gPsZfL/hBuChGdS1BbiuW/5+4Bsj6non8HtzmrfTwJUXWD/zORvxvP4pgy/LzGXOgHcA1wFPDLX9K+Bgt3wQ+PUxtV/wNTmFuv4RsLlb/vVRda3leZ9Sbb8G/PM1PN8znbMV6/8N8KuznrNx7xPTep1dynsEa7mcxV7gUzXwIPDGJFumWVRVnamqR7vl7wAnga3THHPCZj5nK9wIfLOq/niGY36Pqvoy8BcrmvcCR7rlI8AtIzad6iVWRtVVVV+sqvPd3QcZfHdn5sbM2VrMfM5ekiTATwOfntR4a3WB94mpvM4u5SDYCjw3dH+Jl7/hrqXP1CTZAbwNeGjE6h9N8tUk9ye5dlY1Mfim9xeTPJLBpT5WmuucMfjuybj/mPOaM4Crq+oMDP4TA1eN6DPvufs5Bntzo6z2vE/LB7rDVnePOcwxzzn7h8DZqnp6zPqZzNmK94mpvM4u5SBYy+Us1nTJi2lI8nrgd4EPVtW3V6x+lMGhjx8G/h3wX2dRU+ftVXUdgyvE3p7kHSvWz3POXg28B/gvI1bPc87Wap5z92HgPHDPmC6rPe/T8HHgbwN/DzjD4DDMSnObM+C9XHhvYOpztsr7xNjNRrRdcM4u5SBYy+Us5nLJiySvYvDk3lNVn1u5vqq+XVX/t1v+feBVSa6cdl3deM93t+eAzzPYzRw2z8uE/ATwaFWdXblinnPWOfvSIbLu9tyIPvN6ve0H3g38k+oOIq+0hud94qrqbFW9WFXfBf7jmDHnNWebgZ8EPjOuz7TnbMz7xFReZ5dyEKzlchbHgPd1n4S5AfjWS7td09Idd/wEcLKqPjqmzw92/UhyPYPn6c+nWVc31uuSfP9LywxOND6xotvM52zI2N/Q5jVnQ44B+7vl/cB9I/rM/BIrSfYAHwLeU1V/OabPWp73adQ2fG7pH48Zc16Xpflx4OtVtTRq5bTn7ALvE9N5nU3jjPdG+WHwCZdvMDiD/uGu7f3A+7vlMPhDOd8EHgd2z6Cmf8BgN+1rwGPdz7tW1PUB4ASDs/0PAj82o/l6czfmV7vxN8ScdeO+lsEb+w8Mtc1lzhiE0Rngrxn89nUb8DeB48DT3e0VXd83Ab9/odfklOs6xeB48Uuvtf+wsq5xz/sMavtP3WvoawzeqLZshDnr2j/50mtrqO/M5uwC7xNTeZ15iQlJatylfGhIkrQGBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq3P8Hhg1Il4CnufAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = skn.hierarchy.straight_cut(dendrogram, n_clusters=20)\n",
    "ids, counts = np.unique(labels, return_counts=True)\n",
    "print(f'Using max cluster size of {maxClusterSize}')\n",
    "for j,k in zip(ids, counts):\n",
    "    print(f'Cluster ID {j} has {k} items')\n",
    "    \n",
    "plt.bar(np.arange(counts.size), counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Dont try and draw the wikivitals graph its too big.\n",
    "\n",
    "#graph = nx.from_scipy_sparse_matrix(adjacency)\n",
    "#cmap = plt.cm.get_cmap('Spectral')\n",
    "#colors = cmap(np.arange(labels.max()))\n",
    "#nx.draw(graph, node_color=[cmap(i) for i in labels/labels.max()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, a balanced cut:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using max cluster size of 10\n",
      "Cluster ID 0 has 541 items\n",
      "Cluster ID 1 has 649 items\n",
      "Cluster ID 2 has 44 items\n",
      "Cluster ID 3 has 294 items\n",
      "Cluster ID 4 has 709 items\n",
      "Cluster ID 5 has 659 items\n",
      "Cluster ID 6 has 687 items\n",
      "Cluster ID 7 has 584 items\n",
      "Cluster ID 8 has 417 items\n",
      "Cluster ID 9 has 542 items\n",
      "Cluster ID 10 has 590 items\n",
      "Cluster ID 11 has 599 items\n",
      "Cluster ID 12 has 582 items\n",
      "Cluster ID 13 has 219 items\n",
      "Cluster ID 14 has 812 items\n",
      "Cluster ID 15 has 586 items\n",
      "Cluster ID 16 has 586 items\n",
      "Cluster ID 17 has 912 items\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x116e32f50>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASXElEQVR4nO3df5BdZ33f8fenNhAbQrGR7ApJVKYj3MqdieNqXNs0jGfkqX+C3IJBmdqRW3dEwKSoAxMkAiEzYFBaYJx0YhclECuYYqs2qYWxW4waT6aDsSsbAZaFbBE71saKJLsN0IFxbOfbP+5x5mq1K+3u3b27q+f9mtm55z7Pc/Z899HRZ8+ec++5qSokSW34O7NdgCRpeAx9SWqIoS9JDTH0Jakhhr4kNeTE2S7gWBYsWFDLli2b7TIkaV55+OGHn62qhaPb53zoL1u2jB07dsx2GZI0ryT587HaPb0jSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGnLM0E/yxSQHkzza13ZqkvuSPNE9ntLXtzHJ3iR7klzc1/5Pkny/6/vdJJn+H0eSdDQTOdK/BbhkVNsGYHtVLQe2d89JsgJYA5zVrXNTkhO6dW4G1gHLu6/R31OSNMOO+easqvrTJMtGNa8GLuyWtwD3Ax/u2m+rqueBJ5PsBc5N8hTw2qp6ACDJHwFXAvcea/t79uzhwgsvPKztXe96F+973/v46U9/ymWXXXbEOtdeey3XXnstzz77LO985zuP6H/ve9/Lu9/9bvbt28c111xzRP8HP/hB3va2t7Fnzx7e8573HNH/0Y9+lIsuuoidO3eyfv36I/o/9alPccEFF/Ctb32Lj3zkI0f033jjjZx99tl885vf5JOf/OQR/Z///Oc588wz+drXvsZnP/vZI/q/9KUvsXTpUm6//XZuvvnmI/rvuOMOFixYwC233MItt9xyRP8999zDySefzE033cTWrVuP6L///vsB+MxnPsPdd999WN9JJ53Evff2/tk+8YlPsH379sP6X//613PnnXcCsHHjRh544IHD+pcsWcKtt94KwPr169m5c+dh/W9+85vZvHkzAOvWrePxxx8/rP/ss8/mxhtvBODqq69mZGTksP7zzz+fT3/60wC84x3v4Lnnnjusf9WqVXzsYx8D4NJLL+VnP/vZYf1XXHEFH/rQhwCO2O/Afc99b/7uey+b6jn906tqP0D3eFrXvhjY1zdupGtb3C2Pbh9TknVJdiTZ8cILL0yxREnSaJnIJ2d1R/p3V9U/7p7/VVW9rq///1bVKUl+D3igqm7t2r8A3AM8DXy6qi7q2n8J+PWqetuxtr1y5cryNgySNDlJHq6qlaPbp3qkfyDJou4bLwIOdu0jwNK+cUuAZ7r2JWO0S5KGaKqhvw1Y2y2vBe7qa1+T5FVJzqB3wfah7hTQT5Kc171q51f61pEkDckxL+Qm+Qq9i7YLkowAHwc2AVuTXEfv1M1VAFW1K8lW4DHgReD6qnqp+1bvpfdKoJPoXcA95kVcSdL0mtA5/dnkOX1JmrzpPqcvSZqHDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhx7zLpiRpMMs2fH3S6zy16fIZqMQjfUlqiqEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhA4V+kn+fZFeSR5N8JcnPJTk1yX1JnugeT+kbvzHJ3iR7klw8ePmSpMmY8gejJ1kM/DtgRVX9LMlWYA2wAtheVZuSbAA2AB9OsqLrPwt4A/DNJG+uqpcG/ikkHVem8kHiMHMfJn48GfT0zonASUlOBE4GngFWA1u6/i3Ald3yauC2qnq+qp4E9gLnDrh9SdIkTDn0q+ovgM8ATwP7gR9V1TeA06tqfzdmP3Bat8piYF/ftxjp2o6QZF2SHUl2HDp0aKolSpJGmXLod+fqVwNn0Dtd8+okVx9tlTHaaqyBVbW5qlZW1cqFCxdOtURJ0iiDnN65CHiyqg5V1QvAV4ELgANJFgF0jwe78SPA0r71l9A7HSRJGpJBQv9p4LwkJycJsArYDWwD1nZj1gJ3dcvbgDVJXpXkDGA58NAA25ckTdKUX71TVQ8muQN4BHgR+A6wGXgNsDXJdfR+MVzVjd/VvcLnsW789b5yR5KGa8qhD1BVHwc+Pqr5eXpH/WONvwG4YZBtSpKmznfkSlJDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNWSgu2xqfvHDpiV5pC9JDTH0Jakhhr4kNcTQl6SGeCFXs2IqF5W9oCwNziN9SWqIoS9JDTH0Jakhhr4kNcQLuZoU39UrzW8e6UtSQzzSPwaPbCUdTzzSl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0ZKPSTvC7JHUl+kGR3kvOTnJrkviRPdI+n9I3fmGRvkj1JLh68fEnSZAx6pP87wH+vqn8I/AKwG9gAbK+q5cD27jlJVgBrgLOAS4Cbkpww4PYlSZMw5dBP8lrgrcAXAKrqr6vqr4DVwJZu2Bbgym55NXBbVT1fVU8Ce4Fzp7p9SdLkDXKk/ybgEPCHSb6T5A+SvBo4var2A3SPp3XjFwP7+tYf6dokSUMyyG0YTgTOAX6tqh5M8jt0p3LGkTHaasyByTpgHcAb3/jGAUqU2uEtQzQRgxzpjwAjVfVg9/wOer8EDiRZBNA9Huwbv7Rv/SXAM2N946raXFUrq2rlwoULByhRktRvykf6VfWXSfYlObOq9gCrgMe6r7XApu7xrm6VbcB/SfI54A3AcuChQYqXjhcepU+/6ZrT4+3znAe9y+avAV9O8krgz4B/Te+vh61JrgOeBq4CqKpdSbbS+6XwInB9Vb004PbVsOPtP+Nc4C+f499AoV9VO4GVY3StGmf8DcANg2xTkjR1viNXkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGDPrmrDnNN5poGHyTmOYTj/QlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNOXG2C5Bm07INX5/0Ok9tunwGKpGGY+Aj/SQnJPlOkru756cmuS/JE93jKX1jNybZm2RPkosH3bYkaXKm4/TOB4Ddfc83ANurajmwvXtOkhXAGuAs4BLgpiQnTMP2JUkTNFDoJ1kCXA78QV/zamBLt7wFuLKv/baqer6qngT2AucOsn1J0uQMeqR/I/DrwN/0tZ1eVfsBusfTuvbFwL6+cSNd2xGSrEuyI8mOQ4cODViiJOllUw79JFcAB6vq4YmuMkZbjTWwqjZX1cqqWrlw4cKplihJGmWQV++8BXh7ksuAnwNem+RW4ECSRVW1P8ki4GA3fgRY2rf+EuCZAbYvSZqkKR/pV9XGqlpSVcvoXaD9n1V1NbANWNsNWwvc1S1vA9YkeVWSM4DlwENTrlySNGkz8Tr9TcDWJNcBTwNXAVTVriRbgceAF4Hrq+qlGdi+JGkc0xL6VXU/cH+3/BywapxxNwA3TMc2JUmT520YJKkhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNmYlPztIMWLbh61Na76lNl09zJZLmM4/0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDvLWypGnlbcDntikf6SdZmuRPkuxOsivJB7r2U5Pcl+SJ7vGUvnU2JtmbZE+Si6fjB5AkTdwgp3deBD5YVf8IOA+4PskKYAOwvaqWA9u753R9a4CzgEuAm5KcMEjxkqTJmXLoV9X+qnqkW/4JsBtYDKwGtnTDtgBXdsurgduq6vmqehLYC5w71e1LkiZvWi7kJlkG/CLwIHB6Ve2H3i8G4LRu2GJgX99qI13bWN9vXZIdSXYcOnRoOkqUJDENoZ/kNcCdwPqq+vHRho7RVmMNrKrNVbWyqlYuXLhw0BIlSZ2BQj/JK+gF/per6qtd84Eki7r+RcDBrn0EWNq3+hLgmUG2L0manEFevRPgC8DuqvpcX9c2YG23vBa4q699TZJXJTkDWA48NNXtS5Imb5DX6b8FuAb4fpKdXdtHgE3A1iTXAU8DVwFU1a4kW4HH6L3y5/qqemmA7UuSJmnKoV9V/4uxz9MDrBpnnRuAG6a6TUnSYLwNgyQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUkBNnu4AWLNvw9Smt99Smy6e5Ekmt80hfkhpi6EtSQwx9SWqIoS9JDTH0JakhQw/9JJck2ZNkb5INw96+JLVsqKGf5ATg94BLgRXALydZMcwaJKllwz7SPxfYW1V/VlV/DdwGrB5yDZLUrFTV8DaWvBO4pKr+bff8GuCfVtX7R41bB6zrnp4J7JmBchYAz87A950J86XW+VInWOtMmC91Qhu1/v2qWji6cdjvyM0YbUf81qmqzcDmGS0k2VFVK2dyG9NlvtQ6X+oEa50J86VOaLvWYZ/eGQGW9j1fAjwz5BokqVnDDv3/DSxPckaSVwJrgG1DrkGSmjXU0ztV9WKS9wP/AzgB+GJV7RpmDX1m9PTRNJsvtc6XOsFaZ8J8qRMarnWoF3IlSbPLd+RKUkMMfUlqyHEd+se65UN6frfr/16Sc2apzqVJ/iTJ7iS7knxgjDEXJvlRkp3d12/ORq1dLU8l+X5Xx44x+ufKvJ7ZN187k/w4yfpRY2ZtXpN8McnBJI/2tZ2a5L4kT3SPp4yz7tBuZzJOnf8xyQ+6f98/TvK6cdY96r4ypFp/K8lf9P0bXzbOukO9Rcw4td7eV+dTSXaOs+7U57WqjssveheKfwi8CXgl8F1gxagxlwH30nv/wHnAg7NU6yLgnG7554HHx6j1QuDu2Z7XrpangAVH6Z8T8zrG/vCX9N6wMifmFXgrcA7waF/bfwA2dMsbgN8e52c56r49hDr/OXBit/zbY9U5kX1lSLX+FvChCewfQ5vT8Wod1f9Z4Dene16P5yP9idzyYTXwR9XzbeB1SRYNu9Cq2l9Vj3TLPwF2A4uHXcc0mhPzOsoq4IdV9eezXMffqqo/Bf7PqObVwJZueQtw5RirDvV2JmPVWVXfqKoXu6ffpveem1k3zpxOxNBvEXO0WpMEeBfwlene7vEc+ouBfX3PRzgySCcyZqiSLAN+EXhwjO7zk3w3yb1JzhpqYYcr4BtJHu5umTHanJtXeu8JGe8/0FyZV4DTq2o/9A4GgNPGGDPX5vff0PvLbizH2leG5f3dqagvjnPKbK7N6S8BB6rqiXH6pzyvx3PoT+SWDxO6LcSwJHkNcCewvqp+PKr7EXqnJn4B+E/Afxt2fX3eUlXn0Ltb6vVJ3jqqf67N6yuBtwP/dYzuuTSvEzVn5jfJbwAvAl8eZ8ix9pVhuBn4B8DZwH56p01GmzNz2vlljn6UP+V5PZ5DfyK3fJgzt4VI8gp6gf/lqvrq6P6q+nFV/b9u+R7gFUkWDLnMl2t5pns8CPwxvT+N+82Zee1cCjxSVQdGd8ylee0cePlUWPd4cIwxc2J+k6wFrgD+VXUnmkebwL4y46rqQFW9VFV/A/z+ODXMiTkFSHIi8C+B28cbM8i8Hs+hP5FbPmwDfqV7tcl5wI9e/tN6mLrzd18AdlfV58YZ8/e6cSQ5l96/3XPDq/Jv63h1kp9/eZneBb1HRw2bE/PaZ9yjprkyr322AWu75bXAXWOMmfXbmSS5BPgw8Paq+uk4Yyayr8y4UdeT/sU4Ncz6nPa5CPhBVY2M1TnwvM7k1enZ/qL3KpLH6V2V/42u7VeBX+2WQ+9DXX4IfB9YOUt1/jN6f0p+D9jZfV02qtb3A7vovarg28AFs1Trm7oavtvVM2fntavlZHoh/nf72ubEvNL7RbQfeIHekeZ1wOuB7cAT3eOp3dg3APccbd8ecp176Z0Df3l//c+j6xxvX5mFWr/U7Yffoxfki2Z7TsertWu/5eX9s2/stM2rt2GQpIYcz6d3JEmjGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIf8fsLOvLjXbDt4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "maxClustSize=1000\n",
    "labels = balanced_cut(dendrogram, max_cluster_size=maxClustSize)\n",
    "ids, counts = np.unique(labels, return_counts=True)\n",
    "print(f'Using max cluster size of {maxClusterSize}')\n",
    "for j,k in zip(ids, counts):\n",
    "    print(f'Cluster ID {j} has {k} items')\n",
    "    \n",
    "plt.bar(np.arange(counts.size), counts)\n",
    "plt.axhline(maxClustSize, linestyle='--', c='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# test equivalence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000000000000002"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "homogeneity_score(ward_cut_tree_balanced(dendrogram, 1000)[0], balanced_cut(dendrogram, 1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.07 s ± 112 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "ward_cut_tree_balanced(dendrogram, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26.9 ms ± 614 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "balanced_cut(dendrogram, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}