{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 11.1\n",
    "## Air-shower reconstruction\n",
    "Follow the description of a cosmic-ray observatory in Example 11.2 and Fig. 11.2(b).\n",
    "The simulated data contain 9 × 9 detector stations which record traversing particles from the cosmic-ray induced air shower.  \n",
    "Each station measures two quantities, which are stored in the form of a map (2D array) corresponding to the station positions in offset coordinates:\n",
    "- arrival time `T`: time point of detection in seconds\n",
    "- signal `S`: signal strength in arbitrary units\n",
    "\n",
    "The following shower properties need to be reconstructed:\n",
    "- `axis`: x, y, z unit vector of the shower arrival direction\n",
    "- `core`: position of the shower core in meters\n",
    "- `logE`: $\\log_{10} (E / \\mathrm{eV})$, energy of the cosmic ray\n",
    "\n",
    "Reconstruct the properties of the arriving cosmic rays by analyzing their\n",
    "air showers:\n",
    "\n",
    "### Tasks \n",
    "1. Set up a multi-task regression network for simultaneously predicting shower direction, shower core position, and energy. The network should consist of a common part to the three properties, followed by an individual subnetwork for each property. Combine the mean squared errors of the different properties using weights $w_j$.  \n",
    "\n",
    "\n",
    "2. Train the model to the following precisions:\n",
    "- Better than $1.5^\\circ$ angular resolution\n",
    "- Better than $25$ m core position resolution\n",
    "- Better than $10\\%$ relative energy uncertainty $\\left(\\frac{E_\\mathrm{pred} - E_\\mathrm{true}}{E_\\mathrm{true}}\\right)$\n",
    "\n",
    "  Estimate what these requirements mean in terms of the mean squared error loss and adjust the relative weights in the objective function accordingly.  \n",
    "\n",
    "3. Plot and interpret the resulting training curves, both with and without the weights $w_j$ in the objective function.\n",
    "\n",
    "\n",
    "##### Hint: using a GPU for this task may be advantageous!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "keras version 2.4.0\n"
     ]
    }
   ],
   "source": [
    "from tensorflow import keras\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "layers = keras.layers\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import gdown\n",
    "url = \"https://drive.google.com/u/0/uc?export=download&confirm=HgGH&id=1nQDddS36y4AcJ87ocoMjyx46HGueiU6k\"\n",
    "output = 'airshowers.npz'\n",
    "\n",
    "if os.path.exists(output) == False:\n",
    "    gdown.download(url, output, quiet=True)\n",
    "\n",
    "f = np.load(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input 1: arrival times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200000, 9, 9)\n"
     ]
    }
   ],
   "source": [
    "# time map\n",
    "T = f['time']\n",
    "T -= np.nanmean(T)\n",
    "T /= np.nanstd(T)\n",
    "T[np.isnan(T)] = 0\n",
    "\n",
    "print(T.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot four example arrival time maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEYCAYAAAAaryJBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABjsUlEQVR4nO2dd3gc1dWH37O7apbkbsvG3caAAdNsTDMgYyAUmw6BhCRACAFCQkIICalAPkIICUkoCTiU0EIJ1YTeBDbFuGNjuWMbGxvLTdZK2n6+P3ZXqMyuts+udN/nmUfa2Tlzz52Z356ZO/eeK6qKwWAwGAy5xmG3AwaDwWDonpgAZDAYDAZbMAHIYDAYDLZgApDBYDAYbMEEIIPBYDDYgglABoPBYLAFlx2F9u/fX0eOHJnVMhobGykvL89qGcaHr5g/f/42VR1g9d3XppTr9h3Bttt/4n1NVU/KumNdCKObruVDPM1AaroRkQeAacBWVd0/Y85mCVsC0MiRI5k3b15Wy6ipqaG6ujqrZRgfvkJE1sf6btuOAB+8OqTNutI9PuufdacAETkrgc08qvpy1p1JE6ObruVDPM1Ayrr5N3AX8HCavuVENxkLQCLiBOYBm1R1Wqb2ayh8FAgQ7HS7LPEv4AVA4mxzDGBLADK6McQiFd2o6nsiMjIDxedEN5l8AroaqAV6Jmvodrtxu91UVFRQUVGRVbtc2gQCgRa7RMn3OqVynhTFr6GEt88wr6jqJfE2EJFHc+WMBTnVTSFoLVndFEKdjG6syUgAEpGhwKnAzcA1ydi63W5qa2sREVSVcePGJXSSUrHLtY3P56O2trZL1SnZ8wQQAjxqzxOQql6YiW2yQa51UyhaS0Y3hVKnDOqmv4i0boedoaozEtphEuRKN5l6AvobcB1QGWsDEbkMuAygqqqKmpoaAAKBAD6fD6fTSTAYpL6+Hperc7c6s3O73S1lJGqTSjnxbAKBAHV1dRmrUyo2dh6HKIriJ/9yDorIIFXdYqMLfyOHuknEptB00xV/P6LE0M02VZ2Y0A6yRCZ1k3YAEpFoj4v5IlIda7tIlJ4BMHHiRI2+4MvWHYzVS8Rc3/XU1dUxYMAAW+/K7DwOUVTBn3/xB+B+wk8gOccO3SRiU2i66Yq/H1G6g24y8QR0FHCaiJwClAI9ReTRRB/PKioqGDduXNJtpKnY5dqmvr4+qQuuEOqUWlu24NH8G3KmqrYEnwg5102haC0Z3RRKnXKlGxF5HKgm3FS3Efidqt6f1E468yuDukk7AKnq9cD1AJE7uWuTbRtM9sSkY5dLG5fLldf+5cIGwr15/DYHIBEZbrVeVTfk2pdIubbophC0lqxuCqFOudKNql6QdEFxyLZubBkHZOhehBB8OO124yXCmhbCTxyjgBXAfnY6ZTDEojvoJqMBSFVrgJpM7tNQ+OTDE5Cqjm/9WUQOAa60yZ02GN0YrOgOujFPQIasE27LLrLbjTao6gIROcxuPwyGWHQH3ZgAZMg6qoJf7W1KEJHW42wcwCHAFza5YzB0SnfQjQlAhqyjCD5N/lLLcGLF1mNtAoTbtp9Jc58GQ9ZIVTcZJqu6sb12hq5PKPWmhH+TgcSKAKp6Y7r7MBhySRq6yRjZ1k3+Dc4wdDnCL1NdbZaE7FTfA3Zky69IlgGDIS9JVTfZJpO6yYsA5Ha72bJlC263O+t2ubSJJlVMhnyvUyrnSQm3ZbdeiOS0arXYEQziZfrNe3J5DvNZN4VQpwzqJh/ImG5sD6ldMZmgSUbalpAKnlCHpgTbc1qp6r12lp8OJhlpYdUpg7qxnUzqxvYnILfbjYhQWVmJiCR8l5CKXa5tnE5nl6tTsucJ8udOTkROFZHrROS30cUWRzJALs9hPuumUOpkdGON7U9AFRUVqCoNDQ2oalK5nJK1y7VNMBjscnVK9jzBV23ZdiIi9wA9gCnAfcA5wMe2OpUGuTyH+aybQqmT0Y01eRGAuloywVSSKhZKnVJPqphSN+xMJlY8UlUPEJFPVPVGEfkL8EqK+7KdXJ7DfNZNodQpl7rJMFnVje21g66ZTDCVpIrplJWvNhAZUBdK/lLLcGLF5sjfJhHZA9gODM7g/nNOLs9hPuumEOqUS91kmKzqxvbaGbo+0bZsm/mfiPQGbgMWEG7h+JetHhkMcegOujEByJB1Qghem+/kVPX3kX+fEZH/AaWqWm+nTwZDPLqDbmzvBWfo+oRndrSnN08ke287f9TbWkRW2xgMdtMddGOegAxZRxECIduaEh6MTPgWb/Dc/cDBOfHGYEiQ7qAbE4A6IeAPUrvkc0IhZd8Dh1FU1PkhU1VWbK6j0eujvslDrx6lCZX1RdMu1jTUMbKiH8PK+6bret6g9jYl9ALmE19IdTnypdvgafKy7IMVlJaXsM9hY3E4Om9sCYVCLP9kI82NXjzNPkrLinPgaf7SHXSTdu1EZBjhZJFVhF9QzVDVv6e733zgk/nruOnaJwgGQwCIwK/+eB4Tjtgzps2X9W4uv/9ZNu7YzXf3reK4/5vB946bxOXHHx7Txh8Kcv2CZ3l78wqKHE78oSCHDxjF7YeeR6kz/0ZCJ0u4KcGe1l5VHWlLwZ3QlXXz5qPv8vcr/oXD6UBVKe/Zgz+8/EtGjR8R02bt8s385oqHaWr0MO2SvTj/6D/woxvO4LhpB+XO8TyjO+gmE7ULAD9V1X2Bw4EfiMi+yewgH3M5Nexu5jdXP0rD7maaGr00NXppdHu58don2LUjtu0P//0Ca7fuoNnnJ6SKNxDk/nfmUrNsTUybf66o4Z3NK/CGArgDXryhAB/WreXPn76e0TrZYQNfNSW0Xgz26CbbWlv36ef87fsz8DR6adrdTHODh22bdnDd8TcR8AcsbQL+IL+49AG2b91Nc6OPUFDxNPv52++eZ/3qrbbXyQ4b6B66STsAqepmVV0Q+b8BqAWGJGofzZW0ceNGamtrk0qlkaxdMjaz3vzUcn0opNS8ttTyuw3bd7Fm6w6CIW2zvtkf4JHZC2OW9cRn8/CE2orTFwry7PqFqGoMqzDZPg7p2ERRwBdytVm6O3boJhdae2nGG/h9HQON1+Nj4VtLLG0WfLiagC/YYX3AH+DVZ+ZmzLd07YxuMo909gOX1M5ERgLvAfur6u52310GXAZQVVU14YknngAgEAjg8/lwOp0Eg0GKi4txuTo/0J3ZRUceJ2PTmp3b3Wyva7AMAP0GVNK3f2WH9R5/gHV1OwlFbAaUFVHX7AegpMjFmIHW73Vq6zcTslgvwL694o/5yvZxSNRmypQp82MlF+2zz0A97oFz2qx79qh/xty+u5Er3SRik+71suWzrTRYtBCIQ6gaMYDKvh0HZDbUN7N18y5CwbBu+gwsZedWDwCVvcsYNKRPWnVKxS4XuomnGegeuslYSBWRCsIz5f24vYgAVHUGMANg4sSJWl1dDWQvm21NTQ3RMhK1ac2q2i/46aUP4PX426wvKS3ilru/zX4HDe9g4wsEOPrGe2n0+gC48sAh/GPxJopdTi6tPpTq6iMsy3rk/YeYs+0z2oe6cZVVXFkdPxlAto9DOjZRFAiE7O3xLyICfBMYrao3ichwYJCq2poPLpe6ScQm3evl7c2zeeyGe/A0etusLypx8fCau+m/R8ebsLot9Vxy8u0tT05nX7kfz/zjU0rLirjm/87mmOrxadUpFTujmzDZ1k1GApCIFBEW0WOq+mwytvmay2nsuD2YfNw43n+nFk/kKaa0tIiJR45l3wOHWdoUu1z8+ozjuPHZN/EGwmIqcTnpV1nONyfH7q14/fiT+eas+/EGA/g1iAsHRU4nNxx8WkbrlGubKOG2bNuHnP0DCAHHATcBDYSv2UPtcijXusmF1o4553Cev/MVPvtkPZ6mcBAq6VHM2ddMtww+AAMG9eLsi47i+Uc+aNFaSWkRo/cZzJFT478Wy9ffj3T9g+6hm0z0ghPC/cFrVfX2VPaRr7mcrr3xTI6oHsdrLyxAVTlh+sEcc/y+hKtszfQJ4xg5sA+PzV5ID5eHH5x4JOcdPp6K0pKYNnv2HMjzx13Jo2vmsHTXJvbpOYgLxxyWcFfs/M9pBT77X6AepqqHiMjCsE+6U0Rs6+drl26yrTVXkYs/v3MDbz36Hu888T49epYx7fsnMvHEA+PaXXT1iew/YRQvP/UxPSpKuOKX05g6/SBcRZ1fN/n6+5GODXQP3WTiCego4FvAEhFZFFn3S1V9OQP7thWHw8HRU/fl6E7uwtozftgg/njByZHH+MSaaweX9eJn+5+Yipt5T57cyflFxEm4ZQMRGQCWr95yRZfVTXFJESd/dyonf3dqUnYTJ49l4uSxSemmK9MddJN2AFLV2RT41MaG7BO0aTxDK+4AngMGisjNhOc1+bVdzhjdGBKhq+um6/XrM+QdqoI/aG9Tgqo+JiLzgamEf/jPUNVaW50yGOLQHXRjApAh6ygQtL8pAeBLYBbh675MRA6JjsUxGPKN7qAbE4AM2UchqPa2NonI74GLgDXQ0uNdCffuMRjyj26gm24TgBp3NzPz3+/x8Zuf0rt/JWd891gOPGqvuDahkPLqx8t5YfZSVJVpR+7HKYePw+WMf1cy98uN3Fc7l0m7QyxaNJuLxk2gd0lZXJtdvi18tP0ZtjSvYmDpKA7rdzb9SoYmXc98RBECQdvv5M4Dxqiqz25Hujoen58X3l3KG3NW0KO0iHOmHsTRB42O23tUVXm7di1PzFnMkZXCxg8Wcu6h4ynpJPnvp9u2MmPRXNbV72TS4KFceuBEqsqT73GWj3QH3XSLANS4u5kffO1WdtbtxucJj89ZOGsFF18/ndMvOdbSRlW5fsZLvL/kM5ojA+Q+Xfclb81fxd9+eHpMMT2xajE3fvwmzcEA+5eP4J9LP+KJ1Z/wyrSL6VNqHYS2ej7jkXXXEgj5CBFki2c1y+prOH/EHxjaY1wGjoC9qELI/qaEpUBvoPPkYoaU8QeCXHbzk6zbvANPRDeLVn3BOVMO5EfnHxPT7rZX3uPJjz+h2R/goP2GcP/rs5m5qJZHL/s6xS7r9yDvrF/Lla/PxBsMElJl2batPLV8Kf8751sM69krK/XLJd1BN7bXDrKfTPDFh2a1CT4A3mYfD/5hJk1uj6XNp+u+ZHar4APg8QWYv3IjC1ZusrTxBAPcNPctmoOtygkF2eFp4r7a2Dmt3thyD75QMyHC+bCUEH718trmu+PWK0q+J1UECIakzWIDtwALReQ1EZkZXexwJFPkYzLSNz9ewfotO1uCD4DH6+epNxfy5fYGS5svdu3m8TmLaW6VrNTjD7Cmbgevf7rK0kZVuf7d12kOBFpSX/lDIRq8Xv7y8eyM1skOmyhdXTe2PwFlK5VGa+a8sbRN8InidDlYs3Qj4w/vOL3C/BWf4w927O7u8fqZu3wDE/bu2Dy2YmcdDosnI18oyNsbV/Ozg63vADc1L7dcv9W7lpAGcUjsnjC5Sg+SXkoRyYeXqQ8BtwJLsHf8T0bI5TlMxu79xZ/R7PV3WO90Oli0ahNf67dPh+/mr9uE0+EA2iYkbfb5eXf5WqYd2NFma1MjOz0dbx5DKLM3rs9onXJtEyVV3YjIScDfASdwn6r+MemdfEVWdWP7r4Lb7UZEqKysRESSymabqF2fAR0ThwIEAiF69im3/K5XeRnFFu96ilxO+lT2sC6npIxAyPoc9S+1Lgeg2GHdNOekGOnkFKVy/HJl04KChqTNYgNNqnqHqr6jqu9GFzscyQS5PIfJ2PXrVY7TYXV+lV4V1hMz9u5RhlWLtlOEAZXWuikvKo6ZKb53aecTQHZV3UQGjd4NnAzsC1yQ7DQf7ciqbmwPQBUVFagqDQ3hzNPJ5HJK1O6MS6spaTe7osMhDB7RjxF7W2ebnjphLGIhJIdDOPHQvS1thlf2ZlyfgTjbqanM6eK7+8ZOnTSxz3Rc0jZVj0Nd7F9xfNwXt5Da8cuVTWtCIWmz2MAsEblFRI4QkUOiix2OZIJcnsOktFY9HpfFO5uykiImjuuYwBfgiDHDKS0q6jAqt8jp5NxDrRORVhQXc/zIMRQ72pZV5nRx6QGdZ1HowrqZBKxW1bWRjgNPAKcnVWhbsqob25vgcpFM8IAjxvLdX53G/Te/gMPpIBgIMXhkf/7vkSti2lT2KOGuq8/i2n+8SLPPD6oUFzn50xXT6VMZu0fbjCln8d23n2blrm04RChxOvnxQZOpHjI6ps2RA85nl/9Llu1+FydFBNXPyB4TOHHo9zN6HHJtE0UVQvb35olmg209NW3BdsPO12Sko/box+8u/Ro3P/AGihIKKX0qy/jbT8+K2XvU5XTw4HfP4YqHnmdHYxMOgfLiIv7v7BMZNSB2PsQ/TTmJy197gbmbN1IkDnyhEBfufxBfHxc7e3Yqdcq1TZQUdTME+LzV543AYcnupBVZ1Y3tAQhyk0xw+kXHcPy5h7Fm6UYq+/RgxF7x59kBOHDPPXj1z9+jdv3WcPvtiKpOu2APKCtn5qnfYXX9dlbMXcC8aWdQWRw7ESmAQ5xMG3IN1QMvYodvE72LB9OzqH9C9YL8T6oIqTW7ZbItW1WnpGqbr+RjMlKA4yftzTEHj6F23ZeUFRcxdviATp/k9xzYj9evvYTazXWsW7qY9391FsWdzLVTUVzMo9PP5fPdu9jkbmDvvv1j9jRNt065tgljqZv+IjKv1ecZkSk7skK2dZMXAShXlJWXsP9hY5KycToc7D9qUNJl7dmrHxtdRZ0Gn9ZUFPWloiixDNgFRaQtOxlatWWfQPgubq6IzFTVZUnu50JVfVRErrF0LcVM1Ib4FBe5OHBswhO8AiAi7LvHQLauLOo0+LRmWM/eDOvZO0kPCwBr3WzrZEK6TUDr+WKGRtYlRa50060CkMFGkn8CamnLBhCRaFt2UgEIiL7FtuqJkrnpgA2GbJC8buYCY0VkFOHAcz7wjRRKzoluTAAyZB8lFSFlpC1bVe+N/Pumqr7f+jsROSrZ/RkMOSMF3ahqQESuAl4j3HT9gKp+mnTROdKN7W+GDd0DDbVdiLRlt1ouy7ILdya4zmDIGyx007mN6suqupeqjlHVm9N0Iau66VZPQAF/gI0rN1PRpzzm9MDtUVU2bdmFKgwd3LvTl6lR6hob8QYC+IJBip2JpVRvDnrY5t1G3+K+lLusxxp19C8Agc/AUYk4k39XlRMUJGhbW/YRwJHAgHbt2T0J3yEaskBIlbU7dtCjqIg9evZM2G5TYz3eYDi7gdWgbit2+xvZ7q1ncFl/Sp22TXKbeax1kxNypZuMBKAMj7zNCm8/Pos7f3A/wWCQgD/IvofvxW+euoZe/WOLY/W6On592wvU7XAjAn16lfP7n05nnz1j/9DXezz86JWXmLNpIz+qGszPZvyT3xxTzbn77R/TJqQhntzwDG98+TZOcRLUIEcPOJJvj/wGzjhZEELNr8Du3wJ+0CBaNB7pfSfi7JfQMckdYmdbdjFQQfhab92evZvw5Fq2UQi6SYVZ69Zx7cuv0uT3E9QQe/btxz9PP40hvWJrbaN7F5fPeoY1u7fzw+IR/Pa5O/jzEadx9OBRMW28QR9/Xv44729bisvhJKQhLhxxIuePOD4b1bKBlHSTKXKim7QDUCZ6K7nd7pT6ySdqVztnFbd/7x68TV8ldF36/nJ+e/qt/P196yfUZo+PH/72CRoavS3rNm+t5+obnuLpey+jstx6tPUPXnqRuV9swh8KEVLF7fNxQ83bDO/Vi8OGDrO0eXXzG7y59R386sev4TQms+s+oMJVwbnDzrS0Uf+nUP9zwNNq3SLY+T2k/7Mxj0WUVI55qucJSDqJRwbbst8F3hWRf6tq5zlacoRdusm21jbs2sUVz8+kOfBV6qvaujq+8eRTvPO971o+1QRDIS5481E2NzcQUiWkSp2nkcvfe5rXTv0eQyt6W5Z1x8qn+WDbUvwawB/Jv/jY+jcYWNqX46o6HyuZq+OXS91kilzpJhPvgNIaeRvNlbRx40Zqa2uTSg+SqN0zt7+Ir7ltNvGgP8jqRevYuPILS5uaj1YRsMgFFwiGeOv9FZY2m3bvZv7mL/C3S8fTHAgwY/48SxuAl7e8ji/U1j+f+nlty5sx041o47+BtjZCAA2sRv3WCRyjpHLMUz1PYWdBQtJmScgsg23Z+RR8IuRcN7nQ2uOLP+lw/YdU2dHczJzPP7e0+Wjrenb5PC1JRaMEQkEeX73Q0qY56OWdrQvxadscj56Qjyc2vJnROuXapoUUdZNJsq2bTDTBJdRbKfKS+TKAqqoqampqAAgEAvh8PpxOJ8FgkPr6elwJjAHozM7tdreUMWJKFedOOrnDPsQpLFn5Cau/WNnhu+ZdjZwzdZBlh8NA4+fU1OzsaBPw88OqPQhF3hZWFRVzzeDwWIjSYKjFn/ZMbDowZqB5990YaZeC+4KOsPjCCc6VIOHXJa2PQ4v/KRzzVM9TFCn49J8ZJ+e6ScQm3etlj/p6fjik4yBvhwhbly2jZu3aDt/t8nm4omgYIVdYA1WOEq4pDWcO6bWpnppdNR1sAhrkrMaDUAuBOhsdMbWWaJ2MbnJDzjohREbrzgCYOHGiVldXA9nL0FtTU0O0jH+//QTP3fYm/nZZeotLi3jyi39R0btjwsPFyzbyj/97Bk87m9ISF3+47gwmHTSyg02T38+1M/7Z0vxwzeAh3L55E0UOB98+8GAuPcZ67qEblv6BNY0dhTmouIqLD/6OpU3IvQzc9wLeNuuVYhwDZyGOPh2OQ5RcZ/VNsRu2gczqJhGbdK+XJz9Zwj3vvNNmagWAYqeT1048keG9e3ewWdewg+tfug9vKKKb0tHc7llLmdPFDQd+jeoxB3awCWmIr3/wO3b52z5RCMLk/uO5ZP/qDjbJ1MnoJjdkIgCl1VspF/mpzvzRKbx831s0bG8g4A+nfC/pUczXf36GZfABOGDcEPbfew+WLN+ENzK3SUmxi71GVzHxAKsnD+hRVMSPDz+Sv330QUsQcjkcVJaU8L0JE2L6d+HI87ml9s/4Q/6WO7oiKeKi0RfGtJEe30Sb/gOhnUA4SCqlSPnFLcEnFrnOaQX23cmJyJ3EGTinqj/KoTutyblucqG108ftw7/mzmPT7t34gmGtlblcTB+3j2XwARhZ2ZfpI/blpQ21NAfD13Kxw8ke5b04beR+ljYOcXDlnmfylxVP4A2FbRwIJc5iLhp1SkbrlGub1nR13WQiAKXdWynb+al69e/JvQtv44lbn2fOSwvoNaAn5/50OpPPjD2uUUS47Zdn8dxri3jprSWEVDllyv6cdfLBOCzTzYf53oSJjO7Tlxnz51LscHHB/uO58tDDGBhnmuA9K0Zzw36/5PlNL7K+8XOG9NiD0/eYxuiKkbH9c/SC/jNR9wzwvg2OPjjKL4KSkzo9HpDbnFaiIMHOt8sS0ZdvRxFOT/9k5PO5JJ9VIZPYoptsa620qIjnLvwG982dx8srVlJeXMyFBx/EWfvFnxHgj4efyqEDh/HIyvmU+Fx8f98juHTcYZQ4Y/9ETak6hL4lPXl8/Zts9mxnv56j+OaIExjSY0BG62SHDXQP3aQdgDLVWynb9KnqzRW3X8QVt1+UsE1RkZPzpk3gvGmxn16smDp6NFNHj6ampobL2j3Gx2JYj6H8cGzs7NxWiKMv0vMXwC+SsrMFm5oSVPUhABG5ApisGn5rLSL3ALNscYrC0U0qVJaU8JPJR/GTyYkPmHeIcO6YAzl3zIHh5q8DYk/f3ZoDe+/Jgb07TijZZejiusnIOyBVfRl4ORP7MnRN8uBlah/Cg+h2RD5XRNbZhtGNoTO6um66VSYEg03Y25QQ5Y+E57Z/BxDgGOAGWz0yGOLRDXRjApAhJ4jNeadV9UEReYWvujr/XFW32OmTwdAZXV03JhmpIftouCmh9ZJrJJzE73jgQFV9ASgWkUm598RgSJBuoJtu8wTk9wepeXsZcz5cTd9+FZw6/SBGjIzfW0ZVmbNuIzMXL0MVph2wD0eOHt5pQtINjXU8v3EO/Zt9vLRpLscPOogSZ1Fcm0BwF9vcj9PkW0xZ8Tj6V3yTImfis6LmM0JeNCX8g3Bik+OAm4AG4BngUDud6oqENMQn9fNZuPNjShxlHNW/mhHlsaekj7LTs5QNDc/RFBjNF41vM7jHsUicXIgAGtyMNj0RTshbfChSdibiSGX20fyjO+imWwQgr9fPj696hM83bMfj8eNwCC+9uJCf/WIa1cfF7h56y6vv8t8FS1oG1b26bBXTxu/D70+Lnexw1tZl/G7JfwiEgnwzsD/3L3+BJ9bP4t5JP6CHy3p2VG9gPcu3nEZIm1H1UN/8Flt3z2CvqmcpK947vcrnA5oXL1MPU9VDRGQhgKruFJEulDo5PwhpiH+s/jNrGlfgDXkRhDk7ZnPaHucytapjNpIoq3Y9zIqd/yCoPghdyoKtd9GvdAKHD7oDEeuGGvUtQHdeAuoH/OB9F238F/R7Lg8T8qZAN9BNXjTBud1utmzZklyepCTsXv7fIjas34bHEx6wFgopXm+Av/zpZXzegKXNqq3beGr+kjYjupv9fv63pJYlm6ybQAOhIDd/+hTekJ9gJIugJ+RnU/N2nvn8g5j+fb7jBoKhelTDiUUVL0F1s2HHL+PWK0oqxy9XNi2E2i25xx9JAKoAIjLANk8yRC7PYaJ2n9TPbwk+AIriVx8vfPEUDf7dljbe4A6W77yLoHqInpKgNrPdM58tTe9Z2qgqWv9z0CaiA7GhGULbUPcdGa2THTYtdHHd2P4ElK1UPK2pebsWr1WgEVixYjPjD+iYpfq9VesIhjoeZ28gwLurPmP8kI5TMqx2byZoMWuUNxTgrS8/4Vujplj61+CZTcdzqjT65qMajNsMkav0IOmmFMmDpoQ7gOeAgSJyM+GU8r+216XUyeU5TMZu4c65LcGnNU6crGhYxsS+h3f4rq55DoKL9sl1g9rMF41vMri8umNBoW0Q3GzhQQC8bwI3ZqxOubZpoRvoxvYnILfbjYhQWVmJiCSVoTdRux7l1k+MoVCI0lLrdzM9iotwOjseHqfDQXmx9f7KnMUtiUjbU+60bn4DcEis75x0dopSOX65smmN3S9TVfUx4DrgFmAzcIaq/jf3nmSGXJ7DZOzKnGUIHd+RKkpJDA04xdoGHBTFep8jJcTMFCNlMf2LYnSTGNnWje0BqKKiAlWloaEBVU0qP1WidqedPsEy0PTq1YM9x1ZZ2nxt37GW17dDhFP2t34vM7zHAAaV9ekgplJHEWcNOyKmf33Lz0VoK05VF5XFp3ba4SGV45crm9bYLSQRuR8oVdW7VfUuVa0VkRty70lmyOU5TMbuyP7VuKSj1hziYJ9K60kZB5YdARbXuYMihleeYWkjjp5QPAltNzmnUgplF8SpTRijmwTLz7JubG+Cy0WCxMOP3JPTz5zAc8/Mxel0oEBZWRG3/On8mD/wfct7cPu5p/DTp1/GEXl8DqnyxzNPoqqndVkiwq0HfYcfzptBg78ZQSgSF9OGHMpxVQfE9G9I75/j8S/H7VsA6gBClLr2YvTAzqfAKYSkipIfL1O/BkwUkb+o6sORdadRoINR8zUZ6fAeozh9j/N4/osncRJ+deAQJz/Y82cUOaxbG5yOEo4YdBcfbvkhqiGCOHBQzH79fkzvkn1iliW9boMd30aDX6CAEERKpiDl1hnkU61Trm2idAfd2B6AIPsJEkWE711+HGeefShLPtlAz149OOjgEZZNbK05bu8xzL72+3ywdgOqypFjRlBREr8DyNAe/Xn66F+waOdavpi/iieOupZBZfEzVzgcpYytepwmXy0e/wpKi0bTozh2wGpPvidVhLwQ0lZgCvCoiBwGXA2W7T4FQz4mIwU4ruokJvU7ihUNn1LsKGGfyv1jBp8ofUsP4qQRb1HX/BFLPtvFsSNep9jZO66NOPtD/5cQ/wIIfgFF+yGuzrt7RzG6SYis6iYvAlCu6D+gkilTrdO7x6K8pJgTxiWX7NApDib03ZOGoo2dBp/W9CgeR4/icUmVVRDkx52cqGo9MD3ShFAD9LLVoy5MhauSCX06djiIh1OKGdTjGJY7ajoNPlFEBIonAMklDC4IuoFubH8HZOge2N2WDcyM/qOqNwC3Auts8cRgSJCurptu9QRksIk86E6qqr9r9/lF4EWb3DEYOqcb6MYEIEPWEWyd2XG2qk4WkQba9msUQFW1pz2eGQzx6Q66MQEoS/hCQYKqqGqnXamjhDSIP9RIkaMcRyc5sAoNCdmT1ldVJ0f+VtrigMGQBl1dN2kFIBG5DZhOeAjzGuBiVd2VAb8KFk/Qz00LX2Pm+iVc4RzF71++i5smnMIxg8bEtavd9V8Wbb+PQMiD01HM+D7fYv8+30o4eOU1NjYliEjfeN+r6o5432cDoxtDQnQD3aTbCeENYH9VPQBYCVyfyk6ynZ8qlzY/+/gFZq5fijcURFXZ1FTPDz54mqU7rdKGhFlV/z/mb/snvlADIfz4Q40s3v4gy3Y9kRd1Sscmio0vU+cTnt9+vsUyL45dNrFNN4WgtUAg0OXqZHRjTVpPQKr6equPHxHOE5QUuchPlSubOo+bt79YhS/U9rbFG/Rz7/L3ufMI68OzeMcDkUSMXxHEyyfbH2K/PvFHdefjceiAjd1JVXWUPSXHxi7dFIrWfD4ftbW1XapORjfWZPId0CXAk7G+FJHLgMsAqqqqqKmpASAQCODz+XA6nQSDQerr63G5OnerMzu3291SRqI2qZTTmuagn6ucowk6wldNlZRwtSvc9FZaRwd/oji8X6MiRl6rmo3WNon6Z8dxaI8AjmBm27JF5FzCo7HHAZNUtdO7MhHpA4wFSqPrVNU63XLuyJluErGx83qJ2gQCAerq6rJ6XXZX3aRCNnXT6ZEQkTeBjqmf4VeRGfIQkV8BAeCxWPtR1RnADICJEydqdXU1kL07mJqaGqJlJGqTSjmt2e3z8LMX/4Y3FM68fbVrDH8PrMGJcM6wg7h0YrWl3Qvr72eXb02H9T0cVVSPuTwt/+w4Dh3ITlv2UuAs4N5ENhaRSwmP4h4KLAIOBz4kPNFWxslH3SRiY+f1ErWpq6tjwIABtj4BdWHdJEW2ddNpAFLV2LOvASJyETANmKqqSYfrXOSnypVNz+JSvrXnoTy2Zh7NwfAcJQKUuoq4bJ8jY9pN7H8V72z+BUH9Ko29kxIOHfhD2+uUjk1rMt2UoKq1QDKdNK4mPIvjR6o6RUT2Af6QWa++Ih91Uyhaq6+vT/iHulDqlC+6SYGs6ibdXnAnEU7VfayqNqW6n1zkp8qVzXUHHMewit7ct+JDnD4H1YPHct0BUxleETslz5Dyw5i6x20s2HYP9f71VBYN4eB+lzG0PHbQStW/XNsAkTu5Dr+x/UWkdbPZjMjdfrbwqKpHRBCRElVdLiK2TDdrp24KQWsulyspu0KoUwZ1k2uyqpt03wHdBZQAb0TuRD9S1fhtRl0cEeEbYybwjTETqKmp4YrJ1QnZDe4xkVOH35dd52xCVK3GM2xT1Ylx7RJoxkqCjSLSG3ie8PW6E1if5D4yhdGNoVNi6Cb1/aXw3pQs6ybdXnDJZek0dFtSaUrorBkryX2dGfn3BhF5h3BCxVcztf8kfTG6MSREhpvgknpvCtnXjcmEYMg+ChKwvSkh2ptnGNAQWfYHFtjqlMEQiwzrJoX3pkS2z5puTAAy5IRMpxQRkTOBO4EBwEsiskhVvxZn+98DFwFrgeh9pZKlXnAGQyawKxVPS/lZ1o0JQIbsk4UBdar6HPBcEibnAWNU1ZdZTwyGLGGtm7iddzL83hSyrBsTgAxZRwBHwPb+pEuB3oRneDQY8p4YuonbeSeT700jZFU3JgAZso8q2NyUANwCLBSRpUDLgCtVPc0+lwyGOHQD3eRFAHK73SkN1ErFLpc20aSKXalOqQ+os11IDxGezXEJX7VlFzS5PIf5rJtCqFM+6CbZ96YRsqob2wNQV0wmmEpSxUKpU8opRezvBdekqnfY7USmMMlIC6tO+aCbFN6bQpZ1k+50DGnjdrsRESorKxGRhFOWp2KXaxun09nl6pTseYoioVCbxQZmicgtInKEiBwSXexwJBPk8hzms24KpU5GN9bY/gRUUVGBqtLQ0ICqJpXLKVm7XNsEg8EuV6dkzxNERnTbn1Lk4Mjfw1utK9hu2Lk8h/msm0Kpk9GNNXkRgLpaMsGoTTJJFQulTqm2ZWPP3RsAIuIEZqrqX21zIsPk8hzms24KpU5GN9bYHoCgayYTrKhIPqliOmXlqw0Qacu2T0iqGhSRC4AuE4DAJCNNp5xU7YxuMkteBCBDF0fV1ju5CO+LyF2EJ39rjK5UVZOKx5CfdAPdmABkyAl50JZ9UOTvTa3WFew7IEP3oKvrxgQgQ/ZRhYC9Uzuq6hRbHTAYkqUb6Mb2btiGbkIo1HbJMSLSS0RuF5F5keUvItIr544YDMnQxXVjApAh+6hCMNh2yT0PEE4lf15k2Q08aIcjBkNCdAPdZCQAichPRURFpH8m9mfoYijhpoTWS+4Zo6q/U9W1keVGYLQdjkQxujHEpRvoJu0AJCLDgBOBDem7Y+iaqO1NCUCziEyOfhCRo4BmOxyJlG90Y+iErq+bTHRC+CtwHZDKXBNA10wmaJKRtkKxq/mgNZcDD0farwXYQXiiLbuwRTeFoDWTjDRCN9BNWgFIRE4HNqnq4mSneY3SFZMJmmSk7VBFA4HEts0SqroYOFBEekY+77bLF7t0UyhaM8lII3QD3XQagOLNsAf8knAzQqeIyGXAZQBVVVXU1NQAEAgE8Pl8OJ1OgsEg9fX1uFydx8XO7Nxud0sZidqkUk48m0AgQF1dXcbqlIqNncehDUF7B9SJSAlwNjAScEV/+FX1pjhm6ZSXd7pJxKbQdNMVfz/a0MV10+mRiDXDnoiMB0YB0bu4ocACEZmkqlss9jMDmAEwceJEra6uBrJ3B1NTU0O0jERtUiknnk1dXR0DBgyw9a7MzuPQQrQ3j728ANQD82k1sVa2yEfdJGJTaLrpir8fLXQD3aTcBKeqS4CB0c8isg6YqKrbktlPV0wmGLUxyUgjqKJ+e5sSgKGqepLdTtipm0LRmklGGqEb6CYvMiGkmqwvV4kBU7UxyUjDKKD238l9ICLjIwGgS5DLc5jPuimEOhndWJOxAKSqIzO1L0MXQzUfhDQZuEhEPiPclCCAquoBdjpldGOISTfQjajmPtmdiNQB67NcTH8gqWYN40NajFDVAVZfiMirET9asy2XTWIiMsJqvapm+zrMGEY3Xc6HmJqB7qEbWwJQLhCReao60fhgvw+GwiEfrhfjQ/fB5IIzGAwGgy2YAGQwGAwGW+jKAWiG3Q5gfDAUHvlwvRgfugld9h2QwWAwGPKbrvwEZDAYDIY8pssEIBG5TUSWi8gnIvKciPSOsd06EVkiIotEZF6Gyj5JRFaIyGoR+YXF9yUi8mTk+zkiMjIT5bba/zAReUdElonIpyJytcU21SJSH6n3IhH5bSZ9MBQmRjdGN7aiql1iIZzc0RX5/1bg1hjbrQP6Z7BcJ7CG8CRNxcBiYN9221wJ3BP5/3zgyQzXfTBwSOT/SmClhQ/VwP/sPk9mya/F6Mboxs6lyzwBqerrqhpNnPQR4SSPuWASsFrDswX6gCeA09ttczrwUOT/p4GpkmoefgtUdbOqLoj83wDUAkMytX9D18XoxujGTrpMAGrHJcArMb5T4HURmR9JdZ8uQ4DPW33eSMeLuGWbiNjrgX4ZKLsDkWaKg4E5Fl8fISKLReQVEdkvG+UbChqjG6ObnJIXyUgTJd4cK6r6QmSbXwEB4LEYu5msqptEZCDwhogsV9X3suNxbhGRCuAZ4MfaceKoBYRTf7hF5BTgeWBsjl002IDRTXyMbuyjoAKQxphjJYqIXARMA6ZqpAHXYh+bIn+3ishzhJsC0hHSJmBYq89DI+usttkoIi6gF7A9jTI7ICJFhEX0mKo+2/771sJS1ZdF5B8i0l+TnAbAUHgY3cTG6MZeukwTnIicBFwHnKaqTTG2KReRyuj/hF/ALk2z6LnAWBEZJSLFhF+Wzmy3zUzgO5H/zwHejiX0VIi0i98P1Krq7TG2GRRtPxeRSYTPfUbFbCg8jG6MbuykoJ6AOuEuoIRw8wDAR6p6uYjsAdynqqcAVcBzke9dwH9U9dV0ClXVgIhcBbxGuGfPA6r6qYjcBMxT1ZmEL/JHRGQ1sIOw2DLJUcC3gCUisiiy7pfA8IiP9xAW8BUiEgCagfMzKWZDwWJ0Y3RjGyYTgsFgMBhsocs0wRkMBoOhsDAByGAwGAy2YAKQwWAwGGzBBCCDwWAw2IIJQAaDwWCwBROADAaDwWALJgAZDAaDwRZMADIYDAaDLZgAZDAYDAZbMAHIYDAYDLZgApDBYDAYbMGWZKT9+/fXkSNHZrWMxsZGysvLs1qG8eEr5s+fv01VB1h997Up5bp9R7Dt9p94X1PVk7LuWBfC6KZr+RBPM9A9dGNLABo5ciTz5s3Lahk1NTVUV1dntQzjw1eIyPpY323bEWTOa21nei4avKZ/1p3qYhjddC0f4mkGuoduMhaARMQJzAM2qeq0TO3XUPgoil+DnW+YBUTkrAQ286jqy1l3xgKjG0Ms7NRNImRCW5l8AroaqAV6Jmvodrtxu91UVFRQUVGRVbtc2gQCgRa7RMn3OqVynhTFq4GEt88w/wJeACTONscAtgQgcqybQtBasrophDoVoG4SIW1tZSQAichQ4FTgZuCaZGzdbje1tbWICKrKuHHjEjpJqdjl2sbn81FbW9ul6pTseQJQwE8ooW2zwCuqekm8DUTk0Vw5067cnOqmULSWjG4KpU4FqJtESFtbmXoC+hvhaX0r4zhyGXAZQFVVFTU1NQAEAgF8Ph9Op5NgMEh9fT0uV+dudWbndrtbykjUJpVy4tkEAgHq6uoyVqdUbOw8DlEU8Ns08aGqXpiJbbLE38ihbhKxKTTddMXfjyh26iYRMqGttAOQiEwDtqrqfBGpjuPIDGAGwMSJEzX6gi9bdzBWLxFzfddTV1fHgAEDbL0rs/M4RAmp4slDIYnIIFXdYlPZOddNIjaFppuu+PsRJVXdiMhJwN8JT3V+n6r+sd33FwG3AZsiq+5S1fuSLii+DwlpKxNPQEcBp4nIKUAp0FNEHk30rrKiooJx48Yl3Uaail2uberr65O64AqhTqm1ZQt+jddMbBv3E24Cs4Oc66ZQtJaMbgqlTrnSTaRTy93ACcBGYK6IzFTVZe02fVJVr0pq58mRkLbSDkCqej1wPUDkTu7aZJs0kj0x6djl0sblcuW1f7mwgXBTgi8Pxzyrql3BxzbdFILWktVNIdQph7qZBKxW1bUAIvIEcDrQPgBllUS1Zcs4IEP3IoTgUXsvNREZbrVeVTfk2heDIRFi6Ka/iLQeDDYj0kwbZQjweavPG4HDLHZ/togcA6wEfqKqn1tskxDpaCujvwqqWgPUZHKfhsIn/DLV9ieglyKuCOEmr1HACmA/O50CoxuDNTF0s01VJ6a56xeBx1XVKyLfBx4CjktjfylryzwBGbJOuC3b3ktNVce3/iwihwBX2uSOwdApKepmEzCs1eehfNXZILxf1e2tPt4H/CklB7/aX8rasv221ND1UQ03JbReOkNEhonIOyKyTEQ+FZGrM+uTLsC6acJgyAtS0Q0wFxgrIqNEpBg4H5jZegMRGdzq42mEB0JnjGS0ZZ6ADFlHEXzJ38kFgJ+q6gIRqQTmi8gbFr15EkJEWg/0dACHAF+ksi+DIRekohtVDYjIVcBrhLthP6Cqn4rITcA8VZ0J/EhETiOssR3ARen4mY62TAAyZJ1wW7YzORvVzcDmyP8NIlJL+AVrqr15Wg/2DBBut34mxX0ZDFknFd0ARHKvvdxu3W9b/d/SAzNDpKwtE4AMWSfdd0AiMhI4GJiTsg+qN6bsgMFgA/nw7jQR0tFWXrwDcrvdbNmyBbfbnXW7XNpEkyomQ77XKZXzFO5OWtRmIdKdtNVymZWtiFQQvpv6saruTqrgTohVZqGQy3OYz7ophDplUDcFQaLasj28dsVkgiYZaVtUxaopodPupCJSRDj4PKaqzyZUWHLkZXqGRDDJSAurThnUTaGQkLZsfwJyu92ICJWVlYhIwncJqdjl2sbpdHa5OiV7niAyoltdbZbOEBEhnM6jVlVvT7iwJFDVe7Ox31yQy3OYz7oplDrlSjf5QqLasr1GFRUVqCoNDQ2oalK5nJK1y7VNMBjscnVK9jxBuCnBG0q6+eAo4FvAEhFZFFn3y3QmjhORUwkPjiuNrlPVm1Ldn53k8hzms24KpU451I0tpKqtvAhAXS2ZYCpJFQulTqknVUy6F9xsMthEJiL3AD2AKYQH350DfJyp/eeaXJ7DfNZNodQpV7qxg3S0ZXsAgq6ZTDCVpIrplJWvNpA3bdlHquoBIvKJqt4oIn8BXrHbqXTI5TnMZ90UQp0KWDeJkLK28iIAGbo2iuCxvymhOfK3SUT2ALYDg+NsbzDYSp7oJhFS1pYJQIask+qAugzzPxHpTXgirgWE3fqXrR4ZDHHIE90kQsraMgHIkHUUIWCzkFT195F/nxGR/wGlqlpvp08GQzzyQTeJkI62bO+Gbej6hFTwBl1tllwRyczbBlX1thaI1TYGg93YqZtEyIS28qtGOSAYCOJwOggPM0nQJhgCwOlMPF5rCnO5d1Vsng/owciMo/FO+P2EU/0YDHlDnsyjFY+0tZV2ABKRYcDDQBXhYzZDVf+e7n4zzbKPVnLHVQ/w2ZLPKS4t4qRLpvC9P36D4pLYL/m2b2vgb39+hY/nrEFVmXDoaH5y7SkMrOoZ08YfCHLHzNk8PfsTvn1wFXfe/DC//PpxTNhzaDaqVSDY2pTQC5hPfJHU5ciXFgpFNwY7yfsmuLS1lYknoLTT5rvd7pT6ySdq9/mKL/jFSbfgafIC4G328cr9b7P9ix389smfWFcqEORHVzzEtm27CQbDTzPz567lh1c8yCOP/4DiEutDd8Njr/PWotV4/AEUWLN5Oz+4+zkevvZ89hoyIGN1KhQb+KopwQ5UdaQtBXeOLbrJttYyYRPNBdeV6lRoukmETGgr7dqlmzY/F7mcnv7r//B5/W3W+Tx+5ry8kK2fb2fgsH4dbD58fxW7G5pbgg9AKKQ0ur3Menc5U0/cv4PNjoYm3li4Cl8g2LasQJAH35jLLRedkrE6FYJNFAUC+d2UkHPs0E2h5E0zueDCdAfdZDS8xkubH8mOehlAVVUVNTU1AAQCAXw+H06nk2AwSH19PS5X5251Zud2u1vK6H9wBWfvNbXDPhxOB/MXz6VsTWmH73bsaOS0s4Zh9Spn1+511NRs67C+2efnuxMHE4wYDSgv4rLDhgBQWtTc4k+qdUrFpvVxyGY5cVEhEMrrpgRbyZVuErGx83qJ2gQCAerq6rJ6XRrd5AcZC0Cdpc1X1RnADICJEydqdXU1kL07mJqaGqJl3Pnsg7x8/9sE/W2fTIpKXDy88u/0G9ynw/4/+mAVD854nuYmX5v1pWVF/OwX0zm2elwHm/pGD7/71YyWJ6DLDhvCjDmbcDqEUyeN45KIP6nWKd3jkM1y4hECfF1cSKmSS90kYmPn9RK1qaurY8CAAbY+ARnd5IaMBKB00ubnIpfTOT85lTcfm0VzqwBUUlbMMeccZhl8AA49bAz9+1ey+YudBAJf9YLr06ecIyfvZWnTq7yUM47Yn5kffYrHH2hZX+xycckJh2a0ToVgE0URAiF7mxIi2bW/CYxW1ZtEZDgwSFVtyweXa90USt40kwsuTD7oJhHS0VYmesGlnTY/27mcBo8ayN9qbuCf1z7Cpx+spEfPMk6/8kQu+PnpMW2cTgd//8d3uPfuN3mvphZVOPrYvfn+D46nqCj2XcnPz61mUJ9KHntnAQIcMmYIPzvnWEZW9c1onQrFBgDNi7bsfxC+qTwOuAloIPzj3/mdQRawSzeFkDfN5IKLkB+6SYSUtZWJJ6CMp83PBqPGD+dPr/0qKZuePcv42fXT+dn10xO2cTocXHLioVxy4qHU1NRwxbeqk/S06xFC8AVtb0o4TFUPEZGFAKq6U0SKbfSnIHRjsI880U0ipKytTPSCy2jafEPXJGj/nZxfRJyEOxchIgMI37XZgtGNIRHyQDeJkLK2CqJ2hsJGFYIhR5vFBu4AngMGisjNwGzgD3Y4YjAkQp7oJhFS1lb+jnIydCHEdvGo6mMiMh+YSvjJ4wxVrbXVKYMhLvbrJhHS0ZYJQIasowr+/BDSl8Aswtd9mYgcoqoLbPbJYLAkj3STCClpqyADUDAQZM7LC9lQu5Fhew/h8GmH4HR1/rJuRe0XLFqwnp69yjimeh/KKzoOQG3P5h27eWfxalRhyoFj2KNfr05t3D4fr65ZiTY1sWDzFxw8aHCnyU+DGuLDbStY6/6S4T36M3nAOFyOgngB2SmaB3dyIvJ74CJgDZG26sjf4+zyqRAIhZQFC9axauUWqgb1YvLkvSgu7vxnY+WObdRs+IwyVxEnjx5L/x7lndo0B3awwf02nqCDXd7P6F0yqlMbTyDAa5+tYlPDbg4YOIgjhwzHkUSi4XwmH3STCOloq+AC0K663fx48m/Y8eUuvE0+SnoU03tAL/7+/u/pM9A6OIRCys2/e46PP1yN3x+kqMjJP+94g1v+cgH7jY+dJPS/sz7hL8/UtGRDuHPmbK6afhQXTp0Q0+aTL7dw4fP/JajKFf2quOH5pzli6DDuOfV0XA7ri6ne38T3P/4ndZ56PEE/pc5iehb14L7DrqRfSWXiByePCYVs/1E4Dxijqr5OtzQA0Nzs46fX/IcN67fh8wUoLi7iH3e9wR13fps9hliPn1NVfv/+O/xn2ScENYRLHNz8QQ13njCNE0btGbOsDQ01zP7ydwCUBr7By59fxF69zmLigKtj2qyr38nZzz2OJ+DHEwhQ4nKxV5/+PH7aeZQVFcRMop2SB7pJhJS1lRfh1e12s2XLFtxud6fb/vMn/2bL+jqaGzyEgiGaGzxs3bCNu370QEybd978lI8/XI3H4ycYDOHx+Glu8nHDL//bMtVCezbv2M1fnqnB6w/iC4QXrz/IXS++z/qtOy1tVJXLX55Jg89Hk9+PKjQF/HywcQNP1y6N6d8dK15iU9N2moI+QihNQS91nl3cuuy5To8HJHf8cm0D4aaEQNDRZrGBpUBvOwrOFtk+h/957AM+W1tHc7OfYFBpbvaxa1czf/jDzJg2H33xOY/XLsETDOAPhWgOBvAEA/zozf/R6Lf+ffKHGpn95Q0E1UtQvYASVC8r65/jy+aFMcu6+s2X2OlpptHvJ6hKk99P7fat3L2gQ0YjS4xuMkbK2rL9CSjZVBWznv24Q0qdYCDIB8/PRVUtm7pe+d8iPB5/h/Vej5+Vyzczbr8hHb6LNru1JxBU3lq0iktOnNThuxXbt1Hv8XRY3xwI8OSnSzh/vwMs6/TOl0sIaNtAGET5oG45IQ3hkNgXXmEkVRRC9jcl3AIsFJGlgDe6UlVPs8+l1MnFOXz9taX4fIE261SVVSu3sHt3Mz17lnWweW7lMjyBjlpzIMz6fB0nje6YReSLxjmIOL5qvIkQVC9rd79KVVnH6WR2eppZtm0roXYi9QaDPLPyU649bHLMeoHRTYZJWVt5EYBEhMrKShoaGhJIw2490VvcCeBifKexv4q5HhSN0cM9pBpzZEdc92LVKcb61iR//HJn06oihNT2poSHgFuBJdg4/idT5OIcxtNUrO9Ura9aReNqSiy/UzTGqVJVJIbYEpkM0ugmo6SsLdvDa0VFBapKQ0MDqtrpyTli+sQOHQ4cTuHQkw+K+aL/xFMOpLS0Y5twUZGTvfcZbGlTfcAYrHbncjo47iDrtux9+g+gsrjjAOBSl4tz9t3P0gbg2AH74Wz3lCPAob3HxH36geSPXy5tWhMKSpvFBppU9Q5VfUdV340udjiSCXJxDqcev1+HtFMiMHJUf3r16mFpc/pe+9LD1VFrIYXJw0ZY2gzucRghgh3WO6WU0ZUnWdr0LevB2L79OoSgYoeTM/fa19KmNUY3GSVlbeVFABo3bhxDhw5N6PH0B3+/mH579KEs0oOttLyEvoN6c/U/vhfTZuoJ+3PgISNaglBxsZPS0iJ+93/n4HRZH4Ih/Xtx1fSjKCly4nQ4cDqEkiInl37tMEYNss7r5hDh7pOnU15URKkzLNwyl4sJg/bg6/uOj+nf1ftMo6q0N2XOcPAqdRTRt7iC6w84J+6xgOSPXy5toqiChhxtFhuYJSK3iMgRInJIdLHDkUyQi3P4rW8fxbBhfSkrC+umpMRFRUUpv/r1GTFtjhoynNPHjqPM5cIBFDkclDid/Pm4k6gsLrG0KXZWcOTAX+OUEoRwWU5KGF15MlVlsU/R348/lV4lpZRFpjfo4SpiTJ++/OCQw+PWC4xuMkzK2rK9CQ6SS9bXd1BvHlz+N95/7mM21G5i2D5DmHzWpLhTaztdDn5/63l8smgDixeGu2FPmbofvXpb38VFuXDqBI4eP5q3Fq1CFY47cM+YwSfKhMFDmH3R93hx5QpKPt/EfUecyRFDh8Xtht2nuIInjrqG97YuY617C8PLB1I9cD9KnIn15Mn7pIoQs9kyh0RfJLT+dSrobtjZPoc9epRwz4xLmPPRGlau3ExVVS+qp4yjrCx2mi8R4ZbqE/nGfgfw9vq1lBcVc+qYvRlcEb8358ieJzCg7EDWu99k/YZijh12L/1K94lrs2effrz/rct4afUKNjbUc8DAwUwZPgpnjN6m7TG6yRgpaysvAlCyFJcUMeX8o5KyEREOPHgEBx5s3QwQixED+1h2OIhH79IyvnXAQdTs2MWRw4YnZFPkcDF10AFMxbqjQkGjQsjmHjyqOsVWBwoUp9PBkUeN5cijxiZlN37AIMYPGJSUTXnRQPbt8w22Oms6DT5f2RRz3rjYrQsFTR7oJhHS0VZBBiBDAWLTy1QRuVBVHxWRa6y+T3UqBIMhJ+RxJ4RMaCv/w6uhaxBqt+SO6BD8SosltXYRgyFXpKAbETlJRFaIyGoR+YXF9yUi8mTk+zmRKeFTIW1tmScgQ/ZRUJt68KjqvZF/31TV91t/JyLJteMaDLkkBd1EpkW4GzgB2AjMFZGZqrqs1WbfBXaq6p4icj7hLtRfT9q9DGjLPAEZcoKEpM1iA3cmuM5gyBtS0M0kYLWqro2kxnkCaD/18+mEx+4APA1Mlc6SVcYnZW1l5AlIRE4C/g44gftU9Y+Z2G8s5r66kPt/+R82rd7CkDFVXHzzNzjslPi9/j5fv50Zd73J4oXrqago4azzDuOs8w/D4Yh93N3NXu6e+T6vzl2BqnLihL246vTJ9CyPncRUVXl48ULuWzifCyp68e8XnuGXk49lr3794/o3Z9sq/rnqVdY3bmOPsr5cPvZEjh44Lv6BKBRUwKacViJyBHAkMKBdW3VPwterbeRaN6mw7MMVzLjuEdYsXk//Pfpy4W/OYeo3j45r8+XOBv7+3CxmL/mMkmIXZx41nktPnkRxUeyfG18owANr3mbmxrmc4R7N24uf4Kq9T2Zgafzkv6GmF6DxbghtBdc4pPLnSPFBqVQ1/0hNN0OAz1t93ggcFmsbVQ2ISD3QD9iWTEGZ0FbaT0CtHvlOBvYFLhCRzkeCpciHL87jxnP+zJpF6/C4PaxZvJ7fn/cXZj8XO//T1i/r+eGlDzDng1U0N/mo29rAv+97lztuiz37cTAU4rt/eYpnZy+hvtHD7iYvL3zwKd+57Qn8wY6D5qLcPOtd/vTBLDY17Cakyqz16zj7qf+wfteu2HWqW8F1Cx9h+e4vaA76WOPewq8XP86bmxcndEwKgmC7JXcUE26PdtG2jXo30PlAqyyRa92kQu2cVVx3wk18+v4KPG4PG1d+wV+/fy/P3xVbNw1NHi784394Y/5K3B4f23c38cib8/npvS/GLeu6hY/w+LpZ7PC5CWmIt7Z8wsUf3kVjoGNqqyihxodg928huA60Cfzz0R3fQX1dWjf9RWReq+UyG71LW1uZaIJL5JEvLskk65tx3SN4m9omNfQ2+Zhx3SMxbZ5+/CO8Xn+bVCBej5/XX/mEnTusy/xo2Xo2bavHH/jqzZ8/GKJul5v3PllraVPv8fDYkkU0B77Kn6WEU8bfM//jmP7dufIVvKG2+bO8IT93rnwlpk1r8j2pImpfE1xkVPaNwOGqemOr5XZVXZUzRzqSU92kYvPgr/9joTUv//7NkwQD1ncRL3z4KY3NPoKhr8Tm9QeYt3IjqzdZ32CvatjMoh2f4Q19pZsgSmPAy/82zbe0UfWD+w6gud03zaj7b53WDQpWN9tUdWKrZUY7q03AsFafh0bWWW4jIi6gF7A9Oecyo61MNMEl8sgXk2ST9W1atdly/eY1X8ZMRlq7dBOBQMcuJEVFTjas306fvh3LW7mpDq8/0GF9k9fPyo11TD2447iIz3btpMjpxNvuCSmoyqIt1n4DbGi0FuZWTz2BUDDuvECFkFQRiJXCL2eo6np7PehATnWTis2aRdaHzO/zs6tuN/0Gd5ySYcnaLXgsdON0CKs2bWPPIR2bolc3bLbUrSfkZ+muDXx9hMX77NB20I5JTwEIdD4ZZxfWzVxgrIiMIhxozge+0W6bmcB3gA8JP6m8rYkk0IvlYhraylkvuMij4mUAVVVV1NTUABAIBPD5fDidToLBIPX19bhcsd264K/TCPgsLvAiJ++++1X6Ibfb3VLG5OP7ccChPTpkBBUR6ratpqbmsw776x/0cunhQzpk2xWBwY7dLftuTSAU4vv9qlqSIVYVFXPNoPB8Q71KSy1tAC72jccfsqiTOJj93ixLm5YyOzl+rY9DojaplNMZqTz1FMI7kmyTKd0katP6ejn991PxNno77EccwifLFyErOp7TQwfCqCOHdEg8KgLs+pyami0dHQv6uKB57xat9dMyLvLujyD031JBzc4ai9ooBC7Hsm+y9ACnlc1XdFXdRN7pXAW8Rlg3D6jqpyJyEzBPVWcC9wOPiMhqYAfhIGULmQhAiTzyEXlUnAEwceJEra6uBpK/Q2hc9Qb/vP4hvE1fCaOo1MUlf7iA6D4BampqWj6vW1vHVZc+gLfVlAwul4ODJozg0u9PtSzHHwgy/TcPsH13Y0tzggj0Li/lfzefSVmxdZqcF155kTfXrsEbDHLNoKHcvmUjpU4XT517PvsPrLKu06b53LbsBTytmuGKcfLdkcdRvXe1pU2Uzo5f6+OQqE0q5cRFQZJ875Ngd9JCJqe6SdSm9fXycfNCbjr3z22a4YpKXUy/4kSm/MR68HtdvZuzbvg3ja205nQIowf15YkLj7d80lFVvvXhnaxzf0lAQ1zk3Z9/lyylh7OY/x59bcxJGUO7P0abHkP46j2RUoqjz11IyTFpHYtC1Q2Aqr4MvNxu3W9b/e8Bzk1+z5knEwEokUe+mEST9UXTlHd2ck659Hj8Hj8P3/hfmhqaKS0v4byfn8bZV0+PaTNy9ABu/vPX+dufXmbzpp04HA6OOW4ffvzzaTFtilxO/v2z87nhkdeYv3IjAONHDuKmi06KGXwA/nzCyfzfezUtE9ANrezJzcedEDP4AJw6ZAKeoJ8Zq9+gMeChxFHEN4cdxbf36jzDRbLHL5c2bUh+8GnLOxIAEYm+I0kqAInIncRpyFDVHyXtWWbIqW5SsZl08sFc86/LuffaR6jfthtXsYvpV57IpX/4ZkybAb0q+Nc153HTI2+walMdAhy57whuvOjkmPkQRYS7J36XP3z6HO/XhZvPxpYP4jcHnBt3RmCpvBakCG18CPCj0gdHz+s7DT6pHItc2rQhj3PBZUJbaQegWI98yewjmRMjIpzxw1M47Qcn0dzQTFllGY4Ekg8eeMhIHnziSpoavRSXuHC5Ou8lOKhvJfdcfQ4eXwBF4waeKCUuF78/7nh+e+wUZr33Hu+ee17cRKRRzh5+OGcOm0RTwEcPV3Gn0zC0Jt+TKopaNiX0F5F5rT7PaPdCNa13JK2IlnEU4d5mT0Y+n0uSwSyT5Fo3qdocd8HRTDl/Mk27myitKMXp7Fw3+wwbyH9++U2aPD5cTkfc7tdRehWXc+vBF+ILBXj/vVlcPjl2kIsi4kQqr0Errg73gpOKhLQWpUB1k0+kra2MvAOyeuTLNg6Hg/Je5Z1v2I4e5dYp4eNRWpz8YSpyOnGIJCUIhzioKIo9xqiQsWhK2KaqE7Ndrqo+BCAiVwCTVTUQ+XwPEP8FW/Z9y7luUkFEUtNaaeys2bEodrhiTjQXCxEnSPxs24VKKk1wuSIT2jKZEAzZR0FCbZcESOgdSRL0ITxALkpFZJ3BkJ+kphs7SFlbJhecISekIJ603pFY8EfC89a/Q3jC2WOAG9LYn8GQdfI46LQmZW2ZAGTIPpq8kDLxjqTd/h4UkVf46j3Sz1XVok+wwZAnpKAbO0hHW90qAG3esZuFqzbRq7yMw8YNx+XsvAWywetl1obwOKvJw0fQs6Tzd0hBDTF32xp2+5vZ1LSDIT3iz6Ia5cvmVWz3radP8RAGle6T0PsjDW4B38fg6AnFRyGS2CyqOScFIWXyHUkk2eLxwGhVvUlEhovIJFWNnaLCAMC6NVtZXbuZqsG92P+QEQldl9sbm/hg/QbKioo4etQIShIY++IP+VjtXkhz0I07sIsKV+8MeF/gFEAASkdb3SIAqSp/ffo9nqpZjMvpQEQoKXIx45pzGL1Hv5h2r61exU9eewWXw4GqElTl1uNPZPresWdr3NC4jcvn3EdTwMsFnr04b9bfmD50Aj/f97SYwvWHPDz/+W/Y0rwcwQEofUuGc/bwWylxxn75G9r9F2h6EFqCTgn0fQgp2juRw5IzhLy4k/sHYTkfB9wENADPAIfa6VQ+E/AHufnnTzHvg9U4nA4E6D+wJ7fddzF9+sXu1XX/x/O5fdb7Ya0hiMB9Z5/BhKFDYtqsb1zGI+tuBpQx/hP4y/LLOKHqWxw5IPbwiq5OnugmEVLWVl50Qkg1V1Kidu99spZn3luCLxCkyeun0eNjR0MTP7rr+ZasBe3Z1tTET157BU8ggNvno9HvxxMI8PM3XmdzQ4Oljary0/mPsM3bQGPQS4gQvlCAlzct4I3Nn8T0b/bWB9jcXEtAvfi1Gb96qPOs5Z0td8W0Ue970Pww4ANtBG1EdQe681I0gYnkc57TKth2sYHDVPUHEB61qKo7CSdTLFiyfQ6fefQD5n+4Bp83gKfJR3OTjy8+38Gffv1sTJslm7fw19kf4A0GafT5cft8NHh9XPrMC3hb5UhsjT/k45F1N+MNNeENNaMaIqB+3vzyUb5oXpPROqVr1w11kwgpa8v2ABQdKbxx40Zqa2sTPknJ2D317mKafR3zRu1saGLF53WWNq+sWmm5PqghXlq1wvK79Y11bPHsQtuNzWoO+nl6Q+xs3cvq3yCobZM+hgiwYve7MQOkNj0O2jYRowAaagD/kphlQWrHPNXz1OKb/b15/JHsCgogIgMoiAYOa3JxDl96el6b7CEAwWCIT+avo7HBOkv1f5d8is8iW3woFOL9dRssbVa7F2I1njGgfubveCuuj7n4/ci1TWvyQDeJkLK28iIAiQiVlZWISFIXUKJ2zR7rpIUiWAYmCGewDoY6HsNgKESTP4ZN0I8jxiFtCvos1wMEYyRVDBFEY53HUKPlakXCg/LikMoxT/U8RZzKByHdATwHDBSRm4HZwB9s8SQD5OIctg8+X6H4LRKOAjT6fB3yJ0K4daA5hm78IZ/lcHpF8YXaZ7tuSy5+P3Jt00J+6CYRUtaW7QGooqICVaWhoQFVTXjEcDJ2J07cO8ZgUmG/EdYpcqpHjsJpkWGh2OXiuFGjLW32rByEy8KmxOHixMEHxPRvRPkhkXc/rVAYVDIOh8QYeV56KlDWYbUQguKDY5YFqR3zVM9Ti182NyWo6mPAdcAtwGbgDFX9b+49yQy5OIdHTtkHl6vj9Vy1R296W2SQBzhp77H0KOrYESagyhEjhlvajK4YT0g7BrQiKWG/XkfG9TEXvx+5tmmN3bpJhHS0ZXsnhFRzJSVjd+bR+/PSnGWs3byDZq8fp0NwOZ3c8J0TY6YJGduvH98cfwD/WfIJnkjbdanLxVnj9o2Z183lcHLjAedy/cLHCWj4ail1FDGsvB/nDj88pn/VVVfyRfMyAiEvAfXipAino5gTh1wT00Z6nIl6ngX/cqAZxQm4kF43IxI/m0Kuc1qJhhc7EZH7gTtV9e5W625Q1Rvs8yp1cnEOv3PlccydvYrdu5vxNvtxFTlwuZxc939nx7SZuucYDhs2lDmfb6TJ78chQpHTwXXHHk3fHh1vmAAqXL05oepbvPnlowQirQFFUsLoigPYq3JCRuuUjl131E0ipKMt2wMQpJ4rKVG7kiIXD/zs67y1YBWzlnxGv149OGvyeEYOit89+lfHVHPCmD15fnktIVVO33sfDh86LK7N5IH78J/JP+L5z+fS67Nmrt//DI4fNJ5iZ+xD3at4EBePeZBPd73Gl55VDCgZzf69T6LMFXs6YpFi6PsoeF5Hve8gjn5Ij3MR157xD0aEXOa0grxoPvgaMFFE/qKqD0fWnUYBD0bN9jns3beCfz17FW+9tJglC9YzZEQ/TjlrIv0H9oxp4xDh3rNP5+3Va3lt5Soqios554D92a9qYNyyjhwwnZEV+7Jgx9sEtlRy3vCfslflhIRyImb798MOmyh5oJtESFlbeRGAckGRy8lJk/bhpEmxu1BbMWnIUCYNGZqUzbDyfvxwn5Oo2VJD9ZD4zWFRSp2VTOiX3AzRIkVQdipSdmpSdjknxbTyGWYrMAV4VEQOA66GJJOOdUPKepQw7dxJTDt3UsI2DhGOHzuG48eOSaqsPcrGsMeQMdSsqmGfnqZ3fJ7oJhFS1pbt74AM3YM8eJkqqlqvqtOBOqCG8FTEBkPekge6SYSUtWUCkCEn5IGQZkb/ibRN3wqss8UTgyFB8kA3iZCytrpNE5zBPiQPmhJU9XftPr8IvGiTOwZDp+SDbhIhHW2ZAGTICRKypzuPiMxW1cki0kDb0Y4CqKrGfqNuMNiMXbpJhExoK60AJCK3AdMBH7AGuFhVd6Wzz3zCGwjw8KJF/PfTpajCWfvuyyWHHNJpYsW3tyzj4bWzmNzQm/eXvMCle1ZTVRa/SfRLzxo+qPsPWz1r6FcynCP7f4M9esTvMKHBbYQa/4F6a8DRG0ePi5HSaUlNgpcTbMzqq6qTI3/zZsayrq4bQ4bI82zYmdBWuk9AbwDXR1Ln3wpcD/w8zX3mBarKJc8/x8LNm1vGAd055yPeWruGp75+Po4YP/L/XvMeM1a/gyfo5zCt5PmN83lzy1KeOvqHDCi1viHY1FTLk+uvJ6A+QKn3f8mGxk84a9hvGVlxiLV/oV0Et58GoZ2AH4IbCNX/EvEvx9nzZ5k4BBnFLiGJSNy+9qq6I1e+tKLL6saQWfI5AGVCW2l1QlDV16PTsAIfEZ61MmnyMZng3E2bWLxlS0vwgXB6nuXbtjFr/TpLm6aAl3tXhYNPlKCGaAr4eGht7Blq395yDwH10vopNqBe3tzyz5g2oaZHIFQPtE5v0ow2/RsNdf6b2o2SKs4nPHf9fItlXhy7rGGnbvJRa+1tAoFAl6tTAeomEdLWVibfAV0CPJmsUTRZn4igqowbNy6hQVup2CVjs2DzF5bZe5v8fhZ88QXHjhzV4bu17jpcDgfednctfg0yd/vamH596bX+bodvIyENWqbjUe/7gLfjeopQ/zKkZHLM8rJ97NoTTitvT1u2qnY8UflFznSTr1prb+Pz+aitre1SdSo03SRCJrTVaQASkTeBQRZf/UpVX4hs8ysgADwWZz+XAZcBVFVVUVNTA0AgEMDn8+F0OgkGg9TX1+NKYPKqzuzcbndLGYnatGaQx8NPhgzpkFhREAbX7+6wb4CABvmmZ8+WbNj9tZTvBfYFoKKh1NIGYIT3LMtcWIKD97bGeHIKTgPtGGQUB+JqJNwVP/3jkI5NK6fyoilBRPoAY4GWXEWq+l6Wyso73SRiY+f1ErUJBALU1dVl9brM9u9Huv4BeaObREhVW50eCVU9vpOCLwKmAVM11twB4f3MAGYATJw4Uaurq4Hs3cHU1NQQLSNRm9Y0+/0c+a8Z1HvbPmVUFBcz+/TTY86M+sycB1mwYx1+DfK9wL78y7WMUoeLuyZdxCF9R1razN2+g1lbH440w4VxaBEH9zqN6qHVljbq/5Tg9q8TmYIDgJA6wbkPxQNfyNhxSMemNQ7r5Mk5Q0QuJTxCeyiwCDgc+JDwJFoZJx91k4iNnddL1Kauro4BAwbY+gRkdJM46Wgr3V5wJxHOgnqsaidzAMQgX5MJlhUV8cR5X+eq/73Ixt27AagqL+fu6afFnZb7T4dcwK8XPcVH29YgCBXOEq7bb1rM4AMwse+ZNAXqmb/jeQQHIQ0yrnIqU4ZcHNNGivbD0esvhHb/GlUPEERdB1Lc7+6YNqkch3RsWtC8aEq4mvAMjR+p6hQR2QebpmOwSzf5qrX2NvX19Qn/UBdKnQpYN4mQsrbSfQd0F1ACvBHp+vuRql6e7E7yNZng3v3788ZFF7Oxvh4Fhvbs2WkX58qiUv5+6LfZ6W1k7vsf8lb1eRQ5YkypEEFEOLbqYo4YcAEN/q1UuPrFnYo7iqPsa0jp8RBcD9ITcfZPqF6Q26SKeTK1sEdVPSKCiJSo6nIRsWvuctt0k69aa23jcrmSsiuEOhWwbhIhZW2lFYBUNbHUywXO0F7JpwzrU1JOscPVafBpTbGjlH4l1nOmxELECS7r+YnyBlUkaPud3EYR6Q08T/iHfyew3g5HuotuDGmSH7pJhJS1ZTIhGLKPYruQVPXMyL83iMg7hJMlvmqjSwZDfPJAN4mQjrZMMlJDTsiHpIoi0kdEDgAagI3A/vZ4YjAkRj7oJhFS1ZZ5AjJknzy4kxOR3wMXAWuBqJSVLPWCMxjSJg90kwjpaMsEIEPWEfJCSOcBY1TVZ7cjBkMi5IluEiFlbZkmOEP2iXQnbb3YwFKgtx0FGwwpkR+6SYSUtZUXT0ButzulfvKp2OXSJprTqivVKaXxDCgStL0B+xZgoYgspVUOI1U9zT6X0iOX5zCfdVMIdSpg3SRCytqyPQB1xVxOqeS0KpQ6pTSiOz/ash8iPFPjEr5qpy5YTC64wqpTAesmEVLWlu1NcG63GxGhsrISEUk4Y2wqdrm2cTqdXa5OyZ6nKHnQlNCkqneo6juq+m50scORTJDLc5jPuimUOuWDbkSkr4i8ISKrIn/7xNguKCKLIstMq23akbK2bH8CqqioQFVpaGhAVZNKpZGsXa5tgsFgl6tTsucJQFSRgO0PHbNE5BbC89e3biZYYJ9LqZPLc5jPuimUOuWJbn4BvKWqfxSRX0Q+W81D1ayqByWx35S1lRcBqKvlcoraJJPTqlDqlFpbNhCyPQAdHPl7eKt1BdsNO5fnMJ91Uyh1yhPdnA5UR/5/iHDK/ExMhJiytmwPQNA1czlVVCSf0yqdsvLVBrC9LVtEnMBMVf2rbU5kgVyew3zWTSHUKYO66S8i81p9nhHJmJ4IVaq6OfL/FqAqxnalkTICwB9V9flYO0xXW3kRgAxdHFUI2jedo6oGReQCoEsFIEMXx1o321R1YiyTePNQtd21qojEuiscoaqbRGQ08LaILFHVNdYupqctE4AMOSEPevO8LyJ3EZ59tDG6slDfARm6B8nqJt48VCLypYgMVtXNIjIY2BpjH5sif9eKSA3hJjbLABQhZW2ZAGTIPgpkeDyDiNwGTAd8hMVxsaruimNyUOTvTe08K8h3QIZuQOZ1MxP4DvDHyN8X2m8Q6RnXpKpeEekPHAX8qZP9HhT5m7S2TAAy5ADNRieEN4DrVTUgIrcC1xPnhaqqTsm0AwZDdsm4bv4IPCUi3yU8XcJ5ACIyEbhcVS8FxgH3ikiI8DCdP6rqsrhepqEtE4AM2UcVApmdW1hVX2/18SPgnHjbi0gv4HfAMZFV7wI3qWp9Rh0zGDJFhnWjqtuBqRbr5wGXRv7/ABifzH7T0VZGBqKKyE9FRCOPbAZDW6JNCa2XSG+eVstlaZRwCfBKJ9s8QDhV/HmRZTfwYBplpo3RjSEu1rrJR1LWVtpPQCIyDDgR2JDuvgxdmI5NCXF780D8Hj2q+kJkm18R7i76WCcejFHVs1t9vlFEFnVikzWMbgwJYf/4uURIWVuZaIL7K3AdFi+0EqUrJhM0yUhbkWJTQrwePQAichEwDZiqqp11F2oWkcmqOjtiexTQnLRTmcMW3RSC1kwy0ghZaLrOEilrK60AJCKnA5tUdbGIpLSPrphM0CQjbY9moxfcSYR/wI9V1aYETC4HHo60Vwuwg/AkWjnHLt0UitZMMtIomddNlkhZW50GoE4GNv2ScDNCp0Ta+C8DqKqqoqamBoBAIIDP58PpdBIMBqmvr8fl6jwudmbndrtbykjUJpVy4tkEAgHq6uoyVqdUbOw8Di0oaOYHot4FlABvRH7EP1LVy2O6oLoYOFBEekY+7860Q63JR90kYlNouumKvx8tZEc3GScdbXV6JGI1g4jIeGAUEL2LGwosEJFJqrrFYj8zgBkAEydO1OrqaiB7dzA1NTVEy0jUJpVy4tnU1dUxYMAAW+/K7DwOLWSnF9yeyWwvIiXA2cBIwBV98lDVm+KYpUw+6iYRm0LTTVf8/WihQJrg0tFWyk1wqroEGNjKiXXARFXdlsx+umIywaiNSUYaRfPhTu4FoB6YT6uMvbnGTt0UitZMMtIoeaGbREhZW3kxDijVZH25SgyYqo1JRhpBsTUXXIShqnqS3U5kklyew3zWTSHUqYB1kwgpaytjAUhVR2ZqX4auhaoS8tvelPCBiIyPPIHkDUY3hljkiW4SIWVtSee9VzOPiNQRTgWRTfoDSTVrGB/SYoSqDrD6QkRejfjRmm25fCIRkWXAnsBnhJsJhHBS4ANy5UO6GN10OR9iagbyQzeJkI62bAlAuUBE5nU20NH40H0QkRFW61U12z/oBUU+XC/Gh8IiHW3lxTsggyHbmEBjMGSHdLSVkVxwBoPBYDAkS1cOQIlOU5tNjA+GQiMfrhfjQzehy74DMhgMBkN+05WfgAwGg8GQx3SZACQit4nIchH5RESeE5HeMbZbJyJLRGSRiMzLUNknicgKEVktIr+w+L5ERJ6MfD9HREZmotxW+x8mIu+IyDIR+VRErrbYplpE6iP1XiQiv82kD4bCxOjG6MZWVLVLLISTO7oi/98K3Bpju3VA/wyW6wTWAKOBYmAxsG+7ba4E7on8fz7wZIbrPhg4JPJ/JbDSwodq4H92nyez5NdidGN0Y+fSZZ6AVPV1VY0OG/6IcJLHXDAJWK2qa1XVBzwBnN5um9OBhyL/Pw1MlVTz8FugqptVdUHk/wagFhiSqf0bui5GN0Y3dtJlAlA74k3RrMDrIjJf0psGOsoQ4PNWnzfS8SJu2SYi9nqgXwbK7kCkmeJgYI7F10eIyGIReUVE9stG+YaCxujG6CanFNRA1HhzrGjiUzRPVtVNIjKQ8Fwyy1X1vex4nFtEpAJ4BvixdpyTYwHh1B9uETkFeB4Ym2MXDTZgdBMfoxv7KKgApBmYollVN0X+bhWR5wg3BaQjpE3AsFafh0bWWW2zUURcQC9gexpldkBEigiL6DFVfbb9962Fpaovi8g/RKS/JjkNgKHwMLqJjdGNvXSZJjj5aorm0zTGFM0iUi4ildH/Cb+AXZpm0XOBsSIySkSKCb8sndlum5nAdyL/nwO8HUvoqRBpF78fqFXV22NsMyjafi4ikwif+4yK2VB4GN0Y3dhJQT0BdYLlFM0isgdwn6qeAlQBz0W+dwH/UdVX0ylUVQMichXwGuGePQ+o6qcichMwT1VnEr7IHxGR1YTnSz8/nTItOAr4FrBERBZF1v0SGB7x8R7CAr5CRAJAM3B+JsVsKFiMboxubMNkQjAYDAaDLXSZJjiDwWAwFBYmABkMBoPBFkwAMhgMBoMtmABkMBgMBlswAchgMBgMtmACkMFgMBhswQQgg8FgMNiCCUAGg8FgsIX/B2aLKFbVBk0pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nsamples=len(T)\n",
    "random_samples = np.random.choice(nsamples, 4)\n",
    "\n",
    "def rectangular_array(n=9):\n",
    "    \"\"\" Return x,y coordinates for rectangular array with n^2 stations. \"\"\"\n",
    "    n0 = (n - 1) / 2\n",
    "    return (np.mgrid[0:n, 0:n].astype(float) - n0)\n",
    "\n",
    "\n",
    "for i,j in enumerate(random_samples):\n",
    "    plt.subplot(2,2,i+1)\n",
    "    footprint=T[j,...]\n",
    "    xd, yd = rectangular_array()\n",
    "    mask = footprint != 0\n",
    "    mask[5, 5] = True\n",
    "    marker_size = 50 * footprint[mask]\n",
    "    plot = plt.scatter(xd, yd, c='grey', s=10, alpha=0.3, label=\"silent\")\n",
    "    circles = plt.scatter(xd[mask], yd[mask], c=footprint[mask], \n",
    "                          cmap=\"viridis\", alpha=1, label=\"loud\")\n",
    "    cbar = plt.colorbar(circles)\n",
    "    cbar.set_label('normalized time [a.u.]')\n",
    "    plt.grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input 2: signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200000, 9, 9)\n"
     ]
    }
   ],
   "source": [
    "# signal map\n",
    "S = f['signal']\n",
    "S = np.log10(1 + S)\n",
    "S -= np.nanmin(S)\n",
    "S /= np.nanmax(S)\n",
    "S[np.isnan(S)] = 0\n",
    "\n",
    "print(S.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEfCAYAAAAHqhL5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABgoElEQVR4nO2dd5xcZfX/32dmtu+mQAoloYcmnaVJ24ggIE0RBBFB0AiKgoCo4A8RRYpYqGK+gCICAaUY6aAsTQOk0EIoCSSQSGQDyWZn2+zMnN8fz93N7Oz0OzP3zuzzfr3mtTN3nnOfc+/ez5x7n3IeUVUsFovFYik3Aa8dsFgsFsvoxAYgi8VisXiCDUAWi8Vi8QQbgCwWi8XiCTYAWSwWi8UTbACyWCwWiyfYAGSxWCwWT7AByGKxWCyeEPLaAYul1IjI2mxFgA9Vdety+GOxVDrF0pQNQJbRwBJV3TVTARFZUC5nLJYqoCiaEpuKx1LtiMgWqvqu2zIWi8VQLE3ZAGSxWCwWT/CkCW7ChAm62WablbSO7u5umpqaSlpHJfhQLj/mzZu3SlUnpvruUBFdlcoGHlPVQ0vqWBZEZKaqzvDSh2JhdVV9PlSirvLRlCcBaLPNNmPu3LklraO9vZ22traS1lEJPpTLDxFZlu67VcBcSWGjTCihS7nyB68dKBZWV9XnQ4XqKmdNFW0YtogERWSBiDxYrH1aqohAipcPUNV5XvuQCasrS0Z8qKt8NFXMJ6CzgUXAmHwNw+Ew4XCY5uZmmpubS2pXTptoNDpklyt+PqaC7QQIptgey7naoiAiTwEjOj1V9TPl9SQvyqqrStBitemq0HPuB1251VRRApCITAE+D1wGnJuPbTgcZtGiRYgIqsp2222X0z+hELty20QiERYtWlQVx+TGDvDFnRlwfsL7euBYIOqRL1kpt64qRYvVpCtXmgI/6MqVpor1BPQ74AKgJV0BEZkBzACYPHky7e3tAESjUSKRCMFgkFgsRmdnJ6FQdrey2YXD4aE6crUppJ5MNtFolI6OjqIdU6E2yeeiVOc8LYIfhJKqaeB5EXnRE2dy43eUUVeFXEulrCudjR905eV5GMIHunKrKdcBSESOAD5S1Xki0paunKrOBGYCtLa26mAHXqnuulJ1Epb7rqajo4OJEyd6fqeWfC5G6xOQiKyX8DEA7A6M9cidjHihq0KupVLWlc7GD7ry8jwMw2NdudVUMZ6A9gWOEpHDMY9gY0TkL6r61VyMm5ub2W677fJuAy3Ertw2nZ2deV1Qfj4mN3Zp26rLzzxMe7VgmgneA0731KP0lF1XlaLFatJVwZoCv+jKlaZcByBV/THwYwDnTu38XEUySN4n3oVdOW1CoZCv/SvXOQc8v1MDUNXNvfYhV7zSVSVosdp0VbCmwHNdudWUD34WLKMCHw4XBRCRDbz2wWIpGB/qKh9NFdVdVW1X1SOKuU9LFTDYVJD88ge3eO1ANqyuLCnxr65y1pTNhm0pDz64M0uFqn7eax8sloLxoa7y0ZQNQJbS44PhoomIyCRMxz4Aqvq+h+5YLIXhI10VqimfuG+pegpoKhCRQ0XkLRFZLCI/SvH9JiLylJOq5lVnxFim/R0lIu9gRuo8DSwFHsn7WCwWv+BxE5xbTdkAZCk9g3dqeXSWikgQuAE4DNgeOFFEtk8q9hPgHmdhrBOAG7N48nNgb+BtZ/TOQcCcfA7FYvENBeiqBLjSlA1AlvKQv1D2BBar6ruqGgFmAUcnlVHW5UgbC/w3yz4HVPVjICAiAVV9CmjN+RgsFr/hfQBypSlf9AFVQgJEmzTRhV1hE+Y2Bj5I+Lwc2CupzCXA4yLyXaAJ+GyWfa4RkWbgGeAOEfkI6M7bswrB79eF1VXhNoBfJqK60pTnAahSEiCO9qSJbuyAdHdmE0RkbsLnmU5qmVw5EfiTqv5aRPYBbheRHVQ1nqb80UAv8H3gJMxT06V51Fcx+P26sLoq3GYY3rdhudKU5+6Hw2FEhJaWFkSEcDhcMrty2wSDwao5Jjd2GdqqV6lqa8IrMfisAKYmfJ7ibEvkdOAeAFX9D2YUTtrFuFS1W1XjqhpV1dtU9Vqn+aDq8Pt1YXVVuM0QPugDcqspzwNQc3MzqkpXVxeqmld+p3ztym0Ti8Wq5pjc2AGFCOUlYJqIbC4itZhBBrOTyryP6fRERLbDBKCO5B1JDou55VKmkvD7dWF1VbjNMDwKQMXSlOdNcIUm4yvErtw21ZQ00Y1dIW3VqhoVkbOAxxzrW1V1oYhcCsxV1dnAecD/icj3MQMSTlXVEYtjAfuJSHLwSvYweYRdReP368LqqnCbIbztAyqKpjwPQFAZCRALsam2pIlu7Aq5M1PVh4GHk7ZdnPD+DUzW6Gwkj55LRSQ/7/yP368Lq6vCbYbwrg2rKJryRQCyjAI8bOxV1ae9q91iKSEe6apYmrIByFJ6/DFc1GKpLqpAVzYAWcqD58NdLJYqpMJ1VeHuWyoCHwwXTUZExovITt56YbG4oEBdZcuxmFDuWBFREckps0EhmrIByFIefLBuiYi0i8gYZx37+ZgRdL8pvycWS5HIU1c55lhERFqAs4EXsuzPlaZcByARmepkJH5DRBaKyNlu9+k79L+gi0GjeRh1A/2YScK51tMHutDUV0345wlorKquBb4I/FlV9yJ7+h5PGBW6Yg3QlXtxjUPsP8Ba0JUl8qmCKExXueRYBJNk9EqgL8v+XGmqGD8DUeA8Vd0ekxX1O6kiakWiMdCTgC2AXYBtQJfnYHgXZkL+G8Ak4Mkc6noVM/F/H1OffgdSTmmpUPwRgEIisiFwPOD3iafVqyt6gIMx2lgf+CoQy2yi/dDfBpFDQN+FvmkQe6LUjvqf1LqaICJzE14zEixS5VjcOHGXIrIbMFVVH8rBA1eacv0zoKofqup8530XsIikA8pGOBxm5cqV+aWhKNAuP5vfAw9gnmS6UV0GfDmLzXuYDDF9QBwIA1/A3O1l4vPAKswdYT+qfwL+moOP5TgPLu38s3TwpZiJrYtV9SUR2QJ4xxNPsuCVrspzXZwPPAcMAAOo3gf8LrNJ9DLQuRg9xczfyLGm1aDo/vnfBsikq0wprjLvUiQA/AYzyTsXXGmqqKPgRGQzYFeytBsm4u8EiM9g7tYMIjFUX0YkUy2volqDyLqmN1VBZDFps5RrFzC8SUGkh0jkOWrrji/yMfkqGWlZUdW/khDVVfVd4FjvPMqNcumqfNfFMyS27Ij0Eo0+SSiU4Tcv/h+Sm7NV44i+D7J10Y/LzzbDyF9X2XIstgA7AO1ifug2AGaLyFGqmpg4GHCvKUmduSR/xKTkfhq4TM0tTfL3M4AZAJMnT9591qxZAESjUSKRCMFgkFgsRm1tLaFQ9riYzS5Vuvb86/ovJjAknCOtBdkxg00P8BYQJxyeQnPzcsytyo5ATQa7BZgnpkGEeGwjAsENMtjkdkzJ56IU53z69OnzVDVlhG0dKzp375Hb5XHS2hQTEbmOYf/E4ajq90rtQ6GUU1eFXEuF1bWE4S0CQjw+gUBgk/Qm+j7oKhBdpysVCOxMtsfpUpyL4pyH7DbF1pWIhIC3MTkUV2ByLn5FVRemKd8OnJ8cfIqlqaI8AYlIDXAvcEcqkTgOzQRmArS2tmpbWxtQuruu9vZ2BuvI1Wak02uJx3dDEwYF9Ef+TmNjW3obIBJ5gGDwLzz77GXsv/9PiEZ/RF3dwZkPSDtQPYV4XAGhr28TlOdobk6b3DnnY0o+F548AXk7YW7EnVslUG5dFXItFVbX5sTje6DaAyiRyPrEYv+muXlKehP9kHjvzqBreOaFy9l/z58Q5cfUNR2U8ZgK869Mvy8F2gyjNDkWc6EomnIdgMQ8p90CLFLVvIe0+joBoowhEHiV3r57ifR/QqjmYJqasvcD19bOpKfnWGKxTvr7H6Ox8YDsByTHIWzHQPRf9PY2UFNzLM3N6xX/mAq0cWM3NFrHI1T1Nu9qLwwvdFW+62JTAoF36O19jN7eAWprD6G5eWJmE9mQQMMb9PfcRFwnEdFZNDTnko7M34lFXScjLUGOxaTtbWm2F0VTxXgC2hc4GXhNRF52tl3oHGRO+DoBojTS0HAyDQ35+dbY+DmCwfbcgs9QXTtQX78D9fX51eX3pJOAL/qARGQi8EPM/Iehs6yqn/HMqfR4oqvyXRdjaWg4Pj9dyQTqmn5CINhOQ3Nbif3zt80Q3k/odqUp1wFIVZ/DxGKLJT0+CEDAHcDdmCGHZwCnkGL9ID9gdWXJCe915UpT3rtvqX78Mwx7fVW9BRhQ1adV9TTAj08/Fkt2/KErV5qyAchSekqYs0pEjk/IFnBnll0OOH8/FJHPi8iuQPaONovFj/gjw4grTdls2NnQGHQvgHg/NO4AobHZbeIDsPoJGOiF1f+EcdNBRnmsz/PwE3JWHYyZrf2SiMx2FqEbLDMN+DGwr6quFpFJWXb7CxEZi5lkdx0wBvh+fp5ZPGP1u/D2/RDeCN56AKYdAYFR/hPm/c+KK02Nvv+eDgABkByeVT/6Eyz7gQk+EjB/J34NNr8WAnVpbP4Gb88wgav/Enj96xBohO3vgvHTc/AvAoSqK2AVtm7JUM4qABEZzFn1RkKZbwI3qOpqAFX9KNMOVXUwVUgnkMM/w1Iyel6Gzn+A1MC4Y6F+Wvqy0X74x9fgndkmPdWml8E/fgbBOvjSAzA1l0VxqxAfrAfkVlOjJwDF34PeUyH+HFADNadA3TUgaYac/e9mWHo2xHuGb++4HfqXwnaPMiIlwiePw5tfg/jgjO04xLrM67UjYNfnoGXX1PXFFkDfaRB/FWiAmu9B3c9zC5SVQP7xNFXOqr2SymwNICLPY6R4iao+mm6HzoidbwKbkXDtO+3WlnKgCstOh9V3mxs6AvDfS2GDH8JGP01t8+Cp8M4/IJqQdifSBXTBrEPh9AWw3lZp6ouC/gF0PsinQL4DkubmsRLxxyi4gjXliwAUDocLGgefs50OQM9+TgbdONCPDvwZ0T5oSDGcPT7gPPn0jPxOe6HreQjPgZZ9hn+35LyE4JO8z1547yLYKcUoWv0EetqAtc6GbnTgGkRqoe6S9MflUMj5K/k5Tya1UCaISOKEtpn55K3CXL/TgDZMSpFnRGRHVV2TpvzfgWcx2WGzZL+sfHx5Xay+x7yGtBUz+vzfVTDmEGhO0lTnMtPcFkuT8y3aB/+5Aj5/88jvVCF+NNAO9BCP1xOP/5VQzXM53diV6/wVrCnwPADhUlOeB6CyzMqPPQ4aJjHVjdCHRmcheiNI0/DyXc8zPC1OEvEe6LhjeADq/xB6MuXgU/OEFI+ObLceuBuT/HgdQg8a+R2SJQD5O+fX0MGkaypYlSEVT7acVWCeil5Q1QHgPRF5GxOQXkqzz0ZV/WF2hysf314XH10P8e4RmzXei6yaOTIAvZ1lYr5G4c2/pQ5AvIzJYmSCXUD6UF6ht/cxGhoPL94xldlmCB80weFSU57Hz3A4jIjQ0tKCiOScETYvO/2YlAFFAR0pBiOQTFMwFGKdSTY9IDnE81RrCunHpFx2Q7Ofi0LOX1nOeTL5j9Z5CZgmIpuLSC1wApD8a/QA5ukHEZmAaZJ7N8M+HxSRzL88VYJvr4vY6pSbBYWBFNNHoj3mpi0Tsf40X3SS/AutBOjry9hVCJTv/LnSFPhhFJwrTXkegJqbm1FVurq6UNW80n/kbBc8iOQnDACVqSApUoA07uK0T6ch0AQt+w/fVjc1+8CBug0hmKLPKXQoCZOIAYirEJP9Mu+Pws5fWc55IgUMF1XVKDCYs2oRcM9gzioROcop9hjwsYi8ATwF/EBVP86w27MxgukVkbUi0iUiazOUr1h8e12MPRykdsRmlUYYe8TI8pN3gZos6RLWT5ceaxeSLzRBqa31j64K1hT4ZRi2K0153gRXlvxTgY2h9nKI/BhFUQ0iEiTQeNfIgQQAdRvDuINhzWPOqLQkJAQTT0qqoxY2+hasuCF1P1CgEaamWX492Ao1M2DgDyhxVEMgDYSabkl/TA7+zvmVQAFNBdlyVqlJ5X6u88plfy35e1GZ+Pa6mHwurPqj8yRkugyUWqRmMqx/8sjymx8MtWOcQQcpqGmCT6fRlYyDwD8h/kVU3ycen0R/9HaamtMMWCj0mMpsMwzvR8G50pTnAQjKlH+q7hyoORKJPoRIC4S+CJJhTs9Wf4aFB0LfuxB3HosDjSb4bPcYBJtG2mz+C1j7AnTNH97OHWiC9Q6Djc9MX1/9b6HmNCT2T0QmQ+gYkNwSZfk75xeeJyMdRMxKj8l0AsucJ66qwpfXRc0GsP08WHERrJkNEkLWOxE2+llqTUkAjp8Nf2kzAw7iAwn7aoJtj4Ntv5S+PtkNgksRjRMMBWgc+fDl/pg8sAF8oSu3mvJFACobgS2hNselX0LjYKf5sPph6PiLCShjD4ZJp5jvUu6/Dnb+F6x6AFZcB911MP4QmHK2CUCZV7KD4I7mVY34IAABNwK7Aa85n3cEXgfGisiZqvq4Z56NJmo3gc1vz738BrvBjIUw59fwxp1mBNtGe8E+P4Stj8muK6iueXWJeH9YrjTlvft+RoKw3pGwzd2w3YOw0dnpgw+YYZ//eQFm/AUOewveicJhr8E5d8Err5TNbd/hj7ZqMCsM7qqqu6vq7phOgncx2Rau8sQjS26MmQqH/A7O+cj0C506B7b5Qm7Bp1rxh65cacoGoGKxcCFsvTV87nMwezb8738Qi8GHH8Jdd8G++8Juu8H773vtqTd4nzQRYOvElR+dtD7bDmZbsPicjz+Ghx+GNWtg3jxzwzfa8V5XrjRlA1AxmD8f9tkHliyB7u6RwojFoKcHXn3VBKF3k/438Si88wDMaoM/bAK3bA1P/xDWfkBV4I87NYCFIvJ7ETnQed0IvCEidaxLqmjxG0uWwHHHwZQp8JWvwNKl0NYGW2wBt9wyegORP3TlSlOjqw+oEPoWw5r7IN4HzftCy2eGP/avWQMHHwxdaUbpJBKLwerVMH06LF4MNTXQvRLung5dy2EgYQ7A/GtgwXXQ9hvY5YyiH1bZ8cetzqnAt4FznM/PA+djhGJzw/mRl1+GAw+EcBjicejrMzoKh83r7LNhzhyYOXNkc1ykG+b9Hl66Hrr/Bw3rw+5nwh5nQX0OSYUrAe91dSouNFWUACQihwLXYB4Ab1bVK4qx35KgMUwy0ixtxxqDZd+ET+5ybKJmNFvtFNj6X1CzoSl3661GFLkSj8Mnn8ADD8AXj4ZZB8Kad0dOUB2cXNd+LjROhK2PTb2/2FsQuQsYgNrjIbhz7r6UC3/M2EZVe4FfO69k8pwBWHoqSleFoFFA0qfF6e01N3drM0wr6e6GO++E1lb41rfWbe/rhFv3Nql8os60iK4V8NwvYMH/wekvQlO25Ok+xwe6cqsp1/EzIW3+YZhlWU8UkXQzw7xDOyF+CGgdaD3EL8786P7hZfDJ3aB9mGCuZjh232J4+3NOnqk4XH21aV7Lh3AYrrwS3rkfwv9NnR1hkGgvtJ+X2tf+W6BrV+i/DPqvhK5PQ9+V+flSLjxsKhCRe5y/r4nIq8mv8nmSOxWjq0LQHogdA7F65/Ut5yYvibvvzu3mrqcHfvGL4Rp57Hvmxi6aNCcv2gdd/4V/nO7qEHyDR7oqlqaK8QSUS9r8jJQlMaaehMkLFQNiqP4aYVOQFBeiDsBHvzFCGUEUIu9C979h2bjMd2eZeO01mHPl8Ga3dPR+DB/OgY0S8mRpJ/R+F0gUWA/0XQK1XzWTbzPgk2Sk5eJs52+Kqfa+xRNdleW6iJ8J+hhDE1HjtyNsBMGkbNjXXGNu1nJhzRr4z3/g05+G/rXwxj0QSzGJHMxconefgPBKaN6gOMdUZpshvNNVUTRVjACUS9r8tJQvMeY/gXUXpEgP0eg9hGpTBKD+ZanvyAaJD0D3HFjVCqECT2FtLXS+l3v51e8MD0DRfwM1DA9AoASR6D+h9mtpd+WjZKRlQVU/dN6uAnpVNS4iWwPbAo9451lGyq6rsl0X+giJuQ9FeonF7iWYHIA+yGMQjggsW2YC0CfvQLB2+PINyYTqoWNh2gBkk5FmpliaKtsgBBGZAcwAmDx5Mu3t7QBEo1EikQjBYJBYLEZnZyehHH7Us9mFw+GhOgxXkJwPTuPjkUA7I9AB6P0ZJltpKgKwegr0h+GSS0ynaArCU6bQfvXVaXYRgI0C2RMtgplE90ELdCT6GofYJYxMshqAwFiQdWWTz0WpznlGvO8sBXgG2F9ExgOPYxKefhk4KaOVjymmrnKxGamrQuq6GBiea1HjzSO1ePHFMDByIFVKXQWD0NAA7e0m8GxyCWiGjPYSgMV9sKw95df5/76U7pxnxHtdudJUMQJQLmnzcdZ5mQnQ2tqqbW1tQOnuutrb2xmsA4D4YlS/h0gvqkI83kBf5BmamnZPXcEbP4DeNJNHpR52XAofK3zpS2nbqduvvpq2889PvY+GBrjnWHjrrsxPW2ASmB7+LjRvuG6bKnSdg8bfQwabMlRQmUBg7Aqz0uSgH0nnwpMnIO+FAiCq2iMipwM3qupVIvKy106loey6ysVmhK4KqSv+MRo72dEixOP19Ef/QWPT8P3y+9/D3/5m+loTfUilq/p6eOcdM1Rb43DN6aavJx314+G4/0GwJuXXef++FHIeCrQZwh+6cqWpYgSgobT5GIGcAHwlV+OyJcYMfAPRjRiI/IXIQD3Id2lqSrM6KcAmN8A7h4xclC7QBBO/BzWTYQPMBNN//jMnn4cIBuHEE2Gvc2DxvSM7ShORIGx60PDgA6bJofmfSPeX0NirZoyEbEGw5b5hwScVlZKMtASIiOyDuTsbbHv1h2cjKbuuyqfFYxHGMxC5mYEoIN+msSlFhurzz4cHH8xtkM8BB5jgA+bp5jOXw8NnwkAK25pGOODitMGnoGMqo80wvL96XWnKdQBS1aiIDKbNDwK3Js6MzYWyJcaUw6mpO5yaXFbkbd4Xpv0Tlp8DPfOBIATHw0YXw4SE4Z4XXGDmIXSnWFcoHbW18P3vw8QdYPfvw7zfmXVPRvgbhPr14OA/pN5PYBNoeRGJrwRiBLMMPBh2eKMwGSmm4/THwP3O0g5bYJZx8B1e6aps10XgM9TUf4aMt0p77GGGYT/+uBmSnY6mJvh10ijgnb4G4Y/g6f8HiLnJC9WbVvW9zoU9z061p2HYZKQ54UpTRekDSpU2vypo3hu2nQPRNWY4dmjSyKSGBx9sJpY++WRuQ0YbG+FrX4MddjCf9/sFNG0Iz19shmPHBkwdGoeNPw2H3QYtWQJLIP1IHt9QgFBynQcjIscCfwP2UNW56fanqs9g2qwHP78L5JidtvxUra7y4e674fjjTStDT8/wodZNTaY14eGH1+kpkU+fD7t8HV6/EzqXQssU2OHEjCPfKg6PA5BbTdlMCLmQKQGpCPz1r3DkkfDvf2duLmhqgmOOgeuvH26/21mw87dg6WNmZFyoATY9GMZuWqwj8JYCRuskzIM5GDMC7CURme3kmkos14K5C3uhKL5a/EVdnZm0/Z//mDl3Tz9tBvBstRV873vmZm5shqwGjevDnt8tm7tlxQcTUd1iA1AxqK+HRx+Fm26Cq64ymQ76nVE+tbVmqPamm8KFF8JJJ6XOwhCsgS0raZpKnuR/p5brPJifA1cCP3DpocWviJjh1ffdZz63t5sBBxbPn4DcYgNQsQgG4TvfgW9/G559Fl58ETbaCC6/3HSQtrZ67aF3FNZWnXUejLMY1lRVfUhEbACyjC780QfkChuAsqFxGJgHGoaaPSCQrbMwBtstgvVvgfe/BbvcDpM3AN09t/xzsaUQGAeB9Yt0AD4htVAmiEhin81MZ1hxVkQkAPwGkwwxW9nrSD+pC1X1bT+QxUEV1vwblt8CPfvC0ldh41OgpkqSihaKRwGoWJoaXQFIl2EyIjQCR4KkWAI4kch8+Pho0DWY/3QUxlwFzd9Js/8oLD4cws+bND4ag96X4f0ZsPZR2PS29EGo9y4If9fJPReF2s/CmNshML6wY/UT6duqV6lqukfDbPNgWoAdgHYx53QDYLaIHJViIMLg530xedXudj4fRx6pbSxFJD4fog+YJe6DJ0JgWvqyqrDwW/DhHRDrhein4J2LYcmlsNfz0LxNBtt+iN0FsecgsAeETgZpLPrheIK3fUBF0ZQvAlB5csE9D3wOVXUG0kwkEFgAkuYHPt4Fqz5jcq4lsvYCCE2D+kNG2nxyp8kRl5xDLt5tlnRY7xQYc9BIu8iz0PUNYJ2dRp5AOo+C8c9mPi58nPMrkfzv1DLOg1HVTmDC4GcRaQfOTzUKTlVvc8qcCew3uFa9iNwEZD/BFYpvr4vYQxA5DtU+IAADVyH17RBIcy/y0d/hwzshlqCrWI8JRgu+CPunGZ2uMehvg/hrQDfx6J1o5FqCjfNBss/FqNZccNlGl4rIucA3MKljOoDTVHVZYpliacrzAFS+WflfA7oRMQ8h8fh/iUQuprbuutTFe/+aOku19kDXFakDUMe1JtikIt4NHTekDkDdvyQx+AAIEXRgHhJ9E0Lbpj0qX+f8SiRPoaSbByMilwJzVXV2fnsEYDwwBvjE+dzsbKs6fH1dRE4Dep3GgBjQTazvGwQbX05dfulvIZZKVwp9y6DrdWhJMQw7/jjEXweMbUB6icXeo6/7Tuqbv17cYyqjzTDy1FWOo0sXAK1OhoMzMUtrfznNLl1pyvMurHA4jIjQ0tKCiBDOMftt/nYrh30KBAaIx99OXzy6jMELdwSxZam3D/wvswsDy1Nvj6devVa1Jn1dDoWcv/Kdc4fBpoI8lw5W1YdVdWtV3VJVL3O2XZwq+KhqW6Y5QA5XAAtE5E8ichswH/hlbgdRWfj7uvh45CbNkDanL41uwDTh9aexjb9N8qKcgUAv0YHs83nLdf4K1hQUqquh0aWqGgEGR5cOoapPqQ4148zBNH+nw5WmPA9Azc3NqCpdXV2oal7pP/Kz2wnVdf+dWKwekf3TF6/dFaRlxGYlALV7pLap3y5D/UFoSLNYXM2nSXXliPRDaKcM+yzs/JXvnCfg/dLBqOofMSPp7gfuA/YZbEqoNnx9Xcg0VNf1hcY1iEqG63zMLphf2xTE+6EpTQtBYF+SdRWPNxCsOTCzf5Tv/LnSFBSiq1SjSzPNdD+dDNmt3WrK8ya48uUluxuRA1HtAKKofpa6uh+lL15/hMkwEBtckM4gUg8tP0ltM/kC6H5+ZP44AKmFSWnSfzReCP1/MyPtHJQGpP6rENwwtY2Dr3N+DeKT4aJiRit8FthCVS8VkU1EZE9VfdFr34qNr6+L2nuR/gNQHUCJA+MJNfw5ffnNfwgdj47UldTCem3QsElqu2ArhE6H6M3OzWeMeOAYGhqzz7eriFxw6XVV8OjSYbsX+SrQCqSN2G415XkAgjLln5JNQN9GZDHQQCi0aeZh0RKCic/DmjOhbzYQh9AOMO5GqEnR3gww5rMw6Tz439Wgg2sP1Zh9TfktNKSxC02D8f+B8A8g8hwExiENZ0Pj93M6NF/n/AKvR+skciNm/YrPAJcCXcC9QJpH2srGt9dFYHuoX4zE5yAEzZNKppFp4/aET90IC890lu8WCDZCy86w812Z66q9FkLfRuLzILAjNYHMLQoFH1OZbYBSjS41uxb5LHARcKCq9id/n4ArTfkiAJUNqQEyNZMlEZwI6//NrA9EFKQhu81Gl8K4Y2HV76FjDEz8Dkz8NtRnGGYKTnDz69poRcAHT0DAXqq6m4gsAFDV1SJS67VToxIZB8FDcy+/8Skw6Rj4aDYsrIM922FsjvcNgW3Nqxop8uhSABHZFfgDcKiqfpRlf6405Y+fBb8jNbkFn0Ead4ZNboK6aTD1t9mDz2jAB31AwIAzCkgBRGQiI1f0s/iVmrGw8clQOyn34FPt5KkrZ7j04OjSRcA9g6NLReQop9ivMKPZ/ioiL4tIphGnrjQ1up6ALN7gnya4azGdpZNE5DLgS0CaDj2LxecUqKtUWdZV9eKE95/NY3euNGUDkKU8+OBZW1XvEJF5wEEY+R6jqos8dstiKRzvl2NwpSkf/CxYqp7B0ToeN8GJyC1AvareoKrXq+oiEbmk/J5YLEXAB7pyqylX7orIr0TkTRF5VUTuF5FxbvZnqWJ8EICAzwG3icjXErYdla6wV1hdWXLGe1250pRbd58AdlDVnYC3MUuzWujALJG+CDgDMzIxG1FM0+nOmGH1r5TMu7JTYCaEEvARcABwnIjcICIh0s5w9BSrK0t2/KErV5pyFYBU9fHBJHRkT9mQlnA4zMqVK/NLQ1GgXelt+oC9Uf0r0IPqnzDNo9kGhnwL+C3wKqr/RHU/IHWKHnf+FW7jxs4Hd2oAoqqdqnok5i6hHfBdPn8vdVUJWoxGo772r1znHPCDrlxpqpiDEE5jXUrunClnYszy2LyAagciJnuCSD+qCxFZAqQbjh0H/ox5ChpMltpHJHIndXWZB5T4OunkID7JhAAMDSdV1UucztPcZvt6R9l0VSlajEQiLFq0yLf+lS0ZqT905UpTopp2TSFTQORJzForyVykqn93ylyESdnwRU2zQxGZAcwAmDx58u6zZs0CIBqNEolECAaDxGIxamtrCYWyx8VsdoOpLfKxKaSekYSBd4A44fAUmpuXY66S7YFMKeDnM3x9JyEe34BAYCPX/iWfi1Kc8+nTp89LN/u6dWvRuTeM3C6HkNam2vGjrgq5lkpZVzqbaDRKKBQq2m9FITblOg/Vrqus/71sY8JF5FTgCOCgdCJx9jMTmAnQ2tqqbW1tQOnuutrb2xmsI1ebQuoZST/x+LnAUp555nL23/8iVHcjFHqezE2j96B6GyI9qArxeCN9fS/R1JQ5c0Mu/iWfi9H2BCQiz6nqfiLSRXKUB1XVMeX2yY+6KuRaKmVd6Ww6OjqYOHGip09AXp6HITzUVbE05aoJTszCRhdg8gWlyMCZnXImxiyPTR2BwAtEIueh2kQs9k1qa68ie7/cdYhMIRZ7gIGB8cRiv8wafMp3TO7sAE8DkJoONVR1ZHpzH+KVripFi52dnXn9UPv3t8KlpsAzXRVLU277gK7HtCs9YZKiMkdVz8h3J+VMjFkem/HU1t6KSDu1tbmejiBwIcHghQTzHMni26STgxQ4Y1uKsHKjU269TPWo6ieZvvcAz3RVCVoMhUI+1b1vkpGWnGJpylUAUtWt3NhbRhHertw4D9NMkOoxVIEt8vOutFhdWXLGu5aFomjKpuKxlIf8hTK0ciOAiAyu3DgUgFT1qYTyc4CvptqRqm6ed+0WSyXgXRNcUTRlA5Cl9BTWVJBq5ca9MpTPuHLjkCsi4zHj4esHt6nqM3l7Z7F4jU+S/LrRlA1AlvLg8cqNTrlvAGdjJna+DOwN/AezmJbFUnl4PA/Iraa8n8ZkqX7SJ01cpaqtCa/E4JPvyo1HZVm5EYxQ9gCWqep0YFdgTSGHZLF4jg+SkeJSUzYAWcpD/jmrhlZudFZYPIGEWdcwbOXGo3JYuRGgT1X7HNs6VX0T2Cafw7BYfIX3ueBcaco2wVlKTwET5lQ1KiKDKzcGgVsHV24E5qrqbIav3AjwvqpmysS73Mks/QBmiPNqYMSwbYulIvBHKh5XmvJFAAqHwwVNxCrErpw2g0kTq+WY3NgVIpQir9yIqn7BeXuJiDyFSZr4aP6eVQZ+vy6srgq3GcL7BelcacrzAFQpCRBHe9JEN3Y+uVMDhkbsTMWskdEF7IBJxFdV+P26sLoq3GYIn+jKjaY8dz8cDiMitLS0ICI5pyQvxK7cNsFgsGqOyY0d4Ie2akTk58CrwHXAr53X1eX3pPT4/bqwuircZhge68qtpjx/AmpubkZV6erqQlXzyu+Ur125bWKxWNUckxs7v9ypAccDW6pqxGtHSo3frwurq8JthvCHrlxpyhcBqBISII72pIlu7AA/CAXgdWAcZhXHqsbv14XVVeE2w/BeV6405XkAgspIgFiITbUlTSzYzicztoHLgQUi8jowNGcoy8i5isXv14XVVeE2gF905UpTvghAllGA93dqALcBVwKvkX2NdIvF/3ivK1easgHIUh68FwpAj6pe67UTFkvR8F5XrjRlA5Cl9PijqQDgWRG5HJNRIbG5oOqGYVtGAf7QlStN2QBkKQ/e36mByVMFJmHiIIpNRmqpVLzXlStNFSUAich5mLHfE1V1VTH2aakifDBc1Fngbraq/tZbT3LH6sqSEY91VQxNuXZfRKYChwDvu92XpYrxeMKcqsaAE8tba+FYXVlywkNdFUNTxYifvwUuwDx2FUQ4HGblypV5zwIuxK6cNoM5q/LBz8dUsJ0/0sYDPC8i14vI/iKy2+DLE0+y44muKkGL1aarQs+5T3TlSlOumuBE5Ghghaq+4mQjzptKyT812nNWubEDPG+Cc9jF+Xtpwjbf9QF5patK0WI16cqVpsAPutrF+VuQprIGIBF5EtggxVcXARdimgmyIiIzgBkAkydPpr29HYBoNEokEiEYDBKLxejs7CQUyh4Xs9mFw+GhOnK1KaSeTDbRaJSOjo6iHVOhNsnnolTnPC0+6AMCcBbM8gV+1FUh11Ip60pn4wddeXkehvCBrtxqKuuRpkt5LyI7ApsDg3dpU4D5IrKnqq5MsZ+ZwEyA1tZWbWtrA0p319Xe3s5gHbnaFFJPJpuOjg4mTpzo+Z1a8rnw5AnI++GiiMhY4KfAAc6mp4FLVbWz3L74UVeFXEulrCudjR905eV5GIbHunKrqYKb4FT1NWBSgiNLgdZ8R+tUSv6p0Z6zyo2dH+7UHG7F5K463vl8MvBH4IueeZSEl7qqFC1Wk65c5YLzh65cacoX84AqIf+UzVnlzs4HQgGTtffYhM8/E5GXvXKm1Pj9urC6KtxmCO915UpTRXNfVTezcxUsKRmcse3xekBAr4jsN+SWyL5Aryee5IjVlSUt/tCVK02JasGjPAtGRDrIY93wApkAeC1cP/gA5fFjU1WdmOoLEXnU8SGZVap6aGndGubHzsCfMcsGC/AJcKqqvlIuH0qJ1VVV+uBrXbnVlCcBqByIyFxVbR3tPvjJD78gImMAVHWt175UGn64lqwP/qNQTfmiD8hiKQciUgccC2wGhAbn2KjqpRnMLBZLGtxqygYgy2ji70AnMI+EzL0Wi6VgXGmqmgPQTK8dwB8+gH/88Jop5exzqlL8cC1ZH/yDK01VbR+QxZKMiMwErnPm2lgsFpe41ZQNQJZRg4i8AWwFvIdpLhBAVXUnTx2zWCoUt5qqmgAkIr8CjgQiwBLg66q6JkW5pUAXEAOixRjJIiKHAtdgRuHfrKpXJH1fhxmquDvwMfBlVV3qtt6E/U919j8Zkwhwpqpek1SmDdNe+56z6b7R1vkuIpum2q6qpR66XLFYXVldZcK1plS1Kl6Y5I0h5/2VwJVpyi0FJhSx3iBGmFsAtcArwPZJZb4N3OS8PwG4u8jHviGwm/O+BXg7hQ9twINe/5/sq7JeVldWV6V8eZ/IoUio6uOqGnU+zsEkcSwHewKLVfVdVY0As4Cjk8ocDdzmvP8bcJAUmmc/Bar6oTprsKtqF7AI2LhY+7eMXqyurK5KSdUEoCROAx5J850Cj4vIPCeVvVs2Bj5I+LyckRfpUBlHzJ3A+kWoewQishlmnfYXUny9j4i8IiKPiMinSlG/paqxurK6KioVNQw70xoqqvp3p8xFQBS4I81u9lPVFSIyCXhCRN5U1WdK43F5EZFm4F7gHB05I3k+Jq1HWEQOBx4AppXZRYsPsbrKjNVV6aioAKRp1lAZREROBY4ADlKngTbFPlY4fz8Skfsxj/puhLICmJrweYqzLVWZ5SISwuRN+thFnSMQkRqMSO5Q1fuSv08Ujqo+LCI3isgEtYkuRz1WV+mxuiotVdME54yYuQA4SlV70pRpEpGWwfeYDtbXXVb9EjBNRDYXkVpMZ+jspDKzgVOc918C/pVOyIXgtHvfAixS1d+kKbPBYPu4iOyJ+d8XVayW6sPqyuqqlFTUE1AWrgfqMI//AHNU9QwR2QgzhPNwzHDK+53vQ8Cdqvqom0pVNSoiZwGPYUbu3KqqC0XkUmCuqs7GXMS3i8hiTLbYE9zUmYJ9MQtBvZawFseFwCaOjzdhBHqmiEQx6dJPKKZYLVWL1ZXVVcmomnlAFovFYqksqqYJzmKxWCyVhQ1AFovFYvEEG4AsFovF4gk2AFksFovFE2wAslgsFosn2ABksVgsFk+wAchisVgsnmADkMVisVg8oZoyIVgsKRGR5ASSI4oAH6rq1uXwx2KpdIqlKRuALKOBJaq6a6YCIrKgXM5YLFVAUTRlU/FYqh4R2UJV33VbxmKxGIqlKRuALBaLxeIJnjTBTZgwQTfbbLOS1tHd3U1TU1NJ66gEH8rlx7x581ap6sRU3x0qknJxlHnwmKoeWlLHsiAiM1W1GCt4eo7VVfX5UIm6yktTqlqUFyZl+gLgwWxld999dy01Tz31VMnrqAQfVMvjByZFfur/N6gGRr4y2ZTrBezutQ9Z/LO6GsU+VKKu8tFUMYdhnw0sKuL+LNVEIMXLB6jqPK99yILVlSU9PtRVPpoqShOciEwBPg9cBpybr304HCYcDtPc3Exzc3NJ7cppE41Gh+xyxc/HVLCd4AthiMhTwIhOT1X9jAfuZMULXVWCFqtNV4Wecz/oyq2mitUH9DvMsr0t+RqGw2EWLVqEiKCqbLfddjn9EwqxK7dNJBJh0aJFVXFMbuwA05CUJ85y0Nc41jer6hUpyhwPXIIRwSuq+pUMuzw/4X09cCwQzd+zsvE7yqirStFiNenKlaagIF0VGVeach2AROQI4CNVnScibRnKzQBmAEyePJn29nYAotEokUiEYDBILBajs7OTUCi7W9nswuHwUB252hRSTyabaDRKR0dH0Y6pUJvkc1Gqc56WAu7URCQI3AAcDCwHXhKR2ar6RkKZacCPgX1VdbWITMq0zxRNA8+LyIv5eVYevNBVIddSKetKZ+MHXXl5HobwwROQW00V4wloX+AoETkcEwHHiMhfVPWriYVUdSYwE6C1tVXb2tqA0t11tbe3M1hHrjaF1JPJpqOjg4kTJ3p+p5Z8Ljx5AspfKHsCi9WZRyAis4CjgTcSynwTuEFVVwOo6keZdigi6yV5tDswNm/PykPZdVXItVTKutLZ+EFXXp6HYXjfBOdKU64DkKr+GHMXinOndn6ySDLR3NzMdtttl3cbaCF25bbp7OzM64Ly8zG5sUMopKlgY+CDhM/Lgb2SymwNICLPOzVcoqqPZtjnPExTnWCaCd4DTs/bszLgha4qRYvVpKuCNQWF6qrYuNKUL1Lx5H3iXdiV0yYUCvnav3KdcyDdndoEEZmb8Hmmc0efKyFgGtAGTAGeEZEdVXVNqsKqunke+654/H5dWF0VbjOE901wrjRV1ACkqu1AezH3aakC0rdVr1LV1jRWK4CpCZ+nONsSWQ68oKoDwHsi8jYmIL2Us2siG6jqylzLe4HVlSUlPugDSkU+mvKh+5aqJJjilZmXgGkisrmI1AInALOTyjyAefpBRCZgmuTyzed2S57lLRb/kL+uEJFDReQtEVksIj9KU+Z4EXlDRBaKyJ15epWzpnzRBGepcgq4U1PVqIicBTyGkdWtqrpQRC7FzPSe7Xx3iIi8AcSAH6jqx3nW8/n8PLNYfIJPRpcmk4+mbACylIcCnrVV9WHg4aRtFye8V8wEzbwmaTqCqk/Yz/v5e2ex+AAfjC519lOQpmwTnKX0DI7WybOpoOhuiBwlIu9gRuo8DSwFHim/JxZLEUivqwkiMjfhlZgYNNXo0o2T9rw1sLWIPC8ic5wJ4aldcKkp+wRkKQ/+uNX5ObA38KSq7ioi04GchzZbLL4j/8E9uZDP6FJXmvLHz4Kl+vFH0sQBp48oICIBVX0KcCNUi8Vb8tdVrqNLZ6vqgKq+BwyOLk2FK035IgCFw2FWrlxJOBwuuV05bQaTJuaDn4+pYLvBzlLvA9AaEWkGngHuEJFrgG5PPCkDfr8urK4KtwEK1VWxR5e60pTnTXCVkgBxtCdNdGMH+GHGNpjO1l7g+8BJmJQhl3rqUYnw+3VhdVW4zTDy1FUJRpe60pTnT0DhcBgRoaWlBRHJ+S6gELty2wSDwao5Jjd2fnkCUtVuVY2ralRVb1PVa/Mdtl0p+P26sLoq3GaIAnWlqg+r6taquqWqXuZsu9gJPqjhXFXdXlV3VNVZGfblSlOeB6Dm5mZUla6uLlQ1r/xO+dqV2yYWi1XNMbmxAzwNQCLyYDHKVBJ+vy6srgq3GYZHuiqWpjxvgquUBIijPWmiGzsfJE3cT0SS27kTEWD7cjlTDvx+XVhdFW4zhLe6KoqmPA9AUBkJEG3SRHd2Hj9rH51DmUjJvSgzfr8urK4KtxnCO10VRVO+CECWKsfjpImq+rR3tVssJcJDXRVLUzYAWcqDP0bBWSzVRYXrygYgS+nxadp4i6WiqQJdVbj7lorBB8OwExGR8SKyk7deWCwu8ZGuCtGUa3dFZKqIPJWwdsTZbvdpqTL8k4y0XUTGiFnHfj7wfyLym/J7kh2rK0tWfKArt5oqRryMAuep6vaYpHTfERF/DmnVdoifBfEfgy7LzSY+Dwa+DAOHQ+yPoJpDPTHQW4EVoLNys6l2/HGnNlZV1wJfBP6sqnsBn/XEk+xUjq4s3uG9rlxpyrW7qvqhqs533ncBixiZ3jsjZck/Fb8L9PPADaj+Co3vBJpl8cz4ixA9AI3fA/oIGj0LYpdktlEFjkT1u8BKVE8HZmS2ScDv+acKzlvlvVAAQiKyIXA84OuJp17pyuaCqwybIbzXlStNFXUQgohsBuwKvJCrTfnyT30f6HH8jKHaRSRyKbV1f0pvErsY6EEEx64HjV2OBC8GSfes+wKqzyDSs85Gb0f4fyCbFPmY/J/zC/B6wlwil2LyXD2nqi+JyBbAOx77lJVy6crmgqsMmyH8oStXmhItUvOQkxH1aeAyVb0vxfczcB4FJk+evPusWSa9UDQaJRKJEAwGicVi1NbWEgplj4vZ7AZnFq9jARAftg+Nj0EC6bKMA/omKRO7yq6kv9VYCywB4oTDU2huXu6U3RZocHVMhdokn4tSnPPp06fPS7cGSesY0bl7j9wuT5DWxmIop64KuZZKWVc6m2g0SigUKtpvRSE25ToP1a6rojwBiUgNcC9wRyqRAKjqTGAmQGtrq7a1tQGlu+tqb29nsA4A4r9D9VFE+gGIxeoZiF5LfUMbaYk+g8auQKTX7CIeQNmWYP3C9Da6CtWvItJF+9NXc+ABP0B1EoHAUpD69HY5HFOhNsnnwpMnIG9H51wHpL3TUtXvldGdnCm3rgq5lkpZVzqbjo4OJk6c6OkTkJfnYQgPdVUsTbkOQCIiwC3AIlXNe0RR2fJPyZ8Rjkf1SSBINHYu9fXfyGwTvBDRxWh8FhAAmUqw9tEs9UxA+Bfx+HGAENdtCQYeyBp8CjqmMtq4sQO8HnY919PaC8ALXdlccJVhMwzvdFUUTRXjCWhf4GTgNRF52dl2oao+nOsOypJ/SsaAPIpoDAhQVy852ISg5s+IXgf0IkxmqEMoo10rgeB7QDvB4Bu5+efg9/xTBdl53Fatqrd5V3vBeKIrmwuuMmwAT3VVLE25DkCq+hzmVFQGaQcPZLIZi1lnyVIwPpjyLCITgR9isvQOPZKq6mc8cyoNFacrizd4P6HblaZ88LNgqXp8siAdcAdmOPPmwM+ApZglii2WysMfunKlKRuALOXBB5kQgPVV9RZgQFWfVtXTAN89/VgsOeO9rlxpyiYjtZQe/yRNHHD+figinwf+C6znoT8WS+H4Q1euNDW6AtDAJ7D2BQg0wNh9IVCT3SbaBZ88AbFeGPdpaNi89H5WIwUIRUQOBa7B3NfdrKpXpCl3LPA3YA9VzTQ65xciMhY4D7gOGIOZoWwpN33vw+p/mT7Z9T4HtZO89qgy8T4AudJUZQYgVYg+D9F5ENwGag7OPLhAY/DO2bDyFpBas00CMO33MPmE9HbvXweLfwhSAyjoAKx/GOxwBwQzTCrVKEQegXgPRJ6Ems+Y+kYrBYzWEZEgcANwMLAceElEZqvqG0nlWoCzySFLgKoOpgrpBKbn55ElLfFVEH8JZHMIbpul7AAsOhVW3Ye5KASIwsbfgS1/lXmU6cDb0D0TYp+Cng+h4UuONkcpPsiE4FZTlferqAOw9jDoPBS6fwhdx8Ga3SDeld5myfmw8o8Q74PYWvOKroG3TofVT6W2+eg+WPwjiPc6Nl3G/uNHYOEp6euKd8LqnWHtSRBfAZ1fgDUHgvZnObAwJp3SeEx/3mNZylcY+XeW7gksVtV3VTUCzCL1MsA/B64E+rLtUEQmisiFIjJTRG4dfOVxFJZkItdC91ToPRF6doOeo4xG0/HOd2HV/UZL8W6Ih837Fb+HD65Ob9f7D/hoF+i+BuKfwJpvQkcbaLZVn5cBnwbGAbsDb+V3fH7H40EIbjXliwCUVzK+/tth4DlMipx+0DAaewt6f5m6fHQt/Pcm8zSSTLwH3rsotd3in6Sx6YNV/4C+5antei6D2BKgCzNROAzR+dD7f1kO7CRgNrAGWIrqF4DXs9gY/J50ssDROhsDHyR8Xk5SMk4R2Q2YqqoP5ejJ3zHj6Z8EHkp4VSUlvy7ib0P/jzCxvxPoRWNPQuT61OUHVsPK28xN3Yh99cCyyyEeHfmdRmH1qUAvJkk4oN0w8Ar03JXBwQhwAObhuBPVBcTj+2Fu9rLj+2Sk/hgF50pTnjfB5Z2Kov8ekvOzCf3E+/5KoOnyFBW8AlJH2hvkrnkjt8Wj0PNmeh+kFta+CPVTUvh3H5D8tNMD/XdD41np98mjGMEYVKNEIg9SV7dDBhv/J50cInVTwQQRSeyzmemklsmKiASA3wCn5uYAAI2q+sM8ylcsZbkuog+TnF9R6CUWuYNgXYpugGxajEeg/31o2CKpnjdJ1MY6uqHvXmhK1yKxGPhkyEcRJRbrpa9vDo2NmVcMqIhkpFBQE1yR+1ZdacrzJ6BwOIyI0NLSgohkvwsITCTV/LxoLM1E0WCz6QNKu78UKXIkYLIgZCLYlGZ/40dsUgUC2TpZG5NsgvT11WaxKeD8FWjjxi7DndoqVW1NeCUGnxXA1ITPU5xtg7QAOwDtIrIUs2bObBHJlITxQRE5PDenK5uyXBeyPjC8D0YVBqIjNQA4momn/g7Mk04qXclY811ycQLO70E6xrBukJazK4nR05P9vrtcuipYU1DQE1BC3+phmMmjJ6ZaZyqPvlVXmvI8ADU3N6OqdHV1oarZo3/92SRnlY7F64nVnZemgl0glDo4KTUw+Ssjv5AATDyajKdnfFvq7Y0/RJOCCdIAjdkGhvwa1QZUIRarZWBgAsHgSVlsCjh/Bdq4sQMKaSp4CZgmIpuLSC1wAqaNEgBV7VTVCaq6mapuBswBjsoyCu5sjGB6RWStiHSJyNrcD6JyKMt1EfoCUG8CgUNc64iHzk9dvmV3CDSm/g6gaQeonZyinqlQuxs6osGmDprOzODgFOAkVE2dsVgDnZ17U1+/ewYbQ7l05UpT4Ie+VVea8rwJLu9kfDWt0PIXCJ+F6keoNjNQezENY1MEEjCjarb9E7x+tNP2bBK4qtQiNevDZpektpt2FXzyT4h2MuyuLdAA294AgbrUdnVfQppWoN0XA4KyPtJyDdTsl/m4OA2RLYhEHqKvr4VA4Js0N6cQYxJ+Tzo5RJ63OqoaFZGzMKMxgsCtqrpQRC4F5qrq7Mx7SLnPlnxtKpWyXBfSDE1zkL7vobHniMY2Ihr8BY0tn0tTPgDb3gwLvzysf1UBCTTBNr9PX9d69yEfH4kOvA4EUJqQcddCbbZVB25G5GAikbn09GxCbe2pNDdnvwyqOBlpqr7VvRILJPatisgPMu3MraY8D0BQQDK+ui9A7TGIhhFpoj7bEOf1Pgu7PgtLL4E1T0OgDpl0Emz64/TzDxo2h73mw5KfmBFx8QiM2R22vAzWz7LibOPZSMNZEGpHJnyUxxDsNmpr26jN3vI2DL8nnSx0uKiTePPhpG0XpynbltUNI6xkOoFlqinaeCqcslwXgS2h8SEE0xiXdVD0hCNg50dgyY9h7RxAkPUOgS0vh+ad09sFJ8OkF5HoO/D228hGHaZlISsCnEBt7Qm+1VUJkpGWrW/VraZ8EYAKQgQkj+DbshvsmOdNc8NmsMNf8rMZRJy8GKN5/k8i/jgNNwK7Aa85n3fEDDUcKyJnqurjnnk2mhh3AOz+PGgckNwyzA8SmgayIsfgMwpIratVGRaky6dvFWADTN9quuZtV5ryx89CtRGPw5Il0NsLy5Y5oxBGMf4YLgomTciuqrq7qu4O7AK8i5nsepUnHo1mJJBf8Fm2DB54AFavhqefhmjVPbTmR2G6KnbfqitN2QCUifeXwGXnwNG7wOc/Bd/5Asx5Kn1A6eqC3/wGpk6FnXaCt96C7baDadNg5kzoyzpXsnrxPmkiwNaqOrScrZNVYVtVfdcTb0YzkV748E1Y+TZEs0wmfeYZmD4dtt0WTjkFli6FI4+EDTaAn/8cursz21czeerKaRYb7FtdBNwz2LcqIkcV4IErTVVuE1wpGRiAn3wDHrkH4jHzGWDxG/DvJ2GDKXDr47BhwpPsBx/A/vvDRx+ZJx+AWMy8X7IEzj0XbrgBnnoK1htl+S/9kTQRYKGI/B4z8gfgy8AbIlJH8nhdS2lY8yH841L4z5+dZmo1T0JtZ8DhF0Jj0ojVG2+EH/wAepxBC319poWhy8l88stfwl13wXPPWV3lSJH7Vl1pyh8/C+Uisgw++SOsuduk1kmFKnz/BHj0r9Dfty74DNIThqXvwJf2hI8/MtvCYTjgAFi+fF3wSaa7G958Ew4+OHXTgUZh9QPw/g/gv1eYCXnVhD+a4E7FzE48x3m962wbwOaGKz0fLYFLdoTnboZID/R3QX8Y+tbCk9fApbtC16p15R96CM4/f13wSUVfn7nBO+SQ0dnU7b2uTsWFporirogcKiJvichiEflRMfaZFZ0DejPoP7NfeBqH5WfAW9vCf78Ly78Jb2wAa/46suy/n4TnHoO+NIEEzFPRmlXwmwvN5z//2Tz5xDJMeAWIRODtt+HhpFWVB/4Hr24HS06GlVfDip/Cq9vAyhsy769SGByt43ETnKr2quqvVfULzutqVe1R1biq5pkHpfR4oqt80T6I3Ah9R0D/BSb/YcpyCr87DLpXQyzFDVi0H1avgD8kJAc+77z0N3SJRCKmuftf/xr5Xf/jsGpH+F89rNoW+v6e23FVAj7QlVtNuQ5Auc6szUTeuZDil4IehMbPJh47ioHIVzMHoY9/D6tvN2KJd5vEpdoDH5wC/e8ML3vzVdCbQ5tyNAoP3gnhLvjVrzLfpSUSDsOVVw7ftuQk6F9qEjOCSbCofbD8Auh5NYdd+jwXHHh6pyYi9zh/XxORV5Nf5fMkdzzRVb42Goe+g2DgBxB7CB34HdqzA8RTPL2//Qx0/tcZ+ZaGWASWPG+elF56ybQo5O44XJXU5x15DtZ8AWKvA/0QewvtPBH6H8lxlz7PBQee6apYmipGH9DQzFrHocGZtW9ktHLIOxeSvg9cDvSZkdgCqvfT0/sEjY2HpLbpuNoEnORdaRRZdSNs/NvBDfBCmuzYqQjVwJMPwocf5m4DMGeOqUsEIiuh6zmGkiwm+hfvR/53PWyefgh/ReSC874P6Gzn7xGeepEf5dVVITbxdoi/ChhtiQwQ1zADPb+ktvmm4WVfuBP6c7hJ0zgseABejkJ/tgzySTz77PDP4Z8N+TaI0Ets7YUEJx6WcVcVkQvOW10VRVOiLttNReRLwKGq+g3n88nAXqp6VlK5GcAMgMmTJ+8+a5bps4pGo0QiEYLBILFYjNraWkKhTHGxG3gHSGju0gCx+BSCwXV5oQZnFgPQu4C0OaiCY6F2K2c/Cm/Mz/HIgUAAJk2B95ebjtEkwlOm0JzuLm633UwAivdC31vDj2dYHS1Qv3VaF3I5f8PORY42+dY1ffr0eenmHrRuJDp3xsjt8jPS2pQCEWkCelU1LiJbA9sCj6hmWj/AG8qvqwKuJe0AXU6ytlSbkOR1gTqWQM+a3A5+7IYmbqS5scuoq90T0uxEXyN1EtMghHbJ6EK2c5GsqVxsCqnH77pyq6myjYJzZuLOBGhtbdW2tjagkCegVaDHk3hnE4vV0R95hsamPYe2tbe3M1gHb58LfQtG7krqkUk/gcnfWLfxe4dBJMc7r6YxcOVf4Ec/Mu3QSbRffTVt56fIi1VTY+7uREyqnwVfNE1uyf5Ri2x4LkxNcZU55HL+hp2LHG0KrSsl3j8BDfIMsL+IjAcex8yJ+DJmLYyKpGi6ytFm2LUUewX6TmWYFuO1ROUC6prPGL7zex6E569JvdxCIrUNcNyv4b0+uOiilE9BaXXV0gJrE9KQrZlpstAnBEhViAUPIjTxnIxuZDsXyZrKxaaQejLiD1250lQxAlC2mbUZyT//1ATQOzFJBgNAhIHYL4cFnxFs8HNYdvywZjgFRGph/W8OL3vocfDgXWagQTaCQZh+mJnzMzdTDswEAgEzh2FwAl5oLKx/Enx8J+jwDlcJ1MCkb2fcXbXmgisRoqo9InI6cKOqXiUiL3vtVBrKq6tCbII7Q+gUiN6GatQk95Up1DWlSB+279fhqRuzB6B4HPY4HjZZAxdemNXndb4E4YtfTDqgn0HkIUf3UVSDIPWExv066+6qOBdcsXGlqWIEoKGZtRiBnACkyQyamvx/zI4GXY7Ie8AU6uszpWQHxnweNr4BPjzXDHfWKFK7KWxyF4SScsGddh48fm/mUXAAdQ1w8nchFIIf/hBOPTW3CXENDWZ0TyKb3WBWXV3zD2fJcIVALWx1L9RNTbmbRKo1F1wJEBHZB3N3drqzzR+ejaT8uirEpu5GCJ2KxJ9FZAsIHpF6meyNPwVbfRreec6MeEtFbQPs+RVoXt+89thjZL9OWj/qzFy7RELTYP2XofsKGHgRCe0MTT+G0DY57bKCc8GVE1each2A0mUtdrvfrMg4YNfcy693Kow/CfoWmZTwdVulLrfdLnDGhXDT5dCXptO0rh623sGUA/jCF+D3v4d//ztztoPGRjjmGNhnn+HbA3Uw7R7oXwbdL0JwPRhzYPY1iSoJ7+/UwHSc/hi435n9vQWQx6iT8uGZrgohuKd5ZePMe+Hq6Sb7QSTpZq2uCbbaD75647pt110H++6b/causdE8/ey0UwrfNocxf8juW6Xiva5caaoov3CpZtb6EqmBhhQXaTJn/gTGTYCrf2RG5XQ7k1brGgCFQ46FX/wf1DpLMgSD8I9/mED0/PNmSHbi4I5g0NyhHXss3Hpr+vxXdZuaVzXivVBQ1WcwbdaDn98FvuedR5mpGF3lSuNYuHAOvHQPPHoVrHzLaGHT3eDQC2Dno0wT9SA772zmzB1xhJkQnurmrqkJDj/c6Go04rGu3Gqqim6xi8yJZ8Cxp8GTD8D858zAhM23gaNPhvVSNPk1NsKjj5oA9KtfwRNPmO1NTeap59xzzci30Yg/OkstfiBUC/t81bxy4YADzCTTG280qaz6+02Qqq2F/faDCy4wWRDySWpaLVSBrmwAykRtLRx6LLTuaCbJrbcF1GVYAkLEiGLvPWHx0/BuGBY8BFvuBwHvG2s9ZZQfvsUFG25oko7+9KewciUsXGgyj4xNvdLxqKLCdWUDUDriMXj2anjmaoj2mYSJsQjseBwcdhW0bJDCJg5PXAFPXgUo7HwJ/OEUCNXBEb+A/b5V7qPwB1Vwp2YpIp2LYeE1sOJJk5B002Ng+29D00aZ7UIhmDIFFi+2wQeqQlejJwCpQt+/oHs2BMZA88lQm2aCZzwGtx8DS/4FA0kDEV65C955HM6aB2M3Hr7/u74J82eZRItg+o/6u8zr/nNh7Ydw+CVp6gxD523Q/QiENoJxZ0B9FTXZeSgUEbmOwbXYU6Cqvu0HqghiH0J0LgQ3hVCWPtbFd8Kz3zCjUePOXMW1TkA65O+w0WfS2/a/CmtnwsBusPY9aD4BAqN8YTqPdFUsTVV4/MwRjcH/vgArj4a118KaK2DFLuZiTsWCv8CSp0YGHzDzGLpXwb2nD9++5FmYf/e64JNMpMc8GX30zsjvoqvgvR2h4wLofgg6b4H394PVN44sO8Kf2yG2E8Q+BfEb/JkR2PukiXOBeUA9ZvXGd5zXLkCeCzVbhtF9EXyyOXR9FVbvA6v3N7kWU7F6ETz7TYj1rgs+APF+iIbh8aOhb1Vq284/wPK9ofMmiK+Gju/CB7tBrLP4x1QpeKuromjKFwGo5Ikxu/8GvU+CDg7njJpJnx+fDbGPRpZ/+goYyDD0U2Pw3tPQmZAO5J+/Sh98BolHof3akdtXXQzRFQkTZePGv47zTHBKu79ZoGdgVsN9A41fAHpdZh8cRlMyUlW9TVVvA3YC2lT1OlW9DjgII5iqpOSJMSOPQ881QD/oWqAHjb4I3RekLv/ar02wSYfG4M2bR26PfQyrznEmajsTxLUbBt6D1Zdn9xN/JxatxGSkxdKU5wFoMBXF8uXLWbRoUc7/hLzsum5LCD7rUILQ89DwjbEorHoruwPBWvjgxXWfl75AhidSZ98DsDjFEPnwX0m1dpMShO5H0+9Pf0diGhSRHuKxazL7QJnOeSL+WZJ7PDAm4XOzs63qKOR/lbdN7x8xuRnXIUTQvjtTl1/xuAky6Yj1wvv/GLm951FS9xb0QzhNXQmU5VyU0WYIf+jKlaZ8EYBEhJaWFkQkrx/D3O1SD9E0rVUuhm8WNPQzlU0m//LbVzyevQmuPOc8CR+sBwRcASwQkT+JyG3AfOCXnnhSYgr5X+VtI6l/PtInOM5BLyk1lcku+09YWc5FGW2G4b2uXGnK8wDU3NyMqtLV1YWq5pWXLGe7llNBmkZsFmLQmJRNPBiCidtldyAWgal7rfu8+T7ZA1KwBrZO0cna8mVgZPoSIQbNGdLGy/dRXdcJG4vVE42flb68Q1nO+TA/8cOdGqr6R2Av4H7gPmAfpxmh6ijkf5W3Tf3XgcZhm+LxWqKhE1KXn3Jo5uwewQbY5KiR2xsPI+VyJdRBc/acl2U5F2W0GcIHunKrKc9HwZUlMWbTsRC+B3ofQbUb82MfRCbcAMEJI8sf+CN44Mz0/UCBEGwxHcYkDBs96Afw1j9HphhJtjswxeCQCT+D7och+j/QbtP0JrXIpN9BcP0M+zseQYnFfkk8PkAsfib19bkFoNGYjFREBPgssIWqXioim4jInqr6YjbbSqMsiTFrPwuNP4CeK1FqQQeIh1qpHf/b1OV3PBcW3556RVQwQ7K3OX3k9uB4mHADrDqLwSz/Kk1IaDNYL/tCsX5OLFrpyUjdasrzAARlSIwpAZh8D/Q9i/Q8aNbYaf4K1GyZuvwuJ8HC++Gdx0aOhAvUQNMEOPaW4du33A/2OAleuiN1EKpthIN/DBNT1BlcDzZ7DTr/Aj2PIcENYdwMqM8hbVDgywQDXyZIqmeo9IzSZKQ3YnLzfwa4FOgC7gX28NKpUlGWxJhNl0DDt5HofAhsSiiUofVg3DZw4G3w9CmmLyjuLGESrDdPRp97COrT3HCNPQ0a9kY6b4bgesikm6D5OJC60hyXz20Av+jKlaZ8EYDKggg0HGBe2QgE4KS/mvVLnvkV9K81mQziUROcDrkMmieNtPvyTTBxGjz+SzOXSIImc0JdExz5S9j76xnqbITxM8yrGvHBExBmQbfdRGQBgKquFhE7DNstgUlQe2huZbc4DibuAQuvhxWPmcCz2Rdh2xnQODmzbe32MPE3EGqHlja3XlcH3uvKlaZGTwDKl0AQ9j8X9j0HPnnXzFsYtwnUjuxLGkIEDjof2s6Bd5+DxWvg24/AZvsMT7I42ihwxraIHApcg7nPu1lVr0j6/lzgG5gOgg7gNFVdlmGXAyISxBmuKCITSbtUrqVktGwGe18NXO21J5WNPzIhuNKU9+77nUAAJmwFk7bLHHwSCYZgWhs0jIMt9h3dwWeQPDtLnYv6BuAwYHvgRBHZPqnYAqBVVXcC/gZclcWLazGdpZNE5DLgOap0FJxllOD94B5XmrJPQJbSU1hb9Z7AYie9OyIyCzgaeGOwgKomTqqaA2RMsayqd4jIPMxkOQGOUdVFeXtmsfiBAvuAitmy4FZT9tbcUh5S36lNEJG5Ca/EDrCNgQ8SPi93tqXjdOCRTC6IyC1AvareoKrXq+oiEbkk72OxWPyCxy0LbjXlKgCJyK9E5E0ReVVE7heRcW7250s0CvHnIf4E6Ce52UQXwZqTIfoGdJ4JsUzdEoP19ELPtbDmM7D2azCwIEf/+iH+HMSfTpntwTekFsoqVW1NeKVJzpcZEfkq0Ar8KkvRzwG3icjXEralmHjiLaNCV5bikH8T3FDLgqpGgMGWhSFU9SnVobxgc4ApGfbnSlNun4CeAHZwIuXbmKVZqwfthujeEDsUYsdBdEvQ+ZltBhbAx3tA311AL/TeDKt2geh7GeoZgDX7QfePYOAp6L/DfO7PeEMP2gHRHSB2OMSOgug2kLEP3iMKS5q4Apia8HmKs234rkU+C1wEHKWqGRKNAfARcABwnIjcICIhXKXCKBnVrStLcUivq3K2LLjSlKsApKqPq+rgrLJskTIt5UyMmZdN7CJgIRAGOoE1ED06owldFzhPIoM5r6KgXdD9i/Q2kb9D9G2g19kQB3ognGVSaWwGsAwz9H4tqishlttKkxWQjPQlYJqIbO4M6zwBmJ1YQER2Bf6ACT4pssqOQFS1U1WPxLRttwO+W1jGS135VosJNtFo1Nf++SAZaTlbFlxpqpiDEE4D7s7XaDAZn4igqmy33XY5TcoqxC5vG30JGL4OveoKRAdA0kz7HEj1hBSDyL/T1zPwAibIJVUffy9zXfoyiUlMRWJo/LWstx/lPOfGMfK+1VHVqIicBTyGua+7VVUXisilwFxVnY0RRjPwVzMhm/dVNdPj/1AAU9VLnM7T7+fnWdkpm658rcUEm0gkwqJFi3zrXznOOVDoMOx8WxYOzNKy4EpTkj5x4JAjTwIplv/kIlX9u1PmIkyk/KKm2aHzGDgDYPLkybvPmjULgGg0SiQSIRgMEovFqK2tJRTKHhez2Q2mtsjHZiRLnX6fxEMKguySwbFFDGaoDvdOobnBWbJBxkEwTeYFXQWxDxg5fD4EoZ0z+Pe2ebpK3BX1iHxq2Lbkc1GKcz59+vR5qtqayq51a9G5KZY2koNJa1Pt+FFXudgUR1fubKLRKKFQqGi/FYXYlOs8FFtXThPZ25hRayswLQ1fUdWFCWV2xQw+OFRVUyxgVjyy/vdU9bOZvheRU4EjgIPSicTZz0xgJkBra6u2tbUBpbvram9vZ7COXG1GOv0h8YFd0XgnInEUIRK7nYbGtvQ2/RF09TEIvbS/fjVtO5yP0oCs/zzU7Jqmnm74ZGs0/hHiJF2MxesZqL2K+vEZ6tIN0YG9iMf7zN2QBuiPPUpj0/BsD8nnohKegIqJiDynqvuJSBfD7yYEUFUdk8a0ZPhRV7nYFEVXLm06OjqYOHGip09AXp6HITxsWSiWplw1wTnjyS/APKZlWY0tNeVMjJm3jWxIoGYRfb13EImsIVTzORqbsqQ4qjsEGX8vsc4LgACxwB4Ex/0uffABk6l73Fyk+yLi/Y8Q0/WI1v+IhrEnZ/FvG6TmTQZ676G/v5dQzZE0NSePqBzJaEtGqqr7OX9bvPMid7zSla+1mGDT2dmZ1w+1nxOLepGMVFUfBh5O2nZxwvuMN0dOmaJoym0f0PVAHfCEEynnqOoZ+e6knIkx87aR8dQ3nkV9Y/aiQ9QdRnDSYRBqJzgpx+bQ4IYw5tahfsScE4vKBtQ3fi8//xhdyUhFZL1M36vmOr6+bHimK19r0bEJhUK+TRI6WpKRFktTrgKQqm7lxt4yivB2yvM8TDNBqvEZCmxRXncyY3VlyRnvdFUUTdlUPJbS43EfkKpu7l3tFkuJ8FBXxdKUDUCW8uD9uiUAiMh4YBpQP7hNVZ/xziOLxQU+0JUbTdkAZCk9/kgbj4h8AzgbM/fhZWBv4D+YxbQslsrCB7pyqykf/CxYRgXep40HI5Q9gGWqOh3YFVjjiScWSzHwXleuNGWfgCzlwR+3On2q2iciiEidqr4pItt47ZTFUjDe68qVpnwRgMLhcEHj4AuxK6fNYM6qajmmgu38sXY9wHIns/QDmCHOqzHJ9KoSv18XVleF2wB+0ZUrTXkegCol/9Roz1nlxg7ww50aqvoF5+0lIvIUJmniox66VDL8fl1YXRVuMwyPdeVWU57/LITDYUSElpYWRCTnjLCF2JXbJhgMVs0xubEb6iz1vg8IERkvIjthUogvB3bwxpPS4vfrwuqqcJshfKIrN5ry/AmoubkZVaWrqwtVzSu9Rr525baJxWJVc0xu7AA/NBUgIj8HTgXeZV3mV6UKR8H5/bqwuircZhge68qtpnwRgCoh/9Roz1nlxs4Pw0Udjge2VLMSZFXj9+vC6qpwmyH8oStXmvI8AEFl5J8qxKbacla5sfOBUABeB8ZhVnGsevx+XVhdFW4zhPe6cqUpXwQgS5Xjj9E6AJcDC0TkdWBoka3kVPMWS0XgD1250pQNQJby4P2dGsBtwJXAa4xc/c9iqTy815UrTdkAZCk9/mirBuhR1Wu9dsJiKQr+0JUrTdkAZCkP3jcVADwrIpdj1rFPbC6Y751LFosLvNeVK03ZAGQpPf64UwOTpwpMwsRBqnIYtmUU4A9dudJUUQKQiJwHXA1MVNVVxdinpcrwPmtvEJitqr/11pPcsbqyZMVDXRVDU67dF5GpwCHA+273ZaliPJ6xraox4MTy1lo4VleWnPBQV8XQVDHc/S1wAeaxqyDC4TArV67MLw1FgXbltBlMmpgPfj6mgu0Gh4smv8rP8yJyvYjsLyK7Db488SQ7nuiqErRYbboq9Jz7RFeuNOWqCU5EjgZWqOorIqmWBs9OpSRAHO1JE93Y+aStGmAX5++lCdt81wfkla4qRYvVpCtXyUj9oatdnL8FaSprABKRJ4ENUnx1EXAhppkgKyIyA5gBMHnyZNrb2wGIRqNEIhGCwSCxWIzOzk5CoexxMZtdOBweqiNXm0LqyWQTjUbp6Ogo2jEVapN8Lkp1zjPivVBwFszyBX7UVSHXUinrSmfjB115eR6G4X02bFeaynqkqvrZVNtFZEdgc2DwLm0KMF9E9lTVlSn2MxOYCdDa2qptbW1A6e662tvbGawjV5tC6slk09HRwcSJEz2/U0s+F548AXk/XBQRGQv8FDjA2fQ0cKmqdpbbFz/qqpBrqZR1pbPxg668PA9DFKgrETkUuMaxvllVr0j6vg74M7A78DHwZVVdmmZfrjRVcBOcqr4GTEpwZCnQmu9onUpJgDjakya6sQM8v1NzuBWTu+p45/PJwB+BL3rmURJe6qpStFhNunKlKchbV87ItRuAgzFLJ7wkIrNV9Y2EYqcDq1V1KxE5AZPp4MtpdulKU76YB1QJCRBt0kR3dj4JQFuq6rEJn38mIi975Uyp8ft1YXVVuM0Q+etqT2Cxqr4LICKzgKOBxAB0NHCJ8/5vwPUiIqqaakCMK00V7WdBVTezcxUsKfHHaB2AXhHZb8gtkX2BXk88yRGrK0taCtPVxsAHCZ+XO9tSllHVKNAJrJ9mf6405ckT0Lx581aJSM7rhhfIBMBr4frBByiPH5um+2LePB6TABNSfFXuc3MG8Gen3VqATzCLaVUFVldV6UMhuqoXkbkJn2c6fYWlwJWmPAlAqjqx1HWIyFxVbS11PX73wQ9+qOqhXtWdiKq+AuwsImOcz2s9dqmoWF2NLh8K1NUKYGrC5ynOtlRllotICBiLGYyQygdXmvJFH5DFUg6c0T3HApsBocE5Nqp6aQYzi6WaeAmYJiKbYwLNCcBXksrMBk4B/gN8CfhXmv4f15qyAcgymvg7pj17HgmZey2W0YKqRkXkLOAxTI/Rraq6UEQuBeaq6mzgFuB2EVmMaVI7IcMuXWmqmgNQqdo888EPPoB//PCaKX5pDqxg/HAtWR9coKoPAw8nbbs44X0fcFyOu3OlKUnzZGWxVB0iMhO4zplrY7FYXOJWUzYAWUYNIvIGsBXwHqa5QABV1Z08dcxiqVDcaqpqApCI/Ao4EogAS4Cvq+qaFOWWAl1ADIgWYyRLMVNbFFj/VGf/kzGJAGeq6jVJZdow7bXvOZvuG22d7yKSckirqpZ66HLFYnVldZUJ15pS1ap4YZI3hpz3VwJXpim3FJhQxHqDGGFuAdQCrwDbJ5X5NnCT8/4E4O4iH/uGwG7O+xbg7RQ+tAEPev1/sq/KelldWV2V8uWPBClFQFUfVzNrF2AOZnx7ORhKbaGqEWAwtUUiRwO3Oe//BhwkhebZT4GqfqjOGuyq2gUsYuTsZoslb6yurK5KSdUEoCROAx5J850Cj4vIPCeVvVuKndrCFSKyGWad9hdSfL2PiLwiIo+IyKdKUb+lqrG6sroqKhU1DDvTGiqq+nenzEVAFLgjzW72U9UVIjIJeEJE3lTVZ0rjcXkRkWbgXuAcHTkjeT6wqaqGReRw4AFgWpldtPgQq6vMWF2VjooKQJpmDZVBRORU4AjgIHUaaFPsY4Xz9yMRuR/zqO9GKEVNbVEoIlKDEckdqnpf8veJwlHVh0XkRhGZoDbR5ajH6io9VlelpWqa4JwRMxcAR6lqT5oyTSLSMvge08H6usuqh1JbiEgtpjN0dlKZwdQWkCW1RSE47d63AItU9Tdpymww2D4uInti/vdFFaul+rC6sroqJRX1BJSF64E6zOM/wBxVPUNENsIM4TwcM5zyfuf7EHCnqj7qplItfmqLQtgXsxDUawlrcVwIbOL4eBNGoGeKSBSTLv2EYorVUrVYXVldlYyqmQdksVgslsqiaprgLBaLxVJZ2ABksVgsFk+wAchisVgsnmADkMVisVg8wQYgi8VisXiCDUAWi8Vi8QQbgCwWi8XiCTYAWSwWi8UT/j86IiCudEkstQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i,j in enumerate(random_samples):\n",
    "    plt.subplot(2,2,i+1)\n",
    "    footprint=S[j,...]\n",
    "    xd, yd = rectangular_array()\n",
    "    mask = footprint != 0\n",
    "    mask[5, 5] = True\n",
    "    marker_size = 200 * footprint[mask] + 20\n",
    "    plot = plt.scatter(xd, yd, c='grey', s=10, alpha=0.3, label=\"silent\")\n",
    "    circles = plt.scatter(xd[mask], yd[mask], c=footprint[mask], s=marker_size,\n",
    "                          cmap=\"autumn_r\", alpha=1, label=\"loud\")\n",
    "    cbar = plt.colorbar(circles)\n",
    "    cbar.set_label('normalized signals [a.u.]')\n",
    "    plt.grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.stack([T, S], axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis = f['showeraxis']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "core = f['showercore'][:, 0:2]\n",
    "core /= 750"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# energy - log10(E/eV) in range [18.5, 20]\n",
    "logE = f['logE']\n",
    "logE -= 19.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test = np.split(X, [-20000])\n",
    "axis_train, axis_test = np.split(axis, [-20000])\n",
    "core_train, core_test = np.split(core, [-20000])\n",
    "logE_train, logE_test = np.split(logE, [-20000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Model\n",
    "## Task \n",
    "Set up a multi-task regression network for simultaneously predicting shower direction, shower core position, and energy. The network should consist of a common part to the three properties, followed by an individual subnetwork for each property."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input1 = layers.Input(shape=(9, 9, 2))\n",
    "\n",
    "# TODO: define a suitable network consisting of 2 parts:\n",
    "# 1) a common network part (you can try a convolutional stack with ResNet- or\n",
    "#    DenseNet-like shortcuts)\n",
    "#   z = ...\n",
    "# 2) separate network parts for the individual objectives\n",
    "#   z1 = ...\n",
    "#   z2 = ...\n",
    "#   z3 = ...\n",
    "\n",
    "output1 = layers.Dense(3, name='direction')(z1)\n",
    "output2 = layers.Dense(2, name='core')(z2)\n",
    "output3 = layers.Dense(1, name='energy')(z3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task\n",
    "Train the model to the following precisions:\n",
    "- Better than $1.5^\\circ$ angular resolution\n",
    "- Better than $25$ m core position resolution\n",
    "- Better than $10\\%$ relative energy uncertainty $\\left(\\frac{E_\\mathrm{pred} - E_\\mathrm{true}}{E_\\mathrm{true}}\\right)$  \n",
    "  Estimate what these requirements mean in terms of the mean squared error loss and adjust the relative weights in the objective function accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_weights=[1, 1, 1]  # you can give more weight to individual objectives\n",
    "\n",
    "model.compile(\n",
    "    loss=['mse', 'mse', 'mse'],\n",
    "    loss_weights=loss_weights,\n",
    "    optimizer=keras.optimizers.Adam(lr=1e-3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.callbacks import ReduceLROnPlateau, EarlyStopping\n",
    "\n",
    "fit = model.fit(\n",
    "    X_train,\n",
    "    [axis_train, core_train, logE_train],\n",
    "    batch_size=128,\n",
    "    epochs=40,\n",
    "    verbose=2,\n",
    "    validation_split=0.1,\n",
    "    callbacks=[ReduceLROnPlateau(factor=0.67, patience=10, verbose=1),\n",
    "               EarlyStopping(patience=20, verbose=1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot training curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_history(history, weighted=False):\n",
    "    fig, ax = plt.subplots(1)\n",
    "    n = np.arange(len(history['loss']))\n",
    "    for i, s in enumerate(['axis', 'core', 'energy']):\n",
    "        w = loss_weights[i] if weighted else 1\n",
    "        l1 = w * np.array(history['%s_loss' % s])\n",
    "        l2 = w * np.array(history['val_%s_loss' % s])\n",
    "        color = 'C%i' % i\n",
    "        ax.plot(n, l1, c=color, ls='--')\n",
    "        ax.plot(n, l2, c=color, label=s)\n",
    "\n",
    "    ax.plot(n, history['loss'], ls='--', c='k')\n",
    "    ax.plot(n, history['val_loss'], label='sum', c='k')\n",
    "\n",
    "    ax.set_xlabel('Epoch')\n",
    "    ax.set_ylabel('Weighted Loss' if weighted else 'Loss')\n",
    "    ax.legend()\n",
    "    ax.semilogy()\n",
    "    ax.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Unweighted losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_history(fit.history, weighted=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Weighted losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_history(fit.history, weighted=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Study performance of the DNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis_pred, core_pred, logE_pred = model.predict(X_test, batch_size=128, verbose=1)\n",
    "logE_pred = logE_pred[:,0]  # remove keras dummy axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction performance of the shower axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.sum(axis_pred * axis_test, axis=1) / np.sum(axis_pred**2, axis=1)**.5\n",
    "d = np.arccos(np.clip(d, 0, 1)) * 180 / np.pi\n",
    "reso = np.percentile(d, 68)\n",
    "plt.figure()\n",
    "plt.hist(d, bins=np.linspace(0, 3, 41))\n",
    "plt.axvline(reso, color='C1')\n",
    "plt.text(0.95, 0.95, '$\\sigma_{68} = %.2f^\\circ$' % reso, ha='right', va='top', transform=plt.gca().transAxes)\n",
    "plt.xlabel(r'$\\Delta \\alpha$ [deg]')\n",
    "plt.ylabel('#')\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction performance of the shower core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.sum((750 * (core_test - core_pred))**2, axis=1)**.5\n",
    "reso = np.percentile(d, 68)\n",
    "plt.figure()\n",
    "plt.hist(d, bins=np.linspace(0, 40, 41))\n",
    "plt.axvline(reso, color='C1')\n",
    "plt.text(0.95, 0.95, '$\\sigma_{68} = %.2f m$' % reso, ha='right', va='top', transform=plt.gca().transAxes)\n",
    "plt.xlabel('$\\Delta r$ [m]')\n",
    "plt.ylabel('#')\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction performance of the shower energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 10**(logE_pred - logE_test) - 1\n",
    "reso = np.std(d)\n",
    "plt.figure()\n",
    "plt.hist(d, bins=np.linspace(-0.3, 0.3, 41))\n",
    "plt.xlabel('($E_\\mathrm{rec} - E_\\mathrm{true}) / E_\\mathrm{true}$')\n",
    "plt.ylabel('#')\n",
    "plt.text(0.95, 0.95, '$\\sigma_\\mathrm{rel} = %.3f$' % reso, ha='right', va='top', transform=plt.gca().transAxes)\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(19.25 + logE_test, 19.25 + logE_pred, s=2, alpha=0.3)\n",
    "plt.plot([18.5, 20], [18.5, 20], color='black')\n",
    "plt.xlabel('$\\log_{10}(E_\\mathrm{true}/\\mathrm{eV})$')\n",
    "plt.ylabel('$\\log_{10}(E_\\mathrm{rec}/\\mathrm{eV})$')\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}