distance_distributions.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Investigating MO distance distributions in MAC\n",
    "\n",
    "In the proof of concept notebook, we notice that the MO 'energy' distances follow a peculiar, almost discontinuous distribution. The purpose of this notebook is to get to elucidate the origin of this discontinuity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy spectrum contains degeneracies =  False\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "84a0b367f7594c3e8791be549272e912",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz4klEQVR4nO3de3QUZZ7/8U/HJG0MSZsYk6aRmzJe2OAtKBdFUCHgEBSddVQwC7OYFTWwCKxHdH4DsmpcRNxdENFZQZ2jxHUQ9SgbiYggEgQDWQhej4LBkBCE0B0iuZA8vz8cam3CnSTdyfN+nVPnkKpvd3+rxPSHp+qpchljjAAAAGCNiFA3AAAAgNZFAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAy0SGuoG2rLGxUTt37lRcXJxcLleo2wEAACfAGKOqqir5fD5FRNg5FkYAPA07d+5U586dQ90GAAA4BTt27NB5550X6jZCggB4GuLi4iT98hcoPj4+xN0AAIATEQgE1LlzZ+d73EYEwNNw6LRvfHw8ARAAgDbG5su37DzxDQAAYDECIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAoHm4XL8sCHsEQAAAAMsQAAEAACwTGeoGAABAO/Pr08DGhK4PHBUjgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAQCAABOD/f+a3MYAQQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACwTGeoGAABAG+VyhboDnCJGAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAAAtx+VitnAYIgACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWKbZAmBOTo6uuuoqxcXFKTk5WSNHjtTXX38dVGOM0YwZM+Tz+RQTE6NBgwZp69atQTW1tbWaMGGCkpKSFBsbq5tvvlk//vhjUE1lZaUyMzPl8Xjk8XiUmZmpffv2BdWUlJRoxIgRio2NVVJSkiZOnKi6urqgmi1btmjgwIGKiYlRp06dNHPmTBljmuuQAAAAhKVmC4CrVq3SAw88oHXr1ik/P18HDx5Uenq6qqurnZpZs2Zpzpw5mjdvnjZs2CCv16shQ4aoqqrKqZk0aZKWLl2q3NxcrVmzRvv371dGRoYaGhqcmlGjRqmoqEh5eXnKy8tTUVGRMjMzne0NDQ0aPny4qqurtWbNGuXm5mrJkiWaMmWKUxMIBDRkyBD5fD5t2LBBc+fO1ezZszVnzpzmOiQAAADhybSQiooKI8msWrXKGGNMY2Oj8Xq95qmnnnJqampqjMfjMQsWLDDGGLNv3z4TFRVlcnNznZrS0lITERFh8vLyjDHGfPHFF0aSWbdunVNTUFBgJJmvvvrKGGPMsmXLTEREhCktLXVqFi9ebNxut/H7/cYYY+bPn288Ho+pqalxanJycozP5zONjY0ntI9+v99Ict4TAACrSCe+hBG+v41psWsA/X6/JCkxMVGStG3bNpWXlys9Pd2pcbvdGjhwoNauXStJKiwsVH19fVCNz+dTamqqU1NQUCCPx6M+ffo4NX379pXH4wmqSU1Nlc/nc2qGDh2q2tpaFRYWOjUDBw6U2+0Oqtm5c6e2b9/enIcCAAAgrLRIADTGaPLkybr22muVmpoqSSovL5ckpaSkBNWmpKQ428rLyxUdHa2EhIRj1iQnJzf5zOTk5KCawz8nISFB0dHRx6w59POhmsPV1tYqEAgELQAAAG1NiwTA7Oxsbd68WYsXL26yzXXY42CMMU3WHe7wmiPVN0eN+dsEkKP1k5OT40w88Xg86ty58zH7BgAACEfNHgAnTJigd999VytXrtR5553nrPd6vZKajq5VVFQ4I29er1d1dXWqrKw8Zs2uXbuafO7u3buDag7/nMrKStXX1x+zpqKiQlLTUcpDpk2bJr/f7yw7duw4xpEAAKAdOvRsX57v26Y1WwA0xig7O1tvvfWWPvroI3Xv3j1oe/fu3eX1epWfn++sq6ur06pVq9S/f39JUlpamqKiooJqysrKVFxc7NT069dPfr9f69evd2o+++wz+f3+oJri4mKVlZU5NcuXL5fb7VZaWppTs3r16qBbwyxfvlw+n0/dunU74j663W7Fx8cHLQAAAG1Oc80mue+++4zH4zEff/yxKSsrc5aff/7ZqXnqqaeMx+Mxb731ltmyZYu56667TMeOHU0gEHBqxo8fb8477zzz4Ycfmo0bN5obbrjBXHbZZebgwYNOzbBhw8yll15qCgoKTEFBgenVq5fJyMhwth88eNCkpqaaG2+80WzcuNF8+OGH5rzzzjPZ2dlOzb59+0xKSoq56667zJYtW8xbb71l4uPjzezZs094n5lFBACwzsnM/GUWcNhqtv8iko64LFq0yKlpbGw006dPN16v17jdbnPdddeZLVu2BL3PgQMHTHZ2tklMTDQxMTEmIyPDlJSUBNXs2bPHjB492sTFxZm4uDgzevRoU1lZGVTzww8/mOHDh5uYmBiTmJhosrOzg275YowxmzdvNgMGDDBut9t4vV4zY8aME74FjDH8BQIAWIgA2C64jOHRF6cqEAjI4/HI7/dzOhgAYIdTvfYvjOIG3988CxgAAMA6BEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwTGSoGwAAAG2AyxXqDtCMGAEEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAA0PJcLm4lE0YIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAlokMdQMAACCM8fSOdokRQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDKRoW4AAACEIZcr1B2gBTECCAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAluFZwAAA4Bc8/9cajAACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFimWQPg6tWrNWLECPl8PrlcLr399ttB28eOHSuXyxW09O3bN6imtrZWEyZMUFJSkmJjY3XzzTfrxx9/DKqprKxUZmamPB6PPB6PMjMztW/fvqCakpISjRgxQrGxsUpKStLEiRNVV1cXVLNlyxYNHDhQMTEx6tSpk2bOnCljTLMdDwAAgHDUrAGwurpal112mebNm3fUmmHDhqmsrMxZli1bFrR90qRJWrp0qXJzc7VmzRrt379fGRkZamhocGpGjRqloqIi5eXlKS8vT0VFRcrMzHS2NzQ0aPjw4aqurtaaNWuUm5urJUuWaMqUKU5NIBDQkCFD5PP5tGHDBs2dO1ezZ8/WnDlzmvGIAADQBrhcPAbONqaFSDJLly4NWjdmzBhzyy23HPU1+/btM1FRUSY3N9dZV1paaiIiIkxeXp4xxpgvvvjCSDLr1q1zagoKCowk89VXXxljjFm2bJmJiIgwpaWlTs3ixYuN2+02fr/fGGPM/PnzjcfjMTU1NU5NTk6O8fl8prGx8YT20e/3G0nOewIA0CZJrb+EEN/fxrT6NYAff/yxkpOTdeGFFyorK0sVFRXOtsLCQtXX1ys9Pd1Z5/P5lJqaqrVr10qSCgoK5PF41KdPH6emb9++8ng8QTWpqany+XxOzdChQ1VbW6vCwkKnZuDAgXK73UE1O3fu1Pbt24/Ye21trQKBQNACAADQ1rRqALzpppv02muv6aOPPtIzzzyjDRs26IYbblBtba0kqby8XNHR0UpISAh6XUpKisrLy52a5OTkJu+dnJwcVJOSkhK0PSEhQdHR0cesOfTzoZrD5eTkONcdejwede7c+WQPAQAAQMhFtuaH3XHHHc6fU1NT1bt3b3Xt2lXvv/++brvttqO+zhgj16+uTXAd4TqF5qgxf5sAcqTXStK0adM0efJk5+dAIEAIBAAAbU5IbwPTsWNHde3aVd9++60kyev1qq6uTpWVlUF1FRUVzuic1+vVrl27mrzX7t27g2oOH8WrrKxUfX39MWsOnY4+fGTwELfbrfj4+KAFAACgrQlpANyzZ4927Nihjh07SpLS0tIUFRWl/Px8p6asrEzFxcXq37+/JKlfv37y+/1av369U/PZZ5/J7/cH1RQXF6usrMypWb58udxut9LS0pya1atXB90aZvny5fL5fOrWrVuL7TMAAECoNWsA3L9/v4qKilRUVCRJ2rZtm4qKilRSUqL9+/dr6tSpKigo0Pbt2/Xxxx9rxIgRSkpK0q233ipJ8ng8GjdunKZMmaIVK1Zo06ZNuvvuu9WrVy8NHjxYknTJJZdo2LBhysrK0rp167Ru3TplZWUpIyNDF110kSQpPT1dPXv2VGZmpjZt2qQVK1Zo6tSpysrKckbtRo0aJbfbrbFjx6q4uFhLly7Vk08+qcmTJx/1FDAAAEC70JxTileuXGkkNVnGjBljfv75Z5Oenm7OPfdcExUVZbp06WLGjBljSkpKgt7jwIEDJjs72yQmJpqYmBiTkZHRpGbPnj1m9OjRJi4uzsTFxZnRo0ebysrKoJoffvjBDB8+3MTExJjExESTnZ0ddMsXY4zZvHmzGTBggHG73cbr9ZoZM2ac8C1gjGEaOQCgneA2MNZxGcOjL05VIBCQx+OR3+/nekAAQNsTyjNeIYwffH/zLGAAAADrEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsExkqBsAAACtjCdeWY8RQAAAAMsQAAEAACxDAAQAALAMARAAAMAyTAIBAMAGTPzArzACCAAAYBkCIAAAgGUIgAAAAJbhGkAAANozrv3DETACCAAAYBkCIAAAgGUIgAAAAJbhGkAAANojrv3DMTACCAAAYBkCIAAAgGUiNm7cqC1btjgr3nnnHY0cOVKPPPKI6urqQtgaAAAAWkLEvffeq2+++UaS9P333+vOO+/UWWedpTfffFMPPfRQiNsDAAAnxeXi+j8cV8Q333yjyy+/XJL05ptv6rrrrtPrr7+ul19+WUuWLAltdwAAAGh2EcYYNTY2SpI+/PBD/fa3v5Ukde7cWT/99FMoewMAAEALiOjdu7cef/xx/eUvf9GqVas0fPhwSdK2bduUkpIS4vYAAADQ3CKeffZZbdy4UdnZ2Xr00UfVo0cPSdJf//pX9e/fP8TtAQAAoLm5jDHmSBtqamoUGRmpyEjuFX00gUBAHo9Hfr9f8fHxoW4HAGCrtjjp48jxo1Xw/S1FnH/++dqzZ0+TDTU1NbrwwgtD0BIAAABaUsT27dvV0NDQZENtba1+/PHHELQEAACAlhQpSR988IE8Ho+zsqGhQStWrFD37t1D1hgAAABaRqTL5dKYMWOCVkZFRalbt2565plnQtQWAAAAWkpkly5dtGHDBiUlJYW6FwAAcDLa4uQPhIXIbdu2hboHAAAAtKJISVqxYoVWrFihiooK56kghyxcuDAkjQEAAKBlRD722GOaOXOmevfurY4dO8rFcDIAAEC7FrlgwQK9/PLLyszMDHUvAAAAaAURdXV1PPINAADAIhH33HOPXn/99VD3AQAATpTLxQxgnJbImpoavfjii/rwww916aWXKioqKqhgzpw5IWoNAAC0W78OsCF8LrCtIjdv3qzLL79cklRcXBy0kQkhAAAA7U/kypUrQ90DAAAAWlFEqBsAAABA64q8/vrrj3mq96OPPmrFdgAAANDSIg9d/3dIfX29ioqKVFxcrDFjxoSmKwAA0BTX5qOZRD777LNH3DBjxgzt37+/ldsBAABASzvqNYB33303zwEGAABoh44aAAsKCnTmmWe2Zi8AAABoBZG33XZb0ApjjMrKyvT555/r//2//xeitgAAANBSIj0eT9CKiIgIXXTRRZo5c6bS09ND1BYAAABaSuSiRYtC3QMAAABaUeShPxQWFurLL7+Uy+VSz549dcUVV4SyLwAAIHHrF7SIyIqKCt155536+OOPdfbZZ8sYI7/fr+uvv165ubk699xzQ90jAAAAmlHEhAkTFAgEtHXrVu3du1eVlZUqLi5WIBDQxIkTQ90fAAAAmpkrPj7efPjhh7rqqquCNqxfv17p6enat29faDprAwKBgDwej/x+v+Lj40PdDgCgPbLhFLAxrfpxfH9LEY2NjYqKimqyISoqSo2NjSFoCQAAAC0p4oYbbtA///M/a+fOnc7K0tJSPfjgg7rxxhtD2BoAAABaQsS8efNUVVWlbt266YILLlCPHj3UvXt3VVVVae7cuaHuDwAAO7lcdpz+RUhEdu7cWRs3blR+fr6++uorGWPUs2dPDR48ONS9AQAAoAW4/H6/OfwCSL/fr/79+2vBggUaMGBAiFoLf1xECgBoMTaN/jEJpNVFHGnHPR6P7r33Xs2ZMycELQEAAKAlRRxtQ3p6ugoLC1uzFwAAALSCowbAyMhI7d69uzV7AQAAQCs4agDcvHmzOnbs2Jq9AAAAoBVE1NTUNFl54MABTZ8+XRkZGSFoCQAAi3H7F7SCiAsvvFCzZs3SO++8o3fffVf/9m//posuukh79+7Vo48+elJvtnr1ao0YMUI+n08ul0tvv/120HZjjGbMmCGfz6eYmBgNGjRIW7duDaqpra3VhAkTlJSUpNjYWN1888368ccfg2oqKyuVmZkpj8cjj8ejzMzMJo+sKykp0YgRIxQbG6ukpCRNnDhRdXV1QTVbtmzRwIEDFRMTo06dOmnmzJkyrTwTCQAAoLVFpKamatq0abr11ls1cuRIPfLII0pNTdWnn36qlJSUk3qz6upqXXbZZZo3b94Rt8+aNUtz5szRvHnztGHDBnm9Xg0ZMkRVVVVOzaRJk7R06VLl5uZqzZo12r9/vzIyMtTQ0ODUjBo1SkVFRcrLy1NeXp6KioqUmZnpbG9oaNDw4cNVXV2tNWvWKDc3V0uWLNGUKVOcmkAgoCFDhsjn82nDhg2aO3euZs+ezcxnAADQ/hljzN69e8369evNZ599Zvbu3WuagySzdOlS5+fGxkbj9XrNU0895ayrqakxHo/HLFiwwBhjzL59+0xUVJTJzc11akpLS01ERITJy8szxhjzxRdfGElm3bp1Tk1BQYGRZL766itjjDHLli0zERERprS01KlZvHixcbvdxu/3G2OMmT9/vvF4PKampsapycnJMT6fzzQ2Np7QPvr9fiPJeU8AAE7bL3fFs2tpZXx/GxMhSQkJCbrqqqt09dVXKyEhoUWC5rZt21ReXq709HRnndvt1sCBA7V27VpJUmFhoerr64NqfD6fUlNTnZqCggJ5PB716dPHqenbt688Hk9QTWpqqnw+n1MzdOhQ1dbWOre2KSgo0MCBA+V2u4Nqdu7cqe3btzf/AQAAAAgTR50F3NzKy8slqclp5ZSUFGdbeXm5oqOjm4TQw2uSk5ObvH9ycnJQzeGfk5CQoOjo6GPWHPr5UM3hamtrFQgEghYAAIC2ptUC4CGuw2Y2GWOarDvc4TVHqm+OGvO3CSBH6ycnJ8eZeOLxeNS5c+dj9g0AABCOWi0Aer1eSU1H1yoqKpyRN6/Xq7q6OlVWVh6zZteuXU3ef/fu3UE1h39OZWWl6uvrj1lTUVEhqeko5SHTpk2T3+93lh07dhx/xwEAOBHc/gWtqNUCYPfu3eX1epWfn++sq6ur06pVq9S/f39JUlpamqKiooJqysrKVFxc7NT069dPfr9f69evd2o+++wz+f3+oJri4mKVlZU5NcuXL5fb7VZaWppTs3r16qBbwyxfvlw+n0/dunU74j643W7Fx8cHLQAAAG1NswbA/fv3q6ioSEVFRZJ+mfhRVFSkkpISuVwuTZo0SU8++aSWLl2q4uJijR07VmeddZZGjRolSfJ4PBo3bpymTJmiFStWaNOmTbr77rvVq1cvDR48WJJ0ySWXaNiwYcrKytK6deu0bt06ZWVlKSMjQxdddJGkX55j3LNnT2VmZmrTpk1asWKFpk6dqqysLCe0jRo1Sm63W2PHjlVxcbGWLl2qJ598UpMnTz7uKWkAAIA2rTmnFK9cudJIarKMGTPGGPPLrWCmT59uvF6vcbvd5rrrrjNbtmwJeo8DBw6Y7Oxsk5iYaGJiYkxGRoYpKSkJqtmzZ48ZPXq0iYuLM3FxcWb06NGmsrIyqOaHH34ww4cPNzExMSYxMdFkZ2cH3fLFGGM2b95sBgwYYNxut/F6vWbGjBknfAsYY5hGDgBoRqG+FQu3gbGKy5i/zXzASQsEAvJ4PPL7/ZwOBgCcHpvPPrVyFOH7OwSzgAEAABBakaFuAAAAa9k86oeQYgQQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMtwGxgAAFobt39BiDECCAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZZgEDANBamP2LMMEIIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAUMAEBLY/YvwgwjgAAAILRcLkJyKyMAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiG28AAANASmNWKMMYIIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAUMAEBzYvYv2gBGAAEAACxDAAQAALAMARAAAMAyBEAAAADLMAkEAIDmwOQPtCGMAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhlnAAACcDmb/og1iBBAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMswCBgDgZDHzF20cI4AAAACWIQACAABYhgAIAABgGQIgAACAZZgEAgDAiWLyB9oJRgABAAAsQwAEAACwDAEQAADAMlwDCADA8XDtH9oZRgABAAAsQwAEAACwDAEQAADAMgRAAAAAyzAJBACAo2HyB9opRgABAAAsQwAEAACwDAEQAADAMq0aAGfMmCGXyxW0eL1eZ7sxRjNmzJDP51NMTIwGDRqkrVu3Br1HbW2tJkyYoKSkJMXGxurmm2/Wjz/+GFRTWVmpzMxMeTweeTweZWZmat++fUE1JSUlGjFihGJjY5WUlKSJEyeqrq6uxfYdAAAgXLT6CODf/d3fqayszFm2bNnibJs1a5bmzJmjefPmacOGDfJ6vRoyZIiqqqqcmkmTJmnp0qXKzc3VmjVrtH//fmVkZKihocGpGTVqlIqKipSXl6e8vDwVFRUpMzPT2d7Q0KDhw4erurpaa9asUW5urpYsWaIpU6a0zkEAAIQvl+v/FqC9Mq1o+vTp5rLLLjvitsbGRuP1es1TTz3lrKupqTEej8csWLDAGGPMvn37TFRUlMnNzXVqSktLTUREhMnLyzPGGPPFF18YSWbdunVOTUFBgZFkvvrqK2OMMcuWLTMRERGmtLTUqVm8eLFxu93G7/ef8P74/X4j6aReAwAIcxJLqJZWwve3Ma0+Avjtt9/K5/Ope/fuuvPOO/X9999LkrZt26by8nKlp6c7tW63WwMHDtTatWslSYWFhaqvrw+q8fl8Sk1NdWoKCgrk8XjUp08fp6Zv377yeDxBNampqfL5fE7N0KFDVVtbq8LCwpbbeQAAgDDQqvcB7NOnj1599VVdeOGF2rVrlx5//HH1799fW7duVXl5uSQpJSUl6DUpKSn64YcfJEnl5eWKjo5WQkJCk5pDry8vL1dycnKTz05OTg6qOfxzEhISFB0d7dQcSW1trWpra52fA4HAie46AABA2GjVAHjTTTc5f+7Vq5f69eunCy64QK+88or69u0rSXIdds2FMabJusMdXnOk+lOpOVxOTo4ee+yxY/YCAGijuOYPFgnpbWBiY2PVq1cvffvtt85s4MNH4CoqKpzROq/Xq7q6OlVWVh6zZteuXU0+a/fu3UE1h39OZWWl6uvrm4wM/tq0adPk9/udZceOHSe5xwAAAKEX0gBYW1urL7/8Uh07dlT37t3l9XqVn5/vbK+rq9OqVavUv39/SVJaWpqioqKCasrKylRcXOzU9OvXT36/X+vXr3dqPvvsM/n9/qCa4uJilZWVOTXLly+X2+1WWlraUft1u92Kj48PWgAAANqaVj0FPHXqVI0YMUJdunRRRUWFHn/8cQUCAY0ZM0Yul0uTJk3Sk08+qd/85jf6zW9+oyeffFJnnXWWRo0aJUnyeDwaN26cpkyZonPOOUeJiYmaOnWqevXqpcGDB0uSLrnkEg0bNkxZWVl64YUXJEn/9E//pIyMDF100UWSpPT0dPXs2VOZmZl6+umntXfvXk2dOlVZWVmEOgAA0O61agD88ccfddddd+mnn37Sueeeq759+2rdunXq2rWrJOmhhx7SgQMHdP/996uyslJ9+vTR8uXLFRcX57zHs88+q8jISP3+97/XgQMHdOONN+rll1/WGWec4dS89tprmjhxojNb+Oabb9a8efOc7WeccYbef/993X///brmmmsUExOjUaNGafbs2a10JAAAAELHZYwxoW6irQoEAvJ4PPL7/YwcAkBbxeSP8NFKkYTvb54FDAAAYB0CIAAAgGVa9RpAAADCAqd9YTlGAAEAACxDAAQAALAMARAAAMAyXAMIALAH1/4BkhgBBAAAsA4BEAAAwDIEQAAAAMtwDSAAoP3j2j8gCCOAAAAAliEAAgAAWIZTwACA9otTv23Lof9exoS2DwswAggAAGAZAiAAAIBlCIAAAACW4RpAAED7wnV/wHExAggAAGAZAiAAAIBlOAUMAGgfOPULnDBGAAEAACxDAAQAALAMp4ABAG0bp36Bk8YIIAAAgGUYAQQAtE2M/AGnjBFAAAAAyxAAAQAALMMpYABA28FpX6BZMAIIAABgGQIgAACAZTgFDAAIf5z6BZoVI4AAAACWYQQQABC+GPkDWgQjgAAAAJZhBBAAEH4Y+QNaFCOAAAAAliEAAgAAWIZTwACA8MBpX6DVEAABAKFF8ANaHaeAAQAALMMIIAAgNBj5A0KGEUAAAADLMAIIAGhdjPwBIccIIAAAgGUYAQQAtDxG/YCwwgggAACAZRgBBAC0HEb+gLBEAAQAND+CHxDWCIAAgOZD8APaBK4BBAAAsAwjgACA08OoH9DmEAABAKeG4Ae0WQRAAMDJIfgBbR4BEABwYgh+QLtBAAQAHB2hD2iXCIAAgKYIfkC7RgAEAPwfgh9gBQIgAIDgB1iGAAgAtiL0AdYiAAKAbQh+gPUIgADQ3hH4AByGAAgA7RGhD8AxEAABoD0g8AE4CQRAAGhrCHto7w79HTcmtH20YwRAAAh3BD4AzYwACADhgJAHoBURAAGgtRDyAISJiFA3EGrz589X9+7ddeaZZyotLU2ffPJJqFsC0Ba5XMdfACBMWB0A33jjDU2aNEmPPvqoNm3apAEDBuimm25SSUlJqFsDEA5OJNQR7gC0QS5j7J1i06dPH1155ZV6/vnnnXWXXHKJRo4cqZycnOO+PhAIyOPxyO/3Kz4+viVbBXAqCGZA29ZCEYXvb4uvAayrq1NhYaEefvjhoPXp6elau3btEV9TW1ur2tpa52e/3y/pl79IAH7F4wl1BwDagxb6fj30vW3xGJi9AfCnn35SQ0ODUlJSgtanpKSovLz8iK/JycnRY4891mR9586dW6RHAACs1sL/mKyqqpLH0n+wWhsAD3EddorIGNNk3SHTpk3T5MmTnZ8bGxu1d+9enXPOOUd9jfTLvzQ6d+6sHTt2WDvUfDwco2Pj+Bwfx+j4OEbHxzE6vvZwjIwxqqqqks/nC3UrIWNtAExKStIZZ5zRZLSvoqKiyajgIW63W263O2jd2WeffcKfGR8f32b/Z2ktHKNj4/gcH8fo+DhGx8cxOr62foxsHfk7xNpZwNHR0UpLS1N+fn7Q+vz8fPXv3z9EXQEAALQ8a0cAJWny5MnKzMxU79691a9fP7344osqKSnR+PHjQ90aAABAi7E6AN5xxx3as2ePZs6cqbKyMqWmpmrZsmXq2rVrs36O2+3W9OnTm5w+xv/hGB0bx+f4OEbHxzE6Po7R8XGM2ger7wMIAABgI2uvAQQAALAVARAAAMAyBEAAAADLEAABAAAsQwBsQU888YT69++vs84666g3jC4pKdGIESMUGxurpKQkTZw4UXV1da3baJj55ptvdMsttygpKUnx8fG65pprtHLlylC3FXbef/999enTRzExMUpKStJtt90W6pbCUm1trS6//HK5XC4VFRWFup2wsX37do0bN07du3dXTEyMLrjgAk2fPt363z/z589X9+7ddeaZZyotLU2ffPJJqFsKGzk5ObrqqqsUFxen5ORkjRw5Ul9//XWo28IpIgC2oLq6Ot1+++267777jri9oaFBw4cPV3V1tdasWaPc3FwtWbJEU6ZMaeVOw8vw4cN18OBBffTRRyosLNTll1+ujIyMoz6j2UZLlixRZmam/vCHP+h///d/9emnn2rUqFGhbissPfTQQ1Y/7ulovvrqKzU2NuqFF17Q1q1b9eyzz2rBggV65JFHQt1ayLzxxhuaNGmSHn30UW3atEkDBgzQTTfdpJKSklC3FhZWrVqlBx54QOvWrVN+fr4OHjyo9PR0VVdXh7o1nAqDFrdo0SLj8XiarF+2bJmJiIgwpaWlzrrFixcbt9tt/H5/K3YYPnbv3m0kmdWrVzvrAoGAkWQ+/PDDEHYWPurr602nTp3Mf/3Xf4W6lbC3bNkyc/HFF5utW7caSWbTpk2hbimszZo1y3Tv3j3UbYTM1VdfbcaPHx+07uKLLzYPP/xwiDoKbxUVFUaSWbVqVahbwSlgBDCECgoKlJqaGjQ6MXToUNXW1qqwsDCEnYXOOeeco0suuUSvvvqqqqurdfDgQb3wwgtKSUlRWlpaqNsLCxs3blRpaakiIiJ0xRVXqGPHjrrpppu0devWULcWVnbt2qWsrCz95S9/0VlnnRXqdtoEv9+vxMTEULcREnV1dSosLFR6enrQ+vT0dK1duzZEXYU3v98vSdb+nWnrCIAhVF5erpSUlKB1CQkJio6OtvZ0p8vlUn5+vjZt2qS4uDideeaZevbZZ5WXl3fU6yht8/3330uSZsyYoT/+8Y967733lJCQoIEDB2rv3r0h7i48GGM0duxYjR8/Xr179w51O23Cd999p7lz51r7KMyffvpJDQ0NTX4np6SkWPv7+FiMMZo8ebKuvfZapaamhrodnAIC4EmaMWOGXC7XMZfPP//8hN/P5XI1WWeMOeL6tuxEj5sxRvfff7+Sk5P1ySefaP369brllluUkZGhsrKyUO9GizrRY9TY2ChJevTRR/W73/1OaWlpWrRokVwul958880Q70XLOtFjNHfuXAUCAU2bNi3ULbe6U/kdtXPnTg0bNky333677rnnnhB1Hh4O/93bHn8fN4fs7Gxt3rxZixcvDnUrOEVWPwv4VGRnZ+vOO+88Zk23bt1O6L28Xq8+++yzoHWVlZWqr69v8q/Qtu5Ej9tHH32k9957T5WVlYqPj5f0y6y8/Px8vfLKK3r44Ydbo92QONFjVFVVJUnq2bOns97tduv8889v9xern+gxevzxx7Vu3bomzyrt3bu3Ro8erVdeeaUl2wypk/0dtXPnTl1//fXq16+fXnzxxRbuLnwlJSXpjDPOaDLaV1FR0e5+H5+uCRMm6N1339Xq1at13nnnhbodnCIC4ElKSkpSUlJSs7xXv3799MQTT6isrEwdO3aUJC1fvlxut7vdXe92osft559/liRFRAQPTkdERDgjX+3ViR6jtLQ0ud1uff3117r22mslSfX19dq+fbu6du3a0m2G1Ikeo//8z//U448/7vy8c+dODR06VG+88Yb69OnTki2G3Mn8jiotLdX111/vjCIf/v+dTaKjo5WWlqb8/Hzdeuutzvr8/HzdcsstIewsfBhjNGHCBC1dulQff/yxunfvHuqWcBoIgC2opKREe/fuVUlJiRoaGpx7kPXo0UMdOnRQenq6evbsqczMTD399NPau3evpk6dqqysLGf0yzb9+vVTQkKCxowZoz/96U+KiYnRn//8Z23btk3Dhw8PdXthIT4+XuPHj9f06dPVuXNnde3aVU8//bQk6fbbbw9xd+GhS5cuQT936NBBknTBBRcwYvE3O3fu1KBBg9SlSxfNnj1bu3fvdrZ5vd4QdhY6kydPVmZmpnr37u2MiJaUlFh7XeThHnjgAb3++ut65513FBcX54yWejwexcTEhLg7nLRQTkFu78aMGWMkNVlWrlzp1Pzwww9m+PDhJiYmxiQmJprs7GxTU1MTuqbDwIYNG0x6erpJTEw0cXFxpm/fvmbZsmWhbius1NXVmSlTppjk5GQTFxdnBg8ebIqLi0PdVtjatm0bt4E5zKJFi474+8n2r4XnnnvOdO3a1URHR5srr7ySW5z8ytH+vixatCjUreEUuIwxprVDJwAAAELH3gs+AAAALEUABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEcMoGDRqkSZMmNfv77tmzR8nJydq+fXuzvu/f//3fa86cOcetGzRokFwul1wul/MM70Pr28P+jh071tm/t99+u1k/E0DbQAAEEHZycnI0YsQIdevWLWj92rVr5XK5NGzYsKO+9rrrrnPCza+X0aNH609/+pOeeOIJBQKB4/aQlZWlsrIypaamnu7uHNfR9vdoRowYocGDBx9xW0FBgVwulzZu3HjU/f2P//gPlZWVnW7bANowAiCAsHLgwAG99NJLuueee5psW7hwoe666y6tXLlSJSUlTbYbY1RUVKTZs2errKwsaHnhhRd06aWXqlu3bnrttdeO28dZZ50lr9eryMjIZtmvoznW/h7NuHHj9NFHH+mHH35osm3hwoW6/PLLdeWVVx51fz0ej7xe72n3DqDtIgACaBa1tbWaOHGikpOTdeaZZ+raa6/Vhg0bgmqqqqo0evRoxcbGqmPHjnr22WebnFb9n//5H0VGRqpfv35Br62urtYbb7yhSZMm6YYbbtDLL7/cpIdvv/1WVVVVuu666+T1eoOWDh06SJJuvvlmLV68OKz31xijWbNm6fzzz1dMTIwuu+wy/fWvf3W2Z2RkKDk5uckx+Pnnn/XGG29o3Lhxzrrm2l8A7QsBEECzeOihh7RkyRK98sor2rhxo3r06KGhQ4dq7969Ts3kyZP16aef6t1331V+fr4++eQTbdy4Meh9Vq9erd69ezd5/zfeeENer1dXX321Ro8erUWLFskYE1RTWFioyMhIXXrppUft8+qrr9b69etVW1sbtvv7xz/+UYsWLdLzzz+vrVu36sEHH9Tdd9+tVatWSZIiIyP1D//wD3r55ZeDjsGbb76puro6jR49utn3F0D7QgAEcNqqq6v1/PPP6+mnn9ZNN92knj176s9//rNiYmL00ksvSfplNOyVV17R7NmzdeONNyo1NVWLFi1SQ0ND0Htt375dPp+vyWe89NJLTrAZOXKkKioqtGLFiqCajRs3qqGhQeecc446dOjgLFlZWU5Np06dVFtbq/Ly8rDc3+rqas2ZM0cLFy7U0KFDdf7552vs2LG6++679cILLzh1//iP/6jt27fr448/dtYtXLhQt912mxISEpp1fwG0Py17cQsAK3z33Xeqr6/XNddc46yLiorS1VdfrS+//FKS9P3336u+vl5XX321U+PxeHTRRRcFvdeBAwd05plnBq37+uuvtXbtWi1atEiS1KFDB91yyy1auHBh0GSIwsJC3X777XriiSeCXv/rQBQTEyPpl9Ol4bi/X3zxhWpqajRkyJCgurq6Ol1xxRXOzxdffLH69++vhQsX6vrrr9d3332nTz75RMuXLw96XXPsL4D2hwAI4LQdOg3pcrmarD+07lg1v5aUlKTKysqgdS+99JKuuuoqXXjhhc660aNH63e/+50qKyudgLdp0ybNnDlTPXr0OGqvh07RnnvuuSe8f4dryf1tbGyUJL3//vvq1KlTUK3b7Q76edy4ccrOztZzzz2nRYsWqWvXrrrxxhuDappjfwG0P5wCBnDaevTooejoaK1Zs8ZZV19fr88//1yXXHKJJOmCCy5QVFSU1q9f79QEAgF9++23Qe91xRVX6IsvvnB+PnjwoF599VWNGjUqqG7o0KGKi4tzZrh+//332rdvX9Ao2ZEUFxfrvPPOU1JS0qntrFp2f3v27Cm3262SkhL16NEjaOncuXPQa3//+9/rjDPO0Ouvv65XXnlFf/jDH5oEzubYXwDtDyOAAE5bbGys7rvvPv3Lv/yLEhMT1aVLF82aNUs///yzMyM1Li5OY8aMcWqSk5M1ffp0RUREBIWWoUOHatq0ac7I3nvvvaddu3YpNTVVxcXFQZ87YMAAvfTSS8rOzlZhYaEkKSUlpcn1bsnJyYqI+OXfu5988onS09PDdn/j4uI0depUPfjgg2psbNS1116rQCCgtWvXqkOHDhozZozz2g4dOuiOO+7QI488Ir/fr7FjxzbptTn2F0D7QwAE0CyeeuopNTY2KjMzU1VVVerdu7c++OCDoOvv5syZo/HjxysjI0Px8fF66KGHtGPHjqBr4Hr16qXevXvrv//7v3Xvvfc6kyoOvybu1zZu3OjMrv31aWLpl2vzqqqq5Ha7VVNTo6VLl+qDDz4I2/2VpH/9139VcnKycnJy9P333+vss8/WlVdeqUceeaRJH+PGjdNLL72k9PR0denSJWhbc+4vgPbFZQ6/IAUAWkl1dbU6deqkZ555JujedcuWLdPUqVNVXFzsjNw1h+eee07vvPNOk4kShxs0aJAuv/xy/fu//3uzfbYUfvvrcrm0dOlSjRw5stk+E0DbwDWAAFrNpk2btHjxYn333XfauHGjc1uXW265Jajut7/9re69916VlpY26+dHRUVp7ty5J1Q7f/58dejQQVu2bDnlzwvX/R0/frxzY2wAdmIEEECr2bRpk+655x59/fXXio6OVlpamubMmaNevXqFurUgpaWlOnDggCSpS5cuio6OPqX3Cdf9raiocJ4P3LFjR8XGxoa0HwCtjwAIAABgGU4BAwAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBl/j9QtqC0JygnzQAAAABJRU5ErkJggg==",
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz4klEQVR4nO3de3QUZZ7/8U/HJG0MSZsYk6aRmzJe2OAtKBdFUCHgEBSddVQwC7OYFTWwCKxHdH4DsmpcRNxdENFZQZ2jxHUQ9SgbiYggEgQDWQhej4LBkBCE0B0iuZA8vz8cam3CnSTdyfN+nVPnkKpvd3+rxPSHp+qpchljjAAAAGCNiFA3AAAAgNZFAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAy0SGuoG2rLGxUTt37lRcXJxcLleo2wEAACfAGKOqqir5fD5FRNg5FkYAPA07d+5U586dQ90GAAA4BTt27NB5550X6jZCggB4GuLi4iT98hcoPj4+xN0AAIATEQgE1LlzZ+d73EYEwNNw6LRvfHw8ARAAgDbG5su37DzxDQAAYDECIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAoHm4XL8sCHsEQAAAAMsQAAEAACwTGeoGAABAO/Pr08DGhK4PHBUjgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAQCAABOD/f+a3MYAQQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACwTGeoGAABAG+VyhboDnCJGAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAAAtx+VitnAYIgACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWKbZAmBOTo6uuuoqxcXFKTk5WSNHjtTXX38dVGOM0YwZM+Tz+RQTE6NBgwZp69atQTW1tbWaMGGCkpKSFBsbq5tvvlk//vhjUE1lZaUyMzPl8Xjk8XiUmZmpffv2BdWUlJRoxIgRio2NVVJSkiZOnKi6urqgmi1btmjgwIGKiYlRp06dNHPmTBljmuuQAAAAhKVmC4CrVq3SAw88oHXr1ik/P18HDx5Uenq6qqurnZpZs2Zpzpw5mjdvnjZs2CCv16shQ4aoqqrKqZk0aZKWLl2q3NxcrVmzRvv371dGRoYaGhqcmlGjRqmoqEh5eXnKy8tTUVGRMjMzne0NDQ0aPny4qqurtWbNGuXm5mrJkiWaMmWKUxMIBDRkyBD5fD5t2LBBc+fO1ezZszVnzpzmOiQAAADhybSQiooKI8msWrXKGGNMY2Oj8Xq95qmnnnJqampqjMfjMQsWLDDGGLNv3z4TFRVlcnNznZrS0lITERFh8vLyjDHGfPHFF0aSWbdunVNTUFBgJJmvvvrKGGPMsmXLTEREhCktLXVqFi9ebNxut/H7/cYYY+bPn288Ho+pqalxanJycozP5zONjY0ntI9+v99Ict4TAACrSCe+hBG+v41psWsA/X6/JCkxMVGStG3bNpWXlys9Pd2pcbvdGjhwoNauXStJKiwsVH19fVCNz+dTamqqU1NQUCCPx6M+ffo4NX379pXH4wmqSU1Nlc/nc2qGDh2q2tpaFRYWOjUDBw6U2+0Oqtm5c6e2b9/enIcCAAAgrLRIADTGaPLkybr22muVmpoqSSovL5ckpaSkBNWmpKQ428rLyxUdHa2EhIRj1iQnJzf5zOTk5KCawz8nISFB0dHRx6w59POhmsPV1tYqEAgELQAAAG1NiwTA7Oxsbd68WYsXL26yzXXY42CMMU3WHe7wmiPVN0eN+dsEkKP1k5OT40w88Xg86ty58zH7BgAACEfNHgAnTJigd999VytXrtR5553nrPd6vZKajq5VVFQ4I29er1d1dXWqrKw8Zs2uXbuafO7u3buDag7/nMrKStXX1x+zpqKiQlLTUcpDpk2bJr/f7yw7duw4xpEAAKAdOvRsX57v26Y1WwA0xig7O1tvvfWWPvroI3Xv3j1oe/fu3eX1epWfn++sq6ur06pVq9S/f39JUlpamqKiooJqysrKVFxc7NT069dPfr9f69evd2o+++wz+f3+oJri4mKVlZU5NcuXL5fb7VZaWppTs3r16qBbwyxfvlw+n0/dunU74j663W7Fx8cHLQAAAG1Oc80mue+++4zH4zEff/yxKSsrc5aff/7ZqXnqqaeMx+Mxb731ltmyZYu56667TMeOHU0gEHBqxo8fb8477zzz4Ycfmo0bN5obbrjBXHbZZebgwYNOzbBhw8yll15qCgoKTEFBgenVq5fJyMhwth88eNCkpqaaG2+80WzcuNF8+OGH5rzzzjPZ2dlOzb59+0xKSoq56667zJYtW8xbb71l4uPjzezZs094n5lFBACwzsnM/GUWcNhqtv8iko64LFq0yKlpbGw006dPN16v17jdbnPdddeZLVu2BL3PgQMHTHZ2tklMTDQxMTEmIyPDlJSUBNXs2bPHjB492sTFxZm4uDgzevRoU1lZGVTzww8/mOHDh5uYmBiTmJhosrOzg275YowxmzdvNgMGDDBut9t4vV4zY8aME74FjDH8BQIAWIgA2C64jOHRF6cqEAjI4/HI7/dzOhgAYIdTvfYvjOIG3988CxgAAMA6BEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwTGSoGwAAAG2AyxXqDtCMGAEEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAA0PJcLm4lE0YIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAlokMdQMAACCM8fSOdokRQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDKRoW4AAACEIZcr1B2gBTECCAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAluFZwAAA4Bc8/9cajAACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFimWQPg6tWrNWLECPl8PrlcLr399ttB28eOHSuXyxW09O3bN6imtrZWEyZMUFJSkmJjY3XzzTfrxx9/DKqprKxUZmamPB6PPB6PMjMztW/fvqCakpISjRgxQrGxsUpKStLEiRNVV1cXVLNlyxYNHDhQMTEx6tSpk2bOnCljTLMdDwAAgHDUrAGwurpal112mebNm3fUmmHDhqmsrMxZli1bFrR90qRJWrp0qXJzc7VmzRrt379fGRkZamhocGpGjRqloqIi5eXlKS8vT0VFRcrMzHS2NzQ0aPjw4aqurtaaNWuUm5urJUuWaMqUKU5NIBDQkCFD5PP5tGHDBs2dO1ezZ8/WnDlzmvGIAADQBrhcPAbONqaFSDJLly4NWjdmzBhzyy23HPU1+/btM1FRUSY3N9dZV1paaiIiIkxeXp4xxpgvvvjCSDLr1q1zagoKCowk89VXXxljjFm2bJmJiIgwpaWlTs3ixYuN2+02fr/fGGPM/PnzjcfjMTU1NU5NTk6O8fl8prGx8YT20e/3G0nOewIA0CZJrb+EEN/fxrT6NYAff/yxkpOTdeGFFyorK0sVFRXOtsLCQtXX1ys9Pd1Z5/P5lJqaqrVr10qSCgoK5PF41KdPH6emb9++8ng8QTWpqany+XxOzdChQ1VbW6vCwkKnZuDAgXK73UE1O3fu1Pbt24/Ye21trQKBQNACAADQ1rRqALzpppv02muv6aOPPtIzzzyjDRs26IYbblBtba0kqby8XNHR0UpISAh6XUpKisrLy52a5OTkJu+dnJwcVJOSkhK0PSEhQdHR0cesOfTzoZrD5eTkONcdejwede7c+WQPAQAAQMhFtuaH3XHHHc6fU1NT1bt3b3Xt2lXvv/++brvttqO+zhgj16+uTXAd4TqF5qgxf5sAcqTXStK0adM0efJk5+dAIEAIBAAAbU5IbwPTsWNHde3aVd9++60kyev1qq6uTpWVlUF1FRUVzuic1+vVrl27mrzX7t27g2oOH8WrrKxUfX39MWsOnY4+fGTwELfbrfj4+KAFAACgrQlpANyzZ4927Nihjh07SpLS0tIUFRWl/Px8p6asrEzFxcXq37+/JKlfv37y+/1av369U/PZZ5/J7/cH1RQXF6usrMypWb58udxut9LS0pya1atXB90aZvny5fL5fOrWrVuL7TMAAECoNWsA3L9/v4qKilRUVCRJ2rZtm4qKilRSUqL9+/dr6tSpKigo0Pbt2/Xxxx9rxIgRSkpK0q233ipJ8ng8GjdunKZMmaIVK1Zo06ZNuvvuu9WrVy8NHjxYknTJJZdo2LBhysrK0rp167Ru3TplZWUpIyNDF110kSQpPT1dPXv2VGZmpjZt2qQVK1Zo6tSpysrKckbtRo0aJbfbrbFjx6q4uFhLly7Vk08+qcmTJx/1FDAAAEC70JxTileuXGkkNVnGjBljfv75Z5Oenm7OPfdcExUVZbp06WLGjBljSkpKgt7jwIEDJjs72yQmJpqYmBiTkZHRpGbPnj1m9OjRJi4uzsTFxZnRo0ebysrKoJoffvjBDB8+3MTExJjExESTnZ0ddMsXY4zZvHmzGTBggHG73cbr9ZoZM2ac8C1gjGEaOQCgneA2MNZxGcOjL05VIBCQx+OR3+/nekAAQNsTyjNeIYwffH/zLGAAAADrEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsExkqBsAAACtjCdeWY8RQAAAAMsQAAEAACxDAAQAALAMARAAAMAyTAIBAMAGTPzArzACCAAAYBkCIAAAgGUIgAAAAJbhGkAAANozrv3DETACCAAAYBkCIAAAgGUIgAAAAJbhGkAAANojrv3DMTACCAAAYBkCIAAAgGUiNm7cqC1btjgr3nnnHY0cOVKPPPKI6urqQtgaAAAAWkLEvffeq2+++UaS9P333+vOO+/UWWedpTfffFMPPfRQiNsDAAAnxeXi+j8cV8Q333yjyy+/XJL05ptv6rrrrtPrr7+ul19+WUuWLAltdwAAAGh2EcYYNTY2SpI+/PBD/fa3v5Ukde7cWT/99FMoewMAAEALiOjdu7cef/xx/eUvf9GqVas0fPhwSdK2bduUkpIS4vYAAADQ3CKeffZZbdy4UdnZ2Xr00UfVo0cPSdJf//pX9e/fP8TtAQAAoLm5jDHmSBtqamoUGRmpyEjuFX00gUBAHo9Hfr9f8fHxoW4HAGCrtjjp48jxo1Xw/S1FnH/++dqzZ0+TDTU1NbrwwgtD0BIAAABaUsT27dvV0NDQZENtba1+/PHHELQEAACAlhQpSR988IE8Ho+zsqGhQStWrFD37t1D1hgAAABaRqTL5dKYMWOCVkZFRalbt2565plnQtQWAAAAWkpkly5dtGHDBiUlJYW6FwAAcDLa4uQPhIXIbdu2hboHAAAAtKJISVqxYoVWrFihiooK56kghyxcuDAkjQEAAKBlRD722GOaOXOmevfurY4dO8rFcDIAAEC7FrlgwQK9/PLLyszMDHUvAAAAaAURdXV1PPINAADAIhH33HOPXn/99VD3AQAATpTLxQxgnJbImpoavfjii/rwww916aWXKioqKqhgzpw5IWoNAAC0W78OsCF8LrCtIjdv3qzLL79cklRcXBy0kQkhAAAA7U/kypUrQ90DAAAAWlFEqBsAAABA64q8/vrrj3mq96OPPmrFdgAAANDSIg9d/3dIfX29ioqKVFxcrDFjxoSmKwAA0BTX5qOZRD777LNH3DBjxgzt37+/ldsBAABASzvqNYB33303zwEGAABoh44aAAsKCnTmmWe2Zi8AAABoBZG33XZb0ApjjMrKyvT555/r//2//xeitgAAANBSIj0eT9CKiIgIXXTRRZo5c6bS09ND1BYAAABaSuSiRYtC3QMAAABaUeShPxQWFurLL7+Uy+VSz549dcUVV4SyLwAAIHHrF7SIyIqKCt155536+OOPdfbZZ8sYI7/fr+uvv165ubk699xzQ90jAAAAmlHEhAkTFAgEtHXrVu3du1eVlZUqLi5WIBDQxIkTQ90fAAAAmpkrPj7efPjhh7rqqquCNqxfv17p6enat29faDprAwKBgDwej/x+v+Lj40PdDgCgPbLhFLAxrfpxfH9LEY2NjYqKimqyISoqSo2NjSFoCQAAAC0p4oYbbtA///M/a+fOnc7K0tJSPfjgg7rxxhtD2BoAAABaQsS8efNUVVWlbt266YILLlCPHj3UvXt3VVVVae7cuaHuDwAAO7lcdpz+RUhEdu7cWRs3blR+fr6++uorGWPUs2dPDR48ONS9AQAAoAW4/H6/OfwCSL/fr/79+2vBggUaMGBAiFoLf1xECgBoMTaN/jEJpNVFHGnHPR6P7r33Xs2ZMycELQEAAKAlRRxtQ3p6ugoLC1uzFwAAALSCowbAyMhI7d69uzV7AQAAQCs4agDcvHmzOnbs2Jq9AAAAoBVE1NTUNFl54MABTZ8+XRkZGSFoCQAAi3H7F7SCiAsvvFCzZs3SO++8o3fffVf/9m//posuukh79+7Vo48+elJvtnr1ao0YMUI+n08ul0tvv/120HZjjGbMmCGfz6eYmBgNGjRIW7duDaqpra3VhAkTlJSUpNjYWN1888368ccfg2oqKyuVmZkpj8cjj8ejzMzMJo+sKykp0YgRIxQbG6ukpCRNnDhRdXV1QTVbtmzRwIEDFRMTo06dOmnmzJkyrTwTCQAAoLVFpKamatq0abr11ls1cuRIPfLII0pNTdWnn36qlJSUk3qz6upqXXbZZZo3b94Rt8+aNUtz5szRvHnztGHDBnm9Xg0ZMkRVVVVOzaRJk7R06VLl5uZqzZo12r9/vzIyMtTQ0ODUjBo1SkVFRcrLy1NeXp6KioqUmZnpbG9oaNDw4cNVXV2tNWvWKDc3V0uWLNGUKVOcmkAgoCFDhsjn82nDhg2aO3euZs+ezcxnAADQ/hljzN69e8369evNZ599Zvbu3WuagySzdOlS5+fGxkbj9XrNU0895ayrqakxHo/HLFiwwBhjzL59+0xUVJTJzc11akpLS01ERITJy8szxhjzxRdfGElm3bp1Tk1BQYGRZL766itjjDHLli0zERERprS01KlZvHixcbvdxu/3G2OMmT9/vvF4PKampsapycnJMT6fzzQ2Np7QPvr9fiPJeU8AAE7bL3fFs2tpZXx/GxMhSQkJCbrqqqt09dVXKyEhoUWC5rZt21ReXq709HRnndvt1sCBA7V27VpJUmFhoerr64NqfD6fUlNTnZqCggJ5PB716dPHqenbt688Hk9QTWpqqnw+n1MzdOhQ1dbWOre2KSgo0MCBA+V2u4Nqdu7cqe3btzf/AQAAAAgTR50F3NzKy8slqclp5ZSUFGdbeXm5oqOjm4TQw2uSk5ObvH9ycnJQzeGfk5CQoOjo6GPWHPr5UM3hamtrFQgEghYAAIC2ptUC4CGuw2Y2GWOarDvc4TVHqm+OGvO3CSBH6ycnJ8eZeOLxeNS5c+dj9g0AABCOWi0Aer1eSU1H1yoqKpyRN6/Xq7q6OlVWVh6zZteuXU3ef/fu3UE1h39OZWWl6uvrj1lTUVEhqeko5SHTpk2T3+93lh07dhx/xwEAOBHc/gWtqNUCYPfu3eX1epWfn++sq6ur06pVq9S/f39JUlpamqKiooJqysrKVFxc7NT069dPfr9f69evd2o+++wz+f3+oJri4mKVlZU5NcuXL5fb7VZaWppTs3r16qBbwyxfvlw+n0/dunU74j643W7Fx8cHLQAAAG1NswbA/fv3q6ioSEVFRZJ+mfhRVFSkkpISuVwuTZo0SU8++aSWLl2q4uJijR07VmeddZZGjRolSfJ4PBo3bpymTJmiFStWaNOmTbr77rvVq1cvDR48WJJ0ySWXaNiwYcrKytK6deu0bt06ZWVlKSMjQxdddJGkX55j3LNnT2VmZmrTpk1asWKFpk6dqqysLCe0jRo1Sm63W2PHjlVxcbGWLl2qJ598UpMnTz7uKWkAAIA2rTmnFK9cudJIarKMGTPGGPPLrWCmT59uvF6vcbvd5rrrrjNbtmwJeo8DBw6Y7Oxsk5iYaGJiYkxGRoYpKSkJqtmzZ48ZPXq0iYuLM3FxcWb06NGmsrIyqOaHH34ww4cPNzExMSYxMdFkZ2cH3fLFGGM2b95sBgwYYNxut/F6vWbGjBknfAsYY5hGDgBoRqG+FQu3gbGKy5i/zXzASQsEAvJ4PPL7/ZwOBgCcHpvPPrVyFOH7OwSzgAEAABBakaFuAAAAa9k86oeQYgQQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMtwGxgAAFobt39BiDECCAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZZgEDANBamP2LMMEIIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAUMAEBLY/YvwgwjgAAAILRcLkJyKyMAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiG28AAANASmNWKMMYIIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlmAUMAEBzYvYv2gBGAAEAACxDAAQAALAMARAAAMAyBEAAAADLMAkEAIDmwOQPtCGMAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhlnAAACcDmb/og1iBBAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMswCBgDgZDHzF20cI4AAAACWIQACAABYhgAIAABgGQIgAACAZZgEAgDAiWLyB9oJRgABAAAsQwAEAACwDAEQAADAMlwDCADA8XDtH9oZRgABAAAsQwAEAACwDAEQAADAMgRAAAAAyzAJBACAo2HyB9opRgABAAAsQwAEAACwDAEQAADAMq0aAGfMmCGXyxW0eL1eZ7sxRjNmzJDP51NMTIwGDRqkrVu3Br1HbW2tJkyYoKSkJMXGxurmm2/Wjz/+GFRTWVmpzMxMeTweeTweZWZmat++fUE1JSUlGjFihGJjY5WUlKSJEyeqrq6uxfYdAAAgXLT6CODf/d3fqayszFm2bNnibJs1a5bmzJmjefPmacOGDfJ6vRoyZIiqqqqcmkmTJmnp0qXKzc3VmjVrtH//fmVkZKihocGpGTVqlIqKipSXl6e8vDwVFRUpMzPT2d7Q0KDhw4erurpaa9asUW5urpYsWaIpU6a0zkEAAIQvl+v/FqC9Mq1o+vTp5rLLLjvitsbGRuP1es1TTz3lrKupqTEej8csWLDAGGPMvn37TFRUlMnNzXVqSktLTUREhMnLyzPGGPPFF18YSWbdunVOTUFBgZFkvvrqK2OMMcuWLTMRERGmtLTUqVm8eLFxu93G7/ef8P74/X4j6aReAwAIcxJLqJZWwve3Ma0+Avjtt9/K5/Ope/fuuvPOO/X9999LkrZt26by8nKlp6c7tW63WwMHDtTatWslSYWFhaqvrw+q8fl8Sk1NdWoKCgrk8XjUp08fp6Zv377yeDxBNampqfL5fE7N0KFDVVtbq8LCwpbbeQAAgDDQqvcB7NOnj1599VVdeOGF2rVrlx5//HH1799fW7duVXl5uSQpJSUl6DUpKSn64YcfJEnl5eWKjo5WQkJCk5pDry8vL1dycnKTz05OTg6qOfxzEhISFB0d7dQcSW1trWpra52fA4HAie46AABA2GjVAHjTTTc5f+7Vq5f69eunCy64QK+88or69u0rSXIdds2FMabJusMdXnOk+lOpOVxOTo4ee+yxY/YCAGijuOYPFgnpbWBiY2PVq1cvffvtt85s4MNH4CoqKpzROq/Xq7q6OlVWVh6zZteuXU0+a/fu3UE1h39OZWWl6uvrm4wM/tq0adPk9/udZceOHSe5xwAAAKEX0gBYW1urL7/8Uh07dlT37t3l9XqVn5/vbK+rq9OqVavUv39/SVJaWpqioqKCasrKylRcXOzU9OvXT36/X+vXr3dqPvvsM/n9/qCa4uJilZWVOTXLly+X2+1WWlraUft1u92Kj48PWgAAANqaVj0FPHXqVI0YMUJdunRRRUWFHn/8cQUCAY0ZM0Yul0uTJk3Sk08+qd/85jf6zW9+oyeffFJnnXWWRo0aJUnyeDwaN26cpkyZonPOOUeJiYmaOnWqevXqpcGDB0uSLrnkEg0bNkxZWVl64YUXJEn/9E//pIyMDF100UWSpPT0dPXs2VOZmZl6+umntXfvXk2dOlVZWVmEOgAA0O61agD88ccfddddd+mnn37Sueeeq759+2rdunXq2rWrJOmhhx7SgQMHdP/996uyslJ9+vTR8uXLFRcX57zHs88+q8jISP3+97/XgQMHdOONN+rll1/WGWec4dS89tprmjhxojNb+Oabb9a8efOc7WeccYbef/993X///brmmmsUExOjUaNGafbs2a10JAAAAELHZYwxoW6irQoEAvJ4PPL7/YwcAkBbxeSP8NFKkYTvb54FDAAAYB0CIAAAgGVa9RpAAADCAqd9YTlGAAEAACxDAAQAALAMARAAAMAyXAMIALAH1/4BkhgBBAAAsA4BEAAAwDIEQAAAAMtwDSAAoP3j2j8gCCOAAAAAliEAAgAAWIZTwACA9otTv23Lof9exoS2DwswAggAAGAZAiAAAIBlCIAAAACW4RpAAED7wnV/wHExAggAAGAZAiAAAIBlOAUMAGgfOPULnDBGAAEAACxDAAQAALAMp4ABAG0bp36Bk8YIIAAAgGUYAQQAtE2M/AGnjBFAAAAAyxAAAQAALMMpYABA28FpX6BZMAIIAABgGQIgAACAZTgFDAAIf5z6BZoVI4AAAACWYQQQABC+GPkDWgQjgAAAAJZhBBAAEH4Y+QNaFCOAAAAAliEAAgAAWIZTwACA8MBpX6DVEAABAKFF8ANaHaeAAQAALMMIIAAgNBj5A0KGEUAAAADLMAIIAGhdjPwBIccIIAAAgGUYAQQAtDxG/YCwwgggAACAZRgBBAC0HEb+gLBEAAQAND+CHxDWCIAAgOZD8APaBK4BBAAAsAwjgACA08OoH9DmEAABAKeG4Ae0WQRAAMDJIfgBbR4BEABwYgh+QLtBAAQAHB2hD2iXCIAAgKYIfkC7RgAEAPwfgh9gBQIgAIDgB1iGAAgAtiL0AdYiAAKAbQh+gPUIgADQ3hH4AByGAAgA7RGhD8AxEAABoD0g8AE4CQRAAGhrCHto7w79HTcmtH20YwRAAAh3BD4AzYwACADhgJAHoBURAAGgtRDyAISJiFA3EGrz589X9+7ddeaZZyotLU2ffPJJqFsC0Ba5XMdfACBMWB0A33jjDU2aNEmPPvqoNm3apAEDBuimm25SSUlJqFsDEA5OJNQR7gC0QS5j7J1i06dPH1155ZV6/vnnnXWXXHKJRo4cqZycnOO+PhAIyOPxyO/3Kz4+viVbBXAqCGZA29ZCEYXvb4uvAayrq1NhYaEefvjhoPXp6elau3btEV9TW1ur2tpa52e/3y/pl79IAH7F4wl1BwDagxb6fj30vW3xGJi9AfCnn35SQ0ODUlJSgtanpKSovLz8iK/JycnRY4891mR9586dW6RHAACs1sL/mKyqqpLH0n+wWhsAD3EddorIGNNk3SHTpk3T5MmTnZ8bGxu1d+9enXPOOUd9jfTLvzQ6d+6sHTt2WDvUfDwco2Pj+Bwfx+j4OEbHxzE6vvZwjIwxqqqqks/nC3UrIWNtAExKStIZZ5zRZLSvoqKiyajgIW63W263O2jd2WeffcKfGR8f32b/Z2ktHKNj4/gcH8fo+DhGx8cxOr62foxsHfk7xNpZwNHR0UpLS1N+fn7Q+vz8fPXv3z9EXQEAALQ8a0cAJWny5MnKzMxU79691a9fP7344osqKSnR+PHjQ90aAABAi7E6AN5xxx3as2ePZs6cqbKyMqWmpmrZsmXq2rVrs36O2+3W9OnTm5w+xv/hGB0bx+f4OEbHxzE6Po7R8XGM2ger7wMIAABgI2uvAQQAALAVARAAAMAyBEAAAADLEAABAAAsQwBsQU888YT69++vs84666g3jC4pKdGIESMUGxurpKQkTZw4UXV1da3baJj55ptvdMsttygpKUnx8fG65pprtHLlylC3FXbef/999enTRzExMUpKStJtt90W6pbCUm1trS6//HK5XC4VFRWFup2wsX37do0bN07du3dXTEyMLrjgAk2fPt363z/z589X9+7ddeaZZyotLU2ffPJJqFsKGzk5ObrqqqsUFxen5ORkjRw5Ul9//XWo28IpIgC2oLq6Ot1+++267777jri9oaFBw4cPV3V1tdasWaPc3FwtWbJEU6ZMaeVOw8vw4cN18OBBffTRRyosLNTll1+ujIyMoz6j2UZLlixRZmam/vCHP+h///d/9emnn2rUqFGhbissPfTQQ1Y/7ulovvrqKzU2NuqFF17Q1q1b9eyzz2rBggV65JFHQt1ayLzxxhuaNGmSHn30UW3atEkDBgzQTTfdpJKSklC3FhZWrVqlBx54QOvWrVN+fr4OHjyo9PR0VVdXh7o1nAqDFrdo0SLj8XiarF+2bJmJiIgwpaWlzrrFixcbt9tt/H5/K3YYPnbv3m0kmdWrVzvrAoGAkWQ+/PDDEHYWPurr602nTp3Mf/3Xf4W6lbC3bNkyc/HFF5utW7caSWbTpk2hbimszZo1y3Tv3j3UbYTM1VdfbcaPHx+07uKLLzYPP/xwiDoKbxUVFUaSWbVqVahbwSlgBDCECgoKlJqaGjQ6MXToUNXW1qqwsDCEnYXOOeeco0suuUSvvvqqqqurdfDgQb3wwgtKSUlRWlpaqNsLCxs3blRpaakiIiJ0xRVXqGPHjrrpppu0devWULcWVnbt2qWsrCz95S9/0VlnnRXqdtoEv9+vxMTEULcREnV1dSosLFR6enrQ+vT0dK1duzZEXYU3v98vSdb+nWnrCIAhVF5erpSUlKB1CQkJio6OtvZ0p8vlUn5+vjZt2qS4uDideeaZevbZZ5WXl3fU6yht8/3330uSZsyYoT/+8Y967733lJCQoIEDB2rv3r0h7i48GGM0duxYjR8/Xr179w51O23Cd999p7lz51r7KMyffvpJDQ0NTX4np6SkWPv7+FiMMZo8ebKuvfZapaamhrodnAIC4EmaMWOGXC7XMZfPP//8hN/P5XI1WWeMOeL6tuxEj5sxRvfff7+Sk5P1ySefaP369brllluUkZGhsrKyUO9GizrRY9TY2ChJevTRR/W73/1OaWlpWrRokVwul958880Q70XLOtFjNHfuXAUCAU2bNi3ULbe6U/kdtXPnTg0bNky333677rnnnhB1Hh4O/93bHn8fN4fs7Gxt3rxZixcvDnUrOEVWPwv4VGRnZ+vOO+88Zk23bt1O6L28Xq8+++yzoHWVlZWqr69v8q/Qtu5Ej9tHH32k9957T5WVlYqPj5f0y6y8/Px8vfLKK3r44Ydbo92QONFjVFVVJUnq2bOns97tduv8889v9xern+gxevzxx7Vu3bomzyrt3bu3Ro8erVdeeaUl2wypk/0dtXPnTl1//fXq16+fXnzxxRbuLnwlJSXpjDPOaDLaV1FR0e5+H5+uCRMm6N1339Xq1at13nnnhbodnCIC4ElKSkpSUlJSs7xXv3799MQTT6isrEwdO3aUJC1fvlxut7vdXe92osft559/liRFRAQPTkdERDgjX+3ViR6jtLQ0ud1uff3117r22mslSfX19dq+fbu6du3a0m2G1Ikeo//8z//U448/7vy8c+dODR06VG+88Yb69OnTki2G3Mn8jiotLdX111/vjCIf/v+dTaKjo5WWlqb8/Hzdeuutzvr8/HzdcsstIewsfBhjNGHCBC1dulQff/yxunfvHuqWcBoIgC2opKREe/fuVUlJiRoaGpx7kPXo0UMdOnRQenq6evbsqczMTD399NPau3evpk6dqqysLGf0yzb9+vVTQkKCxowZoz/96U+KiYnRn//8Z23btk3Dhw8PdXthIT4+XuPHj9f06dPVuXNnde3aVU8//bQk6fbbbw9xd+GhS5cuQT936NBBknTBBRcwYvE3O3fu1KBBg9SlSxfNnj1bu3fvdrZ5vd4QdhY6kydPVmZmpnr37u2MiJaUlFh7XeThHnjgAb3++ut65513FBcX54yWejwexcTEhLg7nLRQTkFu78aMGWMkNVlWrlzp1Pzwww9m+PDhJiYmxiQmJprs7GxTU1MTuqbDwIYNG0x6erpJTEw0cXFxpm/fvmbZsmWhbius1NXVmSlTppjk5GQTFxdnBg8ebIqLi0PdVtjatm0bt4E5zKJFi474+8n2r4XnnnvOdO3a1URHR5srr7ySW5z8ytH+vixatCjUreEUuIwxprVDJwAAAELH3gs+AAAALEUABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEcMoGDRqkSZMmNfv77tmzR8nJydq+fXuzvu/f//3fa86cOcetGzRokFwul1wul/MM70Pr28P+jh071tm/t99+u1k/E0DbQAAEEHZycnI0YsQIdevWLWj92rVr5XK5NGzYsKO+9rrrrnPCza+X0aNH609/+pOeeOIJBQKB4/aQlZWlsrIypaamnu7uHNfR9vdoRowYocGDBx9xW0FBgVwulzZu3HjU/f2P//gPlZWVnW7bANowAiCAsHLgwAG99NJLuueee5psW7hwoe666y6tXLlSJSUlTbYbY1RUVKTZs2errKwsaHnhhRd06aWXqlu3bnrttdeO28dZZ50lr9eryMjIZtmvoznW/h7NuHHj9NFHH+mHH35osm3hwoW6/PLLdeWVVx51fz0ej7xe72n3DqDtIgACaBa1tbWaOHGikpOTdeaZZ+raa6/Vhg0bgmqqqqo0evRoxcbGqmPHjnr22WebnFb9n//5H0VGRqpfv35Br62urtYbb7yhSZMm6YYbbtDLL7/cpIdvv/1WVVVVuu666+T1eoOWDh06SJJuvvlmLV68OKz31xijWbNm6fzzz1dMTIwuu+wy/fWvf3W2Z2RkKDk5uckx+Pnnn/XGG29o3Lhxzrrm2l8A7QsBEECzeOihh7RkyRK98sor2rhxo3r06KGhQ4dq7969Ts3kyZP16aef6t1331V+fr4++eQTbdy4Meh9Vq9erd69ezd5/zfeeENer1dXX321Ro8erUWLFskYE1RTWFioyMhIXXrppUft8+qrr9b69etVW1sbtvv7xz/+UYsWLdLzzz+vrVu36sEHH9Tdd9+tVatWSZIiIyP1D//wD3r55ZeDjsGbb76puro6jR49utn3F0D7QgAEcNqqq6v1/PPP6+mnn9ZNN92knj176s9//rNiYmL00ksvSfplNOyVV17R7NmzdeONNyo1NVWLFi1SQ0ND0Htt375dPp+vyWe89NJLTrAZOXKkKioqtGLFiqCajRs3qqGhQeecc446dOjgLFlZWU5Np06dVFtbq/Ly8rDc3+rqas2ZM0cLFy7U0KFDdf7552vs2LG6++679cILLzh1//iP/6jt27fr448/dtYtXLhQt912mxISEpp1fwG0Py17cQsAK3z33Xeqr6/XNddc46yLiorS1VdfrS+//FKS9P3336u+vl5XX321U+PxeHTRRRcFvdeBAwd05plnBq37+uuvtXbtWi1atEiS1KFDB91yyy1auHBh0GSIwsJC3X777XriiSeCXv/rQBQTEyPpl9Ol4bi/X3zxhWpqajRkyJCgurq6Ol1xxRXOzxdffLH69++vhQsX6vrrr9d3332nTz75RMuXLw96XXPsL4D2hwAI4LQdOg3pcrmarD+07lg1v5aUlKTKysqgdS+99JKuuuoqXXjhhc660aNH63e/+50qKyudgLdp0ybNnDlTPXr0OGqvh07RnnvuuSe8f4dryf1tbGyUJL3//vvq1KlTUK3b7Q76edy4ccrOztZzzz2nRYsWqWvXrrrxxhuDappjfwG0P5wCBnDaevTooejoaK1Zs8ZZV19fr88//1yXXHKJJOmCCy5QVFSU1q9f79QEAgF9++23Qe91xRVX6IsvvnB+PnjwoF599VWNGjUqqG7o0KGKi4tzZrh+//332rdvX9Ao2ZEUFxfrvPPOU1JS0qntrFp2f3v27Cm3262SkhL16NEjaOncuXPQa3//+9/rjDPO0Ouvv65XXnlFf/jDH5oEzubYXwDtDyOAAE5bbGys7rvvPv3Lv/yLEhMT1aVLF82aNUs///yzMyM1Li5OY8aMcWqSk5M1ffp0RUREBIWWoUOHatq0ac7I3nvvvaddu3YpNTVVxcXFQZ87YMAAvfTSS8rOzlZhYaEkKSUlpcn1bsnJyYqI+OXfu5988onS09PDdn/j4uI0depUPfjgg2psbNS1116rQCCgtWvXqkOHDhozZozz2g4dOuiOO+7QI488Ir/fr7FjxzbptTn2F0D7QwAE0CyeeuopNTY2KjMzU1VVVerdu7c++OCDoOvv5syZo/HjxysjI0Px8fF66KGHtGPHjqBr4Hr16qXevXvrv//7v3Xvvfc6kyoOvybu1zZu3OjMrv31aWLpl2vzqqqq5Ha7VVNTo6VLl+qDDz4I2/2VpH/9139VcnKycnJy9P333+vss8/WlVdeqUceeaRJH+PGjdNLL72k9PR0denSJWhbc+4vgPbFZQ6/IAUAWkl1dbU6deqkZ555JujedcuWLdPUqVNVXFzsjNw1h+eee07vvPNOk4kShxs0aJAuv/xy/fu//3uzfbYUfvvrcrm0dOlSjRw5stk+E0DbwDWAAFrNpk2btHjxYn333XfauHGjc1uXW265Jajut7/9re69916VlpY26+dHRUVp7ty5J1Q7f/58dejQQVu2bDnlzwvX/R0/frxzY2wAdmIEEECr2bRpk+655x59/fXXio6OVlpamubMmaNevXqFurUgpaWlOnDggCSpS5cuio6OPqX3Cdf9raiocJ4P3LFjR8XGxoa0HwCtjwAIAABgGU4BAwAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBl/j9QtqC0JygnzQAAAABJRU5ErkJggg==' width=640.0/>\n",
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "from os import path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from qcnico.molecule import Molecule\n",
    "from qcnico import plt_utils\n",
    "\n",
    "def get_absdiffs(arr):\n",
    "    \"\"\"Compute the absolute value of differences between elements of a 1D array.\n",
    "    To avoid redundancy when histograming; only upper triangular (excluding the main diagonal) is returned.\"\"\"\n",
    "    return np.triu(np.abs(arr[:,None] - arr[None,:]),1)\n",
    "\n",
    "def isdegenerate(energies,tol=27.2114 * 1e-8, return_absdiffs=False):\n",
    "    \"\"\"Returns `True` if `energies` contains the same value more than once (up to precision `tol`).\n",
    "    The default value of `tol` is set to what it is because QCFFPI offers an energy of 1e-8 Ha, and \n",
    "    the energies we are working with have usually been converted to eV (1 Ha = 27.2114 eV).\"\"\"\n",
    "    absdiffs = get_absdiffs(energies)\n",
    "    Nnnz = np.sum(absdiffs > tol)\n",
    "    N = energies.size\n",
    "    isdegen =  Nnnz != (N**2 - N) // 2\n",
    "    if return_absdiffs:\n",
    "        return isdegen, absdiffs\n",
    "    else:\n",
    "        return isdegen\n",
    "\n",
    "\n",
    "datadir = path.expanduser('~/Desktop/simulation_outputs/qcffpi_data/MO_dynamics/300K_initplanar_norotate/')\n",
    "\n",
    "mac = Molecule('frame-80000', datadir)\n",
    "#pos, M = mac.fetch_pos_MOs()\n",
    "mac.fetch_energies(prefix='orb_energy_')\n",
    "ee = mac.energies\n",
    "\n",
    "isdgen, ediffs = isdegenerate(ee, return_absdiffs=True)\n",
    "print('Energy spectrum contains degeneracies = ', isdgen)\n",
    "ediffs = ediffs[ediffs > 0]\n",
    "plt_utils.histogram(np.log(ediffs), nbins=200, xlabel='log($\\Delta E$) [log(eV)]')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import njit\n",
    "\n",
    "@njit\n",
    "def log_distance(ei, ej, mu=0, T=300):\n",
    "    kB = 8.617e-5\n",
    "    return (np.abs(ei-mu) + np.abs(ej-mu) + np.abs(ej-ei)) / (2*kB*T) \n",
    "\n",
    "@njit\n",
    "def dArr_logE(e):\n",
    "    N = e.size\n",
    "    darr = np.zeros(N*(N-1)//2)\n",
    "\n",
    "    k = 0\n",
    "    \n",
    "    for i in range(N):\n",
    "        for j in range(i):\n",
    "            darr[k] = log_distance(e[i], e[j])\n",
    "            k+=1\n",
    "    \n",
    "    return darr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "79d53e2e3f3357e27f4d32bf088d4f1991bddd65b1f38f3a2694c1c60a1b6da0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}