{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Demo: saccade detection\n", "----\n", "\n", "This is a demonstration of how we detected saccades.\n", "\n", "Dependencies:\n", "\n", "- Python 3.5+\n", "- Numpy\n", "- Scipy\n", "- Pandas\n", "- Matplotlib\n", "- sliding1d\n", "\n", "**NOTE**: you need to get content of the `04_formatted` data repository populated beforehand, via:\n", "\n", "```\n", "$ cd ../01_data/04_formatted\n", "$ gin init\n", "$ find . -name \"*.csv\" | xargs gin get-content\n", "```" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "from collections import namedtuple\n", "import math\n", "import numpy as np\n", "import scipy.stats as stats\n", "from scipy.signal import find_peaks\n", "import matplotlib.pyplot as plt\n", "from matplotlib.gridspec import GridSpec\n", "from matplotlib.patches import Rectangle, Polygon\n", "import sliding1d as sliding\n", "\n", "import datareader" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load trials from the dataset:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "trials = datareader.load_trials(\"../01_data/04_formatted\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Browsing original eye-position data\n", "\n", "We use the `left_pupil_normalized_position` and the `right_pupil_normalized_position` columns of the `tracking` table." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "subject = 'MLA007518'\n", "session = 'session2019-09-30-001'\n", "xlim = (9, 17) # in seconds\n", "\n", "in_session = [trial for trial in trials if (trial.subject == subject and trial.session == session)]\n", "trial = in_session[2]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABBM0lEQVR4nO3dd5hU5dn48e8zbXtl2WXZpfciICAiCqJYsKHYS5REjTFqNCbq66ux5E00Gk38JdGo2NDYYu+FYkERpfe69AW2w/adnfL8/njO9h3YMrAzu/fnuvaamTPnnLlnB+be537KUVprhBBCiObYOjoAIYQQoUuShBBCiIAkSQghhAhIkoQQQoiAJEkIIYQIyNHRAbRFSkqK7tu3b0eHIYQQYWX58uUFWuvurTkmLJNE3759WbZsWUeHIYQQYUUptau1x0i5SQghRECSJIQQQgQkSUIIIURAkiSEEEIEJElCCCFEQJIkhBBCBCRJQgghREBhOU9CCCE6ta3zYM8ScMXA+GshMr7DQpEkIYQQoeTLe2Hxk3WP170Dv/gCImIDH7PmLdj2FdgcMOVOSOpT99yyl2DpC20OR8pNQggRKjZ9ahLEMZfAfQVw+euQsxb+fQIEukDc+g/gvV9C1nxY/QZ893jD51e/AWW5kNi7TSFJkhBCiCNh/2p47RJ4OAOemghr3obivXBwD2z+HMryG+6ftxH++zNIGQLTHgC7E4aeYxJG8W4o2g5+H1QUwYGdsH8N/DQb3r0eMsbBb9fCmCvN61QeqDtvYRYMOQuueL1Nb0PKTUIIEUz7VpnWwNq3we6CYTNg3bvw3vUN94tMhDMfhog4cJfAt3+FyAT4xWcQk1K33+Q7zLk++o1JMMW7G54ncwJc9RY4o+C4X8KKV0x5acodJllUFEK3AW1+O5IkhBAiGPx+yN8Ez50K2gcn/Q7G/8KUec74MxRtgz0/mX2T+sJHt8KHN9UdH5sGV73bMEEApA41rYm1b5vHpz0IMammMzs2DTKPA6XMc+mjzOMtX5okUbjdbO82sM1vS5KEEEK0l9bwn/Nhx0Lz+Nc/QNqIuufj081P35PqtvU/BYqzzRe8I8okE4er+fOPuNAkiWEz4KTbDx1L5gRY8qw5d84asy1ZWhLiCNpTVMG2/DJcDhsxLgcjMxKw21RHhyVEaPD74bWLTYI45lIYPqNhgggkOtn8tMSQs+CmnxqOWgpk+Pnw41PwhBWDzWFaLm0kSUIEpLXmq0153PL6Sio9vtrtF43N5PFLRpF9oJKZ/15EaZUXl93Gs9eMY9KAlEOcUYhORmv49HewbYH5C//U+8BmD/7rKGXKTi3R+3iY9THkrjcJImUwOCPb/NKSJEQTOcVVbC8o453l2by3Yi9D0uK4a/oQYiIcPLdwO++uyOanHYXYlKKwvJpfTOrHi4t2sDa7WJKE6FpWvwHLX4Lek2Dq/7YqQazbW8zuogoUMGlACgnRzuDF1W+K+QkCSRJdiM+vyS2pIj0hEqWaLxd9sS6HW15fgddvxmRfdXxv7j1nGNEu809lXJ8kXvh+B4u3FQKmVXHrtIG8+uMuiiqqj84bESIUuMvgy3ug9wnw80/B1vIZBUt2FHHps4trH980dQB3TW9hS+EokyTRBXy0eh8/bS9kUVYBOwsrOCYjgfNGpzOxfzdGZSbW7rcoq4AbX13OiJ7x/M/0oWQkRdE/JaZBQnHabdx48gBuPLlhR1hyjIuiMkkSoovwVsNLZ5khpqf/qVUJwu318b/vrSEzKYrZV4/n8tmLKa3yHsFg20eSRCNaa77ZnE96YiRDewReL0VrTXm1j1W7D7K/uBKXw0aEw0ak005itItuMS425ZSyfNcBusdFcOWE3kS5Wt4ULa3ycO/761i3r5i9Byo5vn83fj6pD6cOTQOgpMpDSaWHcreP/FI3eaVV5Ja48fj8TB6UwuC0OLbmlfHI5xv5cXsRDpsiLT6S04alsWxXEQ9/tgmAif2TGZIWh09rvttaQHKMi3d/PYlIZ+vqqskxLorKJUmITm7fSpj/R8hdB+X5MOEG6HVcq07x1Nfb2JZfzsvXTmB4z3iiXHY8Pv8RCrj9gpIklFLTgX8AduB5rfUjjZ5X1vNnAxXAz7XWK1pybCDFFR6+3JDDzoJyKqp9lLu9eHx+fnvaYPqmxDTc2ecBZW+S7QvL3OSUVHGg3ENhuZuPV+/npx2FtVl9dGYCsyb15YIxGdis0TzVXj+vLN7J3+dtoaLax+EoZfq2Hv5sIwO6xzB9ZDo3TR1Q+yVc5fGxcvdBlIIIh434KCfJ0S6+WJ/DR6v3MbRHHBeNy2TBxlyuf3kZV0/sQ26Jmy/W5wR8zb/P21J7Pz7SwR/OGcbVJ/QhwlH3xV9Q5uZfC7ayZOcB3lu5F4dNEWvt29oEASZJFEqSEJ3dp3eYmdQjLoCRF8OQ6a06fGtuKU9/k8UFY3py8uDuALgcNqq9nThJKKXswFPA6UA2sFQp9ZHWekO93c4CBlk/xwNPA8e38Nim9q8h4a8pbPVcyYv6PJJdPhJcmj2lmiiXnVmT+qI1OIp3k77kIWK3f2aOi0zAF5fBvjJNUZWmwmujkHi26x7M841nk60/pwxJIzMpGo/Pz9KdRTzw1mJWf7iEixM3ExUdz4/FiewptnFxXDqZp15En+7xDE+Px+31U+31U+nxcbCimqKSMvp5dzAmzcHmXdlsyavkzeJh/HPBVl5atIMpg7ujtWbxtkIOVHiafZsJUU4+vXUydpui4pxh3P7fVbzy4y6inXZ+dXJ/BqTEEh1hJyU2grT4SFLjIvD4/Hy8Zj8Vbi8psRFMHpxCalzTkQ0psRH88fyRrfikDy3aZaegzB208wkRkoqzzdIXM/7Z6kP9fs3d760lNsLBfecOr93usttwd/KWxAQgS2u9HUAp9SZwPlD/i/584BWttQZ+VEolKqXSgb4tOLaJAzoa0NzrfJ170tag8jdBtZfyyBieW3Em/16WzjT7Ss6zLUYDC+POIj61N+6SfCryd6N8HtJi7aS6NKO8+zinYim3Oj7AnzIE2ylPQS/zAXr3LIM5v8Hhq2TPge7EHqzkZ6oMnEAVsOhfZlja0HPNMLPB08Fu/UrnPwPf/x2AEdbPzBNu4dtTf8u7y7NZk30Qu01xwoBuzDw2k5gIO26vnwPl1Rys8GBTMKLefIRol4Nnrx7fog/k6oktGEsdZHabwucPsACZEJ2FtxIcbRtO+tqS3SzfdYC/XTKabrERtdtdDnvnbkkAGcCeeo+zMa2Fw+2T0cJjAVBK3QDcAJDUsx8rpv6d0QcXYC/eCaMuh9Rh2LfM57c73wPAZ4tg56DreblqMv/Z4sBfYJp100f0YNakPgzrU28SS0URbPwY28LH4KXpZnGtuB44Fv0DImPhsg/oljaOKq+myu4j0lduluXNXmZWbVzwR3OeiAToPtgkjKwF0HMsnP5/Zvr8on/C4ic5OWMsJ19xUet+w2HAZlP4Aq1SKURn4XW3ac5BTnEVj36+iZMGpnDh2IwGz3X6chPQ3FjKxt8WgfZpybFmo9azgdkA48eP12OnXgBc0GCfyEm3mAWwvFXYY1LoH5XEH4F7vD5Kq7zERzpxOZoZhRCdDONmmRbBW9fAvPvMdpsDLnoBek8kGoiuTf5RMOpS83Pmw1BZZBJG1nwzDX77t2YFx4m/hn6TzSFnPmxqme9cC5//j1lLpftQM01/xIWtGh3RQHkBfHwb7F0BngqzlkvKIOg/FZzRsPVL6DHK/OOOTjZrwDReGyYI7Erhl5aE6My0Bm9Vm1oS93+4Dq/fz0MzRzYZfu6yq06fJLKBXvUeZwL7WriPqwXHtk5iryabIhx2ImJb0Bkb0w1mfQQFW02CaMm0ebsDYlNh6NnmJ5C4NLj2S1j1GhRsMa+x9m0zEWfde3D5a3WLdHkqYdcimPegSUADp5nJOgNOgYh4c/z2b8w5ts41qzyOvMisIFm6H/Yuh02f1L32hg/r7s97wLR27BHmH7vDZd23biNiIakf9J4IvSYc/ndW82uQloTo7LxWn5sj4tD7NfLFuhzmbsjl7rOG0qdbTJPnXQ4bVZ7OnSSWAoOUUv2AvcDlwJWN9vkIuMXqczgeKNZa71dK5bfg2KPLZm/59PfWiukGJ95a99jrhnn3w0/PwEPpJsGV7IfqUisWh2lprHzVLP/bWHQ3M5Fn0m/Ml3oNn9esOOn3mUXDSvZCfIZZoXLp81B5EHxu8/pVxeCrNuO+fdbjygPmtW/8HlKHteit2ZTCH7r/zoVoP2+VuXVEtfiQkioP93+4juHp8Vx/Ur9m93HZbZRUduJ5Elprr1LqFuBLzDDWF7XW65VSN1rPPwN8hhn+moUZAvuLQx3b3pjChiMCpj9iFgPb85Ppxxh8pklSSf2gzySI72m+9HcvNhclqSyC7kNMCSnQGvF2h9mnRs39zPHm51C0Ni2VpybAgv8zV8ZSylzk5NWLoLrc7NfvZDj/qdqOersN6bgWnVttkmh5S+LRzzdRUObm+VnjcdibLyl3hT4JtNafYRJB/W3P1LuvgZtbemyXohSMvcb8BGJ3mL6Nmv6NIx1P9yGQMR42f2bWxk/qY0pYJfth9OXgLoU1b5rJRJe8BJEJUm4SnV9NknC2rCWxdGcRr/20m+tO6tdgZYPGXA471SE8BFYuXyqaN+tjmPx7M7N08xcQn2nGhp//JFz6srmgyrYF8KP5W8AmHdeis/O0vCVhlt5YS0ZiFL87ffAh93XZu0BLQnRCrmiYdj9MuhW0v2kH/mkPmBUwN30Mk38nLQnR+dWWmw4/uunpb7aRlVfGS784jpiIQ3/Nuhw23I2SxLdb8nno0w34NfRMjOJflx8b3FViW0FaEuLQohIDj/DqPxVy1sKcc3Dgkz4J0bnVjm46dJJYsqOIp77OYsbonpwyJPWwp3XZVZO1mxZuyWdHQTm9kqJYuCWf//y4s61Rt5skCdF2Z/zZtDb2/MSZe5+UcpPo3LyV5vYQSWJN9kEun72YaJeD+88bHnC/+kxLwoeu1xI/UFFNalwkL/1iAhmJUWzLL29X6O0h5SbRdjEppt+icDvHr3qVY/UwoHULngkRNrzWApZWn8T+4ko+W5tDlNPOzGMzeHdFNk99nUVcpJN3f30CKbEtGwXVv3ssVR4/a7KLGd0rEYCDFR4SrfJSr+QodhdVBP3ttJQkCdF+Z/wJVr3KMHZ0dCRCHDl+M5fhzWX72LZ6Az9sK2T9vhKgZka1Ji7SwQuzjmNgalyLT3vOqHQe/Gg9767Irk0SByqqSYp2AdA7OZr5G/Pw+3XtatRHk5SbRPtFJqBRxFHW0ZEIceRo02/wn5/28OqPu9lZUFcC8vo1vzl1ICvvO50J/Q6zSkMj8ZFOTh+exrvLs8k+UEFltY/1+0pqWxInDkyhqLyaH3cUBu+9tIK0JET72ey47bHEeyVJiE5Mm+vHnD4ind/+rK6s+uGqvewvruKXk/vXrtrcWnefNZSvN+Xx2JebySmuotrrp0e86fs4Y3gPYiMcvLt8b4dcQ16ShAiKKmc8CdXlHdYkFuKIs1oSdkfDr83zx2Q0t3erZCZFM3NsBq/+uBuAM0ekcdtpgwCIctm5cGwGr/64i+kje3D68LR2v15rSJIQQeF2xJNIGT6tsTW7uK8QYc5vWhJ225H52rz5lIEkRDmxK8V1k/sTF1k3L+K2aYP4YOVefvnKMi4cm8HozEQGpsZy4sAj37KQJCGCwu2MJ0EV4pcJdaKzsv5t2x2tv7xvS6QnRHHnmc0vLtotNoJv7jyFO95ezfsr9/Leir247DZW3H86sYeZrNde0nEtgqLaGU8C5bISrOi8dE1L4sgkicNJjnHx4s+PY+2DZ/LCrPFU+/wsyio44q8rSUIERbUzgQRVLktziE7L7zNDYBv3SRxtsREOpgzuTlyEg2825x3x15MkIYKipiXhC+HVLIVoD5/PakkcoXJTazjtNiYPTuHrTfkNZmofCdInIYLC40zAqXz43WWHv5pfC1V5fKzbW0x+qZvyah9xkQ6O7ZVIanzbLkQvRHv4/D6cgOMIdVy31tQhqXy2NoeN+0sZ3jP+iL1OaLxbEfa8LvOPVFcegKT2J4l5G3L50ycbmixHkBIbwXd3nUKUq+P/mhNdi88bOi0JgKmDuwPw1abcI5okpNwkgsLjSjB3Kg606zzb88v43Vur+OUry6io9tEzIZI/nDOMhXeewsMzj6GgzM3cDTkBjz9YUU1eaRVVHl+LXs/n1+QUV+GVMpk4DL+1LIfdHhp/W6fGRzKhXzJzfthFQZn7iL1OaLxbEfa8rkQAdFXbk8T8Dbnc+OpyvH7NtKGpPP2zcbgcdX/HXJbUi0c+38iirIJmJzA99XUWj325GTAX2BuVmcjsq8eRFqA89e2WfO58ezV5pW6SY1zcNHUA10/u3+b4RedW0yfh6OCO6/ruPXsYlzyzmMmPfs2LPz+OEwZ0C/prhM67FWHNE1FTbjrY5nN8vi6HhCgnb94wkQHdY5vM3LbbFKN7JbJxf2mTY3cVlvPU11mM7pXIxWMzyD5QybMLtzPjye8ZlZnImF6JHNsrkeE940mMdpFXUsUdb6+mstrHHWcM5pXFu3h87mZJEiKg2o5re2iUmwBG90rkg5tP5JoXf+J/31vDTVMHotFoDTUr908blhrwD6WWkCQhgsJntSRUO5LErsJyBqbGMigt8AqaPeIj2ZJblyS01vz7m208PnczdqX42yWjGZgaC0DflBjeX7mXLbmlzNuQC0BPRwlvjN3EngOVzKgoY9b0E+l97Fi0hr/N24Lb6yMiRGrOIrRoawhsKLUkAIb3jOfpn43j2peWcte7a5o8f9meXjx68ag2nz+03q0IW97IZPxa4dq3BPT1pt5To6II9vxkZqxmHgexpsOtyuNjW34Z1V4/ldU+NueUcs6o9EO+Tlp8JAVl1fj8GrtNsXzXAR77cjPj+iTxlwuPqU0QAFdM6M0VE3qjtWZ7QTnZBypZ+fp99FnzBn2Ak5zAgtfgu1hGD7kH6EtplZeIWEkSoqlQLDfVOK5vMovvmcaB8mpsNoVNgUJxx9urWbXnYLvOHXrvVoQnZxQ7dA8GbHoLvhsCU+4w23PXw8szoKJmZqiC0VfAzKf5/dur+XTN/tpT2G2KmcceerG01PgIfH5NUXk13eMieHbhdlwOGy9fO4FYhwafB5QdbHV9GUopBnSPZUD3WEYPrKIiuwefT/uSY9Mc9Ccb3v4FEzc9wlj1O0qrprb4YjGia/Fbazc5QqjcVF9shKPJEh1j+yTx5Fdbqaj2Eu1q29e9JAkRFHalmOW5m3l9/kPUV3+CxU+aZLB3hbk28DUfmR0XPwmrX4dJt5CVW8a4PknccupAXHYbPROj6JcSc8jX6ZUcDcAP2wo4cWAK8zbkcsuUPsR+/zB893fAKsSOvxbOfaLJ8YklWyBzJBcd19fakgk/exf/nJk873qcPRU/Bw4dg+iaapNECLYkAhmVkYBfw8b9JYzr07ah6eHzbkVIs9sgW3dn9/SXGbL/I1gyG378t3ly4k3Q/2Rzv8cx8MQIeP9GfAdv5JhxY1p0sfgaUwZ1Z3h6PHe9uYQH4j/hM9dShi/ZZZ4cdAb0mgDbv4VlL8Kp9zWc2LfhQ8hbD6MuaXjStOHkjriePsv+zJaSAiCp7b8I0Wn5a2dch8/X5qA0U37NyiuTJCE6ls3qg/A44+CEm2Dir/H7NXM35PDakj30en8tv5rSnz7dkuHSV9Bv/5w/63+xJuHVVr2O3ab487kD6f/6NSRWF5EVdyy+Yy/BnpgBY34GDpfp99j5HWz4wLQoAPI3w/u/hshEGHlR0/Mm9zbxF+0GBrXjNyE6q5q1m1xhlCQyk6JxOWxk5bX9gmDh825FSKu5IlfNUuHbCsq56dUVbM4txeWw8d3WAt5auoenrhrLmSNOZ8Xg25m47k/0yHkWvA+ZL/cWGuvKBl8RHHc9A8/5W9MdMsZDXE/45HbYtRiS+sD698EZCTcugvimneP2xF4AOEqy2/DuRVdQU25yhlGSsNsU/VNimLsht1Ut9vpkxrUIipo5DT5rcPZ/l+5he0EZf790NOv/eCYf3nwiEQ4bv3ljJY9+sYlLlg1hGcPpu3E2zL236Qm1hpJ9ULwXqssbPrdzobk96fbmg4mIhRu/g+EXwMaPYOHjZoTVRS80myAAiDGTkOye4ta+ddFF+P0+fFrhsIfXRbXG9UliV2EFVz7/U5uOD5+UKEKa3So3ea0kkVdSRY+ESC4cmwmYST9v3DCRGU8u4ulvttG3Wyzp1y+AH+43/RdRSXDKPeZkfj98eBOsfsM8dkRB/6lgs5sv+90/mJJS/CFGQsWkwKUvtzh+mzMKAOWtat0bF12G3+fHhw2XPbz+tv6/80dyxYTelLu9THy09cdLkhBB0bebGRG0cX8Jx/VNJq/UTWpcw1meozITmXv7FArK3AzrEU9SjAumPwLl+fDto1B50HRsr3sHtn8Dw8+HvpNNa+CgufYvdgdMvQeOazQXo51sLpMkbJIkRADa78OPDUeYJQm7TTEyI6HNx0uSEEHRKzmKjMQoftxeyFXH9+GHbYWcNbJHk/0Gp8UxuP6MarsDLnoeKg/AkmetbRFw1l9hwg0mEUz45RGPv6YlYfNWHvHXEuHJbyUJZ5iVm9pLkoQICqUUE/t3Y+6GHH4xZynAYec81LI74ZoP4cBOqC6D7kPNtqPI4YzAo+3YfNKSEAH4fWFZbmqvrvVuxRE1Y0xPIhw21mQf5JJxmfz+jCEtP1gpSO5nyk1HOUGAmaBdhQub78gtuSzCm9/vQ6PCrtzUXtKSEEFz8uDuLPvD6R0dRps4bDbKcEqfhAhIWy2JrlZualdKVEolK6XmKaW2WrfNTlVVSk1XSm1WSmUppe6ut/1BpdRepdQq6+fs9sQjRFvZbQo3Lux+SRKiedrvt5JE12pJtPfd3g0s0FoPAhZYjxtQStmBp4CzgOHAFUqp4fV2eUJrPcb6+ayd8QjRJnabokpLuUkEpq1ykySJ1jkfqBmM/jJwQTP7TACytNbbtdbVwJvWcUKEDJsyfRJ26bgWgVjlJrtNyk2tkaa13g9g3TY37zsD2FPvcba1rcYtSqk1SqkXA5WrhDjSlDLlJodfWhKieX7tR9O1EgS0IEkopeYrpdY189PS1kBzv1VrPWeeBgYAY4D9QDML8dTGcYNSaplSall+fn4LX1qIlnNLS0Icit+Hn9C8lsSRdNjRTVrr0wI9p5TKVUqla633K6XSgbxmdssGetV7nAnss86dW+9czwGfHCKO2cBsgPHjx+tA+wnRVm7lwuFv+2qZopPz+/CrrtUfAe0vN30EzLLuzwI+bGafpcAgpVQ/pZQLuNw6Diux1JgJrGtnPEK0WTUu7FJuEgHoLlpuau88iUeAt5RS1wG7gUsAlFI9gee11mdrrb1KqVuALwE78KLWer11/F+VUmMw5aedwK/aGY8QbVatIqRPQgTm96GVlJtaRWtdCExrZvs+4Ox6jz8Dmgxv1Vpf3Z7XFyKY3MqFU5KECKCrtiS6XoFNiAA80pIQh6K7ZktCkoQQlmppSYhD8fvxd8GvzK73joUIwK0iseMDb3VHhyJCkfaDjG4Souty28w1JaiWYbCiGX4fWpKEEF2XW1lX0vNUdGwgIiS5/BV47ZGH37GTkSQhhKW6tiUhSUI0FesrptyR2NFhHHWSJISw1CYJaUmIZsT5S6hyJnZ0GEedJAkhLNpmDW/0ezs2EBF6/H7iKcXt6nprkEqSEMKibdZlU32ejg1EhB53MQ78VHfBJCGXLxXCUjtRSloSor7cDbDuHQAqo9MPs3PnI0lCiBo1LQm/tCQ6vaId8OHNUJxtHl84G3pPbLhPaS588zAsnwNAjk4iJ3Xy0Y0zBEi5SQiLtll/M/mkJdGpVR6AOefArkWQOR5K98NLZ8E3j9TtU1EEz0+DFa9AvymUX/05J7ufwBkR3XFxdxBpSQhRoyZJSLmp8yrIgpfPNYnh0v/A8BlwcDfMux+++Qt43ZTEDSDum/tQlQfgZ+/CwGkcPFiJm6+IcnW9tZskSQhhUQ4pN3V6OxeaBHHibSZBACT2Rs98lj0Hquj9/d+JB9bTn8IpT3FS/1OxAU99nQVAv24xHRZ6R5EkIYTF6Ywwd6Ql0XlVFJnbqfc02Pzmilz+sP1KLo4/lqEZ3XijoD9b5lYyadtPnDIklXeWZXPp+EwmDUzpgKA7liQJISwOp/RJdHqVB8AZA86Gy2ss2JhHZrc4/vL7u7DZFFd5/TwxfwtPf7ONH7YVEu2y85tTB3VQ0B1LkoQQFkdtS0LKTZ1WWR7Edq99WFRezVNfZzF/Yy6XjMvEZjMXFXI5bPzP9KHcNHUAfg0RDhuRzq7XHwGSJISo5aztk5CWRKdVlguxabi9Pv42dwvvrcimqLyac45J584zhzTZPS7S2QFBhhZJEkJYnC6XuSMzrjuvsjxIGci/FmQxe+F2ju+XzP+7bBAnDep6fQ0tJUlCCIvLKjf5vB66ZmGhCyjLJSt6NE+uzOLsY3rw76vGdXREIU8m0wlhcVktCY9HrkzXKXmrobKIT3f46NstmjvOaFpeEk1JkhDCUlNu8nqk3NQplecDsN+XwO/OGEL/7rEdHFB4kHKTEBZXRE25yd22ExRtB3cZOCIgZTAoFcToRHvl5e4lFSgijrG9Ezs6nLAhSUIIS4TTjGTxedvQkshZB8+cWPf4iv/CkOlBikwEw568IlKByyYOJjOp663B1FZSbhLCEuVyUK3tbUsS5XnmdtJvzG3p/uAFJoLiQEkZAKP7pnZwJOFFkoQQlkinDR92vC1MEuv2FnPv+2u55/21ZO23lnsYcKq5lUughpxN2aZPIjFe+iJaQ8pNQlginXY82PF7Dz+6ad/BSq57eSm5JW4iHDZU5EYeAohKNjtUlx/RWEXrVHl8rN6ZBy5wNFqSQxyatCSEsEQ67HhpWbnpPz/uIq/UzSe/OYlnrh5HWYXVcnDFmosXSZIIKaVVXlxYM+kdER0bTJiRloQQlppyk997+GU5FmUVML5PEiMzEtBa86Wt3heQK0bKTSGmzO3FhZX87a6ODSbMSEtCCEuk0041DvRhhsCWub2s3VvMCQPMUg5KKRKdfvNkTZKQlkRIKavy4lLSkmgLSRJCWCKddsp0FLbq0kPutz2/DK1heHpc7bZ4V70k4YyWJBFiSt0enDXlJmlJtIokCSEskU4bJUTj8BwuSZgEUH/GbrzDZ+7YpdwUisrq90nYZWXX1pAkIYQl0mmnRMe0IEmUYVPQp1vdhKw4h5SbQlmZ24ujJknYJEm0hiQJISxOu40yonEdJklsyimlV3I0EY66tWKjbF6qcZilOKTcFHJMkqhp7Um5qTXalSSUUslKqXlKqa3WbVKA/V5USuUppda15XghjpYKWwwRvsBJ4p3l2czdkMvJg7s32O7CQzXWl4+Um0LO3oOVRNpqkoS0JFqjvS2Ju4EFWutBwALrcXPmAM0tZNPS44U4KiptsUR4y0DrJs+9uWQ3d72zmsFpsfxycv8GzzmpxlMzolzKTSFn5e6DpMfaQdll4cVWau88ifOBqdb9l4FvgP9pvJPWeqFSqm9bjxfiaCl3dsNe7YPiPaZs9MLp6LI8ynUkP5RfxgkDZvDcNeOJdjX8r+PUXtw1LQkpN4UUt9fHqj0H+W2GAzxSamqt9iaJNK31fgCt9X6lVGtXzmrx8UqpG4AbAHr37t3WeIU4pOzkCZDzHPz3akBD0XYUEEsZf4h8C/f5dzVJEAAu7TZ9EiDlphBz/wfrqfb66RnngINSamqtwyYJpdR8oEczT90b/HAC01rPBmYDjB8/vmktQIggsHUfyr+zZ3BD1XpsNht7x97Fk5tiiPUd5D73E1C0GLqf1eS4VPcu9qhUeoNJEr5qc61sqX/XKciCvA3mfsZYSMg8Ki/79eY8EqKcZCY4YZ98Hq112CShtT4t0HNKqVylVLrVCkgH8lr5+u09XoigGpgWx5+8l/P/8mz0S4lh8w+lRLvszL50ELzzBBRshSGNkoS7lB5V2/mYCzkeTJIAU3KKSjzK7yCEvXEZFGaZ+9EpcPMSiOl2RF8yr6SKvFI3D5w3HEfhxzL8tQ3a23H9ETDLuj8L+PAoHy9EUF1zQh8emjkSn9Zsyy9j8qAUPv7NSZw0cgBEJcH+1VBeYHauPAgl+2HvCmz4Wclgs91pzZ+QklOdqhKTII7/NVz8ElQUwOo3jvjLrttXDMDIjATweaVl1wbt7ZN4BHhLKXUdsBu4BEAp1RN4Xmt9tvX4DUwHdYpSKht4QGv9QqDjhegoTruNq47vw4XHZqKUmWBXK2UwrHsHNn0Ck26FRf8Anxtccfixsco/0OxXvyUhjNz15rb/VHPFvu/+BlnzYNItR/Rl1+8tQSkYlh4Py6olSbRBu5KE1roQmNbM9n3A2fUeX9Ga44XoaFEue9ONM5+FtW/D1w/Bwr9CdDfodRoUbWd+1Fkc3GG1ICRJNJVrTZHqMdLc9j4Blj4HOxdB3xMDH9dOBWVuEqKcxEY4wO+RclMbyFLhQrRUcj848TZTcuo5BibfUTvmfsXnm/Bv22H2k3JTU4XbwBkD8Rnm8eTfw/r34ct74IZvWjV3YWdBOev3lQAwKjOBXsmBr1ddXu0jpmY0ms8rs63bQJKEEK3hiIDLX2uy2W4DX80EPGlJNFW0HZL71yWD+HSYdj98fCtkL4Nexx32FN9tzef/zd/K8l0Harf17x7D/NtPxmZrPslUVHuJrmkV+qrBLl95rSVrNwkRBHabDZ9fo7WWJNGcou2mJQbklVZx97treHR9gnnu4K4WneKF73ewNbeUc45J57Xrj+cP5wxje345P2wrDHhMmdtHTISVGKTc1CaSJIQIArv1F7JfU1dukiRh+H1wYGdtknh7WTZvLt3D5zvMWko/rN5gkushVFR7+WFbIReP68VTV43lxIEp/GxiH5JjXLyyeGfg49xeYiJqWhIeKTe1gSQJIYLAbv1P8vnrtSSkT8LY+LH5Kz5tJFpr3l+5l3F9kljwhwvwYmfVpq3MXri92UPnrs9h8B8+Z8QDX1Lt9XPasLpFGSKddi47rhfzN+ay92Bls8eXV/vqZsj7PFJuagNJEkIEQU1N3B+o3LTwcXh8sPl557oOiLAD5awBwDf8Qh76dCNZeWVcMi4Tu92OPbY7A6Ir+cvnm3hp0Y4mhy7KKsCuFLecMpD7zh3OxP4NJ99ddbxZouf1n5ovWVVUe4mp6ZPwS0uiLSRJCBEENeUmn1+DI8psrJ8kdi82ZZfYVNj8WbOrzHZa2o+2uTjnyR94/vsd/HxSXy4eZ5bkUN0GcKpzHUOTFS82ShIen5+3l2czuEccvz9jCNed1K9JB3VmUjSnDk3jzSV7cHt9TV663O0jOqJeS0L6JFpNkoQQQWC3vrx8WoPNZvol6pebtN/U5I+5xGx3H/rCRp2K9uPV5mJNf714FA/OGIGjpj532oM4y/dzf9oi9hRVUub21h42f0MuFdU+BqXGBjixcc0JfSgsr+bztTlNnit312tJSLmpTSRJCBEENUnC7683DLZ+S8LvA2WDWGutzLLcoxxhx9E+Hx6tuHhcJpeO79XwyV4ToMcohpUvA+ChTzews6CcpTuLeGPpHrNt5shDnv+kgSn0S4lp0oHt82sqPfX6JKTc1CaSVoUIgpok4a1JEo2vKaH95oI3cVaSKM2BlEFHOcqOUVzpxqEVE/omN79Dn0kkLn+ZyQOSeWPJHt5Ysqf2qZnHZjS4TGxzbDbFzyb24U+fbGDd3mKzThNQ6THlp1gpN7WLJAkhgsCmmmlJNC43KVvdqrBVxUc3wA5UUFpJKjZG90psfofuQ1DeSv5zaS+2Vx/D/I25DEqNIz0xksGpcS16jYvHZvLYl5t49cddPHLRKMAMfwWIbjAEVpJEa0m5SYggaNAnAc2Xm2y2ur9k/V66ipKKavwoBgbqW+hmLYxYmEX/7rHcMGUApwxNZWiP+IAzqRtLiHZywZgMPli1l+IKD0Bt/0bdshyywF9bSJIQIghqk8Thyk22mvp410kSWvvR2Gp/R03USxLtcfUJfajy+Hl7uSlXVVSbclPtshx+r5Sb2kCShBBBUDvj2m9taFJusjqua0bX+JsO1+y0/D78HKJFEJduFv/L29Sy8/k85hddUQRz/wAf3waF2xjRM4FxfZJ49cdd+P2a8pqWRP0+CWlJtJr0SQgRBHUd11aWaLbcVL8l4TnKEXYg7cevDvH3qFLQ5wSzKuzJd5m5JPV53fDp72HrPDN02FMOMakQnQz5m8AeARs+ghNv45qJF3Pbf1fz/PfbySl2AzRa4E+SRGtJS0KIIGgw4xqk3FSPKTcdpm/h9D+ZBPDu9eYqdlqbL/53fwmPDYSV/4GyHJMg+p1shs6W5sCZD8OvfzBJef4DnJWUTVp8BA9/tql2cl5MhMO0PLRPhsC2gbQkhAiCuhnX1gZXbPOjm7pgkkD78R/u79G04XDWI/DJ7fDcqZA61Kz5FJUEA0+DUZfBwGlmiY+eY5tef+JXC+GxAbi2z+eT39xJTnEVf3jpY6ZUfU3/9x8Dn7W2k02+8lpLfmNCBEGDBf4AXFZLQmvzhab91uimehfA6Sr8LWhJAIy/1rTAvn/CXLFu9JUw458NS0QZ45o/NjrZXO1u4WN0X/UG3dF86NsLTsDdF3ocA+mjYdiMYLyjLkWShBBBYLeZLNGg3IQGT6VJGDUzrrtoS0Ifqk+ivtGXm5+2OO1BWD7HeqAgMh6GnQe9jjf9QaJNJEkIEQRNWxLWnABPhUkS2tdl+yTQvsOXm4Kh1wTzI4JKOq6FCIKaGdfe+uUmqOu81v5Go5u6UpIw8yREeJJPToggsDc3ugnqkkSXLjfpph3NImxIkhAiCBpcTwIalpvA+qK0m85rZetiScInLYkwJp+cEEHQdKnwmpZEmbmtmXENpjXRpZJEKzquRciRT06IIGiywF9tuclqSdQs8AcmSfi62Ixr+aoJW/LJCREEtsbXk2hSbrJmXINZZK4rrd2k/dInEcYkSQgRBPYm15M4VLnJ3qXKTUprKTeFMfnkhAiCZpcKh0blppqWRFfrk/AhXzXhSz45IYKgyRDY2nJTzTwJXa/c5OhSq8Aq6bgOa/LJCREEdS0Ja4PDZZJB7WS6euUmu6OL9UlIuSmcyScnRBDUzLiuHd0E5kI61fU6rm1ddwgskiTClnxyQgRBXUvCX7fRFVNXbvI3mifRlYbAIkkinMknJ0QQNLmeBNQtFw51C/xBl2tJ2OqX2kTYadcnp5RKVkrNU0pttW6TAuz3olIqTym1rtH2B5VSe5VSq6yfs9sTjxAdxW5vNAQWrKvT1S831U8S0ichwkN7P7m7gQVa60HAAutxc+YA0wM894TWeoz181k74xGiQ9ib65OouTqd1g3r8l2sJaGkTyKstfeTOx942br/MnBBcztprRcCRe18LSFCVk2ftK9+S8IVbSbTaasGVVNuckY1vLRpp+dHSZIIW+395NK01vsBrNvUNpzjFqXUGqsk1Wy5CkApdYNSaplSall+fn5b4xXiiGiyCiyYjuvqirokUZNJYtOgNOcoR9hxbDJPIqwd9pNTSs1XSq1r5uf8ILz+08AAYAywH/hboB211rO11uO11uO7d+8ehJcWIniazLgGMwTWU1HX/1DzRRmXbpJE/dJUJ6aoN5FQhJ3DXr5Ua31aoOeUUrlKqXSt9X6lVDqQ15oX11rn1jvXc8AnrTleiFDRZMY1BC43xfUwQ2PdpeY6zJ2c0n6UTVoS4aq9n9xHwCzr/izgw9YcbCWWGjOBdYH2FSKUNduSqC03NdOSgC5RctJaWy0JSRLhqr2f3CPA6UqprcDp1mOUUj2VUrUjlZRSbwCLgSFKqWyl1HXWU39VSq1VSq0BTgFub2c8QnSIgDOufe66iXO2ei0JgNL9RzHCjuHza5RMpgtrhy03HYrWuhCY1sz2fcDZ9R5fEeD4q9vz+kKEiiZXpoO65cLdJea2ttzUdVoSXr/GLqObwpp8ckIEQc3oJm/jchOA27qmRG25Kc3cdoGWhMfnx4aua0WJsCNJQoggsNkUSjWecV2TJEprdjK3EXHgiusSLQlTbtLSkghj7So3CSHq2JVqNOO6ptxkJYn6w0DjenSRloTGRr0VcEXYkU9OiCCx2VSjBf6slkRJtrmt/9d0XI8u0ZLw+q1yk7QkwpZ8ckIEiV2phvMkIhPM7SfWoD1HZN1zceldoiXh9WnsSuZJhDMpNwkRJHabajhPIv1YuPA5M7rJHgHDzq17Lt5KEp4qcEY2PVkn4fVrXGiUdFyHLUkSQgRJkyRhs8GoS5vfud8UWPQP2PYVDO28K+R7fX4ipdwU1uSTEyJImiSJQ+k7GWxO2PPTkQ2qg9V0XNuk3BS25JMTIkhsjUc3HYojAlKHwf7VRzaoDubzmyShZIG/sCVJQoggsdsazZM4nJ5jTJLoxKvBempGN9klSYQrSRJCBIldtaLcBNDreKgs4vvvv2L9vuLWHRsmvD6NDY1N+iTClnRcCxEkdnvLk0SVx8f96zJ4GBu+uQ9yrucuHHY7WkNSjItHLzqGU4emHeGIjzyvz2/WbpLRTWFL0rsQQdJkxvUhzPlhJ2+tr+BL52mcbF/Dhd2y8fg0Xr8mv9TNdS8vY93e4iMc8ZFnFvjzoaTcFLYkSQgRJLYWjm7y+Pw89XUW04amcs4dL4Ejir9lfsfiu09h4Z2n8M0dU9Eavlwf/jOyvf6aloSzo0MRbSRJQoggaTLjuhm5JVXc895aSqu8XDwuEyJiYfCZsPlT0vfNo3e3aPqmxBAb4aCi2neUIj9yPD5rqXC7VLbDlSQJIYLEabdR7fUHfL6kysOsF5fw9vJsxvZOZPJg61rt51iXdl/7Tu2+kU4blZ7wTxJerx+n8mGTJBG25JMTIkjioxyUVHqbfW7j/hJufn0FuworeGHWeKYNq9cpHZMCx/0SVr4K1eXgiiHSaaeqEySJauuqfJIkwpe0JIQIksQoF8WVnibbX1q0g7P+8R0llV5eu/74hgmixvDzwVsJW+cBdJoksa/QXHApLjqqgyMRbSVJQoggSYhycrCyusn2uetzUQo+u+0kJvbv1vzBfSZBdAps+ACAKKedyiD3SeSWVAX9nIeTlXMQAKdTOq7DlbQBhQiSxGgnBysatiS01qzbV8yVE3qTGneI1V5tdhh2Hqx5C6oriHTaqPIE7t9orb98tpFnF27HblP8YlJfnA4bSdFOPD7T0X7iwBRGZSRgs67VHQxfrNvP/PX7IJKGF1wSYUWShBBBkhDtxO31s7OgnL4p5oJD2QcqKa3yMrxn/OFPMPx8WP4SbFtApDOV0qrm+zdaavWeg7y5dA8Hyqv5Yn0OEQ4bbq+f57/fgcOmGlyP+7EvNzO0Rxzv3TSJaFf7vxayD1Twu7dWMz4zFgoAm3zVhCv55IQIkoxEU3ef+vg3nDWyBx6fZldhOQAjeiYc/gR9T4KoZNjwIZHOG8kvdbc5Fq/Pz7VzlnKgopqMpCguG9+LP10wku+z8hnaI56eiVEUV3hw2BUV1T7+u3Q3j8/dwpNfZXHX9KFtft0aT36VRZXHx18uOAaex7SURFiSJCFEkMwY3ZPeydE89fU2Pl9nJsJNGdydEwemMLIlLQm7E4aeA+s/IK7P9e3quF6dfZDC8mqevPJYzh3Vs3Z7/aU+EqJNP0FMhINbTh3EzsIK/v3NNgrLqvnfs4eSGO1q9etWe/3c+sZKvlifw9Qh3cmIt/oipCURtuSTEyJIlFIc2zuJv182mpe+38mZI9MY2qMFyaG+4RfAyv8wyr2SHzyD2xTH9vwy/rkgC5uCkwamtPi4P84YQaTTxptL9hDlsvPgjBGtfu1P1uzji/U5nDiwG49fMhq8eeYJSRJhSz45IYIsPtLJbacNatvB/aZAZAJjyr6h0jOgxYf5/Zp73l/LltxSVuw+CMC5o9Jb1RqIiXDw5wuOIae4ivkbc3ngvOEo1bKO7O35Zbz6425eXLSDQamxvHrd8ebYIqtfRZJE2JIhsEKEEocLhpzD0OLv8Xla3iexcGs+by7dQ6XHzw1T+vP1HVP51xXHtimEM4b3IPtAJbe9uYr9xZWH3X/vwUqufO4n5vywg/F9krj3nGF1ycVvlcwkSYQt+eSECDXDzydy9euM96/B7z/3sMNSiys8/PHjDaTEuvjg5klEONrXSTxjTE+W7SrirWXZfLxmH7+aMoC7zwrcmf2P+VsorvTw6a2TGZbeqLzmr2lJSMd1uJKWhBChZsAp+JSD42ybcB9iLSi/X/PU11lc9MwP7Cgo5+6zhrU7QYCZ7f3Xi0fzwc0nMiojgWe+3calzy5utiPd6/Mzb0MuZ45Ia5ogQJJEJyBJQohQ44jA44glnopDLvL33sq9PPblZnx+zYPnDTerygbRmF6J/PdXJ3DjyQNYsqOIUQ/O5atNuWhrpVuvz89jczdzoMLD6cN7NH8Sv/RJhDv55IQIQR5nHHFVFQGHwS7dWcS/vtpKRmIUX/3+5BZ3MLdWpNPO3WcN5djeiTzw4XqunbOMuAgHfq0pt5b4OHlwd6YNS23+BFr6JMKdfHJChCCfMy5gS8Lt9fHbN1fh82v+fMHII5Yg6jtzRA+O7ZXIl+tz2JZfjt2miHTaGJ6ewDmj0gMfWNtxLeWmcCVJQogQ5IuIJ04VNWlJ/LS9kMtm/wjAK9dOYErNNSmOgtT4SK4+oW/rDpJyU9iTPgkhQpDfFU88DctNWmse/nwTLruNh2aOZPKglk+U6zCSJMJeu5KEUipZKTVPKbXVuk1qZp9eSqmvlVIblVLrlVK3teZ4IbqkyATiVEWDlWAXZRWyes9BHpwxgquO73Nky0yLn4IXp8Occ2H/6rafR5JE2GtvS+JuYIHWehCwwHrcmBf4vdZ6GDARuFkpNbwVxwvR9USalkT96z88/W0WqXERXDQu48i+ts8L3z4KxXthzxJY/nLbz1XTJyFLhYet9iaJ84Gaf0EvAxc03kFrvV9rvcK6XwpsBDJaerwQXZEtMoFYVUWV28y69vk1i7cVMnNsRlDmQhxS9hKoKoYz/gQDToWseaD14Y9rjsyTCHvtTRJpWuv9YJIBEGAcnKGU6gscC/zU2uOVUjcopZYppZbl5+e3M2whQltEnKm8lpYcAKCovBq/rluO/IjxeeH1y81f/gNOgUGnwcHdULitbeeTclPYO+wnp5SaDzQ3U+be1ryQUioWeBf4rda6pDXHAmitZwOzAcaPH9/GP2uECA/RVpIoPlAIQEGZaVGkxEYcuRct2QdLXwB3MYy9BiITYMA081zWfEgZ2PpzSpIIe4f95LTWpwV6TimVq5RK11rvV0qlA3kB9nNiEsRrWuv36j3VouOF6GpsUYkAlBcXAFBYZq6d3S2m9dd4aLFPboctX0BibzjzL2Zbcj9I6gu7FsHEGxvuv/kLqDoIoy6DQJ3oMk8i7LU3vX8EzAIesW4/bLyDMkMwXgA2aq3/3trjheiSIsw6SCUHiwDYZ63GmhIXoCVRmgPfPwE56yDOavj3PRGOvQbsLfxvnrMWeh0P13wIznplreQBULyn4b5z74Mf/mnur3sPLnrOtDwaqzTlMiITWxaDCDntTRKPAG8ppa4DdgOXACilegLPa63PBk4ErgbWKqVWWcfdo7X+LNDxQnR51hfuvtxcdhaU89bSPfROjqZPcnTdPuWF8MXdUJhlfjwVkD4Gdv0AvmpY9w5smQuXvw62AN2PPz0LK181ZaGSvXDc9Q0TBEBChkkgNbLmmwQx+gqIS4fv/w4vz4CZz0DqsIbHluWBskFMGMzpEM1qV5LQWhcC05rZvg8427r/PdBsWzTQ8UJ0eZGmJRFPBauzD7Jyz0FumjoAh936sq88ALOnQlmuuTb2gFPh5LsgzbqanNZmrsPce+Gx/hDbw5SLxv287jWKtsPnd0HKYOg+BHocAyMvbBpLfCaU54HXDY4IWPxvcMXCGQ9BTDeI7mZe56Nb4eefmmti1CjPM89LuSlsSW+SEKHIKs/EqQo+Xr0Pn18zoV9y3fM7FkLxbrjybRh8RtPjlYITbgZvlWkF5G+Cj2+D5XMgprtJJjV/5V/9gWktBFLzXMk+00dRuBUGTzcJAmDSLVBZBN/9Dd64HC54GuKsa2mX5UFsWvPnFWFBkoQQoSgiDoAkWyWvbs7HYVOM61NvQYKdi8ARBf1PDnwOpWDKHeZ+dQUsftIkl/xNsHWu2T7+ukMnCIB46/nibNi30gyJHXNVw31Ovc+Unj67E16/FK6bZ1oUeRuge+ALFonQJ0lCiFBkd4IzhnHdbfRzxzBpQDeiXfX+u27/GvpMMuWflnBFm3LUyXfB7h/hxTPN9kGnH/7YbtbQ1/kPQGmuuT/myob7KAUTfgmxqfDWNfDCaeCpggM74YRbWhajCEmSJIQIVZHxnJTpZN75jVoLxXuhYIuZy9AWmcfBqX8An6duHsShJGTCsPOgIMu0cC54xgyTbc7w82H6o7DqNdMXcdqDTVsdIqxIkhAiVEXEQ1W9ead+n+mM3rHQPO5/StvOa7PDlDtbvr9ScNmrLd9/4o1N51SIsCVJQohQFZkA2Uvhw5vhxNth73KYd5957phL6kYyCXEESZIQIlQNOgNWvAKr3zQjlA7sgh6j4IZvA897ECLI5F+aEKHq5Dvh9rUwYqa5pkPVQZh0qyQIcVRJS0KIUHfeP0yH8P7V5laIo0iShBChzhVjRhcNO6+jIxFdkLRbhRBCBCRJQgghRECSJIQQQgQkSUIIIURAkiSEEEIEJElCCCFEQJIkhBBCBCRJQgghREBKa93RMbSaUqoU2NzRcbRAClDQ0UG0gMQZPOEQI0icwRYucQ7RWse15oBwnXG9WWs9vqODOByl1DKJM3jCIc5wiBEkzmALpzhbe4yUm4QQQgQkSUIIIURA4ZokZnd0AC0kcQZXOMQZDjGCxBlsnTbOsOy4FkIIcXSEa0tCCCHEUSBJQgghREBhlSSUUrcppdYppdYrpX7b0fHUp5R6USmVp5RaV29bslJqnlJqq3WbFIIxXmL9Pv1KqZAYwhcgzseUUpuUUmuUUu8rpRI7MMSamJqL809WjKuUUnOVUj07MkYrpiZx1nvuDqWUVkqldERsjWJp7vf5oFJqr/X7XKWUOrsjY7Riavb3qZT6jVJqs/X/6a8dFZ8VS3O/y//W+z3uVEqtasm5wiZJKKVGAr8EJgCjgXOVUoM6NqoG5gDTG227G1igtR4ELLAed6Q5NI1xHXAhsPCoRxPYHJrGOQ8YqbUeBWwB/vdoB9WMOTSN8zGt9Sit9RjgE+D+ox1UM+bQNE6UUr2A04HdRzugAObQTJzAE1rrMdbPZ0c5pubMoVGcSqlTgPOBUVrrEcDjHRBXfXNoFKPW+rKa3yPwLvBeS04UNkkCGAb8qLWu0Fp7gW+BmR0cUy2t9UKgqNHm84GXrfsvAxcczZgaay5GrfVGrXVIzV4PEOdc63MH+BHIPOqBNRIgzpJ6D2OADh8ZEuDfJsATwF2EQIxwyDhDSoA4fw08orV2W/vkHfXA6jnU71IppYBLgTdacq5wShLrgClKqW5KqWjgbKBXB8d0OGla6/0A1m1qB8fTWVwLfN7RQQSilHpIKbUHuIrQaEk0oZSaAezVWq/u6Fha4BarhPdiR5dsD2EwMFkp9ZNS6lul1HEdHdAhTAZytdZbW7Jz2CQJrfVG4FFM2eELYDXgPeRBotNRSt2L+dxf6+hYAtFa36u17oWJ8ZaOjqcx64+sewnRBNbI08AAYAywH/hbh0YTmANIAiYCdwJvWX+xh6IraGErAsIoSQBorV/QWo/VWk/BNKValAk7UK5SKh3Auu3QJmi4U0rNAs4FrtLhMcHndeCijg6iGQOAfsBqpdROTOluhVKqR4dG1Qytda7W2qe19gPPYfokQ1E28J42lgB+zKJ/IUUp5cD0Qf63pceEVZJQSqVat70xb7TF2bCDfATMsu7PAj7swFjCmlJqOvA/wAytdUVHxxNIo8EUM4BNHRVLIFrrtVrrVK11X611X8wX3FitdU4Hh9ZEzR9ZlpmYsnMo+gA4FUApNRhwEZqrwp4GbNJaZ7f4CK112PwA3wEbMKWmaR0dT6PY3sA0hz2Y/3TXAd0wo5q2WrfJIRjjTOu+G8gFvgzR32UWsAdYZf08E6Jxvov5IlsDfAxkhGKcjZ7fCaSEYpzAf4C11u/zIyA9RON0Aa9an/0K4NRQi9HaPge4sTXnkmU5hBBCBBRW5SYhhBBHlyQJIYQQAUmSEEIIEZAkCSGEEAFJkhBCCBGQJAkhhBABSZIQQggR0P8HSROfof/qUxIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(trial.tracking.time, trial.tracking.left_pupil_normalized_position)\n", "plt.plot(trial.tracking.time, trial.tracking.right_pupil_normalized_position)\n", "plt.xlim(xlim)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Obtaining pupil velocity\n", "\n", "1. Pupil positions are smoothed using the sliding-window averaging (almost identical to gaussian-kernel smoothing).\n", "2. The smoothed positions are differentiated, and scaled, to obtain the velocity." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def sliding_diff(x):\n", " x = np.array(x)\n", " ret = np.empty(x.size, dtype=np.float64)\n", " ret[1:-1] = (x[2:] - x[:-2])/2\n", " ret[0] = x[1] - x[0]\n", " ret[-1] = x[-1] - x[-2]\n", " return ret\n", "\n", "def smoothing_diff(x, rad=5, num=3):\n", " return sliding_diff(sliding.nanmean(x, rad, num))\n", "\n", "dt = np.diff(trial.tracking.time).mean()\n", "vleft = smoothing_diff(trial.tracking.left_pupil_normalized_position) / dt\n", "vright = smoothing_diff(trial.tracking.right_pupil_normalized_position) / dt" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABNkElEQVR4nO2dd5iU1fm/7zN9tldg6UWQJiAioNh7R4015hujRk2iv5gejakmMUZN02iMFTX23hDEglhRRJFepcOyhe3T5/z+OO/szFZmizs77HNf114z89ZnZmfez/uU8xyltUYQBEHou9hSbYAgCIKQWkQIBEEQ+jgiBIIgCH0cEQJBEIQ+jgiBIAhCH0eEQBAEoY/TLUKglHpQKbVHKbWijfXHKKWqlVJfWH+/7Y7zCoIgCF3H0U3HmQP8G3iknW3e01qf0U3nEwRBELqJbvEItNaLgMruOJYgCILQs3SXR5AMhymllgE7gZ9prVe2tpFS6irgKoDMzMxDxo4d24MmCoIgpDefffZZuda6uCP79JQQLAWGaa3rlFKnAS8Co1vbUGt9L3AvwLRp0/SSJUt6yERBEIT0Rym1paP79EjVkNa6RmtdZz2fCziVUkU9cW5BEAShfXpECJRSA5RSyno+3TpvRU+cWxAEQWifbgkNKaWeAI4BipRS24HfAU4ArfU9wHnA95VSYcAHXKSl7akgCEKvoFuEQGt98T7W/xtTXioIgiD0MmRksSAIQh9HhEAQBKGPI0IgCILQxxEhEASh64QDsORBCPlSbYnQCUQIBEHoOqtfgVd/DO/9PdWWCJ1AhEAQhK4TqDWPNTtSa4fQKUQIBEHoOoEa8xgJptYOoVOIEAiC0HV8e82jjBNNS0QIBEHoOr4q8xj2p9QMoXOIEAiC0HX8VeZRqobSEhECQRC6jt/KEYgQpCUiBIIgdJ1YkjgsQpCOiBAIgtB1YkIgHkFaIkIgCELXESFIa0QIBEHoOpGQeRQhSEtECARB6DqNOQIpH01HRAgEQeg6jR5BQ2rtEDqFCIEgCF0nJgTRcPy5kDaIEAiC0HUSewxJv6G0Q4RAEISuI0KQ1ogQCILQdaJhcHjN80g4tbYIHUaEQBCErhMJgisj/lxIK0QIBEHoGlqbi78z07wWIUg7RAgEQegaUSsUFPMIohIaSjdECARB6BqxclGXeATpigiBIAhdI3bhd8ZyBDKOIN0QIRAEoWu08AhECNINEQJBELpGc48gKkKQbogQCILQNWJCIDmCtKVbhEAp9aBSao9SakUb65VS6g6l1Aal1JdKqandcV5BEHoBjVVDMSGQqqF0o7s8gjnAKe2sPxUYbf1dBfynm84rCEKqEY8g7ekWIdBaLwIq29lkNvCINnwM5CmlSrrj3IIgpJgWVUMiBOlGT+UIBgHbEl5vt5a1QCl1lVJqiVJqSVlZWY8YJwhCF2heNSQDytKOnhIC1coy3dqGWut7tdbTtNbTiouLv2azBEHoMjEh2F88gi+egC0fptqKHsXRQ+fZDgxJeD0Y2NlD5xYE4eukRY4gjctHo1F48Xvm+e+rU2tLD9JTHsHLwLet6qGZQLXWelcPnVsQhK+T/anFRP2eVFuQErrFI1BKPQEcAxQppbYDvwOcAFrre4C5wGnABqABuKw7zisIQi+gxYCyNM4R+Gviz0N+cHpSZ0sP0i1CoLW+eB/rNXBNd5xLEIReRmwksSvLPKazRxBIEIJAbZ8RAhlZLAhC12gMDe0HyWJ/Ql4gWJc6O3oYEQJBELpGY2hoP5iqMtEjCNanzo4eRoRAEISuERMCuxtsjjT3CEQIBEEQOo4VGlqyvQ5td6V399EmHoGEhgRBEJLDEoLvPPIFYezpPY4gUBt/Lh6BIAhCklihoBAOQjjSWwjC/vhzEQJBEIQksS78IRyEtT29cwThBNslNCQIgpAkkSBRbESxEcSR3gPKwv74eAjxCARBEJIkGjIhITA5gnAgxQZ1gXAAPHnWc3+7m+5PiBAIgtA1IiHCygiB8QjSOEcQCVATcRC1uSDkS7U1PYYIgSAIXSMSNJ4AENTpXTUUDfnZURulLpLmnk0HESEQBKFrRIIEtfEI/Dq9B5QF/Q0EcBDAhRaPQBAEIUki4cYcQSCa3h5BOBgggIsATkJBEQJBEITkSPQIor2gfLShsmnzuA4QDfkIaCd+7SIcaOhmw3ovIgSCIHQJHQkS0HYcNkVAO4iGUygEWsOtI+Dhszq3ezhAEAd+XEQC4hEIgiAkRTRsksUFmS5CONCpFIKGCvO464tO7a7DAQI4CeAkKjkCQRCE5IiGAoRwxIUglaGhqq3x54GOjwxWYT9BTGhIksWCIAhJEo2ECOIgP8NFMNUegW9v/HldaYd3t0WDBLTxCHRIBpQJgiAkRTQcJKQtjyDV5aOJQpDYSTRJ7JEgAZz4caFkZLEgCEJyaCtHkJ/pJESKq4b8VfHnnRECHUQ53EYIIjKgTBAEITkiQYI4KMgwOQKVyhYTXfQIHNEgODwEcWIXIRAEQUgOHTFN5/IzTY5ApXJAma8q/ryjQhCNYCcCDhch5cYeFSEQBKGjhAPw+WMQjaTakp4lJgS9wSPwV5l5k6HptJPJYPUWitrdRGwuHCIEgiB0mEW3wUs/gNWvpNqSHkVFYzkCkyy26TBEo6kxJtgA2SXmeUc9glhy2O4mYvfg1EEzQK0PIEIgCN1F9Q7zWLs7tXb0MCoSIqjjOQIgda2oQw3gzQdl74QQGA9AOdxE7W5rWd+oHBIhEIRuw7p7TExY9gFUNGzlCJwEYkKQqsqhYD24MsGd3XEhiCWHHR6ido95LkIgCEKHiCUq+5wQmBxBXqJHkKqEccgHzgwjBqEONo1L8Ai0w/II+sigMhECQeguwlZLgsRa9j6AXRshyHDawe40C1PlEYQawOk1QtDROYetu3/l8qAbPYK+0WaiW4RAKXWKUmqtUmqDUur6VtYfo5SqVkp9Yf39tjvOKwi9ithdcB/qUYPW2KNBtM2BzabA7jLLUzW7Vyw01CkhMOJlc3jAGROCvlE55OjqAZRSduAu4ERgO/CpUuplrfWqZpu+p7U+o6vnE4ReS+yi0ZeEIBpBoYlaAmBzuCFEikNDXnBlddojsLvc4PDGj9cH6A6PYDqwQWu9SWsdBJ4EZnfDcQUhvYiFQ/pIghFoTLBqmxEC1StCQ5kmTxDsWPfRsJUPcLi8KGesaqhveATdIQSDgG0Jr7dby5pzmFJqmVLqdaXUhG44ryD0LhpDQ31nZqvGC6VVbtl4AU2FEGjdNEfQwf9D0G+2tzs92JyWR9BHcgRdDg0BqpVlzUdhLAWGaa3rlFKnAS8Co1s9mFJXAVcBDB06tBvME4QeIlZ+2EcqTYDGC752GI/Abj2mJDQUDoCOgiujUzmCkDUjmcPlweaKhYb6xv+yOzyC7cCQhNeDgZ2JG2ita7TWddbzuYBTKVXU2sG01vdqradpracVFxd3g3mC0EP0YY9ANeYIYkKQAo8g9rk7MzuVI4gJgdOdEReCPhLm6w4h+BQYrZQaoZRyARcBLyduoJQaoJRS1vPp1nkruuHcgtB76JM5AvOeYyEheypDQ41C4DVeQbCuQy0iwjEh8GTgsIQgHJTQUFJorcNKqWuB+YAdeFBrvVIp9T1r/T3AecD3lVJhwAdcpHUfaeIh9B0aq4b6kEcQiZVcNheCFISGgtbnHisf1VHzP4mVgu6DiLW/05OBw20EIBxo6Jb4eW+nW96jFe6Z22zZPQnP/w38uzvOJQi9lsbQUB/yCBJG4wI4UuoRWKGgWPkomPBQkkLQxCNwm/9hKNBAcnunNzKyWBC6i9jFLxLoO62orfcc8wRSKwRWGCdWNQQdKiHVIR9RrfC4vTjcJjQUCfSN0JAIgSB0B9Go6bjpzDCv+0qewPII7NZdt8PdC0JDsXEE0KGEcTToI4ATj8uBKyYEfSRHIEIgCN1BrO1yLCSRynl7e5JYjsBlhMBpPUZSMRArlptxZcT/Dx3I10RDfgI48TrteF0O/NopQiAIQgeIXfhjIYlUTtfYg0StfEjMI3BZQhAK9qxHFI1qwgErDBTrPgodG10c9uHHRYbLjsdpJ4Cz8f3t74gQCEJ3EI4JQd/yCIIBc6F0WyEhpzsmBD3rEVx078fc/cZy86KJEHRgLEHIj1+78LjseJw2/LjQ0mtIEISk6aMeQdC684/F1F0uIwjhYM8J4bbKBj7ZXEl1TbVZ4OqkEIT9+HHhdVoegXaKEAiC0AH6qBDERuO6LE/AbT2GQz3nESzZUglABtY5O+kRqIifIE6cdhsZLoflEUhoSBCEZIld+BuFoG+EhkJWaMjjMR6Bx+0iqlWPJlm3VDSgFAzLUWaGNLuzU0JgCwcI2YxHk+my48fVZ6q/RAgEoTuINZzrYzmCWFLY5THlmh7rTjrag0KwtbKBATkeRuXZ8GkX1Q0hU0IKHROCSICQMr2SvC6TLJY21IIgJE8fDQ2FLSHwWh6B12nHj5NosOfabGzf62NIfgaDsjQNuPloUwXYHaY1dgeqhuxRPxHLI8hwOfBrFzbxCARBaEI0CuveaP0uMXbhd/ctjyBiVQd5vCY34HXZ8eHu0SRrRV2AomwXhc4wAdx8uLHcrOjgnAT2aJCwJQR2myKkXKiICIEgCImsfQ0ePx8W/qXlupg4xDyCaN/wCCIhPyFtJ9NtQioZTnMn3ZNTPFbWBynIdGGP+MGVwQcbYkLQsVbUzmiAiD3eWShsc2OPSGhIEIREdi0zj7uXt1wXaT6OoI8IQThIEAcZLjsAHpfNSrL2jBBEopoqX4iCTDcE63FnZLOxrJ7d1f54K+okcUYDRK2Z1sASgqgIgSAIidTsMo/15S3X9dGqIR0yJZcZbtPI2OQIXNh6SAiqGoJoDQUZTgj5yMrKBjBeQQdnKXPqINoRF4Ko3Y1DhEAQhCbUWhPv1ZW2XBdpFhrqI0IQDQcI4sDrtDwCpx2fdqF6KMlaWW8+54IsN4QayMzKIdvjYOnWvZYQJJkj0Bo3AaJ2b+OiiN0jQiAIQjNiHkHdnpZtpiPNm871jdCQCvvw48ZuM1OXO+02gsrdY9U2FZYQFGa6IFiPcmYwbkAOq3fVdCxHEPZjJ0o09v8DtN2NS/cNQRchEIRkiXkEOgINlU3X9dHyUVvIR1C5mywL2twmcdsDxDyC/AyXyQe4sxhXks2a3bVoZwdyBFbDumhs/AGAw4ODcJ+YW0KEQBCSIdgA/mooHmteB2qarm9eNdRHQkO2sI+AzdtkWVi5cUR7VggKs1zgrwF3DuNKcmgIRqiLupP2CHSg1jxJ9AgcVgVRHxhLIEIgCMlQa4WFisaYx9iFI0YfDQ3ZIr4mJZcAYbsHZw8LQb5LmzyNJ4exJTkAVIScSY8jCDQYYbd7shuXqdgUl31gdLEIgSAkQ40VFmpTCLqnDfVDH3zFpQ9+YtokpAGOsK9JghV6NslaWR8k2+3AFbZCQO5cxvTPQiko9TuMRxCN7vM4vroqAOzenMZlNqf1vvpAB1IRAkFIhphHUHygeWwee+6GqqFgOMofXlnFu+vKePijzZ2zs4dxRP1EnU2FQDu8uHQAtP7az19ZH6QgyxUP1XlyyHA5GFqQwQ6fHdAQrG33GAD+OtPC2p0R9whsbvO/jAY6MLlNmiJCIAhtcNc7G/jza6uIRnWCRzDaPLYZGoqNLA53+HxrdsfzDm+t2dPh/VOBSwdQsfmBLRxuLzZ0j+RJYqOKG4XAbe7oD+yfzfo6K4ndPLHfCrHQkCszt3GZPbMAAF9NK+NG9jNECAShFT7cWM5t89dy33tfsWB1KdTsAHcuZJeYDZoniyNBsDnAZjePnbgIrtxpjnnOwYNYvr2Kal/vDg9FoxqP9qNczYTA0/H5gjtLRX2QggwrUQzgMUIwdkA262tN24tkhCDUYDwCT4IQOLKKAPBVp4codwURAkFohVeW7cTjtJHtdjBvxW6o3g65g8FthQ6ahwvCAbBbFx67q1NCsGOvD7tNcf4hg4lqWLypoovv4uulrC6AlyAub1aT5W6v8Yp6Yk6CyvoARVlu8FkXe08eAGNLciiPWvH+hn1/jtE6c9fvyS1uXObKMUIQ+Bo9gmj06w+fJYMIgdAlAuEIugdiwT1JKBJl3ordnDxhAMeM7cdHGyugehvkDjKzX6FaDw3Znea53dmpqqGdVT4G5Hg4ZHg+HqeNJeu2wqcPGBHqhWyvqMWtQmRkZjdZ7rHi7NVVe7/W82utqaizcgR1ZWZhVj8Apg3LpxLLriSEgIYK6rWbzIT34rZEIVT79QjBIx9tZsyvX+ftNa2MVO8oDZWw+f1O52VECIROs2pnDQfftICrHv2sVTGIRjXXPr6UGTe/yYY9PZhw27oYVjyfVLVIa3ywoZy9DSFOO6iEyYNz2V3jJ1pleQRKgdPborbctFpwmnJGm7NzHkGVj4F5HtwOO4cOL+DoFb+C134Cc86A8NcQb2/v8/HXwAd3wN4tbW5SWmYuvlm5BU2W5+SbO+ndpbu7bmM71PjChKPajCqu3wPKBhmFAPTL8ZBbOMBsWL/v0I6qL6NS55jxCBbZ2XkEtKPRW2iLbZUNvL2mlLnLd/Hal7t4ffkuFq7dQ0Vd25VTWmvuWbiRcFRz/3tfJfFu20FrePxCmHM6LHmwU4dwdM0CoS/z0hc7aAhGWLCqlE++qmTGyMIm65du3curX5pqm7/OW8N9355GtS/EHW+tZ1RxFt+cMTSp82iteeaz7QTDUS6ePrSxnUGr7FkDD55knp+wBY74cYff15OfbKMg08UxBxbz6Vd78eLH5t9rhADA4WkhBJ9tKmWgT3PuvxbxocuJvRNCUFEfZHQ/E2b57sDNzNr2KTXFh5BT9hmsfAEmX9jhY7bKri/hhavNRf68B+DAU1tuM+96+OIx+Px/8IOPqQtFWbh2D4cMy6ck11QJbdtpEui5Bf2a7DpwwEAAdu7exfgkzLnrnQ088tFmvjVjGNcedwBKtfP/TaC83lxoi7LcsL3MiIDN3rh+5rjhVH2SiafsKzxtHcTC7q9kr8phiCt+SSzO9rBdF+Os2mSJpm5y/Iq6AD99ZhkL15a1ekyXw8Zt501i9pRBLdbtqPKxs9qP12lnyZa9BMNRXI5O3peXr4Ptn5jnH9/dqUOIEAidZmNZHYPzvZTVBpi7fFcLIXhrzR4cNsWFhw7hyU+3safGz9/eWMdTS7YB5ody3iGDWz94OADPfReqtzN35G/4xZsm3rxoXRl3XTIVp72NH81bN4ErG4oOMHe0M74Pzn1dBgzbKht4+MPNzFu5mx8cMwq3w87wogyGKeuOMm+YeXR6IWFS860VDeysqGao202dP0xZWJNR7+PXT3xOaY2fX5wylkOG5e/z/DW+EDke00XzyPW3sp3+XF7/S57J+BmexfcTGvsNdlf7GFmUhc2mCEWiPPzhZrSGSw8f3vaFxF9t6ulzBpo+SY9fCDoK3nx4/Rcw+mSwJewb8hmPyp0L5WuJbHyHS+a7Wba9mhyPg/svPRSbgkXL1nE14Mpq6hEUFfcHYOnarzgyHMHtsNMWa3bXcNv8tRRlufnbgnWs21PHFUeMYPLg3H0KQkVdwqjiuj2QWdxk/blTB7N1cT9yt6xlGLB9bwML15Zx9JhihhQ0TXB7AuWU2vOaLCvMcrOEoZy25224ZYgpAph1Hcz8ARurwlz20KeU1vj5+ckHMnNkIZluOwpFVGtqfCFum7+WG55fzpGji01lU4xwgO1LXqOQEBfMmsp/Fm5k9a4aJg9pev6k2f6peZx5DXx8V6cOIaGhPsiGPbXMW7GLQLhrPVQ2ltVz0KBcjh5TzLyVu1skvt5evYfpIwq48siRRLXmJ08v46kl2/juESOYPqKAP7y8kp1VbSQUP7wTVr8Mu75g0ofXMn2wh9+eMZ43VpVy3ZOfE2ktybZrmZk8ZtYP4YTfg6+S91++n+89+hlbK9qvYHlmyTaO/9u7PPThZk4/qIRrjj0AgH7ZHg5UW81G/SeYR4enSb/9hz78CpeKUJibze3nT6Y+rHh39Q7eWl3KV+X1XD7nU/bU7Hukba0/TI7XAQv/gq1yAzUn3EZQufh39RG4d37COb9/gBP+vojz//sRe2r8XP/ccv702mr+PHc11zy+lHCklVBPzU64Yyr8YwK8/0946lvg2wuXPA3H/waqtsLOpU33KV9v3t/pt4Mrm63vPcay7dX8+IQxFGW7ueC/H3H+fz9ioMt6T96mIqes17VV5fx9wbp23/ND72/G47Tx5k+O4qcnjuG1L3dy9l0f8P3/LW39f5xApeURFGa6oWIj5I9osn5cSQ7+rKGoyo38Z+FGTvvXe/z6xRWc+e/3zXwFYMIqgVoKgzuocDe9KbHbFMtdU8yL4UfAkOnw1h8I/HMqc+76M/hreOKqmVxz7AEcMiyfsQNyOHBANuNKcpgxspDfnTmBhmCEd9clhKZCfnjoNGZ+8F0WuH/BecPNe9hY1oXQ6fYlRrRnXQck5001p1uEQCl1ilJqrVJqg1Lq+lbWK6XUHdb6L5VSU5M5biAcZf7K3eyt34eb3Up8OukE5o7P4J4j4L7joXRV4+LSGj8rdlS3/2X0VcGSh2DzB8mdqxmRqObeRRv5+xtrqfV3vlRQa01dIJzUe/58615Ov+N9vve/pXzjPx/ywYZyHnz/K254/kuTFO1AsmlXtY/B+V5On1DMCfWvsXXBXY0NurZVNrC2tJbjxvZjeFEmJ43vz/sbyhmQ4+FHJ47htm8cxOF6Kbffcx8bSpuVYlZtg0W3w7gzKTv3aYZEd/A7z9NcfsQIfn36OOYu380db61vadCi280PYsbVMPwofFlD4IvHmbdyNzfPXd3m+/hsSyXXP7+cQ0fk894vjuWuS6aSafXXdzlsHOreSlg5odCIAw5Po0dQ4w/x9KfbGJJjx+F0cepBJfTPz2FcsYc3f3o0T1w1E38owu9eXtnuZxmKRPGFIowMb4QP/w0Hf4vxR8xm4c+P5ZtX30BEOfn7iCVcf+pYVu2s4bBb3ua5pdu57vjR/P7M8SxYVcqNL6wg1FwMPrkPGsqhZAq8+TvYthjOvhtKJsMBJ5ptvlrUdJ9y6+LdfyKR0SeRt3UBkwdm8f+OO4Bnv3c4/zdzGNceewB/OMkKeTQTArx5ABw20M4jH25pbAPRnIq6AC98sYNvHDyIPI+D/3f8aD698QR+euIY5q3czQPvb2r3MyuLeQReBZUb44P9Ehh7yNEMVaU8OO9jBuVncP+3p+ELRrjp1ZWw7VMjkH8ZjEf72ZU1scX+i/PP4gf9H4NvPgWXPMOKEx5lY52TP3I373IFUxdeBmvntWrf+IGmJfanmxOS5q//HHYs4dm8y7HbFCPe/xk2pdlc3k5PpJAPXvge/GUoLL635fodS2DQwZDdn+iQGe1+Zm3R5dCQUsoO3AWcCGwHPlVKvay1XpWw2anAaOtvBvAf67Fd1pXWcvWjn5HldvC3CyZz8oQBTTco3wAvft9czAdPg+N+TcOgw/ndSyt56Yud5GY4OefgQXz7sGEMzs8gHImytbKBTWX1jOqXxQhPPfzvPPPDjobhoVPg/17k0a0F3PTqKkIRzaA8L5cePowLDx1KrtcZP3dDJTx4cvxHc/Z/YMo3W30fFXUBdlb5Gd0/C4/Vtz1SV8EjTz7Gik27eS96EJ9u3sv/vjujSfw7ULaJ6Nr5eDa/hdq22NyRnvAHGBr/6L7cXsUvn1vO6l01lOR6uHj6UK44YkTjhSyRbZUNXPnIEvrnePje0aP467w1XHL/YgA8TkXx0n9xqOtVHAMnwQUPQ85Aqn0h3l9fzqkTB2BLsM0XjOAPRSnIdHNqxUPMdj4IHwHR7XDqLbxtDYg6bqyJH9/6jclMHbqVkycMIMtlJ2vRT/mv7Unww1N3vcPcI//C1ceOxm23mXAFwMl/4Z31UBc+hcu3Pwn3rOIKbx5qzGn86W04YnQRhw63whJ7VhsP4qifgycXrTUv6qO40P44V4518uS6cqJR3eQ9gBHjvz21gGnZYe751olke5xN1uOv5jj1GWtdk5gQqwpyxnMEzy7ZTn0wwrA8J2jj/mdleDkgyw1WLP2Hx4/mtvlrWbCqlBPH92/1O1LrNwPQjth2j7mwnvSnxnUjhg6FyRdy0LLHOWiAl1PPOoXb1xUzcXABVx01EqUU5XVB/v3OBtaU1vLPC6cwoijTJJg/fxQOPA0u/B98+bSpfBpxlDlwZiEUjDS/n0TK15nEa+EoXo9M5wye409Ta7HZFAWZLv54tnXB/MQSkOZCYHeCO4eZ/aP4tkV44fMdXHFE07t1gMc/3sJV+jl+vGYe/MkPRWMoHHkM1048l8+3FnPnWxs4/+AS8iPlkDO4afgK2FXlo8RWRfHnd5rfb8xjSyBn3HHw3k08fWw1A48/HLfDzg+OOYB731xGcPsNuDwZRI65kdvf2gzDWuZKDizJ5dUv64lGNRvL6rj4TTcDc//Js2fYyd7yJqx8EZ64EA671vzPEsJZdptiXEkO63ZbFWafzYGlj8CRP+POpUcQ6FfAJdtv5/yclXxVkZBH0BrWvGqqnYrHwhu/Nnf9+cNN7mbUcSb0CaYZYumqxlxY2eATgTdavI990R05gunABq31JgCl1JPAbCBRCGYDj2hzy/qxUipPKVWitd7V3oEHuAM8ceVMbpm7iiceewD3RDfHnHYh5A2B5c+iX/4hEZsLDr0Sx9rX4OEzWeKcxZ76I7hxdB5hXw0ffPAJZywaiSOrkBpfmKB1x6QUzC+5n9HBetRlr5u478NnEnzoLNb7zuX0kWdw5JTxPL1kGzfPXcM/31zPKRMGcPCwfAZ4wkx970py927m4SF/YVbZUxzwyk9g4DQc/cY02h+ORPnXW+u5e+FGIlFNtsfBhROzuMI3h/4bn+UyolzmgqAzh19vuZB/vZnPT04cg974Nttf+ytD9pqL9GZdwmrXdA7ZtoyiB09m/Yhv0//sP/LG2mruf2kBR3q/4obxDj6pLebuBdU8+vEWfnPGeM4cEjACUldKbV0dLy+tYHq4mJ9d+m1GDhnEqRMH8MnmSkb3y2LQrgW4n3+WReHJzCpdif2p/yNy2XyufGQJn3xVyd8vmMy5U+Ouc2WDuRsbV78Y16f/4OPc0/iqOsrFi/8DE7/BG6siHNAvi5HFJvmZm+Hk6qNHmZ0/mwNfPglH/ISGQIgLP72TFxb9hAu/vIbHxn1E5tq55keVN4SVO1fwku1bXHbUJNTOpaiKjVxedSMfZt/C9c99ydzrjjQx6Pf+Ds5MkxMAPtxYwV2V07nQ/QTnBl7gvsAZ7KkNMCC3ab7go1WbuKf+R2QrP2rXkPhFEoy4PHgyA6PVPOa6lMbLjMNUDWmteXrJNiYPziXbEYWoNZK12TiCq44ayctf7OS3L61gwsAcBuY1bckAUOsPMVjtYUj5+3D09S0vrqfeYr60y59l2NKHubPfeDjxocYLz89OPpBxJTnc8PyXnHXn+/zq9HFMqnyDCfVlrBx0HkOCUXKmXNzyR1YyOT4FZ4zydZA/nDfW7uWGL/txktvNQTWLgDOablez08TNvU1zBADkDqEguJvJg3N5Zsk2Lp81vEnMPxiOEvzwbn7mfAZGnG5GbO9eDp/eh/r4Lv6bPYTFOgfPP7dBpAYGT4dvv2Smn7Qo2vo6C103Y1sUhKIDYczJrby/KVA8jhGbnwK7+W5cedQICj+8CZdvD/qSt9jiHst/5r3L7SUt8ziHDs/niU+2cuv8tdbYEjsPXDaD7PwMGHMUHPcbmHcDfPRvU1J8xj9MQtmas2JUcSbzV5bC+jfhtZ/CASfgP+KXbHtzAeWTvgG+J7g08Bq/LJ8VP+nSh+GV6+KvnZlw4aMwZIbxYJY8SPXRf+CDDeUcZl9Lvo6YG2FgeeGpwM9bfg77oDuEYBCwLeH1dlre7be2zSCghRAopa4CrgI4pMTOYdvu44WMd7C5PoZ1oNf9CQpHoSo2sFQfyDWBa6n/pD9HDj+LcXvn8N3Qixzl+gCsqrcrnIATvnJN5sth5xI6cDZDi3NZ9fb/GLPtbR7L+g5TIwMp8Lh4YeQdHLLk59zkfBi2PQzhKXxj8KHs6e9m5bYKdq2pxrm8jgNtq8lX5VwXuobFu8fwItfwaPhHbL/nEpYc9wRHjOlPtK6CJc/9Daoa+NbEbzJtwjiWLl/BBcuvoohdPBo9kfzpF3HWwcNwLfgtt265jzXvz8O3DLz123DpPF4ouJzQ2LNYHezHrio/r4brOHH7vzn7q4ep//uTnIKNC5w+CAOb4CjgJ5lulusD8TxXibLF68+zgWtiLx65E065mfxDvmO8LH81PHw94X4T+WXljZxh/4gbd/yDNx7/O598ZS5/S7fubSIEe+uDFFDD4ct+Bf0mED32Nv708MfMzv0M12s/Y/HWn3PVUQe0/LZU74D5N8LIY+C435Bhs0F2Hue8/UfOqf0APgGmfMskvoDVu2sZVVKIOu5Gs7+vCnX3TP7ueozJO37Ov95cz88PdcKKZ1lS8k3WfFnL9BEu/jpvDeHsIUQnfJsDv/gfo9TB7K7xtxCCXcvf5Qhl5Q/e/jNckSAEc38Oysatg+7kjdph8Z+X0wMNleypDbBmdy2/Pn0crAsazxIsIYiH+px2G7edP4lv3reYY29fyJGji7j9/MnkZcQTiDW+MEfZrLmQDzq/5efmzobZ/4ZTb4U1r5k7w/uOhTP+2VhNdPpBAzh4gJPrnvycG59fxkuuO9lICWfMdeF6401+ecpYLm9+Z154AKx62XgPDmOPLlvLZgZx9f8+Y9KgfqiCE2D1K3DKLU3vyivWmwS6vZXLSP5w2PsV508bwq9fXMGKHTUcNDg+aveNz1ZxdeQJKgYeTeFFj8XvpH1VsOY1nGteY2BoM2/6J3PKUbNwLroF3vubyWsA1Oziwl23ssU5kjFXPmzCQq0ll5WCmd+HV34Iq16CCWeTUb2Rb+q5PBU+Btvu/lTUmzr+g4fmtdj9tINKeG7pdu55dyO5XiePXjGdwfkJiWa7E067zYxofu9vsPtLk4i3xPWa/JmMDuShH1+A6jcOznuQrVUBohqG98+FrKsZv+C3uP0r0XoWSmt4/x8waJoJ4e360uQncqwR7aOOI7jyZU5YcixldUF+kjmPH4LZHlhR1cyjTZLuEILWshPNA83JbGMWan0vcC/AtBH5moV/webNJ3rGHdy2OpeMtS8wq3ozc0OXsGHEJdwwbTiL1pXz8aYKoqO/R+UJv2dQcIu5w3dlGWXe8gEjlj3BiI2/g72PwIRzmF76X6pyx3Fb1UlU/eu9xvOfP/U+Jh+ucW16EzYuhGVP0C/ko5/diXa5iGZm4ssZza5D7uCP404gL8NFNKpZMb+aSYt/xiuv387Vr01ljvNWvmkrAyeweT5kncuZu+eh3fUsP2IOJ046MX5neOkrBD95AP+Ch3m/xsXrkdPJmHoBN507tUUoQ+ujWL9kAQ2fP0uOx07W+JnYhk6HrP6wcym29QuYtO1TdvkHc0vFCbwdHMcOXcTBI0v4/clDOCC6xXxhX7nOfD4HnQdv/RHq9+C4+HH+ERjO/z0Q5gzngRyy8S6+Pe1plu0OsaVZsnVvQ5AfO57FEaqF8x7gsOJBjBs2kFvKL+am3XdwjnqX0w46suU/eN71Jo9w5r/iF5WjfgaDpvLmgtd4eEcJ/z31OjJsNrTWrNlVw5mTB8b39+bB0b8g99Uf86uxpdy8UDFp8X84OurgB5sOp2xTPBZ/9yVTcYyYRHjVi9wUnsOuqjOZ0qwyw7ZrKVEUtqN/Ae/+FcrWmotKxUbY/B6c8HvKd0+mriyhltwqH11fahJ840tyYHWwsc8NdmeLPviTBufxyv87gvvf28Rji7fy9JJtXHXUqMb1tf4Qh9tWEMgYgLtwFG3iyoBJ55uLw3NXwAtXwaZ3zPd91csMbCjnGSBWL7nzuH/xaMlhzPnwK2561TjpTcSgcLSZaGfvZigeQ3lNA7ll63kjNIJzpgzi5nMPwrlqB6x7zYSQhhxq9ouEYOvHRtBbI384bHybMycN4I+vruKZz7Y1CoHWmop3/0uW8pM5+89NL+DePDj4Ejj4EnZvKOf/3b+Yf+RP5pzxa2DxPTDzByak9fovsOswzw77Lb/qN7btzwvg4G/BJ/eakOOQ6fDKdShXBm8UXc07z32JUopZBxQyqjirxa4ep53/XTGDr8rrKcp2m6qu5igFx/8W8oaaAYDubDju1xCNUPTR/VzuKKN65Jnknn83eHLYtdWUmw7K88KBlxJ662bOD86lvO47FFd8Cns3Ez36enThGOzN8h4NI08hY908xjq38ONzTiLv1X9SnTGQ3CxTMbWudN8N9lqjO4RgOzAk4fVgYGcntmlJ4Sj43kOQPxybO5ufT9Xc8fZ4fru6lONm9OeB40djt6lW6nQTXheNhuGzTOx4zWuw8C/w3u1QMoW8ix7nDVshC9eUUR8MM2NEIeMHWj/mwQebfRJQgB3Isv5i2GyKSad8F12ziBtWP8EvbM8RdGSz97x55Bf2M0m6ZU/CgINQZ93JpP7NqqttNlwzr2Tg+P/jg8+2c0JRJqdOHNBq+ZxSitGHngSHntTy8zrgBDjgBBQwEPhhMMyZ5fUUZ7vplx27Ex5qXMxHzoKXrjUXvk/vgxnfg0GHMBO4/9LpPDPvav5U+RN+W7iQn4ROMnPAJlBdV89s+4fUjT6bnH7jUMBNsydy/j3VzOZ1fuN5mpzCXze177OHTRz/uN+YC0Uio44jg0m8d99iFq0v55SJJeyq9lPjDzf2l29kyiXw7m1cGXmanMOv5ZSl77Fy1HdZdNFFlNb4eW99GeNKcphm5Q+C03/ArPdu5snSzXBQSZNDeRt2Uu0oJP/Q78Ki20wc/fjfmP+XssGki8iu2ts0me/wQMjXWOlxQL8sc2F0JIaGWo6qHVGUyZ/POYiPNlbw8abKJkJQ4w8xSW3CP2Am7mTq6HNKTKjkzd+b2nGHx4RGBh5sYsy+Sug/kYEHnc9ApZg5soBrH/+cP762ipJcD6dan4MvZzhe4Of/fY53OJRc/zbecoSYMOlQrjp/svkOjjnZlOS+9Qf41vNG6BbeAvVlMOmi1u3rPwHCPnIbtnLKxAG8+PkOfnXaODxOO0s2V3Js3evsKp5JyYCD2nyLM0cWMrQgg+c+28E5Z11v7ug/+CcMOxxWv8wd4Qvx9B+978/KZodz7zUFIf+YCDqCOvse/n7gifxjwToiUc2PTxzT5u5KqcYQZ7sc8h3zl8DKYVfwzf++zz3TZ3Gs1Qep1Kog65/jAW8GZSPPZvb651i5bRvFax5Fu7I5661Cat5YyNNXH9bEi/3T+iHcDNw6aTcl04dQOX8jn0UncJy1fs3u1AnBp8BopdQIYAdwEdA8a/oycK2VP5gBVO8rP9BIwhfFZlP86IQx/OiEtv9pbaIUjDsDxp5uyue8+aAU/YALDh2yz92TOb46579QMAJ73R68x1yPN3axu/B/SR2iX46nsWyxO8hwOZgwMLflCocLLnjUJMcX3Qr9xps7GoujxxRz9Jgr4Mm3cXx0JyMnHsm82gBa60ZxcpetIEc1UDEmnmAbPzCH1354FNtW/oXsd86Hd242sf5NC82FdcVzRqxm/ahVe6cPLyDb42Dh2jJOmVhi5p0Fxg1o2sIAhxuO/Alq7s+4qHQlePOZcMFvwWVneFEmw4sym2zuGXcKvHczWXs+Aw5rXB6NarKC5fiy+5Gf1c/c3S5/Bo690eQwRh4DOSVkuWupD0aIRLVJ5js9EA5QWuPHYVMUZ7utXkPJtZiYMjSPDzY0Ha1aX1/HYFVObfE+7m4TsTvh5D/Dsb8yo5kdrjY3ddht/POiKVx838f86KkvyPI4OHR4AdfOr+UB4LC8KhwD+zGp4SvYAEfMPDx+p+7Ng1P/Ci/9AP41ybTZqNxoBPmA41s/Yclk87jzCy6YdiwvfbGTZz/bzrdmDmP+22/xa1sZwRk3tvv2bDbF2VMGcuc7Gyj1TKb/5Ivgo7vg0wfwF4zl3p2n87fizHaP0Uj/CXDpK/D5IzD8KJh0PrnA789qmWDuTvrnZRLARWmsXBUanxdnmxsH24yr8Gx4Es+Su2Dri6zufyYrNoaAEHcv3MBNs01yfvn2ah5fGeCHReMpKX0X9l5CQbSCtwOjGFFeT0mup/3qo3boshBorcNKqWuB+Zgb5ge11iuVUt+z1t8DzAVOAzYADcBlXT1vp1EKMlpJbnUHrgw48aav59jdTVYxfO99Mxhl0CHx9smJHP87WDuTM7f9nX+Gv0NDMNJYjZRVaUIw9kGHNNlleFEmw48+EWovN678kgdN4tSTC4ddYy6yrcWUMRerGSMK+Nhqtha7uzmwuRAATLvcJBe3fGAuUp5WBM/CVnwgERRZNU1LTivqgxRTSTjDEt+DzjdVaB/8w9TXH2fi0dkeY2+dP0xuhtNKFvvY2xAiL8NpxDESTLrp3PDCTJ5fugN/KNJYRaYqN2FTGme/TtzktPa/awWP0879357Gxfd9zKUPfkJxtps9tQH8OYWcO7SBc8+aBB++a36lRc3utA++xIzc/fIpI3qHXwtTL209Lg8mvGZ3w64vOPyk85gxooDb5q9lQI6HzK/mox0K1/jT92nz7IMHccfbG3j5i51ceeJNULEBQn4WTfgL4Z01jGvuLbbHkEPjoa0eIuaJ704YR1Ja6yc/w9n4v+83aiof6IOYtcm0h7iz5kgmD8ljVHEmLyzdwQ2njsPrsnP3wg1kexzkTzkT3r8VVj4PwCfRsYxeV8bUofl0toddt4ws1lrPxVzsE5fdk/Bck5CrFHoJTm/TKpnmFI+B43/LqDd/x/n2EextOK5RCHKrV1Ghs8ksbMObOvVWU5pYtRVGHm08gVjopB1mjizkzdV72FXtY9WuGoYUeFuWdIJx98+6I5l3CU4Pu9QA8uqb9nTZU+unRO2lLssq6Rx7Bjh+YkYnZxSa19AYF67xh4wQOM04gmpfMJ7wjYSSFoJhhSbZuH1vAwf0MyLnqtoIgLukAx5BJyjMcvPCD2bx9wXr2FRWxy3nDsfz4RiTEwETKswsbv1m6cBTzF8y2J0wYCLsWoZSir9+YxJn/vt9vvvIEp73rCQyYDKOrH77PMyo4iwmD87lxS92cOVRR8J33wTgnee/JMvdwMiiJD2CFOFy2CjKcsUHsAGlNQETFrKw2RT/zf8p9sBDDDv4eF5/u4g/zh7E6P7ZPL90B/NW7uKgQXnMW7mba445APekofD+X833tGAkft8Y3ltfRrQLzR9lZLHQPrOuozZnDBfYF7K3Ph7uyK9dzxo9FLezjfYBdoe5azztVhOOS0IEAA4bZdpUfLSxghU7qpnYWmirE1Tai8gMNg3HVDcEyaMeZU1AgicHTv+bEbAz72gsVcyKeQQBa7IZa2Tx3rogebGxJZHENtTth4YG55sige1746OTvXVm9LK9cGSX3mcyZLod/OaM8Tx02XSOHdvP5OIqNpiV5evj03F2lcGHmvr3cJDhRZk8esUMrj58EFPsm3CMOCLpw5x98CBW7qxpTIQGw1HeXL2HI0cX4Wir1Ugvon+Op4lHsKfGT78EIQAYNHQUV/qu4d7ACTjtijMmDWT68AKGFHh5Zsl27l64AbfDxmWzhkO/cTDxG2bHo3/JkWOK+WhjBR9urDAJ6E7Q+z9FIbUoRe3wE5mq1lNdU9W4ONu/kx2q9YR2Vxg3IIe8DCfzV+5mS0UDEwd1jxDUOwvICjdN4NbX1ZhwTEZefOHBl8APPzf5JAuvJXb+kNWSwyoTrWtoaN0jcLjjU1e2QlGWEcXyurjX4PLtoR4vuJNISnY3hQdAXanpOFq+rvuEYNgs06pi5+cATBmSxw2TfdgiAZPwTZIzJg3EblM8vtiI5eOLt1BWG+ie3F4PUJztpjyhE+nuGj/9s5veGE0bVkCtP8xDH2zm+LH9yc90YbMpzps6hA83VvD80h3838xhFFrfHc69H366FiZfxFFjiqm3mj/OGNm5sLcIgbBvhszArjTssHrSBOvJCu+lzL5v176j2GyKGSMKzCAc4KBuEgKfq5CcZpU8sQnLXVl57e7rdpqfiT9ktW+w5uj1++rJy4h5BPE6/ObjCABTybNuPmz9OEEI4hcHb6CcStvXlLvaF7HWGdsWm2qj7hQCgC3vx5dt/dA8DpmZ9GGKs92cf8hg5ny4mSvmfMrNr6/hyNFFHDOmeN879wKKstyU1xrRj0Q1ZbVNQ0MAJ07oj8MqFf/OrOGNy781cygjizM5sH82PzgmoZDEZoNs02nh8FGF5Fhea5tNHPeBdB8V9olj4BQAPHvXmgVVZmxguX1AG3t0jcNGFjJ/ZSlOu2L6iO65OAY8RWRW+0zfFutCHrCEwJOZ1+6+bXkEkWADWW6HuciHm4eGmuUIlj8Dz18JQObsu/E6CymvjQtBZqiCGnuqhMBKDK96yTzuqy4/WTILoXic6cV15E/Nsm2fGKHJLGx/32b87swJKKV4d+0eThzfnz/Nntjt3ujXRWGWi4r6gDWRjhlM1r/ZwMYcj5Onrp5JaU2AmQldfAuz3Lz546MBWowpipHtcfLStUfwVXkdh48q6pSNIgTCPvEUlFCnPXhrrGRrlXHRK10l7ezVeWZPGcR768s5aUL/xsqKrhJyWxfZ+nLTogQI1Vvz1O7DI/C6jA2+ZkJAyGfWRSOAbj9Z/NkcyBlk8g+v/pjJmf+grC7edyg3XMFWT8umaT1CwQhTfrrsCfN6wOTuO/bwWaZ0OBIy7Si2LW5MwncEr8vOX85te8xBb6Y4y00ooqnxhSmtMeLfPDQEcMiw1m8E2hKAREYUZZr+Up1EQkPCPsl0O9msB5BZb7VjrjaPNa7WG6h1lfxMFw9851AuPDS5iWuSQcdG/SZMMRn2GSGwedovQfQ4mnkE1vwGtmiQDKc9ftFPFAIdbezESiRswmrjZ8M37gel+DavUtUQDx/lRffic6co1OFwm7r/aNgM9Ovg3Xq7DJsFwTqTJyhfb8bwdLJDZroSCwWW1QUak8bNQ0OpRoRA2Cd2m2ILJeQ2WA2canYSwYbP3Tk3NBUoa9L5UENVfGHAan/tbmWcQgJteQRuQmZdLDGcGBqCuEBUbDBJ0wGTTFx3zMkcFvyYWp+1X6CWDPwEvCmMeY+2WlJPOKd7jzvqWONtrHoJNr5tlnUgUbw/kJgTKu2lQiChISEpdtsHkhtYbFz8mp1U2grxutseydrr8Jqkc7C+mtioBHvIGoW5DyHwNOYIrGRxghBkuBzxxHDjyOJYJVHQ5CP2bjavY4O0Rh1H/qqXyGgwTQF17W4UEPJ2f/I9aWb9yCSNW5u2sit4880Yki+fNm0x+k805ap9iKJs830orwuwp8aPTUFRVu/67YhHICRFuXMANqKm9XDNTvaowsY75XTAboV/ggkegSNkhYn2KQSxqqFmHoEKkeFKCA0l9hqCuEDUWm21cqzmeVb7hUE+M9I5uNd0W9GZX0+oLSmcHtOAMMlRyh3isGvMBPK7lsGkbpp3OY2IeQQVdUF21/gpynL3uvEP4hEISVHlKoEgJlFcs5NSXYTXmT5fH0eG8QhCDdXxZeHkPAKX3YZNmcl4zI5GCDwEjRjGjtNWaKhmJyi76RAL0G88GsWQsEm++6p24QZUdgo9gq+TEUeatuJ7N8OhV6Tamh4nP8OFTcVCQy1LR3sD6fNLFlJKjbsE6oCqLVCzk116DF5X77qraQ+nJQTRBCFwhesI4MZtb6WFRQJKKbxOe4JHYO7wTGjI3n5oCIwQZPU3bTGs/evc/RnUUIo/FCFUY8ZMOHNT6BF83Zxyc6otSBl2m6Ig092YI4iNLO9NpM8vWUgpvowSoigoXQmherZH8k18PE1wZ+QQ1YqIPz4/sjPSQMCW3N2Zx2lvJVkcbBoasrcRGvLtbdG7pz5zCEPVHmp8IcK1ZUS0wp2dPsl3oWMUZbkoqw2ypzbQor1Eb0CEQEgKj8dDhSqArR8BsC2c3201/j2B1+2gATfRQLxNrz0SIJikELgcNoLh2MjieI7A02r5aLPQkL+mRXfUQPZQhqpSavwhdN0eKskmJyO5fkxC+lGc7WZXtY/K+iADRAiEdCXT5WAnxY19Y7bpYnM3nCZ4nXZ8uNDB+ExrTh0gbEvu4uu02whFWlYNuR22BCGIhYasYzYKQXV89jKLSO5wilUNtTXVqPoyKnRu67NfCfsFhZkuVu403mhJrgiBkKZkuh2s1/GZ39bqIY2tF9IBt8OGX7tR4XjHT0c0eSFwOWwEG4XA7OMhiMveWtVQzCOwQkOB6hYegSoYBkCwYgsOXwUVOqf1dtvCfkFxwkjirowA/roQIRCSIsNl5+OwaYEQySjGjzu9hMDyCAjFhcAZDRBJNjRkTwgNJXgETocyE79D05HF0NQjaDZ62VlohEDv3YIrUEE5ueR40yfnInSMAwfE///DCkUIhDQlw2VnXngq0bFnUXrULQBpNY7A7bDhw40twSNw6QARe/IeQSAmBDYHUWy4VQinvbXQUIIQRKOmrUUzj8BbZCaQV9Xb8AQrqCQnrYRV6BgHD81rfN7bBpOBlI8KSeJx2qkjg7qzH6S8vB74IK0uXK0JgVMHiTo6kSxWiojdjTscwtUkR9BKaChYZ/oONcsRZBUNJKAdeGs24o40UOfIT5tumkLHGVWcxY2njWPMgOxe+X8WIRCSInb37w9GaLAGVqVTstjtsOPTLmwRIwThSBQPQXSyHoHdRn0w3Pg6rNy4CeJqzyMIB0xYCFp4BG6nk80UUVS1HIAGZ4paUAs9xpVHff2zz3UWCQ0JSRG7+/eFIo319J40EgKnXeHHjd3yCPzhKB4VRHfAI2isGgIiNhcemoeGWskRxBrbtdLhdI+tmIH1qwCocg9qsV4QegoRAiEpmghBGnoESikCyo0jYro/BkIR3HRACBKTxUDI5sKjgthtqv2qoTY8AoByR3xin1pveky7KOyfiBAISRG7+/cF40KQTjkCgJDNgyNqhMAfjuIhhHIkN9y/SY4ACCsXXpsVKgq3kyyOjWR2txSCnZ741IOhVDacE/o8IgRCUiR6BA1WaCidqoYAQjY3zmjcI/AQbBwlvC+aC0FIufAoa5xAe6GhdjyCLVkHA7BSjSbL2/sGGQl9B0kWC0kRCwP5QxH8aeoRBG1enOEAaI0/EMSpIihnBzyChBxBUJkcAZDQdK6V0FBjjqClENTkjuGGyl+xuGEQx2XIYDIhdYhHICRF7KLfEIwni9NNCMI2j5lTIRwgGDCtJpQrSSFoniPAhbfRIwgAKt5dtIlHUGWet5IszvU6eb5+EptC+RT0wtpyoe8gHoGQFLEGcz6rfNRlt/W6yTX2RcQem3S+gZDfCIGt0x6BEzcJoSG7C2L14bGkccwjcHjiyxIoznI3DlIryBAhEFJHev2ShZThTQwNhSKNs3alE5FYYjjkI2R5BHZ3RlL7NvcIArhwKys3EAk1vdDb7KBsxlNopeFcjMQJSgoyRQiE1JF+v2YhJTRJFgfDaTUXQYxEjyASE4JkQ0MOG1FtBqIBBHHiinkE4QDYmn0eDi+E/VafoZb5AYD+uSIEQu9AhEBIinhoKIovFE27iiEAneARhC0hcCTrETjMTyUWHvLjSggNBRob0TXi9JgGd4GaVvMDAP1z3AnPpWpISB1duq1TShUATwHDgc3ABVrrva1stxmoBSJAWGs9rSvnFXoeu03hctisAWXhtEsUA0QThCAaNCOMkxUCp5UPCYajZLgggBOXtkJD4SA4mt3ROzOMELTjESS2Ix6U1/umLxT6Dl31768H3tJa36KUut56/cs2tj1Wa13exfMJKSQ2b68vFElLjwCnddEP+4gEu+gR6EQh8Lf0CByWR+CvgdzWRw27HXb+ddEUAGy23teITOg7dFUIZgPHWM8fBhbSthAIaY7XaW+sGspyp1+OoHHwWMhHxPIIXO7k7sTdCR4BWEJAELSOVw01OZd3nx4BwOwp0mNISD1dzRH011rvArAe+7WxnQbeUEp9ppS6qr0DKqWuUkotUUotKSsr66J5QneS4bLTYPUaSqf5ihuJlYqGGtAxIfB00COwhMAXTZiXONxajiADQg3t5ggEobewz9s6pdSbwIBWVt3YgfPM0lrvVEr1AxYopdZorRe1tqHW+l7gXoBp06bpDpxD+JrxWB6BLxRJq4ZzMZTTismHfOhwLEeQ3GxRzUNDPm0JQdhvCUGzcQJOD/iqzPp2PAJB6A3sUwi01ie0tU4pVaqUKtFa71JKlQB72jjGTutxj1LqBWA60KoQCL0Xr8vKEQQjaZkstrnjyWJCpudQizv5NnA1Cw35tPXTCflN1ZAzr+kOzgwoW2uet9JwThB6E10NDb0MXGo9vxR4qfkGSqlMpVR27DlwErCii+cVUoDXaW9sQ52OoaHYKOJo0AexmcqSHFnstDyC2JwEvqglBI0eQSvJ4trd5rl4BEIvp6tCcAtwolJqPXCi9Rql1ECl1Fxrm/7A+0qpZcAnwGta63ldPK+QAtI9NGR3mXxAJFgf9wiSFQK7qeoJhk20siGWIwgHLCFopXwUK7IpOQKhl9Ol0g+tdQVwfCvLdwKnWc83AZO7ch6hd+B12anxhwhHdVqGhlwuFwHtQAfqsUU6FxqKeQQN+/IIEgVGPAKhlyMji4Wk8TptlNUGAMj2pF/5qNtpw4+LaNCHLezHT0KjuH3gbCYEdZEEjyASaFk+mugFiBAIvRwRAiFpvE57Y7fMbE/69c93O+xGCAIN2CM+gir5tg7NhSDuEfhaH1DmTZiMPrO4S3YLwteNCIGQNImT1Welo0fgsOHTbqKhBuyRAEGVfKM3l8N4DqGIJhrVCUIQaL3FREaCECSKgiD0QkQIhKRJzAukZWjIYcOHC4I+HFE/IVvnPIJgJEoAyyMK+VpvOufNjz+3p99nJfQtRAiEpMnzxsNB2e70Cw15nHb8uNFhH45IgLCt5WQxbZEoBKFEIQjWg47Gp6mMIV6AkEaIEAhJU5gVv9ilrUegXRDy4Yz6CXfCIwhGNKGIxq+tUFBsTuLmoaHs/t1hsiD0COn3axZSRmHCvLppmSNw2qnEhQr5cOkQYXtW0vs2lo+GYx6B9Vn4q81j89BQ3jA45Dsw+qRusFwQvl7S79cspIziBI8gPw3n2I3lCGzhWpw6Gp+xLAmcjcniKMFwQmgoJgTNy0eVgjP/1R1mC8LXjoSGhKQpzo4LgT0N++e7HTb8uFERP24dINoRIWgrWeyvMo9JDkwThN6IeARC0uRluDjmwGLGl6RnywS3045Pu7CHfbix0dCBi7fDEj6TI4gSxIFGofxt5AgEIY0QIRA6xJzLpqfahE5jPAIX9ogfD/b4HMZJoJTCZbc1hoZAEbW7sbeVIxCENEJCQ0KfweQI3NgjAbwE4zOWJYnTrhqTxQBRuzshR5B8Kaog9DbEIxD6DG6HCQ3ZiJChIh3yCMC0ojYegekqqu1uaKiwDp58BZIg9DbEIxD6DE67wq/id+4qyRbU8f1tBCO6cZayqDMTaneZla7kZjoThN6ICIHQZ1BKEUkcTdzBi3csRxCyGu9F3TkQDVsrxSMQ0hcRAqFP0WQ0sSevQ/s67aqxxQSAdmXHV7qz29hLEHo/IgRCnyKcMHZAefM6tK8zVjUUEwJ3QhmthIaENEaEQOhT+B3xO3fl7diEMU67jWBYN05gr2JegLJL+aiQ1ogQCH2KBke8K6gjs2MdQmNVQ6FIs7mI3VlJz3QmCL0REQKhT1Hvil/8XVn57WzZElezHIEtNgWlU8JCQnojQiD0KULOeDjIm91Bj6DJyGJQWUXdapsgpAoRAqFP4XTGJ9TxZnYiR5AwjsCeP8SsiIa6zT5BSAUysljoU2S4HHw/eB3ZqoE/JczBnAxOu61Jiwn7EKvv0vSru9tMQehRRAiEPkWO18Hz0Rk47YpbHR1ziF0O1RgactgUtux+8MvNHR6PIAi9DRECoU+R4zGhoQxXx7/6sRxBKBJtnJ+gyST1gpCmSI5A6FPkWFNsepwd/+obITBzFrs66E0IQm9Gvs1CnyLHazyCqO74viZZbEYWN3oEgrAfIN9moU8RCw11ZqbN2DiCYDiKyy4DyIT9hy4JgVLqfKXUSqVUVCk1rZ3tTlFKrVVKbVBKXd+VcwpCV8jxmtCQw9bxr74joWpIQkPC/kRXv80rgHOBRW1toJSyA3cBpwLjgYuVUuO7eF5B6BQjiky76GGFGR3eN54jkNCQsH/RpaohrfVqMH3e22E6sEFrvcna9klgNrCqK+cWhM4wpn8Wfz5nIqdMGNDhfV12ZXIEYRECYf+iJ77Ng4BtCa+3W8taRSl1lVJqiVJqSVlZ2ddunNC3UEpxyYxhFGZ1fI7h2MW/IRjBKaEhYT9inx6BUupNoLXbpxu11i8lcY7W3IU2aza01vcC9wJMmzatE7UdgvD1ELv41wcjuMUjEPYj9ikEWusTuniO7cCQhNeDgZ1dPKYg9Dgxj6A+ECYzu+MehSD0VnrituZTYLRSaoRSygVcBLzcA+cVhG4lVjJaHwjjcXasT5Eg9Ga6Wj56jlJqO3AY8JpSar61fKBSai6A1joMXAvMB1YDT2utV3bNbEHoeWIeQa0/jFtyBMJ+RFerhl4AXmhl+U7gtITXc4G5XTmXIKSamBDUiUcg7GfIbY0gJElipZB4BML+hHybBSFJEttKiEcg7E+IEAhCkiQOIhOPQNifkG+zICRJEyEQj0DYjxAhEIQkSRSCzsxnIAi9Ffk2C0KSeBPmOHY7xCMQ9h9ECAQhSbLc8Yu/eATC/oR8mwUhSTLd8WE34hEI+xMiBIKQJIlCIB6BsD8h32ZBSJJMV4IQiEcg7EeIEAhCktgTJjrO8TpTaIkgdC8iBILQCWJzHwvC/oAIgSB0ghyPeATC/oMIgSB0gv45nlSbIAjdhvi3gtABnrhyJhvK6prkCwQh3REhEIQOcNioQg4bVZhqMwShW5HQkCAIQh9HhEAQBKGPI0IgCILQxxEhEARB6OOIEAiCIPRxRAgEQRD6OCIEgiAIfRwRAkEQhD6O0lqn2oY2UUrVAmtTbcc+KALKU21EEoid3YvY2b2Ind3HgVrr7I7s0NtHFq/VWk9LtRHtoZRa0tttBLGzuxE7uxexs/tQSi3p6D4SGhIEQejjiBAIgiD0cXq7ENybagOSIB1sBLGzuxE7uxexs/vosI29OlksCIIgfP30do9AEARB+JoRIRAEQejj9EohUEpdp5RaoZRaqZT6UartiaGUelAptUcptSJhWYFSaoFSar31mJ9KGy2bWrPzfOvzjCqlekX5Wxt23qaUWqOU+lIp9YJSKi+FJsZsas3OP1o2fqGUekMpNTCVNlo2tbAzYd3PlFJaKVWUCtsS7Gjts/y9UmqH9Vl+oZQ6LZU2Wja1+lkqpf6fUmqt9Vu6NVX2JdjT2uf5VMJnuVkp9cW+jtPrhEApNRG4EpgOTAbOUEqNTq1VjcwBTmm27HrgLa31aOAt63WqmUNLO1cA5wKLetyatplDSzsXABO11pOAdcANPW1UK8yhpZ23aa0naa2nAK8Cv+1po1phDi3tRCk1BDgR2NrTBrXCHFqxEfiH1nqK9Te3h21qjTk0s1MpdSwwG5iktZ4A3J4Cu5ozh2Z2aq0vjH2WwHPA8/s6SK8TAmAc8LHWukFrHQbeBc5JsU0AaK0XAZXNFs8GHraePwyc3ZM2tUZrdmqtV2ute9Uo7TbsfMP6vwN8DAzuccOa0YadNQkvM4GUV1208f0E+AfwC3q3jb2KNuz8PnCL1jpgbbOnxw1rRnufp1JKARcAT+zrOL1RCFYARymlCpVSGcBpwJAU29Qe/bXWuwCsx34ptmd/4nLg9VQb0RZKqT8rpbYBl9A7PIIWKKXOAnZorZel2pZ9cK0VanuwN4RX22AMcKRSarFS6l2l1KGpNmgfHAmUaq3X72vDXicEWuvVwF8xIYJ5wDIg3O5Own6HUupGzP/9sVTb0hZa6xu11kMwNl6banuaY91I3UgvFakE/gOMAqYAu4C/pdSatnEA+cBM4OfA09Zdd2/lYpLwBqAXCgGA1voBrfVUrfVRGLdnn4qWQkqVUiUA1mPK3cV0Ryl1KXAGcIlOj4EujwPfSLURrTAKGAEsU0ptxoTZliqlBqTUqmZorUu11hGtdRS4D5Mf7I1sB57Xhk+AKKYJXa9DKeXA5ASfSmb7XikESql+1uNQzJtJStVSxMvApdbzS4GXUmhL2qOUOgX4JXCW1roh1fa0RbMChrOANamypS201su11v201sO11sMxF7KpWuvdKTatCbEbKYtzMOHh3siLwHEASqkxgIve24n0BGCN1np7UltrrXvdH/AesAoTFjo+1fYk2PUExnUNYX5UVwCFmGqh9dZjQS+18xzreQAoBeb3Ujs3ANuAL6y/e3qpnc9hLlhfAq8Ag3qjnc3WbwaKepuNwKPAcuuzfBko6Y2fJebC/z/r/74UOK432mktnwN8L9njSIsJQRCEPk6vDA0JgiAIPYcIgSAIQh9HhEAQBKGPI0IgCILQxxEhEARB6OOIEAiCIPRxRAgEQRD6OP8fhR4YerLpIgEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(trial.tracking.time, vleft)\n", "plt.plot(trial.tracking.time, vright)\n", "plt.xlim(xlim)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Perform thresholding\n", "\n", "Here we use standard deviation based thresholding, i.e. we detect saccades if the pupil velocity exceeded `5*std`.\n", "\n", "We consider the velocity to be \"above-threshold\", only if both the left and the right pupil velocities exceeded the threshold." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " left pupil threshold: 0.13 deg/s\n", "right pupil threshold: 0.13 deg/s\n" ] } ], "source": [ "def threshold_values(x, threshold):\n", " ret = np.array(x, copy=True)\n", " ret[np.abs(x) < threshold] = 0\n", " return ret\n", "\n", "threshold_left = vleft[trial.tracking.time.values < 1].std(ddof=1)*5\n", "threshold_right = vright[trial.tracking.time.values < 1].std(ddof=1)*5\n", "\n", "print(f\" left pupil threshold: {threshold_left:.2f} deg/s\")\n", "print(f\"right pupil threshold: {threshold_right:.2f} deg/s\")\n", "\n", "# take above-threshold motions only\n", "vleft_thresholded = threshold_values(vleft, threshold_left)\n", "vright_thresholded = threshold_values(vright, threshold_right)\n", "\n", "# detect as above-threshold position only when both the left and the right pupil velocities\n", "# are above-threhsold.\n", "#\n", "# there are some \"negative spikes\", possibly due to small unsynchronization of pupil motions\n", "# we take both negative and positive spikes here\n", "vavg = np.sqrt(np.abs(vleft_thresholded*vright_thresholded))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqQ0lEQVR4nO2dfbAkZ3Xen9Mz9+4uSEIYrYRYLZZiLxiFIEIWmT+MA3YIEiYokJhIJkCIiUopRJGqOKAUFScpKuU4xHFCLFAUR1nIBwIMtgVZLNsUBmKQ0QokISEEG6FIixRpZX1r996Z6T75o/vtr+m5093TMz2n+/lV3bp35vbMnOmZPn36OR+vqCoIIYR0F69tAwghhCwXOnpCCOk4dPSEENJx6OgJIaTj0NETQkjHGbb1wmeccYaee+65bb08IYSY5JZbbnlEVfdWeUxrjv7cc8/FkSNH2np5QggxiYj836qPoXRDCCEdh46eEEI6Dh09IYR0HDp6QgjpOHT0hBDScejoCSGk49DRE0JIx6GjJ4TM5a4Hn8Sf3fPnbZtBatJawxQhxA4X/4evAQB++GtvgIi0bA2pCiN6QkhpToz8tk0gNaCjJ4SU5uSYjt4idPSEkNKcZERvEjp6QkhpGNHbhI6eEFIaavQ2oaMnhJRmixG9SejoCSGl8QNt2wRSAzp6QsiOqCbOfewHLVpC6kJHTwjZkbGfOPqJz4jeInT0hJAdSUfxE0o3JqGjJ4TsSDqKnwSUbixCR08I2ZFROqKndGOSuY5eRK4TkYdF5I4Z/xcR+YiIHBWR20XkFc2bSQhpi7R0w2SsTcpE9IcAXLTD/y8GcCD6uRzAxxY3ixCyLqSdO8srbTLX0avqVwE8usMmlwD4hIbcBOB0ETm7KQMJIe2SrroZ09GbpAmNfh+A+1O3j0X3TSEil4vIERE5cvz48QZemhCybDJVN5RuTNKEoy9ahaDwtK+q16rqQVU9uHfv3gZemhCybCjd2KcJR38MwP7U7XMAPNDA8xJC1oCMdMOqG5M04ehvAPCOqPrmVQCeUNUHG3heQsgaQOnGPnPXjBWRTwJ4DYAzROQYgH8OYAMAVPUaAIcBvAHAUQAnALxrWcYSQlZPpryS0o1J5jp6Vb1szv8VwHsas4gQslZkNXpG9BZhZywhZEc41Mw+dPSEkB0JAiZjrUNHTwjZEV851Mw6dPSEkB1J185zTLFN6OgJITsSpCN6lleahI6eELIjad/OZKxN6OgJITvikrGbA4919EahoyeE7IhLxm4OvUwFDrEDHT3pFI8+M8J7P/ltPHFi3LYpncElYzcGwqFmRqGjJ53io18+is/f9gA+c8v98zcmpQhSET2rbmxCR086xYmxDwAYekXTs0kdXBQ/9LxMBQ6xAx096RR7NgYAgO0JywCbwgXxm0OP0o1R6OhJp3BxPB19cwQpjZ4RvU3o6EmncBryViThkMVxVTcbA4919EahoyedwkXyW2NG9E3h5JrNoZeZe0PsQEdPOsUocvTbE0b0TZFIN6yjtwodPekUzsEzom+OuGFqYD+if/ipLVx67TfwwOMn2zZlpdDRk04RSzeM6BvDRfHDgZiP6L9w24O46Z5HcfWXj7Ztykqhoyedwkk3PpOGjeGrYuAJBiLmG6Z2bYQur2/VQ3T0pFM46YYLZDSHHwADEXie/REIT21NAIQyVJ/o17slncdJN1zyrjkCVXhe6OytR8InR2EgYD3XUBU6etIpnHTDiL45/EAxEMGgA0PNXO7mxKhfORw6etIpGNE3jx8ovEijt+7ot6NqrJN09ITYxWn01h3SOhG4ZKwn5iUP1zHdt0CAjp50ili64dqmjRGowhOBJwLrilhyxWf8jVSEjp50ChfJ9y1iWyZ+AHgiGHj2r5T6WpVFR086hXNEfTuQl0kQKAYeMPDsLzziAoDxxPb7qEopRy8iF4nI3SJyVESuKvj/c0Tk8yJym4jcKSLvat5UQubjHBGnLDaHr1HVjWe/0ch19o4o3WQRkQGAqwFcDOB8AJeJyPm5zd4D4LuqegGA1wD4DRHZbNhWQuaSRPS2HdI6EXSo6mbS0yu+MhH9hQCOquo9qjoCcD2AS3LbKIBTRUQAnALgUQCTRi0lpASxo+9ZxLZM3AgEz7M/6ybO4VC6mWIfgPRKy8ei+9L8FoCXAHgAwHcAvE9Vp440EblcRI6IyJHjx4/XNJmQ2cQHsnGHtE64hqlhB8orXSQ/ZkQ/RdEqy/lP+/UAbgXwAgAvB/BbInLa1INUr1XVg6p6cO/evRVNJWQ+E0b0jROOQAgjeuuSmPPvLK+c5hiA/anb5yCM3NO8C8DnNOQogB8C+KlmTCSkHGlZgcnY5ohHIIh96SaO6CndTHEzgAMicl6UYL0UwA25be4D8PMAICJnAXgxgHuaNJSQeaSjzb5dmi8TP0CYjO2AdNPX8tvhvA1UdSIiVwK4EcAAwHWqeqeIXBH9/xoAHwJwSES+g1Dq+YCqPrJEuwmZwh3Enthv7FknwhEIYdOUalKFYxEXDLgO6r4w19EDgKoeBnA4d981qb8fAPDXmzWNkGq4aHPXcICTYx+qirAQjCxCOhkLhPvZK0zdrT997ZxmZyzpDG5Vqd3RKkKM6pshiE6YLoq3vF/7Kt3Q0ZPO4A7e3RuD6LZdh7ROpKdXuttWSUf0avh9VIWOnnSGRLoJv9Z9K6FbFumqG3fbKumTf58CATp60hmcA9o1jCL6numwyyIIEC4l6CJ6w+fP9EmqT4EAHT3pDM6x73IafY8uzZeJn5NuLOvbGUffo1p6OnrSGZx2vDlgMrZJ/CBaeCRVdWOVvvZa0NGTzuAO4s0hHX2TxMlY6YJ0EyRloj36ftDRk87g09EvhaI6eqtMAo2T9X36ftDRk86QJGPDr3WfqiqWiR8kQ82ApF/BIn6gvQwE6OhJZ0gi+kHmNlmMILXCFGA7ovcDTaqyevT9oKMnnSHW6JmMbRQ/iBYe6UAdfTaiN5xsqAgdPekMsXSz4aSb/hzIyyTQZHpleNumo1fVnEbfskErhI6edAY/F9HTzzdDmIxFnIy12ojmLkT6GAjQ0ZPO4A7cJBnbnwN5mQQa1dGL7Yg+HwhYlqCqQkdPOoPz630sn1smbv78wHj9eb78lslYQgziIvg+ls8tEz+qurHeGZtc8YVVN9aXRawCHT3pDGyYWg5uKcFhPNTM5n7tc58FHT3pDFPTK3t0IC8Tt5SgG4Fgdb9O4qqs/vVZ0NGTzjAV0RuVGNYNNwLBMx7RB7lkrNUTVh3o6ElnmGqYMloGuG5MJWONnkDzQ++snrDqQEdPOoMr+0vqpPtzIC8TPx6B0I2qG2r0hBjGNfLEDVNGI891w41AsL6UYKLRcwQCIWbpc530Mgm0W3X0u6jRE2KXZHFwV1XRn4htmcTJ2K50xvaw/JaOnnSG6RWm2rSmG6jq1FAzq/s13zBFR0+IQXw/O+uGEf3iuOA9k4w1HtH3MVlPR086g6um7OMY2mXhnLonSEX0Nndsn9croKMnncGfmnVj0yGtE84ZepmqmzYtqk9AjX5nROQiEblbRI6KyFUztnmNiNwqIneKyFeaNZOQ+eQ1+j5dmi8Ll3gdeAIv8hZWG43i8soeavTDeRuIyADA1QBeB+AYgJtF5AZV/W5qm9MBfBTARap6n4icuSR7CZlJ0OMDeVm4fTgQwdCzPVqCDVM7cyGAo6p6j6qOAFwP4JLcNr8E4HOqeh8AqOrDzZpJyHymq276cyAvC6d+eamI3qqDnP5+GNWgalDG0e8DcH/q9rHovjQvAvBcEfkTEblFRN5R9EQicrmIHBGRI8ePH69nMSEz8AOFpJe8M+qQ1gkXvQ8kmV5pVboJchG91VxDHco4eim4L/9JDwH8FQC/AOD1AP6ZiLxo6kGq16rqQVU9uHfv3srGErITfqAYdmBu+joRSzcd6Iztc0Q/V6NHGMHvT90+B8ADBds8oqrPAHhGRL4K4AIA32/ESkJK4Afh2qYDRvSN4ZKxoXRjvTM2W5XVp+9HmYj+ZgAHROQ8EdkEcCmAG3Lb/D6AV4vIUESeBeCnAdzVrKmE7MwkiuhFBJ5Ujzy/dd9jeMtH/xRPnBgvyUJ7ZJOxtk+g6Tp6qfH9sMxcR6+qEwBXArgRofP+tKreKSJXiMgV0TZ3AfgDALcD+CaA31bVO5ZnNiHTuCmLADD0vMrVIR//+r341n2P44/vemgZ5pkkXUfvGZ9emZahhp6YfR91KCPdQFUPAzicu++a3O0PA/hwc6YRUo20ox/UOJDjhiCj0sQyiOvoU5KY1dxHPt/QJ0fPzljSGSaBYhDVAA48iefTV2V77DdplmkyztH4iXCSieg9sxJUHejoSWcIIo0eCA/myknDqL7s6W06ekc+GStiP6Ifel6tHI5l6OhJZ5hkNHqJx9KWxTmwE6NJ47ZZxdWau2h+IGI2Es5E9AOPjp4Qi/hBEDt6r4YGe2IURvKjSX/qq+eRSDfhbc8Ts9JNEEf0oUZv9YRVBzp60hl8Tbpi61RVPBNF8tt09DHJULMo9yFiVrqZpCqIwu9Hfz5nOnrSGfwgiJt66kRsz0Ta/KhPvfFzyEf0oYNs0aAFcI59GJWKMqInxCATP5uMrRrRn4ykmzEj+pg4Cpa0JGZz/2Q1ertXJnWgoyedIdDF6ujHUajKiD4hPY/e/aZGbw86etIZ8lU3VR29c/BMxiakRyAAYWRv9TyYraNnwxQhJkl3xtbRYF2DFR19QhDkI3rbdfSeIJqFxIieEJO4McUAammwru6e0k3CJOfoLXeUhkPvQpdHjZ4Qo0yiMcVAWA5Y1SGNo4i+T5f08/BTnbHhb8tjitM5HLsnrDrQ0ZPO4AeK4cB1cFZ32JMokq87I6eLBDmNfiB2te3MFR81ekJs4qeGmg296i3u42j7qqMTuoyfk24sd8b6gWIwSI9y6M/nTEdPOoMfKKLjuFZ5pYvo+xTpzSNfXjn07GrbkyBIrkw8QY/8PB096Q5TY4orHMlBoHD+q0/a7TzyyVjL1SqZhWkGjOgJMUl+THEVqX2cOuip0SfEK0xlImGb+8fPfz+Mvo860NGTzjBJTa+sOrRqnHLufYr05uGkm+wJ1KaDnExp9DbfRx3o6ElnyDRMVVxhapKqne9TpDePeB59ByLhMIdj/33UgY6edAZfs+VzVeq9sxF9fxzAPJxME08FNVxeOclp9FbfRx3o6EnMI09v419+/k48cWLctim18H2tPaZ4Qo2+kDgZm5leaXP/+H7SGTuoUX5rGTp6EvORL/0A//VP78WffP/htk2pxWSBZJtz7pvDfnVMziPpjA1vD6TGWrxrgp+ebir9unKjoycxbh77scdOtmxJPRYZU+xGFO8eekzGpkhG+yZlq1Yj4UznNCN60lfcmqmPnxi1bEk9FhlT7KK7zeEARgPWpZAfU1y1bHWdSM9C4ggE0lvcmqlPnpy0bEk9fD8b0Ve5NHcR/eZAoPT0MXEdvZNuDK8w5QdBLO15XHiE9JUT0ZqpT5w0mozVrEZfpbHHafQbQw9Pbk3wxv/4NXz7vscWsmd74uOztxzD1thf6HnaxM+NQGh74ZFv/vBR/M4tx2o91p+64rN5wqoDHT2JcRH9CaOOaRIkVTdV56a7iN6dKO740ZO4+sv/ZyF7fvdbP8I//sxt+MQ37l3oedokP9Ss7YVH3vqfvoFf+cxteHq7+lVnVqOndDOFiFwkIneLyFERuWqH7V4pIr6I/O3mTCSrwmn0ViPQdIu7V7He29XRbwySQ8I9V13ue/QEAOCRp23mPIDpMcVDz2utMzZ9grn3kWcqPz49C4kafQ4RGQC4GsDFAM4HcJmInD9ju18HcGPTRpLV8EwUJW0bdPSqmh1TXLEhxlXa7Bomh8T2ZLH94E6cjz5j19FPSTctOsintpIo/sEntio/Pj/dlBp9lgsBHFXVe1R1BOB6AJcUbPdeAJ8FYLMIm8SO6aRBR++O2bot7k6jH6YielehURd3ZfRMDZlhXfADhUTrrAL1FnRpisdS1WCP1agMm/jZ6aaM6LPsA3B/6vax6L4YEdkH4M0ArtnpiUTkchE5IiJHjh8/XtVWskRUNdbot8b2klQuIh/WXFjCafQbg8S5y4KO3p0w6+jJ60JaDgPajegfTxUJpKP7sgS5ERlWh7PVoYyjL/q25/fQvwfwAVXdMRRU1WtV9aCqHty7d29JE8kq2BoHcf24RY0+KBi+FShKl0q6y/i0Rr9oB6hrQDMd0atmrmza7IxNR/FPbVWvDMtMr/Q8qLabWF4lwxLbHAOwP3X7HAAP5LY5COD6KAI6A8AbRGSiqr/XhJFk+aTlGovSjYveB6mGGCBbabETSUSfOPpFT3hdiOiDVEki4BbsaMc5pmcw1Ynos/Pow/smgWJzwaS7Bco4+psBHBCR8wD8CMClAH4pvYGqnuf+FpFDAL5AJ2+L0SR0dM/ZsxFHopYoWtsUiObfDOY/Pq6jT50UtieLSVjbkQT2zLa9/enwg+TkCYR5i7ai4CejKH7oCZ6u4eizSwmGnr4vOv1c6UZVJwCuRFhNcxeAT6vqnSJyhYhcsWwDyWpIO/qRH5g7AJy9Lnp3kVtZmcFdEaQj+kWrbroQ0ftBEJ80gXYXHnHFAmedthtPbVeXboIgu/YtgN7o9GUieqjqYQCHc/cVJl5V9e8tbhZZNSM/PIhO2xN+JbYnPp61WerrsRYULXkHlJ9QOCqoo180KZ129Kq6cHK3DdLdxkD1/oQmcVeaZ5y6q5Z0MwmCTMMUEI7N6APsjCUAEpnitN0bAGBOvnEOPT0CASh/IE8Kqm4WjuijfegHurAM1BZ+gExEP2xxzditsY/dGx5O2z2s3Rk78PKBgM3PpSp09ARAIt04R79lzDHlNfqql+bpefSO7QUj+nQyd9TmgJgFCFLL7wHtNhqdGPnYszHAruGg1tXWJLeUIECNnvQM5+ifvSuUa8bGHX3VZNvY1eF7aY1+cenGddouetJoi/RiHUAY3bdVXnlyHMqJezYHtSqifF/jhri+afR09ARAEnGeujty9MYi0MmUo8/eP/fxBRH9IuWVqoqtsY/n7AmvkKxG9H6uvLLNNWNPRtLN7qFX67MZp8YUx9INNXrSJ5KIPqxFtKYp56tmXERfVk92Gn3aqW1Pgtqz6Ud+gEARO3qL84OAAkdfsRGtSbZGPvZsDrBnc1Cr1yM/vdLd1wfo6AmAxNGfsit0TOYien+xiH4cKDYGMtUGPq4Z8W2Nwv13+rMiR2/sxOkIO2OT2206SKfR796oLt2oKsa5WTdAf9aNpaMnABJp4ZQooh8Zc0zJCIO8Rl/ufUz8IKPPO3Rq2kc5XMT5nD2bAOw6+nxn7KBFbfvk2MeezWHk6KtdbTl/vpFaryC8n46e9AjniE6JNXpbB4Bz6Ol54+H95R4/9sPL+nype10/4By9i+itnTgd6RnuQNKn0EZV4tbYx54ND7s3ogR3hX3qrlAHA2r0pMe4g+bZmzaTsfHCIamFR4DyddJjP8DmwMO0eFOPrTiid9KNTY0+jOiT2+7vNiJ6J93s2QivOqv0ejipaSMn3VCjJ71ilIvorUkNM+voK1TdFA0/WzSiT5Kxtvanw9d8HX0kibUQCYfSTajRA8BWhZNnPoczZMMU6SNJMtZqRO/m0UcR26Caow9L77wp6aYuW6OsdGPtxOnwU+vwAohXaGojot8a+dizMawV0U+tV1BxFpJ16OgJgGlHb01TjoeauaobqR7RF1XdLJ6MdXX0RqWbqYi+PckjjOgTjb5Kd2wyIiObw6FGT3rFyPfhCbBnM4yW7EX0xdMrSzdMBUFmGUHH4snYqOrGqHQTLr+X7YwFVh8JjyYBJoHG5ZVAtXUT8rOQPGr0pI+MJgE2h17ccGStkzNOtg2yybayDVNjP5rSmNNu6roBJyskyVhb+9MR5EYgVM19NIVz6q68EqjWhDbxs9JN1UDAOnT0BEDk6AdePALAmnQTrzA1NZ2wfGfsxsCblm5qRq5uKNzpxqtu8p2xXkVJrClcFVOm6qZGRJ//fnDWDekVIz/A5nCATaMRfbxCVM3yucmMJQfruoGtXERv7cTpmOQWB29Lo3eLjoQafVR1U0Wjz603MGyxeqgN6OgJgFBa2JWSbsYTWwdAHNHHl+YVp1f6ATYKqm4W1ehP2T2EJ3alm4mfbZhqKxJ2Ulj9iD57xefeEqUb0iucRj/wBANPzFWJxCMQah7IY1+xMZTphqkFHP3GQLARyWFmHX0QZBZjqZr7aIqsRu+qbqrX0W/UDASsQ0dPACQaPRAeDNZGIEw3xFQ7kN2sm6bq6E+O/Fhi2DUcmJ1eGUpaqYhe2klipjX63ZuDzH1lSDT6nLRHjZ70iVCjD78OmwPPnKYcl88N6h3I46iOPk/dOvrtiR9LDLssR/R+bs3YtjX6lHRTLaJ3C8vkO6dtfi5VoaMnAMKI3q2GtDn0DCZjswfyoOKBPHGdsbn7a2v00ex0ANi1YdnRB9lkrBtqtmqNfpwkYzcGocRYRaOfaqhjwxTpI9uTALs2koje2lKCSUSfj9hKPn7W9Mqa9pwc+9g9DB29xSskR74aqa2qm61RotEDYWRfqzOWC4+QPpPR6E1G9NkW96oR/TgI4oqjNHXr6E+Og1hL3hwO7Eb0gWbm9Lc1I+bEaAIAsWyze8OrVXWTH4FAjZ70ilC6SSJQayMQ/CCASP2GqSAIm4Gkoc7YcACXXSnMMfGDwoVHVi15uBOlyyN5IpUqf6ZXIGNET3rI9sSPD6INg1LDeEZjT1lnoKqFFTeL1NG76HNzIBgZ7YydBNkkddwZu+JI2L2ayxGIVPtskhXIcksJUqMnfcLV0QMuArV1APg5iaHqLBMFCpccqVt1szVOkrGbQ3snTseshqlVF6s4qcidjAVS6bOZNQKBY4pJr9hOVd1sDMRcMnacqw6pWgaoGjqRxurox0kd/ebAsHQzo2Fq5RF97uU8SdaBLYOryso3TLEzNoWIXCQid4vIURG5quD/bxOR26Ofr4vIBc2bSpZJOqKvGi2tA36uOqTqlEWFhhp9Q52xW2lHP/TMjZQAQtkrUBQmY9uqP48jepFa0g01+hmIyADA1QAuBnA+gMtE5PzcZj8E8FdV9WUAPgTg2qYNJctlO5WMRUX9cx0Yz5AYSidjtTiaX2RMcazRDwcmI/p8SSKQXtBltba46icv9SFVkm5yQ82o0U9zIYCjqnqPqo4AXA/gkvQGqvp1VX0sunkTgHOaNZMsE1XNdMYK6ju4tvDzEkPFcbqhHymoo6+xI1Q1l4y1qdEnJYlpSSz8vepI2L2cs0Qqfkn9/FCzFpdEbIMyjn4fgPtTt49F983ilwF8cRGjyGpx0eauVOla3frxtsivhFT90jysumliKcGRHyBQZJKxFuvox7mSxPTfq05iupdz5a+eSKVPJl6BzHNVO4KhJ70ZgTAssU1xMULRhiKvRejof2bG/y8HcDkAvPCFLyxpIlk2zgk5R1+1dG0dCMsAk7hFJJzCWSkZO+P+qriOzd2pWTcWyyvzq3YBiXSyckcfuZx0RF/FhngEwiAr7zEZm3AMwP7U7XMAPJDfSEReBuC3AVyiqn9e9ESqeq2qHlTVg3v37q1jL1kCo1wziog96WYSZBt7gFC+qVJeGTZMTd9flfSkRcBuw5SrVMmuMBX+XnUgkET00e+KNowLZKiBJ1x4JMXNAA6IyHkisgngUgA3pDcQkRcC+ByAt6vq95s3kyyTfETviZirL85PWQRCPbns+wjUSTeL11e6RTJ2p2YH2dToszPcQ9qSblwdfSK9VLHAz0k3QHX5xzJzpRtVnYjIlQBuBDAAcJ2q3ikiV0T/vwbArwJ4HoCPRh/ERFUPLs9s0iT5iB6wKd3klwKskmuYLd1U3xEnCyL6QKNJkAXzdNaVZGxAWrppxxZFtioqjOjLfzbjXHmlew5rAU1dymj0UNXDAA7n7rsm9fe7Aby7WdPIqnALV7vyyqrR0jqQH74FuAO53OPDEQjNVN2M4+acpNMYCJO0lhy9kzvSEb20pdHnTsRV5UU/CBvq0rOMLOai6mLnW0eWRhzRD5x0U39qY1vk56YDLqIv9/gm3607ubjo0eo6vMkM9+mIfuUaPTTnpKtVhuWrsoCwe9ra97wudPQkdvRuHn3VRNc6UCTdoEplhmabceK7a+yH+DWjp3P+xZpMMC5Ixkqs0a/WFtWsbFT1OzoJpnM4Va74rENHT5IRsANXdWNvBIJb8zVNkeOeRZyMnRpTXH0/OAfkXt/9trVHk8FlGUcfR/Srb5hKJ8qryi6B6tT3wTP4Pa8LHT1JRfSRRg97EX2gySAzh1chonfTK6capmp2xrrXT/+2FtE7J5iJpFuUbtIfTlUn7YbWpZGKg9EsQ0dP4ojeJd1ExNwBoKpTFSFVBl8VOQKgXhSetOvHE7ii+23t1CB3ZZL+e+WRcEFVVJXvqKpOBQJVB6NZho6exA7ISR9iMBkbFGjs1SJ6LWyYqmdLcURvbJdOzYBP/71yjR7Zz7eqky76flQt0bQMHT2ZckwtlUovRFAQ0QPlr0w00m6mpZv6Gr3kNXpjPqVoYmRb7yUIdKqOvsr1VtH3o0pVlnXo6EkiNaQck0WZIZ9IDQ/sKg1TBVU3NWzpikZfJN1I/L9Vz7opqKOvGNEXfT+sfSZ1oaMnU47JYiNJkUYfLiBd8vGYUXVTq7wy/B2367c0NmBR3Hq7WemmnQoizTnq6uMLdOo0bjEXVRc6epKSblwy1l4pYFFEHr6PaiMQpjX6GtINpk+c7jUskZywkvviv1uYXpm3o8qJMwgKNPoK3w/r0NGTOOpNHL1F6UaRK6OPJKhyj88n++L7G4jozWr0yAYA6b/baJjKSDeoU0efvc/ilWtd6OjJVHVFBWl7bQhUp2QXoPyBrDpr4ZF6tgDJ/nQnIGsnz3zjF5Dsn1VXq2j+860o3RRr9ByBQHpELN2kVt+x9vXXovJKr7xDClzE2IBGn69WMavR53I34d8tRfTIfjRV5zFpwRUfRyCQXpFUV4S/wwPA1hEwq3yu0vtoooge6Ug4+7TWnEpeggpvhL/aSMbmryyqfLRF0lyf5tHT0ZOpZKxnULuc2RBT4rHxohZoZs3YfGdsYpetnVrUMJU0f632vQSqufLKaiMQ8o8Pn8NeQFMXOnoyVV1hcaiZG0qWpmwyNq1FNzGPfkqjb0nuWJSihqm4vLJl6aZ6Mrao6oYaPekR05oyStefrwtFGj1K6rhui8JZN7U0+vB3ulwVsBc95iWo9N9tLDySvt6qM71yOhCwd+VaFzp6EjfGpMsrrbFIi3uQkW5yydhaY4pdctvZEb2OsZNncWdsOw1TmKqjrzq9cnpMscBeGXFd6OjJdDK2YkXDOlBYdVOyISaZTTMd1S9UR4/sidOiHAbMGmq2+og+v/BItemVMxqmbH0ktaGjJ6kDOiXdGDsAiuroBeVGIDgH3NSVTL4z1mzDVKFG7/63WlvCZGrOjgWlG45AIL0ir8VaXHknyEV8QPmqikxEP+N/VW0Jny85cYb329unQL7qxp20Vh/RZ5OxVatuioeaWbtyrQsdPSmedWPs+68zIrYqb0PQzDx6zUkeTqu3tk/z3wsg3Rm7WlsU2ZNw2AxX4fGzcjhNGGcAOnoylXQLI+EWDapB8Zqg5SK2Ii3aUa+Ovji5bS2iL666aakzNheRV02kFpdX2vtM6kJHT1IjENw9FQXQNWDWgVzmOI6lGxRU3SxUXpk8r7PREvncTfh3+HvV0l7+iq3qhNXioWZceIT0iHzSzWJ98eyGqYp19Pmqm1q2RM+X64y1pgcXDjVbk4apqjbkrwgAm6M+6kJHT2ZIN7YOgOLyuXIabPpE555hkVb/WZ2xtvZotr8gTRvltzpVdVN1emVxw1RfoKMnBWvG2ktSFV6ao5xcUrTJQispuROnl1whAUljmhWKGqbc7VW/E0VBh26l6ZXFgYC1gKYudPSkYKEMm9JN3WSsRrX2IjJVErnIrJvYMcVNRtWfq01mJanbkDzy5ZFVez2KO6ftdSvXpZSjF5GLRORuETkqIlcV/F9E5CPR/28XkVc0bypZFvnSM4uRTqEGWzLZFjdMASnppv51fXc0epekL4joV15108D0yoKGOmv9InWZ6+hFZADgagAXAzgfwGUicn5us4sBHIh+LgfwsYbtJEukKBq29v3Pt8gDUcRWsWEqJv57kfJKZ4dNjb6ovBIA0EL5rUavmzKhcjK2uKGuAeMMMCyxzYUAjqrqPQAgItcDuATAd1PbXALgExqGADeJyOkicraqPjjrSb//0FN43b/7ygKmk6Z45Ontqe7Hp0cTU5/PyA8KG6a+fd/jc9+HH6Qi+ug5dg08jCYB3nf9rdizMahky+Mnx4ifEImDef/v3I5nbVZ7rjZ5Inof+ZJTT4BP3XwfvnTXQyuz5cEntnDWabvi2yKCHzz8dOnv6P2PncBf3v/czH2eCL59/2Omvud1KePo9wG4P3X7GICfLrHNPgAZRy8ilyOM+HHaC/4CDpx1SlV7yRI4cNYp+Knnnxbf/hsXnI2HntoyJTW86Pmn4uKXnp257+2v+nF88Y6ZsUaGl+57Dl7z4jOxa8PDW16xD2+64AX4wu0P4sRoUsues07bjb2nhI7pJWefhrcePAdPb9d7rjY589TdOPPUXZn7rnztT+K7Dz65UjsOnHUKfuYn98a3L33lfmwOy8trB846Bb/wl16Que/vvurH8dxnbzRm46r44xqPkXkHs4j8IoDXq+q7o9tvB3Chqr43tc3/AvBrqvq/o9tfAvB+Vb1l1vMePHhQjxw5UsNkQgjpLyJyi6oerPKYMsnYYwD2p26fA+CBGtsQQghpgTKO/mYAB0TkPBHZBHApgBty29wA4B1R9c2rADyxkz5PCCFkdczV6FV1IiJXArgRwADAdap6p4hcEf3/GgCHAbwBwFEAJwC8a3kmE0IIqUKZZCxU9TBCZ56+75rU3wrgPc2aRgghpAnYGUsIIR2Hjp4QQjoOHT0hhHQcOnpCCOk4cxumlvbCIk8BuLuVF6/GGQAeaduIEtDOZrFgpwUbAdrZNC9W1VOrPKBU1c2SuLtqd1cbiMgR2tkctLM5LNgI0M6mEZHKIwUo3RBCSMehoyeEkI7TpqO/tsXXrgLtbBba2RwWbARoZ9NUtrO1ZCwhhJDVQOmGEEI6Dh09IYR0nFYcvYi8T0TuEJE7ReQftWFDESJynYg8LCJ3pO77MRH5IxH5QfT7uTs9xyqYYecvRvszEJHWS8Rm2PhhEfletID874rI6S2a6GwqsvNDkY23isgfisgLdnqOVVBkZ+p/vyIiKiJntGFbzpai/fkvRORH0f68VUTe0KaNkU2F+1NE3isid0fH0r9py76UPUX781OpfXmviNw673lW7uhF5KUA/gHCtWgvAPBGETmwajtmcAjARbn7rgLwJVU9AOBL0e22OYRpO+8A8BYAX125NcUcwrSNfwTgpar6MgDfB/BPV21UAYcwbeeHVfVlqvpyAF8A8KurNqqAQ5i2EyKyH8DrANy3aoNmcAgFdgL4TVV9efRzuOD/q+YQcnaKyGsRrn/9MlX9iwD+bQt25TmEnJ2q+nfcvgTwWQCfm/ckbUT0LwFwk6qeUNUJgK8AeHMLdkyhql8F8Gju7ksAfDz6++MA/uYqbSqiyE5VvUtV16bTeIaNfxh95gBwE8KVyFplhp3pBVGfDaD1ioUZ300A+E0A78ca2AjsaOdaMcPOfwjgX6vqdrTNwys3LMdO+1NEBMBbAXxy3vO04ejvAPCzIvI8EXkWwgVL9s95TJuc5VbLin6f2bI9XeHvA/hi20bMQkT+lYjcD+BtWI+IfgoReROAH6nqbW3bUoIrIznsunWQP2fwIgCvFpE/E5GviMgr2zZoDq8G8JCq/mDehit39Kp6F4BfR3gZ/wcAbgMw2fFBpFOIyAcRfub/o21bZqGqH1TV/QhtvLJte/JEQdIHsaYnoRwfA/ATAF4O4EEAv9GqNbMZAngugFcB+CcAPh1FzevKZSgRzQMtJWNV9b+o6itU9WcRXpbMPSO1yEMicjYARL9bv5yzjIi8E8AbAbxNbTRx/E8Af6ttIwr4CQDnAbhNRO5FKIN9S0Se36pVBajqQ6rqq2oA4D8jzM+tI8cAfE5DvgkgQDjobO0QkSHCnNynymzfVtXNmdHvFyI0ttRZqSVuAPDO6O93Avj9Fm0xjYhcBOADAN6kqifatmcWueKANwH4Xlu2zEJVv6OqZ6rquap6LkIn9QpV/X8tmzaFC5Qi3oxQvl1Hfg/AzwGAiLwIwCbWd5rlXwPwPVU9VmprVV35D4CvAfguQtnm59uwYYZdn0R4aTlGeOD8MoDnIay2+UH0+8fW1M43R39vA3gIwI1raONRAPcDuDX6uWZN9+VnETqj2wF8HsC+dbQz9/97AZyxjnYC+G8AvhPtzxsAnL2mdm4C+O/RZ/8tAD+3jnZG9x8CcEXZ5+EIBEII6TjsjCWEkI5DR08IIR2Hjp4QQjoOHT0hhHQcOnpCCOk4dPSEENJx6OgJIaTj/H/bBiSALVBiHgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(trial.tracking.time, vavg)\n", "plt.xlim(xlim)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Event detection\n", "\n", "- Saccade events can be detected as peaks in the threshold velocity plot (above).\n", "- By taking into account the signature of the velocity, we can determine the direction of saccades." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "min. peak distance = 40 samples\n" ] }, { "data": { "text/plain": [ "array([-0.45002551, -0.91814375, -0.57781147, -0.42164273, 1.0871939 ,\n", " 0.42223875, 0.87461791])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def get_sign(x):\n", " return x / np.abs(x)\n", "\n", "# parameter configuration for peak detection\n", "# we care about the peaks only when the peak-to-peak distance exceeds 200 ms\n", "min_distance_seconds = 0.2\n", "min_distance_samples = int(round(min_distance_seconds / dt))\n", "\n", "print(f\"min. peak distance = {min_distance_samples} samples\")\n", "\n", "peaks = find_peaks(vavg, distance=min_distance_samples)[0]\n", "sign = get_sign(vleft_thresholded[peaks]) # leftward if negative, rightward otherwise\n", "events = vavg[peaks] * sign\n", "\n", "events # direction + amplitude of the saccades" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAATiklEQVR4nO3deZCkdX3H8fd3dpeFzXKEXRCQYzkUEEQRPAgKKCBRjIpIGQsFo4mooSSCZxIVAmooqxQJIqKA4FGKAiqIF6JA1ASWWKKcciyHHIHALi4ssLt888fTA+Mw0zPP9NNPP/3wflV1Tc/TPdOfmdn57DPf59dPR2YiSWqfkUEHkCT1hwUvSS1lwUtSS1nwktRSFrwktdTsOh9s4cKFuWjRojofUpKG3pVXXnlfZm5Q9uNqLfhFixaxePHiOh9SkoZeRNw6k49zRCNJLWXBS1JLWfCS1FIWvCS1lAUvSS1lwUtSS1nwktRSFrykCV1394MsXnL/oGOoB7U+0UnS8PjrEy4DYMm/7z/gJJop9+AlqaUseElqKQteklrKgpeklrLgJamlLHhJaikLXpJayoKXpJay4CWppSx4SWopC16SWsqCl6SWsuAlqaUseElqKQteklrKgpeklrLgJamlLHhJaikLXpJayoKX1GorVz/OGb+8hZWrHx90lNpZ8JJa7cxfLeGY86/hzF8tGXSU2lnwklrtwUdWAbD80VUDTlI/C16SWsqCl6SWsuAlqaUseElqKQteklrKgpeklrLgJamlLHhJaikLXpJayoKXpJaa3e3GiJjWfwCZ+fQ7i48kNVzXggdWATmNzzOrgiySpApNVfBbjrm+P/BG4FPArcAWwIeAc/oTTZLUi64Fn5m3jl6PiCOBXTNzaWfTDRGxGFgMfKFvCSVJM1LmIOu6wLxx2+Z1tkuSGmaqEc1YZwIXRcQJwO3AZsB7O9slSQ1TpuA/CNwIvAnYBLgLOAn4Uh9ySZJ6NK2Cj4hZwM+A/TLzlP5GkiRVYbrr3FdTrKiJ/saRJFWlzEHWY4BTImKLiJgVESOjl36FkyTNXJkZ/Jc7b986ZltQPBHKJzpJUsOUKfgtp76LJKkppl3wY5/0JElqvjJ78ETEa4E9gYWMOeCamYdUnEuS1KNpHyCNiI8DX+x8zEHA/wH7AUv7kkyS1JMyK2DeDuybme8DHuu8/RtgUT+CSZJ6U6bg18vM33euPxYRczLzcoqRjSSpYcrM4G+KiB0y82rg98C7I+IB4IH+RJMk9aJMwf8rsKBz/cPAN4D5wHuqDiVJ6l2ZZZIXjrl+ObBNXxJJkipRZhXNZyLidRGxfj8DSZKqUeYg60PAkcAdEfHbiDgxIg6MiA36lE2S1IMyI5qPAkTEXOAlFK/RejrFHN5z0UhSw0y74CNiPrA7xbLIvYDNgR8Dl/QlmSSpJ2VW0TwALAFOBN6emdf1JZEkqRJlCv4Y4GXAPwOviYhLKPbeL8/Mlf0IJ0mauTIz+OPgiZfvewFwAHAhxfx9fl/SSZJmrMwMfn2K+fuewMuBbYErcQYvSY1UZkRzB3A5cCnFcslfZeaKvqSSJPWsTMFvnZl3jd8YERtl5t0VZpIkVaDME52un2T7NVUEkSRVq0zBx1M2RKwDPF5dHElSVaYs+Ii4PSJuA9aKiNvGXoC7gO/2O6Skp68Tf/YHXnr8xYOOMZSmM4N/C8Xe+4XAW8dsT+CezJxsdCNJPfvMT28YdIShNWXBZ+YlABGxMDMf7n8kSVIVyqyiWRUR7wSez7gnNmXmIVWGkiT1rkzBnwXsBJwP3NOfOJKkqpQp+P2ALTNzaZ+ySJIqVGaZ5G3A3H4FkSRVq+sefES8Ysy7ZwHfi4jPMW5Ek5muYZKkhplqRHPaBNs+Oe79BLaqJo4kqSpdCz4zt6wriCSpWmVm8JKkIVLmfPC3U4xjxnuU4lTC5wJfyMxVFWWTJPWgzDLJEylOW3AicDvFi27/I/Bt4H7gKGAz4IMVZ5QkzUCZgn8bsG9m3jm6ISJ+CPwkM3eIiJ8DF2HBS1IjlJnBbwwsH7ftIWCTzvUbgPUqyCRJqkCZgj+fYh38PhGxXUTsA5zT2Q6wG7Ck4nySBEDmRIcA1U2Zgj8M+G/gi8BvgFOBK4B3dW6/Gdi/0nSS1GG/lzftGXxmPgJ8uHOZ6HZfl1VS39jv5U11qoI9MvPSzvVXTHY/T1Ugqd+KEc1TXjlUXUy1B38ysGPn+pnARGvcPVWBpL573F340qY6VcGOABExC9gAWDczH60jmCSNlQ5pSpvWQdbMXA1cDyzobxxJmpgHWcsr80SnrwMXdE4XfAdjjnk4g5ek5ilT8O/uvD163HZn8JL67nF34Usrs0zSUwdLGhj7vTxPFyxpKNjv5VnwkoaCI5ryLHhJQ8F+L8+ClzQcLPjSLHhJXTXlLI4+0ak8C17SUPBUBeVZ8JK6asgOfGP+khgmFrykoWC9l2fBS+qqKcXqDnx5FrykoeCIpjwLXlJXTSnWZqQYLha8pKHQkP9nhooFL6mrpvSqpyooz4KXNBSs9/IseEldNWXHuSnHAoaJBS9pKMy03+OJt9H1fm1U5hWderZsxUrO/+2ddT6kpuHuZY+w0bprDjpGT2b6NSRwT8Vffxu+n2P94Hd3Mntk8PuCF117Dwvnzy39cdfc9WDn7bKnXf9EnX/2zN34WbnxoSfU9niS1Aa3Hv+aKzNz17IfV+se/LOfsTbnHblHnQ+pKXzritv50mW3cOALNuXdew3nS+ueeunNnL34Dv5u90Uc/OLNS33soadfwR+XruDY1+3Ablsv6DnLv11wLZfecC9H7ftsXvXcjXr+fIMXNOPwZq85mvJ1zMyzjp/Zx9Va8HNnj7DNhmvX+ZCawuifvAvmrzG0P5sFna9hg7Xnlv4a1l6z+BXYZL21Kvn61583B4BnrLvm0H4/1R6DH6xJPRqdMvYybYyKj789/Q7nqYkseInqV1hE1f9jSDNgwUtQ+S73iP2uBrDgJaofqbgDryaw4DX0Rsu0l1J1pKI2suA19Co5yFpNlCf4rHo1gQUv4UhF7WTBS/RjFU2ln06aEQteovpCdkSjJrDg9bQ2enDVHW61kQUvQeUN74hGTWDBa+hVUaZPx3OFq/0seA29Xubdo6fLdgavNrLgJamlLHgJT1WgdrLgJao/VYEjGjWBBa+hV8lB1or63XPaqEkseA29KvaWq6rlfh20lWbCgpdwFY3ayYKXgKr24R3RqEkseInq9uDTXXc1iAWvoVfNM1mr5Y68msCCl3CZpNrJgpeobg/eGbyaxILX0KtkmaQzeLWQBS/h2STVTrO73RgRW03nk2TmzdXEkcrzmazSxLoWPHAjkBQjyrF/e45/f1bFuaSh5IhGTdJ1RJOZI5k5KzNHgL8HvglsB6zZefsN4B19Tyn1mTveaqOp9uDHOhZ4Vmau6Lz/h4g4DLgB+ErVwaTp6mWn+cnXZPWZrGqfMgdZR4BF47ZtgeMZDbHRkUriaEXtU2YP/rPAxRFxBnA7sBnwts52STiDV7NMu+Az89MR8TvgIGBn4C7g7Zn5o36FkyTN3LQKPiJmUczan2Ohq2maNPZ2Bq8mmdYMPjNXA6spVs9IjeJURJpYmRn8CcDZEfFJ4A7GrIP3iU4adj6TVW1UpuBP6rzdd9z2xJU0ktQ4ZQ6yet4aSRoilraGXpOOazYpizTtPfiImA28B9gTWMiYU2hn5h7VR5OGz7+8envWmjOL/XfaeNBRpFJ78J8FDgMuBXYBzgE2BC7uQy5pKC2YP5dPHPBc5s72sJQGr0zBvwF4VWZ+DljVeft64OX9CKZ6rLvWHADWmzdnwElmzmWS0sTKrKKZR3GKAoAVETEvM6+LiJ37kEs1OWjXzUjgjbtsOugokipWpuCvBV4IXA4sBo6OiAeBP/YjmOoxayR484s2H3SMnnhgU5pYmYI/guLZrABHAl8A1gbeWXUoSVLvyqyDv2LM9T8A+/QlkSSpEtM+yBoR50bEERHx/D7mkUrzIKs0sTKraH4AvAD4bkTcHxHfj4ijIuKFfcom1cYX/FAblRnRnAacBhARW1DM3j8GzMdz0WiAPMgqTazMM1m3o3gW657AS4G7gS8Cl/QnmiSpF2VW0VwD3AR8CviHzHyoP5Gk+nm6YLVRmRn8IRSnJXg/cGVEnBoRB0fEZv2JJknqRZkZ/NeArwFExDOA9wIn4wxekhqpzAx+Z2Avihn8y4AVwAU4g9eAuUxSmliZEc15wPOA7wMvysxNM/PgzDy1P9Gk/hs9B89G6/pyw2qfMgdZ/zYz/2v8xoh4UWZeXmEmqZRelkm+46Vb8ra/WsTsWb72jdqnzL/qn0yy/UdVBJEGISIsd7XWlHvwETFC8epNEREBf7aebGtgVZ+ySZJ6MJ0RzSogKYp9fJk/Dnyi6lBSGR5klSbWteAj4nnAlhTlfgkw9rVXE7g3M1f0L54kaaam2oO/LDPXAYiIxzLz1hoySaV4LhppYlMV/NKIeA3FaQo2iojRvfk/k5k39yOcJGnmpir4I4ATgC0oVtzcNMF9Ep/JKkmN03V9WGael5nbZOYc4OHMHJngYrlLUgOVWQC8AIplkxGxcZ/ySJIqUqbg50XEN4BHgBsBIuK1EXFcX5JJknpSpuBPAZZRzOMf62z7NfCmqkNJknpX5lw0ewObZObKiEiAzLw3IjbsTzRJUi/K7MEvAxaO3RARmwN3VZpIklSJMgX/ZeCciHg5MBIRuwFnUrwuqySpYcqMaI6nOMD6eWAOcDpwSmZ+rh/BJEm9mepcNK8Yt+kq4PDx98nMi6sOJknqzVR78Kd1uW30DJMJbFVZIklSJboWfGZuWVcQSVK1fCkbSWopC15Dzxf8kCZmwUtSS1nwGnq+4Ic0MQteklrKgpeklrLgJamlLHhJaikLXpJayoKXpJay4CWppSx4SWopC16SWsqC19CbM2uk89antEpjlXlFJ6mR3rXnVjy6cjWH7LZo0FGkRrHgNfTmrTGbj7x6+0HHkBrHEY0ktZQFL0ktZcFLUktZ8JLUUha8JLWUBS9JLWXBS1JLWfCS1FKRmfU9WMSfgOtre8CZWwjcN+gQ02DOag1DzmHICOas2raZuXbZD6r7mazXZ+auNT9maRGx2JzVMWd1hiEjmLNqEbF4Jh/niEaSWsqCl6SWqrvgT6358WbKnNUyZ3WGISOYs2ozylnrQVZJUn0c0UhSS1nwktRSFrwktVQtBR8R20fExRGxLCJujIgD6njcqUTE4RGxOCIejYivjLtt74i4LiIejoifR8QWA4o5ac6IWCMivhMRSyIiI2KvBmZ8SUT8NCLuj4h7I+LbEbFxA3M+p7P9gc7looh4TtNyjrvPxzs/931qjjc2w2Tfz0WdbMvHXD7atJyd2+ZFxMkRcV+noy4dUMxu38+Dx30vH+58f3fp9vn6XvARMRv4HnABsD7wTuBrEfHsfj/2NNwJHAecPnZjRCwEzgU+SpF5MfCt2tM9acKcHf8JvAW4u9ZETzVZxr+kWAGwCNgC+BNwRq3J/txkOe8E3kjx814IfB/4Zr3RnpJnsp85EbE1Rd676gw1ga45gfUyc37ncmyNucbrlvNUip/79p2376sx13gT5szMr4/5Ps4H3gPcDPxPt09WxzNZtwM2AT6bxZKdiyPil8BbKQp0YDLzXICI2BXYdMxNbwCuzsxvd24/GrgvIrbLzOuakjMzHwNO6Ny2uu5cY3XJ+MOx94uIk4BL6k33pC45lwJLO7cFsBrYpv6ET+SZ7N/mqJOADwEn15lrvGnkbITJckbEtsBrgU0z88HO5ivrT1go8f08FDgrp1gGWceIJibZtmMNjz1TOwC/HX0nMx8CbupsV2/2AK4edIjJRMRS4BHgP4BPDjbNxCLiIOCxzLxw0Fmm4daIuCMizuj8Zdw0LwZuBY7pjGh+FxEHDjpUN51x8R7AWVPdt46Cvw74X+ADETEnIl4J7AnMq+GxZ2o+sGzctmVA6ZP96EkRsRPwMeADg84ymcxcD1gXOBz4zWDTPFVEzKf4j+efBhxlKvcBL6QYy+1C8bvz9YEmmtimFDubyygmDYcDZ0bE9gNN1d0hwGWZectUd+x7wWfmSuD1wP4Uc+KjgLOBO/r92D1YDqwzbts6FPNjzUBEbAP8EDgiMy8bdJ5uOn+xnQKcFREbDjrPOMcAX53OL/cgZebyzFycmasy8x6K4nxlRIz/vRq0FcBK4LjMfCwzLwF+DrxysLG6OgQ4czp3rGUVTWZelZl7ZuaCzNwP2Aq4vI7HnqGrgeeNvhMRfwFsTYNHC03W+ZPyIuDYzPzqoPNM0wjFX5nPHHSQcfYG3hsRd0fE3cBmwNkR8aEB55rK6Kx4opHtIF016ABlRMTuFH9pfGc6969rmeROEbFmZznS+4GNga/U8djdRMTsiFgTmAXM6mScDZwH7BgRB3Zu/xhw1SAOsE6Rk4iY27kNYI3ObbX/Ek2WMSKeCVwMfD4zT6k713hdcu4bETtHxKzOXuZngAeAa5uUk6LgdwSe37ncCRwGfL5JOSPixRGxbUSMRMQC4ETgF5k5fvQ50JzApcBtwEc699kd2Av4ccNyjjoUOCczpzdNyMy+X4BPU/yyLKf4M32bOh53GrmOptizGHs5unPbPhTHD1YAvwAWNTTnkgluqz3rZBmBj3euLx97adr3Ejio8/NeDtwLXAjs1LScE9xvCbBP03ICbwZuAR6iWMp5FrBR03J2btsB+HUn6zXAAQ3NuSbFSq+9p/v5PNmYJLWUpyqQpJay4CWppSx4SWopC16SWsqCl6SWsuAlqaUseElqKQteklrq/wEnMUgwDff7jAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# turning events into samples\n", "saccades = np.zeros(vleft.size, dtype=np.float64)\n", "saccades[peaks] = events\n", "\n", "plt.plot(trial.tracking.time, saccades)\n", "plt.xlim(xlim)\n", "plt.yticks((-0.5, 0.5), (\"leftward\", \"rightward\"), rotation=90, va='center')\n", "plt.tick_params(labelsize=12, left=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.10.5" } }, "nbformat": 4, "nbformat_minor": 4 }