exp002-sir-fourier-latent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting BC covid-19 data using an SIR model with varying $\\beta(t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import torch as tc\n",
    "import numpy as np\n",
    "from torchdiffeq import odeint\n",
    "from torch import nn\n",
    "from data import get_cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tc.cuda.is_available()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We retrieve the BC data. We have that $S(t) = S_0 - C(t)$, where $C(t)$ is the number of cumulative cases at time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cases, cum_cases, N = get_cases()\n",
    "\n",
    "I0 = 1\n",
    "R0 = 0\n",
    "S0 = N - I0 - R0\n",
    "y0 = tc.tensor([S0, I0, R0, N]).float()\n",
    "\n",
    "data = S0 - cum_cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcd0077cd30>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGsCAYAAAD+L/ysAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJG0lEQVR4nO3de1xUdf4/8NcZYGYQmJGLcheUq6B4wwt4K29parZtpWVoZRluF9t+1qbVlvvNxXIv6W5Rpol2kbZF0raVtE3wfgdDTUVRQeQiKgwgzMDM+f2BTo5cRy5nmHk9H4/zMM75nOF99mzOq8/5fD5HEEVRBBEREZEFk0ldABEREVFLGFiIiIjI4jGwEBERkcVjYCEiIiKLx8BCREREFo+BhYiIiCweAwsRERFZPAYWIiIisngMLERERGTxGFiIiIjI4lldYNm5cyemT58OHx8fCIKAb7/91qzz33nnHQiC0GBzcnLqmIKJiIioRVYXWKqqqjBgwAD885//vKvzFy1ahMLCQpMtIiICjzzySDtXSkRERK1ldYFlypQpePfdd/HQQw81elyn0+G1116Dr68vnJycMHz4cKSnpxuPOzs7w8vLy7gVFxfj5MmTmDdvXiddAREREd3JXuoCOttTTz2FCxcuIDk5GT4+PkhNTcXkyZORnZ2NkJCQBu3XrFmD0NBQjB49WoJqiYiICLDCHpbmnDt3Dhs3bsQ333yD0aNHIygoCIsWLcKoUaOwbt26Bu21Wi2+/PJL9q4QERFJzKZ6WI4ePQpRFBEaGmqyX6vVwt3dvUH7TZs2oaKiAnPmzOmsEomIiKgRNhVYDAYD7OzscOTIEdjZ2Zkcc3Z2btB+zZo1mDZtGry8vDqrRCIiImqETQWWQYMGQa/Xo6SkpMUxKefPn8eOHTuwZcuWTqqOiIiImmJ1gaWyshJnz541/nz+/HlkZWXBzc0NoaGhmD17NubMmYO//vWvGDRoEEpLS/HTTz+hf//+uP/++43nffbZZ/D29saUKVOkuAwiIiK6jSCKoih1Ee0pPT0d9957b4P9c+fORVJSEmpra/Huu+9iw4YNKCgogLu7O2JiYrB06VL0798fQP2jo4CAAMyZMwfLli3r7EsgIiKiO1hdYCEiIiLrY1PTmomIiKhrYmAhIiIii2c1g24NBgMuX74MFxcXCIIgdTlERETUCqIooqKiAj4+PpDJmu5HsZrAcvnyZfj7+0tdBhEREd2F/Px8+Pn5NXncagKLi4sLgPoLVqlUEldDREREraHRaODv72/8Hm+K1QSWW4+BVCoVAwsREVEX09JwDg66JSIiIovHwEJEREQWj4GFiIiILB4DCxEREVk8BhYiIiKyeAwsREREZPEYWIiIiMjiMbAQERGRxWNgISIiIovHwEJEREQWj4GFiIiILB4DCxEREVk8q3n5YUdZu/s88q/dgCAAAgTcejeTANTvEwQYX9d0W5tb+4RG9uHmObeO/dru188EALmdDK5Ocrh2c4CnSgl/t25QOzp0xmUTERFZFAaWFvzn58vIzCuTugwjtaMD+vuqMTjAFYN7dcegXq4MMUREZPUYWFrw8BA/xAa5QxQBEYAo1u8XUb+jfp94237cbPvrPtxq08ixXz+zft+t4wCgrTPg+o1aXK/SobC8GqWVOpRX12L32VLsPlsKoL5HJqSnM4YEuGJwL1cMCXBFbw+nFl/TTURE1JUIonj712rXpdFooFarUV5eDpVKJXU5HaJKW4fzpVXIyi/D0bzrOHLxOi5evdGgnWs3Bwzv7Y6xYT0wNrQHfLo7SlAtERFRy1r7/c3A0sWVVmpx5OJ1HL14HUfzruPYpXLo6gwmbUI9nTE2tAdigzwwtLcbnBXsWCMiIsvAwGKjdHUGHL9cjt05pUg/XYKs/DIYbrvDdjIBUX5qjA6p730Z6N8ddjI+PiIiImkwsBAAoOyGDjtzSrEnpxT7cq8i75rpIyS1owOGBrphUK/uGOjfHRHeKrg6ySWqloiIbA0DCzXq0vUb2HO2FDvPlGJXzhVoauoatPFUKRDupUK4twv63vyzj4cz5PZctoeIiNoXAwu1qE5vwM8F5Th68Toy88vw86Uy5F+rbrStg52AoB7OCPdyQYinC3y6K+GlcoS3WgkvtRJKB7tOrp6IiKwBAwvdlUptHU4XVeBUkQanCm/+WVSBikZ6Ym7n2s0BPV2UUDs6QOXoAJWjPVRKh19/Vtrf/PPXYypHB7go7CHjGBoiIpvV2u9vThchE84KewwJqF/P5RZRFHG5vAanCuvDy7krlSjW1KCwvAaFZTWortXXrxdzo9bs3ycI9b9TpXSASxOh5tewc/u++jYuSgcOGiYisgEMLNQiQRDg290Rvt0dMb6vp8kxURShqalDUXkNSipqoKmug6amFuXVtdBU3/yzpg7l1bWoqKnfp6mpg6a6Fto6A0QRqKipa7EHpzn1gac+1LgoGws69eHG3UkBDxcF3J3k6OGi4GMsIqIuhIGF2kQQBKgd6x/9hHm5mHVuTa0eFTX1AUdTXXvbP/+679bPFTW/Bp1b+6pr9QDqH2NVautwubzGrN/vJLczBhgPZwXcnRXo4SyHu7MCHs4K9HBRILinM9w4a4qISHIMLCQZpYMdlA526OGiuKvzdXUGY5CpaCLoaGrqg9D1Gzpcq9KhtEKL0koddHoDqnR6VF290ehqwbfr6aJAhI8KkT4qRPqo0d9XDT9XR77+gIioEzGwUJclt5fB/WbPiDlEUUSFtg5XK3UordTWh5ibYeZqlRalFTpcrdKisLwGl65Xo6RCi5LTV5B++orxM3y7O2JksDtGhfTA2JAeUHfjCyiJiDoSZwkRNePWrKmThRqcvFyO4wUanCrSoFb/6782DnYCRgV7YGqUD+6L9ISLkuGFiKi1OK2ZqIPc0NXh0IXr2J1zBRlnruBMcaXxmNJBhvsivfDgIF/E9HHnwF4iohYwsBB1krMllfhvdiE2ZxXg3JUq435nhT1mD++FJ0YEwN+tm4QVEhFZLgYWok4miiKOXSpHypFL2H6yGEWa+llLggCMCvbA48N6YUKEJxzs+IoDIqJbWvv9bdbfnImJiYiKioJKpYJKpUJMTAy2bt3aZPsnn3wSgiA02CIjI41tkpKSGm1TU2PeFFUiqQmCgIH+3fF/D/bDvsXjsGZONEYFe0AUgV05pVjw5VHEJPyEv207jUvXm5+ZREREpsyaJeTn54fly5cjODgYALB+/XrMmDEDmZmZJiHklpUrV2L58uXGn+vq6jBgwAA88sgjJu1UKhVOnz5tsk+pVJpTGpFFEQQBEyI8MSHCExevVuHrQ/n41+FLKK3UYtVPZ/GPHWfx6BB/vDAumI+LiIhaoc2PhNzc3LBixQrMmzevxbbffvstHnroIZw/fx4BAQEA6ntYXn75ZZSVlbWlDD4SIotXqzdg24lifLH/IvblXgUAyARgSn9vxI8JQn8/tcQVEhF1vg55JHQ7vV6P5ORkVFVVISYmplXnrF27FhMmTDCGlVsqKysREBAAPz8/TJs2DZmZmS1+llarhUajMdmILJmDnQxTo7yxcf4IpCyIwegQDxhE4PufCzH9n7vx6jfHcLVSK3WZREQWyezAkp2dDWdnZygUCsTHxyM1NRUREREtnldYWIitW7fimWeeMdkfHh6OpKQkbNmyBRs3boRSqcTIkSORk5PT7OclJCRArVYbN39/f3MvhUgyQwLc8Pm84di6cDRmDPQBAHxz5BLG/TUD/zqUD4PBKsbCExG1G7MfCel0OuTl5aGsrAwpKSlYs2YNMjIyWgwtCQkJ+Otf/4rLly9DLm/63SwGgwGDBw/GmDFjsGrVqibbabVaaLW//teoRqOBv78/HwlRl3Tk4nW89e1xnCys7ykc3tsNa+ZGcxE6IrJ6nTatecKECQgKCsInn3zSZBtRFBEaGopp06bh73//e4uf+eyzz+LSpUvNzkC6E8ewUFdXpzfgsz3n8cGPObih06OXWzf89dEBGBroJnVpREQdpsPHsNwiiqJJT0djMjIycPbs2VYNzBVFEVlZWfD29m5raURdir2dDPPHBCF5/gj4dndE3rUbePSTfVi+9RS0dXqpyyMikpRZgWXJkiXYtWsXLly4gOzsbLzxxhtIT0/H7NmzAQCLFy/GnDlzGpy3du1aDB8+HP369WtwbOnSpfjhhx+Qm5uLrKwszJs3D1lZWYiPj7/LSyLq2qL8umPry6Px28F+EEXg44xzmL/hCOr0BqlLIyKSjFmBpbi4GHFxcQgLC8P48eNx4MABpKWlYeLEiQDqB9bm5eWZnFNeXo6UlJQme1fKysowf/589O3bF5MmTUJBQQF27tyJYcOG3eUlEXV9KqUD/vroAHz8xBAoHWTIOHMFcz47iFLOIiIiG8Wl+Yks3PaTxViYnIkbOj28VEp8OHswhgS4Sl0WEVG76LQxLETUsSZGeGLz8yMR1MMJRZoazPxkH1KOXJK6LCKiTsXAQtQFhHi6YPMLozA1yht1BhGL/n0MyQfzWj6RiMhKMLAQdRHOCnv887FBmBsTAFEEXt+UjQ93nIWVPNUlImoWAwtRFyIIAt55IBK/uycIALDih9NY+t1JroxLRFaPgYWoixEEAa9NDsdb0+pXl07aewF/2Xa6hbOIiLo2BhaiLmreqN54/7dRAICP0s/ho/SzEldERNRxGFiIurBHh/rj/00MBQC8n3Yan2Sck7giIqKOwcBC1MW9OD4EiybVh5aErafwxf6LEldERNT+GFiIrMAL40Lw/L31A3Hf2nwcqZlcp4WIrAsDC5GVWDQpDE/GBkIUgUXf/Iy040VSl0RE1G4YWIishCAI+OO0CDw8xA96g4iXNmZi55krUpdFRNQuGFiIrIhMJmD5Q/1xf38v6PQGzP/8MA5duCZ1WUREbcbAQmRl7O1k+GDmIIwN7YGaWgPmJR1CQVm11GUREbUJAwuRFZLby/DxE0MwwE8NTU0dnlp3EGU3dFKXRUR01xhYiKyUo9wOH84eDE+VAmeKK/HkukOo0tZJXRYR0V1hYCGyYn6u3fD5vOHo3s0BWflleO7zI9DW6aUui4jIbAwsRFYu1NMF654cim5yO+w+W4oXv8pEnd4gdVlERGZhYCGyAYN6ueLTOdGQ28uw7WQxXv33z3zDMxF1KQwsRDZiZLAHPnp8MOxkAlIzCzBv/SGkny6RuiwiolZhYCGyIRMiPPG3RwdAEIAdp6/gyXWHsO/cVanLIiJqEQMLkY2ZMdAXmxbEYmigKwDgtZRjnD1ERBaPgYXIBg3q5YrPnhwK3+6OyL9WjffTTkldEhFRsxhYiGyUi9IBy3/bHwCwft9FJGz9BTd07GkhIsvEwEJkw0aH9MCTsYEAgE8ycjEv6TBEkbOHiMjyMLAQ2bg3p/bFn2ZEwl4mYF/uVezP5csSicjyMLAQ2Th7OxnmxARi1jB/AMCHO86yl4WILA4DCxEBAOaPDoK9TMDus6VI2ntB6nKIiEwwsBARAKCXezcsub8vAGDZ97/g8AU+GiIiy8HAQkRGT40MxPQBPqgziFiYnIWKmlqpSyIiAsDAQkS3EQQBf/5NP/i5OqKgrBqvp2TznUNEZBEYWIjIhIvSAX+fORAOdgK+zy7Eci4qR0QWgIGFiBoYGuiGFQ8PAACs3pmLdXvOS1wREdk6BhYiatSDg3zx2uQwAMCf/nMSW7MLJa6IiGyZWYElMTERUVFRUKlUUKlUiImJwdatW5ts/+STT0IQhAZbZGSkSbuUlBRERERAoVAgIiICqampd3c1RNSuFowNwhMjekEUgYVfZ3HmEBFJxqzA4ufnh+XLl+Pw4cM4fPgwxo0bhxkzZuDEiRONtl+5ciUKCwuNW35+Ptzc3PDII48Y2+zbtw8zZ85EXFwcjh07hri4ODz66KM4cOBA266MiNpMEAQsfaAfJvT1hK7OgHnrDyOnuELqsojIBgliG5e0dHNzw4oVKzBv3rwW23777bd46KGHcP78eQQEBAAAZs6cCY1GY9JTM3nyZLi6umLjxo2trkOj0UCtVqO8vBwqlcr8CyGiJlXr9Hjs0/3Iyi+DazcHrH96GKL8uktdFhFZgdZ+f9/1GBa9Xo/k5GRUVVUhJiamVeesXbsWEyZMMIYVoL6HZdKkSSbt7rvvPuzdu7fZz9JqtdBoNCYbEXUMR7kdPntyKPr7qnH9Ri0eW70fGWeuSF0WEdkQswNLdnY2nJ2doVAoEB8fj9TUVERERLR4XmFhIbZu3YpnnnnGZH9RURE8PT1N9nl6eqKoqKjZz0tISIBarTZu/v7+5l4KEZnBzUmOjfNHYGSwO6p0ejyddAgf7jgLPddpIaJOYHZgCQsLQ1ZWFvbv348FCxZg7ty5OHnyZIvnJSUloXv37njwwQcbHBMEweRnURQb7LvT4sWLUV5ebtzy8/PNug4iMp+zwh7rnhyGhwb7Qm8QseKH03j80/24XFYtdWlEZOXMDixyuRzBwcGIjo5GQkICBgwYgJUrVzZ7jiiK+OyzzxAXFwe5XG5yzMvLq0FvSklJSYNelzspFArjbKVbGxF1PLm9DH99ZAD+8sgAOMntcOD8NUz+YCe+/5nTnomo47R5HRZRFKHVapttk5GRgbNnzzY6MDcmJgbbt2832bdt2zbExsa2tTQi6iCCIODhIX74/qXRGOCnhqamDs9/dRSvfnMMVdo6qcsjIitkb07jJUuWYMqUKfD390dFRQWSk5ORnp6OtLQ0APWPaQoKCrBhwwaT89auXYvhw4ejX79+DT5z4cKFGDNmDN577z3MmDEDmzdvxo8//ojdu3e34bKIqDMEejjh3wti8cGPZ/BR+jl8c+QSdp8txZtTI3B/f68WH+0SEbWWWT0sxcXFiIuLQ1hYGMaPH48DBw4gLS0NEydOBFA/sDYvL8/knPLycqSkpDQ57Tk2NhbJyclYt24doqKikJSUhK+//hrDhw+/y0sios7kYCfDq/eFY+OzI+Dn6ojC8ho8/9VRzFy9H3vPlaKNKycQEQFoh3VYLAXXYSGSXk2tHonp55CYcQ66OgMAYFigGxZOCEFskDt7XIiogdZ+fzOwEFG7u1xWjcT0c/j6UD50+vrgEh3gioUTQjAq2IPBhYiMGFiISHKF5dX4OP0cNh7KN/a4DO7VHQsnhGJMCIMLETGwSF0OEd2mWFODxPRz+OpgnjG4hHu54LXJYbg3rCeDC5ENY2AhIotToqnBxxm52HgwD9W1egBAlJ8aL9wbjIkRngwuRDaIgYWILFb5jVp8lH4WG/ZdNAaXcC8XvDguBFP6eUEmY3AhshUMLERk8a5WarFm93ls2HsBVbr64BLc0xnP3xuE6VE+sLdr89qWRGThGFiIqMsou6HDuj0XsG7PeWhq6lfKDXDvhufvCcZDg30ZXIisGAMLEXU5mppafL7vItbsysX1G7UAgFBPZ7wxNQJjQ3tIXB0RdQQGFiLqsqq0dfjywEV8lH4OZTeDy9jQHnhjal+EerpIXB0RtScGFiLq8spv1OIfP+Vg/b4LqNWLkAnA48N74dVJ4VB3c5C6PCJqBwwsRGQ1LpRWYfnWU0g7UQQA8HBW4I/TIzA9yptToYm6uNZ+f3MkGxFZvEAPJ3wcNwTJ80cguKczSiu1eGljJuZ/fgSamlqpyyOiTsDAQkRdxog+7vj+pVH4/YRQyO1k2H6yGA9+uAdniiukLo2IOhgDCxF1KQp7OyycEIKUBbHwViuRe6UKD/xzNzLOXJG6NCLqQAwsRNQl9fdTY8sLozAq2AM1tQY8s/4QNmcVSF0WEXUQBhYi6rJ6uCjw2ZNDMS3KG7V6EQuTs7B65zlYyVwCIroNAwsRdWlyexlWzRqEp0f2BgD8+b+nsOKH0xJXRUTtjYGFiLo8mUzAH6dH4M2pfQEAH6Wfw4HcqxJXRUTtiYGFiKzGM6P7YNZQfwDACxszUaypkbgiImovDCxEZFXemhaBUE9nXKnQYsEXR6CrM0hdEhG1AwYWIrIqTgp7rI6Lhkppj6N5ZfjLNo5nIbIGDCxEZHUCPZyw4pEBAIDVO3ORdrxQ4oqIqK0YWIjIKt0X6YV5o+pnDr21+QQqtXUSV0REbcHAQkRW6w+TwxHg3g1XKrT44+bjUpdDRG3AwEJEVktuL8OKhwdAJgCbjhZgy7HLUpdERHeJgYWIrNqw3m54cVwIAOCdLSdwtVIrcUVEdDcYWIjI6j1/bzDCPF1wrUqHl7/OQp2eU52JuhoGFiKyenJ7GT6YNRCODnbYlVOKxPRzUpdERGZiYCEim9DXW4Vlv+kHAPjnjrM4mndd4oqIyBwMLERkM34zyBfjw3tCW2fAvKRDKOHS/URdBgMLEdkMQRCw6rFBiPBW4fqNWiz97qTUJRFRKzGwEJFNcVLYY8UjUbCTCfg+uxApRy5JXRIRtQIDCxHZnEgfNV4cFwwAWJyajaz8MmkLIqIWMbAQkU16aVwIJvT1hK7OgPkbDnM8C5GFMyuwJCYmIioqCiqVCiqVCjExMdi6dWuz52i1WrzxxhsICAiAQqFAUFAQPvvsM+PxpKQkCILQYKup4V8eRNRxZDIBf585ACE9nVFSocXzXx1FLddnIbJY9uY09vPzw/LlyxEcXN+Vun79esyYMQOZmZmIjIxs9JxHH30UxcXFWLt2LYKDg1FSUoK6OtOXkKlUKpw+bfoKeKVSaU5pRERmc1E6YPWcaDzwj904dOE63tt6Cm9Oi5C6LCJqhFmBZfr06SY/L1u2DImJidi/f3+jgSUtLQ0ZGRnIzc2Fm5sbACAwMLBBO0EQ4OXlZU4pRETtoreHE/7y6AA89/kRrNl9HoN6uWJqlLfUZRHRHe56DIter0dycjKqqqoQExPTaJstW7YgOjoa77//Pnx9fREaGopFixahurrapF1lZSUCAgLg5+eHadOmITMzs8Xfr9VqodFoTDYiortxX6QXnhvbBwDw2r+P4WxJhcQVEdGdzA4s2dnZcHZ2hkKhQHx8PFJTUxER0XgXam5uLnbv3o3jx48jNTUVH3zwAf7973/j+eefN7YJDw9HUlIStmzZgo0bN0KpVGLkyJHIyclpto6EhASo1Wrj5u/vb+6lEBEZvTopDDF93FGl0yP+i6Oo0ta1fBIRdRpBFEXRnBN0Oh3y8vJQVlaGlJQUrFmzBhkZGY2GlkmTJmHXrl0oKiqCWq0GAGzatAkPP/wwqqqq4Ojo2OAcg8GAwYMHY8yYMVi1alWTdWi1Wmi1v751VaPRwN/fH+Xl5VCpVOZcEhERAOBKhRbT/rELxRotpkV54x+PDYIgCFKXRWTVNBoN1Gp1i9/fZvewyOVyBAcHIzo6GgkJCRgwYABWrlzZaFtvb2/4+voawwoA9O3bF6Io4tKlxhdrkslkGDp0aIs9LAqFwjhb6dZGRNQWPVwU+Gj2YNjLBPzn50J8tueC1CUR0U1tXodFFEWTno7bjRw5EpcvX0ZlZaVx35kzZyCTyeDn59fk52VlZcHbm4PeiKjzDQlwwxtT+wIA3v3+JDZnFUhcEREBZgaWJUuWYNeuXbhw4QKys7PxxhtvID09HbNnzwYALF68GHPmzDG2f/zxx+Hu7o6nnnoKJ0+exM6dO/Hqq6/i6aefNj4OWrp0KX744Qfk5uYiKysL8+bNQ1ZWFuLj49vxMomIWu/J2EDEjQiAKAKv/OsYfjhRJHVJRDbPrMBSXFyMuLg4hIWFYfz48Thw4ADS0tIwceJEAEBhYSHy8vKM7Z2dnbF9+3aUlZUhOjoas2fPxvTp003GppSVlWH+/Pno27cvJk2ahIKCAuzcuRPDhg1rp0skIjKPIAhY+kAkfjvYD3qDiBe/ysTOM1ekLovIppk96NZStXbQDhFRa9XpDXgpORP/zS6C0kGGtXOHYmSwh9RlEVmVDht0S0RkK+ztZPhg5iDcG9YDNbUGzPnsINbsyoWV/HceUZfCwEJE1Ay5vQyJTwzBbwb5Qm8Q8e73v+DFjZmoqKmVujQim8LAQkTUAqWDHf726AAsfSDSOOX5/lW7sD/3qtSlEdkMBhYiolYQBAFzYwPx9XMj4OfqiPxr1Zi1ej/i1h7A4QvXpC6PyOpx0C0RkZkqamqxfOspfH0oH3WG+r9CB/ipMTc2EFP6ecNRbidxhURdR2u/vxlYiIjuUv61G/hwx1lsOloAnd4AAOgmt8O94T0xOdILI/q4o4eLQuIqiSwbAwsRUScprdRi44E8/OtIPvKvmb6N3ketRJRfd0T5qzHArzv6+6mhUjpIVCmR5WFgISLqZKIo4udL5dh6vAg/nSpGTkklGvsb1kulRG8PJwS4d4O7sxyu3eRwd5bDzUkBt25yuDnL4e4kh9KBj5bI+jGwEBFJrFJbh+xL5fj5Uhl+vlSOY5fKcOl6dcsn3uToYAc3p1thRl4fZpzqA82tf3Z3lsPdSQFPlZJjZ6hLYmAhIrJAZTd0yC2twoXSKly6Xo1rVTpcrdLhWpUW16pqb/6pQ63e/L+a3ZzkCHDvhgC3bghwd0KfHk6I9FGjt4cT7GRCB1wNUdsxsBARdVGiKKJSW2cMM9eNocZ0u3XsSoUW1bX6Jj9P6SDD1P4+uDe8B8aE9uAYGrIoDCxERDZCFEVUaOtw6Vo1Ll6twoWrN3DxahVOF1fg5GUNtHUGY1t7mYB+vmq8OC4YY0N7wN6Oy3GRtBhYiIgIeoOIQxeu4ceTxdhxugTnrlQZj/l2d8RDg30xNrQHogPdJKySbBkDCxERNXChtAob9l3Et1kFuFalAwAIAhA3IgC/uycYXmqlxBWSrWFgISKiJtXU6rE5qwD/+bkQu3JKAdQveje1vzdmDeuFIQGuEldItoKBhYiIWiSKInafLcXft5/B0bwyAPVvqH57egRmDe3F2UXU4RhYiIio1QwGERk5V/DPn87iyMXrAIDoAFfMjQ1ETJA7PJz5igHqGAwsRERktlq9ARv2XcTft59BpbYOAOCisMe/4mPQ15t/t1L7a+33N+ezERGRkYOdDPNG9Ubay6Px1MhAqB0dUKGtQ9zag/jfL8VSl0c2jD0sRETUpGJNDR5bvR+5pfXToYN6OGHlrEHo56uWuDKyFuxhISKiNvNUKfHfhaPx9MjeAIBzV6rw3OdH8EuhRuLKyNYwsBARUbOUDnb44/QI7Fs8Dj5qJQrKqjHjwz1YsyvXOM6FqKMxsBARUat4qx3x3YujMD68J3R1Brz7/S8YufwnHC8ol7o0sgEMLERE1GruzgqsmRuNt6ZFwM1JjvLqWsz57CDSjhdJXRpZOQ66JSKiu3KtSofHP92PU0UVAIDRIR6YOdQf06J8JK6MuhIOuiUiog7l5iTH5hdGYsE9QRAEYFdOKV74KhMZZ65IXRpZIQYWIiK6awp7O/xhcjj+8+IoDPCrn+r8dNIhLPv+JHR1BomrI2vCwEJERG0W6aPGF88Mx/39vaA3iPh013nEJPwPq3eek7o0shIMLERE1C5clA74aPYQrI4bAhelPa5W6fDn/57CG6nZuF6lk7o86uI46JaIiNpdtU6PP//3F3y+/yIAoIeLAo8M8cO94T0xNNBN4urIkvDlh0REJLk9Z0vx9pYTOFtSadz3+wmhmD+mDxzldhJWRpaCgYWIiCxCTa0e3x27jM1Zl7H7bCkAoHs3BzwwwAcL7gmCt9pR4gpJSgwsRERkUQwGEZ/vv4g1u3ORf60aQH1wmRMTiJlD/eGjVkIQBImrpM7WIeuwJCYmIioqCiqVCiqVCjExMdi6dWuz52i1WrzxxhsICAiAQqFAUFAQPvvsM5M2KSkpiIiIgEKhQEREBFJTU80pi4iIugCZTMDc2ECkL7oXn8QNQR8PJ5TdqMWq/+Vg5PKfELv8J/zrcD6qdXqpSyULZFYPy3fffQc7OzsEBwcDANavX48VK1YgMzMTkZGRjZ4zY8YMFBcX491330VwcDBKSkpQV1eH2NhYAMC+ffswevRo/N///R9+85vfIDU1FX/84x+xe/duDB8+vNUXwh4WIqKupVZvwLYTxVi/9wIOXrhm3K9S2uO+SC/EBLnjwYG+kMnY62LNOu2RkJubG1asWIF58+Y1OJaWloZZs2YhNzcXbm6NjwqfOXMmNBqNSU/N5MmT4erqio0bN7a6DgYWIqKuq6KmFuv2XMC/Dufj0vVq4/7eHk54eIgfJvT1RHBPZ9gxvFidDl+aX6/XIzk5GVVVVYiJiWm0zZYtWxAdHY33338fvr6+CA0NxaJFi1Bd/ev/Gfft24dJkyaZnHffffdh7969zf5+rVYLjUZjshERUdfkonTAS+NDsPPVe5H01FA8M6o3nOR2OF9ahRU/nMZ9H+zElJU7kXHmCgwGqxh6SWayN/eE7OxsxMTEoKamBs7OzkhNTUVERESjbXNzc7F7924olUqkpqaitLQUv/vd73Dt2jXjOJaioiJ4enqanOfp6Ymioubf/JmQkIClS5eaWz4REVkwmUzAPWE9cU9YT7w8MRRbswux6WgBDl+8hjPFlZj72UH4uTpi1lB/PBLtD0+VUuqSqZOY/UhIp9MhLy8PZWVlSElJwZo1a5CRkdFoaJk0aRJ27dqFoqIiqNX175jYtGkTHn74YVRVVcHR0RFyuRzr16/HY489Zjzvyy+/xLx581BTU9NkHVqtFlqt1vizRqOBv78/HwkREVmhEk0NPtxxFqmZBdDU1AEA7GQCxoX3xMxof9wT1gP2dly8vStq7SMhs3tY5HK5cdBtdHQ0Dh06hJUrV+KTTz5p0Nbb2xu+vr7GsAIAffv2hSiKuHTpEkJCQuDl5dWgN6WkpKRBr8udFAoFFAqFueUTEVEX1FOlxNIZ/bD4/r74b3YhNh7Mw6EL17H9ZDG2nyyGi8Ie8fcEYW5sIJwVZn+1URfQ5jgqiqJJT8ftRo4cicuXL6Oy8tcVDs+cOQOZTAY/Pz8AQExMDLZv325y3rZt24yziIiIiG5ROtjhocF++CY+Ftt/PwbPju4NNyc5KrR1WPHDacQm/A+f7szFpes3pC6V2plZj4SWLFmCKVOmwN/fHxUVFUhOTsby5cuRlpaGiRMnYvHixSgoKMCGDRsAAJWVlejbty9GjBiBpUuXorS0FM888wzGjh2LTz/9FACwd+9ejBkzBsuWLcOMGTOwefNmvPnmm5zWTEREraI3iPj6UD7W7MpFbmkVAEAQgKdieyN+bB/05DgXi9Yhs4SKi4sRFxeHsLAwjB8/HgcOHDCGFQAoLCxEXl6esb2zszO2b9+OsrIyREdHY/bs2Zg+fTpWrVplbBMbG4vk5GSsW7cOUVFRSEpKwtdff21WWCEiIttlJxPw+PBe2P7KWPz5N/0xwL87RBH4bM95jHp/B978Nhs5xRVSl0ltxKX5iYjI6uw4VYJ//JSDo3llAOpDzeIp4ZgTEwi5PQfnWhK+S4iIiGyaKIrYn3sNK/93Bvtz61fSDenpjBWPDMBA/+7SFkdGDCxERESof+nixkN5+Nu2M7hapYNMAGYO7YVlD/bjsv8WoMNXuiUiIuoKZDIBs4cHYPsrY/HgQB8YRGDjwTws+vcx6LlqbpfBwEJERDbBzUmOD2YNwoqHoyATgE1HC/D2luOwkgcNVo+BhYiIbMoj0f748PHBEATgi/152Jx1WeqSqBUYWIiIyOZM6e+NF8eFAABe/fcx/O+XYokropYwsBARkU16aVwwpg/wQa1exDvfneCjIQvHwEJERDbJ3k6G937bH0oHGfKvVeN4gUbqkqgZDCxERGSzusntMaFv/ct2V+/Klbgaag4DCxER2bQF9wQBALZmF6KkokbiaqgpDCxERGTTIn3UGNyrO+oMIr47Vih1OdQEBhYiIrJ5U/p5AwAyzlyRuBJqCgMLERHZvLFhPQAAB3KvoqZWL3E11BgGFiIisnkhPZ3hrVZCW2fAgfPXpC6HGsHAQkRENk8QBIwJqe9lyTjNx0KWiIGFiIgIvz4W2pnDwGKJGFiIiIgAjAz2gL1MwNmSSpy4XC51OXQHBhYiIiIAakcHTO7nBQD4+lC+xNXQnRhYiIiIbpoW5QMA2HfuqsSV0J0YWIiIiG4a3tsNAJBTUomrlVqJq6HbMbAQERHd5OokR7iXCwDgIKc3WxQGFiIiotvc6mXheiyWhYGFiIjoNgP8uwMAzhRXSFsImWBgISIiuo1vd0cAwOWyaokrodsxsBAREd3G51ZgKa+BKIoSV0O3MLAQERHdxlOlhCAAujoDrlbppC6HbmJgISIiuo3cXoYezgoAfCxkSRhYiIiI7uDv1g0AkHfthsSV0C0MLERERHcIcK8PLBevMrBYCgYWIiKiOwS6OwEALpRWSVwJ3cLAQkREdAdjDwsfCVkMBhYiIqI7BNzsYbl4lT0sloKBhYiI6A4BNwfdFmu0qNbpJa6GAAYWIiKiBrp3c4BKaQ+AM4UshVmBJTExEVFRUVCpVFCpVIiJicHWrVubbJ+eng5BEBpsp06dMrZJSkpqtE1NTc3dXxUREVEbCIKAQI/6x0I5JXynkCWwN6exn58fli9fjuDgYADA+vXrMWPGDGRmZiIyMrLJ806fPg2VSmX8uUePHibHVSoVTp8+bbJPqVSaUxoREVG7GuDXHT9fKsfhC9cxLcpH6nJsnlmBZfr06SY/L1u2DImJidi/f3+zgaVnz57o3r17k8cFQYCXl5c5pRAREXWo4X3c8Pn+izh4/prUpRDaMIZFr9cjOTkZVVVViImJabbtoEGD4O3tjfHjx2PHjh0NjldWViIgIAB+fn6YNm0aMjMzW/z9Wq0WGo3GZCMiImovw3q7AQB+KdKg/EatxNWQ2YElOzsbzs7OUCgUiI+PR2pqKiIiIhpt6+3tjdWrVyMlJQWbNm1CWFgYxo8fj507dxrbhIeHIykpCVu2bMHGjRuhVCoxcuRI5OTkNFtHQkIC1Gq1cfP39zf3UoiIiJrU00WJPh5OEEXg8EX2skhNEM18d7ZOp0NeXh7KysqQkpKCNWvWICMjo8nQcqfp06dDEARs2bKl0eMGgwGDBw/GmDFjsGrVqiY/R6vVQqvVGn/WaDTw9/dHeXm5yXgZIiKiu7UwORObsy5j8ZRwPDc2SOpyrJJGo4FarW7x+9usMSwAIJfLjYNuo6OjcejQIaxcuRKffPJJq84fMWIEvvjiiyaPy2QyDB06tMUeFoVCAYVC0frCiYiIzOStdgQAFJZz5qrU2rwOiyiKJj0dLcnMzIS3t3ezn5eVldVsGyIios7gpar/D+NiDQOL1MzqYVmyZAmmTJkCf39/VFRUIDk5Genp6UhLSwMALF68GAUFBdiwYQMA4IMPPkBgYCAiIyOh0+nwxRdfICUlBSkpKcbPXLp0KUaMGIGQkBBoNBqsWrUKWVlZ+PDDD9vxMomIiMznxR4Wi2FWYCkuLkZcXBwKCwuhVqsRFRWFtLQ0TJw4EQBQWFiIvLw8Y3udTodFixahoKAAjo6OiIyMxPfff4/777/f2KasrAzz589HUVER1Go1Bg0ahJ07d2LYsGHtdIlERER3x1tdvyZYYXm1xJWQ2YNuLVVrB+0QERG11pUKLYYu+xGCAPzyp8lQOthJXZLVae33N98lRERE1AQPZzmc5HYQRSCf7xSSFAMLERFRE25/p9D50iqJq7FtDCxERETNuBVYLlxlYJESAwsREVEzervf6mHhIyEpMbAQERE141YPy0X2sEiKgYWIiKgZvT26AQAucAyLpBhYiIiImhF485HQ5fIa1NTqJa7GdjGwEBERNcPNSQ4XZf06qxevchyLVBhYiIiImiEIAnpzarPkGFiIiIhacOuxEKc2S4eBhYiIqAXGtVjYwyIZBhYiIqIW3JopxEdC0mFgISIiakGA+621WDjoVioMLERERC24tdptkaYGN3R1EldjmxhYiIiIWuDqJIeHswIAkH2pXOJqbBMDCxERUSvEBLkDAPacuypxJbaJgYWIiKgVRgXfDCxnSyWuxDYxsBAREbVCbJAHAOBYfhkqtRzH0tkYWIiIiFrB360bfNRK1BlEnCrUSF2OzWFgISIiaiUvtRIAUFqpk7gS28PAQkRE1EruN2cKXa3SSlyJ7WFgISIiaiUPZzkA4Cp7WDodAwsREVEruTvd7GGpZA9LZ2NgISIiaiX3mz0spVXsYelsDCxEREStdGu12+LyGokrsT0MLERERK0U5uUCADhZqIHeIEpcjW1hYCEiImqloB7OcJLb4YZOj7MllVKXY1MYWIiIiFrJTiYg0lcNADhewJcgdiYGFiIiIjOEedY/FjpTUiFxJbaFgYWIiMgMoTfHseQU85FQZ2JgISIiMkNoT2cAwOki9rB0JgYWIiIiM4TefCRUUFbNtzZ3IgYWIiIiM7g6ydHDpX49lpxi9rJ0FrMCS2JiIqKioqBSqaBSqRATE4OtW7c22T49PR2CIDTYTp06ZdIuJSUFERERUCgUiIiIQGpq6t1dDRERUSe4NfCW41g6j1mBxc/PD8uXL8fhw4dx+PBhjBs3DjNmzMCJEyeaPe/06dMoLCw0biEhIcZj+/btw8yZMxEXF4djx44hLi4Ojz76KA4cOHB3V0RERNTBQjxvjmNhD0unEURRbNNSfW5ublixYgXmzZvX4Fh6ejruvfdeXL9+Hd27d2/0/JkzZ0Kj0Zj01EyePBmurq7YuHFjq+vQaDRQq9UoLy+HSqUy+zqIiIhaK/lgHl7flI3RIR74fN5wqcvp0lr7/X3XY1j0ej2Sk5NRVVWFmJiYZtsOGjQI3t7eGD9+PHbs2GFybN++fZg0aZLJvvvuuw979+5t9jO1Wi00Go3JRkRE1BluTW0+cVmDWr1B4mpsg9mBJTs7G87OzlAoFIiPj0dqaioiIiIabevt7Y3Vq1cjJSUFmzZtQlhYGMaPH4+dO3ca2xQVFcHT09PkPE9PTxQVFTVbR0JCAtRqtXHz9/c391KIiIjuSqSPCh7OClyr0uG7Y5elLscmmB1YwsLCkJWVhf3792PBggWYO3cuTp482WTbZ599FoMHD0ZMTAw++ugjTJ06FX/5y19M2gmCYPKzKIoN9t1p8eLFKC8vN275+fnmXgoREdFdUdjb4amRgQCA1Ttz0cbRFdQKZgcWuVyO4OBgREdHIyEhAQMGDMDKlStbff6IESOQk5Nj/NnLy6tBb0pJSUmDXpc7KRQK42ylWxsREVFnmT28FwDgVFEFNNVcj6WjtXkdFlEUodVqW90+MzMT3t7exp9jYmKwfft2kzbbtm1DbGxsW0sjIiLqMN27yeHmJAcAXC6vlrga62dvTuMlS5ZgypQp8Pf3R0VFBZKTk5Geno60tDQA9Y9pCgoKsGHDBgDABx98gMDAQERGRkKn0+GLL75ASkoKUlJSjJ+5cOFCjBkzBu+99x5mzJiBzZs348cff8Tu3bvb8TKJiIjan093Ja5V6XC5rBp9vdnT35HMCizFxcWIi4tDYWEh1Go1oqKikJaWhokTJwIACgsLkZeXZ2yv0+mwaNEiFBQUwNHREZGRkfj+++9x//33G9vExsYiOTkZb775Jt566y0EBQXh66+/xvDhnCZGRESWzVvtiOMFGlwur5G6FKvX5nVYLAXXYSEios72zpYTSNp7Ac+N7YPFU/pKXU6X1OHrsBAREdm68JvrsRw6f03iSqwfAwsREdFdGh3aAwCQlV+G8upaiauxbgwsREREd8m3uyOCejjBIAJ7z5ZKXY5VY2AhIiJqgzE3e1l25jCwdCQGFiIiojYYE1IfWHafvSJxJdaNgYWIiKgNhvV2g4OdgPxr1bh4tUrqcqwWAwsREVEbOCnsMaiXKwBgN8exdBgGFiIiojYaHewBANh2oljiSqwXAwsREVEb3R/lDTuZgIwzV7D3HHtZOgIDCxERURsF9XA2vr35b9vOSFyNdWJgISIiagcL7gkCABzJu46yGzqJq7E+DCxERETtwFvtiFBPZ4gisOfsVanLsToMLERERO3k1posu3K4Jkt7Y2AhIiJqJ7feLbTzzBWIoihxNdaFgYWIiKidDAt0g9xehsvlNTh3hYvItScGFiIionbiKLfDsEA3AED66RKJq7EuDCxERETtaELfngCA77MLJa7EujCwEBERtaP7+3tDEIDMvDJcun5D6nKsBgMLERFRO+qpUmJ47/rHQkl7LkhbjBVhYCEiImpnc2MCAQBrdp/H9z/z0VB7YGAhIiJqZ1P6e2PeqN4AgC/2X5S4GuvAwEJERNQBHhtW/26hI3nXUVOrl7iaro+BhYiIqAME9XCCl0oJXZ0BqZkFUpfT5TGwEBERdQBBEDB/TB8AQGL6ORgMXPm2LRhYiIiIOsisYf5QKe2Rd+0G0s9wIbm2YGAhIiLqIN3k9ng02h8AsH4vB9+2BQMLERFRB4qLCYAgABlnruBYfpnU5XRZDCxEREQdKMDdCRP7egIAHv90P05e1khcUdfEwEJERNTBVjwyAEMDXVGl02P93gtSl9MlMbAQERF1MLWjA34/IRQAsPV4IUortRJX1PUwsBAREXWC4X3cEe7lAk1NHf7831+kLqfLYWAhIiLqBHYyAQkP9QcAbDpagOxL5RJX1LUwsBAREXWSQb1cMWOgDwDg9//KgqamVuKKug4GFiIiok705tQIeKmUOFtSifkbDkNbx/cMtYZZgSUxMRFRUVFQqVRQqVSIiYnB1q1bW3Xunj17YG9vj4EDB5rsT0pKgiAIDbaamhpzSiMiIuoSergo8OmcaDjJ7bA/9xo2HeV7hlrDrMDi5+eH5cuX4/Dhwzh8+DDGjRuHGTNm4MSJE82eV15ejjlz5mD8+PGNHlepVCgsLDTZlEqlOaURERF1Gf391Hj55qyh99JOIfdKpcQVWT6zAsv06dNx//33IzQ0FKGhoVi2bBmcnZ2xf//+Zs977rnn8PjjjyMmJqbR44IgwMvLy2QjIiKyZk+MCECUnxplN2rx3OdHUKmtk7oki3bXY1j0ej2Sk5NRVVXVZBABgHXr1uHcuXN4++23m2xTWVmJgIAA+Pn5Ydq0acjMzGzx92u1Wmg0GpONiIioq3CU22HN3Gh4qhTIKanEc58fxjn2tDTJ7MCSnZ0NZ2dnKBQKxMfHIzU1FREREY22zcnJweuvv44vv/wS9vb2jbYJDw9HUlIStmzZgo0bN0KpVGLkyJHIyclpto6EhASo1Wrj5u/vb+6lEBERSaqnixIfPzEEcnsZ9py9iol/y8C/DudLXZZFEkRRFM05QafTIS8vD2VlZUhJScGaNWuQkZHRILTo9XqMGDEC8+bNQ3x8PADgnXfewbfffousrKwmP99gMGDw4MEYM2YMVq1a1WQ7rVYLrfbXlQI1Gg38/f1RXl4OlUplziURERFJ6lSRBu+nncZPp0rgrLDHj6+MhZfaNsZyajQaqNXqFr+/zQ4sd5owYQKCgoLwySefmOwvKyuDq6sr7OzsjPsMBgNEUYSdnR22bduGcePGNfqZzz77LC5dutTqGUhA6y+YiIjIEhkMIn778V5k5pUhOsAVn88bDke5XcsndnGt/f5u8zosoiia9HTcolKpkJ2djaysLOMWHx+PsLAwZGVlYfjw4U1+XlZWFry9vdtaGhERUZchkwl477dRcFHa4/DF6xi7YgcO5F6VuiyL0fjAkiYsWbIEU6ZMgb+/PyoqKpCcnIz09HSkpaUBABYvXoyCggJs2LABMpkM/fr1Mzm/Z8+eUCqVJvuXLl2KESNGICQkBBqNBqtWrUJWVhY+/PDDdrg8IiKiriPU0wWr46LxzPpDKKnQ4tkNh/GfF0ejl3s3qUuTnFk9LMXFxYiLi0NYWBjGjx+PAwcOIC0tDRMnTgQAFBYWIi8vz6wCysrKMH/+fPTt2xeTJk1CQUEBdu7ciWHDhpn1OURERNYgJsgdu/4wDsE9naGpqcMjn+xFsYaLqbZ5DIul4BgWIiKyJueuVGJe0iFcuHoDvdy6IWVBLHq4KKQuq9112hgWIiIian9BPZyx/ulh8O3uiLxrNzDqvZ+QcuSS1GVJhoGFiIjIQgW4O+HjJ4bAt7sjtHUGLEnNxtG861KXJQkGFiIiIgvW30+NXa/dizGhPaCtM2Du2oPIvlQudVmdjoGFiIjIwslkAj5+YjCGBbqhQluHuesO4ssDF1GnN0hdWqdhYCEiIuoCusntsfbJaPT3VeNalQ5vpB7H65uyoTdYxdyZFjGwEBERdREuSgckzx+BRZNCIQjAv49cwv/7VxYMNhBaGFiIiIi6ECeFPV4YF4IPZg6EvUzAt1mXsfS7E6jW6aUurUMxsBAREXVBMwb6YsUjUQCA9fsuInb5/3DkovXOIGJgISIi6qJ+M8gPf585AL7dHXH9Ri2e3XAYl8uqpS6rQzCwEBERdWG/GeSH7a+MQT9fFa5V6fC7L49CV2d9s4cYWIiIiLq4bnJ7JM4eApXSHln5ZXjz22xo66xrTAsDCxERkRXwd+uGD2YNBAD86/AlzPjnHpRX10pbVDtiYCEiIrIS48I9kTh7MNyc5DhVVIH/+89JqUtqNwwsREREVmRKf298OicaAJBy9BJ+KdRIXFH7YGAhIiKyMkMCXDE1yhuiCLz7/UmIYtdfWI6BhYiIyAr94b5wKOxl2HP2KpIP5UtdTpsxsBAREVmhXu7d8NL4EADA4k3Z2H6yWOKK2oaBhYiIyErNH9MHDwzwAQC89u9jyL92Q+KK7h4DCxERkZVysJPh/YejEOWnxvUbtXjz2+NSl3TXGFiIiIismNLBDqtmDYK9TEDGmSs4cblc6pLuCgMLERGRlQv0cMLkfl4AgM/3XZS4mrvDwEJERGQD5sQEAgC+zSpATnGFtMXcBQYWIiIiGzA00BUj+rihptaA3ybuRUEXe6szAwsREZENEAQB/3hsMMK9XKCpqcOfv/9F6pLMwsBCRERkI3q4KPC3RwdCJgDfZxfi+58LpS6p1RhYiIiIbEiEjwrPju4DoH7Zfr2hayzbz8BCRERkY34/MRTduzmgsLwGO89ckbqcVmFgISIisjFKBzs8NMgPAJCYcQ6amlqJK2oZAwsREZENeny4PxzsBBw8fw0z/rnH4pftZ2AhIiKyQcE9XfDFvOHwVitxvrQKv/loD8pu6KQuq0kMLERERDZqeB93fPXsCHiqFCit1GHt7vNSl9QkBhYiIiIb1tvDCUsfiAQArNtzwWLHszCwEBER2bhJEV4I7umMSm0dvjt2WepyGmVWYElMTERUVBRUKhVUKhViYmKwdevWVp27Z88e2NvbY+DAgQ2OpaSkICIiAgqFAhEREUhNTTWnLCIiImoDmUzArKH+AICVP+bgbInlvWvIrMDi5+eH5cuX4/Dhwzh8+DDGjRuHGTNm4MSJE82eV15ejjlz5mD8+PENju3btw8zZ85EXFwcjh07hri4ODz66KM4cOCAeVdCREREd+3x4b3g5+qIkgotZq3ej7yrljVrSBBFsU1L3Lm5uWHFihWYN29ek21mzZqFkJAQ2NnZ4dtvv0VWVpbx2MyZM6HRaEx6aiZPngxXV1ds3Lix1XVoNBqo1WqUl5dDpVLd1bUQERHZskvXb+DZDUfwS6EGfTycsOl3sejeTd6hv7O13993PYZFr9cjOTkZVVVViImJabLdunXrcO7cObz99tuNHt+3bx8mTZpksu++++7D3r17m/39Wq0WGo3GZCMiIqK75+faDUlPDYVvd0fkllbhgx9zpC7JyOzAkp2dDWdnZygUCsTHxyM1NRURERGNts3JycHrr7+OL7/8Evb29o22KSoqgqenp8k+T09PFBUVNVtHQkIC1Gq1cfP39zf3UoiIiOgOniol3n84CgDw5YGLFrOgnNmBJSwsDFlZWdi/fz8WLFiAuXPn4uTJkw3a6fV6PP7441i6dClCQ0Ob/UxBEEx+FkWxwb47LV68GOXl5cYtPz/f3EshIiKiRowM9sDIYHfU6kXMXXcQ//n5Mto4gqTNGu/2aIZcLkdwcDAAIDo6GocOHcLKlSvxySefmLSrqKjA4cOHkZmZiRdeeAEAYDAYIIoi7O3tsW3bNowbNw5eXl4NelNKSkoa9LrcSaFQQKFQmFs+ERERtcL/zeiH3ybuRe6VKrzwVSb+N6gE7z8cBQc7aVZEafNvFUURWq22wX6VSoXs7GxkZWUZt/j4eGMPzfDhwwEAMTEx2L59u8m527ZtQ2xsbFtLIyIiorvUp4cztv1+LF4aFwx7mYDUzAJ8c/iSZPWY1cOyZMkSTJkyBf7+/qioqEBycjLS09ORlpYGoP4xTUFBATZs2ACZTIZ+/fqZnN+zZ08olUqT/QsXLsSYMWPw3nvvYcaMGdi8eTN+/PFH7N69ux0uj4iIiO5WDxcFXpkUhoG9umP7yWLjWi1SMCuwFBcXIy4uDoWFhVCr1YiKikJaWhomTpwIACgsLEReXp5ZBcTGxiI5ORlvvvkm3nrrLQQFBeHrr7829sAQERGRtMaFe2JcePNDNTpam9dhsRRch4WIiKjr6fB1WIiIiIg6CwMLERERWTwGFiIiIrJ4DCxERERk8RhYiIiIyOIxsBAREZHFY2AhIiIii8fAQkRERBaPgYWIiIgsHgMLERERWTwGFiIiIrJ4DCxERERk8cx6W7Mlu/UOR41GI3ElRERE1Fq3vrdbehez1QSWiooKAIC/v7/ElRAREZG5KioqoFarmzwuiC1Fmi7CYDDg8uXLcHFxgSAI7fa5Go0G/v7+yM/Pb/a112Q5eM+6Ht6zrof3rGux5PsliiIqKirg4+MDmazpkSpW08Mik8ng5+fXYZ+vUqks7iZT83jPuh7es66H96xrsdT71VzPyi0cdEtEREQWj4GFiIiILB4DSwsUCgXefvttKBQKqUuhVuI963p4z7oe3rOuxRrul9UMuiUiIiLrxR4WIiIisngMLERERGTxGFiIiIjI4jGwEBERkcVjYGnBRx99hN69e0OpVGLIkCHYtWuX1CXZpISEBAwdOhQuLi7o2bMnHnzwQZw+fdqkjSiKeOedd+Dj4wNHR0fcc889OHHihEkbrVaLF198ER4eHnBycsIDDzyAS5cudeal2KSEhAQIgoCXX37ZuI/3y/IUFBTgiSeegLu7O7p164aBAwfiyJEjxuO8Z5alrq4Ob775Jnr37g1HR0f06dMHf/rTn2AwGIxtrOqeidSk5ORk0cHBQfz000/FkydPigsXLhSdnJzEixcvSl2azbnvvvvEdevWicePHxezsrLEqVOnir169RIrKyuNbZYvXy66uLiIKSkpYnZ2tjhz5kzR29tb1Gg0xjbx8fGir6+vuH37dvHo0aPivffeKw4YMECsq6uT4rJswsGDB8XAwEAxKipKXLhwoXE/75dluXbtmhgQECA++eST4oEDB8Tz58+LP/74o3j27FljG94zy/Luu++K7u7u4n/+8x/x/Pnz4jfffCM6OzuLH3zwgbGNNd0zBpZmDBs2TIyPjzfZFx4eLr7++usSVUS3lJSUiADEjIwMURRF0WAwiF5eXuLy5cuNbWpqakS1Wi1+/PHHoiiKYllZmejg4CAmJycb2xQUFIgymUxMS0vr3AuwERUVFWJISIi4fft2cezYscbAwvtlef7whz+Io0aNavI475nlmTp1qvj000+b7HvooYfEJ554QhRF67tnfCTUBJ1OhyNHjmDSpEkm+ydNmoS9e/dKVBXdUl5eDgBwc3MDAJw/fx5FRUUm90uhUGDs2LHG+3XkyBHU1taatPHx8UG/fv14TzvI888/j6lTp2LChAkm+3m/LM+WLVsQHR2NRx55BD179sSgQYPw6aefGo/znlmeUaNG4X//+x/OnDkDADh27Bh2796N+++/H4D13TOreflheystLYVer4enp6fJfk9PTxQVFUlUFQH1z2RfeeUVjBo1Cv369QMA4z1p7H5dvHjR2EYul8PV1bVBG97T9pecnIyjR4/i0KFDDY7xflme3NxcJCYm4pVXXsGSJUtw8OBBvPTSS1AoFJgzZw7vmQX6wx/+gPLycoSHh8POzg56vR7Lli3DY489BsD6/j1jYGmBIAgmP4ui2GAfda4XXngBP//8M3bv3t3g2N3cL97T9pefn4+FCxdi27ZtUCqVTbbj/bIcBoMB0dHR+POf/wwAGDRoEE6cOIHExETMmTPH2I73zHJ8/fXX+OKLL/DVV18hMjISWVlZePnll+Hj44O5c+ca21nLPeMjoSZ4eHjAzs6uQcIsKSlpkFap87z44ovYsmULduzYAT8/P+N+Ly8vAGj2fnl5eUGn0+H69etNtqH2ceTIEZSUlGDIkCGwt7eHvb09MjIysGrVKtjb2xv/9+b9shze3t6IiIgw2de3b1/k5eUB4L9jlujVV1/F66+/jlmzZqF///6Ii4vD73//eyQkJACwvnvGwNIEuVyOIUOGYPv27Sb7t2/fjtjYWImqsl2iKOKFF17Apk2b8NNPP6F3794mx3v37g0vLy+T+6XT6ZCRkWG8X0OGDIGDg4NJm8LCQhw/fpz3tJ2NHz8e2dnZyMrKMm7R0dGYPXs2srKy0KdPH94vCzNy5MgGSwWcOXMGAQEBAPjvmCW6ceMGZDLTr3E7OzvjtGaru2cSDfbtEm5Na167dq148uRJ8eWXXxadnJzECxcuSF2azVmwYIGoVqvF9PR0sbCw0LjduHHD2Gb58uWiWq0WN23aJGZnZ4uPPfZYo9P3/Pz8xB9//FE8evSoOG7cOIucvmeNbp8lJIq8X5bm4MGDor29vbhs2TIxJydH/PLLL8Vu3bqJX3zxhbEN75llmTt3rujr62uc1rxp0ybRw8NDfO2114xtrOmeMbC04MMPPxQDAgJEuVwuDh482DiNljoXgEa3devWGdsYDAbx7bffFr28vESFQiGOGTNGzM7ONvmc6upq8YUXXhDd3NxER0dHcdq0aWJeXl4nX41tujOw8H5Znu+++07s16+fqFAoxPDwcHH16tUmx3nPLItGoxEXLlwo9urVS1QqlWKfPn3EN954Q9RqtcY21nTPBFEURSl7eIiIiIhawjEsREREZPEYWIiIiMjiMbAQERGRxWNgISIiIovHwEJEREQWj4GFiIiILB4DCxEREVk8BhYiIiKyeAwsREREZPEYWIiIiMjiMbAQERGRxWNgISIiIov3/wFLw/KbtNAM8gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following corresponds to our model. We model $\\beta(t)$ as a fourier series of order 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ODE(nn.Module):\n",
    "\n",
    "    def __init__(self, y0, data, N):\n",
    "        super(ODE, self).__init__()\n",
    "\n",
    "        #initial parameters\n",
    "        self.gamma = 0.06\n",
    "        self.y0 = y0\n",
    "        self.data = tc.tensor(data).float()\n",
    "        self.N = N\n",
    "        self.T = self.data.numel()\n",
    "        self.tspan = tc.arange(0, self.T).float()\n",
    "\n",
    "        #beta(t) parameters\n",
    "        self.n = 10\n",
    "        self.a0 = nn.Parameter(tc.tensor(0.2))\n",
    "        self.a = nn.Parameter(tc.zeros(self.n - 1))\n",
    "        self.b = nn.Parameter(tc.zeros(self.n - 1))\n",
    "        self.sigmoid = nn.Sigmoid()\n",
    "\n",
    "        #minimum likelihood that can be achieved with the provided loss\n",
    "        self.min_likelihood = (1/N) * (self.data.mean() - (1/len(self.tspan)) * tc.dot(self.data.log(), self.data)) #when pred = data\n",
    "        \n",
    "\n",
    "    def betafn(self, t):\n",
    "        x = tc.full((self.n - 1,), t.item()) * 2 * tc.pi / self.T - self.T / 2 #scaled input t\n",
    "        nspan = tc.arange(1, self.n, dtype=float) #n values\n",
    "        x = nspan * x #element-wise mult.\n",
    "        return self.sigmoid((1/2) * self.a0 + (self.a * tc.cos(x) + self.b * tc.sin(x)).sum().squeeze())\n",
    "        \n",
    "    def forward(self, t, y):\n",
    "        S, I, R, N = y\n",
    "\n",
    "        beta = self.betafn(t)\n",
    "        dS = - beta * S * I / N\n",
    "        dI = beta * S * I / N - self.gamma * I\n",
    "        dR = self.gamma * I\n",
    "        dN = tc.tensor(0.0)\n",
    "\n",
    "        return tc.stack([dS, dI, dR, dN])\n",
    "    \n",
    "    def predict(self):\n",
    "        S, I, R, N = odeint(self, self.y0, self.tspan).T\n",
    "        return S, I, R, N\n",
    "    \n",
    "\n",
    "    def likelihood(self):\n",
    "        S, I, R, _ = self.predict()\n",
    "        return (1/self.N) * (S.mean() - (1/len(self.tspan)) * tc.dot(S.log(), self.data))\n",
    "    \n",
    "    def convergence_criteria(self, current_lik, tol=0.01):\n",
    "        min_lik = self.min_likelihood\n",
    "        return tc.isclose(current_lik, min_lik, atol=tol)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimization routine. Since we know the global minimum of the likelihood, we can establish a convergence criteria. We use adam to optimize objective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss : -8.968239784240723\n",
      "Loss : -14.196749687194824\n",
      "Loss : -14.633702278137207\n",
      "Loss : -14.778605461120605\n",
      "Loss : -14.805987358093262\n",
      "Loss : -14.91203784942627\n",
      "Loss : -14.89608097076416\n",
      "Loss : -14.902745246887207\n",
      "Loss : -15.005561828613281\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb Cell 10\u001b[0m in \u001b[0;36m<cell line: 8>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m \u001b[39mwhile\u001b[39;00m \u001b[39mnot\u001b[39;00m func\u001b[39m.\u001b[39mconvergence_criteria(current_lik\u001b[39m=\u001b[39mcurrent_lik, tol\u001b[39m=\u001b[39m\u001b[39m1e-5\u001b[39m):\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m     optimizer\u001b[39m.\u001b[39mzero_grad()\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=9'>10</a>\u001b[0m     l \u001b[39m=\u001b[39m func\u001b[39m.\u001b[39;49mlikelihood()\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=10'>11</a>\u001b[0m     current_lik \u001b[39m=\u001b[39m l\u001b[39m.\u001b[39mdetach()\u001b[39m.\u001b[39mclone()\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=11'>12</a>\u001b[0m     l\u001b[39m.\u001b[39mbackward()\n",
      "\u001b[1;32m/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb Cell 10\u001b[0m in \u001b[0;36mODE.likelihood\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=46'>47</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mlikelihood\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=47'>48</a>\u001b[0m     S, I, R, _ \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mpredict()\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=48'>49</a>\u001b[0m     \u001b[39mreturn\u001b[39;00m (\u001b[39m1\u001b[39m\u001b[39m/\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mN) \u001b[39m*\u001b[39m (S\u001b[39m.\u001b[39mmean() \u001b[39m-\u001b[39m (\u001b[39m1\u001b[39m\u001b[39m/\u001b[39m\u001b[39mlen\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtspan)) \u001b[39m*\u001b[39m tc\u001b[39m.\u001b[39mdot(S\u001b[39m.\u001b[39mlog(), \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata))\n",
      "\u001b[1;32m/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb Cell 10\u001b[0m in \u001b[0;36mODE.predict\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=41'>42</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mpredict\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=42'>43</a>\u001b[0m     S, I, R, N \u001b[39m=\u001b[39m odeint(\u001b[39mself\u001b[39;49m, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49my0, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtspan)\u001b[39m.\u001b[39mT\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=43'>44</a>\u001b[0m     \u001b[39mreturn\u001b[39;00m S, I, R, N\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/odeint.py:77\u001b[0m, in \u001b[0;36modeint\u001b[0;34m(func, y0, t, rtol, atol, method, options, event_fn)\u001b[0m\n\u001b[1;32m     74\u001b[0m solver \u001b[39m=\u001b[39m SOLVERS[method](func\u001b[39m=\u001b[39mfunc, y0\u001b[39m=\u001b[39my0, rtol\u001b[39m=\u001b[39mrtol, atol\u001b[39m=\u001b[39matol, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39moptions)\n\u001b[1;32m     76\u001b[0m \u001b[39mif\u001b[39;00m event_fn \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m---> 77\u001b[0m     solution \u001b[39m=\u001b[39m solver\u001b[39m.\u001b[39;49mintegrate(t)\n\u001b[1;32m     78\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m     79\u001b[0m     event_t, solution \u001b[39m=\u001b[39m solver\u001b[39m.\u001b[39mintegrate_until_event(t[\u001b[39m0\u001b[39m], event_fn)\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/solvers.py:30\u001b[0m, in \u001b[0;36mAdaptiveStepsizeODESolver.integrate\u001b[0;34m(self, t)\u001b[0m\n\u001b[1;32m     28\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_before_integrate(t)\n\u001b[1;32m     29\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39m1\u001b[39m, \u001b[39mlen\u001b[39m(t)):\n\u001b[0;32m---> 30\u001b[0m     solution[i] \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_advance(t[i])\n\u001b[1;32m     31\u001b[0m \u001b[39mreturn\u001b[39;00m solution\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/rk_common.py:194\u001b[0m, in \u001b[0;36mRKAdaptiveStepsizeODESolver._advance\u001b[0;34m(self, next_t)\u001b[0m\n\u001b[1;32m    192\u001b[0m \u001b[39mwhile\u001b[39;00m next_t \u001b[39m>\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrk_state\u001b[39m.\u001b[39mt1:\n\u001b[1;32m    193\u001b[0m     \u001b[39massert\u001b[39;00m n_steps \u001b[39m<\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mmax_num_steps, \u001b[39m'\u001b[39m\u001b[39mmax_num_steps exceeded (\u001b[39m\u001b[39m{}\u001b[39;00m\u001b[39m>=\u001b[39m\u001b[39m{}\u001b[39;00m\u001b[39m)\u001b[39m\u001b[39m'\u001b[39m\u001b[39m.\u001b[39mformat(n_steps, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mmax_num_steps)\n\u001b[0;32m--> 194\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrk_state \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_adaptive_step(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrk_state)\n\u001b[1;32m    195\u001b[0m     n_steps \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m1\u001b[39m\n\u001b[1;32m    196\u001b[0m \u001b[39mreturn\u001b[39;00m _interp_evaluate(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrk_state\u001b[39m.\u001b[39minterp_coeff, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrk_state\u001b[39m.\u001b[39mt0, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrk_state\u001b[39m.\u001b[39mt1, next_t)\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/rk_common.py:255\u001b[0m, in \u001b[0;36mRKAdaptiveStepsizeODESolver._adaptive_step\u001b[0;34m(self, rk_state)\u001b[0m\n\u001b[1;32m    250\u001b[0m         dt \u001b[39m=\u001b[39m t1 \u001b[39m-\u001b[39m t0\n\u001b[1;32m    252\u001b[0m \u001b[39m# Must be arranged as doing all the step_t handling, then all the jump_t handling, in case we\u001b[39;00m\n\u001b[1;32m    253\u001b[0m \u001b[39m# trigger both. (i.e. interleaving them would be wrong.)\u001b[39;00m\n\u001b[0;32m--> 255\u001b[0m y1, f1, y1_error, k \u001b[39m=\u001b[39m _runge_kutta_step(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfunc, y0, f0, t0, dt, t1, tableau\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtableau)\n\u001b[1;32m    256\u001b[0m \u001b[39m# dtypes:\u001b[39;00m\n\u001b[1;32m    257\u001b[0m \u001b[39m# y1.dtype == self.y0.dtype\u001b[39;00m\n\u001b[1;32m    258\u001b[0m \u001b[39m# f1.dtype == self.y0.dtype\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    263\u001b[0m \u001b[39m#                     Error Ratio                      #\u001b[39;00m\n\u001b[1;32m    264\u001b[0m \u001b[39m########################################################\u001b[39;00m\n\u001b[1;32m    265\u001b[0m error_ratio \u001b[39m=\u001b[39m _compute_error_ratio(y1_error, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrtol, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39matol, y0, y1, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mnorm)\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/rk_common.py:76\u001b[0m, in \u001b[0;36m_runge_kutta_step\u001b[0;34m(func, y0, f0, t0, dt, t1, tableau)\u001b[0m\n\u001b[1;32m     74\u001b[0m         perturb \u001b[39m=\u001b[39m Perturb\u001b[39m.\u001b[39mNONE\n\u001b[1;32m     75\u001b[0m     yi \u001b[39m=\u001b[39m y0 \u001b[39m+\u001b[39m torch\u001b[39m.\u001b[39msum(k[\u001b[39m.\u001b[39m\u001b[39m.\u001b[39m\u001b[39m.\u001b[39m, :i \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m] \u001b[39m*\u001b[39m (beta_i \u001b[39m*\u001b[39m dt), dim\u001b[39m=\u001b[39m\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m)\u001b[39m.\u001b[39mview_as(f0)\n\u001b[0;32m---> 76\u001b[0m     f \u001b[39m=\u001b[39m func(ti, yi, perturb\u001b[39m=\u001b[39;49mperturb)\n\u001b[1;32m     77\u001b[0m     k \u001b[39m=\u001b[39m _UncheckedAssign\u001b[39m.\u001b[39mapply(k, f, (\u001b[39m.\u001b[39m\u001b[39m.\u001b[39m\u001b[39m.\u001b[39m, i \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m))\n\u001b[1;32m     79\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m (tableau\u001b[39m.\u001b[39mc_sol[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m] \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m \u001b[39mand\u001b[39;00m (tableau\u001b[39m.\u001b[39mc_sol[:\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m] \u001b[39m==\u001b[39m tableau\u001b[39m.\u001b[39mbeta[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m])\u001b[39m.\u001b[39mall()):\n\u001b[1;32m     80\u001b[0m     \u001b[39m# This property (true for Dormand-Prince) lets us save a few FLOPs.\u001b[39;00m\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torch/nn/modules/module.py:1102\u001b[0m, in \u001b[0;36mModule._call_impl\u001b[0;34m(self, *input, **kwargs)\u001b[0m\n\u001b[1;32m   1098\u001b[0m \u001b[39m# If we don't have any hooks, we want to skip the rest of the logic in\u001b[39;00m\n\u001b[1;32m   1099\u001b[0m \u001b[39m# this function, and just call forward.\u001b[39;00m\n\u001b[1;32m   1100\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m (\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_backward_hooks \u001b[39mor\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_forward_hooks \u001b[39mor\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_forward_pre_hooks \u001b[39mor\u001b[39;00m _global_backward_hooks\n\u001b[1;32m   1101\u001b[0m         \u001b[39mor\u001b[39;00m _global_forward_hooks \u001b[39mor\u001b[39;00m _global_forward_pre_hooks):\n\u001b[0;32m-> 1102\u001b[0m     \u001b[39mreturn\u001b[39;00m forward_call(\u001b[39m*\u001b[39;49m\u001b[39minput\u001b[39;49m, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m   1103\u001b[0m \u001b[39m# Do not call functions when jit is used\u001b[39;00m\n\u001b[1;32m   1104\u001b[0m full_backward_hooks, non_full_backward_hooks \u001b[39m=\u001b[39m [], []\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torchdiffeq/_impl/misc.py:189\u001b[0m, in \u001b[0;36m_PerturbFunc.forward\u001b[0;34m(self, t, y, perturb)\u001b[0m\n\u001b[1;32m    186\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m    187\u001b[0m     \u001b[39m# Do nothing.\u001b[39;00m\n\u001b[1;32m    188\u001b[0m     \u001b[39mpass\u001b[39;00m\n\u001b[0;32m--> 189\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mbase_func(t, y)\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torch/nn/modules/module.py:1102\u001b[0m, in \u001b[0;36mModule._call_impl\u001b[0;34m(self, *input, **kwargs)\u001b[0m\n\u001b[1;32m   1098\u001b[0m \u001b[39m# If we don't have any hooks, we want to skip the rest of the logic in\u001b[39;00m\n\u001b[1;32m   1099\u001b[0m \u001b[39m# this function, and just call forward.\u001b[39;00m\n\u001b[1;32m   1100\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m (\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_backward_hooks \u001b[39mor\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_forward_hooks \u001b[39mor\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_forward_pre_hooks \u001b[39mor\u001b[39;00m _global_backward_hooks\n\u001b[1;32m   1101\u001b[0m         \u001b[39mor\u001b[39;00m _global_forward_hooks \u001b[39mor\u001b[39;00m _global_forward_pre_hooks):\n\u001b[0;32m-> 1102\u001b[0m     \u001b[39mreturn\u001b[39;00m forward_call(\u001b[39m*\u001b[39;49m\u001b[39minput\u001b[39;49m, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m   1103\u001b[0m \u001b[39m# Do not call functions when jit is used\u001b[39;00m\n\u001b[1;32m   1104\u001b[0m full_backward_hooks, non_full_backward_hooks \u001b[39m=\u001b[39m [], []\n",
      "\u001b[1;32m/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb Cell 10\u001b[0m in \u001b[0;36mODE.forward\u001b[0;34m(self, t, y)\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=30'>31</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mforward\u001b[39m(\u001b[39mself\u001b[39m, t, y):\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=31'>32</a>\u001b[0m     S, I, R, N \u001b[39m=\u001b[39m y\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=33'>34</a>\u001b[0m     beta \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mbetafn(t)\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/william/Projects/ODEAutodiff/ode_model-latent-contact-rate.ipynb#W4sZmlsZQ%3D%3D?line=34'>35</a>\u001b[0m     dS \u001b[39m=\u001b[39m \u001b[39m-\u001b[39m beta \u001b[39m*\u001b[39m S \u001b[39m*\u001b[39m I \u001b[39m/\u001b[39m N\n",
      "File \u001b[0;32m~/opt/miniconda3/envs/ODEAutodiff/lib/python3.10/site-packages/torch/_tensor.py:648\u001b[0m, in \u001b[0;36mTensor.__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    643\u001b[0m \u001b[39mif\u001b[39;00m torch\u001b[39m.\u001b[39m_C\u001b[39m.\u001b[39m_get_tracing_state():\n\u001b[1;32m    644\u001b[0m     warnings\u001b[39m.\u001b[39mwarn(\u001b[39m'\u001b[39m\u001b[39mIterating over a tensor might cause the trace to be incorrect. \u001b[39m\u001b[39m'\u001b[39m\n\u001b[1;32m    645\u001b[0m                   \u001b[39m'\u001b[39m\u001b[39mPassing a tensor of different shape won\u001b[39m\u001b[39m\\'\u001b[39;00m\u001b[39mt change the number of \u001b[39m\u001b[39m'\u001b[39m\n\u001b[1;32m    646\u001b[0m                   \u001b[39m'\u001b[39m\u001b[39miterations executed (and might lead to errors or silently give \u001b[39m\u001b[39m'\u001b[39m\n\u001b[1;32m    647\u001b[0m                   \u001b[39m'\u001b[39m\u001b[39mincorrect results).\u001b[39m\u001b[39m'\u001b[39m, category\u001b[39m=\u001b[39mtorch\u001b[39m.\u001b[39mjit\u001b[39m.\u001b[39mTracerWarning, stacklevel\u001b[39m=\u001b[39m\u001b[39m2\u001b[39m)\n\u001b[0;32m--> 648\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39miter\u001b[39m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49munbind(\u001b[39m0\u001b[39;49m))\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "func = ODE(y0 = y0, data = data, N = N)\n",
    "\n",
    "optimizer = tc.optim.Adam(func.parameters(), lr=0.1)\n",
    "\n",
    "losses = []\n",
    "i = 0\n",
    "current_lik = tc.tensor(0.0)\n",
    "while not func.convergence_criteria(current_lik=current_lik, tol=1e-5):\n",
    "    optimizer.zero_grad()\n",
    "    l = func.likelihood()\n",
    "    current_lik = l.detach().clone()\n",
    "    l.backward()\n",
    "    if i%5==0 : print(\"Loss : {}\".format(current_lik))\n",
    "    optimizer.step()\n",
    "    losses.append(l.clone().detach())\n",
    "    i+=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum Likelihood : -15.92385482788086, Achieved Likelihood : -15.923604965209961, Difference : 0.0002498626708984375\n"
     ]
    }
   ],
   "source": [
    "print(\"Minimum Likelihood : {}, Achieved Likelihood : {}, Difference : {}\".format(func.min_likelihood, current_lik, tc.abs(func.min_likelihood - current_lik)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Loss')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGwCAYAAACpYG+ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4mUlEQVR4nO3de3TU1b3H/c/MJJncJ5AhN4gEDMhNKcWqoNZ4gWLVqlRbsLX42MpBG49oPecp1UKwh4OtB9t69NDaooFVezzVqo+trRZrtCreUFAEiiBKAkm4hdyTyWV+zx/JTAi5MiS/ncm8X2v9FsnMnsk3e40rH/dvXxyWZVkCAACIAE7TBQAAANiF4AMAACIGwQcAAEQMgg8AAIgYBB8AABAxCD4AACBiEHwAAEDEiDJdwFDj9/tVWlqqpKQkORwO0+UAAIB+sCxLNTU1ysrKktPZ87gOwecEpaWlys7ONl0GAAAIQUlJicaMGdPj8wSfEyQlJUlq67jk5GTD1QAAgP6orq5WdnZ28O94Twg+Jwjc3kpOTib4AAAQZvqapsLkZgAAEDEIPgAAIGIQfAAAQMQg+AAAgIhB8AEAABGD4AMAACIGwQcAAEQMgg8AAIgYBB8AABAxCD4AACBiEHwAAEDEIPgAAICIMayCzwcffKA5c+YoJSVFqampWrx4sWpra02XJUk6WN2okop6NTa3mi4FAICINWyCT2lpqS677DLl5ubqnXfe0Ysvvqjt27frpptuMl2aJGnBo2/rwp8VaduBKtOlAAAQsaJMFzBQ/vznPys6OlqPPPKInM62PPfII49oxowZ2rNnj3Jzc7t9nc/nk8/nC35fXV09KPU5HG3/tvqtQXl/AADQt2Ez4uPz+RQTExMMPZIUFxcnSXrjjTd6fN3q1avl8XiCV3Z29qDU52pPPn6L4AMAgCnDJvhccsklKi8v1wMPPKCmpiYdO3ZMP/rRjyRJZWVlPb5u2bJlqqqqCl4lJSWDUp8zEHz8g/L2AACgH4Z88CkoKJDD4ej12rx5s6ZOnar169drzZo1io+PV0ZGhsaPH6/09HS5XK4e39/tdis5ObnTNRicTkZ8AAAwbcjP8cnPz9eCBQt6bZOTkyNJuuGGG3TDDTfo4MGDSkhIkMPh0IMPPqhx48bZUGnvnIE5PgQfAACMGfLBx+v1yuv1ntRr0tPTJUmPPfaYYmNjNWfOnMEo7aS4AiM+TG4GAMCYIR98TsbDDz+s2bNnKzExURs3btS//du/6f7771dKSorp0jrm+JB7AAAwZlgFn3fffVcrVqxQbW2tJk2apF//+te68cYbTZcl6bhbXSQfAACMGVbBZ8OGDaZL6FHgVpfFHB8AAIwZ8qu6hgtH+60uJjcDAGAOwccmLub4AABgHMHHJoENpVnVBQCAOQQfmzg5sgIAAOMIPjYJBB9WdQEAYA7BxyYdq7oMFwIAQAQj+NiEIysAADCP4GMT5vgAAGAewccmweDDHB8AAIwh+NgkeEgpuQcAAGMIPjZxcFYXAADGEXxs0jHiQ/ABAMAUgo9NXExuBgDAOIKPTRyc1QUAgHEEH5u42nuaOT4AAJhD8LEJy9kBADCP4GMTJ8vZAQAwjuBjE46sAADAPIKPTQKruiyCDwAAxhB8bBJY1cXkZgAAzCH42IQjKwAAMI/gY5PAHB82MAQAwByCj02Cq7oY8gEAwBiCj00C+/iwqgsAAHMIPjbpWNVluBAAACIYwccmwX18uNUFAIAxBB+bdOzcTPABAMAUgo9NXA6CDwAAphF8bNKxqstwIQAARDCCj01Y1QUAgHkEH5uwgSEAAOYRfGziYgNDAACMI/jYJHhIKbkHAABjCD42cXGrCwAA4wg+NuGsLgAAzCP42MTJPj4AABhH8LFJcDk7+/gAAGBM2ASfVatWafbs2YqPj1dKSkq3bYqLi3XVVVcpISFBXq9X//qv/6qmpiZ7C+2Bq72nLUZ8AAAwJsp0Af3V1NSk66+/XrNmzdK6deu6PN/a2qorrrhCo0aN0htvvKGjR49q0aJFsixL//3f/22g4s4cbGAIAIBxYRN8Vq5cKUkqLCzs9vm//e1v2rFjh0pKSpSVlSVJWrNmjW666SatWrVKycnJdpXarY6zuoyWAQBARAubW119eeuttzRt2rRg6JGkr3zlK/L5fHr//fd7fJ3P51N1dXWnazA423uaVV0AAJgzbIJPeXm50tPTOz02YsQIxcTEqLy8vMfXrV69Wh6PJ3hlZ2cPSn2s6gIAwDyjwaegoEAOh6PXa/Pmzf1+v8A8muNZltXt4wHLli1TVVVV8CopKQnpd+lL4MiKVkZ8AAAwxugcn/z8fC1YsKDXNjk5Of16r4yMDL3zzjudHjt27Jiam5u7jAQdz+12y+129+tnnIrAiA8DPgAAmGM0+Hi9Xnm93gF5r1mzZmnVqlUqKytTZmampLYJz263WzNnzhyQn3EqnKzqAgDAuLBZ1VVcXKyKigoVFxertbVVW7dulSTl5uYqMTFRc+fO1ZQpU3TjjTfqgQceUEVFhe6++27dcsstxld0SZKTs7oAADAubILP8uXLtX79+uD3M2bMkCQVFRUpLy9PLpdLL7zwgm677Tadf/75iouL0w033KD/+q//MlVyJy7O6gIAwLiwCT6FhYU97uETcNppp+nPf/6zPQWdJCf7+AAAYNywWc4+1DlZ1QUAgHEEH5swxwcAAPMIPjZxsYEhAADGEXxsEjyklFtdAAAYQ/CxSWBVFwM+AACYQ/CxSWCODxsYAgBgDsHHJoFVXczxAQDAHIKPTYL7+PgNFwIAQAQj+NiEVV0AAJhH8LGJs72nWdUFAIA5BB+bcGQFAADmEXxs4mJyMwAAxhF8bMKRFQAAmEfwsYmTnZsBADCO4GOTQPBhwAcAAHMIPjYJzPFhxAcAAHMIPjZxMMcHAADjCD42YVUXAADmEXxswj4+AACYR/CxCau6AAAwj+Bjk8A+PpLkJ/wAAGAEwccmruOSD/N8AAAwg+BjE4ejI/i0EnwAADCC4GOT40d8yD0AAJhB8LGJ6/gRH+b4AABgBMHHJo7jJzcz5AMAgBEEH5t0mtzsN1gIAAARjOBjE6eDVV0AAJhG8LHJ8fv4sKoLAAAzCD42cTgcHFQKAIBhBB8bBVZ2MccHAAAzCD426jiolBEfAABMIPjYyNne2+zjAwCAGQQfGwVGfBjwAQDADIKPjQJzfFjVBQCAGQQfG7GqCwAAswg+Ngrs3uxnjg8AAEaETfBZtWqVZs+erfj4eKWkpHTb5o477tDMmTPldrv1hS98wdb6+qNjVZfhQgAAiFBhE3yampp0/fXX69Zbb+2xjWVZuvnmm/XNb37Txsr6z9k+4sOqLgAAzIgyXUB/rVy5UpJUWFjYY5uHHnpIknT48GF99NFHdpR1Ulzs4wMAgFFhE3wGi8/nk8/nC35fXV09aD/LyeRmAACMCptbXYNl9erV8ng8wSs7O3vQfha3ugAAMMto8CkoKGg/vLPna/PmzYNaw7Jly1RVVRW8SkpKBu1nMbkZAACzjN7qys/P14IFC3ptk5OTM6g1uN1uud3uQf0ZAcHl7NzqAgDACKPBx+v1yuv1mizBVsENDBnyAQDAiLCZ3FxcXKyKigoVFxertbVVW7dulSTl5uYqMTFRkrRnzx7V1taqvLxcDQ0NwTZTpkxRTEyMoco7cGQFAABmhU3wWb58udavXx/8fsaMGZKkoqIi5eXlSZK+973v6bXXXuvS5rPPPhv0W2b9wSGlAACYFTbBp7CwsNc9fCTp1VdftaWWULGqCwAAsyJ+Obud2McHAACzCD42YlUXAABmEXxs5Ajs4+M3XAgAABGK4GMjV/utLlZ1AQBgBsHHRh2rugg+AACYQPCxUceqLsOFAAAQoQg+NnI5mNwMAIBJBB8bOdt7m+ADAIAZBB8bBeb4sIEhAABmEHxs5Aze6jJcCAAAEYrgY6PgBoYkHwAAjCD42IgjKwAAMIvgY6PgHB+CDwAARhB8bMQcHwAAzCL42Ig5PgAAmEXwsZGDOT4AABhF8LGRy8k+PgAAmETwsVHHIaWGCwEAIEIRfGzEqi4AAMwi+NiIfXwAADCL4GMjVnUBAGAWwcdGDvbxAQDAKIKPjVztvc2qLgAAzCD42MgVXNVF8AEAwASCj40crOoCAMAogo+NgpObyT0AABhB8LFRcDk7yQcAACMIPjZycmQFAABGEXxs5GQ5OwAARhF8bOQKBh+SDwAAJhB8bMSRFQAAmEXwsRFzfAAAMIvgYyPm+AAAYBbBx0YcUgoAgFkEHxs5mOMDAIBRBB8buTiyAgAAowg+NnIGDyk1XAgAABGK4GMjVnUBAGBW2ASfVatWafbs2YqPj1dKSkqX5z/88EMtXLhQ2dnZiouL0+TJk/XLX/7S/kJ74WKODwAARkWZLqC/mpqadP3112vWrFlat25dl+fff/99jRo1Sr/73e+UnZ2tTZs2afHixXK5XMrPzzdQcVdOJzs3AwBgUtgEn5UrV0qSCgsLu33+5ptv7vT9+PHj9dZbb+mZZ57pNfj4fD75fL7g99XV1adebA+C+/j4B+1HAACAXoTNra5QVFVVaeTIkb22Wb16tTweT/DKzs4etHqcrOoCAMCoYRt83nrrLf3hD3/Qv/zLv/TabtmyZaqqqgpeJSUlg1aTq723LYIPAABGGA0+BQUFcjgcvV6bN28+6ffdvn27rr76ai1fvlxz5szpta3b7VZycnKna7A4HKzqAgDAJKNzfPLz87VgwYJe2+Tk5JzUe+7YsUOXXHKJbrnlFt17772nUN3A69jA0HAhAABEKKPBx+v1yuv1Dtj7bd++XZdccokWLVqkVatWDdj7DhQnt7oAADAqbFZ1FRcXq6KiQsXFxWptbdXWrVslSbm5uUpMTNT27dt18cUXa+7cubrrrrtUXl4uSXK5XBo1apTByjs4udUFAIBRYRN8li9frvXr1we/nzFjhiSpqKhIeXl5euqpp3T48GE98cQTeuKJJ4Ltxo4dq88//9zucrsVXM7OiA8AAEaEzaquwsJCWZbV5crLy5PUNlG6u+eHSuiRJJeTfXwAADApbILPcODkyAoAAIwi+NiIDQwBADCL4GOjjjk+hgsBACBChRR8SkpKtH///uD37777rpYuXapHH310wAobjjrm+JB8AAAwIaTgc8MNN6ioqEiSVF5erjlz5ujdd9/Vj370I913330DWuBwwunsAACYFVLw+fjjj3XOOedIkv7whz9o2rRp2rRpk37/+9/3eHo6OiY3s48PAABmhBR8mpub5Xa7JUkvv/yyvva1r0mSJk2apLKysoGrbpgJHFnBgA8AAGaEFHymTp2qX/3qV3r99de1ceNGzZs3T5JUWlqq1NTUAS1wOHGwqgsAAKNCCj4//elP9etf/1p5eXlauHChpk+fLkl6/vnng7fA0JWLOT4AABgV0pEVeXl5OnLkiKqrqzVixIjg44sXL1Z8fPyAFTfcBDcwZI4PAABGhDTi09DQIJ/PFww9+/bt0y9+8Qvt2rVLaWlpA1rgcNKxqstwIQAARKiQgs/VV1+tDRs2SJIqKyt17rnnas2aNbrmmmu0du3aAS1wOOF0dgAAzAop+HzwwQe68MILJUlPP/200tPTtW/fPm3YsEEPPfTQgBY4nHSs6iL4AABgQkjBp76+XklJSZKkv/3tb5o/f76cTqfOO+887du3b0ALHE4cgX18CD4AABgRUvDJzc3Vc889p5KSEr300kuaO3euJOnQoUNKTk4e0AKHE3dUW3c3tfgNVwIAQGQKKfgsX75cd999t3JycnTOOedo1qxZktpGf2bMmDGgBQ4n7iiXJMlH8AEAwIiQlrNfd911uuCCC1RWVhbcw0eSLr30Ul177bUDVtxwExvdljMbm1tlWVZwQ0MAAGCPkIKPJGVkZCgjI0P79++Xw+HQ6NGj2bywD4ERH78ltfgtRbsIPgAA2CmkW11+v1/33XefPB6Pxo4dq9NOO00pKSn6yU9+Ir+f2zg9cUd3dHdjc6vBSgAAiEwhjfjcc889Wrdune6//36df/75sixLb775pgoKCtTY2KhVq1YNdJ3DQmBys9Q2zyfJYC0AAESikILP+vXr9dvf/jZ4KrskTZ8+XaNHj9Ztt91G8OmBw+GQO8opX4ufER8AAAwI6VZXRUWFJk2a1OXxSZMmqaKi4pSLGs5io1nZBQCAKSEFn+nTp+vhhx/u8vjDDz+ss84665SLGs4Ct7sY8QEAwH4h3er62c9+piuuuEIvv/yyZs2aJYfDoU2bNqmkpER/+ctfBrrGYSUw4tPYzIgPAAB2C2nE56KLLtInn3yia6+9VpWVlaqoqND8+fO1fft2Pf744wNd47ASGPHxtTDiAwCA3ULexycrK6vLJOYPP/xQ69ev12OPPXbKhQ1XwTk+jPgAAGC7kEZ8EDpGfAAAMIfgYzPm+AAAYA7Bx2aM+AAAYM5JzfGZP39+r89XVlaeSi0RgREfAADMOang4/F4+nz+O9/5zikVNNwx4gMAgDknFXxYqn7q3Iz4AABgDHN8bMaIDwAA5hB8bMYcHwAAzCH42IwRHwAAzCH42IwRHwAAzCH42IzT2QEAMCdsgs+qVas0e/ZsxcfHKyUlpcvzR48e1bx585SVlSW3263s7Gzl5+erurra/mJ7ETyrq4URHwAA7BY2waepqUnXX3+9br311m6fdzqduvrqq/X888/rk08+UWFhoV5++WUtWbLE5kp7x4gPAADmhHw6u91WrlwpSSosLOz2+REjRnQKRWPHjtVtt92mBx54wI7y+o0RHwAAzAmb4HOySktL9cwzz+iiiy7qtZ3P55PP5wt+P9i3xoKruhjxAQDAdmFzq6u/Fi5cqPj4eI0ePVrJycn67W9/22v71atXy+PxBK/s7OxBrY8RHwAAzDEafAoKCuRwOHq9Nm/efFLv+fOf/1wffPCBnnvuOX366ae66667em2/bNkyVVVVBa+SkpJT+ZX6FBvNHB8AAEwxeqsrPz9fCxYs6LVNTk7OSb1nRkaGMjIyNGnSJKWmpurCCy/Uj3/8Y2VmZnbb3u12y+12n9TPOBXuKEZ8AAAwxWjw8Xq98nq9g/b+lmVJUqc5PKYx4gMAgDlhM7m5uLhYFRUVKi4uVmtrq7Zu3SpJys3NVWJiov7yl7/o4MGD+tKXvqTExETt2LFD//7v/67zzz//pEeNBhMjPgAAmBM2wWf58uVav3598PsZM2ZIkoqKipSXl6e4uDj95je/0Z133imfz6fs7GzNnz9fP/zhD02V3K3jR3wsy5LD4TBcEQAAkcNhBe4HQVLbcnaPx6OqqiolJycP+PtX1Tdr+n1/kyTtXnW5ol3DbmEdAAC26+/fb/7q2swd3dHlzPMBAMBeBB+bBTYwlDihHQAAuxF8bOZwODp2b25hxAcAADsRfAzoOKiUER8AAOxE8DGg49gKRnwAALATwccAdzQjPgAAmEDwMSA2ihEfAABMIPgYEBjx8THiAwCArQg+BjDiAwCAGQQfA5jjAwCAGQQfAxjxAQDADIKPAYz4AABgBsHHAEZ8AAAwg+BjACM+AACYQfAxwN0+4sPp7AAA2IvgY0BcTFvwqW8i+AAAYCeCjwEj4qMlSZX1TYYrAQAgshB8DBiZ4JYkHa0j+AAAYCeCjwGpCTGSpAqCDwAAtiL4GDCS4AMAgBEEHwMCwedoXZMsyzJcDQAAkYPgY0BqYlvwaWrxq9bXYrgaAAAiB8HHgPiYKMW2b2LI7S4AAOxD8DEklZVdAADYjuBjSOB2V0UtwQcAALsQfAxhZRcAAPYj+Bhy/MouAABgD4KPIR2bGPoMVwIAQOQg+BgSPLaCOT4AANiG4GNIYHIzt7oAALAPwccQzusCAMB+BB9DWNUFAID9CD6GdGxgyORmAADsQvAxZGT7HJ/GZr/qmzivCwAAOxB8DEmIcSkmqq37WdkFAIA9CD6GOBwOedvn+Ryp5XYXAAB2IPgYlO6JlSSVVzUargQAgMgQNsFn1apVmj17tuLj45WSktJr26NHj2rMmDFyOByqrKy0pb5QZHniJEmlBB8AAGwRNsGnqalJ119/vW699dY+2373u9/VWWedZUNVpyazfcSnrLLBcCUAAESGsAk+K1eu1J133qkzzzyz13Zr165VZWWl7r77bpsqC11mStuITxkjPgAA2CLKdAEDaceOHbrvvvv0zjvvaO/evf16jc/nk8/XMbm4urp6sMrrIisw4lPFiA8AAHYImxGfvvh8Pi1cuFAPPPCATjvttH6/bvXq1fJ4PMErOzt7EKvsjBEfAADsZTT4FBQUyOFw9Hpt3ry5X++1bNkyTZ48Wd/+9rdPqoZly5apqqoqeJWUlITyq4QkMMfnYHWjWlr9tv1cAAAildFbXfn5+VqwYEGvbXJycvr1Xq+88oq2bdump59+WpJkWZYkyev16p577tHKlSu7fZ3b7Zbb7e5/0QPIm+hWlNOhFr+lQzU+ZbWPAAEAgMFhNPh4vV55vd4Bea8//vGPamjomCvz3nvv6eabb9brr7+u008/fUB+xkBzOR1KT47VgcoGlVU1EHwAABhkYTO5ubi4WBUVFSouLlZra6u2bt0qScrNzVViYmKXcHPkyBFJ0uTJk/vc98ekrJS24FNa2aiZY01XAwDA8BY2wWf58uVav3598PsZM2ZIkoqKipSXl2eoqlOX6YmTdIyVXQAA2CBsVnUVFhbKsqwuV0+hJy8vT5ZlDenRHknKTGmb4FxaycouAAAGW9gEn+EqcGwFIz4AAAw+go9hgSXtmz8/ptV/3amjnNQOAMCgIfgYdkZGkhwO6Whdk3792l6t3/S56ZIAABi2CD6GjU1N0J/yL9BV07MkSSXHuOUFAMBgIfgMAdNGe3TppDRJUjnHVwAAMGgIPkNEenL78RU1BB8AAAYLwWeISE9uOzbjICM+AAAMGoLPEBEY8alralWtr8VwNQAADE8EnyEiwR2lJHfbRtrM8wEAYHAQfIaQ9PY9fQ5VE3wAABgMBJ8hJDDPp5zgAwDAoCD4DCHBlV3V7N4MAMBgIPgMIRnB4MOIDwAAg4HgM4SkE3wAABhUBJ8hJBB8mOMDAMDgIPgMIWxiCADA4CL4DCEZgeXsNT7945PDKq3kwFIAAAYSwWcI8Sa65XBILX5L33nsXd30+LuyLMt0WQAADBsEnyEk2uVUlicu+P0nB2v18YFqgxUBADC8EHyGmF8u+IJWXTtNl01OkyT96aNSwxUBADB8EHyGmLNzRupb547VdTPHSJJe+KhMfj+3uwAAGAgEnyEq74w0JbqjdKCyQVtKjpkuBwCAYYHgM0TFRrs0Z0q6JOmZDw7Isixt3HFQxUfrDVcGAED4IvgMYV//Ytvtruc/LNVT7+/XLRs2a+Fv3pavpdVwZQAAhCeCzxA2+/RUjU6JU01ji+599mNJ0oHKBv3u7WLDlQEAEJ4IPkOY0+kITnJuavXL6Wh7/JGiPappbDZYGQAA4YngM8QFgo8k/fu8SRrvTVBFXZMef/Nzc0UBABCmCD5DXPbIeC29bIKuPCtTN83O0R2XTZAkFW76XA1NzPUBAOBkRJkuAH1betnE4NdXnJmp//rbLpVUNOip90t043lj9fruI9r06VF969zTlD0yXjvLqvXc1gP6qKRKmZ5YnTc+VdfNHCNn4F4ZAAARymFxGFQn1dXV8ng8qqqqUnJysulyurXhrc+1/P/bruTYKHmT3Np7uE6SlOiO0uTMJL33edd9f748cZTWXD9do5LcdpcLAMCg6+/fb251haHrZ2YrLcmt6sYW7T1cp9hopyamJ6rW16L3Pj+mKKdDX5martXzz1T+xblyRzn1j08O6/Jfvq7Xdx82XT4AAMYw4nOCcBjxkaTSygZ9fKBKLqdD07NTNCI+Rhve+lwVdU1aeM5pyko5/rDTGt3++y3adbBGknTZ5HTdcekEnTnGE/LPP1jdqLQktxwObp8BAMzr799vgs8JwiX4nKzG5latemGnfvfOPlmW5HI69KOvTtbN5+ecdHj54/v79YOnPtQ540Zq7be+qNREbp8BAMziVhc6iY126SfXTNPGOy/SV8/MUKvf0k/+vEP3PPexTib7Nrf69eDGTyRJ735WoasfeVNHan2DVTYAAAOK4BNhctMS9cgNX9SKq6bI6ZB+/06xVjy/vd8bIj635YAOVDbImxij7JFx2n+sQb99/bNBrhoAgIHBra4TDNdbXd15+v39uvupDyVJTod0fq5Xd86ZqC+eNqJTu5ZWv37y5x16/sNSNTb71dDcqh9ePkmnj0rULRs2Kzk2SpuWXapEd5T2Ha3TJwdrgwesAgBgh/7+/WYfnwh23cwxcjmln2/creKKer2++4he331ESbFRSoiJUl1TizI9sRqZEKO391YEX5eaEKNvnzdW8dEujfcmaO+ROv3hvRJdf/YYfePXb+lgtU/rFp2tSycTfgAAQwsjPieIpBGf4+07WqeHX9mjP36wX/5uPhHRLod++vWzlJYUq3GjEjS6fdXY798p1o+e3abUhBidM26k/vpxuaS2fYM23HyOnb8CACCCDbtVXatWrdILL7ygrVu3KiYmRpWVlV3adLc6ae3atVqyZEm/f06kBp+AmsZmHarxqd7XqrgYp7YUV+qNPUf0jbOzdX6ut0v7xuZWXfs/m7SzrLrLc6/9W57GpibYUTYAIMINu+CzYsUKpaSkaP/+/Vq3bl2Pwefxxx/XvHnzgo95PB7FxcV1aduTSA8+oahubNatv3tfb+45qq9Nz1JVQ7Ne++Sw5n9xtBZ/ebzOSE9ivx8AwKAadnN8Vq5cKUkqLCzstV1KSooyMjL6/b4+n08+X8dy7OrqriMX6F1ybLQK/59z9NH+Sk0fk6KiXYf12ieH9cwHB/TMBwc0NStZt+adrnPHpaq4ok4f7a/S+FGJmjl2hBLdYfMRBAAMA2Ez4hNQWFiopUuX9jjiM3r0aDU2NmrcuHH67ne/q8WLF8vp7HnVfkFBQTBUHY8Rn9D5/ZZ++tI/9e5nFfpnWY0amrs/RT422qnvzMrRN84eo7GpCYp2sbsCACA0w+5WV0Bvwec//uM/dOmllyouLk5///vftXz5ci1btkz33ntvj+/X3YhPdnY2wWeAVNQ16bev79VL28v16eE6JcVGaebYEdpzqFb7jzUE27mcDo1MiNHolDjNPj1VF0zwaubYEXJHuQxWDwAIF2ERfHoabTnee++9p7PPPjv4fW/B50Rr1qzRfffdp6qqqn7XxByfwVPf1KIYl1NRLqcsy9Kruw7rV699qm0HqlTf1HVUyB3l1OmjEnV6WqLGexOUm5aoCemJGudNIBABADoJizk++fn5WrBgQa9tcnJyQn7/8847T9XV1Tp48KDS09lTxrT4mI6Pm8Ph0MWT0nTxpDT5/ZYO1fh0pNanf5bX6M09R/TGniM6XOPTjrJq7ThhxZjTIY0flaizxnj0hewUTR+TokmZSYQhAECfjAYfr9crr7frEumBsmXLFsXGxiolJWXQfgZOndPpUIYnVhmeWE0b7dF1M8fIsiztO1qvTw/Xtl2H6rTncK0+OVijmsYW7TlUqz2HavXMBwckSTEupyZnJesLYzyakJ6knNQEjU2NV1ZKnFxOVpQBANqEzZKa4uJiVVRUqLi4WK2trdq6daskKTc3V4mJifrTn/6k8vJyzZo1S3FxcSoqKtI999yjxYsXy+3m9PBw43A4lONNUI43odMO0JbVNjq0o7RaW0sq9eH+Sn1YUqlj9c36sKTt6+NFuxzKHhmvnNQE5aQmaEJ6oiamJ2pCepKSY6Nt/q0AAKaFzeTmm266SevXr+/yeFFRkfLy8vTiiy9q2bJl2rNnj/x+v8aPH6/vfe97+v73v6+oqP7nO+b4hB/LslRcUa+tJZXatr9Knx2p0+dH61RS0aCmVn+Pr/MmujXem6Acb7zGeRM1zhuv3LQk5aTGK4oVZgAQVsJicvNQRPAZPlr9lsqqGrTvaL0+P1qnvYfr9MnBGu0+WKvy6sYeXxfjcmr8qARNTE/SxPREnT4qUeNHJWpsarxio5lHBABDEcEnRASfyFDd2KzPj9Tps/br8yN12nukTnsO1Xa7wkxqm1Q9ZkS8JqQlKjetbbVZ4OskbpsBgFEEnxARfCKb32/pQGWDPjlYo08O1mr3wRp9eqROew/XqqaxpcfXZSTHakJ6WwjKTUvUhLQkTUhL1IiEGBurB4DIRfAJEcEH3bEsS4drfcHVZXsO1mjP4VrtPlirQzW+Hl+XmhAT3H8od1TbpOoJaYkaleTm/DIAGEAEnxARfHCyqhqa9enhWu05WKvdh2q051Ctdp+wM/WJkmOjdEZGUtuVnqQzMpJ1RnqSPPHcMgOAUBB8QkTwwUCpb2rR3sN12n2obUJ1YO+hfRX1avV3/59derK7PQQlBsPQhPREJlUDQB8IPiEi+GCw+VpagyvM/lleo0/K2/49UNn9CJHTIY1NTdAZ6UmaGBwhYtk9AByP4BMigg9MqWls1icH23an3lXefh2sUUVdU7ftY6Kcyh2V2OmW2cSMJGV5Ypk/BCDiEHxCRPDBUGJZlo7UNgVD0K7yau1qX23W07L7pNgoTc5M1pTAlZWsCemJnGUGYFgj+ISI4INw4Pdb2n+soVMY2lVerb2H69TSzfyhKKdDuWmJHYEoK1mTM5M1kuX2AIYJgk+ICD4IZ00tfu05VKudZdXa2X6y/Y6yalXWN3fbPiM5VlOyOsLQlMxknTYyXk4OdgUQZgg+ISL4YLixLEtlVY1tQai0IwztO1rfbfuEGJcmZXYOQ2dkJLGyDMCQRvAJEcEHkaKmsVm7ymvaglB7INpVXiNfS9eDXV1Oh3JHJWpqVlsYmprl0ZSsZHni2HcIwNBA8AkRwQeRrKXVr8+O1HUKQ9tLq3tcWZY9Mk5TMz2ampWsaaPb/k1LjrW5agAg+ISM4AN0ZlmWDlb7tL20SttLq/XxgbZ/e9p3yJvo1tSs5PbLo2mjk5U9gnlDAAYXwSdEBB+gfyrrm7SjtG1EKBCKPj1cq+42pU5yR2nycWFoalayctMSFc0GjAAGCMEnRAQfIHQNTa3aWd4Whna0h6F/lteoqZt5QzFRTp2RnqRpo5M1pT0MTc5IVlwMk6gBnDyCT4gIPsDAam7169PDtfr4QMfI0M7SatX4Wrq0dTqk8e2TqI8fHUqJZ78hAL0j+ISI4AMMPr/fUsmx+k63yT4+UK0jtb5u249OiesUhKaOTlZGMkdzAOhA8AkRwQcw51B1Y6cwtL20WsUV3e83lJoQE1xaHxghyklNYBI1EKEIPiEi+ABDS1VDs3aWta0mC0ym3nO4Vq3dzKJOiHFpcmb7bbL25fUT0pIUE8UkamC4I/iEiOADDH2Nza3aVV7TaXRoZ1l1t5svRrscmpie1OlW2eTMZCW4owxUDmCwEHxCRPABwlNLq197j9S1BaEDHcvsqxu7TqJ2OKRx3oROt8mmZnk4tBUIYwSfEBF8gOHDstpOse88b6hKB6u7n0Sd6YntHIZGe5TlYRI1EA4IPiEi+ADD35FaX+cwdKBKn/dwaGtKfHTnFWVZHo3zJsjFJGpgSCH4hIjgA0SmmsZm7Syr6bSibPfBGrV0M4k6LtqlyZlJncLQxIxEuaPYfBEwheATIoIPgABfS6t2H6wNnk+2vbRKO8tq1NDc2qVtlNOh3LTE4GGtU7M8mpyZpKRYTrAH7EDwCRHBB0BvWv2WPmufRB1YXv9xaZUq65u7bT9mRJzGj0rU6aMSdPqoxPYrQaOS3MwdAgYQwSdEBB8AJ8uyLJVWNWr7gY7bZDtKq1Ra1djja5LcURo9Ik5jRsRpdEqcxoyI1+j2r0ePiFNqQgzBCDgJBJ8QEXwADJRjdU3ac7hWnx6q1aeHa7X3cJ0+PVyr4or6bk+xP15stFMZybEaleRWWlLbv4ErLfhvrEYmxDDRGhDBJ2QEHwCDzdfSqpKKepUca9CBYw3af6xBByobdOBYvfYfa9Chmu6X23fH6ZBS4mPkiYtWcly0UuKi5Wm/UuKjg4974qKVFBulRHeUEtwd/8ZHuzjmA8NCf/9+s3UpANjMHeVSblqSctOSun3e19KqsspGHaxu1OFanw5V+074t1FHan06WtckvyVV1DWpoq4ppFocDikhJkoJbpcS3FFKag9EHeGo7fGEmCjFRbsUF+NSfPsV1/5Y29ftj0dHKTbGqRiXk1t1GJIIPgAwxLijXMrxJijHm9Bru5ZWv47WNamyvllVDc2qrG9SVUPb19UNzaps/zpw1Ta2qM7Xotr2y29JlqXg91L/R5r64nI6FH9cUIqNDgSmKMXFuLoGppio49oEno8KPp9w3NdxjFLhFBB8ACBMRbmcSk+OVXpy7Em/1rIsNTb7VevrHIbqgv+2dnq8oblVDU2tqm9qUX1T29cdj7U93tDcqubWttkTrX5LNb4W1fi6HhkyEI4PTseHoviYqGB4Cnzd1qZzkOrUrv2WX1yMS+4oRqqGO4IPAEQgh8PRNvIS49KoJPeAvW9zq79LIGpsDnzdEZjavm4JPt6pTXNHuOp4XYvqm1sVmJXa0Nz2PqobsNIldR6pSnC33cpLbJ8bleiOUmJs2+3A4FypwPcntElsvz3IyNTQQ/ABAAyYaJdT0S6nkgdh48bAKFX9cYGpx699bUGpoalt5Or4rwPBK9CmvqlVTS1+SSeMVJ3EJPOeBOZJtQWi6PY5VC4luk+cbO4Kzq3qdp4VIWrAEHxs1Nfw6fEL7GgbPm2Pb0/b4d+2P+1pO/htUwfhfasamoLh6PQeJp4H/PCPH6qmse024Pqbz+21rSSN/X//LEna99MrB7TtzY+/qwR3lP77hi/22fav28oU7XLosikZfbb9+EClXE6HJmd6+mwb6n9HphB8EDEcDke//6M7mbaILCcz/4O5IuHFExfT7//u7//69GDb9Tf33f5/bzlPtb4Wzf1p320npCWqzteiff2o4+//PNSPVm2W/O79fre94qE3+t023DhNF9Bfq1at0uzZsxUfH6+UlJQe2xUWFuqss85SbGysMjIylJ+fb1+RAAB0Y9bpqZozJb1fbTfedZE2Lbu0X23vn3+m7r1icr/azhw7QtPH9D2CI0mjktxKTYjpV9twEzYbGK5YsUIpKSnav3+/1q1bp8rKyi5tHnzwQa1Zs0YPPPCAzj33XDU2Nmrv3r266qqr+v1zBnMDw3AZQqYtt7po23Pb/rSnLW17ak/bwYscw3bn5sLCQi1durRL8Dl27JhGjx6tP/3pT7r00v4l5e6wczMAAOGnv3+/w+ZWV182btwov9+vAwcOaPLkyRozZoy+8Y1vqKSkpNfX+Xw+VVdXd7oAAMDwNGyCz969e+X3+/Wf//mf+sUvfqGnn35aFRUVmjNnjpqaet7KffXq1fJ4PMErOzvbxqoBAICdjAafgoICORyOXq/Nmzf36738fr+am5v10EMP6Stf+YrOO+88/e///q92796toqKiHl+3bNkyVVVVBa++RogAAED4MrqcPT8/XwsWLOi1TU5OTr/eKzMzU5I0ZcqU4GOjRo2S1+tVcXFxj69zu91yuwdu11IAADB0GQ0+Xq9XXq93QN7r/PPPlyTt2rVLY8aMkSRVVFToyJEjGjt27ID8DAAAEN7CZgPD4uJiVVRUqLi4WK2trdq6daskKTc3V4mJiZo4caKuvvpq3XHHHXr00UeVnJysZcuWadKkSbr44ovNFg8AAIaEsAk+y5cv1/r164Pfz5gxQ5JUVFSkvLw8SdKGDRt055136oorrpDT6dRFF12kF198UdHRA39mDAAACD9ht4/PYGMfHwAAwk/E7eMDAADQF4IPAACIGAQfAAAQMQg+AAAgYhB8AABAxAib5ex2CSxy47BSAADCR+Dvdl+L1Qk+J6ipqZEkDisFACAM1dTUyOPx9Pg8+/icwO/3q7S0VElJSXI4HAP2vtXV1crOzlZJSQn7A50i+nLg0JcDh74cOPTlwImkvrQsSzU1NcrKypLT2fNMHkZ8TuB0OoNnfQ2G5OTkYf/hswt9OXDoy4FDXw4c+nLgREpf9jbSE8DkZgAAEDEIPgAAIGIQfGzidru1YsUKud1u06WEPfpy4NCXA4e+HDj05cChL7ticjMAAIgYjPgAAICIQfABAAARg+ADAAAiBsEHAABEDIKPTf7nf/5H48aNU2xsrGbOnKnXX3/ddElDWkFBgRwOR6crIyMj+LxlWSooKFBWVpbi4uKUl5en7du3G6x46PjHP/6hq666SllZWXI4HHruuec6Pd+fvvP5fLr99tvl9XqVkJCgr33ta9q/f7+Nv8XQ0Fdf3nTTTV0+p+edd16nNvSltHr1an3pS19SUlKS0tLSdM0112jXrl2d2vC57J/+9CWfy94RfGzwf//3f1q6dKnuuecebdmyRRdeeKEuv/xyFRcXmy5tSJs6darKysqC17Zt24LP/exnP9ODDz6ohx9+WO+9954yMjI0Z86c4Flrkayurk7Tp0/Xww8/3O3z/em7pUuX6tlnn9WTTz6pN954Q7W1tbryyivV2tpq168xJPTVl5I0b968Tp/Tv/zlL52epy+l1157Td///vf19ttva+PGjWppadHcuXNVV1cXbMPnsn/605cSn8teWRh055xzjrVkyZJOj02aNMn64Q9/aKiioW/FihXW9OnTu33O7/dbGRkZ1v333x98rLGx0fJ4PNavfvUrmyoMD5KsZ599Nvh9f/qusrLSio6Otp588slgmwMHDlhOp9N68cUXbat9qDmxLy3LshYtWmRdffXVPb6GvuzeoUOHLEnWa6+9ZlkWn8tTcWJfWhafy74w4jPImpqa9P7772vu3LmdHp87d642bdpkqKrwsHv3bmVlZWncuHFasGCB9u7dK0n67LPPVF5e3qlP3W63LrroIvq0D/3pu/fff1/Nzc2d2mRlZWnatGn0bzdeffVVpaWlaeLEibrlllt06NCh4HP0ZfeqqqokSSNHjpTE5/JUnNiXAXwue0bwGWRHjhxRa2ur0tPTOz2enp6u8vJyQ1UNfeeee642bNigl156Sb/5zW9UXl6u2bNn6+jRo8F+o09PXn/6rry8XDExMRoxYkSPbdDm8ssv1xNPPKFXXnlFa9as0XvvvadLLrlEPp9PEn3ZHcuydNddd+mCCy7QtGnTJPG5DFV3fSnxuewLp7PbxOFwdPresqwuj6HD5ZdfHvz6zDPP1KxZs3T66adr/fr1wUl69GnoQuk7+rerb37zm8Gvp02bprPPPltjx47VCy+8oPnz5/f4ukjuy/z8fH300Ud64403ujzH5/Lk9NSXfC57x4jPIPN6vXK5XF1S9KFDh7r83w16lpCQoDPPPFO7d+8Oru6iT09ef/ouIyNDTU1NOnbsWI9t0L3MzEyNHTtWu3fvlkRfnuj222/X888/r6KiIo0ZMyb4OJ/Lk9dTX3aHz2VnBJ9BFhMTo5kzZ2rjxo2dHt+4caNmz55tqKrw4/P5tHPnTmVmZmrcuHHKyMjo1KdNTU167bXX6NM+9KfvZs6cqejo6E5tysrK9PHHH9O/fTh69KhKSkqUmZkpib4MsCxL+fn5euaZZ/TKK69o3LhxnZ7nc9l/ffVld/hcnsDMnOrI8uSTT1rR0dHWunXrrB07dlhLly61EhISrM8//9x0aUPWD37wA+vVV1+19u7da7399tvWlVdeaSUlJQX77P7777c8Ho/1zDPPWNu2bbMWLlxoZWZmWtXV1YYrN6+mpsbasmWLtWXLFkuS9eCDD1pbtmyx9u3bZ1lW//puyZIl1pgxY6yXX37Z+uCDD6xLLrnEmj59utXS0mLq1zKit76sqamxfvCDH1ibNm2yPvvsM6uoqMiaNWuWNXr0aPryBLfeeqvl8XisV1991SorKwte9fX1wTZ8Lvunr77kc9k3go9NHnnkEWvs2LFWTEyM9cUvfrHT0kN09c1vftPKzMy0oqOjraysLGv+/PnW9u3bg8/7/X5rxYoVVkZGhuV2u60vf/nL1rZt2wxWPHQUFRVZkrpcixYtsiyrf33X0NBg5efnWyNHjrTi4uKsK6+80iouLjbw25jVW1/W19dbc+fOtUaNGmVFR0dbp512mrVo0aIu/URfWt32oSTr8ccfD7bhc9k/ffUln8u+OSzLsuwbXwIAADCHOT4AACBiEHwAAEDEIPgAAICIQfABAAARg+ADAAAiBsEHAABEDIIPAACIGAQfAAAQMQg+ANAHh8Oh5557znQZAAYAwQfAkHbTTTfJ4XB0uebNm2e6NABhKMp0AQDQl3nz5unxxx/v9Jjb7TZUDYBwxogPgCHP7XYrIyOj0zVixAhJbbeh1q5dq8svv1xxcXEaN26cnnrqqU6v37Ztmy655BLFxcUpNTVVixcvVm1tbac2jz32mKZOnSq3263MzEzl5+d3ev7IkSO69tprFR8frwkTJuj5558f3F8awKAg+AAIez/+8Y/19a9/XR9++KG+/e1va+HChdq5c6ckqb6+XvPmzdOIESP03nvv6amnntLLL7/cKdisXbtW3//+97V48WJt27ZNzz//vHJzczv9jJUrV+ob3/iGPvroI331q1/Vt771LVVUVNj6ewIYAKaPhweA3ixatMhyuVxWQkJCp+u+++6zLMuyJFlLlizp9Jpzzz3XuvXWWy3LsqxHH33UGjFihFVbWxt8/oUXXrCcTqdVXl5uWZZlZWVlWffcc0+PNUiy7r333uD3tbW1lsPhsP76178O2O8JwB7M8QEw5F188cVau3Ztp8dGjhwZ/HrWrFmdnps1a5a2bt0qSdq5c6emT5+uhISE4PPnn3++/H6/du3aJYfDodLSUl166aW91nDWWWcFv05ISFBSUpIOHToU6q8EwBCCD4AhLyEhocutp744HA5JkmVZwa+7axMXF9ev94uOju7yWr/ff1I1ATCPOT4Awt7bb7/d5ftJkyZJkqZMmaKtW7eqrq4u+Pybb74pp9OpiRMnKikpSTk5Ofr73/9ua80AzGDEB8CQ5/P5VF5e3umxqKgoeb1eSdJTTz2ls88+WxdccIGeeOIJvfvuu1q3bp0k6Vvf+pZWrFihRYsWqaCgQIcPH9btt9+uG2+8Uenp6ZKkgoICLVmyRGlpabr88stVU1OjN998U7fffru9vyiAQUfwATDkvfjii8rMzOz02BlnnKF//vOfktpWXD355JO67bbblJGRoSeeeEJTpkyRJMXHx+ull17SHXfcoS996UuKj4/X17/+dT344IPB91q0aJEaGxv185//XHfffbe8Xq+uu+46+35BALZxWJZlmS4CAELlcDj07LPP6pprrjFdCoAwwBwfAAAQMQg+AAAgYjDHB0BY4249gJPBiA8AAIgYBB8AABAxCD4AACBiEHwAAEDEIPgAAICIQfABAAARg+ADAAAiBsEHAABEjP8fSJigJTCRZhMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Losses\n",
    "\n",
    "plt.plot(tc.stack(losses))\n",
    "plt.hlines(tc.full((len(losses),), func.min_likelihood), xmin=0, xmax=len(losses)+10, linestyles='--', color='black', alpha=0.4)\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"Loss\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, '$\\\\beta (t)$')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGzCAYAAAD9pBdvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXZ0lEQVR4nO3deXhTZdoG8Dtbk64ptHRfKGuBspZ9ERCt4j46ioCgDszA4AIyDsowijKO+Lkw6CgoijCOoujIjKKoFGXfCoVC2ZdSWkpLN7q3SZOc74/kpK0t0CVtzjm9f9eV69Lk5ORND22ePO/zPq9KEAQBRERERG6idvcAiIiIqH1jMEJERERuxWCEiIiI3IrBCBEREbkVgxEiIiJyKwYjRERE5FYMRoiIiMitGIwQERGRWzEYISIiIrdiMEJERERuxWCEiIiI3IrBCBG5XEVFBV5++WXExsbCYDAgICAAkydPxtWrV+sct2TJEvTu3Rs2m63O/YIgYMmSJdi+fbvzvtWrVyM8PBzl5eVt8h6IqO0wGCEilxIEAZMnT8by5csxa9Ys/PDDD3jxxRexfv16vPXWW87jLl++jNdffx1LliyBWl33T9GZM2ewePFiZGdnO+979NFH4e3tjddff73N3gsRtQ0GI0TkUtu3b8e3336LDz/8EM888wzGjx+PuXPnIigoqE5W4+2334a/vz/uv//+eudITk4GAMTHxzvv02q1mDVrFt5++21UVFS0/hshojbDYISIXOqrr75Chw4d8Jvf/MZ5344dO3DlyhXcfPPNAACz2YzVq1djypQp9bIi8fHxmDp1KgCgR48eUKlU8PX1hSAImDp1KkpKSvDFF1+03RsiolbHYISIXGrPnj0YNmwYAPtUzJo1a/DAAw9gwoQJuOuuuwAA+/fvR0FBAcaPH1/v+atWrcLgwYMxaNAg7N27F3v37sWuXbugUqkQEhKC2NhYfP/99236noiodWndPQAiUo6qqiocO3YM9957L1599VW8+OKLAICoqCisXbsWKpUKALB3714AwKBBg+qdIz4+HhkZGXj44YcxfPjweo8PGjQIW7ZsacV3QURtjZkRInKZQ4cOwWKxYOjQoZg6dSp++uknvPzyyygtLcVNN92EsrIyAPaMiUqlQmBgYL1zZGZmIjc3t069SG1BQUHIzc2FxWJp1fdCRG2HmREicpmkpCQAwNChQxEYGIguXbogISEBPXr0wOTJk7Fv3z7ccsstqKyshE6ng0ajqXeOhopXazMYDBAEAVVVVfDx8Wm9N0NEbYaZESJymaSkJHTp0qXBjAcAhIWFAQACAwNhNpsb7BmSnJwMLy8vxMbGNniOwsJC6PV6BiJECsJghIhcJikpCaGhoXXuEwQBH330EeLi4tC7d28AcAYa58+fr3eOo0ePIjY2tsGsCQCkpaU5z0NEysBghIhcorCwEOfPn8fu3bsxc+ZM/Pzzz9i4cSPuv/9+7Nq1C6tWrXIeO27cOADAvn376p3H398f58+fxzfffIN9+/bh4sWLzsdsNhuSkpIaXIVDRPLFYISIXEKsF5k+fTp+/PFHTJw4EU8++SS0Wi3279+PESNGOI+NjIzEmDFj8M0339Q7z+LFizFgwABMnToVI0aMwNdff+18bNu2bSguLnb2ISEiZVAJgiC4exBEJH9LlizBK6+8gpKSEhgMhhse//XXX2PSpEm4ePEiwsPDG/Ua06ZNQ1paGnbv3t3S4RKRhDAYISKXuOuuu5Cdne1cDXMjgiBg5MiRiI+Px7vvvnvD48+fP49evXrhl19+wejRo1s6XCKSEE7TEJFLHDhwAEOHDm308SqVCh9++CHCwsLq7drbkIyMDLz77rsMRIgUiJkRIiIicitmRoiIiMitGIwQERGRWzEYISIiIreSxd40NpsNly9fhq+vr3PXTyIiIpI2QRBQWlqKsLAwqNXXzn/IIhi5fPkyIiMj3T0MIiIiaobMzExERERc83FZBCO+vr4A7G/Gz8/PzaMhIiKixigpKUFkZKTzc/xaZBGMiFMzfn5+DEaIiIhk5kYlFixgJSIiIrdiMEJERERuxWCEiIiI3IrBCBEREbkVgxEiIiJyKwYjRERE5FYMRoiIiMitGIwQERGRWzEYISIiIrdiMEJERERuxWCEiIiI3IrBCBEREblVuw5G9qcVYPrHScgtqXL3UIiIiNwi8cQVfLQzDUUVZreNod0GI4Ig4I2fTmPHmTys2Hbe3cMhIiJyi3e3nsMr35/EZ/sz3DaGdhuMqFQqzL+1BwDg86QMmCxWN4+IiIiobZ24XIIjmUXw0KgxaUik28bRboMRABjRNQCBPh4wWWw4eqnY3cMhIiJqU3vTCgAAY7oHItBH77ZxtOtgRKVSYUjnjgCApAuFbh4NERFR2zqUcRUAMCi6g1vH0a6DEQDoG2EEAJy9UurmkRAREbWtY1n2WYH+Ef5uHUe7D0a6dvIBAJzPK3fzSIiIiNqO2WJDZmEFAKB7sI9bx8JgxBmMlEEQBDePhoiIqG1kFFbAJgDeHhoE+bqvXgRgMILoAC+oVUCF2YrcUpO7h0NERNQmLuTbZwRiOnlDpVK5dSztPhjRadQI9jMAAC4XVbp5NERERG0j66p9iiayg5ebR8JgBAAQ5u8JAMguZidWIiJqH644ZgPEL+TuxGAENcEIMyNERNReXHFshcJgRCLCjPYLkcVghIiI2oncEjEz4t7iVaCZwciKFSsQExMDg8GA+Ph47Ny585rHbtu2DSqVqt7t1KlTzR60q4lRYR4LWImIqJ2QdWZk/fr1mDdvHhYtWoTDhw9jzJgxmDhxIjIyrr/BzunTp5Gdne28de/evdmDdrUAHw8AQEGZ+3YsJCIiaks1wYgMMyPLli3DjBkzMHPmTPTq1QvLly9HZGQkVq5ced3nBQUFISQkxHnTaDTXPNZkMqGkpKTOrTUFeNsvREE5MyNERKR8lWYrSqosAIAguWVGzGYzkpOTkZCQUOf+hIQE7Nmz57rPHThwIEJDQzFhwgRs3br1uscuXboURqPReYuMbN2dBDt62zMjheXMjBARkfLlltqzIp46DXz1WjePponBSH5+PqxWK4KDg+vcHxwcjJycnAafExoailWrVuHrr7/Ghg0b0LNnT0yYMAE7duy45ussXLgQxcXFzltmZmZThtlkgT41wYjNxi6sRESkbDnFNVM07m54BgDNCod+PXBBEK75Znr27ImePXs6/3/EiBHIzMzEm2++iZtuuqnB5+j1euj1bTeH1cGRGbEJQFFltTNTQkREpERijxEpTNEATcyMBAYGQqPR1MuC5Obm1suWXM/w4cNx9uzZprx0q9Jp1DB66gAABWWsGyEiImXLdRSvuntPGlGTghEPDw/Ex8cjMTGxzv2JiYkYOXJko89z+PBhhIaGNuWlW51zRQ3rRoiISOHEGslAH2kEI02eppk/fz6mTZuGwYMHY8SIEVi1ahUyMjIwe/ZsAPZ6j6ysLHzyyScAgOXLl6Nz587o06cPzGYzPv30U3z99df4+uuvXftOWijA2wNpeeVc3ktERIpXVFkNAM5ZAXdrcjAyadIkFBQUYMmSJcjOzkZcXBw2bdqE6OhoAEB2dnadniNmsxnPPvsssrKy4OnpiT59+uD777/HHXfc4bp34QLi8t5CLu8lIiKFK6qwf/Hu4CXTYAQA5syZgzlz5jT42Nq1a+v8/4IFC7BgwYLmvEyb6uiYpslnZoSIiBSuqMKeGfH3ksaCDe5N4xDoLdaMMDNCRETKdtUZjEgjM8JgxEFczsuaESIiUrpixzQNMyMSI/YaKXYU9RARESmVmBmRSs0IgxEHP4P9gpQ6evUTEREpUVW1FZXVVgCAvyczI5Li52mv5S2pYmaEiIiUS5wBUKsAX4P796UBGIw4iZmREk7TEBGRgokraYyeOqjV7t+XBmAw4uQrBiNVFggCN8sjIiJluursMSKNKRqAwYiTOE1jtQmoMFvdPBoiIqLW4cyMSKR4FWAw4uSp00DrSFexboSIiJSqiJkR6VKpVPDzFOtGuKKGiIiUydnwTCL70gAMRurwM3BFDRERKVtRpT0zwmkaiarJjDAYISIiZRKz/1LZsRdgMFIHG58REZHSlZnsn3HiKlIpYDBSCxufERGR0pU6PuN89dJoeAYwGKnDV89pGiIiUrayKjEzwmBEkmoyI5ymISIiZSqt4jSNpLElPBERKZ04TePDzIg0OVfTsGaEiIgUqtTEaRpJc07TsOkZEREpkM0m1FpNw2BEkpwFrMyMEBGRAlVUWyHuBSt+5kkBg5FavB3LnMpNzIwQEZHyiPUiWrUKBp10QgDpjEQCfJzBCHftJSIi5SmttaxXpVK5eTQ1GIzU4q3XAGBmhIiIlEkMRqS0kgZgMFKHc5rGbIEgTqoREREpRE33VenUiwAMRuoQgxGbAFRV29w8GiIiIteS4koagMFIHV46jfO/yzhVQ0RECiPF7qsAg5E61GoVvD1YN0JERMrknKZhZkTaateNEBERKYkUN8kDGIzU483lvUREpFAlDEbkgct7iYhIqcR6SB+uppE2bw9O0xARkTKJX7TZZ0TifNgSnoiIFKrcbC9BEBdrSAWDkV/xcgQjZawZISIihalwfNH28mBmRNJ8WDNCREQK5cyM6JkZkTTWjBARkVJVmJkZkQVv1owQEZFCiW0rmBmROB/2GSEiIoUSMyPezIxIm5cjWuTeNEREpCQ2m4AKR82IF1fTSJuYGalgzQgRESlIZXVNxl8sSZAKBiO/IqauuLSXiIiURFyYoVYBeq20Pv6lNRoJYAErEREpUYVYvOqhhUqlcvNo6mIw8itihXEFgxEiIlIQMTPiJbGVNACDkXrEop7ac2tERERyV2GuyYxIDYORXzHoHJkRM4MRIiJSDrH8gJkRGRC70pksNlhtgptHQ0RE5Bo1y3qZGZG82muvqzhVQ0RECiFmRqS2Yy/AYKQevVYNsciYUzVERKQUzsyIxHqMAAxG6lGpVPB01I1UMhghIiKFKDczMyIrzmCE0zRERKQQYp8R1ozIhKeHuKKGvUaIiEgZxD3XpLZjL8BgpEHOXiOcpiEiIoUQv2AzMyITnKYhIiKlKXc2PWNmRBZqpmkYjBARkTJUOJueKSQzsmLFCsTExMBgMCA+Ph47d+5s1PN2794NrVaLAQMGNOdl24yYwmJmhIiIlKJcSe3g169fj3nz5mHRokU4fPgwxowZg4kTJyIjI+O6zysuLsb06dMxYcKEZg+2rXBpLxERKU2FkjbKW7ZsGWbMmIGZM2eiV69eWL58OSIjI7Fy5crrPm/WrFmYMmUKRowY0ezBthVO0xARkdKIS3tlnxkxm81ITk5GQkJCnfsTEhKwZ8+eaz5vzZo1OH/+PBYvXtyo1zGZTCgpKalza0ssYCUiIqUpd66mkXlmJD8/H1arFcHBwXXuDw4ORk5OToPPOXv2LJ5//nl89tln0GobF40tXboURqPReYuMjGzKMFusZmkv+4wQEZEyODMjSilgVYmbtzgIglDvPgCwWq2YMmUKXn75ZfTo0aPR51+4cCGKi4udt8zMzOYMs9k4TUNEREoiCIKk28E3KTwKDAyERqOplwXJzc2tly0BgNLSUhw8eBCHDx/Gk08+CQCw2WwQBAFarRabN2/GzTffXO95er0eer2+KUNzKU7TEBGRkpgsNtgE+3/Lfmmvh4cH4uPjkZiYWOf+xMREjBw5st7xfn5+SE1NRUpKivM2e/Zs9OzZEykpKRg2bFjLRt9K2IGViIiUpNxUU3YgfuGWkiaHR/Pnz8e0adMwePBgjBgxAqtWrUJGRgZmz54NwD7FkpWVhU8++QRqtRpxcXF1nh8UFASDwVDvfinxdFQac5qGiIiUQPw889RpoFHXL6twtyYHI5MmTUJBQQGWLFmC7OxsxMXFYdOmTYiOjgYAZGdn37DniNRxmoaIiJTEWS8iwR4jAKASBEFw9yBupKSkBEajEcXFxfDz82v119t6KhePrz2AvuFGbHxqdKu/HhERUWtKvngVD6zcg6iOXtixYHybvW5jP7+5N00DDDpxNQ2X9hIRkfxVSLjHCMBgpEHixaqqtrl5JERERC3nrBlhMCIfXh7MjBARkXJUOWogmRmRkZppGhawEhGR/NVeTSNFDEYaIEaOJosNNpvk63uJiIiuq9I5TSO9hmcAg5EGedW6WFzeS0REcid+lnnqpPmxL81RuZleW/NjYTBCRERyV8lpGvlRq1U1jc9YN0JERDLnzIxwmkZeuHMvEREpBQtYZYot4YmISCm4tFem2GuEiIiUQvwsMzAYkRdxmoY1I0REJHeVjo7iXpymkRdO0xARkVJUOjIjbAcvMyxgJSIipahZTcNgRFa8OE1DREQKwT4jMuWps6/F5jQNERHJHYMRmfL0sP9oOE1DRERyV8mlvfIkRo8mZkaIiEjmxC/WBmZG5IWraYiISAlsNgEmi2NpLzMj8mJgASsRESlA7S/VXE0jM8yMEBGREtT+HDNoGYzIihiMVDEYISIiGat01ouooVar3DyahjEYuQaDMxixuXkkREREzVezkkbr5pFcG4ORazBwmoaIiBRA6j1GAAYj18SN8oiISAkqak3TSJV0R+ZmrBkhIiIlqOI0jXyJESSnaYiISM4qOE0jX8yMEBGREkh9x16Awcg1sYCViIiUoNJsAcDMiCyJEWRVtQ02m+Dm0RARETWP1DfJAxiMXFPtCFLs6U9ERCQ3ztU0DEbkp/bOhqwbISIiuXJmRjhNIz8atQoeGq6oISIieasys4BV1ri8l4iI5K6m6RmDEVliF1YiIpI7FrDKHHuNEBGR3HFvGpnjzr1ERCR3bHomc2x8RkREcsd28DLnyWCEiIhkjhvlyZyzCysLWImISKZqpmmk+5Ev3ZFJAJf2EhGR3NVM0zAzIksGrqYhIiKZY9MzmWPNCBERyZkgCKioZgGrrDEYISIiOau2CrA6dp5nZkSmWMBKRERyVruDODMjMsWmZ0REJGdiZl+rVsFDK92PfOmOTALY9IyIiOSswmwBIO2sCMBg5LpYM0JERHImh1bwAIOR6xIbxHBpLxERyVEVgxH5c2ZGWMBKREQyJId9aQAGI9elFwtYLQxGiIhIfipl0PAMYDByXcyMEBGRnFXKoOEZwGDkujy5tJeIiGRM/DLtxcyIfIlpLa6mISIiORJrRgxKzIysWLECMTExMBgMiI+Px86dO6957K5duzBq1CgEBATA09MTsbGx+Mc//tHsAbclg5bTNEREJF/il2mpZ0aavJ/w+vXrMW/ePKxYsQKjRo3CBx98gIkTJ+LEiROIioqqd7y3tzeefPJJ9OvXD97e3ti1axdmzZoFb29v/OEPf3DJm2gtBnFpr8UKQRCgUqncPCIiIqLGq1Tqapply5ZhxowZmDlzJnr16oXly5cjMjISK1eubPD4gQMHYvLkyejTpw86d+6MRx55BLfddtt1sylSIV48QQBMFtaNEBGRvNQ0PWty7qFNNSkYMZvNSE5ORkJCQp37ExISsGfPnkad4/Dhw9izZw/Gjh17zWNMJhNKSkrq3Nyh9hwbG58REZHcKLLPSH5+PqxWK4KDg+vcHxwcjJycnOs+NyIiAnq9HoMHD8YTTzyBmTNnXvPYpUuXwmg0Om+RkZFNGabL6DRq6DT2qRkWsRIRkdxUyaRmpFkFrL+unWhMPcXOnTtx8OBBvP/++1i+fDk+//zzax67cOFCFBcXO2+ZmZnNGaZLiEWsXN5LRERyI9aMGCQejDRpEikwMBAajaZeFiQ3N7detuTXYmJiAAB9+/bFlStX8NJLL2Hy5MkNHqvX66HX65sytFZj8NCg1GThihoiIpKdCjEzoqRpGg8PD8THxyMxMbHO/YmJiRg5cmSjzyMIAkwmU1Ne2m24cy8REclVlUzawTe5vHb+/PmYNm0aBg8ejBEjRmDVqlXIyMjA7NmzAdinWLKysvDJJ58AAN577z1ERUUhNjYWgL3vyJtvvomnnnrKhW+j9dR0YWUwQkRE8lJRbQEg/QLWJgcjkyZNQkFBAZYsWYLs7GzExcVh06ZNiI6OBgBkZ2cjIyPDebzNZsPChQtx4cIFaLVadO3aFa+99hpmzZrlunfRisR5Nk7TEBGR3Mhlo7xmLTyeM2cO5syZ0+Bja9eurfP/Tz31lGyyIA0xaGsanxEREcmJYpuetTeezIwQEZFMyaUdPIORG2DNCBERyZWiN8prT7iahoiI5MhmE5xbmTAzInN6MRgxs+kZERHJR+0v0VIvYGUwcgPOaRoWsBIRkYzUDkbEbuJSxWDkBjw97D8iFrASEZGcOFvB69RQq6+/ZYu7MRi5ARawEhGRHNWspGlWF482xWDkBgwsYCUiIhmSS48RgMHIDRmYGSEiIhmqkEn3VYDByA3VLO3lahoiIpIP8Us0MyMKIEaUVSxgJSIiGangNI1ysOkZERHJkfi5xWkaBWABKxERyVGl2QKAmRFFMOgcu/YyGCEiIhmRyyZ5AIORG3LWjDAYISIiGXFuksdgRP6cNSMsYCUiIhlxZkY4TSN/tQtYBUFw82iIiIgap5J9RpRD3LXXJgBmK3uNEBGRPDAYUZDaVchVbHxGREQyUcmmZ8qh06igcex2yCJWIiKSCzEzwtU0CqBSqVjESkREsiNmRgzMjCgDG58REZHcsB28wrDxGRERyU2Vs+mZ1s0juTEGI43A/WmIiEhunJkRD+l/1Et/hBLALqxERCQ3NatpmBlRBGfNiJlLe4mISB7YZ0RhOE1DRERyIggCN8pTGhawEhGRnJitNlht9i1MuLRXIcTMCIMRIiKSg6paZQXMjCiEON/GpmdERCQHFdUWAIBWrYJOI/2PeumPUALY9IyIiORETsWrAIORRmEwQkREciKnTfIABiONUlMzwqW9REQkfcyMKBALWImISE6YGVEgAwtYiYhIRiqYGVEeNj0jIiI5qZJRwzOAwUijsOkZERHJiTMzwmka5WDNCBERyUlNAav0N8kDGIw0CqdpiIhITmoKWOXxMS+PUbqZs4CVwQgREcmAmBnxYmZEOZyZETP7jBARkfSJNSNy2CQPYDDSKOLFNDEzQkREMlDJ1TTKw5oRIiKSk0qzfaM8rqZREPFiWmwCqq2cqiEiImkrd0zTeOtZM6IYBo+aHxOzI0REJHUVjsyIt56ZEcXw0KihUtn/u4ot4YmISOLKTWx6pjgqlYo79xIRkWzUZEY4TaMoLGIlIiK5qDBzNY0iGRiMEBGRTFSwgFWZxG2YK1kzQkREEldusk/TMDOiMM6dey0MRoiISLosVhtMFnt9ozfbwSuLs4CVmREiIpKwilrlBF5c2qssrBkhIiI5qHAs69WqVfDQyONjXh6jlACupiEiIjkoN9fUi6jEJlkS16xgZMWKFYiJiYHBYEB8fDx27tx5zWM3bNiAW2+9FZ06dYKfnx9GjBiBn376qdkDdhcWsBIRkRyImRG5rKQBmhGMrF+/HvPmzcOiRYtw+PBhjBkzBhMnTkRGRkaDx+/YsQO33norNm3ahOTkZIwfPx533303Dh8+3OLBtyWD1rFzr4VNz4iISLrEzIinTFbSAM0IRpYtW4YZM2Zg5syZ6NWrF5YvX47IyEisXLmyweOXL1+OBQsWYMiQIejevTteffVVdO/eHRs3bmzx4NsSMyNERCQHzu6rMllJAzQxGDGbzUhOTkZCQkKd+xMSErBnz55GncNms6G0tBQdO3a85jEmkwklJSV1bu7GAlYiIpIDcV8aufQYAZoYjOTn58NqtSI4OLjO/cHBwcjJyWnUOd566y2Ul5fjoYceuuYxS5cuhdFodN4iIyObMsxWwQJWIiKSg0qZdV8FmlnA+uvqXEEQGlWx+/nnn+Oll17C+vXrERQUdM3jFi5ciOLiYuctMzOzOcN0KWfTM4VN0wiCgA2HLuHJdYfwRVIGBEFw95CIiKgFaq+mkYsmhU2BgYHQaDT1siC5ubn1siW/tn79esyYMQNfffUVbrnlluseq9frodfrmzK0VifWjCitA+u/9qTjpY0nAADfHc1GbqkJT0/o7uZRERFRczn3pVFqzYiHhwfi4+ORmJhY5/7ExESMHDnyms/7/PPP8dhjj2HdunW48847mzdSN3PWjCgoM1JYbsbrP50GAAyO7gAAeOfns0jLK3PnsIiIqAWc+9LIpPsq0Ixpmvnz5+Ojjz7Cxx9/jJMnT+KZZ55BRkYGZs+eDcA+xTJ9+nTn8Z9//jmmT5+Ot956C8OHD0dOTg5ycnJQXFzsunfRBpRYM7L+QCYqzFb0CfPDV7NHYHzPTrDYBKzZne7uoRERUTMpPjMCAJMmTcLy5cuxZMkSDBgwADt27MCmTZsQHR0NAMjOzq7Tc+SDDz6AxWLBE088gdDQUOdt7ty5rnsXbaAmGFFOn5Efj2UDAKYMi4JKpcKM0V0AAP9LyVJUBoiIqD2RY2akWWHTnDlzMGfOnAYfW7t2bZ3/37ZtW3NeQnIMCtsoL6uoEkcuFUOlAhJ6hwAARnYNQGRHT2QWVmJTajYeiI9w8yiJiKip2kVmpL3y9HCsplFIAWvicXsR8pDojujkay8WVqtVeGCQPQD56XjjlmoTEZG0yHE1DYORRlJaAevetAIAwM296i6xvqWXfVXUrnP5MCkk8CIiak8qnE3PmBlRHCUVsAqCgEMZRQCAeMcqGlGfMD8E+epRYbZif1qhG0ZHREQt4cyMyKhmhMFIIzn7jCggGLl0tRJ5pSZo1Sr0DTfWeUylUmF8T3u2ZOvpXHcMj4iIWoA1Iwom7tpbbRVgscp7Rc2hjKsA7FkQcfqpttHdAwGAmREiIhmqYM2IctXeirnKIu9g5LBjimZgVIcGHx8WY9/E8GROCYorqttqWERE5AJizYji96Zpj/Tamh+V3ItYT2Tbd0HuF2Fs8PEgPwNiAr0hCMDBi8yOEBHJhSAIzpoRb2ZGlEelUjmLWOVcNyIIAs5cKQUA9AzxveZxYnZk/wUGI0REcmGy2GBz7HfqxcyIMok798p5RU1eqQlFFdVQq4CunXyuedxQRzCSfPFqWw2NiIhaSOy+CtSsApUDBiNNoITMyKkce1akc6B3g8WrogGR/gCAY1nFqJZ5wS4RUXshrqTx1GmgUavcPJrGYzDSBAYP+Tc+c07RBF97igYAOgd4w8+ghcliw2lHAENERNLmrBeRUY8RgMFIkyih8ZkYWFyvXgSwt4bv78iOHLlU1MqjIiIiVyh3rKTxlFHxKsBgpEmUME1zLq8MANA96PrBCAD0j/AHABzJLGrFERERkauINSNyangGMBhpEoMCMiMXCyoAAJ0DvW54rDMzklncmkMiIiIXEYMRXwODEcWq2SxPngWdxZXVKCw3AwCiA7xveHz/SHsfkjO5pSirVaFNRETSVOr4W+0jo2W9AIORJhELgsRWu3KT4ciKBProG/UPNcjXgHB/TwgCkHqJ2REiIqkrq3IEIwadm0fSNAxGmkDcjlksEJKb9IJyAEDngBtP0YjE7AiLWImIpK+MmRHl85F5ZuSiIxhpzBSNiEWsRETyUcaaEeUTMyNyrZ9IF4tXm5QZ8QfAYISISA5Kq5gZUbyamhF5TtOINSPRgY3PjMSFG6FSAZeLq5BXamqtoRERkQtwmqYdELdjLpdtZqTpNSM+eq1zD5vUrKLWGBYREblIWVU1AMCH0zTKJTaRKZdhzUiF2YJcR2YjumPjMyMA0C/CXsR6lCtqiIgkzVkzwsyIcnk52uvKcTWN2OzM30sHo1fTlnz1C2cwQkQkB86aEWZGlEucg5PjaprmrKQR9XMUsR69VAxBEFw5LCIiciHWjLQDXnr59hm5dLUSABDZwbPJz+0d6getWoX8MhOyi6tcPTQiInIRLu1tB7zFaRoZZkayiuzBSHgzghGDToMewfaN9ThVQ0QkTYIg1HRg1bMDq2KJq2kqZJgZuSwGI/5ND0aA2kWsRa4aEhERuZDJYoPFZp9KZ82IgomracxWG8wWeW2Wd7nIPr0SZmxuMOIPAEjNYmaEiEiKajfk9HJs7CoXDEaawEtfc3HlVsQqTtOEtTgzwiJWIiIpKqvVfVWtVrl5NE3DYKQJdBo1PLT2H1m5jLqwVpqtKCw3A2hezQgA9Aj2hYdGjeLKamQUVrhyeERE5AJyXUkDMBhpMmcRq4y6sF4utmdFfPRa+DVzHtFDq0avMD8AwBEWsRIRSY5ce4wADEaaTNwsT1bBiHOKxgCVqvmpO7H5WSqLWImIJIeZkXakpvGZfKZpsq62bCWNSKwbYWaEiEh6Sh370sitxwjAYKTJxCLWMllmRloajPgDAI5nFcNqYxErEZGUFFfagxE/T3n1GAEYjDSZuLxXTqtpssRlvS0MRroF+cBTp0G52Yq0vDJXDI2IiFykpNL+uWRkMKJ83nr5bZbX0oZnIo1ahbhwexErO7ESEUmLMzNiYDCieN4yLGBtSSv4XxOnatiJlYhIWkqqxGka1owonrdeXsGIzSYgu9g1NSMA0N+xg+/hzKIWn4uIiFynxJEZ4TRNOyCu3y6VSTCSX2ZCtVWAWgUE++pbfL7B0R0AAMcvl8gmICMiag84TdOOiEumxOYyUidO0YT4GaDVtPxyh/l7IsxogNUm4AizI0REklFSxQLWdsPXEXGK67mlzpX1IqL4zh0BAAcvXnXZOYmIqGVKuLS3/fB11IzIpc+Iq3qM1CZO1TAYISKSDmcwwqZnyie3aZrLLuoxUlu8Ixg5dPEqm58REUmA1SY4axk5TdMO1EzTyCMYyXJRj5HaYkN84e2hQZnJgtM5pS47LxERNU/t0gFfFrAqn7g3jWyCERftS1ObVqPGIEd2ZF9agcvOS0REzSN2X/XUaeChld9Hu/xG7GY10zTyKGC97MIeI7WN6hYIANh9Lt+l5yUioqYrlnGPEYDBSJOJ67dNFhvMFpubR3N95SYLiirs/0DD/A0uPfdoRzCyL60A1VZp/xyIiJROzt1XAQYjTSbuTQNIf0WN2HnVz6B1+Rxi71A/+HvpUG62st8IEV2XhV9YWp2cu68CDEaaTKtRw8vDHpBIfarm0tXWmaIBALVahVFd7dmRXZyqIaJfMVts+OfPZzH6/35Bt0U/YMJb2/Dd0csQBK7Aaw1y7r4KMBhpFrks7xWX9bqyeLU21o0QUUNKq6rx8Kq9eCvxjPNL0fm8cjy57jBe+f4kA5JWUMTMSPsjlxU1rdHwrLYx3e3ByKGMIhRVmFvlNYhIXqqtNsz810EcyiiC0VOHtx7sj70Lb8bTN3cDAKzedQGrdqS5eZTKc9XxN7iDt4ebR9I8DEaaQS4t4bNaORiJ7OiF2BBfWG0Cfj6Z2yqvQUTysnzLGey/UAhfvRafzRyGB+IjEGr0xPyEnnjxrt4AgDc3n8aJyyVuHqmyXC23ByMdGYy0H+I0jdQLWFtjX5pfS+gTAgD46XhOq70GEcnDsaxirNh2HgDw2gP9EBdurPP446M6I6F3MKqtAhb+N5XTNS5UWG7/ctzBqx0FIytWrEBMTAwMBgPi4+Oxc+fOax6bnZ2NKVOmoGfPnlCr1Zg3b15zxyoZcqkZaY2GZ792W59gAMCOs3moNFtb5TWqqq3YfDwH7/5yFiu2ncPuc/mwsQ09kaQIgoCXvj0OQQDu7h+GO/uF1jtGpVLhld/EwctDgyOZRfg+NdsNI1UmcZqmo3c7qRlZv3495s2bh0WLFuHw4cMYM2YMJk6ciIyMjAaPN5lM6NSpExYtWoT+/fu3eMBS4KuX/jSN1SYgp6R1C1gB+xLfiA6eqKq2YfuZPJeeWxAEfLb/Ika+9gv+8O9kvLn5DF7/8TSmfrQfCct3IPlioUtfj4ia79sjl3Hw4lV46jT4yx2x1zwuyNeAWTd1BQAs23yGXyxcRJymaTeZkWXLlmHGjBmYOXMmevXqheXLlyMyMhIrV65s8PjOnTvj7bffxvTp02E0Ghs8Rm6cmREJT9NcKamC1SZAp1EhyFffaq+jUqlwu2OqZuPRyy47b4XZgt9/koxF/z2GwnIzQo0GPDAoAvf0D4OvQYtzuWV48P29+PJgpstek4iax2oT8I/EMwCAJ8Z3Rajx+l+AZoyJga9Bi7T8cvx8ivVmrlCTGWkHwYjZbEZycjISEhLq3J+QkIA9e/a4bFAmkwklJSV1blLi51g6JTaZkSKxXiTEaIBarWrV17pvYDgAIPH4FZesqjFZrJj172RsOXkFHho1XrirN3YsGI+3HuqPdyYPxO7nb8Y9/cNgE4AF/zmK/x6+1OLXJKLm25SajfSCCnTw0uF3o2NueLyPXoupw6IBAB9yZU2LWW2Cc2mvf3vIjOTn58NqtSI4OLjO/cHBwcjJcV0B49KlS2E0Gp23yMhIl53bFfy97MFIsYSDkcutsFvvtcSFG9Er1A9mqw3fpLQsO2Kx2jDvixTsPJsPLw8NPv/DMMwYHQOdpuafqp9Bh7cfHoDHRnYGADz3n1QczrjaotclouYRBMFZtPrYyBh4eTSuHfnjozpDp1EhKb2Qv78tVFxZDbEWWPx8kptmFbCqVHW/aQuCUO++lli4cCGKi4udt8xMaaXixaYy4r4vUnTJWbzq1Sav99DgCADAv/ddbPYcsM0mYOGGVPxwLAceGjVWTRuM+OiODR6rUqnw4l29kdA7GGarDfPWp6BcwtNmREq17UweTmaXwNtDg0dHRjf6ecF+BtzT355V/Wx/wzWH1DiFjnoRP4O2zhc3OWnSqAMDA6HRaOplQXJzc+tlS1pCr9fDz8+vzk1KxDSYlIORmsyIazfIu5bfxkfAV2+v5dh2pulzwIIg4JXvT+Kr5EtQq4B3Jg/EaEdTtWtRq1V486H+CDMacLGgAkt/ONnc4RNRM63cas+KTBkW1eQpginDogAA3x/Ndm70Rk0n93oRoInBiIeHB+Lj45GYmFjn/sTERIwcOdKlA5Myf0/pT9O0RY+R2nwNOjw81D6d9t7W803uH/DOz+fw8e4LAIDXf9sft8eFNOp5fgYd3njQvkrr030ZSLrAFTZEbeVgeiGS0gvhoVFj5pguTX7+oCh/dAvyQWW1FRuPuK4Avr0RMyNy7b4KNGOaZv78+fjoo4/w8ccf4+TJk3jmmWeQkZGB2bNnA7BPsUyfPr3Oc1JSUpCSkoKysjLk5eUhJSUFJ06ccM07cANxTk7KLdCzWnGTvGuZMboLDDo1ki9exeYTVxr9vDW7L+AfW+yV+C/e1Ru/jY9o0uuO6haIyY5A6KVvj8PKpYJEbUKsFXkgPhzBfk3PwqpUKjw8xP67++UBaU3Hy4n4WdRRpsWrQDOCkUmTJmH58uVYsmQJBgwYgB07dmDTpk2IjrbPFWZnZ9frOTJw4EAMHDgQycnJWLduHQYOHIg77rjDNe/ADfw97Re83GyF2SK9rbEFQWjTAlZRiNGAmaPt345e+f5EozrUfpGUgZc32gPTuRO6N6oSvyHPJvSEr0GLE9klXO5L1AZOZpfgl1O5UKvg7BvSHL8ZGA6dRoUjl4rZIr6ZxO6rcl1JAzSzgHXOnDlIT0+HyWRCcnIybrrpJudja9euxbZt2+ocLwhCvVt6enpLxu1WvgYtxHpdKU7VFFdWo9zRDbUtMyMAMHtcV4T7eyKzsBIvf3v8utM1n+67iOc3pAIAfjcqBvNu6d7s1w3w0eOZW3oAAN746bQkrwuRkqx0ZEXu6BuKzoHezT5PgI8et/a21xzyi0TzFJSZAAABPu0sGGnv1GqVc0VNcaX0pmrElTSBPh4w6DRt+to+ei3eeqg/VCrgq+RLWJZ4pl5AYrbYsPibY/jr/44BsC/xe+GuXi1ekTVtRDS6BfmgsNyMFVvPtehcRFIhCAIOpBfio51peOfns9hw6JLzw8ddLhaU4ztHk8M/jmt+VkT04GD7VM03KVmSzDZLXZ7j30NrNrhsbY1bEE71+HvqUFRRLckVNe6YoqlteJcAvHR3Hyz+9jj++cs5pGYV43ejYhDkp0dKRhFW7UhDWn45AGDeLd0xd0J3lywN12nUWHRHLzy+9gDW7EnHtBHRiOjQNkubiVrDjjN5WPLdCZzLLatzv1atwqQhkVhwWyyMbugr8cGONNgEYFzPTugT1vLO2mO6BaKTrx55pSb8ciq30QXsZJdbYg9GOjEYaX+MXh5AQYUkgxFxJU1bT9HU9ujIzlCrVfjbxhPYdjoP207X3bcm0McDS+/v50zPusq4np0woksA9qYV4K3NZ/CPSQNcen6itmCzCXhj82nnVIiPXosx3QNh9NQhNasYxy+X4LP9GfjlVC5WPhKPAZH+bTa2KyVV+M9Be9fjOeO6ueScWo0a9w8Mxwc70vD1oUsMRppIzIx08mEw0u6Iy3uLJFib4O7MiGja8GiM6NIRH+64gH0XClBSWY2oAG9MjAvB1GFR8DW4/hudSqXCX+7ohbvf3YX/Hs7CjNEx9bYxJ5IyQRCw6H+p+DzJXj/x6Iho/Om2nvCr9fuyL60ACzek4kJ+OSav2of3p8VjbI9ObTK+D7anwWy1YUjnDhga03BTwuZ4ID4CH+xIw9ZTuSgoMyFAxh+sbS3XsSlqkJ98f2asGWkmKS/vlUJmRNQtyBf/99t+2P7n8Tj8YgK+eWIUZo/t2iqBiKhvhBH3DQgDALy66WSTe54QudPfvjuJz5MyoVYBbz3YHy/fG1cnEAHsU6EbnxqNMd0DUVltxcx/HcDuc/mtPrbc0ip8tv8iAODpCc0vOG9Ij2Bf9A03wmITWrytRHtSVW1FSZV95WInn7ZpctkaGIw0k5Qbn7mjx4jU/CmhJzw0auw5X4BtZ/Ju/AQiCVi3P6NO878HrtNzx0evxepHh2BiXAiqrQJm/TsZx7KKW3V8H+5Ig8liw8Aof4zudv0Oyc0h9hj6+hA3v2ysfMcUjYdWDT9P+U52MBhpJqOEW8JnFFYAAKI6tt/izciOXnhsVGcAwGubTrERGkneoYyrWPytfYXZn2/r2ajmfx5aNZY/PAAjugSgzGTBY2sOINPx++9q+WUmfLrP3kPqaRcVnf/aPf3DoNOocPxyCU5ms+dIY+SW1tSLtMY1aSsMRppJzIwUSmyapqSqGlcdAVJUQPsNRgDgiXHdYPTU4fSVUnydzG9aJF2lVdWY+8VhVFsFTIwLwZwmLJfVazX4YHo8eoX6Ib/MhN+tPdAq+7z88+ezqKy2on+EEeNaqT6lg7cHJsTai9r5O9s4eaXyX0kDMBhpNrG5jLvX+/9aRoH9W1GAtwd89PJN2bmC0UuHp262V/u/lXgaFWbu6kvS9LfvTiCzsBLh/p54/bf9mvwN18+gw5rHhiDYT4+zuWV44rNDqLa6rl/H+bwy5866z90e26rfwMWpqf+lXHbpe1AqMTMi5x4jAIORZhOXUOWXSSsz4pyiaedZEZG914gnrpSYsHrnBXcPh6ien47n4MuDl6BSAcse6t/s4u4QowGrHx0CT50GO8/mY/ENOiA3xdJNp2CxCZgQG4SRrVArUtu4np0Q4O2B/DITdrDe64aYGWnnAn3FYERamZGLjsxIdDuuF6lNr9Xgz7f1BAC8v/285K4XtQ2L1YaqaqvkVlblllZhoWNLhD/c1AXDugS06Hxx4Ua8M3kgVCp7MezqXS0PwH88lo0tJ69Ao1Zh4R2xLT7fjeg0atw7IBwAC1kbI6/UsazXV74raQAGI80W6MiMFFVUSyqVmFFo72waFdD8vSKU5u5+YegXYUS52Yq3t5x193CojWQWVmDZ5tNI+Md29PjrD4h94Uf0fWkzHvloP748kOn2tuOCIOC5/xxFYbkZvUL9MP/WHi457629g7Hojl4AgL9vOonNx3Oafa7iimq88M1xAMDssV3QLcjXJWO8kQfi7cHIlhO5kmyfICVZRfZgJNSfwUi75O+pg0ZtnzctkNBUDVfS1KdWq7Bwov2P87qkDJzPK7vBM0jOSquq8cp3JzD+zW1455dzOHOlDOJiqjKTBbvO5WPB10dx0+tb8e2Ry27LlqxLysDW03nw0KixfNIA6LWu20dqxugYTB0WBUEA5n6R0qwlv4Ig4C//TUVeqQldOnnjqZtd21fkevqEGdEr1A9mqw0bj7DnyPVkXbX/zY+QeSsHBiPNpFarEOBtL2KVUurfOU3DmpE6RnQNwITYIFhtAl774ZS7h0Ot5FhWMe54Zyc+2nUBFpuAEV0C8PbDA7B34c1IfSkBP84bgwW390SQrx45JVV4+vPD+OOnh1pl9cn1nM8rwyvfnQQALLi9J3qGuDbjoFKp8NI9fZxN0Wb86wDSHftBNdaHO9PwfWo2dBoV3nywf5tvuvnAIHt25D+Hstr0deVEEARJNblsCQYjLRDoI626EbPF5mwFz5qR+p6fGAuNWoXEE1fwQ2q2u4dDLvZDajbuX7kHmYWViOjgiTWPD8HnfxiOeweEI9ToCV+DDrEhfpgzrht2Pjcez9zSA1q1Cj8ez8F97+1GWhtlzMwWG+Z+cRiV1VaM6haA342KaZXX0WnUeG/qIPQI9sGVEhMe/GAvTuU0rnfH90eznUH7i3f3waCoDq0yxuu5b2A4tGoVjmQWtXozN7m6WlGNqmr7dGOIkdM07VZNEas0pmmyiiphEwCDTi37yurW0D3YF7Nu6gIAWPS/Y5IJIqnl/nc4C0+sOwSzxYYJsUH4/qkxGN8z6JrH67UazL2lO/47ZxRCjQak5ZXjvvd2Y8/51m+p/tbm0ziWVQJ/Lx3eenAA1OrWWybrZ9Dh05nDEBvii7xSE+5fsQf/PXz9otDP9l/E018chk0AJg+NxCPDolptfNcT6KPHxL6hAIB/7Ulvs9cVBAE/pGZj+sdJGLn0Z9z81jYs2XhCcm0cgJpu25189W2euXI1BiMtEOgjrWma2vUicu7E15rm3tIdsSG+KCw349mvjsDGzqyyt/5ABp75MgU2AXgwPgKrpg+G0atxy2P7Rhjx7ZOjER/dASVVFjz6cRL+d7j1pgV2n8vHqp1pAID/e6Bfm3ybDfI1YP0fRmBk1wBUmK14Zv0RTFu9HwfSC531MoIg4FhWMWb+6yAW/fcYrDYBv42PwCv39XXr35LHRnYGAHxz5HKbBANFFWY8tuYA/vjZIew4k4fLxVVIyyvHx7svYMKy7Ui6UNjqY2gKpUzRANy1t0WcvUZKJRKMFDhW0nTkSppr0Ws1WPbQAPxmxW5sO52Ht38+i2dctIrh1wRBQLbjj1mF2QKtRoWojl6IDvCGTiOf7wHVVhtOXC7BlZIqVFls6OSjR9cgb0ksJfxkbzpedKz2eGR4FJbcE9fkTEMnXz0+mzkMf/rqCL4/mo1561OQVVSJOeO6uvSDOLOwAk99fhiCAEweGoXb+oS47Nw3YvTS4d8zhuG9refwzs9nsfNsPnaezUdHbw+EGg3ILzPhSon975hOo8L8W3ti9tgubv9SMyjKH/0ijDh6qRjr9mfgKRdvzldbdnElpq9OwtncMnho1Zg5Oga39A5GbokJb/98FiezSzD94/1Y9/vhbpm2aogYjMi9eBVgMNIiUqsZOZ9nD0a6dGIwcj29w/zw6m/64k9fHcHbP59FTKA37hsY7pJzC4KAvWkF+CE1B1tOXkF2cVW9Y3z1WoyPDcJDgyMxqluA2//gN0QQBOw9X4B/77uI7WfyUGG21jsmNsQXvxkYjklDIuHv2KupLX24Iw1/32QvAp0xOgZ/vbNXs3+WBp0G/3x4IML9PbFqRxre+Ok0sooqseSePtC6IHAsN1nw+08OorDcjL7hRrx4V+8Wn7OpNGoVnp7QHfcNCMe7W8/iu6PZKCw3o7DcPs2s06hwe1wonr65G7oHt80S3htRqVT43agYzFufgrV70jFjTAy8PFz/sXW13IzJq/YhvaACwX56rH18KHqF+jkfH9ezE/7w72TsOJOHP36ajO+eGiOJqfDLzsyI+78YtBSDkRYI9LX/Ac6TTDBiL8Dr1snHzSORvgfiI3DscjHW7E7Hn746AoNOg9vjmv9Ntaraiv8ezsKa3Rdw5kpNIaRWrUJUgBeMnjqYqm1ILyhHqcmCb49cxrdHLqNvuBEv3NUbQ2M6uuJtucSxrGIs2XgCSek1KWl/Lx2iA7xh0KqRW2rChfxynMopxdIfTmH5lrP4/ZgYzBrbFd5ttAXBP38+i7cSzwAAnhjfFc8m9GxxUKdWq/CXO3oh3N8TL208jnX7M5BdVIl3pwxq0fuqttow94sUnMopRaCPHqumx8PTw33z+1EBXnj9t/2x5N44nLlSioJyM/wMWvQONbp1XNdyV79QLEs8g4zCCqzbn4GZY7q49PwmixWzPk1GekEFwv09sX7WcER0qLsAwKDTYMXUQbj33V04n1eOxd8ew4qp8S4dR3OINSPhzIy0b8F+9mi0oW+/7nA+1/4h2DWImZHGeOHO3iiuqMaGw1mY81kyXryrNx4d2blJH2pFFWb8e+9FrN2TjgLHN0wvDw3u6R+G2/qEYETXgDqFZTabgMOZRfgmJQtfHsxEalYxHvpgLx4YFIG/3tkLHbzbPsMgqqq24v9+PIW1e9IhOAqhfxsfgYeHRKF3qF+d6Y+CMhMST1zBJ3sv4kR2Cd755RzWJWXixbt74+5+oa2W7REEAf/342m8v/08AOBPt/Zweer+0ZGdEWI04OnPD2Pr6Tw8vGofVj82uFnTUlabgPlfHsGWk1fgoVXjg2mDEGqUxgeHQadBvwh/dw/jhrQaNZ4c3w0Lvj6K97en4ZHh0S4r1hQEAX/ZcAxJFwrhq9dizeND6gUiIh+9Fu9OGYS7/7kLm1JzsOXEFdzSO9gl42iuC47l2tGB8v+brxKk1h+5ASUlJTAajSguLoafn9+Nn9BGLhaUY+wb22DQqXFyye1uTbeXmyzos/gnAMDhF25164eanFisNizckIqvHDuEJvQOxkv39LlhQVh6fjk+3XcR65IynFMYER088djIznhoSCT8GrG/SEGZCW9uPo3PkzIBACF+Bix/eACGt7AleHOcvVKKpz4/jFM5pQDsW7k/PzH2hj8HQRDw0/EcvPbDKaQ7etzc2jsYr9wX5wzWXcVqE/DX/6U6f16L7uiF39/k2m/JtR3KuIqZ/7JPrUR08MTax4c0qQNpucmCuV8cxpaTudBpVPhgWjxujnXvh5dcVVttGPfGNmQVVeLPt/XEE+O7ueS87209hzd+Og2NWoWPHxuCsY3Yjfi1H07h/e3n0SXQG5ufuckl03jNYbUJ6PXijzBbbNj+53GIlmjX7cZ+fjMYaQGTxYqef/0RgPsDgNRLxbj73V0I8PZA8gu3um0cciQIgrNOwGIT4KFV4+5+YZgYF4K4cCOMnjpUmC3IKKzAgfRC/HT8CpIvXnU+v1eoH2aP7YI7+4Y26w/ToYyrePbLI0jLL4dKBTw1vhuentC9zf7IfXkgEy98cwwmiw0B3h5488H+GB977WWxDTFbbFi57Tze3XoW1VYBvgYtXrirNx6Mj3BJkF5ptuJPX6VgU2oO1Crg1d/0xcNDW3/JaXp+OR5bk4T0ggp4edj3OZo+orOz+/K1HL9cjD99eQSnckqh16rxzuSBbVqwqkT/O5yFeetT4OWhwS9/GtfilUjfH83GE+sOAQD+dl8cpg2PbtTzSquqMfaNbSgsN2Pp/X0xuQ3+HTYks7ACY17fCg+NGif/dvsN/026C4ORNjLk71uQV2rCd0+NRly40W3jEH9Rh8Z0xJezRrhtHHJ2KqcEL/zvGA6kX73hsWoVMKZ7Jzw+qjPG9ujU4g/ccpMFL288ji8P2jM0Qzt3xNuTB7RqSr/SbMWL3xxzZoXGdA/EWw/1b9EqmVM5JVjwn6M4esnepGp4l4545b6+6BbU/DqmjIIKzPo0GSezS6DTqPD2wwNxh6P/RFsoKDPhiXWHsC/NXkPTtZM3Zo/titvjQurtsHshvxyrd6Xh86RMWG0CAn08sGr6YMmsvpAzQRDw2/f3IvniVdzTPwzvTB7Y7HMlX7yKKR/ug8liw+OjOmPx3X2a9PzVuy7gb9+dQLCfHtueHe+WWpttp3Px2JoD6B7kg8T5Y9v89RuLwUgbuffdXThyqRirpsUjwY3ffN786TTe3XoOk4dGYun9/dw2DrkTBAEpmUX48uAlHEwvxLm8Moi/IZ189egfYcSIroG4q1+oy6chAOCblCws+u8xlJks6OClw7JJA67bvKu50vPLMfvTZJzKKYVaBfwpoSf+OLarSxpwWaw2rN51Af/YcgZV1TboNCrMHtsVfxzXtUkrIQRBwIZDWXh543GUVFkQ6OOBd6cMcss0ls0mYF1SBt746TSKK+2t4z20avQJ80OYvyfMFhvOXil1TlUBwB19Q/DyPXGSWHWhFMey7BlgQQDemzIId/ZrelCalleGB1buwdWKakyIDcKq6YObnFUwWay4+c3t9hVX9/bB9BGdmzyOllqz+wJe3ngCt/UJxgfTBrf56zdWYz+/WcDaQmH+njhyqdi5xMpdxDbPPSWyJE+uVCoVBkZ1wEDHN1mrTUBltRU6jcqlG5ldy70DwtEvwh9PrjuE45dL8PiaA5g1tgueTejpst4km1Kz8dx/jqLUZEGAtwfemTwQo7oFuuTcgL3gcNbYrrijbygWf3scv5zKxT9/OYfP9mdg5pgYTB0WDaPn9Wtqki4UYlniaWc2YmCUP1ZMdV/xp1qtwiPDo3HvgDD8e99F/OfgJaTll+NwRhEOZxTVHKcCxvbohN/f1AUju7ruZ0p2ceFGzBnXFe9tPY+FG45iQJR/k1aSXCmpwqNrknC1ohr9I4z455SBzZre0Gs1mD22C1745jhW7UjDlKFRbV47kuZs5aCM1ZMMRlpILPC77OYVNSez7YWHtdfGU8tp1Cr4tNFyVVFMoDe+/uNILN10Ev/aexEfbE/D/rRCvP7bfujRgmCzuKIaL208jv86OowOju6Ad6cMarUuoJEdvbD60cH46XgOXt10ChmFFXj9x9NYvuUsJsQGYWTXAPQK9UOAjx4Wqw1ZRZU4klmMH45lOwtp9Vo15t3SAzPHxEiiUZyvQYc547rhj2O74nxeOc5cKUVOcRUMOg1CjQbEd+7QqOJlar55t/TArnMFOJJZhJn/Ooj1s4Y36md+6WoFpn60H5mFlYjq6IXVjw1pUc+SBwdHYvmWs7h0tRLfp2bj3gGu6VXUWGn59tWTMQpYSQMwGGkxMRjJcmNmpLiy2vn6sSEMRpTAoNPg5XvjMLxLABZ8fRQpmUW44+2d+MNNXTB7XNcmfeBZbQK+OpiJNzefQX6ZCWoVMHtsVzxza49W/4BXqeyNtG7pFYxvUi5j1Y40nL5Sih+O5eCHYznXfJ6HRo0HB0dgzvhukuyhoFKp0C3Ip0W1MNQ8Oo0a704eiN+s2GPviro6CR8/NgQdr7OA4GB6IZ5YdwhXSkyI7OiJz2YOczatbC6DToPHR3XGm5vPYOW287inf1ibragUBMH5BbQlX1CkhMFIC4U7Ot+JzWfc4bTjW2SY0dDoPTlIHib2DUX/SH+8+M1xbDl5BSu2ncen+y7isVExmDI06rpZjUqzFRuPXMZHu9Kcjdi6dPLGmw/2b/OCSq1GjQfiI3D/oHAcv1yCLSftK5IuFlSgoMwED60aHb090CfMiFHdAnB7n1D+W6ZriuzohbWPD8Ejq/cjJbMId72zE3+/vy/G/aqYvLSqGh9sT8PK7edhtQnoFuSDf88Y6rLpvmnDO2PltvM4lVOKbafzmrwKrbmyi6tQWG6GRq1CbAiDEQKcDXIyCytucGTrEetFYjlFo0hh/p74cHo8Np+4gjd+Oo1zuWV45+ez+OcvZxEf1QHxnTsgJsAbRk8dqixWZBZW4uilIuw6l+/cXtzoqcPcCd3xyPBoeGjdN92hUqkQF25068ozUoa4cCP+M3sEZvzrIC4WVODxNQcQG+KL4V0C4KPXIi2/DNtO12xlcE//MCy9v69LuwQbvXSYMiwKH+68gJXbzrdZMHL8sv1vfvcgH9nv1itiMNJC4nxdQbkZRRVmt+zRccLxD1MpETLVp1KpcFufENzaKxg/HMvB2j0XcCD9Kg5etN+uJbKjJx4ZFo2Hh0Qx00CK0y3IF5ueHoNliWfw6b6LOJVT6qw3qjnGB88m9MRtfYJbZRplxuguWLsnHUnphTiYXojBnVt/a4djWfal833ClBPUMxhpIW+9FqFGA7KLq3A+rxzx0W0fjKRkFgEA+kUo5x8mNUytVuHOfqG4s18osooqsftcPlIyi5BTXIWSymoYdBoE+ekRG+KLUd0C0TvUT5Ib8RG5irfe3mDvyfHdsPV0Lk7nlKKq2opgowHDYgIwKMq/VX8HQowG3D8wAusPZmLFtvP4+LHWD0bEzEifMOVkwxmMuECXTt6OYKQM8dFtOxdfbrLgzBX7N4EBkWys1J6E+3viocGReGhwpLuHQuR2Hbw9cP+gCLe89uxxXfFVciZ+OZWLk9klrbqqURAEHL1UBACKmu50/1o5BejqWOctrvtuS8eyimET7PuatNYSTSIiuraYQG9MdHQFXrntfKu+1oX8cuSW2ou+lZQNZzDiAl0cdSNpeWU3ONL1xCmaAZH+bf7aRERk98exXQEA3x29jIsFrffF1NkIMNJfMcWrAIMRl+jq6DVw3g3BiNj9sT+DESIit4kLN2Jsj06wCcAHO9Ja7XX2pRUAgFu2RWhNDEZcoLtjW/H0ggpUOpaRtQWbTcD+C/Z/mENjWr9oioiIrm3OOHt25D8HLyGnFbpy22wCgxG6tmA/PTr56mG1CTh+ubjNXvdUTimuVlTD20OjqLlDIiI5GhrTEUM7d4TZasPyLWdcfv4jl4qQW2qCt4cGA6P8XX5+d2Iw4gIqlQr9HcHAkUttF4zsOZ8PABgS01ES+3YQEbVnKpUKz03sCQD48mAmzuWW3uAZTfPjcfsWCuNjgxRVLwIwGHGZfhH+AOBcctUW9py3p+tGdlVWuo6ISK7iozsioXcwbAKwdNMpl51XEAT85NjPaWJcqMvOKxUMRlxEnCY52kaZkQqzBbvP2TMjY7p3apPXJCKiG1tweyy0ahV+PpWLH6+zIWRTHEi/ivSCChh0aozrqby/+QxGXKS/IzNiXwPu+sKlX9txJg8miw2RHT3ZBp6ISEK6Bflg1tguAICXvj2O0qrqFp/zX3vTAQD3DQh36f46UsFgxEU6eHsgLtzedW/X2fxWf73Nx68AAG7rHcJ230REEvPUzd3ROcALOSVVWPzNcQiC0OxzXbpa4ZyimTYi2lVDlBQGIy50k2O6ZMeZvFZ9nUqzFYkn7MFIQp+QVn0tIiJqOoNOg9d/2x8atQobDmfhiwOZzT7XW5vPwGITMKpbgKI2x6uNwYgLje3hCEbO5sNma34UfCObUrNRarIgqqMXBrfxXjhERNQ4Q2M64tkE++qaxd8cd9b5NcWB9EL8LyULAPD87b1cOj4pYTDiQoOiO8DXoEVhuRn7HM3IWsMXBzIAAJOGREKt5hQNEZFUzbqpCybGhcBstWHmvw5if1rjPxuulpsx74sUCALwYHwE+iq4nxSDERfSadS4u38YAOCrg5da5TUOZVzFgfSr0KpVeMBNO1QSEVHjqNUqLH94AMZ0D0RltRWPrN6PdfszblhDUlRhxqNrkpBVVImojl5YfE+fNhqxezAYcbEH4+0Bwg/HslFc2fIK6l/7589nAQD3DwrnLr1ERDKg12qwatpgTIwLQbVVwF/+m4qpH+1HSmZRvaBEEARsPZ2LO9/ZhaOXitHR2wOrHx0MHwWuoKlN2e/ODQZE+iM2xBenckqxetcFzL+1h8vOveNMHraezoNaBcwZ181l5yUiotbl6aHBiqmD8P72NPxjyxnsOV+A+97bjS6B3hgQ6Y9OfnoUlVcjKb0QF/Ltu/5GdfTCqunx6B6s/PYNKqEl643aSElJCYxGI4qLi+Hn5+fu4dzQptRszPnsELw9NNi+YDwCffQtPmeZyYI739mJiwUVeHxUZyy+W9kpOyIipcosrMCyxDPYlJoNk8VW73GDTo1pw6Px5M3dYfTUuWGErtPYz28GI63AZhNw73u7kZpVjFt7B2PVtPgW9QKx2QQ89flhfJ+ajVCjAYnzxyo+ZUdEpHSlVdXYfS4f5/PKUVBmho9Bi14hvhjTo5Ni/sY39vNbGe9WYtRqFf7vgX64773dSDxxBe/8fA5zb+nerHNZbQJe+OYYvk/Nhk6jwrtTBinmHykRUXvma9DhdgXuM9McLGBtJb3D/PDXu+xrwv+x5Qxe3ngc5gbScdeTVVSJRz9Owrr9GVCpgNd/2w/x7CtCREQKw6/YrWj6iM4orbLgjZ9OY83udPxyKhczR8fg7v5h8PfyaPA5VpuAo5eK8OXBS9hw6BJMFhs8dRq88WA/3NUvrI3fARERUetrVs3IihUr8MYbbyA7Oxt9+vTB8uXLMWbMmGsev337dsyfPx/Hjx9HWFgYFixYgNmzZzf69eRWM/Jrm4/nYNH/jiGv1AQAUKmA7kE+6BLogw7eHgAEVFXbkHW1EqdySlBSZXE+d1hMR/z9N3HoFqT8amoiIlKWVqsZWb9+PebNm4cVK1Zg1KhR+OCDDzBx4kScOHECUVFR9Y6/cOEC7rjjDvz+97/Hp59+it27d2POnDno1KkTHnjggaa+vCwl9AnBqG6B+OJAJv6TfAkns0tw5koZzlwpa/B4X70WY3t2wtRh0RjepSM3wiMiIkVrcmZk2LBhGDRoEFauXOm8r1evXrjvvvuwdOnSesc/99xz+Pbbb3Hy5EnnfbNnz8aRI0ewd+/eRr2m3DMjv5ZbWoXUS8W4XFSJwvJqaNT27q0hRgO6dvJBbIgvtBqW8xARkby1SmbEbDYjOTkZzz//fJ37ExISsGfPngafs3fvXiQkJNS577bbbsPq1atRXV0Nna7+GmqTyQSTyVTnzShJkK8BE3qxeyoRERHQxNU0+fn5sFqtCA4OrnN/cHAwcnJyGnxOTk5Og8dbLBbk5ze8g+HSpUthNBqdt8jIyKYMk4iIiGSkWXMBv65hEAThunUNDR3f0P2ihQsXori42HnLzMxszjCJiIhIBpo0TRMYGAiNRlMvC5Kbm1sv+yEKCQlp8HitVouAgIAGn6PX66HXt7yFOhEREUlfkzIjHh4eiI+PR2JiYp37ExMTMXLkyAafM2LEiHrHb968GYMHD26wXoSIiIjalyZP08yfPx8fffQRPv74Y5w8eRLPPPMMMjIynH1DFi5ciOnTpzuPnz17Ni5evIj58+fj5MmT+Pjjj7F69Wo8++yzrnsXREREJFtN7jMyadIkFBQUYMmSJcjOzkZcXBw2bdqE6OhoAEB2djYyMjKcx8fExGDTpk145pln8N577yEsLAzvvPNOu+kxQkRERNfHXXuJiIioVTT285udtYiIiMitGIwQERGRWzEYISIiIrdiMEJERERuxWCEiIiI3IrBCBEREblVk/uMuIO4+lhpu/cSEREpmfi5faMuIrIIRkpLSwGAu/cSERHJUGlpKYxG4zUfl0XTM5vNhsuXL8PX1/e6uwM3VUlJCSIjI5GZmclmajLBayY/vGbyw2smP1K9ZoIgoLS0FGFhYVCrr10ZIovMiFqtRkRERKud38/PT1IXj26M10x+eM3kh9dMfqR4za6XERGxgJWIiIjcisEIERERuVW7Dkb0ej0WL14MvV7v7qFQI/GayQ+vmfzwmsmP3K+ZLApYiYiISLnadWaEiIiI3I/BCBEREbkVgxEiIiJyKwYjRERE5FYMRoiIiMit2nUwsmLFCsTExMBgMCA+Ph47d+5095DapaVLl2LIkCHw9fVFUFAQ7rvvPpw+fbrOMYIg4KWXXkJYWBg8PT0xbtw4HD9+vM4xJpMJTz31FAIDA+Ht7Y177rkHly5dasu30i4tXboUKpUK8+bNc97H6yU9WVlZeOSRRxAQEAAvLy8MGDAAycnJzsd5zaTFYrHgr3/9K2JiYuDp6YkuXbpgyZIlsNlszmMUdc2EduqLL74QdDqd8OGHHwonTpwQ5s6dK3h7ewsXL15099Dandtuu01Ys2aNcOzYMSElJUW48847haioKKGsrMx5zGuvvSb4+voKX3/9tZCamipMmjRJCA0NFUpKSpzHzJ49WwgPDxcSExOFQ4cOCePHjxf69+8vWCwWd7ytdiEpKUno3Lmz0K9fP2Hu3LnO+3m9pKWwsFCIjo4WHnvsMWH//v3ChQsXhC1btgjnzp1zHsNrJi2vvPKKEBAQIHz33XfChQsXhK+++krw8fERli9f7jxGSdes3QYjQ4cOFWbPnl3nvtjYWOH5559304hIlJubKwAQtm/fLgiCINhsNiEkJER47bXXnMdUVVUJRqNReP/99wVBEISioiJBp9MJX3zxhfOYrKwsQa1WCz/++GPbvoF2orS0VOjevbuQmJgojB071hmM8HpJz3PPPSeMHj36mo/zmknPnXfeKfzud7+rc9/9998vPPLII4IgKO+atctpGrPZjOTkZCQkJNS5PyEhAXv27HHTqEhUXFwMAOjYsSMA4MKFC8jJyalzvfR6PcaOHeu8XsnJyaiurq5zTFhYGOLi4nhNW8kTTzyBO++8E7fcckud+3m9pOfbb7/F4MGD8eCDDyIoKAgDBw7Ehx9+6Hyc10x6Ro8ejZ9//hlnzpwBABw5cgS7du3CHXfcAUB510wWu/a6Wn5+PqxWK4KDg+vcHxwcjJycHDeNigD7HOj8+fMxevRoxMXFAYDzmjR0vS5evOg8xsPDAx06dKh3DK+p633xxRc4dOgQDhw4UO8xXi/pSUtLw8qVKzF//nz85S9/QVJSEp5++mno9XpMnz6d10yCnnvuORQXFyM2NhYajQZWqxV///vfMXnyZADK+z1rl8GISKVS1fl/QRDq3Udt68knn8TRo0exa9eueo8153rxmrpeZmYm5s6di82bN8NgMFzzOF4v6bDZbBg8eDBeffVVAMDAgQNx/PhxrFy5EtOnT3cex2smHevXr8enn36KdevWoU+fPkhJScG8efMQFhaGRx991HmcUq5Zu5ymCQwMhEajqRcZ5ubm1osyqe089dRT+Pbbb7F161ZEREQ47w8JCQGA616vkJAQmM1mXL169ZrHkGskJycjNzcX8fHx0Gq10Gq12L59O9555x1otVrnz5vXSzpCQ0PRu3fvOvf16tULGRkZAPg7JkV//vOf8fzzz+Phhx9G3759MW3aNDzzzDNYunQpAOVds3YZjHh4eCA+Ph6JiYl17k9MTMTIkSPdNKr2SxAEPPnkk9iwYQN++eUXxMTE1Hk8JiYGISEhda6X2WzG9u3bndcrPj4eOp2uzjHZ2dk4duwYr6mLTZgwAampqUhJSXHeBg8ejKlTpyIlJQVdunTh9ZKYUaNG1Vsuf+bMGURHRwPg75gUVVRUQK2u+xGt0WicS3sVd83cVDjrduLS3tWrVwsnTpwQ5s2bJ3h7ewvp6enuHlq788c//lEwGo3Ctm3bhOzsbOetoqLCecxrr70mGI1GYcOGDUJqaqowefLkBpewRURECFu2bBEOHTok3HzzzZJcwqZEtVfTCAKvl9QkJSUJWq1W+Pvf/y6cPXtW+OyzzwQvLy/h008/dR7DayYtjz76qBAeHu5c2rthwwYhMDBQWLBggfMYJV2zdhuMCIIgvPfee0J0dLTg4eEhDBo0yLmUlNoWgAZva9ascR5js9mExYsXCyEhIYJerxduuukmITU1tc55KisrhSeffFLo2LGj4OnpKdx1111CRkZGG7+b9unXwQivl/Rs3LhRiIuLE/R6vRAbGyusWrWqzuO8ZtJSUlIizJ07V4iKihIMBoPQpUsXYdGiRYLJZHIeo6RrphIEQXBnZoaIiIjat3ZZM0JERETSwWCEiIiI3IrBCBEREbkVgxEiIiJyKwYjRERE5FYMRoiIiMitGIwQERGRWzEYISIiIrdiMEJERERuxWCEiIiI3IrBCBEREbnV/wNpTwivNkVZTgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Beta\n",
    "\n",
    "plt.plot([func.betafn(tc.tensor([t])).detach().clone() for t in func.tspan])\n",
    "plt.title(r\"$\\beta (t)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fccff45f3a0>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEKklEQVR4nO3deXxU9b3/8ddsmewDCYQkJLLvCCKLRlERBYvW5bZVqxZw6e/WFlBrtRW8j6qtXqze3tvFNvW6YL0ieFGw2KupVASEFgUEQVD2nSxAyJ5MMjPn98fJTBZISEKSM2Hez8djHknOnDPzHc+l877f7+f7/doMwzAQERERsYjd6gaIiIhIZFMYEREREUspjIiIiIilFEZERETEUgojIiIiYimFEREREbGUwoiIiIhYSmFERERELKUwIiIiIpZSGBERERFLdakwsmbNGm688UbS09Ox2Wy8++67rbr+ySefxGaznfaIi4vrmAaLiIjIWXWpMFJeXs7o0aN54YUX2nT9I488Qm5uboPH8OHDufXWW9u5pSIiItJSXSqMTJs2jaeffppvfetbZ3y+urqan/70p/Tu3Zu4uDguueQSVq1aFXo+Pj6e1NTU0CM/P58dO3Zw3333ddInEBERkcacVjegPd1zzz0cOHCAxYsXk56ezrJly/jGN77Btm3bGDRo0Gnnv/zyywwePJgrrrjCgtaKiIgIdLGekebs3buXRYsWsWTJEq644goGDBjAI488wsSJE1mwYMFp53u9XhYuXKheEREREYudNz0jn3/+OYZhMHjw4AbHvV4vycnJp52/dOlSSktLmTFjRmc1UURERM7gvAkjgUAAh8PBpk2bcDgcDZ6Lj48/7fyXX36Zb37zm6SmpnZWE0VEROQMzpswMmbMGPx+PwUFBWetAdm/fz8ff/wxy5cv76TWiYiISFO6VBgpKytjz549ob/379/Pli1bSEpKYvDgwdx1113MmDGDX//614wZM4YTJ06wcuVKLrzwQq6//vrQda+++ippaWlMmzbNio8hIiIi9dgMwzCsbkRLrVq1iquvvvq04zNnzuS1116jpqaGp59+mtdff52jR4+SnJxMVlYWTz31FBdeeCFgDuf06dOHGTNm8Mwzz3T2RxAREZFGulQYERERkfPPeTO1V0RERLomhRERERGxVJcoYA0EAhw7doyEhARsNpvVzREREZEWMAyD0tJS0tPTsdub7v/oEmHk2LFjZGZmWt0MERERaYPDhw+TkZHR5PNdIowkJCQA5odJTEy0uDUiIiLSEiUlJWRmZoa+x5vSJcJIcGgmMTFRYURERKSLOVuJhQpYRURExFIKIyIiImIphRERERGxlMKIiIiIWKpVYSQ7O5tRo0aFCkmzsrL44IMPmr1m4cKFjB49mtjYWNLS0rjnnns4efLkOTVaREREzh+tCiMZGRk8++yzbNy4kY0bNzJ58mRuvvlmtm/ffsbz165dy4wZM7jvvvvYvn07S5YsYcOGDXz/+99vl8aLiIhI19eqqb033nhjg7+feeYZsrOzWb9+PSNGjDjt/PXr19O3b18eeOABAPr168cPfvADnnvuuXNosoiIiJxP2lwz4vf7Wbx4MeXl5WRlZZ3xnMsuu4wjR47w/vvvYxgG+fn5vP3229xwww3NvrbX66WkpKTBQ0RERM5PrQ4j27ZtIz4+Hrfbzf3338+yZcsYPnz4Gc+97LLLWLhwIbfffjtRUVGkpqbSrVs3fv/73zf7HvPnz8fj8YQeWgpeRETk/GUzDMNozQXV1dUcOnSIoqIi3nnnHV5++WVWr159xkCyY8cOrr32Wn784x9z3XXXkZuby6OPPsr48eN55ZVXmnwPr9eL1+sN/R1cTra4uFgrsIqIiHQRJSUleDyes35/tzqMNHbttdcyYMAAXnzxxdOemz59OlVVVSxZsiR0bO3atVxxxRUcO3aMtLS0Fr1HSz9Mq9VUQUkuVBaCzQE1leCrrnveGQWumPY/7oyChF6QmA6u6Pb7PCIiImGkpd/f57w3jWEYDXox6quoqMDpbPgWDocjdJ2lig7D5/8Duz+E0mNQXQGGH/w+83mbDewOsNnBCLTf8YAfHFEQnwIZE2DwN6B7XwUbERGJWK0KI/PmzWPatGlkZmZSWlrK4sWLWbVqFTk5OQDMnTuXo0eP8vrrrwPm7Jv/9//+H9nZ2aFhmoceeogJEyaQnp7e/p+mpWqqYPcKCr74AEdpPjGBIlx4CW7jYwMMAH/Dy1pyPPR7E8cNbOD3QmEpFO7Dv/VtzNtgAIHanzbMch577bGmjhsEbE7KXckc81zEvh6TqXHFYg+YQchmsxGwu6hxROPyV2IP+Bp8xoDdhc8ZjdNfBY4oKuL6YsR4cNhtuOx2HHYbTocNl8P8PdrlIC7KQUyUg7goJ7Fu82eMy4Hd3vwmSCIiIk1pVRjJz89n+vTp5Obm4vF4GDVqFDk5OUyZMgWA3NxcDh06FDr/7rvvprS0lBdeeIGf/OQndOvWjcmTJ/OrX/2qfT9Fa1WXQ2ku5RVV+GvAjh074KDuKz/4M6ilx42zHA9gYA8dM3BQA9SEXq/+ufWd6XgAcOAl0VdOQuVhBuR9gB8bTgxsBDCw4cOOgQ0bRui42b665wzslOHmQKAXH/jH8xV9qcRNDNW46rWtBmeTx/2uWHo4/CREGyS4ncRFuYiJiSYqOo5uUT48LoOEaAexUS7i3Q66JSTQ3ZOAO1DdcUNjMd0gOhFik9VrJCISxs65ZqQztHvNSE0VbHmTk+tew1aWT4zvFM4z9Yw083t9bTkO1MaCur8bBx1acNzehrbS6LkA9tpgEwwwUfixYyeAo0F4MaOMDQOHYWALNB14WhqE7A4Dl9lfRKjnx27HYQO74cMBbRsac7jAnQjd+kGvYZB5CcT1Mp9va+BxOMGTYYYcERE5q06rGemSXNEwaArJpXm1NSNG59SMBPzg80LtV7rDHCzp5A9/JvXbYBDFmWuAGrCB4WhZ6+sHpOC7RQWHuvx1z4UCkr/uPMNR+3ztsWCvEo3OPY2/GqorofQ4HF4PGxfU3oc23j9HlNnDkpgJg6+HlKFnDi+q4RERabXI7BkJqqmC0jyoKDS/cKorzS+xoOCXTXseLzkCez6CoxugrMAMKDbOIdgE+01sNP3N3BpN9am0r2CFTONjoZ8BMIwzDI3VhhN7wz87gd2ccWXU/jd2uE6/RzabGVoSUqHvROg/Cdzd615CQUVEIkynTe3tDB0WRqwSDEEleea37LkGm7wtUH7C/P/QW9K7U/85A/BVEeytIUx6a+qHFYPTC4dtgGEzP1owqHReMGmODbMlNqhf1NtUUFFAEZHzmMJIpGhr7079545sgi/fhhO7obqsNpO085TmMwWhVr6WHzt+v59AwA8E8PtpUHsSBbhs4LBhdqsEezHCIFyZGgWVYECJT4ELsmDYTWadCyigiMh5QWFEWqeyCAr3Q1VpxwxPNRWE2vhaBrCroJQ1+0r5aHcxlVVluPBhs8HFmUl8c3AMI6s2Ydv/iRnWfLV1MG0JQtjAX0O9PhraZyirfsm03Xx/pwtwQHyPunVoeg7RrCAR6ZIURiRiVPsCfLgjjzfWH2T9vsLQ8aE9XHx/dAzf6Gcn3lZvjKctQahwL3yxCAp2grcEMBqGF1+N+bPBZO629MgE15KpHTZzuCC6OyQNqJsVlDzAfGhWj4iEOYURiUh7Ckp5Y/0h3t50hDKv2bsR73by7Yt7M+OyvgzoGd/2F6/fewSnh5fyAjj8GeR/Caf2g7fU7JExjHMIKsE6Hod5jd0O7m7QcyiMvgN6DYWY7hrSEZGwpDAiEa20qoZlm4/y538cYO/x8tDxKwb1YGZWX64emoKjI1aNramCylNQWWyGl2APS+OgUlUCNRW0bTZU7QRnuwOiPRCXDGljYMBk9ZqISFhRGBHB3ANp3Z6T/PmfB/j7V/kE/689o3sM0y/tw+3jM+kWG9U5jWkcVBrPhqquqO1JqaHlM5vqbSVgd5pDOj2HwfCbIXkgRMcrnIiIZRRGRBo5XFjBG58e5K0NhymqMJezj3LauWZoCjdflM6kISlEuzp5gnDj2VCFB2HXB+Y6NKX55vNGcCZSS4Z2guHEDk4nRCVA8iAY+R3om6VCWBHpVAojIk2oqvGzfMsxXvvHAXbkloSOJ0Q7mTYylWkj08gakNz5wSSo/jo05QWw/2PY/wkUHwNfBXVFrn6an9UTXBrOAJsLYpPqCmEHXAN9L1OPiYh0KIURkbMwDIOvckv5yxdHWb7lGLnFVaHn4qIcXDWkJ1OG92LykF54Yl3WNbR+OMn7Anb/zSykLck1d4AOLVhX/59y412I6hXC2gFHDKQOh2HfggvGaShHRDqEwohIKwQCBhsOFPLe1mP8fUcBeSV1wcRhtzGhbxJXD+3JpCEpDEqJx2brgOLXlqosgpJjUHwEdiw3994pyTVX4G20BWKd+mEluL2iATjNHhPVmYhIB1AYEWmjQMBg29FiVuzIZ8WOfHbmlzZ4vne3GK4a0pOrh6Rw2YBk4twW7jdZv9ek5AjsWw1HPoWio+Aro65HpP6CbY17UZqpM8m42Aw5zigFFBFpNYURkXZy8GQ5H39dwMc7j/PPfSep9tX1OEQ57Ezol8SkIWavyYCecdb2mjQe0tnxFzi51+xNCQQXcmuqxwQa1JngBFdU7cqwcZDUH4bdDOmj6nYrBi1dLyJNUhgR6QCV1X7W7zvJxzsL+HhnAYcLKxs8n5kUw6TBKVw9tCdZ/XsQE2Xx9n3BhdqKDpuFsHs/hlNHgOD0YWjYY1L/7yAHZkipMX93uc2VZ4PP1V+6vntfM6T4aoOPNgIUiWgKIyIdzDAM9p0we01W7zrOp/sKqfbX6zVx2rm0fzJX1/aa9OsRZ2Frqes12bPS3Bjx+C4zrDRYbK25XpPGGi1db3eay9djmHv7NLVTcX2qTxE5rymMiHSycq+Pf+w9yaqdBazaeZyjRQ17TfomxzJpSAqThvTk0v4WTh2Guh6Tk3ubqTOB09c2OVtAaazRTsX12e3g9kCPoXXFs0HqURE5LyiMiFjIMAx2F5SxamcBH399nA0HCvEF6v6pRbvsZPVP5uqhKUwanMIFybHWNfZMdSaF+8wVYQM14KuibvpwS1aFbSkbZvFsAHPX4tqwYrOB0w0JaZAxHvpdBXG91Isi0gUpjIiEkdKqGtbtqes1qT91GMwZOpf2T+bS/klc2j+ZzCQLw0mw18TnNQtVj2wyh3VO7IbqstpOE5tZF9KqpetbI1inUhtUnA6ISoQeQ7RBoEgXojAiEqYMw+DrvFJW7TzOxzsL2HTwFP5Aw3+GvbvFMLZPd0ZleLiwt4cRvT3EWzmFuP6OxcHdiusvXV9WAAF/7ShOAHw1rdypuCXDP7XDWo03CBx2k1aTFQlTCiMiXUS518fGg6dYv+8kn+47ydYjxQ2GdMAcuRjQM54hqQkM7BnPgJR4BvaMp3/POGtrT+oP8dgwQ0p1pblbceOdir2lZm8LgGGA0XiX4qbCypmCirPumCsO0kbC6O9B79HgyVAwEQkTCiMiXVS518emg6fYeqSIrUeK2Xa0uMFS9fXZbJDuiaF39xgyuseQ0S2GjO6x9O4eQ+9uMaR1i8bttDis1N+p2F9d91x5QV3xbElubW0KjYKKjYazfeqrN5MHmzmTJyENuvU110MZOEnDOCIWUxgROY8cL/Xy5dFidheUsregnD3Hy9hTUEZxZU2z19ls0CPeTWpiNL0So0nzRJPqqfu9V6L5t2VDQPV7VuoHlcI95lL3BduhoqhebUpTGwQGZ+rUFsTaHOBJN6cUD7sJPL3VYyJiAYURkfOcYRicLK/m4Mlyjpyq5MipSo4W1f48VcHRokqqalpWVJrgdtLLUy+g1IaU0E9PNEmxUdgbT8/tSPXrVI5/1YINAuu3LVir4qidJpwKngvUYyLSyRRGRCJcMKzkFVeZj5K6n/klVeQWV5FfXEWp19ei13M5bKQk1AaW2qBSv3cl2PsS5bSf/cXaovEGgQfXmmujGD6a7jHBfM7mMAtqgz0mF2SZK8amjVIwEelACiMi0iJlXh95xWZAaRxagj9PlHlp6f9SJMdF1YUTTzQ94t10j3XRPTaKbrEuusVG0b32Z2K0s217+dRfTXbnX83F20qOmjN6QoHkTEM59XYudrohMQMuuERrmYh0EIUREWk3Nf4ABaXeBr0s9XtX8krMR/1NBFvCYbfhiXHhiXGRGO0kMcZFYrSLxBhn7c96x2uf89R7zu20Y6sqNntMjm2GzQshdyvUlNe+w9kKYO3U7VjcaC2T5P51e+yAgopIGyiMiEinMgyDUxU1tWGlkrxiL3klVRSWezlVUUNRRTWnys2fRZU1VFQ3FRJaLsphDwWXhBgXvaKqGB34ivEV67igaiexvkKi/SXY8dO4/8WODdsZNwXE3AjQ6TL32IGzL10f3MW4fnhp6nhLrnE4VXAr5wWFEREJa1U1foorazhVUU1JpY+SyhpKqmpqf9b/22f+rP97ZQ2BZv6Xy001PTlFImUM4jC3OlYzzH6QOKpCE4IdEMoe0LD89cwDR+baJo7GS9fbHWZ4MQJ14aWp4y25xhEFscmQmAmDr4eUoa0PPAo1EiYURkTkvGUYBuXV/oaBpZkw46s4RWrRF4yv/ITh/q/pThFxeLHVzsixY4QGbCC0CH1oTk79nXkAcNT9bnB6kGl87EzsTZ4TLLit7TlyuFoXeJoLNaCgIp1KYUREpLHawteqU8coP7oN59fLcZ/cid1bhB0fBgYB7NgINAgfBoRKT+oHk8ahwzjL8ca/h0aFaDhBuakg07rl6xqFGruzZb0vCinSjhRGRETO5kxrmZw6COXHzR2L/b7QirANK1yaL9Q1Gv3eIJicoVSmuZ15mgov9a9t9HQLNNH70lRIiekG0Ynmc5oGLa2gMCIi0lrBtUz8NXV77EDTS9e3oWbEsNkI2B2AjUDAgIAPjLqtBO34sBFoGEIaBZimwovNcfo5rQsrZwgpDhe4E6FbP+g1DDIvMadBO6MgoZfWaZFmKYyIiLSnppauDw5v1A8vzR0/2zWFe+GLRVCwE7wlgIFhdxCoDS+G3wcYBLDhrx1QChDA6TdqQ4wNHxBV+wzUm8TcqGelflBpuoYluAx/sPQXcyjH6Tb3AsoYX7dOS/3PZvghprvCSoTrkDCSnZ1NdnY2Bw4cAGDEiBH8/Oc/Z9q0aU1e4/V6+cUvfsEbb7xBXl4eGRkZPP7449x7773t/mFERM4L9YePoEWBp9RbwfEiL/llXgrK/FQf38PQgv+jV/UBYvwVOJoJKqFw4qjrTWkcTE7vVXFQF1BqZxgFe4TsTnDGQGx3SBsDAyaDJ/PMM3+0fst5rUPCyHvvvYfD4WDgQHOO/Z///Geef/55Nm/ezIgRI854zc0330x+fj5PP/00AwcOpKCgAJ/Px2WXXdbuH0ZEROqpDTWBylIKq6o5XuIlryJAVf4u+h9ZTnLlfmKoBH+g2ZACNKhXOfvQj43QpoXBjQudUacPW9ntzS80p16WLq/ThmmSkpJ4/vnnue+++057Licnh+9+97vs27ePpKSkNr+HwoiISDur1/tSbRgcrzAoProTz64lxJfsxe0vw8DAwI/d30RAgVAasTf8sxWaWGiuNb0sql8JWx0eRvx+P0uWLGHmzJls3ryZ4cOHn3bOj370I3bt2sW4ceP4n//5H+Li4rjpppv45S9/SUxMTJOv7fV68Xq9DT5MZmamwoiISEdrPERUVUTVrpUE9q/BXp4P/upQxQp+c1AntE6Lo67npOkalNY6Sy9LwG/OAopPMTdAHHaTWXALCihhoKVhxNnaF962bRtZWVlUVVURHx/PsmXLzhhEAPbt28fatWuJjo5m2bJlnDhxgh/96EcUFhby6quvNvke8+fP56mnnmpt00RE5FzFdIPeYxocih489bTiXaNwD+Wbl2Kc2IGruhQ/AQJ+8GHDQQCnw48TsGHD0eTE5ZYwgJp6f/qgpv5O0zazTYXlUHgQti4xe1hwQHwPM6CMvhMyxiqUhLFW94xUV1dz6NAhioqKeOedd3j55ZdZvXr1GQPJ1KlT+eSTT8jLy8Pj8QCwdOlSvvOd71BeXt5k74h6RkREuoB6vSg+w2B7bglr95Xy+d4DXO37lAmOr0lxFBNNAJfLgaN+zUgggDlnOdiPUtvz0S6Cr2mYPSmxyTB0Glx0F6ReqFDSiTqtZuTaa69lwIABvPjii6c9N3PmTNatW8eePXtCx7766iuGDx/Orl27GDRoUIveQzUjIiJdR5nXx7sb9rJs9UZ8pXm48JGe5OFfrxnByGRX3YygphaaA7NmxGY3h2ECPurWvm3pV1b9lVbs5lBOUh8YfAOMuRN6DDzL9dIeOmyYpjHDMBr0YtR3+eWXs2TJEsrKyoiPjwdg165d2O12MjIyzvWtRUQkDMW7nXxv4hBuvXQgb204zH+t2MXGwhqWLynh3sv78bNpQ3A7HTDgShh9+5kXmgOzNsRbAjuWw+H15lCRv7puQbmAH3xeGk5IDq4Q1yi0BKqh8ABsfgPyt8O1P4fUM88Clc7Xqp6RefPmMW3aNDIzMyktLWXx4sU8++yz5OTkMGXKFObOncvRo0d5/fXXASgrK2PYsGFceumlPPXUU5w4cYLvf//7XHXVVbz00kstbqR6RkREuq7C8mrmv/8VSzYdAWBYWiIv3DmGAT3jW/YCjRecq7/uSskR2PMR5G2B8hNQXWEGFKOGhhORg6vi1s7QuWACTPmlAkkH65Cekfz8fKZPn05ubi4ej4dRo0aFgghAbm4uhw4dCp0fHx/PihUrmDNnDuPGjSM5OZnbbruNp59+uo0fS0REupqkuCiev3U03xiZyqNvb+Wr3BJueWEdv7tjDFcPTTn7C7iiIamv+TiTYd80w0pFoRk4Cg/C13+Fg2ugvLB2mIe6HQ99fji8AVb+EqY+rSGbMKDl4EVEpNMUlFQx+83NfHagEJsNfvaNofzgyv7YbO0zETikpgqObYYvl8LelVB0xJyJg722HsUBMR4Ycxdc+aiKWjtIS7+/7U0+IyIi0s5SEqN54/uXcMeECzAMePaDr/nxW1uoqjnDdsbnwhUNfbJg6i/hphfggkvNxdOChbF2O9RUwO4VZg2JWEphREREOlWU086//8tIfnnzCBx2G+9uOcadL63nRNmZJ0OcE1c09M2Cb/y7WSfijDYDiWGYoaTsBHz9vtmTIpZRGBERkU5ns9mYntWX1++dQGK0k88PFXHzC+vYmVfaMW+YOsIsWM0YBw6XOdXX5jD3vdm/Sr0jFlMYERERy1w+sAfLZl1O3+RYjhZV8u3sf/Dx1wUd82apI+CKn0D3fuBy186ucUFprnpHLKYwIiIilhrQM55lP7qcS/snUeb1cd+fN7Bg3X46ZH5FxlgYeI05vdcRBdRAwIC9H8HRze3/ftIiCiMiImK57nFRvH7vJdw2LoOAAU+9t4N5y77smMLWIdMgPtX82xED9trekc9eghN7mr9eOoTCiIiIhIUop51ffXsU864fis0Giz47xC1/WMfu/HauI0m90OwdiUkEDPBVgMMNJ3bCtiUarrGAwoiIiIQNm83Gv145gNfumUCP+Ci+zivlxhfW8tq6/fgD7TRs44qGUbdBj6EQ8JqrudpsUFUEOz/QcI0FFEZERCTsXDW4J+8/eAVXDOpBVU2AJ9/bwb/8cR1fHi1unzfoMRDG3QcJqeYeN9XlUFUORQdg3W81XNPJFEZERCQspSRE8+d7JvDLW0aSEO1k65FibnphLT/53y84dLLi3N8gYyz0n2z+7vdjbrJnh9yt8MUiDdd0IoUREREJW3a7jemX9uGjn1zFjaPTCRjwzudHmPzrVcxduvXc6kmCxayxPYCAWTdic4KvXCuzdjKFERERCXspCdH8/o4xvDvrcq4Y1ANfwGDRZ4eZ8l9ruOvl9Xy4PY8af6D1L5x6IQy8FqJizb1r/F7AppVZO5k2yhMRkS7ns/2FvLJ2Hyt25BOsa02Ki+KGC9O4+aJ0Lr6gO3Z7CzffO7EHch6Hw+sBm7kYmt0J3TJg2nPmcI60SUu/vxVGRESkyzpyqoL/WX+QdzYd4URZdeh4724xXDsshSnDU7mkfxIux1kGAvavg7/Ng7K82n1rogAfDP4GXP2YWegqraYwIiIiEcPnD/CPvSd5d8tR/vZlHuXVdYulJUQ7uXpIClOG9+KqIT1JjHad/gI1VbD6Odj6FmAzZ9f4KsHphowJcNmD0H9i532g84TCiIiIRKSqGj+f7D7Bih15fPRVASfL63pMXA4bl/ZP5roRqUwd3ouUxOi6Cw9vgA8eg5JcqCkFA3Mjvag4SBoA3/w19BrR+R+oC1MYERGRiOcPGGw5fIoPd+SzYkc++46Xh56z2WBMZjeuG5HKdSNS6etxwCf/CTuWQtlJCNSYJ9miwGGDC2+Da35uzsKRFlEYERERaWTv8TI+3J5PzvY8vjhc1OC5Ib0S+O6Aam4peZNuueuweUvMYlZ/NRg2iOkOkx6Fi+5SIGkhhREREZFm5BVX8eGOPP62PY/1+wpDy833JZcH4/7GlMB6Yowq7HYbNrsdfNXgToCR/wJZs81VXKVZCiMiIiItVFRRzd+/KiDny1zW7DqB21/CI/a3mOL4jDhHDdH4cNjsOAhgc8ZCxsVw+UPQJ0u9JM1o6fe3sxPbJCIiEpa6xUbxnbEZfGdsBiVVNaz8qoDNmyD/4H4G+I9Qgw0/AXCAs6Yc+4F/4Di+C9vAa+DSH0KqClvPhXpGREREmlDx5f/h+79HiaoswO/3YwPsGIANHGBgpzw2k/JxD5Jy6W3ExMZZ3eSwomEaERGR9rDzA/joafynDhCoqcAP4Kc2kphfoaXE8KF/LJtjriDQczDRSb3xJCbSI95Nj/go4qKcxEY5iHWbP2NcDtwuO067HafDhtNuM3+321q+cmwXoDAiIiLSXvK3w6bXYcdyKD+J3/ADAQwC+LFhYIDfTjlu8o0kNgf68ffAJWxkCKXEt+qt7DZCIcVht+Fy2HHYbUQ57EQ57XU/6/0e7bLTPTaK5PgokuLc9E2OZUS6h1SPtfUsCiMiIiLtqaYKduXAqufg5B4IeDH3mzW/Rv2AUdtXYgOqiOFwVH/+7r6WnfYBHPIncdIXQ2WNn4pqH9W+QGhfnY4yPC2Rf72yPzeNTrekx0VhREREpCPkb4d//AG+/itUl5lrkBCofYAZUAKYkcQGDhckpEG3vjD4ekgZCtHxkDyAgNuDL2DgCwTMn/7a3/0G/oBBjT9Q+9P8vdofwFsToNrvp9oXwOsLUO0LUOULcKq8msLyao6XetlTUMbugtJQ2LlmaAp/uOtiol2OTv1PpTAiIiLSUWqq4Kv3YO1/wqkj5j42ho9gL0mdYG+EE7PvBHBGQVQi9BgCo++A5P7m8YRekJjeblOFC8urWfTZIX730W68vgC3js3g+VtHt8trt5TCiIiISEfL3w5f/C8c/hQKvgJvKWbosFEXTOr/HlTbQ2Gzg8tt/h3fw9yUb/A3oOcQiE6E2ORzDifr9pzge698imHA2/dnMa5v0jm9XmsojIiIiHSGmiqoPAUFO+HLt+HgWig6WttTYidYTdI0Gw2GdhwuiO5ubs7XaxhkXgKezHPqOXnsna0s3nCYKcN78dKMcW35lG2iMCIiItLZaqqgNA/2rISdf4WTe6HkKASCgaQlX7m1tSY4gIC5WZ8rBhJSoe9E6D/JDCfJAyCmW4uatTu/lCn/tQaXw8bGx6fgiXW18QO2jlZgFRER6WyuaEjqCxPuhQu/BSXH4Nhm+GKR2XNSWQRGDQ17QwKNXiQYWgJ1f1aXwcn9cPIAfLEY3N0gqT+M/A4MnHTWHpNBvRIY3CueXfllrNpVwM0X9W7fz32OFEZEREQ6Qkw389FrOAy9AQr3Q1UpHP8Kdv8NTh2EsgKoqawXUIIzcxozAJ/5q88PvhNQngdHNkC3DOh/FQz9JmSMa7K3ZOLAnuzKL2PjgVMKIyIiIhEnphv0HmP+PuBKGH272Wvir4HCg7DrAzi6AUrzzaEeozZ4NBtMDDPEnDoIm9+Erz8wZ+gMv+WMvSXj+nbn1XX72XCgsCM/aZsojIiIiHS2YK8JQPpoGDLVrDUpyYPyAtj/Mez/BEpyzZ6T0FJq9YZvgNCQTqAGyk9C+Vo49E/wpJu1JRP+NbSJ3+hM8/32FJRR4w/gctg757O2gMKIiIiI1YK1Jkl9zb8H1wsnJUfg8GeQ/6VZEFt+ArNnpPGwTm1gMaqh6BB88ZZ53cUzIW006QkpdIvyU1Tt4ODJCgamtG6Z+o7UqliUnZ3NqFGjSExMJDExkaysLD744IMWXbtu3TqcTicXXXRRW9opIiISOYLhpO+lMOo7MPUX8J1XYMZf4PrnITMLohIwv8Yd1C2uVq/XxF8Dx7+GD/8N3rwd28Jbeda1gCl8ytG9O8zhoDDRqqm97733Hg6Hg4EDBwLw5z//meeff57NmzczYsSIJq8rLi7m4osvZuDAgeTn57Nly5ZWNVJTe0VEROqpLIJD680pxLlbIO9LcxXY0NTh+l/twR4UgxoMqvwuvHFp9BhxTZumCbdGp60zkpSUxPPPP899993X5Dnf/e53GTRoEA6Hg3fffVdhREREpD0EF1zbvxY2LYDcrVBTXvtk4yEcc/k1nx9sDjtROMDpNKcJ9xwMlz8Egya3a/Na+v3d5uoVv9/P4sWLKS8vJysrq8nzFixYwN69e3niiSda/Nper5eSkpIGDxEREWnEFQ2JaTD6VrjjTfjOyzDyVkhIB7uT0LLzjdgIADXg80JloVlb8sEjsOejTm1+UKsLWLdt20ZWVhZVVVXEx8ezbNkyhg8ffsZzd+/ezWOPPcYnn3yC09nyt5o/fz5PPfVUa5smIiISuWK6wZBvmEMvxzbD1rfgwFo4dQgC3tqTGu+TY6tdf81hroGyeSH0HtshQzbNaXUYGTJkCFu2bKGoqIh33nmHmTNnsnr16tMCid/v58477+Spp55i8ODBrXqPuXPn8vDDD4f+LikpITMzs7VNFRERiTyuaOiTBeljzBk5x7aa65gc+ZRA0VEMfxV+bLiCwcTAXHI+pjv4qqDiZKeHkXOuGbn22msZMGAAL774YoPjRUVFdO/eHYejrosoEAhgGAYOh4MPP/yQyZNbNjalmhEREZFzULtnTtGBTSx5ezHD7Qe51FOMo7LQDCIJaeDJgB6D4Non2y2MdNreNIZh4PV6TzuemJjItm3bGhz74x//yMqVK3n77bfp16/fub61iIiItETtVOF4zwX8x//a8QTKyPnOaJKOb4R9K8EIQHwKjL6z03tFoJVhZN68eUybNo3MzExKS0tZvHgxq1atIicnBzCHV44ePcrrr7+O3W5n5MiRDa5PSUkhOjr6tOMiIiLS8ZwOO+7oWAqqoiiM60dSv1Hmhn4VJyE22ZIgAq0MI/n5+UyfPp3c3Fw8Hg+jRo0iJyeHKVOmAJCbm8uhQ4c6pKEiIiJy7rrHRVFS5eNURY15oP7S9BY555qRzqCaERERkfZxyx/WseVwES9OH8t1I1I79L06fJ0RERER6Xq6x7oAKKqotrgldRRGREREIkj32CiAumGaMKAwIiIiEkG6hcKIekZERETEAvHR5tyVsiqfxS2pozAiIiISQRLcZhgp9yqMiIiIiAXiasNImddvcUvqKIyIiIhEkDi3uU1LmVcFrCIiImKBhOjgMI16RkRERMQCcVGqGREREREL1dWMKIyIiIiIBYLDNAojIiIiYolgz0hFtZ9AIDy2p1MYERERiSDxtWEEoKw6PHpHFEZEREQiiNtpx24zf6+qDo8ZNQojIiIiEcRmsxHtMtcaqaoJWNwak8KIiIhIhAmGkcoa9YyIiIiIBaKd5td/lcKIiIiIWKFumEZhRERERCzgDoYRn2pGRERExAIxLvPrv1KzaURERMQKwWEar09hRERERCygmhERERGxVLQrOJtGNSMiIiJigWin1hkRERERC0VHaZhGRERELBTsGdEwjYiIiFiirmZEPSMiIiJiAU3tFREREUu5nZpNIyIiIhaKqg0j1X6FEREREbGAy2F+/ddobxoRERGxQrBnpEY9IyIiImKFKIeGaURERMRCoZ4Rn2FxS0wKIyIiIhEmWDPi7Yo9I9nZ2YwaNYrExEQSExPJysrigw8+aPL8pUuXMmXKFHr27Bk6/29/+9s5N1pERETarq5npAuGkYyMDJ599lk2btzIxo0bmTx5MjfffDPbt28/4/lr1qxhypQpvP/++2zatImrr76aG2+8kc2bN7dL40VERKT1XA4bED4FrDbDMM5pwCgpKYnnn3+e++67r0Xnjxgxgttvv52f//znLX6PkpISPB4PxcXFJCYmtrWpIiIiAmw6WMi3s/9Jn+RYVj96dYe9T0u/v51tfQO/38+SJUsoLy8nKyurRdcEAgFKS0tJSkpq9jyv14vX6w39XVJS0tZmioiISCNdfp2Rbdu2ER8fj9vt5v7772fZsmUMHz68Rdf++te/pry8nNtuu63Z8+bPn4/H4wk9MjMzW9tMERERaYKrq0/tHTJkCFu2bGH9+vX88Ic/ZObMmezYseOs1y1atIgnn3ySt956i5SUlGbPnTt3LsXFxaHH4cOHW9tMERERaUJoOfgw6Rlp9TBNVFQUAwcOBGDcuHFs2LCB3/72t7z44otNXvPWW29x3333sWTJEq699tqzvofb7cbtdre2aSIiItICwUXPavznyTojhmE0qO9obNGiRdx99928+eab3HDDDef6diIiInKOwm2jvFb1jMybN49p06aRmZlJaWkpixcvZtWqVeTk5ADm8MrRo0d5/fXXATOIzJgxg9/+9rdceuml5OXlARATE4PH42nnjyIiIiItEawZ8QcM/AEDh91maXta1TOSn5/P9OnTGTJkCNdccw2ffvopOTk5TJkyBYDc3FwOHToUOv/FF1/E5/Mxa9Ys0tLSQo8HH3ywfT+FiIiItFiwZwTCY62RVvWMvPLKK80+/9prrzX4e9WqVa1tj4iIiHSw4KJnYA7VRLscFrZGe9OIiIhEHJe9Xs9IGMyoURgRERGJMHa7LdQ7Eg5FrAojIiIiEahuFVbrp/cqjIiIiESguum9fotbojAiIiISkZz28Fn4TGFEREQkAgVrRvwBhRERERGxQHChs3BYZ0RhREREJALVX4XVagojIiIiEaiuZ0RhRERERCzgrA0jvoCGaURERMQCTkcwjKhnRERERCwQnNrr0zCNiIiIWCE4TOPXMI2IiIhYIThMowJWERERsYSm9oqIiIiltOiZiIiIWCpYwKqeEREREbFEsIC1RmFERERErBAsYPVrmEZERESsULcCq3pGRERExALO2tk0mtorIiIiltCiZyIiImIpLXomIiIiltLUXhEREbFU3dReDdOIiIiIBYIFrH4N04iIiIgVNLVXRERELBUsYPVpmEZERESsEOoZ0TCNiIiIWCFYM6JhGhEREbFEXc+IhmlERETEAtq1V0RERCzl0NReERERsZLLrtk0IiIiYiFHV11nJDs7m1GjRpGYmEhiYiJZWVl88MEHzV6zevVqxo4dS3R0NP379+dPf/rTOTVYREREzl1wnZEutzdNRkYGzz77LBs3bmTjxo1MnjyZm2++me3bt5/x/P3793P99ddzxRVXsHnzZubNm8cDDzzAO++80y6NFxERkbax28InjDhbc/KNN97Y4O9nnnmG7Oxs1q9fz4gRI047/09/+hMXXHABv/nNbwAYNmwYGzdu5D/+4z/49re/3fZWi4iIyDkJDtOEQxhpc82I3+9n8eLFlJeXk5WVdcZz/vnPfzJ16tQGx6677jo2btxITU1Nk6/t9XopKSlp8BAREZH2E5zaGzC6YBjZtm0b8fHxuN1u7r//fpYtW8bw4cPPeG5eXh69evVqcKxXr174fD5OnDjR5HvMnz8fj8cTemRmZra2mSIiItKM4DBNlytgBRgyZAhbtmxh/fr1/PCHP2TmzJns2LGjyfNttR82yKhNYI2P1zd37lyKi4tDj8OHD7e2mSIiItKM4DBNIAzCSKtqRgCioqIYOHAgAOPGjWPDhg389re/5cUXXzzt3NTUVPLy8hocKygowOl0kpyc3OR7uN1u3G53a5smIiIiLdRlp/aeiWEYeL3eMz6XlZXFihUrGhz78MMPGTduHC6X61zfWkRERNqoyxawzps3j08++YQDBw6wbds2Hn/8cVatWsVdd90FmMMrM2bMCJ1///33c/DgQR5++GG++uorXn31VV555RUeeeSR9v0UIiIi0ioOW/gUsLZqmCY/P5/p06eTm5uLx+Nh1KhR5OTkMGXKFAByc3M5dOhQ6Px+/frx/vvv8+Mf/5g//OEPpKen87vf/U7TekVERCwWTsM0rQojr7zySrPPv/baa6cdu+qqq/j8889b1SgRERHpWOFUwKq9aURERCKQPVgzEgbDNAojIiIiESi46JnfrzAiIiIiFgjtTaOeEREREbFC3dReixuCwoiIiEhECg3TBKxPIwojIiIiEcjeVRc9ExERkfND3aJnFjcEhREREZGIVLfomYZpRERExAJ1i55Z3BAURkRERCKSU4ueiYiIiJXqF7AaFgcShREREZEIFCxgBeuLWBVGREREIpDDURdGrJ7eqzAiIiISger3jCiMiIiISKcLzqYB64tYFUZEREQiUIMwop4RERER6WwaphERERFL2dUzIiIiIlYLLnwWUM2IiIiIWMEe2p9GYUREREQsENq5V2FERERErOCstyS8lRRGREREIpSGaURERMRSDhWwioiIiJUcGqYRERERKwULWBVGRERExBLqGRERERFLhcKIakZERETECuoZEREREUsFt6dRGBERERFLhKb2KoyIiIiIFew21YyIiIiIheoWPbO2HQojIiIiEcqujfJERETESvauuBz8/PnzGT9+PAkJCaSkpHDLLbewc+fOs163cOFCRo8eTWxsLGlpadxzzz2cPHmyzY0WERGRc9clZ9OsXr2aWbNmsX79elasWIHP52Pq1KmUl5c3ec3atWuZMWMG9913H9u3b2fJkiVs2LCB73//++fceBEREWm74HLwVveMOFtzck5OToO/FyxYQEpKCps2beLKK6884zXr16+nb9++PPDAAwD069ePH/zgBzz33HNtbLKIiIi0h1DNSFcuYC0uLgYgKSmpyXMuu+wyjhw5wvvvv49hGOTn5/P2229zww03NHmN1+ulpKSkwUNERETal702BXSpYZr6DMPg4YcfZuLEiYwcObLJ8y677DIWLlzI7bffTlRUFKmpqXTr1o3f//73TV4zf/58PB5P6JGZmdnWZoqIiEgT7GEyTNPmMDJ79my2bt3KokWLmj1vx44dPPDAA/z85z9n06ZN5OTksH//fu6///4mr5k7dy7FxcWhx+HDh9vaTBEREWmCI0xm07SqZiRozpw5LF++nDVr1pCRkdHsufPnz+fyyy/n0UcfBWDUqFHExcVxxRVX8PTTT5OWlnbaNW63G7fb3ZamiYiISAvVrTNicTtac7JhGMyePZulS5eycuVK+vXrd9ZrKioqsNsbvo3D4Qi9noiIiFgjNLW3Kw3TzJo1izfeeIM333yThIQE8vLyyMvLo7KyMnTO3LlzmTFjRujvG2+8kaVLl5Kdnc2+fftYt24dDzzwABMmTCA9Pb39PomIiIi0SrhslNeqYZrs7GwAJk2a1OD4ggULuPvuuwHIzc3l0KFDoefuvvtuSktLeeGFF/jJT35Ct27dmDx5Mr/61a/OreUiIiJyTmxhMrW3VWGkJcMqr7322mnH5syZw5w5c1rzViIiItLBHNq1V0RERKwULOm0uoZTYURERCRCBWfTdNlFz0RERKRrCxawKoyIiIiIJYI9I1avtKEwIiIiEqHsKmAVERERKwUXPbN6OXiFERERkQgVLoueKYyIiIhEqHBZ9ExhREREJEI5alOAZtOIiIiIJRyhnhGFEREREbGATWFERERErFS36Jm17VAYERERiVDBqb3am0ZEREQsYddy8CIiImIlrcAqIiIilnJobxoRERGxkoZpRERExFLam0ZEREQspUXPRERExFIaphERERFL2bVRnoiIiFgpVDOinhERERGxQnA5eNWMiIiIiCVsoUXPrG2HwoiIiEiEcmhqr4iIiFgpNEyjmhERERGxQmiYRmFERERErFBXwGptOxRGREREIpSWgxcRERFL2bUcvIiIiFjJrpoRERERsZIWPRMRERFL2UNTey1uh7VvLyIiIlYJFrD61TMiIiIiVnDU1owYCiMiIiJihS656Nn8+fMZP348CQkJpKSkcMstt7Bz586zXuf1enn88cfp06cPbrebAQMG8Oqrr7a50SIiInLuwmXRM2drTl69ejWzZs1i/Pjx+Hw+Hn/8caZOncqOHTuIi4tr8rrbbruN/Px8XnnlFQYOHEhBQQE+n++cGy8iIiJtFy6LnrUqjOTk5DT4e8GCBaSkpLBp0yauvPLKJq9ZvXo1+/btIykpCYC+ffu2rbUiIiLSboKzabrUME1jxcXFAKGQcSbLly9n3LhxPPfcc/Tu3ZvBgwfzyCOPUFlZ2eQ1Xq+XkpKSBg8RERFpX/ZQAau17WhVz0h9hmHw8MMPM3HiREaOHNnkefv27WPt2rVER0ezbNkyTpw4wY9+9CMKCwubrBuZP38+Tz31VFubJiIiIi0QLsM0be4ZmT17Nlu3bmXRokXNnhcIBLDZbCxcuJAJEyZw/fXX85//+Z+89tprTfaOzJ07l+Li4tDj8OHDbW2miIiINCFc9qZpU8/InDlzWL58OWvWrCEjI6PZc9PS0ujduzcejyd0bNiwYRiGwZEjRxg0aNBp17jdbtxud1uaJiIiIi1UF0YsbkdrTjYMg9mzZ7N06VJWrlxJv379znrN5ZdfzrFjxygrKwsd27VrF3a7/axBRkRERDpOlxymmTVrFm+88QZvvvkmCQkJ5OXlkZeX12C4Ze7cucyYMSP095133klycjL33HMPO3bsYM2aNTz66KPce++9xMTEtN8nERERkVYJzqaxuoC1VWEkOzub4uJiJk2aRFpaWujx1ltvhc7Jzc3l0KFDob/j4+NZsWIFRUVFjBs3jrvuuosbb7yR3/3ud+33KURERKTVwqVnpFU1Iy1Zu/6111477djQoUNZsWJFa95KREREOpgtTApYtTeNiIhIhAoVsAYsboe1by8iIiJWCQ7TaNdeERERsUSXnNorIiIi5w9bmBSwKoyIiIhEKIddPSMiIiJioXBZDl5hREREJEKFyzojCiMiIiIRKrTOiMXjNAojIiIiESo4TNOlloMXERGR84eGaURERMRSWmdERERELKV1RkRERMRSqhkRERERS2mdEREREbGUvTYFKIyIiIiIJVTAKiIiIpYKhhGwduEzhREREZEIZa/LIpYO1SiMiIiIRChb/Z4RC4dqFEZEREQilHpGRERExFL1a0asnFCjMCIiIhKhGhSwqmdEREREOptNwzQiIiJiJYddBawiIiJioYY1I+oZERERkU7WcDaNhe2w7q1FRETESvXXGfFrBVYRERGxQrB3RMM0IiIiYolw2CxPYURERCSC1YUR9YyIiIiIBYJlIwojIiIiYolgz4iWgxcRERFLBBc+U8+IiIiIWKJumMa6NiiMiIiIRLAuV8A6f/58xo8fT0JCAikpKdxyyy3s3LmzxdevW7cOp9PJRRdd1Np2ioiISAfocuuMrF69mlmzZrF+/XpWrFiBz+dj6tSplJeXn/Xa4uJiZsyYwTXXXNPmxoqIiEj7Cod1RpytOTknJ6fB3wsWLCAlJYVNmzZx5ZVXNnvtD37wA+68804cDgfvvvtuqxsqIiIi7S+4JHyXXQ6+uLgYgKSkpGbPW7BgAXv37uWJJ55o0et6vV5KSkoaPERERKT92bvyOiOGYfDwww8zceJERo4c2eR5u3fv5rHHHmPhwoU4nS3riJk/fz4ejyf0yMzMbGszRUREpBldep2R2bNns3XrVhYtWtTkOX6/nzvvvJOnnnqKwYMHt/i1586dS3Fxcehx+PDhtjZTREREmhEOPSOtqhkJmjNnDsuXL2fNmjVkZGQ0eV5paSkbN25k8+bNzJ49G4BAIIBhGDidTj788EMmT5582nVutxu3292WpomIiEgr2O1drIDVMAzmzJnDsmXLWLVqFf369Wv2/MTERLZt29bg2B//+EdWrlzJ22+/fdbrRUREpGOFwzojrQojs2bN4s033+Qvf/kLCQkJ5OXlAeDxeIiJiQHMIZajR4/y+uuvY7fbT6snSUlJITo6utk6ExEREekcXW6dkezsbIqLi5k0aRJpaWmhx1tvvRU6Jzc3l0OHDrV7Q0VERKT9hcM6IzbDyijUQiUlJXg8HoqLi0lMTLS6OSIiIueNa369ir3Hy3nrXy/lkv7J7fraLf3+1t40IiIiESwcekYURkRERCJY3TojXaRmRERERM4vtVkEv8KIiIiIWEHDNCIiImIph936dUYURkRERCJYl1tnRERERM4vtuAwTcC6NiiMiIiIRLBw2ChPYURERCSCqYBVRERELKV1RkRERMRSttAwjXVtUBgRERGJYHXDNNalEadl7ywiIiKWu2VMOmP7dGdgSrxlbVAYERERiWC3j7/A6iZomEZERESspTAiIiIillIYEREREUspjIiIiIilFEZERETEUgojIiIiYimFEREREbGUwoiIiIhYSmFERERELKUwIiIiIpZSGBERERFLKYyIiIiIpRRGRERExFJdYtdewzAAKCkpsbglIiIi0lLB7+3g93hTukQYKS0tBSAzM9PiloiIiEhrlZaW4vF4mnzeZpwtroSBQCDAsWPHSEhIwGaztdvrlpSUkJmZyeHDh0lMTGy315WOo3vW9eiedT26Z11PuN4zwzAoLS0lPT0du73pypAu0TNit9vJyMjosNdPTEwMq5snZ6d71vXonnU9umddTzjes+Z6RIJUwCoiIiKWUhgRERERS0V0GHG73TzxxBO43W6rmyItpHvW9eiedT26Z11PV79nXaKAVURERM5fEd0zIiIiItZTGBERERFLKYyIiIiIpRRGRERExFIRHUb++Mc/0q9fP6Kjoxk7diyffPKJ1U2KSPPnz2f8+PEkJCSQkpLCLbfcws6dOxucYxgGTz75JOnp6cTExDBp0iS2b9/e4Byv18ucOXPo0aMHcXFx3HTTTRw5cqQzP0pEmj9/PjabjYceeih0TPcr/Bw9epTvfe97JCcnExsby0UXXcSmTZtCz+uehRefz8e//du/0a9fP2JiYujfvz+/+MUvCAQCoXPOq3tmRKjFixcbLpfLeOmll4wdO3YYDz74oBEXF2ccPHjQ6qZFnOuuu85YsGCB8eWXXxpbtmwxbrjhBuOCCy4wysrKQuc8++yzRkJCgvHOO+8Y27ZtM26//XYjLS3NKCkpCZ1z//33G7179zZWrFhhfP7558bVV19tjB492vD5fFZ8rIjw2WefGX379jVGjRplPPjgg6Hjul/hpbCw0OjTp49x9913G59++qmxf/9+4+9//7uxZ8+e0Dm6Z+Hl6aefNpKTk42//vWvxv79+40lS5YY8fHxxm9+85vQOefTPYvYMDJhwgTj/vvvb3Bs6NChxmOPPWZRiySooKDAAIzVq1cbhmEYgUDASE1NNZ599tnQOVVVVYbH4zH+9Kc/GYZhGEVFRYbL5TIWL14cOufo0aOG3W43cnJyOvcDRIjS0lJj0KBBxooVK4yrrroqFEZ0v8LPz372M2PixIlNPq97Fn5uuOEG4957721w7Fvf+pbxve99zzCM8++eReQwTXV1NZs2bWLq1KkNjk+dOpV//OMfFrVKgoqLiwFISkoCYP/+/eTl5TW4X263m6uuuip0vzZt2kRNTU2Dc9LT0xk5cqTuaQeZNWsWN9xwA9dee22D47pf4Wf58uWMGzeOW2+9lZSUFMaMGcNLL70Uel73LPxMnDiRjz76iF27dgHwxRdfsHbtWq6//nrg/LtnXWKjvPZ24sQJ/H4/vXr1anC8V69e5OXlWdQqAXMM9OGHH2bixImMHDkSIHRPznS/Dh48GDonKiqK7t27n3aO7mn7W7x4MZ9//jkbNmw47Tndr/Czb98+srOzefjhh5k3bx6fffYZDzzwAG63mxkzZuiehaGf/exnFBcXM3ToUBwOB36/n2eeeYY77rgDOP/+nUVkGAmy2WwN/jYM47Rj0rlmz57N1q1bWbt27WnPteV+6Z62v8OHD/Pggw/y4YcfEh0d3eR5ul/hIxAIMG7cOP793/8dgDFjxrB9+3ays7OZMWNG6Dzds/Dx1ltv8cYbb/Dmm28yYsQItmzZwkMPPUR6ejozZ84MnXe+3LOIHKbp0aMHDofjtGRYUFBwWsqUzjNnzhyWL1/Oxx9/TEZGRuh4amoqQLP3KzU1lerqak6dOtXkOdI+Nm3aREFBAWPHjsXpdOJ0Olm9ejW/+93vcDqdof/eul/hIy0tjeHDhzc4NmzYMA4dOgTo31g4evTRR3nsscf47ne/y4UXXsj06dP58Y9/zPz584Hz755FZBiJiopi7NixrFixosHxFStWcNlll1nUqshlGAazZ89m6dKlrFy5kn79+jV4vl+/fqSmpja4X9XV1axevTp0v8aOHYvL5WpwTm5uLl9++aXuaTu75ppr2LZtG1u2bAk9xo0bx1133cWWLVvo37+/7leYufzyy0+bLr9r1y769OkD6N9YOKqoqMBub/gV7XA4QlN7z7t7ZlHhrOWCU3tfeeUVY8eOHcZDDz1kxMXFGQcOHLC6aRHnhz/8oeHxeIxVq1YZubm5oUdFRUXonGeffdbweDzG0qVLjW3bthl33HHHGaewZWRkGH//+9+Nzz//3Jg8eXJYTmE7H9WfTWMYul/h5rPPPjOcTqfxzDPPGLt37zYWLlxoxMbGGm+88UboHN2z8DJz5kyjd+/eoam9S5cuNXr06GH89Kc/DZ1zPt2ziA0jhmEYf/jDH4w+ffoYUVFRxsUXXxyaSiqdCzjjY8GCBaFzAoGA8cQTTxipqamG2+02rrzySmPbtm0NXqeystKYPXu2kZSUZMTExBjf/OY3jUOHDnXyp4lMjcOI7lf4ee+994yRI0cabrfbGDp0qPHf//3fDZ7XPQsvJSUlxoMPPmhccMEFRnR0tNG/f3/j8ccfN7xeb+ic8+me2QzDMKzsmREREZHIFpE1IyIiIhI+FEZERETEUgojIiIiYimFEREREbGUwoiIiIhYSmFERERELKUwIiIiIpZSGBERERFLKYyIiIiIpRRGRERExFIKIyIiImIphRERERGx1P8HTRuRErLfwloAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Fit\n",
    "\n",
    "S, I, R, _ = func.predict()\n",
    "\n",
    "plt.plot(S.detach())\n",
    "plt.plot(data, '.', alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 ('ODEAutodiff')",
   "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.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "9d9b1ab54ce76eb5c20cb46fa74b57f4d5689344267895880bd96eda7b19c259"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}