time_series_analysis_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## a useful function to generate a data list for further analysis\n",
    "import os,json,sys\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from sorting_data_pairwise_comparison import calculate_relative_position,generate_points_within_rectangles\n",
    "import scipy.stats as st\n",
    "\n",
    "##need to add this additional cell because useful tools are in another folder. Need to integrate these two folders one day\n",
    "current_working_directory = Path.cwd()\n",
    "parent_dir = current_working_directory.resolve().parents[0]\n",
    "sys.path.insert(0, str(parent_dir) + \"\\\\utilities\")\n",
    "from useful_tools import select_animals_gpt,find_file\n",
    "from data_cleaning import findLongestConseqSubseq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Session 0.2: Load analysis methods in python dictionary form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "json_file = \"./analysis_methods_dictionary.json\"\n",
    "with open(json_file, \"r\") as f:\n",
    "    analysis_methods = json.loads(f.read())\n",
    "    \n",
    "#Put the folder of your Unity folder below\n",
    "#thisDataset =\"D:/MatrexVR_Swarm_Data/RunData\"\n",
    "#thisDataset =\"D:/MatrexVR_blackbackground_Data/RunData\"\n",
    "thisDataset =\"D:/MatrexVR_grass1_Data/RunData\"\n",
    "#thisDataset =\"D:/MatrexVR_2024_Data/RunData\"\n",
    "#parameter name means independent variable in the experiment\n",
    "variable_name='mu'\n",
    "#variable_name='initial_position'\n",
    "#variable_name='agent_speed'\n",
    "#check trace in trial 115 from VR1_2024-11-16_155242_score_full, maybe there is a jump"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Session 1.0: select animals based on condition and return which a directory list and a list of vr rig number to specify which animal to analyse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the path to your Excel file\n",
    "dir_list = []\n",
    "file_type=\".h5\"\n",
    "using_google_sheet=True\n",
    "sheet_name = \"Unity_MatrexVR\"\n",
    "experiment_name=analysis_methods.get(\"experiment_name\")\n",
    "# if type(thisDataset) == str:\n",
    "#     thisDataset = Path(thisDataset)\n",
    "if analysis_methods.get(\"load_individual_data\") == True:\n",
    "    if using_google_sheet==True:\n",
    "        # database_id = \"1UL4eEUrQMapx9xz11-IyOSlPBcep3I9vBJ2uGgVudb8\"\n",
    "        #         #https://docs.google.com/spreadsheets/d/1UL4eEUrQMapx9xz11-IyOSlPBcep3I9vBJ2uGgVudb8/edit?usp=sharing\n",
    "        # url = f\"https://docs.google.com/spreadsheets/d/{database_id}/gviz/tq?tqx=out:csv&sheet={sheet_name}\"\n",
    "        database_id = \"1UL4eEUrQMapx9xz11-IyOSlPBcep1I9vBJ2uGgVudb8\"\n",
    "                #https://docs.google.com/spreadsheets/d/1UL4eEUrQMapx9xz11-IyOSlPBcep1I9vBJ2uGgVudb8/edit?usp=sharing\n",
    "        url = f\"https://docs.google.com/spreadsheets/d/{database_id}/gviz/tq?tqx=out:csv&sheet={sheet_name}\"\n",
    "        #df = pd.read_excel(url, engine='openpyxl')## use this function if the file is not google sheet but uploaded excel file\n",
    "        df = pd.read_csv(url)\n",
    "    else:\n",
    "        excel_file_path = \"Z:/DATA/experiment_trackball_Optomotor/Locusts Management.xlsx\"\n",
    "        print(f\"using a database {excel_file_path} from the server but this file might be outdated\")\n",
    "        # Create a 'with' statement to open and read the Excel file\n",
    "        with pd.ExcelFile(excel_file_path) as xls:\n",
    "            # Read the Excel sheet into a DataFrame with the sheet name (folder name)\n",
    "            df = pd.read_excel(xls, sheet_name)\n",
    "        ##list up the conditions and answers as strings for input argument to select animal. One condition must pair with one answer\n",
    "    if analysis_methods.get(\"select_animals_by_condition\") == True:\n",
    "       #animal_of_interest=select_animals_gpt(df,\"Independent variable (list up all of them in the experiment)\",\"gregarious_leader_grass\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        #animal_of_interest=select_animals_gpt(df,\"Independent variable (list up all of them in the experiment)\",\"gregarious_leader_black\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        animal_of_interest=select_animals_gpt(df,\"Independent variable1\",variable_name,\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        #animal_of_interest=select_animals_gpt(df,\"Independent variable1\",variable_name,\"Independent variable2\",\"choice_vr_locust_black_locust\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        #animal_of_interest=select_animals_gpt(df,\"Independent variable1\",variable_name,\"Independent variable2\",\"choice_vr_locust_sta_black_locust\",\"Background\",\"grass1\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        \n",
    "        #animal_of_interest=select_animals_gpt(df,\"Independent variable1\",variable_name,\"Independent variable2\",\"closed_loop_sta_black_locust_open_loop_sta_black_locust\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "        #animal_of_interest=select_animals_gpt(df,\"Independent variable1\",variable_name,\"Independent variable2\",\"sta_black_locust_2dir_3_initial_position\",\"Excluding this animal from analysis (Usually when animals die or molt, T/F)\",\"F\")\n",
    "    else:\n",
    "        animal_of_interest=df\n",
    "    folder_name=animal_of_interest[\"folder name\"].values\n",
    "\n",
    "    dir_tile=np.tile(thisDataset, (len(folder_name), 1))\n",
    "    vr_no=animal_of_interest[\"VR number\"].values\n",
    "    dir_list = [''.join([x[0], '/', y]) for x,y in zip(dir_tile,folder_name)]\n",
    "    #dir_dict = itertools.zip_longest(dir_list, vr_no.tolist())\n",
    "    dir_dict = zip(dir_list, vr_no.tolist())\n",
    "else:\n",
    "    for root, dirs, files in os.walk(thisDataset):\n",
    "        for folder in dirs:\n",
    "            folder_path=os.path.join(root,folder)\n",
    "            if any(name.endswith(file_type) for name in os.listdir(folder_path)):\n",
    "                dir_list.append(folder_path.replace(\"\\\\\", \"/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['D:/MatrexVR_grass1_Data/RunData/20240824_143943',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_143943',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_143943',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_163554',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_163554',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_184143',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_184143',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240824_184143',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_132616',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_151648',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_151648',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_170748',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_170748',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240825_170748',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_131447',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_150826',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_150826',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_150826',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_174030',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_174030',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_174030',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240826_174030',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_151537',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_151537',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_151537',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_151537',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_171158',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_171158',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_171158',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_171158',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_190731',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_190731',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_190731',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240904_190731',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_134125',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_134125',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_134125',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_134125',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_155137',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_155137',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_155137',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_155137',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_174813',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_174813',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_174813',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240905_174813',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_142802',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_142802',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_142802',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_170446',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_170446',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_190839',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_190839',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240907_190839',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_125638',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_125638',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_125638',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_125638',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_150715',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_150715',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_150715',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_150715',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_174232',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_174232',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_174232',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_193840',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_193840',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_193840',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20240908_193840',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_162147',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_162147',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_162147',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_162147',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_182305',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_182305',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_182305',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_201706',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_201706',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241012_201706',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_144457',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_144457',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_144457',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_164304',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_164304',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_164304',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_164304',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_184451',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_184451',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_184451',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241013_184451',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_134432',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_134432',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_134432',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_134432',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_155538',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_155538',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_155538',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_155538',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_175759',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_175759',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_175759',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241014_175759',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_123054',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_123054',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_123054',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_123054',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_144142',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_144142',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_144142',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_144142',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_165438',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_165438',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241110_165438',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_122609',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_122609',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_122609',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_143727',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_143727',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_143727',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_143727',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_164808',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_164808',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_164808',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241111_164808',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_123949',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_123949',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_123949',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_123949',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_150308',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_150308',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_150308',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_171922',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241112_171922',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_134457',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_134457',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_155210',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_155210',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_155210',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_180011',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_180011',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_180011',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241116_180011',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_130954',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_130954',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_130954',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_152042',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_152042',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_152042',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_152042',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_173052',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_173052',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_173052',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241117_173052',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_130021',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_130021',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_130021',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_130021',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_151102',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_151102',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_151102',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_151102',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_172046',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_172046',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_172046',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241118_172046',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_131605',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_131605',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_131605',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_131605',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_152639',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_152639',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_152639',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_174120',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241201_174120',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_131919',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_131919',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_131919',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_131919',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_152804',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_152804',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_174716',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241202_174716',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_132141',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_132141',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_132141',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_132141',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_160127',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_160127',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_160127',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_190928',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_190928',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_190928',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241203_190928',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_141818',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_141818',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_141818',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_141818',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_160907',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_160907',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_160907',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_175920',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_175920',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_175920',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241207_175920',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_134852',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_134852',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_134852',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_134852',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_155704',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_155704',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_155704',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_155704',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_180456',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_180456',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241225_180456',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_134945',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_134945',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_134945',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_155750',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_155750',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_155750',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_155750',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_182032',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_182032',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_182032',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241226_182032',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_131107',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_131107',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_131107',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_175417',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_175417',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20241227_175417',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20250103_140602',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20250103_140602',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20250103_140602',\n",
       " 'D:/MatrexVR_grass1_Data/RunData/20250103_140602']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_list=dir_list[32:]\n",
    "vr_no=vr_no[32:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_list=dir_list[:32]\n",
    "vr_no=vr_no[:32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(dir_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Session 1.1: introduce helper functions to make plot and calculate speed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 2.0: Pool animal's response together according to some criteria (from none criteria to criteria that can define the follow behaviour)\n",
    "Output1: a list 'follow_proportion_across_animals' showing the proportion of 'follow' time for each animal (across trials)\n",
    "\n",
    "Output2: a list 'relative_pos_all_animals' showing relative position between virtual and focal locusts across time. 1st and 2nd columns shown relative x and y, 3rd columns shown virtual animal's moving direction. 4th column shown the timestamp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_methods.update({\"follow_locustVR_criteria\": True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'experiment_name': 'choice',\n",
       " 'overwrite_curated_dataset': True,\n",
       " 'save_output': True,\n",
       " 'time_series_analysis': True,\n",
       " 'filtering_method': 'sg_filter',\n",
       " 'plotting_trajectory': False,\n",
       " 'plotting_event_distribution': True,\n",
       " 'distribution_with_entire_body': False,\n",
       " 'load_individual_data': True,\n",
       " 'select_animals_by_condition': True,\n",
       " 'active_trials_only': True,\n",
       " 'align_with_isi_onset': False,\n",
       " 'extract_follow_epoches': True,\n",
       " 'follow_locustVR_criteria': True,\n",
       " 'graph_colour_code': ['r', 'b', 'g', 'k', 'c', 'y', 'm', 'r'],\n",
       " 'follow_within_distance': 50,\n",
       " 'camera_fps': 100,\n",
       " 'trackball_radius_cm': 0.5,\n",
       " 'monitor_fps': 60,\n",
       " 'body_length': 4,\n",
       " 'growth_condition': 'G',\n",
       " 'analysis_window': [-10, 10]}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "analysis_methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relative_pos_all_animals=[]\n",
    "simulated_relative_pos_all_animals=[]\n",
    "trial_evaluation_across_animals=[]\n",
    "raster_across_animals=[]\n",
    "animal_id=0\n",
    "for this_dir,this_vr in zip(dir_list,vr_no):\n",
    "    if Path(this_dir).is_dir()==False:\n",
    "        continue\n",
    "    agent_pattern = f\"VR{this_vr}*agent_full.h5\"\n",
    "    xy_pattern = f\"VR{this_vr}*XY_full.h5\"\n",
    "    summary_pattern = f\"VR{this_vr}*score_full.h5\"\n",
    "    agent_file = find_file(Path(this_dir), agent_pattern)\n",
    "    focal_animal_file = find_file(Path(this_dir), xy_pattern)\n",
    "    summary_file = find_file(Path(this_dir), summary_pattern)\n",
    "    relative_pos,trial_evaluation_list,raster_pd,num_unfilled_gap,simulated_relative_pos=calculate_relative_position(summary_file,focal_animal_file,agent_file,analysis_methods)\n",
    "    if animal_id==0:\n",
    "        largest_unfilled_gap=num_unfilled_gap\n",
    "    elif num_unfilled_gap>largest_unfilled_gap:\n",
    "        largest_unfilled_gap=num_unfilled_gap\n",
    "    else:\n",
    "        pass\n",
    "    relative_pos_all_animals.append(relative_pos)\n",
    "    simulated_relative_pos_all_animals.append(simulated_relative_pos)\n",
    "    trial_evaluation=pd.concat(trial_evaluation_list)\n",
    "    trial_evaluation.insert(0, 'animal_id',np.repeat(animal_id,trial_evaluation.shape[0]))\n",
    "    trial_evaluation_across_animals.append(trial_evaluation)\n",
    "    raster_pd.insert(0, 'animal_id', np.repeat(animal_id,raster_pd.shape[0]))\n",
    "    raster_across_animals.append(raster_pd)\n",
    "    animal_id=animal_id+1\n",
    "analysis_methods['largest_unfilled_gap']= largest_unfilled_gap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_evaluation=pd.concat(trial_evaluation_across_animals)\n",
    "##for grass background 1 agent\n",
    "# fair_follower_threshold=0.2\n",
    "# good_follower_threshold=0.3\n",
    "##for grass background 2 agent\n",
    "# fair_follower_threshold=0.1667\n",
    "# good_follower_threshold=0.25\n",
    "##for blackground\n",
    "# fair_follower_threshold=0.05\n",
    "# good_follower_threshold=0.25\n",
    "follow_time_tbt=all_evaluation['num_follow_epochs']/all_evaluation['number_frames']\n",
    "all_evaluation[\"follow_ratio_previous_trial\"]=pd.concat([pd.Series(np.nan),follow_time_tbt[:-1]],ignore_index=True).to_list()\n",
    "all_evaluation[\"follow_ratio_next_trial\"]=pd.concat([follow_time_tbt[1:],pd.Series(np.nan)],ignore_index=True).to_list()\n",
    "all_evaluation[\"follow_ratio_this_trial\"]=follow_time_tbt\n",
    "all_evaluation.reset_index(drop=True, inplace=True)\n",
    "all_evaluation.loc[all_evaluation.index[all_evaluation['trial_id']==0].tolist(),'follow_ratio_previous_trial']=np.nan\n",
    "all_evaluation.loc[all_evaluation.index[all_evaluation['trial_id']==all_evaluation['trial_id'].max()].tolist(),'follow_ratio_next_trial']=np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['red', 'blue', 'green']\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5), tight_layout=True)\n",
    "ax1, ax2= axes.flatten()\n",
    "i=0\n",
    "for keys, grp in all_evaluation.groupby(['duration_ISI']):\n",
    "    print(keys)\n",
    "    this_color=colors[i]\n",
    "    ax1.scatter(grp['follow_ratio_previous_trial'],grp['follow_ratio_this_trial'],c=this_color)\n",
    "    ax2.scatter(grp['follow_ratio_this_trial'],grp['follow_ratio_next_trial'],c=this_color)\n",
    "    i=i+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_evaluation\n",
    "all_evaluation.to_csv(\"all_evaluation.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1D histogram to show the frequency of following response\n",
    "fig, axes = plt.subplots(\n",
    "        nrows=2, ncols=2, figsize=(9,10), tight_layout=True\n",
    "    )\n",
    "ax1, ax2, ax3, ax4 = axes.flatten()\n",
    "follow_time_aba=all_evaluation.groupby(['animal_id'])['num_follow_epochs'].sum()/all_evaluation.groupby(['animal_id'])['number_frames'].sum()\n",
    "simulated_time_aba=all_evaluation.groupby(['animal_id'])['num_chance_epochs'].sum()/all_evaluation.groupby(['animal_id'])['number_frames'].sum()\n",
    "simulated_time_tbt=all_evaluation['num_chance_epochs']/all_evaluation['number_frames']\n",
    "print(\"first 1/3 best of followers:\", np.quantile(follow_time_aba, 0.66))\n",
    "print(\"middle 1/3 best of followers:\", np.quantile(follow_time_aba, 0.33))\n",
    "print(\"first 1/3 best of follow epochs:\", np.quantile(follow_time_tbt, 0.66))\n",
    "print(\"middle 1/3 best of follow epochs:\", np.quantile(follow_time_tbt, 0.33))\n",
    "print(\"first 1/3 best of simulated followers:\", np.quantile(simulated_time_aba, 0.66))\n",
    "print(\"first 1/3 best of simulated epochs:\", np.quantile(simulated_time_tbt, 0.66))\n",
    "fair_follower_threshold=np.quantile(follow_time_aba, 0.33)\n",
    "good_follower_threshold=np.quantile(follow_time_aba, 0.66)\n",
    "ax1.hist(follow_time_aba,color='red')\n",
    "ax1.hist(simulated_time_aba,color='tab:gray',alpha=0.5)\n",
    "ax1.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,0.5),yticks=[0,20],ylim=(0,20),title='proportion of time following aba')\n",
    "ax2.hist(follow_time_tbt,color='red',bins=20)\n",
    "ax2.hist(simulated_time_tbt,color='tab:gray',alpha=0.5,bins=20)\n",
    "ax2.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,1),title='proportion of time following tbt')\n",
    "follower_of_interest=all_evaluation.groupby(['animal_id'])['num_follow_epochs'].sum()/all_evaluation.groupby(['animal_id'])['number_frames'].sum()<good_follower_threshold\n",
    "rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "ax3.hist(all_evaluation[rows_of_follower.values]['num_follow_epochs']/all_evaluation[rows_of_follower.values]['number_frames'],density=True)\n",
    "ax3.set(xticks=[0,0.25,0.5,0.75,1],xlim=(0,1),title='proportion of time from animals below threshold')\n",
    "follower_of_interest=all_evaluation.groupby(['animal_id'])['num_follow_epochs'].sum()/all_evaluation.groupby(['animal_id'])['number_frames'].sum()>good_follower_threshold\n",
    "rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "ax4.hist(all_evaluation[rows_of_follower.values]['num_follow_epochs']/all_evaluation[rows_of_follower.values]['number_frames'],density=True)\n",
    "ax4.set(xticks=[0,0.25,0.5,0.75,1],xlim=(0,1),title='proportion of time from animals above threshold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##plot the relationship between travel distance and proportion of following time\n",
    "# print(\"first 1/3 best of followers:\", np.quantile(follow_time_aba, 0.66))\n",
    "# print(\"middle 1/3 best of followers:\", np.quantile(follow_time_aba, 0.33))\n",
    "camera_fps=analysis_methods.get(\"camera_fps\")\n",
    "fair_follower_threshold=np.quantile(follow_time_aba, 0.33)\n",
    "good_follower_threshold=np.quantile(follow_time_aba, 0.66)\n",
    "glocust_list=[]\n",
    "exp_agent_list=[]\n",
    "fig, axes = plt.subplots(\n",
    "        nrows=2, ncols=2, figsize=(9,10), tight_layout=True\n",
    "    )\n",
    "ax1, ax2, ax3, ax4 = axes.flatten()\n",
    "ax1.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,0.5))\n",
    "ax2.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,1))\n",
    "for keys, this_data in all_evaluation.groupby(['animal_id']):\n",
    "    p_follow=this_data['num_follow_epochs'].sum()/this_data['number_frames'].sum()\n",
    "    if p_follow>good_follower_threshold:\n",
    "        color='k'\n",
    "    elif (p_follow>fair_follower_threshold) and (p_follow<good_follower_threshold):\n",
    "        color='k'\n",
    "    else:\n",
    "        color='k'\n",
    "    ax1.scatter(this_data['num_follow_epochs'].sum()/this_data['number_frames'].sum(),this_data['travel_distance'].sum()/camera_fps,c=color,s=8)\n",
    "    ax3.scatter(this_data['num_follow_epochs'].sum()/this_data['number_frames'].sum(),this_data['turning_distance'].sum()/camera_fps,c=color,s=8)\n",
    "    # if p_follow>good_follower_threshold:\n",
    "    #     pass\n",
    "    # elif (p_follow>fair_follower_threshold) and (p_follow<good_follower_threshold):\n",
    "    #     pass\n",
    "    # else:\n",
    "    #     continue\n",
    "\n",
    "    for this_object in this_data['object'].unique():\n",
    "        these_num_follow_epochs=this_data['num_follow_epochs'][this_data['object']==this_object]\n",
    "        these_num_frames=this_data['number_frames'][this_data['object']==this_object]\n",
    "        these_travel_distance=this_data['travel_distance'][this_data['object']==this_object]\n",
    "        these_turning_distance=this_data['turning_distance'][this_data['object']==this_object]\n",
    "        if this_object=='mov_glocust' or this_object== \"LeaderLocust\":\n",
    "            # ax2.scatter(these_num_follow_epochs/these_num_frames,these_travel_distance,c=color,s=8)\n",
    "            # ax4.scatter(these_num_follow_epochs/these_num_frames,these_turning_distance,c=color,s=8)\n",
    "            ax2.scatter(these_num_follow_epochs/these_num_frames,these_travel_distance,c='b',s=8)\n",
    "            ax4.scatter(these_num_follow_epochs/these_num_frames,these_turning_distance,c='b',s=8)\n",
    "            follow_time_ratio_tbt=these_num_follow_epochs/these_num_frames\n",
    "            glocust_list.append(follow_time_ratio_tbt.values)\n",
    "        else:\n",
    "            # ax2.scatter(these_num_follow_epochs/these_num_frames,these_travel_distance,edgecolors=color,marker='o', facecolors='none',s=8)\n",
    "            # ax4.scatter(these_num_follow_epochs/these_num_frames,these_turning_distance,edgecolors=color,marker='o', facecolors='none',s=8)\n",
    "            ax2.scatter(these_num_follow_epochs/these_num_frames,these_travel_distance,c='r',s=8)\n",
    "            ax4.scatter(these_num_follow_epochs/these_num_frames,these_turning_distance,c='r',s=8)\n",
    "            follow_time_ratio_tbt=these_num_follow_epochs/these_num_frames\n",
    "            exp_agent_list.append(follow_time_ratio_tbt.values)\n",
    "ax3.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,0.5))\n",
    "ax4.set(xticks=[0,0.25,0.5,0.75,1],xticklabels=(['0', '25', '50', '75', '100']),xlim=(0,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_evaluation['batchs'] = pd.cut(all_evaluation['animal_id'], bins=3, labels=['1st', '2nd', '3rd'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\neuroLaptop\\AppData\\Local\\Temp\\ipykernel_3292\\2344363941.py:19: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.\n",
      "  for keys, grp in all_evaluation.groupby(['batchs']):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAPdCAYAAAA3faEtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjpElEQVR4nO39f5RkdX0n/j+LaegCtrqXn93UOuL0ZL4qgpFgwg/zCWwahhgJyfbGUTEuRpJgMOgkIkhMssT9ZCYQBKIkRl0jBJYQEpvEk49RsPdIQoiKiHsEXTU0UbCYmTU76S50qoeB+/3jOg3NDDAzdHd133o8zqlTvG+9u+p1L3V486z3ve9bK4qiCAAAAJWwX7cLAAAAYP4IeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUSF+3C3i6J554Iq1WK41GI7VardvlALBEFEWRdrudZrOZ/far/m+UxkMAnm5Px8IlF/JarVZWrlzZ7TIAWKIeeuihvOAFL+h2GQvOeAjAM3musXDJhbxGo5GkLHxgYKDL1QCwVExPT2flypWz40TVGQ8BeLo9HQuXXMjbeUrKwMCAQQ2AXfTKqYvGQwCeyXONhdW/qAEAAKCHCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4APePv//7v8zM/8zNpNpup1Wr567/+6zmvF0WRyy67LM1mMwceeGBOO+203H///XP6zMzM5MILL8zhhx+egw8+OGeffXYefvjhOX22bt2aN73pTRkcHMzg4GDe9KY35d/+7d8WeO8AoCTkVUink0xOls8A7Op73/tefviHfzjXXnvtbl+/4oorctVVV+Xaa6/N3XffneHh4Zxxxhlpt9uzfdavX59bb701N998c+688848+uijOeuss/L444/P9jnnnHPy5S9/OZ/61KfyqU99Kl/+8pfzpje9acH3b7F1dnQyuXUynR0GHoClpFYURdHtIp5qeno6g4ODmZqaysDAQLfLWTYmJpKxsWR6OhkYSMbHk9HRblcFMH/me3yo1Wq59dZb83M/93NJylm8ZrOZ9evX55JLLklSztoNDQ3l8ssvz/nnn5+pqakcccQRueGGG/K6170uSdJqtbJy5cp88pOfzJlnnpmvfe1rOeaYY/K5z30uJ554YpLkc5/7XE4++eT87//9v/PiF794t/XMzMxkZmZmzv6uXLlyyY6HE5MTGbtlLNMz0xnoH8j4uvGMjhh4ABbSno6FZvIqoNMpA97OH5rb7bJtRg9gzz344IPZtGlT1q5dO7utv78/p556au66664kyT333JPHHntsTp9ms5ljjz12ts8//dM/ZXBwcDbgJclJJ52UwcHB2T67s3HjxtnTOwcHB7Ny5cr53sV509nRydgtY2nPlANPe6adsVvGzOgBLBFCXgW0WuUM3s452aIo261Wd+sCWE42bdqUJBkaGpqzfWhoaPa1TZs25YADDsghhxzyrH2OPPLIXd7/yCOPnO2zO5deemmmpqZmHw899NDz2p+F1Gq3Mj0znSLlwFOkyPTMdFptAw/AUtDX7QJ4/prN8hTNdrsMeLVa0miU2wHYO7VabU67KIpdtj3d0/vsrv9zvU9/f3/6+/v3struaDaaGegfSHumnSJFaqml0d9Is2HgAVgKzORVQL1eXoPXaJTtRqNs1+vdrQtgORkeHk6SXWbbtmzZMju7Nzw8nO3bt2fr1q3P2mfz5s27vP//+T//Z5dZwuWq3lfP+LrxNPrLgafR38j4uvHU+ww8AEuBkFcRo6PJ5s3JAw+UzxZdAdg7q1atyvDwcG6//fbZbdu3b88dd9yRU045JUlywgknZP/995/T55FHHsl999032+fkk0/O1NRUvvCFL8z2+fznP5+pqanZPlUwOjKazRdtzgNvfyCbL9ps0RWAJcTpmhVSrycjI92uAmDpevTRR/PP//zPs+0HH3wwX/7yl3PooYfmhS98YdavX58NGzZkzZo1WbNmTTZs2JCDDjoo55xzTpJkcHAw5513Xt75znfmsMMOy6GHHpqLLrooxx13XE4//fQkyUtf+tL81E/9VH75l385H/rQh5Ikv/Irv5KzzjrrGVfWXK7qffWMHGLgAVhqhDwAesYXv/jF/Mf/+B9n27/xG7+RJDn33HNz3XXX5eKLL862bdtywQUXZOvWrTnxxBNz2223pbHzfPgkV199dfr6+rJu3bps27Yto6Ojue6667JixYrZPv/jf/yPvP3tb59dhfPss89+xnvzAcB8c588AJaFXhsfem1/AXhu7pMHAADQg4Q8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfJ4UqeTTE6WzwAAwLIk5FGamEiGhpLVq8vniYluVwQAAOwDIY9y5m5sLGm3y3a7XbbN6AEAwLIj5JG0Wsn0dFIUZbsoynar1d26AACAvSbkkTSbycBAUquV7VqtbDeb3a0LAADYa0IeSb2ejI8njUbZbjTKdr3e3boAAIC91tftAlgiRkeTzZvLUzSbTQEPAACWKSGPJ9XrychIt6sAAACeB6drAgAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5LHmdTjI5WT4DAADPTshjSZuYSIaGktWry+eJiW5XBAAAS5uQx5LV6SRjY0m7Xbbb7bJtRg8AAJ6ZkMeS1Wol09NJUZTtoijbrVZ36wIAgKVMyGPJajaTgYGkVivbtVrZbja7WxcAACxlQh5LVr2ejI8njUbZbjTKdr3e3boAAGAp6+t2AfBsRkeTzZvLUzSbTQEPAACei5DHklevJyMj3a4CAACWB6drAgAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAT9Fut7N+/focffTROfDAA3PKKafk7rvvnn29KIpcdtllaTabOfDAA3Paaafl/vvvn/MeMzMzufDCC3P44Yfn4IMPztlnn52HH354sXcFgB4l5AHAU/zSL/1Sbr/99txwww35yle+krVr1+b000/Pd77znSTJFVdckauuuirXXntt7r777gwPD+eMM85Iu92efY/169fn1ltvzc0335w777wzjz76aM4666w8/vjj3dotAHrIXoW8HTt25Ld+67eyatWqHHjggRkZGcl73/vePPHEE7N99uQXTgBYirZt25aPf/zjueKKK/ITP/ET+aEf+qFcdtllWbVqVT74wQ+mKIpcc801ec973pOxsbEce+yxuf766/P9738/N910U5JkamoqH/3oR/O+970vp59+eo4//vjceOON+cpXvpLPfOYzz/jZMzMzmZ6envMAgH2xVyHv8ssvz5/8yZ/k2muvzde+9rVcccUV+YM/+IN84AMfmO2zJ79wAsBStGPHjjz++OOp1+tzth944IG588478+CDD2bTpk1Zu3bt7Gv9/f059dRTc9dddyVJ7rnnnjz22GNz+jSbzRx77LGzfXZn48aNGRwcnH2sXLlynvcOgF6xVyHvn/7pn/KzP/uzec1rXpMXvehF+fmf//msXbs2X/ziF5Nkj37hfDq/XAKwVDQajZx88sn5b//tv6XVauXxxx/PjTfemM9//vN55JFHsmnTpiTJ0NDQnL8bGhqafW3Tpk054IADcsghhzxjn9259NJLMzU1Nft46KGH5nnvAOgVexXyfvzHfzwTExP5xje+kST5X//rf+XOO+/MT//0TyfJHv3C+XR+uQRgKbnhhhtSFEX+w3/4D+nv78/73//+nHPOOVmxYsVsn1qtNudviqLYZdvTPVef/v7+DAwMzHkAwL7Yq5B3ySWX5A1veENe8pKXZP/998/xxx+f9evX5w1veEOS7NEvnE/nl0sAlpLVq1fnjjvuyKOPPpqHHnooX/jCF/LYY49l1apVGR4eTpJdxrQtW7bMjn3Dw8PZvn17tm7d+ox9AGAh7VXI+4u/+IvceOONuemmm/KlL30p119/fa688spcf/31c/rtzS+cfrkEYCk6+OCDc9RRR2Xr1q359Kc/nZ/92Z+dDXq33377bL/t27fnjjvuyCmnnJIkOeGEE7L//vvP6fPII4/kvvvum+0DAAupb286v+td78q73/3uvP71r0+SHHfccfnWt76VjRs35txzz53zC+dRRx01+3d+vQRgufj0pz+doijy4he/OP/8z/+cd73rXXnxi1+cX/zFX0ytVsv69euzYcOGrFmzJmvWrMmGDRty0EEH5ZxzzkmSDA4O5rzzzss73/nOHHbYYTn00ENz0UUX5bjjjsvpp5/e5b0DoBfsVcj7/ve/n/32mzv5t2LFitlbKDz1F87jjz8+yZO/cF5++eXzVDIALJypqalceumlefjhh3PooYfmP//n/5zf+73fy/77758kufjii7Nt27ZccMEF2bp1a0488cTcdtttaTQas+9x9dVXp6+vL+vWrcu2bdsyOjqa6667bs51fQCwUGpFURR72vnNb35zPvOZz+RDH/pQXvayl+Xee+/Nr/zKr+Qtb3nLbIi7/PLLs3HjxnzsYx+b/YXzs5/9bL7+9a/PGQCfyfT0dAYHBzM1NeXUTQBm9dr40Gv7C8Bz29OxYa9m8j7wgQ/kt3/7t3PBBRdky5YtaTabOf/88/M7v/M7s3325BdOAAAAFsZezeQtBr9cArA7vTY+9Nr+AvDc9nRs2KvVNQEAAFjahDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyGPJ6nSSycnyGQAA2DNCHkvSxEQyNJSsXl0+T0x0uyIAAFgehDyWnE4nGRtL2u2y3W6XbTN6AADw3IQ8lpxWK5meToqibBdF2W61ulsXAAAsB0IeS06zmQwMJLVa2a7Vynaz2d26AABgORDyWHLq9WR8PGk0ynajUbbr9e7WBQAAy0FftwuA3RkdTTZvLk/RbDYFPAAA2FNCHktWvZ6MjHS7CgAAWF6crgkAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQt1R1OsnkZPkMAACwh4S8pWhiIhkaSlavLp8nJrpdEQAAsEwIefugs6OTya2T6exYgFm2TicZG0va7bLdbpdtM3oAAMAeEPL20sTkRIauHMrq96/O0JVDmZic51m2ViuZnk6KomwXRdluteb3cwAAgEoS8vZCZ0cnY7eMpT1TzrK1Z9oZu2Vsfmf0ms1kYCCp1cp2rVa2m835+wwAAKCyhLy90Gq3Mj0znSLlLFuRItMz02m153GWrV5PxseTRqNsNxplu16fv88AAAAqq6/bBSwnzUYzA/0Dac+0U6RILbU0+htpNuZ5lm10NNm8uTxFs9kU8AAAgD1mJm8v1PvqGV83nkZ/OcvW6G9kfN146n17EML29pYI9XoyMiLgAQAAe8VM3l4aHRnN5os2p9Vupdlo7lnAm5goV8icni6vrxsfL2frAAAA5pmZvH1Q76tn5JCRPZ/Bc0sEAABgkQh5C80tEQAAgEUk5C00t0QAAAAWkZC30NwSAQAAWEQWXlkMbokAAAAsEiFvsey8JQIAAMACcromAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBwFPs2LEjv/Vbv5VVq1blwAMPzMjISN773vfmiSeemO1TFEUuu+yyNJvNHHjggTnttNNy//33z3mfmZmZXHjhhTn88MNz8MEH5+yzz87DDz+82LsDQA8S8gDgKS6//PL8yZ/8Sa699tp87WtfyxVXXJE/+IM/yAc+8IHZPldccUWuuuqqXHvttbn77rszPDycM844I+12e7bP+vXrc+utt+bmm2/OnXfemUcffTRnnXVWHn/88W7sFgA9pFYURdHtIp5qeno6g4ODmZqaysDAQLfLAWCJWKzx4ayzzsrQ0FA++tGPzm77z//5P+eggw7KDTfckKIo0mw2s379+lxyySVJylm7oaGhXH755Tn//PMzNTWVI444IjfccENe97rXJUlarVZWrlyZT37ykznzzDN3+dyZmZnMzMzM2d+VK1caDwGYtadjoZk8AHiKH//xH8/ExES+8Y1vJEn+1//6X7nzzjvz0z/900mSBx98MJs2bcratWtn/6a/vz+nnnpq7rrrriTJPffck8cee2xOn2azmWOPPXa2z9Nt3Lgxg4ODs4+VK1cu1C4CUHF93S4AAJaSSy65JFNTU3nJS16SFStW5PHHH8/v/d7v5Q1veEOSZNOmTUmSoaGhOX83NDSUb33rW7N9DjjggBxyyCG79Nn590936aWX5jd+4zdm2ztn8gBgbwl5APAUf/EXf5Ebb7wxN910U172spfly1/+ctavX59ms5lzzz13tl+tVpvzd0VR7LLt6Z6tT39/f/r7+5//DgDQ84Q8AHiKd73rXXn3u9+d17/+9UmS4447Lt/61reycePGnHvuuRkeHk5SztYdddRRs3+3ZcuW2dm94eHhbN++PVu3bp0zm7dly5accsopi7g3APQi1+QBwFN8//vfz377zR0eV6xYMXsLhVWrVmV4eDi333777Ovbt2/PHXfcMRvgTjjhhOy///5z+jzyyCO57777hDwAFpyZPAB4ip/5mZ/J7/3e7+WFL3xhXvayl+Xee+/NVVddlbe85S1JytM0169fnw0bNmTNmjVZs2ZNNmzYkIMOOijnnHNOkmRwcDDnnXde3vnOd+awww7LoYcemosuuijHHXdcTj/99G7uHgA9QMgDgKf4wAc+kN/+7d/OBRdckC1btqTZbOb888/P7/zO78z2ufjii7Nt27ZccMEF2bp1a0488cTcdtttaTQas32uvvrq9PX1Zd26ddm2bVtGR0dz3XXXZcWKFd3YLQB6iPvkAbAs9Nr40Gv7C8Bzc588AACAHiTkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkLqLOjk8mtk+ns6HS7FAAAoEcIeQtkYnIiQ1cOZfX7V2foyqFMTE50uyQAAKAHCHkLoLOjk7FbxtKeaSdJ2jPtjN0yZkYPAABYcELeAmi1W5memU6RIklSpMj0zHRa7VaXKwMAAKpOyFsAzUYzA/0DqaWWJKmlloH+gTQbzS5XBgAAVJ2QtwDqffWMrxtPo7+RJGn0NzK+bjz1vnqXKwMAAKqur9sFVNXoyGg2X7Q5rXYrzUZTwAMAABaFkLeA6n31jBwy0u0yAACAHuJ0TQAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKiQvQ553/nOd/ILv/ALOeyww3LQQQflFa94Re65557Z14uiyGWXXZZms5kDDzwwp512Wu6///55LRoAAIDd26uQt3Xr1rzqVa/K/vvvn7/7u7/LV7/61bzvfe/Lv//3/362zxVXXJGrrroq1157be6+++4MDw/njDPOSLvdnu/aAQAAeJq+vel8+eWXZ+XKlfnYxz42u+1FL3rR7D8XRZFrrrkm73nPezI2NpYkuf766zM0NJSbbrop559//vxUDQAAwG7t1UzeJz7xibzyla/Ma1/72hx55JE5/vjj85GPfGT29QcffDCbNm3K2rVrZ7f19/fn1FNPzV133bXb95yZmcn09PScBwAAAPtmr0Le5ORkPvjBD2bNmjX59Kc/nbe+9a15+9vfnj/7sz9LkmzatClJMjQ0NOfvhoaGZl97uo0bN2ZwcHD2sXLlyn3ZDwAAALKXIe+JJ57Ij/zIj2TDhg05/vjjc/755+eXf/mX88EPfnBOv1qtNqddFMUu23a69NJLMzU1Nft46KGH9nIXAAAA2GmvQt5RRx2VY445Zs62l770pfn2t7+dJBkeHk6SXWbttmzZssvs3k79/f0ZGBiY8wAAAGDf7FXIe9WrXpWvf/3rc7Z94xvfyNFHH50kWbVqVYaHh3P77bfPvr59+/bccccdOeWUU+ahXAAAAJ7NXq2u+eu//us55ZRTsmHDhqxbty5f+MIX8uEPfzgf/vCHk5Snaa5fvz4bNmzImjVrsmbNmmzYsCEHHXRQzjnnnAXZAQAAAJ60VyHvR3/0R3Prrbfm0ksvzXvf+96sWrUq11xzTd74xjfO9rn44ouzbdu2XHDBBdm6dWtOPPHE3HbbbWk0GvNePAAAAHPViqIoul3EU01PT2dwcDBTU1OuzwNgVq+ND722vwA8tz0dG/bqmjwAAACWNiEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8pawzo5OJrdOprOj0+1SAACAZULIW6ImJicydOVQVr9/dYauHMrE5ES3SwIAAJYBIW8J6uzoZOyWsbRn2kmS9kw7Y7eMmdEDAACek5C3BLXarUzPTKdIkSQpUmR6ZjqtdqvLlQEAAEudkLcENRvNDPQPpJZakqSWWgb6B9JsNLtcGQDQDZ1OMjlZPgM8FyFvCar31TO+bjyN/kaSpNHfyPi68dT76l2uDKD6XvSiF6VWq+3yeNvb3pYkKYoil112WZrNZg488MCcdtppuf/+++e8x8zMTC688MIcfvjhOfjgg3P22Wfn4Ycf7sbuUAETE8nQULJ6dfk84TJ94DkIeUvU6MhoNl+0OQ+8/YFsvmhzRkdGu10SQE+4++6788gjj8w+br/99iTJa1/72iTJFVdckauuuirXXntt7r777gwPD+eMM85Iu92efY/169fn1ltvzc0335w777wzjz76aM4666w8/vjjXdknlq9OJxkbS3Z+vdrtsm1GD3g2taIoim4X8VTT09MZHBzM1NRUBgYGul0OAEtEt8aH9evX52//9m/zzW9+M0nSbDazfv36XHLJJUnKWbuhoaFcfvnlOf/88zM1NZUjjjgiN9xwQ173utclSVqtVlauXJlPfvKTOfPMM3f7OTMzM5mZmZltT09PZ+XKlcbDHjc5Wc7gPd0DDyQjI4tfD9BdezoWmskDgGewffv23HjjjXnLW96SWq2WBx98MJs2bcratWtn+/T39+fUU0/NXXfdlSS555578thjj83p02w2c+yxx8722Z2NGzdmcHBw9rFy5cqF2zGWjWYzGRhIauVl+qnVynbTZfrAsxDyqsDV2AAL4q//+q/zb//2b3nzm9+cJNm0aVOSZGhoaE6/oaGh2dc2bdqUAw44IIcccsgz9tmdSy+9NFNTU7OPhx56aB73hOWqXk/Gx5NGeZl+Go2yXXeZPvAs+rpdAM/TxER5cv70dPnT3vh4Mur6PYD58NGPfjSvfvWr03zatElt57TKDxRFscu2p3uuPv39/env79/3Yqms0dFk8+ak1Spn8AQ84LmYyVvOXI0NsGC+9a1v5TOf+Ux+6Zd+aXbb8PBwkuwyI7dly5bZ2b3h4eFs3749W7dufcY+sLfq9fIaPAEP2BNC3nLWapUzeDvXzimKst1y03SA5+tjH/tYjjzyyLzmNa+Z3bZq1aoMDw/PrriZlNft3XHHHTnllFOSJCeccEL233//OX0eeeSR3HfffbN9AGAhOV1zOdt5NXa7XQa8Wq08Wd/V2ADPyxNPPJGPfexjOffcc9PX9+RQWavVsn79+mzYsCFr1qzJmjVrsmHDhhx00EE555xzkiSDg4M577zz8s53vjOHHXZYDj300Fx00UU57rjjcvrpp3drlwDoIULecrbzauyd1+S5GhtgXnzmM5/Jt7/97bzlLW/Z5bWLL74427ZtywUXXJCtW7fmxBNPzG233ZbGzpUxklx99dXp6+vLunXrsm3btoyOjua6667LihUrFnM3AOhR7pNXBZ2Oq7GByuu18aHX9heA57anY4OZvCrYeTU2AADQ8yy8AgAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AwE6dTjI5WT4DLFNCHgBAkkxMJENDyerV5fPERLcrAtgnQh6wV/zIDVRSp5OMjSXtdtlut8u2/9gBy5CQB+wxP3IDldVqJdPTSVGU7aIo261Wd+sC2AdCHrBH/MgNVFqzmQwMJLVa2a7Vynaz2d26APaBkAfsET9yA5VWryfj40mjUbYbjbJdr3e3LoB90NftAoDlYeeP3O12GfBqtfL/gfzIDVTG6GiyeXP561WzKeABy5aZPGCP+JEb6An1ejIy4j9uwLJmJg/YY37kBgBY+oQ8YK/s/JEbAIClyemaAAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQB1dPpJJOT5TMAQI8R8oBqmZhIhoaS1avL54mJblcEALCohDygOjqdZGwsabfLdrtdts3oAQA9RMgDqqPVSqank6Io20VRtlut7tYFALCIhDygOprNZGAgqdXKdq1WtpvN7tYFALCIhDygOur1ZHw8aTTKdqNRtuv17tYFALCI+rpdAMC8Gh1NNm8uT9FsNgU8AKDnCHlA9dTrychIt6sAAOgKp2sCAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeXdPpJJOT5TMAADA/hDy6YmIiGRpKVq8unycmul0RAMCe8UM1S52Qx6LrdJKxsaTdLtvtdtn2H0oAYKnzQzXLgZDHomu1kunppCjKdlGU7Varu3UBADwbP1SzXPR1uwB6xNRUMj6e/N//m/9wwikZaJyU9qO1FEVSqyWNRtJsdrtIAIBntvOH6p2e+kP1yEj36oKnE/JYeH/7t8kb3pB8//tJvZ7+738/D/zYT+dlX/urbGkfmEajzH/1ercLBQB4Zs1mMjBQzuD5oZqlzOmaLKzNm5N165Kf/Mnk4YfL/yp+/OM5/L7P5ju/+Ft54IGyy+hotwsFAHh29Xr5w3SjUbb9UM1SZSaPhXXTTckTTyQf+1hy6KHltrGx5POfT9+HPpSRq/8g2c9vDQDA8jA6Wv5A3WqVM3gCHkuR/7tmYW3alBx11JMBb6djjimv05uZ6U5dAAD7qF4vr8ET8FiqhDwW1gknJP/yL8mXvvTktqIoz2146Uv91xEAAOaZ0zVZWD/3c8nLXpacdVbym7+ZvOhFyY03Jp/4RPI//kd5xTIAADBvhDwW1gEHJJ/5TPK2tyXveEd5fd4LX5j86Z8m55zT7eoAAKByhDwW3vBw8vGPl9fgTU+XVymvWNHtqgAAoJKEPBbP4GD5AAAAFoyFVwAAACpEyAMAAKgQIa9qOp1kcrJ8BgAAeo6QVyUTE8nQULJ6dfk8MdHtigAAgEUm5FVFp5OMjSXtdtlut8u2GT0AAOgpQl5VtFrl7QmKomwXRdlutbpbFwDsK5cgAOwTIa8qms1kYCCp1cp2rVa2m83u1gUA+8IlCAD7TMirino9GR9PGo2y3WiU7Xq9u3UBwN5yCQLA8+Jm6FUyOpps3lyeotlsCngALE87L0HY6amXIIyMdK8ugGXCTF7V1OvlACjgAbBcuQQB4HkR8gCApcUlCADPi9M1AYClxyUIAPvMTN4i6OzoZHLrZDo7XDAOAHvMJQgA+0TIW2ATkxMZunIoq9+/OkNXDmVi0hLQAADAwhHyFlBnRydjt4ylPVMuAd2eaWfsljEzegAAwIIR8hZQq93K9Mx0ihRJkiJFpmem02q3ulwZAABQVULeAmo2mhnoH0gt5RLQtdQy0D+QZsMS0AAAwMIQ8hZQva+e8XXjafSXS0A3+hsZXzeeep8LyAEAgIXhFgoLbHRkNJsv2pxWu5VmoyngAQAAC0rIWwT1vnpGDhnpdhkAAEAPcLomAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkA8DTf+c538gu/8As57LDDctBBB+UVr3hF7rnnntnXi6LIZZddlmazmQMPPDCnnXZa7r///jnvMTMzkwsvvDCHH354Dj744Jx99tl5+OGHF3tXAOhBQh4APMXWrVvzqle9Kvvvv3/+7u/+Ll/96lfzvve9L//+3//72T5XXHFFrrrqqlx77bW5++67Mzw8nDPOOCPtdnu2z/r163Prrbfm5ptvzp133plHH300Z511Vh5//PEu7BUAvaRWFEXR7SKeanp6OoODg5mamsrAwEC3ywFgiVis8eHd7353/vEf/zH/8A//sNvXi6JIs9nM+vXrc8kllyQpZ+2GhoZy+eWX5/zzz8/U1FSOOOKI3HDDDXnd616XJGm1Wlm5cmU++clP5swzz9zlfWdmZjIzMzPbnp6ezsqVK42HAMza07HQTB4APMUnPvGJvPKVr8xrX/vaHHnkkTn++OPzkY98ZPb1Bx98MJs2bcratWtnt/X39+fUU0/NXXfdlSS555578thjj83p02w2c+yxx872ebqNGzdmcHBw9rFy5coF2kMAqk7I41kVRXLXXcmVVyYf/Wjyb//W7YoAFtbk5GQ++MEPZs2aNfn0pz+dt771rXn729+eP/uzP0uSbNq0KUkyNDQ05++GhoZmX9u0aVMOOOCAHHLIIc/Y5+kuvfTSTE1NzT4eeuih+d41AHpEX7cLYOn6/veTn//55O/+Ljn44GTbtuQd70j+/M+Tn/mZblcHsDCeeOKJvPKVr8yGDRuSJMcff3zuv//+fPCDH8x/+S//ZbZfrVab83dFUeyy7emerU9/f3/6+/ufZ/UAYCaPZ/Hbv5189rPJxz+eTE8n3/lOcsYZybp1yTP8EA2w7B111FE55phj5mx76Utfmm9/+9tJkuHh4STZZUZuy5Yts7N7w8PD2b59e7Zu3fqMfQBgoQh57NYTTyR/+qfJhRcmY2PJfvslw8PltiS56abu1gewUF71qlfl61//+pxt3/jGN3L00UcnSVatWpXh4eHcfvvts69v3749d9xxR0455ZQkyQknnJD9999/Tp9HHnkk991332wfAFgoTtdkt7ZvL6+/e+lL524/5JDkqKPM5AHV9eu//us55ZRTsmHDhqxbty5f+MIX8uEPfzgf/vCHk5Snaa5fvz4bNmzImjVrsmbNmmzYsCEHHXRQzjnnnCTJ4OBgzjvvvLzzne/MYYcdlkMPPTQXXXRRjjvuuJx++und3D0AeoCQx2719yfHHJOMjyfnnpvsvITkS19KHnwweeUru1sfwEL50R/90dx666259NJL8973vjerVq3KNddckze+8Y2zfS6++OJs27YtF1xwQbZu3ZoTTzwxt912WxqNxmyfq6++On19fVm3bl22bduW0dHRXHfddVmxYkU3dguAHuI+eTyjm29O3vCG5LWvTX7hF5JvfSvZsCE57LAy7B1wQLcrBHpJr40PVdrfzo5OWu1Wmo1m6n31bpcDsGzt6dhgJo9n9PrXl6dt/tf/mvzlXyYrViT/6T8lH/iAgAfAnpmYnMjYLWOZnpnOQP9AxteNZ3RktNtlAVSahVd4Vv/lvyQPPJA89FDyr/9ahr0fLCwHAM+qs6OTsVvG0p5pJ0naM+2M3TKWzo5OlysDqDYhj+e0337JC16QDA52uxIAlpNWu5XpmekUKa8MKVJkemY6rXary5UBVJuQBwAsiGajmYH+gdRSrt5VSy0D/QNpNppdrgyg2oQ8AGBB1PvqGV83nkZ/uepoo7+R8XXjFl8BWGAWXgEAFszoyGg2X7TZ6poAi0jIAwAWVL2vnpFDRrpdBkDPcLomAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIe9IhOJ5mcLJ8BAKguIQ96wMREMjSUrF5dPk9MdLsiAAAWipAHFdfpJGNjSbtdttvtsm1GDwCgmoQ8qLhWK5meToqibBdF2W61ulsXAAALQ8iDims2k4GBpFYr27Va2W42u1sXAAALQ8iDiqvXk/HxpNEo241G2a7Xu1sXAAALo6/bBQALb3Q02by5PEWz2RTwAACqTMiDHlGvJyMj3a4CAICF5nRNAACAChHyAFhYnU4yOem+HQCwSIQ8ABbOxEQyNJSsXl0+T0x0uyIAqDwhD4CF0ekkY2NJu1222+2ybUYPABbU8wp5GzduTK1Wy/r162e3FUWRyy67LM1mMwceeGBOO+203H///c+3TgCWm1YrmZ5OiqJsF0XZbrW6WxcAVNw+h7y77747H/7wh/Pyl798zvYrrrgiV111Va699trcfffdGR4ezhlnnJH2zl9yAegNzWYyMJDUamW7VivbzWZ36wKAitunkPfoo4/mjW98Yz7ykY/kkEMOmd1eFEWuueaavOc978nY2FiOPfbYXH/99fn+97+fm266abfvNTMzk+np6TkPACqgXk/Gx5NGo2w3GmXbjRoBYEHtU8h729velte85jU5/fTT52x/8MEHs2nTpqxdu3Z2W39/f0499dTcddddu32vjRs3ZnBwcPaxcuXKfSkJgKVodDTZvDl54IHyeXS02xUBQOXtdci7+eab86UvfSkbN27c5bVNmzYlSYaGhuZsHxoamn3t6S699NJMTU3NPh566KG9LQnLkwNLWb2ejIyYwQOARbJXIe+hhx7KO97xjtx4442pP8tgXdt5/cUPFEWxy7ad+vv7MzAwMOfBXrA8OQAA8BR7FfLuueeebNmyJSeccEL6+vrS19eXO+64I+9///vT19c3O4P39Fm7LVu27DK7xzywPDkAAPA0exXyRkdH85WvfCVf/vKXZx+vfOUr88Y3vjFf/vKXMzIykuHh4dx+++2zf7N9+/bccccdOeWUU+a9+J5neXIAAOBp+vamc6PRyLHHHjtn28EHH5zDDjtsdvv69euzYcOGrFmzJmvWrMmGDRty0EEH5Zxzzpm/qintXJ683S4DXq1Wrl5neXIAgF10OuVv4c2my4Sptud1M/Tdufjii7N+/fpccMEFeeUrX5nvfOc7ue2229LYuYQ288fy5AAAe8QyBvSSWlHsPNdvaZiens7g4GCmpqYswrKn/CwF9IBeGx96bX9hIXU6ZbB7+slPmzf7XyeWlz0dG+Z9Jo8usDw5AMAzsowBvUbIAwCg0nYuY7Dzjl61Wtm2jAFVJeQBAFBpljGg1+zV6poAALAcjY6W1+BZxoBeIOQBANATdi5jAFXndE0AAIAKEfIAAAAqRMgDAACoECGPZ9TZ0cnk1sl0dnS6XQoAALCHhDx2a2JyIkNXDmX1+1dn6MqhTExOdLskAABgDwh57KKzo5OxW8bSnmknSdoz7YzdMmZGDwAAlgEhj1202q1Mz0ynSJEkKVJkemY6rXary5UBAADPRchjF81GMwP9A6mlliSppZaB/oE0G80uVwYAADwXIY9d1PvqGV83nkZ/I0nS6G9kfN146n31LlcGAAA8l75uF8DSNDoyms0XbU6r3Uqz0RTwAABgmRDyeEb1vnpGDhnpdhkAAMBecLomAABAhQh5AABQRZ1OMjlZPtNThDwAAKiaiYlkaChZvbp8npjodkUsIiEPAACqpNNJxsaSdrtst9tl24xezxDylimz7wAA7FarlUxPJ0VRtouibLda3a2LRSPkLUNm3wEAeEbNZjIwkNRqZbtWK9vNZnfrYtEIecuM2XcAAJ5VvZ6MjyeNRtluNMp23X2Pe4X75C0zO2ffd3rq7PuIW9oBAJAko6PJ5s3l/yQ2mwJejzGTt8yYfQcAYI/U6+UsgIDXc4S8ZcbsOwAA8GycrrkMmX0HgGfR6RgkgZ5mJm+ZMvsOALthCWoAIQ8AqAhLUAMkEfIAgKpwA2iAJEIeAFAVlqAGSCLkAQBVYQlqgCRW1wQAlrDOjk5a7VaajWbqfXsQ1ixBDWAmDwBYmiYmJzJ05VBWv391hq4cysTkHq6UaQlqoMcJeQDAktPZ0cnYLWNpz5QrZbZn2hm7ZSydHVbKBHguQh4AsOS02q1Mz0ynSLlSZpEi0zPTabWtlAnwXIQ8AGDJaTaaGegfSC3lSpm11DLQP5Bmw0qZAM9FyAMAlpx6Xz3j68bT6C9Xymz0NzK+bnzPFl8B6HFW1wQAlqTRkdFsvmjz3q2uCYCQBwAsXfW+ekYOGel2GQDLitM1AQAAKkTIAwAAqBAhDwAAoEKEPAAAgAoR8gAAACpEyJsnRVHkT+/905z430/MC656Qc666az8/bf+vttlAQAAPUbImyfvvO2dOe8T52Xo4KG8+RVvznfa38l/vP4/5q//9193uzQAAKCHCHnzYHLrZK753DW54vQr8ok3fCL/70/+v/niL38xP/VDP5V33f6uPFE80e0SAQCAHiHkzYPPTH4mtVotb/uxt81uW7Hfivzaj/5a/vn//nMmt052sToAAKCXCHnzoH9Ff54onsj3tn9vzvb29vbs6wAAAItByJsHZ/3/zkq9r553f+bd2fHEjiTJv37/X7PhHzbkpBeclJWDK7tcIQAA0Cv6ul1AFRx20GH545/+45z3ifPy6Qc+nWOPPDZ3fvvO9Pf1Z+K/THS7PAAAoIeYyZsnv3j8L+ZL538pP/eSn8uB+x+Yd53yrtz3q/flFcOv6HZpAABADzGTN49eMfyKXPvT13a7DAAAoIeZyQMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AHiKyy67LLVabc5jeHh49vWiKHLZZZel2WzmwAMPzGmnnZb7779/znvMzMzkwgsvzOGHH56DDz44Z599dh5++OHF3hUAepSQBwBP87KXvSyPPPLI7OMrX/nK7GtXXHFFrrrqqlx77bW5++67Mzw8nDPOOCPtdnu2z/r163Prrbfm5ptvzp133plHH300Z511Vh5//PFu7A4APcbN0AHgafr6+ubM3u1UFEWuueaavOc978nY2FiS5Prrr8/Q0FBuuummnH/++ZmamspHP/rR3HDDDTn99NOTJDfeeGNWrlyZz3zmMznzzDN3+5kzMzOZmZmZbU9PTy/AngHQC8zkzZdOJ5mcLJ8BWNa++c1vptlsZtWqVXn961+fycnJJMmDDz6YTZs2Ze3atbN9+/v7c+qpp+auu+5Kktxzzz157LHH5vRpNps59thjZ/vszsaNGzM4ODj7WLly5QLtHQBVJ+TNh4mJZGgoWb26fJ6Y6HZFAOyjE088MX/2Z3+WT3/60/nIRz6STZs25ZRTTsm//uu/ZtOmTUmSoaGhOX8zNDQ0+9qmTZtywAEH5JBDDnnGPrtz6aWXZmpqavbx0EMPzfOeAdArnK75fHU6ydhYsvNajHa7bG/enNTr3a0NgL326le/evafjzvuuJx88slZvXp1rr/++px00klJklqtNudviqLYZdvTPVef/v7+9Pf3P4/Kl6hOJ2m1kmbTuAiwSMzkPV+tVjI9nRRF2S6Kst1qdbcuAObFwQcfnOOOOy7f/OY3Z6/Te/qM3JYtW2Zn94aHh7N9+/Zs3br1Gfv0DGe6AHSFkPd8NZvJwECy89fZWq1sN5vdrQuAeTEzM5Ovfe1rOeqoo7Jq1aoMDw/n9ttvn319+/btueOOO3LKKackSU444YTsv//+c/o88sgjue+++2b79IRnOtPFtesAC07Ie77q9WR8PGk0ynajUbaf4ZSUzo5OJrdOprPDIAewFF100UW544478uCDD+bzn/98fv7nfz7T09M599xzU6vVsn79+mzYsCG33npr7rvvvrz5zW/OQQcdlHPOOSdJMjg4mPPOOy/vfOc7MzExkXvvvTe/8Au/kOOOO252tc2e4EwXgK5xTd58GB0tr8F7jmsOJiYnMnbLWKZnpjPQP5DxdeMZHRld5GIBeDYPP/xw3vCGN+S73/1ujjjiiJx00kn53Oc+l6OPPjpJcvHFF2fbtm254IILsnXr1px44om57bbb0tj5Y1+Sq6++On19fVm3bl22bduW0dHRXHfddVmxYkW3dmvx7TzTpd0uA16tVv4Q6kwXgAVXK4qdP7EtDdPT0xkcHMzU1FQGBga6Xc686ezoZOjKobRn2ilSpJZaGv2NbL5oc+p9LkQHeC5VHR+eSSX2d2KiPEVzeroMfOPj5Q+jAOyTPR0bzOQtkla7lemZJ29sW6TI9Mx0Wu1WRg4Z6WJlALBA9vBMFwDml2vyFkmz0cxA/0BqKRdoqaWWgf6BNBtOWwGgwur1ZGREwANYRELeIqn31TO+bjyN/vKajUZ/I+Prxp2qCQAAzCunay6i0ZHRbL5oc1rtVpqNpoAHAADMOyFvkdX76q7BAwAAFozTNQEAACpEyAMAKqPTSSYny+cl59vfTt797uTVr07OOy/5/Oe7XRFQUUIeAFAJExPJ0FCyenX5PDHR7Yqe4u67k+OOSz70oaS/P/nsZ5OTTkr+5E+6XRlQQUIeALDsdTrlfdfb7bLdbpftJTOj96u/mqxZk/zLvyR//dfJN7+ZvPWtyTvekXz3u92uDqgYIQ8AWPZarWR6OimKsl0UZbvV6m5dSZJvfSu5557k0kuTwcFy2377Jf/tvyXbtyd/+7fdrQ+oHCEPAFj2ms1kYCCp1cp2rVa2m83u1pUkeeyx8vnAA+dur9fLQne+DjBPhDwAYNmr15Px8aTRKNuNRtmuL4Vb0o6MlKdq/uEfJjt2PLn9mmvKkHfmmV0rDagm98kDACphdDTZvLk8RbPZXCIBLylPzbz66uTnfi55+cuTn/7p5MtfLleGefe7kxe+sNsVAhVjJg8AqIx6vZw4WzIBb6fXvCb5h39Ijjkm+Zu/KWf0brop2bCh25UBFWQmDwBgMZx0UvJXf9XtKoAeYCYPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhD56u00kmJ8tnAABYZoQ8eKqJiWRoKFm9unyemOh2RQDLylRnKn/4uT/Mur9cl7f+7VvzTw/9U7dLAug5Qt4iM0m0hHU6ydhY0m6X7Xa7bPuXBbBHHpp6KMd/6Pi86/Z3Zcv3tuS2B27LKX96Sjb+w8ZulwbQU4S8RWSSaIlrtZLp6aQoynZRlO1Wq7t1ASwT77r9Xdn++PZ8/de+ns+++bP557f/cy798Uvzm//zN/P173692+UB9Awhb5GYJFoGms1kYCCp1cp2rVa2m83u1gWwDDz2+GMZ/9p43nHiO7LqkFVJkv1q++V3Tv2dDPQP5C+/+pddrhCgdwh5i8Qk0TJQryfj40mjUbYbjbJdr3e3LoBl4PHi8Tz2xGMZrA/O2d6/oj8H7X9Qtj22rUuVAfQeIW+RmCRaJkZHk82bkwceKJ9HR7tdEcCyUO+r5yeO/ol85EsfycyOmdntH//ax7Pp0U1Zu3ptF6sD6C193S6gV+ycJBobK2fwTBItYfV6MjLS7SoAlp0NP7kho382mh/+kx/Oa495bf5l6l9y830352df/LP5iaN/otvlAfQMM3mLyCQRAFX2qhe+Kv/4ln/My458Wf77vf89X3rkS/m9n/y93PLaW1LbeSoLAAvOTN4iM0kEQJWd0DwhH1/38W6XAdDTzOQBAABUiJAHAABQIUIeAABAhQh5AABLVKeTTE6Wz7314cDzIeQBACxBExPJ0FCyenX5PDHRKx8OPF9CHgDAEtPplPfWbbfLdrtdthdlUq2rHw7MByEPAGAPbNmSfPvbSVHM7/vu7qzIViuZnn7ys4qibLda8/vZu9XVDwfmg5AHAPAs/vf/Tn7yJ8uzFo8+OnnpS5NPfGJ+3vuZzopsNpOBgWTnPeRrtbLdbM7P5z6rrn44MB+EPACgGhZgoZDvfjc57bTkkUeS669P/uZvklWrkv/0n5LPfvb5vfeznRVZryfj40mjUb7WaJTtev35feYe6eqHA/Ohr9sFAAA8bxMTZUKani5nncbHk9HR5/22f/qnydRUcu+9yVFHldte85rkx34s+f3fLwPgvtp5VuROTz0rcmSkLH/z5rLdbC5yxurqhwPPl5k8AGB5W8CFQu65JznppCcDXpKsWJH87M8mX/rS83vvPTkrsl4vA19XMlZXPxx4PoQ8AGB5W8CFQo46Kvn615PHHpu7/b775ga/feGsSGChCHkAwPK2gAuFnHdesmlT8qu/mvzrvyYzM8kf/3HyV3+V/MqvPO+3nz0r8oEHyud5OMMUQMgDAJa5BZwSO+645L//9+SGG5Ijjyyz49velvzSL5XBbz44KxKYbxZeAQDmVWdHJ612K81GM/W+RUouC7hQyFvekvzMz5Qra37/+8kZZ5S3UQBYqoQ8AGDeTExOZOyWsUzPTGegfyDj68YzOrJI5yDunBJbAEccUc7eASwHTtcEAOZFZ0cnY7eMpT1TrnLZnmln7JaxdHbM333rAHhuQh4AMC9a7VamZ6ZTpFzlskiR6ZnptNrPf5VLAPackAcAzItmo5mB/oHUUq5yWUstA/0DaTae/yqXAOw5IQ8AmBf1vnrG142n0V+uctnob2R83fjiLb4CQBILr/SeTmdBVh4DgCQZHRnN5os2L/7qmgDMMpPXSyYmkqGhZPXq8nliotsVAVBB9b56Rg4ZEfAAukTI6xWdTjI2lrTLFc/SbpftjhXPAACgSoS8XtFqJdPTSVGueJaiKNstK54BsPA6OzqZ3DrpdgoAi0DI6xXNZjIwkNTKFc9Sq5XtphXPAFhYE5MTGbpyKKvfvzpDVw5lYtLlAgALScjrFfV6Mj6eNMoVz9JolG2LrwCwgNwgHWDxWV2zl4yOJps3W10TgEWz8wbpOz31Bukjh4x0sTKA6jKT12vq9WRkRMADYFHs9gbp+/+7NPc/tMuVAVSXkAcALJhdbpA+U2T8vz+a+n842q18ABaI0zUBgAU1OjKazb/2rbSOWZnmI4+mviNJ7Qe38tm82dklAPPMTN5C6nSSyUn3ogOg59W3/N+MPPSDgJe4lQ/AAhLyFsrERDI0lKxeXT47JQWAXuZWPgCLRshbCJ1OeQpKu1wuOu0fnJJiRg+AXuVWPgCLxjV5C6HVKk9B2empp6SMWC4agB7lVj4Ai8JM3kJwSgoA7J5b+QAsOCFvITglBQAA6BKnay4Up6QAAABdIOQtpJ2npAAAACwSp2sCAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAAAspE4nmZwsnxeBkAcAALBQJiaSoaFk9eryeWJiwT9SyAMAAFgInU4yNpa022W73S7bCzyjJ+QBAAAshFYrmZ5OiqJsF0XZbrUW9GOFPAAAgIXQbCYDA0mtVrZrtbLdbC7oxwp5AMCytM/rGCzyAghAD6vXk/HxpNEo241G2a7XF/RjhTwAYNnZ53UMurAAAtDjRkeTzZuTBx4on0dHF/wjhTwAYFnZ53UMurQAAkDq9WRkZMFn8HYS8gCAZWWf1zHo0gIIAItNyAMAlpV9XsegSwsgACw2IQ8AWFb2eR2DLi2AALDY+rpdAADA3tq5jkGrVU7E7XFO2+c/BFg+hDwAYFnauY7B4v0hwPLgdE0AAIAKEfIAAAAqRMgDAJ5Rp5NMTrqVHMByIuQBALs1MZEMDSWrV5fPExPdrgiAPSHkzbPOjk4mt06ms8NPngAsX51OMjaWtNtlu90u22b0AJY+IW8eTUxOZOjKoax+/+oMXTmUiUk/eQKwPLVayfR0UhRluyjKdqvV3boAeG5C3jzp7Ohk7JaxtGfKnzzbM+2M3TJmRg+AZanZTAYGklqtbNdqZbvZ7G5dADw3IW+etNqtTM9Mp0j5k2eRItMz02m1/eQJwPJTryfj40mjUbYbjbLt3uEAS5+boc+TZqOZgf6BtGfaKVKklloa/Y00G37yBGDp6+zopNVupdlopt5XJrnR0WTz5vIUzWZTwANYLszkzZN6Xz3j68bT6C9/8mz0NzK+bnx2oASAperZrimv15OREQEPYDkR8vbCc62cOToyms0Xbc4Db38gmy/anNGR0UWuEID5tHHjxtRqtaxfv352W1EUueyyy9JsNnPggQfmtNNOy/333z/n72ZmZnLhhRfm8MMPz8EHH5yzzz47Dz/88CJXv2dcUw5QPULeHtrTlTPrffWMHDJiBg9gmbv77rvz4Q9/OC9/+cvnbL/iiity1VVX5dprr83dd9+d4eHhnHHGGWnvvNdAkvXr1+fWW2/NzTffnDvvvDOPPvpozjrrrDz++OOLvRvPaW+vKXerIIClT8jbA37lBOgtjz76aN74xjfmIx/5SA455JDZ7UVR5Jprrsl73vOejI2N5dhjj83111+f73//+7npppuSJFNTU/noRz+a973vfTn99NNz/PHH58Ybb8xXvvKVfOYzn+nWLj2jndeU11Iuo1lLLQP9A7u9ptytggCWByFvD1g5E6C3vO1tb8trXvOanH766XO2P/jgg9m0aVPWrl07u62/vz+nnnpq7rrrriTJPffck8cee2xOn2azmWOPPXa2z+7MzMxkenp6zmMx7Ok15X7wBFg+rK65B6ycCdA7br755nzpS1/K3XffvctrmzZtSpIMDQ3N2T40NJRvfetbs30OOOCAOTOAO/vs/Pvd2bhxY373d3/3+Za/T3ZeU/701TWfaucPnjs99QfPkUNGFrNcAJ7DXs3kbdy4MT/6oz+aRqORI488Mj/3cz+Xr3/963P67MkF6cuNlTMBesNDDz2Ud7zjHbnxxhtTf5blJGs77xD+A0VR7LLt6Z6rz6WXXpqpqanZx0MPPbR3xT9Pz3VN+d6c1glAd+1VyLvjjjvytre9LZ/73Ody++23Z8eOHVm7dm2+973vzfbZkwvSlyMrZwJU3z333JMtW7bkhBNOSF9fX/r6+nLHHXfk/e9/f/r6+mZn8J4+I7dly5bZ14aHh7N9+/Zs3br1GfvsTn9/fwYGBuY8lhI/eAIsH7WiKIp9/eP/83/+T4488sjccccd+Ymf+IkURZFms5n169fnkksuSVJeYzA0NJTLL788559//nO+5/T0dAYHBzM1NbXkBjgAumcxxod2uz172uVOv/iLv5iXvOQlueSSS/Kyl70szWYzv/7rv56LL744SbJ9+/YceeSRs+Pc1NRUjjjiiNx4441Zt25dkuSRRx7JC17wgnzyk5/MmWeeuUe1LNXxcHc3TQdgcezp2PC8rsmbmppKkhx66KFJnvuC9N2FvJmZmczMzMwpHAC6odFo5Nhjj52z7eCDD85hhx02u339+vXZsGFD1qxZkzVr1mTDhg056KCDcs455yRJBgcHc9555+Wd73xnDjvssBx66KG56KKLctxxx+2ykMtytPO0TgCWrn0OeUVR5Dd+4zfy4z/+47MD355ckP503bzQHAD21sUXX5xt27blggsuyNatW3PiiSfmtttuS6PRmO1z9dVXp6+vL+vWrcu2bdsyOjqa6667LitWrOhi5QD0in0+XfNtb3tb/r//7//LnXfemRe84AVJkrvuuiuvetWr0mq1ctRRR832/eVf/uU89NBD+dSnPrXL++xuJm/lypVL7vQUALprqZ6+uFB6bX8BeG4LerrmhRdemE984hP5+7//+9mAl5QXmyfljN5TQ96zXWze39+f/v7+fSkDAACAp9mr1TWLosiv/dqvZXx8PP/zf/7PrFq1as7rq1atyvDwcG6//fbZbdu3b88dd9yRU045ZX4qBgCWpk4nmZwsnyuis6OTya2TbvoOLCt7FfLe9ra35cYbb8xNN92URqORTZs2ZdOmTdm2bVuS8r5BOy9Iv/XWW3PfffflzW9+85wL0gGACpqYSIaGktWry+eJiW5X9LxNTE5k6MqhrH7/6gxdOZSJyeW/T0Bv2Ktr8p7pJq4f+9jH8uY3vzlJOdv3u7/7u/nQhz40e0H6H/3RH+2yWtkzcQ0CALvTa+PDstrfTqcMdu12UhRJrZY0Gsnmzcmz3FR+Kevs6GToyqG0Z9opUqSWWhr9jWy+aLNbRwBdsyDX5O1JHqzVarnsssty2WWX7c1bAwDLVauVPPUWSEVRtlutZGR53m6h1W5leubJfSpSZHpmOq12yy0kgCVvr07XBADYRbOZDAyUM3hJ+TwwUG5fppqNZgb6B1JLuU+11DLQP5BmY/nuE9A7hDwA4Pmp15Px8fIUzaR8Hh9ftqdqJuVN38fXjafRX+5To7+R8XXjTtUEloV9vhk6AMCs0dHyGrxWq5zBW8YBb6fRkdFsvmhzWu1Wmo2mgAcsG0IeADA/6vVlew3eM6n31V2DByw7TtcEAACoECEPAACgQoQ8AACAChHyAAAAKkTIAwAAqBAhDwAAoEKEPABgjs6OTia3Tqazo9PtUgDYB0IeADBrYnIiQ1cOZfX7V2foyqFMTE50uyQA9pKQBwAkKWfwxm4ZS3umnSRpz7QzdsuYGT2AZUbIAwCSJK12K9Mz0ylSJEmKFJmemU6r3epyZQDsDSEPAEiSNBvNDPQPpJZakqSWWgb6B9JsNLtcGQB7Q8gDAJIk9b56xteNp9HfSJI0+hsZXzeeel+9y5UBsDf6ul0AALB0jI6MZvNFm9Nqt9JsNAU8gGXITB4AMEe9r56RQ0YWLOC5RQPAwhLyAIBF4xYNAAtPyAMAFoVbNAAsDiEPAFgUbtEAsDiEPABgUbhFA8DiEPIAgEXhFg0Ai8MtFACAReMWDQALT8gDABbVzls0ALAwnK4JAFRDp5NMTpbPAD1MyAMAlr+JiWRoKFm9unyecP89oHcJeQDA8tbpJGNjSbu8/17a7bJtRg/oUUIeALC8tVrJ9HRSlPffS1GU7Zb77wG9ScgDAJa3ZjMZGEhq5f33UquV7ab77wG9ScgDAJa3ej0ZH08a5f330miU7brbMwC9yS0UAIDlb3Q02by5PEWz2RTwgJ4m5AEA1VCvJyPuvwfgdE0AAIAKEfIAAAAqRMgDAACoECEPAACgQoQ8AACAChHyAIBF1ekkk5PlMwDzT8gDABbNxEQyNJSsXl0+T0x0uyKA6hHyAIBF0ekkY2NJu1222+2ybUYPYH4JeQDAomi1kunppCjKdlGU7Varu3UBVI2QBwAsimYzGRhIarWyXauV7Wazu3UBVI2QBwAsino9GR9PGo2y3WiU7Xq9u3UBVE1ftwsAAHrH6GiyeXN5imazKeABLAQhDwBYVPV6MjLS7SoAqsvpmgAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkPcsOp1kcrJ8BgAAWA6EvGcwMZEMDSWrV5fPExPdrggAAOC5CXm70ekkY2NJu1222+2ybUYPAABY6oS83Wi1kunppCjKdlGU7Varu3UBAAA8FyFvN5rNZGAgqdXKdq1WtpvN7tYFAADwXIS83ajXk/HxpNEo241G2a7Xu1sXAADAc+nrdgFL1ehosnlzeYpmsyngAQAAy4OQ9yzq9WRkpNtVAAAA7DmnawIAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCXpLOjk4mt06ms6PT7VIAAACel54PeROTExm6ciir3786Q1cOZWJyotslAQAA7LOeDnmdHZ2M3TKW9kw7SdKeaWfsljEzegAAwLLV0yGv1W5lemY6RYokSZEi0zPTabVbXa4MAABg3/R0yGs2mhnoH0gttSRJLbUM9A+k2Wh2uTIAAIB909Mhr95Xz/i68TT6G0mSRn8j4+vGU++rd7kyAFgYnU4yOVk+A1BNfd0uoNtGR0az+aLNabVbaTaaAh4AlTUxkYyNJdPTycBAMj6ejI52uyoA5ltPz+TtVO+rZ+SQEQEPgMrqdMqA1y7XGku7XbbN6AFUj5AHAD2g1Spn8IpyrbEURdluWWsMoHKEPADoAc1meYpmrVxrLLVa2W5aawygcoQ8AOgB9Xp5DV6jXGssjUbZrrtSAaByen7hFQDoFaOjyebN5SmazaaAB1BVQh4A9JB6PRkZ6XYVACwkp2sCAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgA8xQc/+MG8/OUvz8DAQAYGBnLyySfn7/7u72ZfL4oil112WZrNZg488MCcdtppuf/+++e8x8zMTC688MIcfvjhOfjgg3P22Wfn4YcfXuxdAaBHCXkA8BQveMEL8vu///v54he/mC9+8Yv5yZ/8yfzsz/7sbJC74oorctVVV+Xaa6/N3XffneHh4Zxxxhlpt9uz77F+/frceuutufnmm3PnnXfm0UcfzVlnnZXHH3+8W7sFQA+pFUVRdLuIp5qens7g4GCmpqYyMDDQ7XIAWCK6OT4ceuih+YM/+IO85S1vSbPZzPr163PJJZckKWfthoaGcvnll+f888/P1NRUjjjiiNxwww153etelyRptVpZuXJlPvnJT+bMM8/co880HgLwdHs6NpjJA4Bn8Pjjj+fmm2/O9773vZx88sl58MEHs2nTpqxdu3a2T39/f0499dTcddddSZJ77rknjz322Jw+zWYzxx577Gyf3ZmZmcn09PScBwDsCyHvKTo7OpncOpnOjk63SwGgi77yla/k3/27f5f+/v689a1vza233ppjjjkmmzZtSpIMDQ3N6T80NDT72qZNm3LAAQfkkEMOecY+u7Nx48YMDg7OPlauXDnPewVArxDyfmBiciJDVw5l9ftXZ+jKoUxMTnS7JAC65MUvfnG+/OUv53Of+1x+9Vd/Neeee26++tWvzr5eq9Xm9C+KYpdtT/dcfS699NJMTU3NPh566KHntxMA9CwhL+UM3tgtY2nPlBfNt2faGbtlzIweQI864IAD8kM/9EN55StfmY0bN+aHf/iH84d/+IcZHh5Okl1m5LZs2TI7uzc8PJzt27dn69atz9hnd/r7+2dX9Nz5AIB9IeQlabVbmZ6ZTpFyDZoiRaZnptNqt7pcGQBLQVEUmZmZyapVqzI8PJzbb7999rXt27fnjjvuyCmnnJIkOeGEE7L//vvP6fPII4/kvvvum+0DAAupr9sFLAXNRjMD/QNpz7RTpEgttTT6G2k2mt0uDYBF9pu/+Zt59atfnZUrV6bdbufmm2/OZz/72XzqU59KrVbL+vXrs2HDhqxZsyZr1qzJhg0bctBBB+Wcc85JkgwODua8887LO9/5zhx22GE59NBDc9FFF+W4447L6aef3uW9A6AXCHlJ6n31jK8bz9gtY5memU6jv5HxdeOp99W7XRoAi2zz5s1505velEceeSSDg4N5+ctfnk996lM544wzkiQXX3xxtm3blgsuuCBbt27NiSeemNtuuy2NRmP2Pa6++ur09fVl3bp12bZtW0ZHR3PddddlxYoV3dotAHqI++Q9RWdHJ612K81GU8ADWGJ67b5xvba/ADy3PR0bzOQ9Rb2vnpFDRrpdBgAAwD6z8AoAAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAV0lMhr9NJJifLZwAAgCrqmZA3MZEMDSWrV5fPExPdrggAAGD+9UTI63SSsbGk3S7b7XbZNqMHAABUTU+EvFYrmZ5OiqJsF0XZbrW6WxcAAMB864mQ12wmAwNJrVa2a7Wy3Wx2ty4AAID51hMhr15PxseTRqNsNxplu17vbl0AAADzra/bBSyW0dFk8+byFM1mU8ADAACqqWdCXlIGu5GRblcBAACwcHridE0A6FXuEQvQe4Q8AKgo94gF6E1CHgBUkHvEAvSuSoe8zo5OJrdOprPDiAZAb3GPWIDeVdmQNzE5kaErh7L6/aszdOVQJiadowJA73CPWIDeVcmQ19nRydgtY2nPlOeotGfaGbtlzIweAD3DPWIBelclb6HQarcyPTM92y5SZHpmOq12KyOHuIcCAL3BPWIBelMlZ/KajWYG+gdSS3mOSi21DPQPpNlwjgoAvWXnPWIFPIDeUcmQV++rZ3zdeBr95Tkqjf5GxteNp95nhAMAAKqtkqdrJsnoyGg2X7Q5rXYrzUZTwAMAAHpCZUNeUs7ouQYPAADoJZU8XRMAAKBXVTbkdTrJ5GT5DAC9qrOjk8mtk24jBNBDKhnyJiaSoaFk9eryecJ90AHoQROTExm6ciir3786Q1cOZWLSgAjQCyoX8jqdZGwsaZf3QU+7XbbN6AHQSzo7Ohm7ZSztmXJAbM+0M3bLmBk9gB5QuZDXaiXT00lRlO2iKNutVnfrAoDF1Gq3Mj0znSLlgFikyPTMdFptAyJA1VUu5DWbycBAUivvg55arWw33QcdgB7SbDQz0D+QWsoBsZZaBvoH0mwYEAGqrnIhr15PxseTRnkf9DQaZbvuNnkA9JB6Xz3j68bT6C8HxEZ/I+Prxt03FqAHVPI+eaOjyebN5SmazaaAB0BvGh0ZzeaLNqfVbqXZaAp4AD2ikiEvKYPdiPugA9Dj6n31jBxiQAToJZU7XRMAAKCXCXkAAAAVIuQBAABUyIKFvD/+4z/OqlWrUq/Xc8IJJ+Qf/uEfFuqjAAAA+IEFCXl/8Rd/kfXr1+c973lP7r333vw//8//k1e/+tX59re/vRAfBwAAwA8sSMi76qqrct555+WXfumX8tKXvjTXXHNNVq5cmQ9+8IML8XEAAAD8wLyHvO3bt+eee+7J2rVr52xfu3Zt7rrrrl36z8zMZHp6es4DAACAfTPvIe+73/1uHn/88QwNDc3ZPjQ0lE2bNu3Sf+PGjRkcHJx9rFy5cr5LAgAA6BkLtvBKrVab0y6KYpdtSXLppZdmampq9vHQQw8tVEkAAACV1zffb3j44YdnxYoVu8zabdmyZZfZvSTp7+9Pf3//fJcBAADQk+Z9Ju+AAw7ICSeckNtvv33O9ttvvz2nnHLKfH8cAAAATzHvM3lJ8hu/8Rt505velFe+8pU5+eST8+EPfzjf/va389a3vnUhPg4AAIAfWJCQ97rXvS7/+q//mve+97155JFHcuyxx+aTn/xkjj766IX4OAAAAH5gQUJeklxwwQW54IILFurtAQAA2I0FW10TAACAxSfkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIUIeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAVIuQBAABUiJAHAABQIX3dLuDpiqJIkkxPT3e5EgCWkp3jws5xouqMhwA83Z6OhUsu5LXb7STJypUru1wJAEtRu93O4OBgt8tYcMZDAJ7Jc42FtWKJ/ST6xBNPpNVqpdFopFar7fP7TE9PZ+XKlXnooYcyMDAwjxUuT47HkxyLuRyPuRyPJy21Y1EURdrtdprNZvbbr/pXGxgP559jMZfj8STHYi7HY66ldDz2dCxccjN5++23X17wghfM2/sNDAx0/V/GUuJ4PMmxmMvxmMvxeNJSOha9MIO3k/Fw4TgWczkeT3Is5nI85loqx2NPxsLq/xQKAADQQ4Q8AACACqlsyOvv789//a//Nf39/d0uZUlwPJ7kWMzleMzleDzJsagG/x6f5FjM5Xg8ybGYy/GYazkejyW38AoAAAD7rrIzeQAAAL1IyAMAAKgQIQ8AAKBChDwAAIAKEfIAAAAqZFmFvD/+4z/OqlWrUq/Xc8IJJ+Qf/uEfnrX/HXfckRNOOCH1ej0jIyP5kz/5k136fPzjH88xxxyT/v7+HHPMMbn11lsXqvx5Nd/H4rrrrkutVtvl0el0FnI35s3eHI9HHnkk55xzTl784hdnv/32y/r163fbrxe+G3tyLHrpuzE+Pp4zzjgjRxxxRAYGBnLyySfn05/+9C79lut3I5n/47Hcvx/LkbFwLuPhk4yFcxkPn2QsnKsnxsJimbj55puL/fffv/jIRz5SfPWrXy3e8Y53FAcffHDxrW99a7f9Jycni4MOOqh4xzveUXz1q18tPvKRjxT7779/8Vd/9Vezfe66665ixYoVxYYNG4qvfe1rxYYNG4q+vr7ic5/73GLt1j5ZiGPxsY99rBgYGCgeeeSROY/lYG+Px4MPPli8/e1vL66//vriFa94RfGOd7xjlz698t3Yk2PRS9+Nd7zjHcXll19efOELXyi+8Y1vFJdeemmx//77F1/60pdm+yzX70ZRLMzxWM7fj+XIWDiX8fBJxsK5jIdPMhbO1Stj4bIJeT/2Yz9WvPWtb52z7SUveUnx7ne/e7f9L7744uIlL3nJnG3nn39+cdJJJ822161bV/zUT/3UnD5nnnlm8frXv36eql4YC3EsPvaxjxWDg4PzXuti2Nvj8VSnnnrqbv9D3ivfjad6pmPRq9+NnY455pjid3/3d2fby/W7URQLczyW8/djOTIWzmU8fJKxcC7j4ZOMhXP1yli4LE7X3L59e+65556sXbt2zva1a9fmrrvu2u3f/NM//dMu/c8888x88YtfzGOPPfasfZ7pPZeChToWSfLoo4/m6KOPzgte8IKcddZZuffee+d/B+bZvhyPPdEr34091avfjSeeeCLtdjuHHnro7Lbl+N1IFu54JMvz+7EcGQvnMh4+yVg4l/HwScbCuXppLFwWIe+73/1uHn/88QwNDc3ZPjQ0lE2bNu32bzZt2rTb/jt27Mh3v/vdZ+3zTO+5FCzUsXjJS16S6667Lp/4xCfy53/+56nX63nVq16Vb37zmwuzI/NkX47HnuiV78ae6OXvxvve975873vfy7p162a3LcfvRrJwx2O5fj+WI2PhXMbDJxkL5zIePslYOFcvjYV9XfvkfVCr1ea0i6LYZdtz9X/69r19z6Vivo/FSSedlJNOOmn29Ve96lX5kR/5kXzgAx/I+9///vkqe8EsxL/HXvluPJde/W78+Z//eS677LL8zd/8TY488sh5ec+lYL6Px3L/fixHxsK5jIdPMhbOZTx8krFwrl4YC5dFyDv88MOzYsWKXRL2li1bdkniOw0PD++2f19fXw477LBn7fNM77kULNSxeLr99tsvP/qjP7qkf51K9u147Ile+W7si174bvzFX/xFzjvvvPzlX/5lTj/99DmvLcfvRrJwx+Pplsv3YzkyFs5lPHySsXAu4+GTjIVz9dJYuCxO1zzggANywgkn5Pbbb5+z/fbbb88pp5yy2785+eSTd+l/22235ZWvfGX233//Z+3zTO+5FCzUsXi6oijy5S9/OUcdddT8FL5A9uV47Ile+W7si6p/N/78z/88b37zm3PTTTflNa95zS6vL8fvRrJwx+Pplsv3YzkyFs5lPHySsXAu4+GTjIVz9dRYuDjruzx/O5c7/ehHP1p89atfLdavX18cfPDBxb/8y78URVEU7373u4s3velNs/13LpP867/+68VXv/rV4qMf/eguyyT/4z/+Y7FixYri93//94uvfe1rxe///u8vi+VfF+JYXHbZZcWnPvWp4oEHHijuvffe4hd/8ReLvr6+4vOf//yi79/e2tvjURRFce+99xb33ntvccIJJxTnnHNOce+99xb333//7Ou98t0oiuc+Fr303bjpppuKvr6+4o/+6I/mLIH8b//2b7N9lut3oygW5ngs5+/HcmQsnMt4+CRj4VzGwycZC+fqlbFw2YS8oiiKP/qjPyqOPvro4oADDih+5Ed+pLjjjjtmXzv33HOLU089dU7/z372s8Xxxx9fHHDAAcWLXvSi4oMf/OAu7/mXf/mXxYtf/OJi//33L17ykpcUH//4xxd6N+bFfB+L9evXFy984QuLAw44oDjiiCOKtWvXFnfddddi7Mq82NvjkWSXx9FHHz2nT698N57rWPTSd+PUU0/d7fE499xz57zncv1uFMX8H4/l/v1YjoyFcxkPn2QsnMt4+CRj4Vy9MBbWiuIHVxwDAACw7C2La/IAAADYM0IeAABAhQh5AAAAFSLkAQAAVIiQBwAAUCFCHgAAQIUIeQAAABUi5AEAAFSIkAcAAFAhQh4AAECFCHkAAAAV8v8HDejBO2KH2e8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x1000 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = ['red', 'blue', 'green']\n",
    "camera_fps=analysis_methods.get(\"camera_fps\")\n",
    "fig, axes = plt.subplots(\n",
    "        nrows=1, ncols=2, figsize=(9,10), tight_layout=True\n",
    "    )\n",
    "ax1, ax2 = axes.flatten()\n",
    "for keys, grp in all_evaluation.groupby(['animal_id']):\n",
    "    if grp['batchs'].values[0]=='1st':\n",
    "        this_color=colors[0]\n",
    "    elif grp['batchs'].values[0]=='2nd':\n",
    "        this_color=colors[1]\n",
    "    else:\n",
    "        this_color=colors[2]\n",
    "    travel_distance=grp['travel_distance'].sum()/camera_fps\n",
    "    turning_distance=grp['turning_distance'].sum()/camera_fps\n",
    "    num_follow_epochs=grp['num_follow_epochs'].sum()/grp['number_frames'].sum()\n",
    "    ax1.scatter(grp['num_follow_epochs'].sum()/grp['number_frames'].sum(),grp['travel_distance'].sum()/camera_fps,c=this_color,s=8)\n",
    "    ax2.scatter(grp['num_follow_epochs'].sum()/grp['number_frames'].sum(),grp['turning_distance'].sum()/camera_fps,c=this_color,s=8)\n",
    "for keys, grp in all_evaluation.groupby(['batchs']):\n",
    "    if grp['batchs'].values[0]=='1st':\n",
    "        this_color=colors[0]\n",
    "    elif grp['batchs'].values[0]=='2nd':\n",
    "        this_color=colors[1]\n",
    "    else:\n",
    "        this_color=colors[2]\n",
    "    ax1.scatter(grp['num_follow_epochs'].sum()/grp['number_frames'].sum(),grp['travel_distance'].sum()/camera_fps/grp['animal_id'].unique().shape[0],edgecolors=this_color,s=20,facecolors='none')\n",
    "    ax2.scatter(grp['num_follow_epochs'].sum()/grp['number_frames'].sum(),grp['turning_distance'].sum()/camera_fps/grp['animal_id'].unique().shape[0],edgecolors=this_color,s=20,facecolors='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "#chop the data into 3 groups depends on time. Then we can analyse the 1st, 2nd and 3rd of the data\n",
    "all_evaluation['time_group'] = pd.cut(all_evaluation['trial_id'], bins=3, labels=['1st', '2nd', '3rd'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_follow_list=[]\n",
    "travel_distance_list=[]\n",
    "for keys, grp in all_evaluation.groupby(['animal_id','time_group']):\n",
    "    num_follow_list.append(grp[\"num_follow_epochs\"].sum())\n",
    "    travel_distance_list.append(grp[\"travel_distance\"].sum())\n",
    "num_follow_arr=np.array(num_follow_list).reshape(-1,3)\n",
    "travel_distance_arr=np.array(travel_distance_list).reshape(-1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(\n",
    "        nrows=1, ncols=2, figsize=(9,10), tight_layout=True\n",
    "    )\n",
    "ax1, ax2 = axes.flatten()\n",
    "ax1.plot(np.transpose(num_follow_arr), 'o-')\n",
    "ax1.plot(np.mean(num_follow_arr,axis=0),'k',linewidth=5)\n",
    "ax1.plot(np.median(num_follow_arr,axis=0),'k',linewidth=5,alpha=0.6)\n",
    "ax2.plot(np.transpose(travel_distance_arr), 'o-')\n",
    "ax2.plot(np.mean(travel_distance_arr,axis=0),'k',linewidth=5)\n",
    "ax2.plot(np.median(travel_distance_arr,axis=0),'k',linewidth=5,alpha=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarise_as_csv(input_list,agent_type='glocust'):\n",
    "    df=pd.DataFrame(np.hstack(input_list))\n",
    "    df.insert(0, 'agent_type', np.repeat(agent_type,df.shape[0]).tolist())\n",
    "    df.insert(0, 'ID', np.repeat(np.arange(len(input_list)),int(df.shape[0]/len(input_list))).tolist())\n",
    "    df.to_csv(f'{agent_type}.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summarise_as_csv(exp_agent_list,'exp_agent')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summarise_as_csv(exp_agent_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Section 2.1: Using 1D histogram to plot proportion of time and 2D histogram to plot relative position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_trials=pd.concat(relative_pos_all_animals)\n",
    "all_simulated_trials=pd.concat(simulated_relative_pos_all_animals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extract_follow_epoches=analysis_methods.get(\"extract_follow_epoches\",True)\n",
    "distribution_with_entire_body=analysis_methods.get(\"distribution_with_entire_body\",True)\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(all_trials[\"type\"].unique())*len(all_trials[\"degree\"].unique()),figsize=(18, 6),tight_layout=True,sharex=True, sharey=True)\n",
    "i=0\n",
    "if extract_follow_epoches:\n",
    "    xlimit=(-5,40)\n",
    "    ylimit=(-15,15)\n",
    "else:\n",
    "    xlimit=(-20,100)\n",
    "    ylimit=(-45,45)\n",
    "for keys, grp in all_trials.groupby(['type','degree']):\n",
    "        print(keys)\n",
    "        if distribution_with_entire_body:\n",
    "            body_points=generate_points_within_rectangles(grp['x'].values,grp['y'].values,1,4,2,21)\n",
    "            axes[i].hist2d(body_points[:,0],body_points[:,1],bins=1000)\n",
    "        else:\n",
    "            axes[i].hist2d(grp['x'].values,grp['y'].values,bins=400)\n",
    "        axes[i].set(\n",
    "        yticks=[ylimit[0],0,ylimit[1]],\n",
    "        xticks=[-5,0,4,20,40],\n",
    "        xlim=xlimit,ylim=ylimit,title=f'agent:{keys[0]},deg:{int(keys[1])}',adjustable='box', aspect='equal')\n",
    "        #rect = patches.Rectangle((0,-1), 6, 2, linewidth=1, edgecolor='red',linestyle=\"--\",facecolor='none')\n",
    "        #cir1 = patches.Circle((4,0), radius=0.5, linewidth=0.5, edgecolor='yellow')\n",
    "        cir2 = patches.Circle((0,0), radius=0.4, linewidth=0.4, edgecolor='red',facecolor='red')\n",
    "        #axes[i].add_patch(cir1)\n",
    "        axes[i].add_patch(cir2)\n",
    "        i=i+1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=len(all_simulated_trials[\"type\"].unique())*len(all_simulated_trials[\"degree\"].unique()),figsize=(18, 6),tight_layout=True,sharex=True, sharey=True)\n",
    "i=0\n",
    "if extract_follow_epoches:\n",
    "    xlimit=(-5,40)\n",
    "    ylimit=(-15,15)\n",
    "else:\n",
    "    xlimit=(-20,100)\n",
    "    ylimit=(-45,45)\n",
    "for keys, grp in all_simulated_trials.groupby(['type','degree']):\n",
    "        if distribution_with_entire_body:\n",
    "            body_points=generate_points_within_rectangles(grp['x'].values,grp['y'].values,1,4,2,21)\n",
    "            axes[i].hist2d(body_points[:,0],body_points[:,1],bins=1000)\n",
    "        else:\n",
    "            axes[i].hist2d(grp['x'].values,grp['y'].values,bins=400)\n",
    "        axes[i].set(\n",
    "        yticks=[ylimit[0],0,ylimit[1]],\n",
    "        xticks=[-5,0,4,20,40],\n",
    "        xlim=xlimit,ylim=ylimit,title=f'agent:{keys[0]},deg:{int(keys[1])}',adjustable='box', aspect='equal')\n",
    "        #rect = patches.Rectangle((0,-1), 6, 2, linewidth=1, edgecolor='red',linestyle=\"--\",facecolor='none')\n",
    "        #cir1 = patches.Circle((4,0), radius=0.5, linewidth=0.5, edgecolor='yellow')\n",
    "        cir2 = patches.Circle((0,0), radius=0.4, linewidth=0.4, edgecolor='red',facecolor='red')\n",
    "        #axes[i].add_patch(cir1)\n",
    "        axes[i].add_patch(cir2)\n",
    "        i=i+1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##plot the distribution of follow epochs in a trial\n",
    "for keys, grp in all_trials.groupby(['type','degree']):\n",
    "    sim_grp=all_simulated_trials[(all_simulated_trials['type']==keys[0])&(all_simulated_trials['degree']==keys[1])]\n",
    "    print(keys)\n",
    "    fig, axes = plt.subplots(\n",
    "        nrows=1, ncols=2, figsize=(9,5), tight_layout=True\n",
    "    )\n",
    "    ax, ax2 = axes.flatten()\n",
    "    ax.hist(grp['ts'].values,bins=100,density=False,color='r')\n",
    "    ax.hist(sim_grp['ts'].values,bins=100,density=False,color='tab:gray',alpha=0.3)\n",
    "    ax2.hist(grp['ts'].values,bins=100,color='r',density=True,histtype=\"step\",cumulative=True,label=\"Cumulative histogram\")\n",
    "    ax2.hist(sim_grp['ts'].values,bins=100,color='tab:gray',alpha=0.3,density=True,histtype=\"step\",cumulative=True,label=\"Cumulative histogram\")\n",
    "    ax.set(xlim=(0,60),title=f'agent:{keys[0]},deg:{int(keys[1])}')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "#https://stackoverflow.com/questions/78384537/fitting-2d-histograms-with-2d-gaussians\n",
    "def twoD_Gaussian(xy, amplitude, xo, yo, sigma_x, sigma_y, rho):\n",
    "    x, y = xy\n",
    "    xo = float(xo)\n",
    "    yo = float(yo)\n",
    "    a = 1 / sigma_x ** 2\n",
    "    b = rho / sigma_x / sigma_y\n",
    "    c = 1 / sigma_y ** 2\n",
    "    exponent=  - (1 / (2 * (1 - rho ** 2)))\\\n",
    "                 * (a * (x - xo) ** 2 - 2 * b * (x - xo) * (y - yo) + c * (y - yo) ** 2)\n",
    "    g = amplitude * np.exp(exponent)\n",
    "    return g.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbr_bins = 100\n",
    "B, xedges, yedges = np.histogram2d(grp['x'].values,grp['y'].values, bins=nbr_bins)\n",
    "extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]\n",
    "bin_centers_bx = (xedges[:-1] + xedges[1:]) / 2.0\n",
    "bin_centers_by = (yedges[:-1] + yedges[1:]) / 2.0\n",
    "X, Y = np.meshgrid(bin_centers_bx, bin_centers_by)\n",
    "data = B.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = [0, 4]\n",
    "p0 = (B.max(),mean[0],mean[1],1.7,1.0,0.9)\n",
    "coeff, var_matrix = curve_fit(twoD_Gaussian, (X, Y), data, p0=p0)\n",
    "print('hist fit', coeff)\n",
    "print('var_matrix', var_matrix)\n",
    "data_fitted_hist = twoD_Gaussian((X, Y), *coeff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "xlimit=(-30,30)\n",
    "ylimit=(-30,30)\n",
    "ax[0].imshow(data.reshape((nbr_bins, nbr_bins)), origin=\"lower\",\n",
    "             extent=extent, interpolation=\"none\")\n",
    "ax[1].imshow(data.reshape((nbr_bins, nbr_bins)), origin=\"lower\",\n",
    "             extent=extent, interpolation=\"none\")\n",
    "ax[1].contour(X, Y, data_fitted_hist.reshape(nbr_bins, nbr_bins), 2, colors=\"w\",linestyles='dashed',linewidths=0.5)\n",
    "ax[0].set(\n",
    "        yticks=[-30,0.0,30],\n",
    "        xticks=[-30,0.0,30],\n",
    "        xlim=xlimit,ylim=ylimit)\n",
    "ax[1].set(\n",
    "        yticks=[-30,0.0,30],\n",
    "        xticks=[-30,0.0,30],\n",
    "        xlim=xlimit,ylim=ylimit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_methods.update({\"distribution_with_entire_body\": True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##plot the spatial distribution of virtual locusts in a 2D histogram\n",
    "\n",
    "distribution_with_entire_body=analysis_methods.get(\"distribution_with_entire_body\",True)\n",
    "xlimit=(-30,30)\n",
    "ylimit=(-30,30)\n",
    "for keys, grp in all_trials.groupby(['type']):\n",
    "    fig, ax = plt.subplots(dpi=300, figsize=(2,2))\n",
    "    if distribution_with_entire_body:\n",
    "        body_points=generate_points_within_rectangles(grp['x'].values,grp['y'].values,1,4,2,21)\n",
    "        ax.hist2d(body_points[:,0],body_points[:,1],bins=1000)\n",
    "    else:\n",
    "        ax.hist2d(grp['x'].values,grp['y'].values,bins=100)\n",
    "    ax.set(adjustable='box', aspect='equal')\n",
    "    ax.set(\n",
    "        yticks=[-30,0.0,30],\n",
    "        xticks=[-30,0.0,30],\n",
    "        xticklabels=(['-0.3', '0.0', '0.3']),\n",
    "        yticklabels=(['-0.3', '0.0', '0.3']),\n",
    "        xlim=xlimit,ylim=ylimit,title=f'agent:{keys[0]}')\n",
    "    '''\n",
    "    the position of virtual locust if their butt is at the origin\n",
    "    rect = patches.Rectangle((0,-1), 6, 2, linewidth=1, edgecolor='red',linestyle=\"--\",facecolor='none')\n",
    "    ax.add_patch(rect)\n",
    "    '''\n",
    "    cir = patches.Circle((0,0), radius=0.5, linewidth=0.5, edgecolor='red')\n",
    "    ax.add_patch(cir)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##plot the temporal distribution of follow epochs in a trial\n",
    "for keys, grp in all_trials.groupby(['type']):\n",
    "    print(keys)\n",
    "    sim_grp=all_simulated_trials[all_simulated_trials['type']==keys[0]]\n",
    "    fig, axes = plt.subplots(\n",
    "        nrows=1, ncols=2, figsize=(9,5), tight_layout=True\n",
    "    )\n",
    "    ax, ax2 = axes.flatten()\n",
    "    ax.hist(grp['ts'].values,bins=100,density=False,color='r')\n",
    "    ax.hist(sim_grp['ts'].values,bins=100,density=False,color='tab:gray',alpha=0.5)\n",
    "    ax2.hist(grp['ts'].values,bins=100,density=True,color='r',histtype=\"step\",cumulative=True,label=\"Cumulative histogram\")\n",
    "    ax2.hist(sim_grp['ts'].values,bins=100,color='tab:gray',alpha=0.5,density=True,histtype=\"step\",cumulative=True,label=\"Cumulative histogram\")\n",
    "    ax.set(xlim=(0,60))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Section 2.2: save data into Mat file if needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##save the raw data into mat file for matlab\n",
    "from scipy.io import savemat\n",
    "all_raw=pd.concat(relative_pos_all_animals)\n",
    "#all_raw['trial_id']=all_raw['trial_id'].astype(int)\n",
    "data_dict = {name: col.values for name, col in all_raw.items()}\n",
    "summary_file_name = Path(thisDataset) /\"time_series_curated.mat\"\n",
    "savemat(summary_file_name, data_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 3.0, analysis walking behavours before and after the presence of stimulus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fix_data_type(all_trials):\n",
    "    all_trials['id'] = all_trials['id'].astype(int)\n",
    "    all_trials['mu'] = all_trials['mu'].astype(int)\n",
    "    all_trials['velocity'] = all_trials['velocity'].astype(float)\n",
    "    all_trials['omega'] = all_trials['omega'].astype(float)\n",
    "    all_trials['normalised_v'] = all_trials['normalised_v'].astype(float)\n",
    "    all_trials['normalised_omega'] = all_trials['normalised_omega'].astype(float)\n",
    "    if 'density' in all_trials.columns:\n",
    "        all_trials['density'] = all_trials['density'].astype(int)\n",
    "    return all_trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "if type(raster_across_animals)==list:\n",
    "    all_trials=pd.concat(raster_across_animals)\n",
    "    all_trials=fix_data_type(all_trials)\n",
    "else:\n",
    "    all_trials=fix_data_type(raster_across_animals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_baseline_distribution(all_trials,analysis_methods,metrics_name='velocity',duration_for_baseline=3):\n",
    "    monitor_fps=analysis_methods.get(\"monitor_fps\")\n",
    "    align_with_isi_onset=analysis_methods.get(\"align_with_isi_onset\",False)\n",
    "    these_baselines=[]\n",
    "    for keys, this_data in all_trials.groupby(['animal_id','id']):\n",
    "        #print(this_data['object'][1])\n",
    "        if align_with_isi_onset:\n",
    "            if this_data['object'][1]!='empty_trial':\n",
    "            #if keys[1]%2==0:#for the Swarm scene to only use stim trials to get the baseline\n",
    "                this_metrics=this_data[metrics_name].values\n",
    "                #print(this_metrics)\n",
    "                these_baselines.append(np.mean(this_metrics[-duration_for_baseline*monitor_fps:]))\n",
    "        \n",
    "        else:\n",
    "            if this_data['object'][1]=='empty_trial':\n",
    "            #if keys[1]%2!=0:#for the Swarm scene to only use ISI trials\n",
    "                this_metrics=this_data[metrics_name].values\n",
    "                these_baselines.append(np.mean(this_metrics[-duration_for_baseline*monitor_fps:]))\n",
    "    fig, axes = plt.subplots(\n",
    "        nrows=1, ncols=1, figsize=(9,9), tight_layout=True\n",
    "    )\n",
    "    axes.hist(np.vstack(these_baselines),bins=500)\n",
    "    if metrics_name=='velocity':\n",
    "        axes.set(xlim=(0,10),ylim=(0, 250))\n",
    "    else:\n",
    "        #axes.set(xlim=(-0.001,0.001),ylim=(0, 60))#if used rad\n",
    "        axes.set(xlim=(-3,3),ylim=(0, 400))\n",
    "\n",
    "    plt.minorticks_on()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics_name='omega'\n",
    "check_baseline_distribution(all_trials,analysis_methods,metrics_name,3)\n",
    "# check_baseline_distribution(all_trials,analysis_methods)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_trials(analysis_methods,all_trials,metrics_name='velocity',metrics_name2='normalised_v',walk_threshold=1,duration_for_baseline=3):\n",
    "    monitor_fps=analysis_methods.get(\"monitor_fps\")\n",
    "    movement_trial_boolean=[]\n",
    "    these_metrics=[]\n",
    "    these_normalised_metrics=[]\n",
    "    for keys, this_data in all_trials.groupby(['animal_id','id']):\n",
    "        this_metrics=this_data[metrics_name].values\n",
    "        baseline_metrics=np.mean(this_metrics[-duration_for_baseline*monitor_fps:])\n",
    "        movement_trial_boolean.append(abs(baseline_metrics)>walk_threshold)\n",
    "        these_metrics.append(this_metrics)\n",
    "        these_normalised_metrics.append(this_data[metrics_name2].values)\n",
    "    return movement_trial_boolean,these_metrics,these_normalised_metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for simplicity, classifying trial only based on stationary and walk trails\n",
    "metrics_name='omega'\n",
    "_,these_metrics,these_normalised_metrics=split_trials(analysis_methods,all_trials,metrics_name,'normalised_omega',1)#1 degree or 0.0002 rad\n",
    "movement_trial_boolean,_,_=split_trials(analysis_methods,all_trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "#movement_trial_boolean,these_metrics,these_normalised_metrics=split_trials(analysis_methods,all_trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_trial_index(movement_trial_boolean,analysis_methods):\n",
    "    align_with_isi_onset=analysis_methods.get(\"align_with_isi_onset\",False)\n",
    "#int((trial_id)/2) means the number of stimulus trial\n",
    "    after_movement_ith_trial=[]\n",
    "    after_no_movement_ith_trial=[]\n",
    "    if align_with_isi_onset:\n",
    "        after_movement_ith_trial=[i+1 for i, x in enumerate(movement_trial_boolean[1::2]) if x and i % int((trial_id)/2) != (trial_id)/2-1]\n",
    "        after_no_movement_ith_trial=[i+1 for i, x in enumerate(movement_trial_boolean[1::2]) if x==False and i % int((trial_id)/2) != (trial_id)/2-1]\n",
    "    else:\n",
    "        after_movement_ith_trial=[i for i, x in enumerate(movement_trial_boolean[::2]) if x]\n",
    "        after_no_movement_ith_trial=[i for i, x in enumerate(movement_trial_boolean[::2]) if x==False]\n",
    "    return after_movement_ith_trial,after_no_movement_ith_trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "after_movement_ith_trial,after_no_movement_ith_trial=extract_trial_index(movement_trial_boolean,analysis_methods)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fill_between_range(data,mean_data,using_confidence_interval=True):\n",
    "\n",
    "    ##to plot distribution with 95% confidence interval with t distribution (since the sample is usually not big)\n",
    "    if using_confidence_interval:\n",
    "        confidence_level = 0.95\n",
    "        #cl95=st.t.interval(confidence=0.95, df=len(data)-1, loc=np.mean(data), scale=st.sem(data)) \n",
    "        cl95=st.norm.interval(confidence_level,loc=mean_data,scale=st.sem(data))\n",
    "        dif_y1=cl95[0][:]\n",
    "        dif_y2=cl95[1][:]\n",
    "    else:\n",
    "        sem_response = np.std(data, axis=0, ddof=1) / np.sqrt(data.shape[0])\n",
    "        dif_y1=mean_data + sem_response\n",
    "        dif_y2=mean_data - sem_response\n",
    "    return dif_y1,dif_y2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_visual_evoked_behaviour(these_metrics,these_normalised_metrics,after_movement_ith_trial,after_no_movement_ith_trial,analysis_methods,metrics_name='velocity',row_of_interest=None,type_key=\"\",variable_values=None):\n",
    "    exp_name=analysis_methods.get('experiment_name')\n",
    "    number_frame_scene_changing=analysis_methods.get(\"largest_unfilled_gap\",12)\n",
    "    number_frame_scene_changing=10#set an arbitrary value to escape from the effect of missing value when plotting the histogram\n",
    "    align_with_isi_onset=analysis_methods.get(\"align_with_isi_onset\",False)\n",
    "    save_output=analysis_methods.get(\"save_output\",False)\n",
    "    analysis_window=analysis_methods.get(\"analysis_window\")\n",
    "    monitor_fps=analysis_methods.get(\"monitor_fps\")\n",
    "    camera_fps=analysis_methods.get(\"camera_fps\")\n",
    "    all_animals=False\n",
    "    tmp=np.vstack(these_metrics)\n",
    "    tmp3=np.vstack(these_normalised_metrics)\n",
    "    if align_with_isi_onset:\n",
    "       stim_evoked_metrics=tmp[::2]\n",
    "       stim_evoked_norm_metrics=tmp3[::2]\n",
    "    else:\n",
    "       stim_evoked_metrics=tmp[1::2]\n",
    "       stim_evoked_norm_metrics=tmp3[1::2]\n",
    "\n",
    "    fig, axes = plt.subplots(\n",
    "            nrows=4, ncols=2, figsize=(9,10), tight_layout=True\n",
    "        )\n",
    "    ax1, ax2, ax3, ax4,ax5,ax6,ax7,ax8 = axes.flatten()\n",
    "    plt.rcParams.update(plt.rcParamsDefault)\n",
    "    plt.rcParams.update({'font.size': 8})\n",
    "    # Set the axis line width to 2\n",
    "    plt.rcParams['ytick.major.width'] = 1\n",
    "    plt.rcParams['xtick.major.width'] = 1\n",
    "    plt.rcParams['axes.linewidth'] = 1\n",
    "    cmap = plt.get_cmap('viridis')\n",
    "    if all_animals==False and type(row_of_interest)==pd.Series:\n",
    "        animal_interest_stationary=[]\n",
    "        animal_interest_move=[]\n",
    "        for i in np.where(row_of_interest)[0].tolist():\n",
    "            if i in after_no_movement_ith_trial:\n",
    "                animal_interest_stationary.append(i)\n",
    "            else:\n",
    "                animal_interest_move.append(i)\n",
    "        p3=stim_evoked_metrics[animal_interest_stationary,:]\n",
    "        p4=stim_evoked_metrics[animal_interest_move,:]\n",
    "        p5=stim_evoked_norm_metrics[animal_interest_stationary,:]\n",
    "        p6=stim_evoked_norm_metrics[animal_interest_move,:]\n",
    "\n",
    "    else:\n",
    "        p3=stim_evoked_metrics[after_no_movement_ith_trial,:]\n",
    "        p4=stim_evoked_metrics[after_movement_ith_trial,:]\n",
    "        p5=stim_evoked_norm_metrics[after_no_movement_ith_trial,:]\n",
    "        p6=stim_evoked_norm_metrics[after_movement_ith_trial,:]\n",
    "    x=np.arange(0,p3.shape[1])\n",
    "    p1=np.nancumsum(p3, axis=1)/camera_fps\n",
    "    p2=np.nancumsum(p4, axis=1)/camera_fps\n",
    "    ax1.plot(np.transpose(p1),linewidth=0.1)\n",
    "    mean_p1=np.nanmean(p1,axis=0)\n",
    "    ax1.plot(mean_p1,'k',linewidth=1)\n",
    "    #ax1.plot(np.nanmedian(p1,axis=0),'k--',linewidth=0.5)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p1,mean_p1)\n",
    "    ax1.fill_between(x,dif_y1,dif_y2, alpha=0.4,color='k')\n",
    "    ax2.plot(np.transpose(p2),linewidth=0.1)\n",
    "    mean_p2=np.nanmean(p2,axis=0)\n",
    "\n",
    "    ax2.plot(mean_p2,'k',linewidth=1)\n",
    "    #ax2.plot(np.nanmedian(p2,axis=0),'k--',linewidth=0.5)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p2,mean_p2)\n",
    "    ax2.fill_between(x,dif_y1,dif_y2, alpha=0.4,color='k')\n",
    "    ax3.plot(np.transpose(p3),linewidth=0.1)\n",
    "    mean_p3=np.nanmean(p3,axis=0)\n",
    "    ax3.plot(mean_p3,'k',linewidth=1)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p3,mean_p3)\n",
    "    ax3.fill_between(x,dif_y1,dif_y2, alpha=0.4,color='k')\n",
    "    #ax3.plot(np.nanmedian(p3,axis=0),'k--',linewidth=0.5)\n",
    "    mean_p4=np.nanmean(p4,axis=0)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p4,mean_p4)\n",
    "    ax4.plot(np.transpose(p4),linewidth=0.1)\n",
    "    ax4.plot(mean_p4,'k',linewidth=1)\n",
    "    ax4.fill_between(x,dif_y1,dif_y2, alpha=0.4,color='k')\n",
    "    #ax4.plot(np.nanmedian(p4,axis=0),'k--',linewidth=0.5)\n",
    "    mean_p5=np.nanmean(p5,axis=0)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p5,mean_p5)\n",
    "    ax5.plot(np.transpose(p5),linewidth=0.1)\n",
    "    ax5.plot(mean_p5,'k',linewidth=1)\n",
    "    ax5.fill_between(x,dif_y1,dif_y2, alpha=0.4,color='k')\n",
    "    #ax5.plot(np.nanmedian(p5,axis=0),'k--',linewidth=0.5)\n",
    "    mean_p6=np.nanmean(p6,axis=0)\n",
    "    dif_y1,dif_y2=get_fill_between_range(p6,mean_p6)\n",
    "    ax6.plot(np.transpose(p6),linewidth=0.1)\n",
    "    ax6.plot(mean_p6,'k',linewidth=1)\n",
    "    ax6.fill_between(x,dif_y1,dif_y2,alpha=0.4,color='k')\n",
    "\n",
    "    #ax6.plot(np.nanmedian(p6,axis=0),'k--',linewidth=0.5)\n",
    "    if metrics_name=='velocity' and all_animals==False:\n",
    "        ylimit=10\n",
    "        ylimit_log=100\n",
    "        ax1.set_ylim([0,5*ylimit])\n",
    "        ax2.set_ylim([0,5*ylimit])\n",
    "        ax3.set_ylim([0,1*ylimit])\n",
    "        ax4.set_ylim([0,1*ylimit])\n",
    "    elif all_animals==False:\n",
    "        ylimit=80\n",
    "        ylimit_log=1000\n",
    "        # ax1.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax1.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax2.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax2.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax3.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax3.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax4.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax4.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        ax1.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax2.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax3.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax4.set_ylim([-1*ylimit,1*ylimit])\n",
    "    elif metrics_name=='velocity':\n",
    "        ylimit=10\n",
    "        ylimit_log=100\n",
    "        ax1.set_ylim([0,5*ylimit])\n",
    "        ax2.set_ylim([0,5*ylimit])\n",
    "        ax3.set_ylim([0,1*ylimit])\n",
    "        ax4.set_ylim([0,1*ylimit])\n",
    "    else:\n",
    "        ylimit=80\n",
    "        ylimit_log=1000\n",
    "        # ax1.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax1.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax2.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax2.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax3.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax3.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        # ax4.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))\n",
    "        # ax4.yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))\n",
    "        ax1.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax2.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax3.set_ylim([-1*ylimit,1*ylimit])\n",
    "        ax4.set_ylim([-1*ylimit,1*ylimit])\n",
    "    ax1.set(\n",
    "        ylabel=f\"sum of {metrics_name}\",\n",
    "        #xlabel=\"frame\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    ax2.set(\n",
    "        ylabel=f\"sum of {metrics_name}\",\n",
    "        #xlabel=\"frame\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    ax3.set(\n",
    "        ylabel=metrics_name,\n",
    "        #xlabel=\"frame\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    ax4.set(\n",
    "        ylabel=metrics_name,\n",
    "        #xlabel=\"frame\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    ax5.set(\n",
    "        ylabel=\"normalised values (ratio)\",\n",
    "        xlabel=\"frame\",\n",
    "        ylim=[1/ylimit_log,ylimit_log],\n",
    "    )\n",
    "    ax6.set(\n",
    "        ylabel=\"normalised values (ratio)\",\n",
    "        xlabel=\"frame\",\n",
    "        ylim=[1/ylimit_log,ylimit_log],\n",
    "    )\n",
    "    # if metrics_name=='velocity':\n",
    "    #     peak_distribution_stationary=np.argmax(p3, axis=1)\n",
    "    #     peak_distribution_move=np.argmax(p4, axis=1)\n",
    "    # else:\n",
    "    peak_distribution_stationary=np.argmax(abs(p3[:,number_frame_scene_changing:]), axis=1)\n",
    "    peak_distribution_stationary=peak_distribution_stationary+number_frame_scene_changing\n",
    "    \n",
    "    #print(peak_distribution_stationary.shape)\n",
    "    peak_distribution_move=np.argmax(abs(p4[:,number_frame_scene_changing:]), axis=1)\n",
    "    peak_distribution_move=peak_distribution_move+number_frame_scene_changing\n",
    "    #print(peak_distribution_move.shape)\n",
    "    ax5.set_yscale('log')\n",
    "    ax6.set_yscale('log')\n",
    "    ax7.hist(abs(peak_distribution_stationary),bins=10)\n",
    "    ax7.set(\n",
    "        ylabel=\"Count\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    ax8.hist(abs(peak_distribution_move),bins=10)\n",
    "    ax8.set(\n",
    "        ylabel=\"Count\",\n",
    "        xlabel=\"Time (s)\",\n",
    "        xticks=[0,analysis_window[1]*monitor_fps],\n",
    "        xticklabels=(['0', str(analysis_window[1])]),\n",
    "    )\n",
    "    if save_output==True:\n",
    "        fig_name=f\"ts_plot_{exp_name}_{variable_values}_{metrics_name}_{type_key}.png\"\n",
    "        fig.savefig(fig_name)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_methods.update({\"save_output\": True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "##criteria to sort out data from good followers\n",
    "# fair_follower_threshold=0.20\n",
    "#fair_follower_threshold=0.1667\n",
    "# good_follower_threshold=0.3\n",
    "#data from grass1 secene\n",
    "# first 1/3 best of followers: 0.04542389707280022\n",
    "# middle 1/3 best of followers: 0.024365961537554066\n",
    "use_follower_threshold=True\n",
    "if use_follower_threshold==True:\n",
    "    p_follow=all_evaluation.groupby(['animal_id'])['num_follow_epochs'].sum()/all_evaluation.groupby(['animal_id'])['number_frames'].sum()\n",
    "    #follower_of_interest=p_follow>good_follower_threshold\n",
    "    #follower_of_interest=(p_follow>fair_follower_threshold) & (p_follow<good_follower_threshold)\n",
    "    follower_of_interest=p_follow>0.025\n",
    "    #rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "else:\n",
    "    #good_trial_threshold=0.125\n",
    "    good_trial_threshold=0.125\n",
    "    long_following_epochs=all_evaluation['num_follow_epochs']/all_evaluation['number_frames']>good_trial_threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if all_evaluation['object'].unique().shape[0]>1:\n",
    "for type_key in all_evaluation['radial_distance'].unique():\n",
    "    print(type_key)\n",
    "    rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "    for this_condition in all_evaluation['mu'].unique():\n",
    "        print(this_condition)\n",
    "        if use_follower_threshold==True:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['radial_distance'].reset_index(drop=True)==type_key) &(rows_of_follower.reset_index(drop=True))\n",
    "        else:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['radial_distance'].reset_index(drop=True)==type_key)&(long_following_epochs.reset_index(drop=True))\n",
    "        plot_visual_evoked_behaviour(these_metrics,these_normalised_metrics,after_movement_ith_trial,after_no_movement_ith_trial,analysis_methods,metrics_name,row_of_interest,type_key,this_condition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if all_evaluation['object'].unique().shape[0]>1:\n",
    "for type_key in all_evaluation['object'].unique():\n",
    "    print(type_key)\n",
    "    rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "    for this_condition in all_evaluation['mu'].unique():\n",
    "        print(this_condition)\n",
    "        if use_follower_threshold==True:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['object'].reset_index(drop=True)==type_key) &(rows_of_follower.reset_index(drop=True))\n",
    "        else:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['object'].reset_index(drop=True)==type_key)&(long_following_epochs.reset_index(drop=True))\n",
    "        plot_visual_evoked_behaviour(these_metrics,these_normalised_metrics,after_movement_ith_trial,after_no_movement_ith_trial,analysis_methods,metrics_name,row_of_interest,type_key,this_condition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movement_trial_boolean,these_metrics,these_normalised_metrics=split_trials(analysis_methods,all_trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics_name='velocity'\n",
    "for type_key in all_evaluation['object'].unique():\n",
    "    print(type_key)\n",
    "    rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "    for this_condition in all_evaluation['mu'].unique():\n",
    "        print(this_condition)\n",
    "        if use_follower_threshold==True:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['object'].reset_index(drop=True)==type_key)&(rows_of_follower.reset_index(drop=True))\n",
    "        else:\n",
    "            row_of_interest=(all_evaluation['mu'].reset_index(drop=True)==this_condition)&(all_evaluation['object'].reset_index(drop=True)==type_key)&(long_following_epochs.reset_index(drop=True))\n",
    "        plot_visual_evoked_behaviour(these_metrics,these_normalised_metrics,after_movement_ith_trial,after_no_movement_ith_trial,analysis_methods,metrics_name,row_of_interest,type_key,this_condition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics_name='velocity'\n",
    "for type_key in all_evaluation['object'].unique():\n",
    "    print(type_key)\n",
    "    if use_follower_threshold==True:\n",
    "        rows_of_follower=follower_of_interest.repeat(int(all_evaluation.shape[0]/follower_of_interest.shape[0]))\n",
    "        row_of_interest=(all_evaluation['object'].reset_index(drop=True)==type_key)&rows_of_follower.reset_index(drop=True)\n",
    "    else:\n",
    "        row_of_interest=(all_evaluation['object'].reset_index(drop=True)==type_key)&(long_following_epochs.reset_index(drop=True))\n",
    "    #print(row_of_interest[row_of_interest==True].index.to_list())\n",
    "    plot_visual_evoked_behaviour(these_metrics,these_normalised_metrics,after_movement_ith_trial,after_no_movement_ith_trial,analysis_methods,'velocity',row_of_interest,type_key)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "unity_analysis",
   "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.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}