Projects/ABM_DA/experiments/grand_central_terminal_data/Real_Data_correct_trails.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Grand Central Terminal - Trajectories\n",
    "    author: P. Ternes\n",
    "    created: 31/05/2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data assimilation experiments using real data are conducted using information from a video sequence, recorded at Grand Central Terminal (GCT). The original video is a 33:20-minutes video sequence shoot at 25fps with a resolution of 720x480.\n",
    "\n",
    "The Grand Terminal data is available on the website: http://www.ee.cuhk.edu.hk/~xgwang/grandcentral.html\n",
    "\n",
    "A large analysis of the real data can be found at [`Grand Central Terminal Callibration notebook`](gcs_callibration.ipynb).\n",
    "\n",
    "A explanation of the StationSim new version can be found at [`Grand Central Terminal Model notebook`](StationSim_GrandCentral_version.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The trajectories were determined through a KLT keypoint tracker on GCT video. Many tracking failures occur due to frequent occlusions in crowded scenes, and most trajectories are highly fragmented with large portions of missing observations.\n",
    "\n",
    "When we use the particle filter, this fragmentation can be a big problem. \n",
    "\n",
    "**In this notebook we have corrected the fragmentation of trajectories.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading original data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import loadmat\n",
    "\n",
    "#open and read original file\n",
    "file_name = \"gcs_real_trajectories.mat\"\n",
    "file = loadmat(file_name)\n",
    "trajectories = file[\"trks\"][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Organize the data\n",
    "\n",
    "To facilitate future analysis, the trajectories will be reorganized in a list called \"data\". Each line in this list will have 4 elements: pedestrian ID, position $x$, position $y$, and time. The positions $x$ and $y$ are  undistorted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "for i in range (len(trajectories)):\n",
    "    x, y, t = trajectories[i]\n",
    "    x, y, t = x.flatten(), y.flatten(), t.flatten()\n",
    "    for j in range (len(x)):\n",
    "        values = [i, x[j], y[j], t[j]]  # indices: ['pedestrianID', 'x', 'y', 'time']\n",
    "        data.append(values)\n",
    "\n",
    "        \n",
    "from operator import itemgetter\n",
    "data3 = []\n",
    "for i in range (len(data)):\n",
    "    data[i][2] = -1.* data[i][2] + 455.   # correct upside down\n",
    "    data3.append(data[i])\n",
    "    ratio = data[i][2]/431\n",
    "    data3[i][1] = data[i][1]-347.5 #centre x axis to camera\n",
    "    data3[i][1] = (data3[i][1] * (1 + ratio)) + 355.1403712296984 # new x position\n",
    "    data3[i][2] = data[i][2] * (1 + ratio)  # new y position\n",
    "\n",
    "#sort by frame\n",
    "data1 = sorted(data3, key=itemgetter(3))\n",
    "\n",
    "#sort by pedestrian\n",
    "data2 = sorted(data3, key=itemgetter(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select 3000 frames from the frame 20000 to 23000.\n",
    "\n",
    "Just data from this frame interval are analysed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from operator import itemgetter\n",
    "data4 = []\n",
    "for i in range (len(data3)):\n",
    "    if (data3[i][3] >= 20000 and data3[i][3] < 23000):\n",
    "        data4.append(data3[i])\n",
    "        \n",
    "#sort by frame\n",
    "data1 = sorted(data4, key=itemgetter(3))\n",
    "\n",
    "#sort by pedestrian\n",
    "data2 = sorted(data4, key=itemgetter(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the data from 3000 frames in files\n",
    "\n",
    "### By agent\n",
    "\n",
    "Each file conteins information about one agent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "directory = 'agents'\n",
    "if not(os.path.exists(directory)):\n",
    "    os.mkdir(directory)\n",
    "\n",
    "agent = data2[0][0]\n",
    "save_file = open(directory+'/agent_'+ str(agent) +'.dat', 'w')\n",
    "print('# time', 'x', 'y', file=save_file)\n",
    "for i in range (len(data2)):\n",
    "    if (data2[i][0] == agent):\n",
    "        print(data2[i][3], data2[i][1], data2[i][2], file=save_file)\n",
    "    else:\n",
    "        save_file.close()\n",
    "        agent = data2[i][0]\n",
    "        save_file = open(directory+'/agent_'+ str(agent) +'.dat', 'w')\n",
    "        print(data2[i][3], data2[i][1], data2[i][2], file=save_file)\n",
    "save_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Remove some informations\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove pedestrians that never enter the environment\n",
    "\n",
    "Remove pedestrians that never enter the environment. It also moves pedestrians that have complete trajectories to a different folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "before remove/move, the number of file is:  1049\n",
      "agents/agent_21303.dat\n",
      "agents/agent_21405.dat\n",
      "agents/agent_21408.dat\n",
      "agents/agent_21414.dat\n",
      "agents/agent_21532.dat\n",
      "agents/agent_21559.dat\n",
      "agents/agent_21571.dat\n",
      "agents/agent_21575.dat\n",
      "agents/agent_21589.dat\n",
      "agents/agent_21592.dat\n",
      "agents/agent_21605.dat\n",
      "agents/agent_21618.dat\n",
      "agents/agent_21624.dat\n",
      "agents/agent_21625.dat\n",
      "agents/agent_21633.dat\n",
      "agents/agent_21638.dat\n",
      "agents/agent_21658.dat\n",
      "agents/agent_21675.dat\n",
      "agents/agent_21680.dat\n",
      "agents/agent_21691.dat\n",
      "agents/agent_21694.dat\n",
      "agents/agent_21695.dat\n",
      "agents/agent_21702.dat\n",
      "agents/agent_21738.dat\n",
      "agents/agent_21741.dat\n",
      "agents/agent_21742.dat\n",
      "agents/agent_21743.dat\n",
      "agents/agent_21744.dat\n",
      "agents/agent_21756.dat\n",
      "agents/agent_21758.dat\n",
      "agents/agent_21768.dat\n",
      "agents/agent_21780.dat\n",
      "agents/agent_21788.dat\n",
      "agents/agent_21831.dat\n",
      "agents/agent_21838.dat\n",
      "agents/agent_21840.dat\n",
      "agents/agent_21852.dat\n",
      "agents/agent_21854.dat\n",
      "agents/agent_21865.dat\n",
      "agents/agent_21866.dat\n",
      "agents/agent_21868.dat\n",
      "agents/agent_21870.dat\n",
      "agents/agent_21871.dat\n",
      "agents/agent_21877.dat\n",
      "agents/agent_21884.dat\n",
      "agents/agent_21885.dat\n",
      "agents/agent_21888.dat\n",
      "agents/agent_21900.dat\n",
      "agents/agent_21904.dat\n",
      "agents/agent_21919.dat\n",
      "agents/agent_21933.dat\n",
      "agents/agent_21950.dat\n",
      "agents/agent_21963.dat\n",
      "agents/agent_21989.dat\n",
      "agents/agent_21991.dat\n",
      "agents/agent_21996.dat\n",
      "agents/agent_21999.dat\n",
      "agents/agent_22013.dat\n",
      "agents/agent_22018.dat\n",
      "agents/agent_22033.dat\n",
      "agents/agent_22057.dat\n",
      "agents/agent_22067.dat\n",
      "agents/agent_22069.dat\n",
      "agents/agent_22097.dat\n",
      "agents/agent_22101.dat\n",
      "agents/agent_22113.dat\n",
      "agents/agent_22115.dat\n",
      "agents/agent_22116.dat\n",
      "agents/agent_22117.dat\n",
      "agents/agent_22134.dat\n",
      "agents/agent_22145.dat\n",
      "agents/agent_22178.dat\n",
      "agents/agent_22189.dat\n",
      "agents/agent_22201.dat\n",
      "agents/agent_22203.dat\n",
      "agents/agent_22230.dat\n",
      "agents/agent_22232.dat\n",
      "agents/agent_22237.dat\n",
      "agents/agent_22242.dat\n",
      "agents/agent_22243.dat\n",
      "agents/agent_22246.dat\n",
      "agents/agent_22257.dat\n",
      "agents/agent_22268.dat\n",
      "agents/agent_22285.dat\n",
      "agents/agent_22287.dat\n",
      "agents/agent_22291.dat\n",
      "agents/agent_22306.dat\n",
      "agents/agent_22319.dat\n",
      "agents/agent_22320.dat\n",
      "agents/agent_22326.dat\n",
      "agents/agent_22330.dat\n",
      "agents/agent_22332.dat\n",
      "agents/agent_22338.dat\n",
      "agents/agent_22345.dat\n",
      "agents/agent_22350.dat\n",
      "agents/agent_22360.dat\n",
      "agents/agent_22365.dat\n",
      "agents/agent_22408.dat\n",
      "agents/agent_22433.dat\n",
      "agents/agent_22434.dat\n",
      "agents/agent_22435.dat\n",
      "agents/agent_22442.dat\n",
      "agents/agent_22446.dat\n",
      "agents/agent_22449.dat\n",
      "agents/agent_22461.dat\n",
      "agents/agent_22509.dat\n",
      "agents/agent_22542.dat\n",
      "agents/agent_22550.dat\n",
      "agents/agent_22555.dat\n",
      "agents/agent_22557.dat\n",
      "After remove/move, the number of file is:  919\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import shutil\n",
    "    \n",
    "def is_within_bounds(boundaries, loc):\n",
    "    return all(boundaries[0] < loc) and\\\n",
    "           all(loc < boundaries[1])\n",
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]\n",
    "    deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "boundaries = np.array([[0, 0], [740, 700]])\n",
    "agentsCompleteTrails = []\n",
    "\n",
    "directory = 'agents/'\n",
    "files = os.listdir(directory)\n",
    "number_files = len(files)\n",
    "print('before remove/move, the number of file is: ', number_files)\n",
    "\n",
    "directory1 = 'agents_complete_trails/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "\n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    locI = [x[0], y[0]]\n",
    "    if not (is_within_bounds(boundaries, locI)):\n",
    "        locF = [x[-1], y[-1]]\n",
    "        if not (is_within_bounds(boundaries, locF)):\n",
    "            dist = normaldist(locI, locF)\n",
    "            if (dist<130.): # 130 defined by eye\n",
    "                # Remove agents that never enter the environment.\n",
    "                os.remove(file_name)\n",
    "                print(file_name)\n",
    "            else:\n",
    "                # moving agents that have complete trajectory to another folder\n",
    "                shutil.move(file_name, directory1 + agent_file)\n",
    "\n",
    "files = os.listdir(directory)\n",
    "number_files = len(files)\n",
    "print('After remove/move, the number of file is: ', number_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove out-of-boundaries positions\n",
    "\n",
    "Rewrite the position files with only the positions within the station limits.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def is_within_bounds(boundaries, loc):\n",
    "    return all(boundaries[0] < loc) and\\\n",
    "           all(loc < boundaries[1])\n",
    "\n",
    "boundaries = np.array([[0, 0], [740, 700]])\n",
    "agentsCompleteTrails = []\n",
    "\n",
    "directory = 'agents/'\n",
    "files = os.listdir(directory)\n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    save_file = open(directory + agent_file, 'w')\n",
    "    print('# time', 'x', 'y', file=save_file)\n",
    "    for i in range(len(x)):\n",
    "        loc = [x[i], y[i]]\n",
    "        if (is_within_bounds(boundaries, loc)):\n",
    "            print(t[i], x[i], y[i], file=save_file)\n",
    "    save_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove pedestrians with very short time life\n",
    "\n",
    "If the pedestrian exists only for 1 second (25 frames) or less, remove!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "agents/agent_21502.dat\n",
      "agents/agent_21403.dat\n",
      "agents/agent_21517.dat\n",
      "agents/agent_21464.dat\n",
      "agents/agent_21521.dat\n",
      "agents/agent_21546.dat\n",
      "agents/agent_21465.dat\n",
      "agents/agent_21561.dat\n",
      "agents/agent_21581.dat\n",
      "agents/agent_21585.dat\n",
      "agents/agent_21586.dat\n",
      "agents/agent_21484.dat\n",
      "agents/agent_21591.dat\n",
      "agents/agent_21603.dat\n",
      "agents/agent_21609.dat\n",
      "agents/agent_21643.dat\n",
      "agents/agent_21665.dat\n",
      "agents/agent_21679.dat\n",
      "agents/agent_21689.dat\n",
      "agents/agent_21724.dat\n",
      "agents/agent_21726.dat\n",
      "agents/agent_21733.dat\n",
      "agents/agent_21747.dat\n",
      "agents/agent_21784.dat\n",
      "agents/agent_21789.dat\n",
      "agents/agent_21791.dat\n",
      "agents/agent_21797.dat\n",
      "agents/agent_21812.dat\n",
      "agents/agent_21818.dat\n",
      "agents/agent_21824.dat\n",
      "agents/agent_21859.dat\n",
      "agents/agent_21862.dat\n",
      "agents/agent_21864.dat\n",
      "agents/agent_21892.dat\n",
      "agents/agent_21911.dat\n",
      "agents/agent_21856.dat\n",
      "agents/agent_21912.dat\n",
      "agents/agent_21914.dat\n",
      "agents/agent_21935.dat\n",
      "agents/agent_21941.dat\n",
      "agents/agent_21948.dat\n",
      "agents/agent_21956.dat\n",
      "agents/agent_21990.dat\n",
      "agents/agent_21994.dat\n",
      "agents/agent_22011.dat\n",
      "agents/agent_22024.dat\n",
      "agents/agent_22027.dat\n",
      "agents/agent_22051.dat\n",
      "agents/agent_22052.dat\n",
      "agents/agent_21975.dat\n",
      "agents/agent_22058.dat\n",
      "agents/agent_22062.dat\n",
      "agents/agent_22074.dat\n",
      "agents/agent_21978.dat\n",
      "agents/agent_22075.dat\n",
      "agents/agent_22154.dat\n",
      "agents/agent_22155.dat\n",
      "agents/agent_22161.dat\n",
      "agents/agent_22138.dat\n",
      "agents/agent_22166.dat\n",
      "agents/agent_22179.dat\n",
      "agents/agent_22182.dat\n",
      "agents/agent_22184.dat\n",
      "agents/agent_22191.dat\n",
      "agents/agent_22199.dat\n",
      "agents/agent_22222.dat\n",
      "agents/agent_22223.dat\n",
      "agents/agent_22133.dat\n",
      "agents/agent_22150.dat\n",
      "agents/agent_22224.dat\n",
      "agents/agent_22228.dat\n",
      "agents/agent_22256.dat\n",
      "agents/agent_22262.dat\n",
      "agents/agent_22284.dat\n",
      "agents/agent_22295.dat\n",
      "agents/agent_22299.dat\n",
      "agents/agent_22313.dat\n",
      "agents/agent_22340.dat\n",
      "agents/agent_22375.dat\n",
      "agents/agent_22387.dat\n",
      "agents/agent_22398.dat\n",
      "agents/agent_22270.dat\n",
      "agents/agent_22399.dat\n",
      "agents/agent_22404.dat\n",
      "agents/agent_22432.dat\n",
      "agents/agent_22457.dat\n",
      "agents/agent_22473.dat\n",
      "agents/agent_22475.dat\n",
      "agents/agent_22496.dat\n",
      "agents/agent_22498.dat\n",
      "agents/agent_22499.dat\n",
      "agents/agent_22436.dat\n",
      "agents/agent_22501.dat\n",
      "agents/agent_22503.dat\n",
      "agents/agent_22514.dat\n",
      "agents/agent_22516.dat\n",
      "agents/agent_22520.dat\n",
      "agents/agent_22522.dat\n",
      "agents/agent_22527.dat\n",
      "agents/agent_22530.dat\n",
      "agents/agent_22451.dat\n",
      "agents/agent_22531.dat\n",
      "agents/agent_22540.dat\n",
      "agents/agent_22582.dat\n",
      "agents/agent_22583.dat\n",
      "agents/agent_22559.dat\n",
      "agents/agent_22576.dat\n",
      "agents/agent_22579.dat\n",
      "agents/agent_22580.dat\n",
      "agents/agent_22581.dat\n",
      "agents/agent_22584.dat\n"
     ]
    }
   ],
   "source": [
    "directory = 'agents/' \n",
    "files = os.listdir(directory)\n",
    "\n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    try:\n",
    "        if(len(t)<25):\n",
    "            os.remove(file_name)\n",
    "            print(file_name)\n",
    "    except TypeError:\n",
    "        '''\n",
    "        This error occurs when len(t)=1. In\n",
    "        this case, the data is read as a float instead of an\n",
    "        array.\n",
    "        '''\n",
    "        print(file_name)\n",
    "        os.remove(file_name)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Put together some trajectories\n",
    "\n",
    "If a pedestrian has his/her trajectory closed in the middle of the station, it is likely that the rest of his/her trajectory has been associated with one or more pedestrians with a different ID. To identify which trajectories represent the movement of pedestrian $i$, we will initially look for trajectories that have started at most <font color=\"red\"><b> 25 frames</b></font> (1 seconds) after the trajectory of pedestrian $i$ has ended, within a radius of <font color=\"red\"><b> 25 * 0.87 pixels</b></font>  from its final position. The radius is associated with the mean pedestrian velocity determined over all trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import shutil\n",
    "\n",
    "def is_within_bounds(boundaries, loc):\n",
    "    return all(boundaries[0] < loc) and\\\n",
    "           all(loc < boundaries[1])\n",
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]\n",
    "    deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "boundaries = np.array([[0, 0], [740, 700]])\n",
    "\n",
    "directory = 'agents/'\n",
    "files = os.listdir(directory)\n",
    "number_files = len(files)\n",
    "print(number_files)\n",
    "\n",
    "time_REF = 1 # time in seconds!\n",
    "radius_REF = time_REF * 25 * 0.8719127839925644\n",
    "\n",
    "for agent_file1 in files:\n",
    "    try:\n",
    "        file_name1 = directory + agent_file1\n",
    "        t1, x1, y1 = np.loadtxt(file_name1,unpack=True)\n",
    "    except OSError:\n",
    "        pass      \n",
    "    for agent_file2 in files:\n",
    "        file_name2 = directory + agent_file2\n",
    "        try:\n",
    "            t2, x2, y2 = np.loadtxt(file_name2,unpack=True)\n",
    "            deltat = t2[0] - t1[-1]\n",
    "            if(deltat>0 and deltat<time_REF):\n",
    "                dist = normaldist([x1[-1], y1[-1]], [x2[0], y2[0]])\n",
    "                if (dist<radius_REF):\n",
    "                    save_file = open(file_name1, 'w')\n",
    "                    # save the file2 inside the file1\n",
    "                    for i in range (len(t1)):\n",
    "                        print(t1[i], x1[i], y1[i], file=save_file)\n",
    "                    for i in range (len(t2)):\n",
    "                        print(t2[i], x2[i], y2[i], file=save_file)\n",
    "                    save_file.close()\n",
    "                    # remove the file2\n",
    "                    os.remove(file_name2)\n",
    "                    print(file_name1, file_name2)\n",
    "        except OSError:\n",
    "            pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Move pedestrians that have complete trajectories\n",
    "\n",
    "To prevent finished trajectories from being associated with new trajectories, move pedestrians who have trajectory to another folder (**agents_complete_trails/**). By doing this, the loop from the previous step will also become smaller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import shutil\n",
    "    \n",
    "def is_within_bounds(boundaries, loc):\n",
    "    return all(boundaries[0] < loc) and\\\n",
    "           all(loc < boundaries[1])\n",
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]\n",
    "    deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "boundaries = np.array([[20, 20], [720, 680]])\n",
    "agentsCompleteTrails = []\n",
    "\n",
    "directory = 'agents/'\n",
    "files = os.listdir(directory)\n",
    "number_files = len(files)\n",
    "print('before remove/move, the number of file is: ', number_files)\n",
    "\n",
    "directory1 = 'agents_complete_trails/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "\n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    locI = [x[0], y[0]]\n",
    "    if not (is_within_bounds(boundaries, locI)):\n",
    "        locF = [x[-1], y[-1]]\n",
    "        if not (is_within_bounds(boundaries, locF)):\n",
    "            plt.plot(x,y)\n",
    "            print(file_name)\n",
    "            # moving agents that have complete trajectory to another folder\n",
    "            shutil.move(file_name, directory1 + agent_file)\n",
    "\n",
    "files = os.listdir(directory)\n",
    "number_files = len(files)\n",
    "print('After remove/move, the number of file is: ', number_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Repeat steps 2 and 3\n",
    "\n",
    "As we are deleting files by running the loops, some comparisons between trajectories are lost during the execution of the loop. To ensure that all possible trajectories are grouped, repeat steps 2 and 3 until no file is changed. It is a manual and laborious approach, but suitable for that moment.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Plot the complete trajectories to check\n",
    "\n",
    "The approach taken in step 2 to join pedestrian with different IDs is relatively simple. In this approach, the possibility of more than one pedestrian being started in the same radius is not taken into account and no comparative analysis is performed. In a preliminary analysis we found that the cases where the trajectories are selected incorrectly represent a small percentage. However, it is recommended that the new trajectories are checked systematically to avoid errors. To do that, plot the trajectories and look the results.\n",
    "\n",
    "If something looks wrong, split the file into two files and move them to the 'agents' folder so they can be processed again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_station():\n",
    "    # plotting the station walls\n",
    "    a = [0, 0, 740, 740, 0]; b = [-190, 700, 700, -190, -190];\n",
    "    a = [0, 0, 740, 740, 0]; b = [0, 700, 700, 0, 0]; plt.plot(a,b,'-', linewidth=1,color='k')\n",
    "\n",
    "    # plotting the station gate 0\n",
    "    a = [-150, 0, 0, -150, -150]; b = [400, 400, 150, 150, 400];\n",
    "    plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(-84, 252, '0', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 1\n",
    "    a = [0, 250, 250, 0, 0]; b = [850, 850, 700, 700, 850];\n",
    "    plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(112, 770, '1', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 2\n",
    "    a = [455, 700, 700, 455, 455]; b = [850, 850, 700, 700, 850];\n",
    "    plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(560, 770, '2', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 3\n",
    "    a = [740, 860, 860, 740, 740];\n",
    "    b = [700, 700, 610, 610, 700]; plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(784, 630, '3', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 4\n",
    "    a = [740, 860, 860, 740, 740]; b = [550, 550, 400, 400, 550];\n",
    "    plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(784, 455.0, '4', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 5\n",
    "    a = [740, 860, 860, 740, 740]; b = [340, 340, 190, 190, 340];\n",
    "    plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(784, 252, '5', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 6\n",
    "    a = [740, 860, 860, 740, 740];\n",
    "    b = [130, 130, 0, 0, 130]; plt.plot(a, b,'-',linewidth=2,color='b'); plt.text(784, 42,  '6', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 7\n",
    "    a = [555, 740, 740, 555, 555]; b = [0, 0, -70, -70, 0]; plt.plot(a, b, linewidth=2,color='b')\n",
    "    plt.text(630, -49.0, '7', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 8\n",
    "    a = [370, 555, 556, 370, 370]; b = [0, 0, -70, -70, 0]; plt.plot(a, b, linewidth=2,color='b')\n",
    "    plt.text(448,-49, '8', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 9\n",
    "    a = [185, 370, 370, 185, 185]; b = [0, 0, -70, -70, 0]; plt.plot(a, b, linewidth=2,color='b')\n",
    "    plt.text(266, -49.0, '9', fontsize=12)\n",
    "\n",
    "    # plotting the station gate 10\n",
    "    a = [0, 185, 185, 0, 0]; b = [0, 0, -70, -70, 0]; plt.plot(a, b, linewidth=2,color='b')\n",
    "    plt.text(63.0, -49.0, '10', fontsize=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [],
   "source": [
    "directory = 'agents_complete_trails/'\n",
    "\n",
    "files = os.listdir(directory)\n",
    "\n",
    "directory1 = 'trails/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "    \n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    fig, ax = plt.subplots(num=None,figsize=(7.4, 7),dpi=128,facecolor='w',edgecolor='k')\n",
    "    plot_station()\n",
    "    plt.plot(x,y, 'k')\n",
    "    plt.plot(x[0],y[0], 'o', ms=10, color='r')\n",
    "    \n",
    "    file_name =  directory1 + agent_file + '.png'\n",
    "    plt.savefig(file_name)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Repeat steps 2, 3, 4 and 5.\n",
    "\n",
    "After performing the steps above until all possible files are changed and moved to another folder, it is necessary change the radius of comparison. Increase the time by 1 second and repeat steps 2 to 5. Then increase the time again by 1 second and repeat the entire procedure. Do this until all trajectories are selected.\n",
    "\n",
    "It is important to note that as the radius increases, the probability of selecting an incorrect trajectory increases. The bigger the radius, the more important step 5. On the other hand, the number of files will be smaller for the tests using larger radius, which will facilitate the manual part of the process\n",
    "\n",
    "It is also important to note that in the end, there will be some short trajectories that could not be related to any pedestrian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7: Gates in and out\n",
    "\n",
    "To determine the starting and ending gate for each pedestrian, we will calculate the distance between the starting/ending position and the gates. The minimum distance is related to the desired gate. The more complete the trajectory, the better the gate selection. However, not all trajectories are complete. Therefore, after determining the initial and final gates using this code, check the figures in the **trails folder**. If any gate was determined incorrectly, open the **gates.dat** file and manually change the gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]\n",
    "    deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "# gates locations!\n",
    "gates_locations =\\\n",
    "    np.array([[0, 275],  # gate 0\n",
    "              [125, 700],   # gate 1\n",
    "              [577.5 , 700],  # gate 2\n",
    "              [740, 655],  # gate 3\n",
    "              [740, 475],  # gate 4\n",
    "              [740, 265],   # gate 5\n",
    "              [740, 65],   # gate 6\n",
    "              [647.5, 0],   # gate 7\n",
    "              [462.5, 0],   # gate 8\n",
    "              [277.5, 0],   # gate 9\n",
    "              [92.5, 0]])   # gate 10\n",
    "\n",
    "n_gates = len(gates_locations)\n",
    "gates = []\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "directory = 'agents_complete_trails/'\n",
    "files = os.listdir(directory)\n",
    "directory1 = 'trails/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "save_file = open('gates.dat', 'w')\n",
    "yy = []\n",
    "j = 0\n",
    "for agent_file in files:\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    fig, ax = plt.subplots(num=None,figsize=(7.4, 7),dpi=128,facecolor='w',edgecolor='k')\n",
    "    plot_station()\n",
    "    plt.plot(x,y, 'k')\n",
    "    plt.plot(x[0],y[0], 'o', ms=10, color='r')\n",
    "    \n",
    "    gate_in = 0\n",
    "    gate_out = 0\n",
    "    min_dist_in = 700\n",
    "    min_dist_out = 700\n",
    "    for i in range(n_gates):\n",
    "        dist_in = normaldist(gates_locations[i], [x[0], y[0]])\n",
    "        dist_out = normaldist(gates_locations[i], [x[-1], y[-1]])\n",
    "        if (dist_in < min_dist_in):\n",
    "            gate_in = i\n",
    "            min_dist_in = dist_in\n",
    "        if (dist_out < min_dist_out):\n",
    "            gate_out = i\n",
    "            min_dist_out = dist_out\n",
    "    print(j, gate_in, gate_out, file=save_file)        \n",
    "    plt.plot(gates_locations[gate_in][0], gates_locations[gate_in][1], 'o', ms=10, color='b')\n",
    "    plt.plot(gates_locations[gate_out][0], gates_locations[gate_out][1], 's', ms=10, color='b')\n",
    "    file_name =  directory1 +str(j)+ '.png'\n",
    "    plt.savefig(file_name)\n",
    "    plt.close()\n",
    "    j += 1\n",
    "save_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 8: Data correction\n",
    "\n",
    "It is necessary to make 3 types of corrections in the data:\n",
    "\n",
    "* Delete data outside station at the entrance;\n",
    "* Delete data outside station at the exit;\n",
    "* Complete the trajectory at the intervals that are missing;\n",
    "\n",
    "Final data will be saved in the pedestrian_data folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "gates_locations = np.array([[0, 275], [125, 700], [577.5 , 700], [740, 655], [740, 475], [740, 265], [740, 65], \n",
    "                            [647.5, 0], [462.5, 0], [277.5, 0], [92.5, 0]])\n",
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]; deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "def is_inside(gate, x, y, tolerance):\n",
    "        if(gates_locations[gate][0] == 0):\n",
    "            if(x <= 0+tolerance):\n",
    "                return False\n",
    "            else:\n",
    "                return True\n",
    "        elif(gates_locations[gate][0] == 740):\n",
    "            if(x>=740-tolerance):\n",
    "                return False\n",
    "            else:\n",
    "                return True\n",
    "        elif(gates_locations[gate][1] == 0):\n",
    "            if(y <= 0+tolerance):\n",
    "                return False\n",
    "            else:\n",
    "                return True\n",
    "        else:\n",
    "            if(y >= 700-tolerance):\n",
    "                return False\n",
    "            else:\n",
    "                return True\n",
    "\n",
    "def new_position(position_desire, initial_x, initial_y, speed):\n",
    "    direction = (position_desire - [initial_x, initial_y])/\\\n",
    "                normaldist(gates_locations[_in], [initial_x, initial_y])\n",
    "    noise = np.random.normal(2.5, 0.5, [2])\n",
    "    initial_x = initial_x + noise[0]*direction[0]*speed # timestep = 1.0\n",
    "    initial_y = initial_y + noise[1]*direction[1]*speed # timestep = 1.0\n",
    "    return initial_x, initial_y\n",
    "    \n",
    "    \n",
    "directory = 'agents_complete_trails/'\n",
    "files = os.listdir(directory)\n",
    "\n",
    "directory1 = 'pedestrian_data/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "\n",
    "id_, gate_in, gate_out  = np.loadtxt('gates.dat',unpack=True)\n",
    "\n",
    "i = 0\n",
    "for agent_file in files:\n",
    "    # open trajectory file\n",
    "    file_name = directory + agent_file\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    \n",
    "    _in = int(gate_in[i]); _out = int(gate_out[i])\n",
    "    # excluding values outside station\n",
    "    while not (is_inside(_in, x[0], y[0], 1)):\n",
    "        x = np.delete(x, [0]); y = np.delete(y, [0]); t = np.delete(t, [0])\n",
    "    while not (is_inside(_out, x[len(x)-1], y[len(y)-1], 1)):\n",
    "        x = np.delete(x, [len(x)-1]); y = np.delete(y, [len(y)-1]); t = np.delete(t, [len(t)-1])\n",
    "\n",
    "    # getting some usefull informations\n",
    "    initial_time = t[0]; initial_x = x[0]; initial_y = y[0]\n",
    "    final_time = t[-1]; final_x = x[-1]; final_y = y[-1]\n",
    "    initial_new_lines = []; final_new_lines = []\n",
    "    speed = normaldist([x[0], y[0]], [x[-1], y[-1]]) / (t[-1] - t[0])\n",
    "\n",
    "    # open a new file to store the complete trails.\n",
    "    save_file = open(directory1 + 'agent_' + str(i)+'.dat', 'w')\n",
    "\n",
    "    # estimating the beginning of the trajectory\n",
    "    while(is_inside(_in, initial_x, initial_y, 1.05*7)):\n",
    "        initial_x, initial_y = new_position(gates_locations[_in], initial_x, initial_y, speed)\n",
    "        initial_time -= 1.0;\n",
    "        initial_new_lines.insert(0,[initial_time, initial_x, initial_y])\n",
    "\n",
    "    # saving the beginning of the trajectory\n",
    "    for line in initial_new_lines:\n",
    "        print(line[0], line[1], line[2], file=save_file)\n",
    "\n",
    "    # estimating and saving the intermediate parts of the trajectory\n",
    "    tpass= t[0]-1\n",
    "    for j in range(len(t)-1):\n",
    "        if(t[j]>tpass+1):\n",
    "            initial_x = x[j-1]\n",
    "            initial_y = y[j-1]\n",
    "            speed1 = normaldist([x[j-1], y[j-1]], [x[j], y[j]]) / (t[j] - t[j-1])\n",
    "            for k in range (int(tpass), int(t[j]), 1):\n",
    "                tpass+=1\n",
    "                initial_x, initial_y = new_position(np.array([x[j], y[j]]), initial_x, initial_y, speed1)\n",
    "                print(tpass, initial_x, initial_y, file=save_file)\n",
    "        else:\n",
    "            tpass+=1\n",
    "            print(t[j], x[j], y[j], file=save_file)\n",
    "            \n",
    "    # estimating the end of the trajectory\n",
    "    while(is_inside(_out, final_x, final_y, 1.05*7)):\n",
    "        final_x, final_y = new_position(gates_locations[_out], final_x, final_y, speed)\n",
    "        final_time += 1.0; #initial_x = location[0]; initial_y = location[1]\n",
    "        final_new_lines.append([final_time, final_x, final_y])\n",
    "\n",
    "    # saving the end of the trajectory\n",
    "    for line in final_new_lines:\n",
    "        print(line[0], line[1], line[2], file=save_file)\n",
    "    save_file.close()\n",
    "\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Step: Organize the data to use in the StationSim_gcs model\n",
    "\n",
    "Create a file called **activation.dat** and put the pedestrian ID, time activation, gate in, gate out, and speed informations. Also creates files separeted by frames, where each frame file has the pedestrian positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "def normaldist(locI, locF):\n",
    "    deltaX = locI[0] - locF[0]\n",
    "    deltaY = locI[1] - locF[1]      \n",
    "    return np.sqrt(deltaX*deltaX + deltaY*deltaY)\n",
    "\n",
    "directory = 'pedestrian_data/'\n",
    "files = os.listdir(directory)\n",
    "n_files = len(files)\n",
    "id_, gate_in, gate_out  = np.loadtxt('gates.dat',unpack=True)\n",
    "directory1 = 'gcs_final_real_data/'\n",
    "if not(os.path.exists(directory1)):\n",
    "    os.mkdir(directory1)\n",
    "\n",
    "save_file1 = open(directory1 + 'activation.dat', 'w')\n",
    "print('#', 'pedestrianID', 'time_activation', 'gate_in', 'gate_out', 'speed', file=save_file1)\n",
    "data = []\n",
    "for i in range(n_files):\n",
    "    # open trajectory file\n",
    "    file_name = directory + 'agent_' + str(i) + '.dat'\n",
    "    t, x, y = np.loadtxt(file_name,unpack=True)\n",
    "    for j in range (len(x)):\n",
    "        values = [int(i), x[j], y[j], t[j]-19620]  # indices: ['pedestrianID', 'x', 'y', 'time', gate_in, gate_out]\n",
    "        data.append(values)\n",
    "    speed = normaldist([x[0], y[0]], [x[-1], y[-1]]) / (t[-1] - t[0])\n",
    "    print(int(i), t[0]-19620, int(gate_in[i]), int(gate_out[i]), speed, file=save_file1)\n",
    "save_file1.close()\n",
    "\n",
    "from operator import itemgetter\n",
    "#sort by frame\n",
    "data1 = sorted(data, key=itemgetter(3))\n",
    "\n",
    "frame = data1[0][3]\n",
    "print(frame)\n",
    "save_file = open(directory1+'frame_'+ str(frame) +'.dat', 'w')\n",
    "print('#', 'pedestrianID', 'x', 'y', file=save_file)\n",
    "\n",
    "for i in range (len(data1)):\n",
    "    if (data1[i][3] == frame):\n",
    "        print(data1[i][0], data1[i][1], data1[i][2], file=save_file)\n",
    "    else:\n",
    "        save_file.close()\n",
    "        frame = data1[i][3]\n",
    "        save_file = open(directory1+'/frame_'+ str(frame) +'.dat', 'w')\n",
    "        print('#', 'pedestrianID', 'x', 'y', file=save_file)\n",
    "        print(data1[i][0], data1[i][1], data1[i][2], file=save_file)\n",
    "save_file.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove unnecessary files and folders\n",
    "\n",
    "For experiments using the particle filter only files in the **gcs_final_real_data** folder are needed. All other files and folders can be deleted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra-Step Re-write traces in CSV format\n",
    "\n",
    "It can be useful to save the traces as a csv file wht columns 'individual', 'time', 'x', 'y'. Do that now. Note that this relies on the files created earlier, stored in the [GCT_final_real_data](GCT_final_real_data) folder:\n",
    "    - `GCT_final_real_data/activation.dat` - the IDs for each individual and the times that they enter the scene\n",
    "    - `GCT_final_real_data/frames/frame_x.0.dat` - locations of individuals in each frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import re\n",
    "import tqdm  # for progress bar on for loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>#</th>\n",
       "      <th>pedestrianID</th>\n",
       "      <th>time_activation</th>\n",
       "      <th>gate_in</th>\n",
       "      <th>gate_out</th>\n",
       "      <th>speed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>2273.0</td>\n",
       "      <td>8</td>\n",
       "      <td>5</td>\n",
       "      <td>1.417886</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>2401.0</td>\n",
       "      <td>9</td>\n",
       "      <td>3</td>\n",
       "      <td>1.130044</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>2411.0</td>\n",
       "      <td>9</td>\n",
       "      <td>3</td>\n",
       "      <td>1.120553</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>2564.0</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>0.554143</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>2591.0</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>0.600241</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>269</th>\n",
       "      <td>269</td>\n",
       "      <td>3753.0</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "      <td>0.579463</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>270</th>\n",
       "      <td>270</td>\n",
       "      <td>3720.0</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>1.025881</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>271</th>\n",
       "      <td>271</td>\n",
       "      <td>3849.0</td>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "      <td>1.406983</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>272</th>\n",
       "      <td>272</td>\n",
       "      <td>3948.0</td>\n",
       "      <td>9</td>\n",
       "      <td>5</td>\n",
       "      <td>1.408405</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>273</th>\n",
       "      <td>273</td>\n",
       "      <td>3957.0</td>\n",
       "      <td>9</td>\n",
       "      <td>6</td>\n",
       "      <td>1.335276</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>274 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       #  pedestrianID  time_activation  gate_in  gate_out  speed\n",
       "0      0        2273.0                8        5  1.417886    NaN\n",
       "1      1        2401.0                9        3  1.130044    NaN\n",
       "2      2        2411.0                9        3  1.120553    NaN\n",
       "3      3        2564.0                9        1  0.554143    NaN\n",
       "4      4        2591.0                9        1  0.600241    NaN\n",
       "..   ...           ...              ...      ...       ...    ...\n",
       "269  269        3753.0                9        0  0.579463    NaN\n",
       "270  270        3720.0                8        0  1.025881    NaN\n",
       "271  271        3849.0                8        2  1.406983    NaN\n",
       "272  272        3948.0                9        5  1.408405    NaN\n",
       "273  273        3957.0                9        6  1.335276    NaN\n",
       "\n",
       "[274 rows x 6 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "activation = pd.read_csv(os.path.join(\"GCT_final_real_data\",\"activation.dat\"), sep=\" \")\n",
    "activation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reading files: 100%|██████████| 5688/5688 [00:01<00:00, 3734.83it/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>individual</th>\n",
       "      <th>time</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>233</td>\n",
       "      <td>0</td>\n",
       "      <td>130.2375</td>\n",
       "      <td>692.7607</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>233</td>\n",
       "      <td>1</td>\n",
       "      <td>130.3354</td>\n",
       "      <td>692.6439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>233</td>\n",
       "      <td>2</td>\n",
       "      <td>130.4934</td>\n",
       "      <td>692.4644</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>233</td>\n",
       "      <td>3</td>\n",
       "      <td>130.5529</td>\n",
       "      <td>692.3266</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>233</td>\n",
       "      <td>4</td>\n",
       "      <td>130.6455</td>\n",
       "      <td>692.1638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226385</th>\n",
       "      <td>246</td>\n",
       "      <td>5683</td>\n",
       "      <td>572.5485</td>\n",
       "      <td>692.1653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226386</th>\n",
       "      <td>246</td>\n",
       "      <td>5684</td>\n",
       "      <td>572.6494</td>\n",
       "      <td>692.3114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226387</th>\n",
       "      <td>246</td>\n",
       "      <td>5685</td>\n",
       "      <td>572.7741</td>\n",
       "      <td>692.3983</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226388</th>\n",
       "      <td>246</td>\n",
       "      <td>5686</td>\n",
       "      <td>572.8418</td>\n",
       "      <td>692.4915</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226389</th>\n",
       "      <td>246</td>\n",
       "      <td>5687</td>\n",
       "      <td>573.0393</td>\n",
       "      <td>692.8700</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>226390 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        individual  time         x         y\n",
       "0              233     0  130.2375  692.7607\n",
       "1              233     1  130.3354  692.6439\n",
       "2              233     2  130.4934  692.4644\n",
       "3              233     3  130.5529  692.3266\n",
       "4              233     4  130.6455  692.1638\n",
       "...            ...   ...       ...       ...\n",
       "226385         246  5683  572.5485  692.1653\n",
       "226386         246  5684  572.6494  692.3114\n",
       "226387         246  5685  572.7741  692.3983\n",
       "226388         246  5686  572.8418  692.4915\n",
       "226389         246  5687  573.0393  692.8700\n",
       "\n",
       "[226390 rows x 4 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get all the filenames of the frames\n",
    "frames_files = [ filename for filename in os.listdir(os.path.join(\"GCT_final_real_data\",\"frames\")) if filename.endswith(\".dat\") ]\n",
    "\n",
    "# A regular expression to get the frame number from the filename\n",
    "regex = \"frame_(\\d+).*\"\n",
    "\n",
    "# Lists to store columns in the dataframe\n",
    "individuals = []\n",
    "times = [] \n",
    "xx = []\n",
    "yy = []\n",
    "\n",
    "# Read the trajectories\n",
    "for file_no, filename in tqdm.tqdm(enumerate(frames_files), desc=\"Reading files\", total=len(frames_files)):\n",
    "    \n",
    "    frame = int(re.search(regex, filename).group(1))  # Get frame from the filename\n",
    "    #print(filename, frame)\n",
    "    \n",
    "    with open(os.path.join(\"GCT_final_real_data\",\"frames\", filename)) as f:\n",
    "        line_count = 0\n",
    "        #print(f\"Reading file {file_no} / {len(frames_files)}\")\n",
    "        assert f.readline().startswith('#')  # First line should be a comment\n",
    "\n",
    "        while True:\n",
    "            line = f.readline()  # Now read remaining lines\n",
    "            if not line: \n",
    "                break\n",
    "            try: \n",
    "                l = line.strip().split(\" \")\n",
    "                individuals.append(int(l[0]))\n",
    "                xx.append(round(float(l[1]),4))\n",
    "                yy.append(round(float(l[2]),4))\n",
    "                times.append(frame)\n",
    "                line_count +=1\n",
    "            except ValueError as e:\n",
    "                print(f\"Exception reading file {filename}, line {line_count}: '{line}'\")\n",
    "                raise e\n",
    "\n",
    "traces = pd.DataFrame({'individual':individuals, 'time': times, 'x':xx, 'y':yy})\n",
    "traces = traces.sort_values(by=[\"time\"]).reset_index(drop=True)\n",
    "traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "traces.to_csv(os.path.join(\"GCT_final_real_data\", \"grand_central_terminal_traces.csv\"), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a couple of traces to check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAHSCAYAAAAubIVMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5DU9Z3v+9eb4YfIQnB04JIBdhSnzI8LgjsVhuJWihN3IuImzGbjogdWNuWRc288p0LYeBwDtcR7ZSXHFGGtk+seWLMHj2yCmqTlBII7ZaS21gI2Y2aEuOoFXIQZCEOCiIvIz/f9o79t2qFn+tvd3+7+dvfzUTXV3d/+dPebXx9e/fl8vp+vubsAAABQuGHlLgAAAKBaEKwAAAAiQrACAACICMEKAAAgIgQrAACAiBCsAAAAIjK83AVI0nXXXedNTU3lLgNACb3yyiu/cfeGctdRKPovoPYM1X/FIlg1NTWpq6ur3GUAKCEze7vcNUSB/guoPUP1X0wFAgAARIRgBQAAEBGCFQAAQEQIVgAAABEhWAEAAESEYAUAABCRrMHKzG4ys560n9NmttzM6s2s08z2B7fXBO3NzB43swNmttfMbin+LwMAAKD8sgYrd3/T3We6+0xJfyDpfUk/kdQh6UV3b5b0YvBYkm6X1Bz8LJP0RDEKBwAAiJtcNwi9VdJBd3/bzBZKmhcc3yRpp6QHJS2U9JS7u6TdZjbezCa5+7GIagZia1Vin36w54guuavOTHfPnqJH2qeXuywAyIg+K3q5Bqu7JP0guD8xFZbc/ZiZTQiON0o6kvaa3uAYwQpVbVVin57effjDx5fcP3xMRwUgThLdfVq+pecjx+izohF68bqZjZT0RUnPZmua4ZhneL9lZtZlZl0nTpwIWwYQW+mhKsxxVC76L1SqRHefbnho2xWhKt0P9hwZ9Dlkl8uI1e2Sfunux4PHx1NTfGY2SVJ/cLxX0pS0102WdHTgm7n7BkkbJKmlpeWK4AUAcUX/hUo0Y/UOnT53KWu7S85f6ULkst3C3frdNKAkbZW0NLi/VNLzacfvCc4ObJX0LuurAAAoj1WJfWrq2BYqVElSnWWaeEJYoUaszOxqSW2S/mPa4bWSnjGzeyUdlnRncHy7pAWSDih5BuFXIqsWAACEkuju09e39Fy5FieLu2dPyd4IgwoVrNz9fUnXDjj2WyXPEhzY1iXdH0l1QJVoW7dTnSvmlbsMADVi9ppOHX/vfM6vM7FwvVC5nhUIIA/7+8+UuwQANaBt3c6C+pt/XXtHhNXUJoIVAAAVbvHGXXr54Mm8Xz93Wr023zcnwopqF8EKiMjcafUFdWwAkI+mjm15v3bi2JHas7ItwmpAsAIisvm+OQV1cACQi0Kn/Q4x7VcUBCsAACpIodN+zRPGcDJNERGsAACoAAMvm5WrcaPqtPfh+RFWhEwIVgAAxFyho1RLWqeyjUKJEKwAAIixQtZScbZf6RGsAACIoXw3+ZSk9Ytmqn1WY8QVIQyCFQAAMVLIWqqr6kxvrFkQcUXIBcEKiMjsNZ3lLgFABWOTz+pAsAIiku+QPYDaVsiUn0SgihuCFQAAZZDo7tPyLT0FvQebfMYPwQoogSWtU8tdAoAYmbF6h06fu5T367kUTXwRrIASYP8YAClc26+6EayACCS6+8pdAoCYK2QtlUn6V6b9KsKwMI3MbLyZPWdmb5jZ62Y2x8zqzazTzPYHt9cEbc3MHjezA2a218xuKe4vASi/bzz7arlLABBjTR3b8g5VzRPGEKoqSNgRq7+WtMPdv2xmIyVdLembkl5097Vm1iGpQ9KDkm6X1Bz8zJb0RHALVK2Ll73cJQCIoXwXqA+TtI5NPitS1mBlZuMkfVbSn0uSu5+XdN7MFkqaFzTbJGmnksFqoaSn3N0l7Q5Guya5+7HIqwcAIIYS3X36+pYe5fqViw0+K1+YEasbJJ2Q9HdmdrOkVyR9TdLEVFhy92NmNiFo3yjpSNrre4NjHwlWZrZM0jJJmjqVM6YAVA76Lwwl37VULEyvDmHWWA2XdIukJ9x9lqQzSk77DcYyHLsitLv7BndvcfeWhoaGUMUCQBzQf2Ew+a6lWtI6lVBVJcKMWPVK6nX3PcHj55QMVsdTU3xmNklSf1r7KWmvnyzpaFQFAwAQN/lejma4SQceZWF6Nck6YuXuv5Z0xMxuCg7dKulfJG2VtDQ4tlTS88H9rZLuCc4ObJX0LuurAADVKNHdp6aObXmFqvWLZhKqqlDYswL/s6TNwRmBb0n6ipKh7Bkzu1fSYUl3Bm23S1og6YCk94O2AABUlUL2peJSNNUrVLBy9x5JLRmeujVDW5d0f4F1AQAQS23rdmp//5m8XssC9erHzutAgVYl9pW7BAAlcn3Htpy3UJCkudPqtfm+OZHXg/ghWAEFenr34SGfHzeqrkSVACiWfDf6lJJrqdjos3YQrIAi2/vw/HKXAKAAqxL7sn6BGgxrqWoPwQooABdfBqpLortPDzzbowuXy10JKhXBCihAtqmBJa3syg1Uinz3ogLSEayAPLWt25m1zSPt04tfCICCMEqFKBGsgDxlO916RJgLRgEom0LWTgGDIVgBeQgzWvXYnTOLXwiAnJVqhIqF67WJ79RAHrKNVl1VZ5xeDcRQatuEYoaqcaPqCFU1jGAFFMEbaxaUuwQAGeS7F1VYzRPGsMVKjWMqEMjRjNU7yl0CgJgZUWd67Ms3M1INghWQi1WJfTp97tKQbdhiAagtjeNH64HbbiJUQRLBCshJmDOI2GIBqA1c/w+ZEKyAkMLsss4WC0C8DTfpYj5XUQ4QppAN/w0AIT34o71Z27DFAhBvI4bn/98eoQphMGIFhHTu4tDnZy9pncoaCyDmzuaxz8LoEcP06Jdm8O8boYQKVmZ2SNJ7ki5JuujuLWZWL2mLpCZJhyT9qbu/Y2Ym6a8lLZD0vqQ/d/dfRl86UDqz13RmbcPaKiDeViX25dSeQIV85DJi9e/c/Tdpjzskvejua82sI3j8oKTbJTUHP7MlPRHcAhVpVWKfjr93fsg2zRPGlKgaAPkKe/kapvxQiELWWC2UtCm4v0lSe9rxpzxpt6TxZjapgM8ByipMZ9y5Yl7xCwFQVCOGSesXzSRUoSBhR6xc0j+YmUv67+6+QdJEdz8mSe5+zMwmBG0bJR1Je21vcOxY+hua2TJJyyRp6lT2/UHl4kzA2kP/VXmyndV7zdUj1P2Xny9RNahmYf9LmOvutyg5zXe/mX12iLaW4dgVJ7e6+wZ3b3H3loaGhpBlAKX1iZXbs7bhTMDaQ/9Veb619bUhn1/9hU+XqBJUu1DByt2PBrf9kn4i6TOSjqem+ILb/qB5r6QpaS+fLOloVAUDpZLo7tMHl4be8Gb9opksbAUqwKmzFwZ9btTwYfw7RmSyBiszG2NmY1P3JX1e0q8kbZW0NGi2VNLzwf2tku6xpFZJ76amDIFK8vUQF2ulMwYq37f/ZEa5S0AVCbPGaqKknyR3UdBwSX/v7jvM7BeSnjGzeyUdlnRn0H67klstHFByu4WvRF41UGSLN+66cv56gLnT6ktSC4DCZFtfxRckRClrsHL3tyTdnOH4byXdmuG4S7o/kuqAMnn54MmsbThzCKgMD/04+1UTgKhwPhOQh/WLWLAOVIp8dlsH8kWwAga4vmNb1jZMHQDV4ZqrR5S7BFQZghWQZsbqHVnXVi1pZd8ioFqwzQKiRrACAos37tLpc5eytuOagED1YPQZUSNYAUqGqjAL1jkTEKgs2c4IBKJGsELNS3T3hQpVEmcCApUm247rQNQIVqh533j21VDtOBMQqDxD7bgOFAPBCjUt0d2ni5ezLVfn0jVANWJqH8VAsEJNWx7isjVLWqcSqoAqxNQ+ioFghZo1Y/WOrG3GjarjLECgQrFwHeVAsEJNCru1wt6H55egGgDF8MCz2UekgagRrFBzwp4FyPoLoHIluvvElWxQDgQr1JwwZwGOG1XH+guggmXbZoFL2aBYCFaoOWHOAmQKEKhs2bZZ4FI2KBaCFWrK7DWdWduwXxVQ/TjTF8VCsELNWLxxl46/d37INhPHjqTDBaocF1JHMYUOVmZWZ2bdZvbT4PH1ZrbHzPab2RYzGxkcHxU8PhA831Sc0oHchFmwvmdlWwkqAVBObKGCYsplxOprkl5Pe/xtSd9192ZJ70i6Nzh+r6R33P1GSd8N2gFlFWbPqvGjWcwKVAP2r0I5hQpWZjZZ0h2S/jZ4bJI+J+m5oMkmSe3B/YXBYwXP3xq0B8qibd3OUHtWfeuLLGYFqsFjL7xZ7hJQw8KOWK2X9F8kpXYFuVbSKXe/GDzulZRamNIo6YgkBc+/G7T/CDNbZmZdZtZ14sSJPMsHstvffyZrm7nT6llbhdDov+Lt6Kmz5S4BNSxrsDKzP5LU7+6vpB/O0NRDPPe7A+4b3L3F3VsaGhpCFQsUC3tWIRf0X/H2sSzT+m3rdpamENSkMCNWcyV90cwOSfqhklOA6yWNN7PhQZvJko4G93slTZGk4PmPScq+ahgogsUbd2Vtw/YKQHXJtvhkf/8ZNXVsYy0WiiJrsHL3h9x9srs3SbpL0s/dfbGklyR9OWi2VNLzwf2twWMFz//c3bPvyAgUQbYzAZe0TmUKEKgyp94fenPQlOVbekJ9+QJyMTx7k0E9KOmHZvaIpG5JTwbHn5T0P83sgJIjVXcVViKQnzDfRjntGqg+Hx8/Wn0h11m9fPCkmjq2feSYSfruopl86UJecgpW7r5T0s7g/luSPpOhzQeS7oygNiBvqxL79PTuw0O2YZNAoDo9cNtNWr6lJ+/Xu/Th6wlXyBU7r6PqhAlVEqNVQLVqn9Wo9YtmFvwf3EM/3htJPagtBCtUnTChCkB1a5/VqLfW3qHmCWPyfo+zFy6zwB05I1ihqoTtBJkGBGpD54p5mjutPu/XryhgShG1iWCFqvLgj7IP3U8cO5JpQKCGbL5vjg6tvUMTx47M+bWXxSVykBuCFarKuYuXh3x+4tiRXGgZqFF7VrbltfaqkIXwqD0EK1SNMPvREKqA2pZae3Uox/VXs9d0FrEqVBOCFapCorsv1GagAJDSuWKe1i+aqfFZLoEjScffO8+UIEIhWKHiJbr7Qg3Vs64KwEDtsxrVs/rzoRa4MyWIMAhWqHgrnsne2Y3gbzqAIYS9EDtTgsiG/25Q8S6HuBLlY3dyoWUAQwszanX8vfMlqASVjGCFijZj9Y6sbeZOq+eyFACy2nzfHI0bVZe1HRduxlAIVqhYs9d06vS5S0O2aZ4wJvQQPwDsfXh+1jbZTpRBbSNYoSIluvtCDcl3rphX/GIAVJX1i7IvHbi+Y1sJKkElIlihIn09xNk5o1mxDiAP7bMaNXyYDdnGJbWt21mSelBZ+J8HFWVVYp+aOrYpxHp1PfqlGUWvB0B1+s6dN2dts7//jFYl9pWgGlSSrMHKzK4ys382s1fN7DUzezg4fr2Z7TGz/Wa2xcxGBsdHBY8PBM83FfeXgFqxKrFPT+8+HKrtktapLFgHkLf2WY2hri0Ytk9C7QgzYnVO0ufc/WZJMyXNN7NWSd+W9F13b5b0jqR7g/b3SnrH3W+U9N2gHVCwsB3Y3Gn1bAYKoGBhL4HFWYJIlzVYedK/BQ9HBD8u6XOSnguOb5LUHtxfGDxW8PytZjb0ZDWQxSdWbg/VbuLYkZwFCCAyYRayc5Yg0oVaY2VmdWbWI6lfUqekg5JOufvFoEmvpNS8S6OkI5IUPP+upGujLBq1ZfaaTn1wKfuqqnGj6rjIMoBIhZ0SZCE7UkIFK3e/5O4zJU2W9BlJn8zULLjNNDp1xf+KZrbMzLrMrOvEiRNh60WNWbxxV6htFcaNqgu1/wwQBfqv2hLmC9v+/jMlqASVIKezAt39lKSdkloljTez4cFTkyUdDe73SpoiScHzH5N0xTipu29w9xZ3b2loaMivelStRHefbujYFnqInVCFUqL/qj1LWqdmbcN1BCGFOyuwwczGB/dHS/pDSa9LeknSl4NmSyU9H9zfGjxW8PzP3T3M2fGApGSoWr6lR5dDtg/T4QFAIR5pn66r6oZeLsx1BCGFG7GaJOklM9sr6ReSOt39p5IelLTCzA4ouYbqyaD9k5KuDY6vkNQRfdmoZt949tXQbZe0TuUMQAAl8caaBVnbNLEje80bnq2Bu++VNCvD8beUXG818PgHku6MpDrUnER3ny5eDjfASagCUGrNE8ZkXU/Vtm4nl9OqYey8jlgJc6kaib2qAJRHmMDEQvbaRrBCbNz4ULhL1SxpncpeVQDKpnnCmHKXgBgjWCEWZq/p1MUQqWr9opmMVAEoK6b5MBSCFcou7F5VXP8PQKVgw9DaRbBCWc1e0xlqr6rmCWMYqQIQG3VZrtS2v/+MEt19JaoGcUKwQtmEHamSGHoHEC93z56Stc3ykCfjoLoQrFAWie6+0Luqh7kIKgCUEiPoGAzBCiWX2lk9DNZVAYgrrvqATAhWKLmwe1WxASiAOKN/QiYEK5RU27qdofeqotMCAFQaghVKZvHGXaF2JGZXdQBApSJYoSQWb9wVarH6xLEj2VUdQFWYO62+3CWgDAhWKLpczgDcs7KtyNUAQGnwJbE2EaxQdCueCbdYnW0VAACVjmCFolq8cZcuh1itzrYKACrN7DWd5S4BMUSwQtGEXVfFYnUAlWb2ms7QV45AbckarMxsipm9ZGavm9lrZva14Hi9mXWa2f7g9prguJnZ42Z2wMz2mtktxf5FIH7CrqtqnjCGdQgAKkrbup2EKgwqzIjVRUl/4e6flNQq6X4z+5SkDkkvunuzpBeDx5J0u6Tm4GeZpCcirxqxF3ZdFdcABFBJZqzeEWrbmGuuHlGCahBHWYOVux9z918G99+T9LqkRkkLJW0Kmm2S1B7cXyjpKU/aLWm8mU2KvHLE1qrEvlDrqlisDqCSLN64S6fPXQrVdvUXPl3kahBXw3NpbGZNkmZJ2iNporsfk5Lhy8wmBM0aJR1Je1lvcOxYocUi/trW7Qz1bY7F6gAqTdhtY+ZOq6d/q2GhF6+b2e9J+pGk5e5+eqimGY5dMX5hZsvMrMvMuk6cOBG2DMRY2FDVPGEMi9VR0ei/as+M1TtCtZs7rZ51ozUuVLAysxFKhqrN7v7j4PDx1BRfcNsfHO+VNCXt5ZMlHR34nu6+wd1b3L2loaEh3/oRE4nuvlChSmJdFSof/VdtSXT3hZoC5GQcSOHOCjRJT0p63d3XpT21VdLS4P5SSc+nHb8nODuwVdK7qSlDVK8Hf7Q3VDvWVQGoNMu3ZD8ZZ+60er40QlK4NVZzJf2ZpH1mlvrb9U1JayU9Y2b3Sjos6c7gue2SFkg6IOl9SV+JtGLE0rmLl7O2YV0VgEoTZgrwmqtHMFKFD2UNVu7+T8q8bkqSbs3Q3iXdX2BdqCCLN+7K2oZNQAFUmrZ1O0NNAXIGINKx8zoKlu1MGdYdAKg0YdeNjqwzRuLxEQQrFCTMtbJYdwCg0oRZVyVJ//XLNxe5ElSanPaxAtLd+NA2XcyyEejoEWR3AJVjVWKfnt59OFRb1o0iE/7XQ15mr+nMGqok6dEvzSh+MQAQgVxDFetGkQnBCnkJcwHSiWNH8m0OQMUIG6rY5BhDIVghZ6sS+0K127OyrciVAEA0wqwXlaSr6ox1oxgSwQo52xziWx0bgQKoFG3rdoYahZekN9YsKHI1qHQsXkfOsi2tWr9oJlOAAGIvlzVVEl8YEQ4jVojc917aX+4SAGBQie4+3fjN7TmHKr4wIgyCFSK3v/9M6PUKAFBKbet2avmWHl28HOK05gDbKiAXTAUiZ6NHDNPZC0NfG/D4e+fV1LGNb3kAyi7R3aeVP9mnM+ezX55mIM4ARK4IVsjZo1+aEXpX4lQ7whWAUkp09+nh//Wa3nn/Qt7v0TxhDGcAImdMBSJn7bMaNXHsyNDtv/Hsq0WsBgA+alVin5Zv6SFUoSwIVsjLnpVtGjeqLlTbi5ddM1bvKHJFAGpdPovSM1nSOpVQhbwxFYi87X14vhZv3KWXD57M2vb0uUtq6timudPqtfm+OSWoDkAtCdsXDYXL1CAKjFihIJvvm5PT3i4vHzzJ6BWAyCS6+3TDQ9sKClVLWqfq0No7CFWIRNZgZWbfN7N+M/tV2rF6M+s0s/3B7TXBcTOzx83sgJntNbNbilk84qF9VqMOrb0jdPvT5y6pbd3O4hUEoCaktk7IYeeEj5g7rZ5AhciFGbH6H5LmDzjWIelFd2+W9GLwWJJul9Qc/CyT9EQ0ZaIS5DJytb//jBLdfUWsBkC1SnT36fqObdrffyav16dGqFiWgGLIusbK3f/RzJoGHF4oaV5wf5OknZIeDI4/5e4uabeZjTezSe5+LKqCEV/tsxr1vZf2h+7slm/p0bNdh+ncAITWtm5nXoFq4tiRXBgeJZHvGquJqbAU3E4IjjdKOpLWrjc4hhrRuWKe5k6rD93+5YMn1dSxjdErAEMqZJRqSetUQhVKJuqzAi3DsYyz32a2TMnpQk2dOjXiMlBOqRGoXL5ZLt/So663T7LWARWB/qu08h2l4ixklEO+I1bHzWySJAW3/cHxXklT0tpNlnQ00xu4+wZ3b3H3loaGhjzLQJzlOnr19O7DjFyhItB/lc7sNZ15har1i2YSqlAW+QarrZKWBveXSno+7fg9wdmBrZLeZX1Vbdt835ycdmn/eshL5QCobqmpv+Pvnc/pdakz/biMFsol61Sgmf1AyYXq15lZr6TVktZKesbM7pV0WNKdQfPtkhZIOiDpfUlfKULNqDB7VrbpEyu364NL2c+JdomLNwM1LNHdpwd/tFfnLg59ofeBrqozvbFmQZGqAsKz5Al85dXS0uJdXV3lLgNFlus6Cc7iqW5m9oq7t5S7jkLRf0VnVWJfXpej4bp+KLWh+i92XkfJdK6Yl9NeV8ffO69PrNxexIoAxEG+1/gbN6pOh9beQahCrBCsUFK57tL+wSXX4o27ilgRgHJKdPdp+ZYeXcxx+/S50+q19+GBe1cD5UewQlnken1BANUnFapykRql4ow/xBXBCmXRPqtRS1rD7//T1LGtiNUAKKVEd59uWvWzvEIVo1SIO4IVyuaR9uk5jVwRroDKlxqlyvWsvyWtUwlVqAgEK5RVas1V2L2uCFdA5VqV2JfzKFXzhDE6tPYOrsqAikGwQizsWdmm5gljQrVlMTtQeRZv3JXzWX9LWqdyxh8qDsEKsdG5Yp6GZ7ra5AAsZgcqS9u6nTn/u12/aCajVKhIBCvEyoFHw23FwJQgUBnyuYAyV15AJSNYIXbC7nNFuALia1Vin5o6tuUUqlLrqQhVqGQEK8RS2HA1Y/WOIlcCIFf5rKdav2gm66lQFbJehBkolxHDpAtZzsg+fe5SaYoBEMqM1Tty/nfJ1B+qCSNWiK3H7gy3xxVTgkD5pab+cglVTP2hGhGsEFvtsxpDbyDKxZqB8mlbtzPnqb/mCWOY+kNVIlgh1sJetJmLNQOll88CdSl5AWVCFaoVa6xQEeZOq8+6Dw77WwGlkeju04pnenTZc38t66lQ7YoyYmVm883sTTM7YGYdxfgM1JawV7Jn1AoonkR3n2785nYt35J7qJo7rZ71VKgJkQcrM6uT9D1Jt0v6lKS7zexTUX8Oak+Y9VYvHzypRHdfCaoBakvq4skX8ximWr9oZugvR0ClK8aI1WckHXD3t9z9vKQfSlpYhM9BjWmf1aglrVOztnvshTdLUA1QW7619bWcXzNx7EhGqVBzihGsGiUdSXvcGxwDChbm2mFHT50tQSVAbTl19kJO7ZsnjNGelW1FqgaIr2IEq0yX0b1i7NjMlplZl5l1nThxoghloFplG7X6+PjRJaoEtYr+a3AjhrGLOmpbMYJVr6QpaY8nSzo6sJG7b3D3FndvaWhoKEIZqFaPtE/XVXWZ8nvSA7fdVMJqUItqsf+65uoRWdusXzRT+/+KqT/UtmIEq19Iajaz681spKS7JG0twueghr2xZoEmjh15xXFO5QaKY/UXPq0Rg3yhWdI6lbVUQCDyfazc/aKZ/SdJL0iqk/R9d8991SOQBes3gNJJhabHXnhTR0+d1cfHj9YDt91EmAIGKMoGoe6+XRLXGAGAKtI+q5EgBWTBJW0AAAAiQrACAACICMEKAAAgIgQrAACAiBCsAAAAIkKwAgAAiAjBCgAAICLmfsVl/EpfhNkJSW9Luk7Sb8pczkBxrEmKZ11xrEmKZ11xrEkqbV2/7+4Vfz2YmPdfUjzrimNNUjzrimNNUjzrikX/FYtglWJmXe7eUu460sWxJimedcWxJimedcWxJim+dVWCuP7exbGuONYkxbOuONYkxbOuuNTEVCAAAEBECFYAAAARiVuw2lDuAjKIY01SPOuKY01SPOuKY01SfOuqBHH9vYtjXXGsSYpnXXGsSYpnXbGoKVZrrAAAACpZ3EasAAAAKhbBCgAAICIEKwAAgIgQrAAAACJCsAIAAIgIwQoAACAiBCsAAICIEKwAAAAiQrACAACICMEKAAAgIgQrAACAiBCsAAAAIkKwAgAAiAjBCgAAICIEKwAAgIgQrAAAACJCsAIAAIgIwQoAACAiBCsAAICIEKwAAAAiQrACAACICMEKAAAgIgQrAACAiBCsAAAAIkKwAgAAiAjBCgAAICIEKwAAgIgQrAAAACJCsAIAAIgIwQoAACAiBCsAAICIEKwAAAAiQrACAACIyPByFyBJ1113nTc1NZW7DAAl9Morr/zG3RvKXUeh6L+A2jNU/xWLYNXU1KSurq5ylwGghMzs7XLXEAX6L6D2DNV/MRUIAAAQEYIVAABARAhWAAAAESFYAQAARIRgBQAAEBGCFQAAQESyBiszu8nMetJ+TpvZcjOrN7NOM9sf3F4TtDcze9zMDpjZXjO7pfi/DAAAgPLLuo+Vu78paaYkmVmdpD5JP5HUIelFd19rZh3B4wcl3S6pOfiZLemJ4BYAAFSIxRt36eWDJz98PHdavTbfN6eMFVWGXKcCb5V00N3flrRQ0qbg+CZJ7cH9hZKe8qTdksab2aRIqgUAAEU3e03nR0KVJLNH3HMAAB8oSURBVL188KQWb9xVpooqR647r98l6QfB/YnufkyS3P2YmU0IjjdKOpL2mt7g2LFCCk1JdPdpxZYeXU471jxhjDpXzIvi7QEAqGmrEvt0/L3zGZ8bGLZwpdAjVmY2UtIXJT2brWmGY57h/ZaZWZeZdZ04cSJUDYnuPi0fEKokaX//Gc1e0xnqPQCgUPn0X0CleHr34XKXUNFymQq8XdIv3f148Ph4aoovuO0PjvdKmpL2usmSjg58M3ff4O4t7t7S0BDuOqyPvfDmoM8df++82tbtDPU+AFCIfPovALUhl2B1t343DShJWyUtDe4vlfR82vF7grMDWyW9m5oyLNTRU2eHfJ6RKwAA8scaqsKFClZmdrWkNkk/Tju8VlKbme0PnlsbHN8u6S1JByRtlPTVqIr9+PjRWdscf+884QoAgDxkW0O1pHVqiSqpXKGClbu/7+7Xuvu7acd+6+63untzcHsyOO7ufr+7T3P36e7eFVWxD9x2U6h2x987r6aObUp090X10QAAVLUwo1WPtE8vQSWVraJ2Xm+f1ai50+pDt1++pUerEvuKWBEAAJUv0/YKA00cO7JE1VS2igpWkrT5vjk5/eE+vfswI1cAAAxi9prOQbdXSLdnZVsJqql8FRespOQf7lV1mXZ1yGz5lh4W5AEAkEGYUIXwKjJYSdIbaxZo3Ki60O1fPnhSM1bvKGJFAABUlrCDDixaD69ig5Uk7X14fk5rrk6fu6Smjm2MXgEAalqiu0+f/ssdoXZSb54whkXrOajoYCUl11ytXzQzp9e8fPCkruesQQBADVqV2KflW3p05vylrG0njh3JJeNyVPHBSkqeLXho7R0aHn7ZlVysvQIA1JZEd19Ol6xhwXruqiJYpRx4NLdwJSVHr7gUDgCgFizf0hO6Leuq8lNVwUpKhqtc99rY33+GcAUAqGqfWLk9dNtxo+pYV5WnqgtWUnLoMpdF7VIyXDEtCACoRk0d2/TBJQ/V9qo6096H5xe5oupVlcFKyn9ROyNXAIBq0tSxLXTbudPq9caaBUWspvpVbbCSfreoPZfRq/39Z3TjQ+H/EgIAEFe5XNbt0No7tPm+OUWspjZUdbBK2XzfnJzOGrzoyYTPdQYBAJUs7BmAuS6fweBqIlil5Lqw/endhwlXAICKNHtNZ6h2c6fVM1IVoZoKVlJyYXvzhDGh2+ey3wcAAHER5hqATP9Fr+aClSR1rpiX03UGucYgAKCShLmyyKG1d5SgktoTKliZ2Xgze87M3jCz181sjpnVm1mnme0Pbq8J2pqZPW5mB8xsr5ndUtxfQn72Pjw/9MjV6XOXCFcAgIqRbSPQXPd7RHhhR6z+WtIOd/+EpJslvS6pQ9KL7t4s6cXgsSTdLqk5+Fkm6YlIK45Q54p5obdkOH3uUk6bqwEAUA5h1gZzqZriyRqszGycpM9KelKS3P28u5+StFDSpqDZJkntwf2Fkp7ypN2SxpvZpMgrj0hqS4YwJwx+cMnZigEAEGvZ1gYPH5bjtd+QkzAjVjdIOiHp78ys28z+1szGSJro7sckKbidELRvlHQk7fW9wbFY+27IkauLzporAEDl+s6dN5e7hKoWJlgNl3SLpCfcfZakM/rdtF8mmaLwFfvom9kyM+sys64TJ06EKraY2mc1hr7g5Olzl9ihHahhceu/gJRs04DDLfn/HYonTLDqldTr7nuCx88pGbSOp6b4gtv+tPZT0l4/WdLRgW/q7hvcvcXdWxoaGvKtP1KPtE8PvUka1xYEalcc+y9Ayj4NeOBRzgQstqzByt1/LemImd0UHLpV0r9I2ippaXBsqaTng/tbJd0TnB3YKund1JRhJdh835zQ4YprCwIA4oIv+/EQ9qzA/yxps5ntlTRT0l9JWiupzcz2S2oLHkvSdklvSTogaaOkr0ZacQlsvm9O6GnB/f1nOFsQAFB2Lx88OeTzbLFQGsPDNHL3HkktGZ66NUNbl3R/gXWV3SPt0yWF23n9g0uupo5tWr9oJnPXAICSCzNaxRYLpVGTO6+H9Uj79Jwuf7N8Sw9DsQCAkss2WsUOC6VDsMqic8W8nIZPWXcFAIibdX8abkshFI5gFcKelW05XVtwf/8ZwhUAoCTCzJSwTKV0CFYh7X14fk4jV/v7z2j2ms4iVgQAQPZpwLAnYyEaBKsc7FnZltNf0OPvnWeXdgBA0YQZrUqdjIXSIFjl6JH26aEv3Cwld2m/vmObEt19RawKAFBrViX2ZR2tCrsvI6JDsMpD6sLNYacGXckzBsNccRwAgGxmr+kMtR3Q5vvmlKAapCNYFWDPyjZdVRf+HNandx9m5AoAUJC2dTt1/L3z5S4DgyBYFWhYjpuDLN/SQ7gCAOQl0d2n/f1nQrVl0Xp5hNp5HYN7/8LlnF+zfEuPlm/p0TCT/v3sqSwsBACEsnxLT6h2w4cZ/7eUCcGqjC57cnowfZ581PBh+vafzGDPEQDAR+SyP+J37ry5eIVgSASrmDl38fKHI1pXjximUSPqdOr9C/r4+NF64LabCFwAUIMWb9wVegqQ69aWF2usYuz9C5f1zvsX5JL6Tp3V8i09+vRf7mCNFgDUkER3X9ZtFSSpecIYHVp7B6GqzBixKtDcafWh/sJH5cz5Sx+OaKVb0spaLQCIm0R3n1b+ZJ/OnL9U1M+5qs7UuWJeUT8D4RCsCrT5vjmavaaz7Ke+DlyrNWZkndb88XS+uQCoeaUKN+X0xpoF5S4BAYJVBPasbFNTx7Zyl/ERqZGtrrdPMpIFoCokuvv0F8/06JKXu5J4YVuFeAm1xsrMDpnZPjPrMbOu4Fi9mXWa2f7g9prguJnZ42Z2wMz2mtktxfwFxEUul7kppad3H1ZTxzY1dWxj53cAFSvR3aflWwhVA11Vx7YKcZPLiNW/c/ffpD3ukPSiu681s47g8YOSbpfUHPzMlvREcFvVUlNu33j2VV28HM9/+enTheNG1Wnvw/PLXBEAhPPYC2+Wu4RYYgowfgo5K3ChpE3B/U2S2tOOP+VJuyWNN7NJBXxOxWif1agDf7WgIoZlT5+7pBmrd5S7DAAI5eips+UuAQglbLBySf9gZq+Y2bLg2ER3PyZJwe2E4HijpCNpr+0NjtWMR9qn69DaO7Skdapyu+BNZuNG1UXwLlc6fa56F3ICqC4fHz+63CXETlyXoNS6sFOBc939qJlNkNRpZm8M0TZTlrhibiwIaMskaerU+I/w5OOR9ukfmfuesXpHXmHm9LlLOrT2jo8cS3T36VtbX9OpsxcKrhNAbmqh/4qbB267KfTlXOLIJC0uYFucRHefHnvhTR09dZYNo2MuVLBy96PBbb+Z/UTSZyQdN7NJ7n4smOrrD5r3SpqS9vLJko5meM8NkjZIUktLSzwXJUVs78PztXjjrkj2vWqf1XjFP6pViX0f2XIBQHHUYv9Vbqn+rtThau60em2+b05JPzOTTH0+4inrVKCZjTGzsan7kj4v6VeStkpaGjRbKun54P5WSfcEZwe2Sno3NWWI5L5XxVqD9Uj7dK1fNFOjR7ChPoDq0z6r8cNlFqXy8sGTaurYpsUbd5XsM1HZzH3oL1tmdoOknwQPh0v6e3dfY2bXSnpG0lRJhyXd6e4nzcwk/TdJ8yW9L+kr7t411Ge0tLR4V9eQTarSqsQ+/f2ew8p2EmEh35iy7a81cIoRKBUze8XdW8pdR6Fqtf+Ki0R3nx54tkcXLpevBq58UXuG6r+yBqtSoGMafBqv0GFoghXiimCFYkp09+mhH+/V2RInrmuuHqHVX/g003ZVbqj+i53XY2LgQvdSqIRtIQAgH+lrkpIha5/OXij+mdDvvH9By7f06Nmuw7FYm4XSYzFOFUt09w35PEPXAGpB+6xGPfql6Wos4ZYNLx88qZtW/SxrP4zqw4hVFavkU5MBIEoDz6pLbV/Qd+qs6sx0qQjLYs5dvKzlW3o+0hczVVj9CFYAgJoz1PYFxdy65p33L2jFMz0f1oDqw1QgAABp0q+eUQyXPTmjcOM3tzNVWIUYsQIAIIOhTiqKYrPni5ddy7f06Hsv7VfninkFvRfig2BVw2785nZ9586bGY4GgByln/FXaMja33/miq1xxoys05o/nk7/XIEIVlVsZJ3p/KXBF2Smvi1JzPUDQL423zcn8nVZZ85f+nDh+6jhw/TtP5mRsZ+evaZTx987H+o9h5n072ezmWmxsUFoFUt09+V1ZiBnraAU2CAU1Sb9TMM4a54whqnHArHzeg278ZvbdTHbNXOGsH7RTEnSt7a+plNnL0gieCEaBCtUu0R3n1Zs6VEZr7YTSlwuNF1Jhuq/OCuwyn3nzpsLen1qKDoVqqTk6cIPPPfqh2ezJLr7NHftz3V9xzbNXftzznIBACWXWLy19g7NnVZf7lKG9PLBk5yhGCFGrGpA27qd2t9/JvL3bRw/Wg/cdpMeePZVXUgbFRsxzPQYi+KRBSNWqEX5LtEoh6HWdtU6RqxqXOeKeZo4dmTk73v01Fl9a+trHwlVknThsuvrz/Tw7QcABmif1aj1i2Zq/OgR5S4lq9TO8Tc8tI3+PAeMWNWQKPZdSdc4fvSQizRN0uJWzkBBZoxYAVcq5q7vUVpS4307i9fxoSinBdcvmslZh8gbwQrILNHdd8USi7AGLkQv9gL6Wg1YBCt8RBThavSIYXr9/7lds/7vf9A771/I/oIMOBOlthGsgMGltm44euqsPh6sZy3ky2ipRsJq5YtzJMHKzOokdUnqc/c/MrPrJf1QUr2kX0r6M3c/b2ajJD0l6Q8k/VbSInc/NNR70zGV3mD/aMNOF65fNFPtsxqT36yee1UXhtiINIzGCDoOVBaCFVAepVhAP3yYVfWVPaIKVisktUgaFwSrZyT92N1/aGZ/I+lVd3/CzL4qaYa7/59mdpekP3b3RUO9Nx1TvGT7Rzdw6DfR3feRfa6iVCvffmoRwQoon+R0Y48ulGiTrWqboSg4WJnZZEmbJK2RtELSFySdkPS/uftFM5sj6VvufpuZvRDc32VmwyX9WlKDD/FBdEzxtCqxT5t3H1bqDy7btasGti8mTgOufAQrIJ6KtUXPuFF12vvw/MjftxyiCFbPSXpU0lhJ35D055J2u/uNwfNTJP3M3f93M/uVpPnu3hs8d1DSbHf/zYD3XCZpmSRNnTr1D95+++08f3mIk1Jf0uHqEcP0V18iYFWiSg5W9F+oBcVel1XJo1gFBSsz+yNJC9z9q2Y2T8lg9RVJuwYEq+3uPt3MXpN024Bg9Rl3/+1gn8E3vupUzCnCMLjgaLxVcrBKR/+FWlCsqcNKXYs1VP81PMTr50r6opktkHSVpHGS1ksab2bD3f2ipMmSjgbteyVNkdQbTAV+TFJ0myehYrTPavzwH0s5dhu+7NLTuw9/5BsXe2sBQO5S/Xmiu08P/6/X8j4bfKCLl13Lt/So6+2TVdMv57TdQmrEKli8/qykH6UtXt/r7v+vmd0vaXra4vUvufufDvW+fOOrHcWau88Va7TKjxEroDpEuVdWpXz5jWwfqwHB6gb9bruFbklL3P2cmV0l6X9KmqXkSNVd7v7WUO9Lx1SbEt19eujHe3W2VKelZMHUYWkRrIDqEuV0YfOEMepcMa/wNyoSNghFRSl34KrUOf9KQ7ACqlNUfXicF7cTrFDRop7Tz1e27SaQG4IVUN2iCFhx/aJLsEJVKfWWDkNhx/j8EayA2lLIGtu4jV4RrFAz4nRleEa4hkawAmpP2MumZRKndVcEK9SUgddBbLp2dN7/kKNCyLoSwQqoTVFMEZa7TyVYoebFafpwMMMsufdWrUwvEqwASIXNNKxfNLMsfSXBChhCIUPTxVbNF6EmWAFIKaQfLke4IlgBWZT78ju5qobARbACkK7QNbKlXOBOsALyVO49tXK1pAJ2LE4hWAHI5MZvbtfFy4Vlk6tHDNNffal4V9cgWAERqrTRLUmqM9Pds6fEKnQRrABkkuju0wPPvaoLl6LJJ8UYySJYAUVUaaNaUjyCFsEKwGCKecKRSfpugeuyCFZAGVTayFapT18mWAHIxYzVO3T63KXI3q+QRe8EKyCm4rShaa4KXUBPsAKQq0J2b8/k0No78nodwQqoEJWw31Ym+axhIFgByEeiu08PPNujKFZfEKyAGlUJ04q5hiuCFYBCFboPYb6XyRmq/xqedzUASqZ9VuOgU25xWTwf101WAVSv9C9z+YSsKKcVUwhWQIUbLHRFOVwOAHGXClnl7vuyBiszu0rSP0oaFbR/zt1Xm9n1kn4oqV7SLyX9mbufN7NRkp6S9AeSfitpkbsfKlL9AAaRKXDF+fI9ABCFgX1fU8e2IdvPWL1Dex+eH9nnhxmxOifpc+7+b2Y2QtI/mdnPJK2Q9F13/6GZ/Y2keyU9Edy+4+43mtldkr4taVFkFQPI22BroKIIXHOn1Rf0egAohuYJY4ac8otyCwdJGpatgSf9W/BwRPDjkj4n6bng+CZJ7cH9hcFjBc/famYWWcUAIrf5vjk6tPaOD3/WL5qpxvGjJSU308umlNfoAoBchFmcnm1UKxeh1liZWZ2kVyTdKOl7kg5KOuXuF4MmvZJS426Nko5IkrtfNLN3JV0r6TcD3nOZpGWSNHXq1MJ+FQAiNdRiedB/ARhc1hErSXL3S+4+U9JkSZ+R9MlMzYLbTF9wr9jTwd03uHuLu7c0NDSErRcAyo7+C6gsS1qzfwFKdPdF8lmhglWKu5+StFNSq6TxZpYa8Zos6Whwv1fSFEkKnv+YJFbLAgCAsghzXdTHXngzks/KGqzMrMHMxgf3R0v6Q0mvS3pJ0peDZkslPR/c3xo8VvD8zz0Ou5ACAICalW3U6mhEV7wIM2I1SdJLZrZX0i8kdbr7TyU9KGmFmR1Qcg3Vk0H7JyVdGxxfIakjkkoBAADylG3U6uPBCTuFyrp43d33SpqV4fhbSq63Gnj8A0l3RlIdAABARA6tvUPXd2y7YuH36BF1euC2myL5jJzWWAEAAFSyf03bUsYkNY4frUe/ND2yM6G5pA0AAKgpxdxShhErAACAiBCsAAAAIkKwAgAAiAjBCgAAICIEKwAAgIgQrAAAACJCsAIAAIgIwQoAACAiBCsAAICIEKwAAAAiQrACAACICMEKAAAgIgQrAACAiBCsAAAAIpI1WJnZFDN7ycxeN7PXzOxrwfF6M+s0s/3B7TXBcTOzx83sgJntNbNbiv2LAAAAiIMwI1YXJf2Fu39SUquk+83sU5I6JL3o7s2SXgweS9LtkpqDn2WSnoi8agAAgBjKGqzc/Zi7/zK4/56k1yU1SlooaVPQbJOk9uD+QklPedJuSePNbFLklQMAAMRMTmuszKxJ0ixJeyRNdPdjUjJ8SZoQNGuUdCTtZb3BsYHvtczMusys68SJE7lXDgBlQv8FYDChg5WZ/Z6kH0la7u6nh2qa4ZhfccB9g7u3uHtLQ0ND2DIAoOzovwAMJlSwMrMRSoaqze7+4+Dw8dQUX3DbHxzvlTQl7eWTJR2NplwAAID4CnNWoEl6UtLr7r4u7amtkpYG95dKej7t+D3B2YGtkt5NTRkCAABUs+Eh2syV9GeS9plZT3Dsm5LWSnrGzO6VdFjSncFz2yUtkHRA0vuSvhJpxQAAADGVNVi5+z8p87opSbo1Q3uXdH+BdQEAAFQcdl4HAACICMEKAAAgIgQrAACAiBCsAAAAIkKwAgAAiAjBCgAAICIEKwAAgIgQrAAAACJCsAIAAIgIwQoAACAiBCsAAICIEKwAAAAiQrACAACICMEKAAAgIgQrAACAiGQNVmb2fTPrN7NfpR2rN7NOM9sf3F4THDcze9zMDpjZXjO7pZjFAwAAxEmYEav/IWn+gGMdkl5092ZJLwaPJel2Sc3BzzJJT0RTJgAAQPxlDVbu/o+STg44vFDSpuD+Jkntacef8qTdksab2aSoigUAAIizfNdYTXT3Y5IU3E4IjjdKOpLWrjc4BgAAUPWiXrxuGY55xoZmy8ysy8y6Tpw4EXEZAFA89F8ABpNvsDqemuILbvuD472SpqS1myzpaKY3cPcN7t7i7i0NDQ15lgEApUf/BWAw+QarrZKWBveXSno+7fg9wdmBrZLeTU0ZAgAAVLvh2RqY2Q8kzZN0nZn1Slotaa2kZ8zsXkmHJd0ZNN8uaYGkA5Lel/SVItQMAAAQS1mDlbvfPchTt2Zo65LuL7QoAACASsTO6wAAABEhWAEAAESEYAUAABARghUAAEBECFYAAAARIVgBAABEhGAFAAAQEYIVAABARAhWAAAAESFYAQAARIRgBQAAEBGCFQAAQEQIVgAAABEhWAEAAESEYAUAABARghUAAEBEihKszGy+mb1pZgfMrKMYnwEAABA3kQcrM6uT9D1Jt0v6lKS7zexTUX8OAABA3BRjxOozkg64+1vufl7SDyUtLMLnAAAAxEoxglWjpCNpj3uDYwAAAFWtGMHKMhzzKxqZLTOzLjPrOnHiRBHKAIDioP8CMJhiBKteSVPSHk+WdHRgI3ff4O4t7t7S0NBQhDIAoDjovwAMphjB6heSms3sejMbKekuSVuL8DkAAACxMjzqN3T3i2b2nyS9IKlO0vfd/bWoPwcAACBuIg9WkuTu2yVtL8Z7AwAAxBU7rwMAAESEYAUAABARghUAAEBECFYAAAARIVgBAABEhGAFAAAQEYIVAABARMz9isv4lb4IsxOS3pZ0naTflLmcgeJYkxTPuuJYkxTPuuJYk1Taun7f3Sv+ejAx77+keNYVx5qkeNYVx5qkeNYVi/4rFsEqxcy63L2l3HWki2NNUjzrimNNUjzrimNNUnzrqgRx/b2LY11xrEmKZ11xrEmKZ11xqYmpQAAAgIgQrAAAACISt2C1odwFZBDHmqR41hXHmqR41hXHmqT41lUJ4vp7F8e64liTFM+64liTFM+6YlFTrNZYAQAAVLK4jVgBAABUrJIGKzP7vpn1m9mv0o7Vm1mnme0Pbq8JjpuZPW5mB8xsr5ndUqSappjZS2b2upm9ZmZfK3ddZnaVmf2zmb0a1PRwcPx6M9sT1LTFzEYGx0cFjw8EzzdFXdOA+urMrNvMfhqHuszskJntM7MeM+sKjpX171XwWePN7DkzeyP4+zWnzH+vbgp+j1I/p81seRx+ryoB/VdOdcW2D4tb/xV8Vuz6sLj1X8HnVEYf5u4l+5H0WUm3SPpV2rH/KqkjuN8h6dvB/QWSfibJJLVK2lOkmiZJuiW4P1bS/yfpU+WsK3jv3wvuj5C0J/isZyTdFRz/G0n/V3D/q5L+Jrh/l6QtRf5zXCHp7yX9NHhc1rokHZJ03YBjZf17FXzWJkn/Ibg/UtL4ONQVfF6dpF9L+v241BT3H/qvnOqKbR8Wt/4reP/Y9WFx7r+Cz4xtH1b0D8jwm9E0oGN6U9Kk4P4kSW8G9/+7pLsztStyfc9LaotLXZKulvRLSbOV3PhseHB8jqQXgvsvSJoT3B8etLMi1TNZ0ouSPifpp8Ff2LLWNUinVNY/P0njJP3rwF9vuetKe//PS3o5TjVVwg/9V141xaYPi2P/Fbx/rPqwuPdfwWfEtg+Lwxqrie5+TJKC2wnB8UZJR9La9QbHiiYY6p2l5LerstYVDFf3SOqX1CnpoKRT7n4xw+d+WFPw/LuSro26psB6Sf9F0uXg8bUxqMsl/YOZvWJmy4Jj5f57dYOkE5L+Lph2+FszGxODulLukvSD4H5caqpEsfm9i1P/FdQTxz4sjv2XFL8+LO79lxTjPiwOwWowluGYF+3DzH5P0o8kLXf300M1zXAs8rrc/ZK7z1TyG9ZnJH1yiM8tSU1m9keS+t39lfTD5a5L0lx3v0XS7ZLuN7PPDtG2VDUNV3La6Al3nyXpjJJD1OWuS8Eaki9KejZb0wzHivZvsMrUdP8lxa8Pi3H/JcWvD4tt/yXFvw+LQ7A6bmaTJCm47Q+O90qaktZusqSjxSjAzEYo2Sltdvcfx6UuSXL3U5J2Kjk/PN7Mhmf43A9rCp7/mKSTRShnrqQvmtkhST9Ucjh9fbnrcvejwW2/pJ8o2YmX+8+vV1Kvu+8JHj+nZEdV7rqkZOf9S3c/HjyOQ02Vquy/d3Huv6RY9WGx7L+kWPZhce6/pJj3YXEIVlslLQ3uL1VyjUDq+D3Bqv5WSe+mhvqiZGYm6UlJr7v7ujjUZWYNZjY+uD9a0h9Kel3SS5K+PEhNqVq/LOnnHkwoR8ndH3L3ye7epOQw7M/dfXE56zKzMWY2NnVfyXn3X6nMf6/c/deSjpjZTcGhWyX9S7nrCtyt3w2hpz673DVVKvqvzHXFrg+LY/8lxbMPi3n/JcW9Dyv2Iq70n+A34pikC0omyXuVnLN+UdL+4LY+aGuSvqfkvPw+SS1Fqun/UHJocK+knuBnQTnrkjRDUndQ068k/WVw/AZJ/yzpgJJDoKOC41cFjw8Ez99Qgj/LefrdWTVlqyv47FeDn9ckrQyOl/XvVfBZMyV1BX+OCUnXlLsuJRcS/1bSx9KOlf33qhJ+6L9yqivWfVhc+q+0z49dHxbH/iv4rNj3Yey8DgAAEJE4TAUCAABUBYIVAABARAhWAAAAESFYAQAARIRgBQAAEBGCFQAAQEQIVgAAABEhWAEAAETk/wdmau1DJAheIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(2,2, sharex=True, sharey=True, figsize=(10,8))\n",
    "\n",
    "for i, individual in enumerate(random.choices(traces.individual.unique(), k=4)):\n",
    "    ax = axes.flatten()[i]\n",
    "    trace = traces.loc[traces.individual==individual, [\"x\", \"y\"]]\n",
    "    ax.scatter(x=trace.x, y=trace.y)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}