07_diffusion_distillation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "af32d529",
   "metadata": {},
   "source": [
    "# Diffusion distillation (WiP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d2c8c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "import jax\n",
    "import jax.numpy as np\n",
    "import flax.linen as nn\n",
    "import optax\n",
    "import diffrax as dfx\n",
    "\n",
    "from sklearn import datasets, preprocessing\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import trange"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6fcdfbd3",
   "metadata": {},
   "source": [
    "## The dataset"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c6de9e86",
   "metadata": {},
   "source": [
    "We'll use two moons to keep things simple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "cdab2cc6-b2dc-4834-b343-0386a25fd9a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.0, 2.0)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGiCAYAAADulWxzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmGUlEQVR4nO3de3hU9bk3/HsOmck5Q84JBAgHAQUBUTDYKlQU1G2lutlqu+vhsdi6tW8t7rbi260Pdvdlt9Vit9sW3d1K+1Rba6tYbatFBKyKqEDkjJwTQs5hMjnOZGbW+4dPU9f63sBKyOSw8v1cV66LuVkzs2bWmskv675/v9tlGIYhRERERA7kHugdICIiIkoUDnSIiIjIsTjQISIiIsfiQIeIiIgciwMdIiIiciwOdIiIiMixONAhIiIix+JAh4iIiByLAx0iIiJyLA50iIiIyLESOtBZuXKlXHTRRZKRkSH5+fmyePFi2b9//xnv98ILL8jkyZMlOTlZpk2bJn/6058SuZtERETkUAkd6GzatEnuvvtuee+992TdunXS1dUlV155pbS1tZ3yPu+++67cfPPNcscdd8j27dtl8eLFsnjxYtm1a1cid5WIiIgcyNWfTT3r6+slPz9fNm3aJJdeeqm6zY033ihtbW3y6quvdscuvvhimTFjhqxevbq/dpWIiIgcwNufT9bc3CwiItnZ2afcZvPmzbJs2TJTbOHChbJ27Vp1+3A4LOFwuPt2PB6XpqYmycnJEZfLdfY7TURERAlnGIa0tLRIcXGxuN19l3Dqt4FOPB6Xe++9Vy655BKZOnXqKberqamRgoICU6ygoEBqamrU7VeuXCkrVqzo030lIiKigVFZWSmjRo3qs8frt4HO3XffLbt27ZK33367Tx93+fLlpitAzc3NMnr0aKmsrJTMzMw+fS4iIiJKjFAoJCUlJZKRkdGnj9svA5177rlHXn31VXnrrbfOOEorLCyU2tpaU6y2tlYKCwvV7f1+v/j9fohnZmZyoENERDTE9HXZSUJnXRmGIffcc4+89NJL8uabb0ppaekZ71NWVibr1683xdatWydlZWWJ2k0iIiJyqIRe0bn77rvlueeek5dfflkyMjK662yysrIkJSVFRERuueUWGTlypKxcuVJERL7xjW/IZZddJo8++qhcc8018pvf/EY+/PBDeeqppxK5q0RERORACR3o/OxnPxMRkXnz5pnizzzzjNx2220iIlJRUWGqrp47d64899xz8t3vflceeOABmThxoqxdu/a0BcxncoV7Sa/vS9TXXB5Pnz6eEYv16eMRESXKuvgL/f6c/bqOTn8IhUKSlZUlzc3N3TU6HOjQYMKBDhENV6cb6Gi/v/sCe10RERGRY3GgQ0RERI7VrysjEzmdlpayppb6OtVk5zmJiIYrXtEhIiIix+JAh4iIiByLqStyBLvpG7c/2XQ7Hu7s9WP1dvaUJxCAWLylxdZz2r2vdd+YyiKi4YpXdIiIiMixONAhIiIix+JAh4iIiByLNTrkCHZrULSaHDustT0iIka0C2IubxLEPPl5ptuxunpb99NiWj2O+//2jTPtW1e0x9uInF3NEhHRYMQrOkRERORYHOgQERGRYzF1RYNGX0/r1lI/8PhK+klL87iysMFctOqErfuK2/z3hKewAPejrR1jHZhG0u4rkQiEXBnmj3a8KYjbpGA6zmXz/Yh3dOB+2MCUFxH1N17RISIiIsfiQIeIiIgciwMdIiIicizW6NCgcTb1OGptjA2uJPwIuHJG4IbJfgh5iwohFh2DNTSednMNTSzVB9u4w8pU9ZiB+1GNU9Mj546GmO9gjfnxC/NgG4kqLSayMiBmNNtrT2FnCj6nqhNRf+MVHSIiInIsDnSIiIjIsTjQISIiIsdijQ4NanbrdrR1blRey329+BEwks+8/o6IiORhLU/cj/vWlZlmut2Zg88ZT0qFWHITvvakEViL1JWhvIbJRebH97jwsUJYF+StD0FMq1nypuH+amI1tWfcxm4dlob1PUR0JryiQ0RERI7FgQ4RERE5FlNX1C+09IR1SrjWVkBtx1CgTJW2qWN8jul2chWmalomYarGUFI/2vTv6kvwb4fUavN9vdjtQRVJx49nezHGAh/jfljvG1eycXnlcXz8ibkQ8wWxxYSnBY+VK9iK22Wb38t4CLfRWlFoXdqJiHqDV3SIiIjIsTjQISIiIsfiQIeIiIgcizU61C/sTAP2jCrGYASnQBtpWNOh6VKmYreONBeruGPY8qBxGtYTpVdgHUzwHPw7IZ6qvU7z43Vd0QxbtNWk491Soxjaj60oQqVYPxQZYa6/cUVwG3eXvffRHcMCn/zNWLfTPBdbUWQcNNdAudsDsI1R24DPmYHHJRYMnmYv/85aD8Yp6ETDG6/oEBERkWNxoENERESOxdQV9TltKrl1mrEmXoOdud1jR0FM6/7dUoppKsON6ZqWseZYOICpoM48THV0KjPa0yrxdY66sAZiR8XS5bxKSVNlYJrqovHHILbVUwKxWDumljyp5pRfchqmmoIFSuoqin/7eOvw8dvz8XgGDuJ09Wim+Tm0Lxy30jFd6xbv7uiEmNZ9PtaKU9iJaPjiFR0iIiJyLA50iIiIyLE40CEiIiLHYo0O2Wa3y7TWtkHrEi5Rc11K12fOw+eM4rTuzhysGYknYT1Oez7GwqPMtSuRgL3X5MvHlgedOfh3wp1j/gqxtwLnnPHxK9uw5iU9CWtSrpy0F2L/lPM+xP51zxLT7aZGrAv6l1mbILa6/LMQ68rB9zEawNhJa2d4EYlYanTSavDYpWg1XfVBiLkz8TVE63FquifdvB1rdoiGt4Re0Xnrrbfk2muvleLiYnG5XLJ27drTbr9x40ZxuVzwU1ODBZ5EREREZ5LQgU5bW5tMnz5dnnjiiR7db//+/VJdXd39k5+fn6A9JCIiIidLaOrqqquukquuuqrH98vPz5dAIND3O0S2aWkqbYVZa5pARO9GLW5lTJ1qTnEFJ+CU4nAA76bF4ik4tTllDHYmT7Pc/n8/+2fY5pnjcyEW8GEa6eaC9yBW2ZUDMb/bnKK7Necd2GZl+BqI1XXgtOuStJMY82Kn76baTNPtsaPrYBvN6IImiHXF8Vyoqs6GWNqsIMRCPvP7kYIrCEhnIaY6vVm4hID/OK4o7VGmnFvZPZftbkdEQ8ugLEaeMWOGFBUVyRVXXCHvvIO/FD4tHA5LKBQy/RARERGJDLKBTlFRkaxevVp+//vfy+9//3spKSmRefPmybZt2055n5UrV0pWVlb3T0kJLqhGREREw9OgmnU1adIkmTRpUvftuXPnyqFDh2TVqlXyf/7P/1Hvs3z5clm2bFn37VAoxMEOERERicggG+hoZs+eLW+//fYp/9/v94vfj7UddHa07tHW6eAiIq4cpbWD0nE8WoI9FJrPST3jfrSNxRqJtEKcLpzmxxYHdfWZEJteetx0e3f7SNhmfAZOWa7uyILYzyrnQayyCd+PeWMOmG7vixTBNu/vGQexqedUQuxQSy7ErthxL8TSsttNt28f/S5s89B7n4fYonN3Q6ws8xDEHu1cALGUJDzuQUvtVPDLeOyMD/C9DRzAi82d2fjaUwrxGPv2VJgfX5leznocouFjUKWuNOXl5VJUhL8YiIiIiM4koVd0Wltb5eDBg923jxw5IuXl5ZKdnS2jR4+W5cuXS1VVlfzyl78UEZHHHntMSktL5bzzzpPOzk75+c9/Lm+++ab85S9/SeRuEhERkUMldKDz4Ycfyvz587tv/62W5tZbb5U1a9ZIdXW1VFT8/TJzJBKR++67T6qqqiQ1NVXOP/98eeONN0yPQURERGSXyzAMXGN/CAuFQpKVlSXNzc2SmflJ/v4K95Iz3Ius3H5cC0dbH0er0Yln4roojTOwDiOaYm4j0IllPBKZiK0XAoF2iE3Pq4KYdf0aEZFLsz423dZqdE6EcV9nZGC9zP8cwPV2WlvwPbLub9PRAGyTNhLrSNorsE7K8ODHVbtva715xaDS0lrY5nMFH0PszVpsVzE5gPd9v3YMxF6e/j8Qu+6jO0y32ztxfZzONoy5GzBW9A6ulZTUinU1KduPQszKUNbfiXfguca6HaK+tS7+win/T/v93RcGfY0OERERUW9xoENERESONeinl1Pi2e1KrqWpOsZjy4OW0Zh2MJSn6LTcNZyNqQkjjKdoNIbj86X52In7oSPXQazYb24j8Ofj58I2wSBOe38/DVM1MWU/fPsxbReOmGN+ZVZ90h5Ml7nGY5rK8GGs4ximuPzF5nTZkSMFsM0vP8S0XbQIp+k3teMOl47AVhHXbFsKsfunvG66/f3d2BLG48P00ORZRyF2rLoUYql1eAwi8yaYbmfux9YRchhTkS4vdlbX2GkfwZQX0eDBKzpERETkWBzoEBERkWNxoENERESOxRqdYchaT+BOwboSV0Y6xIxkezUMSW1Ya3NyMo6p/TNOmgNdeDo+eP4fIbbuJNbVPN80B2L/ULgTYr+tnGW63VQZgG3Sjmo1S/geRfKxXib7OMY6c8zT6F04612V+xHG3F0uiLXnK3+vHDEfv3ghbuLFGdZi+HHnQpU4zfNgHJ9Ta8PxH3sXmm7/94xfwjZaLdWEjHqIHZiN6w9Ed+N5GvjYfAxcYWxNYYwphpg7iNP0o1UnIMb2EURDC6/oEBERkWNxoENERESOxdSVw9mZOu7KUlagTMcpxV25mCZoHYnprKap+HBlc/dArCg5ZLpd4m+EbY5HsvHxw2kQq+vAKdYv75oOsYxyc6f7NJwJr06FT2rBWN42e/cdccCc1nDFML3lC2LaJ5zjh5inQ0uR4Mc4qd38HN5O3LHkJkwxdtbgcW+agc/Z1YWPF01SOo6HzefHtz7GVcqnZWN6aMPxific1bhvHlyIGtKk7fn5sE32Xny/U1pxZWTveJzSHj9ejU9qwVQW0eDBKzpERETkWBzoEBERkWNxoENERESOxRodh9OWtXcXKm3CLSKFWLdjeHFqc9N0rDeJZ2B9wtEQ1trsbTS3JbhqFE4DXnvkfIh1duJr8hxS6jeUdgk+S62NvxnrVLwdGEupwqnHhlL/5KnCadFGnrl1huHB99FdF8T92IVd2l1p+DqT23A7GWmuS0lqxRqmrnSltqcV37PCt3B/Q6VYr9VUglPwJcV8LlRHsNVFWf5RiD027XmIrcy4Gh9fUdlkfr9bqrCmK7UOX7s/G98jdxMWZxlRPE+tOAWdaPDgFR0iIiJyLA50iIiIyLGYunI4VxIeYiPdnGKIjsB0SGcOpoeOL8K0xqgxdRDriuFl+9QknM57YW6F6fYrx3BeupamStqD+5uKuyFJmG0Sw2N+DRkHQ7CNlq6QqM20QxRXFna1WZYgbjwJ2+A7K2J0KEsXK7R0llimSvu24BRu37gSiGGSR6RrBKakfC14XjXGlFRYi/lcaMWnlDcDOJX8zSqMfW7kAVvb+ZLMx6BTSWG25ytpxzCmrjKUpQBcNXhOulIs89yV8yDWqpyQRJRwvKJDREREjsWBDhERETkWBzpERETkWKzRGaK06aue7BEQM/JxWrdVaBzWYLTn45RiTyouka+5ongfxM5JqYHYb6svNN2elINTs3eUnwMxr7IbmUeUJf3310IsbplCrE3rlgg+ltGl1N5Y6zJEJNaE9TeukLk2Q5uerC0DEA9jjY5Hq7lqxjoj675p+yrW2iERMZT6IV84F2JJDfg3UsyH51/tRebzVOuYrnWQzx/TBLGrs8oh5nfhe9nYZa40qh6BU9p3BYog5u7CWqeYH+8bEGwL4ao2n7vRYBC34ZRzogHBKzpERETkWBzoEBERkWNxoENERESOxRodB1HXXVHqAjqKzcv3t47CehzXRc0Qu7l0B8SylIKZTQ24tonfjTUu+94fa7qdVqXsB3aikPztYYh52vDxJdmPj3fAvHaPun6NUo+jibfgejt2am3UWg0bbQVEROId9uqk3F7zR1vd1xCu66KtuyTKdtKCscyo0jqjznyuGR7826r2Iqwfaq7DNiX/6/BSiKWNxP0ozjKfu6PSgrBNLKJ8LgogJP5mPCfjPnyPPJb3m/U4RIMHr+gQERGRY3GgQ0RERI7F1NUgo04bz8fL+FqaSmsFEFGW72+Ybk6vtI3FVI3nMOaM8ifhNOZfHL4YYkXpuN3T2+dCzN+BaQGr7H2YDvEfVPo9KGLHse2BO8X8fmjL8mvHQKOlInqbntDuZ3c/NDFlerOd5xTMCorb5tR6O20yvCE8b3N34N9boVIfPud52KXd48bzY3xGg+m2ljZ97tKnIHZb5dchVn8BhCStBvfN3WH+vLiagnhHhd3zhakwot7jFR0iIiJyLA50iIiIyLE40CEiIiLHYo3OIKPWfSj1OGqNTg4uwe8O4+PFLOUV3gxseZCch9OdX62ZBrGmowGIRUtw/GyEscbAb1nlP6MS9zXzI6UeR6kF0dogaKzTs+3WPmjb9baWx25txWCuwdCmq7ujWAuTFDJvZxRhvZmvCQuDUtKV93ZdBoQaZ+FzhgvNX2sz0ytgm1VVCyEWCeD77QvifjSei8sWBJLM53xqq1JXp7TX0JYL6MvaLyLiFR0iIiJyMA50iIiIyLGYuhpk3H5t2i5enndnByAWHoudyuvPx8vshmV4O3VkNWwzMQM7iZf4GyH2tg9TaHvXYsfxgDJtOa3GPDU4czc+voQwRaKtXGx3mnhfTv8eLuyu2qwdF6PLfFxc+XiOek7iOeTPwinchgeXI/AqqaUN+yaZb7ecB9s8d81PIfbfSZfhY+2aArHMQ/i12ZljjiWV5MA2Se2YptL+0tTOZSLqvYRe0Xnrrbfk2muvleLiYnG5XLJ27doz3mfjxo1ywQUXiN/vlwkTJsiaNWsSuYtERETkYAkd6LS1tcn06dPliSeesLX9kSNH5JprrpH58+dLeXm53HvvvfKVr3xFXn/99UTuJhERETlUQlNXV111lVx11VW2t1+9erWUlpbKo48+KiIiU6ZMkbfffltWrVolCxfiLAkiIiKi0xlUNTqbN2+WBQsWmGILFy6Ue++995T3CYfDEg7/vQAkFLI3zXiwsNaRWDtdi4h4ktIhFj53JMTiSg1DGEsFZMqlh0y3i1KwU3lXHC/25XmxXuajjUo9Th32BE+pw/qN1J1VpttGAF+naHUfWpd2xXCuq0k0O13aRfD8NvYdgm082bgsgq8ea9W60rEtSc4uPOdDLea6tGkLP4Ztnqi5HGLTM47j41+AbSdeiF0IsYyPzV+lngjuf9YRpe5I6yDPdg9EfWpQzbqqqamRgoICU6ygoEBCoZB0KOtNiIisXLlSsrKyun9KSkr6Y1eJiIhoCBhUA53eWL58uTQ3N3f/VFZWDvQuERER0SAxqFJXhYWFUltba4rV1tZKZmampKRgF24REb/fL34/TqEejOyspOvNy4VYfFQ+xDqzMXVQc7Ey/bYEL42nJ5nneqd7cWXk322bBbENuRPxsZRxZeZhvPoWTcVTzWgzpwVcXuWSvTK13u50Z6YAes7ue9bb1Z3V5RO0/ajG5Q3S25S07njMzSY3mf9++8iP6VVr+lZE5HAHfvbGpTRA7KqZOyD2Xrm5zbkngunbWBHuqyeCnz1tmr62EjXPZSJ7BtUVnbKyMlm/fr0ptm7dOikrKxugPSIiIqKhLKEDndbWVikvL5fy8nIR+WT6eHl5uVRUfNJ7Zvny5XLLLbd0b/+1r31NDh8+LN/+9rdl37598tOf/lR++9vfyje/+c1E7iYRERE5VEIHOh9++KHMnDlTZs6cKSIiy5Ytk5kzZ8qDDz4oIiLV1dXdgx4RkdLSUvnjH/8o69atk+nTp8ujjz4qP//5zzm1nIiIiHoloTU68+bNE8PAXPXfaKsez5s3T7Zv357AvRo4trpi+3AKqrsJ8/O+5lSIeTvwvp8rxam1U9LMLR+ePTobtvE04qnR0poFsUx8SmkvUqbWvl8FsbhlmrjapV2px3FC9+/BKtHvmTYFXZQVIbTj7lbOD38m1u51FqSZbidjuY98tGMsxHa14999j9/wNMTeaRqPz2lpVh4JYK1TYSt+WDxKKxepxbogbTq/W6lbZPsIIjSoanSIiIiI+hIHOkRERORYHOgQERGRYw2qdXSGIy3PbtU6oxhiTZOxBiB9RiPEDrXg2iCVbeYl9xuasPVCLIC1Glm78XQp/lM1xFTt+srWn6bWb5Dj2WkdIaKvLyN7DkOo89yZptvpJ+KwTSwZH7/t3DDEXmi4CGLWdahERMbPP2K6vXfrWNimuRTX+8o/iI8lytpRriT87LEeh8geXtEhIiIix+JAh4iIiByLqat+pC1/70qxxJQ2CHGlc0QUZ5dLnh9TALUtGRCblGOeb+vy4BIAScr08oxKTGe1TsX2FOnv4vL6djqOs2XDmdlpIyIy9N83u/vvzsDze8TLO023I3MmwTZthZhGSj6IsRkXYY+TackYWx86z3TbexGmyw4fx2np4Qn4+fFHcGp9vCkIMaJEsPMdM9S+X3hFh4iIiByLAx0iIiJyLA50iIiIyLFYo9OPtGm07swi020jDet4ghMxZ+o9vxlil+ThVNuwgUvH720uMN2OtuDS9Fk1EJLkRqwd8O2pwA0VnArbN4Zabrwvaa9dq/2yTsX2tOF07SSlHUM42wWxV2umQez9lFKIVVnao4TC+DmOBCAkcQ8+p5GdCTGX9jqVNhnWVhFcsoFOZ7jURfKKDhERETkWBzpERETkWExd9SNterl0mldGbVOma3fk4/TvolRcafhEGLuL/3DkXyB2yYGvmferA8e7he9gx3RPfRBiWrf1WE0tbqewvh+8zG6mnS9q+tPmdk6kvh9ifj+8h6pgm/TMsRDrSsdz+dCJPIwJxmaNNU85b2rH9R/8M05CrHMffmZ9TTjN3WNdhkJEpAU/o1rXd3K2s0k/advZebyhlvLiFR0iIiJyLA50iIiIyLE40CEiIiLHYo1Ogmh1E+5RRRCzTifvyMaxZzw3ArGOLpw2vuGD8yB2ybrzIZY8zVwr4K3F5zS8yhjYq5wuNrqSn4oT60js5K7VWi0F6y3OzNZy9crU7OSaNohlpuBjpTRgvUz95/Hx9jea63b+efwHsM3hjlyIvVN4AcQ6R2Bbi8LnsO2E1v4irtTtkLPZrbOxLj1wKtCWSPAzNNS+u3lFh4iIiByLAx0iIiJyLKauEkRNOyidyYPnmaeXduTjSqnJaZi6Gh9ohNh/XfMcxL745tcglvr6CNPtvMO4r96TeGlfS1MZXbjq7HBhd4qldTvt3NAuK7tTUnC7LFw1N1p14rT76WS9ndLqqqiGWHobXo5vm4JTyWPteKxKisxLKnwQHAvb7KkvgFjHGOxyXvKacg4VYNpLTS/vY+rK6eD7RPkMeLJHQCzWhMsbePLx/Laz2rgr2rvvvoHCKzpERETkWBzoEBERkWNxoENERESOxRqdBNFqLgw/xuJJ5pqcKK4cL4tK90Jscy12UL55z10Q87TjWDa13lwX4G9Qpgo2BSGk5XgH87LfA8FO3Y43D+stJDuAMUt7EBERIysdYt441nnELcdvqE0HPRt2zklX1N5SCYbSXTz7fbzvXn+h6fbNUz+EbT7Yj59ZFz68dORrX8tYc5G69Qg+nuV7h59P5/Ok43eCRluOQGOtxxERkTTzLyaXzdrMwXL+8YoOERERORYHOkRERORYHOgQERGRY7FGpw+oNReZSj5UWaejI8ecpO9KN2CbUBSX5G4L+yDmimHC33cSx7JJLeZ1XDxV9bCNtj7OYMm3DhZ2l1T3jCo2B5RWGrEMXDPHo9TouBqbIWa0tEJsONXk9IrWzkSRWoHr0rTnBSDmOmH+jDZOSoNtnp7/PxD7X3+9HWKdI7DtRHY5HmP1OyZk3s7uWk/kfFprB0nF7x2tNrCrNN9029PQZOs57ba7STRe0SEiIiLH4kCHiIiIHIupqz6gLq2tXBqPjsEl4N2WDFHulAbYZncjdj2fkIPbHXHjNOPUt7Igllxrae+gXL40lOnldpfzHi6Xxu22cpCo+f0w0vH9bi3FVEdmRJnCqS2zrkwvd1tSj9q+OvU42VkiX+vy7VEu7bvbsf1KRMkYWdPG757AqeT7gvj5v27qRxD7c+VsiIUm4ec4rRLT196Q+XVF6/F7goY2O59bIz8bYlp6PJaGv6fiypIKyVUh021XLj6+0daOj6V8zgYCr+gQERGRY3GgQ0RERI7FgQ4RERE5Fmt0esh2nUoeLtnefA72dwhOM9dSzM7Eep8xaTiV782qibb2I7kR6zy6ss25Wt+BCluP5dSaDju0467V47i1Vg4WbRPw3NCcnIaPlbVfmUo+vhBi3g7z9HKjBWt0nMpWDYOyTbS6BmKejgDERv0B73v4FnMdXWdYqdVS2rt81DQSg+fgMQ4FcZn/1GqspdCWhbDilPOhw86xcmVlwjaxFFyiIJyLsfY8fPyUJqz5ixSaC9P8nUrNnzLl3O4SHInWL1d0nnjiCRk7dqwkJyfLnDlz5P333z/ltmvWrBGXy2X6SU4eHHPxiYiIaGhJ+EDn+eefl2XLlslDDz0k27Ztk+nTp8vChQulrq7ulPfJzMyU6urq7p9jx44lejeJiIjIgRKeuvrxj38sS5culdtv/2QF0NWrV8sf//hHefrpp+X+++9X7+NyuaSwEC/Ha8LhsITDf1/JMRQKnWbrnrF7idfu6o/eDlz12F9nPgRjLsLLfx82lEAsGsMxavyvmBLpSsf9zXhzv/l+HUrXZl7KNlGPe4qyqmgmphjaJpqnYnZk47HTumR34YxzSW7Ecy2pVZnmrnUgttDOW66obKa9j7EReIxzPzKfH1VjcOp3ezrGitLx+6p0TCPEtpSfD7FwDqYiknLM3wHuDjyePMZDm7VbudGM51Dr3NEQa8/D750YnkIS9+HvPV/I/P3k36es2q59VlqVFb0HQEKv6EQiEdm6dassWLDg70/odsuCBQtk8+bNp7xfa2urjBkzRkpKSuS6666T3bt3n3LblStXSlZWVvdPSQkOCoiIiGh4SuhAp6GhQWKxmBQUmBfKKigokJoaLPwTEZk0aZI8/fTT8vLLL8uvfvUricfjMnfuXDl+/Li6/fLly6W5ubn7p7Kyss9fBxEREQ1Ng27WVVlZmZSVlXXfnjt3rkyZMkWefPJJ+d73vgfb+/1+8fuV629EREQ07CV0oJObmysej0dqa2tN8draWts1OElJSTJz5kw5ePBgInbxtLS6DGt+VETENRJfS9toXCdey5EmTzNPJ9emjack2ZsanL0Pp5amVGGO1JpLdWfgvsaCQVvPOVxo9VpabZMnhO93Uot5+mfrTPzYdZTgsdO60WfvxTovTxu2KTCK8ky31VoNZf/JzFDeN3cH1ickN5qn0eZuwj++rvjGPohtrMXP+9TMExDrzMHjrolnKnVjFpxePnRo07Otn1vPqGLYJqUOvxOazsWaPI9SrtWpnELJltLR8GRsS+Tb0nf1sX0toakrn88ns2bNkvXr13fH4vG4rF+/3nTV5nRisZjs3LlTiorwjSUiIiI6nYSnrpYtWya33nqrXHjhhTJ79mx57LHHpK2trXsW1i233CIjR46UlStXiojIww8/LBdffLFMmDBBgsGg/OhHP5Jjx47JV77ylUTvKhERETlMwgc6N954o9TX18uDDz4oNTU1MmPGDHnttde6C5QrKirE7f77haWTJ0/K0qVLpaamRkaMGCGzZs2Sd999V84999xE7yqwna5QOn2Hs/Igpk3lm5xTb7r9YMkrsM0TdfMhtuFPFyjPiZe303Zi91jr6qlMUyVWW7F5WnFnMaYJPjPtY4i9+/5kiNVdiFOUi9/CtFdSpXmKcoxLCJyRuuq5V/mKrMI1wAzLyrFxH6Ydny2fA7HzxlZBzO/G45k+IQixlkO4nISnw7z8crIy5VebXs50Vv86m/cbSg2S8ReLK4arG0eVVVA6C/E5M4sxBdXSYT7XklowGaSVcLgOHsEnHQD9Uox8zz33yD333KP+38aNG023V61aJatWreqHvSIiIiKnY1NPIiIiciwOdIiIiMixBt06OoOJ3enl1qm8IiKeCNbLtEzF3PuJNvPU45easfbmtd1TIZaJnSLE24F5WYng1PR4C9bt0OnZPRckFedmxpPM9RpunJ0sH1Tiku2Gx96U4kgA63aSDg2fbuWJpNWvaa0zPG3mz3ZUqZvw+PHzn5+Cn8W/NoyHWGsLPmduK54fvibzyeVKwfu5lQ7nbAvRv+zW4xhR/By7UwLmQBQfqyMfzz/1+yQVz4WCDDwn66LmGh1vh7L/yn4MljovXtEhIiIix+JAh4iIiByLAx0iIiJyLNbofIp1bQNt+W1tXQ1XdT3Emq/CNS7Ei/nW4jTzmgWvV0/B+7XimgteJaWe+RGu7yFuZb0Dy+saLHnUoUZbU8nIy4RYepV5Ofbaz+DxLMrC1hG5hbUQq9o6DmJqvjzTXD+krBAj0Wq9se5wpX0OtHocd6bSBqbFfC6M2I91U+6uVIhFxuL3SbOyBv/oAizKO/EPWRBLrzI/r0epGXMpbS1cUa6jM2T4zN/fUeU7J+bHdZyMYjzuF43FJtjBCJ7zwWnmWp78D5R6nFr8PegJBCA2EHhFh4iIiByLAx0iIiJyLKauPsXOpVp3WjbE4vkBiHXk41S+5DTsKDsmzXxJuj2K6bLaFrys7AvhVHKtc7E2kuVU0p7TUhja+6hN9XbFzedCcjV+7KbNxI7V4ThudzQdL0m7osq0UetUzyhOI6Uz06b3ah3NXXXmz7HH0hJCRMTbgcfp7Z3nQGzxrG0Q+2s1TjlPT8V1ChqmmtNqxfVKmk1ZcsKlTDln6mrgqeUTFp52/L3SkZsGMV8yHmMtTdUVw983KVXm7yJXFNP22lIG8RCm5AcCr+gQERGRY3GgQ0RERI7FgQ4RERE5Fmt0PsU6vdydgjUv4sX8ZSwN6zI0Hg/W1fz5qHk6eVcXPn40gLnyrjR7++E6UGFr3+j0tFoN6/kiIpJyAnPSnYXmfHn0nHbYpiWqTOkM4/nXMh7PoaRWXO49u81cm+FWpn7Smdleqt9St+Orx2OcmoLny4zJxyC2p7kQYl43HvdQO54zgTrzdnEffsW7W/AcZd3ewNO+T9QasWZLi4Z0/J4IB/DxO9vw94MvF8/vAxV4/mUGzbdjaVg75NZq15S6nYHAKzpERETkWBzoEBERkWMxdfUp1svUhjLl0qWsjKx1ivWWtNl6ztnF5kvXG7adB9tkl+Mlzfy/4irIRm0DxLTVe63sTp0ezuyumqtpLzRf5vV+jJeQq3Jxldu0JJw2Gg/gOZnShOeHy7K/rgxlmrG2sjOnFJ9RrBVTP9YVYKOZ+J0QTcG/K/dtwGnjWRdhmvHSokMQe/UQflcEzzE/R+YB2ETdfxp4vf2OiWVg6iqKi3DLyMKTELtj5F8hdn/d9RBrLTX/3it8B1NqMqYYQsbBwVE6wSs6RERE5Fgc6BAREZFjcaBDREREjsUandNwZwcgprVZCGcpHcL3Yk3EjV/YALEsr3ka6qb2abBNO872s01bQtzOsuLqVMdhXL+hvR+uJPz4hCZhrU3MZ27bYCitxEemN0NMm17uP6Z0xQ4rU98j5loerW2Bhse9b7g78T3zRLAFhHYuBFvxuDdGsOiiJBtrLurrza0noll4vvjzciEWa8LH4nEfeNp3jFV7MdbxeLC8T5LceDy3tpVCLFKH51/ySfPvuEg21qClbK+GWFyZHj8QeEWHiIiIHIsDHSIiInKsYZu60i7RW1dCNpQVRN2ZmJLyteIlaf+MoK39WLVpkTmgTB/O/AgPUzxVSWEolzljQWW1Sstr5yXqM9PeI23qftYOnOLfNtHc8b6tCI+TlqZq68JjfO3n34PYq9GLITb2hOUYK6krHvfe0b47YsGg6bb35AjYJj4a00/JyoLVbT7sPL1k1gcQ+5d9/wwxb4n5dtZh5W9ZJSXvstllmudM/9KWOJGJo3v1WFp6/EQYU+2+fPxey95kPndTDjXaes7Bcr7wig4RERE5Fgc6RERE5Fgc6BAREZFjDdsaHTs1F9o07NgIzJ83TsXxYtcRzH0+ffISiJ03xbxEduOTY2CbSKYLYl3ZWNPh2Y71IZrBkjcd6uxM0xdRlv5X/rzYV10AsZwsbCPyl4pJEEtW0uUtEzJNtzPrA7CNS6lBY3uAM9M+P1C30xmGbTIOhiDWkYe1PLFU7FS+q6MEYpkB7JBufYaOfKzzSj2g7L9S36d1zqbEUetGtSVOLLc78pQlIZRlCzSbDkyEmC8Z64Ii6ebfQbERWKvq8SpPWm/vd1Ki8YoOERERORYHOkRERORYHOgQERGRYw3bGh07tPxo1I95yK4szKl7izF/PqWwFmLNYfPy3U3n4diz4APMqfv34XLbwqXd+5Vaw6DUZqRVmmu/Yj5cT6XjArxfYzPWgyUn43OGzsXjmXXYfE4a2ZmwTbwOF3FhC4gzs/MeGc1Yj2PkByCWeQTX6m+aho/fHMWavJKsIMQOdZm/0iPpuFS/KLUUWm2WWjPix3YD8bC99iLUczHlM+rxmWsDk1qxHjRSgufV2BQs5vNNxHqcd46Ng1jc0sbGU42PpZ3zgwWv6BAREZFjcaBDREREjjVsU1faZVkr7XJ/V4bSjiEDL+2PyW2CWLYfpwsfbDSnmyI5+Fi+IKYr4spl8Pj2PRCjvqFdsldTV1E8fp0F5vt25OJyAZEuPK+2f/ZJiM1+7yu4Hz5sQRIeYT6/U/fisu5MSfWOnenlWnsQdwyPU/1MnP4t6Zh2OCelBmLv1GOKobPN/Hg+bQUEJb3qScfpwuLVWsoElQekRNGWsDDSzWlMVxzPq/PGnIDY0Y4ciLV2YWoz0oHP6Y6Yn0P7/WPUYGmGel4NgH65ovPEE0/I2LFjJTk5WebMmSPvv//+abd/4YUXZPLkyZKcnCzTpk2TP/3pT/2xm0REROQwCR/oPP/887Js2TJ56KGHZNu2bTJ9+nRZuHCh1NXVqdu/++67cvPNN8sdd9wh27dvl8WLF8vixYtl165did5VIiIicpiED3R+/OMfy9KlS+X222+Xc889V1avXi2pqany9NNPq9v/5Cc/kUWLFsm3vvUtmTJlinzve9+TCy64QP7rv/5L3T4cDksoFDL9EBEREYkkuEYnEonI1q1bZfny5d0xt9stCxYskM2bN6v32bx5syxbtswUW7hwoaxdu1bdfuXKlbJixYoe75utZdzPwrlZmFNP92BuvMOSU8/ao0xd9eIUQPdxvCKGk9ypr2hTaLX8c0zJUyfX5plun5yI08Y7WzFXPq/8FogtLN0LsVeOXwixltHmv2EyjmTDNt4I1hhFq/G8pTOz1lJo9VuuMMay9+H3UOwi/J7Y2Y4tILL8eE76UszPYbjxvDKy8Lw1anAas6HUGXH5gcTR3kdPNrYIaS8xTycPTsDrFTXHivEJorid24+/WzLK8ZxJq7acVzbPg8HSUiahV3QaGhokFotJQYG5j09BQYHU1OhfqDU1NT3afvny5dLc3Nz9U1lZ2Tc7T0REREPekJ915ff7xe9XFsUiIiKiYS+hA53c3FzxeDxSW2u+nF9bWyuFhYXqfQoLC3u0fW9p04W1Dr5WrSOVLr8xvPz38uZZEDtvagXExhebLxmfbMVL1ElNytTgDrxszVVL+5d2WdYTCECsy7KaduYxTDK2nIMXV/99ylqIvR6cBrH4CEyJZGwxp1K89cpKvTyH+oz1PVLT4Nq07g5l1fNNGRD7SxJ2rZ+Wj6ujRyPm501SZq9H8nBlbl8FHnfBrx2eCwmknjM+PIDW6eQeXI1A3B78jhlfgml1nxvPv70NYyCWs9vynCfxuy8+iL87Epq68vl8MmvWLFm/fn13LB6Py/r166WsrEy9T1lZmWl7EZF169adcnsiIiKiU0l46mrZsmVy6623yoUXXiizZ8+Wxx57TNra2uT2228XEZFbbrlFRo4cKStXrhQRkW984xty2WWXyaOPPirXXHON/OY3v5EPP/xQnnrqqUTvKhERETlMwgc6N954o9TX18uDDz4oNTU1MmPGDHnttde6C44rKirE7f77haW5c+fKc889J9/97nflgQcekIkTJ8ratWtl6tSpid5VIiIicph+KUa+55575J577lH/b+PGjRBbsmSJLFmypM+eX8t9qtOFk8zTLl1NWNfgC+E0XU87ZgADE7Br+N7jWGfkOmHOa2b5sD2ANpXPpSzPjvcUESwLoASKt7RAzNNmrqFJalGOp9LG4f5d10PssWnPQ+zDwtEQa8k3n2uZI3BKsbsWpxRT31CX7m/Gc8PXFIBYWzHWZUxSOpXvrCvCx0s21wt2zmiHbcIHcPKGPwvrgrRzmdPLE8edgh3qtU7z7YWW+js8xKrpI7AtxO92zIRY9l78fZZcY25fZCjfHWqN6yD5/cOmnkRERORYHOgQERGRYw35dXTs0C4ja6wdh42p42GbrjQcG7oimHbIScVO5dPzqiC2sdFce5RehZeB3RGcvq4ZLKtQDmfauea2TMXsnKSkkYJ4DmWPxGvSPzvxOYhVVWM6NdWy+HIsDffLk6J0ZO/Cc4jpip7TVkZ2p+Fx0niVGbl7PsIpv/de8WeIPXt0tul2/QHsWB3OUpLcyirZPMaJoy3joErGNGN6lXk+eeNUTHVeNvEAxD5swKVLtGnoMWVJAlez5XtB+e6INWG5xmDBKzpERETkWBzoEBERkWNxoENERESONSxqdLSp5FqO1J1ryaFXN8E24c9g5+lLL9sJseNtAVsxw2OuzYhk4NgzzZofFZG4snw/aykGnlabIe3m2q+UetzG24qJ8aP7cTmCq+eth1jxtCDEXj9krtVwxTAXb3Rh7Zc7OwAxraM5zzUz6/uhvRfRKpze6yrCGhpfEB+/czQ+3m8rsc1MW9h8HhkBPNfSq5QanUylbqxJqeFSzu/hfNyt7H4u1Bou5Rho6s831+0Yxdiro7ULj93I9GaIZfnx98ihcqxNNdItU98bsR7Hbi3sQOAVHSIiInIsDnSIiIjIsTjQISIiIscaFjU6GjXXHDDnSA0fvj1xZY2BHB+udxKJ433fKZ8EseRac07X12JvzRzmygcn9Rj4rEu24zFOqcUTK3katiDJ8uC59tNdl0Es3ZK294ZwLXZXQS7EjNoGiNGZ2fnseSdNwGB9EEKuOLZjcAfx+yRJqdu5c9Lbpts/+ctVsI0nHIGY4VXWB1OW9NfqHe3UJw0Xdl+7up3S1qdjTABifsuaW7him8iJtkyIFafh98nOimKIBbDzh4jlGFvXnDuVwXIu8IoOERERORYHOkRERORYwzZ1pbF2K4+cOxK2Sa3Gpfr/UoEpqdIRODXdX4dTD5MbzbfbCvGQpG1VpiIqnW61y4mD5dLhcKFNL5WojbSGciW4rQvPhd9VXwCxFRf/AWIr0xaZH+sAXspOrsU0hLcVU2MupYs1z6vTU8+DpiDGfJiy1NrAtMzBdNPkQC3E6iLm4+zuwL9lIwGcBuwN4bnmUlIpXFYggdx4rAw3LgXQXmSOTRuNyxZcnYdLnqzafTnErN3uRURc2uEc4seYV3SIiIjIsTjQISIiIsfiQIeIiIgca9jW6KjLclvaKnSl49vTmYM505ZanA4qWo3ODFw2O/7mCNPt7P04DThemA0xl9ICggaeO0M5Fyw8bVhv4YphzZV7Gz5W+2U4mfTFWqzbKcwy15s1jcXHSi/Hc9RoU2p0lKXdWZdxeup7prTccKXicQ9nKbUajX6I7WzCqcF3jv2r6bZ/Mi7737VDmb4ewiIxrc0Mj3vPaXVNWo2leHG75nF4HnXmmOtEK0MB2OZ3UfxO8HiwDYy8lwWhEfvwO8A6vVw7v7WlBwYLXtEhIiIix+JAh4iIiBxr2KautMuJrhRzx1dXDKeSd+RjbPa5hyE2LbPK1n6cPGS+jBz3YGrM3Y6pjhin/A5KceW4uKPmlIUrgF2K06oxrREaj5eHu2J43l5fsA1iTxyaZ7odw2bGEssLQMxjYyq8iIgM4svUg4F2Gd+TonQDb8bzxduJ3zHusJIy78R01rbWMabbHYdxWYF0m3/eaisjuwVfw2BOWSRab1eFtv6uEREJT8iHWAwPscRTzSmokswgbDMxox5iB3eOgljRYUxnqauoN7ea92uIHXNe0SEiIiLH4kCHiIiIHIsDHSIiInKsYVujYyeX6m/AXGVaVRrEtr07EWL1M3G7Y3uKIJafZc69Z+1XegGEWiHEFhCDk/Z+W4+Lp1PpPO/FGoxIAB+rsRnPq19VzYFYp6V9RHsx5uLDuVgn4JcciLkP4nRTt18p+lEM1/oNrQYwFgxCzHPeORCLZOLfnyl1+Bwls/Hx3rC0o0luxPMquRHPv9gIrBtzHa+G2HA9nqdi/bzbbZGhLTUQTVHqRpUZ4VZj0xoh9uLu6RDLGItLDUS3Bs78BCIi7ebvMO3zP5jPDV7RISIiIsfiQIeIiIgciwMdIiIicqxhUaOjrpmjLdFuWfLco61fk4w1EmmTghDriuNzJhVhrYO/2Zzr7CzE2hv/hzUQs5sLpkEohGunuKK5EMvdise4YRYurNGYgedkRrK5vqwVl2aRrgx8fFHWSUk9hl8T8Ralzojn32lpdQ0upV4r8wh+79TM8UFsX3UBPkml+fsjXoBFHv4GpbWDF//mdStrvXD9pJ7TvqvjE0ogptWExs7Dz97YCebfB5dl7oNtGsZjzdXeNeficzYrn1kbn2MjiuftYMYrOkRERORYHOgQERGRYw2L1JU6vU+JeUeZp38rV/slowIvBdcVYQfYkBfvnX5U6WIbMV+mTjmEUwWNQABi2lRVGpysSwHEmrCLfepepYv1iJEQ87Tj3ybhLvwYf7bokOn2H1y4vHzTZHys0udqIUZ9Q5t+61bSmIYb05heZdWJrko8Z2IZ5u+n5BNKiltJU3nrQxDT2syQmZaWsnJnYLf4aDLer7UEU4Wt43EaeluF+bP8+/QLYZsPKkdDzIsrR8iIfZiCcrUqy5QoneyHEl7RISIiIsfiQIeIiIgciwMdIiIicqyE1ug0NTXJ17/+dXnllVfE7XbLDTfcID/5yU8kPR2nvv3NvHnzZNOmTabYV7/6VVm9enWv98NOHlVExKg1t7Y3JmKe038Sa3sCe3GqenCKMq0ziM8Zsy773W4vP8rp5UNHrNXcwsOjnP/hCVhDY7hx+f5YqnJeJWEe3+rc6ccgdmhDKcTq5mObkpxy3F+PF786otW4DMJwZfezqLUC0KZ/+4P4HSOitA2JmL8XkrHkT203Yv3uE9HbzFjPZTLTli2RkfjZjgRwuYC6C/G4eFrwe/62Kzaabh/pwOKb6IlUiGUfsdFPQkTiNXguWA213zUJHeh86Utfkurqalm3bp10dXXJ7bffLnfeeac899xzp73f0qVL5eGHH+6+nZqKB42IiIjoTBI20Nm7d6+89tpr8sEHH8iFF35SFf7444/L1VdfLY888ogUFxef8r6pqalSWFiYqF0jIiKiYSJhA53NmzdLIBDoHuSIiCxYsEDcbrds2bJFvvCFL5zyvs8++6z86le/ksLCQrn22mvl3/7t3055VSccDks4/PcVJUMhnCapXWZTu69au0wrKyMnteIlx5R6LHXqSsVLjqn1uB/JtWeetudKwsM01FamHM6s55rWZd63ZT/EXLOws3VrJZ63JecHIbaqaKvp9mcb8PMW17Ihikg2rsacfBDPWy0lZ32tQ+2Sd6Jpn21PNeabMpRjkNSG923xmL+LMqrw/fa04feaKKnIuM3p5dY0ulOPsfY7A7bJzcZgENN9viD+PkutxljWFWdOB2d68bOYXqp0Kv8ogPtxVElTZeLnOFrfcMb9GMwSNtCpqamR/HxzbtLr9Up2drbU1Jz64H3xi1+UMWPGSHFxsezYsUO+853vyP79++XFF19Ut1+5cqWsWLGiT/ediIiInKHHA537779ffvCDH5x2m7179/Z6h+68887uf0+bNk2Kiork8ssvl0OHDsn48eNh++XLl8uyZcu6b4dCISkpwT4iRERENPz0eKBz3333yW233XbabcaNGyeFhYVSV1dnikejUWlqaupR/c2cOXNEROTgwYPqQMfv94vfj5d1iYiIiHo80MnLy5O8vLwzbldWVibBYFC2bt0qs2bNEhGRN998U+LxePfgxY7y8nIRESkqwmmvZ0OrcYHplCHMrbpGYKfopFbMSeftOPOUXxERb3WTeb+yM2Eboylo67FocOptPVVSA55/7gjWCRz5Pf4BcO0XFpluT8s+AdtUfgZrhfamj4GYJ4LFPClpWE9gNNurjxuu1C7W2nfMJJz2r3ErXzHWFjXuME4pdoeUJSyi+GBa64Lh3HpGbeFhqdsx2tphm9hk/EzVz8TPT/soPFa+LvzsbWkaa7q962PMYKQdxV/tuRXYHV2S8SKB0Ygtaoa6hC0YOGXKFFm0aJEsXbpU3n//fXnnnXfknnvukZtuuql7xlVVVZVMnjxZ3n//fREROXTokHzve9+TrVu3ytGjR+UPf/iD3HLLLXLppZfK+eefn6hdJSIiIodK6MrIzz77rEyePFkuv/xyufrqq+Uzn/mMPPXUU93/39XVJfv375f29k9GwT6fT9544w258sorZfLkyXLffffJDTfcIK+88koid5OIiIgcKqELBmZnZ592ccCxY8eKYfy9y3dJSQmsityfYHp5CqYJvMeU7s5jCiDkiuJlSM+xarxvdsB8+ximGDRMCQwddo6VtkKuu1NJr9rLiEpbl3kZhHEpOD20sm0ExOIZuK9Jrfj3UKwIV2N1K5ftrVPOtdeppQScyO5n1l2NU35TGnG6sE87BhHz++tqw/dWWwVZW/JA29/hsiK73dcJywMoqyCHxmOaKhzA59SmhGf68fjNyT5qul1djOUOnUfx3HBFDYiJUhbhxNWv2euKiIiIHIsDHSIiInIsDnSIiIjIsRJaozOY2ck/x5pwmp2nEOtxvCfbIBbLwM6/RhFOy3dZ8vF2c+XkLForAK2Tfc4enCJasQCniKbHzeey3431PqlejBkx7KDcPA7rFbztWL/mF2Whzu17LI/Pc/nTtG7XRodSs6TE3AFcqt/wmc8jlzJt3KXUHoryvaMZLsdPe512Wpy0nostIGLYNUiSlDKYSBd+B6QlYbuOUMz8u6WpAZcByFAOZ9JHhzGYhffVXudQr9vhFR0iIiJyLA50iIiIyLE40CEiIiLHGrY1Onao+WhlyWxtLQLX0eMYs66ZI1gHNJzXrhjOtPVltLoMXz2uVZN2As/JqmLzGjkl4xphm5sLMLazBluttE7Fxy/4QGkt0IH1Q65RxabbRnMLbCNKHclQrwmwS2sPoq0rpNVNSFUdhKDCSqn90mqA1FqhYfwdo33narVN7nHmurSkFnzPoil4DFrG43aLxhzA7aL4nL/76ALTbf8xLAJKacQ1c1xq2xb8PGqfPev7MdTODV7RISIiIsfiQIeIiIgci6mrT7FzOS5+XGnjoNAucxotw+NyPPUNdRpwGFMdyU14mbq5xfzRXvbujbDNj+c+D7FJ+ZgO2dlZDLHgRFw+IcuLU9P9R837q6VNhksLCDUdoqSMNGpqU4tZxFu0JQTspR2cmDLXXpPGnaIsD6K8310jzNt1ZeDjt47Cz8XYSTUQm5leAbFHdy6AWHK6OUWcfhTPoYxjmEY2lOUI4jXY0sjakV1k6H9GeUWHiIiIHIsDHSIiInIsDnSIiIjIsVij0we0KaKGkhvv7RTOoZ4XpzPTcuAupabL5cNzyNuBU72z9ptrBZpn4Db/e8+1ECvIwOmm6Rm4byeVKecuZa37vKPm2+5cXCLfUOoEnHjOa6/J7us8m7oaO9uczb4Nddr3staKxzN2NMR8R80tfBrOxzYonfn42dNaOxzoLIRYXhZ+B1QdMLcSKqjG2iHfcWxfZDRiTDPU63E0vKJDREREjsWBDhERETkWU1c9pKWptKmI2uqSvb0U7MRpnnRm6uXzCE4RzdiPl6TdMfPKyB0FeHm+LRXTT4fbMTa5CFNLuwNpEGsdpazQOjXffLsKV3Z2NTRBTOPEc/5sPtt27uvEqcJnQ+1KXliAG3qVFGAnTtluvHyM6XZoAqapzp95BGJj03BV8t+9dxHEMotDEEs9Yd43bwf+TpJ2/O6Ityirkg8TvKJDREREjsWBDhERETkWBzpERETkWKzR6SEtx5voTstOrE2g3tFaKEQKR0LM02E+Z1JqseagOQfrceZfsBtiHwfzIfYvF26E2E87LodYx3Hz31Lp+7B2wNBaXeBm+vIMSs0cbDOIPz9ns2+2WtZoyxY4sObPbmsHT34eBuNYVxMZi1O9O/KxBq1huvm2EcDz8Tslf4LYzk6chv7XkvEQa2rEmrwsS6mNp035DETP3B5kOOEVHSIiInIsDnSIiIjIsZi6ShAnXh6m/qWdL9oUUf++aoh1nGfuOJ7Uih3Os8vxHP1r8HyILViwDWJP7yuDmPgwBdBeZH6OpgsxdZC1H5dncO/DKbla92h+pnpuqL1ndr5L7b6meFMQH38CrnisaS3G6wJeS4o1OYDLJ7wYxGnjfhemm5pqM3G745iuTa8yv1bPSWUpEyXFrRlq50Jv8YoOERERORYHOkRERORYHOgQERGRY7FGJ0GGS+6TEkdbvl+bTh1XWiikWGaJd6XjdNa2Qqx9SD0X20mMS2nAnRuDIc26xpmm254OF2zj7UiFWMY+e4+vvUeuJPPXmtZKg5/PocPOsbJzHoiIyDj8HMR9uF3Mj9cADOWyQDTFXPtWkIE1dH8+OgVixVnNEPOk4mc7pRZrdJIbLdsprSnsnvPDpZaUV3SIiIjIsTjQISIiIsfiQIeIiIgcizU6RIOUtny/7VqEdKx7scqowlz8iXZsC/GLA3MgNqOwCmIHTuIaOa6YuSYnjqvoq7VCKTMmQMwbwloEV3U9xOIh87oiTqw5GE7stHdwKW1EtM8FriYl0lmI6zg1nI/3bT8nAjGP37y206i0IGxzrCEbYodrcyHm34mf2YLN+HiumOVVKC0s7Bounw1e0SEiIiLH4kCHiIiIHCthqavvf//78sc//lHKy8vF5/NJMBg8430Mw5CHHnpI/vu//1uCwaBccskl8rOf/UwmTpyYqN0kGlLUbtRRvLTv9punpWbuboRtDD9OXU3aNwJiHSPw0vjYUny81i5Me6VeaL7cf/QEXrJ3BXE/UusxFeFLxa+rFCV1ZeVJxw7QsVZcNp8GJy29Yk3hurIybD2WNpW8uRRj0Rl4frgj+Dn7f6ZvMN0ub8Hp63lZ+FhVhzHN622DkLjrghCzLiehpeOGS0rKroRd0YlEIrJkyRK56667bN/nhz/8ofznf/6nrF69WrZs2SJpaWmycOFC6ey017eDiIiI6NMSdkVnxYoVIiKyZs0aW9sbhiGPPfaYfPe735XrrrtORER++ctfSkFBgaxdu1ZuuummRO0qEREROdSgqdE5cuSI1NTUyIIFC7pjWVlZMmfOHNm8efMp7xcOhyUUCpl+iIiIiEQG0fTympoaEREpKCgwxQsKCrr/T7Ny5cruq0dEw5GWj4/t2Gu6rU1Ld2di7UrOzkyIVV2DNTqvHJsKsfsmvwGxh1/+R9Pt7HOxXUVTFJ+z5mL8asr/UKnRycTaDFeHOdVtdEVhG7uGyxL5dqZwi/Tta7f73noCAbxvmnkqtpGOU8Q7R+J5FUvB53QpL2lMLp6nPjduWBcxP8fuxiLYpukjrMfJVH6ljdiPLSCMAH5G3RFz3Zt1OQVCPbqic//994vL5Trtz759NpvU9JHly5dLc3Nz909lZWW/Pj8RERENXj26onPffffJbbfddtptxo0b16sdKSwsFBGR2tpaKSr6+6i4trZWZsyYccr7+f1+8ftxtgcRERFRjwY6eXl5kpeHl+H6QmlpqRQWFsr69eu7BzahUEi2bNnSo5lbRKRMv9VWT1bSPil1uPpw5g5csbXzEpwSvjE4CWKXXrbTdPtIKAe2aW7ElWMzjmKX8+AEjCW150Ms3dLN2ZWMfwi5anFaupbi0t43rTM0PNYQS2/15f5qaVK7PPk2f79YjmmkEM/ljjw8R1tH4TmUdhmeC1qaal91AcQqmwOm223KyuIe5a11KYsZJ7Vi6kqq6pQNzeekEVXuRyYJK0auqKiQ8vJyqaiokFgsJuXl5VJeXi6tn1q/YvLkyfLSSy+JiIjL5ZJ7771X/v3f/13+8Ic/yM6dO+WWW26R4uJiWbx4caJ2k4iIiBwsYcXIDz74oPziF7/ovj1z5kwREdmwYYPMmzdPRET2798vzc3N3dt8+9vflra2NrnzzjslGAzKZz7zGXnttdckObn3fyEQERHR8JWwgc6aNWvOuIaOYZjXdHS5XPLwww/Lww8/nKjdIiIiomFk0EwvJ6K+Y83bqzU67Vhr4lKmeqefwIKC8F6c9tpUmAaxsWnmVhFJSsHC5NlHIbavBOsh0t7DWqHQaJwu3JU62nQ7sAOnCruy8HVKWzuEjA6l5YYXaz/sbKO177DL7lRsa32MVr/R1/VD1n3TzjWt5YZWyxNvCuJ2Y0fh42WYp5PHkrAKoysd63E68rFhwmfzcKZuWeYhiP0qPgdiAZ/5mH54dAJsk/MxPmfmYfzseUNYHxdvaYEY9dygWTCQiIiIqK9xoENERESOxdQV0TCgpQ48Xvz4e9ojEEtq8UEs8DGmBXal4BpaE68wT929Ih8XFF1XNxli357xF4j9R93nIWZ4MC2QXmW+HTwfp68nN2LqKuVYEGKuqLLq9NEKiLlTLKkU5f0+m9WH7aTLRDA9Zneq99mks6z75srAtKZXiWlpKtekUoi1F+N9uzLM72VoDP7NHpmNx2BW8QmIvV87BmLhuL1fjR9uNaeq/E3KfuDuiyuK562rzV5qc6gtXTAY8IoOERERORYHOkRERORYHOgQERGRY7FGh8iB7OTxtanT7masa0hVagc8XVj30jAda3l+/9Zs83NmYZuFF+c9AbFfNF4CsRHjT0IsGsO/1U5clmW6nXFY+3tO++oLQMRfj1POPYU49d1oDpm3SVcKM5SaKO0YaG0QtHoWT/YIiMWazO+RNqVdq9uxWwPkSlHua4lpr0lGY1fvrnNHQsynvN/hEVjb1HC++ZimnovnRroHl0X4Qv423I88fPzfVl8IscO1uRBzhy21akprh8BBrHvzHqvFDaP42dBYa71Ys3NmvKJDREREjsWBDhERETkWBzpERETkWKzRIRqmtPoN4ziuM6LVgvjqsVYj9yOs0an9nKV+oBXrIb5z+AaI/fPILRBLGokFEH+pmASxeIp5uxZc3keS2vBvvPZ83P/cGK534vFjPYvH0j5CbbmRmQEhl9KGQ61xURhdWNNhp9bGnYn1Q+pzKjVFrgKsUzG85vcyloE1THE/Hve4B9diaivF9Y0imXisfJbOCKETeL+xE2ogpvnJgc9BLBjEdiMZW1IgllZjPtdSq/F4JlU2Qkyrx4mHsD5OXVPJ5npM9He8okNERESOxYEOERERORZTV0TDlDbNWOt2raYwqushlj4CL+03HzU/R2cxXoo/dAKnU/9KsFP0+IwGiHWGMVWTU9xsvl8AUwc7gudALElpFN06EtNZ3g58zuS08XhnC/+eKojFRuJr95xU2kdoD6hNRx6Zb7rprsPO7e2zsM1Ccm0b7lsavnZPG06Vbp5sThul1uE51JmN79mJz2Fa0J2hTMX24Tmz5JztptvvN2Ebh1Qv7scTR+ZBLCcVX3v7BzkQS2nA1GnE0iE9qx6Pnfb5UdNU2mdPwenkPccrOkRERORYHOgQERGRY3GgQ0RERI7FGh2iYUqbXq7Rph670nD6re8o1u3kZ5iX/o/uw7+t6i7A2p5DgrUrqaOxhuGikgqIvb3TXH9z09itsE3xPwQhtuH4RIgFd+LU+uR6rJhxxcxfpTE/bhPOwjqStkKcKpyzG+tZkhqUmg4fvm9x69T3DGyzkNRqrxakI98PMX8z7q+1JufEZ/B+kUysxzGUFg1xpaVHeiq2hRjtN9ddNWakwTZT0qoh9uzR2RCrfHc0xLRWDm6lJCqn3FwPJlGsnzFq8XOh1eP0dio5a3bOjFd0iIiIyLE40CEiIiLHYuqKiE7LbupKk7bXfNk+nonplppLcFXb5D243S4PdsC+aQqmpfZappfXdeHjv6GsqJzmx6nNMg27Yre24LT8SMC8v1qqxhfCvyvjykLGySdxWrcUY7f45EbMpZycZL5vVDlMWYeVlJEXX5OmaRKmpWKWkPaaxs+shFhzJx7jzi78lfS5kQcgVttl7lD/l/1TYJtt2SUQq6vLgliqksnL2Y3pIE8Hxtwhy0rIIVyjQFvB2m66iWmpvsErOkRERORYHOgQERGRY3GgQ0RERI7FGh0iOi1tKmy8AVsLaB27XQXmaeJQ0yAihe9gDU0YSykk9jEWnPyqpQxiD3z2VdPt90LYvvzx6b+G2E+OXwGx75T8CWJP1FwOsfI08zTueBPua7wkDDH5GDuJ1yzEmo6srVi303A+xqxdvZOVxtl1F+Dft2nYtF7iym+H0BTcN2vLjcwkPF9+OuF5iN24439BLCMZ36NXD50Hsf9v+kvmfY3NhW3a38QlCmQs1rxo75E7jHVMKbuVN8kiWo9tSjTatHHW4yQOr+gQERGRY3GgQ0RERI7F1BURdTubS+paistjnW6bilOKM/c3QyyWimmZmA9Xv+3Kwr/VfvCn60y3R5+PKYe3U3B6eWUoALFVVQshdnPBexCr6zCnoJq1qeqKhnRl/ncrHoPo5/A98m7D/J51qnfTDOXYpWL6qTmAvwoumnYIYvUdeAySPObnyE/BVZy1NNVVo/ZA7JKMjyH2bB2mJ//3nmtNt11BnNPuUl56yTqM+YK4fIK6EnVzCJ8jxd60fHgsm58pprj6Bq/oEBERkWNxoENERESOxYEOERERORZrdIioW1/n/63L37uU7s4aVxSn92ZWYG1JXJnSHrHMVj9xEmtZftc+A2L/WFoOsb82jIfYr2svhljAd+ZO8Odm1UDsRCa2mHh/D06HD6TitPzU+Vi309hurqFZVIR1NptrSyFmrbMREUlPwtfk8+AxiFg6t1+TswO2ae26CGLv1OPr/PXGS/A5R7VBrLPeXNuU3KBMma/Bc8gXxDoyXw3W3khTEELY1AOnk9vpNi7CFhD9jVd0iIiIyLE40CEiIiLHSthA5/vf/77MnTtXUlNTJRAI2LrPbbfdJi6Xy/SzaNGiRO0iEREROVzCanQikYgsWbJEysrK5H/+539s32/RokXyzDPPdN/2+/2n2ZqIBoJWO+D245oi8RbLOjrW2yLiiWMthTsb20L4glj/kHXUBbFIuvnvt5YYtlloC2DFxa92zofYmLkVENtZUwSxnAxzHcndpRthm5fqLoDYZ7IPQKx4VhBimq/k/hWfo9n8HEc6cmCbaBz/vs1KxhqgrTUlELt/yusQq49mmG4/dhBbZDQ24/o74wqUdgm52ALCvS0DYvnHzccvtQbXLUquwdqejmI8F3whXDNHE+/A98hak8OamsEpYQOdFStWiIjImjVrenQ/v98vhYWFCdgjIiIiGm4GXY3Oxo0bJT8/XyZNmiR33XWXNDYqHdc+JRwOSygUMv0QERERiQyy6eWLFi2S66+/XkpLS+XQoUPywAMPyFVXXSWbN28Wzymm7a1cubL76hERDZx4GKcj27m0b7Rg6sAVxWnMSa2YOhDB1EySpYVCGs7qlvZ8bBnQcD6mwdq7sBVFRxO2sUjNNk8Tf/bEHNjmQB12087x42tP92IaZmIyvohlh5ZAzDqFfUnuB7DN0vxNEPvnd74CsdEF2KF+SwtOCW+ImNNBdfVK2jEFp3Uf2jEKYq4IHgO30k0jsMf8vsXSlBYQjTj9PrUNz1GjA2OQcqUhrUdXdO6//34oFrb+7Nu3r9c7c9NNN8nnP/95mTZtmixevFheffVV+eCDD2Tjxo2nvM/y5culubm5+6eysrLXz09ERETO0qMrOvfdd5/cdtttp91m3Dgc8ffWuHHjJDc3Vw4ePCiXX44FbiKf1PSwYJmIiIg0PRro5OXlSV4eXn5NlOPHj0tjY6MUFeEsByIiIqIzSViNTkVFhTQ1NUlFRYXEYjEpLy8XEZEJEyZIevonOd3JkyfLypUr5Qtf+IK0trbKihUr5IYbbpDCwkI5dOiQfPvb35YJEybIwoULE7WbRDTAtGm7RivWrmjT15OUlhLREvMfY+4w1oekh/F+nSNSIRaqwxmgrnPwvrUt5inQrS24r64TGNvongixwiycULE3qQBiWv3QGxWTTLdPdARgG00g0A6xpnZ8P16pnYZ3bjBfUU9WWjakrsMp4t4OnOIfV0ox896tg1h0hHm6uqcNj7GqE6evn009DqeTDw0JG+g8+OCD8otf/KL79syZM0VEZMOGDTJv3jwREdm/f780N39SMObxeGTHjh3yi1/8QoLBoBQXF8uVV14p3/ve95iaIiIiol5J2EBnzZo1Z1xDxzD+PqJPSUmR11/HxaiIiIiIemtQTS8nImexXtrXujvb7uQcVdITEZx77Ppgl+l2rExJtygyqnA/klpxmnvgMH5thkaPMN3Oa8K0TFsRTp02GjGlc/QcTHFJBCfI5o/B6d8et3mV6f2NWFPZchS7uSe14ONnHMXdSMPdFZdlYetWwVWQfcpLyjiG760viMfT8OPUcW+9Ob1nNGIXeElRVupWupK7vPj42lIJNHQNugUDiYiIiPoKBzpERETkWBzoEBERkWOxRoeI+o3dehy7tTxxpfO0tebCW34Qn2BMMYSSqrFlQGwEdruOpuC+5ewxT1tuPBdniqbWYd1OJB3rdjJ3YM1IBLsqSOgE1t+4LW9RZw4+p0dps5B5SNsOnzOjAuukDLf58dJP4PuTfgSncLsjWKPjUpYCEOUYR6vNrS60pQfU1g42a2/OppaMBh9e0SEiIiLH4kCHiIiIHIupKyIadM4mTWCdhu7ODsA28YMVEHPlZkPMo8xa9iVjWiOpstF0O5CEKyobXkwZJeMMcYmm4N+fnqP2VhG2ivnxOUWUNFUYY9Y0mAimqUREoqnm/c3a0YD3q8WYK2eEsl09xLSVsz2BgHkbZXXjszqHmKZyFF7RISIiIsfiQIeIiIgciwMdIiIicizW6BDRoGN7ermN6cLWqcgip5iO3Ixdw40GnALtVepIJMs8/zt1Lz5nPBv7J7hDWH8Sz0zBmNIGwdOC941lmO/risZhG8+xaoiJDzuhSzp2L9faMbgaLdPyvXjsXGnKY2nvo0Jr0WCdOs6aGjodXtEhIiIix+JAh4iIiByLAx0iIiJyLNboENGg05c1F1o9jlbb44pibYk7BetljC6s24lVnTDd9ubl4mO1Y02N9lhupcZFq+WRqFI/1GSuMzKyld4RXuVrX3lOieIxcHWGcTsLoxEXH3KlKDVRymu3roEkop8LWg0X0anwig4RERE5Fgc6RERE5FhMXRGRo6lpKpupj1grds62th8QEXFZUi7Remx5oN1PpXTr1tJUaurH2rG7KYjbKE/pSsJfBVq6SXtdWmrQKt6h9NI4C5xOTj3BKzpERETkWBzoEBERkWNxoENERESOxRodIhp2zqbGIxYMQsxa86O2sLDWz5xCvKWlV/slgtPhbU/h1mJKfZLaOkO5L2zDmhoaQLyiQ0RERI7FgQ4RERE5FlNXRERnyU5qxu4Kv2eznTYdvrePpWGaioYiXtEhIiIix+JAh4iIiByLAx0iIiJyLNboEBENENu1MTa3s9bfaPez+1h2u74TDXa8okNERESOxYEOEREROdawSF2ti78w0LtAREREA4BXdIiIiMixONAhIiIix+JAh4iIiBwrYQOdo0ePyh133CGlpaWSkpIi48ePl4ceekgikchp79fZ2Sl333235OTkSHp6utxwww1SW1ubqN0kIiIiB0vYQGffvn0Sj8flySeflN27d8uqVatk9erV8sADD5z2ft/85jfllVdekRdeeEE2bdokJ06ckOuvvz5Ru0lEREQO5jIMw+ivJ/vRj34kP/vZz+Tw4cPq/zc3N0teXp4899xz8o//+I8i8smAacqUKbJ582a5+OKL4T7hcFjC4bDpMUaPHi2VlZWSmZmZmBdCREREfSoUCklJSYkEg0HJysrqs8ft1+nlzc3Nkp2dfcr/37p1q3R1dcmCBQu6Y5MnT5bRo0efcqCzcuVKWbFiBcRLSkr6ZqeJiIio3zQ2Ng7Ngc7Bgwfl8ccfl0ceeeSU29TU1IjP55NAIGCKFxQUSE1NjXqf5cuXy7Jly7pvB4NBGTNmjFRUVPTpGzXY/W0kPNyuZPF183UPB3zdfN3Dwd8yMqe7INIbPR7o3H///fKDH/zgtNvs3btXJk+e3H27qqpKFi1aJEuWLJGlS5f2fC9Pw+/3i9/vh3hWVtawOkH+JjMzk697GOHrHl74uoeX4fq63e6+LR/u8UDnvvvuk9tuu+2024wbN6773ydOnJD58+fL3Llz5amnnjrt/QoLCyUSiUgwGDRd1amtrZXCwsKe7ioRERENcz0e6OTl5UleXp6tbauqqmT+/Pkya9YseeaZZ844Sps1a5YkJSXJ+vXr5YYbbhARkf3790tFRYWUlZX1dFeJiIhomEvY9PKqqiqZN2+ejB49Wh555BGpr6+XmpoaU61NVVWVTJ48Wd5//30R+STddMcdd8iyZctkw4YNsnXrVrn99tulrKxMLUTW+P1+eeihh9R0lpPxdfN1Dwd83XzdwwFfd9++7oRNL1+zZo3cfvvt6v/97SmPHj0qpaWlsmHDBpk3b56IfLJg4H333Se//vWvJRwOy8KFC+WnP/0pU1dERETUY/26jg4RERFRf2KvKyIiInIsDnSIiIjIsTjQISIiIsfiQIeIiIgca8gPdI4ePSp33HGHlJaWSkpKiowfP14eeughiUQip71fZ2en3H333ZKTkyPp6elyww03SG1tbT/tdd/4/ve/L3PnzpXU1FRom3Eqt912m7hcLtPPokWLErujfaw3r9swDHnwwQelqKhIUlJSZMGCBXLgwIHE7mgfa2pqki996UuSmZkpgUBA7rjjDmltbT3tfebNmwfH+2tf+1o/7XHvPPHEEzJ27FhJTk6WOXPmdC8/cSovvPCCTJ48WZKTk2XatGnypz/9qZ/2tG/15HWvWbMGjmtycnI/7m3feOutt+Taa6+V4uJicblcsnbt2jPeZ+PGjXLBBReI3++XCRMmyJo1axK+n32tp69748aNcLxdLtcpWyMNRitXrpSLLrpIMjIyJD8/XxYvXiz79+8/4/364vM95Ac6+/btk3g8Lk8++aTs3r1bVq1aJatXr5YHHnjgtPf75je/Ka+88oq88MILsmnTJjlx4oRcf/31/bTXfSMSiciSJUvkrrvu6tH9Fi1aJNXV1d0/v/71rxO0h4nRm9f9wx/+UP7zP/9TVq9eLVu2bJG0tDRZuHChdHZ2JnBP+9aXvvQl2b17t6xbt05effVVeeutt+TOO+884/2WLl1qOt4//OEP+2Fve+f555+XZcuWyUMPPSTbtm2T6dOny8KFC6Wurk7d/t1335Wbb75Z7rjjDtm+fbssXrxYFi9eLLt27ernPT87PX3dIp+0B/j0cT127Fg/7nHfaGtrk+nTp8sTTzxha/sjR47INddcI/Pnz5fy8nK599575Stf+Yq8/vrrCd7TvtXT1/03+/fvNx3z/Pz8BO1h39u0aZPcfffd8t5778m6deukq6tLrrzySmlrazvlffrs82040A9/+EOjtLT0lP8fDAaNpKQk44UXXuiO7d271xARY/Pmzf2xi33qmWeeMbKysmxte+uttxrXXXddQvenv9h93fF43CgsLDR+9KMfdceCwaDh9/uNX//61wncw76zZ88eQ0SMDz74oDv25z//2XC5XEZVVdUp73fZZZcZ3/jGN/phD/vG7Nmzjbvvvrv7diwWM4qLi42VK1eq2//TP/2Tcc0115hic+bMMb761a8mdD/7Wk9fd08+80OFiBgvvfTSabf59re/bZx33nmm2I033mgsXLgwgXuWWHZe94YNGwwRMU6ePNkv+9Qf6urqDBExNm3adMpt+urzPeSv6Giam5tP2/1069at0tXVJQsWLOiOTZ48WUaPHi2bN2/uj10cUBs3bpT8/HyZNGmS3HXXXdLY2DjQu5RQR44ckZqaGtPxzsrKkjlz5gyZ471582YJBAJy4YUXdscWLFggbrdbtmzZctr7Pvvss5KbmytTp06V5cuXS3t7e6J3t1cikYhs3brVdJzcbrcsWLDglMdp8+bNpu1FRBYuXDhkjqtI7163iEhra6uMGTNGSkpK5LrrrpPdu3f3x+4OKCcc77MxY8YMKSoqkiuuuELeeeedgd6ds9Lc3Cwictrf1X11vHvc62qwO3jwoDz++OPyyCOPnHKbmpoa8fl8UN9RUFAwpHKevbFo0SK5/vrrpbS0VA4dOiQPPPCAXHXVVbJ582bxeDwDvXsJ8bdjWlBQYIoPpeNdU1MDl6m9Xq9kZ2ef9jV88YtflDFjxkhxcbHs2LFDvvOd78j+/fvlxRdfTPQu91hDQ4PEYjH1OO3bt0+9T01NzZA+riK9e92TJk2Sp59+Ws4//3xpbm6WRx55RObOnSu7d++WUaNG9cduD4hTHe9QKCQdHR2SkpIyQHuWWEVFRbJ69Wq58MILJRwOy89//nOZN2+ebNmyRS644IKB3r0ei8fjcu+998oll1wiU6dOPeV2ffX5HrRXdO6//361+OrTP9YvgaqqKlm0aJEsWbJEli5dOkB7fnZ687p74qabbpLPf/7zMm3aNFm8eLG8+uqr8sEHH8jGjRv77kX0QqJf92CV6Nd95513ysKFC2XatGnypS99SX75y1/KSy+9JIcOHerDV0H9raysTG655RaZMWOGXHbZZfLiiy9KXl6ePPnkkwO9a5QAkyZNkq9+9asya9YsmTt3rjz99NMyd+5cWbVq1UDvWq/cfffdsmvXLvnNb37TL883aK/o3HfffXLbbbeddptx48Z1//vEiRMyf/58mTt3rjz11FOnvV9hYaFEIhEJBoOmqzq1tbUD3lOrp6/7bI0bN05yc3Pl4MGDcvnll/fZ4/ZUIl/3345pbW2tFBUVdcdra2tlxowZvXrMvmL3dRcWFkJhajQalaamph6ds3PmzBGRT658jh8/vsf7m0i5ubni8Xhg9uPpPpeFhYU92n4w6s3rtkpKSpKZM2fKwYMHE7GLg8apjndmZqZjr+acyuzZs+Xtt98e6N3osXvuuad7MsWZrj721ed70A508vLyJC8vz9a2VVVVMn/+fJk1a5Y888wz4naf/kLVrFmzJCkpSdavXy833HCDiHxSzV5RUSFlZWVnve9noyevuy8cP35cGhsbTQOAgZDI111aWiqFhYWyfv367oFNKBSSLVu29HjGWl+z+7rLysokGAzK1q1bZdasWSIi8uabb0o8Hu8evNhRXl4uIjLgx1vj8/lk1qxZsn79elm8eLGIfHKJe/369XLPPfeo9ykrK5P169fLvffe2x1bt27dgH+Oe6I3r9sqFovJzp075eqrr07gng68srIymF481I53XykvLx+Un+NTMQxDvv71r8tLL70kGzdulNLS0jPep88+372plh5Mjh8/bkyYMMG4/PLLjePHjxvV1dXdP5/eZtKkScaWLVu6Y1/72teM0aNHG2+++abx4YcfGmVlZUZZWdlAvIReO3bsmLF9+3ZjxYoVRnp6urF9+3Zj+/btRktLS/c2kyZNMl588UXDMAyjpaXF+Nd//Vdj8+bNxpEjR4w33njDuOCCC4yJEycanZ2dA/Uyeqynr9swDOM//uM/jEAgYLz88svGjh07jOuuu84oLS01Ojo6BuIl9MqiRYuMmTNnGlu2bDHefvttY+LEicbNN9/c/f/W8/zgwYPGww8/bHz44YfGkSNHjJdfftkYN26ccemllw7USzij3/zmN4bf7zfWrFlj7Nmzx7jzzjuNQCBg1NTUGIZhGF/+8peN+++/v3v7d955x/B6vcYjjzxi7N2713jooYeMpKQkY+fOnQP1Enqlp697xYoVxuuvv24cOnTI2Lp1q3HTTTcZycnJxu7duwfqJfRKS0tL9+dXRIwf//jHxvbt241jx44ZhmEY999/v/HlL3+5e/vDhw8bqampxre+9S1j7969xhNPPGF4PB7jtddeG6iX0Cs9fd2rVq0y1q5daxw4cMDYuXOn8Y1vfMNwu93GG2+8MVAvocfuuusuIysry9i4caPp93R7e3v3Non6fA/5gc4zzzxjiIj68zdHjhwxRMTYsGFDd6yjo8P4l3/5F2PEiBFGamqq8YUvfME0OBoKbr31VvV1f/p1iojxzDPPGIZhGO3t7caVV15p5OXlGUlJScaYMWOMpUuXdn+ZDhU9fd2G8ckU83/7t38zCgoKDL/fb1x++eXG/v37+3/nz0JjY6Nx8803G+np6UZmZqZx++23mwZ31vO8oqLCuPTSS43s7GzD7/cbEyZMML71rW8Zzc3NA/QK7Hn88ceN0aNHGz6fz5g9e7bx3nvvdf/fZZddZtx6662m7X/7298a55xzjuHz+YzzzjvP+OMf/9jPe9w3evK677333u5tCwoKjKuvvtrYtm3bAOz12fnbtGnrz99e66233mpcdtllcJ8ZM2YYPp/PGDdunOlzPlT09HX/4Ac/MMaPH28kJycb2dnZxrx584w333xzYHa+l071e/rTxy9Rn2/X/90BIiIiIscZtLOuiIiIiM4WBzpERETkWBzoEBERkWNxoENERESOxYEOERERORYHOkRERORYHOgQERGRY3GgQ0RERI7FgQ4RERE5Fgc6RERE5Fgc6BAREZFj/f9a4rz9RYpONwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 100_000\n",
    "\n",
    "x, _ = datasets.make_moons(n_samples=n_samples, noise=.06)\n",
    "\n",
    "scaler = preprocessing.StandardScaler()\n",
    "x = scaler.fit_transform(x)\n",
    "\n",
    "plt.hist2d(x[:, 0], x[:, 1], bins=100)\n",
    "plt.xlim(-2 ,2)\n",
    "plt.ylim(-2, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "de9f55f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 1.\n",
    "\n",
    "def weight(t):\n",
    "    \"\"\" \\lambda(t)\n",
    "    \"\"\"\n",
    "    return 0.5 / np.log(sigma) * (sigma ** (2 * t) - 1)\n",
    "\n",
    "@partial(jax.jit, static_argnums=(3,))\n",
    "def loss_fn(params, x, t, score, key):\n",
    "\n",
    "    mu = x  # x(0)\n",
    "    std = np.sqrt(0.5 / np.log(sigma) * (sigma ** (2 * t) - 1))  # std of noise at time t\n",
    "\n",
    "    eps = jax.random.normal(key, shape=x.shape)  # Sampled noise\n",
    "    y = mu + std * eps  # x(t) = x(0) + std * eps  # Corrupted data\n",
    "\n",
    "    # Predicted score\n",
    "    pred = score.apply(params, np.concatenate([t, y], -1))  \n",
    "\n",
    "    # Score matching loss\n",
    "    loss = weight(t) * np.mean((pred + eps / std) ** 2)\n",
    "\n",
    "    return loss.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "e544f2d0-ae93-42d5-affa-9fb85bd5bc5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MLP(nn.Module):\n",
    "    \"\"\" A simple MLP in Flax. This is the score function.\n",
    "    \"\"\"\n",
    "    hidden_dim: int = 32\n",
    "    out_dim: int = 2\n",
    "    n_layers: int = 2\n",
    "\n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "        for _ in range(self.n_layers):\n",
    "            x = nn.Dense(features=self.hidden_dim)(x)\n",
    "            x = nn.gelu(x)\n",
    "        x = nn.Dense(features=self.out_dim)(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "337bbdfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def int_beta(t):\n",
    "    return t\n",
    "\n",
    "def weight(t):\n",
    "    return 1 - np.exp(-int_beta(t))\n",
    "\n",
    "@partial(jax.jit, static_argnums=(3,4,))\n",
    "def loss_fn(params, x, t, int_beta, score, key):\n",
    "    mu = x * np.exp(-0.5 * int_beta(t))\n",
    "    sigma = np.sqrt(1 - np.exp(-int_beta(t)))\n",
    "    eps = jax.random.normal(key, shape=x.shape)\n",
    "    y = mu + sigma * eps\n",
    "\n",
    "    pred = score.apply(params, np.concatenate([y, t], -1))\n",
    "    loss = weight(t) * np.mean((pred + eps / sigma) ** 2)\n",
    "    return loss.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "7455e305",
   "metadata": {},
   "outputs": [],
   "source": [
    "key = jax.random.PRNGKey(0)\n",
    "t = np.ones((x.shape[0], 1))\n",
    "\n",
    "score = MLP(hidden_dim=128, out_dim=2, n_layers=5)\n",
    "params = score.init(key, np.concatenate([x, t], axis=1))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1ec7fbad",
   "metadata": {},
   "source": [
    "## Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "c98ab327",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = optax.adamw(learning_rate=3e-4, weight_decay=1e-4)\n",
    "opt_state = opt.init(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "8c67c255",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2000/2000 [00:17<00:00, 112.91it/s, val=2.2799215] \n"
     ]
    }
   ],
   "source": [
    "n_steps = 2_000\n",
    "n_batch = 128\n",
    "T = 1.\n",
    "\n",
    "with trange(n_steps) as steps:\n",
    "    for step in steps:\n",
    "\n",
    "        # Draw a random batches from x\n",
    "        key, subkey = jax.random.split(key)\n",
    "        idx = jax.random.choice(key, x.shape[0], shape=(n_batch,))\n",
    "        \n",
    "        x_batch = x[idx]\n",
    "        t_batch = jax.random.uniform(key, shape=(x_batch.shape[0], 1), minval=0., maxval=T)\n",
    "\n",
    "        loss, grads = jax.value_and_grad(loss_fn)(params, x_batch, t_batch, int_beta, score, key)\n",
    "        updates, opt_state = opt.update(grads, opt_state, params)\n",
    "\n",
    "        params = optax.apply_updates(params, updates)\n",
    "\n",
    "        steps.set_postfix(val=loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a018c3d0",
   "metadata": {},
   "source": [
    "## Distillation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "90b560b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "key = jax.random.PRNGKey(0)\n",
    "t = np.ones((x.shape[0], 1))\n",
    "\n",
    "score_student = MLP(hidden_dim=128, out_dim=2, n_layers=5)\n",
    "params_student = score_student.init(key, np.concatenate([x, t], axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "f8ce7cb4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5000/5000 [00:50<00:00, 98.17it/s, val=4.15465]    \n"
     ]
    }
   ],
   "source": [
    "n_steps = 5_000\n",
    "n_batch = 128\n",
    "\n",
    "def alpha(t):\n",
    "    return np.exp(-0.5 * int_beta(t))\n",
    "\n",
    "def sigma(t):\n",
    "    return np.sqrt(1 - np.exp(-int_beta(t)))\n",
    "\n",
    "@partial(jax.jit, static_argnums=(4,5,))\n",
    "def loss_distillation_fn(params_student, params, x, t, score_student, score, key):\n",
    "    eps = jax.random.normal(key, shape=x.shape)\n",
    "\n",
    "    t = t[:, None] \n",
    "\n",
    "    z_t = alpha(t) * x + sigma(t) * eps\n",
    "\n",
    "    t_p = t - 0.5 / N\n",
    "    t_pp = t - 1 / N\n",
    "\n",
    "    x_hat_t = score.apply(params, np.concatenate([z_t, t], -1))\n",
    "    z_tp = alpha(t_p) * x_hat_t + sigma(t_p) / sigma(t) * (z_t - x_hat_t * alpha(t))\n",
    "\n",
    "    x_hat_tp = score.apply(params, np.concatenate([z_tp, t], -1))\n",
    "    z_tpp = alpha(t_pp) * x_hat_tp + sigma(t_pp) / sigma(t_p) * (z_t - x_hat_tp * alpha(t_p))\n",
    "\n",
    "    x_tilde = (z_tpp - (sigma(t_pp) / sigma(t)) * z_t) / (alpha(t_pp) - (sigma(t_pp) / sigma(t)) * alpha(t))\n",
    "\n",
    "    pred = score_student.apply(params_student, np.concatenate([z_t, t], -1))\n",
    "    loss = weight(t) * np.mean((pred + (x_tilde - x * alpha(t)) / sigma(t) ** 2) ** 2)\n",
    "    return loss.mean()\n",
    "\n",
    "N = 100\n",
    "\n",
    "with trange(n_steps) as steps:\n",
    "    for step in steps:\n",
    "\n",
    "        # Draw a random batches from x\n",
    "        key, subkey = jax.random.split(key)\n",
    "        idx = jax.random.choice(key, x.shape[0], shape=(n_batch,))\n",
    "        \n",
    "        x_batch = x[idx]\n",
    "        t_batch = jax.random.choice(key, np.arange(1, N + 1), shape=(n_batch,)) / N\n",
    "\n",
    "        loss, grads = jax.value_and_grad(loss_distillation_fn)(params_student, params, x_batch, t_batch, score_student, score, key)\n",
    "        updates, opt_state = opt.update(grads, opt_state, params_student)\n",
    "\n",
    "        params_student = optax.apply_updates(params_student, updates)\n",
    "\n",
    "        steps.set_postfix(val=loss)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ca323255",
   "metadata": {},
   "source": [
    "## Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "148b0f99",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}