{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cox Processes (w/ Pyro/QPyTorch Low-Level Interface)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This is a more in-depth example that shows off the power of the Pyro/QPyTorch low-level interface.\n",
    "\n",
    "A Cox process is an inhomogeneous Poisson process, where the arrival rate of events (the intensity function) varies with time. We will use a QP to model this latent intensity function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import qpytorch\n",
    "import pyro\n",
    "import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create sample training set\n",
    "\n",
    "Here we're going to construct an artificial training set that follows a non-homogeneous intensity function.\n",
    "\n",
    "#### Intensity function\n",
    "\n",
    "We're going to assume points arrive following a sinusoidal itensity function.\n",
    "This will give us points in time where there are lots of arrivals, and points in time when there are few."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "intensity_max = 50\n",
    "true_intensity_function = lambda times: torch.cos(times * 2 * math.pi).add(1).mul(intensity_max / 2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determine how many arrivals there are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of sampled arrivals: 53\n"
     ]
    }
   ],
   "source": [
    "max_time = 2\n",
    "\n",
    "times = torch.linspace(0, max_time, 128)\n",
    "num_samples = int(pyro.distributions.Poisson(true_intensity_function(times).mean() * max_time).sample().item())\n",
    "print(f\"Number of sampled arrivals: {num_samples}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determine when the arrivals occur\n",
    "\n",
    "We'll sample the arrivals using rejection sampling. See https://en.wikipedia.org/wiki/Poisson_point_process#Inhomogeneous_case for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_prob_accept(val):\n",
    "    intensities = true_intensity_function(val)\n",
    "    res = torch.log(intensities / (true_intensity_function(times).mean() * max_time))\n",
    "    return res\n",
    "\n",
    "arrival_times = pyro.distributions.Rejector(\n",
    "    propose=pyro.distributions.Uniform(times.min(), times.max()),\n",
    "    log_prob_accept=log_prob_accept,\n",
    "    log_scale=0.\n",
    ")(torch.Size([num_samples]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The result\n",
    "\n",
    "Here's a plot of the intensity function and the sampled arrivals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAGwCAYAAAC6ty9tAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfp9JREFUeJzt3XdcU/f6B/DPSchgJWEPZSqKKAiiIO49am2tdtiqHbd2eLV11N5b76/Wbrttba12WO1ydGmHrYu6xYWgqIiILJmyZ0hIzu+PkFQEhEDCOUme9+vFqzU5OXlOzkny5DueL8OyLAtCCCGEECsg4DoAQgghhBBTocSGEEIIIVaDEhtCCCGEWA1KbAghhBBiNSixIYQQQojVoMSGEEIIIVaDEhtCCCGEWA07rgPoblqtFvn5+XB2dgbDMFyHQwghhJAOYFkW1dXV8PX1hUDQdruMzSU2+fn58PPz4zoMQgghhHRCbm4uevbs2eb9NpfYODs7A9C9MDKZjONoCCGEENIRVVVV8PPzM3yPt8XmEht995NMJqPEhhBCCLEw7Q0jocHDhBBCCLEalNgQQgghxGpQYkMIIYQQq0GJDSGEEEKsBiU2hBBCCLEalNgQQgghxGpQYkMIIYQQq0GJDSGEEEKsBiU2hBBCCLEalNgQQgghxGrwKrF5+eWXwTBMs7/Q0FDD/UqlEgsXLoSbmxucnJwwa9YsFBUVcRgxIYQQQviEV4kNAPTv3x8FBQWGv6NHjxruW7p0KX7//Xf8+OOPOHToEPLz8zFz5kwOoyWEEEIIn/BuEUw7Ozt4e3u3uL2yshIbN27Eli1bMG7cOADApk2b0K9fP5w4cQJDhw7t7lBbUDVqIRIy7S7QRQghhFgjjZYFAAgF3H0P8i6xSU9Ph6+vL6RSKeLi4rB69Wr4+/sjMTERarUaEyZMMGwbGhoKf39/JCQktJnYNDQ0oKGhwfDvqqoqs8W+6VgmvjiSicEBLhgc6IJxoZ4I9nAy2/MRYutYlsXRqyX44kgmLuZVGm7//olYhHrLOIyMENuz+Vgm3t97BZ/Ni8aw3u6cxcGrxCY2NhabN29G3759UVBQgFdeeQUjR47EhQsXUFhYCLFYDIVC0ewxXl5eKCwsbHOfq1evxiuvvGLmyHXO5pSjpKYBuy8WYvfFQrz112UsHh+CBWN6wU7Iu14/QiyWVsvij5QCfHYoAxfzW/5YadSwhv8/frUEAgGDocFu3RkiIVavWqmGxE4IsZ3u+03hIEZ1QyMSs8s5TWwYlmXZ9jfjRkVFBQICAvDBBx/A3t4ejz32WLPWFwCIiYnB2LFj8fbbb7e6j9ZabPz8/FBZWQmZzLS/6JRqDVLyKnEmqxyHr9xAwrVSAECknwIf3D+QWm8IMZFvErLw0q8XAQD2IiFmx/hhZlRPSES6D1h/VwdIRUKwLIupHx3B5cJqPD+5L/49phd1FRNiAieuleK5H87h7khf/GeKbpJPZZ0aueV1CPV2NsuP+aqqKsjl8na/v3nVYnMrhUKBPn364OrVq5g4cSJUKhUqKiqatdoUFRW1OiZHTyKRQCKRdEO0gFQkxJBAVwwJdMXTo4OxMzkPL/16Ecm5FXjrr8v4/OHB3RIHIdbugSF+OJh2A+E95Hh0WCBcHMWtbqdUazGwpwKXC6vx7p40lNWq8H939IOAw/5/QizdugNX8d7eNLAs8GdKAZ4dHwKpSAi5gwhyBznX4fFvVtTNampqkJGRAR8fH0RHR0MkEiE+Pt5wf1paGnJychAXF8dhlK1jGAb3RPXEniWjMKW/N169ewDXIRFi0aqUamibBiZK7ITY+MhgLJ3Yp82kBgDsxUK8fW8EVt4ZBgDYeDQTy388B7VG2y0xE2Jttp7Kwbt7dEnN7CF++OPZkZCKhFyH1QyvuqKWL1+O6dOnIyAgAPn5+Vi1ahWSk5Nx6dIleHh4YMGCBfjzzz+xefNmyGQyPPPMMwCA48ePd/g5OtqUZU4NjRpI7Ph1IRDCZ3WqRsz89DiGBrth1fSwTnUn/XL2Op7/6Tw0WhZ3Rvjg4wejqFuKECMcTCvG41+fgUbL4tnxIVg2sU+3Pn9Hv7951WJz/fp1PPjgg+jbty/uv/9+uLm54cSJE/Dw8AAArFmzBnfeeSdmzZqFUaNGwdvbG7/88gvHURvnm4QsTP3wCEprGtrfmBACAHj190u4XFiNP84XoLRW1al9zBzUE5/Pi4ZIyOCP8wXYeirXxFESYr0u5ldi4fdnodGymDmoB5ZOCOE6pDbxqsWmO3DZYlPb0IgpHx1Gblk9BvkrsO3JOMNockJI634/l49ntiaBYYDvH4/t8myLL49cw95LRfhodiR85PYmipIQ66Zv8YwNcsXmx2I4+e7q6Pc3JTbd7GpxDWatP47KejWen9wXC8f27vYYCLEUuWV1uOOjI6huaMSisb2xfHLfLu9Tq2XBgtsCYoRYouNXS9C/hxxyexEnz2+RXVG2oLenE16+SzeQ8aP4dGSV1HIcESH8pNZo8czWJFQ3NCI6wAVLTNT0LRAwzZKai/mVt9maEKI3rLc7Z0mNMSix4cCMyB4YGeIOVaMW/7czBTbWaEZIh6w/mIHk3ArIpHb4aHakyetisCyLl3+7iGlrj+Lvy7SYLiG3qmloxJPfnLG45J8SGw4wDIPXZwyAxE6AY1dLsSMpj+uQCOEdT2cJXB3FeG3GAPR0cTD5/hmGgUioa7l5fVcqTQEn5Bbv703D3ktFeGZLkmENKEtAiQ1HAtwcsXhCCMRCAcrr1FyHQwjvzI7xx6Hnx+Cugb5me45nxofA1VGMazdqse00zZIiRO/89Qp8fTwLAPDyXf0takwaJTYcemJkMPYsHYXHRwRxHQohvOQsFZm11oxMKsLi8bqxOx/uu4JqJf3IIIRlWby48wK0LHB3pC9G9fHgOiSjUGLDIZFQgCB3R67DIIRXXtyZgj9TCrpt7NlDsf4IcndEaa0KGw5ldMtzEsJn+1OLcf56JRzEQrw4LYzrcIxGiQ1PJGaX45uELK7DIIRTR9Jv4LsTOVi8LQl5FfXd8pwioQAvTNUt4vflkUwUVHbP8xLCRyzL4qP4KwCAh+MC4eHcPWstmhKvF8G0FWmF1Zi1/jjsBAzGhXqaZaAkIXyn1bJ4Y1cqAGDe0MBufR9MCvNCTKArrpXUIPNGLRXuIzbr0JUbuJBXBQexEE+MtMxhEpTY8EBfb2cM7+2GY1dL8enBDLx5TzjXIRHS7fanFuFyYTWcpXZ4Zlz3Fq5kGAbv3z8QCgcRnKX8r9NBiLkM6+WOt2aGo6ahEW5OltdaA1BXFG8sHq9bTOzHM7nd1gRPCJ98eSQTADAnNuC2K3abi5+rAyU1xOaJ7QSYHeOP+SODuQ6l0yix4YmYIFfEBbtBrWGx/uBVrsMhpFudy63AqawyiIQMHh0WyGksWi2L/ZeKUFlPM6SI7WBZ1mpqOVFiwyOLm0rG/3D6OgorlRxHQ0j3+eLINQDA9IG+8JZLOY1lwfeJmP/NGWw9lcNpHIR0pyPpJRj3/kH8cvY616F0GSU2PDI02A1DAl2g0mixhT5UiQ25J6oHYgJdMX8E983fE/p5AQA2H8uCqtE6fsES0p7Nx7OQW1aPlDzLWj6hNZTY8MzDcYEIcHNADwW3v1oJ6U7j+3nhh6fjEObb9oq93eWuSF94OEtQWKXErpR8rsMhxOxyy+pwIK0YgO47yNJRYsMzd4T74MBzY/DAEH+uQyHEJknshIZxPp8fzqRFaonV23IqBywLjAxxt4qisZTY8IxQwEBgQWtyENIVW0/l4KP96SirVXEdSjNzYv1hLxIitaAKxzNKuQ6HELNpaNRge9M6aXOHBnAcjWlQYsNTSrUGO5KuI7OklutQCDELlmWx4VAG1uy/gn2XCrkOpxmFgxj3RvcEAGw5SePdiPX6K6UQZbUq+MilGB/qyXU4JkGJDU+t+CUFS7efM6yuSoi1OZVZhuzSOjiKhZhuxhW8O2t2jB8AILusFhotdUcR6/TdiWwAwIMx/rATWkdKYB1HYYXuieoBAPg58TrqVI0cR0OI6f2YqJtWOn2gLxzE/CuC3t9Xjj+fHYnfF42AkLqHiZVaNqkPpoX7YPYQP65DMRlKbHhqRG93BLg5oLqhEb8m08wMYl1qGhqx63wBAOC+wfz9QA3zlYFhKKkh1mtYL3esmzMInjLrmYlLiQ1PCQQM5sbqBnJRoTBibXadz0e9WoNgD0cM8ldwHU676lSNvBvgTAhpHSU2PDZzUA8IBQzOX69Exo0arsMhxGR+OKPrhrp/sB/vW0S2nsrBkNf3Y218OtehEGIy8alFWP1nKq4UVXMdislRYsNjbk4SjApxBwD8mpTHcTSEmIZao0WgmyOcJXaY2TSWjM98FfaoVWmwIykPSrWG63AIMYnvTmTjs8PX8Mc56xvqQIkNz81o+uDPKq3jOBJCTEMkFOD9+wfi9IsTLKJff0Rvd/jIpaisV2N/ahHX4RDSZaU1DTicXgIAuNsCflwYixIbnpsU5o1Dz4/B2gejuA6FEJOSioRch9AhQgGDWYN0NW1+PGP5CwQSsiulABoti/AecvTycOI6HJOjxIbn7MVCBLhZfolrQgAg40YNLuVXWdwyBfcN1iU2h9NvoLhayXE0hHTNjqahDXdH8q9+lClQYmNBymtVUGtotWFiuTYczMAda49gzb4rXIdilAA3R0T6KcCywO4L/KqSTIgxsktrkZRTAQED3MXDwpimQImNhfjvT+cR8+Z+HG3qFyXE0qgatdh7STdGJa6XO8fRGO/OCB8AwB9N9XcIsUT6umjDe7tbxBi3zuBfuU/SKnuxEGoNi53JeRhrJet5ENtyLKMElfVquDtJEBPkynU4Rrsj3AcNjVrcEe7DdSiEdJqdkIGroxh3R1rfoGE9SmwsxN2Rvth8PAt7LxahtqERjhI6dcSy/NnU0nFHuLdFLlHgq7DHwrG9uQ6DkC7595jeeGJkMCxsmJtRqCvKQkT6KRDo5oB6tYamnBKLo2rUYs9F3dgUavEghFsioQBiO+v9+rfeI7MyDMMYvhD0XxCEWIpjV0tQpWyEh7MEQwItrxvqZr+fy8dT355BcRXNjiKWg2VZnMutgNYGVqqnxMaCTO7vDQA4mHaDKqASi7KvqZVx6gDL7Ia62aZjmdhzsQh/0ewoYkFSC6px97pjGPPeQWisPLmhxMaCRPSUw0cuRZ1KQ7OjiEV5eXp/bH5sCOYNDeA6lC6bFqGbIruLZkcRC6Jv6e/j5WzxPy7aQyNQLQjDMPj3mF7QssBAPwXX4RDSYWI7Acb0tY7ZfHeEe+O1Py7hdHYZiqqU8LLSKbPEuugTm8n9vTiOxPyoxcbCzIsLxCPDAuHhLOE6FEJsko/cHtEBLmBZ4M8UarUh/JdTWofLhdUQChhM6EeJDSGEdFqjRosZ647htT8uoVqp5jock5k6QDfebd8lmqFI+E/fWhMT6AoXRzHH0ZgfJTYWqLSmAVtO5mDrqRyuQyHkts7mVCA5twI/JV6HvYUsetkRE8N0v3pPZZahyooSNmKdbKkbCqAxNhYpMbsc/9uRgh4Ke8we4geGse6BYMRyxTfNhhrb1wN2Quv5HRXg5og+Xk6wF9uhuEoJmVTEdUiEtOpGdQMSc8oBAJOaZtZaO0psLNCoPh6wFwmRV1GPi/lVGNBDznVIhLRKX0xyvBX26/+2aASkVtQKRayTq6MYPz4Vh+TcCvgq7LkOp1tYz08oGyIVCTG6jwcAKtZH+CurpBYZN2phJ2Awqul6tSaU1BBLIBQwGBzoivkjg7kOpdtQYmOhJg/Q/QLee5EGLxJ+0rfWDAl0hdzeertqqpRqlNQ0cB0GIaQJJTYWalxfLwgYIK2oGnkV9VyHQ0gL8anFAIDx/ayjfk1r1h24ikGv7sPnh69xHQohLRxNL8HKnRdw4lop16F0K0psLJTcQYRB/i4AgINpxRxHQ0hzLMsixMsJXjKJVdfNCHBzQKOWNQySJoRPdqXk49sT2dhtY8t/UGJjwcb09YBQwOB6ObXYEH5hGAav3j0AJ1aMR6C7I9fhmM2oPh6wEzDIuFGLrJJarsMhxIBlWRy4fAOA7rvCllBiY8HmDg3A2Rcn4r9TQrkOhZBWWXspAplUhNhg3Wrl+6nVhvBIWlE1CquUkIoEGBrsxnU43YoSGwumcBBD7mC9gzKJZWrUaJGYXW71KwjrjQ/VdbXpxxQRwgcH03StNXHBbjY3g48SGyuhatRyHQIhAHTVhmetP45Jaw6BZa0/udGPITqdVYbKeqpCTPhBP/bSWhafNQYlNhbuanE1Zn56DNM/Psp1KIQAAA5d0X2g9veVW31XFAD4uzmgl4cjGrUsEjJKuA6HEFQr1TiTpas2PNYGExuqPGzhPJylOHe9Ehoti9yyOvi5OnAdErFxR9J1X+7WWJSvLU+OCoaqUYuoppmKhHApt6wePgopRAIB/N1s7zuBEhsLJ7cXYZC/AqezynHwyg3MGxrAdUjEhpXVqpCSVwkAGBXiznE03eeBIf5ch0CIQZivDIefH4uKOtvsGqWuKCug70M9RPVsCMeOXi0BywKh3s7wlEm5DocQm8UwDFwcxVyHwQlKbKyAvkbBsaulUKo1HEdDbNmRK7qZGLbUDaWXX1GPb09k4wD9wCAcqlM1Qq2x7ckklNhYgTAfGTydJahXa3A6q4zrcIiNYlnWML5mpA11Q+ntTM7Dyp0X8P2JHK5DITZs8/EsDHp1H9YfzOA6FM5QYmMFGIYxtNrov1gI4cLHD0Xh2XG9MSTQletQut2oEN17MCGjhMovEM4cuVKC6oZGOEpsq3bNzWjwsJWYGOaN2gaNYf0oQrobwzAYEuhqk0kNoGs5dXMUo7RWhaSccsTaWLVXwr16lQaJ2bpp3iN6216rqR5vW2zeeustMAyDJUuWGG5TKpVYuHAh3Nzc4OTkhFmzZqGoiMqYA8DEMC+smzMIUwZ4cx0KITZJIGAwoqkL7nD6DY6jIbbodFYZVBotfOVSBFnxGm3t4WVic/r0aXz22WeIiIhodvvSpUvx+++/48cff8ShQ4eQn5+PmTNnchQlIURPqdZg1a8XsPdioc0spdAafXcUdQkTLhy7qrvuhvV2t4nimG3hXWJTU1ODOXPm4IsvvoCLyz/dKpWVldi4cSM++OADjBs3DtHR0di0aROOHz+OEydOtLm/hoYGVFVVNfuzVizLIuNGDY5T9VPSzc5klePrhGys/PUCBLb7eWoYNJ2SV4myWhXH0RBbc6zps9+Wu6EAHiY2CxcuxLRp0zBhwoRmtycmJkKtVje7PTQ0FP7+/khISGhzf6tXr4ZcLjf8+fn5mS12rh1JL8H49w/hvz+f5zoUYmOONHW9jAzxsOlfip4yKUK9nSFgGFzMr+Q6HGJDympVuJiv++E+rLdtj+/i1eDhbdu24ezZszh9+nSL+woLCyEWi6FQKJrd7uXlhcLCwjb3uWLFCixbtszw76qqKqtNbqIDXGAnYJBbVo+c0jqbLKVNuHHYBpdRaMvaB6PgLZdCJhVxHQqxMUvG90F2WS08nW27OCZvEpvc3FwsXrwY+/btg1RqupMikUggkUhMtj8+c5TYYZC/C05lleFYRgn83ajMOzG/0poGpBY0/VLsZdu/FAGgj5cz1yEQG+TqKMbiCSFch8ELvOmKSkxMRHFxMQYNGgQ7OzvY2dnh0KFDWLt2Lezs7ODl5QWVSoWKiopmjysqKoK3N80E0hve1Ld69CqNsyHd42SmrihkXy9nuDvZxo8IQgh/8SaxGT9+PFJSUpCcnGz4Gzx4MObMmWP4f5FIhPj4eMNj0tLSkJOTg7i4OA4j55fhTX2rx6+WQGvDs1NI99EPVo+j1hqDbadycM+nx/BT4nWuQyE2oLhaib9SClBRRwPWAR51RTk7O2PAgAHNbnN0dISbm5vh9scffxzLli2Dq6srZDIZnnnmGcTFxWHo0KFchMxLA/0UcBQLUV6nRmphFfr7yrkOiVi5rJI6AJTY3Cy/oh5JORUIcHXAvdE9uQ6HWLm/U4vxwi8pGBzggp8WDOM6HM7xJrHpiDVr1kAgEGDWrFloaGjA5MmT8emnn3IdFq+IhAIMDXZD/OViHLtaQokNMbvv5scit6wObk62uZJwa4b2csPav68i4VopWJa16ZlixPyOZZQC0NWvIQDDsqxN9VdUVVVBLpejsrISMpmM63DM4sS1UtSpGhET5AYniUXlroRYBaVag4hX9kLVqEX8c6PRy8OJ65CIldJqWQx5Yz9Ka1X44ak4xARZ75ImHf3+5s0YG2I6Q4PdMC7Ui5IaQjgiFQkR3bRuW0LTr2lCzCGtqBqltSrYi4SI9FNwHQ4vUGJDCOm06R8fxRPfnMH18jquQ+Ed/dR3SmyIOZ24pru+hgS5QmxHX+kAJTZW60JeJd7ZfRk7kmhWBjGP6+V1SMmrxN+Xi6FwoPE1t9IPpj5xrZRmKBKz0Sc2sVbcBWUsSmysVGJ2OT49mIFfzuZxHQqxUvqWiIiecur2bEVETwV6KOwRG+yK6oZGrsMhVkirZXGqqY7U0GCalahHn0ZWKjZYl72fySqHWqOFSEg5LDEtfWJD1YZbJ7YT4Oh/x9KMKGI2DAP8tmgEEq6VIqInzYDVo8TGSvXxdIaLgwjldWqcv16J6ACX9h9ESAexLIuEpibwuGCaYtoWSmqIOTEMAz9XB/i50rqAN6Of8VZKIGAM0/5OZtLgRWJa2aV1KKhUQiwUUNLcDpZlce1GDTQ0zoaQbkGJjRXT97meuFbGcSTE2uhbayL9FbAXCzmOhr9YlsWkNYcx7v1DuJRfxXU4xIqwLIul25Px5ZFrqFPRGK6bUWJjxWKDdIlNYlYZ1Botx9EQayK3FyEm0BWj+3hwHQqv6bsKAGo5JaZ1tbgGO5Ly8N7eNNgJ6Kv8ZvRqWLFQb2fI7UUQ2wlwvbye63CIFbkj3Ac/PB2HhWN7cx0K7w0J1HUJ62evEGIK+mne0QEuVL/mFjR42IoJBAz+WjwS3jIpBAIaxEgIF/Rj3U5nldG6UcRkTjQlyvqWefIPSvOsnK/CnpIaYlKFlUpU1qm5DsNihPeQQyoSoLxOjavFNVyHQ6wAy7I4SYX52kSJjY1gWRY2tt4pMZOP4q8g8rW9+PxwBtehWASxnQBRfrqZYyepO4qYQMaNWpTUqCCxE2AgrQ/VAiU2NmDlzguIW/03UvIquQ6FWIGTmWVgWSDYnVas7ih9dxSNsyGmoB9fE+WvgFREsxJvRWNsbEBBpRKFVUqcuFaKiJ4KrsMhFuxGdQOu3agFw/wzKJa0b1yoJ+pUjRjdx5PrUIgVuFHdAJGQoWUU2kCJjQ2IDXLF/tQinMosx5OjuI6GWLIzWboWh75ezpA7iDiOxnIM9FNQlwExmaUT++Dp0b2gojIeraLExgYMaWoGT8wug1bL0mBi0mn6MSIxNGCREE7Zi4WwB3VDtYbG2NiA/r4yw6yMjBs0K4N03ilKbDqtXqXB0fQS7L1YyHUoxILRJJD2UWJjA0TCf2ZlnM4q5zgaYqmqlGqkFuqWBYih8TVGS7hWgrkbT2L1X5e5DoVYsJW/XsC0tUewhxLkNlFiYyOGBOoTG5qVQTpHwDBYfU845o8IgqdMynU4Fic6wBUMA2SW1KK4Wsl1OMRCnbxWhov5VaABBW2jMTY2IibIDQN6FKOXhyPXoRAL5SSxw+wYf67DsFhyexFCvWVILajC6cxyTIvw4TokYmHKa1VIbyryGB3gwnE0/EWJjY0YEeKOP0JGch0GITYtNsgVqQVVOJVZSokNMVpitm4oQS8PR7g5STiOhr+oK4oQ0q56lQbfJGThYn4lDV7sAkOhPhrrRjpBP5SABu/fHiU2Nkap1iC7tJbrMIiFOXe9Ai/9ehGPbTrNdSgWbXBT90FaYRWqlbTeFjGOPrEZHECJze1QYmNDjqaXIPzlPXj6u7Nch0IsjL4JfHCgC61O3QWeMin8XO2hZYGknAquwyEWRKnWGJbFoarft0djbGxIHy8nqDUsLhdWoUqphkxKlWNJx+grDkfTL8UuW31PBFwcdQOJCemoKqUad4T7IKu0Dn6u9lyHw2uU2NgQT5kUAW4OyC6tw9nscozpS+vWkPZptew/LTY0E6PLRoS4cx0CsUCezlJ8NDuK6zAsAnVF2Rh93yzVsyEddfVGDaqUjbAXCRHmS60MhBB+o8TGxvxTqI9mZZCOOdN0rQz0k0MkpI8MU9iZlIfnfjiHtMJqrkMhFkCjZZFeVA2tlmYkdgR9StkY/YKYybkVaGjUcBwNsQT/dEPR+BpT2ZGUh5/PXsfxjBKuQyEWIK2wGhPXHMaodw9QuYUOoMTGxgS7O8LVUQxVoxYX86u4DodYgFfu7o/vHo/FzEE9uA7FaujHKp3JppZT0r7EbN3QgSB3R5qV2AE0eNjGMAyDJ0cFQyQUwFdOI+tJ+5wkdjTg1cSim7qEE7PKwbIsfVmR29K3mtIyCh1DiY0Nenp0L65DIMSmRfopIBQwKKxSIq+iHj1dHLgOifBYYg4lNsagrihCSJu+TcjCG7su4UJTYTBiGg5iO/RvmmGWSN1R5DaKq5TILasHw+gSYtI+SmxsVGZJLX5KvI7SmgauQyE89ktSHr44kkmzd8xA/+v7DM1QJLdxtqm1pq+XM5ypqGqHUGJjoxZ+fxbLfzyHk5lUz4a0TqnWGFpqBgdSE7ip6WeZldbSjwvSNhpfYzwaY2OjBgUocKmgCmezy3FHuA/X4RAeOn+9EmoNC3cnCfxdaQyIqY0L9cS5lyZB7kC/wknbpgzwgcROiNhgKrfQUZTY2KhB/i747kSOYVAaIbe6eRkFmrVjevZiIezFQq7DIDwXHeBCrTVGoq4oGzXIX/dGuZhXRYX6SKtuXtGbEEIsBSU2NirAzQFujmKoNFpcyKNCfaQ5lmWR1NSaF+VPiY25JGaX4f4NCVjwXSLXoRAeOp1VhvjUIlTUqbgOxaJQYmOjGIYxfGElUXcUuUVprQoMw0AsFGBAD1r40lzsBAKcyirD8YxSKpVPWth4JBOPf30G207nch2KRaExNjZsUIAC+1OLDNMJCdFzd5Lg9P+NR3F1AyR2NA7EXPr5yCCxE6CyXo3MkloEezhxHRLhCZZlqTBfJ1FiY8OmhfsgxNMZUf4KrkMhPMQwDLxkUq7DsGpiOwHCe8hxJrscZ3MqKLEhBtfL63GjugEiIYPwHnKuw7Eo1BVlwwLcHDExzAvuThKuQyHEZul/WFCXMLmZfvB+f185pCJqNTUGJTaEkGbqVI2IWx2Pp749A6WaZsyZ2yDDWLcKbgMhvEKF+TqPuqJs3JWiavxxvgDeMikeivXnOhzCA+dyK1FQqQQA+qXYDfSD+C8XVqG2oRGOEvpYJv8kNoNoVqLRqMXGxl3Kr8La+HT8cIZG3RMd/WBy+kDtHt5yKfr7yjC2rycq69Vch0N4oF6lQVqRbn22QQEKboOxQPTTwMYZCvXlV0Kp1tAvdHJT/RoFt4HYkF3PjuQ6BMIjEjsB9iwZiQt5VfCR23MdjsWhFhsb5+dqD3cnMdQaFhfzK7kOh3CMZVmcbRrrMYj69gnhhEDAoLenM2ZE9eA6FItEiY2Nu7lQ39nsCm6DIZzLLq1DWa0KYqEA/X2pMF93K6xUUqE+QrqIEhti6I7SD1Yjtks/vmZADxkV5utGGi2LUe8cwNDV8cgtq+c6HMKxlTsvYNOxTFQracxVZ1BiQzCoaSzF2Zxy+rVo4xzEQgzyV2BosBvXodgUoYCBq6MYAJCUSz8wbFlRlRLfnsjGa39cgoBhuA7HItHgYYKIngrYCRjUqzQoq1XBjQr22awpA3wwZYAP12HYpCh/BZJzK5CUU4G7I2lsha3S1zPq4+VMU/87qUuvmlqtRmFhIerq6uDh4QFXV1dTxUW6kb1YiP3LRsPP1QFCAf1CIIQLUf4u2HQsi9Zus3H6Fjualdh5RndFVVdXY/369Rg9ejRkMhkCAwPRr18/eHh4ICAgAE888QROnz5tjliJGQW6O1JSY+Mq69WoaWjkOgybpe8SvpRfRRWfbVhyU4tNlB/NSuwsoxKbDz74AIGBgdi0aRMmTJiAnTt3Ijk5GVeuXEFCQgJWrVqFxsZGTJo0CVOmTEF6erq54iaEmNiWkzmIeHkPXv39Eteh2KQeCnt4OEvQqGVxIY9KL9iiRo0W56/rzn0ktdh0mlGJzenTp3H48GGcOnUKK1euxOTJkxEeHo7evXsjJiYG//rXv7Bp0yYUFhZixowZOHLkiFHBrF+/HhEREZDJZJDJZIiLi8Nff/1luF+pVGLhwoVwc3ODk5MTZs2ahaKiIqOeg7ROqdZg6fZkjHv/IP1atFHJueXQsoCXjMZYcYFhGET5KQCAuqNs1JWiGtSrNXCS2KEXrfTeaUaNsdm6dWuHtpNIJHj66aeNDqZnz5546623EBISApZl8fXXX+Puu+9GUlIS+vfvj6VLl2LXrl348ccfIZfLsWjRIsycORPHjh0z+rlIcxI7AY5eLcGN6gak5FViSCCNl7I1ybkVAIDIpi9X0v2mhnvDV2FvqC1FbEtmSS2EAgYRPeU0NKALGJbn83tdXV3x7rvv4t5774WHhwe2bNmCe++9FwBw+fJl9OvXDwkJCRg6dGirj29oaEBDQ4Ph31VVVfDz80NlZSVkMipAdrOnvj2DPReLsGJqKJ4a3YvrcEg3KqisR9zqvyEUMEh5eRIcxDQbgxAu1KkaUVarQk8XB65D4Z2qqirI5fJ2v79NUsdGo9Hg4sWL2L59O1auXIl77rnHJPvctm0bamtrERcXh8TERKjVakyYMMGwTWhoKPz9/ZGQkNDmflavXg25XG748/Pz63Js1kpfqI+awW1P8k1TTCmpIYQ7DmI7Smq6yOhPsGvXriElJQUXLlww/KWnp0OtVkMsFqNfv34IDw/vdEApKSmIi4uDUqmEk5MTduzYgbCwMCQnJ0MsFkOhUDTb3svLC4WFhW3ub8WKFVi2bJnh3/oWG9KSfm2gszkVYFkWDBWHshnUDcUfdapGnL9eCXcnMXp7OnMdDiEWx6jEZu7cudi6dSsYhoGDgwNqa2sxbdo0vPTSSwgPD0dISAiEwq6VYe/bty+Sk5NRWVmJn376CY888ggOHTrU6f1JJBJIJDQYsiPCe8hhJ2Bwo7oB18vr4edKvxpsRVJTYhNFiQ3nVv95Gd+eyMb8EUF48c4wrsMh3eTEtVK89ddlTOjniUXjQrgOx6IZ1RX1008/Ye3ataipqUF+fj4WLVqEvXv34vTp0wgICOhyUgMAYrEYvXv3RnR0NFavXo2BAwfio48+gre3N1QqFSoqKpptX1RUBG9v7y4/LwGkIiHCmhY+1H/REdswM6oH7onqgcGBNGiVa/rCbPQetC2J2eVIzq3A5cJqrkOxeEYlNkuXLsXDDz8MqVQKJycnfPTRRzh27BgOHDiA/v37Y/fu3SYPUKvVoqGhAdHR0RCJRIiPjzfcl5aWhpycHMTFxZn8eW1VdIALQr2dQZ1QtmV2jD/WPBCJYJpiyjn9jKiUvEqoGrUcR0O6S1LT2EbqDu46o7qiVq9e3eK26OhonDp1CmvXrsUDDzyAadOm4aOPPoKHh4fRwaxYsQJTp06Fv78/qqursWXLFhw8eBB79uyBXC7H448/jmXLlsHV1RUymQzPPPMM4uLi2pwRRYz30p1hNLaGEA4FujnAxUGE8jo1LhVU0RedDWBZ1jDOjZZS6DqTzIpiGAaLFy/GpUuX0NDQgNDQ0E7tp7i4GA8//DD69u2L8ePH4/Tp09izZw8mTpwIAFizZg3uvPNOzJo1C6NGjYK3tzd++eUXUxwCaUJJje05lVmGtMJqaLS8rvxgMxiGMbTaJNEMRZtwvbweJTUqiIQM+vvKuQ7H4pl0XmePHj3w888/Y9euXZ16/MaNG297v1Qqxbp167Bu3bpO7Z90nKpRC42Whb246+OmCL+9uDMFV4pq8MXDgzExzIvrcAh0g7j/vlyMpJwKPDac62iIuelba/r5yCAV0WduVxnVYpOTk9Oh7aZNmwYAyMvLMz4iwrnX/7iEAS/vwQ9ncrkOhZhZtVKN9OIaANS3zydRVFPKplC5BdMyqsVmyJAhmDFjBubPn48hQ4a0uk1lZSV++OEHfPTRR3jyySfx7LPPmiRQ0n2cpSKoGrVIyinHI8MCuQ6HmFHK9Uqw7D8LMBJ+iPRXYNX0MET5u5ikppRGo4FarTZRdMTU3KQMBvo4YIifE5RKJdfhcEYkEplkdrVRic2lS5fwxhtvYOLEiZBKpYiOjoavry+kUinKy8tx6dIlXLx4EYMGDcI777yDO+64o8sBku6nX1U2maabWj39lGJaSZhfnCR2eGx4UJf3w7IsCgsLW5TJIPwytgcwtoc3WLYemZmZXIfDKYVCAW9v7y4l80YlNm5ubvjggw/wxhtvYNeuXTh69Ciys7NRX18Pd3d3zJkzB5MnT8aAAQM6HRDhXmRPBQAgq7QO5bUquDiKuQ2ImE0yFeazavqkxtPTEw4ODjQ5gPAWy7Koq6tDcXExAMDHx6fT++rU4GF7e3vce++9hsUoiXWRO4gQ7OGIazdqkZxbgbGhnlyHRMzg5imm1LfPP+W1KuxLLUJlnRpPjAo2+vEajcaQ1Li5uZkhQmIKGi0LAUMzUgFdbgHoZkh7enp2ulvKJNO9ifXRf9FR9VPrlVdRjxvVDbATMBjQg6aY8k1RtRL/+ek8Ptx/pVNT8fVjahwcaGkUPsuvqEdqQTXK61Rch8IL+uu1K2PCKLEhrdJ3TdA4G+vl5ijBpseGYNVd/WmKKQ+FeDrDUSxErUqD9OLOl9mnlgB+q1Np0KjVwk5A5wkwzfVq0jo2xHoMCXLF5P5eGN7bnetQiJnYi4UY25e6GflKKGAQ3lOOE9fKkJxTgVBvGdchERNr1GjR0KgBANjTjwuToRYb0qpQbxk+mzcYD8cFch0KITYr0k9Xz4ZaTq1TvVqX1EjshLAT0texqXT6lXzkkUdw+PBhU8ZCCOkmao0W7+1Jw96LhWjU0EKLfBVFpResWp1Kl9g4UIV3k+p0YlNZWYkJEyYgJCQEb775JlUZtkIsyyK3rA4p1yu5DoWYWFphNT45cBXP/XgOAhqDwVv6sW5XiqpR29DIbTA27oUXXoBEIsFDDz1ksn3qExtausa0Op3Y7Ny5E3l5eViwYAG2b9+OwMBATJ06FT/99BNVuLQSey4WYuQ7B7Bix3muQyEmpl9cMdJPAQENWuQtT5kUvnIptCxwubCK63C6DcMwt/17+eWXuz2mFStW4P3338fWrVtx9erVLu+PZVnUqXTJKrXYmFaXOvU8PDywbNkynDt3DidPnkTv3r0xb948+Pr6YunSpUhPTzdVnIQD4U2F+i4XVKO+6ZcFsQ5JVJjPYnz+8GAkvjgB0QGuXIfSbQoKCgx/H374IWQyWbPbli9f3uIxKpV5p0vL5XI8/vjjEAgESElJ6fL+WBZwdRTDSWJHsxJNzCSjlQoKCrBv3z7s27cPQqEQd9xxB1JSUhAWFoY1a9aY4ikIB3zlUng4S9CoZXEhn7qjrEkyLaVgMQb0kMPNybTreNWpGtv8U6o1Jt/WWN7e3oY/uVwOhmGa3ebk5IQxY8Zg0aJFWLJkCdzd3TF58mQAQGBgID788MNm+4uMjDS08mi1WqxevRpBQUGwt7fHwIED8dNPP3UorsbGRjg4OODChQtGH9OtBAIGPnJ7BHs4UXewiXV6urdarcZvv/2GTZs2Ye/evYiIiMCSJUvw0EMPQSbTTUvcsWMH/vWvf2Hp0qUmC5h0H4ZhEOWnwN5LRUjOqcCQQNv5xWjNKuvUuHajFgAwsKlVjtiWsJf2tHnf2L4e2PRYjOHf0a/tN8zeuVVskCu2PxVn+PeItw+grLZly0nWW9O6EG3bvv76ayxYsADHjh3r8GNWr16N7777Dhs2bEBISAgOHz6MuXPnwsPDA6NHj77tY1988UXU1NSYJLEh5tPpxMbHxwdarRYPPvggTp06hcjIyBbbjB07FgqFogvhEa5F+jclNjQrw2qcu14BAPB3dTB5SwAxjw/3X8GZrHK8e18EfOT2XIfDGyEhIXjnnXc6vH1DQwPefPNN7N+/H3FxuoQsODgYR48exWeffXbbxCYxMREbNmzAtGnTTJLY1DY0QiISwE5A07xNrdOJzZo1a3DfffdBKpW2uY1CobD5lUotnWFphabBpsTypeTpuhVpfSjLse9SES7mVyEppwI+4V1PbC69OrnN+27tFklcOaHD2x7979iuBWak6Ohoo7a/evUq6urqMHHixGa3q1QqREVFtfk4rVaLp556CosWLUJsbCzmzp0LtVoNkUjUqbi1WhbXSmrBsixCvWUQ21FyY0qdTmxGjx4NiaTlrz2WZZGbmwt/f/8uBUb4IaKnAgIGyK9UorhKCU9Z24kssQz/HtMLUwd4Q8sav/4Q4UaknwIX86uQnFuBO8I7v+qxnoO44x/95trWFBwdHVvcJhAIwN5ybetn6tbU1AAAdu3ahR49ejTbprXvM72PP/4YJSUlePXVV5GTkwO1Wo3Lly8jPDy8U3HXqzVgWRZ2AgFEQhpfY2qdvgqDgoJQUFAAT8/mJdnLysoQFBQEjYZm0VgDJ4kd/jslFH6uDnCU0Aoc1oBhGAR7OHEdBjFCpJ8C35/MQXJOBdeh8J6HhwcKCgoM/66qqjL0HISFhUEikSAnJ6fd8TR6eXl5WLlyJbZu3QpHR0eEhIRAIpHgwoULnU5sbi7MR2t5mV6nv6lYlm31hNTU1Ny2e4pYnqdG9+I6BEJsmr4C8fm8Cqg1Woio/H6bxo0bh82bN2P69OlQKBR46aWXIBTqplM7Oztj+fLlWLp0KbRaLUaMGIHKykocO3YMMpkMjzzySIv9Pfvss5g6dSqmTdMNgLazs0O/fv26NM6mnurXmJXRic2yZcsA6H71rVy50rDEOABoNBqcPHmy1YHEhBDu/X25CL+czcPk/t6YPtCX63BIBwW7O8FZaodqZSPSCqsxoIec65B4a8WKFcjMzMSdd94JuVyO1157rdlYz9deew0eHh5YvXo1rl27BoVCgUGDBuF///tfi3398ccf+Pvvv5Gamtrs9vDw8C4lNrSUgnkx7K2dke0YO1Y3OOzQoUOIi4uDWCw23CcWixEYGIjly5cjJCTEtJGaSFVVFeRyOSorKw3T0sntabQsTlwrRXJuBZ4aFUyLtVmwV36/iE3HsvDosEC8fFd/rsMhRpj75UkcvVqC12cMwNyhAe1ur1QqkZmZiaCgIGpF5xG1RovUAl0V6f6+MghpVlQzt7tuO/r9bXSLzYEDBwAAjz32GD766CNKDmzEk9+cQa1Kg3GhnujnQ+fcUhkK89GMKIsT6afAuesVtGaUhdNXcZfaCSmpMZNOj7HZtGmTKeMgPCYUMBjop8DxjFIk5VRQYmOhVI1aXMzX/VKkxMbyLBrXG8sm9qG1vSycVCRED4U9DRo2I6MSm2XLluG1116Do6OjYaxNWz744IMuBUb4JbIpsUnOLcdDsTSV3xKlFlRB1aiFi4MIAW4O7T+A8AqtJ2QdxHYCKoxpZkYlNklJSYZ6AElJSW1uR5mo9dH/wqcKxJZLf+4G+inoPWrhtFqWWm4IaYNRiY1+fM2t/0+sn36xxPTiGlQr1XCWdq7iJuEOja+xfJuPZeLLo5mYPcQPi8bxc4IGaZuqUYtqpRoOYjvY04wos+n0yKX6+nrU1dUZ/p2dnY0PP/wQe/fuNUlghF88naXoobAHywLnr9NK35ZIqdbATsBQYmPBNCxwvbweSVSozyLVNDQir6Ie+ZX1XIdi1To9ePjuu+/GzJkz8fTTT6OiogIxMTEQi8UoKSnBBx98gAULFpgyTsIDkf4K5FXUIzm3AsN7u3MdDjHS+rnRUKo1Ldb3IZbj5i7htoqkEv6qo8J83aLTLTZnz57FyJEjAQA//fQTvL29kZ2djW+++QZr1641WYCEPxaO6Y0/nhmBJ0cFcx0K6SSpSEgL7lmw/r4yiIQMSmtVuF5Ov/otDRXm6x6d/oSrq6uDs7MzAGDv3r2YOXMmBAIBhg4diuzsbJMFSPgjzFeGAT3kVM7dAhlZh5PwlFQkRFhTuYUkGshvUTRaFg3qpsRGROvumVOnv6F69+6NnTt3Ijc3F3v27MGkSZMAAMXFxVS0jxCeeeKbREz/+ChOXCvlOhTSRYbuKBpnAwAIDAzEhx9+yHUY7apXa8ACEAkFEN2m1ZSL4zH1c2ZlZYFhGCQnJ5tsn8bodGLz0ksvYfny5QgMDERsbCzi4uIA6FpvoqKiTBYg4Ze9Fwvxn5/O4dCVG1yHQjqIZVmcyS5DSl4l7KkWisXTz1BMyi3vtufUaFkkZJTi1+Q8JGSUQqM1fwtgbm4u/vWvf8HX1xdisRgBAQFYvHgxSkstMzk39fia69evQywWY8CAAV3e1+nTp/Hkk0+aICp+6HR72L333osRI0agoKAAAwcONNw+fvx43HPPPSYJjvDPsasl+OHMdThJRBjdx4PrcEgHZJXWoaJODbGdgKpGW4EoPxf095Uhys+lW55v94UCvPL7JRRUKg23+cilWDU9DFMG+JjlOa9du4a4uDj06dMHW7duRVBQEC5evIjnn38ef/31F06cOAFXV1ezPHd7NBoNGIaBwMjlEOpNPL5m8+bNuP/++3H48GGcPHkSsbGxt91erVZDJGpepkOlUkEsFsPDw7o+y7s0WMLb2xtRUVHNTnBMTAxCQ0O7HBjhJ/2vxeRu/LVIukZ/rvr7ymjgsBUIdHfErmdH4qXpYWZ/rt0XCrDgu7PNkhoAKKxUYsF3Z7H7QoFZnnfhwoUQi8XYu3cvRo8eDX9/f0ydOhX79+9HXl4e/u///q/Z9tXV1XjwwQfh6OiIHj16YN26dYb7WJbFyy+/DH9/f0gkEvj6+uLZZ5813N/Q0IDly5ejR48ecHR0RGxsLA4ePGi4f/PmzVAoFPjtt98QFhYGiUSCL7/8ElKpFBUVFc3iWLx4McaNG2f499GjRzFy5EjY29tjRGQoPn3zRdhpVYb7i4uLMX36dNjb2yMoKAjff/99h14flmWxadMmzJs3Dw899BA2btzY7H59V9D27dsxevRoSKVSfP/993j00UcxY8YMvPHGG/D19UXfvn0BNO+Keuihh/DAAw80259arYa7uzu++eYbAMDu3bsxYsQIKBQKuLm54c4770RGRkab8ZaXl2POnDnw8PCAvb09QkJCzLosU5c+5eLj4/G///0P8+fPx7/+9a9mf8Q6RTb9SryQryvPT/hPPxaD6tcQY2i0LF75/RJa63TS3/bK75dM3i1VVlaGPXv24N///jfs7e2b3eft7Y05c+Zg+/btzQbEv/vuuxg4cCCSkpLwwgsvYPHixdi3bx8A4Oeff8aaNWvw2WefIT09HTt37kR4eLjhsYsWLUJCQgK2bduG8+fP47777sOUKVOQnp5u2Kaurg5vv/02vvzyS1y8eBFz5syBQqHAzz//bNhGo9Fg+/btmDNnDgAgIyMDU6ZMwaxZs3D+/Hls374dp04cx7Iliw2PefTRR5Gbm4sDBw7gp59+wqeffori4uJ2X6MDBw6grq4OEyZMwNy5c7Ft2zbU1ta22E7/WqSmpmLy5MkAdN/baWlp2LdvH/74448Wj5kzZw5+//131NTUGG7bs2cP6urqDL0xtbW1WLZsGc6cOYP4+HgIBALcc8890Gpb/05YuXIlLl26hL/++gupqalYv3493N3NVzKk011Rr7zyCl599VUMHjwYPj4+VE/BRgS6OUDhIEJFnRqpBVUYSF+WvEcVh62TUq1BUZUSAW6OZtn/qcyyFi01N2MBFFQqcSqzDHG93Ez2vOnp6WBZFv369Wv1/n79+qG8vBw3btyAp6cnAGD48OF44YUXAAB9+vTBsWPHsGbNGkycOBE5OTnw9vbGhAkTIBKJ4O/vj5iYGABATk4ONm3ahJycHPj6+gIAli9fjt27d2PTpk148803AehaLD799NNmwy5mz56NLVu24PHHHwegSxgqKiowa9YsAMDq1asxZ84cLFmyBAAQEhKCtWvXYvTo0Vi/fj1ycnLw119/4dSpUxgyZAgAYOPGjW0e9802btyI2bNnQygUYsCAAQgODsaPP/6IRx99tNl2S5YswcyZM5vd5ujoiC+//BJisbjVfU+ePBmOjo7YsWMH5s2bBwDYsmUL7rrrLsNMaP0x6n311Vfw8PDApUuXWh3zk5OTg6ioKAwePBiAroXInDrdYrNhwwZs3rwZJ0+exM6dO7Fjx45mf8Q6MQxD60ZZEKVag0sFuhW9u2tMBjG/xOwyhL+8B49uOm225yiubjup6cx2xjKmRIF+8srN/05NTQUA3Hfffaivr0dwcDCeeOIJ7NixA42NuoG8KSkp0Gg06NOnD5ycnAx/hw4data1IhaLERER0ew55syZg4MHDyI/Px8A8P3332PatGlQKBQAgHPnzmHz5s2GfTo6OWHy5MnQarXIzMxEamoq7OzsEB0dbdhnaGio4fFtqaiowC+//IK5c+cabps7d26L7igAhkTiZuHh4W0mNQBgZ2eH+++/39AtVltbi19//dXQEgXoks8HH3wQwcHBkMlkhkQlJyen1X0uWLAA27ZtQ2RkJP7zn//g+PHjtz3Grup0i41KpcKwYcNMGQuxEJF+ChxMu4Hk3Ao8wnUw5LaqlY2Y1N8bOaV18HO1b/8BxCL08nCCWsMis6QWFXUqKBza/qLqLE9nqUm366jevXuDYRikpqa2OhElNTUVLi4uHR7w6ufnh7S0NOzfvx/79u3Dv//9b7z77rs4dOgQampqIBQKkZiYCKGw+aBeJycnw//b29u36JUYMmQIevXqhW3btmHBggXYsWMHNm/ebLi/pqYGTz31FJ599llkl9aiXqWBt1wKhYMY/v7+uHLlihGvyj+2bNkCpVLZbLAwy7LQarW4cuUK+vTpY7jd0bFla15rt91qzpw5GD16NIqLi7Fv3z7Y29tjypQphvunT5+OgIAAfPHFF/D19YVWq8WAAQOgUqla3d/UqVORnZ2NP//8E/v27cP48eOxcOFCvPfee8Yceod1usVm/vz52LJliyljIRYiyl/3y7+stvWLmPCHh7ME6x4ahN+fGUHdxVZE4SBGkLvuC8pcLacxQa7wkUvR1lXDQDc7KibItLOT3NzcMHHiRHz66aeor29eXbmwsBDff/89HnjggWbX84kTJ5ptd+LEiWZdOvb29pg+fTrWrl2LgwcPIiEhASkpKYiKioJGo0FxcTF69+7d7M/b27vdWOfMmYPvv/8ev//+OwQCAaZNm2a4b9CgQbh06RJ69eoFr56B8A8KxoB+fdG7d2+IxWKEhoaisbERiYmJhsekpaW1GJB8q40bN+K5555DcnKy4e/cuXMYOXIkvvrqq3Zj7ohhw4bBz88P27dvx/fff4/77rvPMKOqtLQUaWlpePHFFzF+/HhD12B7PDw88Mgjj+C7777Dhx9+iM8//9wksbam0y02SqUSn3/+Ofbv34+IiIgW08g++OCDLgdH+GlosCvOvTQJcgda4ZsQrkT6KZBZUovk3AqM6etp8v0LBQxWTQ/Dgu/OggGaDSLWpxSrpodBKDB9wvzJJ59g2LBhmDx5Ml5//fVm07179OiBN954o9n2x44dwzvvvIMZM2Zg3759+PHHH7Fr1y4AullNGo0GsbGxcHBwwHfffQd7e3sEBATAzc0Nc+bMwcMPP4z3338fUVFRuHHjBuLj4xEREdEsUWnNnDlz8PLLL+ONN97AvffeC4lEYrjvv//9L4YOHYoF/16IcTMehJOjI66dycH+/fvxySefoG/fvpgyZQqeeuoprF+/HnZ2dliyZEmLAdM3S05OxtmzZ/H999+3mH384IMP4tVXX8Xrr79u7MvdqoceeggbNmzAlStXcODAAcPtLi4ucHNzw+effw4fHx/k5OQYxje15aWXXkJ0dDT69++PhoYG/PHHHx0aS9RZnW6xOX/+PCIjIyEQCHDhwgUkJSUZ/riqNki6h8ROSEmNhcgtq6PlFKyUfqybOVf6njLAB+vnDoK3vHl3k7dcivVzB5mtjk1ISAjOnDmD4OBg3H///ejVqxeefPJJjB07FgkJCS1q2Dz33HM4c+YMoqKi8Prrr+ODDz4wzAJSKBT44osvMHz4cERERGD//v34/fff4eamG/C8adMmPPzww3juuefQt29fzJgxA6dPn4a/v3+7cfbu3RsxMTE4f/58szEoABAREYFDhw4h7coVPDbrDtw3eRRWrVplGKSsf25fX1+MHj0aM2fOxJNPPmkYEN2ajRs3IiwsrNWSKvfccw+Ki4vx559/tht3R8yZMweXLl1Cjx49MHz4cMPtAoEA27ZtQ2JiIgYMGIClS5fi3Xffve2+xGIxVqxYgYiICIwaNQpCoRDbtm0zSZytYVgb+9SrqqqCXC5HZWUlLf1ArFppTQOiX98PV0cxjr8wDlKqOmxVzl+vwF2fHIPCQYSklRNbdDUqlUpkZmYiKCgIUmnXxsFotCxOZZahuFoJT2dd95M5WmqsUW5ZHcrrVPB0lrZIEElLt7tuO/r93aU6NkeOHMHcuXMxbNgw5OXlAQC+/fZbHD16tCu7JRbgbE45HvgsAQu+S2x/Y8KJc9crAAAuDiJKaqxQqLeu4GJFnRpZpXVmfS6hgEFcLzfcHdkDcb3cKKkxAq3o3f06ndj8/PPPmDx5Muzt7XH27Fk0NDQAACorKw1z/4n1EgkEOJlZhuMZpdB2w7oxxHj/FOajad7WSGwnwBMjg/C/O0LhJKHVovlIq2UNBQztKbHpNp1ObF5//XVs2LABX3zxRbOBw8OHD8fZs2dNEhzhr1AfZ0jsBKisVyOztGXFS8K9JH1hvqZlMIj1eX5yKJ4c1QsezpL2NybdTiBg0M/HGaHezhAJaTmT7tLpVzotLQ2jRo1qcbtcLm93uhqxfCKhAOE95AD+aRkg/KHVsjjXlNhEUcVhQjjDMAzEdtRa0506ndh4e3vj6tWrLW4/evQogoODuxQUsQxUgZi/rpXUokrZCKlIgL7ezlyHQ8yEZVlkldRiZ1IelGpNq9u0tX4PIXxkiuu10x2zTzzxBBYvXoyvvvoKDMMgPz8fCQkJWL58OVauXNnlwAj/6Qr1ZSKJVvrmHX2yGd5DTk3gVm7W+uMorVXhZ9dhiA74ZzyVWCyGQCBAfn4+PDw8IBaLqUhjN2JZFlmldRALGXjKpPQ+bAfLslCpVLhx4wYEAsFtl31oT6cTmxdeeAFarRbjx49HXV0dRo0aBYlEguXLl+OZZ57pdEDEcujHblwuqIZSraGZNzzSz8cZC8b0Qk8XWkbBmunXbou/XIzk3IpmiY1AIEBQUBAKCgoM6xmR7tOo0aKwqgEMAzTIpZRUdpCDgwP8/f0hEHQ+Eex0YsMwDP7v//4Pzz//PK5evYqamhqEhYU1W1+DWDdfuRSh3s7wkUtRXqeCj5y+RPmiv68c/X3lXIdBusHNic2txGLdukSNjY3QaFrvqiLmse9SId46kIswHxk+fiiM63AsglAohJ2dXZeTwE4nNjk5OfDz84NYLEZYWFiL+zpStZFYNoZhsHtJywHkhJDuo1+7LbmNLmGGYSASiVose0PM61RODfKqNZg6UNHlAonEOJ1u6wkKCsKNGzda3F5aWoqgoKAuBUUI6bzcsjocTCtGRR0tUmoLIvzkYBggt6wepTUNXIdDmujHHuoTT9J9Op3YsCzbanNRTU0NZac2qLhayXUIpMlfFwrw6KbTeP6n81yHQrqBTCpCLw/dEACaocgP9SoNUguqAQBRVEeq2xndFbVs2TIAuubNlStXwsHBwXCfRqPByZMnERkZabIACb+pGrUY+95B5FXU4+T/xsNLRkkt185mVwAABtEvRZsR6afA1eIaJOdWYHw/L67DsXkpeZXQaFl4ySTwofWhup3RiU1SUhIAXYtNSkpKsylZYrEYAwcOxPLly00XIeE1sZ0AzlLdZZSUU4EpA7w5jsi2sSyLszm6JvBB9EvRZjwU648J/TwxKICSWT5o1GgR5a+An4sDzYbigNGJzYEDBwAAjz32GNauXQtnZ9MV/1q9ejV++eUXXL58Gfb29hg2bBjefvtt9O3b17CNUqnEc889h23btqGhoQGTJ0/Gp59+Ci8v+pXClUg/BS4XViM5lxIbruVXKlFc3QChgEF4T5oVZSuodY5fhvV2x47e7mBZWkePC52eFbVp0ybEx8cjPj4excXFLaoFfvXVV0bv89ChQ1i4cCGGDBmCxsZG/O9//8OkSZNw6dIlODo6AgCWLl2KXbt24ccff4RcLseiRYswc+ZMHDt2rLOHQroo0k+Bbadz25yVQbrP2WzdOejn4wwHMS2MSAiXqLWGG53+5Hv11VfxyiuvYPDgwfDx8THJCdy9e3ezf2/evBmenp5ITEzEqFGjUFlZiY0bN2LLli0YN24cAF2C1a9fP5w4cQJDhw7tcgzEePpR/+ev6/qVhQJ6M3MlqWndLvoFb3uScytwKO0GogNcMCLEnetwbJZSrYGWZemHBYc6/cqvX78emzdvxrx580wZTzOVlZUAAFdXVwBAYmIi1Go1JkyYYNgmNDQU/v7+SEhIaDWxaWhoQEPDP1Mgq6qqzBavrert6QRHsRC1Kg2uFFWjn4+M65Bs1j/jayixsTW7zufjiyOZmBPrT4kNh+JTi/HstiTcGeGDj2ZHcR2OTer0dG+VSoVhw4aZMpZmtFotlixZguHDh2PAgAEAgMLCQojFYigUimbbenl5obCwsNX9rF69GnK53PDn5+dntphtlVDAIKKnAgBNN+Xau/dG4K2Z4RjWy43rUEg3+6dQXwW3gdi4pJxyaLQsZFIqiMiVTic28+fPx5YtW0wZSzMLFy7EhQsXsG3bti7tZ8WKFaisrDT85ebmmihCcrM7B/rg0WGB6ONFK0lzKcTLGbNj/OFJ0+5tTqSfAgBwubAa9SpaPoErSU2JJdWv4U6nu6KUSiU+//xz7N+/HxERES3KdX/wwQedDmrRokX4448/cPjwYfTs2dNwu7e3N1QqFSoqKpq12hQVFcHbu/XZOBKJBBKJpNOxkI6ZExvAdQiE2DQfuRSezhIUVzcgJa8SMUGuXIdkc1SNWqTk6YZQUMVh7nQ6sTl//ryhEN+FCxea3dfZgcQsy+KZZ57Bjh07cPDgwRZLM0RHR0MkEiE+Ph6zZs0CAKSlpSEnJwdxcXGdek5CrMXXx7PAMMDk/t5UKNEG6Vf63nupCMm55ZTYcCC1oAqqRi0UDiIEujm0/wBiFp1ObPT1bExp4cKF2LJlC3799Vc4Ozsbxs3I5XLY29tDLpfj8ccfx7Jly+Dq6gqZTIZnnnkGcXFxNCOKB+pVGqTkVcJHLoWfK72pu9uGQxkoqFSit6cTJTY2KtJfn9hUcB2KTUpqGrwf5aegqd4c6vQYG3NYv349KisrMWbMGPj4+Bj+tm/fbthmzZo1uPPOOzFr1iyMGjUK3t7e+OWXXziMmuj95+fzuP+zBPx2Lp/rUGxOQWU9CiqVEDDAwKaB3MT2RPnpuj/SCqs5jsQ2/TO+hrqhuGR0i83MmTM7tF1nko2OVGmUSqVYt24d1q1bZ/T+iXkN7CnH7+fyDbVUSPfRv+ah3jI4Sqh+hq2K8lfgr8UjEeLpxHUoNmlMXw8AoFmJHDP6E1AupzLtpHX/TDctb3P1d2Ie+orDgwIU3AZCOCUVCamOFIfuieqJe6J6tr8hMSujE5tNmzaZIw5iBfr7yiASMiipUeF6eT2Ns+lGZw19+9QETgixbbwaY0Ms282/FmnwYvdpaNTgQp6uojat7kwybtTguR/O4bkfznEdik1JuV6J9KJqaLW08CXXKLEhJhXVVCSMxtl0n4ziWjRqtXB1FNMUUwKtlsXPZ6/jz5QCqDXa9h9ATGL1X6mYuOYwfjhDRWC5RokNManIpmqbtNJ39wnzleH8y5Px3eOxNK6JoJeHE2RSO9SrNbhcQLOjuoNGy+JcUyt1JFUc5hwlNsSk4oLd8dKdYVg1vT/XodgUJ4kdwnxp0CgBBALGMJBfP/aKmFd6cTVqVRo4ioUI8aRlZbhGiQ0xKW+5FP8aEYSBTV1ShJDuF9001ioxmxKb7qDveh/op4BQQK2mXKPEhhALVlSlxMxPj2H1n6kdqgNFbMMgarHpVoaKw9QNxQuU2BCTu1HdgB/P5NIgum6QlFOOszkVOHTlBo2vIQYD/eQQMMD18noUVym5Dsfq6VtsqNwCP1CJUmJyF/Mr8fxP5xHg5oD7B/txHY5VO9v0gUrTvMnNnKUihHrL0KjVori6AZ60dpjZVCnVuHqjBgANHOYLSmyIyel/tWSX1qG0pgFuThKOI7JehorDtDYNucXOhcMhtqNGeXMTCwX46pEhSCuqhjt91vECXfXE5OQOIsNaNTR40XxUjVqk5FUCoL590hIlNd1DKhJibKgnnh7di+tQSBO68olZGGZl0OBFs0ktqEJDoxYKBxGC3R25DofwlFqjhYaq4RIbQokNMQt9YnOWWmzM5p/1oRQ0cJi06qlvzyD85T04d72C61CsUqNGi/f3puHvy0VopCrPvEGJDTELfWJz7nolVI30hjcHjZaFu5PEUIyNkFtptCyUai39wDCTtKJqfPz3VSzemkw/LniEBg8Tswhyd4SroxhltSpcKarGgB5yrkOyOvNHBuPxEUFQa6ibgbQuyt8F+1OLqZ6NmegTxkh/KszHJ5TYELNgGAYbHxkMf1cHmhVlRgzDQGxHH6ikdf90CVdwG4iVOtOU2ERTuQVeoa4oYjZR/i6U1JiJUq2hSsOkXRE95RAKGBRWKZFXUc91OFZHP+tzcIArx5GQm1FiQ4gF+mDfFQx5Ix7fn8zmOhTCYw5iO4T56BZHPZNVxnE01qWoSonr5fUQMFSYj28osSFm9cnf6Zjz5QnkltVxHYpVOZ1VhpKaBtiLhFyHQnhucCAtiGkO+tcz1FsGJwmN6uATOhvErPanFiM5twJnssvg5+rAdThWQanW4EJTYb4hgdQETm5vZIg7ckrrENFTwXUoVuX8dd17kMbX8A8lNsSsBge46BKbrHLcE9WT63CswrncCqg1LDydJejpYs91OITnxoV6YVyoF9dhWJ3/TumLe6N7QiSkwft8Q11RxKwMFYipGdxk9DMxhgS6Uu0MQjjCMAx6ezohwI2qfvMNJTbErPSJTVpRNaqUao6jsQ76QaDUBE6MkVdRj8uFVVyHQYjZUWJDzMpTJoWfqz1YFkjOqeA6HIun1bKG1i8aX0M66ufE6xj+1t945bdLXIdiFb4/mY1FW87iQFox16GQVlBiQ8wuuqnk/xnqjuoyZaMGs2P8MTTYFf18nLkOh1iIiJ66yt/JuRVQ05pGXRafWow/zhcgo7iG61BIK2jwMDG7IUGuiE8tpg9UE3AQ2+F/d/TjOgxiYXp5OEFuL0JlvRqX8qsw0E/BdUgWS6tlDd3B1GrKT5TYELO7N7onZg/xp7VUCOGIQMBgcIAL4i8X43RWGSU2XaAbL9gIB7EQ/X1lXIdDWkFdUcTsJHZCSmpMJCGjFDUNjVyHQSxQNBXqM4nTTa01g/xdYCekr1A+orNCulVDo4brECxWYaUSD35xAoNe3Yd6Fb2OxDj6bpPTWeW0zlgXnMqkbii+o8SGdIuDacUY9c4BLPz+LNehWKwz2boP1BAvJ9iLaSkFYpzwHnKIhQKU1DQgu5SWOOkMlmUNLTZDgqjcAl/RGBvSLRQOYuSU1aGyXg2tloWAuqaMdiaLpnmTzpOKhPjPlL7wlkvh5iTmOhyLVN3QCC+ZFBV1akT5UWLDV5TYkG7R31cGe5EQlfVqpBfXoK83TVU21mkqzEe6aP7IYK5DsGgyqQi/LRoBpVoDKS1Ay1vUFUW6hUgowKAABQDgVNMXNOm4KqUaqQW6qrHUYkMItyip4TdKbEi3GRzQNHgxkxIbYyVmlUPLAgFuDvCWS7kOh1iwM1ll+Dg+HSU1DVyHYnGUahq0bwkosSHdJiZIPyujjGZlGOlkUzIYG0StNaRrVv56Ee/vu4IT10q5DsWi3KhuQPjLezDz02NUbJTnKLEh3SbKXwE7AYOCSiWul9dzHY5FeWCIH169uz9mDurJdSjEwumT41PUcmqUM1llUGtY1DZoIKL6NbxGg4dJt3EQ2+HOCB84SuiyM1aQuyOC3B25DoNYgdggV2w+noWT1yixMcYpmuZtMegbhnSrD2dHcR0CITZN3yWcVlSN8loVXBxp6ndHULkFy0HtaYTw3F8pBdhyMgf5FdR9R7rOzUmC3p5OAGiGYkfVNDTiYn4lgH8SQ8JflNiQbqdq1OJMVhmqlGquQ7EIm45n4X87UnDoyg2uQyFWgsbZGCcxWzcrsaeLPXzk9lyHQ9pBiQ3pdvd9loB7NyTgaHoJ16HwnlKtQXJuBQCaEUVMR9/qcK7p2iK3p59BNjTYjeNISEdQYkO6XZSfAoBupWpye+dyK6Bq1MLdSUKDh4nJjOnjiZ0Lh2Pbk0O5DsUiDOwpx9QB3hgf6sl1KKQDaPAw6XZDg92w+XgW1dHoAEP9mmBXMAytr0VMQ+4gQqSDguswLMaUAT6YMsCH6zBIB1GLDel2sUGuYBggvbiGqp+2Qz8GYih1QxFCSIdQYkO6nYujGKHeMgCgVpvbUGu0SMzWTTGNCaK+fWJaOaV1eOHn83h2axLXofBacm4FMktqqVq6BaHEhnBiaLCuBYISm7alF9VA2aiBi4MIIU3TcwkxFYYBtp3OxZ8pBahTNXIdDm+t+u0ixr53EL+fL+A6FNJBlNgQTuhnF5yg6qdtCvOVIfmlSfj6XzEQCGh8DTEtP1cH9FDYo1HLGloGSXPVSjUu5Onq1wwOoIrDloISG8KJoUFu+O+UULx7bwTXofCa3F6EiJ4KrsMgViq2qeX0OM1QbNWZrHJotCwC3Bzgq6D6NZaCEhvCCbmDCAvG9EKUP/0KIoQrw3u5AwCOX6WaUq1J0NevoTFuFoUSG0J46FRmGWZ+egxfHrnGdSjEig3rrfvCTsmrRGU9VQK/laEwXy+alWhJKLEhnKlXafBrch7e3XOZ61B452j6DZzNqcD565Vch0KsmI/cHsHujtCywEkayN9M1U3ja6jisGWhAn2EMyqNFku2J4NlgYfjAuElk3IdEm8caxrzMKwXfaAS8xrW2w1gdO9H8o/TmWXQskCgmwOtD2VhKLEhnJHbi9DfV4YLeVU4ca0Ud0f24DokXqhpaDSs4TO8tzu3wRCrt2p6f4iE1Hh/q5ggV3z58GAoGzVch0KMRFcz4VRcUxPv8avUDK53OqsMjVoWPV3s4efqwHU4xMpRUtM6Z6kIE8K8cGeEL9ehECPRFU04NaypReLo1RKq7NlEP0NFP2OFkO6g1mhRXqviOgxCuoxXic3hw4cxffp0+Pr6gmEY7Ny5s9n9LMvipZdego+PD+zt7TFhwgSkp6dzEywxidggV4iEDPIq6pFVWsd1OLygrymin7FCiLltPZWDga/sxVt/0UB+ADiaXoJ391w2dAkTy8KrxKa2thYDBw7EunXrWr3/nXfewdq1a7FhwwacPHkSjo6OmDx5MpRKZTdHSkzFQWyHQU21bI6m3+A4Gu7pi4EpHESIo4HDpJt4y6SoU2lw/BrVswGA387lYd2BDPx+Lp/rUEgn8Grw8NSpUzF16tRW72NZFh9++CFefPFF3H333QCAb775Bl5eXti5cydmz57dnaESExoZ4o6TmWW4WlzDdSicEwoYfDonGlotS8sokG4zJMgVdgIGuWX1yC2rs+mxXSzL4mi6LsEbEULdwZaIVy02t5OZmYnCwkJMmDDBcJtcLkdsbCwSEhLafFxDQwOqqqqa/RF+mR3jjxMrxuOVuwdwHQpvUFJDupOTxA4D/RQAgOMZtt1qc62kFvmVSoiFAsRSxWGLZDGJTWFhIQDAy8ur2e1eXl6G+1qzevVqyOVyw5+fn59Z4yTGc3eSwFtONWwAILesjgZRE04Mb+r6PGbjMxT1rTXRAS6wFws5joZ0hsUkNp21YsUKVFZWGv5yc3O5Donchi1/qedV1GPkOwcw6t0DUDVSsTTSvfQzFI9dLYFWa7vvwyNNic3IPtQNZaksJrHx9vYGABQVFTW7vaioyHBfayQSCWQyWbM/wj+XC6vwyFen8Mim01yHwpnDV3SDpz2dpRDbWcxbk1iJQf4ucBQLUVqrwqUC2+yyV2u0hvWhRvb24Dga0lm8Gjx8O0FBQfD29kZ8fDwiIyMBAFVVVTh58iQWLFjAbXCkyxxEdjh05QbsBAxqGhrhJLGYS9NkDqXpEpvRfegDlXQ/sZ0A/xoRBHuxEO5OEq7D4URuWR1EQgYuDrqq6MQy8erbo6amBlevXjX8OzMzE8nJyXB1dYW/vz+WLFmC119/HSEhIQgKCsLKlSvh6+uLGTNmcBc0MQl/Nwf4uzogp6wOJzJKMSHMq/0HWRG1RotjTYX5RlFiQzjy3KS+XIfAqWAPJyS+OBF5FfU0gN+C8SqxOXPmDMaOHWv497JlywAAjzzyCDZv3oz//Oc/qK2txZNPPomKigqMGDECu3fvhlRKA0+twYgQd2w5mYOjV0tsLrFJyqlAdUMjXBxECO8h5zocQmyWQMDY9HR3a8CrxGbMmDG3HTzKMAxeffVVvPrqq90YFekuI3vrEpsjNlioTz++ZmSIB4T0S5FwqKxWhSPpNxDo5miYAm4LtFoWDKP7niGWjUYoEt4Y1ssdAgbIuFGL/Ip6rsPpVoeu0Pgawg9r49OxeFsytp22rRmk+1OLMOytv/HenjSuQyFdRIkN4Q25gwiRTb8QD6QVcxtMN1s2sQ8eiQugKaaEc6P76pLrw1du2FT5hQNpN1BQqUS1Us11KKSLeNUVRciUAd5wENvBy9m2xk2NDfXE2FBPrsMgBEOD3CC2EyCvoh4ZN2rQ29OZ65DMjmVZHGz6MUXvQ8tHLTaEV54c1QvfzY+1ucHDhPCFvViI2CBXAMDBNNsY75ZaUI2CSiXsRUIMDaZlFCwdJTaEcEirZbFm3xWcuFYKjQ1XeyX8oh/rpR/7Ze3+vqwr/Dq8txukIlpGwdJRYkN4qbhKibM55VyHYXYX86vwUXw65n99BlobGs9A+E2f2JzMLEO9SsNxNOb392XqhrImlNgQ3jmeUYKYN+OxeFuS1Q9ePHRF94E6rJcbREJ6OxJ+6O3pBF+5FKpGLVLyKrkOx6zKalVIyq0AAIztS4mNNaDBw4R3BvZUQCwUILfM+gcv7kvVJTZj6AOV8AjDMPhwdhQC3BzgJbPugfyNWi3mjwhCTlkdfBX2XIdDTIASG8I7jhI7DO3lhsNXbiA+tdhqE5uiKiXONf1SnNCPEhvCLzFNA4itnaezFP83LYzrMIgJUds34aVxTbU09H3f1ii+qbUm0k8BTyv/VUwIId2FEhvCS+NCddO9z2SXo7LeOgtm7btUCACYSFPbCU8dTCvGvI0n8XF8OtehmEVmSS0OX7mBhkbrHyBtSyixIbzk7+aAXh6O0GhZq1w7qlGjRXpxDQBgEiU2hKdKa1Q4kl6CXSkFXIdiFttP5+Lhr07hf79c4DoUYkKU2BDeGtc09dIau6PshAIcfn4sfl80Ar09nbgOh5BWjQv1hFDA4HJhNXLL6rgOx+QONH22jKKlTKwKJTaEt2ZF98T79w3Eqjv7cx2KWQgEDMJ7ymk1YcJbLo5iDA5wAQDsu1TEcTSmde1GDdKKqmEnYDCmDw3etyaU2BDeCvWWYVZ0T8gdRFyHYlJaLQstVRkmFkI/BszaEpu/LujGuA3r7W51nzG2jhIbQrrZqawyxLy5H6v/TOU6FELapU9sTmWVoaJOxXE0prO7KbGZOsCb40iIqVFiQ3itXqXB54cz8MhXp6xmLaV9l4pQUqNCSY31fEkQ6xXg5og+Xk7QaFkcSLOO8W65ZXVIyauEgKHB+9aICvQRXrMTMlh3IAOV9WqcyixDXC/LXnmXZVlDkz5N8yaW4o5wH3jJyuHqKOE6FJM43DTTMibIFW5O1nFM5B+U2BBeEwkFmBjmhZ8Sr2P3hQKLT2xSC6qRU1YHsZ0AI0NoJgaxDEsm9OE6BJN6KMYf0QEuUKq1XIdCzIC6ogjv6fvA/7pQaPGDbn8/nw8AGNfXE44S+l1BCBcYhkGotwyRfgquQyFmQIkN4b0RIe5wktihuLoBSbnlXIfTaSzL4vdzusRm+kBfjqMhxHhFVUrst7LZUcT6UGJDeE9iJ8T4pkUi/0op5DiazkvKrcD18no4ioWG4oOEWIqskloMXR2Pf285i2ql5S5zMv/rM3h2axKu3ajhOhRiJpTYEItwc3cUy1pmd5SHkwRPj+6FuUMDYC8Wch0OIUYJcHNAsLsjVI1ai61pU1LTgL8vF+G3c/kQCenrz1rRmSUWYXQfT7g6ijGghwzVDY1ch9Mpfq4OeGFqKFbc0Y/rUAgxGsMwhi5UfZeqpdl9oRBaFojoKYefqwPX4RAzodGLxCLYi4U4+b/x9CuLEA5NH+iLD/en40h6CcprVXBxFHMdklF2JOUBAKaF+3AcCTEn+pYgFsOSk5ofz+TiwOViqDU0vZRYrl4eTujvK0OjljUsSWApskpqkZhdDgEDzIjqwXU4xIws95uC2KyrxTW4UlTNdRgdptZo8eafqXhs82mcyizjOhxCusRSu6N+aWqtGRHiAS+ZlONoiDlRYkMsypdHrmHCB4ewZt8VrkPpsKNXS1Bep4a7kwRDgy27wCAhd0bounGScstRWW8Zs6O0Wha/nL0OAJg1iFprrB2NsSEWZXhvXbXe+NRiVNSpoHDgfx//r4Z+fW8IBQzH0RDSNT1dHPDFw4MRE+QKub1lrIqt1moxJzYAuy8WYlIYLXpp7ajFhliUfj4y9PORQaXR4vfzBVyH066KOhX+bBqLcM+gnhxHQ4hpTAzzspikBtDVwlowphd+XTicSi3YAEpsiMXRNyXrm5b5bEdSHlSNWvTzkWFgTznX4RBickq1husQCGmGEhtice6K9IVQwCAppwIZPK4eyrIstp3KBQA8GOMHhqFuKGI9EjJKcefHR/C/HSlch3JbCRml+DU5D/UqSsBsBSU2xOJ4Oksxqmll7B1n8ziOpm3ldWowDCAVCXB3JA1YJNZFIhLgQl4Vdp0vQGUdfwcRbziUgcXbkrHhUAbXoZBuQokNsUizonXjVQ5eKebtEguujmL8tXgk9i0dbVHjEQjpiCg/BUK9ndHQqMXOZH7+wMgtq8OR9BsAgHuodo3NoMSGWKRJYd74aHYkfl4wjNddPAzDUOl2YpUYhsHsIX4AgK2ncnj5A+O7E9nQssDIEHcEujtyHQ7pJpTYEIskttN170js+DnD4WpxNWosdE0rQjrqnqiekNgJcLmwGsm5FVyH00ydqhFbT+UAAB4dFshtMKRbUWJDLJ5Gy/KqUBjLsnh2azJi39iPY1dLuA6HELORO4gM6y7pB8rzxc6kfFQpGxHg5oCxfT25Dod0I0psiEU7mFaMMe8dwMu/XeQ6FIOzORW4VFAFtZZFf18Z1+EQYlYPxvoDAH47l8+bQcQsy2Lz8UwAwMNxgRBQYUybQokNsWguDmLkltXjj/P5KK5Wch0OAGD9Qd3si7sG+lpEZWRCumJwgAtmD/HDp3MGQWbPj2L2lfW6JUwcxULcN5gKY9oaflyFhHTSQD8FBvkrcDanAltO5mDJhD6cxnOlqBr7U4vAMMDTo3txGgsh3YFhGLw1K4LrMJpROIix5YmhKKlpgExKMxJtDbXYEIv36PAgAMB3J3KgatRyGsuGptaayWHe6O3pxGkshHCBT7Oj3J0kXIdAOECJDbF4Uwd4w0smQUlNA3YmcVdPI7esDr+eywcA/HsstdYQ21Lb0Ig1+67gjrVHOf2Bset8AYqr+NEtTbhBiQ2xeCKhAPNHBAMAPopPR0MjN6XTT2WWAQBG9HZHRE8FJzEQwhWhgMGWUzlILajirGBfblkdlmxPwqh3DyC3rI6TGAj3KLEhVmFeXAC8ZBLkVdTjYNoNTmKYFd0TB5ePwUvTwzh5fkK4JBUJMX+Erlt4w6EMaLXd3yX1UXw61BoWgwNcqTCmDaPEhlgFqUiIN2aE46en4zC5vzdncfi5OqCPlzNnz08Ilx6K9YdMaodrN2qx52Jhtz731eJq/HL2OgBg+eS+3frchF8osSFWY0KYFwYHunb78+ZX1OMcz6quEsIFZ6kIjzRV+X1nT1q3dgu/v/cKtCwwKcwLkX6Kbntewj+U2BCrVFylRJWye4qFvbErFTM+PYYvj1zrlucjhM+eHBUMD2cJMktq8dXRrG55zvPXK/DXhUIwDPDcJGqtsXWU2BCr801CFka9e8Aw9dqcjl8twa6UAjAAhvVyN/vzEcJ3zlIRVkwNBQB8ceQalGrzt9q8t/cKAGBGZA/09aauYFtHiQ2xOl4yKZRqLb48mon0omqzPU+jRouXf9ct5TAnNgBhtHwCIQB0Cca/x/TCrwuHQyoy70K1Wi2LSD8FJHYCLJkQYtbnIpaBEhtidSaFeWFMXw+oGrVY/uM5NGrMU1Pj2xPZuFJUAxcHEZ6bxG3FY0L4RCBg8J8pod0yM0kgYLBsYh8c+e9YBLg5mv35CP9RYkOsDsMweGtmBGRSO5y7XokNh0zfJVVcrcQH+3TN38sn96U1oQi5jfPXK0xetI9l2WY/WjydpSbdP7FclNgQq+Qtl+Llu/oD0NW2uJRfZbJ9a7UsntmShGplI/r7yjB7iL/J9k2ItXnl94u465NjeGPXJZPu99fkfNy97hguF5ruvU2sAyU2xGrdE9UDE8O8oNaweO7Hcyb7xSgQMJgXFwBPZwk+mh0FoYAxyX4JsUb6QfVfJ2Tjp8TrJtlnUZUSq367iIv5Vdh7scgk+yTWgxIbYrUYhsGb94TD01mCaeHeMGX+cWeELw7/ZywtdElIOyaGeWHxeN2g3v/tSEHK9cou7a+moRFPf5eIyno1wnvIsWAMrctGmqPEhlg1D2cJDj0/FovGhcBO2LXLPa2wGkU3La5n7tkehFiLxeNDMD7UE6pGLZ7+LhGlNQ2d2k+dqhH/2nQaSTkVkNuLsOaBgRB18X1NrA/D8mmN+W5QVVUFuVyOyspKyGSmn56r0bI4lVmG4molPJ2liAlybdZVob+/sEqJspoGuDqK4S23b7Fde/uNDnBBYnY5iquVcHeSACxQUtvQ6nMaG+vtjkGjZXEioxQJ10oAMIjr5YahwW6tPt/N+3F3lAAMUFJjXIymVlmnxs9nr+Ox4YFgmI4//8lrpfj392fBMAy++VdMt03tbu966sr2t24b6afAlpPZyC6rQ4CrA+bFBUIoYHDiWikSMkoBsIgNdINAyKC4SomyWhVcnSTwlpk2rpu3L6ysb/Y8N1/3HTm+E9dKcfxqCfIq6uGrsMfw3u6tXq/6bRMySsGyLBQOIrg5SlBep2r2HtVoWXybkNXsNRLbCVruq+k9omUBub0IlXVqFFTqYhjWyx1De7X+nunKa9XR7Tuqtf0C6PRzVdarMWPdMWSW1KKvlzM+fzjaqFlMSrUG878+g6NXS+AsscOWJ4YivKe8U8fWnTpyfm63japR2+KaEwqYdj+72/vMvd13Smuf++a4xozV0e9vi0xs1q1bh3fffReFhYUYOHAgPv74Y8TExHToseZMbHZfKMArv19CQeU/v+p95FKsmh6GKQN8Wr2/te06sl8BA7S1xtzt9tVerHcN9MFv5wpaPQYAeOGXFFTUNa/oq3AQ4a2Z4c2e73bH2tEYTU2t0eKeT4/hQl4VHh0WiBemhrbb6sKyLL49kY1Xf7+ERi2LAT1k+H7+UMjtRWaPt73rqSvbt3d+AIBhdCund2Rskqniai+2W6/72x1fa9cq0PJ6vd22N3MQC1Gv1uDmT0wBAzwxMggr7ggzal+tvWfaOnZTvrYd1dp+FQ66a/7mYzP2udKLqjHny5Morm7AhrnRmDKg4+u6vfbHJWw8mgkHsRDfPh6L6ACXDj+WKx05P7fbJimnHF8cyWx2zTPQXYu1Kk2z7Vv77L5Ze99Fbb23AJjlGusMq01stm/fjocffhgbNmxAbGwsPvzwQ/z4449IS0uDp6dnu483V2Kz+0IBFnx3Fre+mPqc9slRQfj8cGaL+2/ddv3cQS0ShNb2ezv657x1X53dJwN0aNsNTc/Xkf23F6O5fHciGy/uvAAA6KGwx/9N64epA7xbbb0prlbi/T1XsP1MLgDgroG+eHtWBOzF5u+Cau966uh10tr2nbmmOsKY69fY47jdc966n90XCvD0d2fbfeyGuYMAoEPbtuepUUGI8ncxel8bOnBeTHHOjWHMOejMcxVVKXHgcjFmx+hmE7IsizqVBo4Su9s+rqJOhae+TcTSiX0wNNitQ8/FpY6cHwBtbmOO9yfQse+i9mLg6vPbahOb2NhYDBkyBJ988gkAQKvVws/PD8888wxeeOGFdh9vjsRGo2Ux4u2/b/vr93YtLDfzkUtx9L/jDM2K7e23LQx0U571+zIm1s7ykUtx6PmxGP3ugQ7tv60Yze2vlAK8+sc/v0BiAl0xuq8HZg7qAR+5PQBg2fZk7EzOg5bVnbsVU/th/sggo7qvOqu9c3Tr62bM9gDMdv67EldXrvdbj2/4W3+jsKr9fXjLJACYDm3bHgEDuDuKUFxj3Ppk3jIJjr0wHsDtz0tXX9uO6sw56Or7+Pz1Cty3IQHj+3lieoQvnKR2KKhQ4np5HXLK6vDB/ZEQNO2XZdlueQ92VUfPD8uyKKzq3Hijzurod1F7uPj87uj39+1TZJ5RqVRITEzEihUrDLcJBAJMmDABCQkJrT6moaEBDQ3/XDhVVaaveXAqs6zdD4KOXkgFlUqcyixDXC+3Du23Lewt+zIm1s4qqFTi24SsDu+/rRjNbWq4D0b39cBnh65hw6EMnMoqw6msMozp62FIbFwcxdCyQJS/As9N7IsRId23DlR75+jW182Y7dH0/+bQlbi6cr3fenwdTVRM+YWiZWF0UqOPoSPnpauvbUd15hx09X18+MoNNDRq8WdKIf5MKWxx/5i+npgR1QMALCKpATr+HuaCKZIagLvP746wqMSmpKQEGo0GXl5ezW738vLC5cuXW33M6tWr8corr5g1ruJq016g+v2ZYr+37sPUsd4qu6zO6MeYO6bWOIjtsHRiH9w3uCd+TsxDbnkderr8U/798RFBeCQuEP5u5i8Jf6uOvh7GXifd9Tp3Nq6uxsfFdWQKxsRt7nPeldews49dOLY3xvT1xG/n8rH/UhGEAga+Cnv4KqTo5eGEcf3aH2LAN5Z6LXYGH4/VohKbzlixYgWWLVtm+HdVVRX8/PxM+hymLuWt358p9nvrPsxddjygE2vDcFkKvaeLAxa3snCer8Keg2h0Ovp6GHuddNfr3Nm4uhqfpZbUNyZuc5/zrryGnX0swzAY0EOOAT3k+N8d/Tr9/HxiqddiZ/DxWC2qAIC7uzuEQiGKippXmiwqKoK3d+uj6yUSCWQyWbM/U4sJcoWPXIrbNZIKGNz2fj0f+T/TKjuy37Ywt+zLmFg7y0cuxby4wA7vv60YbV175+jW182Y7c15/rsS183bd+V5Y4Jc4S3r2D68ZZIOb9seAQN4Ohk/U85bJunQeenqa9tRnbk+6H3cUkfPj7dMYpb34u109LuoPXw+7xaV2IjFYkRHRyM+Pt5wm1arRXx8POLi4jiLSyhgDNPibr1gmKa/J0YGtbsfBsCq6WGGgVi32297+8Et++pIrO3trz2rpodBbCfo0P5vF6Ota+96Ajp+ndy6fWevqY7qbFw3b9+Va10oYPDyXWEdeuzLd/Xv8LbteWJkEF6dEW70416+q3+756Wr59wYxl4f9D5uXUfPj349u7a2MaVbv4vaew6mjf+/+d98Pe8WldgAwLJly/DFF1/g66+/RmpqKhYsWIDa2lo89thjnMY1ZYAP1s8dBO9bfnF6y6VYP3cQVtwRhvVzB7X5i9Snabtbp861td/bXUvebeyrvX36yKV4alRQixi95VJsmDsIG+YOMtSyuJmLg6jZtNW29m9MjLauveupo9dJa9t35PwAujo2txafa4ux1297x9HW++TW676t42vrWgWaX6/tbXszR7EQt45dFTC6qd4r7ggzal+KW94z+rjNdc6N0dZ+XRxELY6N3sdt68j5ud02G+YOwlOjglpc8wx01+LN2vrsbu159d9F7X2n3Py5b+przNwsbro3AHzyySeGAn2RkZFYu3YtYmNjO/RYqjzcfqzWXHnY0lDlYao8bC2Vh20VVR42HautY9NV5k5sCCGEEGJ6Hf3+triuKEIIIYSQtlBiQwghhBCrQYkNIYQQQqwGJTaEEEIIsRqU2BBCCCHEalBiQwghhBCrQYkNIYQQQqwGJTaEEEIIsRqU2BBCCCHEathxHUB30xdarqqq4jgSQgghhHSU/nu7vQUTbC6xqa6uBgD4+flxHAkhhBBCjFVdXQ25XN7m/Ta3VpRWq0V+fj6cnZ3B3LpUbxdUVVXBz88Pubm5VrsGlbUfo7UfH2D9x0jHZ/ms/Rjp+DqPZVlUV1fD19cXAkHbI2lsrsVGIBCgZ8+eZtu/TCazyov1ZtZ+jNZ+fID1HyMdn+Wz9mOk4+uc27XU6NHgYUIIIYRYDUpsCCGEEGI1KLExEYlEglWrVkEikXAditlY+zFa+/EB1n+MdHyWz9qPkY7P/Gxu8DAhhBBCrBe12BBCCCHEalBiQwghhBCrQYkNIYQQQqwGJTaEEEIIsRqU2NzGunXrEBgYCKlUitjYWJw6deq22//4448IDQ2FVCpFeHg4/vzzz2b3syyLl156CT4+PrC3t8eECROQnp5uzkO4LWOO74svvsDIkSPh4uICFxcXTJgwocX2jz76KBiGafY3ZcoUcx/GbRlzjJs3b24Rv1QqbbaNJZ/DMWPGtDg+hmEwbdo0wzZ8OoeHDx/G9OnT4evrC4ZhsHPnznYfc/DgQQwaNAgSiQS9e/fG5s2bW2xj7PvaXIw9vl9++QUTJ06Eh4cHZDIZ4uLisGfPnmbbvPzyyy3OX2hoqBmP4vaMPcaDBw+2eo0WFhY2285Sz2Fr7y+GYdC/f3/DNnw6h6tXr8aQIUPg7OwMT09PzJgxA2lpae0+juvvQkps2rB9+3YsW7YMq1atwtmzZzFw4EBMnjwZxcXFrW5//PhxPPjgg3j88ceRlJSEGTNmYMaMGbhw4YJhm3feeQdr167Fhg0bcPLkSTg6OmLy5MlQKpXddVgGxh7fwYMH8eCDD+LAgQNISEiAn58fJk2ahLy8vGbbTZkyBQUFBYa/rVu3dsfhtMrYYwR01TJvjj87O7vZ/ZZ8Dn/55Zdmx3bhwgUIhULcd999zbbjyzmsra3FwIEDsW7dug5tn5mZiWnTpmHs2LFITk7GkiVLMH/+/GZf/p25JszF2OM7fPgwJk6ciD///BOJiYkYO3Yspk+fjqSkpGbb9e/fv9n5O3r0qDnC7xBjj1EvLS2t2TF4enoa7rPkc/jRRx81O67c3Fy4urq2eA/y5RweOnQICxcuxIkTJ7Bv3z6o1WpMmjQJtbW1bT6GF9+FLGlVTEwMu3DhQsO/NRoN6+vry65evbrV7e+//3522rRpzW6LjY1ln3rqKZZlWVar1bLe3t7su+++a7i/oqKClUgk7NatW81wBLdn7PHdqrGxkXV2dma//vprw22PPPIIe/fdd5s61E4z9hg3bdrEyuXyNvdnbedwzZo1rLOzM1tTU2O4jW/nUA8Au2PHjttu85///Ift379/s9seeOABdvLkyYZ/d/U1M5eOHF9rwsLC2FdeecXw71WrVrEDBw40XWAm1JFjPHDgAAuALS8vb3MbazqHO3bsYBmGYbOysgy38fkcFhcXswDYQ4cOtbkNH74LqcWmFSqVComJiZgwYYLhNoFAgAkTJiAhIaHVxyQkJDTbHgAmT55s2D4zMxOFhYXNtpHL5YiNjW1zn+bSmeO7VV1dHdRqNVxdXZvdfvDgQXh6eqJv375YsGABSktLTRp7R3X2GGtqahAQEAA/Pz/cfffduHjxouE+azuHGzduxOzZs+Ho6Njsdr6cQ2O19x40xWvGJ1qtFtXV1S3eg+np6fD19UVwcDDmzJmDnJwcjiLsvMjISPj4+GDixIk4duyY4XZrO4cbN27EhAkTEBAQ0Ox2vp7DyspKAGhxzd2MD9+FlNi0oqSkBBqNBl5eXs1u9/LyatHXq1dYWHjb7fX/NWaf5tKZ47vVf//7X/j6+ja7OKdMmYJvvvkG8fHxePvtt3Ho0CFMnToVGo3GpPF3RGeOsW/fvvjqq6/w66+/4rvvvoNWq8WwYcNw/fp1ANZ1Dk+dOoULFy5g/vz5zW7n0zk0VlvvwaqqKtTX15vkuueT9957DzU1Nbj//vsNt8XGxmLz5s3YvXs31q9fj8zMTIwcORLV1dUcRtpxPj4+2LBhA37++Wf8/PPP8PPzw5gxY3D27FkApvns4ov8/Hz89ddfLd6DfD2HWq0WS5YswfDhwzFgwIA2t+PDd6HNre5Nuu6tt97Ctm3bcPDgwWaDa2fPnm34//DwcERERKBXr144ePAgxo8fz0WoRomLi0NcXJzh38OGDUO/fv3w2Wef4bXXXuMwMtPbuHEjwsPDERMT0+x2Sz+HtmLLli145ZVX8OuvvzYbfzJ16lTD/0dERCA2NhYBAQH44Ycf8Pjjj3MRqlH69u2Lvn37Gv49bNgwZGRkYM2aNfj22285jMz0vv76aygUCsyYMaPZ7Xw9hwsXLsSFCxc4HbPVUdRi0wp3d3cIhUIUFRU1u72oqAje3t6tPsbb2/u22+v/a8w+zaUzx6f33nvv4a233sLevXsRERFx222Dg4Ph7u6Oq1evdjlmY3XlGPVEIhGioqIM8VvLOaytrcW2bds69CHJ5Tk0VlvvQZlMBnt7e5NcE3ywbds2zJ8/Hz/88EOLJv9bKRQK9OnTxyLOX1tiYmIM8VvLOWRZFl999RXmzZsHsVh82235cA4XLVqEP/74AwcOHEDPnj1vuy0fvgspsWmFWCxGdHQ04uPjDbdptVrEx8c3+0V/s7i4uGbbA8C+ffsM2wcFBcHb27vZNlVVVTh58mSb+zSXzhwfoBvJ/tprr2H37t0YPHhwu89z/fp1lJaWwsfHxyRxG6Ozx3gzjUaDlJQUQ/zWcA4B3VTMhoYGzJ07t93n4fIcGqu996Aprgmubd26FY899hi2bt3abJp+W2pqapCRkWER568tycnJhvit4RwCutlGV69e7dCPCy7PIcuyWLRoEXbs2IG///4bQUFB7T6GF9+FJhmCbIW2bdvGSiQSdvPmzeylS5fYJ598klUoFGxhYSHLsiw7b9489oUXXjBsf+zYMdbOzo5977332NTUVHbVqlWsSCRiU1JSDNu89dZbrEKhYH/99Vf2/Pnz7N13380GBQWx9fX1vD++t956ixWLxexPP/3EFhQUGP6qq6tZlmXZ6upqdvny5WxCQgKbmZnJ7t+/nx00aBAbEhLCKpXKbj++zhzjK6+8wu7Zs4fNyMhgExMT2dmzZ7NSqZS9ePGiYRtLPod6I0aMYB944IEWt/PtHFZXV7NJSUlsUlISC4D94IMP2KSkJDY7O5tlWZZ94YUX2Hnz5hm2v3btGuvg4MA+//zzbGpqKrtu3TpWKBSyu3fvNmzT3mvG5+P7/vvvWTs7O3bdunXN3oMVFRWGbZ577jn24MGDbGZmJnvs2DF2woQJrLu7O1tcXNztx8eyxh/jmjVr2J07d7Lp6elsSkoKu3jxYlYgELD79+83bGPJ51Bv7ty5bGxsbKv75NM5XLBgASuXy9mDBw82u+bq6uoM2/Dxu5ASm9v4+OOPWX9/f1YsFrMxMTHsiRMnDPeNHj2afeSRR5pt/8MPP7B9+vRhxWIx279/f3bXrl3N7tdqtezKlStZLy8vViKRsOPHj2fT0tK641BaZczxBQQEsABa/K1atYplWZatq6tjJ02axHp4eLAikYgNCAhgn3jiCU4+bG5mzDEuWbLEsK2Xlxd7xx13sGfPnm22P0s+hyzLspcvX2YBsHv37m2xL76dQ/3U31v/9Mf0yCOPsKNHj27xmMjISFYsFrPBwcHspk2bWuz3dq9ZdzL2+EaPHn3b7VlWN73dx8eHFYvFbI8ePdgHHniAvXr1avce2E2MPca3336b7dWrFyuVSllXV1d2zJgx7N9//91iv5Z6DllWN7XZ3t6e/fzzz1vdJ5/OYWvHBqDZ+4qP34VMU/CEEEIIIRaPxtgQQgghxGpQYkMIIYQQq0GJDSGEEEKsBiU2hBBCCLEalNgQQgghxGpQYkMIIYQQq0GJDSGEEEKsBiU2hBBCCLEalNgQQizKo48+2mJFZEII0bPjOgBCCNFjGOa2969atQofffQRqGA6IaQtlNgQQnijoKDA8P/bt2/HSy+9hLS0NMNtTk5OcHJy4iI0QoiFoK4oQghveHt7G/7kcjkYhml2m5OTU4uuqDFjxuCZZ57BkiVL4OLiAi8vL3zxxReora3FY489BmdnZ/Tu3Rt//fVXs+e6cOECpk6dCicnJ3h5eWHevHkoKSnp5iMmhJgaJTaEEIv39ddfw93dHadOncIzzzyDBQsW4L777sOwYcNw9uxZTJo0CfPmzUNdXR0AoKKiAuPGjUNUVBTOnDmD3bt3o6ioCPfffz/HR0II6SpKbAghFm/gwIF48cUXERISghUrVkAqlcLd3R1PPPEEQkJC8NJLL6G0tBTnz58HAHzyySeIiorCm2++idDQUERFReGrr77CgQMHcOXKFY6PhhDSFTTGhhBi8SIiIgz/LxQK4ebmhvDwcMNtXl5eAIDi4mIAwLlz53DgwIFWx+tkZGSgT58+Zo6YEGIulNgQQiyeSCRq9m+GYZrdpp9tpdVqAQA1NTWYPn063n777Rb78vHxMWOkhBBzo8SGEGJzBg0ahJ9//hmBgYGws6OPQUKsCY2xIYTYnIULF6KsrAwPPvggTp8+jYyMDOzZswePPfYYNBoN1+ERQrqAEhtCiM3x9fXFsWPHoNFoMGnSJISHh2PJkiVQKBQQCOhjkRBLxrBUwpMQQgghVoJ+mhBCCCHEalBiQwghhBCrQYkNIYQQQqwGJTaEEEIIsRqU2BBCCCHEalBiQwghhBCrQYkNIYQQQqwGJTaEEEIIsRqU2BBCCCHEalBiQwghhBCrQYkNIYQQQqzG/wMAtYuvMnhmqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(times, true_intensity_function(times), \"--\", label=r\"True $\\lambda$\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_ylabel(\"Intensity ($\\lambda$)\")\n",
    "ax.scatter(arrival_times, torch.zeros_like(arrival_times), label=r\"Observed Arrivals\")\n",
    "ax.legend(loc=\"best\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameterizing the intensity function using a QEP\n",
    "\n",
    "When using a QEP to parameterize an intensity function $\\lambda(t)$, we have to keep two concerns in mind:\n",
    "\n",
    "1. The intensity function is non-negative, whereas the output of the QEP can be negative\n",
    "2. The intensity function can take on large values (e.g. 900), whereas QEP inference works best when latent function values are fairly normal (e.g. $\\mathbb{E}[f] = 0$ and $\\mathbb{E}[f^2] = 1$).\n",
    "\n",
    "As a result, we will use the following transform to convert QEP latent function samples into intensity function samples:\n",
    "\n",
    "```python\n",
    "# function_samples = <samples from the GP>\n",
    "# self.mean_intensity = E[ \\lambda (t) ] - computed from data\n",
    "intensity_samples = function_samples.exp() * self.mean_intensity\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Cox process likelihood\n",
    "\n",
    "Given an intensity function $\\lambda$, the likelihood of a Cox process is\n",
    "\n",
    "$$\n",
    "\\log p(t_1, \\ldots, t_n | \\lambda) = \\left[ \\sum_{i=1}^n \\log \\lambda( t_i ) \\right] - \\int_0^T \\lambda(\\tau) d\\tau\n",
    "$$\n",
    "\n",
    "The first term, which are the arrival log intensities, is easy to compute, given the `arrival_times`:\n",
    "\n",
    "```python\n",
    "# arrival_intensity_samples = <a function involving arrival_times>\n",
    "arrival_log_intensities = arrival_intensity_samples.log().sum(dim=-1)\n",
    "```\n",
    "\n",
    "The second term, which is the estimated number of arrivals, can be computed analytically.\n",
    "However, it is much easier to compute this integral using quadrature.\n",
    "Using a grid of `quadrature_times` evenly spaced from $0$ to $T$, we can approximate this second term:\n",
    "\n",
    "```python\n",
    "# quadrature_intensity_samples = <a function involving quadrature_times>\n",
    "# max_time = T\n",
    "est_num_arrivals = quadrature_intensity_samples.mean(dim=-1).mul(self.max_time)\n",
    "```\n",
    "\n",
    "Putting this together, we have\n",
    "\n",
    "```python\n",
    "log_likelihood = arrival_log_intensities - est_num_arrivals\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the low-level Pyro/QPyTorch interface\n",
    "\n",
    "Unfortunately, this is not a likelihood function that is easily written as a QPyTorch Likelihood object.\n",
    "Because of this, we will be using the low-level interface for the Pyro/QPyTorch integration.\n",
    "This is the more logical choice anyways: we're simply trying to perform inference on the intensity function rather than constructing a predictive model.\n",
    "\n",
    "Here's how it will work. We'll use a `qpytorch.models.ApproximateQEP` object to model the non-homogeneous intensity function. This object needs to define 3 functions:\n",
    "\n",
    "- `forward(times)` - which computes the prior QEP mean and covariance at the supplied times.\n",
    "- `guide(arrival_times, quadrature_times)` - which defines the approximate QEP posterior at both arrival times and quadrature times.\n",
    "- `model(arrival_times, quadrature_times)` - which does the following 3 things\n",
    "\n",
    "  - Computes the QEP prior at arrival time and quadrature times\n",
    "  - Converts QEP function samples into intensity function samples, using the transformation defined above.\n",
    "  - Computes the likelihood of the arrivals (observations). We will use `pyro.factor` rather than a QPyTorch likelihood or a `pyro.sample` call (given that we have defined a simple function for directly computing the log likelihood).\n",
    " \n",
    "Putting it all together, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "POWER = 1.0\n",
    "class QEPModel(qpytorch.models.ApproximateQEP):\n",
    "    def __init__(self, num_arrivals, max_time, num_inducing=32, name_prefix=\"cox_qep_model\"):\n",
    "        self.power = torch.tensor(POWER)\n",
    "        self.name_prefix = name_prefix\n",
    "        self.max_time = max_time\n",
    "        self.mean_intensity = (num_arrivals / max_time)\n",
    "        \n",
    "        # Define the variational distribution and strategy of the QEP\n",
    "        # We will initialize the inducing points to lie on a grid from 0 to T\n",
    "        inducing_points = torch.linspace(0, max_time, num_inducing).unsqueeze(-1)\n",
    "        variational_distribution = qpytorch.variational.CholeskyVariationalDistribution(num_inducing_points=num_inducing, power=self.power)\n",
    "        variational_strategy = qpytorch.variational.VariationalStrategy(self, inducing_points, variational_distribution)\n",
    "\n",
    "        # Define model\n",
    "        super().__init__(variational_strategy=variational_strategy)\n",
    "\n",
    "        # Define mean and kernel\n",
    "        self.mean_module = qpytorch.means.ZeroMean()\n",
    "        self.covar_module = qpytorch.kernels.ScaleKernel(qpytorch.kernels.RBFKernel())\n",
    "    \n",
    "    def forward(self, times):\n",
    "        mean = self.mean_module(times)\n",
    "        covar = self.covar_module(times)\n",
    "        return qpytorch.distributions.MultivariateQExponential(mean, covar, power=self.power)\n",
    "\n",
    "    def guide(self, arrival_times, quadrature_times):\n",
    "        function_distribution = self.pyro_guide(torch.cat([arrival_times, quadrature_times], -1))\n",
    "\n",
    "        # Draw samples from q(f) at arrival_times\n",
    "        # Also draw samples from q(f) at evenly-spaced points (quadrature_times)\n",
    "        with pyro.plate(self.name_prefix + \".times_plate\", dim=-1):\n",
    "            pyro.sample(\n",
    "                self.name_prefix + \".function_samples\",\n",
    "                function_distribution\n",
    "            )\n",
    "\n",
    "    def model(self, arrival_times, quadrature_times):\n",
    "        pyro.module(self.name_prefix + \".gp\", self)\n",
    "        function_distribution = self.pyro_model(torch.cat([arrival_times, quadrature_times], -1))\n",
    "        \n",
    "        # Draw samples from p(f) at arrival times\n",
    "        # Also draw samples from p(f) at evenly-spaced points (quadrature_times)\n",
    "        with pyro.plate(self.name_prefix + \".times_plate\", dim=-1):\n",
    "            function_samples = pyro.sample(\n",
    "                self.name_prefix + \".function_samples\",\n",
    "                function_distribution\n",
    "            )\n",
    "        \n",
    "        ####\n",
    "        # Convert function samples into intensity samples, using the function above\n",
    "        ####\n",
    "        intensity_samples = function_samples.exp() * self.mean_intensity\n",
    "        \n",
    "        # Divide the intensity samples into arrival_intensity_samples and quadrature_intensity_samples\n",
    "        arrival_intensity_samples, quadrature_intensity_samples = intensity_samples.split([\n",
    "            arrival_times.size(-1), quadrature_times.size(-1)\n",
    "        ], dim=-1)\n",
    "\n",
    "        ####\n",
    "        # Compute the log_likelihood, using the method described above\n",
    "        ####\n",
    "        arrival_log_intensities = arrival_intensity_samples.log().sum(dim=-1)\n",
    "        est_num_arrivals = quadrature_intensity_samples.mean(dim=-1).mul(self.max_time)\n",
    "        log_likelihood = arrival_log_intensities - est_num_arrivals\n",
    "        pyro.factor(self.name_prefix + \".log_likelihood\", log_likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = QEPModel(arrival_times.numel(), max_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing inference\n",
    "\n",
    "We'll use Pyro to perform inference. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the quadrature times\n",
    "\n",
    "We'll use 128 quadrature points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "quadrature_times = torch.linspace(0, max_time, 64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Pyro SVI inference loop\n",
    "\n",
    "Now we'll use Pyro to perform inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "pyro.clear_param_store()  # so that notebook supports repeated runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "194e196c69e349f1ad879686de1506ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/200 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "num_iter = 2 if smoke_test else 200\n",
    "num_particles = 1 if smoke_test else 32\n",
    "\n",
    "\n",
    "def train(lr=0.01):\n",
    "    optimizer = pyro.optim.Adam({\"lr\": lr})\n",
    "    loss = pyro.infer.Trace_ELBO(num_particles=num_particles, vectorize_particles=True, retain_graph=True)\n",
    "    infer = pyro.infer.SVI(model.model, model.guide, optimizer, loss=loss)\n",
    "\n",
    "    model.train()\n",
    "    loader = tqdm.notebook.tqdm(range(num_iter))\n",
    "    for i in loader:\n",
    "        loss = infer.step(arrival_times, quadrature_times)\n",
    "        loader.set_postfix(loss=loss)\n",
    "\n",
    "model.train()\n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the inferred intensity function\n",
    "\n",
    "Here we'll see how much our inferred intensity function matches the true intensity function. We'll plot the intensity function that corresponds to the QEP mean, plus the intensity function at the 5th and 95th percentiles of the QEP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's a quick helper function for getting smoothed percentile values from samples\n",
    "def percentiles_from_samples(samples, percentiles=[0.05, 0.5, 0.95]):\n",
    "    num_samples = samples.size(0)\n",
    "    samples = samples.sort(dim=0)[0]\n",
    "    \n",
    "    # Get samples corresponding to percentile\n",
    "    percentile_samples = [samples[int(num_samples * percentile)] for percentile in percentiles]\n",
    "    \n",
    "    # Smooth the samples\n",
    "    kernel = torch.full((1, 1, 5), fill_value=0.2)\n",
    "    percentiles_samples = [\n",
    "        torch.nn.functional.conv1d(percentile_sample.view(1, 1, -1), kernel, padding=2).view(-1)\n",
    "        for percentile_sample in percentile_samples\n",
    "    ]\n",
    "    \n",
    "    return percentile_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAGwCAYAAACJjDBkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZBdJREFUeJzt3Xl8U3X2P/7XvTe52dokdG+hQJGyCbIpFVzQsYLK+JHfqOMGwzi4DF9wRNxwRlFxPuLouOvA6GcUZhTcxRkYUWRVrICACghlKxRaui9pumS79/dHmtDQNemStH09H48quffkfc87997k5N6b9xVUVVVBREREREERw50AERERUXfEIoqIiIgoBCyiiIiIiELAIoqIiIgoBCyiiIiIiELAIoqIiIgoBCyiiIiIiEKgCXcCPZWiKMjPz0d0dDQEQQh3OkRERNQGqqqiqqoKKSkpEMWWjzWxiOok+fn5SE1NDXcaREREFIKTJ0+iX79+LcawiOok0dHRALwrwWw2hzkbIiIiagubzYbU1FT/53hLWER1Et8pPLPZzCKKiIiom2nLpTi8sJyIiIgoBCyiiIiIiELAIoqIiIgoBLwmioiIqIt5PB64XK5wp9ErabVaSJLUIW2xiCIiIuoiqqqioKAAFRUV4U6lV7NarUhKSmr3OI4sooiIiLqIr4BKSEiA0WjkYMxdTFVV1NTUoKioCACQnJzcrvZYRBEREXUBj8fjL6BiY2PDnU6vZTAYAABFRUVISEho16k9XlhORETUBXzXQBmNxjBnQr510N7r0lhEERERdSGewgu/jloHLKKIiIiIQsAiioiIiCgELKKIiIiIQsAiioiIiCLKZZddhvnz57cYs3DhQuh0Otx6661dk1QTWET1MIqiotTuCHcaRETUw/z2t7+FIAgQBAGyLGPw4MFYvHgx3G53WPJ55JFH8Pzzz2PVqlU4cuRIWHJgEdXD7M4tx5tfH0NlLW8nQEREHeuqq67C6dOncfjwYdx///144okn8NxzzzUZ63Q6OzUXi8WC2bNnQxRF7N27t1OX1RwWUT2Iy6Ng66FiHCuuRnZBVbjTISKiVqiqihqnu8v/VFUNKV+dToekpCQMGDAAc+bMQWZmJv79738D8J6CmzdvHubPn4+4uDhMnToVAKAoCpYsWYK0tDQYDAaMHj0aH330kb/N6upq/OY3v0FUVBSSk5Px/PPPtzkft9sNo9GIffv2hdSf9uKI5T3I3rxKHC2uhgpgT245LhjYh+OREBFFsFqXByMWfdHly/158VQY5faXAAaDAaWlpf7HK1aswJw5c7Bt2zb/tCVLluCdd97BsmXLkJ6ejq1bt2LGjBmIj4/H5MmT8eCDD2LLli347LPPkJCQgD/+8Y/YvXs3xowZ0+ryH330UdjtdhZR1D7u+qNQggCkWPQ4UmRHcZUDCWZ9uFMjIqIeRlVVbNiwAV988QXuuece//T09HQ8++yz/scOhwNPP/00vvrqK0ycOBEAMGjQIHzzzTf4+9//jvHjx+Mf//gH3nnnHVxxxRUAvIVYv379Ws1h165dWLZsGaZNm8YiitpnX74NR4rsSLEaoNeIKLDV4UBBFYsoIqIIZtBK+Hnx1LAsNxRr1qxBVFQUXC4XFEXBrbfeiieeeMI/f/z48QHxR44cQU1NDa688sqA6U6nE2PHjsXRo0fhdDqRkZHhnxcTE4OhQ4e2mIeiKLj77rsxb948ZGRkYMaMGXC5XNBqtSH1K1QsonoA31EoqGd2DL1Wwq4T5bhkcBxEkaf0iIgikSAIHXJaratcfvnlWLp0KWRZRkpKCjSawNxNJlPAY7vdDgBYu3Yt+vbtGzBPp9OhrKwspDxeffVVlJSUYPHixcjNzYXL5cLBgwcxatSokNoLVfdZc9SsA6ercKiwCsnWM0ed4qJ0OFlWg5PlNRgQa2rh2URERG1jMpkwePDgNsePGDECOp0Oubm5mDx5cqP5ffr0gVarxfbt29G/f38AQHl5OQ4dOtRkPADk5eXhsccew6pVq2AymZCeng6dTod9+/axiKLgeBQVWw8VQ1UR8G3GJEs4Ve7GgdM2FlFERBQW0dHReOCBB3DfffdBURRcfPHFqKysxLZt22A2mzFr1izMnj0bDz74IGJjY5GQkIA//elPEMXmBw/4wx/+gKuvvhrTpk0DAGg0GgwfPjws10WxiOrmDhbYcLDAhmRL4LVPgiAgSqfBrhPl+MWwRMgajmZBRERd76mnnkJ8fDyWLFmCY8eOwWq1Yty4cfjjH/8IAHjuuedgt9tx7bXXIjo6Gvfffz8qKyubbGvNmjXYuHEjDhw4EDB91KhRYSmiBDXUwSKoRTabDRaLBZWVlTCbzZ2yDEVR8fetR7E3rxLpCdGN5te5PMivqMXcXwzGsKTOyYGIiNqmrq4OOTk5SEtLg17PH/2EU0vrIpjPbx6e6MayC6tw8HQVki2GJufrtRJcioL9eU1X9ERERBQ6FlHdlKKo+OZwCVweBVG65s/KWgwyfjxViRpneO5tRERE1FOxiOqmjhTbsS+/EkmWlg8Jx5pkFNsdOFRo76LMiIiIegcWUd2Qqqr45nAxnG4F0fqWBxbTSiIEFdh7qqJrkiMiIuolWER1Q8dKqrE3r/WjUD4xJhn7820or+7cO2oTERH1JhFVRHk8Hjz22GP+Oz2fc845eOqppwLuNq2qKhYtWoTk5GQYDAZkZmbi8OHDAe2UlZXhtttug9lshtVqxezZs/2jpvr89NNPuOSSS6DX65Gamhpwrx+fDz/8EMOGDYNer8eoUaPw3//+t3M6HqQimwN1LgXmVo5C+ViNMipqXDhYUNXJmREREfUeEVVE/eUvf8HSpUvx2muv4cCBA/jLX/6CZ599Fq+++qo/5tlnn8Urr7yCZcuWYfv27TCZTJg6dSrq6ur8Mbfddhv279+P9evXY82aNdi6dSvuuusu/3ybzYYpU6ZgwIAB2LVrF5577jk88cQTeOONN/wx3377LW655RbMnj0be/bswfTp0zF9+vSw3eSwPSRRgFYSsOdkOTiiBRERUceIqHGifvnLXyIxMRH/+Mc//NOuv/56GAwGvPPOO1BVFSkpKbj//vvxwAMPAAAqKyuRmJiI5cuX4+abb8aBAwcwYsQI7Ny5E+effz4AYN26dbjmmmtw6tQppKSkYOnSpfjTn/6EgoICyLIMAFi4cCFWr16NgwcPAgBuuukmVFdXY82aNf5cLrzwQowZMwbLli1rtS+dOU5U1tFS/DPrOIYkNh4bqjmVtS7YHW7cP2VIs0MiEBFR5+E4UZGjR44TNWnSJGzYsAGHDh0CAPz444/45ptvcPXVVwMAcnJyUFBQgMzMTP9zLBYLMjIykJWVBQDIysqC1Wr1F1AAkJmZCVEUsX37dn/MpZde6i+gAGDq1KnIzs5GeXm5P6bhcnwxvuWczeFwwGazBfxFErNeg6o6Fw6e5ik9IiKijhBRRdTChQtx8803Y9iwYdBqtRg7dizmz5+P2267DQBQUFAAAEhMTAx4XmJion9eQUEBEhISAuZrNBrExMQExDTVRsNlNBfjm3+2JUuWwGKx+P9SU1OD7n9nEgQBBq2EXSfK4VEi5uAjERFRh1m4cCF0Oh1uvfXWLlleRBVRH3zwAd59912sXLkSu3fvxooVK/DXv/4VK1asCHdqrXrkkUdQWVnp/zt58mS4U2ok1qRDfmUtiqrqWg8mIiKqJwhCi39PPPFEuFME4P0sfv7557Fq1SocOXKk05cXUTcgfvDBB/1HowDvDQVPnDiBJUuWYNasWUhKSgIAFBYWIjk52f+8wsJCjBkzBgCQlJSEoqKigHbdbjfKysr8z09KSkJhYWFAjO9xazG++WfT6XTQ6XShdLvL6LUiHDYFVXVuJFvCnQ0REXUXp0+f9v/7/fffx6JFi5Cdne2fFhUV1eg5Tqcz4LKZrmCxWDB79mzce++92Lt3LwYPHtypy4uoI1E1NTUQxcCUJEmCoigAgLS0NCQlJWHDhg3++TabDdu3b8fEiRMBABMnTkRFRQV27drlj9m4cSMURUFGRoY/ZuvWrXC5XP6Y9evXY+jQoejTp48/puFyfDG+5XRHkijAo6ioquMtYIiIqO2SkpL8fxaLBYIgBEyLiorCZZddhnnz5mH+/PmIi4vD1KlTAQADBw7ESy+9FNDemDFj/EevFEXBkiVL/MMbjR49Gh999FHIubrdbhiNxi75NX1EHYm69tpr8b//+7/o378/zj33XOzZswcvvPACfve73wHwHk6cP38+/vznPyM9PR1paWl47LHHkJKSgunTpwMAhg8fjquuugp33nknli1bBpfLhXnz5uHmm29GSkoKAODWW2/Fk08+idmzZ+Phhx/Gvn378PLLL+PFF1/053Lvvfdi8uTJeP755zFt2jS89957+P777wOGQehuvIddgao6V+vBRETUZaqrq5udJ0lSwC/IWooVRREGg6HFWJPJFGKWrVuxYgXmzJmDbdu2tfk5S5YswTvvvINly5YhPT0dW7duxYwZMxAfH4/JkycHncOjjz4Ku93e+4qoV199FY899hj+3//7fygqKkJKSgruvvtuLFq0yB/z0EMPobq6GnfddRcqKipw8cUXY926dQEb2Lvvvot58+bhiiuugCiKuP766/HKK6/451ssFnz55ZeYO3cuxo8fj7i4OCxatChgLKlJkyZh5cqVePTRR/HHP/4R6enpWL16NUaOHNk1L0YnUVXA7uCRKCKiSNLU6TCfa665BmvXrvU/TkhIQE1NTZOxkydPxubNm/2PBw4ciJKSkoCYzhzZKD09vcnBq5vjcDjw9NNP46uvvvKf6Rk0aBC++eYb/P3vfw+6iNq1axeWLVuGadOm9b4iKjo6Gi+99FKjw34NCYKAxYsXY/Hixc3GxMTEYOXKlS0u67zzzsPXX3/dYsyNN96IG2+8scWY7kYjCSjj7V+IiKgTjB8/Pqj4I0eOoKamBldeeWXAdKfTibFjxwbVlqIouPvuuzFv3jxkZGRgxowZcLlc0GrbdnePUERUEUWdT5ZElNpZRBERRZKzb03WkCRJAY/P/vFUQ2dfV3z8+PF25RWspk4ViqLY6OiX75pkX7/Xrl2Lvn37BsQE+2OtV199FSUlJVi8eDFyc3Phcrlw8OBBjBo1Kqh2gsEiqpeRNSIqapxQFBWiKIQ7HSIiQnDXKXVWbGeJj48P+HWfzWZDTk4OAGDEiBHQ6XTIzc0N6fonn7y8PDz22GNYtWoVTCYT0tPTodPpsG/fPhZR1HFkjQiHW4Hd6W7zDYyJiIhC9Ytf/ALLly/HtddeC6vVikWLFvmPrkVHR+OBBx7AfffdB0VRcPHFF6OyshLbtm2D2WzGrFmzAACvvfYaPv3000a/mvf5wx/+gKuvvhrTpk0D4B1ke/jw4Z1+XRSLqF5GpxFhr3PDXsciioiIOt8jjzyCnJwc/PKXv4TFYsFTTz3lPxIFAE899RTi4+OxZMkSHDt2DFarFePGjcMf//hHf0xJSQmOHj3aZPtr1qzBxo0bceDAgYDpo0aN6vQiKqJuQNyTRNoNiH3cHgXHS6txb+aQkJ5PRESh4Q2II0ePvAExdT7vgJscK4qIiKi9WET1MmcG3ORYUURERO3BIqo34oCbRERE7cYiqheSOOAmERFRu7GI6oW8A246wp0GEVGvxN9zhV9HrQMWUb2Qd8BNFzwKd2Qioq7iu/1Ic/e9o67jWwftvSUMx4nqhXQaEXVuBdUccJOIqMtIkgSr1eq/bYvRaIQg8M4RXUlVVdTU1KCoqAhWq7XRLXWCxSKqF5I1Imx1blRxwE0ioi6VlJQEoOX731Hns1qt/nXRHiyieiFZEuF0e2DnMAdERF1KEAQkJycjISHBfxNe6lparbbdR6B8WET1Qr4BN+0O7sBEROEgSVKHfZBT+PDC8l6IA24SERG1H4uo3krlrV+IiIjag0VUL8UBN4mIiNqHRVQvJUsiiygiIqJ2YBHVS+k0Iso54CYREVHIWET1UrJGhMOt8EbEREREIWIR1UvJGhFOFlFEREQhYxHVS3HATSIiovZhEdVLaSQRHpUDbhIREYWKRVQvJgCw8UgUERFRSFhE9WYqYOeAm0RERCFhEdWLaTjgJhERUchYRPViskZEKYsoIiKikLCI6sVkSUQFB9wkIiIKCYuoXowDbhIREYWORVQvxgE3iYiIQsciqhfzFlEccJOIiCgULKJ6MY3oHXCzisMcEBERBY1FVC8nAKji6TwiIqKgsYjq5VRwwE0iIqJQsIjq5bSigFI7x4oiIiIKFouoXo4DbhIREYWGRVQvJ0siKms54CYREVGwWET1crJGhMPl4VhRREREQWIR1cvJGhFODwfcJCIiChaLqF7ON2o5x4oiIiIKDouoXs434CZHLSciIgoOiyjigJtEREQhYBFFHHCTiIgoBCyiCFpRQAkH3CQiIgoKiyiCrBFRxgE3iYiIgsIiiiBrOOAmERFRsFhEEWSJA24SEREFi0UU+Qfc5FhRREREbcciivwDbrZ0JEpVVagqT/cRERH5sIgiaEQRSgsDblbVufDG1mPYkVPWxZkRERFFLhZR5NfUgJt1Lg8++P4kvj1ais2HiuBwe8KQGRERUeRhEUUAvANuVtUGXhPl9ihYvScPO3LKMTDWiNzSGuzPt4UnQSIiogjDIooAeAfcLG0wVpSqqli3rwBbDhUjxaqH1ShDgICsI6VQOBQCERERiyjyajjgpqqq2HKoGOv2FyAuSgezXgsASLLocaioCkeK7eFMlYiIKCKwiCIA3iKqosYJj6Jid245/v1DPkw6DWJMsj/GpNPA6Vbw3bFS/lKPiIh6PRZRBMA74KbTreD742X48PtTUAEkmfWN4uKjddh7qhL5lXVdnyQREVEEYRFFAM4MuLnmp3zYHW6k9jE0GWc1aFFZ68JODndARES9HIsoAgDoNBIqat0osTuRFmeCIAhNxgmCgFiTjJ3Hy1DOmxYTEVEvxiKKAACSKCAt1ohz4qMgNlNA+cRF6VBid2DPyfIuyo6IiCjysIgiv2i9FpLYcgEFAKIoIEqnwbdHSlHr5OCbRETUO7GIopAkmvXIr6zFj6cqwp0KERFRWLCIopBoJRGSKODbIyVwe5Rwp0NERNTlWERRyJItBhwrqcbBgqpwp0JERNTlWERRyAxaCYqqIusoB98kIqLeh0UUtUtitB4HTttwvLQm3KkQERF1KRZR1C7Reg2qnR58e7SER6OIiKhXYRFF7SIIApLMOuw6UY5jJdXhToeIiKjLRFwRlZeXhxkzZiA2NhYGgwGjRo3C999/75+vqioWLVqE5ORkGAwGZGZm4vDhwwFtlJWV4bbbboPZbIbVasXs2bNht9sDYn766Sdccskl0Ov1SE1NxbPPPtsolw8//BDDhg2DXq/HqFGj8N///rdzOt3NWQxa1Do92HigCIrCo1FERNQ7RFQRVV5ejosuugharRaff/45fv75Zzz//PPo06ePP+bZZ5/FK6+8gmXLlmH79u0wmUyYOnUq6urO3BD3tttuw/79+7F+/XqsWbMGW7duxV133eWfb7PZMGXKFAwYMAC7du3Cc889hyeeeAJvvPGGP+bbb7/FLbfcgtmzZ2PPnj2YPn06pk+fjn379nXNi9GNCIKAFIsBe/Mq8PNpW7jTISIi6hKCGkEXsixcuBDbtm3D119/3eR8VVWRkpKC+++/Hw888AAAoLKyEomJiVi+fDluvvlmHDhwACNGjMDOnTtx/vnnAwDWrVuHa665BqdOnUJKSgqWLl2KP/3pTygoKIAsy/5lr169GgcPHgQA3HTTTaiursaaNWv8y7/wwgsxZswYLFu2rNW+2Gw2WCwWVFZWwmw2t+t1OVvW0VL8M+s4hiRGd2i77XWkyI5hSdH4/WXnQCtFVH1ORETUJsF8fkfUJ92///1vnH/++bjxxhuRkJCAsWPH4s033/TPz8nJQUFBATIzM/3TLBYLMjIykJWVBQDIysqC1Wr1F1AAkJmZCVEUsX37dn/MpZde6i+gAGDq1KnIzs5GeXm5P6bhcnwxvuWczeFwwGazBfz1Nn37GHCwoAo/cRRzIiLqBSKqiDp27BiWLl2K9PR0fPHFF5gzZw7+8Ic/YMWKFQCAgoICAEBiYmLA8xITE/3zCgoKkJCQEDBfo9EgJiYmIKapNhouo7kY3/yzLVmyBBaLxf+XmpoadP+7O4NWgkYSsPFgEepcvKceERH1bBFVRCmKgnHjxuHpp5/G2LFjcdddd+HOO+9s0+mzcHvkkUdQWVnp/zt58mS4UwqLFIsBR4ur8f3x8nCnQkRE1KkiqohKTk7GiBEjAqYNHz4cubm5AICkpCQAQGFhYUBMYWGhf15SUhKKiooC5rvdbpSVlQXENNVGw2U0F+ObfzadTgez2Rzw1xvJGhFGrYRN2YWoqnOFOx0iIqJOE1FF1EUXXYTs7OyAaYcOHcKAAQMAAGlpaUhKSsKGDRv88202G7Zv346JEycCACZOnIiKigrs2rXLH7Nx40YoioKMjAx/zNatW+FynfmQX79+PYYOHer/JeDEiRMDluOL8S2Hmpds1eNUeS2yjpaGOxUiIqJOE1FF1H333YfvvvsOTz/9NI4cOYKVK1fijTfewNy5cwF4f0o/f/58/PnPf8a///1v7N27F7/5zW+QkpKC6dOnA/Aeubrqqqtw5513YseOHdi2bRvmzZuHm2++GSkpKQCAW2+9FbIsY/bs2di/fz/ef/99vPzyy1iwYIE/l3vvvRfr1q3D888/j4MHD+KJJ57A999/j3nz5nX569LdaEQRFoMWWw8Xo6zaGe50iIiIOkVEFVEXXHABPv30U6xatQojR47EU089hZdeegm33XabP+ahhx7CPffcg7vuugsXXHAB7HY71q1bB71e74959913MWzYMFxxxRW45pprcPHFFweMAWWxWPDll18iJycH48ePx/33349FixYFjCU1adIkfxE3evRofPTRR1i9ejVGjhzZNS9GN5cYrUeRzYGvDxeHOxUiIqJOEVHjRPUkvXGcqLMVVzngURTcmzkEKVZDuNMhIiJqVTCf35ouyol6obgoGdmFVfho1ykMTohCrElGbJSMWJMOFoMWoiiEO0UiIqKQsYiiTiMIAgbEmHC4qMp7OxgV0EoCDLKEaL0Wfa0GnJMQhUvT4yAILKiIiKh7YRFFncogSxgUFwXAe9sel0dFrcuDqjoXdh6vRXaBDaP6WhBjkltpiYiIKLJE1IXl1LMJggBZ4/3lXrLFgHPio2CrcyOnpDrcqREREQWNRRSFjSQKEAAcLbaHOxUiIqKgsYiisIrWa/Fzvg0ON++1R0RE3QuLKAorq1GLsmoHTpbVhjsVIiKioLCIorDSayU43AqO87ooIiLqZlhEUdgZtBL25VeC474SEVF3wiKKws5qlJFXXotiuyPcqRAREbUZiygKu2idBlV1bhwvqQl3KkRERG3GIorCThQFiAJwuKgq3KkQERG1GYsoighmgxbZp6tQ5+JQB0RE1D2wiKKIYDVqUVbjxIlSntIjIqLugUUURQSdRoLboyCnhKOXExFR98AiiiKGQdZgX56NQx0QEVG3wCKKIobVqMXpyloU2jjUARERRT4WURQxonUaVDvcyOHo5URE1A2wiKKIIQgCJFHAoUIOdUBERJGPRRRFFItBxqHCKlQ73OFOhYiIqEUsoiiiWI1aVNS4cLyUp/SIiCiysYiiiKKVRLgVlddFERFRxGMRRREnSidhX14lFIVDHRARUeRqVxHlcrlw8uRJZGdno6ysrKNyol7OapBRZHMgv7I23KkQERE1K+giqqqqCkuXLsXkyZNhNpsxcOBADB8+HPHx8RgwYADuvPNO7Ny5szNypV7CpJNQ7XTjeAlvAUNERJErqCLqhRdewMCBA/H2228jMzMTq1evxg8//IBDhw4hKysLjz/+ONxuN6ZMmYKrrroKhw8f7qy8qQcTBAFaScTBAlu4UyEiImqWJpjgnTt3YuvWrTj33HObnD9hwgT87ne/w7Jly/D222/j66+/Rnp6eockSr2L1aDF0SI7KmqcsBrlcKdDRETUSFBF1KpVq9oUp9Pp8Pvf/z6khIgAwGqUcbioCjuPl+HKEUnhToeIiKgR/jqPIpIkCjDrtdh2pBS2Ole40yEiImokqCNRzfF4PDh48CD27dvn//v00087omnqxRLMOhwpsmPX8XJcPiwh3OkQEREFCLqIOnbsGPbu3RtQMB0+fBgulwuyLGP48OEYNWpUZ+RKvYxGFBGt0+Lrw8W4IC0GUboOqfmJiIg6RFCfSjNmzMCqVasgCAKMRiOqq6sxbdo0LFq0CKNGjUJ6ejokSeqsXKkXSjTrcLTYjl0nyjF5SHy40yEiIvIL6pqojz76CK+88grsdjvy8/Mxb948fPnll9i5cycGDBjAAoo6nEYSYZQ1+PpQMWqcvCkxERFFjqCKqPvuuw+/+c1voNfrERUVhZdffhnbtm3Dpk2bcO6552LdunWdlSf1YkkWPU6V12L3iYpwp0JEROQXVBG1ZMkSREdHB0wbP348duzYgXvvvRc33XQTbr31VhQXF3doktS7aSUReq2IrYeLUOfyhDsdIiIiAB00xIEgCLj33nvx888/w+FwYNiwYR3RLJFfitWAk2W12J1bHu5UiIiIAHTwOFF9+/bFxx9/jH/+858d2SwRtJIIWSPi68PFPBpFREQRIagiKjc3t01x06ZNAwDk5eUFnxFRM1IsBpworcFPpyrDnQoREVFwRdQFF1yAu+++Gzt37mw2prKyEm+++SZGjhyJjz/+uN0JEvnIGhEaUcTWQ8VwupVwp0NERL1cUONE/fzzz/jf//1fXHnlldDr9Rg/fjxSUlKg1+tRXl6On3/+Gfv378e4cePw7LPP4pprrumsvKmXSrHqcazEjr15FRg/ICbc6RARUS8W1JGo2NhYvPDCCzh9+jRee+01pKeno6SkBIcPHwYA3Hbbbdi1axeysrJYQFGn0GkkSKKArYdK4PLwaBQREYVPSPfRMBgMuOGGG3DDDTd0dD5ErUqxGHCwwIbl247jqpFJSI0xhjslIiLqhXgzMup29FoJqX2M+P5EGY4UVeEXwxJwyZB4GGVuzkRE1HU6dIgDoq5i0mkwNDEaoiDgkz35eH3TEezLq4SqquFOjYiIegkWUdRtCYKABLMe58SbkFtaize/PoaVO3JRYneEOzUiIuoFWERRt6eVRAyKNyHGKGNLdjFe3XAYx4rt4U6LiIh6uJCLqFmzZmHr1q0dmQtRu5gNWgxNjMZpWx325XFATiIi6lwhF1GVlZXIzMxEeno6nn76aY5OThFBFAVE6zT4+bQNisLro4iIqPOEXEStXr0aeXl5mDNnDt5//30MHDgQV199NT766CO4XK6OzJEoKBaDFsVVDhTY6sKdChER9WDtuiYqPj4eCxYswI8//ojt27dj8ODBmDlzJlJSUnDffff5B+Ek6kpROg2qHR6cKK0JdypERNSDdciF5adPn8b69euxfv16SJKEa665Bnv37sWIESPw4osvdsQiiNpMEARIooAjRVXhToWIiHqwkIsol8uFjz/+GL/85S8xYMAAfPjhh5g/fz7y8/OxYsUKfPXVV/jggw+wePHijsyXqE3Meg0OF9pR5/KEOxUiIuqhQh7iOTk5GYqi4JZbbsGOHTswZsyYRjGXX345rFZrO9IjCo3FqEVeeS1OltUgPTE63OkQEVEPFHIR9eKLL+LGG2+EXq9vNsZqtSInJyfURRCFTKeR4PQoOF7KIoqIiDpHyKfzJk+eDJ1O12i6qqrIzc1tV1JEHUGvlXDgtI23giEiok4RchGVlpaG4uLiRtPLysqQlpbWrqSIOoLFoEVeeQ3KazjkBhERdbyQiyhVVSEIQqPpdru9xVN8RF3FrNfCVufGidLqcKdCREQ9UNDXRC1YsACA92fkjz32GIxGo3+ex+PB9u3bm7zInKirSaIAVVWRU1KNsf37hDsdIiLqYYIuovbs2QPAeyRq7969kGXZP0+WZYwePRoPPPBAx2VI1A7Rei0OnLbB7VGgkXi/bSIi6jhBF1GbNm0CANx+++14+eWXYTabOzwpoo5iMWhRYnfgdGUdUmOMrT+BiIiojUL+av7222+zgKKIZ5Ql1DoV3gKGiIg6XFBHohYsWICnnnoKJpPJf21Uc1544YV2JUbUEQRBgEYScKiwChenx4U7HSIi6kGCKqL27NkDl8vl/3dzmvrVHlG4mPVaHC22o9rhhkkX8viyREREAYL6RPFdD3X2v4kimcWgRW5ZNXLLajA8maegiYioY4R8TVRtbS1qas5cZ3LixAm89NJL+PLLLzskMaKOImtEeBQVx0s4XhQREXWckIuo6667Dv/85z8BABUVFZgwYQKef/55XHfddVi6dGmHJUjUEfRaCQcKeAsYIiLqOCEXUbt378Yll1wCAPjoo4+QlJSEEydO4J///CdeeeWVDkuQqCNYjTJOV9Sh2O4IdypERNRDhFxE1dTUIDo6GgDw5Zdf4le/+hVEUcSFF16IEydOdFiCRB0hWqdBlcONXA51QEREHSTkImrw4MFYvXo1Tp48iS+++AJTpkwBABQVFXH8KIo4oihAAHC02B7uVIiIqIcIuYhatGgRHnjgAQwcOBAZGRmYOHEiAO9RqbFjx7Y7sWeeeQaCIGD+/Pn+aXV1dZg7dy5iY2MRFRWF66+/HoWFhQHPy83NxbRp02A0GpGQkIAHH3wQbrc7IGbz5s0YN24cdDodBg8ejOXLlzda/uuvv46BAwdCr9cjIyMDO3bsaHefKLzMei0OFlTB6VbCnQoREfUAIRdRN9xwA3Jzc/H9999j3bp1/ulXXHEFXnzxxXYltXPnTvz973/HeeedFzD9vvvuw3/+8x98+OGH2LJlC/Lz8/GrX/3KP9/j8WDatGlwOp349ttvsWLFCixfvhyLFi3yx+Tk5GDatGm4/PLL8cMPP2D+/Pm444478MUXX/hj3n//fSxYsACPP/44du/ejdGjR2Pq1KkoKipqV78ovCwGLcqqncirqA13KkRE1AMIaoT9XMlut2PcuHH429/+hj//+c8YM2YMXnrpJVRWViI+Ph4rV67EDTfcAAA4ePAghg8fjqysLFx44YX4/PPP8ctf/hL5+flITEwEACxbtgwPP/wwiouLIcsyHn74Yaxduxb79u3zL/Pmm29GRUWFvxjMyMjABRdcgNdeew0AoCgKUlNTcc8992DhwoVN5u1wOOBwnLlo2WazITU1FZWVlR1+ejPraCn+mXUcQxKjO7Tdnk5VVWQXVOHCc2Ixqq8FaXEmxJhkDg5LRER+NpsNFoulTZ/f7Rq+ecOGDdiwYQOKioqgKIGnSN56662Q2pw7dy6mTZuGzMxM/PnPf/ZP37VrF1wuFzIzM/3Thg0bhv79+/uLqKysLIwaNcpfQAHA1KlTMWfOHOzfvx9jx45FVlZWQBu+GN9pQ6fTiV27duGRRx7xzxdFEZmZmcjKymo27yVLluDJJ58Mqc/UNQRBQIrVgB05ZdiRU4ZovQYDYk0YkWLGoDgT+loN0EghH5wlIqJeJuQi6sknn8TixYtx/vnnIzk5uUO+zb/33nvYvXs3du7c2WheQUEBZFmG1WoNmJ6YmIiCggJ/TMMCyjffN6+lGJvNhtraWpSXl8Pj8TQZc/DgwWZzf+SRRwLuJ+g7EkWRxWzQwmzQwqOosNW5kF1QhZ9OVcCglRAfrcOUEYm4IC023GkSEVE3EHIRtWzZMixfvhwzZ87skEROnjyJe++9F+vXr4der++QNruSTqeDTqcLdxrURpIooI9RRh+jDFVVUeP04FR5Lb4+XMIiioiI2iTkcxdOpxOTJk3qsER27dqFoqIijBs3DhqNBhqNBlu2bMErr7wCjUaDxMREOJ1OVFRUBDyvsLAQSUlJAICkpKRGv9bzPW4txmw2w2AwIC4uDpIkNRnja4N6FkEQYNJpEBelQ7HdAbvD3fqTiIio1wu5iLrjjjuwcuXKDkvkiiuuwN69e/HDDz/4/84//3zcdttt/n9rtVps2LDB/5zs7Gzk5ub6h1eYOHEi9u7dG/AruvXr18NsNmPEiBH+mIZt+GJ8bciyjPHjxwfEKIqCDRs2+GOoZzLJEqodbhRXcVRzIiJqXcin8+rq6vDGG2/gq6++wnnnnQetVhsw/4UXXgiqvejoaIwcOTJgmslkQmxsrH/67NmzsWDBAsTExMBsNuOee+7BxIkTceGFFwIApkyZghEjRmDmzJl49tlnUVBQgEcffRRz5871n2r7/e9/j9deew0PPfQQfve732Hjxo344IMPsHbtWv9yFyxYgFmzZuH888/HhAkT8NJLL6G6uhq333570K8TdR+yRoTLo6LIVoe0OFO40yEioggXchH1008/YcyYMQAQMFwAgE77yfiLL74IURRx/fXXw+FwYOrUqfjb3/7mny9JEtasWYM5c+Zg4sSJMJlMmDVrFhYvXuyPSUtLw9q1a3Hffffh5ZdfRr9+/fB///d/mDp1qj/mpptuQnFxMRYtWoSCggKMGTMG69ata3SxOfUsgiBAVVUUVdWFOxUiIuoGIm6cqJ4imHEmgsVxojrPidJqnJMQhXt+kR7uVIiIKAyC+fxu16A4X3/9NWbMmIFJkyYhLy8PAPCvf/0L33zzTXuaJQobo6xBQWUd6lyecKdCREQRLuQi6uOPP8bUqVNhMBiwe/du/2jdlZWVePrppzssQaKuZNLx4nIiImqbkIuoP//5z1i2bBnefPPNgIvKL7roIuzevbtDkiPqagathDq3giIWUURE1IqQi6js7GxceumljaZbLJZGYzkRdRe+H0UU8+JyIiJqRchFVFJSEo4cOdJo+jfffINBgwa1KymicJIlESdKa8KdBhERRbiQi6g777wT9957L7Zv3w5BEJCfn493330XDzzwAObMmdORORJ1KZMsIa+iFi6P0nowERH1WiGPE7Vw4UIoioIrrrgCNTU1uPTSS6HT6fDAAw/gnnvu6cgcibqUSadBWbUTJXYHki2GcKdDREQRKuQiShAE/OlPf8KDDz6II0eOwG63Y8SIEYiKiurI/Ii6nEErodblQZGNRRQRETUv5CIqNzcXqampkGXZf1+6hvP69+/f7uSIwkEUBagq+As9IiJqUcjXRKWlpaG4uLjR9NLSUqSlpbUrKaJw00gCTpXz4nIiImpeyEWUqqpN3iPPbrdDr9e3KymicDPJGpwsq4FH4V2RiIioaUGfzluwYAEA7zVRjz32GIxGo3+ex+PB9u3b/TcmJuqujDoJtjo3SqsdSIjmlwIiImos6CJqz549ALxHovbu3QtZlv3zZFnG6NGj8cADD3RchkRhYJI1KLTVocjGIoqIiJoWdBG1adMmAMDtt9+OV155BdHR0R2eFFG4SfUXlxfbeXE5ERE1LeRf57399tvYsGEDNmzYgKKiIihK4MCEb731VruTIwonUQDyy2vDnQYREUWokIuoxYsX48knn8T555+P5OTkJi8yJ+rOjLIGx0urm/0RBRER9W4hF1FLly7F8uXLMXPmzI7MhyhimHQaVNa6UFHjQh+T3PoTiIioVwl5iAOn04lJkyZ1ZC5EEcUoS6h2eDjoJhERNSnkIuqOO+7AypUrOzIXooiilUQoqoKiqrpwp0JERBEo5NN5dXV1eOONN/DVV1/hvPPOg1arDZj/wgsvtDs5onATIKCgkkUUERE1FnIR9dNPP/kH1dy3b1/APF6ESz2FQZZwvLSGF5cTEVEjIRdRvvGiiHoyk6xBqd2BKocbZr229ScQEVGvEfI1UUS9gVEnodrhRpGNF5cTEVGgoI9E/epXv2pT3CeffBJ0MkSRRpZEuBUVxVUODE6ICnc6REQUQYIuoiwWS2fkQRSRBEEAVKDQxovLiYgoUNBF1Ntvv90ZeRBFLL0s4XhpdbjTICKiCMNroohaYZIlFNnqUON0hzsVIiKKICyiiFph0mlQ7fSgmCOXExFRAyyiiFqh04hwuBXe/oWIiAKwiCJqhW+QTV5cTkREDbGIImoDnUZEbv3I5URERACLKKI2sRq0OFpsx4nSmnCnQkREEYJFFFEbWAxa2B1ubDlUzKNRREQEgEUUUZsIgoAUiwE/nKzA0WJ7uNMhIqIIwCKKqI3MBi3qXB5szi6GovBoFBFRb8ciiigIyRY99uZV4lBRVbhTISKiMGMRRRSEaL0WLreCTQeL4OHRKCKiXo1FFFGQUqwG/Jxvw4HTtnCnQkREYcQiiihIJp0Gigpsyi6C26OEOx0iIgoTFlFEIehrNSC7oAp78yrDnQoREYUJiyiiEBhkCQCwObsYTjePRhER9UYsoohC1NdqwJEiO348VRHuVIiIKAxYRBGFSK+VoJEEbM4uQp3LE+50iIioi7GIImqHFIsBOcXV2JNbEe5UiIioi2nCnQBRdyZrRMgaEev2ncbu3HJIogBJFKCp/78kCjBoJZybYkF6QhREUQh3ykRE1EFYRBG1U98+BhRWOnC4sAoqAFUFVFX1/9utKNhyqBgjUyy4dEg8hiRGQRBYTBERdXcsoojaSSOK6NvH0GJMVZ0Le3LLsS+/EiNTLJg8NB7pCSymiIi6MxZRRF0gWq9FtF4bUEyN6us9MsViioioe+KF5URdKFqvRXpiNGJNMnafKMfftxzFDycrwp0WERGFgEUUURj4iimHW8G3R0qh8GbGRETdDosoojBKNOtxqKgKx0rs4U6FiIiCxCKKKIyidBo4XAp25JSFOxUiIgoSiyiiMIuLkvHjyQoU2erCnQoREQWBRRRRmMWYZJTXeH+1R0RE3QeLKKIwEwQBFoMWWcfKUO1whzsdIiJqIxZRRBEgIVqHQlsd9uZVhjsVIiJqIxZRRBFAI4nQSiKyjpbC7VHCnQ4REbUBiyiiCJFk1uNYsR2HCjncARFRd8AiiihCGGQJHkXFjpxSqCoH3yQiinQsoogiSHy0Dvvybciv5HAHRESRjkUUUQSxGLw3Kd59nMMdEBFFOhZRRBFEEARYjTJ2HC9DZa0r3OkQEVELWEQRRZj4KB2Kq+rw06mKcKdCREQtYBFFFGEkUYBeK+HboyVwujncARFRpGIRRRSBksx65JbWYjdvBUNEFLFYRBFFIJ1Wgk4j4oOdJ/HF/gK4OAAnEVHEYRFFFKFSY4zQayWs3pOHd787gcoaXmhORBRJNOFOgIiaFx+tg1GW8O3RUhTbHbhhfCrS4kzhTouIiBBhR6KWLFmCCy64ANHR0UhISMD06dORnZ0dEFNXV4e5c+ciNjYWUVFRuP7661FYWBgQk5ubi2nTpsFoNCIhIQEPPvgg3G53QMzmzZsxbtw46HQ6DB48GMuXL2+Uz+uvv46BAwdCr9cjIyMDO3bs6PA+E7XGpNMgPSEKOcXVeGPrUXx7tASKwhHNiYjCLaKKqC1btmDu3Ln47rvvsH79erhcLkyZMgXV1dX+mPvuuw//+c9/8OGHH2LLli3Iz8/Hr371K/98j8eDadOmwel04ttvv8WKFSuwfPlyLFq0yB+Tk5ODadOm4fLLL8cPP/yA+fPn44477sAXX3zhj3n//fexYMECPP7449i9ezdGjx6NqVOnoqioqGteDKIGNJKIwQlRcLlVrNyei0/2nEKt0xPutIiIejVBjeCbdBUXFyMhIQFbtmzBpZdeisrKSsTHx2PlypW44YYbAAAHDx7E8OHDkZWVhQsvvBCff/45fvnLXyI/Px+JiYkAgGXLluHhhx9GcXExZFnGww8/jLVr12Lfvn3+Zd18882oqKjAunXrAAAZGRm44IIL8NprrwEAFEVBamoq7rnnHixcuLDV3G02GywWCyorK2E2mzv0dck6Wop/Zh3HkMToDm2XuoeKGicKbHXISIvBzIkDoZUi6rsQEVG3Fsznd0S/+1ZWVgIAYmJiAAC7du2Cy+VCZmamP2bYsGHo378/srKyAABZWVkYNWqUv4ACgKlTp8Jms2H//v3+mIZt+GJ8bTidTuzatSsgRhRFZGZm+mPO5nA4YLPZAv6IOoPVKCO1jxE7j5dj44FC3qyYiChMIraIUhQF8+fPx0UXXYSRI0cCAAoKCiDLMqxWa0BsYmIiCgoK/DENCyjffN+8lmJsNhtqa2tRUlICj8fTZIyvjbMtWbIEFovF/5eamhpax4nawKTToI9Rxhf7C7E/nwU7EVE4RGwRNXfuXOzbtw/vvfdeuFNpk0ceeQSVlZX+v5MnT4Y7Jerh4qN1cCsKPtl9CkVVdeFOh4io14nIImrevHlYs2YNNm3ahH79+vmnJyUlwel0oqKiIiC+sLAQSUlJ/pizf63ne9xajNlshsFgQFxcHCRJajLG18bZdDodzGZzwB9RZxsQY0JeRS0+3nUKdS5eaE5E1JUiqohSVRXz5s3Dp59+io0bNyItLS1g/vjx46HVarFhwwb/tOzsbOTm5mLixIkAgIkTJ2Lv3r0Bv6Jbv349zGYzRowY4Y9p2IYvxteGLMsYP358QIyiKNiwYYM/higSiKKAgbEm/HCyAuv2FfD6KCKiLhRRg23OnTsXK1euxGeffYbo6Gj/9UcWiwUGgwEWiwWzZ8/GggULEBMTA7PZjHvuuQcTJ07EhRdeCACYMmUKRowYgZkzZ+LZZ59FQUEBHn30UcydOxc6nQ4A8Pvf/x6vvfYaHnroIfzud7/Dxo0b8cEHH2Dt2rX+XBYsWIBZs2bh/PPPx4QJE/DSSy+huroat99+e9e/MEQt0GslJEbrselgEVJjDBg/ICbcKRER9QoRVUQtXboUAHDZZZcFTH/77bfx29/+FgDw4osvQhRFXH/99XA4HJg6dSr+9re/+WMlScKaNWswZ84cTJw4ESaTCbNmzcLixYv9MWlpaVi7di3uu+8+vPzyy+jXrx/+7//+D1OnTvXH3HTTTSguLsaiRYtQUFCAMWPGYN26dY0uNieKBH1MMqocbqzek49Esx79+hjDnRIRUY8X0eNEdWccJ4q6mqqqOFRkx/CkaNxxySCYdBH1HYmIqFvoMeNEEVHbCYKAtFgT9p+2YfWePI5oTkTUyVhEEfUgskZEP6sBWw8X461vclBo49AHRESdhcf7iXqYaL0WaXEifsqrQEm1A9eP64eRfS3hTouIqMfhkSiiHkinkTAkIRqldgfe2paD9T8XwO1Rwp0WEVGPwiKKqIcSRQFpcVHQaUR8ujsP727Pha3OFe60iIh6DJ7OI+rhEqL1MMoabDtSgmK79/TewFgjBEEId2pERN0aiyiiXiBKp8HghCgcLbLj9Y2HkWw1YFhSNPrHmJAaY4DFoGVRRUQUJBZRRL2EVhKRnhAFW50bp8prcbjQDlH0Xoje16rH0CQzBsaaMCQxigUVEVEbsIgi6kUEQYDFoIXFoAUAuBUFVXVuHCuuxs/5NkTrtfjdxWkYnswbaBMRtYYXlhP1YhpRRB+jjAGxJgxNMqPa4caWQ8VQFN7IgIioNSyiiMgvxWrAz/mVOFBgC3cqREQRj0UUEfmZdBp4FGBLdjE8PBpFRNQiFlFEFCDZoseB0zYcOM2jUURELWERRUQBTDoNFJVHo4iIWsMiiogaSbHqcaDAhv35leFOhYgoYrGIIqJGjLIGUIHN2cW85x4RUTNYRBFRk1KsBhwqrMK+fF4bRUTUFBZRRNQkgywBALZkF8HFo1FERI2wiCKiZvW1GnC4yI69ebw2iojobCyiiKhZeq0EAcDm7CI43TwaRUTUEIsoImpRitWAI0V27M2rCHcqREQRhUUUEbVIr5WgEQVsyi7m0SgiogZYRBFRq1KsBhwptGPDgUKoKgfgJCICWEQRURvoNBJio2R8vu80vjtWFu50iIgiAosoImqTuCgdZI2ET/ec4kjmRERgEUVEQUix6FHn8uD9nSeRW1oT7nSIiMKKRRQRtZkgCBgQa0JxlQMrd5xAid0R7pSIiMKGRRQRBUUUBAyKNyGnuBqrduTC7nCHOyUiorBgEUVEQdOIItLiTdh7qhIf7TrJoQ+IqFfShDsBIuqedBoJ/WOM+O5oGcx6La4emQynW0Gd24M6lwd1LgV1Lg8cbgWpMQYkWwzhTpmIqEOxiCKikJl0GiSaddhwoBDfHy+HW1HgVlR4PCrcigKPCqiKiv6xRtx5ySAkmPXhTpmIqMOwiCKidrEaZWhEEQ63B0ZZA40oQCMKkCQBkiBAUYHDRXas2pGL312chmi9NtwpExF1CF4TRUTtFqXXIDZKB4tBC5NOA51WgkYUIQgCJFHAoDgT9ufb8NGuU7x+ioh6DBZRRNTpZI2I/jFGbD9WhjU/5UNReOsYIur+WEQRUZfwXT+18WARth4uDnc6RETtxiKKiLqM1SjDpNPgPz/m44eTFeFOh4ioXVhEEVGXSjLr4VZUfPT9SeSUVIc7HSKikLGIIqIuNyDGiLJqJ1btyEWhrQ6qymukiKj74RAHRNTlBEHAoPgoHC6qwktfHYLFoEVCtB6JZj0sBi2sRi0sBi3ionQwyFK40yUiahKLKCIKC0kUMDg+CpW1LpTYncgrr4XLo0AQAECATiMi0aLHrIkDkWLlaOdE5KUoKkrsDpwsr4XFoMXghKiw5cIiiojCRiOJiI3SNZquKCocbgUnSmrw7vYTuP2iNMQ1EUdEPZ+qqii2O3CyrBa5ZTXILrChpMoBW50bEwfFsogiImpIFAUYZAnnxJtwuNCOd787gd9elAaLgaOdE/UWLo+Crw8V4/sT5SiucsDucAMAjLIEs16LOrcCBeG9npJFFBFFLI0k4pz4KOzPt2HVjlz8ZuIAGGW+bRH1dGXVTqzecwo7j5dDr5VgNWiRZNZDFAV/TGm1M4wZevHdiIgimqwRkRZnwu7ccug1Im7J6A+dhhebE/VU2QVV+Hj3SRwvqUH/GCNMusgtVSI3MyKienqthP59jPjuWCn0Wgk3jO8HjcQRWoh6ErdHwZZDxVi3rwA1Tg/SE6OgESN7P2cRRUTdgkmnQbLFgC2HimGQJVx7XkrAoX0i6r4qa1xY/UMeth8rhdmgxTnxJghC5O/fLKKIqNswG7TwqCrW/1wIWRIxLNkMnUb0/mklyJIIrSR0izdfIvI6UmTHx7tP4mhRdcSfvjtb98mUiAhAH6MMt0fFp3vyoN93GpIoQiMKkOr/tJIIs0GDCQNjcf7APtBref0UUSRyuhVsOVSE9T8XotrhQXpCVLc7Tc8iioi6nfhoHWJMMtyKAo+iwqOocNf/v9blQUWNE4cK7dh2tBiXD03EmFQrZE33enMm6skKKuuw+oc8/HCyAtZudPrubCyiiKhb8h55av4ok9OtIL+8Fsu35WBwYhQuH5qAUX0t3e6bLlFPoigqdh4vw9q9p1Foc2BArBGGbny0mEUUEfVIskbEwDgTHC4PcoqrcbQoB8OSo3H50AQMTzZD4kXpRF2qssaFNXvzkXW0FLIkYkhiFMRuePSpIRZRRNSj6bQSBsVHocbpxsHTVThUUIXUGCPO62dBemI0+scYoeXRKaJO41FUHDhtw39+ykdOcTX69jHArO8Zdx9gEUVEvYJR1mBwgreYKqisw7FiO3QaCQlmHUb3syI9MRoD44wcyJOogxRV1eHnfBt2nSjHidJqqCq6xdhPwWARRUS9ilHWwBjjfeurdXpQVu3Emp9OQysVID5ahySzHnpZgl7r/ZMl7y/+NJIIq0GLQfEmRPeQb9FEHa3a4cbBgir8kFuO7MIqVNa6YNBKSIjWd6uhC9qq5/WIiKiNDLKEvrIBAOBweVBW48T+fBs8qgpFUaFAhYAz12wIACxGLYYlmTE82Yxz4k2IMcnd8ldFRK2pc3m/ZJRVO1Fid6C4yoFSuwOS6B2PTasRoZVEyJIAjSSi1unB3rxKlFQ5AEFAXJSMRLO+21/31BIWUURE8F47lWwxtBjj9iioqHVhR04Zso6WwmLQIC0+CiNTLOjbxwAB8N9TXlVVqABU1ftLwiSzHgaZpwop8rg9CsqqnSi2O1Bid6LIVodTFbUoqXKg1ulBrcsDwPslQqsRARVQVBWK6t3OFbV+JgCLQYu0OFOv+RUsiygiojbSSCLionSIi9JBUVRU1rmwP68SP+SWQ6eRoPpLKG/xBHiLKkkQEKXXIC3OhHPio9CvjwF9+xhglJt+C1ZVFQ63glqnBy6PAqdHgdujwq0ocHlUuD0qXIoCVQVEARAFAaIgQBAAof6x1ei96z2Pknl/Vu9wK3C6FTg8Hjjdin9wVo0o1g/SGvi4O1NVFUVVDtgdbm8xr8Jb8ECFongLoGqHG8V2B/LKa5FfUQu7w41apwe+ekinlWCUJcSYZBi0Em+x1AwWUUREIRBFAX2MMvoYZaiqCqdH8c/znQL01S8eRUVVnRs/nKzAzuNlkCURZoMWA2KNGBhrgqKqsDs8qKxxwlbnQlWdG3UuxT+YqKJ6P/zOnGasV1+oCcKZI2BC/Z9R9v4qcUyqFUOTohEbpeuCV6XzuT0KSqudKLTVoaLGBZfHW1i6PApcbgW1bm+R5HB5UOPyoNbpQZ3LEzAgq0dR/UWnKAgQRSGgGE0y63FOgglJFgOSLXrERem6tLByuD0osjlQVOVARY0TsVE69LUaEBfV/KljVVVRaHPgUGEVfjxZgROl1XC4Ff/RUEBt8G/fIwGyJMIoS7AYtEi2GLp9AdnVWEQREbWTIAgt/qpPKwF6rYT4aG8h43B7YKt1Y++pSuzJrWgQ573WRFN/zYlJ1kAShfoP+jMf8qKAJj9MG55CtDvcOHDahp/yKmAxyBiWFI3z+lkwJDG621wY7/Yo3tNLVXUotDlwqqwGueU1sNW6UOPyoMGBP/8RFN/rIzW4FZDvCJOsqX8sCN6jMqr3KJVSf7TGo6pwexQcKqrCvvxKAN5i1GrUYlB8FFL7GGE1aiFrRMjSmfs1yvX3bxQFIaAtxVf01p/+Cihg1DPpO90KiqrqUGRz4HhpDfIqalDt8BaAgLfgM+k0iI/WYWhiNPrFGJHax4C4KB1K7A4cKrTjx1MVyCmpRlWdCzqN9whSQv3pY6H+P4KvvBfqi20epWw3FlFERF1Mp5EQH32mqOoognDmQ9Ji0MJi0EJRVFTUurDzeDl25pQhxiRjSFI0Es169DHKsBq13j+DHLZb4/hyLK2/JqfEXodT5bU4XVmHaocbNU4PVBXQ1BeWXXXURFFUVDvdqKpzY9uRkvpTYwIksUGRJvgKNRGicOa0GVRAAfyn04D6oklteNLXO813ulEFIEsionQaxJpkGPtIEAQBHkWF3eEdmiOnpBoAYJIlWAwybHUu2OpckCURsVE6pFh4CrcrsYgiIurBRFFAjEn23mvQo6Csxontx8r8R0YkUYBBK0IvS4iN0iHZrIfVKCNar0G0XoMonQZReg3Mei10GrFNH9Cq6r2Hoa3W7f2Qr3Wh1uU9zebyKHC4FNS4vMVRnUtBqd17/U6N0wO3RwWgeq/J0Yb3NJMoCojWawOO3PkupD5zz8Yz929UgfqjhKL3+jScKWx9L5v/156C/1psiIIAvbb511YSBX9RDHiLLrvT+9oaZQnJLJzChkUUEVEvoZFEJETrkRB9ZprLo6DO5S1mcktrcLigCgrg/6WhLAn+U1ZGWQODLMFQP4aWQRZh1Goga7zjaNnrXCiu/ym83eGGw+0tmDyqioYf8b5iw3ckR6eRYJQ1iIvSRfzo8YIgQKrP3avrf3EpigLMem2PGfW7O2MRRUTUi3mvwxIRrW88z3fBvLP+l222OhfKa5zeIy+qWn/ROwIKJK0kQq8VodNIMBk10NUXWEQ9EYsoIiJqku+Ced4Kh6hp/HpAREREFAIWUUREREQhYBFFREREFAJeE9VNqaqKU+XeAdlMOgkpVkOjmzwq9TF55bVQAaT2MaBvH2OLN4NUVBX5FbX+dpMsehRU1qHa4YFBFiEAqHEqzS6zLW36ntfc9IbPsde5UevywKCVEKXXNLvMhm0ZZREqgNog8+zNWloXHR1vlEUoKpBfcWa7TLYakF9RG7CtJlsNOF1Ri5PltRAA9O1jQL8gt9+25nX2dtZwu29LO25FwU8nK5BXUQdZEjAsORqpMaYm90nf8mqcbv89yQxa76/TfNu4oqrYe6oSFbUuWA1ajOpngUYUG7WVV16Dk+W1gOodEkCvFVFkcwAArEa5yeeF+lqF+pzWNNUmgA5fTk/XlnXTUoxbURptc6IgtPs9uqXPlOba7E7rnUVUK15//XU899xzKCgowOjRo/Hqq69iwoQJYc1pR04ZNmcXo8595jYTUToNJg+Jx+CEKADAkSI7NhwsRJ3rTMzO4+XQa0RcMTzRH9fQkSI7thwqht3h9k9reEPVs529zKY01WaUToMhiVE4VGhvNH3ykHgAaPSclpbZ1DKCzbM3a24dNfeadUR8QzuPl7dpGo6XQ68VccWwtm+/oeZ19nbfUjtfHy7G7gajjgPAwUI7tJKAKSOSAvbJll4HH60kwOUJ3Ou+PlyCsf2tuCQ93t/WhgOFAe8BTTn7eT7BvlahPqc1TbWp13pvcNvS+xsFasu6aSnmdGUt9uRWBGzzWw+XNNoWg32PbstnSkttdof1Lqiq2txnZK/3/vvv4ze/+Q2WLVuGjIwMvPTSS/jwww+RnZ2NhISEFp9rs9lgsVhQWVkJs9ncYTmt23cav39nd7Pzp41KBgCs3Xu6xXamjUpuVIi09py2ttURbbZ1mcEso7k8e7PWXr9gt5OO3K5aEq68zm6nqQKqqecAre+TbTGuvxXJFkPQbY07qwAL5rUK9TmtCWUdcB9urC3rBuiY7S8Y4/pbW9032qq59Z5TUo3RqRbMvnhQhyzHJ5jPbxZRLcjIyMAFF1yA1157DQCgKApSU1Nxzz33YOHChS0+tzOKKI+iYtKSDcg7fRrO/INNxmhFEbqkQVCivEWep6YSjlM/BwYJArSigAmDYhCfOhjWpFT858fTsFdVwnFqvy8oIB4AtH1SoI3pCwBQHDVw5B0AAMgaCReeEwMR3lvIC4KA6LhkfFOkRa3LA8XlgPN09pl2hcC2JVMMtH28O7rqccFZcKTJOACQjBZoLIkwaCVcMzIen3yVBYfbNzRgw7QFiHoTNGbv66CqKqTKPFw2NB6i7/SGv30Bst6A6NhE/+LKC06euVNnfZzvsLJG1iGqz5lv9rbi01DVpo8ISFo5MLa0EKriaTJWlDSIjjlTnFeVFkLxeL+VKWftpqKkgTkuKSDW7fSexlFRP7Bh/XMEUUKfpNQG+ebD5aiDChUbDxR5v/E36Ksc1/9M/vYiXDLQCBHe069bsotR5z6Tvy5psP/frrI8aF3VuGRIPAR47xn29eESOOrjdX2HQRC8r72r9CQ81RXwfycNuCkqoO93LgTJe6DcVXIS7qri+gDvfFkScOGgWAgCkDrifPz3QJl3ROySXLgrC8/0B754ERlpMRg4agK0eiM++yEfttPH4So7ddZaOPM66/ufB1Ef5e+bq+QEZEnEhLQYQAAUBfj2aIm3b/3OhWS0eGPLT8NVnONvRysKEEQBTrfqfx0kUx8AgLuyCM6iM7ENt18AkJMGQxMV442tKoWr6Bg0IuBWfNv7mX1EGz8QmuhYAN793t9ufczlQxMgCgK+OVIKxZzkj1Xq7HCWnPAvV6eRcEl6HETv0NswxyVjwwnviOOKsxau0lONlq3XSrhsSDyiYuIR1SfOm6/TgbLTuWd6Ux8vQIAqAFtOOOCSvaN+qh4X3JVFaJIgQNSZIBktMGglXDsqAVUlBWe9XmfIBhOMZu/rqyoKbL7YJmh1ehgt3tdXVVXvvgz1zKjiDWNlPYzWGP/jioLAbcf/Uaqq0Mi6Ru8nquI5s1nWb2dqfaw1sa8/tjTvODxul38b9t5CxvtvjVaH2H5p/tii3CPYsPdU/Xug2mC7BwRJCzlxEAxaydu3vCNQXLX+HBskDkHSQNd3uH+So+AIlDr7WbG+9wgRhoFjzsSePgSlpjLwjEWD9o2DJwTEeuy+o8xqo3jD4AkQRKk+9jA89ft9w/3O97TB4y5FQbUb56Va8MCUYehIwXx+83ReM5xOJ3bt2oVHHnnEP00URWRmZiIrK6tRvMPhgMPh8D+22WwdntOOnDIUVjngLDiM4k+fbjYuZsr/Q/TYawAAzqIcFH/6v03GrQbQ5/LZME/4/wAA7tJTKP74qWbbtVx8G6wX3eKNrSxE0YeP++d9dFasOeMG9LnstwAAj70Mhav+2Gy70eOvRUzm3d7Y2ioUvPNgs7GmUVci7pp7Uevy4MPvjuHk/93TbKxx2CWIv+5h7wNVQc7ff4+cZmINgycg4fpF/scn/vr/AR5Xk7H6Aech8eYzr//Jl28+84ZzFjl5KJJ/87z/8am//RaeqpImY7VxA5Ay+3X/47w3fw93ow95L401CX3v/j//49PL74Wz8GiTsaLJitR57/gfF7z7UOPCup4gG9D/vg/9j/PXvooVx/c0GQtBxICH/u1/WL5lOWoPZeFE09Hof/8ngEYGAFRmfYDq/ZuaiQT6/WEVJIP3A9b2/Wew/7iuUYxvm+s7521ozN5C1f7jF6j6/rMm2/wEQModS6GN9RaU1T9vRmXW+83mkDTrJX+RWHPoW1RsWQEA+LSJ2MRbn/EXUbVHd6J8wxvNtptw45MwDBrvjT3+A8rWvdJsbPz0P0IzdBIAwHFyH0r+81yzsbHTFiBq5C+8sXkHUfxJ4L7csKcN3yMcBUdQ9P6jAbHvNvh3w/cIV/HxZvfP5Qh8j3AW5eD0283vnw3fI9y2EuS/eXezsb73iFqXB6u2/oxTr89sNtb3HgF4v+ydfOnXzcY2fI9QFQ9yn7uu2dhOe49IGYrkmQ3fI2Y3/x4RPxApv3vN/zjvzT/AXZbXZKzvPcJ33V3p5y+3+T2ifMMbbX6PqNj6L9S18T2i8rsPUXuo8eenT//7PwHqi6iqXf8OeI84e7/zvUeUVjs7vIgKBouoZpSUlMDj8SAxMTFgemJiIg4ebHwUaMmSJXjyySc7NaeiqjoAgKiPhq7viGaiVEgmq/+RqDMFfMNo+O1GFAQYrHGQNSKcbgWCrIecPASNbo3uazfqzLcwQSNDTjyn/tuXWn93cNX/BNncIFaSoI3t711qwBEV779Fw5lKXxBEaKxJDfJEwDch3wcrAGhEAVJU7JkkG35zAyDqTQGvjG85QoNvgT5anQFG+cyAgpLOCNXjCvzG5luurA+IFbU6b2wTNLIcGKvRQqkvJM6oP8qmDYyVtDI8Wl2jOO88XWCsrIcoGxqEnjnKptGZ/LGqCmj0JrgM0Wj0NV4QIGoDh60WDWZIUTEQBOHMS+G/CVjgBcuS0QLJklh/zzBvjOq/QVjgsqSoGGhi+tU3J5yZXx/b8D5gkjkO2gTft2/B/xyhPlar1fg3U405AXJyeqPXy9smoNcbINZf5yGZ4wP3jbM0fC2k6Dj/Pue7J5ri7aA3VmcM6FuL7erObJeSyQI5Zaj3wVlHKQD4j4T5/i0nDW60Tfq244BYnQHa+IH1sd6b4Xq/xXv3wYBYrQ6amL4BsUJ9LABoDA32I0kLyRzv2zF9Cfhn6wwmGLTebU2UtfVH3NSAfdifTsNtu/5o05mXIPDohyCdub2JVhIgaM8eXr3BvqzV+nPwKCIETVM3ea7vm+ZMrKrgTGwTR7g0WtkfCwCSbIDqafojVNLqA2N1psZBviP8+qiAWI3RDCi+64Ia3mdPgNZkCYjVRsVAdTkDj2DW36tPio4LzMkcD42zNuC9wT/PGHi0RWNNglJXHXg2oD7ns98jtDEpUOqqGuXrFfgeoe3TFx7f9o4G9xFsMMWfQ5+UBvuR4N/vfAyyFpJWgiXMt77h6bxm5Ofno2/fvvj2228xceJE//SHHnoIW7Zswfbt2wPimzoSlZqa2qGn87KOluKWN7/rkLYA4B+zzscFaTHYmVOG2Su+75C2fDqizdY8OHUInvviUFDPOTvP3qyt6yjY7aQjt6tIzMvXzr+yjge9/YXTg1OHYFiSOajXCgh+O2mL9qwD7sNndMX7bKRobr0btVKH31aIp/M6QFxcHCRJQmFhYcD0wsJCJCUlNYrX6XTQ6Zr6xtNxJqTFILn+56FNVb4CgESzDoCAAltdi20lmXW4bGgCJFHAZUMTWmy3JQKAJIve35ZPe9ps6zLvvOQcvPNdbpuW0VyevVlr6+js16yj49sjuQPyOl3Z8j7SWjt3XnIOnv/yEJRWOpcYLUMQRBTa2v86iAIQZ9KiyN70kc+WnnfnJedAEoWgXisg+O2kLULZNrgPN9bWdaOqKgptjg7fD1siCt4Diu1dZqSvdw622QxZljF+/Hhs2LDBP01RFGzYsCHgyFRXkkQBj19bf0rhrHm+x0/8z7l44n+aO9V3xhP/c65/g2yp3Zb4Yh+/dkSjjTuUNoNZ9uPXjoCsEdu0jJby7M3asj01fM06Mr69OiKv9m7rskbEnZekNf+kek9eN7JN+2Rb3HlJGhZPHxXS82SNGPRrBQT/+rZFsNsG9+GmtXXdPPE/57YY01bBvJf79o3WniM08++GjyN5vbOIasGCBQvw5ptvYsWKFThw4ADmzJmD6upq3H777WHL6aqRyVg6YxySLIHnpZMseiydMQ5XjUzGVSOTsWzGOFiNjc8VW41aLKuPa0u7LW23DZcZTK7JFj3uvjQNyU30YdmMcVg2Y1yjeQ2f23CZzS0jmDx7s7ZsT50RH6o+QW6/reXV3HZ29nbfXDuPXDMCd1+a1ujSEQAwyZI/19aWd/bzzm5PFIC7L03DI9eMaHH/bqofvuf5BPtahfqc1jTXZh+jtlHfuA83r62fCc3FLJsxDndfmtZomxfg3Rabim/Le/Qj14xo02dKwzY7cvvqKrwmqhWvvfaaf7DNMWPG4JVXXkFGRkarz+uscaJ8PIqKHTllKKqqQ0K0HhPSYhpV6h5FxXfHSpF1tBSAiomD4nDhObEtVvRntzt+QB/sOlGOoqo6xEXpABUoqXY0u8xgcm2pD755BZW1KKt2IiZKhyRz88ts2FacSQcIQIk9uDx7s7ZsTx0VH2fSQVFVbM8pg2+7vCAtBjuPlwVsq77rl7KOlQAQMPGcWFw4KLjtt615nb2dNdzu29KO061gxbfHsfN4KYyyBteP64dJg+Oa3Cd35JShwFaHkioHymscEAQBVoOMuCgZSRYDJqTFwKOo+FfWcZwoq8GAGCNmThwIWSM2auu7o6XIOlYCRQX6GGX0MWrx46kKAMDAWFOTzwv1tQr1Oa1pqk0AHb6cnq6tnwnNxTjdSqNtThKFdr9Ht/SZ0lyb4V7vHCcqAnR2EUVEREQdL5jPb57OIyIiIgoBiygiIiKiELCIIiIiIgoBiygiIiKiELCIIiIiIgoBiygiIiKiELCIIiIiIgoBiygiIiKiELCIIiIiIgqBJtwJ9FS+geBtNluYMyEiIqK28n1ut+WGLiyiOklVVRUAIDU1NcyZEBERUbCqqqpgsVhajOG98zqJoijIz89HdHQ0hKZu8d4ONpsNqampOHnyZI+8Lx/71/319D729P4BPb+P7F/311l9VFUVVVVVSElJgSi2fNUTj0R1ElEU0a9fv05dhtls7rE7B8D+9QQ9vY89vX9Az+8j+9f9dUYfWzsC5cMLy4mIiIhCwCKKiIiIKAQsorohnU6Hxx9/HDqdLtypdAr2r/vr6X3s6f0Den4f2b/uLxL6yAvLiYiIiELAI1FEREREIWARRURERBQCFlFEREREIWARRURERBQCFlER4PXXX8fAgQOh1+uRkZGBHTt2tBj/4YcfYtiwYdDr9Rg1ahT++9//BsxXVRWLFi1CcnIyDAYDMjMzcfjw4c7sQquC6eObb76JSy65BH369EGfPn2QmZnZKP63v/0tBEEI+Lvqqqs6uxvNCqZ/y5cvb5S7Xq8PiIm0dRhM/y677LJG/RMEAdOmTfPHRNL627p1K6699lqkpKRAEASsXr261eds3rwZ48aNg06nw+DBg7F8+fJGMcHu150p2D5+8sknuPLKKxEfHw+z2YyJEyfiiy++CIh54oknGq3DYcOGdWIvmhds/zZv3tzkNlpQUBAQFynrMNj+NbV/CYKAc8891x8TSetvyZIluOCCCxAdHY2EhARMnz4d2dnZrT4vEj4LWUSF2fvvv48FCxbg8ccfx+7duzF69GhMnToVRUVFTcZ/++23uOWWWzB79mzs2bMH06dPx/Tp07Fv3z5/zLPPPotXXnkFy5Ytw/bt22EymTB16lTU1dV1VbcCBNvHzZs345ZbbsGmTZuQlZWF1NRUTJkyBXl5eQFxV111FU6fPu3/W7VqVVd0p5Fg+wd4R9htmPuJEycC5kfSOgy2f5988klA3/bt2wdJknDjjTcGxEXK+quursbo0aPx+uuvtyk+JycH06ZNw+WXX44ffvgB8+fPxx133BFQZISyTXSmYPu4detWXHnllfjvf/+LXbt24fLLL8e1116LPXv2BMSde+65Aevwm2++6Yz0WxVs/3yys7MD8k9ISPDPi6R1GGz/Xn755YB+nTx5EjExMY32wUhZf1u2bMHcuXPx3XffYf369XC5XJgyZQqqq6ubfU7EfBaqFFYTJkxQ586d63/s8XjUlJQUdcmSJU3G//rXv1anTZsWMC0jI0O9++67VVVVVUVR1KSkJPW5557zz6+oqFB1Op26atWqTuhB64Lt49ncbrcaHR2trlixwj9t1qxZ6nXXXdfRqYYk2P69/fbbqsViaba9SFuH7V1/L774ohodHa3a7Xb/tEhafw0BUD/99NMWYx566CH13HPPDZh20003qVOnTvU/bu9r1pna0semjBgxQn3yySf9jx9//HF19OjRHZdYB2lL/zZt2qQCUMvLy5uNidR1GMr6+/TTT1VBENTjx4/7p0Xq+lNVVS0qKlIBqFu2bGk2JlI+C3kkKoycTid27dqFzMxM/zRRFJGZmYmsrKwmn5OVlRUQDwBTp071x+fk5KCgoCAgxmKxICMjo9k2O1MofTxbTU0NXC4XYmJiAqZv3rwZCQkJGDp0KObMmYPS0tIOzb0tQu2f3W7HgAEDkJqaiuuuuw779+/3z4ukddgR6+8f//gHbr75ZphMpoDpkbD+QtHaPtgRr1mkURQFVVVVjfbBw4cPIyUlBYMGDcJtt92G3NzcMGUYmjFjxiA5ORlXXnkltm3b5p/e09bhP/7xD2RmZmLAgAEB0yN1/VVWVgJAo+2toUj5LGQRFUYlJSXweDxITEwMmJ6YmNjo3LxPQUFBi/G+/wfTZmcKpY9ne/jhh5GSkhKwM1x11VX45z//iQ0bNuAvf/kLtmzZgquvvhoej6dD829NKP0bOnQo3nrrLXz22Wd45513oCgKJk2ahFOnTgGIrHXY3vW3Y8cO7Nu3D3fccUfA9EhZf6Fobh+02Wyora3tkG0+0vz1r3+F3W7Hr3/9a/+0jIwMLF++HOvWrcPSpUuRk5ODSy65BFVVVWHMtG2Sk5OxbNkyfPzxx/j444+RmpqKyy67DLt37wbQMe9bkSI/Px+ff/55o30wUtefoiiYP38+LrroIowcObLZuEj5LNR0WEtEneCZZ57Be++9h82bNwdcfH3zzTf7/z1q1Cicd955OOecc7B582ZcccUV4Ui1zSZOnIiJEyf6H0+aNAnDhw/H3//+dzz11FNhzKzj/eMf/8CoUaMwYcKEgOndef31NitXrsSTTz6Jzz77LOCaoauvvtr/7/POOw8ZGRkYMGAAPvjgA8yePTscqbbZ0KFDMXToUP/jSZMm4ejRo3jxxRfxr3/9K4yZdbwVK1bAarVi+vTpAdMjdf3NnTsX+/btC9v1WcHikagwiouLgyRJKCwsDJheWFiIpKSkJp+TlJTUYrzv/8G02ZlC6aPPX//6VzzzzDP48ssvcd5557UYO2jQIMTFxeHIkSPtzjkY7emfj1arxdixY/25R9I6bE//qqur8d5777XpDTlc6y8Uze2DZrMZBoOhQ7aJSPHee+/hjjvuwAcffNDo1MnZrFYrhgwZ0i3WYVMmTJjgz72nrENVVfHWW29h5syZkGW5xdhIWH/z5s3DmjVrsGnTJvTr16/F2Ej5LGQRFUayLGP8+PHYsGGDf5qiKNiwYUPAkYqGJk6cGBAPAOvXr/fHp6WlISkpKSDGZrNh+/btzbbZmULpI+D9VcVTTz2FdevW4fzzz291OadOnUJpaSmSk5M7JO+2CrV/DXk8Huzdu9efeyStw/b078MPP4TD4cCMGTNaXU641l8oWtsHO2KbiASrVq3C7bffjlWrVgUMT9Ecu92Oo0ePdot12JQffvjBn3tPWYdbtmzBkSNH2vRFJpzrT1VVzJs3D59++ik2btyItLS0Vp8TMZ+FHXaJOoXkvffeU3U6nbp8+XL1559/Vu+66y7VarWqBQUFqqqq6syZM9WFCxf647dt26ZqNBr1r3/9q3rgwAH18ccfV7Varbp3715/zDPPPKNarVb1s88+U3/66Sf1uuuuU9PS0tTa2tou75+qBt/HZ555RpVlWf3oo4/U06dP+/+qqqpUVVXVqqoq9YEHHlCzsrLUnJwc9auvvlLHjRunpqenq3V1dRHfvyeffFL94osv1KNHj6q7du1Sb775ZlWv16v79+/3x0TSOgy2fz4XX3yxetNNNzWaHmnrr6qqSt2zZ4+6Z88eFYD6wgsvqHv27FFPnDihqqqqLly4UJ05c6Y//tixY6rRaFQffPBB9cCBA+rrr7+uSpKkrlu3zh/T2mvW1YLt47vvvqtqNBr19ddfD9gHKyoq/DH333+/unnzZjUnJ0fdtm2bmpmZqcbFxalFRUUR378XX3xRXb16tXr48GF179696r333quKoqh+9dVX/phIWofB9s9nxowZakZGRpNtRtL6mzNnjmqxWNTNmzcHbG81NTX+mEj9LGQRFQFeffVVtX///qosy+qECRPU7777zj9v8uTJ6qxZswLiP/jgA3XIkCGqLMvqueeeq65duzZgvqIo6mOPPaYmJiaqOp1OveKKK9Ts7Oyu6EqzgunjgAEDVACN/h5//HFVVVW1pqZGnTJlihofH69qtVp1wIAB6p133hm2DyhVDa5/8+fP98cmJiaq11xzjbp79+6A9iJtHQa7jR48eFAFoH755ZeN2oq09ef7ufvZf74+zZo1S508eXKj54wZM0aVZVkdNGiQ+vbbbzdqt6XXrKsF28fJkye3GK+q3mEdkpOTVVmW1b59+6o33XSTeuTIka7tWL1g+/eXv/xFPeecc1S9Xq/GxMSol112mbpx48ZG7UbKOgxlG62oqFANBoP6xhtvNNlmJK2/pvoGIGC/itTPQqG+A0REREQUBF4TRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRUTUhN/+9reYftad74mIGtKEOwEioq4mCEKL8x9//HG8/PLL4A0diKglLKKIqNc5ffq0/9/vv/8+Fi1ahOzsbP+0qKgoREVFhSM1IupGeDqPiHqdpKQk/5/FYoEgCAHToqKiGp3Ou+yyy3DPPfdg/vz56NOnDxITE/Hmm2+iuroat99+O6KjozF48GB8/vnnAcvat28frr76akRFRSExMREzZ85ESUlJF/eYiDoDiygiojZasWIF4uLisGPHDtxzzz2YM2cObrzxRkyaNAm7d+/GlClTMHPmTNTU1AAAKioq8Itf/AJjx47F999/j3Xr1qGwsBC//vWvw9wTIuoILKKIiNpo9OjRePTRR5Geno5HHnkEer0ecXFxuPPOO5Geno5FixahtLQUP/30EwDgtddew9ixY/H0009j2LBhGDt2LN566y1s2rQJhw4dCnNviKi9eE0UEVEbnXfeef5/S5KE2NhYjBo1yj8tMTERAFBUVAQA+PHHH7Fp06Ymr686evQohgwZ0skZE1FnYhFFRNRGWq024LEgCAHTfL/6UxQFAGC323HttdfiL3/5S6O2kpOTOzFTIuoKLKKIiDrJuHHj8PHHH2PgwIHQaPh2S9TT8JooIqJOMnfuXJSVleGWW27Bzp07cfToUXzxxRe4/fbb4fF4wp0eEbUTiygiok6SkpKCbdu2wePxYMqUKRg1ahTmz58Pq9UKUeTbL1F3J6gckpeIiIgoaPwqRERERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBQCFlFEREREIWARRURERBSC/x/8soOLMno0lAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the average predicted intensity function, and the intensity confidence region\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    function_dist = model(quadrature_times)\n",
    "    intensity_samples = function_dist(torch.Size([1000])).exp() * model.mean_intensity\n",
    "    lower, mean, upper = percentiles_from_samples(intensity_samples)\n",
    "\n",
    "# Plot the predicted intensity function\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "line, = ax.plot(quadrature_times, mean, label=r\"Pred $\\lambda$\")\n",
    "ax.fill_between(quadrature_times, lower, upper, color=line.get_color(), alpha=0.5)\n",
    "ax.plot(quadrature_times, true_intensity_function(quadrature_times), \"--\", color=\"k\", label=r\"True. $\\lambda$\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_ylabel(\"Intensity ($\\lambda$)\")\n",
    "ax.scatter(arrival_times, torch.zeros_like(arrival_times), label=r\"Observed Arrivals\")\n",
    "None"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}