{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from collections import namedtuple\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import scipy.stats as stats\n",
    "from sklearn.linear_model import Perceptron\n",
    "import sliding1d as sliding\n",
    "\n",
    "import datareader as dr\n",
    "import saccades as sc\n",
    "import trial_envelope as envelope\n",
    "from kw_dunn import kw_dunn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials  = dr.load_trials(\"../01_data/04_formatted\")\n",
    "\n",
    "# representative session\n",
    "subject = \"MLA007518\"\n",
    "session = \"session2019-09-30-001\"\n",
    "\n",
    "in_session = [trial for trial in trials if (trial.subject == subject and trial.session == session)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "figroot = Path(\"../05_figures/prediction/perceptrons\")\n",
    "saved   = False\n",
    "\n",
    "labelsize = 12\n",
    "\n",
    "def savefig(fig, basepath, root=None, dpi=600, as_png=True, as_svg=True):\n",
    "    if saved == False:\n",
    "        return\n",
    "    if root is None:\n",
    "        root = figroot\n",
    "    base = root / Path(basepath)\n",
    "    if not base.parent.exists():\n",
    "        base.parent.mkdir(parents=True)\n",
    "    is_saved = False\n",
    "    if as_png == True:\n",
    "        fig.savefig(str(base) + \".png\", dpi=dpi)\n",
    "        is_saved = True\n",
    "    if as_svg == True:\n",
    "        fig.savefig(str(base) + \".svg\")\n",
    "        is_saved = True\n",
    "    if saved == True:\n",
    "        print(f\"--> saved to: {base}\")\n",
    "    else:\n",
    "        print(f\"--> ***not saved: {base}\")\n",
    "\n",
    "def psign(p, append=True):\n",
    "    if p < 0.001:\n",
    "        return \"***\"\n",
    "    elif p < 0.01:\n",
    "        return \"**\"\n",
    "    elif p < 0.05:\n",
    "        return \"*\"\n",
    "    else:\n",
    "        return \", NS\" if append else \"NS\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single session: collecting data for prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def baseline(vec, timerange):\n",
    "    return np.nanmedian(vec[timerange])\n",
    "\n",
    "def get_weights(vec):\n",
    "    return np.ones(vec.size) / vec.size\n",
    "\n",
    "def in_camel(txt):\n",
    "    return txt[0].upper() + txt[1:]\n",
    "\n",
    "class TimeRange(namedtuple(\"_TimeRange\", (\"fromFrame\", \"toFrame\"))):\n",
    "    @classmethod\n",
    "    def from_range(cls, timerange):\n",
    "        return cls(*[int(v) for v in timerange])\n",
    "    \n",
    "    @classmethod\n",
    "    def from_row(cls, trial_row):\n",
    "        return cls(int(trial_row.FromFrame), int(trial_row.ToFrame))\n",
    "    \n",
    "    def as_slice(self):\n",
    "        return slice(self.fromFrame - 1, self.toFrame)\n",
    "\n",
    "def time_normalized(vec, num=100):\n",
    "    size     = vec.size\n",
    "    origtime = np.linspace(0, 1, size, endpoint=True)\n",
    "    desttime = np.linspace(0, 1, num, endpoint=True)\n",
    "    return np.interp(desttime, origtime, vec)\n",
    "\n",
    "def before_start(trial):\n",
    "    \"\"\"return the time range tuple corresponding to the initial AtEnd phase.\"\"\"\n",
    "    return TimeRange.from_row(trial.states[trial.states.State == \"AtEnd\"].iloc[0])\n",
    "\n",
    "def count_epochs(trial,\n",
    "                 pattern=\"AtEnd <Backward> Left\",\n",
    "                 **kwargs):\n",
    "    return [[1] for rng in trial.get_timeranges(pattern)]\n",
    "\n",
    "def setpoint_asymmetry_from_trial(trial,\n",
    "                                  pattern=\"AtEnd <Backward> Left\",\n",
    "                                  radius_sample=10):\n",
    "    ranges = [TimeRange.from_range(rng) for rng in trial.get_timeranges(pattern)]\n",
    "    if len(ranges) == 0:\n",
    "        return []\n",
    "    else:\n",
    "        sides      = (\"left\", \"right\")\n",
    "        envelopes  = dict((side, envelope.whisker(trial, side=side, radius_sample=radius_sample)) \\\n",
    "                          for side in sides)\n",
    "        base_range = before_start(trial)\n",
    "        baselines  = dict((side, baseline(envelopes[side].raw, base_range.as_slice())) \\\n",
    "                         for side in sides)\n",
    "        subtracted = dict((side, envelopes[side].bottom - baselines[side]) \\\n",
    "                         for side in sides)\n",
    "        asymmetry  = subtracted[\"right\"] - subtracted[\"left\"]\n",
    "        return [asymmetry[rng.as_slice()] for rng in ranges]\n",
    "\n",
    "def eye_position_from_trial(trial,\n",
    "                            pattern=\"AtEnd <Backward> Left\",\n",
    "                            radius_sample=10):\n",
    "    ranges = [TimeRange.from_range(rng) for rng in trial.get_timeranges(pattern)]\n",
    "    if len(ranges) == 0:\n",
    "        return []\n",
    "    else:\n",
    "        sides = (\"left\", \"right\")\n",
    "        raw   = dict((side, envelope.interpolate(trial.tracking[f\"{side}_pupil_normalized_position\"])) \\\n",
    "                          for side in sides)\n",
    "        if all(np.isnan(raw[\"left\"])):\n",
    "            return []\n",
    "        smo   = dict((side, sliding.nanmean(raw[side], rad=radius_sample)) \\\n",
    "                     for side in sides)\n",
    "        avg   = (smo[\"left\"] + smo[\"right\"]) / 2\n",
    "        return [avg[rng.as_slice()] for rng in ranges]\n",
    "\n",
    "def bin_trace(trace, nbin=5, reduce=None):\n",
    "    \"\"\"bin single-trial traces based on time, and return the list of split traces.\"\"\"\n",
    "    if reduce is None:\n",
    "        reduce = lambda trace: trace\n",
    "    t = np.linspace(0, 1, num=trace.size, endpoint=False)\n",
    "    borders = np.linspace(0, 1, num=nbin+1, endpoint=True)\n",
    "    return [reduce(trace[np.logical_and(t>=start, t<stop)]) for start, stop in zip(borders[:-1], borders[1:])]\n",
    "\n",
    "class NoEpochFoundError(ValueError):\n",
    "    def __init__(self, msg):\n",
    "        super().__init__(msg)\n",
    "\n",
    "def collect_binned(trials, fun, nbin=5, reduce=None, merge=None, **kwargs):\n",
    "    epochs = sum([fun(trial, **kwargs) for trial in trials], [])\n",
    "    if len(epochs) == 0:\n",
    "        raise NoEpochFoundError(f\"no epoch found from {len(trials)} trials ({fun.__name__}, {kwargs.get('pattern', '(no pattern info)')})\")\n",
    "    # print((f\"{len(epochs)} epochs found from {len(trials)} trials ({fun.__name__}, {kwargs.get('pattern', '(no pattern info)')})\"))\n",
    "    ret = []\n",
    "    for ibin in range(nbin):\n",
    "        ret.append([])\n",
    "    for epoch in epochs:\n",
    "        for ibin, binned in enumerate(bin_trace(epoch, nbin=nbin, reduce=reduce)):\n",
    "            ret[ibin].append(binned)\n",
    "    if merge is not None:\n",
    "        for ibin in range(nbin):\n",
    "            ret[ibin] = merge(ret[ibin])\n",
    "    return ret\n",
    "\n",
    "def compute_binned_points(trials, base_pattern=\"AtEnd <Backward> {direction}\", nbin=5):\n",
    "    binned = dict()\n",
    "    for key in ('whisker', 'eye'):\n",
    "        binned[key] = []\n",
    "        for ibin in range(nbin):\n",
    "            binned[key].append({})\n",
    "\n",
    "    for direction in (\"Left\", \"Right\"):\n",
    "        pattern = base_pattern.format(direction=direction)\n",
    "        for ibin, vals in enumerate(collect_binned(trials,\n",
    "                                                   setpoint_asymmetry_from_trial,\n",
    "                                                   pattern=pattern,\n",
    "                                                   merge=np.concatenate)):\n",
    "            binned['whisker'][ibin][direction.lower()] = vals\n",
    "        \n",
    "        for ibin, vals in enumerate(collect_binned(trials,\n",
    "                                                  eye_position_from_trial,\n",
    "                                                  pattern=pattern,\n",
    "                                                  merge=np.concatenate)):\n",
    "\n",
    "            binned[\"eye\"][ibin][direction.lower()] = vals\n",
    "    \n",
    "    return binned\n",
    "\n",
    "def colored_hist2d(bin_whisker, bin_eye, nbin=10,\n",
    "                 range_whisker=(-1, 1), range_eye=(-0.2, 0.2),\n",
    "                 cmap=dict(left='bwr', right='bwr_r')):\n",
    "    H = {}\n",
    "    for direction in ('left', 'right'):\n",
    "        h, xedges, yedges = np.histogram2d(bin_whisker[direction],\n",
    "                                           bin_eye[direction],\n",
    "                                           bins=nbin,\n",
    "                                           range=(range_whisker, range_eye))\n",
    "        H[direction] = plt.get_cmap(cmap[direction])(h / h.max()*2 + 0.5)\n",
    "    H = (H['left'] + H['right'])/2\n",
    "    return H, xedges, yedges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAB1CAYAAACf45KrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIn0lEQVR4nO3dT2wU5x3G8efBlENKK4V6SVL+1FRCjciligxK20ukConmYo7JCaWVfAmtKrWVLFWqxa0NinLi4gMql5RbKamoUsSlhwbJBgW1VKV2ERQLGoyQiriEVvn14EEsxt7ZXb87+874+5FGuzvz8puXx8Y/Zndm7IgQAAAYrk3DngAAAKAhAwCQBRoyAAAZoCEDAJABGjIAABmgIQMAkIHNw9z56OhojI2NDXMKjXPjxg3du3fP7evIeTAuXbp0LyJaj1+T82CQczXIuRorc2431IY8Njam2dnZYU6hcfbv3//MOnIejE2bNt1sf03Og0HO1SDnaqzM+altVU4EAACsjoYMAEAGaMgAAGSAhgwAQAZoyAAAZICGDABABmjIAABkgIYMAEAGaMgAAGSAhgwAQAZoyAAAZICGDABABmjIAABkgIYMAEAGaMgAAGSAhgwAQAZ6bsi2T9q+a/uvbeu22T5ve754fD7tNAEAaLZ+jpB/LenQinVTki5ExF5JF4rXAACgSz035Ij4k6T7K1ZPSDpVPD8l6fD6pgUAwMaS6jPkFyLijiQVj9vXGmh70vac7bmlpaVEuwcAoN4qP6krImYiYjwixlutVtW7BwAgS6ka8qe2X5Kk4vFuoroAAGwIqRryWUlHiudHJP0uUV0AADaEfi57+o2kjyV9w/ai7R9I+qWkg7bnJR0sXgMAgC5t7vUPRMRba2z67jrnAgAbz61bnbfv2lXNPJru8uXO2z/8sLzG9HSauayBO3UBAJABGjIAABmgIQMAkAEaMgAAGaAhAwCQARoyAAAZoCEDAJCBnq9D3qiO/exh6Zjp41srmEl9HfvJg/JB779fOmT65vc7D9jg120eO1Y+ZlpdDEphwNdtDlLZt+KRI523S9K2bV3s6OTJrubTUY1zrsyrr3be/txz1cyjA46QAQDIAA0ZAIAM0JABAMgADRkAgAwkbci2D9m+ZnvB9lTK2gAANFmyhmx7RNIJSd+TtE/SW7b3paoPAECTpTxCPiBpISKuR8QjSaclTSSsDwBAY6VsyDsktf9iz8ViHQAAKJHyxiBeZV08M8ielDQpSbt37064+/51c9MPvfde+ZjjXJzfUTc3/SDCdesuwy4Gld1hZKK+b4DNz5ePeVByH5uubvqRAv8o0ij7fn777Wrm0UHKI+RFSe23SNop6fbKQRExExHjETHearUS7h4AgPpK2ZBnJe21vcf2FklvSjqbsD4AAI2V7C3riPif7aOSPpI0IulkRFxNVR8AgCZL+sslIuKcpHMpawIAsBFwpy4AADJAQwYAIAM0ZAAAMkBDBgAgA0lP6qqtrVtLh0x/zsX5Zcquu+/mBgf3f1i+n8puyFBTMzPlY+7cKR8z3eAbUnzwwfprXLxYPua117oo1OCcS38mqMK/fg1y5ggZAIAM0JABAMgADRkAgAzQkAEAyAANGQCADNCQAQDIAA0ZAIAMcB2yanF5Wi2QYx4mJ4c9g2bg+7ncmTPDnkGz9HSEbPtl2x/b/sz2T1dsO2T7mu0F21NppwkAQLP1eoR8X9KPJB1uX2l7RNIJSQclLUqatX02Iv6WYpIAADRdT0fIEXE3ImYl/XfFpgOSFiLiekQ8knRa0kSiOQIA0HipTuraIelW2+vFYt0zbE/anrM9t7S0lGj3AADUW6qG7FXWxWoDI2ImIsYjYrzVaiXaPQAA9VbakG2/Y/uTYvnqGsMWJe1qe71T0u0UEwQAYCMobcgRcSIivlksazXZWUl7be+xvUXSm5LOppwoAABN1tNZ1rZflDQn6cuSPrf9Y0n7IuKB7aOSPpI0IulkRFxNPVkAAJqqp4YcEf/W8tvRq207J+lcikkBwCBx0480rlzpvJ2ce8OtMwEAyAANGQCADNCQAQDIAA0ZAIAM0JABAMgADRkAgAzQkAEAyECvv34RAABJXGecGkfIAABkgIYMAEAGaMgAAGSAhgwAQAZoyAAAZICGDABABmjIAABkgIYMAEAGHBHD27m9JOlm4rKjku7VoOag6n4tIlrtK2qU86DqDmquT2VNzuRcUd065yzVK5NKfkY/NtSGPAi25yJiPPeag6xbhTplQs7V1CXnaurWOWepXplUnTVvWQMAkAEaMgAAGWhiQ56pSc1B1q1CnTIh52rqknM1deucs1SvTCrNunGfIQMAUEdNPEIGAKB2at2Qbb9u+z+2PymWX7RtO2T7mu0F21Pr2EeqOjds/6WY51yxbpvt87bni8fn+60/SHXKuahF1mvXJ2eRc1XIuUcRUdtF0uuSfr/K+hFJ/5T0dUlbJF2RtK+P+knqFLVuSBpdse5dSVPF8ylJvxp2pnXPmazJmZzzWci5t6XWR8gdHJC0EBHXI+KRpNOSJoZYZy0Tkk4Vz09JOpywdhXqkrNE1qlqlCFncu4GOa+iCQ35W7av2P6D7VeKdTsk3Wobs1is61WqOpIUkv5o+5LtyWLdCxFxR5KKx+191q5CXXKWyHot5Pw0cq4GOXdp8yCLV+Cylm9D9tD2G5LOSNoryauM7ed08lR1JOk7EXHb9nZJ523/vc86w1CnnCWyXgs5P0HO1SDnHtTuCNn2O49PEJC0NSIeSlJEnJP0BdujWv6f0q62P7ZT0u0+dpeqjiLidvF4V9Jvtfx2y6e2X5Kk4vFuP7UHoa45F3Mk69WRMzkPHDmvwyA/oB70IulFPbmW+oCkf2n5f02bJV2XtEdPPux/pY/6qep8UdKX2p7/WdIhScf19AkD7w470zrnTNbkTM55LeTc4zyG/QVb5xf7qKSrxRfioqRvt217Q9I/tHwW3s/XsY9119HyWYBXiuXq4zqSviLpgqT54nHbsDOtc85kTc7knNdCzr0t3KkLAIAM1O4zZAAAmoiGDABABmjIAABkgIYMAEAGaMgAAGSAhgwAQAZoyAAAZICGDABABv4P4FMASYjNRHYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "binned = compute_binned_points(in_session)\n",
    "\n",
    "fig, axes = plt.subplots(1, 5, sharex=True, sharey=True, figsize=(8, 2))\n",
    "for ibin, ax in enumerate(axes):\n",
    "    H, xedges, yedges = colored_hist2d(binned['whisker'][ibin], binned['eye'][ibin])\n",
    "    ax.imshow(H.transpose((1,0,2)), aspect=5, origin='lower', extent=(xedges.min(), xedges.max(), yedges.min(), yedges.max()))\n",
    "    ax.tick_params(left=(ibin == 0), labelsize=10)\n",
    "plt.xticks((-0.5, 0, 0.5), (\"-50\", \"0\", \"50\"))\n",
    "plt.yticks((-0.1, 0, 0.1), (\"-10\", \"0\", \"10\"))\n",
    "plt.subplots_adjust(wspace=.05)\n",
    "savefig(fig, f\"2d-hist-representative\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAB4AAACPCAYAAADtCmtEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAABkklEQVR4nO2cPXLCMBSE93lSUbkIdbgEV0rOlFyJS5A6KaholSJNMrHlQU/SKuzbGQpGA5/lXZ5+sG0pJTA0UagBlgA/5Bofd7t0mOfiLz9fLvi8Xu1m8GGecXp+KQYf315X27JgAMDi8fqVB5t9v0qV+ex2jxt1eRs8tQn+oB4DPo+LwdRwGc1jVqq1wgWDr8ul4TIAkyfV602jnmommDJImPkGCV/lKufmNLDHguGiDYvhcS+wx+NiMHfOpeex1EoiwtUVHAWkksLjv9LbEaB5rDfn0qtcUUDuHBzbTfU0MJjzO3ZWLle4XCVzvSkq1xKZAI7tpq5gsVkms3LRhkUxj1lgzcoV4eoCpoYrPO4G5swymatF2oRezGMobkVEuLqBYyuiDjc8XuLSwE24I/81QLtGQCxcsRVRT6N6zALHarEvuAl31B5TrzYWDFcUkDoa9FQDagUEBt/9Tv9y0SbmMVgrCc3tpnJuTr5wJfw+sJ/vN54SYLnnCJjZB4D3zYNb11NKaX8zuKX07qoPcDd9AfnAHrsHy6rzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 36x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# color bar, left\n",
    "grad = np.linspace(0, 1, 100).reshape((-1,1))\n",
    "null = np.zeros((100, 1))\n",
    "img  = (plt.get_cmap('bwr')(grad/2 + 0.5) + plt.get_cmap('bwr_r')(null/2 + 0.5))/2\n",
    "fig  = plt.figure(figsize=(0.5, 2))\n",
    "plt.imshow(img, aspect=.08, origin='lower')\n",
    "plt.tick_params(bottom=False, labelbottom=False, left=False, labelleft=False)\n",
    "plt.subplots_adjust(bottom=.05, top=.95, left=.05, right=.95)\n",
    "savefig(fig, f\"2d-hist-representative_cb-left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAB4AAACPCAYAAADtCmtEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAABo0lEQVR4nO2cMVLDMBBF/2aoUrmAmlyCK8GZ4Eq5BNRQpEorCjcwwDrRWv7L/P0zKTIe50na75W0sW2tNTC0o1ALLAG+8Q7u97dtmg7dP346veJ8/rCrwdN0wNPjsRv8/PLw5zEXDAD4tb1xuWCz+dMr79zFHg/q8DJ4N8j3OWMMxGLcDaaaayc11DOZBOaYC7EOe+f6PTbAIgkklDKlYgyQJgmz2CQRylxa1zEXPIZ7wbTIMlfFeDswY5XJXdDLxVhrJ1Hm2hZcCWQlVYx/gvUSSD/XVc6h5oL1Mlc/11XOoQbUXF3lpjWVGDyGe4G5IpmrF4xgyoyZSy9zjeFWuSkTeAw3aY+5mYs2LQa4nnLGmAYWzVxjuJl7rGeufq6rnEPNBVNKEWDuFmkLeqkYA5KliDLXZuAx3MwLer1SxBhu4stJ749r2j0CWuaqUsSayhljGrh2i5uCtW6Bpd5trGeuSiBrKedQz2QSmLb0iTzv9D83bVoxBkg7Cc1yU05zNXxv2NfvC28JMO89Amb2DuBtqXGO7ltrd1eDR0rvqfoCb6ZP+wQeu1BpzccAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 36x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# color bar, right\n",
    "grad = np.linspace(0, 1, 100).reshape((-1,1))\n",
    "null = np.zeros((100, 1))\n",
    "img  = (plt.get_cmap('bwr_r')(grad/2 + 0.5) + plt.get_cmap('bwr')(null/2 + 0.5))/2\n",
    "fig  = plt.figure(figsize=(0.5, 2))\n",
    "plt.imshow(img, aspect=.08, origin='lower')\n",
    "plt.tick_params(bottom=False, labelbottom=False, left=False, labelleft=False)\n",
    "plt.subplots_adjust(bottom=.05, top=.95, left=.05, right=.95)\n",
    "savefig(fig, f\"2d-hist-representative_cb-right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple sessions: fidelity testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_design_matrix(bin_whisker, bin_eye, shuffle_whisker=False, shuffle_eye=False, random_state=539167):\n",
    "    \n",
    "    def _shuffled(x):\n",
    "        np.random.seed(random_state)\n",
    "        x = x.copy()\n",
    "        np.random.shuffle(x)\n",
    "        return x\n",
    "        \n",
    "    x_whisker = np.concatenate([bin_whisker['left'], bin_whisker['right']])\n",
    "    x_eye     = np.concatenate([bin_eye['left'], bin_whisker['right']])\n",
    "    \n",
    "    if shuffle_whisker == True:\n",
    "        x_whisker = _shuffled(x_whisker)\n",
    "    if shuffle_eye == True:\n",
    "        x_eye = _shuffled(x_eye)\n",
    "    \n",
    "    X = np.stack([x_whisker, x_eye], axis=-1)\n",
    "    y = np.concatenate([np.ones(bin_whisker['left'].size, dtype=bool),\n",
    "                        np.zeros(bin_whisker['right'].size, dtype=bool)])\n",
    "    return X, y\n",
    "\n",
    "def perceptron_fidelity(X, y, random_state=539167):\n",
    "    y_hat  = Perceptron(random_state=random_state).fit(X, y).predict(X)\n",
    "    return np.count_nonzero(y == y_hat) / y.size\n",
    "\n",
    "def binned_fidelity_for_session(trials, nbin=5, random_state=539167):\n",
    "    points_binned = compute_binned_points(trials, nbin=nbin)\n",
    "\n",
    "    raw_fidelity = {}\n",
    "    for shuffle_whisker in (False, True):\n",
    "        for shuffle_eye in (False, True):\n",
    "            binned = []\n",
    "            for ibin in range(nbin):\n",
    "                X, y = create_design_matrix(points_binned['whisker'][ibin], points_binned['eye'][ibin],\n",
    "                                            shuffle_whisker=shuffle_whisker, shuffle_eye=shuffle_eye,\n",
    "                                            random_state=random_state)\n",
    "                binned.append(perceptron_fidelity(X, y, random_state=random_state))\n",
    "            raw_fidelity[shuffle_whisker, shuffle_eye] = np.array(binned)\n",
    "\n",
    "    fidelity = {\n",
    "        'full': raw_fidelity[False, False],\n",
    "        'null': raw_fidelity[True, True],\n",
    "        'whisker_unique': raw_fidelity[False, False] - raw_fidelity[True, False],\n",
    "        'eye_unique':     raw_fidelity[False, False] - raw_fidelity[False, True],\n",
    "    }\n",
    "    return fidelity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***in MLA007518, session2019-10-04-001: no epoch found from 5 trials (eye_position_from_trial, AtEnd <Backward> Left)\n",
      "***in MLA005757, session2019-07-18-001: no epoch found from 5 trials (eye_position_from_trial, AtEnd <Backward> Left)\n",
      "***in MLA007519, session2019-10-01-001: no epoch found from 5 trials (eye_position_from_trial, AtEnd <Backward> Left)\n",
      "***in MLA007518, session2019-10-07-001: no epoch found from 4 trials (setpoint_asymmetry_from_trial, AtEnd <Backward> Right)\n",
      "***in MLA007519, session2019-09-30-001: no epoch found from 3 trials (setpoint_asymmetry_from_trial, AtEnd <Backward> Right)\n",
      "***in MLA007519, session2019-10-02-001: no epoch found from 4 trials (eye_position_from_trial, AtEnd <Backward> Left)\n",
      "***in MLA007519, session2019-09-29-001: no epoch found from 3 trials (eye_position_from_trial, AtEnd <Backward> Left)\n",
      "***in MLA005757, session2019-07-16-001: no epoch found from 5 trials (setpoint_asymmetry_from_trial, AtEnd <Backward> Right)\n"
     ]
    }
   ],
   "source": [
    "nbin = 5\n",
    "data = []\n",
    "for subject, session in set((trial.subject, trial.session) for trial in trials):\n",
    "    try:\n",
    "        in_session = [trial for trial in trials if (trial.subject == subject and trial.session == session)]\n",
    "        data.append(binned_fidelity_for_session(in_session, nbin=nbin))\n",
    "    except NoEpochFoundError as e:\n",
    "        print(f\"***in {subject}, {session}: {e}\")\n",
    "\n",
    "data = dict((key, np.stack([item[key] for item in data], axis=-1)) for key in ('full', 'null', 'whisker_unique', 'eye_unique'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN0AAAC3CAYAAACWq1smAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5gUlEQVR4nO29e5gjV3ng/XtVukvdavVlem6eGc94ZnxLAK8NSfAHxHEcB0ySjRNuCSH7xbBZNskmhPBkvycbICRPyG5uu7k6JFy/AIHAAvHyBduAAxhI7BiIr+PxjOfmmemrulstqVSlqvf746hqpL6qu9XqmXH9nkdPS0el0lut89Y5570dUVUiIiJ6R2yrBYiIeL4RKV1ERI+JlC4iosdEShcR0WMipYuI6DGR0kVE9JhI6S4SRERF5Kpl3vspEbm3g3N8UER+u/vSRXSTSOk2ARH5ryLy+QVtR5dpe91q51PVv1XV27otZy8QkRMiUhOR+ZbHn261XFtJpHSbw1eAl4qIBSAi24EEcMOCtquax15UBDKu8P7PisgH13DKV6tqvuXxCxuT8NImUrrN4SGMkr2w+fplwJeBIwvajqnq2ZbP3doc/Uoi8mciIhB28q81n4uI/JGIjIvIrIj8m4hcv1AAEekTkS+LyP9qfuZqEblPRKZF5IiIvKbl2A+KyF+IyOdFpAJ8f3f/HYsRkVRTlu9qadvWHBVHmq/vEJFvi8iMiHxdRL57s+XqBZHSbQKq6gD/jFEsmn+/CnxtQdvCUe4O4CbgBcBrgB9a4vS3NT97CBgAXgtMtR4gIkPAF4EHVfWXgCxwH/BRYBvweuDPReS6lo+9AfgdoK8p56aiqnXg48BPtzS/HrhfVSdE5Abg/cB/BIaAu4HPiUhqs2XbbCKl2zz+iQsK9n9hlO6rC9r+acFn3quqM6p6CjMyvnCJ87oYxbgaEFV9UlXPtby/s3neT6rqbzTb7gBOqOoHVLWhqo8AnwJ+ouVzn1XVB1XVV1V7Hde7Ep9pjlbB483N9g8BbxCRoB++EfhI8/mbgbtV9Z9V1VPVDwF14Hu6LFvPiW+1AJcxXwH+s4gUgRFVPSoiY8CHmm3Xs3ikO9/yvArkF55UVb/UNET8GbBHRP438HZVnWse8ipgHvjLlo/tBV4iIjMtbXEudHCA0ytdjIj8OWY0BEgCcRH5sebrU6q60tTvx1T1/iWu5Z+b09mXi8g5zBr3cy0yv0lEfrHlI0nMTeWSJhrpNo9vAAXgLcCDAE3FONtsO6uqz67nxKr6v1T13wHXYaaZv9by9vuAfwQ+LyK5Zttp4J9UdaDlkVfV/9R62lW+863BZ4G3Ah9tOddG1lofwkwx3wj8fcsoexr4nQUyZ1X1Yxv4rouCSOk2CVWtAQ8Db8NMKwO+1mxbl9VSRG4SkZeISAKoADbgLTjsFzBGm3tEJAPcAxwSkTeKSKL5uElErlmPDF3mI8C/xyjeh1va3wf8fPNaRURyIvIqEenbEim7SKR0m8s/YQwXrYaJrzbb1usq6Md0yBJwEmNE+f3WA9QkSb4FM1p8FrMOvA14HWakPQ/8HtAro8Q/LPDT/e8WWc8Aj2BG2q+2tD+MWdf9KeZanwF+tkfybioSJbFGbDUi8n7MdPs3Vj34MiAypERsKSKyD/hx4EVbLErPiKaXEVuGiLwHeAz4H+s1Kl2KRNPLiIgeE410ERE9JlK6iIgec0kr3e23364YU3P0iB4Xw6MjLmmlm5yc3GoRIiLWzCWtdBERlyKR0kVE9JhI6SIiekykdBERPaYjpRORHxGRKGQsIqILdDrSvQc4JyJ/KiIv2SxhROT9zdofj23Wd0REbDUdKZ2qvgC4FagBn2oWtvmNZrBqN/kgcHuXzxlxEaO+Mn3vNKf++ymm751G/Y7dXZesDB1PGVX1O8B3ROQdwA8AfwC8W0QexBSN+Ziq+hsRRlW/sgmKHHGRor7yxBueYPrz02Hb4CsHufaj1yIxuWxlWNM6TUQOYDJ8fxrwgd8ETmEyle/EpGhcMjjjDl55YdJ1b5GUkN6d3pLvVlXskzb2Cdv8mh0cjzb/+rQ914aa5x6op2a08Eynxm9vwzft5UfKTH5qEvU0jOeY+LsJvvLpryDW2ju8LhUUstyg1WxXX6HRbLPAylpMf36a0v0lBm8bXLMMndCR0onIf8bUsLgK+ATwRlX9Zsv7nwLGN0XCxbK8BZMVzZ49ezZ0rsRQgsrRCl5pixRPwS25uFMuqZ0pENqCisLOTbNz6IWOj3JBUdYwG1JfcSYc3HEXfLAKFiLS+eo+hjleLrwmBhIT025x4W/M/EVALEFigopSfbpK+Rtl5v91Hm2oOVdLedvEjgTJnUlzvWhHNwRadXTB/1Aw3ysxCf/H2jD/NHfCxT3nmlO0KPr8d+a3VumAH8ZMJz/brOnYhqpWRaQno5yq/hXwVwA33njjhibfYgnFlxa7IlcnqCrenIc75eI7PiKClbfwGz72szbpK9PE83HTMYLOHcN02maHkdgSbbLyqOA7PrVnatinbFAo3FwgezhLLB6j8mSFzIEMseTmeY9UlflvzzP2/44x8YkJ3GmXeH+cwi0FZr8yaxQ2dkGRD/3JIYq3FY2iBsrqK37dx7f98K86uviGE4NYKkYsHSOWiiEpwa/7NKYb+HWjvVbeIjGYMKPavdM8/hOPL5I5/4JFhdi6RqdK94CqfnJho4i8TVX/EEBVV93g4vnGkkrWb5HanSKWjOHZHo2pBuk9adJXpKk+XUVzSnIoueHv9ioetWdq1M/VQSBzZYbirUVi8QvKVT9XJzGU2DSFqz9XZ/zvxhn7yBjVo1UkLhRuLrD9Z7eTPpDGOe8QS8aYfWDWfECh8IoCmYMZ6ifr7SdrUSYra5EoJpCkLLrhqK80Zhs0phu40675n/dZpHamiKUWX2fx1iKDrxxctKYr3rp5N+OOklhFZE5V+5don1bVro3BIvIx4BXAMDAGvFNV/2a542+88UZ9+OGHu/X1G0ZV8coe7mS7ki3s2I35BvUzdWKJGPFinPrZuhl5ErFwNErvXds6T1VpzDSwT9g4Yw6SEDJXZkjtSbUpWoBne9gnbHJX55Y42/rxKh6T/zDJ2IfHmP3qLKpK7vocwz86zPY3bScxnGDmqzPEB+Lkrs1RPV6l/PUy9mmbvhf1MXj7ILF0bNXRO7xuz1y3O+2inplCxgtx4oNxYonObibqK6X7S8x/Z578C/IUby2u14jS0YdWVDoRuaX59B5MEdPWk+4H/puq7l2PdN1gq5VuWSUbTCx5V3VLLs45BytnRrtgDeE7PtUj1XB66U65OOMO2cPZFX989RV30sU+beNOuFhZi9TelLmrr9Lh5h+bJ3dNbl0Gi6XkmP3qLOc/cp6pz07h1TxSu1IMvnKQ4VcP0/99/VhpC/uMTeWJCrlrc6R2pZh/dB6xhPTe5rS6A3zXp1Fq0JhpoL4ilhAvxs3I14Vr2SAdCbDalQajTApTVz5AMWXcfnHRJy5jllOy1K6lpy4BzriDO+ESL8bJXpNddBePJWPkrs9Re7qGX/RJbksSy8aoPF4hczCDlb5gZfBdH2fMGEKcKYd4Pk5yV5L8d+c7vrPbJ+02pV8v1aerjP3tGOMfG6d+to6Vtxj6sSGKtxTpe3Ef6b1pYvEYXtWj/GgZv+LTf1M/iaIxYKmrZK7MrKhwft3HnXbx5jxUNZwdZK7K9Myt0G1WVDpVvRJARD6sqj/TG5EuHhYqGUC8P76qkoG5+9fP1vHmPBIjCXLXrTyNExGyh7PYZ2xqx2tk9mfM9OupKlbRQm0zjWrMNLD6LZKjSXLfletY0QIa8w20oSQGEmv6XIA75TL+9+OMfXiM+e+YkWrgFQPs+X/2kP/uPImRBKldKWP8UKX2bM3coBBy1+eI5+M4Ew7Ocw59N/QR72vvgl7NrIG9eWNRjqViJAYTJLcnO55yXux0NKY/XxQuVLIpF7/uX1iEd6Bk4Tk8xT5t49d8UjtTa/bBpXencUoOs9+YJV6M4zs+81+ZJ16Mk7smR+763LoNH+or9rM2uevXto7z6z7T905z/sPnKd1XQhtK7roc+9+7n8LLCsSsGPGhOMnRC4rRmGtgn7JJ7kziPeeFI7Zf9yk/XKb/e/uJ95vu55ZcnDFjFLcyFvGh+Jb5LnvBskonIk+q6jXN56dZxhukqhtzlm0hyyrZMpaulfAdH/uUjTaU9BVprNyK+ypeWLx/e578C/MM3DIQWtz8mo/netSfqJO7NsfIj4+Eo5wk1n+3D0bQTkYMVaX8r2Vj5v/kBI3ZBsmRJDvfupPR140SH4jjzXsktydJDF4YNdVTasdrxJIx0vvT2MdsslcbI5GqMn3/NH0v7iMxkMCrGWNOvD/edYPOxcxKI92bW57/9LJHXcI45x1Q1qVkAV7Vo366DjFI70l3dJ620KOmA7dwc4F979mHlTOGmNy1OSQhVJ+u4pW9tnVe9lB2zaOdO+0SS8babgYLFb94a9GY+T9mzPy1Z2vEUjGGXjXE6M+MUnhpAeecg+/4xItx0nvaRyNn0sEdc0nvT6MNpf5sndx1uXDtNfPlGXLX5UgUE1SfqZop9aHsxWAA6SmXdN3LrbReNuYa1J+rE0vHSF2xtFkemuFPnpowqYZ5PvbRMY796jETFRFEWwhkD2dJ7kiGjt1Y2jx81yeWjJEYSSAJMUo4miQ+EG9zBAfHS0ra2oiDc9oh/8K8cRgnBZS2mEN1FWvAojHXQEQofG+B0Z8ZZfjHhpGEUD9l/GapPak2ww4Y407tWM2sd3emcKdbrK/NUXXuX+dAILktiVf2SF+ZXnSey4CNWS9F5Lc6OYGq/manEl3MqBqFwCNUjlZFCZ67k6ZDWTljzFBHsY/by4diWcY4MP/wPHPfnGP2a7PUnqmZaIoYSPzC7yRJId4fNxEUMw0TfRE8aj5e1TOyOCaOsS0ca6Vr83WxpS8G3lwz/K0puz/uM/qmUa5815Wk96Tb/InpK9NLGm3qZ+s0ZhtkrsoQS8Son6/jV/226WLl6QrOeYf03jRW3rqs12udsNL08oqeSbFF1M/VL3S8ID4w3owbjDdDkBKCpAV3wjUjzK4kfTf2rbgu8hs+84/MU/piiZkvzVB+uIx6ipWzKNxcIP+iPJOfmVykCAf++4EV4/28mokySR9IIyLYJ2zcGZfUjhRa17YQqeCvM+bg2z6xdKztmNIXS8x+dfbCtTevN3dNDqvPovJEBStnkT249PTPq3nYx20Sowly1xgFs0/ZEIPM/kx4XOXJCuWHyhReViCzL7PoPM9HllU6Vf0PvRRkK0jtSMGO5d9XX6mfqYcGg+yB7LLH1k7UmPnSDKUvlZh5YMYos0D+RXl2v203xVuKZK/LUj9TJ3MogzZ0zaFHVsYyboQjVZI7kuSuydGYNaNR5nBm0RTXr/vUnq0tMlI4Ew6J4QTlh8oLLhjihTi+7S/pTwQzI6ifquM7PtmrLyhk7VgNK29GfzA3ntkHZ2lMNxh85SDJ4Y2Htl0udJza09xA8CeAUVX9BRE5DKRU9d82Tbotwnd96qdNx0rtSi0yGIBZ0818ZcYo2hdL2M+aDURTu1OM/PsRBm4ZYODlA6FlrzFr1oC5a41h4dqPXruu0COxhNy1OWrP1vAqHundaWKZGNUnq2T2Z9oMJbVnamQOXxhdGvMN6ifrJIYT7LhrB6UvldqMOcUfLLLjrh3LytEoN7BP2qSvSBMvmK6jquYmMJokUUwYpTxTp3q0SmI4QfqKdKRwC+g09vIngT/HbA7/BlXtF5EbMRvT37rJMi5Ltw0pnt20RCqkrkhhZVoiQZaZMsayMQZeNkDxB4oM3DJgIiUWjBDOuENjtkH24PIj5Xpwxh0apQaZg0axqk9XSQwlSI4ksU/ZYTS975gsBkkJ6T3pUKnUV85/4Dz2GZvC9xaWVXz1jRtALCG9Lx1en3pK9akqqb2pMHwtCLC2+iy8Oc9YLy8Tp3YHbDz2MjxI5Eng9ar6bREpqWqxuf3uWVUd2aCg66ZbStdqMEjtuRC3uNKUsfgDRRPudFPfiuZ7+3QzgHmJ0bIbeBWP2vGaCZhOxrBP2zTKDWKJGJkDGeyTNuoo6f2LDSHVI1USo4kVo1PcaZf62boZRbPt4WjVp6pG4X2Mv20gHo6a3pxnUobW6Yq5ROlK7GXANuA7zefa8veS9Teor0x+dpLZr82Sf2Geba/bhl/zmf7C9PJTxu8fYOAVA23O4JVYuM7ZDKycRfaaLLUjNVJXpEjtTmE/aJz0vu1fyNFbwGoKF7gBrLxF/vr23LLQoHNVmvrpOmKZEDZn3DEWYIHEtqWDviM6V7p/xWSOt27E/jrgX7ouUQ9odU4HfrRjv3YM3/bBJ5wy7nrrrmWnjKudv3qkaqI1iuuLcVwLsXiM7LVZasdq2P9imyDsPSm8Wa/NJRGwmsLVz9dpTDeWHKmCdZ3VZ1E/WSe9z/jb6ufrqKvEB+PGPfA8dwusRKdK90vAvSLyc0BORL4AHMJsHn/JUbq/xNQ9U/hVPxyr3UmX0Z8aZcebd9D/4v51xzf6jea068rMqqFg9hmbxkyDeCFuov43sPYJo/GrHqk9KdJXpGEXoaUzUP7KUxVSO1KhIaQVz266AYZNRMxC3GmX6tNVYrkYiWIiPIczbtwS6SvSJnVnjbGdzzc6DXh+SkSuBu7A5NadBu5R1fnNFG6zmP/2vCkR0PRNBeUP0gfSDNw8sO7z+nWf6tNVE6a1wtQqGC2CgOjGbIPqE1Xig3HjxlgD6in2CTuMWhl+9TBe2TPhYoezFyyd8x5exVtW4exTNr7tkzm02PUAJm6z+nSV/AvybTI6Ew5exSNzZYbKUyYV6XlkOFkXaynBV8UUJbrkyb/QrFFi6fbOlT2YpfJUBREhMZIgXox33IEa8yZrO3ft8omhQTCwxCV0HYDxjcULcZxJh/nH5kluT3ZkZq8/Z6JB0vvM2ip7tfGtxfvjZA9nTWLsvjSZKzPMfH3GxF72tY++XtUYYpZzjfiu8bdJQhi8bbDNuulMmmpqmf0ZU/qhmLgcQ7u6zkphYF+lA0OJqr6sqxL1gOXqYmx77bawCI474VI9UgUgMZwgMZRYVgGDWMOVzOPOmIM7aYKBW10RrSSHjbLVz9epPF4huSu55LrLLbnUn6uT2pkitytnwtL6rLYOH0vEyF2Xo/p0lfmJefLflUfiEo6AgaXTt/22G0CAqlI/XafypMn0Tl/RrpDutEtjpkH2qixe1cMre2QPddclcrmyrMtARN7U8vIA8H8DHwJOAnuANwHvV9V3braQy7ERl0GndTHUV9wpUyYPIDGYIDGcCI8NYg1bQ59aCcOlhhNrsmKqKvXnTJhaao/xg7WmwiR3mtw13zFWxiAUa+E5qkeqoCYZNLM/g9/wqTxaMRElB7NLWmKdSQfnvFmnZQ5kFhmD3JKLO+mSPZhFfaXyeKUtm+B5TFf9dN8Efk5VH29puxajdN+zbhE3SK+zDFTV5LxNumEuXnI0ueS0TH2z1lJPTQ7bOtNXAsd07RnjEshd3T59rTxeMeuwBT64QOGCNZw74+I85xAvxmnMNkDMtLZ1feZVjFJbAxbejGcskwuMQe6MizvuhqNa9WiV1M7Uqkaj5wldVbpZTPiX3dKWAc6p6sB6JdwoW5Hao2oqGVeeNmu/IDsh3h8nsc0Ux3GmHepn6yZSIxcn+B+rKp7nUa/XcRyHRqOxyrc1q1BXTT6dV/JQX03pgoTgjrvE0rFFhhFVNeFeIwmsnCkm63s+teM1vFmP/PV5YukYzqSD1pTEjgTuWZdkOknuihzeCY/MVZlF67NGuYFzzgkVzplw0Iau2fizEgun562vu/XeJtJVpfscUAX+G3AGk4HwLqBPVV+9fhk3xkaVrlwuU6/XVz+wFQX7uE1yR5J4/oKhxSt7OOcd6mdMZ89cncFTD8dx8DyTySAixGIxUqkUyWSSRGJ5H54zbYoPJXckw7IG0MxQP22bTHdLFq2jVJXaMzWS25Nh/RG3ZEan1L4UEhdqx2tmmjyYMIo05pDanaJer1M+VsbabYUWzHg8TiqVIubEcMdcMgeMddKzPepn6mSv6u46bmF/XOn1et9bjoGBAeLxDe0I11WlG8TEXv44pgB2A/g08IuqOrkBITdEr0e6MPTpqkxoDPF933TWZ8s4FYfkFUm0bqahCRJkChlyuzuva+JVPOyTNvHi8u4DVaX8cJlYzhReDeIpwynlzhTx/riZmh4z2d8Lp8D2KRv1NUy3aZQb1E/VTXZBy9rMdV0q0xXmT86TujIV3jgazzYY/K5B4slo28IWuqd04cEiMWAEmNjoDj3doJdK16g1mH1yFmufhadeeOf0qz46ruT35skOL06H8SpmBPQdHytrkdyeXNKH57vNoOR4M6h4BaNE7XiNxHCCeH/cxI2eqhPLx0zWwa60aS8bF8bCzINW3CmXuTNz1FI1KGNGwgXy+zWf+tk66f0XAp2rJ6r4fT5e3MP3TTcQERKJxKoj+MVMr0a6lVwG+1T1RPP5/gVv58JIc9Xj6xRwyymXy+RyOWKxC0qgqjiOEz5UFW/ewz3v0n9NP+lMmkQiYZzSx20kKaT/XXrZdYOVs8gcMKOJV/WoP1fHr5uk0uR2U5qhfqqOVzMO5tXiFRtzTSNIc8oZz8exrrGYe2gOFLwBD3fabA6Su37lCH8pCG7VZZhh0lcvNgZ5VQ97zGbkxSPhedwpl/xofpELQVWp1+vU63Vc1w3bE4kE6XSaZHL5EnrVapVqtbridW82tVqNfD6/UaXriJW+4VGgr/n8GYzPbuF/TWnbb+XSIpFIcPLkSfL5fNghgjt2KpUin8/jTrk06g0yL7kQaVE/V6dRahif2xqcwVbWCl0Lnt0cAWs+qd0p0n2rxyqqr9in7LYamqomvSZ3dQ5JCnMPzRGLx1aNDFFVpqam2LZrW9tNJ8CredjP2m3JrL5jMtGXChETEdLpNOn0An+e62LbNvPz8+HswLIsUqkUlmUxNzdHOp1meHh41evfDFSVUqlEoVBYJPtmsZLSheXSVfWyDRcfGRmhXC4zOjq6qPPVn6vjuz7Zw8ZY4FVNR0xsWzo2cS1YaWvN5Qtqz9RCQwZcULjUrhRe2cOb9xj4vgHEklUd7FNTUwwODi6tcLZnfH8tTnNVs8VV9uq1GU4SicSi6Waj0WB8fBzbtsnn89TrdWKxGLlcb2M2G40GU1NTDAwMkEp1zwK7Gisp3QmgH0BE7t/KZNXNIrizpdNpJiYmGBkZCTthUH4usy8TFmlV1UWGhl7hTJpiSIEBJ1C4xLZEaDENbg4Aqe0pkqNJ6s/Vcc46oYMdYG5ujkwms+Tay7NN2s7CKBX7hCnFvlzVs06pVCpUKhWGh4dJJi8EC0xPTxOLxchkelNHpVqtUqlU2n7zXrGS0lVF5HrgSeDFYm6vi3rbxWBQ2SjxeJyhoSEmJycZHh7GPmqHoV9BdEYnWQObhe/6uONuOLoGCicpMWb8g4ud49Cc8u1Oh9PS+pk6jIKHR3//ok2YTE2Vo4sVzp1xTY7cOkuxg5lmlkolstks27ZtW/T+4OAgk5OToUtlswimk5ZlMTKyNfnXK6n4uzH5cg6Qw7gJ3JZH8PqyIB6PM9A/wKlvniKxI4GVt6g8WUEdJX99fksjLmpHa2FJBlUTduVXfayUKVS02n4GEhMy+zIk9iWYPj5NaiIV7s0Q4DvNDIlr2qt/+a6P85yz7kpevu8zPT1NuVxmZGSEfH75zRaHhoaYnZ1tM8R0E8/zmJiYIJfLUSgUNuU7OmGlamB/ISLvA7YDTwHX9UyqHuE4TjjF8R0f56jD9uu3c/74eQYzg2YE2eB0aqPUnzNTx6AsefmhMqpK/gX5NRlxVJXSbIldL9oFDTNdJAbpfWlQk3eXvSa76HprT19Q+LUyPz9PrVZjYGCgIzeCiDAyMsL4+DjDw8NYVvdudLVajfn5eYaHh3s+nVzIarv2NIAzIvIiVT3ZI5l6RqPRYH5+nkK6YNJbdqaon6ozvHuYsl8ma21t1LxXM8aR7OEsvu8zc/8MyZ3JdRX7mZqaolgsmg6XhOyhbFh2QR1dUuHsUybyZq0JvY7jMDMzQy6XW/MUTkQYHh5mcnKyK+stVWVmZoZYLLZl08mFdJrEenSzBdkKstksjbkGZ586y2BxEK9yoXqV5VjhGm8rkjJVNcweaMw3KH2pRN+NfaR3rt2sHRhOWg0X0KyjuczGHY25Bupqx/VgwEwlS6VS2MHX+3+zLCtcY2/kPJ7nMTU1RX9/f8/cAZ1w2boCOsEZd/CP+qQzaexB24RTNX/gZDJJoVBgcnKyo7i9bhPUl7TP2Mx+dZaBlw+sS+Fs28bzvDWZ49Uzhpf0/s6/b25ujqmpKQqFAsViccM3qng8zsDAAJOT64syrNVqTE1NMTQ0dFEpHDzPlc4re2YrqheNkM6nmZmZaXs/mUzS39/fc8VrlBv4NRPY7DznULi5QKKwdsuh53nMzc1RLK5t0/rq01VTTr0DxanX64yPj5NIJBgZGelqREcymaSvr4+pqak1fa5UKlGv19m2bVtX14Xd4nmtdMkdSeKDppPkcjni8fgixUulUvT39zM1NdUTxVNfmXt4Dl991Fdy35VbtFtpR+dpRpysFunheV5btEj9OVMBerVwNM/zmJycpFarMTIysmn+tSDKZeHvspxM4+PjZDIZBgYGNkWebtDRr9nMMng78EKgzeZ7KZZrCPDrJsg42Dc8n89TLpeZnZ1tMymnUqmOO/GG5Gn4lO4rkTmcQW0ltT+14n7cKzE9Pc3AwMCKhohqtcr8/Dz5fN7cVGwlNZ+i/+rFPrwAVWVubg7XdSkWiz0ZSXK5HOVymXK5TF9f35LH2LbN3NwcQ0NDF+Xo1kqnv+hHgRSmMNHWRqZ2kUQxQaKYoDHXoHqkSiwTI39Fnvnq/CLFC9YFgXGl27jTrtmP4Kqmwu1Zv8IF8YwLDScBge8smUyGjupMOsPco3M0djdCX1Y22269DTp2f39/z/1cfX19zMzMUK1WF8k1MzODqi7pdL8Y6fRX/T5gRFXXmPF5aRDvjxPvj+NVPKpHq8QTcerFetjBAtLpdDjiDQ0NdeW7g5w3iQtWwULdjSmcbds0Go0lI07AGBjK5TKDg4Nt66/a0Rr5Q3msjGUiXqpVJiYmSCQSZLNZ5ubmwnXbVpXYGxgYYGpqilgsRjqdDq2TfX19PQsf6wad/rL/BuwGjm2iLFuOlTMmdM/20FNKuVrG2+VRHLlgiAh+3JUULyjL4Pv+kn/DXDzbp3aqRt/+PqxJi1gjRvqqpcugd0JgOFnKH6WqTE9PE4/HF40I9fN1rMKFuE4RCUe6qakpTp06RV9fH319K+/L1wuGhoaYmJjAdV1qtdq6p5NBCle9Xmd2dpYTJ07g+z67d+9m//79m3qdnf66XwL+UUQ+AJxvfUNV3991qbYYK22ROZghbscZe2KMueNzFK8sQppQcWq1GseOHVtymhVkV1uWhWVZxGIxEolE2CYiOBMO7rzLyPeMUD5dZurMFJlrMghCTtfu/G5dcy78bDAtLBaLbZEhjTlTqkEsWVTNLBgRC4UCw8PDJhqmXA6nrq3pUBtBVTl+/Djnz59n+/btHXX4RCLB2NgYe/fuXVXhVJVGoxHWpfF9H8dxsG0by7JIJBI8+OCDPPvss+FnDh48yJ133rlpitdpuYYvL/OWquot3RWpczZUgk+Vxx57jHPnzjE8PMwVV1yxqJhNoCTzc/P4Ez45K0dmd4ZU0WRYV6tV6vX6mk3ytWfNdDJ9RRqv5lH6YomBVwwQz8ep1WpUKhXAOO8zmc4qJgfTrNZ1XBDcG4vFKBQKiEioaNpoFlMaTbRFojQaDUqlUmi1XYogP05EKBQK63YTqCqf+tSnOHr0QuzFSh3e932mpqbI5XKkUinGx8cZHBwMb4SO41Cr1ajVariu2zbDaDQa4SxDREJlPX/+PN/+9reBC785wGte8xoOHDiw1kvqfrmGi431Kl3wYz995GnwAYF9A/u47eBty3bwOXsOC4vEXAKtmo0yrD6LmlvD9Vz60/0rl+bVCzu7xgfjxHIx6o061bkq/S/tJ9G/2A9Xq9WwbTvMNWtVqFY55+fn2/LRRMTUbWmOVDEnhjvumqplfXETy9lUtDBB1feZnZ1FVUMFXfg/W/g6cDc0Gg0ymUw49V6pKJDv+9i2Tblc5vjx43zrW98KK6wFFAoFEolEqDDBI6ieFrwOPrdUJbCFMiz3uwbnAojFYuEN5GUvexkvfelLl/zMCnR1qyxEpAi8GtgFPAf8g6qW1irVxcDx48c5evQonu8hIogIJ2ZPcHLwJAcPHiSTybTfvQUGGWRmZiY0LDjnHRrTDfIjedyci9twl/QNBT92o9Kg9FSJ5A2mfJ6IMJgZZFdm15KdRFXDEdT3/bBymWVZ9PX1hfLV63WSyWR4bBBr2Kg0KNpFmIdYPkbfgb4wG2Hhd1UqFarVKoVCYcmyCiu9Dkz41WqVWq1GPB4nk8mErpfZ2VlmZmYolUrh6yCLIFCqhecMRtBgKt5oNFBVMpkMlmWFFcqSySSxWAzHcchms+G5LMsKbwLB9D6Y6rdO+WOxGGfOnOGee+5Z9Ltt3759UVu36NRP973A/8FkG5zEbCTyxyLyKlX9xqZJt0mcP38+vEsGPxTAP977jzzy7UcYGBhgcHCQwcFBtm3bxrZt20gkEhQHi5RKJWp2jdyuHOwytR/llOCnfGaYaev8juNQrVaxJ20apQYDVw+Qy+eW9J2t1LFjsVh43iBI23VdLMvCcRxGR0cREWqlGuPHxulL9tFf7Cd5MLlk2k8gl+u6nD59mpmZGfbu3duxocTzPMrlMjMzM+GjVCqFihXUOwnOFY/HKRQK9PX1MTo6Sn9/P0NDQziOw3333dd2blXl5ptvZs+ePaEVdWBggP7+/jbZGo0G1Wq1rZ7N7t271zzVHRgY4MiRI4umuPv3LywL1D06XdP9M/BHqvrxlrbXAm9X1Zs6+iKR24H/iamp8teq+t4F778C+CwQrGg/raq/tdI51zu9PHbsGJ/4xIW9UAIFPHToEI7jMDU1xfz8hQ2JgjjAwcFBRkdHyefz7Nixg6GhofCz9ek6E8cmqFEjvyuPoqaA0aRiWRbJHSuXVBcR4vF4WN4gHl958xJV5ezZs1iehTvhUq/VyQ5kGT00ipVsNy4EShZMz5LJJOl0ms997nNLrqfAjH4zMzPMzs5SKpVCxZqZmWmLYAFzU+jr6wvjLovFYjhFDK4llUqRSqXIZrO4rku9XqfRaHDvvfdy8uTJ8FoDGYKk16GhIeJxU7A3WK+pmv9pNpsNE16DqXgn6+tWy2VQfOr06dNUKhV27ty5EetlV+teloCh1ixxEbGASVVd9Sqbxz4N/CCmWO1DmO2Un2g55hUYJb6jE8Fh42u6o0ePhp1n79693HbbhTVdvV6nVCoxPT0dPkqlErZth+dIJBL09/dTLBYZGhpieHiYvkQfyYoJlsaG9O40yaHlK2G1ytRoNHBdF9d1wylVQGABTSQSxNwYE89MkLWyxPNxKqkK6Vw6VKrgbt+qZMGU2bZtbNvm6NGj3H///W3rKVWlv78/VIhWgsTPgYEBisUiAwMD4aOvr29RRbX5+XnK5XJovKjX6+G0MZAluObjx48zMTHB8PAwe/bsoVKp4HkemUwG27bxfT8sfNRqWFo4JQ0+12pRDtaQgWElODaZTJJKpUKrMtD2fJ10dU13FLPz6kdb2n6Szv12LwaeCcr1icjHgR8FnljxU5uEiHDnnXdy7NgxxsbGljVV79q1q+11EG84NjZGqVTi3LlzlMtlTp8+zbFjx8Jzp5IpCskCo0OjbBvbxo5DO9i2Y1tbtPtSpvKlivgEuBWXyukKtWqN2cYs2qdMN6apTdbM1LJWC8vfBXfwQIHr9Xo4EgQK1rqeav2/pFIprrvuulCxCoVCOGoth+/7VKvVcCT0PC9U9GCkC5zZMzMzTE5OksvlwvKH+/btY//+/biuy9jYGLFYjGw2G65rg+8O1t/B84UVndPpdDhLCaJWLMsKA6cXKlRgDApuDL0ovwedK90vA/eIyC9h1nT7gIOYtV0n7MJsJBlwBnjJEsd9r4h8BziLGfUeX+KYriAi9PX1hYv1M2fOhCbj4AGEtRxVNYyE2LFjBzt37uS6665jZmaGVCoV5pKNnR5j/Nw4M/UZHjv1GP4J34zrGBdAcbDI4OAgZ8+eZWpqKuxEO3bs4IYbbghL1tXrdWqVGrUZoyyO5+Cog123wynWwusJ2uLxeFj0NR6Phw7xQqFAf38/uVyOmZkZvvnNby660dx0003s3r07tOrNz5uQuKA8fKDMjuMAhKNzoCjB/zT43zmOE35HYB2Mx41rpFQqkUwmyefzoSV0dHQ0VLjgEUwBg0fw/YHfLbj+RCJBJpMJ15TBTS4YaVfqC0sZszaN1otZ6oEZMvcD24CfBt7R/Du42mdbzvGTmHVc8PqNwJ8sOKYfyDefvxI4usy53gI8DDxcKBQUY6hXQB9++GF9+OGH29re+c53qqrqjh07wrYbbrhBVVXvuuuutmOffvpp/fCHP9zW9t73vlfHx8fb2m677TadmJjQ2267ra194okJ/cPf+sO2tve973163333tbXdeOON+q53vatNpnw+r+95z3v0lltuaTv2V375V/Ttb397W9uP/MiP6Cc+8QkdHBwM2w4fPqyPPPKI3nnnnW3HfvnLX9Y/+ZM/aWt729vepl//+tfb2g4fPqzve9/79OUvf3lb+5NPPqnvfve729o+/vGP66OPPtrWdtddd6nv+3rDDTeEbTt27FBV1Xe+852r/k7veMc71Pf9Rb+Tbdv6pje9qe3YJ554Qj/5yU+2td19992qpoOEj1e+8pWqqnrHHXe0tauq3n333W1tn/nMZ/TMmTNtbW9+85tVVTu+pk71odM1XQWzWci6Kn81rZ/vUtUfar7+r5ir/90VPnMCuFFX2Ctho2XVgztgrVYLN/kITM2dzu1VFa/hcfJfTlLYXSA7kg1HCc/zmJ2dDUcBALfu8o0HvsETx59AFzj29m7by1V7ryI7kiWWuGDmDoKXz507x44dO+jv7w9H49bRuXXUWu53DUaKcrkc1p8slUoUi0W2bdvW5rcKRspYzNRnsW0zyooI2Wx22YDqgFb5Fsq58HmwJ0S9Xl927bUW44aqMjExsSgKZ5Pp6pruW8AhjMtgPTwEHBSRKzE+vtcBb2g9QES2A2OqqiLyYkyu39qyF9dIpVIJjQfrTQfxbR/7GZu9L9rLTHUmtKqBmXqlUikmJibC6PhEKsGB6w9w5NSR8D6pKAhcee2VjGwfCddAiUQCy7Ko1+tMTU1x6NChZVNb1kLgSwNT+s627dDvtfD/4DgO5XI5XF/t2rWr484fOLBb/7Y6ulvbgtJ7q6UjdcpmFjnaKJ0q3QOY2MsPYtZm4W1UO4i9VNWGiPwC8AWMy+D9qvq4iPx88/2/BH4C+E8i0gBqwOt0kyfZK5WD6wR3ymx7nL46jeM6WJYVlmlPJpPh2mr37t2Uy2USiQS5XI6hoSGOHDnCM888A2I6yIEDB7j55pvbOnRQMs5xHIaHh2k0GovKFwRrmdb121pGhKVKoQd+ONd1SSaTDAwMrKvTBsqzVR0+ULyFhYS3mudt7OV6aN3QsXKyAoLZhbRZmz8YKSYnJ8Pojlamp6dxHCeMqpiYmGBiYmJJ62kQJZLL5bBtm8HBwSVlCkaNwM2w0NUQ+P8CH2AwXVyIqlKpVKjVaqHfbbXp46VCo9Fgenq6F2lJUezlegmsdAv9VZZlkYgn8E56ZHdkSQwtvVYI1hMDAwMkk8kwDhJMpeNCobAoETOgNcE0l8sxNTW1oc4S+MJafYDL/eZrCbC+1HAch9nZ2c0uw9e9NZ2IfEtVX7RE+8OqeuNaJbtYCBbvgRk8IFhf5HK5tkW4Z3sm2fOqfJh7thQiwuDgICdPniSTyZDP59s265ienkZEFiVetiaYWpbFxMQEQ0NDG1KC1unncomeukTQ8OVGa5GjbiUgr5dO13RXLWwQ8yttXoBaD6jVaoC5w69m4XKnXZzzZpuo1rLjCwn2WhMRdu/ezezs7CJr6ODgYFjhKpPJoLo4wTQoZ9eL9dDlrnABrQ76rSxctKLSiciHm0+TLc8D9gGb5rzuBZ3Wggy2Cl5ue6xGo8Hc3FwYutQ6OgUL+YWm66GhIaampsJpbOv75XI5NJVHdJdOihxtNquNdMeWea7Ag8Anuy7RRYT6Zk+2xHCC5HC7UaHV8BCPx+nv718yjKjVgrawLsnQ0BDVarUtgj6Y7m71FOhyJihyVKlUer4nHqy+l8G7AUTkm6r6hd6IdHEQrN8yBzJY2QtTvFa/Vae1+qVZnz9Yo7UqXqtBJXCmXyw19y9ngiJHQfBBL+l0L4MviMgtwOuBnZjYyI+r6hc3U7itwi25OOcurN+CJNJgl59wI441ENT3X0rxwIycwSL/+bLG2mqCIkexWKyn7pGOeo6IvA34ODCNSWadAj4qIr+6ibJtCfZpm8ZMg9y1OWzHZnJyklKpRDqdZmRkxJQ/WKeTNVC8qampMOwsINj3+mKKnHg+MDw8bDLtF6QybSadOsefA35IVR9rabsOuE9Vd26ifCvSTT9dsH6LDcSwU6Z2ZDerXrXi+36YP2ZZVpgUulUL++c7C3+PDdDdGinAMwteHwcuXc96C57tMfXoFLpdSSQT9PctbRTpFrFYLNyDLUgajQwnW0fwe/RqWt9pz3oX8Dci8i5MLtwVwH8D3iki4VxrvVkIW4XrukyfmaY+VmfouiGy+c52qukGlmUxPDwcliSI2Fp6Oa3vdHrZqkxK+zAavFZV7emCZKPTy9lnZ7F8i/yBjQU+R0Q06er08soNCHLR0r+nf8XokoiIzaBTl8FJgOZUclRVz22qVD0iUriIraBTl8GAiHwUsGkaVETkR0TktzdTuIiIy5FOHU5/CcwCewGn2fYN4LWbIVRExOVMp2u6HwB2qqorIgqgqhMicmnswhcRcRHR6Ug3C7RtPyoie4DLYm0XEdFLOlW6vwY+JSLfD8Sa1b0+hJl2RkRErIFOp5e/hzGi/BmQAN4P3I3ZmyAiImINdOoyUOCPm4+IiIgN0KnL4NdF5KYFbS8WkXdsjlgREZcvna7p/guLN/t4ArPHQURExBroVOmSgLugzQF6m3IbEXEZ0KnS/Svw1gVtPw880l1xIiIufzq1Xv4KcJ+IvBFToOgqYBSzyWNERMQaWFXpmvUta5gNRO7A5NJ9GrhHVedX+mxERMRiVlW65i46j2K2yvr4asdHRESsTKdrumCrrIiIiA3Sk62yIiIiLtCp0r0UeBZ4+YJ2xYSERUREdEinYWDfv9mCREQ8X+i4aqqIDInIG0Xk15qvd4rI7s0TLSLi8qTT2MuXA0eAnwJ+s9l8EPiLTZIrIuKypdOR7o+B16rq7UBQf/qfgRdvhlAREZcznSrdvpbNQgLLpcPaKkRHRETQudI9ISI/tKDtVuDRLssTEXHZ0+lI9avAPSLyf4CMiNwNvBr40U2TLCLiMqWjkU5Vvwl8N2a74/djfHYvVtWHNlG2iIjLktX2HM8CvwFcj0nj+V1VrfdCsIiIy5XVRro/xUwjnwJ+Avj9TZcoIuIyZzWl+2HgNlV9R/P5HZsvUkTE5c1qSpcLNgtR1dNAYfNFioi4vFnNehlvFpiVZV6jql/aLOEiIi5HVlO6cdqzCKYWvFZgf7eFioi4nFlR6VR1X4/kAEBEbsdUjbaAv1bV9/by+yMiekHHWQabjYhYmLLtPwxcC7xeRK7dWqkiIrrPRaN0mODpZ1T1uKo6wMeJIl4iLkMuJqXbhSkFEXCm2RYRcVlxMWUJLLUBuC46SOQtwFuaL+si8timSrU6w8BkJEMkA/CYql6/2kEXk9KdwdTUDNgNnF14kKr+FfBXACLysKre2BvxliaSIZKh9fs7Oe5iml4+BBwUkStFJAm8DvjcFssUEdF1LpqRTlUbIvILwBcwLoP3q+rjWyxWRETXuWiUDkBVPw98fg0f+avNkmUNRDIYIhk6/H4xm6xGRET0iotpTRcR8bzgklU6EbldRI6IyDMi8utbLU/E8xcReb+IjHfqvroklS4KGYu4yPggcHunB1+SSkcUMhZxEaGqXwGmOz3+UlW6KGQs4pLlUlW6jkLGIiIuRi5VpesoZCwi4mLkUlW6KGQs4pLlklQ6VW0AQcjYk8AnopCxiK1CRD4GfAM4LCJnROTnVjw+ikiJiOgtl+RIFxFxKRMpXUREj4mULiKix0RKFxHRYyKli4joMZHSPQ8QkQdE5K7m858SkXu7fP59IqIi0lFStIjsEZH5ZuD6845I6bqAiJwQkTERybW03SUiD2yhWEuiqn+rqrf18jub/59bW2Q4pap5VfV6KcfFQqR03SMO/JeNnkQM0e9yGRP9uN3jfwBvF5GBpd4Uke8TkYdEZLb59/ta3ntARH5HRB4EqsD+5nTtrSJyVETKIvIeETkgIt8QkTkR+UQzBA4RKYrIPSIyISKl5vPdy8jxsyLytebzdzSnecHDFZEPNt8riMjfiMg5EXlORH47mA6KiCUivy8ikyJyHHjVcv8UEfkIsAf4h+Z3vGPhdLR5/b8tIl9vHvMPIjIkIn/bvNaHRGRfyzmvFpH7RGS6mcj8mo5/pYsBVY0eG3wAJ4BbgU8Dv91suwt4oPl8ECgBb8SMiK9vvh5qvv8AcAq4rvl+ApM18Tmgv9leB76I2SWpADwBvKn5+SHgTiAL9AGfBD7TIt8DwF3N5z8LfG2Ja7gCEzT+yubrzwB3AzlgG/AvwH9svvfzmN15r2he25eb8sZX+v+0vN7XenxTvmeAAy3X9nTzfxoHPgx8oHlsDpPW9R+a792AKTB73Vb3g04f0UjXXX4T+EURGVnQ/irgqKp+RFUbqvoxTKd9dcsxH1TVx5vvu82231PVOTVxpY8B96pJ3J0F/j/gRQCqOqWqn1LVqqqWgd8BXt6p0CKSwSjZ/1TVz4vIKCYr/5dVtaKq48AfYQLLAV4D/LGqnlbVaeB3O/2uFfiAqh5rubZjqnq/mjjbTwbXitkN+ISqfqD5v3oE+BRme+5LgouqBN+ljqo+JiL3AL+OCcQO2AmcXHD4SdoTb0+zmLGW57UlXm8HEJEsRiluB4rN9/tExNLOjBV/AxxR1d9rvt6LGW3PiYSpi7EWGXcukHfhta2H1a413yLbS0RkpuX9OPCRLsjQEyKl6z7vBB4B/qCl7Syms7SyB/jHltcbiTz/VeAw8BJVPS8iLwS+xdLJvm00izodBm5uaT6Nmc4ON0eahZyjPZ9xzypf082o+tPAP6nqD3bxnD0lml52GVV9Bvg74Jdamj8PHBKRN4hIXEReiymodE+XvrYPMxrMiMggRvFXRUR+uCnnj6lqreUazgH3An8gIv0iEmsacYIp6yeAXxKR3SJSxIzsKzFG93bsvQfzv3yjiCSaj5tE5JounX/TiZRuc/gtzIIfMGsuzFrkVzFbSL8DuENVu7XDzB8DGYxB4Zu0j6Ar8VpgBHiyxYL5l833fgZIYowaJeDvgR3N996HyWX8DmZU//Qq3/O7wG+IyIyIvL1D2ZakuWa9DbO+PAucB34PSG3kvL0kyqeLiOgx0UgXEdFjIqWLiOgxkdJFRPSYSOkiInpMpHQRET0mUrqIiB4TKV1ERI+JlC4iosdEShcR0WP+f/Rmg1kzFYY8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(nbin)*0.2 + 0.1\n",
    "xticks = (0, 0.2, 0.4, 0.6, 0.8, 1)\n",
    "xticklabels = (\"0\", \"\", \"\", \"\", \"\", \"1\")\n",
    "marker = 12\n",
    "\n",
    "fig = plt.figure(figsize=(3, 2.5))\n",
    "plt.plot(x, data['full'], \"-\", color='m', lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['full'], -1), \".-\", color='m', lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "plt.plot(x, data['null'], \"-\", color=\"gray\", lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['null'],-1), \".-\", color=\"gray\", lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "plt.hlines(0.5, 0, 1, color=\"k\", linestyle=\"--\", linewidth=1)\n",
    "\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1.1)\n",
    "plt.xticks(xticks, xticklabels)\n",
    "plt.yticks((0, 0.5, 1), (\"0\", \"0.5\", \"1\"))\n",
    "for side in (\"top\", \"right\"):\n",
    "    plt.gca().spines[side].set_visible(False)\n",
    "plt.xlabel(\"Normalized time\", fontsize=labelsize)\n",
    "plt.ylabel(\"Perceptron fidelity\", fontsize=labelsize)\n",
    "plt.title(\"Whisker + Eye\", fontsize=labelsize)\n",
    "plt.subplots_adjust(left=.22, right=.98, bottom=.22)\n",
    "savefig(fig, f\"fidelity-fullnull-allsessions\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOEAAAC3CAYAAADkZDmQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAviElEQVR4nO2deXycd3ngv8/cI2k0um3ZlizLlm3ZbJKasGkLBBYSCAn02AXKTViO5WjLblvY7ge6YUn4UHZLoQe0pQcBlh5cW2g2hcDSJKQQGgMhjWzL8iXJto6RZjT3Pc/+8b4zGd2j0YxGlt/v5zMfzfze3/t7n3f0PvO7nkNUFQsLi8Zha7QAFhbXO5YSWlg0GEsJLSwajKWEFhYNxlJCC4sGYymhhUWDsZSwDoiIisihVY69XkQeqqCN+0XkvtpLZ7HdsJRwHUTkv4nIg0vKxlYpe8167anqF1X1JbWW83pGRC6JyG2NlqNaLCVcn0eB54qIHUBEdgNO4MSSskNm3W1FUcbrGRFxNFqGtbCUcH2ewFC6m8zPtwL/BIwuKTuvqlfLzrvN7B1DIvIpEREAEblbRB4z34uIfEJEZkUkLCJPicizlgogIj4R+ScR+UPznKMi8m0RCYrIqIi8uqzu/SLyJyLyoIjEgX+3QntvEZHTIhIVkQsi8p/KjnWJyAMismC2/z0RsYnI+0Tkq0va+SMR+aT5/mERuU9Evi8iMRH5BxHpFJEvikhERJ4QkYGyc1VE3m1+R1ERuVdEDorID8z6XxIRV1n9l4vIk6Zc3xeRG8zyLwD9wD+Y132/iAyY7b9VRCaA74rI/xWRX1si/1Mi8ktLv58tR1Wt1zovDKX7L+b7Pwb+I/CRJWV/VVZfgQeANowHJADcYR67G3jMfP9S4EdmPQGGgV7z2P3AfUAn8C/AfWZ5MzAJvAVwACeAOeB42Xlh4LkYP7KeFe7nLuCgec0XAAnghHnso8CfYvzwOIHnm/V6gTjQZtZzALPAs83PDwPnzHb9wCngLHCbWffzwGeXfEffAFqB40Aa+H/AYNn5bzbrnjCvdQtgB94MXALc5vFLwG1lbQ+Y7X/e/L68wKuBH5bVuRGYB1yNfr6snrAyHsHo7cB4KL9nvsrLHllyzu+q6oKqTmAo8U0rtJsFfMBRQFT1tKpOlR3fY7b7ZVX9oFn2cuCSqn5WVXOq+mPgq8Ary877uqr+s6oWVDW19KKq+n9V9bwaPAI8ZN5DUaZeYL+qZlX1e2a9KYzh9qvMencAc6r6o7KmP2u2Gwb+EWN08B1VzQFfBn5miSgfU9WIqo4ATwMPqeqFsvOL9d8O/Jmq/lBV86r6OQyl/dkVvtNyPqSqcVVNAl8HhkRkyDz2RuDvVDWzTht1x1LCyngUeJ6ItAPdqjoGfB/4ebPsWSyfD06XvU8ALUsbVdXvYvSinwJmROQzItJaVuUujF/xPy0r2w/cYg7LFkRkAXg9sLuszuRaNyMiLxORx83h5gJwJ9BlHv5fGD3aQ+ZQ9bfLTv0c8Abz/RuALyxpeqbsfXKFz0u/g0rr7wd+c8k992H8SK1F6XtQ1TTwJeANImIDXruC/A3BUsLK+AHGEOkdwD8DqGoEuGqWXVXVi9U0rKp/qKrPxhiSHQbeV3b4z4FvAg+KSLNZNgk8oqptZa8WVX1XebOrXU9E3Bg95+8Bu1S1DXgQY8iJqkZV9TdVdRB4BfAbIvJi8/S/B24w560vB75YzT1XwSTwkSX33KSqf2MeX+1+l5Z/DuMH68VAQlV/UCd5N4SlhBVgDmdOAr+BMQwt8phZVtWqqIg8R0RuEREnxnwrBeSXVPtVjEWgB0TEizHXPCwibxQRp/l6jogMV3hZF+DGmKfmRORlQGnLxFwAOWQuJEVMefIA5tD2K8BfA/9iDrW3gj8H3ml+VyIizSJyl4j4zOMzGHPJNTGVrgB8nG3SC4KlhBvhEaAHQ/GKfM8sq3ZrohXjAQsB4xgLBb9XXkGNVYR3YPQGX8eYs70EeA1GTzwNfAxDsdZFVaPAr2MMzULA6zAWSIoMAd8BYhgjgE+r6sNlxz8H/Bu28CFW1ZMY88I/xpD5HMYCV5GPAh80h6q/tU5zn8eQ/3/XQdSqEHOlyMKiIkSkHzgD7DaH5NcUIvIm4B2q+rxGy1LE6gktKsZc0PgN4G+vUQVsAt4NfKbRspRjKaFFRZgLQxHgduCeBouzYUTkpRjz4BmMOe22wRqOWlg0GKsntLBoMJYSWlg0mB2lhHfccYdibNBaL+u1HV4VsaOUcG5urtEiWFhsmB2lhBYW1yKWElpYNBhLCS0sGoylhBYWDcZSQguLBmMpoYVFg7GU0MKiwVhKaGHRYCwltLBoMJYSWlg0GEsJLSwajKWEFhYNxlJCC4sGYymhhUWDaYgSisgdZiKTc0siPJfXeaGZAGRERJaGmLew2DFsecooMVJ1fQojYNBl4AkR+Yaqniqr0wZ8GiOJyoSI9Gy1nBYWW0UjesJ/C5wzE39kgL8FfnFJndcBXytGeFbV2S2W0cJiy2iEEu5lccKSy2ZZOYeBdjPn3Y/MgK0WFjuSRmQwlRXKlsbjcADPxkjc4QV+ICKPq+rZZY2JvAMjTDz9/f1VC6WqLKQWaPe2V92GhUU1NKInvIyR1qrIPoycCkvrfNPMLTeHkevhxpUaU9XPqOrNqnpzd3d31UKJCMFkkGAyWHUbFhbV0AglfAIjWeMBMx3ya1ickASMxCfPFxGHGbr8FuB0vQU72HGQucQc4VS43peysCix5cNRVc2JyK8C38JIffxXqjoiIu80j/+pqp4WkW8CT2GksvoLVX16K+Q73HmY04HT2G12WlzL8npaWNScHRUG/+abb9aTJ09uuh1V5VTgFAfaD9DkbKqBZBbXKSutfyzDsphZARFhuHuYC6ELpHLLUr5bWNQUSwlXwSY2hruGGZsfI5PPNFocix2MpYRrYLfZGe4eZnRulFwh12hxLHYolhKug8Pm4EjXEU4HTpMvLE0nb2GxeSwlrACX3cVQ5xCn505T0EKjxbHYYVhKWCEeh4fB9kHOzJ1hJ60oWzQeSwk3QJOzib7WPkbnRxstisUOwlLCDeJz+9jdspux+bFGi2KxQ7CUsAraPG10eDu4ELrQaFEsdgCWElZJZ1MnLa4WxhfGGy2KxTVORUooIr8gIo1we9rW9DT34LQ7uRK50mhRakImn+GpmaesofYWU2lPeC8wJSJ/LCK31FOga409vj0UtMB0bLrRomyKeCbO2fmzDHcN0+5t52LoYqNFum6oSAlV9UbgNiAJfNUM0vRBERmop3BbTTqXruq8Pn8fyWySQDxQY4m2hmAyyGRkkmPdx3DanXQ1deF1epkMT65/ssWmqXhOqKo/VdX3YTjkvgd4FXBeRB4VkdeLyDU9v1RVzgXPVW0neqD9AOF0eFNOwblCjmAySCwTI5vPVt3ORrgavUo4FeZo11FsZf/C3S27EZFrvoe/FtiQK5OIHATeYL4KwOeBCeDdwJSq/vt6CFkpm3VlyhVynA6c5njP8UUP5EYYnRul19dLq7t1Q+clsgkuhC6wq3kXmXyGdD69qiLaxIbb4cZld+G2m38dbhy2jU3bL4Qu4HV46fX1rlrn0sIlmp3NdDdXH7XgOqYiV6aK/msi8h7gjcAh4EvAG1X18bLjXwWu+YhoDpuDgx0HOTt/lqNdR6tq43DnYc7MncHmt1XsFDyXmGMuMcex7mMVKX9BC6Rz6ZKyxrNx0rn0qkbmdpt9mbJeCl1ib+vedWPqDLQNMDY/hsPmsOLv1ImKekIReQD4HPB1M0zhSnVeoqoP1Vi+DVErp975xDyxTIz9bfurOr/oFDzYPojX6V2z7vjCOCJCv7/6IFXrkS/kSecNpY2lY5yeO83ult247K4V6zc5m+jwdtDsai6VnZk7w17fXnxuX93k3IHUricEHlbVLy+7gshvqOrvAzRaAWtJZ1MniWyCQDxQ1TCs6BR8KnCKoY4h3A73sjoFLXB2/izdTd10NnXWQuxVsdvsNNmMCAGhVIgXDrwQp925av14Jl5arAFjhNDl7eLSwiUOdhy0og3UmEp7woiqLpvkiEhQVTvqIlkV1KonLDI6N8re1r1Vx5rJF/KcCpziaNfRRQ99Opfm7PzZLX2gQ8kQ07FpjnQd2fB8N5vPspBaIJgMMjo/yoB/gM6mTjq8Hev29Nc5m+8JReRFxXoi8u+WNDoIRKuT7drgcOdhRgIjHOk8smbPsRp2m52jXUc5M3eG4e5hHDYHkXSEyfBk6fNWMBWdIplLMtw9XNX5TruT7uZuupu7OdRxiJHZEbwOL4FEgEQ2ARjuXu2edvwef9WLWtcra/aEIlLcse3HWAUtosA08LuqujRcYcOodU8IRi8wOj/K8e7jiFT0w7aMdC7NWHCMdk87qVyKwfbBqtvaKBdDF3E73Ozx7alZm9l8ljNzZzjWfQy7zQ4Y1jahZIhwOkxBCwhCq7uVdm87HoenZte+xqjon1zpcPTzqlqzUPQicgfwBxghD/9CVX93yfFfxLDSKQA54D+r6mPrtVsPJQRjjnQ5cpkjXUeqOr+4UBOIB7h14NYt6SmKc86e5h46vLWfMRR/WFZb0VVVIukIoVSoFCzLbXfT4e3A5/ZdL71l7ZSwlphZmc5SlpUJeO2SrEwtQFxVVURuAL6kquvuGdRLCQEC8QDJXHLDq5i5Qo7RuVH6/H3Yxc54eJzhruG69oTF3vtA24FFK5y1JpFNcGnhUsX3k86lCaVChFNhFMUmNvxuP+3e9lVXaq9xNjcnFJHTqjpsvp9keb4IAFR1o2vrpaxMZtvFrEwlJVTVWFn95tWuvZV0N3czvjDOXGKOrqauis5JZBOcD57ncOfh0grpvtZ9nJ0/W3Wvuh7JbJJzwXMc6TpS9we73Mm5kn1Vt8PN7pbd7G7ZDRi9dTgV5nLkcslSyevwlnrL64W1VgbeXvb+DTW85kpZmZYZhYvILwMfBXqAu2p4/RVRVYLJ4JrbBfvb9nNm7gxeh3fdHqa4Ab/U+qbV3Uq+kGdsfoyhzqGayQ8QToW5Er2yKYufjeJz+9hV2MW54DkOdRza0Lk2sdHubV9kBJDMJgkkAszEZzbc3rXKqv+p8jmYqj6y2quKa1aSlQlV/T/mEPSXMOaHKzcm8g4ROSkiJwOB6g2oRYRMPsOpwKk1A/4e7jzMhdCFNUMgToQniGfiy+wxixQfvFo6BU/HpgkkAhVb3dSSdm87fre/Jp4XXqeXfn8/nd7O6yaez6pzQhH5cCUNqOp/39AFRX4O+JCqvtT8/N/Mdj66xjkXgeeYGZpWpRZzwlwhx4XQBVx2F/v9+1ec62TyGc7On122YlpcDOlq6qpoyDoTmyGdT2/aWubSwiUcNgf7Wvdtqp3NMh2bJlfI1UyOaDrKRHiC4e7ha3UhZ9Nh8PsqfG2UdbMyicghMZ9uETkBuID5Kq61IVQVh83B4c7DdHg7GAmMsJBaWFbPZXcZNpXBZ5xfM/kMI7Mj9Pv7K54z7mrZhV3sVTsFqyqjc6P4XL6GKyAYnheqWjPPC5/bx2D7IKcCp3Z08OWGJIQRkTuBT/JMVqaPlGdlEpH/CrwJyGL4ML6v3lsUqspIYITB9sGSFYuqcjlymUQ2wWD74LIN+9n4LJl8hlZ3K5PhSY50HalqA34iPIHL7iotWFRCrpDjzNwZBtoGtl32qIuhi/jcvop/jNYjk88wOje6JYtNNaa2WxQiMgy8Etilqr8qIkcAt6o+Vb2MtWWzw9HicLLT27nIZjSTz3A+eJ5Wdyt7Wxdn9v7h5R9it9l5du+zN7XtsJEHN5VLMTY/tq0fyrH5Mbqbu2nztNWkvaKb2aGOQ9eSqVztsjKJyKswsuXuxeihAHzA71cl2jbFJjaOdh0llUstWjRx2V0Mdw/jdXp5evZpYplYyQm4399Ps7OZZC65qWsfaD9AKBliLrHmtJdIOsL54HmOdR/btgoIMNQ5xHRsmmi6NpaNDpuD4z3HuRC6ULM2twuVznY/DNyuqu8EigkZfsoqKayvdfr8fcaccHZk0Vykw9vBse5jXI1e5ZvnvkmHp4NeXy9Huo5wPnh+0/OWoc4hktkkY/NjK4bbn43PMhObWWQutp050nmEifBEyb50s9jExrHuY1yJXiGUDNWkze1ApUrYg6F08Mx2grINNtHrRZunjaHOIc7MnSGWecZ2IJVLkc6l+bl9P8dUbIrZ+Cw2sXG48zCjc6ObXlLv8/exx7eHkdkR4pl4qXwiPEE6l2aoc2jL7E43S3mex2rj96zU5tGuowSTQWbj17wfOVC5Ev4Iw7O+nNcA/1JbcbYXLruL493HuRq9ynRsmvnEPBPhCY73HKfN28bxnuMUtMCpwCkUpc/fx7nguU1ft9nVzPEe47oTCxOcnT+L1+Glz1/NYnRjKeZ5PDt/tqZxcw52HCSZTe6IcJOVGnAfBR4CLgI/CzwMHAZeoqrbJkhlPW1Hn7jyBPFMnBcMvGBZT5Qv5LkQuoDD5sDj8FDQwrIFnGrIF/I8fvlxnDYnP9P7M1W5U1VKQQsE4gFCqRBeh7FhXssedyXPi1pwNXqVbD5bdRSEOlPz1dEm4OXAfgyzsweW2Hg2nHooYfkGvNvuZjw8vqp/YTQdZTw8TjafZaBtYFMxWYqOv4c7D2MTG2fnz7K7ZXdNvfBzhRyz8VnCqTA2sdHd3E27p51oxtgk3+tbPwbNRkjlUpwLnqu5VU8gHiCcDm9HM7ft6UVRT2qthMX9qcH2wZKtaDaf5ez8Wfr8fatGVLscucyT00/ygv0vqMoQuWgpcrTr6KJe43LkMslskoMdB6t+iLP5LNOxaWKZGA6bg57mHvwe/7J6qsqV6BWi6SiD7YMrhuioho16XlRKKBliJj7Dkc4j22nOvDklFJHvUcHCi6reujG56kctlbDoAb/SBryqciF0AY/Ds+qwM5VL8dD5h7hx140bGioVh4RDHSsvwMQzcS6ELiz6YViPVC7FdGyaVC6F0+ZkV8uuijf4s/ls6V5rNUSNpqNciV6pOqLdasQyMcYXxreTmdumlfDNZR8PAv8RI+LaOIan/ZsxrF3u2ZyctaNWSljsKQ62H1zzoZuNzxJKhhjqHFrxn57Kpfjx1I/xuXwMtA2s2SumcikmwhOlh30tClrgfPA8Hodn1cWaRDbBdGyaTD6D2264EG1mkzuSjtR0iBpKhphPztd8CFl05drK8CFrUFPP+seBt6rqSFnZMQwl/NmqRawxm1XCYg/ndXorDgdR7JmGOodWDOMQToWZjc/isDnIFXIMtg8uGmKmcinGF8ax2+zs9+/f0OLLfGKe6dg0hzsP47Q7iaajzMRnyBVyNDmb2NW8q2bDyCKXI5drNkQNxAPEs3EG2gZqI5xJcRpR7sfZIGqqhGEMc7VUWZkXI+p2W7US1prN2o6eCpxiX+u+FedIa5Ev5BmdH2V3y+4VQ0lMRacAY+/x4sJFOr2dtLpbGQ+P47Q56ff3V73yGYgHeOLqE7R72un395cyRdWTWg5Rp6JT5DVfcwP0om1tuS1wA6ipEn4DSAC/g+GE2wd8CPCp6iuql7G21CIM/maGMJcWLmET24rDyfPB83Q1deG0O/nJ1E+IZqI8t++5G164KTofzyXmUJQ2TxvdTd1MxaY2vWizUYpD1D2+PZuKYzMZnsRpd27IgL0SClrgdOD0motodaamStgBfBr49xieDznga8Cvrefjt5XUc5+wUuYT88zGZznceXjRsDOWjvHdS9/lcOdhhjoMj/ri3uL+tv1rKk5BC8wl5krJZjq9nXQ2dS47Z6OLNpl8hmg6SjQTLVm0KIqqlobkldin1mKIWmvPiyKqagRZbu6uS8Crdaj9FoWZeakbCKiuYNzYYLaDEsIzXg4HOw5S0AKT4UnjoW7Zw+j86KIN61gmxqWFS8t6k1whV9r/EqS0h7fe0K+gBS6ELuC2u0uLNtl8lmgmSjQdXRQ1wGl34nP58Ll9y+aziWyCqegUmXyGNk8bu1p2rflDkc1nubhwcU1n6PWotedFOeeD52lxtbCrZVfN216DTa+ODqjqJfP94GoNFAM2bQe2ixICLKQWeGziMfa07OGm3ptKD3DRQ+NY97FF9S9HLhPLxGh2NhPPxrGLfdU9vNXIFXKlnm0qOkUgHmCgfQCvw4vP7aPV3VpVDNCF1AIzsRkUXTeE4maHqGfmzrCvdV9dfCQnwhPYxV4Ta6YK2bQSRlXVZ74vYOwZLm1UVXXbmPNvByUsBltqcbWwr3UfVyJXyBayDLY/8zsWSoYIpUKLysDoTdL59LoPYK6QI5aJEU1HF3koOGwOfG4fPpfRs+UKOcaCY/Q099RkmFc0bQsmgzhsDvb49qw67L0SuUIkHdnwEFVVGZ0fxWFz0OntpM3TVtPN96noFOl8uuYrsquw+eGoiOxW1WsmS2S1SljQAt+58B2+feHbHO8+zvP7n4/T7sQmNmxiwy5246/NvurnSDrC1ejVkuNv+dAtmAxyKXSJA+0HUJRsPstEZAItKG3ethU9Lxw2By67i4IWyBay5PK50sNot9lLw0ivw7vuQ7qepU1BC0TSEYLJ4LLQgy2ulhXbz+azTMWmiGfiq84fNzJEzRVyRqo3M+Wby+4inUsTThsxSt12N11NXTUJhRiIB1hILdQ82t0K1EQJFyWCEZGvNToR6FpUo4QFLfC6r76OB8ceLJXdcegO7v+l+1FVClogr3my+SypXIpMPlP6m8lnCMQDTMWm8Dq9dHg6VnzQig/+VGyKPl8frZ5WHDYHVyJX2N2yG5/bh4hgExsFLaCqpYfSaXeWlH4pHocHr9OL1+HF6/SuubJbXLQ50H4Aj8NDKBliIbVAXvPYxEaru5U2TxuPjj/Kk9NPMtw1zE27byKRTaAoTpuTDm/HirkmktkkV6NXyeQz+D1+I8svUsqfOBef4+LCRdq97TQ7V+45HTZHKfGpy+4inokbztModrHjcXjIF/Kk8sactsnZRHdTd9UGCAupBaaiUxztOlpPM7eaKGFpSGp+3lZZmJZSjRI+dP4hXvmlV5YeGEEQEW7Zews37LqB3pZeelt62du6lz5/H53eTlx2F9FMlPnEPD3NPfS19mGzrb8tUPTG97l9paBIpwKnONRxqNTr2cRW0UOhqqTzaZLZJMlckmQ2ucyp2Ca20h5ZKpcilUtxOXIZr9PL8e7jtHnaSgtEK/0Y3Tl0J3/9H/4am9jI5rMEEgEC8QCpXIpcIYfX6cVj9yxaBY6kI0Z0AIHdzbvZ1bwLj9ODy+5iLj5HMpfc8BA1X8gTTocJJUNkC9lF5cXvq9XdSldT14aiDWyBmVtdesIdp4T/85//J/c9eh95zZMv5I2eyFyAKPZ+5RQtUfr9/Qx1DtHX2kefv49+fz99rX3rriKCMS+JZWIc6jhEXvObTtG9lFgmRjAZJJqOks6nKRQKi3rKhdQCgXiA3b7dpd7lsUuPcd9j95WGxsWH+5Mv/SS37jfMg4tDZLfDjdPmJJlLEs/EyasRbMHn8tHZ1Fly5yqfP/b6emlxtVQUUrISktkkwWSQWCZGQQukcimyhSxehxe3w02bp41Ob+e6blPFleyl6etqRE2UMIER/brY2N9jhKwvNa6q361axPWEWydxzFI20xMu5Suv/gq3D97OXGKOycgk/zrzr4zOjxJOhQmmgkyGJ5kITywLiei0O9nr21tSyn5/P/ta95UUda9vL26Hu+T2dKTzSOnBvBK9wpPTT3LT7pu4bfC2ipQyX8hzNXqV8bARoj+SjpDOG6mzE5kEC6kFQiljIai4IBRMBgkmg8QzcUSkNHQs72WKeOweBjsG6WvtM+6j+Nf/zGef24eqlpS/GG+nOK90291Mx6cXzR/TuTTj4fFNb/SXfw/F3jKdTxNJR0hlU7S4W2h1tdLV3EWbp23F77To61gHM7eaKOEl1vGkUNUDGxKrQipJHLOUWs0Jy4dhM7EZ5hJzdDd309Pcs+z8WCbGZHiSyYihlMX3k5FJJsOTTMWmli287GrZVXqYm5xNHOk6wpdHvsyT009iFzuK8pw9z+Fdz3mXoUTJZ5QomAwym5glmAgSToeJpCPkNY+s8v922V34PX58LmOLotXdSoe3g57mHnY1G722x+Ehno3z4Uc+jIigqijGfPiuobuwiY3p2LQR5TseKPV8Rfwe/4pK2t3UTZOzqdQLO21OPA4PiWyCXCGH3+Mnl88Rz8ZLQ9TiItlGf4xWIplNlr63hdQCsUyMFlcLXU1d9Pp68bl8pZ44X8hzeu40g+2DeByeWslQP39CEdkF3I1h1H14ww1Udo0NR+re7OroT6d/yo27b+TFB15MIBEozfmqSZldJJPPMBWdMhTUVNTLkcuL3heHjUvxOoyHV1G8Dq+R78/TTmdTZ2m/rt1jhNRvcbXgtrtx2p14HV78Hj9tnja6mrrwu/00OZtWHfolsgnG5se499F7eej8M1nPy3+MwJiHLqQWGAuOcTF0kenYtOFJkgoRSAS4Gr3KZHhy2ejAbrPT29LLvtZ9pejkPc09tHvaafO0lQzXu5q6+MB3P7DqD+JmKWiBhdQCs/FZZuOzxDKxkofJQNsAXqeXU4FT3PPwPXz7/LdrIUNdLGbuwEgUcyfwfeCF9donFJFXAneo6tvMz28EblHVX13tnM3uE8Yz8dJDtat5V91zyYPxYNz7yL38wQ//YFGENZvYeO2/eS3vvvndDLQN0OHtKPVSiWyCSDpCNBMt9Vouu6vU01UTClFVGQuO8cPLP2QyMskNu27gRQdetO6DVxwGLqQWStsb2UKWZDZJJBXhauwqVyJXuBy5zJXoFSYjk0xFp0qLSOVmci6bi+n4tDEfLXt+7zh0x6ZTBaxG+Z5rXvME4gF+NPUjbGLDZXeV7v8rr/4KLzn4ko02X9mEV1XXfPGMsfYEcA7DiHsf0ALk1zu/2hfwKox5YPHzG4E/WqHeO4CTwEm/31+MAKeAnjx5Uk+ePLmo7J577lFV1d7e3lLZiRMnNJ/P6w0vu2FR3U986xP63k++d1HZ3R+8W+//yf2Lym689Ua9/yf364233rio/P6f3K93f/DuRWXv/eR79RPf+sSishtedoN67vOo7JFnyn3oPd+9R1/1rlctqvuBL3xAP/TFDy0qe/2vvV6/ff7b2tnTWSo78qwj+vjk4/rLr//lRXUfeeoR/fT//vSisv/x8f+hZwJnKr6nN33gTYvK3vnxd+p9D9y3qOzEXSf0/d96v/Yc6imVtXS26Ece/Yi+6O4XLar7vPuepzf9zk2LyngByodQWp4pkz2i7nvdarvZtqiu630udbzesajM8QsOdd3rWlQmR0RdH3apHJZF5c4PO9X+C/bF138t6n6/e1HZ29/+dlVVPXHiRKmst7dXVVXvueeeZc9epc/6enPC3wHeB3wdw3fwn8qONQOROvaEWzYcLedi6CIislUWFcDieWlxdfbE7hO86znv4ubem+nz9625H1bQAvlCvrTCm9c8uUJu1fcr/c9n4jMksgn2+9c2JgdKhgoOmwO72LHb7NjF/Gy+Lz8uIqRyKRZSC0TSkVI67WZXM+2e9tI2yhd++gXe84/vWdQLKspfvOIveG7/c8nms+QKuUUjhnwhX9qzLb5yhRx2ux2nOEtyOGwO3HY3za5mmp3NeJwe3HZjX7J8mL7WQl29esL1/HbyGGERXUB9ndSWU0ocA1zBCLH4unpf9ED7ARZSC4zMjjDUObQlUa5tYuP3X/r73D54O4FEgBO9J3jxgRczE5/h6dmnKVBYMxSETWzY7DacVfyL8oU8Y0FjiX6lhada4XF4FiUIVVXi2TihVIgLoQuML4xzsOMgt+y9hccvP1467/n7n8+RriPMJ+ax2Ww4xLHoR8LlcOFz+3A73Ljt7tKGf7VzyNsGb+POoTuXzUtvG7ytyjtfn3XnhCKyG3gL8FYMRfws8FcYWZLq1hOa116WOGat+rW0Hc3ms4wFx1Z11K0VRSVo97SvaOFf0AI/mfoJM/EZXjjwwpo6qEbTUS4tXFo1KsBWMJ+YZzI8icvhwi52/G4/37n4HU4HTnPj7hu5ffB2Y764CcXaKEsX6rbV6qiIvBh4G/ALGD3V83e6AfdEeIJs3jDArrV5UzgVZjIyyaGOQ+sqQSQV4dHxRznQfqDq7E/lXI5cJpVLrRtHp14UrYVm47P0+fs2HNrjGqGuWxQdGIGe3qqqz9pwA3WiXl4URfecSpSlElSVSwuXABhoG6hYCfKFPP86+68kc0l6W3rp9/dv+Bc6V8hxdv5szTwrqiGYCPLY5GPs9+/neM/x7RCQqV5YcUdrSb6QLwUB3sy+YTEAbl9r34Zj2RSZCE8QTUcREdo8bfS29FakyMWed6hjqCEBkKLpKE/NPEUyl+TW/ltxObZvVqkaUZOFGQsTu83OcPcwVyJXSl7zG+2FZmIzhFIhhruGNxUKvt/fTygZYio2hdPm5FTgFLtadq3Zs40vjFPQwrIU31tBOBXmcuQy88l5DrQduCZzatQTSwk3yN7WvcQzcUZmRypOWJkv5DkXPEeru7VmAW/bvcbS/tn5swy2DxLNRBmZHVkW1KgYMXyPb09NQ9pXQjAZZCo6VXJDOtF7YttlFd4OWMPRKilogbH5sZL/3GoU55MH2w/WJcNsUY42Txs9zT1MRiaJZ4xYnqlciqvRq6W4pFtFIB5gNj5rWPkgLKQXGOoYuiZyKtYYa064FUzHpgmnwitG4R5fGCeveQ60Haj7EPBK5AqpXIrB9kHyhTzfm/geAM/rf96WKWAxfVx3czddTV2cC56jxdVScSDlHYg1J9wKdrfsxu/2MzI7Ugo1mM6lORc8t6VDwL2te4mkIzw5/SQAz97zbNx2NxdCF3DanQy0DdRln01VmYpNEUqG2NWyi+M9x0lkE4u+D4u1sZSwBnidXo73HOd88DzJbBKn3VmTvbyNkslnsIkNRREEt8PNka4jxDNxRudG8bl97PXtrUmvXNACVyJXiGai9Lb0sqfH6O2molNEM9GaOinvdCwlrCGKlh7w1fz76nJdVc6HjAQxN+6+sfQ5momWIqINdw+zkFrgVODUqr6RlVCMoxrPxkt+g8Xy4hz5cGddvNt2LNacsAbEMjEuhi5ysOMgTc6m0l5gv7+/7uHXi+EZDrQfWLbyOB2bJpKOLEuzNhufJRAPsLd1b8WBdvOFPBPhCVK5FH3+vkXXKgaRqnS1+DrCWpjZClYza1Mzw5PT7qybL1xx33GtlcfiD8TS0A1qJgGNpCPs9+9fde6WzWeNeyxk6ff3L7NdvRK5QiKb2NIcGNcQlhLWk0w+w9j8GL2+3jUNvOcT88zEZxjqGKrZKmVBC0bUNpePXl/vuvVzhRyjc6Mr9nwFLXBp4VIpxXdRUTP5TGmDf3/b/mXmekXD805v56YsiHY41upovZhLzJWSvqynWJ1NnfjcPs7On93Q8G81ikO/4tC3Ehw2B8d7jnMxdJFoOrrIYsUmNgbbB0uBeouxTwWh39+/onlbMUhVrWxpr3esnnADFDPkep3eqvLplRttV8PV6NVSqMRqh35ziTnmEnMc7jy8anZhu9hX/XGZDE+SyWfq4lWyA6noC7IG8RVSNFXb49tTdULLgbYB2jxtjMyOlOKxVEK+kOfM3BnsYl9VeSqlq6mL/f79jMyOLIupCobz7UoKmCvkOBU4RZOziYMdjXF/2qlYw9EKmAxPksqlarL31eZpo9nZXLHDcD0cb4v7mqNzo/Q096wb0KrR3hc7HasnXINsPlv69V/JLK1anHYnx7qPEcvEOB88v2LMFzCUfyY+w7N6nlXzuZdNbAx3DxPPxkvD5JUYXxhnPjnP8e7jlgLWCUsJVyGYDHJ2/ixDHUN1C33Y7++nu7mbkcDioWFx6Od1ejnUcaiuQ79+fz9+t59TgVPkC88E9c3ms4zMjuBz+6z5X52xhqNLKM90e7zneN2v1+puZbhruOQw7LK7tnzoV3SLOj13mgNtB8jkMw3xvrhesZSwjEQ2wfng+RWtT+pJucNwPBtviOOt2+HmWPcxzgfP47Q7t+QHyMLAUkITVWUqOtVQw+MtTOO8IjaxbUXiTIslbPnTJgZ/KCLnROQpETmxSr2/FJGfmnW+IiJ17ZpExDK9smgIjXjiXgYMma93AH+ySr3/oqo3quoNGCH4V81BYWFxLdMIJfxF4PNq8DjQJiLLDCBVNQJGzwl4WSdFm4XFtUojlHAvMFn2+bJZtgwR+SwwDRwF/qj+ollYbD2NUMKVlv1W7OVU9S3AHuA08CsrNibyDhE5KSInA4FA7aS0sNgitkQJReQ9IvKkiDwJXMVIt1Zkn1m2IqqaB/4O+A+rHP+Mqt6sqjd3d1suNRbXHluihKr6KVW9SVVvwsh7/yZzlfRngbCqTpXXN48dKr4HXgGc2QpZLSy2mkbsEz6Iken3HJDAyPgEgIg8iJFwZhr4nIi0Ygxffwq8a+tFtbCoP1uuhGpYK79nlWN3ln187tZIZGHRWKydaQuLBmMpoYVFg7GU0MKiwVhKaGHRYCwltLBoMDsq2pqIRIHRBovRBcxdx9e3ZHgGTyXp5HeaP+Goqt7cSAFE5GQjZWj09S0ZFstQST1rOGph0WAsJbSwaDA7TQk/02gBaLwMjb4+WDIUqUiGHbUwY2FxLbLTekILi2uOHaGEInKHiIyawaN+u9HyWFzfiMhficisiDxdSf1rXglFxA58CiOA1DHgtSJyrLFSWVzn3A/cUWnla14JgX8LnFPVC6qaAf4WI5iUhUVDUNVHgWCl9XeCElYcOMrCYjuyE5Sw4sBRFhbbkZ2ghJfZQOAoC4vtxk5QwieAIRE5ICIu4DXANxosk4VFxVzzSqiqOYwQ+d/CiE/6JVUdaaxUFtczIvI3wA+AIyJyWUTeumZ9y2LGwqKxXPM9oYXFtY6lhBYWDcZSQguLBmMpoYVFg7GU0MKiwVhKeJ0hIg+LyNvM968XkYdq3P6AiKiIVBS/SET6RSRmGuJfl1hKWGNE5JKIzIhIc1nZ20Tk4QaKtSKq+kVVfclWXtP8fm4rk2FCVVvMFHjXJZYS1gcH8N7NNmKmiLP+Rzsc6x9cH/4X8Fsi0rbSQRH5eRF5QkTC5t+fLzv2sIh8RET+GSN13KA5vHu3iIyJSFRE7hWRgyLyAxGJiMiXTJM9RKRdRB4QkYCIhMz3+1aR424Recx8/35zWFh8ZUXkfvOYX0T+UkSmROSKiNxXHD6KiF1Efk9E5kTkAnDXal+KiHwB6Af+wbzG+5cOX837v09Evm/W+QcR6RSRL5r3+oSIDJS1eVREvi0iQdOx+9UV/5e2C6pqvWr4Ai4BtwFfA+4zy94GPGy+7wBCwBsxeszXmp87zeMPAxPAcfO4E8Mr5BtAq1meBv4fMAj4gVPAm83zOzGyGjcBPuDLwN+Xyfcw8Dbz/d3AYyvcQx+GEfyd5ue/B/4MaAZ6gH8B/pN57J0YCVz7zHv7J1Nex1rfT9nngfL6pnzngINl93bW/E4dwOeBz5p1mzHc2N5iHjuBEfD3eKOfg428rJ6wfvx34NdEZGkO77uAMVX9gqrmVPVvMB7iV5TVuV9VR8zjWbPsY6oaUcMu9mngITUcmcPAPwI/A6Cq86r6VVVNqGoU+AjwgkqFFhEvhtL9gao+KCK7MKIW/GdVjavqLPAJDEN5gFcDn1TVSVUNAh+t9Fpr8FlVPV92b+dV9Ttq2Al/uXivwMuBS6r6WfO7+jHwVeCVNZBhy9hpEbi3Dar6tIg8APw2hmF5kT3A+JLq4yx2RJ5kOTNl75MrfN4NICJNGEpyB9BuHveJiF0rW/z4S4xI5h8zP+/H6I2njMzlgDGNKcq4Z4m8S++tGta715Yy2W4RkYWy4w7gCzWQYcuwlLC+3AP8GPh4WdlVjIennH7gm2WfN2NV/5vAEeAWVZ0WkZuAn7Cy8/MizCBZR4DnlRVPYgx/u8yeaClTLPbn7F/nMrX0GJgEHlHV22vY5pZjDUfriKqeA/4O+PWy4geBwyLyOhFxiMivYASoeqBGl/Vh9BYLItKB8UOwLiLyMlPOX1LVZNk9TAEPAR8XkVYRsZmLQsUh7peAXxeRfSLSjtHzr8UMxly2FjyA8V2+UUSc5us5IjJco/a3BEsJ68+HMRYQAGPOhjGX+U1gHng/8HJVrVUGoU8CXowFisdZ3MOuxa8A3cDpshXSPzWPvQlwYSyShICvAL3msT/H8OX8KUav/7V1rvNR4IMisiAiv1WhbCtiznlfgjE/vQpMAx8D3Jtpd6ux/AktLBqM1RNaWDQYSwktLBqMpYQWFg3GUkILiwZjKaGFRYOxlNDCosFYSmhh0WAsJbSwaDCWElpYNJj/D+2tjybUVjbzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(nbin)*0.2 + 0.1\n",
    "xticks = (0, 0.2, 0.4, 0.6, 0.8, 1)\n",
    "xticklabels = (\"0\", \"\", \"\", \"\", \"\", \"1\")\n",
    "marker = 12\n",
    "\n",
    "fig = plt.figure(figsize=(3, 2.5))\n",
    "plt.plot(x, data['whisker_unique'], \"-\", color='g', lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['whisker_unique'], -1), \".-\", color='g', lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "# plt.plot(x, data['eye_unique'], \"-\", color=\"firebrick\", lw=.5, alpha=.3)\n",
    "# plt.plot(x, np.median(data['eye_unique'],-1), \".-\", color=\"firebrick\", lw=1.5, alpha=.95, \n",
    "#          ms=marker, mec='none')\n",
    "plt.hlines(0, 0, 1, color=\"k\", linestyle=\"--\", linewidth=1)\n",
    "\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(-0.4, 0.8)\n",
    "plt.xticks(xticks, xticklabels)\n",
    "plt.yticks((-0.3, 0, 0.3, 0.6), (\"-0.3\", \"0\", \"0.3\", \"0.6\"))\n",
    "for side in (\"top\", \"right\"):\n",
    "    plt.gca().spines[side].set_visible(False)\n",
    "plt.xlabel(\"Normalized time\", fontsize=labelsize)\n",
    "plt.ylabel(\"∆Fidelity\", fontsize=labelsize)\n",
    "plt.title(\"Whisker asymmetry\", fontsize=labelsize)\n",
    "plt.subplots_adjust(left=.22, right=.98, bottom=.22)\n",
    "savefig(fig, f\"fidelity-whisker-unique-allsessions\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOEAAAC3CAYAAADkZDmQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqTElEQVR4nO2dd3RkeXXnP7dyUCmHltRStzpNz8AExgMzgEljMEMyThgGTDLBLAzYXmwffA7Baw+MWS9eOEvGZHsJBgzD7GBjGw/JpJ4BhglM55bUklpZqpze3T9+r6pLaoWSVFJJ6t/nnDqq+r1X79339L71S/d3r6gqFoulfnjqbYDFcrljRWix1BkrQoulzlgRWix1xorQYqkzVoQWS52xIrRsCBH5hoi8fIXtHxaRt22lTTsNsfOE9UdEzgJdQLGi+FOqelt9LFofIvIK4NWq+qv1tmUn4au3AZYyz1PVf6+3EZatxzZHtzEiEhSRaRG5uqKsU0TSItLhfn6uiPxMRGZF5L9E5JoVjqci8iYROS0ikyLytyLicbd5ROStInJORMZF5DMi0uRuC4nIP4jIlHuen4hIl7vtHhF5tYhcCXwYeLyIJERk1t3+KRG5vcKG14jISfe67hSRnkX2vU5ETojIjIh8QESkpjd1G2JFuI1R1SzweeD3K4pvBf5dVSdE5HrgE8AfAm3AR4A7RSS4wmF/C7gBuB54PvAHbvkr3NfTgANAA/B+d9vLgSagzz3P64D0Ilsfdst/oKoNqtq8+MQicjNwB/B7QDdwzr2+Sp4LPBa41t3vmStcy67AinD78FW3lim9XuOWfxp4canGAl4KfNZ9/xrgI6r6I1UtquqngSxw0wrnebeqTqvqIPBejKgBXgL8naqeVtUE8BfAi0TEB+Qx4jvknudeVZ1fxzW+BPiEqt7n/sD8Babm3F+xz9+o6qxr338C163jPDsKK8Ltw2+qanPF62MAqvojIAk8RUSOAoeAO93v7APeXCleTG3Vs8TxSwxVvD9XsW+P+7lymw8zYPRZ4F+Bz4vIiIj8TxHxr+MaF5zDFfsU0Fuxz1jF+xSmRt7VWBHuDD6NaZK+FPiSqmbc8iHgnYvEG1HVz61wrL6K9/3AiPt+BCPqym0F4IKq5lX1f6jqVcATME3Gly1x7NWG2hecQ0SimBr2/Crf29VYEe4MPovpy/0+8JmK8o8BrxORG8UQFZHniEhshWP9mYi0iEgf8EfAF9zyzwF/IiIDItIAvAv4gqoWRORpInK1iHiBeUzztLjEsS8Ae0UksMy5/y/wShG5zu23vgv4kaqereYm7FasCLcPX3dHFUuvfy5tUNVh4D5MTfPdivJjmH7h+4EZ4CRmcGUlvgbcC/wM+H/Ax93yT2DE/h3gDJAB3uhu2wN8CSPAh4FvA/+wxLG/BTwIjInI5OKNqvofwNuALwOjwEHgRavYu+uxk/U7BBH5BDCiqm/dwDEUOKyqJ2tnmWWj2Mn6HYA7evjbwGPqbIplE7DN0W2OiPw18ADwt6p6pt72WGqPbY5aLHXG1oQWS52xIrRY6syuEuEtt9yimGF8+7Kv7fCqil0lwsnJS6amLJZtz64SocWyE7EitFjqjBWhxVJnrAgtljpjRWix1BkrQoulzlgRWix1xorQYqkzVoQWS52xIrRY6owVocVSZ6wILZY6Y0VosdQZK0KLpc7URYQicouIPOImBnnLMvs81U108qCIfHurbbRYtootj7bmBpD9APAMYBj4iYjcqaoPVezTDHwQuEVVB0Wkc6vttFi2inrUhI8DTrqJR3KYrDzPX7TPi4GvuElBUNXxLbbRYtky6iHCXhYmJRlmYUIQgCNAi5v77l4RWSrvgcWyK6hH8N+lkj4ujsfhA34F+DUgDPxARH6oqscvOZjIa4HXAvT399fYVItl86lHTTjMwsxAe7mYGahyn39R1aSqTmLyI1y71MFU9aOqeoOq3tDR0bEpBlssm0k9RPgT4LCb/SeASQhy56J9vgY8SUR8IhIBbsQkIrFYdh1b3hx1U23dhkk66cVkbn1QRF7nbv+wqj4sIv8C3A84wN+r6gObbBeaz+MJLJfVy2LZHHZVGPwbbrhBjx07tq7vqirx48eJHTmCyFLdVotlzVT1IFmPGRcRIdzTQ3pkcffUYtlcrAgr8MdiONksxWy23qZYLiOsCBcR6e8nNTS0+o4WS42wIlyEx+fDH4uRnZ6utymWywQrwiUIdXWRnZhAHafepmw56ZERcjMz9TbjssKKcBkifX2khofrbcaWUnT7w9mJCbRYrLc5lw1WhMvgi0TAcSik0/U2ZctIDQ4S7e8nsm8fyXPn6m3OZYMV4QpE+vtJXyaDNNnJSQItLYjXizcYxBsKkZudrbdZlwVWhCsgHg+BtjYyExP1NmVT0WKR3PQ0wfb2clm4p4fMhQu2WboFWBGuQrCtjfzsLE6hUG9TNo3k4CCRffsuKY/u20dycLAOFl1eVCVCEfkNEanHsqdtwW6eO8zH43j8frzB4CXbvKEQnmCQ3NxcHSy7fKi2JvxrYFRE3i8iN26mQdsRbzCIx++nkEjU25SaoqqkR0YI9y5eU32RSE8PmbGxy3K6ZquoSoSqei3wdCANfNkN0vRWEdm/mcZtJ8K9vaTOn2c3ObynR0cJd3ev6rAe6e+3o6WbSNV9QlX9uar+GWZB7huAFwCnROQ7IvISEdnV/UsRIbRnD5nR0XqbUhOK2SxOJoO/sXHVfX3hMJ5AgHw8vgWWXX6sSTgichB4O/AhIOS+/xhwG/Clmlu3zQg0NVFIp3FyuXqbsmFSQ0NE1hAOpLTCxDZLa0+1AzNvEJEfAj8CuoCXquoVqvpOVf0sJhbMr2+induGaH8/yR0+SJOdmsLf2IjHV/1Ym4gYLyI7Wlpzqq0JnwW8B+hR1der6g8rN6pqCvjtWhu3HfH4/fgikR07YqiOQ25qilDn2kO5+iIRxOu1zdIaU60I71HVf3LjhJYRkf9eeq+q36ypZduYcHe3GTGs8SBN8tw54idObOrgT2pwkHBf3+o7LkN4717bLK0x1Yrw7cuUv7VWhuw0Ir29pGvo4J0eGcEXjRLu7SVx8uSmCLGQTILHgy8cXvcxRITI3r27dt60HqzYKRCRm0v7icjTWBgz4wBw2bZLfA0NZCYmKGazS050r4Xs5CQKZbexcHc3iVOnaDh4sKbxblLDw8QOH97wcXzRKNmpKQqJBL6GhhpYdnmzWs/84+7fIPCJinIFxoA3boZRO4Vofz+J06c39GDn5+cpJBJE9+8vl/kaGgh1dZE4fZqGAwdqIsT06Cihri7EU5uZpEhfnw2MVSNW/I+o6oCqDgD/WHrvvg6o6hNUdXG80KpYLSuTiDxfRO53szIdE5FfXc95NhvxevE3N5OdnFzX9wupFJnx8QUCLOGPxQh1dpI8c2aDVoKTy1FIpQg0N2/4WCVEhEhvr22W1oBqPWZqlguiIivTs4CrgFtF5KpFu/0HcK2qXgf8AfD3tTp/rQl1dJCbnl7zagMnlyM1NETDgQPL7uOPxQi2t5PYoBCTQ0NENyFFgK+hAVRNX9OybpYVoYg8XPF+SEQGl3qt45yrZmVS1YReHJmIcmmuim1FuK9vTTWCFoskzpwxfb5Vmof+xkYCLS3rFmJudhZfNIrH71/X91ejFIFgN7nzbTUr9QlfU/H+92t4zqWyMl3iFC4ivwXcAXQCz6nh+WuOLxwm6/FQSKXMivwVUFXiJ08S3b+/6snyQHMzqJI8e3bJpuuy53IcMhcu0HjFFVV/Z62Ix2O8ac6fJ7J376adZzez7FOgqt+reF/LTLnVZGVCVf8Z+GcReTJmFcfTlzzYNsnKFNm7l/iJE6s+8InTp4n09a15RDXQ0oKqknRDUFRDaniYyAbmBKvFH4uRm5qq6kfIcinLilBE/qqaA6jqcnOIy1FNVqbK439HRA6KSLuboWnx9o8CHwUTBn+NttQM8XgIdXaSuXCBUFfXkvskz50j2N6+7gc12NoKqsbvcxVxFdJpUN0yUUT6+4mfOGFHS9fBSh2Svipfa2XVrEwickjc/6SIXA8EgKl1nGtLCbS0kI/HcfL5S7alR0fxRiIEmpo2dI5gWxvecHjVSHDpKoRaS8rNUptGYM2s1Bx95WacsJqsTMDvAC8TkTxmDeMLdYf0/CPuIE3lqGd2chJVJVyj/InB9nYyExOkRkaI9PRcsj1z4QKB9vaazQlWiz8WIzs5SSGd3pBXzuVG1VmZRORK4HeBLlW9TUSuAIKqev9mGrgWNpKVqZakzp/H39iIPxYjPz9Pbnp6TQMq1ZIZH0cLBcIVQnTyeZLnzhE7dKjm56sGdZyq+saXCbXLyiQiL8Bky+0FSnOGMeDv1mXaLqfULMuvMBlfC0KdneDxkB4bK5dV01/cTMTjIbRnDynbLK2aatsrfwU8Q1VfB5RmpX/OMimsL3dEhGB7OzPHjq04GV8Lwnv2gCqZCxfIz8/jDYU27Mu6UQJNTSa7VSZTVzt2CtWKsBMjOrg4naBs80n0eqHFItnJSUJ79iw5SFNrwt3dFAsF5h58kFB396afrxqiNi5N1VQrwnuBly4qexHw49qas/NRVRKnThHdv5+GgYEt860UxyHQ3r5uP9ZaI14voa4uO1paBdWK8E3A7SLybSAqIv+KmUD/k02zbIeSOH2acG9vOUyiPxbb9CxHxUwGp1Cg8fBhnGx22wgx0NxMMZOxSVdXoVoH7l8CRzGO128FPglcraonNtG2HUdycNBMxkej5bJQV5cZxdzEleipCi+ayN69FFIpslPbY1o1um8fKdssXZGqI/24cWS+uIm27GjSo6N4w+ElJ+NLTs6bsZIhMzFBoLUV8XrLZaX+mHg8BFpaan7OtSBeL8GODtJjY2YQyXIJK7mtfZcqBl5U9ck1tWgHstpkvC8SIeumWavlJLZTKJCfnV1yUXF03z6SZ8+CSE3XEa6HQEsL8VOnahKFYDeyUk1YuYbvIGZd36eBc0A/8HIWrra/LMnH4+TjcRoGBlbcL9LXR+LUKWJHjlyyrZBKkZuaopjLrcnvMjU8TLCzk8Tp08vuM3v//fibmvDHYlUfV3w+wt3dNV3+FN23j+Tp00tefyW5uTmy4+N4QyHCPT0LavjdSlUeM27M0Vep6oMVZVdhXM5u2kT71sRWe8wU0mnSa4jbUqoxg21t5GZmyM/Noar4wmECbW1rqiXy8Tj5+XkiK+SRKJE4fZpgR0dVQnQKBZxslsyFC+CunvcEAlXbtRLZqSm0UFjSwT07PU1uchI8HtN/dp/Lkhg3az3kJlPVL2q1fcIrgVOLys5gBmsuS5xcjtTgYNUCLKTTFDMZ5o8fJ9rfT7Ctjej+/evy7ywlclmtVinRcOAA8VPm37eSELVYJHHqFOLzIR4PwbY2UufPg2p5xHcjBNvaiJ86RaClpSzs7OQkuelpxBVZKb6Ok8uRHhmhkE4TP3ECbzhMuKdnVzZnq60J7wRSwNu4uBTpL4GYqj5vMw1cCxupCVWV5JkzRPbuXfWXX4tF4idP0nDw4LILc9VxFtR23nCYYFsbYMIbrtZ8XYlK39RqKc1fhvfsWTJCmqqSOHGCyL59eINBnHyezNgYxUwGf3MzxVQKp1Ag3N29oeVRTj5P4swZ/I2NFObn8QSDONksvsZGAm1tZC9coJhOl/uyvqYmsq43kFMo4G9oMGIMhdZtwxZS05rwFcAHgQcxKx8KwFeATVlpUS+80SjJwUF8DQ3LjuRVTsYvFmAxkyE7NUUxk0E8HvxNTUvWdh6fb93hAovZLE4utyYBgnGlazh40Aixu3vBNApA8syZBTWNx+8n0teHqpKdnMTJ5fCGw2QnJ0nncoS6utZsgzoOmYkJ8nNzFN1BKk8wSLCri8zoKOmhIULd3eXBq9zMDKlz5xCvl+jAAIV4nOzUFHMPPYS/sZFwT8+uWERc9SoKADfzUgcwoarbLgTzRmvCie9/n1BnJ8G2NrKTk0T6+y8ZzYyXapNo1NR2s7OmtnMcvKGQWe+3yq+0qhI/fnxdKw3iJ04QHRhYUx6JxedOnDxJuLe3/ACnhofxRiJm0fAK5ONxMuPjIIKIoIUCwc7OVddIarFIenSUYiaDLxKhkEiQnZ4mduQIuelpPH7/iv2+BbVyUxPi95drR19TE5GenjX/IGwRVdWEy4pQRPar6ln3/bJeyKq6/NDcFrPh5ujZs8z98pd4vF7CPT04uRyBlhai+/YhHg/JwUE8gQBaLJraTgR/czOB5uY19+1yc3MUUynCa/D1zE5NoY5DaIPrEtVxjBD7+igkEmixuKY5PCeXIz02RjGbNU+ZKoH29ktE7BQKpEdGcPJ5fA0NFObm8DU2Ij4f2fFx8rOztN1005ruXW5ujuzkJOLx4GtsJD8zQ3ZyEl8sRmTv3g0vmq4xGxZhXFVj7nsHM2e4+KCqqttmDHm9IlTHYeL732f+oYdovOoqAq2txjs9lzMDBzMziN9PuLub6N69BNraajLfFz91imhfX1Wjj6VBk2oHY1Y9nuMw89Of4ovFaFznMVWVrNu8dHI5xOsl0NaGv7HR+Iyq4mtoID87i9dtehcTCQKtreWFyaUVJ2vFKRTIjI9TTCbxhEI4+TzZ8XE8wSDRfftWrdW3iI31CVU1JiJ7VHVMVXdtAlB1HO774z9m/NsXY1l1POlJ9N96K/7GRqJ9feVwEumhITzBIOLx4PH5NjxsHu3rIzk0ROzgwVX3TQ4Orimf4GoUMxl8sZiJG5pK4eRyRkyFAmCWIwVaWlacpxMRQp2dhDo7ycfjpIaHmT52DPF6Cff24vH5zNxnMEgxlSLU1bUgEkCoo4P4iRP4m5rWfC89Pl/5WPl4nOzERFnMiZMnSQDRgYENtxq2ghX7hCIyr6qNFZ+/oqrbNgXaemrC8e9+l3tvuw2nUEDzedM08nq59t3vpunIEfzNzRSSSRoGBsgnEsz94hc4hQICeBsaCHZ0EOroWPfQeXpkBG80umIzqpBIkJuZqdli3dz8PHMPPGAWBTsOyeFhGq+8klB7O4iACIV4nNzMzMW+bnv7stdYSKdJj4yYH6ZgkPgjj5hkqoUCofZ2mq+5ZtkBFCeXIzk4WJNIAFoskhkfp5BIUMznKcTjOIUCsQMHCO3ZU48AVBtrjsLCJqn7eVpVt0U9vxTrEeHJj3yEEx/8IE4+bzLwuvdD/H4aDh4k2NpK580303z11TQdPYqvoYHs5CTZiQl8sRi5uTlyMzN4QyFCHR2E9uxZkyBLgzTLRSlbbfuqxy8Wyc/Pl4f41XHIjI3RfM01eMNhCokE+bk5E4qxtxdvKGQC+Zbug9cLPh+ayZj+iMdDsL0dfyxGIZEgPTZmVoxEIuSnpsxkuyoevx+Pz2dGi9Npwt3d5dpxMZnxceNj6k7h1IJCMlkWZG5+HvJ5ogMDRPr6tjL2Tk1EuLgm3HUiLNWEAKiiqqjj0P3MZ5IZGyM1PEx2asoIQITovn00PepRxK64gkBrK01XXUXDoUPkZ2fJjI+Tm5xEAgEjyK6uqobQC4kEudnZJYPnVlNTLjhWOk1+bq4cml48HpOVNxSimEox9+CDBLu68AWDiM+HPxYzUyUiJhPUwMCCPqpTKFBIJEytks/j5POkx8bITk4SaGkxuQ7zeRTweL34m5oIdXYueNDz8/MkBwfJTk8T2bPHjO4uan5udNR3OdRxyE5OkpmYIDM2hlMs0rB/P7GDB7fCJa4mIkxhol+XDvZVTMj68sFV9VvrNnE140RuAd6HmZv8e1X9m5X2X48Il+oTdj75yRx6/euJHT6Mx+cjPTbG2Le+RXZ8nPmHHiJ+4sSCNXuhzk5TU15zDY1HjxLp6yv7gyJCqKNj2UnyEolF83TgeuUMDy8bIqNcy83N4bi5MLzhsKnNikUz8uk4iAgev5/szAyR/n6Cyzh0O4XCxXnEWOySmjc3M0N2YgJvQwMKJE+dIn3hAk4+bybR9+wxx29tXbKPV8xmSZ47R3p0lGBrK7EjR8rXW8xmyytNCskk4vXiCQTwBAI1a0aWms3Jc+coJBI0DAzQeOWVNRd+BTUR4VlWWUnhZm2qOW7imOPAMzBeOj8BblXVh5b7zoZHRx9+mMajRwl1dhIdGLhEEIkzZ4gdOoR4vWSnp5m9/35mfvYz5h96iPlHHin7PooIoZ4eU0seOGDm0lpazFxce/uS81pOoUDyzJkFbnDxkyeJ7ttXfqALqZSp5VIpc488HlNrieBks+U1i95QyNRw0Wj51z41MoI3GFyxyaeOw8R3v8v0vfcS2buX5sc8Bl84jDoOTiaDv6UF8fvJjIyUszxV/rhk5+ZIDw6Sm51FHYdAS0vZDl8shjcUMvOLjkPy3DnmHn7YZHfatw9fKGSmHoJB0kNDxI8fJzowQPM11ywpQk8ggCcYxBMImOZwIFB1M7PkgBB/5BFS588T6e+n9TGPKc/vLh4t73jiE9fbhN24CJf9kkgXxovmVapamzHzS8/xeOAvVfWZ7ue/AFDVO5b7zkbnCTMXLpCbmyO6d+8lHiXg/pKfPUvs8OEl/ymp0VHG/u3fSA8PkxwcJHHqFJmxMTO5DQS7uoju20eoq4twTw8t111H06Mfjb+xEVQZuftuZu+/n9iRIzQePUoxncbf0ABeb1ncUnJwdvGGwxcFt8yDUvJ4CS8Ro7R8/YtbBKq0Pf7xHH7DG8qjnOnRUVCl4cABGhb9SC3GKRTITU+b2nN62vS3HcfY7zqtBzs7Ea+XxNmzCND06Efz0ze/mZl77zUDREDnU57C9e9974JrU1WcXM68slmK7l+WeJbF5zMirRBsZTO0mM0y//DDzD38MIHmZtqe8AQeeNvbFraMlrChSmovQrd5+Brg2cB/AU/drHlCEfld4BZVfbX7+aXAjap623Lf2egqirlf/hLN5/FGIngDAYIdHZd4vxTSadJDQzQcPrxsMyk9NkYhmSTa309ubo6Zn/6UuV/8gvjJkyROnyYzMmKaF6r4GxsJ9faSHR8ne+ECeDygSqS/n/4XvxiKRRMsynHA6y2HzQBTE0pJoO5LHcekaXP/FpJJComEyWXhllX+Le2XOHuW0W98w1xA6ZkQofOpT8UXiaDFIr5YzPh65nLlVHAev9/UcB5Puc+4OLiVx++/aJcqHp8PCQTMteVyph/u+uPOP/CA+6WLD3z7k5+84g/ISqjjwKJrZakoBx4PxUzGTG+cOmVaGX5/WbC/8v730/mkJ6319NW1o7U0GLHMi4vO2oPASYwT916gASiu9v31voAXYPqBpc8vBf7PEvu9FjgGHGtqaipFgFNAjx07pseOHVtQ9o53vENVVbu7u8tl119/vTqOoy9/wQsW7Hv/176mn3nnOxeUve/223XmF79YUHbLzTdr/NQpveXmmxeUj33rW/qet7xlQdnnP/hBve8LX1hQ9qyODr3zwAE9GAiUy1q9Xv36FVfore3tC/Z978CAvndgYEHZre3teteVV2qrz1cuOxgM6l1Hj+ozF92TTx86pG/bu3dB2W3d3fr1I0cWlD02EtE7DxzQx0ajC8rvvvZafWNv74Kyt/f362cWff+W1la9+9pr9VA4fPGafD69+9pr9cWdnQvv6aFD+r7DhxeUvai5We88cEBbvd6L1xQK6V1XXqnPbG5e/Zr27NG7rrxy4TVFo3rX0aOXXNPXr7hC39DVtaDsrV1d+ulDhxaUveY1r1FV1euvv75c1t3draqq73jHOy559qp91lfrE74N+DPga5i1g/9ZsS0KzG9iTbilzdGVcAoFclNT5BOJsquaeL3kpqdXXA3h5PMkh4ZIDw0RaG8v1xj+WMz0kdzm3PEPfpCTH/pQuaYAQISum2+m62lPW1DT4dZ+pVqw/PL5zGCG+1lVSY+OEt23DxwHp1gs135OPk8xlSovIvY3NJAcHubhd70LoNzcFY+Ha++4g+arr16QBNUbCuFraFjQ5yx9Lz8/X55fLIXX8Dc2LtuUq7y308eO8eDtt5dHokv34rr3vIe2xz0OzefNNEuhYGrblVpxIuZ++P2Iz4e6tW4xl0NUy3Z7/H68kQi+hga8weDC0fIKNrMmXG1YqIgJBhUAtnpVZTlxDHAeE2LxxVtsA2C8M0JdXYS6uspO27mpKXJzc0zfey/N115rmn7xOIVEojxS6fH5CLgOxpmxMeNJsnfvJc3Y5quvvticLBbNAytC/wtfuJ5/fDkUffcznmGagsWicTSfnzfi8PnM1InjUEgmcYpF/M3NNF93HTP33VcWQdtNN9H+hCfgj8Wq8mgRj4eA60sLZvQ2NzND8tw51HHM/WhtxdfQUL4Hlfc2un8/4/fcw8T3vlcWWNtNN9H0qEcZRwq/H384jLjeSouFraoU02kKySTFZBItFCi6HkCeUMisvFglYWrHE59I51OeckmfsOOJT1zT/2AtrNonFJE9mCVLr8II8ZOYsBZTbGJN6J772cB7uZg45p0r7b/RgZnUuXMXM86ucF+K2az5R6dSpMbGyE9PE2hrI9DcTKi7G3/l3GBJcK57WObCBYKdnfij0fI2LRZ56I47mPrxj8srFDYwGMD88eMm3KC7BtDJZsHrNX0x93ilAR28XrLj4zhuXzh+/LhZ4XHllRsZFVwSJ5cjNztrpk5U8QaDRpQV92vBSPUyNmixSCGVMv8DNwUcACIXa+lIZN3zgNXYUCWbMjDza8Crgd/A1FRP2g0O3CW0WCy7bYmI+WVNpcjH4xTdaQFwf1UrRiQz4+OoW5vkJicp5nJ4fD6CbW2XzA2WmolONkukr688R7XRf3wpzOHcQw+ZGiMWw9/YiNfvvzhFEA4jIji5HJmJCYrptPH06eysWQiLtVDMZMhNT5tcipiAWIHW1nIz3cnnKaRSFJPJhSH1Rcw1RSJ4I5Etzz61BjZ1iqIVE+jpVar66DUfYJPYqAjz8TiFeHzBP7xUY3ij0RUnjdMjI4jPZ/wxMQ9QdnKSQip1sW/U1FQ+RjGbJXXuXHlFwVopptMkz583K9FzObyhkFn53tlJdN++S0Z1S6H58/G4+YHo6tp26csKySS5mRmzRMrt0/miUbzRaHmOcYexeSLcrmxUhJnxcXyx2LofztTw8IIwFiXUccyc2dwcAP7GRoJurNDMxAT52Vki/f3Lzrs5+Tz5+XlSIyPkpqeNU3UgQKS3l+CePeWJ7mIut2CVgqqa887OLvD5tGwZVoT1IHn2bHmh71KoajlnoRaL+CIR/M3NZMbG8IRChNrby5HUSn9RxReNGtF1dFzS/MrH42QnJ8sjtfn5+XJkt0BLi1mStPNqkd1ATWPMWKokun8/idOnEa93yVpHRMxaPdchu5BKlQWTn50lff58efQyvGcPTUePrthfK2YyZEZHCfX2kjx71gRDamw00QAug5iduwErwk0gOjBQ9rpYyv2tEl8kgs9drFt0/T+rbQ7nEwlm7r2XUHc3+bk5s1RoZ8bnvKyxItwEypHNTpwwfb0qw/NVsw7RKRRMSIl4nPTICK3XX298Ty07lm07trvTEREaDh0iee7chlODldbExU+eJDU0hD8WQzwe2m680QpwF2Brwk1EPB5ihw6ZQMEHDqy5qZibnS2nOAu2tdFw8CAiQvLcOUKdndtuisGyPqwINxnxemk4cIDE6dMrRuwuUUgkyExMmMn/piYaBgYWjIamR0fxRiK2BtxFWBFuAR6/nwZ3sGa5tYilhCi+hoYFnjSL91HHWTYFm2VnYvuEW4QnECC6fz/xEycWLMrNTk4SP34cHIfYkSMmEvUSAszH4+Tn5qrKwmTZWdiacAvxBoPlPIW+piYKs7ME2tpWDehbzGbJjI7SUGUGKMvOwopwC1HV8nKnfDxO6/XXr/odp1AwITUOHbJeL7sU2xzdAlSV9NgYiRMn8ASDtDzmMWYe8cyZlb/nOOUwhNb7Zfdia8JNpBRot5BMmuBOFUlXAk1NaLFIcnCQ6DLh7RNnzphBmjosM7JsHVaEm4A6jslalEwS2rNn2SBFwdZWcBxSIyMLVj+AyT0RbGvbFfn3LCtjm6M1RItFUsPDJE6fJtDYSOzw4VWXDgXb2xGPh/TYWLksPTaGNxRadiWGZXdha8IaoMUi6ZERirncutJJh/fsITUyQmZiAo/Ph+bzJry85bLAinADOIUC6fPnzYr2np4NuZFFenpIDQ1RKBQ2lM/esvOwIlwHTj5P+vx5473S01P1KonVqFXqM8vOwopwDTi5HKnz5wEuSd5isayXLRehmBnn92FC6aeAV6jqfUvs93HgBkyIgOPufomttLVEMZs16Z+BSG+vnTKw1JR61ITPAg67rxuBD7l/F/MnqjoPICJ/B9wGrJgardYUs1kTbsLjIbJ3r121btkU6iHC5wOfURNh6oci0iwi3ao6WrlThQAFCMPKKdpqSWX650hfnxWfZVOphwh7gaGKz8Nu2ejiHUXkk5hm60PAmzfbsEIqRXp0FI/fb/ICbl7ySIulTD0m65fyQl6yllPVVwI9wMPAC5c8mMhrReSYiBybmJhYt1GqSm5mhob9+4n291sBWraMLRGhiLxBRH4mIj8DRjDp1krsdcuWRFWLwBeA31lm+0dV9QZVvaFjA4tdRYRIb691lLZsOVsiQlX9gKpep6rXYfLev0wMNwFzi/uD7rZDpffA84BfboWtFstWU482192Yft5JzBTFK0sbRORuTMKZMeDTItKIab7+HPhvW2+qxbL5bLkI3VHRNyyz7dkVHzcvIZzFso2wqygsljpjRWix1BkrQoulzlgRWix1xorQYqkzuypJqIjEgUfqbEY7MHkZn9/acJFQNenkd5tv1iOqekM9DRCRY/W0od7ntzYstKGa/Wxz1GKpM1aEFkud2W0i/Gi9DaD+NtT7/GBtKFGVDbtqYMZi2YnstprQYtlx7AoRisgtIvKIiJwUkbfU2x7L5Y2IfEJExkXkgWr23/EiFBEv8AFMAKmrgFtF5Kr6WmW5zPkUcEu1O+94EQKPA06q6mlVzQGfxwSTsljqgqp+B5iudv/dIMLlAkdZLDuC3SDCqgNHWSzbkd0gwmHWEDjKYtlu7AYR/gQ4LCIDIhIAXgTcWWebLJaq2fEiVNUCJkT+v2Lik35RVR+sr1WWyxkR+RzwA+AKERkWkVetuL/1mLFY6suOrwktlp2OFaHFUmesCC2WOmNFaLHUGStCi6XOWBFeZojIPSLyavf9S0TkmzU+/n4RURGpKn6RiPSLSMJ1xL8ssSKsMSJyVkQuiEi0ouzVInJPHc1aElX9R1X99a08p3t/nl5hw6CqNrgp8C5LrAg3Bx/wRxs9iJsizv6Pdjn2H7w5/C3wpyLSvNRGEXmCiPxERObcv0+o2HaPiLxTRL6PSR13wG3evV5ETohIXET+WkQOisgPRGReRL7ouuwhIi0icpeITIjIjPt+7zJ2vEJEvue+/3O3WVh65UXkU+62JhH5uIiMish5Ebm91HwUEa+I/C8RmRSR08BzlrspIvJZoB/4unuOP1/cfHWv/3YR+S93n6+LSJuI/KN7rT8Rkf0VxzwqIv8mItPuwu7fq/q/tF1QVfuq4Qs4Czwd+Apwu1v2auAe930rMAO8FFNj3up+bnO33wMMAo9yt/sxq0LuBBrd8izwH8ABoAl4CHi5+/02TFbjCBAD/gn4aoV99wCvdt+/AvjeEtfQh3GCf7b7+avAR4Ao0An8GPhDd9vrMAlc+9xr+0/XXt9K96fi8/7K/V37TgIHK67tuHtPfcBngE+6+0Yxy9he6W67HhPw91H1fg7W8rI14ebxduCNIrI4h/dzgBOq+llVLajq5zAP8fMq9vmUqj7obs+7Ze9W1Xk1frEPAN9Us5B5DvgG8BgAVZ1S1S+rakpV48A7gadUa7SIhDGie5+q3i0iXZioBX+sqklVHQf+N8ZRHuD3gPeq6pCqTgN3VHuuFfikqp6quLZTqvrvavyE/6l0rcBzgbOq+kn3Xt0HfBn43RrYsGXstgjc2wZVfUBE7gLegnEsL9EDnFu0+zkWLkQe4lIuVLxPL/F5D4CIRDAiuQVocbfHRMSr1Q1+fBwTyfzd7ud9mNp41GQuB0w3pmRjzyJ7F1/beljtWhsqbLtRRGYrtvuAz9bAhi3DinBzeQdwH/CeirIRzMNTST/wLxWfN+JV/2bgCuBGVR0TkeuAn7L04ucFuEGyrgB+taJ4CNP8bXdrosWMsnA9Z/8qp6nlioEh4Nuq+owaHnPLsc3RTURVTwJfAN5UUXw3cEREXiwiPhF5ISZA1V01Om0MU1vMikgr5odgVUTkWa6dv6mq6YprGAW+CbxHRBpFxOMOCpWauF8E3iQie0WkBVPzr8QFTF+2FtyFuZcvFRG/+3qsiFxZo+NvCVaEm89fYQYQANNnw/Rl3gxMAX8OPFdVa5VB6L1AGDNA8UMW1rAr8UKgA3i4YoT0w+62lwEBzCDJDPAloNvd9jHMWs6fY2r9r6xynjuAt4rIrIj8aZW2LYnb5/11TP90BBgD3g0EN3LcrcauJ7RY6oytCS2WOmNFaLHUGStCi6XOWBFaLHXGitBiqTNWhBZLnbEitFjqjBWhxVJnrAgtljrz/wH34KC+uM7JbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 216x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(nbin)*0.2 + 0.1\n",
    "xticks = (0, 0.2, 0.4, 0.6, 0.8, 1)\n",
    "xticklabels = (\"0\", \"\", \"\", \"\", \"\", \"1\")\n",
    "marker = 12\n",
    "\n",
    "fig = plt.figure(figsize=(3, 2.5))\n",
    "# plt.plot(x, data['whisker_unique'], \"-\", color='g', lw=.5, alpha=.3)\n",
    "# plt.plot(x, np.median(data['whisker_unique'], -1), \".-\", color='g', lw=1.5, alpha=.95, \n",
    "#          ms=marker, mec='none')\n",
    "plt.plot(x, data['eye_unique'], \"-\", color=\"firebrick\", lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['eye_unique'],-1), \".-\", color=\"firebrick\", lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "plt.hlines(0, 0, 1, color=\"k\", linestyle=\"--\", linewidth=1)\n",
    "\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(-0.4, 0.8)\n",
    "plt.xticks(xticks, xticklabels)\n",
    "plt.yticks((-0.3, 0, 0.3, 0.6), (\"-0.3\", \"0\", \"0.3\", \"0.6\"))\n",
    "for side in (\"top\", \"right\"):\n",
    "    plt.gca().spines[side].set_visible(False)\n",
    "plt.xlabel(\"Normalized time\", fontsize=labelsize)\n",
    "plt.ylabel(\"∆Fidelity\", fontsize=labelsize)\n",
    "plt.title(\"Eye position\", fontsize=labelsize)\n",
    "plt.subplots_adjust(left=.22, right=.98, bottom=.22)\n",
    "savefig(fig, f\"fidelity-eye-unique-allsessions\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOEAAAC3CAYAAADkZDmQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8lUlEQVR4nO29eZgkaV3g//lGRN5nVdZ9V/VZ3cDMDoygKKCCjPe64AKuoMixKqOL97Eq6/VDdtXFR3C9UVwfFQEFWUVQuRVhmGEG+u7qrvvKrLzvyIj390dkZdeR1V3dXdXZ3ROf58mnMt6IjPhGVnzzfd/v+z1EKYWLi0vn0DotgIvL0x1XCV1cOoyrhC4uHcZVQheXDuMqoYtLh3GV0MWlw7hK6HKoiMiLRGSx03LczbhKeI8jIrMiUheRnh3tXxQRJSITHRLNZZ+4Snh/cBV41eaGiDwTCNzKiUTEOCihXPaHq4T3B38GvGbL9vcA797cEBGfiPy6iMyLyJqI/K6IBJr7XiQiiyLyUyKyCrxLRAIi8qcikhGRcyLyk1uHlCIyJCLvE5GkiFwVkR/esi8gIn/S/OxZ4OFDv/t7HFcJ7w8+C0RFZFpEdOAVwP/dsv9twHHgQeAoMAz8wpb9A0A3MA68EXgLMAFMAS8BvnvzQBHRgL8Dnmye5+uBN4vIS5uHvAU40ny9FOcHweV6KKXc1z38AmaBFwM/B7wVeAT4KGAACpgESsCRLZ/5SuBq8/2LgDrg37L/CvDSLduvBxab758LzO+Q4WeAd2357CNb9r1x87Puq/3LHf/fP/wZ8EkcpXv3lvZeIAh8QUQ22wTQtxyTVEpVt2wPAQtbtre+HweGRCS7pU0HPrXHZ+du6i6ehrhKeJ+glJoTkavANwGv27IrBVSA00qppb0+vmN7BRgBzja3R7fsW8DpRY/tca6V5vFnmttj+7uDpy/unPD+4nXA1ymlSlvabOAPgP8tIn0AIjK8ZQ7XjvcAPyMiXSIyDDy6Zd/ngHzTkBMQEV1EniEiD7f57AjwQwd1c/crrhLeRyilZpRSj7XZ9VPAZeCzIpIH/gk4cZ1T/RKwiLP08U/Ae4Fa8xoW8K04Rp6rOD3tHwKx5md/EWcIehX4CM4w2eU6SHPy7OKyJyLyA8ArlVIv7LQs9yNuT+iyCxEZFJHni4gmIieAHwP+ptNy3a+4hhmXdniB38OxtGaBvwR+p5MC3c+4w1EXlw7jDkddXDqMq4QuLh3mvlLCRx55ROEsPLsv93U3vPbFfaWEqVSq0yK4uNw095USurjci7hK6OLSYVwldHHpMK4Surh0GFcJXVw6jKuELi4dxlVCF5cO4yqhi0uHcZXQxaXDuEro4tJhXCV0cekwrhK6uHQYVwldXDqMq4QuLh2mI0ooIo+IyAURuSwiP73HMS9qlvc6IyKfuNMyurjcKe54oqdmwZJ34hQaWQQ+LyIfVEqd3XJMHCex0CNKqfnNpLUuLvcjnegJvwK4rJS6opSq42Ty+vYdx3wX8H6l1DyAUmr9Dsvo4nLH6IQSDrO9YMhis20rx4EuEfm4iHxBRF6Di8t9Sifyjkqbtp35OAzg2Ti17wLAv4nIZ5VSF3edTOSNOOW3GBu79dojSimy1Sxdga5bPoeLy63QiZ5wke1VfkaA5TbHfFgpVVJKpXBKfj3Q7mRKqd9XSj1HKfWc3t7eWxZKREhX0qQr6Vs+h4vLrdAJJfw8cExEJkXEC7wS+OCOYz4AfI2IGCISxClMee6wBTvSfYRUOUWumjvsS7m4tLjjw1GlVENEHgX+Eae45B8rpc6IyPc39/+uUuqciHwYeAqntNcfKqW+fMhyoUyT44njnEueQ9d0wt7wYV7SxQW4z9LgP+c5z1GPPdauMtiNUUpRuHiRyPHjAJxNnmWya5KgJ3iQIro8vWhn/9iF6zHTREQIDA1RWV5GRJjuneZK5grVRvXGH3ZxuQ1cJdyCJxLBrtWwajU00ZjumebSxiXqVr3Torncx7hKuIPg2BjlBWcZU9d0pnunuZC6QMNudFgyl/sVVwl3oBkGnkiEWtpZqjA0gxM9JziXPIdlWx2WzuV+xFXCNvj7+6klkyjbBsCrezmWOMa51DlsZXdYusOlsrxMPZPptBhPK1wl3IPg6CjlxcXWtt/wM9U1xfnUee4ni/JWrOZ8uJZMoiy3179TuEq4B0YwCLZNo1JptQU9QUajo1zYuNBByQ6P8vw8obExguPjlObmOi3O0wZXCa9DcGyMysLCtraIL8JAeIBLG5c6JNXhUEul8HZ1IbqO7vOh+/3Us9lOi/W0wFXC6yCahjeRoJpMbmuP++N0B7q5krnSIckOFmVZ1NNpfD09rbbA0BDVtTV3WHoHcJXwBvgSCcxsFruxfYkiEUwQ9oaZy977w7bS/DzB8fFd7aHxcUrz8x2Q6OnFvpRQRL5NRDoR9nRXsHXtcCt9oT48uoel/FIHpDoYzEIBzeNB9/moW3WeWnuqNdTW/X40n496znVoP0z22xP+MrAiIu8QkecepkB3I7rPh+bx0CgWd+0bigxhK5vV4moHJLs9lFJUlpcJDA9Tqpe4uHGR6Z5pugJdXM1cBSA4NER1dbW1XONy8OxLCZVSDwAvBirA+5pJmn5ORCYOU7g7Ta1R23NfYHiY8tJS2+WJ0dgoFbNCspRs88m7l8rKCoHBQTLVDAv5BU71nsKje+gJ9hDwBFjIOb1/cGzMtZYeIvueEyqlnlRK/QROQO6bgO8EZkTkkyLyX0Tknp5fKqW4nL68p5+oiOAfGKC6stJ2/2TXJLla7raCght2g3QlTbFexLTMWz7PfrBqNexqlaQUyVVznOw5ibblXzgQHkBEWC2uYgQCaF4vZqFwqDI9XbmpeZ6IHAG+u/mygV8A5oFHgZcB/+mgBbxTiEjLPe103+ltD+Qm3liMQiqFXa+jeb279h/tPsqF1AUMzSDqi97U9ctmmSuZK/SH+smZOWpWbU9F1ETDZ/jw6l58evOv4cPQ9v/vLC8ssB5TBBEmuybbHjMSHWE2O0uylKRnaIjCxYsYx44h2j39e3vXsa94QhF5E/Bq4CjwHuDdSqnPbtkfBNaVUh2Ngr2deMJNymaZ+dw8J3tOtt1vmyal+XkiR4603a+U4nzqPKOx0X0HBafKKVLlFMcTx9sq/y4ZlE2tUaNu1alZzb+N2p5O5rqmb1NWO5tnfmOW8aln7iunzqWNS/QEe4goH7X1dUITE/u6L5f9xRPuVwk/BPwp8IFmmsJ2x3yDUuojNyXiAXMQSgiwUd6gWC8yHt9ttgdnLqUHg3hjsbb7lVKcTZ5lqmuKgCdw3WvNZecQEcZit56k6kZYttVS1kIlz9nHP8bQMx7Cq+/uzcHxDOoOdBPyhlpt51PnGY4Mo69n8cTjeCKRQ5P3PmJfSrjf8cvHlVJ/vesKIj+qlPpNgE4r4EGSCCYom2WSpSS9od3JowKDg+QvXMATjSKy+3veDAo+mzzLse5j+AzfrmNsZXNx4yK9wV4SwQQApbk57Hqd8NGjbc97q+iaTlBzMgQkr5znhc/9DvzhvYfLpXqJdCXNQt4xzBiaQU+gh9nsLFP9U1SuLrrD0gNkvz1hXim1678mImmlVPehSHYLHFRPuMmF1AWGo8Nth5WNYpF6JkNwdLTNJx0s2+Js8iwne07i0T2t9lqjxsWNixzpPtJKn1FZXkbzetGDQSpLSweuiJlKhuX1q4zTTXh84qY+a1om2WqWdCXNhY0LjHn6iVaEwePPvGFP/zTn9ntCEfm6zeNE5Gt3nHQKuK/NZccTxzmTPMOJxIltSgRghMNUk0msWg3dt7unA6cHOtlzkvOp80z3TmNoBvlanoXcQmsbHL9NBS23scDgIMWZGcJHjhyIIq4UVqg0KoxWQ4SO3fyw16N76A310hvq5Wj3Uc6sn8FvwmpyjlpzROvVvXT5u4j5Y/ua17pc47o9oYhcbb4dw7GCbqKAVeDXlFI70xV2jIPuCcHpBS5sXOB07+ldCqEsi+KVK0SOHbvuOWqNGpfSl+jyd1FtVJnqmmqdy8znqafTu4wdZqFAdX2d8NTUbSni1cxVfIaPrqKgBwJ44/FbPldLNsvkXPIcoxmD+MlpRIS6VSdTyZCr5bCVjSBEfVG6Al34Df9tX/Me5UANM+9WSh1YKnoReQT4LZyUh3+olPq1Hfu/HcdLxwYawJuVUp++0XkPQwnBmSMt5hc50XNi175qMomIbHN+3smmoSZZSvKCiRe0eopGuUxleZnI0aNtP2cWCtSSScJTUzct8+acsy/UR1wPU1pY2NOieyvUGjUuLDzJlNbXdnirlCJfy5OpZlrJsny6j+5ANxFf5OnSWx6cEh4kzapMF9lSlQl41Y6qTGGgpJRSIvIs4D1KqfZrBls4LCUESJaSVBqVtlbMwsWLztBR13fta9gNLqQuMBobRRedudwc0z3TKNOkePUqkRsYOMx8ntrGBuHJ9mt5bT/T7L0n45OEvCEKMzOExsbQPJ4bf/gmKJtlLn353zh17Ll4wjdejqk1amSqGXLVHAqFJhoxX4yuQNeeltp7nNubE4rIOaXUdPP9ArvrRQCglLrZSUarKlPz3JtVmVpKqJTa6qQZ2uvad5LeUC9z2TlS5RQ9we29XmB0lPLCwq4hZdksM5Oe4XjieMtCOhId4cL6OYayuqO4N7AweqJRlG1TvHp1X4pYMStcTl/mRM8JvLqXejaLEQoduAKCs5QxefI5nH3y0zzruS+94bDZZ/gYCA8wEB4AnN46V82xmF9seSoFjECrt3y6cD3DzBu2vP/uA7xmu6pMu5zCReQ7gLcCfcA3H+D126KUIl1Jt5YL2jEeH+d86jwBI7BtDc0IBKhpGo1y2YnI59oC/E7vm4g3QnitxEpvkBPG/laIvPE4KEVpdva6C+W5ao6lwlLrmsq2qa6tET2xexh9UEQDMYaOnObcmc9w6hlffVOf1USjK9C1zWGgYlZIlpOsldY42t1+mH6/sefP8NY5mFLqE3u9buGa+6nKhFLqb5pD0P+IMz9sfzKRN4rIYyLyWDJ56w7Um8aFs8mz1034ezxxnCuZK7u8U4IjI61wp/ncPKV6aZc/JkDxyhUGjp4mEeu/qaBgb1cXRjS6Z3zfanGVZDnJqd5TrWuWFxevu4RyUPT2jhLRg8wsn73xwTcg4AkwFhsjEUjc1/l8trLnnFBEfmk/J1BK/cJNXVDkK4H/oZR6aXP7Z5rneet1PnMVeLhZoWlPDmJO2LAbXMlcwat7GY+Ntx1i1a06Fzcu7rKYVtMbXF47z8D4iV1DVnAW4z3xeMvTZq24Rs2q3ZS3TG1jA6tc3qZcs9lZDM1gJDpy7T4qFcfFrE2w7mGgbJurT34Wz9Qoo7GDUfxCrcB8bp7p3ul71ZBz22nwR/f5ulluWJVJRI5K8+kWkYcAL7BxC9e6KZRSGJrB8cRxugPdnEmeIVvN7jrOq3uZiE9wKX0tz0zdqnPJXKZfYnR7druztXN16w/3o4t+U0HBvkQCPRCgvLiIUooLqQtEvJFtCghQWVi4I73gJqJpjB59JvXl1QOLrYz4Ikx1TXE2efa+Tr7ckYIwIvJNwNu5VpXpV7dWZRKRnwJeA5g4MYw/cdhLFEopziTPMNU11fJiUUqxmF+kbJaZ6pratWC/XlqnbtWJ+qIs5BY40XMCMS3H42XLskItlcKq1wkODbW99nxuHq/ubRks9kNpbYVLa+c5evLhXR491bU1xOPB133nnZmKV6+y5q8Si/W2HQ3cCnWrzoXUhZax6R7iYJcoRGQaeDnQr5R6VEROAD6l1FO3LuPBcrvD0c21tUQgsc1ntG7VmUnPEPVFGY5ur+z974v/jq7pPHvw2a2haXlpCU80iicS2XMxfidXM1eJ+CL7enCrjSqXNi4xYXdhKI3AFuW2TZPS3Nyea4+HjbJtCpcusdaj0RvqJe6PH8h5G3aDc8lzHO0+ei+5yh1cVSYR+U6carnDOD0UQAT4zVsS7S5FE42TPSepNqrbjCZe3ct07zQBT4Avr3+ZYr3YCgIei40R8oSoNK7lJ92s7mSWy1T3Gfoz2TVJppIhVb7utJd8Lc9MeoZTvaeIDI6AplFZvTb8K9/hYehORNPwDwwwXAuxWlylUDsYz0ZDMzjdd5ormSsHds67hf3Odn8JeIlS6vuBzRx4T7JHCet7ndHYqDMnXD+zbS7SHejmVO8plgvLfPjyh+n2dzMYGeREzwlm0jOtYzc9aDKPPXZT3i7HEseomBUubVxqm25/vbTOWnGNU72n0DXHMSAwMABKUV1bw8zn0f3+PX1Z7xTeWAy7VuNoaJz53Dxls3wg59VE41TvKZYKS2Qq90+q/v0qYR+O0sG15QTFXbCIfljE/XGOJY5xPnWeYv2a70C1UaXWqPGVI1/JSnGF9dI6mmgcTxznQuqCU/HXsqilUvgHBrDNm0tTMRobZSgyxJn1M5TqpVb7fG6eWqPGscSxXRbbwOAgVqNB7swZ/IODt3fjB0RobIzy/HyrzuP18vfcDCLCyZ6TpCtp1kvrB3LOTrNfJfwCTmT9Vl4JfO5gxbm78OpeTveeZrmwzGpxlY3yBvO5eU73nSYeiHO67zS2sjmbPItCMRob5dLGJYozM4QmJghPTrZNlXgjQt4Qp/uc685n57m4cZGAEbiu6V9sG29PD7XU9YezdwrRdaewzsoq0z3TXNy4eKB5c450H6FiVu7pdJOb7FcJfxj4FRH5BBASkX/EWUD/kUOT7C5BRDieOM5CboEvrX2JE4kTuxIinUicYCG3wEZ5A2M5TSYqrTSJnkjklqocaaIx1TXFQn6BbCV7XQOHVa1iNxpEjx3DrtVuSRFtZbNWXON86jxz2bkDWST3xuNY1SqYjVZI10GWlxuPO+u493oC5puxjgaBbwHGcdzOPrTDx7PjHIYD96bFtCfYg0/3MZebaxtfCLB++Rwr5LBDPibiEy13rPyFCzd01N7JZuDvZt6ZixsXGQgPtHWr2+lAXpqfxwiF8CX2dsEDx+K4XlonV82hiWPN7PJ3Uag7i+TDkeF95aC5HsqyKM7MEDl+nGqjyuX05W1ePQdBspQkV8vdjW5ud2cUxWFy0Eq4uT411TXV8hU1LZOLGxcZjY1uy6hWWVlBDAN/by+L+UW+uPpFXjj+QiK+CI1ymVoqRWhsf54xm54iJ3tOtgwwAIv5RSpmhSPdR1oP8V6hVKW5OTzRKN6u7UpkWiarxVWK9SKGZtAX6iPm3+1coJRiqbBEoVZgqmuqbYqO/VLPZLBqNQIDA5TNMrPZWaZ7pg88c8BaaY0TiRMHet7b5PaUUEQ+xT4ML0qpF9ycXIfHQSrhZgT8iZ4Tu1IJKqW4krmC3/AzHB1uuxhfbVT5yMxHeKD/Acbj45RmZ/H192MErr/GlSwlyVQzHOvebYABJ7bxSuaKk0RK81FqhkO1ozQ7iycexw77WS2uUm1U8Wge+sP9+84EZ1pm617HYmO3/IAXZmYIjoyg+3wUagWWCkt7ZrS7VYr1InPZubvJze22lfB7tmweAb4PJ+PaHE6k/ffgeLu85fbkPDgOSgk3e4ojXddPL7FeWieVXGRUxYlM7l6KqDaqPL7yOBFvhPHIKLK4TuT48V3HNcplCmvLLKbn8Bk+hiLtPWs2sZXNfG4eWUszMvkM9DY5UCtmhVQ5RXlpEX+sm/6+cfye/UW4i2EQGBzcFv6Ur+Vva4hqNxqUrlxp3X+mkmGjsrFrCFnP5aitr6P7/QSGhtrGaF6PzVCurelDOsiBRtZ/FnidUurMlrZTOEr4vFsW8YC5XSXc7OECnsANFQEcJ+n01UusJYRjiWNt0zjkqjnWS+sYmkE1lWI8Okqgt496JoOZy1Ezq6w00ni7u5nsO9Z2rtkOs1BgfW2ObNQxHHl0D4VagbXSGg27QdATpD/Uj8/wUbxyBV9v777SFNqNBnatRnVtDUQIDg9vS3S8mF+85SFqbWMD1Wjg7+8HnF6/ZJaYiE9QS6epp1KgOSFYNJ/LTWW8mXjIzWnE1jjODnGgSpjDcVerbmkLACtKqfitSnjQ3K7v6NnkWUaiI23nSDux6/VWZLyN4sLGBQbCA3QHdvtrrhRWsCpVgiWLmS9+mp6paboHx1gmg9fwMRYb27fybcpauHiRyPHjpMopPr/8ebr8XYzFxlqVonZSmJnB39d3XUVUlkXh8mXEMBBNw5dIUNvYAKUIDA+3nABuZ4hamJkhNDraUuz52TPUMmkG4yMo08SIRvH392PX61SWl7FqNbBt9ECAwNDQvh0RGnaD86nz23yBO8CBKuEHgTLw8zhBuKPA/wAiSqlvvXUZD5bbVcLczCWiYxNtU9xvO7b5sIaPHEHbEpg7m51FE42x2BjKtlu9nVKKpdo6/cNHMHQPTzz1L5T7Qjx/9Pk3HUGulGJ55gw5o44WDhL3x+kN9rJSXNlltNn5ueLMDIGBAYw2qSiUUhQvXSI4Po7u82GbJtXVVaxqFU88jlUuYzcaBAYHW4HLm0PUochQ2x+fdtjNtB6eaJRGPo/m87G0MYc/1sXw2Elqa2tYlQqI4I3HMWIxak1vILvRwBMOO8rov/HQ2la2k5BqhxHtDnKgStgN/A5OrQkdJ/nS+4EfulGM353kdpWwur5Oo1DACIcdd7A9jtv6sG7FqlZZXbxMKrPMRPcU/q5uvPE4omkUa0X+ZfZfOJ44znDJj7eri3nTGaaOx8eva0iwlU2qnCJdSWPX6oSzJiOn/sOuz2w12myN/N8m+8yMo0ihkJORu1agUC9QmJnB6OlGQgGUUq0huUfzUEulMLNZ9EAAZdvY9Tr+/v5Wr7rfIaqybSqrq1SWltC8XoxAACMaxROPc/nC5wkZQQanTrWMV/VMhtrGhrPwPzBAo1Bw4ikrFTzRKIGhodYPwp7XVMpJshzq3fcPxQFy8EsUzcpLvUBSqTbOjR3mdpUw+ZnP4O/rc4ZhqRTBsbFd1szCZm8SCjm9XTbr9Ha2je73O5V9Dad+w5HuI9jKZiG34DzU4SEubFxgumea8uUZoidOUKwXmc3O7upNGnajtf4lSGsNr3j5MqHJyW098FZsZXMlcwWf7mt52JiW6SharUDFrFC/uoBnsA9fOELEG8FI5glGu7aFPpXNMiuFFepWnbg/Tn+4H6tYorq+DiKICKrRwNfXhzcWw7RMrmavtg2GVpZFZWUFq1rFCAZpFIvU0mkix49TT6fRPB4CQ0PM5GfbRl5s65VjMcTjafWORixGcGjohvPdmfQMYW+Y/nD/vp6HA+K2raMTSqnZ5vs9vZA3EzbdDdyuEpZmZ8mdP4+m6wSGhrDrdbxdXYTGxxFNozQ/j+b1oiwLq1pFRJxI+WZvt5VsNcun5z/NUHiIBwcfvLau14zQOOodxiqXCTR9PRfzixTrRUKeECWzhC76rjW82sYGyrbx9+5OzQ+O4m72bCuFFZKlJBNdEwSMABFfhKgvit/wO4mjLl8mMDpKo1hEWdaePf/mvawV11CoVgrFyuoqVq3mPGVK4e3pwdfdvW2IGvdEqSwvY5smRjhMI5fDiEYRw6C2vo6ZzZJ43vO2fXfnU+cZiY7suYRSz+WopVKIpmFEo5iZDLVUCiMSITgysmd9EHD8b3XRd4WjHSK3rYQFpVSk+d7GWTPceVKllLo5G/IhcqtKqGyb5Gc+Q/7sWaKnTuHt7na80+t1aqkU9UwG8XgIDA4SGhnBm0jsud63mWwp7A0zEh1hKb+EaZtMdV37HctUMmSqGXrTapuRwrRMalat7QO41fOkYTco1osUaoVtEQqGZhDxOb2b3/DTsBtcSl+iL9S3K05R2TaZJ57AiESItlk2aYetbJKlJOlKGkMzGAwPYuTKmLkcdr2O6DreRAJPNMrsxSco1gpMDJ1EK1bQm/NQq1jE292Nr6enraOBUo6Ry9AMEoEEcX+8reHHbjSorq9jlUpofj+2aVJbX0fz+QiNj+8Z0LxSWKFm1ZiIT+zrnm+T2x+OisiAUuqeqQN9K0qobJsvvPnNrHzsn2nYDScX5lc9l56XfRueWBTd66OxkcVcXkX3eAhPTRIYGMLf24vH50cTDV108rU8y4XlVuDv1vlaupJmNjPLZNckCoVpmczn57FrJsF0Be/4yC65DM3Aq3uxlY1pm5SvzuHp70Xz+9A1nYg3QsQXIWAEbmidbOdp0yiXKS8uglL4R4bJFTbYSC5QrzlxkcGuHvoGJ4kEYm3Pb1omK8UVSvUSAU+AXiI0VtYpzc4iuk5geBhbg+XGBh50Bv19BAYGdg0bC5cuEZqYwNalVeKtbtXx6l5qjRq5mpOj1Kf7nPJsbQxZm0mSN0t619NOodbQ5GTbUUOylCRbzXIscf3M6QfAgSjhtkIwIvJ+pdRdWwj0VpRw7ZOf4DPf/3rshoVugxInMPWB//k2QkePoscimKUixsggpdwG+S+foVorY9kWDb9B1m+x6isTCEbo9ne3fWA3H/yV4gqjkVGi/iiGZrCUX6K3ZBCJ92DEomiiYSsbpRQNu0HdquPRPUi5hsoX8QxtHzL6DT8BT4CAESDgCVx3cXrTaDPZNYlerrP0hc9SiXmxLQtrNUniGc+id2iSTy18midXn2Q6OMG0b5xKrQQ+L/6eXnriA21rTRTzaeYvP0lDFJFIN4HFLFalQr1eQeuOUZ8cYL6+Slegi5Bnu8FImSb1xRVCR6ZahU+9updSveQET6PQRcdv+LFsi6rlrJIFPUF6g73bouyVZTnGtWIRyzRpFArYjQaRqSn8AwPb/jfZapaVwgone04eppvbgShha0ja3L6rqjDt5FaU8O/f+uMU//R96LZgWNe+NUsXSv1hJBZGPXyK+DOfyciDz6W3ZwSyBbKri+S1Gl0NH111L0YwiL+3F//AwJ5rWZvR+BFfhIHwgJPXZv0MAxvQffJUKyv1NqPGljXBne01q0bFrFBpVKiYlV3JkDTRCGg+rEKRcnaDarXMan4Fb7bMM776G+mKD6DKFcxcjvzMZX7x3G/zL0ufBARNwQsmXsBvfuPb0Txe6qUimWqadC2LGfFjB734TDA2ihh+PxLwY2VyFM0ymUoaPDp90QHiFQPdtAkPD1OK+ahS32VFra6vI7q+p8O5ZVvkajkylQymbW5r3/y+or4oPcGeVg6aRqnUUsh6Pg+mSWhykuDoaGsOegfc3A6lJ7zvlPCdv/tD9L79Qy0nWVGgKXjiiE2sCIm8EKtoKAGFIhlTJAd91EYSxPqHCZ88wcCpZzEi3fRWfYTLCt3ndxSyv7+tCX2lsEKxXuRo91EsZXFm9nGOeAYIj+528K4sL6OHQtc1OGylUamQSy2zseEMFU0aqFCAYCiK1BrULlymEvGQsYv0xgZQAR81v/CZuX/lT//f/2Q5atLQaT3cb3/xb/L8xLOxSyW0hsJQgkpmUdkcdiRIoy+G1TDRbEUwEKW3f4zY4ChKaM0fpVQlljHRClW8/b2sxxU+X3CbFbVw6dJ1rb67vhezQrqSplgvYiubaqOKaZsEjAA+w0fcHycRSKAh1FIpqskk1dVVbMsiPDFBpBl1spmvZ2f5ugPiQJSwjJP9evNkf4uTsr51cqXUv9yyiDcS7gaFY3ZyK0r4j5c+zCd+8HU8uHBtgf6LI3W+5id/mW/42leTrme5evmLXPnwB9lYmsU3lyS4nMdXqGPjuFdlgjZXehpc6W2w1Av6QC/j/gHG7G76I30khicZmjzF+OAJhiPD+AzHiXkzLKphN7jw1KdZC9Z5Mn2GBwce5MVTLwazQXlxcc8UGcqyqGYzLC1fYim3QKaSoSB1arqN1ahTK+QpVHLk6nmyVoHaRpq5UIkVcqQraUr1EiKC4MzHrIbJaM5DKmRR8tgg4Nf9THVPMRodZVLvZ8Lupq93lIHwID1pi2jBxosHTziM6o5S7Q3TiPoc/9NmSnuf7mO1tEqxmEVW00RzDYhHSPd5GU6M0x3oxqrVKC8uEhobo1EqIbqO5vWieb37Gi5u7S1rVo18LU/VrBL2hYl6o/SEnHLf1ZVVSnNzNIpFwpOTRKenscQpcX4Ibm4HooSz3CCSQim1/0olN8F+Csfs5FaU0FY23/XXr2LhE//EWNpgvqvBM6afzzve8Fd4/AHWimukyikS3jjBtSKRo0cRXaeWTpN96inWvvA5kk89QfnSDCqTwxawUeSjGnOJBpfDRTIBi6LfpmpALmBj9Cfo7R1jJDpC0BPkRM8J3vul97B+4Qyr3QqF4uGhh/nBgZeR6dLJmE51o3w2RSmTJJ9LkqvkyDUKZBoFLGx8loY0/1Omrih7FRWPwuPxEvPHGKmG8AcjeLq76A500xfqoz/UjyYafsNPySzxyx//JZ615OP4qsF6pMGF3jrPP/YiNCVk80nm7RSr1Q26igqfKRR9ikzIxheNOkpqDHC0FmeQGL3BHuI9QwTjCQLRON5oFG8gTMAToFQrUplfQJ9dxbJtrMFuJvqOorJ58Hr43Bf+gZUvP8HgyWfyVV//qm3hXJtoXi+az+cUVm3+3blMVDErZKoZMpUM2WqWYr1I2BsmEUiQqHlRs8tUlpYIjo0Re+BZXCxedYbKmpePvvcdLD75OUYe+Ape8vJH0fVbcgY/vHhCEekHvhfHqXt/9u2bv8ZNZ+q+1SUKW9l8dOajnD33r5z0T/B1z3sZGSmzUd6gL9TXSn9o1WqUZmf3DNAtr6yw+tGPUllcpDQ/T3FmhurqKgqFjcKKhyn1BsiFFKvBOl/qLvB4aI01VaBu1Xj+bIhjKR9zXXWuxGtE8dPwaVgaeHU/IW+ISCBGty9GIpgg5osRjnQT6e4jFO3G5/Hj0T0EjAAxf4y4P05PsIdAoYFP6QSH26+Plc0yF5MX+OSbXs/AuXSrPXWsm1f86h9jGB6spi9nqV5kvVtYCtdYrW+wXlonU82QLCdZLiyzkFsgW82i2RCtaoSqQqym0e/poS+QIBKIEfPFiMR6CPb2EwnFGcjqeDSD/gefzccefQPdV69lIlib7uEH3/tv25RAKYVdrzuvWg2r+Zc2z7IYhqOkPh94DAqqQrLqyF2sF/FaGvHlEtGlAqFEL+vHu/nUTzy67XtoJ8M+ORSPmUdwCsV8E/CvwIsOa51QRF4OPKKUen1z+9XAc5VSj+71mduNolj70hOkCusUtTq9sQH6R4/t8lG0KhUqi0uEju4d5lRdW8MqlgiMjdLI58k+8UXyXzpDaWaG0tWrVFeaqz5KYUQjrIcaVNbX6C7pqGZ3thppkP7ak7xg5Gvo8kQJaL7WA2VqiqpZoWbXUbozdzMML0FfEJ/mQ1cCto2yLBqlEo1iEW9XF6rZtvXv5nHFq1dZ/vA/OImqlEJDEE3oeeEL0IMhlGVhhMNoPi+qbqKsphHIMDA9Qs2uY5l1lNnANus0rAamZVJqlClaZUr1IuVqkbJZIt8oY4qNBuiWMw/3YjCU0ZlIac4a7ZavduF4BL33FiP8bRClnPGcrcBWrREDgMLGskzqtolWt4gky/Ru2CigoauWHNFf/G888oo33+zV92d23fzS93pxzVl7HriM48Q9AoQB60afv9UX8J0488DN7VcDv93muDcCjwGPxWKxzQxwClCPPfaYeuyxx7a1veUtb1FKKTU4ONhqe+ihh1Sj0VAPveDEtmN/++0/qn700e/c1vbaR1+u3v3e/29b24MPT6t3f+jX1YMPT29r/79/+LPqtd/3LdvafuRnv1f99v96dFvb13fH1QemJtURr7fV1qXr6m+OHVWvTCS2HfsbE2PqNybGt7W9MpFQHzx5QnUbeqvtiM+nPnjiuPqGHd/Jnxw9on5+ZHhb25sGBtQHjx/b1vZwMKg+MDWlHg4Ft7X/7TOm1Q8MDmxr+9mRIfXHR6e2tb2kK6bef/q4mvL7rt2Toau/OX1CvaJ3+z29bWJE/drEyPZ7ijvfSZd+7Z4mfV7118cm1NdHw9uO/b3JEfVTg33b2t7Yl1DvPTaxre3ZwYB679EJ9exgYFv7+45OqP+6Q6b/3t+vfn9ydFvbG97wBqWUUg899FCrbXBwUCml1Fve8pZdz95+n/UbzQl/HvgJ4AM4sYMf27IvBOQPsSe8Y8PRrVzNXEVEtnlU2I0G9Y0NzGKx5aomuk49nb5uzUDbNCktLFBZWMDb04Pu9yOahicSwYhEWksZF975Ts684+0oywLl/HyKJqivOM3oi19Kd7AHn9f5LJqG6Dqy+VfXQROUpmHroJrvTcukurKKd2QA02pgmyaWZdJo1LHNBqpSQdVN0DS0UJCNq5eo/Z+/QERzegwBROj/72/GP30MrGsJmvRAAE84jDccwTC86JqOLjo6GqpYxs4V0BQYHi/+rgTeWIyaXSdbzZKv5VvltEPeEDEjgpGv0CgW+X9/+w60P/jArv5Dfvw1PPzIKzFrFcx6Dds0UY0GKIVlW5i2iWk1X7ZJw26gGwYe3Zkv6oYXAw2PJfiVTsAI4vcF8GpefP4gRiiEEQ6j+3x8+K/eTv4tv7Xrf3mYPeGNBrkWTkY2L3DwVSavT6twDLCEk2Lxuw77opNdk2SrWc6sn+FY4hhe3YtmGPj7+/H397ectusbG9RzOdJf+ALxBx5whn6FAo1iEbv5wGqGgbfpYFxdXXU8SUZGdg1ju571LHwePw29AbaNphsYojP53d/H/NEgkdgYR24iFcRmKvrws78KzeNBWZbjaJ7Po2wbMQxn6cS2qRXyzKavcLTvhaz++yUyjz/uKLsIiec9j2e99GV4IpH9B9WGEtD0kVaWRT2ToTQ3h7JtIoZBonuwFUpVMktkq1kykmHOmmP4Zf+Jxz72z4xcvpY/bHUqwvMe/mrShXU0jxdP0IfuiSCG80Okazo+3ecs9GtePKaNVE1UudLyoAHQ/H6MYPCGBVNf8vJH+Z2//HP6z10LDlqb7uFlL99zFnTb3HBOKCIDwGuB1+Eo4ruAP8apknRoPWHz2rsKx1zv+Nt14C7PXUv1ZzbqzOXm6An27PLqt2o1GqUSVrlMeXUVM53Gm0jgjcfxDw7i2bo2uKlwStEol6mureHr68MTCoEIlm0xm75C7h3vpvj4k60Ihb4XvpCH3v52lMATK0+wVlrjRRMv2leAav7iRSfdYDMG0K7VQNfRmsG64PRmVZ8wV1hgzIrhsTX0YJDCxYsULl4kOj1N7/Off1MZ4m6EXa9Tz2Ydp3Gl0H0+in5YNlN4DS+66EQ9ET7xgd8n9eWnGHjWQ3zdd/wAQV/I+TFsLqgry6JRLjv/g0rlmkFGBN3vxwiHMYLBm06NsYllNfjoe9/B0pOfZ/iBh+8u66iIfD3weuDbcHqqr7kfHLg3UZblKM1mqI5SzK5epJbLMuzpafVgmt/vDClDIUTTqK6voywLTzxOvZn0STMMfInErgBapRSVlRXsWo1Gb5TF8gpHu4/i07yOE/m5c20VIF/N88m5TzLZNdk2+VSjXKa2sUHu7Fmkme/UE42iezzOcCsSQQ84fqZ2vc7slaeolPJM9B0j0N9/w0Dmg0YpxZcXnmB95SoDni6GwsP4I1G83d2tYbptmjTKZaxSyclfuomIc0/BIHoweKA/FAfMoS5RdOMkenqdUuoZN32CQ+J2ldAsFGgUCtv+4XogQMUHS2Zqzzwy4Hi2iGHg7+sDnAeolkrRKJcRTcPb1eXEwjWV+8raReqLy0xOPBP/jnSFe2HZFl9a/xKVRoU+o4uegkZ9PYlVr6P7/U7ke18fofHxXVbdzdT81VyWueICA+Mn6U/sdhy/E6TLaT698GnGY+Oc7jvd+kFplEqt9Igigug6RiiEHgo58+m7J5XhfnHzjt4s1fV1jEikbZiSZVutJMBby6Ztpby4iB4I7PKBVLZNPZ2mnstRb9RYsDeYGDtNPNRNNZnEzGYJjo3t6XNqmyZmPk95eZl6Os16YY2aNDCG+kmMTDLSM0l9Y2NX2kWllHPdbBbRNGphDysqy7HuYx1JgFSoFXhq7SkqjQovGHsBXuOeqjV4K7hKeBgs5Zcom+U9c7ls5vr0xuO79q0V10hX0kwa/TSyOWftLRjEE49TXV1F8/vx9/RgFgqY+XzrL0phhEIEh4fx9fYimkamkmGluEJvsJfV9VniVZ3R6YcAMPN5aqkUSim8XV1OKo3cPLaymYhP3PEeJVfNsZhfZKOywWR88sDKad8DuEp4WGyGBe1VsHJnikHLtricvkzUF2Uwsr1q0mZ2bts0nXliqdSyXnq7uvDG43vO12qNGueXv8RAwYPZEyG1OEN/oJeuniF8iQSi662M4UORodtOaX+zpCtpVgorjltcvcRU99S+kw7fJxzIEoVLGzarJl3auETMH9tV5jo0OUlxZgbRNMqGxXxuniNdR9oqrBEMYjTT41u1Gsq2b5ilexOtatJ/pcBqoE6XF575zBeyVFnlan2DCTtCtV5lubDcykt6p0iWkqyX1ukOdJMIJMjWspzuO93WB9TF7Qlvm9XiKrlqjmOJY9uGp0opLj7+CfShfo4MHFzgqN1oUEsmMQsFKsvLdD/0EJ5olKX8EtVGlamuKSzb4lPznwLgq8e++o4p4Gb5uN6QU6/+cvoyYW94X4mU71Pc4eidYjP1+maqwVqjxuX0ZQZDAxiLKUITE7dVPbdl2MlmnfR/vb1UVlcJDA9v6zU3S2mDU7/Pp/uYzc7i0T1MxCcOJXBVKcVKcYVMJUN/uJ+eYA9ls8xMembP1ItPI9zh6J0i4Alwuu80M+kZKmYFj+5preWpo3EnUfDU1E2lcgeoZ7NOBmzAl0g45c9EKM3N4e/r2zVsrVt1NNFQKATBZ/g40XOCUr3EhdQFIr4Iw5HhA+mVbWWzlF+iUC8wGB5kqM/p7VYKKxTqBU73nb5birLc9bhKeIAoVOsBl+aPoOg64akpileu7MrY3Y5GsUg1mXQW/2MxwpOT2xajKysr6MEgnui1jNJKKWYyM/gNPw8MPNDaLtQLDEWGCHlDTPdOk61mOZs8S2+ol75Q3y3d42Ye1ZJZYiQ60rJ02spuzZGPJw4luu2+xR2OHgDFepGrmasc6T5C0BNsFcMci4210q9vrV3RzsNjsyCKEQ7j6+trq6y1dBqrUtkWF7iZnmGya3KX5XG1uEq+lt9VZm29tE6ylGQ4OnzdCsBbsWzHwFRtVBmNjW671o2sxU9j3DnhnWA+N49pOXlFdyZiupK5gkf3MBa7Zv3cGRRcS6Wop9OtXJx7YRYK1FKpbVEba8W1Vi3DvSyPmz8QO1M3KOUUAc3X8ozHxvecu5mW6dyjbTIWG9vlu3qjddOnOa4SHiZ1q86ljUsMRgavW+Ngo7zBWmmNY91O2bNGuUxlaQkjFqORzeJNJG5Y1tqq1SjPzRE+5vRotrKdrG3eyK51x3Y07AYXUhfa9ny2spnNzmJaJhPxiZai1q06c9k5bGUzHh/f5a5n2RaX0pdIBBJ7ehC5uIaZQyNVTrFeWt/X+lsimCDii7QWzAOFOo1iEbNQoPuhh254LbvRcHrPo0cdo0xz6Lc59N0PhmZwuu80VzNXKdQK2zxWNNGY6ppq1ZLYzH0qCGOxsbbubZtJqo52H93Tl9Zl/7g94U1gK5uZ9AwBT4CR6P6dn5VSVNfWmFs6j9HTzdHxZzl+pDcICm7FBU5Oonm9LBeWW6kSb3XolyqnSJVTHE8cb3uOaqOKLvqePy4LuQXqVn3X8NulLfv6gtxB/D4p1UucWT/DUGRo3wqobJvK8jLFy5cxgkFOPvtF9PSPcWb9DIQDeGIxSvPze36+ePUqwdFRlKFzPnUeXfQ9lWe/9AR7GI+Nc2b9DNVGddd+v+Fvq4ANu8HZ5FmCniBHuq9fRtzl5nCHo/tgIbdAtVHd99rXZh0+q1TCPzBAYEtkQ9wfJ+QJcSl9iYHwACE7SHl5eVv0A0Bpfh5fIkFFt5hNnr1uGNXNsrmueSF1gb5QH4ng9eekuWqOhfxCx6Iv7nfc4eh1MC2TS+lL9If6b/igAtvq8O1VEXcrm5bVwVoQEWmVJ6usriKaRtJXo2bVONJ1eD3P1uiKdsxl57CUxWR80u39bh7XMHM7bEYA7Mf4oizLqa9er28rJ30jxmJj5Gt5rjTmGa4GkWQSzTAwaxUWgxX6PH2HHvYzFhsjU8lwNnmWE4kTraWOzeiLG1l/XW4fVwl3sLXS7em+09c/ttGgsrTkRLQPDe07+mErUV+U6Z5pLm5cJJwuYChho1vjWNedG/p1BboIeoKcS51jMj5J3ap3JPri6YqrhFvYdDxu532yFds0qSwtoWybwNDQrlQSN4uu6Uz3TrPkW6Jmm5zeUW76TuAzfJzqPcVMegaP7rnhD5DLweEqYROlFCuFlesaX+x6nfLSEoCjfLcRGdGOO1jGuS2aaHeicKbLDu64EorzE/9bOKn0y8D3KqUeb3PcHwHPwZncXmweV9x53AHKxZHuI233WbUaleVlAILDw3c8M5nL/U0nesJvBI41X88F/k/z705+RCmVBxCR3wQeBa5bGu2gsWo1KktLiKYRHBm56VAkF5f90Akl/Hbg3cpZG/msiMRFZFAptbL1oC0KKEAArl+i7SBpVCpUlpfRDIPg6KirfC6HSieUcBhY2LK92Gxb2XmgiLwLZ9h6FvixwxasUS5TWVlB83gIjY/vu2qsi8vt0Am3tXZmv7a9nFLqtcAQcA54RduTibxRRB4TkceSyeQtC6WUop7JEJ6YIDQ25iqgyx3jjiihiLxJRL4oIl8ElnHKrW0y0mxri1LKAv4KeNke+39fKfUcpdRzentvPaRGRAgOD99y/QIXl1vljiihUuqdSqkHlVIP4tS9f404PA/I7ZwPNvcd3XwPfCtw/k7I6uJyp+nEmOvvceZ5l3GWKF67uUNE/h6n4Mwq8KciEsUZvj4J/MCdF9XF5fC540rYtIq+aY9937Rl8/l3RiIXl87ixhO6uHQYVwldXDqMq4QuLh3GVUIXlw7jKqGLS4e5r9JbiEgBuNBhMXqA1NP4+q4M1/Dvp5z8/eabdUEp9ZxOCiAij3VShk5f35Vhuwz7Oc4djrq4dBhXCV1cOsz9poS/32kB6LwMnb4+uDJssi8Z7ivDjIvLvcj91hO6uNxz3BdKKCKPiMgFEbksIj/daXlcnt6IyB+LyLqIfHk/x9/zSigiOvBOnARSp4BXicipzkrl8jTnT4BH9nvwPa+EwFcAl5VSV5RSdeAvcZJJubh0BKXUJ4H0fo+/H5Rwr8RRLi73BPeDEu47cZSLy93I/aCEi9xE4igXl7uN+0EJPw8cE5FJEfECrwQ+2GGZXFz2zT2vhEqpBk6K/H/EyU/6HqXUmc5K5fJ0RkT+Avg34ISILIrI6657vOsx4+LSWe75ntDF5V7HVUIXlw7jKqGLS4dxldDFpcO4Suji0mFcJXyaISIfF5HXN9//FxH5yAGff0JElIjsK3+RiIyJSLHpiP+0xFXCA0ZEZkVkTURCW9peLyIf76BYbVFK/blS6hvu5DWb38+Lt8gwr5QKN0vgPS1xlfBwMID/drsnaZaIc/9H9znuP/hw+F/Aj4tIvN1OEfkqEfm8iOSaf79qy76Pi8ivishncErHTTWHdz8oIpdEpCAivywiR0Tk30QkLyLvabrsISJdIvIhEUmKSKb5fmQPOb5XRD7dfP+TzWHh5ssUkT9p7ouJyB+JyIqILInIr2wOH0VEF5FfF5GUiFwBvnmvL0VE/gwYA/6ueY2f3Dl8bd7/r4jIvzaP+TsRSYjInzfv9fMiMrHlnCdF5KMikm4Gdv/nff+X7haUUu7rAF/ALPBi4P3ArzTbXg98vPm+G8gAr8bpMV/V3E40938cmAdON/d7cKJCPghEm+014J+BKSAGnAW+p/n5BE5V4yAQAf4a+Nst8n0ceH3z/fcCn25zD6M4TvDf1Nz+W+D3gBDQB3wO+K/Nfd+PU8B1tHlvH2vKa1zv+9myPbH1+KZ8l4EjW+7tYvM7NYB3A+9qHhvCCWN7bXPfQzgJf093+jm4mZfbEx4evwD8kIjsrOH9zcAlpdSfKaUaSqm/wHmIv3XLMX+ilDrT3G82296mlMorxy/2y8BHlBPInAP+AfgPAEqpDaXU+5RSZaVUAfhV4IX7FVpEAjhK91tKqb8XkX6crAVvVkqVlFLrwP/GcZQH+M/A25VSC0qpNPDW/V7rOrxLKTWz5d5mlFL/pBw/4b/evFfgW4BZpdS7mt/V48D7gJcfgAx3jPstA/ddg1LqyyLyIeCncRzLNxkC5nYcPsf2QOQFdrO25X2lzfYAgIgEcZTkEaCruT8iIrran/Hjj3Aymb+tuT2O0xuvOJXLAWcasynj0A55d97brXCjew1vke25IpLdst8A/uwAZLhjuEp4uLwFeBz4jS1tyzgPz1bGgA9v2b4dr/ofA04Az1VKrYrIg8ATtA9+3kYzSdYJ4Ku3NC/gDH97mj3RTlbYHs85doPLHGTEwALwCaXUSw7wnHccdzh6iCilLgN/Bfzwlua/B46LyHeJiCEir8BJUPWhA7psBKe3yIpIN84PwQ0RkW9syvkflVKVLfewAnwE+A0RiYqI1jQKbQ5x3wP8sIiMiEgXTs9/PdZw5rIHwYdwvstXi4in+XpYRKYP6Px3BFcJD59fwjEgAM6cDWcu82PABvCTwLcopQ6qgtDbgQCOgeKzbO9hr8crgF7g3BYL6e82970G8OIYSTLAe4HB5r4/wInlfBKn13//Da7zVuDnRCQrIj++T9na0pzzfgPO/HQZWAXeBvhu57x3Gjee0MWlw7g9oYtLh3GV0MWlw7hK6OLSYVwldHHpMK4Surh0GFcJXVw6jKuELi4dxlVCF5cO4yqhi0uH+f8B6twp5QTG+M8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(nbin)*0.2 + 0.1\n",
    "xticks = (0, 0.2, 0.4, 0.6, 0.8, 1)\n",
    "xticklabels = (\"0\", \"\", \"\", \"\", \"\", \"1\")\n",
    "marker = 12\n",
    "\n",
    "fig = plt.figure(figsize=(3, 2.5))\n",
    "plt.plot(x, data['whisker_unique'], \"-\", color='g', lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['whisker_unique'], -1), \".-\", color='g', lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "plt.plot(x, data['eye_unique'], \"-\", color=\"firebrick\", lw=.5, alpha=.3)\n",
    "plt.plot(x, np.median(data['eye_unique'],-1), \".-\", color=\"firebrick\", lw=1.5, alpha=.95, \n",
    "         ms=marker, mec='none')\n",
    "plt.hlines(0, 0, 1, color=\"k\", linestyle=\"--\", linewidth=1)\n",
    "\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(-0.4, 0.8)\n",
    "plt.xticks(xticks, xticklabels)\n",
    "plt.yticks((-0.3, 0, 0.3, 0.6), (\"-0.3\", \"0\", \"0.3\", \"0.6\"))\n",
    "for side in (\"top\", \"right\"):\n",
    "    plt.gca().spines[side].set_visible(False)\n",
    "plt.xlabel(\"Normalized time\", fontsize=labelsize)\n",
    "plt.ylabel(\"∆Fidelity\", fontsize=labelsize)\n",
    "plt.title(\"Merged\", fontsize=labelsize)\n",
    "plt.subplots_adjust(left=.22, right=.98, bottom=.22)\n",
    "savefig(fig, f\"fidelity-merged-unique-allsessions\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "groups = [\n",
    "    # full model\n",
    "    data['full'][0],\n",
    "    data['full'][1],\n",
    "    data['full'][2],\n",
    "    data['full'][3],\n",
    "    data['full'][4],\n",
    "    \n",
    "    # null model\n",
    "    data['null'][0],\n",
    "    data['null'][1],\n",
    "    data['null'][2],\n",
    "    data['null'][3],\n",
    "    data['null'][4],\n",
    "    \n",
    "    # non-whisker model\n",
    "    data['full'][0] - data['whisker_unique'][0],\n",
    "    data['full'][1] - data['whisker_unique'][1],\n",
    "    data['full'][2] - data['whisker_unique'][2],\n",
    "    data['full'][3] - data['whisker_unique'][3],\n",
    "    data['full'][4] - data['whisker_unique'][4],\n",
    "    \n",
    "    # non-eye model\n",
    "    data['full'][0] - data['eye_unique'][0],\n",
    "    data['full'][1] - data['eye_unique'][1],\n",
    "    data['full'][2] - data['eye_unique'][2],\n",
    "    data['full'][3] - data['eye_unique'][3],\n",
    "    data['full'][4] - data['eye_unique'][4],\n",
    "]\n",
    "\n",
    "pairs = [\n",
    "    # full vs null\n",
    "    (0, 5),\n",
    "    (1, 6),\n",
    "    (2, 7),\n",
    "    (3, 8),\n",
    "    (4, 9),\n",
    "    \n",
    "    # full vs non-whisker\n",
    "    (0, 10),\n",
    "    (1, 11),\n",
    "    (2, 12),\n",
    "    (3, 13),\n",
    "    (4, 14),\n",
    "    \n",
    "    # full vs non-eye\n",
    "    (0, 15),\n",
    "    (1, 16),\n",
    "    (2, 17),\n",
    "    (3, 18),\n",
    "    (4, 19),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Item1</th>\n",
       "      <th>Item2</th>\n",
       "      <th>P</th>\n",
       "      <th>Sign</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(all)</td>\n",
       "      <td>(all)</td>\n",
       "      <td>4.090074e-24</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Full-1</td>\n",
       "      <td>Null-1</td>\n",
       "      <td>3.673581e-01</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Full-2</td>\n",
       "      <td>Null-2</td>\n",
       "      <td>3.791739e-02</td>\n",
       "      <td>*</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Full-3</td>\n",
       "      <td>Null-3</td>\n",
       "      <td>9.748826e-03</td>\n",
       "      <td>**</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Full-4</td>\n",
       "      <td>Null-4</td>\n",
       "      <td>2.873942e-05</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Full-5</td>\n",
       "      <td>Null-5</td>\n",
       "      <td>4.411170e-06</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Full-1</td>\n",
       "      <td>Non-whisker-1</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Full-2</td>\n",
       "      <td>Non-whisker-2</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Full-3</td>\n",
       "      <td>Non-whisker-3</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Full-4</td>\n",
       "      <td>Non-whisker-4</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Full-5</td>\n",
       "      <td>Non-whisker-5</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Full-1</td>\n",
       "      <td>Non-eye-1</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Full-2</td>\n",
       "      <td>Non-eye-2</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Full-3</td>\n",
       "      <td>Non-eye-3</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Full-4</td>\n",
       "      <td>Non-eye-4</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Full-5</td>\n",
       "      <td>Non-eye-5</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>NS</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Item1          Item2             P Sign\n",
       "0    (all)          (all)  4.090074e-24  ***\n",
       "1   Full-1         Null-1  3.673581e-01   NS\n",
       "2   Full-2         Null-2  3.791739e-02    *\n",
       "3   Full-3         Null-3  9.748826e-03   **\n",
       "4   Full-4         Null-4  2.873942e-05  ***\n",
       "5   Full-5         Null-5  4.411170e-06  ***\n",
       "6   Full-1  Non-whisker-1  1.000000e+00   NS\n",
       "7   Full-2  Non-whisker-2  1.000000e+00   NS\n",
       "8   Full-3  Non-whisker-3  1.000000e+00   NS\n",
       "9   Full-4  Non-whisker-4  1.000000e+00   NS\n",
       "10  Full-5  Non-whisker-5  1.000000e+00   NS\n",
       "11  Full-1      Non-eye-1  1.000000e+00   NS\n",
       "12  Full-2      Non-eye-2  1.000000e+00   NS\n",
       "13  Full-3      Non-eye-3  1.000000e+00   NS\n",
       "14  Full-4      Non-eye-4  1.000000e+00   NS\n",
       "15  Full-5      Non-eye-5  1.000000e+00   NS"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = kw_dunn(groups, pairs)\n",
    "\n",
    "kwtab = {\n",
    "    'Item1': ['(all)'],\n",
    "    'Item2': ['(all)'],\n",
    "    'P': [results.p_omnibus,],\n",
    "    'Sign': [psign(results.p_omnibus, append=False)],\n",
    "}\n",
    "\n",
    "flag = ['Full', 'Null', 'Non-whisker', 'Non-eye']\n",
    "\n",
    "for (item1, item2), dunn in results.pairwise.items():\n",
    "    flag1 = flag[item1 // 5] + \"-\" + str(item1 % 5 + 1)\n",
    "    flag2 = flag[item2 // 5] + \"-\" + str(item2 % 5 + 1)\n",
    "    kwtab['Item1'].append(flag1)\n",
    "    kwtab['Item2'].append(flag2)\n",
    "    kwtab['P'].append(dunn.p_corrected)\n",
    "    kwtab['Sign'].append(psign(dunn.p_corrected, append=False))\n",
    "\n",
    "kwtab = pd.DataFrame(kwtab)\n",
    "if saved == True:\n",
    "    kwtab.to_csv(str(figroot / \"KW-Dunn.csv\"), header=True, index=False)\n",
    "kwtab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}