diff --git a/.gitignore b/.gitignore index 7290c69..a790567 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,5 @@ data/*.csv docs/sphinx/_build/ docs/sphinx/auto_examples/ -scripts/*.pdf \ No newline at end of file +scripts/*.pdf +docs/figures_notebooks/Figure_1_C.ipynb diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..26d3352 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/.idea/APE_paper.iml b/.idea/APE_paper.iml new file mode 100644 index 0000000..7a1bce4 --- /dev/null +++ b/.idea/APE_paper.iml @@ -0,0 +1,12 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/Project_Default.xml b/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 0000000..3dce9c6 --- /dev/null +++ b/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,12 @@ + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..e7edcdc --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,10 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..186c2c8 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..35eb1dd --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/README.md b/README.md index e6acae9..7bbe773 100644 --- a/README.md +++ b/README.md @@ -97,3 +97,5 @@ git clone https://github.com/HernandoMV/APE_paper.git

+PLACEHOLDER FOR YVONNE FIGURE + diff --git a/data/README.md b/data/README.md deleted file mode 100644 index 65556e2..0000000 --- a/data/README.md +++ /dev/null @@ -1,4 +0,0 @@ -# This directory will be filled with data when running the notebooks. -#### If there are problems fetching the data automatically, get it from here: - -https://zenodo.org/record/7261639#.Y1_Sm9LP2Xk \ No newline at end of file diff --git a/data/TS20_20230512_RTC_restructured_data.pkl b/data/TS20_20230512_RTC_restructured_data.pkl new file mode 100644 index 0000000..02a7306 Binary files /dev/null and b/data/TS20_20230512_RTC_restructured_data.pkl differ diff --git a/data/TS20_20230512_RTC_smoothed_signal.npy b/data/TS20_20230512_RTC_smoothed_signal.npy new file mode 100644 index 0000000..0374000 Binary files /dev/null and b/data/TS20_20230512_RTC_smoothed_signal.npy differ diff --git a/docs/figures_notebooks/Figure_S5_TU.ipynb b/docs/figures_notebooks/Figure_S5_TU.ipynb new file mode 100644 index 0000000..95e6a31 --- /dev/null +++ b/docs/figures_notebooks/Figure_S5_TU.ipynb @@ -0,0 +1,381 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "0a0160a8-61a0-4d9b-92bb-8cd1e0ed71b2", + "metadata": {}, + "outputs": [], + "source": [ + "# This notebook compares the dopamine signal aligned to movement during the CoT task to the\n", + "# dopamine signal evoked by the same high and low frequency sounds played while the mouse is freely moving." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75f5c411", + "metadata": {}, + "outputs": [], + "source": [ + "# run this on Colab\n", + "!rm -rf APE_paper/\n", + "!git clone https://github.com/HernandoMV/APE_paper.git\n", + "%cd APE_paper/docs/figures_notebooks\n", + "!git checkout YvonneJohansson" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "3e40afcb-ef99-4802-ae0a-5bbceacd3790", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "import pandas as pd\n", + "import urllib.request\n", + "from os.path import exists\n", + "import numpy as np\n", + "import pickle\n", + "from scipy.signal import decimate\n", + "\n", + "import warnings\n", + "warnings.filterwarnings('ignore')" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "88c23d07-da61-4490-b9de-1b3e4835b6b2", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import os\n", + "sys.path.insert(1, os.path.dirname(os.path.dirname(os.path.abspath(os.curdir))))\n", + "from scripts import YJ_analysis_utils as yj_utils\n", + "from scripts import yj_plotting as yj_plot" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "1e0c1c17-8a9c-4ae0-80dc-8ac0dfdb1f9f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Downloading data...\n" + ] + } + ], + "source": [ + "# Supplementary Figure 5T & 5U:\n", + "\n", + "# Get dataset information:\n", + "dataset_name = 'DataOverview_SF5TU.csv'\n", + "zenodo = \"https://zenodo.org/records/10803926/files/\"\n", + "url = zenodo + dataset_name\n", + "\n", + "#dataset_name = 'DataOverview_SF5TU.csv'\n", + "dataset_path = '../../data/' + dataset_name\n", + "\n", + "if not exists(dataset_path):\n", + " print('Downloading data...')\n", + " urllib.request.urlretrieve(url, dataset_path)\n", + "else:\n", + " print('DataOverview already in directory')\n", + "\n", + "#print(dataset_path)\n", + "info = pd.read_csv(dataset_path)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "a2748c53-a8c0-4e5d-af85-1609cf0ec66d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Unnamed: 0AnimalIDDatefiber_sideprotocol1protocol2
01TS320230203right2ACRTC
12TS2020230512leftpsychometricRTC
23TS2120230510leftsilenceRTC
34TS2620230929right2ACRTC
45TS3320231106rightSORRTC
\n", + "
" + ], + "text/plain": [ + " Unnamed: 0 AnimalID Date fiber_side protocol1 protocol2\n", + "0 1 TS3 20230203 right 2AC RTC\n", + "1 2 TS20 20230512 left psychometric RTC\n", + "2 3 TS21 20230510 left silence RTC\n", + "3 4 TS26 20230929 right 2AC RTC\n", + "4 5 TS33 20231106 right SOR RTC" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "info\n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "b0e7f4c1-cd09-4ac3-bc65-304cded04237", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Downloading TS3_20230203_smoothed_signal.npy\n", + "Downloading TS3_20230203_restructured_data.pkl\n", + "Downloading TS3_20230203_RTC_smoothed_signal.npy\n", + "Downloading TS3_20230203_RTC_restructured_data.pkl\n", + "Downloading TS20_20230512_smoothed_signal.npy\n", + "Downloading TS20_20230512_restructured_data.pkl\n", + "Downloading TS21_20230510_smoothed_signal.npy\n", + "Downloading TS21_20230510_restructured_data.pkl\n", + "Downloading TS21_20230510_RTC_smoothed_signal.npy\n", + "Downloading TS21_20230510_RTC_restructured_data.pkl\n", + "Downloading TS26_20230929_smoothed_signal.npy\n", + "Downloading TS26_20230929_restructured_data.pkl\n", + "Downloading TS26_20230929_RTC_smoothed_signal.npy\n", + "Downloading TS26_20230929_RTC_restructured_data.pkl\n", + "Downloading TS33_20231106_smoothed_signal.npy\n", + "Downloading TS33_20231106_restructured_data.pkl\n", + "Downloading TS33_20231106_RTC_smoothed_signal.npy\n", + "Downloading TS33_20231106_RTC_restructured_data.pkl\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAC+CAYAAAAWcCl5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABSb0lEQVR4nO29Z3Rc2XWg+51buVDIOUeCAJFIkGAOIMEmqZbUCm6pZY9tBVsaazxj95v33tgeL1uyZ2kshzf2yHJqy7JlSR5FS261ZCYws5s5ggkkmEDkHCqH835UEQRAAATJKhTC+daqhap7T927q1Dn7rvPTkJKiUKhUCgU4UKLtgAKhUKhWFwoxaJQKBSKsKIUi0KhUCjCilIsCoVCoQgrSrEoFAqFIqwoxaJQKBSKsBJVxSKE+LoQolsI0TTN/nohxJAQ4lLo8ftzLaNCMd9Q80Yx39FH+fz/BHwV+OcZxhyXUn5gbsRRKBYE/4SaN4p5TFQtFinlMaA/mjIoFAsNNW8U852F4GPZIIS4LIT4dyFERbSFUSgWCGreKKJGtJfCnsUFIF9KOSqEeBX4MbBs8iAhxOeAzwGsWLFi9bVr1+ZUyLnE4XAAYLVaoyzJokJEW4Awo+aNYi6Ydt7Ma4tFSjkspRwNPf8ZYBBCpEwx7i0p5Rop5RqLxTLncioU8wk1bxTRZl4rFiFEhhBChJ6vJShvX3Slii4ejwePxxNtMRTzGDVvFNEmqkthQoj/A9QDKUKIR8AXAAOAlPJvgdeBzwshfIAT+IRc4uWYzWZztEVQRBk1bxTznagqFinlzz9j/1cJhlUqQvh8vmiLoIgyat4o5jvz3XmvmEQgEIi2CAqFQjEjSrEsMNRSmEKhmO8oxbLAeBxubDQaoyyJQvE0bQODfO/sBQ5cu4nT48FiNPJKRRkfr6slOzEh2uIp5gilWBYYSqEo5iunW+7xhR//FG8ggD+0ZOvweHjnchP7rl7nDz78ftYVF0ZZSsVcMK/DjRULn87OTj7xiU9QXFzM6tWrefXVV3nrrbf4wAeer4zVr/7qr3L9+vUISal4WdoGBvnCj3+Ky+cbUyqP8QcCuHw+vvDjn9I2MBgdARVzilIsC4yFlMcipeQjH/kI9fX1tLS0cP78ef7oj/6Irq6u5z7W1772NVasWBEBKRXh4HtnL+B9RmCJNxDg+2cvzJFEimiiFMsCw2KxsFCypA8fPozBYODXfu3XxrbV1NSwZcsWRkdHef311ykrK+M//If/wOM0i8bGRlatWkVVVRWf+cxncLvdANTX13Pu3DkA9u7dS21tLTU1NTQ0NABgt9v5zGc+w9q1a1m1ahX/9m//Nsefdmlz4NrNCZZKfMAHk1Jn/IEAB67dnGvRFFFA+VgWGI+tFZPJNOv31P/xX0REliO/9eaM+5uamli9evWU+y5evMi1a9fIyspi06ZNnDx5kjVr1vCpT32KxsZGSktL+eVf/mX+5m/+hjfffHKenp4ePvvZz3Ls2DEKCwvp7w8W+f3Sl77Ejh07+PrXv87g4CBr165l586dxMTEhOvjKmbAOcmKlkCu30OrfuLv1LFArG3Fy6EsFkVUWLt2LTk5OWiaxsqVK7l//z63bt2isLCQ0tJSAD75yU9y7NixCe87deoUW7dupbAw6AROSkoCYP/+/Xz5y19m5cqV1NfX43K5ePjw4dx+qCWMZVJQybCmxy40snzuCdutKvhkSaAslgXGi0SFPcuyiBQVFRX84Ac/mHLfeItLp9O9dEUBKSU//OEPWb58+UsdR/FivFJRxjuXmyYsh/XrDKT5PaT5PXTrjOg0jVcqyqIopWKuUBbLAsPpdOJ0OqMtxqzYsWMHbrebt956a2zblStXOH78+JTjly9fzv3797lz5w4A3/zmN9m2bduEMevXr+fYsWPcu3cPYGwpbPfu3fzlX/7lmK/m4sWLYf88iun5eF0tBu3py0m3zoheSpL8XgyaxsfqaqMgnWKuUYplgWE0GhdMLosQgh/96EccPHiQ4uJiKioq+J3f+R0yMjKmHG82m/nHf/xHPvaxj1FVVYWmaRMc/wCpqam89dZbfPSjH6WmpoY33ngDgN/7vd/D6/VSXV1NRUUFv/d7vxfxz6d4QnZiAn/w4fdj1uvRTVIwXUYLCQLe3FSnkiSXCGKxFT1ds2aNfBw9tBhRjb4iwmJr9PXchGvetA0M8v1Q5r3D48Eayrx/fc0qbl++REVFBWlpaWGQWDEPmHbeKMWywBgcHAQgISEhqnIsMpRimYN5I6XkwIED1NXVkZiYGNFzKeaEhdlBUvE0VqtVWSuKBYkQgoaGBk6dOsXo6Gi0xVFEEKVYFhgulwuXyxVtMRSKF0Kn09HQ0MDRo0fV73gRoxTLAkPTNLQpom8UioWC0Whkx44dHDp0CK/XG21xFBFAXaEWGHq9Hr1epR8pFjYWi4XNmzfT2NiomtctQpRiWWCopTDFYiEuLo66ujoOHTrEYgsiWuooxbLAWEh5LBCs4VVRUUF1dTUrV67k9OnTETvXF7/4Rf7sz/4sYsdXhJ/k5GQqKio4duyYUi6LCLWmoogY7733Hu+88w4XLlzAZDLR29u7YEr+K+aOzMxM3G43p0+fZv369dEWRxEGlMWywFhI/Vg6OjpISUkZqwuWkpJCVlbWtKXxCwoK6O3tBeDcuXPU19cDQUvkM5/5DPX19RQVFfGVr3xl7Bxf+tKXKC0tZfPmzdy6dWtuP6DiKdoGBvnz/Yd49c//mu1//Be8+ud/zZ/vP/TMBl8FBQUkJCSoUjyLBGWxLDBepAz8o9+eujbXy5Lz5S0z7t+1axd/+Id/SGlpKTt37uSNN95g3bp1zyyNPxU3b97k8OHDjIyMsHz5cj7/+c9z5coVvvOd73Dp0iV8Ph+1tbXTlulXRJ6XbU1cVlbG5cuXuXHjBuXl5XMltiICRNViEUJ8XQjRLYRomma/EEJ8RQhxRwhxRQix5CvYORyOsbIu8x2bzcb58+d56623SE1N5Y033uDv/u7vnlkafyre//73YzKZSElJIS0tja6uLo4fP85HPvIRrFYrcXFxvPbaa5H+SIppCFdr4pqaGoaHh8eKjCoWJtG2WP4J+Crwz9Psfx+wLPRYB/xN6O+SRafTPfd7nmVZRBKdTkd9fT319fVUVVXxV3/1V9OO1ev1Y6GnkyPfwl1mXxFenqc18Zu7dsw4bt26dRw7dgyTyURWVlY4xVTMEVG1WKSUx4D+GYZ8CPhnGeQUkCCEyJwb6eYnCymP5datW9y+fXvs9aVLlyguLp62NH5BQQHnz58H4Ic//OEzj79161Z+/OMf43Q6GRkZ4Sc/+UkEPoViNkxuTZzi96J7idbEW7Zs4fr162M+N8XCYr4777OB1nGvH4W2TUAI8TkhxDkhxLmenp45Ey4aLKQ8ltHRUT75yU+yYsUKqquruX79Ol/+8penLY3/hS98gd/8zd9kzZo1s7LMamtreeONN6ipqeF973sfdXV1kf5IimmY3Jp4SNNR6HMRF5hoWc62NbEQgh07dnDu3DmGhobCJqdiboh6dWMhRAHwjpSycop97wBfllKeCL1uBH5LSjltGdbFXt1Ylc2PCKq68UvOm1f//K+fVhpSkun3YkDSqjMihSDGaOSn/9d/mvVxfT4f+/fvp76+Xv3m5x8LtrpxG5A77nVOaJtCoZhHvFJR9lSDL4SgQ2+kV9NT4nMRg3zu1sR6vZ6GhgYOHz68YMLsFfNfsbwN/HIoOmw9MCSl7Ii2UNFkujwWt8+HRzm0FVFiutbEAA5Nxx29mcyAl1Vxz291mEwm6uvrOXjwoAraWCBEO9z4/wDvAcuFEI+EEL8ihPg1IcTjfrQ/A+4Cd4C/B2ZvQy9SbDYbNpttwjav38/Prl3j369fx+f3R0kyxVJmptbEOk3DZDDw2Y9/jJzUVPbv3//cfsKYmBg2btyoilYuEKIaXiSl/Pln7JfAr8+ROAuCxw2SxneQ7Bgawu7xgJS0DgxQmJISJekUS5l1xYX8w2d+ccrWxB+rqx3rd5+ZmcmRI0dYsWIFeXl5sz5+QkICtbW1HD58mB07diDEkneNzVtm5bwXQmwGlkkp/1EIkQrYpJTzMoNpsTvv7XY7MDED/+Tduzzs70cTghSbjYbly6Ml3kIlIlcoNW+mR0rJxYsXcTgcbNiw4bnysx49esT9+/fZvHlzBCVUzIIXd94LIb4A/BbwO6FNBuBb4ZFL8bxM1eira3gYk16P2WCga3hYLYfNA9S8mRkhBLW1tZSVlbF37176+2dKZ5tITk4OmZmZnD17NoISKl6G2fhYPgK8BtgBpJTtQGwkhVJMj9vtHivaCODx+XB4veg1DU0IJDDgdEZPQMVj1LyZBSkpKezZs4empiYuX74869L5xcXFWK1Wrly5EmEJFS/CbBSLJ+TrkABCiOevgqgIGyaTaUJ5kxG3GwFj680S6Astlymiipo3s0Sn07F161bi4uLYv38/zlneGFVUVODz+Whubo6whIrnZTaK5XtCiL8jWE7ls8BBghFaiigQCAQmRMWMjrNeAAyaRufw8FyLpXgaNW+ek8LCQrZs2cLRo0d58ODBrN5TW1tLb2/vrMcr5oYZo8JE8Db4u0AZMAwsB35fSnlgDmRTTIHX653wesTlmrB8YNTr6RkdRUqpomaihJo3L47VamX37t1cunSJhw8fsnHjxmc69jds2MDRo0cxmUxkZGTMkaSKmZhRsUgppRDiZ1LKKkBNinnA5ByWQacT/Thnvk4I7D4fTq8X6wJqYbyYUPPm5RBCsGrVKvr6+ti7dy8bNmwgKSlpxvHbtm3jwIEDGI3GGccq5obZLIVdEEKo6n7zhNHR0bFcFoBhlwv9uDs6IQQCGFIO/Gij5s1LkpycPObYv3Tp0oyOfSEEDQ0NnD59mpGRkTmUUjEVs1Es64D3hBAtoWZbV4UQKhQjShiNRozjLBG72z3BYoFgjkD/AmkGtohR8yYMPHbsJyQkPNOxr9Pp2LlzJ8eOHZt1AIAiMswm8353xKVQvBD+QAC3z4dJr8fp8TLsdJIWF4tBp6N7dJSKaAu4tFHzJowUFBSQnp7OsWPHWL58OQUFBVOOMxgMNDQ00NjYyK5duzAYDHMrqAKYhWKRUj4QQtQAj9sQHpdSXo6sWIrpeFyA0mq14vR6EUJgd3v45onTODxeqnOz2F6xnF7lwI8qS3XetA0M8r1QSRenx4MlVNLl4+NKurwoFouFXbt2ceXKFY4fPz6tY99sNrN161YOHjzIrl27XqjrquLlmE3m/W8C3wbSQo9vCSH+S6QFU0yN2WzGbDYD4PJ6EcCZlvs4PMFosSut7fQOj+D2+XCpSrBRYynOm9Mt9/iVr3+Ldy434fB4kAQbe71zuYlf+fq3ON3y8tVshBDU1NRQUVHBvn376Ovrm3JcbGws69ato7GxcdZJl4rwMRsfy68A66SUvy+l/H1gPfDZyIqlmA6fzzdWOtzl8+GXklsd3QBjd4RNjzoQBCPGFFFjSc2btoFBvvDjnwZ/k5OqD/sDAVw+H1/48U9pGxgMy/mSkpLYvXs3169f5+LFi1Mqj6SkJKqrqzly5IhSLnPMbBSLAMYXn/KjOu5FDb/fjz9UC8zl9TIwasfh8RBrNrGjohSA5s5u/IEAAyoDP5osqXnzvbMX8D6jnL03EOD7Zy+E7Zw6nY4tW7aQlJTEvn37xrqrjicjI4OioiJOnToVtvMqns1sFMs/AqeFEF8UQnwROAX8Q0SlUkyL1Woda9Fq93joHQkqj8yEeFJjbcRZzDg9XobsTrpU2GU0WVLz5sC1mxMsFaMMwCQrwR8IcODazbCfOz8/n23btnH8+HHu3Xt6uS0/P5/k5GQuXAifUlPMzDMVi5TyfwGfBvpDj09LKf8iwnIppsFut4+Vzh91u+kdCea0ZCTEIYSgICUZgPbBobEMfMXcs9TmjXNSV9OYgJ8Sn4v4wEQ/nyNC7YUfO/aHh4c5evTomFX/mNLSUvR6PdeuXYvI+RUTmY3zfj1wW0r5FSnlV4AWIcS6yIummIrxeSx2t5u+0aCSSY8PFs4tSA1mHT/o7cfj9wcbgCnmnKU2byyTqjwM6Azc0ZsxSkmJ10lcSMFEshrEY8d+VVUVe/fupbe3d8L+6upqHA4HLS0tEZNBEWQ2S2F/A4yOez0a2qaIMk6vl0F70EGfFGr8lZechBCCjsFhvD6fcuBHjyU1b16pKHuqJTFC0BNSMGYZYJnPxY78nIjLkpSUxJ49e7h58yYXLlyYYLXX1dXR0dHBo0ePIi7HUmZWzns57j8jpQwQ5ZbGSxmPxzOWy9Jvd+D2+TDodMSYgneCJoOezIQ4pJS0DwzRo/ws0WJJzZuP19VimKxYHiME3TojraYYVmcEe953dXVFVB6dTsfmzZtJSUl5yrG/adMmbt26RU9PT0RlWMrMRrHcFUL8hhDCEHr8JnA30oIppuZxHovP72fAHpwsCTGWCYmQ+SnB5bDOwWHah4aiIqdiac2b7MQE/uDD78es1z9lueg0DbNezx985APUb9xAQ0MDbW1tHDhwIOIX97y8POrr6zlx4gR37wa/fiEE27dv58KFCwwODkb0/EuV2SiWXwM2Am3AI4I1kD4XSaEU0/M4j8Xt9zPsCC5zJcZYJ4x57MBv7RtgwOnEo1oVR4MlN2/WFRfyP19/jbykhAnb85IS+J+vv8a64kIgaE3U1tayfft2Hjx4wMGDB5+rNfHzYjab2bVrF3a7naNHj+Lz+dA0jYaGBt59992xYBhF+JhNVFi3lPITUso0KWW6lPIXpJTdcyGc4mkeN/py+3wMOV0AJFonKpaM+FhMej2DDiejThd9o6NTHUoRQZbivDndco///oO3edg/OGH7w/5B/vsP3n4q816v17NmzRq2bdvG7du3aWxsjKgFUVVVRXV1Nfv376enpwe9Xs/OnTs5cuTIhHbfipdnNlFhfyKEiAuZ841CiB4hxC+G4+RCiD1CiFtCiDtCiN+eYv+nQue7FHr8ajjOu5B5vBTm9vkYcQUVS4LVgj8QYMjlwu52o2kaucmJALT1D9KuOkrOOZGcN/ORl8m8NxgMrFu3ji1btnDz5k0OHTrEcIR+s4mJiezZs4fm5mbOnz+PwWBg+/btNDY2jlW0ULw8s1kK2yWlHAY+ANwHSoD/92VPLITQAX8FvA9YAfy8EGLFFEO/K6VcGXp87WXPu9BxOBw4HEGnvd0VvMuymU2Mut1UZmZi0Otx+3xjfpaOwWHu9/WpfJa5JyLzZr4Sjsx7o9HI+vXr2bRpE01NTRw+fDgivVU0TWPTpk2kpaWxb98+pJRs2rSJgwcPTmj7rXhxZqNYHkeyvB/4vpQyXN7gtcAdKeVdKaUH+A7woTAde9HyOI/F7fVidwejwyxGA0a9nqqsLNbk5uLyeikIKZZH/QPY3W4GVH+WuSZS82ZeMjnzPtXvxTyF5TKbzHuTycTGjRvZsGEDly9f5ujRoxHxg+Tm5lJfX8+7775Lb28va9as4dChQ+omLAzMRrG8I4S4CawGGoUQqYArDOfOBlrHvX4U2jaZnws1SvqBECJ3qgMJIT4nhDgnhDi3VEII7W43jpBi0ek0CpOT0WkaOYmJGPV6bGYT8VYLLq+PvlE79yLoHFVMSaTmzbxkcuZ9v6YnOeClyOvCGngSPPI8mfdms5nNmzdTV1fH+fPnI9LAy2w288orr+BwOGhqaqK0tJTjx4+H9RxLkdk473+bYHTLGimlF3Awd5bFT4ACKWU1wd7h35hGxreklGuklGtSU1PnSLTo8DiPpc9uxxcIYNDpMOr1ZMXHA8HQztK0NOxeL/nJQaule2iE293deFV02JwR5Xkz50zOvPcLQZvexD29ifiAn+JQ9v2LZN5brVa2bt1KbW0tp0+f5sSJE7hc4dXRVVVVrFy5kqamJmw2G6dPnw7r8Zcas7FYkFL2Syn9Qoi3pJR2KWVnGM7dBoy3QHJC28aft09K+Thc42sE7/6WNBaLBYvFQtdwcO3ZZjYhgeRQ5j1AYXIySEl+StCB39o3gM/v5+6kEheKyBKheTMvmTLzHpBC0KE3cldvxgpsjjHS0tLyQstNNpuN+vp6ampqeO+993j33XfHkoXDQUJCAnv27MHlctHT08OlS5fCduylxqwUyzjWhPHcZ4FlQohCIYQR+ATw9vgBQojMcS9fA26E8fwLkjGLJVQjzGYyYjUYMI9rwRpnNpNgtZIeH4cmBG0DQ0gpudzWhkdFvkSDcM6becmMmfcEFcyg0cwvvf5zBAIB9u/fz40bN15IwcTGxrJ9+3YqKio4ceIEp06dwuv1voz4Y2iaxsaNG1m5ciVXr17l8uVF3/QzIjyvYglbHL6U0gf8Z2AfQYXxPSnlNSHEHwohXgsN+w0hxDUhxGXgN4BPhev8C53HWfcWk5EUm23CPiEEy9LSQEBhWjJSSm539uDx+bjWuWhvmuczizp/BWaZef/h95OTlMiyZcvYvXs3VquV/fv3c+nSpaeqEc+G+Ph4duzYwfLlyzl27BhnzpwJW8hwTk4OH//4x7l8+TInT54MyzGXErPJY9kSCg1GSrkntK02HCeXUv5MSlkqpSyWUn4ptO33pZRvh57/jpSyQkpZI6XcLqUMfzOHBYbRaESn1zMcWmO2GAykTlIsADkJCQghqMgOGn3XHnUQYzRyo6ODUZUMFnEiOW/ma/7X48z73KTECdtzkxInZN4/Jj8/n927d5ORkcGhQ4c4c+bMC1keiYmJNDQ0UFxczJEjRzh//vwLKarJmEwmfumXfon29nZ+8pOfqDyX52A2Fss+4JAQIm3ctiWfTxItnE4nw6OjYxFhNrOJOLP5qXExRiMpMTFkJsYTYzLSN2qnNZQRfbmt7anxirATkXkzn/O/Hmfet/YPTNje2j8wZeb9YzIyMnjllVcoKSnh+PHjvPvuuy/knE9OTmbnzp3k5uZy6NAhLl68+NIKRgjB66+/jqZpvP3223R3L3rjMyzMRrHcAv4UOCqE2BjatmhbrM53jEYjUtPGclisJiM2k2nKsSWpqfgDAWoL8gA4deceMUYj9/r6VF5L5InUvJmX+V/h6HmflJTEjh07qKys5NSpUxw9epTRFyhHlJaWxiuvvEJmZiaNjY1cuXLlpRIfhRC8733vw2KxcO3aNc6ePatyXZ7BbBSLlFK+Q9B5/lUhxH8G1LcaRdw+35jFMpNiyUlMRAhBdW4mZoOe9oEhbnf1oKGsljkgUvNmXuZ/hbPnfVxcHPX19dTV1XH58mUOHTrEwMDAM983mYyMDHbt2kVycjIHDhygqanphRWCpmns3LkTh8NBUlISe/fufSGlt1SYVT8WACnlbWALsBWojqRQiunxeDzYnU7sIT9JaqxtyjBPAJNeT1FKCl4p2bK8BIDD15sxaDoeDQ6qJmCRJZrzZs7zv6bKvI8L+Cb0vX/envdWq5VNmzaxefNmmpubOXjw4AstRWVnZ7N7927i4+NfKhrNYDDQ0NDAzZs32bx5M6dOneL27dvPfZylwGwSJFeNe26XUn4cKIqoVIppsVqt+DQNl9eHJgQZcXEzjl+enk4gEKAyJ5PMhDjsbg9n7j5AALci3GxpKRPBeTMv878mZ973aXqMUlLkc1PgcxEbUjIv0vPeaDSybt066uvrefToEfv376e1tfXZb5xEbm4uu3fvxmKxsG/fPpqbm59bwZjNZrZt28bx48fZvn07Xq+Xw4cPhy3cebHwvOHGAEgpH4ZbEMXscLlcPOrtA4I1whLHJUZORZLVSmZ8PE6vlx0rlgNw4X4r0h+gpbcXt4p0mTPCNG/mZf7X5Mz7gBD06gzcNZh5qDNhCimZEumlvb39hSwGvV5PbW0tO3fuZHh4mH379r1QsmVBQQG7d+9Gp9O90DFsNhsbNmygsbGR8vJyamtrOXDgQMS7Yi4kXkixKKKHpmn0hJIjY0wm4qeICJtMdXY23kCA9PhYlmem4Q8EONVyn0AgwANVQ2xBMV/zv6bLvIcnSuaBycqKmpUMDg5y6NAhjh07RkdHx3MrBk3TqKioYNeuXUgp2b9/P9evX38uB70QguLiYnbv3o3f72ffvn3cv39/1u9PTExk5cqVHDlyhLi4OPbs2cO9e/eUYz/EbPJYCmezTTE36PV6ekPJkTEm46xqL6XExJAdH8+ox8Om0mIAbrR34g9IbnR2qokQASI5b+Zj/tezMu8BDJrGx9fXsWLFChoaGli/fj39/f00NjZy7NgxOp/ztyiEoKSkhN27d2Oz2Th48OBzhxgLISgtLWX37t04nU727t0762W2tLQ0SkpKePfdd9E0jfXr15Odnc3evXsjUu5/ITEbi+WHU2z7QbgFUcwOl8tF30gwGmW2ikUIQV1+PoJg3ktJempQqbR1MOJy0atas0aCJTVvZpt5n52YMLbdaDRSUVHBzp07Wb9+Pb29vWNKpqur67mUTF5eHrt27SIrK2ss2fJ56ogJISgvL2f37t1jy2xts4iczM3NJS0tjXPnzgGQlZVFQ0MDZ86c4datW7M+/2JDP90OIUQZUAHECyE+Om5XHPDs9ZclhC8Q4EpbG0NOJ6tyc0mwWCJ2LqPRyIjnSaixdVyNsJmINZtZW1DAybt3WZmfw52uHpoedVCVm0VzV9eU2fuK52cpz5t1xYX8w2d+ke+fvcCBazdxeDxYjUZeqSjjY3W1E5TKZIxGI5WVlVRWVuJ2u2lububq1auYTCaWL19OWlratO8dT3p6Oq+88gr9/f2cOHECk8lEbW0tllnOSSEEFRUVlJeX09TUxLVr16iuriYjI2Pa9yxbtoympiaampqorKzEaDTS0NDAjRs3OHToEFu2bMEwy3m6WJhWsQDLCXa/SwA+OG77CPDZCMq04Dj/8CHNXV3odToO3rzJ+ysrsUToh+Tz+xkNdY6MtZgx6HSzfm9RcjL3+/rQCUGs2cyIy8Wgw8l90U9tXl7EZF5iLOl5k52YwJu7dvDmrh0vfAyTyURVVRVVVVW4XK4JSqasrIyUlJRnHuNxsuXIyAhnzpwBoLa2ltjY2FnJoGka1dXVVFZWcvnyZa5cucLKlSunVXCVlZWcO3eO27dvs2zZMgDKy8vJzs7mwIED1NbWzqicFhviWeamEGKDlPK9OZLnpVmzZo18bJbOBUNOJz9paiLWZEITgmGXi2WpqawtKIjI+dq6unjzB2/TM2rnjfW1fH7b1ud6/6DTyTtXr3L5wSPO3n1IdW42a5cVUJWVRU32VHl2S4KwV5JQ8yb8uFwubt26RU9PD2azmfLycpKTk2f1XqfTyYULF3C73axcuZKkpKTnOrff7+fSpUv09/dTW1s77XlPnjxJXl4eublPIsKllJw5cya4JF1XhxCLpnDJtB9kJovlMXeEEP8dKBg/Xkr5mZeXa+Fzp7cXAWihH4vNZOJ2Tw8rMjOnzYh/GTSjEUcoZj491NzreUiwWMhLSmLA7uTs3Yc0d3axtbyE6x0dLEtNfaFGTIopUfMmzJjNZmpqaoCgorh58yYXL17EYrFQXl4+o7KwWCxs2rQJj8fDpUuXGBoaorq6mvT09FmdW6fTsXr1anw+HxcvXuTixYusXr2axMSJBTc3btzIkSNHMBqNY8cWQrBu3Tra29vZu3cvmzdvnrXltFCZjWL5N+A4cBBQLQjHEZCSu729E/wcjxXM7Z4eVuXkhP2cA8MjY0lm2QnPr1gAytLTedDfT7Ithr5RO4/6BkmLt3Hq/n3qly0b+wyKl0LNmwhisVhYtSqYg+pwOLh58yYXLlzAarVSVlY2rZIxGo2sXbsWv9/PlStXuHz5MuXl5RMsjJnQ6/XU1dXh9Xo5f/48drudNWvWEB+6yRNCUF9fz4EDB6irq5ugeLKyskhJSeHEiRNkZWVRVlb2kt/C/GU2isUqpfytiEuyABlwOPD4fJgn5ZJYjUZudXVRlZmJ/jl8ILOh3+lESjAb9CQ9IzlyOlJtNqxGI6WZabx3+x432jspSqugbXCQd+/eZbXyt4QDNW/mCKvVSm1tsCOBw+Hgxo0bY0qmvLz8KasCghbIqlWrCAQC3Lhxg3379lFcXExxcfGslqoMBgPr16/H7XZz/vx53G43a9asITY2FiEEDQ0N7Nu3j61bt2IbFxhjNBrZsWMHN2/eXNSO/dkolneEEK9KKX8WcWkWGJ3Dw1NWFdRrGvZAgPahIfKecy33WXSNhRqbXvjiL4SgJDWVnuER3rt9j5auHrx+P/FmMw/6+7nf10es2UxeUhLFKSlPleX3BwLc6+ujpbd3rCPlqNuNTtMoTE6mPCMjIsuACww1b6KA1Wpl9epgBRu73c6NGzcYGhoiJiaGsrIyEhISJox/nGy5YsUKWlpa2L9/Pzk5OZSXl6M9Iy8HgoEGGzduxOVycf78ebxeL3V1dcTExLBz504OHDhAQ0PDUzefZWVlZGdnc/DgQVauXElmZuY0Z1iYTOu8F0KMEKzGKoAYwA14Q6+llHLmIlVRYi6dkHuvXw9GVtmdnGxuYdTtYVlGKuuKC3D7fCRYrewuLw/rOf9y3wF+eOkaecmJ/PnPvz6h1/1MeLvseLscGNKtGNJj6B0dZd+NG/z0YhPtA0O8r2YFK0JNwaSUeP1+XD4fEqjOyqIiMxOdpuH0ejnR0kLn8DBmvX5s2Uyv0xEI1YISQE1ODmXp6dNmY88zwrb2p+bN/GR0dHRMydhsNsrLy8eWrybT2trKjRs3SE1Npbq6Gt1zrDo4nU7OnTuHlJI1a9YghODQoUPs2rVrSstESjmWrb927dqF5th/fue9lHJBe5eklHSNjPCgvx9/IEBBcjJZL+DsBgg4vAQcPnQJJoQ+eKH0+v302e243F6+e+o8vlA5id6RUbqGRnhtVSU9o6MMOZ3EhzGvpd8ZbIBkNRlnbbF4Oux4Hg4Hnz8YQZdoJikmBr1OR2lmOu0DQ9xs7xpTLEIIjHo9Rr0efyDA5bY2bnR2Em+xMOBwIIF4s/mpSaAJQZzZjC8Q4EJrK3d7e1mdl0dmXNy8mzABpw9vtwNTfniv8wt93ixWbDYbdXV1AIyMjHDjxg2Gh4eJjY2lvLycuHHFXHNzc8nNzaWrq4vDhw8TGxvLqlWrMM4isMVisbBlyxbsdjtnz55Fp9NRV1dHY2Mju3btesoKEkKwdu1aOjs72bt3L5s2bZogy0LlmUth07RTHQIehOoWzSsCUvKwv59rHR302u3EDks8MRotvb2szc9n+SyjQB7jt3txXetDSolm0WOpTEFoIniBlZKD127iCwRYnplGWVYG+65c5253L8eb71JbmMvd3l5WzdIxOBseN+iKMZkwz0KxBJw+hg8+wNfjRBg1jDmx6DOsmPLiyE1MxO5yownB/d5+RlwuYieZ7DpNI8FiwRcIMOp2YzEYnmmF6EPvcXq9NN66RZrNxoaioik7XUYLT9sovj4n+hQLupjwr3EvtHmzlIiNjWXt2rUADA8Pc+PGDUZGRoiLi6O8vHwsYis9PZ309HQGBgaeO9kyJiaGbdu2MTo6yrlz58YKXu7Zs2fKm6yMjAx27tzJ8ePHyczMpDzMKx1zzWx8LH8N1AJXQ6+rgCaCmcWfl1Luj5RwL8Kgw8HxlhZGnC4utbTS0T+EUaejKCsVnz9AamwsSVbrrI/naR0ZKy0RcPrw9TkxpFrptdtpHxymY3AYs8HArqpyjHo95tXVfP/0Rc7fe0h+ShK3tG4qs7KeK5FxOnx+P8MhiyUhxvJU9Ja3007A5cOYG4fQBff1/6AZ17W+J2MejeLrcZDyK1XkJiRwt7eXkvRUmju7OXf3IdtXlE55br2moX/OZS2LwYBZr6ff4WDv9evsKi+PaFWC2SL9EvvZTnw9jqBiKU6IxGkW1LxZqsTFxbFu3TogqGSuX7/OyMgI8fHxlJeXY7PZSExMfOFkS5vNRn19PcPDwxw5coRvf/vbfOxjH8M0hQ/SYDCwY8cObt26RWNjI1u2bJmVlTQfmc2Voh1YJaVcLaVcDawE7gKvAH8SQdleCE0Irj5s43unLmB75OFPesv4j705XLnbyo/OXGLftetP1SCS3qmrogYcXnyDLqT/yX5fX/DC3jE0xO2OYJnsVQU5GPVBHZ2TlMj6kmCtwcamm7i8XtoGB2ctv/QFcN8fwttpf0pOh9eL0xuMXJ1cgkX6AnhaR/B2OXA19yP9EkdTb1CpaGBdnYa5Ihk0gbtliNF328aOsa44H4BLDx4x+GiirMInJzRrel6EENhMJqSUHLl9G+9L9iAPB66WAZyXe/C22xn4bsTqOS2oeaN4omR27txJaWkpTU1NHDx4kLNnzzI6OkpsbCzbtm1j7dq1XLlyhcbGRvpnWR08Li6O1157jZUrV/Ld736X06dPT9vDZfny5axdu5bGxkba29vD+RHnjNlYLKVSymuPX0gprwshyqSUd+fbujlAv93BqTv30UvBrzkLMUuNMl8sn3UV8hXRwreOn6YyK4uq7CwAXPeH8D4cJmZ9FppxolUxeqaTkaOPkE4fukQT5uXBCC+/10/n8DAPeoPtUsuzJpZqWFucz+3ObnpGRmlqbSc9Lo6CWWYIu+8N4esPKi//iAdTcQJCC37PDo9nrHNkRvzEdVhfjxP7hS78A270KRY8j0YYfbcDAEt1KqbiBKQ3gDBoOC/1MNzYSuaaDJJjYhjV6ajNzOJCRzvfu3aFXb5SMrNTsHk1Ytu9+E2CoWwDaAIpJXe7e7ne3onH58eo02HU65CAy+slMcZKSVoqWYnxE0x+q9HIkMvFlfZ2VodxaRCCSlUG5Nj/b9jlYsDhIM5sJnEK69R+pnPsuX/4+RtPzZIFNW8UE4mPj2f9+vUADA4OcvXqVex2O4mJiZSVlT2VbFlVVTWrki2VlZXo9Xp6eno4evQoCQkJrFy5Er1+4qU4NjaWPXv2cPbsWR4+fMjatWtnFaU2X5iNYrkmhPgb4Duh128A14UQJoLRLi+MEGIP8L8BHfA1KeWXJ+03Af9MsANeH/CGlPL+TMf0+v3EW8z8QlIp5iYNGZrDa10JbI5J44S3mz/+8V7+4bO/hOhx0f/P1wk4fPgGPCR+qHjsOPYL3QwffAC+4N26f8CN/VQHxvxYPNlmWvsG8Ph8pNhiSIyZePHSaRo7Kkr57qkLXLz/iOL0VAadTmxujYEfNqOLNZL4kWUIw8Qfit/uxdNhJzDiQbMZ8PW7CLj6MGbb0CeZsbvd2MeSIxMmvHf40EM894IOev+gG/edQQD0aVYSXy9FiNCyHhLPwxH8/S6GGx9SsDqJ862tfCCxAFe/k2GHk3+9dg2uQ6rRSp45ljxzLAnuGG6LoKIccMzc0vjc3YckWC1U5WZRkZNJTMjsjzUaudHZSVFy8pQX/BfF1TJIwO6FigSudnRwZ1z/9pKUFNbk548t4wXcftzNwRsCS00K5rLZKfwXIGLzRjG3JCQksGHDBgAGBga4fPkyDoeDpKQkKisrMZlMXL16dSzZMi8vb8bjlZWV4Xa7MRqNpKSkcPjwYZKTk6mpqZkQgTbZsb958+YF49ifjWL5FPCfgDdDr08C/w/BybH9RU8shNABf0VwaeARcFYI8baU8vq4Yb8CDEgpS4QQnwD+mOAEnRarychH165ixbkA4MO6Mg2/w4vn1gCfsudyI2aIR6PDfPd7R9h5P5aAI+hHtZ9ux7YxE0OqFf+Ih6Gf3gWfxJBjw7oyDfeDYVzX+/A8GMF1qJUHWtAELsmYuld4TlIiZVnp3Gzv4vSde2zIziP9B0GLAgBNkPjhEgIOH8KoIfQa9jOdQWXmlyDAWBiPZUUSjn4nxrxY2kf68QckRr2OlNgnS2HSF8B1KyiPpSYVX58TX5cDfaqF5E9VooUi2UwF8RiybOCHobdbsJ/tJGNrDZbRAHH98NvdRWgB6NV5OGjsZsTtwzAcYK/xHo7uJ0tYcWYTNfk5pMTa8Pr9uL3B79Bk0NMxOMyt9i4GHU6O32rhxK0WcpMTqSvKpyA1GZ0QnLl/n1fKyyf4iKRf4u2yo5n16JOmd/J3j4zQ0tvLsMuFPxBA3++h9KAdnUdyr0lHa7lhrG6blJLbPT14/H42FxejCYH9TAfSE0CXYCJ2R15EHPchPkUE5o0iuiQmJrJx40YA+vv7uXjxIk6nk+TkZDZv3syDBw/Yt28fRUVFlJSUTBsNWVNTw5kzZzCZTLzyyit0dXXR2NhIeno6VVVVE6yTjIwMXnnlFU6cOEF6ejorVqyYk8/6MjxTsUgpncD/F3pMZvQlzr0WuCOlvAsghPgO8CFgvGL5EPDF0PMfAF8VQgg5Q+XMGJOJFIuVxO7g3XvstmyEEPS0jWIe9fJHgQr2a51UX9cIBLzoEk1oFj3edjsjh1tJ+vhyBn92j4DdixZrxLo6HaEJzCUJ6OKM2E+2I64PYU4MLldVmpLwDLlwhnxxBrckflAiU0xsLSuhpauXru5BPD+8h39AhzDqkB4/jvNd6JPMaKELm3/IzejxNvBLhFWPdPjw3B3Cc3coeGAB5AR9PXH6iaHGzmt9SE8ALc6IeVkCLE9EM+sxL0t8yirSjDpsGzKxv9cejBQ70U2qU5B814MWACkgxW/kE84n5Wg+6M3kLxLvEWs2UxObQmlMEn6rjpG0py/KyzPT2bq8mHvdfVy/28adoX4e9g3wsG+AjcsKWV9SSPfoKFfb26nMzKTXbufeww7M54dIfOhFbzOS0JCPpTQJzfTk7s3j93Ph4UNu9/Sg0zQMmhZs0nTJid4d/DkU3vCjizPitvmRGgSMAptTz0BfL00mM5WpGdjPBf1i1to0jOkxyEBkmpxFcN4o5glJSUls2rQJgL6+Pi5duoTT6SQtLQ2fz8f+/fvJzs5mxYoVUy5jrV27lmPHjmE2m8nKymLXrl10dHRw4MABsrKyqKioGHufwWBg+/btNDc3LwjH/kz9WL4npfy4EOIqPJ1gLqWsfslzZwPjW7U9AtZNN0ZK6RNCDAHJQO+0cgOb9OmMegbRJZoxZtiQUhJbn8PIkUdYR718mGC+RpfRQ/YvV5Dg0tHzd1dwXOnFWBCH80pwKSX+fQUYc2PxPBhB+gMY0qyYyhJx3xzg04O5dCV5yLYbMLglZTWFxCTFMHCjmy7vIIN3HMRl6vlwUhFr7xjJDOiQBkH6r9cwcvQRjgvdOJt6iVmXiad9FMe5LvBL9OlWYtZnEhjx4Ljcg7/PFUqtg8JWjW+xmr9PbcUsnlx0HReCF0tzSQLWVeljEWHTfkdCYNuaw+APb+O41E1Ogg6DW+KxCLoqTJiHAlj7/CDANBIg2WPgd91ldBcYIXQHphsNoHcG8FmenjCaprHaFcd2acVv8XPPM8z33Pd59/Y9dJrG6sI8rrS10dT6iKTOAIU3/NiGgj8xOeqk68fNaBtSSaxMh3gjHcPDNLW34/J6iTeb6RkZ5XxbJ6O9I9S2FQAanRYvGU4DuWecjGTocSbo8MQEZY0dEfib79Nm6EJ0OUATxGwM+tge+6/CxRzMG8U8JDk5eUzJ9Pb2cvPmTXQ6HT09Pezdu5eMjAyqq6uf8qVs2bKFQ4cOjS2LZWZmkpmZyaNHjzhw4AC5ubmUl5ePWT6lpaVkZWXR2NhIdXU12fO0IvlMFstvhv5+YC4EeRmEEJ8DPgeQk5tLoCm0hl6Z/Hg/pqIE0Gl420bp7OznX/qbOW0c4LWbOhrKy0jKjMHXYWfwX+8AYCpJIKY2mPOisxlxNfcTcPnRSuLpuNdFptvIn/SVE+j3oZklOuMw/jQfxnsOUq+NkCpBNvlZSXBN9IHOwdmqAG+mxxC/uwDHlWBU0ujJdnzdwdwUQ24sMbVpoAl08SZit+Yg/QE0mxH3gyEcJ4PO+M/25GLs90IsBDx+XC2DANi2ZD9TqTwmZlUaQz+7R2DIg3EoeAXsLzAgdYKRTD2OZB0iADqPJKPJhXk4QMIjH4O5T6yUmF4/QzliTNkAEJBkXnaRdsuNFlo9S8PGWir4jqWNn95qwWI0UmdLI/Wmm6QHPjQ/BDQYyjZg6/ZhtAdwnOriQncvjngNtwlsAR3Jfo3G/hZOPHoAwEedmRjRuGgY5C9MLfxHfyEbPUnEt/uIb/cxagpg1HQYnY+v78HvWVeeQL90Y3VLYozGcCdvLph5o4gMKSkpbN68GYDu7m5u3brFw4cPuXnzJvn5+dTV1Y1ZG0IItm/fzv79+9mwYcNYNYCcnBxycnLGltYKCwspLS0N3hTabOzZs4dz587x8OFD1q1bN+8c+zNl3neE/j6YvE8IcRLY9JLnbgPGhwflhLZNNeaREEIPxBN04k+W9S3gLYA1q9dI57WgQROz5kkypCHVird9FGNuLLm5sTiuP8DbLTl95z7xVgvFxQbyO0KD9YKEcY58zaLHsiIFT+swfT2D/G3sff6bp4gYqUeTgNOP/d129GlWfF3Bi5cwaOANENDgdoqHP/bchFaNj/T1k5+cROyWHEYOt44pFXNFMjHrMjAXJiClxNfnwtdtB6nDXJaEpySGo9eus20wWFDPfrAV8y/acFzqAZ9En2bBmDn7LpBCrxGzOp3RE8Gv3JFhYDBO4s43EjAICEgMTklsl4/eEiNpNz3Edfhwx2g4koLWkt4VwNbtw5GkR+8KkHTXQ8odDwZX8ELujNPwWjV0Hom138/PO3Po0dxcuNTCB3wmkhzB4zgSNLrLTAwUGrF1+Sg8bsc6Cvm3A7jjwDwYwOAO+ruXkcJ6g4HThW4+0BL0bw2ssJLjSOKve+5xzNjLam8C6zyJxLkNgMSlC+DMMAYDFzTJw+xRXDdvEgAsej0/F6qSGw7mYN4oFhBpaWmkpaUhpaSnp4dLly7xzW9+k/j4eHbu3ElCQgKaprFz5072799PfX091nGBLfn5+eTl5XHv3j327dtHSUnJWKHMuro6urq6xjL2pytREw1m47yfipnDHmbHWWCZEKKQoAL5BPALk8a8DXwSeA94HTg0k38FIDDqCTpnk80Y0p/U0RIGDWNeHJ4Hw+iTzPzC9nVc++E7POofpHtwBEN6PDE5OlIGBbHrszCkToxaEgYNU1ECZwfv0eIf4YsJt/hPKVWUFWXDvRHczYNjSsWyMhVzUQIiRo+WaeF+fwcphx7SPjDEVw8d4cs/92Fit+cS8Prx97nQZ8RgrUrBmBVUDAIwZsZgyLCCDC7X2B0ezpuH2UZQsbhvD9Lz1lWkO2gWWGtm17p1PHE78/ANBn1F+jwrF3XdmHUSHQI0gTdGMJqqQ/glg3kGEh96SW7xYOoTeI0wnKShuTXym1zEdfmDShbwGqGv0IA9PWgBISC5xUPiAy//xV6EF4kRDYfm52TqCM0pHmINMZT60yHTQPtKC5lXXRidEqMz+Pm8IoAdP3FSz0pvPCubg+eyxwq0IhsfMlXjcHto6xvgfredG45OLMOSPp+Lm/pRimJT+EBpOfFDYEjRYwg59/2BqXOYIkQ45o1iASKEIC0tjV27diGl5P79+7z99tv4fD6qqqqoqamhoaGBgwcPsnPnzgkJlEIIioqKKCwspKWlhX379rF8+XIKCwvHWjGfPHmS1NRUKioqovgpn/CiiuWlPZ4hn8l/BvYRDDf+upTymhDiD4FzUsq3gX8AvimEuAP0E1Q+M+IfCYbjxjU8PYcNaVYMaUGFUeGNYU1RHieb77K/6QYfX1dLa5WJvi4vGSU6YqV8KrPd4/fz7u27ACRkxGPMjyehLBV/Tjz6RDOeDjuGjBisFckYC+PH8io2JRTRPTTMVw8e4/y9Vn5y5SqvVlViLkrAn+pBs+rHlMp4hBBjZd4GHQ5uixH+1HabjyYXUtJhxNsW8gHrBTHrnr/tqWbWE7M6Hf+wB2OWjVqjhTMPHqCFaoXpNQ1njMBgBldSAM2hEd8bIG4g+O9P6vaD8CNC12Zfphl3ZRzt6QE6XXYC8slFe7BSR9loAFufHyOCS/oh/jbmHqMefzCVsB0O37jNiuwM1mRlQ4UJW7ePAZeTd/xtnA/0odfr+IXlVdSf1qEPFUXJ2F3C5kIDJ+7exWIysiI3C3LBYA9aUzeHernddZvmrh6+5XTxgVWVJArD2Pcb7rYGzyAykQKKBYUQgsLCQgoLC3E6nRw7dox/+Zd/ITk5eayV8Z49e57yxwghxiyW5uZm9u7dS3l5Ofn5+dTX13P79m0OHjzIli1bxhSTr8/JyPE2HBe7kW4/wqTDuiqN2C3Z6JMjVwVjpurGH53uPcDfSimnjrONMrWVK+WhP/lXEl4teubYY3fu8E/H3+Nudx8GnY4P1lZRFBPPgPCSERfHluLiCfW4bnZ28lvf/zeGHE721Kzg05s2kB4q6+BqGcTX60Qz64L1xHSTclQCAd78zg+42tpOYWoyP7d2FfWp+XB/FPPyRPQJM9fROnjrFl/6t39HSviz8i0ss6Yw+m470ukjdlsu8XsKnv/LAnxDbnyddkzLEhGaoGd0lHt9ffTZ7Yy6XEjA5BMU9hrIjo9HN+rHP+jG1+fE2zoCEvSpFuJfLcRS/iQnJCAlLq8XXyDAkNPJhdZW/N0OilsgIOBsjpN2nwN/IBCs79bbT2v/4Nj7c2Lj8Lh9dHuCVmCKLYb3razAbDRQFIih+LIPY0YM8R8oQmiC1oEBjty+jc1kGstZEX6JCECPy86Pzl9hKJR/Y9LrEUKQYLVQlJbM//zoh8JZ3XhBzpvFXt14vuP1erl48SJtbW1IKenr62Pjxo2UlpZO269FSsmNGzd4+PAhlZWV5OTkYLfbOX78OJWVlSTbrfR98/pYLt4E9ILkX1qBZflLtfWYdt7MpFj+caYjSik//TISRYo1a9bIs6fOjFUhnok+u52fXL3Ke833uNXRhRCCjcsKqSvMw+71YjOZqF+2jFizmc7hYf7l9FneudiExWjgV7dv4qM1NWOOXxmQeNtG0SWZp82NaO7u5vPf+A7+QIDX1lSzIjOD7TlFmOKeXZzx70+c5NsnzxJrNPInm3ZStroEb7sdT8coMavSnworDjeetlE8j0YmbAu4fEhvgJhVac+8+/H4/Ry7fZvBtkGssWb8lqcthd6RUS7ef8T19g58oTI6FqOBtUX5rMzPxev3YdTrebWyErqcaCbdhPPe6uri9IMHxI5TLo9xeb0cvt7MzfYuApN+80d+681wKpYFO2+UYok+fr+fpqYm7t27x+joKNnZ2fj9frKzsykpKZlSyQQCAa5du0Z7ezvV1dVkZGRw9tgp+g7epcKbizbd9V8nyPivq1/GcnmhsvnzcgLMhtkoFYDkmBjyk5IwlOtIjLFw6s59Tjbf5V5PH6/WVOD0ePjJ1asIIQgEAlx+GHR0l2VlUJOdPSGaSGgCY+7MRelK09LYUFLIieYWzt15QIothqahHlbHzVzixOP30zkUzMtJsJjRZZgROg1jbuwzzxkuDJkx+HocBNxPEiU1sx5dumlWP0yjTkd9aSlnTQ+409OD0aPDqNejE2Lse0yJtfFKVRlbyopp7RvEbNCTlRiPTtPw+f14AgF2FBdj1OlgiqXD0rSgn+lcaytSyrHCmUadDrPBwPtqKthVVY7H5ycgA3QMDnOve9rI9RdiIc8bRfTR6XTU1NRQXV3NyZMnaW5uZtOmTVitVk6ePInP5yM3N5fi4uKxpTJN06iqqqKiooIrV65w9epV8jvisPiSOG24TaUvl1g5xRz1S4YOPCD5E+FvkTxTHst/nemNUsr/FXZposDqvDw6mppYW1xATlIie69cp31giH8+cZqdlWWUZQYjyy4/bKO1bwCTXs/K/GzyX7Az5Kc3r+fcvQe0Dw7RPTiMpmkUJCXN2LBr2OkcW8aJt1lJiEJZB6EJjHlxuG4PPNmm1zAVzF4WvaaxvqCAnPh4bvX0MOBwYPf5QEpsJtNYOX6zwcCycRUNfH4/ox4PGwoLSbFNH/kmhGB5ejrZCQk87O+nz+Fg1OWi/3GrAaMRnaZhMQbPU5KeSkl6eFemlsq8UUQWIQSbN28mIyOD27dvo2ka2dnZlJeX097ezokTJ/D7/eTm5lJUVIRer0fTNFauXInf7+fIF37AkM5BqS+TO7pO4qWVIv/TLUOcV3pgLhUL8PhWeDlQRzBCC+CDwJmwSxIl4sxmNhcVcezOHTIT4vjlzevYf/UGd7p6+Nmla8Glk0CA+73BkinrSgrYUFj4wmXwC1NSWL+siCPXmznRfJePrVvFuYcP2VVWNm0+xbDLxYA9eHGMN5swRSlmXZ9kxpgXh7d9FF28CWOuDc30fPEfQghyk5LIDSlmu8fD3Z4errS3o9c0rJOyiZ1eL16/n01FRRSlpMzqHDaTiRXjWr06vV5udHZyo7NzwjhBUNmEmSUxbxRzQ0lJCW63G7/fT0JCAo2NjWN5MjqdjtbWVo4fP04gECAvL4+ioiJ0Oh3LfVn4CdCs68CHH6/00ydGSJ7chy5CQZEzLYX9AYAQ4hhQK6UcCb3+IvDTyIgTHfKSkqhftozjLS0EpOSDqyq5+qidI9dvcze0VKLXNDaUFrK2KJ/i1Be/y9WE4PXVK7l4v5VBh5PbnT0ECNbASp/GEukZHR2zWFJs1pcqY/+yGDNjMGbOrh3ybIgxGqnKziY7MZFjd+4w5HQSYzKhE4IRtxuzwUDD8uWzbsE8FRaDgdrcXMrS0+l3ONCEICAlrQMD3OsN+1LYkpk3irmhoqKCCxcu4HA42LVrF93d3Rw5cgSbzcaqVavIy8sjEAjQ2trK0aNHkVISo42QFUik3J+NDz/Nug6EDJsr8ZnM5nYzHRhfW9wT2raoyElM5INVVZxvbeVhfz9lWRkUpCRzt7sXg05HXkoiQtPYVFz80n3ci1JSWFucz4GrNzl15z6FaSlcaW9nZ2zslFZL1/Awg/agYilMT5uQQLVYSLJaeX9FBddDloUvECAvMZG6/PxZt2B+FlajcYJFlJOQwMqcnBne8VIsiXmjmBtqa2t57733ePDgAfn5+ezcuZPBwUHeffdd9Ho9q1evJj8/n/z8fAKBAOcPvcNF/X0kkoxAAmX+LLSp2m9FSNfMRrH8M3BGCPGj0OsPA/8UGXGii81kYmtxMR2pqZy4cwdNp1GTn4PT48EvJVuKi8NS7t1qNLK+qICm1nY6Bodp6exBpwn6HY6n7sx9fj+PBgbxBQJYjQaSDAbsdjsJk8rmLwYMOh012dlUZWXhDwTC0nXzWYRLaU3Bkpk3irlh/fr1HD16FJPJREZGBgkJCWzfvn2s/XEgEKC2tpa4uDgyZSKZ/kQCBOjUBrmgv0eRP50kOclHGS3FIqX8khDi34EtoU2fllJejIw40UcIQVZ8PK9WVnLm/n06R0ZIs9mozc0l6SWWYyazLD2dyrxsOgaHufzwEaXZ6TS1t7Nt2bIJ4wadTvpG7QAk2WJIioub11VNw4EmBNrcJi6GnaU2bxSRRwjBtm3bOHjwIEajkaSQn9Jms7F161ZcLhcXLlzA6XSSIX3EY0VDIyuQRFZgmmCjufaxjEdKeQG4EBkR5ic2k4kdy5dH7PiZcXHkpyQRZzEz5HTROxTMERlwOCZYRX12O70jQcWSEmsba5qlmP8sxXmjiCxCCHbs2MH+/fvZvHkzsbFPnPFms5mNGzfi9Xo5cviH3KKdIn86KZMd9uOJUBzQ/CqJuYQw6HQUpaSwIjtYiuXig0dowOW2iXU42wYHGQhZLMmxMRgBjydi7XQVCsU8R6fTsXPnTo4dO4bT+XQ3V4PBQLnMYY2viAFtlNP623RoA8ipKgpFqE22UixRpCg5mdLMdPSaxoPefrw+P60DA3SPBK0Xn99P58gIvSHFkmKLITEuDrP52Zn6CoVi8WIwGGhoaODQoUN4vVN0uvZLNDSW+TNZ6yvBjY8BYZ9yXCRQiiWKpNpsxJiMLM96koRp0uk4ff8+/kCAjuFhRl1uHG4PRr2O7MQEAn4/Pp8vypIrFIpoYzab2bp1KwcPHsTv90/YJ8Z1XxUICgKpTzvuJ40LJ0qxRBGdplGUkkJpZrAUybVHHeg1HYNOJ2fu3+dCayv9Y9aKjRSbjUAgQGBuS70rFIp5SmxsLOvWraOxsZHxdR+tq9KefXXXQuMigFIsUSY/KYmUWBtZifG4fT5utHcSZzZzt68Ph8dDz3CwNH5aQiwpMTGYzWa1FKZQKMZISkqiurp6LDkSIHZL9lMV1icjdBqxWyLT2lgpliiTYrNh0Omoygn2YL/04BGCYKkZm8lEW6iUfGZ8HHEWCw6HA0eo9pVCoVAAZGRkUFRUxHvvvQeAPtlC0i+WB6ueT77Ka8HGhUm/WB6xnixKsUQZTQgKk5PJTk7AajTSMzJKa3+w0KPb66N7eAQBpMbHEmc2YzQaF30ei0KheH7y8vJISUnhwoVghLtleRLpb9YSszYz6EsRQZ9KzNpM0t+sfdleLDPyoh0kFWGkIDmZ5u5uVubn8O7tu5xsvkvu+kTu9fQSkJLsxHhsJhMWgwGnctwrFIppKC0t5erVq1y7do2Kigr0yRYSP1xC4odL5lQOpVjmAckxMZj0eqrzsrj4oJX2gSHudPVwo70LgPyUZFJsNoQQYzksi7FemEKheHmqqqo4e/YsLS0t5CdkRaU1sVoKmwdoQlCSmoo/EGB9SSEAb1+4yt3uXnSaRmFaMhmhDFuLxYLFErkfhEKhWPjU1dXx8EoLV//iMPYzHchQgz7p9mM/00HXX1zAeas/YudXimWeUJCcjBSClXnZFKY+6R2/tjgfq8lIcqjBlcfjUZn3CoViRnx9TpZdsPJAdjMgJyVGBkB6A/R/6wa+vqcz98OBWgqbJ8SbzSRaLDg8Hj60upqb7V0Y9TqK01IY8XjCUlVZoVAsDUaOt4FfUhso4qz+Div8OU+1J5b+ACPH2yLif1EWyzxBCEF5RgZuvx+dplGRk8myjDQ8fj/JMTHBPu+gosIUCsUzcVzshgBoCNb4imnSt+Jk0kpHIDQuAijFMo/ITUxEr2n4xmXWu30+CpKehAU6nc4pC88pFArFYx77VAB0aKzxFuEXT1fskB7/U9vCQVQUixAiSQhxQAhxO/Q3cZpxfiHEpdDj7anGLCYMOh0rMjKwu90AwSxaIcge19RLWSwKheJZTK4BZkCPTT5dsUMYF1etsN8GGqWUy4DG0OupcEopV4Yer82deNFjeXo6Bp0Oj9+Pw+slzWYjTpVwUSgUz8FSrRX2IeAboeffINi2VQGY9Ho2Fhbi8nqRUrI6L2/CfhUVplAonkW0a4VFKyosXUrZEXreCaRPM84shDgH+IAvSyl/PNUgIcTngM9BsKzBQic3KYn3WyxoQjxlrajESIVC8Swe1wrr/9YNpD8wsQWxFlQqkawVFjHFIoQ4CGRMset3x7+QUkohxHTdZvKllG1CiCLgkBDiqpSyZfIgKeVbwFsAa9asiUznmjkmYZokSJfLBaD8LAqFYkYe1woby7z3+BHGucm8j5hikVLunG6fEKJLCJEppewQQmQCU8a8SSnbQn/vCiGOAKuApxTLUkLTVCCfQqGYHdGqFRatq9TbwCdDzz8J/NvkAUKIRCGEKfQ8BdgEXJ8zCecper0evV7ltSoUivlLtBTLl4FXhBC3gZ2h1wgh1gghvhYaUw6cE0JcBg4T9LEsecXicrnGlsMUCoViPhKVW18pZR/QMMX2c8Cvhp6/C1TNsWjzHuVbUSgU8x21YK9QKBSKsKIUywJD5bEoFIr5jvICLzBiYmKiLYJCoVDMiLJYFhgOhwOHwxFtMRQKhWJalMWywNDpIlM0TqFQKMKFUiwLDJXDolAo5jtqKWyBofJYli6q3YRioaAUywJD9WNZ0qh2E4oFgVIsCsXCQbWbUCwIlGJZYKg8liXNc7WbEEKcEkJ8eG5EUyieIKRcFFXmxxBCjAC3oi1HhEkBeqMtRASZ68/XK6XcM4fnm5ZntJv4hpQyYdzYASnlU34WIUT2+HYTQMNU7SbG9zEClhP+ebPYf6cLjXD/P6adN4tRsZyTUq6JthyRZLF/xsX++V4UIcQtoH5cu4kjUsrlz3jPPwHvSCl/MBcyTjq3+j/OI+by/6GWwhSKhYNqN6FYECjFolAsHFS7CcWCYDFm270VbQHmgMX+GRf753shFmC7CfV/nF/M2f9j0flYFAqFQhFd1FKYQqFQKMLKolQsQog/FULcFEJcEUL8SAiREG2ZwoEQYo8Q4pYQ4o4QYrqs6wWLECJXCHFYCHFdCHFNCPGb0ZZpKSOEyBBCfEcI0SKEOC+E+JkQ4nNCiHee8zhfE0KsiJScSw0hxO+G5seVUNmedRE81xeFEP/P875vUSoW4ABQKaWsBpqB34myPC+NEEIH/BXwPmAF8POLcLL6gP9bSrkCWA/8+iL8jAsCIYQAfkQwpLlYSrma4DyaLilzWqSUv6oCCMKDEGID8AGgNnR92wm0Rleqp1mUikVKuV9K6Qu9PAXkRFOeMLEWuCOlvCul9ADfIVjiY9EgpeyQUl4IPR8BbgDZ0ZVqybId8Eop//bxBinlZeA4YBNC/CC0KvDtkBJCCNEghLgohLgqhPj6uLDnI0KINaHne4QQF4QQl4UQjaFtMaHxZ0LvX1S/6zCTSTAx0Q0gpeyVUrbP8N3fD4WdP44ePBJ6/sXQuCNCiLtCiN94fIKQRdQshDhBMHH2uVmUimUSnwH+PdpChIFsJt6ZPGIRX3SFEAXAKuB0lEVZqlQC56fZtwp4k6DlXARsEkKYgX8C3pBSVhGMOP38+DcJIVKBvwd+TkpZA3wstOt3gUNSyrUEFdqfCiFUq9Sp2Q/khi78fy2E2Dab734ayoDdBG9avyCEMAghVgOfAFYCrwJ1LyLkglUsQoiDQoimKR4fGjfmdwkur3w7epIqnhchhA34IfCmlHI42vIonuKMlPKRlDIAXAIKCN7Z3pNSNofGfAPYOul964FjUsp7AFLK/tD2XcBvCyEuAUcAM5AXQfkXLFLKUWA1wVI8PcB3gf/Is7/7qfiplNItpewFugkuc24BfiSldITm3gu1XViweSxSyp0z7RdCfIrgWmSDXBwx1W1A7rjXOaFtiwohhIGgUvm2lPJfoy3PEuYa8Po0+9zjnvt5+euIIGjFLPYaf2FBSuknqICPCCGuAr8+w3AfTwwI86R94f4/jrFgLZaZEELsAf4b8JqUcrE0iD8LLBNCFAohjATN1UXVxCm0Vv8PwA0p5f+KtjxLnEOAKVSoEgAhRDXBO9qpuAUUCCFKQq9/CTg6acwpYKsQojB0vKTQ9n3Afxnnq1kVno+w+BBCLBdCLBu3aSXQwvTf/X2CFg7Az83iFMeADwshLEKIWOCDLyLnolQswFeBWOBAKBzvb5/1hvlOKBjhPxOchDeA70kpr0VXqrCzieCk2CGedEB8NdpCLUVCVv5HgJ2hcONrwB8RLNc/1XgX8Gng+6G76ADwt5PG9BBcwvnXUMmZ74Z2/Q/AAFwJned/ROAjLRZswDdCIflXCPq5fpvpv/s/AP63EOIcQatkRkLBM98FLhP0TZ99ESFV5r1CoVAowspitVgUCoVCESWUYlEoFApFWFGKRaFQKBRhRSkWhUKhUIQVpVgUCoVCEVaUYlEoFApFWFGKZZ4ihEgel8vRKYRoCz0fFUL8dYTO+aYQ4pdn2P8BIcQfRuLcCoVi8aDyWBYAQogvAqNSyj+L4Dn0wAWC5bh904wRoTGbFlFFA4VCEWaUxbLAEELUi1CjpVDp628IIY4LIR4IIT4qhPiTUOnsvaG6WwghVgshjopgs6Z9QojMKQ69A7jwWKkIIX7jcXavEOI7MJaNfYRgDTaFQqGYEqVYFj7FBJXCa8C3gMOh0tlO4P0h5fKXwOuhZk1fB740xXE2MbFM+m8Dq0LNhH5t3PZzTF8vSqFQKBZudWPFGP8upfSGagTpgL2h7Vd5Us68kmDdNEJjOqY4TibBGmSPuQJ8WwjxY+DH47Z3A1nhE1+hUCw2lGJZ+DzuJBcQQnjHtQgIEPz/CuCalHLDM47jZGJZ7fcT7OnwQeB3hRBVoWUyc2isQqFQTIlaClv83AJSRbBXNqEucRVTjLsBlITGaECulPIw8FtAPMGqqgClQFPEpVYoFAsWpVgWOVJKD8GGTX8cKlV+Cdg4xdB/50nXOR3wrdDy2kXgK1LKwdC+7cBPIymzQqFY2KhwY8UYQogfAf9NSnl7mv3pwL9IKRvmVjKFQrGQUIpFMYYQYjmQLqU8Ns3+OsArpbw0p4IpFIoFhVIsCoVCoQgryseiUCgUirCiFItCoVAowopSLAqFQqEIK0qxKBQKhSKsKMWiUCgUirDy/wNxIm7a8tjZCwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# load mouse data from zenodo and run basic analysis of the CoT session\n", + "\n", + "x_range = [-2, 3]\n", + "y_range = [-0.75, 1.5]\n", + "nr_mice = len(info['AnimalID'].unique())\n", + "for m, mouse in enumerate(info['AnimalID']):\n", + " date = str(info[info['AnimalID']==mouse]['Date'].values[0])\n", + " fiber_side = info[info['AnimalID']==mouse]['fiber_side'].values[0]\n", + " protocol = info[info['AnimalID']==mouse]['protocol1'].values[0]\n", + " dataset_names = []\n", + " # CoT session\n", + " CoT_trial_data = mouse + '_' + date + '_restructured_data.pkl'\n", + " CoT_photo_data = mouse + '_' + date + '_smoothed_signal.npy' \n", + " # RTC session \n", + " RTC_trial_data = mouse + '_' + date + '_RTC_restructured_data.pkl'\n", + " RTC_photo_data = mouse + '_' + date + '_RTC_smoothed_signal.npy' \n", + " dataset_names.extend([CoT_photo_data,CoT_trial_data, RTC_photo_data, RTC_trial_data])\n", + "\n", + "\n", + " for i, dataset_name in enumerate(dataset_names):\n", + " url = zenodo + dataset_name\n", + " dataset_path = '../../data/' + dataset_name\n", + " \n", + " if not exists(dataset_path):\n", + " print('Downloading ' + dataset_name)\n", + " urllib.request.urlretrieve(url, dataset_path)\n", + " #else:\n", + " # print(dataset_name + ' already in directory')\n", + "\n", + " \n", + " \n", + " if i == 0:\n", + " photometry_data = np.load(dataset_path) \n", + " \n", + " if i == 1: \n", + " trial_data = pd.read_pickle(dataset_path)\n", + " session_traces = yj_utils.SessionData(mouse, date, fiber_side, protocol, trial_data, photometry_data)\n", + " file_name = mouse + '_' + date + '_aligned_traces.p'\n", + " dataset_path = '../../data/' + file_name \n", + " pickle.dump(session_traces, open(dataset_path, \"wb\"))\n", + " if i == 2: \n", + " RTC_photometry_data = np.load(dataset_path) \n", + " \n", + " if i == 3:\n", + " RTC_trial_data = pd.read_pickle(dataset_path)\n", + " RTC_data = yj_utils.ZScoredTraces_RTC(RTC_trial_data, RTC_photometry_data, x_range)\n", + "\n", + " \n", + " if session_traces.protocol == 'SOR':\n", + " APE_aligned_data = decimate(session_traces.SOR_choice.contra_data.mean_trace, 10)\n", + " APE_time = decimate(session_traces.SOR_choice.contra_data.time_points, 10)\n", + " APE_sem_traces = decimate(session_traces.SOR_choice.contra_data.sorted_traces,10)\n", + " else:\n", + " APE_aligned_data = decimate(session_traces.choice.contra_data.mean_trace, 10)\n", + " APE_time = decimate(session_traces.choice.contra_data.time_points,10)\n", + " APE_sem_traces = decimate(session_traces.choice.contra_data.sorted_traces,10)\n", + " RTC_aligned_data = decimate(RTC_data.mean_trace, 10)\n", + " RTC_time = decimate(RTC_data.time_points, 10)\n", + "\n", + " if m == 0:\n", + " APE_traces = np.zeros((nr_mice, len(APE_aligned_data)))\n", + " RTC_traces = np.zeros((nr_mice, len(RTC_aligned_data)))\n", + " APE_sem_traces_upper = np.zeros((nr_mice, len(APE_aligned_data)))\n", + " APE_sem_traces_lower = np.zeros((nr_mice, len(APE_aligned_data)))\n", + " RTC_sem_traces_upper = np.zeros((nr_mice, len(RTC_aligned_data)))\n", + " RTC_sem_traces_lower = np.zeros((nr_mice, len(RTC_aligned_data)))\n", + " APE_peak_values = []\n", + " RTC_peak_values = []\n", + "\n", + " APE_traces[m,:] = APE_aligned_data\n", + " APE_sem_traces_lower[m,:], APE_sem_traces_upper[i,:] = yj_utils.calculate_error_bars(APE_aligned_data, APE_sem_traces,\n", + " error_bar_method='sem')\n", + " RTC_traces[m,:] = RTC_aligned_data\n", + " RTC_sem_traces = decimate(RTC_data.sorted_traces,10)\n", + " RTC_sem_traces_lower[m,:], RTC_sem_traces_upper[m,:] = yj_utils.calculate_error_bars(RTC_aligned_data, RTC_sem_traces,\n", + " error_bar_method='sem')\n", + " # get the peak values: # APE_time: 16000 datapoints, half: 8000 datapoints = time 0, only consider time after 0\n", + " start_inx = 8000\n", + " APE_range = APE_aligned_data[start_inx:start_inx+8000]\n", + " APE_time_range = APE_time[start_inx:start_inx+8000]\n", + " RTC_range = RTC_aligned_data[start_inx:start_inx+8000]\n", + "\n", + " APE_peak_index = np.argmax(APE_range) # from time 0 to 8s\n", + " APE_peak_time = APE_time_range[APE_peak_index]\n", + " APE_peak_value = APE_range[APE_peak_index]\n", + " RTC_peak_value = RTC_range[APE_peak_index]\n", + " APE_peak_values.append(APE_peak_value)\n", + " RTC_peak_values.append(RTC_peak_value)\n", + "\n", + "# calculate mean and sem across mice:\n", + "APE_mean_trace = np.mean(APE_traces, axis=0)\n", + "RTC_mean_trace = np.mean(RTC_traces, axis=0)\n", + "APE_sem_trace = np.std(APE_traces, axis=0)/np.sqrt(nr_mice)\n", + "RTC_sem_trace = np.std(RTC_traces, axis=0)/np.sqrt(nr_mice)\n", + "\n", + "figure = yj_plot.plot_SF5TU(APE_mean_trace, RTC_mean_trace, APE_sem_trace, RTC_sem_trace, APE_peak_values, RTC_peak_values, APE_time, RTC_time)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1314abc-abea-4346-8fe6-3c62f32e79f3", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/YJ_analysis_utils.py b/scripts/YJ_analysis_utils.py new file mode 100644 index 0000000..5670a29 --- /dev/null +++ b/scripts/YJ_analysis_utils.py @@ -0,0 +1,431 @@ + +from scipy import stats +import numpy as np + + +class SessionData(object): + def __init__(self, mouse, date, fiber_side, protocol, trial_data, photometry_data): + self.mouse = mouse + self.date = date + self.fiber_side = fiber_side + self.protocol = protocol + self.choice = None + self.cue = None + self.reward = None + + if protocol != 'SOR': + self.choice = ChoiceAlignedData(self, trial_data, photometry_data) + #self.cue = CueAlignedData(self,trial_data, photometry_data, save_traces=True) + #self.reward = RewardAlignedData(self, trial_data, photometry_data, save_traces=True) + elif self.protocol == 'SOR': + self.SOR_choice = SORChoiceAlignedData(self, trial_data, photometry_data) + + + +class ChoiceAlignedData(object): + """ + Traces for standard analysis aligned to choice (=movement from center port to left or right port) + """ + + def __init__(self, session_data, trial_data, photometry_data): + + # "RESPONSE": RIGHT = 2, LEFT = 1: hence ipsi and contra need to be assigned accordingly: + fiber_options = np.array(['left', 'right']) # left = (0+1) = 1; right = (1+1) == 2 + ipsi_fiber_side_numeric = (np.where(fiber_options == session_data.fiber_side)[0] + 1)[0] # if fiber on right ipsi = 2; if fiber on left ipsi = 1 + contra_fiber_side_numeric = (np.where(fiber_options != session_data.fiber_side)[0] + 1)[0] # if fiber on right contra = 1, if fiber on left contra = 2 + + params = {'state_type_of_interest': 5, + 'outcome': 2, # 2 = doesn't matter for choice aligned data + 'no_repeats': 1, # 1 = no repeats allowed + 'last_response': 0, # the previous trial doesn't matter + 'align_to': 'Time start', + 'instance': -1, # -1 = last instance; 1 = first instance + 'plot_range': [-6, 6], + 'first_choice_correct': 2, # 2 = doesnt matter + 'SOR': 0, # 0 = nonSOR; 2 = doesnt matter, 1 = SOR + 'psycho': 0, # only Trial type 1 and 7 (no intermediate values / psychometric sounds) + 'LRO': 0, # 0 = nonLRO; + 'LargeRewards': 0, # 1 = LR + 'Omissions': 0, # 1 = Omission + 'Silence': 0, # 1 = Silence + 'cue': None} + + self.ipsi_data = ZScoredTraces(trial_data, photometry_data, params, ipsi_fiber_side_numeric, ipsi_fiber_side_numeric) + self.contra_data = ZScoredTraces(trial_data, photometry_data, params, contra_fiber_side_numeric, contra_fiber_side_numeric) + +class SORChoiceAlignedData(object): + """ + Traces for SOR analysis: aligned to movement for trials when cue has been played on return already + """ + + def __init__(self, session_data, trial_data, photometry_data): + fiber_options = np.array(['left', 'right']) # left = (0+1) = 1; right = (1+1) == 2 + contra_fiber_side_numeric = (np.where(fiber_options != session_data.fiber_side)[0] + 1)[0] # if fiber on right contra = 1, if fiber on left contra = 2 + + params = {'state_type_of_interest': 5, + 'outcome': 2, # 2 = doesn't matter for choice aligned data + 'no_repeats': 1, + 'last_response': 0, # doesnt matter for choice aligned data + 'align_to': 'Time start', + 'instance': -1, # last instance + 'plot_range': [-6, 6], + 'first_choice_correct': 2, + 'SOR': 1, + 'psycho': 0, + 'LRO': 0, + 'LargeRewards': 0, + 'Omissions': 0, + 'Silence': 0, + 'cue': None} + + self.contra_data = ZScoredTraces(trial_data, photometry_data, params, contra_fiber_side_numeric, contra_fiber_side_numeric) + # no ipsi data for SOR trials as ipsi trials were classic 2AC trials. + + + + +class ZScoredTraces(object): + def __init__(self, trial_data, dff, params, response, first_choice): + self.params = HeatMapParams(params, response, first_choice) + self.time_points, self.mean_trace, self.sorted_traces, self.reaction_times, self.state_name, self.title, self.sorted_next_poke, self.trial_nums, self.event_times, self.outcome_times = find_and_z_score_traces( + trial_data, dff, self.params, sort=False) + +class HeatMapParams(object): + def __init__(self, params, response, first_choice): + self.state = params['state_type_of_interest'] + self.outcome = params['outcome'] + self.response = response + self.last_response = params['last_response'] + self.align_to = params['align_to'] + self.other_time_point = np.array(['Time start', 'Time end'])[np.where(np.array(['Time start', 'Time end']) != params['align_to'])] + self.instance = params['instance'] + self.plot_range = params['plot_range'] + self.no_repeats = params['no_repeats'] + self.first_choice_correct = params['first_choice_correct'] + self.first_choice = first_choice + self.cue = params['cue'] + self.SOR = params['SOR'] + self.psycho = params['psycho'] + self.LRO = params['LRO'] + self.LR = params['LargeRewards'] + self.O = params['Omissions'] + self.S = params['Silence'] + + +def find_and_z_score_traces(trial_data, dff, params, norm_window=8, sort=False, get_photometry_data=True): + response_names = ['both left and right', 'left', 'right'] + outcome_names = ['incorrect', 'correct', 'both correct and incorrect'] + title = '' + + # Filter trial data according to selection of special trials (e.g. SOR, LRO, psychometric etc) + if params.SOR == 0: + events_of_int = getNonSORtrials(trial_data) + elif params.SOR == 1: + events_of_int = getSORtrials(trial_data) + elif params.SOR == 2: + events_of_int = trial_data + if params.psycho == 0: + events_of_int = getNonPsychotrials(events_of_int) + if params.LRO == 0: + events_of_int = getNonLROtrials(events_of_int) + if params.S == 0: + events_of_int = getNonSilenceTrials(events_of_int) + + + + # 1) State type (e.g. corresp. State name = CueDelay, WaitforResponse...) + events_of_int = events_of_int.loc[(events_of_int['State type'] == params.state)] # State type = number of state of interest + state_name = events_of_int['State name'].values[0] + title = title + 'State type = ' + str(params.state) + ' = state_name ' + state_name + ';' + # -------------- + + # 2) Response, trials to the left or to the right side + if params.response != 0: # 0 = don't care, 1 = left, 2 = right, selection of ipsi an contra side depends on fiber side + events_of_int = events_of_int.loc[events_of_int['Response'] == params.response] + title = title + ' Response = ' + str(params.response) + ';' + else: + print('Response = 0, so both left and right responses are considered') + # -------------- + + # 3) First and last response: + if params.first_choice != 0: + events_of_int = events_of_int.loc[events_of_int['First response'] == params.first_choice] + title = title + ' 1st response = ' + str(params.first_choice) + ';' + if params.last_response != 0: + events_of_int = events_of_int.loc[events_of_int['Last response'] == params.last_response] + title = title + ' last response = ' + str(params.last_response) + ';' + # -------------- + + # 4) Outcome: + if not params.outcome == 2: # 2 = outcome doesn't matter + events_of_int = events_of_int.loc[events_of_int['Trial outcome'] == params.outcome] + title = title + ' Outcome = ' + str(params.outcome) + ';' + # -------------- + + # 5) Cues / Sounds: + if params.cue == 'high': + events_of_int = events_of_int.loc[events_of_int['Trial type'] == 7] + title = title + ' Cue = high;' + elif params.cue == 'low': + events_of_int = events_of_int.loc[events_of_int['Trial type'] == 1] + title = title + ' Cue = low;' + # -------------- + + # 6) Instance in State & Repeats: + if params.instance == -1: # Last time in State + events_of_int = events_of_int.loc[ + (events_of_int['Instance in state'] / events_of_int['Max times in state'] == 1)] + title = title + ' instance (' + str(params.instance) + ') last time in state (no matter the repetitions);' + elif params.instance == 1: # First time in State + events_of_int = events_of_int.loc[(events_of_int['Instance in state'] == 1)] + title = title + ' instance (' + str(params.instance) + ') first time in state;' + if params.no_repeats == 1: + events_of_int = events_of_int.loc[events_of_int['Max times in state'] == 1] + title = title + ' no repetitions allowed (' + str(params.no_repeats) + ')' + # -------------- + + # 7) First choice directly in/correct? + if params.first_choice_correct == 1: # only first choice correct + events_of_int = events_of_int.loc[ + (events_of_int['First choice correct'] == 1)] + title = title + ' 1st choice correct (' + str(params.first_choice_correct) + ') only' + elif params.first_choice_correct == -1: # only first choice incorrect + events_of_int = events_of_int.loc[np.logical_or( + (events_of_int['First choice correct'] == 0), (events_of_int['First choice correct'].isnull()))] + title = title + ' 1st choice incorrect (' + str(params.first_choice_correct) + ') only' + if events_of_int['State type'].isin([5.5]).any(): # first incorrect choice? + events_of_int = events_of_int.loc[events_of_int['First choice correct'].isnull()] + elif params.first_choice_correct == 0: # first choice incorrect + events_of_int = events_of_int.loc[(events_of_int['First choice correct'] == 0)] + # -------------- + + event_times = events_of_int[params.align_to].values # start or end of state of interest time points + trial_nums = events_of_int['Trial num'].values + trial_starts = events_of_int['Trial start'].values + trial_ends = events_of_int['Trial end'].values + + other_event = np.asarray(np.squeeze(events_of_int[params.other_time_point].values) - np.squeeze(events_of_int[params.align_to].values)) + # for ex. time end - time start of state of interest + + last_trial = np.max(trial_data['Trial num']) # absolutely last trial in session + last_trial_num = events_of_int['Trial num'].unique()[-1] # last trial that is considered in analysis meeting params requirements + events_reset_index = events_of_int.reset_index(drop=True) + last_trial_event_index = events_reset_index.loc[(events_reset_index['Trial num'] == last_trial_num)].index + # index of the last event in the last trial that is considered in analysis meeting params requirements + next_centre_poke = get_next_centre_poke(trial_data, events_of_int, last_trial_num == last_trial) + trial_starts = get_first_poke(trial_data, events_of_int) + absolute_outcome_times = get_outcome_time(trial_data, events_of_int) + relative_outcome_times = absolute_outcome_times - event_times + + if get_photometry_data == True: + next_centre_poke[last_trial_event_index] = events_reset_index[params.align_to].values[ + last_trial_event_index] + 1 # so that you can find reward peak + + next_centre_poke_norm = next_centre_poke - event_times + event_photo_traces = get_photometry_around_event(event_times, dff, pre_window=norm_window, + post_window=norm_window) + norm_traces = stats.zscore(event_photo_traces.T, axis=0) + + + if other_event.size == 1: + print('Only one event for ' + title + ' so no sorting') + sort = False + elif len(other_event) != norm_traces.shape[1]: + other_event = other_event[:norm_traces.shape[1]] + print('Mismatch between #events and #other_event') + if sort: + arr1inds = other_event.argsort() + sorted_other_event = other_event[arr1inds[::-1]] + sorted_traces = norm_traces.T[arr1inds[::-1]] + sorted_next_poke = next_centre_poke_norm[arr1inds[::-1]] + else: + sorted_other_event = other_event + sorted_traces = norm_traces.T + sorted_next_poke = next_centre_poke_norm + + time_points = np.linspace(-norm_window, norm_window, norm_traces.shape[0], endpoint=True, retstep=False, + dtype=None, axis=0) + mean_trace = np.mean(sorted_traces, axis=0) + + return time_points, mean_trace, sorted_traces, sorted_other_event, state_name, title, sorted_next_poke, trial_nums, event_times, relative_outcome_times + + + + + +def getSORtrials(trial_data): + trials_2AC = trial_data[trial_data['State name'] == 'WaitForPoke'] + trial_num_2AC = trials_2AC['Trial num'].unique() + trials_SOR = trial_data + for trialnum in trial_num_2AC: + trials_SOR = trials_SOR[trials_SOR['Trial num']!=trialnum] + return(trials_SOR) + +def getNonSORtrials(trial_data): + trialnums_trial_data = trial_data['Trial num'].unique() # all trial numbers + events_2AC = trial_data[trial_data['State name'] == 'WaitForPoke'] # 2AC WaitForPoke events + trialnums_2AC = events_2AC['Trial num'].unique() # 2AC trial numbers + trial_data_2AC = trial_data + for trial_all in trialnums_trial_data: + include = 0 + for trial_2AC in trialnums_2AC: + if trial_all == trial_2AC: + include = 1 + if include == 0: + trial_data_2AC = trial_data_2AC[trial_data_2AC['Trial num'] != trial_all] + return trial_data_2AC + +def getNonSilenceTrials(trial_data): + total_trials = len(trial_data['Trial num'].unique()) + trials_2AC = trial_data[trial_data['State name'] == 'WaitForPoke'] + trial_num_2AC = len(trials_2AC['Trial num'].unique()) + if trial_num_2AC == total_trials: # This is not a Sound-On-Return session + trial_data_NonSilence = trial_data + trial_data_NonSilence = trial_data_NonSilence[trial_data_NonSilence['Sound type'] != 1] + new_trial_num = trial_data_NonSilence['Trial num'].unique() + #if len(new_trial_num) < total_trials: + #print('non-silent trials ' + str(len(new_trial_num)) + ' out of ' + str(total_trials) + ' trials') + #print(trial_data_NonSilence['Sound type'].unique()) + else: + trial_data_NonSilence = trial_data + return trial_data_NonSilence + + +def getNonPsychotrials(trial_data): + for sound_number in [2, 3, 4, 5, 6]: + trial_data = trial_data[trial_data['Trial type'] != sound_number] # Trial type == 1 or 7 == classic high and low frequency + return trial_data + +def getNonLROtrials(trial_data): + LR_trial_numbers = trial_data[(trial_data['State type'] == 12) | (trial_data['State type'] == 13)][ + 'Trial num'].values # 13 = RightLargeReward, 12 = LeftLargeReward + O_trial_numbers = trial_data[(trial_data['State type'] == 10) & (trial_data['State name'] == 'Omission')]['Trial num'].values # 10 = Omission or REturnCuePlay + LRO_trial_numbers = np.concatenate((LR_trial_numbers, O_trial_numbers)) + #NonLRO_trial_numbers = trial_data[~trial_data['Trial num'].isin(LRO_trial_numbers)]['Trial num'].values + NonLRO_trial_data = trial_data + for trial_num in trial_data['Trial num'].unique(): + if trial_num in LRO_trial_numbers: + NonLRO_trial_data = NonLRO_trial_data[NonLRO_trial_data['Trial num'] != trial_num] + return NonLRO_trial_data + +def get_next_centre_poke(trial_data, events_of_int, last_trial): + ''' + This function returns the time of the first centre poke in the subsequent trial for each event of interest. + + last_trial is a boolean that is true if the last trial in the session is included in the events of interest + ''' + + next_centre_poke_times = np.zeros(events_of_int.shape[0]) + events_of_int = events_of_int.reset_index(drop=True) + for i, event in events_of_int.iterrows(): + trial_num = event['Trial num'] + if trial_num == trial_data['Trial num'].values[-1]: + next_centre_poke_times[i] = events_of_int['Trial end'].values[i] + 2 + else: + next_trial_events = trial_data.loc[(trial_data['Trial num'] == trial_num + 1)] + wait_for_poke_state = next_trial_events.loc[(next_trial_events['State type'] == 2)] # wait for pokes + + if(len(wait_for_poke_state) > 0): # Classic 2AC: + wait_for_pokes = next_trial_events.loc[(next_trial_events['State type'] == 2)] # wait for pokes + next_wait_for_poke = wait_for_pokes.loc[(wait_for_pokes['Instance in state'] == 1)] # first wait for poke + next_centre_poke_times[i] = next_wait_for_poke['Time end'].values[0] # time of first wait for poke ending + elif len(wait_for_poke_state) == 0: # SOR: (SOR trials don't have WaitForPoke state) + wait_for_pokes = next_trial_events.loc[(next_trial_events['State type'] == 3)] # CueDelay + next_wait_for_poke = wait_for_pokes.loc[(wait_for_pokes['Instance in state'] == 1)] # first wait for poke + next_centre_poke_times[i] = next_wait_for_poke['Time start'].values[0] # start time of first poke + + + if last_trial: # last trial in events of interest == last trial in session, last_tial is true or false, not a number + next_centre_poke_times[-1] = events_of_int['Trial end'].values[-1] + 2 + else: # last trial in events of interest != last trial in session + event = events_of_int.tail(1) + trial_num = event['Trial num'].values[0] + next_trial_events = trial_data.loc[(trial_data['Trial num'] == trial_num + 1)] + + wait_for_poke_state = next_trial_events.loc[(next_trial_events['State type'] == 2)] # wait for pokes + if (len(wait_for_poke_state) > 0): # Classic 2AC: + wait_for_pokes = next_trial_events.loc[(next_trial_events['State type'] == 2)] + next_wait_for_poke = wait_for_pokes.loc[(wait_for_pokes['Instance in state'] == 1)] + next_centre_poke_times[-1] = next_wait_for_poke['Time end'].values[0] # end time of wait for poke + elif len(wait_for_poke_state) == 0: # SOR: (SOR trials don't have WaitForPoke state) + wait_for_pokes = next_trial_events.loc[(next_trial_events['State type'] == 3)] # CueDelay + next_wait_for_poke = wait_for_pokes.loc[(wait_for_pokes['Instance in state'] == 1)] # first wait for poke + next_centre_poke_times[i] = next_wait_for_poke['Time start'].values[0] # start time of first poke + return next_centre_poke_times + + +def get_first_poke(trial_data, events_of_int): # get first poke in each trial of events of interest + trial_numbers = events_of_int['Trial num'].unique() + next_centre_poke_times = np.zeros(events_of_int.shape[0]) + events_of_int = events_of_int.reset_index(drop=True) + for trial_num in trial_numbers: + event_indx_for_that_trial = events_of_int.loc[(events_of_int['Trial num'] == trial_num)].index + trial_events = trial_data.loc[(trial_data['Trial num'] == trial_num)] + wait_for_pokes = trial_events.loc[(trial_events['State type'] == 2)] + if len(wait_for_pokes) > 0: #Classic 2AC: + next_wait_for_poke = wait_for_pokes.loc[(wait_for_pokes['Instance in state'] == 1)] + #next_centre_poke_times[event_indx_for_that_trial] = next_wait_for_poke['Time end'].values[0]-1 #why -1 in FG code? + next_centre_poke_times[event_indx_for_that_trial] = next_wait_for_poke['Time end'].values[0] + + elif len(wait_for_pokes) == 0: #SOR: (SOR trials don't have WaitForPoke state) + next_wait_for_poke = trial_events.loc[(trial_events['State type'] == 3) & (trial_events['Instance in state'] == 1)] #First CueDelay + next_centre_poke_times[event_indx_for_that_trial] = next_wait_for_poke['Time start'].values[0] + + return next_centre_poke_times + +def get_outcome_time(trial_data, events_of_int): # returns the time of the outcome of the current trial, indep of rewarded or punished + trial_numbers = events_of_int['Trial num'].values + outcome_times = [] + for event_trial_num in range(len(trial_numbers)): + trial_num = trial_numbers[event_trial_num] + other_trial_events = trial_data.loc[(trial_data['Trial num'] == trial_num)] + choices = other_trial_events.loc[(other_trial_events['State type'] == 5)] # 5 is the state type for choices / wait for response + max_times_in_state_choices = choices['Max times in state'].unique() # all values in max times in state available for this trial an state type + choice = choices.loc[(choices['Instance in state'] == max_times_in_state_choices)] # last time wait for response + outcome_times.append(choice['Time end'].values[0]) + return outcome_times + + +def get_photometry_around_event(all_trial_event_times, demodulated_trace, pre_window=5, post_window=5, sample_rate=10000): + num_events = len(all_trial_event_times) + event_photo_traces = np.zeros((num_events, sample_rate*(pre_window + post_window))) + for event_num, event_time in enumerate(all_trial_event_times): + plot_start = int(round(event_time*sample_rate)) - pre_window*sample_rate + plot_end = int(round(event_time*sample_rate)) + post_window*sample_rate + event_photo_traces[event_num, :] = demodulated_trace[plot_start:plot_end] + return event_photo_traces + + +class ZScoredTraces_RTC(object): + def __init__(self, trial_data, df, x_range): + events_of_int = trial_data.loc[(trial_data['State type'] == 2)] # cue = state type 2 + event_times = events_of_int['Time start'].values + event_photo_traces = get_photometry_around_event(event_times, df, pre_window=5, post_window=5) + norm_traces = stats.zscore(event_photo_traces.T, axis=0) + sorted_traces = norm_traces.T + x_vals = np.linspace(x_range[0],x_range[1], norm_traces.shape[0], endpoint=True, retstep=False, dtype=None, axis=0) + y_vals = np.mean(sorted_traces, axis=0) + self.sorted_traces = sorted_traces + self.mean_trace = y_vals + self.time_points = x_vals + self.params = RTC_params(x_range) + self.reaction_times = None + self.events_of_int = events_of_int + + +class RTC_params(object): + def __init__(self, x_range): + self.plot_range = x_range + +def calculate_error_bars(mean_trace, data, error_bar_method='sem'): + if error_bar_method == 'sem': + sem = stats.sem(data, axis=0) + lower_bound = mean_trace - sem + upper_bound = mean_trace + sem + elif error_bar_method == 'ci': + lower_bound, upper_bound = bootstrap(data, n_boot=1000, ci=68) + return lower_bound, upper_bound + +def test2(a, b): + c = a * b + return c \ No newline at end of file diff --git a/scripts/__init__.py b/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/scripts/yj_plotting.py b/scripts/yj_plotting.py new file mode 100644 index 0000000..4ff9d7d --- /dev/null +++ b/scripts/yj_plotting.py @@ -0,0 +1,50 @@ + +import matplotlib.pyplot as plt +import numpy as np + +def plot_SF5TU(APE_mean_trace, RTC_mean_trace, APE_sem_trace, RTC_sem_trace, APE_peak_values, RTC_peak_values, APE_time, RTC_time): + x_range = [-2, 3] + y_range = [-0.75, 1.5] + plt.rcParams["figure.figsize"] = (3, 6) + plt.rcParams['axes.spines.right'] = False + plt.rcParams['axes.spines.top'] = False + + fig, ax = plt.subplots(1, 2 , figsize=(6, 3)) # width, height + + # Plot with average traces: + ax[0].axvline(0, color='#808080', linewidth=0.25, linestyle='dashdot') + ax[0].plot(APE_time, APE_mean_trace, lw=2, color='#3F888F', label = 'Choice') + ax[0].fill_between(APE_time, APE_mean_trace - APE_sem_trace, APE_mean_trace + APE_sem_trace, color='#7FB5B5', + linewidth=1, alpha=0.6) + + ax[0].plot(RTC_time, RTC_mean_trace, lw=2, color='#e377c2', label = 'Sound') + ax[0].fill_between(RTC_time, RTC_mean_trace - RTC_sem_trace, RTC_mean_trace + RTC_sem_trace, facecolor='#e377c2', + linewidth=1, alpha=0.4) + ax[0].legend(loc='upper right', frameon=False) + ax[0].set_ylim(y_range) + ax[0].set_ylabel('dLight z-score') + ax[0].set_xlabel('Time (s)') + ax[0].set_xlim(x_range) + ax[0].set_ylim(y_range) + ax[0].yaxis.set_ticks([-0.5, 0, 0.5, 1, 1.5]) + ax[0].xaxis.set_ticks([-2, 0, 2]) + plt.tight_layout() + + # dotplot with peak values: + for i in range(0, len(APE_peak_values)): + x_val = [0, 1] + y_val = [APE_peak_values[i], RTC_peak_values[i]] + ax[1].plot(x_val, y_val, color='grey', linewidth=0.5) + ax[1].scatter(0, APE_peak_values[i], color='#3F888F', s=100, alpha=1) + ax[1].scatter(1, RTC_peak_values[i], color='#e377c2', s=100, alpha=1) + + x_text_values = ['Choice', 'Sound'] + ax[1].set_xticks([0, 1]) + ax[1].set_ylabel('dLight z-score') + ax[1].set_xlim(-0.2, 1.2) + ax[1].set_xticklabels(x_text_values) + ax[1].set_ylim(-0.5, 1.5) + ax[1].yaxis.set_ticks([-0.5, 0.5, 1.5]) + fig.tight_layout(pad=2) + + return fig \ No newline at end of file