{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates the CUDA-accelerated likelihood function from the heron package.\n",
    "This likelihood function incorporates the variance from the waveform models in heron as well as the variance from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "\n",
    "from heron.models.torchbased import HeronCUDA\n",
    "from heron.matched import CudaLikelihood\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "generator = HeronCUDA()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a smple demonstration, inject a signal into Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = 5e-19*torch.randn(1000)\n",
    "signal = generator.time_domain_waveform({'mass ratio': 0.7}, times=np.linspace(-0.02, 0.02, 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from elk.waveform import Timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "detection = Timeseries(data=(torch.tensor(signal[0].data) + noise), times=signal[0].times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbfd62a7630>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD+CAYAAAAqP/5ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOxdd7wcVfX/3tm+r/eSl7yE1EcqISShl4BKRyCRJiCIyE9ARUTEAop0UFFQsNMUEVS6QAKBhNBSSH3p7ZV9++q+sn137u+PKTszOzM7215J5vv5JG+n3HvP3Llz7rnnnkIopTBhwoQJE2MfzEgTYMKECRMmcgOToZswYcLEIQKToZswYcLEIQKToZswYcLEIQKToZswYcLEIQKToZswYcLEIYIRZeiEkHpCyHpCSIgQYs2g/JmEkO2EkNWSc0cSQj7k/92dW4pNmDBhYvRipCX0XgBLAHycYfmPAcxVnPsmgB9SSo8HsJgQUpoFfSZMmDAxZpC2VJxLUEpDAEKEEAAA4X78DsB0AEEAV1BK+3TK9/HlpKd3ACghhFj443DuKTdhwoSJ0YeRltCVOAfAQUrpaQAeAydtp4t3APwGHGP/iFIazCF9JkyYMDFqMaISugqaAFxCCPkiONo+IoTUAnhecV8HpfQSjTruBrAMwDoALxFCJlJK9+eLYBMmTJgYLRhtDH0HgKcppY8AACHERimNAjgljToIgF5KKUsI6QdQlHsyTZgwYWL0YaStXGyEkOXgNjbfAtABYCIh5F1CyLsAzkxRfgFffhYhZDkhxAngAQDPEEJWAYhQSjfn+TFMmDBhYlSAmNEWTZgwYeLQwGjbFDVhwoQJExnCZOgmTJgwcYhgxDZFV6xYYep6TJgwYSIDLFmyhKidH1ErlyVLlmRUrrm5GU1NTTmmJnuYdKUHk670MVppM+lKD9nQtWLFCs1rpsrFhAkTJg4RmAzdhAkTJg4RmAzdhAkTJg4RmAzdhAkTJg4RZMXQCSGLCCFrCCGrCSG/Ulyr5z0+1xBCTs+OTBMmTJgwkQrZSugHAJxGKT0BQDUhZLbk2u0AfgLgCwB+nGU7JkyYMGEiBbJi6JTSDj6mOQBEAcQll2cDWEMpHQIwSAgpzqYtEyZMmDChj5zYoRNC5gCoopRuk5y20ESgmH4ApQAGpOVmzZol/l66dCmWLVtmqL1QKITm5uasaM4HTLrSw6FAl/3Nf8H23mvwP/hUnqnicCj02XDicKMra4ZOCCkHl4xCyY1Zye9iAD5l2S1btmTU5kg7C9BoBGwoCEtRiez8SNOlBZOu9JAOXS23vgIAw/Ych0KfDScORbra29s1r2W7KWoF8CyAWymlHYrLmwghxxJCCgAUU0oHkmsYm+j62XfRfklmXq4mEnA9dDv6n/n9SJNhwsQhg2w3RZcCOAbAg4SQlTwD/y1/7UEA9wBYDuDeLNsZVQhv+CSn9cV7uzH0zis5rXMswOJtw8Dzfx5pMkyYOGSQlcqFUvoPAP9QnP6Iv9YK4LRs6j9c0PWz7yK6uxmuo4+DpbxypMkxYcLEGIXpWDQKwPb1AABoPDbClJgwYWIsw2ToowFENRKmCRMmTKQFk6GPBgj83IwQb8KEiSxgMvRRAZ6jU1b/NhMmTJjQgcnQRwMElcthnrCbxmNovfikw9Lix4SJXMBk6KMBhH8NhztDD4VAgwH4/vDISJNiwsSYhMnQRwPMPVEOh/mEZsJEtjAZehagOWZAua5vzEF4ftPqx4SJjGAy9GygYMCMpwW+v/wmbcZMTB26CRMmcgCToWcDhVWK64n7MPjS02AH+9OrR9Chm3aLYxaH/erKxKiAydCzgfIb5j09B196Jr16RDt0kylwGIMqF/PdmRgFMBl6NlDajfMf9eCLGcbGPsx5Ah0jHRDZswP+lf+TnzQZuolRgJwkuDhsofyI2Qwdg0wdOgeWe34yyjdFvTdfDgAoOOVLkrOH+bszMSpgSuhZgLKKj1jCkNMKtMXr0KMH9yC4dk0uSBubEFY8o5yhmzChhnhfD8LbPh9RGkwJPSsoGXpCQqfRGIjFYPfyDKznvtsBAONfX5sT6sYcxvIKZZhJj+zdCRoMwDFz3vA2bEIT3u9eiXiXd0S/X5OhZwMlA5IeRyOA05lRtfHBfhDGAqagMAviRidY/5D2c41phj68tHtvugzAYTz5j0LEu7wjTYKpcskKik1RItGh01jUcDVKnXH7JUvQ9pVTs6NtFCK04WO0LTsFoc8/Vb9hLDN0U4duYhQg25yi9YSQ9YSQEJ9fVHrtb4SQT/jUdJdlR+Yohc43TKPGGbqqznhMMzd1hLdw+kUtPSMVdejDRVEOcQi+LxNjD9mqXHoBLAHwH43rl1NKd2fZxqhFvNurqT5IR0I/bDYBmRTWPOIm89jrD0rpGKTaxKGGrCR0SmmIUtqndRnA04SQVwkhjdm0M1rR8X9f0b4YjaRRkzorGHrrv4ge2JMeUaMZwsSltA4SMYal3DFMuoDIrm3wXHch2MDQSJNiIkPkc1P0e5TSXkLICQAeAXCx8oZZs2aJv5cuXYply5YZqjgUCqG5uTlXdKYNqUwupUN6fu+unWADxpi6KxyGReV8329+AQAYejhNz1MFRrq/BNi6e+AA0N3VifbmZrG/BNpITycKAMTj8RGlN1V/iXRv2yb+3rG9GXBktgmeCW3KvssFnH94CNb2g9j91uuIz5iTEV2jDcNJVzrvJF905Y2hU0p7+b+rCSH3q92zZcuWjOpubm5GU1NTFtRlhxbJbykd0vMTG8bDMcMYjR0uJ/QUNNk+60j3l4CBjWvQD6CiogKlTU1ifwm0xTxF8ACwWCwjSm+q/hLonjF9Olr539OnTwfjcg8bbcq+ywW6CgsRAjB+fANcadY7WsaYEsNJVzrvJBu62tvbNa/lzcqFEFLM/50OwJevdkYLQhs+lh9/tiqN0sa1r6x/CDRTj9SRRgqPWDrWwudK38OhsCk6VvrdhCaytXKxEUKWA5gL4C1CyMmEkB/xl58jhKwG8CcAt2dJ56hH/9//KDuO9/WKvyml6H/mCUTbW5TFOBj8kOL9fWhbdgoGnv9TxnSOKCSZmVSjE44xT1HKxqVH3P+xGGg8rl5grOBQmJxGECMZeTMrlQulNArgdMXp9/lr52ZT95iD0qpFYqMe7/Ji4Pk/IbDqbdT94d9gA36Et2wAGxiC+/glxhl6Xw8AILh6BUou+0bOSB82iFElWXWmMdYYiUxC5/60nr8YjnkLUX3P70aGpmwgvp8RpWLsg9IRE0pMT9EcgUZjimPJhqigaQiFAAC9v/45gh+uAABED+xJ+92PlaiESWB4CZ2l6sxb0/pllIKqq1zCWo5Tox5jY2WUb4Q2fgZLaTlsjZMzq2AEBRPTUzRHUNqd00iCoQsxXYQleqz9oHgt3u3FYfMh8SoXqgw7LIA/z/p6R1XCiFhXBzzXXYhYZ4f8QvwQ06GLOJSeJX103XGDvklySpgMfcwi3tvN/9CT0HmGLehWpR8/O3LLs+GG+JiUQm3QS5m4dNIbafjffhmx9oMYfOlptJy9QDw/ZjentWCGcc4NRnClaTL0LNH105sAJLv6yxi68IGoqhnihw1DT6TaA2JtKgxb2j/xOEIbPk5MmCMJ/gONKTe1JauykVKDxTytqW8yisNlHOYdJkMfsxCW4Umu/pFkhi63ihAuURw+KhfuOWNtB9BxQ8KJbOj1F7kfUoZOKbp+fCM6bx8Fm7+Cisgid/9qv/LMxO9LlgwnRSK67/5e7is1JfTsYOrQxwYie3ei8wcaDCamrXIRVQlq5mxsXJ+fMyqvaKx+bzxDj/fIpe6+3/F+ZxLdOhsKAsixBJohxPen9i5GGGklUjFc6VgdYKMDY9Zs8XCD74mHEN66QXE2YX8sOytTufCMSk3nylJ9RiFZBo/21GwpkYohSnSPlGfoZBjc6VOCjo3UeFnjUH++4cIIzoejT+QYzWC0B3yylUtYcqCtcgFlQfREdKImoY9NCYqIz5I6OBcN+LkyDld+iTICqj5pjwrkYSyMzdE1iqBlxTUMMBl6OlBjrgIU0RVjbQcR2b2dO2C1JfRUlhLEImlzrEtQKciX9gUbDHBFRoWEztGVF/VGtsgh9yWmlUtuYOrQxwhUGao2l/J++wqw/qHEC2bZ5Gw9qbzKpJPIWP/QxOdMft6ht/4rY040xDN0p2MYCEsBga7RKKHnEmNdYBg1MBn62EAGAz7e7ZVtknT96P/kN1DWOEMf6xCeReVxB198Wr4pGuR16PbRIKHzKpfRKKHnA1kKDjQaOfRs9NOBqUMfI1BlvCneHiH6OjU2BUOX6u3HugQlLunVLlL0PPgj8Si6fxdXxKIWKX74wAYDib2PUSmh55J75CaYS+sFx6Hv8fuyJycXiEXRcvYCDL72wvC1aapcxgZIptKyzgum8RQ69ENJQuehmhGHUsQ7PeJh4N03uB8jOInFB/vRdvFJGHr5HwAOAwk9B10trEb9/9PKSqmPWHcnp37LEKEt69Fy9gLEOtoAAITfXB/4x58zrjNdaIa2GAYcetwiB4js3y03OxSQqR2yniswZfXVKqp26GNTly5408a9yQH6NS1IRpChs/2K7IqjUULPx1jIps4s6em+82b0/eYXiPdnlkLB//YrAIDQ5nXcCV71IzMuyDdMCX30IN7TBe+3LkHf7x5Ivpgpc1HO2DKPSDYNx6KxrXLx/V6lT3kIdudKjOgKRfG+D30JXYhXn0UdWTKzuDCJppNkXU4AAMm4EXT5w+kUZurQRw/YILdES3YgSlhepAPKsvqeYymDc0mvUcXfQwfsYL/6BR3b//xDwdBHoYSeU69EkgsdurGy/pX/Q3jb58kkMNyeScZqCzHrlXAsMPTh3IsxPUVHDcRQtypu+uHN69OvUG3TU3JM2bi+3C1laGNU1ZIVRpGEroyoqcRIh/ztvv+HsFbXovSab2dXUVYqF2O39T70YwDA+NfXyi8I4z1TKxkl7YLKZVgl9DGqciGE1BNC1hNCQoQQq+LaLELIakLIh4SQ9FKIjyQEqwoFQ9f8WFO9vHg8eXAqglDpqRWk1xI0jG3VS1oYQYaudPVPKaGPxIcsaTK46h0MvvRM5nXlwrEo2z5g1L8/4+3zf8VnGQmVyxhl6AB6ASwB8LHKtbsBXApgGf97bEDMqqMYUGpu+wbASfp6KpcU4XMZNceiw0hSH1GViwKpmEyeba+jrftBvG2Ks7lUueSikuzoESTpzO3YFUKPYJAwnILBWA3ORSkNAQhpBC0qo5S2AAAhpDSbdoYX3LMkqVy0BliqjVI2BrBZ6O/I4a1yGU1mm+xACsuLDCd9o+i4/mIUAMApa1PdmjaG3ngRwQ/eAZCl6ihrCV2xmZlp+/x3Q8TQx6nHUWT/btBg+vtkmjSMAPKpQ09pnjFr1izx99KlS7Fs2TK125IQCoXQ3NycEVG2d/4L664tCP7fj1WvE18vCgDEIxF5G5EwClXuZ+MsmrdtU70GAPv37gVsDrgV9AssPhgMgTJDmi8iGo2KdDBtB+AGEAlH0nr+bPorl9DqIz0M+f3oHmbahf4iPZ0cAzWI7ZJxkI/+ltYt/I5Kxmk2bRc+fr/4u72tDbE06xDHWDRiiA6te9zRGBgAe3fvAjukbvmkB0d/P2wA2tvbEWtuRjQUhBtAOBJN2S+Ft35VdpxuPwrPtGvXTtBO/cQs+fom88nQpdOU6nS7ZcuWjCpubm5GU1NTRmVb+JemVT7W1QEPuJW+9B424IdysQsAjIXBjKlToBW1u3H8eBCnC52Scw6bDYI21uVyIrJ9kya9NrtDpCNiJ/ACsNttaEzj+bPpr1yiJfUtSSgsKsKkYaZd6K9YRzE8qW8XMX3aVHGM5KO/hf5ramoSf9vsdrEt6fVM6waA+vp6FKRZh9BnbChkqA+0aO1wuRAFMGniRNgnTU2LBgDoLi5CEED9uHEoaGrCjtb9AACny4WJKZ5JOT7T7Ueh/JQpU2CtrNG9N5tvsr092Y9DQD7Xs72EkAZCSD2AgTy2k1sI05BRlQug7oQkIB5PWoJJ74/39ejTo2a1eDiBCCqw2OiPD5LC6xcAOr51CYYy9KJUxWhzLMp2kGatcuH+iGpgdgTMFsfqpighxEYIWQ5gLoC3CCEnE0KEgBx3AvgngH8B+Gl2ZA4PQhs+xuC/OSuBpIQVGvpR6h9C+1fPVL0G8M4oisEpjZVOw2FlEW0chjp0EILue76P1vMWo+eBO4a3bY3+LrntPjw76Ut4cuoFGLIm4rUbmXCi+3ej77f35IzE/DD0LCbOUaZDF57lcDFbzHZTNArgdMXp9/lrmwAcn039w42uH9+YOEiyctGR0DW8HAHwkr7iBUsSSqvZu0sh3XAeyRgRIwZCEFzzHgAguHr5sDattTn4ZHch/t14GgCg316E27bypoKHyvvJJmt9jswWtQQoAbGONtBYFLaGiertmxK6CRmUDDxDiYHG40kfiExFk8oyYoxbuQRWL0dk/+6RJiMzqPR3j70Y/94XxOntn+Cyvf/Dx1WzsbewHgAQ/HTVcFOYF2QnOOTGbDGViajn2vPRcf3FOu0rGLrEyqXrzpvzu9ozGfroR8b6W5ZN+kBk6enScaAYg2boPffdDu+3LtG9hwXBjuIJWFs+AwGLIqGFot8j+3blmkQdJHf0BzVHIUaBCw++hy+1fwR7PIoVdccAAPoeHQl3i+HRoYc2rUXL2QtSB83KlhzRDt3Yd5EUuVPh+k/4eqQql9DaNQh88HZ2dOpB0n+xrg70/uaeYQsbYTJ0A2ADfkQzlDI5HbpilEteuIy5q+LQ3hX12Qpx19zr8MP5N+LeOdfgm4tvx+dlCesGZfwc742XDh9xKqqHT+eegyPL7agN9aIwFsRRvdvxWcWRSW/G/96biPd0IbTxs4yapvEY2FAIlGXR89CPMfif5wzTmDVUhJfBfz8LAIjs2KxfVjK2Wb9KmGStYjFu05uk6SnadskSRUX8X6XKhfdniPfqmxOqIdbtRXh7ZhZ5fb+5B/63/pOcqSxPMGO5GEDXnTcjsm1jZoXZFJ6iqSDj52NQRNdB0GLH3XOuRbu7Ctft/A8aAp34y5Rzcf+sq3D3509g6mBrWkwh51BIqoNWF3YOUVwzM7EROr93Bz6pmo2DBTVo9HvF870P/wRMaQVYX09yvBID6PrhDVyAOIsFiMcRWPk/DRKNrRwje3aAOJ2wjWtMfbNKnUKikVR7PtI+a1t2iuFnbz1/MRxHLUp4BhtduSbdp6Fy4SX09q9+yVi9EniuvQCIRQ0/y0haY5kSugFkzMzBSffxnq7cEDIGdeh6ePqIs7C/sA7f3/oMzmz/CLN9e3DXxj+iKBrAYzOWIUosYPkEBSMDeX9vKz0CFMDR1Ym0eLP7uJXb9pJJSaVZXwqTVB2I0T7TYKB68N58OTq+cZGhe6ma1C/GOOJUB+HmTQitT474QbMQNsIbPhEl6YxDFSdtigoqlyw2RdMO5avWB8Pz7Y55hh7ZvxvR9kxcVoYHvt8/iN5HsrDaJGNX5SK1ElHqEHcVNeCtccfh7NYPMb93BwDAUlWDkqgf1+/6D1oKavFO/SJAz8Y/z1BKWltKJ8NhZXBkRUJCrwn1ojgyhJ1F41PXN+psxrXqVJEwFUHrOm+9Bl0/uTH5vhE2W6SaVi76rC6n70ZaVy4CnqWBMc/Qvd+6BB3XfTnpfLyvB4P/eVbzReV1UySnkDD0fOhL8wmJdNnzkDzUwvMTv4CiqB+X7k+8B+LgJN+je5oxo38f/jv+ZEQMOOsMF7YXN2J2fTHstsRnQwBMGziIncUTUleQj486H8t7FToJw4eVTtVelo8oxkNPw1ig/7knxWxYAu3iV2M0wUUuNy2lfaAIesIG/PD99bdgDuTH8mvMM3Qt9DzwQ/j+9GtED+xRvy5JSDxWkM1ydiQglcqlNuQ7iidgQ8UMXHBwJVxxyaYwv9wmAC468B66nWX4uCh99++cQcLYosSCA4V1aKopTorcN23gINoKauC3OpU1aNaXDcT0ajmsUwY1pq0RVjoZ2vS0nL0g9Xdn0GxRioG//xFDb74kPykE51KxclEi1ulB6wXHGm4vJdTeCX+KhoIYfPEpMO0Hc9eeBIcsQ4+niow3VjCW7dA19KCvNZyAgmgAZ7atkV+QhMo9qncHqsI+LK+cm08K9SFhbK0F1YgxVsyoKUpiDpMHuUg++3h7dE3k6P113X59okqexlyqDNSk8MSmaHYx4QPvv6VfnslMhy5ai2Wgcgl+8kFabaUmhn8nLJvsJS5a3eQnLPQhw9CTlmjxEXD5zTeEwTpGGLua7a3PVoBPKmfh1I51cLLyzSZpqFwGFKf3bMTGksnocoxU9OVEP+8tHAcAmF5TlPQxNvo7AAAtbv2ATHnVd0sYB42Ewep5L/OI9/WADYdU6lShUytPgJGyMG75QQy0o/psym8jD56iAo8Jff4pYp0d2vfxNPTce5tk41igz2TohiDEYBGRambOM1MMWux4fdxxuHfW1fjRvBvw6IyvYHX1XMSyie89Nvh4AipS1nu1CxBjrPiC55Pk+xWD/EQft1n6SeXMvJCXEpL+3ls0Dq5YCOPLXEljqiwygIJoAAcLatWryae5qVClhGF6vnEh2i46MWXR9iu+iK47btCpNIFEasZUOnQNhm4wH2/w4/f5dtQZerh5k/qzUcUPgQ6JykUziJ5RXsDX1fWj/0PHDUtT3h78aGVyG+KEkx/We8gw9JhHHtxW9DQb1uSwHLYXN+LmY27Fn6deAI+7EhbKYkP5dPzyyMtx48LbsK7ySOOVqVi5GFleD772AlrOXpDsSTeMSApwBuCd+kU40rcXDYHO5AIKhj4u3o8JQx58UjUr+d7hgMTa42BBLSb4O8AQkkQnATDB79Vk6HldWYnL+wQDjHdx9vDxfh/YFAkbIttVHIUkdMb7fZx0LejQ2XhGoRzSNj/lGXrM04qWsxcgtIXL5xvZs0OjAJX9ESc4iYTeedt16dGgbEEyyejFb9J0ApPSkycJ/dBxLFJmJBmJ5LAAtpQegbvnXIvKUD/uWf84mgYOcOSAYH3FDDxzxJm4Z9bVOLflA1y553VYUkltBIh1dyK6fxeI1WaYjqFXngeQmWdczqCQ0HcXjUeHqxIXH1ihfr9SamEsWNy9BS82LoHPVoDS6PDapPueelz83equxoIePiGBinQ1wd+BD6vngkIlmwtlATB5MltMVrkIaL/sdFgqqlD/9Jvp1clbU8V7u0VHHOd8ftMwHtMM5RDa+Jlm9FAa4CcWo98jP0EJnraBFa/DOWs+mAKNNCmKSVPsaZGhE0R2bjXWthYMWsIElr+Kiu/eqXqNSujJBw4ZCT2JcQuzaZ5mQjV4XBV4cOaVqA724r4Nj4nMHOB0wgt6mvHw2kdxZseneHX8SXhk5hWGVDCdt1yF7ju/LX68Gin/5BgF6hmlhL66ei6sbAyLujU+LOU7tFiwqHsLWMJgfcUM2aVYVweiB/flklwZKKWcowuAQasb/fYijONXFWpCwviAF0M2N/rsRSqVKf7mkk7BlFVDR610aot5WnWkXLFWrqzEsCC0/iPuio71SdcdN6D7Z99RrzHMSbTE7lC9nnS/0I640ub6nClQ6V95Sf4Pb76Yw5V6OqaUyYWVOnRT5aIPxQsTl6AGpSKZKZgSBhhoHASPzuAklx9t/iuKo+pLXRuN4/rWt/G13a/i46rZ+O2Mr4DVyc4bO7A38VGKYzUNzpCnCY0NBcEODereI7VUYEHwUdVszOvdiYKYykYckicqYrFg4pAHpZFBbCybJrvmufocQ3rMjCFRYbQWVAMAxvs7BUKTbq8LcCuhDldFcl2i6iYPHF1kfMY2HT1fvwDemy/Xv0lghhaVBbxKOzQaQcGP9NUZUt27ofHLP5c4YfHfN3FoTAhUoY6k8okuJyv1TL1XZfQoNm1zjEOIoaurXIwy9JietGeg899oOB47Sxpx7a6XURvq1b+ZYXBu6ypcsfcNrKo5Cs8eoZ0gQ4q8LNkzRMc3LkTbV07Vv0kioe8snoBuZxmO79IJo6DsZ4sVBMCcvl3YVDZlWBcd0tVFq5tj6A0BPlaLCnOoC3IM3eOqVKlMW4fOBvzZSX6swPhy52AkxoexJEu1arTGPG0gatYysoKS79HAOBbb4cv533gRkT3bNcv2P/P7RP3gJpm+Jx4CGeznzutJxBp1Dr3xIlrOXiChKQfOR6aEbgxJ2eENDnA24MfQ6y8CVr3tBH2GPmh14Z8Tz8C83h04qXND6kZ5hvDlgyvxxbaP8N8Jpxiz5BAZQ+pb8w0j8WmkTPHjqlmwsjEc071Nu4DiHQq2z3P6dqHfXqS96ZgPSBhXq7sa9ngEVSFeBaEywVeFfbCwcVWGLkqZylSELIu2pSej7/H7sqfXYLhZQxBVRCrfkApDT5lGEZA8uzGGLrYjuTe8aV3qsvzl4IcrMPTqP2Ff+Tp3QrnHJm3Kpy6ADf7n7+o0ZYBEl/Kb2KNRQieE/IoQsooQ8qji/N8IIZ8QQlYSQi7LjkSDUEhNwmxKKUW0ZT/CGgG2fH94BH2/ux/hrZ9n3PR/JpyKoMWBq/a8noL1y2klAK7Z/QomD7TgtzOWocderF9ulGTEGXztBfF3cO2H2jdKJJp1FU2Y5dsDd1wnXLByo4hf8s/hA2BtLBs+r1GpNNbqrsa4QBcY/rNMEh4AWCiL6lBveioX/uP2v/Nq9gSnKaHrrvaEayomimpSqqEgZFJhxIiEzvLtJD2XdlkuRC0voSuYr57KZfCFv2pRIT+KxzNfJStXaaNtU5QQMh9AIaX0RAB2Qsgxilsup5SeQin9u0rx3INhMPDC3xDr4M0XJSqXjm9ejM7vX5tUhPUPJXa+9ZZTOtJPj70Yrzccj1O860UHk1SQRn6z0Ti+t+05xIgFT067UF/4HiVJkn2/f1D83X3ntxEXlrUKCBK6x1WBNnc1jhasRDSgZJSCVU9luB91gS5sK+UiGkb5TO55hYQhtLur5GaWGsyhLtiTQuVirOmYp9UolYkm0pXQdRl6simkCBUplU0Z0z/RHqXUmGASVyS4VoEAACAASURBVFElEaIeCZJH769/rq3eykTFoaxDJT9w2hjFnqKLAbzD/14OQBoMgQJ4mhDyKiGkMYs2DCOyYwv6n3oMvj/9mjvBpt6I8t5ydSLWi8oGUHzAh8jenbrtvtFwPOLEgqX708h3yTCoe+oNUc1TG+rFpfvextrKI7Gmao5msYz1pKFgbqRADRDFuoQdHMDAC38F5cOOrivnLFSO7tmuX5FCaiG2xDuZ0b8f24snggKy1GNZ6Z/1INhBEwZdzjLUBiVSqMbHWBvshsdVmTziFBt0ifPq79Pz9Qs0ydIczelKjnpjSahKrW/THIMtZy9A769/nohDZFCHHms9gMCqd+B/898K2rTLEodDvB7KQTpA5eRB4/H0+1ksrBgDo1CHXgpggP/dzx8L+B6l9DgADwB4JIs2DEOQ2phCzk7ViJVLTCLpERUduu/JR8AODSSdFxBibHinbhEWdm9JvREqBcPAWlkNS0mZeOqcttWYNNiGpyafjTCjoc/PcFPG8eJf0PvrnyGcrR2uQfQ98RD6n3pc/KjWVTRhnN9roI8InEcfJx6xocRG24yBAxiwFyZJwKkzPmUGQbXQ5SwDSxg5Q9eU0LsRsjrQb1PYSou20VT1vBEMWt34wfwbcfmJd+O9mvmya/4Vr6UvOepJyVoTEDSk9hTP4X/nFZmApSdlCwh88DZ67v8hYtIgVoRAd5nDWLSvZ8SIFWVisazVnjTPm6LZOBb1AxCUvsUARKNVSmkv/3c1IeR+rQpmzUp4AC5duhTLli0z1HAoFEJzM7d8Fz6dGEvBAPD1D8Db3IyCeBwEwJ49e1DA39Pc3AwNtwT4+vqgdNvxf7YaPTOOgku1BOfGPmRz47yW9KSBaCyG5uZmuKNRcUa1UBZf2/Mqfjrvm3i94QRceHBlUrm2lhY4AUSjEfH5teCOhMGAe347v2m1f+cOsPHsB5KyD3fs3Am4C8RjZ3cnrAD6PO2IWRzYWnoEzmrT0bXz8AeDCF31HRSu44J2hX19Yv9M7+ds+ncUN6I+mHCW2rFlM1CYYu8hTYRCIezesQMFADqcnE68JsT1YXNzMxCNqI6jOp7pe9yVKO1PeOju2LEdcBcC/kGxnLQeSlnZ+1Sr+4lpF2J/YR0m+L14fMZSHDHUJmZI6v3lXfD/8BEUqJQDILPUELB92zYZLdI2u7s60d7cDObAbrgV5fo3roXU9qV52zZYPR6kiDOJgwcOwAVulblje7O8H6D+zEp4vV6wMWh+j5FQCNTqh5rFeW9vL+yKc6najkYiMol3f/NW0B3NYp8ov0FpPco+bW1tRby5GcyBvXADiMSiKb/hTJANQ/8IwPUAXgBwOoC/CRcIIcWU0gFCyHRIGL0SW7ZklqevubkZTU1NCHz4LgS5ye52I+brQWlJMcqbmtDCz8hHNDZCSAzW1NQErVQYJUVFUFqOE/8gKnZtTDoPcHbVrzWcgKkDBzFd4kAEAJU//SXYgB+9D/9EtS2704nGpia0MRZI5/tZvr04pnsrXppwGs5o/xRFMXnL9XV16AVgs9nQ1NSk8SQcPHY7YgAmT56MNl6ibJw4Ec4U5YxA2YfTpk2FpahEPO4uKUEQQLHbhffKpiDGWBNeljooKCzEJMk7chQWIsob0zQEOlEQDWB7yUSc6k34DExtbERo7WrYpjTBMS03MV+am5sxubIMHUjYldcGudVFU1MTaDQCNS13DX+P11mOpv794vlpU6fCUlyKeL8P7fy5pqYmsKEQ2sDZ30vfp7J/W93V+Kh6DpbuX45zWlfj/xbdhn9OPAO3bX1WvGfK5MnwpPGM06dNgxAsQ/ldVFRUoLSpCWEagTJAg2W/PEn3jOnT4G/Zib4U7U0YPx5d4AwBpG1Pq66EpaJK87uUoqamBta6Bmj5PtvtNlhcbqit2crLyqAMgiH0uVbbNqsV0vVIwXO/k63Yld+gtB5lnzaMq4e7qQlhRNEJwO5wYFqG32J7e7vmtYzFNUrpegAhQsgqAHEABwkhQrDj5wghqwH8CcDtmbaRCj333pY4ELOFs4lg90hjSa6xqRRp3qR6fl3FDHjcVTi3ZZVMg2ypHQfXopNQcKqObblo35u8DLx87/8QtDrxWsMJKjRmsNwbjqQYyuWsEGo1HMKG8hlwx4KYIWFwWkiyRJCEOmBAMX3gILaXyLdkaDSCvsfvR+d3r8qIdC0IKpcOVwXs8QjKIhLVm8ZyuTLMsbUupyI6pFZwLoPL99XVc8FQFme2rUFRLIDTPZ/h08qZ6JV4paa9v2LkfqP3GLlPtlmZ6If2K435YAAACEH3z76rfdnp1swZkJl1iryMnvrVeJWjd1MUlNJvU0pPpJTeRCntoJTew58/l1J6An8tMzE8bQjJZVlZsHoqcXgIrl2jLJSAxuZavEcliBSAVxtOQmWoD8d2p8iCLoHjqEUcpYKVi8ogmxDwYnHXZrzecDz8FsVCNh1nKZmrcZ6ZutBUNIKYtz2RdSYawaayKZjl2wOrIeal2BRVMPjpA/vRUlCLIatk0Z0vyx9+PHhd5agJ9copk3yM1roG8beDjaE0MoguRxlk0LK8UBzHur2yfQMBH1XNRlP/PpRGORnztI7PwBILPpX5LuRwU1SwlTaw4Rx4/21jMVKkfZAnIcMxfZbOt5FBm7kkU6hL9BQdfZuiowuKpLACpAy9+86bNYtreYGpBRvaW1iPLWWTcVbbGljS2CRxH38a90NgVBplLzrwLgJWF94cJ8+ikrZpGiAb4EprlNyBa6PnoR/Dc815Yob61rYueF0Voh15SiilFsXx1AFuEStNJJET7z0VCPV6nRWiKkWEZKIpu/EO2aWqUB86nQqGrjkRy489V52NtovkKzOPqwItBbVY3JWQixoCnagLdMkZepoSqK5Erwg9q4feX90F/9svp25PSp/SwSpHHtA0FtPuh0zayGIDVIgOKamMJ2MUS+ijCnwHKQeHVvS3JKRh/vZqw4lwxsM4Qy2mtwT2mfPkJwQVgjA7a4yxyUNtmN/TjNcaTpRbvGQgjXJ2v2kXywiy+M8APrdwFilGGXr0oHq6QAFHDHGa1z1FCak4q/gaeojHQcFJ6DILFyhizig+zKpQH7qUDF2AkqkYkFS3lXC293P7EuazBMCi7q3YUjoZAYvDcF1yWvTGEs98crn6kbQX3q5QY+aKoeuNhQzayGai6frBN9QvjGZP0dEEoiGh+5582FB5o/bMvfZifFg9F6d5PlMNMiX92F2LTpJf423diUTfr4XzWz7AgL0Qq6slk0ImNtesROWSr8iTGmZum8qmoCw8IEYpTAUhjrcIBb3F0QAqQ33YUzQu0XQuk/tKQOMx+OyFCFkcqA1pe0JSRQyTqpAP3c5SWcA1zfdsgGE0l0xCUdSPcQF5qIWjencgxlhFhp+uNOm9+Qrti8LkkEsbf8mzdt8lj8gYePcNQ1X4/vhL/RtiMW3hJSMVeg4lIeU3kqew3ocMQ09KOcVDSweeBIOD981xxyJOGJzTulr1umxWT/J85KVtRntTtOCMcwEAs3x7MN7fgTfGHZ+IA5HJElBaJm8al+TnYEGwuXQK5vTtQtHZGUZFVJmAJg+2Ym9hQkLXzEKTLeJx0WRRKaHLblPEAakK9yHK2NBvlxoRquvQqcZ5KZpLJmJG//6kVzdt4CCsbAzbSo8AAAy984r2s6jRrfNdRHY3Y+ClZ/ITH0YFvb+6y1gdKVYMNK5nJ54/HfrAC39FZHcKpzmxzvwKV2OOoUcP7oN1w0eIdSsGJM88lct+oxC8GvUQYmx4q/5YLOzeasiRKCluuXAsBApSG6BCmFAAZ7V+iH1F47C9uFH7/lQwGDsjK6jUf7CgBgP2Qsz27VaN2mcIqgy9DR53ZWLDOJJHhs6bLCbp0AFUP/JXVP7sUTGAmHg+xFu6SDdGWYnbuxSSY7Xl/YDNDY+7KslCqOCM8+Bgo5g60IKtPENnc5jIJLx5Hfr/8igCH2okIskEwxGHSE+HnhGM1dX/1OPwfltnxQMgPtCPzjtuQLSdN2Y0N0U5BD9bBedzv4PnqrPkF7Kc8IxIeu/WHYMhmxvnt2hnCdfTrwrHovWG2qpAwiBO9q5HQTSA1wUTRpUIdEAKtYPsQ8pfsg+l/nIzH0hrTp82Qw9/8aK02zliiLMA38urXdTMUmNdxmLqSEEplacYi8XgdVWA8EG3lHDMmA3XguPhPuVMFF+SiBNUJTB0mR5danEkazTxW0UaFjZ/Jw/K0ysSJ2flM7N/L/YUjUPQ4sjLEt6oKsQQhiH0c/Cj93LbfpplrOMnaV4Lf/4Jwhs/S6iATQmdA9HKPJLleJHarqshDoJXG07E9P79mKFwJNKE8iMTXiI/O6tt4khN9ZxsFEs6PsPHVbPQay9S1cVGdjWj9fzFCH6moQIapoBe/ZJ0bQCnP68LdKEy3K8andA+bSZoeZVunVKrnIIzLwSQYG57NBh6YPVyeK4+B6ENHyMd9D74I7SetyhxIh6D11mOinA/bFRb9UAsFriXnCMeCwxdZukiWrkoCssYevIAPlBQBwBo9CtchvgJ8kjfXrDEwq3ghjnVYroYrlj+rF8j6cowWLkoV2tSMII3szBxmwydh1anZavvS6Fy+bRyJryuCpzHS+fW+vGp60x6aYLKhX8GNQldMWGd0f4pWGLBypqjJc9IQaMRsMEAws1cWOCQIoxtUuaWLNFy9gL0Pnq3+kVKEdq4VjyMEQZbS45IWLeo0FDwhfNSD2rJdftULrF2SdSPylCfqEdXZpqKbOfM+yJ75R6NqRD44G3ZMY3HRBt0AUyZSmhcyFdl7ngYBdGAXELXdCySqFxUxu/+wjqUhQdQosilKgg10wcOglAWu4onaAs6owXDlZxFQ4AJbVqrel4X6ZKs84xJwpupcuGhMXCzlUT1yrMgeLFxCWqD3VjI58N0LT4FltpxmmUA6KhcdD4+haQ1LtiFJt8+vFt3DCItiZVBx42Xou3ik1IzRdkmbXZSgZa9sVL62l00HiGrA3P6dvHXNfo2JUOX/JQETztiqE2U0P1vvCSnJVc2mvE4Op1lqOb15+Xf+xlqf6sRCVrxzqrCPhlD15xcdWyzAWB/YT0mDqm4efMCgSseRkOgE7uKx2eVM1OaOzRvGKZY/tqhh5PHGqUUMa+2G72huO3SSVln7CVpAEwJnYPmsibbAaMj4X9SORP7isZh2f7lsAgvTWuJK3lRSaoG4ZrO8ljt+ZZ0fIp2dxXWrRW8UglirQq1j6a51nB8SPLGN5VNAaEsZvl4u3JVdRYxMKgl1yXhjScNtqPDVYEQIw+n1nHjpTmzIhjcthl99mJxk9PRNBcWDQldKW1Vh3rR5ZC4/0vTr0kh856Uj78osaDVXY2JSnUL5Gq5qQMt2FU0HjQLlUv7padnXNYwRjh9oirfoBSea87TLGNITSS5J3Zgr+ZtwVXvyE+YZos8NFUu2TGuJAbJIw6C5yd9AeP8XpzolaSXYxhVJlrwpS9L7tFgKjovkykpTzp3XOcmuGIhrKhbyJ8xIHVLXY2H2cplc9kUTBpqR1GMy/SuuuFMiMy5wtowEaU3/CDpHvGn5L03+j2ghEGLIiVddN8utaJpw9K8EftXLAcljMjQdT9AxXuuCnESutgrGhK6jGEosgO1u6sQY6xoHFIJuSWRxqcOtmDAXggvSRXv0DjihMH+glr4lGGAs8HI8nONsL+peIYRhp4h3zEldA5aEnq+Nv/eqV+EloJaXLr/7YR0DvWUVkMPP4Pii66UnFG8NAMZyK3VdXDOl7v8O9koTujciDVVcxKegQpoSROyZeAwSElhxoYdxY0y71Athi6VbMtuugNF53D26q5jT+ELSuiVMLGJPJPbX1iXXG8OnpHp8YqbmqKFi96qSiGhV4X6ELQ6EzFntLw4JcxAqZZqKagBAPUsWJJvYAofDmEHzU0I4V1FDbhx4fdxyzG34Nrjfownp34ZUZID/fxIp09U4w+pvGsNSeiZkWN6igpQySwEIC9Bmny2Ajw36UuY3bcLx3YpgnARJrUVoJIJGPESI4Bt4pSk00s8nyJssXOeo1RRQKMeWZuAMZ1gNJoyS1NyocTP5pKJiDFWzJYx9NQ2/lJU/OBe1D/9pswcUDqRV4f64IyHUySNzuKDoVTUgQsSOtEad4Cq+z8gNV3UskOX/FaM33ZXFQhlURtMti+XCgSNfg9sbBQ7WK1o6MbR4SzH3XO+DgC4qfmfOKttDd4adywemvlVxLO2Cx5xEV3lVApDihQ0x7o6spDQTZULgDzq0JXVAfjrlPMQtthx3a6Xk4ezER2YchYWBpDeBhZR1y1PHWzBeH8HVtQlJytIUKx2OqFy8S9/VTWanxR9Tz4E702XIdaZji13ou3NZVNgZWNo6t+XuKph409CiXjvUhNFYrPDUlEF57yFiZslm6IMKMb7O0SzvpyDUnid5bCwcZQLYXN1VS7Jm6IA0M3r0amR4FxKhu6uREW4Hw5WxcdA0p6VsjhisA07YspUFOmBAni06RIAwJ0b/4hTvetw7e5XcO2u/2Jt5ZH418Ts9OzDZbaoCTWBLxVNKa5HdmzJfKIyJXQeWswwxxL6ypqjsarmKFx8YIU8QbBIh4EXonhp1IDKBZSqMg8C4NSOddhV3Ig2e8KCIskbVayH/yvpF//bL8P351/rkhzZzq1E0or9LBnUm0snY+pAC5ysRCpXk9AJARnSsBmW3ubgdcMKiWbiUAf2F9Zq5+/M8IOhfLzuLmcZKsO+RDRN3VVVssoFUDoXAfpmi/Lx63FVol4Rv0WEYrUwZbAFu+IuxLOQ+j6snosdJRNx5Z7XZV7QZ7etwSkda/Fi42nYK4lymTZGmKGrqmRT8YwUNBOL1VA6PfXCJkMHMDw69OaSiXhy2oWY1bcHFx14V50OwmgzU8k9MogSuk63x+OaL/tE7wYQymJlqUqmE7040JJrrE87LknGEFYAFif2FjVw7v7Syxo6dBKQMHSN91r+7Z/AUlENxiWXQBv9HgzZCtBn19AdZ/rBsCzAUs5kUcLY9E1NlUHE/LDHIwnnIi0JXaYOk+jTAbS7KlGnom7haJGPn6kDLQjDgoO83j1dUAD/alyCCUMenNqRbK99ze5XURzx4/HpSzOfNEZah672faScZFJct1hT36MFk6HzyJOVi4Dm4kbcO+tqVIZ9uHXbM7KNUBmMqFwYdQldjNdiT97gpPGYJvOoiAxgTt9uvF/alKAqlRm6UoLIUnenuXSmFNtKJ4ElDGb3ycPgWiqqk24nIIiefBYc8xaiaNnXYJ8xW7Va98lfRP3Tb4DY5CaKjeLGqEKPnq0kyEvoXmd5wsIlBZQTN4EQRrdUTpOeNCcZv4M2N/w2tyx3qgyKsTdtgEukvLNogiF6ldhYNhUtBbU4v+UD1fFeGAvi67tfxr6icXirfnFGbUjfi9/qxAuNp+Omhbfi4pPvw1ePvwv3z7wSG8qm5c8Yhu9fS3VCTRfv0VgBCUhFjMWS+XjLkxHHmGPomptT2WbjBvBW/WL8bO43UBwN4Keb/oTiqFo2UR6G7Ej1rVxqn3xJWYCLy2LRrvtk73p47aXYXjJRUVBDP6vslywlg9Bnq5Pr4NveXDoF9ngU0xShEcq+dTscs+WZ6kEIaFEJqu/5HUqv+lbK1Y5yIhKsPzT16Jk+J2URYQGfo1hUnQBIy2wR4J2LhABdGp6isslR8oG3u7iQCJOapqo2p/wGakK9KI4MYWdxZgz9tYYTUBoZxAmdn2vec2zXZszt3Ym/T/piZuaM/KNuK5mIbx/zPTw/6QuoDPlwwcH3cVzXZuwsacTdc7+OX8y+Bl6tePJZgOWdpwokYRpCmz7TJzmlysWSMd8hBoIBZoKsGToh5FeEkFWEkEcV52cRQlYTQj4khMzJth0RmhJ6ZjNllFjwceUs3Db/Jjw57UI09e/DvRseTymd6erBBSRZucg3Ra3VKlYasaiuFH3y+AI42ChW1vAMUmBcmhoXxQUjun8dsIP9mtc2l03GjP59SbFPGHcBXCecIb85XTIUDLowFkRFyJds6ZKlhE5ZFt0slx9eGAOl198Kxq1jRaLyvgRbdBlNSVYuEjWLhKF7XFxikIn1lSj/7l3J7THJK4Kpgy3YlQFD99kK8Hn5NJzm+Uw/Zg2Ar+96GRHGhqcnn6V5nxYoZbGxbArumvsNOOIRPLDuN7hz059wxb7/4YadL+EPH92Dr+1+FdtLJuJ7C76Djypnpd1GUpvgYuqsrWjCezXzsaJ2AT6KFGFfQR1nipnK+irVWLJak787q441lLRqg/eli6xqJYTMB1BIKT2REPJ7QsgxlFJh2rsbwKUAWAC/A3B+dqRyCFOCflsBKAhYQkDBO6jYCsA6LWD5Y4rENeFcyGLHgK0AA7YCdLgqsK9wHLaXNCJgdaEm2IMbt/8Tp3asM8ZrDKgulFInFZxHdCRwGouJkwWxO2TBp9ynnAl2aADHBg9gTdUcXLv7lZSJOfqf+wOIJAmDWqCstKAi+VJK0c84caCwHpfvfVO9WJYTiZqE3Oj3qNiiZ69y6WI59Y4Qx8W18ET9MqoMvQ8D9kKEGRsGX/47Km69Wzd8rkxCd1eCoXHUOoyr+6YOHMT68unwW5woiOtbMknxUdUcsMSCEzs3pLx3XLAL57e8j5cal+B0z6c40kDibwHb/QQPzLoK9YEu/OLzJ1DIO50JsFIW57auwjHdW/HLIy/DQ7OuxBfb1uDqPa/DwaYnzXY6SrGibiFW1C1Ar0ORsLsHwDHfhZWNYeJBFkdMuxBTBw5i6kALGgKdYKTjR0eKDjM2HPBTdO7rRXP9seh1FCPM2BG1OxFjKWxsDHY2Cns8CgcbhSsehisWhisegiseRnlRshoyF8h2mlgMQPBpXQ7gWAACQy+jlLYAACGkVKVsRnipJYrfH39n1vVY2DgaAl4c17kJi7q3Yl7fzrTyg2ZktkgNWLnEYyKDKDx3GQZfekbSJgEYBqcO7cbKgilYVz4DZ6Zg6NFd2+TLsKw3Y5LLs4P9+NwXB8Yh4e6fVEwjDILhZpPvb/R3YFPZVMQIk5SEOt7lRazTA2t1mqaNLItOyu1tiCqXNGLOCKgKC5YupXC89yYqblUJbCbl11SucqkJ9sLKENW61fZYpg0cBCUMdhc3YK7RHK4APqg5ChOGPGj0e1PfDC7f7fs18/HHqV/Gw+seNfTNtLqr8JO9LhRHfPjppj8nMXMpakO9uGfD7/HcEV/CK+NPxtbSyfhO8z9whFpMGwlihMHaiiPxTt1CfF4+DQBwVO9OLD3wLhqHPCiNDIKAgr3wGuxZvhz7CutxYOoirK6ei7f5fQFXLIQpgy2YOORBSXQI7lgILGEQYWzotRejx1GCbmcJOp3l6LcXAR/0A+gHpn0ZDI3DEY/CTmOwsHHEiAURiw0RxgZWZcK/2hvCcSl7Ln1ky9BLAQgBDPoBSNOQS59C9YuYNSuxrFq6dCmWLVuWssH6aA+u2/kWJ39TCgJwvxkGTDyaOEcpGP4eBtygc8SjKI4OoSTqR1l4QHeJmQrezk7YIhHZQ4ZCITQ3N4vHlvZ2uKRlPB44APT6fOjg71NqI73tHoAQOAD0dPfALrnWPzAAEvBjel8PykoG8H7tfMzn6/T5+uCVtO1W0JaoYxCdkvuUcIVCsADYt3cv2HCifwQ62z3tcED+QtsffwCbS2fwH0RycKTm5mZYvV5IndPb29sRKquX9ZceSEcblEqPxiEPYowVbe4qkSH19vbCDmDw389g8N/PYOjhZ5LqUoPwfDu2b0dnzAorG0NZmDPd3L13H2ivtqoJlMI5cz6sWxOJgYXJoNtRioZAF5qbm8G0H4Rgq9Pc3AziTTzT3t27xWseN2fh0uFlAYcTSqf+dkVfApzKBQB2Fk8wzNA7HaXYUTIRl2msqtTgZKO4dvcreGDWVXhj3HE4VyNzl4BuRwl+PufroCyLn276U8KuXwc2GsfVe17HvN6deGzGMtx29E04zbMWF7S8L9sopgD2Fo7Dh9VzsLJ2AXz2IlSEfFh6YAWWeD4T/QGkCA+0o65zI07o3IjIZAesq19Gu7sSu4omYFfxeOwqnoC36o9FxCLfhHfEI6gM+1AR8mFh9zZUhXpRevKpqKgqQ+Njd6A0MgQGFGxRCRiJWpICiBELglYHAhYnglYHBifNhKtkqeGxnw6yZej9AAS7sWIA0h6U+cGpFd6yZYvaaV0cEQ+h888fJV+w2YF8pSNTQW19PQYdDkjdPpxOJ5qaEiaFge42CEaCxOlCTVUlfADKKypRxt/Xoqi3qrwMxOFEP4DyslIMMYy4HC8tK0PcwiDuH8CJnRvwxrjj4SjnbORLS0pQLmm73W6H2nRVUlqCiiYVs0ceHU4XogAmTZoE++TpADiVisCm68c1QJnuwW2zYUvxZBzZv09VYmtqasJQy05IdyXqxzXgoKK/9BAtdEHp6iTECT9QUCcy9LLSMkiDzRqtX3gP06dORc97e1AV8okWH1OnT4elvFK/ggf/gM47bkB4I7dArQpxn4KgR29qakLEYYEgBzc1NSHqtovPNGniRHjBfTQeVyVm9e1BXf1MEIczqb/HjR8PpfFpQSyEcX4vdqVh6SLkqz2hc6PhMgCwsHsr5vdsx/MTv4DjOzdpMulBqxs/n/N1+K0uPHpEAJXvpZdVaV7fLvzqs1/ihYmn4636Y7G8fhHqA12oDfYgyljQUlQPn7UAFjaO+b3bcUb7Jziqd4e2ZRqAmro6kUlVlJVhEBQNgS40BLpwqjcRijnM2OC3umChcdhpDM5YOEkqrZ53Jax1DWiPJMxvbU4n4hKGTsBNULZoQDSycDono7vUbXhsKtHerr1ayXZT9CMAS/jfpwOQZhXoJYQ0EELqAaThpaIPzcxCGTBz2yR1KwJDMKSLTgysmt/+XeJYpGPTLNGhxChByQAAIABJREFUg1L5cp8wIAwDyrI4uWM9YowVKwc5WY1GIwhv07ZSSI9uBSS6XqLiyerZux9t7mqZu39Ss9mqXFR08OMCXbCysRQhANIDpZzKpUqapUjP7V8D5ZEBMDSuERddOEzWoffaixG22FEf7NJ+VxrjZ9pgC3YWTzC8i7C6Zh6m9R8wlE5RCgLg2t0vI8ZY8MS0C1XDAgxZXbhr7nXwusrxwy1/w1RnZsm8i2JBXLv7VTzx8f24avdraPB70W8vRJixY77/AG7Y8SL+suZu/HDLU1jQu12XmQOQq0p1Mn052CjKI1wseheNqasYWDbpnRKrTe3OYUNWDJ1Suh5AiBCyCkAcwEFCyI/4y3cC+CeAfwH4aVZUStvMpblPNiEsJWWrH/4Lap98Mfke/mW7TjgdtvrxqrFcym+5S14kHktcj8flHzVhABDEWvZhkt+DxiEP/rePkw4C772Jzu9/HdGD+6AHvc3JgX8/i+g+lTguKWKqbymbDABJDkUyZLspqsLcrJTFuECXbGM0/Pkn2bVDKbqpQ2blpJeJRkGk+MtCWZSHB8QwurFOD+J9EnWBhsVLu5szWawLdmuPT43z0wYOcpEXnckRO5VocVdjf2E9TtQxVVRD/T+W8/T14Ko9r2Nt5ZF4YvpFiEnej8dVgZ/Mux4tBTX4wZanMMu3N3OPSh7lkQGc3/oBbt/6NB5a9xvcv+FxfMf7Ls7wfIqimI55sRKSvjPMSzSFIJpsT67D0O3Ts7fcSYWsbWcopd9WnLqHP78JwPHZ1q+ErWFi7irLJsuLZGAwxSWwjWsEtHRiPBMsPHspoi37Ubzsa+KlgiXnoPeXd4nHzqMWIdrCM2VKQSwMqCBIEMiU1yd71+PpyWej3VUp6hYFe1ttkyttxtqvGRZAn6FvLp2CwqgfjUM68V+yzKijZac+wd+BZolNfqxdqcRKD6FIDP2wy81WMzQxqw71iRK652vnyq7FWvbLmIFgASWYLNYHtBm61gQj2P9vL5mYUupeVT0PDGVxXFd66hapDfxZbWvQbyvEvyaejt1F43Fc10b0OErxfs182NgY7tj8V8zjk5z4fv9AWu0YQgYCmbTvdHPxSsswjLrcr3JSV0IX2s5jGIQx51hkrWvIWV2GbMl1y6aRLQicPXbFLXeBKSzSLOKce4yonqCsQkJnGNmxEArgg5qjJE2msDrI5JmV0R1lWVq4hBYzfXuhMez5YkRxmK7KRZ3uRr8H3c4y+K25iQfuGeTMRGVu/4YldDmqQn3y3KJSMARUYk4q+Ci0uythj0dRIeRiVesnjcmxcagDhVE/NpdO1qWLglO3zPLtQVlkSPdeUqDYtlf0xaX738ZtW54GBfCPSV/C+zXzsbhrMx5e+2uRmecNOua/mpD2nUpOX/Uy6u1QqqJysWkz9OFIE5gf6/axAsWLsjfNQaR5k7GyhEEqm2dxWZ2utkFQT7BUxsCVH3hFZACz+3bj/Zr5+Mr+d7hmUi1tMzFblE4SCtVJu6sK3c4yXHTwPd0qkiaxtM0W1T+qCfyq4GBBLZrSsIvWQnu/wNClXqIZMvSwD72OYsQJk7xZzFJ03npN4pgfKx5XJWqD3dzkmGaCFAYUs3x7sLlsCii0h93uovHocFXiwgP67wzgxpx0RKlNbou7t2Bx9xaEGSusbDy1HjtHyIhBSsadUQldU+VCabIApbeaE9+bKaHnBcoBUXCGdjqqJKQh6ZJ0ObowgCgr/6gJSXI1PsW7Hl5Xhah2CKzm3AJoQJ5YWFZHupBJIfLygs3vvN4dulU4jz4O5d+9C07BSSdnEroQAiA3G6Ot+zmVTbVRt38dVIX6wBILelUCiCljcQvHsqBcRH0VqLdimNO3G93OMlF1o4aVNfNhY6M4tnuz5j0ilJOKROVS9YvH4Zh9tHjsYGOwuoyF8S256luG7tNFBmOZZCSha0U0pckSut5qTlS5GGs2ExzWDF35otSCZWnBkMpFvDlD5kVZuXUIw4BG5NY8i7s2wxkL471aLk66/42XED24Vzv8bQZWLlS5KSp5ng3l01AX6EJNqlAJhKDg9HNAbHbd+zTLazDVyrAP7lgwZwx9/8efwsZGUSoxRTOsHlLcJtii9044MvnepMiLFHEQeF0VCVtrQiB8/UyJRHUjMCVCUHbTj2TVCJmiNpclJ0kBuFAXH1bPxTHd21AQS+1Ryg7I7e+l78F51CI4jpwru86UcJvAxJUi4Ybk+S1VmUWJzMhiKxMduqaEjmSG7nCi+PLr1evJUHWXDg45hl55p368bxmUUfLsaTAbI4wpw9jcCR26XEK3lJbL9a7gHD2O69qENVVzEOaTJvc/+6R23ZlYm0gGbWjDx+JxlFiwtXQy5vWmk+FIqCt7T1Ghlgn+DhzMUbKLroJKVIV8+vsBBlHJO7b01kxKvqjMlsNyJo4xxoq6AMfQ5cxzMWy8X4BUd6xcVdYFu1ET7MFnFYlJRBrJckP5dAzYC3GKdz3yAcbp4uk6N8WdEmQYjoJkokOXjiOjVi6aEnqyDh0ASi67TqMek6GnDdfCE+A+PY3BJAGxGZfQlfG5VUHTY151T70OALDP5Bw+3Cd9QRzsTHklii66UhbbRcCpHesQtDrxMR/QiPXrJY7IjqH733hJPG4umYSwxY55fWmmrMuEDB21R+NQBw6oJbvIAB2WQtmGaDYQnIu8VGXDVrnXQSna3byFS5AP66qlcpH6KagE6jq2azM2lU0Rc5o6j0kYm71XezSKI0MpVWQZQ7DySEdNlWl8oQwYZCZWLlrt0GhEhaFrD2xiWrlkiAw7jNiM7xETp4Sh6yWXMIC6p95A3VNvwFrJLT1t9eMx/vW1cC04XhzshWecB2K1qjL0pv59qAn24L1aTp/J9uuoP7LWoSfwefk0WNmYdvwWvbpyEMtFwAR/BwJWF3ocJenVqYJOS4E8bG4WcLBRFEeG0EmSGbpyM42ycVHvLerQZZKhRF8rYTBq6qBjuzYjxljFCR6UK+N1luGzypk43fNpUuybXIHwm4Ipg8BJndVSrBptUzQ8KiWTRsFZF4m/Xcedql2Z1A7dqA5dY+x13/UdRFsPGLpX3rbJ0NOCpaIqw4LpMHRnailT/AD1b7RWVsNaqRF9TSjLaM/uDChO9q7H5rIp6HaUIN6nI2FmZLaoPgA3lE/DjP4DcMXT8NJNc9UiQGAQxOlCzW//LrsmDQGQDYIWOwYsLsOJLTQhYbhVoT50yiL68FA4pFD/ENpdVXDHgigVTAmVZovCXJhCFztlsAUThjx4veF4CeugeHn8ySCU4sy2NbL7rQ2NBh7KGAQ77JTms7JC+mPBMeso1fPSScMxM3FPxR0PalcmlbbTsEPXQvedNxuqI6ntPOGQZOgll12HoouuTH2jUoeehgOJoCvUhci7cpBuSmDsGsz1lI51oITBypr5+vlAcySh99qLcKCwHvP60lu6M0WcFC3mCjVckH9XhIH9iGmouvsxcYNrAh/HJduN0Q5nBQCgNphhmj6hbyUMtyrch06SmqH3PPgjLihXoDsx1cWi6u87xaRMAJzbugoHCuuxrnwGAKDFWYm36xfhdM8nqJDEXql/5n8o+vIVqZ5MhvLv/wI1j/0DgIrHqyAUpYgCSqVSqvA8SvWRqwC1T76I0qtvUq9EIihJma7eJrb0mlRCL1p6NezTZqoVSU8lpNc2Py7ymTD7kGToxGaH+5QvpV8wLQndBfeJXwAAMMUa0YFz8eIEKYIfsFqST22oF0f69mJl7QJ96cgwQ+dop/EY4n3JDG5DObdBl96GKFB63S0ovf5WOI8+Nq1ySoHeOX+xGDBLTHahTEcnAesfSvkhdbiyZOgCqVKGHvKhi3ElL7JVUpC1uatk0QST90KMr25O8m7AuEAnnph+IVYH3Xh45hVwx0K4dP/bCmL16yq+5OtJ5wpO+RLsGnGQBLVlqjj9ilIAAOcCuWM5jUU5z3AtQUsq8RpdeVrUJXRH0xwxVWLFHQ8mzGvTqRtA2Q23pW7bZOgZIBMbVavV0MsrvOAyMCVlKL78Gxj3wkpYNBg6zdSiQwLBlZ/6/UKlmji1Yy3a3VXYWagdcS/dBBe+PzyCjm9enHT+s4ojURHyYVKKONVKMO4CFJ13SfqeoiLUO6DR36Epocf7+9C27BQMPP9nDL39MsJbNogMJ7B6OYbefAkgJGcMXWarHepDmFgxYJOb8Skn3TBjRbejFPWBLtE7kx0cVKhcEvsP7iVno+qe32mSYKNx3LLtOcSIFXf11KLHUYLvb302Oa0io+GNKiDdsAcGJfQkGqCiSuItgTTHivRbNTiurTX14m9ZLBeLRfQFIAwjX60brLv8ez/Xjb+f2Mw2VF1GOGwYev2zb6ncozi2WA1J6WXX3QJCCBeDXekaLUWmG4AqYAND8jpVcFzXZjjiEbxbO1/znnRpCXz4btK5MGPFxvJpOKZnWxZTVbpQUTlJfk7we9DmrpYFiRLA+rg9hcD7b6Hv0bvR+YPr4PvDwwCAnvtuR99j94kMvTgylFbGH1VKJYxJ0Md3KgJmsT65nr7DVQFKGNQHu2Ep4e5lgxLnMIUvQMUtP4Nz3kJdOiYNefDYpw/h3goPHvvkIcxW27xOpb5JdV0Z0kGwcmFT6Kdlj0OEH5r3qNMm2RxW0qll5uouRNW9T3DVSyR0wkgSPjOMjA8YNvVNJQwy+XfMH9sMPY1QlZayiqRzjML5gVitaenRjSIXKnQROuoUVzyMxV2b8WHVHIQYjb4xSAzb70P3fberepxuKpuKsMWOhd1bU9Yj6MyzBlFj6Infjf4OxBirmGBZBuGdSnSmgVXLlQ2gw1WBmlDm0jlR0aELFivtCs/Nnvtvlx17eLrrA10oPP9SFJxxLoouvEJdQk8DhbEgFjoDmPy9H4MpTf4GUq7Y0txEF76ftFQuAmNO2jRMI4yFUYZOSEL3Lt0UleQdAGHkqwWDEnrK2FCZ2M2niTHJ0INfvQmu45eg4vsqab10ICxRLRXVqLzz17COnyi7TizW3O5EZ2jRoV6Xsk51nOH5BAGrCx9Wz9W9T6w2EkbL2Qvgf0+etWbghb8guHp5kiMTAHxaORPuWBAzfVyyKqZUPVxr/fMrUPfXVw3RkQrE6YLrpDNQeZfUcSzRF2JMFxU9uhAhUOZlq5wYGY6h12WrbpG0B3DqG0IpPG79BBkJG/RuMIVFKP/OnbAoJkNRhadgVu6Tv5iSJveJp6Pm108lX0ilcknBpIR9ifgR01H5s0cTQlYqhi5bcfBN6SXiBuBcoEjaloHKBYQkHPfiGgydYUCkRg85ktATDN/UocsQn7sQlXc8kDbzFdzOicMJ18ITkm+w5EdCz4mILjCgFGOhqX8/xvs78L96jU1HxYQQ59UR/U89Jr9N44OMg2BtRROO6t0hpvDTcue3FJUYc8AyAEIIKn9wH5xzFiRolDxLQ6ATDI0n6dHbv3auKDHTiDS6obwfooTTYSv15w3/+TB9YiXSnY3GUR3q1Y2tAnBBzsrCA3DFw9pjMGM+wBVUDWYlYUKqDnkSRll6ww80W4hNmw3XguPh4p2YbBPTSB4jtGGxwFJdlwgbIHm/9c++hcofPywrJtNXJzFTjW9OIqHHvYn9H8JYxH2NJFWqUT6T6j5hojc3RTWgxyglnVZ6/a2y+7VC8OZa5SKY5pEUkocRJKQJ/cFAAHyx/WPsKR6P3UUGQg3zgzCJgWsw9J0ljei3F2Fh97bEyWGIUZEKNhrHuEBXki16vNMj/pauNpQWL53WIlDCoEbC0Mu+dXta8X1EKNz664LdYuIKLbS7K0UPUaK1j8PTnPGGsko5wjCiOlK6YSi9LqDonKUpm3CfeAbGvbAyEabACFkOro+J3YG6v7yC6vs4Hbcs3ktZRZLg4Jy/WJVOXTCMOuNlmMQkTwgYd4KhG+3vlPsNph26PtId2PZpM1F63S2o+MG96jcY3BQ1CveJZ6DkmptR8tUbsq9M0PcZmN1P7lgPRzyCt+sWpbxXHIQKBq7lFr2qeh7s8QgW9CQYeuYWK7kFF9NFxdKF7zIalnjZKlQugoWLVOVSeFaydY8RKNVUdYFueFwVmlMxBdDqruaSWgDqYzBF1ihDUGM4hIFzwfGo/OkvUfyVr6HijgdSl0nVTEGhgXGauF5w2lkovvQ6lFzxTbmOWwXORScl2pHuiyXp0DUqUEmhCIATSoSJmGHkceCN9ndKHfooNlskhBQRQl4lhHxICEny4iGE7CCErOT/qYSbywHSmZXBSeBFF1ymaZlCGCanEdGIxYLii6405oSUAgKDNZLKqyAewoneDVhVcxT8lhQOPPzgSgoXoAweBSBOGKypmoMFPc1y79ARYuhKE7HGoQ50usoRtCilapU+U5wSGHq9PR37aSW4fqAhOUOvD3YjYHWh36Y+7vrsxRiyFYger9oqF413n4pBSK03lGAYEELgWnQSvzIgSdczQzJN7iXnqN9qsaLkiusTOnQdfXilxAtUJrEbzVnLW6clnWYYWT8x2TJ0tdXQaGboAK4D8Pz/t3fmYXJU5cL/nV5mevYls89kMpNJIBOGECBA2EK4CZsi4ieJiFdFUfS64fL5XMQNxYt47/XDC254UbjyIV5w/YQrn7JoEjZB8Eqkk0DWmUz2ZJLMvp37R1d3V1dXdVd113TXdM7veeaZ7lOn6rx1qurtU+95z/sCK4APCCGMhtQDUsqV2t9rybu7gV1FYrOeP9FdaY7BZpdXbJpcolzW/zxj/iKeal5m2JK4vzEud6zcxOTy1+oFHCsqd5yHcqao++I3E763D0WTXRjCsZqZjwwj9D0lcwhNjtG6KoMFacZDG+LtNGumlH6LidEd2kRuLIWf7h4sWhiJY1JywSqcLCwyw9SjJY13SMZZvUyUVsnyi2KfixctsWwz1Y+ICARovOtB5nzhXxOzA9mcuBTCZ77ATPgoPu0sAAJ1jckrmW2YS9KaVLRorraDgmVANgp9OfB7KeUU8N/AIsP2WiHEOiHEPUKYRCdyg1QX0UHwn/jx/AmjI0tbZh6QDkwuAF2Du1k8sI3ftF2Y6Jtt3N3qeCaBi9Y3LqV0coQzDm0ybMnPCN1v8K6xiukiTVZlGs97d2kDrcP7qbzm3dnLNScxLk/UlLK71DxeTzT0b3t0hK6fVG2dR9ujL1J6/qr4Dpm+EZm4zSWZy5K+2/ceScBUZ+rCAZ+xPLYa05gAJp0Jr6jrZErPXZkwv5H8Y2U9QpcjIyblUHntDTTf9yiB5jaDcha0PmR0czU/dsrNRRE16GqiewPZaKxqIBoU4qj2Xc8FUsrDQohbgBuBu4wH6OmJZ8Fes2YNa9eutdXw6Ogo4XAYf19fQtij6fomwlqiZl//TqL+FXv27mPSJIFz0cGD6F8rNm/ZQsnEBNFL2du/2yysUqwNK7ncJPriNzg6xsFwmNLxcdu/wlf3/pHbT30fz9SfxkX7XwHg0MGD7NHJKA7tx2zKdnx0NKGdoUCIZ+uXcNG+l2PeLbG6E+YypesLN/pLb8SoHx0gNDnGToPr4tY3Xk86x+mpKcLhcGz/vtIGega2snnnrliZU9lCQ4MEgKHLryH0YMRFVvr8NIweoXhqPPnNQWNnWRO1YwNUTEYUzY6+PqZNTGWlY2P4gK3btiIH40qp+NgxUq3IOHjgIP3hMIyNYjT6GM/R+Ezt2bePkEVdiD9DkxMTCdv9vbuSnp3e3f2xsnA4TPGUJAj07+lPeD7Fvt2x65XqGojjR2P1dvTuij3v4XCYMsxV+uYtW/Bv3Zwk2/btO5ge136FDh3B378nVmd0bIzNr7+e1HdGdvX2MhWKZKcya3//4cMUA6NDgzOiK8CGQhdCNBExrejZS0SJVwKj2v8BfQUpZTTk3y+BT5kde+PGjQ7FjRAOh+nu7mZ05CgHdOXz7n809nm8yMc+7XNzSwvl3ckhOAf+VIc+Wsai7m4OdSxgZHckJGZ7R2fC8aN0mxxLL5ebRHPYt3/udvzVtfQHAtiy8grBGYc2MXdoL79qv4gV+19BALU11dToZJzoL2evye5Bny+hnXUNpzPuL+LS/heS6hYVF2P2EpmuL9zor17dZx+SzsF+tpW3JtSZP29e7F6IIjT5eoERfzGHQtW0Du+PldmR38iB8gpGgdbO+RzX8tP6ggHk2BjtQ3vYWZbsRQKws7w5bm4BOrsWxEwtevYEi5gEurq6IjFONA5WVmAy5oxRV1dHVXc306Oj7DZsM57jyNH9HNR9b2lt47BFXYg/Q4FgkAW67SPHDiQcB2DuvHmxsu7ubg5VVTIMtLS0Uqa/JytKYvdkqmswPTQYO5+Ozvns1+3TazFaPnnRIqbqatl7350J5Z1dXQnxaUaGB2KyhkLFtC9alNR3RtrndRDS5O012d7UPo8jQLHfTygUyvje7++3DreRdrAnpdyrs4VH/64FngNWCSH8wFIg9h4uhCgSQkTfh84HHATMdoDt1GB26/mo/eQX49894I4XJWZeSGVy0U+mBYL4kFzVu46d5S38peYk830sbOj6cgn8rmU5Xcf76BqM3NY1H7nZfL88EQ2fOn+wj+3lLUzpX8FNTS7xsj7NpbBteH9yPScI3b/oZdLMdh2De9hR3pxkiZgUPvpKG2LmIiB9/JQMTS62MvxoxxZl5QTmzSd0zoVpdrDA5D5Nat8qibrdlZl6G7qhTyyzBgkfwbkd1N9hyOqVFMIgkLjN1lxCGpdizUQ0kyaXbGzo9wLvAtYDP5JSjgshLhdCvBmoAZ4TQqwD3gJYRxHKhpQ3dgZuXpr/aXRhQy78Rp2SKpKi0IVCiN6QK/a9wpzRAX7aeWmkR3TKbfSVF5js3WHejm4icVNVBzvLm1mtG52LkN4kkH+3xUBLxOd+wfE+xv1FMSUNWCj0+P0RtW23DZu9jzlAHy45ukglEFfog8HSpCQcvWWNTPoCdAzGFbq1p5W5wkibhDymOBMf99JVbzY5WORYwdZ5NH/3YcvAc0nHTio2C/trOC+L8NK2J2J1Xi7GfSrXvo/qf/hHRj50Mw3fvF93cE2U0kQDSpLd3nANnE4Om71hxSZaZ3BSNGMbupTyGHCloexx3dcUEaJcwq7Dv92JHeNFM3xv+fFvTTMG5RSLB6jk/FVUrn0f+27SYltrD09QTnHtjt/xnUVreaGuh9W6H4QDX0iReV2n0H85dyUVE0NcpMtDGWiysWgph0TdObuOR94gtla0MU+Lk26amUbKmNLpK23APz1F04jRSOCM6eORhMq+yqr4ddJ+ZDu0Efj28hbqxuKJlzdVdgCw6Jgu842VQs82lITufvbPqWfOp79iUskiUJZTbMVxtzgf2891ilguRBZC9YXDCL/eQSK+KjVVm8ZJUadp8uq/+m0merex/7Px8MNRN0uvernkn5QrRW3WS6iW2n3KP6fecpVpzrB4q6v9+OcJ6myA+lHeyn0v0za0jwfnX86UzUQy0ZtuZ1kjL9Ut5s19zxCajr8qBjt1WeXzuLCo+sOfjZy39kPVPHyA0OQYW8t118lihD76YmRZf19pA00jB7NOyxaNG++vqSN6oXxaqsL5g/0EpicJV3Uk7LO5ah41Y8cS0t6lfTM09Le0GxPAznWKVsnWnGlmckkaoVtEI83EVTKFwpUJN71mUkoXKiAhOJfFYqQU+CoqKV68NLEFbUWs7eTUGTC7FbrthUUZ2hwzXlQxg1gpHWOQJd0N6ZfTvGv74+wubeD/TdhLzye10K0PdV5GaHIsKW2Zr7Sc+q9p8V/yqNAr3vIOmr79UGzZth/J/MHdbK1Io9CB4z//MQC9ZU3MzdZ+DjEF5auZE9dVmmmqWE6x4Hgvf6uen7DL5sp5nHxsR6I6sejPqDuksIgyanSXNJI4YEkR68QVHGRaSmozAxlSPav6eaJYSkfj27hhhG50WbajC9K4FEeT0HvVhp53UtoOEzrX/ihC2xD55wEbesXb30P51dfFC3Qyl158RbxckPBgGEdDZx/8G6cf2sR9E3PpP5rKJyLO36rm86e6Ht6+6ykqJoeTK0TbyL8JPeGHrut4HzvKm2MTo9aLpyYZ9hezp7SOzsF0Pgzpqb/tbqo/8MnIymDtR0REg5PJaRYPbGdbeRsj/sir94HGLvaVzKH76I7EA1kojzm3fIPRd36YQIN5Io/qG26ykMzJykRh+J8h0ftUP9I1npfV85fBACzl4Gs6eYRufLaNobQxTIq6Ed5CaAuL/CmSYGTLrFboboW1NBLUwuomzKLnier3f4KaD35aVxJ/CILtXbpyw01nnNQBPrzlFwjgtsc3MZUmhMCYL8APTrqautEjXNm33qKWFzS5hk45dGkTozG/b6tQrtPT7CyPPFyduknJTAm2z9fl59QCaUVH01LSM7CVKZ+fjdWR6/ZiRWS0vuyQwR/ZaoReVcPkmeebbrNLqoiJqdq2xDIcgXY43YpL4yAjaipKSpLh0F6tHdxaxIQRuhYGRKcTGv7l3qTE8jPhECGCRdR94V9pSJFpKltmt0JPoVACLXPjtRzepHW33kXdV+9GhNwJ/eoqVno4zYQuQP3YAO/f9igv7TrCd9en9iR9YP6b6C1r4iObf0axVfaZmJte/hW7fjXoIm3EG51wtDK5yOlptpVHfMM7j2c/Qk88uJaTVfd6fcrANsonhtjQELGt/rHmFOYO7U2Owe5UoaUZgPub4n75xYu1OPmpEkBAwqMVaG6j/Ep7i/6MQiUsoU8aoccaNQrhsC1SDu70OX/NkpAYbd3G7bbFSVdRQMm5K2O5cGeCWa3QU4Wl9ZWVx2NHOLxB/FXVlJx5bhaBiWYO02XsgPHHzcr1bdWuZ7iqdowf/2kXP2+/2LTOY63n819tF3Bl73qWHnndWhir1GH5QGdyqR8bYM7oAOGqzsgmMy8XgKlJdpS3UDk+SO34MfM6GcujaasJ3USynGJpD9YlAAAYFklEQVTFvld4tn4Jj7aez+tlLVzW/3zSrm6Fxw20dVB/+/cp08c5T3dPm7TdfO+vUic/NiM6n6APTGe8Jy1NLu6O0Is6FiSVpRuBZxSkL11Yjhw8J97TWA4o6lhgHQpXTyavcHh0UjRhOJac9SVGijg0Hyvfy2WLGnhw/hV8q/ta9odqABgIlnHPwqv54cK3cvbBjbx322P2RDLcqKEzlruWqcg2OhOSIJLo47XqjkgPWZlcpqbYVt5K52C/++8Y0RH6ROIE2Nt3PUXx9Dg/WvhW5g7vZ9WeF5P3devBl9OETluW8AMRm+yzasLp4EcLYSvLKxPKQ8vOo+ikxVS97xPxQ9v1csngeU2Z29d0B5uhbh2QPnb+zCt070SfypDSFZdy6Bu3pK6UaT/qbrT6276domIOscipmWSHTHHDBoqKuPWSBVQ9/iA/b/871jWeQfnEMEOBEAJ4S+863rPtv/Cnc+OzePh9ldUps5/PCAZZu49uZ0PjUg4UV1Nn8VYzPh1Z2GM9R5CFOJo8gZa5TGzfEiuvGR/kjpe/w59rF7FiIEzxtInHg0sRDkuWr3R+bIcKvfzKtfgqq9nV0J7YTGk5jXf+mMk9fbrCxGNXrrmesVdfonjJWQYZnT+wPsMPihFRUhbz3IrskG6ErleN6eUJdiygqHtJ6ko5GKHPeoWeithqtQxH6FF8VTUJ2VHyioVCTzrHFCOM4Q1PUHb527h2x+9ZvedPbGg4jX2hWmrGj3PegVdtL4EXFh4RCd43OcK4MjFqR/9b9Xy6LbxcXi+qY9IX4ORju2ZAoMi/imvew8gzTyZsahveH+lji9R92d6vUarem7xwLL0pwZkZTfj9lF18BVgFmkrhWFC0sJvWh55MKs9kUjTdeTXd/SBjm16N13d5hF5z0xfTmspykQimoBV6dImtlbdKyXkXc/yR+633N8s0n2+sZDGO0FPckOObXuXYw/cBUDd2lKt712Umi8kNGuxYQMmy7DwxMsLQL/OG9lI1fpy/1J7MZZq/uZGNExH77iKj26CL8lj5jAPWcXQyXDdh6x6IKTIrP/TMmrYUKdX5W5HpupFQCcU9p5tuCzS3JS4KTBfXRt+XKRRx0eLTaLj9e5Z5dS2POUN40UjsGlOHI0u5Iyv3kik+6RTmPvaS9QE8qNBLL7wk9lmfjstpgKPjmkLPCt0Afc7NX8/+eNlgMKv4kCw9vIVXak9iZIv56PG1qk5ah/ZRNTFkuj0roiYg2/7RcTJy27OL3Yk7l5SPr6om+djpyPD8236+nvqvJEXpNiedyUUvQ5rcxbaUebrjuERhK/RDEdOB0cfULtHRRaC1PU3N3FHz8Vto+cnvmfvYSwnhPpNTcOXw0gpBoHVe7HNeMFGOZx7exGCwjDcq5yZtmxQ+NlV1sNhkdD50y520/N//n5080RF6KoXu4sIax6RI0QbuuaIKIQi0OHx+cnDvpjPRxJbpp8PRmi2l0LOi5iM3E2huw2fIbGMXX1k5dV/+FnVfvjN95Rwh/AH8+lFPbIMxHkYuFKtLqwpdwCwK5WmHt+CT0/y51phMC8JVnQwHSjjjsDH7Esjaupj3RuYCaf9ddsEzb8uBVrGbf9RV4sdsefB3NOvyFphhO8NYNqS5LqKomDmx/KXW8jhzb1QKPStKL1xN872/ysr9sOTsC9KHEPUArqy2c97ozLdhFxNFVDE5wikDW9nQuDRpIPWnusUEpyc47fCWpP1cEijyz+Z1CHbpfnQyvXY2roe/oZmyN709/SDFxWsbNXmKQBB/dS2BevPQBfG2c3Dv2lgJ6q9JPxCs/d+32W8zF79TM9+EIi84HOXUfeUumr73sLM2PLWwyHxkuXLvn9lbUscmXZTDKQQv1PWw5MgbCREkXRUnagKy2Te+kviq5JkcoQqfj9qPfs50sQ0wIyP0Obd8g5pPfolAo3nGpiRyMEJ35HFiVTVYZBlXx6JR+3UzRCn0QsXhKMdXWkqwfX76iomNaP8Ewfb5hM66gNqbvuTwGC5hMcF47oFXKZ0c4Tdt8cw7r9SezMFQDRfvTTEhni2xNTNpHuLoK3sqF9R84KLu8VfXUn7JVflpPBusQhNkSC7emj1w5yiyQR8DPQGnN0+aYF2pEAJEIED9rd8yzdSSEzSFWH/HPYSWnRcrDk1PcGXfBp6vP5UtFXOZQvCfHZdQN3aUsw6+NpMCRf6lUwZmr/45mf+wwAsOXV5Zoa0tLvJZhBhx7FfuZT90IcQVwJ3AQSnlBSbb3wV8FDgMXKdlOFK4TP3Xv5+4Gi+KQ6WQNKkYCDoIxJ//EVVMfpOH5qredTzZdBb/3PMe5g3uYWvlXD5z6A8Epa1025lh0+QifD4khiQVLkyK+qqznNTN4zXNxQKcKD4zBwONopNPoer6j1Hm6O0iv2TzU/g8cJrZBiFEEPgwsAJ4APhQFu0oUuCvqKL4pFOSyh2/3hlMFrYSCntiOKcRjT9uct6lU2Pc8up9lE6O8rfqLq7b9jgXT7kcXTEJm5OiZl4SLii0xjvvz3BPD13TGab+69+n6e6fWG4XQlC55vp4gvakCg4b9LINXUp5REpplWBzIfCqlHISeAI4N9N27BKcb5HV/kTFYoReuuJS0/KkCSufjZc3q+BK+SAqi8XreufQHu568Zs8tP7zXLPrKcvXaNeImrBMuqb+9u/FPsthbVFTirg8aTGpn/NYOrOQ0JJlGa9RiVBAJpc0VANRE8tR7XsSPT09sc9r1qxh7Vp7MZdHR0cJ62NH3PZ9CAQ5YhVPIg3ROG3hDPe3lCsPRM9lcGjY9OIePXYM42LsyZOX8PqRY3DkWGz/aVLfruFwGN+O7ZQCIxmet5v9VTI0hB/YsXMnRYODljd29JyOT0wl9QNEzssNuUonJvABb7yxFeNPx87ePoyR9keGR4iO1VO1bSZbsXZNd/ftJhp9PFP5/Tt3UgIMDQ87Oobb977dZzJdPTflirYlQyWI0RGmp6dtywfw+pYtyIoq1+XSk1ahCyGagJ8aivdKKa9NsdtRIBr+rBIYMKu0ceNGOzImEQ6H6e52b/KtV/uf7THdlisToudSXlHBqMn2qspKjMnkKpua6dTkju7vDwaZTpGprru7mzE5zn6gpKSUjgzO283+2hsKMQF0zJ/PsfJy03PXU93YhNmC/+7ublfk2u33MQ0sPOkk+g3bOjo6MIY/KykJMa6TwQoz2Q5WVDACtLa1Ek2Vkan8o2PHOQCUlZXF7gk75OuZTFfPTbmibQUqqpgaHcHn99uWD2DhwoWxBWvZyNXfb7yj4qRV6FLKvcBKh21uAXqEEH5gNRF7uyIHhM44l9GXn3P2emdW15angYfsrTEbuvV5l1+5hsFHH4nUSxu72iWEYM7NX2d6eIjB3/6Ciddfm0EvFheP6wUzmkfx19YxdWAvFWve62zHHMSEysbLZRlwBxHF/QRwJRHF75dSPiaE+HdgPXAEuM7yQApXKVl+EaMvP+cwJ6LJw+vGwoscEmhoYmLrJkRxieXvTMn5q3Km0H2lZUwfOQRCxAKqlV50OXJkiMm9JhOyXvlt9FAgOq9SfMrpNP6f+/MthikZK3Qp5UtERt96Htdtf4CIh4sil0RHf06UrFldOwrd5YUX2VD7qVsZXfEswfZO60r6zD1FNiPkZUj9V+9m5LmnE2LC+EIhCIWY3Gf9yuwd8n9NvUrV9clx5m2Rgx9Lj3jwK9wjqtDNL600GwqaKGQhfNT/03dpuPM/rJvykJeLr6zc0oMnhl7OTHJGOiDQ1ErF2/4+vRwxsnnY3VQU3hmhl668PN8imCJSpHdMhTEJy0xQ2AkuTkSitm+/L25PT4epDV0QWnp2Uk5MPaY/Dh7GrZCwWTNTP4BuHjbPXdX2qEmu1VmOryJ1mjxX2pjxFhS5RRfPuv62u6l85wfT72Jqb9eOkyKFWLymRxRlOhIUaT5XQpo8dl6xXXtEDCFETleM2qHopMUUn3pGRvvOfewlfMWh9BWzRI3QC43oQxAdqdtRFClTlaXAK0rILl7RDzPl5eLm5fCYMrWi8VsPML5tc27autM8laGXUAodIjfvbFNOFsRGyzGFbOO8zEbotiZFvWNDt4e9PJEzj0nbHrn/cmHndZOihd35CwjnQZRCB1of+WO+RXAfB/lQzbKupHzdjR07+t2hbPnCK1H8vGxDj6bOmzU/0go9SqGTmFxgthOPOujTvtsZoZsoulRBpZLqe+vhF5bXU5h+zDmmsdvdGRlX/K93M/zHLPOhAl67prOV5nt/xfTwYM7aUwq90EiXnNjMa9HMDSuFnTde35uv5zUfvZmizoUcfeB7iRsENP/HY8ihQYafezo/wgFycjK5zKWurL7hJqpvuMmdgymyJtDcltP2PPIOqnANYxzuTCdFdSO06hs/k3hjGut77PXcX1lN5bU3JG8QgkBdI8F5Xfn1zJlOVuhZMcvs3oqZQyn0QsOY6CFDG7rerFLx1nfSfO+vkrfNsklRr9iFzUborihlF84v6qYqZjq8sGJGUAq9wJDRONxO3BZNzDOplF/Mb116e1Y0eeWoR+SccnmE7iLFS5ZR9b6PU/uxW/ItiiIDlEIvNKQhc48the7QbTE2ovf2q/6cf7w9sUCYuy1WvuP9OZIoQqApYr4qu+xt1H35W5FCj5hNhBBUXvPenKxqVLiPmhQtNDIwg5iaXFJ4uSTV98jANy0WcgY7c5vtKtDcRusjf0SUlDLxxqasj1e59v2MvfZXik890wXpFLMZpdALjeikaNYrRVNo6ZjJxZlo+cdiYVEefpCSU+Bl3plFC7tpfdANV0XFbEeZXAoNw6RozC/dQHHPGfhq6yJVU8Ry0dP47YciH7T6cpZNilrKmU/5Z0vfKWYFSqEXGLEBuWGEXnbJVTR9/2exeuVvvobGf/khwXldpmFKS5avSCqLJoUQ/sixi3tOp/j0c6i58TMunsEMYhmcywNKVUpaHnicph/8It+SKGYxSqEXGtNTQHJEv+C8LoJzOyi/7K0AFC1aQqCplabv/mdCEoYoldd+wOTYUXNOZITuKw7R8LXvEGyf7+IJuIuvqib+xWo0nMeQAFH3wEDbPPy1dQRb2/Mmi2L2k/GdLIS4QgixSQixwWL7ZiHEH7S/xZmLqHDEtMEMYjCLhE5fztzHXiLQ0JTyMGYrTaXmbmc6iepRmu75WUxhCwu7efSjv7Eld4JpBFvmUveVu6j9xBdz3rai8MhmUvR54DTgSYvtB6SUK7M4viIjLPzQ3bDVRkfogfQx0r2Cv6KKQGMLk3v6UtrQW37yBKI4R4mjDZQsOy8v7SoKj2xyih6BlAtQaoUQ64AwcJOUcjTTthQOsFr674JCD3YsoPzq66i4cm3Wx8opsXO38nIR+KuqcyqSQjETzKTb4gVSysNCiFuAG4G7jBV6enpin9esWcPatfYUxejoKOFw2C05XcMLcgX376cYOHT4MHvCYYoOH6YI2LdvP302ZCvX/luexwVXsHfgOAxkf5656q/SiQl8wNbt25BDkXFFtJ8Aevv6mKqMy+GF62iFV2VTcjljpuRKq9CFEE3ATw3Fe6WU16baT0p5WPv4S+BTZnU2btxoR8YkwuEw3d3eC2rvBbmOvrKeY0BdfQNV3d0c7ezi2LNP0LxgIWU2ZNt3ylKmjxzKyXnkqr/2FBczCXR1zifY3gnAsdde5Ki2fe7cdkp0cnjhOlrhVdmUXM7IRq7+/n7LbWkVupRyL7DSSYNCiCJASCnHgPOBrU72V2ROUdeiyH8ti0vlmuvZPzFN6UWX2dq/8Z/vnTHZ8kfUvGKxeMcDXosKhRtkbHIRQiwD7gB6hBBPAFcSUfx+4CXgt0KIQeAI8PfZi6qwQ8k5K2i+79GYF4sIBpk8e4V1fPQTgXSRJ1Ml81AoZhHZTIq+BKw2FD+u+5xZemxF1qRzSTzhiK2atVLoaoiuKAzU0ERR8Jh7Yok02xWK2YdS6IrCJ6qwTXN5omzoioJBKXRF4WNmQ7eM66JQzF6UQlecAKTzclEKXVEYKIWuKHzS6HOl0BWFglLoisLHxOSSaHFRCl1RGCiFrjgBUCYXxYmBUuiKgsevZWYiWKQrtQjUpVDMYlROUUXBU/uZrzDy7NMUdSww3S6Ul4uiQFAjdEXB46+oovyyqxMLhRqhKwoPpdAVCqXPFQWCUuiKE5KSc5KTYCsUsx2l0BUnJIHmNopPjcSPkxOTeZZGoXAHpdAVJyxC83qRE+N5lkShcAel0BUnLlE3RqXQFQWCUuiKExY1QlcUGkqhK05YRJGm0MeVQlcUBhkrdCHEjUKI57W/60y2/50Q4jkhxNNCiLbsxFQo3Kfk7AsBCFosOFIoZhvZrBT9nZTyB0KIIPA88BPD9i8ClwKLgc8BH82iLYXCdUovvITQmefiKy3PtygKhStkPEKXUu7QPk5qfzGEEKXAiJTyuJTyBeCUjCVUKGYQpcwVhYQbsVw+DPzaUFYNHNN995vt2NPTE/u8Zs0a1q5da6vB0dFRwuGwMylzgJLLGUou53hVNiWXM2ZKrrQKXQjRBPzUULxXSnmtEOIc4E2AIVAGR4FK3fcps2Nv3LjRgahxwuEw3d3dGe07kyi5nKHkco5XZVNyOSMbufr7+y23pVXoUsq9wEpjuRCiFfgmcJWUcsqwz5AQokQIUU7Ehv6aQ5kVCoVC4ZBsTC5fAhqBX4hItLorgJOBM6WUPwT+Cfg9MAq8N0s5FQqFQpGGbCZFPySlXCilXKn9jUgp/6Ipc6SUT0gpz5VSXiyl3OWeyPDwww+7eTjXUHI5Q8nlHK/KpuRyxkzJNSsXFj3yyCP5FsEUJZczlFzO8apsSi5nzJRcs1KhKxQKhSIZIaVF4twZ5sknn8xPwwqFQjHLWbVqlWlalrwpdIVCoVC4izK5KBQKRYGgFLpCoVAUCJ5U6EKICiHEb4QQzwgh3mOy/V1CiGeFEI8KISq1+k8KIdZpZRVaPdcjPjqVTSv7tRBiQAixWldvsxDiD9rfYg/JlVQvD3KZlbnSX0KIO4UQ64UQ/2Yo7xFCbNDkXOKkzA1ckOt+IcQLWv8kRT/NkVyfF0L0CyG+lqqeR+TyQn/do33foCtrEUI8pd3/q83aSImU0nN/wKeBdxGJAbMOKNJtCwLriSyKegfwWSAENGvbPwh8XPv8NFABnAN8Jx+yaeXNwK3Aal3dDfnsMzO5rOrl+FpayZp1fwFnAP+uff4ecJZu2y+BuUAr8GsnZR6R635ggcv3lFO5GoGLga+lqucRubzQX53a/4XAz7XPdwHnA+XAH5zK4MkROrAc+L2MhBT4b2CRbttC4FUp5STwBHCulHJUSrlH2z4BTImZi/joSDYAnWx6arU3inuEECGPyGVaL8dyWcngRn8tJ7J6GZLPr0ZK2Sul3E0kuJyTsmxxQy4J/Fh7G5qXD7mklPs0OUhVzyNyeaG/tmvbJojHuzoVeFZKOQgcd/qW7FWFro/WeJTEG8Fym4jEjvkQkdjstiI+5ko2Ey6QUq4AdgI3ekQuJ/LPlFxW9d3or1Sy6J8F4bAsW9yQ6zNSyvOAbxCJsZQPuczwQn+Z4aX++jqRkTmAX2pDdZNjpMWN8LkZIywiORKP1jiq/R/QbddHcoxtE0II4EfA56WUA0KICWxEfMyFbGZIKQ9rH38JfMojctmWfwblMpUh0/4ykOr89CO4aYdl2ZK1XNH+kVJuEELckSe5zPBCfyUL5ZH+EkJ8EnhNSrnBuM3kGGnJq0KX1pEcPw2sEkI8DCwFNuk2bwF6hBB+YDWRbEkAXwWekVI+pR07q4iPLstmPEYRkTUAY0TsZVu9IJeDejMpV1JZNv1l4Dkib3APa8e+X7ftsIhMnE8TH2XZLcuWrOUSQlRKKY8JIU7GoRJwUS4zvNBfSXihv4QQlwLnEZkrivJXIcS5wF+BSimlsz5zc1LAxcmFSuBR4Fngeq3scuDN2ud3a9seA6qAFmAc+IP29w9avdVaJz8NtOdDNhmf6NgGvEzEXNCofV5HJDlIhRfksqqXB7mM19e1/gL+jcik691AE5E3OoAlwDPa31InZS7dV9nK9Rtgg3aMnjzJdQPwZ2A7mhOCR/rLTC4v9Ndm4EUiOuserawNeIqI3rrUaftqpahCoVAUCF6dFFUoFAqFQ5RCVygUigJBKXSFQqEoEJRCVygUigJBKXSFQqEoEJRCVygUigJBKXSFQqEoEJRCVygUigLhfwAC1YIwsEQqnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(detection.times, detection.data)\n",
    "plt.plot(signal[0].times, signal[0].data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the likelihood object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = CudaLikelihood(generator, \n",
    "                   data = detection,\n",
    "                   psd=noise.cuda().rfft(1),\n",
    "                  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%timeit\n",
    "masses = np.linspace(0.5,1.0,200)\n",
    "likes = [l({'mass ratio': m}) for m in masses]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7738693467336684"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "masses[np.argmax(likes)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbfd65c8b38>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD+CAYAAAAeRj9FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhcZdn48e9zZp/J0mxN2pSmtAVSmi6UtpSlUloqi4CKtqCggq+KuPBDX1RURBSFF1BBRfHldUNUEBQEkb2l2NKWpQslNN33pEmTNHtmOzPP74+ZDG1pumSWM5Pcn+vq1VnPuZ9OcveZ+zyL0lojhBAi+xlWByCEEOLYSMIWQogcIQlbCCFyhCRsIYTIEZKwhRAiR0jCFkKIHJHWhK2UGqmUWq2UCiil7P285kKl1JL4n71KqY/EH5+vlFocf/z0dMYphBC5QKVzHLZSyg14gCeB87XW5lFe/zowD4gADwNXaK0jaQtQCCFySFp72FrrgNa6re++inkg3nP+t1Kq6IDnxgJNWutu4EwgCjynlHpYKeVLZ5xCCJELMl3DvgTYpbWeC9wPfPGA5y4n1hMHKAdGABcBy4HrMhmkEEJko0wn7AnAlUqpJcB3geIDnrsUeDp+uwNYFi+HLI6/TwghhrTDXghMo43An7TWPwVQSjnif1cAIa11a/x1bwJfiN+eCmzPcJxCCJF10j1KxKGUehmYArwANAJj4jXsxcRKHgAfBp7qe5/Wuhl4VSn1H+Ba4DfpjFMIIXJBWkeJCCGESB2ZOCOEEDlCErYQQuSItF50XLRokdRbhBDiOM2bN08d7vG0jxKZN2/egN5XV1fHhAlDazSftHnwG2rtBWnz8Vq0aFG/z0lJRAghcoQkbCGEyBGSsIUQIkdIwhZCiBwhCVsIIXKEJGwhhMgRkrCFECJHSMIW4hiFd++g/uoLCG3daHUoYoiShC3EMfKveIVoWytdTz1CuGE3TV/7NMH1a9GmSdeTfyHS0Xb0gwiRhKRmOiqlzgDuJbad15ta66+lJCohskDLj25CeX0U3/h9lGEQWPsmAP6lLxFpaSK0aT37f3477hnn0P3kX1BOF3kf+rjFUYvBLNmp6TuBuVrrgFLqL0qpSVrrd1IRmBBWCu/ahn/FEgBsw4opuOo6guvX4poyk+DbbxB8+03cM84h8OYyuvfsBMBsbrQwYjEUJJWwtdYH/oSGie12LkTO63n5GbDZ8M7+IF3/eJhodxeEQ+R/5JNoM0S0rZXS797N/vvvJLxjC9G2FiLNTVaHLQa5lCz+pJSaDJRprdcf+lxNTU3i9oIFC1i4cOExHTMQCFBXV5eK8HKGtDlLRCJ4X3qaaPUU9l10BZ6d2+h54Z9ow8YOhxeu+CLoKO1btsIFC0BrPA/cQWjnNvYdpS1Z2d40kzanTtIJWylVTGwH9MNm4tra2gEdV1b4Ghqysc09S55nf2c7ZR++Em9NDeZt99J0w1U4Rp/I6NOmHfY9rVVjCda9zZijtCUb25tu0ubj09DQ0O9zyV50tAN/Bm46pDwiRE4J1K6m/Tc/wTPzHLqe/AvOkyfimTkbAHvpcMrvexhls/X7fltZOZGlTehI5IivEyIZyfawFwAzgLuVUgDf1lqvSDoqITIssGoF4e2bCG/fhH3kCZTedh/K4Ug8bx9eccT328sqIBIh2r4fW0lZusMVQ1SyFx0fAR5JUSxCpF3PkufpXfRvSn/4C+KdDADM+l3YK0dTfNPt2MsqsBUWHddxbWXlseM0N0nCFmkjE2fEkBJYtZzA6hVEWvcd9Li5dzf2ESfgOnkitqKS4z6urTSWsCPxoX2RtlaCtWuSD1iIA0jCFkNKpGkvAOHNdUT9vYQ216G1xty7B/vIEwZ8XHtZrGQSaYkl7LYH7mbfd68n2tOdfNBCxEnCFkOKuS92BT60ZQOdj/yWpv++BnP3DrS/F/uIUQM+rsrLR7k9mM1NRFqbY5NuTJPAmpVE2lrxr3w1RS0QQ1naN+EVIlto0yTS2gxAaOsGzL27IRKh5+V/ASTVw1ZKYSurINLcSPcL/4RoBOX24H/zNXpfW4T/Py9R8eA/UtIOMXRJwhZDRqSlCaJRcDgJrnsTHQwC0LPo3wDYRww8YQPYy8oJbV5PcP063NNmYeQV4F/xCtrfC0Dv4udg+rnJNUIMaVlXEtGRCD1LnsfYtdXqUMQgEdxQS/cL/8SM16/d02YlkrWtopJoeysYNuzDRyR1HvuoMUSam1A2GwVXXYd7xtnonm4wbDjGnULPK8+B1km3RwxdWZewMQza7r8D++rXrI5EDBJdT/yJtl/dSTjeCfCePRcAe9VYfHMuBMA2vOKgcdcDUfip66l44HFGPPRvXNWTcJ9+Fths+OZ9iPzLPkGkqR5jx+bkGiOGtKxL2Eop7CNPwJCFdESKhHfvgEiE3iUvgFJ4Zn4A5XThnTUHV01sqrkjifp1H8PrwzH6xMT4blthEeU/e4hhn/86nrPOA4cTe+2qpM8jhq6srGHbR47GWP82AMHaNdjKyrGXj7Q4KpGLdMTEbNgFQGjDOmyl5Rj5BZTf/1fsZRXoaBTsduwjR6fl/M7x1YnbjpEnYDbvTct5xNCQdT1siP1gq7YWdDhE82030vnIb60OSeQoc+8eMM3EfVt5rE7tqKxCOV0Ybg9lt99PwcJr0x6LfeRojJbYN8eOP/+vDPUTxy0rE7Z95GhUNEpg9Uq0vycxFEuI4xXevQOIXWgEsJe9/8Kie/L0jEwnt48chWrdhw4F6Xzs9/Qsfjbt5xSDS5Ym7Fg9sfc/LwKxab5CDIS5ezsAvos+BoC9PLmRIMmwjzgBFTEJrHkdIpH3TY8X4miyM2FXxuqJ/tf/A0CkrcXKcEQOC+/ejq2kDM/0s3BNmYEr3tO2QqIj8tpiAEnY4rhl5UVHo2AY2u2B+ISDaEcbOmKibFkZrshi4d07sI86EeV0MfyOByyNpW9iTl/tOtLaIutni+OSlT1spRTR0thiOsrlBq2JdrRbHJXIJf43l9Hx1wcx9+zAccIYq8MBwFY6HG13oHu6Yg9EI0Q72qwNSuSUrEzYANH4cpWuKTMAiOyXsog4NsG6dbT86Bt0/uVBtL8XR9VYq0MCQBkGujh2cdOIr7dttsh8g3TREZPQ1o2YTf1vuZVrsjZh63gPu2+bJqlji2MR6Win5Uc3YS8rp+LXf6P4ptvxzv2Q1WEl9HVE3NPPApARUGlitjTR8JkP0XTDVey7+Tr0IFkSIGsTtjn1DPIuWYB76kxAetiif1prQju2oLWm5+Wnibbvp+Tbd+GoGofvvIsw3B6rQ0zoS9ieGecA8QWpRMr1vvoC0bZWvPM+RGTfXsxd26wOKSWyNmFHK0ZRdP23EuNjozK0T/QjWLuapi9fSfe//07P8//EOXEqznGnWB3WYUVHj0f58nBPPQPsdulhp4n/tcU4xlVT+KnrY/dXHX6rWbOxHh0KZjK0pGRtwu6jnC6UL1/GYot+Bd+Jrc/R/uBPMBt2kXfhRy2OqH/mlJlU/uVFjPwCbMVl0sNOkg4F6Vn8b/QBs1nNliZCG2vxnjMPe1kF9qqxBFYtJ9rdhdlYn3hdpK2VxusX0vmPh60IfUCyPmED2IpLpYYt+hWqewdbWTnK4UT58vGcPc/qkPqnFMrhBMBWUiY97CT5V77K/p9+n66n3tsL3L9sEQCe+KqMnmlnEaxdQ+MNV9N046fR4RAAPYueQYeCBN9ZnfnABygnBjbbikqkhy0OS0ejBDe+g3f2fDxnngdmGMPltjqsY2IrLSe8daPVYeS0cHwma+cj/4etpIye558k+M4qHOOqcVRWAeA+/Uy6nvxzbKKSGSaweiXumbPpeeEpAEKb3s2Z8fC508Nubab5thtpveeWnKo5ifQy9+xA93Tjqp6MZ/pZeGblzo4usR72vkEzgsEK4T07MPIL0abJ/ntuwWxqoODqL1J2232J17gmn07B1V+k/Gd/xMgroHfpSwRrV2M27MJ12hlof08i8We7nOhhG0UlRJoaiMTHU0Za91F2+6+SXnA+WZGuDjr+cD8FCz6T1AauYuCCdesAcE6YZHEkx89eOhwdDKC7u1D5BVaHk5PMPTtxnjwR3wUfIdrZju/8SxIlpz7KZqfwE58DwHPWefQufZnQpncxCgoZ9l830vSVTxDa8A7OMeOtaMJxyY0edlEpAM6JUyn68s0E31lN4O03LI1Jh4K03n4TPS88Se/SlyyNZSgLbXgHI78Qe/zrby6xxbckM5saMPfuYd93vkSkU2b0HisdjWLW78I+qgrv2XPJu+jy9yXrQ3lnz4+tANqyj9Lv/RTHmPEYBcMIbViXoaiTkxMJ21E1FuVyUXTdN/CedzEYNkLxDQ6soKNRWn92G8F316A8XkKb11sWy1CmtSbw9ls4J0xO7PKSSxyjxgCxOqz/rdcIvv0GwbXWdkRySaRlHzoYwDHq2P+zdk2ZTt5lV1J62324Tp2KUgpn9SSCG2rTGGnq5ETC9kw/m8pHF+McdwqGx4tj7MkELUzYHX/4Jf6lL1H42RvwzDiH0Oa62KSNJc9LDymDwts2EWmqz6m69YHsI0eDYcPcvZ3wji0ABDe9a3FU2S+0bRP7bv4CwfVrgdjmx8dK2ewUXXcT7snTE4+5TqnB3L2daHdXqkNNuZxI2BAbj93HdeoUQptqDxp72Z9UJ9LQ5jq6nniYvA8tIP/yT+E46VQizY0EVi1n/z230PX3h1JyHnF0/uWvgGHkbMJWDgf2kaMI79pOeGdsg+CQJOyjCqxaTvCd1XT88X7gvW8qA+WcMBmA4Mbs72XnTMI+kOvUKehgkNDWDUdNxOHtm9l/zy10P/PYYZ+P7G+h+7knjvlKvf+t10ApCq66LvZ16qQJAHT8/hdAbK3jvmNFA35C22WX7FTRkQjh+l2J+70rFuOaOBVbfCGlXOQ44UTCu7clEnZ4Sx06cvSOyFBm1u8EINLciPL6MIpKkjqe86RTQSlCG99JRXhplZMJ2zlhCgD7f3YbDZ84P5ZE+9G3CUKonxpVz0tP03b/HfS+8ly/xwiuX0v9J84ntGUDwTWv4xhXja1wWCyW+BTo8M6tKI+XSGM94W2b0FrTeue3aLrhqoOSjBg4//JXaLzuY4S2bCBcvxNz5zY8Z821Oqyk2E84EXPPTnRvD66aaehgkPCOrVaHldXC9bswCmK/f45RVUlfvzC8PhxV4wjVScJOC3vpcGzlIzH37EB5fHT+9f/67SEnEvbG2sO+JrxnBwDtv/850d7u9z2vTZP9999BtLOdjj/9muCGdYn9AQEMbx72+EWPwmu+CoaBf/kr9Lz4FIG3lkM0Sve//kawbh2t9/5AFrFKQnjPDtCa7mf/TvfTj4LdnvMJ23HCiYnbvvmXAVIWORpzz048Z87Bc+Yc3DNmp+SYzupJBDfWoqPRlBwvXXJiHPbhFH/t+xAOYzY10Hb/HQTXvoH7tDMOeo3Zso/w5vXYK0dj1u/CbNiVmP2UeM2eHdjKKoi0NNH+2/so+up3D/ofu+vpRzF3bsNZPYnAquUA7zuP85RJRPa34Jt/Kf5lL9P594fADOOadDq2snJ6Xnqa3lefJ9rZQXDN65TcfAeuU6em6V9m8Io0x9bd6F3yPDoaxTfvEuylwy2OKjkHJmzPrHMx8gsJrHsT34UfzcmRL+kW6eog2tmOvbKKghtuSdlxndWT6Hn+Scz6XVmz4cXh5GQPG8A96XTc02bhO/8SbCXDabnzZlrvueWgtRn6tmLqW7ErdMhFBa014T078Zwxm/yPf4aeF/5J1+MPEelsR2uNjph0/eNh3NNmUfLtu8BmQ7ncuOIXKfoMu/arDL/ndxguNwVXXItn5jkUfup6Sr59F/kfuQod8ENUU/Kdu8FmY983PkfrPd+TWuVxirQ0oXx56GAAIhEKFlxjdUhJs8eTg62sHCMvH8/s+fj/8xJtv7oTHYlYG1wWMuPlxeMZyncsXNXxC49ZPh47Z3vYfZTDSen3fkrXv/6Gf9nLmA27KLvzQfzLF9Pxu/twnHgynrPmxsZLb6jFd8Bi9tG2VnRvD/bKKvIuWYjZsIuOh+6n46H78V38MTwzZxNtbyXv4m9jLx1OwYJr0OHQ+wbn24pKsMUvfLhPm4X7tPdKJrbCYRR95ds4xlXjOnki7mmz6Hj4AbqfegTfBy/DKBjG/p98j9If/Dwz/2A5zGxpxD3pdLQZxj5y9KCYXWq4PdjKK3GMjvW0i774DQy3h64nHsZVMw3fnAstjjC7mHtiFxztKU7Y9srRGMNK6HnuCXxzP5S164rkfMIGcJ40gZKv30bvmXNo/dFNNHzyfHQwgGviaZR85y6UzYbz5IkEVq+g65nH8M6ej62wiHD8arNj1BiUYVDyjR/jP3dprAb97D8I1a3DKCjEPf1s4L2e+vHKu+hjiduGx0vhp66n+9+PE1i9Em2GCe/YQs9L/4KpZyf/jzGIRZqbcE+aTtH137Q6lJQq/c5dGL48AJTNRuG1X6V36Uv0Ln5WEvYhwvU7wWbDXl6Z0uMqw2DY578WG5r71CMUXH51So+fKjlbEjkc75lzGPbFb+CeOZvim26n7Me/xjasGAD39LMx9+6h/YG76Yqvf3vo/9bK4cB79lyKvvJtjOJSwts34/3ABSlfs8TweHFVTyawegX+FbGyTe/if4MsAtSvaG83urcHW1m51aGknHN89UHfFpRh4Jt7MYE1K2X51UOY9TuxV1Si7Knva3rPvQDPrHPpfPgBor09KT9+KiSdsJVS9yqlliqlsuI7ff6lV1B68534zrvooERbcPnVVD7+auxqcHyGVLh+J8rlwlZ6cBIwPF6GXXsDGAa+D16Wljjdp52RmKnnPHUKZsNujJ0yZrs/fRccD/2sBivv3IshGqX1Z9+n6evXJL4NDmVaa8I7t6Zt3RilFN55l6BDwUStPNsklbCVUtOAPK31bMCplJqRmrDSw/D6cJ06NTaVPBzC3LMT+8jRKOP9/wy+uRcz8i8v4hxXnZZYXH1DA5Wi+Mbvo1wuHEcYTz7UmfGEbS+rsDiSzHCMGoNz4lSCa98gtLGWQD9bXA0lXU88HBvSd8YH0naOvm86ZuOetJ0jGcn2sGcBfUvVvQycmeTx0s45YTKYYUJbNsS+Xh3hf2tbfHB+WuIYV42RX4jzlBoclaNxTZmJbduGtJ0v1/VtpTUYSyL9Kf3uTxjxu6cw8gsTMyGHomDdOtp+ey8df/wVnrPn4rvgI2k7l70iVhs392Znwk62EDQM6NuOuAOYeOgLampqErcXLFjAwoULj+nAgUCAurq6JMN7P2X34APq//QbHA276Zl2Ns1pOM+xsF31JbQ3j/11dTiGleFsXkbdmtWQRbt8p9uxfs7OundxKMXmfS3QmrsLbA3k59pTWkF4w7s0WfRzmqykfpejUXzf+xKYJpFTJtF84RU0b0hvx8abX0hLXS31Sfx7pyt/JZuwO4C+ldcLgPf9JtXWDmxBlbq6OiZMmDDwyI5g74hRsO4NjGHFjP3M9RheX1rOc1QHtM/f1UrLi09wogPcaWp3NjrWz7n1+QjB4lIm1OTeRgUHGsjP9f4JNfS++gLV1dU5OZkmmd9ls7GevcEARTfcQl4ae9YHahpVhfJ3MzaJ38Nk2tzQ0NDvc8mWRFYAfTueng+sTPJ4GdG3OlfBJz5nXbI+hHN8rFYe3pybvah0izQ3DZkLjodyVI1D93QPyREjfUtHZHL2oX3EqMFZw9ZarwYCSqmlQERrnROrr/vmX4b3Ax8k74KPWh1Kgq24lGhBEaEtUsc+lI5GMRvrh27CHj0WAHPXtqO8cvAxd+8AwF45JmPntFdUxjZHiO+unk2SHsyotf5/qQgkk9yTpx+0gHm2iJ5wIqEtdbEForTGVlJmdUhZoeuffyXSVI/7ys9aHYolHFXjgNiKkAcuPDYUhPfswCgoTKyOmQn2EaNAa8zGhqxbV2RQTZzJdZHKMZj1O9n7hctp+Z+brQ4nK4S2bqDjofvxnHleYjW7ocZWWIRRWER4KPaw9+w8rh1lUsFeER/al4UjRSRhZ5Fo1XjQGh0OEd6++Zg3VRjMOv78vxheH0X/75acvOCWKo7RY4fE0L5g3To6Hv5N4mc/vGdn0jvKHK9sHostCTuLRE6uoeyO3zDsmq+i/b1E21qtDslS4V3bCbyxlLxLFmLLL7Q6HEs5TzqV0NaNRAMBq0NJq54Xn6Lz0d8SXPM60a5Oou2tGe9hG8OKUR5vVs52lISdTZTCPWU6jjHxmmX8CvlQ1fXkn1FOF3mXHNvY/cHMNWV6bMJXnXWbT2eCuXc3AB1/fdCSESIQm6LuPKWGwJqVWfctVxJ2FuqbfZmN/8NnSnjnVnoWPYNv/mU5vWdjqrhOnQo2G4F1b1kdSlqZDbtRvjxCdeto/919AGlbO+RIvGfNjW16kmXXDSRhZyFbaTnK6RqyCVtHo+z/xY8xfHkUXHWd1eFkBcPrw3nyRIKDOGFHA34irc3kX3ol7mmzMPfsxDFmPPaKkRmPxXPmHFCK3uWvZPzcRzIo1sMebJRhYB95wkErtB1u44TBqveVZwltWEfxf/8go8O5sp17ygw6H/sj0d5uDG+e1eGkXN+oDEfVWAo/9UVLY7EVl+KsnoR/+WIKP/E5S2M5kPSws5S9sgqzfidmyz5af3Ybey4/B/+KJVaHlRGBVSuwlQzHe97FVoeSVVyTp0M0QrB2jdWhpIXZEKtf20eOtjiSGM9Zcwlv24QZX3gsG0jCzlL2ytGYjfU03/Jlev/zIsrpouc/L1gdVkYEN9birK4Z0sP4Dsc1YTLK5U5sBj3YmA2xEqB9ZHZs/da3tHLfRifZQBJ2lnJUVkEkgrl7O2W3/gzv7PkEVq1Em4Nz495IVwehrRuIdLQRaazHeUrN0d80xCinC9fUmfjfWJZ1oxdSwdy7B2NYcdaUe/pq52ZT/4sxZZok7Cxlj2/KmnfpFbinzcI94xx0TxfBQTqsq+tvf6Dp69ckyj6SsA/PM3M2kX17MXdtw2xqQEejVoeUMmbDLuwjT7A6jARb6XCw2TAb660OJUESdpZynnQqpT/4BcP+60YgtqUYdjuBN5dZHFl6hHdvA9Ok4w+/BMOGc/zQWWL2ePRtCL3/V3ey97OX0f3MYxZHlDpmw27sI7InYSubHfvwEVk141ESdpZSSuGZflZiX0rD68NVMw3/6/8ZnF+H4xecot2dOKrGYQyhTRyOh710OI5xpxB6N7YvaWB1TqxofFSBt98i0tqMI4t62AC2ikoi0sMWA+E953zMPTsJbXzX6lBSSpsmZmMDzurYOuXOU963cZE4QMGCa/DOuwTveRcRfHcNOhKxOqSk9C59meZbvoz9hBPxnn+p1eEcxF5RKSURMTDecz+IcrnpefEpwru20fPyM4Oit2021UM0Qt6FH6Xgk5+XqehH4Z09n5Kv34Z7+tno3h7C2zdbHVJSup/9O/YRlZT/7I/YS4dbHc5B7BWVRDs7iPZ2Wx0KIBNncorhzcMzez69/3kB/2uLiHZ3Eli9guIbb0U5XVaHN2BmfXz8beVofPOzq4eVzVw10wAIvrMqsWNRrtHhMKGN7+C78PKs2f3pQPby+Ka8jQ04x55scTTSw845eR/8MNrfi/J4yP/4Z+h99QVaf/K9nO5pJ8bfVmbHhIlcYS8djn3EKIK1q60O5biEtm6g+dYbaLnjm4S2bkAHg7gmTrU6rMNK7KKeJWUR6WHnGOepUyj+7x/iqpmGfXgFRsEwOn7/czof/V1WTaE9Hn0L/hgFMg39eLlqpuFfsQRtmih79v86m82NNH39WkCDaUK8o+GaeJq1gfUj2xK29LBzjFIK39yLsQ+vACD/8qvxzJ5P51//j2jAb3F0AxNu2IV95GiZ2TgAnjPnJEpjucDcvQPMMKXf/QnK68O//BXso6qwDSu2OrTDMvILUL78rBkpIgk7xyml8J49D6KRnF3dz2zYnXXDuXKF+/SzMAqG0bP4WUKb6+h87A9ZXR6L7G8BYmtc92355po4zcqQjurQhdisJAl7EOir/fbVgnNKOESkuVHq1wOk7Ha8H/gg/pWv0nzrV+h46FcE175udVj96kvYRlEJ+ZdegXK58cw42+KojiybtmeThD0I9K1uFs6iRWqOlW3jOxCNJsZgi+PnnXsxhEMQ1RiFRXT9869Wh9SvSFsLyuPDcHuwjxhF5aOLYmtPZzHHmPFE21qJdLRbHYok7MHAcLuxlZVn1apix8q+dgVGwTDcU2ZYHUrOcp48kcLPfIXSH/6CvEsWEnhrOeFd260O67Ai+1uwFZck7ufCcFRHVXzLvl3W97IlYQ8S9soqzIbcSthRfy/2d9fgOef8nBjhkK2UUhQsvAbXKTXkfejjYHfQ8/K/rA7rsCJtLdiKS60O47gkEnYWlEUkYQ8SjlFVhPfszOoLTofyr3gFFQ7hPfeDVocyaNgKi3CMGU9o28bEY1przJZ9Fkb1nmhbK7ai3ErYtpIylC9fErZIHXtlFbq3h0hjPYEsn0gRDfhp++297L/vh0RLK2IbzIqUcYwZT3j7lsT93ldfYO+1l2ZF0o7sb8HIsR62UgrHmHGSsEXq9O0s3fz9G2j+1hcS++Nlm9COLTTdcBXd//wrvvMvxf/lW1CG/BimknPMeKLtrUTa9wPgX7YoNuxz725L44r29qAD/pwriUCsLBLesdXyb7DymzJIOPqG9sXHYvtXZedEirZf3kG0p5uyH/+a4htuQecXWh3SoOMYMx6A8I4t6HCIQHyYX6S12ZqAtCZYt45IW2xIn62o5ChvyD6OqnHoni7r/g3jJGEPErayCpTHi2viadjKKwmsXkE0GMiqxB2u30VowzryP3qVjApJI8eJJwGxhB18dy3a3wtApNWakohtax37bvos3c89EbufYzVsyJ4Lj5KwBwllszH87t9S+v17cZ9+JsF1b9H+vz+h5davEtxQa3V4APS+8hwohXfORVaHMqjZhhVjDCsmvGMz/jeXgcOJcrks6x0a8eFwPS8+HYsvF0sioyVhixRzjj0Zw5eHZ9ostL+Xnhf+CYB/5ZKMnK5ul1sAABOKSURBVD+0dcNBtXNz316avvm52N6DWtPzyrO4pszIujWPByPHmPEE312Lf+lLuCdPx1Zabl3Cjk/r1j1dsfs5mLBthcMwikokYYvUc02ZDjYbRlEJzupJiY1tAXoWPXPQ/VRq+fG3aP/dzxP3u597gtC7a+l58SmCtauJNNbjm3txWs4tDuYYMx6zYTfRQICCq6/DVjLcupJI/Q5U31rXDidGXoElcSTLMWa85QlbZisMQoY3j6Iv3Yx9VBXh7Ztp/809hPfsoHfZy3Q+/Buw2yn/6R9Tuuh9tLebSFM9ptcLgI5GYyUQoGfJ84Trd6F8+XjOPj9l5xT98549j9CGWoq++h2cY8ZjKykjGN8HMpOivd0YLU3kfezTdD31CLaikpxdldFRNY6e5/6BjkRQNpslMUgPe5DKu/CjuGum4Zl1LgD7vvl5Oh/+Dd4PfBBbYTGtd30npcuxhnduA2I7c2itCa5fS6S5Efe0WUQa6/EvfQnf+ZdguN0pO6fon+vUKZT/9Pc44yNGbMVlRPY3o6PRjMYR3hbbvsxVMw3PmXMSI1hykaNqHDoYxGxqsCwG6WEPcvayCtwzziHS3Ejh1dfhu+CjBNa+QcutX8W/Ygm+81JzAbDvq6L29xDt7qR38XMot4fir91Gw2cvg3CIvIs/lpJzieNnKx0Opkm0sz2ja0/3zbh0jq/GPW0W5GbnGjh4pIhVywEPOGErpb4AfDZ+9xda6+xdImyIK7vtvoPuu6fOxCgoJLBmZcoTNkCkaS+Bt9/AffqZ2IpLybvgw0S7OnGMGpOSc4njZyspA2JD+zKasLdsIJpfmJMjQw7lGD0WgPDOLWDRCoPJ9LBf1Fo/qJRyACsBSdg5QtlsuKeeQWD1SrTWKakphnduRbk96ICf8PbNRBrrccYXqC+6/ltJH18kx1YSG5kTaW2BcZk7b3jbRqIjqzJ3wjQyPF5s5ZWWXngccA1ba70jftOM/xE5xDVtFtG2VsLbN6fkeOGdW3GfNgsA/5tLgfcmcAjrHdjDThcdjRJ4ZxXd//47ofgsy/CubUQH0Tcre8VIIhauyZKKi45fBJ5KwXFEBvUl18CalYnHdDRK1zOPEek8voXaI+37ibbvxzVxCsrrIxCfXek48eTUBSySYisqBaXSOha7d8nzNN98HW2//p/YyKQdWyESITJIetgARl4B0a5Oy85/1JKIUqoCePSQhxu11lcqpc4ALgY+0t/7a2pqErcXLFjAwoULjymwQCBAXV3dMb12sMh0mz0Vo2j7x8M0dfVgzjwXY9dWvA/czb63VxG6/JpjO0g4hOP1V3EBew0nzsJibHt3oz1eNrfsh9a2I759qH3OVrbXm1dA69ZNNMTPb+zdjfvBuzBnnkvo/A+Dw5nU8V1LXsTmyydSczp61WvsWvIiLsBfWjFoPmNXJIqtff9R25Ouz/moCVtr3QjMOfRxpVQl8FPgMq11pL/319YObFp0XV0dEyZMGNB7c1Wm2xz8+m20/e9PMB7/HSXjxhPxd9AOuN5axpgvfTOxSI+OROj6x5/wzJqDY/SJifdHe7ppuvHTmA27MAqGMW7uBbStex3/3t24x1Uz+tRTjxrDUPucrWzvvhNPItq0h7Hx87ctfZbuni6ci56mMNBN6XfuHvCxtdbs3bUF52ln4Jv3IVpeX4Jv1TIiXh/OEaMGzWfcfkIVXW8tpbq6+ojXfpL5nBsa+h82mExJ5FagHHhCKbVEKeVJ4ljCAq6JUym/708YBcPwv7GU4Pq3MfIL0eEQnX/7PdGAH6017Q/+lI6HfkXPS7G1IFp/dltsPetf/hizsZ6Sm+9k5J+exZZfiK18JCD162zkmTmb8M6tmHv3oCMmva++iGfWuXjnXkywbt1R36+1pu039xDc8M77nos0NRBpbsQ16XRcE08Dw8DcvR3nuFNgEC2fa+QXgGmiUziH4XgMeJSI1vq6VAYirKEMA/e0WQRWrwQU7ulngdZ0/+tvdD/zGEZ+AdHODgDMpnqivT30Lnom8f6CT30R7+z5ifv24SMASdjZyDPrXNp/ey/+N5biGH0i0fZWvHMuil0YXPws0WAAw9X/xKbQ5vV0/+tvEI3iqp500HOBdW8B4J58OoYvD+f4akKb1uMYe0pa25RpfdPqo12dGB5vxs8vE2cE7tPPpHfJ80Bshpxv/mV4z72A0NaNRPY3YyspI/Tu25iN9YlF8PMXXouRl0/+R6466FiOqnGgFM5Tat53HmEt+4hR2EePxf/aYgJr81C+PDwzzkaHggBE9jVinDCm3/f3rUET3rPjfc8F31mFMawY+wmxkplr0nRCm9bHetiDiBFfvz3a1QHDKzJ+fknYIjFiBMA5YQrK4cQzczaembMTj7c9cBfBDeswG2IJ2zt7Ps6x7x8F4po6kxG/fzrR0xbZxTPrXLoe+wMAhZ+6HuV0YY+XsczGPTiOlLCXLwYgvDu2I3twQy220uEop5PAW8txTZmRqOt6zppL93P/wFUzDfZ3pLFFmWXkx3vY3daMFJGELbAVleAYV425d3diNteh7BWV6N4eQpvWx+6PGHXY1ymlJFlnMd/5lxKqe5v8yz+V+A/ZXlEJgNnUgH/lq4Q2vUv+FZ89qDwS3rUdc89ObOWVRJrqiexvYd/NX8DIL8BRNY5obzcFCz6TeL2ruoZRj78auzOoEvYBPWwLSMIWAAy79qtE9jf3uwqZrSKWoANrXscoKrGkfieS56gczfD/efCgx4yiEpTTRaSxgZ6XnyG8eT3+t5ZTdvsvsRUWAe+tqZ5/+dW0P3AXPS8/A+EQ0c4Ogmtep/AzX8E5LnWrP2arA2vYVpCELQBwn3bGEZ/v+9oc3r4J50TZ5XwwUUphKx9JaMdmwls34po6k+DaN+h56WkKPh7rNQffWYWjahzu02YC0P3Ck2AYlP/k9wTXryXvkmObX5HrEiURSdgim9krRr53e4Q1K5WJ9LFXVBJYtRyiUQo+9mna9jUSXP82ADpiEqxbh/e8i2PlE7udSGM9jnGn4DxpAs6TBscY62NhuNwop4totzUlkcEzQFKkleHNwygYBoCjn/q1yF328pEQjYJh4KyehPPUKYTq1qG1Jrx9C9rfi2viVJTNjn3kaCC2xvVQZORbNz1dErY4Zn0Xp+wWrQUs0qfvs3WMq8bw+nCdOoVoZztm/U6C764BYhOtAByjquL3T7MmWItZuZ6IJGxxzBIJW0oig44t/tn2JWXXhCkABNe/TfDdtdiGj8BeFht37BgzHgwj8dqhxsgvtKwkIjVscczsI0aBUv0O6RO5yzlmPNhseGacA4B9VBVGXgGBN5cRXP/2QRel8z/8SdxTZmZ0I4RsYuQXYO7dY8m5JWGLY5Z36RU4T6nByMu3OhSRYvYRo6h8dDFGfHdzZRg4J0zGv/wVILapbx8jLx9XzdAsh4C1JRFJ2OKY2YpK8JzxAavDEGnSl6z7FHzy8zira/DOno+jcvCsaZ0sKYkIIbKO6+SJuE6eaHUYWcfIL0AHg0ddLCst587o2YQQIsf1TZ7R3V2ZP3fGzyiEEDmsb3p6xIL1RCRhCyHEcUisJ9Jx5O3v0nLujJ9RCCFymGPsyaBUYup+JknCFkKI42ArLMIxrprA6hUZP7ckbCGEOE7uabMIbagl2tOd0fNKwhZCiOPkPv1MiEYIrH0jo+eVhC2EEMfJVT0Z5fFlvCwiCVsIIY6TsttxTZpGsHZNRs8rCVsIIQbAOfZkzIbdiV3nM0ESthBCDIBjzEkQjSR2kc8ESdhCCDEAjjHjAQhv35Kxc0rCFkKIAbCPHIVyugjt2Jyxc0rCFkKIAVA2O/bRYwnvkB62EEJkPceY8VISEUKIXOAcM55oeyuR9v0ZOZ8kbCGEGKDEhccMlUUkYQshxADZSssBpIcthBDZTnm8AGh/T0bOJwlbCCEGyPDGEna0tzcz58vIWYQQYhBS7r4etiRsIYTIasowUG4PUSmJCCFE9lMeb+70sJVSTymlfpSKYIQQItcYHl9u9LCVUpMBT4piEUKInKM8XnSOXHS8Afh1KgIRQohcZHgzVxKxD/SNSqlqoBloP9LrampqErcXLFjAwoULj+n4gUCAurq6gYaXk6TNg99Qay8M/ja7zSiqo+2gNqarzUdN2EqpCuDRQx5uBDqBW4HqI72/trZ2QIHV1dUxYcKEAb03V0mbB7+h1l4Y/G1uHV5OqL2FMQe0MZk2NzQ09PvcURO21roRmHPo40qpF4A/AsVAiVLqJa31qwOKUAghcpTyeIlme0lEa30BgFJqDnC+JGshxFBkeHwZm5o+4ITdR2u9BFiSdCRCCJGDlNeLDgbRERNlSzqlHpFMnBFCiCQYiQWg/Ok/V9rPIIQQg5jy+AAyMnlGErYQQiThvR52+i88SsIWQogkJHrYGZjtKAlbCCGSoLyZ28RAErYQQiShrySSibHYkrCFECIJfSUR6WELIUSWkx62EELkiEQPWy46CiFEdlMuFxiGjMMWQohsp5TK2DZhkrCFECJJsW3CJGELIUTWi/WwpSQihBBZz5CSiBBC5Abl8cnUdCGEyAWG10t411Za7/0BZsu+9J0nbUcWQoghwjvnQuwjR9P78r8IrFqetvNIwhZCiCR5zzmfstt/CYAOBtJ2HknYQgiRAsrpBiRhCyFE1lNOJwA6FEzbOSRhCyFECijDQDld0sMWQohcoFxudFB62EIIkfWU0yUlESGEyAWxHraURIQQIuspl/SwhRAiJ0gPWwghckRslIj0sIUQIuspp4uo9LCFECL7KZdbathCCJELlEsmzgghRE6QiTNCCJEjZGq6EELkiFgNWxK2EEJkPeVyQSQCETMtx5eELYQQKWK4YmtiEw6n5/hpOaoQQgxBfZsYqDQN7bMP9I1KKQO4G5gKtGmtF6QsKiGEyEEq0cMOpeX4A07YwMeBOq31TakKRgghcplyuWJ/pylhJ1MSuQSYqJRaopT6fKoCEkKIXJXNPexyYBnwTeBlpdTTWuumQ19UU1OTuL1gwQIWLlx4TAcPBALU1dUlEV7ukTYPfkOtvTC02mxrbMIDhLq709LmoyZspVQF8OghDzcCHcCrWmtTKbUCGA+8L2HX1tYOKLC6ujomTJgwoPfmKmnz4DfU2gtDq81BwuwDXAacPMA2NzQ09PvcURO21roRmHPo40qpG4HJwEZgEvCrAUUnhBCDhHLGatiEsq+G/TvgE0qp14A3tNZ7UhSTEELkJOWOD+vLthq21roLuDyFsQghRE5L9LCzcJSIEEKIAyQmzkjCFkKI7JYY1peFNWwhhBAHUE5n7G/pYQshRHZThhGrY4fTs5aIJGwhhEihWMKW1fqEECLrKZcbJT1sIYTIfsolPWwhhMgJyuVO23rYkrCFECKFpIYthBA5QmrYQgiRI6SHLYQQOUJq2EIIkSNio0RkpqMQQmQ95XJn5RZhQgghDlFw5efYN/nMtBxbErYQQqSQvXQ4uqw1LceWkogQQuQISdhCCJEjJGELIUSOkIQthBA5ImsT9mOPPWZ1CBknbR78hlp7QdqcSlmbsB9//HGrQ8g4afPgN9TaC9LmVMrahC2EEOJgSmudtoMvWrQofQcXQohBat68eepwj6c1YQshhEgdKYkIIUSOkIQthBA5IisStlLqXqXUUqXUzw95/I9KqdeVUkuUUp+0Kr50OEKbi5VSjymlFiulvmtVfOlwhDY/Gv+MVyil1loVXzococ0LlFJvxH++P2xVfKl2hPbOV0qtVEq9opSqtiq+dFBKjVRKrVZKBZRS9kOeq1FKLVNKvaaUmpzsuSxP2EqpaUCe1no24FRKzTjkJVdpredorf9qQXhpcZQ2fx+4VWs9V2v9Y2siTL0jtVlrfaXWeg5wN/CMRSGm3FE+568Bc+J/vp756FLvKO29FZgHfBL4gRXxpdF+Ym1beZjnbgc+ASyM306K5QkbmAW8FL/9MnDguoQa+JNS6l9KqaqMR5Y+R2pzDfCdeE8kPWs0WuNIbe7zUeCJjEWUfkdq81bAB+QBnRmOK12O+BlrrXu01nuBcZkOLJ201gGtdVs/TxdprXdrreuBYcmeKxsS9jDe+4Ht4OBG/bfW+izgLuCnmQ4sjY7U5rOAO4ErgXsyHFc6HanNKKUcwCSt9epMB5ZGR2rzk8AaYC3wywzHlS5H+4zL4+WQCZkOzEIH5tjDDtU7HtmwHnYHUBC/XQC09z2htd4f/3uZUup/LIgtXfptM7BJa10HoJSKZjqwNDpSmyFWGliSwXgy4UhtvhU4NX77WeDFDMaVLkdq7zeBR4GdwGsZjstKB46bTvr3ORt62CuI1X8AzueAOpBSqiD+9ym8/xc8l/XbZmCTUmqEUspHdvyHmipHajPEyiFPZjSi9DtSm4NAL9ADODMcV7r0216t9Qqt9XnAj4E6C2Kzyn6l1Cil1EhSUPqyPGHHvwIHlFJLgQiw64DREX9RSi0DfgvcbFWMqXaUNn8feARYDPzIohBT7khtVkopYvXOZRaGmHJH+ZwfINbTXA48aFGIKXWUz/i7SqlXiJX7fmhhmCmnlHIopV4GpgAvKKXOPeT3+W/A48S+VSV3LpnpKIQQucHyHrYQQohjIwlbCCFyhCRsIYTIEZKwhRAiR0jCFkKIHCEJWwghcoQkbCGEyBGSsIUQIkf8f+nP6a+xOcWfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.axvline([0.7], color='k')\n",
    "plt.plot(masses, likes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(-2.0416e+36, device='cuda:0', dtype=torch.float64)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l({'mass ratio': 0.2101})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(1.7698e+35, device='cuda:0', dtype=torch.float64)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l({'mass ratio': 0.92})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(1.0497e+35, device='cuda:0', dtype=torch.float64)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l({'mass ratio': 1.0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 1000\n",
    "\n",
    "times = np.linspace(-2, 2, 1000)\n",
    "q = 0.9#masses[np.argmax(likes)]\n",
    "\n",
    "points = np.ones((1000,2))\n",
    "points[:,0] = times\n",
    "points[:,1] = points[:,1] * np.log(q) * 100\n",
    "\n",
    "test_x = torch.tensor(points).float()\n",
    "test_x = test_x.cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "invalid argument 0: Sizes of tensors must match except in dimension 0. Got 2 and 8 in dimension 1 at /pytorch/aten/src/THC/generic/THCTensorMath.cu:71",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-45-b245066f9339>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m#%%timeit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mno_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgpytorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfast_pred_var\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0mf_preds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtest_x\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0my_preds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlikelihood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf_preds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.virtualenvs/gaston/sandbox/local/lib/python3.7/site-packages/gpytorch/models/exact_gp.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    307\u001b[0m                     \u001b[0mtrain_input\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtrain_input\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mbatch_shape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mtrain_input\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    308\u001b[0m                     \u001b[0minput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mbatch_shape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 309\u001b[0;31m                 \u001b[0mfull_inputs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrain_input\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    310\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    311\u001b[0m             \u001b[0;31m# Get the joint distribution for training/test data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: invalid argument 0: Sizes of tensors must match except in dimension 0. Got 2 and 8 in dimension 1 at /pytorch/aten/src/THC/generic/THCTensorMath.cu:71"
     ]
    }
   ],
   "source": [
    "#%%timeit\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    f_preds = model(test_x)\n",
    "    y_preds = likelihood(f_preds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with torch.no_grad():\n",
    "    lower, upper = y_preds.confidence_region()\n",
    "    f, ax = plt.subplots(1,1, dpi=300)\n",
    "    ax.plot(test_x[:,0].cpu().numpy(), f_preds.mean.cpu().numpy()/1e21)\n",
    "    #samples = y_preds.sample_n(100)\n",
    "    #for sample in samples:\n",
    "    #    ax.plot(test_x[:,0].numpy(), sample.numpy()[:,0], alpha=0.01, color='r')\n",
    "    #    ax.plot(test_x[:,0].numpy(), sample.numpy()[:,1], alpha=0.01, color='b')\n",
    "    #ax.plot(training_x[:,0].cpu(), training_y.cpu()/1e21, '+')\n",
    "    ax.fill_between(test_x[:,0].cpu().numpy(), lower.cpu().numpy()/1e21, upper.cpu().numpy()/1e21, alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generator.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We use SGD here, rather than Adam. Emperically, we find that SGD is better for variational regression\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},\n",
    "], lr=0.1)\n",
    "\n",
    "# Our loss object. We're using the VariationalELBO\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "\n",
    "epochs_iter = tqdm.tqdm_notebook(range(training_iterations), desc=\"Epoch\")\n",
    "for i in epochs_iter:\n",
    "    optimizer.zero_grad()\n",
    "    # Output from model\n",
    "    output = model(training_x)\n",
    "    # Calc loss and backprop gradients\n",
    "    loss = -mll(output, training_y).cuda()\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (sandbox)",
   "language": "python",
   "name": "venv-sandbox"
  },
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}