exp001-sir-constants-fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import torch as tc\n",
    "from torchdiffeq import odeint\n",
    "from torch import nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = 0\n",
    "tmax = 200 #days\n",
    "\n",
    "tspan = tc.arange(tmin, tmax).float()\n",
    "\n",
    "N = 10**6\n",
    "I0 = 1\n",
    "R0 = 0\n",
    "S0 = N - I0 - R0\n",
    "\n",
    "y0 = tc.tensor([S0, I0, R0, N]).float()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ODE(nn.Module):\n",
    "\n",
    "    def __init__(self, init_params, y0):\n",
    "        super(ODE, self).__init__()\n",
    "        self.init_params = init_params\n",
    "        self.y0 = y0\n",
    "        self.params = nn.Parameter(init_params)\n",
    "\n",
    "    def forward(self, t, y):\n",
    "        S, I, R, N = y\n",
    "\n",
    "        dS = - self.params[0] * S * I / N\n",
    "        dI = self.params[0] * S * I / N - self.params[1] * I\n",
    "        dR = self.params[1] * I\n",
    "        dN = tc.tensor(0.0)\n",
    "\n",
    "        return tc.stack([dS, dI, dR, dN])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_data(model, params, y0, tspan):\n",
    "    f = model(init_params=params, y0=y0)\n",
    "    with tc.no_grad():\n",
    "        S, I, R, N = odeint(f, y0, tspan).T\n",
    "        S = tc.distributions.Poisson(S).sample()\n",
    "        I = tc.distributions.Poisson(I).sample()\n",
    "        R = tc.distributions.Poisson(R).sample()\n",
    "\n",
    "        data = tc.distributions.Poisson(f.params[0] * S * I / N).sample()\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc068930c70>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcOElEQVR4nO3dfYxc1Znn8e9vDYaRA2MzbYzHxmNDN5bwaEVMibDqkCVDsA2axYSNskarsWeGjQc5KJMZr2ZhI41RZqUlmXGiQQpkTGJhr7K8ZGKGlpbEdqxM0FpLcNtxjA043RBebNnuNoQAToDYPPtHnTLXTVe73fVyb1X9PlKpbz11q+r07er73HOec28pIjAzs872b/JugJmZ5c/JwMzMnAzMzMzJwMzMcDIwMzPgrLwbMFFdXV0xd+7cvJthZtZSdu7ceTQipo+Mt2wymDt3Lv39/Xk3w8yspUh6ebS4h4nMzMzJwMzMnAzMzAwnAzMzw8nAzMxwMjAzM5wMzMyMcSQDSeslDUnam4k9Iml3ur0kaXeKz5X0m8xj38w850pJz0galHSvJKX4BZK2ShpIP6c14Pc0s1FsHzzKLfdtZ/vg0bybYjkbT8/gQWBJNhAR/ykiroiIK4DvAZsyD79QeSwibs/E7wc+B/SkW+U17wS2RUQPsC3dN7MmWLtlP7teeYOVG/udEDrcaZNBRDwJvD7aY+no/rPAQ2O9hqSZwPkR8VSUv01nI3BzengpsCEtb8jEzazBVi+az5TJkzj23gnW9O11L6GD1VozuAY4EhEDmdg8ST+V9GNJ16TYLOBAZp0DKQYwIyIOpeXDwIxqbyZppaR+Sf3Dw8M1Nt3Meru7WLe8xMI5UxFi1ytvsHbL/rybZTmo9dpEt3Jqr+AQMCciXpN0JfAvkhaM98UiIiRV/R7OiFgHrAMolUr+vk6zCdo+eJS1W/azetF8eru76O3uOiVmnWfCyUDSWcAtwJWVWES8C7yblndKegG4DDgIzM48fXaKARyRNDMiDqXhpKGJtsnMTm/74FFWbuzn2HsnWLtlP73dXQAnk4J1plqGiT4FPB8RJ4d/JE2XNCktX0K5UPxiGgZ6U9LVqc6wHHg8Pa0PWJGWV2TiZlZHlZlDd/ft49h7J5gyedKHegGeXdS5xjO19CHg/wHzJR2QdFt6aBkfLhx/AtiTppr+M3B7RFSKz6uAbwGDwAvA91P8HuB6SQOUE8w9E/91zKyaysyhIFg4Zyrrlpc+1BOorONicudReXJP6ymVSuHvMzAbn+2DR1nTtxch7r5pQdXhoErd4K13jjMw9DYL50xl06reJrfWGknSzogojYy37JfbmNn4rd2yn8GhYyycM3XMuoCLyZ3LycCsA6xeNP+Mdu4uJnceJwOzDuCdu52OL1RnZqPyzKLO4mRgZqPyzKLO4mRg1sZqObpfvWi+L1PRQVwzMGtjlaP77JnG4+WZRZ3FycCsjZ3pLKLRuPjcGZwMzNpYPXfkIy9uZ+3FNQMzG5fskJO1HycDMxuXSkHZtYP25GEiMxsX1w7am3sGZjZuPhGtfTkZmNm4uW7QvpwMzNpUI47iXTdoX64ZmLWpWk44q8Z1g/blnoFZm/JRvJ0JJwOzNlMZHgLYtKrXR/I2Lk4GZm3GRV6biNMmA0nrJQ1J2puJ3S3poKTd6XZj5rG7JA1K2i9pcSa+JMUGJd2Zic+T9JMUf0TS5Hr+gmadphnDQ55i2n7G0zN4EFgySvzrEXFFuj0BIOlyYBmwID3nPkmTJE0CvgHcAFwO3JrWBfhKeq1u4JfAbbX8QmadrFnXD3Lvo/2cNhlExJPA6+N8vaXAwxHxbkT8AhgErkq3wYh4MSLeAx4GlkoS8EfAP6fnbwBuPrNfwcwqmrWTdnG6/dRSM7hD0p40jDQtxWYBr2bWOZBi1eK/B7wREcdHxEclaaWkfkn9w8PDNTTdrD01ayfd293l4nSbmWgyuB+4FLgCOASsrVeDxhIR6yKiFBGl6dOnN+MtzVpKM3fSrhu0lwklg4g4EhEnIuJ94AHKw0AAB4GLM6vOTrFq8deAqZLOGhE3s4Jz3aC9TCgZSJqZuftpoDLTqA9YJukcSfOAHuBpYAfQk2YOTaZcZO6LiAB+BHwmPX8F8PhE2mRmzeW6QXs57eUoJD0EXAt0SToArAGulXQFEMBLwF8ARMQ+SY8CzwLHgc9HxIn0OncAm4FJwPqI2Jfe4r8BD0v6H8BPgW/X65cz6xR5fAuZL03RXlQ+OG89pVIp+vv7826GWSHcct92dr3yBgvnTGXTqt68m2MFJmlnRJRGxn0Gslkb8JCN1cpXLTVrA3kO2eQxRGX1556BmdXEs4rag5OBmdXEQ1TtwcNEZlYTzypqD+4ZmJmZk4GZmTkZmLU8XyPI6sHJwKzFeTaP1YOTgVmLK8JsHvdOWp9nE5m1uCLM5sn2TvJui02MewZmVrMi9E6sNk4GZlazypfqAB4ualFOBmZWNy5mty4nA7MWVcSirYeLWpcLyGYtqohF2yIUs21i3DMwa1E+Crd6cs/ArEX5KNzqyT0DMzM7fTKQtF7SkKS9mdjfS3pe0h5Jj0mamuJzJf1G0u50+2bmOVdKekbSoKR7JSnFL5C0VdJA+jmtAb+nmZmNYTw9gweBJSNiW4E/jIh/C/wcuCvz2AsRcUW63Z6J3w98DuhJt8pr3glsi4geYFu6b2ZVFHEWkbW+0yaDiHgSeH1EbEtEHE93nwJmj/UakmYC50fEUxERwEbg5vTwUmBDWt6QiZvZKIo+l9/JqjXVo2bw58D3M/fnSfqppB9LuibFZgEHMuscSDGAGRFxKC0fBmZUeyNJKyX1S+ofHh6uQ9PNWk/RZxEVPVnZ6GqaTSTpS8Bx4DspdAiYExGvSboS+BdJC8b7ehERkmKMx9cB6wBKpVLV9czaWdFnEa1eNJ+1W/YXNlnZ6CacDCT9KfDHwHVp6IeIeBd4Ny3vlPQCcBlwkFOHkmanGMARSTMj4lAaThqaaJvMLH9FT1Y2ugkNE0laAvwNcFNE/DoTny5pUlq+hHKh+MU0DPSmpKvTLKLlwOPpaX3AirS8IhM3M7MmOW3PQNJDwLVAl6QDwBrKs4fOAbamGaJPpZlDnwC+LOm3wPvA7RFRKT6vojwz6Xco1xgqdYZ7gEcl3Qa8DHy2Lr+ZmZmNm9IIT8splUrR39+fdzPMrIrtg0dP1g48bFQcknZGRGlk3Gcgm1lDeFZRa3EyMGshrTSHv+hTYO1UvlCdWQsp4mWrq/GsotbinoFZC/HRtjWKewZmLcRH29Yo7hmYmZmTgZk1TisVvDudk4FZC2jVnaqnl7YOJwOzFtCqO1UXvFuHC8hmLaBVrwTqgnfrcDIwawHeqVqjeZjIzMycDMys8Vq1AN5JnAzMrOEqBfCVG/udEArKycDMGm71ovlMmTyJY++daLkZUZ3CycDMGq63u4t1y0ueZlpgnk1kZk3hGVHF5p6BWcG5+GrN4GRgVnCtevaxtZZxJQNJ6yUNSdqbiV0gaaukgfRzWopL0r2SBiXtkbQw85wVaf0BSSsy8SslPZOec68k1fOXNGtlvqSDNcN4ewYPAktGxO4EtkVED7At3Qe4AehJt5XA/VBOHsAa4GPAVcCaSgJJ63wu87yR72XWsXq7u9i0qtfj7dZQ40oGEfEk8PqI8FJgQ1reANyciW+MsqeAqZJmAouBrRHxekT8EtgKLEmPnR8RT0VEABszr2VmZk1QS81gRkQcSsuHgRlpeRbwama9Ayk2VvzAKPEPkbRSUr+k/uHh4RqabmZmWXUpIKcj+qjHa53mfdZFRCkiStOnT2/025lZA3h2VDHVkgyOpCEe0s+hFD8IXJxZb3aKjRWfPUrcrKO1607Ts6OKqZZk0AdUZgStAB7PxJenWUVXA79Kw0mbgUWSpqXC8SJgc3rsTUlXp1lEyzOvZdax2nWn6dlRxTSuM5AlPQRcC3RJOkB5VtA9wKOSbgNeBj6bVn8CuBEYBH4N/BlARLwu6e+AHWm9L0dEpSi9ivKMpd8Bvp9uZh2tVb/Q5nR8JnIxqTzc33pKpVL09/fn3Qwzs5YiaWdElEbGfQayWQG1a73AisvJwKyA2rVeUOFkVzxOBmYF1O5F1nZPdq3Il7A2K6B2L7K2a3G8lTkZmFnTtXuya0UeJjKz3Lh2UBxOBmaWG9cOisPJwMxy0+6F8lbimoGZ5ca1g+Jwz8DMzJwMzIrEBVXLi5OBWYG4oGp5cTIwKxAXVC0vTgZmBbF98OjJs3I7qajqobFicDIwK4hOHSLq1N+7aJwMzAqiU4eIOvX3Lhp/uY2ZFUKnDpM1m7/cxswKzcNF+XIyMLNC8HBRviZ8OQpJ84FHMqFLgL8FpgKfA4ZT/L9HxBPpOXcBtwEngC9ExOYUXwL8IzAJ+FZE3DPRdplZa/KlKfI14WQQEfuBKwAkTQIOAo8BfwZ8PSL+Ibu+pMuBZcAC4PeBH0q6LD38DeB64ACwQ1JfRDw70baZmdmZqdeF6q4DXoiIlyVVW2cp8HBEvAv8QtIgcFV6bDAiXgSQ9HBa18nAzKxJ6lUzWAY8lLl/h6Q9ktZLmpZis4BXM+scSLFq8Q+RtFJSv6T+4eHh0VYxazk+6cqKoOZkIGkycBPw3RS6H7iU8hDSIWBtre9RERHrIqIUEaXp06fX62XNcuVZNKdycsxHPXoGNwC7IuIIQEQciYgTEfE+8AAfDAUdBC7OPG92ilWLm3UEz6I5lZNjPuqRDG4lM0QkaWbmsU8De9NyH7BM0jmS5gE9wNPADqBH0rzUy1iW1jVrez7R6sNWL5pP94VTeOud4+4dNFFNyUDSFMqzgDZlwl+V9IykPcAngb8CiIh9wKOUC8M/AD6fehDHgTuAzcBzwKNpXbO256PgD+vt7uL8c89mYOhtb5cm8uUozHLknsHovF0ap9rlKJwMzMw6iK9NZGYtxzOLmsfJwMwKyzWV5nEyMLPC8rTb5qnX5SjMzOrOF69rHvcMzMzMycDMzJwMzHLhWTJWNE4GZjnwLJnxc+JsDicDsxx4lsz4OXE2h2cTmeXAs2TGb/Wi+ScvTWGN42Rg1mS+7s6ZceJsDg8TmTWZhz0mxrWDxnIyMGsy1wsmxkm0sTxMZNZkHvaYGNcOGsvJwMxagpNoY3mYyMxahusGjeNkYNYk3pHVznWDxnEyMGsS78hqVym+L15wkRNrndWcDCS9JOkZSbsl9afYBZK2ShpIP6eluCTdK2lQ0h5JCzOvsyKtPyBpRa3tMisazyKqXW93F5tW9bJ532En1jqrV8/gkxFxReZ7Ne8EtkVED7At3Qe4AehJt5XA/VBOHsAa4GPAVcCaSgIxaxeVHZmLoLVbvWg+3RdO4a13jrt3UCeNGiZaCmxIyxuAmzPxjVH2FDBV0kxgMbA1Il6PiF8CW4ElDWqbWdO5XlBfvd1dnH/u2QwMve3eQZ3UIxkEsEXSTkkrU2xGRBxKy4eBGWl5FvBq5rkHUqxa/BSSVkrql9Q/PDxch6abNYfrBfXnYbf6qsd5Bh+PiIOSLgS2Sno++2BEhKSow/sQEeuAdQClUqkur2nWDD5hqv583kF91ZwMIuJg+jkk6THKY/5HJM2MiENpGGgorX4QuDjz9NkpdhC4dkT8X2ttm1lReMdlRVfTMJGkKZLOqywDi4C9QB9QmRG0Ang8LfcBy9OsoquBX6XhpM3AIknTUuF4UYqZmY3J9Zj6qLVmMAP4v5J+BjwN/J+I+AFwD3C9pAHgU+k+wBPAi8Ag8ACwCiAiXgf+DtiRbl9OMbOW5h1V47keUx+KaM2h91KpFP39/Xk3w2xMt9y3nV2vvMHCOVPZtKo37+a0JX8/xJmRtDNzGsBJvlCdWQO5cNx4rsfUh5OBWQN5R9Uc7h3UztcmMrOWV6kbrOnb6xrNBDkZmDWIi8fNUzkBTcjF5AnyMJFZg2RnuXjoorEqw3HZ4SI7M04GZg3i4nHzuUYzcR4mMquzyvAQ4KuU5sDDcxPjZGBWZz4JKl8uJk+Mk4FZnflqmvlyMXliXDMwqzOPW+fLxeSJcTIwqxOf+FQslb9BpWfgv8nYPExkVieuFRSP/ybj52RgVieuFRSP/ybj56uWmlnb8xDeB6pdtdQ9A7M68Nz2YvNw0ek5GZjVgXc2xbZ60Xy6L5zCW+8cd8KuwsnArA48Nl1svd1dnH/u2QwMve2EXYWTgVkNfOmJ1lFJ2IsXXOQhvVE4GZjVwMNDraO3u4tNq3rZvO+wL1cxigknA0kXS/qRpGcl7ZP0lyl+t6SDknan242Z59wlaVDSfkmLM/ElKTYo6c7afiWz5tg+eJQ33/ktPRd+xMNDLcSXqxhdLWcgHwdWR8QuSecBOyVtTY99PSL+IbuypMuBZcAC4PeBH0q6LD38DeB64ACwQ1JfRDxbQ9vMGm7tlv0MDh1j4ZypHh5qIb5cxegmnAwi4hBwKC2/Jek5YNYYT1kKPBwR7wK/kDQIXJUeG4yIFwEkPZzWdTKwQvP3FbQ2X67iVHWpGUiaC3wU+EkK3SFpj6T1kqal2Czg1czTDqRYtfho77NSUr+k/uHh4Xo03eyMuWjcPny56w/UnAwkfQT4HvDFiHgTuB+4FLiCcs9hba3vURER6yKiFBGl6dOn1+tlzc6Ii8btw/WDD9SUDCSdTTkRfCciNgFExJGIOBER7wMP8MFQ0EHg4szTZ6dYtbhZ4bho3F4qM4zuvmlBx5+UVstsIgHfBp6LiK9l4jMzq30a2JuW+4Blks6RNA/oAZ4GdgA9kuZJmky5yNw30XaZNVKlaHzeuWd5eKiNZE9K69Qho1p6Br3AnwB/NGIa6VclPSNpD/BJ4K8AImIf8CjlwvAPgM+nHsRx4A5gM/Ac8Gha16wwKnWCxQsu8pnGbWrkkNHKjf0dlRB81VKzcbjlvu3seuUNFs6ZyqZVvXk3xxpo++BRVm7s59h7J9ry7+2rlppNkOsEnaW3u4t1y0sd1wN0MjA7DdcJOk+lsAx0TP3AycDsNHxF0s5VmUa8cmM///TjF9o6MTgZmFXhk8ts9aL5TJk8iWPvneDebQNtfS6Ck4HZKCpFxHb+57fTy9YPvnBdz8keYjt+s10tF6ozazuVi5e99c5xjr13gimTJ3l4qMNVLmwH8Iezfvfk56PyRTnt0mN0z8AsozJGHAQL50xl3fJS2/yzW+2yn492O2PZPQOzJDuF9O6bFjgJ2Idkr1SbvcgdgFBLf26cDKyjbR88ypq+vQgB+PsJbEzZISPglCEjgDV9ezn/3LNZvWh+y32GPExkHa1yDsHA0Nsnh4ZcI7DxGHmRu54LP9LSl7Lw5Siso2V7Bq3cxbdiyF7KovvCKSd7CcDJ4aW8P2PVLkfhZGAdxwnAGik7I21g6G26L5zCoTfeKcy1jqolA9cMrGNk/0kHh44BtNXUQCuGkd+xnJ2mvHjBRXzqa/9ayAMRJwNrayMLxJUjte4LpyDk+oA1zMikUJmBVDkQKdosJCcDa2vZf77uC6ecLBDn/Y9nnWPkDKSRByeVGJQTw2eunM3mfYeb/jl1zcDaTrY38JkrZ/Pdna8W5ujLrGK0Xitw8lpIUyZP4gvX9dT98+sCsrWlShd88YKLTv7TwAf/WEUo2JmdzsgDmHu3DZxMCMfeOwGUe7ZQ+7CSC8jWNrIJoPJPs//wW6f807gmYK2kt7uLH/71tSfvV66BVO0gpxETH5wMrLBGO+rPHjVVEkCjutNmecnWGf7i318KnNp7aMRBTmGGiSQtAf4RmAR8KyLuGWt9DxO1vuzOfvO+w1V3+tmu8sjx1DwKbWatrNA1A0mTgJ8D1wMHgB3ArRHxbLXnOBnkIztNLjttbuSOfDzLI3f2Y+30s89zAjCbuKIng38H3B0Ri9P9uwAi4n9We85Ek0EtO6/xLI92hNtO7wPlccvKTnqso/fxLGeP8Ef+Pt7pm9Vf0ZPBZ4AlEfFf0v0/AT4WEXeMWG8lsBJgzpw5V7788stn/F633LedXa+8MeGd13iWG/naeb9P9tT6sY7ex5tovLM3a662SAZZ7hnk8z5337QA4JSxfu/QzVpH0ZNB04aJzMw6WbVkUJTvM9gB9EiaJ2kysAzoy7lNZmYdoxDnGUTEcUl3AJspTy1dHxH7cm6WmVnHKEQyAIiIJ4An8m6HmVknKsowkZmZ5cjJwMzMnAzMzMzJwMzMKMh5BhMhaRg481OQy7qAo3VsTr0UtV1Q3La5XWfG7TpzRW3bRNv1BxExfWSwZZNBLST1j3bSRd6K2i4obtvcrjPjdp25orat3u3yMJGZmTkZmJlZ5yaDdXk3oIqitguK2za368y4XWeuqG2ra7s6smZgZman6tSegZmZZTgZmJlZ5yUDSUsk7Zc0KOnOHNtxsaQfSXpW0j5Jf5nid0s6KGl3ut2YQ9tekvRMev/+FLtA0lZJA+nntCa3aX5mm+yW9KakL+a1vSStlzQkaW8mNuo2Utm96TO3R9LCJrfr7yU9n977MUlTU3yupN9ktt03m9yuqn87SXel7bVf0uImt+uRTJtekrQ7xZu5vartHxr3GYuIjrlRvjz2C8AlwGTgZ8DlObVlJrAwLZ8H/By4HLgb+K85b6eXgK4Rsa8Cd6blO4Gv5Px3PAz8QV7bC/gEsBDYe7ptBNwIfB8QcDXwkya3axFwVlr+SqZdc7Pr5bC9Rv3bpf+DnwHnAPPS/+ykZrVrxONrgb/NYXtV2z807DPWaT2Dq4DBiHgxIt4DHgaW5tGQiDgUEbvS8lvAc8CsPNoyTkuBDWl5A3Bzfk3hOuCFiJjoGeg1i4gngddHhKtto6XAxih7CpgqaWaz2hURWyLieLr7FDC7Ee99pu0aw1Lg4Yh4NyJ+AQxS/t9tarskCfgs8FAj3nssY+wfGvYZ67RkMAt4NXP/AAXYAUuaC3wU+EkK3ZG6euubPRyTBLBF0k5JK1NsRkQcSsuHgRk5tKtiGaf+g+a9vSqqbaMife7+nPIRZMU8ST+V9GNJ1+TQntH+dkXZXtcARyJiIBNr+vYasX9o2Ges05JB4Uj6CPA94IsR8SZwP3ApcAVwiHI3tdk+HhELgRuAz0v6RPbBKPdLc5mTrPLXot4EfDeFirC9PiTPbVSNpC8Bx4HvpNAhYE5EfBT4a+B/Szq/iU0q5N8u41ZOPeho+vYaZf9wUr0/Y52WDA4CF2fuz06xXEg6m/If+jsRsQkgIo5ExImIeB94gAZ1j8cSEQfTzyHgsdSGI5VuZ/o51Ox2JTcAuyLiSGpj7tsro9o2yv1zJ+lPgT8G/nPaiZCGYV5Lyzspj81f1qw2jfG3K8L2Ogu4BXikEmv29hpt/0ADP2Odlgx2AD2S5qUjzGVAXx4NSeOR3waei4ivZeLZcb5PA3tHPrfB7Zoi6bzKMuXi417K22lFWm0F8Hgz25VxytFa3ttrhGrbqA9YnmZ8XA38KtPVbzhJS4C/AW6KiF9n4tMlTUrLlwA9wItNbFe1v10fsEzSOZLmpXY93ax2JZ8Cno+IA5VAM7dXtf0DjfyMNaMyXqQb5ar7zyln9S/l2I6PU+7i7QF2p9uNwP8CnknxPmBmk9t1CeWZHD8D9lW2EfB7wDZgAPghcEEO22wK8Brwu5lYLtuLckI6BPyW8vjsbdW2EeUZHt9In7lngFKT2zVIeTy58jn7Zlr3P6a/8W5gF/Afmtyuqn874Etpe+0Hbmhmu1L8QeD2Ees2c3tV2z807DPmy1GYmVnHDROZmdkonAzMzMzJwMzMnAzMzAwnAzMzw8nAzMxwMjAzM+D/A6w+TA7BiGdgAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = simulate_data(ODE, tc.tensor([0.2, 0.1]), y0, tspan)\n",
    "plt.plot(data, '.', markersize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(func, data, y0, t):\n",
    "    S, I, R, _ = odeint(func, y0, t, method='euler').T\n",
    "\n",
    "    pred = func.params[0] * S * I / y0[-1]\n",
    "\n",
    "    return (1/y0[-1]) * (pred.mean() - (1/len(t)) * tc.dot(pred.log(), data))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "func = ODE(init_params = tc.tensor([0.5, 0.3]), y0=y0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss : -0.0020960418041795492 , Params : tensor([0.5000, 0.3000]), Grad : tensor([ 0.3176, -0.2695]) \n",
      "Loss : -0.007476808037608862 , Params : tensor([0.4809, 0.2973]), Grad : tensor([ 0.3228, -0.2817]) \n",
      "Loss : -0.02929786406457424 , Params : tensor([0.4316, 0.3209]), Grad : tensor([ 0.2106, -0.1960]) \n",
      "Loss : -0.024280589073896408 , Params : tensor([0.3959, 0.3374]), Grad : tensor([-0.3713,  0.3633]) \n",
      "Loss : -0.0308226365596056 , Params : tensor([0.4099, 0.3070]), Grad : tensor([ 0.1318, -0.1211]) \n",
      "Loss : -0.028831271454691887 , Params : tensor([0.4089, 0.2942]), Grad : tensor([ 0.2279, -0.2117]) \n",
      "Loss : -0.03132360428571701 , Params : tensor([0.3950, 0.2956]), Grad : tensor([ 0.0836, -0.0751]) \n",
      "Loss : -0.031207874417304993 , Params : tensor([0.3823, 0.2965]), Grad : tensor([-0.1059,  0.1060]) \n",
      "Loss : -0.0309747327119112 , Params : tensor([0.3756, 0.2918]), Grad : tensor([-0.1367,  0.1355]) \n",
      "Loss : -0.03136257454752922 , Params : tensor([0.3717, 0.2844]), Grad : tensor([-0.0910,  0.0916]) \n",
      "Loss : -0.03161302208900452 , Params : tensor([0.3678, 0.2770]), Grad : tensor([-0.0437,  0.0462]) \n",
      "Loss : -0.03170012682676315 , Params : tensor([0.3632, 0.2702]), Grad : tensor([-0.0167,  0.0204]) \n",
      "Loss : -0.031736645847558975 , Params : tensor([0.3579, 0.2639]), Grad : tensor([-0.0063,  0.0104]) \n",
      "Loss : -0.031763751059770584 , Params : tensor([0.3522, 0.2578]), Grad : tensor([-0.0050,  0.0091]) \n",
      "Loss : -0.03178839385509491 , Params : tensor([0.3463, 0.2518]), Grad : tensor([-0.0079,  0.0118]) \n",
      "Loss : -0.03181163966655731 , Params : tensor([0.3403, 0.2457]), Grad : tensor([-0.0116,  0.0152]) \n",
      "Loss : -0.03183569386601448 , Params : tensor([0.3342, 0.2395]), Grad : tensor([-0.0136,  0.0171]) \n",
      "Loss : -0.031863003969192505 , Params : tensor([0.3282, 0.2330]), Grad : tensor([-0.0124,  0.0158]) \n",
      "Loss : -0.031892772763967514 , Params : tensor([0.3222, 0.2264]), Grad : tensor([-0.0074,  0.0109]) \n",
      "Loss : -0.031920675188302994 , Params : tensor([0.3163, 0.2197]), Grad : tensor([-0.0004,  0.0042]) \n",
      "Loss : -0.031944602727890015 , Params : tensor([0.3104, 0.2130]), Grad : tensor([ 0.0053, -0.0014]) \n",
      "Loss : -0.03196736052632332 , Params : tensor([0.3044, 0.2064]), Grad : tensor([ 0.0069, -0.0030]) \n",
      "Loss : -0.03199036419391632 , Params : tensor([0.2982, 0.2001]), Grad : tensor([ 0.0042, -0.0006]) \n",
      "Loss : -0.03201162442564964 , Params : tensor([0.2920, 0.1939]), Grad : tensor([0.0004, 0.0029]) \n",
      "Loss : -0.03203127905726433 , Params : tensor([0.2861, 0.1876]), Grad : tensor([-0.0008,  0.0039]) \n",
      "Loss : -0.03204990550875664 , Params : tensor([0.2804, 0.1814]), Grad : tensor([0.0010, 0.0019]) \n",
      "Loss : -0.03206680715084076 , Params : tensor([0.2748, 0.1754]), Grad : tensor([0.0027, 0.0002]) \n",
      "Loss : -0.03208228945732117 , Params : tensor([0.2693, 0.1696]), Grad : tensor([0.0020, 0.0007]) \n",
      "Loss : -0.03209617733955383 , Params : tensor([0.2640, 0.1640]), Grad : tensor([0.0005, 0.0019]) \n",
      "Loss : -0.03210853040218353 , Params : tensor([0.2590, 0.1587]), Grad : tensor([0.0007, 0.0015]) \n",
      "Loss : -0.03211936354637146 , Params : tensor([0.2542, 0.1535]), Grad : tensor([0.0015, 0.0006]) \n",
      "Loss : -0.03212875500321388 , Params : tensor([0.2497, 0.1487]), Grad : tensor([0.0011, 0.0008]) \n",
      "Loss : -0.03213677182793617 , Params : tensor([0.2455, 0.1442]), Grad : tensor([0.0006, 0.0012]) \n",
      "Loss : -0.03214354068040848 , Params : tensor([0.2415, 0.1400]), Grad : tensor([0.0008, 0.0008]) \n",
      "Loss : -0.03214915469288826 , Params : tensor([0.2379, 0.1361]), Grad : tensor([0.0009, 0.0006]) \n",
      "Loss : -0.03215375915169716 , Params : tensor([0.2345, 0.1325]), Grad : tensor([0.0005, 0.0007]) \n",
      "Loss : -0.03215748071670532 , Params : tensor([0.2314, 0.1292]), Grad : tensor([0.0005, 0.0006]) \n",
      "Loss : -0.0321604385972023 , Params : tensor([0.2286, 0.1262]), Grad : tensor([0.0005, 0.0004]) \n",
      "Loss : -0.03216277062892914 , Params : tensor([0.2261, 0.1235]), Grad : tensor([0.0003, 0.0005]) \n",
      "Loss : -0.032164569944143295 , Params : tensor([0.2239, 0.1211]), Grad : tensor([0.0004, 0.0003]) \n",
      "Loss : -0.032165948301553726 , Params : tensor([0.2219, 0.1190]), Grad : tensor([0.0003, 0.0003]) \n",
      "Loss : -0.032166987657547 , Params : tensor([0.2202, 0.1171]), Grad : tensor([0.0002, 0.0003]) \n",
      "Loss : -0.03216775879263878 , Params : tensor([0.2186, 0.1155]), Grad : tensor([0.0002, 0.0002]) \n",
      "Loss : -0.032168325036764145 , Params : tensor([0.2173, 0.1140]), Grad : tensor([0.0002, 0.0002]) \n",
      "Loss : -0.032168734818696976 , Params : tensor([0.2162, 0.1128]), Grad : tensor([0.0002, 0.0001]) \n",
      "Loss : -0.032169025391340256 , Params : tensor([0.2152, 0.1117]), Grad : tensor([0.0001, 0.0001]) \n",
      "Loss : -0.03216923400759697 , Params : tensor([0.2143, 0.1108]), Grad : tensor([0.0001, 0.0001]) \n",
      "Loss : -0.03216937929391861 , Params : tensor([0.2136, 0.1101]), Grad : tensor([8.7368e-05, 8.8200e-05]) \n",
      "Loss : -0.03216947987675667 , Params : tensor([0.2130, 0.1094]), Grad : tensor([7.0838e-05, 7.2052e-05]) \n",
      "Loss : -0.03216954320669174 , Params : tensor([0.2125, 0.1089]), Grad : tensor([5.7016e-05, 5.8495e-05]) \n",
      "Loss : -0.032169587910175323 , Params : tensor([0.2121, 0.1084]), Grad : tensor([4.5736e-05, 4.6985e-05]) \n",
      "Loss : -0.03216962143778801 , Params : tensor([0.2117, 0.1081]), Grad : tensor([3.6817e-05, 3.7105e-05]) \n",
      "Loss : -0.0321696400642395 , Params : tensor([0.2115, 0.1078]), Grad : tensor([2.9075e-05, 2.9367e-05]) \n",
      "Loss : -0.0321696512401104 , Params : tensor([0.2112, 0.1075]), Grad : tensor([2.2087e-05, 2.3712e-05]) \n",
      "Loss : -0.032169658690690994 , Params : tensor([0.2110, 0.1073]), Grad : tensor([1.8022e-05, 1.7653e-05]) \n",
      "Loss : -0.03216966614127159 , Params : tensor([0.2109, 0.1071]), Grad : tensor([1.3290e-05, 1.4153e-05]) \n",
      "Loss : -0.03216966986656189 , Params : tensor([0.2108, 0.1070]), Grad : tensor([1.0965e-05, 1.0026e-05]) \n",
      "Loss : -0.03216967359185219 , Params : tensor([0.2107, 0.1069]), Grad : tensor([6.9383e-06, 8.7983e-06]) \n"
     ]
    }
   ],
   "source": [
    "func = ODE(init_params = tc.tensor([0.5, 0.3]), y0=y0)\n",
    "\n",
    "optimizer = tc.optim.Adam(func.parameters(), lr=0.1)\n",
    "\n",
    "losses = []\n",
    "tol = 1e-5\n",
    "g = tc.tensor(1.0)\n",
    "i = 0\n",
    "while g > tol:\n",
    "    optimizer.zero_grad()\n",
    "    l = likelihood(func, data, y0, tspan)\n",
    "    \n",
    "    l.backward()\n",
    "    if i%5==0 : print(\"Loss : {} , Params : {}, Grad : {} \".format(l, func.params.clone().detach(), func.params.grad))\n",
    "    optimizer.step()\n",
    "    g = func.params.grad.square().sum().sqrt()\n",
    "    losses.append(l.clone().detach())\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc099b2da60>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArHUlEQVR4nO3deXxV9Z3/8dfnLllIwhIIIYAIKi4sYiGC1Na6sLhVXIqtdSq209r+xv5+v+m0M9XxN0Or06mdrtPaDkWLw3Sdjq1CayuNqHWpVIJFZNOgghDDGiBhyf75/XHPDRe84SZcyE1y38/HI497z/d+77mfo4E33+/33HPM3REREemKUKYLEBGR3kfhISIiXabwEBGRLlN4iIhIlyk8RESkyyKZLqA7DBkyxEePHp3pMkREepVVq1btdveSZK9lRXiMHj2aysrKTJchItKrmNmWjl7TtJWIiHSZwkNERLpM4SEiIl2m8BARkS5TeIiISJcpPEREpMsUHiIi0mUKj0549vVdvL3nUKbLEBHpMRQenfB/f/EXFr3wVqbLEBHpMRQenXCwqZWG5tZMlyEi0mMoPFJwd5pa2mhu1R0XRUTi0goPMys2swozqwoeB3XQb17Qp8rM5iW0P2Fmr5jZOjNbYGbhoH1u0NZmZuUJ/Ueb2WEzWx38LEin/s5obGkDoKWt7VR/lIhIr5HuyOMuYLm7jwWWB9tHMbNiYD4wDZgKzE8ImZvdfRIwASgB5gbta4EbgWeTfOYb7n5B8POZNOtPqbE5FhrNrQoPEZG4dMNjDrA4eL4YuD5Jn9lAhbvXuvteoAK4EsDd64I+ESAH8KB9g7u/lmZtJ0VDS2ytQ9NWIiJHpBsepe5eEzzfDpQm6TMC2JqwvS1oA8DMlgE7gXrgkU585hgz+4uZ/dHM3n9iZXdefOTRopGHiEi7lPfzMLMngWFJXronccPd3cy6/M9zd59tZnnAT4HLiY1MOlIDjHL3PWY2BXjMzMYnjGAS674DuANg1KhRXS2rXXzk0dKmkYeISFzK8HD3GR29ZmY7zKzM3WvMrIzYCOJY1cClCdsjgWeO+YwGM1tCbBqsw/Bw90agMXi+yszeAM4G3nWnJ3dfCCwEKC8vP+G/+eMjj6YWjTxEROLSnbZaCsTPnpoHLEnSZxkwy8wGBQvls4BlZlYYBA5mFgGuATYe78PMrCThjKwzgLHAm2kew3Fp5CEi8m7phsf9wEwzqwJmBNuYWbmZPQTg7rXAfcDK4OfeoK0AWGpma4DVxEYtC4L332Bm24DpwOPBugjAJcAaM1tNbH3kM8G+ThmteYiIvFta9zB39z3AFUnaK4FPJmwvAhYd02cHcGEH+30UeDRJ+6+AX6VTc1fFv1nepLOtRETa6RvmKbR/SVAjDxGRdgqPFOIjD615iIgcofBIIT7y0NlWIiJHKDxSODLyUHiIiMQpPFI4suahaSsRkTiFRwpHzrbSyENEJE7hkYJGHiIi76bwSEFrHiIi76bwSCE+8mhuddw1+hARAYVHSo0J9y7Xdz1ERGIUHik0Jny/Q+seIiIxCo8UGhJGHs1a9xARARQeKSWOPCo313Lx/U+x71BTBisSEck8hUcKiSOPddV1VO87zKadBzJYkYhI5ik8UkgcedQ3tgBQs78hU+WIiPQICo8UEkce9Q2x8Niu8BCRLKfwSKGxpY1o2AA4qJGHiAig8EipobmVwtzYDRfrG5oB2F53GIDH/lLN7G8/qy8PikjWUXik0NjSRmFeLDwOHDPy2LC9jtd21OuiiSKSdRQeKcRGHlHg3Wsejc1tQR+Fh4hkl7TCw8yKzazCzKqCx0Ed9JsX9Kkys3kJ7U+Y2Stmts7MFphZOGifG7S1mVn5Mfs638xeDF5/1czy0jmG43F3GlvaKMo9euSxs76R1jbncFNsMT3xEiYiItkg3ZHHXcBydx8LLA+2j2JmxcB8YBowFZifEDI3u/skYAJQAswN2tcCNwLPHrOvCPAT4DPuPh64FGhO8xg6FD9NtyA3DBwJj9Y2Z/eBRhpaYqGhkYeIZJt0w2MOsDh4vhi4Pkmf2UCFu9e6+16gArgSwN3rgj4RIAfwoH2Du7+WZF+zgDXu/krQb4+7n7J/9sfDozAvNm11IJi2gti6R3zkEQ8REZFskW54lLp7TfB8O1CapM8IYGvC9ragDQAzWwbsBOqBR1J83tmAm9kyM3vZzP6ho45mdoeZVZpZ5a5duzpxKO8Wn46Kn23V0ubkRWP/ybbvP8zh5vjIQ+EhItklZXiY2ZNmtjbJz5zEfh47X7XL56y6+2ygDMgFLk/RPQK8D7g1eLzBzK7oYL8L3b3c3ctLSkq6WhZwZORRFJxtBXB6cQEQWzTXgrmIZKtIqg7uPqOj18xsh5mVuXuNmZURG0Ecq5rY2kTcSOCZYz6jwcyWEJsGqzhOOduAZ919d/D5vwMmE1tvOeniI4qCnCP/mYYU5fDaDjjY1KqRh4hkrXSnrZYC8bOn5gFLkvRZBswys0HBQvksYJmZFQaBE18IvwbYmOLzlgETzaxf8J4PAOvTPIYOHVnzOBIe/fOiRELGwcYWhYeIZK10w+N+YKaZVQEzgm3MrNzMHgJw91rgPmBl8HNv0FYALDWzNcBqYqOWBcH7bzCzbcB04PFgXYRgwf1bwX5WAy+7++NpHkOHCnMjfHDScEYP7tfelhcN0y8nzKGm1vbQaGjRtJWIZJeU01bH4+57gHetObh7JfDJhO1FwKJj+uwALuxgv48Cj3bw2k+Ina57yo0eUsD3bnkPu+ob29tyIyEKciMcbGw5Eh4aeYhIltE3zDshJ3zkP1NeNExBboRDTa36kqCIZC2FRydEgqvqAuRGQxTkhDnQ2NI+XaWzrUQk2yg8OiExPPIiYfrlRNh/uJnWttiZyZq2EpFso/DohGjo2GmrMLUHj9zHXN8wF5Fso/DohFDICIdio4+8aIh+ORH2HDiyiK5pKxHJNgqPTooE4ZEbiY08DjYdGW1o2kpEso3Co5PiZ1zFRx6JNPIQkWyj8Oik+KJ5XjRMQU74qNe05iEi2Ubh0UmRxJFH7tEjD33PQ0SyjcKjk9qnrSJHjzzCIWuftmpobmX/oVN2byoRkR5D4dFJ8Wmr3Gj4qDWPgflRGppbaWppY+KXlnH1d5/LVIkiIt1G4dFJR862CrXflhZgYL8oDS2tfPMPr9Hc6lTvO0xbm/Olpet4Y9eBTJUrInJKKTw6Kdq+5nHMyKNfDg3NbfzXi1uAWLhsr2vgP/+0mac2JLu9iYhI76fw6KRowoJ54shjUL8ou+obOdzcSv+8CI0tbe3fPq9r0PqHiPRNCo9OSjxV99iRx/7DsZA4c2ghAO/sOwxAfUNLN1cpItI9FB6dFL++Vex7HrHwCNnR9zc/s+To8Kg7rJGHiPRNCo9OikaOLJj3C6at8qJh8qJHprDaw2N/AwB1GnmISB+l8OikSChEOGREw6H2kUd+NExeJDE8CgCojo88tOYhIn2UwqOTomEjL3Jk0dwsPvKItfXLCTNsQB6gNQ8R6fsUHp0UDYfap6jMjIKcCHnRI21Di3IpyosCWvMQkb4vrfAws2IzqzCzquBxUAf95gV9qsxsXkL7E2b2ipmtM7MFZhYO2ucGbW1mVp7Q/1YzW53w02ZmF6RzDJ0VSQgPiI008nOOjDyG9s+jMLjm1c762L0+6jVtJSJ9VLojj7uA5e4+FlgebB/FzIqB+cA0YCowPyFkbnb3ScAEoASYG7SvBW4Enk3cl7v/1N0vcPcLgI8Bb7n76jSPoVNmjivlpskj2rcLcyOxNY+jRh6RoM5Yn/rGFtqCW9WKiPQlkdRdjmsOcGnwfDHwDPDFY/rMBircvRbAzCqAK4Gfu3tdQh05gAO4+4ag7/E++xbgF2nW32nXTRp+1HZhXoT8nAi5kXh45JEbCRENG82tscBwh4NNLe3TWSIifUW64VHq7jXB8+1AaZI+I4CtCdvbgjYAzGwZsRHJ74FHuvDZHyYWXkmZ2R3AHQCjRo3qwm47556rzyMvGmbvodi3yUv752JmFOVFj7q/eV2DwkNE+p6U01Zm9qSZrU3yc9Rf3O7uBCOHrnD32UAZkAtc3pn3mNk04JC7rz3Ofhe6e7m7l5eUlHS1rJSmnTGYSacNPDJt1T8XoH3dI37Pc617iEhflHLk4e4zOnrNzHaYWZm715hZGZDsSoDVHJnaAhhJbHor8TMazGwJsZFERSfq/gjw8070O+VGDy6gtH8uE0cMBI5847xsQB7b9h6m7nAL7s66d+oYP7x/qqk4EZFeId0F86VA/OypecCSJH2WAbPMbFCwUD4LWGZmhUHgYGYR4BpgY6oPNLMQcDPduN5xPMMG5PHnf5zBWcF1reIjjxED84HYyON7T23i2u89z5/fqs1YnSIiJ1O64XE/MNPMqoAZwTZmVm5mDwEEC+X3ASuDn3uDtgJgqZmtAVYTG7UsCN5/g5ltA6YDjwfrInGXAFvd/c00az8l4usbIwbFwuPxNTV8q+J1AKr3xu710aozsESkl0trwdzd9wBXJGmvBD6ZsL0IWHRMnx3AhR3s91Hg0Q5eewa46ISLPsXi01Yjg5HHY6urGdgvyr5Dzew+0Mh9j6/nte31/OxTPfYQRERS0jfMT7J4eMRHHm0Ol4wtITcSYs/BJtZV17Fxe30mSxQRSZvC4ySLr3kMLsglN7gW1rQzihlSmMvu+kZ21DdQe7CJ5ta2TJYpIpIWhcdJFl/z6J8fbX8+bcxghhTmsOtAI9uDy7UnfhdERKS3UXicZIXBtFVRXoT++RGGFOZyZkkBQwpzeXPXQRpbYiOOXcH1r0REeqN0v2Eux5g0cgDnlBYxclA+F55eTP/8CGbGkMJcqvcd+RrM7gMKDxHpvRQeJ9n5Iwey7HOXAPC1D53f3j64MOeofrsPaNpKRHovTVt1kyGFuUdta+QhIr2ZwqObDCk6Eh45kZDWPESkV9O0VTcZUhCbthrUL3YW1u4DjTS1tJETUX6LSO+jv7m6SXzkUdo/jyGFOby2vZ733PsHFj77RoYrExHpOoVHN4mveQwbkMeQwlw2bq/nYFMr//bEa6zeui+zxYmIdJHCo5sMzI8SDhnD+ue1j0LKBuRRXJDDd5dXZbg6EZGuUXh0k1DI+LuZZ/OhKSMpCUYhs8cP471nDmZjTR3V+w5z1b8/x9baQxmuVEQkNYVHN7rzsrMoH11MSTDymHFeKWcPK+Kd/Q08vuYdNtTUsWrL3gxXKSKSms62yoCrJgzjUFML088cTENzKwCPrNoGoJGHiPQKCo8MGFyYyx2XnAnAOcOKAHh9xwEAtu5VeIhIz6dpqwwbMTCf/Gi4fXtr7eEMViMi0jkKjwwLhYyxpbH7nw/Ij2rkISK9gsKjBzi7NDZ1NWtcKTX7G/jyb9bxycUrM1yViEjH0goPMys2swozqwoeB3XQb17Qp8rM5iW0P2Fmr5jZOjNbYGbhoH1u0NZmZuUJ/aNmttjMXjWzDWZ2dzr19xQ3l5/Gp94/himnD6K1zfnJii0s37iTPbp4ooj0UOmOPO4Clrv7WGB5sH0UMysG5gPTgKnA/ISQudndJwETgBJgbtC+FrgRePaY3c0Fct19IjAF+LSZjU7zGDJu6phi7rlmHKcV9wOgudVxh+c37c5wZSIiyaUbHnOAxcHzxcD1SfrMBircvdbd9wIVwJUA7l4X9IkAOYAH7Rvc/bUk+3KgwMwiQD7QBNQl6dcrjRyUD0BeNMTAflH++NquDFckIpJcuuFR6u41wfPtQGmSPiOArQnb24I2AMxsGbATqAceSfF5jwAHgRrgbeAb7l6brKOZ3WFmlWZWuWtX7/hLePjAfMIh45KxJbx/bAnPVu2mrc0zXZaIyLukDA8ze9LM1ib5mZPYz92dYOTQFe4+GygDcoHLU3SfCrQCw4ExwOfN7IwO9rvQ3cvdvbykpKSrZWVENBzi3246n7+ffQ4zzhvK7gON/GH9DmoPNlHX0Jzp8kRE2qX8kqC7z+joNTPbYWZl7l5jZmXERhDHqgYuTdgeCTxzzGc0mNkSYtNgFccp56PAE+7eDOw0sxeAcuDNVMfRW9w0ZSQAY4YU8MBTm/jyb9ZR39DCJWcP4Qe3TslwdSIiMelOWy0F4mdPzQOWJOmzDJhlZoOChfJZwDIzKwwCh2AN4xpgY4rPe5tgdGJmBcBFnXhPrxQJh7j76nOp2d/AgcYWVrxZS2xwJyKSeemGx/3ATDOrAmYE25hZuZk9BBCsSdwHrAx+7g3aCoClZrYGWE1s1LIgeP8NZrYNmA48HqyLAHwfKDSzdcG+Hnb3NWkeQ4912TlD+dmnpnHXVedSe7BJ3z4XkR7DsuFfs+Xl5V5ZWZnpMk7Y2ur9XPu95/nsZWfx+7U1/GjehYweUpDpskSkjzOzVe5enuw1fcO8Fzh3WBF50RDff2YTb+w6yPKNyZaWRES6j8KjF4iEQ0wcMYD4IPGlt/ZktiARyXoKj17iojMGU5gb4fJzh/LSW1o8F5HMUnj0Ep+9/Cye+sIHuHLCMPYeambl5r269pWIZIzCo5fIjYQZWpTHtDHFANz8wxe57oEXaGxpzXBlIpKNFB69zKjiflw9cRhXnDuU6n2H+dmf3850SSKShRQevYyZ8YNbp/DQvHKmnzGY7z+9iQONLZkuS0SyjMKjlzIzvnjVuew+0MQDT23KdDkikmVSXttKeq4LThvIh6aM5EfPv0lTSxvXTipj8qik9+MSETmpNPLo5b545bmcM6yIH6/YzD2Prs10OSKSJRQevVxJUS6//d/v5+6rzmNDTR2bdtZnuiQRyQIKjz7i2vPLCBn8+MUt/OaVd2jVTaRE5BTSmkcfMbR/HhedMZjFL25h8YtbaGlr44b3jMx0WSLSR2nk0Yd8ftbZfPzi0YwYmM+vX67OdDki0odp5NGHTDm9mCmnF1OUG+GBpzexfX8DwwbkZbosEemDNPLog26YPJI2hznff55vV7ye6XJEpA9SePRBY4YUcO+c8YweXMC/L6/i5bf3ZrokEeljFB591G3TR7Po9gspKcrlK49voE1nX4nISaTw6MMKciP8/exzWLVlL/+8dK3uASIiJ40WzPu4uVNG8uaugyz44xsAfPm6CYRDluGqRKS3S2vkYWbFZlZhZlXBY9ILK5nZvKBPlZnNS2h/wsxeMbN1ZrbAzMJB+9ygrc3MyhP655jZw2b2avC+S9OpPxuYGV+88hw+/YEz+MmKt7nvt+szXZKI9AHpTlvdBSx397HA8mD7KGZWDMwHpgFTgfkJIXOzu08CJgAlwNygfS1wI/DsMbv7FIC7TwRmAt80M029pWBm3H3Vedw2/XQWv7iZtdX7M12SiPRy6f7FOwdYHDxfDFyfpM9soMLda919L1ABXAng7nVBnwiQA3jQvsHdX0uyr3HAU0GfncA+oDxJP0niC7PPYXBBDl/4n1d0BpaIpCXd8Ch195rg+XagNEmfEcDWhO1tQRsAZrYM2AnUA4+k+LxXgOvMLGJmY4ApwGnJOprZHWZWaWaVu3bt6tTB9HX986J89cbz2V7XwI0/+BPPVem/i4icmJThYWZPmtnaJD9zEvt57FSeLp/O4+6zgTIgF7g8RfdFxMKnEvgO8Ccg6U283X2hu5e7e3lJSUlXy+qzZo4r5YUvXs5pxfl87YmNOoVXRE5IyvBw9xnuPiHJzxJgh5mVAQSPO5PsopqjRwcjg7bEz2gAlhCbBjteLS3u/jl3v8Dd5wADAX2FuosKciN8bsbZrK2u42tPbGTfoaZMlyQivUy601ZLgfjZU/OIBcCxlgGzzGxQsFA+C1hmZoUJwRMBrgE2Hu/DzKyfmRUEz2cCLe6u04dOwJwLRnD1xGH88Nk3mbvgRZpa2jJdkoj0IumGx/3ATDOrAmYE25hZuZk9BODutcB9wMrg596grQBYamZrgNXERi0LgvffYGbbgOnA48G6CMBQ4GUz2wB8EfhYmvVnrXDI+MGtU1jwV1Oo2nmA/3pxc6ZLEpFexLLhW8fl5eVeWVmZ6TJ6rI8//BIvvVXLN+ZO4qqJZZkuR0R6CDNb5e5Jz2jVdySEr9wwkTNKCvlfP31ZIxAR6RSFhzB8YD6//pv3ctk5JfzL4xtY/05d6jeJSFZTeAgA0XCIb8ydxMD8KJ9cvJKttYcyXZKI9GAKD2k3uDCXhz9+IQebWrnlwRW8s+9wpksSkR5K4SFHGT98AD/+66nsP9TMRx9cwa76xkyXJCI9kMJD3uX8kQNZ/NdT2V7XwGd+sorGlqRf4heRLKbwkKQmjxrEt26+gFVb9vLpH6/iYGNLpksSkR5E4SEdunpiGV+9cSLPVe1m3qKX9C10EWmn8JDjumXqKL794Quo3LKXf/3dBt3KVkQA3YZWOuG6ScP5y9t7efiFzRxuauXe68eTGwlnuiwRySCFh3TK/7tmHIW5Eb731CZCIfjqjednuiQRySCFh3RKOGR8ftY5tLnz/aff4JzSIm6/eEymyxKRDNGah3TJ3808hyvOHcqXfrOee3+zXmsgIllK4SFdEg4ZP/zYFG5/72gWvfAWC599M9MliUgGaNpKuiwSDjH/g+PYdaCR+5/YSH5OmNumj850WSLSjRQeckLMjG98aBKNza3885J11B5s4m9nnJ3pskSkm2jaSk5Yfk6YH36snJsmj+Q7T1bx4xVbMl2SiHQTjTwkLeGQcf9NE9l7qIl/emwtBxpa+F+XnpnpskTkFNPIQ9IWDYdY8FdT+OCk4XztiY18q+J1nYUl0selFR5mVmxmFWZWFTwO6qDfvKBPlZnNS2h/wsxeMbN1ZrbAzMJB+9fNbKOZrTGzR81sYMJ77jazTWb2mpnNTqd+OXlyIiG+8+ELmDtlJN9dXsU3/6AAEenL0h153AUsd/exwPJg+yhmVgzMB6YBU4H5CSFzs7tPAiYAJcDcoL0CmODu5wOvA3cH+xoHfAQYD1wJ/CAeOJJ54ZDxtZvO5yMXnsYDT2/iq7/fqAAR6aPSDY85wOLg+WLg+iR9ZgMV7l7r7nuJBcOVAO4ev1l2BMgBPGj/g7vHrwG+AhiZ8Hm/cPdGd38L2EQskKSHCIWMf71hIh+76HQWPvsmd//6VVrbFCAifU264VHq7jXB8+1AaZI+I4CtCdvbgjYAzGwZsBOoBx5J8v5PAL/vzL4SmdkdZlZpZpW7du3qxKHIyRIKGffOGc9nLzuLX6zcyv/5+V90QymRPibl2VZm9iQwLMlL9yRuuLubWZf/ienus80sD/gpcDmxkUn8s+8BWoLXurrfhcBCgPLycv3Tt5uZGV+YfQ4D8qN85XcbqGto5ocfm0K/HJ3gJ9IXpPyT7O4zOnrNzHaYWZm715hZGbERxLGqgUsTtkcCzxzzGQ1mtoTYtFRFsO/bgWuBK/zIxHk1cNox+6pOdQySOZ+65AwG5Ee569dr+KuH/szDt09lQL9opssSkTSlO221FIifPTUPWJKkzzJglpkNChbKZwHLzKwwCBzMLAJcA2wMtq8E/gG4zt0PHfN5HzGzXDMbA4wFXkrzGOQUu/nC0/j+RyfzavV+PrzwRXbWN2S6JBFJU7rhcT8w08yqgBnBNmZWbmYPAbh7LXAfsDL4uTdoKwCWmtkaYDWxUcuCYL8PAEVAhZmtNrMFwb7WAb8E1gNPAHe6uybTe4GrJpax6PYL2bLnEHMXvMjW2kOp3yQiPZZlw6mU5eXlXllZmekyBFi1ZS8ff/gl8nPCLP7EVM4d1j/TJYlIB8xslbuXJ3tN3zCXbjXl9EH88jPTAfjQf7zIM68lWyYTkZ5O4SHd7txh/Xnszos5rbgff724UhdUFOmFFB6SEWUD8vmfz0znA2eX8E+PreW+367XlwlFehGFh2RMYW6EB28r5/b3juZHz7/Fp39cycHGltRvFJGMU3hIRoVDxpeuG8+9c8bz9Gu7uOk//kT1vsOZLktEUlB4SI9w2/TRLLr9Qqr3HmbOAy+wakttpksSkeNQeEiP8YGzS3j0zvdSmBvmIwtX8N8r3850SSLSAYWH9ChnDS1iyZ3v46IzBvPFX73KPz22lqaWtkyXJSLHUHhIjzOgX5SHb7+QOy45gx+v2MItD65gR50uaSLSkyg8pEeKhEP849Xn8cBH38OGmjqu/d7zrNysdRCRnkLhIT3atecP59G/uZiCnDC3LFzBgj++QZu+DyKScQoP6fHOGVbEks++j5njSrn/9xuZ9/BLujKvSIYpPKRXGJAf5Qe3TuarN05k5eZarvrOczyt62KJZIzCQ3oNM+OWqaP4zWffR0lRLh9/eCX3/Xa9bnErkgEKD+l1xpYW8didFzNv+un86Pm3uPEHf+KNXQcyXZZIVlF4SK+UFw3z5TkTePC2ct7Zd5ir//05Hnz2TV1cUaSbKDykV5s5rpRlf3sJ7x9bwld+t4Gbf/iiRiEi3UDhIb3e0P55PHjbFL794Uls2nlAoxCRbqDwkD7BzLjhPSOp+NwlvH/sEL7yuw3c9B9/Ym31/kyXJtInKTykT4mNQsr5zocvYGvtIa574Hm+/Jt11Dc0Z7o0kT4lrfAws2IzqzCzquBxUAf95gV9qsxsXkL7E2b2ipmtM7MFZhYO2r9uZhvNbI2ZPWpmA4P2wWb2tJkdMLMH0qld+i4z4/r3jOCpz1/KR6eN4j//tJkrvvlHfvPKO7hrKkvkZEh35HEXsNzdxwLLg+2jmFkxMB+YBkwF5ieEzM3uPgmYAJQAc4P2CmCCu58PvA7cHbQ3AP8EfCHNuiULDOgX5V+un8hjf3Mxpf3z+N8//wu3PvRn1r9Tl+nSRHq9dMNjDrA4eL4YuD5Jn9lAhbvXuvteYsFwJYC7x/8UR4AcwIP2P7h7/H6kK4CRQftBd3+eWIiIdMqk0wby2J0Xc++c8ayvqeOa7z3H3//PK2zfr18jkROVbniUuntN8Hw7UJqkzwhga8L2tqANADNbBuwE6oFHkrz/E8Dvu1qYmd1hZpVmVrlr166uvl36mHDIuG36aP74hcv45PvGsGT1O1z2jWf4VsXrum+6yAlIGR5m9qSZrU3yMyexn8cmk7s8oezus4EyIBe4/JjPvgdoAX56Avtd6O7l7l5eUlLS1bdLHzWgX5R7rhnHk3/3Aa44byjfXV7FB77+DA899yYNzbrMiUhnpQwPd5/h7hOS/CwBdphZGUDwmOxKddXAaQnbI4O2xM9oAJYQmwYj2N/twLXAra5VTjnJRg3uxwMfncyv/+a9nF1ayL88voH3/9vT/Oj5txQiIp2Q7rTVUiB+9tQ8YgFwrGXALDMbFCyUzwKWmVlhQvBEgGuAjcH2lcA/ANe5+6E0axTp0ORRg/jZpy7iv++4iLNKCrnvt+t5/789zaLn3+JQk6azRDpi6fyj3swGA78ERgFbiJ09VWtm5cBn3P2TQb9PAP8YvO0r7v6wmZUCvyU2XRUCngY+5+4tZrYpaN8TvGeFu38m2NdmoD+xBfZ9wCx3X3+8OsvLy72ysvKEj1Oyx4o39/CdJ19nxZu1DMiPcuu0Udz+3tEM7Z+X6dJEup2ZrXL38qSvZcOMkMJDumrVlloefPYtlq3fTiRkXDdpBJ+6ZAznDuuf6dJEus3xwiPS3cWI9AZTTi9myseK2bLnIIuef4tfVm7jVy9vY+roYm6ZdhpXTSgjLxrOdJkiGaORh0gn7DvUxH+v3MrPX3qbzXsOMSA/yo2TR/DRqaMYW1qU6fJETglNWyk85CRpa3NWvLmHn730NsvWbae51Tl/5ACumzScD04aTqnWRqQPUXgoPOQU2HOgkV+/XM2SV6pZW12HGUw/YzBzLhjOrHHDGFSQk+kSRdKi8FB4yCn2xq4DLF39DktWV7N5zyHCIaP89EHMHFfKzHGlnD64INMlinSZwkPhId3E3Xm1ej8V63dQsX4HG7fXAzB2aCGXnzuUi88awoWji8nP0WK79HwKD4WHZMjW2kNUrN/Bkxt2sHJzLc2tTk44xOTTB/K+s4Yw/czBjB8+QGduSY+k8FB4SA9wqKmFlZv38sKm3bywaTfrgkvDR8PGuOEDmDxqIJNHDWLy6YMYPiAPM8twxZLtFB4KD+mB9hxopHLLXl5+ey9/2bKPNdX7aGhuA2BoUS7jh/fn3LL+nDusiHFl/RkzpIBIWDf/lO6jLwmK9ECDC3OZPX4Ys8cPA6C5tY2NNfW8/PZeVm/dx4aaOp6r2k1LW+wfeDmREGOHFjJ2aCGjhxQwZkgBowcXMHpIAQPyo5k8FMlCCg+RHiIaDjFx5AAmjhzQfrXRppY23th1gI3b69hYU8/6mjpWbt7LklfeIXHSoLggh9MH92P4gHzKBuQxbEAewwfmxx4H5FNSlEs4pGkwOXkUHiI9WE4kxHll/TmvrD+850h7Q3Mrb9ce4q3dB9m8+yCb9xxky55DbKipY/nGHe3TX3HhkDGkMIfiglwGF+RQHPwMLsihuDD22D8/Sv+8KEV5EYqCx6imyaQDCg+RXigvGubs0iLOTnJpFHdn/+Fm3tnXwPa6w7yzr4Ga/YfZVd9I7cEm9hxsYuveQ9QeaKI+xV0UcyMhivKi9M+LUJQXIT8nTF40TH409hj7CbVv5wfbucF2NGREwiEiYSMaCh7DRqT9eYhIKHgM2qNhIxQyQmaEDEJmWPAYb9PJBJmn8BDpY8yMgf1yGNgvh3HDj38V4MaW1ligHGiirqGZ+oYW6htaOBB/3thCfUMzdUF7Q1Mrew400dDcSkNLK4eb2mhsbuVwc2v72kx3CYeOBEkoIVyOBE08eGLP4xJzx7B3tSfGUkch1d43xb4S95e4r6P2mqTvyXTZOUP5f9eOO+n7VXiIZLHcSJiyAfmUDchPe1/NrW2xUGmOPTa2tNLc6rS0Os1tbbS0Oi2tbTS3BY+tTkvQ3tzaRktCe5s77tDmTlvw6AnP29oSnjvBa4l94685rW0Qv0N24jrRUc+Tvc7x+9JhX39Xe8f7fXffk61sYPr/b5NReIjISRENh4iGQxTp2pBZQathIiLSZQoPERHpMoWHiIh0WVrhYWbFZlZhZlXB46AO+s0L+lSZ2byE9ifM7BUzW2dmC8wsHLR/3cw2mtkaM3vUzAYG7TPNbJWZvRo8Xp5O/SIicmLSHXncBSx397HA8mD7KGZWDMwHpgFTgfkJIXOzu08CJgAlwNygvQKY4O7nA68Ddwftu4EPuvtEYB7w4zTrFxGRE5BueMwBFgfPFwPXJ+kzG6hw91p330ssGK4EcPe6oE8EyCE4W83d/+Du8W8vrQBGBu1/cfd3gvZ1QL6Z5aZ5DCIi0kXphkepu9cEz7cDpUn6jAC2JmxvC9oAMLNlwE6gHngkyfs/Afw+SftNwMvu3ngCdYuISBpSfs/DzJ4EhiV56Z7EDXd3M+vy91zcfbaZ5QE/BS4nNjKJf/Y9QEvwWmJN44GvAbOOU/cdwB0Ao0aN6mpZIiJyHCnDw91ndPSame0wszJ3rzGzMmIjiGNVA5cmbI8EnjnmMxrMbAmxabCKYN+3A9cCV3jC1zDNbCTwKHCbu79xnLoXAguD9+wysy3HOcxUhhBbb8lWOn4dv44/O53e0Qtp3QzKzL4O7HH3+83sLqDY3f/hmD7FwCpgctD0MjAFaAKKguCJEBtdPOfuD5jZlcC3gA+4+66EfQ0E/gh82d1/fcKFd5GZVXZ0Q5RsoOPX8ev4s/f4O5Lumsf9wEwzqwJmBNuYWbmZPQTg7rXAfcDK4OfeoK0AWGpma4DVxEYtC4L9PgAUARVmttrM4u2fBc4C/jloX21mQ9M8BhER6aKsuA1turL9Xx46fh2/jj97j78j+oZ55yzMdAEZpuPPbjp+eReNPEREpMs08hARkS5TeIiISJcpPI7DzK40s9fMbFNwKnKfZ2abgwtPrjazyqCtUxfA7K3MbJGZ7TSztQltSY/ZYr4b/E6sMbPJHe+5d+jg+L9kZtUJZzVenfDa3cHxv2ZmszNT9cljZqeZ2dNmtj64SOv/Ddqz5nfgRCg8OhBc4ff7wFXAOOAWMzv5NwLumS5z9wsSzjBJeQHMXu4/Ca63lqCjY74KGBv83AH8RzfVeCr9J+8+foBvB78HF7j77wCCPwMfAcYH7/lB/GrYvVgL8Hl3HwdcBNwZHGc2/Q50mcKjY1OBTe7+prs3Ab8g9g34bNSZC2D2Wu7+LFB7THNHxzwH+C+PWQEMDK6u0Gt1cPwdmQP8wt0b3f0tYBOxPyu9lrvXuPvLwfN6YAOx6+9lze/AiVB4dOy4F3Tswxz4Q3C/lDuCts5cALOv6eiYs+n34rPBtMyihKnKPn38ZjYaeA/wZ/Q7cFwKDznW+9x9MrGh+Z1mdknii8F1xrLq/O5sPGZiUzFnAhcANcA3M1pNNzCzQuBXwN8m3C4CyNrfgeNSeHSsGjgtYXtk0NanuXt18LiT2AUopwI74sPy41wAs6/p6Jiz4vfC3Xe4e6u7twEPcmRqqk8ev5lFiQXHTxOum5fVvwOpKDw6thIYa2ZjzCyH2CLh0gzXdEqZWYGZFcWfE7vk/Vpixx2/ffA8YElmKuxWHR3zUuC24Iybi4D9CVMbfcYxc/g3EPs9gNjxf8TMcs1sDLFF45e6u76TycwM+BGwwd2/lfBSVv8OpJLykuzZyt1bzOyzwDIgDCxy93UZLutUKwUejf1ZIgL8zN2fMLOVwC/N7K+BLcDNGazxpDOznxO7bcAQM9tG7LbJ95P8mH8HXE1sofgQ8PFuL/gk6+D4LzWzC4hN1WwGPg3g7uvM7JfAemJnKd3p7q0ZKPtkuhj4GPCqma0O2v6RLPodOBG6PImIiHSZpq1ERKTLFB4iItJlCg8REekyhYeIiHSZwkNERLpM4SEiIl2m8BARkS77/3uForfZSOAZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(losses[50:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc099c76520>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1/klEQVR4nO3deXhb93ng++8LgCDBHeC+SaQkardWWnK8ynZiy05sOZtjJ9M4jSdOWvu26cy0dW7vNJnO7TxpOpneyTOZ5DqNGyfXseM43pJ4UxyvsSWL2iiREk2KIiUuokgCFBeAG/C7f+DQgWlSCwkSAPF+ngcPD344B+flIXhe/JbzO2KMQSmlVHKzxToApZRSsafJQCmllCYDpZRSmgyUUkqhyUAppRTgiHUAs5Wfn28qKytjHYZSSiWU/fv39xpjCqaWJ2wyqKyspLa2NtZhKKVUQhGRtunKtZlIKaWUJgOllFKaDJRSSqHJQCmlFJoMlFJKoclAKaUUmgyUUkqRwNcZKKWio/HMIL871k2G0841KwtYXpAZ65BUDGgyUCpJTQRD/Odn63l83ykmb2vitNv4+seq+bPrliMisQ1QLShNBkolIWMM//nZeh579xRfvqqK+69fTmA8yH97/hjfebERY+D+61fEOky1gLTPQKkk9LM9bTz27in+fMdy/v62teRlplLuTuf7n9/CHZtK+e8vN/LKse5Yh6kWkCYDpZJM79Ao//xSI9dU5/PXN6/6wGsiwrc/vYHVxdn8n08fwT82EaMo1ULTZKBUkvnOi8cJjAX55m3rkJF+6DsBAd/7r6el2PmHXevoHhjloTdaYheoWlAXTAYi8rCInBWRoxFlvxCRQ9ajVUQOWeWVIhKIeO2HEdtsFZEjItIsIt8Tq3dKRDwisltEmqyf7nn4PZVSQGvvML/c386XrqxkRdY41D8Djc/D3h+Br/X99S6v9PDxy0r4f19v4ezASMziVQvnYmoGPwF2RhYYYz5njNlkjNkE/Ap4KuLlE5OvGWO+FlH+A+ArQLX1mHzPB4FXjDHVwCvWc6XUPPi3P5wkxWbjvuuWgd8LY0PhJNBzHGofgY4D79cS/vrmVYxMBPnpO9POeKwWmQsmA2PMG4B3utesb/d3Ao+d7z1EpATINsbsMcYY4KfAHdbLu4BHrOVHIsqVUlF0zj/OE7Xt3L6plMKsNEj3wFgA/D5wpoO3Cd57OVxbCPiozM/go2uKeHRvGyPjwViHr+bZXPsMrgG6jTFNEWVVInJQRF4XkWussjKgPWKddqsMoMgY02UtnwGKZtqZiNwnIrUiUtvT0zPH0JVKLr+oPUVgPMi9V1eFv/37vbBqJxSshvQ8sKdC/nIwofBrwL1XV+Hzj/P0wY4YR6/m21yvM7ibD9YKuoAlxpg+EdkKPCMi6y72zYwxRkTMeV5/CHgIoKamZsb1lFIfZIzhyf3tbF3qZk2aF/Y+Ac40cGbBxjvDtYO2t2HwTLi2YLMDsL3Kw9qSbH72Tht3b1sS499CzadZ1wxExAF8CvjFZJkxZtQY02ct7wdOACuBDqA8YvNyqwyg22pGmmxOOjvbmJRS06vvHOBM9xnureoP9w30HAdfG4wNQigIZZthzcfDicCZBidehYAPEeFzl1fQ0DXA8TMDsf411DyaSzPRR4Hjxpj3m39EpEBE7NbyMsIdxS1WM9CAiFxh9TN8EXjW2uw54B5r+Z6IcqVUlLyw7xi3O95lh+0AeJvDfQR+H4yNhPsOIJwUMvIgrzrcsdx5CAI+bttYisMmPH1Am4oWs4sZWvoY8A6wSkTaReRe66W7+HDH8bVAnTXU9Enga8aYyc7nPwf+FWgmXGN4wSr/NvAxEWkinGC+PftfRyk1VTBkOFx3iG35Y6TnlYf7Blx54b6CjXeCyxrNne4BsUFfM3Qdhu4jUP8MHtswO1YV8vTBDoIhbZ1drC7YZ2CMuXuG8i9NU/YrwkNNp1u/Flg/TXkfcOOF4lBKzc7BxpNUjxxhc6kPTr4JBWtgxQ3gqfpjIoDw8ro7wjUCTLiGMNAJfi+f3lLG7451886JPq6uzo/Rb6Lmk16BrNQit7e+iZDNSf6WXZBfDStuhLItH0wEk1xuKN0U7lge6AzXFNI9XL+6EFeKnRfruz68jVoUNBkotYgZY3i+ZYzKwixcwSHILAJP5fk3mqwhlG6Bko1AeIqKHasKeLm+m5A2FS1KmgyUWsSazg5R77WRueWzsGxH+CQ/XY1gOl2HoPPg+xeh7VxfzNnBUQ6e7p+3eFXsaDJQahHb3RCehvraDdWQt/ziE4HfG774LLv0/ZFF11c6SbELL9WfmceIVaxoMlBqEas91sJNxUMUpQQubcNpRhZlN/+GG6tSeVmTwaKkyUCpRWqwv4fSzhf5jLvp/aaeizbZb1C0Hko2hEcWmRAfW2Kntc9Pa+/wfIWtYkSTgVKL1OHjJyAUYknlyg/MN3TRpo4sGh/hI2UpZDPEa406UcBio8lAqUXq9fYQDoeD5ann3h8ieskiRxYBpf7j3JN9iL3H9KY3i81cJ6pTSsWp3SdHuWzpraRULw0ngovtPJ7K5YY0L6SkQXYpq4o7eamljcBYEJfTHt2gVcxozUCpRei0109rn58tq6subRTRTCY7lAc6qS7K4cxEBnta+qITrIoLmgyUWoQmT9RXrYjS1BGTzUXLdrD0mi8wYs/m7RO90XlvFRe0mUipRWhPixdPhpPqwszovalVu0jze7m63M7bJ7RmsJhozUCpRWhPSx/XL01BvC2XNqT0fAK+8BDVltf4fOYB2rs66fePRee9VcxpMlBqkTnt9TPY38OnU/dBy2uXfo3BTCKuSl6Rn46bQfa0XOJwVRW3NBkotcjsaenDI4Msy0sPTycxm2sMphPRiVyel0HAkcs72m+waGifgVKLzOGmNkpTxylM5wPTUM/ZZCey34sj3cOqo03ab7CIaM1AqcUk4MN98rfcmt8V/ucu3XxpM5VeiMv9/lDVK5fn0XR2iJ7B0ei8t4opTQZKLSK+3jN4hwLklVSGLxJLy4leIogU8HFd3gDZDPGOXm+wKFzMPZAfFpGzInI0ouxbItIhIoesx60Rr31DRJpFpFFEbo4o32mVNYvIgxHlVSKy1yr/hYg4o/kLKpVMjnYN45ZB1jrPRq95aCprVNHK4Vo+lbaPg40no78PteAupmbwE2DnNOX/YozZZD2eBxCRtcBdwDprm/8tInYRsQPfB24B1gJ3W+sC/JP1XisAH3DvXH4hpZJWwMdQw27GxElJBrD8+vmpFVijiuy55SzPS6fxZGv096EW3AWTgTHmDeBihyLsAh43xowaY04CzcA269FsjGkxxowBjwO7RESAG4Anre0fAe64tF9BKQWA38vpvkGMpxpnVgGEgvOzn4hRRVUFWRz1Oejov8T7Jai4M5c+gwdEpM5qRpr8+lEGnI5Yp90qm6k8D+g3xkxMKZ+WiNwnIrUiUtvT0zOH0JVafMacbk73j7Ipd3j+mojgAzOZlq69AoB3dFRRwpttMvgBsBzYBHQB341WQOdjjHnIGFNjjKkpKChYiF0qlTAa+m08O76NjNUfje4Iopl0HaIy0Mhn0mo52tw2v/tS825WycAY022MCRpjQsCPCDcDAXQAFRGrlltlM5X3Abki4phSrpS6RPvbfAyQyZr1m+c/EVj9BrbcMqo8aZw4pckg0c0qGYhIScTTTwKTI42eA+4SkVQRqQKqgXeBfUC1NXLISbiT+TljjAFeBT5jbX8P8OxsYlIqqQV8nG46wuqcIEXZafO/v4h+g4r8LA732ekb0usNEtkFr0AWkceAHUC+iLQD3wR2iMgmwACtwFcBjDH1IvIE0ABMAPcbY4LW+zwAvATYgYeNMfXWLv4WeFxE/m/gIPDjaP1ySiWFgA9z9Glc7Q3ck58FgWvmv2YQcTVybo6NgboG9rX62Lm+eH73q+bNBZOBMebuaYpnPGEbY/4R+Mdpyp8Hnp+mvIU/NjMppS6V30u/f5T3Atlc604NN+HMdzKA8D5cbtbkBHE6jlPb6tVkkMD0CmSlEl26h1O+EYrxsiQvc/5GEc0gdXyAnSXDNLRov0Ei02SgVKJzuXkhdAX77Rsp+sjdC1MrmGRdjfyJzONUnd2N/5zOYpqoNBkotQi80xkkp3w1jsyFrRVMjioqrlgOoRDHT+jUFIlKk4FSCW5sIsSxrkE2VeQu/M6tUUXVaQOExMbes7LwMaio0PsZKJXgjp8ZYCwYYkN57sLv3BpV5PJ7ad67hNPtQf5s4aNQUaA1A6US3OH2cwBsKM+JTQDWPQ7WLKvgwCkfE8FQbOJQc6LJQKkE19R6mo3pfZSnxfair8srPfjHgjR0DcQ0DjU7mgyUSmQBH57W3/KpnCak4Zno3Ph+lnFckdNPNkPsa41RDGpONBkolcAC587SOxggp3hp9G58f8lBhIeXFvS8wxcyD+j1BglKk4FSCayh30nQ2FjlGpjfaavPxxpeSnYplXlptJ46RXjaMZVINBkolcAO9hh+G9pO0aZbFmba6ulETFq3JC+TpqFU2n16s5tEo8lAqQR2uP0cmTn5eJasiU0igD9OWrdsB3kbduKRQeqaWmMTi5o1TQZKJaqAj75TDVxRYo91JOGEkO5h+eA+bnAeY/zIr2LXma1mRZOBUoko4GP44JNUndvHbc698XHi9XuxY8gsqOB033BsOrPVrGkyUCoR+b10eoc5g4eKnLT4OPFafQcbc/x0DowxYM+OdUTqEmgyUCoRpXto7w9PW13qzojNKKKprL6DnLUf47fB7Rw8G+uA1KXQZKBUInK5eWZ8O6fc20jf/OnYdR5P5XKzev1mhiST/a1xUFtRF02TgVIJak9XkPwla+MnEVgyUx2sKclm/6k46MdQF+2CyUBEHhaRsyJyNKLsn0XkuIjUicjTIpJrlVeKSEBEDlmPH0Zss1VEjohIs4h8T0TEKveIyG4RabJ+xtcnW6k41D0wQvfAaOwmp7uAmqVuDp7q10nrEsjF1Ax+AuycUrYbWG+M2QC8B3wj4rUTxphN1uNrEeU/AL4CVFuPyfd8EHjFGFMNvGI9V0rNJOCjqeEQ2QzFZtrqCwn4uCbvHI6xcxw/MxjraNRFumAyMMa8AXinlL1sjJmwnu4Bys/3HiJSAmQbY/aY8HXqPwXusF7eBTxiLT8SUa6UmsqaByjQ+Aq3Ofayzh1n37yt+LYG6/i4ba9efJZAotFn8GXghYjnVSJyUEReF5FrrLIyoD1inXarDKDIGNNlLZ8BiqIQk1KLkzUP0JHBDEqyU0kb7491RB9kxecursTtstHUqpPWJYo5JQMR+TtgAnjUKuoClhhjNgP/Afi5iFz0YGOr1jDjDFcicp+I1IpIbU9PzxwiVypBpXswCIHe05S60+NjSGmkiHmKyvOyeLtTJ6xLFLNOBiLyJeATwBeskzjGmFFjTJ+1vB84AawEOvhgU1K5VQbQbTUjTTYnzTg62RjzkDGmxhhTU1BQMNvQlUpcLjftpTeze3QNobV3xN1Iosh5ilh7B40DDjr7ddK6RDCrZCAiO4G/AW43xvgjygtExG4tLyPcUdxiNQMNiMgV1iiiLwLPWps9B9xjLd8TUa6UmsbBXqHVlLCmakmsQ5medRvMDdWVANS26RDTRHAxQ0sfA94BVolIu4jcC/wvIAvYPWUI6bVAnYgcAp4EvmaMmex8/nPgX4FmwjWGyX6GbwMfE5Em4KPWc6XUDOpO95PqsLGyKCvWoZzX6pwga5xnqW9ujXUo6iI4LrSCMebuaYp/PMO6vwJ+NcNrtcD6acr7gBsvFIdSKqyu/RzrSrNJscfxNaMBH47jz/IZdzPDJ5ogEMMpttVFieNPk1JqquCwj8HO42wrjvN/XWtUUV5JFWfP+Rnu14mK4l2cf6KUUu8L+OjZ+zhbJg5zE+/Ex7TVM7FGFa3NHCJobBzxXrARQsWYJgOlEoXfS4d3KDxtdW5qfExbPRNrVFHplp08b7aztyvOLo5TH6LpWqlEke7hdP8IS1POkZdVHn/XGEzlcpPpclNSdJbatjhOXArQmoFSicPl5hf+GvqLr8S2/pMJ0yG71Zq0LhjSC9DimSYDpRJBwMfY2fd4r3uQoqp1CZMICPi4Ln8A22g/jTppXVzTZKBUvLMmf+s++CI3mT1sLYx1QBfJirvGmrTuyAmdpyieaTJQKt5ZwzTfC2RjI8QGTzDWEV2c9yetW0pOmo2mk62xjkidhyYDpeKdNUzT29VKRpqTgsLSWEd0cay4ZbCLirxM/qCT1sU1HU2kVLyzhmk+9+oLuMtLkPQE6S+YnLTO7yU0NMSxl9rpHhihKDst1pGpaWjNQKkEMCCZvNmXw8rK895HKv5Yk9ZdtmIpALWtcXyhXJLTZKBUAjjSfg6AjRW5sQ1klta5Q2xNOUlnw9vxfeV0EtNmIqUSwKHT/QBsKMuNaRyzEvCRUvcYX09/meAJ4OAZ2Pz5xBkemyS0ZqBUAqhr76cqP4Oc9JRYh3Lp/F4YGyQzx0N7wMGYvz++p9JIUpoMlEoAh0+fY2N5TqzDmJ10DzizKEkbJ8P4afenxP9UGklIk4FSce7MuRHODIywoTw31qHMjssNmz9P5vV/yb+ZW3kh5QZtIopDmgyUinP1LW1UShdbCiXWocyey01m1TZMyRbeODUe62jUNDQZKBXPAj4m6p7iGttR1vpeSfiRONurPBw83c/IeIJcRZ1ELioZiMjDInJWRI5GlHlEZLeINFk/3Va5iMj3RKRZROpEZEvENvdY6zeJyD0R5VtF5Ii1zfdEJIG/AikVRX4vnd5hHO5ynDYSvuN1W5WHsYkQddZQWRU/LrZm8BNg55SyB4FXjDHVwCvWc4BbgGrrcR/wAwgnD+CbwHZgG/DNyQRirfOViO2m7kuppBRK89B2boSNOX4QW8J3vG6rCse/t6UvxpGoqS4qGRhj3gCmfiXZBTxiLT8C3BFR/lMTtgfIFZES4GZgtzHGa4zxAbuBndZr2caYPcYYA/w04r2USmqt/hSeGrmclOrrw1M7JHjHa266k9XFWbzbmtg1nMVoLn0GRcaYLmv5DFBkLZcBpyPWa7fKzlfePk25UknvcHs/A2SyfPXGhE8Ek7ZXedjf5mM8qLfCjCdR6UC2vtHP+5SEInKfiNSKSG1PT898706p2Ar4ON10lGJngOrCrFhHEzVXljkoHG+noeVUrENREeaSDLqtJh6sn2et8g6gImK9cqvsfOXl05R/iDHmIWNMjTGmpqCgYA6hKxXnrBvDpLS9zp/mHMI+2h/riKIj4OPKkTe5So7Sv++JhB8dtZjMJRk8B0yOCLoHeDai/IvWqKIrgHNWc9JLwE0i4rY6jm8CXrJeGxCRK6xRRF+MeC+lkpPfy/jEBIcGMljiTk34UUTv83vJSrVBTimtPYOL5/daBC5qojoReQzYAeSLSDvhUUHfBp4QkXuBNuBOa/XngVuBZsAP/CmAMcYrIv8V2Get9w/GmMlPwp8THrHkAl6wHkolr3QP3YNj5Ad7Kc1bmvCjiN5n3fCmxj3C4Y4AwTQP9ljHpICLTAbGmLtneOnGadY1wP0zvM/DwMPTlNcC6y8mFqWSxeFgJYdNkD+7/DOLpvN48oY3WSMNPNXSzqd8woaMWAelQK9AVir+WP0F/pP7+IirnbJcV6wjii6Xm40rl+ORQWqPn4x1NMqiyUCpeGPdSP5Av4ulnjQksMja1QM+Ck69wB3ZTZj6p7UTOU5oMlAq3qR7GBwLYRvsosKTuXj6CyZZya6wbBmnewcZG9Rh4vFAk4FS8cblZl/6tfzBrCd762cXT3/BJKsTeVPuMKNB4VCv3nAxHmgyUCoO7T0TosNWypplS2IdSvRZncgVW2/lLbOeI00ntKkoDmgyUCoOHWjzsa40h7SURTrw0uUmy1PEXXktTDT9Huqf0YQQY5oMlIozk1M8b126yJqHpvJ7WZGfRq0vncDYmF6AFmOaDJSKMw1dA4xOhNiyZJEng3QPywtzKDR9tPaNLL6O8gSjyUCpOHOgLdxcsmVpbmwDmW8uN+VXf559tg38emL74usoTzCaDJSKJwEfp5qPsCp7gpKcRXax2TTSsvMoqFzL79vGYh1K0tNkoFS8sK48zjj9BvfkHkqaDtUrl+dz/MwgPYOjsQ4lqWkyUCpe+L30D49wPJBNpWcRzVR6AddU5wPwVrNefBZLmgyUihc2O2fOdFBJFxV5WUnTobrebdiU3sfehpZYh5LUNBkoFQ8CPjjxKq0Dhiz7OEWbbk2ODtWAD9uxZ/hCQQsZTb8mOJwcTWPxSJOBUvHAmq9nT7+bbE8RTluS3B/Y+r0rKlcwMjbO8RadxTRWNBkoFQ/SPfjHDRP97SwtSJ4mosl5ii7LGiJVxjnYfCppOs7jjSYDpeKBy83e9Gt5K7Se7K13JkcTEbw/T1FG1XbKcl0MNO/VqSliRJOBUnHi7Y4JOu1lXLZiaaxDWVguN6RlU1Wcx4F+F0OB0aQZSRVPNBkoFSf2nvSyqSJ38U5Odz7pHqqLsykyXt7r8SdPM1kcmXUyEJFVInIo4jEgIl8XkW+JSEdE+a0R23xDRJpFpFFEbo4o32mVNYvIg3P9pZRKNEOjExztOMf2qiQ9CbrclF35eRpcW/iFvyZ5msniyKyTgTGm0RizyRizCdgK+IGnrZf/ZfI1Y8zzACKyFrgLWAfsBP63iNhFxA58H7gFWAvcba2rVHII+Kg/coBMM8T2qrxYRxMztgw3VSs38FLLCMGQiXU4SSdazUQ3AieMMW3nWWcX8LgxZtQYcxJoBrZZj2ZjTIsxZgx43FpXqcXPmoLCd/RlbnPsZUthrAOKretWFdDvH+fQ6f5Yh5J0opUM7gIei3j+gIjUicjDIjJZ3ysDTkes026VzVT+ISJyn4jUikhtT49euq4WAWuc/cH+dCpyU0mfOBfriGJqx8pCPLZhag/s0xFFC2zOyUBEnMDtwC+toh8Ay4FNQBfw3bnuY5Ix5iFjTI0xpqagoCBab6tU7KR7GAtCoPcUS/Ozk77jNEeGuL/wKMGG5zF7HwJfa6xDShrRqBncAhwwxnQDGGO6jTFBY0wI+BHhZiCADqAiYrtyq2ymcqWSQgNVHAguI3PrZ7Tj1O9lQ6GNjMBpBk8fhcNPaA1hgUQjGdxNRBORiJREvPZJ4Ki1/Bxwl4ikikgVUA28C+wDqkWkyqpl3GWtq9TiZvUX9Da+wwZbKxsrcmMdUeyle1idl0KuGaLVnwrONL3mYIE45rKxiGQAHwO+GlH8HRHZBBigdfI1Y0y9iDwBNAATwP3GmKD1Pg8ALwF24GFjTP1c4lIqIVj9BXu9LipyU8kODgBJ3vzpcpO17QsMHu5k3zknG5xJNDVHjM0pGRhjhoG8KWV/cp71/xH4x2nKnween0ssSiUcaz6iQO8pVqwq0pPeJHclo5vv5We/P8DHCz5C8WTNINmb0OaZXoGsVKy43LyddjVvhdbj2XaXnuwi7Ni0Eq/JomXPr6HlNZ2vaAFoMlAqhl5tG6cnpZwN1Uk2H9EFLC/IZHNeiGOd/ZBdCiakfQfzTJOBUjH0ZlMvH1meR4pd/xWn2rp2OS19AYZ72kBs2ow2z/QTqFSMtPUNc8rr55rqJO80nsGOjSv59cR2Xp9YD+vu0Ga0eabJQKkYebOpF4CrrRvCqw+6rCwHd14hP29O0USwADQZKBULAR+N9QdZnRNkWX5GrKOJSyLCro2lvH2il7MDI7EOZ9HTZKDUQgv4CB55ipTTb/KVvMPISH+sI4pbt28qJdMM8caePTqaaJ5pMlBqofm9tHuHaRvLobogXUfJnMeKrAn+vecwZw+/pMNL55kmA6UWWrqHph4/JeJleZFOTndefi+by7PZ73PROxjQxDmPNBkoFQPP9xbiL9hIxmadnO680j1cVuGmRLwc7hjQxDmPNBkotZACPs7V/hJXTx23FZyNdTTxz+Umt+ZOBkuv5l+7qjDDfdpUNE80GSi1kPxejnf1cwYPa4sztdnjYrjcXL1xFZW+vXQefEH7DuaJJgOlFlK6h4YzQ6zJGKIw26XNHhfp5koHTpvh9S6nTk0xTzQZKLWAhm1ZfL/7MlyrbkTWf1L7Cy5Sdl4xq8tyOdHyHmMhNInOA00GSi2gN5t66A2ms2VzjSaCS+Fys+y6f8fe0aW8NTTtLdLVHGkyUGqhBHwcOLif8rRRLq/URHCptlV62JHVSfOhN7TfYB5oMlBqIVhXHZsTr/IXxUdxjJ2LdUQJRwJeti3NZW9fGj29fdB5SBNCFGkyUGoh+L2c6humZSyHNUUZ2gE6G+keNi/NY5l00fneXug+ojWEKJpzMhCRVhE5IiKHRKTWKvOIyG4RabJ+uq1yEZHviUiziNSJyJaI97nHWr9JRO6Za1xKxRVrFFG5zUd1Sa52gM6Gy03W1s+QsWQTu71FjKcXwVA3eFtjHdmiEK2awfXGmE3GmBrr+YPAK8aYauAV6znALUC19bgP+AGEkwfwTWA7sA345mQCUWoxCIXgyY48pGwLaRs/pZ3Hs+Vys+2qGxkeM5ypfRZ6m6Dtba0dRMF8NRPtAh6xlh8B7ogo/6kJ2wPkikgJcDOw2xjjNcb4gN3AznmKTamFFfDR9uajlPqPsauoJ9bRJLwr1i6nI3sTb3pzMSt3QkqaNrtFQTSSgQFeFpH9InKfVVZkjOmyls8ARdZyGXA6Ytt2q2ym8g8QkftEpFZEant69J9KJQi/l7rTXrz2fNaW6FXHc2WzCTdcuZ0j/Wmc6uzSW2JGSTSSwdXGmC2Em4DuF5FrI180xhjCCWPOjDEPGWNqjDE1BQV6q0CVGMZT3dR1DLKjZJzUlBQ9cUXB7Ves5U3n1fy4cyksvz6cYLWpaE7mnAyMMR3Wz7PA04Tb/Lut5h+sn5MzcnUAFRGbl1tlM5UrlfDe6pjglyM1lG+9Ve/lGyXpTge3X7GW55pG6Dv8ErS8piOL5mhOyUBEMkQka3IZuAk4CjwHTI4Iugd41lp+DviiNaroCuCc1Zz0EnCTiLitjuObrDKlElvAx9vvvkt2WgqXb71cE0EU3XNlJfm2Id4+cRayS3XOojlyzHH7IuBpEZl8r58bY14UkX3AEyJyL9AG3Gmt/zxwK9AM+IE/BTDGeEXkvwL7rPX+wRijf1WV2AI+xg4/xURTPX9V7sY5fj04NBlES1F2GtvXr6T22DFuXHOKdKdTm+DmYE7JwBjTAmycprwPuHGacgPcP8N7PQw8PJd4lIor1nTVpyZy+Xfl2eFvrVoziKovXLeRuw53sq7HxZ015bEOJ6HpFchKzZd0D4faB6hOG6CyIEu/tc6DtaXZXLMin8aDbxJo26f9BnOgyUCpeXJmYIQnu/IoWXsldp2uet48cIWHkbFxXulwwNiQzlk0S3PtM1BKTSfg4+jLj7BBzvCJgtRYR7OorVleyaqSHI7VH+QmE8CJAV+rjty6RFozUGoeTAz2sr+1l7ySKvLSHTrKZT653NR8/MscGivjHX8Z5FXryKJZ0GSgVLQFfOxrbGUkEODmJUG9QnYBrF22lMKVV/DqST/+ruMw3Ac2e6zDSiiaDJSKpoAP6p+h8cAbZKc6WLn5Om2uWCD33byV3aNr2d/cCc40OPGq9h1cAk0GSkWT34t3aIS3zqawoaoIR3quJoIFsqYkm1vWFvC7tgl6nEu0M/kSaTJQKprSPbzb1k8xXrZW5Wnz0AL79zfXYLDz2ttvQ9dhvQHOJdBkoFQUBezZ/LfWVZhlO8ituVNrBQusuKiE8qvu5pkzbk6mrNDO5EugyUCpaAn4eOnNP9AfGGfXDddqIoiRL964ib6sdTx37BzBnve0M/kiaTJQKhoCPoJHn6Zlz2/4av4RLi+SWEeUtNJS7PzFJy7nCV81dSfPaGfyRdJkoFQ0+L00dPioH87k2uVuJKDNErF0y/piti/N4dfNY3jTlmpn8kXQZKBUFISws//YCbZm9LG23K0dxzEmIvzV7R9hIiS89PpbmE7tTL4QTQZKzVXAx7E/PMPpQcN1yzKxr7hB+wviQEVZKatu/BN+0+OhLrgUskpgqBu8rbEOLS7p3ERKzVFoqI/Xjnfjz6pk9TIPhIKxDklZ7r52Ay/Xn+HFY79hZfAELqcdnFngqdSEPYXWDJSai4CPPzScoLd/gM+utGO3O7SJKI7YbMI377yKvROreKE7k1D1TiCk/QfT0JqBUrMV8DFx5Cn2vdVIYZaTjdtvgLwq/cYZZ5YVZHLnTdey78UWig8f5cqMDsBAdwMs/Qh49G8GWjNQava8JznS+B4nhx3sWFeB3ZWjJ5U49blrLmN89S5+eNxFm3NFuP+g8yC897J2KltmnQxEpEJEXhWRBhGpF5G/tMq/JSIdInLIetwasc03RKRZRBpF5OaI8p1WWbOIPDi3X0mpBRDwMdz8Jp0njnBndj2r81O1eSiOiQh/f+eV9OWu4+cHvQx3NgIG8pfrFcqWuTQTTQD/0RhzQESygP0istt67V+MMf89cmURWQvcBawDSoHfichK6+XvAx8D2oF9IvKcMaZhDrEpNb/8Xn7fPMDu8U38X1udyNIrtVYQ57LTUvinL1zLl34wwlBzgP+yKRfH4BkYC+gVysyhZmCM6TLGHLCWB4FjQNl5NtkFPG6MGTXGnASagW3Wo9kY02KMGQMet9ZVKj4FfDSd6uRQyxk+WumisGRpeHSKinvry3L42zuu4NHTeXz31ArMmD98hfKx30LHgaRuLopKn4GIVAKbgb1W0QMiUiciD4vI5NelMuB0xGbtVtlM5dPt5z4RqRWR2p6enmiErtSlCfgIHnmKl3c/T4bTwfU37tT7FSSYz9ZU8MD1K3ihrp03OkLaf2CZczIQkUzgV8DXjTEDwA+A5cAmoAv47lz3MckY85AxpsYYU1NQUBCtt1Xq4nlPcqC+keP9wsc2LCErx6OJIAH9x5tWsn3dSl6s7+F4/UHAQHZxUl+UNqdkICIphBPBo8aYpwCMMd3GmKAxJgT8iHAzEEAHUBGxeblVNlO5UvEl4KP76Kt0tRzhC+4G1helaadxghIR/svnrqKr/Gb+tq6E90w5nHwLepug+ZWkbDKay2giAX4MHDPG/I+I8pKI1T4JHLWWnwPuEpFUEakCqoF3gX1AtYhUiYiTcCfzc7ONS6l5EfAxdmo/zx5sZ49jG5dtqNFO4wSXlmLne1++AUo38zeH8jlhiqHqGug5DvXPwt4fga811mEumLnUDK4C/gS4Ycow0u+IyBERqQOuB/4KwBhTDzwBNAAvAvdbNYgJ4AHgJcKd0E9Y6yoVH6z7Gr/0+1fIH2rkTzekkpFXpp3Gi0B2Wgo/vXc7aYUr+LfDAVpamiE4CoG+cFI4/ETS1BDEGBPrGGalpqbG1NbWxjoMlQw6DlD36i/5UYPwqRV2rr/uBijdpLWCRaTfP8ZXH3qF0Z5mvr2mjdWpXkh3Q2YhlG9bVH9vEdlvjKmZWq5XICt1PgEf7XW/p7XxEHflNnDt6uJFdWJQYbnpTh766kdJW1rDvUfXsneoAJNRGO5DaH833GTUVQd9JxZtTUHnJlJqJgEffU17ear2FH2p2/nbTQ7s2k+waOW4Unjky9v4T79M4yuHM/i6fYg/KYeUobMw2A2dB6DyyvCsp4twOLHWDJSaTsDH8MEneeI3L1I50cxXt6aT7tF+gsUu1WHnf35uE3dft4H/pyGTX9aeZri/F1LSwlcpOzMX7fQVWjNQKlLAB34vQz1tPPPWAZoCHm7Zdi2lKy/X5qEkYbMJ37hlDZeV5fC9J0fwNb3NxzcUU5nSAcM9YATGhsNNRumL5zoTTQZKTbJGDY36z1H3h5fxDI/z18sKKcm7UhNBEvrEhlJWF3+av/5ZPk/WdvG5DZv4crARpx14+3tQsmFRNRlpM5FSk7wnGe7r4KfvnKZnOEjFxh2ULLsMtJ8gaa0ozOTR/+Nmrtl+BY8f7uP7e3y0jWRAKATYwlcsdx1ZFB3LWjNQyc1qFsJmp6/hNfbv3UPF+CibKj0Ul2X98RaJKmmlOx38w6717Fnh4vWnm/jZG/V8Iq+P1WNvk+awQ/cxKF4bbj5atROcGQnZfKTXGajkZTULMTZEe1sTzzaOcJpi/mKTULrhRsgpS8h/ajV/Bvt7+NkrB3hr/2Gucr7H1hVLqLEdx5GaDuMjEJqAss3hxLDxTnBXxjrkD5npOgOtGajk5fcSHBmg7kgdXR2tbHIKn9u+kvy8PCi5TJOA+pCs3AL+/NM3s/Pytbz69I/Y3dDFWOpZKguzKct3Yw+Ohq9NCIWg9hFYe1vC3FZTk4FKPgEfeFs509PNH96px9HfTn5hBTXrVpNatVU7i9UFLVtSQdVX/op9Dc08+kY9Fe1vktuXwk25o5S4A6RkusHbFJ4WOzULiteHf8ZxYtBkoJKHlQTG39vNifpaGrsH6bYVcf3qy1izfEm4f0ATgbpIku5mW83lXL61htcOX8kvXj/Ik+093NZbx9ZS4bJMBxnZxdD8ezj5BqRmQsEaWHFDXCYFTQYqOQR8mKNP09zSzNnGvbSPplOcl8cNa5eTddknIC1b+wfUrIgI129ayY6N1exr9fHoa+t4sekIl0sjazv3cVVGP/mF5aTYBdr3wehguJN51U4IBQETF8lBk4Fa3AI+TN9J6o4d5d2DR6jrd/Hp1BRuWGKnIC8Vsj3h0UKaBNQciQjbqjxsq7qedt92frP3GC8c2M8J7yjrulupygiRn+0iz9WLY+wEnN4LThfYU8M1horLw01J6Z5wkljgLyc6mkgtPlZz0PCgj8Z9r+BtrWNkJIArNYX8VVeyvjwXe8mGuG/DVYkvFDLsO97CH96t5ejJdjaPH6LK3kN+toulaQFyC4pxpWXAQBfYU8DhBGwfHKoa5drDTKOJNBmoxBZxnQB+H8GRAU4deg1v22HOnfMhE+OMpRexvDSPypISHNU7tF9AxcTYRIiDdYcYfPf/o6l7iPLxVoQQbqchzyW48sopdIFr7Bxk5oWHqo4Hpq89zCEx6NBStXhMJoCxYWh8kZGJEN6TB+gZGKF/cIjQ+BheRx5lOUWsyRwiO8OOpAYhu1ATgYoZp8PG9i1bYE0VNwz10dgToP5EK2+c6iStu5bKc6dJMRO4U8ZJyxgjKyOTfNsIGSlZpDrSoO3tcJ9DaiaUboHNn4/qZ1mTgYpfU771gwGbg6Ejv6Gj34/pOIDXH+LsiJATGmDIkYMnO58VaYPk5ThxpGVAQU1Uvk0pFTUuNzaXmzUFsGbtZQCMD32W5sYjNHQN8Fq3n8Ku39Pvm2C1zYe0t5PlCJHjBJNVQlaGjbJsHy6/V5OBWiQiT/ah4IdO+mPHfkvf8BgTnUcZGhlncGSc8ZEh+sccjOIkyxbAmZ7DKo+hKMNDTlYO9pRUKNimCUAllJRMD2u2XseayYLATfT3neFYt5+Ozk7aes6S0XMQz9kTBEN9uCvW4Er3RDWGuEkGIrIT+J+AHfhXY8y3YxySms7kCXxypMN0397TPeddDhnDudEQ48deYHh0AoeviTOOJaT3H2dkPEhgdAIz4WdgPHzSdzPAgC2LrNQUitLsFBbmk5caIic1G3v5lnnraFMqZlxucsvdfKQc2Lo1XBbwEew9SfdAgNzKtVH/jMdFMhARO/B94GNAO7BPRJ4zxjREfWezOHld0vLUb7jRfO8o7ccYQzDNTcjfRzBkmEhzY4b7CBmYcOZiAt5weWouoeAEQWyYQD8mMICjt57gRJDQmJ9u91Zyzu5jwoTIOvceoVCIUMhgQuOMh2xMhEIEgxOMBcPLoeAEo8bGxEQIF6MEjJNRnKQyxjHS2GLzEXRmke50UpA6RkF6PtmOCbLsTtIzsrDbBPAsiknBlLpkLjf2Cjel8/T2cZEMgG1AszGmBUBEHgd2AdFNBgEfTz/6vzhzto9VE40IBjA4TIgJsQMGuwkRFDvGWg6JHWMMNoKExIExBjtBQjgIGYPDKg8ZAwgnbEtZFmpFBDDgIMgEdgxY2820bCOEYCcYPvkiOJhgAjsg2JkgaC0b4D0qqOYUNgzGCI6I7ewECRrbB7czIASZMDaMAYeE4+I8y0aE46ElrLadRgjhNOMYBK9kk2OGyZKXaTMZ739775csxEChrR+f3YPDJhQxSCAlD2eqDU/IS8BZQIrDRlZoADLycDFOmt1wfXk1Lu8EQggQIHf64XUxGoOt1GIXL8mgDDgd8bwd2D51JRG5D7gPYMmSJZe+F78Xj8tBIDeX/MEgo/Z0jIGcoJeBlDwwkD3Rx2BKHgZr2ZEH7y9nYzDvL4fLexlIyUaA1OAQE+nplPph1J4BQNZ4L4Mp+eF1x3sZdOZ/oFwQMid6GHLkg7U87MjHYLOWc0GEzPEehhy5YBNSJ4ZwZuRQOmxjzJGJAbLGehlyhk+OmWO9DFvLGWO9DDs9CJA53os/NQ8BMsZ7CaR6piwL6WO9jKTlYwNSJoZYnleMx3eWoCML40glPdCBpNgwaWXYx4awubKxT4xiD6Vic2Zis4E95AZ7OHERzI5Ydv1xmfQPf8OfWtPRk75SCyZeksFFMcY8BDwE4esMLvkN0j1ct6oYxoagvTh8L9OpJ6lg2iUup0ac4LJYV1wM3X0R7506/bofWC6bYbl0+mUyWF1cCN09EfspiVj3YpaLI5aLIpYLI/aTR2lxNtjc4f3YbVBxTXjZmRaenTG/evpv7xfb1DX1ZB+HU/4qlQziJRl0ABURz8utsuhyucO3qPN7Yc1t89eWv+FzC9NnEIv9eKrCxzJyFNB0394jT+ozLSul4kZcXIEsIg7gPeBGwklgH/B5Y0z9TNvoFchKKXXp4voKZGPMhIg8ALxEeGjpw+dLBEoppaIrLpIBgDHmeeD5WMehlFLJyBbrAJRSSsWeJgOllFKaDJRSSmkyUEopRZwMLZ0NEekB2ma5eT7QG8VwoiVe44L4jU3jujQa16WL19hmG9dSY0zB1MKETQZzISK1042zjbV4jQviNzaN69JoXJcuXmOLdlzaTKSUUkqTgVJKqeRNBg/FOoAZxGtcEL+xaVyXRuO6dPEaW1TjSso+A6WUUh+UrDUDpZRSETQZKKWUSr5kICI7RaRRRJpF5MEYxlEhIq+KSIOI1IvIX1rl3xKRDhE5ZD1ujUFsrSJyxNp/rVXmEZHdItJk/VzQ24+JyKqIY3JIRAZE5OuxOl4i8rCInBWRoxFl0x4jCfue9ZmrE5EtCxzXP4vIcWvfT4tIrlVeKSKBiGP3wwWOa8a/nYh8wzpejSJy8wLH9YuImFpF5JBVvpDHa6bzw/x9xowxSfMgPD32CWAZ4AQOA2tjFEsJsMVaziJ8P4e1wLeA/xTj49QK5E8p+w7woLX8IPBPMf47ngGWxup4AdcCW4CjFzpGwK3AC4RvH3cFsHeB47oJcFjL/xQRV2XkejE4XtP+7az/g8NAKlBl/c/aFyquKa9/F/j7GByvmc4P8/YZS7aawTag2RjTYowZAx4HdsUiEGNMlzHmgLU8CBwjfC/oeLULeMRafgS4I3ahcCNwwhgz2yvQ58wY8wbgnVI80zHaBfzUhO0BckWkZKHiMsa8bIyZsJ7uIXwnwQU1w/GayS7gcWPMqDHmJNBM+H93QeMSEQHuBB6bj32fz3nOD/P2GUu2ZFAGnI543k4cnIBFpBLYDOy1ih6wqnoPL3RzjMUAL4vIfhG5zyorMsZ0WctngKIYxDXpLj74Dxrr4zVppmMUT5+7LxP+BjmpSkQOisjrInJNDOKZ7m8XL8frGqDbGNMUUbbgx2vK+WHePmPJlgzijohkAr8Cvm6MGQB+ACwHNgFdhKupC+1qY8wW4BbgfhG5NvJFE66XxmRMsog4gduBX1pF8XC8PiSWx2gmIvJ3wATwqFXUBSwxxmwG/gPwcxHJXsCQ4vJvF+FuPvilY8GP1zTnh/dF+zOWbMmgA6iIeF5ulcWEiKQQ/kM/aox5CsAY022MCRpjQsCPmKfq8fkYYzqsn2eBp60YuierndbPswsdl+UW4IAxptuKMebHK8JMxyjmnzsR+RLwCeAL1kkEqxmmz1reT7htfuVCxXSev108HC8H8CngF5NlC328pjs/MI+fsWRLBvuAahGpsr5h3gU8F4tArPbIHwPHjDH/I6I8sp3vk8DRqdvOc1wZIpI1uUy48/Eo4eN0j7XaPcCzCxlXhA98W4v18ZpipmP0HPBFa8THFcC5iKr+vBORncDfALcbY/wR5QUiYreWlwHVQMsCxjXT3+454C4RSRWRKiuudxcqLstHgePGmPbJgoU8XjOdH5jPz9hC9IzH04Nwr/t7hLP638UwjqsJV/HqgEPW41bgZ8ARq/w5oGSB41pGeCTHYaB+8hgBecArQBPwO8ATg2OWAfQBORFlMTlehBNSFzBOuH323pmOEeERHt+3PnNHgJoFjquZcHvy5Ofsh9a6n7b+xoeAA8BtCxzXjH874O+s49UI3LKQcVnlPwG+NmXdhTxeM50f5u0zptNRKKWUSrpmIqWUUtPQZKCUUkqTgVJKKU0GSiml0GSglFIKTQZKKaXQZKCUUgr4/wFKo0OidZRvYQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S, I, R, _ = odeint(func, y0, tspan, method='euler').T\n",
    "pred = func.params[0] * S * I / y0[-1]\n",
    "plt.plot(pred.detach())\n",
    "plt.plot(data.clone().detach(), '.', 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
}