{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "859f09ac-1e0f-40d1-87b3-3f3e16c00272",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "********************************************************************************\n",
      "WARNING: Imported VTK version (9.3) does not match the one used\n",
      "         to build the TVTK classes (9.2). This may cause problems.\n",
      "         Please rebuild TVTK.\n",
      "********************************************************************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%gui qt\n",
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "from mayavi import mlab\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from kifmm_py import (\n",
    "    KiFmm,\n",
    "    LaplaceKernel,\n",
    "    SingleNodeTree,\n",
    "    EvalType,\n",
    "    BlasFieldTranslation,\n",
    "    FftFieldTranslation,\n",
    "    read_stl_triangle_mesh_vertices,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bca79974-3942-4ea4-9325-5314dfb064eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "(sources, faces) = read_stl_triangle_mesh_vertices(\"battleship.STL\")\n",
    "x = sources[:, 0]\n",
    "y = sources[:, 1]\n",
    "z = sources[:, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a744f8cc-2d34-4d84-9aec-cabea34cfaa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.view(azimuth=40, elevation=70, distance=\"auto\", focalpoint=\"auto\")\n",
    "fig = mlab.figure(\n",
    "    size=(960, 1080), bgcolor=(1, 1, 1)\n",
    ")  # This sets the window size for rendering\n",
    "plot = mlab.triangular_mesh(\n",
    "    x, y, z, faces, color=(0.5, 0.5, 0.5), representation=\"points\", figure=fig\n",
    ")\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5ba2c3af-949d-44a5-b798-2e2ec22c2445",
   "metadata": {},
   "outputs": [],
   "source": [
    "dtype = np.float64\n",
    "\n",
    "# # Set FMM Parameters\n",
    "expansion_order = np.array([6])\n",
    "# n_vec = 1\n",
    "n_crit = 15\n",
    "n_sources = sources.shape[0]\n",
    "n_targets = sources.shape[0]\n",
    "\n",
    "sources = sources.flatten().astype(dtype)\n",
    "charges = np.random.rand(n_sources).astype(dtype)\n",
    "\n",
    "# Set Random charges\n",
    "# charges = np.random.rand(n_sources).astype(dtype)\n",
    "\n",
    "kernel = LaplaceKernel(dtype, EvalType.Value)\n",
    "\n",
    "field_translation = BlasFieldTranslation(\n",
    "    kernel,\n",
    "    svd_threshold=dtype(1e-16),\n",
    "    surface_diff=2,\n",
    "    random=True,\n",
    "    n_components=100,\n",
    "    n_oversamples=10,\n",
    ")\n",
    "\n",
    "tree = SingleNodeTree(sources, sources, charges, n_crit=n_crit, prune_empty=True)\n",
    "\n",
    "fmm = KiFmm(expansion_order, tree, field_translation, timed=True)\n",
    "\n",
    "fmm.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3e52e08-a2c0-473a-9c43-8d5a9fb47b12",
   "metadata": {},
   "source": [
    "# Plot Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f615beab-0210-48e1-8170-dfa360d8d5bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.view(azimuth=40, elevation=70, distance=\"auto\", focalpoint=\"auto\")\n",
    "fig = mlab.figure(\n",
    "    size=(960, 1080), bgcolor=(1, 1, 1)\n",
    ")  # This sets the window size for rendering\n",
    "solution = mlab.triangular_mesh(\n",
    "    x,\n",
    "    y,\n",
    "    z,\n",
    "    faces,\n",
    "    scalars=np.log(fmm.all_potentials_u[0][:, 0]),\n",
    "    representation=\"surface\",\n",
    ")\n",
    "\n",
    "colorbar = mlab.colorbar(solution, title=\"Log-scaled values\", orientation=\"vertical\")\n",
    "colorbar.title_text_property.font_size = 10\n",
    "colorbar.label_text_property.font_size = 8\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88525acc-aae5-4ccf-a069-0229429584ef",
   "metadata": {},
   "source": [
    "# Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "91893002-cd58-49e6-b421-0c4d3f6d7c1f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stats = defaultdict(int)\n",
    "\n",
    "for leaf in fmm.leaves_target_tree:\n",
    "    found = fmm.coordinates_target_tree(leaf)\n",
    "    n = found.shape[0]\n",
    "    stats[n] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6ae7f253-fc1b-4ee2-8d5c-ced227dba336",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAJICAYAAADxUwLTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABB6klEQVR4nO3deXwV5d3///fJRjZWgbDaIBDBFqqlYZG6UQHLjRYqAlbRUpG7rdy9laBW6hdFocX7vsUHVlzAKovYIqCI2MItVKkiixEFgbDEGknYIWwh20kyvz+4z/nlZL3OOuecvJ7/xMxcM/OZmSsxb851zTgsy7IEAAAAAAZi7C4AAAAAQOQgQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjMXZXUC4++KLL2RZluLj4+0uBQAAAAgKp9Mph8Oha665ptG2fALRCMuyFM4v67YsS+Xl5WFdY2OOHTtmdwmIENHQ3wFT9Hc0JfR3+3nzNy+fQDTC9clDnz59bK6kbsXFxcrJyVGPHj2UnJxsdzk++d3vfqf333/f7jIQAaKhvwOm6O9oSujv9vvqq6+M2/IJBAAAAABjBAgAAAAAxggQAAAAAIwRIAAAAAAYI0AAAAAAMEaAAAAAAGCMAAEAAADAGAECAAAAgDECBAAAAABjBAgAAAAAxggQAAAAAIwRIAAAAAAYI0AAAAAAMEaAAAAAAGCMAAEAAADAGAECAAAAgDECBAAAAABjBAgAAAAAxggQAAAAAIwRIAAAAAAYI0AAAAAAMEaAAAAAAGCMAAEAAADAGAECAAAAgDECBNAEjBo1yu4SAABAlCBAAE2A0+m0uwQAABAlCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADG4uwuoD7Z2dmaN2+eli5d6rF8w4YN2r17tzp27KgDBw5owoQJSk9Pd6/Pzc3VihUrdMUVV+jQoUPq37+/brjhhhBXDwAAAESnsAsQH330kT777DOdOHFChw8f9li3Z88eLVu2TK+//rok6cyZM5owYYKWL1+ulJQUFRcX68EHH9Rf/vIXNW/eXJZl6Z577lGbNm3Up08fO04HAAAAiCphN4Tpxhtv1MMPP6xBgwbVWrdw4UINGTLE/X3r1q2Vnp6utWvXSpJWr16tbt26qXnz5pIkh8OhoUOHatGiRSGpHQAAAIh2YRcgGrJt2zZ17NjRY1nHjh21devWetd36tTJvR4AAACAfyImQBQVFamwsFDJyckey5OTk3Xo0CFJUn5+vpKSkjzWJyUl6dSpUyouLg5ZrQAAAEC0Crs5EPUpKiqSJMXGxnosj42Nda+7cOGC4uLiaq13rasZPkxZlhW2AaSkpMTjaySqqqoK2+sbLaLlGkdDfwdM0d/RlNDf7WdZlhwOh1HbiAkQLpZl1fq++rK61vvL6XQqJyfH7/0EU15ent0l+MThcKi4uFj79u0LyL0KR64fRrvOLxqvcV393e7rDARLpP5+B3xBf7dXQkKCUbuICRCpqamSpIqKCo/lTqdTKSkp7jZ1rZfkbuOL+Ph49ejRw+ftg6mkpER5eXlKT0+vNXwrUiQnJ6tXr152lxHVouUaR0N/B0zR39GU0N/tl5uba9w2ogJE69at3cOVXIqKitS1a1dJUteuXetc36pVK3cA8YXD4fB5+FOoJCUlhX2N9YmJiYnY2iNFtF3jSO7vgLfo72hK6O/2MR2+JEXQJGpJGjhwoAoKCjyW5efnKzMzs8H1/fv3D1mNAAAAQDQL2wBRVVWlqqoqj2X33XefNmzY4B7jfPLkSRUUFGj06NGSpNtuu03ffvutTp8+7d7Hxo0bdf/994e2eAAAACBKhd0Qph07dmjTpk3atGmTjh8/rpkzZyotLU2/+tWv1KdPH/3yl7/Us88+q86dOysnJ0cvvPCCe3hSamqqXnjhBc2fP189e/bUoUOHdP/996tv3742nxUAAAAQHcIuQPTt21ff//739eCDD8rhcMiyLPdEaEkaNmyYhg0bVu/2PXv21IwZM0JRKgAAANDkhF2AqPkeB4fDYfxIKQAAAADBFbZzIAAAAACEHwIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMBYnN0F+MKyLK1atUolJSWqrKzUsWPHNHbsWF1xxRWSpOzsbG3cuFHp6enKzc3Vrbfeqr59+9pcNQAAABD5IjJALF26VNdff73S09MlSaWlpXrggQf00ksv6cyZM5o9e7ZWrFihuLg4lZWVady4cXr55ZfVoUMHewsHAAAAIlxEDmHasGGDunbt6v4+MTFRGRkZOnjwoJYsWaIBAwYoLu5SNmrWrJkGDBigN998065yAQAAgKgRkQGisrJSTz75pMrLy93LDh06pC5dumjbtm3q2LGjR/uOHTtq69atoS4TAAAAiDoRGSAmTpyot956S7fffrt27dqlV199VSNHjlTLli2Vn5+v5ORkj/bJyck6dOiQTdUCAAAA0SMi50DcfPPNeu6555SVlaU77rhDd9xxhyZOnChJunDhgmJjYz3ax8bGqqioyOfjWZal4uJiv2oOlpKSEo+vkaiqqipsr2+0iJZrHA39HTBFf0dTQn+3n2VZcjgcRm0jMkAcPXpUmzdv1vLlyzV//nytWLFCR48e1cKFCyVdugDVWZZVa5k3nE6ncnJy/Ko52PLy8uwuwScOh0PFxcXat2+fX/conLl+GO06v2i8xnX1d7uvMxAskfr7HfAF/d1eCQkJRu0iLkBYlqVp06bp+eef12WXXaZXXnlFK1as0MyZM7V27VqlpqaqoqLCYxun06mUlBSfjxkfH68ePXr4W3pQlJSUKC8vT+np6UpKSrK7HJ8kJyerV69edpcR1aLlGkdDfwdM0d/RlNDf7Zebm2vcNuICxMGDB9WiRQtddtll7mV33HGHTp48qS+//FJdu3atNVypqKhIXbp08fmYDoej1ryKcJOUlBT2NdYnJiYmYmuPFNF2jSO5vwPeor+jKaG/28d0+JIUgQGiqqrK4+lLLhkZGaqsrNTAgQNVUFDgsS4/P1/9+/cPVYkAAABA1Iq4pzBdeeWVKi8v15dffuleVlFRofXr1+uOO+7QXXfdpe3bt6u0tFTSpY/Etm/frnvuucemigEAAIDoEXGfQDgcDr344ot6+eWXtXHjRjVr1kwXL17Ub37zG/ebpp988knNnTtX3bp1U25urmbNmqVOnTrZXDkAAAAQ+SIuQEhS8+bN9fDDD9e7PjMzU5mZmSGsCAAAAGgaIm4IEwAAAAD7ECAAAAAAGCNAAAAAADBGgAAAAABgjAABAAAAwBgBAgAAAIAxAgQAAAAAYwQIAAAAAMYIEAAAAACMESAAAAAAGCNAAAAAADBGgAAAAABgjAABAAAAwBgBAgAAAIAxAgQAAAAAYwQIAAAAAMYIEAAAAACMESAAAAAAGCNAAAAAADBGgAAAAABgjAABAAAAwBgBAgAAAIAxAgQAAAAAYwQIAAAAAMYIEAAAAACMESAAAAAAGCNAAAAAADBGgAAAAABgjAABwG+jRo2yuwQAABAiBAgAfnM6nXaXAAAAQoQAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYCzO3x0cO3ZMHTp08Fh28OBBrV+/Xi1bttTVV1+tPn36+HsYAAAAAGHA708gnnzyyVrLevbsqSlTpmj06NE6duyYXn75ZX8PAwAAACAMBHUIU2pqqoYOHapTp04F8zAAAAAAQsTrIUxHjhzx+L6srExHjx6VZVm12paUlOirr77SgQMHfK8QAAAAQNjwOkBs2LBBhw8f1saNG3X48GFJ0pAhQ+psa1mWWrdurXnz5vlXJQAAAICw4HWAuOeeeyRJWVlZmjZtmk6cOKGpU6fWaudwOJScnKyePXsqISHB/0oBAAAA2M7npzAlJCTo97//vZ544gn1798/kDUBAAAACFN+TaJOS0vTSy+91Gi7r7/+2p/DAAAAAAgTfj+FyeFwNNrmueee8/cwAAAAAMKA3y+Sk6R169bp888/V0lJiaqqqjyeyFRVVaUtW7YE4jAAAAAAbOZ3gHj++ef10ksvqXXr1kpMTKy1vry8XMXFxf4eBgAAAEAY8DtA7N27Vx9++KE6dOhQb5uJEyf6exgAAAAAYcDvANGnT58Gw4Mk3X333f4eppby8nItWLBADodDKSkpio+P11133SVJys7O1saNG5Wenq7c3Fzdeuut6tu3b8BrAAAAAJoavwOEZVkqLy9v8F0PsbGx/h7GQ2VlpaZMmaJ7771XgwcPVkFBgX7605+qX79+at26tWbPnq0VK1YoLi5OZWVlGjdunF5++eVGgw4AAACAhvn9FKa7775bf/7zn3XgwIF62yxZssTfw3h48803FRcXp8GDB0uS2rZtq7vuukudO3fWkiVLNGDAAMXFXcpGzZo104ABA/Tmm28GtAYAAACgKfL7E4isrCw5nU4tWLBAktSmTRuP9eXl5Tp16pS/h/HwxhtvuN+ILUmJiYnut2Fv27ZNt956q0f7jh076m9/+1tAawAAAACaIr8DxP79+/Wd73xHt9xyS53ry8rK9OGHH/p7GLfjx48rLy9PSUlJWrRokWJjY3Xw4EFNnjxZXbp0UX5+vpKTkz22SU5O1qFDhwJWAwAAANBU+R0g2rZtqyVLljQ4z+HXv/61v4dxO3bsmCTpyy+/1FNPPSVJ2rVrlyZPnqw1a9bowoULtWqJjY1VUVGRz8e0LCtsH0VbUlLi8TUSVVVVhe31jRbBvsahuofR0N8BU/R3NCX0d/tZlmX0gmgpAAFizpw5jU6SnjZtmr+HcausrJQk9e7d272sb9++Onr0qHuYUvUX2bm+r7nMG06nUzk5OT5vHwp5eXl2l+ATh8Oh4uJi7du3z697FM5cP4x2nV+wr7Ed97Cu/m73dTYRCTUi/ETq73fAF/R3ezX0UKTq/A4QvXr1anB9ZWWlli1bphkzZvh7KElSamqqJCktLc1jecuWLbV3716lpqaqoqLCY53T6VRKSorPx4yPj1ePHj183j6YSkpKlJeXp/T0dCUlJdldjk+Sk5Mb7UfwT7CvcajuYTT0d8AU/R1NCf3dfrm5ucZt/Q4Qx48fd38qUJe9e/fq008/9fcwbpdffrliY2PldDo9lldVVSkmJkZdu3atNVypqKhIXbp08fmYDoej1ryKcJOUlBT2NdYnJiYmYmuPFMG+xqG+h5Hc3wFv0d/RlNDf7WM6fEkKQIAYN26cjh8/3mCbFi1a+HsYt8TERF1zzTU6evSoe5llWTp37py+973vyeFwqKCgwGOb/Px89e/fP2A1AAAAAE2V3wGiVatWGjZsmMecBEm6ePGicnJyVFJSopEjR/p7GA+TJ0/WwoULde+998rhcOjTTz9Venq6hg8frquvvlr333+/SktLlZiYqJKSEm3fvl2vvfZaQGsAAAAAmqKABIhp06bVO+ni66+/1vbt2/09jIcbbrhBhYWFmj59utq1a6cTJ05owYIFio2NVadOnfTkk09q7ty56tatm3JzczVr1ix16tQpoDUAAAAATZHfAeL1119vcMxU9+7d9d577/l7mFpGjx6t0aNH17kuMzNTmZmZAT8mAAAA0NTF+LsDkwkXPLIQAAAAiA5+fwLxxRdf1HpsqotlWcrLy9P+/fv9PQwAAACAMOB3gMjKytLRo0fr/ZShc+fO+tOf/uTvYQAAAACEgYBMon700UfVu3fvWm+kTkxM1GWXXebvIQAAAACECb8DRGZmpoYNG+bVyycAAAAARCa/A8Rjjz0WiDoAAAAARAC/A4TL7t27tWzZMu3du1eVlZXq2bOnxo0bp4EDBwbqEAAAAABs5vdjXCVpwYIFGjt2rP73f/9XFy5cUHFxsT788ENNnDhRzz33XCAOAQAAACAM+P0JxPr167Vjxw6tWbNGPXr08Fi3e/du/fGPf9TatWs1cuRIfw8FAAAAwGZ+fwLx8ccfa/78+bXCgyR973vf06uvvqr169f7exgAAAAAYcDvANG+fftaj2+tLikpSWlpaf4eBgAAAEAY8DtAlJSUNNqmtLTU38MAAAAACAN+Bwin06kvvvii3vWbNm1STExA5moDAAAAsJnfk6h/85vf6Oc//7m++93vqk+fPmrVqpUk6dSpU8rOzlZOTo5WrFjh72EAAAAAhAG/A0SbNm20ePFizZgxQ3PmzPFYd/XVV2vp0qVq3769v4cBAAAAEAYC8iK5tLQ0vfLKKyooKNDevXtVUVGhnj17qmfPnoHYPQAAAIAw4VWAmDJlinJyctzfOxwODR48WDNnzpQkdenSRbt27dLmzZv1zTffKC0tTS1atAhsxQAAAABs41WAmDt3rp566imtW7dO//mf/6lbbrlF7dq182gzYsQIjRgxQidPntQLL7ygkSNHqm/fvgEtGgAAAIA9vAoQlmWptLRU77//fqPvdmjXrp2mT5+uJ554QmlpabwLAgAAAIgCXj1fdenSpZo8ebJXYeCBBx7QwoULvS4MQNMzatQou0sAAACN8CpAHDhwQBkZGV4doH379jp9+rRX2wBompxOp90lAACARngVIIqLi306SHl5uU/bAQAAAAgvXgWIs2fP+nSQ8+fP+7QdAAAAgPDiVYBITEzU3r17vTrAnj171KxZM6+2AQAAABCevAoQI0eO1FNPPWU8Trm8vFxPP/20brvtNp+KAwAAABBevAoQP/3pT+V0OjV58mQdPny4wbZHjhzR5MmTVVlZSYAAAAAAooRX74FwOByaP3++xo8fr+HDh+tHP/qRBgwYoK5duyolJUUXL15UQUGBtm7dqk8++UTt2rXTX//612DVDgAAACDEvAoQktShQwetXLlSWVlZ+uijj7Rp06ZabSzLUv/+/fXss8/WelM1AAAAgMjldYCQpLZt22rx4sXasmWL/v73vysnJ0fnz59XixYt1Lt3b91yyy269tprA10rAAAAAJv5FCBcBg0apEGDBgWqFgAAAABhzqtJ1AAAAACaNgIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjMXZXUAgPPLII8rKylJaWpokacOGDdq9e7c6duyoAwcOaMKECUpPT7e3SAAAACAKRHyA2LRpk95991399re/lSTt2bNHy5Yt0+uvvy5JOnPmjCZMmKDly5crJSXFzlIBAACAiBfRQ5guXryogwcPeixbuHChhgwZ4v6+devWSk9P19q1a0NdHgAAABB1IjpALF++XHfeeafHsm3btqljx44eyzp27KitW7eGsjQAAAAgKkVsgNi1a5d69OjhMSypqKhIhYWFSk5O9mibnJysQ4cOhbpEAAAAIOpE5ByIiooKbd26VZMnT/ZYXlRUJEmKjY31WB4bG+te5wvLslRcXOzz9sFUUlLi8TUSVVVVhe31jRbBvsaB2n9j+4mG/g6Yor+jKaG/28+yLDkcDqO2ERkg3n77bY0ZM6be9ZZl1fq+5jJvOJ1O5eTk+Lx9KOTl5dldgk8cDoeKi4u1b98+v+5ROHP9MNp1fsG+xoHavzf7qau/e3ud7bgvdveFUGgK5xhqkfr7HfAF/d1eCQkJRu0iLkB8++23atGihdq0aVNrXWpqqqRLn1BU53Q6/XoCU3x8vHr06OHz9sFUUlKivLw8paenKykpye5yfJKcnKxevXrZXUZUC/Y1DtT+G9tPNPR3wBT9HU0J/d1+ubm5xm0jLkB88sknysvL0+effy5JKi8vlyTNnz9fGRkZat26da3hSkVFReratavPx3Q4HLXmVYSbpKSksK+xPjExMRFbe6QI9jUO1P5N9xPJ/R3wFv0dTQn93T6mw5ekCAwQd911l8f3BQUF+utf/6oHHnhAXbp00c6dO1VQUODRJj8/XzfeeGMIqwQAAACiU8Q+hcnFNc7W9fW+++7Thg0b3N+fPHlSBQUFGj16tG01AgAAANEi4j6BqG7VqlXasGGDJGnOnDm6+eabNXr0aP3yl7/Us88+q86dOysnJ0cvvPCCe34EAAAAAN9FdIC49dZbNXr0aDkcDlmW5Z48PWzYMA0bNszm6gAAAIDoE9EBovqjphwOh/GjpwAAAAD4JuLnQAAAAAAIHQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECES1UaNG2V0CEHb4uQAA+IMAgajmdDrtLgEIO/xcAAD8QYAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAANGrUKLtLAABECAIEAEBOp9PuEgAAEYIAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAIkFGjRtldAgAAQUeAAIAAcTqddpcAAEDQESAAAAAAGCNAAAAAADBGgAAAAABgjAABAAAAwBgBAgAAAIAxAgQAAAAAYwQIAAAAAMYIEAAAAACMESAAAAAAGCNAAAAAADBGgAAAAABgjAABNFGjRo2yuwQEGfcYABAMBAigiXI6nXaXgCDjHgMAgoEAAQAAAMAYAQIAAACAMQIEEOEY546mgH4OAOGDAAFEOMa5oymgnwNA+CBAAAAAADBGgAAAAABgjAABAAAAwBgBAhATNBEZ6KcAgHBAgADEBE1EBvopACAcxNldgC9KS0u1aNEilZWV6dChQ4qNjdWjjz6qyy67TJKUnZ2tjRs3Kj09Xbm5ubr11lvVt29fm6sGAAAAIl9EBojXXntNd955p1q3bi1JmjFjhiZOnKiVK1fqzJkzmj17tlasWKG4uDiVlZVp3Lhxevnll9WhQwebKwcAAAAiW8QNYSorK9PixYu1efNm97KJEydq//792rRpk5YsWaIBAwYoLu5SNmrWrJkGDBigN998066SEWCMA7cH1x0AAEgRGCAcDoeSk5N16tQp97JOnTpJkr799ltt27ZNHTt29NimY8eO2rp1a0jrRPAwDtweXHcAACBF4BCmhIQEffjhhx7LDh8+LEm6/PLLlZ+fr+TkZI/1ycnJOnToUMhqBAAAAKJVxAWIurz//vvq3LmzbrzxRl24cEGxsbEe62NjY1VUVOTz/i3LUnFxsb9lBkVJSYnH10hUVVXl1fX1pr1pW29rCCcmtdfVJpDXPVDXr7H9hHt/D3Y/CtQ9C0R/CLVwry8Ywr2/A4FEf7efZVlyOBxGbSM+QBw/flwrV67USy+9pISEBEmXLkB1lmXVWuYNp9OpnJwcv+r0letGNlZ/Xl5eCKoJPIfDoeLiYu3bt8/oHnnT3rSttzV4y/Qe+rrvxmqvq00gr3ugrp83+6mrv3t7nb1pb9LW9F54U6O3+zdpH4j+0NhxpcD292D/jIa7SP39DviC/m4v19/SjYnoAFFRUaHHH39c//Vf/6WrrrpKkpSamqqKigqPdk6nUykpKT4fJz4+Xj169PCr1mApKSlRXl6e0tPTlZSUZHc5PklOTlavXr2C0t60rbc1hBOT2utqE8jrHqjr19h+wr2/B7sfBeqeBaI/hFq41xcM4d7fgUCiv9svNzfXuG1EB4hnn31WkydPVmZmpiSpsLBQXbt2rTVcqaioSF26dPH5OK6J2+EsKSkp7GusT0xMjFe1e9PetK23NYQTk9rrahPI6x6o62e6n3Dt78HuR4G6Z4HoD6EW7vUFU7j2dyAY6O/2MR2+JEXgU5hcVq1apRtuuMEdHoqLi7Vx40YNHDhQBQUFHm3z8/PVv39/O8oEAAAAokpEfgKxZcsWbd26VYMGDdLbb78ty7J08OBB3XTTTRo8eLDuv/9+lZaWKjExUSUlJdq+fbtee+01u8sGAAAAIl7EBYiioiI9+OCDOnv2rNasWeOxbuLEiUpLS9OTTz6puXPnqlu3bsrNzdWsWbPc74oAAJgZNWqUVq9ebXcZAIAwE3EBIjU1Vdu2bWuwTWZmpntoEwDAN7w8EABQl4idAwEAAAAg9AgQAAAAAIwRIABIujTePZjtAQBAdCBAAJDk/Xh3xscDANA0ESAAAAAAGCNAAAAAADBGgAAQUE1hbkRTOEcAAOpDgAAQUE1hbkRTOEcAAOpDgAAAAABgjAABAAAAwBgBAgAAAIAxAkQTx2TQ8MR9QaSgr9aPawMgWhEgmjgmg4Yn7gsiBX21flwbANGKAAEAAADAGAECAAAAgDECBGCI8cwIBfqZvQJx/bmHAKIdAQIwxHhmhAL9zF6BuP7cQwDRjgABAAAAwBgBAgAAAIAxAgRgI8ZKI5Tob8ERyOtac18Oh0Pjx48P2P4BIBAIEICNGCuNUKK/BUcgr2td++K+AQg3BAgAAAAAxggQAAAAAIwRIAAAAAAYI0AAgE2Y1Bw+uBcAYI4AAQA2YXJs+OBeAIA5AgQAAAAAYwQIAAAAAMYIEEAUYjx3ePHmfsTFxQWvkChhd/8OxvHtPicA8AYBAohCjOcOL9yPwLL7egbj+HafEwB4gwABAAAAwBgBAgAAAIAxAgSAOgV7nLedY76jYbx5XecQDuc1ZsyYRtv4W2dj80RMajARiDrtuCfh0A8ARDcCBIA6BXuct51jvqNhvHld5xAO52VSQ7DrDNT+A7EfO+5JOPQDANGNAAEAAADAGAECAAAAgDECBAAAAABjBAgEBJP2AN+Z/vz4+pK5YP18BmqyMsIfv+NDh2uNSECAQEAwaQ/wXaRMKg7VfhF+uNehw7VGJCBAAAAAADBGgAAAAABgjAABIKJ5M16YscW+82W+g6/XO9TbhYK3tdXX3rW8of35Olcm3IXz/QWaGgIEgIjmzXhhxhb7zpdr5+v1DvV2oeBtbfW1dy0P53MNlqZ4zkC4IkAAAAAAMEaAAAAAAGCMAAEg7FUf++xwODR+/HivtolEkV5/Qxobox8XF+fTnIu6tgvkdQzEvsLh3RnVzyOa+1m04V4hnBAgAIS9mmOfTcZCR/p46Uiv31+BmgcRyOsYiH2Fw32tXkM41AMz3CuEEwIEAAAAAGMECAAAAADGCBAA4Ac7xiWHwzj6+mqIi4ur85qE4jrZ8f6D6tfBdY7V52KEqiaTPjFmzBiPdmPGjKmzPm/7F2Pzwwf3AqESnW+bAYAQsWNccjiMhW6ohrrWhUPNwVDffIJQn68v84Iae9dEII+N0OBeIFT4BAIAAACAMQIEAAAAAGMECAAAAADGCBAA4KNwmMzsj+oTLuub/BxMrus3fvz4Ro/tmuxb/SWCNScE+6LmcUN1HRqrNy4uTuPHj5fD4VB8fHytbbydnF1z24Zexlj9Gphci+rta07Urq7mvny9zt5sF4p7GaqfG2/uOZOpEWxMogYAH0X6hMVgvnTNm+M7nU7jYzc0admXJx7ZNeE7kJOefTleY/uqfm+82b83k+sD9bLAQLX1VTj+HgjHmhBd+AQCAAAAgDECBAAAAABjBAgAEcM1JjzS1DeuvvoLxwIpmOOf7Zoz0NC18vYajhkzpsGaExMT69xnY3MHTNR1Ht7ss75hWjVfXldff6trWUPXouY519e2vmvmUt9L62pqrBbTbRqai+ELX/p4Q9uEao6CXXMhfD2uw+Hw+2cMocEcCAARI5LH9QbqpV3+HCtY+w7VfalvnoMvLz6rqKhotI03y709fqD3WXM/gZwnYLpfX9f50q6xbQLdJwNVlz/784VdvzMDOV8H4YlPIAAAAAAYI0AAAAAAMEaAABD1IvGZ6MGYW1Bz/L03jz1taPx/zVpdxzEZl2/C5J0JJtuNGTPGfQ7ezGcwrduX61v9+piM2a/ZLlh9u77rU/Oc6mvTWF2m51B9LoPJNarZD03b+sK1fWJiose+6ttv9XdlVD+X6u9Bca1r6D0cvsztcG1j+o4OoDHMgQAQ9SJ1TG24zY/w5Rn/vo7Lr7mvQLzjwel0yrIsr+sJ5nsHqp+bv+/CCDRf3lVhuq2v8yq82W9j/L12Db2TpKH2dX1taF0g6m7sPR2R+jsS9uETCAAAAADGCBAAAAAAjBEgAAAIU5H43hMA0Y8AAQARzJfJjw29jKz6hE5v9tfYBOLq39c1ydol0C/W82XuRPVt65skPH78eI/JvfUdq6Fzbay+pKQkTZ8+vcEaG3rpnbfHa+ja13UtvL1XdU22r7lPf+6X6X5cNVRv09jLBevavi4N3euG1OxDJpPVXX2w5s+d6bFNtqs5mbuhyd3e/t4wmXiO8EWAAIAIFowXZvkz2dnbtsF6qVqgNDTp1NtJz4F+GVljbXyZ0O3N+kCcjx33u76JyYGYVO/Pvfb2mK6afb2mJtsFc3J3qB4GgOAgQAAAAAAwRoAAAAAAYIz3QABAAAR67H4omb5Qraa4uDhVVFR41d7kWK7x36NGjdLq1at9qsvbcwrkePy6rovJtaqvTX21BWO+iMlL+7w9bvXa69retWzlypUec0YaqqG+muqai9NYva7jVb/2rnpq1m+67/r6QF3Hakhdc1zqOkfXz4pr3w3Ngan5c1X95YrVt3Gdf83t62tb/b5Ur7vmfrzpP9Vr9fX3QV37CiTXfoO1/3BFgACAAIjkMbyhrD1UL/oK1nwDu9Q35jxQk4692ae/cysaGz9f84V/3h6j5jmYzO2oec6mL4nz9h74cs/qq6+xczSdG+PrXARfX+rn6zEC+eK/QArWvJdwxxAmAAAAAMai9hOIBQsWKCYmRs2aNdORI0c0ZcoUpaSk2F0WAAAAENGi8hOIxYsXq7y8XJMmTdKECRP0ox/9SI888ojdZQFA0ARjDkYgxtmb7qO+cf71vYfB9Hn7rjb1zYnwZb5Eze193S7Y82YaO/dg1FJ9jHz1d2XUbNPQPIbG6qv+ffXx+9W/1ldTzeU152e4am7onRWu8f31jfmvq8665nV4c87eauhY1etv6P7U1ab6dajrGDXf51DXsprbVr8HDodDSUlJ7vdbVN+++n+7aqq+7+rXuub7KgL5s1bfuzCq19fQey3qamOyXbiJugBRVVWlhQsXaujQoe5l1157rbKzs/XNN9/YWBkABE+wx/eGYh++jOM2fU+CybP7Qy0UxzV9v0EwavH2nSLe3KNAjTuva78m7ygx7auN1RnKfllX/b6808Lbmn2ZK1HffWjofRTVl5m8t8JXJu/C8Pa6RuI8iqgLEP/617908uRJderUyb3M4XCoQ4cO2rp1q42VAQAAAJEv6gLEoUOHJElJSUkey5OSktzrAAAAAPjGYXn7rLQw9+677+qRRx7R/v37PZbffffd6tatm55++mmv9rdjxw5ZlqX4+PhAlhkwlmW5n/nscDi83v7EiRNq376933X4sx9vt/WmvWlbk3aBulYm+/T3HGsu86WNN9/X998N1WpSU83lx48fV5s2bVRYWKi0tDSjbar/d82vvqq+H0m1/ruu9nWtcy2vvq6ufTd23IaW1dWmZruTJ0+6H5vpcDhkWVatfZ48eVKSPNa52rq+Vt++oa+u/fjaxt/tfam1vnOufr1c19HfWqtfe9dyk1obqqOxery9Hi7Buj8mbfw515p92tdjNHZdqv+s+XvvGmtT17rq5+rLta+uvmtW8xh11VHz91hdNbdr104VFRU6c+aMamrsd5xrvya/H/3V2O/Xxv4f01jtgf47wxtOp1MOh0M/+MEPGm0btQFi3759Hn9Q33XXXerWrZtmzZrl1f6++OKLsA4QAAAAgL9cAeKaa65ptG3UPcY1NTVV0qWLkJCQ4F7udDp9eoyryUUEAAAAmoqomwNx+eWXS5KKioo8lhcVFalr1652lAQAAABEjagLEN27d1e7du1UUFDgXlZRUaGjR48qMzPTxsoAAACAyBd1ASImJkaTJk3SunXr3Mv++c9/asCAAbryyittrAwAAACIfFE3iVq69DK5BQsWKCYmRomJicrLy9PUqVPd8yMAAAAA+CYqAwQAAACA4Ii6IUwAAAAAgocAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwFmd3AQi+jRs3qri4WLfeeqvdpQBBs379ehUWFuqzzz7T7bffrsGDB9tdEhA0W7du1YULF7RlyxbddNNNuu666+wuCQi6v/3tb2rfvr1++MMf2l1Kk8cnEFGssrJSS5Ys0Ztvvimn02l3OUDQ5OXlKTExUXfeeaceffRRTZ06VUVFRXaXBQRFUVGRHn74YQ0dOlRjx47V9OnT7S4JCLrS0lK9/vrrqqystLsUiAAR1WJjY3XPPffo6quvtrsUIKgOHjyol156SZKUlpam+Ph4ffPNNzZXBQRHamqqli5dKkmKiYlRVVWVzRUBwffFF1+oe/fudpeB/8MQpgiRnZ2tefPmuf+nUd2CBQsUExOjZs2a6ciRI5oyZYpSUlJsqBLwny99/frrr9cVV1whSSopKdH58+fVsWPHUJcOeM3X3+3p6emSpE8//VT33ntvKEsGfOZrfz9+/Ljatm0b6nLRAAJEmPvoo4/02Wef6cSJEzp8+HCt9YsXL1Z5ebmmTJkiSdq8ebMeeeQRzZ8/P9SlAn7xp683a9bM/S9Ta9eu1e23387/bBDW/P3dblmW/v73v2vnzp16/PHHQ1o74C1/+3tOTo5uvPHGUJaMRjCEKczdeOONevjhhzVo0KBa66qqqrRw4UINHTrUvezaa69VdnY2wzcQcQLR1wsLC/XJJ5/oscceC0nNgK/87e8Oh0MjRoxQVlaWxo4dqwsXLoSsdsBb/vT33bt368orrwxluTDAJxAR7F//+pdOnjypTp06uZc5HA516NBBW7duVbdu3WysDggck75eVVWlxYsXa/bs2ZKks2fPqlWrVjZVDPiusf5eWlqqb775RiNGjFCXLl1UVlamLVu2aNiwYTZWDfimsf7erFkzlZWVSZK++eYb/eMf/1CbNm3Us2dPu0qGCBAR7dChQ5KkpKQkj+VJSUnudUA0MOnr77zzjsaNG6f4+HhlZ2era9euBAhEpMb6e1xcnAoLCyVdejLNxYsXmVyKiNVYf3/00Ufdy3bu3KkhQ4YQHsIAASKCuT6yjovzvI1xcXHuR1iuXbtW27dvV6tWrdSqVSsNGTIk5HUC/mqsr+/atUuPP/64+2k0CQkJ2rlzZ8jrBAKhsf4+atQorVy5Um+99ZYOHjyoJ554ggCBiGXyt4wkffzxx9qzZ48SExPVqVMnde3aNaR1whMBIgpYliWHw+HxvWVZkqThw4dr5MiRsixLFRUVdpUIBER9fb1v377KycmxsTIg8Orr7/Hx8brzzjttrAwIvIb+lpGkgQMHas2aNbWWwx5Moo5gqampklTrJXFOp9P96LP4+HhJl8YTuv4biDQmfR2IFvR3NCWm/T0+Pl4Oh0MxMTGKjY0NaY2ojQARwS6//HJJqvXG3aKiIj7aQ1Shr6Mpob+jKaG/RyYCRATr3r272rVrp4KCAveyiooKHT16VJmZmTZWBgQWfR1NCf0dTQn9PTIRICJEVVWVe4KoS0xMjCZNmqR169a5l/3zn//UgAEDeGYyIhZ9HU0J/R1NCf09ejgsZqKEtR07dmjTpk3atGmT9u/fr/HjxystLU2/+tWvJF36YXS9/j0xMVF5eXmaOnWqe0whECno62hK6O9oSujv0YcAEeYqKirck4YcDocsy5LT6VRCQoLdpQEBRV9HU0J/R1NCf48+BAgAAAAAxpgDAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAARzel0qry8XNUfKmhZlsrLy1VRUWFjZQAQneLsLgAAwt2ePXs0d+5cff311zp69KhSUlL0gx/8QLGxsaqqqlJRUZFKSko0fPhwTZw4UYmJiT4dZ9GiRZo/f75efPFFZWZmBvgs/Ddt2jR9++232rVrlyTpmmuuUcuWLSVJZWVlOnXqlK666ir9+7//u7p37x6yurKzs7V8+XLt3LlTR44cUXx8vAYMGKBWrVpp0qRJ6t27d8hqcTqdeuWVV5Sfn6/i4mK1bNlSs2bNqrNtzX6VnJysfv36KTY2VpJUWVmpEydO6Nprr9Wvf/1r97UGANtZAAAjBw4csDIyMqzf/e53tdYVFBRYI0aMsG677TarsLDQp/0/+OCDVkZGhrVs2TK/6jx//rz1/PPP+7WPhvzsZz+zevXqZVVVVXksLysrs/7whz9Yffr0sf7xj38E7fj12bRpk5WRkWFlZWWF/Ngus2fPtn7+859blZWV1tSpU62+ffta58+fb3AbV7969NFHa607ffq0dcstt1g//vGPrWPHjgWrbADwCkOYAMBQUlKSJCkmpvavzs6dO+uhhx7Svn37NGfOHJ/2P2PGDM2fP19jx471q87Dhw/7tX1jkpOT5XA45HA4PJYnJCRo2rRpSk1N1cMPP6zz588HtY6amjVrJkmKi7Pvw/X33ntPw4YNU0xMjP77v/9bn376qZo3b97gNq5+5frkobo2bdpo6tSpys/P1x/+8Ieg1AwA3iJAAECAXHnllZKkDRs2+LR969atdfPNN/v9B/Df/vY3v7b3R3x8vK644gpduHBB27Zts60OuxQWFroDQ0xMjFJSUvze58CBAyVJ//znP/3eFwAEAgECAAIsPj7eluOWlJRo5cqV+vOf/2zL8V2s/5vMbNd1iDYXLlyQdOnTCG9VVVWpqqqqznWVlZXuNq7/BgATTKIGgAD58ssvJUmjR4/2WH7x4kW99NJLKi0tVYsWLXTq1Cn17t1b48ePdw8DWr16td5++20VFhbqgQce0E9+8hNJ0pQpU7R3716dOXNGa9eu1ZYtW3T8+HGdOnVK+fn5mjVrljp06CBJ+stf/qJt27YpLi5Oa9eu1e7duyVJQ4cO1ZgxYyRJBw4c0OrVq9W5c2dVVVXp4sWL+sc//qG33norINegpKREBw4cUIcOHdz/ci5d+iN10aJF2rVrl77zne/o22+/1aBBgzRu3DhJ0oIFC7Rq1SoVFhZq7ty5ys/PV15enrZv367XXnvNpz+eTXz99dd6/vnn3dfw1KlTmjZtmjp27Ohus2fPHn355Zc6d+6cjhw5oosXL2ratGnq3Lmzu81rr73m/sRl8eLFWr9+vb773e/qt7/9rd81rl27VpI0fvz4WutycnK0ePFitW/f3l3/hAkT3BPHP/74Y7366qv67LPPZFmWbrnlFs2bN09nzpzR8OHDde7cOXXp0kW/+MUvNGHCBL9rBdBE2D0JAwAiRX5+vpWRkWFNnz691rqdO3daN954o/Uf//EfVllZmXt5UVGRNXr0aOv999/3aD9jxgzr//2//+f+vqSkxNq2bZuVkZFhrVq1yqPtsmXLrIyMDOuVV16xcnJy3Mvvvfdea9KkSbVquemmm+qcRF1aWmrdcccdVklJiXvZuXPnrB/+8IcGZ///u/vuu63evXvXWn7mzBlr6tSp1g033OBRp2VZ1vTp060pU6ZYFRUVlmVdmnA9cuRI67333nO3yc3NtTIyMqzZs2dbO3futNavX29lZGRYGzduNKpr69at9U5Grsv+/futgQMHWjt37nQve/fdd60RI0a46zx9+rR1zTXXWG+++aa7zauvvmpdd9111oULF2rts67715iG+tUnn3xiXX311dacOXOsyspKj3Wffvqp9ZOf/MQ6efKke9nJkyet4cOHW59++qlH20WLFlkZGRnWjh073MtmzZplzZs3r9ZkeABoDJ9AAICXduzYodmzZ0u6NFzn7Nmzio+P11NPPaXrrrvOo+2zzz6ruLg4jRgxwmN5VlaWBg0apMGDB2v48OFKTEys99GtrkeiXrx4Ub169XIvv+qqq/TXv/7VuO6vv/5ap0+f9hha1KJFC/Xr1894Hy6VlZWaNWuW+xOU0tJSnTt3Tv369dPMmTOVmprqbrtp0yatXLlS7733nnuicEJCgkaNGqXnn39eI0eO9DjPvLw89e3bV71799b8+fN1/fXXe12fienTpyszM1N9+/Z1Lxs5cqRmzpyp999/X7fddpuKiooUGxurkpISd5sJEyZo7ty5eu+993TnnXcGrJ79+/drwYIFki71q6NHj2r79u3KysrSmDFjPCbvFxcXa9q0aZoyZYratm3rXt62bVtNmDBB06ZN0wcffKDk5GRJ0r333qstW7bo6aef1ltvvaWzZ8/KsqyAfEICoOkhQACAl3r27Knf//73jbYrLS3VypUr63yqUosWLXTFFVfojTfe0PDhwyWp1lONanJN0nZJSkry+MO2MZ07d9bJkyd13333aezYserfv7/atm2rrKws431U99BDDxlNEl61apUSExOVkZHhsfyqq67St99+q7Nnz6pVq1bu5X369JF0aQ7FzTff7FNtjTlw4IC++uor97V3iYmJUa9evbRr1y7ddtttuvzyy7V9+3aPe5OQkKC2bdvq1KlTAa3pyiuv1OTJkz2WVVZWas6cORo5cqTmzJmjH/7wh5Kkv//97zp16pRH+HH5/ve/r1OnTmndunX62c9+5l7+9NNPa+TIkXrxxRdVWFjo830HAAIEAARJbm6uysrKPP44rq5Vq1bas2eP8f4uu+yyWsvqmyBbl5YtW+p//ud/9Pjjj+uhhx6Sw+HQ9773PT3yyCPG+/BFbm6u4uPj3f+67nLx4kVde+21Ki0t9Vjerl27oNWybt06DRkyRLm5uZKkr776qlZdaWlpHv+q73A4lJOTo48++kilpaVq2bKlV8HNH7GxsXr44Ye1ceNGTZo0SatXr1Z6err27t0rSXX2LdeyvXv3egSIdu3aacaMGZo6dar++Mc/Nvp4WQCoDwECAILE9WQb6/+eSlSTZVkhefrNgQMH3P/6P2zYMA0aNEhbt25Vdna2PvjgA9133316++231bNnz6Ac3/U405r/ul6fhISEoNQhXZrofsstt7iHA1133XW644476m1fVlam6dOna/fu3XrmmWd09dVXS5LeeOONoNVYU0JCgn70ox9p+fLlWrVqlbKyslRRUSGp7r7lWuZqU1337t3VoUMHvf766xo5cmRQrzWA6MVjXAEgSLp166bY2Nh6h7qcPn1aPXr0CHod69evlyRlZ2drzZo1at68uYYOHarHHntM69atU79+/bRmzZqgHf/KK6/U6dOn5XQ6a62rqKioN2AFw5EjR9w1SdKxY8fqbOeqde7cuVq/fr0WLlzoDg815efnB77QGlzvBrl48aIkucPe6dOna7UtLCz0aONSXl6uN954Q++8847OnTun5557LpglA4hiBAgACJIWLVroJz/5iXbs2FFr3enTp5WXl1fnozn9lZKS4jHExjV+v7KyUh988IFH24SEhKDNM3AZP368nE5nnS+WW7JkifsP3mBzPQ5XuhTuBgwYoM2bN9dqV1JSoldeeUWStHnzZvXq1UuXX365e31FRYVOnDjh/v7dd98Nat3l5eXuOl2T9EeMGKGUlBRlZ2fXav/555+refPmtSbuv/jii5o8ebLatGmjWbNmadGiRfrss8+CWjuA6ESAAABDZWVlklRrzH5DHn/8cVVUVGjZsmXuZZWVlXrqqac0YsQI9/sZpEt/KEqq9S/1ruWurzWX1xwG1b9/f3311VceNbt88MEHtf5o/Pzzz2tNJm6Ia581912fzMxM/epXv9Izzzyjs2fPupcfPHhQpaWl7rkdrvPx5vpW5wpNdQ0L27dvn6ZOneoxv+KPf/yjjhw5ouXLl7uXVVVV6U9/+pN7WFPv3r2Vn5/vfpmbJC1fvlxXXXWVzp8/r7Nnz6p169aSfOsf1euuaz5LWVmZnnjiCeXl5WnMmDG66aabJF16qdwzzzyjxYsXq6CgwN3+0KFDWrRokZ555hl3XZK0dOlSFRQUuIPQ9ddfrx//+Md69NFHde7cOa/qBQCHFcrPjgEgAu3fv18LFy7Uvn37dPDgQcXHx+u6665Tt27djCYgnz9/XvPnz9fZs2fVvHlzHT9+XJmZmbr77rvdY/GXLVumd999Vzt37lRaWppuuukmzZw5U08++aQ2b96sQ4cO6Tvf+Y5uuukmTZkyRY899pg+++wznT17Vv369dO//du/6a677pIknTt3TjNnzpRlWWrZsqV+/etfKy0tTTt27NDHH3+sdu3a6fTp03I4HDp37pwGDRqkIUOGNHoes2bNUn5+vjZt2iTLsvTd735X3bt31wMPPKD09PRGt1+7dq3WrFmjLl26KCYmRh06dNDEiRMVGxvrcf6tW7dWv379NHbsWN1www2N7vfjjz/W0qVLtXfvXp08eVLJyckaMGCAYmJiVFlZqSNHjujAgQOSpIkTJ+p3v/ude9sTJ05o3rx57iBTUVGhsWPHuh+Xe/HiRc2dO1d5eXnq06ePqqqqdO211yo5OVkPPfSQBg8erKysLK1du1bvvPOOvvrqK1122WX6wQ9+oF/84hfupybVZe/evZo3b54OHjyow4cPKzU1VZmZme4+UVxcrLy8PLVv31533nlnrRcUStLOnTv1+uuvq1WrVqqqqtLZs2c1adIk99OZVq9erZdfflnffPONUlJStHbtWnXq1En79u3TPffco3PnzqlDhw76/ve/r+eff77Raw0AEgECABpVWVmpyspKxcXFKSYmRpZlqaqqShUVFWrWrJnd5YVMaWmp4uPjFRMTI4fD4b4GsbGx7vc7BILr+koK6H5DxVW/w+HweHeDHVyfxkTidQQQvggQAAAAAIwxBwIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjBEgAAAAABgjQAAAAAAwRoAAAAAAYIwAAQAAAMAYAQIAAACAMQIEAAAAAGMECAAAAADGCBAAAAAAjP1/TB1tc6KgTa4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# # Plot the histogram using seaborn\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "\n",
    "# Convert the dictionary to a pandas DataFrame\n",
    "data = []\n",
    "for points, occurrences in stats.items():\n",
    "    data.extend([points] * occurrences)\n",
    "\n",
    "df = pd.DataFrame(data, columns=[\"Points\"])\n",
    "\n",
    "hist_plot = sns.histplot(\n",
    "    df, x=\"Points\", bins=len(stats), color=\"black\", log_scale=(True, False)\n",
    ")\n",
    "\n",
    "# Load CMU font\n",
    "import matplotlib.font_manager as fm\n",
    "\n",
    "font_path = \"/Users/sri/Downloads/cmu/cmunrm.ttf\"\n",
    "font_properties = fm.FontProperties(fname=font_path)\n",
    "\n",
    "hist_plot.set_xlabel(\"Points Per Leaf Box\", fontsize=15, fontproperties=font_properties)\n",
    "hist_plot.set_ylabel(\"Count\", fontsize=15, fontproperties=font_properties)\n",
    "\n",
    "hist_plot.set_ybound(lower=None, upper=117)\n",
    "# Set font properties for tick labels\n",
    "for label in hist_plot.get_xticklabels():\n",
    "    label.set_fontproperties(font_properties)\n",
    "for label in hist_plot.get_yticklabels():\n",
    "    label.set_fontproperties(font_properties)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"stats.svg\", format=\"svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d00829b-4de1-4b75-8cce-571390d86b7e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2678b7d-7d6c-4c12-9a53-7b9693631619",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}