Poisson-weighted-convolutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An infinite sum of Poisson-weighted convolutions\n",
    "\n",
    "Kyle Cranmer, Aug 2018\n",
    "\n",
    "If viewing on GitHub, this looks better with nbviewer: [click here](http://nbviewer.jupyter.org/github/cranmer/poisson-convolution-sum/blob/master/Poisson-weighted-convolutions.ipynb)\n",
    "\n",
    "Consider a variable $x = \\sum_{i=1}^n z_i$, where the $z_i \\sim p(z|\\theta)$ are iid and $n \\sim Pois(\\nu)$. For each value of $n$, the distribution $p(x|n)$ is given by repeated convolutions\n",
    "\\begin{equation}\n",
    "{p(x|n, \\theta)} = p(z|\\theta)\\underbrace{\\star \\dots \\star}_\\text{n times} p(z|\\theta)\n",
    "\\end{equation}\n",
    "and for a Poisson mean of $\\nu$ the expected distribution for $x$ is  \n",
    "\\begin{equation}\n",
    "{p(x|\\nu, \\theta)} = \\sum_{n=0}^\\infty e^{-\\nu} \\frac{\\nu^n}{n!} p(x|n, \\theta)\n",
    "\\end{equation}\n",
    "\n",
    "The infinite tower of Poisson-weighted convolutions can be computed efficiently with a nice trick documented in *Analytic Confidence Level Calculations using the Likelihood Ratio and Fourier Transform* by Hongbo Hu and Jason Nielsen https://arxiv.org/pdf/physics/9906010.pdf. See also [this old paper](https://arxiv.org/abs/physics/0312050) and this [code](http://phystat.org/phystat/packages/0703002.1.html) for a C++ implementation.\n",
    "\n",
    "First we take advantage of the convolution theorem relating convolutions to multiplication in the Fourier domain, denoted with a bar.\n",
    "\\begin{equation}\n",
    "\\overline{p(x|n, \\theta)} = \\overline{p(z|\\theta)\\underbrace{\\star \\dots \\star}_\\text{n times} z|\\theta} = \\left[\\overline{p(z|\\theta)} \\right]^n \n",
    "\\end{equation}\n",
    "and then we can compress the infinite sum into the following exponential\n",
    "\\begin{equation}\n",
    "\\overline{p(x|\\nu, \\theta)} = \\sum_{n=0}^\\infty e^{-\\nu} \\frac{\\nu^n}{n!}  \\overline{p(x|n)}  =  \\exp\\left( \\nu \\left[ \\overline{p(z|\\theta)}-1 \\right] \\right)\n",
    "\\end{equation}\n",
    "\n",
    "In the code below, we will reimplement the technique using pytorch. An interesting feature of implementing this in pytorch is that we can backprop through the inverse FFT, the exponentiation, the multplication, the subtraction, and forward FFT to calculate the gradient\n",
    "\\begin{equation}\n",
    "\\nabla_\\theta\\, p(x|\\nu, \\theta)\n",
    "\\end{equation}\n",
    "and \n",
    "\\begin{equation}\n",
    "\\nabla_\\nu \\, p(x|\\nu, \\theta)\n",
    "\\end{equation}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classes for Poisson convolutions with complex FFT "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def complex_exp_pois(z, nu):\n",
    "    real1, imag1 = z.t()\n",
    "    real1=real1-1.\n",
    "    real_out = torch.exp(real1*nu)*torch.cos(nu*imag1)\n",
    "    imag_out = torch.exp(real1*nu)*torch.sin(nu*imag1)\n",
    "    return torch.stack([real_out, imag_out], dim = -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This defines p(z|\\theta)\n",
    "class base_dist():\n",
    "    def __init__(self):\n",
    "        self.theta = torch.ones(1,1, requires_grad=True, dtype=torch.double)\n",
    "        self.sigma = torch.tensor(0.051,requires_grad=True, dtype=torch.double)        \n",
    "    def forward(self, x):\n",
    "        max_ = torch.max(x)\n",
    "        min_ = torch.min(x)\n",
    "        nbins = x.size()[1]\n",
    "        norm = torch.tensor(max_/nbins/np.sqrt(2.*np.pi), dtype=torch.double)/self.sigma\n",
    "        y = torch.empty(1,nbins,2, dtype=torch.double)\n",
    "        y[0,:,0] = (torch.exp(-(x-self.theta)*(x-self.theta)/2./self.sigma/self.sigma))*norm\n",
    "        y[0,:,1]= torch.zeros(nbins)\n",
    "        return y\n",
    "    def parameters(self):\n",
    "        return [self.theta, self.sigma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This does Poisson convolutions of the base distribution p(z|\\theta)\n",
    "class pois_conv():\n",
    "    def __init__(self):\n",
    "        self.nu = torch.tensor(3.,requires_grad=True, dtype=torch.double)\n",
    "        self.base_dist = base_dist()\n",
    "    def forward(self, x):\n",
    "        y = self.base_dist.forward(x)\n",
    "        yfft = torch.fft(y,1,normalized=False)\n",
    "        y2fft = complex_exp_pois(yfft[0], nu=self.nu).view(-1,nbins,2)\n",
    "        y2 = torch.ifft(y2fft,1, normalized=False)\n",
    "        return y2\n",
    "    def parameters(self):\n",
    "        temp = list(self.base_dist.parameters())\n",
    "        temp.append(self.nu)\n",
    "        return temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check with simple Gaussian $p(z|\\theta)$ base distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_=0\n",
    "max_=15\n",
    "nbins=3000\n",
    "\n",
    "x = torch.arange(min_, max_,(max_-min_)/nbins, dtype=torch.double)\n",
    "x=x.view(-1,nbins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[tensor([[ 1.]], dtype=torch.float64), tensor(1.00000e-02 *\n",
       "        5.1000, dtype=torch.float64), tensor(3., dtype=torch.float64)]"
      ]
     },
     "execution_count": 239,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y2_ = pois_conv()\n",
    "y2_.parameters()\n",
    "y_ = y2_.base_dist\n",
    "y2_.forward(x)\n",
    "y2_.parameters()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 2, 6]"
      ]
     },
     "execution_count": 244,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#quick helper to for chekcing Poisson probabilities\n",
    "def myfact(xin):\n",
    "    return [np.math.factorial(x) for x in xin]\n",
    "myfact(np.arange(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Poisson probs\n",
    "nu = y2_.nu\n",
    "xint = np.arange(min_, max_)\n",
    "pois = np.exp(-nu.item())*np.power(nu.item(),xint)/myfact(xint)\n",
    "scaled_pois = pois\n",
    "#heights of Gaussian bumps will be lower by width of Gaussian, which scales like sqrt(x()\n",
    "scaled_pois[1:]= scaled_pois[1:]/np.sqrt(xint[1:])/np.sqrt(2*np.pi)/y_.sigma.item()*max_/nbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXmYVOWZ8P27a+19YZFdwYggKIsiMcEkKjHgRIPj6BvNFaOj+Zw4yTdvJhMTnXmTvPqOX3ScjNnMzJjRyCQmaIxDeCcqWUjG0bih4AKKIgh0g9DQ9F7VtT3fH+ec6urqWk6dU0hD37/r4qLq1HPuOqcbnvvcuxhjUBRFUZRiBI72BSiKoiijG1UUiqIoSklUUSiKoiglUUWhKIqilEQVhaIoilISVRSKoihKSVRRKIqiKCVRRaEoiqKURBWFoiiKUpLQ0b6AajBhwgQzc+bMo30ZiqIoxxQvvvjiQWPMxHLrjgtFMXPmTDZu3Hi0L0NRFOWYQkR2uVmnridFURSlJKooFEVRlJKoolAURVFKclzEKBRFeW9IJpO0tbURj8eP9qUoFVBTU8P06dMJh8OezldFoSiKa9ra2mhsbGTmzJmIyNG+HMUFxhgOHTpEW1sbs2bN8iTDletJRFaKyDYR2S4iNxf4PCoiD9mfPyciM3M+u8U+vk1EVuSdFxSRTSLynznHZtky3rJlRjzdmaIoVScejzN+/HhVEscQIsL48eN9WYFlFYWIBIF7gIuAecBVIjIvb9n1wGFjzCnA3cCd9rnzgCuB+cBK4Ae2PIf/CbyeJ+tO4G5jzGzgsC1bUZRRgiqJYw+/vzM3FsVSYLsxZocxJgGsAVblrVkFrLZfPwIsF+vKVgFrjDGDxpidwHZbHiIyHfg48G+OEPucC2wZ2DIv9XJjiqIoSnVwoyimAXty3rfZxwquMcakgG5gfJlzvw18BcjkfD4e6LJlFPsuRVGUo8bmzZt57LHHKj7vvPPOO2YLg90oikI2i3G5puBxEbkYOGCMedHDd1kLRW4QkY0isrGjo6PQEkVRlKrjVVEcy7hRFG3AjJz304G9xdaISAhoBjpLnLsM+ISIvIPlyrpARH4CHARabBnFvgsAY8y9xpglxpglEyeWbVWiKMpxwr//+7+zYMECFi5cyNVXX82uXbtYvnw5CxYsYPny5ezevRuAa6+9lhtvvJHzzz+fk08+mf/6r//iuuuu47TTTuPaa6/NymtoaOBv/uZvOPPMM1m+fDnOg2euBXDw4EFmzpxJIpHg61//Og899BCLFi3ioYceor+/n+uuu46zzz6bxYsX88tf/hKAWCzGlVdeyYIFC/jkJz9JLBZ7b39QVcRNeuwLwGwRmQW0YwWnP5W3Zh1wDfAMcDmwwRhjRGQd8FMR+SdgKjAbeN4Y8wxwC4CInAd82Rjzafv9720Za2yZv/R1h4qiHBFu/b9b2Lq3p6oy501t4huXzC/6+ZYtW7j99tt5+umnmTBhAp2dnVxzzTV85jOf4ZprruH+++/nr/7qr1i7di0Ahw8fZsOGDaxbt45LLrmEp59+mn/7t3/j7LPPZvPmzSxatIj+/n7OPPNMvvWtb3Hbbbdx66238v3vf7/g90ciEW677TY2btyYXfO3f/u3XHDBBdx///10dXWxdOlSPvrRj/Kv//qv1NXV8corr/DKK69w5plnVvVn9V5S1qKw4wVfANZjZSg9bIzZIiK3icgn7GX3AeNFZDvwJeBm+9wtwMPAVuAJ4PPGmHSZr/wq8CVb1nhbtqIoChs2bODyyy9nwoQJAIwbN45nnnmGT33Kena9+uqreeqpp7LrL7nkEkSEM844g0mTJnHGGWcQCASYP38+77zzDgCBQIBPfvKTAHz6058edr4bfv3rX3PHHXewaNEizjvvPOLxOLt37+bJJ5/k05/+NAALFixgwYIFfm//qOGq4M4Y8xjwWN6xr+e8jgNXFDn3duD2ErL/APwh5/0O7MwoRVFGL6We/I8UxpiyqZ65n0ejUcBSBs5r530qlRpxbu75oVCITMbKtSlVg2CM4Re/+AVz5swpeS3HMtrrSVGUY4bly5fz8MMPc+jQIQA6Ozv54Ac/yJo1awB48MEHOffccyuSmclkeOQRKyP/pz/9afb8mTNn8uKLVr6N8zlAY2Mjvb292fcrVqzge9/7HsZYeTebNm0C4MMf/jAPPvggAK+99hqvvPJKxfc7WlBFoSjKMcP8+fP5u7/7Oz7ykY+wcOFCvvSlL/Hd736XH/3oRyxYsIAf//jHfOc736lIZn19PVu2bOGss85iw4YNfP3rlrPky1/+Mv/8z//MBz/4QQ4ePJhdf/7557N169ZsMPtrX/sayWSSBQsWcPrpp/O1r30NgBtvvJG+vj4WLFjAP/zDP7B06bHrKBFHCx7LLFmyxByr+cmKcizx+uuvc9pppx3ty6gqDQ0N9PX1He3LOOIU+t2JyIvGmCXlzlWLQlEURSmJKgpFUcY0Y8Ga8IsqCkVRFKUkqigURVGUkqiiUBRFUUqiikJRFEUpiSoKRVHGNA0NDQDs3buXyy+/vOTab3/72wwMDFQk/w9/+AMXX3yx5+tzw8yZM4fVelQbVRSKohx3pNPlWsqNZOrUqcMqsAvhRVEcD6iiUBTliLF2UzvL7tjArJt/xbI7NrB2U7svee+88w5z587lmmuuYcGCBVx++eXZjXvmzJncdtttnHvuufz85z/n7bffZuXKlZx11ll86EMf4o033gBg586dfOADH+Dss8/OVlE7sk8//XTAUjRf/vKXOeOMM1iwYAHf+973+O53v8vevXs5//zzOf/88wGrIeAHPvABzjzzTK644opsqu0TTzzB3LlzOffcc3n00UcL3ssDDzzAqlWrWLlyJXPmzOHWW2/NfvaTn/yEpUuXsmjRIv7iL/4iq/huvPFGlixZwvz58/nGN74xQmYsFmPlypX88Ic/9PVzHoEx5pj/c9ZZZxlFUY48W7dudb32P15qM3P/1+PmpK/+Z/bP3P/1uPmPl9o8f//OnTsNYJ566iljjDF//ud/bu666y5jjDEnnXSSufPOO7NrL7jgAvPmm28aY4x59tlnzfnnn2+MMeaSSy4xq1evNsYY8/3vf9/U19dnZc+fP98YY8wPfvADc9lll5lkMmmMMebQoUPZ7+jo6DDGGNPR0WE+9KEPmb6+PmOMMXfccYe59dZbTSwWM9OnTzdvvvmmyWQy5oorrjAf//jHR9zLj370IzN58mRz8OBBMzAwYObPn29eeOEFs3XrVnPxxRebRCJhjDHmxhtvzF6vcx2pVMp85CMfMS+//HL2unbu3GmWL1+eXZtPod8dsNG42GPVolAU5Yhw1/ptxJLDXUCxZJq71m/zJXfGjBksW7YMGNkW3GkX3tfXxx//+EeuuOKK7FP5vn37AHj66ae56qqrAKsteSF++9vf8rnPfY5QyGqwPW7cuBFrnn32WbZu3cqyZctYtGgRq1evZteuXbzxxhvMmjWL2bNnIyLZVuOFuPDCCxk/fjy1tbVcdtllPPXUU/zud7/jxRdf5Oyzz2bRokX87ne/Y8eOHQA8/PDDnHnmmSxevJgtW7awdevWrKxVq1bx53/+53zmM59x/bN0i6s244qiKJWyt6vwRLdix92S37o79319fT1gdYRtaWlh8+bNrmTkY1y0MzfGcOGFF/Kzn/1s2PHNmze7bi9e6F6MMVxzzTV885vfHPbZzp07+cd//EdeeOEFWltbufbaa4e1P1+2bBmPP/44n/rUp6re3lwtCkVRjghTW2orOu6W3bt388wzzwDws5/9rGBb8aamJmbNmsXPf/5zwNrUX375ZcDaUHPbkhfiYx/7GP/yL/+SnVnR2dkJDG8xfs455/D000+zfft2AAYGBnjzzTeZO3cuO3fu5O23385eYzF+85vf0NnZSSwWY+3atSxbtozly5fzyCOPcODAgex379q1i56eHurr62lubmb//v08/vjjw2TddtttjB8/nr/8y78s9yOsGFUUiqIcEW5aMYfacHDYsdpwkJtWjBzwUwmnnXYaq1evZsGCBXR2dnLjjTcWXPfggw9y3333sXDhQubPn5+dZf2d73yHe+65h7PPPpvu7u6C5372s5/lxBNPzM7m/ulPfwrADTfcwEUXXcT555/PxIkTeeCBB7jqqqtYsGAB55xzDm+88QY1NTXce++9fPzjH+fcc8/lpJNOKnov5557LldffTWLFi3iz/7sz1iyZAnz5s3j7//+7/nYxz7GggULuPDCC9m3bx8LFy5k8eLFzJ8/n+uuuy7rfsvl29/+NvF4nK985SuV/lhL4qrNuIisBL4DBIF/M8bckfd5FPh34CzgEPBJY8w79me3ANcDaeCvjDHrRaQGeBKIYrm/HjHGfMNe/wDwEcD5DV5rjClsP9pom3FFeW+otM342k3t3LV+G3u7YkxtqeWmFXO4dPE0z9//zjvvcPHFF/Paa695ljFaeOCBB4bN3j7S+GkzXjZGISJB4B7gQqANeEFE1hljtuYsux44bIw5RUSuBO4EPiki84ArgfnAVOC3InIqMAhcYIzpE5Ew8JSIPG6MedaWd5MxpnRCs6Ioo55LF0/zpRiU0YEb19NSYLsxZocxJgGsAVblrVkFrLZfPwIsFyuasgpYY4wZNMbsBLYDS+3MLKe3b9j+c+xPUFIU5Ygyc+bM48KaALj22mvfM2vCL24UxTRgT877NvtYwTXGmBSW22h8qXNFJCgim4EDwG+MMc/lrLtdRF4Rkbttt5aiKKMEN+5qZXTh93fmRlEUyrPK/9Zia4qea4xJG2MWAdOBpSJyuv35LcBc4GxgHPDVghclcoOIbBSRjR0dHeXvQlEU39TU1HDo0CFVFscQxhgOHTpETU2NZxlu6ijagBk576cDe4usaRORENAMdLo51xjTJSJ/AFYCrxlj9tkfDYrIj4AvF7ooY8y9wL1gBbNd3IeiKD6ZPn06bW1t6MPZsUVNTQ3Tp0/3fL4bRfECMFtEZgHtWMHpT+WtWQdcAzwDXA5sMMYYEVkH/FRE/gkrmD0beF5EJgJJW0nUAh/FCoAjIlOMMfvsGMelwPHhkFSU44BwOMysWbOO9mUo7zFlFYUxJiUiXwDWY6XH3m+M2SIit2H1CVkH3Af8WES2Y1kSV9rnbhGRh4GtQAr4vDEmLSJTgNV2RlUAeNgY85/2Vz5oKxIBNgOfq+YNK4qiKJXhqo5itKN1FIqiKJXjto5CK7MVRVGUkqiiUBRFUUqiikJRFEUpiSoKRVEUpSSqKBRFUZSSqKJQFEVRSqKKQlEURSmJKgpFURSlJKooFEVRlJKoolAURVFKoopCURRFKYkqCkVRFKUkqigURVGUkqiiUBRFUUqiikJRFEUpiSoKRVEUpSSqKBRFUZSSuFIUIrJSRLaJyHYRubnA51ERecj+/DkRmZnz2S328W0issI+ViMiz4vIyyKyRURuzVk/y5bxli0z4v82FUVRFK+UVRT2XOt7gIuAecBVIjIvb9n1wGFjzCnA3cCd9rnzsOZnzwdWAj+w5Q0CFxhjFgKLgJUico4t607gbmPMbOCwLVtRFEU5SrixKJYC240xO4wxCWANsCpvzSpgtf36EWC5iIh9fI0xZtAYsxPYDiw1Fn32+rD9x9jnXGDLwJZ5qcd7UxRFUaqAG0UxDdiT877NPlZwjTEmBXQD40udKyJBEdkMHAB+Y4x5zj6ny5ZR7LsURVGU9xA3ikIKHDMu1xQ91xiTNsYsAqYDS0XkdJffZX2hyA0islFENnZ0dBS9eEVRFMUfbhRFGzAj5/10YG+xNSISApqBTjfnGmO6gD9gxTAOAi22jGLf5Zx3rzFmiTFmycSJE13chqIoiuIFN4riBWC2nY0UwQpOr8tbsw64xn59ObDBGGPs41faWVGzgNnA8yIyUURaAESkFvgo8IZ9zu9tGdgyf+n99hRFURS/hMotMMakROQLwHogCNxvjNkiIrcBG40x64D7gB+LyHYsS+JK+9wtIvIwsBVIAZ83xqRFZAqw2s6ACgAPG2P+0/7KrwJrROTvgU22bEVRFOUoIdZD/LHNkiVLzMaNG4/2ZSiKohxTiMiLxpgl5dZpZbaiKIpSElUUiqIoSknKxigU5WiwdlM7d63fxt6uGFNbarlpxRwuXawlNYpyNFBFoYw61m5q55ZHXyWWTAPQ3hXjlkdfBVBloShHAXU9KaOOu9ZvyyoJh1gyzV3rtx2lK1KUsY0qCmXUsbcrVtFxRVGOLKoolFHH1Jbaio4rinJkUUWhjDpuWjGHUGB426/acJCbVsw5SlekKGMbVRTKqOPSxdNYMrM1+35aSy3fvOwMDWQrylFCs56UUcn4higAkVCAp2++4ChfjaKMbdSiUEYlvXFrJEkilSGRyhzlq1GUsY0qCmVU0htPZl/3D6ZKrFQU5UijrielKlS7ktqxKAD6Eyla6yPVuExFUTygikLxzZGopO6NJ2mMhugdTDGQSJc/QVGUI4a6nhTfHIlK6t54iknNNUB1XE9rN7Wz7I4NzLr5Vyy7YwNrN7X7lqkoYwVVFIpvql1Jnc4YBhJpJjVZmU/9g/4sCsfiae+KYRiyeFRZKIo7VFEovql2JfVgylIM4+ttRZHwZ1Fo7yhF8YcrRSEiK0Vkm4hsF5GbC3weFZGH7M+fE5GZOZ/dYh/fJiIr7GMzROT3IvK6iGwRkf+Zs/5/i0i7iGy2//yJ/9tUjiQ3rZhDUKpXSe2kw46zA9gDPhWF9o5SFH+UVRT2XOt7gIuAecBVIjIvb9n1wGFjzCnA3cCd9rnzsOZnzwdWAj+w5aWAvzHGnAacA3w+T+bdxphF9p/HfN2hcsS5dPE0Tp/WlH3vt5J60FYUzbVhwL/rSXtHKYo/3FgUS4HtxpgdxpgEsAZYlbdmFbDafv0IsFxExD6+xhgzaIzZCWwHlhpj9hljXgIwxvQCrwNHrT+DBjr901hjbeqnTmrg6Zsv8JUa61gUTbaiGPRZcHfTijmEtXeUonjGjaKYBuzJed/GyE09u8YYkwK6gfFuzrXdVIuB53IOf0FEXhGR+0WklSOIBjqrQ3csOexvPzgxiqaa0LD3Xrl08TTOm3NC9v3U5hrtHaUoFeBGUUiBY8blmpLnikgD8Avgi8aYHvvwPwPvAxYB+4BvFbwokRtEZKOIbOzo6Ch9ByXQQGd16LErqatR8+BYEI6VEk/6b+HRUhfOvn70L5epklCUCnCjKNqAGTnvpwN7i60RkRDQDHSWOldEwlhK4kFjzKPOAmPMfmNM2hiTAX6I5foagTHmXmPMEmPMkokTJ7q4jcJooLM6OJXUsSooCsf1FA0HiIQCvi0KgIN9g9nXue1BFEUpjxtF8QIwW0RmiUgEKzi9Lm/NOuAa+/XlwAZjjLGPX2lnRc0CZgPP2/GL+4DXjTH/lCtIRKbkvP1T4LVKb6oSNNBZHQZtqyyVMb6b+DkWRTQUIBoKMFgFiyK3JUiPKgpFqYiyisKOOXwBWI8VdH7YGLNFRG4TkU/Yy+4DxovIduBLwM32uVuAh4GtwBPA540xaWAZcDVwQYE02H8QkVdF5BXgfOCvq3WzhbhpxRxqw8FhxzTQWTmJdIZoyPrnlO/Kq1hWjqKoCQerYlH0DaaY2GjVZVQjjqIoYwlXvZ7sFNXH8o59Ped1HLiiyLm3A7fnHXuKwvELjDFXu7mmauH4qv/2P15lIJFmWhUa2o01MhlDMm2Y1BRhf88gsUQ6m9rqBceiiASDVbMo+gZTTG2ppaN3kJ6YdqNVlErQymwsZfGJhVOZ1BT1ndo5FkmkrY28pbY6BXK5MYqacJB4FSyK/sEU01qs3lHqelKUylBFkYPJz+VSXJEtkLMzi3y7ntLW+ZFg9WIUfYMppjRbcacedT0pSkWoorCRgo4wxQ2OBdDqKAqfmU+OYqiWRTGYSpNMG1rrwoQCQr+2LVeUilBFofjGcT05cQm/tRSOvGpZFE4LkIZoiNpIsCopvIoyllBFkYN6nrzhpMY2ZQvkqmNRRELVsSiceRb10RB1qigUpWJUUWRR35NXHAvAqaR23vuVFw1VJ+vJqaFoiIaoi4QY8KnIFGWsoYoiBw1meyORbblh92byubE7Fko4KHYdhU/XU2LIoqgJB4n5zMpSlLGGKgobDWZ7J19R+LUoBu3iPREhGgr4dmX12a6nhhrL9aQzuBWlMlRRKL4ZHGFR+K/MjthV3tWwKAbsYHZ9xI5RqOtJUSpCFcUw1PfkhUSqujGKwdRQO5BqWBTO+TXhALVhDWYrSqWoorBRz5N3RloUPoPZqQzRkNV/KxoKMJjKYHwEkIaaDAapVdeTolSMKoocNJjtDceCqIuECEh1LArH9RS1Gzb6kZlrUajrSVEqRxWFjQazvePEJKIhZ36EX4sinXU9RYIB+5gPRZFyFEWQ2nBIXU+KUiGqKBTfZCupQwGioaDveRSJYRaF9bcf5ZMt4AsG7KynlC9XlqKMNVRR5KBbhzdy50dUYyLdYCqTtSSqZVFEQgECAaE2EiRj/CkeRRlrqKKwEQ1ne8bZxCNVmkiXSGWylkS1LIoa20JxhlT5zaRSlLGEq8FFYwV1R3hjaNCQbVFUIZjtZFBFgnYw24+iSKWzQfHaiPV3LJmmxcc1rt3Uzl3rt7G3K8ZUHXalHOe4sihEZKWIbBOR7SJyc4HPoyLykP35cyIyM+ezW+zj20RkhX1shoj8XkReF5EtIvI/c9aPE5HfiMhb9t+t/m/TzT2+F98yOli7qZ1ld2xg1s2/YtkdG1i7qd2XvEQqQ0AgFLRiFNVOjwV8ubPiyQw14aG6DPCXwrt2Uzu3PPoq7V0xDNDeFeOWR1/1/XNUlNFKWUUhIkHgHuAiYB5wlYjMy1t2PXDYGHMKcDdwp33uPOBKYD6wEviBLS8F/I0x5jTgHODzOTJvBn5njJkN/M5+r1SJI7HJWfOyrY09EgpUpSmgE8x2/vYVo0imqbGvr8a2LPy4su5av21Eim0smeau9ds8y1SU0Ywbi2IpsN0Ys8MYkwDWAKvy1qwCVtuvHwGWi4jYx9cYYwaNMTuB7cBSY8w+Y8xLAMaYXuB1YFoBWauBS73dWuWMBcfTkdjkhmUphQK+W3gMJtPD5IHPGEUqk1UQjjw/MYq9XbGKjivKsY4bRTEN2JPzvo2hTX3EGmNMCugGxrs513ZTLQaesw9NMsbss2XtA04odFEicoOIbBSRjR0dHS5uozRjxfN0JDa5wdTwjb0aFkW0yhbFUEsQ/xbF1Jbaio4ryrGOG0VRaA/Nf/gutqbkuSLSAPwC+KIxpsfFtQwJMeZeY8wSY8ySiRMnVnJqCZlVETOqORKbXG46azWyngaTuRaKs7H7iVGksxZFTdh/zOOmFXMIBYb/064NB7lpxRzPMhVlNONGUbQBM3LeTwf2FlsjIiGgGegsda6IhLGUxIPGmEdz1uwXkSn2minAAbc34wcZI9Hsm1bMyaaIOvjd5BKp4RZAddqMD8U8oBqup+GKJ+5DmV26eBpnnjiUYzG+PsI3LztDs56U4xY3iuIFYLaIzBKRCFZwel3emnXANfbry4ENxso1XQdcaWdFzQJmA8/b8Yv7gNeNMf9UQtY1wC8rvSmlOJcunsY3LzuDcNBSjLXhoO9NbniMIujrad0YMyLmAf4UheV6smMUVbAowFJgExqiAHx15VxVEspxTVlFYcccvgCsxwo6P2yM2SIit4nIJ+xl9wHjRWQ78CXsTCVjzBbgYWAr8ATweWNMGlgGXA1cICKb7T9/Ysu6A7hQRN4CLrTfvyeMlTqKSxdP48RxdQAsnTXO9yY3LKYQDPiKJwyNQR2uKPzFKIYK+JzsJ7/usa5YglMnNQDQ0TfoS5aijHZcFdwZYx4DHss79vWc13HgiiLn3g7cnnfsKYrEj40xh4Dlbq5L8Y7zhN4VS/qXlRyezupLUaTyFYX/4POwrCdbYcR9WhRdA0lOPaGRhmiIjl5VFMrxjbbwyGFs2BMWzsZbjfnR+XUP/qqoh9qB5P7tS+awrCf/BXcA3bEkzXVhWurC9FRB2SrKaEYVhc0YiWVnceoIqjHEJ7eS2m+b8XyLYiiY7SPrKZWb9eTfQkmlM/TGUzTXhmmsCdMTV0WhHN+ooshlDJkUQxZFdRRFbrfXVMaQyXj7YSbyLIpgQAgFxLNFkc4YkmmTjU041+mn4K4nbllhLbVhmmpC2feKcryiimIM4mQWAVWZ9pZbcJd1FXlMkR1qMDiUwhv1YaU4logTmwgEhEjQn9XTNZAAoKUuQmNNmF5VFMpxjioKm7HUZjy322ssmfad7VXNdNZ81xNY41C9WhROvUTNMHn+ZmZ02zGJ5rowTbUhjVEoxz2qKHIYK54nJ5DbUhfGGH/FZzA8PdZvOquzgUdyNnbLAvC2sQ/mjEF1iIaCvu7ZyRRrrg3TVBOmV2MUynGOKgqbsRTMdjbP1roIAAM+M58GUwW6vXp0PRW2KLxnUjkKwXE9ObL9WBSOq6mpJkRTTYjewZTnmIyiHAuooshhrBTcxXMsCvAfpyioKDxbFMOD2YCvmIITtHaC2WD1e/IToxgYtBRFfTREY41llfVXIc1YUUYrqijGIPkWhZ/MJycwHs1mPfmbSFdIUfixKBx5+a4nP63Q++2fV104RFOtVbOqmU/K8YwqCpsx5HnKbp7VsCiSacsKq5ZFMdTCY2hjr4ZFMcz15NOicIoUayNBGmusn6HGKZTjGVUUOYwNx9PQ5tlca21y1enNNLzbayLtMfjsbOy5FkXIT9aTIy/H9eRzXGt/Ip2dD14ftSyK/kG1KJTjF1UUNmMrmG1tko01oWHvPclKDs9ScgraPKfHpkcGsyM+gs9DrqfhFoWfXk8DgylqI5biqbf/rkaFu6KMVlRR5DBGYtnZTbfJtij8ZAA5G3vVgtnJAjEKHwV32WB2OK+Az4dFMZBIZxVEXcSxKFRRKMcvqijGIE7Wk+N68rNpJnKK98B/HUW+4nFe+1U8+a4sP8pxIJGmznY51Ucdi0JdT8rxiyoKm7Ey4Q5yLAo7EOtnIl227iFc3ToKR/GAs7H7K+DLtShqwgGsExXFAAAgAElEQVRfBXf9iRR1tkXhuKD61fWkHMeoosjBjJFwtvOU7aR2+rEoBvM2dudvP5XZwYAQCubHKHy28MhPj/UVo0hnFUW97Xoa0GC2chyjisJm7NgTQ357x6Lws2lWe35Ebidah2goQMLjNcYLZlH5LLhLprIKwpk/rsFs5XjGlaIQkZUisk1EtovIzQU+j4rIQ/bnz4nIzJzPbrGPbxORFTnH7xeRAyLyWp6s/y0i7QVGpB5xxk4w27EowsPeeyG/LXg1usfm1jyA3+6xGYIBIZyjfGrCQeI+miEODKazLqdAQKiLBDVGoRzXlFUUIhIE7gEuAuYBV4nIvLxl1wOHjTGnAHcDd9rnzgOuBOYDK4Ef2PIAHrCPFeJuY8wi+89jRdYoHskqipoqKIq8dNYjZlGkM5429ngyPaxzrCMvYyDlsT9Tf2LIogCoiwQ1RqEc17ixKJYC240xO4wxCWANsCpvzSpgtf36EWC5WNHhVcAaY8ygMWYnsN2WhzHmSaCzCvdQHcaQ7ymeTBMOSra2oBp1FNG8wUB+2oznWxSRUABjhqrAKyGeShPNiU/AUODd6zVaWU9DMusiIY1RKMc1bhTFNGBPzvs2+1jBNcaYFNANjHd5biG+ICKv2O6p1kILROQGEdkoIhs7OjpciCzPGPE8We6dUBAR8VXMBgXqKHwqisGCFkVw2HdVJC+ZGWFROIFtL1PujDF2HcVwi0JjFMrxjBtFUehZO39PLbbGzbn5/DPwPmARsA/4VqFFxph7jTFLjDFLJk6cWEZkecbW4KJ01prwW3yWn84aCAjhoPfRpVYn2uEWQHZutoeNPZ7KDMt4An/DlQZTGdIZk41RQHUUxdpN7Sy7YwOzbv4Vy+7YwNpN7b7kKUo1caMo2oAZOe+nA3uLrRGRENCM5VZyc+4wjDH7jTFpY0wG+CG2q+o9YYyYFPFkJvuU7qdGAUYGs8FSGn7SY6MFYgrg1aJID7s2S14w+1mlOAqhPkdR1EdDvtqMr93Uzi2Pvkp7VwwDtHfFuOXRV1VZKKMGN4riBWC2iMwSkQhWcHpd3pp1wDX268uBDcaKPK4DrrSzomYBs4HnS32ZiEzJefunwGvF1iresFxPQxaFn6aAjpLJ783ktSlg7ljVXHngrd6jkEXhWFNeiu6c7CanMhtsi8JHC4+71m8b0cE3lkxz1/ptnmUqSjUJlVtgjEmJyBeA9UAQuN8Ys0VEbgM2GmPWAfcBPxaR7ViWxJX2uVtE5GFgK5ACPm+MSQOIyM+A84AJItIGfMMYcx/wDyKyCOv5/h3gL6p5w8UYQ4XZxJNDAV6/86MLWRR+ur0m0hkaosP/WWYtAC+uomQhC8WR592iqMu1KCIhBpLeLYq9XbGKjivKe01ZRQFgp6g+lnfs6zmv48AVRc69Hbi9wPGriqy/2s01HQnGTGV2jkXhZ9YDHJneTOPrC1sUXmTGUxla7HoRh6gPi8JpJ54bzK71aVFMbamlvYBSmNpS61mmolQTrcy2GUMGxbCn7GjYX4wiv4UHOK4n7xbFyJiCE3yufDMuZFFks548yIsVsih8xihuWjGHUGD4v8DacJCbVszxLFNRqokqihzGSmV2rt/eT3sMGIop5DZV9B/MLpz15EXmYKEYhY9gdnYMal56bDxpZUN54dLF05g7uTGrLOojQb552RlcuthNJrmiHHlcuZ7GAmMpRjGYTBNtjAKWoujzUSw2mEpn52U7+GniV6wy2/ouD66nZHrY0CKoVjB7eIzC+cwZjVopBvjwqRPZ0znAyRPrVUkoowq1KMYgVvXzkEXht46iUJaSnzqKQpXZzmee5IXys568F9w5A4ryYxTgrzHg/p5BJjVFmdxcw7s9g57lKMqRQBVFDmPE8zSs/5HfltuFFEXUT4yiRGW2l+ssbFF4VxQFLYqoP0WRTGc41D/ICY01TG6q4d1uzXZSRhfqerIZzZXZaze1c9f6beztijG1pZabVszx5ZrIfWr323K7UPDZT4yimOJxPquUQhaFIy/uQV42PTacW5ntjEP15sI72DeIMTCpqYaMMXT0DpJKZ4bN5FCUo4kqihy8tp0+kjhVu05BllO1C3hWFrmbp1VH4acpYGZEVpFX11M6Y0hlTNGNvdLrTKatAHM1LYr+RIpIKDBsE3cyoPKL5txywHY1ndAYJWMMGQMH+xJMbq7xJE9Rqo0+stiM1mD2kajazXXH+CmOgyIWhUfXU6Hivdz3lV5nPK+zrUPQ7kflJZgdS6SHte8A/xbF/p44YFkU4+sjABweSHiSpShHArUoRjnVrtpNpTPDntqjfrvHFogpeHU9Dc23LhajqEymsz5fHlgpsl6D2bmpseA/RrG/17YomqLZeozOflUUyuhBLYocRp/jqXh1rteq3fxBQ06MwqvbrZpZT0N9o6pTR1HMogCn0NBbMLsu36II+7Moumyl0FoXYZxtUaiiUEYTqihsRqnniZtWzMnOZXbwU7XruFuyBXfhoOehQFC8QM6TokiObDAIlqsoFJCKN/as4ilkUYS9pQX3J9LDGgLCUAaU1xhFVyxJfSRIJBSgtU5dT8roQxVFDqMwls2li6fxzcvOIOhU7Ub9Ve06m21ur6fc45XLK2xRDHppCe5cW4GN3YvyKWVR1ISDHlt4pIZlPMFQTUW/x35PXQNJWmwF0VJnFeypRaGMJlRROIzWaDaWsnDqHlbMm+wvNTY5/Cnb71jQgr2Z7BhFpe6sYq4n61jl2VklYxThgMemgOlsTCJXlshQjUWldMcSNNuNC8PBAE01IQ6rolBGEaoojgESqUy2x1BXLOlLlvMUXZMTzAafM64LtPCAygcN5Vs7+TKralF4DGZbMYrhricRoS7sfcqdZVEMtf4YVx+hc8Df71lRqokqimOAWM4G5PWp1WGEReHMoz4CLTcqlVksRmEdqzz4PJgsZVF4VRTpEcFssAYZef3ddMVGKgq1KJTRhCoKm9HreGLYUBy/s5nz3Tt+WnhD8fRY5zNP1xYeuRF7qc0YslAKxSi8uZ4GEulh87Id6iNBXzGK5tpI9n1rXURjFMqowpWiEJGVIrJNRLaLyM0FPo+KyEP258+JyMycz26xj28TkRU5x+8XkQMi8lqerHEi8hsRecv+u9X77VXOaKzOzlUOXlMwHZyn6Jr8GIXHxoC5DQYdIo6VUkXXk5fmhfESFkXUQzDbGEN/IjViAh9AbSTkSYkbY+iOJYZZFC11Ebp9uhgVpZqUVRQiEgTuAS4C5gFXici8vGXXA4eNMacAdwN32ufOwxqLOh9YCfzAlgfwgH0sn5uB3xljZgO/s98fcUZxLDvremquDR8Bi8L7mFFjjJ0eW53eTIXmbzv4sigKWCheFY8xjIhRgGVReHE9DSTSJNNm2BS+lrqwpscqowo3FsVSYLsxZocxJgGsAVblrVkFrLZfPwIsF2uSzSpgjTFm0BizE9huy8MY8yTWfO18cmWtBi6t4H58MwoNimx+/oSGiG+LIv+p3Y/rKZWx+hIV6vUEPmIUVdzYgWzGWC5eYhRO1XR+1hNYMYp+D0rcSU7ItSha66wHAj8V84pSTdwoimnAnpz3bfaxgmuMMSmgGxjv8tx8Jhlj9tmy9gEnuLjG4xrHopjQEPVtUYwouMtOe6t+JXXl6aylsp6CFddmlLIovGQ9OXOxC1kUdeEgMQ8WRZdtOeTGKJrtmopuzXxSRgluFEUhp0z+c3exNW7O9YSI3CAiG0VkY0dHh395ozicPZCjKFIZ46uJ3wiLwkcdhTNKtGjWU8Ube6msp0DFo0tLWxSBituMZy2KgllP3oLZjjJorh1uUYD/VGhFqRZuFEUbMCPn/XRgb7E1IhICmrHcSm7OzWe/iEyxZU0BDhRaZIy51xizxBizZOLEiS5uwx2j0PNEzM56mtBgPWn6SZHNd+/4cT0V29ijfrOeCmQpeY1RhAJScK5DTThIOmNIViBzaGhRoRiFt/TY7gKupxbbuuhSi0IZJbhRFC8As0VklohEsILT6/LWrAOusV9fDmwwVvrQOuBKOytqFjAbeL7M9+XKugb4pYtr9M3oDmZbm9mEBmvOtRdfuEN8RIzCezC72k38BpNpRCAcHPnL8BqjKGSdwFAmVCX3PTQGtYBFEQlWLUbhvPYT0F67qZ1ld2xg1s2/YtkdG1i7qd2zLEUpqyjsmMMXgPXA68DDxpgtInKbiHzCXnYfMF5EtgNfws5UMsZsAR4GtgJPAJ83xqQBRORnwDPAHBFpE5HrbVl3ABeKyFvAhfb794zRmR5rPamOtxXFgI+Adn5RWzae4KH4rFhMwbOiSFkbuxTQ2l7Gq1pzN0Zu6uBteJHjAixUR1EXCZFIZUhVeI2O1dCSE6NwFIXXGIUz7Kq9K4ZhaNiVKgvFK67mURhjHgMeyzv29ZzXceCKIufeDtxe4PhVRdYfApa7ua6xghPMHm+7nvxYFPmbcdRjPAFGVnk7+IlRFHI7gV2ZXaEyc+61EE4Lk8oUhROjKOB6cmZSJNM0VTDCtCuWIBIKDKv1aPHZQbbUsCs/fcKUsYtWZtuMYs8TsaTla2+qsZ40Y35cT8nhdQ/ZGEU1s548xyhG1mRkZXqwKGLJNDUFnv5hSLlVUp3tKOi6AumxjpUxUGFAu3sgSXNteJgVVR8JEg6K52B2tYddKYoqijxGn+NpqG3E0GxmH66nvErqUDBAMCAeYxRVdj0lR/aNcvAyYCmWSI+Y5eHgyfU0WMKisI9VGtDuGkhms5wcRITm2kg2dbZSqj3sSlFUUdiM7mC21Ygu+9Tqy/WULjBq1Ns41KF4R5E6iiq6niLBAMZYRX5uKTSNzsFRFJXct2NRFFI+dR5/N92x5LDUWIfWurDnrKebVswhmPcP2s+wK0VRRZHHKIxlE0taT8bOBuVLUSRHbsZeZj1A8Qly0aC3jrSlXE9e6j1iyQy1BZ7+Yai2ohLX08CgpXgCgZFPFfVRb+NQu4ooCj9tPC5dPI1JzdHszzIg+Bp2pSiqKI4BLNdTKPvU6qU9tkPh3kxBjzGKI5P1VAgvcY9C0+gcoh5cT/2JdMGqbMiJUVT4u+mJDe8c69BSF/FsUcSTaQ70DHL9ubP4uz85jYyBD59avVojZeyhisKmUErmaCGWtJ5k67J+cD/B7JGbcTTs0fVU9TqKEllPHlxFxVqCw1AdRUUWRSJVsM8T5MQoKgxmdw0kClsUtWHPHWR3HRoglTHMmdzInMmNALzxbo8nWYoCqihGYEZhONsJyjobvD9FMXLz9Ox6Sha2KIIBsQPklaazposGs71ZFCUUhYf02P7B4haFY+31VxDMTqatyYUFYxT1Ec+uJye7aXprHSdPrAdg96EBT7IUBVRRHBM4T8aBgFDrsfmcQzyVzm6SDtb0uOrFKMDa2CsfXTpyCFL2Gj3FKNJFXU/ZrKeKLJRUwapsyAlmVxCj6ClQle3QXBsmnsx4cjO22YpiWkstU5prCQWE3Z2qKBTvqKLIY7QGs52NqC4SHFFMVQnxZGZEtXLUwzxqGNq0C23uXuoe4qniFkClFoUxZtjPLR+nFqKSmpT+RLpgnycYCmZXEqNw6iSKBbPBW7+nvV0xwkHhhMYowYAwrbVWFYXiC1UUxwC59QC1kaAv11MsMdK94z1GUbzpXsSD8ik2j9q6xspiFM6QoWIFd05Moa8CC6BkcDwUICCVxSi6SyiKVrs6uytWufup/XCMyc012eysE8fVsUcVheIDVRQ2oziWPczXbrme/NZRVMn1lMxkA9f5eHE9WQqx8BO7Y1G4vU7H6iq2sQcDQk04UJHS7R9MF6zKBisZoi4SqihGkVUUBVxPzsS7w/3eLIppOcV1M8bVsUsVheIDVRSjHGMMA8khi6LOp0URT2ZGFIxFgpV3ZoXS6azRcKCigjtjDAOJFLWR6sQosi3BiwSfwbIqKql7KFXAZ31XZUq80CwKB6ffU7cXi6IrNqwKe0ZrHV0DyYqsJ0XJRRWFzWgdXJRMG9IZk92gaivcjPKJJQtUZvtwPZWqpK7EokikM2SKzKN25IH7GIXzMyrmegIrruBWURhj6I2naKwZuakPk1eJojgCMYpkOsP+njjTcxTF9Fbrdfth7fWkeEMVRR6jLZjtbHhOhXFdJOQ5mJ1MZ0hnTIGsJ2/psbFkZoTSycqscCZ1rER7DKh8foRjdRVzPUFlMyTiyQypjKGxpriFUhcJVmShuIlRHK5QUbzbHSdjhvd1mmYrirbD6n5SvOGqzbhSGWs3tXPX+m3stV0AN62Y47l9woDdADAbzA4HPU+4czbuasUoYolU0RYZdRXGUrIbe9Gsp8ragmRjFCUsioYKLIreuLVhl7IoGqIh+uLufzddA0m7U2zhUa2RUKDiYLZTQ+EoByBrXbRr91jFI2pR2FQrmF3toTH5G6gf11N2hnShpoAerJRSWUqVxlJKDQWC3Lbg7mTGysgDa6SpW4uix1YATSUsiqbaMD1x9xZAsYaAYAXHW+vCdFUYzHaUQa5FMaEhSiQUUNeT4hlVFHn4rcwuNTTGC/kbXl0kWHE/IYeiFkXYm+uplKKorbDew7m2Yq6nbIt1lxt7OcUD1twHtwVyQxZFCUVRE84W0bmhO5akuW5knyeHltqId4siR1EEAsK0llraVFEoHnGlKERkpYhsE5HtInJzgc+jIvKQ/flzIjIz57Nb7OPbRGRFOZki8oCI7BSRzfafRf5u0R3VCmVXe2hMdS2KwoqiJhQklTEkKx0MVEJR1EdCFbnIhu6zSEFbhXUPjkupoUiBHFQWzO61LYpSrqem2lB2nRu6Ywmaa4tfX3NduOIYRXtXnPH1kRG/4+mttdmKbS/oDO6xTVlFISJB4B7gImAecJWIzMtbdj1w2BhzCnA3cKd97jzgSmA+sBL4gYgEXci8yRizyP6z2dcdVojfYHa1h8YMpXkOxSgGU1ZQulKGXE/DNxGn2rjStNuBZKpkN9XKXE+p7HmFCASEukjQtaJwXECN0RJZShUEs3uzrqcSiqImTO9gyvXvppTrCayZFJXOzW7vig2LTzhMa6ml3WMwW2dwK24siqXAdmPMDmNMAlgDrMpbswpYbb9+BFguVjvWVcAaY8ygMWYnsN2W50bmMclNK+aMcJ/4GRozlA3kZD05U+4qtyqcvkb5MYpK3Tq511ZsY6+0piDuIvhciQXgKJSGEq6i+qh7q8eV68ne9N0GtDv7E4yrjxb93KvraWrzSEUxvbWWg30JT72jqu1OVY493CiKacCenPdt9rGCa4wxKaAbGF/i3HIybxeRV0TkbhEp/j+pijjBbL/ZsZcunsbtf3p69n1DNORraMxI15O3kZuQU1uQb1F46HzqXFux9NO6iOXOcpulNFAmPRbsrCK3iiJuFccFCwwZcqiPhkimjasakiHXU6kYhfWZm4B2JmPo7E8wvr5EjKLecj25Hf9qjKH9cBGLIpsiW7n7SWdwK24URaH/afn/coutqfQ4wC3AXOBsYBzw1YIXJXKDiGwUkY0dHR2Flhw1Lpw3Kfv64gVTfE0WG8h70nY25nii8uBzsYCx4z6qxAIo13SvtkKZ5dJjAeqj7usUrOK40tnfTidYN/2ZeuJJRArPy3ZwLAo3cyS6YkkyBsaVUhS1ERKpjOuZGV0DSWLJdEE35/TWOsBbiqzO4FbcKIo2YEbO++nA3mJrRCQENAOdJc4tKtMYs89YDAI/wnJTjcAYc68xZokxZsnEif6nd1WzMjt3o/A6fMYhlue7H5qk5sGiKBLMdjbMSorFnKZ7ResoKrzOcjEK6zpD9Ltsutc3mCoZyIah2IwbK6U3bskrNAbVwYlfuLEoOvsHARjfUFxRtNrV2W7nUrQXyHhycI55SZH90oWnjjgWDQZ0BvcYwo2ieAGYLSKzRCSCFZxel7dmHXCN/fpyYIOx7OV1wJV2VtQsYDbwfCmZIjLF/luAS4HX/Nxgpbg180uRqxy8jrN0yM8GyioKD5lPziabv4F6kZkfZM+nrkKZffFU2Sf2SlxPvYMpGkoEniFnKp2La+yJJ0sGssHKegLoiZW/xkN91uY/vlSMosI2HqUUxaSmGkIB8VSd7Qw/as1pXnjBaSfoDO4xRNnKbGNMSkS+AKwHgsD9xpgtInIbsNEYsw64D/ixiGzHsiSutM/dIiIPA1uBFPB5Y0waoJBM+ysfFJGJWO6pzcDnqne77w3ORlEfCWZnDngllkgTCQWyvnbH9eQlRdaxGPLHedZ7yHoqV6dQW+F19rh4Yrd6KbkPPpcqjoOhQHevCwugeyBZcMBQLpVZFJaiKOl6clqNu7UoDo+synYIBoQpLTWeXE+bdncB8MQXP8ykphou/t5/c6jf2/Q95djEVQsPY8xjwGN5x76e8zoOXFHk3NuB293ItI9f4Oaaqk21gtkwZFHMGFdXUQFWIfKL2rzEExycp/H8lFZnU68kmO24sYpZAJUqn954quwTe0NNBVlP8RSTm2pKrmmt4Im9cyBRclOHoRiFm9/5QXujLeV6yloULv8NtXfFqAkHhj355+K16G7Tni6mNtcwyf55LnvfBO5/eifx5MiW9crxiVZmHwGcjeKk8XWunwaLkZ9ZNBSj8GZR1IZHZgI5m3o1ezPVVphJ1RtPlg0+V+R6chHMbql1hgOV34gP9yeyT/jFaIyGEHGnKDpt11NrmcpscB+jaDs8wIzWOqRIP5rprXWeYhQv7TrM4pNas+8Xn9hKMm3Yuq+nYlnKsYkqiiOA43o4aXw9/Ym0pzGjDrFkaph7pzZb81B5MLs/kcoqhVy8pMc6tQKF5MFQHKSSuodywef6SMjq4uqigrw3nqShRLEdWOmn4M6109mfYFwZ11MgILTWRVy5ZTr7B2msCRUd/ASVxyjaDseyLcULMa2llv298Yr+PR7oidPeFWPxjJbsscUnWq9f3tPlWo5ybKOKIo9qtBnvjlmplE5Q0U/mk+V6GtpAHevCazC7ocCENi9jPB1lWKyyuLmCVFFwmc4adRRa6escTKXpT6QZV196Y2+MhggGpOxGnEpn6ImnaC3jegIYXx/hYN9g2XUH+xJMaChdIlQTDlITDrj+GVqKoq7o59NbazEG9nW7typesuMTi08csigmNdUwuamGzaooxgyqKGyKmete6IlZGTLOE6EfRdGXt4HW+qjM7h8sbFGIiN2bqfKhO01FehVVriiSJfsowVCxW7ngszM+tNzGLiK01IbLunYc11S5GAVYnVqdjKZSvNsTLxtDARhX507x9MSTdMeSpS0KDwOMNu05TCQY4PRpTcOOL5zRrBbFGEIVRT5Vsiiaa8O+xlk69OT57p2nf6/B7GKuoroKitlgyA9fLABdEw4SCbl/Gu6Np0q224Ahf345C8DJKCpV9ezQXBcuG6M4bMsrF6MAmNAYdbWxv9sdZ0pzeUUxpaWWd7vjZde1dVqb/4xxJSyKFuuzSgLam3Z3MW9q04hJhgtntPDOoQHPMThtMnhsoYriCNAdS9JUG6Kl1ts4y1zys4FExB5e5MGiSBSPAzTXhiuyfHriSYJ2o75iNNe6a7udyRi6Y8nsz6sYzhN9Z5kYgPN5qUCxQ2tdpOxm58QcxrmQN74+UtaiyGQM+3viTHajKJprXLXKcOojSlkUk5trCAiuu8gm0xleaevizBy3k8MiO2axyYNVoU0Gjz1UUdhUc2J2TzxlWxT+FUVPLJlNu3Soq6CZXS6l4gCVNqBzrKZSLju3yqcnniSVMYwv47NvdasoBsrXKDi01IbL/n7291hP9JOayrcdm9gYpXcwVbL53sH+QVIZ48qimNZSy97ueNlC0N2djqIoblFEQgEmNdW4Lrrb9m4v8WQmG7zOZcH0FgICm3dXrii0yeCxhyqKPPwOLoIc11MF6ZeFSKUz9CfSIzb3Sp/+HTr7E0Wfspvrym+YufTEUmUL2txe5yGXriLnib6conBcRW6Czy11kbL37bh+3FgAzj2Ucj858ia5iFFMbaklkcqUzaR6u6Oflrpw0RoKh+mtta5jFC/tPgxQUFE0REOcOqnRk0WhTQaPPVRR2FQxlk3XQILm2giNNVZevddgtlMzkB8HaK0LZwO2bkmmM/TGU8UVRYXKp7uApeNVppsqZbAK2gJSvq7AkVfOlQUwoTFCR99gySf2d3vi1EeCZYPtYFkUAB29xRXFvgoUj2N1lNtE3z7QxykTG8omZUxrqXVdnf3SrsOc0Bgt2BIErEyozbsPk6lwNkqxbC9tMjh6UUWRh9/0WGMMXQNJWuvCBAJCU02Ybo8BP6cVSP6G3FwbcV2E5eA8NbcWSRl144LJxZqlUHpjd21R9LlTFMGA0OKiTuFA7yDj6iOEguX/eU9pqiGRypS0Ut7tdhdPgKHNrtRmvPuQ5fo5sUTgOV9eOUWxvaOPU05oKCtvemsd+7rjrmpRNu46zJKZrUWVz+ITW+iJp9hxsL+srFymthT+WV59zkkVyVHeO1RRVJn+RJpUxmTjEy0usmqK4dQq5Lt4Wusqdz05AdtiFkVLXZhYMu1qNgPAgd44E8vFFOoi2QrkUhyyO6mWqysAS5mUk7m3K1b0KTifKfa6fSUyi951GXiGoWDyns7iG/uOg/2Mq4+4yqKaYcccnBhEITr7E3T2J1wpihnjaklnTFmrYn9PnLbDMc46aVzRNWfaLqlNtovKDW2HB3i1vZvlc09gWkstAkxuqiEcFHZ7nMCnHHlc9XoaC1TL8+RsyE58otIn9Vw6i6RlttZXblE4s5eLxyiGUk8nNZXu35POGA72JTihTHB3cnOU/kS6bI3E/p5BAuIu+DyxIcqB3tLponu7YtmOp+XIde2cPq254Jq2wzHOO9VdK/tGu35mT4lNb+fBPmaOL29NgBU7mtAQ5a39fUXXvLm/F4D3uVAU86ZY9/haew8njS/+M9r4jrX5n3XSyIwnh8N0gngAABILSURBVJMnNNBYE+Kl3V1csWRG0XW5rP7jO4gI/+fS04e5mr7yyMv8x0vtfHXl3JLjYYuxdlM7d63fZk34a6nlphVztLttFVGLIg+/oWxHKTgWRXNdxHOMwvFzn9A4fENurg0TT2YqGmvpzD8o1gF1kv0dbnL2O/sTpDOGExpLP2U7wdpyMtsODzC5qaZkOwuH6a2lG9sZY9jbFWNKgXGghXDWvdtT+Bp740k6egc5eWL5TdhhRmtdyWvcebCfWRPcyzt1UgNvHiiuKJzCtwVFFN0wWZMbiAQDvNJeOgj9/M5D1IQDzJ/aVHRNICAsOamVZ3ccKvu9YBV8rnlhDxedPnlEPOIzH5hJLJnmkRfbXMnKRdNtjzyqKGyqVZk9pCisp2OvGUoAHXbmzMQ8ReEmsyafcgHUqVkXTPlAZzEFlo9TeVxsE3ZoKzK+sxAzxtVxoHewqJI8PJCkP5EuWU+Qy/j6CJFQgD1FXDs7Oiz/u1sLBaxmkDs6Cm/sB/sG2d8zyJzJlSiKRt7a31s0aLx5Txcnjqsrm14MEA0FmTulkdfau4uuMcbw+20dLHvfBMJl4jwfOXUiOw/2s+tQ+TjFo5va6Y2nuO7cWSM+O31aM2ed1MrqP75DusLguKbbHnlUUeThd3CR4293+gxZridvweyO3kHqIsER1dReWjG0H7ZaUBdLQR1ywbioArbdKlPKxAEcpbSvjMz2Mj2KcnEUQDEfu+OGmT2p0ZW8QEA4dVIDb7zbW/Dzt+wn+fdVYFHMm9pE2+EY3QVcjq+02U//00emnBZjzuRGBhJpdhVQZsYYNu3uyhbAuWHB9GZe3tNNskhAe8fBfnZ3DnDe3BPKyjpvjrXmD9tKjyM2xvDA0ztZML15WIPBXK4/dxa7Owf4zdZ3y35vLsX+LXiZvaEURhVFlXH+cToujRY78FxpCiFYT/eF+gE5m2olrRjabd9tMctpXH2EaCjgyqJwslzKPWVPbaklEgzwdpGna7DcEXu7Y8ws4S/PZdYEa932Iq6YbfaGP3eyO0VhrW0qqihe3tNFYzTEyRPcWxSnT7VcQFv2jnxq37ynGxGKxkMKcfZMK07w/M6RLp439/fxbk+cc04e71reuadMpG8wxYu7CgehH3tlHwAXuFAUMyfUc/LEen716r6S6za8cYC3O/q59oMzi/4bXDF/MjPG1XLvkzvKfq9De1dsRNt8h9pw0FV2VyG0xchwVFHYVKuOYm9XjJa6cNYKaK4NkzHWaM5K2dHRn90Yc5naUoMIJQOm+bzd0VdyMxYRThpfx9sd5V0I2w/0MaEhWnbQUDgY4H0nFH9aB3jj3R6MsZ7C3XDalCaCASnqOnm5rYvWunBZt1gucyc30tE7WDBIvnHXYRbMaC45eS8fRwlsLLARP/lmBwumNZdtqZ7L+yY2MKEhwtPbRyqK376+H4Dlp5Xf1B3OnT2BcFD4nX1uLpmM4ZGX2vjAyeNdZ45dftZ0nt/ZWdTdlskYvvXrNzlpfB2XLJxaVE4wIHz23JN5aXcXv3/jQNnvPdAb51M/fJZwUIjmuchCASGWTPOVX7xS8UOaxjxG4kpRiMhKEdkmIttF5OYCn0dF5CH78+dEZGbOZ7fYx7eJyIpyMu052s+JyFu2zPKpMFXEbzB7T2eMqTmBVMf3X8wHXox0xrDjYH/BTJZoKMjU5lpXmzpAPJnm7Y7+koFJgDOmtfBKW1dZ99tLuw+zcLq7J+LTJjeyZW930f+sz++0NlO38mrCQWaf0FDwadgYwx+3H+KD75tQUczp/bOsp/H/fvPgsONthwd4fV8P557iLuPJYVx9hEUzWrKbeK68l9u6WH7apIrkiQgXzD2B376+f1jjxnTGsOaF3SydNc5VlbdDQzTEBXNP4JEX20bEev7z1X3sOjTAp95/omt5l585nUgwwL/819sFP//xs7vYuq+HL114atmYx1VLT+TkCfXc+n+3lGx82TWQ4DP3Pc+BnkF++v+cw52XL8im205rqeUfr1jIX3/0VB59qZ2v/fK1ilzK//DEGwVjHnc8/oZrGYU4lq2UsopCRILAPcBFwDzgKhGZl7fseuCwMeYU4G7gTvvceVjzs+cDK4EfiEiwjMw7gbuNMbOBw7bsI041DApjDFv29gxze8ybYm3OhdwQpXilrYtEKlPURbHoxBY2vtPp6j/Aczs7SWdMwVYMw2TOaOZgX6JkAdXuQwPs6Ojn/ScXz6/P5dzZEzjYl2BzW+Esmye2vMtpU5o4oYKN7oK5J/Dczs4Rwfxndhzi3Z4457twmeQyf2oTk5tqWLt5+H/cnz2/G4A/OWNyRfIALjp9Mq+0dWdjEgD3PrkDAf7srOkVy/sfS2YwkEjz42d3ZY89+Nwu9nTGuG7ZzIrlXbdsFocHksPcPAd649z+q63MndzIx8+Y4lrWCU01fOYDJ/HzF9v477eGxyqe3n6Q2x97nfPmTOQTJawJh0gowP+59HR2dQ7wxYc2FazreWt/L5f94I/s6Ojn3s+cxZkntnLp4mk8ffMF7Lzj4zx98wVcungaf7X8FD73kffx4HO7+Zufv1y2P9pr7d18/ZevsbdIlt67PXE+/9OXeG7HoYpjmWs3tXPzL145Zq0UN/bvUmC7MWYHgIisAVYBW3PWrAL+t/36EeD7Yj3SrQLWGGMGgZ0ist2WRyGZIvI6cAHwKXvNalvuP3u6Ow/86pV9nNAYpS4aot4OJDdEQ9mgcjQUQEQwxmCMZYGkM4a+wRQb37E2r7NnDW2iJ46rY2JjlF9u3suF8yZTF7FGkQZECMjIbCtjDP2JNA8+t5tQQDj3lAkFr/Mjp07kV6/s44nX3mXF/MkFXSOZjKEnnuSHT+6gsSbEB99XWJbDx+ZP5hvrtnDPhu38f5edMWweciqdoW8wxZ1PvEEwICVdCLksP20SdZEg//TrN/neVYuz6bnxZIY1L+zm5T1d3LZqvitZDpedOZ1/+a+3+ca6Ldz2ifnUR0PsPNjP13+5hYmNUS5e4H6TAyugfe2ymdzx+Bv88MkdrJg/meff6eSHT+7k4gVTStYbFONT7z+Rf31yB19cs5mvrJzL1n09/Pszu7j6nJNcu3RyWTJzHB89bRLf+vU2Yok0qUyGf/2vHXxo9gRWzK9ckb3/5PF8YuFU/uk3b3KgN860ljp+8uwuumNJfnTt0opcbQB/feGpPLX9INev3sh1y2bxvon1vLS7i4c37uHkCfXc/T8Wubbylp0yga99fB63/edWVn3/aa7+wEmcNK6ewwMJ/rCtg3Uvt9MQDfGTz76fpbOKP7CICF9dOYeacIDv/O4t/vutg3xq6Ym8/+RxnNAYJZEyvNsTY/Oebta/9i7b9vcSCQWoDQcLzntpiIZ46q2D/OqVfcyd3MglC6eyaEYLU1tqaakNkzaGRCrDgd5B9nXF2HN4gF2HBtjdOcAf3z40Ipsrlkxzy6Ov8s4hy8U8c3w9Exuj1Nv7TcYY0hlDKmOIJ9L0DqboH0zRF0/RNzj056OnTTri7U+knGYUkcuBlcaYz9rvrwbeb4z5Qs6a1+w1bfb7t4H3Y23yzxpjfmIfvw943D5thMyc9afYx2cAjxtjTi91jUuWLDEbN26s4LZH8ust73LDj18su06kdJuPKc01rP/rDw/z3//wyR3c/tjrReXlKo1kOpOV/9lzZ/G/Ls433iziyTQf/+5/83ZHPwGx/LsigtgyM4ZhIy///tLT+bSLFgl3rX+De35vuRAiwQAiVp+o3H/jN62Yw+fPP6WsLIefPrebv/2PVwEIB4VUxmTv8UOzJ3D/tWeXdUnkc8/vt49If2yIhrj36rP4YBHlWopEKsNf/Hgjv8/J3lk4vZnV1y11VUFdiOd3dvK5n7yYLZz8+BlT+Nb/WDhMAVdCdyzJF9dsyl7jivmT+McrFrrqQVWIeDLNrf93Cz/f2EYqY5g/tYn/70/PYGEFGVS5HOob5Gu/fI0nXnuXjLFmp1x25nRu+ZO5ZeNZhVi/5V3uePwNduZYuA3REH+6eBpf/OhsV+nADhvf6eTu377JH98+VPD/79kzW7lk4VQ+sXAqf9jWwS2PvjpMWdSGg3zzsjNYMX8y615u58fP7uK19vIzw1vqwpw0ro6X24p7FMrtKeX40bVnV2xFD323vGiMWVJ2nQtFcQWwIm9TX2qM+X9z1myx1+QqiqXAbcAzeYriMSyX1wiZOetzFcVjxpgzClzXDcAN9ts5QDWSpicAB8uuOnqM9uuD0X+Nen3+GBPXF6htGhdsGDdNgqGISacS6b7O9kysp9OrvPDEmWdIMBQBSA90E6yzXMomnUokO9551e/1/v/t3VuI3OUdxvHvg1uJ0ap4oagpGkVUEG1EwSr0wgOUKtELLwTbCgremNYWRSuCFyIqVKSFQkU8BGxUJCouXqjBCL3oQUM8hhR7ocRNo1EEKYLoyuPF+65ONrv/nZ1l5v2nPh9Ydmd29j/Pzun3P72/dwVOsL3kQbhhdj3NAIPj89cA/13kNjOSpoAjgE+X+NuFrv8EOFLSlO3ZRe4LANsPAg8OkX9okrYNU11b6Xs+6H/G5FuZ5Fs5SdtmP9vb64zzDbO9/xpwSj0b6WDKwenpebeZBq6pP18JbHXZVJkGrqpnRa0FTgFeXWyZ9W9eqcugLvO50f+9iIhYqSW3KGzPStoAvAgcBDxie4ekO4FttqeBh4HH6sHqTykf/NTbPUU58D0L3GD7a4CFllnv8lbgSUl3Aa/XZUdERCNLHqP4PpF0fd2l1Ut9zwf9z5h8K5N8K3cgZJwvhSIiIjqlhUdERHRKoaiWalPSkqQfSXpF0k5JOyTd2DrTQuqo+9clPd86y3ySjpS0WdK/6+P4k9aZBkn6XX1u35H0hKThh6qPL9MjkvbWcVJz1x0laUttsbNF0uIzG7XJ94f6HL8l6VlJow0KGVO+gd/dLMmSlj/op4EUCoZuU9LSLHCT7dOB84AbepZvzo3AwiML2/sT8ILt04Cz6FFOSccDvwHOqYNLD6KeENLYRkrrnUG/B16uLXZerpdb2cj++bYAZ9g+E3gXuG3SoQZsZP98c+PDLgF2TTrQqFIoim/blNj+EphrU9ILtvfY3l5//h/lQ65X8zxKWgNcCjzUOst8kg4Hfko9g872l7a7p3ibvCngkDoOaTWLjB+aJNt/o5zFOOhySmsd6vcrJhpqwEL5bL9Ux2AB/JMyFquJRR4/KP3wbmHlPUgnJoWiOB74YODyDD37IJ6j0pl3HfCvtkn280fKi3+0CQDG6yTgY+DRumvsIUnLb+A0JrZ3A/dR1jD3AJ/ZfqltqkUdY3sPlBUYYLTeEZNxLd+1DOoFSeuB3bbfbJ1lOVIoioW6lfWu2ks6DHga+K3tpRvNTIiky4C9tpdultXGFHA28Bfb64DPabvLZB91P//lwFrgOOBQSb9om+rAJul2yi7bTa2zzJG0GrgduKN1luVKoSiGaVPSlKQfUIrEJtvPtM4zzwXAeknvU3bbXSjpr20j7WMGmLE9txW2mVI4+uJi4D3bH9v+CngGOL9xpsV8JOlYgPp96RmGJkzSNcBlwNXu1/n/J1NWBt6s75U1wHZJy2//O2EpFMUwbUqaUenP/DCw0/b9rfPMZ/s222tsn0h57Lba7s0ase0PgQ8knVqvuoh92+S3tgs4T9Lq+lxfRI8Ots8z2K6ndy12JP2M0t1hve3lzRY2Zrbftn207RPre2UGOLu+PnsthYLSpgSYaymyE3hqoKVIH1wA/JKypv5G/fp561AHmF8DmyS9BfwYuLtxnm/VLZ3NwHbgbcr7svnIXUlPAP8ATpU0I+k64F7gEkn/oZy5c2/P8v0Z+CGwpb5PHuhZvgNSRmZHRESnbFFERESnFIqIiOiUQhEREZ1SKCIiolMKRUREdEqhiIiITikUERHRKYUiYgwknVvnRFgl6dA618QZrXNFjCID7iLGRNJdwCrgEEqvqXsaR4oYSQpFxJjUvmGvAV8A59v+unGkiJFk11PE+BwFHEbpPdR8atOIUWWLImJMJE1T2q6vBY61vaFxpIiRTLUOEPH/SNKvgFnbj9c52f8u6ULbW1tni1iubFFERESnHKOIiIhOKRQREdEphSIiIjqlUERERKcUioiI6JRCERERnVIoIiKiUwpFRER0+gbQMP2qFJsr7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    print(y2_.nu.item())\n",
    "    #plt.plot(x[0].detach().numpy(),y_.forward(x)[0,:,0].detach().numpy())\n",
    "    plt.plot(x[0].detach().numpy(),y2_.forward(x)[0,:,0].detach().numpy(), label=\"computed\")\n",
    "    plt.scatter(xint,scaled_pois, label=\"predicted peak\")\n",
    "    plt.ylim(0,np.max(1.5*y2_.forward(x)[0,1:,0].detach().numpy()))\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do a fit\n",
    "\n",
    "First make some dummy data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [],
   "source": [
    "y2_ = pois_conv()\n",
    "with torch.no_grad():\n",
    "    y2_.nu+=1.5\n",
    "    y2_.base_dist.theta+=.1\n",
    "    y2_.base_dist.sigma+=.1\n",
    "data = torch.zeros_like(y2_.forward(x))\n",
    "data.data = y2_.forward(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "t=0\tloss=0.003986922902300565\n",
      "\tparam  1.0\n",
      "\tparam  0.051\n",
      "\tparam  3.0\n",
      "\n",
      "t=500\tloss=0.0004082261412487773\n",
      "\tparam  1.1069657856157538\n",
      "\tparam  0.1769350213826769\n",
      "\tparam  3.5400895045643637\n",
      "\n",
      "t=1000\tloss=0.00019254864368328299\n",
      "\tparam  1.104256439223183\n",
      "\tparam  0.1674518635502764\n",
      "\tparam  3.769002476383715\n",
      "\n",
      "t=1500\tloss=0.0001108094209823346\n",
      "\tparam  1.1031418506638224\n",
      "\tparam  0.16309889590260718\n",
      "\tparam  3.91099116864832\n",
      "\n",
      "t=2000\tloss=7.033956699939845e-05\n",
      "\tparam  1.1024834462038822\n",
      "\tparam  0.16046594784319412\n",
      "\tparam  4.011160320625154\n",
      "\n",
      "t=2500\tloss=4.735130257572868e-05\n",
      "\tparam  1.102034468115289\n",
      "\tparam  0.15866951450029707\n",
      "\tparam  4.086748654745574\n",
      "\n",
      "t=3000\tloss=3.314869426685245e-05\n",
      "\tparam  1.1017038863350115\n",
      "\tparam  0.15735672225680727\n",
      "\tparam  4.146203055164765\n",
      "\n",
      "t=3500\tloss=2.3863101343512058e-05\n",
      "\tparam  1.1014485810067902\n",
      "\tparam  0.15635343805652147\n",
      "\tparam  4.1942969840307684\n",
      "\n",
      "t=4000\tloss=1.7540769419754663e-05\n",
      "\tparam  1.1012449280000847\n",
      "\tparam  0.15556214994448553\n",
      "\tparam  4.233993118894525\n",
      "\n",
      "t=4500\tloss=1.3103148307666452e-05\n",
      "\tparam  1.1010786615022556\n",
      "\tparam  0.15492333853285223\n",
      "\tparam  4.267256955428157\n",
      "\n",
      "t=5000\tloss=9.914245638328102e-06\n",
      "\tparam  1.1009405412434996\n",
      "\tparam  0.1543983108783803\n",
      "\tparam  4.295459111782398\n",
      "\n",
      "t=5500\tloss=7.57951268948578e-06\n",
      "\tparam  1.1008242563348958\n",
      "\tparam  0.1539606685597788\n",
      "\tparam  4.319593106709681\n",
      "\n",
      "t=6000\tloss=5.844156926365047e-06\n",
      "\tparam  1.1007253141406788\n",
      "\tparam  0.15359169879646456\n",
      "\tparam  4.340401711692256\n",
      "\n",
      "t=6500\tloss=4.538209163846984e-06\n",
      "\tparam  1.1006404082358734\n",
      "\tparam  0.15327771812946442\n",
      "\tparam  4.358454400061147\n",
      "\n",
      "t=7000\tloss=3.545217693286112e-06\n",
      "\tparam  1.1005670380541883\n",
      "\tparam  0.15300845972835223\n",
      "\tparam  4.374196991355995\n",
      "\n",
      "t=7500\tloss=2.7836023312572676e-06\n",
      "\tparam  1.1005032692037167\n",
      "\tparam  0.15277605152118925\n",
      "\tparam  4.387984677341635\n",
      "\n",
      "t=8000\tloss=2.1951300970023443e-06\n",
      "\tparam  1.1004475765038195\n",
      "\tparam  0.15257434506746897\n",
      "\tparam  4.4001046931404355\n",
      "\n",
      "t=8500\tloss=1.7375649295138236e-06\n",
      "\tparam  1.1003987377854716\n",
      "\tparam  0.1523984611320626\n",
      "\tparam  4.410792304345434\n",
      "\n",
      "t=9000\tloss=1.3798512568451279e-06\n",
      "\tparam  1.1003557599943536\n",
      "\tparam  0.15224447378635297\n",
      "\tparam  4.4202423467314675\n",
      "\n",
      "t=9500\tloss=1.0988831508716698e-06\n",
      "\tparam  1.100317826502317\n",
      "\tparam  0.15210918571107085\n",
      "\tparam  4.428617727693528\n",
      "\n",
      "t=10000\tloss=8.772923137453809e-07\n",
      "\tparam  1.1002842587278634\n",
      "\tparam  0.15198996510843235\n",
      "\tparam  4.4360558034693085\n"
     ]
    }
   ],
   "source": [
    "y2_ = pois_conv()\n",
    "prefit_prediction = y2_.forward(x)\n",
    "\n",
    "# do the fit\n",
    "learning_rate = 0.5\n",
    "for t in range(10001):\n",
    "    out = y2_.forward(x)\n",
    "    loss = (data-out).pow(2).sum()\n",
    "    if t%500==0: print(f\"\\nt={t}\\tloss={loss}\")\n",
    "    loss.backward()\n",
    "    with torch.no_grad():\n",
    "        for param in y2_.parameters():\n",
    "            if t%500==0: print(f\"\\tparam \", param.item())\n",
    "            param -= learning_rate*param.grad\n",
    "            param.grad.zero_()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEKCAYAAAA8QgPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXmUXGWd8P/53trXru5OZ09IgLAkLBEQHQFFwSPMj4FRcQivjviOIwPqvOf8Zjni8qqvc/DI6Iy/13FhdHDQGYdFHTE6LKMiioJAkC0sgSxAOiTppNfq2qvu8/vj3qqu5dZe1elOns85OVTf5XufCunne7+7KKXQaDQajaZdjCO9AI1Go9EsTrQC0Wg0Gk1HaAWi0Wg0mo7QCkSj0Wg0HaEViEaj0Wg6QisQjUaj0XSEViAajUaj6QitQDQajUbTEVqBaDQajaYj3Ed6Af1kyZIlat26dW3fl82kyEy+iPItJTq4oqs1ZAtZMoUMYW+YvJknnU8T8oQwpDPdncglcIkLv9sPQM7MdSUzmUsiIgTcgYpnGGJUHGsVU5kkcgn8bj8ew+O4Zo1Gs7B5/PHHDyulRppdd1QrkHXr1rFt27a273t1z7PsvvMi0id8mEuu/GRXa3h15lV2T+/mglUXcCh1iBcmXuDc5ecS9AQ7kvfI/keIeCNsHN4IwIHEga5kPnbgMQLuAKctOa3hsVZJ5pI8euBRNg5vZGlwKQC/2/87BrwDnDp8atvyNBrN/CMir7RynXZhLTIU/e9dJgid9khzWp8g87JujUYzv2gFMo8I0nM5vZDZq3UBJcVTLVMrEI3m6EMrkD5T3DhFerdJ9xuRzi0Gp+/bSwWl0WgWDkd1DORopNft95VS9Ht/FxG0AaJZaORyOUZHR0mn00d6KUcMv9/P6tWr8Xg8Hd2vFcgipOLt3v7cSxdRNzGLevdpF5ZmoTE6OkokEmHdunWLykPQK5RSjI+PMzo6yvr16zuSoV1YDej1234/NvtuUai+x0C6CcprNP0inU4zPDx8TCoPsPaj4eHhriwwrUCc6OFGXy+o3LG8qg2/+LmXG7RI9xt+9S/lQlKaGk2RY1V5FOn2+2sFcowzL2nBx/gvqUZztKIVyDzSL2uhaxlV1lE31lJdF5a2QDSapnz2s5/lS1/6Ut3zd911F88999w8rqgxWoHMEz17C6+zD3cc9O51VhcOCqQHLjGNRqMVyDHHYnjzdlJuug5Eo5kfbrzxRk4++WQuvvhiduzYAcC3vvUtXv/613PmmWfy7ne/m2QyyUMPPcTWrVv527/9WzZv3syuXbscr5tPdBrvPNLPWECv0257zWJQpJpjl//zk2d57rWZnsrcuDLKZ/5oU8NrHn/8cW6//XaeeOIJ8vk8Z511FmeffTbvete7+NCHPgTApz71KW655Rb+8i//kssvv5zLLruMK6+8EoBYLOZ43XyhFcgRoJvNtB8bcU0MRATTNDuSpV1VGk3rPPjgg7zzne8kGLQaoV5++eUAbN++nU996lNMTU0xOzvLO97xDsf7W72uX2gF0mf6suEv4KwmRxdWFwpJo5kPmlkK/cTp9/kDH/gAd911F2eeeSa33norDzzwgOO9rV7XL3QMxIFeb9A9LdSrUkhdZ3Y53NaLrCmdhaXRNOfNb34zP/rRj0ilUsTjcX7yk58AEI/HWbFiBblcju9973ul6yORCPF4vPRzvevmi5YUiIhcIiI7RGSniNzgcN4nInfY5x8RkXVl5z5uH98hIu8oO/5tERkTke1Vsu4QkSftPy+LyJP28XUikio7d3OnX/pI08vNebGgXVsaTS1nnXUWV111FZs3b+bd7343F1xwAQB/93d/xxve8Abe/va3c8opp5Su37JlC1/84hd53etex65du+peN180dWGJiAv4GvB2YBR4TES2KqXKc8k+CEwqpU4UkS3ATcBVIrIR2AJsAlYCPxeRk5RSBeBW4KvAd8ufp5S6quzZ/wBMl53epZTa3P7XPHrp1lpSqBoZXc0D6WMdSMEsoFC4De151Rw9fPKTn+STn6wdXHf99dfXHDvvvPMq0nivv/56x+vmi1YskHOBnUqp3UqpLHA7cEXVNVcA37E//wC4SKxd6QrgdqVURim1B9hpy0Mp9Wtgot5D7fv/BLitje+z8FALO2ZRj56m8fbo+z8x9gTbDrY/YVKj0fSHVhTIKmBv2c+j9jHHa5RSeSyrYbjFe+txAXBQKfVS2bH1IvKEiPxKRC5oUc6CoReV6Er1uPmhQzPFXle392q9s7lZ0vm0dodpNAuEVnwBTr/91b/B9a5p5d56XE2l9bEfWKuUGheRs4G7RGSTUqoieVtErgWuBVi7dm2Lj1q89KM9Sjc4rkO6X1/BLJQ+Z80sPpevK3kajaZ7WrFARoE1ZT+vBl6rd42IuIEBLPdUK/fWYMt4F3BH8ZjtBhu3Pz8O7AJOqr5XKfVNpdQ5SqlzRkZGmn65hiyQTbmfOFk0XfXC6tM8kIKaUyDlykSj0Rw5WlEgjwEbRGS9iHixguJbq67ZClxjf74SuF9Zr5xbgS12ltZ6YAPwaAvPvBh4QSk1WjwgIiN2QB8ROd6WtbsFWW0j0rvs5vKNsxeuoeqgd8kCWWBpsr1uZZI386XPOTPXtTyNRtM9TV1YSqm8iHwUuA9wAd9WSj0rIp8DtimltgK3AP8mIjuxLI8t9r3PisidwHNAHviInYGFiNwGXAgsEZFR4DNKqVvsx26hNnj+ZuBzIpIHCsB1Sqm6QXhNazgqnh64nCrFdZ+FlVdzCqRcmWg0miNHS/mQSqm7gburjn267HMaeE+de28EbnQ4fnWD533A4dgPgR+2st6jnepOt13Lc0rj7bK7b6+78Za7rbQC0RwtfOUrX+Eb3/gGBw4c4GMf+xg33HADd911FyeddBIbN2480strik6od6KHabfzFdzu5g2/H5Xy1S6snsZAlI6BaI4Ovv71r3PPPfdUzCS/6667uOyyyxaFAtGtTOaBnleOL76ykq4ptzq0AtEcDVx33XXs3r2byy+/nC9/+ct89KMfdWzZvpDRFsg80qs6kF6ilKpRSF25sBwGStknuqJcgZhKN2bU9Jh7boADz/RW5vLT4dIv1D198803c++99/LLX/6Sn/70pwC86U1vqmnZvpDRFkgjFlEab18USxf31cRAutQg5evRCkSjWRhoC6TPOPWa6lreAq5ELymKMhG9iIGYzCkNrUA0PaeBpaCpj7ZAHOhX76p+1GwU19rLIHqv26/3QuEVLRC34dYxEM1RTXXL9oWMViCLCCf3UC/k9ZSSAVK5xm6fZWIiCIYYC6Zti0bTD6pbti9ktAurAQututuJXrQd6bWLDfoTAxERXOLSFojmqOHll18GrMmCH/jAB4Dalu0LGW2B9BmnVia9bJXeNzoecOhcB9L1cuyeXYYYOgai0SwQtAKZB/o5QbAXsp2C6D21vnrQGsXExBADQwxtgWg0CwStQI4SOtmg+xFL6HWcplxu0YWlYyAazcJAK5A+08vNrudB9DpFfz0pJOxxKxNTzQXRtQWi0SwMtAJxoNdvz73OSKqQ3Yt4SA+/rmNdSS8q8FElF5aOgWg0CwOtQBYRdduE0Flgvt49Qhfdcx1mwPeiVkUpS4H0orOvRqPpDVqBaPrSjbfXFK0OA6OiKl2j0VhkMhkuvvhiNm/ezB133MGf//mfl9KBP//5z/flmboOZD6w9+depfFWtwmBDt1DdW7pdRZWxRo71FXahaU5VikUCrhcrqbXPfHEE+RyOZ588kkArrrqqtK5z3/+83ziE5/o+dq0BXIM05dCwgYz1rtRSqay0niLz9BoFjsvv/wyp5xyCtdccw1nnHEGV155JclkEoB169bxuc99jvPPP5/vf//77Nq1i0suuYSzzz6bCy64gBdeeKFC1tjYGO973/t48sknS23gL7zwQrZt28YNN9xAKpVi8+bNvPe97+3pd9AWiKaGbqvbaxSSzJ3rSm6xlcki6BCgWVzsnNzJbG62pzLDnjAnDp7Y8JodO3Zwyy23cN555/Fnf/ZnfP3rX+dv/uZvAPD7/fzmN78B4KKLLuLmm29mw4YNPPLII3z4wx/m/vvvL8lZunQp//Iv/8KXvvSlUmv4Il/4whf46le/WrJMeolWII3o8Ztu17GGOn2mrFOdB9H73d23J3KVwjC0C0tzdLFmzRrOO+88AN73vvfxla98paRAii6o2dlZHnroId7znrmp4ZlMZv4X60BLCkRELgH+L+AC/kUp9YWq8z7gu8DZwDhwlVLqZfvcx4EPAgXgfyml7rOPfxu4DBhTSp1WJuuzwIeAQ/ahT9gz2evKWsgshrflmg2/i8pxp/tkzgTpGFOZePDoLCxNX2hmKfSLehmLAKFQCADTNInFYn2xILqlaQxERFzA14BLgY3A1SJSPaz3g8CkUupE4MvATfa9G4EtwCbgEuDrtjyAW+1jTnxZKbXZ/lNUHo1k9Rajv3UgneJYpNfN7I4+bcT15ot068IyxMDAcmFpJaI5Gnj11Vd5+OGHAbjttts4//zza66JRqOsX7+e73//+4D1e/vUU0+19RyPx0Mul+t+wVW0EkQ/F9iplNqtlMoCtwNXVF1zBfAd+/MPgIvE2jWuAG5XSmWUUnuAnbY8lFK/BibaWGtdWYuFfvbE6oaat6BuK8cdRuRC90F0hFIgXbuxNEcDp556Kt/5znc444wzmJiY4Prrr3e87nvf+x633HILZ555Jps2beLHP/5xW8+59tprOeOMM45IEH0VsLfs51HgDfWuUUrlRWQaGLaP/67q3lUtPPOjIvJ+YBvw10qpyS5kHVGc3pQ7dg85tUrvQTv3XtKoEr1buQZGSdmZmLjojwGq0cwXhmFw88031xwvtnkvsn79eu69996Gsi688EIuvPDC0s8PPPBA6fNNN93ETTfd1M1SHWnFAnH67a/eeepd08q91XwDOAHYDOwH/qGNdSAi14rINhHZdujQIYdb5p/5aL/ey2aK3dSWFJsetvO8Viim8RYtkEUQWtJojnpaUSCjwJqyn1cDr9W7RkTcwACWe6qVeytQSh1UShWUUibwLebcVC3JUkp9Uyl1jlLqnJGRkSZfbXHS715d3dBoxnq3rUzAqkQHdDW6ZtGzbt06tm/ffqSX0RWtKJDHgA0isl5EvFiB7K1V12wFrrE/Xwncr6zf+K3AFhHxich6YAPwaKOHiciKsh/fCRT/htuW1Sm93J57OlDKKcOpD9ZNLzZ8J3rRTLHkwtIxEE0PONaTMbr9/k1jIHZM46PAfVhpvN9WSj0rIp8DtimltgK3AP8mIjuxLI8t9r3PisidwHNAHviIUlYvbhG5DbgQWCIio8BnlFK3AH8vIpuxnBQvA3/RTFb/6M0/rvkIni+UdOF+xUCKLqxedPbVaMAq1BsfH2d4eHh+pnwuMJRSjI+P4/f7O5bRUh2InUp7d9WxT5d9TgPvqb7PPncjcKPD8avrXP+nDdbhKKv3LMx/TL1uPVJPXlcbfoNuvN1QjK3oLCxNr1i9ejWjo6MslFjpkcDv97N69eqO79eV6IucXqTI1lMYnTQ/bLSOXrQyKc/C0mi6wePxsH79+iO9jEWNbqY4j3Trfum126ZpFlYn7VEaVKJ3uv5SK/eyLCztwtJojjxagfSZfscmehET6Gk33j5kYZXPAillYWkXlkZzxNEKZB5YqGm3dTf0LrrnOnbjbfa8FmQW19WvDDGNRtM+WoE0YFFsUl1u9v2g1wqzVANS5sLSFohGc+TRCqQRvdhfHWQsNMVUN+22kyGHDpXoPYuBoNN4NZqFhFYg80C/0247ljfP3Xg7pfx7awtEo1k4aAXiwGIqKurFrI16dRu9cov1KvusPI13oVlxGs2xiFYgfcaplUnHsso20sVEt+t1SuPVFohGc+TRCkRTNwbSaR1Iv7KwRESn8Wo0CwitQBYhFRMJe1z0Vy6zE/pZByJIV1lnGo2mt2gFMo90HQvodXdchwFVFec7nQfSQ3nl9xVH2nYjS6PR9A6tQBY5vQgq97wSvYdBeZjre1XRykRbIBrNEUcrkIZ0v0n18k25XxZINb3OdOo2iF7+d6jngWg0CwetQBzodRpvcQNdTOnB3dDzSnTmXFjF/2oXlkZz5NEK5AjQ7UTC8g262/nl3Zyve0+V/uhlN97qYxqN5sihFUifWci++pYGSr30M/j8anj81pZl9joLq1pxGmIs6L9XjeZYQSuQeaDnLrHyNN5+V2Y/+I+QjcOvvggtWBCNRtp2m31W/K4GhrZANJoFgFYg80gvpgf2g7oxi1wKc++jTMkAzIzCoR1NZfUjNlHjwpKF93eo0RyLtKRAROQSEdkhIjtF5AaH8z4RucM+/4iIrCs793H7+A4ReUfZ8W+LyJiIbK+S9UUReUFEnhaRH4lIzD6+TkRSIvKk/efmTr/0fNKPja7f80VKVs2BZzBUnm9m7f9tB55uKqsf80DKu/GCDqJrNAuFpgpERFzA14BLgY3A1SKyseqyDwKTSqkTgS8DN9n3bgS2AJuAS4Cv2/IAbrWPVfMz4DSl1BnAi8DHy87tUkpttv9c19pX7IaFkcZapFHleD8q0XMHngPgJ+abyCg38VeebEluv7rxFsUaGHomukazAGjFAjkX2KmU2q2UygK3A1dUXXMF8B378w+Ai8TaNa4AbldKZZRSe4CdtjyUUr8GJqofppT6b6VU3v7xd8DqNr/TgqVXsyx6boHU2eAPj75IRnl454VvZLdaSWrfdsfrKnD4ak7fW5kmv//Pf+TZ+25pLrJYiW7/cxURbYFoNAuAVhTIKmBv2c+j9jHHa+zNfxoYbvHeRvwZcE/Zz+tF5AkR+ZWIXNCGnLYwSkbSwmK+CwkL4y+zV43wR5tX86paimtmr+P11TLrVaKX8/QvbuOsp/8Pmx7+K/Y+/UBDmdUxEEErEI1mIdCKAnF6Pa3+7a13TSv3Oj9U5JNAHviefWg/sFYp9Trgr4D/EJGow33Xisg2Edl26NChVh7Vd/pdQCgiHU8PhPoWjTu+l70s44SRMDO+FYRSrzXNxHLKwio/V8R88j+YUUEyysOh39zaVCaUZWGJdmFpNAuBVhTIKLCm7OfVwGv1rhERNzCA5Z5q5d4aROQa4DLgvcre5Ww32Lj9+XFgF3BS9b1KqW8qpc5RSp0zMjLSwtfrLxVvyl12ku31RMLGD1MMpF9jxr8KwxAy4dX4VRqSNV7HGpq1hy/kc5yU2MbzQxfztP8clh/+XZOlqIoiQu3C0mgWBq0okMeADSKyXkS8WEHxrVXXbAWusT9fCdxvb/xbgS12ltZ6YAPwaKOHicglwMeAy5VSybLjI8UAvIgcb8va3cL6jwl6NT2wRHYWLymyUVv/x9Za/516pbFMh429Ogay54UnCJHGtf48Ustfz0pzP4mJ/XVlmsqsqb7XabwazZGnqQKxYxofBe4DngfuVEo9KyKfE5HL7ctuAYZFZCeWe+kG+95ngTuB54B7gY8opQoAInIb8DBwsoiMisgHbVlfBSLAz6rSdd8MPC0iT2EF6q9TSjV/Hda0jSCo1BQARmQFAN5hS4Gkx5soEKc03irLa3LX4wAsPen1BNafC8Brzz/cskxDdCGhRrMQcLdykVLqbuDuqmOfLvucBt5T594bgRsdjl9d5/oT6xz/IfDDVta70Oh3Gm+zc81wcollZydQgDe2HICBpZYlMn1oH/5m8pp8X3XgaTLKw4oTTscIL4MHIL73WeBdjtebyixlYIFWIBrNQkFXojdgofrZG/au6kSeQ91GcmYCEIJDKwFYsnQVBSUkJxqHsBo1Uyx6nXxTuxl1rcLj8bJy5SoOqwHUoRfqy6yyQBbbTHiN5mhFKxAHej1gqSS3aiPtRlYvaKQgs7OWdzAybGVdL4sFmSCCGT/YWKbDGquD6AOZfcQDVnmPYQj7PWuJxuuHs0xl1gTRtQWi0Rx5tAJZhDhZDL2cSCgI+eQkaeViZGQpAEvCPg6rGEaieWp0o268hYLJ8sIBstG1pfPTwbUM5/bVlVedGqxdWBrNwkArkD6jVP3eUEeauoWECKSmGWeAZQNWxMPvcTFhDOJJ1VcgzepKAMYP7MUvORhcVzpWiK5mSE1TyCQd76n+O9RZWBrNwkArkHmkV914exaUp/6Gb2SmmZQBgt65PItZ9xCB7HhzeXViNEopJva9CIBv5ITSeffQOgDG9+10lOvkwlqo8SmN5lhCK5BFRK+bKZbfX/GzCJ7cDHH3YMXxtG8JkfxE3Wr0ehZIuQsrddBSEpEVG0rnw8uOB+orkBoXlp4HotEsCLQC6TONWnscaRq9xQcKcRKeSgWSDy7BSw7S0w3lOn3fopLLTbwKwMiq9aVzg6uszO3U2J6666yxQLQLS6M54mgF0pDeblK9mh7o5CLqaH55VZv0OXkQMJNkPQOVx4N2a5jE4cbyHCitOX6AKcJEInNtzEZWHEdWuVB1qtxNZdZMYdQuLI3myKMVSJ9RqncWyLy9dWcTGCgK3spelZ7IEgAycedAej2FVHYB7uRBJo3hisMBn4eDsgRP3DkTq7qViYE1E10rEY3myKIViAMLNWuqSN8nEqZnrA/+SMVx/4ClQOITY86CSvqjfjfeUOYQCe9wzbkp1zC+tLNcRa0Lq3hco9EcObQC6TONRry2LaveG3en4uuIU8UYh7/SheWPWi6s9ExjC8Tp+xbdbAOFCdL+pTXnE94Rwtk6rrGqNN5iWxMdSNdojixagcwjvZpI2CvqpfGm4pMAuAKVLqxQzNr4s/HGMRDHILoIpllgWE2SCy6rOZ8NLGXQdM7wMql0YXXbFl+j0fQGrUAWEY02zF5WoqdnLAXiCcUqjkcGhsgrg/xs4ybI9VxYyZlxPFJAhZfXnFOR5QRJk0/VZnhVZ2EVPy8URazRHKtoBTIP9DuNVxB4+bdw3yehkGv5vnpKJ5uwFIgvVOnCioV8TBFGJZ2LCRtt6IKUZn5IdEXNefeAdWzy4Kij3OogOmgXlkZzpNEKZB7pVxovZgH188/Cw1+FZ3/UvrwqBZdPWLNAvOFKCyQW9DClwhipSUc5zWIgyWmrEaN3cGXNef+g1bRxeuzVmnMmtZXooBWIRnOk0QqkET1wkTi1N+81MrEH8inrh12/bPm+ehZDITmFCYQilQrE4zKYkQjubB0F0ujvSyBnB9+Dw6trTodHrGPJidpU3moLZKEWZmo0xxpagTghC/uvpSbtdnIPCsgOnwqjj7Usp57FYKanSeMjHPDW3JNwR/Fm269EByjYBYjhoVoLJDpiDawqTNeOtq1O4y1+1haIRnNkWdg75VGAUyuTTl1Ydd/wJ/aQET//OrYBNbkHCvmO5JdIz5BUPiI+V82pjDtGIN9YgThRHJMbVwGGY5Ga80NDwySVD+K1CsSpEh10FpZGc6TRCmSe6YX7pdpiyE2O8kphCTvN5YiZh+m9Lcmp16fLnZkhhY+Ap1aBZL0xQoWZuvKc1lc8JulpJolUdPgt4vO4OSwx3InKYsKi0qyIgaBjIBrNQkArkEVEvTfuwuwYY8R42bTTYyfqT/drBVd2hrQEHGM3eV8MH1nI1s7uaDQPRBBcuRnixkDNuSIzxiDeTGWGV1FJVA+UKn+eRqM5MrSkQETkEhHZISI7ReQGh/M+EbnDPv+IiKwrO/dx+/gOEXlH2fFvi8iYiGyvkjUkIj8TkZfs/w7ax0VEvmLLelpEzur0S88nfe/GqxTe1DiHiXHG6WdYx+o0JXRYnKO14M3FyboCzgorMGT9N1VbC9LIAgHw5mdIuGOO5wASniECuUq5JmaNzFIMBGcL5PFXJvnaL3cynWo9pVmj0bRPUwUiIi7ga8ClwEbgahHZWHXZB4FJpdSJwJeBm+x7NwJbgE3AJcDXbXkAt9rHqrkB+IVSagPwC/tn7OdvsP9cC3yjta+4sBCRrmeiVyik1CRulSXrH2bFyrWYSkhPvNbVGgOFWXKugOMbvoQsBZKfra0FadiNF8FfiJPx1lcgWd8wkcJUtdDS/TXPc1jf3okkV3/rd3zxvh389Z1P1X2WRqPpnlYskHOBnUqp3UqpLHA7cEXVNVcA37E//wC4SKxXxiuA25VSGaXUHmCnLQ+l1K8Bp5LmclnfAf647Ph3lcXvgJiI1Fak9ZTepPH2tTnjtFV4p0JLOW5pjHEizI7Xny9esbY63y9kzlJwBx3Pe2wFkpiubWfS0IUlQtCcJe8bqrueQnCYATVTkQRQtDIcK9Ed1vcvv97FjcY/81zoOtw7fsILB5zjNRqNpntaUSCrgPKo7Kh9zPEapVQemAaGW7y3mmVKqf22rP1AsfNeS7JE5FoR2SYi2w4dqj+/+0jSbRZWxQYd348ARmQJ64aDjKlB8lOtWSCOreYLOfwqhekOOt7ji1odeZMOCqQh+Qxe8hQCtZ14S+sJLsVAVfTaKsZAWknjNU1F6pkf8x7jlwQky997vsl/P7GrvXVqNJqWaUWBOL0+V++A9a5p5d5WaUmWUuqbSqlzlFLnjIyMdPik3lkM/U41zUwfQBT4o0tYPuBnTMVwJQ90LjA9gwDKE3T8P1Vs6Z6p01ARnGMgucQUCjBC9RWIO2q9K0yPzylAJ6VZrynl9tem+cPsf5MMrET+9C6ikoRnflD3eRqNpjtaUSCjwJqyn1cD1a+4pWtExA0MYLmnWrm3moNF15T932JeZyeyFgTVm1+n2UNOA5viE5ayiMSGCfvcTBhD+NOtWQdOreZVahIBxBdyVH6hAUspOzVUbOTCSsYnUALuSH2l7otZWWTxw2UKxCEwX68O5Jmdr3C+8QzmpnfB2jcy7V/FxvhvdTBdo+kTrSiQx4ANIrJeRLxYQfGtVddsBa6xP18J3K+s3WQrsMXO0lqPFQB/tMnzymVdA/y47Pj77WysNwLTRVfXsUxqaoyccjE0EEVESHqXEMpNgNm8RsJJQaTiE5ZLzBt2PD8QiZBRHgqJBllYDgokY3f49TZQIMHB5fZ3mrOgGrmwqhXx7Eu/wS0m4dMuBRESa9/GecazPPPKwbrP1Gg0ndNUgdgxjY8C9wHPA3cqpZ5g5d4uAAAgAElEQVQVkc+JyOX2ZbcAwyKyE/gr7MwppdSzwJ3Ac8C9wEeUUgUAEbkNeBg4WURGReSDtqwvAG8XkZeAt9s/A9wN7MYKxH8L+HBX37whPQx6V+3BItJTt5Y5e4gZQgwEfQAUAkMYmJCeanKncwwkNWNlV7n8Ycd7YkEfU4TAoaFiozTerO3CCg7VzgIpMrDEyonITc8VEzrVgZQKCavSeEMHHyOPG1ZaGd6xTRcTkCyvPd/snUWj0XRCbUmwA0qpu7E28PJjny77nAbeU+feG4EbHY5fXef6ceAih+MK+Egr6z1qcUhpleQ4MwRZGbCyo43gMMxgbfDB+hlPpftrhklNMAR4AmFHV1vE7+aACmO0oKDKKSQtBRIdrK9AhoaXkVMuzNna0bZOLqzyIPp0MsfJ2Wc5HDuV5V4rASC4/g3Ws/c+BryrrfVqNJrm6Er0PlNtbQi9tUDc6QlmVajUdsQTsYLc1JnZUUOVsZCdtWIg3kDEcZ2GIcwaEdzZWgVSshYcLBCVmiavDGKD9V1YIb+HCaJIci57ruTCKvunWvxcruBeOjDF6bKH7PKz5wRGVzDpHmHp9DN1n6nRaDpHK5BG9KhVRs9mojvEGPzZCdKuSEkReCJWllNmpnkg3alKPm/HNnzBaF1Fl3JF8Obq11c41oGkZ4gTZCBY2+G3nGkjhic1p/yclJJTEP3AKy/glxzhNWdUyBuPbmJdbiepbKHhczUaTftoBeJAL+v+HC2QLhVT+WYayk+R8cx1t/VFrTf85FTzwLFTkWMhOUVGuQkGAnXvS3sG8OdrFUijLCxXdppZCWEYjf9yZz1D+LNzAfri31+zOpDkqGVlxNZtrhS4dCPr5AC7XmuzbkWj0TRFK5B5oFe9sGosgmwSn0qT9wyUNu9gzKqlSHdogZCeYoYgIZ+nbsVOzjPg2JG3URDdl4uTcDkH5ise7x0ikp8L0NdTSoYYFX8frsMvWMeXnVJxXXjt6bhEcWD3002frdFo2kMrkPmmC11Ss5kWGxr6oqVrogND5JSLXLyzKnxJzxBXIYJed/3uv74YfjKQSzdeXxm+/CwZd6jp8/P+YQbUVMl96JTGW/08gIH4Sxz2rARv5TOWHG9ZJEULRaPR9A6tQBrQiwiIY7uQLim94aetwU4u/1zAezDkY5IIpkOdRs3aHCwQd3aaWSOMSP126SpgN0SsysRqlBwQMmfJeWoHSdXIDo3gJ0shM1shs9qqMcQoKZdkNs/a/CvMRE+skece2UAON57xHU2frdFo2kMrkHmmmyys6vsys9YG7glGS5v9UMjLpAojqRaysBzauXtycVJGuKHSEzs9uJ6SqnFhmQXCKkHBG3W8vhwjbMVwpg9bDSHrWTWClBTI3sOzrJMDFIY21Ap0eTjoWcNgQvfE0mh6jVYgfaaXKbvVFsOs3dDQG5jbmAcCHiaJ4ErXFvq1sjZvPk7aHWlY8OiyFUh1Q0Wn6YEAheQkBgrlb65AvANWNfrMYavJgFM3XqhUUmP7duGVAr6lDgoEmA0fx0huH6apB1BpNL1EK5CG9D6NV6S7LKxyWem43R4kNJdy6zKEhBHBk20+t9zJvRYsxMl5og0tJW/YShVOzlTGWepdHx+3Owa3oEACxXYmk1Y7k3q1JeUurNn9LwEwsMpZgeRjx7OaMcama6coajSaztEKZDFRtT/nEpYC8Qcrx8Rm3RG8+XgL4qoEKkVQJcjbrqZ6is4/YCmQ9Eylm6yUclv1zyo+cRABXMHmCiQybCmQzPSBolCgjgvLtk7yhy33VHSFswLxLd2AT/Lsf/Wlps/XaDStoxWIA70eAFXTjbdHMZB80oqB+EOVCiTvjeIvzLYkr+K75pJ4yKN80YYxkKDd0j036zy/vPrW1JSlQLyhwaZrGhyxRrzkZ6x2JvVcWIYYJeXinn6FHG5kYLWjzIFVJwEwNfpC0+drNJrW0QpkEVFtEZjJKZLKRzjgrzhX8EYJqFTFZL96VCgKO6uLQKxhunEkOkReGbUdeetYC9mZQ4gCb7j+ONsiA+EQ0yqEJCz3WL003vIgejixl3HPcjBcODF8nDWBOXtoZ9PnazSa1tEKpM8opWo2424C6+Wbs0pPM0OQsL+yZkMV60Iyjce5VleiZ+KWRWEEYg0tpVjIyxRhVFVH3nrt3POzhxAgGG2uQESESYnhTh0uyXSyhgwxMDFRSjGce414YE3NNUVc0RWk8OGe3N30+RqNpnW0AplnuqkJqXY5GZlpZlSQkLeqqbLfuU6j2XoSM5ZF4bZdTfViIAMBj2UlpJzrQGqGVCXGSSkvoYC/6XoAZtyD+DOHS2twcikWkxEOxzOs4SD5gePqCxThsGcV4eSrLT1fo9G0hlYgfaamF5ZIr5K7cGVnSEgYl1HZ1kMCVkzETDbOxKpeW3EWiCc01DAO5HEZxCWCO1OlQOrUbLhS48QJ43a1pjyT3iWE8pYyM5XpWIVetJD27x8lIilcw8c3lBkPrWUku6/rPmQajWYOrUDmgV72wiqX5cnFSTn0l3IHLQskOdO4mLDGhWWPqfVHBpsG+xOuCJ5cpYKqZ4F4MpMkjeZ9sIpk/UuIFSz3mKlM5+6+YsVApva9CEB4xUkNZeZj61nNQcbjqdKxXCHH4ZRusqjRdIpWII3o09tqx1lYVXUbvnyctCtcU1viCVuFful443Ym1evIzVqbdjC6pKkCSbujNR1568UrfLkp0u5oy2//ZnCEMElUNgE4Z8UZGCilSB+0AuNDa05uKNMzcgJeKbD/1blA+q7pXWw/vJ2pNodjaTQaC61AHOjZ/I4GzQU7pkyUv2D1l6re7P0RS4EULYrG4uYEFpKWAgkP2JMMG+z3We8AweqOvA6tUQDChSlybucBVY5rCltTC2cn9td1YZW68U7uwUTwjzR2YUVXWgpmet9cT6wp2wU33ULRpUajqUUrkCNAT/zwShFSs6Wiv3ICUUsBFAsNG62jXIGYqWlmlZ+BUKCpEi14Y4RUEgq5OXl1FETEnKHgGyg9sxmemN3O5NAoJo1dWN74XiaNIfA0DtAPr7UUSGZsridWtpAFIJnTFeoaTSe0pEBE5BIR2SEiO0XkBofzPhG5wz7/iIisKzv3cfv4DhF5RzOZIvKgiDxp/3lNRO6yj18oItNl5z7NAscpJtB1Flbx/mwCFybK3pjL9+5wOEZBCYVkY9dMdVaXpKeYJkTE765YvxNmwC4KTM+9vTu1RlHZBAEyKN9QU5lF/IMrAEhMHADl3Mq9aHUNZF5jyreyBZlryODBmNwDWMqjWEeSKWSa3q/RaGpxN7tARFzA14C3A6PAYyKyVSn1XNllHwQmlVInisgW4CbgKhHZCGwBNgErgZ+LSDHa6ShTKXVB2bN/CPy47DkPKqUu6/TLLhS6iYGUsDdu5RuoaXwYDXqYIYSZai+N18jOkCCE21U7c7zmPluBqOQEErIq02sq24Hk1CFCgGEH9uu5pMoJD1vV6Nmp/Si1qW4dSL5QYE3hIFOhc5p8S8AwOOReQSBhpfIWlYYgJUtEo9G0RysWyLnATqXUbqVUFrgduKLqmiuA79iffwBcJNZOcgVwu1Iqo5TaA+y05TWVKSIR4G3AXZ19tYVJt/GV4v15O15hVPXBAqtOY0YFId28kLB8b/ZkZ0jY2VLNLCWXXStS3g/LKWNqZtzqqusOWf2zWlGesZEVmEoozBzApE4arwjj8STLGUfFGtSAlK8lsJrhzCgAOdNyvQU9QbKmViAaTSe0okBWAXvLfh61jzleo5TKA9PAcIN7W5H5TuAXSqnyXfAPROQpEblHRDY5LVZErhWRbSKy7dChzqby9QqnIHq3vbCKspLFor9gbdV4yOtmhiCuTHt1IMVW7kDDdu4A3rBldSSqOvJWK8jEpDWb3Wd38G3lqw+Fg0wQgdmx+oWECNOH9+EShWd4fXOhQDa6nlXqIMlMjpwduwl5QuTNfMV8dY1G0xqtKBCnV9HqbaDeNe0eL+dq4Layn38PHKeUOhP4J+pYJkqpbyqlzlFKnTMyMuJ0SRss3KKzYoquJzRYOy/cEJJGGHeusQUClcrNn4+TLSqQJorOF7U78pbNBHFK483YTREDA3Nurma4XQaTEsOTOtSwlUli4jUAIitOaCoTwLVkPUHJsG/0ZfLK6hMW8lgjcIsWiUajaZ1WFMgoUN5oaDXwWr1rRMQNDAATDe5tKFNEhrHcXP9VPKaUmlFKzdqf7wY8IrKkhfW3Tc/SeB16Q3Uju3wzLSoQn13zUR2vSLvCeHKNW7pXb84hcy6rq9nckkC0tiOvk7WQm7EUTHRwWemZrTDjHsKfOVy/El2E7IzV8n1odeMiwiIhu9hwYu+OkgXid1nZW3mzeeNJjUZTSSsK5DFgg4isFxEvVlB8a9U1W4Fr7M9XAvcra/fZCmyxs7TWAxuAR1uQ+R7gp0qpdPGAiCy34yqIyLn22luY29o+vW7nXk3HabxldRa5hBUgD0SHHN1NaXcUf6GJAinf8M0CIZKYLYydBatWxFRCvqwjr5O1YCYOkVcG0ZizoqtH0jNMKDfRsJWJzI6Rx4V/0LmNezUja04BIHVwJzkzh9tw43F5AK1ANJpOaJqFpZTKi8hHgfsAF/BtpdSzIvI5YJtSaitwC/BvIrITy/LYYt/7rIjcCTwH5IGPKKUKAE4yyx67BfhC1VKuBK4XkTyQAraoxdLYqEf6qFxJ5FNWfCMYHcQKOVWS90QIZBrPBKnY8O2sLtNuxNgsiB4LBZghiJmcqzVxUiCSmmCKCCGfF5KtWyAZ/xJiqUledUgNBsuF5ckcZswYYaWr6T9jACLLjyePARO7yJk5PIYHt2Hd24oLyzQVL43NctxwEL/HuXW8RnMs0dJvnu0yurvq2KfLPqexrAane28EbmxFZtm5Cx2OfRX4aivrXSg4bZa9qkovzgIZCIdIqdpYR8EbxU8G8llwe+stsEQ2MYUXq5X73On6m30s6OGACiPlLd0dKtHdqXGmjShhaZ4aXI4ZGsE3mcXMziLe2j5aBgbB3CRT/hU0rwKxcXk4ZCzFF3+1ZIG4xfoVKJiFhrcmMnne/+1HefyVSdYMBbjtQ29k9WCw1SdrNEcluhJ9vnn2LtToox3dWlGoZ88Cifo9ILUbc6nAMF0/E6tcQSTsYHgxPbdZDMTvcTEjYVxlHXmdFI4vO0XCNVBSLC23M4lYMZPUzCFHpWsqiJmTpEOtua+KTPlXM5gepWAWLAViWyDFoHo9/v7eF3ji1Un+6q1rmUpk+dgPn9adfTXHPFqB9JPqCX37n4LffQ3u+6RlGbQtbi5mIfYskKC3jisl0JoCKcortnJvZexskYQRwZNtXIkeyE+R9sxlirW66foHLbtiZvKgYwxkYmqaAZKYA/UHSTmRjhzH8sJ+svkcbilTIA1iIKOTSe54ZDc/WP5v/K+Hz+cXsRt5fuceHnu5casYjeZoRyuQ+WTPr+e214PbuxLlzs6QNKxOvE4pty7bFdWsoWJ1Vpc3PFg63sxaSLuj+Mo68pqYNS6sSGGanG+wbQsksnQdYM9Td0hqmDxgVZR7Rta1JK/E0HpikuDQ+Bguw4UhhlXV3kCB3PnYXv7CuIuzJu+BM7YwkniRr/j/mVsf2tPeszWaowytQBrRIxdFaQM89IIlFoGx5xrcUWc5ZUFqT35uFoiTi6c4VTDVoKV7ucWQtRsvBuz6DqGxCwsg4xkgWNXSvcJaMAtEVZxCYLjiO7TC8AqrODAXH8Nw+GeaPPQKAOHlrRURFgkss1N5D72CSyzrzS3uugqkYCp+se0ZPuLeCqe9G971z8jb/jfn8wSzL9zPbEZnb2mOXbQCcaDXdSBFsvufZ7u5jhxu1Njz7csrS7v1lRX9lZ8vUrQk0g2GSpW7sPKJ4iwQW4E0qUQHyHsHCKlZMM2a5wOk4+MYopDgcEmxtOrCWjoc47CKYiYOO7qw8hNWI4PBFY3buFcztMZSIMnxvbgMW4EY9RXI06NTvC1xD15ycOEnrIPnfoicb4j3cQ+/fGGsredrNEcTWoHMA8W3fHN8F/vMEcZUlIT9Bt0pAXsWCODoHvKVZoI09tOX1paasuo1orGG15dj+mMYKLBbplSn8U7bfbBckfY7AvjcLg7JElyZcUeF7preSxY33oFlbcldsvpkCoA5c2DOAjHcdYPoD744xv9w30/uuAthyYnWQbcP19nv5yLXE/zmqfZfBDSaowWtQPpIxRt8PoM/N8VhooyrKNnxVzuSKVbKFSGVoOCJ1L0uZFsS+WR9BVLRTDFlZXUNBL1zz2lG0B48ZafyVleixyfsPljRpSU3VDuZSxOeEXzZKUcLJJgYZdK1BNWmsWj4ghyQYVypsZICMcSo2wvrteceYoVM4Dnrf1TKOf3duDAJ77kX09TZWJpjE61A5ou49TbujYwwzgCexP62RZQUUi6JmwKmnarrlOEUDkfJKlfTmSBFjMwUccJ43fY/CYfU4Jp7gpabLBMfL62vXPGkpiwFEhxcWlIsJq03LUwEVhDOTzvGQEayr5HyjXTUBPGQbyXBzOFSBlY9F9ZMOsfKQ7/GxIANb688ufx0ZoNreEv+YZ7b37znmEZzNKIVyDwgCMStvk3BwaUkvIOEsoeh0GEAtpiaG6hSIGUWz0DQywwhVIOZIOUbvjs7w6wxV7DXShaWx+6wm7JrSKrTeLN2p97I0JybqR0LJB9egZ8Mkq2sqI8nEqxQh8mHlrcsq5xEaDUDhYmSYnKJi4KqLSR8eNc4F8oTJEZeN2dtFRFBTv0j3mg8x6Mv7q25V6M5FtAKpI+Ub5a5KatXpCe8DNMXw8CEZHutvIouomJqrlGs9XBw40QDbmZUEKlTB1Ldat6XmyblmnOJtaJAfBGroWLKbphYPVDKjFsB5tiSlW2n8QIwYHX4z0xU9u7c//IOPGLiiq3syALJxdYyQIKZGatXmCGGYyX6E88+zxnGHgKbLnWUEzr1YrxSYPy5X7W9Bo3maEArkEb0MI03MW4NMgrERnDZ0/lIHm5wl8NybIshOW0pEFdwsOZ8EZ/bRVxCuLLO7pXqcbv+fJxMWVZXo7hAkeCAZYEUXVhQNeEwMcakChMKBjtq4eJZYg2KmjnwcsXxqb3PIyiCw6s7UiBFuROj1nx0t+F2tEBk58+s86dc4ixo7R+QFw9Lxh6moOMgmmMQrUD6SPmGnprcT1a5iAwM4w1ZCqQw29nAq2JtR7FqvN7mnDLCuJu0dC/eGzTjZD1V0w2b7IlBe8ZH3m7pXjPhMHWISWOwVOxYuqZFBtecDFjdc8vJHnwRQ8HgsrUdDecKrbSyqVKvvQhYLixTmRVre/lwgjNTj5LwL4NlpzkL8gaZHD6Lc9XT7DjQ+O9Zozka0QrEAWkys7sTcvFxpgkzFPLhtVNsE3aWUqsUYwzF1Fx/pKrtSNVemnZHKirFq2WV/UBIzVLwzaXwtlIHEgsHiatAqSNvdRA9kD1M3GWtsVQH0saGv2rFMqYJoMZfrjjuHX+eSQbwhgY66kc1vPZk8hiY+60G0MV6kPJU3t/u2Mf5xjMUTng7NKgL8mx4K5uMV9i+c3fb69BoFjtagfSR8s2tkJxgSoVZEvHjK1aJT7epQOwYQ65YNW4ronrxhaw7ir/g3NK9YthVJo4bE+Wfs0BaiYHEgl6mCZXSeKtnoodzEyS8lTO/2tnwwz4XhxkmMFuZ8jyUeIkx//GlZ7aN28W4a4TI1EsApXTeclkHn/klYUkTOf0PG4oaOOVCAOIv/qb9dWg0ixytQOYJSU0xTZio300gEsNUQna6syrmgv3GX4xBOGVhgTUTJGTOOsZyStcK5OyhUFLWyr2VViYhr4tpFcawO/Iq1FzNhlLEzEmyfkuBdBJEN5VJ3LeUofRo6Vg2k2Ft/lWSsZPbllekYBaI+1ewOr0LpVRJgRQD6bmCyZL9D5ATL3L8hQ1lyaqzyImH4IHH2l6HRrPY0QpkHhAEd3aKlDuCiDAQ9DJFqKMYiCCYqWmSykc0FGp4rekbwEMecqnak2WdgmftNFx3WSfeVlxYIsKsEcFTrEQvUzgqEydAhkJoaek51dc0QylFNricAXMGEtYaR3c9hVfyeFac1ra8IgVVgIHjWM4Eh8f217iwfv/KJOer3zO19FzwNv47xu1jPHoaJ2e3MxZPN75WoznK0Aqkj5RnOvnzM6Tdloso4vcwoaKlTbFlecUgdXEWSMAqhKu7Ofvrt3QvX1ty2gqCu0NzTQ+bzQMpknJH8eXmWpkULZD44X2WzOjykrzy57aCiYkMn4ABTL74MACTLz0CwPDxZ5euaZe8mSe0zAqk79uxbc6FZff0euqpxzneOEDkjMtaE3jcGzlN9vDkrteaX6vRHEVoBTIPCEKwMEPOO4AghHwuJojiSjdutV5NcfM1MjPECeFzWxtf3eaPtkuq4FBMWJKFQXrGWkexf1Zxza1s9hnvIKF8bQxk6pClQHyx5SV55c9tBaUUw2tOxlQw/pKlQOSVh5hUEdactLn0zHYpqALLj9sEwOzuR2osEHnxHgD8G53rP6oZPvUteKXAgecfanstGs1iRiuQBnSb2V96gy/kCKpUKcvJMISkK4I7W3/YUz0sd9g0SaPWtVK9ORddUkmHjrzFjVdEyM5allAxLbf4nIrvUIdsYBkRFYec5b4pWiCz45YCCS9Z1Za86u+zemSI3eZq3PseA6VYObmN3cHT8Lg9xYvapmAWCA+MsNdYTXD/oxVjbQ9Mpzkj8RCHQxtgcF1L8jzr3giAa/SR9hej0SxiWlIgInKJiOwQkZ0icoPDeZ+I3GGff0RE1pWd+7h9fIeIvKOZTBG5VUT2iMiT9p/N9nERka/Y1z8tImd188WbfN+eyCm5iTJWjYDyx0quoYw7ii/XXg+lYpqsNx8n7aptpFi9OXtL2V61CqT82mIr93C5AmnR5WSG7XYi8f2WBWLfl5m03DmxEXvkrMx9h1YxlYnXbbAjdBarpx9n4pn/ZrkaY2bN20qKqhMXVkFZ42zHhs7ixPR2Cnmz9LwHn3qBc2QHnNI4+6qCwCCHAsezOv4U6Vzj2eoazdFEUwUiIi7ga8ClwEbgahHZWHXZB4FJpdSJwJeBm+x7NwJbgE3AJcDXRcTVgsy/VUpttv88aR+7FNhg/7kW+EYnX/iIkLZcSFJWOZ73Rgma7RWfFVuZ+PJxMu7KvlVOFGeCNGrpLghmcrKmlXurFoMrZo2eLVbal9YyvY+08rBkqXW+k268RSvJPOUK3BSI/ecWUsrLuguunlMgHbiw8mYel7hwrz+fqCR56WnLcsirPBOP/ycuUQyf9c62ZGZWnsvr5CWeeqW99jQazWKmFQvkXGCnUmq3UioL3A5cUXXNFcB37M8/AC4S61X0CuB2pVRGKbUH2GnLa0VmNVcA31UWvwNiIrKihfW3TSdtN5wovm0n7RiDOzRcii0UvDH8Kg2FXFvyBCFoJsiVtR2pZy0U60SyDgqkFAMRA9LTTBEmardybySzGv+Q5aKaHnt5Th7gnt3HAVmC3+uuON7Ohl9UNm8676380HwLBibfH7iG9WusOeitxmmqKagCLsPFhvP+mLwymHn8PwEYnUzw+sn/YiKwHlm5uS2Zg6e8magk2fP8422vR6NZrLSiQFYB5e1GR+1jjtcopfLANDDc4N5mMm+03VRfFhFfG+tARK4VkW0isu3Qoc5ahfSK4gaYshWIJzwXpDaLNRcNuuXWCrT+E1SzFHyVRX/lzysSjlnPyzu0dC+PgbgyU8xKCJcxpzhbDXqHl6y1vsbEvor7Qqn9TLiXzsnroJ178drVsRCh99zMJ06+l7d+4HMVMjtN43WLm2BsGS+FXseJY/9NOpvnFw/ez1nGTtyvv6Zh9bkToRPPByC757dtr0ejWay0okCcfpOqf2vrXdPucYCPA6cArweGgI+1sQ6UUt9USp2jlDpnZKT9SXj9IGt3zw1E7RiDmivaU6nWM7EUCsmlrKpxX7R0vJ61EAmFSSkvZqq+BQLgyU6TMCpjKiWZTTbo4eER0spDdmpfxX0DuTHivsp264YYbQW9y5XcJaev5PNX/wFrhoIV8jpSIGahlLob+IMPsZoxfvvT77Fxz7+ScA0QfdOftS2T2FqmPEtZOrFND5jSHDO0okBGgTVlP68GqhPeS9eIiBsYACYa3FtXplJqv+2mygD/iuXuanUdC5LcbLFyfE6hue35Etl4ewqE4mwMf+3o2WoFMhDwWK1GnOpA1JwLy5ubqQnKt2qBLI0GOKgGkVlrQJaIQD7LoDlBJriyRmYnLiyngVIleR0G0Q3DkrnuvKsYHTib/yf1X5xq7MF16efn6mfaQYTppW/gbPUsLx7UA6Y0xwatKJDHgA0isl5EvFhB8a1V12wFrrE/Xwncr6zf/q3AFjtLaz1WAPzRRjKLcQ07hvLHwPayZ7zfzsZ6IzCtlGp/rF9bdNBnqYziBlhITlFQQiQ2VKrw9tgB7sRM68WESimyCXuGRXBOgZQaFVa9jQe9LuIqiJGpny4sCP5CnIwnWnPcSWY10YCbMRnClbDashgYZCZHMVBIbE3FtYYYbW34FXEap7V34MIqtisppu5iGKz+8I+Rc6+j8IdfxH/O+9qSV074lLcyIjO8uH1bxzI0msWEu9kFSqm8iHwUuA9wAd9WSj0rIp8DtimltgK3AP8mIjuxLI8t9r3PisidwHNAHviIUtbgBSeZ9iO/JyIjWC6rJ4Hr7ON3A3+IFYhPAv+z62/fZ0pv76lJpgkRC/mYsGPmvog9zW+mPQskPWvFMzyh2oypakSEhBEm4jATpOQeQggWZslXK5AWg+giwqRnOSszz5V+Pjz6IquAwMhxNdd2koVVL626XYsGKM39KBYPAuCL4Nn0xxTKj3XA0Ka3wS8gs/NX8PaLupKl0SwGmioQAKXU3VgbeLGZ2jUAACAASURBVPmxT5d9TgPvqXPvjcCNrci0j7+tjhwFfKSV9XZNafPsjlIdSHqGKRVmbcDDZM6yQAJ2zUVbLiylyCYsa8JT3rfKViBOm2nKFWEo16CViVkgwix5f+XI1nYK/+KhtQzO/AryWau31r4dAERWnVIjs91K9HrWB9gxkDb/LxVnnxdjIEVchstxKmE7yOA6JtzLGBp7pJRyrdEczehK9HnAnZ0mLhHcLqPU5bZYtJdPtF43oFDkk5YLy1/edqSBtZBxR/Dla1u6lxRDegoDhRkYrjjfTu+qQmy95ZiaPWiNhz30IinlZfnq4yuua2XKYTnV7eGr6SSIXnx+tQJxi/NUwrYQYXr5G3mduZ3nX2u/y4BGs9jQCqSPFDc3T3aGpB2kLsZABsJBZpW/NIypVVm5pOWOCkaH615TTs4TIWA6KBBbMaTtNicSqlIgbVggnpENKCA9sQ8RwTO1m5dZwUgkUHFduxu+idnQAoH2CwkdXVj2z0XrpBuGN13EkMyy/fd6Pojm6EcrkD5S3KT9+ThptxVjKG7MsYCHKcKOGVKNZKmMpUBCLVogBW+UkEqAWbnRFjfy2Ukr+O2JLK04344FEll1CqZAcnwvghCb3cWYby2GUWs9tNvKpJEbqKcuLHF1b4EA0dMvxURQO+5t6frn98/w/m8/ymmfuY/L/ulBHtjR2YwYjeZIoBXIPBAsxMl57dRQsTbvoNdFnFBpGFMz5lxOsySUj2h47u2+kZvH9MdwYUKmMpBe3HhTU1axpT9WpUDa6J67ZsVyXlPDmBN7kNQkI4WDzMQ21VzXiQureqPvRh6UWSB9iIEAEF7K/vBpbJr5DeOzmYaX/vrFQ7z7Gw/x7L5p3n3GMKlMnv9562Pc9cS+7teh0cwDWoH0G9MkpBLk7U68xUByMUPK3SDFtpxS0DsbZ5oQEb+ndK6Ru0nZsQ2VrIy1lFxidhpxZHBZxflWCwkBThgJ8Yw6nsDMbtKvPGHdt7K212UnCqSRC6uTVialGEi1C0tcKFRnI3KrcG/6I04z9nD/o0/UvWbHgTjX//vjvDP6Ar8b+gz/55mL+Llcz8eXPsLHfvgULxzQtSSahY9WIA3pLg9LKQW5BAYK5beypso3+5Qrgi/f2kZR3ChdmRkSVLUdaeBukpAVrE9NVbpGitcWpyJGl1RWjbcTA3G7DA6ENxLKTZF57LtWE8WT31S7ljY3/KYKpIM6kKILq1QHYlM91rYblp17JQCJx293PD+ZyPLn332MP/H8mhtnP4PHzMJbPoYMrufa6f/L/+v5EZ/60faOquw1mvlEKxAHetnO3czMIoBhd+It3+yz7ij+fIsdee29xJubIe6qrJRutNl77fYpiYkDleLsa1VyghkVZDgarjjf7gCo2bUXAjB84EF+rTaz+YTaPped1IE0dGHRhQurygJxG5ZCKQ6V6orhE9gfO5u3xO/miaruvPmCyV/e9gSnzfyWT5s3w/Fvhet+A2/9BHzgv2Dze7lO3Ulk7/385Ok+18lqNF2iFUgfUUqRSUwjgCtUaYEA5HwDBB0ypBxl2Ru5Lxcvjcatd005wZjlmkpOO1sgrtQEU0Tweyo31HZH0L7utE3cVXgTLxTW8vMVHyLorS0xanfDbxZEb2VuezUFVUCQGsumlxYIwOAFH2K9cZCf/eS2CqX5xft2kNv1IF/1fsXq+HvVv4PHb500DLjsy6hlm/gH37f49i+e0n21NAsarUD6TDphuag8YcsSKC/6K3gH8JMpTfNrRHGjDBfiZLyVfbAauXkiQ5YlkK2jQLyZKWZc9Xs/tWoxbF47wPbBi9mS/yyXXeRYC9p21lRBFXofRDcLNdYHzFkkvcjEAvCf+S5m/Su4dOxbfOvXO8kVTL503w4eevDnfCfwj7iG1sN7fwC+SssPtw+5/J8YUlNcPHEb9z17wPkBGs0CQCuQPqJQpbqNUife8sl8AbuaPN08E0spBUpVBOSLNAp4xwYGSCsP5mxla/vixuvPTpJ0D9bcV+qv1eKGL6K44dJT+PlfXcibT3LugizSfjPFRhZIJ6m3BVVwVLhFRdWLWhAA3D6Cf/h3nG68jPmzz3Lq/76H3//qLu4IfAFfdAn86X9CcMj53lVnY572Hj7kvoc7fvlYb9aj0fQBrUD6TC4VR4BAbEnNuWJcJDfbvBrdVCbkkrgpoAKVG37JqnFoVDgc8TFOFJWsapli64VIfpqst7azbyczzN0ugxUDgbrn2w2iN7NA2lVIYFkg1QF0mIuB9MoCATBOvxLzrA9wnfunPBP4C/7D+3mCQyuRa34CA6sb3/vWj+MhzxvH7uCJV1srNtVo5hutQPqIUgozYymQqK1ASm+/Ctwh6w00Md1aR9582pIlwdarxoNeN1NEcaWr0nixLJqomibnq30TbjcGUtzIm/WuatsCaVDj4hJX2wokr/LOLqwex0AAEMH4o/8P3vUtAmdcAe/4PFz7AAwe1+xOGD6Bwql/zJ+6fs6dDz7TuzVpND1EK5AGSA/ilyozS1wFiIWsN/Nya8ETKabYNp+cqFCkZqcQ5hTP3ELnrnFi1jWAL1NpgSgUKjOLjzz58PKae9q1QIpv7r1Mu201BtKLzK6iUulJFlY5InDGn8AVX4M/+Ah4Qy3f6nnLXxOSNMte+C5jM83jZBrNfKMVSB9RKCSbYJoQ0YBV+Ff+Ru23B0xl4s0tkLmMLoUnUhljkDkN4kjKE8Ofq4yzFFSB9Oy49Q8g4qBA2rRAmg1/Kp5rywKhcTfeoiJoR2bezJfcVU6yemqBdMvy00gedxF/atzLHQ+9eKRXo9HUoBWIA9KkgV+rKBRGbpZZwhWFf2BtuCFbgeTizWMgCkUuMWPHU9prO5LxDhEuVFW8K0hNWYrLP7iy5p7ixt3q5tyKC6udtFulVNNCwnbXCPWD6IYYVifhLmIgyVySJ8ee5HCq9SFhzQi+9a8ZljjJR28lk19Ayk2jQSuQvqKUwp1PkiqbN17+Zh8biJFVLgottHRXKPKpGVAQilVZIE2shYJ/iCApyKVKx0xlkrU78YaXrKm5p7jJtrqhFjf7pnUbLbqbWlFIJauhjU3fNOsXJ7qku468L8+8zFRmihcnX+xdFflxb2Jmyeu4Or+VrY+/2huZGk2P0AqkjygUnkKSjGeuzqI8tjAY9jFFpDZDykmWHZBXSkoBeSeZTphhq5hQxedqCkxMCgkru2doWa0CaeYWq3lGk9kdYLmwFKqlzbXVoHz5ta2QV/n6CsTovCOvUoqJ9ASGGGQLWeK5FjsMNEOEyEV/zVrjEDse+HddWKhZUGgF0md8ZpKsx7lQL+R1MUUYI93iTJDMLLMEGQxXpso2s0AkarmoUhOjc7KUgsQE0+r/b+/MoySrygT/+957sUfua1XWvldRgCwKiguK2tgojLYcwdbWUQf1wNg69rTaPS6tYuNpR7THhfaAuHQpreAoMqIgasu+71VWUQu15Vq5x768O3+8iKzIzBcRL6IyKxK5v3PqZORbvvdFZMX97v22G6a7vXwdiNc9zL3s3VFLXKX43EpB9FqL/2xlYyvbNQZSfFa9zRTj2Tg5O8fa5rUATHpskukF2XwR09G1vDX+U/5Tt3vXLCG0AVlMbJuQSpEvKfwrHURFhJjRjC9d3YAopTAz00wRIeSfP6hWqrEIdTg1B1NDB4+rpmzM9AQj0u4qrx4XVrUeYsUAuxeZxWB2tf1Ais/2QlGmWxovOLUg9bqwUnknS6o12Irf8BPLeGtR4wnDIHT+xzjFOMj9d/504eRqNCeIJwMiIheKyG4R2Ssin3Q5HxCR/yicf1BE1pSc+1Th+G4R+YtqMkVkR+H4MyLyXRHxFY6fLyKTIvJE4d9nWGRq7bM0Fzs9hYGC4HEDMlPhXXDjJK0W/NnqHXltbKzMFFPm/NUCVI4vtPSsASBx7PAsecH0OFPW/J0NS3WtJWZRabUAYBiz33slip/9QsZAiim6PsPnev5ENpVKFuJLQStI1B8lll1AAwJYL7mMeKCLV4/8iHueW7ggvUZzIlQ1ICJiAt8E3gRsAy4XkW1zLns/MK6U2gBcC3y5cO824DLgFOBC4FsiYlaRuQPYApwKhIAPlDznbqXUSwr/Pl/PGz6ZZBNO6qwRctl+tjBApn0thOdmSLmglMKfmybl0nYEmNlr3Y3uzi7iKkBu4vhGRUopWrJjTAbmZ2B5kemmXzUXVi0Dvpe6klpXIDOt3Mu5sE5gU6lELoHP8OEzfIStMMlcsvpNtWAFCLzySl5pPsutt/5Ux0I0SwIvK5CXAXuVUvuVUhngJuCSOddcAny/8Ppm4AJxfA+XADcppdJKqQPA3oK8sjKVUr9SBYCHgMo9HxaBhWrnnphyguNWdH6ld9GA5AJtRO0pqDJQKxRhe5pMwH3FYIhRNl7R0xJiULXDVP/MMTsTp1lNk4qU/3hrSWstdrmtRC0DfnEgd2s7Uo88KL+dbRFLrLoLCZO5JCHLiU2FfCFsZZPOV96RsFascz5IItTLeye/zY4H9i+obI2mHrwYkD7gcMnvRwrHXK9RSuWASaCjwr1VZRZcV+8GSjeXfrmIPCkit4vI/D1TlxiJ6VEMFP6m44P+3IwpO9SGnxxk4hVl2dkUIZXGDs3vqQWVXVh+y2DU7CSQGJo5lp7ox0AhFdpq1FK3YSu7bGyhSE0GpMy+HfXKK5VZKYhe7wpklgEp/ExmF3gV4g8T+ssvsc04yKFff43dgwuU6aXR1IkXA+I2rZw7qpS7ptbjpXwL+KNS6u7C748Bq5VSpwP/B/i5q7IiV4jIIyLyyMhI9RYhi0kq5gTHg80lBmTO6sYoNEasVguSnXYGf4l0u56v1mdq2t9NU+b4BkXjQwcRINyzvuw9tfSuUlTuW1WUB7WtQCpmYdUaA6niwrIMi7zK15yJVVxtzDMgC+3GAmT720ivfT3/0/gRV1//Y/YOL2ysRaOphfL+geMcAUoLBVYA/WWuOSIiFtACjFW5t6xMEfks0AV8sHhMKTVV8vpXIvItEelUSs2KKCqlvgN8B+Dss89uqKM4Mz1KCIi2H9+db+4KxIw6xiU+MUxze/nVwNSoM/j7Wua3HYFCm5AKKbfTkdW0j90F6WkINBEfO4IA7Ss2lr3HwHsQPW/n8Zv+itfUMuCX27t8ln5Fg2TXGAMp4xab6chr5zFM7wmKqZyTgVU0HEEziCHGohgQRAi8/Ttkv3Ue18a/yPu+Aeee+yrOXtNO3lbsG4mxa2CK3YPTHBxNgMC6zghvPKWXd5+7mq6mwMLrpHnR4sWAPAxsFJG1wFGcoPg751xzK/Ae4H7g7cDvlFJKRG4FfiQiXwWWAxtx4hpSTqaIfAD4C+ACpY5PBUWkFxgqyH0Zzuqpegl3A8nGjxFQQltHz8yxuW1H/IWGirGJEZoryIqNOwYkVMaAVKvyzrdvgDHIj+zFXHEG+fFDJFWQzctXlb2nlnbpClU9C6uGFUgxFrGgWVjFGEgZo1TMzsraWXyme6aWG4lcAnAysMD53IJmcHEMCECkA99//SWtN17EzYl/4Ib7LuS6u88kTogmEmxuSvKepjRrVscJ5sZ5LtXKjb/fwvfuXcPnLj6Ft5150sOKmj9TqhoQpVRORK4CfgOYwHeVUs+KyOeBR5RStwI3AD8Ukb04K4/LCvc+KyI/AXYCOeBKpZxvu5vMwiOvAw4C9xfcPT8rZFy9HfiwiOSAJHCZWrB+EbOp5orxTHKcOGE6moLHZc8ppgsV+mElq7R0T0w4LqymTvcvfzV3U2T5NtgLY4eeoWvFGQRihzhCD+dEyq8aTDG9FxJ6qAOppWGhFxeWiNTkZsvZuZmeV24UVyBZO+tJXpGioQj7wjPHQlZoxrAsCp0bMT98D+btn+CDO3/Oh6xfHj+XoTC1Egi2cHZqgsv98FjgpXz0p+9k9+DL+eSbtixYsojmxYuXFQhKqV8Bv5pz7DMlr1PApWXuvRq42ovMwnFXnZRS3wC+4UXfBeME7ZOZHmeKCD4Xd0jR9kXanBVFdmpo3jWl5GIjmEBbp3vabbWBtHvNVvL/KUwf3UWXUjSn+9kbOGvBeleVa1I4Vx54rEQv9NaqJtMQw3PmVF7ly8Y/4LgBqbWYMJlLYhnWrPqSoBVkIl19p8kTItoNl96ITF8DA09CLuW0i492Q6Qbwh1gWhAbhsd/yBl3X8udoX/gI/d8iH/KXcpn37JNGxHNCaEr0RcRKzM5q5EizB9Emzt6yCoTe6ry3tf52DDTBOluc2+LUm23v/W9HexTy7EGnsAePYBfpci1rq34zFpm914KCWtxOeVUrqrxgELg22PmVM7OVUwLLnVh1UJpBlaRoBUkr/Jk87XJqoumHtj0Rth2MWy4AHpPdY6Zhfca7YZXfRy58kH8y07h3/zXwoPXcc2v/7T4umn+rNEGxIWFmpUF81NkfXMMyBz3WGdTkGO0QKzyCsSXGGaCFkI+9wGw2mDfEvbxjLWd7onHGX36N9giWH1nV3xmLft3eKpEryHo7UUe1NZBt1pso2hA6lmBzDUgxd8X1Y1VKy19yHv/H2rLm/mc7wdk7/km3/7DvkZrpXkBow3IIhK2p8n5Z4fG52ZhBSyTMWnDSlQ2IMH0MAlfa1nj5mW1MNnzMoJ2gra7P8sxFWX9xu0VrzfE8OxuqqUOxGsvLC8GxGf4PAfRs3a2bBsTqM+FZSubdC49z4CELSceUszQWjL4gsil30NteQuf8f2QoTu/xo8f0m3iNfWhDchikY7hJ4vyz2494paJNG11EEpXSChTiub8KBn//Ir2Il7iFZHTLmZUNWHaaR60zmRVR7Ti9V6zsLy0HSmeF8RTYL7c3uVzMcX07HLK5isbEBGpSR44BkKhXF1YsDi1ICeM6UMuvRF7s7MS2fmLr/Krpweq36fRzEEbkEUiPTmEAoxQ66zjboHkZLCLpmx5A6KS4/jJoEJdZa/xMjBfeMZa3sfn+GT2/eS3XlwxoAzemwt6aTtSxBDDU8yi3NazcykW/3nBS3quz/TVtAIpGoi5BsQQA7/pX5oGBMD0YVx6I/mNF/IF343c+x9f4e7nGlt4q3nhoQ1IRerPwhoffB4lYM3Zv9xtBZILddGsJiHvPnBNDOzHFrCa3KvQwdteFs1BH1dfcSmnXXwlb9zeV3WG77W1h5fNn2Zkety0KWfnKq4WSuV5GfCzdtbZ4KuKTEtqa+lezoCA48Yqtnlfklh+zHf8gOy61/MF6wZu/+H/5tGD1Tc302iKaAOySEwPPY8NRObsN17cE2PWYN/Ug4EiP+0eB5kc2EceCLYucz1fxEvK7fa+Fv7qrD5Mw/DUfr2Woj8vLief4W2GX9MKxIORK2ZDVTMgPtNHxs5UlVckmUtiiulahR+0ggvfD2uhsQL4Lt9BbvWr+aJxHT+//mp+9tiR6vdpNGgDsmikR53AZGvXbAMy48IqGeytVueayZINn0qJD+7DRmjpLF81Xstuel5jFpZYKFTVAboWF5bXTZs8G5CCjtVkFuMa1QxIwAyQyddmQIrxjrmErBAZO1N3g8aThi+I/69vIrfmfL5gXk/6Z1fxke/fzfPHKjf41Gg8FRK+2BAPrphq2JNHGCdMR+vsGIjbVrHBTqceY2pwP+2bXzFf1shuJgizvMe9jQl4Xy1AbUFvKKTUUn51MSPP8ODCErPqAJ2zc467Saq7sGb6V6k8VoX/zjMGpEoMxG/4yeQzKKU8pXMncgma/E2u54KmY1hS+RQRI1JVVkPxh/H/zS3k7/oCl997La/d/yT/eu1bGV5zMeduWcUpy1tY3RGmtzmIYejiQ42DNiCLhD92lHFpocs3e8ByS2XtWLEJgOSQe05+eGofu8weNvnKD36mmCjUTAV3JYp1GNVcWDNprSqHj/LPLs7+vcQsLMOqGliu1jXXVUc7R8As3yiwaLQqXQPgN/0zK5pqxqaYwtsT7nE9H/Id78ob8dVuQGxlMxAfwBST3kj5ycOCYZiYb/gcbLmIjtv+ji8N3UDy6A7uO7SVe+0N7FC9HJZlZFtW09beyYrWMNuWN3Pehg7Wd0V1VfuLEG1AFolwapCpQJvrYD53q9iVvd0cU83YYwdcZXWmDjLVeoa3xoIeOsnWswKpRHF272XA91L4V4sB8RtO7CGTz1QcpFP5FILMXF+OooFJ59NVDUgyl3Q2+rLCrudPtK37n8b+xHBieEbG2pbKnQMWjJUvxfehP8ChBwg9/VNes/f3XDBRshd7AqaSzRw60s2uJ5bzTfsUDrWdwxtedhpvP2sFnVHd8ffFgjYgi4Ft05XtJ9Gy3nWWP7foL+Q32Wf04J+eX9CVjx2jRU2RjS6vuGLwuloA7yuG4vO8xBcE8RQD8VL4VwxiV2sPX3pNNbdYOpfGb/qrzpKL8tL5NFEq18nMZGD55mdggfNevay43BhPjTOcGGZ182qSuSSHpw/TE+6Z1bBxURGB1S+H1S93BolMAsb2wdh+GH+e5rEDbB8/wNb+p7g09Ufs+L9x92+388W7Xkv72X/Ff3vNJpa1uH8umj8ftAFZBFJjhwmRJt/c5zpgubUImQj0sTE1vzdR/877WAlYPRs8r0Cq4TWgXBpfqEQx4O3FhWEaTrA/b+fLZm0VjYEXA1JcJVTLnErn01XdV6XP9BJIT2SdNiVuKbxF6m3rfmT6CH7Dz+rm1WTtLMeSxzgaO8rGtvL7tywq/rDTY6v31FmHTduGoacxdv2SVzy2g9fEvs6RR3dw/cMXkT39XXzgtdtZ1XGSjJ7mpKOzsBaB4QNPAeBrX+G6anDbKjbRvJ6e/CAqNTXr+OTeB7GV0LZys6fNlbzWWED1tNuZ5oJVGgJ6zZiC4y6nStXexb3EvQz4PsOHIUZVHdP5NAGruryAGUAQT/Ub09lpAmagoiGO+qPEq2xXPJdkLsloapTeaC+GGATMAF2hLgbjgzX36Vp0DAOWnQ6v+1/4/sezcPlNdPWt49Pm9/nY02/ltq99mH/acRf7RvTOiX+OaANSiTrrCKcPO1ubRLpWus7K3Sq8zb7TARjd99is49bgY+yX5XS0NXlqRe7VgJhiVo2BeO1OW63H1CyZhRVD0Ui4kc6nsQzLU2FiUc9K8pRSpPPpmayoShhiELACnpogxjKxshlYRaK+KBk7U1G/uQzEBxCE5ZHjKeB9TX3kVZ6hKj3TGophwOY3EbjiTnjfHYQ2vIoPm7/gU3su5al/fQdXX/d9bnvyKKnsEk9r1nhGu7AWAbv/CYZUCz3d3a4rELfq6e7N58DjMPzcw3Secr5zMJ9jxdSTPBh5JaZQMa21VheW1wwnQaq7h3JpIn5vWUbFQbySiyiTz3hafczItIIVB+hkLomtbM+ZUGErPOOeKkfOzpHMJatmRxUNzHRmmkCo+nuylc1AbICOUMes+pJmfzNRX5T+WD990T4P76LBrDqH4LtvgrH95O/5Fhc9sYO3Dt7D4Vu62HHLOYx1nUvTpvNYu2I5G7ujrGoPY9WwjbBmaaANiAsnmo7YMf4UzwW24LfcZ/lu1dMb121gRLWgjjwyc2xo1z30ECez9nxCVM5KqqWTrNcVg4hgGVZV91Aqn6LdKN/osZTSLKdyJHNJT6uFIiErxFiqfAuO4mqiUqyilLAVZiI9UbEWZCrjuBqjvsqB9qgviiBMZaboDHVWffZIYoSsnWV5dP7GYX3RPnaP72YyPUlLwH1fmCVH+zpCF38F3vhp7J2/JPLIT3hv/+2Yo7dh3yf008F+exn3sIyp0CpyLavwda6jadl6VnZ3sqYzwoq2kOumbJrGow3IApOeHGR5/iiP97yJCLimjVpikbZnD6BBv8WDwTM4ffR+sPNgmIw88GNalcWqs1/PCP0V00pLs4eq6pif3368ktxKLqysncVWtqf4AjguLEOMsnoqpUjmkrQHvRkkKFR85zNlYzGxjON/97oCifgi2MommUuWzXqaSE1giFF1IDcNkyZ/E+OpcfAw5vfH+wlZIdoCbfPOdYe72Te5j6Oxoy8cA1Ik2IJx5rtoP/NdkInD0UfJ7r+PyNFn2T66l3Ni9xDIJGAE598uGFatHFTdPKF6mAgsJ922gWDfdnrXnsLWFZ2sbAvrosYGow3IArP/3lvYCoS3vRZwr3y2DMt1G9bMujfSuusPDD/9O9o3vpTVR37JQ8Hz2NrTzsix/opZSYYY+A2/ZwPSGmiteh04g3OlLCIvmUhuMsu5iGp1N8FxwxDPxl0H1snMJBFfxHOgvyhjMj1Z1oCMpkZp9jd7ktkR6uDA5AEy+UzFv+FkepLJ9CTrW9e7x84Mk95wL/3xfs9ZZaUMxAboj/eTtbO0BlpZ1bTq5KUFl+KPwNpXE1j7ambegVKQGIPx51Fj+0kM7UOG97Fm7ACbp/cQzdyDMaJgBLKPm+xXy7hDVjLZtAHVvY3mlaeyasMpbFzWQsCq3pNNszB4+kaJyIXA1wETuF4pdc2c8wHgB8BZwCjwDqXU84VznwLeD+SBjyilflNJpoisBW4C2oHHgHcrpTKVnrGUMJ69haOqi02nnc1zk39ydRUV+0HNdZFsfd1ljO68huTtn2b/vavYRBzjvKtmNiWq5tbxm9UNSDqfJmfnPA/4RfdQOXdOLOvM7qu5ckqJ+qKMp8ddz02mJwHH5++V4rVurp2cnWMyPcmySOVGlKWEfWH8hp+x9BjLovPvm85ME8/GPafUdgQdAzKUGGJl08qy1x2YPIDf8FfUdXl0Of3xfp6ffJ7N7Zs9PT+dT7N7bDdjqTGivihNviaGE8MMJ4ZZ27yWFU0rGl9FLgKRDoh0ICvOIgLMmkJkUzD6HJmBnUw8/yTRgWd52cQe2mP3QwzYD6k/+NinlnPUv5bpYC8q3IkR6cRq6sbX3EWwqZ1wcwdNLW20RgK0hvwEfUbj3/sLmKoGRERM4JvAG4AjwMMicqtSOFH6dQAAC3NJREFUamfJZe8HxpVSG0TkMuDLwDtEZBtwGXAKsBz4rYhsKtxTTuaXgWuVUjeJyHUF2d8u94wT/QAWkkOP3cHm+KP8dsWH2Go4MQ63GWfICjltMPLpWYHSFV0d/GLj3/OW5z6NkdrFba3v4qJXvYE943uwDKtqXUTY5/juKzGdmQaomj1UpOjOiWXdM44mUhP4TX/ZhoJuNPmbGEoMkcgm5s2AjyWP4Tf9Nc2M/aafiC/CaGqUVc2zG06OJEawlU13uHwrfDc6w50MxgddVw0Hpw5iiOFZZtQfpSXQwuHpwyyLLHNdtRyeOsxEeoJNbZsqrmrCvjB90T4OTx+mI9RRNa4yGB9k78RebGWzsW0jyyPLERHS+TTPjT/Hvsl9HEsdY3Pb5rpXI1k7iyXe6oDqxheE3lPx955K9xklX/tMHHt4N6MHnmDq0NMER3ZxVmwnrbE/YsTcOyjYSogRYpQwMcIkjCgpM0rWaiLnb8YONEOgGSPUghVpwx9pI9jURripg2hrO02tHTSFQ9rw4G0F8jJgr1JqP4CI3ARcApQakEuAzxVe3wx8Q5xP9xLgJqVUGjggInsL8nCTKSK7gNcB7yxc8/2C3G+Xe4by0sO8XjJxho4+DxS75xYeVfJa2TbZVJzhnXez9sl/4QC9nPq2jzOWHsUyLFc3Q3H2H8vG5g28b3nnVfz+3peQSyV4/fkXkFM5xtPjnmbkLYEWhhPDTGemXQf7vJ1nIDaAZVieDUgxFtEf62dT26aZL42tbKYz04ymRmua3QN0hbvYN7GPg1MH2dS2aaa4cCQxwmhqlNXNq2uSB9Ab6WXfxD76Y/0zvakmM5Psn9xPk7+p5pjBiugKBuOD7BrbxcbWjfhNp8ni0dhRjiWPsa5lnefUZYB1Let4fPhxnjn2DOta1xG2wiiceM9AbICB+ABdoS7X4Plc1jSvYSI9wc7RnaxsWkl3uJugGcQQg5zKkclnmEpPMRAfYCozRbO/mS3tW2YZiIAZYHvn9hkD88jQI/RF++gIdhD2hWfSqJVSM1sWp/IpUrkUyVySeDZOIpcgkU3MpI4HzABRX5SIP0LEihDxRfCb/rIp2bayyav8TGGprWxs7Fltfoq7WBbvF5GZYyKCWAGMvpfQ2XcGXaUDum1DagIVHyE1OUxyYojk9DiZ2Di5xDh2chKVnMLITNKamcKXHSWQPUgoHSMylcCoksOfUAGmCROTCCkJkRcfecNH3vBjiw/b9GMbPpThR5l+MP1g+cEMOK8LvyvTj5h+xAqA5cewAhiWD8MwsSwTw7AwTAsxDAzTRAwLkeOvDdPEMMzC785PQ0wMyyIYjtLaVj1x40TwYkD6gMMlvx8Bzil3jVIqJyKTQEfh+ANz7i3mILrJ7AAmlJoJEJReX+4Zxzy8h5oQEXIYtE3ezu5/v33muJr56T7zuCPSgzrvoyxPOLa1XIpnS6AFv+HnmWPPzPtiKaUw1yhMAjwweM/M8Q2tG6rq3R3q5sDkAR4denSe3NLK93Ut6zzXWPhNv+M2ifUzGB+cZUCK5+fO+qsRMAOsbF7JoalDDCWGEI4XVkZ90YpunnL0RfsYTY6yZ3wPe8b3zNJ/a/vWmuWFfWE2tm5kz/geHhp8aN6zatWxJdDClvYt7Bnfw2NDs2t9BKEv2sf61vWeZJmGyWldp7FnfA8Hpw5ycMp9G4CQFWJT2yaWRZaVnS33RnppC7axd3wvR6aPcHj6+Ney9O8yF7/pJ2yF6Yn0EDJD5FRuxrCMTY3Nu09cvjPlZNdL0ajMIwj0djv/vGDbkEuSTcVJJiZJx2NkUjGyiRj5dAw7HUNlkhi5BFY+iZVPISqH2CkMNY2h8hgqh6nyGLkcVjaPSR6fymHK/FVR6aewkBUyxwJn8raP3l79whPAiwFx+5839y9f7ppyx91Gr0rXe9UDEbkCuKLwa0xEdrvc55VOajJQQ8D7TuBxNVOjfg1hqeu41PUDreNCsNT1gwXX8dfwsbrdbJ7cAF4MyBGgdLq1Augvc80REbFwEhbHqtzrdvwY0CoiVmEVUnp9uWfMQin1HeA7Ht5XVUTkEaXU2QshazFY6vrB0tdxqesHWseFYKnrBy8MHefixY/xMLBRRNaKiB8nKH7rnGtuBd5TeP124HeF2MStwGUiEihkV20EHions3DP7wsyKMj8RZVnaDQajaYBVF2BFOINVwG/wUm5/a5S6lkR+TzwiFLqVuAG4IeFIPkYjkGgcN1PcALuOeBKpZyIm5vMwiM/AdwkIl8EHi/IptwzNBqNRtMYRE/iyyMiVxRcYkuSpa4fLH0dl7p+oHVcCJa6fvDC0HEu2oBoNBqNpi50hzKNRqPR1IU2IC6IyIUisltE9orIJxutz1xEZKWI/F5EdonIsyLyt43WyQ0RMUXkcRG5rdG6uCEirSJys4j8qfBZvrzROpUiIh8r/H2fEZEfi4j3cv/F0+m7IjIsIs+UHGsXkTtF5LnCz/mdIBuv478U/s5Picj/FRFvzeBOoo4l5/5ORJSILG4V4AKgDcgcSlq3vAnYBlxeaMmylMgBH1dKbQXOBa5cgjoC/C2wq9FKVODrwK+VUluA01lCuopIH/AR4Gyl1HacZJOlkDjyPeDCOcc+CdyllNoI3FX4vZF8j/k63glsV0qdBuwBPnWylZrD95ivIyKyEqfF06GTrVA9aAMyn5nWLUqpDE5jx0sarNMslFIDSqnHCq+ncQa+JbXLkIisAC4Crm+0Lm6ISDPwagpZfkqpjFKqciOxk48FhAp1T2Hm11+ddJRSf2R+/dUlOG2HKPz8LydVqTm46aiUuqOkw8UDODVmDaPM5whwLfD31L0f6slFG5D5uLVuWVKDcykisgY4A3iwsZrM42s4XwT3jnaNZx3OzhM3Ftxs14uI9x7yi4xS6ijwFZyZ6AAwqZS6o7FalaVHKTUAzuQGqK1z5cnnfcDi9vioAxG5GDiqlHqy0bp4RRuQ+XhqmbIUEJEocAvwUaXUVKP1KSIibwaGlVKPNlqXCljAmcC3lVJnAHEa73qZoRBHuARYi9PJOiIi72qsVi98ROQfcVzAOxqtSykiEgb+EfhMo3WpBW1A5uOldUvDEREfjvHYoZT6WaP1mcN5wMUi8jyOC/B1IvLvjVVpHkeAI0qp4srtZhyDslR4PXBAKTWilMoCPwNe0WCdyjEkIssACj+HG6yPKyLyHuDNwF8vwS4W63EmC08WvjcrgMdExL0j6xJBG5D5eGnd0lAKrfJvAHYppb7aaH3mopT6lFJqhVJqDc7n9zul1JKaPSulBoHDIlLclekCZm9R0GgOAeeKSLjw976AJRTkn0Npm6HS9kNLBnE2sPsEcLFSyn07zAailHpaKdWtlFpT+N4cAc4s/D9dsmgDModCoK3YZmUX8JOSNitLhfOAd+PM7J8o/PvLRiv1AuS/AztE5CngJcCXGqzPDIWV0c04u3I+jfNdbXiVsoj8GLgf2CwiR0Tk/cA1wBtE5DmcDKJrKslokI7fAJqAOwvfl+uWoI4vOHQlukaj0WjqQq9ANBqNRlMX2oBoNBqNpi60AdFoNBpNXWgDotFoNJq60AZEo9FoNHWhDYhGo9Fo6kIbEI1Go9HUhTYgGs1JREReWtiTIigikcJ+H9sbrZdGUw+6kFCjOcmIyBeBIBDC6cf1zw1WSaOpC21ANJqTTKHH2sNACniFUirfYJU0mrrQLiyN5uTTDkRxejM1fJtajaZe9ApEoznJiMitOG3u1wLLlFJXNVgljaYurEYroNG8mBCRvwFySqkfiYgJ3Ccir1NK/a7Rumk0taJXIBqNRqOpCx0D0Wg0Gk1daAOi0Wg0mrrQBkSj0Wg0daENiEaj0WjqQhsQjUaj0dSFNiAajUajqQttQDQajUZTF9qAaDQajaYu/j/3jRqTdReSZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    plt.plot(x[0].detach().numpy(),data[0,:,0].detach().numpy(), label=\"data\")\n",
    "    plt.plot(x[0].detach().numpy(),y2_.forward(x)[0,:,0].detach().numpy(), label=\"fit\")    \n",
    "    plt.plot(x[0].detach().numpy(),prefit_prediction[0,:,0].detach().numpy(), alpha=.3, label=\"pre fit\")    \n",
    "    plt.ylim(0,np.max(1.5*data[0,1:,0].detach().numpy()))\n",
    "    plt.xlabel('x')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## It works!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}