From ba60fa949be99e9bdba3e104aead402509de0db1 Mon Sep 17 00:00:00 2001 From: sfmig <33267254+sfmig@users.noreply.github.com> Date: Thu, 28 Nov 2024 18:59:11 +0000 Subject: [PATCH] Add notebook to explore trajectories of escape clips (#257) * add movement as dependency * Rename existing movement notebook * Add notebooks dependencies * Rename old movement notebook * Add notebook to generate plot of escape clips * Update path * Add plot saving * Add as a script (to run in interactive job in the cluster) * Small edits to script * Update notebook * Add notebook for visualising trajectories with movement * Add last version of .py vscode notebook (to keep it in history) * Delete duplicate * Remove movement as dependency and `[notebooks]` as an extra --- MANIFEST.in | 1 + notebooks/notebook_movement_escapes.ipynb | 401 ++++++++++++++++++ ...y => notebook_movement_trajectories_gt.py} | 6 +- scripts/escape_trajectory_plots.py | 138 ++++++ 4 files changed, 543 insertions(+), 3 deletions(-) create mode 100644 notebooks/notebook_movement_escapes.ipynb rename notebooks/{notebook_trajectories_movement.py => notebook_movement_trajectories_gt.py} (97%) create mode 100644 scripts/escape_trajectory_plots.py diff --git a/MANIFEST.in b/MANIFEST.in index 5c94df88..e2de5f30 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -6,6 +6,7 @@ recursive-include guides *.md recursive-include crabs/tracker *.md recursive-include bash_scripts *.sh recursive-include notebooks *.py +recursive-include notebooks *.ipynb recursive-include scripts *.py recursive-include crabs *.yaml recursive-include guides *.png diff --git a/notebooks/notebook_movement_escapes.ipynb b/notebooks/notebook_movement_escapes.ipynb new file mode 100644 index 00000000..1bf1b20e --- /dev/null +++ b/notebooks/notebook_movement_escapes.ipynb @@ -0,0 +1,401 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inspect crab escape trajectories using `movement`.\n", + "\n", + "### Requirements\n", + "1. Create and activate a conda environment with `movement` by following the \n", + " instructions at https://movement.neuroinformatics.dev/getting_started/installation.html\n", + "\n", + "2. Install some additional dependencies on the conda environment by running:\n", + " ```\n", + " pip install ipykernel ipympl\n", + " ```\n", + "\n", + "3. Mount the zoo directory in ceph following the guide at\n", + " https://howto.neuroinformatics.dev/programming/Mount-ceph-ubuntu-temp.html\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import required packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7abaed3-a297-40fe-8d2d-a52dea9f6254", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from movement.io import load_bboxes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Uncomment and run the following line to enable interactive plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93cfc22f-f64d-4b7f-8f76-fa15e37f5931", + "metadata": {}, + "outputs": [], + "source": [ + "# %matplotlib widget" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Set input and output data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79c9824a-2960-40f8-8937-ca3255942e59", + "metadata": {}, + "outputs": [], + "source": [ + "# Ensure the input data points to the directory containing the \n", + "# csv files in ceph\n", + "input_data = Path(\n", + " \"/ceph/zoo/users/sminano/escape_clips_tracking_output_slurm_5699097\"\n", + " \n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Set output directory for figures" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2944aabe-b53e-4faf-903a-4c9fb676da38", + "metadata": {}, + "outputs": [], + "source": [ + "output_figures_dir = \"path/to/directory/where/figures/will/be/saved\" \n", + "# replace with actual path\n", + "\n", + "# Create output directory if it doesnt exist\n", + "if not output_figures_dir.exists():\n", + " output_figures_dir.mkdir(parents=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "List all .csv files in the input directory" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "d147111e-69cd-47c0-be94-aa7b1e64d0ed", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "234\n" + ] + } + ], + "source": [ + "list_csv_files = [\n", + " x\n", + " for x in input_data.iterdir()\n", + " if x.is_file() and x.name.endswith(\"_tracks.csv\")\n", + "]\n", + "list_csv_files.sort()\n", + "print(len(list_csv_files))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Read first file as movement dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c7fdb4b3-e130-4f4d-b499-b3836f1452ed", + "metadata": {}, + "outputs": [], + "source": [ + "csv_file = list_csv_files[0]\n", + "\n", + "ds = load_bboxes.from_via_tracks_file(\n", + " csv_file, fps=None, use_frame_numbers_from_file=False\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Print summary metrics of dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "697a1fcb-1f14-4f95-83ec-40696da9a710", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "04.09.2023-01-Right-Spontaneous1_tracks.csv\n", + "Number of frames: 187\n", + "Number of individuals: 105\n", + " Size: 789kB\n", + "Dimensions: (time: 187, individuals: 105, space: 2)\n", + "Coordinates:\n", + " * time (time) int64 1kB 0 1 2 3 4 5 6 ... 180 181 182 183 184 185 186\n", + " * individuals (individuals) " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(1, 1)\n", + "for ind_idx in range(ds.sizes[\"individuals\"]):\n", + " # plot trajectories\n", + " ax.scatter(\n", + " x=ds.position[:, ind_idx, 0], # nframes, nindividuals, x\n", + " y=ds.position[:, ind_idx, 1],\n", + " s=1,\n", + " color=list_colors[ind_idx % len(list_colors)],\n", + " )\n", + " # add ID at first frame with non-nan x-coord\n", + " if flag_plot_id:\n", + " start_frame = ds.time[~ds.position.isnull()[:, ind_idx, 0]][0].item()\n", + " ax.text(\n", + " x=ds.position[start_frame, ind_idx, 0],\n", + " y=ds.position[start_frame, ind_idx, 1],\n", + " s=ds.individuals[ind_idx].item().split(\"_\")[1],\n", + " fontsize=8,\n", + " color=list_colors[ind_idx % len(list_colors)],\n", + " )\n", + "\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_xlim(-150, 4200) # frame size: 4096x2160\n", + "ax.set_ylim(-150, 2250) # frame size: 4096x2160\n", + "ax.set_xlabel(\"x (pixels)\")\n", + "ax.set_ylabel(\"y (pixels)\")\n", + "ax.set_title(Path(ds.source_file).stem)\n", + "ax.invert_yaxis()\n", + "\n", + "# Save plot as png\n", + "plt.savefig(\n", + " output_figures_dir / f\"{Path(ds.source_file).stem}_tracks.png\",\n", + " dpi=300,\n", + " bbox_inches=\"tight\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot histogram of trajectories' lengths" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "1f54d085-e579-4a5b-abea-d2cf00d8930b", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAHFCAYAAAAUpjivAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWo0lEQVR4nO3deVhUZf8G8HvYdxAEBhQRExXFHVNRAfdwDzNNSyk119xK3FLBVJLKJdc0Q0vRyr1UFDfSzERRM/cFRFPCTAFF9uf3hz/O6zhsA4MznO7Pdc31Ns955pzvOWd45/Y5m0IIIUBEREQkUwa6LoCIiIioIjHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMexQiR4/fowJEybA1dUVZmZmaNKkCTZv3lxkfyEE/Pz8oFAoMHbs2FIv58CBA2jdujUsLCxQtWpVBAcHIyUlRa3f1atX0bdvX1SpUgUWFhZo2bIldu3aVaplHDp0CO+99x7q1asHS0tLVKtWDb1798bp06cL7R8fH49OnTrBysoKdnZ2CAoKws2bN9Xq+eijj9C8eXPY2dnB3t4ebdq0wZYtWwpdx86dO8PV1RWmpqZwcnJChw4dsGfPnlLVr0ldALB48WIEBQXBw8MDCoUCAQEBGi1n3bp1UCgU0svIyAguLi4YMGAArl27ptY/ICBA42UUCA4OhpWVVYn9MjIyEBoaiiNHjmg0/0uXLuGdd95BrVq1YGZmhqpVq6JZs2YYO3Ys0tLSylSzth0/fhyhoaF49OiRrkvRufT0dISEhKBLly5wdHSEQqFAaGhomealD9u1tN9vqhgMO1SioKAgrF+/HrNnz8bevXvRokULvPXWW4iKiiq0//Lly3H9+nWNlhEbG4vAwEA4Oztj586dWLJkCQ4cOICOHTsiKytL6peYmIjWrVvjypUrWLVqFX788Uc4OjqiT58+2Lp1a4nLWblyJRITEzF+/Hjs2bMHS5YsQUpKClq1aoVDhw6p9L18+TICAgKQnZ2NH374Ad988w2uXr2Kdu3a4f79+1K//fv3Y/fu3ejbty9+/PFHbNy4EZ6enujXrx/mzJmjMs8HDx6gQYMGWLRoEfbv34+vvvoKxsbG6N69OzZs2FCqbVXaugBg1apVuHXrFjp06ABHR8dSzb8wkZGR+O2333DgwAGMHTsWu3btQtu2bfHw4UOVfitWrMCKFSvKvJzSyMjIQFhYmEZh58yZM2jevDkuXryIWbNmITo6GqtWrUL37t2xb98+/PvvvxVXsAaOHz+OsLAwhh08+1tZvXo1srKy0KdPn3LNi9uVIIiKsXv3bgFAREVFqbR37txZuLq6itzcXJX2hIQEYWVlJbZt2yYAiDFjxpRqOS1atBD169cXOTk5Utuvv/4qAIgVK1ZIbSNGjBBmZmbizp07Ultubq7w8vISbm5uIi8vr9jl/P3332pt6enpwtnZWXTs2FGlvV+/fqJq1aoiNTVVaktMTBTGxsYiJCREart//77Iz89Xm2/37t2FhYWFyMzMLLam7OxsUa1aNdGuXbti+2lalxBCZXs0aNBA+Pv7l2oZBSIjIwUAERcXp9IeFhYmAIhvvvlGo/kVZ8iQIcLS0rLEfvfv3xcAxOzZs0s978GDBwtLS0uRlpZW6PTC9p8ufPbZZwKASEhI0HUpOpefny/tl7Ls8+dpul0zMjLKtJzilPb7TRWDIztUrO3bt8PKygr9+vVTaX/33Xdx9+5d/P777yrt77//Pjp37ozXX3+91Mv466+/EBcXh3feeQdGRkZSu6+vL+rUqYPt27dLbb/++isaN26MatWqSW2GhoYIDAzE7du3cfLkyWKX5eTkpNZmZWWF+vXr4/bt21Jbbm4ufv75Z/Tt2xc2NjZSu7u7O9q3b69SU9WqVaFQKNTm++qrryIjI6PEUQNjY2PY2dmprHtRNKkLAAwMKuZP3MfHBwDw999/q7QXdhjrzp07eOONN2BtbQ07OzsMGjQIcXFxUCgUWLdundq8r1+/jm7dusHKygpubm748MMPpdG9xMREaYQqLCxMOrwWHBxcbL0PHjyAjY1NkYcRnt9/AQEB8Pb2xtGjR9GqVSuYm5ujWrVqmDlzJvLy8lQ+9++//2L06NGoVq0aTExMUKtWLcyYMUNlNLJg/mPHjsV3330HLy8vWFhYoHHjxvj555+lPqGhoZg8eTIASIcdFQqFNIL1/fffo0uXLnBxcYG5uTm8vLwwdepUPHnyRGVZBYdLituOBbKzszF37lzUq1cPpqamcHR0xLvvvqs2Qpifn4+IiAipn5OTEwYPHow7d+6o9KtZs2ah++LF70V+fj7mzp2LunXrwtzcHHZ2dmjUqBGWLFmiss0K+7vSVEnbtWbNmujRowe2bduGpk2bwszMDGFhYQCejVL7+fnByckJlpaWaNiwISIiIpCTk6O2nOjoaHTs2BG2trawsLCAl5cXwsPDi63t119/RdWqVdGjRw9pPx46dAgBAQFwcHCAubk5atSogb59+yIjI6Pc2+K/rOT/d6X/tD///BNeXl5qP8SNGjWSpvv6+gIAvv76a5w8eRIXL17UeBnPz/PF5fz666/S++zsbNjb26v1MzU1BQD88ccfaNWqlUbLT01NRXx8PDp06CC13bhxA0+fPi2yppiYGGRmZsLMzKzI+R4+fBiOjo6FBqz8/Hzk5+cjJSUFX331Fa5evYoFCxaUWKs26tKGhIQEAECdOnWK7ffkyRO0b98e//77LxYsWIDatWsjOjoa/fv3L7R/Tk4OevXqhaFDh+LDDz/EL7/8gk8++QS2traYNWsWXFxcEB0djddeew1Dhw7FsGHDAKDEQ3StW7fG7t27MWjQIIwYMQKvvvoqzM3Ni+yfnJyMAQMGYOrUqZgzZw52796NuXPn4uHDh1i2bBkAIDMzE+3bt8eNGzcQFhaGRo0a4ejRowgPD8fZs2exe/dulXnu3r0bcXFxmDNnDqysrBAREYHXX38dV65cQa1atTBs2DD8+++/WLp0KbZt2wYXFxcAQP369QEA165dQ7du3TBhwgRYWlri8uXLWLBgAU6ePKl2CLak7Qg8+w727t0bR48eRUhICHx9fXHr1i3Mnj0bAQEBOHXqlLSNRo0ahdWrV2Ps2LHo0aMHEhMTMXPmTBw5cgTx8fGoWrVqsdv/RREREQgNDcXHH38MPz8/5OTk4PLlyxVymKmk7Qo8Owfu0qVL+Pjjj+Hh4QFLS0sAz/7eBg4cCA8PD5iYmODcuXOYN28eLl++jG+++Ub6/Nq1azF8+HD4+/tj1apVcHJywtWrV6X/byvMDz/8gMGDB+O9997D0qVLYWhoiMTERHTv3h3t2rXDN998Azs7O/z111+Ijo5GdnY2LCwstL59/jN0PbRE+s3T01N07dpVrf3u3bsCgJg/f74QQog7d+4IW1tb8dVXX0l9UMrDWBs3bhQAxG+//aY27f333xcmJibS+z59+gg7OzuRnp6u0q9du3Yq9Whi0KBBwsjISJw6dUpqKziEtmnTJrX+8+fPFwDE3bt3i5znmjVrBACxZMmSQqd37dpVABAAhI2Njdi2bVupai1PXeU5jHXixAmRk5Mj0tPTRXR0tFAqlcLPz0/lsKMQQvj7+6ssY/ny5QKA2Lt3r0q/ESNGCAAiMjJSahsyZIgAIH744QeVvt26dRN169aV3pflkEZmZqbo06ePtM0NDQ1F06ZNxYwZM0RKSoraOgAQO3fuVGkfPny4MDAwELdu3RJCCLFq1apC612wYIEAIPbv3y+1ARDOzs4qh9GSk5OFgYGBCA8Pl9pKe7glPz9f5OTkiNjYWAFAnDt3TppW2u24adMmAUBs3bpVpV9cXJzK4eNLly4JAGL06NEq/X7//XcBQEyfPl1qc3d3F0OGDFGr98XvRY8ePUSTJk2KXcfnVeRhLHd3d2FoaCiuXLlS7Dzy8vJETk6O+Pbbb4WhoaH4999/hRDPDoPb2NiItm3bFns49PnDWJ9++qkwNDQUCxYsUOmzZcsWAUCcPXtWwzWkkvAwFpWouKHkgmkjR45E48aNMXz4cK0v5/n2sWPHIjU1FYMHD8bNmzfx999/Y+bMmTh+/DgAzQ/bzJw5Exs3bsSiRYvQvHnzUtdU3LS9e/dizJgxeOONN/DBBx8U2mfp0qU4efIkdu7cia5du6J///7YtGmTND0/Px+5ubnS68XDJ2WpqyhCCJVl5ebmqvVp1aoVjI2NYW1tjddeew1VqlTBzp07Szz0FhsbK33meW+99VaRtffs2VOlrVGjRrh161ap1uXF9RBCAHg28rd9+3ZcvHgRixYtwoABA3D//n3MmzcPXl5euHLlisp8rK2t0atXL5W2gQMHIj8/H7/88guAZ4cbLC0t8cYbb6j0KziMc/DgQZX29u3bw9raWnrv7OwMJyenUq/bzZs3MXDgQCiVShgaGsLY2Bj+/v4Anl1p9rzSbMeff/4ZdnZ26Nmzp8o2a9KkCZRKpXSY5/DhwyrrVeDVV1+Fl5eX2nqWxquvvopz585h9OjR2Ldvn86vhmvUqFGho5RnzpxBr1694ODgIG3zwYMHIy8vD1evXgXw7OTntLQ0jB49usS/PSEERowYgdmzZyMqKgohISEq05s0aQITExO8//77WL9+faFXWFLZMOxQsRwcHPDgwQO19oLzUOzt7bFlyxZER0cjIiICqampePTokTQcnZ2djUePHhV6jPv5ZQAocjnPH7bq2LEjIiMj8csvv+CVV16BUqnEtm3b8MknnwCAyrk8JQkLC8PcuXMxb948tUvkS6pJoVDAzs5Obdq+ffsQFBSEzp07Y+PGjUX+n5+npydatGiBXr164YcffkDHjh0xZswY5OfnAwDee+89GBsbS6+OHTuWq67ixMbGqizL2NgYiYmJKn2+/fZbxMXF4dChQxgxYgQuXbpUZGB53oMHD+Ds7KzWXlgbAFhYWKgdgjM1NUVmZmaJy0pMTFRbj9jYWJU+Xl5emDBhAjZs2ICkpCQsXLgQDx48wMyZM0usT6lUSutU8L9KpVJtHzs5OcHIyEhtHxXsuxfX7enTpyWu2+PHj9GuXTv8/vvvmDt3Lo4cOYK4uDhs27YNANTmUZrt+Pfff+PRo0cwMTFR227Jycn4559/VNa34PDP81xdXQv9LpZk2rRp+Pzzz3HixAkEBgbCwcEBHTt2xKlTpzSelzYUtm5JSUlo164d/vrrLyxZsgRHjx5FXFwcli9fDuB/27zg/Kbq1auXuJzs7Gx8//33aNCgAQIDA9Wmv/LKKzhw4ACcnJwwZswYvPLKK3jllVdUzmWisuE5O1Sshg0bYtOmTcjNzVX5V/z58+cBAN7e3ti/fz9yc3MLPVdmzZo1WLNmDbZv317k5aPe3t7SPLt166Yy7fz589L0AkOGDMGgQYNw7do1GBsbo3bt2ggPD4dCoUC7du1KtV5hYWEIDQ1FaGgopk+frjb9lVdegbm5ubSeL9ZUu3ZttR+Tffv2oU+fPvD398fWrVthYmJSqlqAZ//SjY6Oxv379+Hs7IzQ0FCVAFYwIlCWukrSvHlzxMXFqbS5urqqvPfy8pJOSm7fvj3y8vLw9ddfY8uWLWojG89zcHAo9KTx5ORkjWosDVdXV7X1qFu3bpH9FQoFJk6ciDlz5qidW/HiidfA/2ouCC0ODg74/fffIYRQCTwpKSnIzc3V+DyW4hw6dAh3797FkSNHpNEcAOU6x6Vq1apwcHBAdHR0odMLvnMF63vv3j21H/S7d++qrKeZmZnaSdAA8M8//6j0MzIywqRJkzBp0iQ8evQIBw4cwPTp09G1a1fcvn37pZ+bUtg/Snbs2IEnT55g27ZtcHd3l9rPnj2r0q/gfLEXT9YujKmpKQ4fPoyuXbuiU6dOiI6ORpUqVVT6tGvXDu3atUNeXh5OnTqFpUuXYsKECXB2dsaAAQPKsHYEcGSHSvD666/j8ePHavewWb9+PVxdXdGyZUsEBwfj8OHDai8A6NOnDw4fPoy2bdsWuYxq1arh1VdfxYYNG1QO15w4cQJXrlxBUFCQ2meMjIzg5eWF2rVrIzU1FatXr0bv3r1V/k+pKJ988ol0cuTs2bML7WNkZISePXti27ZtSE9Pl9qTkpJw+PBhtZr279+PPn36oG3bttixY4d0wnRpCCEQGxsLOzs76YelZs2a8PHxkV4FP9qa1lUa1tbWKsvy8fEpMahFRESgSpUqmDVrljQaVRh/f3+kp6dj7969Ku3F3ZSyJAXb9sXRDBMTE7X1KPjBvnfvXqHzunv3LtLS0tTCXXp6utqNKqOiomBgYAA/Pz8Az0YZHz9+jB07dqj0+/bbb6Xp2lq3gh/jF79XX331lcbLKNCjRw88ePAAeXl5atvt+e9cwYn7L94HKi4uDpcuXVJZz5o1a+KPP/5Q6Xf16lW1w4TPs7OzwxtvvIExY8bg33//VRtV1IaitmtxCtvmQgisWbNGpZ+vry9sbW2xatUq6bBpcZo2bYrY2FjcuXMHAQEBhd44FXh2lWnLli2lkaT4+PhS107qOLJDxQoMDETnzp0xatQopKWloXbt2ti0aROio6OxYcMGGBoaombNmqhZs2ahn69WrZrapchGRkbw9/dXOda/YMECdO7cGf369cPo0aORkpKCqVOnwtvbG++++67ULyUlBV988QXatGkDa2trXL58GRERETAwMJD+T6HAnDlzMGfOHBw8eFD61/AXX3yBWbNm4bXXXkP37t1x4sQJlc88PzoVFhaGFi1aoEePHpg6dSoyMzMxa9YsVK1aFR9++KHU79ixY+jTpw+USiWmT5+u9i+/+vXrS5eJ9+7dG40bN0aTJk3g4OCAu3fvYt26dYiNjcXy5ctLdfl5aesCgFOnTkk/HmlpaRBCSHd2btGiRanCYWGqVKmCadOmISQkBFFRUXj77bcL7TdkyBAsWrQIb7/9NubOnYvatWtj79692LdvH4CyXRpvbW0Nd3d37Ny5Ex07doS9vT2qVq1a5HcQeHZLhEePHqFv377w9vaGoaEhLl++jEWLFsHAwABTpkxR6e/g4IBRo0YhKSkJderUwZ49e7BmzRqMGjUKNWrUAAAMHjwYy5cvx5AhQ5CYmIiGDRvi2LFjmD9/Prp164ZOnTppvG4NGzYEACxZsgRDhgyBsbEx6tatC19fX1SpUgUjR47E7NmzYWxsjI0bN+LcuXMaL6PAgAEDsHHjRnTr1g3jx4/Hq6++CmNjY9y5cweHDx9G79698frrr6Nu3bp4//33sXTpUhgYGCAwMFC6GsvNzQ0TJ06U5vnOO+/g7bffxujRo9G3b1/cunULERERalfL9ezZE97e3vDx8YGjoyNu3bqFxYsXw93dHZ6enlK/vXv34smTJ1Kwv3jxovT97datW6lHgIrars+fQ/Wizp07w8TEBG+99RZCQkKQmZmJlStXqt1I08rKCl988QWGDRuGTp06Yfjw4XB2dsb169dx7tw56eq953l5eeHo0aPo1KkT/Pz8cODAAVSvXh2rVq3CoUOH0L17d9SoUQOZmZnSVV9l+T7Rc3R3bjRVFunp6WLcuHFCqVQKExMT0ahRo0KvBnoRirgaC0ChVwXt379ftGrVSpiZmQl7e3sxePBgtZsAPnjwQHTp0kU4OjoKY2NjUaNGDfHBBx+I+/fvq81v9uzZAoA4fPiw1FZwpU1RrxedOnVKdOzYUVhYWAgbGxvRp08fcf369UKXU9Tr+eUvWLBAtGjRQlSpUkUYGhoKBwcH0bVrV/Hzzz+XsDU1r0uI/12ZU9jr+SuhilLUTQWFEOLp06eiRo0awtPTU7q55ItX3QghRFJSkggKChJWVlbC2tpa9O3bV+zZs0ftiqeibrpWsH2fd+DAAdG0aVNhamoqABR6BdDz9u3bJ9577z1Rv359YWtrK4yMjISLi4sICgpSuwrQ399fNGjQQBw5ckT4+PgIU1NT4eLiIqZPn6529dmDBw/EyJEjhYuLizAyMhLu7u5i2rRpajeSLOpvobCrl6ZNmyZcXV2FgYGByvfn+PHjonXr1sLCwkI4OjqKYcOGifj4+EKvaivtdszJyRGff/65aNy4sTAzMxNWVlaiXr16YsSIEeLatWtSv7y8PLFgwQJRp04dYWxsLKpWrSrefvttcfv2bZX55efni4iICFGrVi1hZmYmfHx8xKFDh9S+F1988YXw9fUVVatWFSYmJqJGjRpi6NChIjExUW37FPX91fTGi0VtV3d3d9G9e/dCP/PTTz9J26ZatWpi8uTJYu/evWp/10IIsWfPHuHv7y8sLS2FhYWFqF+/vsrVVoXtlzt37oh69eqJmjVrihs3bojffvtNvP7668Ld3V2YmpoKBwcH4e/vL3bt2qXRupI6hRClGHcjItKi+fPn4+OPP0ZSUlKpTux8mQICAvDPP/8Ue48UIqpceBiLiCpUwTB+vXr1kJOTg0OHDuHLL7/E22+/rXdBh4jkiWGHiCqUhYUFFi1ahMTERGRlZaFGjRqYMmUKPv74Y12XRpVcwZ3Ii1Oa8+BI/ngYi4iIKqXQ0FDpOVZFSUhIKPbkdfpvYNghIqJK6e7du7h7926xfRo1aqTRPa9Inhh2iIiISNZ4U0EiIiKSNZ65hWcnud29exfW1tYaP0SRiIiIdEMIgfT0dLi6uhZ7k1KGHTw77uvm5qbrMoiIiKgMbt++XeytLBh28L8H3t2+fVu6rT8RERHpt7S0NLi5uRX76A+AYQfA/x74ZmNjw7BDRERUyZR0CgpPUCYiIiJZY9ghIiIiWWPYISIiIlnjOTtERKRzeXl5yMnJ0XUZpGeMjY1haGhY7vkw7BARkc4IIZCcnIxHjx7puhTSU3Z2dlAqleW6Dx7DDhER6UxB0HFycoKFhQVv7EoSIQQyMjKQkpICAHBxcSnzvBh2iIhIJ/Ly8qSg4+DgoOtySA+Zm5sDAFJSUuDk5FTmQ1o8QZmIiHSi4BwdCwsLHVdC+qzg+1Gec7oYdoiISKd46IqKo43vB8MOERERyRrDDhERkZ4KDQ1FkyZNpPfBwcHo06dPueapjXlUNjxBmYiI9M6imKsvbVkTO9fR+DPBwcFYv349AMDIyAhubm4ICgpCWFgYLC0ttV2iZMmSJRBClKpvYmIiPDw8cObMGZXApMk85IJhh4iIqAxee+01REZGIicnB0ePHsWwYcPw5MkTrFy5UqVfTk4OjI2NtbJMW1tbvZhHZcOwU9GePNF1BURE+ikrC8jPB/Lynr2eJ/JfXh0vLrs0hICpiQmUjo4AgIH9++PwoUPYsWMHnB0dsWPXLowbOxZz589HYmIi8rKzkZaWhslTpmDHzp3IzMyET/PmWPTFF2jcuLE0208XLMCiJUuQkZGBN/v1g2PVqio1Br/3Hh49eoQd27YBAPLz8/HZ559jzdq1uH37NpydnTFi+HDMmD4dHh4eAICmTZsCAPz9/HDk0CEEDx36bB47dgAAsrKyMHnyZGzevBlpaWnw8fHBokWL0KJFCwDAkSNH0L59exw4cABTpkzBxYsX0aRJE0RGRqJu3bqabzsdYNipaFZWuq6AiEg/ubsDq1YBT5+qT7v3+OXVcaYMy3rwAHj8GDhzRmoyf/wYOU+fAsnJuH71Kn5YuxZbP/kEhgYGwJkz6D58OOxtbLDn889ha2WFr7ZtQ8cOHXB161bY29rih5gYzA4NxfKQELRr0gTf7d2LL1esQC1X1/8t54XlTlu6FGt27MCiiRPRtkkT3PvnH1xOTATOnMHJdevwanAwDixfjga1asHE2Fil3gIhISHYunUr1q9fD3d3d0RERKBr1664fv067O3tpX4zZszAF198AUdHR4wcORLvvfcefv31V823nQ4w7BAREZXTyQsXEBUdjY7/PxqSnZOD7+bMgWOVKgCAQ3FxOH/9OlL274epiQkA4PMJE7AjNhZbDh7E+0FBWLxpE97r1QvD/v/k4bmjRuHAyZPIzMoqdJnpT55gyebNWDZ5Mob06AEAeKV6dbT9//NzCpbtYGsLZcEI0QsKDrutW7cOgYGBAIA1a9YgJiYGa9euxeTJk6W+8+bNg7+/PwBg6tSp6N69OzIzM2FmZlbWzfbSMOxUtMcv8V8nRESVSVYWcO8eULMm8OIP5r/XXl4dTT01/4yDA36OjoZVQAByc3ORk5OD3r16YemKFVixciXca9aEY4cOUvfTBw/i8dOncOjSRWU2T58+xY2cHKBpU1y6fRsjJ00C/v+wEwC07tABh2Nj/9fm4AAYGj7rf/IksrKz0TE4GPj/Q1Yq/j/soF494LkTlJ9348YN5OTkoE2bNlKbsbExXn31VVy6dEmlb6NGjaT/Lnh0Q0pKCmrUqFHS1tI5hp2KVoFn5RMRVWqGhoCBwbP/ffExAIqXeGeUsjyCQKFA+/btsXLlShgbG8PV1fV/JyEbGDy7Iuu5+ebjWUA4cuSI2qzs7Oz+17dgexQwMFCtUaF49jI0hHnBaRKFbb/nP1PUdEC6KuvFG/cJIdTanj/JumBafv5LPLeqHHifHSIiojKwtLRE7dq14e7uXuLVVs2aNUNycjKMjIxQu3ZtlVfV/z/E5OXlhRMnTqh87sX3z/P09IS5uTkOHjxY6HST/z9cllfMCdi1a9eGiYkJjh07JrXl5OTg1KlT8PLyKnadKhOO7BAREVWwTp06oXXr1ujTpw8WLFiAunXr4u7du9izZw/69OkDHx8fjB8/HkOGDIGPjw/atm2LjRs34sKFC6hVq1ah8zQzM8OUKVMQEhICExMTtGnTBvfv38eFCxcwdOhQODk5wdzcHNHR0ahevTrMzMzULju3tLTEqFGjMHnyZNjb26NGjRqIiIhARkYGhg4d+jI2zUvBsENERFTBFAoF9uzZgxkzZuC9997D/fv3oVQq4efnB2dnZwBA//79cePGDUyZMgWZmZno27cvRo0ahX379hU535kzZ8LIyAizZs3C3bt34eLigpEjRwJ4drPDL7/8EnPmzMGsWbPQrl27Qg+jffrpp8jPz8c777yD9PR0+Pj4YN++fahScM6PDCjEf+02ioVIS0uDra0tUlNTYWNjo+tyiIj+EzIzM5GQkAAPD49KcUUP6UZx35PS/n7znB0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYYdIiIikjWdhp1ffvkFPXv2hKurKxQKhfQE1gJCCISGhsLV1RXm5uYICAjAhQsXVPpkZWXhgw8+QNWqVWFpaYlevXrhzp07L3EtiIiISJ/pNOw8efIEjRs3xrJlywqdHhERgYULF2LZsmWIi4uDUqlE586dkZ6eLvWZMGECtm/fjs2bN+PYsWN4/PgxevToUewdI4mIiOi/Q6c3FQwMDJSesvoiIQQWL16MGTNmICgoCACwfv16ODs7IyoqCiNGjEBqairWrl2L7777Dp06dQIAbNiwAW5ubjhw4AC6du360taFiIiI9JPenrOTkJCA5ORkdHnuCbGmpqbw9/fH8ePHAQCnT59GTk6OSh9XV1d4e3tLfQqTlZWFtLQ0lRcRERHJk96GneTkZACQbqNdwNnZWZqWnJwMExMTtVtaP9+nMOHh4bC1tZVebm5uWq6eiIhIVXBwMPr06VPu+Tx/jmtiYiIUCgXOnj1bps8XpizzLI2aNWti8eLFWp1naen9s7FK89j5F5XUZ9q0aZg0aZL0Pi0tjYGHiIgq1JIlS6DtJzS5ubnh3r170pPTS+PevXuyeu5VaejtyI5SqQQAtRGalJQUabRHqVQiOzsbDx8+LLJPYUxNTWFjY6PyIiIiqki2traws7PT6jwNDQ2hVCphZFT6sQulUglTU1Ot1qHv9DbseHh4QKlUIiYmRmrLzs5GbGwsfH19AQDNmzeHsbGxSp979+7hzz//lPoQERFpW0BAAMaNG4eQkBDY29tDqVQiNDS02M+8eBirNPO4du0a/Pz8YGZmhvr166v83gGqh5zy8/NRvXp1rFq1SqVPfHw8FAoFbt68CUD9MNbJkyfRtGlTmJmZwcfHB2fOnFH5/Lp169RC2o4dO1SOoNy4cQO9e/eGs7MzrKys0KJFCxw4cKDY7REaGooaNWrA1NQUrq6uGDduXLH9y0Onh7EeP36M69evS+8TEhJw9uxZ2Nvbo0aNGpgwYQLmz58PT09PeHp6Yv78+bCwsMDAgQMBPEvJQ4cOxYcffggHBwfY29vjo48+QsOGDaWrs4iIqBJ68uTlLcvSskwfW79+PSZNmoTff/8dv/32G4KDg9GmTRt07txZK/PIz89HUFAQqlatihMnTiAtLQ0TJkwocl4GBgYYMGAANm7ciJEjR0rtUVFRaN26NWrVqqX2mSdPnqBHjx7o0KEDNmzYgISEBIwfP16j7QA8+z3v1q0b5s6dCzMzM6xfvx49e/bElStXUKNGDbX+W7ZswaJFi7B582Y0aNAAycnJOHfunMbLLTWhQ4cPHxYA1F5DhgwRQgiRn58vZs+eLZRKpTA1NRV+fn7i/PnzKvN4+vSpGDt2rLC3txfm5uaiR48eIikpSaM6UlNTBQCRmpqqrVUjIqISPH36VFy8eFE8ffpUfSLw8l5l4O/vL9q2bavS1qJFCzFlypQiPzNkyBDRu3fvUs9j3759wtDQUNy+fVuavnfvXgFAbN++XQghREJCggAgzpw5I4QQIj4+XigUCpGYmCiEECIvL09Uq1ZNLF++XJrH85//6quvhL29vXjy5Ik0feXKlSrzjIyMFLa2tip1bt++XZQUIerXry+WLl0qvXd3dxeLFi0SQgjxxRdfiDp16ojs7Oxi5yFE8d+T0v5+6/QwVkBAAIQQaq9169YBeDbUFhoainv37iEzMxOxsbHw9vZWmYeZmRmWLl2KBw8eICMjAz/99BNPNiYiogrXqFEjlfcuLi5ISUnR2jwuXbqEGjVqoHr16tL01q1bFzu/pk2bol69eti0aRMAIDY2FikpKXjzzTcL7X/p0iU0btwYFhYWpV5GYZ48eYKQkBDUr18fdnZ2sLKywuXLl5GUlFRo/379+uHp06eoVasWhg8fju3btyM3N1fj5ZaW3p6zQ0RE/2GPH7+8VxkZGxurvFcoFMjPz9faPEQhV26VdDUyAAwaNAhRUVEAnh3C6tq1a5FXaxW2jBcZGBio9cvJyVF5P3nyZGzduhXz5s3D0aNHcfbsWTRs2BDZ2dmFztPNzQ1XrlzB8uXLYW5ujtGjR8PPz09tvtqi95eeExHRf1AZz6ORk/r16yMpKQl3796Fq6srAOC3334r8XMDBw7Exx9/jNOnT2PLli1YuXJlscv47rvv8PTpU5ibmwMATpw4odLH0dER6enpePLkCSz/f7+8eA+eo0ePIjg4GK+//jqAZ+fwJCYmFlunubk5evXqhV69emHMmDGoV68ezp8/j2bNmpW4jpriyA4REZEe6tSpE+rWrYvBgwfj3LlzOHr0KGbMmFHi5zw8PODr64uhQ4ciNzcXvXv3LrLvwIEDYWBggKFDh+LixYvYs2cPPv/8c5U+LVu2hIWFBaZPn47r168jKipKOt2kQO3atbFt2zacPXsW586dw8CBA4sd5Vq3bh3Wrl2LP//8Ezdv3sR3330Hc3NzuLu7l7h+ZcGwQ0REpIcMDAywfft2ZGVl4dVXX8WwYcMwb968Un120KBBOHfuHIKCgqQRm8JYWVnhp59+wsWLF9G0aVPMmDEDCxYsUOljb2+PDRs2YM+ePWjYsCE2bdqkdon8okWLUKVKFfj6+qJnz57o2rVrsSM0dnZ2WLNmDdq0aYNGjRrh4MGD+Omnn+Dg4FCq9dOUQpTmgJ3MpaWlwdbWFqmpqbzBIBHRS5KZmYmEhAR4eHjAzMxM1+WQnirue1La32+O7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQEZFO8ToZKo42vh8MO0REpBMFdw/OyMjQcSWkzwq+Hy/ebVoTvIMyERHphKGhIezs7KRnQVlYWJTqcQj03yCEQEZGBlJSUmBnZwdDQ8Myz4thh4iIdEapVAKAxg/QpP8OOzs76XtSVgw7RESkMwqFAi4uLnBycqqwh0BS5WVsbFyuEZ0CDDtERKRzhoaGWvlRIyoMT1AmIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIlnT67CTm5uLjz/+GB4eHjA3N0etWrUwZ84c5OfnS32EEAgNDYWrqyvMzc0REBCACxcu6LBqIiIi0id6HXYWLFiAVatWYdmyZbh06RIiIiLw2WefYenSpVKfiIgILFy4EMuWLUNcXByUSiU6d+6M9PR0HVZORERE+kKvw85vv/2G3r17o3v37qhZsybeeOMNdOnSBadOnQLwbFRn8eLFmDFjBoKCguDt7Y3169cjIyMDUVFROq6eiIiI9IFeh522bdvi4MGDuHr1KgDg3LlzOHbsGLp16wYASEhIQHJyMrp06SJ9xtTUFP7+/jh+/LhOaiYiIiL9YqTrAoozZcoUpKamol69ejA0NEReXh7mzZuHt956CwCQnJwMAHB2dlb5nLOzM27dulXkfLOyspCVlSW9T0tLq4DqiYiISB/o9cjO999/jw0bNiAqKgrx8fFYv349Pv/8c6xfv16ln0KhUHkvhFBre154eDhsbW2ll5ubW4XUT0RERLqn12Fn8uTJmDp1KgYMGICGDRvinXfewcSJExEeHg4AUCqVAP43wlMgJSVFbbTnedOmTUNqaqr0un37dsWtBBEREemUXoedjIwMGBiolmhoaChdeu7h4QGlUomYmBhpenZ2NmJjY+Hr61vkfE1NTWFjY6PyIiIiInnS63N2evbsiXnz5qFGjRpo0KABzpw5g4ULF+K9994D8Ozw1YQJEzB//nx4enrC09MT8+fPh4WFBQYOHKjj6omIiEgf6HXYWbp0KWbOnInRo0cjJSUFrq6uGDFiBGbNmiX1CQkJwdOnTzF69Gg8fPgQLVu2xP79+2Ftba3DyomIiEhfKIQQQtdF6FpaWhpsbW2RmprKQ1pERESVRGl/v/X6nB0iIiKi8mLYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWdPrZ2MREREBwKKYq7ougcphYuc6Ol0+R3aIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWNA478fHxOH/+vPR+586d6NOnD6ZPn47s7GytFkdERERUXhqHnREjRuDq1asAgJs3b2LAgAGwsLDAjz/+iJCQEK0XSERERFQeGoedq1evokmTJgCAH3/8EX5+foiKisK6deuwdetWbddHREREVC4ahx0hBPLz8wEABw4cQLdu3QAAbm5u+Oeff7RbHREREVE5aRx2fHx8MHfuXHz33XeIjY1F9+7dAQAJCQlwdnbWeoFERERE5aFx2Fm8eDHi4+MxduxYzJgxA7Vr1wYAbNmyBb6+vlovkIiIiKg8jDT9QKNGjVSuxirw2WefwdDQUCtFEREREWlLme6z8+jRI3z99deYNm0a/v33XwDAxYsXkZKSotXiiIiIiMpL45GdP/74Ax07doSdnR0SExMxfPhw2NvbY/v27bh16xa+/fbbiqiTiIiIqEw0HtmZNGkS3n33XVy7dg1mZmZSe2BgIH755RetFkdERERUXhqHnbi4OIwYMUKtvVq1akhOTtZKUURERETaonHYMTMzQ1pamlr7lStX4OjoqJWiiIiIiLRF47DTu3dvzJkzBzk5OQAAhUKBpKQkTJ06FX379tV6gURERETloXHY+fzzz3H//n04OTnh6dOn8Pf3R+3atWFtbY158+ZVRI1EREREZabx1Vg2NjY4duwYDh06hPj4eOTn56NZs2bo1KlTRdRHREREVC4ah50CHTp0QIcOHbRZCxEREZHWlSrsfPnll3j//fdhZmaGL7/8sti+48aN00phRERERNqgEEKIkjp5eHjg1KlTcHBwgIeHR9EzUyhw8+ZNrRb4MqSlpcHW1hapqamwsbHRdTlERPSCRTFXdV0ClcPEznUqZL6l/f0u1chOQkJCof9NREREpO80uhorJycHtWrVwsWLFyuqHiIiIiKt0ijsGBsbIysrCwqFoqLqISIiItIqje+z88EHH2DBggXIzc2tiHqIiIiItErjS89///13HDx4EPv370fDhg1haWmpMn3btm1aK46IiIiovDQOO3Z2dnwsBBEREVUaGoedyMjIiqijSH/99RemTJmCvXv34unTp6hTpw7Wrl2L5s2bAwCEEAgLC8Pq1avx8OFDtGzZEsuXL0eDBg1eap1ERESknzQ+Z6fA/fv3cezYMfz666+4f/++NmuSPHz4EG3atIGxsTH27t2Lixcv4osvvoCdnZ3UJyIiAgsXLsSyZcsQFxcHpVKJzp07Iz09vUJqIiIiospF45GdJ0+e4IMPPsC3336L/Px8AIChoSEGDx6MpUuXwsLCQmvFLViwAG5ubiqjSTVr1pT+WwiBxYsXY8aMGQgKCgIArF+/Hs7OzoiKisKIESO0VgsRERFVThqP7EyaNAmxsbH46aef8OjRIzx69Ag7d+5EbGwsPvzwQ60Wt2vXLvj4+KBfv35wcnJC06ZNsWbNGml6QkICkpOT0aVLF6nN1NQU/v7+OH78uFZrISIiospJ47CzdetWrF27FoGBgbCxsYGNjQ26deuGNWvWYMuWLVot7ubNm1i5ciU8PT2xb98+jBw5EuPGjcO3334LAEhOTgYAODs7q3zO2dlZmlaYrKwspKWlqbyIiIhInjQ+jJWRkaEWLgDAyckJGRkZWimqQH5+Pnx8fDB//nwAQNOmTXHhwgWsXLkSgwcPlvq9eJNDIUSxNz4MDw9HWFiYVmslIiIi/aTxyE7r1q0xe/ZsZGZmSm1Pnz5FWFgYWrdurdXiXFxcUL9+fZU2Ly8vJCUlAQCUSiUAqI3ipKSkFBrICkybNg2pqanS6/bt21qtm4iIiPSHxiM7ixcvRmBgIKpXr47GjRtDoVDg7NmzMDMzw759+7RaXJs2bXDlyhWVtqtXr8Ld3R3As6exK5VKxMTEoGnTpgCA7OxsxMbGYsGCBUXO19TUFKamplqtlYiIiPSTxmGnYcOGuHbtGjZs2IDLly9DCIEBAwZg0KBBMDc312pxEydOhK+vL+bPn48333wTJ0+exOrVq7F69WoAzw5fTZgwAfPnz4enpyc8PT0xf/58WFhYYODAgVqthYiIiConjcPOL7/8Al9fXwwfPlylPTc3F7/88gv8/Py0VlyLFi2wfft2TJs2DXPmzIGHhwcWL16MQYMGSX1CQkLw9OlTjB49Wrqp4P79+2Ftba21OoiIiKjyUgghhCYfMDQ0xL179+Dk5KTS/uDBAzg5OSEvL0+rBb4MaWlpsLW1RWpqKmxsbHRdDhERvWBRzFVdl0DlMLFznQqZb2l/vzU+QbmoK50ePHig9lBQIiIiIl0r9WGsgjsUKxQKBAcHq5zgm5eXhz/++AO+vr7ar5CIiIioHEoddmxtbQE8G9mxtrZWORnZxMQErVq1UjuPh4iIiEjXSh12Cp5PVbNmTUyePFmrz8AiIiIiqigan7MzePBg/PXXX2rt165dQ2JiojZqIiIiItIajcNOcHBwoQ/Z/P333xEcHKyNmoiIiIi0RuOwc+bMGbRp00atvVWrVjh79qw2aiIiIiLSGo3DjkKhQHp6ulp7ampqpbzHDhEREcmbxmGnXbt2CA8PVwk2eXl5CA8PR9u2bbVaHBEREVF5afy4iIiICPj5+aFu3bpo164dAODo0aNIS0vDoUOHtF4gERERUXloPLJTv359/PHHH3jzzTeRkpKC9PR0DB48GJcvX4a3t3dF1EhERERUZhqP7ACAq6sr5s+fr+1aiIiIiLRO45Ed4Nlhq7fffhu+vr7SPXe+++47HDt2TKvFEREREZWXxmFn69at6Nq1K8zNzREfH4+srCwAQHp6Okd7iIiISO9oHHbmzp2LVatWYc2aNTA2NpbafX19ER8fr9XiiIiIiMpL47Bz5coV+Pn5qbXb2Njg0aNH2qiJiIiISGs0DjsuLi64fv26WvuxY8dQq1YtrRRFREREpC0ah50RI0Zg/Pjx+P3336FQKHD37l1s3LgRH330EUaPHl0RNRIRERGVmcaXnoeEhCA1NRXt27dHZmYm/Pz8YGpqio8++ghjx46tiBqJiIiIyqxM99mZN28eZsyYgYsXLyI/Px/169eHlZWVtmsjIiIiKrcyhR0AsLCwgI+PjzZrISIiItK6UoWdoKAgrFu3DjY2NggKCiq2r5WVFRo0aICRI0fC1tZWK0USERERlVWpwo6trS0UCoX038XJysrCqlWr8Ouvv2LXrl3lr5CIiIioHEoVdiIjIwv976JcvHgRLVq0KHtVRERERFpSpmdjlaRu3bo4fvx4RcyaiIiISCNlOkE5Li4OP/74I5KSkpCdna0ybdu2bTA0NETjxo21UiARERFReWg8srN582a0adMGFy9exPbt25GTk4OLFy/i0KFDPCGZiIiI9I7GYWf+/PlYtGgRfv75Z5iYmGDJkiW4dOkS3nzzTdSoUaMiaiQiIiIqM43Dzo0bN9C9e3cAgKmpKZ48eQKFQoGJEydi9erVWi+QiIiIqDw0Djv29vZIT08HAFSrVg1//vknAODRo0fIyMjQbnVERERE5aTxCcrt2rVDTEwMGjZsiDfffBPjx4/HoUOHEBMTg44dO1ZEjURERERlpnHYWbZsGTIzMwEA06ZNg7GxMY4dO4agoCDMnDlT6wUSERERlYdGYSc3Nxc//fQTunbtCgAwMDBASEgIQkJCKqQ4IiIiovLS6JwdIyMjjBo1CllZWRVVDxEREZFWaXyCcsuWLXHmzJmKqIWIiIhI6zQ+Z2f06NH48MMPcefOHTRv3hyWlpYq0xs1aqS14oiIiIjKS+Ow079/fwDAuHHjpDaFQgEhBBQKBfLy8rRXHREREVE5aRx2EhISKqIOIiIiogqhcdi5desWfH19YWSk+tHc3FwcP34c7u7uWiuOiIiIqLw0PkG5ffv2+Pfff9XaU1NT0b59e60URURERKQtGoedgnNzXvTgwQO1k5WJiIiIdK3Uh7GCgoIAPDsZOTg4GKamptK0vLw8/PHHH/D19dV+hURERETlUOqwY2trC+DZyI61tTXMzc2laSYmJmjVqhWGDx+u/QqJiIiIyqHUYScyMhIAULNmTXz00Uc8ZEVERESVgsZXY82ePbsi6iAiIiKqEBqfoExERERUmTDsEBERkawx7BAREZGsMewQERGRrGl8gjIAHDx4EAcPHkRKSgry8/NVpn3zzTdaKYyIiIhIGzQOO2FhYZgzZw58fHzg4uJS6N2UiYiIiPSFxmFn1apVWLduHd55552KqIeIiIhIqzQ+Zyc7O5uPhSAiIqJKQ+OwM2zYMERFRVVELURERERap/FhrMzMTKxevRoHDhxAo0aNYGxsrDJ94cKFWiuOiIiIqLw0Djt//PEHmjRpAgD4888/VabxZGUiIiLSNxqHncOHD1dEHaUSHh6O6dOnY/z48Vi8eDGAZ09hDwsLw+rVq/Hw4UO0bNkSy5cvR4MGDXRWJxEREemPSnNTwbi4OKxevRqNGjVSaY+IiMDChQuxbNkyxMXFQalUonPnzkhPT9dRpURERKRPKkXYefz4MQYNGoQ1a9agSpUqUrsQAosXL8aMGTMQFBQEb29vrF+/HhkZGTyJmoiIiABUkrAzZswYdO/eHZ06dVJpT0hIQHJyMrp06SK1mZqawt/fH8ePH3/ZZRIREZEeKtPjIl6mzZs3Iz4+HnFxcWrTkpOTAQDOzs4q7c7Ozrh161aR88zKykJWVpb0Pi0tTUvVEhERkb7R65Gd27dvY/z48diwYQPMzMyK7PfiVWBCiGKvDAsPD4etra30cnNz01rNREREpF/0OuycPn0aKSkpaN68OYyMjGBkZITY2Fh8+eWXMDIykkZ0CkZ4CqSkpKiN9jxv2rRpSE1NlV63b9+u0PUgIiIi3dHrw1gdO3bE+fPnVdreffdd1KtXD1OmTEGtWrWgVCoRExODpk2bAnj2OIvY2FgsWLCgyPmamprC1NS0QmsnIiIi/aDXYcfa2hre3t4qbZaWlnBwcJDaJ0yYgPnz58PT0xOenp6YP38+LCwsMHDgQF2UTERERHpGr8NOaYSEhODp06cYPXq0dFPB/fv3w9raWtelERERkR5QCCGErovQtbS0NNja2iI1NRU2Nja6LoeIiF6wKOaqrkugcpjYuU6FzLe0v996fYIyERERUXkx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGs6XXYCQ8PR4sWLWBtbQ0nJyf06dMHV65cUekjhEBoaChcXV1hbm6OgIAAXLhwQUcVExERkb7R67ATGxuLMWPG4MSJE4iJiUFubi66dOmCJ0+eSH0iIiKwcOFCLFu2DHFxcVAqlejcuTPS09N1WDkRERHpCyNdF1Cc6OholfeRkZFwcnLC6dOn4efnByEEFi9ejBkzZiAoKAgAsH79ejg7OyMqKgojRozQRdlERESkR/R6ZOdFqampAAB7e3sAQEJCApKTk9GlSxepj6mpKfz9/XH8+HGd1EhERET6Ra9Hdp4nhMCkSZPQtm1beHt7AwCSk5MBAM7Ozip9nZ2dcevWrSLnlZWVhaysLOl9WlpaBVRMRERE+qDSjOyMHTsWf/zxBzZt2qQ2TaFQqLwXQqi1PS88PBy2trbSy83NTev1EhERkX6oFGHngw8+wK5du3D48GFUr15dalcqlQD+N8JTICUlRW2053nTpk1Damqq9Lp9+3bFFE5EREQ6p9dhRwiBsWPHYtu2bTh06BA8PDxUpnt4eECpVCImJkZqy87ORmxsLHx9fYucr6mpKWxsbFReREREJE96fc7OmDFjEBUVhZ07d8La2loawbG1tYW5uTkUCgUmTJiA+fPnw9PTE56enpg/fz4sLCwwcOBAHVdPRERE+kCvw87KlSsBAAEBASrtkZGRCA4OBgCEhITg6dOnGD16NB4+fIiWLVti//79sLa2fsnVEhERkT7S67AjhCixj0KhQGhoKEJDQyu+ICIiIqp09PqcHSIiIqLyYtghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWWPYISIiIllj2CEiIiJZY9ghIiIiWTPSdQFytyjmqq5LoHKa2LmOrksgIqJy4MgOERERyRrDDhEREckaD2MRkezxcDLRfxtHdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNaMdF0Akb5bFHNV1yUQEVE5cGSHiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZE02YWfFihXw8PCAmZkZmjdvjqNHj+q6JCIiItIDsgg733//PSZMmIAZM2bgzJkzaNeuHQIDA5GUlKTr0oiIiEjHZBF2Fi5ciKFDh2LYsGHw8vLC4sWL4ebmhpUrV+q6NCIiItKxSh92srOzcfr0aXTp0kWlvUuXLjh+/LiOqiIiIiJ9UelvKvjPP/8gLy8Pzs7OKu3Ozs5ITk4u9DNZWVnIysqS3qempgIA0tLStF5f5pPHWp8nERFRZVIRv6/Pz1cIUWy/Sh92CigUCpX3Qgi1tgLh4eEICwtTa3dzc6uQ2oiIiP7Lplfw/NPT02Fra1vk9EofdqpWrQpDQ0O1UZyUlBS10Z4C06ZNw6RJk6T3+fn5+Pfff+Hg4FBkQCqLtLQ0uLm54fbt27CxsdHafPWJ3NdR7usHyH8duX6Vn9zXketXdkIIpKenw9XVtdh+lT7smJiYoHnz5oiJicHrr78utcfExKB3796FfsbU1BSmpqYqbXZ2dhVWo42NjSy/wM+T+zrKff0A+a8j16/yk/s6cv3KprgRnQKVPuwAwKRJk/DOO+/Ax8cHrVu3xurVq5GUlISRI0fqujQiIiLSMVmEnf79++PBgweYM2cO7t27B29vb+zZswfu7u66Lo2IiIh0TBZhBwBGjx6N0aNH67oMFaamppg9e7baITM5kfs6yn39APmvI9ev8pP7OnL9Kp5ClHS9FhEREVElVulvKkhERERUHIYdIiIikjWGHSIiIpI1hh0iIiKSNYadCrRixQp4eHjAzMwMzZs3x9GjR3VdUpmEh4ejRYsWsLa2hpOTE/r06YMrV66o9AkODoZCoVB5tWrVSkcVayY0NFStdqVSKU0XQiA0NBSurq4wNzdHQEAALly4oMOKNVezZk21dVQoFBgzZgyAyrf/fvnlF/Ts2ROurq5QKBTYsWOHyvTS7LOsrCx88MEHqFq1KiwtLdGrVy/cuXPnJa5F8Ypbx5ycHEyZMgUNGzaEpaUlXF1dMXjwYNy9e1dlHgEBAWr7dcCAAS95TQpX0j4szXdSn/dhSetX2N+jQqHAZ599JvXR5/1Xmt8Fffo7ZNipIN9//z0mTJiAGTNm4MyZM2jXrh0CAwORlJSk69I0FhsbizFjxuDEiROIiYlBbm4uunTpgidPnqj0e+2113Dv3j3ptWfPHh1VrLkGDRqo1H7+/HlpWkREBBYuXIhly5YhLi4OSqUSnTt3Rnp6ug4r1kxcXJzK+sXExAAA+vXrJ/WpTPvvyZMnaNy4MZYtW1bo9NLsswkTJmD79u3YvHkzjh07hsePH6NHjx7Iy8t7WatRrOLWMSMjA/Hx8Zg5cybi4+Oxbds2XL16Fb169VLrO3z4cJX9+tVXX72M8ktU0j4ESv5O6vM+LGn9nl+ve/fu4ZtvvoFCoUDfvn1V+unr/ivN74Je/R0KqhCvvvqqGDlypEpbvXr1xNSpU3VUkfakpKQIACI2NlZqGzJkiOjdu7fuiiqH2bNni8aNGxc6LT8/XyiVSvHpp59KbZmZmcLW1lasWrXqJVWofePHjxevvPKKyM/PF0JU7v0HQGzfvl16X5p99ujRI2FsbCw2b94s9fnrr7+EgYGBiI6Ofmm1l9aL61iYkydPCgDi1q1bUpu/v78YP358xRanBYWtX0nfycq0D0uz/3r37i06dOig0lZZ9p8Q6r8L+vZ3yJGdCpCdnY3Tp0+jS5cuKu1dunTB8ePHdVSV9qSmpgIA7O3tVdqPHDkCJycn1KlTB8OHD0dKSoouyiuTa9euwdXVFR4eHhgwYABu3rwJAEhISEBycrLKvjQ1NYW/v3+l3ZfZ2dnYsGED3nvvPZUH31bm/fe80uyz06dPIycnR6WPq6srvL29K+1+TU1NhUKhUHvO38aNG1G1alU0aNAAH330UaUakSzuOymnffj3339j9+7dGDp0qNq0yrL/Xvxd0Le/Q9ncQVmf/PPPP8jLy1N76rqzs7Pa09krGyEEJk2ahLZt28Lb21tqDwwMRL9+/eDu7o6EhATMnDkTHTp0wOnTp/X+rqAtW7bEt99+izp16uDvv//G3Llz4evriwsXLkj7q7B9eevWLV2UW247duzAo0ePEBwcLLVV5v33otLss+TkZJiYmKBKlSpqfSrj32hmZiamTp2KgQMHqjxocdCgQfDw8IBSqcSff/6JadOm4dy5c9JhTH1W0ndSTvtw/fr1sLa2RlBQkEp7Zdl/hf0u6NvfIcNOBXr+X83Asy/Ei22VzdixY/HHH3/g2LFjKu39+/eX/tvb2xs+Pj5wd3fH7t271f6A9U1gYKD03w0bNkTr1q3xyiuvYP369dIJkXLal2vXrkVgYCBcXV2ltsq8/4pSln1WGfdrTk4OBgwYgPz8fKxYsUJl2vDhw6X/9vb2hqenJ3x8fBAfH49mzZq97FI1UtbvZGXch9988w0GDRoEMzMzlfbKsv+K+l0A9OfvkIexKkDVqlVhaGiolkxTUlLUUm5l8sEHH2DXrl04fPgwqlevXmxfFxcXuLu749q1ay+pOu2xtLREw4YNce3aNemqLLnsy1u3buHAgQMYNmxYsf0q8/4rzT5TKpXIzs7Gw4cPi+xTGeTk5ODNN99EQkICYmJiVEZ1CtOsWTMYGxtXyv364ndSLvvw6NGjuHLlSol/k4B+7r+ifhf07e+QYacCmJiYoHnz5mpDjTExMfD19dVRVWUnhMDYsWOxbds2HDp0CB4eHiV+5sGDB7h9+zZcXFxeQoXalZWVhUuXLsHFxUUaQn5+X2ZnZyM2NrZS7svIyEg4OTmhe/fuxfarzPuvNPusefPmMDY2Vulz7949/Pnnn5VmvxYEnWvXruHAgQNwcHAo8TMXLlxATk5OpdyvL34n5bAPgWcjrc2bN0fjxo1L7KtP+6+k3wW9+zvU6unOJNm8ebMwNjYWa9euFRcvXhQTJkwQlpaWIjExUdelaWzUqFHC1tZWHDlyRNy7d096ZWRkCCGESE9PFx9++KE4fvy4SEhIEIcPHxatW7cW1apVE2lpaTquvmQffvihOHLkiLh586Y4ceKE6NGjh7C2tpb21aeffipsbW3Ftm3bxPnz58Vbb70lXFxcKsW6PS8vL0/UqFFDTJkyRaW9Mu6/9PR0cebMGXHmzBkBQCxcuFCcOXNGuhKpNPts5MiRonr16uLAgQMiPj5edOjQQTRu3Fjk5ubqarVUFLeOOTk5olevXqJ69eri7NmzKn+XWVlZQgghrl+/LsLCwkRcXJxISEgQu3fvFvXq1RNNmzbVi3Usbv1K+53U531Y0ndUCCFSU1OFhYWFWLlypdrn9X3/lfS7IIR+/R0y7FSg5cuXC3d3d2FiYiKaNWumcql2ZQKg0FdkZKQQQoiMjAzRpUsX4ejoKIyNjUWNGjXEkCFDRFJSkm4LL6X+/fsLFxcXYWxsLFxdXUVQUJC4cOGCND0/P1/Mnj1bKJVKYWpqKvz8/MT58+d1WHHZ7Nu3TwAQV65cUWmvjPvv8OHDhX4nhwwZIoQo3T57+vSpGDt2rLC3txfm5uaiR48eerXOxa1jQkJCkX+Xhw8fFkIIkZSUJPz8/IS9vb0wMTERr7zyihg3bpx48OCBblfs/xW3fqX9TurzPizpOyqEEF999ZUwNzcXjx49Uvu8vu+/kn4XhNCvv0PF/xdNREREJEs8Z4eIiIhkjWGHiIiIZI1hh4iIiGSNYYeIiIhkjWGHiIiIZI1hh4iIiGSNYYeIiIhkjWGH6D8iNDQUzs7OUCgU2LFjh67LeanWrVsHOzu7Evv9F7cN0X8Bww7Rf8ClS5cQFhaGr776Cvfu3VN50vt/Qf/+/XH16lXpfWhoKJo0aaK7gnQgODgYffr0UXmvUCigUChgbGwMZ2dndO7cGd988w3y8/N1VyhRBWDYIfoPuHHjBgCgd+/eUCqVMDU1VeuTnZ39sst6aczNzeHk5KTrMvTOa6+9hnv37iExMRF79+5F+/btMX78ePTo0QO5ubm6Lo9Iaxh2iCqZgIAAjBs3DiEhIbC3t4dSqURoaGiR/UNDQ9GzZ08AgIGBARQKBYD//Us/PDwcrq6uqFOnDgBgw4YN8PHxgbW1NZRKJQYOHIiUlBRpfkeOHIFCocC+ffvQtGlTmJubo0OHDkhJScHevXvh5eUFGxsbvPXWW8jIyJA+J4RAREQEatWqBXNzczRu3BhbtmyRpj98+BCDBg2Co6MjzM3N4enpicjIyELX6aeffoKdnZ00AnH27FkoFApMnjxZ6jNixAi89dZbAFQPY61btw5hYWE4d+6cNLKxbt066XP//PMPXn/9dVhYWMDT0xO7du0qdn+sWLECnp6eMDMzg7OzM9544w1pWnR0NNq2bQs7Ozs4ODigR48eUvAEgMTERCgUCvzwww9o164dzM3N0aJFC1y9ehVxcXHw8fGBlZUVXnvtNdy/f19luZGRkfDy8oKZmRnq1auHFStWFFtnYUxNTaFUKlGtWjU0a9YM06dPx86dO7F3716VbUJU6Wn9aVtEVKH8/f2FjY2NCA0NFVevXhXr168XCoVC7N+/v9D+6enpIjIyUgCQnkwshBBDhgwRVlZW4p133hF//vmn9IC+tWvXij179ogbN26I3377TbRq1UoEBgZK8yt4wGGrVq3EsWPHRHx8vKhdu7bw9/cXXbp0EfHx8eKXX34RDg4O4tNPP5U+N336dFGvXj0RHR0tbty4ISIjI4Wpqak4cuSIEEKIMWPGiCZNmkhPeY6JiRG7du0qdJ0ePXokDAwMxKlTp4QQQixevFhUrVpVtGjRQupTp04d6WnSkZGRwtbWVgjx7MGnH374oWjQoIHak5oBiOrVq4uoqChx7do1MW7cOGFlZVXkwxfj4uKEoaGhiIqKEomJiSI+Pl4sWbJEmr5lyxaxdetWcfXqVXHmzBnRs2dP0bBhQ5GXlyeEENIDPQu2y8WLF0WrVq1Es2bNREBAgMr2HTlypDTf1atXCxcXF7F161Zx8+ZNsXXrVmFvby/WrVtXaJ0F+7t3795Fvn9e48aNVfY5UWXHsENUyfj7+4u2bduqtLVo0UJMmTKlyM9s375dvPhvmyFDhghnZ2eRlZVV7PJOnjwpAIj09HQhxP/CzoEDB6Q+4eHhAoC4ceOG1DZixAjRtWtXIYQQjx8/FmZmZuL48eMq8x46dKh46623hBBC9OzZU7z77rvF1vK8Zs2aic8//1wIIUSfPn3EvHnzhImJiUhLSxP37t0TAMSlS5eEEKphRwghZs+eLRo3bqw2TwDi448/lt4/fvxYKBQKsXfv3kJr2Lp1q7CxsRFpaWmlqjklJUUAkIJlQdj5+uuvpT6bNm0SAMTBgweltvDwcFG3bl3pvZubm4iKilKZ9yeffCJat25d5LI1CTv9+/cXXl5epVonosqAh7GIKqFGjRqpvHdxcVE51FRaDRs2hImJiUrbmTNn0Lt3b7i7u8Pa2hoBAQEAgKSkpCJrcHZ2hoWFBWrVqqXSVlDTxYsXkZmZic6dO8PKykp6ffvtt9JhnVGjRmHz5s1o0qQJQkJCcPz48WJrDwgIwJEjRyCEwNGjR9G7d294e3vj2LFjOHz4MJydnVGvXj2Nt8nz62VpaQlra+sit23nzp3h7u6OWrVq4Z133sHGjRtVDt3duHEDAwcORK1atWBjYwMPDw8AJW9L4Nm+eb6toIb79+/j9u3bGDp0qMq2nDt3rsohsvIQQkiHO4nkwEjXBRCR5oyNjVXeKxSKMl1BY2lpqfL+yZMn6NKlC7p06YINGzbA0dERSUlJ6Nq1q9oJzM/XUHBFT1E1Ffzv7t27Ua1aNZV+BSdLBwYG4tatW9i9ezcOHDiAjh07YsyYMfj8888LrT0gIABr167FuXPnYGBggPr168Pf3x+xsbF4+PAh/P39Nd4eL67Xi+vxImtra8THx+PIkSPYv38/Zs2ahdDQUMTFxcHOzg49e/aEm5sb1qxZA1dXV+Tn58Pb27vEbVlY24vbcs2aNWjZsqXKfAwNDcu0zi+6dOmSFMyI5IBhh4gkly9fxj///INPP/0Ubm5uAIBTp06Ve77169eHqakpkpKSig0hjo6OCA4ORnBwMNq1a4fJkycXGXb8/PyQnp6OxYsXw9/fHwqFAv7+/ggPD8fDhw8xfvz4IpdjYmKCvLy8cq8XABgZGaFTp07o1KkTZs+eDTs7Oxw6dAj+/v64dOkSvvrqK7Rr1w4AcOzYsXIvz9nZGdWqVcPNmzcxaNCgcs/vRYcOHcL58+cxceJErc+bSFcYdohIUqNGDZiYmGDp0qUYOXIk/vzzT3zyySflnq+1tTU++ugjTJw4Efn5+Wjbti3S0tJw/PhxWFlZYciQIZg1axaaN2+OBg0aICsrCz///DO8vLyKnKetrS2aNGmCDRs2YMmSJQCeBaB+/fohJydHOvxWmJo1ayIhIQFnz55F9erVYW1tXejl+CX5+eefcfPmTfj5+aFKlSrYs2cP8vPzUbduXVSpUgUODg5YvXo1XFxckJSUhKlTp2q8jMKEhoZi3LhxsLGxQWBgILKysnDq1Ck8fPgQkyZNKvV8srKykJycjLy8PPz999+Ijo5GeHg4evTogcGDB2ulViJ9wHN2iEji6OiIdevW4ccff0T9+vXx6aefFjmyoqlPPvkEs2bNQnh4OLy8vNC1a1f89NNP0uESExMTTJs2DY0aNYKfnx8MDQ2xefPmYufZvn175OXlScGmSpUqqF+/PhwdHYsNSn379sVrr72G9u3bw9HREZs2bSrTOtnZ2WHbtm3o0KEDvLy8sGrVKmzatAkNGjSAgYEBNm/ejNOnT8Pb2xsTJ07EZ599VqblvGjYsGH4+uuvsW7dOjRs2BD+/v5Yt26dxoeeoqOj4eLigpo1a+K1117D4cOH8eWXX2Lnzp1aOyRGpA8UQgih6yKIiIiIKgpHdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNb+D6h2K1/Yq7o6AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(1, 1)\n", + "ax.hist(\n", + " non_nan_frames_per_ID.values(),\n", + " bins=np.arange(0, len(ds.time) + 50, 50),\n", + " alpha=0.5,\n", + " label=\"Prediction\",\n", + ")\n", + "ax.set_xlabel(\"n frames with same ID\")\n", + "ax.set_ylabel(\"n trajectories\")\n", + "ax.hlines(\n", + " y=len(ds.individuals),\n", + " xmin=0,\n", + " xmax=len(ds.time),\n", + " color=\"red\",\n", + " label=\"n individuals\",\n", + ")\n", + "ax.legend()\n", + "ax.set_title(Path(ds.source_file).stem)\n", + "\n", + "# Save plot as png\n", + "plt.savefig(\n", + " output_figures_dir / f\"{Path(ds.source_file).stem}_histogram.png\",\n", + " dpi=300,\n", + " bbox_inches=\"tight\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "movement-env", + "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.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/notebook_trajectories_movement.py b/notebooks/notebook_movement_trajectories_gt.py similarity index 97% rename from notebooks/notebook_trajectories_movement.py rename to notebooks/notebook_movement_trajectories_gt.py index 3acdebc3..8883ea81 100644 --- a/notebooks/notebook_trajectories_movement.py +++ b/notebooks/notebook_movement_trajectories_gt.py @@ -1,4 +1,4 @@ -"""Inspect crab trajectories using movement""" +"""Inspect crab groundtruth clip using movement""" # %% import itertools @@ -22,8 +22,8 @@ # load ground truth (corrected) groundtruth_csv_corrected = ( - "/Users/sofia/arc/project_Zoo_crabs/escape_clips/" - "04.09.2023-04-Right_RE_test_corrected_ST_csv_SM_corrected.csv" + "/Users/sofia/arc/project_Zoo_crabs/tracking_groundtruth_generation/" + "04.09.2023-04-Right_RE_test_corrected_ST_SM_20241029_113207.csv" ) diff --git a/scripts/escape_trajectory_plots.py b/scripts/escape_trajectory_plots.py new file mode 100644 index 00000000..a0f1dfc6 --- /dev/null +++ b/scripts/escape_trajectory_plots.py @@ -0,0 +1,138 @@ +"""Generate trajectory plots for escape clips.""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from movement.io import load_bboxes + + +def main(input_data, output_figures_dir): + """Read input files as movement datasets and generate plots.""" + # Create a directory if it doesnt exist + if not output_figures_dir.exists(): + output_figures_dir.mkdir(parents=True) + + # List all csv files in the input directory + list_csv_files = [ + x + for x in input_data.iterdir() + if x.is_file() and x.name.endswith("_tracks.csv") + ] + list_csv_files.sort() + print(len(list_csv_files)) + + # Prepare plots + # select whether to plot ID at first frame + flag_plot_id = False + + # Define colors - ideally more than max n individuals + # so that we don't have repetitions + list_colors = ( + plt.get_cmap("Pastel1").colors # 9 colors + + plt.get_cmap("Pastel2").colors # 8 colors + + plt.get_cmap("Paired").colors # 12 colors + + plt.get_cmap("Accent").colors # 8 colors + + plt.get_cmap("Dark2").colors # 8 colors + + plt.get_cmap("Set1").colors # 9 colors + + plt.get_cmap("Set3").colors # 12 colors + + plt.get_cmap("tab20b").colors # 10 colors + + plt.get_cmap("tab20c").colors # 20 colors + ) # 96 colors + + # loop thru escape clip files + for csv_file in list_csv_files: + # Create movement ds + ds = load_bboxes.from_via_tracks_file( + csv_file, fps=None, use_frame_numbers_from_file=False + ) + + # Print summary metrics + print(Path(ds.source_file).name) + print(f"Number of frames: {ds.sizes['time']}") + print(f"Number of individuals: {ds.sizes['individuals']}") + print(ds) + print("--------------------") + + # Compute number of frames with non-nan position per ID + non_nan_frames_per_ID = {} + for ind, _id_str in enumerate(ds.individuals): + non_nan_frames_per_ID[ind] = ( + len(ds.time) + - ds.position[:, ind, :].isnull().any(axis=1).sum().item() + ) + + # Plot trajectories per individual + fig, ax = plt.subplots(1, 1) + for ind_idx in range(ds.sizes["individuals"]): + # plot trajectories + ax.scatter( + x=ds.position[:, ind_idx, 0], # nframes, nindividuals, x + y=ds.position[:, ind_idx, 1], + s=1, + color=list_colors[ind_idx % len(list_colors)], + ) + # add ID at first frame with non-nan x-coord + if flag_plot_id: + start_frame = ds.time[~ds.position.isnull()[:, ind_idx, 0]][ + 0 + ].item() + ax.text( + x=ds.position[start_frame, ind_idx, 0], + y=ds.position[start_frame, ind_idx, 1], + s=ds.individuals[ind_idx].item().split("_")[1], + fontsize=8, + color=list_colors[ind_idx % len(list_colors)], + ) + + ax.set_aspect("equal") + ax.set_xlim(-150, 4200) # frame size: 4096x2160 + ax.set_ylim(-150, 2250) # frame size: 4096x2160 + ax.set_xlabel("x (pixels)") + ax.set_ylabel("y (pixels)") + ax.set_title(Path(ds.source_file).stem) + ax.invert_yaxis() + + # Save plot as png + plt.savefig( + output_figures_dir / f"{Path(ds.source_file).stem}_tracks.png", + dpi=300, + bbox_inches="tight", + ) + + # Plot histogram of trajectories' lengths + fig, ax = plt.subplots(1, 1) + ax.hist( + non_nan_frames_per_ID.values(), + bins=np.arange(0, len(ds.time) + 50, 50), + alpha=0.5, + label="Prediction", + ) + ax.set_xlabel("n frames with same ID") + ax.set_ylabel("n trajectories") + ax.hlines( + y=len(ds.individuals), + xmin=0, + xmax=len(ds.time), + color="red", + label="n individuals", + ) + ax.legend() + ax.set_title(Path(ds.source_file).stem) + + # Save plot as png + plt.savefig( + output_figures_dir / f"{Path(ds.source_file).stem}_histogram.png", + dpi=300, + bbox_inches="tight", + ) + + +if __name__ == "__main__": + input_data = Path( + "/home/sminano/swc/project_crabs/escape_clips_tracking_output_slurm_5699097" + ) + + output_figures_dir = input_data / "figures" + + main(input_data, output_figures_dir)