{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*PyPlr* and Pupil Core\n",
    "==========================\n",
    "\n",
    "<img src=\"../img/Pupil-Labs-Glasses.png\" alt=\"drawing\" width=\"500\"/>\n",
    "\n",
    "*PyPlr* was designed to work with [Pupil Core](https://pupil-labs.com/products/core/)—an affordable, open-source, versatile, research-grade eye tracking platform with high sampling rates, precise model-based 3D estimation of pupil size, and many other features which make it well-suited to our application (see [Kassner et al., 2014](https://arxiv.org/abs/1405.0006), for a detailed overview of the system). In particular, we leverage real-time data streaming with the forward facing scene camera to timestamp the onset of light stimuli with good temporal accuracy, opening the door to integration with virtually any light source given a suitable geometry.\n",
    "\n",
    "The best place to start learning more about Pupil Core is on the [Pupil Labs website](https://docs.pupil-labs.com/core/), but the features most relevant to *PyPlr* are:\n",
    "\n",
    "* [Pupil Capture](https://docs.pupil-labs.com/core/software/pupil-capture/#world-window): Software for interfacing with a Pupil Core headset\n",
    "* [Pupil Player](https://docs.pupil-labs.com/core/software/pupil-player/#load-a-recording): Software for visualising and exporting data\n",
    "* [Pupil Core Network API](https://docs.pupil-labs.com/developer/core/network-api/): Fast and reliable real-time communication and data streaming with [ZeroMQ](https://zeromq.org/), an open source universal messaging library, and [MessagePack](https://msgpack.org/index.html), a binary format for computer data interchange\n",
    "\n",
    "*PyPlr*'s `pupil.py` module\n",
    "---------------------------\n",
    "\n",
    "*PyPlr* has a module called `pupil.py` which facilitates working with the Pupil Core Network API by wrapping all of the tricky ZeroMQ and MessagePack stuff into a single device class. The `PupilCore()` device class has a `.command(...)` method giving convenient access to all of the commands available via [Pupil Remote](https://docs.pupil-labs.com/developer/core/network-api/#pupil-remote). With Pupil Capture already running, we can make a short recording as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'OK'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from time import sleep\n",
    "from pyplr.pupil import PupilCore\n",
    "\n",
    "p = PupilCore()\n",
    "\n",
    "p.command('R our_recording')\n",
    "\n",
    "sleep(10)\n",
    "\n",
    "p.command('r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Annotations and notifications\n",
    "-----------------------------\n",
    "\n",
    "To extract experimental events and calculate time-critical PLR parameters (e.g., constriction latency, time-to-peak constriction) we need a reliable indication in the pupil data of the time at which a light stimulus was administered. The Pupil Labs [Annotation Capture](https://docs.pupil-labs.com/core/software/pupil-capture/#annotations) plugin helps us with this. The plugin allows timestamps to be marked with a label manually via keypress or programmatically via the Network API in a process that is analogous to sending a 'trigger', 'message', or 'event marker'. With *PyPlr*'s `PupilCore()` interface, annotations can be generated programmatically with `.new_annotation(...)` and sent with `.send_annotation(...)`. It is important to make sure that the Annotation Capture plugin has been enabled. You can do this manually in the Pupil Capture GUI or programmatically by sending a [notification message](https://docs.pupil-labs.com/developer/core/network-api/#notification-message), which is a special kind of message that the Pupil software uses to coordinate all activities. The following example shows how to enable the Annotation Capture plugin programmatically and then send an annotation with the label `'my_event'` halfway through a 10 second recording."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'OK'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = PupilCore()\n",
    "\n",
    "p.annotation_capture_plugin(should='start') \n",
    "\n",
    "p.command('R our_recording')\n",
    "\n",
    "sleep(5.)\n",
    "\n",
    "annotation = p.new_annotation(\n",
    "    label='my_event', \n",
    "    custom_fields={'whatever':'info','you':'want'})\n",
    "\n",
    "p.send_annotation(annotation)\n",
    "\n",
    "sleep(5.)\n",
    "\n",
    "p.command('r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the recording is finished, we can open it with Pupil Player and use the [Annotation Player](https://docs.pupil-labs.com/core/software/pupil-player/#annotation-export) plugin to view and export the annotations to CSV format. Any custom labels assigned to the annotation will be included as a column in the exported CSV file. By default, the timestamp of an annotation made with the `.new_annotation(...)` method is set with `.get_corrected_pupil_time(...)`, which gives the current pupil time (corrected for transmission delay); but this can be overridden at a later point if desired. \n",
    "\n",
    "Notifications can be used for many things, but there is no single exhaustive document. One way to find out what you can manipulate with a notification is to open [the codebase](https://github.com/pupil-labs/pupil) and search for `.notify_all(` and `def on_notify(`. Alternatively, if you just want to access the pupil detector properties, there is a handy method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'subject': 'pupil_detector.properties.1.Detector2DPlugin',\n",
      " 'topic': 'notify.pupil_detector.properties.1.Detector2DPlugin',\n",
      " 'values': {'blur_size': 5,\n",
      "            'canny_aperture': 5,\n",
      "            'canny_ration': 2,\n",
      "            'canny_treshold': 160,\n",
      "            'coarse_detection': True,\n",
      "            'coarse_filter_max': 280,\n",
      "            'coarse_filter_min': 128,\n",
      "            'contour_size_min': 5,\n",
      "            'ellipse_roundness_ratio': 0.1,\n",
      "            'ellipse_true_support_min_dist': 2.5,\n",
      "            'final_perimeter_ratio_range_max': 1.2,\n",
      "            'final_perimeter_ratio_range_min': 0.6,\n",
      "            'initial_ellipse_fit_treshhold': 1.8,\n",
      "            'intensity_range': 23,\n",
      "            'pupil_size_max': 100,\n",
      "            'pupil_size_min': 10,\n",
      "            'strong_area_ratio_range_max': 1.1,\n",
      "            'strong_area_ratio_range_min': 0.6,\n",
      "            'strong_perimeter_ratio_range_max': 1.1,\n",
      "            'strong_perimeter_ratio_range_min': 0.8,\n",
      "            'support_pixel_ratio_exponent': 2.0}}\n"
     ]
    }
   ],
   "source": [
    "from pprint import pprint\n",
    "\n",
    "p = PupilCore()\n",
    "\n",
    "properties = p.get_pupil_detector_properties(\n",
    "    detector_name='Detector2DPlugin', \n",
    "    eye_id=1)\n",
    "pprint(properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows some of the pupil detector properties that can be modified with notifications. Bare in mind that for most use cases it will be best to verify manually in Pupil Capture that all of your most important settings are as they should be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting pupil data in real-time\n",
    "-------------------------------\n",
    "\n",
    "The Pupil Capture continuously generates data from the camera frames it receives from a Pupil Core headset and makes them available via the [IPC backbone](https://docs.pupil-labs.com/developer/core/network-api/#reading-from-the-ipc-backbone). `PupilCore()` has a `.pupil_grabber(...)` method which simplifies access to these data, empowering users to design lightweight applications that bypass the record-load-export routine of the Pupil Player software. Just specify the topic of interest and how long you want to spend grabbing data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Grabbing 10 seconds of pupil.1.3d\n",
      "PupilGrabber done grabbing 10 seconds of pupil.1.3d\n"
     ]
    }
   ],
   "source": [
    "p = PupilCore()\n",
    "pgr_future = p.pupil_grabber(topic='pupil.1.3d', seconds=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of note, `.pupil_grabber(...)` does its work in a thread using Python's `concurrent.futures` framework which means the grabbed data can be accessed via a call to the `.result()` method of a returned `Future` object once the work is done:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'id': 1,\n",
       " 'topic': 'pupil.1.3d',\n",
       " 'method': 'pye3d 0.0.6 real-time',\n",
       " 'norm_pos': [0.5360164625892022, 0.6350228570342871],\n",
       " 'diameter': 32.04137148235498,\n",
       " 'confidence': 1.0,\n",
       " 'timestamp': 55693.943421,\n",
       " 'sphere': {'center': [6.045138284134944,\n",
       "   -0.6276101187813217,\n",
       "   49.29280081519474],\n",
       "  'radius': 10.392304845413264},\n",
       " 'projected_sphere': {'center': [130.62958410928272, 92.31126323451176],\n",
       "  'axes': [142.37459004874788, 142.37459004874788],\n",
       "  'angle': 0.0},\n",
       " 'circle_3d': {'center': [0.7818853071011702,\n",
       "   -3.159317419301394,\n",
       "   40.696951453773934],\n",
       "  'normal': [-0.5064567538505914, -0.24361364857743406, -0.8271359904549626],\n",
       "  'radius': 2.0321335156335523},\n",
       " 'diameter_3d': 3.6276776361286016,\n",
       " 'ellipse': {'center': [102.91516081712683, 70.07561144941687],\n",
       "  'axes': [27.30486055019297, 32.04137148235498],\n",
       "  'angle': 32.95347595832888},\n",
       " 'location': [102.91516081712683, 70.07561144941687],\n",
       " 'model_confidence': 1.0,\n",
       " 'theta': 1.8168863457712132,\n",
       " 'phi': -2.120212108874702}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pgr_future.result()\n",
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data are returned as a list of dictionaries, with each dictionary representing a single [data point](https://docs.pupil-labs.com/developer/core/overview/#timing-data-conventions). With the `unpack_data_pandas(...)` helper function from `pyplr.utils`, we can organise the whole lot and inspect the pupil timecourse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAEGCAYAAABM/fUaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABmvElEQVR4nO3dd3xUVfrH8c9J770QCAFC6L2DAgKiYu+9rGtbXbtuc4tr2eLqrvpz19V1Xctad21r74IIUqQX6T1Aeu/t/P6YySQhhQCZTCb5vl+vvJi599x7n4EwM8895zzHWGsRERERERHpCXw8HYCIiIiIiEhnUQIkIiIiIiI9hhIgERERERHpMZQAiYiIiIhIj6EESEREREREegw/TwdwpOLi4mz//v09HYaIiIiIiHRRK1euzLHWxre0z+sSoP79+7NixQpPhyEiIiIiIl2UMWZPa/s0BE5ERERERHoMJUAiIiIiItJjKAESEREREZEeQwmQiIiIiIj0GEqARERERESkx1ACJCIiIiIiPYYSIBERERER6TGUAImIiHSgTzdmkFFY4ekwRESkFUqAREREOkh5VS0/emkllz+7FIDvdufx9Nc7qKiu5ZmFO6iqqfNwhCIi4ufpAERERLqy5xfv4v73v+f7B04h0M8XXx9DbZ1l8fYcxveLZu2+AkYnR5KeX46fjwFgR3Yp1bV1XPj0EgC2Zhbz9qr9BAf4ceXUfvzfF9vwMTAoMYz8smounZxCUUU1VTV1xIUFevLlioh0e0qARERE2vDk/O0ALN6ey/X/XsG10wfwr0W7AOgVEURGUcNwtxMGx7seX/viCtfjt1ftB6C4opqckkoe+2Jrk2ucN74Pxz/0FcUVNex+6HS3vRYREdEQOBERkTYZ4+jVeXeNI4mpT36AJskPwNdbs12PFzofj0+Jcm2rrrF8ujGj2TXS88sprqhpsu2DdQf458Kdxxa8iIg0owRIRESkDXV1FoAFW7Jb3H/uuD5tHv+jEwa6HmcWV7B0Zx5RIf5Nhrqd/NhC1+Pq2jq+3prNLa+u5vcfbTqW0EVEpAUaAiciItLIV5szufXV1YxOjiIlJoTc0ioASiqb9tCcPbY37645wPCkCGYPTWDh1mzeXJmOMfCr04bRJyqYsqpapgyIcR1zsKCcjQeKmDEontlD4rnrv2sBqHUmWQAZhRX84LnlrufWWlcvlIiIHDslQCIiIo38+p0NlFbVsmRnLkt25gIQHuhHcaMEKDzIjz+cO4raOstJwxPpHxdKVLA/b65MZ1SfSK6bkdrknIMSwtiWVcJ8Zy/SiUMTmhU7+OkpQ3jk0y3MeHh+k+3FlTVEBPm746WKiPRIGgInIiLSSHRoQLNtp47qBcA5Y3sDcP2MVEID/fjbZePpHxcKwIR+0cwYFMfDF4xudvznd53Al3efQFJkEH2igjl1VC8igx1JzZDEcL5/4BSmDYxtMZ6c4soOeV0iIuKgHiAREZFGQgMaPhp/MK0fX2/NZlhSBAAzB8fzl4vG4tPCiLTQQD9eunZKq+cdGB/G4p/PodZa/H19GNknggfPHsEZo3sTEuBHfKMeoZF9ItiwvwiAAwUVpMaHddCrExERJUAiIiKN5JY6elyeuXICJ49w9PxU19aRFBnEKSN6HdN8HB8fgw+O440xXDmtv2tf35gQnrxsPCcMieepBdtdCdAV/1rGK9dNIa+0ijNGJ2k+kIjIMVICJCIiAtTU1uHn60NmUSU/PL6/K/kB8Pf1Yd7IJLfHcPpoxzV8fZqOUL/82WUA1FnL2WPbrjonIiJt0xwgERHp8ZbvyiPtVx/z2cYMSipr6BMV7NF4RvZ2DLn7xalDm2x/bfleDhaWeyIkEZFuw+0JkDHG1xiz2hjzQRttzjfGWGPMRHfHIyIiPYO1li83ZVJRXXvYtl9uzgTghpdWAo6CBp508ohefHrHTG48YSBXH9cfgFlD4lm6M49pf/yK15fv9Wh8IiLerDN6gG4HWl3JzRgT7myzrBNiERGRbuzdNft59LMtAKzam8+1L67gqkZr6rRkT24pe3LKmmwb1SfSbTG215Be4QD89szhLPvliZwwON61769fbae6ts5ToYmIeDW3zgEyxiQDpwO/B+5qpdmDwJ+An7ozFhER6f5uf30NAGNTolizrxBwDG9raTHRv3y2hb15Zby75gAAKTEhpOeX8dNThuLn23VGiBtjSIwIalIJbn9BOVP+8CU/nzeEiyeleDA6ERHv4+4iCI8DPwPCW9ppjBkP9LXWfmiMUQIkIiJHpbCsmjdW7nM9v+aFFU32ZxVXkhAeSJ2FrZnFfLT+IH/9anuTNtcc359LJqcQ5O/bKTEfqVTnekOnj04iPb+ctfsK+Plb67lgQl98W6rLLSIiLXJbAmSMOQPIstauNMbMamG/D/AocHU7znUDcANASorudImISFMPf7qZV5a1Pi/mxW938/cFO1rcFxbox4jeEV06+QFHmeznfziJyf1jeHHJbtbuKwAcCV39OkUiInJ47uzjPx44yxizG3gdmGOMebnR/nBgJLDA2WYq8F5LhRCstc9YaydaayfGx8cfultERHq41ubDzBri+MxoLfm578zhbLj/FP7zo2ldOvmpN3tIAqGBfswanODatnpvgecCEhHxQm5LgKy191hrk621/YFLgK+stVc02l9orY2z1vZ3tlkKnGWtXdHyGUVERFoWEtB8QEPfmGD+euk46keHXTt9QLM2vT1c7vpoDe8dwfbfn0pogC9bM4s9HY6IiFfp9FmexpgHjDFndfZ1RUSk+8orrQIgNT7Ute1P540mPMjfVdAgKTKIe04dyvS0OH5y8mAAkqNDOj/YDuLn60NaYrgSIBGRI+TuIggAWGsXAAucj+9tpc2szohFRES8X3FFNZ9uzOT88X2oqq0jr7SKEb0jePvHxzHk158AEB0a4GhsHX8kRgRx5pje/OiEgVhrOWl4L1epaW81OCGMzzdlsvFAISN6e750t4iIN+g6dT5FRETa6Y8fb+Ynb6zlo/UZDPn1JyzankPvqGAC/Rrm8cTUJ0BOCeGBrsfGGK9PfgDOHtuHyuo6zv37t3y1ORNrradDEhHp8pQAiYiIV9iWWcy8xxcy7/GFvOqs+Hbzq6tc++PCHAnPUGdiExXiD4B1dgElRgR1ZridYvqgOJ67ehJVNXVc88IK3lyZ7umQRES6vE4ZAiciInIsXlu+l3veXt9kW7C/L+XVtYAj6blqWn8AXrp2Chv2F7p6g0b2iWT13gISIgLpjgbENcx72pKh+UAiIoejBEhERLq00soa7n9/Y7PtL107mV++s569eWV8csdM1/b48EBmD20oE/2vH0xi44HCFivFdQeNh/bV1GkInIjI4XTPTwMRES+zLr2AYUkR+PtqZHK9DfsLAfjNuxuorKnjRyekMmdIAvvyy3n2m52M7RvFe7dMp7Km5TWA6sWEBjBjUPddQ86nvs438MK3u7luxgCvrm4nIt6lqqaOX72znlvmpNEvNvTwB3QB+qQVEfGw7VklnPW3xTzy6RbXtqKKajYdLHI9t9ZSWlnjifA85oy/LuKMvy5i9d4CHj5/NPecOowpqbFcMCGZT+6YiZ+vD0H+vkQG+3s6VI97+doprsdvrdzvwUhEpKdZs6+AN1am8/O31nk6lHZTD5CIiAe8v/YAsaEBHJcWx968UgAWb88ht6QSgAueXsKunFLevHEaT3+9gy82ZeHva5g9JIFTR/Xim6053HBCKkN7RXjyZRwzay3GGNfjxdtzOT4tlrKqWlebm2YN5MKJfT0VoleYPiiOT+6YwbzHv6GypvbwB4iIdJCi8moAlu7Mo6yqxiuGG3f9CEVEupkd2SXc+tpqADY/OI+P1mcAsPFAERN+9wX+vobqWsdcjgueXuI6rrrW8tn3mXz2fSYA/eNCvToBWr03n3P//i1v//g4xqdE88G6g9z62mpevGYy4UGOj6dnrpzAySN6eThS7zC0VwQxoQEUOr+MiIi4Q2FZNZHOKpsHCsq57t8rXPv+9PFm7j97pKdCazcNgRMR6WRbG1XqeuLLbc1KF1fXWi6ckMyMQXEAnD4qiUU/n91kmBNAXmmV+4N1o8XbcwD4ZIMjAXx7lePvYX16AY99vhVjYEQfLe55JCKD/Smq6FlDJUXE/erXGPtyUyZjHviM5bvyKK+q5ZoXvmvS7roZqZ4I74ipB0hExM0qqmv5x9c7uXRyXxIigtieVeLa9/cFO1o85pELx7A7p5RnF+3kpycPJTLEn+ToEJ67eiLXvOC425ZdXMnO7BK+3ZHLBROSCfL3bfFcXVV9vEXl1WQXV7JwmyMh+vNnWwH42bwh9IkK9lh83igi2F89QCLSobKKK5j8+y/522XjWLE7H4CL/rGEXhFBZBZXuNp9cOt0+sZ4RwEWJUAiIm726OdbeWbhTnJKKrn75MEs2ZlLn6hgjIH0/HLOHtub00Yl8d7aA0xPiyM21LGgZ/+4UH53zqgm55ozNJGPb5/Bfe9tZG9eGec99S0FZdVUVNd6zZ23evXjxvPLqnhp6R5q6yyTB8SwfFceQ3uFc+30AR6O0PtEBPm5/l5FRI5URXUt7605wAUTkl0VJtfuc1Tk/Okb67hsSoqrbUZRQ/IzPS2OkV7UY68ESETEjay1rNrjuGP20tI9vLR0DwC/PG0ofaJC+PX/1nP3SUNIiQ3hlHbOdRmWFEFCRBDvrz3g2vbmynSvS4CySxxD+D7dmMkXm7I4fXQSj188FmshwE8jtI9GZLA/+/PLPR2GiHiR7OJKLJaE8CCe+HIbf1+wg8gQf9dn0o5sx6iF8upa/rVoV7PjX71+CscNjOvUmI+VEiARETepqK7l1P/7hl05pc32XTc9FR8fw2mjermqoB2J+sUvQwJ8OXdcHz5cf/CY4+1s9RXvAKKC/fn9OSO1DtIxitQQOBFph53ZJaTnlzNzcDyTfv8FAB/eNp2tmY45qvvyyrjj9dVcMKEvWzOLSQgPZHxKNJ9szHCd44zRSXyw7iBpCWEeeQ3HQgmQiEgH2JZZzLasEk4blQRAXZ3lxW93u5Kf0cmRnDA4noTwQJIig11DC44m+QGY0C+afy3aRVlVLbFhgRSWV1NbZ/H1ObrzdZbiimrmPvo1/3fJOA4UljM1NYYLJ/RlxuA4okICPB2e16ufA1RVU6deNBFp1Zy/fA3Ak5eNd207/YlFDEkMB+B3H24C4H9rHCMNZg+J594zh7sSoPvPGsGlk1O4YWYqCeFBnRl6h1ACJCLSAU56bCEAu/54GsYY/r5gu2syP8CZo3tz/cyOG6I2d1giAD86IZWYEH+shcLyamJCu3YSsWF/EZlFlVzyzFIArp8xgPMnJHs4qu5jamosTy3YwX++28uV0/p7OhwR6QLyS6vw8TFsyyzm1eV7efj80a59N7+6qknbLZkNVUpT40LZ6byJNzU1lt5RwXx02wyiQ/1JinQUqBmdHOX+F+AGSoBERDpQen45Ly3dwxsr9gEQHuhHcWUNQQEdW6EtwM+Hbb8/FT8fw3vOuUB5pVVdPgGqrbOux8bAXScN8WA03c/MQXH0iQpm5Z58JUAiPcT2rGJ8fXwYEBfabF9dnWXe/y0kt6SKAD8fyqpqmTM0oUmbl66dzL8W7WLBlmwGJYRx8+w0Xv9uL/efNZLMogp+/Moq5o10zAca3tt7155rTAmQiEgH+uPHm1wLm9510mCC/X35/UebiAr27/Br1c+XiXYOHcsv6/rrAhWUN8TYLyaE4A5ODHs6YwxJkUFkFlUevrGIdAtzH3WMQNj90Omubc9+s5MVu/O5bEqK6/2gpqoWgD84h7cBDO0VzoxB8cwYFE9WcQUxIQH4+fpwzrg+AAzpFc76+04+6uHaXZUSIBGRY7Rhf6Hr8bKdea7HZ4xOol9sKMnRwa67Z+5Q3+uT38UXRl2fXsjafQWu5/1buFspxy4xMohNB4o8HYaIeEhtnXXN4SmvdiQ94UF+FFfUMHNwPAu3ZpOWEMYdcwcxLTXWdVxrc3m6W/IDSoBERI7IvrwywoP8CPDzoaSihoSIIM746yLX/tzSKnx9DPefNYLUeEdlnFOdhRHcJSrE0bvU1XuAzvzboibP+8cqAXKHxPAgFhRleToMEfGATzZkcOPLK13Pv96aTVJkEP/90TSyiiuoqK5j6c5cHjpvFBP7x3gwUs9SAiQicgRmPDyf3pFB9IsNZcnOXLb9/tRmbT65fQaDnJV0OkN9D1Beadctf7x4e06zbd1lLHlX0ysykNKqWoorqgkP6vihlyLSNVXX1vHTN9e6nocG+FJaVUtaQhh9Y0LoGxMCwIb7TunxVSKVAImItFP9BP4DhRUcKHSsgP3595kAPHbxGAJ8fekXG9KpyQ9AsL8vgX4+FHTRHqDaOsvlzy5rtn10svesGu5NejmrM+0vKGdoLyVAIt2ZtQ2FZcbe/xmlznk+AAMTwliXXsj4lOgmx/T05AdAfwMiIu2UU9J8YvmPX1mFv69hYr8YTh+dxMg+nf+l3hhDdEgAeV10DlDjoXmNk560eO9bPM8b1K/jsSWj+DAtRcTbVdbUuR43Tn4Arjl+AJHB/vzw+P6dHFXXpx4gEZF2ynD2+jR2wuB4Lp3c1zW0wFOiQwN4f90BbjtxkMdjOVRuSUMC9NszRzAsKZz8smr8fHUPzh1S4x1zq25/fQ3DkiIY3Mk9kiLSeUoqa5ptu/3EQQxKDOOM0b1d1dykKSVAIiLtsHBrNlc9t7zJto33n0JoYNd4G40I8qOiuo4Ln17C0l+e6Olwmsgtbeg5SwgPJCTAj5CArvH31h35N0os52/OUgIk0o2VOhOgqakxXD8jlTdWpHP7iYPw8el+lds6km6/iYi0w1ur0ps8v21OWpdJfgAKyhwFEDKKmvdSeVr90Lw75w7ucr1T3dX8n8wCYFtWiWcDERG3Kq10DHu7+rj+nDgskaevnKDkpx2UAImIHMZNL6/k3TUHiAsLpP5zJSGi5fUSPCWzuCHxWb4rr42Wna8+AbpsSoqHI+k5BsSFMnNwPN9rPSCRbmf13nw2HijEWktplaMHqCvdkPMG+tsSETnEgYJyAHpHOappfbwhA4ApA2L4w7mj+Nv8bVwwIdlj8bVk3ohevP7dPgAu+scStv/+1C4zx6Z+DlB0iCqSdaYhiWEs35VLXZ3VHWGRbuTcv3/rejxjUByAhhUfoa7x6Sgi0kWUVdVw3ENfcdxDXzXbFxniT2SIP786fThB/r4eiK51D54zknkjermej7zvU4/FsmpvPlc9t5zpf/qKDfsLeXt1OgPjQ7tMQtZTpMaHUVFdx4HCck+HIiJu8s02xxprYeoBOiL6NBIRaWT4vQ2Jw8yH51NY1rC4aGV1XUuHdAn+vj6EBDQkZRXVddTV2TaOcJ/z/v4tC7dmk55fziXPLGVfXjl3zB3skVh6stQ4RzW4HdmlHo5ERDpKS+/rxkCf6GAPROO92kyAjDHTjDFPGmPWGWOyjTF7jTEfGWNuNsZoBTsR6VYqa5quobA3r4yXl+1xPY/q4kO46seC18vtAusC1ZdoHZakSmSdbWCCY50lzQMS6T4Ky6ubbUuNC1UP0BFqNQEyxnwMXAd8CswDkoDhwK+BIOBdY8xZnRGkiEhnOFDQvILao59vBeDccX24++Su3Yvxi1OHMbRXQ6Kxv6DrDH1S9bfOFxcWyJjkSP63er+nQxGRDtJ4WYF6ozywALe3a6sH6Epr7bXW2vestQestTXW2hJr7Spr7V+stbOAb9s4XkTEq6TnlzXbVltnSY4O5o/njeryk0wHxIXyyR0z+fj2GQDsz+/8BKjxkMEzx/R2PQ7061pzpnqKU0clsSWzuMm/i4h0bda2PHz5/bUHmPvowibbbjtxELdriPERa/XT3Fqb0/i5MSaicXtrbd6hbUREvNm+vKYJw3nj+xAdEsCFE5O7XNGDtiQ7x4Lvyun8NWC2ZhUD8LtzRnLhxGTOG9en2dA86Twpzp639IIyIkN0l1ikq5v6hy85dVQvfnvmCAA+2XCQYUkR9IsN5dbXVjdrf9dJSn6OxmFvZxpjfgTcD1QA9SmpBVLdGJeIiNt9siGDJ77cxms3TGXzwSLuf38jwf6+/P2K8QT5+TJ5QAy+Xlg+ODzIn4HxoazeW9Dp1964vxCAk4YnEujny+yhCZ0egzSoL+V+oKCCEb2VAIl0ZRXVtWQUVfD84t1cPKkvGYUV3PjyKmYNiWfygJgmbf951UTXkg1y5NoznuMnwEj19ohId/P4F1vZnFHMmPs/c20bGB/K7CHe/6V9XEo0b65MZ/muvGYfnO608UARcWEBJIQHdto1pXV9nAnQ/haGd4pI15JT0jC/56GPN7tuYJRX1fLwJ1tc+0b1ieSk4YmdHl930p4y2DsAvXOKSLdRVlXDL99ZT3Zx88mkZ4/t44GIOt6lk1MID/Lj7jfWUF3bUL77nwt38tyiXW655pPzt/PGynTG9o3CGO/rOeuO4sICCPTzYZ8H5oOJSPtYa1mXXtDkMymjsIJdzhL2GUUVBDuHYb954zReu2GqR+LsTtrTA3QP8K0xZhng+pex1t7mtqhERNwkq7iC8/7+LenOL4Q/mpnKvJG96B0VjL+vD5HBXbvUdXtN6BfNE5eM44cvfMd/vtvHqD6RLNiSzWNfOKraXTN9QIddy1rLL9/ZwGvL9wIwqX/n9ThJ24wxjE+J5stNmfz69GFKTEU8qLKmlqLyGuIb9ZBba/l0YwY3vrzKtW1MciR788oocBYv2ZPr6If4w7mjmKj31w7RngToH8BXwHrgiFcBNMb4AiuA/dbaMw7ZdxeOUts1QDZwjbV2T/OziIgcmcLyaiKD/ckorODhTzaTXVLJiUMTWLgtx5X8gKNS2chuWkJ01pB4JveP4b73NlLjhkVRt2YWc+2L3/Gn80e7kh+AU0cmdfi15OidPbY3v3h7PduzShiUqPWYRDzl5ldW88WmTHb98TSMMVhrmff4N81KWw/vHcna9EKgmt6RQRworHBuj/BA1N1Te4bA+Vtr77LWPm+tfbH+5wiucTuwqZV9q4GJ1trRwJvAw0dwXhERsoorePTzrdQ4h3llFFZw5b+WMeb+z/hmWzZPzt/O26v38822HO57/3u+2pzlOvaGmandNvkBx93/p64Yz1ljehPg1/Tt/tBFX49UdW0d9767gX155Vz2z2Wu7c9cOYGUWK3505XUJz3pGgYn0iFue201z36z84iP+2JTJgBFFY7KmGvTC9mSWUxOSdNFqxsnOtfNaKg5NkIJUIdpTwL0sTHmBmNMkjEmpv6nPSc3xiQDpwPPtrTfWjvfWls/v2gpkNyuqEWkRyutrKG0sobq2joe/mQLT3y5jQVbsimuqObKfy3jm22Omi3XvriCj9YfBGBQQhj/vmYy549PZmivcG48YSC/PG2YJ19Gp4gNC+TRi8fywa3Tm2zPKmo+/6kttY16kPbmljH7zwtYujOvWTsteNr11BdCSFfFKJGjUlxRzWOfb6WqxnGj7b21B/jdh63d23d4btEu/vhRy22yiytYtTefc55cDIC/b8PQ1D+eN4pk5/9ZgCum9gMgMtgff9/2fG2X9mjPELhLnX/e02hbe8tgPw78DGhPn/u1wMct7TDG3ADcAJCSktKOU4lIdzblD19SUllDbGgAUwfGArB8dx7LduWyPbuEm2cP5J1V+zlQWEFuTRWXT0nhl6cNIzTQj5mD4z0cvWf0jW6amOzNK3OtF3S4eSG///B7nl20i5eumcLy3Xl8f6CoWW/Cg+eM5PnFuxgQF9qxgcsxSwgPxN/XqGSuiNOXmzJJiQlp95DQz7/P5P++3MaMQXGMSm7fqIFPNmSwel8+t8xJIzyo6dzSxz7fxtKdua7n790ynfLqWqKC/UmND6OwrJpTR/bi1FFJBPj5sOSeOfj5KPnpSIdNgKy1RzVT1hhzBpBlrV1pjJl1mLZXABOBE1qJ4RngGYCJEyd2/EB2EfEalTW1lFQ6hg/kllaxJ9dRJeeZhY7hCKnxofz0lKFcOz2VPbmlhAf5MSAuzCvX8+lIwQFNF3L9+VvrKCir5swxvRnXN4pVe/N56PzR1NVZfA75u/rnN46qcVf8q2Go27nj+vDO6v2u55dPTuFK551K6Vp8fAy9IoPYryFwIlhruf31NcwYFMeTl42nzlr8DtOzstNZja2oovqwvefVtXXsyysju6SS6lrLN9tyGNornCe+3OZq86FzZMKpI3tx75nDSYoMbnKOyBB/nrpiguv5ofvl2B02nTTG+BpjzjLG3GaMuav+px3nPh44yxizG3gdmGOMebmF888FfgWcZa09sjEZItLjbMssafJ8w/6iJs+jnFXcYkIDGJcSTVpCeI9PfuqdN85R4vsvF44hPb+cksoaXlu+l5+9tY7Xv9vHZxszSP3lR2zPcvwdV1TX8tXmTNcQqt6RQa5zTUuN5cu7G+5ZHZo0SdeSFh/G6n35LNuZ67qBINITHCgo562V6a7nmUWVlFTWsPFAETe+vJK0X7U4+IiCsirO/ftiNmcUsSPb8Z5YVF5DVgvLJzT25sp05vzla3blOJKmH7+yijl/+Zr/rTnQrO2gxHAlNx7SniFw7wMVHGEVOGvtPTiHzTl7gH5irb2icRtjzDgcVebmWWuzDj2HiMihVu/NB+C88X14e5WjB2JaaixLnMMJ3FDsrNv40wWjOW98MtMHxXH3G2ub7f/zZ46F9uY++nWzfTeeMJCfzxvCgHs+AmBKagz9YkP54q4TmqwzJF3TWWN7c+d/1nLxM0uZOyyRZ38w0dMhiXSKHzy3nG1ZJZw8IpHwIH9XMrM3r4y9eY5p6NZajDHkllSSW1rF4MRwlu7MY/XeApbuyHX1ABVXVJNZVNHk/N/tzuOxz7cyPiWaW+aksX5/oWufj2n5M8nf11Bda+kTFdR8p3SK9iRAyc4qbR3CGPMAsMJa+x7wCBAGvOEcg77XWntWR11LRLqHrKIK3l1zgAn9o/nNuxvpExXMXy4cw/CkCNbvL+SOuYM55fGFDEkM54/njfJ0uF2Wv68P0wfFAY71JP782RZOGBzvGspWf8eyJSkxIU3mCqU4ix2kJYS5MWLpKBP7NdQu+v5AYRstRbqPpTtz2ebs0X5+8W5eXrqnxR6covIaIkP8Offv37I3r4x/XjWR1fscN9u2ZpW43htzSqp44qvtruM+WHeAW15dDcC3O3JJPGSo6U9PGUp8eCA/aXTD6eKJffnDeaP4aP1BThulJQM8pT0J0MfGmJOttZ8d7UWstQuABc7H9zbaPvdozyki3Vt+aRXRoQH897t9/OytdQDUf/8+bVQvjDFNyoNu/d2pngjTa102JYXLpqRQU1vHrCHx3P76GqprHbcqb5uTxonDEkmMCGLe/y2koKyamYMdidPzV0+isqZOC2p6mcYLL/r7aTK19AyXPLPU9fjRz7e22m5vXhnLV+W5eoSu//cK175XlzWscfbmynSyiytdPTj1yU+93JLKJjeSokP8OX98H1cCdNuJg7hscgq+PoYzx/Q+thcnx6Q9CdBS4B1jjA9QDRjAWmtVjFxEOtyO7BKu+tdy9heU88IPJ7mSn74xwezLK+fccX34ySlDPBxl9+Hn68PZY/vw9/k72JJZzG/OGM610xtq37zz4+Opqa0j2VlFbvbQBE+FKscgyL+hCIaf5mtJDxca4EtwgK9r/Z1fvrO+ydC11uwvKKdPVDCXTUnhkU8dQ4a/uGsmcx9dCDjmpKbnl3H1cf2ps5Yzx/TGGMPcYYmkJYRx10mD3fei5Ii0JwF6FJgGrLfWanS9iLjVfe9tZL+zXO97axsmjX502wy+3prNrCEJBPr5tna4HKXBvcLZklnMyEMW2lNZ6+5Ha4lIT9DW3ERjDGkJYeSUONYya0/yU++ccb0JdVbVvGhiMmkJDaW06xc6nT00gRMaLbmgOXddT3veBfcBG5T8iEhHueHfK/jnwqaraC/alsODH3zvWsQUHAtuArx4zWTCg/w5Y3RvwgLbc99GjtSY5EiC/H2arEAu3VOdPs6lB9jnHM7WktNHJTF5QGyL+1744ST++6Np+PoYXr1uCs9fPYkPb2tYSPqcsX0oq64FICokAIAl98zh9+eOdLUZnxLVAa9A3Kk93yR2AguMMR8Drplj1tpH3RaViHQ7FdW1WAvbs0r47PtMPvs+k+tnOubwWGubrDFTb8UexyRUVcpxv6um9efUUUnNFuyT7qN/bAi7c8soLK/2dCgiblVTW8cH6w62uG/ygBgePGckVbV1pOeXuaqJ1jthcDzGGHb84bQm24ckOnrJByWG866zpHW0MwFKigzmool9Sc8vp1dEkN5HvUB7EqBdzp8A54+IyBGb9cgCMosr+MnJjvk7/WJDXPtW7yto81itk+B+AX4+rvV+pHt6/9bp/OZ/G/hoQ4ar7K9Id/T01ztaLXrQNzqEAD8fAvx8ePSisezILmXtvgJ+c8ZwRvSOaPX/xds/Po5aZ+/pNdMHcKCgnCumprj2+/v68PN5Qzv+xYhbHDYBstbe3xmBiEj3luFcO2GHsyRpaWUtLy3ZzWmjkljqXMOnsZtnD+TJ+TuICwsgVMPeRI5ZeJA/Q5Mi+N+aA/z1q+1cNiWFuLDAwx8o4kUqqmtZvL3hM+WmWQPZmV3CpxszW2z/t0vH8eWmTH5wXP82bwo0/hyKCQ3g0YvHdljM0vla/VZhjPkn8IS1dn0L+0KBi4FKa+0rboxPRLxUjXMCqt8hE67fdq45k1NSyW/e3cjz3+5mZ3YpqfGhvH79VCb/4UsAbp0ziOKKGi6dnIKIdIzYUMdAjkc/38r8LVm88+PjPRyRSMepqqlj3AOfU+6cowPws1OGYIzhrZXp3P3GWpKjm/Z0940J4erjBxx6Kunm2rqt+iTwG2PMKGADkA0EAYOACOA5QMmPiLjsLyhnb24Z0wbGcv7TSzDA/24+nopGH0aHql9h+4IJya61Si6dnEKQvy8PnD2y1eNE5Mg1Xg9oX155Gy1FvM/i7TlNkh/A1atz3vg+hAT4Mnd4oidCky6m1QTIWrsGuMgYEwZMBJKAcmCTtXZL54QnIt7k5Ee/prSqli/vPoG1znk93x8o4tBRBVdMTeHlpXs5f3wyKTEhTB8UywTnSvWbH5xHgMr0irhF4yFv+WVV1NTWNeulFfFGi7fn8MMXvmt1vzGGU0cldWJE0pW1Zw5QCbDA/aGIiLfZmV3CgLhQ1x220irHnbevt2S72pz2xDdNjkmNCyU+zFHVLT48kNvnDmqyv/GCjSLSsRr3ANXWWTKKKlyL3Ip4s5XOqqGN/e4cjSKQlum2j4gcVllVDf9atIvaOkcFnD98tImHPt7MnL98zbPf7AKaLjq3fJdjcbngFpIZP1/DueP60CcqmMunaH6PSGeKCW1azPX3H27yUCQiHSejsKLFqm9XTO3ngWjEG6i0kogc1uNfbOOZhTvpFRHE6aOTeKbRIqbfbM9h3sheXPyPJa5ty3blEhcWwIKfzmbRtmxufHmVa5+fjw8psSEs/sWcTn0NIuIo1QuOmxPl1bV8ujGDiupa9byKV3v66x1Nnl86uS/Hp8V5KBrxBm0mQMYYX+BP1tqfdFI8ItIF5ZZUAVBcUU1xRdNFFKtr6vjjx5s4UFjh2pZfVs3wpAjCAv2Ymtp0tW1/P3U8i3jSZ3fOJDY0gO9253Pjyyv5/mAR41OiPR2WyBFbu6+AfyzcwUfrMzh5eCLhQf4s2ZHDH88b7enQpItrMwGy1tYaY6Z3VjAi0jX5OIsYvPbdPt5be6DJvtKqGrZnl7ieT0uNZcnOXAbEhQIQFdJ0yI2/jxZfFPGkwYnhAIxLiQJg9d4CJUDiFerqLN9sz2FvXhkXjE/mxpdXctB58+2h80c3G+Ip0pr2DIFbbYx5D3gDKK3faK19221RiYjH/fbdDZRX1/LwBWNwTv1xVXZrbF16YZPnr14/hZV78ukd1bDWwmMXjyG3pIrffbiJ1PhQd4YtIu2UGBFEn6hgVu/NB7QOinRtNbV1XPmv5SxxLpz9+feZruQnNS5UyY8ckfYkQEFALtB4wL4FlACJdFN1dZYXl+wBHAuS5pVWttn+jNFJfLDuIOAoNTqxf0yT/eeOSwYgKTKYOUMT3BCxiByNsSlRrN5b4OkwRA7rm205LNmZy5ljevP+2gMs3JrNoIQwtmWVUFlTd/gTiDTSnjLYP+yMQETE8z5af5D1+wvZllns2vaf7/a57rIdau6wRH56yhD6x4Vw8+y0JpXgWnL6aK3BINKVTB0Qw4frDvLd7jwmHXLjQqQrWbozlwBfHx4+fzS3zE5j/pYsLpnUl5MeW8hvzhju6fDEyxhrbdsNjBkMPAUkWmtHGmNGA2dZa3/XGQEeauLEiXbFihWeuLRIt3bFs8tYtD2nzTZTU2NYutNR4vqsMb2576wRGnYg4sXKqmqY/ecFWAuf3jGTaP1/li7qnCcX4+9reOPG4zwdingJY8xKa+3Elva1pxzTP4F7gGoAa+064JKOC09EPC2rqKJZ8vPEpeOaPH/0ojHcNmdQk/1KfkS8W0iAH49dPJas4krmb8nydDgizezLK+PEvyxgzb4CRvWJ8nQ40k20JwEKsdYuP2RbjTuCERH3yiyqYGuj4W31PcCLdzTv+RnirBRV74TB8c0quomI95s6IJbY0ADmb8n2dCgiTRSUVTHj4fnsyHbU4BqaFH6YI0Tapz0JUI4xZiCOwgcYYy4ADro1KhFxi1mPLODkxxYCsDO7hAH3fMSry/ay0/nhctucNMID/egXG0JqfCjv3zKdP184hm9/MYfYsECiQ/09Gb6IuIGPj+GM0Ul8vP4ge3PLPB2OiMv+gvImzwclhHkoEulu2lMF7mbgGWCoMWY/sAu43K1RiYhblFfXAo6en7XpBQD88p31ACRFBnHXyUO4Y+5gfJxr9YxKjmRUcqTr+GhnD9CsIfGdGLWIuNvNs9P4z4p9PPr5Fh6/ZNzhDxDpBFnFjgqkz1w5gaziSsb2jfJsQNJttKcHyFpr5wLxwFBr7fR2HiciHlZVU8dVzy1nxe68JttLq2rJKmpa2rqPc90enzYWKg3y9+XD26bz5GXjOz5YEfGYhIggLpmUwkfrMyit1Ch36RqynZ9Tw5IiuGJqP4zRQtrSMdqTyLwFYK0ttdbWTx54030hiUhH2ZpZzMKt2Szb1TQByiupYk9e06EuA+PbN7RgRO9IQgPb03ksIt7k5OGJVNXWsWRHrqdDEQEc81YB4sMDPRyJdDetfosxxgwFRgCRxpjzGu2KwLE4qoh0cfUFD0oOuaObW1rJ3twyxiRH8s6Pj2fpzlzSNLZapEeb0D8agA0HCpk7PNHD0UhP9+2OHP7y+VbAMfpApCO1dRt3CHAGEAWc2Wh7MXC9G2MSEafq2jr8fY9+xOkWZwJUXFHdZHtuSRUbDhRyyvBe+PgYjkuLO6Y4RcT7Bfr5EhnsT35pladDkTZUVNd2+4Qgt6SSO/+zBoBJzsRcpCO1mgBZa98F3jXGTLPWLunEmER6tJLKGjbsL+SZhTv5anMWb9w4rdUV2mtq66izEODnSJIKyqqalKrekuHsAapo2gO0am8+BWXVjE2Jcs+LEBGvFBMaQF5Z9eEbikes3VfA2U8u5sVrJnPC4O5ZjKaiupYz/rqIzKJKXr1uCpMGtPz5J3Is2nNrOdcY86UxZgOAMWa0MebXbo5LpMe65oXvuOSZpXy12bEo4YI2Fie85JmljLzvUwBW781n7AOf89nGDNf+rRn1PUA11NZZ1/b6c49uVOFNRCQmNIC80srDNxSP+M5Z0Ga+8z08q6iC0/7vG9fNru7g4w0HOVhYwYNnj+C4tLhjGgUh0pr2/Fb9E7gHqAaw1q4DLnFnUCI92fJDChY0yltcKqpreWXZHlbsyaeqpg5rLRv2FwLwwTrHMl2F5dUcKKxwPV6zL991/Gbnh2VqnOb9iEiD6JAA8krVA9SVZBVVsDO7BABfZ5XO2jrL84t38fO31vH9wSIe+XQz32zz/oVs6+osTy3YQVpCGJdP6efpcKQba08ppxBr7fJDSg+qRqZIJ6mvglPPWsvQ33zSZFtBWTWr9hYAkFNSycYDhZz+xCIAwgP9WLEnn/OfajqSNSE8kOCA7j2OXESOTEyoP+v3F/DBugNEBvszY1D3HGblTU54ZAHl1bWs/e3J/Oe7fQC8tHRPkzZfbMrii01ZbH5wntfODzr/qW8praxha2YJD503qs0lGUSOVXsSoBxjzEDAAhhjLgAOujUqEXFJz3OshJ1fWoWPj+FPn2xu1mbcg5+7Hq9LL+SdVfsBmJ4WR2SwPx+ub/5fNiUmxE0Ri4i3CvTzJbOoklteXU2Anw9bf3eqp0Pq8eoXsL74H0tcvfetOVBQTmo7lzToSmrrLCv3NIxSmJoa68FopCdozxC4m4F/AEONMfuBO4Cb3BmUSE/me8hdr+3ZJVhrGffg54y5/zNeXbbXta+lNeFKKmt4dtEuRvaJ4KVrJ5Oe33S9n1F9HPN+UmKVAIlIUyN6RwCOhZGraur43+r9lFfVejiqrquwvJrHPt9KVU2dW85vbcMY6EOTH18fQ0xoABP6RXPz7IEAfLox0y1xuNveRuvSBfj50E+fT+Jmh02ArLU7rbVzgXhgqLV2urV2t9sjE+mByqtqXcUKXrluCr89czh5pVUMuOejFts/cPZIAI4bGMv8n8xiwU9mkRwdDMDUAbEYY6iqbTqJ6N4zh3PaqF7cPDvNja9ERLzRxZP6su6+k3nkwtEA3PGfNTz99Q4PR9V13fWfNfzfl9tYujOXFbvzKOzgCnrZJa0XpNj0wDyW3nMi/7lhKpdMSgHgT59sZtPBog6NoTNsbhTz6aOSMC3d3RPpQIcdAmeMiQKuAvoDfvW/lNba29wZmEhP8dDHm3n66x3sfuh0/r1kNwD/vGoix6fFtdjD09jghDD+ceUEjhsYS3iQPwAv/HAyC7ZkcdkUxwfi3y4bx8Kt2YzqE4mvj2FcSnSrZbVFpGczxhAR5M+UAbH8+vRh/PmzLSzensOdJw32dGhdjrWWL53V2HJLK7nqubUMjA/ly7tnddg19jmHQB9qdHKka/kDgF6RDevTb84oYlhSRIfF4G4b9hfyxaYsgvx9uO/MEZwzro+nQ5IeoD1zgD4ClgLrAff08Yr0YPV3Vx/5dDNPzt/BmORI5g5LAGB0chRxYYHkHHIX8KKJydTUWcamRBHo13TCa1pCGGkJDWPAB8aHMdALx4SLiOf4+hium5FKTkkVT3+9g//7Yhu3zx3k6bC6lO1ZJa7H32zLAWBHdimFZdVEhvh3yDXynIvS/mBaPzZlFFNWVUPf6BCeumJCk3b+vj7MGZrAV5uz2HywGMZ1yOXdbvH2HC5/dhkAZ4xO4pLJKR6OSHqK9iRAQdbau9weiUgP8e32HBbvyGF/fjl9GxUieHK+IxFKjQ9zdf+HBfqx4tdz+WRDBkt35nL1cf0J8vdtcrdPRMRdrp0+gHdWp/Pq8j3cOidNlbka+Wh9w5prC7c2lKDeklnM5A5avDPfmQBdNyO1yedFS567ehJn/nURS3flYa31imFkzy/eDTjm/dw6Rwm2dJ72JEAvGWOuBz4AXLehrbV5rR8iIq25zHm3qzV9ooKbbZs3shfzRvZyV0giIi2KDw/knlOHccd/1rBkZy7Hp8V5OqQu4cn523nsi60cnxbLtztyySmpcu1buSeftIQwYkIDjvk6+WWO80a381wXTEjmt+9tZF16IWP6Rh3z9d1tbXoB543rw+/OHUlIQHu+kop0jPZUgasCHgGWACudPyvaewFjjK8xZrUx5oMW9gUaY/5jjNlujFlmjOnf3vOKdFehgfoQEJGuY97IXsSFBfDCt7s9HUqXsWK34x7wr08fTrjzPfv4tFhCA3z50yebGf/g55RX1WKt5dvtOU2quR2J/LJqAnx9CG3nmm0zBzvWbdqRXXKYlp6XWVRBdnElo5IjlfxIp2tPAnQ3kGat7W+tHeD8ST2Ca9wObGpl37VAvrU2DXgM+NMRnFfEK4UdJsFJCA/spEhERA4vyN+XE4cmsnJP/lF/ke9u9uSVcerIXgxLinDN9+kfG0paYrirzQfrDvDmynQue3YZA+75iA37C9s8Z0V1rWsoXV5pFWv2FZBdXElUiH+7h7PFhTl6ig6dN9oVffa9o2T3hH7RHo5EeqL2JEDbgbLDtmqBMSYZOB14tpUmZwMvOh+/CZxovGHQqkgbvtqcyR8+ajnnr6mtcy1q15KnLh/PuaqAIyJdzPDeEeSVVpFRVNGjk6CyqhreWLGPndml9IsNBaCy2lEfakBcKKlxoa62Gw8UNVnf5g8fbaKmto4N+wux1nLvuxv4Zls2+aVVVFTXcttrq7nqueW8uTKd2X9ewDlPLuatVenUHcHfd1igH0H+PmQXd/0E6L/f7WN0cqRrbTqRztSePsdSYI0xZj5N5wC1pwz248DPgPBW9vcB9jnPV2OMKQRigZzGjYwxNwA3AKSkqEKIdG3XvOAYIfqLeUObTRjenVvmWuenJaeOSnJrbCIiR2O4c4HUaX/8imB/Xz6/aybJ0T1nsUprLb//cBP/XbGPoooaAAbEOV6/v6/jXvL0QXHszCl1HbPpYBGDEhsqcH67I5e0X30MwFXT+vHvJXv495I9za71kzfWAvCzeUN4+JMtTeYXHY4xxlk5tP3HeEJdnWVbVjGXT+nnFcUapPtpTwL0P+fPETHGnAFkWWtXGmNmHenxjVlrnwGeAZg4cWLPvfUkXqWgvJqY0ADKqmr4eH0GE/pF8+nGjBbbhgf68cI1kzo5QhGR9hnXN4qkyCAOFlZQXl3LHa+v4bUbprq+/Hd3RRU1PLtoFwCzh8QzaUCM64bVM1dNoLSylqG9IpieFsery/aSGh/K9weL8PNt+ct9S4lPYz+fN5SbZg0ks7DisNXfDtXS0gldTWZxBRXVdfRv1GMm0pkOmwBZa188XJtWHA+cZYw5DQgCIowxL1trr2jUZj/QF0g3xvgBkUDuUV5PxCMqqmsJ8velvKqWN1buc23PKq7gQEE5Z/x1kWtbeKAf01JjOXFYAr/7sGGY3BOXjmNCPy1OKiJdk5+vD/93yTgu+scSAFbsyee6F1fw+MVj212hzFtZa3nRWQDiqmn9eODskU32j+jdMITrtFFJfPuLOSzalsPP3lrH4u25XDYlhVeX7XW1efiC0Tzx5TbS81te5BTgplkDAbj/kGu1R1xYIOn5RzVzodPscvaUDYhVAiSecdhbN8aYQcaYN40x3xtjdtb/HO44a+091tpka21/4BLgq0OSH4D3gB84H1/gbKMeHuly0vPLePHb3c3Gvu/MLmHobz7h6a93MOzeT7j33Y2ufQ99vJl/Oe8Y1iuurOHBc0Zy7fQBbP3dqa7t0weptKyIdG0j+ziGwQ3t5RjV/vXWbB784HtPhtQpPv8+k0c/3wrQruUIekcFM6F/w8T+u08azBTnukBDe4Vz6shefP3T2fzk5MHNjvX1Mbx4zeRjindYUjhbM4vZl9c1k6BdOaX83xfbABgQrwRIPKM9Q+CeB36Lo0rbbOCHtK94QouMMQ8AK6y17wH/wrHO0HYgD0eiJOJxC7Zk8bsPNzGpfzS/OWM4f/hoEx+tz2B47wj8fAz7C8o5Y3Rv1jur+jz08eYWzpFNZHDz1cDTEhxjwgP8DJ/fOZOSypoeM4xERLxXSIAf3/5iDlU1dcz68wIASqtqPBtUJ1i6s2HZw4Tw9i1CnRoXytxhCZw1tg+xYYE8+4OJZBVXMjC+YU7QRZP68ufPtrqe940J5su7ZhHgd2yfB5dP6cffF+zgteV7+dm8ocd0Lne4+vnl7Ml1JGe9tai3eEh7EqBga+2Xxhhjrd0D3GeMWQnc296LWGsXAAucj+9ttL0CuPCIIhbpBDe8tJKqmjq2Z5UwJDHclci8vHQP7645AMAbK9I5bmBsm+cpLK92PT53XB/GJDetdjMosbX6ICIiXU/vqGAqGlWy9PPp3jdv6uosn33fMHczvp3LFBhjePYHDfM6w4P8CQ9qekMsITyIL+8+geziSi55ZimzBiccc/ID0CsyiOMGxvL2qv38aOZAV5nuruD+9ze6kh9ABRDEY9qTAFUaY3yAbcaYW3DM2wk7zDEi3q3RSLf73m8Y4vHV5izX46+3Zjf5InD22N6u5Cgs0I+SSsed0YHxoYQH+fPYxWPdG7OISCcI8m9YlLO1Sf7dxep9BU3m6kQEdeyCnQPjwxgYH8ZbN01jZAeWg752+gCueeE7Zj4yn79eOs61QKqnPb94t+vxv34w0XOBSI/XnlsNtwMhwG3ABOBKGubtiHQ7lTW1VNfVtbiv2Fn+tH6x0mW7GoZGTOrfUMTg7R8f53r8wNkj+d/Nx7sjVBERj6pfA6e72njAMcz58YvHcs+pQ93WYzGhXwyBfr6Hb9hOs4YkMG9kLwrLq7nqueUddt5jYa0lwDnc+8u7T+DEYYkejkh6svZUgfvO+bAEx/wfkW7nYGE5S3fmcu64ZOdq53DWmN6MT4lq0gNU77M7ZzL2gc8BOG5gLAPjw7h0cgpbMorZnVvKoISGTtIxfaM662WIiHSq/LKuvd7M0bLWYoxhZ3YpoQG+nD22t9cN17rphDQ+Wp9BeGDH9lodjZ3ZJWQWVVJVW8dvzhjeZC6UiCe0+r/CGPO4tfYOY8z7NBkQ5GCtPcutkYl0op+9uY5vtuWwK6eMJ750VKc5e2xvThyWyNaskiYlTGNCA4gKCWBCv2hySyr5++XjiQpxlIF98JymJUtDAnwJ6wIfPiIiHWliv2hW7MmnoKz68I29zO6cUmb9eQF3nTSYF77dTf/YEK9LfgBGJUdy8cS+LNiadfjGbrQlo5hTHl/oet7vCNc1EnGHtr6ZveT888+dEYhIZ/t6azb9YkLoHxfqmstTn/xAw2TXM0f3bpIAJUY4qta8eeO0Nj8UF/18tpIfEemW3rzpOH7+5jrmb/Hsl2t3WLknH8BV+jotwXt7K6JC/D2apGYVVzRJfkDLPkjX0Oq3M2vtSuefX3deOCJHZ09uKb2jggEwwO2vr+GGmaltDj/7wXPL8TGw84+ntzj2Oi7MkQAdWvUnJcZ5ncPcEUyO1l0uEem+YsMCyCutoq7O4uPjfT0krcksrmjy/A/njvJQJMcuItifypo614Ldne3QxV5/esoQj8Qhcqi2hsCtp4Whb/WstaPdEpHIEcoorOCERxbQLzakSXnNZbvyWPHruS0eU17l6PGpc/6GHyxsviJ3bJhjWFt8WNMEqPGq3yIiPVVCeCA1dZa8sirXDSNv987qdB7+ZIvr+ZQBMSREeO9aNVHOEtiF5dUeSTyyiyubPE9oZxlxEXdra3zOGc4/b3b+WT8k7graSIxEOlt6viPpaZz8AOSUVLIju4SB8WF8sO4Axw2MIybUkdQ0flO+8l/LOFjY9I4f4OoVigzx57GLxzAwPoznF+/mh8f3d9MrERHxHvXDgbOKKrtNAnT/IUVv6m+EeauoYEf8BWXVrn+vzlT/WRsT6ugt9OZkUrqXtobA7QEwxpxkrR3XaNfPjTGrgF+4OziR9sgpqWx13+X/XMZJwxN5aeke5g5LZGL/aDYeKGLF7oby1d9sywFgXEoUq/cW8LtzRnLo6LZzxyUDaC0fERGnhAhH0vOnTzbzwg8neWWhgHpfbc5k4dYcCsqq+fXpwxiXEsX5Ty0hOsS7E6D6RbwziyoY0qvzF97OLq7EmIYEKCq46yzKKj1be2ZoG2PM8dbaxc4nx9G+9YNEOsWBgua9N/Uyiip4aekeAL7YlMkXmzJbbJccHcx/bphGaWUN0aHe/YEnItIZEsIdd/O/3prN9weLvHp48DUvrHA9HhAXyoR+Mbx2/VSGJ0V4MKpj1yvSkaTe//5Gvrx7VqdfP6u4kpiQAH48ayB3/Xct/WI1N1a6hvYkQNcCzxlj6t/ZCoBr3BaRSDtYaxn7wOfcOieNjEbD1xLCA8lydrnHhQU26x3qFxtCoJ8PWzNLmmyfOTieAD8fAvyU/IiItEfjAjGLtuV4dQLU2IC4UACmDYz1cCTHLi0hnJOHJ7J4e06nXndLRjH3vruBZbvyGJ0cyXnjkzlvfHKnxiDSlsP25FhrV1prxwBjgDHW2rHW2lXuD02kdZsziiksr+YPH21if0FDAYOIRt3rvzp9qOtxgJ/jV/28cclcOa1/k3OdPz6Ze04dioiItF+Qvy//vmYycWGBrqHE3UHfbrZOzZi+UZRW1bqWe+gMj3+xlWW7HEPNf3vmiE67rkh7tXsom7W20Fpb6M5gRNqrfu2J4b0j2JHd0JsTEdTQqTkoIbzZ9vjwQGIbDXHbcP8p/OWiMYQHaVyyiMiRmjk4nrPH9mbR9hzeWZ3u6XCOmp+zjHd4oB/+vt1rlH/9Z15uaVWnXbN+WQqACf2iO+26Iu3Vvf6XS49QWlnDc4t2AeDn48PunIbqb40TmT6N3oDry13Hhwe6KsEZgxYqFRE5RqeM6AXAr9/Z4OFIjp6fryMB+t8tx3s4ko5X/5mXV9J5CVB9b9OL10zutGuKHAklQOJ1Xvh2NzklVfSJCmZdegFVtXWEOxOZxkPg6tc/AKiuqQOa9gCFaDE2EZFjNnlADNdNH0BZdecOs+oo5VW1VFTX8bN5QxgYH+bpcDpcrLNEeU5p6xVTO1pBeTWp8aGcMDi+064pciTaWgj1vLYOtNa+3fHhiLTupSW7mdAvhgVbshiXEsXI3pGuCm8XTuzLc4t3UVhezWd3ziS3pKpJSdaq2oYEKMg5H2hM36hOfw0iIt3R2JQo7CLYkV3idcUQ8sscPSMxXl7yujWxHugBKiyrdpXgFumK2hr/c2Yb+yygBEg6TXlVLb95dyPGOMZonzW2t6uHZ8agOH552lAslrnDEhmcGA6JjuMePGckoQG+vLpsLyv25BMbGkCQvy/PXjWRSf1jPPiKRES6j/o5l1syir02AWo8aqA7iQ8PxMfApoNFnXbNwvJq4rx8EVnp3tpaCPWHnRmISFvqCx1YC0UVNQxJDGfjAceb+awhCfj5+rRYaebKqf0AmDM0ga2ZJQQ5h73NHZ7YSZGLiHR/aQlhhAf6sWJPvteVOy6uqAEgopsWwwkN9OO0UUm8uGQ3MwfHM7MThqUVlFeRltD9hhNK99HWELgrrLUvG2Puamm/tfZR94Ul0mBPbinvrtnfZNvgxHAm9o/hQGEFF048/IdtVEgAkweox0dExB18fQwT+0ezbGeup0M5YvUJUHeuBnrvGcP5ZEMGb69K75wESEPgpItrawhcqPPP8DbaiLjdGU8soriypsm2wYnhRIcG8G9VmBER6RJG9onk663ZVNfWeVUp6eKKagDCg7pvVdCEiCBGJ0eSXeK+Qgi7c0q58rllPHPlRIoraposOSHS1bQ1BO4fzj/v77xwRBwyiyrYnlVCYkSgK/npFRFERlEFCeGBROuNVUSkS+kbHUKdhYzCCq9aTLShB6j7JkAAcWGB7M4tddv5311zgH155dzz9noABmoInHRhh/3fboxJBf4PmIqj+MES4E5r7U43xyY90I7sEkID/Dj7yUVkFjW9U/XydVO4/t8rvOqDVUSkp0iOdqy9ti+/zKvepxt6gLr3kK348EBW7Ml32/nLqh2J5Jp9BQDdsqS4dB/tud3xKvAkcK7z+SXAa8AUdwUlPdeJf/m6xe03zRpIWkIYT10xnmCt3yMi0uUkRzuSnvT8cg9HcmSKK2oI9PMhwM97hu0djbiwQPJKq9w2RHFXdtPepX6x3pMES8/TngQoxFr7UqPnLxtjfuqugKTnqq2zre6rH0s8tFdEZ4UjIiJHICkqiAA/HxZsyeKiiX09HU67rNidx6aM4m7f+wOOHiCAvNIqEiOCOvz8u3NLmTsskVNGJFJcUeOquirSFbUnAfrYGPML4HUcQ+AuBj4yxsQAWGvz3BifdDMV1bX8ff52aq3lp6cMdW1/ZdkevtqU1epxlTV1nRGeiIgcJX9fH26cmcoTX21nZ3YJqV18CFRWcQUXPL0EgNS40MO09n6p8Y7X+PXWbLckqAcLKjhuYBwXeknyKz1bexKgi5x//uiQ7ZfgSIhSOzQi6baeWrCDP32y2fW8PgHaX1DOr97Z0OaxU1TCWkSky7twYl+e+Go787dkd/kE6I0V6a7HYd28AALAtNRYhvYK5z/f7evwBKiksobiyhq39CyJuMNh/8dbawd0RiDS/TVOfgCenL+dgfFhbNhf2OZxH9w6nZF9vGtlcRGRnqhvTAgpMSGs2pPPtdO79teHr7dkux4PT+r+w6uNMUxNjeW/K/ZRW2fx9TEddu6MwgoAkiKVAIl3aE8VuKta2m6t/XfHhyPdVckh6/gAPPLpFgJ8faiqbRjeFuTvQ0V10+FuqiQjIuI9ekUGuXW9mY5QVlXDqr35/OiEVE4YFM+U1FhPh9QphveOoKyqlj25pR3aQ5dZ5EiA1AMk3qI9ZUAmNfqZAdwHnOXGmKQb2nywqMXt9clPZLA/r98wlcU/n9OsTXCAJlKKiHiLuLAAcrt4ArTpYDE1dZYJKdEclxbXob0hXdkwZyGhLRnFHXbOmto6Djp7gHqpB0i8RHuGwN3a+LkxJgpHQQSRFlVU1wK4KsBYa3lrVXqLbWNCA8grrWLygBimOu/A9YkKZn9BOVcf158ThyV0TtAiItIhYkMDyS3N9XQYbfr+gGPo9YgeNry6fq2m+oTlWFz0jyUE+fuyYnceZVW1+BjoHaUESLzD0cz6KwW69sBe8agZD8+nrLKGjQ/MA+Chjzfz2vJ9JEUGce8ZwwkL8uMPH23m9FG9uGBCX5btyuWUEb1cx39650xqauuICgnw1EsQEZGjFBsWQEFZtdvWm+kI6/cXEhXiT+8e1mMRFeJPoJ8PGUXHlgC9tHQPy3c1LQJ81pjeBPppxIZ4h/bMAXofR7U3cAyZGw78151BiXfLLnYMfVi6M5dfvLWO3blljOkbxYs/nORKaj6+Pd7V/uyxfZocHxbY/avxiIh0V7FhjvVmbnp5Jc/+YJKHo2muuKKaj9dnMG1gLMb0jKFv9YwxJEUGsS+vjMyiiqOas1NaWcNv/tdQufXyKSlcOjlFC5+KV2nPN80/N3pcA+yx1rY8nkl6PGsbFjP9ems2u3PLAJg9JF49OiIiPUCIc/jzF22s7eZJ1724guLKGqYPivN0KB4RGezPxxsy+HhDBmvuPemIP5sb9x4lhAdy24mDVPxAvE6rfdPGmCBjzB3AhcBQYLG1drGSH2lLbmmV6/FeZ/IDKo0pItJTzB7aMHez8U2xrmJLpqMAwPnjkz0ciWeEB/m7Hi/ZceRzteorvg3tFc6Se05U8iNeqa3BuS8CE4H1wKnAXzolIvFq6fnlrscfrj/IsKQIzhrTm3kjkzwYlYiIdJaY0ADuOdWx0HVpVa2Ho2nq3TX7KSir5ubZAwntocOt/3DuKP59zWTCAv246ZVVPDl/+xEdX58APXn5+B5TPU+6n7YSoOHW2iustf8ALsBRAlukTX/9cluT51HB/jxx6Tgig/1bOUJERLqb+HDHPKCz/rrIw5E0qKuz3P76GgCSIoM9G4wHpcSGMHNwPBP6RQOONfkyjqAqXEahY55vL/X8iBdrKwGqrn9grW2+iuVhOIfQLTfGrDXGbDTG3N9CmxRjzHxjzGpjzDpjzGlHeh3pGkora9iZXcKCrdn0bzQR8s6TBnswKhER8YT6BGhnTmmXGQaXWdzwJV9f3h1D2Oo9/+2udh2zN7eMBVuyCA/y67E9aNI9tPXbO8YYU796pQGCnc8NYK21EYc5dyUwx1pbYozxBxYZYz621i5t1ObXwH+ttU8ZY4YDHwH9j+qViEf98p31vLvmAABzhiby3OJdDO0VzuQBMR6OTEREOltUcMPE+vLqWkICPP9leVd2qetxckzP7QGqd9uJg0iKDGL57jxeXbqXW2anNZkf1JLznlpMTkkVM3poAQnpPlrtAbLW+lprI5w/4dZav0aPD5f8YB1KnE/9nT+H3gayQP25IoEDR/EaxIP25JZy2T+XupIfgLEpUQCaGCki0kMNTQrHzzk/pKTiiAeRuMWuXEcC9PK1Uxja67BfY7q90EA/rj5+ADedkEZxZQ1vr9p/2GNyShyFjs45ZPkKEW/j1lsyxhhfYCWQBjxprV12SJP7gM+MMbcCocBcd8YjHaOyppZfvbOBRdtyKKuqocj54Xbt9AGcPbY3w5IiuGJqCjfPTvNwpCIi4gn+vj785aIx3P76GooqakjoAvnGruxSAv18OG5grKdD6VJGJUeSlhDGX7/azoxBcfz2vY389swRpCWENWvbJyqY1PhQzhuvBEi8m1uXaLbW1lprxwLJwGRjzMhDmlwKvGCtTQZOA14yxjSLyRhzgzFmhTFmRXZ2tjtDlnZYl17ImyvTySiqcCU/AKeNSmJ0chT+vj787pxRPXqSqYhITxfhHE5VXFF9mJadY3duKf1jQ/FR5bJmZg6KJ6ekkjl/+ZpvtuXw9Nc7WmxXUFbF4MTwHreArHQ/bk2A6llrC4D5wLxDdl0L/NfZZgkQBDQbWGqtfcZaO9FaOzE+Pt7N0crh7Gm0vk9jo/pEdnIkIiLSVYUFOQaZ/OqdDdTWeb4Qws6cUgbEhXo6jC7p5tkDiQ5pmP9TVN48aa2qqaO0qrZJOxFv5bYEyBgTb4yJcj4OBk4CNh/SbC9worPNMBwJkLp4uri9uaX4GBiW1DCmYUhiOAF+nZJPi4iIFwh3JkDfHyzi+wNFh2ntXpU1tezLK2NAvBKglsSGBfLeLdM5aXgivSODWJde6NpXW2epqK6loNwx/ycyJKC104h4DXd+Y00C5htj1gHfAZ9baz8wxjxgjDnL2eZu4HpjzFrgNeBq21XqZUqLsooreOKr7USFBHDhBMcq2it/PZdP7tAyUSIi0qBxRbEtmcUejATeWrmf6lqr+T9t6BsTwj+vmsi1M1LJKKogy7ng6a2vrWLobz6hoMzRK6QeIOkO3FYEwVq7DhjXwvZ7Gz3+HjjeXTFIx/tsYyYAE/tF88Pj+3PZlBSC/H09HJWIiHQ19T1AABsPFHKB86aZJyzcmk1KTAjT01S++XDGJDuGs//kzXX8/pyRfLQ+A4D9+eVA0xLnIt5KY5akRen5ZS0uXrclo5jwID/+ceUEjDFKfkREpEVhjdb+aW3uaGfZX1BO/7hQTd5vhxG9HQnQwq3ZzHh4vmv7419sxRjoqzWUpBtQAiTNbNhfyPQ/zef17/bx3KJdfLs9x7VvS0YxQ1QBRkREDsPHx/Dq9VMY2SeC3JJKj8Xxt6+2sX5/IX2i9MW9PYIDfPn8zpnNtq9NL+SGGan0i9U8KvF+nl+aWbqcnTmOxeLeWb2f5bvyAIgK8cfXGHJLq7h+xgBPhiciIl7iuIFxDE4MZ9nOPI9cv6a2jj9/thWApEgtzt1egxLDufq4/rzw7W7XtqgQf35x6lDPBSXSgdQDJM1kFzvu1NUnPwAFZdXEhQUyOjmSW2YP8lRoIiLiZeLDAskuqWxxWLW7rW1UzcxX6/8ckfvOGsFTl4/HGLhtThrz756l0R/SbagHSJpJz28+VnveiF48feUED0QjIiLeLC4skKqaOoora1yLo3aWVXvyAbhiagpXTuvXqdfuDk4dlcT235+m5FG6HfUA9RC1dZa6OsvSnbnsLyhvtd3GA4U8v3g3Ab4+fHX3CdwyOw2Ay6emdFaoIiLSjcSFO6qG7cnpvEIIlTW1nPDIfH7/0SZ6Rwbxu3NGdXry1V0o+ZHuSD1APcR9723klWV7qLMwODGMz+48ocn+vy/Yzsb9RUxNjQHgoknJpMaHcedJg5kzLIHxKdGeCFtERLxccnQIAD97ax0f3945a8ZlFla6Ks8Nd1Y1ExGppwSoBygsr+alpXtcz7dmlvD+2gOkxIQwODEcf1/Dw59sAeDD9QfxMXD/WSMBx50fJT8iInK0JvaLZvKAGHZml3baNbOKHYt4psaFcp0K94jIIZQAdVPWWmrrLJszijnjr4sA+MuFY3h20S42HSzi1tdWt3psnVWXt4iIdAxjDMcNjGX5rjyqauoI8HP/6Pv6Yj5/vWyca10bEZF6SoC6qVteXc1Xm7OY4hzSBjBjUByT+scw85H5bRwpIiLSsepLUGcVV7iGxLlTtnPdoYRwlb4WkeZUBMGL1dZZSiprWtz34fqDlFfXUlhe7doWHx5ISmwIj140psVj7j1juFviFBGRni0xwpGIZBRWdMr1sooq8TEQExrQKdcTEe+iBMgLWWtZvTef3763gZG//ZSa2rpW267eW+B6XF+/f1L/GMICm3b+jUmO5JrpjnHSs4bEd3zQIiLSYyVFBgO0WYW0I2UXVxIbFqjh3CLSIg2B80IvLd3Dve9udD3fkV1KWkKY642+qKK62TEPnj3C9bhvTAjr7zuZM/+2iA37i1j567lEhTjukq2/72QC/Xzd/ApERKQn6R8XQligH0t25HL22D5uv97BogrXsDsRkUMpAfJCzy/efcjzXbz+3T7eumka41OiXZV2nrh0HGv2FjBnaALTB8U1OcYYw79+MImVe/KJDQt0bQ/XOgkiItLBAv18mTM0gU82ZnDvmcMJCXDv148DBeWkxYe59Roi4r2UAHmZ9PwyduU0LSX6+nf7ADj/qSX8aGYqPs6eoKkDYjhrTO9Wz5UYEcRpo5LcF6yIiIjTD47rx3trD/Df7/Zx9fHuK01treVAQTkzB2k4t4i0TAmQF6iureOd1fs5a0xvvtyU1WTf6aOSyCqu4Lvd+QD8Y+FOAIYnRZAQoe5/ERHpGib0iyEpMoh1+wvdep2i8hrKqmrpHaXPQBFpmRIgL7BgSzY/e3MdK3bnYTDEhAaQV1oFwG0nDmJIr3DGPfAZ+WUNc39UyEBERLqa/rGhbDpYTGVNrdvmm6YXlAENhRdERA6lKnBeYF+e4838ndX72ZRRxLCkcNe+QQmOMc4f3z6zyTHHDWw650dERMTT+sYEs+lgEWPu/4yP1x90yzV25zg+M/vHuX+9IRHxTuoB8gK7cx1zfqprLevSC/nRzFRuP3Ew2cWVrvk+vSKDmP+TWXy1OYvwQD+OT4v1ZMgiIiLNhDqXYKioruOmV1ax64+nuZZoOBaPfLqZrzZn8/HtM9iZXQLAgLjQYz6viHRPSoC6uPzSKv69ZA8Bvj5UOdf7uW5GKvHhgc3aDogL5drp7ptYKiIicixunp1GVnElH65z9P4cLKygd9SxD1V7cv4OADZnFLErp5SkyCC3V5oTEe+ld4cu7v11BwC4eFJfpqTGMD4lusXkR0REpKuLCwvkycvGc+GELK5+/jv+/NkWBsaHcdMJA10jGo5UXZ11PX5h8W725JXRP1a9PyLSOiVAXcRbK9PJL6viuhmpALy+fC+ZRZVU19bhY+D+s0Yc9YeDiIhIVzJjUDyjkyN5e9V+AMb1jeK4tKObu7ov3zHnJ8jfhzdWphMbGsCk/jEdFquIdD9KgLqIu99YC+BKgH7x9nrXvsSIQCU/IiLSbfj6GM4d14d16Y6S2NuzS446ATpYWAHAD48fwFMLdpBVXElEsBb1FpHWqQqcF/D31T+TiIh0L3FhDcO56xOho1FU7lgC4oTB8fg6bxZGBOv+roi0Tt+su5jaOku+c42fegWN1vcRERHpDmLDAlyPDxaWH/V5iipqAOgdGUxsqOOckeoBEpE2KAFygw37C7ngqW/JLak8bNuiimpeWbbH9fxvX23nxpdXAjBvRC8ASipr3BOoiIiIh8Q36gHKLalqo2XrLnjqW+59dwPg6PUJDnAsrqoESETaogSog1hr+fz7THbllHLGXxexYk8+a/YVHPa4t1am86t3NrieP/bFVpbtyiMxIpA/XzQGgFtmp7krbBEREY9oPAQut/TIE6C6OsuKPfmUVdUCEBboR7C/IwGKCFICJCKt0yDZDrJhfxHX/3tFk237Cw7fpb8nt6zZtocvGM201FjCAv3Y/dDpHRajiIhIV9G4lya/tApr7REtilpQ3nR4uJ+vD0H1CZB6gESkDeoB6iDr9he4HvsYCPDzIT3/8AlQS20unJBM35iQjgxPRESkS/HxMbx/y3RunZNGTZ2lqLwGa+3hD3TKK20+zLy+B0h1U0WkLUqAOsh3u/Jcj/18fEiODmZfXvPenUOl55cxa0g8D5w9AoBLJvU9ojtgIiIi3mpUciSp8Y5FSyf9/gtm/3lBk4VN29LSvKGLJ/UFYECcFkIVkdZpCFwH+HprNv9bc6DJtuTokMMOgauurSM9v5wpA2K4alp/jhsYS2pcmDtDFRER6VLq5wJV1daxO7eM9PxyUmIPPwoir4V5Q+eM68Ppo5O0fISItEnvEB3g6QU7SIoM4p9XTXRtS4oIIsO5OFu9hVuz+WxjBuCo/nb8Q19RUlnDrCEJAKQlhGvBUxER6VGmpsZy7rg+DE+KAODNVensbWF+bD1rLQu3ZpPjrLQ6d1git504yLVfyY+IHI56gI5RUUU1y3blcvPsNPpEBbu2J0YEklNSSU1tHTV1lozCCq56bjkAT142Hn9fQ1ZxJaeN6sWsIfGeCl9ERMSj/H19eOzisVRU1zL0N5/wxJfb+GxjBp/cMbPF9v9bs587/7OW6BB/jIG/Xz6eAD8lPSLSfkqAjtL2rGJ+9uY6JvaPoc7CtIGxxDVa1C0xMog6CwcLK7j6+eXsyC517bv51VXMHhJPgJ8Pj140VnN+RESkxwvy9+UP547il++sb3F4G0BNbR3PfrMLgPyyamY5P0tFRI6EEqCjFB8WxKq9BazaW0CfqGDGp0Tj6xy+dvroJHpFBAEw4+H5LR4/f0s241OiXCU7RUREerrLpqSwZl8+X2/NBuCd1emM7B1JcnQIFdW1LNyWzcYDRQ3tJ6d4KlQR8WJKgI5SZIg/t584iI83HOSV66a6Epml95xITGgAWzOLXW0n949h+e68ZucY0iu80+IVERHxBvHhgeSUVFFdW8ed/1kLQHigH9GhAZw3vk+TttMHxXkiRBHxckqAjsGdJw3mjrmDmgxh6xXp6PkZlhTBT08ZwpjkKCb0i+bNVensyi7lucW7XG0Hxqvim4iISGMJ4UHU1lk2H2y4kVhcWUNxZQ07sktJCA/k4QtGsy+vjJAAfY0RkSPntncOY0wQsBAIdF7nTWvtb1todxFwH2CBtdbay9wVkzu0Nn/H18dw8+w01/Mrp/ajts42SYDqS3+KiIiIQ3y447Pxyfnbm+17f+0BRvSOcFVPFRE5Gu68dVIJzLHWlhhj/IFFxpiPrbVL6xsYYwYB9wDHW2vzjTHd+h3N18fw3a/mEhzgyzur0jljdJKnQxIREelS0hIcoyM+cS4bcahdOaUtbhcRaS+3lU6xDiXOp/7On0OXd74eeNJam+88Jstd8XQV8eGBhAX6ceW0/vhprQIREZEmBieG89FtM5psmzIghqevmADAWWN6eyIsEelGjLWH5iQdeHJjfIGVQBqOROfnh+z/H7AVOB7wBe6z1n7SwnluAG4ASElJmbBnzx63xSwiIiKet3RnLt/tyuMHx/fH38eH4ABfCsurCQnw1WKnInJYxpiV1tqJLe5zZwLUKIAo4B3gVmvthkbbPwCqgYuAZBxzhkZZawtaO9fEiRPtihUr3BqviIiIiIh4r7YSoE65heJMaOYD8w7ZlQ68Z62tttbuwtEbNKgzYhIRERERkZ7HbQmQMSbe2fODMSYYOAnYfEiz/wGznG3igMHATnfFJCIiIiIiPZs7q8AlAS865wH5AP+11n5gjHkAWGGtfQ/4FDjZGPM9UAv81Fqb68aYRERERESkB+uUOUAdSXOARERERESkLR6fAyQiIiIiItIVKAESEREREZEeQwmQiIiIiIj0GEqARERERESkx/C6IgjGmGxgj6fjaEMckOPpIMTr6PdGjoZ+b+RI6XdGjoZ+b+RoePr3pp+1Nr6lHV6XAHV1xpgVrVWcEGmNfm/kaOj3Ro6UfmfkaOj3Ro5GV/690RA4ERERERHpMZQAiYiIiIhIj6EEqOM94+kAxCvp90aOhn5v5Ejpd0aOhn5v5Gh02d8bzQESEREREZEeQz1AIiIiIiLSYygBEhERERGRHkMJkBsYY241xmw2xmw0xjzs6XjEexhj7jbGWGNMnKdjka7PGPOI871mnTHmHWNMlKdjkq7LGDPPGLPFGLPdGPMLT8cjXZ8xpq8xZr4x5nvnd5rbPR2TeA9jjK8xZrUx5gNPx3IoJUAdzBgzGzgbGGOtHQH82cMhiZcwxvQFTgb2ejoW8RqfAyOttaOBrcA9Ho5HuihjjC/wJHAqMBy41Bgz3LNRiReoAe621g4HpgI36/dGjsDtwCZPB9ESJUAd7ybgIWttJYC1NsvD8Yj3eAz4GaDKJNIu1trPrLU1zqdLgWRPxiNd2mRgu7V2p7W2Cngdx806kVZZaw9aa1c5Hxfj+DLbx7NRiTcwxiQDpwPPejqWligB6niDgRnGmGXGmK+NMZM8HZB0fcaYs4H91tq1no5FvNY1wMeeDkK6rD7AvkbP09EXWTkCxpj+wDhgmYdDEe/wOI6bunUejqNFfp4OwBsZY74AerWw61c4/k5jcHQVTwL+a4xJtao33uMd5vfmlziGv4k00dbvjbX2XWebX+EYqvJKZ8YmIj2DMSYMeAu4w1pb5Ol4pGszxpwBZFlrVxpjZnk4nBYpAToK1tq5re0zxtwEvO1MeJYbY+qAOCC7s+KTrqm13xtjzChgALDWGAOOYUyrjDGTrbUZnRiidEFtvd8AGGOuBs4ATtSNFmnDfqBvo+fJzm0ibTLG+ONIfl6x1r7t6XjEKxwPnGWMOQ0IAiKMMS9ba6/wcFwuWgi1gxljbgR6W2vvNcYMBr4EUvTFRNrLGLMbmGitzfF0LNK1GWPmAY8CJ1hrdZNFWmWM8cNRKONEHInPd8Bl1tqNHg1MujTjuCv3IpBnrb3Dw+GIF3L2AP3EWnuGh0NpQnOAOt5zQKoxZgOOSaY/UPIjIm7yNyAc+NwYs8YY87SnA5KuyVks4xbgUxwT2f+r5Efa4XjgSmCO8z1mjfOuvohXUw+QiIiIiIj0GOoBEhERERGRHkMJkIiIiIiI9BhKgEREREREpMdQAiQiIiIiIj2GEiARERERETlmxpj7jDH7D1c10Biz2xiz3tlmRaPt/2l07G5jzJp2XjfCGJNujPlbe9orARIREYwxtc4PnA3GmDeMMSFHeZ4bjTFXOR+/YIy5oIU2Vxtjejd6/qwxZvjRR3/YmH7prnO3cU1jjPnKGBPRRpvXjTGDOjMuEZGOYoyZZYx5oYVdj1lrxzp/PmrjFLOdbSbWb7DWXlx/LI4FeNu7+O6DwML2xq4ESEREAMqdHzojgSrgxqM5ibX2aWvtvw/T7GrAlQBZa6+z1n5/NNdrp05PgIDTgLXW2qI22jwF/KyT4hER8RrORXgvAl5zPvc1xjxijPnOGLPOGPOjRm0nAInAZ+09vxIgERE51DdAmvPu3gf1G40xfzPGXO18vNsY87BzCMNyY0yac/t9xpiftHZiZ4/QROAVZ49TsDFmgTFmonN/ifNDbqMx5gtjzGTn/p3GmLOcbVr8IDTGJBljFjbqyZphjHkICHZue8XZ7n/GmJXOa9zQKLb2XPtqY8y7zu3bjDG/beWlXg686zwm1BjzoTFmrTOuixv9Pc81xvgd4b+PiEhXdovzvfk5Y0x0K20s8JnzvfiGFvbPADKttducz68FCq21k4BJwPXGmAHGGB/gL0CrnzstUQIkIiIuzi/jpwLr29G80Fo7Cvgb8Hh7zm+tfRNYAVzu7HEqP6RJKPCVtXYEUAz8DjgJOBd4wNmmxQ9C4DLgU+fQiTHAGmvtL2jo3brcefw11toJOBKx24wxsUdwbYDJwPnAaODC+uTtEMcDK52P5wEHrLVjnD1snzj/LuqA7c5YRUS8gjFmmXNuzrPAWY3m7JyCo2d7IDAWOIgjOWnJdGvteByfNzcbY2Yesv9SnL0/TicDVzmvuwyIBQYBPwY+stamH8lr0F0nEREBZy+J8/E3wL+A4w5zzGuN/nysg+Kowpkg4EjCKq211caY9UB/5/aTgdGN5hdF4vgg/A54zhjjD/zPWrumlWvcZow51/m4r/PY3HZeG+Bza20ugDHmbWA6jqSusRhrbXGjc/3FGPMn4ANr7TeN2mXhGA64EhERL2CtnQKOOUDA1dbaq1tqZ4z5J/BBS/ustfudf2YZY97BcWNpofM4P+A8YELj0wG3Wms/PeQaPwBmGGN+DIQBAcaYEufNr1YpARIREXD2kjTeYIypoelIgaBDjrGtPD4W1dba+nPVAZXg6C1pNFSsxQ9CAOddxNOBF4wxjx46H8n5gT0XmGatLTPGLKDhdbXn2tD8tbb02muMMT7W2jpr7VZjzHgc84J+Z4z50lpb36MUBBzaCyYi4pWMMUnW2oPOp+cCG1poEwr4WGuLnY9Ppmkv+1xg8yG9Op8CNxljvnLemBoM7G/Us49ziPbEwyU/oCFwIiLSuj3AcGNMoDEmCjjxkP0XN/pzyRGctxgIP4a46j8I/QGMMYOd82z64Rgz/k8cQzPGO9tX17fF0VuU70x+hgJTj+L6JxljYowxwcA5wOIW2mwBUp3x9QbKrLUvA480igtgMC18QRAR8VL1c0PXAbOBO8HxPmiMqa8IlwgsMsasBZYDH1prP2l0jktoOvwNHO/p3wOrjDEbgH9wDB056gESEZEWWWv3GWP+i+ML+i5g9SFNop0fcpU4xmu31wvA08aYcmDaUYT2LI4haauMMQbIxpGIzAJ+aoypBkqAq5ztnwHWGWNWAdcANxpjNuFIUpYexfWX4yjPmgy8bK09dPgbwIfOeLYDo4BHjDF1QDVwE4AxJhFHz1vGUcQgIuJR1toFwIJDtl3ZStsDOHrBsdbupI25jy0NqXPOmfwlbVT1tNa+gOPz5bBMQ2+/iIhI+xhjduMYapDj6Vg6U6MhFrccpl0S8G9r7UlttLkTKLLW/qtjoxQRkbZoCJyIiEgHc46B/6dpYyFUoAB4sXMiEhGReuoBEhERERGRHkM9QCIiIiIi0mMoARIRERERkR5DCZCIiIiIiPQYSoBERERERKTHUAIkIiIiIiI9xv8DT8e/R1qEHzMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from pyplr.utils import unpack_data_pandas\n",
    "\n",
    "data = unpack_data_pandas(data, cols=['timestamp','diameter_3d'])\n",
    "ax = data['diameter_3d'].plot(figsize=(14,4))\n",
    "ax.set_ylabel('Pupil diameter (mm)')\n",
    "ax.set_xlabel('Pupil timestamp (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Timestamping a light stimulus\n",
    "-----------------------------\n",
    "\n",
    "The obvious way to timestamp a light stimulus would be to control the light source programatically and send an annotation as close as possible to when we change the status of the light:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```Python\n",
    "p = PupilCore()\n",
    "p.command('R')\n",
    "sleep(5.)\n",
    "annotation = p.new_annotation('LIGHT_ON')\n",
    "my_light.on() # turn hypothetical light source on\n",
    "p.send_annotation(annotation)\n",
    "sleep(1.)\n",
    "my_light.off() # now turn it off\n",
    "sleep(5.)\n",
    "p.command('r')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But in reality this will be problematic as the light source will have a latency of its own, which is difficult to reference. In fact, our own light source takes commands via generic HTTP requests and has a variable response time on the order of a few hundred milliseconds. Given that we may want to calculate latency to the onset of pupil constriction after a light stimulus, which is typically around 200-300 ms, this variable latency is far from ideal.\n",
    "\n",
    "To solve the issue and to make it easy to integrate *PyPlr* and Pupil Core with any light source, we developed a method called `.light_stamper(...)`. This method uses real-time data from the forward facing World Camera to timestamp light onsets based on the average RGB value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Waiting for a light to stamp...\n",
      "Light stamped on frame.world at 55986.582772\n",
      "(True, 55986.582772)\n"
     ]
    }
   ],
   "source": [
    "p = PupilCore()\n",
    "\n",
    "p.annotation_capture_plugin(should='start') \n",
    "\n",
    "p.notify({'subject':'frame_publishing.set_format',\n",
    "          'format':'bgr'})\n",
    "\n",
    "p.command('R our_recording')\n",
    "\n",
    "sleep(5.)\n",
    "\n",
    "annotation = p.new_annotation(label='LIGHT_ON')\n",
    "\n",
    "timeout = 10\n",
    "\n",
    "lst_future = p.light_stamper(\n",
    "    annotation=annotation,\n",
    "    threshold=15,\n",
    "    timeout=timeout,\n",
    "    topic='frame.world')\n",
    "\n",
    "sleep(timeout)\n",
    "\n",
    "p.command('r')\n",
    "\n",
    "print(lst_future.result())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like `.pupil_grabber(...)`, this method runs in its own thread with `concurrent.futures` so it does not block the flow of execution. The underlying algorithm simply keeps track of the two most recent frames from the World Camera and sends an annotation with the timestamp linked to the first frame where the average RGB value difference exceeds a given threshold. To work properly, the `.light_stamper(...)` requires a suitable stimulus geometry, an appropriately tuned threshold value, and the following settings in Pupil Capture:\n",
    "\n",
    "* **Auto Exposure Mode** of the relevant camera must be set to **Manual** \n",
    "* **Frame Publisher Format** must be set to **BGR**\n",
    "* **Annotation Captre Plugin** must be enabled\n",
    "\n",
    "In our testing, `.light_stamper(...)` flawlessly captures the first frame where a light becomes visible, as verified using Pupil Player and the Annotation Player plugin. Timestamping accuracy therefore is limited only by frame rate and [how well the Pupil software is able to synchronise the clocks of the eye and world cameras](06b_pupil_core_timing_analysis.ipynb). \n",
    "\n",
    "`PupilCore(...)` also has some other cool features, like a `.fixation_trigger(...)` method that allows you to wait for a fixation that satisfies criteria relating to dispersion, duration and location of gaze samples. Check out the code for more insights."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}