{ "cells": [ { "cell_type": "markdown", "id": "d0e5f760", "metadata": {}, "source": [ "# Missing data\n", "\n", "Pupil and gaze data have gaps and bad samples for many reasons: the eye\n", "tracker loses the pupil during a blink, the gaze leaves the screen, the\n", "two eyes point to different places, or a part of the recording is simply\n", "not part of the analysis. `pypillometry` keeps track of *why* a sample is\n", "bad, and a *mask policy* decides which samples the analyses leave out.\n", "\n", "This tutorial goes through these steps with a real recording. The\n", "concepts are explained in the user guide ([Missing data](../user_guide/missing_data.rst))." ] }, { "cell_type": "code", "execution_count": 1, "id": "76a76342", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:35.084963Z", "iopub.status.busy": "2026-11-09T20:27:35.084769Z", "iopub.status.idle": "2026-11-09T20:27:36.189982Z", "shell.execute_reply": "2026-11-09T20:27:36.188736Z" } }, "outputs": [], "source": [ "import pypillometry as pp\n", "from pypillometry import u\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "pp.logging_set_level(\"WARNING\")" ] }, { "cell_type": "markdown", "id": "7771a5b3", "metadata": {}, "source": [ "We use the example dataset `rlmw_002`: a 20-minute recording of both eyes\n", "(pupil and gaze, 500 Hz) from a reinforcement-learning task. It is\n", "downloaded from OSF on first use (see [Example datasets](example_datasets.ipynb))." ] }, { "cell_type": "code", "execution_count": 2, "id": "44cb66d8", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:36.194270Z", "iopub.status.busy": "2026-11-09T20:27:36.193763Z", "iopub.status.idle": "2026-11-09T20:27:38.291215Z", "shell.execute_reply": "2026-11-09T20:27:38.289626Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "

EyeData — rlmw_002

\n", "\n", "
\n", "
\n", "
Samples
\n", "
597,301
\n", "
\n", "
\n", "
Sampling rate
\n", "
500.0 Hz
\n", "
\n", "
\n", "
Duration
\n", "
19.91 min (start 0.00, end 19.91)
\n", "
\n", "
\n", "
Eyes
\n", "
left, right
\n", "
\n", "
\n", "
Variables
\n", "
left_x, left_y, left_pupil, right_x, right_y, right_pupil
\n", "
\n", "
\n", "
Events
\n", "
1320
\n", "
\n", "
\n", "
Size
\n", "
39.9MiB
\n", "
\n", "
\n", "\n", "\n", "
\n", " Signal overview\n", "
\n", "
\n", "
left_pupil
\n", " \"left_pupil\n", "
\n", "
\n", "
left_x
\n", " \"left_x\n", "
\n", "
\n", "
left_y
\n", " \"left_y\n", "
\n", "
\n", "
right_pupil
\n", " \"right_pupil\n", "
\n", "
\n", "
right_x
\n", " \"right_x\n", "
\n", "
\n", "
right_y
\n", " \"right_y\n", "
\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", " Data glimpse\n", " 
EyeDataDict(vars=6,n=597301,shape=(597301,)): \n",
       "  left_x (float64): 655.6, 655.7, 655.0, 654.5, 655.0...\n",
       "  left_y (float64): 599.9, 598.9, 597.6, 597.8, 597.8...\n",
       "  left_pupil (float64): 1121.0, 1122.0, 1124.0, 1126.0, 1126.0...\n",
       "  right_x (float64): 773.1, 773.8, 774.8, 776.3, 776.7...\n",
       "  right_y (float64): 594.2, 593.9, 596.5, 597.4, 597.3...\n",
       "  right_pupil (float64): 949.0, 951.0, 956.0, 959.0, 960.0...\n",
       "
\n", "
\n", "\n", "
\n", " Recent history\n", "
    \n", "
  1. reset_time()
    2. \n", "
\n", "
\n", "
\n" ], "text/plain": [ "EyeData(rlmw_002, 39.9MiB):\n", " n : 597301\n", " sampling_rate : 500.0\n", " eyes : ['left', 'right']\n", " data : ['left_x', 'left_y', 'left_pupil', 'right_x', 'right_y', 'right_pupil']\n", " nevents : 1320\n", " experimental_setup: {'screen_resolution': (1280, 1024), 'physical_screen_size_mm': (300.0, 200.0), 'eye_to_screen_distance_mm': 600.0, 'camera_position': 'not set'}\n", " calibration : not available\n", " duration_minutes : 19.91\n", " start_min : 0.0\n", " end_min : 19.91\n", " params : {}\n", " nblinks : {}\n", " blinks : {}\n", " glimpse : EyeDataDict(vars=6,n=597301,shape=(597301,)): \n", " left_x (float64): 655.6, 655.7, 655.0, 654.5, 655.0...\n", " left_y (float64): 599.9, 598.9, 597.6, 597.8, 597.8...\n", " left_pupil (float64): 1121.0, 1122.0, 1124.0, 1126.0, 1126.0...\n", " right_x (float64): 773.1, 773.8, 774.8, 776.3, 776.7...\n", " right_y (float64): 594.2, 593.9, 596.5, 597.4, 597.3...\n", " right_pupil (float64): 949.0, 951.0, 956.0, 959.0, 960.0...\n", "\n", " History:\n", " *\n", " └ reset_time()" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = pp.get_example_data(\"rlmw_002\", cache_dir=\"example_data\")\n", "d" ] }, { "cell_type": "markdown", "id": "1271294a", "metadata": {}, "source": [ "## Missing samples in the raw data\n", "\n", "Samples that the eye tracker did not record are NaN (Eyelink codes them as\n", "0; they are converted to NaN when the data are read). NaN is the only kind\n", "of missing data that is stored in the values themselves.\n", "`d.mask(eye, variable)` returns `True` for each unusable sample:" ] }, { "cell_type": "code", "execution_count": 3, "id": "77f21c90", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:38.297064Z", "iopub.status.busy": "2026-11-09T20:27:38.296457Z", "iopub.status.idle": "2026-11-09T20:27:38.325531Z", "shell.execute_reply": "2026-11-09T20:27:38.323990Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "left_x 3.13% NaN, 3.13% masked\n", "left_y 3.13% NaN, 3.13% masked\n", "left_pupil 3.13% NaN, 3.13% masked\n", "right_x 2.07% NaN, 2.07% masked\n", "right_y 2.07% NaN, 2.07% masked\n", "right_pupil 2.07% NaN, 2.07% masked\n" ] } ], "source": [ "for key in d.data:\n", " print(f\"{key:12s} {np.isnan(d.data[key]).mean():6.2%} NaN, \"\n", " f\"{d.mask(key).mean():6.2%} masked\")" ] }, { "cell_type": "markdown", "id": "d19aa190", "metadata": {}, "source": [ "## Marking bad samples: mask layers\n", "\n", "The detection methods do not change the values. They mark the samples\n", "they find in a *mask layer* of the channels they check:\n", "\n", "- `d.pupil.blinks_detect()` sets the `\"blink\"` layer of the pupil;\n", "- `d.gaze.mask_offscreen()` sets the `\"offscreen\"` layer of x and y (it\n", " needs the screen size, which is stored with this dataset);\n", "- `d.gaze.mask_divergences()` sets the `\"divergence\"` layer of x and y of\n", " both eyes." ] }, { "cell_type": "code", "execution_count": 4, "id": "1abc0353", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:38.329359Z", "iopub.status.busy": "2026-11-09T20:27:38.328952Z", "iopub.status.idle": "2026-11-09T20:27:39.512131Z", "shell.execute_reply": "2026-11-09T20:27:39.511058Z" } }, "outputs": [ { "data": { "text/plain": [ "['blink', 'offscreen', 'divergence']" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.pupil.blinks_detect()\n", "d.gaze.mask_offscreen()\n", "d.gaze.mask_divergences()\n", "d.data.layer_names" ] }, { "cell_type": "markdown", "id": "2baf1d88", "metadata": {}, "source": [ "To see what is masked where, we write a small helper that shows, for each\n", "channel, the percentage of samples that are NaN, the percentage in each\n", "layer, and the percentage that is masked under the `\"default\"` and the\n", "`\"strict\"` policy (explained below). A sample can be in several layers at\n", "once (a blink is usually also NaN), so the columns do not add up." ] }, { "cell_type": "code", "execution_count": 5, "id": "650e0706", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:39.516479Z", "iopub.status.busy": "2026-11-09T20:27:39.515850Z", "iopub.status.idle": "2026-11-09T20:27:39.574354Z", "shell.execute_reply": "2026-11-09T20:27:39.573188Z" } }, "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", "
NaNblinkoffscreendivergencemasked (default)masked (strict)
left_x3.130.000.970.974.864.86
left_y3.130.000.970.974.864.86
left_pupil3.135.840.000.005.845.84
right_x2.070.001.040.974.064.06
right_y2.070.001.040.974.064.06
right_pupil2.074.080.000.004.094.09
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence masked (default) \\\n", "left_x 3.13 0.00 0.97 0.97 4.86 \n", "left_y 3.13 0.00 0.97 0.97 4.86 \n", "left_pupil 3.13 5.84 0.00 0.00 5.84 \n", "right_x 2.07 0.00 1.04 0.97 4.06 \n", "right_y 2.07 0.00 1.04 0.97 4.06 \n", "right_pupil 2.07 4.08 0.00 0.00 4.09 \n", "\n", " masked (strict) \n", "left_x 4.86 \n", "left_y 4.86 \n", "left_pupil 5.84 \n", "right_x 4.06 \n", "right_y 4.06 \n", "right_pupil 4.09 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def layer_table(d):\n", " rows = {}\n", " for key in d.data:\n", " row = {\"NaN\": np.isnan(d.data[key]).mean()}\n", " for name in d.data.layer_names:\n", " if name == \"excluded\": # one layer for the whole recording\n", " row[name] = d.data.get_rec_layer(name).mean()\n", " else:\n", " row[name] = d.data.get_layer(key, name).mean()\n", " row[\"masked (default)\"] = d.mask(key).mean()\n", " row[\"masked (strict)\"] = d.mask(key, policy=\"strict\").mean()\n", " rows[key] = row\n", " return (100 * pd.DataFrame(rows).T).round(2)\n", "\n", "layer_table(d)" ] }, { "cell_type": "markdown", "id": "246d6b6b", "metadata": {}, "source": [ "Each layer belongs to the channels it was detected on: blinks are on the\n", "pupil only, off-screen samples and divergences on x and y only. The gaze\n", "channels are still masked during blinks, because the eye tracker records\n", "no gaze then (NaN).\n", "\n", "The layers can be read back as `Intervals` (in ms), one per channel:" ] }, { "cell_type": "code", "execution_count": 6, "id": "44c9633a", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:39.577669Z", "iopub.status.busy": "2026-11-09T20:27:39.577366Z", "iopub.status.idle": "2026-11-09T20:27:39.612499Z", "shell.execute_reply": "2026-11-09T20:27:39.610849Z" } }, "outputs": [ { "data": { "text/plain": [ "Intervals 'left_pupil_blinks' | 321 intervals (total duration: 69762.00 ms), 217.33 ms +/- 136.83 ms, [22.00 ms, 984.00 ms]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "blinks = d.get_blinks()\n", "blinks[\"left_pupil\"]" ] }, { "cell_type": "code", "execution_count": 7, "id": "f69d6605", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:39.615743Z", "iopub.status.busy": "2026-11-09T20:27:39.615364Z", "iopub.status.idle": "2026-11-09T20:27:39.642724Z", "shell.execute_reply": "2026-11-09T20:27:39.641648Z" } }, "outputs": [ { "data": { "text/plain": [ "{'left_x': Intervals 'left_x_offscreen' | 383 intervals (total duration: 11582.00 ms), 30.24 ms +/- 79.67 ms, [2.00 ms, 544.00 ms],\n", " 'right_x': Intervals 'right_x_offscreen' | 405 intervals (total duration: 12448.00 ms), 30.74 ms +/- 77.91 ms, [2.00 ms, 542.00 ms]}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.get_layer_intervals(\"offscreen\", variables=\"x\")" ] }, { "cell_type": "markdown", "id": "3dc464e9", "metadata": {}, "source": [ "With these, we can draw where each layer is set in each channel, here for\n", "30 seconds of the recording (the grey bars are NaN samples).\n", "`d.plot.plot_masked()` gives a similar overview of the masked stretches\n", "for the whole recording." ] }, { "cell_type": "code", "execution_count": 8, "id": "44e5d7f5", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:39.645562Z", "iopub.status.busy": "2026-11-09T20:27:39.645334Z", "iopub.status.idle": "2026-11-09T20:27:40.632206Z", "shell.execute_reply": "2026-11-09T20:27:40.631234Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/gAAAEmCAYAAADWcDvhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAASS5JREFUeJzt3Ql4U1X6+PG3UFpooaWUshWoLFX2XRGUAVwQx6qAsgsuCAoCoqKCuwMOCOLgwqLgoAIKyOjgDxVBWQQFlVU2RRRZCkKp0EILLYX8n/fwTyYpTRdIs9x8P8+TJ8nNzb0n996c5D1riM1mswkAAAAAAAhoJXydAAAAAAAAcOkI8AEAAAAAsAACfAAAAAAALIAAHwAAAAAACyDABwAAAADAAgjwAQAAAACwAAJ8AAAAAAAsINTXCQhW586dk4MHD0q5cuUkJCTE18kBAAAA4CM2m01OnDgh1apVkxIlqIPFxSPA9xEN7mvUqOGr3QMAAADwM/v375fq1av7OhkIYAT4PqI19/YvcVRUlK+SAQAAAMDH0tPTTeWfPUYALhYBvo/Ym+VrcE+ADwAAAICuu7hUdPAAAAAAAMACCPABAAAAALAAAnwAAAAAACyAPvgAAAAAEADTbGdnZ/s6GfCBsLCwQk+fSIAPAAAAAH5MA/s9e/aYIB/Bp0SJElKrVi0T6BeEAB8AAAAA/JTNZpNDhw5JyZIlzVR6ha3JhTVooc7BgwfNNVCzZs0CZ1ogwAf8zM6dO2XBggUuy4YPHy4xMTESbJKTk2XmzJkuy1q2bClJSUk+SxMAoPicPHlSJk2alO86ffr0kcTERE6Dl6xevVqWL1/OOfChnJwcyczMlGrVqklERIQvkwIfiYuLM0G+XgulSpXKd12KfwAAAADAT509e9bcF6Z5NqzJfu7t10J+CPABAAAAwM8V1DQb1lWUc0+ADwAAAACABYTYdNQGeF16erpER0dLWlqaREVFcQYAAACAIJVfbHD69Gkzgr6Ool66dGnH8uTjp+RYhvemzYuJDJP48mUKvf4999xj+ozPmTPH7TpDhw6VAwcOyH//+988n3tqP4HO3TWQFwbZAwAAAIAAosH9da+slKwc702bFx5aQpaP7FCkIB/eR4APAAAAAAFEa+69Gdwr3Z/utzgD/DfffLPYth0s6IMPAAAAACiWqS8HDRok8fHxpuvBwIEDTXNzd7SJfpcuXVya33ft2tUsr1GjhlSoUEH69+8vGRkZbrexbt06iY2NlQkTJkgwIsAHAAAAAHjcokWLpGzZsrJ161ZZuXKlfP311/LUU08VaRvaH79mzZqybds2WbNmjdnOpEmT8lx3yZIlctNNN5ng/oknnpBg5JEA/88//5SjR48W6T0HDx6U1NRUCVS5P/PFHAMAAAAAsKrq1avLxIkTTc17ixYt5J///KdMmzZNTp06VehtXHnllSZY10EIGzRoIL1795bVq1dfsN68efOkR48e8u6778qAAQMkWHkkwH/wwQflhRdeKNJ7tGmFnuz8HDp0yG+D5tyf+WKOAQAAAABYVfPmzaVkyZIuwbp9RPjCSkxMdHkeExMjx44dc1m2YsUKueuuu2ThwoWmSX8w80iAX7VqVYmLixNP0z4aY8eOFX9UXJ8ZAAAAAKwgJCTEK9uoV6+e1K5dW2bOnGmmzAtmoUVtVh8eHm4GLThy5Ig5eNWqVZPnn39eQkMv3NThw4dNn4vIyEjzuESJEnkGxSdOnJBz586ZZhd2uv3MzEwzF+Tu3bvNsoSEBClVqlSh0qfvO3v2rGkO4mz//v1Srlw5KV++vEs6ndNWmO24+8wAAAAAAJFNmzaZWMpei//jjz+aOEvnc/d05evs2bOlQ4cO0rdvX5k7d27QxmolitqsXkdBbNWqlTRp0kQef/zxPJuna8Dcrl070+dCg/LrrrtO+vTpI6NHj3bZngbubdq0Mc0uKlWqJH//+98d/TF04ITvv/9ePvnkE+ncubO5aXBemPS1bNlS6tata7Z5xx13mIICu1tvvdWU7Dh7+OGHXdJWmO3QJB8AAAAA3NP4TWPGv/76SzZu3GgG2NM4qkwZz0+1pxXPK1asMPvRIF8LFoJRkZvof/rpp/Lyyy+bQeW0ZCQvI0aMkKysLBPoax/6bt26yfLlyy9Y77PPPjMjHOq2/vjjD1PCM336dPOa7qNjx45magQtCNCbNrsoiBYIPP3005KSkiK///672ebFNPP31Hbs9Hikp6e73AAAAADAqm6//XbTWrtRo0bSvn17E9+NGzeu2PYXHx9vRtkP5iC/yO0WtAb8+uuvd/u6Bq4LFiyQxYsXO5q167yFb7zxxgXraq241vTbm1XolAZ6Mi6FNsvQAgWl0ynoiItaUqQjNvpiO3Z6Ib/44osX9V4AAAAAsIuJDJPw0BKSlXPOawdF96f7LSwdzb4gb775Zr7P89rGqFGjzM3dOvHx8fLrr79KsCpygF+nTp18X9cREbU/fcOGDV2W65QGuWkTfmfaX//AgQNyKbR0yFnjxo3NKIvaLCR3P3pvbMdOuwA8+uijLgUhNWrUKPJ2AAAAAAS3+PJlZPnIDnIsI9tr+9TgXvcLiwX4ztMc5KV06dKOJunOdDoEb8i9H/tze7ryGoUxr5EWC9pOUelgEnoDAAAAgEulwTYBN4plmjxnOiJiVFSUfPPNNy7B8YYNG4q8LQ2IizrNwerVq8Vmszmer1q1ygziFxERYZ5XrFjR9Pm309YGW7duLfJ2AAAAAADwJx6fOyAsLMw0Rdcm6TrtnfZff+WVV0zT9qLOg3jFFVfI559/Lps3bzbN9wuaJk9pfwsdmVFHwdf36Wj8U6dOdbx+8803m/7w2sde+/1PmzZNfvvtN8dYAIXdDgAAAAAAARvg64AFWgOemwbKzvPbP/vss6Yp/0svvWTmnL/llltMLb5z8/a8tqXbcB7pUAsKdE76u+++WzIyMmTp0qUFjqQ/bNgwM+2C3mvt/6uvvmre7/y6juSoA95FRkZK165dZciQIaYAoSjbyf2Zcz8HLsWeoxmSkXW+9cqxzGypHhMhtSpGSrAfC7vI8NCgPR7AxXx/+M4g0K5bpfm8Pj5wLFNiIs4P7MW17LtzYj8Pwf6/BPB3ITbndugeokG6c1/9M2fOmMH5dLRDDaaLyw033CCtWrWS8ePH+8V28qOD7GkLh7S0NNOlAXD+Ee34ysoLDsiKkR2C7sfU3bEI1uMBXMr3h+8MAu26nT3gKun3zg8XrMO17LtzwnkQn8QGWlGqA5lrV+iLHQ8Mga0o14DHm+irr776Sr788kvp3r27CfYnT54sp06dkp49e170NrWQYO/evW5f1+b7gBXkrq0uaLmV5feZg/F4AEWR+zvCdwaBwPk6TT2Z9+jgXMvexf8SILAUS4DfqVMn2bVrl+mHr03rmzZtaua3j42NvehtHjp0SDp37uz29eXLl7vtQlBUntoOAAAAAAABHeDrYHrad11vnqKD9e3evTvfdd577z2P7MtT2wEAAAAAIGCnyQMAAAAAoLDdu/v3729mO1u7dq3p4v3666+b7t7aglunNYePa/ABAAAAAMXo+H6RzFTvHeKIWJHyNTy6SR1jTWdc08HN69evbwZmf+edd2TChAnyr3/9y8zIVtSp1oMdAT4AAAAABFpw/2ZLkZws7+0zNFxk6AaPBvkbNmwwQfwjjzziWLZu3Tozq5nW4KPoCPABAAAAIJBozb03g3ul+9P9FiHA//77702N/OHDhyUxMVGGDx9uxlZTb731lkydOlVOnDjhMpj6li1bJDQ01Cxr27atPPfcc3Ls2DGZMmWKeU0Hbu/du7e0b9/e8Z7ffvtNpk+fLr/++qtcfvnl8uijj0qVKlXMa5MmTTLN/OvWrSsfffSRlClTxqRJLVmyRObPn2/S0LBhQ3n44YelQoUKLp+hoHXuvPNOGTx4sEnbDz/8YKY5HDp0qDRp0kR8gT74AAAAAACP0lnO2rVrJ2XLlpU+ffqYAdNbtGhhZkdTGqDfdNNNEhkZKSNGjHDcrrjiChNI6+Nbb73VrKv3K1euNLX6V111lYwdO1ZWrFhhXtN++82aNZPk5GQT+FeuXFm6dOniSMfWrVvN+m+88Ybcdtttcs8995jl48aNk0GDBpk06ft+//13s520tDTHewuzjo4hoOk6ePCgY5r4a665Rv7880+fXFHU4AMAAAAAPOqJJ56Q++67T1599VXzvEePHmb69JdeeknefPNNqVevnqnlDgsLc6nBf/fdd02hgH3Z6dOn5dtvvzXTrjdv3tws0+2mp6ebx9oqQAsAPvjgA8c2Bg4c6JKWiIgI+fzzz6V06dLm+YEDB0zLAO0iYK9pv+OOO0zhgbYEePLJJwu1jp0WYLzyyivmcdeuXc2+vvjiC7n33nu9flUR4AN+JjI8tEjLrSy/zxyMxwMoitzfEb4zCATO12ls2bAC10Hx438JLkZ2drZs2rRJXnzxRccyHSxPa9C//PLLIm1Lg3ItGNBm93rTVgHly5c3TeE1yF+/fr08++yzLu+Jiopyea5BuT24V9oaQJvtP/XUU+a5zWYzN62F3759e6HXsdOuBHYlSpQw3RDsLRW8jRwS8DO1KkbKipEdJCMrxzw/lpkt1WMizPJgPxbOfzaC8XgAF/v94TuDQLtunR8fOJYpMRHng32uZd+dE/t5COb/JSg8bcKuwXF0dLTLcn3+119/XVRzf63119H2tSXA9ddfL2+//bbj9ZiYmHzfX65cuQvSp33xta98bva++4VZx05bITjTwgxfTe9HgA/4IX40ORYAeQmClfNvoD7mN9H3OA8oqri4ONPMXge9u/baax3Ld+3aJbVq1Sry9nRQO20ur7fU1FTTd/+FF14wg/RpH/4dO3aYmv3C0jRkZGSY/v7u0lOYdfwRg+wBAAAAADzqrrvuksmTJ8vx48fN8507d8q8efOkf//+RdpOSkqK6ZdvrxHXUfQ14NfB7HS0fd3ehAkTZN++fY73vP/++/luU6fhq127tgwbNswE8XZr1qwx/f0Lu44/IsAHAAAAAHiUDqanTfJ1ejrto96qVSvp1auX9OvXr0jb0Rr6VatWSbVq1czI+3Xq1JH9+/fL008/bV6fOHGiGeW+fv36Zj8JCQlmQL78aJP6zz77zAykp/3ltfb/sssuk2eeecaMwl/YdfxRiE1HCoDX6YAQesFr347cg0AAAAAACB75xQY6ivyePXtMM3HHQHHH94u82fL83PTeEhouMnSDSPkaRXqbNp/XKeMSExOlRg3X9+pAdD///LN07NjRsUznk9eaeZ0qz5k2zddtaX97fU37uTvbu3evOU7apL5q1aou0+RpyOtuXnptWXD48GFTcJA7fYVZ5+uvv5ZGjRq5BP3r1q2TSpUqmRYAnpDnNeAGAb6PEOADAAAAKCg2cBvcaZCfmeq9AxgRW+TgHuL1AJ9B9gA/oyWECxYscFmm83sWNDqoFSUnJ8vMmTNdlrVs2VKSkpJ8liYAQPE5efKkTJo0Kd91dL5prQmEd6xevdqMYM458EMabBNwIxf64AMIKCeyTsje9L2+Tgbc0HPz3cHvOEcAis3hzMPkMQDgBgE+4GeOZB65YNmc7XMkJTPFJ+nxNx//9rEkfZLEnzs/De713Dyw7AHOkR/QPGPq5qnyc+rP5p48BIEg9VTBzY2fWfMMeYyP/5cs/m0xeQrgpwjwAT+TlcdgKbN3zpaUUwT4zjLO/G+6EvjnOeEc+ZbmGdO2TJPdabvNPXkIAkHq6cL3JyaP8Y688o5Fvy0iTwH8FH3wAT9TvW51WXjZQvN4ftJ8c79w8fnnwSY+Pl66D+0uPRf3DPpjAQDBoHREafMbqL9/DWIbXPD6jtQdQfub6CvVm1SXKWlTZFy7cZJUO4lzAPg5avABAAAAALAAAnwAAAAAACyAAB8AAAAAAAsgwAcAAAAAFLtvv/1W/vvf/zqer1y5Uv7v//6PI+9BDLIHAAAAAAHm0MlDcizrmNf2FxMeI1XLVr2kbWgwv3nzZunSpYt5vnDhQjlw4IDceuutHkolAi7A7927t1SrVk0mTZpU6PfoBdOiRQt58cUXizVtAAAAAOCN4D7pv0mSfTbbawc7rGSYLO6y+JKDfGcdO3aU48ePe2x7CMAa/FOnTklWVpbH39O9e3epXbu2vPzyy5eYQgDFJu0ABxcAgiW/z8oWCS8nElvH16kB/I7W3HszuFe6P91vUQL8Y8eOyeLFi6VkyZLSpk2bC16PjY2VsLAw83j79u2yevVqefDBB13W2bRpk6xfv14GDhxonp89e1aWL18uv/32m9SoUcMUEkRERDjWX7p0qdhsNmnSpIksW7bMPL777rvNa7t27TLdAipXrix/+9vfzLp16tSRVq1aOd5/4sQJ+fLLL+Wvv/6SevXqmfWc/fvf/5bWrVtLaGiofP/991KmTBm5+eabpWzZsi7rnTx50uz/6NGjcvXVV0vjxo1dXi9oP0ET4M+bN09KlCjhFwUHQHGILBX5v8cHf5KI9bNkcN1uElcmLigPeNxZmwyOaiRxf+2TzPl9RWpUu+A4wT/kPiecI9/SPGNw08FSN7quuQ/WPAQBeN0m9pS4uT1Fzp47v/CqQSLtHjMPIzfMdqxLHuMdmoe0im0idX/9RuRchMStfjmo/5eg8Hbv3i3XXnutVK9eXRo1aiSjR4+W8uXLS9Wq/ysgcG6iX7p0aRk8eLAJhps1a+ZY54UXXjDBswb4Gix37txZcnJy5KqrrpK5c+fKww8/bALpWrVqmfU/+OADWbdunVlH95+YmGiWz58/X/r16yc33XSTREZGyqOPPmriv6FDhzoC/B9++MGkpWnTpnLZZZfJK6+8YrarhRSlSpUy6/zjH/8wBQSZmZkmrbqvp556yhRE2IP8tWvXmm4IWgCh25oyZYr06dNHnnjiiULvx3IBvn5gLU3Rko/PP//cnGS9v/fee12a6OtJGTlypLk4oqKi5JZbbpEjR46Yi+fNN990bC89Pd1cMHrytdRHm/qPHTvWFBZoKdGSJUvM45kzZ5r1N27cKJdffnmeaXvvvffMhfbrr7+akhs7TZsWFGghBHCxEqISZHHbiSLz+0rCZftF9q6TIS0HiEQE5w9pXM4ZGbLlc5E6XUVycmTx/oMiPeea4wQ/vHa7Lpbkk8kSXzaec+RjcRFxMqTZEPO4Xmw9XycHKPx1e1mSyNmJ/1v4w9sizfqahwnfTZXFfeeLVK5PHuMlmn/MavG4yNvtRUrHSVyQ/y9B4T355JOOGE7jrD/++EMaNGjgEuA709hPA+Y5c+Y4Anyt3f7iiy9k0aJF5vnw4cNNwK5BfEhIiFl23333yWOPPSYff/yxY1t79+41LQK0hbbSGG3EiBEmhtNgXP3nP/+RO++80/GeM2fOmFbd2q3b3opAY83mzZvLtGnTzL7t9PNovKjBuG47ISHBFCAMGDBAsrOzzXY0wJ8+fbpJp7Yi0IKAou7HUqPo64F67bXXpH79+ubg6QnIq6ZdT9Cnn34qCxYsME0stCRFS3JOnz7tsr233npLGjZsaJp9vPPOO/LGG2+Y9ZTu58YbbzSlQn/++ae52Ut68tKtWzdJSUkxF6tdWlqaOak9e/YshqOBYJMQWUUScnJ8nQy/pMdFjw/8N8hvW60tf7wBFO9vJIW8gF87d+6cqY0eNGiQo/W11lRrZWx+7rrrLvnwww/N+5XGeDExMSZW08BZg3it1NUacY3nXn/9dRM86+j8ztq1a+cI7pU28dcYz7n5/x133CGVKlVyPP/uu+9k3759pum8VhTr9jWG1G4EubevAbq9pl2b6Gvze+0yoNasWSPJycny7LPPOgoh9N7eRaEo+7FUDb667rrrTGmMOxrEaymH1rrrSVR6gPKaaiEpKck0v1BaaqTPtf+FNtMIDw83/UL0JOXuO5GXcuXKSY8ePUz/i9tuu80s01IkXa7bzYsWSjgXTGiLAgAAAACwGq0M1YA8Pj7eZbk219+5c6fb92llqda0a5ymsaBWyPbq1cu0mtagWeMpHZRPm//bRUdHm4IBZ5WcAnd18OBB0wWgQoUKLsu1ZbiddhXQwgjdj7OWLVvKFVdc4bJMCxmcaRypn1cdOnTIpFc/a16Ksh/LBfg6MEJ+fv/9d1Ojr4McOB9cbd6QW+6DpSUkuQ9qUdx///3Svn170x1ALyAN9rWwwF2fiXHjxjGKPwAAAADL01hLg1wN9J1p7JSfihUrmj7yGthrk32t0X711VfNa9oFWwNj7cqdO6AvSKVKlUzlsHb/dq7QdU6fpllbDmj3b3fBeWHodrT/v44XoJ8nr9c9sZ+Aa6KvtGY9P/amG1r77iz3c5XXwHzanONitW3b1hQavP/++7Jt2zbT7EP7f7ijg0poM377bf/+/Re9bwAAAADwVxrc66jw2srZeUR95y7O7vTt29d0z9YKVO02feWVV5rlOjDe9ddfb5rl22vL7bS/fX5atGhhavqdx0qzN6W30wH5tDuADnjnTAsG7M3vC0O3o4UR2o3AmX0bntpPQNbgF0RHGtQacx2xUPt02IP+n376yZz8ol6E9gKDwtJBFLR7gDb50FYEOmhEfoUVBRVYAAAAAIAV6PTjGuTr+GU6WryOV6ZBdkFuv/1203df3//000+7vKaD1mmcp83Zu3btaoLiFStWmJH1x4wZ43abul8d1O6hhx4ysaIWFmjhg9am2yuCtWZ/9uzZpiv2L7/8Itdcc41pbv/VV1+ZAd61RUFh6Ha0cEIHdd+6dasZMFBHzdfK4YkTJ3psP5YM8PXE3HPPPfLMM8+YgQ20icP48ePNCI1FpVMYaE28NqdwHhk/P/3795dRo0aZrgI6UB8AAAAAQMzUc1u2bDGBtFbKasWojorvXGuuc9hrn3pnOqe9Bro7duxwzF9vpwPnaW29zqCmr2vT+xkzZrhMq9epU6c8K251Oj2NGXVWNZ3mbvny5WbwPudCBx0EUGdK0+1rOnV6Py1kcO6rr622dbkzLWxwnh1An+tYA7odbb39wAMPuAwwWJj9BGWAr7Rpgx7kevXqmcEOtM/GDTfcIGFhYUXazrBhw0wpihYa6AWY3zR5dlrioydPB/XTwR8AAAAAoLjFhMdIWMkwyT7r2lS9OOn+dL9FoU3sn3/+ebev60j2edEafHe0AEArWt3p06eP28HtdOA+vakNGzaYimHn8dyUBtn5TVX33HPPXbBMA/i8Wps//vjjbrdT0H4sF+DrtAp59ZvXfhPOyzWo15IPndve3vdep8PT0pj8tqWDNTiX7Ggwv3nzZlODr009NNAvDJ3eQOdPzD2SInBJylURaT9KJL6lSMI1InFBPIe1/ViUr/m/ZeHlfJkiAEBxyp3HXzXo/G+B0t8D+2N4D/9L/E7VslVlcZfFcizrmNf2qcG97jdQ/fDDD2aUfo0TtVZ91qxZpsu1c+2/FYTYLmWkOT+watUq2bt3r+nboR9F+zVMmDDBNNlwnvuwuC4Snc/wxx9/NAM3FIVOk6fNQfTionAAKKTU/z/wSOyl9U0CAARAfp914nywT56PIJBfbKCVj3v27DE1wjrVGy7ejz/+aJroa9x49dVXF3ncNl8pyjXgtzX4haUDLOhI9o888ohkZGSY/gtffvnlJQX32qc+vyn6dGCGm2++2TTp0P7/RQ3uAVwk/uQBQHAgvwdQDK688krHqPxWFfA1+M60yX1ezfqLSg+JFha4o833T506ZUbFz2tKvsKgBh8AAABAQbEBNfg4HUw1+M48EdyrkJAQM31BfnRwBwAAAAAA/IVnImIAAAAAQLGxUMNrFOO5t1QNPgAAAABYiU7hrS2MU1JSJC4uzjxGcAX3KSkp5rzrtVAQAnwAAAAA8FM65lf16tXNPO46yDeCT0hIiLkGCjP+GwE+AAAAAPgxHR8sMTFRzpw54+ukwAe05r6wg7sT4AMAAACAn9MA72Jn8ELwYJA9AAAAAAAsgAAfAAAAAAALIMAHAAAAAMACCPABAAAAALAAAnwAAAAAACyAAB8AAAAAAAsgwAcAAAAAwAJCfZ0AAK527twpCxYscFk2fPhwiYmJCbpDlZycLDNnznRZ1rJlS0lKSvJZmgAAxefkyZMyadKkfNfp06ePJCYmchq8ZPXq1bJ8+XLOARAgqMEHEFBOZJ2Qvel7fZ0MuKHn5ruD33GOABSbw5mHyWMAwA0CfMDPHMk8csGyOdvnSEpmik/S428+/u1jSfokiT93fhrc67l5YNkDnCM/oHnG1M1T5efUn809eQgCQeqp1ALXeWbNM+QxPv5fsvi3xeQpgJ8iwAf8TFZO1gXLZu+cLSmnCPCdZZzJ8OJZwcWcE86Rb2meMW3LNNmdttvck4cgEKSeLjjAtyOP8Y688o5Fvy0iTwH8FH3wAT9TvW51WXjZQvN4ftJ8c79w8fnnwSY+Pl66D+0uPRf3DPpjAQDBoHREafMbqL9/DWIbXPD6jtQdQfub6CvVm1SXKWlTZFy7cZJUO4lzAPg5avABAAAAALAAAnwAAAAAACyAAB8AAAAAAAsgwAcAAAAAwAII8AEAAAAAsICAC/CffPJJmTRpUpHeM3z4cJk6dWqxpQkAAAAAAF8LuGnyfvnlF8nIKNr81zt27JCIiIh81xk5cqRUr15dRowYcYkpBFBs0g5wcAEgWPL7rGyR8HIisXV8nRoACBgBF+BPmDBBwsLCPL7dn3/+WXJycjy+XaCoIktF/u/xwZ8kYv0sGVy3m8SViQvKgxl31iaDoxpJ3F/7JHN+X5Ea1S44TvAPuc8J58i3NM8Y3HSw1I2ua+6DNQ9BAF63iT0lbm5PkbPnzi+8apBIu8fMw8gNsx3rksd4h+YhrWKbSN1fvxE5FyFxq18O6v8lgL8LsdlsNvFD2qz+8ssvN0H3559/LomJiTJlyhTTRL9SpUry2GPnM3r11ltvycKFCyUqKkpuueUW2bdvn4SHh8vo0aPN6zfccIO0aNFCEhISZNmyZXL27Fnp3bu39OnTx7w+ZswYefXVV6V06dJmHTV//nzH49xSUlKkR48eMnToULnjjjvMsiNHjphlDz/8sHTt2rXAz5eeni7R0dGSlpZm0g042/vrEpH5fSWh3SiRFWNFus0QadIjOA/Swc0ib7c/fww+Hih7Q0NFes6VhMTOvk4Z8rA3fa8kn0yW+LLxkhCVdx4KAIXK950NWnX+/u32srfvfJHK9cljfHFOOj7D/5JiQmwAy9fga7P6GTNmSK9eveSJJ56QatWq5dlEXwPzF154QV5++WWpWbOm/Otf/5Jvv/1W+vbt67K9N954Q+6++2556KGHZOfOneZxhQoVpHPnzibQ//LLL12a6FeuXNlt2uLi4kwQf//998uVV15p9nvvvfeagoPbbrut2I4JgkdCZBURWpTkfWz0uOjxgV/SoJ7AHkCx5jP6G0ABIgAEVoCvGjRoILNmzXL7+pkzZ+Sf//ynuQ0ePNgsa9++vQnUc9NAfPr06ebxjTfeKF9//bUsWrTIBPh16tSR8uXLS5UqVeTqq68uVNqGDRsmX3zxhfTr188E+1qosGXLFilZsmSe62dlZZmbcykdAAAAAABBMYp+27Zt83197969kpqaaprg25UtW1Zat26dZ4DvTAsB/vzzz4tOW0hIiLz77rum7/4jjzxiRul316RfjRs3zjTJt99q1Khx0fsGAAAAACCgAvzIyPwH0Tpx4kSe62mQn1upUqUuCNAvdfgBHexPbyVKlJC6devmu66OB6D97e23/fv3X9K+AQAAAAAImAC/ILVq1TL3u3btclme+3lhaJBe1ID/gQcekKpVq8p9991n+vznN32fDvqng+k53wAAAAAA8JSADvC137wOavfSSy/JqVOnzLJ58+bJtm3birwtHVQvOTm50Otr83wd3X/OnDny2muvSWhoqBlBHwAAAAAAX/DrQfYK4/XXX5ebb75Z4uPjpWLFilKmTBnTB18D7qLQUfD//ve/S6NGjUwT//ymyfv999/NNH46Yr9O5ac++OADM0CfTtNXmGnyAAAAAAAIigBfp7XLqy/9hAkTTL93Ow3CtcZep74rV66cmbJOA3xtOp/ftkaOHCnZ2dkuA/ppDb4G79rUPr9p8nT/Ogq/88B9zZs3l82bN19yv37AKFdFpP0okfiWIgnXiMTVC94DYz8W5Wv+b1l4OV+mCABQnHLn8VcNOv9boPT3wP4Y3sP/EiBghNgCPCLduHGjCapbtmxpni9YsEB69eol69evlxYtWoi/0mnydDR9HXCP/vhAIaX+dv4+tg6HDACsnt9nnTgf7JPnIwgQG8DyNfiFpTX1Xbp0kaNHj8rZs2fNtHlvvfXWJQX3Bw4ckDvvvNPt6wsXLjTT7AHwMv7kAUBwIL8HgOCswXcOyrVpfe3atS+YEq+osrKyZNOmTW5f1+b4Oir+paCUDgAAAACxATwp4Gvw7TxZo67Buw6YBwAAAABAoAjoafIAAAAAAMB5BPgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAUQ4AMAAAAAYAEE+AAAAAAAWAABPgAAAAAAFkCADwAAAACABRDgAwAAAABgAQT4AAAAAABYAAE+AAAAAAAWQIAPAAAAAIAFEOADAAAAAGABob5OAABXS5culbVr17osa9y4sXTr1i3oDtXKlStl1apVFyzv1KmTtGnTxidpQuHOV/v27aVDhw4cLh/Jzs6WcePGuSwbPXq0hIWFcU7gt06ePCmTJk0q1Lp9+vSRxMTEYk9TsEtOTpaZM2fm+drw4cMlJibG62kCkD9q8AEAAAAAsAACfAAAAAAALIAAHwAAAAAACwix2Ww2XyciGKWnp0t0dLSkpaVJVFSUr5MDAAAAwEeIDeAp1OADAAAAAGABBPgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAUUS4C/fPlyGTNmTJHft3//fhk7dqwMGTJEvvrqK/Fn06dPl/nz57t9DgAAAABAwAf4P/30U5GD3ezsbLn22mtl586d0rhxY4mLi5OpU6fKRx99JP5oyZIlsnr1arfPAQAAAADwplDxE7/++qvs27dPXn/9dYmNjTXLnn76aalbt650795d/M3gwYMlMjLS18kAAAAAAMC7Af7p06dl7ty5sm3bNqlUqZLcfvvt0qBBA/PaypUrTWCvnnzySQkNDZWmTZvK1q1bTdD/4IMPmtdefPFFqVy5stt9vPbaa5KQkGDmlddtnj17Vu666y6pX7++Y52XXnpJWrduLTfccINj2TvvvCPh4eFm3cJuZ8+ePRITE1MMRwoAAAAAAD8dZO+vv/6SVq1amWb78fHxkpKSIm3btpX//ve/5nVtjp+YmGgea/P8Zs2amQA7OjravKbP9aZBeH7+7//+T4YOHSqPPvqoREREyO+//y4tWrSQH374wbGONvnfvHmzy/uWLVtmAvmibIcm+QAAAACAoKvBf+aZZ6RWrVomcLbT2vuHH35YunTpIg0bNpSuXbvKhAkTZMCAAVK2bFmzjvbB1yb69hr8wrYU0L7w5cqVM8+15l1bBaxYsaJIafbUduyysrLMzS49Pf2itgMAAAAAgM8CfA3sNcDXWnGbzWZuR48eNc3vU1NTHX3uPeHvf/+7IyhXvXv3NoUIOTk5pum/t7djN27cONPFAAAAAACAgA3wNZi/7rrrpFGjRi7LdVlYWJhH95W7sKBixYomKNduAtr339vbsRs9erRp8u9cg1+jRo0ibwcAAAAAAJ8F+NWqVTNBcVGa2quQkJAi72v//v0uz7WVgPbd1778SgsUdEo+Z8ePH3d0CyjsdopK31vQGAIAAAAAAPj1IHt9+/aVGTNmyO7dux3LtDbcuU9+XnSU+mPHjhVpX59//rkZFM++D+3Hf9tttzkKC7RP/7p16xzr67pr1qwp8nYAAAAAAAjKQfZ27dplRsLX6elKlSolGzdulF69esmtt97q9n362j333GOmqdMa9oKmyVP16tWTDh06yDXXXCPbt283ffx1Gjy7kSNHSseOHc06VatWNemoWbNmkbcDAAAAAIDlA/zrr7/eNMu302bx8+bNkx07dsj69etNsD558mQzZZ5d7dq1Zdq0aVK6dGnHsu7du5tAW6e1y8jIKFQT906dOsmwYcPk+++/NyPzd+7c2cxnb6fT3f3888/yzTffSGRkpKmZ37Bhgyl0KMp2Bg8ebN7v7jl8a8/RDMnIynH7emR4qNSqGBkw6fd1egs6nsWZTn87Fijc+eI8+Qfn7w/nBIFyzR44likxEe7HaOJa9v054RwA/ivEpkPaW4S2DmjVqpWMHz/eL7aTHx1kLzo6WtLS0lwKDuCZH6KOr6wscL0VIzv4ZaDoLv2+Sm9hj2dxpNPfjgWKdr44T76V1/eHcwJ/VpTfG65l358TzoFnERsgoJroe0JKSoo8++yzbl8fM2aMV9MD/1WYmuairOdt7tLlq/QWdb+eTKe/HQvkL/d54Tz5Vl7Hn3MCf1aU65Nr2TvyO86cA8A/BUyAr83ztQ9/fq+PGDHioke5d+ap7QAAAAAA4C0BE+BrM/aCptlLSkryyL48tR0AAAAAACw1TR4AAAAAACheBPgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAUQ4AMAAAAAYAEE+LCcyPBQj67nbe7S5av0FnW/nkynvx0L5C/3eeE8+VZex59zAn9WlOuTa9k78jvOnAPAP4XYbDabrxMRjNLT0yU6OlrS0tLMFIDwrD1HMyQjKyffH6VaFSMDJv2+Tm9Bx7M40+lvxwKFO1+cJ//g/P3hnCBQrtkDxzIlJiLM7Tpcy74/J5wDzyM2gKcQ4PsIX2IAAAAAxAbwJJroAwAAAABgAQT4AAAAAABYAAE+AAAAAAAWQIAPAAAAAIAFEOADAAAAAGABTCYN+JmlS5fK2rVrXZY1btxYunXrJsFm5cqVsmrVqguWd+rUSdq0aeOTNKFw56t9+/bSoUMHDpePZGdny7hx41yWjR49WsLC3E89BvjayZMnZdKkSYVat0+fPpKYmFjsaQp2ycnJMnPmzDxfGz58uMTExHg9TQDyRw0+gIDz1+m/fJ0EFOBY1jGOEYBiczjzMEcXAPJAgA/4mbSstAuWbUvdJimZKT5Jjz+asmWK7E3f6+tkIJ+gftb2WZwjP3M086ivkwB4zDNrniGP8bHUU6m+TgKAPBDgA37mzNkzFyzbnrpdUk4R4DvLOJPhxbOCi7l2OUf+5ehpAnxYC3mMb9GaDvBP9MEH/Eyzds1kXNr5vrPzk+ab+56Le0ow0j7clRpXMp8/2I9FIGjUupGMSRnj62RAxPS17z60u/m+jGs3TkavHi3dS3Xn2MCvlS1b1nHdap7fILbBBevsSN3B74AXxcfHS8t+LU0eonlJUu0kxznoHk2eAvgjavABAAAAALAAAnwAAAAAACyAAB8AAAAAAAsgwAcAAAAAwAICPsDfunWrLFiwoMjvy87OlsWLF8vbb78tW7ZsKZa0AQAAAADgLQE/iv7XX38tM2fOlB49ehT6PefOnZP27dvLmTNnpGXLllK7dm354osvJDo6Wtq2bVus6QUAAAAAoDgEfIB/MXbv3i3r1q2T5ORkqVatmlmWlJQkdevWJcCHz0WeTHV5HFGmggyOaiRxZ20SjOLKxMngpoPNfWbaXsfyyFKRPk0XLpT7nJhrOZYj5SuaZ2jeUffMuaDOQxCY16276zUyO/N/j/kd8Iq6oeWlVVicuS/MOQLgW5YN8Dds2CDbtm2TSpUqyTXXXCNRUVFm+fbt2+X99983jz/66CMJDw+XmjVryr59++TUqVMyffp081qvXr2kfPnzGVlu69evlz179kj37q7zf2qhwcGDB6Vbt27F/vlgXQmlysni/QcdjyXnjAzZ8rlI6yclGMVFxMmQZkPOPzmefP7Y9JwrCVEJvk4actFzsrjtREn+T3+Jzzl7/vqFz8TZ8446XaVeEOchCNDr1s31mmAL5XfAy+qFlpVZv2wQaV+2UOcIgG8FfB/8vPrW33HHHdK1a1dZtmyZTJw4UerVq2eCcpWSkiK//vqro//+5s2bZe/evZKWlmZe0+d6y8rKcrsPbdrfu3dv0wLA2YMPPujYD3ApEnJyzA1ujk1kFQ6Nn9Jz0/Z0FtcvgOLLZ/gdAIDgqcEfP3687N+/X3bt2iWlS5c2y1588UUZOHCgbNq0STp06GCWf/LJJzJ58mQpW/Z8aeRnn31mmujrsoK0adPGFBq899578tRTTzlaDOhgfdoqIC9aYOBcaJCenu6hTwwAAAAAgAVr8OfOnSuXXXaZzJkzR2bMmGFGyddafa2VP3HihMf2M2DAAJk1a5bj+b///W9p166dJCYm5rn+uHHjzCB+9luNGjU8lhYAAAAAACwX4Gtf+tTUVNNUXmvVN27caJ4/8MADJtD3lP79+5uWAt98842cPn1aPvzwQ7nvvvvcrj969GjTDcB+0/cCAAAAAOAplmuiHxMTY6bAe+6554p1P7GxsdKlSxdTc6998XNyci4YdM+ZDuanNwAAAAAAioPlavBvvfVWeeedd+TkyZMuy7VPfn50lP3MzP9NvVIY2q9/4cKF8uabb0rPnj0lMpJpuwAAAAAAvmG5Gvx//vOfsnbtWmnWrJncddddUqpUKfn2228lIiLCBOP5DZw3duxYufzyy83Ae/lNk2d33XXXSeXKleW7776TV155pRg+DQAAAAAAQRLgN2nSxNSeOzed1/73GszrvQbrI0aMkE6dOjnW0aBc++Rr8G83ZMgQ07xfR9rPyMgw0+wVJCQkRDp37izLly83BQQAAAAAAPhKwAf4WouuN2dhYWHSp08fc8tLrVq1ZPr06S7LSpYsaWr89VZYZ8+elcWLF5sCBMBjwsu5Pg6LFGk/SqQcc7+bY8CxCKzrF77/vsTV43sD6+Tz/A54X+5jzjkA/FqIzWaz+ToR/kbnqP/ggw/cvq4FB4sWLZIvvvhCVq1aJb/88otpKVDUfeh0eTqivvb/B1yk/nb+PrYOBwaBh+sXAIAiITaApwR8DX5xyMrKks2bN7t9/Y477pAdO3aYuey1eX5Rg3ugQAT2CGRcvwAAAD5BDb6PUEoHAAAAgNgAnmS5afIAAAAAAAhGBPgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAUQ4AMAAAAAYAEE+AAAAAAAWAABPgAAAAAAFkCADwAAAACABRDgAwAAAABgAQT4AAAAAABYAAE+AAAAAAAWQIAPAAAAAIAFEOADAAAAAGABob5OAABXS5culbVr17osa9y4sXTr1i3oDtXKlStl1apVFyzv1KmTtGnTxidpQuHOV/v27aVDhw4cLh/Jzs6WcePGuSwbPXq0hIWFcU7gt06ePCmTJk0q1Lp9+vSRxMTEYk9TsEtOTpaZM2fm+drw4cMlJibG62kCkD9q8AEEnL9O/+XrJKAAx7KOcYwAFJvDmYc5ugCQBwJ8wM+kZaVdsGxb6jZJyUzxSXr80ZQtU2Rv+l5fJwP5BPWzts/iHPmZo5lHfZ0EwGOeWfMMeYyPpZ5K9XUSAOSBAB/wM2fOnrlg2fbU7ZJyigDfWcaZDC+eFVzMtcs58i9HTxPgw1rIY3yL1nSAf6IPPuBnmrVrJuPSzvednZ8039z3XNxTgpH24a7UuJL5/MF+LAJBo9aNZEzKGF8nAyKmr333od3N92Vcu3EyevVo6V6qO8cGfq1s2bKO61bz/AaxDS5YZ0fqDn4HvCg+Pl5a9mtp8hDNS5JqJznOQfdo8hTAH1GDDwAAAACABRDgAwAAAABgAQT4AAAAAABYAAE+AAAAAAAWEPAB/sGDB2XTpk0X9d4//vhDvvvuO0lOTvZ4ugAAAAAA8KaAD/AXLFgg/fr1K/L7+vfvL1dddZU8/vjjsnr1atmxY4fs2bOnWNIIAAAAAEBxC8pp8vbt2yezZ8+WXbt2SWJiolmWlJQkdevWlcmTJ/s6eQhykSdTXR5HlKkgg6MaSdxZmwSjuDJxMrjpYHOfmbbXsTyyVKRP04UL5T4n5lqO5Uj5iuYZmnfUPXMuqPMQBOZ16+56jczO/N9jfge8om5oeWkVFmfuC3OOAPhWwNfgu5Oeni4bNmyQ/fv3uyw/dOiQLFmyxNFEf82aNbJ9+3Y5duyYae6vz/V26tQpt9tOS0sz65w4ccJl+fr162X37t3F9IkQLBJKlZPF+w+amz6OyzkjQ7Z8bu6DUVxEnAxpNsTcJ9hCzx+bthMlISrB10lDLnpO9Ny8deiw4/qF79jzjnpZp4M6D0FgKeg3j98B76sXWlZm/bLB3Ktg/18C+DtL1uA///zzpia+Vq1acuTIEUlISJB58+aZ+7Vr18qUKVMc65UoUULat28vO3fuNAG/Bvnqgw8+kJo1a+a5/YiICBkxYoTUqVNH5s+fb5Z9/PHH0rt3b/n++++9+ElhVQk5Ob5Ogn8fm8gqvk4G3EiIrCIJp7M4PgCKDb8DABBEAf6MGTNM8/tt27ZJjRo1JCcnR+6//34ZOHCgLF26VLp16ybVqlWTNm3amOdly54vjdyyZUuhm+iXKlXKFAA0b95c3n33XbnxxhvN9l966SVp1qxZnu/JysoyN+cWBgAAAAAAeIrlAvypU6fKTTfdZEbGP3DggNhsNrnmmmvkwQcflNOnT0vp0qU9sp/LL7/cFAYMGzZMGjZsaIL9xx57zO3648aNkxdffNEj+wYAAAAAwPIBvg6cp/3nt27d6rJca+yPHz8uVap4rmmv1tpPmzZNfvzxR/n9998lJCTE7bqjR4+WRx991KUGX1sYAAAAAADgCZYL8MuUKSMDBgww098VNx2sT7sCxMbGmq4BY8eOdbtueHi4uQEAAAAAUBwsN4p+hw4dTP/4c+fOuSw/efJkgQUDZ84UfjTQlJQUuffee+XZZ5+VOXPmyMsvv2xG1gcAAAAAwBcsV4Ovfd3btm0rnTt3lkGDBpkB8b799lvZtGmTLFu2zO37GjdubAJ1rZXXgfdatmxpgn53dOC+2rVry1NPPSUlS5aUoUOHSr9+/cxgfVFRUcX06QAAAAAAsGgNfnx8vLRo0cLxPDEx0QTZrVq1krffflvee+890+9+0aJFjnU0ANeB9zQwt3vkkUfMCPv/+te/ZNSoUaaG3p1Vq1aZ/vxaIGDfxvjx481Ae7o/AAAAAAC8LcSmw8zD63SQvejoaElLS6PGH65SfxN54/8XWg3bKBIWKbJ+lkire0XKBfn87yf+5FgE2vUbW8fXKQpe9u9LvVtEfv6MPATWyOf5HfD9OeEcFAtiA3gKAX4edDq99evXuz1o2jrgUqfb40uMAoMkRXCEQMT1CwBAkRAbwFMs1wffE1JTU00zfXfmz59vugYAxYbAHoGM6xcAAMAnCPDzoME7I+IDAAAAAAJJwA+yBwAAAAAACPABAAAAALAEavABAAAAALAAAnwAAAAAACyAAB8AAAAAAAtgFH0fsdlsjjkvAQAAAAQve0xgjxGAi0WA7yOpqanmvkaNGr5KAgAAAAA/ixGio6N9nQwEMAJ8H6lQoYK537dvH1/iIC6p1QKe/fv3S1RUlK+TAx/gGgDXALgGwO8BVFpamtSsWdMRIwAXiwDfR0qUOD/8gZbQEdwFNz3/XAPBjWsAXAPgGgC/B3COEYCLxRUEAAAAAIAFEOADAAAAAGABBPg+Eh4eLs8//7y5R3DiGgDXALgGwDUA/hOAawCeFGJjLgYAAAAAAAIeNfgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAWE+joBVpWVlSU//vhjnq/Vrl1bqlWrZh5v3bpV0tLSXF6Pi4uTK664wivpRPHTYS727Nkjx44dkxo1akilSpXyXG///v1y+PBhufzyy818yLCW5ORkOXTokFSuXNlcB850+W+//XbBe6699lovphDFLSUlRfbt2yfVq1c310Fe0tPTZdeuXSafqFmzJifFYvT86ne9QoUKkpCQ4PLamTNn5Pvvv7/gPfXr15fY2FgvphKepP/xtm3bZv77Va1a1W3e8Mcff5hrwt1/hMKsA/+NCTZs2CBVqlQx10FR8369hjReyK1ly5ZSpkyZYks3ApgOsgfPO3TokO2aa65xuTVs2NCmh3zu3LmO9XR5jRo1XNZ7/vnnOSUWsWXLFlu9evVslStXtrVo0cIWERFh69atmy0zM9OxzqlTp8yyMmXKmHX1/vXXX/dpuuE5hw8ftrVr184WHR1ta9mypa18+fK2q6++2pacnOxY54033jDnPXeeAWtIT0+3de3a1Va2bFlbs2bNbFFRUbb777/flpOT47LelClTHPmA5hW33367S16BwHX27Fnb8OHDzXlt2rSpLTY21tapUydbWlqay/8G/Y/QvHlzl3xgxYoVPk07Ls6ePXtsAwcOtFWpUsVWsmRJk8/nZeTIkbbw8HBbgwYNzP2wYcNs586dK/I68D+pqam2xx9/3FatWjXz3X/ooYfyXK+gvH/ZsmUmb8j9H2Hv3r1e/DQIJAT4XvTUU0+ZP/kZGRmOZQT01ta2bVvbTTfdZDtz5ox5rpmxXgMTJ050rDNq1Chb9erVbQcPHjTPP/nkE5ORr1u3zmfphucMGjTIlpiY6Pgjr8Fe/fr1bXfffbdjHf3jd8UVV3DYLWrw4MG2WrVqmcIepff63Dkf+PHHH20hISG2jz/+2BHsab6gfw4R+PQPfGRkpG3btm3muf5579ixo0s+YA/wt27d6sOUwlOWLFlie+utt2wnTpwwv/t5Bfhz5swxgd2GDRsclQIa4L3zzjtFWgf+6aeffrKNHz/eduTIEVvr1q3zDPALk/fbA3ygsOiD7yVnz56V9957T/r27SsREREXNMvR5vwHDhzwVnLgJdqk7sorr5TQ0PO9YbTZlTbP1eV2s2bNkvvvv9/RdK9Lly7SqFEjsxyBT891w4YNHd0uypUrJ02aNHG5Bux5hDbB27lzp2mqC+v4+uuvpXv37o5mtXqvz99++23HOvp91+uka9eu5rk25dR8QZczm601roH27dubc6y0We29994rH374oZw4ccJlXe3GsXHjxgu67yGw3HTTTTJo0CApW7as23X+/e9/yy233CItWrQwz/W34bbbbjPLi7IO/FPjxo3lySefNF1v3SlK3v/LL7/ITz/9JKdOnSr2tCOwEeB7yZIlS0wf3IEDB17w2vTp081y/YI3bdpUNm/e7K1koZj94x//MJn0O++8I1999ZU88cQTkpmZKUOGDDGvHzx40PS7135Uzq666irZtGkT58cC9JyvXbtWXnnlFfMn/9VXX5WVK1fK008/7bKe9svt0aOHdO7cWSpWrGjyBViD9rfW/N+ZPv/1118dwZ1+3/PKB44ePUrhr0WuAc3vc18D2dnZpn+2s/vuu0/69+9vgoK7775bTp486eXUwlvcfe+df/8Lsw4CV1Hy/ptvvlnuvPNOiYmJkaeeeorCX7jFIHteogFeq1atpFmzZi7LH3zwQenWrZup1c/IyJB+/fqZGlz9wc+v1BeB4brrrjMZ9ejRo03N/e+//y6jRo1yDKDy119/mfvcAyjpc/trCGz6ne/Zs6eMHTtW6tSpY66B3r17u+QF2mLj559/NgMsKi0U0j/5OhhPp06dfJh6eMKIESNM663ExERp06aNKfDRQl+l33Nt1aH3eeUD9nVyD8yIwKKFuq1btzY1unfccYcp3LEX4tnz+vDwcFmwYIFp3WGvrdPfkMcee0zeeustn6Yfnqe1s8ePH8/ze68VATowW1hYWIHr6HWDwFWYvF8HZf3222+lbdu2ZrlWEmhlgLb8HDZsmE/SDf9GDb4XHDlyRBYvXpxn7f1dd93laLIfGRlpavf27t1r/gAi8H+8tbRVaSmsNrnUJtivvfaavPTSS2Z5qVKlzP3p06dd3qvNr/SHHYFPv/erVq0yox/rKLr6/V6/fr0J4O06dOjgCO6VNt3VYGD+/Pk+SjU8SQt4vvzyS9NKY8KECabptbboUPYRkDUvyCsfUOQFga958+bm+x8SEmLOvdbaacG/8zWgtXL24F7pbDqPPPII+YBF6bWg3ffcfe81TyjMOghshcn7tam/Pbi3/2fo1auXzJs3z8upRaAgwPeC999/33xJtdauIPapk3I350Tg0eaYGtRrgGfPpLUkVvvOffrpp47nJUqUyLP5LlNkWYMW7mntbfny5c1z7YuvLXXs10B+eQH5gHVcf/315rdAu+poQa5Oh6RNsO398nXqq7zyAf2DT+29NeifdK2JX7ZsmQnutXuWsvfLd5cPaIEQzfStSX/n8/rea4s//W9Q2HUQuC427+c/AvJDzuAF+kOuJW3aDDN3Cd25c+dcli1dutTRZBeB3+dSf3xz96HS+e7tA65o6w0tlXUO9rSrhgYBN954o9fTDM/Tc53fNWA/57kH3tT5sMkHrCH3gEgarM2ePdv0s7bT7/vy5ctdroVFixbJ1VdfTXctC8jJyXEZPFNbeGkTfR2IzV6wnzsfsP8n0ACALnvWpN97LQS2D6am9/p/wPn3vzDrIHAVJu/PnTdo7KBj+vAfAe7QB7+YaZ8Z7Vv77rvvXvDanj175J577pEBAwZIrVq1TPNtbbqtNf3aXx+BTZtdau29DoSiI6Rrf2r9s6Z9bz/77DPHeto3WzN47aev/XPfeOMNU6unfTUR+IYPH2760OqfeJ1RQZvm6jkeM2aMY51bb71V/va3v5nXNbjXGl5ttvfoo4/6NO3wDM3b9XxrtwzNC7SJthbwvPDCC451dNTkKVOmyO233276VGoBz8cff+wo9EVg0z/o2mXroYceMq15ZsyYYbpsrFmzxrHO66+/Ljt27DAjpmuFgP7J1ya4c+fO9WnacXG0IM8+aLJ+73X8FT3fOohqvXr1zHIdYV27YmmrLu3K85///Mf8N1y4cKFjO4VZB/5JA/HvvvvOPNYBVQ8dOmSuAe2Sq912Cpv3Dx061OQbOhOHFvBoxaGO0cHYHHAnROfKc/sqLtnkyZPNl1sHzsmL/phPmzbNNNfUwTKSkpLMCJmwBv1R12a52v82NTXV1MRoZq4ls7kLgjSD1yab2oxTB+LTqVJgDVqgo3/Q9Mddz6t+x/XH3E5/+KdOnWquA+3OodMh6Q+6fWo9BD5tlq1BnQZ6HTt2NIOu5Z4yVcdrGT9+vGzZssUU8uk67dq181ma4Vk6vZUW3mn3LS3Me/jhhx1dNOw0ePvkk0/M74UOyqgD8TZo0IBTEYA0ANMKnNz0++9cwKv//yZOnGgKfLSyZ+TIkVK/fn2X9xRmHfhn6628WlpohY/+Nyxs3q+tf3RaRP0d0Zk3tFuPVh7Yp1cGciPABwAAAADAAuiDDwAAAACABRDgAwAAAABgAQT4AAAAAABYAAE+AAAAAAAWQIAPAAAAAIAFEOADAAAAAGABBPgAAAAAAFgAAT4AAEXw0UcfycGDB312zDIyMmThwoUFrpednS0LFiyQnJwcr6QLAAD4HgE+AABuaID8559/uizr16+fbNy40WfHbOzYsbJkyZIC1wsLC5NZs2bJ1KlTvZIuAADgeyE2m83m60QAAOCPQkNDZfHixdK5c2fHsv79+8sjjzwizZs393p6Dh8+LJdddpls375dateuXeD633zzjXTr1k32798vZcqU8UoaAQCA74T6cN8AAPitRYsWiZaBr1q1So4fPy6RkZFy6623mlvlypXNOvr6/Pnz5brrrpPTp0/Ltm3bJDY2Vlq3bm1e37Vrl+zcuVPq1KkjjRo1umAfp06dkrVr15r7Jk2aSI0aNfJN04wZM+Sqq65yCe7PnTsnP/74oxw5ckQaNGhg9mXXrl07KVeunEnjPffc48GjAwAA/BEBPgAAefjiiy9MAK8B+J49eyQuLs4E99pEX/vAV6tWTc6ePSu9e/eWv/3tb3Lo0CG5/PLLZeXKlXLbbbdJlSpV5PPPPzcB94oVK+S5556TUaNGObavBQc9e/aUWrVqSYUKFeS7776ToUOHypgxY9yeD21NcPPNNzuep6WlSceOHeXEiRPSsGFD+fnnn6VDhw4yffp083pISIh5ru8jwAcAwPoI8AEAyIMGyTNnzjRBuXMT/bxoLfnXX3/taNKvBQHaNF6b0pcsWVLmzZsn9913n2naHx4ebgJzfX3KlCnSq1cvsw0tRGjatKlcf/31JijPy6ZNm+Sxxx5zPNeChvT0dNNKoFSpUmbZJ5984vKexo0by2uvvcY5BgAgCDDIHgAAl2jAgAEmuFdt2rQx9/fff78J7u3LtBm+9oW3N//Pysoy79FR+fWmzewTEhJMbX9eNJDXkfFjYmIcy7Rfvdbea+GAXdeuXV3ep+sfPXqUcwwAQBCgBh8AgEvkHHRrDb27ZdpPX/3xxx8muM893Z02s69Zs2ae+9AxALTJvU6TZ9ejRw9ZvXq1tGjRwjT1v+GGG2Tw4MGmq4Cdrq8tDAAAgPUR4AMA4GVRUVGmf/+HH35ogvbC0NYAOriec229FhJMmzZNJk+eLN9//73pUtCyZUv55ZdfzBgBSte/4ooriu2zAAAA/0ETfQAA3Chbtqyj1t2TOnXqZJrWa4DvTJvtp6amun2f9s//9ttvHc8PHjzoaCGgA/3pKPtaY6+j+dvp+lqzDwAArI8afAAA3GjVqpW8+eabcvLkSYmOjjaD53mCTmf3wgsvyL333is//PCDGQhPa9r/85//yLvvvmum2svLwIEDpX379qY/vrYC0L78un6XLl2katWq8umnn5qp9nQqPXtXgI0bN5o+/gAAwPqowQcAwI3Zs2ebAfK+/PJLWbp0qaPfe3x8/Pkf0RIlzFR3lSpVcrxHR7PXZRUrVnQZDE+XlS9f3rFMp81bvny5aXq/Zs0aiYiIMPtp3bp1vgUOGuBrTb3S/vY6Er8G/NoX/+qrr5b169c79qOFEzqtnwb9AADA+kJs2gkQAAAEhF27dsnrr79ugvf8nDlzRh544AEZP368SwEEAACwLgJ8AAAAAAAsgCb6AAAAAABYAAE+AAAAAAAWQIAPAAAAAIAFEOADAAAAAGABBPgAAAAAAFgAAT4AAAAAABZAgA8AAAAAgAUQ4AMAAAAAYAEE+AAAAAAAWAABPgAAAAAAFkCADwAAAACABL7/ByfRrLpC84qHAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def plot_layers(d, start, end, layers=(\"blink\", \"offscreen\", \"divergence\")):\n", " fig, ax = plt.subplots(figsize=(10, 3))\n", " keys = list(d.data)\n", " sel = (d.tx >= start) & (d.tx < end)\n", " for row, key in enumerate(keys):\n", " for col, name in enumerate(layers):\n", " runs = d.get_layer_intervals(name, units=\"sec\")[key]\n", " ax.broken_barh([(on, off - on) for on, off in runs if off > start / 1000 and on < end / 1000],\n", " (row - 0.4 + 0.25 * col, 0.25), color=f\"C{col}\",\n", " label=name if row == 0 else None)\n", " nan = np.isnan(d.data[key][sel])\n", " ax.scatter(d.tx[sel][nan] / 1000, np.full(nan.sum(), row + 0.35), marker=\"|\", color=\"0.5\", s=10)\n", " ax.set_yticks(range(len(keys)), keys)\n", " ax.set_xlim(start / 1000, end / 1000)\n", " ax.set_xlabel(\"time (s)\")\n", " ax.legend(loc=\"upper left\", bbox_to_anchor=(1, 1))\n", "\n", "plot_layers(d, 75_000, 105_000)" ] }, { "cell_type": "markdown", "id": "73c16379", "metadata": {}, "source": [ "Most blinks are NaN in all channels, and the gaze around them is often\n", "off the screen or divergent as the eyelid closes. Some blinks (here at\n", "78.5 s and 82 s in the left eye) were found from the velocity of the\n", "pupil size although no sample is missing: the layer records the detection,\n", "not only the gaps." ] }, { "cell_type": "markdown", "id": "b09eb98b", "metadata": {}, "source": [ "## Re-running a detector replaces its layer\n", "\n", "Each detector only writes its own layer and replaces it when it runs\n", "again. Trying other parameters therefore never adds up masks from earlier\n", "runs, and `clear_layer()` undoes a detection:" ] }, { "cell_type": "code", "execution_count": 9, "id": "b92a6b6f", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:40.636009Z", "iopub.status.busy": "2026-11-09T20:27:40.635630Z", "iopub.status.idle": "2026-11-09T20:27:40.836070Z", "shell.execute_reply": "2026-11-09T20:27:40.834766Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'left_pupil': 321, 'right_pupil': 307}\n", "{'left_pupil': 253, 'right_pupil': 231}\n", "{'left_pupil': 321, 'right_pupil': 307}\n" ] } ], "source": [ "def nblinks(d):\n", " return {ch: stats[\"n\"] for ch, stats in d.blink_stats().items()}\n", "\n", "n_default = nblinks(d)\n", "d.pupil.blinks_detect(vel_onset=-20, vel_offset=20) # less sensitive\n", "n_other = nblinks(d)\n", "d.pupil.blinks_detect() # back to the defaults\n", "print(n_default, n_other, nblinks(d), sep=\"\\n\")" ] }, { "cell_type": "code", "execution_count": 10, "id": "5265be78", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:40.839774Z", "iopub.status.busy": "2026-11-09T20:27:40.839511Z", "iopub.status.idle": "2026-11-09T20:27:40.892558Z", "shell.execute_reply": "2026-11-09T20:27:40.891528Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(0)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d2 = d.clear_layer(\"divergence\", inplace=False)\n", "d2.data.get_layer(\"left_x\", \"divergence\").sum()" ] }, { "cell_type": "markdown", "id": "2ff7b9cc", "metadata": {}, "source": [ "## Excluded periods\n", "\n", "Some parts of a recording are not part of the analysis, e.g. the time\n", "before the task starts or breaks. They are not *missing* data: blink\n", "detection and the quality report ignore them instead of counting them as\n", "gaps. `exclude_intervals()` and `keep_intervals()` store them in the\n", "recording-wide `\"excluded\"` layer, which masks all channels, including\n", "channels added later. Here we keep the time from 1 s before the first event\n", "to 3 s after the last one:" ] }, { "cell_type": "code", "execution_count": 11, "id": "03ac6768", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:40.895844Z", "iopub.status.busy": "2026-11-09T20:27:40.895471Z", "iopub.status.idle": "2026-11-09T20:27:40.919514Z", "shell.execute_reply": "2026-11-09T20:27:40.918455Z" } }, "outputs": [ { "data": { "text/plain": [ "Intervals 'excluded' | 1 intervals (total duration: 4.18 sec), 4.18 sec +/- 0.00 sec, [4.18 sec, 4.18 sec]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "first, last = d.events.onsets[0], d.events.onsets[-1]\n", "d.keep_intervals([(first - 1000, last + 3000)])\n", "d.get_excluded(\"sec\")" ] }, { "cell_type": "code", "execution_count": 12, "id": "830ff7b1", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:40.922497Z", "iopub.status.busy": "2026-11-09T20:27:40.922014Z", "iopub.status.idle": "2026-11-09T20:27:40.982752Z", "shell.execute_reply": "2026-11-09T20:27:40.981820Z" } }, "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", "
NaNblinkoffscreendivergenceexcludedmasked (default)masked (strict)
left_x3.130.000.970.970.355.135.13
left_y3.130.000.970.970.355.135.13
left_pupil3.135.840.000.000.356.196.19
right_x2.070.001.040.970.354.334.33
right_y2.070.001.040.970.354.334.33
right_pupil2.074.080.000.000.354.444.44
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence excluded masked (default) \\\n", "left_x 3.13 0.00 0.97 0.97 0.35 5.13 \n", "left_y 3.13 0.00 0.97 0.97 0.35 5.13 \n", "left_pupil 3.13 5.84 0.00 0.00 0.35 6.19 \n", "right_x 2.07 0.00 1.04 0.97 0.35 4.33 \n", "right_y 2.07 0.00 1.04 0.97 0.35 4.33 \n", "right_pupil 2.07 4.08 0.00 0.00 0.35 4.44 \n", "\n", " masked (strict) \n", "left_x 5.13 \n", "left_y 5.13 \n", "left_pupil 6.19 \n", "right_x 4.33 \n", "right_y 4.33 \n", "right_pupil 4.44 " ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "layer_table(d)" ] }, { "cell_type": "markdown", "id": "16f821b7", "metadata": {}, "source": [ "## Mask policies\n", "\n", "A *mask policy* decides which layers make a sample unusable. All analysis\n", "and processing methods take a `policy=` argument; there are three\n", "predefined policies:\n", "\n", "| policy | masked samples |\n", "|---|---|\n", "| `\"default\"` | NaN, blinks, off-screen, divergences, user masks, excluded periods; **interpolated samples are used** |\n", "| `\"strict\"` | the same, and interpolated samples |\n", "| `\"raw\"` | nothing, not even NaN |\n", "\n", "`d.mask()` gives the mask, `d.values()` the values with NaN where the\n", "mask is set, and `d[\"left_pupil\"]` a numpy masked array under the default\n", "policy:" ] }, { "cell_type": "code", "execution_count": 13, "id": "27ad12c0", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:40.986253Z", "iopub.status.busy": "2026-11-09T20:27:40.985864Z", "iopub.status.idle": "2026-11-09T20:27:40.997326Z", "shell.execute_reply": "2026-11-09T20:27:40.996623Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "default 6.19% masked\n", "strict 6.19% masked\n", "raw 0.00% masked\n" ] } ], "source": [ "for policy in [\"default\", \"strict\", \"raw\"]:\n", " print(f\"{policy:8s} {d.mask('left_pupil', policy=policy).mean():6.2%} masked\")" ] }, { "cell_type": "code", "execution_count": 14, "id": "9cc0052e", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.000510Z", "iopub.status.busy": "2026-11-09T20:27:41.000313Z", "iopub.status.idle": "2026-11-09T20:27:41.007204Z", "shell.execute_reply": "2026-11-09T20:27:41.005718Z" } }, "outputs": [ { "data": { "text/plain": [ "('MaskedArray', np.float64(0.0619))" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = d[\"left_pupil\"]\n", "type(x).__name__, x.mask.mean().round(4)" ] }, { "cell_type": "markdown", "id": "020bcd99", "metadata": {}, "source": [ "A policy can also be a list of layers (NaN is always included), e.g. to\n", "look only at the blinks:" ] }, { "cell_type": "code", "execution_count": 15, "id": "0f6935b7", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.011120Z", "iopub.status.busy": "2026-11-09T20:27:41.010688Z", "iopub.status.idle": "2026-11-09T20:27:41.018697Z", "shell.execute_reply": "2026-11-09T20:27:41.017734Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(0.0584495924165538)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.mask(\"left_pupil\", policy=[\"blink\"]).mean()" ] }, { "cell_type": "markdown", "id": "a73e995a", "metadata": {}, "source": [ "## Interpolating blinks\n", "\n", "`d.pupil.blinks_interpolate()` replaces the pupil size during the blinks\n", "by a smooth curve. The interpolated samples get the `\"interpolated\"`\n", "layer; the `\"blink\"` layer stays, so it is always known which samples\n", "were reconstructed. Under the default policy, the interpolated samples are\n", "used, under the strict policy they are not:" ] }, { "cell_type": "code", "execution_count": 16, "id": "b8e14bec", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.021440Z", "iopub.status.busy": "2026-11-09T20:27:41.021022Z", "iopub.status.idle": "2026-11-09T20:27:41.315148Z", "shell.execute_reply": "2026-11-09T20:27:41.314044Z" } }, "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", "
NaNblinkoffscreendivergenceexcludedinterpolatedmasked (default)masked (strict)
left_x3.130.000.970.970.350.005.135.13
left_y3.130.000.970.970.350.005.135.13
left_pupil0.005.840.000.000.357.490.357.84
right_x2.070.001.040.970.350.004.334.33
right_y2.070.001.040.970.350.004.334.33
right_pupil0.004.080.000.000.355.660.356.01
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence excluded interpolated \\\n", "left_x 3.13 0.00 0.97 0.97 0.35 0.00 \n", "left_y 3.13 0.00 0.97 0.97 0.35 0.00 \n", "left_pupil 0.00 5.84 0.00 0.00 0.35 7.49 \n", "right_x 2.07 0.00 1.04 0.97 0.35 0.00 \n", "right_y 2.07 0.00 1.04 0.97 0.35 0.00 \n", "right_pupil 0.00 4.08 0.00 0.00 0.35 5.66 \n", "\n", " masked (default) masked (strict) \n", "left_x 5.13 5.13 \n", "left_y 5.13 5.13 \n", "left_pupil 0.35 7.84 \n", "right_x 4.33 4.33 \n", "right_y 4.33 4.33 \n", "right_pupil 0.35 6.01 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.pupil.blinks_interpolate()\n", "layer_table(d)" ] }, { "cell_type": "markdown", "id": "ef230897", "metadata": {}, "source": [ "After interpolation, no pupil sample is NaN any more, and the pupil is\n", "only masked in the excluded periods (default policy). The strict policy\n", "still masks the interpolated samples. Note that the interpolated stretches\n", "are a bit longer than the blinks: the method extends each blink by a\n", "margin, because the samples next to a blink are distorted as well." ] }, { "cell_type": "code", "execution_count": 17, "id": "34baf351", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.317813Z", "iopub.status.busy": "2026-11-09T20:27:41.317532Z", "iopub.status.idle": "2026-11-09T20:27:41.579923Z", "shell.execute_reply": "2026-11-09T20:27:41.578999Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1wAAAFNCAYAAAAU8G7TAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAm2xJREFUeJzt3Qd4FFXXB/CTvuk9JCEBEjokgPTemwhYEUVFfbHAKzYsCIKin7xgFxF7QSlWVFR670167zWEJKSX3fTvuXd2dme2JJvN9v3/nmefzO5OJpPNZnbO3HPP8aiurq4mAAAAAAAAsDhPy28SAAAAAAAAEHABAAAAAABYEUa4AAAAAAAArAQBFwAAAAAAgJUg4AIAAAAAALASBFwAAAAAAABWgoALAAAAAADAShBwAQAAAAAAWIm3tTbsiqqqquj69esUHBxMHh4e9t4dAAAAAACwk+rqaiosLKT4+Hjy9DQ+joWAqw5YsJWYmGiJvw8AAAAAALiAq1evUkJCgtHnEXDVARvZEl/UkJCQ+v91wOFVKZVUfOAgefj48Ft9VZeX81tgx1vI09/fIvvoqq+17mtlzt8CrzcAmHsMMpWjHGfq+3nlKL+Hu7H0eUZd4G9efwUFBXwwRowRHDbgys/Pp0WLFtGpU6fo2WefpebNm+sN1a1du5Z27txJ3t7e1Lt3bxowYIDedtLT02nx4sWUkZFBqampNG7cOPLReeOask5NxDRCFmwh4HIPVT4+5BUYSJ7s5udX/+2VllJVcTEFhoTgA62W11r3tTLnb4HXGwCsfbx3lONMfT+vHOX3cDeWPs+o08/G39xiaptqZNeiGV9++SW1bt2adu3aRQsWLKC0tDS9OVO33HILffTRR+Tl5UXFxcV011130RNPPCFb78yZMzyA2rRpE48wZ8+eTUOHDqXKyso6rQMAAAAAAGBJdh3h6tGjB509e5Zyc3Np6dKlBqNFNvrFAiXR4MGDaciQIfTUU09R+/bt+WNTp06llJQUWrFiBf+eCRMmUNOmTWnJkiU0fvx4k9cBAAAAAACwJLuOcLFAKjAw0OjzLDCSBltM27Zt+dcbN27wr+Xl5bRq1SqeHigO57FJa/3796fly5ebvA4AAAAAAICl2X0OV119++23FBAQQJ07d+b3r1y5QqWlpZSUlCRbLzk5mXbs2GHyOoaw72E36cQ4AAB3xdK8y8rK7L0b4GDYXGiW9g8AAC4QcK1fv55mzZpF8+fPp8jISP6YUqnkX3Wrg7CiFiUlJSavY8icOXPojTfesPjvAQDgbFigdfHiRR50AegKCwuj2NhY9KgEAHDmgGvbtm10xx130MyZM2nixImax4OCgvjXvLw82fpsXphYSdCUdQyZNm0aTZkyRa/0IwCAO2HVYlmVVzaKwY6BNTV3BPd7b7ALl5mZmfx+XFycvXcJAMDhOEXAxdL+RowYQS+88AK99tprsucaNWrEA6qTJ0/S8OHDNY+z+23atDF5HUP8/Pz4DQDAnVVUVPCT6vj4eJ7SDSDlry4hzoKumJgYpBcCAOhw+MuUrP8WC5LYSJOh9D52pXXMmDG0cOFCTergoUOH+PeNHTvW5HUAAMAwsX2Gr68vXiIwSAzEWZEqAABwoBGuPXv28LLvrL8WM2/ePPrtt9/4aBa7FRUV0a233soP5NnZ2TR58mTN9z7wwAO8rDwzd+5cXnGwU6dO1KFDB1q9ejU98sgjNGrUKM36pqwDAADmN3YE94X3BtTFqWs3afpvB+mpgc1pULsmePHA5dk14AoNDaVWrVrxZRYIiaKiovhXb29v3qDY2PeKWArDwYMHac2aNZSRkUHPPPMMde/eXba+KesAAIDr+PPPP3mRj7vuuoucnb1/l5UrV/KLoPfee69dfj64luGf7OFfJyw9TpcQcIEbsGvAxYItMeAyRKFQyEa1asLmWo0ePbre6wAAgGtgQQqbf1aXIOX333/naeisSJMjsffvwgIu1v8SARfUx/9+30tf7s2SPZZTWEIRwZgbCq7N4edwAQAAmIMFGnUdEWJBCgtuXIEr/S7g/MrKK/SCLebRr7ZTfrGKftlxih77chOVqNDrD1yPU1QpBAAAqK9ffvmFzwlu3bo1bdy4kVQqFZ8n3KxZM/48m9t75MgRPirEej4yjz32GCUkJPCCSyx4uXz5MjVu3JhGjhwp6+0obrt58+a0atUqXrnvySef1DzOfsaWLVv4du68806+DambN2/y7WdlZfHMD5aNUVMzYdaXcvv27Xw5PDycOnfuTL169dI8X5/fhbly5QoP2FimSb9+/fDmg3rLKVIZfPxwZjkt33eeXlt9id//aOUhmn5XV7zi4FIwwgUAAC6JBRUsaBCx4IdVvGXBDGsLwoKudu3a8fm9NWENn9u2bUtff/0179/4/fff8/vscem2n3/+eb5tFsiw/lTSx1khqBMnTvBAiH3vvn37NN/LAqOWLVvS0qVLeWl11gKlb9++vNm0Kc6ePcuDOPZzamPK77J//37+GNtX9jqxfV+7dq1J+wJgzM2CEqPPZeQLFaSZq7naZQBXgREuAACoExZMiKXibY2N+tSnIh4r/HDq1ClN03vWduTTTz+lr776ii8vXryYF2wSR4WYoUOH8oq2rJKuiFW5nTZtGv3000+yUapz585RZGSk7GdevXqVTp8+rRnVeuihh3hwJI5QsSJOrIrusmXL+P1XX32VB2CfffYZPfvsswZ/j8GDB/ObiG2PfQ+b99y0adN6/S4sKL399tv594v32bZZcApgruxCwyNczIJdGZrlozeMB2YAzgoBFwAA1AkLttjohz2wQIIFEeYaOHCgJthiWJsQNqJjTF5eHk/fS0xMpLfeeosHm+yWn5/P+zlKDRo0SC/YYgYMGCBLIZwwYQJ/rLCwkP8uW7du5UUppJV67777bp6aaCzgYo4ePUo7duzgo2KsgiFL/zt27BgPuMz9XViDa7Y/bD0R23dpcAdgjtxibcC19YXe1Pd94YKDrmuFVXiBweUgpRAAANxGUFCQ7D4LeGpq1svaiLCghDV9ZlUCWbDJgpv27dvzOVpSYksTXbGxsbL7cXFx/Ov169f5jW1ffEwUHx9PaWlpRvfrzTffpD59+tDOnTv5XDSGzddiaYL1+V3S09Nr3GcAc+UUlfKvKdE+1Cha29rHWIENqWOXM6nlqyto6pId+AOAU8IIFwAA1Dmtj4002UNNhSQsQTddMTo6WjN6dc8995i1TVZOXUoMalhQxQI+9jPZYyzwka7TsGFDg9tjgdKcOXP4/DCWHsiwAOp///tfvX8XMbBi+8zmcUn3h42gAZgrt0SYkxjsJ/wPB3oTFcvjKo1vNx2niUO1/w8TF++n0kqin4/m0bRiFYUF4r0IzgUjXAAAUCfsRJ4FCva41Wf+lilCQ0N5qp8oIiKCp/+99957mpEkMcD5999/Tdrmpk2beCEN0TfffEM9e/bklQFZNUNWIIM9JsrJyeHzuVgFRUPYiBwrqCENgNj3647UmfO7sIqKbOTsu+++0zzP9l2aYghgjnx1wBWiEAKuvyb3pP8b3oQ6xvrqrfvOxmt05FImZReU0Myfd8nSDLMLUVQDnA9GuAAA7GzR5uPk5elB4/q2sfeuuD02x0ssasECFlZKnQUzw4YN4yNQLAhihTdYwYv//ve/vBx7bVgpdlYUg1UwPHPmDJ8jtWHDBs3zH3/8MQ+E2MgTK0zx119/8SIVkyZNMrg9FmiNHz+e38aNG8fTEjdv3qyXLmnu7/L+++/z/WGjmKxE/fLly/m8L4D6yFcKFwRCFT78a9PYcH77ZKv2YoSIhVcv/XKQYoK8adsVebGNnEIl/z4AZ4KACwDAjvKKVTRT3X/mtk7JFIpUGYs2PmZzlET33nsvH8HRDUrEPlwMa5TM5mLt3buXByNMUlISL0axYsUKXuGQBUUzZsygJk2a1LhtESs48dxzz/GgiH0Pqz4o/V62PVbF8I8//uB9uNgIlG4fLt3f5dtvv+WBGSs136ZNG1qwYAEtXLiQFwGp7+/SpUsXOn78uKYPFyvowcrGi9sAMEe+UsgfDAsQAi6RqkJooaDrdE4FXcnXzznECBc4IwRcAAA2duJqFv2x9wI9PiiFqtT9mpibhUpZwLXj5FU6cOkmPTWsPS+IAHXDghQpFhTpYgGXLpbix25SrNAE63VljKFtS7GRInYzhgVGjz/+uMm/C0utZKXb2U3Eyrdb4ncRKxNK+3qxZtEA9VHAJmHxgEueQlguKUrYItyLzuRqW06UGShYuO54Og3vaLgSJ4CjQsAFAGBjzy7dT2dzK+lGwX56aaR2YnhekTx15oHvj/CvDcMD6K7uLfF3AgCndTRTmMMVHmg84EoI86Mzudo+XJUGBr+WHc+n9624nwDWgEumAAA2cvDCDfp6/REebDHbLhVSabn2am55ZZXBssjT/j7Hv+YXq2j+yoOUll2Av5mTYCNfI0eOtPduANhVRWWVJrCKCvaXPdcxzo9/jQ30oOxi4y0apEpKtev9suMUbTqqPw8MwJFghAsAwEbu/FLeYJdlE5ZXStJnKrTLBy5kaJZZJs7u02n0/fZztOpsES3Zn06PdGtID/VtTYEK/Qpf4DhqSzUEcAdXsvI1y/3ayguwfPxQd/pi/XF6qE8L+mLDCTqcmWdwG36eRKXqoO3nHafJ29ODAvy86eW/z/PHLs3VNhcHcDQIuAAA7KSqmo1kaYOsH3acp8SoYAoNUNB93x2Srfvbvos82GJuFFfT3I3X6HhaPs1/VD4/BwDA0WQVCGmCob5Evj7yU88GYUH02j3d+PK02zvRyRtb6XCm/kjX/R0iafHBbGI1Nt5Yqz+iVaQsoyB/XIACx4SUQgAAOwZcKknAtf5CCQ34cCcNeHej3rq/HdNeIRatOqPtsQQA4Khy1PNTA31r7qPHiga9OMxwcZmXR3fkwZYx6bmoogmOCwEXAIAdA65iyVwEUb4wt7xW7OTj992nLb9jAAAWdDVbCIaCfbWtDozx8Zavw/okvzOqKQUofCmkhgGszPzi+u8ogJUgpRAAwAakPZQ0j/E0GOOTxHVLJBsy5U+hoAaqGAKAo9p7MYd/9fOueYSLkRYSOvHGUArw0/bt+ntybxr9yXaDF6Wu3CykXpbaYQALwwgXAIANKMv0G3iyohkXs4ynBS6Z2IcitW25jFpzLL2+uwcAYBUnr96kw+lKvhxkwghXiKRPlzTYYhrHhNLO6UMMft9rqy5STqHwczYfu0Lnb+TWc88BLAcBFwCADfy664z+AdiDaPlRbTVCqQBvoujQQPrukS61bvvoDW3fGnA+paWltHv3biorK3OI7QBY0riv91CWUph81bVJeK3r35IcS490jKRpAxMMPm+sMisrOz/k/Y28/cYji4/SoI921nPPASwHARcAgJWxE4BZa/SraqlqyBZsFSWcVLRrEkMjmgfVuP3rRdV0AxPGndbVq1epR48elJmZ6RDbAbBkKnVuqfa+btNjY2bd252eHKptCq/r5f4NDT6erSLafgoj/uB4EHABAFjZqTR5aktqtDZN5mK+MLfr83vllbleG52qWX55ZHvqnaig929vZvRnIOACAEdTUipPpY4MNiFH2gT/Hd6BHuxgeLSsRJK+rTKQyg1gDwi4AACszIvlDkr8+vRAvXXYXIXByQGa+x2SYzXLTRqE0eKnBtHdPVrS7a1CDP6MGqolu62zZ8/SxYsX+VX2U6dO0dGjR2XPs+eOHz/OU/GMuXLlCl+nUtKgWsS2e+LECTp27JjB5/fs2UP5+flUUlJChw4domvXrmmey8nJ4fujVApzTgxh22Q/+9y5cwaLrpi6HQB7OXJZPto6MNVyzYnfvLc73d02lP7bo4Hs8X8vaxsnl5QivRYcAwIuAAArK1TJKxEqfL0NBlxvjelCPRMUNP/ulka3FeinnXTer7G/wcpeIJg2bRo99thj1L59e7rzzjvp3Xff5Y+fOXOGOnbsSD179qSxY8dSTEwMLViwQPaysQCtU6dOlJKSQmPGjKHk5GTatm2b5vmdO3dSUlISDR48mG699VZKTEyk9evXy7bRq1cvevbZZ/l6EyZMoNWrV/PHP/jgA4qPj6d77rmHGjZsSO+8847en+yvv/6ihIQEuuOOO2j48OHUpEkT2rJli2wdU7YDYE+nr2uDnxmDE/WKYNSHp6cnvf9Qb3r59s6yx/ela4OsYp1jL4C9oCw8AICVLTt4Xe+xoU0Dae15bd+YQD8fig0PoqWTB9W4LYWPNuD6ftJA6vLGSj4hXdpA2dqqq6tJaacAz9/Hizw8ai8tLWJBysaNG6lv3778PisoMWLECBo/fjzNmDGDn7QdPHiQ+vTpw4MwNgeKjXixIKpr1660Y8cOUigUlJaWxpcZlUpF999/P40cOZI++eQTvj9Tp06lcePG8dGokBDtKOS+ffv4CFh0dDS/f+TIEXrppZd4QHXbbbfxEbABAwbI9vnkyZP0wAMP0LJly2jo0KH8sS+++ILuvfdePmrHtm/KdgDsTazCykbvHxvczmo/h83pemdzmknVYQHsAQEXAEA9VFRW0fErWZTaOJqfvBui8NY+/sEdwjysAJ3yyIEK0678NgiRz4EQ4y9lme2u5LJgq81ra8geTrw5jAIMjBAaM2jQIE2wxbBRJpbaxwKZ/fv3awLIdu3a0YoVK3jAtXLlSp5KyIIlFmwxbASJBTzMunXr6Pr16/S///1PE/zNmjWLPv30Ux4APfjgg5qfN3nyZE2wxSxcuJD/DBYkMaGhofTKK6/wkTbRV199RW3atKGIiAi+D2z/OnToQAUFBbR3714+qmbKdgDs7WqukOqaGK4djbcGNqfrpwPpdKVAnnprqLE8gD0g4AIAqIcXF++gP08W0NO9YumFUZ0MrpOvLkf4xrDGmgbFl3NUsnWCTAy4Hu7fljadzqLOjcP4/TCFF10vqqAfdl7k8yP8fITDOpvzcykzn5Jjay/D7MoaN5bPGTl9+jQPkp577jm9dcXgio0isQArKirK4DbPnz/PU/lYkCPy9/fnqYPsuZp+Pnu+ZUt5ymirVq309vHSpUs8WJNiqZFiyXdTtgNgb2nqDsXJ0cFW/1kVlfozWeeuOEY/P62dDwtgLwi4AADqgQVbzFe7b9D4vsU8NTBAp0/M+Twh4GLPidIK5FdeI4K1BTNqwuZ//fT0YL3Hd10rpVEfrKe1U4fz+1MWCYHg1AENadKwDmTptD420mQP7GfXhe6oIwuqvL29adeuXUZTEwMCAvhokjHBwcFUVFSk9zh7jD1X089nzxcXa1NJGd37bB87d+7MR9xq2ofatgNg79H/s7kspc+DmsZqL05YS5mBujJ70owXxAGwJRTNAACwADaI1WXOZuo1Z53scTbSJB5om8cJo1LM7Ntba5Yf6RRp9s/tnawdwTqTW6kXCM7fpj+vob5YoMLS+uxxq8v8LUP69evHAyOW+qeLVRNkWAoimxO1efNm2fNs7hbDgiFWHVBMSRSLbLBRKfZcTdjzW7dulTUnXrt2rWyd/v3787lnLG1Ris0tE6shmrIdAHsa+cE6qibh/7VlfITVf16YH05pwXHh3QkAYEHSJp9MobKMxAuvSQ20V3m7tYjTLPvolI2vC0OpiOzKsqgK9eJl2Fwtlqr38MMP0/vvv88Lanz77bc8yBKrALJ1Hn/8cT5n67PPPuNztqZPn05vvfWWJrWPFchgz//yyy/0+++/01133cULbUjnixnCtuvl5UV33303n0/24Ycf0kcffSRbh1VWZHO4WBGMRYsW8eqHbD02j6uwsNDk7QDYC2/FkK0tWBEZYtoIfn1MGdKc4gM96IV+8Vb/WQB1hZRCAAALKyuvIF/1XKrcIpXm6laQJNUw0E+77O1l/rWvID/9w/g9H2vLk6unj7mlFi1aUFiYdlRRNH/+fOrduzf99ttv9Pfff1Pz5s1pzpw5vIy7iFUFZEUp/vjjDz6yNHDgQJoyZYrm+e+++44+/vhjXuCCFbVgVQWlzzPdu3fX+/mBgYG8vPybb77Jy7izeVf//PMPvfDCC+Tn56eZD8ZG11iw9+OPP1J5eTkPtlhgJW7PlO0A2ENBiYqeX7TL5j93RKem/MbsvZhD266oKABnueAg8FYEALCwIlUZRagDrrxiIeDy95bP5/GSBFn1CbhC/OXzxZhDGfL5YR/9c4CeG9mR3A2rImgMq+ZXU0U/lrr46KOP8pshvr6+9OKLL/KbMdu3bzf4OOupxUbVpHbv3q03j4wFT+xmjCnbAbClrPxi6vfuZiqxczX26aPa060L9lB5lTqt20gFWQBbwTsQAMDCrmQV8A95JrdYyDEM8PGwWCGIulY3/Gh7utnbBwAw1aGLmXYPtphG0UIvPBZw5amPwQD2hIALAMDC7vhyPzWdvpIvZxeqjAZcd7QOoUbBnvRgX/PLeYcE6I9wAQDYw0318c7e/CW9+sYsEOZmAtgTUgoBAKyAVec6dOEGvfjXOeFga6AwxkcP96n3z/HUqdpXWq5/eTnQtBZfAAAuEXBJUwhZW46z13OouQ0qJQIYgxEuAAArmfLzQb1eXJZ2S3IDCpfUSHjws02a5dRoIdIK8K5fKXUAAFNkF8nT94J9iN66tYndX7yPVh+19y6Am0PABQBgJRfyteXZ+zX2t8rP8PPxpvUvDNDc33dd25fp9dtT+VdVBWrDA4D1rT2TI7t/+I1b6cF+be3+0jePkTckB7A1BFwAADYw4/YOVtu2QjJfQSo0QBj6UlpoEjsrfw6A9wYYUqwqo+tF8mOEo1QHLJP0JgSwB8f4TwAAcHHWnD/ARrkMCQsUAi42wGVobpepWINdpqxMO3oGIFVSUsK/+vhgwqC72nL8KjkSaSZ1oUreKgPA7Ypm5Ofn06JFi+jUqVP07LPP8gaUdV3n559/5g0gpeLj42n69Omyx9LT02nx4sWUkZFBqampNG7cOHw4AIDVDUoKsOr2jfXxEke4mPziUooJM++Q7+3tzftCZWVl8WOmo1y1Bvtjo54s2MrMzORNmcXgHNyPbsXUF/rGkz1N6hlL83fc4MuFKgeoVQ9uza4B15dffkmzZs2iAQMG0NKlS+mee+7RC6ZMWWfTpk20Z88emjBhguaxyMhI2Tpnzpyhnj17UteuXalbt240e/ZsWrhwIa1fvx4fEABgVc8Nt/4chie7RdMXe7Jkj/lKRr6WbD9Fz4/sZNa2WRPguLg4unjxIl2+fLne+wquhwVbsbGx9t4NsKOKSiGdkBXx2fXqMKOpzrbywqhOlFW4g346koeAC+zOrv8NPXr0oLNnz1Jubi4Ppsxdh2natClNnjzZ6PNTp06llJQUWrFiBT95YMEZ+54lS5bQ+PHjLfL7AIB7G9DEnzZdUtKbw5tQcoMQeviHIzSkaSClNo6x+s+OCam5KMe87TfIx+sgTb71FrO27+vryy92Ia0QdLFRT4xswZRfhUqALJPP3sGWKJYfF/Now8USyswrppiwQHvvErgpu/5HsLQ+hgVT9VmHOXfuHL388ssUGhpKffr0ob59+2qeKy8vp1WrVtHHH3/Mgy0mISGB+vfvT8uXL0fABQBmY8l14nTsh3om00cPxVJooILfPzgzmoIUtmlMHBUk/EzRrc2D+Nf4IA/NRPb3tlw3O+BiWCqhQiH/OQAAJaXllKOuCC9+dQTB/to5hV9vPE7T7+pq1/0B9+USifgsiGrQoAFPabh+/ToNHz6cJk2apHn+ypUrVFpaSklJSbLvS05O5qNnxrDvKSgokN0AAOTHH+2yl6enJthiQgIUNpvvFB2qHeGaf3dLWvCo0FR527ThNvn5AOC+VGUVDt8X7Mu9WXTy6k277g+4L5cIuFhxjDVr1vCvCxYs4KNZn3/+Oa1bt44/r1Qq+dfgYHkfhpCQEE1lJUPmzJnDR8zEW2JiopV/EwBwNtJK6R2SrJ86aEyrhhHk7yU0Gh2U2kgT6Hl5edKvE7SjWpUojwwAFlafKqjWdE+3prL7b/x5yG77Au7NJQIu3UCoX79+/LHt27fz+0FBQmpNXl6ebD2WpsiCLmOmTZvGKySKt6tXHavkKQDYn5hO+M+kLrLRLVuLCA6grS/3py0vD6QAnTTGFpKS9EoHvRINAM5LJQm47msXRo4iOTacEoK1p7rnctDaAuzDJQIuQ9jE7ooK4QDQqFEjHnSdPHlStg6736ZNG6Pb8PPz4wGZ9AYAYGi0KDK45qIVthAdGkgRBvYjUBKAFalwwgEAlqUqq9Qszx7bw6Fe3iBf7anuTWU1NXllBf2287Rd9wncj9MHXCyo2rBhg+yxH374gffaGjZsGL/PUmvGjBnDy8CL6YWHDh2inTt30tixY+2y3wDg/CqqqmTztxwV69MlNgEtLkUDUACwLPG4wlKaWRqzI6mqkuR9q7341zm77Au4L7tWKWS9s1hD4+LiYn5/3rx59Ntvv9GIESP4zZR1WMGM999/n6f/tW3blhfI2LVrF59/Ja1UOHfuXF6VsFOnTtShQwdavXo1PfLIIzRq1Cg7/fYA4EojXN5ekuoZDqhCfc6x/WQaNY0Nt/fuAIALWXdUmHIR6Ot4x0H9cEuQnlNIcRHyuf0ALhlwsUIUrVq14sssEBJFRUWZvA7r/bFy5Uo6cuQIHTx4kMLDw2nx4sW8SadUTEwMf54V12CjX8888wx1797d6r8jALguZxnhknp9zWV6eECKvXcDAFzIqRtFekWEHMUDXRNo1hr9hu2PfrOTVr8kZEIBuHTAxQIpMZiqzzpMu3bt+K0mbE7W6NGj67yfAACGVFRWy9L2AADcTVVVFW25LEzXeLhrPDma8f3aUNuECJr++xE6m6uda+bn4FkJ4FpwhgAAYKZKJxzhYqXjAQAs5cCFDM1y39aOF3CxefxdmsfT0ol9qEW49gAY6IuDIdiOc5whAAA4oMoq5xnheqGvcCLUOMyuiQ0A4GIuZRbwrz6eRCmN7deL0JQqrp0baatNF5ZqR7sArM2xzxAAAJxkhMvTwbNTkqKFyeHFZdp9BgCor8wCIZ0wJVre/88Rqcq1QVZRWSXlF6to0jdbaPWBC3bdL3B9CLgAAOo5wuVB1TxtxZGFBfnxr8XlCLgAwHKyi0uFY4y/44+e923ZQLN8Mb+K2v/fBlp1togm/iLv0wpgaY59hgAA4AQjXJ4eDj68RUThgULAVYI2XABgQd/su8m/hgU4fsA1ukszmjE40d67AW4IARcAQD2rFDp6OiETEeTPv6oq2X5jlAsA6u/Y5UzN8v6rQml4R8YyER4b3I4UXoarLQJYCwIuAID6jnA5wSsYrg64mIsZeXbdFwBwDdI5UT2TtAUpHF0bA/PNCpRltP7wRZr1y25clAKLc4bzBACwscrKKnpr2R76ZccpvPYmBFxOkFFICl9tus8rv+63674AgGsoLa/QLL8xphs5iwpJhVlRVn4xPfbjCVp4IJs+WXXILvsFrsvxE24BwOb+2X+evuZ5+Tfp3l61Nx53V2JmnjOkFEodulFm710AABegLBMCrgYBHuTn4zynlNfy9SezZuaXaJb3XMq18R6Bq8MIFwDo2XZG28iyTHIFE4wVzXCuVyYuCId+AKg/ZZmQUujr5VwHwb7J+umPNwtUmuUD6dplAEvApy4A6PntWL5muUiFsna1lYV3loCrY6wwb+FaYRWVqDDKBQCWGeFytoDrjXu60PN94mSPZRdpgyxzeyKzC5Tv/fUvbTxyqb67CC4GARcA1Gju8gO0/1w6XiUDxGp/zhJw3ZKovar75u/77LovAOA6RTP8vJ3kIKgWGqigZ2/rSAsfTNU8tuVMVr23u2z3WfpkZwZNWHrMpPW3Hr9C89YeowpUSHR5CLgAQE+Y0LKJ++VYHj3900G8Sk5epZAJ9ddW5tp9qcCu+wIAzo0VmVi056pTjnCJ+qc0Ij/1AXzLZWW9y8Sn5wnzwKrJg4pNyCJ46qejNG9XJv1yTWgeDa7LWc4TAMCGFDpXK68XVdPiLcfxNzCaUugcJxsNwwM0y5cK0HMGAMz31Pe76GyuOMLlvKeTSWGGi32UVdT9GJmn1AZZby/fTzdya+5NVqjO2D+QY2YOo4UrTn61/ghdzEDBEGtw3v8QALCaCgPH/hmrLtHaQxfxqjvxHK7kBqGy++k5hXbbFwBwbnuva0dlgv0MdBJ2Ep893LXWkvem+uFAjnb5YA5N+3lfje1XHAkLEGevv0p3fbrT3rvikhBwAYCeimr9HiXMTkn1QnC+KoW3JMfK7mcVyFNoAADMEarwcdoXrnG0/EKU6Kcdp+u97W2XtaXmdRWXakfDHGGAcMMZIVjMRXajVTjAnxgAHAXLOf9h83HKM3LAdZLMOZupqHSuES4mKVR72C8pRQVKAKi/2FCF076Mnp6GT4XnbLxGu0+n1WlbzcPlI301jWEVStIPy6ur6YayhDKLi+lmSQnlqpRUUFpKxeVlpKqo4EU1qoxcCLWUMvXnGViH83SpAwCr+2jFQfqKNzw2rNLKB3xnI34AejpRxHVn+1j6YOt1vpxXjEuZAGCeAG+iEnXW3YSBbZz6ZRybGkY/H83Te/xads1zsHT56BQPUWedcxl5RbTx6BUa3aUZBSp8qaBEG3AVV1bRyTz9n6+LzRf29fIiPy8v9VdvCvBhN18K9PHhj3uYeWUUAZd1IeACAL2UAlHvRAUdTFdRcYX5k4hdWbk6B9+ZCnQ9OSRVE3AVSK6wAgDUhbrnMW14ridFBGsL8jijJwe1oZ+P6s9dOpaWS/fUYTulFfoXJVnBqSYxIbRgw2nada2UTl7PpzfH9qAiSRXDViGeFO7rQ9UewoVNNj+4srpKvVxF1ZKLfGzEi90M8fb0pFA/PwrzU1Cowo9C/RQmF3US/55gHQi4AEAzQfhCvjygevu+LvTNppP07b/CqBcCLsN9uJwp4PLz8abUaB86mlVOeSUY4QKAumPpyGJsERYg6SPipJJjw0nhRaTSCToW7s+mWWNM347SwEVJVnBKihXTCA/cT2czhNGzuEAPeq5ZMHkGBpKnn/5rWV1dzQMuFnix1MLSykoqq6zkX0srKqikopyKy8pJWVHOn89WKvlNDMAi/f0pOiCQovz9yctI+iRTViXvD9a3bSPTf3GoFeZwAQD37l/79V4JX28venHkLZr7WUUYETFUpdDLiVIKmQbBwgT3Taez6HJmvr13BwCcTIHkYk2wpL+fM9v2cn8amKQ/UldQojLp+3MKlbyFCrP44XY1rjtv+w1aeVYIuIJ8az4VZymCbJTKx8uL/H18KEyhoJjAQEoMCaFmERHULqYB9UhIoP6Nm1CXuHhqHhFBMQGB5OPpyQOwjOJiOpaVSduvXaXT2dlUVGb4c1waKx64aHxqAZgHARcAcN8amLvl7+tDAQpfah0pDIZvu2LaB4/7jXA5V8CVGO7Pv7L0ln4fbLf37gCAkwZcPp5Evj6ukSwVHRpI3z45QO/x+z7dYtL3j563WbPMRv2kBYpq/LlsMpwFsKAsxM+PGoWEUmpMDPVJbESdYuP4fTa3iwVf1woLaM/1NDqccYPyVfLPc+lulJTVvSQ+1Mw1/ksAoF6KlGUGqyn5qmvV3ijCwdeQCrEPl5NduooIcv4UIACwn5wi4WTd3w3OIk/crODVCmNCA3jqoTHXirSfohFB/pQQ6ksX82u/SNmhYTBL6CNLYyNjbDSM3ZqFh1OOUklpRYWUVVJCN5VKfovw86Nkbx8KJHk65ckb6NFoaU52mgAA1pCea/jgKl65/N/trTSP/baz/r1JXIWzjnDlS6pjMSpczQSAOsz3nfLLYb6s8HauY5+57vvuEA38aCflFBrvqyUVFqSgqbelUqe4mtMtBzTxp2eGWL/CIwu+IgMChPTDhgkUFxRE7C+XU1pKW/KK6PtNh6hcctWVZbOk5yDosiQEXAAgK0/LDE4OoEc6RWrut26oXf5up3wCsDtjKRqMl5MdSW/qzMXLN3GOAgDA9B930/ViYXQ/zM/JDn4mmH93S6PPjZSkDdYkwM+HUhrH0LJnhxhdZ3zHCPpu4kDy9bbtMGGAjw+1iYqm7g0TKMpPQQduEr2xVV6hmFlx4CJdzMiltOwCm+6fq3K9/xQAqLN8yQRoFmx9/cQAmjWmu+Yxfz/tBwI6cekXzTC17K6jeHlkB9n9P/eet9u+AIDzqKysoj9OagvtPNk3iVzNqC7NeKl7Q8SiGMayHQx5slu0ZrlRiHDazQYGJw5KIXtigVdqRASlqoso6Xpr/VUa8tFOGvzBNiqRlLAH87hB9i0A1CZX0gA3LED/4Buk0KZFBNZSUckd53A52whXfCSbM6C17GA6PTm0vd32BwCcw+SF2zTNfFmRhTu6NidX1DQ2nObclkTTVlw0qUT+3/uMX7SaentnOp+1hWKC/eh/9xsO5OypEZ+IZ3ieNiv9X1FJlJZTSM3jtZkuUHcIuACAruUUa16FqaO0ZeBFgZKAKzYYBRdElZXCmYe3k5WF19U4QmHvXQAAJ7BKXcqc+eL+VPJytqtNdXB/nzYUFx5Ejyw+WuN6bV5fW+Pznp6ePGvEUQX51P75lV2oItcMrW3Hdf9TAMDkCdAfbkvny22ivHlpXGP55kyhChUL9eZwOVlKoa5154vpszWH7L0bAODAdIvrNAgz/FnhSnq0jKeOsb78s1H047aT5EoUJlwwfPWPmoNOqB0CLgA3dy5dO1k2u0RSF1aHj/pK5ubLSiorR9DlzI2Pmc/u1VaeZN7elGa3fQEAx3cpM8/tAi4/H2/6/bkh9NWjPTSPfbip9jRDZ9JQIQ8FLs29TW+d83nGzw3ANAi4ANzczQJthbpydYqcIapy7QH3WjbKxcrKwjthwHVrx6bUNMzL3rsBAE6iSFUuux8a6D6pyP6+2hGuUEllxip1loPUO6OakjNhKY+6dTP+ntTFXrvjshBwAbi54lLth6i/j/FDQrUkFiuQVDV0Z9qiGc4XcDGz70y19y4AgJO4WaAkdxUmCS6bRmqXlTpplqym1L295NkDzuDDO1pQfKAHzbuzBb+f2jiGL4uxpQdVGwwuwXQIuADcnPSq5ftj2hldr0Iy+pWn07fL3VMKnbVoRveWDUk8dwituT8nALi5jzee1SxP6RtP7oSNAt3XLowvrz5XrAk+8iUVfpkZQ5uQMxrYJoF2zhxBt3fTlsZgy3tfHcSXq8mDSkoxlaA+EHABuDkx4EqN9uEn4MZUSoa4pH273Jmzj3Axs0e35l/zy4hW7kc/LgB7c9Q5siduavfrmRH61WxdXU6x9uJkdqEw2pdVUMK/spEgNvdpfP+25EqkLWGWbnetYiG2hoALwM0Vq69aKWopDfvUEO0HSaESI1yyohlOXKUwOsRfs7xwx8UaG3gCgHXtPp1Gqa+vobeW7XGol1p6MaZ9jOFGua6uSnLRcd6qw/T2n/vo9i/+5fcVLtpkSVr2f/7Wa/Vumu3OEHABuLki9RyuAJ+aCygkx4bz8rhMgVI+edpdOWvjY6mokADNcp6ygjq+sYoe+3KTXfcJwF3N/PMolVYRfb3vJjmS73Zc0Czf3dF4JoQrK63QBgyLD+XSZ7sztc+5QRG/wnLDRUJM8faf+6jta6v4BQV35cSnCQBgCZmFQnpgeEDtVy0D2Yxg3osLAZezl4UXJUQGa5bP5FZSQRnR+gsl/IMVk6QBbEt9SHE4+65rsxqigt2nOqGpfxuVCwdc0o+3W2atMutz4bPdmfw1YhcU3FW9Aq7Kykq6evWq5fYGAGyqSFlGf5wo4MtxobV/iEYECCNcGZJS8u7MFQIuljLycn/9K9Z3zVtPA+eu0Wt2CgDurXsL9yqYIZo20j2rui58UPt7s7m+eTqFQmrz++7TDn9BwWEDrpKSEpo4cSIFBgZSo0aNNI8/9NBDdPSo+0avAM5m/ZFLmuUQ/9pHuOLDhKAsqwhzuFyhSqGoQah2HpfoUEY5XSqooo1HL9tlnwDcfZ6Qoygo0V5gu7ttKEUE6x8v3EFK4xg6P/tWcje9WiXI7peaWNTlP19sot5vraIpf56TNVA+dOEGuSOzAq6ZM2fS4cOHae3atbLH77nnHpo1a1adtpWfn0+ffPIJTZ48mc6ePWv2Ounp6fTuu+/Siy++SN9//z2Vl5ebtQ6AOwlSaIMsUz7nQ9UjXNuuqKhYhaBLrNzo7enc2dkt4sONPleI+XoANuOA8Rb9uVdbMGNkh0RyZ9IiElJPdI0mV/6dpYFmj3e20uXM/Fq/b+PFErpWpJ9++MUm7YiXOzHrLOHXX3+lH374gfr27St7vFevXrRu3TqTt/Pll19S69atadeuXbRgwQJKS0sza50zZ85Qamoqbdq0iYKDg2n27Nk0dOhQnvJYl3UA3E2ZZBKwKUIkAdofewxf/HDLlEIv5x7hahprPOD654j7TnIGsLUqBz7OMe0au25gYa7FD7ej6Xd1JXcKNN/844DRdd9Z/i+1enWF0efzlO6Zpm5WwJWZmUnx8UIOr4ekHDIbMapLANOjRw8+YvX222/Xa52pU6dSSkoKrVixgl5//XXauHEj7dy5k5YsWVKndQDcjVIyP+fu7s1qXd9bctBVluNihaukFAb4GU8nZaOZGySppwDgXiNcxaXaXo2RkqqmIAgP9HO7l2LDxRL6Z582VVDqn2NZNRYR8XbuhBCzmfVrs8Bl/fr1egHX559/Tp06dTJ5O2zEic0Dq886LMhbtWoVjRs3TrMvCQkJ1L9/f1q+fLnJ6wC4c8CVFOpJDcKCal2/rEJ7FHX3fk2sp8imS0LzSy8nTymszYSlx6nv7FX0y45T9t4VALcJuByhb1F6TiG9t+U6Xw7zr7l1iLt4pFOk7H5ksHsEoQ92kGdCTF522uB8LmWF/lWDsalh9GQ3YXQ0T+meF2vNOkt444036OGHH6YZM2bw+yzVb9SoUTxNj40e2dKVK1eotLSUkpKSZI8nJydr5nuZso4h7HsKCgpkNwBX8uM+oZFhuYmf69JRrbc3pdFX64+Quzp4McPgSKGzaxZm+KTqSmEVvfy3di4HAFheFVXXuTiBNS1Yqy2EFuTnot1962jWmO609YXepPAiahjkSTGh7hFwPTmord5jBy9oPweZ5XvPUpZSP+CKCPSllIQIvpzryjX0LR1w3XbbbbRs2TI+ryooKIheeeUVHoysWbOGBg0aRLakVApXmNm8LKmQkBBeTdHUdQyZM2cOhYaGam6Jie49WRRcy5FLmXT8pvCBfq3QtIhrePvGsvuz17tvW4ibBdpjR54LFJYY0ESoPDaxbxIdmjmI3hgm/1uLUCYewDYjXPa8kFNWXkFjPl7HG/yK/H1ceyS/LhpFh9Lu6YNo3YuDyNPFMxxEidGhtPeV/tQjQdtC5tyNPNk6z/5+RrPcoYE2Vb1Nw3Dy8/HSmxPoTsx+l7DAasOGDbyCYGFhIW3ZssXmwRbDAj4mL0/+R8/NzeUBlanrGDJt2jT++4k39BwDV/Lp+hNmHXBvb2X8f8adFEiCLOeewSX46vH+tOn5nnRPz5YUFqigJtGG/86tXhN6c32y6iB1nLWCmryygj78Z3+t22fNMh/6dAON/2wjGioDGCGtY1SinjtlD5uPX5U1O2ZUpqZCuAl2nAxQCJV73UVMWCC1jNFO80nLNT5o8fPkgZrl7i3iyNdbCLjqWKvLvQMuNg+KBVc5OTkGn7Ml1geMBVQnT56UPc7ut2nTxuR1DPHz8+MBmfQG4CryVdqrp13iTf/QaBLlHukTtSlQak9GXOF6HSuIktRAm6Of2jjK6Lod3ljD53XkqNvzzNt+gwded320jvKLVXT/JxvorWV7ZN+z/shlXoBj62Ulnde5KgoAgrJK7dFEZaeUwgUbjtMTP+lfkFMi4ALWAqpbsuZ1+Gx3Jq0+cEEzKirl5+NNx14fQrun9qPo0EDyVVfLqMAIV92wFMIuXbrQsWPH7PoGZEO5Y8aMoYULF2pSBw8dOsQrEI4dO9bkdQDcyfQfd9Kua9pu8W/c2cHk7+0gORH3doWhHTMVqrRXn8d0lc8PdQURwQH0fJ84g88ZS8E/cKOMXvt1L+26pqKv992UjWRN//Ok3U8kARxdmeR/S2WnSrDvb5fPyxGp3HVoAvQaQPdrrG1+PfEX4dh+98cb9V6pIH9fig0XsszEES7JNQW3YnZK4bZt23jJdnb7888/zdrGnj17eDNj1kiZmTdvHr+/cuXKOq0zd+5cXomQVUhklQgHDhxIjzzyCC/kUZd1ANwBqy649LA2L59pHB1q8vcPSG1M3z+YKmyrmujn7afcOqUwJdqHerZKIFckDa5NtfxUoWZ54abjfOSL3W6qJFfuXajICIAlSQeRHO3/ZFw3w/M6wf38d1BL2f3FW47T0SztRcgTs4bofY+feg6Xu8btZpecUSgUtHjxYnr33Xf56NFrr72mqVpoKlaIolWrVnxZWk4+KiqqTuvExMTQwYMHedGOjIwMeuaZZ6h79+6yn2XKOgDuoKBEO7LFPNAhnALrmIfeuzULMITqVdNXnKexvYX/UXdSVCqcDDWL0l7pczXSv7M53lx3xeDjrlTVEcBSWEqWtKK2PQKuSsmotNRdbUJoVJfaezWCe2geJy8RP2OVvFejobltfuIIFwIu87z00ku8V9b9999PR4/W7YOZBVJiMFWfdcT5VqNHj673OgCurlAy94hpHWf66JahrvPumh5QVCqk+wT5uW5vGvZ3ZoXJLD11Y9JPR+n4/zWy7EYBnFyJToBVKul7aCsFOjHeHa1DaPQtidSvLf5fQZ5yvu6Z7rRo2xn64aB+PQdDfMWiGW56zmBWSqGXl/wEY/jw4Tz1r64BFwDY3ndb5MVjWsbLr1SZas5twrylaH/3nMhVrJ5sEazQlr51RabObw6oQ75EcblQtRDAlljVv8e/2kw/bpMfAx1FsUq3KqDtA64snaCvUUQADWzXRHaRDYBpHh9Jr9/TTe/F6C+Z3yXlK5nD5Y7Hf7P+gyoq9Ie5W7RoQfv376+xkTAA2M9rP+/i5bsX7s+WPZ7aOMas7SU3EEbGWJNDVoTDXQ6gbA7c099t5dX23CHgqmkE88CrA8hfff2tbxNtqWBT5BWrSxwC2MjQd9fTuvPFNG2FUFXN0RRLCvEwpXYIuKo9qig1QvtP37WZeZ8P4B7ELAjm1uZB9Pm9rWnBo70NruunnsPFlLnhRC6Ltg0PCAigZs2Q4wvgaM7fyDU67K/wNe8wEOinDTRYEY5/Tq6if2cOI18fix5WHM6iLSfo79PawhAh/q4dcOmmF/15soAvs55dLK1kw5S+dDY9l9Jyimj1uYsGvy/UlyhffvGeMvJK+PcD2EJlZRVdK3Lsk7xinb5bZXZIKVRRJXWMIjqao3+cBzBk4/O96WJmPvWtJe3UT3JuUFpeYfa5h7My+bdNSUnhX1kZeHHZGHuXigeAmq+cWkKgzshOQRnR5+uO0jMjbnHpl/9shhBwiJrEuE9/vkFt4jQBV4i/H/8aHxnMbzWlafl5sbRT+VDZpEX7aMGDnejT9SfpP31bUKdmhkvQA1hCVkGJ3kg16z3nSE5ek18UK7XxKEBldTXlV1aQr2TWiL+bnRRD3SVGh/JbbXzVfbjEwkmhdUuKcHom/yc99thjBpcBwPEV6cwNEK34bxeztxnop1+F6GS6PBhxRYXqhtEeVE0f3dWSerZsSO4iuUGI5ndn/VWkvDwNz+XrmaCgy3llegHXpYIquu3TfXx5+6UDdPjN26y23wAZecWyFyErv5jiIoId6oWZvvKiXVMK88vK+H+pmCLGBGCECyzEW3KB4+jlLE1/LndhcsD13HPPGVwGAOcc4eqdqKC2jczPz9c94Waqql2//FC+OuB6olsM3d61Obm6XokK2nFVRYOTA/j75bN7W1F4oEKWHsKkNIpkyat8uXGIJ/0+uR+tP3KFRnVpSoPfWVfjz9BNNwSwtMz8Er0AzNECLl22nueSWyq0DPH3ZdkLwnEOKYVgKZ6e2qq39mrqbU9mjafn5ubSt99+q7n/448/UufOnWns2LGUnS2fkA8A9ldoYITrvXH61YXqQiGZACuqVJe023nqGr2waDst23WaF+twtAae9VGgDrjCA4WUOlc3f3xP+r/hTeijh3rx+7d2bErdDYzqtUmM1iyzwa7IkADen41dIb9erA3EFa5bRR8c1OXMfHr8pxOyxzadSCNHV27jhkW5ZULAFRKiDUT9McIFFtQ4xMttezGaFXBNnTqVPDyE9JGsrCyeYjhgwADKzMykF1980dL7CABmSs8ppJk/76LT6fmyx+9uG1rv4Xx2tUrXmSyh8twry47SsuP59MLyc7xYx/xVh8hV5KuEK3MRbhJwRQT700P92xoc0dRbVyF8HdqKjXZpxQe5Z+sAcAy3zd+u99i87Tf0yrDbW4cGPnYLuMorK6mwXMiECAsOqfHCGoC5fL09DPaccwdmBVx//fWXpoHwqlWrqHfv3vTuu+/SDz/8wO8DgGOYuHAXLTqYQ1/syeL3YwI86JWBCTR7bHer/DyVuqPhlUL5icK5rCJyFTdLhN8tyY2KZZjq14k9acbgRHppdGfZ4/+7Q1toaURLw6/b7tNpNOePvVRW7n4fxGBdRUZqBr25TJhD6CiOZwk7KvZSL7NhwJWrEi6WBXh6Usv4iBovrAGYS6EunIERLhMplUqqVs/V2LhxIw0ZMoQvh4WFUUmJPE8aAOzncKb8TCOrpIomDm1vsXKsfRqphzTUWEYhSyfUxXrfHLmUSa5AWaEd+QG5prHh9NjgdnrV36JDtO+TSYPb8PmDuu777hC/MPDFuqN4WcFidPsDtorQHvvWnsmj//2+l7J1KhjaS4U6JTtQvYvlNpzDlaMS+gqGe3tTi7hwev/2ZrTwwVSb/XxwD37qES5lGeZwmaRbt270wgsv8BGt3377jUaOHMkf//fff6lLF/OrngGAdTVS509byhcT+tILfePpozuF4hHKimq6oVMNjGHnEaM/d6yryeZgc9HUg3gUYkKKHYivlTb9snF0KL14a1ujL83if9PxsoHF/PfbrbL7Syb20SznlhJ9uTeLnvp+p0O84uJUjeZRfrYf4VIKI1xhXkK0d3ePltQ/pea+SgDmjnCp3LBohlmXuT/55BOaMGECbdq0iV577TVq1aoVf/ztt9+mV1991dL7CAAW8uIQy1bVYwURnh5xC13KyGMdqqi4guizLYab37qC/BLhpIQJDXCPOVyWwHq0/G9EEoUF+vHG2B2SY+n1oY3pjbWX9dbNKKmmH7edoNziMmoeG0pDOiTZZZ/B+aVlF9Dqc/ILQKyYi67daUKxCHs6k5bNL0wxIQp2alZK5ZW2qfqqqqigkgr1/C1v9N0C6/FDwFU3LVq0oG3btuk9vnr1aov9UQDA8sKDrBMkhAVpU8TO5rrulauCEuHEjJW2ZYEDmG5c3zay+48OTKHoYAVNXnZab91pK8SgPY2WhwZQ+6QGeKmhzs5cz5XdX/xwO/71k3ta0uTf9N939jT6092a5RA/b01zZlu4qRRSKkN9fMlbPcoGYA1+3uoqhW44woXZkABueLCztGCFe6TXnbvBRvJw4LSUEZ2SeRPlmkxast9iPw/cy/kMbXXWuEAP6t06kS+P7NxMU1HTUaiLn3L+vsKpma1GuHKUwvytCD+M2oN1+avf2yrWjMvNIOACcCO+Vgq4vHSKJLiqZ347xb+Wut9nhVWwCmgKr5qvqF8vcv1m2mAdO89r+4LqvstytNnB3C87TtGhy0I1V1uorKyiRZuP09nr+r1LfdTH03Ixx9CKWLN6MeCKVCDgAutSqM9BVLUUhMkpLKFvNxyl/GKdf1Qn5h5nSQBuyFCPGT8r9lQZ2jRQs/xy/4Y0qXsMzburhWyd137ZRRdu5NLiLcep1AnLf7vhRTmr49NVAKxg40Vt9cEuidpmvoa8/Pd5umuh7SpkLtl2kmauvkQjPtmtV0nRVz3PxRYphfkqFVVWV5OPpycFecv7gAFYmr+vcA5SaiTgYnMZ3/h1N3WcvYneXHeFBryzwWX+CPioA3BBJaXl1O/t9TYNuHo3j6a154UJ6iH+PvTf4R34ckJEEN399QG+/MOBHPrhgFAR7GZhKT03sqPV9gecg7+3B+WW1nwlv0hZZlLjZQBRTqEwasPE+HvQuw/0MOnFKa2sIls0fPh+91XNRZwHP91keITLBimF2eLoln8Ar5KI8WSwJoX6HMRYSuGEhXvpqqSPZ04p8bYNhordOBuMcAG4oHkrD9JNpf5Hp8KKhR6igrWTIsIkFfw6NYszuP5H29Pp2GXn6s2VGi1cAR6c7PwHf0fhzyqQqIX4EnWN109revDzzWZte92hi/yKKbifnafTNMtThzbVK3Lzw0OpBt9r10utP/LORq7O52knbe28prJbSqE24EJfQbBdwFVmZPRWGmyJcopcI63Q5LOvlJQUkzd67Ngxc/cHACzgn+OGTzKtOcIVE6oNQkID5KMR/+kcRd/+e1Pve0Z+to8u/O9WPpfHmfRoGmXvXXAZPpI5XK2i/ChAEoCJDmWU89YD32w+Sff1aEptG8XUut3959Lp8Z9O8OVLc2+TpdpuPHqFBrVrzNsagGuSViEMD9QPrPq2bUTN48KpxzvaPl1B3tWUVqGk+IpgCrRiAYmSUv10bxFrNnziWo5NRrhYOfiicmFfIljAVeF8ad7gXAJ8hbBDVVHN52dtPn6Vht+SRH4+3nzeliEqJ5x+UK+A67HHHrPungCAxZQZ+aBmBzVriZYM+YcGykuAvXZPN/ru33+oWm/qOtHPO07T/X1akzMQT4DEXiJQf6eytR+m93VNpMW7rxhcr/+HO/jXRQdzZAGUMb/uvaBZZvMFxff+M9/voA0XS+jWI2n02YR+FvgNwNGwghRSxhr46h6n2MASG3c6mZdHnQKEFDtrUJUZLondMMiT76tYDbXCyiNcYrGMEF8/8vXyoioEXGBlAeqWByyl8KmFO2j7VRXdeyaT3nmwF/V+W55aK50i4QpMPvt67rnnrLsnAGARJaoyyjKQTmj1gEsywhWk0B85eKxLNH21T3+Ua/2JG04TcBWr886lKZNgOfHhQRTH8gpvCFfdn+wWTV/sqXvluIISFf10RDhpZU5dy9b08mLBFrPqbJEsEAPXkS2Zv8UYG0HXHeH09iBiOQD55WV0rbCQEkNCrLJ/yjLDV+yHtAiXVZO1dsCVre6/hXRCsJVGUULxmpvKKrpUIKQK/nIsj14tVlGJkYEspZELFM4Gl2kBXMys3/ZplgO85ScT3lYs385OXlhlwgc7hFNyrHDiIBVspOiBs/TZVJVVaPLLrdVA2h0tUTej7dbQj7o2j6OZd3aiAU0C6Kv72tCztwqFV3SV1ZJiotvwdsIP++mbDUf1SnDP+FnbbBZcx6ZjQkEKJrQOtVbY3KkkhTDqdT43h6fcWYOxCq1TRnawWdEMoRy8cMKLgAtspXl8BP9apDNo9fk64xVCiwxUXHaLOVxsflZt87kwhwvAftjVItHEnvH0wdbrfFldjdWqpt7RxehzDcMNF5rYn2Y4b9vRTPtpl2Y5KSbMrvviSnq1TqRLc4WGtExseBB9N3FAjd+TV1xKMWH6H1+svPbvu8/SG6vOyR5nBWT+b90VInaT+PVYPr1b798AHAlrOzH1n/Oa+2+OlLem0MVaWLyzWSiwweKceB9fyqqs4qNcZ3KyqV2MMDJqSapy/Sv2Sx9pTyEBCllZ+EorjnDlqlRUUV1Fvp5eFIKGx2AjYYEKCvYhKtQJuD7bbbyA1rTlp3mzcrecw4X5XADO4b6eLTQBVy39Za2uRbx21OuNYY3p9TWX+XKYwjkG2v84UaBZjg3X9hwD62IjsxU65515xSqKCdP/G3zwzwH6ZGcG/iRu7OgVbdpyQrAn3d61eY3rs/YVYsDVINCbz9tqERpK/97MoqySEj7PiReUsCD2/hWxkd3/DmhBPVslaB7z06QUktXcLBFaeEQF+FttrhqAIWEKTyqsQ1NLFpxl5hUbPOa7/BwuzOdyz6uGoQF+LtELwZ1I51J52vkztW1iFI1oHsT7vDzUrw1FBilo8rLTVj2psBZnq6rozL4f344+XHuKZo5uR+O/+5fyy4h2n71BLRpGatILT1/P4e8vBFvAPqdEZbqRuhFf39+GPt98jt4a1Zbo2nkK8vGhhsEhdK2wgI9ydY1vSJ4WDErEMteNQzzp56cH6z3v62PdOVzV1dU8mGSiA5z7JBacT5Av+/ys2wf/XZ9spe0zbiVnVq/ZwmfPnqWTJ0/y5datW1Pz5jVfSQLndCUrnwZ+tJPnwv/8RDdqER+BE04H1iTEky4VVFGbKG9ZQQB7j3CxIOVTSVU4cT6ZNdNmwDVSDtmN8ffxoPyyanpt9SVqERdGnZvF0eiPNvBKh2zU1BzrD1+kwe2TLLzXYC/lFdp0vX5NTSt6wf7+7FalVFLxNeGx5LAwyiguouLyckqzcAGN3OIyyYmnPm3RDLKKwrIyKq2s5EFkuHrOGoCtBBl539fkWpETXpnVYdZl2qysLBo5ciS1aNGC7r77bn5jy6NHj6abN/WrkIFzW3dEmPfAriwP/2QPfbjioL13CWpwo1g4MD03uAV5SYpkeNp7iEuHOE/BxIvQDmNce/2CIGAbvpKrBptOpNHEb7ZoysqLKap19cvey5SWrU0XBedWKgm4XlQXoTCHj5cXJYcJ/+sX8nKpvNJyldJyi0v51xCFV4296aw1wnVTPbrFimV4YbQebCzQhAnl97VzvXnSZgVcjz/+OA+69u3bRyqVit/YckZGBj3xxBOW30uwq1KdCb6f7bxht32Bmm07foVU6j9XiE5VQAeLtzSVuJxhhIsVYxBfvycGtrH37rgt1ixTlFtcTusv1Fxw5fWhtY96rT1fTP3f30bXswstso9gX2XqgCsp1JMahAXVa1sNg4MpyMeXKqqqeNBlqWPJ2cwivhyiMJxkJGYmWOtilJhOGOWPKQJge4F+xgMulpkz57Zkmjuul8sFXWYFXGvXrqWlS5dS586dycvLi9/Y8pIlS2j16tWW30uwq2u58pMaxz89dl8PLdKWVg0JkAdc9k4p1OWjTpuxYuVjiykprSAxLgwPQgqOvZRKzkCl1TiNeXRgCr1/u1DdKj7Ig7rECf8Tc25LovEdhfLEDJu/fTJNXjIenFO5+j3iZYE5V7yARoTwPmF9uYrK6l+ees4f+3gPuJoCLmumFCoryqmoXPg9ogIQcIHtBRt53zMrXxym6cv50shbZM9dyqj9mO9yc7hiYmIoMFB/omVQUBA1aGD5EqpgP3P+2EtLD8uv7LETZFZliZX3BOeYPM5YctK3ZUe4yOH9sOVErf3EwPpC/Dwpv6z2Nwwb3ZjcP5kv392jJb/p+nT1Idn9QqVr9Hpxd1nqpsdBfpYpbBPu70/RAQF8VIgV0LilQWy9qvpJm7/HhvrXnG5thWNjZrFwATXMT0G+XjboFQKgI8jPtNCDFWm7NPc2avLKCn7/THouNWngvKNeZh2Rxo0bR08//TTl5ORoHmPL7LH777/fkvsHdlRQoqIv9mQZfK7H/zbwMp3guOLCgxx6hEusxOUMI1xi2WgGFQrt57mBQhBVm03TbjUYZEn1ahknu59XgoDLFby3RWiDkaO03Jyr5uER5EkevHfVTaXl+gamJmpHWaXEizrs2oKxJsnmyiwWPrdjDFw0B7CFBpILDQ92qH1OdPNw4VzhQmaBe6YU/vbbbxQfH09t27alNm3a8GX2GHuue/fumhs4n8rKKh5s9Xt7g9F12GfZoUvGG9WBneYZqZe3TOmtFxhI6mc4BGca4eoYi1EtR8CCqA3P9ZQ9tu2F3rL7o1oGm7Qt3V5qBUqdTpzgsEpUZfx4p6tCcjDJsGBVM38fH2oUKlQpPJuTQ1XV5l8lig3UXvkalGp4jqE0e2TQO+vIUpTl5VRQJhTsQMAF9hITog24BrSJp54JNWdL5SqF/+W5G6/R/nPp5FYphXfccQe/gWsa/eE6On6z9qtqT/x0gi51QDllR6Eqq9B0tjCU9uZoKYXiPAVnGOEK8BGCwye7Rdt7V9xe01j5FdHE6FDZ/f8ONq2oSXRIAPVI8KNd14QT0AKkFDqF9JxC6vfeVmrfwJd+fXaI7LnLWfma5c/ua2vRn9s4NIyuFxWRsqKCrhbk8/vmiA3yphvF5TR9UIKsiqyUv6/21OxaYRXlF6so1AIp/BnqZsesFLwf0gnBToIl/UFZP8437upAD3+zh+7vJM86MFTwa+up69SpmeH1XDLgmjFjhuX3BBwCu2poSrAlHQ0z9qEBtiUdcQzRmb/FeDtYmULNCFc9rhbb6n9i+1WVQwat7oq9laXFLdmowY1i4YHWiVGmbcPTk36cPJie/m4r/X26kHJL7DPCxRo3+0r65UHNlu44w1Pt9qXrp4C+889hzfLAdk0s+lJ6e3pSs/BwOnHzJl3My6PYwCDy8/Y2u9JmgCSo0qWbnXA2PZf3nLNUOmEDpBOCHQVKAq6IYH9qFB1KO2eOMLr+8wOTaNqKC3w5o0D4LHZGOFMGmYuZ2iuEusRqX1L5Jc775nc1X246o9dUmGFXUsP9iF4fbdkrvpaaw1VNHrJUIEdTpNKe2F28iXmLjuCviV0oQkH0VA+hSNO8se35e3zqgIZmlyhedjyfjyTY0omrWdRu1hp6YdF2m/5cZ1bTJY8156z7/8mCrBBfP36R6LyZZeJLWElMnVGs2ny/7axZP0v+c8t5w2P2+kUHYP4W2E9qo2hKDPbkJeAbRtSeAs6qFo5JETIZfjqSR8evOOd0FrMuqzVrpn/iLXXu3Dlz9wfs7MRVw6WRk0M9aVTnpvTCcvnfNrtQRRHBKC3rCFIahtDmy0KFLqknhrTnN0fDGotKr/I7as3LIsncnnJnyH90AymNY+jArNs097u1aEgH36h7sKX7Nz1yKZP6tG1EtjJ/7QneN48Fe+/b7Ke62LxVGzbuFcvE/3sjndKLiig+KJjCFKYfuS5n5tOVAiHgCqylUttdbULo9xNCkQA2Aju/nvt+o0goRR+u8Ed1QrCrAIUvbXv11jp9j7+kWfKHq4/T10/EkLMx60j14osvym5TpkyhUaNGUXp6Oo0ZM8byewk2k1Wgf8L+TK9YWj91GE97Yb1rArWjwTTk4110McMyDSGhficeKnWD6gFNnCMA9lOPcDE3HLjiZYFSmOPDPDfcsUYJof6kBRBYH7vZy/bor2OgQEO9fmZVFb9tviCcBDOzftlt8Z/jiqSXPMrUddPF102sBM+unFtLqEJBcUFCBdgTN7Oosg5/s52nhQqKTMfk2BrX/WB8H/L3kldpM1d1dTUPEBlx3wGcSX6JdqpLqCQl0ZmYdVSaOHGiwcf79etH3377bX33CezozXVXZPe/uq8NDZEUxnjz3h40654q6v7WasosET765q85xj8cwH7FMoa+t05z5TTU3znmg0hTan7bc55eHOqYwUyRShjhCvYhSm3sfFfVoGa658usT9Krd8sr4g1+bz21iFTQwkkD6/1y/rLjFL3893m+rK7Fwi08kE1/Hl9Fv0/qSck6hUFAS5p+rCqvoAVrDtPCvTfou4c7UsNgT7qQX0X/6WnZ+VuGysTnKJW8gAZLLWwWWHsQ8/HKg/TBViHgah/jQzFhtaf1PdY9lubvuEFR0qucZmDl7FWVFeTt4cl7igE4mxL1BWWmQjqB14lYdCx+yJAhtH07ctFdSZSkfKeIpXAE+mgz6Z3zre86dpy6pgm2mCD1nBRHx95HEeraHvkOXCFObIir8EbBDFdUWqHfr6nj6yto3ooDfPmf/RfoelG1wXRdc0z9W5uWrZ7Oo5FXSjR7ubwhMxgPuH7ffZYHJIXlRP/311EqURekCAnwtXo6dKtIoTjL1YICylbVPPcvK79YE2wxpgZQIf7CesVl9espll5UqCmW4WXDFEwAS3lmqLb6bIH6Iqizseil8G3btpG/v+HO6cbk5+fTokWL6NSpU/Tss89S8+bN9dZhqYqLFy+mjIwMSk1N5Y2XfXy0B6yff/6Z/2wp1hds+vTpddqOuyuWFAdgIhVELeINN2bMUfdFYP44UUAfWn3vwJgfd1+U3b+9k3Wv7lrSgKahtOxEPqlY2TEHJY5wIeByTZMGtabV5/6VPZZTSvThtnR69jbWJ057SYk1ofWrZ0VBVumypqmApeo0OTDy+kgCLmlGBntNS8qFFzbUygEXExUQQA2DgymtsJBO5OVSx8BAMjZmtXT7adn9Aa1MGykPVbf3KK7H8bGiqooyS4RmzXHBSCcE55TaOIam9I3nFy7y2cRXJ2TWpY6RI0fq3Xr16sW/PvXUUyZv58svv6TWrVvTrl27aMGCBZSWlqa3zpkzZ3hwtGnTJgoODqbZs2fT0KFDqbJS+4Kz53bs2EGtWrXS3JKSkuq8HXeXlS8clJn5d7ekna8Oo0CF4Q+uYue8wOCSciS5zVEKD2prYllsR6DwFQ5BSkm6gKPJLxHmcPljhMsltU8SKh2aIs8CVQxrGydVSPMMQebZhVvphwM5Bl+VC7nlJFb2jwqxTdpci4hICvXzo4rqajpSXExKAyNdrHUKC95Fc25Lpgf7mZY+HRYopACcz6ukOz5cy7dlTrEMNk8x0MeHV1gEcFbRwUKBmgM3ymjfWe2IsdtVKQwPD6f/+7//o4EDTc9x79GjB509e5Zyc3Np6dKlBteZOnUqpaSk0IoVK3iFoAkTJlDTpk1pyZIlNH78eM167LHJkycb/VmmbsedHb6Uxb+yz/tRXWquRPl4txj6bHem3apFgZZSnZf03x4N6JlbO5CiDuWG7c1fXThD/B0cTUZeEc1YdYkvhyicI1UT6o79aY1dNBWL0TDd5m7hBQyWPzeIAvzMy45oHeVDR7PKa/2fAH3LTwmpcaIYfw/KVAqjWqoK0jR+jzVhfpQlsNHK1OgY2p9+nZSVlbR7/37q1LkzhYVpmyKfSrsp+57bOskvBtdEOlJ3KKOczt/IpRYNI+tULONqoVDpMD44mJ/7ADiruHDt//XU3w7Txmnx5EzMOjP76KOPLPLD2YgTwwIuQ8rLy2nVqlX08ccfaw4UCQkJ1L9/f1q+fLksUGKl6F9++WUKDQ2lPn36UN++fc3ajjtbfuga/2pKH9qpd3Sh+3vlU9/3t2uu/KI8vH2o1ClITWOCnSrYkjb/FH8HRzNv1RG9/jngPliZ+DfWXpY9dja3kvadTad+KeaVj48J8iHSCbiWPdaR7v76gKYE+N+vrKDP721Fwzs2rcfeuz6les4WI/3vNJaZYQ2s+XGHyCg6dDOL1l5U0fNbttNDnRrQS3d04c8XStpKsD99SIDpZeQjguTrnrlet4ArR6Xi/be8PDx4CXsAZxYvCbgc9JShRvUakmCjU3/99Re/sWVLu3LlCpWWluqlByYnJ8t+HguiGjRowK8qXb9+nYYPH06TJk2q83Z0se8pKCiQ3VzZ1ktCSmGPRNM+EFh3cLE+Q3qutrwx2Fa+SjjyRAbXbf6kI/DyFC6A7LqmoisGUortrUQyWX1Y62i77gtYj7E5VW/9pQ24pQp15ruaW4Ze1DohklpFyC+WPPPbKbN/hisyVDKfFcvQ1aeR7Tv6Kby86JbAILpZ7ksF5R70+1EhW4TJKdKmGaoPdyaL1OlxeSFLPsJXm2sF+fxrXFAweSMDBZxcQmSwU1cqNCvgysrK4vO1WrRoQXfffTe/seXRo0fTzZvy4fP6UCqFqlBszpVUSEgIlagngTKsOMaaNWv4VzYXjI1mff7557Ru3bo6bUfXnDlz+IiZeEtMTCRXxUofixcLH+qZbPL3BfsKnyCZ+Zap4AV1/7vllgp/uBZxzldKevdF7ej2wROnHDYg1G28CK7lpQEJmmazUnuva3uwSU1dftrsnlmGgjs2Mq17PszqJFhizpirkKZ21iQm2D7zlHw8POg/fVrw5fTiaj7HhB2f//ur+cc1cQ6XKMNAn0xjisrK6Kb63CcxRP6+BnDWhsmippF+7hFwPf744zzo2rdvH6lYfweVii+z6n9PPPGExXYuSN2gLy8vT/Y4S0FkwZJINxBi/cDYY2KJelO3o2vatGm8iqJ4u3r1KrmqgxczNMv925oeWIYrhLfQvPXyKkxgG2fTc6maPPgclFjJcLuzeGqQcILCfHacnYw61lWrCsnZcUmp4xb2gPp5Ykh72vZCb3rvwV40OLn2ggusaND6I/JUQ1NVGbgyy+a/Pj9Y+78gevdvIc0QiP7eJ/Quq01koO3SCXW1idfO3Zr++xH6ddcZ2fNedRzhYu+LwzMHUZc44XdS1qE8/KV84XyH9d0KQDVmcBEPdhAuLPt5O1/NALP2eO3atbzIRefOncnLy4vf2DIrQLF69WqL7VyjRo14sHTy5EnZ4+x+mzbamvyGlJWVUUVFRb224+fnxwMy6c1VlaqvHob4EvnWoexxHPsG9YTefFyNtVmj43s/Xk9v/7mPLmUJaa5R/p5OWbSkd2ttcJ+lIipysKqh+ZL5F7d1bGzXfQHrSowO5f9D88b3Mvj87a2CZSfMT/503KLpi6zB/KdjWskeu5yNzAHR1H9MC7jCbFAS3pjQAD/ZXL89F7Jlz08ZUPeWHaGBCmodK1w03nnZcEoh6xl3z7x1VFIqHK+KK8opo7iYLyeFaoNAAGcXqm6VsP6C8ew0R2XWGVpMTAwFBupfTWdBDZtLZSnsw2/MmDG0cOFCTVrgoUOHaOfOnTR27Fh+nwVVGzZskH3fDz/8wEfbhg0bZvJ23J3Y/DNYXabbVB8+0F2zPOOXvRbfL9D33A/beaoTqxK5/5KQwhvm7/zpbmwSbFGVYwVceaz0GRE90jGSmsebPlkdnJexggudm0RSQrD2+MhGlnMKlRaZwyUa0akpzbtTO9KVo9S2fAAtVinSmEA/xygcxD5KV57Vzm1e9VQ3eqi/aeXgdYkj/wFGWlOwsvP/ppfRpxuFi8qXCoXALMo/gIL9nC/1CsCYrEJtmndOYYnrB1ysYfDTTz9NOTnafhhsmT12//33m7ydPXv28FLuM2fO5PfnzZvH769cuVKzzty5c3mVwU6dOvGfy8rOP/LIIzRq1ChNwYz333+funbtSo8++igNGjSIJk6cyOdfSSsV1rYdd1emrsBW11HayJAASggSvukCrsbaxJU87QFH7EkT5gIly9mVf0cb4WIjt0xKgvPNjwPzPdpJP7iOC9dPNfxgxUGzUwoHJgUYDB76pbjuXGFLiQ2Wl+Qf2VI7P7tLM8td9DXHooeE6su6vYpb16M/Yq/mwu90Ib+KysqNB+FZRWWUV1FBmep+YMnhGN0C1zI4JV5WtdOZeJubUrh//376+++/eS8r1uvhwoULvKofC2g2btyoWXf37t1Gt8MKUbAmxQz7PlFUVJRsNO3gwYO8KAYbtXrmmWeoe3ftqApLZ2QB2pEjR/h6rB/Y4sWLKS4uTvazatuOuxNHuHzqWkaJpXoMa05PLztN2UrHOll2J+EB5vUEcgTsxHPjxRIaluBYAZe0KEJCpJDSA+5hxl1d6bv9q2SPdUxuQGM75tM7m7XVNBcfyqX/9M+l5NjwOqcUdm4cTs8Pb0tJMfKT4rBABT3SKZIW7s8263jsioqU8qqQQTqjWC+PbE8j22dTg7AAapNo32qiLQ2UbV89uVu9thmk0B7f/9l/nu7q3tLo5/h5lTDq2jAomILR6BhczNAOSaT49QTvm3juRh51b9mQXDrguuOOO/itvliwJQZcNWFzqVgFxJq0a9eO3+q7HXdVpg64vOs6q5cFs6FCOXJluWMVPHBVxbqXTu08Uby+PvtPH9p98grlXz5BxVVVPOXK00HmyolaxEfYdV/Atry8PCkp1JMu5gv/a+ue6cH7DD45pB2lJEbQ9D+O07Ui4bmxX+yifa+PqHNKITvWpjaOMbhOcjQbscmmw5nltHL/eZ5q6M6KJGX42WB+gE7FUFYumrUpcQRRIfL2HG2jvKlVgvmjW0xSjPZ3W7LnitGAq6CoiNjb0tvDg5LDMSoPrqlbgj9tuaykizedqx2RWQHXjBkzLL8nYFel6i5y5lxRFScKq6e7gJUVGgi4ghSOMW/BHH4+3tQ3NYnWXj9D5RUVfMJ3KNm+l46uYvUEdCZQcoUZ3MPfzw2in3eepsGpjaix+oSXBWJ92zaiButP07UiIQjw9/EwK03Vv4biRH7e2oCClRW/5OYBl1gMgtn5ykD6aNUhzf0hTQMdqmCQ7r60VRe8qG9BF9H+9DJeoIoV09CjvubZIjSUfL2cP80cwJBH+zSj29qV2D19uK4c5ygFDjHC5WPGCBdLgWFKq4iOXMq0+L6BnPo8T6ZbU8NXyp0Fm4sZou6TV1huoJupHaw7rC377eOFQ6W7CfL3pQmDUjXBlrF2Aa2i/c1KU72SbfzqrJ+Pc73fHlywgZq8soJGvLeGdpy8arWAi41uRQT7U5Cf9gJIYrhjN3wP9rf8xTBpgaq0bKFSLcPeldHePhSjcOzXBKA++qc0ont7taKkBs41iutcR3WwmnL1CUR9Ai7msR/+teh+gRxrpGlggIv6tG3k9C9VqIMFXDNWXtQsO9IVdLC/ke1iNcvllaY3QFZK0lRLTGzkG+7gReZYELn9qlCk4cTNCnrg+yMW/xni6yZmEkrncDliP56P7myuV8a6vqRNuYskPQFf/eVf2Zz2lv7+/AIWADgWxztSgV37cJmTUqjw1X74ZZZU05Nfb6YLN5yreoyzePP3feSqxBGuIgcJuEw/jQZ3858BKdRZ3Yx20yUl3f/JBsrKF/oeGfPV+iP0mmRkYkBreWEnqfZNtIUf1NneDktlYuBYHyWl6oBL/fkULEnxlaZfOgrpnKsQf8ukI8+6pwuJVeE3XSrRjJZey9emPPx3SHvyQrAF4JAQcIHsKq05I1xiHr1ozblievufw3hlLWzZrtP005E8vmzmn8lJAq4KXvnU3sRy3ZO6O3e6Jlgem8vVOEKbtrXrmoq+3yL0QNK1+sAFemnxdpq9/iotO5GveXxw+ySj22epMm+PFOZtFZVXU2UdRtFsTVVm/QskW06l868lFcJxIVAyZzXcAQsGSRsgWyrgCglQ0FO9tCOrYvq+tNBVx0bOlWIF4E4QcAFXJhbNMHOuSps4bR8UMegCy3ph+TlZpTMxEIgJcI3oKzAggFiIU0XVVOwAo1xnc4Ur90m8YhyA4WOm6JOdGfTucm16FxuBmPXrbpr4y0n69Zg20GICTJjWM6pLU22D5aK6N1i2lWKV9f9Xv/1XaPBeos7IDJGk6YUFOl7OZVSItmebJa8dPTlY6PHFfLxOCPATwxzv9wcAfQi4gFOq00L8zZysbSitJKfQcU8SnM3u09reP+JJ2AujOtGc25Loj6d6kytg8w4C1ZW1Csu0zZ3tXRK+puIG4L6GSRpwihbsytCMRn2y+jDvpWWITlVzgwL8fMhPfTjOKighR1VkIOCSFgepr1JJo18xUA2WBFyRwfavaGqo4IqoWazlmg8HKHwpxl+4wMZ6F7ImyDHBws9qFem8lWoB3AECLpAFTAof8/Lhi9U59lIbj17Bq2sh932nLYMs8vbypPv7tKGGkdrJ1M4uyFN4/xWVGSjFaEObj2krFLZpiDQd0DeySzNNQCTFGtO+snQHfbD1utGXzdfEnOAg9Xl7Vr7jXrwyNPr2977zFtt+saQH18+Pd9ELuKKCHbMi37LHOvLiGe2TLFu6uqJKO2TWYuYaWnJImC/dKkY7qgYAjgcBF/DJ3stPFfJXwt/MgKtdon5j2Bf/0qbAgflYzxVdg5Jc88M1SDPCZd+A62ah9jUffovxuTbg3g7PGka3tZD3WXr29zOauZbG+IvVD2oR7Ct8RP+27xI5qpwi/dHoE9dz6c3f9tDZ6zkWS1n0oGpqmyg0EI6UBFlJDSw3gmRJnZrF0R3dWlh8u95GCluZO/8aAGwDY9BA037WVs7yNyXXxYB7erSgAmUZvbXe8j1Y3N2mY9qRwplDGlF+SRnd19PyH+SOFnDZs3BGbrEQ8HWJ80VJeKixQuuUEe1oxZmddXqVTE3dzmfNDdmI0elCmldV5ZDvxdxi/YDriz1Z/Os/J7Jp72sjLFISngUa4u/fICyIjx4FsrTLGhpIu6L+TUPpl2P6AX2ohYpzAIB1ON7RG2zuULo2JSRAUuK9LtgH4WOD29HUAQ0tuGcgnSPBmn6yRqxTRnWi+EjXLOQQ6OlJ7DptRVUVlVZav9y0MUfThCIHIZJqaACGNI0N53Mp6yLc37QLW9JLDo46j0sMuLrG+1HfxvL0PtYmpL7+PZ/Bv0p6TXNs9GhIB/cbfX7l9lsMPh4W4HjVGgFACwEXkPRiq7kjXKInBrfDK2pBbPL54j3X+HKLSNe/gunp4UEB3kKQU1hhv0qF1/KFk0ik6YAp2FzKumgYZlqhB29JT6WCEvsWkjEmr0QYDQ5ReNHQNpadr8RMVzcgl0xdcmsRwYbTySMcsFojAGgh4AJSSOYTRAT51bs/DVjO/y3bS6dyhJSa4jLH7cVjScE+PnZvgFyuvpzeNl7bwBTA3B6FA5oEyPq5tUswrRBLj8ba+WE5RfpzOR3BoWuFmpS2tomRsufiAus3r+joZaHXFNQuLlzbCxMAHA/OjoEuFwgn8tH+HnRbJ6H3S328MayxTV/Vv/edoz/3nCFXHN36TlJW+nye/VLsbClIE3DpV760lVyV8D+RYqAYDEBdvHFXJ/pu4gCaekcX+v7BVPpvjwb0YN/WJn3v/93bVbO864yQWudoLucJI1wNQhQUESRPKaxvGYc/9gmjW/b4XHE2/dom2nsXAKAGCLjc3Nnr2ZpUjSkDkywyAflWSVW3L9cdlpX1tUaFxaeXnabn/jhL17OFK62uIs/AZHR3EOTtY9eUQhbo5pcK/xS4agz1FRakTR/sl9KIXr69s8nFL8ICFZQUKqyb46DHgyL1v2m/1nEUHSIPuK4Xm58HyPqZiQ2PmYf61S1t053EBmgLigCAY8J/qJu7nFWgWb69S/1Ht3SbPv5vwzWa+uNuspbM/GLN8mfrj9GZNMONRp2R9HdjWrtJY0txhIsVzSivtn0aZVp2IZVXCWWoHbXkNDie1ZO7GW1gXB+dE4UCOcuP52iaKjsKdnFCzHROiAzmjXlfH9pYlj7JLuqZ4+Q1bbD1+b2tEFBILHwwlW5trk03TYl1zF5kAKCFgMvNXcgUAq6mYV78w9IS2AlGsOQcY/056408FSq1oyCLDubQbQusF9xZ+kTl992nadmu07T9pOFS+kt3npXd//zh7uQOfDw9SaEunFFkhxPMtJwi/jXE18PtSk6D+VolCD2iLE3sjZhfRvS/P/eRIymX/H/6egv7+ejAFJ4+KbqWLfw/1dXTS/drlod3tMzFQFfRP6URfTahH/32WEe6NyWM3rnfcLAPAI4DZxNubvs54epjmIllik0VF+RFhbnCnCOVFaceFemkK7KRiUsZedTEwUcm/jp4iab8o21meuDVKIqQNPPMzCumHw5om4aem30rebtRQZJgX19SVVRQkR1Kw2fkCeW3Q/3c5/UGy1j0UCr9sf8yTR19C8396yD1aBpd721WSvrRfbsvi2beTQ5Dpe6Rxfj7ykfyGod48vnBBeoqhsyFG7n0574L1LVZDPVubXzOUUVlFV3Md6zRPEfUuVkcvwGA48MZhRsrK6+gvdeEHlzRgZYtOf7y8Jay+9lW6iEj9qiSuvfzujUhtYd524RS76LsQvnrc+WmNtWzUYinWwVbYsDFFFXZPuBKzy+pU68kAFGfto3og/F9eGPeD8f3oXt7tar3i1NeoQ24qsmDCkocp1phabn2/9PX29Ngc+e9F4UmyMyI+Tvp4x036MHvj9CP204a3e7TC7dZZX8BAOzFvc7iQObrDcc0o0+dm1i2Gtvg9kkUIBk/ffCLbTYLuDKVjt2wpby6mi7nyyvw/bX/kuxEasMxbUD2xUPa9Bx3EeQrtCcotsMI15VsYe5cXAgaiYL9RQfLW3VkFWgb1dubSl1JlM139NVJvz2rbmex5FAuHwnLK1bJsh2mrbhA+8+l621z6/ErtOqsNg1ROh8MAMBZIeByY6uOa3uchFuhaeLdqdog7mS2dUp8F5Ua3m5OoeOclOg6p1JSkLc8KJy/4wY9+e0OvnwxI5c+2y38bVKifah1onXmhjjDCFdxVRVV2jjouporvHcaRRhuMApgS/8d1k5WIMGRjm3iCJe3p34B+OQw7Qjx3OX/0gOfbdFb58iVbJ4+yFp7/PPvOb68/fQNzfMKL5LNBwMAcFYIuNxYlWRuQJRk/pClPD+ig+y+NSpsFUiKZhircOVIMpVKyiwvp/+0ImoXJT9J2XVNGOEaOV+bEhmgTstxN35eXuTjIfzuhcXyao3WVlgqnERG1rMJOIAlBCp8eYEE1ieRKTZykckeyiqEffEy0HDrjdtTNMsL92fT8Zv6+30xq5DeWf4vb+0x+bfT9NGKAxSs0Ka3zxjaxFq7DgBgU+55NgdUWl4h+wDs1TrB4q+KtAgE892mY3pzyHafTuNXNc1VoDTc4ys91zpzxuqDlTk/U5DPl2/t0IT+enGE3jpXs/KpWHJeMqSNe6bTeHh4UJA6Ramg0Lb91YrVda5D/C07rxGgPnzVUY2yzD796QxRlQkXJwxdF+rZKoGGNQus8ft/OJhDi/Zr53h9ty+D0vO1I3h3dW9hyd0FALAbBFxu6rkfhPQ15sV+8VYryjD+Fm1a4bKD6Xr7cN93h2jWr+aXci9UCdHJmJRQ+uRubaGOC1mO1QS5urqaTt7MovKqKgr09KRmSUJz6Gd7xxotlsGM613/SffO3o/LlgEXK9d/Lld4T4UGYIQLHIevOkOvxIFGuMSy8Dr1MjRax4bUug3pR09xOdGBa8L/+xNdo+vdwwwAwFEg4HJDrOS4dFLyY4NSrfazJgxoo1mu0qllsVK9D4sP5dY7/YuNRozs0oyGNBWuqP54IEMz32HPmTR+Im1Pl/PzKVup5P9wrfwDyMtT+Nd7dEBb2XpKSZllMZ3IXQWLAVeReX18zPHHnrO8EhwTEuC+rz04Hl91ZKJ7jLAnVQ1zuJhgA6PEbF7W1/drPxeahvvInjulnu+bFK2dtwYA4OwQcLmhY1e0KRzvjGpKCl/rtWNrHBNKcYHCh3FhWSVPHzx7Xej9ZQmF6qu9of7CyXG7hqH8a26p8Pzdn2ymsd8eoh+2nOBVsi5n5tOZNMv9fFPkqVR0IU8IKpuHhlKQl3YyeViggu5JEfaZef3vU5rlBzqEkzsL8hZOxAqLivgIoS0s2XNFs4wRLnAk4ihSeaXjVGE9oC75bmgOF/NAn9Z6jy1+tCMNTG1M3urvySrRBpDSKoa3dUq29O4CANgNAi43NGP5CZt+qH14b3v+tbismh79YhMN+Xg3Ldp8vN7bZaNWl3KF+Qyh6iupwzs00jy/7tBFTfPMT7ZcoQ7/t4H6fbCdhs7fzXvAsCac1lZaUUHHsjKJnSLFBgZSnL9+5bv3HuxNyaHCv2JakXYkbmz3ZuTOAry9+QGKVSkstlHhjObR2r9PmwT3qw4JjkscRaqw82i9IQVlhoNAdjFvfEd5y5GOyQ3I09OTxNaP14v0vzdSwUaYFdbZWQAAO0DA5YauF1fbNGVNTCvJLyPadkWoxDdrzaV6b/fDFQcpS91zK1Sd/hXsr5138/hP2sDypkr+oc56wAz8aCe9u/xfg9suVpVRWrZ8PlVdVVZV0ZHMTF4sI8DHh1pGRvFiEDWlC0m1dcNy8FLstQr0FEYDCwrq97cwVZl6TgobdfRys2bT4Ni81AFXeYXjBFxiemO3BONVbif016YPtor05sGW+HlgTASajgOAi8EZBVidNAgS6WbFmDPHivWukqbmMSHq1EJTLdglzPXSNeKDDdTnvW10Ki2HzC6SkX2TCspKydvTk9rHNOBfjblSIJ+X8UK/eJzws7RCdfplfr5Q3dHaitRzAoP8rJdmC1C/ES7HSSlUqYM/hY82TdpQWjlLXW8e7kUf3NdR83hwDfUwwv3x/wcArgUBF1hdQmRwrev837K99foZYeoRrgCFL594XR8s+LtcUMWLfGw8Ja+saKrzebmUUVzMyy+0i47hI1w1kUxjqHESursJUo8y2WqESwy4pL2AAByBeEwQKwM6AqW6aIailn6B9/ZqReumDqc2idGaxz6+V9unixkhae58d0fLtykBALAnBFxu7O422mIN1sRSSBKCa36rfbe/7oUsotSNQJmUxtp+Vadm31an7RQpy3gxDZU6PUZaBSxfVfeeN5fy8nhVQqZVZBSF+9etqXSTEE+6q5t7z98SBXlqR7hsUTijRN1XKFiBK+zgWLzVlSnq07fQ0krLhX3xM1YXvgYDUhtT60jt/1njSO38yQg0HQcAF4OAyw2J6fGP9Gth82aylhToI5yAfHRn8xr7iLFg78lu2iurulLeWMeLaXR8cw1vxlxQoi5xSERf/ZtD229q79fmSn4+H91imoVHUHxw7aN7UizNZvP0W6lBGEoiM4HqlMKysjIqLath0oeFlKhPIEMNpMEC2JOXh+OOcPnXkFJYkxUvDKEeCQrq3tCPmkt6doUHoWAGALgWBFxuprKyipTq0rthNmzs2jCk9hQtFuywYhWmUpZXG5239UwvbUPhuXel0LQ7u9LMIdoKhsbS+i5k5NFry/bLHv/1au37xEZfzufm0tlcYc5Xk9Awahxq+gjinNuSqEGAB30+rp3J3+MuJ5lBAQE2a4BcIr6nApBSCA46wuVAc7iK1SPC/ma2FmHZDz9OHkQ/PT2YooO1mQBRkmUAAFeAgMvNfLHuiGY5woYfanPv6cBTAEe1ND7i89KSXZQyay19tuaQSdtUVognx/qB45RRnejS3Nv4rXfrRP7YhEGptOeVfhQtSUXUdTEzn9adl5cgry2Traq6mk5lZ9Ol/Dx+PzksjN/q4v4+bWjPayOol3pfQStYPUpoiwbISvUIV12LrwDYrGiGg4xw7T6dRruuldZrhEsqWPI/Fx2q3z4DAMCZIeByM+9sTtMsm3tV0hxsjtW/r4+g+Y/2pWZhhj+cl58qoGryoHlbtftYE/XFVQqqQ4EDlqq37/URPBAzZNIv2sbDogCxQ6cBqooK2p+eTteLhNGXlhGRlBQWbrT8O9RdqBhw2WKESz19D02PwdE4WpXC91Zr2274+3pZ5DOiVYQ3dY7ztUm7EgAAW0LA5cbEfii29tadqTU+r1IHUrVRZ3+Rwse8wJGVXq+NB1XThOZeBlMI0woLaU9amqz0e0KIdh4CWEZwUJBNAq7le85q3nsIuMDR+KhTCst1e2rYSVyINrNApZ7LVR9sHu7ql4fRb88Oqfe2AAAcDUpxuZHScp3a43bSvWVDuvC/OL78xNdbaP2FErPmookXehVmjtQ9Naw9rTqeSaUV1eTj6UGncuSvz91NPWlAXCVVUbUs0MpWKulCXi4Vqos4hPj6UUp0NPnXUvodzBOiHuEqUSqpwtuHrHXt++XlZzTL4UGYQwKO2fi4wkECLj9JKfi8EusXtAEAcGYIuNzIias3Nct+no4xuvb1EwP4nIRmr66q0/erJMGjuQEX24d/pghXU1fuv0CTl53WPPdIp0h6qGM0nTx5kiqqqymntJTyS0ooo7iIlBUVmoIOyWHhlBgSghRCK/Lz9SWFQkEqlYqKKyutFnCVSqbGBKDxMTgYH3Ul1kobtEcwRZFKewwe37eVXfcFAMDRIaXQjZRK0j7WPNuLHAVLJWlqZF6XMWLPrPqkFIpBF7uN7NKMdrzUh4J8iHonKmjWmO6aoPBmRQUdzhGKYrBgiwVajUJCqWdCIjUKDUWwZQMh6lTNoqr6py4Z0ynO1+7ptgDOUjRDbBI+qXsMNYq2TU9HAABnZfcRLtbQdNGiRXTq1Cl69tlnqXnz5nrrpKen0+LFiykjI4NSU1Np3Lhx5KOTvmWpdVzZpKWHNH24mjSoWxU9a3uyTxN6+e/zmvtNXlnBA6CGkSG1jnD5mtF00xD2sw68NkzT0ys8PJyfeFdVVZHCy4tCFQqK9PenmIBA8sIJuU2FhoZSZmYmFVVaL+DyVf/dJ/dsYLWfAWAu8bjkKEUzitRVi0IDUOACAKA2dr2M++WXX1Lr1q1p165dtGDBAkpL069Od+bMGR4cbdq0iZeHnj17Ng0dOpQqJSdellrHlWXkFVGuun+v2IfLkfRtk0AKnUGuZXu0AZgucZI2KyBoydEIXx9vzfbCwsJoSN++1Cc4hHrENKCU6BiKCwpGsGXPES4r/r+q1CXh61L1EsDWRTMcZQ6X2Mw+FC0UAAAce4SrR48edPbsWcrNzaWlS5caXGfq1KmUkpJCK1as4KlbEyZMoKZNm9KSJUto/PjxFl3HlR28kEGOLDY8iPbPHEJXbxbQ8E/28MeUkrRBXcrScv7VQoNbRnl5eZEnSrw7TMBVXMWKpVRb/EoRa7h9MEOY+B9ow3YJAKZilVAdaYSLFRticIECAMDBR7jYiFNgYKDR58vLy2nVqlU89U/sa5SQkED9+/en5cuXW3QdVyed4OyoWO+VVglRfA4V89nuTGoxfQUt2nzc6AiXpFAWuLCAgADy9vLi9SJL1EVLLOmvfec0y8H+GOECx+MtjnA5SMDFiglJR94AAMA4hz5dvXLlCpWWllJSUpLs8eTkZD4yZsl1DGHfU1BQILs5q/c2GE/PczTS8wmWtTJz9SWjJe6tPcIFjoFdKBHLwxeWC6OblhQkabRa6SAntACG5nA5yvtTrN3h413/pscAAK7OoU9XlUol/8rmXOmmF5WUlFh0HUPmzJnDJ+uLt8TERHJWN0u0H9JPdI0mR6Y+rzCp4iLrnwXuQQy4iqwQcHnozCcEcNSy8BWOUaSQxKlkYqojAAAY59BHyqCgIP41Ly9P9jib8yXO6bDUOoZMmzaNV1EUb1evXiVnxOanqNPt6e9JXWj6XV3JkRmaM3XXR+tk90sr1EUzcHHVbYSo/4+LKiwfcInvpxh/D4oONZ7mDGD/gMtRRrjUKYVIMwAAcO6Aq1GjRjxYYs1npdj9Nm3aWHQdQ/z8/HhAJr05owPnb2gq+rVNjCJHZyjgOnCjjPadva65jxEu9x7hqrZw81fx/RTBeiYAOCBHa3wsjnD54qoXAIBzB1ysPPeYMWNo4cKFmrTAQ4cO0c6dO2ns2LEWXceVHbuWw7+yUS5naOh6dyfDqZtjvjnIe2IxpeoS3kgpdB9BgYE89Y9N1ldZuHBGsboipi+7KgHggAL9hOqZRaVVjjWHy5QccAAAN2fXI+WePXto8uTJNHPmTH5/3rx5/P7KlSs168ydO5dXGezUqROvMjhw4EB65JFHaNSoURZfxxWVlVfQ25uE/mZd4p2jQeXILs3or4ldKD5I/+T3dFqOLAUMJ8jug10sCFRfMCgsE0q4W0pesbC9ED+McIFjSogUUmrzHSXgUo9wIaUQAKB2dm04wwpRtGrVii+zQEgUFaVNe4uJiaGDBw/SmjVrKCMjg5555hnq3r27bDuWWscV/blXW+46RH2F1Bm0axJDq54fSL/vOUdvrL2sefzRhXtp+/ThmoALI1zuJcjLi4qqqnjAFVNDS4m6ylOqAy6F8/yPgHvxU6fuifNxHSbgwggXAECt7Hp2wYItMeCqCZtLNXr0aJus48qOZwrplM4iNFBBjw5MoQ83X6YC9YDGjeJqWrrtpHYOF3rAuF3AReXlVFRWatHt5iuFlMIQfwRc4JjEkSQxlc9xAi6MCgMA1AbJ1y5ozh97adI3W6iisoqk86vH3tKAnNH6Kf1lDY43ncqk97YIBTRwddW9BHp6WSWlMF8pVD4MD3COtFtwP2JxCjHQsSc2l1YsloiiGQAAtcPlXBdzI7eIvtiTxZeHv7uWzuWJqXdET996CzmjmLBAUtfI4DZd0vZO88UIl/uNcLE5fJWVVFZZSb4WuLp+MSOXtl1R8eVQBFzgoMSRJEcY4SqX7ATmcAEA1A4jXC5my/FrmmUx2GKahnmTN3Ltwcl5e3iQv/rEs8hCo1xLt5/RLDdr4JytH8D1iSNJLNQRq7XaS7mk+zJGuAAAaoeAy8VkFwlX6nUVS4eInFDvRIXBx7NLLFseHBxfkI8P/1pooXlchaXCeyg20IMGpja2yDYBLE06klQmCXjsoUzSlgFp3QAAtUPA5WJyig2fhDp75t2nj/Yy+HjrBgE23xewryBvMeCyzAhXqfrktVeTEKfoUwfuSTqSJA147J1S6OuDmQkAALXB2YWLyS0WJv/rClU4dyWpkAAFNQ/X/x16N3fOQiBgvmAfCwdc6tFfhbQyC4CD8ZEEXNKUPnsoU1eJ9aBqpKoDAJgAl6ZcTJ662pquW1OcPzAplZTnGpMSSrc0iqBRXZradZ/AfimFJeXlVFlVRV71HJUSR7gUPs59UQLcow8XU6buQ2jvES4vDydPnQAAsBEEXC4mT2U41eSJwank7MokHT/v7NyEerZKsOv+gH34enqSj6cnlVdVUVF5GYX6GZ7fV9eAS3pCC+BovLw8eUpKlUPM4ap0iVR1AABbQQ6NiylQGb7y6QpzU3o2CdYuI9hyWx4eHhTs68eXC0vrn1YoBlz+GOECBycWmrXUCFdpeQXv27jvrNDXsM4jXM7/sQIAYBM4XLqQ+SsP0tlc4YN4QBN/6hgrNHH1dZG/8qx7utCL/eJpx0t97L0rYGchfsJ7u6C01GKpqgpfjHCBYxNHlMorLRNwfbTiIO/bOOabg3X6vgp1wOWJES4AAJMgpdBFnLx6k97fqr1K+czQNpTaOIY+W3OYBrRtSK6AFc6Y7KTNm8GyhDTCfMorNdwGoS6U6qIZgb44HIKTBFwWSik8mlZg1vchpRAAoG5whuEisgpKZPdjQgN49ainRyBAAdcT6iekFCorKqisspJ81c2QzSGOCjeJQdNjcGw867WCqEhlmQqdVdppsXWClEIAgLpBwOUidp/LlN2PCw+y274AWJuPlxcF+vhQcXk55ZeqKDog0KztHLus/b9p1TDCgnsIYHmR/l6UV1pJm05cpxNpufSfASm8mEZt/txzhgL9fKhpbBj9ue8Cf6xAVU47r5k3QiyOsHmjSiEAgEkQcLmIT3dlaJZnDmnkEkUyAGob5eIBl6rU7IDrYqY2pSoiGE20wbEJ/RQr6et9N/l9Tw8PmjBIXoG2SFlG3246Ti1iQ2l4x2TaevwKPffHWf5cmB9RXv2nPZKyTKiG64MyhQAAJsFZuYthvVt1P4ABXJFYDp6NcJlr22nhQkVqtNDbC8CR6Taw/+fIDb11vtl4jD7Yep0m/nKSVx8cv+io5jlLBFvSgAt1ZgAATIOAy8W8PqyJvXcBwCZCFcI8roKyMqqqrvtklJLScvrlWJ56W6hQCI4vRCG/MHAwQ38uV0aB9gJEXasPmiqnWIjc/FAXHgDAJAi4XIS3unpVj5bx9t4VAJsI8PYhb09PHmwVltW9iMCby/Zqlh/pnWzhvQOwPH8Tenwoy80rGV+pLvVuioOXc4UFlIUHADAJAi4XwK7UV6gv8IcFCFf9AdyhAXKYulphvqruaYVrTwujW8zA1MYW3TcAaxjYpvYLaiVl5gVc+SU1/w9dzy6k5XvPUlVVFV3MEdZViFf6AACgRgi4XEBBiTYxPzRQmNcC4E7zuMzpx9UpQSi0ERPggSIz4BQMVdIsLRfmUx25lEl7zqRp+srV5ttxbSkxWHsKMPPXfUbXvXAjl3q+u5We/f0MfbL6MJ3MFn5mh4ZopQAAYAoEXC4gt0g42fTzJN57C8BdhPv786+5KhVV13EeV1q+cKHi7nbRVtk3AEsLC9TPYLhZUEIVlVU0+vN9NPbbQ3Qys/aLD1MHNKSB7ZrQtldv1Ty24kyR0fUf/ma3ZvmHvemaZVQpBAAwDc7OnRy7uvn1plN8WYEi/+Bmgn19ycvDgyqqquo8j+uM+ip9gB/+ccA5BCl89R67kVtMV7LyNfezlLVfeJg0rEONo1lbjl2RPXa1UDtqdlOl3T4u8AEAmAYBl5Ob/uNuWnZC+LA1M3UfwGmxPkThCnGUS1mn7w1UF3xr3iDUGrsGYHGG+itevlkgSysXDWhSt75yrFAnm5818KOd9PDio7Tj5FX++KUM7VxHXb7eOIUAADAFjpZOTgy2GCUCLnBD4f7CPK4cZe2pVIcu3KCcQiWVlVdQvnpArGFkkLV3EcBqZvxzjlQGKhN+N3EA7X2lP/Vt7E/z725J4epsxFCdQbKkUOE0IDXGl7aeuKZ5fNfZDP5/0v/DHUZ/duMo/O8AAJgCuTQupG0U/pzgfiLUI1yscAYrEc9GvQxhV+wf+P4IxQd5UK8m2sn+Ch/834DzKqlgqeXygMtP3VYuJiyQfpg0kC8nRAbR3BXH6OVb28rWHd8tkd5Ye5lO3SyjRxZrmyTvuphLlxfvrPFn396lueV+EQAAF4YzDSfXKMSTrhQI+fWvjU6x9+4A2Fygjw/5eHpSeVUVLw8vFtLQtXTXBf71elE1/XpMOzIc6CdvJgvgTCL8iEor5AFXqYFsh1uSY+nnp2P1Hk+IEKp1FpbLH9+fXsZvNfFCkSYAAJMgpdCJsXz7bKUQbL03uhl1a9HQ3rsEYJd+XJH+wnyVm8oSo+tVVBkuJhAWhFYK4LyC/Dxpyylt5cC6ah4Xbtb33dkGJeEBAEyFgMuJfbjiIBWrr0qydBEAdxUVIARcWSUlRsvDG2vRGoARLnAiM4c0ovYxPvwiG3OtoIqWHMo1e3vxEaZ9drw2pJHs/rgeTc3+mQAA7gYBlxObv+OGZrlJDCqtgfuK9PfnAZWyooJKynVyo9QMhWGtI5FVDc5lwqBUWj5lKMWrUwFNa3NsnK+BOYz3pOh/ntzRVQjwpCmKAABgGgRcLiI6pG4lgAFcibenp6Y8vLG0QkMDX789PcDauwZgFQG+xi8W6FYirKsgA73pAhXauY6RCvTgAgCoCwRcLgKTl8HdRavTCjNLDAdcqgr5WMD/RiRRoIFGsgDOQBoAMVH+HvTZva14mffPxrWv07a6xqtrxqsF62yb8ZOMhOn8KwEAQC0QcDmp7eqmlEwAsqIANAFXQWmpwbTCEp3O4N6osAZOLNBPfrHgh/90oVs7NqVN026lnq0S6rSt757oK7sf6l9z5c5AH2MzIgEAwBAEXE7quZ+1/VL+2yvervsC4Aj8vL35XC7mRlGR3vNFZfLL8hFB8qv6AM4kQDIKNXVAQ2qTGG32tthIb5AkxgoPlP9vdI4Tgrt3RjXlfew+uLduI2gAAO4OAZeTkl6cj9D5cARwV7GBQsW19OIivWqFp7IrZPf7tkm06b4BWGsOlyVGa+ffm0INAjzogzuaUYROq4QfnuzHv97bqxXtnDGCurdECxIAgLpAMpqTKqvUnkyG40o9gCat0NvDk1QVFZStVFKElxd/vLRcHmx9ck9Lg9XZAJyF9P0bHlD/i24DUhvTntTGfHmHJGWdjXwFYK4jAEC94IzDCZWVV1BuqfZ+VLCQRgXg7rw8PalhcDBdLsiny/n5FBERwR+/kVfMv7JxgLOzb0WRGXAJfRop6Fp+GY3olGzR7XZuFkfJoccoNtiXlk4eZNFtAwC4IwRcTujo5SzZ/QgEXAAaCSEhdKUgn/JKVZRbWkqso9Cu85n8uWBfVPQE17Hov9YJhlhFwo3TbrXKtgEA3BHmcDmhrAKl7H4kAi4ADYW3Nx/lYs4U5FNVdTXdLBKGhPPL8EIBAACAbSHgckLDOyZTl3ihahTruRLij15CAFLJYeHk4+lJJRUVdEalpKJSYQ5Xv8ZIvwUAAADbQkqhk/r1mSFUUVlFnh5Enp6ImwGkfLy8qG10NB3KyKCM8nK6lplNRB4U4CsU0QAAAACwFacIuPbt20fr1q2jvLw86t69O915553k4aFtvPjzzz/Ttm3bZN8THx9P06dPlz2Wnp5OixcvpoyMDEpNTaVx48aRj0/NDR4dGRq3AhgX6R9AKWHhdCovl0rVPY/DgjDCBQAAALbl8EMjb7/9NvXv35+ysrIoNDSUpk2bRvfdd59snU2bNtGOHTuoVatWmltSUpJsnTNnzvAgi60bHBxMs2fPpqFDh1JlpfpMDABcTrS/P3UPDqE3bu9Am5/rQa+MusXeuwQAAABuxqFHuNiI1quvvkqff/45PfbYY/yxxx9/nJo0aUIrV66kESNGaNZt2rQpTZ482ei2pk6dSikpKbRixQo+OjZhwgT+PUuWLKHx48fb5PcBANvz8vCg+OhI8vTH6BYAAADYnkOPcF26dImPQHXs2FHzWExMDDVs2JB+//132brnzp2jl19+mY9cbd26VfZceXk5rVq1iqcQiqmICQkJfORs+fLlNvptAAAAAADA3Th0wNWsWTNSKBS0du1azWNnz56ly5cv8xRBEQuiGjRoQGFhYXT9+nUaPnw4TZo0SfP8lStXqLS0VC/NMDk5mW/PGPY9BQUFshsAAAAAAIBLpBQGBQXRggUL6KmnnuKjVlFRUbR9+3Y+F0up1PaiYsUxEhMTNffvvfdePnp111130ZAhQzTrsrlbUiEhIVRSUmL058+ZM4feeOMNq/xuAAAAAADg+hx6hIv5z3/+w0ezHn74YRo4cCAPuOLi4njwJZIGW0y/fv34Y2xdMXAT54RJ5ebm8qDLGFagIz8/X3O7evWqhX87AAAAAABwZQ49wiViwdPYsWP5cllZGe3cuZOmTJlS4/ew9SoqhGanjRo14kHXyZMnebqhiN1v06aN0W34+fnxGwAAAAAAgEuOcO3du1cTODFvvfUW/ypWLWTPbdiwQfY9P/zwA++1NWzYMH6fNQYeM2YMLVy4UJNeeOjQIR64iYEcAAAAAACA241wsaIWLK2wa9eudPr0aV6N8I8//qDY2FhNwYz333+fp/+1bduWF8jYtWsXn3/Vt29fzXbmzp3L53V16tSJOnToQKtXr6ZHHnmERo0aZcffDgAAAAAAXJnDB1wPPPAA9ezZkxfNuP3222nQoEGaOVmMl5cX78l15MgROnjwIIWHh9PixYv5PC8pVk6ePb9mzRo++vXMM89Q9+7d7fAbAQAAAACAu3D4gIth5dx1S7rrateuHb/VhM3HGj16tNn7UV1dzb+iPLz7qFIqqbi4mDzKysjDx6fe26suL+e3yoIC8iwvt8g+uuprrftamfO3wOsNAOYeg0zlKMeZ+n5eOcrv4W4sfZ5RF/ib158YE4gxglMHXI6isLDQYFVEAAAAAABw3xghNDTU6PMe1bWFZKBRVVXFGyuzfl5s7hg45pUGFhCzEv41lfwHwHsGcJwBW8FnE+A945pYGMWCrfj4eF6kzxiMcNUBeyETEhIs8fcBK2PBFgIuwHsGcJwBR4LPJsB7xvXUNLLlNGXhAQAAAAAAnBUCLgAAAAAAACtBwAUuhVWifP311/lXALxnAMcZcAT4bAK8Z9wbimYAAAAAAABYCUa4AAAAAAAArAQBFwAAAAAAgJUg4AIAAAAAALAS9OECp1BZWUmnT58mLy8vSkpKIl9fX6PrpqWl0cWLF6lly5YUHR0te66srIwuXLhAKpWKkpOT0avLhVVUVPD3jI+PD3/PsK9SBw8epOLiYtljcXFx1LRpU71tnTx5kr9nUlJS9LYDrqO8vJy/Z1iBA/ae8faWf0T++++//H0gxXozNmnSRPZYSUkJnTlzhh9/GjZsaJN9B/tgnynsPePv78/fM+wzSmrv3r18HalGjRrxmyHs8+n69euUmppqUm8fcD7s/XDq1CkKDAzkxw7d98zu3bv555dU48aNKTExUW9b+fn5dP78eX4+ExYWZvV9h3qoBnBwH3zwQXVcXFx1ampqdZMmTaqjo6OrlyxZYnDdgoKC6hYtWlSzt/aiRYtkz/3+++/VDRo0qE5OTq5u3759tb+/f/XUqVNt9FuALc2dO7c6Nja2ul27dtWNGjXiy7/88otsnbZt21YnJSVV9+rVS3N79913ZetcunSJbyMqKoq/99j7Z9OmTfhjuqBZs2ZVx8TE8GNDQkJCdcOGDav//vtv2TqNGzeubtasmew988knn2iev3HjRvWjjz5aHRoaWt2hQ4fqiIiI6h49elSfP3/eDr8RWFNlZWX1jBkz+OcR+1vHx8fzY83q1atl67FjBvtMkr5nvvzyS4PbTEtL4+uzzy8cZ1xPWVlZ9SuvvMI/T9h7hp3XsM+VDRs2yNZjx49WrVrJ3jPff/+9bJ3y8vLqp59+mp/H3HLLLfy99+abb9r4N4K6QMAFDu+dd96pzsvLk9338vKqzsrK0lt33LhxPIjSDbjYgS4wMFAWYK1fvx4fbC4ccLHgW8Q+iHx8fKrz8/NlAdeHH35Y43Z69+5dPWjQIP7+YV544YXqyMhI2XbA+VVUVFS//fbb1UVFRfx+VVVV9csvv1wdEBBQrVKpZAHXV199ZXQ7+/btq/7uu+/4yRDDtjdw4EAedIFrKSkpqX7//ff5V/E98+yzz1aHhITw95OIBVC6F/+MBXADBgzQfH4h4HI97HNj3rx51UqlUvM3nzhxIg/A2PtHGnD9+uuvNW5r8uTJPMg/c+YMv8/ec5999pmVfwOoD8zhAof30ksvyVIr+vfvz1MMMzIyZOt9++23dPbsWZo1a5beNliKD7v16NFD81jPnj3Jw8ODsrKyrPwbgK1NnTqVgoODZe8Zli6m+7fOzc2lffv2UXp6ut422Htp+/btNG3aNE0aIVtmKRx//fWXDX4LsBWW0vPyyy/zFB+GHRfYe4YdM/Ly8mTrZmdn89RC3eMP07lzZ3rkkUc0qYhse+PHj6c9e/bopQiBc2MphFOmTOFfxfdMv379qKCggAoLC2XrsuMOe89kZmYa3d5bb73F3y8TJ060+r6DfYSEhNAzzzxDCoWC3/f09KS+fftSTk6OXqoye6+w94yh8xOWcvr555/z90zz5s01xzC8dxwb5nCBU7h69SpdvnyZbty4Qe+88w498MAD1LZtW83zLB+anQyzE2TdeRcMC9heeeUVmjFjBj/xDgoKoq+++op69+5No0aNsvFvA7Zw5coVfmMfTnPnzqX//Oc/evOz3n//fVq+fDnPgW/Tpg19//331KpVK80cL6ZTp06a9SMjI3muPHvuwQcfxB/SxbC5n9euXePzQNnJDDs5atCggWyd2bNn048//kjnzp2jjh070sKFC/l7whgW0LP5F4aOS+D82LGDHWPYZ9T//d//8cBddy4NuwjI5nexizjdunWj7777jr8nROxz68svv6QDBw7wIB9cG3sfsHOZS5cu0ZtvvkmvvvqqJnAXscfYPD82F5QFZeyCsjgfdMuWLfwCzsiRI/lnHAvYmjVrxs9rwIHVa3wMwEZY/nLPnj2rExMT+RyKjRs3ap5jKT9s3oWY6sPSeQzN4Tpy5Eh1ly5deFpQSkoKn9dT27A9OC82T4K9Z9h8nJYtW1Zv3bpV9vzChQs16WIsZXXIkCHVrVu3ri4tLeWPsfQMlroqTfVgWHrYI488YsPfBGyFpfuwv684Z3TPnj2y57/99ltNeunNmzer+/TpU92xY0dZCpkUSwvz9vbm3weuic37ZO8Z9nnC5tLs379f9vw333yjSTHNzMys7tatW3X37t01x5Xs7Gw+/2blypX8/sWLF5FS6OJmz57N3wNszig7Jzl8+LDseXYuIx5T0tPT+TGmb9++muffe+89PkXiiSee4J9v7HyGzeXCHC7HhoALnM78+fP5SczJkyf5/ZkzZ/IDzrZt2/ht8+bN/AOLPX7w4EHNyVFYWFj166+/rtkOOwFn2xE/6MA1sRMbdlLk6+tbY/EC9l5h75u9e/fy+19//TW/LwZgIjbZmX3QgWu/Z1gRDXZSwwoZGMOON+w9cuLECb3n/v33Xz4Xg837A9fH5uNMmzaNz+FigZUx69at4++ZCxcu8PusyAqb5yd+frGLgOz5jz/+uPr48eM2/A3A1lhQNWXKlOrw8PDqnJwco+ux4j3sPXH9+nXNhSF2//HHH+fvO4YVa/Hw8NAr2gKOA3O4wOk89dRTvCz8+vXr+X2WDy2mDLIbG4pnfv75Z5o3b54mZYPNxWDfK+rTpw+1a9eO/v77bzv9JmALbG7Fc889R1VVVbRx40aj64mpYyydjBFTfli6kBS7b6ykM7jOe+bFF1/kbQO2bt1q8ntGxFLDhgwZQo8++ii99957Vt9fsD82H+eFF17gc7h27Nhh8nuGfXaVlpZqPr/efvtt/jhLMWQpzuC62LwrNg+QzSVm8zxNfc+IbSgef/xx/r5jhg0bxtNWt23bZpN9h7pDUjk4NJbPznKb2QmQiB102ONsPg0zffp0fhOx3GZW5GDmzJmaeTZiPy42P0NcZuuxPGrdXl3g3Nh7IyAgQPYYm//H/t7ie8bQ+2rt2rX8qzg3kBVYYZPYWYEMNpeH2bVrF5/MzE6mwXWwwEosmCFic7QY8T1jaB32nmEnPK1bt9Y8dujQIf7+eOihh+jDDz+0yf6Dc71n2Hw+ca6o7nuEzethJ87z58/nhVvAddTnPcMuMosFMticLnahmZ0LdenShT/Gim6wgj44n3FcCLjAoR0/fpyef/55evjhh/mHEJsgygodtG/fnu68806Tt9O1a1c+WZkV23jttdd4taCvv/6aV5Ni2wbXwUYX2CgnO+FlVwLZCcy7777LP5jYJGPm8OHDvPolqyjHRqtYNag5c+bQpEmTNB9q7EOPvVdYMM8+7MLDw/l22fuOvZ/AdbARcDYSNW7cON5clJ0EseI87MRm4MCBmmCbFcxg7ys2eZ01J2WjEWxUQ5zMzia4Dx48mDfIHjNmDN+uiL1namrYDs6FZViwSnFjx47lza/Z3569H1iwzYoxMZs3b6aPPvqIv6/i4+P5yBc7FrECT1FRUfb+FcDG/vnnH1q0aBE/NrBjxsmTJ/l7hn0uiYHT6tWr+bnJfffdx0e22Ag7O+d5/fXXNdWaWVEWdkGZXQhkI6rsvfTpp5/yohko5uS4PFheob13AqAmR44c4ekV7AONHVhY6V12ouzn52dwfVYynq3DTpaHDh2qebyoqIgWLFjAT5TY1SB2hfHpp5+uscIYOCdWRZBVoWQnzuyK34ABA3h5bukJr3QddsJ01113aQIyqSVLlvD0VJb2w06+WXqisfceOC+W0sOqx7Gqc7GxsTxwYhdopNUF9+7dS9988w0P4llgxk62paOd7Eo0qzpmyB9//IGrzy6GBVAs7Y9Vt4yLi+NpXexEmaWKiXbu3MnfV+xiIbu4w4IvdjwyhmVd3HPPPXyE65ZbbrHRbwK2wioMsqCLZV2wIHzEiBE8ABNTA6XrsMqX7KIhC6LYFAhdrFqq+NnELkKz9MSYmBj8MR0UAi4AAAAAAAArQdEMAAAAAAAAK0HABQAAAAAAYCUIuAAAAAAAAKwEARcAAAAAAICVIOACAAAAAACwEgRcAAAAAAAAVoKACwAAAAAAwEoQcAEAgEsrKyujn376ifLy8uy2D1lZWbRixYpa18vPz+dNkgEAwHUg4AIAAJehUql4cMUCF1FBQQHdf//9dOnSJbvt15QpU2jfvn21rhcSEkKzZs2i33//3Sb7BQAA1oeACwAAXAYbxWLB1dWrVzWP+fn50dixYyk8PNwu+3T06FEeQD3//PO1ruvh4UFTp06l6dOnU3V1tU32DwAArAsBFwAAuAQWoCxfvpwvr169mo90bdiwgQdcd9xxB4WGhvLniouL+XNFRUV07tw5+vvvv+nw4cOa7Rw5coRvhz1nCBsxW7NmDb9lZGTUul8LFiyQ/XymvLyctm3bxtMMr127Jlv/zjvvpLS0NL7vAADg/LztvQMAAACWCrhWrVrFl1mwwgKcVq1aUfv27fmo18GDB6lDhw58PhW7369fP8rJyaGEhARav349Pf3005Senk4nTpyguLg42rRpE3355Zc0fvx4zc9YtmwZPf7445SamkoBAQG0a9cumj17Nj311FNG94sFVSxNUMQCLPaz2fcnJyfTsWPH6OGHH6bXXnuNP+/v70/du3enf/75hwYPHow3BwCAk0PABQAALsHT05M+//xzPjr17rvvUkpKCn/85s2bBtdv3rw5D6hYGh/7vkmTJtHkyZPp0KFD/Pn33nuPZsyYoQm4Ll68yAMjFkCxgIn5999/qU+fPjRo0CAe3OnKzMzkAVabNm00j33zzTcUGxtLO3bs4Perqqr4KJsUC+j27NljsdcGAADsBymFAADglp544gkebDE9evTQPCZij7G5YEqlkt//8ccfKTIyko+Q/frrr/TLL7/QhQsXKCIigqcHGiIGe9L5Y2wEKzs7m6cNioHi7bffLvs+tr6xQBEAAJwLRrgAAMAtSYMgNs/L2GOlpaU8SGJVDlmJ+d9++022HTbCFRUVZfBnBAUFaeaNidhIGhtFYyNsbORryJAhfGStYcOGmnXY+sHBwRb7XQEAwH4QcAEAAJiAlWxnARkruGGq+Ph4HnSxdMROnTrxx1ggtXTpUh5UsbTC+fPnU+fOnXmRjsDAQL4OW79ly5b4uwAAuACkFAIAgMsQR5RYPy5LGz58OJ08eZK2bNkie7ywsJDfDPH29qa+fftq5msxYiohC66GDh1KH3/8Md24cYMuX76sWWfnzp0omAEA4CIwwgUAAC4VcLGRobfffpuXV2/QoAGvUmgJrGIgm+M1cuRIXtGwadOmdPr0afrjjz9o7dq1RlMAWVVDljLIinB4eXnRZ599xqsbjhgxgsLCwmjJkiXUrl07atGiBV+fBWes9PyYMWMsst8AAGBfGOECAACXwsqps3LrK1eu5MUsdBsfs5Eldl8aILF0QfYYK9UuYsUw2GO+vr6ax7744gs+h4v18GKjUCygY8FTUlKS0f0ZPXo0RUdHa+Z+vfXWWzRz5kxegn737t08MGT7yUbDGJZi+Nxzz2lG6wAAwLl5VKOVPQAAgFWxwIoFgG+++WaN6+Xm5tLzzz/PR8FYoQ4AAHB+CLgAAAAAAACsBCmFAAAAAAAAVoKACwAAAAAAwEoQcAEAAAAAAFgJAi4AAAAAAAArQcAFAAAAAABgJQi4AAAAAAAArAQBFwAAAAAAgJUg4AIAAAAAALASBFwAAAAAAABWgoALAAAAAADAShBwAQAAAAAAkHX8Pxze/TcL0PVYAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(10, 3.5))\n", "sel = (d.tx > 347_000) & (d.tx < 357_000)\n", "t = d.tx[sel] / 1000\n", "ax.plot(t, d.get_values(\"left_pupil\", policy=\"raw\")[sel], color=\"0.7\", label=\"interpolated\")\n", "ax.plot(t, d.get_values(\"left_pupil\", policy=\"strict\")[sel], color=\"C0\", label=\"recorded\")\n", "for on, off in d.get_blinks(eyes=\"left\", variables=\"pupil\", units=\"sec\")[\"left_pupil\"]:\n", " if t[0] < on < t[-1]:\n", " ax.axvspan(on, off, color=\"C3\", alpha=0.2)\n", "ax.set_xlabel(\"time (s)\"); ax.set_ylabel(\"pupil size\"); ax.legend();" ] }, { "cell_type": "markdown", "id": "58534cbf", "metadata": {}, "source": [ "The red stretches are the detected blinks, the grey line the values that\n", "the default policy uses there.\n", "\n", "## Pupil samples during off-screen gaze\n", "\n", "The off-screen and divergence layers mask only the gaze channels. Whether\n", "the pupil size is still valid when the gaze is off the screen depends on\n", "the setup; to mask the pupil there as well, copy the off-screen periods to\n", "the pupil with `mask_intervals()`. They go into the pupil's `\"user\"`\n", "layer:" ] }, { "cell_type": "code", "execution_count": 18, "id": "7ac53bb9", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.582802Z", "iopub.status.busy": "2026-11-09T20:27:41.582444Z", "iopub.status.idle": "2026-11-09T20:27:41.720109Z", "shell.execute_reply": "2026-11-09T20:27:41.718478Z" } }, "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", "
NaNblinkoffscreendivergenceexcludedinterpolatedusermasked (default)masked (strict)
left_pupil0.05.840.00.00.357.490.971.048.53
right_pupil0.04.080.00.00.355.661.041.046.70
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence excluded interpolated user \\\n", "left_pupil 0.0 5.84 0.0 0.0 0.35 7.49 0.97 \n", "right_pupil 0.0 4.08 0.0 0.0 0.35 5.66 1.04 \n", "\n", " masked (default) masked (strict) \n", "left_pupil 1.04 8.53 \n", "right_pupil 1.04 6.70 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d3 = d.copy()\n", "off = d3.get_layer_intervals(\"offscreen\", variables=\"x\")\n", "d3.mask_intervals({\"left_pupil\": off[\"left_x\"], \"right_pupil\": off[\"right_x\"]})\n", "layer_table(d3).loc[[\"left_pupil\", \"right_pupil\"]]" ] }, { "cell_type": "markdown", "id": "1f7f3e93", "metadata": {}, "source": [ "The pupil is now in the `\"user\"` layer for about 1% of the samples, but\n", "the default policy masks fewer of them. The reason: under the default\n", "policy, an **interpolated sample is always usable**, whatever other layer\n", "of the channel is set and in whichever order the layers were set\n", "(only excluded periods mask it). Off-screen periods that fall into an\n", "interpolated blink therefore stay usable:" ] }, { "cell_type": "code", "execution_count": 19, "id": "75bd78e7", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.723391Z", "iopub.status.busy": "2026-11-09T20:27:41.723036Z", "iopub.status.idle": "2026-11-09T20:27:41.731587Z", "shell.execute_reply": "2026-11-09T20:27:41.730139Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "user layer: 0.97%, of that interpolated: 0.28%\n" ] } ], "source": [ "user = d3.data.get_layer(\"left_pupil\", \"user\")\n", "interp = d3.data.get_layer(\"left_pupil\", \"interpolated\")\n", "print(f\"user layer: {user.mean():.2%}, of that interpolated: {(user & interp).mean():.2%}\")" ] }, { "cell_type": "markdown", "id": "c014c50c", "metadata": {}, "source": [ "To leave these samples out as well, use the strict policy (which leaves\n", "out all interpolated samples), or a policy that lists the layers to mask\n", "without the exception for interpolated samples, e.g.\n", "`policy=[\"user\", \"excluded\"]`.\n", "\n", "## Filtering\n", "\n", "The filters read the usable samples only. Unusable samples are bridged\n", "linearly for the filter and keep their values, so neither blinks nor\n", "excluded periods leak into their neighbours:" ] }, { "cell_type": "code", "execution_count": 20, "id": "4ffbf083", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.734568Z", "iopub.status.busy": "2026-11-09T20:27:41.734225Z", "iopub.status.idle": "2026-11-09T20:27:41.940560Z", "shell.execute_reply": "2026-11-09T20:27:41.939387Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "masked samples unchanged: True\n", "usable samples filtered: True\n" ] } ], "source": [ "before = d.data[\"left_pupil\"].copy()\n", "d.pupil.lowpass_filter(cutoff=4)\n", "masked = d.mask(\"left_pupil\")\n", "print(\"masked samples unchanged:\", np.array_equal(before[masked], d.data[\"left_pupil\"][masked]))\n", "print(\"usable samples filtered: \", not np.allclose(before[~masked], d.data[\"left_pupil\"][~masked]))" ] }, { "cell_type": "markdown", "id": "aa4955e7", "metadata": {}, "source": [ "## Event-related analyses\n", "\n", "Segments around events (`get_segments()`) and statistics per event\n", "(`stat_per_event()`) use the policy as well. With the default policy, the\n", "interpolated samples are part of each trial; with the strict policy, they\n", "are masked. This lets us find the trials that consist mostly of\n", "reconstructed data:" ] }, { "cell_type": "code", "execution_count": 21, "id": "05aa248c", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:41.943344Z", "iopub.status.busy": "2026-11-09T20:27:41.943012Z", "iopub.status.idle": "2026-11-09T20:27:42.003842Z", "shell.execute_reply": "2026-11-09T20:27:42.002169Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "

SegmentedEyeData — S_

\n", "\n", "
\n", "
\n", "
Epochs
\n", "
S_
\n", "
\n", "
\n", "
Segments
\n", "
300
\n", "
\n", "
\n", "
Time points
\n", "
1,500
\n", "
\n", "
\n", "
Window
\n", "
-500 to 2498 ms
\n", "
\n", "
\n", "
Missing
\n", "
0.0%
\n", "
\n", "
\n", "
Sampling rate
\n", "
500 Hz
\n", "
\n", "
\n", "\n", "\n", "
\n", " Channels (mean over segments; red line = event onset)\n", "
\n", "
\n", "
left_pupil (0.0% missing)
\n", " \"left_pupil\n", "
\n", "
\n", "
\n", "\n", "\n", "
\n", " Epochs\n", " 
Epochs 'S_' | 300 epochs, window (-500, 2500) ms
\n", "
\n", "
" ], "text/plain": [ "SegmentedEyeData(S_): 300 segments x 1500 time points (-500 to 2498 ms, 500 Hz)\n", " channels: left_pupil (0.0% missing)\n", " epochs: Epochs 'S_' | 300 epochs, window (-500, 2500) ms" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stim = d.get_epochs(\"S_\", window=(-500, 2500) * u.ms)\n", "seg = d.get_segments(stim, \"left_pupil\")\n", "seg_strict = d.get_segments(stim, \"left_pupil\", policy=\"strict\")\n", "seg" ] }, { "cell_type": "code", "execution_count": 22, "id": "1f1d3607", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:42.007004Z", "iopub.status.busy": "2026-11-09T20:27:42.006607Z", "iopub.status.idle": "2026-11-09T20:27:42.158619Z", "shell.execute_reply": "2026-11-09T20:27:42.157248Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAEmCAYAAAAk6TltAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJsRJREFUeJzt3Qu0jPX++PHPdtuuewsR2WwicjnIJbcfkrsIJylOKCdlZSc6QqVQQrdVIc6qSCWX1JFrLB23UOQeTu73+/2+3Z7/+nzXmvnPbLO32WN/zZ79vF9rzZqZZ5555vt855lnPs/3GuU4jiMAAAAWZLKxUQAAAAINAABgFSUaAADAGgINAABgDYEGAACwhkADAABYQ6ABAACsySIucOPGDTl48KDkyZNHoqKiwp0cAAAihg63de7cOSlSpIhkypT68glXBBoaZMTFxYU7GQAARKx9+/ZJ0aJFU/0+VwQaWpLhyaSYmJhwJwcAgIhx9uxZc7Hu+S9NLVcEGp7qEg0yCDQAAEi9UJse0BgUAABYQ6ABAACsIdAAAADWEGgAAABrCDQAAIA1BBoAAMAaAg0AAGANgQYAALDGFQN22RLff3ZI79s9vGWapwUAgPSIEg0AAGANgQYAALCGQAMAAFhDoAEAAKwh0AAAANYQaAAAAGsINAAAgDUEGgAAwBoCDQAAYA2BBgAAsIZAAwAAWEOgAQAArCHQAAAA1hBoAAAAawg0AACANQQaAADAGgINAABgDYEGAACwhkADAABYQ6ABAACsIdAAAADWZJEwW7RokcyZM0dOnTolxYoVk86dO0vx4sX91tmxY4eMGzdOjhw5IhUrVpTu3btLjhw5wpZmAAAQASUa7733nrRo0UKyZs0qNWrUkPXr10u5cuVk3bp13nU2bNggVapUkZ07d5ogQwOO+vXry5UrV8KZdAAAEIQox3EcCZPSpUtL27ZtTcChNCllypSRNm3aeJdpIHL9+nWZN2+eeX706FFT4vHpp5/Kc889F9TnnD17VmJjY+XMmTMSExOTZumP7z87pPftHt4yzdIAAIBNt/sfGtYSjfj4eNm/f7/3+fnz5+X06dNSsmRJ81xLLRYsWCBPPPGEd52CBQtKw4YNZdasWWFJMwAAiJBAY8KECSa4qF69ujz++ONSuXJlSUhI8JZU7N27V65evXpTmw19rlUpyUlMTDQRmO8NAAC4LNBYvHixrFixwpRQNGnSxLTFmDhxouzZs8e8funSJXOfO3duv/flyZPH+1ogw4YNM8U8nltcXJzlPQEAAOkq0Lh8+bLpPTJgwAAZMWKEeTxt2jRTNdK/f3+zjgYJSnuk+Dp58qT3tUB0m1qX5Lnt27fP8t4AAIB0FWgcP37cVJuUL1/eb7k+37Vrl3msJRF58+aVP//802+djRs3mh4oyYmOjjYNVnxvAADARYHGvffea0ovfvzxR++yc+fOyfz5800VioqKipInn3xSvvzyS/OaWr58uaxcuVI6duwYrqQDAID0HmhoEKGNQX/44QfTCLRdu3Zy//33mxKMoUOHetd79913TZuMChUqmK6uTZs2ld69e5s2HQAAIH0L6zgaSksqVq9eLSdOnDC9SapWrWqCEF86jsayZcu8I4OWLVs2VZ/BOBoAAITmdv9Dwz4EuZZWNGjQIMV1MmfOLPXq1btjaQIAAGmDSdUAAIA1BBoAAMAaAg0AAGANgQYAALCGQAMAAFhDoAEAAKwh0AAAANYQaAAAAGsINAAAgDUEGgAAwBoCDQAAYA2BBgAAsIZAAwAAWEOgAQAArCHQAAAA1hBoAAAAawg0AACANQQaAADAGgINAABgDYEGAACwhkADAABYQ6ABAACsIdAAAADWEGgAAABrCDQAAIA1BBoAAMAaAg0AAGANgQYAALCGQAMAAFhDoAEAAKwh0AAAANYQaAAAAGsINAAAgDUEGgAAwBoCDQAAYA2BBgAAyNiBRmJiovz222+yYcMGcRwn4DpbtmyRRYsWyZEjR+54+gAAQIQGGt9++60ULlxYnn/+eXOrV6+eXzBx4cIFadq0qdSuXVv69u0r8fHxMnz48LCmGQAABCeLhNGCBQukS5cu8t1330mHDh3MspUrV8qpU6ekUKFC5vmbb74pW7dulW3btkmBAgVk3rx50qxZM6lbt665AQCA9CusJRpDhgyR5s2be4MMVaNGDSlbtqz3+ddffy3dunUzQYbS0o0qVarIhAkTwpJmAAAQAYHG5cuXZcWKFdKiRQs5duyYLF682JRa+Nq/f78cP37cBBa+9Pn69etTbPNx9uxZvxsAAHBRoHHixAm5du2aaQRauXJlGThwoNSqVUvq1KnjbaNx+vRpc58vXz6/9+bPn99UryRn2LBhEhsb673FxcVZ3hsAAJCuAo1s2bKZew00Nm3aJEuWLJFdu3aZ4OKVV17xW+fSpUt+77148aL3tUAGDBggZ86c8d727dtndV8AAEA6awyqbS7y5MkjrVu3lrx585pl+rxt27amcajSkojMmTPfFCholYr2PklOdHS0uQEAAJeWaERFRZneI1qK4Uufa3dXlSNHDtPddfr06d7XtYRCe6voewEAQPqWJdy9TrRdhlaVaNuM33//XaZOnSozZszwrvPuu+9KgwYNJCEhwaw7ZswYU9KhPVEAAED6FtburdqNddWqVaZRqHZj1bYY+ly7vHrUrFnT9E7RniRTpkyR+vXry7JlyyRnzpzhTDoAAAhClJPcmN8ZiHZv1d4nWu0SExOTZtuN7z87pPftHt4yzdIAAEB6/g8N+xDkAAAg4yLQAAAA1hBoAAAAawg0AACANQQaAAAg/Qca586dExd0YAEAALYDjS1btki/fv28z3v16mW6vBQpUkTWrVsXyiYBAEAGFFKg0adPH2ncuLF5vHXrVvnyyy9l1qxZ8tRTT8mrr76a1mkEAABuGoJ8+fLlMm3aNPN43rx58thjj0nLli3NMOIlSpRI6zQCAAA3lWhkz55dDh8+bB7PnTtXGjZsaB5fuXIlxenbAQCAu4RUoqGlF+3atZOqVavK0qVL5auvvjLLf/nlF2+VCiIfQ6wDAMJSojFq1Chp3bq1mQxt5syZUrBgQbNcg4633nrrthMFAABcXKKhM6e+/fbbNy3/7LPP0iJNAADAbYHG8ePHg95ogQIFQk0PAABwY6Bx9913B71RBu4CAACpCjTWrl1LjgEAADuBRuXKlVO3ZQAA4HpMqgYAANJXrxM1adIk+f7772Xv3r2mm6sv5jsBAAAhl2h89NFH8vLLL0u5cuVk9erV8vjjj0uxYsVk/fr1Urt2bXIWAACEHmiMHTvWlGa888475vkbb7whM2bMMAGIlnAAAACEHGjs2rVLatas6Z335Ny5c+bxM888I4sXLyZnAQBA6IGGtsnwTJ6mVSZafaIOHDggWbKE3OwDAABkMLcdFXTp0kU6dOhgJlNbtGiRtG3bNm1SBgAA3BloHDp0yPt4wIABUqhQIVmxYoVpIJqQkJCW6QMAAG4LNO655x7v46ioKOnWrZu5AQAA3Nakajph2q0mWGNSNQAAENKkajph2q0mWGNSNQAAEPKkakywBgAArE2qdvr0aWnQoEGwbwUAAC4V0jgajRo1onoEAADYCTTi4+Nl69atobwVAAC4SEiBhs5t0rVrV/ntt9/k7NmzcvnyZb8bAABAyONo6JwmqlatWgFfp9cJAAAIOdBYunQpuQcAAOwEGqNGjZLJkycHfO3JJ5+UunXrhrJZAACQwYTURmPKlCnJVplMnTr1dtMEAADcGGjo+Bl6833suZ08eVJmzZolhQsXDikhQ4cOlbx585pJ2pIaN26cVKpUycyxorPEMmAYAAAZsOrkrrvuCvjYI1OmTDJixIiQ2nyMHz9e8ufPL5cuXfJ7bdKkSdKjRw/58ssvTePT999/Xxo2bCibN28OOagBAADpMNBYtWqVua9evbr3sUfWrFklLi5O8uXLl6oEaElI586d5ZtvvpGePXve9PqwYcNML5d//OMf5vno0aNl+vTpMmbMGBkyZEiqPgsAAKTjQKNatWrmfteuXWbQrrTw7LPPSqdOnQI2INUqmY0bN8rAgQO9yzJnzmxKNH799dc0+XwAAJDOep2kVZAxcuRIOXjwoEybNi3g6/qaKlSokN/yggULypo1a5LdbmJiorl56KBiAAAgQnqdpIUNGzbI4MGDZeLEiZIlS8rxjrb98KXrpzQomFa3xMbGem9apQMAAFwUaCxZskTOnDlj2ntobxO9aTWJtr3Qx9evXzclF+r48eN+7z169Kj3tUC054pu23Pbt2+f9f0BAADpKNDo3r27HDt2THbv3u29lStXzrTZ0MfaFqNAgQJSsmTJm0Yi1SDloYceSnbb0dHREhMT43cDAAAuCjSyZcvmLcnw3DS40CBBH3v06tVLvvjiC1m2bJlcuXLFVIscPnxYnn/++XAlHQAA2GwMeiclJCSYqpMWLVrIhQsXTENU7d5aunTpcCcNAABEUqChVSRJG4ZGRUWZ8TK04ahOQZ8jR46wpS9SxfefHe4kAABcKl0FGnny5En2NQ04CDIAAIgsYWujAQAAMj4CDQAAYA2BBgAAsIZAAwAAWEOgAQAA3NHrBO7uTrt7eMs0TwsAILwo0QAAANYQaAAAAGsINAAAgDUEGgAAwBoCDQAAYA2BBgAAsIZAAwAAWEOgAQAArCHQAAAA1hBoAAAAawg0AACANQQaAADAGgINAABgDYEGAACwhkADAABYQ6ABAACsIdAAAADWEGgAAABrCDQAAIA1BBoAAMAaAg0AAGANgQYAALCGQAMAAFhDoAEAAKwh0AAAANYQaAAAAGsINAAAgDUEGgAAwBoCDQAAYA2BBgAAsCaLhJHjODJ//nxZvny5ZMmSRerWrSsPP/zwTesdOnRIvv32Wzly5IhUrFhROnbsKFmzZg1LmgEAQASUaNy4cUOqVKkiH3/8sWTOnFkuXLgg7dq1k+7du/utt3XrVhNcLFy4UPLkySNDhw6VJk2ayPXr18OVdAAAkN5LNKKiouSbb74xQYRHo0aNpHHjxvLiiy9KpUqVzLJ+/fpJhQoVZPbs2eY93bp1k/vuu08mTpwonTt3DlfyAQBAei7R0KDBN8hQ5cuXN/eHDx8291evXpW5c+eaqhJdXxUtWlQaNGggP/30UxhSDQAAIqaNRlLjxo2TnDlzSrVq1czzvXv3SmJiopQoUcJvvZIlS8qyZcuS3Y6+R28eZ8+etZhqAACQ7nudLFiwQAYNGiQffvih5M+f3yy7dOmSude2Gb5iYmLk4sWLyW5r2LBhEhsb673FxcVZTj0AAEi3gcbSpUulTZs2MnDgQHnhhRe8y3Pnzm3uT58+7bf+qVOnTLCRnAEDBsiZM2e8t3379llMPQAASLeBhlaBtGjRQl555RV58803/V4rVqyYCTa2bNnit1yflytXLtltRkdHm0DE9wYAAFwWaOj4Gc2aNZM+ffrI4MGDb3o9U6ZM0r59e/nqq6+81Sjr1q0z7+vQoUMYUgwAACKiMej58+elefPmpvHniRMnpGfPnt7XOnXqJLVq1TKPhw8fbnqZVK1aVSpXriw///yzdO3aVVq1ahWupAMAgPQeaOhIoDr4ViDagNOjYMGCsnbtWpk3b54ZGfSll16SmjVr3sGU4k6J7z87pPftHt4yzdMCAIjwQCN79ux+pRgp0TYXrVu3tp4mAACQwRqDAgCAjItAAwAAWEOgAQAArCHQAAAA1hBoAAAAawg0AACANQQaAADAGgINAABgDYEGAACwhkADAABYQ6ABAACsIdAAAADWEGgAAABrCDQAAACBBgAAiDyUaAAAAGsINAAAgDVZ7G0auDPi+88O6X27h7dM87QAAPxRogEAAKwh0AAAANZQdQLXCrXKRVHtAgDBoUQDAABYQ4mGS67AgTtxrFHSAyApSjQAAIA1BBoAAMAaqk6AO4gqCQBuQ4kGAACwhhINIAOjATGAcKNEAwAAWEOgAQAArCHQAAAA1hBoAAAAa2gMCkQAGnUCiFSUaAAAAGso0QBCQAkDAASHEg0AAGBNxAQa+/btkz/++EPOnj0b7qQAAICMUnVy+fJl6dSpk8ydO1eKFy8ue/bskREjRkhCQoJEKordgbT5TTAtPZD+pftAY/DgwbJy5UrZsWOHFC5cWKZPny5t27aVGjVqyEMPPRTu5AEAgEgONMaPHy89evQwQYZq06aNVKhQwSwn0ADSl0gprbvT6aTkJXl8Fxm/lC9dBxoHDx6UI0eOSNWqVf2Wa2nG2rVrk31fYmKiuXmcOXPG3Kd1+44biRfTdHsAUifU3/Sd/u3Stix5fBdpmy82jjXPNh3HyXiBxsmTJ819/vz5/Zbrc89rgQwbNsxUuSQVFxdnIZUAwiX248jI+0hJpxtk9O8i1uL+nTt3TmJjYzNWoJE1a1Zvg1Bfly5dkmzZsiX7vgEDBkifPn28z2/cuGECEw1QoqKi0izC08BFe8PExMSkyTYzGvKIfOJY4veW3nBeSn0eaUmGBhlFihSRUKTrQEN3NFOmTHLgwAG/5fq8WLFiyb4vOjra3HzlzZvXShr1SyDQII84lu4Mfm/kEcdSeH5voZRkRMQ4Gjlz5pTatWvLjBkzvMsuXLggCxYskMaNG4c1bQAAIMJLNNQ777xjggqtDqlVq5aMHDlSChYsKN27dw930gAAQCSXaKj69evLwoULzUBdn3zyiZQvX15+/fVXyZ07d1jTpVUzb7311k1VNCCPOJb4vYUD5yTyKb0eS1FOqP1VAAAAIr1EAwAARC4CDQAAYA2BBgAAcG+vk/RIBy7566+/pECBAhIfHx/u5KQL2u143bp1ZuyT5MY4OXXqlGzfvl3uvffekAd+ifTjRvdf5+255557Aq5z/fp12bRpkxlYThs+6zgybqMNv3XagBIlSkiePHkCrrN//345fPiwlCpVytoYOZFgzZo1ZgDDOnXq3PTaxYsXZcuWLSZ/7rvvPnGTDRs23DQUt/ZWvP/++29aV8/lev7S35tbG/dv377dTNvxwAMPBDznnDhxQnbu3GnO78mdu1KkjUERvM8//9zJmTOnU6ZMGSdXrlxO8+bNnfPnz7s2Cw8ePOgkJCQ4hQsXdqKjo5233nor4HpDhw41r5crV87Jnj278/TTTztXr1513GD37t1Ohw4dnLx58zqVK1d2YmJinEceecQ5dOiQ33pr1qxxihcv7tx7773OPffc49x3333Oxo0bHbeYPXu2OT5KlSrllC9f3smRI4fTr18/v3USExNNXuprDzzwgDmW3n//fceNNL8yZ85sbklNnjzZHGelS5d28uTJ49SvX985efKk4xYPPfSQU6xYMadOnTre25AhQ246d1WtWtW56667nJIlSzr58+d35syZ47jJ6tWrnQoVKjiFChUyeVGxYsWbzjmvvfaa99yt9927d3euX7+eqs8h0EiFDRs2OJkyZXK+++478/zo0aNOfHy8+aN1q2XLljmffPKJc+rUKfPHGCjQ+Pnnn83JcOHCheb5jh07nHz58jkjRoxw3ED3e8qUKd4fp+aV/qhbtWrlXefKlSvmZNelSxfz/MaNG0779u2dsmXLpvpHHclB/K5du/yOLb0WmjlzpnfZoEGDTBC2d+9e81z/GKKiopzFixc7bnLgwAGnaNGiTq9evW4KNHbu3Olky5bN+eyzz8zzM2fOmD8JDe7dFGi8/fbbKa7TrFkzp3bt2s6lS5fMcz13aXB27Ngxxw0OHjxozsM9e/Z0rl27ZpZt27bNmTdvnnedadOmmWNpxYoV5vnmzZtN4Dp69OhUfRaBRir06dPHXG35Gj58uBMbG+v9otwsuUDjiSeecBo0aOC3TA9uLRVyqw8++MBcSXnMnz/f/Klu377du2zdunVm2dKlSx030gBLSy7Gjh3rXaZXqf379/dbr1q1at4AzS358vDDDzsff/yxCc6SBhp65V6wYEG/AFXzUK9G3VL6qoHGv/71L2flypXO/v37AwZqGqBOnz7du0zzRo+3MWPGOG7Qt29fE7RrKWFyWrRo4Tz66KN+y7p27WoulFLDfRXAt0Gnpg80Zb3WJ2v9FVKXb1u3bjX1yG60atUq077AN490LgHfuvRKlSqZyQP1NbfQ35IOyDd37lzp1KmTqVPv0KGDeU0nRty7d2/AY8lNeTR06FBzXLz00ksBX9e8qFKlil9du+aR1sFv3rxZ3OKzzz6T5557TsqVKycPPvigabfhoe3J9ELb91jKlSuXaaPglmPpl19+kaZNm5r2YNrWZ9euXWYC0mDO3ZqXSddNCY1BU0FPdHog+vJMYZ/StPVu55k515c+1x+6NhDVOW3c5D//+Y9MnTpVfvrppxTzSOkyNx1bGrD379/fND47dOiQjBgxwtvY05MPgY4lt+SRBmGjR482fwDJzUSteaENrt18nnrxxRfl8ccflxw5csj58+dN0NqmTRv5888/zfmGY0nk4MGDpsG1BmIaZOnv7e6775ZJkyZJxYoVUzx3X716NVVTxlOikcpp6wNNWa9Smrbe7ci3/0+H09eT3rvvviutWrVKMY88x5ebji29Etc/U+0tMWfOHOndu7dMmDDBm0cq0G/QLXnUsWNH6dKli+zYscPkk/YWUPpY/zgUvzeRp59+2gQZSqer+PDDD80V+++//+7NI7cfS1mzZpXZs2ebix4t4dEp4UuWLGmOsbT+zyPQSIXixYsHnLJepTRtvdsll296ZRHoKj6jWrx4sQkuXnvtNXPVnjSPjh8/LleuXPEu0y53WpXg1mNLZ27WbptajaK0S7Se+AIdS27JI93PZcuWmeNHbz/88IPpEq2PdbniPHWzQoUKmXvPsaN55PvcjcdSfHy8qQbR4N4TODzzzDOm1EdLmlM6lnRoB08gFwwCjVTQWWSXLFliTv4eWvytX5Sb/jBDybeff/7ZFLf55tsjjzzimnEili5dKi1btpR+/frJG2+8cdPrmheaP5pPHjNmzDD507BhQ3GDpO11rl27Jrt37/b+tjTI0EkWNV889Gpr3rx55hhzAy258L3p8ZQ5c2bzuH379mYdzYs//vjDFIX7/t60mNwN42noFXfS9gPz58839xUqVDD32u4gX758fsfSxo0bTamHW46lpk2bmmPEd7ozHZ8me/bs3vFrNC+01EODWd9jKdV5FHKTVRfSblDaTaxevXqmtfKbb75pWnzPnTvXcavLly+bXhF60/EfnnnmGfNYuwJ7aHexIkWKOI899pgzY8YM0+NExz/QPtxu8Mcffzi5c+c23Vk9eeW5aTdWjx49ephW4F9//bUzfvx4069fezq5gfba0t/Whx9+aLpDT5061WnatKnpmfPXX39519Nudtrd7pVXXnF++uknM46N9kTRLsNuFKjXieZl9erVze3HH390hg0b5mTJksV0sXaD9evXOzVq1DA9bbSr5nvvvWfGsEnavVe7aOp5SLvn6/Gm47I0adLEcYvTp087JUqUMPmivzntbaPdXV9//XXvOtpjp0CBAqbnoJ67//nPf5pz2aZNm1L1WczemkpavK0N1LQxll5pvfDCC/Lwww+LW+nojNroKikt5Rk5cqT3ufYW0Hz73//+Zxqq9erV66bWzBnV999/L5988knA1/773/966zr1qmHMmDHmCkIb+rVu3Vq6d+/umlKfI0eOyKhRo2T16tWmWFZ73fTo0cM0UPOl9ex6bOnVmI7mqFf1SRs/usWsWbPkgw8+kEWLFvkt11LX9957z+SVNqbt1q2bNG/eXNxCi//Hjh1rerZplZv+ltq1axfwtzlx4kRTmlavXj3p06ePqxqnHzlyxBwn69evN78zbTDr6eXl20Bb19m2bZupVtI88jQWDRaBBgAAsMYdl0oAACAsCDQAAIA1BBoAAMAaAg0AAGANgQYAALCGQAMAAFhDoAEAAKwh0EC6onN9TJ48WU6fPi2RLtz7osMw6+frTKiRxnbe2dq+zk+j29UZQ8MllH0L9j3Hjh0zA8qFSgcS09mL4S4EGkhXzp49K0899ZSZ4yJYOt+FniR956BJD8K9L/rnoZ+vI/pFmlDy7k5vP9B3pX/Eul0dMTctpea4CGXfgn2Pjgq5atUq73OdKlznmtHZdpPSKcZ95xJRMTExMmjQIPnxxx+DThsiH4EG0pXo6GgzBO5dd90V9Hv0KkxPkjrNcaTLSPuS0Y6zcH5XqfmstNi3QHTSMQ0QevfubZ7rDJ+VK1eW4cOHS4MGDeTzzz/3W//ll1+Wv/76y2+ZDq2vQ8brDMa+k3khY8sS7gQASU+SOt5+bGystyh65syZ8uijj8rRo0fNlVPRokXNPBhKT1Y6m6DSmU91jgMds19nQ/Vcqa1YscI81pOiZ7ropNvWKzmdF0HnX9EpkD3LDx48aE6WOqVyoPH99cp15cqVJt06rbln1sPkTJ061VRp6IybOgWzzgmjs5Le7r54aOnF5s2bzSydpUqVuuXBpTPG/vbbb2bbmqeat8GkNWn+aT7od6NzIXi+mw0bNpjZMHU+Et+0+L7vVvmb1K3yIKX9uZ3jLKnkvivf2VF1joiUthPM95nSZ9WsWTPg8Vu4cGG/fQvmuwzG6NGj/bY7fvx4851Nnz7dzE6scx4999xz3tlSN23aJOPGjbtpO23btpXnn39efvnlF2nUqFGq0oAIZWFSOCBkOtOrHpZr1641z3ft2mWe68ynOrtny5YtzeyBCQkJ5vXr16+bWWF1nWbNmjkdOnRw3nrrLfPatGnTzOyfOtuuvhYbG+uMGjXK+1mebT/66KNOmTJlnPbt25sZVT3LdWbQUqVKmRkdc+bM6fTq1csvrSNHjnRy5MjhNGjQwKlSpYqZ+XDhwoXJ7ovq2LGjSePf//53p3Tp0maf9u7de9v7ooYPH+5ER0c7jzzyiFOhQgWTft2WzngayL59+5ySJUuadVu3bm0eDx48OKi0+uZf/fr1nYoVK5rPy5o1q5lx9qmnnnIqVapk0qppmjBhwk3vSyl/A+XdrfLgVvtzO8dZUsl9V8FuJ5jvM9jPSnr8pva4Sy6/kypatKjzxRdfeJ+/8MIL3vzVWZz1/WfOnHEuXLhgPiOl2ZkbNWp00+8JGReBBtKV5P4Ann32We+U6osWLXKioqKcPXv2mOeHDh0y62zcuNG7nZ07dzq5cuUy63qsWrXKTAu9ZcsWv23rFMh6MvfwLNeT4ZUrV8yy33//3cmUKZOzZMkS83zbtm3mT3Xy5Mne9+mfSXx8vDnpBtqXpPQz9cTftWtX77JQ92Xr1q1muvCZM2ea55pX+mefUqAxaNAgp3bt2n7pmT59etBp9eSTTh3t+W50qmld1rNnT+9677//vhMXF5eq/E2ad8HkQWr2J5TjLKlA31Uw2wlmX1LzWUmP31COu1u958iRI+b15cuXe5fpdOI9evTwpiVbtmxm2xpo9u3b1wQcOk37hg0bbtpe7969/b4rZGxUnSAiaFGr1u+qunXrmqnTtahYi+oDmTRpkuTPn980ztOpoD31wfny5TPFvGXLlvWu++KLLwacil3rmD3FyzVq1DD10FoE/X//93/yww8/SFxcnN+Uyq+//rqZvlyn5tYpp5OjVRtarK49E7QIfMmSJSnuezD7oukpXbq0KUZXmld9+/Y1702OTsWuPVIOHDhgplnXPHjsscdSnVadyt7z3dSqVcu7zEOXaVouXbpkPjOY/A0lD4LZn7Q+zkLZTmqOzWAkd/ze7nHn6/jx4+bet93HE088YY7zMmXKyJw5c0wbkrVr15qp65ctWybVq1c3+7l9+3Z59dVXzfftodvxbBMZH4EGIoKehD20njlLliymJX5ytM5ae11MmzbNb7n+iWkbDF9apx2IthvwVaJECdmzZ495rPclS5b0e13r2HPlyuVdJyn9o23VqpU5GetJWOu6tXGftglISTD7snfv3oDpTUmPHj1k3bp1JkApV66cNG7cWHr27Gn+pFOTVt8/H237kNyyxMREv0AjpfwNJQ9S2h9bx1ko20nNsRmM5I7f2z3ufOXOndvblsXjb3/7m2mfob1hNEhMSEgw+6BtObQLq6ZrwYIFsn79ehOQvPTSS96ASLdzq/ZMyDgINJAhaTc6/bPTk+CteK48k9JW9Umfe/4I9F5LLnzpH+nFixeT/bP45ptvTGNN/TP1nLhHjRpluvvd7r7oleOaNWtSTH9SeqL/7rvvzElfr0C1NKZatWrmCnTixIkhpTU1UsrfUPIgpf3RADASj81gJHf83u5x56tIkSLmvdq4VxucemiAoTelvU+00as28NQSrAcffNAs12Ua7GiDYd2O0u1oSQjcge6tiHiek6fvlWezZs1Mi//Fixf7rav9/vUWDL1a89DxC7SVfJ06dbzF4XqltmPHDu86eoWaPXt27wk2KT3R6tW1J71KqzzSYl80PRpo+HZ/vNVYBVrFoPRPuEmTJvLpp5+aNOofUjBpvV0p5W9SweRBSvtjQ6DvKhihHJuhfpZKi+9SS2S0VEIDuEA0mBs7dqx89NFH3sDXM5aIVpFcu3bNr5Rr+fLl9DhxEUo0EPH0BKpXRyNGjDBd57QKQ6+qtJ2AtlnQIl3tdqjdKLVIV7veBVNsq1eCeoLUOnMdI0DbZHTu3Nm8pn9kLVq0MMXzWvesf5T6+QMHDky2m6Km5e233zbrV6hQwaRFi/q1aP1296Vp06bmj0DX1XX0JK/dD1MyZswY071S9yNv3rymFEOLw++///6g0nq7UsrfpILJg5T2x4ZA35Vv99bkhHJshvpZKq2+S+26qlVRH3zwwU3v1f3REg1PlZF2ddWSDw08Fi1aJO3atfNWm2mwol1727dvn6rPR+SiRAPpStLBhvTqVJ8nPflqQzTfundtgKZtJrRRmjaoU//+979NKYM2ftMrKD056x+Rp+1Cctv20G1pUa+Ok6GNCnW7vmMP6FWhNgDVk7aOBzFlyhQZMGBAsvuiJR26Df1z1ZNty5YtTamDnoR9hbIvSsdU0D+D1atXm5IVfV0/P7nqiHfeeccERocOHTJjT+gfmH6eXr0Gk9ZA+afVArosZ86c3mX656PLsmXLFnT+Bhp06lZ5kNL+pNVxllTS7yrY7QTzfYb6WaEcd8EM8tW6dWvTiDRp2xIdL0Xbfjz55JPeZTpuyty5c02Vjb729ddfe1/TKi0NenxLWJCxRWnXk3AnAkhPtLGenvD1JBnMoFcgf91CAzgNdIYMGRLS+7Udjo4sqqVPvg2DkbFRdQIACIqORqq3UGmJyVdffUVuuwxVJ0ASt6pSwe0hfwF3oeoEAABYQ4kGAACwhkADAABYQ6ABAACsIdAAAADWEGgAAABrCDQAAIA1BBoAAMAaAg0AAGANgQYAABBb/h9Hchf4nAIF7QAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "recon = np.ma.getmaskarray(seg_strict[\"left_pupil\"]).mean(axis=0) # one value per trial\n", "fig, ax = plt.subplots(figsize=(6, 3))\n", "ax.hist(100 * recon, bins=30)\n", "ax.set_xlabel(\"interpolated samples in the trial (%)\"); ax.set_ylabel(\"trials\");" ] }, { "cell_type": "markdown", "id": "ae2e44d9", "metadata": {}, "source": [ "`Epochs.keep()` selects the trials to keep, e.g. those with less than 25%\n", "interpolated samples:" ] }, { "cell_type": "code", "execution_count": 23, "id": "3dbcf252", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:42.161243Z", "iopub.status.busy": "2026-11-09T20:27:42.160895Z", "iopub.status.idle": "2026-11-09T20:27:42.166065Z", "shell.execute_reply": "2026-11-09T20:27:42.165105Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "300 trials, 284 with less than 25% interpolated samples\n" ] } ], "source": [ "good = stim.keep(recon < 0.25)\n", "print(len(stim), \"trials,\", len(good), \"with less than 25% interpolated samples\")" ] }, { "cell_type": "code", "execution_count": 24, "id": "f36fc380", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:42.168234Z", "iopub.status.busy": "2026-11-09T20:27:42.167778Z", "iopub.status.idle": "2026-11-09T20:27:45.416583Z", "shell.execute_reply": "2026-11-09T20:27:45.415444Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "largest difference in a trial mean: 48.1 (median 3.13)\n" ] } ], "source": [ "mean_default = d.stat_per_event(stim, variables=\"pupil\")[\"left_pupil\"]\n", "mean_strict = d.stat_per_event(stim, variables=\"pupil\", policy=\"strict\")[\"left_pupil\"]\n", "diff = np.abs(np.asarray(mean_default) - np.asarray(mean_strict))\n", "print(f\"largest difference in a trial mean: {np.nanmax(diff):.1f} \"\n", " f\"(median {np.nanmedian(diff):.2f})\")" ] }, { "cell_type": "markdown", "id": "e73e2b00", "metadata": {}, "source": [ "## Downsampling\n", "\n", "`downsample()` averages the usable samples in each new sampling interval.\n", "The layers are carried over: a new sample is in a layer if the layer\n", "covers more than half of its interval. The proportion of samples in each\n", "layer therefore stays about the same:" ] }, { "cell_type": "code", "execution_count": 25, "id": "64e0d2d6", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:45.419624Z", "iopub.status.busy": "2026-11-09T20:27:45.419443Z", "iopub.status.idle": "2026-11-09T20:27:45.653259Z", "shell.execute_reply": "2026-11-09T20:27:45.652566Z" } }, "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", "
500 Hz50 Hz
NaN0.000.00
blink5.845.78
offscreen0.000.00
divergence0.000.00
excluded0.350.35
interpolated7.497.45
masked (default)0.350.35
masked (strict)7.847.80
\n", "
" ], "text/plain": [ " 500 Hz 50 Hz\n", "NaN 0.00 0.00\n", "blink 5.84 5.78\n", "offscreen 0.00 0.00\n", "divergence 0.00 0.00\n", "excluded 0.35 0.35\n", "interpolated 7.49 7.45\n", "masked (default) 0.35 0.35\n", "masked (strict) 7.84 7.80" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dd = d.downsample(50, inplace=False)\n", "pd.DataFrame({\"500 Hz\": layer_table(d).loc[\"left_pupil\"],\n", " \"50 Hz\": layer_table(dd).loc[\"left_pupil\"]})" ] }, { "cell_type": "markdown", "id": "08a19154", "metadata": {}, "source": [ "## Combining the eyes\n", "\n", "`d.pupil.merge_eyes()` computes `mean_pupil` from the usable samples of\n", "both eyes, recorded samples first: where one eye was interpolated and the\n", "other one recorded, the recorded value is used (after the offset\n", "correction). The merged channel gets the layers that the eyes used for a\n", "sample have in common: it is `\"interpolated\"` only where both eyes were\n", "interpolated, and it is masked where neither eye is usable:" ] }, { "cell_type": "code", "execution_count": 26, "id": "41725cd5", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:45.656986Z", "iopub.status.busy": "2026-11-09T20:27:45.656799Z", "iopub.status.idle": "2026-11-09T20:27:45.913647Z", "shell.execute_reply": "2026-11-09T20:27:45.912636Z" } }, "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", "
NaNblinkoffscreendivergenceexcludedinterpolatedmasked (default)masked (strict)
left_pupil0.05.840.00.00.357.490.357.84
right_pupil0.04.080.00.00.355.660.356.01
mean_pupil0.03.330.00.00.354.840.355.19
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence excluded interpolated \\\n", "left_pupil 0.0 5.84 0.0 0.0 0.35 7.49 \n", "right_pupil 0.0 4.08 0.0 0.0 0.35 5.66 \n", "mean_pupil 0.0 3.33 0.0 0.0 0.35 4.84 \n", "\n", " masked (default) masked (strict) \n", "left_pupil 0.35 7.84 \n", "right_pupil 0.35 6.01 \n", "mean_pupil 0.35 5.19 " ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dm = d.pupil.merge_eyes(inplace=False)\n", "layer_table(dm).loc[[\"left_pupil\", \"right_pupil\", \"mean_pupil\"]]" ] }, { "cell_type": "markdown", "id": "d70dd48d", "metadata": {}, "source": [ "Because `mean_pupil` keeps the `\"interpolated\"` layer where both eyes\n", "were interpolated, the strict policy leaves these samples out in the\n", "following steps (column \"masked (strict)\"). Merging with\n", "`policy=\"strict\"` gives the same channel here:" ] }, { "cell_type": "code", "execution_count": 27, "id": "9dd659ed", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:45.916557Z", "iopub.status.busy": "2026-11-09T20:27:45.916356Z", "iopub.status.idle": "2026-11-09T20:27:46.210098Z", "shell.execute_reply": "2026-11-09T20:27:46.209182Z" } }, "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", "
NaNblinkoffscreendivergenceexcludedinterpolatedmasked (default)masked (strict)
mean_pupil0.03.330.00.00.354.840.355.19
\n", "
" ], "text/plain": [ " NaN blink offscreen divergence excluded interpolated \\\n", "mean_pupil 0.0 3.33 0.0 0.0 0.35 4.84 \n", "\n", " masked (default) masked (strict) \n", "mean_pupil 0.35 5.19 " ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dm = d.pupil.merge_eyes(policy=\"strict\", inplace=False)\n", "layer_table(dm).loc[[\"mean_pupil\"]]" ] }, { "cell_type": "markdown", "id": "5e6eb5a1", "metadata": {}, "source": [ "## Missing data in the quality report\n", "\n", "The quality report counts missing data in the *values* (NaN) and reports\n", "the mask layers in its \"Artifacts\" section; excluded periods are left out\n", "of all measures. Interpolation fills the NaN samples, so after\n", "preprocessing the pupil shows no missing data any more. Run the report on\n", "the raw data to judge the recording, and after preprocessing to see what\n", "was masked and interpolated:" ] }, { "cell_type": "code", "execution_count": 28, "id": "e8b5e454", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:46.212205Z", "iopub.status.busy": "2026-11-09T20:27:46.211984Z", "iopub.status.idle": "2026-11-09T20:27:47.937140Z", "shell.execute_reply": "2026-11-09T20:27:47.936215Z" } }, "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", "
rawpreprocessed
pupil.left.missing_pct3.1388820.0
pupil.left.dropout_rate15.4230760.0
gaze.left.missing_pct3.1388823.138882
artifacts.left.pupil_masked_pctNaN0.0
artifacts.left.gaze_masked_pctNaN4.795098
\n", "
" ], "text/plain": [ " raw preprocessed\n", "pupil.left.missing_pct 3.138882 0.0\n", "pupil.left.dropout_rate 15.423076 0.0\n", "gaze.left.missing_pct 3.138882 3.138882\n", "artifacts.left.pupil_masked_pct NaN 0.0\n", "artifacts.left.gaze_masked_pct NaN 4.795098" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raw = pp.get_example_data(\"rlmw_002\", cache_dir=\"example_data\")\n", "raw.keep_intervals([(first - 1000, last + 3000)])\n", "keys = [\"pupil.left.missing_pct\", \"pupil.left.dropout_rate\", \"gaze.left.missing_pct\",\n", " \"artifacts.left.pupil_masked_pct\", \"artifacts.left.gaze_masked_pct\"]\n", "pd.DataFrame({\"raw\": pd.Series(raw.quality_report().to_dict(flat=True)).reindex(keys),\n", " \"preprocessed\": pd.Series(d.quality_report().to_dict(flat=True)).reindex(keys)})" ] }, { "cell_type": "code", "execution_count": 29, "id": "917bbd25", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:47.939533Z", "iopub.status.busy": "2026-11-09T20:27:47.939323Z", "iopub.status.idle": "2026-11-09T20:27:48.307475Z", "shell.execute_reply": "2026-11-09T20:27:48.306566Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "
\n", "

Quality report — rlmw_002\n", " (EyeData)

\n", " good\n", "
\n", "\n", "
Recording
\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", "
duration19.9 min
included19.8 min
sampling rate500 Hz
samples597301
events1320
time gaps0
\n", "
Pupil
\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", "
leftright
missing data0.0% good0.0% good
dropouts0 /min good0 /min good
longest dropout0 s good0 s good
noise0.027% good0.032% good
noise (SD)0.283 0.286
jumps0 /min good0 /min good
\n", "
Gaze
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
leftright
missing data3.1% good2.1% good
off-screen1.0% good1.1% good
precision0.0185 deg good0.0186 deg good
\n", "
Binocular
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
left/right
pupil correlation0.943 good
gaze disparity0.418 deg good
\n", "
Blinks
\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", "
leftright
number321 307
blink rate16.2 /min good15.5 /min good
mean duration217 ms 159 ms
time in blinks5.9% good4.1% good
long (>500 ms)14 10
\n", "
Artifacts
\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", "
leftright
pupil masked0.0% good0.0% good
interpolated7.5% good5.7% good
gaze masked4.8% good4.0% good
masked off-screen0.97% 1.0%
masked divergence0.89% 0.89%
other artifacts0.0% 0.0%
\n", "\n", "\n", "
\n", " Thresholds\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pupil.missing_pctgood <= 10, ok <= 20, poor <= 40, else bad
pupil.dropout_rategood <= 30, ok <= 60, poor <= 150, else bad
pupil.longest_dropoutgood <= 2, ok <= 5, poor <= 30, else bad
pupil.noisegood <= 0.5, ok <= 1, poor <= 2, else bad
pupil.jump_rategood <= 1, ok <= 10, poor <= 100, else bad
gaze.missing_pctgood <= 10, ok <= 20, poor <= 40, else bad
gaze.dropout_rategood <= 30, ok <= 60, poor <= 150, else bad
gaze.longest_dropoutgood <= 2, ok <= 5, poor <= 30, else bad
gaze.offscreen_pctgood <= 5, ok <= 15, poor <= 30, else bad
gaze.precision_deggood <= 0.05, ok <= 0.1, poor <= 0.2, else bad
binocular.pupil_correlationgood >= 0.9, ok >= 0.8, poor >= 0.6, else bad
binocular.gaze_disparity_deggood <= 1, ok <= 2, poor <= 3, else bad
blinks.rategood <= 20, ok <= 35, poor <= 70, else bad
blinks.time_pctgood <= 10, ok <= 20, poor <= 35, else bad
artifacts.pupil_masked_pctgood <= 10, ok <= 20, poor <= 35, else bad
artifacts.gaze_masked_pctgood <= 10, ok <= 20, poor <= 35, else bad
artifacts.interpolated_pctgood <= 10, ok <= 20, poor <= 30, else bad
calibration.mean_errorgood <= 0.5, ok <= 1, poor <= 1.5, else bad
calibration.max_errorgood <= 1, ok <= 1.5, poor <= 2, else bad
\n", "
The overall label is the worst label of all measures. Hover over a measure for a description.
\n", "
\n", "
" ], "text/plain": [ "Quality report: rlmw_002 (EyeData)\n", "Overall: \u001b[32mGOOD\u001b[0m\n", "\n", "Recording\n", " duration 19.9 min\n", " included 19.8 min\n", " sampling rate 500 Hz\n", " samples 597301\n", " events 1320\n", " time gaps 0\n", "\n", "Pupil left right\n", " missing data 0.0% \u001b[32mgood\u001b[0m 0.0% \u001b[32mgood\u001b[0m\n", " dropouts 0 /min \u001b[32mgood\u001b[0m 0 /min \u001b[32mgood\u001b[0m\n", " longest dropout 0 s \u001b[32mgood\u001b[0m 0 s \u001b[32mgood\u001b[0m\n", " noise 0.027% \u001b[32mgood\u001b[0m 0.032% \u001b[32mgood\u001b[0m\n", " noise (SD) 0.283 0.286\n", " jumps 0 /min \u001b[32mgood\u001b[0m 0 /min \u001b[32mgood\u001b[0m\n", "\n", "Gaze left right\n", " missing data 3.1% \u001b[32mgood\u001b[0m 2.1% \u001b[32mgood\u001b[0m\n", " off-screen 1.0% \u001b[32mgood\u001b[0m 1.1% \u001b[32mgood\u001b[0m\n", " precision 0.0185 deg \u001b[32mgood\u001b[0m 0.0186 deg \u001b[32mgood\u001b[0m\n", "\n", "Binocular left/right\n", " pupil correlation 0.943 \u001b[32mgood\u001b[0m\n", " gaze disparity 0.418 deg \u001b[32mgood\u001b[0m\n", "\n", "Blinks left right\n", " number 321 307\n", " blink rate 16.2 /min \u001b[32mgood\u001b[0m 15.5 /min \u001b[32mgood\u001b[0m\n", " mean duration 217 ms 159 ms\n", " time in blinks 5.9% \u001b[32mgood\u001b[0m 4.1% \u001b[32mgood\u001b[0m\n", " long (>500 ms) 14 10\n", "\n", "Artifacts left right\n", " pupil masked 0.0% \u001b[32mgood\u001b[0m 0.0% \u001b[32mgood\u001b[0m\n", " interpolated 7.5% \u001b[32mgood\u001b[0m 5.7% \u001b[32mgood\u001b[0m\n", " gaze masked 4.8% \u001b[32mgood\u001b[0m 4.0% \u001b[32mgood\u001b[0m\n", " masked off-screen 0.97% 1.0%\n", " masked divergence 0.89% 0.89%\n", " other artifacts 0.0% 0.0%" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.quality_report()" ] }, { "cell_type": "markdown", "id": "2aa7cb76", "metadata": {}, "source": [ "## Saving\n", "\n", "The layers are saved with the object (`write_file()`, `write_file()`),\n", "so a saved dataset keeps the information on why each sample was masked:" ] }, { "cell_type": "code", "execution_count": 30, "id": "e7758da0", "metadata": { "execution": { "iopub.execute_input": "2026-11-09T20:27:48.311041Z", "iopub.status.busy": "2026-11-09T20:27:48.310705Z", "iopub.status.idle": "2026-11-09T20:27:48.879985Z", "shell.execute_reply": "2026-11-09T20:27:48.878999Z" } }, "outputs": [ { "data": { "text/plain": [ "(['blink', 'offscreen', 'divergence', 'excluded', 'interpolated'], True)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import tempfile, os\n", "fname = os.path.join(tempfile.mkdtemp(), \"rlmw_002.pypil\")\n", "d.write_file(fname)\n", "d_loaded = pp.EyeData.from_file(fname)\n", "d_loaded.data.layer_names, np.array_equal(d_loaded.mask(\"left_pupil\", policy=\"strict\"),\n", " d.mask(\"left_pupil\", policy=\"strict\"))" ] }, { "cell_type": "markdown", "id": "33c04399", "metadata": {}, "source": [ "## Summary\n", "\n", "- Samples the eye tracker did not record are NaN in the values. All other\n", " reasons are mask layers, which record *why* a sample is bad. Each\n", " detector writes its own layer on the channels it checks and replaces it\n", " when it runs again.\n", "- Excluded periods are one layer for the whole recording; they mask every\n", " channel but are not counted as missing data.\n", "- The mask policy (`\"default\"`, `\"strict\"`, `\"raw\"`, or a list of layers)\n", " decides which samples the analyses use. By default, interpolated\n", " samples are used, even where another layer is set.\n", "- To mask the pupil where the gaze is bad, copy the gaze layers with\n", " `mask_intervals()`.\n", "- `merge_eyes()` prefers recorded samples and keeps the layers the eyes\n", " have in common, so `policy=\"strict\"` still leaves out samples that were\n", " interpolated in both eyes." ] } ], "metadata": { "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.15.1" } }, "nbformat": 4, "nbformat_minor": 5 }