{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "1f32dffd", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "from collections import namedtuple\n", "import math\n", "import numpy as np\n", "import pandas as pd\n", "import scipy.stats as sstats\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 as reader\n", "import saccades" ] }, { "cell_type": "code", "execution_count": 2, "id": "78d28ed6", "metadata": {}, "outputs": [], "source": [ "dataset_root = \"../01_data/04_formatted\"\n", "figdir = Path(\"../05_figures/saccades\")\n", "saved = False" ] }, { "cell_type": "markdown", "id": "df086eed", "metadata": {}, "source": [ "## Preparation of trials\n", "\n", "- scaling whisker positions\n", "- computing whisker asymmetry\n", "- detect saccades" ] }, { "cell_type": "code", "execution_count": 3, "id": "f02fb5fc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "## Number of trials/epochs used for analysis\n", "\n", "- all annotated trials: 113 trials out of 21 sessions from 4 subjects (876 epochs)\n", "- trials with eye positions: 91 trials out of 17 sessions from 4 subjects (728 epochs)\n" ] } ], "source": [ "def get_info(trials):\n", " sessions = set((trial.subject, trial.session) for trial in trials)\n", " subjects = set(trial.subject for trial in trials)\n", " epochs = sum(trial.states.shape[0] for trial in trials)\n", " return f\"{len(trials)} trials out of {len(sessions)} sessions from {len(subjects)} subjects ({epochs} epochs)\"\n", "\n", "alltrials = reader.load_trials(dataset_root)\n", "trials = [trial for trial in alltrials if trial.has_eyedata()]\n", "\n", "print(\"## Number of trials/epochs used for analysis\\n\")\n", "print(f\"- all annotated trials: {get_info(alltrials)}\")\n", "print(f\"- trials with eye positions: {get_info(trials)}\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "1a1b2c22", "metadata": {}, "outputs": [], "source": [ "def normalize(vec, subtract_min=True):\n", " m = np.nanmin(vec)\n", " M = np.nanmax(vec)\n", " den = vec - m if subtract_min == True else vec\n", " return den / (M - m)" ] }, { "cell_type": "code", "execution_count": 5, "id": "0142d165", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
timeleft_whisker_angle_degleft_whisker_radius_pxleft_pupil_normalized_positionleft_pupil_normalized_diameterright_whisker_angle_degright_whisker_radius_pxright_pupil_normalized_positionright_pupil_normalized_diameterleft_whisker_normalizedright_whisker_normalizedleftward_whisker_asymmetryrightward_whisker_asymmetrysaccades
00.0003.513870271.1515560.0486740.08670213.644373243.7318270.0569160.0978010.0260680.0933150.067246-0.0672460.0
10.0054.533231272.7853000.0486740.08670214.096692247.6674800.0569160.0978010.0336310.0964080.062778-0.0627780.0
20.0105.875672273.4176050.0486740.08670215.003054248.1331200.0569160.0978010.0435900.1026070.059017-0.0590170.0
30.0156.574637271.1895450.0486740.08670215.494432246.5431580.0569160.0978010.0487750.1059680.057192-0.0571920.0
40.0206.493624271.0029090.0486740.08670215.364230246.8022740.0569160.0978010.0481740.1050770.056903-0.0569030.0
\n", "
" ], "text/plain": [ " time left_whisker_angle_deg left_whisker_radius_px \\\n", "0 0.000 3.513870 271.151556 \n", "1 0.005 4.533231 272.785300 \n", "2 0.010 5.875672 273.417605 \n", "3 0.015 6.574637 271.189545 \n", "4 0.020 6.493624 271.002909 \n", "\n", " left_pupil_normalized_position left_pupil_normalized_diameter \\\n", "0 0.048674 0.086702 \n", "1 0.048674 0.086702 \n", "2 0.048674 0.086702 \n", "3 0.048674 0.086702 \n", "4 0.048674 0.086702 \n", "\n", " right_whisker_angle_deg right_whisker_radius_px \\\n", "0 13.644373 243.731827 \n", "1 14.096692 247.667480 \n", "2 15.003054 248.133120 \n", "3 15.494432 246.543158 \n", "4 15.364230 246.802274 \n", "\n", " right_pupil_normalized_position right_pupil_normalized_diameter \\\n", "0 0.056916 0.097801 \n", "1 0.056916 0.097801 \n", "2 0.056916 0.097801 \n", "3 0.056916 0.097801 \n", "4 0.056916 0.097801 \n", "\n", " left_whisker_normalized right_whisker_normalized \\\n", "0 0.026068 0.093315 \n", "1 0.033631 0.096408 \n", "2 0.043590 0.102607 \n", "3 0.048775 0.105968 \n", "4 0.048174 0.105077 \n", "\n", " leftward_whisker_asymmetry rightward_whisker_asymmetry saccades \n", "0 0.067246 -0.067246 0.0 \n", "1 0.062778 -0.062778 0.0 \n", "2 0.059017 -0.059017 0.0 \n", "3 0.057192 -0.057192 0.0 \n", "4 0.056903 -0.056903 0.0 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for trial in trials:\n", " for side in (\"left\", \"right\"):\n", " trial.tracking[f\"{side}_whisker_normalized\"] = normalize(trial.tracking[f\"{side}_whisker_angle_deg\"], subtract_min=False)\n", " trial.tracking['leftward_whisker_asymmetry'] = trial.tracking.right_whisker_normalized - trial.tracking.left_whisker_normalized\n", " trial.tracking['rightward_whisker_asymmetry'] = trial.tracking.left_whisker_normalized - trial.tracking.right_whisker_normalized\n", " trial.tracking[\"saccades\"] = saccades.detect(trial.tracking.time,\n", " trial.tracking.left_pupil_normalized_position,\n", " trial.tracking.right_pupil_normalized_position)\n", " \n", "trials[0].tracking.head()" ] }, { "cell_type": "markdown", "id": "eeb5221b", "metadata": {}, "source": [ "## Extracting whisker asymmetry around saccades" ] }, { "cell_type": "code", "execution_count": 6, "id": "0641f98b", "metadata": {}, "outputs": [], "source": [ "def extract_traces_around_saccades(trials,\n", " pattern,\n", " trace,\n", " detect='left',\n", " threshold=1e-3,\n", " window_radius=500,\n", " baseline_radius=10,\n", " shuffle=None):\n", " \"\"\"returns a list of 'epochs', with each epoch being dict consisting of keys:\n", " \n", " session: tuple of (subject, date)\n", " baseline: a float value (peri-saccade median)\n", " trace: a np.ndarray of size (window_radius * 2) + 1 (peri-saccade trace)\n", " \"\"\"\n", " pattern = str(pattern)\n", " trace = str(trace)\n", " detect = str(detect).lower()\n", " threshold = abs(threshold)\n", " window_radius = int(window_radius)\n", " baseline_radius = int(baseline_radius)\n", " \n", " if shuffle is not None:\n", " np.random.seed(shuffle)\n", " \n", " if detect == 'left':\n", " matcher = lambda val: val > threshold\n", " elif detect == 'right':\n", " matcher = lambda val: val < -threshold\n", " elif detect == 'both':\n", " matcher = lambda val: np.abs(val) > threshold\n", " else:\n", " raise ValueError(detect)\n", "\n", " epochs = []\n", "\n", " for i, trial in enumerate(trials):\n", " # detect epochs that match the given pattern\n", " valid = np.zeros(trial.tracking.shape[0], dtype=bool)\n", " for start, stop in trial.get_timeranges(pattern):\n", " valid[start-1:stop] = True\n", " if np.count_nonzero(valid) == 0:\n", " continue\n", " \n", " # detect saccade events\n", " sacc = matcher(np.array(trial.tracking.saccades.values, copy=False))\n", " if shuffle is not None:\n", " extract = sacc[valid]\n", " np.random.shuffle(extract)\n", " sacc[valid] = extract\n", " evts = np.where(np.logical_and(valid, sacc))[0]\n", "\n", " # extract time windows around events\n", " for evt in evts:\n", " if evt < window_radius:\n", " continue\n", " window = trial.tracking[trace][(evt-window_radius):(evt+window_radius+1)]\n", " base = np.nanmedian(trial.tracking[trace][(evt-baseline_radius):(evt+baseline_radius+1)])\n", " epochs.append(dict(session=(trial.subject, trial.session), \n", " trial_id=i,\n", " baseline=base,\n", " trace=window))\n", " return epochs" ] }, { "cell_type": "code", "execution_count": 7, "id": "c4b83a9f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "leftward: 120 traces / 63 trials / 16 sessions / 4 animals\n", "rightward: 94 traces / 51 trials / 16 sessions / 4 animals\n" ] } ], "source": [ "pat = ''\n", "attr = '{toward}ward_whisker_asymmetry'\n", "radius = 120 # in samples -- 600 ms at 200 Hz\n", "\n", "epochs = {}\n", "for toward, side in (('left', 'right'), ('right', 'left')):\n", " _epochs = {}\n", " \n", " # true traces\n", " _epochs['data'] = extract_traces_around_saccades(trials,\n", " pattern=pat,\n", " trace=attr.format(toward=toward),\n", " detect=toward,\n", " window_radius=radius)\n", " \n", " # traces using shuffled saccade timings\n", " _epochs['rand'] = extract_traces_around_saccades(trials,\n", " pattern=pat,\n", " trace=attr.format(toward=toward),\n", " detect=toward,\n", " window_radius=radius,\n", " shuffle=539167)\n", " _epochs['N'] = len(_epochs['data'])\n", " n_traces = _epochs['N']\n", " n_trials = len(set(epoch['trial_id'] for epoch in _epochs['data']))\n", " n_sessions = len(set(epoch['session'] for epoch in _epochs['data']))\n", " n_animals = len(set(epoch['session'][0] for epoch in _epochs['data']))\n", " epochs[toward] = _epochs\n", " print(f\"{toward}ward: {n_traces} traces / {n_trials} trials / {n_sessions} sessions / {n_animals} animals\")\n", "\n", "# data merged\n", "dmerged = np.stack([epoch['trace'] for epoch in epochs['left']['data']] + \\\n", " [epoch['trace'] for epoch in epochs['right']['data']], axis=-1)\n", "\n", "# shuffled merged\n", "rmerged = np.stack([epoch['trace'] for epoch in epochs['left']['rand']] + \\\n", " [epoch['trace'] for epoch in epochs['right']['rand']], axis=-1)" ] }, { "cell_type": "code", "execution_count": 8, "id": "44aae4d5", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALMAAAECCAYAAAC19eaXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABAxElEQVR4nO2deZgcVdW439PrrNlDQhYStrBGQEAIS4isYRFQ1gQ1qBhUQEEQFZFNAf1k+z5REUFRCCGAbAYQBISAkADRHwRkS1iyh6yz91JV5/fHre6pnultprtnS73P089U36q699b06dvnnnvuOaKq+PgMBAK93QEfn3LhC7PPgMEXZp8Bgy/MPgMGX5h9Bgy+MPsMGPqkMIvIRBFRETnJU7a0yHt/JiKfiMgznrJ9RORfIrJARJ4TkR3cchGRX4vIiyIyX0SGZanveREZl6e90SLyioj8U0Q+JyJTu/a0xSEi40Tk+QLXnCgii9znOavDuUkikhSRQ4po6xYRGZnn/NkicnmW8pNFZLtC9XcHEfmyiFyV75o+Kcwu7wI/FhHp4n2/BT7foWwNMF1VpwI3AFe75ccANap6KHA/cGk3+vl54GlV/TywO1AWYRaRLn027vW/Ao4EDgfOF5Ehnkt+CrxQRD1BVb1QVdd3pX2Xk4GyCLOIBLt6T6gcDVeIVcD7wEnAI8XepKprRGRih7K1nrcJwHKPpwHz3eO/Ad/KV7eIXA8cBESAa4H/AFcCUREZC3wOqBeRI4FZwP+q6skiciMwRlVniMhvgDlAFLgC8xlsAs5Q1Zj7C3Q/MMX9ZZrrXvtmgUcfAaxX1Sa3r++7/XlaRD4HrAXsHM91FTARGAbMFZFzgS9jBoG7gbHAK8CpqrqTe9tkEfkrsAvwbWAjMB3Yy32G+cBQVb1FRF4H7lDV20TkVeBg4ALgBKAeeEJVr3Q/twcwA1lSRG4A7gLWu/V/mPc/oKp97uX+Y58BxgELAQGWuuemAM9neR3e8f4s9dYCi4Dd3fe3A9PcYwHezXLP824/pgO3uWU1wBvuPWcDl7vl6WP3/WIgCDyB+bII8BoQBmo91/0S+Kp7/DEwxT2+CPixe3wW8Hye/5m4QjAWGAQsA850zz0GDHcF45As914F/D7LM3/J88wHAx97nvMR9/gg4EH3OF2/e//fgKFu+3OBUcDfU5+F+zcAvIwZ0SdiBHeQe+5Rz//iD8BV+eSmL4/MqOpKEVmM+flKlb2CGVG7hIiEgXnA9ar6X7d4EzDEPR4MbM5TxWTgMI/eGsUISD4WA0dhRpUV7vF6VU26evzP3XpGAY3uPTbmCwwwCXjQPV4EfDNXQ6qqIjIbM5K2AEuA1SJyPPC6qm4soLG9nKVsZ8yXL9W+1/dhsft3OVn+D+5nty1wBEaov+Ae/9O95BQROcetcwdgPObX+C1VTf0vdgZe9bSfc+4CfVtnTnE98MPUGxGZ4k7KOr4Oz1WBq0/egxlNHvGcegE4zj0+jvw65dsY3Xiaqk4DPqOqGzpckyBTdXsOo58/7x7/zD0G+Alwpaoehhm5UpKm2u4w8wGwn3u8v+d56kVkaMcOquoCVT0cmAnUYQRgb2CaiPwd82W6QUQmZHm+bCrI0g7te78NXsFOlXd8/teBH2Ce/U3ge7QL888wc5bPAx956vD2o2P7eenTIzOkv+GvY37mC47MInI+cCawm2vROBfYBzgeGCUiXwaWqOoFwFPACSLyImZk/GqefjyR+iJhPsiVwFc6XPYvzMRrT+B8zId4r9uftcBewHnutfcBd4rIe0AD7SOzlz8A94vIUcBbnvIzMbrmTR2e/X8wH7qFUU/iGN3+Wvf8XRjd9ZNcz9mBR4DTROQFzAgdL3D9fOAaEXlHVc8FngWOVdVlIvIsRk9OjegPYf5f7wLNOeq7DPijiGwEOg4cnZD2QcCnvyAi/wtcq6qf9kBbYVctOhjzBTmh0m12F1+YffLiWixGYHT7c1X1jV7uUk58YfYZMPSHCaCPT1H4wuwzYPCF2WfA0G+Fefr06YoxkfmvMr6uuuqq3my/JPqtMG/YUNDs6NMNZs+e3dtd6Db9Vph9KsPtt9/e213oNr4w+2QwefLk3u5Ct/GF2SeDurq63u5Ct/GF2SeDV155pdfafl6e7+pGjAx8YfbJYMaMGT3Wlt1mYzVZ3qJIKfX5wuyTwfz58wtfVCacVgen1fEWVZdSny/MPhnE44W8PLuP2orVaEZidRQn4aB2xo6XaCn193l/Zp+eZebMmRWr24k7RoBVsVvt9DKJWgpmgO7yJlYv/sjsk8Fdd91VsbqduINailqaoV5Ymy2cmJPnzuLwR2afDPbdd9+K1Ku24sQdcMBpM0KdPmcptpV143iX8IXZp+KorSQ+TaTVCruldMHNhq9m+GSwePHiwhd1EbU0042oQvtBfGH2yeDss88ua32qito9s5vJF2afDO69995u36uOoo4RXLXNsd1s+8Ls0ztEo8WbejvuH9WkmpcqVpNlTHExpyhhdhyHGLEu99eLL8w+GZxwQnGRBOxWG2tzxlK0MbklXPNbQtPC7bQVNrvZts17vFfSoklRwiwi40XkBDes6AkiMr6URn36LnPnzs17PrVqZzfbOIlMIfXakdVyTXFQ1IRPVXmbt2u722/IY5pzY7Od6752wIRKasJE0tlJRD4CbgNuV9VEKZ3w6TtMmTIl73m7xUYCkmEnBiOMTsxBwoImXb05Wbyu7DgObbSVZCrOd/MbmPBS5wKLVDVtHHRj534OE5nyP8AepXTCp+/Q3GwiZamtSDDTIzO5JWlW6jzF6qgRbrt94tdxxC4GVcXBqZgL6DRVPV9VX/YKstuwraqvqOr5dCMip0/fZcmSJQBYDR30YXV1XyXlR2HKndRKiFvggCa6br1wtPTl7JzCnC+OmYiMSEW01+5FWPfpo8yePdt4tMUcrGaL5MYkTtLpvPCRwhXitFAXSVNzE8lkMm0RKUdkra6mGpgqIp9gwpOuF5HTSu6BT59BVfn973+f1nXtRhsn7mBtsnLqv171olgsyyKeiJuX63LqOBUcmQFEpOPs8kpgqqqOAQ4Dbim5Bz59BrWV4UOGt1shPOUdy7zngAzVIxuWZdHW1gZAMpnEsiwsy8J2zNDeEyPzAhE5xfM+CYwWkRAmirlvxRhAaFI59KBDswpuPmF2rMILI7Ztk7SS5tixcRyHpGWEOnW+VAoJ8+HAESLyhIjsiImC/n+YwNg3At8ouQc+fQa1lL/e/9fsKkWOkVcTitPiFDTD2Y6dFtjU39QI7TgO8UTpO1zy2vVUtQH4jput6C+YpDlT3YjsPv0cdYwTUCBsxjS1lalTupb5TS3FEaeT3dmLbdsZKoVXP44n4mzasqlnJoCu1eJDjI68AXhFRI4tuWWfXic9uUs5B1nKmnVrulxPoVE5mUzS2taaFuKOk71EojzaaqEJ4OkYAV6CSen1FiaRzVdF5KF8mUt9+j5qmZE5Para8MGHH5S9nZSOnFIvymFTzkahkflm4POqui0mwc11qrpWVWcAv8NkSfLpp6QmbVajlZ7EfeOs8k+DUiOxqhp1owyTvWwUEuY4xoIBxmSe9tFT1X8AB1akVz4VRx1F4656kVCSn5qP+c45d5a9La/wxhPxstiUs1HIseObwDwRqQE+pUM6Xt/BqP/ixLOb00ZvM7rsbXmFubklV5a00ilkzXgW+EzFWvfpNXLZhffdq2u7sx3HwbIsIpH2yFqWZREKtYuWZbf7eVRKxYA8aoaI7FVMBcVe59PHyCFTj//j8aKrSCaTxOPGtJZCVUkmkxnvKynAXvKNzL8RkUZMLuYXVHV16oSbE/kwTEbTOqBrxkmfXifXyHzk1COLrmNL4xYi4QiO4+A4DoFAAMdxSCQTVFebsHE9JciQ32vuEOC3GJ/lpSLSJCKrRaQJk9P5TOBWVfUFuR+SS5i7YpqzLIu2WFv6GIza0dLaktaNW1pbSuxp8RTSmecD891dJzsDQ4DNwAeqauW716fvoqo5V+w+WVlcWm3HcbwBD7FsC6vNIhQ0ItXc0kxNdU1a2HuCorapqGoS+G+F++LTQ2gyh28yZLUzO45DW6yN2pp2J8rUSJwimUySSCbS1ziOQ2NTY8XMcNnwd2dvheTb1pTNztzQ2EBjU2NGWUddOJFMYFlWxtJ0a1triT3tGn6sua2QfNuaJoybkD5Ojb4pVUFVcTcYZZjbgLQFI5HsvaUHf2TeCsnnGLTzDjunj1Mumim8KkNHNaNQeU9QbNyME12HfJ8BQL79es8seCZ9nLSSORc8Oo7M3aV1QSvWGqtjOohuUayA/gy4U0TmAXer6qKSW/bpFVRzT/7A2JljsRiRSATLsggG24PZp0bmRCKRsTDSXeJL4qy/qHz7oYsamVV1L+BIoA34q4i8JyKXi8jEsvXEp2fIMwCqKv969V9p3+N4Ik4s3h7/LeW6WS7bccsT5bVBF60zq+obqvoDYDxwHnAasExEFojIWSLi69/9gHwqRjKZ5NMNn6b9jzsuRcdiMRKJRFm2OAHE344j0ZLivmTQJT3Y3Qf4ZfflAFcAy4HzgVOAL5WtZz6VIc/InLSSzPjSjPTydEdi8RjxRLwsW5w0oSTeTzDozEG0vdJGcmnpakuxE8DzRWQhsAgYBXxFVXdR1WtV9W7gCODoknvjU3HyjcyWZTH3obnpndNZ7y+DIAMkPkhAEiJ7RBgzbwx1p5We5rjYkXk6Zjf2o9l8mFW1VUT8Ubk/kGVkTrlsJq0k22+3Papatn15uUi8beqP7mmi2AaHlZQ1DShuQ2sQ4xn3WD5nfFV9uuTe+FQM76bVFKloQrZtp3dQjxo5Cqj84kf8v3ECQwMERxsh7hFhdoMmbk9G7Eef/obVYHVyMGpsNkvUtmMTi8dwHIeFixcC5VMnvKiltL3ShlpK4u0E0T2i6RXF4PDShblYNeNq4DYRuRJYicdSqVqhrbY+ZcWJOdhN7flFUk70KYtFygQ3/fDpJbeV+DBBeGIYCWSOf1tu20Ljn9p9PGqOqkkfhyaUviZXrDntDowj/oeYkFxJwKJ9s2vRiMjOIhITkXs8ZUeIyLsi0ioi/xSRCfnq8Ckex3JwLBOK1m620yOzNyyW4zjpRZBFi7u/HmatsWhb1Maa09aw5qw1bLrRBHdxmh2aH2um6d4mAAKDA9QeW0v9afXpe8Pbhbvdbopivw7bl9xSO78BXku9EZERwEPAOcDfMKuN8/B3fpeFjsvEbW1txBPx9J49y7IyLBebGzZ3u621s9Zib3RDb72fJPl+kubHmglvHyaxxOjg287dlvBOnUdtCZWuxRY7Mp+mqp90fGFsy0UjImcCW4BnPcVfAt5W1QdUNQZcBewlIrt2pW6f7DgxxyRdd7FsszskZa2IJ+IZS9MzvjQjfayqbLh8A00PNqXLkh8lafl7SyfPO1VNC3JGebMaQQ7C8CuGE5kU6STI5aJYYb4iR/nlxTYkIoOAa4CLO5zaA5NyAgBVbQGW4aeWKBnHcmPAeQbn1Mpeyq2zta01Y7I396H2BD3xN+O0PNnCpus3maxRTQ5rZ69lw082sOGqDaitrPvOOhrvbcTekCnI4xeMZ9xz7QGvRt8xmrqTSrcl5yOvmiEih7uHQRH5PJkWjR0wCXuK5WfAnaq6IjWDdakDOnqbNGASAXXsz2xgNsB2223Xhaa3PrxBXrx4owt5/6bYZadd0sexRe1+GYn3EjiNDs4mh9DYEK1PtbL8qeXp65wmU29w2yChbUIEas04OfTioVhrLSK7R6g0hXTm1LaDKuCPnnIF1gEXFNOIiOyNcVTaJ8vpZmBQh7JBZPmiqOrtwO0A++23X8+k/ewH2C02wdpM05bT5mA3d/7ZL7SNqba6fWtU8uN29SP5YRJ7vamv/sx6Nt+YqVs33N5AoD7A2IfHIuH2wWrQzI4fbeUotKF1ewAR+YuqfrWEdqYBE4Hl7qhchxntd8ekX5uVutCN1r8j8HYJ7W1VWE0WEpH20LSuPTnbDmwtkJTv30v+zdSDzIZ762OLqgOqiC2KsfHqjelrovu0557cbtF2bPzZRlrmt1B9WHWGIPc0xbqAflVEwiJyqIicAUbosqSJyMXtGAHd233dBjwOHAM8DOwpIqeISBVGP39TVd/t0pNspdgtdqcMT9YmK+fu63wjs73R5ti9jmXLH7ZgrbZILk8S3iFMeJd2s1n9zHrCE837QbMGISEhvL15X31wdTkeqdsUZZoTkcmYiJ9xTPqHeZggMLOAMwrdr6qtQHp3o4g0A7FUpio31cStwD0YZ6Yzu/QUAxjvvrtsODEn/TdQHQDpnLIhmUxi2zZVVVU5hdlpclh59Eoe4zFO53RjsWhTwhPCDP3eUJYfaPTjYRcPA2D8C+ORGtOvQWcNIjwhTPW0fiDMmPC1V6jq3SKSUpZeAP7QnUZV9aoO758BfFNcFuwm2wgppNUIL6kR2Ik7WI0WwbpM3Tm1j6+xuZFIJJJTmONLjJ9Gwk1TY31sFlXCE8NIWBg8ezASaf9SBera+yJhoebzNfQ2xQrzHphRE9ylbFVtEZHe/SoOcFTdHNVtJitqZJtIp/NevTiVs9pLapnatu28AVmsNUZ4jyQzPFdoohGRIecO6fZz9BTF2pk/BjLCQ7p5TpaWu0M+bjbUuJNWIVLR7a1GCzvWbqHolAHKztx57TgOiUQivUCSV5hXWxCA5w94npE3jqT2uFrqZ9YTHFG6A1BPUezI/FPgcRG5DYiIyI8xsZq/WbGebcU4cQeNK04yU1jtZhuJCMGqIOoodkOm6U1tTX8BwN14aiUzNqLmwlptERoTYr8z96Nmag0103pfbegqxVoz5gPHAiMxuvIE4Eu+D3Nl0IRit9nZg7W4smo1WNlNb56ROeXWWQinySG2KEZkUuUXNipJ0X53qvpv4DsV7IuPi5Nwcufdc5TWhlaCbdl//lM7q+tq6/JufdK4EqgyY1nri604DQ71Z9WzZNESjph6RFmeo6cp1jQXAmZgVvAyFthVdXYF+rXVorbmDZ+FA/GWOKFEKO35llqSFhHi8TiWbVFXW5czNnLTvU1svskYpaoPqwYbpEaI7hnl1G1OLe8D9SDFTgDvAX6EGS/WdXj5lBG7pXBw7mRbMmNbUywWo7Gp0WRzsk02J1UtaIYDaHuhjbaX2ohOjiIh4dG/P1r6Q/QSXdnQOl5Vu+JY5NMNUu6a3sWS1Aibii5ktVnYYoO7/hqLx7Bsiy0NW9LunPliviWXJ6k6qIpt/ncbVp+2Gutji+pDjJU1Go7mvK+vU6ww/xcYRte85HyKJCW46rS7a8ZiMQKBAPFEPD3SRqNRwqEwiWQCVaWltYXamloSyUQ6UU6KXOFkVRVruUXVPlVIQNjmhm3YdNMmak8w34zDDz086339gWKF+cvAHSLyNB1UC1X9S9l7tZXhtJqlaK/d2LItmhrM2BEMBlHVjFBZYKLTRyPRrLpxrhBazhYHbVNCY81HH94+zKhfj0qff+ypx7jw3AtLfaS8PP3Pp3lv6Xtc8M2inC6LplhhPhs4FBiKiTeXQjEJ4n26iRM3O0HsVju9TTiZTGbs/sg1kbNtu8tx31JO9MGR2a0hn5382S7VVyyr164mEokwYtgIvnWJSSf5rbO/RThc+t6/FMUK8/eAfVT1nbK17AN4Npl6DBhbGrcUHWWzy8Lsbm3KtbLX0laZhDqHHH8I4VCYt156K1327tJ3WfT6IvbcbU/232f/ktsoVpjXYWLK+ZQRJ+50WpL27pSuBGlhzhF05b2l73HsEceWtc2U6TBpJXn1P6+myy+54pJ0dquzzzy75HaKNc3dDMwRkQNFZAfvq+QebMVkM8OlogxVrM2UMOcIuuLd0FouGhob0sf3P3I/NdVmqfyDDz9g+LDhANx1310lt1OsMP8GOBF4GeNclHoVnzTOBzCbTFNkW46udKoxZ6ODRAWpze4j7d3QWgpbGrbwwssvEE/EWbl6Zbr8qX8+xSlfOCUt0KefdDp/vvXPZWmzWN+MQI5X/3Gp6iM4LR5hzrIbpJIx3tQxYWRDE0I5Hf6HDh5alrb+cPcf+NoFX+O7P/oul/38snR5IpHg9JNOT5sOd5u0G4dOOZRfX//rktv0A4T3IIn1CexW2+ycdjqnY0hFF6oUm2/cTOzVGNG9ci+MHLDvAWVpKzUa/+OFf/DWu+2TvtqaWvbYtT2KxA4TjKZ6/NHHl9xmsfGZtxORO0Xk3yLyvvdVcg+2EtTWdDJJJ+FkVTEqnamp+TGTArjuhNzxK/7+3N/L0tb6Ddlzldxw9Q0Z7yduN7Es7UHx1owHgHcxm017Ln/sAMIbsFATCq6CZts2mxs2E41EOy2KlBMn5qAxZfA3BqdjImfjwH3LExVt3Yb2tbV5d8xj5x12prq6mmjEtP3nW//MglcWpHXnclCsMO8KTPEjfpoJXCDUde0sJcxtsTaqA9XpOlI7QSod3LtpThM4ENkzv8/yuvXl8R1bv2E9J04/kS8c84WsNuRDpxzKoVMOLUtbKYr9VP6G2Y291ZNrC3+x91mWRXNTc9qJvlzJbvK2rUrTA01Ufa4q7VCUi4+Wf1Ryexs3b6S5pZndJu1WvG90GWZvxY7M3wVeFpFldPbN+Hrp3eg/qKXYMZtgVdcMOSkfZdu2sWyL+rp64vF4j+SXTixJYK+3GfKdIQWDFpbDzvzw/IcBOOyg4se/QKR0aS62hj8BNvAOsKrDa+vCcUNftdo4STORc5IOyS25V+2cuJNOvm47dnqi19zS3CNd3nL7FgLDAkWFAyhkZ35v6Xu8uPDFTuWqyuI3FqOqLHl3CePHjmfXnYuPHhGIli7MxY7MhwNjfH9mY6t14g5Om0OwPkggGsBusTvvlE5dr0q8KU5AA7S0tqSj1be0tvSYipF4O0HNUTUE6gsLzDYjtsl7/tgzzFL3h4s/BGDFqhXMfWguu+y0CxddfhE3/+pmNjZsZOTwkUX1TyKCBM2rVIoV5jeB4fj+zMbf2JVbtY292GlrDwlAgIwFiZaWFho3NjK8fjixWCztp9DU3DP/SnuDjdPoEN6hOO+0ybtPLuq6VIaqmefOZNWaVey7l4lE8a9F/2Ljpo1MHD+xcCViAtt0DPrYXYoV5ueAp0XkT3TWmf+Y/ZaBg5NwCEQCne3DtvtKXRd3SJCguqZ9ktXa2koikWDtp2szwsdWcnHES/Ijo/4UK8zPLniWybsVFuj1G9czdPBQVq0xmubiNxYD8OCDDwKw3377IUHJujiUQsJCIBJI6/EJEoX3jOWhWGE+BKMfd0xcqWSGuh1wqKPYLWbVztpkdN30zhA7M7ZFfHOcNtqorqk2fsqufpwawXuD5DIjzJEdigsjcOiBuc1l3l+TRYsXpUffidtNZJsR27C5cTMfLDXuOkOHDiVYF8RJGJVMwpIRBiEQDRCsC2aEHGuksSR3waKEWVU/X0oj/ZlUaFjvjumW1hbqaus6hcOyEhaxZIxYa4yG1Q3UVNdgtVp5U/xWmuRHSQKDAgSGFzfB+uiTj9IqQ0fue+i+9PH3f/p9zv/G+QDc/du7GbvtWNa1rGPK1CkAbNmyhWC1CVZDAMLDwiQ3JdGkIiHJCMGb4r/8tyRn6mKXs29yA4ZvVajjCqudaV/O5RCfSCZItCRoWNtAa2MrseYYdlNJv5wloaok3koQ3iGcN5Kol5VrVmYtX/LOEq7/3+szym6981bq6+oZM3oMAOPGjWPp0qXMnDmTb3/72wBI0FUlRAhEA0hYCNWHspoIV7KypPX8Yu0hYeApEXlLRH4oIuMK3tEPceJOeyZTW7GbbLMTxG4fgR3HSWcz7Ug8EQeF1i3GdtwTNuR8ND/QTOK9hImNUSS57Mxr160F4PKLL+fDxR+y2y67AbDjxB3bvygBqK6u5le/+hUTJ040RdFA2uwmIVeww+WxXnSkWBfQC4AxmNgZewPviMgzIvJVEals1pUeJBWwECC5OWl0ZSszTW9qP15HV81UkMKOZb1Jw10NRPeLMugrxadiyGVnTjnYHznVRAmdMsWoE+PHjU9fk01ARYRgtbFWSMgdnQNSkQj7RVuqVdVW1fmqOgOTo28kcBewVkTuEJGxZe9dL5DSjb2TFS9pYe4guJV2qu8qGlfsdTZV+1cVrWIAjNs280f3nQ/e4cQvn8iKVSsAGDxoMATg/O+ez/FHH88Zp56RTttUaHUxEAqkV/q60qdiKTrWnJv67DRM2IHPAH/FxJ5bjkmH9qRb3q9x4u0ZTbORGpE7LnjEYpXzeOsO1qdGDQqN7loa3+0nZOYv/cM9f+Ctd97irXeMT3J9XT0SEkYNG8WtN95KIBxArc4RS3uDYieAD2JMc1/C5CMZo6qzVfVfqroC+D7lzeLaa6ilWJvzRANyN5umArOkjntiNS8f1lqL5VOXE/uP+VKtv8T4E4dGdU2YvUvVEhYmTMjM/BwIBNLqhASM7hsaEipLhtVSKfZJFwLnq+rabCdV1RGRUdnO9UdyqRhggrOkjy2LcDjcJ1SM1n+0oi1K492NxBbFSC41X7rg6K6trn32M5/l3Q/eZddJuxKsCRJPZn5JJSzpFTsJSnrFsxyOQqVS7ATwBq8gi8jnRWRqh2tyTt1FJOruVPlERJpE5D8icqznfL9JBO+d1LW0ttDU3FRRp/piSY3IifcSNPyhfTd0V9WMcy8+l+POPA4JC07QobG5kZEj2v0sgjXBtC+2BNttxeVwFCqVYtWMF0TkYPf4h8B9wFwRuSz/nWlCwAqMT/RgTCT++0VkoicR/E8x8exex2Sz6nFUtdPELpFIpPXhVKreFK1trTQ1N1Xcsb4QsdditL1gfh3ste127TGPjemy1SAV6uDN/77JpF0mce+99zJ48GAAopFoxgjsNbFlSx7U0xT7td0To2qASf0wDZNZ9V/AdYVudvNhX+Upmi8iH2HypAzHTQQPICJXARtEZNeezgUYj8eJxWNEIhEamxqJhCMkkgksy8qbdqy3sDfbNP21iYbfmZE4NDaEtapdDQptW/yofN9D92VsdTrpjJPSx4MGD2LJkiUZ+jL0DQH2UuzTBgAVkR0BSYXpEpFu7Ut39etJmCys36ZDInh3E8AemH2HPUpqgheLxzIicHZMmN4XWH/xeuJvtOu0tV+opeE2I9jj/jmuoKmssamReCLOyOEjuezazj+yo0aNYt26dQwZMoRhw4aVt/MVoNiv1kuYpJM3YDKq4gr2hq42KCJhYA7wZ3fkrcMkfveSMxG8iLwuIq+vX59992+pJC3jb2xZFslkMq1abGnYQmNTY0Xa7C5eQa46qIrBXx2cfl+M7/IJM0/ggKMP4KVFL2WUf+Mb3+C5557jqKOOAmDQoJ7Lf10KXYkCejGwHviVW7Yr8L9daUxEAsDdQAI43y3uU4ngHcfJGSKrL43Mqbx9KSKTIkhUGPfUOKx1VsFFiVgslo5t8ds7f5suHzZsGNdccw0Au+1mlqynTp3auYI+SLFecxuByzqUPd6VhsT8d+8ERgHHqWrK3e9t+lgi+ObWntnOVAptizLNgYFBZiQOjggWzN3X0NiQDlgIsHDxwvRxXV27d8KMGTM48MAD2XXX/pE8tysrgHtjYjSPIL2ACap6RZFV/A7YDThSVb2fxMPAr9z82Y/TBxLB97Z1ohhiC2MEhgZwNptJaTa1wrZtYwMOmHMNjQ1898ffzVgYGTN6DKvXrgagvr6eq666Kn0uGo32G0GG4k1zszGWi8OBHwKTMWrHTkXePwE4F+OktFZEmt3XWW4y+FOAa4HNwAH4ieDzogkltiiWETYgmzDv/LmdufDyC9Pvb/jNDZ02o3527/bg4otfX8yqVf13j3KxE8BLgemq+kWgzf17KlDUzgBV/URVRVWrVLXO85rjnn9GVXdV1WpVnaaqH3fnYbYG2l5to/nxZpxGh9pja9PlHYU55aI6/6n5gBmln3zmSfaZvE9GiKz9DtwvfVxbV8vmzZsr2f2KUqyasY2qpr7SjogEVPVJEZlTqY75dKb1+VbWX2ysOMGRQar2r6Lui3U0P9zcaaVv4+aN6WNV5dqbrmXTlk1c86NrOPbIY/lw5YcMHz2cE044IUO1OOuss3rkWSpBscK8UkQmuiPm+8BJIrIBY5XwqQCtC1qRoNA4r5GR144kuSLJxmvaBbR6ajUSEIb9ZBhDvjOkUyR8b5iteQ/P46777uKMk89g+hHTERF+9MMfpf2ML7roIkaNMq41c+bM4eKLL+6BJyw/xQrz/2Ambx8D1wAPAhFMpCOfMpP8JMn6i9rt6CumrSC0XQinwWHwtweTfD/J4K8bm7KIZE3p4I3CmVoQ+f63v5+eDHqXuS+55JL0cX+a8HWkWNPcXZ7jJ92Vv4iq9n0bVj9Dk8rqL63uVG4tt6g/s54h5wzpdO6h+Q8RCoU4cfqJ6bKPV3wMwMnHncwjTzzChPETMhyGcq0O1tbWZi3vD3TNpcpFVRP4KkZFSCzN/W+NTu4citayLC650oys48aMY8L4CaxZu4bf/fF3TN59MheeeyHNLc18faYnJKDkFubFixczbdq0kp6ht+iWMPuUF7WV5r82U/uFWhJvG2Ee89AY1p23jtpja6k/vZ7YqzFqjsiMFff+svc58cvto/GpXzuVcWPGpVf2LjnvErYbtx2333R7xn35fDZOOumknOf6Or4w9wFan21l0y83Ya2z0DZFaoXQdiHG/m1selm67vjO+4Zvvu3mTgs83mQ4O26/Y/YG8xhkn3nmGXbaqajlgz5HQTuziARE5HARKS4kjk9OrLWW2QHeZMJ8qZqEOW0LzYJo412NNM1rIjTGJNDJ5l+xpWFL2kdk+UqTmnHK/lP4xRW/6HTtTttnF8p8I3N/WP3MRUFhdqPlP+rqyT7dxN5ks+r4Vaz71jpWTFtBwx8baJnfwpoZa2h5NDOoTGib9h9Mb6rhxqZGPnv4Zzn/h+ezfuN63nn/HS745gX86dd/4tQvnEoomPlDO3RIDg/dPJ/6F7/4xa4/XB+h2BXABSJSnmQXWympiV38deOR1/qPVmKL3O1WIRj2w2FEP2MmeHaDEeCFry9k1ym78vy/ngfaR+Inn30yXXbkYUcSCUcIBALpHSEA9995f86+5BuZ778/9319nWJ15k+AJ0XkUcz2p7QvZBccjbZqUtE40++XJUkuSzJn/BzmbpjLKR+fQvOoZi7gAqo+WwXA8y8/j23b3PibG5l28LSM0FmPPPEIVdEqdp+0e7qsrc2oK9dffj377b0fOckzhE2eXFxI275IscJcDTziHg/I0Fyl4LQ4aEIJDAmw+ebN1B5dm5HRyWly0tE4M4jAHSvuAOAv8/4CwEV3XcTqwGree/09lvx3CWCsFm1tbenwsQCvvPYKn9njMwSD7QsmqXBghxx4SN7+FtqB0l8pdtHka5XuSH9m9amrsT+1GTt/LE1zmmia08SExWaDeeyNGOu+bpaWqw+rJjQmRO3xtUR3i6KWEjoolBG+4NCzM0PKbj9hez765COuueEaPvjwAyLhSHo3zK47Za7WXfmDK3nh5RcYu23+4FL5hHnJkiUcfXTHyMX9g674M++G8ZQbparni8guQFRV36xY7/oBaiv2p27Irnfa58jWWovQ6BBtz7e7btefVk/1lHa3zYSTyBDkbNxw9Q383x/+j3mPmA3re+66J+PHjufJZ5/k1BNPzbh21pmzmHXmrGzVZJJHzTj99NML399HKdaf+TRgATAW+KpbXA/cVKF+9Rvib7VvsVr/g3Z/iFXHr8JaYxF7rT2mRtUBVRn3vr8sM8Htrjvvyve//X0Axo8dzx233ME+k/fhmGnHpK857xvncf1Pr2feHfPy68V5yDcyP/zww92qsy9Q7Mh8DXCUqv4/ETnDLXsD2Ksy3eo/tDyRaVaLfiZK/E0j4KtOMDru4G8Opv70et58503+/ca/mXXmLP40909ce9O1Gfeef875HHfkceyy0y5M3n0yo7cZDRj/itVrV/OV07+S9q/IliiyKCR7tM4UkUj/XU6QYjZpishGYISqqohsUtVhIhICVqtq/vREFWK//fbT119/vax1tmxp4dMPPy14XXJ5koa7Ghh6/lDWnr2W4Igg8TfiDDp7EEO+NYTlBy5PXxscHmTbedsSHBrk4OMOZs26NUzacVLGqHzP7+5hw6YNTD9iOpFwZYVJQkJkRO42li5d2msrgGPHjh2rqp29rIqk2JF5MfAV4C+esjOBV7vbcH+gbWEbmlBqpmb6RDTOaaTl0RZiC2PY62zqz6xn1B2j0j/fYx4eQ3B4EHuzzQP/fADnGYczTj6DjZuMP/L7y95nv733Y9iQYZx20mlM2X9KRUK8ZqOQJePRRx8d8P7M38Vkm/oGUCsiT2GCuPTPaW8RxN6I8el5ZpQecd0Iao6qQVsViQhtL7uhsNbZEITqQ4yj/MbNG1m+cjnr1q/jqdufQkR45IlHALjiF5nm+Osvvz6370QFKRStc999s+cz6Q8Ua5p7V0R2BU4A5mMWTuYzgN1AN169keDIICg03NlA00NNWCsswtuHsVfbDPvRMOJvxRl8zmDC401aspt+exMPPPoAQFYrxawzZnHqiafy7IJn2WHiDj36PGkKBAVtaSkpR06vUpQwi8gtqnohcL+nLAo8BkyvTNd6B7vBpuXvLVifWAy9cCh2g03jn9ojGdnrbOpn1FN/mnl5eWnRS2khvu4n19HS1sLI4SM54egTeH/Z+0zacRKBQIA9dt2jR5/JSyE149133+X444/vod6Ul2LVjM+IyNWqeiWAiFRjYlysqFjPeom5F8xlwtsT2I7tiO4VxVrtRqAfH8JaYY5rjuycg3rFqhWsWLWCcChM0kpyyIGHMG5M+2JpV/JIV5JCiXH684bWYh2NTgKmi8j33YhDTwHLMGG7BgwrVqzgirevYBazsLGJ7B6h6oAqIntEGPnLkYQmmO9+ZJdMa8B7S9/j8usuB+CXV/6S2bNmF1yF6y0KCfOcOf13w32xwcabgGOBrwH/D3hDVb+pfSn4Whm458/3pI8/+cEnSEgIDg2y7V+2JbJLhBG/H8EN+9/Ai//JDKRy+XWX8+LCFxkzegwnHXsSP/ruj3rMOtElCtiYwWRW7a/kVDNE5Josxa8CxwObU+cHitfcxd+5mJt+dxP7sA8fVH3AU+8/xefJTEz7+KLHefy1x1myegnPP/Y8YAIQLvnvEg4+4GBu/vnNfVOIAyaWckbe7xykUqL1R/LpzONzlD+Z51y/ZacVOzGTmZwy9BRu3+t2Xlz4Yjouc22N2bH84GMmyfnyVct5dsGz1NXW8es7fk0imWD2V2czYtiI3nyErEhQCA0NEQgFsJoKJ0B94okn0tE/+xs5hXlr85T7ynVf4bjjjyOxXYKpL0/lqeef4jd3/oZbfn8LZ5x8BmfPOJtF/17E/vvsz2v/eY3Lr7s8HWjlpGNP4pAD8rtd9haBqkA6B0mgqrBW2Z9H5mIdjXZPZZMSkToRuVpErhCRztP6fkrd5DpGnDmC0GgTf2L4sOHc9LubqK6q5r6H7+OY044hGAhy889v5pjDj2Hd+nWEgiGu/MGV/OqqX/VN9YLMYC/FpG1Yt25dwWv6KsWa5u4FzgDWYaLn7wLEgN9jlrkHFLU1tfz51j/zwKMPcPLxJxOLxXjmhWc49MBDGTN6DNtvZ1Ie3nHLHUw9qG8H4u5qjuoPP/ywQj2pPMU6Gm1R1SFuwPC1mHwjbcBHW6OjUWtbK6//5/U+L8gAkW0iXdpZsnbtWkaPHl3BHuWmVEejYu3McRGpBz4HrFDVDUAcqMp/28CkprqmXwhyvshFuejPduauqBnPYRzyb3XLPgt8VIlO+ZSHYiZ8Hdlmm175oS0LxToaXSQiRwNJVf2nW+wAF1WsZz4lE6rvesCqvfbqv/stiv7qqurTHkFGVV9X1ecq0y2fUpGIdGsX9j/+8Y8K9KZnyLcC+HdVne4ev4gnVoYXVe0HyuPWR3hIuFv39Zc0adnI9zvk3VVyR6U7IiLDMKnVjsYky/yxqt5b6XYHIhLs3qgMxjS3//7d3F/Yy+RbAbzXc/znHujLbzDO/qMwWakeF5E3VLXX8gH2V7oz8UuxcuXKwhf1UboSN+NojJBlxFYth6OR61Z6CrCnG43/JRF5DLMg86NS69/aCNYW2E6ShwHvzywitwL3APtinIy8r3IwCbBV1RtI4g3M4oy3HxXPnX3L72/hlt/fAsDhXzycDz/5kCXvLOHEs0xQ72tvupY77jZa14HHHMi69etY+PpCZsyeAcBlP7+MuQ/NBWDyoZNpbmnm2QXPcs6F5wDwvcu+x6NPPgrADvuarVOPPvko37vsewCcc+E5PLvgWZpbmpl8qIn7NvehuVz2c5OXZMbsGSx8fSHr1q/jwGNMLMs77r6Da2+6lmBNkBO/fCJL3lrCsmXLOOQQ4y9y4403cuONNwJwyCGHsGzZMt58802mTzebhK6++mpuu+02wOjMa9eu5eWXX+bUU02QmUsvvZR77jHusZMmTaK5uZmnn36aWbNMwJnzzjsvHW9j7Fjjx/3www9z3nnnATBr1iyefvppmpubmTRpEgD33HMPl156KQCnnnoqL7/8crc/sxRdCTWwt6pWZGeJiBwKPKCqoz1l3wTOUtVp2e7pzRXAPoe7OBIaFMJutQkP7d7kD2DevHmcccYZhS+sAD0VamAjsKW7jRRB0cngfTojYSE8OAwBEKs0h6cdduiljbZlIKeaISI7pF7AjcAcEZniLXfPlYP3gZCI7Owp24teTAbfZyhCNgOhgLFgiBTlGZePBQsWlHR/b5JvZF6apeyEDu+VgpvXC6OqLSLyEHCNiJyDmWieBBxUat39nWBtELvFzmHlN3jdPL3H3eGoo44q6f7eJOfXWFUDRbxKFmQP38HEgf4UmAt82zfLQSASyAjcIkHpNFpnCHOJftVvvPFGSff3Jnl1ZhE5HXhBVSvusa2qm4CTK91OvyEAqBFUCQlqKSjpSZ4Td8x1QnonSTn49NN+OAF2KTQB/Dmwo4gsw4S0fQFYoKqfVLxnWznB6iBqa1oPFjEjciAawEk6ZnlJi9s9AhCNRnEch2QySwR/DwPWzqyqk4AxwE8wzvgXA8tE5BMRudvVb30qgISF0GAz1khECFQHCETdvXyRAOFhYQJVASRSnFpRVVVVVLharz9zV8LbBoPBXt86VkzqtHWq+oCqXqCqewMjMEvPR2G2TfmUirgBDV1ZCEQCBKvahSMQChAIBwhE2oU5EA4QHhIuerWvWGH2muY65tGuqsq+F0NECIVCXRL+VEL6clLQzuxuldobmOq+DgJWY+LOvZj7Tp+iEON3HKwJktiUQC3tkkXCOxoGg0Ecx6HjQlgkEiEYDBYlbKNGjep0n23bBINBqquricViGdeLCEOHDsWyTMLOeDzescpOpPpR7gSaeb8eIjIfWIVJ9zACuB2YpKr7qup3VfWBsvZmK0NCYkbhGjO6SlAI1YUKhp3NxaBBgxgyZEin8upqk0clGAyms1PV1NQwcuTITqrBK6+8YvoiknF9KBTKyGyVIhqNUlVVRU1NTdbz2QiHw+lryxmpv9BYvwtmr99HmNhyS91QXT5lQMJCsL5dAAKRgIlzEe3eT3BqJO0ooN6f9JSqEA6HCYfD1NRkRos47rjjACOk3r9ewYb2X4RwOJxuo1hh9l7bQXVpLaqCXPXmO6mqOwMHYqIY7Qv8VURWisg8ETlfRPYupfGtFQmJmdSFAxlmtWB1EAkU9kUOhTprh6nRLhgMUlNTw6BB7d4BXmGORCJEo9H0iNhxZEyNzKkRPiXsqZFZRKipqcl6f7Z+ZSMlzIFAgOrqau+Xry3ffQXrLXRBlgngnpi0EJe7fwcupU7Oc9yfUi2663ecUhu8pEa4QCBAKBTKmLx5R8xIJEJ1dXV6RI1Go4hIWig3b95MKBRKfwFSunZKqAcNGsSgQYPSdabqSV07fPjwgv0PBAJEo9H0ly81+gOF44floTsTwEOAIcDrwB9LabyvIyEBxSxYdINgdRC71c4sDJh6g1XtAhYOh4lGo7S2tuI4Tva63IkYGGFubm7OmOilhEtECIfDaQuDZVkZI3MwGMz4aQ8EAml9N5FIcNZZZ3UaYYcPH54ePVNfktQvQEerRDQaTbcbiUSyTvJSX7j6+vp0XSKCqtqdLu4ChSaAjwObgJcwvhJvYhLzDFbVg1X1slIa7+tISIwtt4N1IdcEzWvzDQ0OZS4zB4WauhpqhtR0qq+6uppoNJrXXOUdDTvqr6myFKnRMiWUHevt+L6+vj59z5w5c7wjpel7FvtxdXV1RuJ5L6l2s/2CeOtLfamCwWDRKko+CtXwInAt8Jqq5l86GoAEwkYdcMIOyc3tjx+IBrCtzEGktq6WNruNYDiIjU2wOohjmVE2WBs0iyCBEFItncxnVVVVaT3SsizC4XCnlbrUT31LSwsiQl1dHU1NTenROtvkLKWXFnxOt+1QKMTkyZM7TQqzkU/4QqEQIkJVVRUtLS1YVqb2kO3LVezkMW+f8p1U1V+U3EJ/xv2fp0dS1wk+9T4QCKABJaABautribfEGTFqBJtjm41AkjQTvYixUIQ0RLg2TNAOphPheEcl74hVU1NDQ0ND+ppgMJhh0koJXDKZpK2tLasweHXfQgSDQYYMGcLo0aNLXskLh8PU1tamn80rzCLSqU/lGpnLvwzTz6mKVqVNYyl1QgJG3QiEXb9htzwaiRKtjxIKhghXhYnWRakabEbZujqzVTJYHUxfH4wGqamrydBZU9cBGZOq2tpaAoEAkUgkPTnrqF7U1NQwePBg6urqsgpgsQslQHoSuHDhwuL+UXmorq5uf/4OX7L6+vpOwhwOh8tiby796zDAGDJiCE1WE23r29IRNEWE6LAoiS0JRo8fTSAaYHXzaiJVEWSQoFElPCTMkOQQJGhWxFI2X63S9KQsXGP00kgkkprwZHyIHW23tbW11NbWplWSUCiUYT1I4f1CeAmHw51WAwsxY8aMLl2fi5TAdhxxs/2CBAKBnEvlXcEXZg/hqjDBSJAap4ZwOEybZcyeVVVViAh21CZSa4Rvm/HbENIQiVACqskwb6UmULW1tTiOk77fq06kfCW8whkOh6mrq0t/4KnZfgoR6eQvkY/UEnRXmD9/PrvsskuX7inUBy+V8MlI4Quzh0iVEcZqu5qamhqq4lVs2bKFcDhMIBDAqm/X/arqzUiiydwjX3V1NYFAIKsKMHjw4E4fbHV1dUHhq6QwAEX5VnSFaDTK8OHD2bRpE6paloleLnydOQve1amUnprrJz5XeaqeXJOpSgtld5k5c2ZZ6xMRotFoenGmks/dN/+jfQivMGebcfe2D2+5ueuuuypSbzQaZciQIb6a0ZukfBJSS7ADnUolgk+pXJXEF+YCeE1J5bCFbq30hFrlqxkF6Ku6baVYvLj/+o5tXZ+UT0HOPvvs3u5Ct/GF2SeDe+/tvyGxfWH2yaA/T3J9YfbJ4IQTOkZg6z/4wuyTwdy5c3u7C93GF2afDPpzIviigo33RURkPdBfwoSNwCQd6ut1VqreYuvckMpw1h36rTD3J0TkdVXdr6/XWal6K9XXjvhqhs+AwRdmnwGDL8w9w+39pM5K1Vupvmbg68w+AwZ/ZPYZMPjC7DNg8IW5jIjIcyKiIhLylA0TkYdFpMXNODCzwz1HiMi7ItIqIv8UkQm90O+8fewv+MJcJkTkLLJvdvAmuD8L+J2I7OHeMwJ4CPgpMAwTv29ej3S4yD72K1TVf5X4AgZjEnMeiMnYF3LLazFCMslz7d3AL9zj2cDLnnO1mLCuu/Zg3/P2sT+9/JG5PFwH/A5Y26G8UIL7Pdz3gEnuiQnq3pOjYqE+9ht8YS4REdkPOBj4dZbTdUBDh7IGoL7I8z1BX+hDWfB3aHYRVzdOZdl6ERgOfE9VrSxhBwoluC90vifoC30oC/7I3EVUdY6q1qlqHTAD2A+YJyJrgdfcy1aKyKEUTnD/tvseABGpBXb0nO8JCvWx3+CvAJaAm1VglKdoPPAqMA5Yr6oJEbkPMylMJbh/AjhIVd8WkZHAUuDrwOPA1cBhqnpgzz0F5OtjT/ajZHp7BjqQXsBEPNYMt2wY8AjQAiwHZna450jgXYwV43lgYi/0O28f+8vLH5l9Bgy+zuwzYPCF2WfA4Auzz4DBF2afAYMvzD4DBl+YfQYMvjBXEBF5W0Sm9XY/OiIiZ4vIS73dj3Lj+2aUgIg0e97WAHEglbr1XFXtd55n/RlfmEtAjX8GACLyMXCOqj7Tez1K9yWkqlbhKwcWvppRQUTkYxE50j2+SkQeEJF7RKRJRJaIyCQR+bGIfCoiK0TkaM+9g0XkThFZIyKrROTnIpI175hb94Nu3Y3A2V25f6DgC3PP8gXMLo6hwH+ApzCfwVjgGtpdSwH+DFjATsA+wNEYR6BcnAQ8CAwB5nTj/n6PL8w9y4uq+pSrAjwAjMRsT0oC9wETRWSIiIwCjgUuVNUWVf0UuBk4M0/dr6jqI6rqYPyRu3p/v8fXmXuWdZ7jNkzUS9vzHszOjzFAGFjjcfgPACvy1O09N6Eb9/d7fGHum6zAWEZGdGEi53V/7M79/R5fzeiDqOoa4GngRhEZJCIBEdlRRA7rifv7K74w912+CkSA/wKbMZO7bStxv7u4c1ZJve0D+M75PgMGf2T2GTD4wuwzYPCF2WfA4Auzz4DBF+Z+TqnunCJyrojcUsYuISJRN0zvNuWstxBbjTB7nX58DCISAS4HflXOelU1DvwR+GE56y3EViPMPlk5CXhXVVdVoO57gVki0mOZ5bcKYRaRu4HtgL+JSLOIXOqWn+guGGwRkedFZLcc94uI3Oy6ajaIyJsisqd77ngR+Y+INLpunFd1uPcQEXnZbWOFiJztlleLyI1upPoGEXlJRKrdcw+IyFq3fIE38LeIDBeRx9z2XsXEpvO2t6uI/ENENonIeyJyep5/zbHAC557J4qJ/P81t6+bReRbIrK/+8xbRORWz/U7icgLbj83iEg6ULqqrsQs1vRcqLHeDqnUgyGoPgaO9LyfhAlHdRTGKedSTNy3SJZ7jwEWY9wrBdgN2NY9Nw2YjBkYPoNxJjrZPbcdJprmDLeN4cDe7rnfYMJxjQWCwEFA1D33dUxI2ShwC/D/PH25D7gfEyR8T2AV8JJ7rhbjl/E1jN/NZzFpfvfI8T95DTjN834ixsfjNqAK4zYaw4Tu2sbt66eYeHgAc4GfuM9eBRzSof7HgO/22Gfc20LWi8L8U+B+z/uAKxjTstx7OO2R8QMF2rkFuNk9/jHwcJZrAhgvub2K6PcQV8AGu0KfxBNZHxPoPCXMZ2DcTL33/x64MkfdHwDTPe9TwjzWU7YROMPz/q8Y11KAv2By/I3LUf8c4Iqe+oy3CjUjB2PwJJJX4we8AjP6ZKCqzwG3YkbTdSJyu4gMAhCRA8Qk1lkvIg3AtzCJz8FEBV2Wpe0RmJGs0zkRCYrIL0Rkmbtr5GPPPSMxI67XlfMTz/EE4ABXHdgiIlswOUpG5/gfbCZ7UPGOrqod36e2i12K+aV61VXXvt6hnnpgS462y87WJMwdnVBWYz58IB2edjxmdO58s+r/qeq+mPQIk4AfuKfuxfycjlfVwZif6JQT8Qo66LQuGzA/39nOzcRMzI7EjMYTU10E1mN2j4z3XL+d53gF8IKqDvG86lT129meCXjTfZZuoaprVfWbqjoGOBf4rYjs5LlkNzxpLirN1iTM64AdPO/vB44Xk7osDFyM8QF+ueON7gToAPe6Fowgppzq64FNqhoTkc9hhDHFHOBIETldRELu5G1v91fgj8BNIjLGHY2nuDP/ercfGzE7vq9LVabGkf8h4CoRqRGR3YFZnvbmA5NE5CsiEnZf++ea2GLiMHfbLVREThORce7bzZgBw3bPjcWEyl3Y3fq7TG/qsT35wox2yzE/e5e4ZV/EuEg2YGb1uSZKR2BGsWbMqDoHqHPPnYr5qW/CCNOtwD2eew8FFgGNmJFzlltejdGvV7ntL3DL6oBH3fo+wbhyKrCTe99It51GTGDzn+HqzO75XTCBy9djvhDP4U46szxX2P2fjNFMndkbX3olnnkEcA9wuXv8P27/mzEq02zPdT8AburJz9h3Ad3KEZHZwO6qemEZ64xi1IupavYf9gi+MPsMGLYmndlngOMLs8+AwRdmnwGDL8w+AwZfmH0GDL4w+wwYfGH2GTD4wuwzYPj/U/dti0G72nsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(2.5, 3.5))\n", "for lab, val, col in (('Data', dmerged, 'm'), ('Shuffle', rmerged, 'k')):\n", " m = np.median(val, axis=-1)\n", " l = np.percentile(val, 25, axis=-1)\n", " u = np.percentile(val, 75, axis=-1)\n", " x = (np.arange(m.size) - (m.size // 2 + 1))*5\n", " plt.plot(x, m, '-', c=col, lw=1.5, label=lab)\n", " plt.fill_between(x, l, u, color=col, alpha=.1, lw=0)\n", "plt.vlines(0, -1, 1, color='k', linewidth=.5, linestyle='dashed')\n", "plt.hlines(0, -600, 600, color='k', linewidth=1, linestyle='dotted')\n", "plt.title(f\"N={epochs['left']['N']} leftward, {epochs['right']['N']} rightward\", fontsize=9)\n", "plt.xlim(-800, 300)\n", "plt.ylim(-0.15, 0.65)\n", "plt.xticks((-400, -200, 0, 200), (\"-400\", \"\", \"0\", \"\"))\n", "plt.yticks((0, .2, .4), (\"0\", \"20\", \"40\"))\n", "plt.xlabel(\"Time rel.\\nto saccade (ms)\", fontsize=12)\n", "plt.ylabel(\"Whisker asymmetry (%)\", fontsize=12)\n", "plt.tick_params(labelsize=12)\n", "# plt.legend(fontsize=12, frameon=False, loc='upper right')\n", "for side in (\"top\", \"right\"):\n", " plt.gca().spines[side].set_visible(False)\n", "plt.subplots_adjust(left=.3, bottom=.22, top=.93, right=.95)\n", "\n", "if saved == True:\n", " figpath = figdir / \"saccade_traces.png\"\n", " if not figdir.exists():\n", " figdir.mkdir(parents=True)\n", " fig.savefig(str(figpath), dpi=300)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "id": "24280930", "metadata": {}, "outputs": [], "source": [ "class BoxSummary(namedtuple('_BoxSummary', ('n', 'q1', 'q2', 'q3', 'min', 'max', 'outliers'))):\n", " @classmethod\n", " def from_values(cls, values):\n", " values = values[~np.isnan(values)]\n", " n = values.size\n", " q1 = np.percentile(values, 25)\n", " q2 = np.median(values)\n", " q3 = np.percentile(values, 75)\n", " iqr = q3 - q1\n", " lower = q1 - 1.5 * iqr\n", " upper = q3 + 1.5 * iqr\n", " vmin = values[values >= lower].min()\n", " vmax = values[values <= upper].max()\n", " outliers = values[np.logical_or(values < vmin, values > vmax)]\n", " return cls(n, q1, q2, q3, vmin, vmax, outliers)\n", "\n", "def plot_box_around(x, boxsummary, ax=None, color='gray', mcolor='w', ecolor='k', ocolor='k',\n", " boxwidth=.5, mlinewidth=1.5, elinewidth=1, marker='o', markersize=2, boxalpha=.6):\n", " v = boxsummary\n", " if ax is None:\n", " ax = plt.gca()\n", " ax.add_patch(Rectangle((x-boxwidth/2,v.q1), boxwidth, v.q3-v.q1, color=color, linewidth=0, alpha=boxalpha))\n", " ax.hlines(v.q2, x-boxwidth*.45, x+boxwidth*.45, color=mcolor, linewidth=mlinewidth)\n", " ax.vlines(x, v.min, v.q1, color=ecolor, linewidth=elinewidth)\n", " ax.vlines(x, v.q3, v.max, color=ecolor, linewidth=elinewidth)\n", " if v.outliers.size > 0:\n", " ax.plot((x,)*v.outliers.size, v.outliers, marker, color=ocolor, markersize=markersize)" ] }, { "cell_type": "code", "execution_count": 10, "id": "4e51ee34", "metadata": {}, "outputs": [], "source": [ "def plabel(p):\n", " if p < 0.001:\n", " return \"***\"\n", " elif p < 0.01:\n", " return \"**\"\n", " elif p < 0.05:\n", " return \"*\"\n", " else:\n", " return \"NS\"" ] }, { "cell_type": "code", "execution_count": 11, "id": "da18aedf", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAADtCAYAAAC/HcCwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAApeklEQVR4nO2deZgV1Zn/P18xjUsjhCUKKuAGKuJGq1EMghqCisuEaAQ0wYliYtQYnTHRQSGaqIOGMb+skhiMCZCoaAaIJkxGECVopplEGSMSF0BZDLZKbBWa4Pv741RD9eUudatv3Vt97/k8Tz1ddarqfN9Tt+utsx+ZGR6Pp3bZpdIGeDyeyuKdgMdT43gn4PHUON4JeDw1jncCHk+N452Ax1PjeCeQIJL6SzJJ54bCXop4762SVkv6fSjsGElLJC2W9LikA4NwSfqupCclzZfUPUt8iyTtl0dvH0lLJS2UdLykYcWlNhqS9pO0qMA150h6JkjP+IxzAyRtlXRyBK27JfXKc36CpElZws+T1LdQ/HGQdJGkKUnEHRfvBJJnBXCDJBV53w+AERlh64FRZjYMuAv4RhD+KWAPM/sE8ABwfQw7RwALzGwEcDhQEicgqaj/seD6O4HTgVOBKyV1C11yE/BEhHg6mdk1ZraxGP2A84CSOAFJnUoRT5LsWmkDaoC1wErgXODXUW8ys/WS+meEbQgdtgD/CPaHA/OD/XnAF/PFLel24CSgDvgW8CdgMtBZ0r7A8UAXSacDnwe+Y2bnSfo20MfMxkr6PjAT6AzcjPtfegv4rJltDnI8DwAnBjmh2cG1zxVIek9go5m9G9i6MrBngaTjgQ3AthzpmgL0B7oDsyVdDlyEc54/B/YFlgKfMbODg9sGS5oDDAS+BDQBo4CjgjTMBz5qZndLagR+YmY/kvRHYChwFTAa6AI8amaTg9/tQdwHYKuku4D7gI1B/K8UeAZlxecEysNtwNfDuQFJJwZZ9Mzt1EKRSdoT9/LeGQR1B94O9t8JjnPdOwr3T30KcFoQzzrgDuBeM7sUmBbsDzez1cD+wRftMKA+SMfxwP8AfzSzEUEuZAVwQSC1KzAvyFl8AXjKzE4HlhVI3kagp6R9Je0FnBxKz6TAznxsMbNzzGx2KOxc4O9BmufR9uPX2czGABOBr5jZX4DfAleZ2fnAfwOnSfpo8JxOkbQ38JaZbQXuMbNTgROAT4aKEf2BL5vZPwO3B3GfBWwpYH/Z8TmBMmBmr0tahstmtoYtxX3Bi0LSR4BfAbcH/7DgvsDdgv2u7HAI2RiM+0deFBx3BnoUkF0GfBL3FXst2N9oZluDeopvBvHsDfw9uGcb8HSwPwB4KNh/Brgsl5CZmaSJuC/3e8ByYJ2ks4BGM2sqULL6Q5awQ3AOq1U/3Fe+1SmtIctzCH673jiHOQ84O9hfGFwyRtKlQZwHAvvjcn//Z2atz+IQ4I8h/Zx1M5XA5wTKx+3A11oP4uQEgvLyL4Bfm9mvQ6eeAM4M9s8kf5n5eVzZf7iZDQeONLM3M65poe0H4nFc/cOiYP/WYB/g34DJwVd2LtD6hprtGJjyV6Ah2D8ulJ4uwRe2DWa2OPi6jgPqcS/O0cBwSb/FOaG7JPXLkr5sRYWXMvTDXiTsEFrDM9PfCPwrLu3PAV9hhxO4FVcnMwJ4NRRH2I5M/VThcwJlIviiNOLKmwVzApKuBC4EDgtaCC4HjgHOAvaWdBGw3MyuAn4HjJb0JO5L/Lk8djza6oBwL8DrwMUZly3BVcgdAVyJ++efFdizATgK+HJw7S+BeyW9CGxiR04gzI+BByR9Evi/UPiFuLL0tIy0T8W9LP8AbjCzLbhiy7eC8/fhyuarc6Uzg18D50t6ApcjKJQlnw/cIukFM7scVyQ4w8xelvTfuHqA1hzEw7jntQJozhHfjcBPJTUBmQ634siPIvRUCknfAb5lZn8rg9ZHguLLUJxjGZ20ZkfBOwFPTRC0APTE1V1cbmbPVtik1OCdgMdT4/iKQY+nxvFOwOOpcbwT8HhqnJpxAqNGjTJck1gqtilTpnhdr1vOLSc1UzHY0NBgjY2NlTZjO+vWraNPnz5e1+uWi5zdLGsmJ5A2pk+f7nW9birwTqBC9OhRqLu+1/W65cE7gYhI4rrrrtt+fNdddzFlyhQAXnzxRYYPH87RRx/NYYcdxsSJEwvGN3z48LJrRtXNpLXI2KodpwhZSDeXRvg4jh21phsHP3YgIp07d+bhhx/mhhtuoGfPnm3OXX311Xz1q1/l3HPPBWD58uUF45szZw6DBw8uq2ZU3UxmzpzJunXr2Lx5M1OnTqVPnz5cdNFFRcVRSDebRuYxULQdtaYbi1aPU+3bkCFDrD3sueeedtttt9mNN95oZmZ33nmnTZ482czMBg8ebI2NjUXFt3DhwrJrRtXNxqxZs0ySzZ49O9b9UXQzNbJpFmtHrenmIee74VsHIlJfX8+6des48sgjefbZZ/nxj39Mc3MzU6ZMYcaMGVxzzTWcdNJJjBw5kksuuYRu3brljW/WrFmMGzeurJpRdbPd8/rrr/PWW2/RvXt39ttvv1hx5Lsnm0bmMVC0HbWmm4fckzDk8xDVtJUiJ2BmdtNNN9ktt9zS5qtsZrZ27Vq799577ZxzzrGBAwfa5s2b88YXvrdcmlF1M/nwww/b3Nt6XAyFdHNphI/j2FFrunnI+W5U/OUs11YqJ9DU1GT9+vWzKVOm5PyhBw0aVDCrvnbt2rJrRtVNAq9bcXK+G751oEi6d+/OBRdcwL333rs97Le//S1bt24FYMOGDTQ1NbHvvvvmjaeYduRSaRarW0q8bnrxTiAG1113HW++uWOCmAULFnDEEUdw1FFH8alPfYo777yTffbZJ28cvXv3LrtmHN1S4XXTi68YrBDLli1jyJAhXtfrlgvfbThtzJs3z+t63VTgnUCFGDlypNf1ukWRVD1DJCcgaX9Jo4N11EZL2j8Ra2qIlStXel2vWxRJOYGc3YaDRS4uD7YDcXOnv4ubIvpgSa8CPwKmm1lLItZVMatWrfK6XjcV5KwYlPQXdsw3/4yZbQud64Rbhmo8MMLMBpXB1naRtorBWhvn7nXbT0NDA+34H45VMTjczK40sz+EHQCAmW0zs6VmdiUxltLy1F77dVp110xdw6s3v8qaqWtKFmfUa9JCTidgeRaEkNSzdXFNi7f0c83Tv39/r5sC3bcWvMXqW1fz1oK3IscZ5QWvVHrjUOza8cMkrcatx7ZR0vnJmFX9DBgwwOt63VSQ1wkES2CHmQwMM7M+wCnA3XGFJV0pqVHSlmBtufC50yStkPS+pIXhhSfl+HdJTcE2Nbzkd0dhwYIFXjcFut1HdqffTf3oPjLnau6J6KaJQpOKLJZ0m5nNCY63AvtIWotbXrk9rQLrgG/iVnTdvTVQUk/cIo+X4paCvhW3FPfHg0sm4pb4Pgo3i+p/Aa/gWio6DGeffbbXTYFu3+v7VkR3zdQ1bGveRqf6TonZEJVCxYFTgdMkPSrpINzyzP8Pt/Lst4EvxBU2s4fNLa/dlHHq08DzZvagmW0GpgBHSTo0OP954Ntm9rqZrQ3smBDXjkpRqZYKr5sO3W3N21h962q2NWdbSb285HUCZrbJzK7AvYj3A5/BFQf2MLMjzOzxBGwaBGxfLNLM3gNeDsJ3Oh/sp76JMpP169d73XYSpYIurentVN+JprOa6FTfqUwW5aZgxWBQ3n4FVwfwJrBU0hkJ2lSPW+c+zCZcJ6Vs5zcB9dnqBSRNDOodGjduTFcjRtSJQWtFN6mmurSmt+/1fbl9w+0VLwpA4YrBC3Av/nJgFfB/wJnA5yQ9LGm/BGxqBvbKCNsL11sx2/m9gGbL0uvJzKabWYOZNfTq1SsBU+OT1nbzSukmlT1Oa3rTRKGcwH/gegT2Bs4CbjOzDWY2FvghMDcBm57HVfoB21soDgrCdzof7D9PB6PWmq4K6Xaq70S/m/qVPHuc1vSmiUKtA1twLQLgauI3t54ws/+S9ERcYUm7BvqdgE6SdgP+ATwC3ClpDPAb4GbgOTNbEdx6P3CtpEcDm64DvhvXjkpRqSWq0qqbVLY4relNE4VyApcBv5L0Eq4J7ivhk+0cODQJ+AD4OnBRsD8p6IE4BvgW8DZwAnBh6L57cE2Hy3HFk98EYR2KRYsWeV2vmwry5gTM7L+BI5MQNrMpuFaHbOd+Dxya45wB1wdbh2XMmDFeNyW606dPL3kFYqXSG4ecOQFJR+U6F+c6T1uq6QsVpRIsim4SlWlp1k0L+YoD3w86CY2V1KaAI6m3pAuDcnmHK4+ngaamzD5SHVc3yksURTeJl7GannNS5CwOmNnJkkYDXwTulbSNHZOKCPg98D0ze7QsllYZaW2/9rodWzcOhXoMzjez0UBXXAXdZ3CTiXQzs/O8A4hPrbVfe930EmlVYjPbCvwlYVtqimJXBva6Xjcp/NLkJWbJkiW0tLRQV1fH0KFDc15XX19fRqu8bq3oxsFPOV5iWlpaWLx4MS0t+btQLF26tEwWFa+bRFY2zemtJt04eCdQYurq6hg2bBh1dXV5rxs7dmzBuJJ4GdOsmwS1phuHqOsOnBN08/UUYOjQoYwYMSJvUQBg/vz5BeNK4mWMopsEXje9RM0J3Aqsl/Q9SSckaVAaSeJl3LJlS8nj9LpeNw6RnICZHQWcjuvfP0fSi5ImSeqfpHFpIQknMG7cuJLH6XW9bhwi1wmY2bNm9q/A/sCXgfOBlyUtljRekq9fKIL77rvP63rdVFBUOT+YZ/CiYPsQN8x3DXAlbuTfp0ttYLVSqWWrvW5168YhkhOQdCXuxT8YeAC42MyeDp2fA+RcrMTj8aSXqFn4UbhZffuY2RVhBwBgZu/jcwFFsWzZMq/rdVNBlIlGO+Em95ybbxIRM+s4qy0kTJSKxAkTJlSNbhS8bnop6ASCxUgPIM+qpp62RHkZZ82aVTW6UfC66SVqceAbwI8k9ZPUSdIurVuSxlUznTt39rpeNxVEbR34SfD34lCYcBN9Vn71hA7I6NGjva7XTQVRv+QHBNuBoa312BOD2bNne12vmwqiOoHzzWx15obrG+CJwYknnuh1vW4qiOoEbs4RPqlUhtQazc3NXtfrpoJCy5CdKulU3OIgI1qPg+1SdiwN5imS5cuXe12vmwoKVQzeG/zdDfhpKNyAN4CrkjCqFqi1CTC9bnopNNHoAWZ2ADCzdT/YDjSzE80sibUIa4JamwDT66aXqBONfk7SR4CP47oO/ypYKBQzey9JA6uVHj16eN0a033x8hfbHJ+zxzk7hQ28Z2CidmUj6gCiwbgViLcA+wG/Ak4BPg98NjHrqpjhw4d73QrpluNljJLezx6WjlcnauvAD4GbzexQdqxS/ARwciJW1QBz5szxuinRTeJlrFR64xDVCQwCfhHsG2wvBuyehFG1QJq+jF63enTjENUJrALazJIg6XjgpVIbVEokTZTUKKlx48aNlTanDevWrfO6XjcVRHUCNwG/kfQNoE7SDcCDpLyzkJlNN7MGM2vo1atXpc1pw8qVK72u100FUVsH5ks6A7gUVxfQD/i0mXWcmRNSRq21X2fTLUcFXZrSWyzz5s1rc3zCCSe0CTv77LPbrQHFTTT6v8GsQmeZ2Re9A2gftdZ+HUU3iQq6UqR33rx5bbbWlzG8JaGbyahRo0oeJ0RvItwVGAscg5tlaDtm1nG6RqWI3r17e90OqhvlZaxUeuMQdT6BXwCDgcdw3YU97aShocHret1UENUJjAL2NzM/YKhEzJs3ryLTUndk3UJlZNi5nJym9O5/7f5573lt2mtJmpSTqE7gL0B3/KjBkjFy5Eiv206iZMuz6ZbjZazUc45DVCdwEfATSQvIKA6Y2f0lt6oGWLlyJSeddFIqdMtRS5+m9FZKt1Jf+kJEdQITgE8AH8WtR9iKAd4JxGDVqlU7hZXjZcymm0kStfRRdJMgm245XsZKpTcOUZ3AV4BjzOyFJI2pJaK0IyfxMnbkdnOvmwxR+wm8gVtz0FMi0txeX4i0tJtHodZ04xA1J/AfwExJd5Cx5qCZvVJyq2qA/v37tzuOOLXlpdDNJEoFXTbdclTQJZHeNOvGIaoT+H7w95yMcL/uQEwGDBhQ8jijvIxJ6EbB66aXqGMHamqloXJU0C1YsKAitdZp0i1HBV2a0ptWouYEappCFXRxsuWlGvxRLF63unXjEHXsQF9gMtnHDnScfE+ZiJItb2xsrEhPNq/bcXTLtYBJ1JzAg8AK3CIkHxS41hOB9evXe90a080sQk6fPj0VTYlRncChwIlm9mGSxtQStdZ+7XV3ppATWLp0ad77yz2fwDzc7MKeElFr7ddp0h14z8A22xNDntgpLAndtBI1J3A18AdJL7Pz2IF/LrlVNUCamq7K0V6fpvRmUuiLXKhsnu2LXXVNhMAMYBvwAr5OoCT06dPH63rdVBDVCZyKW3nIDyUuEYsWLarItNTZdMvRXl+K9Mb5IpdCt1DZPBuV+n3jENUJPAf0wM8nUDLGjBmzU1g5suXZdMuB100vUZ3A48ACSTPYuU7gp9lv8eRj0aJFDB482OsWQdwvckdNb7mI6gROBtYCmdOlGG2XLPdEpKmpaaewcmTLs+kWS5xseTbdcrSblyK9HUk3DlHHDoxI2pA0UY5seZrbryulm4QTSHN600LUbsPTgPvN7M/JmtMxifNlnD59OlOmTOmQunGy5aXQjUOt6cYhanHgI8DvJG0Efg7MNLPXkzOrspQjW16p8qLXrW7dOEQtDlwl6RrgDGA8MEnSM7j5BR82s+bkTEw/cb6M9fX1hS9KqW4cvG56KWYZsm1mNt/MxgIfB3oB9wEbJP1E0r4J2ViVxHmBva7XTYLITkDSXpK+IGkhsBh4BjcD8WFAM251orIhqbukRyS9J2m1pHHl1G8vY8eO9bpeNxVEcgKSHsI1EX4a+BGu9+BEM1tiZq8B1wIHJGdmVr4PtAB744ooP5Q0qMw2xGb+/Ple1+umgqgVg08DV5rZhmwnzexDSXuXzqz8SNoTGAMcEdRHPCVpLnAx8PVy2dEetmzZ4nW9biqIWjF4V/hY0ghgm5ktDl3zfolty8eAQH9lKOxZOtBw53HjKlN6yaZbjk47aUpvNevGwswKbsATwNBg/2u4rsNrgRuj3F/qDVcXsSEj7DJgUUbYRKARaOzatavhejgaYI2NjdbY2NgmbPLkyWZm1muPXtvDDu95uK2YuMLOP/T8NteuXbvW5s6d2ybsiiuu2CnsuOOOs7lz59ro0aN30rrnnnvahP3gUz+wJ8Y/0Sbs/EPPtxUTV9jhPQ/fHta7d28zM7vwwgvbXDtt2jSbNm1a1jT17t17e9ixxx5rZmaXXXZZwTTdc889Zu5h7pSm4447rk343Llz7YorrtgpbO3atW3CLrvsMjMzO/bYY9ukaciQITZ58uSCv9OFF15oc+fOte7du28PO+igg2zu3Lk7penaa6+NlKbRo0ebme30O5nZTmmaNGmSzZgxI1KazCxSmrL9Tq1pGjlyZJtrZ8yYYZMmTSoqTfneJwUvS14kNQEfM7Ntkl4CzsZVBi4xs74FIygxko4JtPcIhV0HDDezrNOtNDQ0WGNjY6T4M2cWzkb465ltsY1Msq2WmxlWKd1MGhoayPesvG46dIucWUi5TkRtHdgFMEkHATKzF4IKwY8WY0UJWQnsKumQUNhRwPMVssfj6bBErRh8Cvge0Bt4BCBwCG8mZFdezOw9SQ8Dt0i6FDgaOBfoGBO9A8uWLdvJk5ejbJ5NtxyUQjfz/ijPpyOnt1wUsyrxdcBG4M4g7FDgOwnYFJUrcCMY/wY0AV8ys5LkBMrxMk6YMKHgNZXSLUScl7EUupmUSzct6U2KSMUBM2sysxvNbHJrF2Ez+42Z3Z2odfltesvMzjOzPc2sr5nNSkoriUkjZ81KzNyy60Z5Pj696SXyCkSSjsbVyvckVMlgZjeX3qyORZwvRefOnTusbhy8bvHE+X3jELXH4ERgCW6uwa8Bg3HFg4NLblEVEOVLMXr06KrRjYLXbT9JTWMetXXgemCUmf0T8EHw9zPA1kSsqgFmz57tdb1uKojqBD5mZk8G+x9K2sXMHsP1F/DEoFzrzHnd2tKNQ9Q6gdcl9TezVbg2+nMlvYkbwOOJQXNzZaZg8LrVrRuHqDmBqbghwwC3AL/AzUD8jSSMqgWWL1/udb1uKog6gOi+0P5jkj4K1FmNzyjUHmptAsw06yZhW0eaaDTypCJhzKzFO4D2kaYFOmtdN4kXtlK6cYjlBDztp0ePHl63xnXT4gQidxbylJZKrVPnddOhu2bqGrY1b6NTfSf6Xl/2gbhtKJgTkLSLpFMl1ZXDoFphzpw5XreGdbc1b2P1ravZ1rytTBblpqATMLMPgf80M98cWELS+oXyuuXR7VTfiX439aNTfafyGJSHqMWBxZI+bmZPJ2pNDbFu3TqvW8O6lS4ChInqBFYDj0n6T+A13PRFgB9AFJeVK1cWvsjrJq6bVNm8UumNQ1QnsDvw62B/v2RMqS3S3G5eS7qtZfN+N/Urq+5zZz633fkc+eiRJdUulqjzCVySa0vawGolze3mtaT7buO7dP1EV95tfLesuu+/8D6bntzE+y+Uc5Lu7BSzAtFhkm6S9L3geKCkyrqwDkzv3r29bgZJ5BYK6XZp6MKmJzfRpaFLWXXplPG3gkRdmvx84AfAHGAccCXQBbgDOD0x66qYhoYGr5tBEk6gkG5SOYFCun0m9tleHKg0UXMCtwCfNLMvAq0Nm8/iZvj1xCDKNNZeN3ndpHIChXT7Xt+XA245IBWtBFErBj+Ge+lhR8uAhfY9RTJy5EivmwLdpNrrC+l2qB6DActw6/yFuRD4Y2nNqR3S2mRWa7pJfZEL6aapx2DUnMDVwAJJXwD2lPQ73HqAlXHvVcCqVau8bg3rdrgeg2a2QtKhwGhgPq7D0Hz8zEKxSWu7udctj26liwBhos42fLeZvW9mD5jZnWb2S9wko3OTNa96SWu7OVRmfP2aqWt49eZXWTN1TVl1k6JSunGIWidwpKTtU4lJ2h14DLc6sScG/fv3T61uEk6gkG5SZeQ0P+e0ELVO4Fzg95I2AffgHMCLuKW/PTEYMGCA1w2RVBk5relNE1G7Db8LnAFcAvwZeNbMLrMo65p7srJgwQKvGyKpWvq0pjdN5MwJSLolS/AfgbOAt1vP+1GE8ajUirVp1U2q3byQblIDeTrKisSQPyewf5ZtF1xRIBzmiUFjY6PXDZFUnUAU3U1Pbiq77nNnPsefhv2J5858rqS6cciZE6iGEYLBGooTAfr2TU+TDMD69eu9boik6gSi6Hb9RNey67Y6n66f6FpS3ThEHUB0ONBkZm9Iqgf+FTeG4C4zq/xYyByY2XRgOkBDQ0Oq6i/S2n5dKd2k2s0L6SY1lr+QblLOJw5RmwhnAd2C/buAYcCJuJYCTwzS2n7t2+vLo3vko0dyzOJjKj6hCERvIuxvZi9KEvBPwCDgA+DVxCxLEUl8PdPadJXUTDtpTW+ldDviAKItkroAxwOvmdmbwBZgt8QsSxFJOIE+ffqUPM5S6CZVNk9reiulm6YBRMUUBx4HfgbcF4QdS43kBIolitNYtGhRSeIplkK6SbXXR0lvEqRVtyMOIPqqpJHAVjNbGAR/CHw1Mcs6MFFe3jFjxpQknmKJopsEXrctlS4ChIk8x6CZLQg5AMys0cweT8as6ietXyiv27F145Cvx+BvzWxUsP8kOWYRMrNhCdlW1TQ1NXldr5sK8hUH7g/t/yRpQ2qNtLbXe92OrRuHnMUBM5sV2v9Zrq08ZlYfaW2/9rodWzcOkZcmDyoGjwbqw+F+AFE8Bg8e7HW9biqI2m34e8AFwEIgtd2EOxL19fWFL6qAbrGdWJYsWcLYsWNZsmQJQ4cOja2bFLWmG4eorQNjgSFm9lm/DFlpWLp0aSp1i+3E0tLSQnNzMy0t+aebTCK9UcrdaX3OaSJqcaAJeCdBO6qGmTNn0tLSQl1dHePHj8953dixY0uuHeWlKKRbbCeWuro6hg0bRl1dXbt041CK9CZFpXTjkK+J8MDQ4beBmZJuJ2NeQTN7JSHbOiQtLS2sWbOm4NDl+fPnM3DgwJJqR3kpCukW24klXxGgGN2kqDXdOOTLCbyUJWx0xrGRiiUV00NdXR19+/Yt+GXcsmVLmSxKVnfJkiXbcz75HEK1pDftunHIN6lI5N6Enh3kKwKEGTduXMKWlEe3paWFxYsXM2xY/j5j1ZJeiJbjqlR645D3RZd0gaS9y2VMLXHfffdVhW7UOoFqSS9EcwKVSm8cClUMfhM4SNLLwGLgCWCxma1O3LIqZ8iQIVWhG7VOoFrSm3bdOOR1AmY2IMgJDAu264AZktYSOAUz812KPWUnal2EpzAFy/1m9oaZPWhmV5nZ0UBP4PvAJ/HTi8Vm2bJlXrcdtNZFFOqfUC3pTZKC/QSCKcWOZkdu4CRgHfAA8GSSxlUzEyZM8LrtIGpdRLWkN0kKVQzOB9YC03A5gOnAADMbYmZXm9mDZbCxKpk1a1bhi7xuToYOHcqIESMKFgWqJb1JUqg4MBA3l+CrwMvAS8GSZJ520rlzZ6/rdVNBoYrBQzIqBq+R1BNYgisKPGVmf07cyipk9OjMfldetxiiVgxWS3qTJE7F4BHAMmBS8NcTg9mzZ3vddhC1YrDUujNnzmTGjBnMnDmzrLpJEqdi8GTcQiSNwE8TtK2qOfHEE71uO4haMVhq3ahjQ5J4zknNVpTXCUj6Da41oA54BtdZ6HvAUjPbnIhFNUJzc7PXbQdR+waUWjfq2JAknnNSTqBQceBJ3FLk3czsVDP7hpkt9A6g/SxfvtzrdkDd8ePHc8kllxQcI1Kp9MZBZqlapzMxGhoarFLLcmdj3bp1FVkdx+tWt24elOuEHylYIWptAkyv2z6WLFnCwoULWbJkSUnjhQo5AUlHSPqdpDcl7ZQVkdRd0iOS3pO0WtK4jPOnSVoh6X1JCyWVdvXMMtCjRw+v63UjE7U1JA6VyglsxXU7/kKO898HWoC9gfHADyUNAgj6KTwM3AR0x7VS/Cppg0vN8OHDva7XjUzU1pA4VMQJmNmLZnYv8HzmOUl7AmOAm8ys2cyeAuYCFweXfBp4Pui7sBmYAhwl6dDyWF8a5syZ43W9bmSidpOOQxrrBAYA28xsZSjsWWBQsD8oOAbAzN7DdWkeRAeiWr5QXjddunFIoxOoBzZlhG0CukQ8vx1JEyU1SmrcuHFjyQ1tD+vWrfO6XjcVlMUJSBovqTnYHitweTOwV0bYXsC7Ec9vx8ymm1mDmTX06tUrjumJsXLlysIXeV2vWwYq2k9A0sHAX81MobA9gbeBQWb21yDsfmCdmX1d0kTg82Y2NHT9RuBYM1uRS8v3E/C6taCbh3T1E5BjN1x3ZCTtJqkzbC/jPwzcImlPSUOBc4GfB7c/AhwhaUwQx83Ac/kcQBqplvZrr5su3ThUqk6gH/ABO1oHPgBeDJ2/Atgd+BswG/iSmT0PYGYbca0H38LlGE4ALiyP2aWjd+/eXtfrpoLIqxKXEjNbRZ7siZm9BZyX5/zvgQ7VJJhJQ0OD1/W6qSCNrQM1wbx587yu100FNTOASNJGINt6CT2BN8tsjtf1uuXmTTMble1EzTiBXEhqNLOy5928rtdNC7444PHUON4JeDw1jncCbi0Fr+t1q0W3aGq+TsDjqXV8TsDjqXG8E/B4apyacQKSHpdkknYNhVXVNGaSPi9pmaS/S3pd0tRqTm8+CqXVs4OacAKSxpO9i3S1TWO2B3ANrqPKCcBpwL+EzldbevORM62eDMysqjegK7AS+DhgwK5B+J64f5IBoWt/DtwR7E8E/hA6tyduoNOhlU5TEWm/FphXK+nNsD1nWv3WdquFnMBtwA+BDRnhtTCN2TB2jNSshfS2UiitnhBV7QQkNQBDge9mOV2yaczSiKRLgAbgriCoqtObQTWlJXGqyglkmcbsB8BXzOwfWS4v2TRmlSLXtG2SzgPuAM4ws9ZBLB0+vUVQTWlJnKpyAmY208zqzaweGIv7Ev5K0gbgf4LLXpf0CVw9wa6SDglFcRQ7ss/PB8fA9mnMDiLLNOmVIpxeMzsDQNIo4MfA2WYWXhCvw6e3CAql1ROm0pUSSW24SUv2CW3H4SoG9wXqgmt+iZu5aE9csWETbm5DgF7B8RhgN+Dfgacrna4CaT4VaAKG5ThfVekt8CxyptVvGc+q0gaU8Z+iP6HWgSCsO/Br4D1gDTAu457TgRW4WvJFQP9Kp6NAGhcC/8Blh1u3x6o1vQWeRd60+m3H5scOeDw1TlXVCXg8nuLxTsDjqXG8E/B4ahzvBDyeGsc7AY+nxvFOwOOpcbwTSAhJz0saXmk7MpE0QdJTlbaj0khaJen0mPd2lvQXSfuU2KarJd1Ryjij4J1ATEJ99pslfSjpg9DxeDMbZGaLKm1ne/AOIycTgcVmljkytb1MBy6S9LESx5sX7wRiYjv67NfjeqSdHQqbWSm7wjMJeRLjcnaskl0yzGwz8BjwuVLHnQ/vBBIinN2UNEXSg5J+IeldScslDZB0g6S/SXpN0sjQvV0l3StpvaS1kr4pqVMOnSmSHgri/jswoZj789h/GPAj4MQgd/NOyLb7JW0Mpu2aJCnr/5Gk4yU1BtOdvSFpWujcg5I2SNokaXF41h9Ju0v6dhD/JklPSdo9OHeypD9Ieid4bhOC8LMk/SnQek3SlAxbLg7ia5L0bxnndpH0dUkvB+cfkNQ9R5r64gZWPRMKu0/SDyQ9FjyrJZL2kXS3pLflpmw7JnT914Lf5V1JL0o6LSSxCDgr9y9TerwTKB9n474eHwX+BPwO9/z3BW4B7gld+zPcGICDgWOAkcCleeI+F3gI6AbMjHH/TpjZC8AXgaVB7qZbcOq7uNmaDgROwX21LskRzXeA75jZXrgX54HQuceAQ4CPAf8b2N3KXcAQ4CTcGIDrgQ+DF/CxwIZewNHAn4N73gts6YZ7ib4UDKlG0uG4iWUuBvoAPYD9QnpX41bBPiU4/zZuerJsDAZesZ2Hp18ATMJN7bYFWBqkqyfut5kW2DIQuBI4zsy6AJ8CVoXieYHQaM6yUOnBC9WwBT/i6bnCgCnAf4XOnY0b3NMpOO6CG9zUDTcn3hZg99D1Y4GFObSn4Mqnrcd57wcmAE9FTFeba4FOQdyHh8IuBxbluH8x8A2gZwGdbkH6u+Ic4wfAUVmuuwF4JKLtdwP/EezfDPwydK51+rHW3+cF4LTQ+d7AVkKDzULnxpMxuhK4D/hx6Pgq4IXQ8WDgnWD/YOBvuMFaH8kS/yG4WZHK9v/rcwLl443Q/ge4VWK3hY7BzYjTD/gIsD7I8r6DyyXkqyx6LbQf5/6o9ATqaLu682pcbiYbX8BN9bVC0v9IGg0gqZOkO4Ls99/Z8SXsGWy74aY2y2T/HOFIOkFuhuSNkjbhcjE9g9N9CD0jc1OnNYVu7wc8EnpeLwDbcA41k7fJPkNR5u+beVwfaL+Emwx2CvA3Sb+U1Cd0bRd2nhUpUbwTSB+v4b62Pc2sW7DtZWb55scLDwWNc3+UeMEttb0V99K00hdYm/Vms7+a2VicA/p34CG5yUrG4Yowp+O+/v2DWxRobMYVHzJ5LUc4wCxgLrC/mXXF1WcoOLce50CciLQHrkgQjveM0PPqZma7mVm2dD0HHNieClgzm2VmJ+Oeo+GeTSuHEZrrsRx4J5AyzGw9sAD4tqS9gkqrgySdUo77M3gD2E9SXRD3Nly5/luSusitS3At8ItsN0u6SFIvM/sQeCcI3ob72m3BfY33wE0G22r/h8BPgWmS+gS5hhMldcbVG5wu6QJJu0rqIeno4NYuwFtmtlnS8ThH08pDwOigUrEOVwcT/t//UZCmfoHdvSSdmy1NZvY68Ffg+ALPLiuSBko6NUjPZlwuYVvoklNw9R5lwzuBdPI5XLb7L7js50O4cmrJ75fr1DQ+RzyP46bk2iCpda7Cq3CVcK8AT+G+wD/Ncf8o4HlJzbhKwgvNNYPdjytGrA1sfDrjvn8BluOmhHsL96XcxczWAGcC1wXhf2ZHJdoVwC2S3sXVAWyvhDSz54EvB7auD57J6yG97+ByEQuC+5/GrduQi3twlYxx6Iyb//FN3AzYHwNuBJC0W5C+n8WMOxZ+UhGPp0iCr/ifcJWJ60sY71W44sz1pYozkq53Ah5PbeOLAx5PjeOdgMdT43gn4PHUON4JeDw1jncCHk+N452Ax1PjeCfg8dQ43gl4PDXO/wfZv49KQjcCRAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
GroupDataShuffle
RangeLeftRightPIQR_minQ1Q2Q3IQR_maxOutliersIQR_minQ1Q2Q3IQR_maxOutliers
0-400 – -300 ms120942.788900e-01-0.530916-0.0493020.1224880.3059170.778465[]-0.561890-0.0771810.1048670.2814170.698002[-0.7038,-0.7033,0.8216]
1-300 – -200 ms120941.243208e-01-0.560208-0.0339210.1573500.3723990.811325[]-0.693004-0.0845410.1063820.3309920.693507[]
2-200 – -100 ms120944.421350e-02-0.4020240.0207610.2004680.3792750.829799[-0.6771,-0.5341]-0.658452-0.0571550.1336070.3472310.913788[-0.7372,-0.6775]
3-100 – +0 ms120949.214959e-06-0.3119360.1386470.3115170.4741350.908852[-0.4840,-0.4222,0.9837]-0.681709-0.0410690.1799100.4099111.011172[-0.8646,-0.7956]
4+0 – +100 ms120945.216176e-08-0.1687440.2365230.3860470.5422400.919303[-0.4836,-0.4795,-0.3390,-0.3100,-0.2868,1.061...-0.738590-0.0530440.2016860.4564451.093199[-0.8822,-0.8696]
5+100 – +200 ms120941.409669e-06-0.1830120.2556200.4031900.5578781.011172[-0.5468,-0.4975,-0.4861,-0.3602,-0.3033,-0.29...-0.868905-0.0766790.2628970.4771931.149776[-0.9441]
\n", "
" ], "text/plain": [ " Group Data \\\n", " Range Left Right P IQR_min Q1 Q2 \n", "0 -400 – -300 ms 120 94 2.788900e-01 -0.530916 -0.049302 0.122488 \n", "1 -300 – -200 ms 120 94 1.243208e-01 -0.560208 -0.033921 0.157350 \n", "2 -200 – -100 ms 120 94 4.421350e-02 -0.402024 0.020761 0.200468 \n", "3 -100 – +0 ms 120 94 9.214959e-06 -0.311936 0.138647 0.311517 \n", "4 +0 – +100 ms 120 94 5.216176e-08 -0.168744 0.236523 0.386047 \n", "5 +100 – +200 ms 120 94 1.409669e-06 -0.183012 0.255620 0.403190 \n", "\n", " \\\n", " Q3 IQR_max Outliers \n", "0 0.305917 0.778465 [] \n", "1 0.372399 0.811325 [] \n", "2 0.379275 0.829799 [-0.6771,-0.5341] \n", "3 0.474135 0.908852 [-0.4840,-0.4222,0.9837] \n", "4 0.542240 0.919303 [-0.4836,-0.4795,-0.3390,-0.3100,-0.2868,1.061... \n", "5 0.557878 1.011172 [-0.5468,-0.4975,-0.4861,-0.3602,-0.3033,-0.29... \n", "\n", " Shuffle \n", " IQR_min Q1 Q2 Q3 IQR_max Outliers \n", "0 -0.561890 -0.077181 0.104867 0.281417 0.698002 [-0.7038,-0.7033,0.8216] \n", "1 -0.693004 -0.084541 0.106382 0.330992 0.693507 [] \n", "2 -0.658452 -0.057155 0.133607 0.347231 0.913788 [-0.7372,-0.6775] \n", "3 -0.681709 -0.041069 0.179910 0.409911 1.011172 [-0.8646,-0.7956] \n", "4 -0.738590 -0.053044 0.201686 0.456445 1.093199 [-0.8822,-0.8696] \n", "5 -0.868905 -0.076679 0.262897 0.477193 1.149776 [-0.9441] " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grouping = [\n", " ('-400 – -300 ms', slice(40, 60)),\n", " ('-300 – -200 ms', slice(60, 80)),\n", " ('-200 – -100 ms', slice(80, 100)),\n", " ('-100 – +0 ms', slice(100, 120)),\n", " (' +0 – +100 ms', slice(120, 140)),\n", " ('+100 – +200 ms', slice(140, 160))\n", "]\n", "\n", "fig = plt.figure(figsize=(3.5, 3))\n", "ax = plt.gca()\n", "sep = np.arange(1, len(grouping)) - 0.5\n", "\n", "tab = {\n", " ('Group', 'Range'): [],\n", " ('Group', 'Left'): [],\n", " ('Group', 'Right'): [],\n", " ('Group', 'P'): [],\n", " ('Data', 'IQR_min'): [],\n", " ('Data', 'Q1'): [],\n", " ('Data', 'Q2'): [],\n", " ('Data', 'Q3'): [],\n", " ('Data', 'IQR_max'): [],\n", " ('Data', 'Outliers'): [],\n", " ('Shuffle', 'IQR_min'): [],\n", " ('Shuffle', 'Q1'): [],\n", " ('Shuffle', 'Q2'): [],\n", " ('Shuffle', 'Q3'): [],\n", " ('Shuffle', 'IQR_max'): [],\n", " ('Shuffle', 'Outliers'): [],\n", "}\n", "\n", "for base, (group_name, rng) in enumerate(grouping):\n", " baseoffset = base\n", " tab['Group', 'Range'].append(group_name)\n", " tab['Group', 'Left'].append(epochs['left']['N'])\n", " tab['Group', 'Right'].append(epochs['right']['N'])\n", " \n", " _data = np.median(dmerged[rng], axis=0)\n", " _rand = np.median(rmerged[rng], axis=0)\n", " for shift, lab, vals, col in ((-1, 'Data', _data, 'm'),\n", " (1, 'Shuffle', _rand, 'gray')):\n", " box = BoxSummary.from_values(vals)\n", " plot_box_around(baseoffset + 0.2*shift, box, marker='x', color=col, ocolor=col, boxwidth=.3, elinewidth=.8)\n", " tab[lab, 'Q1'].append(box.q1)\n", " tab[lab, 'Q2'].append(box.q2)\n", " tab[lab, 'Q3'].append(box.q3)\n", " tab[lab, 'IQR_min'].append(box.min)\n", " tab[lab, 'IQR_max'].append(box.max)\n", " tab[lab, 'Outliers'].append('[' + ','.join(f\"{v:.4f}\" for v in sorted(box.outliers)) + ']')\n", " \n", " p = sstats.mannwhitneyu(_data, _rand).pvalue\n", " plt.text(baseoffset, 1.15, plabel(p), fontsize=10, ha='center', va='bottom')\n", " tab['Group', 'P'].append(p)\n", "\n", "plt.vlines(sep, -2, 2, color='k', linewidth=.5, linestyle='dashed')\n", "plt.hlines(0, -0.7, 6.7, color='k', linewidth=1, linestyle='dashed')\n", "plt.xlim(-0.7, 5.7)\n", "plt.ylim(-1.1, 1.35)\n", "plt.xticks((-0.5, 0.5, 1.5, 2.5, 3.5, 4.5), (\"-400\", \"\", \"-200\", \"\", \"0\", \"\"))\n", "plt.yticks((-1, -0.5, 0, 0.5, 1), (\"-100\", \"\", \"0\", \"\", \"100\",))\n", "plt.xlabel(\"Time rel. to saccade (ms)\", fontsize=12)\n", "plt.ylabel(\"Whisker asymmetry (%)\", fontsize=12)\n", "plt.title(f\"N={epochs['left']['N']} leftward, {epochs['right']['N']} rightward\", fontsize=9)\n", "for side in (\"top\", \"right\"):\n", " ax.spines[side].set_visible(False)\n", "ax.tick_params(labelsize=12)\n", "plt.subplots_adjust(left=.22, right=.98, bottom=.18, top=.98)\n", "if saved == True:\n", " figpath = figdir / \"saccade_binned.png\"\n", " if not figdir.exists():\n", " figdir.mkdir(parents=True)\n", " fig.savefig(str(figpath), dpi=300)\n", "plt.show()\n", "\n", "tab = pd.DataFrame(data=tab)\n", "if saved == True:\n", " tabpath = figdir / \"saccade_binned.csv\"\n", " tab.to_csv(str(tabpath), header=True, index=False)\n", "tab" ] }, { "cell_type": "code", "execution_count": null, "id": "6f78d945", "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": 5 }