{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cox Processes (w/ Pyro/GPyTorch Low-Level Interface)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This is a more in-depth example that shows off the power of the Pyro/GPyTorch 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 GP to model this latent intensity function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\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: 58\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": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGn0lEQVR4nO3dd3xb5b348c8jWd57xvGI7cTZO85OyCIhQNibMksLLS0NtJcy29L7K9De9EIXLfQWaBgNCQESdghZZBM7y9nOcBI7TrziPSRLz+8PycZJPCRb0tF43q+XX5aOjs756uhIX51nCikliqIoiv/SaR2AoiiKoi2VCBRFUfycSgSKoih+TiUCRVEUP6cSgaIoip8L0DqAnoiPj5cZGRlah6EoiuJV8vLyyqWUCRcv98pEkJGRQW5urtZhKIqieBUhxMmOlquiIUVRFD+nEoGiKIqfU4lAURTFz6lEoCiK4udUIlAURfFzbm01JIQoBGoBM9AipcwRQsQCS4EMoBC4VUp53p1xKYqi+DMtrghmSSlHSylzbPefBNZIKbOBNbb7iqIoipt4QtHQdcBi2+3FwPWu3NlHu4rILaykwdjiyt0ofkRKSW5hJUt3nGJZ7mmtw1F8UFWDkU/2nOFEeb1Ltu/uDmUS+EoIIYHXpJT/BJKklCW2x88CSR09UQjxIPAgQHp6eo923mQy8/j7e2mxSAIDdDx0WRYPzxxASKC+R9tT/JvZIvnP9pO8tfUkBaV1AATqddyakwbAl/tKGNwnkoz4MC3DVLyY2SLR6wTHyup4ZMkunrpyMA/N6O/0/Qh3TkwjhEiRUhYLIRKB1cAjwMdSyuh265yXUsZ0tZ2cnBzZ057F52qa2FtUzSd7zvDxnjOkxoTw59tHM65fbI+2p/iv9749xZMf5jMqNYrvTezHlAFx6ISgb3QIUkouW7SOqnoTf75jNLMHd/j7RlE61Gg088xH+QQZ9Lx44wiaTGaOltYxMCmCwICeF+QIIfLaFcu3cesVgZSy2Pa/VAjxETABOCeESJZSlgghkoFSV8aQFBnM3KHBzB2axJ0T0/m/b46TFhPqyl0qPurmcanEhwdx+dBLv+SFECz54SQeejuPBxbn8l/zBvGTWQM0iFLxNo1GM9//9w62n6jgp7MGIKUk2KBneEqUy/bptjoCIUSYECKi9TYwD9gHfAzca1vtXmClu2KalBXH6/eNJzEymBazhf1nqt21a8VLldc1c/fr2zleVkeAXtdhEmiVGhPKBz+ewjUj+7Jo1WHe2dbhMC+K0qZ9EvjfW0fx83mDEEK4fL/uvCJIAj6yvagA4D9Syi+FEDuAZUKIB4CTwK1ujKnNn74u4PVNJ/jkkWkMSAzXIgTFwzW3mHno7Tz2FVdT22RfY4Ngg56XbxtNbZOJo7Z6BEXpzC8/2NuWBG4Yk+q2/bq1jsBZelNH0JlzNU3M/9M3JEeF8NFPphAUoCqQlQv91/t7WJ5XxCt3juXqkckOPdfYYulV2a7iH7YcLedYeT13T+rnku13VkegzkybpMhgFt08igMlNfzPl4e1DkfxMKsPnGN5XhGPzB7gcBIA2pLArlPnefqjfLzxB5jiOsYWCwBTBsS7LAl0RSWCdi4fmsQ9k/vx+qYTbD5arnU4igd579tTDEqK4JHZ2b3azqGztfxn+ynezy1yUmSKtzNbJLe8uoVXNxzTLAaVCC7y9FVDmDM4kQCd6ytoFO/x6t3jePP+8b0u3rktJ40JmbE8//lBymqbnRSd4s2W7jjNnqJqUqJDNItBJYKLBBv0vH7feCZmxWkdiuIBis43UNtkwqDX0dcJH1SdTvDCDSNoNJr53WcHnBCh4s1qmkz871eHGZ8Rw4IeFDk6i0oEnahtMvHS6iNUN5q0DkXR0FMf5nPdK5uxWJxXpj8gMZwfXpbJyt1nOFhS47TtKt7nlXVHqag38qsFQ93STLQzXjlnsTucqmzgr2sLaDS28MzVQ7UOR9HAtuMVbCwo55mrhqBzclHhg9P7kxQZTFaCGn7CX9U0mXhn60luGpvKyNRoTWNRiaATw/pGceOYVN7aepKHZvQnPjxI65AUN5JS8tJXR0iMCOLuyc5vxREVauCeyRlO367iPSKDDXz6s+mEBWnfVF0VDXXh4Vn9MZotvLn5hNahKG62saCcbwsreWT2AIINrvugrtxdzC+X73HZ9hXP1FrUmBkfRmJEsMbRqETQpf4J4Vw5vA9vbT1JbZOqK/An6w+XkRIdwq3j01y6n3M1TSzLLWLXKTUXkz958YuD/GBxLmYn1j31hkoE3Xh45gAmZcVRY+eQAopv+PU1Q/l84XSX9zD/3sR+RAQH8MbmQpfuR/Ec5+uNvLv9FOFBevQe0kxdJYJuDE+J4v/uydG0ja/iXnXN1qQfFWJw+b7CggK4fXwan+eXUFLd6PL9KdpbvLWQBqOZH8/0nNFoVSKw0/GyOvacrtI6DMXFzlY3Mf53X7NiV7Hb9nnP5AyklLy1VY1O6uuMLRbe2XaSWYMSGNQnQutw2qhEYAcpJQ++ncdvPt6vdSiKi729rZDmFjPj+nU5N5JTpcWG8sjsbMZnuG+fija+OnCW8jqjx7UYU4nADkII7pyQzu7TVewrVnMW+Komk5n/bD/F3KFJpMW6d7Kix+YOVLOY+YHJWXH89tphXDYwQetQLqASgZ1uGptKsEHHu9vV5buvWrX/LOcbTNyr0a+1czVNLPn2lCb7VtwjLjyIe6dkeEwlcSuVCOwUFWrgmpF9Wbn7DDWqKalPWpZ7mvTYUCZpNM7UF/klPPVhvpopz0e9tbWQT/ac0TqMDqlE4IC7JvWjxSJVpbGPevGGkfzhppFOH07CXteNTiFQr1NDVPugRqOZP646zFcHzmkdSodUInDAqLRodjxzOdOzPat8T3GO9LhQJvfXbtTZmLBA5g1LYsXuYppbzJrFoTjfVwfOUtPUwh0u7qDYUyoROKi1bbmn9AhUes9skfxy+R7yTmrfu/fWnDSqGkys9tBfjkrPrNhVTN+oYM2KHbujEoGDTGYLN/59My+tVtNZ+orNR8tZllvkER26pg6IJy02hCPn1ET3vqKstplvCsq5bkyKZsWO3VGJwEEGvY7wYAMrdp1x6hj1inY+2lVMVIiBuUO1b76p1wlWPzaDn88dqHUoipOcq2liYFIEN45J0TqUTqlE0AM3jkmhuKqRHYWVWoei9FKTycxX+89y5fA+Lh9XyF6to52azBaNI1GcYXhKFF8snE52kuf0JL6YSgQ9MG9YEqGBej5y4zAEimusO1RKvdHMNaP6ah3KBZ75KJ87/rlN6zCUXqpqMFLf7PkDVqpE0AOhgQHMH96Hz/JLVOsOL2eWkpx+MUzMjNU6lAukxISQe/I8pysbtA5F6YV/bTzBpBfWeHwyUImgh+6e1I+nrxqCVNUEXm3ByL4s//EUAvSe9VG4ZqT1CuWz/BKNI1F6SkrJZ/kljEyLIizIsyeD9Kyz34uMSY/hjgnpLp29SnGtczVNGFs8sxw+LTaU0WnRHtsTVene4XO1nCiv58rhyVqH0i2VCHqhoq6Zd7ef9NgvE6VrT36wl+te2ax1GJ26ZlRf9p+p4ViZakrqjT7PP4sQcMWwPlqH0i3Pvl7xcDtPVfHMR/tIjQllhoeNJqh0rabJxKaj5dw3JUPrUDp1zUjrL8loN0yQozjfF/klTMiIJSEiSOtQuqUSQS9Mz44nPCiAL/JLVCLwMusOlWIyS+YP99xfa4mRwTwwLVPrMJQe+vPtY2jyksYkqmioF4INeuYMSWTV/rOqzbeXWbX/LPHhQYxJ8+zJYOqaW/hwZxHnapq0DkVx0NC+kYxN9+zzq5VKBL101YhkzjeY2H5cdS7zFk0mM+sPlzF3aJLHdvlvda6miZ8v28MXqvWQV1m06hB5J73nO8HtiUAIoRdC7BJCfGq7nymE2C6EOCqEWCqECHR3TL0xY2ACIQY9+WrmMq8RqNfx9gMTeGBahtahdKt/QjgDEsNZtV8NQuctTlc28Mq6Y+w6VaV1KHbT4opgIXCw3f0/AC9LKQcA54EHNIipx4INerY9PYcfz+yvdSiKnXQ6wbh+sQxI9Nwu/+1dMSyJbwsrOV9v1DoUxQ5fH7Qm7cuHaD92lb3cmgiEEKnA1cC/bPcFMBtYbltlMXC9O2NyhijVqsNrmC2SFz8/yKGzNVqHYrd5Q/tgtkjWHCrVOhTFDl8fPMeAxHAy4sO0DsVu7r4i+BPwS6C1ZjUOqJJStva/LgI6HKJPCPGgECJXCJFbVlbm8kAd0WK28MC/d/DahmNah6J0I+/keV775jhHS72nbf7I1CiSo4LVzHheoKbJWl/oTVcD4MZEIIRYAJRKKfN68nwp5T+llDlSypyEBM9qqhmg11FRb+RzVaHn8dYcPEeATnhVc18hBF8uvIz/d/1wrUNRunGyvIGEiCDmDk3UOhSHuPOKYCpwrRCiEHgPa5HQn4FoIURrf4ZUwCuH9Jw7NIk9RdWUqmZ+Hm3toVImZMYSEexdxXlRod4Vr78akRrFlidne02z0VZuSwRSyqeklKlSygzgdmCtlPJ7wDrgZttq9wIr3RWTM80ZYv0FoMpxPdfpygYKSuuYPdi7fq21+sWyPfxxlZoZz1NZLBKLRSKEwFr96T08oR/BE8DPhRBHsdYZvK5xPD0yKCmC1JgQ1hxUzfw8VWFFPTGhBmZ5aSI432Dk4z1nkGrIW4+0/UQlE15Ywz4vbEquyRATUsr1wHrb7ePABC3icCYhBPdNyaBFTV/psaZnJ5D77Fw8vA9Zp2YNTmTtoVKOldUzIDFc63CUi6w/XEp1o9GrWgu1UmMNOdEPpmdpHYLSCeslu3VOYG81e3Aiv8I6TpJKBJ5n/eEyxmfEEu7hcw90xBOKhnxKg7GFo6W1WoehXGTd4VJmLFrv1UM6p0SHMLhPBGtVPZTHKalu5PC5WmYO8p7WaO2pROBkP3l3Jw+93aMWsooLrT1USnldM6kxIVqH0it3TEhnbL9oVU/gYTYctvZtmjHQO+ufVCJwsmnZCRwrq1dzzXoQKSUbjpQxdUA8QQHePaPcvVMyePyKwV7XKsXXDe0byY9m9GdgkncW2alE4GStHZXWH/Gs3s/+7ER5PUXnG7nMizqRdcXYYuFkRb3WYSjtjEyN5skrvTdBq0TgZP0TwkiNCWm7VFS0940tKc/I9o1E8MiSndz35g6tw1BsTlU0sPt0FRYvbjGoEoGTCSGYOSiBLcfKafaS2Yl83ZDkSB68LIv0uFCtQ3GKiZlxtqscVfzoCf7z7Slu/scW6o0t3a/sobyvnZMX+P7UTG4fn45Bp/KsJ5iYFcfErDitw3CaywbGA7CpoJzbJ6RrHI2ysaCMsf1ivG7YkvbUN5ULZCWEMzwlyuNnv/IHRecbOHy21qda2fRPCKdPZDAbC8q1DsXvVdYb2X+mhukD4rUOpVdUInCRbccr+Md6NSy11t7Zdoqr/7KReqPvFNMJIZiWHc/mY+WYvbhc2hdsOWZNxlOzVSJQOrDlaDmLVh2ipsmkdSh+7Zsj1st2b+zt2ZUHpmXyz7tzUNec2tpyrIKIoABGpkRpHUqvqETgIlMHxGORsO1Yhdah+K3yumYOlNR41dwD9hqSHMmEzFhV/KixXy8YytKHJhOg9+6vUu+O3oONSY8hxKBn81FVjquVLbYkPNXLy287k3eykre2Fmodhl8LNugZ2jdS6zB6TSUCFwkM0DEhM5bN6opAM1uPlRMRHMBwH/igduSr/ef4f58eoMGLmy16s68PnGPRqkM0mby//kklAheaNiCe+uYWalU9gSaevHIIb31/gtdftndmWnY8JrNk+/FKrUPxSyt2F7M8r4igAO8/v7z/FXiw+6ZmsOXJ2V7dvtibRYUYGONlUwY6YnxGLIEBOtWMVAMWi2TLsQqmDoj32mEl2lOJwIUMep1PnCTeaMORMl5Zd9Sne3cHG/RMyIhl01E1nIm7HTxbQ2W9kan9faP+SSUCF3t720mu/dsmn+rQ5A0+2lnEm5sLCfTRYqFW07LjKattVs2U3ay1EYivNETw7U+JBzDoBHuLqjla6r0TongbKSWbj1UwpX+cz1+R3Tclg9xn5xKpih/dqsFoZlRaNH2igrUOxSlUInCx1l8Mqhmp+xwrq6Ostpkp/X1nfKHOBBv0Xj39prd69PKBrHh4itZhOI1KBC6WFhtKemwom46qZqTusvmob/cfuNj7uae56R9bVPGjm7QON+1LV5sqEbjB1AHxbD9eQYvZonUofqGirpn+CWGkxfrGsNPdsUhJ3snzXj0fszf509dHuOrPG33q8+xbA7B4qCuH90Gvg3qjmagQlXtd7efzBvHo5QO1DsNtJtmG2N56rIIBiREaR+P7Nh0tJ8ig86n+Kb7zSjzYZQMT+N31I4gKURV67uJPY/Ckx4bSNyqYbapjmcvVNJnYU1TNNB8rdlSJwE3MFqku3d3gXxuPc/0rm32i27+9hBBMyopj2/EKVU/gYtuPV2K2SJ+rf3I4EQghwoQQelcE48teXn2E+X/6hkYfGhffE31TUE59cwvBBv86RecNS2LGoASfmnfBE205Vk6wQceY9GitQ3GqbusIhBA64Hbge8B4oBkIEkKUA58Br0kpj7o0Sh8wrl8MJrNk16nzTPGxXxOewthiYceJSm7NSdU6FLebPzyZ+cOTtQ7D503MjCMuLJCgAN/6oWHPFcE6oD/wFNBHSpkmpUwEpgHbgD8IIe5yYYw+IScjBp2wzlymuMbu01U0msx+m2illJTWNGkdhk+bP7wPP52drXUYTmdPq6HLpZSX9F+XUlYCHwAfCCFULWg3IoINDE+JYtsJVaHnKpuPlqMT37Wi8TfPfbyflXvOsPPZuX5VWe4uxVWNWCyS1JgQn+pDAHZcEXSUBHqyjgITM2PZfarKryoy3Wlwnwjun5rpt62zRqRGU9Vg4vC5Wq1D8UlvbDrBnJc2YPSh/gOtHOpHIIRIA4YBw4ERwDApZY6dzw0GvgGCbPtdLqX8jRAiE3gPiAPygLullEZH4vIWt+akMS07QQ0J4CJXjkjmyhH+W04+uf93/QmGJPvmZDxa2n6igrHp0T5XPwB2XBEIIR4SQmwRQlQBR4AfAOHAx8CdDuyrGZgtpRwFjAbmCyEmAX8AXpZSDgDOAw849Aq8SHZSBDMGJmDwoY4onqKstpmy2matw9BUSnQI6bGhqh7KBaobTRw4U8PETN8sdrTnG+kp4DFgHPApEAy8IaX8QEp5xN4dSavWhvQG258EZgPLbcsXA9fbu01vtLeoimU7Tmsdhs95e2shk15cQ12zf0/bOCkrlu0nKtvGw1GcI7ewEov03fonexLBAinldinlMSnlLcArwCdCiMdsTUvtJoTQCyF2A6XAauAYUCWlbP30FgEpnTz3QSFErhAit6zMeyfiWLn7DM+u3KfqCZxs24lKhvWNJDzIv0dN+d7EfvzhphFYVMcyp9p+opJAve/1H2hlT2XxvovufwFMAGKBzY7sTEppllKOBlJt2xjswHP/KaXMkVLmJCQkOLJbjzIpKw5ji4U9p6u0DsVnNJnM7D5VxcTMWK1D0dyotGjmD0/2qXFwPMEPp2fxxn3jfbajoj11BJfUbEopm6WUvwLu7Wydrkgpq7D2T5gMRAshWn/GpQLFjmzL20zIiEUI1LgwTrTrVBVGs8VnL9sddbCkhi/yS7QOw6ckRAQxLdt3+6fYcx29TgjxAbBSSnmqdaEQIhBIFUI8g/VL/d9dbUQIkQCYpJRVQogQYC7WiuJ1wM1YWw7dC6zsyQvxFlGhBob0iWT7iQrA9zqmaGHb8QqEgJwMdUUAsHhLIZ/llzBvWJ+2Fmomk4mioiKamlSHM0cZWywYWyyEBunReUn/geDgYFJTUzEY7GtKbU8imA98H1hia+pZhbXCWA98BfxJSrnLju0kA4tt4xTpgGVSyk+FEAeA94QQvwN2Aa/bFbkXm5QVx4rdxZgtUjUldYJbx6cxJDnCb/sPXGxy/zje23GaA2dqGJEaBUBRURERERFkZGT4XGcoVyupbqS8zsiQ5Eiv6KgnpaSiooKioiIyMzPtek63iUBK2QT8Hfi7rQdxPNBoK95xJLi9wJgOlh/HWl/gNxZens1TVw1WScBJUqJDSIkO0ToMj9E2P8Hx8rZE0NTUpJJAD9U3mwk16L0iCYB1NNq4uDgcaVTjUI2SlNIkpSxxNAkoF4oKMai+BE5ScK6WJd+e8vtmo+0lRQaTFR/G1mMX9idQScBxZouk0WgmLMi7Kokdfa/Vt5FGXt90gmdX5Gsdhtf7PP8sT3+Uj9msmku2NzErjj1F1ao/QS81GFuQSMJ8vFmySgQaOV3ZwPK8IowtvjduiTttP1HBkD6RRIWq+oH2fj53IJuemOU1xRmeqrnFgk4IQgNVIgBACPGIECLGlcH4k0lZcTSZLOwtqtI6FK/V3GJm56nzTMxSrYUulhAR5FFfXhUVFYwePZrRo0fTp08fUlJS2u4bjc4dWiw/P59+/frxj3/8o9fbig8PYmhypM/X5zlyRZAE7BBCLBNCzHe074ByoQm2zk9qXJie21tUTZNJ9R/ozL83n2DRqkNahwFAXFwcu3fvZvfu3fzoRz/isccea7sfGBjYtp6UEould1fJI0aM4L333uOtt97qbdiAf8x/bXcikFI+i7Xh++vAfUCBEOIFIUR/F8Xm02LDAhncJ4Ltan6CHjt0thYhrJ30lEsdLKnlnW2nPL6eoLCwkEGDBnHPPfcwfPhwNm7cyPDhw9se/+Mf/8hzzz0HwDvvvMOECRMYPXo0Dz30EGZzx0O1JCYmsn///l7FVdfUwvGyOppbfH84GIeuHaWUUghxFjgLtAAxwHIhxGop5S9dEaAvmzesD8XnG7UOw2vdPakf14xMJjo0sPuV/dDErFiW5p7m4NmaS37x3fba1kvWXzAymbsnZ9BoNHPfm99e8vjN41K5JSeNynojP34n74LHlj40uVexFhQUsHjxYiZNmkRhYWGH6xw8eJClS5eyefNmDAYDDz/8MO+++y733HPPJes++eSTNDc3c/LkSfr169ejmOqaW6hvNhOg8/2qVLsTgRBiIXAPUA78C3hcSmmyDTxXAKhE4KCfzx2odQheTyWBzk20FZltP17JZA8vPevXrx+TJk3qcp01a9aQl5fH+PHjAWhsbCQxMfGS9b744gvq6+u5+uqr2b9/f48TQX1zCyGBep+vHwDHrghigRullCfbL5RSWoQQC5wbln9pMpl9djArV9lbVMWfvy7g2QVDyYwP0zocj5QSHUJabAjbT1QwOe7CY9TVL/iQQH2Xj8eGBfb6CuBiYWHfxRcQEHBBPUHrsBhSSu69915efPHFTrfT1NTEE088wccff8ybb77Jvn37uOqqqxyOx2KRNJjMxIf7xw8NR655gi9OAkKIPwBIKQ86NSo/cv+b3/Lg23ndr6hcYNPRctYcKiUy2HNaxniiOYOTCPGyHxlJSUmUlpZSUVFBc3Mzn376KQBz5sxh+fLllJaWAlBZWcnJkxd8JfG73/2Oe+65h4yMDEaMGMG+ffsu2b49GowtSCkJ86CWV67kSCKY28GyK50ViL9Kiw0lt7ASkw/Og+pK245XMjApnLjwIK1D8WjPXTuMP91+ycguHs1gMPDrX/+aCRMmMHfuXAYPto5WP3ToUH73u98xb948Ro4cydy5cykp+W6U1cOHD7N69WoeffRRgF4lAoQgPCjA63oU95SQ3UxgIYT4MfAwkIV1IplWEcBmKeVdrguvYzk5OTI3N9fdu3WJz/aW8JP/7OSjh6cwJl1107BHi9nCqN9+xU3jUvnv64Z3/wSFAwcOMnToEK3DUNzo4MGDDBly4XsuhMjraJ55e64I/gNcg3WO4mva/Y3TIgn4mtbOUGp+AvvtO1NDvdHss/PHOtv3/72D8w3O7bTlyyxSYu5lXwZvY88MZdVSykIp5R1SypPt/tQ3lxPEhweRnRhum59AsUezyczY9Oi2TnlK16JDDTSbzHR39a9YNRjNHDhT41cDGXZbEyKE2CSlnCaEqMU62TxAa3sqKaWMdFl0fuInswao0UgdMDErjg8fnqp1GF5jUmYcZlMZzS0W1TrNDvXNLUgg2OA/n0l75iOYZvsf4fpw/NP1Y1K0DsFrmC0Sk1l9oTliYlYsRw6XUd/coo6bHeqbWwgx6P2iI1krRwadu0UIEWG7/awQ4kMhhHc1R/BgBedq2VdcrXUYHm//mWpGPvcVGwvsn3TD36XHhhKgE9T7UVFHT1mkpMFo9vlhpy/mSMr7lZSyVggxDbgc65hDr7omLP/zo3fy+N+vDmsdhsfbfrwSo9nCoCR1gWovIQThwQFEqqk8u9VoNGORvj//wMUcSQStIy9dDfxTSvkZ4B/d7txgUlYcOwrP06L6E3Rp2/EKsuLDSIwM1joUrxIeFKD5cBxFRUVcd911ZGdn079/fxYuXNg2BPW///1vfvrTn2oaH0BggI6+0SGEBVqL0MLDwztdd8WKFQghOHTIsRFep0yZ0qPYCgsLLxiMz5kcSQTFQojXgNuBz4UQQQ4+X+nCpKw46ppb2H+mRutQPJbZIvm2sFLNP9BDxhYzzSZtRtKUUnLjjTdy/fXXU1BQwJEjR6irq+OZZ55x2T5bWhwvCjPodcSHBxFgR+ONJUuWMG3aNJYsWWLX/lvvb9myxeG4XM2RL/JbgVXAPNucxTHA464Iyh99159ANSPtzMGSGmqbWtT8Az10rKyeczXNdq27YlcxU3+/lswnP2Pq79eyYldxr/a9du1agoODuf/++wHQ6/W8/PLLvPHGGzQ0NABw+vRpZs6cSXZ2Nr/97W8B2gaPGzVqFMOHD2fp0qUA5OXlMWPGDMaNG8cVV1zR1sN45syZPProo+Tk5PD888/Tr1+/tnGL6uvrSUtLw2QycezYMebPn8+4ceOYPn06hw4dwiIlew4cZtKkyYwYMYJnn32209dTV1fHpk2beP3113nvvffalq9fv57p06dz7bXXMnTo0Evuw3dXGbfffjufffZZ23Pvu+8+li9fTmFhIdOnT2fs2LGMHTu2w8Sxf//+tuG4R44cSUFBQc/eGBtHCsLMQDBwixCi/fO+6lUECgCJEcFkJYSx/UQlD81QUzx0JD48iF/OH8Tk/ioR9ERYYAD1tjF0uppXasWuYp76MJ9G29VDcVUjT31onV+7py3c9u/fz7hx4y5YFhkZSXp6OkePHgXg22+/Zd++fYSGhjJ+/HiuvvpqTp48Sd++fdu+MKurqzGZTDzyyCOsXLmShIQEli5dyjPPPMMbb7wBgNFopHXkgZ07d7JhwwZmzZrFp59+yhVXXIHBYODBBx/k1VdfJTs7m+3bt/Pwww/zyRdfsXDho9z3gx/yox98n1deeaXT17Ny5Urmz5/PwIEDiYuLIy8vr+317dy5k3379pGZmcn69esvuN/ebbfdxrJly7j66qsxGo2sWbOGf/zjH0gpWb16NcHBwRQUFHDHHXdw8UgKr776KgsXLuR73/seRqOx03kZ7OVIIlgJVAE7Aft+VigOeeXOsfSNCtE6DI/VJyqYh2cO0DoMrxUWpKeq0YixxUJQF81IF6063JYEWjWazCxaddilTZ3nzp1LXJw1yd94441s2rSJq666il/84hc88cQTLFiwgOnTp7Nv3z727dvH3LnW4c/MZjPJyclt27ntttsuuL106VJmzZrFe++9x8MPP0xdXR1btmzhlltuaVuvubmZ+uYWdudu54tPVgBw991388QTT3QY65IlS1i4cCFg/WW/ZMmStkQwYcKEC770L77f6sorr2ThwoU0Nzfz5ZdfctlllxESEkJ1dTU//elP2b17N3q9niNHjlzy3MmTJ/P8889TVFTEjTfeSHZ2tl3HuDOOJIJUKeX8Xu1N6dKQZNU3rzMWi2TNoVImZMYSpVq/9EhrS5h6Y0uXieBMVceTJXW23B5Dhw5l+fLlFyyrqanh1KlTDBgwgJ07d15ylSKEYODAgezcuZPPP/+cZ599ljlz5nDDDTcwbNgwtm69dHIduHBI62uvvZann36ayspK8vLymD17NvX19URHR7N79+4Lnne8rA4hwBDQdV+LyspK1q5dS35+PkIIzGYzQggWLVp0yf47ut8qODiYmTNnsmrVKpYuXcrtt98OwMsvv0xSUhJ79uzBYrEQHHxpw4g777yTiRMn8tlnn3HVVVfx2muvMXv27C7j7oojdQRbhBAjerwnpVsWi+SVdUf5Ir+k+5X9zKGztfzwrVzWHDyndSheKyhAR4BOR31z18UIfaM7virtbLk95syZQ0NDQ9s8wmazmV/84hfcd999hIaGArB69WoqKytpbGxkxYoVTJ06lTNnzhAaGspdd93F448/zs6dOxk0aBBlZWVticBkMnU6LWV4eDjjx49n4cKFLFiwAL1eT2RkJJmZmbz//vuAtSJ71+7dNBjNjJ84ua3M/9133+1wm8uXL+fuu+/m5MmTFBYWcvr0aTIzM9m4caPDx+W2227jzTffZOPGjcyfb/2dXV1dTXJyMjqdjrfffrvDYp/jx4+TlZXFz372M6677jr27t3r8L7bcyQRTAN2CiEOCyH2CiHyhRC927tyAZ1O8MHOIt7PK9I6FI/TOhbTRFVR3GNCCPrFhdInquumt49fMeiSOQxCDHoev2JQr/b90Ucf8f7775Odnc3AgQMJDg7mhRdeaFtnwoQJ3HTTTYwcOZKbbrqJnJwc8vPz2ypFf/vb3/Lss88SGBjI8uXLeeKJJxg1ahSjR4/usiXObbfdxjvvvHNBkdG7777L66+/zqhRoxg2bBgffrQCi5Qs+t+XeOWVVxgxYgTFxR1XkC9ZsoQbbrjhgmU33XRTp62HujJv3jw2bNjA5ZdfTmCgtXnvww8/zOLFixk1ahSHDh3q8Ipi2bJlDB8+nNGjR7Nv374Op+t0RLfDULetKESH871dPFmNO/jSMNQXe+rDfD7dc4bdv5nnF1Pk2euht3M5UFLDxl/2/PLXn3U0JHFXVuwqZtGqw5ypaqRvdAiPXzHI54dCMbaY0et0PvO5c2QYakfqCE4B3wOypJT/LYRIB/oAbk8EvmxSVixLvj3FgTM1jEiN0jocj2CxSL49UcmcIUlah+L1LFJSWW8kOEBHeHDndS3Xj0nx+S/+iwV2UzfgyxwpGvo7MBm4w3a/Fui8fZXSI61t5FV/gu8cLavjfINJ9R9wAgGU1jRxvsGkdSgewyIlpyob/HosJkeuCCZKKccKIXYBSCnPCyHUEBNOlhQZzOA+EVSqiUTaZCeGs/qxy0iMUMNK9JYQgrCgAL/+0rtYo9FMVYORqBD/Gl+oPUdeuUkIocc2J4EQIgFQA+O4wOc/m47OR8opnUEIQbYaZK7XWjuShQUFUN1owthi9uvikFatSdGXJqp3dBIiR4qG/gJ8BCQKIZ4HNgEv2vtkIUSaEGKdEOKAEGK/EGKhbXmsEGK1EKLA9t/vJ+5VSeA7Ukp+s3IfOwrVhHi9ERwcTEVFBVLKti+87pqR+os62zwN9owv5A2klFRUVHTY/6AzdqdAKeW7Qog8YA7WosbrpZQHHYivBfiFlHKnbV6DPCHEauA+YI2U8vdCiCeBJ4GOu/P5iUajmTv/tY1rR/Xl/qmX9kj0J0dL61i89STD+kYxPkMNNtdTqampFBUVUVZWhpRQXt1IY1kAkV1UGPsDKSUl1U2EBgZwsNJ3jkVwcDCpqal2r293IhBC/EFK+QRwqINl3ZJSlgAlttu1QoiDQApwHTDTttpiYD1+nghCAvVUNZjYfLTc7xPBthPWKwE14mjvGAyGC4Y5yDKZ1WxlQNH5Bn7+VS6PzR3I5CF9tA5HM45cC83tYNmVPdmpECIDGANsB5JsSQLgLNBhG0EhxINCiFwhRG5Zme/PTjUpK5btJyoxW/x7wvHtxytIjgomPTZU61B8ikoCVqkxoXz56GVcMcx/kwDYkQiEED8WQuQDg2w9ilv/TgAO9ywWQoQDHwCPSikvGHxfWms4Ovzmk1L+U0qZI6XMSUhIcHS3XmdiZhy1TS0cLPHf+QmklGw7XsnEzNguR8tUHFdR18z3/rWNz/18OBOLn//QamXPFcF/gGuAj23/W//GSSnvcmRnQggD1iTwrpTyQ9vic0KIZNvjyUCpI9v0VWp+AqisNxIWpFf9B1wgOjSQvUXVfj33c3OLmQkvfM3b21Sf2G7rCKSU1UA133Uk6xFh/Un3OnBQSvlSu4c+Bu4Ffm/7v7I3+/EVyVEh3DEhza+LROLCg9jw+Cz1q80F9DrBhIxYth/339ZYe4uqKa8zkhAepHUomnOksjgIuAnIaP88KeV/27mJqcDdQL4QYrdt2dNYE8AyIcQDWIeruNXemHzdizeO1DoEj6Ca07rGpKw41hwqpbSmyS/ngN52rAIhYGKmaojg6MQ01UAePZiYRkq5CWuz047McXR7/qK8rpkAndB84nF3k1Jy5Z83cktOGg9M8++WU67SWvy49XgF1432r3GFwPq6B/eJJCbMvz5bHVET03iw0tomJjy/hl8tGOp3X4ZHS+s4dLaW0EDVusVVhvWNYs7gRCKCfadHrb2aW8zknTzPnRPTtQ7FI6iJaTxYYkQw/eJC/bLCePPRcgCm9o/XOBLfpdcJXr9vPLMH+9+ors0tFh68LIsrhyd3v7IfcOSnwDTgfiHEcaxFQwJri09VkO1CkzLj+HL/WSwW6Vdl5ZuPVZAWG0J6nP9WlrtLTZMJnRCEB/nPlUFksIFfzOv5RDu+xpErgvnAAKwdyxYAV9v+Ky40qX8s1Y0mDp71n/4ELWYL245XqKsBNzhd2cCY/17NJ3vOaB2KW+0tqqLRqMZaamVPh7JaIUQNsA/It/3fB+y3/VdcqLUN/dZj/lM81Ggyc1tOGleNUJftrpYaE0JcWCCbbEVx/qDJZObmV7fyx68Oax2Kx7CnH4Ea/1dDyVEh/Pn20UzwoyZuEcEGnl0wVOsw/IIQgqkD4tlwpMxvih9zC89jbLEwbYC64mzlG+Ou+rjrRqeQHBWidRhuc/hsLSazmurCXab0j6Oy3sihs7Vah+IWG4+WYdALv/px1R2VCLxAfXMLS3ec4rAffFCbTGau/dsm/ufLQ92vrDjFVNsv4y3H/KN4aPPRcsakxxDmR5Xj3VGJwAuYpeTpj/b5RYXezpPnaW6xqPGF3KhvdAh/uGkEc4f6fjPSynoj+8/UMF0VC11AJQIvEBlsYFRqlF9U6G0+Vm4dB0ddtrvVbePT6RcXpnUYLhcVYmDlT6Zy0zj7J23xByoReIlpA+LZW1RFdYNJ61BcatPRCkanRRPh5zNnuVuDsYWVu4s5VlandSgupdcJRqZG0zfaf+rc7KESgZeYlp2ARVrHR/FV1Y0m8ouqmNpfFQu5m7HFwqNLd/PpHt+dn0BKyaJVh9h56rzWoXgclQi8xOi0aMIC9T5dYRwaqOc/P5zEzePStA7F70SHBjKsbySbfbjC+FRlA6+sO8a+4mqtQ/E4qtrcSwQG6Njy1ByiQny3yMSg16lKYg1N7R/PG5tP0GBsITTQ974aNhZYk5zqP3ApdUXgRXw5CQC8tuEYe4uqtA7Db00ZEI/JLNlR6JtFJ5uPltM3KpjMeN+vFHeUSgRepKbJxA8W57Jyd7HWoThdaU0TL35xiC1+NJSGpxmfEYNBL9hzukrrUJzObJFsOVbBtOx4Nf91B3zv+s+HRQQFsLeoimCDzucmEmlNAGqgOe2EBgaw+YnZPjlb2ZmqRoICdG2d55QLqUTgRYQQTBsQz3ofHBdm09FyokIMDO0bqXUofs0XkwBAWmwo25+eg5r+umOqaMjLTB0QT2W90aeGpZZSsrGgjGkD4tH7UHLzRlUNRh5ZsovVB85pHYrTCSHU+dUJlQi8zLRs66XtpgLfaeZ3tqaJJpOFGYMStA7F70UEG9hYUMaq/We1DsVpGo1mZi5a5xdDtPSUKhryMkmRwVwzqi9x4UFah+I0yVEh7PzVXFosasRRrel11uLHDUfKkFL6RMXqt4WVFFY0EOnjre56Q10ReKG/3jGGm31srBS9ThAUoCaq9wQzBiZQVtvMwRLf6Ly4qaCMQL2OCRlq/KrOqETgpZpMZs7XG7UOo9fqm1uY9/IGvvKhoghvN2OgtYhuw5EyjSNxjm+OlJOTEUNIoPqh0RmVCLyQyWxh0otr+Ovao1qH0mvbjldw5FydX02c7ukSI4O5akQfIoK9/z0prmrk8LlaZg1K1DoUj6YSgRcy6HWMSo1m/eFSrUPptQ1HyggN1DMuI0brUJR2/v69cdw1qZ/WYfSaxSK5c2I6c4aoRNAVlQi81KxBCRwvr6ewvF7rUHplw5EyJmfFqfoBD9Ritnj9sOdpsaG8cMMIshLCtQ7Fo6lE4KVm2i51vfmq4HhZHScrGlSzUQ9ksUim/886fu/FU4Y2t5jZV1yNRfUi65ZKBF4qIz6MrPgw1h323go9i4QbxqQwe7C6bPc0Op1gVGo06w6VIqV3fpF+e6KSBX/dxIYC7/2MuIv31wb5sd9cO4y4sECtw+ixAYnhvHzbaK3DUDoxZ0giX+4/y/4zNQxPidI6HIetO1RGYICOSZlqaPPuqCsCLzZjYIJXfkDB2mz0aGmt1/7a9AczByUiBKw95J3Fj+sPlzI5K041G7WD2xKBEOINIUSpEGJfu2WxQojVQogC23/VdMRBGwvKvHJY6rWHSrn8pW/YU6Rmi/JUCRFBjEqNZo0XJoIT5fUcL69npqp/sos7rwj+Dcy/aNmTwBopZTawxnZfccC7207xwucHva5CbM3Bc8SGBTLCS69o/MXjVwzi6SsHax2Gw762DZo3d2iSxpF4B7clAinlN0DlRYuvAxbbbi8GrndXPL5i3rAkztU0k+9F87C2mC2sP1LGzEEJajRIDzd1QDwTvXD60Lsn9+M/P5hIakyo1qF4Ba3rCJKklCW222eBTtO3EOJBIUSuECK3rEy1Amg1e3Aiep3gqwPeM0TDrtNVVDWYmDNY/VrzBjtPnWdZ7mmtw3BIsEHPFDUJjd20TgRtpLXWsNPyDSnlP6WUOVLKnIQEVe7XKjo0kPEZMV41fvzqA+cw6AXTB6oPqjf4cGcRv1m5nyaTWetQ7PL1gXO8vPqI18TrCbROBOeEEMkAtv/eVyvlAeYO7UNVg4mqBu8YhO6nswew+P4JRAarYYG9wfxhyTSazF4zCN17O07zfu5pggK0/nrzHlofqY+Be2237wVWahiL17prUjrbnppDdKh39CmIDDaoy3YvMjErlqgQA6v2eX7xY6PRzKajZcwdmuQTcym4izubjy4BtgKDhBBFQogHgN8Dc4UQBcDltvuKg4IC9Oh0wiva5C/bcZp/bTyudRiKAwx6HZcPSeLrg+cwtnj25EEbC8poMlmYO7SP1qF4FXe2GrpDSpkspTRIKVOllK9LKSuklHOklNlSysullBe3KlLstO5QKVN+v5bS2iatQ+nS/208ztcHvac+Q7G6cngfhBCc8PBBDr/cd5bI4AAmZKpJaByhddGQ4iQpMSGUVDd59OX70dI6CkrruHJ4stahKA6aOSiB3GcvZ1CfCK1D6ZqA60anEKjqBxyixhryEQOTIshODOfTvSXcPTlD63A61Doh+rxhqtmotwnQW79YW4sfPbX8/aVbR3tFEamnUWnTh1w9MplvCysprfHM4qEv9pUwOi2a5KgQrUNReuDQ2Rpm/XE92094Zglu69wJnpqkPJlKBD7k6hHJSAlfeGDxUKPRTFSIgWtH9dU6FKWH0mNDOVfTzMd7zmgdyiWaTGam/mEtf1tboHUoXkklAh+SnRTBT2cNYGSq543fExKo590fTOL+qRlah6L0UGhgAPOGJfF5fonHtR5af7iMuuYWRqVFax2KV1KJwMf81xWDGJPuWYO4Sikpq20G1GW7t7t+dApVDSaP61z2WX4JMaEGJnvhuEieQCUCH3TgTA1bjpVrHUabPUXVTHzha9Z54XDGyoWmZccTGxboUUOf1zW38PWBc1w5IrmtUltxjGo15IOeWZFPQ7OZLx+d7hG/wFfuLiZAp2NsP8+6UlEcZ9DreGL+IOLCgrQOpc3n+SU0mszcNDZV61C8lkqfPuimsakcPlfLvuIarUPBbJF8sqeE2YMTiQpRYwv5gtvGp3O5B43zP29oEotuHsnY9GitQ/FaKhH4oGtG9iUwQMcHO4u0DoVNR8spr2vmutGqtZAvKTrf4DFDU0eHBnJLTppHXP16K5UIfFBUqIF5Q5NYsbuY5hZth+JduuMUsWGBzB6SqGkcinMtzyvil8v3crqyQfM43vv2lKYx+AKVCHzUTeNSMbZYOHK2TtM4fr1gGH+5fQxBAWoCcV9ya04aOgFLd2h3VWCxSP685gif7PW8fg3eRiUCH3VZdgI7nrmcERr3KegTFcy0bDXktK/pGx3CzEGJLMs9TYtZmz4F209UcrqyUVUSO4FKBD5KrxOEBQUgpdSkeMhikfxi2R62Hqtw+74V97h9fBqltc2s1ahZ8FtbC4kONXDVCDWIYW+pRODDTGYL1/5tMy99dcTt+952vIIPdhZxtqbR7ftW3GP24ERSokM4Vub+oanPVDXy1YFz3JaTRrBBFTv2lkoEPsyg15EeG8p7O07TaHTvVcE7208SGRyghpz2YQF6HWv/awY/ntnf7fuurDcyPCWKuyb1c/u+fZFKBD7unsn9qG408fEe9/UEPVXRwJf7znLHxHT1a83HtTYCaB1CxF2Gp0Sx8idTSYsNdet+fZVKBD5uQmYsg/tEsHjLSbeN0/76puPodYL7p2S6ZX+Ktl7bcIwZi9ZR1WB0y/6OnKt12778hUoEPk4Iwb1TMjhQUuO2ittBfSJ58LIs+kQFu2V/irZmDEqgwWjm3e2ub89vsUh+tmQXd72+3eX78idqrCE/cOPYFAx6HTkZ7pnH9c6J6W7Zj+IZBveJ5LKBCby5uZAfTM90aZ+R1QfPcehsLS/dOspl+/BH6orADwQF6Ll5XKrL53FtMLawLPc0TSZtezMr7vfQZVmU1zXzQZ7r6qKklPxlTQEZcaFqgiMnU4nAjyzPK+JnS3a5bPtvbi7kl8v3sv9Mtcv2oXimKf3jGJ0WzVtbC11WF7XmYCn7z9Twk1kD1HDTTqaKhvxITaOJj/ec4c6J6Uxy8gQelfVGXl1/jLlDkxjXzz1FUIrnEELw8m2jiQ0NdNngbztOVpIeG8r1Y1Jcsn1/ptKqH7lzYjpJkUE8/9lBzBbn/mr769oC6o0t/PKKQU7druI9MuPDiAo1YLFIlxQPPnXlED772TQM6mrA6dQR9SPBBj1PXzWE/OJqljhxxMbTlQ28s+0kt4xLIzspwmnbVbxPc4uZa1/ZxB9XHXbaNmuaTBw5VwtARLCa08IVVCLwM9eO6sukrFgWrTpMbZPJKdtsMJoZkhzJY3MHOmV7ivcKCtAzNDmSf28pZF+xc+qK/rS6gAV/2URpTZNTtqdcSiUCPyOE4HfXD+cvd4xx2q+rQX0iWPmTqarfgALA01cNISYskP96f0+vBzw8cKaGxVsLuSUnlcRIdX65ikoEfmhAYgQzBiYA9KqH5unKBp77eD91zS1qdiilTXRoIL+/cQSHztbylzUFPd5OXXMLjyzZSUyogcdV3ZNLqUTgx1buLmb6/6zjaGmtw89tMpl5bOlulucVqe7+yiXmDEnilnGprDtUhrHF8fkKpJQ8sXwvhRUN/PWOsUSHBrogSqWVaj7qxyZmxhEUoOOht/P46CdTibSzqMhskSx8bxe5J8/z1zvGkBqjBv5SLvXctcMI0IsedWSUEgYkhvNE2iAm93duU2flUuqKwI/1iQrmL3eM4WRFA7e9ts2uyjgpJc+uyGfV/nP8esFQrlE9PJVOhAUFEBSgp8HYws+X7eZEeffzFkgpKaluRKcTPDZ3IA9e5v4hrv2RR1wRCCHmA38G9MC/pJS/d8d+V+wqZtGqw5ypaqRvdAiPXzGorbNK62PFVY3ohcAsJSm2dQB++8l+zjdYW90IQELbeq3/dQJam+sLYf2Vk3LRfjqLZdbgBNYdKuswtt68rotN6R/PG/eN50fv5HHTq1v4cuFlhAV1flqU1TXz2d4SHp7Zn+9P63h0UUf23xMrdhXz3Mf7qWq0Hv+gAB3GFkvbezApK4b9Z2rbHo8JNXD1yGQ+3VPStgwgOsTAc9cOA+jyvb74tcCF73+oQYfJbMHUrgQkUG+dIe58g+mSbXZ0LLo6364fk8KzK/J5d/spWjvthhp0jEmPZtvx85jb9eTtah+d7dMZ51xXSmuaWXeolDUHS3n+huEsGNnxjwdji4UXPj/IR7uK+exn03p8pemK86+7bXZ1LKNCDAgBVQ0m+kaHkBEXwpZjlbS+a2GBep6/YQTABed1qEGHEIL6dnOJxIQa+M01w5zeqU64a2jiTgMQQg8cAeYCRcAO4A4p5YHOnpOTkyNzc3N7td8Vu4p56sN8Gtt1fAkx6HnxRusbcvFjrQw6gQV61SGrdT/tk05n++vsOZ3p6nV19dy9RVXknTzP/VOtX+6nKhpIiQlBrxOcqmjg28JKbhyTgk4nKCyvp19caIcVxD3dv71W7Crm8ff3YHJShzidsCaPjrZn0AkQYDLLC5b15v3v6Fh09f6HGPSMTY9i87HKXu3jYs485+xRWF7PwqW72XO6igUjk7lzYjpT+lvnsq5tMvHp3hL+tvYoxVWNfH9qJs9ePQSdzvEGCK44/7rbpj3Hsjs6269Je2pTDHrBoptH9ej1CCHypJQ5lyz3gEQwGXhOSnmF7f5TAFLKFzt7jjMSwdTfr6W46tJpFFOiQwA6fMyZUqJD2Pzk7C5j6eo5nenqdXX33FZbj1Vwx/9tQ68TRIUYqKy3Vgb/7c4xnf6ac+b+e7J9b3LxsXDFa+rueDvznLOXyWzhb2uP8uqGY4zPiOWdH0xESsmgZ7/EaLYwKi2axy7PZuagxB7vwxXnX3fb1OKc7Onr6SwReELRUApwut39ImDixSsJIR4EHgRIT+/9MMdnOnnjOlvubO33Y+8+7VnPGa+rf0IY/+/64ZyrbqK8rpnspAguH5JIv7gwt+y/J9v3Jhe/Ble8pu626cxzzl4GvY7H5g7kRzP6U15nndFMCMGvrxlKRlwYUwfE9boZsivOv+62qcU56ex9ekIisIuU8p/AP8F6RdDb7fWNDukwi/d10xVB6366iqWr53S1Tlevyx6JkcHc3cO5YJ2x/55s35tcfCxc8Zq6O97OPOccFRKov2CKSWfOO+yK86+7bWpxTjr7ffGEVkPFQFq7+6m2ZS71+BWDCLloPt0Qg57HrxjU4WOtDDqBvgdllx3tp6tYuntOZ7p6Xe7g6v0/fsUga9m9k+gEnW7PoBMY9OKSZb15/zs6Fl29/yEGPVP7Ozaaqz3H25nnnCdxxfnX3TbtOZbd0Qn7v4wNeuH098UTEsEOIFsIkSmECARuBz529U6vH5PCizeOICU6BIG1zK218qf9Y2CtTMS2zqJbRvG/t4wiJvS7NvetXwut67X+b/990XrF234/XcVy16T0DmPrzetyB1fv//oxKSy6ZRTRId8d/6AA3QXvwdT+sRc8HhNq4K5J6RcsA2uroZduHc2iW0Z1+l4vunnUBa+lo/c/1KDDcNEnKVAv2tZpv82OjkVX59uLN47g3R9O5q5J6bQvNQk16JjaP7Zt3Vb2Hm9nnnOexBXnX3fb7O5YRocYiAk1tD02tX8s7d+1sEA9L906mpduG33BORpq0BEWeGGCiQk19LiiuCuaVxYDCCGuAv6EtfnoG1LK57ta3xmVxYqiKP7GkyuLkVJ+DnyudRyKoij+yBOKhhRFURQNqUSgKIri51QiUBRF8XMqESiKovg5j2g15CghRBlwsodPjwfKnRiOs6i4HKPicoyKyzG+Glc/KWXCxQu9MhH0hhAit6PmU1pTcTlGxeUYFZdj/C0uVTSkKIri51QiUBRF8XP+mAj+qXUAnVBxOUbF5RgVl2P8Ki6/qyNQFEVRLuSPVwSKoihKOyoRKIqi+DmfSgRCiPlCiMNCiKNCiCc7eDxICLHU9vh2IURGu8eesi0/LIS4ws1x/VwIcUAIsVcIsUYI0a/dY2YhxG7bn1OH57YjrvuEEGXt9v+Ddo/dK4QosP3d6+a4Xm4X0xEhRFW7x1xyvIQQbwghSoUQ+zp5XAgh/mKLea8QYmy7x1x5rLqL63u2ePKFEFuEEKPaPVZoW75bCOHU4XztiGumEKK63Xv163aPdfn+uziux9vFtM92PsXaHnPl8UoTQqyzfQ/sF0Is7GAd151jUkqf+MM6hPUxIAsIBPYAQy9a52HgVdvt24GltttDbesHAZm27ejdGNcsINR2+8etcdnu12l4vO4D/tbBc2OB47b/MbbbMe6K66L1H8E6dLmrj9dlwFhgXyePXwV8gXV6iknAdlcfKzvjmtK6P+DK1rhs9wuBeI2O10zg096+/86O66J1rwHWuul4JQNjbbcjgCMdfB5ddo750hXBBOColPK4lNIIvAdcd9E61wGLbbeXA3OEEMK2/D0pZbOU8gRw1LY9t8QlpVwnpWyw3d2GdZY2V7PneHXmCmC1lLJSSnkeWA3M1yiuO4AlTtp3p6SU3wCVXaxyHfCWtNoGRAshknHtseo2LinlFtt+wX3nlj3HqzO9OS+dHZdbzi0AKWWJlHKn7XYtcBDrfO7tuewc86VEkAKcbne/iEsPZNs6UsoWoBqIs/O5royrvQewZv1WwUKIXCHENiHE9U6KyZG4brJdhi4XQrROKeoRx8tWhJYJrG232FXHqzudxe3KY+Woi88tCXwlhMgTQjyoQTyThRB7hBBfCCGG2ZZ5xPESQoRi/TL9oN1itxwvYS2yHgNsv+ghl51jHjExjWIlhLgLyAFmtFvcT0pZLITIAtYKIfKllMfcFNInwBIpZbMQ4iGsV1Oz3bRve9wOLJdSmtst0/J4eSwhxCysiWBau8XTbMcqEVgthDhk+8XsDjuxvld1wjpD4Qog2037tsc1wGYpZfurB5cfLyFEONbk86iUssaZ2+6KL10RFANp7e6n2pZ1uI4QIgCIAirsfK4r40IIcTnwDHCtlLK5dbmUstj2/ziwHusvBbfEJaWsaBfLv4Bx9j7XlXG1czsXXbq78Hh1p7O4XXms7CKEGIn1/btOSlnRurzdsSoFPsJ5xaHdklLWSCnrbLc/BwxCiHg84HjZdHVuueR4CSEMWJPAu1LKDztYxXXnmCsqPrT4w3p1cxxrUUFrJdOwi9b5CRdWFi+z3R7GhZXFx3FeZbE9cY3BWkGWfdHyGCDIdjseKMBJFWd2xpXc7vYNwDb5XeXUCVt8Mbbbse6Ky7beYKyVd8Idx8u2zQw6r/y8mgsr8r519bGyM650rHVeUy5aHgZEtLu9BZjvxrj6tL53WL9QT9mOnV3vv6visj0ehbUeIcxdx8v22t8C/tTFOi47x5x2cD3hD2ut+hGsX6rP2Jb9N9Zf2QDBwPu2D8a3QFa75z5je95h4Eo3x/U1cA7Ybfv72LZ8CpBv+zDkAw+4Oa4Xgf22/a8DBrd77vdtx/EocL8747Ldfw74/UXPc9nxwvrrsAQwYS2DfQD4EfAj2+MCeMUWcz6Q46Zj1V1c/wLOtzu3cm3Ls2zHaY/tPX7GzXH9tN25tY12iaqj999dcdnWuQ9r45H2z3P18ZqGtQ5ib7v36ip3nWNqiAlFURQ/50t1BIqiKEoPqESgKIri51QiUBRF8XMqESiKovg5lQgURVH8nEoEitIFIURcu9Eozwohim2364QQf9c6PkVxBtV8VFHsJIR4Duvopn/UOhZFcSZ1RaAoPWAbT/9T2+3nhBCLhRAbhRAnhRA3CiH+xzZ2/Ze2oQMQQowTQmywDVq2yjZypKJoTiUCRXGO/lgH5LsWeAdYJ6UcATQCV9uSwV+Bm6WU44A3gOe1ClZR2lOjjyqKc3whpTQJIfKxTq7ypW15PtaxbQYBw7GOWoltnRIN4lSUS6hEoCjO0QwgpbQIIUzyu8o3C9bPmQD2SyknaxWgonRGFQ0pinscBhKEEJPBOuRwu8lYFEVTKhEoihtI67SLNwN/EELswTq65BRNg1IUG9V8VFEUxc+pKwJFURQ/pxKBoiiKn1OJQFEUxc+pRKAoiuLnVCJQFEXxcyoRKIqi+DmVCBRFUfzc/wfGV1YbBwgRUQAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "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 GP\n",
    "\n",
    "When using a GP 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 GP can be negative\n",
    "2. The intensity function can take on large values (e.g. 900), whereas GP 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 GP 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/GPyTorch interface\n",
    "\n",
    "Unfortunately, this is not a likelihood function that is easily written as a GPyTorch Likelihood object.\n",
    "Because of this, we will be using the low-level interface for the Pyro/GPyTorch 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 `gpytorch.models.ApproximateGP` object to model the non-homogeneous intensity function. This object needs to define 3 functions:\n",
    "\n",
    "- `forward(times)` - which computes the prior GP mean and covariance at the supplied times.\n",
    "- `guide(arrival_times, quadrature_times)` - which defines the approximate GP posterior at both arrival times and quadrature times.\n",
    "- `model(arrival_times, quadrature_times)` - which does the following 3 things\n",
    "\n",
    "  - Computes the GP prior at arrival time and quadrature times\n",
    "  - Converts GP 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 GPyTorch 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": [
    "class GPModel(gpytorch.models.ApproximateGP):\n",
    "    def __init__(self, num_arrivals, max_time, num_inducing=32, name_prefix=\"cox_gp_model\"):\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 GP\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 = gpytorch.variational.CholeskyVariationalDistribution(num_inducing_points=num_inducing)\n",
    "        variational_strategy = gpytorch.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 = gpytorch.means.ZeroMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "    \n",
    "    def forward(self, times):\n",
    "        mean = self.mean_module(times)\n",
    "        covar = self.covar_module(times)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean, covar)\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 = GPModel(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,
   "outputs": [],
   "source": [
    "pyro.clear_param_store()  # so that notebook supports repeated runs"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "  0%|          | 0/200 [00:00<?, ?it/s]",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "416f73ab71c84e2c8f247b26feceef51"
      },
      "application/json": {
       "n": 0,
       "total": 200,
       "elapsed": 0.012999296188354492,
       "ncols": null,
       "nrows": null,
       "prefix": "",
       "ascii": false,
       "unit": "it",
       "unit_scale": false,
       "rate": null,
       "bar_format": null,
       "postfix": null,
       "unit_divisor": 1000,
       "initial": 0,
       "colour": null
      }
     },
     "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 GP mean, plus the intensity function at the 5th and 95th percentiles of the GP."
   ]
  },
  {
   "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": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABzWUlEQVR4nO2dd5xcVfm4nzO9bO89vfdCSAgl9I5UQQEFkSIKCspPsSD6xQqKgAqi0hTpiPTekxBI72XTtvedmZ3ezu+PmZ3sbN9kZ7ad5/OZZOfeM/e+c+fe+97zViGlRKFQKBRjF81QC6BQKBSKoUUpAoVCoRjjKEWgUCgUYxylCBQKhWKMoxSBQqFQjHF0Qy3A4ZCTkyPHjx8/1GIoFArFiGLdunVNUsrczstHpCIYP348a9euHWoxFAqFYkQhhDjY3XJlGlIoFIoxjlIECoVCMcZRikChUCjGOCPSR6BQKBSBQICqqiq8Xu9QizLsMJlMlJSUoNfr+zVeKQKFQjEiqaqqIjU1lfHjxyOEGGpxhg1SSpqbm6mqqmLChAn9+owyDSkUihGJ1+slOztbKYFOCCHIzs4e0ExJKQKFQjFiUUqgewZ6XJQiUCgUijGOUgQKhUIxQvD4g/iD4UHfrlIECoVCMUJw+UMEw0oRKBQKxbBCq9Uyf/58Zs+ezSWXXILb7T7sbd15553cc889XZZv2bKFcePG8Y+/PXQkovaIUgQKhUJxBJjNZjZu3MjWrVsxGAw89FD8zVpKSfgIn+LnzJnDk/95iuee/s8RbacnlCJQKBSKQeK4446jvLycAwcOMG3aNL72ta8xe/ZsKisr+fe//82SJUuYP38+119/PaFQCIBf/epXTJ06lWOPPZZdu3b1uO2snBx279yRELlVQplCoRjx/OKVbWyvcQzqNmcWpfHzc2f1e3wwGOSNN97gjDPOAGDPnj08/vjjLF26lB07dvDMM8+wcuVK9Ho9N954I08++SSzZs3i6aefZuPGjQSDQRYuXMiiRYu63f5Pfvxj/H4fFQcPMn3KpEH5ju0oRaBQKBRHgMfjYf78+UBkRnDNNddQU1PDuHHjWLp0KQDvvfce69at46ijjop9Ji8vj5aWFi644AIsFgsA5513Xrf7eOONN3A6naw49XS2b9+uFIFCoVB0ZiBP7oNNu4+gM1arNfa3lJKvf/3r/OY3v4kb86c//anP7Xu9Xn74wx/yzyef499PPMb2bdu48EvnHqnYcSgfgUKhUCSYk08+meeff56GhgYAWlpaOHjwIMcffzwvvfQSHo+HtrY2XnnllS6fveuuu7jyyispLC5lxqzZbN++ddDlUzMChUKhSDAzZ87krrvu4rTTTiMcDqPX6/nLX/7C0qVLufTSS5k3bx55eXkx01E7u3bt4p133uHDjz+hyRVk+oxZPHRf1/DSI0VIKQd9o4lm8eLFUrWqVCjGNjt27GDGjBlDLUZS8ARCNLX5AMhOMWAx9P0M393xEUKsk1Iu7jxWmYYUCoVimBMMDX42cUeUIlAoFIphji8YRpPASqtKESgUCsUwxx8Mo0ng3VopAoVCoRjGhMOSUFgiUDMChUKhGJMEw2EEkMgePElVBEKIW4QQ24QQW4UQTwkhTEKICUKINUKIciHEM0IIQzJlUigUiuFMMCxJdGxn0hSBEKIYuBlYLKWcDWiBy4DfAfdKKScDrcA1yZJJoVAohjv+YDiBRqEIyTYN6QCzEEIHWIBa4CTg+ej6x4HzkyyTQqFQDFsijuLEqoKkKQIpZTVwD1BBRAHYgXWATUoZjA6rAoq7+7wQ4johxFohxNrGxsZkiKxQKBQ90tzczPz585k/fz4FBQUUFxfH3vv9/kHZh5QSfyhMZz3Q3qjmwQcfHJT9JNM0lAl8CZgAFAFW4Iz+fl5K+bCUcrGUcnFubm6CpFQoFIr+kZ2dzcaNG9m4cSM33HADt9xyS+y9wXDI1XkkjWlCYYmUIDp5iufMmcPTTz/NE088cUTfoZ1kmoZOAfZLKRullAHgRWA5kBE1FQGUANVJlEmh6Bf+YJjGNh976tvYUmXH7Q/2/SHFmKVzY5pPPvmE2bNnx9bfc8893HnnnQA9NqyBiKO4J6NQXl4e27ZtGxR5k1l0rgJYKoSwAB7gZGAt8AFwMfA08HXgf0mUSaHokXBY8uKGKrbXOGjzBmm/IsNh0GsFx03JYenEbLJTjEMrqAKAFStWdFn25S9/mRtvvBG3281ZZ53VZf1VV13FVVddRVNTExdffHHcug8//PCI5OnYmObAgQPdjumpYc3XvvY1IFJaoqeIoR/96Ef4fD4OHjzIuHHjjkjWpCkCKeUaIcTzwHogCGwAHgZeA54WQtwVXfbPZMmkUPTG2oOtrCpvJi/NSGG6KW56HgiF+WhXIx/sbGReaTrLJ+dQkmnBoFOpOYoIHRvT9ERPDWva8XXjH4BIoxqXy8XZZ5/Ntm3bRo4iAJBS/hz4eafF+4AlyZRDoegLm9vPfzdUk5tqxKjTdlmv12oozDATDkt21DrYXGVHCJiQY2VmYRrjsq0UZ5rRa5ViSBa9PcFbLJZe1+fk5BzxDKAzHRvT6HS6OD+B1+sFem5Y006gmxpD7Y1qXn75ZR599FG2bt3a7WxnIKizVKHohJSSlzZUI6XEpO+qBDqi0QhyU00UZZjJSzVR7/Dx6uZaHvhgD4+tOsBILPOuGHzy8/NpaGigubkZn8/Hq6++CvTcsAYgLCWBkOziKP7db37N1772NcaPH8+cOXPYuvXIG9UoRaBQdGJrtZ3NVXbyUgdm+9dqBOlmPUUZZorTzeyocbC9dnAbqitGJnq9njvuuIMlS5Zw6qmnMn36dCC+Yc3cuXM59dRTqa2tBeDss86irrY2rrTE3j27ef+9d/ne974HMGiKQDWmUSg64PQF+d2bOzFqNViNR2Y5bfMG0AjBbWdM69a8pDgyRntjGo8/SJPTH+d38gfDqjHNYLDuYEvCmzwoRiZSSl7dVIMvEDpiJQCQatLT6vazqrx5EKRTjDU8gXBCC811ZMwpgtc217K12j7UYiiGIXsanKzZ30J+qmnQtpmXauKtbXW0ugYn01QxNpBS4vEH0Sa4tEQ7Y04RePwhXttSq2YFijjCYcl/11eRYdYPal0Xg06DlPDG1tpB26biECPRtN0fAiFJWHLYXckGelzGnCIAaGzzsaVKzQoUh9hV30Z9m480s37Qt52XamTtwVYONLkGfdtjGZPJRHNz86hUBr5gqO9BPSClpLm5GZOp/zPbpOYRDBcyLQZe21LL7JJ0FeetQErJW9vqSB0Ev0B3aDSCFKOOF9dX8d1TpiZtuj/aKSkpoaqqitFYhNLhDRAO0yWZLBSWtBh1fSYumkwmSkpK+r2/MakIrEYdtXYPmyttLBqfNdTiKIaYfU0uKlrclGSYE7aPDLOeapuHdQdbWTJBnXODgV6vZ8KECUMtxqDT5g3wi1e2UZRu7pJDUGPzcOWyYmaUZAzqPsfs43CmxcDrW2sJKF/BmEZKyTvb67HotV0uusFECEG21cgrm2tUwTpFrxxsdoMUCT0fOzNmFYHVqMPmDrCp0jbUoiiGkKpWD3vq28iyJr5DqtmgxeMP8fHu0WfKUAweW6rtSa9ZNWYVAURnBVtq8QfVrGCs8v7OBow6TdKevvJSjby3o4Empy8p+1OMLEJhydZqOxmWwQ9a6I0xrQisRh12T5ANla1DLYpiCKh3eNlcZUtqGWm9VoNWI3hjiwonVXSlxubBHwwnPYhlTCsCgCyrnje31ClfwRjkw12N6LWaw47VPlxyU41sqLCxX4WTKjqxp75tSPY75hWBxaDD4Q2wTWUbjymanT7WHmghZwiaymhEJJz0v+urCIVHXwy84vDZUGkjzZRcsxAoRQBAmlnPO9vrCauLcsyw7mDEHDhUMf0ZFj1VrZEQZoUCwO4OUGvzYjUmv0ChUgRAmklPncPL3kbnUIuiSAJSStYebE26Q64jQgiyrAZe2liNN3D4WaSK0cOBZheIrv0HkoFSBFEsBh3vbK8flenqinianH5aXH7MfTSdSTRWow6nL8jmKtuQyqEYHmyusvdZrjxR9yelCKJkWvTsbXRS1eoZalEUCaa8ITLzG4onr86kmfR8sqdJPYCMcYKhMDtqHWT0Uuuqzu7lre31VLW6B33/ShFEEUJg0Gn4cFfDUIuiSDDrD7aSkqC6QgMl1aSjxuahzuEdalEUQ0hlq4dAKIyuh7BRbyDE61trCYUlGebBT35UiqADOVYjmyptNLapZJ/RitMXZH+zi1TT8FAEQgi0GsH6gyqXZSyzucqGpoe7sZSSd3fU4/IFWToxi5QEnLtKEXRAoxFoNIKVe5uGWhRFgtjf6IIjqPOeCLKsBlbva1a5LGOUUFiy9kALmZbun/Q3VdnZ2+hi+aQcsq2JCXdWiqATuSlGVpc30eYNDLUoigSwqcqGUT+8TnujTosnEI75LhRji4oWN95AuFtHcb3Dy6d7mpiQY2VBWUbCZBheV8QwQKfVEJLw+f6WoRZFMcgEQmG21fTukBsqzDoNq/eq3sZjkZ7MQr5giDe21mE2aDl1Zn5CgxuUIuiGHKuBD3Y1qGJ0o4yKFjfBXhxyfdHi8vPR7kZWljcNeux/psXAjlo7do+aiY4lgqFwt2YhKSXv7WjA4Q1w5uyChIc6Dw+P2TDDqNfS5PSzq87BnEFuAKEYOrbXOBjoQ5WUkooWNxsqbRxsdqMVgpCUbKm2s3h8JvNLMg5bsXREoxFIBFuqbBw7JfeIt6cYGbSbhbKs8Tf6apuHPQ1Olk3KpiiBDZPaUYqgB6xGLR/uamR2cfqwiDdXHBlSSjZUtPbokOuO6lYPH+xuoNkZST47ekIWc0vScftDrCxvYmV5M5sq7SybmM30wtQjdkBnmCM5Bcsn56hzbowQMQt1/a231zowaDUsKM1IihxKEfRAulnPgWYXtXZvUjSyIrHUObw4vAGKMyz9Gu/xh3htSy16reCUGXlMy0+NPflbDDq+NL+YqlY3n5Y38c6OenbXt3HWnMIjaihiMWipsXuoavVQmtU/ORUjl2AozNqDrWR1ejjxByOBA1PzU5NWjlr5CHpACIFOo+GzfcqBNxrYXd/GQJ6xP9zdgC8Y4tx5RcwqSu/W/FOSaeHSxaWcOC2XihY3L6yvOqI2lO05BWsPqECFscDBFjfeQKjLw8PeRieBkGRGQVrSZFGKoBdyUgx8vr9F9ZgdBaw/aCO1n+V9yxuc7K53cvSE7D7LVAshmFuSwTlzC2lx+Xl2bdUROXyzrUY+P9CiCtGNAbZU2dB1Ey60vdZBullPUYYpabIoRdALOq2GYDjM5krVq2AkY3cHqG719KushMcf4v2dDeSlGlk0LrPf+5iYm8KFC4vxBUI8u7aShsMsGaHXagiEpMopGOW0m4U6+6wcngBVrR5mFKSq5vXDiQyzgfd3NaheBSOYgZT3bTcJnTozf8C9CgrTzVyyuBStRvD8+irq7IenDEw6DV8o89Copiez0M66SIeyGYXJMwuBUgR9YjXqaHb62KfaCo5YdtQ6MGj7jsMeiEmoJ7KsBr68qBSTXsub2+oOKxclw2Jge40Dp0+ZJEcrmyq7moWklOyodVCSYSYtyUmPShH0A6NOy8pyVX9oJCKlZEedgzRz72ahwzUJdUeKScfpMwuwewJ8sqdxwJ9vn4nsrnMckRyK4UkwFGZdN2ahWrsXmyeQ9NkAKEXQL7KtBrZU2Wl1+YdaFMUAaXL6cflCfTb8WHuwBW8wxCkzBm4S6o7iTDOLxmWytcbBvqaB2/stBi2r9ynz0Gik1u7FHwx3MQvtqHWg0wgm56UkXaakKgIhRIYQ4nkhxE4hxA4hxDIhRJYQ4h0hxJ7o/0f2OJYAIlVJ4XNltx1xVLa4oI+mL+01iCbnppCbOnjVHZdOzCInxcC72xsGHHmWbtazv8mlHj5GIQe6MTMHQ2F21zuZkpdyRLkoh0uy93gf8KaUcjowD9gB/Ah4T0o5BXgv+n7YkW018uGuBlzKbjui2F7bhrGPOi07a9vwBcPMH+QsTp1Gw+mzCvAHw7y/s2FAXcgiju2IzVgxuthW48BiiD8n9za68IfCQ2IWgiQqAiFEOnA88E8AKaVfSmkDvgQ8Hh32OHB+smQaCAZdJKxPVSUdOYTDkp11DtJ6yR+QUrKpykZuqpHC9MGP285JMXLM5Gz2NrrYPsCbeppJz8q9zaqN5SjCHwyzr8nZJadlR62DVJOOksyhqWKQzBnBBKAReFQIsUEI8Q8hhBXIl1LWRsfUAflJlGlA5FgNvLezXiX7jBAanT68ga622I5UtXpodvmZX5KRsLjtBaUZlGSa+Wh3I44B9LlIMeqod3hoUB3zRg01Ng9SEueH8gZCVLS4mZaf3NyBjiRTEeiAhcCDUsoFgItOZiAZefTp9vFHCHGdEGKtEGJtY+PAIzEGA6Nei9cfZt1BNSsYCVQ09x3yu7HShlmvZWp+4hx0QghOnZFPWDKgkiVCCASCLdW2hMmmSC77uzknK1vdSGBCjjX5AkVJpiKoAqqklGui758nohjqhRCFANH/u+0eL6V8WEq5WEq5ODd36Mr0ZlkNvL2tXvUqGAHsqG3D1MtswO4JsL/JxezitEEpJd0baWY980rS2VHbRpOz/0/4GRY9q/c2q4TGUcLWKnuXDPfKFg8GrYb8tOSVlOhM0hSBlLIOqBRCTIsuOhnYDrwMfD267OvA/5Il0+FgNmhx+kJsqrQNtSiKXgiFJbvq23pNzNlcZQMBc4rTkyLT4vFZGHQaVg2gE5nFoMPuCVJt8yRQMkUyaDcBdW4+X9HipjjTPChhy4dLsqOGbgKeFEJsBuYDvwZ+C5wqhNgDnBJ9P6zJtOh5c1stQdVsfNhS7/DiD4V7LOPbMWS0v8XojhSzXsvicZnsb3JR3dr/G7tWCDaqB48RT1X0N+/Yt8LuCWD3BCgb4rLjSVUEUsqNUfPOXCnl+VLKVills5TyZCnlFCnlKVLKYW+Atxp1tLoCbK1WxeiGKxXNbmQvenpnXSRkdF6SGn+0M780A6tRy8q9Tf2OBsq06vl8fzMhZR4a0exrdHZxBle2uAEoHaJooXZUZvFhkmHR8+a2OnVxDlO213aN1W5HSsmmykjIaFECQkZ7Q6/VsHRCNrV2b7/rVxl1WryBMHWHWdFUMTzYUm0n1djVLGQ1asmy9r9zXiJQiuAwSTXpaWjzsb1GzQqGG8FQmD0NbaSauq8vVGPz0uzyM69kaNqQzixMI9OiZ1X5QJzAkn2NqjT1SMXlC1Jr92I1Hno4kVJS2eqmLNMy5K1JB6wIhBBWIUTfpRzHABlmPa9uVr6C4Uadw0swJHuMBNrT0IZOI5iSl5pkySJoNIJjJuXQ4vazvZ+F5axGPVuq1EPHSKWq1YMQxN3w2/NcempL6gmEaBtA3smR0KciEEJohBBfFUK8JoRoAHYCtUKI7UKIu4UQkxMv5vAk1aSnyeljg3LkDSsONrt7XCelpLzRybhsy5DUdGlnUq6VgjQTa/a19OtBItWoY3+zSyUzjlDKG9roHBRUEfUP9OQobnb6cHiD+JLwm/fnSvgAmATcDhRIKUullHnAscBnwO+EEFckUMZhTZbVwKubatQFOozYVmPv0T9Qa/fi8oWGpMJjR4QQLJ2YhdMXZFd9W5/jNRqBlIeci4qRxdYaB6nG+Oi0yhYP2VYD1m465zl9QbKsBi5eWEJ9m5dwgsuM9EcRnCKl/D8p5WYpD8VhSClbpJQvSCkvAp5JnIjDG4tBh9MXVDWIhgmBUJi9ja4eQ0LLG5xohRjSLM52yrIs5KQY2FBh61cEkVYDu/uhNBTDC4c3QGObN+7hJBgKU23z9GgWsrn9nDGrgKMnZrGoLIv6BAcK9KkIpJR9Gqn6M2Y0k5ti5M1tdaoy6TCgzu4lHJbdJue0m4XKsi199ifoCSkldk9gUPxCQggWlmXS7PJzsB9P+mkmPZuq7KoI3QijqiWSP9DRP1Bj9xIKy27NQh5/iFSTjjnR+lcXLCwm1aTD7kncbXZARlIhRKkQ4gwhxA+EEI8LIdYmSrCRhFGvJRAM86nqYjbkVLd6ui9WBdS3+WjzBg/bLBQOS6ptHqxGLc0uP9U2dyRx7QjKjUzNT8Vq1LK+orXPsWa9lhanH5t7TD93jTh217eh7dSWsrLFjUZAcUbX/IEWl5/TZhbEfFhWo46vLRtPmzdAIJyYwJT+OIuvF0KsEkLYgN3AN4EUIqUhvpoQqUYgualG3t/RgM2tGokMJbvr2zD30H+gvMGJRsDEwzALBUJhqmwejpmUzS2nTOWXX5rN9cdPYvG4TBzeAJUt7sOaJWg1gvklGVS2eGjso8qoEAIE/Zo9KIYP22scpHcyVVa0uClIM3UJWPAFQpgMGhZ2apc6LtvK2XMLcXoTY3Xoz4zgduAWYBHwKmACHon6B3YnRKoRiF6rQUrJ+zu7rZmnSAJSSvY0OLsU9WpfV97gpDTLgqmPRjWd8QRC1No9nDevkAsXlqDTajDoNEzJT+XixaX84rxZnD23kBq797CCBmYXp6PXin7NCow6jcpoH0E4fUFa3D5M+kO3Wm8gREObr1uzUJPLz8nT87s9R0+YmseySdk9lk05EvqzxXOklGuklHullJcAfwFeEULcIoRQCWkdyEszsXpvM80DqC6pGDyaXX68gVC3YaFNTj92T4DJuQMzCzm8AVpdfq46ZgInTs/vNvFHp9Vw8ox8rlw6jmaXH+cAfUUmvZZZhensrm/r84kv3axnR61DZbSPEOrsXgQi7ryJlZXopAj8wTA6jWDJhKxut6XVCK4+ZgLT8gc//6U/zuKtnd6/ASwBsoCVgy7RCKbdQakKhA0NtTYvsgcPwZ6GNoSASQNQBB5/CLcvxLdPnNyvmkQLx2VywwmTcPuDtA7QRDi/LAMpYWOVrddxeq0GfzBMrV1VIx0J1NjcXc7IilZ3t2Wnm1w+VkzL7TactJ1I//TBz0Luj4+gy16llD4p5c+Ilo/ubsxYJdNiYGV5k6ofPwTsbXSi03Q9pdtNRiUZZsw95Bd0JhSWNLb5+MqSUsYPwKcwOS+Fm0+egk4jaHb2Xxmkm/VMzktha7W9X87n/Y39q1OkGFr2NDi75LRUtXi6lJ0OhsIIYNmknCRLGKFfCWVCiJuEEGUdFwohDECJEOJxDvUTGPOYDVrs3gAH+tEdSzG47K5v69Y/0OyKRNoMJFqozuHl2CnZh1WdtDDdzHdOmgLIAYUULyzLxBcMs62P+lUpRh2b+pg5KIYeKSX7m1xx56TbH8TmCXSJFmp1B5hfmkF6L/0zEkl/FMEZQAh4SghREy0tsQ/YA3wF+JOU8rEEyjjiMGg1rDvYt+NPMXi4/UEa2nzdZhSXN0SKtfXXLNTi8lGQZuSceUWHXQwsy2rgK0vKaHb5+m3PL0g3UZhuYmNl7wlmKUYdB5vdePwqm3040+oO4AvE98SotUcSwwo7Vb31B8MJsf33l/74CLxSyr9KKZcD44h0FlsopRwnpbxWSrkh4VKOMLIsBtYebFFlJ5JIjc2LgG5v3OUNToozzL3aXtvxBkL4Q5Irl40/7KSzdmYWpbF8cs6AykfPL83A4Q32GiKq0QgkkV63iuFLnb3r715r96IRkJdqjFsuNFA0hD0J+r4yOhDNIK5NkCxJ4d1Hf09x6Tgmzl5E0aTpaLUDOgT9QqfVEAzBzloH88sy+/6A4ojp6aZoc/tpdvk5YWrffa7DYUlDm48rlo4bcP/Yzz//nA8//DD2XqvVcvbZZ3PO3CmUNzhpdfvJtPRdc35SbgpmvZat1XbGZ/fsm9BqBLvq2pg6hE+Rit6pbHXT+bmk1uYhL9UUVxk3GI5EC+Wldn/O1dbWsnLlSlatWsXFF1/MMcccM+iyDv5dcBjjdrvZveY91r5WB4DBaKZs+lyWn3c58447fVD3FelC1awUQZLYVdvW7RP//mjzl/7UFqpr87J0YjYLyzL6HNvc3MyDDz7I1VdfTXFxMZ988gk//OEP48b8v//3/6ioqOCKpeP407t7sBrCfVY81WoEM4vSWF/RitMb7NLftp00k44tVXbOmVs45LXsFd1T3tCG1XDo9wuFJfVtPuaWxPfIdnqDjM+2xjmPpZQIIZBSMnfuXJqamjCZTEybNi0himBM5QFYLBZufOgtfvKvD7jy9j+y5IyLcLQ0QgJqt2SY9exvcqmcgiQQDIU50OyKa/rRzoFmN5kWfZ9OOF8whF4jOHde7zfWiooKbr75ZsrKyvjZz37GG2+8AcDNN9+My+WKvaqrq3n88ccpLi6mJNNC44f/4tN3XyfcjxIBs4vSkJJencZmvZYWtz+h9WcUh08oLKlo9sQ9nDS2RfxFnf0DLn+IGYWRmZ2UkkcffZSlS5fi8XgQQvDwww+zevVq7HY7119/fULk7feMQAhxE/BvKeWI94Jm5BaSfWIRC048m3A4jCYacrj6tWewpmcy99jTjngfQggEsLnKzonT8454e4qeaYheYJ1DRwOhMNWtHuaWpvfwyUM0Of2cObsAi6HnS2LVqlWcc845OJ1OvvrVr/KDH/yA2bNnA6DX69HrDykbi8XCFVdEqrN7PB5Wvf0/ysvLqd76GZff8gs02p79DxkWA2VZFrbWODhqfFa3cePtyqqyxUNGP0xOiuTS5PQRlvHFD2uiuR+F6fG+AAGUZFqora3luuuu49VXX+WEE07AbrdjNpu54IILEi7vQGYE+cAXQohno4XnRsV8tF0JhMNh1r73Px775U38+zffx+PsX+eo3siw6Pm0vFHlFCSYGlv3heYqW9yEpOzV1g6RiA29VrB0Ynav4x566CFycnLYvn07jz32WEwJ9IXZbGbnzp1853vfZ8Pbz/Ho/92M39e7A3lOcTpOX7DXMGS9VrCrnx3OFMmlzt41ubHW7iXNpIsLJ5UyMuqz915n9uzZvPvuu/zpT3/i/fffp6CgIGny9lsRSCl/CkwB/glcBewRQvxaCDEpQbIlFY1Gw42/f5wzvnYzGz9+k7//9Dr83iPL3rQYdNg9wVgnIkVi2FPfhrEb2/uBZjd6raAoo3fHb6PTy0nT83qMKnK7I7/fww8/zKpVq5g8eeBN+bRaLQ/cew9X/+AXbFv9Ho/8/MZeQ0Qn5FixGLRs6aWuULpJz7YahypLPQzZ3+SKm6FKKam1eyjslD/g9ofItWj4w92/Y8qUKWzcuJHvfve7sQfUZDGgvcnIGVcXfQWBTOB5IcTvEyBb0tHq9Jx2xbe58sd/5OCOjfzr17f0y6bbGzqNUDkFCaQ9azi10028PZmnLMvSbbZxO/5gGL1GwzHdZHRKKfn5z3/OUUcdhc1mw2QykZNzZJmff7jzR5x7y+845rzLe/VFaDWC2UXpHGh24+jBD2DUa2nzBWl2qYq3w43yTsUP27xBXL5QF/9AmzfIzJJsPvjgAz766COmTZuWbFGBASgCIcR3hRDrgN8TqTE0R0r5LSJVSS9KkHxDwrzjTufC79zB1IXLj1gzZ1kNfHGgBV9Q5RQkArsngMMb7BKN0178rS+zUKPTx0kzup8N3H777fzyl79k6dKlpKQMTmvLTKuB675+ObmzIpEfW1a+i6O5+4q1s4rSANhW05v5R6r2lcMMbyBEvcMbV86kPZGsqIN/4OCOTfzv3v9HYYogLS0No9HYZVvJYiB3uSzgQinl6VLK59q7kkXbV56TEOmGkOXnfpXjzr8SAFtT/WFvR6/VEAyH2adqwySEnhLJ2m3rvSmC9mqP3dV3eeGFF/jd737HDTfcwD/+8Q90usGLtD5+ai4aIWhpauTJ393GE7++lVCoaymKNLOe8dkWttXYe8xONuq0bK9VfoLhRHtfCU2Hc7LW7kGvFWRbI479+oq9/P2n11G3dxupYug7Gw5EEZiklAc7LhBC/A5ASrljUKUaRtTu381vv3EGH//38cPehkGrVeahBHGg2dUlaQfgQJObnBRDj3H4EJ0NTM/rUp9o9+7dXH311Rx99NHcd999gx6nn2bSc/KMPDy6FC757i/Yt+ULXn/kj92OnVOcjssfiuVDdLetHbUOFZAwjKi1dw1eqLF7KUgzodEIbE31PPzjb6LRafn6Lx5m6oSSIZGzIwNRBKd2s+zMwRJkuJJfNolpi4/lpQd/zY7PPzqsbWRa9WyptqmSEwlgV10bqcb4HAFfMESN3dPrbCAQap8NdI0UslgsnHDCCTz77LMYDIkJzVw+OQejTsOs48/mmHO+wgfP/ZMtK9/tMm58jpUUo46tPeQUGHQavIEwjSpfZdhQ3uCMC17wB8M0tfkoTDcTDPj55x034Hba+fJP/sqS+TOHRUJgf8pQf0sIsQWYJoTY3OG1H9iceBGHFo1WyxU/uoeC8VN55t6f4m4beHconUZDMCRjxc8Ug4M3EKLG5sHSKZGsotmNlPRaPrrJGan9ntqhhaCUknA4TElJCa+88gplZWU9fv5IsRh0nDm7gCanj/Nv+DGl0+bw1N0/pK21OW6cRgimF6RS0eLG7e/ehCCRVKhqt8OGvY3xFUfrHV4kUJhhoqmmAkdLI5f/8G4yy6YNaaG5jvRnRvAf4FwiPYrP7fBaJKW8IoGyDRv0BiNf+cFvcLY289KDvzqsbZj0WtYdbBlkycY29Q4viHhbLMD+ZhdGnYbCHuoFhaUkHIYlE+JnAw8//DDnnnsuTmdyFPZRE7JINenwSQ1f/+l9XHDjT0nN7DpDmVaQipSR2vbdYdHr2Kb8BMMCpy+I3eOPmxHEKo6mmSgYN5kfP/o2s5edjACKumlePxT0p/qoXUp5QEr5FSnlwQ6vMXVXK506m1O++i2MlhTCoYGbeDItBrbVOFTp4EGkptXTJYZeSsnBZjfjsi09dnJqcfmZWZRGlvWQ2WfDhg3cfPPNhEIhLJauvWQTgVGn5dy5RbS4/WTmFXHUaZEMUpcj3p+Uk2Ik22pgV11bt9tJM+nYXedU7SuHAXV2b6SqgIjPKM40wBdvPE0oGMBothAMh9FqRJcqpENFf0xDn0b/bxNCOKKvtvb3iRdx+HD6ld/hou/c0Wt5gJ7QagThMMo8NIjsaXBi1sc7ehvafLj9oV79A75AmGOnHIoUCoVCXHvttWRlZfHvf/87qck888syKc4wY4vmCuzbspb/u/xE9mxYHTduakEqtXZvtzkFOq2GQEi1rxwO1NjccaXLpJTU2b04Vj3Ji3/+JQe2bwTA5Yuco7oENKI/HPozIzg2+n+qlDIt+kptf594EYcP7Vr+wI6NvPrPewb8eZNew1plHhoUpJSUNzq7RPy0h42Oy+7+qd7jD5Fm1sU1qXn44YdZt24d99577xEnjA0UrUZwwYJinL4gYSkpmTqblMxsXvjzLwkGDiWKtduSd9d3PyuQRGZCiqGlc2vKFpcfx8Ft7H//GZae+WUmzT0KAJcvyLTC4eEfgIEllF0ihEiN/v1TIcSLQogFiRNt+CClxOb2x6beezas5v1n/s7Gj98Y0HYyLAa21zgG1L5Q0T02dwCPP9QlkexAk5v8NGOPxeNa3X6Om5IbKwYWDoe5//77Oemkk7j00ksTLnd3TMixMq8knaY2HwajiQu//TMaKvfx0YuPxcakm/UUpJnY1YMisBp0bO2lHIUi8bSbJTuWnq5ssNH02r2k5RRw3vUdypQLKMtKjgmyPwxkXvIzKWWbEOJY4BQiNYceSoxYw4PID+viuXVVPL76IA99tJcXN1SRevRFFEyayfP339nFntsbWo0gLCV7eriYFf2ntpvuT25/kDqHt0ezULsiX9ihR4RGo2H16tU89thjQxbGJ4TgrDlFBKUkGAoz8+gVzF52Mu/8+6+0NtTExk0rSKXJ6e+2tHmqSc++RheB0JGVRFEcPg5vEHenh5NPnvkrwdYavvKDX2OyRGahUkqkjM8yHmoGogjavZxnAw9LKV8DRmX9WyklFS1unltXxUsba2jzBlk+OZvZRem4/SHWHLAjTrgRd5uDF//5wIC2bTHo+OKASi47UrpLJGs3jfTUhKbF5WdOSTrplkjIaGVlJYFAgIyMDEpLSxMqb1/kpho5cVoe9dGs1PNv/AlCI9i59tPYmCl5KQjodlbQ/pBRY1N+gqGiwRHJcu+IbuISJp19PVMXLIstc/tD5KUa40pQDDUDyZuvFkL8DTgN+J0Qwsgoa2wjpeRgi5sv9rdQY/eSYtRx4rRcZhalxRUuc/uDVLcW8Oym09nw1jNMO+lilsyb2a99ZJj17Kpvw+kLdrFvK/rP7vq2Lsdvf5MLi0HbYySGLxiKFZcLBoOcffbZlJWV8eqrryZc3v6wYlouq/c24QmEyMov5iePv0dKRlZsvdWoozTLwu56J8smZneZwUhgX5OLcX3UV1Ikhlq7J85R7PGH8OdM5Zijl8WNc/qCHF2cxXBiIDfyLwNvAadJKW1EKo/eNtAdCiG0QogNQohXo+8nCCHWCCHKhRDPCCGSPsuQUXPNU19U8r+NNTi8QVZMy+Xrx4xjbklGl+qVFoOOKfmpXPe9HzHhrOtZVQ+fH2jpVzngSEijZHcPoYCKvvEHw9S0erq0ATzY7GZCjrVbE4/LFyTTYmBidLbwwAMPsGXLFq655pqkyd0XFoOOc+YV0RQ1/bQrgcrdWwj4I8um5adi9wSod3RjHjJG2lcqhoa9ja6Yo3j/tvU8de/PCHudXUxAwbBkYs7gFDEcLAZqGjIBlwgh7gCuA5Yexj6/C3SsTfQ74F4p5WSgFUjalSmlZFddG//67CCvb60jEAxzyow8rjpmPPO6UQCdyc0v4Fvf+Q7Ti7JYvbeZ93c29Kvmi8Wg4/P9zX2OU3RPvcOLFMTlCdTaPfhD4R79AzZPgBXTctFoBLW1tfz85z/nzDPP5Pzzz0+S1P1jUVkmBWmmWJhofcVe/nTTJbz/7N8BmJRnjTWu70yKSddrBrIicbT7Ey1GLVJKXn74d5Sv/RCNVkdeWvwMVUCXctRDzUAUwf+A84j0IXB1ePUbIUQJER/DP6LvBXAS8Hx0yOPA+QPZ5uHiC4Z4a1s9b26rQ6sRnDm7gCuXjWNWUXpce7m+0Gk0FDRvIPDW3WyptvPallrCfcwM0s16yhtdtHlVv9nDofMUHCJmIa0Q3UZihMISjYB5pRkA3HXXXXg8noQUlDtSdFoN58wtwu4NIKUkv2wSs5efyofP/ROnrQWjTsv4bAu7G9q6nGeaaHtU1Qgp+bT5gji9IQxaDZs+eZODOzZScsrV5Geno++QKxAMhdFpBdkpwyORrJ2BKIISKeVlUsrfSyn/0P4a4P7+BPw/oD20IRuwSSnbH2GqgOLuPiiEuE4IsVYIsbaxsXGAu42n3uHlqc8r2V3fxrKJ2XxlSRlT81O7lCroLz6Pi5qNHzHesZl9TS42VNh6Ha8RAiQqeugw2V3vxNQpbHR/k4viTHOXcFKAZpeP+aUZpJr0BINBPvvsM6655hqmTJmSLJEHxPSCVIozzLHG9Gdd9T38Pm9sVjAtPxW3P0RVa1fHsEYj2KnKTSSdBocPoYFQIMBr//wDBeOnEp58fJf+xC5/iHFZ1gE9bCaDgSiCVUKIOYe7IyHEOUCDlHLd4XxeSvmwlHKxlHJxbm7uYckQDkt21Dl4YX0VobDkokUlLJmQddgKoJ2jTr2AwglT2fa/h5iQqWfV3qZYTfKesBi1fK6ihwaMlJK9Dc648tI2t59Wd4DxPSSRBUKSJRMi9nadTsfnn3/OPfcMPCEwWWg0gnPmFtHmC8ZmBYtOOpdPX34SR3MDE3KsGLSabs1DGWY9m6vsqn1lkqmze5ESVr7yH5prKzn+ilsIoeliAnL5gkzNH17+ARiYIjgWWC+E2BWtPrpFCDGQ6qPLgfOEEAeAp4mYhO4DMoQQ7Vd1CVA9gG32Gykl33lqPRsr7YzPtnL50WUUD1LBJ41Wy3nX/pDm2kos5e9j1mt5c1sdwV5iutPNesobnMo8NEAcniBtviCGDtPtA72EjfqCISwGLeOzrdTV1eFwONBqtYPWcSxRTM1PYXyWFZs7cn6cdsV30BtNVO7Zhk6rYWKulb2NToKdWqma9Foc3iBNTtW+Mpnsa3Ji1muZs/wUzr7m+5gmLAS69wWUDqNEsnYGogjOBCYTCR89l0hXsnP7+2Ep5e1SyhIp5XjgMuB9KeXlwAfAxdFhXyfiixh0hBCcOC2PxeMyOXN2ASb94MbwTlt8LNMWH8sHzzzECRNSaHH5WVnes0NYmYcOj1qHp0tHsv1NLjItejIsXQPOWlx+lk7MRqfV8P3vf5/Zs2cTCAx/5SuE4Ky5hTijs4KcojJ+/p+PmbX0RACm5qfiC4ap6K6shISDqix1UjnQ5MJq0JJVUMLJl15Hrd1DilHXtcy5lBQMM0cxDEwRVADHAV+PdiqTQP4gyPBD4FYhRDkRn8E/B2Gb3XLJ4tJIUk6CHITnfvM2vvqD3zK5OJf5pRlsrLL1ekFajFqVXDZAKprdcYlk/mCY6lZPt7OB9gtvXmkGW7du5amnnuLyyy9Hr9d3GTscmZRrZXJ+5KECwGA0IaWkeu8OyrIsmHQadtd3LWJoNmjZrMJIk4bTF6S+sZGnfnsL9RV7gUjme+fZgC8YJt2sj1MOw4WBKIK/AsuAr0TftwF/OZydSik/lFKeE/17n5RyiZRyspTyEinliG21VDRxOrOWnYQQguWTssm2Gnh7e32PpafTzXr2KPPQgNhV1xbXaL6y1U1Iym7DRl3+ELkpJorSTfz85z8nNTWV224bcOrLkCGE4KzZhXgCoZjN/4Pn/sm9376IltoKJuelsK/J2aWsRJpZx+76NlVuIkk0OLysf/0/bP7kLWQ4jNMbpM0b7KIInL4gU/KGp0lyIIrgaCnltwEvgJSylVFaYuJICIdCvP7ovXz+xrOcPqsAXyDMh7sauh0bcVJL9nTzVKfoSiAUprJTItn+JhcGrabbBh92T4BjJmezfv16XnzxRW699VaysoZXRmdfjMu2MK0gLWbzX3zKl9Dq9bz97z8zNT+VQEhyoFM/Y51GQ0hKqruJKlIMPuWV9ax//SnmHns6BeOnxMqBd44Y8gVDTB4FiiAghNASMQkhhMjlUBioIopGq+XA9g289e8/k2GExeMz2d3gpM7RtUgatNceUqWp+0NDmw+JjIXeSRm5CY7LtnQJx2t/gp5dlM7LL79MVlYW3/ve95It8hEjRCTHxRuMzArSsnI59ktXsP79V9Daq7AYtN1XJJVQ3qj8T8ng4b89iN/j5JSv3gBEzEJajSA3tXMimaBwmHQk68xAFMH9wH+BPCHEr4BPgd8kRKoRzqmX30hbSyNr3nyehWWZmPVaVpU3dTs23RQxDzlVaeo+qbXFdyRraPPh8oe69Q/YPAEm5ljJtBr4xS9+wdatW0lPT0+muINGaZaFafmpMV/BSV/+JgazhfeeeoipeakcaHbjC8abH9NMejZVKj9BonE6nbz+1D+YvmQFJZMj9cZq7V7y04xxDyfhsEQAeanDz1EMA1AEUsoniSSD/QaoBc6XUj6bKMFGMpPnHc2EWQt5/5m/o5FBjhqfSWWrp9uMT1V7qP/saXBi1B2K9mo3iXTXhMblC7JsUjatrRFnfGFhYXKETBAnTc/DE4jc7K1pmRxz9mUc2LGRCZk6QmHJvsZ485DVqKXO7lH+pwTj8gWYd9qlnH7Ft4FI5nBDm7ebRLIgRT0kPA4HBtKY5ndSyp1Syr9IKf8spdwhhPhdIoUbqQghOOWr38LWWMvad//HnJJ0Uk06VpY3dZvoYzHoVOeyPmgvDJjawVG8v9lFQZqpSxOaYDiMTqtBY6+lsLCQF198MdniDjqTclPISTHGbuynXn4jtz/yJqW5GaSadF06l0Ui44TqWpZgPNLAcZfdyLjpc4HILDUsu+YPOH1BpuYPn45knRmIejq1m2VnDpYgo43pi4/jxEuuoXTqbHQaDUsnZtPQ5qO8satjON2kZ0+9Mg/1RpPTT5s3iDGa/+HyBal3+JiQ29Us1OoKMLc4nT/ff28kgmv58mSLO+hoNIJTZ+Zj90TOEZMlBZ3eQCgYYGK6looWd2zG0I5OK9ihyk0kjGeffZZnnnuBUIforPaGSZ0VgZQ9t08dDvSnef23hBBbgGnRjOL2135gIJnFYwohBOde+/8onjQDiNSPybIaWL23uUuFUo1GIJV5qFcqml10PGr7o2ahid34B/yhEOMsfp544gmuvvpq8vMHI91l6JlTko7ZoIn5A3weN7+5+nQaP32GsIS9DfEPGelmPVur7f2qiKsYGD6fj1tuuYXH//5XTB1mpLV2D+lmfbetUgvShqd/APo3I/gPkQzil6P/t78WSSmvSKBso4KGqv28/ui9ICXHTMqm1R1ge13XpzQVPdQ7W2vsmDtkg+9rcpFm0pFtjY9g9gVDmPRaXvrX3wkGg/zgBz9ItqgJw6jTcsLU3FgoqdFsoWTqbDa+/Syp2kCX6CGjTovLH6Khj7pXioHzxBNPUFNTw9KLro2FM0spqbF5Keo0G/AHwxj1GrKswzfavk9FIKW0SykPSCm/IqU82OGl7lr9oHLXFt596iF2fvExE3OsFKSZWLOvpUsdovbaQw7l3OtCMBRmV52TtGihuUAoTEWLm4k5XbPEW1wBFhRZeeSf/+CSSy5h4sSJQyFywlgyIRsBsRpDJ335WjxOB5pd71HV6sHVjXlxXzfmSMXhEw6H+eMf/8j8BQvJnLIQkz5yG7V7AngCoS4lJFy+IJO6OVeHEwNxFhuFEF8VQvxYCHFH+yuRwo0G5p9wJhk5BXz4/CMIIThmUjZOX5DN1fGhfRohkFKqEsLdUGv3EghFHMAAlS1uQmHZrX8gLMMsmVzAunXr+PWvf51sURNOulnPUeOzYrOCcdPnMnne0ez/4FlkMNDFaZxi1LGx0jYEko5eXn/9dXbu3MnV138HjUYTu8HXxfwD8RFDbn+IKcPYUQwDb0zzJY6gMc1YRKvTc9wFX6N80xqq9myjNMtCSaaZ9RWthDrZblNMelbvVZ3LOrOvydnpfSSbuHP1WJcvSKbZQEmmmQkTJoy62UA7x07JIRgKxyLQTrr0WtpaGtBVretSeyjVqGN/k6vHMieKgRMOh1mxYgWLTzwT2WFiX23zYNBqyE6JNwEJAUUZw9c/AANvTHPpETamGZMsPevLGC1WPnz+EQAWj8vE5Qt1qSefZtJR2eqh2alsuh3ZXGWPNaqXUrK/ycX4brKJ7Z4Ato1vcPbZZ+NwjN6ZVVGGmcn5qbRGS1RPW3QsN979BEeffBZ1Dm+soQ0caud5QFUjHTTOO+88PvjgAw62+uLyAqpaPRRnmuP6m0gpkXSdJQw3ktaYZixjtqZy3JeuxJqRhZSSsiwLOSkG1lW0xuUVtE8xtyvzUAyPP0RFizvWiKbe4cPtD3UxC0kpCYXCvPjYQzQ2NpKaOryn4kfKydPzYr2JhRBMnnc00wrSALo8YGg1gu016pwaDN599108nkiG+56GQ34rpzeIzROgpNMs1ekLUpxhxmwY3LL3g00yG9OMac66+hYu+NZPEEIghGBRWSYtLn+sqUo76eaIeUh1mIpQ1epGSmJPWfuanAhBl2qjDm8Q55417Ntbzm233TasHXODweTcFPJSjXGZw2tffhzHy79iV11b3PmTYdGzqcqmwkiPkIqKCs444wz+7//+D7sngMMbjM0IqmyR67gkM14RtHmDzCpKS7qsA2UgiuAMIo1pTiXSlObs6P+KfiKlZO/mL/C6nEzJTyXFqGPdwfh+BFaDlnqHl3qHMg8B7K5vo0MzMvY1uShKN3dpLOT0BtnwxpOUlpZy4YUXJlnK5BNJMCvoZAbS0LpjNXX7d8R1KDPqtHj8oR4LHyr6x3333QfADTfcQI3NC8jYA0dVqweDTkNOatem9BNzh2fF0Y70J6GsTQjhALYCW6L/bwW2Rf8fs0gpB/TkXndgD3/5wRV89sZzaDWCBWUZVNs8sWgDiEzzNUKwtcaWAIlHHluq7aRFG3k4PAGanX4mdjILhcKSxoM7Wbv6U2666SZ0uq7JPKOR2cXppBh1sYzio8+8BIPZQtsXL3VTkVR2SThT9B+73c7f//53vvzlL1NWVsaBZmecL6Cq1UNxRrx/IBy9N3SeJQxH+pNHkCqlTOvmlSqlHP5zngTh9gepbHVTY/dSY/NQa4+8qlvdPRb6KpwwlUlzl/Dxfx8nFAwwuygdg07Duor4WUGGRd9tBvJYw+b20+T0xxLJ9kWziTtXG211+1kyZyq//vWv+eY3v5l0OYcKg07DyTPyY1VJzdZUlpx2Ea6dn7Ct/EDcQ0qKUc/6TueZov/84x//oK2tje9///sA7KhtI9UYeUBp8wawewJdbvhOX5DSLMugt8VNBMOzFN4wp9nlo80b5BvLJ/DbC+fw03NmcsupU7n+hEl8ZUkZeq2GylZ3zJnXkRUXXY2tsZZNn7yFQadhbnE6exuc2NyHpvIWgw67O0C1bWw3FolUaz00/d7X5CTToiezU29iXyDESfMmcfvtt5OZmTkEkg4di8ZlotOIWDey4y/4GjIcom71y9R0mGmmRiPSVDXSw2P16tWsWLGCRYsW4Q2EqLN7sRgjN/j2BkBdFMEI8Q+AUgQDIiwl1TY36WY9t546lTklGei0GtLNegrTzUzKTWHR+CxuO30aX11SRiAUpqrVja9DMbAZR68gr2QCH734GFJK5pdmoBGC9RW2uH1pNBo2V8UvG2tsq3Fg0EYuNl8wRHWrp4u9NRAKs+Pjl9nw4etDIeKQYzXqOH5KLk3RkOOcojLOufaHpE5fHhc9JIRAgKpGepg899xz/O9//wMiiWOCQwEMVTYPRp2G3JROjWhE19nrcEUpgn7iD0Zu6ovKsrjppCnk9VJASqfVsHh8Fj86cwbnLyjG5gngiDr1NBoNx51/JS11VThaGrAadUwvTGV7rSNuBpFl1bNmf0uXpLOxQjgs2V5jJ90cmX4fbHYTll2LzNW3tvHpk/fx9FP/GQoxhwXLJmVHwmej58pJF1/NjNlzKW9wxp0/eq2GLdWqWc1AsdvtCCFIS4s83Ve2RM7Fdtr9Ax0j1drNusUjwD8AShH0C38wTJ3Dw0ULS7hsSWm/bX4mvZbjpuTy7RMn4w6EYnVglpxxMT/794ekZ0eqYi4syyQUlmyuOnSRGnVa3P7gmE0Eqm/z4gmEY+F55Q1OzHptlzouWz96FYethVtuuWUoxBwWZFoNLBqXFZsVAOQGG6l59xEqOpw/GZZINdKx+nBxOKxZs4bCwkLef//92LJd9W0xs1Bv/oFxWda4RkrDGaUI+iAYClNn93LRwhKOnZJ7WPHppVkWrlk+nlZ3AG8ghN5gxGA0EQ6F8HncZFkNTMyxsqnKFrP1Ahi0mjFbcuJg0yEThi8YYl+Ti6n5KXFRGS5fgI1v/Id58+axYsWKIZBy+HD81FwCHcpOhBr34vjsOVZ9/EFsjF6rwR8MUzPGfU8D4b777kOv13PUUUcBkQi1vY3OWIOkqph/IL7XQJsvyKzikeEfAKUIeiUYDlNj93L2vAKWT845om1NLUjjiqVlNLT58AfDBPw+fvfNs3jnPw8CsHBcJt5AOC6rONtqZGNFa9yT3lhhfUUr1mg25p6oiWN6QfyFtXHNJzRUlPO9731v1CeQ9UVRhpmpBYf6Gi9acTbGtCx2vPt03MMFwJ7umt0rulBdXc1zzz3HN77xjVimepPTRzAsYwUQq1oj/oGcbuoLdU56HM4oRdADoXCktvgpM/I5eXr+oNxoFpRlctHCYuocHoROT8GEqax+/Rn8Xg9F6SYK0kxsqLDF4o81GoFGI1jZQ+P70YrdE2B/k5O0qH9gZ20bmRY9+WmHnHGhsESGQiw/7ni+8pWvDJWow4qTO/Q11hkMLDrjUtx717Ju46F0nzSTCiPtLw8++CChUIibbroptqzG5kF2MK1VtbopyYz3D4SijeqLMkaGfwCUIuiWSHSQh+Om5HDm7IJBfdpcPjmH02YVUGPzctz5V+Jps7PuvZcRQrBwXAZ2TyAu8Sc3xciq8uYxFfa3u64tVlbC7omE0U4vTIv7HVpcfs4752w+/fgjjMau2ZxjkYk5KWSnGGO+qNMvuhKh1fPJ//4VGxNpau/D7h4759PhEAgE+Nvf/sZ5550XV8V2T70TQ9Tu74iWmehiFvIGmJBjHbaN6rtj5EiaRGptXo4an8WX5hfHqjcOFkIITp9ZwKyiNFLGzaF48kw+eelfSCmZlJtCulkfV4xOp9UQlpLP94+dPkBr9reQEs0m3hnt5ja9IL6I3L4tXzC/UCmAjmg0gpOm5WHzRMxDqZnZTD7+AtzCFOuHLYQAIbv0LVDEo9fr+fjjj+N6WkgZOW6p0UJz7Xk+Xcuhh5hdlJ48YQcBpQg60eT0UZRh4sKFxV3KHA8WGo3gooUlaDQalp57BXUH91C+aQ0aIVhYlkG9wxetZRIhO8XAB7sa8AZGf015uzvAgWgbykijnjZKMsyxMhMAzS0t/Pe33+HPv1Z9kToztzQdnUYT64B32Xd/RsbxX4spVIhkGa/eN7bMjYfDjBkzmDlzZuy9wxvE5g5gjD7pV7a6Mem79w+MG0H+AVCKIA6XL4iU8LVjxic8LTzTauDiRSUULjyZq+/8K5PmRKISZhamYdZr48pOGHVavP4Qm8ZAp6mddQ4QkSfXOocXmyfA9ML42cCnrz+P3+vhW9/61hBJOXyxGHQcPeFQB7NMi4GCNCOfr1pJMBgxB6WZdFQ0e8ZkEEJ/+PDDD7nkkkuora2NW15r9yA0h8rFV3eTPxAMh9FqBIXDvBFNZ5QiiBIMhWl2+fnasnHkpCTH5LCwLIP543PJm70cTTSDVqfVMLcknf1NrlgECECGxcDb2+u69DoebazZ3xJrQrOjtg2dRjA571A2cTAQZMMbT7HsmOUsXLhwqMQc1iydlE0wfCiUNL15O+WP/oBV778FtJuHYFuNSi7rjnvvvZcPP/ywS7mSg81u2m/5PfsHgkzMTUGvHVm31pElbYKQUlJj93LWnAKmFyYv9lcIwUULi9EIwetP/IU3n3gAgHklGeg0Ii66w2rUYXMH2DGKm9a0uvxUNEfMQsFwmN31bUzMjU/K+fyTd3E01nDrLd8bOkGHOQVpJibkWLFFs9lPWHESuvR8PnnpidiYDLOeVeWq70Vn9u7dyyuvvML111+PyRT/VL+j1hF7SInUwepaX8jjD42Y+kIdUYoAqHN4mV2czknT85O+7wyLgUsWlVBdsY8Pn38Ej9OB2aBlZmEaO2vbYhEgAKkmPW9vrx+1F+/OujaINu450OTGFwwzo5Ni3rvpM4qKSzn//POHRsgRgBCCFdPyYueO2WRgwvEX0ly+iYo92wGwGLQ0u3xxhekU8MADD6DT6bjxxhvjlvuCIWpsHqyGiCLY1+Qi1aQj22roso3xI6S+UEfGvCJodftJt+i57KjShDmH+2J+WQZf/vp1+L1u1rz5AgALyjIII+N8BWkmHdWtHvY2js668p/ta45lbO6sc2AxaCnrMPX2BEKc/62fsH79ujHTc+BwmVaQitWoiwUYnHjeZQi9iTefifTNbu97MRb8Tv3F4XDwyCOPcOmll1JUVBS3rt7ui4Q0Ryu9VrS4mZSTEucf8ARCpJp0FPZSh2y4MqYVgccfwhcM843lE7Aah+7GIoTge5edTsmMhXz80hOEQkEyLAamF6Sypcoee7ITQpBi0vHqptpR16ug2emjqtVDqkmHxx9if5OL6QWpceG79a0OjpuSS35e7hBKOjLQazWcMDWX5qifaWppPhnzTuHAls8JBiLLsiwG1uxrVrWHooTDYW6++eZu61YdaD708FXR4iYUll0aJLW6/CyZkDXoIefJYMwqgmAoTKPTx+VLyihMH/oMwAyLgW9cfyO2hhq2rY4UuFoyPouQjJ8VZJj1HGxxjzpfwY5aB4KIsttZ5yAsiSspYW9p5m/Xn8quT14dOiFHGIvHZSGIZLpqNIJll95IwTUPESDiczHqtbjGcGHDzmRkZHDXXXd1G4SwsdIe8w/sbXRi1Gm6ZA5LJLNGWP5AO2NSEUgktXYvp8/KZ25pxlCLE+O2665kwWkXY8kqBOhxVpBh0fPypppRE0EkpWTN/hbSzDrCUrKpyk5huoncDv1f33/pP/hcbRy/fOkQSjqySLfomVuSTnM0THTexBKkVs/O6lZCocj5pNNoWH9QlZz49NNPeeWVVwiHu15TTl+QihY3qabI+Xmgyc34bGucKdkbCJFm0ndJLhspjElFUG/3MqsojVNnFgy1KHGkmA08+NBDmAonxZZ1NytIM+lpcvpGjX232eWn1uYhxahjX6MLuyfAgg4K2u/1sf7NpznhxJPjEnwUfXPslFx80aqkualG0kN2nrntAta/H5lZZVkNrDvYii84+pMVe+NnP/sZ3/nOd7pVBAeiLVKFENTavHgCoS5moRaXn6PGj0yzECRREQghSoUQHwghtgshtgkhvhtdniWEeEcIsSf6f2J7DQrISjFy2ZKyIXMO98aiskyw1/Lxa88D3c8KIJIo9MrmmlFxAW+psseihTZUtJJm0jGpQyeyT995GZetidt/eNsQSjkyGZ9tIT/1UImJudMmIvUm3nvuEaSU6LUagmFJ+RhubL9x40Y+/PBDbrrppm6DELbXOtBrD7VL1QgYlx2fPyCB2cUj0ywEyZ0RBIHvSylnAkuBbwshZgI/At6TUk4B3ou+TxjT8lOH3DncGzqthvrP/sfLf/459uYGoPtZgdWoo80b5IsRXoPIGwjx7o56cqwG6hxeauzeSPvOqJKWUvL5q08yZdoMTjvttCGWduQhhODE6bnYozkFM4vSyVhyPg0HdlG+8TMATDrNiD+PjoT77rsPq9XKNddc02VdOCzZUmUjw2xASsm+RhelmZa43JaIWUg3Ys1CkERFIKWslVKuj/7dBuwAioEvAY9Hhz0OnJ9IOa5aPqFLl6vhxh0/vBUZDvHef/8N9DwryEkx8sbWurgWlyONtQda8AVDGPVaNlS0YtBqmNkhIafF5efbdz3Ao//8+5jvOXC4zCnOQK/VEAiFMem1LDrpXDSWdN5//lEgcn5tq3GMqQq37TQ0NPCf//yHr3/9610yiQFqHV680U55re4ANk+ACZ2jhdyBEW0WgiHyEQghxgMLgDVAvpSyvahHHdBtVpcQ4johxFohxNrGxsbkCDpETJ06lZNPP5N1bz6LzxtJ+OluVmDSa/EFw3y6Z2QWEPMFQ7y9vZ4sq5E2b4A9DU5mFafFnraklLgDIS4/9WiWL18+xNKOXMwGLcsmZcfqDy2YkE/qgrPZ9cVHNFTtR6sRSGBDhW1I5RwK9u7dS3FxMTfffHP36xsOVWltz9/p3Dc7LOWINgvBECgCIUQK8ALwPSllXAykjKTMdhvULKV8WEq5WEq5ODd39MeR337b9/E4Wvnw9ReBQ7OCzVV2nN5DM4C8VCPv7WjA5vb3tKlhy8YKGy5fELNey6ZKO0iYX5IRW7+3fDfv3PcDcNQPnZCjhKMnZBMKH3IaTznhAsZ/5RdkF5YBkZyCD3Y2jJpItP6ybNkyysvLmTZtWrfrO4aN7mt0kZdqJLVDJVxvNIlsJJuFIMmKQAihJ6IEnpRSvhhdXC+EKIyuLwQakinTcOXEE0/kqCVLcdmaYhfn0ROyAVi599AMQK/VIAS8sL56RJWeCITCvLm1jiyrAX8wzJYaO5PzUmJdyQA+fekJdq/7NNYmUHH45KcZmZiTEqs/tHjGBGTZIipbIzNOs0GLwxtgd/3YcRrv2rULr9eLRtP9bbBj2KjLF6TO4e2aRDYKzEKQ3KghAfwT2CGl/GOHVS8DX4/+/XXgf8mSaTgjhOCz1Sv50e0/pr4tEgeebtazsCyDnXVt1NoPNSDPSzWytdrGxhE0td9cacPhDWAx6Nhe68AfDLOw7JCNtq6hgR0fv8oVl19OXl7eEEo6OhBCcMK03Jg/aVJuCma94KV//IGP/xspRmc16vhgZ8OIeqA4XEKhEOeddx4XXnhhj2MiYaMSIQT7oyGkE3NS4sZIKZkzws1CkNwZwXLgSuAkIcTG6Oss4LfAqUKIPcAp0fcKQKPRcMK0PGxV5fiiNWMWj8vCatDy0e7G2AUrhCA3xcRz66pGRAvCYCjMG1vryLQYCEvJxkobhemmOCf+yleeJuDzcuuttw6hpKOLaQWpWAyR+kNajWBuSSYN+3fy7tN/I+j3k2HWs6/JSe0YKET38ssvs3v3bq666qoex0TCRiO3yPYicx2b0HgDIVKMI98sBMmNGvpUSimklHOllPOjr9ellM1SypOllFOklKdIKcduHFs3vPLic/zz1ovZtP4LAAw6Dcsn51Dv8LGj7pAjy2zQEpaSFzdUDfsnum01dlrdfqxGHbvq2rokkDXanWx+51lOO+10Zs2aNXSCjjL0Wg3HT82J9bmYXZxO+pIv4WxtYsOHryGEQK/VsGrvyAw+GAh33303EyZM6HFG0DFstKcic63uAEeN0NpCnRmTmcUjiXPPPZfMzEw2v/6vWFLQ9IJUCtJMrCxvwh885NzLSzWyucrOxmGccRwKS97YWke6WY8vGOLT8iYK0kyx5jNhKWlze7n6mm9y++0JTSkZkywelwVEbnQpRh2zjzoOQ+44PnzhUaSUZFsNfLG/ZVSHkq5cuZLVq1dz66239ljFtmPY6J56J6GwjGuQJKUkFA4zfxiVqDkSlCIY5qSkpPCtb32LHZ+9x77yPUgZsVmeMDUXtz/EFwcOTaAiJiIjzw9jE9GmShsNbT5STXo+29eC2x9ixbTc2JNWY5uPpVOLue/3v2bFihVDK+woJMNiYG5JOk2uiN9pfmkmqUddSO3+Xez44mN0Wg0hObpDSV988UWysrK4+uqrexzTMWx0S7WdTIueog7tJ22eABNzUigc5jlJ/UUpghHATTfdhF6vZ/d7T2OL3uAL0k3MKEhlQ4UtLnTUbNASDEv+u6Fq2JWqbnH5eW5dJTlWY6RWUpWN2cVp5Efrt/uD4UiZ5H1ruq35ohgcTpmZTyAYJhyWFGWYKFtyKjmLzyIjN1J7a7SHkt5zzz2sW7cOq7XnBjLtYaONbT7qHF7mFKfHmYVcviAnTs8bNUmOShGMAAoKCvja177GvnUfY3d7Yj6AYybnoNHAJ3ua4vwC+VET0Zvb6oaNvyAUljzzRQVIgUmv4cNdjRi1Go6ZlBMb09DmYcPzD3DXz38ybOQejRSmm1k4LpNGpw8hBIsn5GI9+UZc1kgzlvZQ0l31bX1saeTh80W+8/jx43sc4+oQNrql2o5WI+I65bn9QdLNeqbmp/S4jZGGUgQjhF//+tfs2b2LBeNzaIiGk6YYdRw9IZt9TS621hzKzRNCUJhh4p1tdawsHx6Ov1V7m9hd30ZuqoHd9U6qbR6OmZSDWR/JInb6gjTu/ILd2zbzwx/+EK1W28cWFUfCKTPzCYbChMKSaQWpZFr0fPj5Jt5/7p/A6Awlra+vp7i4mOeee67Xcdtq7IAkEJLsrHMwNT8Fk/7Q+djqDnDi9Dx0I6xBfW+Mnm8yysnNzSUlJYVz5xQgg8FYPPjCsgzGZVn4aHcjDW2Hwv50Gg0F6WZeWF/NhoqhrTdfa/fw8sYaCtLMBEKST8obyUs1Mqs48pQlpaTV5Wfnm/+iqKiIK6+8ckjlHQvkpZo4emI2jW0+NEKwbGI2tZs+5tW//57K3VsjoaSNrlHVFvX++++npaWFefPm9TgmHJa8t6OBdLOBXXVtBELxeQLBUBiNECwoS2yR5GSjFMEIwmazsfyoBfg3vUKT00846jg+fVYBZr2W17fUxZWlNug05KUa+fdnFeyqG5qOZv5gmCfXVGDUaTDoNKzZ34zLF3EQa6L21SanH9Gwmy9Wf8oPfvADjEZjH1tVDAYnz8hHIgmGwkzOS2H8seejMVp57+mHEUKQatLx0oaaUdHKsqWlhT//+c9ceOGFTJ06tcdxexudNDl9WA1atlTbyUkxUNChB3GT08/RE7NiZSdGC0oRjCAyMjIYN24c//77n5mVa6Tecag8wJmzC3B4A7y7I346b9JrybToeWTlASpb3EmX+d0d9dTaPGSnGKlu9bCx0sbMwrRYe1BfMEQwHOaoYjNLlizh2muvTbqMY5Usq4Hlk3NivoJjZ5aRsvAcNq98m/qKvWRYDNTYPWwexuHI/eVPf/oTDoeDO+64o9dxH+1uxKzXUu/w0ej0xTmJpZQEw2GOmZSdDJGTilIEI4yf//znNDU1UbfmZcx6bawAXVGGmeWTcihvcLKpyh73GatRh0Wv5S8flLO12pYUOaWUrN7bxDvb6ylIM2H3BHhtSy1pZj3HTcmJjal3+LhwQTGXXXgua9asISVl9DjgRgIrpkbKdwRCYSbkWJl04sUIrYH3nvk7EIkgenlTDd7AyG2A5Ha7uf/++7n44ouZO3duj+MaHF521jrItBrYUm1HrxVxfbNt7gCTclOGRY/zwUYpghHGMcccw2mnncZ9f7yHC+bk0OL2x6buC8symJBj5ZM9jdR1KhOQZtaTatLxz08P8PLGagIJDA0MhsK8tLGGZ9dWkZ9mJCQlr2yqISwl580rijne6h1e5pWkY9uzFo/H08dWFYkg3aJnxbQ8Gtsis4Lj50widfG5uIUZKSVWow6HN8Bn+5qHWtTDxmKx8Mknn/Cb3/ym13Gr9zWj0Qj8wTC769uYlp+KQXfoFunyR0JGRyNKEYxA2mcFq9/6L8dNyaEuaiISQnDazHysRh0vb6qhKdq0vB2LIVIX5aPdjTz44d5YqYHBxOUL8sjK/Xy6p5HiDDN6rYY3t9bR4vZz1pxCMi2RWi1t3gBmg4556R7OPvvsPi9SReI4bkoOOo3AFwxRlmVh1nk3wNFXxh4wclONvLW1bkRmG7ebSefMmcPkyZN7HOf2B1m9t5mcFCM769oIhiVzStLj1kdCRkdnJVylCEYgxxxzDO+88w7XXXcdZ80pJMtqiN3UTXotF8wvRqsRvLCuKuZHaEerERRnmKmze/nD27vYWm0btMSzBoeX+9/fw94GJ8UZZrQawcryJg40uzlhai5lWZE+r8FQGJsnwBVLx/H7X/0fJpOJG2+8cVBkUAycVJOes+cW0hCdFSyblI3LH+Lt996npb4ao05LSEre3znyKsTfeeedXH755YRCvZu2NlS0EgyF0WkEW6rs5KcZyUs95CRucfk5dWb+sOxzPhgoRTBCOeWUU9BqtRh1Gq45dgLBsIzVIsq0Grh4UQkGnYYX11dTY4s3uwghyE01YtZreeTTA/z5g3L2NToPO2a82enjtc013PP2Lty+IAXpZoQQbKuxs77CxtzidOZ1aDhT5/Byyox8HFW7efrpp7nlllsoKCg47GOhOHKWTsymKMNMq9tPSaaFIlOA9/50Cy//M1IxPjfFyCe7m+JClIc7LS0t3HvvvQQCgV7zUkLhiJLLsBjYXuugxe1nQemh8FCb209+momjxmclQ+whQSmCEczzzz/PrFmzsGqCXHPsBGxufyx8NN2s5+JFJViMWl7aWN1txJDVqKMk00xjm48/v1/OPz7ZT1Vr/yKLpJQcaHLxxOoD/Ob1HXy0q5GcFCNZVmOk4fz+Ft7d0UBpppnjpx7qKFfv8FKWZeHUmfncfvvtZGdnc9tttw3OAVEcNjqthi8vKsXpDRIKS05fOIX0xeey+aPXqN67E51Wg04reG1z7bArXdITf/zjH2lra+szUmhPQxs2dwCdRrCyvJnCdFMsazgsJW3eIBcuLBlVCWSdGb3fbAxQXFzMjh07+Otf/8rkvBQuW1JGncNLMFqnJ9Wk5+KFJaSZ9PxvUw2769u6PPULIciyGijJNLO/ycUf397N3W/t5D9rDvLx7ga21zios3upbHGzpcrO+zvq+ddnB/jtGzt54P097Kh1UJhupjDqD/AFQryyuZbV+5qZVpDKufOKYtPpeoeXrBQDVy2fgNftoq2tjR//+Mekp4/8xh6jgbJsC8dOyaG+zUuaWc8ZX70BjcHCsw/9Hmivbmvj02GSrd4bzc3N3H///VxyySXMnj2717Hv72jAYtCxZn8LnkCIE6bGF0FcUJbJpNye6xKNBkZXVsQYY9myZZx++un8/ve/59prr2XxuEwa23y8s72O0kwLQgisRh0XLSrh5Y01vLG1jl05VlZMy43ruwqHzEVhKfH4Q2yvccQqUAoBUh5qJm3WazHpNRRnmOOKbjW2+XhtSy1t3gArpuYyt+RQDHad3UtumpHrjp9ImkkP6Fm5cqUqLjfMOH1WARsrbbj9QY6aVsqqE79CxZv/YOu6z5i9aCkFaWb+t7GakkwzE3OHb6hvv2cD9W3saXBiNWgjRRCL4osghiWcPbdw1BSX6wk1Ixjh/Pa3v6WlpYVf/vKXCCE4Y1YBi8qyqLYdKk5n1mu5ZFEJx03OoaLFzb8+O8imShvhbnwCmqjyyE4xUpRhpijDTGF65P/i6CvLasBi0MUujmAozKZKG8+urSQYDnPRwhLmlWZ0UAIe8tONXB9VAuvXr6e+vh4hhKopNMywGnVcuKCYJmck+ODyq6/DkDeBz7btQ0qJQachzaznsVUHhm2pc4Dvfve7PProo73OBtz+IE99XkGGWcfHe5rQaTUs65As1uj0cvqsfLKshh63MVoQI7Go1OLFi+XatWuHWoxhw/XXX88jjzzCgQMHKC4uxhcM8cSqg2yvdVCUboqzbdo9Ad7f2UBFi5uCNBNHTchkXJb1sKIhnL4gm6tsbKm24w2EKck0c8asAqzR9HspJXUOL8UZZq45biIpRh3BYJBZs2aRkZHBmjVrBu0YKAYPKSUPf7yPg80uclNNrD/YwiflzZw2Mz9WhbPd13Pt8RNj7RyHC+09O/ri+XWVfLavBW8gxKubazlham6s0YzTG0Qi+eGZ0zHqRs/DihBinZRyceflw+sXVBwWd911F++88w7FxcUAGHVavnHsBM6YXUCN3YsrGk0EESfy+fOLOH1mPnZPgFc21fL3T/bxzvZ6Dja7enUEhqWk1e2nvMHJW9vqeHTlfr440EpxhpmLFhZz4YLimBIIhsNU2zyUZVn4ZlQJADz22GPs3r2bH//4xwk8IoojQQjBBQuLCYYlvmCIBWWZFKTqefWlF3C4I7kpealG9jY6eXNr3RBLG88777zDsmXLqKys7HXc7vo2VpY3k2M18PHuRrKthlhxORk9zy9aWDKqlEBvqBnBKMPv92MwHJrK7qi186/VFSAiIYAdCYUlFS1u9tS3sbfRhT8UxqjTkGLUYdBpMOo0kQtBQKvLT7PrUBazXiuYVZjOvNJ0MizxU2eXL0izy8fJM/I5fVZB7InR5XIxbdo0SktLWbVq1ai3u4501uxr5unPKynJMvPFx+/yzK++w8QLf8B137wGg05DKCyptnn42rJxw6IaZzAYZN68efh8PrZt29Zj8UK3P8jdb+5CCNhR28bqfc1csKA4ludSZ/cwOS+Vbx43YdSdoz3NCJSzeBTxhz/8gUceeYQNGzbElMGMwnRuPW0q/1p9kKpWN/lpptiNWasRTMixMiHHSjAU5mCLmwNNLjyBEL5gGLc/RKs7QFhKsiwG5pWkk2U1kJ1iJNtq6GISkFLS2OZHo4Frj5vIzKL4aKA77riD6upqnnnmmVF3gY1GlkzIYn+zi7UHWlly/Cl89Ow8DrzxMP+duYyLjp2JThOpbvuvzw4SDMshj7N/6KGH2L59Oy+99FKPSkBKyauba3H6gwSCYT7b38yUvJSYEmhx+Um3GLhsSemYOkfVjGAU8cYbb3DWWWdxzz338P3vfz9unS8Y4qNdjby7owFBpGzAYGZJuv1Bmp1+JualcPmSMjI7OdiklFx99dWYzWYefPDBQduvIrF4AyEeeH8Pdk+AQGMF93zrfMzTj+Poq+/grNmFaKKlKeodPr40vygu9DKZNDc3M2XKFBYtWsTbb7/doww7ax387eN9GHUaXlhfFcu3Meq0uHxBXP4g3ztlaixyaLTR04xAKYJRxtlnn82nn37K7t27yc/P77K+1eXn7e31fL6/GaNOS06K4bAvXCklNk8Aly9Sh+Wk6XksnZjda+JNKBRSkUIjjAaHlz++s5s0k54Pn/oz7zz5V/Iu+QULl6/g1Bn5CCEIhMLU2r2cOiOfM2YXoElyKYaf/vSn/Pa3v2Xjxo09RgpVtrj528d7CQYlr26pBeDSxaWkmHT4g2Ea2rxcf/xEpnaoODraUIpgjLBr1y5mz57NFVdcwaOPPtrjuKpWN69sqol0oIr2EU4z6/uMAAmFJS5/EJcvkoE6IcfKSdPzmZqf0qMC+M9//sPs2bN7LQGsGN5srrTx6Kr95Ft0/OtXN5N/7JfZpytjfmkGx0/JQQhBMBymxublmEnZXLCgOKmZuF6vl48//pjTTjut2/Xba+w8tuoAOo2Gt7bX0eYJcsniEnJSjBFfR6uHixeXsHxyTrefHy0oRTCGuP3227n33nvZuXNnr026IzZ9HxUtbrbXONhV34Y/eCjBS0b/kUTC8QSRPIOSTDOT8lKYV5pBUbqp1xlFeXk5c+bM4ZxzzumzV6xi+CKjpcQ/3NVIcaYZAXy8pynWaGjFtFz0Wg3hsKTG7mFctpULFxZTkmlJqFwejwefz0dGRkaPcq8sb+KF9dVkWvS8s6Oe6lYPX5ofcQ5LGXF4Hzslh/PnF496v4BSBGMIv9/P3r17mTFjxoA+FwpL6h1e7J4AwZDEHwoTDIXxB8OkW/Tkp5nIthr6/aQnpeSUU05h7dq17Nixg6KiosP5OophQiAU5tm1law90EquCd558s9oJx1DucwnO8XA2dEy41JKWlz+WLmGU2bmYzEkJi7l29/+Nm+88QabN2/u0tQoFJa8ujmivNLNej7a3cj+JhenzsxnZmFabAYzrySdK5aOG9W1hNpRUUNjCIPBEFMC77zzDieddFK/7PJajYhlEw8Gjz32GO+//z4PPvigUgKjAL1Ww2VHlZFm0vPGur2se+9l0tav4vw7H+PdXc089XkFp8zIZ2p+KtkpRoLhMB/vaWLtwVYuWFDMvJKMQfUdvPzyy/z1r3/l+9//fpwSkFJysNnNm9tq2VXnJBhVYL5AmBVTc5lZmIY3EKKhzcfps/I5dWbBqC0v3V/UjGAU8/HHH3PCCSfwm9/8hh/96EdJ3ff+/ftZsGABc+bM4aOPPkKjGf1PW2MFKSWf7Gni7r89zqt/vI0TL7mGFVfewhtb66i1e5lbks4xE7MxRjvRuXxBWlx+slMMrJiWy9ySjC61rgZKTU0Nc+fOpaysjNWrV2M0Rqre7mty8fa2OsobnOg0gh21bWyutpNtNXD6rAJyU43YPQHc/iBfXVLG/GGQ/5BMlGloDCKl5NJLL+W///0vq1at4qijjkravoPBIHfeeSfXXnst48aNS9p+FcljY0UrV117PZvefo6v/OC3LDzlfFaWN7Gh0oZBq2FOSToLSjNi2eYuXxCbJ4AAFo/P5OgJ2ZRkmgdskgmHw5x22mmsXr2a9evXUzhuIgeb3Hyws4EDLS4sei1t3gCf7Gmmxe1nfmkGyydlo9UIGtp8mPWRzPvSrMT6L4YjShGMUVpbW5k3bx5Go5H169eTmprYVnuBQIDW1lby8kZnb1dFPDtrWjn7rLOo2buD2x97j9S0VBrbfKw90MKeBicajWBGYSoLyzJjbUqD4TAtLj+BYBiNRkNZlpmp+amUZlnITzNhMUQaLnV23EopCYQkdU0tXHrJRSw95TyKlp5NvcOLQKDXCursXjZW2mlx+0kx6jhlRh7jsq04PAHsngAzi9K4ZHEp6eYjm5GMVJQiGMN88sknrFixgvPPP5/nn38+oZERN954I6+88gqbN28mM3NsTbvHKrUNjbz+2Tb2BLIIhyV5aSa0GoHN7WfdwVZ21LYRkpLsFAMTo5ns+WkmNELEhSO334oOlTvXYDXq0Gs1uHxBnL4QoVAYjUYgw2GEJlIOpc0b4ECzm63VdrzBMLmpRhaUZjAlP4VgSNLo9JGbauTCBSVMzU8Z9ZFBvaGcxWOY4447jvvvv5+cnJyEXgR/+ctfePDBB7ntttuUEhhDFOblcs15K7B7Atzxx4c4WLoIizWFHKuBk2fks3RiNjvr2jjQ5GLtwVa+ONCKWa+lONNMlsVApkVPhtVAllUfK/ImpSQYjrw8/hBajaBmy2o+fv7vnPbdP9Aa0lNj89DQ5ovVv5qUa2VBaSZFGSa8gTD1Dh9GnYaLF5Zw1ISsYVcldTihZgRjkHZ/gV4/eNPjt99+m7POOoszzzyTl156SWUPj0F27NjB7NmzOfWMs/jGHQ+wvspBMBQmxaQnzRTpX+ENhKhocbOvyUWd3YvDG6DjLUirEWiFiPwfffmCIVp3fUHDi3dhyCkj79K70FtSyUs1UZRhivbMMGHQaWhx+QmGJKkmHcdOzuGoCVlH7JgeTagZgQKIRPOsWLGCc845h2eeeWZQlMF7773Hueeey4wZM/jPf/6jlMAYZcaMGdx7771897vfJeB18+gTT9IUNPLR7kaqbW4kYNRqKc20MCUvYqIJhSV2T4BWt59Wtx9vIEwoLONett1fsOe/vyKreCIX/uRBsrKzyYnWsnL6grgDIWzuAELAgtJMlkzMYnz24fXYGKuoGcEY5IEHHuDmm2/mggsu4Omnn44rW304uFwufvSjH/GLX/yCrKyhrUCpGHoee+wxbrjhBvLy8njxxRdZtGgRDW0+Klvc7Gt0Ud7opMXlRxNtgRqWEo0Q6LUadNpIBns7+zat5qm7vkN2yUQu+unfsKamR1qnAnqNhnHZFibnpVCcaaY00xKLUFJ0j5oRKGLcdNNNANx8883MmzeP++67r8caLT1RWVnJT37yE/7617+SkpLCAw88kAhRFSOQq666ijlz5nDRRRdx8OBBFi9eTH6aifw0E4ujpard/iANDh9OXxBPIITd48fuDuDwBglLGesyZp4+lXnLVvCDu/7I+KJ8Uox6Ukw60s16Mi36Me34HUyGxYxACHEGcB+gBf4hpfxtb+MHa0bw0oZq7n5rFzU2D0UZZm47fRrnLyiOW1dt86AVgpCUFEfHAPzilW20Rnu2CiJPKO3j2v/XCGhv+NXeAL640356kuXE6bl8sLOxW9mO5Ht15NVXX+XWW2/lhhtu4NZbb+3Xtm02Gy+88AI/+tGP8Pl8vP322yxduvSw9n+4vLShmjtf3obNEzn+Rp0GfzAc+w2WTsxkW01bbH2mRc/Zcwt5dVNtbBlAhlnPnefNAuj1t+78XSD+97foNQRCYQKHyjRh0EZ6P7e6A1222d2x6O18O39BMT99aQtPrqmI2dMteg0LyjL4bF8roQ7XcG/76Gmfg3HOdYfH48FsjmSpv/nmmxQUFDBv3rxeb94+n4/77ruPDz74gNdff/2wbvSJOP/62mZvxzLdrEcIsLkDFGWYGZ9tZtXellh0lNWg5VcXzAGIO68t+kgIrcsfiu0n06Ln5+fOOuzvM2zDR4UQWmA3cCpQBXwBfEVKub2nzwyGInhpQzW3v7gFT+DQQTbrtfzmwsgP0nldO3qNIAyxSIXDoX0/HZVOT/vr6TM90dv36u6zfr8fIQR6vZ4nn3ySTz75hOOPP57S0lLKysooKipCr9cTCoW45JJLeO211/D7/cydO5dnn32WadOmHdH+B8pLG6q57blNBI7g+HdEIyLKo7vt6TUCBARCMm7Zkfz+3R2L3n5/s17LwrJ0Vu5tOaJ9dGYwz7m+qKqqYurUqXg8HqZNm8Zll13GpZdeGiuD0tzczKpVq9i7dy/3338/+/fv59xzz+Xf//43aWkDKwmdiPOvr23251j2hSb6NBnuc2SkO+DdF887rO8znHsWLwHKpZT7pJR+4GngS4ne6d1v7eryw3kCIe5+a1e369oJRB1YR0L7fnqTpa/P9ERv36s7DAZDzGG8d+9e/vGPf3D55Zdz/PHHM378eEpLSwHQarWkpqbyne98h88//5yNGzd2UQKHs/+BcvdbuwZNCUBkxtbT9gJhGacE2pcdye/f3bHo7ff3BEIDUgI97aMzg3nO9UVJSQkVFRX87W9/o6ioiF/+8pfMnDmTF154AYDNmzdz3nnnccstt5CSksI777zDyy+/PGAlAIk5//raZn+OZV+E+6kEIPJgMljXUzvDwUdQDHTsNF0FHN15kBDiOuA6gLKysiPeaY3NM6Dlg03H/fR3n/0ZdyTf64477uD73/8+lZWVVFRUUFlZid/vj61//PHHE7r//pCs3yeRdP4OifhOfW1zMM+5/pCTk8N1113HddddR21tLc899xw5OZHa/4sXL2bNmjWUlZWRn59/RHb/RJx/fW1zKM7Jwd7ncFAE/UJK+TDwMERMQ0e6vaIMM9XdHMz2ypvdrRtMOlb47EmW3j7T25jevldfWK1Wpk+fzvTp0/s1frD3f7jbH0l0PhaJ+E59He/BPOcGSmFhITfffHPsfWpqKkuWLBmUbSfi/Otrm0NxTg727zIcTEPVQGmH9yXRZQnlttOnYdbHx7ub9VpuO31at+va0UeTXI6E9v30Jktfn+mJ3r5XMkj0/m87fVrEdj9IaAQ9bk+vidSv6bzsSH7/7o5Fb7+/Wa9l+aSBheT253gP5jk3nEjE+dfXNvtzLPtCI/p/M9ZrxaD/LsNBEXwBTBFCTBBCGIDLgJcTvdPzFxTzmwvnUJwR6bZUnGGOOX86roOIM5HomLsvmccfLplHpuVQIlb7baF9XPv/He8X7bPdjvvpTZYrlpZ1K9uRfK9kkOj9n7+gmLsvmUdGh6JhRp0m7jdYPikrbn2mRc8VS8vilkEkauiPX57P3ZfM6/G3vvvieXHfpbvf36LXoO90JRm0Ijam4za7Oxa9nW+/uXAOT167jCuWltHRYmLRa1g+KSs2tp3+Hu/BPOeGE4k4//raZl/HMqM91DW6bvmkrLhcCatByx+/PJ8/Xjo/7hy16DVYDfEKJtOiP2xHcW8MedQQgBDiLOBPRMJHH5FS/qq38SqhTKFQKAbOsE4ok1K+Drw+1HIoFArFWGQ4mIYUCoVCMYQoRaBQKBRjHKUIFAqFYoyjFIFCoVCMcYZF1NBAEUI0AgcP8+M5QNMgijNYKLkGhpJrYCi5BsZolWuclDK388IRqQiOBCHE2u7Cp4YaJdfAUHINDCXXwBhrcinTkEKhUIxxlCJQKBSKMc5YVAQPD7UAPaDkGhhKroGh5BoYY0quMecjUCgUCkU8Y3FGoFAoFIoOKEWgUCgUY5xRpQiEEGcIIXYJIcqFED/qZr1RCPFMdP0aIcT4Dutujy7fJYQ4Pcly3SqE2C6E2CyEeE8IMa7DupAQYmP0Najlufsh11VCiMYO+/9mh3VfF0Lsib6+nmS57u0g024hhK3DuoQcLyHEI0KIBiHE1h7WCyHE/VGZNwshFnZYl8hj1Zdcl0fl2SKEWCWEmNdh3YHo8o1CiEEt59sPuVYIIewdfqs7Oqzr9fdPsFy3dZBpa/R8yoquS+TxKhVCfBC9D2wTQny3mzGJO8eklKPiRaSE9V5gImAANgEzO425EXgo+vdlwDPRv2dGxxuBCdHtaJMo14mAJfr3t9rlir53DuHxugr4czefzQL2Rf/PjP6dmSy5Oo2/iUjp8kQfr+OBhcDWHtafBbxBpD3FUmBNoo9VP+U6pn1/wJntckXfHwByhuh4rQBePdLff7Dl6jT2XOD9JB2vQmBh9O9UYHc312PCzrHRNCNYApRLKfdJKf3A08CXOo35EtDeePd54GQhhIguf1pK6ZNS7gfKo9tLilxSyg+klO7o28+IdGlLNP05Xj1xOvCOlLJFStkKvAOcMURyfQV4apD23SNSyo+B3rrIfwl4Qkb4DMgQQhSS2GPVp1xSylXR/ULyzq3+HK+eOJLzcrDlSsq5BSClrJVSro/+3QbsINLPvSMJO8dGkyIoBio7vK+i64GMjZFSBgE7kN3PzyZSro5cQ0Trt2MSQqwVQnwmhDh/kGQaiFwXRaehzwsh2luKDovjFTWhTQDe77A4UcerL3qSO5HHaqB0Prck8LYQYp0Q4rohkGeZEGKTEOINIcSs6LJhcbyEEBYiN9MXOixOyvESEZP1AmBNp1UJO8eGRWMaRQQhxBXAYuCEDovHSSmrhRATgfeFEFuklHuTJNIrwFNSSp8Q4nois6mTkrTv/nAZ8LyUMtRh2VAer2GLEOJEIorg2A6Lj40eqzzgHSHEzugTczJYT+S3copIh8KXgClJ2nd/OBdYKaXsOHtI+PESQqQQUT7fk1I6BnPbvTGaZgTVQGmH9yXRZd2OEULogHSguZ+fTaRcCCFOAX4CnCel9LUvl1JWR//fB3xI5EkhKXJJKZs7yPIPYFF/P5tIuTpwGZ2m7gk8Xn3Rk9yJPFb9Qggxl8jv9yUpZXP78g7HqgH4L4NnDu0TKaVDSumM/v06oBdC5DAMjleU3s6thBwvIYSeiBJ4Ukr5YjdDEneOJcLxMRQvIrObfURMBe1OplmdxnybeGfxs9G/ZxHvLN7H4DmL+yPXAiIOsimdlmcCxujfOcAeBslx1k+5Cjv8fQHwmTzknNoflS8z+ndWsuSKjptOxHknknG8otscT8/Oz7OJd+R9nuhj1U+5yoj4vI7ptNwKpHb4exVwRhLlKmj/7YjcUCuix65fv3+i5IquTyfiR7Am63hFv/sTwJ96GZOwc2zQDu5weBHxqu8mclP9SXTZL4k8ZQOYgOeiF8bnwMQOn/1J9HO7gDOTLNe7QD2wMfp6Obr8GGBL9GLYAlyTZLl+A2yL7v8DYHqHz34jehzLgauTKVf0/Z3Abzt9LmHHi8jTYS0QIGKDvQa4Abghul4Af4nKvAVYnKRj1Zdc/wBaO5xba6PLJ0aP06bob/yTJMv1nQ7n1md0UFTd/f7Jkis65ioiwSMdP5fo43UsER/E5g6/1VnJOsdUiQmFQqEY44wmH4FCoVAoDgOlCBQKhWKMoxSBQqFQjHGUIlAoFIoxjlIECoVCMcZRikCh6AUhRHaHapR1Qojq6N9OIcRfh1o+hWIwUOGjCkU/EULcSaS66T1DLYtCMZioGYFCcRhE6+m/Gv37TiHE40KIT4QQB4UQFwohfh+tXf9mtHQAQohFQoiPokXL3opWjlQohhylCBSKwWESkYJ85wH/Bj6QUs4BPMDZUWXwAHCxlHIR8Ajwq6ESVqHoiKo+qlAMDm9IKQNCiC1Emqu8GV2+hUhtm2nAbCJVK4mOqR0CORWKLihFoFAMDj4AKWVYCBGQh5xvYSLXmQC2SSmXDZWACkVPKNOQQpEcdgG5QohlECk53KEZi0IxpChFoFAkARlpu3gx8DshxCYi1SWPGVKhFIooKnxUoVAoxjhqRqBQKBRjHKUIFAqFYoyjFIFCoVCMcZQiUCgUijGOUgQKhUIxxlGKQKFQKMY4ShEoFArFGOf/A1Tpoyig9YhVAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}