CSV_to_stats-medial-bias.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "0ec00416",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "import cv2\n",
    "import pandas as pd\n",
    "import pyautogui as pa\n",
    "from scipy.stats import binom_test\n",
    "from matplotlib.markers import TICKDOWN\n",
    "import datetime\n",
    "from scipy.stats import norm, kstest, sem\n",
    "\n",
    "from lib.auto_GUI.auto_GUI_base import AutoGUIBase\n",
    "from lib.auto_GUI.auto_PhotoZ import AutoPhotoZ\n",
    "from lib.analysis.laminar_dist import *\n",
    "####################################################################################\n",
    "#    Aggregates SNR (.dat), latency (.dat), and ROI photos (png)\n",
    "#            scraped from PhotoZ\n",
    "#    Given ROI files in format slice#_rois.dat \n",
    "# 8/23/23 We want to add:\n",
    "#     Take into account the RELATIVE STIM LOCATION within every barrel. Does velocity depend on distance from stimulating location?\n",
    "#            New column: distance (um) from stim loc to relevant barrel pair boundary (1:2)\n",
    "#            New column: distance (um) from stim loc to relevant barrel pair boundary (2:3)\n",
    "#\n",
    "#     METHOD\n",
    "#         To calculate stim location to boundary, take these files produced by MovieMaker2.ipynb:\n",
    "#             electrode.dat\n",
    "#             corners_layer.dat\n",
    "#             corners_barrel.dat\n",
    "#         Calculate: Distance to lower barrel, distance to higher barrel,\n",
    "#                    Distance to lower layer, distance to higher layer\n",
    "#                       (if applicable for each)\n",
    "####################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "5b23cb4f",
   "metadata": {
    "scrolled": true
   },
   "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>Date</th>\n",
       "      <th>Slice/Loc/Rec</th>\n",
       "      <th>Slice Edge</th>\n",
       "      <th>Most Medial Barrel</th>\n",
       "      <th>Medial Dir</th>\n",
       "      <th>Slice Location</th>\n",
       "      <th>Stim Location</th>\n",
       "      <th>Actual Stim Time (ms)</th>\n",
       "      <th>Stim Barrel</th>\n",
       "      <th>L4 barrel 1 Latency</th>\n",
       "      <th>...</th>\n",
       "      <th>Distance to Barrel 2</th>\n",
       "      <th>Distance to Layer n1</th>\n",
       "      <th>Distance to Layer n2</th>\n",
       "      <th>Distance to Barrel A</th>\n",
       "      <th>Distance to Barrel B</th>\n",
       "      <th>Distance to Layer nA</th>\n",
       "      <th>Distance to Layer nB</th>\n",
       "      <th>Actual latency barrel 1</th>\n",
       "      <th>Actual latency barrel 2</th>\n",
       "      <th>Actual latency barrel 3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>78</th>\n",
       "      <td>9/7/2023</td>\n",
       "      <td>01_01_06</td>\n",
       "      <td>down</td>\n",
       "      <td>1.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4 Barrel 3</td>\n",
       "      <td>48.2821</td>\n",
       "      <td>3.0</td>\n",
       "      <td>53.0759</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>122.707482</td>\n",
       "      <td>162.175515</td>\n",
       "      <td>299.607853</td>\n",
       "      <td>56.903902</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>79</th>\n",
       "      <td>9/7/2023</td>\n",
       "      <td>01_02_03</td>\n",
       "      <td>down</td>\n",
       "      <td>1.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.2821</td>\n",
       "      <td>2.0</td>\n",
       "      <td>51.3798</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>100.126320</td>\n",
       "      <td>32.581904</td>\n",
       "      <td>128.666263</td>\n",
       "      <td>142.752381</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2.5640</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>80</th>\n",
       "      <td>9/7/2023</td>\n",
       "      <td>02_01_05</td>\n",
       "      <td>up/right</td>\n",
       "      <td>3.0</td>\n",
       "      <td>up/left</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.2821</td>\n",
       "      <td>3.0</td>\n",
       "      <td>50.1073</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>167.903664</td>\n",
       "      <td>41.159660</td>\n",
       "      <td>203.602330</td>\n",
       "      <td>40.800000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>81</th>\n",
       "      <td>9/7/2023</td>\n",
       "      <td>02_02_01</td>\n",
       "      <td>up/right</td>\n",
       "      <td>3.0</td>\n",
       "      <td>up/left</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.2821</td>\n",
       "      <td>2.0</td>\n",
       "      <td>46.5406</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>116.515753</td>\n",
       "      <td>24.081049</td>\n",
       "      <td>38.274314</td>\n",
       "      <td>99.841938</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>94</th>\n",
       "      <td>9/28/2023</td>\n",
       "      <td>01_01_06</td>\n",
       "      <td>down</td>\n",
       "      <td>3.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.2821</td>\n",
       "      <td>1.0</td>\n",
       "      <td>53.8574</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>40.480239</td>\n",
       "      <td>171.167934</td>\n",
       "      <td>225.909278</td>\n",
       "      <td>88.456883</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>169</th>\n",
       "      <td>10/19/2023</td>\n",
       "      <td>05_06_02</td>\n",
       "      <td>down/right</td>\n",
       "      <td>1.0</td>\n",
       "      <td>down/left</td>\n",
       "      <td>Anterior</td>\n",
       "      <td>L2/3</td>\n",
       "      <td>48.4228</td>\n",
       "      <td>1.0</td>\n",
       "      <td>52.5402</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>24.997803</td>\n",
       "      <td>181.431965</td>\n",
       "      <td>314.744870</td>\n",
       "      <td>103.797327</td>\n",
       "      <td>4.1174</td>\n",
       "      <td>2.8533</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>10/19/2023</td>\n",
       "      <td>06_01_01</td>\n",
       "      <td>down</td>\n",
       "      <td>3.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.5290</td>\n",
       "      <td>2.0</td>\n",
       "      <td>50.6833</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>118.898651</td>\n",
       "      <td>202.376373</td>\n",
       "      <td>288.667605</td>\n",
       "      <td>69.296479</td>\n",
       "      <td>2.1543</td>\n",
       "      <td>1.8040</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>171</th>\n",
       "      <td>10/19/2023</td>\n",
       "      <td>06_02_02</td>\n",
       "      <td>down</td>\n",
       "      <td>3.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L4</td>\n",
       "      <td>48.4103</td>\n",
       "      <td>3.0</td>\n",
       "      <td>50.3022</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>151.005549</td>\n",
       "      <td>215.305522</td>\n",
       "      <td>75.027164</td>\n",
       "      <td>306.897033</td>\n",
       "      <td>1.8919</td>\n",
       "      <td>2.6575</td>\n",
       "      <td>3.6943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>172</th>\n",
       "      <td>10/19/2023</td>\n",
       "      <td>06_03_01</td>\n",
       "      <td>down</td>\n",
       "      <td>3.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L2/3</td>\n",
       "      <td>48.5911</td>\n",
       "      <td>3.0</td>\n",
       "      <td>52.3050</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>322.923386</td>\n",
       "      <td>22.555042</td>\n",
       "      <td>70.935246</td>\n",
       "      <td>350.093066</td>\n",
       "      <td>3.7139</td>\n",
       "      <td>3.1466</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>173</th>\n",
       "      <td>10/19/2023</td>\n",
       "      <td>06_04_01</td>\n",
       "      <td>down</td>\n",
       "      <td>3.0</td>\n",
       "      <td>right</td>\n",
       "      <td>Posterior</td>\n",
       "      <td>L2/3</td>\n",
       "      <td>48.4633</td>\n",
       "      <td>2.0</td>\n",
       "      <td>53.4492</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>222.331315</td>\n",
       "      <td>55.298944</td>\n",
       "      <td>239.536169</td>\n",
       "      <td>171.112560</td>\n",
       "      <td>4.9859</td>\n",
       "      <td>2.5805</td>\n",
       "      <td>3.1151</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>69 rows × 40 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           Date Slice/Loc/Rec  Slice Edge  Most Medial Barrel Medial Dir  \\\n",
       "78     9/7/2023      01_01_06        down                 1.0      right   \n",
       "79     9/7/2023      01_02_03        down                 1.0      right   \n",
       "80     9/7/2023      02_01_05    up/right                 3.0    up/left   \n",
       "81     9/7/2023      02_02_01    up/right                 3.0    up/left   \n",
       "94    9/28/2023      01_01_06        down                 3.0      right   \n",
       "..          ...           ...         ...                 ...        ...   \n",
       "169  10/19/2023      05_06_02  down/right                 1.0  down/left   \n",
       "170  10/19/2023      06_01_01        down                 3.0      right   \n",
       "171  10/19/2023      06_02_02        down                 3.0      right   \n",
       "172  10/19/2023      06_03_01        down                 3.0      right   \n",
       "173  10/19/2023      06_04_01        down                 3.0      right   \n",
       "\n",
       "    Slice Location Stim Location  Actual Stim Time (ms)  Stim Barrel  \\\n",
       "78       Posterior   L4 Barrel 3                48.2821          3.0   \n",
       "79       Posterior            L4                48.2821          2.0   \n",
       "80       Posterior            L4                48.2821          3.0   \n",
       "81       Posterior            L4                48.2821          2.0   \n",
       "94       Posterior            L4                48.2821          1.0   \n",
       "..             ...           ...                    ...          ...   \n",
       "169       Anterior          L2/3                48.4228          1.0   \n",
       "170      Posterior            L4                48.5290          2.0   \n",
       "171      Posterior            L4                48.4103          3.0   \n",
       "172      Posterior          L2/3                48.5911          3.0   \n",
       "173      Posterior          L2/3                48.4633          2.0   \n",
       "\n",
       "     L4 barrel 1 Latency  ...  Distance to Barrel 2  Distance to Layer n1  \\\n",
       "78               53.0759  ...                   NaN                   NaN   \n",
       "79               51.3798  ...                   NaN                   NaN   \n",
       "80               50.1073  ...                   NaN                   NaN   \n",
       "81               46.5406  ...                   NaN                   NaN   \n",
       "94               53.8574  ...                   NaN                   NaN   \n",
       "..                   ...  ...                   ...                   ...   \n",
       "169              52.5402  ...                   NaN                   NaN   \n",
       "170              50.6833  ...                   NaN                   NaN   \n",
       "171              50.3022  ...                   NaN                   NaN   \n",
       "172              52.3050  ...                   NaN                   NaN   \n",
       "173              53.4492  ...                   NaN                   NaN   \n",
       "\n",
       "     Distance to Layer n2  Distance to Barrel A  Distance to Barrel B  \\\n",
       "78                    NaN            122.707482            162.175515   \n",
       "79                    NaN            100.126320             32.581904   \n",
       "80                    NaN            167.903664             41.159660   \n",
       "81                    NaN            116.515753             24.081049   \n",
       "94                    NaN             40.480239            171.167934   \n",
       "..                    ...                   ...                   ...   \n",
       "169                   NaN             24.997803            181.431965   \n",
       "170                   NaN            118.898651            202.376373   \n",
       "171                   NaN            151.005549            215.305522   \n",
       "172                   NaN            322.923386             22.555042   \n",
       "173                   NaN            222.331315             55.298944   \n",
       "\n",
       "    Distance to Layer nA  Distance to Layer nB  Actual latency barrel 1  \\\n",
       "78            299.607853             56.903902                      NaN   \n",
       "79            128.666263            142.752381                      NaN   \n",
       "80            203.602330             40.800000                      NaN   \n",
       "81             38.274314             99.841938                      NaN   \n",
       "94            225.909278             88.456883                      NaN   \n",
       "..                   ...                   ...                      ...   \n",
       "169           314.744870            103.797327                   4.1174   \n",
       "170           288.667605             69.296479                   2.1543   \n",
       "171            75.027164            306.897033                   1.8919   \n",
       "172            70.935246            350.093066                   3.7139   \n",
       "173           239.536169            171.112560                   4.9859   \n",
       "\n",
       "     Actual latency barrel 2  Actual latency barrel 3  \n",
       "78                       NaN                      NaN  \n",
       "79                    2.5640                      NaN  \n",
       "80                       NaN                      NaN  \n",
       "81                       NaN                      NaN  \n",
       "94                       NaN                      NaN  \n",
       "..                       ...                      ...  \n",
       "169                   2.8533                      NaN  \n",
       "170                   1.8040                      NaN  \n",
       "171                   2.6575                   3.6943  \n",
       "172                   3.1466                      NaN  \n",
       "173                   2.5805                   3.1151  \n",
       "\n",
       "[69 rows x 40 columns]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load data, not from Drive for Desktop since PhotoZ complains about loading zda from Drive\n",
    "data_dir = \"C:/Users/jjudge3/Desktop/Data/mm_full_pipeline_targets/\"\n",
    "# data_dir = 'C:/Users/jjudg/Downloads/'\n",
    "cell_type_prefix = 'scnn1a'\n",
    "csv_prefix = cell_type_prefix + '_all_slices-final.csv'\n",
    "\n",
    "# load the all-data CSV dir\n",
    "df = pd.read_csv(data_dir + csv_prefix)\n",
    "df = pd.DataFrame(df)\n",
    "df = df[df[\"Exclude Reason\"].isnull()]\n",
    "\n",
    "# filter to only rows with a numerical Medial entry\n",
    "df = df[pd.to_numeric(df['Most Medial Barrel'], errors='coerce').notnull()]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "84586b97",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Barrel pairs in L2/3\n",
      "Medial\n",
      "\t 9 crossing out of 12 in layer L2/3\n",
      "Lateral\n",
      "\t 8 crossing out of 12 in layer L2/3\n",
      "\n",
      "Barrel pairs in L4\n",
      "Medial\n",
      "\t 10 crossing out of 21 in layer L4\n",
      "Lateral\n",
      "\t 2 crossing out of 16 in layer L4\n",
      "\n",
      "Barrel pairs in L5\n",
      "Medial\n",
      "\t 1 crossing out of 2 in layer L5\n",
      "Lateral\n",
      "\t 0 crossing out of 2 in layer L5\n"
     ]
    }
   ],
   "source": [
    "# figure out if each barrel crossing is in the medial direction\n",
    "data = {}\n",
    "for stim_loc in ['L2/3', 'L4', 'L5']:\n",
    "    df2 = df[df['Stim Location'].str.contains(stim_loc)]\n",
    "    print(\"\\nBarrel pairs in\", stim_loc)\n",
    "    \n",
    "    df_medial3 = df2[df2[\"Most Medial Barrel\"] == 3]\n",
    "    df_medial1 = df2[df2[\"Most Medial Barrel\"] == 1]\n",
    "    # if most medial barrel is 3, then:\n",
    "    #    - if stim barrel is 3, crossing 2/3 is lateral.\n",
    "    #    - if stim barrel is 2, crossing 2/3 is medial, crossing 1/2 is lateral\n",
    "    #    - if stim barrel is 1, crossing 1/2 is medial\n",
    "    # if most medial barrel is 1, then:\n",
    "    #    - if stim barrel is 1, crossing 1/2 is lateral.\n",
    "    #    - if stim barrel is 2, crossing 1/2 is medial, crossing 2/3 is lateral\n",
    "    #    - if stim barrel is 3, crossing 2/3 is medial\n",
    "    \n",
    "    k_medial_cross = len(df_medial1[(df_medial1['# crossings 2/3'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 3)])\n",
    "    k_lateral_cross = len(df_medial3[(df_medial3['# crossings 2/3'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 3)])\n",
    "    k_medial_cross += len(df_medial3[(df_medial3['# crossings 2/3'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)])\n",
    "    k_lateral_cross += len(df_medial3[(df_medial3['# crossings 1/2'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)])\n",
    "    k_medial_cross += len(df_medial3[(df_medial3['# crossings 1/2'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 1)])\n",
    "    k_lateral_cross += len(df_medial1[(df_medial1['# crossings 1/2'] > 0)\n",
    "                                   & (df_medial1['Stim Barrel'] == 1)])\n",
    "    k_medial_cross += len(df_medial1[(df_medial1['# crossings 1/2'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)])\n",
    "    k_lateral_cross += len(df_medial1[(df_medial1['# crossings 2/3'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)])\n",
    "    \n",
    "    n_medial = k_medial_cross\n",
    "    n_lateral = k_lateral_cross\n",
    "    n_medial += len(df_medial1[(df_medial1['# no crossings 2/3'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 3)])\n",
    "    n_lateral += len(df_medial3[(df_medial3['# no crossings 2/3'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 3)])\n",
    "    n_medial += len(df_medial3[(df_medial3['# no crossings 2/3'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)])\n",
    "    n_lateral += len(df_medial3[(df_medial3['# no crossings 1/2'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)])\n",
    "    n_medial += len(df_medial3[(df_medial3['# no crossings 1/2'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 1)])\n",
    "    n_lateral += len(df_medial1[(df_medial1['# no crossings 1/2'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 1)])\n",
    "    n_medial += len(df_medial1[(df_medial1['# no crossings 1/2'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)])\n",
    "    n_lateral += len(df_medial1[(df_medial1['# no crossings 2/3'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)])\n",
    "    \n",
    "    print(\"Medial\")\n",
    "    print(\"\\t\", k_medial_cross, \"crossing out of\", n_medial, \"in layer\", stim_loc)\n",
    "\n",
    "    print(\"Lateral\")\n",
    "    print(\"\\t\", k_lateral_cross, \"crossing out of\", n_lateral, \"in layer\", stim_loc)\n",
    "    \n",
    "    data[stim_loc] = {\"Medial\": {}, \"Lateral\":{} }\n",
    "    data[stim_loc][\"Medial\"]['n_barrel_cross'] = k_medial_cross\n",
    "    data[stim_loc][\"Medial\"]['total'] = n_medial\n",
    "    data[stim_loc][\"Lateral\"]['n_barrel_cross'] = k_lateral_cross\n",
    "    data[stim_loc][\"Lateral\"]['total'] = n_lateral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "db8a7042",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0 12\n",
      "3.1622776601683795 21\n",
      "2.8284271247461903 12\n",
      "1.4142135623730951 16\n",
      "{'# crossing medial': [25.0, 15.058465048420855], '# crossing lateral': [23.570226039551585, 8.838834764831844]}\n",
      "{'# crossing medial': [75.0, 47.61904761904762], '# crossing lateral': [66.66666666666667, 12.5]}\n",
      "Binomial test pval for dir medial, layer L2/3 : 0.14599609375000003\n",
      "Binomial test pval for dir medial, layer L4 : 1.0\n",
      "Binomial test pval for dir lateral, layer L2/3 : 0.3876953125000001\n",
      "Binomial test pval for dir lateral, layer L4 : 0.004180908203124999\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEiCAYAAADjxEWuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtR0lEQVR4nO3deZwU1bn/8c8XRAFBRXHBjVGDC7uIW+JPxiVqYlSMxmhiBDdyMW65xmi8iZpEc40ajRqXuEYTo3CNGqNxRQlRUVkEQcGQKAIKCLiLqMDz+6POYDP09DQz09M9zPf9evVrutbzVHfXPFXnVJ1SRGBmZlZbm3IHYGZmlckJwszM8nKCMDOzvJwgzMwsLycIMzPLywnCzMzycoKwvCTtIOlFSR9KOn01l/1/kl4tVWw55VRJCklrlbqs5ibpBkk/K3ccDVGq70XSw5KGNMF61tjfTVNrlQlC0kxJn0j6KOe1ebnjqjA/BkZHROeIuLr2REm9JD0m6V1J70maIOnrABHxz4jYodkjXjXG3O95vqTbJHUqYrnRkk6qNS4kfal00a4sIv4rIn7ZkGVT/EvSdr8vaYykPk0dY0M19HuJiK9FxO3NEWOxJF0o6U/ljqNUWmWCSA6JiE45r7dyJ/rogu7AywWm/w14HNgU2AQ4HfigGeJaXYdERCdgALAr8NPmDqBMv6VT03ZvBIwG/tiQleSLvYm2p0m/l9a6v5Z6u1tzglhFOkr8gaQZwIw07huSJqWj5Gcl9c2Zf2dJE1M1zAhJd0u6KE0bKunpPOv/Unq/jqTLJc1KR1E3SOqQplVLmiPpLElvS5or6fic9XSQ9BtJb6QjxKfTuIcknVarzJckDa5jew+V9HLattGSdkrjnwT2AX6XjvK2r7VcV2Ab4KaI+Cy9nomIp3Pjz5l/QE511f+lz+qi3HkLbOvBadkPJM2WdGFRX2YtEfEm8DDQW1IXSQ9KWpDOgB6UtGUq72Lg/+Vs++8kjUmrmZzGfTvNW+i3MVPSOZJeAj6W9KX0/Q9J3/lCSf9TV7yS/lDsZ1TPdi8F7gZ65qx7N0ljU9xz0zaunTN9pf0gp/xzJM0DbpPURtK5kv4jaZGkkZI2LCamWvEV9b2kuFac2aX96xlJV0p6B7gwfcb/SPvEQkkj6in+BElvpc/grLTezSQtlrRRTrm7pJjarc625Xw+H0p6RdLhafw6kt5RzlmdpE2UnVVtnIZX57dVuiQREa3uBcwE9s8zPsiOijcEOpAd3bwN7A60BYakZdcB1gbeAH4ItAOOBD4HLkrrGgo8nWf9X0rvfws8kMrqTHZE/r9pWjWwFPhFWvfXgcVAlzT9WrKjwi1SXF9OMR0FPJ9TXj9gEbB2nm3dHvgY+Goq48fAv2vmTes/qY7PT2QJ9EFgMLBprenVwJz0vuZzOiOV803gs5zPqb5trQb6kB3M9AXmA4PTtKr0ma5V3/cMbEV2RvRLsqPqI4CO6bP/P+D+nOVW2fbc7y4N1/nbyCl7Uiq3Q06sN6XhfsCnwE51xP6HYj+jPMuuiD99/hcDY3Km7wLsAayV4poGnFlgP6gp/9dkv7MOwJnAc8CWadzvgbua63sh27+WAqel7egA3AX8D9lvpT2wVx3l18R3F7Au2e9rQU5MfweG58x/JXBNHeu6EPhTHdO+BWye4vk22f7WLU27Dvh1zrxnAH9ryG+rpP8rS7nySn2lD/gj4L30uj9nx9g3Z77rgV/WWvZVYBCwN/AWoJxpz1JEgiD7B/sxsF3OtD2B19P7auCT3B0s/WD2SD+2T4B+ebZrHeAdoEcavhy4ro7P4GfAyJzhNsCbQHUaXrEz1rH8lsDvgP8Ay4ExOeVW80WC2DutN/dzepqV//nl3dY6yv0tcGV6X0X9/4hqvuc3yHbKVXYooD/wbs7wKtvOqgmizt9GTtkn5EyriXXLnHEvAEfXEfsfGvEZjSZLIO+RJeP3gf0KfJdnAvfV2tbc/aA6rad9zrhpuesEupEdINUknZJ+L2T716xa898B3Jj7GddRfk18O+aMuxS4Jb3/NvBMet8WmAfsVse6LqSOBJFn3knAYen97sBsoE0aHg8c1ZDfVilfrbmKaXBEbJBeg3PGz8553x04K53mvSfpPbKsvXl6vRnpG0veKLLsjcmOkibkrPeRNL7GosiqB2osBjoBXcmOjv5Te6UR8SkwEjhWUhvgGOque948N96IWE627VsUswERMSciTo2I7cg+p4/JdtB85dT+nGbXmqeubUXS7pKeSqf47wP/RfYZFKvme+4eEadExCeSOkr6vbIqug/IktsGktquxnoL/Tbq2k7I/tmssp1FqPMzqsPpEbEB2W/lG8A9NdUUkrZP1Tfz0vb/ilU/09qxL4iIJTnD3YH7crZ9GrCMrE2qGE3xvdSO8cdkB18vKKs6PaGeGHKXf4Mvvru/Aj0lbUt2hv1+RLxQ5HatIOm4nGqi94DepM85Ip4n22cGSdqR7MDxgbRoQ39bTa41J4i61P5HdnFOItkgIjpGxF3AXGALScqZf+uc9x+TJQEgq9vMmbaQ7IiwV85614+s0a4+C4ElwHZ1TL8d+C6wH7A4IsbWMd9bZD/EmvhE9iN8s4gYVhIRs8mqvXrnmZzvc9pqNVb/Z7IdZ6uIWB+4geyfQGOcBewA7B4R65Gd5ZCz3si71MoK/TZqFLOekoqI5RHxT7LqwwPS6OuB6WRnfOsB57HqZ1o79trDs4Gv1dr+9pG1KTRUfd9LbSvFFBHzIuLkiNgc+D5wnQpfeZb7O9yabJ8gJcKRZPvR92hAA7+k7mTViacCG6VkPbXWttwOHJvKuCcnAVfMb8sJorCbgP9KR7GStK6yRtPOwFiyOtDTJa0l6ZvAbjnLTgZ6SeovqT3ZqSiw4mj9JuBKSZsASNpC0oH1BZSWvRW4QtLmktpK2lPSOmn6WLIqn99Q+Ic9EjhY0n6p8e0ssjrxZ+uLITUm/jw1CrZR1mh9AlmddG1jyY4sT02f02Gs/DnVpzPwTkQskbQb8J3VWLbQOj8B3lPWsHpBrenzgW3rGVfot1FRJO1J1khdc1VaZ7Irzj5KR6/DG7DaG4CL0z9CJG2cvtvGqO97KUjSt/RFo/a7ZP9ElxVY5GfprKUXcDyQ26h9B1k11qFAfZextpHUPue1DlnbRpC1baDsooLaB1B/BA4nSxK5Z98V89tygiggIsYDJ5PVtb9LdhQ2NE37jKzBdWia9m3g3pxl/0XWqPgEWYPuSlc0Aeek9T2XTqefIDt6KsaPgCnAOLI2h1+z8nd5B1nDW50/7Ih4leyHeQ3ZWckhZJceflZE+Z+R1eM+QfaPZipZchmap5yaz+lEsjrnY8katz8tohyAU4BfSPoQOJ8ssTXWb8kaNReSJbVHak2/CjhS2ZU0NfeAXAjcnk75jyr026gQNVdhfUT2j+inEfFwmvYjskT7Idk/o/qu9snnKrIzu8fSd/McWb16Y/yWwt9LfXYFnk/b/ABwRkS8XmD+f5B9b6OAyyPisZoJEfEM2YHWxIiYWU+5x5AltprXfyLiFbKDtLFkBxd9gGdyF4qIOcBEskTyz5zxFfPb0spVw9YYkv5A1jjb7Nfa14rjOGBYROxVzjjqIul54IaIuK3csZjVRdnl3n+OiJtLWMatwFvl/p9Rl1Z5c8maTFJHsqPu68odSw1Jg8iuwlhIVq/bl9U/OjRrNpJ2JbvctLHVZoXKqCI7u965VGU0lquY1iCpDWMB2Sntn8scTq4dyNpk3idr6zgyIuaWNySz/CTdTlZ9emZEfFiiMn5JVjV7WT3VYGXlKiYzM8vLZxBmZpaXE4SZmeXVohupu3btGlVVVeUOw8ysRZkwYcLCiNi4vvladIKoqqpi/Pjx5Q7DzKxFkVRUt0CuYjIzs7ycIMzMLC8nCDMzy6tFt0GYWdP5/PPPmTNnDkuWLKl/ZmvRHn/88T6TJ09+DZi6dOnSk3bZZZe3883nBGFmAMyZM4fOnTtTVVXFyr2z25pm2bJlS3v27PnuggULes6bN+9msl5rV1GyKiZJtyp7fu7UnHEbSnpc0oz0t0vOtJ9I+rekV4vp9trMmtaSJUvYaKONnBxaiTZt2sTGG2/8Pvmf45LNU8Ly/wAcVGvcucCoiOhB1sXuuQCSegJHA73SMtcVeIqUmZWIk0Pr0qZNm6BAHihZgoiIMWTPKsh1GNlTlEh/B+eMvzsiPk0dV/2b1XuojJmtQX7yk58wevRo7r//fi655JKyxnL++efzxBNPlDWGfIYOHco999wDwEknncQrr7xScP7q6urVvm+sudsgNq3pxTMi5tY8TY3sOci5TyObQx3PRpY0DBgGsPXWW+ebpWyqq6sBGD16dFnjMGsKVec+1KTrm3nJwUXP+/zzz3P++edz3nnnceSRRzaovIggImjTpnHHwb/4xS8atXxzuPnm0jyyolIuc813Xpu3m9mIuDEiBkbEwI03rvdOcTNrQc4++2z69u3LuHHj2HPPPbn55psZPnx43n/S8+fP5/DDD6dfv37069ePZ599lpkzZ7LTTjtxyimnMGDAAGbPns3ZZ59N79696dOnDyNGZA/Pmzt3LnvvvTf9+/end+/e/POf/2TZsmUMHTp0xbxXXnklsPKRelVVFRdccAEDBgygT58+TJ8+HYAFCxbw1a9+lQEDBvD973+f7t27s3DhwlVi7tSpE+eccw677LIL+++/Py+88ALV1dVsu+22PPDAAwAsW7aMs88+m1133ZW+ffvy+9//HsgS3qmnnkrPnj05+OCDefvtLy48yj07GD58OAMHDqRXr15ccMFqPbV1Fc2dIOZL6gaQ/tZs4RxWfoD4lqQHiJtZ63HZZZdx8803M3ToUMaNG0ffvn156aWXOP/881eZ9/TTT2fQoEFMnjyZiRMn0qtXLwBeffVVjjvuOF588UXGjx/PpEmTmDx5Mk888QRnn302c+fO5c9//jMHHnjgimn9+/dn0qRJvPnmm0ydOpUpU6Zw/PHH542xa9euTJw4keHDh3P55ZcD8POf/5x9992XiRMncvjhhzNr1qy8y3788cdUV1czYcIEOnfuzE9/+lMef/xx7rvvvhXbeMstt7D++uszbtw4xo0bx0033cTrr7/Offfdx6uvvsqUKVO46aabePbZ/I+Pv/jiixk/fjwvvfQS//jHP3jppZdW+3uo0dwJ4gFgSHo/BPhrzvijJa0jaRugB/BCM8dmZhXgxRdfpH///kyfPp2ePXvWOd+TTz7J8OHDAWjbti3rr78+AN27d2ePPfYA4Omnn+aYY46hbdu2bLrppgwaNIhx48ax6667ctttt3HhhRcyZcoUOnfuzLbbbstrr73GaaedxiOPPMJ6662Xt9xvfvObAOyyyy7MnDlzRTlHH300AAcddBBdunTJu+zaa6/NQQdl1+706dOHQYMG0a5dO/r06bNiXY899hh33HEH/fv3Z/fdd2fRokXMmDGDMWPGrNiWzTffnH333TdvGSNHjmTAgAHsvPPOvPzyy/W2TRRSsjYISXcB1UBXSXOAC4BLgJGSTgRmAd8CiIiXJY0EXgGWAj+IiGWlis3MKs+kSZMYOnQoc+bMoWvXrixevJiIoH///owdO5YOHToUtZ511113xfu6Hoi29957M2bMGB566CG+973vcfbZZ3PccccxefJkHn30Ua699lpGjhzJrbfeusqy66yzDpAlpaVLlxYsp7Z27dqtuFKsTZs2K9bVpk2bldZ1zTXXcOCBK1/t//e//73eq8xef/11Lr/8csaNG0eXLl0YOnRoo258LOVVTMdERLeIaBcRW0bELRGxKCL2i4ge6e87OfNfHBHbRcQOEfFwqeIys8pUU82z/fbb88orr7Dvvvvy6KOPMmnSpLzJYb/99uP6668Hsnr7Dz74YJV59t57b0aMGMGyZctYsGABY8aMYbfdduONN95gk0024eSTT+bEE09k4sSJLFy4kOXLl3PEEUfwy1/+kokTJxYd+1577cXIkSOB7Azg3XffbeCnAAceeCDXX389n3/+OQD/+te/+Pjjj9l77725++67WbZsGXPnzuWpp55aZdkPPviAddddl/XXX5/58+fz8MON+1fqO6nNrGIsWLCALl260KZNm3qrmK666iqGDRvGLbfcQtu2bbn++uvp1q3bSvMcfvjhjB07ln79+iGJSy+9lM0224zbb7+dyy67jHbt2tGpUyfuuOMO3nzzTY4//niWL18OwP/+7/8WHfcFF1zAMcccw4gRIxg0aBDdunWjc+fODfoMTjrpJGbOnMmAAQOICDbeeGPuv/9+Dj/8cJ588kn69OnD9ttvz6BBg1ZZtl+/fuy888706tWLbbfdlq985SsNiqFGi34m9cCBA6OSngfhy1ytJZs2bRo77bRTucNokT799FPatm3LWmutxdixYxk+fDiTJk0qd1h1mjp16uLevXtPA5g8eXLXfv36VeWbz2cQZmaNNGvWLI466iiWL1/O2muvzU033VTukJqEE4SZWSP16NGDF198sdxhNLlKuVHOzMwqjBOEmZnl5QRhZmZ5OUGYmVleThBmVnHW1O6+q6qq8nbil+tXv/pVk5RVI7ezwdXlq5jMLK9pOzbtPRE7TZ9W9LytubvvX/3qV5x33nmrtcyyZcto27bpn7HmMwgzqxhrenffuQYPHswuu+xCr169uPHGGwE499xz+eSTT+jfvz/f/e53AfjTn/7EbrvtRv/+/fn+97/PsmVZN3WdOnXi/PPPZ/fdd2fs2LH84he/YNddd6V3794MGzas6P6hCnGCMLOKsaZ3953r1ltvZcKECYwfP56rr76aRYsWcckll9ChQwcmTZrEnXfeybRp0xgxYgTPPPMMkyZNom3bttx5551A1nV47969ef7559lrr7049dRTGTduHFOnTuWTTz7hwQcfbOjXsIKrmMysoqxOd9933HEH8EV33++++27R3X2fcMIJfP755wwePJj+/fuv1N33wQcfzAEHHJC33Nzuvu+9994V5dx3331A4e6+c1199dUrlpk9ezYzZsxgo402WmmeUaNGMWHCBHbddVcAPvnkEzbZZJMV23zEEUesmPepp57i0ksvZfHixbzzzjv06tWLQw45pN44CnGCMLOK0Bq6+64xevRonnjiCcaOHUvHjh2prq7O2y13RDBkyJC8HQe2b99+RbvDkiVLOOWUUxg/fjxbbbUVF154YaO6+a7hKqZWrLq6ekUHg2bl1pq6+37//ffp0qULHTt2ZPr06Tz33HMrprVr125FV9/77bcf99xzz4rHi77zzju88cYbq6yvJhl07dqVjz76qMFXLdXmMwgzqxitpbvvgw46iBtuuIG+ffuyww47rKgSAxg2bBh9+/ZlwIAB3HnnnVx00UUccMABLF++nHbt2nHttdfSvXv3lda3wQYbcPLJJ9OnTx+qqqpWVEk1lrv7bkItrbvvlhavlZa7+244d/dtZmZ5ubtvMzPLy919m5lZq+IEYWYrtOQ2SVt9y5cvF7C8rulOEGYGZNfVL1q0yEmilVi+fLkWLFiwPjC1rnncBmFmAGy55ZbMmTOHBQsWlDsUK7F58+attWzZsi7A1KVLl55U13xOEGYGZDdobbPNNuUOw5pBz549p0TEwPrmcxWTmZnl5QRhZmZ5OUGYmVleThBmZpaXE4SZmeXlBGFmZnnVmyAkfUXSuun9sZKukNS9vuXMzKxlK+YM4npgsaR+wI+BN4A7ShqVmZmVXTEJYmlk994fBlwVEVcBdT8JowiSfijpZUlTJd0lqb2kDSU9LmlG+lv/Q13NzKxkikkQH0r6CXAs8JCktkC7hhYoaQvgdGBgRPQG2gJHA+cCoyKiBzAqDZuZWZkUkyC+DXwKnBgR84AtgMsaWe5aQAdJawEdgbfIzlBuT9NvBwY3sgwzM2uEevtiSknhipzhWTSiDSIi3pR0OTAL+AR4LCIek7RpRMxN88yVtElDyzAzs8Yr5iqmDyV9UOs1W9J9krZd3QJT28JhwDbA5sC6ko5djeWHSRovabx7nTQzK51ienO9gqwK6M+AyNoLNgNeBW4FqlezzP2B1yNiAYCke4EvA/MldUtnD92At/MtHBE3AjcCDBw40B3Xm5mVSDFtEAdFxO8j4sOI+CD9g/56RIwAGnKl0SxgD0kdJQnYD5gGPAAMSfMMAf7agHWbmVkTKeYMYrmko4B70vCROdNW+wg+Ip6XdA8wEVgKvEh2RtAJGCnpRLIk8q3VXbeZmTWdYhLEd4GrgOvIEsJzwLGSOgCnNqTQiLgAuKDW6E/JzibMzKwCFHMV02vAIXVMfrppwzEzs0pRb4KQtDFwMlCVO39EnFC6sMzMrNyKqWL6K/BP4AlgWWnDMTOzSlFMgugYEeeUPBIzM6soxVzm+qCkr5c8EjMzqyjFnEGcAZwn6VPgc7Kb5SIi1itpZM2g6tyHmnR9815bVJL1Asy85OAmX6eZWSHFXMXUqK69zcysZaozQUjaMSKmSxqQb3pETCxdWGZmVm6FziD+GxgG/CbPtAD2LUlEZmZWEepMEBExLP3dp/nCMTOzSlFMd9/fktQ5vf+ppHsl7Vz60MzMrJyKucz1ZxHxoaS9gAPJnvZ2Q2nDMjOzcismQdTcPX0wcH1E/BVYu3QhmZlZJSgmQbwp6ffAUcDfJa1T5HJmZtaCFfOP/ijgUbIHB70HbAicXcqgzMys/Iq5k7ob8FBEfCqpGugL3FHKoMzMrPyKOYP4C7BM0peAW4BtyJ5PbWZma7BiEsTyiFgKfBP4bUT8kOyswszM1mDFJIjPJR0DHAc8mMa1K11IZmZWCYpJEMcDewIXR8TrkrYB/lTasMzMrNzqTRAR8QrwI2CKpN7AnIi4pOSRmZlZWRXzTOpqsrunZ5I9C2IrSUMiYkxJIzMzs7Iq5jLX3wAHRMSrAJK2B+4CdillYLayaTvu1OTrXDzrjZKte6fp05p8nWbWvIppg2hXkxwAIuJfuJHazGyNV8wZxARJtwB/TMPfBSaULiQzM6sExSSI/wJ+AJxO1gYxBriulEGZmVn5FUwQktoAEyKiN3BF84RkZmaVoGAbREQsByZL2rqZ4jEzswpRbGd9L0t6Afi4ZmREHFqyqMzMrOyKSRA/L3kUZmZWcepMEKn31k0j4h+1xu8NvFnqwMzM8qmurgZg9OjRZY2jNSjUBvFb4MM84xenaWZmtgYrlCCqIuKl2iMjYjxQ1ZhCJW0g6R5J0yVNk7SnpA0lPS5pRvrbpTFlmJlZ4xRKEO0LTOvQyHKvAh6JiB2BfsA04FxgVET0AEalYTMzK5NCCWKcpJNrj5R0Io24k1rSesDeZE+nIyI+S8+6PoysU0DS38ENLcPMzBqv0FVMZwL3ScrtWmMgsDZweCPK3BZYANwmqV9a9xlkDeJzASJirqRNGlGGmZk1Up0JIiLmA1+WtA/QO41+KCKebIIyBwCnRcTzkq5iNaqTJA0DhgFsvbXv3zMzK5V674OIiKeAp5qwzDlkDx16Pg3fQ5Yg5kvqls4eugFv1xHPjcCNAAMHDowmjMvMzHIU0913k4qIecBsSTukUfsBrwAPAEPSuCHAX5s7NjMz+0Ixd1KXwmnAnZLWBl4je+51G2BkagSfBXyrTLGZmRn19+baFng0IvZvykIjYhJZg3dt+zVlOWZm1nD19ea6DFgsaf1misfMzCpEMVVMS4Apkh5n5d5cTy9ZVGZmVnbFJIiH0svMzFqRYi5zvb2+eczMbM1TqLvvkRFxlKQpwCr3G0RE35JGZmZmZVXoDOKM9PcbzRGImZlVlkJdbdT0i/RG84VjZmaVot47qSXtIWmcpI8kfSZpmaQPmiM4MzMrn2K62vgdcAwwg+w5ECcB15QyKDMzK7+iutqIiH9LaptunLtN0rMljsvMzMqsmASxOPWZNEnSpcBcYN3ShmVmZuVWTBXTsWm+U8nupN4KOKKUQZmZWfkVug+iB3A5sB0wBfhRRPy8uQIzM7PyKnQGcSvwINnZwkTcMG1m1qoUaoPoHBE3pfeXSZrYHAGZmVllKJQg2kvaGVAa7pA7HBFOGGZma7BCCWIucEXO8Lyc4QD2LVVQLdVm37mk3CGYmTWZQl1t7NOcgZiZWWUp5jJXMzNrhZwgzMwsLycIMzPLq9CNcgMKLeirmMysPlXnNv3Tiue9tqgk6555ycFNur41QaGrmH5TYJqvYjIzW8P5KiYzM8urmAcGdZT0U0k3puEekvwYUjOzNVwxjdS3AZ8BX07Dc4CLShaRmZlVhGISxHYRcSnwOUBEfMIX3W+YmdkaqpgE8ZmkDmQN00jaDvi0pFGZmVnZFfNEuQuAR4CtJN0JfAUYWsqgzMys/AomCEltgC7AN4E9yKqWzoiIhc0Qm5mZlVHBBBERyyWdGhEjgaa/48XK6vatu5c7BDOrYMW0QTwu6UeStpK0Yc2r5JGZmVlZFdMGcUL6+4OccQFs2/ThmJlZpSimDeLciBjR1AVLaguMB96MiG+ks5IRQBUwEzgqIt5t6nLNzKw4BauYImI5K585NKUzgGk5w+cCoyKiBzAqDZuZWZmUpQ1C0pbAwcDNOaMPA25P728HBjemDDMza5xytUH8Fvgx0Dln3KYRMRcgIuZK2iTfgpKGAcMAtt5660aEYGZmhdSbICJim6YsMHX093ZETJBUvbrLR8SNwI0AAwcOjKaMzczMvlDMGQSSegM9gfY14yLijgaW+RXgUElfT+tbT9KfgPmSuqWzh27A2w1cv5mZNYFiuvu+ALgmvfYBLgUObWiBEfGTiNgyIqqAo4EnI+JY4AFgSJptCPDXhpZhZmaNV0wj9ZHAfsC8iDge6AesU4JYLgG+KmkG8NU0bGZmZVJMFdMnqcuNpZLWI6v6aZKb5CJiNDA6vV9ElojMzKwCFJMgxkvaALgJmAB8BLxQyqDMzKz8irmK6ZT09gZJjwDrRcRLpQ3LzMzKrWAbhKS1JCm93woYCLRtjsDMzKy86kwQkk4ma294I70fRdZgfbekc5opPjMzK5NCVUxnAtuR3e08DegeEQsldQTGAb8ufXhmZlYuhRLEZ6k31Xcl/bvmKXIRsVjSZ80TnpmZlUuhBNFB0s5k1VBrp/dKr/YFljMzszVAoQQxF7givZ+X875m2MzM1mB1JoiI2Kc5AzEzs8pSTFcbZmbWChXVm6uZWaXY7Dvupq25FLoP4ivpbyk65jMzswpXqIrp6vR3bHMEYmZmlaVQFdPnkm4DtpB0de2JEXF66cIyM7NyK5QgvgHsD+xL1ourmZm1IoUuc11I1u/StIiY3IwxmZlZBSjmMtdFku6T9Lak+ZL+ImnLkkdmZmZlVUyCuI3sedGbA1sAf0vjzMxsDVZMgtgkIm6LiKXp9Qdg4xLHZWZmZVZMglgg6VhJbdPrWGBRqQMzM7PyKiZBnAAcRdZB31yyhwadUMqgzMys/Ip5JvUs4NBmiMXMzCqIO+szM7O8nCDMzCwvJwgzM8ur6AQhaQ9JT0p6RtLgEsZkZmYVoM5GakmbRUTuo0X/m6yxWsCzwP2lDc3MzMqp0FVMN0iaAFwWEUuA94DvAMuBD5ohNjMzK6M6q5giYjAwCXhQ0veAM8mSQ0dgcOlDMzOzcirYBhERfwMOBDYA7gVejYirI2JBM8RmZmZlVOiRo4dKehp4EpgKHA0cLukuSds1V4BmZlYehdogLgL2BDoAf4+I3YD/ltQDuJgsYZiZ2RqqUBXT+2RJ4Gjg7ZqRETEjIhqcHCRtJekpSdMkvSzpjDR+Q0mPS5qR/nZpaBlmZtZ4hRLE4WQN0kvJrl5qKkuBsyJiJ2AP4AeSegLnAqMiogcwKg2bmbVY1dXVVFdXlzuMBqvvkaPXNHWBETGXrFdYIuJDSdPIHkR0GFCdZrsdGA2c09Tlm5lZccra1YakKmBn4Hlg05Q8apLIJmUMzcys1StbgpDUCfgLcGZEFH3jnaRhksZLGr9gga+2NTMrlbIkCEntyJLDnRFxbxo9X1K3NL0bOQ3juSLixogYGBEDN97YTz41MyuVZk8QkgTcAkyLiCtyJj0ADEnvhwB/be7YzMzsC/U+Ua4EvgJ8D5giaVIadx5wCTBS0onALOBbZYjNzMySZk8QEfE0WY+w+ezXnLGYmVnd/MAgMzPLywnCzMzycoIwM7O8nCDMzCwvJwgzM8vLCcLMzPJygjAzs7ycIMzMLC8nCDMzy8sJwszM8ipHX0xmZhVn2o47Nfk6F896o2Tr3mn6tCZfZ20+gzAzs7ycIMzMLC8nCDMzy8sJwszM8nKCMDOzvJwgzMwsLycIMzPLywnCzMzycoIwM7O8nCDMzCwvJwgzM8vLCcLMzPJygjAzs7ycIMzMLC8nCDMzy8sJwszM8nKCMDOzvJwgzMwsLz9y1MysRG7funu5Q2gUn0GYmVleThBmZpaXE4SZmeVVcQlC0kGSXpX0b0nnljseM7PWqqIShKS2wLXA14CewDGSepY3KjOz1qmiEgSwG/DviHgtIj4D7gYOK3NMZmatUqVd5roFMDtneA6we+4MkoYBw9LgR5JebabYyirPaVRXYGGzB1IsqdwRmK2WOqoqKnc/a9w+VtT1t5WWIPJtcaw0EHEjcGPzhFO5JI2PiIHljsNsTdba97NKq2KaA2yVM7wl8FaZYjEza9UqLUGMA3pI2kbS2sDRwANljsnMrFWqqCqmiFgq6VTgUaAtcGtEvFzmsCpVq69mM2sGrXo/U0TUP5eZmbU6lVbFZGZmFcIJwszM8nKCqECSPsoz7r8lvSLpJUmjJHWvNf0RSVtIukXS5DTfPZI6NV/kZi1Xvv0uZ9qRkkJSq7rk1Qmi5XgRGBgRfYF7gEtrJkjqAGwYEW8CP4yIfmm+WcCpZYnWbA0hqTNwOvB8uWNpbk4QLUREPBURi9Pgc2T3iNSoBkan+T4AkCSgA7VuNDSz1fZLsgOyJeUOpLk5QbRMJwIP5wx/DXikZkDSbcA8YEfgmuYNzWzNIWlnYKuIeLDcsZSDE0QLI+lYYCBwWc7orwBP1wxExPHA5sA04NvNGqDZGkJSG+BK4Kxyx1IuThAtiKT9gf8BDo2IT9O4bYHZqffbFSJiGTACOKLZAzVbM3QGegOjJc0E9gAeaE0N1U4QLUQ61f09WXJ4O2fSiuolZb5U8x44BJje3LGarQki4v2I6BoRVRFRRdb2d2hEjC9zaM2morrasBU6SpqTM3wF8HWgE/B/2f9+ZkXEocBBwGlpPgG3S1ovvZ8MDG+2qM1atlX2u4i4omzRVAB3tdGCSVoHeKY1d0dsZqXjBGFmZnm5DcLMzPJygigDSRtJekrSR5J+lzO+o6SHJE2X9LKkS5oxph0ljZX0qaQf5YzfKsU6LcV0RnPFZNYYkr4qaYKkKenvvjnTLpY0u1D3GiWKKe9+lqZtkLrHmZ72tz2bM7Z83EhdHkuAn5FdQte71rTLI+Kp9MCkUZK+FhEPr7KGpvcOWXcCg2uNXwqcFRETU5cDEyQ9HhGvNENMZo2xEDgkIt6S1JvsOTNbpGl/A34HzGjmmOrazwCuAh6JiCPT/t+xOQPLx2cQDSSpKmX5m9KR9WOpT6R6RcTHEfE0tW7dj4jFEfFUev8ZMJGVu9SoL6ZqSaNzjkLuTJe7FhPT2xExDvi81vi5ETExvf+Q7Oa7LfKswqzJNXI/ezEiah5Z/DLQPl3YQUQ8FxFzGxhTk+9n6crDvYFb0nyfRcR7DYmvKTlBNE4P4NqI6AW8R7opTdLZkibleV1d7IolbUB2H8Oo1YxpZ+BMoCewLdld1ki6so6Yzl2NmKrS+ltdp2VWVk2xnx0BvFhzg2kTaOr9bFtgAXCbpBcl3Sxp3SaKtcFcxdQ4r0fEpPR+AlAFEBGXsXJXGKtF0lrAXcDVEfHaai7+QkTMSeuZlGJ6OiJ+2NB40ro6AX8BzqzpENCsmTRqP5PUC/g1cEATxtTU+9lawADgtIh4XtJVwLlkVdFl4wTROLlHI8vIek9F0tnAd/PMPyYiTi9ivTcCMyLit00Q01oppiuBffLMf3dEFGwMl9SOLDncGRH3NiAms8Zo8H4maUvgPuC4iPhPCWNq7H42B5gTETVn5/eQJYiycoIogcacQUi6CFgfOKmJY2rQkU2qW70FmNba7yq1ylLffpaqaR8CfhIRzzRTTA3azyJiXrqqaoeIeBXYDyj7hSBOEGWirPOv9YC1JQ0mO/39gKwzvunAxNTu9buIuLkZ4tkMGJ9iWi7pTLL61b7A94Ap6VQa4LyI+HupYzJrpFOBLwE/k1RTVXNARLwt6VLgO3zRvcbNEXFhqQOqaz9L1banAXemK5heA44vdTz18Z3UZmaWl69iMjOzvJwgzMwsLycIMzPLywnCzMzycoIwM7O8nCCsRZH0P6lPnpdSFwa7p/FnSuqYM9/f03XwTVFmwR4/Uy+cp+QMby7pniYqe7Ra0TOQrbI4QViLkbo//gYwICL6AvsDs9PkM8np/TIivt6MnZ1tAKxIEBHxVkQc2UxlNwlJbcsdg1UeJwhrSboBC2s6XIuIhakr59OBzYGnJD0F2Y2Ikrqm3kCnp87PpqaeN/eX9IykGZJ2S/NfqJWfgzE1dU5IzrhOkkZJmqjsGQOHpUmXANulM5rLUplT0zLtJd2W5n9R0j5p/FBJ90p6JMVxabEfQlr/P1McEyV9OY3/Y05MpG09VFLbFNe4dOb1/TS9WtmzPv4MTFmtb8Jah4jwy68W8QI6AZOAfwHXAYNyps0EutYeJutEbSnQh+yAaAJwKyDgMOD+NP+FwI9ylp8KVKX3H6W/awHrpfddgX+n9VQBU3OWXTEMnAXclt7vCMwC2gNDye6WXT8NvwFslWebRwMDa43rCLRP73sA49P7QTnbsz7weop5GPDTNH4dsjt5twGqgY+Bbcr93fpVmS+fQViLEREfAbuQ/cNbAIyQNLSIRV+PiCkRsZzsuQCjIiLIjpqrViMEAb+S9BLwBNlzMTatZ5m9gD+m+KeTJYLt07RREfF+RCwh63ene5FxtANukjQF+D+yLlGIiH8AX5K0CXAM8JeIWErWjctxqauU54GNyBILZL2Svl5kudbKuC8ma1EiYhnZUfXo9A9yCPCHehbL7Xlzec7wcr7YB5aycpVr+zzr+S6wMbBLRHye+tPKN1+uQg+SydsjaBF+CMwH+pHFnPvgqT+mOI8GTsiJ4bSIeHSlwKRqsjMIs7x8BmEthqQdJPXIGdWf7Igc4EOgcyNWP5OsP34kDSCrgqltfeDtlBz24Ysj/kJljyF1SS1pe2Br4NVGxFkTx9x0RvQ9ILeB+Q9kDfZExMtp3KPA8NRtO5K2r4SH0Vjl8xmEtSSdgGvS5atLydoAhqVpNwIPS5obEfn646/PX/iiGmYcWTtHbXcCf5M0nqwtZDpARCxKjd5TgYeBa3OWuQ64IZ3tLAWGRsSnKu4JlTUeklTziMqxwHnAXyR9C3iKnLOAiJgvaRpwf87yN5NVpU1M3bcvIP8zkc1W4t5czdYg6V6QKWSXAr9f7nisZXMVk9kaQtL+ZGc11zg5WFPwGYSZmeXlMwgzM8vLCcLMzPJygjAzs7ycIMzMLC8nCDMzy8sJwszM8vr/C5u9ykwW7BUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def significance_bar(start,end,height,displaystring,linewidth = 1.2,markersize = 8,boxpad  =0.3,fontsize = 15,color = 'k'):\n",
    "    # draw a line with downticks at the ends\n",
    "    plt.plot([start,end],[height]*2,'-',color = color,lw=linewidth,marker = TICKDOWN,markeredgewidth=linewidth,markersize = markersize)\n",
    "    # draw the text with a bounding box covering up the line\n",
    "    plt.text(0.5*(start+end),height,displaystring,ha = 'center',va='center',bbox=dict(facecolor='1.', edgecolor='none',boxstyle='Square,pad='+str(boxpad)),size = fontsize)\n",
    "\n",
    "n_L23 = {'medial' : data['L2/3'][\"Medial\"]['total'],\n",
    "         'lateral' : data['L2/3'][\"Lateral\"]['total']}\n",
    "n_L4 = {'medial' : data['L4'][\"Medial\"]['total'],\n",
    "         'lateral' : data['L4'][\"Lateral\"]['total']}\n",
    "n_L5 = {'medial' : data['L5'][\"Medial\"]['total'],\n",
    "        'lateral' : data['L5'][\"Lateral\"]['total']}\n",
    "\n",
    "plt.clf()\n",
    "i = 0\n",
    "tick_labels = []\n",
    "colors = {\"# crossing medial\": 'tab:blue', \n",
    "          \"# crossing lateral\": 'tab:red'}\n",
    "# in format # no crossing, # crossing, # no signal; l2/3, L4, L5\n",
    "\n",
    "length_of_previous_sl = i\n",
    "\n",
    "all_layers = ['L2/3', 'L4']\n",
    "totals = {'medial' : [], 'lateral': [] }\n",
    "for med in ['medial', 'lateral']:\n",
    "    if n_L23[med] > 0:\n",
    "        totals[med].append(n_L23[med])\n",
    "    if n_L4[med] > 0:\n",
    "        totals[med].append(n_L4[med])\n",
    "    if n_L5[med] > 0:\n",
    "        totals[med].append(n_L5[med])\n",
    "\n",
    "\n",
    "# in columns L23, L4, L5\n",
    "crossing_data = {\n",
    "    '# crossing medial': [data[x]['Medial']['n_barrel_cross'] for x in all_layers],\n",
    "    '# crossing lateral': [data[x]['Lateral']['n_barrel_cross'] for x in all_layers]\n",
    "}\n",
    "errbars = {}\n",
    "for k in crossing_data:\n",
    "    errbars[k] = [np.sqrt(x) for x in crossing_data[k]]\n",
    "    med = 'lateral'\n",
    "    if 'medial' in k:\n",
    "        med = 'medial'\n",
    "    for j in range(len(errbars[k])):\n",
    "        print(errbars[k][j], totals[med][j])\n",
    "        errbars[k][j] /= totals[med][j] / 100\n",
    "        crossing_data[k][j] /= totals[med][j] / 100 \n",
    "\n",
    "labels = [x for x in crossing_data.keys()]\n",
    "#tick_labels += [all_layers['medial'][j] + \"\\n(n=\"+str(totals['medial'][j]) + \")\"\n",
    "#                            for j in range(len(all_layers))]\n",
    "#tick_labels += [all_layers['lateral'][j] + \"\\n(n=\"+str(totals['lateral'][j]) + \")\"\n",
    "#                            for j in range(len(all_layers))]\n",
    "\n",
    "tick_labels = ['L2/3\\nn=' + str(n_L23['medial'])\n",
    "                   + \"   n=\" + str(n_L23['lateral']),\n",
    "               'L4\\nn=' + str(n_L4['medial'])\n",
    "                   + \"   n=\" + str(n_L4['lateral'])]\n",
    "#,\n",
    "#              'L5\\nn=' + str(n_L5['medial'])\n",
    "#                   + \"   n=\" + str(n_L5['lateral'])]\n",
    "\n",
    "print(errbars)\n",
    "print(crossing_data)\n",
    "\n",
    "# This is a binomial experiment of n independent trials.\n",
    "# We use the Binomial case of goodness of fit test to a \n",
    "# chi-squared distribution\n",
    "pvals = []\n",
    "for med in ['medial', 'lateral']:\n",
    "    for i in range(2):\n",
    "        layer = ['L2/3', 'L4', 'L5'][i]\n",
    "        n = [n_L23[med], n_L4[med], n_L5[med]][i]\n",
    "        cr = crossing_data['# crossing ' + med][i]\n",
    "        k = n * cr / 100 # number of crosses\n",
    "\n",
    "        null_hyp_p = 0.5\n",
    "\n",
    "        pval = binom_test(k, n, null_hyp_p, alternative='two-sided')\n",
    "        print(\"Binomial test pval for dir \" + med + \", layer\", \n",
    "              layer, \":\", pval)\n",
    "        pvals.append(pval)\n",
    "    \n",
    "##########################\n",
    "for nm in crossing_data:\n",
    "    x = np.arange(len(all_layers))\n",
    "    # else:\n",
    "    plt.bar(x + (i-1) * 0.2, \n",
    "        crossing_data[nm],\n",
    "        0.2,\n",
    "        label=(nm if length_of_previous_sl == 0 else None),\n",
    "        color=colors[nm])\n",
    "    plt.errorbar(x + (i-1) * 0.2, \n",
    "                crossing_data[nm],\n",
    "                yerr=errbars[nm],\n",
    "                fmt='none',\n",
    "                color='black')\n",
    "    i += 1\n",
    "\n",
    "#for i in range(3):\n",
    "    #if pvals[i] < 0.001:\n",
    "        #significance_bar(i-.5, i+0.25, 20, \"P<0.001\", fontsize=10)\n",
    "    \n",
    "plt.title(\"Frequency of Signal Pattern in Barrel Pairs by Layer\")\n",
    "plt.xlabel(\"Stimulation Layer\")\n",
    "plt.ylabel(\"% of Barrel Pair Crossings\")\n",
    "#plt.xticks(rotation='vertical')\n",
    "plt.xticks(np.arange(0.12,len(tick_labels)+0.12), tick_labels)\n",
    "plt.legend(bbox_to_anchor=(0.82, 1.0), loc='upper center')\n",
    "plt.savefig(data_dir + \"results_count-medial_oct2023.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "d5226216",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEuCAYAAACESglMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqD0lEQVR4nO3deZxcVZn/8c83YUlY0umQQMDQhC0SyDSLkUUiJP5wwwE06IyMokAcHJwBF6KzoEMjMqMSVHBhyBBEUCOKEYHBbTABwwCSADaBRFQMTYCERDodlgRC8vz+uLeSm04v1d1VXVW3vu/Xq159665Pn7791Klzzz1XEYGZmeXPkEoHYGZm5eEEb2aWU07wZmY55QRvZpZTTvBmZjnlBG9mllNO8DboJJ0naZWkFyXtUel4ykXSWZIWlmG/L0o6oNT7LSdJ/yXpc5WOo944wVcJScslndSP7RZI+kg5YioHSTsCXwHeFhG7RcRfOi0fLykk7dDLfn7d23rp8pfShLhG0lxJI0vyiwxQ5vd8MX0tl/QvxWybltsTfTze2yXdLekFSasl3SXp1P5F33cR8Q8RcelgHc8STvA22PYChgGP9ncHkj4A9PgBkHF4ROwGHAA0Ai39POZ2x+vtQ6hII9P4zgD+XdI7BrIzSUO7mPde4EfADcA4kr/BvwOndLOPUvxeVgWc4KucpEZJt6e1rvZ0ely67DLgzcA30lrgN9L5h0j6laTnJf1e0t9k9ne9pG9K+p+0Nne/pAMzyw/LbLtK0r9JGivp5WxziqQ3pDHt2EXMO0v6mqRn0tfX0nkTgN+nq62V9Ot+lEcDcDHwmb5sFxHrgFuBQzP7OlvS0rQcnpD00cyyqZJWSPpnSSuBb0tqkXSzpO9KWgecJalB0hxJz0p6WtIXukqyRcR3L8mH3iRJR0u6V9LadL/fkLRTJraQdFA6fb2kqyXdIeklYFqn8hLJN6ZLI+LaiOiIiM0RcVdE/H26zlmS7pH0VUnPAy3p73VD+jd+UtJnJQ1J1z8o/QbQkX4zuqlwrHQfz6XLWiVNysT5hU5le2G67rOSzs7EvIek2yStk/RAWqYLezuGdSEi/KqCF7AcOKmL+XsApwO7ALuT1MRuySxfAHwk835X4CngbJJa7lHAGuCwdPn1wPPA0eny7wE/SJftDjwLXEhSy94dOCZddgdwXuY4XwW+3s3v8nngPmBPYAzwfyQJBmA8EMAO3Wzb2/JvAp/sbb103QAOSqcbgV8Cn88sfxdwICDgROBl4Kh02VTgNeBLwM7AcJLa/0bg3SSVo+HALcA1abnvCfwW+Gi6j7OAhb39nunxj0+P//+ANwDHpsvGA0uBT3Tze10PdKTbDwGGdTrOIen6+/dQTmelv+v56TGHk9T2f5qeA+OBx4EZ6fpzgYsKxwOmpPPfDiwGRqa/00Rg70ycX+hUtp8HdgROTn/3xnT5D9LXLiQfyE8VyrGnY/jVxd+20gH4lf4huknwXax3BNCeeb+AbRP83wK/6bTNNcDF6fT1wLWZZScDy9LpM4CHujnu3wL3pNNDgZXA0d2s+yfg5Mz7twPL0+ktia2bbbtdDkwGHmZr4ismwa8D1gKbgGXA63pY/xbg4+n0VOBVMgmTJMHfnXm/F/AKMDwz7wxgfjp9Fr0n+LVAO0kSv6CbdT8B/KTT75VN8Df08Dsdn64/rId1zgLaMu+Hpr/XoZl5HwUWpNM3ALOBcZ328xaSD4JjgSGdll3Ptgl+ffZvBzyXbjeU5EP09ZllX2Brgu/2GH5t/3ITTZWTtIuka9KvyeuAu4GRPTQD7Acck369XytpLfABYGxmnZWZ6ZeB3dLpfUmSc1d+ChyqpPfGW4GOiPhtN+vuAzyZef9kOq/f0uaBb5Ek4Nf6sOlRETGSpKZ5NfAbScPSfb5T0n1pc9Rakg+70ZltV0fEhk77eyozvR9JDfTZTFlfQ1KTL9boiGiMiIkRcVUa1wQlTXEr07/5f3SKq7OnelhWuIi9dy9xZPcxGtiJ7f+Gr0unP0NSe/6tpEclnQMQEb8GvkHyLWuVpNmSRnQXV6e/Y+E8HEPyAZ6NZ8t0H49R95zgq9+FwOtJmkpGACek85X+7Dwc6FPAXRExMvPaLSLOK+JYT5E0WWwnTXQ/JPmwOBO4sYf9PEOS/Aqa0nkDMYKkBn9T2ib+QDp/haQ397ZxRGwErgX2J2nn3hn4MTAL2Cv9ELiDreUK25dt53lPkdR0R2fKekREHNa3X207V5N82zg4/Zv/W6e4eoqps9+ncZ7eyzGz+1hDUovu/Dd8GiAiVkbE30fEPiQ1+28VrglExFUR8QbgMGAC8OlejtvZapLmm3GZeftuE+jAj1E3nOCry46ShmVeO5C0ga4nuSg5iuQCY9Yqkh4iBbcDEySdKWnH9PVGSROLOP7twFhJn1ByUXR3Scdklt9A8nX+VOC7PexnLvBZSWMkjSbpsdHT+l3ZOVsWwAsk3wKOSF8np+u9Abi/t52l33jOJinLJ0hqqDuTJhRJ7wTe1pcAI+JZknb9KySNkDRE0oGSTuzLfrqwO0nT0ouSDgGK+XDuLsYAPgV8TslF5UKcUyTN7mabTSQf5pel58B+6T6+CyDpfUov9JM0LwWwKT3PjlFy4f0lYANJ01hf4t0EzCO50LtL+vt/qLC8FMeoJ07w1eUOkgRUeLUAXyO56LWG5MLlzzttcyXwXiU9bK6KiBdIEtX7SWrNK9l6obBH6bZvJek+txL4A5leGRFxD7AZeDAilvewqy8Ai4BW4BHgwXReX7zItmUxLa05royIlSSJGWBVRLzaw35+J+lFkkT0YeA9EfF8+rteQJLI2oG/I+ll01cfIvmweCzdz8303hzSm5lpPC8A/w3cNJCdRcTNJNdQziE5J1aR/D1+2sNm55Mk0CeAhcD3gevSZW8E7k/L9VaSZrM/k3zL+m+ScniSpHloVj9C/ieggeQcvJGkwvBKuqxUx6gLSi9cmBVFSdfG70fEtZWOxeqDpC8BYyPiw5WOpda4Bm9Fk/RGkm6XA6pRmvVEyX0czWmf96OBGcBPKh1XLfIda1YUSd8h6f/98bR5w6xcdidpltmHpPvkFfTcnGTdcBONmVlOuYnGzCynnODNzHKqqtrgR48eHePHj690GGZmNWPx4sVrImJMV8uqKsGPHz+eRYsWVToMM7OaIenJ7pa5icbMLKec4M3McsoJ3swsp5zgzcxyygnezCynnODNzDJaV7Zy+udOp3G/RoYMHcJBhxzE3LlzKx1Wv1RVN0kzs0pqXdnKP375H2md28ppnz6NEQeP4PGHHufT/5I8U+SMM86ocIR94wRvZpaat2wej817jHd/5t3sf+T+AAx941D23n1vLrvssppL8G6iMTNLtXW00f5UO01/1bRlXsOwBtQkli5dWsHI+scJ3sws1dTQROO+jbQ90rZlXseGDqItmDixmKdeVhcneDOz1PRDpnPo9EO55cu38MSDT7DmhTUsfWApd37tTi666KJKh9dnboM3M0s1j23mm5/5Jpfsegm3XXkbHSs6OODgA7j8i5fXXPs7VNkDPyZPnhwebMzMrHiSFkfE5K6WuYnGzCynnODNzHLKCd7MLKec4M3McsoJ3swsp5zgzcxyygnezCynnODNzHLKCd7MLKec4M3McsoJ3swsp5zgzcxyygnezCynnODNzHLKCd7MLKec4M3McsoJ3swsp5zgzcxyquwJXtJQSQ9Jur3cxzIzs60Gowb/cWDpIBzHzMwyyprgJY0D3gVcW87jmJnZ9spdg/8a8Blgc5mPY2ZmnZQtwUv6a+C5iFjcy3rnSlokadHq1avLFY6ZWd0pZw3+eOBUScuBHwBvkfTdzitFxOyImBwRk8eMGVPGcMzM6kvZEnxE/GtEjIuI8cD7gV9HxAfLdTwzM9uW+8GbmeXUDoNxkIhYACwYjGOZmVnCNXgzs5xygjczyykneDOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznHKCNzPLKSd4M7OccoI3M8spJ3gzs5xygjczyykneDOznHKCNzPLqbIleEnDJP1W0u8kPSrpknIdy8zMtrdDbytIGgIcDuwDrAcejYhVRez7FeAtEfGipB2BhZJ+FhH3DShiMzMrSrcJXtKBwD8DJwF/AFYDw4AJkl4GrgG+ExGbu9o+IgJ4MX27Y/qK0oVuZmY96akG/wXgauCjabLeQtKewN8BZwLf6W4HkoYCi4GDgG9GxP0DjtjMzIrSbYKPiDN6WPYc8LXedh4Rm4AjJI0EfiJpUkQsya4j6VzgXICmpqbiojYzs171epFV0vsk7Z5Of1bSPElH9eUgEbEWWAC8o4tlsyNickRMHjNmTF92a2ZmPSimF83nIuIFSVOAt5M0yVzd20aSxqQ1dyQNJ2nLXzaAWM3MrA+KSfCb0p/vAq6OiJ8COxWx3d7AfEmtwAPAryLi9v6FaWZmfdVrN0ngaUnXkNTAvyRpZ4r4YIiIVuDIAcZnZmb9VEwN/m+AXwDvSNvSRwGfLmdQZmY2cL3W4CPiZUnzgX0zF1fXlDcsMzMbqGLuZL0UOAv4E1tvVArgLeULy8zMBqqYNvi/AQ6MiFfLHYyZmZVOMW3wS4CRZY7DzMxKrJga/H8CD0laQjKAGAARcWrZojIzswErJsF/B/gS8AjQ5cBiZmZWfYpJ8Gsi4qqyR2JmZiVVTIJfLOk/gVvZtonmwbJFZWZmA1ZMgi/cjXpsZp67SZqZVblibnSaNhiBmJlZaXXbTVLSB9PH9XW3/MB0hEkzM6tCPdXg9yDpHrmY5KlMhUf2HQScSDJcwb+UPUIzM+uXnp7odKWkb5C0tR8PNJM8dHspcGZEtA1OiGZm1h89tsGnj9z7VfoyM7MaUsxQBWZmVoOc4M3McqqY4YKHpk011oXWVpg3D9raoKkJpk+H5uZKR2VmVlwN/o+SLpd0aNmjqTGtrTBrFrS3w7hxyc9Zs5L5ZmaVVkyCbwYeB66VdJ+kcyWNKHNcNWHePGhsTF5Dhmydnjev0pGZmYEiove1CitLJwBzScaHvxm4NCL+WKpgJk+eHIsWLSrV7vpMUkn315eyNTPrD0mLI2JyV8uKaoMH3gWcDYwHrgC+B7wZuAOYULJIK6yYhCxpy3otLUmzTGPj1uWF9y0t5YnRzKxYxTTR/AE4Dbg8Io6MiK9ExKqIuBn4eXnDq27TpycJvb0dNm/eOj19eqUjMzMrbjTJ5oh4sasFEXFBieOpKc3NMHPmtr1oZsxwLxozqw7FJPhvSvp4RKwFkNQIXBER55Q1shrR3OyEbmbVqaheNIXkDhAR7WwdI97MzKpUMQl+SFprB0DSKIqr+ZuZWQUVk6ivAP5P0s3p+/cBl5UvJDMzK4Vinuh0Qzom/DRAwPSIeKzskZmZ2YAU29SyDGgvrC+pyePBm5lVt2JudDofuBhYBWwiqcUHyRAGZmZWpYqpwX8ceH1E/KXcwZiZWekU04vmKaCj3IGYmVlpFVODfwJYIOl/gFcKMyPiK2WLyszMBqyYBN+WvnZKX0WRtC9wAzAW2AzMjogr+xOkmZn1XTHdJC8BkLRrRLzUh32/BlwYEQ9K2h1YLOlX7mJpZjY4em2Dl3ScpMeApen7wyV9q7ftIuLZiHgwnX4h3f51A4zXzMyKVMxF1q8Bbwf+AhARvwNO6MtBJI0nGb/m/r6FZ2Zm/VVMgicinuo0q+iHcEvaDfgx8ImIWNfF8nMlLZK0aPXq1cXu1szMelHMRdanJL0JCEk7AReQNtf0RtKOJMn9exHR5ZNKI2I2MBuSR/YVFbWZ2SBqXdnKvGXzaOtoo6mhiemHTKd5bPXf61lMDf4fgH8kaT9fARwBfKy3jZQ84HQOsNRdKs2sVrWubGXWvbNoX9/OuBHjaF/fzqx7Z9G6srXSofWqmAT/+oj4QETsFRF7RsQHgYlFbHc8cCbwFkkPp6+TBxStmdkgm7dsHo3DGmkc3sgQDaFxeCONwxqZt6zLRomqUkwTzdeBo4qYt42IWEgybo2ZWVVJGhgG7hIuASCiOluXu03wko4D3gSMkfSpzKIRwNByB2ZmVi7FJmRJXDz/YtrXt9M4fMtzj7a8b5naUqYIS6OnJpqdgN1IPgR2z7zWAe8tf2hmZpU3/ZDptG9op319O5tjM+3r22nf0M70Q6ZXOrReqbdPMkn7RcSTgxHM5MmTY9GiRYNxqH6TVLVfx8ystAr/79Xci0bS4oiY3NWyYtrgX5Z0OXAYMKwwMyLeUqL4zMyqWvPY5qpJ6H1RTC+a75E80Wl/4BJgOfBAGWMyM+u3UaNGIakkL6Ak+xk1alRFyqKYGvweETFH0scj4i7gLkl3lTswM7P+aG9vr7pm1FL12umrYhL8xvTns5LeBTwDjCtfSGZmlVfN7e7FKqaJ5guSGoALgZnAtcAnyhmUmVkl1fLdq1m9JviIuD0iOiJiSURMi4g3AAcOQmxmZhVRy3evZhU1mmQXPtX7KmZmtamto42GYQ3bzGsY1kBbR1uFIuqf/iZ4D0FgZrnV1NBEx4aObeZ1bOigqaGpQhH1T38TfHVdojYzK6Favns1q9sEL+kFSeu6eL0A7DOIMZqZDarmsc3MPG4mjcMbWbFuBY3DG5l53Mya60XTbTfJiNh9MAMxM6smtXr3alZ/m2jMzKzKOcGbmeWUE7yZWU45wZuZ5ZQTvJlZTjnBm5nllBO8mVlOOcGb1ajWVmhpgXPOgdNPn8tBB01i6NChTJo0iblz51Y6PKsCxYwHb2ZVprUVZs2CxkZYu3Yud955EcccM4ebbprCunULmTFjBgBnnHFGhSO1SnINvh+yNaeWluS92WCaNy9J7o2NsHDhZbz73XOYOHEat922I9OmTWPOnDlcdtlllQ7TKswJvh9mzYL2dhg3Lvk5a5aTfHdaV7bSsqCFc356Di0LWrY8MGHu3LlMmuQmhf5qa4OGdDTbNWuW0tQ0hYaGZD7AlClTWLp0aeUCtKrgJpp+KNScCtOQ1Kiaa3vYipIrPBWncVjjNk/FOWzlYVxz+TXMmTOHKVOmsHChmxT6qqkpqVw0NsLo0RNpa1vIyJHTaEpHs124cCETJ06sbJBWca7B90NDw/bv22rrOQCDorun4lzx5SuYM2cO06ZNY8cd3aTQH9OnJwm+vR2mTLmIW26ZwdKl8znllI3Mnz+fGTNmcNFFF1U6TKsw1+D7oaNja8298L6ptp4DMCjaOtoYN2Lb57M3DGtgTdsapkyZss18Nyn0TXMzzJyZfHN84YUzOOkk+N3vzufoo5cyceJELrvsMn8bMif4/mhvT342NCTJvb0d0hYGy2hqaKJ9fTuNw7d+GnZs6GB002gWLlzItGnTtsx3k0LfNTdnmwXPSF9mW7mJph9mzkxq8CtWJD9nznT7e1e6eyrOhZ+5kBkzZjB//nw2bnSTglm5uAZfhNbW5Ktwtp29paVi4dSMwlNx5i2bR1tHG00NTcw4cgbNY5tpGtnE+eefz9KlblIwKxdFVM/jVSdPnhyLFi2qdBhbtLbC1VfDr34Fe+wBRxwBs2eLM88M19rNqpQkqimvQXljkrQ4IiZ3tcxNNN0o3Cn44IMwalQy7777kp+NjUmN3sysmpWtBi/pOuCvgeciYlIx25SzBj9q1CjaC1dHq0RjYyPPP/98pcOwnMk2KTY1JV0q6+nbpmvwW5WzDf564BvADWU8RtHa29v7VMDnnJPcqXr33bB+PQwfDhGwbh2ccEJSix9oO7ykge3ArJPsGDXZO63dpFifytZEExF3AzVbPW1qSrpATpwIGzYkSX79ethpp+SfZvr0Skdotr3sGDVDhmyddpNifXIvmm5Mn761JnTssfDww/D883DSSfCxj7k2ZIOrFN/2Lrlk63S1NWFYeVQ8wUs6FzgXoKmKbgfN3inY1gbvfGf9tWVa9Sg2IUviggtimzutC2PWuGtv/al4go+I2cBsSC6yVjicbWx7p6BZbfCd1lbgbpJmOeM7ra2gbDV4SXOBqcBoSSuAiyNiTrmOZ2YJf/O0grIl+IjwfedmVlNaV7ZuM7TG9EOm0zy2dj8t3URjZsbWB9S0r2/f5gE1haeQ1SIneDMzun9AzbxltXsTQd0MNlZvty9bffI5BbQ09L5OJbR0lGW3lRqqwMxs0OmSdf36kGtZ0LLdA2oK71umtgwsJokY2C76xU00ZmZsfUDN42seZ/7y+fzw0R+yYPkCJo0paqzEquQavFkdqPcRJovRPLaZUyecyqW/uZSNmzYyZpcxjBsxjlsfv5UJe0yoyd40TvBmOecRJou3ZPUSpu43dbtmmnnL5tVkgncTjVnOeYTJ4rV1tNEwbNuLtA3DGmjraOtmi+rmGrxZhZX6YTR9GXkyO8JkZ/X4QJqmhqbtLrR2bOigqaF6BkLsCyd4swrr68No+qqlZeuIkluP2fsIk/X4QJrph0xn1r2zgKTm3rGhg/YN7cw4sjZHa3MTjVnOTZ+eJPT2dti8eeu0H1qzveaxzcw8biaNwxtZsW4FjcMbmXnczJpsfwff6FRR1RiTDb7BOA/604umVs/Paow7j89kzRV3MyuNvA3mVCl9PR89wmR9chNNEQrdzNrbt+1m1lq7YxBVRB4Hc6oEn4+9k1RVr8bsBZBB5ARfBHczK408DuZUCT4fexYRJXuVan+V6o1UN000cfGIfg9C1ALQ+QO4Mbuwf+LiEf3fuAb0qbsePfTXS1Vbu2qltLUlNfeshoZkvllW3ST4/g5ABP3vZtZrTBUagGiwdC7vrgZzumTaJVw8/+IBD+ZUT5qatj8fOzqS+QW+ZmTgJpqiuJtZaRQGc2pf387m2Ez7+vYt8614vZ2PbqMvjcK1oXN+eg4tC1pq8lqRE3wRmpv9IONS6KqPcWG+Fa+389Ft9ANX6BAA1HSHAPeDr6BqjGmwuQyo3gdUQNkeUlHtCs2JV518FRfPvxgo3djwpeZ+8GZVbCDXh7qSvWa0ahUsXQrPPQd77glXXVX8N8+8XyPqSVtHG+NGbHsluxYHHXMTjVnOFNro//AHuOceWLsWdtgB9tnHbfHFampoomPDtt9eanHQMSf4IrS2JrWic85JfvofxKpZoY3+6afhpZeSBP/KK/DMM7Bpk9vii1HoEABs6RDQvqG95joEOMH3wj0SrBY1N8PIkbDrrrDXXjB2LKxfD488Ag8/XOnoql+hQwBQ04OOuQ2+F9keCbD157x59duLppLjl/ekHscv78natUkvmuHDk/fDhyc1+bVrKxlV7Sgk8+tOu67CkfSfE3wvfNfg9so9fnl/1fL45YMdezGHq9T4KVY6TvC9KOauQbOBKOWHZbbbaUsLPP540vbe0ZFUTPbZByZMGNgd2FY73Abfi+7uGpw0yRderbpNn570njn8cDjllOTnDjv4Dux64hp8Lwo9ErLjerz5zXDrrX5KvVWH7LgzADffDEuWJO933TVpd1+xIjl3Z8zwOdqVvD6nwHey9kNJbySpovIvVl/jHqx/nlotz4Eo9PJqbEyaYC69VOy/f3DssXDQQUnTTHu7Kx89KQxL0DiscZvnsM48biaH73141Z9TvpM1Ve4LWUuXJl+Di1UPF7Gy/zzZMT0KXdDyWGsaTJ17eQGMGJG0u0+Y4F5fxcg+pwDY8jMPzymomxr8QGRrSVddJY4+Oli3Lvn6O3Ro0v1s/frk5+GHD3wY4WrXl5pyV0MEt69v59VNr/LSxpe6rDX1N8nnuQZf6spJXsupWHkqT9fgB6hzLenII+Guu+DJJ2HixCS5b9iQzK/3LpSddTemx22/v40T9juhy1qTa/Hb6y6BdH5WwYIFST/3kSNh6tRkXimeXZA32fLsrhJSjQOL9ZVr8Bl5+lQvK49+WDU6t8H/6U9w7724Db4PemqDr4XKRk81+LImeEnvAK4EhgLXRsQXe1q/0gm+O9090enVV5OxPgr/XPXyz9SXppDCP8+mzZt4ZNUjDNEQXn7tZVa9uIohGsIBjQeww5Ad2PDaBibtOYmD9zi437WmPDfR9KTz05smTdrai8ZPcypOLfeiqUiClzQUeBx4K7ACeAA4IyIe626bak3wnWtJ2UQO9fdotP70orng5xfw3EvPscuOu/D0uqdZv3E9r256ld122o2JYybSsaGDUHDju290G7xZH1SqDf5o4I8R8UQaxA+A04BuE3y16qovfLY/cd4TeldK1Zy1nvWsZvWW94ef14duSJ3UQ68ks74oZ4J/HfBU5v0K4JgyHq+smpvrM5F3pT+15JYFLfzsDz8DYNVLq9i4aSMAz8x8hhO/fSKH73V4Li5qmVWTcib4rqp422UGSecC5wI0eYCXXCm2ln/X2XdxF3cBcAmXdLuem1/M+qacCX4FsG/m/Tjgmc4rRcRsYDYkbfBljMcGWeeEfN7t5/Hgsw/y6qZXGTpkKAAvb3yZPXfdk6vecVXNXNQyqxXlTPAPAAdL2h94Gng/8HdlPJ5VufMmn1fT3dHMak3ZRpOMiNeAfwJ+ASwFfhgRj5breFb9Ck/JaRzeWNNPyTGrFWW9kzUi7gDuKOcxrLY0j212QjcbJB4P3swsp5zgzcxyygnezCynnODNzHLKCd7MLKeqarhgSauBJysdRy9GA2sqHUSOuDxLy+VZWrVQnvtFxJiuFlRVgq8FkhZ1N3Kb9Z3Ls7RcnqVV6+XpJhozs5xygjczyykn+L6bXekAcsblWVouz9Kq6fJ0G7yZWU65Bm9mllNO8BmSXuxi3qckPSapVdKdkvbrtPznkl4naY6k36Xr3Sxpt8GLvPq4LMuvqzLOLHuvpJBUsz1ABlsey9MJvncPAZMjohm4GfhyYYGk4cCoiHga+GREHJ6u10YyVLJty2U5CCTtDlwA3F/pWPKglsvTCb4XETE/Il5O395H8mSqgqnAgnS9dQBKnlM3nC4eT1jvXJaD5lKSD88NlQ4kJ2q2PJ3g+2YG8LPM+3cCPy+8kfRtYCVwCPD1wQ2t5rgsy0DSkcC+EXF7pWPJg1ovTyf4Ikn6IDAZuDwz+3hgYeFNRJwN7EPyBKu/HdQAa4jLsjwkDQG+ClxY6VjyIA/l6QRfBEknARcBp0bEK+m8A4CnIuLV7LoRsQm4CTh90AOtAS7LstodmAQskLQcOBa4tdYuDFaRmi9PJ/hepF/RriFJSM9lFm1pUlDioMI0cAqwbLBjrXYuy/KKiI6IGB0R4yNiPMl1jlMjYlGFQ6tJeSjPsj6TtQbtImlF5v1XgJOB3YAfJfmGtog4FXgHcH66noDvSBqRTv8OOG/Qoq5OLsvy266MI+IrFYum9uWuPH0naz9I2hm4p5ZHmasWLkuz8nGCNzPLKbfBm5nlVN0l+J5uR+5i3amS3lTCY4+XtKRU+6sW6S3cN2be7yBptaQ+9R2WtKDQQ0HSHZJG9rL+ckmj+xV0FfK5WXr1fm7WXYLvo6lAn/6JJNXjheuXgEnpcAMAbwWeHsgOI+LkiFg70MBybCo+N4tR1+emEzwg6RRJ90t6SNL/StpL0njgH4BPSnpY0psljZH0Y0kPpK/j0+1bJM2W9EvghrQ29BtJD6avktW0qtjPgHel02cAcwsLJO0q6bq0zB6SdFo6f7ikH6SDit1EMixBYZstNSBJt0haLOlRSecO3q9UeT43S6J+z82IqKsX8GIX8xrZesH5I8AV6XQLMDOz3veBKel0E7A0s95iYHj6fhdgWDp9MLAonR4PLKl0GZSjTIHCAGLDgIdJapi3p8v/A/hgOj0SeBzYFfgUcF06vxl4jWQwMoDlwOh0elT6cziwBNij8zp5ePnc9LlZ6lc9fmXryjjgJkl7AzsBf+5mvZOAQ9M+3AAjlIw0B3BrRKxPp3cEviHpCGATMKEsUVeRiGhNa5ZnAHd0Wvw24FRJM9P3w0iS0AnAVZntW7vZ/QWS3pNO70uSmP5SwvCrmc/NAarnc9MJPvF1kpsabpU0laTW05UhwHGZfxYA0n+qlzKzPgmsAg5Pt6m5Uej66VZgFkkNaY/MfAGnR8Tvsyun5dZjP93073ESSbm/LGkByT9hvfC5WRp1eW66DT7RwNYLLx/OzH+BZDyKgl+SGZs8rQV1t79nI2IzcCYwtGSRVrfrgM9HxCOd5v8COF/pf42SIQsA7gY+kM6bRPJVuLMGoD39BzqEZDyQeuJzszTq8tysxwS/i6QVmdenSGpFP5L0G2BNZt3bgPcULmSRDPo/Ob3w8hjJha6ufAv4sKT7SL4Cv9TNerkSESsi4souFl1K0jTQqqQr3qXp/KuB3dKvv58BftvFtj8HdkjXuZRkPJC88rlZJvV6bvpOVjOznKrHGryZWV1wgjczyykneDOznKr7BJ/esXaXpH73JpDUKOkn6QWu36ZX3ZG0k6S7VUe3iJeoPD+dXjx8WNISSZskjaq38vS5WVr1eG7WfYIHzgHmRfJ4uP76N+DhiGgGPgRcCRDJI+jupL6eKTrg8oyIyyPiiIg4AvhX4K6IeL4Oy9PnZmnV3bnpBJ/0df2pktH5Fki6WdIySd8r9I0twqEkf1wiYhkwXtJe6bJb0mPUi1KUZ9Y2Y4dQX+Xpc7O06u7crKqvE4NN0k7AARGxXMmtzEcChwHPAPcAxwMLJX0VmNbFLn4QEV8keazc9HTdo4H9SG4xX0UyPsUby/27VIMSlmdhf7uQPM7vnzLr1EV5+twsrXo9N+s6wQOjgbWZ97+NiBUAkh4mGYBpYUR8spf9fBG4Mt3mEeAhksGJiIhNkl6VtHtEvFDS6KtPqcqz4BSSx/k9X5hRR+Xpc7O06vLcrPcEv55tx454JTO9ibR8evtUj4h1wNnpuiIZECo7KNTO1MeYHyUpz8z797PtV+CCeihPn5ulVZfnZl0n+IholzRUUo8DBPX2qa7k6S4vpxdaPgLcnf5jIWkPYHVEbCxR2FWrVOUJIKkBOBH4YKf5dVGePjdLq17PzbpO8KlfAlNIv7b200SShylsAh4DZmSWTWP7IUrzrBTlCfAe4JcR0XmslHoqT5+bpVV352bdj0WjZPS4T0XEmWXa/zzgXzsPR5pXLs/ScVmWVj2WZ913k4yIh4D5A7n5oTvplftbqukPXm4uz9JxWZZWPZZn3dfgzczyqu5r8GZmeeUEb2aWU07wZmY55QRvFSPpIkmPKhnp8GFJx6TzP5HeCl5Y7460P3cpjvliL8tHSvpY5v0+km4u0bEXSJpcin2ZFcMJ3ipC0nHAXwNHpSMdngQ8lS7+BLAlwUfEyRGxdpBCGwlsSfAR8UxEvHeQjl0S5eglYrXJCd4qZW9gTUS8AhARayLiGUkXAPuQdGebDyBpuaTRksano/9dq2Qs7u9JOknSPZL+kA6mhaQWSTMLB0rXHZ89uKTdJN0p6UFJj0g6LV30ReDA9BvF5ekxl6TbDJP07XT9hyRNS+efJWmepJ+ncXy52EJI9/+bNI4HJb0pnX9jJibS3/XU9G7MyyU9kH7z+Wi6fKqk+ZK+TzLmjBlEhF9+DfoL2A14GHgc+BZwYmbZcmB05/ckA0K9BvwVSeVkMXAdIOA0kn7IAC3AzMz2S4Dx6fSL6c8dgBHp9Gjgj+l+xgNLMttueQ9cCHw7nT4EaCMZ3+Qs4AmgIX3/JLBvF7/zAmByp3m7AMPS6YOBRen0iZnfp4Fk/JgdgHOBz6bzdwYWAfsDU4GXgP0r/bf1q3persFbRUTEi8AbSBLWauAmSWcVsemfI+KRiNgMPArcGRFBUmsd34cQBPyHpFbgf4HXAXv1vAlTgBvT+JeRJPIJ6bI7I6IjIjaQDAmwX5Fx7Aj8t6RHgB+RjN9ORNwFHCRpT5Jxx38cEa8BbwM+lI6AeD+wB8kHAyQjJP4Zs5THorGKieTJOguABWmC+zBwfS+bZUcB3Jx5v5mt5/NrbNv82NUAUx8AxgBviIiNkpZ3s15WTw+F6HJ0wiJ8kmRs9sNJYs6ORHhjGuf7SZ5GVIjh/Ij4xTaBSVNJavBmW7gGbxUh6fWSDs7MOoKkRgzwArD7AHa/HDgqPc5RJE0YnTUAz6XJfRpba9w9Hftu0if2SJoANAEDvTW9AXg2/UZyJpC9QHo9yQVnIuLRdN4vgPMk7ViIQ9KuA4zBcso1eKuU3YCvp90fXyNpAz83XTYb+JmkZyOiq7G5e/NjtjZjPEDSzt/Z94DbJC0iuRawDCAi/pJetF0C/Az4ZmabbwH/lX7beA04KyJeUd+e9vY/kgrDyd5L8szUH0t6HzCfTC08IlZJWkryKLiCa0maoh5UcuDVwLv7EoDVD49FY1al0nsBHiHpStpR6Xis9riJxqwKSTqJ5FvF153crb9cgzczyynX4M3McsoJ3swsp5zgzcxyygnezCynnODNzHLKCd7MLKf+PwgOWKWXQQzBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# latency plots\n",
    "plt.clf()\n",
    "i = 1\n",
    "xtick_labels = []\n",
    "df_datasets = []\n",
    "for stim_loc in ['L2/3', 'L4']:\n",
    "    df2 = df[df['Stim Location'].str.contains(stim_loc)]\n",
    "    df_medial3 = df2[df2[\"Most Medial Barrel\"] == 3]\n",
    "    df_medial1 = df2[df2[\"Most Medial Barrel\"] == 1]\n",
    "    medial_lat = []\n",
    "    lateral_lat = []\n",
    "    medial_lat.append(df_medial1[(df_medial1['abs(delta 2)'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 3)]['abs(delta 2)'])\n",
    "    lateral_lat.append(df_medial3[(df_medial3['abs(delta 2)'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 3)]['abs(delta 2)'])\n",
    "    medial_lat.append(df_medial3[(df_medial3['abs(delta 2)'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)]['abs(delta 2)'])\n",
    "    lateral_lat.append(df_medial3[(df_medial3['abs(delta 1)'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 2)]['abs(delta 1)'])\n",
    "    medial_lat.append(df_medial3[(df_medial3['abs(delta 1)'] > 0)\n",
    "                                   & (df_medial3['Stim Barrel'] == 1)]['abs(delta 1)'])\n",
    "    lateral_lat.append(df_medial1[(df_medial1['abs(delta 1)'] > 0)\n",
    "                                   & (df_medial1['Stim Barrel'] == 1)]['abs(delta 1)'])\n",
    "    medial_lat.append(df_medial1[(df_medial1['abs(delta 1)'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)]['abs(delta 1)'])\n",
    "    lateral_lat.append(df_medial1[(df_medial1['abs(delta 2)'] > 0) \n",
    "                                   & (df_medial1['Stim Barrel'] == 2)]['abs(delta 2)'])\n",
    "    \n",
    "    medial_lat = pd.concat(medial_lat)\n",
    "    lateral_lat = pd.concat(lateral_lat)\n",
    "    \n",
    "    #dfd1_ant = df2_ant['abs(delta 1)'].dropna()\n",
    "    #dfd2_ant = df2_ant['abs(delta 2)'].dropna()\n",
    "    \n",
    "    #df2_post = df2[df2['Slice Location'] == \"Posterior\"]\n",
    "    #dfd1_post = df2_post['abs(delta 1)'].dropna()\n",
    "    #dfd2_post = df2_post['abs(delta 2)'].dropna()\n",
    "    \n",
    "    #dfds_ant = np.concatenate([dfd1_ant, dfd2_ant])\n",
    "    #dfds_post = np.concatenate([dfd1_post, dfd2_post])\n",
    "    \n",
    "    jittered_xs_1 = np.random.normal(i, 0.04, lateral_lat.shape[0])\n",
    "    jittered_xs_2 = np.random.normal(i+1, 0.04, medial_lat.shape[0])\n",
    "    plt.scatter(jittered_xs_1, lateral_lat, alpha=0.5, color='blue')\n",
    "    plt.scatter(jittered_xs_2, medial_lat, alpha=0.5, color='green')\n",
    "\n",
    "    xtick_labels.append(stim_loc + \"\\nLateral\\n(n=\" + str(len(lateral_lat)) +\")\")\n",
    "    df_datasets.append(lateral_lat)\n",
    "    xtick_labels.append(stim_loc + \"\\nMedial\\n(n=\" + str(len(medial_lat)) +\")\")\n",
    "    df_datasets.append(medial_lat)\n",
    "    i += 2\n",
    "\n",
    "\"\"\"# calculate significance\n",
    "for sl1 in ['ant', 'post']:\n",
    "    for sl2 in ['ant', 'post']:\n",
    "        for la1 in ['L2/3', 'L4', 'L5']:\n",
    "            for la2 in ['L2/3', 'L4', 'L5']:\n",
    "                if sl1 == sl2 and la1 == la2:\n",
    "                    continue\n",
    "                elif len(dists[sl1][la1]) < 3 or len(dists[sl2][la2]) < 3:\n",
    "                    continue\n",
    "                else:\n",
    "                    p_val = kstest(dists[sl1][la1], dists[sl2][la2]).pvalue\n",
    "                    if p_val < 0.05:\n",
    "                        print('KS test for ' + sl1 + \" layer \" + la1 + \" and\", \n",
    "                              sl2 + \" layer \" + la2,\n",
    "                             \"is significant with p-value \", p_val)\n",
    "\"\"\"\n",
    "# significance: L2/3 Rostral and L4 Rostral, and \n",
    "#               L2/3 Rostral and L4 Caudal\n",
    "#significance_bar(0.75, 2.25, 3.35, \"P<0.05\", fontsize=10)\n",
    "\n",
    "plt.boxplot(df_datasets)\n",
    "plt.xticks(range(1, len(xtick_labels)+1), xtick_labels)\n",
    "plt.title(\"Latency of L4 Barrel Pair Crossings\")\n",
    "plt.xlabel(\"Stimulation Layer\")\n",
    "plt.ylabel(\"Latency (ms)\")\n",
    "plt.savefig(data_dir + cell_type_prefix + \"_results_latency_medial-bias_\" + str(datetime.date.today()) +\".png\", bbox_inches=\"tight\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "0ff4e0ce",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KS test p-value  0.3620106185909979\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEiCAYAAADziMk3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAr+klEQVR4nO3dfZyUdb3/8dcbBBcFlhFQVNrIjMA4qymZlqewYzfHyk6b3XA8nkxOdne8KalfZSZm9uic6M7SytTMbqiTbmpmx+5EozQFqgUF9US6oqKgwyK3cfP5/XFdA7Pr7uy1szszy+z7+XjMg2uu28/MDp/5zvf63igiMDOz+jOs1gGYmVllOMGbmdUpJ3gzszrlBG9mVqec4M3M6pQTvJlZnXKCt6qT9AFJT0raKGl8reOpFElnSFpUgfNulHTYQJ+3kiR9U9KFtY5jqHGCHyQkPSzppDKOWyjpPyoRUyVIGgF8CXhdRIyOiKe7bJ8iKSTt08t5ftvbfun2TWlCXCdpgaRxA/JC+qnodW5MHw9L+niWY9P3bVUfr/d6SXdKelbSWkl3SDqlvOj7LiLeHxGXVOt6lnCCt2o7CGgA7iv3BJJOA0p+ARQ5MiJGA4cBOWBemdd8zvV6+xLKaFwa32zg05Le0J+TSRrezbpTgZ8A1wGTSf4Gnwbe3MM5BuJ12SDgBD/IScpJuiUtdeXT5cnptkuBfwS+npYCv56unybpV5KekfSApHcUne9aSZdL+nlamvujpBcWbX9J0bFPSvqkpEmSNhdXp0g6Jo1pRDcx7yvpK5IeTx9fSddNBR5Id1sv6bdlvB+NwEXAx/pyXERsAG4Gjig613skrUjfh1WS3le0bZak1ZL+n6Q1wHckzZN0vaTvS9oAnCGpUdLVkp6Q9Jikz3aXZDPEdxfJl94MScdKukvS+vS8X5c0sii2kHR4unytpG9IulXSJuDELu+XSH4xXRIRV0VER0Tsiog7IuK96T5nSPq9pC9LegaYl76u69K/8SOSPiVpWLr/4ekvgI70l9GPC9dKz/FUuq1N0oyiOD/b5b09P933CUnvKYp5vKSfSdog6d70PV3U2zWsGxHhxyB4AA8DJ3WzfjzwNmA/YAxJSezGou0Lgf8oer4/8CjwHpJS7tHAOuAl6fZrgWeAY9PtPwB+lG4bAzwBnE9Syh4DvDzddivwgaLrfBn4Wg+v5TPA3cCBwETgDyQJBmAKEMA+PRzb2/bLgQ/3tl+6bwCHp8s54JfAZ4q2vxF4ISDg1cBm4Oh02yxgB/BfwL7AKJLS/3bgX0gKR6OAG4Fvpe/7gcA9wPvSc5wBLOrtdabXf2V6/X8CjgGOS7dNAVYA5/Xwuq4FOtLjhwENXa4zLd3/BSXepzPS13p2es1RJKX9m9LPwBTgQWBOuv8C4ILC9YAT0vWvB5YA49LXNB04uCjOz3Z5bz8DjABOTl97Lt3+o/SxH8kX8qOF97HUNfzo5m9b6wD8SP8QPST4bvY7CsgXPV9I5wT/TuB3XY75FnBRunwtcFXRtpOBlenybOBPPVz3ncDv0+XhwBrg2B72/StwctHz1wMPp8u7E1sPx/a4HZgJ/Jk9iS9Lgt8ArAd2AiuBQ0vsfyNwbro8C/g7RQmTJMHfWfT8IGAbMKpo3Wzg9nT5DHpP8OuBPEkSP6eHfc8DftrldRUn+OtKvKZXpvs3lNjnDKC96Pnw9HUdUbTufcDCdPk64EpgcpfzvIbki+A4YFiXbdfSOcFvKf7bAU+lxw0n+RJ9cdG2z7Inwfd4DT+e+3AVzSAnaT9J30p/Jm8A7gTGlagGeD7w8vTn/XpJ64HTgElF+6wpWt4MjE6Xn0eSnLtzE3CEktYbrwU6IuKeHvY9BHik6Pkj6bqypdUDV5Ak4B19OPToiBhHUtL8BvA7SQ3pOf9Z0t1pddR6ki+7CUXHro2IrV3O92jR8vNJSqBPFL3X3yIpyWc1ISJyETE9Ii5L45qqpCpuTfo3/1yXuLp6tMS2wk3sg3uJo/gcE4CRPPdveGi6/DGS0vM9ku6TdCZARPwW+DrJr6wnJV0paWxPcXX5OxY+hxNJvsCL49m93MdrDHlO8IPf+cCLSapKxgKvStcr/bfrcKCPAndExLiix+iI+ECGaz1KUmXxHGmi+x+SL4vTge+VOM/jJMmvoCld1x9jSUrwP07rxO9N16+W9I+9HRwR24GrgBeQ1HPvC9wAzAcOSr8EbmXP+wrPfW+7rnuUpKQ7oei9HhsRL+nbS3uOb5D82nhR+jf/ZJe4SsXU1QNpnG/r5ZrF51hHUoru+jd8DCAi1kTEeyPiEJKS/RWFewIRcVlEHAO8BJgKfLSX63a1lqT6ZnLRuud1CrT/1xgynOAHlxGSGooe+5DUgW4huSl5AMkNxmJPkrQQKbgFmCrpdEkj0sfLJE3PcP1bgEmSzlNyU3SMpJcXbb+O5Of8KcD3S5xnAfApSRMlTSBpsVFq/+7sW/xeAM+S/Ao4Kn2cnO53DPDH3k6W/uJ5D8l7uYqkhLovaUKR9M/A6/oSYEQ8QVKv/0VJYyUNk/RCSa/uy3m6MYakammjpGlAli/nnmIM4CPAhUpuKhfiPEHSlT0cs5Pky/zS9DPw/PQc3weQ9HalN/pJqpcC2Jl+zl6u5Mb7JmArSdVYX+LdCbSS3OjdL339/17YPhDXGEqc4AeXW0kSUOExD/gKyU2vdSQ3Lv+3yzFfBU5V0sLmsoh4liRRvYuk1LyGPTcKS0qPfS1J87k1wEMUtcqIiN8Du4ClEfFwiVN9FlgMtAHLgKXpur7YSOf34sS05LgmItaQJGaAJyPi7yXO8xdJG0kS0buBt0bEM+lrPYckkeWBfyVpZdNX/07yZXF/ep7r6b06pDdz03ieBb4N/Lg/J4uI60nuoZxJ8pl4kuTvcVOJw84mSaCrgEXAD4Fr0m0vA/6Yvq83k1Sb/Y3kV9a3Sd6HR0iqh+aXEfJ/Ao0kn8HvkRQYtqXbBuoaQ4LSGxdmmShp2vjDiLiq1rHY0CDpv4BJEfHuWseyt3EJ3jKT9DKSZpf9KlGalaKkH0dz2ub9WGAO8NNax7U3co81y0TSd0naf5+bVm+YVcoYkmqZQ0iaT36R0tVJ1gNX0ZiZ1SlX0ZiZ1SkneDOzOjWo6uAnTJgQU6ZMqXUYZmZ7jSVLlqyLiIndbRtUCX7KlCksXry41mGYme01JD3S0zZX0ZiZ1SkneDOzOuUEb2ZWp5zgzczqlBO8mVmdcoI3s6pZsGABM2bMYPjw4Rz4ggOZdf4s5i2cR9uatlqHVpcGVTNJM6tfCxYs4IILLuCTX/gkC3ctZPNfN3P7ZbfTsE8Dq/KrmHv8XJonNdc6zLriEryZVcWll17K1Vdfzerxqxk/ejzNxzVzykdPYen/LCXXkKN1ZWutQ6w7TvBmVhUrVqzghBNOoL2jncaGRgCa/qGJde3raGxopL2jvcYR1h8neDOriunTp7No0SKaGpvo2NoBQPuydiY0TaBjawdNjU01jrD+OMGbWVVccMEFzJkzh8lPT+bpjU/TdncbN3/hZo5+x9Hkt+ZpmdZS6xDrjm+ymllVzJ49G4BLL7qUFStWML5pPEfOPpLj/vk4Wqa1+AZrBQyqCT9mzpwZHmzMzCw7SUsiYmZ321xFY2ZWp5zgzczqlBO8mVmdcoI3M6tTTvBmZnXKCd7MrE45wZuZ1SkneDOzOuUEb2ZWp5zgzczqlBO8mVmdcoI3M6tTvY4mKelA4JXAIcAWYDmwOCJ2VTg2MzPrhx4TvKQTgY8DBwB/Ap4CGoB/AV4o6XrgixGxoQpxmplZH5UqwZ8MvDcinjOPlqR9gDcBrwVuqFBsZmbWDz0m+Ij4aIltO4AbKxGQmZkNjF5vsko6V9JYJa6WtFTS66oRnJmZlS9LK5oz03r21wETgfcAn69oVGZm1m9ZErzSf08GvhMRfyla1/vB0nBJf5J0SzkBmplZebIk+CWSfkmS4G+TNAboSxPJc4EV5QRnZmbly5Lg55A0l3xZRGwGRpJU0/RK0mTgjcBVZUdoZmZl6bWjU0TskrQDeFXaPLKgLcP5vwJ8DBhTXnhmZlauLD1ZrwGagfvYUzUTQGsvx70JeCoilkiaVWK/s4CzAJqamjIFbWZmves1wQPHRcQRZZz7lcApkk4m6QE7VtL3I+LfineKiCuBKwFmzpwZZVzHzMy6kaUO/i5JfU7wEfGJiJgcEVOAdwG/7ZrczcyscrKU4L9LkuTXANtImkhGRDRXNDIzM+uXLAn+GuB0YBl9ax65W0QsBBaWc6yZmZUnS4Jvj4ibKx6JmZkNqCwJfqWkHwI/I6miASAiSraiMTOz2sqS4EeRJPbiAcZ6bSZpZma1laWjU6Zeq2ZmNrj02ExS0qckHVBi+2vSzkxmZjYIlSrBLwN+JmkrsBRYS9Jh6UXAUcCvgc9VOkAzMytPqRmdbgJukvQikl6pBwMbgO8DZ0XEluqEaGZm5chSB/8Q8FAVYjEzswGUZagCMzPbCznBm5nVKSd4M7M61WuClzRV0m8kLU+fN0v6VOVDMzOz/shSgv828AlgO0BEtJEM/2tmZoNYlgS/X0Tc02XdjkoEY2ZmAydLgl8n6YUk488g6VTgiYpGZWZm/ZZlsLEPkUypN03SY8DfAM/MZGY2yGXp6LQKOEnS/sCwiHi28mGZmVl/ZWlF8zlJ4yJiU0Q8Kykn6bPVCM7MzMqXpQ7+nyNifeFJROSBkysWkZmZDYgsCX64pH0LTySNAvYtsb+ZmQ0CWW6yfh/4jaTvkLSkORP4bkWjMjOzfstyk/W/JS0D/gkQcElE3FbxyMzMrF+ylOCJiF8Av6hwLGZmNoCytKJpkfSQpA5JGyQ9K2lDNYIzM7PyZSnB/zfw5ohYUelgzMxs4GRpRfOkk7uZ2d4nSwl+saQfAzcC2worI6K1UkGZmVn/ZUnwY4HNwOuK1gXgBG9mNohlaSb5nmoEYmZmA8szOpmZ1SnP6GRmVqc8o5OZWZ3yjE5mZnWq3BmdTqtoVGZm1m8lE7yk4cAHIqLPMzpJagDuJBlaeB/g+oi4qL8Bm5lZNiUTfETslHRMurypj+feBrwmIjZKGgEskvSLiLi7zFjNzKwPslTR/EnSzcBPgN1JvreerBERwMb06Yj0EWXGaWZmfZQlwR8APA28pmhdpp6saRXPEuBw4PKI+GM5QZqZWd9VtCdrROwEjpI0DvippBkRsbx4H0lnAWcBNDU1lXspMzProio9WdNJuxcCb+hm25URMTMiZk6cOLEvpzUzsxIq1pNV0sS05F6YqPskYGXZkZqZWZ9kqYPfLyLukVS8LktP1oOB76b18MOA/4mIW8qI0czMypAlwZfVkzUt6b+0f+GZmVm53JPVzKxO9ZjgJZ0bEV8FDi6nJ6uZmdVWqZusheaRX4OkJ6uTu5nZ3qNUFc0KSQ8DEyW1Fa0XSUfV5opGZmZm/dJjgo+I2ZImAbcBp1QvJDMzGwil6uB/ExH/JOm2iHikmkGZmVn/laqiOVjSq4E3S1pAUjWzW0QsrWhkZmbWL6US/KeBjwOTgS912RZ0HnzMzMwGmVJ18NcD10u6MCIuqWJMZmY2AErVwU+LiJXAzyUd3XW7q2jMzAa3UlU0HyEZxveL3WxzFY2Z2SBXqormrPTfE6sXjpmZDZTeJt0eD/wrMC1dtQL4YUQ8U+nAzMysf3ocqkDSdGA5cAzwIPAQ8DJguaRpPR1nZmaDQ6kS/CXAuRHxP8UrJb0NuBR4WyUDMzOz/ik12Ng/dE3uABFxAzCjciGZmdlAKJXgN5W5zczMBoFSCf5ASR/p5nE+4NmxB6kFCxZw+OEzGDZsOLncDN72tgW0tfV+nJnVn1J18N8GxvSw7aoKxGL9tGDBAj760Qs44oirecMbTmDDhkXcdNMcnnoKLr98Ns0e4NlsSFFE1DqG3WbOnBmLFy+udRh7rRkzZnDMMV9j3LgTyeWSdX/72+3ccsvZnH32cubNq2l4ZlYBkpZExMzutpWqorG9zIoVK5BOoLFxz7qmphPI51fQ3l67uMysNpzg68j06dOJWERHx5517e2LyOWm09RUu7jMrDZ6TfCShlcjEOu/Cy64gN/8Zg4rVtzOunXbWbXqdm68cQ5HHHEBLS21js7Mqq3kUAWp/5N0PfCdiLi/0gFZ+WbPng3AhReeza9/vYLGxumcdNKlXHSRb7CaDUW93mSVNAZ4F/AekhL/NcCPImLDQAfjm6xmZn3Tr5usEfFsRHw7Il4BfAy4CHhC0nclHT7AsZqZ2QDJVAcv6RRJPwW+SjI+/GHAz4BbKxyfmZmVKUsd/EPA7cAXIuIPReuvl/SqyoRlZmb9lSXBN0fExu42RMQ5AxyPmZkNkCzt4C+XNK7wRFJO0jWVC8nMzAZClgTfHBHrC08iIg+8tGIRmZnZgMiS4IdJyhWeSDqAbFU7ZmZWQ1kS9ReBP6SdnQDeTjKjk5mZDWK9JviIuE7SEuBEQECLe7SamQ1+WataVgL5wv6SmiKi5PiEkp4HXAdMAnYBV0bEV/sRq5mZ9UGvCV7S2SS9V58EdpKU4gPobXSTHcD5EbE0He5giaRfufRvZlYdWUrw5wIvjoin+3LiiHgCeCJdflbSCuBQwAnezKwKsrSieRTo6HWvEiRNIWla+cdutp0labGkxWvXru3PZczMrEiWEvwqYKGknwPbCisj4ktZLiBpNHADcF53I1BGxJXAlZCMJpnlnGZm1rssCb49fYxMH5lJGkGS3H8QEa19D8/M9nZta9poXdlKe0c7TY1NtExroXmSJyiohsyTbkvaPyI2ZT6xJOC7wDMRcV6WYzwevFl9aVvTxvy75pNryNHY0EjH1g7yW/PMPX6uk/wA6dd48JKOl3Q/sCJ9fqSkKzJc95XA6cBrJP05fZzcl8DNbO/WurKVXEOO3KgcwzSM3KgcuYYcrSv9g74aslTRfAV4PXAzQET8JcswwRGxiKRJpZkNQcmP+J5dzMXdrs9aq2C9y9TRKSIe7fLH2lmZcMysXkQE8xbOI78lT27U7uGsyG/Jc9nJlzmRV0GmZpKSXgGEpJGS5pJW15iZldIyrYX81jz5LXl2xS7yW/Lkt+ZrHdaQkSXBvx/4EEknpdXAUcAHKxiTmdWJ5knNzD1+LrlROVZvWE1uVI65x8+tdVhDRpYqmhdHxGnFKyS9Evh9ZUIys3rSPKnZLWZqJEsJ/msZ15mZ2SDSYwle0vHAK4CJkj5StGksMLzSgZmZdccdp7IrVYIfCYwm+RIYU/TYAJxa+dDMzDordJzKb8kzeexk8lvyzL9rPm1r2mod2qDUYwk+Iu4A7pB0bUQ8UsWYzMy6VdxxCtj9b+vKVpfiu5HlJutmSV8AXgI0FFZGxGsqFpWZWTfaO9qZPHZyp3WNDY20d5Scf2jIynKT9QckMzq9ALgYeBi4t4IxmZl1q6mxiY6tnUcv79jaQVNjU40iGtyyJPjxEXE1sD0i7oiIM4HjKhyXmdlz9NRxqmVaS61DG5SyJPjt6b9PSHqjpJcCk0sdYGZWCT11nHL9e/ey1MF/VlIjcD5J+/exwHmVDMrMrCfuOJVdrwk+Im5JFzuAEwEknVfBmMzMbABkqaLpzkd638XMzGqp3ATvcd7NzAa5chO8B3I2MxvkSo1F8yzdJ3IBoyoWkZmZDYhSQxWMqWYgZmY2sMqtojEzs0HOCd7MrE5lmnTb9i5tbdDaCu3t0NQELS3Q7H4hZkOOS/B1pq0N5s+HfB4mT07+nT8/WW9mQ0uvJfgeWtN0AIuB8yNiVSUCs/K0tkIulzxgz7+trUkp3qV7GyzmLZznWZkqLEsVzZeAx4EfkjSRfBcwCXgAuAaYVangrHdS9j5nF1+8Z/nCC2N36X7uXCd5q57C7EtdZ2XyoGEDL0uCf0NEvLzo+ZWS7o6Iz0j6ZKUCs2wiOv+4mjcvqZa57DJx0UXJtnx+T0m+eLlr6d6sGlpXtgJ4VqYqyFIHv0vSOyQNSx/vKNrmHq2DTEtLksQBdu1KlvP5ZH17OzQ2dt6/sTFZb1Yt3c2+5FmZKiNLgj8NOB14CngyXf43SaOA/6xgbFaG5uakygVg9eqklF6ogmlqgo7Ok+HQ0ZGsN6uW7mZf8qxMldFrgo+IVRHx5oiYEBET0+X/i4gtEbGoGkFa3xSqW665JqmyKTwvlO7z+eeW7s2qpTD7kmdlqrwsrWgmAu8FphTvn07dZ3uRQum+uBXNnDmuf7fqKtSz50bldreimfPSOa5/rwB1vUn3nB2kPwC/A5YAOwvrI+KGgQ5m5syZsXjx4oE+bd0q1eRRSm6yujmkDUaSntNAwMojaUlEzOxuW5ZWNPtFxP8b4JisnwodmnK5zh2aCvXv8NzOTm4OaZVwwAEHkC/c2e+DvjTxzeVyPPPMM32+xlCXJcHfIunkiLi14tFYZqU6NBX85S/JTdTGRjj0UDeHtMrI5/MVL4335cvA9siS4M8FPilpG7CdpLNTRMTYUgdJugZ4E/BURMzod6TWSXt7UjovVmjyWCjobNkCY8cm/y5bBps2VT9OM6udLJNulzsu/LXA14HryjzeipTzM/ieezqXepYtg1IFIf8MNqsvpWZ0mhYRKyUd3d32iFha6sQRcaekKf2Mz1JdfwYX18E3NiZVMfl8Us/+6U/DqlUwahQ0NMDWrUkp/rDD4MYbe76Gfwab1ZdSJfiPAGcBX+xmWwCvqUhElkmpJo9HHQX77QePP76nDv6FL4SpU2sdtZlVU6kp+85K/z2xkgFIOovki4Qmd6nsk+bm7m+atrQkpfsjj+xcuneHJrOhpdeerJLeLmlMuvwpSa2SXjpQAUTElRExMyJmTpw4caBOO6QVSve53HOHKzCzoSNLK5oLI+Inkk4AXg/MB74JvLz0YVZLPZXuzWzoyDLYWKH36huBb0TETcDI3g6StAC4C3ixpNWS5pQfppkNRW1r2pi3cB6QTBBSGEvessmS4B+T9C3gHcCtkvbNclxEzI6IgyNiRERMjoir+xusmQ0dbWvamH/XfPJbkubBhYlBnOSzy5Lg3wHcRjLxx3rgAOCjlQzKzKx1ZSu5hlyniUFyDbndE4ZY77Ik+IOBn0fEQ5JmAW8H7qlkUGZm7R3tNDZ0nqHGE4P0TZYEfwOwU9LhwNXAC0jmZzUzq5imxiY6tnaeocYTg/RNpin7ImIH0AJ8JSI+TFKqNzOrmJZpLeS35jvVwXtikL7J0kxyu6TZwL8Db07XjahcSGY21LStaaN1ZevuCUBaprXQPKmZucfP7TRJtycG6ZssE34cAbwfuCsiFkh6AfDOiPj8QAfjCT96NlATJPQ2SYgnYbC+6u/nptBaJteQo7GhkY6tHeS35pl7/NzdydyfzZ6VmvAjS3PH+4G5wDJJM4DVlUjuVnmFAcq6TgTS5lZnVkPFrWWGaZhbywygLEMVzAIeAi4HrgAelPSqyoZllVA8SciwYXuWW/3/yGrIrWUqJ0sd/BeB10XEAwCSpgILgGMqGZgNvFKThJjVSlNjE/kt+d3t3cGtZQZKllY0IwrJHSAiHsQ3WfdKTU3JyJLFOjqS9Wa1UtxaZlfscmuZAZQlwS+RdLWkWenj28CSSgdmA6+lJal3z+dh1649yx5G2Gqp0FomNyrH6g2ryY3KdbrBauXL0opmX+BDwAkk87HeCVwREdsGOhi3oulZX1oRlGop41Y0NtCq8bnxZ7NnpVrRlEzwkoYBbdWaNNsJvmelPuDFSXvffeHRR5MZnLpO5dfb8MH+T2TlyPK56amd+0BeY6gqu5lkROwC/iLJtbSDVNemj0uXwl//Ctu2uaWMDQ7Fo0JOHjvZo0JWUZZWNAcD90m6B9hUWBkRp1QsKsusuOkjwN//DmPGwMqVMGlSss4tZayWuhsVsrDe9eyVlSXBX1zxKKwsbW1w000QAePGwfTpSTLfvLlzaxm3lLFaau9oZ/LYzu1z3c69OnpM8OnokQdFxB1d1r8KeKzSgVlphaqZkSOTBL9lC/zhD/DiF8OyZTB2bNJSplAHP8fzaVkNtK1pY1V+FXevvpsD9z+Q6ROmc9Dog9zOvUp6vMkq6RbgkxHR1mX9TOCiiHhztwf2g2+y9qzrTaZ585LEvW0b3HUXNDQkiV6CAw+EQw9Nqmu6tpTpyzXMMpnX2Ps+A3Kdjt73GYLKakUjaXlPrWckLYuIfxjAGAEn+JL8n8gGqZ4KBvMWztvdQ3XNxjWsXLeSpzY9xYH7H8hlb7jMrWgGSKkEX6oOvqHEtlH9C8n6Shdv6LYEn9vTu3v383nzyryGRJR5rFlXxXXvk0ZPYtLoSeyKXazesNo3V6ukVDPJeyW9t+tKSXNwT9aaK7dXaltb8gVw5pnJvx5J0irFMzLVXqkS/HnATyWdxp6EPhMYCby1wnFZL5qbk85Lxb1S58zpua69rQ2uuAJ+/vPkC2H0aLjvPli8GD73uWx19GZ90TKthfl3zQfoNM77nJf6jn+19FiCj4gnI+IVJM0kH04fF0fE8RGxpjrhWSnNzUmJvakpSfKtrd2XyAstbhYtSlrb7NiRlPY3b046RV1xRfVjt/rnMWZqr9d28BFxO3B7FWKxPiok7lyu8wQeXYclKHSGWrcuGcpgxAjYvh2efTZpcXP33bV7DVbfmic1O6HXUJbRJG2QyjqBR3t70gEKkqaUAPvsA1u3JstS9WI2s+pxgt+LFSfugu6GJSiMAz95ctI2fvv2pJpm+PCkFH/ccdWL2cyqxwl+L5Z1Ao9Ci5upU+GAA5Lkvnlzsnz44fCBD1QvZjOrnl7Hg68md3TqmapQjzJyZI5PfOKZzD1fzcDjwddauR2dbBDJMhb8yJHw2GNw2GFw2WXi2GODp5+Gk06CD36w+6RdfKO2sbHnG7VmldTf8eKtey7B15F58+Chh+CBB2DJEnHIIUEuB/vvnwxC1l3SrkSPWBtaeitdF5L3n5/4M+u3rWfcvuM46uCjdifxwnjxuYZcp/byxU0qXYLvWdkTftje5Y474Fe/gvvvT55v3w5PPZXUy/c06UfWG7Vm5Sgk74eefohV+VWs37KeVetX8eC6B3dP+lE8XvwwDSM3KkeuIUfrSs9S019O8HWirS1J7BFJk0mATZtg586kOWRPSTvrjVqzchSS92PPPsaoEaMYN2oco/YZxeMbH9+dxNs72mls6FzK8HjxA8N18HWitRX22y9pHTOs6Gt761aYMKHnpN3SktS5Q+c5XD1+vPVFOY0A7uCOXvf5Dt8BIFdch2iZuQ5+L1du65riv3vxjdq+jB9vVkrbmjZO/+npCJHflmfb9m0EwYH7H8jyDy3nnFvPITcqt3vMmh07d/DYxsdYu2ktI4aN4MJXXcipLzm11i9j0KtZKxpJbwC+CgwHroqIz1fyekNRIVEXbrAuX55M/rFjR9Ki5u9/h3e+M2nr3lPSbm52QreB17qylZdMfAn3rbuPEcNGsO7v69ixawcbt28E2D3wWPOkZk6ZegqX/O4Stu/czsT9JjJ57GRufvBmpo6f6tY0/VCxBC9pOHA58FpgNcnwwzdHxP2VuuZQVqhqmTEDVq+GtWvhoIPgwgvhVBeCrAbaO9o5fPzhbN+1nTvb72Tk8JHsM2wfdu7aCcApU0/ZnbyXr13OrOfP2j0hN0B+S94Tc/dTJUvwxwL/FxGrACT9CHgL4ARfAcXDB++7L5x4oqtarLaaGpvIb8mzdvNapjROYdSIUWzZvoVRI0ZxB3ewfO1yTiUpfXhi7sqoZII/FHi06Plq4OUVvN6Q56oWG0wuPvHiXrdfzHP3uej2iwBPDjIQKpngu7v795w7upLOAs4CaHLbPLO6ERG0rWnjnP89Z/dcrNMnTOeg0Qftnqt13qx5AJ06O+2KXZ4cZIBUsh38auB5Rc8nA4933SkiroyImRExc+LEiRUMx8yqrXlSM5e94TJmHjKTIw86kon7TyS/JU9+a56WaS2d9vPkIAOvYs0kJe0DPAj8E/AYcC/wrxFxX0/HuJmkWX3yWDOVU5NmkhGxQ9J/AreRNJO8plRyN7P65ZmdaqOi7eAj4lbg1kpew8zMuuexaMzM6pQTvJlZnXKCNzOrU07wZmZ1alCNJilpLfBIreOoExOAdbUOwqwH/nwOnOdHRLediAZVgreBI2lxT21jzWrNn8/qcBWNmVmdcoI3M6tTTvD168paB2BWgj+fVeA6eDOzOuUSvJlZnXKCNzOrU07wg5SkjX3Yd5akVwzgtadIWj5Q57P6Iykkfa/o+T6S1kq6pY/nWShpZrp8q6Rxvez/sKQJZQU9BDnB14dZQJ8SfDpev1m5NgEzJI1Kn7+WZN6HskXEyRGxvr+B2R5O8HsRSW+W9EdJf5L0a0kHSZoCvB/4sKQ/S/pHSRMl3SDp3vTxyvT4eZKulPRL4Lq0pP47SUvTx4D9CrAh4RfAG9Pl2cCCwgZJ+0u6Jv38/UnSW9L1oyT9SFKbpB8Do4qO2V06l3SjpCWS7kun9bQyuBS3d1kEHBcRIek/gI9FxPmSvglsjIj5AJJ+CHw5IhZJaiKZdGV6eo5jgBMiYouk/YDXRsRWSS8i+Q/q3oWW1Y+AT6fVMs3ANcA/ptsuAH4bEWem1S73SPo18D5gc0Q0S2oGlvZw7jMj4pn0F8K9km6IiKcr+mrqkBP83mUy8GNJBwMjgb/1sN9JwBHS7nnPx0oaky7fHBFb0uURwNclHQXsBKZWJGqrSxHRlv6CnM1zJ/Z5HXCKpLnp8wagCXgVcFnR8W09nP4cSW9Nl58HvAhwgu8jJ/i9y9eAL0XEzZJmAfN62G8YcHxRIgcgTfibilZ9GHgSODI9ZuvAhmtDwM3AfJL7QOOL1gt4W0Q8ULxz+hks2fkm/WyfRPIZ3ixpIckXhPWR6+D3Lo3suZH17qL1zwJjip7/EvjPwpO0hN7T+Z6IiF3A6SRz55r1xTXAZyJiWZf1twFnK83okl6arr8TOC1dN4OkaqerRiCfJvdpwHEViXwIcIIfvPaTtLro8RGSEvtPJP2OzkOt/gx4a+EmK3AOMDO9kXU/yU3Y7lwBvFvS3STVM5t62M+sWxGxOiK+2s2mS0iqANvSJreXpOu/AYxOq2Y+BtzTzbH/C+yT7nMJcPfARz40eKgCM7M65RK8mVmdcoI3M6tTTvBmZnXKCX4QS3v93SGp7NYtkqZJukvStqI2yYVt50panvYWPK9o/XxJr+lH6FbnKv3ZTLcPT3vB3lK0zp/NPnCCH9zOBFojYmc/zvEMSaua+cUr0yZq7wWOJWkH/6a0Nysk7e0/3o9rWv2r2GezyLnAii7r/NnsAyf4we004CYlo0UulHS9pJWSflBoX9ybiHgqIu4FtnfZNB24OyI2R8QO4A7grekxjwDjJU0awNdi9aWSn00kTSYZ5+aqLsf4s9kHTvCDlKSRwGER8XC66qXAecARwGFAYQCxL6ft37s+eivlLAdeJWl8OibNySRdwguWFq5hVqwKn02Ar5C0k9/VzTZ/NjPyUAWD1wRgfdHzeyJiNYCkPwNTgEUR8eFyTh4RKyT9F/ArYCPwF2BH0S5PAYeUc26rexX9bEp6E/BURCxJhy3oyp/NjJzgB68tdB5/Y1vR8k7Sv52kLwMndnP8jyLi86UuEBFXA1en5/kcsLpoc0Mag1lXlf5svpJkoLKT0+uMlfT9iPi3dLs/mxk5wQ9SEZFPWxGUHGSp3FISgKQDI+KpdEjhFuD4os1TgZ+Ue26rX5X+bEbEJ4BPwO6Bx+YWJXfwZzMzJ/jB7ZfACXSuOumT9GbUYmAssCttDnlERGwAbpA0nuQm14ciIp8eMwI4PD3OrDuV/mz2dIw/m33gsWgGsXQEvo9ExOlVvu5bgaMj4sJqXtf2Hv5s7h3cimYQi4g/Abf3pzNJmfYBvljla9pexJ/NvYNL8GZmdcoleDOzOuUEb2ZWp5zgrWYkXZAOdNaW9nB8ebr+vLR3bWG/WyWNG6Brbuxl+zhJHyx6foik6wfo2gslzRyIc5ll4QRvNSHpeOBNJC0imkkmWX403XwesDvBR8TJEbG+SqGNA3Yn+Ih4PCJOrdK1B0QNbnzaIOUEb7VyMLAuIrYBRMS6iHhc0jkk3dBvl3Q7gKSHJU2QNCUd0OqqdJjjH0g6SdLvJT0k6dh0/3nFw8+m+04pvrik0ZJ+I2mppGWS3pJu+jzwwvQXxRfSay5Pj2mQ9J10/z9JOjFdf4akVkn/m8bx31nfhPT8v0vjWCrpFen67xXFRPpaT0k7GH1B0r3pL5/3pdtnSbpd0g+BrhNg21AVEX74UfUHMBr4M/AgyeTfry7a9jAwoetzkjFOdgD/QFI4WQJcAwh4C3Bjuv88kt6PheOXA1PS5Y3pv/sAY9PlCcD/peeZAiwvOnb3c+B84Dvp8jSgnaTb/BnAKqAxff4I8LxuXvNCYGaXdfsBDenyi4DF6fKri15PI/C3NOazgE+l6/cl6fDzAmAWyaTpL6j139aPwfNwCd5qIiI2AseQJKy1wI8lnZHh0L9FxLKI2AXcB/wmIoKk1DqlDyEI+JykNuDXwKHAQb0ccwLwvTT+lSSJfGq67TcR0RERW4H7gednjGME8G1Jy0i63x+Rnv8O4HBJBwKzgRsiGdb5dcC/p4N6/REYT/LFAMmgX3/LeF0bAjxUgdVMJJNFLAQWpgnu3cC1vRxWPLDVrqLnu9jzed5B5+rH7sZMOQ2YCBwTEdslPdzDfsVKjXPe7YBbGXwYeJJk0pVhwNaibd9L43wXyQQbhRjOjojbOgWWjNmyKeM1bYhwCd5qQtKLtWcGKYCjSErEAM8CY/px+oeBo9PrHE1ShdFVI8mQtNvTuvRCibvUte8kSbhImgo0AQ/0I85CHE+kv0hOB4pvkF5LcsOZiLgvXXcb8IF0TBYkTZW0fz9jsDrlErzVymjga2nzxx0kdeBnpduuBH4h6YmI6G642d7cwJ5qjHtJ6vm7+gHwM0mLSe4FrASIiKfTm7bLgV8AlxcdcwXwzfTXxg7gjIjYpmwTGBX8XFJhBqO7gE+SDPr2duB2ikrhEfGkpBXAjUXHX0VSFbVUyYXXAv/SlwBs6PBQBWaDVNoXYBlJU9KOWsdjex9X0ZgNQpJOIvlV8TUndyuXS/BmZnXKJXgzszrlBG9mVqec4M3M6pQTvJlZnXKCNzOrU07wZmZ16v8DymZYVbWCsXUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Lump layers together\n",
    "# latency plots\n",
    "plt.clf()\n",
    "i = 1\n",
    "xtick_labels = []\n",
    "df_datasets = []\n",
    "df2 = df\n",
    "df_medial3 = df2[df2[\"Most Medial Barrel\"] == 3]\n",
    "df_medial1 = df2[df2[\"Most Medial Barrel\"] == 1]\n",
    "    \n",
    "medial_lat = []\n",
    "lateral_lat = []\n",
    "medial_lat.append(df_medial1[(df_medial1['abs(delta 2)'] > 0) \n",
    "                               & (df_medial1['Stim Barrel'] == 3)]['abs(delta 2)'])\n",
    "lateral_lat.append(df_medial3[(df_medial3['abs(delta 2)'] > 0)\n",
    "                               & (df_medial3['Stim Barrel'] == 3)]['abs(delta 2)'])\n",
    "medial_lat.append(df_medial3[(df_medial3['abs(delta 2)'] > 0)\n",
    "                               & (df_medial3['Stim Barrel'] == 2)]['abs(delta 2)'])\n",
    "lateral_lat.append(df_medial3[(df_medial3['abs(delta 1)'] > 0)\n",
    "                               & (df_medial3['Stim Barrel'] == 2)]['abs(delta 1)'])\n",
    "medial_lat.append(df_medial3[(df_medial3['abs(delta 1)'] > 0)\n",
    "                               & (df_medial3['Stim Barrel'] == 1)]['abs(delta 1)'])\n",
    "lateral_lat.append(df_medial1[(df_medial1['abs(delta 1)'] > 0)\n",
    "                               & (df_medial1['Stim Barrel'] == 1)]['abs(delta 1)'])\n",
    "medial_lat.append(df_medial1[(df_medial1['abs(delta 1)'] > 0) \n",
    "                               & (df_medial1['Stim Barrel'] == 2)]['abs(delta 1)'])\n",
    "lateral_lat.append(df_medial1[(df_medial1['abs(delta 2)'] > 0) \n",
    "                               & (df_medial1['Stim Barrel'] == 2)]['abs(delta 2)'])\n",
    "\n",
    "medial_lat = pd.concat(medial_lat)\n",
    "lateral_lat = pd.concat(lateral_lat)\n",
    "\n",
    "#dfd1_ant = df2_ant['abs(delta 1)'].dropna()\n",
    "#dfd2_ant = df2_ant['abs(delta 2)'].dropna()\n",
    "\n",
    "#df2_post = df2[df2['Slice Location'] == \"Posterior\"]\n",
    "#dfd1_post = df2_post['abs(delta 1)'].dropna()\n",
    "#dfd2_post = df2_post['abs(delta 2)'].dropna()\n",
    "\n",
    "#dfds_ant = np.concatenate([dfd1_ant, dfd2_ant])\n",
    "#dfds_post = np.concatenate([dfd1_post, dfd2_post])\n",
    "\n",
    "jittered_xs_1 = np.random.normal(i, 0.04, lateral_lat.shape[0])\n",
    "jittered_xs_2 = np.random.normal(i+1, 0.04, medial_lat.shape[0])\n",
    "plt.scatter(jittered_xs_1, lateral_lat, alpha=0.5, color='blue')\n",
    "plt.scatter(jittered_xs_2, medial_lat, alpha=0.5, color='green')\n",
    "\n",
    "xtick_labels.append(\"Lateral\\n(n=\" + str(len(lateral_lat)) +\")\")\n",
    "df_datasets.append(lateral_lat)\n",
    "xtick_labels.append(\"Medial\\n(n=\" + str(len(medial_lat)) +\")\")\n",
    "df_datasets.append(medial_lat)\n",
    "\n",
    "p_val = kstest(medial_lat, lateral_lat).pvalue\n",
    "print(\"KS test p-value \", p_val)\n",
    "\n",
    "# significance: L2/3 Rostral and L4 Rostral, and \n",
    "#               L2/3 Rostral and L4 Caudal\n",
    "#significance_bar(0.75, 2.25, 3.35, \"P<0.05\", fontsize=10)\n",
    "\n",
    "plt.boxplot(df_datasets)\n",
    "plt.xticks(range(1, len(xtick_labels)+1), xtick_labels)\n",
    "plt.title(\"Latency of L4 Barrel Pair Crossings\")\n",
    "plt.xlabel(\"Stimulation Layer\")\n",
    "plt.ylabel(\"Crossing Latency Difference (ms)\")\n",
    "plt.savefig(data_dir + cell_type_prefix + \"_results_latency_medial-bias_\" + str(datetime.date.today()) +\".png\", bbox_inches=\"tight\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "6a525873",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6907000000000001 +/- 0.6193717033600074\n",
      "1.028032 +/- 1.0011922707332492\n",
      "KS test p-value  0.4475725833721871\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEiCAYAAADziMk3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqU0lEQVR4nO3de3xcdZ3/8dcnadqktEkDbe0lhEKhFMQApaIsrBRlV0QX3IiuqOxCu4uXVUSpl98iS6rib12rIooIchPUqovRsizeViksikDLJRSK4GJJK/TKNL3QC20/+8f3THIynSQzSU5mcvJ+Ph7n0TPn+p3J6We+8z3f8/mauyMiIulTUeoCiIhIMhTgRURSSgFeRCSlFOBFRFJKAV5EJKUU4EVEUkoBXoacmX3QzNab2XYzO6TU5UmKmV1oZvcncNztZnbEYB83SWb2LTO7otTlGGkU4MuEma02szP7sd8yM/vHJMqUBDOrAr4C/LW7j3P3zTnrZ5iZm9moPo7zm762i9bviALiJjNbYmYTBuWNDFDsfW6PptVm9ulC9o0+t+eKPN+bzew+M9tmZhvN7F4zO6d/pS+eu3/A3T83VOeTQAFehtqrgGrgyf4ewMzeC/T6BRBzvLuPA44A6oGWfp7zgPP19SVUoAlR+c4H/tXMzhrIwcysMs+y84D/AG4DGgh/g38F/qaHYwzG+5IyoABf5sys3szuimpdmWi+IVp3FfCXwDeiWuA3ouWzzexXZvaSmf3BzN4VO96tZnatmf1XVJt70Mxmxta/OrbvejP7FzObYmYvx5tTzOykqExVeco8xsyuNrMXounqaNks4A/RZlvM7Df9+DzqgCuBTxazn7tvBe4Ejo0d6yIzWxV9Ds+Z2ftj6+aZ2Voz+5SZrQNuMbMWM7vDzL5rZluBC82szsxuMrMXzezPZvb5fEG2gPI9QPjSO87MTjazB8xsS3Tcb5jZ6FjZ3MyOjOZvNbPrzOxuM9sBnJHzeRnhF9Pn3P1Gd+9w9/3ufq+7/1O0zYVm9lsz+6qZvQS0RO/rtuhv/LyZfcbMKqLtj4x+AXREv4x+mD1XdIwN0bo2MzsuVs7P53y2l0XbvmhmF8XKfIiZ/aeZbTWzh6PP9P6+ziF5uLumMpiA1cCZeZYfArwDGAuMJ9TEfhpbvwz4x9jrg4A1wEWEWu4cYBPw6mj9rcBLwMnR+u8BP4jWjQdeBC4j1LLHA6+L1t0NfDB2nq8CX+/hvXwW+D0wGZgE/I4QYABmAA6M6mHfvtZfC3ysr+2ibR04MpqvB34JfDa2/q3ATMCA04GXgTnRunnAXuCLwBighlD7fwV4O6FyVAP8FLg++twnAw8B74+OcSFwf1/vMzr/qdH53wScBLw+WjcDWAVc2sP7uhXoiPavAKpzzjM72v7wXj6nC6P3+pHonDWE2v7S6BqYATwDLIi2XwJcnj0fcFq0/M3ACmBC9J6OAabGyvn5nM/2s0AVcHb03uuj9T+IprGEL+Q12c+xt3NoyvO3LXUBNEV/iB4CfJ7tTgAysdfL6B7g/w74n5x9rgeujOZvBW6MrTsbeDqaPx94tIfz/h3w22i+ElgHnNzDtv8LnB17/WZgdTTfGdh62LfH9cBc4DG6Al8hAX4rsAXYBzwNTO9l+58CH43m5wF7iAVMQoC/L/b6VcBuoCa27Hzgnmj+QvoO8FuADCGIX9LDtpcCP8l5X/EAf1sv7+nUaPvqXra5EGiPva6M3texsWXvB5ZF87cBNwANOcd5I+GL4PVARc66W+ke4HfG/3bAhmi/SsKX6NGxdZ+nK8D3eA5NB05qoilzZjbWzK6PfiZvBe4DJvTSDHAY8Lro5/0WM9sCvBeYEttmXWz+ZWBcNH8oITjnsxQ41kLvjb8COtz9oR62nQY8H3v9fLSs36LmgW8SAvDeInad4+4TCDXN64D/MbPq6JhvMbPfR81RWwhfdhNj+2509105x1sTmz+MUAN9MfZZX0+oyRdqorvXu/sx7n5NVK5ZFpri1kV/8y/klCvXml7WZW9iT+2jHPFjTARGc+DfcHo0/0lC7fkhM3vSzOYDuPtvgG8QfmWtN7MbzKy2p3Ll/B2z1+Ekwhd4vDyd80WeY8RTgC9/lwFHE5pKaoE3RMst+jc3Hega4F53nxCbxrn7Bws41xpCk8UBokD3I8KXxQXA7b0c5wVC8MtqjJYNRC2hBv/DqE384Wj5WjP7y752dvdXgBuBwwnt3GOAHwOLgVdFXwJ30/W5woGfbe6yNYSa7sTYZ13r7q8u7q0d4DrCr42jor/5v+SUq7cy5fpDVM539HHO+DE2EWrRuX/DPwO4+zp3/yd3n0ao2X8ze0/A3a9x95OAVwOzgE/0cd5cGwnNNw2xZYd2K+jAzzFiKMCXlyozq45NowhtoDsJNyUPJtxgjFtP6CGSdRcwy8wuMLOqaHqtmR1TwPnvAqaY2aUWboqON7PXxdbfRvg5fw7w3V6OswT4jJlNMrOJhB4bvW2fz5j4ZwFsI/wKOCGazo62Owl4sK+DRb94LiJ8ls8RaqhjiAKKmb0F+OtiCujuLxLa9b9sZrVmVmFmM83s9GKOk8d4QtPSdjObDRTy5dxTGR34OHCFhZvK2XKeZmY39LDPPsKX+VXRNXBYdIzvApjZOy260U9oXnJgX3Sdvc7CjfcdwC5C01gx5d0HtBJu9I6N3v/fZ9cPxjlGEgX48nI3IQBlpxbgasJNr02EG5c/z9nna8B5FnrYXOPu2wiB6t2EWvM6um4U9ira968I3efWAc8S65Xh7r8F9gOPuPvqXg71eWA50AY8ATwSLSvGdrp/FmdENcd17r6OEJgB1rv7nl6O87iZbScEon8A/tbdX4re6yWEQJYB3kPoZVOsvyd8WTwVHecO+m4O6cvCqDzbgG8DPxzIwdz9DsI9lPmEa2I94e+xtJfdPkIIoM8B9wPfB26O1r0WeDD6XO8kNJv9ifAr69uEz+F5QvPQ4n4U+cNAHeEavJ1QYdgdrRusc4wIFt24ECmIha6N33f3G0tdFhkZzOyLwBR3/4dSl2W4UQ1eCmZmryV0uxxQjVKkNxae42iK+ryfDCwAflLqcg1HemJNCmJm3yH0//5o1LwhkpTxhGaZaYTuk1+m9+Yk6YGaaEREUkpNNCIiKaUALyKSUmXVBj9x4kSfMWNGqYshIjJsrFixYpO7T8q3rqwC/IwZM1i+fHmpiyEiMmyY2fM9rVMTjYhISinAi4iklAK8iEhKKcCLiKSUAryISEopwItI4trWtdGyrIV5l81j8uGTqays5LjjjmPJkiWlLlqqlVU3SRFJn7Z1bSx+YDEbfr+Bx5c8zhmXnMHYmWOZVzGPyz9xOQDnn39+iUuZTqrBi0iiWp9upb66nkd+9AjnfOIcml7fxCHjDmHtIWu56aabuOqqq0pdxNRSgBeRRLV3tFNXXcem9k00vqYRgLrqOto72jnttNNYtWpViUuYXgrwIpKoxrpGOnZ1MLFxIu1PtAPQsauDxrpG7r//fo45ppDRJKU/FOBFJFHNs5vJ7Mow511zuPNLd9L2+zY2b99Mw+YGFixYwOWXX17qIqaWbrKKSKKapjSx8JSFtNa3smvvLu6//n42t2/msWMe46qrrtIN1gQpwItI4pqmNNE0pQnmEcZnkiGhJhoRkZRSgBcRSSkFeBGRlFKAFxFJKQV4EZGUUoAXEUkpBXgRkZRSgBcRSSkFeBGRlFKAFxFJKQV4EZGUUoAXEUkpBXgRkZRSgBcRSSkFeBGRlFKAFxFJKQV4EZGUUoAXEUmpxAO8mVWa2aNmdlfS5xIRkS5DUYP/KLBqCM4jIiIxiQZ4M2sA3grcmOR5RETkQEnX4K8GPgnsT/g8IiKSI7EAb2ZvAza4+4o+trvYzJab2fKNGzcmVRwRkREnyRr8qcA5ZrYa+AHwRjP7bu5G7n6Du89197mTJk1KsDgiIiNLYgHe3f+fuze4+wzg3cBv3P19SZ1PRES6Uz94EZGUGjUUJ3H3ZcCyoTiXiIgEqsGLiKSUAryISEopwIuIpJQCvIhISinAi4iklAK8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISinAi4iklAK8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISinAi4iklAK8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISinAi4iklAK8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISinAi4iklAK8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISiUW4M2s2sweMrPHzexJM1uU1LlERORAo/rawMwqgOOBacBO4El3X1/AsXcDb3T37WZWBdxvZj9z998PqMQiIlKQHgO8mc0EPgWcCTwLbASqgVlm9jJwPfAdd9+fb393d2B79LIqmnzwii4iIr3prQb/eeA64P1RsO5kZpOB9wAXAN/p6QBmVgmsAI4ErnX3BwdcYhERKUiPAd7dz+9l3Qbg6r4O7u77gBPMbALwEzM7zt1Xxrcxs4uBiwEaGxsLK7WIiPSpz5usZvZOMxsfzX/GzFrNbE4xJ3H3LcAy4Kw8625w97nuPnfSpEnFHFZERHpRSC+aK9x9m5mdBryZ0CRzXV87mdmkqOaOmdUQ2vKfHkBZRUSkCIUE+H3Rv28FrnP3pcDoAvabCtxjZm3Aw8Cv3P2u/hVTRESK1Wc3SeDPZnY9oQb+RTMbQwFfDO7eBpw4wPKJiEg/FVKDfxfwC+CsqC39YOATSRZKREQGrpCa+MvAPUBNdHN1KrAp6YKJSPosWbKEI2cfSUVlBfWH1fOOK95B27q2UhcrtQp5kvVzwIXA/9L1oJIDb0yuWCKSNkuWLOETn/4Exy44lrNOPIutz25l6ZeWsmHHBq795LU0TWkqdRFTp5A2+HcBM919T9KFEZH0uuqqq3jTpW9iwuwJ1NfUM3HORN7+ybdz19fuovWcVgX4BBTSBr8SmJBwOUQk5VatWoU1GnXVdZ3LGl/TSGZNhvaO9hKWLL0KqcH/f+BRM1tJSCAGgLufk1ipRCR1jjnmGLzd6RjbQX1NPQDtT7RTf2g9jXV6ij0JhQT47wBfBJ4A8iYWExHpy+WXX97ZBj/rxFmdbfBN5zfRPLu51MVLpUIC/CZ3vybxkohIqp1/fkhvdcWiK/jvRf9NXUMdZ/7jmVz5z1eq/T0hlpMo8sANzL5CaJq5k+5NNI8MdmHmzp3ry5cvH+zDioiklpmtcPe5+dYVUoPPPo36+tgydZMUESlzfQZ4dz9jKAoiIiKDq8dukmb2vmi4vp7Wz4wyTIqISBnqrQZ/CKF75ArCqEzZIfuOBE4npCv4dOIlFBGRfultRKevmdk3CG3tpwJNhEG3VwEXuLueTBARKWO9tsFHQ+79KppERGQYKSRVgYiIDEMK8CIiKVVIuuDKqKlGhrG2NmhthfZ2aGyE5mZo0sODIqlWSA3+j2b2JTM7NvHSSCLa2mDxYshkoKEh/Lt4cVguIulVSIBvAp4BbjSz35vZxWZWm3C5ZBC1tkJ9fZgqKrrmW1tLXTIRSVIhT7JuA74NfNvM3gAsAb5qZncAn3P3PyZcRumFmfVrvyuucNrV0VUk1QpqgwfeClwEzAC+DHwP+EvgbmBWguWTPvSULM7MOte1tIRmmfr6rvWZTGiLF5H0KqSJ5lngXOBL7n6iu3/F3de7+x3Az5MtngyG5uYQ0DMZ2L+/a75ZKbhFUq2QbJJN7r493wp3v2SQyyMJaGqChQu796JZsEC9aETSrpAAf62ZfdTdtwCYWT3wZXefn2jJZFA1NSmgi4w0BfWiyQZ3AHfP0JUjXkREylQhAb4iqrUDYGYHU1jNX0RESqiQQP1l4HdRt0iAdwJXJVckERmO2ta10fp0K+0d7TTWNdI8u1ljrZZYnzV4d78NOA9YD2wAmt399qQLJiLDR9u6NhY/sJjMzgwNtQ1kdmZY/MBi2tbpcelSKrSp5Wkgk93ezBqVD15EslqfbqW+up76mtCam/239elW1eJLqJAHnT4CXEmowe8DjDDotv5qIgJAe0c7DbUN3ZZdc/Y1ACxiUVHH6unhPSleITX4jwJHu/vmpAsjIsNTY10jmZ2Zzpo7wCV3X0J9TT0t81oO2D7+pLUkp5BeNGuAjqQLIiLDV/PsZjK7MmR2Ztjv+8nszJDZlaF5th6XLqVCavDPAcvM7L+A3dmF7v6VxEolIsNK05QmFp6ysFsvmgUnLlD7e4kVEuDbo2l0NBXEzA4FbgOmAPuBG9z9a/0ppIiUv6YpTQroZaaQdMGLAMzsIHffUcSx9wKXufsjZjYeWGFmv3L3p/pZVhERKUKfbfBmdoqZPQWsil4fb2bf7Gs/d3/R3R+J5rdF+08fYHlFRKRAhdxkvRp4M7AZwN0fB95QzEnMbAYhf82DedZdbGbLzWz5xo0bizmsiIj0opAAj7uvyVlU8CDcZjYO+DFwqbtvzXPsG9x9rrvPnTRpUqGHFRGRPhRyk3WNmf0F4GY2GriEqLmmL2ZWRQju33N3jQAqMkLl5qmRoVFIDf4DwD8T2s/XAicAH+prJwuDhd4ErFKXSpGRK1+emuxySVYhAf5od3+vu7/K3Se7+/uAYwrY71TgAuCNZvZYNJ09oNKKyLATz1NTYRXd8tRIsgppovk6MKeAZd24+/2EvDUiMoLly1OTXS7J6jHAm9kpwF8Ak8zs47FVtUBl0gUTkXTIl6cmu1yS1VsTzWhgHOFLYHxs2krIDy8i0qd8eWqyyyVZPdbg3f1e4F4zu9Xdnx/CMolIiuTLU5NdLsmyvlJ2mtkk4JPAq4Hq7HJ3f+NgF2bu3Lm+fPnywT7siKR0rFLOdH0OHjNb4e5z860rpBfN9wgjOh0OLAJWAw8PWulERCQRhfSiOcTdbzKzj8aabe5NumBSmLY2aG2F9nZobITmZmjSL18ZJjRQd7IKqcG/Ev37opm91cxOBA7s8yRDrq0NFi+GTAYaGsK/ixeH5SLlTgN1J6+QGvznzawOuIzQ/70WuDTJQklhWluhvj5M0PVvq54fkWFAA3Unr88avLvf5e4d7r7S3c9w95OAmUNQNulDezvU1XVfVlcXlouUu/aOduqqu1/AddV1egBqEBWUTTKPj/e9iSStsRE6ckbL7egIy0XKXWNdIx27ul/AHbs69ADUIOpvgFcKgjLQ3Bza3TMZ2L+/a75Zz4/IMKCBupPX3wCvDqxloKkJFi4Mbe9r14Z/Fy5ULxoZHrIPQNXX1LN261rqa+pZeMpCtb8Pot5y0WwjfyA3oCaxEklRmpoU0GX40kDdyeotVcH4oSyIlI760oukU3+baCQl1JdeJL0U4Ee4eF/6ioquefWlFxn+FOBHOPWlF0kvBfgRTn3pRdJLAX6EU196kfRSgE+h7A3S+fOhpaX3G6bqSy+SXoUkG5NhJNsrBrr3iuktaKsvvZQrpRMeGNXgUybbKwbUK0bKW8uyFuYvnU/Lspa8KYKVTnjgFOBTptBeMW1tofmmkGYckcGUDdB9Be54OuEKq6C+pp766npan1ZtpVBqoilDA3mytLERnnkmzC9dGoL7tGkwa1b34y9eHGr2hTbjiAyWbIDuKw98e0c7DbXdxxZSOuHiKMCXmYEG3+OOg9tuC/Pjx8OWLeGLIt4rpreBQhTgpVgHH3wwmUym6P0WnbHowGUcuAygenw1n7rzU0onXCQF+DIz0OC7ciWccgr86U+wbRtMmACvfnVYft55YZv29vDlEaeHm6S/MpkM7oUnmG1Z1kJmZ4Y9+/awatMqOnZ1MLpyNHOmzuG6t13XuV22Db6+up5rzr6mM53wghMXJPE2Uklt8GVmoE+WtrfDzGi8rXPPhXnzwuv4/nq4SUqpeXYzz2WeY9nqZby852WqKqvYunsra7au6dYOH08nDCidcD8owJeZgQbfQvbXw01SSk1TmpheO53aMbW8sv8VxlaN5fQZpzOzfuYBN1CbpjTRMq8FgJZ5LQruRbJiflolbe7cub58+fJSF6Ok4m3wdXUhOGcysHTpwWzdWnw7ZzHq6+t56aWXEj2HpI+ZFdVEAzB/6XwaahuosK465n7fz9qta7n53JsH5RwjhZmtcPe5+dapDb7MZJ8sjfeiWbAAbr+98HbO/vbCMdNIjDI0GusayezMdDa/gMZjTYICfBkq9MnSngJ5dsquv/pqDeQh5aV5djOLHwiPXNdV19Gxq0M3UBOgNvhhqq+BOjSQh5Qzjcc6NFSDH6b66k6pvu5S7goZjzWbiwZC90rloilOYjV4M7vZzDaY2cqkzjGS9dWdUgN5yHAXz0UDKBdNPyTZRHMrcFaCxx/R+uoOqb7uMtzFc9EAykXTD4kFeHe/D1Cfu4T01Zddfd1luGvvaKeuuvvPUOWiKY7a4IepnrpTZtvX+1ovkpTByuGurpQDV/IAb2YXAxcDNKr9oCh9dafUQB4y1OL5Y7KpgC//zeVMr53Onn17igr48a6UgHLR9EPJu0m6+w3uPtfd506aNKnUxRGRAcjN4b5n3x7++NIfefTFR4setEO5aAau5AFeRNIjt9181aZVjB89nj379vRr0A7lohmYxJpozGwJMA+YaGZrgSvd/aakzicipZfbbt6xq4OqyqpuQT/fjVKNvZqMJHvRnO/uU929yt0bFNxF0q95djOZXRkyOzPs9/2MrhzNtt3bmD1xduc2uTdKNfZqctREIyKDJjcFwZypc5h58EzGVI5hv+/vvFHaPLurv67GXk1OyXvRSPEGMmarSNJyUxDkNr8sOHGBxl4dIgrww4wGzJbhpq+cM/F2+/Xb17Nq0yo27NjA5IMm07auTW3xA6AmmmEmnkSsoqJrvlW/ZmWYyrbbP7v5WX7b/lu27NzCqIpRTBs3TW3xA6Qa/DDT04DZjz0GLS1qtpHS6W9PmGy7/SU/v4S9vpfJYycze+JspoybQmZnRm3xA6Ah+4aJ7JBlLS2hWaa+6+ltnnkGnnwyDLAdH+av2GYbDYsm/dJS1/c2g3Kejr63GYE0ZF+KNDeHNnfoCuZPPgnHHafc71IatmgrV95z5QF5Y7Kvsw8qZcXTGcRHczqo6iBGV47Oe4xFZyzCux9GCqA2+GEmm0Ssvh7Wrg3/Hn44zJzZfTvlfpehVEzmx566RTrerQ99vi6VUhzV4Ieh3CRi+ZptlPtdhlIxmR976ha5dutaFp6ysNculVIcBfgUyNdsk8mE9MAiQ6GYQbR7+zIoZBg/KZyaaFIgX7ON+sXLUCpmEO3cdAZqikmOetEME/3t4ZJ96vWxx2DLFpgwAU44IX83SvWikf7oz3VTbJdKXZs9Uy+aFCokXUH2qde9e+G558KDUS+9BAcdpKdfpbQKbYrJfhEAtCxrUZbJIqmJZhjKBu5Mpnu6gracB/6yT72+8ALU1ITae00N/PnPevpVyl88yySgLJP9oAA/DBWarqC9veuma3V1WFZdHV6rG6WUu3h3SkBZJvtBTTTDhF9Z2/nEYAtAfc4G2dctXYtujnqoffbcHg7a2H17v7J2YIUUGUTtHe1UVVSxbPUyAJatXsbRhxytLJNFUIAfJmzR1s6bTPn6vWdft7R0LYu3wa9cGWr7+/fDa14DlZUHtsGbmZ4WlJKK33xtW9/G5h2bmTxuMgA7X9nJfc/fx7wZ80pbyGFETTTDUHNzCOiZTAjY2fnmnF5m2e6Ts2bBEUeENvgjjoCjjtINVik/uSM77d67m5d2vcSOPTu6beeoN02hVIMfhrKBO96LZsGC/AE796lXkXKV2+Y+qmIUh9YeyvZXtgNQU1XDCVNOYM++PaUs5rCiAD9MKXBL2uSmMKirruPlPS9TVVnF8zzPvBnzyOzMMLVmaglLObwowItISeQ+7DSmcgwduzo6a/DHTDyGZauXUTsm3PzPPvGaL/2B5KcAn2Iau1XKVTxlcFVFFT979me8sO0FaqpqeO201zLz4JmMrhzNkQcfyfTa6TzEQ9TX1Cv5WJGUqmCYKPZR7fjYrYUOAqLHwaU/+nPdtCxrIbMzw+59u3lg7QNUV1bj7rz8ystUVlRy+ITDOWHqCZ1Prura7JlSFaRMITXz+MNQoEFAJFlmNqjHa6ONpSxlEYsAqK/PffBDCqEAP8zEa+bxNAW5NfOexm7V06sy2PpTszYzLrn7Eu57/j5qx9RiZux8ZSc1VTW84bA3sHbrWm4+9+YESjuyKMAPM62tsG8fPP54V8qB6dMPrJk3NmoQEClvmV0ZRleOZufenRjGrr27mDN1To8DhUjx9KBTmWtr63o6taUF7r0XnngCdu6E2trw7xNPhHTAcYU+DCVSKgtPWciJU0/kpZ0vAfD6htczunK0csMPItXgy1i8OQZCgH7qqfBE6sEHh2U1NbB7d8j1HtfUBOecA9/4RsgeOX06fPjDan+XwVVsXve4pilNfOtt3+p2jKk1U9VTZhApwJexfDdKx47tGrijuhp27Qo19AkTuu/b1gZ33gnHHw9veENonrnzzpC2QEFeBkO8q2NDbUNnOt+eRnKKB/Ls63he+Oy6bLZIBfmBU4AvY/lulDY2huU1NfDii6H2XlERhur7wAdgz56wzfr1A+9Fo3700pt86Xyzy3ODc+6XAdD5ZZCdL/SLQgqnNvgy1tgYat5x06eHtvdp02DMmFBzr6yEzZvhF78IN19/9CP4wQ9C00xcMb1oCh1UREau9o526qrrui2rq67Lm843/mVQYSHsZHO7565T3vfBoxp8GWtuDkE1a9Girr7GPQXa1au75n/84zDluuWW/PvG+xqrH730pbGukczOTGfNHeixB0xunhno/mXQ2zrpPwX4MhbPGnnRRd6tmWT+/FCzrqiApUtDDfuVV0IXytmz4d57jdGjndpaOPbYsG2+HPA9UT966Uvz7GYWPxBqIHXVdXTs6uC5zHNM3zed+Uvnd7vpGv8yWLd9HQB3PHUHkw+azDETj+mWgwZ6/qKQ4ijAl7meskbG+7nX1YU2+KqqEPCfeSZskx3gY+PGcGP2iisKr32rH730pWlKEwtPWdh5c3RM5RgcZ0zlGCYfNLlbW3r2y2Djjo2s3LgSgFE2iunjp7Nm6xoM4wiO6PyiUFKxwZFoG7yZnWVmfzCzP5rZp5M810gT7+d+9NHgDtu3w44dIaDHnX46zJsXRnXqz/HVj1560jSliZZ5Ldx87s28atyrmFk/M29bevbL4IXtL7B3/14ATm08laMOOYqZ9TOZXjud+pp61m5dS31NvW6wDpLEko2ZWSXwDPBXwFrgYeB8d3+qp32UbKw48V4u27bBL38ZHnzavx/27TNqapyqqjCC09lnh1r+zUU8/a1eNFKM+Uvn01Db0HkTddEZi/p1HCUVK06pko2dDPzR3Z+LCvED4FygxwAvxcltvjn66FCL7+gINfm6uq4ulP1pXtGgIlKM3JuuV95zZefrlnktQFcWyXh7e+42MniSbKKZDqyJvV4bLZNBZGad0zPPGC+8YOzYEXrbrFsXXq9fb1xzjbFoUde2IoOteXYzmV0ZMjsz7Pf9nQN0xNMOFLKNDJ4kA3y+KHLAby8zu9jMlpvZ8o25jcfSJ3fvnN7/fmf2bOfkk505c5xp05y6Ouekk5zHH/du24oMtmw7e29t6YVsI4MnySaatcChsdcNwAu5G7n7DcANENrgEyxP6n3oQ6E5ZuPG0F2yoQEmTYIvfEFNLTI04qkHBrKNDI4kA/zDwFFmdjjwZ+DdwHsSPN+I19QUgrlujIoIJBjg3X2vmX0Y+AVQCdzs7k8mdT4JdGNURLISfdDJ3e8G7k7yHCIikp+SjYmIpJQCvIhISinAi4iklAK8iEhKJZaLpj/MbCPwfKnLkRITgU2lLoRID3R9Dp7D3H1SvhVlFeBl8JjZ8p4SEImUmq7PoaEmGhGRlFKAFxFJKQX49Lqh1AUQ6YWuzyGgNngRkZRSDV5EJKUU4EVEUkoBfpgxs+15ln3czJ4yszYz+7WZHZaz/udmNt3MbjKzx6Pt7jCzcUNXchmp8l2zsXXnmZmbmbpMJkABPh0eBea6exNwB/Dv2RVmVgMc7O5/Bj7m7sdH27UDHy5JaUUAMxsPXAI8WOqypJUCfAq4+z3u/nL08veE0bOy5gHLou22AlgYlLWGPEMoigyhzxEqI7tKXZC0UoBPnwXAz2Kv3wL8PPvCzG4B1gGzga8PbdFEAjM7ETjU3e8qdVnSTAE+RczsfcBc4EuxxacC92dfuPtFwDRgFfB3Q1pAEcDMKoCvApeVuixppwCfEmZ2JnA5cI67746WHQGscfc98W3dfR/wQ+AdQ15QERgPHAcsM7PVwOuBO3WjdfApwKdA9HP3ekJw3xBb1dk8Y8GR2Xngb4Cnh7qsIu7e4e4T3X2Gu88g3Dc6x92Xl7hoqZPomKySiLFmtjb2+ivA2cA44D9C7Kbd3c8BzgI+Em1nwHfMrDaafxz44JCVWkayA65Zd/9KyUozgihVQUqZ2Rjgt0rJKjJyKcCLiKSU2uBFRFJKAV5EJKUU4MuYmdWY2b1mVjmAY7w3yj3TZma/M7PjY+smRDlpnjazVWZ2SrR8sZm9cTDeg6TTEFybq83sCTN7zMyWx5br2iyCetGUt/lAa9Rvvb/+BJzu7hkzewthoIXXReu+Bvzc3c8zs9HA2Gj514FvA78ZwHkl3ZK+NgHOcPfcgbl1bRZBN1nLmJn9DngPMANoIYxCfxywAnifF/nHM7N6YKW7T4+6Sz4OHJHvOGa2Aniru68b0JuQVEry2oxeryYk0MsN8Lo2i6AmmjIV1aiPcPfV0aITgUuBY4EjCCkIMLOvRj9jc6dP5zlsPE/NEcBG4BYze9TMbjSzg2LbPpI9h0jcEFybEBLh/dLMVpjZxTnb6toskJpoytdEYEvs9UPuvhbAzB4j1Jzud/ePFXIwMzuD8J/otGjRKGAO8BF3f9DMvgZ8GrgiWr+BkLNGJFfS1ybAqe7+gplNBn5lZk+7+33ROl2bBVKAL187gerY692x+X1Efzsz+ypwRp79f+Du/xZt0wTcCLzF3TdH69cCa909m4v7DkKAz6qOyiCSK+lrE3d/Ifp3g5n9BDgZyAZ4XZsFUoAvU9GNp0ozq+5ju15rSWbWCLQCF7j7M7H91pnZGjM72t3/ALwJeCq26yzgP/r/DiStkr42o6bCCnffFs3/NfDZ2K66NgukAF/efkn42bp3AMf4V+AQ4JtRnpq9sfQFHwG+F7WpPgdcBGBmVcCRgJI/SU+SvDZfBfwkWjYK+L67Z5Pm6dosgnrRlLEoS+TH3f2CIT7v3wJz3P2KPjeWEUnX5vCgXjRlzN0fBe4ZyMMk/TQK+PIQn1OGEV2bw4Nq8CIiKaUavIhISinAi4iklAK8lIyZXW5mT0bJph4zs9dFyy81s7Gx7e42swmDdM7tfayfYGYfir2eZmZ3DNK5l2ncURlKCvBSElHmyrcRekQ0AWcCa6LVl9KV+Ax3P9vdtwxR0SYAnQHe3V9w9/OG6NyDogQ3PqVMKcBLqUwFNrn7bgB33xQ9mn4J4TH0e8zsHuhMHTvRzGZEqY1vNLOVZvY9MzvTzH5rZs+a2cnR9i1mtjB7omjbGfGTm9k4M/u1mT0SpaU9N1r1b8DM6BfFl6Jzroz2qTazW6LtH40escfMLjSzVjP7eVSOfy/0Q4iO/z9ROR4xs7+Ilt8eKxPRez0nesDoS2b2cPTL5/3R+nlmdo+ZfR94oqi/hKSXu2vSNOQTYZDwx4BngG8S0sZm160GJua+JuQ42Qu8hlA5WQHcTBhE/Fzgp9H2LcDC2P4rgRnR/Pbo31FAbTQ/EfhjdJwZhKyG2X07XwOXAbdE87OBdsJj8xcSHhSri14/Dxya5z0vI2RIjC8bC1RH80cBy6P502Pvp46QWncUcDHwmWj5GMIDP4cD84AdwOGl/ttqKp9JNXgpCXffDpxECFgbgR+a2YUF7Pond3/C3fcDTwK/dncn1FpnFFEEA75gZm3AfwPTCU9Q9uY04Pao/E8TAvmsaN2v3b3D3XcRUj4cVmA5qoBvm9kThMfvj42Ofy9wZJRs63zgx+6+l/DY/t9HSb0eJDwJelR0rIfc/U8FnldGAKUqkJLxMFjEMmBZFOD+Abi1j93iia32x17vp+t63kv35sd8OVPeC0wCTnL3VyzkH+81twrhS6GQcnUm3CrAx4D1wPGEMu+Krbs9Kue7CQNsZMvwEXf/RbeCmc0j1OBFOqkGLyVhZkeb2VGxRScQasQA24DxAzj8akIqZMxsDqEJI1cdsCEK7mfQVePu7dz3EQIuZjYLaAT+MIByZsvxYvSL5AIgfoP0VsINZ9z9yWjZL4APRjlZMLNZ1j2Pv0gn1eClVMYBX4+6P+4ltIFnB3a4AfiZmb3o7vnSzfblx3Q1YzxMaOfP9T3gPy2M9/kY8DSAu2+ObtquJAxAcW1sn28C34p+bewFLnT33VFSrEL9l5m9Es0/APwL8GMzeydwD7FauLuvN7NVwE9j+99IaIp6xMKJNwJvL6YAMnIoVYFImYqeBXiC0JW0o9TlkeFHTTQiZcjMziT8qvi6grv0l2rwIiIppRq8iEhKKcCLiKSUAryISEopwIuIpJQCvIhISinAi4ik1P8BfQJG91TcNHAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now do ALL latencies and see if there are any multimodal distributions\n",
    "plt.clf()\n",
    "i = 1\n",
    "xtick_labels = []\n",
    "df_datasets = []\n",
    "df = pd.read_csv(data_dir + csv_prefix)\n",
    "df = pd.DataFrame(df)\n",
    "df = df[df[\"Exclude Reason\"].isnull()]\n",
    "for stim_loc in ['L2/3', 'L4']:\n",
    "    df2 = df[df['Stim Location'].str.contains(stim_loc)]\n",
    "    \n",
    "    df_d1 = df2[df2['abs(delta 1)'] > 0]['abs(delta 1)']\n",
    "    df_d2 = df2[df2['abs(delta 2)'] > 0]['abs(delta 2)']\n",
    "    df_lat = pd.concat([\n",
    "        df_d1, df_d2\n",
    "    ])\n",
    "    jittered_xs_1 = np.random.normal(i, 0.04, df_lat.shape[0])\n",
    "    plt.scatter(jittered_xs_1, df_lat, alpha=0.5, color={\"L2/3\" : 'blue', \"L4\": \"green\"}[stim_loc])\n",
    "\n",
    "    xtick_labels.append(stim_loc + \"\\n(n=\" + str(len(df_lat)) +\")\")\n",
    "    df_datasets.append(df_lat)\n",
    "    i += 1\n",
    "    print(np.mean(df_lat), \"+/-\", np.std(df_lat))\n",
    "\n",
    "p_val = kstest(df_datasets[0], df_datasets[1]).pvalue\n",
    "print(\"KS test p-value \", p_val)\n",
    "# significance: L2/3 Rostral and L4 Rostral, and \n",
    "#               L2/3 Rostral and L4 Caudal\n",
    "#significance_bar(0.75, 2.25, 3.35, \"P<0.05\", fontsize=10)\n",
    "plt.boxplot(df_datasets)\n",
    "plt.xticks(range(1, len(xtick_labels)+1), xtick_labels)\n",
    "plt.title(\"Latency of L4 Barrel Pair Crossings\")\n",
    "plt.xlabel(\"Stimulation Layer\")\n",
    "plt.ylabel(\"Latency (ms)\")\n",
    "plt.savefig(data_dir + cell_type_prefix + \"_results_latency_medial-bias_\" + str(datetime.date.today()) +\".png\", bbox_inches=\"tight\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "a5f8e9e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEuCAYAAACESglMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABHGElEQVR4nO2de3xU5Zn4v0/u4ZIw3JUYEVBBaQBlW7BUobqt9qI2pe2yXWsLu3bbXXtFd7esv8Z1sd2WtlbdttuKWlulF5oq2tpaFbB4R4QIJfWKIVwDDOGWBJI8vz/ec5KTYUImyUzm9nw/n/OZOZc555l3znnmfZ/3uYiqYhiGYWQeOckWwDAMw0gMpuANwzAyFFPwhmEYGYopeMMwjAzFFLxhGEaGYgreMAwjQzEFbyQcEakSkZ8nW45kIyLFIvKwiDSKyK8TfK33iMhfE3Tur4nIXYk4txFfTMFnASKyRkTCIlIY4/GfFpF1iZbLu9ZcEanvw+fGi4iKSF4i5EoQ84ExwAhV/Vh3B3ntryLy8VhP7B0/yV9X1T+r6rn9Ezf676Oqt6rqP/b33EbiMQWf4YjIeOA9gAJXJlearOdM4FVVbe3huGuBA96rYfQZU/CZz6eA54B7iVAYInKGiFSLSIOI7BeRO0VkCvAjYLaIHBGRg96xa0TkHwOf7dLLF5Hvi8h2ETkkIi+JyHv6K7iIfFBEXvbOuV1EqgK7n/JeD3pyzvY+s1BEtnojlj+KyJmB86mI/LOIvObt/18RkcD+f/I+e1hE/iIiF4jIDSLymwi57hCR27qReYrXVgdFZIuIXOltvxn4f8AnPHkXdfP5M4FLgOuA94vImMC+XM888oYn40veb+i3xSbv3J8I9rxF5N9FZGXEdb4vIrd77z8T+N5vishnve2DgUeB073zHhGR0yNNbiJypfddD3rffUpg3zYRWSwiNZ5p6pciUhTtuxsJQFVtyeAFeB34PHAhcAIY423PBTYB3wMGA0XAHG/fp4F1EedZA/xjYL3LMcA/ACOAPOCrwG6gyNtXBfy8G/nmAvWn2PcOXEekAtgDXO3tG48bleQFjr/a+75TPDn+E3gmsF+BR4BhQDnQAFzu7fsYsAP4G0CASbge92nAUWCYd1wesBe4MIq8+d71vwYUAO8FDgPn9tQOgXPcBLzgvX8F+Epg3w3etnM9GafhzD3+d5sUrV2973EMKAn89ruAWd76B4GJ3jkv8Y69oLvfJ/g9gHO89vlb7/vf6LVBgbd/G/ACcDowHNgK/HOyn4tsWawHn8GIyBzcw/0rVX0JeAP4e2/3O3EP3Q2qelRVm1W1z3Z3Vf25qu5X1VZV/Q5QiFNEfUZV16jqK6rarqo1wAqcAuqOzwLfUNWt6swgtwLTg7144JuqelBV64DVwHRv+z8C31LVF9Xxuqq+raq7cKMF32Z+ObDPa89IZgFDvGscV9UncX8oC3rxtT8FPOC9f4Cuo65/BP5TVf/qybhJVff3dEJVfRvYgPsDBPfHc0xVn/P2/05V3/DOuRZ4DGfWi4VPAL9T1T+p6glgGVAMXBQ45nZV3amqB4CH6WxzI8GYgs9srgUeU9V93npQYZwBvK0924NjQkS+6g3zGz2zTikwsp/nfJeIrPZMSI3AP/dwzjOB73umgoM4O7YA4wLH7A68P4ZTyODa441uzvtT3AgF7/Vn3Rx3OrBdVdsD296OuH63iMi7gbOAX3ibHgDeISLTY5CxJx6g84/m7+n8E0FErhCR50TkgNduHyD23+503HcEwPvu24mtzY0EYwo+QxGRYuDjwCUisltEdgNfBqaJyDTcQ1gu0b1QoqUYPQoMCqyPDVzrPcC/edcLqeowoBGnXPvDA8Aq4AxVLcXNDfjnjCbjduCzqjossBSr6jMxXGs7zkwRjQeBChGZCnwIuL+b43YCZ4hI8Lkqx5l+YuFa3Pfb6P1ez3vbPxWDjD3xa2CuiJQBH8FT8OI8q36D63mP8X6733Pqdg6yE/fHinc+wf0RxfqdjQRiCj5zuRpoA87DDYmn42zTf8YpjBdwdthvishgESnyepDgbN1lIlIQON9GoFJEBolzxwtOEg4FWnE27TwR+X9ASW+E9a4fXMQ77wFVbRaRd9JpXsK7VjswIbDtR8B/iMj53jlLRaRbd8QI7gIWi8iF4pjkm3ZUtRlYiVOKL3jmnWg8j/sjvFFE8kVkLvBhOnvkp/z+uD/I6+j8vaYD1wOf9P6I7wJuEZGzPRkrRGSEd4o9dG2LLqhqA24e5R7gLVXd6u0qwJnTGoBWEbkCeF/go3uAESJS2s2pfwV8UEQuFZF83PxLCxDLn6qRYEzBZy7XAveoap2q7vYX4E7gk7ge2odxk4l1QD3OngrwJLAF2C0ivnnne8Bx3AP/U7r2Yv+I87Z4FTdcb8b1NmNlHNAUsUzETQ7/l4gcxnmg/Mr/gKoeA5YCT3smmVmq+lvgf4BfiMghYDNwRSwCqOqvvfM9gJsYfRA3KejzU9yEb3fmGVT1OM4V9QpgH/AD4FOqWhuDCFfjvvd9Eb/Xctyk6OXAd3Ft8BhwyNtX7H2+Cvip1xbd+c8/AFxGwDyjqoeBL3jnDeP+RFcF9tfi5j7e9M59esR3/ivObHWH950/DHzYawsjyYiqFfwwjJ4QkXKgFhirqoeSLY9hxIL14A2jBzyb+leAX5hyN9KJdArzNowBxwv22YMzPV2eZHEMo1eYicYwDCNDMRONYRhGhmIK3jAMI0NJKRv8yJEjdfz48ckWwzAMI2146aWX9qnqqGj7UkrBjx8/nvXr1ydbDMMwjLRBRN7ubp+ZaAzDMDIUU/CGYRgZiil4wzCMDCWlbPDROHHiBPX19TQ3NydblIRQVFREWVkZ+fn5yRbFMIwMI+UVfH19PUOHDmX8+PGI9Df7bGqhquzfv5/6+nrOOuusZItjGEaGkfImmubmZkaMGJFxyh1ARBgxYkTGjk5OxYoVK5g6dSq5ublMmjyJj970URY+tJCqNVXU7K5JtniGkRGkvIIHMlK5+2Tyd+uOFStWsGTJEu644w5e2PYCE/5+Ao/f9TgHXzxIuCnMsmeXmZI3jDiQFgo+2QwZcnKFse9+97ucd955VFRUcOmll/L2211dUS+//HJ27NjBokWLmDZtGhUVFcyfP58jR44MlNgpy9KlS1m+fDnz5s3j4TceZsrfTOHqG69m3f3rCBWHCBWFqK6tTraYaUVwRDR69FTmzl1BVRXU2P9kVmMKvo/MmDGD9evXU1NTw/z587nxxhs79jU1NXHgwAHGjRvH9773PTZt2kRNTQ3l5eXceeedSZQ6Ndi6dStz5swBoK6xjtKiUsrfUc6+OldbpLSolLrG7oomGZH4I6IvfekOFixoZs6cO9i0aQnPPbeCZctMyWczpuD7yLx58xg0yJUonTVrFvX19R371qxZw9y5cwEoKXGV61SVpqamrDTJRDJlyhTWrVsHQHlpOY3NjdS9UsfIclfnubG5kfLS8mSKmFb4I6L6+nmMGJFPRcU8rrxyORs2LCUUgmobDGUtGafga2qgqgoWLmTAhqjLly/niis6K8M9+uijXH55Z+rwz3zmM4wdO5ba2lquv/76xAuU4ixZsoRFixaxevVqPjzxw2x9cSsPfutB5nxyDuGmMOHmMJWTK5MtZtrgj4jq6qDUq5xaXj6Hffu2UloKdTYYyloySsHX1MCyZRAOQ1mZe030EPXnP/8569ev54YbbujY9vTTT3eYIADuuecedu7cyZQpU/jlL3+ZOGHShAULFrB06VKuv/563jn+nbz5wJtc9o+XMexvhhEqDrF49mIqxlYkW8y0wR8RlZdDY6PbVle3jpEjp9DYCOU2GMpaUt4PvjdUV0Mo5BbofK2uhooE6IvHH3+cpUuXsnbtWgoLCwF48803OeOMMygoKOhybG5uLp/4xCf49re/zWc+85n4C5NmLFiwgAULFiRbjIzAHxF97WvLef31OWzfvo7VqxfxrnctJRyGRYuSLaGRLBKq4EVkGHAXMBVQYKGqPpuo69XVuZ57kEQNUV9++WU++9nP8oc//IHRo0d3bA+aZ1SVN954g0mTJqGqPPzww0yePDn+whhZjf9HuXTp9WzdupURI6YwbdpSZs1aQGVlYjo3RnqQ6B7894E/qOp8ESkABiXyYuXlzizj99yBuAxRjx07Rlngn+MrX/kKv//97zly5Agf+9jHvGuXs2rVKv7whz9wxx13AE7BX3vttRw6dAhVZdq0afzwhz/snzCGEQUbERnRSJiCF5ES4GLg0wCqehw4nqjrAVRWOps7uJ57YyNxGaK2t7eftO0rX/nKSdtaWlrYtWsXftGSnJwcnn766f5d3DAMo48kcpJ1AtAA3CMiL4vIXV6F+oRRUQGLF7sefH29e128eOCGqIWFhVawxDCMlCGRJpo84ALgelV9XkS+D/w7cFPwIBG5DrgOnJmjv1RUmM3RMAwDEtuDrwfqVfV5b30lTuF3QVV/rKozVXXmqFFRywoahmEYfSBhCl5VdwPbReRcb9OlwF8SdT3DMAyjK4kOdLoeuF9EaoDpwK0Jvp5hGEafqdldQ9WaKuZ+dS6jzxpNbm4uU6dOZcWKFckWrU8k1E1SVTcCMxN5DcMwjHhQs7uGZc8uY+9ze9m0YhPzvjCPQRMHMTdnLktuWAKQdq6oGZWqIFFESxfss3LlSkTEvGcMI82prq0mVBRiw682cOUNV1Ixq4IRQ0ZQP6Ke5cuXs3Tp0mSL2GtMwfeDw4cPc/vtt/Oud70r2aIYhtFP/NTV++r2Uf4O59Hnp66eM2cOW7duTbKEvccUfD+46aabuPHGGykqKkq2KIZh9BM/dfXI8pHUveLym/ipq9etW8eUKVOSLGHvyTgF70+SJLq+58svv8z27dv50Ic+lJDzG4YxsFROriTcHOaCj1/Aqm+voua5GvYf2U/Z/jIWLVrEkiVLki1ir8mobJL+JEmoKERZSVlHfc94p59tb2/ny1/+Mvfee2/czmkYRnKpGFvB4tmLqQ5V09zazLr/W8f+uv1snLKRpUuXpt0EK2SYgvcnSULFLtuY/1pdWx1XBX/48GE2b97cUbVp9+7dXHnllaxatYqZM81pyDDSlYqxFU5XzAW+k2xp+k9GKfi6xjrKSrrmC05Efc/S0lL27dvXsT537lyWLVtmyt0wjJQio2zw/iRJkHjU9/TTBfvLd7/73X6dzzAMYyDIqB585eRKlj3r8gWXFpXS2NxIuDnMohn9yxccLV1wkDVr1vTr/IZhGIkgo3rw/iRJqDhE/aF6q+9pGEZWk1E9eAhMkhiGYWQ5GdWDNwzDMDpJCwWvqskWIWFk8nczDCO5pLyCLyoqYv/+/RmpCFWV/fv3W6oDwzASQsrb4MvKyqivr6ehoSHZoiSEoqIiysrKej7QMAyjl6S8gs/Pz+ess85KthiGYRhpR8qbaAzDMIy+YQreMAwjQzEFbxiGkaGYgjcMw8hQTMEbhmFkKKbgDcMwMhRT8IZhGBmKKXjDMIwMxRS8YRhGhmIK3jAMI0MxBW8YhpGhmII3DMPIUEzBG4ZhZCim4A3DMDKUhKYLFpFtwGGgDWhV1ZmJvJ5hGIbRyUDkg5+nqvsG4DqGYRhxo2Z3DdW11dQ11lFeWk7l5EoqxlYkW6xeYSYawzCMCGp217Ds2WWEm8KUlZQRbgqz7Nll1OyuSbZovSLRCl6Bx0TkJRG5LtoBInKdiKwXkfWZWpbPMIz0orq2mlBRiFBxiBzJIVQcIlQUorq2Otmi9QpJZDFrETldVXeKyGjgT8D1qvpUd8fPnDlT169fnzB5DMMwAEQkrudLpB7tCRF5qbv5zYTa4FV1p/e6V0R+C7wT6FbBG4ZhDAQ9KeSqNVWEm8Lc/oHb+frqrwMQbgoTKg5RNbdqACSMDwkz0YjIYBEZ6r8H3gdsTtT1DMMw4kXl5ErCzWEA2rWdcFOYcHOYysmVSZasdyTSBj8GWCcim4AXgN+p6h8SeD3DMIy4UDG2gsWzFwNQf6ieUHGIxbMXp50XTUJt8L3FbPCGYaQSIpJU+3osnMoGb26ShmEYGYopeMMwjAzFFLxhGEaGMhCpCgzDiDM1NVBdDXV1UF4OlZVQkV7zf8YAYD14w0gzampg2TIIh6GszL0uW+a2G0YQU/CGkWZUV0Mo5JacnM731ekVRW8MAGaiMYw0INbQ+ptvju18qe76Z8QHU/CGkQYEFXJVlTPLhEKd+/31m29Ofb9tY+AwE40xYNTsrqFqTRULH1pI1ZqqtEu9mipUVjqFHg5De3vn+8r0iqI3BgBT8EbCqdldw+ce+RyVv6rk0dceJT8nP23za6cCFRWweLHrsdfXu9fFi82LxjgZM9EYCcUvnPDXfX9leNFwAJ7b8Ryzy2Z35Nf283tkQgWdgaKiwhT6QJHO96XlojESRibl3E4X0iF3SjohIlxTfQ2hohClRaU0NjcSbg6nVOKxpOWDN7IbVWXhQwspKynjqbefoulEE8X5xagqh1oOcfGZF3P7B25HVTvyb4eKO2cO0zH/tpF5+JWdgI7X4MgzlTEbvJFQykvLaWxuZMrIKTS3NtN0oomm1iYKcgs68m0D1DXWUVpU2uWzpUWl1DXWDbTIhtGFdL4ve1TwIpIjIjNE5IMi8l4RGTMQghmZgV84oSC3gFllswA40HSAGafN6Mi3DZ1/BEEamxspLy0fUHkNI5J0vi+7NdGIyETg34DLgNeABqAIOEdEjgH/B/xUVdsHQlAjPfELJ/iTVFecfUXUSarKyZUse3YZQBdb56IZi5IhtmF0eHj97rXfMaJ4BNPHTqcoryit7stuJ1lFZAXwQ+DPGnGQV0T774Gwqv40XsLYJGv2EZwUTGdvhVTBJln7h38Pbty1kbca36LmczVct+o6Nu7eyP6m/fzthL/lczM/l1L3ZZ8mWVV1wSn27QVu679ohtFJxdiKlHpwjOzCd+kNFYUIt4QRdV5gOZLDFWdf0THpn073aCw2+I8Fimf/p4hUi8gFiRfNMAxj4Kiure7wmDnUfKhjcnXrvq1Aek2u+sTiJnmTqv5aROYA7weW4Uw370qoZEbGEc0EY/SPyLzwRt+pa6yjrKQMcMq86UQT0DnJmk6Tqz6xuEm2ea8fBH6oqg8BBYkTychE/OFvuClMWUkZ4aYwX3vyawCWm6aPRMsL7283ek/Qk2vyyMk0tzUDUFJYQrgpTLg5nHadklgU/A4R+T/g48DvRaQwxs8ZRgfB4W+O5NDS1sIbB94A6FD4lpumd0TLC+9vN3qP79IbbgozevBozh95PtAZ6JRK0auxEouJ5uPA5cAyVT0oIqcBNyRWLCPTCA5/AWr31TK0cCjgJrHSLUIwFaircz33aNuN3hPp0nvOyHMAeHDBg8kVrB/0qOBV9ZiIrAbOCEyu7kusWEamUV5a3iUVQWNzI/k5+V2OScdJrGRSXn5yXnh/u9E3Ij25bibGCiopSixeNLcANcDtwHe8ZVmC5TIyjODwt13bKcgt4PDxw12OScdJrGQSLS+8v90wIDZb+seBiao6V1Xnect7Ey2YkVn4w99QcYj6Q/XMOG0Gk4ZPAqBd29N2EiuZRMsL7283DIjNBr8ZGAbsTawoRqYTOfyt2V3DIzzCw399GEU7ctUYsROZFz7WmqxGdhBLD/4bwMsi8kcRWeUviRbMyB4uPvNirjz3SgpzC82TxkgJ/PKSQFq78MbSg/8p8D/AK4AlFjPiRnWt8+dL11zbRmYQGYA3ddRUVr26ilCRux99F95MdZPcp6q39/UCIpILrAd2qOqH+noeI/OI5jFjnjTGQBLMP+PHY9zy51s4f+T5GdHxiMVE85KIfENEZovIBf7Si2t8EdjaR/lShhUrVjB16lRyc3OZNGkqH/3oChYuhKoqixzsK9E8ZsyTxhhIIgPwQsUhTrSdYMeRHV2OS9eORywKfgYwC7iVXrpJikgZLsXBXX0VMBVYsWIFS5Ys4Y477uCFF5qZMOEOHn98CQcPriAcduHipuR7j+8x47tOmieNMdBEqyQ2atAoGo42dNmWrh2PHhV8wDVyXh/cJG8DbiTNbfdLly5l+fLlzJs3j4cfzmfKlHlcffVy1q1b2hEqbuHhvccf7vquk+kaDm6kL9EqiZWVlJGfk0+4yQUWpHPH41QVnf4BeKC7ik1exafTVHVdN/s/BOxV1ZdEZO4prnMdcB1AeYqG4G3dupU5c+YAneHhJSVz2OenES218PD+YEW1jWQRrZJYbk4uN118E5sbNgOuA7JoxqK07HicapJ1BM498iXgJTpL9k0CLsGlK/j3U3z+3cCVIvIB73MlIvJzVf2H4EGq+mPgx+AqOvX1iySSKVOmsG7dOubNm9cRHn7w4DpGjpwCQGOjhYcbqUNkCuHKSgt+6o7I/DPlpeUdynw+87mZm9O6A9JtyT7o8IB5L05ZnwY04SZMH1XVmPusXg9+cU9eNKlass+3wS9fvpySkjn8x3+s4/nnF3HppUs544wFhMMuotAeolMTLR/8tNOmWYm5OCIiXHONEgq5kWVjI3Z/9oN0KIHYp5J9AKraBvzJW7KWBQtc9cLrr7+erVu3ctZZU7jssqWUli4gFIJFi7Lr4Rk+fDhhP/FJP/ATOYlIv88FEAqFOHDgQFzOlc7480L+e3A9+my6Rw1HLH7w/UZV1wBrBuJaiWLBggUdij7bCYfDve7VVK2p6pJNEuhYj9cQOF5/FOlOaenJ6zZHlJ1Y4Q5jQIjmjpauvsWpTmPjyes2R5SdxJIuOHcgBDEym2juaOnqW5zqRKYQDoezK4Xw8OHDEZG4LEBczjN8+PCktEUsJprXRWQlcI+q/iXRAhmZSTR3tHBzmEUzFiVZssxj8eKuXjTZNkfUFxNiJMEUBsH7ta9xGskyH8ai4CuAvwPuEpEc4G7gF6p6KKGSGRnFqdzRjPgSmULY6D3BFAaQvvloYinZdxj4CfATEbkYWAF8z+vV36KqrydYxqQR9CdubFzBpk1LeeutrUyZMoUlS5bYpGsvicwHHySaC2U6PUhG5lCzu4aHah9CVRlWPIwpI6cwZsiYtJwz6lHBezb4DwKfAcbjctHcD7wH+D1wTgLlSxorV8Itt8CJE6C6grfe+iJ5eYNRVY4ePcoXv/hFAFPy/cBX6ht3beStg29x/qjzmTRiUlqnZzXSG980U5BbgKI0nWjime3PcNEZF1GQW5B2c0axeNG8BlwFfFtVZ6jqd1V1j6quBP6QWPGSQ02NU+4iMGoUvPHGjbS15TJr1t0sWdLC3XffTW5uLjfeeGOyRU1b/AfJz/PRdKKJJ7Y9wS+3/JJNezbR2tbakS/eMAYK3zQz47QZtLS1AFCYW8iGXRvSMh9NTDZ4VT0SbYeqfiHO8qQE1dWu5z5qlFPyJ07Uc+65j3HixDx27IB58+Zx33338b73vS/ZoqYtQRvn7iO7aWxxHjbNJ5ppOtHE5obNHDtxLMlSGtlGXWMdZSVl5EgOs8tmU7uvloPNBwHSckQZi4L/XxH5oqoeBBCREPAdVV2YUMmSSF2dU+7NzVBc7Lbl50NDA8ybl1zZMgX/QQI6ekpFeUW0tLXQpm3sObKHvUf3UrWmyuzxRsLxzYUbdm1gy94tXHDaBYwdMpaxQ8Z2BOSl4z0Yi4mmwlfuAKoaxuWIz1jKy2HcOKfgm5qgoKCM2tpraWlZzYc/fILVq1dz7bXXUlZWlmxR05agX3xRXhGqSktrC4Kw7eA2WttbKS0s7bDHp2tNTCP1WbllJdf89hp+teVX5Eoue4/uZc22New6vCutUwVDbAo+x+u1AyAiwxmgFAfJorIS8vLg/POhqAhGjvwWqq0cPryQd76ziIULF9La2sq3vvWtZIuadvjFjDfu2siabWt4dd+rjBkyhlBxCBXlRNsJCnMLOX3o6Zw29DRCxSFCRSGzxxsJoWZ3Dbf8+RYEYdSgUeTl5JGfm09eTh4v7HiBUHGIK8+5kuraahY+tDDtCnDHouC/AzwjIreIyC3AM0BGa7aKChcscs45MHEi/NM/LWDp0u8zevRgAAYPHsz3v/9986DpJcGJ1YqxFUwdPZUt+7aQozkoyuhBo2nTNsCZbSaPnAxYSgMjcVTXVnOi7QSlRaWICMX5xYSKQgwfNJwZp82gcnIlq15dRbgp3FGzNZ1GlLH4wd/n5YSfBwhQmQ0RrScHiyzg3/7NFHp/iAweOXvE2YwcNJKWthYGFw6m4WgDhXmFtLW3gYLgov8spYGRKOoa6xg1aBTNrc0U57sJt6K8IhqONjBv/Ly0D3iKNdlYLVANPAQcERF72oxe013Csefqn2NiaCJXnH0FV0++mjFDxzAofxB/afhL2ttAjdSmvLScspIymlud95aq0tjcSH5OPpWTK9M+SV4sycauB/bgcsI/AvzOezWMXtFdwjFBOh6isUPGMrtsNqXFpew8stPqtPZATQ1UVcHChe7V6B2VkyvJzcll6uiprud+rAEV5aaLb6JibEXaJ8k7ZUUnABF5HXiXqu5PtDCpWtHJ6Epfq9x0l8BpcP5gCnIL+p0rPh2q78STmhpYtowu1Ztuv13YtEmzOhdNb++DyDQZU0dNZXPDZuoa6yjILWDHoR1MCE3oV9KxRN6bfa7o5LEdaOzxqAzCalomhu4SjgGWabIPVFdb9aZ4EMyRFOyElJWU0djciKIcbztO/aH6tEuSF4uCfxNYIyK/A1r8jar63YRJlUSCvaKyMpdLe9kyq2kZRL9eAlWlPR8YhQpv6WDN9wC4L9rBmx7qvVxZRF2du0ejbTf6RrRJ1YlMjGvlsYEkFgVf5y0F3pLRWK+oZ+TmQ3EZbiYi57ZW9VustKG83HVAQiHYvRtqa932N990HRW7X3tPMMLaJ50mVSOJxU3yZgARGayqRxMvUvLoKSn/zTf37nzZZA/uC+nugpZsKivd6LKhATZvhhzPZWLcOLf9yivddjM1xk55aflJtYPTaVI1kljSBc8GlgNDgHIRmQZ8VlU/n2jhBhpVpaqqs1fk469XVWXfRF4iybTe0kDjB+R94QvQ2gqjR7vtZ58Nr73mMqLOnWumxt6QaZXHYvGDvw14P7AfQFU3ARcnUKakUllpNS0HinR3QUsFKipgwgSYP98pc5/6epcRNRRyPXvf7FhtGR9Oie8IECoOUX+oPu3ddGPKKaOq2yPMF22JESf5+L2ibK5pGW+6q9ZUObmSrz35NRrqG2hpa6Ewt5BRg0dx64xbky1yWhG0xfs0NLiMqEFKS20CNhYivWp+uP6HPFf/HIoyq2wWn5/5+bRR+DG5SYrIRYCKSAHwBWBrYsVKLlbTMn5Eup0FqzVBZzoCPKtXx7oREzU1boL18cdhxAi3LRx26a3Hjet6bGOj+zMwTiZaJwRgyZNLeP3A6wwtGAoCa7etpf5QPbe+99a0UPKxKPh/Br4PjAPqgceAjLO/G4nhVBOpABNCE7jw9As7jg83hW2SNUaCLr2XXgobN7rtx4/DTTfBqlVO2ftBUOGwG40aXemuEzI4fzB7j+6lpLCkI0+NIDQcbUibezQWBX+uqn4yuEFE3g08nRiRjEyip4lUm2TtO5EuvaedBi+8AGPGOJv8Oedkp6mxt3EaFXQTh3EqdtR2xHDELFMSiEXB3wFcEMM2wziJntzOMsklbaDpKdApW02NscRpBHvta7etpSCvgJa2FmaXzWbskLG0azu/2PwLjh4/yvG24wwuGMzI4pHk5uQCcMXZV/Q+jUbsh8eNbhW85x55ETBKRL4S2FUC5CZaMCMz6MntLJNc0gaaaJOrAAUFzqXX/N+7J2g6HFY8jKYTTRTlFlG7r5axQ8byYv2L7G/aT0trC23tbbS1t3Hk+BGG5g/l/DHnp012026TjYnIJcBcnA3+R4Fdh4GHVfW1eAuTDsnGzA++/8mcgjVWT7Uv0XKlCsOHDyccDidbjJMIhUIcOHAg2WL0mljug4UPLeworr3nyB6e2f4MhbmFHG8/znkjz+PR1x9l9ODRDCkYwvZD2znccpjC3EImjZjE3Vfe3et7NFnJxmLJJnmmqr6dEMkiMAWfHqRqG6SqXD3RH7n9xHgbN8LBg86jprQUZsyAsWPdMcFAvYGSK5nEInfVmqou5sE9R/awYdcGWtpaGFY0jO2N2ykrKeuIbm860URRXhETh0/k7qvuTohMfeVUCj6WQKdjIvJtEfm9iDzpLzFctEhEXhCRTSKyRUR6GehvGEZPVFQ4E0xJCUyb5moIq8KzzzplD+b/Ho3KyZWEm8OEm8K0azsFuQWcO/Jc7rv6PiaEJlBW6oqA+Pi54tNtfigWBX8/rqLTWcDNwDbgxRg+1wK8V1WnAdOBy0VkVt/ENAyjO4LeNMOGgYhT9H7yMfN/P5lTRayWl5Yzbsg4mtsiqjzl5qeN7d0nFi+aEaq6XES+qKprgbUisranD6kbjxzxVvO9Jf3Ge3TND++v26SVkSoEvWmmTIFnnoHCQmey8VNtmP/7yURGrFbXVnPb87dRkFvA4eOHOX/k+ew4soOGow3k5+Rz03tuSgvf9yCxKPgT3usuEfkgsBOI4px1MiKSC7wETAL+V1Wf75OUSSQyPzxY0qb+EM9JVcMR9KYZMwYuugg2bHD7QqHs8X/vK90V+SgtKqUov4h54+el7X0ai4L/bxEpBb6K838vAb4Uy8lVtQ2YLiLDgN+KyFRV3Rw8RkSuA64DKE/BcWRkMAl0Jm2yh6Z3nCptQTo+PKmCnzYYnL29oADOPdc6IbGSaUU+gvRog1fVR1S1UVU3q+o8Vb0QmNibi6jqQWANcHmUfT9W1ZmqOnNUZHakFKCuzj00QWzSqm8EH6QcySFUHCJUFOpIW2D0DT9BXijkRpybNsGhQ64TUlOTbOlSn7rGuo6i7z6ZElEdyyRrNL7S0wEiMsrruSMixcBluMnatKK83E1S7dkDa9a4bX/8o7NxGr0j+CDtObKHNdvWsHbbWh6qfYia3aaJ+kOkN01FRWcOeFPypyaT01b3VcHHkvLvNGC1iNTgvG7+pKqP9PF6SaOy0pVAW7MGjh1z2w4dgu3b7cHpLf6D5AeWNJ1ooiCvgILcApY9u8yUfD8JmhMtB3x0anbXULWmioUPLaRqTRU1u2tOcpkMN4UJN4fTzmMmGn1V8D16w6hqjarOUNUKVZ2qqv/Vx2sllYoKl3a1pMQVUAC45BKYONEenN7iP0gbdm2gMNcNgVraWphx2gwz1fSCmhoXtHT11a7Ix9VXu/WNG82ceCr8OaBwU7jLHBDQxWXyeNtxBuUP4rbnb+v4E0hXulXwInJYRA5FWQ4Dpw+gjEnn+HF4//vhqqvc+tix9uD0Bd/3uKWthePtxynOL+5I7pQpNs9E43t1vfaaG1kePOheX30V3noLXn+96/HmA9/JqeaA/AI0hbmF/OnNP/HyrpfJz8nv+BNIVyXfrReNqg4dSEFSmWhJnezB6RsVYyu4evLVlkWyj/hmmE2boLjYLU1NsHMnTJ3qimyPGmU54CMq0PXIzXQNtH+DN3iBFzrWf8bP+iVPKDIj3AARU8m+bCfohgYWPNJfMq2w8UDiBzU1NjqzIbio1cZGuPhiOHrU/QFkWw74IN3lfInMPwNd01WHm8I89fZTlBSWICLsP7afzf+ymWk/mIaI8NOrf5p27rx9tcFnFUE3NHCv5mPcdzKtsPFA4nt1lZZCs5cqpbm5s8c+fbqzx999t3u1e7STU02m+h5epUWlNLc2c+T4EfYc2QOQ1o4ApuBjpKKiMxufPTj9x7d5lpeWU9dYR3Vtddo9PMmgstKNHseNc6aZgwfd6+mnu+2V6e/4kTB6yj/T2NzIlJFTaG5tZtfhXR1mnnR2BDATjTFgBNMUFOYWsv3QdiaGJlpUay/wR5PV1c4cc/CgSzB2zjlW2CMWgvlngvhmw1BRiFlls6jeWk1reytAlypP6eYIYAreGBAi0xT88fU/cqjlUEfRhWAxblPwpyZaKT4/Id5tt1kVp77g9+79Dsj5o8/n9CGns4IVjB3iEuunoyOAKXhjQIjM93G87ThDC4dSu88FN9fuq+Vg80GAtE3slCwiE+L5Eaw2T9Q7gr37lVtWcsufbwFg9VurKSspIzcnN+0cAcwGbwwIkfk+SotKQWHXkV08W/+si2rNKaAwtzAtJ7OSiUWwxpea3TWsenUV5488H4CGYw1s3ruZaaOnUV1b3SUKNtWxHnyM+ENgcJOs2T4E7q2fcW94gRf65HecLF/jZFJTAw895Ko4DRvm8sGPGWOBeP0hcrT58fM/zqv7XuW+V+5j7plz02rOyBR8FGpq4Ac/gOeec9Vxzj7b5aGZ6OXQzPYhcF9qSwZt8L7v+5vhN3nkk48w7BvDGFQwiLGDx5KXk0dTaxMThk3gwQUPxl/4DMI3zRQUOAW/f39nb760FC64INkSphe+E8D9r9zP6UNPZ8rIKR37dhzZwYm2E11SCkPqzxmZiSaCmhr42tdg7Vr34OTnw5NPwpYt0NLijrEhcO+J5qK29L1LATiz9EzOHn42QwuHUpxfTA45HGw5mFyB0wBfmc+Y4bxpdu1y2w8ftoR4vSWYp+b0IafT2NTIk2+50tMP1T7EX/b+heLc4i6fSYf0GtaDj6C6GhoaXJRgsfd7irh8NLWBZMc2BO493bmotWt7R9X65tZm2rWdYYXDBl7ANMOPas3Jcffr0aPQ1uaWSy5xKa2tME1sBM0y5406jyfefIKGYw0A5Ofko6ocaD7A7iO708qrxhR8BHV1rqcezMo3aJB7eBoDKaODuWiCNVvNRS12/EmqI8ePcOj4IYpyixg7dCyThk/i7BFnJ1m6gUO/XgJVpT0fGMHdAd1S9aFTHFjV61MDnlxZQl1jHWUlribnmCFjKCkq4eiJoxzkIIMKBvHes97LK3tf4eVdL/P+Se9Pm/QapuAjKC935pjm5s4efEEB7N0L9fVu/bXXIDfX5fkwF7W+4Q+JAYYUDiGHHNq1nXFDx5Gbk5sRubhjRW4+1Kd5jZUrnTlxzx5ob3f3qQgMHw6XXurWQ6HOCOxeyyWC9vGzqU5vnATWfmZtl/VoScj68vsNBKbgI6ishPXr4Y033MTVsWPuARo+3KUJ3rfPZey76abO9AXBmq3+qw2NT40/JAa46IyLqN1Xy96je9lxeAe3X357Sk9cpQI1NbBqlTPDDBniRphHj8IZZzhvmg0bXF1WS4gXnUiFHOkE4AfizR0/lzFDxgCdicnSqU6rTbJGUFEBt97qbJjHjzuFfuaZ8OEPuwVckYXNXulwq9naN4J+8WOHjGXu+LnMP28+E0ITTLnHgD/BmpcHkya5JGOTJrmiNMePOzOjjSJjJ9IJYMZpM5g0fBIFuQVpXeXJevBRqKiAH/3IvV+4sHMiyyeowC1XfN8oLy0n3BTusi0dJq1SBX+CtbTUJRsrLoYRI5zX18UXu/vRlHvviHQCCOZOKi8tZ9GMRWnX+TAF3wM9KfBgrvhsL7IQKzW7a9h9ZDePv/k4ALsO76IorygtJq1SBf++nDwZnn3WbVN1dne7/+JDd15f6YSZaHrAT88aDruJLOialjWYK76+3nLF94Rv6yzMLeTSsy4F4Im3nuB42/GUjwpMJfxi8C+/7Hrw9fXOD/6CC+z+MzqxHnwPBNOz+maZyAcoWnY/IzqRYeAAHzz7g4SKQ6bce4k/TzhokOtYjBoFn/uc3YtGJ1mj4IcPH044HO75wBiYNi0+eVhCoRAHDhyIy7nShaC/sU86RASmGtXVLnXGzJmd28Jh894yupI1Cj4cDqecr2oiE3alKv7kqhXc7h/+JGsQ894yIjEbfA/U1Dhf94UL3avl9ugfkXUxgbR0P0s2fm3WIOa9ZURiCv4U+FGq4XDXKFVT8n0n0t8YsMnVPhA5+e+/t5qsRhBJJbPFzJkzdf369Qk5t4j02kRTVXWyi6S/7vfm+5ODpi8yZRrWBn1vg+D956cpaGmJXz4k+23SAxF5SVVnRtuXNTb4vhC0c27ZAi++6NKwFhTA0KGwaZPloDGSh3+f/fCH8Kc/uUCn6dPtXjQ6MRPNKfDtnFu2wOOPO3/joiLnnvbf/+3SCluZNCNZ+CbEDRtcriRwRWpaWuxeNBym4E+Bb+d85hmnxI8dc+u5uS7n9tatXY83LwZjIPHz0Rw/7lIVFBe7Dkhtrd2LhsMU/Cnwg5yOHHGZ+gAGD3YVc44ehZ07YffuzuPNi8EYSPxEd6WlLr01OAXf2Gj3ouFImIIXkTNEZLWIbBWRLSLyxURdK5FUVLhe0tChbmlqctvz8pypZu1aFyJuXgzGQOObEKdMcQq+qcktfj4auxeNRPbgW4GvquoUYBbwLyJyXgKvlxB8l8jGRmdzV4XWVrfk5blJ17VrLQeNMfD4JsSCApg1y207cMDVaLV70YAEKnhV3aWqG7z3h4GtwLhEXS8R+JNYY8a4QgrgekonTrgiC2PGwLhxTtlbmT5joAkmutuzx9ngp051hWkMAwbIBi8i44EZwPMDcb14EaxaP2QIjB7tbJyFhe5hGjnS+R6PGGEeC0ZyqKhwnYuSEpg2za1bQJ7hk3AFLyJDgN8AX1LVQ1H2Xyci60VkfUNDQ6LF6RX+JNbYsTB7Npx2mvOeUXWKf9cuV5+1tRU2bky2tEa24ndEzGXXiCShgU4iko9T7veratTbTVV/DPwYXCRrIuWJhWB04JtvOpPMOec4Jf+hD7ne0a5drmariOvN19c7z5qaGjPTGAPPxo3uvjx0yHVIRo1yReJ37nT7zXyYvSTSi0aA5cBWVf1uoq4TTyJzz4wb5wJHXn3V5ft49VUXRJKb63pKeXnOHt/S4l5/+MNkfwMj26ipgbfeck4AJSWwf78Lytu1C04/3cw12U4iTTTvBq4B3isiG73lAwm8Xr+JHOqefbYzzezc6XrpO3e6Yty5uU7ht7Q4k01+vus1Pfdcsr+BkW1UV7uJVVU32jx0yL33YzQ2bXL3qJlrspOsSTZGVWlizttfqhp7PiaDsYRWqVsXIBsL0qQjlmwMkJsP9ahIusseefy4i1zdsAHefrszqrWgoDPgqbgY5s+HH/2oFzKJoFW9/ioZR9Waqo7K9ZWTK7MudXB//uAi79mbbxZGjlQGD4bx4922gwdh2DBYs6Z/chrph6UqCNBdjm0/uOn4cTcMzslx21pa3DJ4sFPwn/98sr9BelGz2xmGw01hykrKCDeFWfbsso7tRs9EKwrf2urcelVdZGt7u1PwRvaRNT34WIgssF1eDosWwW23wY4drseen+8mVHNy3IOTk+P84+fMMU+F7qjZXUN1bfVJvfTqWmcY9sv3+a/VtdVZ14uPhe7qD1x5Jdx5p7tHAc47z92njY3Oq2biROcJZmQfpuAjqKg4WVGXl7sJ1GPHXM/o6FHnIpmb27lY7z06K7es5JanbuFE+wlGDR5F84lmlj27jMWzF0cttG0FuKPje3hF1h+48kpYtcoFOV18Mdxyi/OgmTXLrTc2Wl6abMZMNDFQWel6RIcOdUaytrS44W9bm8vFbb33k6nZXcMtf74FQRg1yCn3Lfu20NrWSnVtddRC21aAOzrdBTPdeWfX7dDV88tyJGU3puBjoKICbrrJKfmDB51yLylxD8/YsS7Bk/kZn0x1bTUn2k5QWlSKiFCcX0xRbhE7juygrrGuo9C2X4A73BS2Atzd4EdVByktdWaZ5mY3gfrQQ277kCEwYQLcfbebhDXlnr1kjZtkPNzxVq6Ef/1XN9k6dKhT8jk5cP75zsZZVTXwMqUy8Xb/y+S26onuPLyeftqZDEtK3OjyqaeEyZOVSy7pnUeXkb6Ym2ScmD8ffv7zrmHhkye7SVarnnMyX1/9dV7b/xqb926mKK+IorwiGpsbUZSffeRnNpHaCyornc0d3H3n29YnTIBXXjn5+BR1rTcGGDPR9JLp091y1VUwd64z0Vj1nOhUTq4kNyeXqaOnUpRXRMOxBlSUmy6+yZR7LwmmBg7a1ocMcZOpxcWu0wFuvaUlufIaqYH14HtBTY3Lux2sYF9U5HpSixYlW7rUo2JsBYtnL6a6tprCvELmnTUvKwOZ4kV3Hl7hsOtsgEs8VlTkMp8ahin4boj0OZ461bmjhUJw6aXuQXriCbjsMvNSOBUVYytMoSeQSNMNWIfD6MQUfBSCPseHDsFPfgL79rkHaPp0N8na2OjcJNes6aygY0reGGgig/PAOhxGJ+ZFEwXfY2HnTpd6taCgM0tfe7sr1Xf8uDv2xAm44goX7NTbByvTvWiMgae7e6q7KFgj/TmVF41NskbB9zl+8UWn3IuLXXBTW5tLLtbQ4F7z8ty+HTusgo6RukTWObAc8dlDVplo4uGX3drqXv1qOf77tWvd+5tvjv1coaBTs2EkiGAULHS+VldbLz7TyRoF3xtTiN/jeeIJp9Dz8txrSws0NgoiLh1rbq4z1Qwa5BI6zZhhwSVG8ok0x2zceLIiLy212I1swEw0UfAnri66yCUYa2115ftyc93+ceNc2oLmZrfv2DGorYUnn3TRroaRTCLNMW+95WoIB7HYjezAFHw3VFTAr38N99wDU6a4otr5+W7fxIkuUjAnpzNicNIk15O/5RazbRrJJTIp2fnnw+bNJ9c5sAyTmU/WmGh6S3CYO3euexiqq52Nfe5c1wM6cMCZbYqKXG4aVTcBa7ZNI5lEJiWbNMmNMkOhrnUO7B7NfEzBR+FUubfBrZeUuPJ9+fkwcqTb3tzsim+bbdNIBv7IceVKlx9p8uTOVBrTp/c+GZ6R/piCj0J3XgebN3euh0LOhbK01JXsa2pyCn7iRLNtGgOP3ykB5xRw8CA884yLwM7Ls8jWbMUUfBTq6lzPffduN3na2NiZ/x06e0IrVzqbe0OD67lPnOgeJrNtGgON3ykBePe7YetW2LvXufDefruZY7IVU/BRKC+HV1+FLVucfb2kxCn5gwfd/qB9fvZsN9Ha0mIRgkbi6E0Mx49+1Hns1q2unF8kFkGdHZiCj0JlJVxzjVPcRUXO9KLqhruR9nk/L7fl/zASSU8KubuCIKGQ2d6zGXOTjEJFBZx1lvN7/+tfOydThwxx+6PVxrQ0BUYyqazsdH80V0jDxxR8N/i1VnNzXe991y547DG3L1ptTPOcMZJJdwVBbFSZ3ZiJphsOHID9+13PvbXVpSRob3f7Ghu7DoUtKtBIBaIVBDGyG+vBd8Nrrzlf4uZmlxK4oMC5QwK8+aYNhQ3DSH2sB98Nqs63fdiwzhQFJ044ZT5unEUFGoaR+mRNwY9YiEc64SCp1LaGYWQmpyr4kbAevIjcDXwI2KuqUxN1nXgSVMg1NTB/fmel+txcZ6KZORPOOcdczwzDSH0SaYO/F7g8gedPKBUVcOutrjzf8OEuYdPMmRapahhG+pCwHryqPiUi4xN1/oFg/nzXW7daloZhpCM2ydoD5npmGEa6knQ3SRG5TkTWi8j6hoaGZItjGIaRMSRdwavqj1V1pqrOHDVqVLLFMQzDyBiSruANwzCMxJAwBS8iK4BngXNFpF5ErOSAYRjGAJJIL5oFiTq3YRiG0TMpFckqIg3A28mWowdGAvuSLUQGYe0ZX6w940s6tOeZqhp1AjOlFHw6ICLruwsLNnqPtWd8sfaML+nenjbJahiGkaGYgjcMw8hQTMH3nh8nW4AMw9ozvlh7xpe0bk+zwRuGYWQo1oM3DMPIUDJKwYvIEhHZIiI1IrJRRN4lIl8SkUGBY34vIsPidL0jPewfJiKf78N5q0Rkcd8liw/Wnn1HRFREvhNYXywiVT185koR+fcejpkrIo90s2+biIzsk8AxXj9ZWHv2jYxR8CIyG1dg5AJVrQAuA7YDXwI6FJKqfkBVDw6QWMOAXiukVMDas9+0AJW9URCqukpVv5lAmbpFRPKSef0YsPbsAxmj4IHTgH2q2gKgqvuA+cDpwGoRWQ2d/8oiMl5EakXkLhHZLCL3i8hlIvK0iLwmIu/0ju/S+/OOHR+8sIgMEZEnRGSDiLwiIld5u74JTPR6v9/2jr1BRF70esU3B86xRET+KiKPA+cmqpF6gbVn/2jFTdB9OXKHiIwSkd94cr8oIu/2tn9aRO703k8Ukee8/f8VMboZIiIrvfa+X6RLrckbROQFb5nknetMrz1rvNdyb/u9IvJd77f8n4jr3ysit4vIMyLypojM97bniMgPxI3sHhE3gpufgPaLxNqzL6hqRizAEGAj8CrwA+ASb/s2YGTguG246LTxuJvmHbg/upeAuwEBrgIe9I6vAhYHPr8ZGO+9P+K95gEl3vuRwOveecYDmwOffR/uJhXvmo8AFwMXAq/gesYl3ucXx7uNrD0HtP2OeNfeBpQCi4Eqb98DwBzvfTmw1Xv/aeBO7/0jwALv/T8H2mYu0AiUed/52cC5tgFLvPefAh7x3j8MXOu9Xxj4Le71rpMb5fr3Ar/2rnEe8Lq3fT7we2/7WCAMzLf2TM32zJiCH6p6REQuBN4DzAN+KT3bv95S1VcARGQL8ISqqoi8glMmsSLArSJyMdAOjAPGRDnufd7ysrc+BDgbGAr8VlWPebKs6sW1E4K1Z/9R1UMich/wBaApsOsy4LxAR7FERIZGfHw2cLX3/gFgWWDfC6paDyAiG3Ftu87btyLw+r3AufxCkz8DvhU4169Vta2br/CgqrYDfxERv/3neJ9pB3Z7vdUBwdqz92SMggfwGnYNsMZTKtf28JGWwPv2wHo7nW3TSldTVlGU83wSGAVcqKonRGRbN8cJ8A1V/b8uG0W+BKScv6q1Z1y4DdgA3BPYlgPMVtWgkqKrZeCUBNu5ja7PsXbznm62H43xOhLxmixuw9ozZjLGBi8i54rI2YFN03GJyw7jenR9ZRtwgXeNC4CzohxTCuz1lNE84Exve+S1/wgsFJEh3vnGicho4CngIyJS7PU8PtwPeeOCtWd8UNUDwK+AYLrsx4B/9VdEZHqUjz4HfNR7/3e9uOQnAq/Peu+fCZzjk3T2TvvCOuCjnu14DM7EMWBYe/aOTOrBDwHuEOey14qzu14HLAAeFZFdqjqvD+f9DfApb+j2Is4mHcn9wMMish5nt64FUNX94iYZNwOPquoNIjIFeNbrXRwB/kFVN4jIL73Pvg38uQ9yxhtrz/jxHQIKCGdi+F8RqcE9g0/h7MJBvgT8XES+CvwOZyeOhUIReR7XefNTdn8BuFtEbgAagM/05Ut4/Aa4FDd38irwfC9kixfWnjFikayGkYKIizVo8uYw/g43QXhVT58bCERkiDdHMwJ4AXi3qu5OtlynIlvbM5N68IaRSVwI3Om57B3EeWukCo94I7sC4JZUV+4eWdme1oM3DMPIUDJmkjVWRKRNXKDMJnGBNBcl8FqnDL3PBKw944e1ZXyx9sxOE02Tqk4HEJH3A98ALonlg97wTjyf1ZPWsxRrz/hhbRlfsr49s64HH0EJLnKs2/B4cSH4W0XkBzj/2/dErJ8h3YTL+4jIaSLylNeb2Cwi7xnA7ziQWHvGD2vL+JKd7ZnoEONUW3CBDBtxrneNuGAaOHV4fDswy9sXuR41XN7b54dDf5XOkOdcYGiy28HaM/UWa0trz3gv2W6imQ3cJyJTOXV4/Nuq+lzgHMH17sLlnwoc/yLObzYfF668Me7fKnlYe8YPa8v4kvXtmY0KvgNVfVZc+tFRwAfoPjw+Mvw4uB41XD7iOk95N9MHgZ+JyLdV9b54fY9Uwdozflhbxpdsbc+stsGLyGTcMGo/3YfH90R34fLB65zpnfsnwHK8UP1Mw9ozflhbxpdsbc9s7MEXiwuTB/ePfK2qtolI1PD4nlDVxyRKuDywN3DYXFxe6RPe/k/F4XukCtae8cPaMr5kfXtaoJNhGEaGktUmGsMwjEzGFLxhGEaGYgreMAwjQ8lKBS+uEMRaEcmNw7n+RlzOC7+IbpG4Ar2bxBXSDRaCXiYi7+3vNVOJeLSliHzSiwysEVeUeFpg35e9dtwsIitEpMjbnnFtCXFrz6u8ttwoIutFZE5g390isldcTv3gZ6w9uz9H+t6fyY42S8YC/AvwxTicJxd4Elc0d763TYAh3vt8XAJ/PxLuTOCxZH//VGtL4CIg5L2/Anjeez8OeAso9tZ/BXw6U9syju05hE4HigqgNrDvYpzr3uaIz1h7dn+OtL0/s7IHjyuz9ZCIzBWRNSKyUkRqReR+kdgLOQLX4yqydLhJqcPPLJfvLertexsYISJj4/M1UoJ+t6WqPqOqYW/1OVyFe588nLtbHjAI2Ol9JhPbEuLTnkfU0zLAYAI1Q1X1KeBAlM9Ye3ZDOt+fWafgRaQAmKCq27xNM3DlvM4DJgDv9o77njfEjVz+3ds/DvgI8KMo18j1/G/3An9S1ecDuzf410h34tWWESwCHgVQ1R3AMqAO2AU0qupjgWMzpi0hvu0pIh8RkVpcebpYi1tYe2bY/ZmNgU4jcRVdfF5Q1XoATymPB9ap6pd7OM9twL+pC5zoskNV24Dp4qq0/FZEpqqqb/PcC5zev6+QMsSrLfE+Mw/3AM3x1kPAVbjC3AeBX4vIP6jqz72PZFJbQhzbU1V/i7v3LgZuAS6L4frWnqcgHe/PbFTwTXTmnQBoCbxvw2sTEfkeEK2o9C9U9ZvATOAXnnIfCXxARFpV9UH/QFU9KCJrgMtxRXXxrt0Ul2+SfOLVlohIBXAXcIWq7vf2Xwa8paoN3jHVOHuo/wBlUltCHNvTR11ulIkiMlJV9/VwfWvPrqT9/Zl1Cl5Vw54JpaiH4075r66qZ/nvReRe4BFVfVBERgEnPOVejLsJ/ifw0XOAX/f5C6QQ8WpLESkHqoFrVPXVwK46YJZ4BZNx1efXB/ZnTFtCXNtzEvCGqqqIXICr9bn/VJ/xsPaMQjrfn1mn4D0eww2zWhNw7tOAn4pzy8oBfqWqjwCISyE6ia43QboTj7b8f8AI4AfeiKhVVWeq6vMishJny2zFpWn9MWRsW0J82vOjwKfE5UNpAj7hT7qKyApcvpSRIlIPfF1Vl1t7npK0vT+zMheNiMwAvqKq1wzwdT8CXKCqNw3kdROJtWV8sfaML9nenlnnRQOgqi8DqyUOgU69JA/4zgBfM6FYW8YXa8/4ku3tmZU9eMMwjGwgK3vwhmEY2YApeMMwjAzFFLxhGEaGYgreSBoiskRcJj4/8+G7vO1f8nyL/eN+70UFx+OaR3rYP0xEPh9YP91zhYvHtdeIyMx4nMswYsEUvJEURGQ28CGcK1kFLiBsu7f7S7jETQCo6gdU9eAAiTYM6FDwqrpTVecP0LXjQhI8RowUxRS8kSxOA/apaguAqu5T1Z0i8gVc/o7VIrIaQES2ichIERkvLhPgXeLyb98vIpeJyNMi8pqIvNM7vkpEFvsX8o4dH7y4iAwRkSdEZIOIvCIiV3m7vglM9EYU3/auudn7TJGI3OMd/7KXmwQR+bSIVIvIHzw5vhVrI3jn/7MnxwYRucjb/rOATHjf9UovMvPbIvKiN/L5rLd/roisFpEHgFd69UsYmctA5ye2xRbVjpzlG4FXgR8AlwT2bQNGRq7jkkO1Au/AdU5eAu7G5eC/CnjQO74KWBz4/GZgvPf+iPeaB5R470cCr3vnGU8gV3pwHfgqcI/3fjIuVL0I+DTwJlDqrb8NnBHlO68BZkZsGwQUee/PBtZ77y8JfJ9SXN7xPOA64D+97YW4SMmzcNGpR4Gzkv3b2pI6i/XgjaSgLmf+hTiF1QD8UkQ+HcNH31LVV1S1HdgCPKGqiuu1ju+FCALcKiI1wOO44g1jevjMHOBnnvy1OEV+jrfvCVVtVNVm4C+4gg+xkA/8REReweUtOc87/1pgkoiMBhYAv1HVVuB9uDQEG3HFZEbg/hjAZUt8K8brGllAtuaiMVIAdWmV1wBrPAV3LXBvDx8LZgRsD6y303k/t9LV/Bgt2dQngVHAhap6QkS2dXNckFMViIiaqTAGvgzsAabhZG4O7PuZJ+ff0ZnTXYDrVfWPXQQTmYvrwRtGB9aDN5KCiJwrImcHNk3H9YgBDgND+3H6bbiydHjZFM+KckwpsNdT7vPo7HGf6tpP4RQuInIOUA78tR9y+nLs8kYk1+DKQPrci5twRlW3eNv+CHzOS2aFiJwjIoP7KYORoVgP3kgWQ4A7PPfHVpwN/Dpv34+BR0Vkl6pGy9PdE7+h04zxIs7OH8n9wMMish43F1ALoKr7vUnbzbjKPf8b+MwPgB95o41WXP3NFulVlUd+52V5BHgW+BrwGxH5GLCaQC9cVfeIyFbgwcDn78KZojaIu3ADcHVvBDCyB8tFYxgpihcL8ArOlbQx2fIY6YeZaAwjBRGRy3CjijtMuRt9xXrwhmEYGYr14A3DMDIUU/CGYRgZiil4wzCMDMUUvGEYRoZiCt4wDCNDMQVvGIaRofx/hjTL/eFsyYUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L4 home\t 3.3740107142857148 +/- 0.1903230577341382 n = 28\n",
      "L2/3 home\t 2.0320232558139533 +/- 0.10822365744798969 n = 43\n",
      "L4 neigh\t 3.4806678571428575 +/- 0.15143091996001262 n = 28\n",
      "L2/3 neigh\t 2.6919000000000004 +/- 0.1894841070874011 n = 31\n"
     ]
    }
   ],
   "source": [
    "# Actual latency of neighboring barrel activation\n",
    "# To exclude direct stim, from each row, exclude the lowest latency if it is less than 1ms\n",
    "# separate by within barrel latency and crossing latency\n",
    "# latency plots\n",
    "plt.clf()\n",
    "i = 1\n",
    "xtick_labels = []\n",
    "df_datasets = []\n",
    "df2 = df\n",
    "\n",
    "act_lat_23_home = []\n",
    "act_lat_4_home = []\n",
    "act_lat_23_neigh = []\n",
    "act_lat_4_neigh = []\n",
    "df_23 = df[df['Stim Location'].str.contains(\"2/3\")]\n",
    "df_4 = df[df['Stim Location'].str.contains(\"4\")]\n",
    "\n",
    "# which is home barrel?\n",
    "for home_barrel in [1,2,3]:\n",
    "    df_home_23 = df_23[df_23['Stim Barrel'] == home_barrel]\n",
    "    df_home_4 = df_4[df_4['Stim Barrel'] == home_barrel]\n",
    "    home_snr_label = 'Actual latency barrel ' + str(home_barrel)\n",
    "    neigh_barrels = [1,2,3]\n",
    "    neigh_barrels.remove(home_barrel)\n",
    "    neigh_snr_labels = ['Actual latency barrel ' + str(neigh_b)\n",
    "                         for neigh_b in neigh_barrels]\n",
    "    \n",
    "    act_lat_23_home.append(df_home_23[(df_home_23[home_snr_label] >= 1)][home_snr_label])\n",
    "    for nsl in neigh_snr_labels:\n",
    "        act_lat_23_neigh.append(df_home_23[(df_home_23[nsl] >= 1)][nsl])\n",
    "        \n",
    "    act_lat_4_home.append(df_home_4[(df_home_4[home_snr_label] >= 1)][home_snr_label])\n",
    "    for nsl in neigh_snr_labels:\n",
    "        act_lat_4_neigh.append(df_home_4[(df_home_4[nsl] >= 1)][nsl])\n",
    "           \n",
    "act_lat_23_home = pd.concat(act_lat_23_home)\n",
    "act_lat_4_home = pd.concat(act_lat_4_home)\n",
    "act_lat_23_neigh = pd.concat(act_lat_23_neigh)\n",
    "act_lat_4_neigh = pd.concat(act_lat_4_neigh)\n",
    "\n",
    "jittered_xs_1 = np.random.normal(i, 0.04, act_lat_4_home.shape[0])\n",
    "jittered_xs_2 = np.random.normal(i+1, 0.04, act_lat_23_home.shape[0])\n",
    "jittered_xs_3 = np.random.normal(i+2, 0.04, act_lat_4_neigh.shape[0])\n",
    "jittered_xs_4 = np.random.normal(i+3, 0.04, act_lat_23_neigh.shape[0])\n",
    "plt.scatter(jittered_xs_1, act_lat_4_home, alpha=0.5, color='blue', label=\"L2/3\")\n",
    "plt.scatter(jittered_xs_2, act_lat_23_home, alpha=0.5, color='green', label=\"L4\")\n",
    "plt.scatter(jittered_xs_3, act_lat_4_neigh, alpha=0.5, color='blue')\n",
    "plt.scatter(jittered_xs_4, act_lat_23_neigh, alpha=0.5, color='green')\n",
    "\n",
    "xtick_labels.append(\"Stimulated\\nBarrels\\n(n=\" + str(len(act_lat_4_home)) +\")\")\n",
    "df_datasets.append(act_lat_4_home)\n",
    "xtick_labels.append(\"Stimulated\\nBarrels\\n(n=\" + str(len(act_lat_23_home)) +\")\")\n",
    "df_datasets.append(act_lat_23_home)\n",
    "xtick_labels.append(\"Neighboring\\nBarrels\\n(n=\" + str(len(act_lat_4_neigh)) +\")\")\n",
    "df_datasets.append(act_lat_4_neigh)\n",
    "xtick_labels.append(\"Neighboring\\nBarrels\\n(n=\" + str(len(act_lat_23_neigh)) +\")\")\n",
    "df_datasets.append(act_lat_23_neigh)\n",
    "\n",
    "#p_val = kstest(act_lat_23, act_lat_4).pvalue\n",
    "#print(\"KS test p-value \", p_val)\n",
    "\n",
    "plt.boxplot(df_datasets)\n",
    "plt.xticks(range(1, len(xtick_labels)+1), xtick_labels)\n",
    "plt.title(\"Actual Latency of Activation\")\n",
    "plt.xlabel(\"Stimulation Layer\")\n",
    "plt.ylabel(\"Latency (ms)\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print(\"L4 home\\t\", np.mean(act_lat_23_home), \"+/-\", sem(act_lat_23_home), \"n =\", len(act_lat_23_home))\n",
    "print(\"L2/3 home\\t\", np.mean(act_lat_4_home), \"+/-\", sem(act_lat_4_home), \"n =\", len(act_lat_4_home))\n",
    "print(\"L4 neigh\\t\", np.mean(act_lat_23_neigh), \"+/-\", sem(act_lat_23_neigh), \"n =\", len(act_lat_23_neigh))\n",
    "print(\"L2/3 neigh\\t\", np.mean(act_lat_4_neigh), \"+/-\", sem(act_lat_4_neigh), \"n =\", len(act_lat_4_neigh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1c2476e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}