{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise session, week 11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VgIIA5l2fkWd"
   },
   "source": [
    "**Course**: [_Systèmes dynamiques en biologie_](https://moodle.epfl.ch/course/info.php?id=14291) (BIO-341)\n",
    "\n",
    "**Professor**: _Felix Naef_ , _Julian Shillcock_\n",
    "\n",
    "SSV, BA5, 2025"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "executionInfo": {
     "elapsed": 382,
     "status": "ok",
     "timestamp": 1638526734073,
     "user": {
      "displayName": "Lorenzo Talamanca",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GhUF0b2-P5yq7rSXU0BUszliuOAWwm2O82n9owJyh4=s64",
      "userId": "06843078895857810674"
     },
     "user_tz": -60
    },
    "id": "xv0WBenAfkWl"
   },
   "outputs": [],
   "source": [
    "# import important libraries\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import set_matplotlib_formats\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# set_matplotlib_formats(\"png\", \"pdf\")\n",
    "matplotlib.rc(\"image\", cmap=\"RdBu\")\n",
    "import matplotlib.animation as animation\n",
    "from random import random\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Lz2eBiVHfkWl"
   },
   "source": [
    "## A. Organization of the circadian oscillator network in the brain\n",
    "  \n",
    "In mammals, circadian timing is regulated by the master clock, which is composed of about 20000 neurons located in the suprachiasmatic nucleus (SCN) in the hypothalamus. The neurons in the SCN are synchronized with each other, in part, by neurotransmitters that are secreted by the neurons. Here, we are interested in how different neuron connectivities affect collective synchronization and oscillatory behavior. For simplicity, we will consider a theoretical model of the SCN in which coupled oscillators are arranged in a 2-dimensional space (grid). We can construct a symmetric adjacency matrix (coupling matrix) $\\mathbf{A}$ where $\\mathbf{A_{ij}}=\\mathbf{A_{ji}}$ specifies the coupling between oscillator $i$ and oscillator $j$. Here, we consider the simplest case where $\\mathbf{A_{ij}}$ is either 0 or 1. This coupling is modulated by a strength $K$. \n",
    "\n",
    "For each neuron $i$, $\\dot{\\theta}_i = f_i + K  \\sum_{j=1}^{N}\\mathbf{A_{ij}} \\sin (\\theta_j-\\theta_i)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1GERYChgfkWm"
   },
   "source": [
    "We provide you with two functions: \n",
    "- heatmap: plot a $n \\times n$ grid ($N := n^2$ being the number of neurons), which represents the SCN. Each cell shows the value contained in the matrix _datagrid_. The values are plotted from bottom left to top right, row by row.\n",
    "- value_to_size: create the grid on which the oscillators are. This function is called by heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "executionInfo": {
     "elapsed": 3,
     "status": "ok",
     "timestamp": 1638526734454,
     "user": {
      "displayName": "Lorenzo Talamanca",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GhUF0b2-P5yq7rSXU0BUszliuOAWwm2O82n9owJyh4=s64",
      "userId": "06843078895857810674"
     },
     "user_tz": -60
    },
    "id": "0Pnn7iJ1fkWm"
   },
   "outputs": [],
   "source": [
    "# create a grid on which the oscillators are\n",
    "\n",
    "def value_to_size(val, size_min, size_max, size_scale):\n",
    "    val = abs(val)  # to handle negative  numbers\n",
    "    val_position = (val - size_min) * 0.99 / (\n",
    "        size_max - size_min\n",
    "    ) + 0.01  # position of value in the input range, relative to the length of the input range\n",
    "    return val_position * size_scale\n",
    "\n",
    "\n",
    "def heatmap(datagrid, marker, ax):\n",
    "    x = list(range(datagrid.shape[1]))\n",
    "    y = list(range(datagrid.shape[0]))\n",
    "\n",
    "    size_min, size_max = 0, 1\n",
    "    size_scale = 500\n",
    "\n",
    "    # dot the scatter + 2 extra invisible points to normalize colors\n",
    "    scat = ax.scatter(\n",
    "        x=x * len(y) + [-100, -100],\n",
    "        y=[v for v in y for p in x] + [-100, -100],\n",
    "        marker=marker,\n",
    "        s=[value_to_size(v, size_min, size_max, size_scale) for v in datagrid.flatten()]\n",
    "        + [1, 1],\n",
    "        c=[v for v in datagrid.flatten()] + [-1, 1],\n",
    "        cmap=\"RdBu\",\n",
    "    )\n",
    "\n",
    "    ax.set_xticks(x)\n",
    "    ax.set_xticklabels(x)\n",
    "    ax.set_yticks(y)\n",
    "    ax.set_yticklabels(y)\n",
    "\n",
    "    ax.grid(False, \"major\")\n",
    "    ax.grid(True, \"minor\")\n",
    "\n",
    "    ax.set_xticks([t + 0.5 for t in ax.get_xticks()], minor=True)\n",
    "    ax.set_yticks([t + 0.5 for t in ax.get_yticks()], minor=True)\n",
    "\n",
    "    ax.set_xlim([-0.5, max(x) + 0.5])\n",
    "    ax.set_ylim([-0.5, max(y) + 0.5])\n",
    "    ax.set_facecolor(\"#F1F1F1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ywcR40TifkWm"
   },
   "source": [
    "Below is an example of how to use the functions to plot the phase of each neuron in each cell of the grid: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 592
    },
    "executionInfo": {
     "elapsed": 849,
     "status": "ok",
     "timestamp": 1638526735301,
     "user": {
      "displayName": "Lorenzo Talamanca",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GhUF0b2-P5yq7rSXU0BUszliuOAWwm2O82n9owJyh4=s64",
      "userId": "06843078895857810674"
     },
     "user_tz": -60
    },
    "id": "GRpz-5RCfkWm",
    "outputId": "ae715686-1266-4d07-9528-18efc028638c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100,)\n",
      "(10, 10)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy0AAAMtCAYAAACSEJ4GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAC0jklEQVR4nOzdd5hU5d3/8c85Z7Y3lt57ExRUEEXsNfbeMNaYaGJ6TDHleUxMYn5PEpMYicZeEDXW2LtgQ6QIIr33vpWtM+fcvz9mF0FZdmd2Z+ce5v26rrnUnTPuly+n3J/7lHGMMUYAAAAAYCk32QUAAAAAwL4QWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArBZq718YBIE2btyogoICOY7T3r8eAAAAgCWMMaqsrFTPnj3luk2fT2m30DJp0iRNmjRJ9fX1WrFiRXv9WgAAAACWW7dunXr37t3k+44xxrRjPSovL1eHDh302WefqaCgoD1/dcryfV+ffPKJxo0bJ8/zkl1OyqBvsaNn8aFvsaNn8aFvsaNn8aFvsaNn8SkrK9MhhxyisrIyFRUVNblcu18e1nhJWHFxsQoLC9v716ekSCSi3NxcFRcXKxRq97+ylEXfYkfP4kPfYkfP4kPfYkfP4kPfYkfPWqe520a4ER8AAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYLZTsApLN+BHVL56n8OqlMvV1cvMKlDlyjDJ69U92aVYzxijYul6molRyXblde8st6JDsslKCMUYKfDle2m9+AAAALZLWo6baeR+r+s1nZaqrJNeVjCRHqv34HXndeiv/3CsV6tYr2WVaJ7LkU9XPflemdNtuP3Xk9R+mjHEny+vSM2m12coYI21ZLbN0trRppWQCmYwsacBBcoaMkVNQnOwS7VWzU+6GxVJVNCCb4l4yPQZJXkayK7ObCeTsLJXjh2Wy8mRyCpJdEQAAcUvb0FIz411Vv/70Fz8Igug/TfQf/taNKn/wryq6+scKde/d/gVaqn76awrPmbaXd4z8NUvkr1um7DOultdncLvXZitjjMyct6RlsyXHkUzDShauk5bNlln+qTThXDm9hiS3UNsYI3f5J3JWz5PkqHFWwdm8Qlo6XcFBJ8p06ZvkIi1kjNxtqxTauEROpG7Xj4O8YkX6HCiT3ymJxdnNra1UqGyD3JpKyXXk53dWpKiHTCgr2aXZzRh54Wo5QURBKEtBKDvZFQHYD8V8T0tlZaV++MMfql+/fsrJydGRRx6pmTNnJqK2hIls3bhnYNkbE0jhsCr/c4+MCdqnMMtFls9vIrA0MEbyA9W+8qhM9c72K8x2Sz6JBhbpi8DSqOFSMfPBczJlW9u/Nou5y2fKWT1PjiRHZo9/yg/Lnfe6VLopuUVayNu4SBlrP5N2CyyS5FSVKmPJB3IqtyepMruFyjYqe91ceZXb5EZq5dbXKFSyTtlr5sipY3/WFK++Srkly5VTvk7ZlZuUW7pa2aVr5ASRZJdmNxMoVL9TWdXblVVTIi9S89XjA5pGr9JSzKHluuuu05tvvqlHH31U8+fP1ymnnKKTTjpJGzZsSER9CVE7c1r0crDmmEBB2Q6FVyxOfFEpIPzptOiZgn0ykh9WeNGsdqnJdsaPyCyc3rJlF3+S4GpSSF21nDXRwLI3jiQZyVtGz3bn1FYqtGlp9N+//J4kGaPQ6k854H+JU1eljG0rov+++88lKYgoa+NCerYXjl+v7PL1cgJ/j597kRpll6+jZ00xgbJrdiijvlJuEJbr1ymrtkwZ9ZXJrsx6rvGV5VcrK6hRpl8bnWBG2ogptNTU1OiZZ57R//3f/+mYY47R4MGDdcstt2jAgAG66667ElVjmzImUN28GV9cDtYcx1XdvI8TW1QKCEq3Kti6oWUHIWMUWchgUpK0cblUX9v8ciaQ1iyUCdc1v2wacDYt23WpZpPLyMgp3yJVlbVLTanA3ba68VzUXjmS3LoqOTt3tF9RKSBUvrHJ9xxJbqRObnVp+xWUIjJqSqW9rHGOJC9SJy9cnYSq7BcK18gJIg1nj78IyqFwFWeo9sExgTKC6DEy2rdAmUEd4biFPM+T25IJe4vFVH0kEpHv+8rO3vN61ZycHH3wwQd7/UxdXZ0qKir2eCWTqauVIuEYPhAoKOdgFVSWxbS8qUru37M1qspbcHaqgQmk2qrE1pMinOryr54qaHJZ1rVGbk2lnGbSnpHk1DCjuzu3duc+VzejaNjDnrxwTZN9M5LcSAsmbNKQ6+99csqR5Pr17VtMCnFN9Ixe4zrnSHKbm92CJMlzXXmuq5DnyWnpmMRCMYWWgoICjR8/Xrfeeqs2btwo3/c1efJkzZgxQ5s27f3a8ttuu01FRUW7Xn369GmTwuMV12NmQ2n7vIJdnFCMT2pyvcQUkmpcL7ZZIPoWFUsfUnzmqC0Z12v2EO5IrGdf1oK+sZ59lXGaGTI69Gyv9tUXetYkI+crIZnI0jKmYRxijNn176ko5q3j0UcflTFGvXr1UlZWlu644w5NnDhRnrf3g+DNN9+s8vLyXa9169a1uujWcDIy5XXurpZP4zrK6DMooTWlArdzD6mlgc9x5fYckNiCUkXXfi1fNq9Iyi1MXC0pxHTqI6cFO1bjhWSKurVDRakh6NC92WWMHAVFXduhmtTh53dufpk8nrr2ZZHsfe+vIlk8ZntvIhm5ex18G8eV7/GkuqYEjqeg4Vxy49HBd0Itv5ohjQXGKBwOKxxJ7csPYw4tgwYN0rRp07Rz506tW7dOn3zyicLhsAYM2PsgNSsrS4WFhXu8ki173HGKJZ9njZmQsFpShZOZrdDwQ1s2C2QCZRw0PvFFpQCnQxepc68W7VSdIYem9GnbtmQ695bJzt/n/RlGjkyv4VKsZwH3Y0HH3lIoq8m9m5EUdOotZfBI2t1FCrvJNNE3IylS2F2Gnn1FJKtIfkbuHoPIxn/W53WRcblKYW8CL1P1mYV7rG/GcVWb3ZEB+L44jurdbEWcDPlOSPVupiJuZrKrShmpe37lC3Gfh8zLy1OPHj1UWlqq119/Xeecc05b1pVQWaPGyS3u3IIBuKOsMUfJK+SL/yQp49DjpMysZk5tO3L7DJHXl+8caeSM/Vr0LFVTByPHkTr1kAYf2r6F2cxx5Y8+RfI8mb30zciRCjopGHRYEoqzmOspPPRIKZS514Gkye+sSN/RSSrOYq6nut6jFWRHzwzs3rdIUQ+Fu3K2fa8cR7VFvaMBxcuMninIyFVNYW+FczkztS+RzDzV5HVTXXaxarM7qja3qwxfmNs8x5HvZijiZipwCMXpJua/8ddff13GGA0bNkzLly/XT3/6Uw0bNkzXXHNNIupLCCczS4Vf/54qHrlDQUXJV+85aPgCwMwRByvv1IuSU6SF3MJi5Zz3LdW+8IBMdeWeX5TouJIJ5PYZouxTJ8rhutxdnA5dpBO/LjP9Balix26hz0T713uYnHGnxX7f0P6usLP8w8+Xu2KWtGXVrhvMTShTpvcIBQMPkTjIf4XJLVL9gSfJ275GbskGGb9eJitfQZf+0cvH2Db3ymRkqa7PwXJqd8qtq4wG59xiKcRM7j45rsK5nQgp8XBc+XwRJ9BiMYeW8vJy3XzzzVq/fr06duyoCy64QH/4wx+UkZFagwevuLOKvvUL1c56T7Uzp8ns/OIJRKFe/ZU97jhljjyUwfeXuJ26K+eKnyqy/DNFFnyioHyHHMeV27O/Mg4aL7dHfy5x2gunuJt02nXStvUyG5dLflhOdr7Ub4Sc/A7JLs9eeR0UjDpJqq+RqiuiA+6CjtxI3pxQpvzuQ+R354xnrEx2vvzs/GSXAQD4kphDy8UXX6yLL744EbW0OzcnV7lHf005E05WULpdJlwvJzefy8Ga4YQylDF8jDKGj0l2KSnFcRypax85XZP7BL2UlJkTfQEAgLTEBYGSHNeT14knEAEAAAA24tonAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYLVQsn5xJBJRJBJJ1q9PKY19ol+xoW+xo2fxoW+xo2fxoW+xo2fxoW+xo2fx8X2/Rcs5xhiT4FokSZMmTdKkSZPk+76WLl2qKVOmKDc3tz1+NQAAAAALVVdXa+LEiSovL1dhYWGTy7VbaGlUUVGhoqIirVq1ap+F4QuRSETTp0/X+PHjFQol7eRYyqFvsaNn8aFvsaNn8aFvsaNn8aFvsaNn8SktLdXgwYObDS1J62goFOIvNEb0LD70LXb0LD70LXb0LD70LXb0LD70LXb0LDae57VoOW7EBwAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFgtlOwCks0Yo02z56t0xRpFamqVWZCvXuNGq7BPz2SXZjW/PqytcxeotrRcbshTQe8e6jhsULLLsl64ukZVGzbLr6lVZmGB8vr0kOt5yS7LeiYIFNlZLSfkycvJluM4yS4JAAC0o7QNLeHqGi186mV9eu/jKlm2as83HUcDTzpKB193qfodewQDpN3s3LhFCyc/q88fflq1JWV7vNdp5FCN+salGnLuqQrlZCenQEuVLlqmJY88rRVPvSy/rm7Xz3N7dNXwqy/W4EvOVnan4iRWaB9jjLbP+FQrH31KG994TyYSkSRlFBVqwGXnasDE85TH5MJXGN/X9vc+1rrHnlHFgqUK6uoUKshXlxMmqM/E85U/ZGCyS7RSEI6oZOoH2v7Guwpv2yEn5CmnXx91u+As5Y8YluzyrBVEIip7/yNVL16ioK5OXkGBOhx5uPIOGJ7s0qxmIhFVLfhckZJSOSFPmb17K2cA22ZzTBDILytVUF8vLy9fXkFBsktCO0rL0LJz8zY9c/F3tGPJCmlvgcQYrXrnI618832NuvICnXDbz+WG0rJVe1j1+jS98a1fKIhEZILgK+/vWLRc7/74d/r0zod01hP/UkGfHkmo0i6B72vmLX/VkoeekuN5Mr6/x/vVm7Zqzp8mae5f79Exd/5efU87PkmV2qV22w5Nv+4mlX628Ct9C5dXaNk9k7X0349o0NWXatSvvi+Hs1WSpNKZczX/J7eodtOWPfoWqdyp9Y8/r3WTn1HHIw/TqNt/q0xC8i6bn35Ba++8T+GSUsl1pYb9W/nsedr81H+Vd8BQDf7fnxFedhNEIto8+Qltfeo5RcrL5YQ8SU706oUHH1XOkEHqee2VKj56QrJLtYq/c6d2vPiCdrz8svzy8j3eyx40SJ3OPkcdjjtOjsvV+7vzKytUMfVtlb/1qvzSkl0/zx4yTEUnn668sYfLYZz2FUFtjWpmfajqT6ZF++a6yujdX3kTTlLW8INSbj1LrWrbQE1puZ485xsqWb46+gNj9rpc48H+s0ef1Vs//aNME8uli9VvvqdXr/mJ/HB4r4FF0hcH+jUb9MzZ16hqy7Z2rNA+xhh99JPfacnDT0f/+0uBZbcFFYTDmnr9z7X6pbfasUI71W7boXfPvUZlC5ZI2nvfTBBIRlrx4BOa+eNbml4n08j29z/WrCu/q9qG7e7LfWv879IZc/TxBd9Q3faSr/w/0tGaf96rFb/7czSwSLv2Y5Kkhp5VLVmuz678jso+mZOECu0T1NVp2U2/1Mb7H1akYeBtIn70bGhDz2qWr9SKm/9Xm6f8J5mlWqV+61at+NGPtPWJJ74SWCSpduVKbfjb7Vp/+1+bPl6kodrlS7X2p99TydNT9ggsje9t+dfftOGP/yN/Z2WSKrRTzeezteWW76v82UcUXr9GQVWlgspy1S2Zr5L7/qptf/6lIiWpNU5Lu9Dyzs1/UsW6TS3fIRijz6c8r6UvvJnYwixWs71Ur3/rF9H/aEF4M76vmm0leuu7v0lwZXZb9thzWvnMKy3qWeMyH3z/N6pcuyHBldltxnduVu2W7S3eRte/8IaWP/hkgquyW836TZr7nV/I+MGeg+69ML6vus1bNfeGn6X9ZMyW51/W+nsfaX7BIJCJhLXoez9X7fqNiS/Mcqv++GdVzvl03/u2hvfW/+selbz1bjtVZi+/qkqrf/0r1W/d2vQ22tCz8vfe08a7727H6uxVv36dNv6/3yqoqdn7+tbws7qVy7XpL39QUF/fzhXaqXbBpyp98A6Z+rqvvtmw/kW2bdb2O26VX17aztXFL6bQEolE9Otf/1oDBgxQTk6OBg4cqN/97ncKUmSWs2rrdi194c2YZzAc19Wce6YkqCr7LXriv/Lr61s2+G5gfF8bPpipkqWrml94P2SM0YK7H9375YdNf0gmCLR08rOJK8xyJfMWaseseTFvo8vueVRBwz0v6Wjt5KdkwpEWb6PG91U+b4HKZn+W4MrsZYJAa//1QMs/EBgF9WFtfOzpxBWVAqpXrFTp21OloOXHg/V335f2Z0NLXn9N9Zs3S0EL9m3GqPS1V1W3fn3iC7Pc9scekAmHJdPM+hMEqlu5XJXvvdM+hVnM1Ner9LF/N388CAIFOytU/sLj7VNYG4gptPy///f/dPfdd+vOO+/UokWL9H//93/685//rH/+85+Jqq9NzZ/8nBTHxKIJAm2a9Zm2LVja9kVZLvB9zb//yZgOUI0cz9OCh59KQFX22/zRLFWuWR9T0JMk4wda+tiz8mv3MjuSBlZOfiau+1Nqt+7Q5nc+TEBF9vNra7X+yf/GPhnjeVo3OX0H4GUffqL6zVtj+5Dva8tzL8mvrklMUSlg23MvSDFuo/Wbt6hi5uwEVWQ/EwQqefHF2I4HrquSV19JXFEpILxlk2oWzG/27PEXHJW/+Wran0GumTtDpra6ZQsHgWrnfSK/8quXK9ooptAyffp0nXPOOTrjjDPUv39/XXjhhTrllFM0a9asRNXXpla++X7csz2O62r1ux+1cUX2K1u+WlWxHtgbGN/Xqtemtm1BKWL9Wx803Jgau3DFTm379PM2rig1bHpjalzXcjuep01vv5+AiuxXNusz+TtbeIDajfF9bXljWtoe4He8817Mg29JCqprVD4zfe9tKXlr6q77VlrM81Q69b2E1JMKapYvV3j79tg+FAQqmzo1IfWkip0zPoo+GKPFjMKbNqh+w7qE1ZQKqmd9EPNVHjVzZySuoDYUU2g56qij9Pbbb2vp0ugZh3nz5umDDz7Q6aef3uRn6urqVFFRsccrWb78iN5YOJ6r2rL0u8mrtrR1f191FenXM0mqLyuP66zeF59P3naSLMYYhSur4vus76dlz6SGdS1OJhxWUJeeZ/XCZeUtu1RnLyJpuq4ZY+RXxbGN+r4irVhPU93ebrpv0ecqK9N2UkGS/PKy2Abfu38ujfnlpbGd1XNcBSlypiWm58P9/Oc/V3l5uYYPHy7P8+T7vv7whz/osssua/Izt912m37729+2utC24GZmxv9hI3mZ6fc4PS8ro3Wfz2jd51OVm5UZ1852j8+nGcdx5IS86L0ZsXIduem6rrXyz52ujwl1MzIkx23+Wvm9cDLTc12TolcdxHw21HHkpOn2KcW/jTmhUHp/T1wovnXGifNz+wvHi3V9M3F8JjliOtPy5JNPavLkyZoyZYrmzJmjhx9+WH/5y1/08MMPN/mZm2++WeXl5bte69Yl77Rd8YA+cX+fQxCJqDANv8iuoFf3+AffrqPCfr3atqAUUdCvd6tuPC3om559y+vbK671zXFc5ffrnYCK7Jc7oE/cn83u0S1tv4Mqp3/fpHw2lTmOo+x+fWPfRh1HOf37JaaoFJDVt09cPcvqm57rWaOsfv3juhQxs2d6Hj8bZfQfFNtldUGgjD4DEldQG4optPz0pz/VL37xC1166aU66KCDdMUVV+hHP/qRbrvttiY/k5WVpcLCwj1eyXLg5efG/ezzUE62hp51UhtXZL/crp3V78QJ8YW9wGjEFRe0fVEpYOD5pyueqOe4rjofcqCKBvdv65JSwsDLz4/rcybw1e/is9q4mtRQMHSQCg86IMZrvyW5jvpcnp7bpyR1O++MmB+UIddR3vAhyj9gaGKKSgFdzj879r7JqPOZpyWknlSQ0amzCg47LLZt1Bh1OvPMxBWVAvLGHiE3N7flH3Bd5Y8bL68geeNMG+QfeWIMDy+Q3KJiZQ0flcCK2k5MR7nq6mq5X9roPM9LmUceDzjpKOV17xLz5xzP08jLzlZmfl4CqrLfgVdfHFfYy8jL0ZBzT01ARfbL7dZZfU87PuawZ4JAw6+5JEFV2a/vBWfIzYht5t/xXHU9+vC0PdMiSX2vujimg5QUDci9LkzfQVFWj24qPvbI2G7GD4x6TEzfoCdJnU45UW5OTsvPHHiuOhx1pDK7xn7s3Z90Ovuclm+jrisvP19FRx2d2KIs52ZmquiUM1q+rhmjolPTd5/WKKN3f2UOPbDFIbng5HPkxDrplSQxVXnWWWfpD3/4g15++WWtXr1azz33nG6//Xadd955iaqvTbmep6N++d2YPuO4rkI5WRpz/dcTVJX9+h4/Xj2OOFSOF9tKPe6n31ZGbk6CqrLf6B9/S14M97Y4nqdOo0eo3xknJrgye2UWFujAn9/Y8g840XtZDvx5bNv1/qb76Sep4/ixUgzb6NCffVeZnYoTWJX9Bv78BwoV5Lesb66rDhMOV9cz03MippGXm6uB//vL6H80t2/zPGUUF6vvj7+f+MIslz96tLpc0oIJKdeVXFd9f/1rudnZiS/McsXnXKjcUYe06Dja+arrlD1wcDtUZb+OV96oULde0fv29iH36FOUO/74dqqq9WIahf7zn//UhRdeqO985zs64IADdNNNN+n666/Xrbfemqj62tzIS87SETd9q0XLOp4rNzND503+hzr0T99ZXMd1dfpDt6vj8MEtTuOjr79co741McGV2a3D0IE64aG/ycvKbDbwOZ6rosH9dOLDf5eXxjf5StKgay7VsO9eE/2PfRyoHM+Tl52l8fffrg4j0vdyHUlyM0I6eNKf1OHQUdGeNdW3hvVw0A++qX7XXNqOFdopu1cPHfTAP5XRsbjpnjXs84onjNPwv96atg8u2F2Ho8Zr0B9uifZib31r6FlWj+4aftc/lNm5UztXaKeul39d3a6+Onp278t9a+iZV1CoAb//g/JGHtj+BVrI8Tx1/8HPVHTqGVLj+tbYu8aedShWt+/8SEUnpPeEwu7c3Dx1/t6vlXfsqXKyGsKv4+7qnde5mzpc8g0VnXt5Sj3sIaa9b0FBgf7+97/r73//e4LKaR9H/vQG5XfrrPd+d4fqK3dGn4ay22lbx/NkfF8dBvTR6Xf9Ud1GHZDEau2QVVSg8/97v9775f/T0mdelUwgs9sXTjb2MKtDkcb99HodeM3FKbUhJEr38WN02n8f0Kzf/k2bP5q1a92SpOhNL47czAwNuvAMjf31D5SRppcg7s5xHI38yQ0qGj5Ei+98QBWLl+/2nTfOrm21x8nH6IAfflNFwwYlr1iLhAryNPahO7Tu8ee09pH/qGbtBsl15TiOTOBLRup4xBj1/8ZEdT76iGSXa43cwQN0yLOPaMszL2rT48+ofsu2Pd4vGH2gelx2vjqffFzcD3LZHxUfM0Gjnpmi7S+/qq3P/Ffh7Tt2vZd3wHB1vfBcFR97VOue2rmfcRxHXS64UMUnnazSN99U2TtvK1xSKifkKatPH3U6/QwVHHFE2j4JsSlOKKTOE69W8TkXqvKDaapbvkRBXZ28/ALljRmn3IPHsG3uhZudo6KzL1PB185X3YJPo49Cdj1l9OyjzEHDU3KMlrZTRqOuvFAHXHSmlr7wpuY9+JRKVqyWX1unjPxc9ZlwmA6+9hL1Hn9oSv6lJkpGXq5O/MdvNf7X39eix/+rlS+/o5odpXJDIRX07akRl52jgaefEL0kCrt0HDFUpzx5l8pXrNayKc+rbOlKRaqqldmhSD0mjNWgC89UZlFBssu0Tu8zTlSv009Q6byF2vj6VNWVRNe13N491Pe805TTvWuyS7SOm5mhflddrL5XXqTSGXNUuXi5gro6hfLz1XHCYcrrH/+TxvZnGUWF6n3t5ep11aXauWCxwiVlckIhZffpqZx+9KwpGR2L1eOKiep++aWKlFc0DCTzFMrPT3ZpVgsVFanLhReqy4UXJruUlOLl5avDqWdIp56R7FJSipuZpZxD9o+JqrQNLZKUkZOtkZecpZGXpOdTh+KV26WTxnz/Wo35/rXJLiWlFA3qr7G/+WGyy0gpjuOo48Ej1fHgkckuJaU4jqOOR4xRxyPGJLuUlOJ4ngpGsa7FynFdZRR3SHYZAPZzqfG4AAAAAABpi9ACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsFooWb84EokoEokk69enlMY+0a/Y0LfY0bP40LfY0bP40LfY0bP40LfY0bP4+L7fouUcY4xJcC2SpEmTJmnSpEnyfV9Lly7VlClTlJub2x6/GgAAAICFqqurNXHiRJWXl6uwsLDJ5dottDSqqKhQUVGRVq1atc/C8IVIJKLp06dr/PjxCoWSdnIs5dC32NGz+NC32NGz+NC32NGz+NC32NGz+JSWlmrw4MHNhpakdTQUCvEXGiN6Fh/6Fjt6Fh/6Fjt6Fh/6Fjt6Fh/6Fjt6FhvP81q0HDfiAwAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKwWSnYByRaprdOSV6dqx4o1ClfXKruoQH2OOFh9xh0sx3GSXZ6VjDFaP3u+lrzxvmpKy+WGQirq2U2jLviaCnt2S3Z51tq5vUSfPvGSdqxco/rqWmUXFqj/EQdrxJknKpSZkezyrBT4vpa8/ZGWvPG+qkvL5YY8FfXopkMvPVPdhg9KdnnWKtuwWbOnvKDtKxrXtXz1P+IQjT7/VGXm5iS7PCsFvq/Fb37QsK5VRNe1nt106CVnqPuIIckuz1qVW7Zr9uMvaPuKtaqvqVVOYb76HXGIRp1zskJZmckuz0rGGK2ePmfP/VrPbjr4wtNV3K9XssuzVn1VtRa+8Ja2N47XCvPV94hDNPDYwxmv7UPZ2g1a8daH0fFaRkiFPbpq6GnHKzM/N9mlxSxtQ0vl5m2a8e8pmvPw06otr5Qb8iTHkQkCGT9Ql2EDNe5bE3Xw5efIy2BAKUlBEGjuEy/qo39N1uYFS+V6nuRGdxTGD/T6LX/XAacfr6O/f7X6HDYqydXaY8ui5Xr39vs0/7nXZYJAjudJgZFcR9PvmaLcjh10+LUX65jvXaXsooJkl2uFSH1YH979mD64e4rKN2yWGwrJmECOouvbu3+7XwPGH6rjfvwNHXDK0Umu1h7rP12ot//8by18dZrkSI6caN9cVzMffU4v/PxPGnfVBTrhJ9cpr1Nxssu1QqQ+rA/+9ag+uHuKKjZtja5rQSA5jhxHmvr3B9Tv8IN1/I+/oRFfOzbZ5Vpj0+dL9fZf7tH8F96SjOS4jowxclxX0+//j/770z/qiGsu1rE/uEa5xYXJLtcKge9r5iPP6sO7Jmvb0lXRcUd0Q5UJAr3x+zs17KSjdOwPr9GACWOTXa41ytZv0vRJj2rO5OdUX1UtN6Nh6GqMgoiv4v69dfg3L9PYay5SRnZWcou1yMqp0/XJ3Y9pxdsfSpJcz5VMdD3MyM3R6MvP1WHfvEwdB/ZNcqUtl5aXh22ev1j/PvZiTZ/0sGrLKyVJQcRXEI7I+IEkadvSVXr5pt/rsYtvVF3FzmSWa4VIXb3+c90v9OyN/6Mti5ZJiq74QTgS7VsQyASBFr86Vfd87SrNevTZJFdshyVvvq87j79M8599XUHElwlMtGcNvZOk6pIyTb39Pv3r5K+rfMPmJFecfLWVVbr/gm/r5f/5265+BJHothn4vgLflySt/mSuHrzoRr3zl3uTWa41Pnv+DU06+eta9Np7uyZfAr9hnYtEe1a3s1of3j1Fdxx3mXasWp/kipOvtmKn7jv/Br16yz9UsWmrpIZ1LQhkfH9X39bO/EwPXfI9vf3ne5JZrjUWvjpVdxx/mT5/8W0ZP7rvDyJ+dJ1r2K/VlFVq6h0P6s4TL1fp2o1Jrjj5wrV1euzKH+v5H92qbctWS2oYd0QiX4w9jNGydz7SPWd+Q9PveyK5BVti/azPdNfRF+mT+55QfVW1JO0adzRun6VrNui1X/9FD5/7TVWXliezXCuYINA7v/uHHr/w21r57nTJmF0Br/H4Ga6u0ewH/qN7j71oV6hJBWkXWnYsX6OHz/6makrLdwWUvTJGMtKaD2fpict/IL8+3H5FWiYIAj3znd9owX/fkCSZwDS9rO/LBIGe//5v9dkzr7ZXiVZaPX2OHpn4A0Xq6nftKJpigkA7VqzVfed8UzVlFe1UoX38cFiPXvEjrfxwVnQb3IfG7fe1W/+pD+6a3B7lWWvp2x/psWt/Jn+3g1JTAt9X+YYt+veZ39DObTvaqUL7ROrDevjyH2rVR7NlmlvXgui69vrv79R7kx5tj/KsteL9mXrk8h/J323Q2BTjBypZvV73nP1NVZWUtU+BFgqCQE9e9wstenVa9Af7WN8C35eM0Qs3/VGzp/y3nSq005aFy/TQud9SXWXVvvdrDYPy9bM/1+SLb1S4prb9irTQ1D/eqel3PChJMvvom/H96IT05d/X2o9mt1d5rZJ2oeW5G36p+qrqfQeW3Rg/0NrpczT9X+l7oFrw/Bua/+xr+wwrX+E4eubG/1F1mh6oAt/X49f+bNfsWUs/s2PVer3xhzsTXJ29Zjz0jJZNndHi7bPRizf/WSVrNiSoKrtF6uo15bpfRC85jGFdq9i0VS//z98SXJ29Zjz4lFa8PzPmde3lX/1FJavT8yyVH4loyjd+Hg15MaxrpWs36rXf3pHg6uw19z8va8FLb+8Kvy317Pd/q8qt6TmxYIzR09f9XJG6+hb3zfi+Nn66QB/848EEV2evdR9/qo/+/kDLPxAYmcDo6at/okhdfeIKayMxhZb+/fvLcZyvvG688cZE1demNs5dqI2fLthn8twbY4xm3vt4szOY+6uP/j1FjhdjvjVGQX1Yc9J0pmjx6++pYtPWmA9Sxvc1e/JzqqusSlBl9jLG6IO7H4vvw66jGQ893bYFpYj5L76l6pKyZs8WfFng+5r79KtpOQPe2nXt4wfTc11b+MpUVW7ZHvN+LfB9zX78BdU0XI6dbj68a7IcN/Y5YhMEaXup9boZc7V18YrYx2tBoE/uf0J+OD2vjpl5z5To/cYxMEGgmpIyLXnp7QRV1XZi2opmzpypTZs27Xq9+eabkqSLLrooIcW1tVkP/KfhxrfYVW7epmVvftDGFdlv84JlWvfJvJhnI6XowODjex5XEOMBbn8w/Z7HYw96DcK19Zrz5IttXJH9Vn4wS9uXr2nxDO7ujB/o4weeSomZorb24d3xr2tBOKJZk59v24JSwIr3PtGOleviXtdmPPiUwrV1CajMbh/GM4HVIFIf1uzHX2jjiuy3Ye5CbZy3KOagJ0UHk9PvSc8J00/ufzLu8Vr1jjIteumdNq7IfpWbt2nxy2/Htb44rquZ9z6egKraVkx7ny5duqh79+67Xi+99JIGDRqkY49NjSeqLHvz/WavwW2KG/K08p2P2rgi+62YOj2uGaJGZes2qTTNLtsJgkAr3vskrqAnSY6kFVNntG1RKWDZ1I/jPkhJUk1ZhTZ+vqQNK7JfuLZOa2fGN6kgRScWlr4zvY2rsl+r17XySm1Ks3Ut8H2t/HB23OuaJC179+M2rCg1LHt3etxBT4o+Unr78jVtWFFqWPbWB60ar614N/32a2s+iP1y10YmCLRh1meqr6pp46raVtyPPK6vr9fkyZP14x//eJ/Px66rq1Nd3RczUhUVybvJuDVPATOBUU15+t0gXVNWIcdz45olalSbZjeW1++sblW/jDGq2lHahhWlhpqyCqmVz9qvLUuvy09q2+Bym6qSNF3X1Lp1Ld0emFFbvjOuM1O7pPF+zXVd+a0Ie+m2rklq1SXSgR+kZc9qyxuOoa3YTmvLK5SZZ+93ecUd/59//nmVlZXp6quv3udyt912m4qKina9+vTpE++vbDU31IqvpXEchTLT74uyvKxMqRXHKUny0qxvXht8UWQ6Pmvey8xo5TCybXqfStriO6TScr+WmdHafJx2398VymK/Fo9QZmZrD6Fp+SWdXkb84zXHcdK0Zxmtm1iQ/cfQuEPL/fffr9NOO009e/bc53I333yzysvLd73WrVsX769stQ59erRiJteoqE+PNq0nFXTo3UNBJBL35x3HUWGPLm1Ykf0ysrOUU1wU9+fdkKcOffe9Xe2Pivv0aPW12x16p9c2ml2Ur4zc7Lg/73qeOvbv3YYVpYYOvXsoaMXMd+P/I51k5OYouzA/7s+7IU/Fabhf69Cne9yXOUmSHEeFPbq2XUEporBX97g/6zhOWo7Xilq5fWXkZCung91fBBtXaFmzZo3eeustXXfddc0um5WVpcLCwj1eyXLolRfE/VkTGI265Kw2rCY1jDjzRGXkxDcocj1Pw087TrkdO7RtUSngsCvPj//m6IivMRPPaeOK7HfwhafLceLrmeO56nf4weo0IL0G4K7n6bDLz437/ozA93XY5ee2bVEp4JCLTo/76jDHddX3sNHqPCh1vkW6LTiO07Bfi3Ndi/gae3n67dcOOucUheKcvXY8T0NPnKCCbp3buCr7jbnyfDlufBtp4Ps65LL0W9f6Hz1O+d3jmyR2PU+jJp7TuiuS2kFcI4QHH3xQXbt21RlnnNHW9STUqEvPjOvUl+N5GnzikSru1ysBVdktKz9XY75+bsyP0JOiO47Dv3lpAqqy3+HXXBTb99o0chx1PWCQ+o4b3fZFWa6gaycddO7JcQ3AjR9owvUTE1CV/Y74xsVxz+R27N9bg44d18YV2a+wexcddNZJce3XTBBowg2XJaAq+x1x7UUxP4JWkuQ46jKkvwYcOabti7JcdlGBDr3s7PjWNd/X+G+l57p2yOXnxhWQHddV/6PGqvOQ/m1flOVcz9PY6y6NK+wFvq8x19j/JOCYQ0sQBHrwwQd11VVXKWR5Ivuy7MICHfndq+KYYTM6+iffTERJKWH8DZdHrwGPYUNwPU+9Dj1QA49JvwGRFB0Mjr7wtNifvGaMTvr5t/f5cIv92bHfvzrme6jckKdOA/rowLNOTEhNtut+wGAdcNpxcZ3ZO+ln18ttxdMBU9mxP7gm5s80Xk530FknJaAi+3UZ1E8HnXtKfPu1X9yQtvu1Cd+5IjoZE8Of3/U89ThwmIaeeGQCK7NXXqdijbvu0pgv6TfG6NibvpWgqux3yBXnK6djcUyBz3FdDT/rJHUZPjiBlbWNmI9Wb731ltauXatrr702EfUk3HG/+LZGnH1yyxZ2HMlxdPYdv1WfcQcntC6bdRrYVxMfvV2O67XoYOV4ngp7ddMVT9yRtgMiSbrgjlvU+5CRMR3gT/zFt3XQuackrijL9T54hC759+93bXvNcT1POUWFuu65f8d9Ccb+4LJ7/qhuwwfFFFyO/u6VaXm5TqM+h47URf/6XUzrWnZhvr753L/T8ibfRpf863fqNXp4TOvaCTd9U4dceHoCq7Jb16EDNPGhv8hxnRZN/rmep4JunXX1U5PiOkOzvzjltz/SkJOPallwaVjmzL/8UgOPOTzBldkrt1OxLnvqX8rIzW5RcHFcVz0OGaGzJ93aDtW1XswjylNOOUXGGA0dOjQR9SSc47o6/94/6Ygbr5TjunsdUDb+LKsgT5c8crtGX5p+97J82ZATJ+ia5+/ZdZPW3vrWuHPtM/Yg3fDWZOV36dSuNdomIydb171wrw447ThJavKyJ8d15XqezrztZzrpF99uxwrtdMhFZ+jKyX+LPmmoiQFl47rWaWAffe+dx9LuXpYvyy7M17dffUiDj46e2WxyXfNcOY6jU3/9XZ35+5+0Z4lWGnPpWbri4b8olJUZPQuwl3Wt8cDfsX8vfe/dx9VpYPKegGmDzLxcXf/SAxp20lGS1OSg2nFdOZ6rM279sU799Xfbs0QrHXDacbr22buVVRB9mMHewktjL3scNEzfeWdKWt6AvzsvFNJlj/5NY666QHKcvQflhm02IydLF977Jx12zcXtXKV9uh80XNe8MVkdG+6721t4cbzomb8R556irz93rzJy7X3M8e5S6/quNuJ6nk753Y91xA2X69PJz2nW/f9R1faSXe93HzVc4755mUacc3LcN6HvjwZMGKOfLXxTC154S9PvmaL1s+bves8NhXTQeafq8OsuUZ/DRqXtZQBflpmXqyse+7s2zF2ojx/4jz594kX59eFd7xd066zx37pMY684XwVd0zvk7e7AM0/Qr5e8rTlPvKgP7p6iHSvX7vH+oGPHacL1EzX85KPSeiZydzlFBbru+X9rzYy5+vCexzX/+Tf3eBpbTodCjf/GxTr86gvT8ilOTTnonJM16JjDNfvx/+rDu6eo5Etfhjvo6LGa8K2JGn7q0fJS7JLoRMnKz9W1/7lT62Z/ro/ue0Jzn3pFfviLp0zmd+mo8d+8VOOuPF9FaT7w3t3gY4/QLxe/rc+efV0f/vsxbfps8a73HM/ViDNP0PjrLtWAo8ZyDG3gZWTo7Nt/o6N/cI1mPfS0Zj309B7fT9Vl6AAdfv1EjbrwDGXl5yaxUrt0HjJA13/4rNZ+OEuz7ntCS1+fpqBhG80uLtLBl5+rQ6+6UMUDUmsSxjGmlQ91jlFFRYWKioq0bt26pD5JbHfGGIWraxWprVVWQb51z6mORCJ6//33dfTRR1t1H1FNeYVqSivkhkLK69TBuoBnY9/qq6pVsWW7wtW1yi7MV1GvblYNum3smTFGJas3qLq0XF5GSIXdO1t3Fs/GvlWXVqh8w2bV19Qou7BAnQb0seoSOht7FgSBSlatV3VZubxQSIXdu1j35CYb+1ZbWaWKTVtVX12jnMICdejbw6qAZ2PPJGnn9hJVl0T3a/mdOyqrIC/ZJe3Bxr4FQaDa8spdx1B61jLGGEVqauWGQtaNcSWppKREAwYMUHl5+T6zgT0dTSLHcZSZl2P1t4DaKKeoUDlFdgTPVJGZl6vOA9PrUamt5TiOOg3onfaXgMUqt7hQucVsn7FwXTftHmXcFrIL8pRdMCDZZaSc/M4dld+5Y7LLSCmu6yq3uEhqxXehpSPHcVLmErB9Sd+7pAEAAACkBEILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYLZSsXxyJRBSJRJL161NKY5/oV2zoW+zoWXzoW+zoWXzoW+zoWXzoW+zoWXx832/Rco4xxiS4FknSpEmTNGnSJPm+r6VLl2rKlCnKzc1tj18NAAAAwELV1dWaOHGiysvLVVhY2ORy7RZaGlVUVKioqEirVq3aZ2H4QiQS0fTp0zV+/HiFQkk7OZZy6Fvs6Fl86Fvs6Fl86Fvs6Fl86Fvs6Fl8SktLNXjw4GZDS9I6GgqF+AuNET2LD32LHT2LD32LHT2LD32LHT2LD32LHT2Ljed5LVqOG/EBAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAVgsluwAAAGxijJEkOY6T5EpSQxCJqK6kTH5tnTIK85VZVEjvmmGCQKVLV6q2pExuRkj5Pbsrv1f3ZJdlvR2LV6hk+WpFamqUWVCgHmNHKbdzcbLLstrOLdu1+t2PVFNaLjfkqaBndw08cYJC2VnJLi1mhJYGxhgpEpGTkZHsUlKCMUb1mzYqUloix/OU2b2nQh06JLssqxljVL1okaqXLFVQWyuvoECF4w5TZteuyS7NWsYY7fzsc+14a6rCJaVyQiFl9+mlLmd+TVnduyW7PGtVzl+kjY8/p8rPFymoqVWoQ5E6HXukul90lrK6dUl2eVaqWLJCax57Rhtfflvh8go5nqucXj3U79Jz1efCM5XZsUOyS7RO1fpNWj7lOa2Y/Kzqyyt2/bxw6EANveYS9T/3awrl5iSxQvvUlpZr6ZMvauGDT2rn+k17vNfjyDEacfXF6nfKsXIzGJ418uvDWvbCm5p7/xPaPOfzPd5zQ56GnH2KRl97sXoeNjpJFdpp/Yy5mnPfFC196W0ZP5DjutGxrjHKKirQ6Csv0CHXXqLCFArLab1VmEhE1Z9+osp3X1fdyqWSMXIyMpV76OEqOP4UZfYfzGzRl5hIWOXvTVPpqy+pbs3qL95wHOWPO0IdTztTuSNGJq0+GxljVPL6m9r69NOqW7NWchzJdSXflxxHhYePU7fLJypv+LBkl2qV0g+ma92d/1bNqjWS50kNs9+StOG+h9XhqCPU78ffU3bPHkms0i5Vy1Zq8c9uVdWipdGe+X70jXUbtXPBEq3514PqetYpGvK/N8ljMClJCldW6dMf/Y+2vvuhHM+TaeiZifiqXrNei/5vkhb/9W4N+9G3NOj6KzgmKLpPW3DHA/r8b/fIcRyZINjj/YplqzTr5ts074//1FH3/J+6HTk2SZXaZdP02Xrj6h8rXFW9x/6s0eYZn2rTR7PV6aDh+trkO5TTuWMSqrRLTUmZ/vv1H2jz7Ply3K/e0RBEfC194Q0tefZVjf3eNZrwyxv3ulw6Mcbo47/frw9um9SwT4tun7tvp3XllZr5r0c096GndP7kf6jP+DHJKjcmafs3G9mxTRt/91Ntv++OXYFFkky4XlUzP9TmP/1GOx66SyYSSXKl9vCrq7T21v/V5rvvVN3aNXu+aYx2zpyhtbf8StufemLX5RXpzgSB1t3+N6376+2qW7uu4Yfmi8GkMar4ZKaW/fBHKp32XvIKtczW517U0p/8UjWr10Z/4PtSEHzxMkZlH83QgqtuUPWKVckt1hKVny/Wp5d8S1VLl0d/0LiONWro3dYX39C8K78rv6q6/Yu0TLiySh9der22TpsuSbsCyx6MkYlEtPjP/9KiP93ZzhXaad4f/6nPb/93tDdfCiySdh1Pw1XVmvr172nzB5+0c4X22TzjU7162Y2KVNfsNbBI2jW4LFm4TC+ee53qyir2uly6qK+q1jPnf0tb5i6UpL2va4pOMEjSrH8+qA9u/Ue71WerGXc8oA9umySpiX1aA+MHClfV6KmLv6NNny5or/JaJS1Di19Rrs1/vkWRbVuiP/jyDqRhw6ia8b62P3wXA3BFV/z1//dH1SxZ3PCDvfSkoW/bn3pCpS+/0I7V2WvTgw+r5LU3ov/R1HoUBJIfaM0f/6Sdn81vv+IsVfbxJ1r1p9uj/7Gvbc8PFKnaqcXfu0nh8vL2Kc5S4dJyzb/uxwpq6yR/7wf2XYJAOxcu1eKf39o+xVls3s9/r8qlK3ftu5qz8r7HtP75VxNcld3Wvfy2Ft8zuWULB4FMEOiDb/5UtTtKE1uYxep3VumNq38k4wdNDrx3Z3xflWvX6/2f/aEdqrPXe7f8TTuWrtznwPvLZv/rUa16+8MEVmW3zfMW6v0/tHxyxQSBgnBEz1/9YwUpMEkfc2jZsGGDvv71r6tTp07Kzc3VwQcfrNmzZyeitoQpe/Ep+eWlzR+ojFH1Jx+qdsHcdqnLZpUzpqtm4YIWH9y3TXlU/s7KBFdlt3BJibY+9VRMn9l4/wMJqiZ1rP/3g9FL6FrCDxQuKdW2/76S2KIst/mZlxQpr2jx9qkg0I633lP1itUJrctmVWvWa/Pr77a8Z5LkOFp+18NpPZG16O5Hope3tpQxitTWatV/0ncia8Wzr6q+YmeLAksj4wda/co72rlhcwIrs1dtWYUWPfHCrrNPLeV4nubeMyVBVdnv0/uflON5MX3GBIF2btqqFW++n6Cq2k5MoaW0tFQTJkxQRkaGXn31VS1cuFB//etf1SGFbsAOaqpVNX1ayw9UrquKd19PbFEpoPTVl1o+kFR0pqh86jsJrMh+O159TQpiGNwEgaoXLlLNqvS93KlqyTJVLVy87zMsX2aMtjz1XEyzcfsT4/vaOPnp2HomSZ6njU88n5CaUsGaKc/J8WKctzNGO5evVunseYkpynIl8xer5LNFsQU9SQqMlj70HwVpuI0aY7TggSdjOn42chxHix97LgFV2W/Rf16SH4595t/4vtZM+1hlq9cloCq71ZSWa+Ezr8Z1LHQ8T3PufzIBVbWtmPbY/+///T/16dNHDz74oMaNG6f+/fvrxBNP1KBBgxJVX5urnjdbJhxu+QeCQLWfz1VQXZW4oiwXKS2JXhYW40Cy/L1piSsqBZS+MzWugWTZNPtnOxKl5N33ozeQx6h+67Yv7uVIM1VLV6hu89bYP+j72vZa+k4sbHrt3ZhncSXJCXnanAIzkomw8Z0PY57FbVSzeZsqlq5s44rsV7Vxi8qWrYr9WKDoDPia16e2fVEpYNVbH8TVs0Zr3p3ehtWkhvXT5yiIZXy7G+P7Wvv+J9ZfIhZTaHnhhRc0duxYXXTRReratasOOeQQ3Xvvvfv8TF1dnSoqKvZ4JZNfURbbqe3Gz1Wm7w1xkfKyOD+XvtcwS63pW/renxEpL5fifDhTpCw9+xYujf/PHSlP30s4w+Vx7tNNKz6b4urLyiU3/qenpeON5a39M9e1YvtOZTWtuAfK8dy0XNdqW7tfMkZ1lXZP0Mc0el+5cqXuuusuDRkyRK+//rpuuOEGff/739cjjzzS5Gduu+02FRUV7Xr16dOn1UW3hhPKiCu9p/P3tzih+P7s8X5uf+GEYn+iuOM4ab2uuRkZcuJMLenaNzcz/j93On8XRNx/dsdpVc9TWWv/3F4abqNeVmarPu9mtu7zqcprzbpmTFpuo22xfdm+jcYUWoIg0KGHHqo//vGPOuSQQ3T99dfrm9/8pu66664mP3PzzTervLx812vduuReZ5jZd0DMocXNy5dX1CExBaWAjC5d5WTF+M2prqvsAQMTU1CKyBk4MOazeiYSUU7/fgmqyH7Z/fvK+HGcnnak7D69276gFJDdt3dc18vLdZQzoG/bF5Qi8gb0i+usuwl85fVLz3WtoH+fXY+XjZnjKK9P+n2nUm73LnJC8V1S53iuigYkd6I3WToM7Bf3pYjGD1TUt1cbV2S/wj49W/X5rKICZeTZ/f1dMe2xe/TooREjRuzxswMOOEBr165t8jNZWVkqLCzc45VMWYOGKqN7DCuz46rg2JPleGk8I5mVpQ7HnxTbAT4IVHzqaYkrKgV0PvusmG9YdbKy1OH44xJSTyrodMqJcjJinFn0PHU48oi0/ab3rK6d1fG4I2O/Fygw6nnZ+YkpKgX0m3hu7DeUS3JcV73OTc99W58zTpIXx8y/43nqcfyRyknDbTQzP08Dzz4lrgG48QMNv/y8BFRlv5ETz4n74SpZRQUacMoxbVyR/XqNG60O/eObxHI8T6O+fr71X54bU2iZMGGClixZssfPli5dqn79Umdm2HEcFZ52TksXlhMKKf+YkxJbVAroEEsAcV1l9uyl3ANHJa6gFFB42FhldO3a8rDnuup02tfk5dg905FIofw8dTnra7EFZN9Xt4vS88DeqOflF3z1yyT3xXHk5eWqy+knJq4oy/U47QRlFBbE9BnH89TzjJOU1ak4QVXZLbMwX/0vPD32R6r6voZedXGCqrLfiKsuimsAnt2pWP1OPa7tC0oBvY8cow4D+8Y8AHc8VwddeaFCrbwsLxU5jqNDr7ssrs+aINDoKy9o44raXkyh5Uc/+pE+/vhj/fGPf9Ty5cs1ZcoU3XPPPbrxxhsTVV9C5B1+tApPPnPfCzmO5Hrq8p2fKFTcqX0Ks1hWr97qccN3m1/QdeXm5Kr3z39lfWJPNMfzNOB/fh29t6W5QbjrKrt/P/W49up2qc1mfa7/hrJ69mjZmQPHUZdzzlDREYclvjCLFU8Yp65nndLyA7wxGnLrz+XlZCe2MIt5WVk68Lc/bfHyjusqVJCnYT++PoFV2W/EjdcoozBfTgsnFhzXVfdjjlD3Yw5PcGX26jrmIPU99dgW96zRuF99P23vO3McR0f95vsxXc7veK6yigp18HWXJrAyux146Vkq6NE1pokFx3U08pIzVZwClyLGtAUddthheu655/T444/rwAMP1K233qq///3vuvzyyxNVX0I4jqMOF1yujhO/8cW9Kp4XfTXsVDL7D1L3m/5XOSNGJ69QyxQdd4J6/uincvPyoz/YfYDU0LesPn3V/7Y/K7NH666t3F/kDhuqwX/9s7yChhndLw8qG3YseQeO1OC//F9an2VpFCoq1Ii7/qac/g33W+xtIN7Qty7nnqkBP/sRAdlxNPQPv1SXrx0f/UFTg6OGfdzQP9ysrqdzBrnX2afowFtukhzt8ztbHM9TqKhARzzyT+X2Tu99W16v7jru0X8qoyC/+YGR66jToQdpwt1/innAvj9xHEcnTPqDuo07uMVPXxv78+9o6CVnJbYwyw0+/QQd87uftGhZx/OUkZuj8564U/nd0+8yxEZZhQW66Km7lFXYgu1TkhxHfSaM1Sl//nXii2sDMUf4M888U2ee2cxZihTgOI4Kjj1Z+UedoJrPP1Xd8iUy9XVy8wqUe8g4ZfZJnUve2lPh+AnKHztOO2dMV/l77yq8Y4ecUEhZffup+OSvKXvI0LQfQH5Z3vBhGvnYIyp7731te+6/qlmxQgoCOZkZKho/Xp3POVt5B46kb7vJ7NpFBz54t3a8/a62PPmsqhYv/eJNz1PH449RtwvPUcHBo+hbAzczQ8P/+lt1OulYbZz8lCo+/XyP952Qpy6nn6ReV16sggOHJ6lK+/S/4kLlDx6g5f96SNs/mim57q4BtvF9uZkZ6n3+GRrynauU07N7kqu1Q8eDhuuUlx7W3Nv+qQ2vTdWuufAgkON5Mr6vjMJ8Db7iQh34g+ta/QSt/UEoJ1unTblTM2+7U4sefUZ+XX30jYYzCY19y+vZTYfd/F0NPj8975v6skOvv1y5XTpp2m/+oprtJXI8d4/vV2rsW9eDhunUSb9Xx8H9k1esJToNGaAr3nhML1z3U22Zt2hXj3bXuI87aOI5OulPN7fuaW3tKD3PO+7G8Tzljh6r3NFjk11KynAzMlR41DEqPCr9bnSLl5uVpY4nn6SOJ58kY4xMJCLX8kcLJpublakup5+qLqefqrpNmxUuKZWTkaGsbl0VKkruAz1s5biuup5xkrqecZKqlq1U1eLl8mtrFSooUIfDD1FGcYdkl2ilzuPHqPP4Mdq5aq02vzFN9SVlcjNCyundQz3POEkZBfnJLtE6+X176ai7/qSaLdu08j8vqnzpSvm1tcooLFC3I8eq7xknycuO8amT+zkvK1NH3PJjHfqTb2n5M69q9avvqGZbidzMDBX07aVhl5ylXseNlxvnU7P2V8PP/5qGnHWiVr3xnuY9+JRKlq5UpLZOmQX56nvMOI26+iJ1Gz2i+f9RGunQr5eueOMxbf50gT598D9a9sq7qt9ZJcd1lNe1s0Zdfp5GXXG+Cnp0TXapMUn70AK0t3T/LpZ4ZPXorqwezHLHIm/IQOUNSe/Hjscqf0BfDb7+imSXkVJyunXRyO9dm+wyUkpmQb5GXH2RRlx9UbJLSRleRoYGn3GiBp+Rvg8PiZXjOOpx6IHqceiB0j+jN9vLcVL66gRCCwAAALAf2x/uK0v9PwEAAACA/RqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArBZK1i+ORCKKRCLJ+vUppbFP9Cs29C129Cw+9C129Cw+9C129Cw+9C129Cw+vu+3aDnHGGMSXIskadKkSZo0aZJ839fSpUs1ZcoU5ebmtsevBgAAAGCh6upqTZw4UeXl5SosLGxyuXYLLY0qKipUVFSkVatW7bMwfCESiWj69OkaP368QqGknRxLOfQtdvQsPvQtdvQsPvQtdvQsPvQtdvQsPqWlpRo8eHCzoSVpHQ2FQvyFxoiexYe+xY6exYe+xY6exYe+xY6exYe+xY6excbzvBYtx434AAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqhZBcAAACQbkxdjUxVuRzXkwqK5XgMyZrjb92gyJqlUrheTm6BQkMOkptXkOyyrGWCQPXLF6p29ofyy7bL8UIK9eqnnHHHKdSle7LLixlbCADsrwI/+vJCksOJ9WaZQE7JRjnVZZLjynToLpPfMdlV2a2mUs7az6VNy6VIvZRTINNnpNRrWHS9w1cE2zfK//xDmTWLJGOiP8zMljtsjLyR4+Vk5yW3QAtFNq5W7bvPK9i8TnKc6CsIpHeeU2jYwco+/hy5ufnJLtMq4U3rVP7YvxSUbIvu/00Q/fnaFar58C1ljjhEhRd9Q25WdpIrbTn2KEB7Mya683Dc6I4X+xYECtWUyq3fKclVJKeDgqx8ercPTk25QltXyK3YKkeScUPyO/ZWpMtAKZSZ7PKs5OxYr9DyGXLCtTKOIxnJkVFQ0FmRYROkbAZEX7FhsZxP35CMkaPo4NvUVMop2SgtmS4z/gIpvzjJRdolWL1AkWnPRv+jMbBIUn2tgs8/UrBivjJOu1pOAX1rFFm9RNXP3b9r0B09hjb0zgSKLJmrqk1rlDfx+wSXBpHN61X27z/JhMPRHzT2ToqGPUn1i+aq7P6/qPibP5OTkRrHhZim3m655RY5jrPHq3v31Du9hDZiAnmRWnl+3Z47X+ydCZRZvV25ZWuUV75WOeXrFKqroHf74IRrlLP5c2WWrlGoaodCVduUs32ZsrYv37XjxZ7cnTuUufxjuRXb1BjrnCAib/tqZa74ODobjj04JRsUWjRNCtdG/3u3QbhTuUMZ896Q6muSWaJ9tq2VM+e16NkpfbEPcxpeqquS89HTUn1tsiq0TlC6JRpYTLDnILKRMVJNpcJvPibD/k2SZGqrVf3CQw09a+JYaQKZilLVvPZEu9ZmK2OMyp+4JxpY9raefbGgIhvWqGrqK+1XXCvFfL3AyJEjtWnTpl2v+fPnJ6KudmbkysiT2WPni6Z5kTpl15UpI1KtzHCVsuvK5AR+ssuylzHK3rlZobrKLwZDxldW9Y5ocMFXmUDZ25bLCSK7BkKNg3CvrlKZ5euTWJylgkAZa+dJe9mXOZKc+hqFNi1JSmnWMkahFTMlfbF+7c6RkcK18tYvbN+6LOcsmS7J2WvPpGjwU12VtI6+NfIXfNz8QsZIFTtk1i9LfEEpoH7BLCkcbn5yzwTyVy1SULa9fQqzWHjNcvlbN+47sDQyRjUzpspEIokvrA3EHFpCoZC6d+++69WlS5dE1NWuPBllONHQkkFwaZ4JlBGp2mMQKRllhKuSV5Pl3EitvEjdHgf4xn/PrC1r2c4lzXg15XKDcBMDSSlUtUMKUmNH217cyq1y/PqmB5Iy8so3Sn64XeuymVO+RU5dVZM9k6J9c7csj94fBGlnqZzSTS06VjprPmuHguxnImGZlfNbtq93HPlL5yS+qBQQnj9DaumYzHGjISfN1X46XXJbPrw31TtVvyI1JhdiDi3Lli1Tz549NWDAAF166aVauXLlPpevq6tTRUXFHi/beA3/bLxE3iO07JPbMPO9O0eSZyJc6tQEL1Lb5FrlmEAug8ivcOurZfYxlHRk5Ia59GR3bm3lPnsmNVz6VF/dThXZz6mpaNEe3/EjUrgu4fWkhOryFi3mSFK1fcf8pKitbnnoNUZmZ2li60kRQWVZTMubnS1bN/dnQUVZzJdPBxVlCamlrcUUWg4//HA98sgjev3113Xvvfdq8+bNOvLII7Vjx44mP3PbbbepqKho16tPnz6tLrqtGe051m7uoI+994e40jTTzJObmns/Lbmuml2r6NsejOupRVui4zW/TLpwvJbv8V36Jim2PsQw47tf82Jbd3j8cZQTS98c8cQ6SU4oQ02N05oUykhILW0tpr3JaaedpgsuuEAHHXSQTjrpJL388suSpIcffrjJz9x8880qLy/f9Vq3bl3rKk6AyB4XOUlcALBvgRtS4Lh7DI2MpIiXzROdmuBnRh9h+eXhpJHke1kyXmrsMNpTJKdDk7tdIynwMhVk5LRnSdYLCrru81BlJAWZuTJZPFK1UVDco9mYZyQF+R2ljKz2KMl+HbrJtGBwaBxH6ty3HQpKAdl5UkELH5/tOHJ6DEhsPSnC6z2w5ZNTQSCvZ/+E1pMKQn36x5xZMnr1S0gtba1VUyB5eXk66KCDtGxZ0zeMZWVlqbCwcI+XbYwc1ctRnXEUlquY/7bTjeOoLrNQvpsZPUslR5FQjiIhBpBNMW5I9bmdo/+++8txVZfXOZmlWctk5Cic22mvQU+S6jv0JiR/icnOl1/YtelLESVFug6ib7vLylXQqW8zlyJKQa8D2q8m24Uypb4HRkPJPjjGyAw4uH1qspzjOPIOGNeyhY2RN2xsYgtKEZkHH9Xyez6zcpQxdFRiC0oBOYce1fJ9vOMoo/9Qhbr2TGxRbaRVoaWurk6LFi1Sjx492qqeJNrztnI0w3EVzsxXbXZH1WYXRwMLA6F9imQVqKawl8JZRYpk5qs+p6NqCnvLeKnxfPRkqC/uq3BhD5ndLmcKQtmq6zRIfk6H5BVmsXDvUQryO0lqfIaYsyskh7sNVVDcK6n12cgfMk4mt2hXnxo1Bhm/xzAFnVNjJrK9mGFHSHnF+wwupv9oqbN9l4Qnizv0UDkduzd7rHRHHcX3tDTw+gyS129oi8YX2Uef3nBpVHpz8wuUc8QJLVvYSHknnJXYgtpQTBf/3XTTTTrrrLPUt29fbd26Vb///e9VUVGhq666KlH1AfsV42UqnMs3bLeY4yhc2EPhgm5yInXRbyn3MgnI++KFFO4/VpHqMnnlm+X4YQWZufKLe0mZnA3dq1CWIqNPkbtpmbxNS6OP6pVkCjsr0vMAmU6c1fuKjGyZCRfL+XyqzMYl0bMqik79mYxsmSGHSQMPpW+7cUIZCp16pSLv/kdm8+o9vqW8sU/uqKPlHXxc0mq0jeM4yj37KlU//6D8dcujfdr9JmTXlYJAWUedpszRRyavUMvkn3aR/PIS1S9o4il0jisZo4Lzr1Tm4NQ5ixxTaFm/fr0uu+wybd++XV26dNERRxyhjz/+WP36MQMFIIEcV4b7V1rOcWTyihXJY7a2xbwMBb1HKOg9QvIj0YM6N5HvW2a2zKFfk0YeI7N1jeTXy2TnS13789CCJjhZOcr42lUKtm9QsHSOTMUOyfXkdu0bPROTW5DsEq3jZGYr98JvKbJ8geo/fV/++oan1nohZQw/VJmHTJDXrXdyi7SM43kquuwG1c58T9Ufvil/+5Y93s8cMlK5x56mzAFDk1RhfGIKLU88wbeNAgD2czyBKDZZuVKf1JmttYHbuZfczlyq2VKO6ylj6ChlDB0lE/hSJCJlZMrhTF6THNdVzuHHKXvcsYpsWCO/bIecUEihbr3lFXdKdnlxYc8MAACAlOC4npTJmbyWchxHGb37K6N3/2SX0mqc+wYAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWCyXrF0ciEUUikWT9+pTS2Cf6FRv6Fjt6Fh/6Fjt6Fh/6Fjt6Fh/6Fjt6Fh/f91u0nGOMMQmuRZI0adIkTZo0Sb7va+nSpZoyZYpyc3Pb41cDAAAAsFB1dbUmTpyo8vJyFRYWNrlcu4WWRhUVFSoqKtKqVav2WRi+EIlENH36dI0fP16hUNJOjqUc+hY7ehYf+hY7ehYf+hY7ehYf+hY7ehaf0tJSDR48uNnQkrSOhkIh/kJjRM/iQ99iR8/iQ99iR8/iQ99iR8/iQ99iR89i43lei5bjRnwAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgtVCyC0i2spqIpq0q1cKtVaoNByrKDumIvkU6rHeBMjwyXVPqIoE+31KlkpqIPNdRn6IsDeqYLddxkl2a1YwxKq31FfYD5WV6ys/0kl0SAACA9dI2tASB0ZPzt+r1pTtkjGQafu5Imrtppx6d4+qb43rp0F4FySzTOpHA6LVlJfpobYXqfSO3IaMERuqYE9LpQzvq4B75yS3SQpHA6NNNOzV7U5Uq6vxdP+9TmKlxvQo0qGN2EquzV2CMNlSGtbikRmW1vlxH6pmfqeGdslWcnba7rxapDgfaXBVWJJByMxx1z8tQyGVSoTnGGEVM9DIEj34BgDXS8qhvjNEDszfqvVXlX32v4Z9V4UD/+HCdvndkb43tXdi+BVoqEhg9MHuzlu2o2dWnwHzxfklNRJPnbVVlna+j+xclpUYb1fuB/rNguzZWhr/y3vqKeq2r2KEj+xToqL6sZ7ur8wO9tbpCW6sjcvTFtllZUqvFJbUa1SVHh3bLlcPZvT1EAqN5W2u0qSoiSbt65zm1GtYxSwOKMunZXgTGaGNVRJt2RlQfRH9WkOGod0GGOmZzRnRf6nyjinoj30iZrlSY6RCQWygw0W2UTRJoXquuf7rtttvkOI5++MMftlE57WP2hsq9Bpa9uXvGBlXX+80vmAbeXlG6R2Bpyn8X79Da8tp2qSkVvLasTJv2ElikLwbiH62r1KJt1e1XlOWMMXp7dYW2VUcH3ruvc43//tm2Gi3cwXq2u8AYzdhYvSuwSF/0yzfSwh11WlFWn5ziLBYYo0Ul9VpT8UVgkaTKcPTnu/cTXzDGaEu1r9WVgUrqjMrrjbbVGq2sCLQz3NyRIr3V+tK2OmlLnbS5TiqplyK0rFnGRPtUHzjy6VfaiTu0zJw5U/fcc49GjRrVlvW0izeXlaglk0BGUtg3+mBNywLO/iwSGH2wpqLZwCJJriN9uKYi4TWlgvLaiBa3IOhJ0vT1lTKGvbAkbaoKa0t1pNm+zd1SrUhAzxpt3BlWad2+J1mWlNSpzg/2uUy62VLtq6yu6Z6sLA+rNkLPvqy83mhvGdhI2lAVKMyocq/qfKk0vGdIqQukHXV7Xr2APRkjVQeuagJPdcZRdeCp1nfEYbN5RpJxXBnHadF4xFZxhZadO3fq8ssv17333qvi4uK2rimhSmvCWrStusU7BiPp/VWlCa0pFSzeVq2aFh60AxO9LyjMwEifb61WS8/6b6+OaGvV3s/IpJulJbUt6lt9YLSugjMHjdaUN7/+GEnrmzjzl6427mz+TMqWas64784Yo5K6fR9Iy+pTeXiUOJVNrG6BJFazptUZR1+MKqJHiLBc0bJ9M5LkhiTHjf7TTd3LXeMKLTfeeKPOOOMMnXTSSc0uW1dXp4qKij1eyVRWE/tp/tI4PrO/KauNrQe+karqCS0Vzcx6t3b5/VVFnd+i2SBHUiWXb+6yM9x8LxxJO9k2dzHGqLYFZwSqwvRsd0ZScy2p4UzLVxgj7evKOTbNpkWMI31lOsvIN9wQtE+O88VLkhw3Zc+2xHwj/hNPPKE5c+Zo5syZLVr+tttu029/+9uYC0uUDC/2lZsbCuPrAU/eifbAcdTi09f0LCqWPvCY7S9Ee9H8yhbHbnC/1pKusZ7tqSXd4EsD9m5f6xs9a1q0b0ZfXvvYMpuxx02hqRpXomLaPtatW6cf/OAHmjx5srKzW/aI1ptvvlnl5eW7XuvWrYur0LbSNS9T2aGW/7FdRxrcKTeBFaWG/h1ieyRvcU5I+ZnsfvsUZrb4UkTXkXrkZyS2oBTRKz+zRQciI6knPdule16o2b4ZSV3z6Fkjx3HUKbv5fVWnHPZnu3McRwXNrEaFmQwnv8xxpJx9XJ2zr/fSXaYbaM/IF/1nyEntgXjiGSnwvwgsgZ+yQS+mvfDs2bO1detWjRkzRqFQSKFQSNOmTdMdd9yhUCgk3//qpQlZWVkqLCzc45VMmSFXxw3s0KIb8aXo/RknDUmt+3YSoXtBpvp3yGrRiu5IOqpvIY9VlTS0U45yWhCSHUkHdM5RTgZHLEka2oLvrXEkdckJqWNOWj65fa/6FWbu84yBIyk35KgLI6M99Gom+GZ7jjrx2OOv6JTtNnlMyPaij4zGVxWEpNBeWpPrRR8Zjb3LcKRsJ2hY54w8Sblu0OLxXLpyJDkmkIKIFESUyrfix7R5nHjiiZo/f77mzp276zV27Fhdfvnlmjt3rjwvNXbqJw/pqAzPaXYA7jrS0M45GtaZMy2SdPrQjs0+S951omdZxvHdNpKilzmdPKjDPpdxJGWHXB3dj541ys1wdWSvpr+k1JGU4To6qjdfZLq7wixPo7tEA9+XN9XGnh3Wg++2+bL8TFfDizP3GPw0/mu252hkp0wuD9uLLM9R3wJXubvNGziSijId9cl3Wc+a4DpS50ypQ0b0zEquJ3XKlIoy+L6W5mS4RvleoAIvUK4XcKlrDPZ2R1CqiWmKsqCgQAceeOAeP8vLy1OnTp2+8nObdcnL1E+O7qu/vrdW4cDs9fIdR1Kvwiz9YEIfdrwNBnbM0ZUHd9Ojc7coMHu/JrdDdkg3HNZDORlMFzUa3jlHgSnWa8vL9ng8b+NJ7qJsTxcc0EmFWZwx2N3QjtnK9BzN2lylyi/dndo9L0Pje+WrKCs1JkraU5/CTOVnulpZVq/NVdHHRofc6FmYAUWxXR6bTjrleDosK1tbq31VhQM5jlSc7aljFoPvfcn2HPXJ9xQJol8umeFy/09LNF4mxklPoOXSdpQ0vEuefnvyQL20eLs+Xluxx2CyMMvTiYM76rShnZTN4HsPB3bL08+P7qPp6yr08brKXY9B7paXoQn9ijSmZ76yGBR9xYguuRpUnK2F26q1orRW9b5RQaankV1zNaBDFoOiJvQvylK/wkxtrgqrvM6X6zjqlpdBWGlGcXZIY7qHZEx0UsZ1xDrWAiHXUc/8tD0stkrIddJ3QAGgXbR6HzN16tQ2KCM5ehZm6Vvjeumy0d20sqRm10BycOdcnhi2Dx1zM3TGsE46fWhH1ftGnuvQrxbICrk6pEe+DunBJU2xcBxHPfIzRdti5zgOl08AAPYLTIxIKsgKaXSPgmSXkXIcx1HW3u4mBAAAANoQ1/EAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAq4WS9YsjkYgikUiyfn1KaewT/YoNfYsdPYsPfYsdPYsPfYsdPYsPfYsdPYuP7/stWs4xxpgE1yJJmjRpkiZNmiTf97V06VJNmTJFubm57fGrAQAAAFiourpaEydOVHl5uQoLC5tcrt1CS6OKigoVFRVp1apV+ywMX4hEIpo+fbrGjx+vUChpJ8dSDn2LHT2LD32LHT2LD32LHT2LD32LHT2LT2lpqQYPHtxsaElaR0OhEH+hMaJn8aFvsaNn8aFvsaNn8aFvsaNn8aFvsaNnsfE8r0XLcSM+AAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqhZBeQbEFg9NH81VqxYYdq6iIqyM3SESP7alDvTskuzWrVtfX6ZOF6lVbWyPNc9epSqEOH9pTjOMkuzWollTVas7lcNfURFeZmaUjvjsrK8JJdlvV21tSrrKpeIddRx4JsZdKzFjPGsF0CAFJe2oaWiqpaTXljrh56aZbWbyuXI8lxHQWBkSQdcWBfXXPmYTrtiGFyXQ74jVZtLNHDr36qJ97+TNW14T3e69e9g64+/VBddPxBKszLSlKF9jHGaObiTXr0zfl6debKXeuYJBXmZumyE0do4gkj1adrYRKrtE/ED/TuvPV69O3FmrF4866fZ4ZcnTN+oC4/YbgO6NsxiRXaqbbe16uz1+nx91Zq2cYKhSOBcrJCOuqArrrs2EE6bEhnQsxe7KyN6JU5m/TavM3aURkNyP265Oq8cb105LDO8jgO7FVFTVivztuiJZsqVRsOVJAd0pFDO2nCkI4KeVzM0ZTK2og+XLZDO6rCCrmO+nbM0WEDixViPdun6npfa3bUqDYSqDDbU/9OuWybLRQJjBwppfuVlqFlzeZSTfyfKVq7pVzGRAeQRpLZbTA5c+E6ffz5Wp1zzAjd/oOzlJWRlq3awzNTF+imO1+RJPm79arR2s1l+t2D7+jfz3+ix265REM4W6W6sK+f3v22Xvp4ubzdQnGjiuo63ffyXN338lz9/tpjdcnxI5JUqV3WbKnQdX9/W2u3Vn5l0qA+EujZD1foqfeX64xx/XXbtRM4W9Xgrbkb9KtHZ2tnbUSOIzXs3lRdF9Hbn23SG3M3akiPQv3zhvHq0zkvucVaIgiM7nlrpR6atlr14UBS9HggSWu2V+ndBdvUrShLv7lwhCYM65y8Qi2zszaif76xQi/P3ayIb+Q6UiDJdRw9P3uTOuVn6ppj+unCcZyB393Gslo9+tFavbFgm8K+kedE17fASJ3yMnT+mJ66ZFwvZYYIfLtbW1Kj5+dt0RsLt6ne/+I42iU/U+eM7qbTD+yi/CzGaV+2saJOry3arreW7tDOel+S1LsoS6eP6KLjB3dUXmZqHTvTbqvYUlKp83/xiNZv/SKw7E3joPzF9xfpB3974SuDzXTz3LQF+tEdL8sPzF4Di9QQ/Iy0raxKF/zyMa3ZXNauNdom4gf6zt9f08szlkvae9Br/LkfGN1831RNeXtBe5ZopbVbK3XxH17Rhu07JWmv215jL1+duVo33PGOwpGgXWu00Qsz1uqH985QVW1E0heBpVFjz1ZuqdSl//eu1mzd2d4lWicIjP7nPwv077dWqi4cRPdhu73vN6xWW8vr9P0HPtVrczfv7X+Tdsqq6nXdfXP0wpxNCvtGRpJvoutc43q2Y2e9/vLKMv3l5WX7PNamkyWbd+q6Bz/Va/O3Ktww8PZNNLBI0o6qsO57f41++Ph8VddFklipXd5bVqLrp3yuVz/fukdgkaRtO+t1/4frdMOUz7WxrDZJFdrHGKNnPtuiG55aqBcWbN0VWCRpfXmd7pm+Xtc9uUCLtqTWcSDtQsvP7nxF28uqmhxAfllgjF76YJGefGtegiuz1/qt5frJna+2eHk/MKqsrtO3//LftD5Y3fvyXE2du+Yrg8d9+c2D07Ro7fbEFWW5IDC6/h9vq7Im3KJtNDDS9EWbdOcL6bt9StKidWX61aOzvzLo3hs/MKqsCev6SR8q7Kd32Htw6mq9PGdTs8s1zoT/+onPtXB9ReILs1gQGP1kyudau71aLTmMPj1zoyZ/uC7xhVluS0WdfvzE56qu9+Xvo2/GSAs3Vuo3zy1O6+Nno1lryvWH15ZHJ/eaaIdRNCT/9NnFKqsO732hNPP851v18MyNktTkdloT9vWbV5drxfbqdqysdWIKLXfddZdGjRqlwsJCFRYWavz48Xr11ZYPZpNt7eYyvTNreYsDSyPHke578ZO03YE89sY8BTH+2f3A6POVWzRveXrOTIYjvh54dV6zA8gvcx1Hk9/8PCE1pYIPF27Uys0VMW2jxkiPvr1YNWk8M/nIO8sUyxU4fmC0bnuVps5vfsC+v6oN+3rw3VUxf+7haavbvpgU8snKUn2+vmKfA+8ve/C9NardbaY3HT0xY72q6yItCnqBkWauLtOcNeWJL8xixhjdOXV18zMxip6x2lFVr2c+Tc8xx+62V9XroYbAsi+Bid7ncucHa9uhqrYRU2jp3bu3/vSnP2nWrFmaNWuWTjjhBJ1zzjlasCA1LmmZ/PqcuG6qN0Zasmab5izZkICq7FYXjuixN+bGdXmc57p65NU5CajKfm9/ulo7Kmpi/pwfGD3z/hJVVNcloCr7TX5nSVw3CVbVhvXyJ7EPQPcHpTvr9Mrs9TFPxriuoylTVySoKvu9MW+LqupiG0j7gdHb87dqe2V6bp+S9PSMDYr1/vqqOl9vLtiamIJSQHW9r5c/2xJT0PNcR8/NaX7guT+bt75SG8vrWjz5Fxjp5c+3qj7NLxd+ffEOtfQoGhhpxY4aLduWGmdbYtr1nHXWWTr99NM1dOhQDR06VH/4wx+Un5+vjz/+OFH1tam3PlkW84G9Uchz9e6c9DvAL1i5VWU747tO1A8CvTlzeRtXlBre/XRN3E/oqA/7+mRR+h2sgsDo/fkb4tpGHUeaOm99Aqqy30eLtyoSy2ioQRAYfbJs+657YNLNO59vjensVCM/MPpgcXpewhmOBPpg2Q7FelWh40jvLNiWmKJSwKxVpaoNx9Y0PzD6YFlJWl/COXXZDnkxbqM763zNTfNLON9ZtqNFZ/QaeY703sqSxBXUhuJ+1ILv+3rqqadUVVWl8ePHN7lcXV2d6uq+mJWqqEjeyhTv4LtReSs/n4rKdsZ+tmB3O2vq26iS1FJeVRd3QJaksp3pN5O7s7Zl97HsjTFSSWX6bZ+SVLazXo5adAXFXlXU1CsvO/2eulOysy6m+80aea5UXpWe181X1Ebi6pkx0XsO0lVZTXwTA4GJPhq5Y15mG1eUGkqrwzGdnWqU7ve1lMc4EWUU/zra3mK+EX/+/PnKz89XVlaWbrjhBj333HMaMaLpx7TedtttKioq2vXq06dPqwpujYxWPELQkZQZSq1Hw7WFjFb+mdP1Of2ZIU9uKx7zmY5fntjaR3ymY88kKTPDjTuwSK3ve6qKd38emNYdS1JZZqzT3rvJStOeSVJGK74XIyNNj6GSlOm5cZ0NDbViPd0fxPpdP04cn0mWmLeGYcOGae7cufr444/17W9/W1dddZUWLlzY5PI333yzysvLd73WrUveU0T6dS+O+4siI4FR765FbVyR/fp0i//P7Ejq1SU9vzCxT9fCuHa2jfqm4RdNZmeG1KkgO67Peq6jft3Sr2eS1K9LftyfzcsOqSg3PWdxB3bLi+sSTmOk/l1yE1CR/fKyQuqQmxHz5zzX0cCu6fu9QAPiXF+K8zKUn5WekzGSNLBzfH0bEOfn9heDOucqll1bYKRBnVKjZzGHlszMTA0ePFhjx47VbbfdptGjR+sf//hHk8tnZWXtetpY4ytZJp5ySNzftxJyHZ17zMg2rsh+/bsX67ADesUX9hzp66ce3OY1pYILjx0e370Zkgb26KDRg7q2fVEp4JLjhsZ1hsoPjC4+ZkgCKrLf2MGd1bNjbotvvGzkuY4uPLJ/2p4NveDw3nFto12LsnT4kPT84lzXdXTBYT1jGhBJ0e3zvLE9E1NUChjWPV+Du+bFNJHlOtL5h/ZI6y/mPHVE55j2a64jHdA9XwNSZACeKGeM6BLTPS0hz9HxQzomrqA21OqjlTFmj3tWbHbakcPUIT/2mVzPdXTOsSNVXJieG8LVp4+JK+xleJ4uPP6gBFRkvwHdO2jCyN5xhb2rTh2VtgeqS44ZoljvznAd6cB+HTWyX/oOJC8/bpBiTS1+YHTJ0QMTU1QKGNazQCP7FMY0AHcd6ZIj+8T9kI39wTljesS0vOtIw3vma3jPggRVZD/HcXTR2J4x3Q/kOI7OGN09cUWlgI55mTp+aKcWb6OBkS44JL17JkmH9y1S94LMFvXNkXT68M7Ky0yNM3oxhZZf/vKXev/997V69WrNnz9fv/rVrzR16lRdfvnliaqvTWVlhPTjicfE9BnHkUIhT985v+mHDezvvnb4EI0c0DXmA/W3zz88rpC4v7jpksPlOk6Lx5Ke62hAjw46/+hhCa3LZt075umbpx3Y4uUdRQ/uv7jksMQVlQIuPmqAhvYsimkbvfbkoerXNf5Ly/YHvz7/AGWE3BYd3D3X0ZAeBbpsQt/EF2axbkXZuun0lp3VdB0pK8PT/553QIKrst+pB3XVCQc0f+ag8f2bzxiizvnpeenm7r57XH/17ZjTom30vIO76ZgUOWOQSJ7r6JZTB6kgK7TPvjmSRvcq0JWHpc5Z0JhCy5YtW3TFFVdo2LBhOvHEEzVjxgy99tprOvnkkxNVX5u75oyx+sZZLRvguI6jkOvq3psv0NC+XRJcmb0yQp4e/vWF6t215YOiC48bqR9dPCHBldlt9KBuuvP7p8jzmh8Uea6j7h3z9fAvzlJeduzXjO9PfnjeIbro6MHNLuc6jjzX0e3XH6PDhnVrh8rslZsV0r9vnKAB3Qr2eXld41sXTeivH52dfpe7ftnwXoWa9I1DlZ3pNblvcxpeQ3vk61/XHaqcFJmRTKQLxvXSTacPliM1+Z0tjqSC7JDuunp0Wt/P0sh1HP36zKE6c3R0X7W3e8UdRS/V+d9zhumUkel5ifCX5WV5uv3CAzSufwdJ+sqx1JGU4Tm6Znxv3XB0ek8o7K5nUbb+evYwje1TKEfRvnmus2s/lx1ydf6obvqfUwal1MMeYnrO5f3335+oOtqN4zi65bqT1b1Tgf4y5T3V10ckR3uctvVcR35g1L1Tgf75k3N0+Eg2hK7F+Xrx/12hmya9qjc/WS63oUeNnIYe5mSFdOP5R+h7F45P20ucdnfK2IF67Ffn6JaH3tOitTt2rVuNGgeYp4wdqN9dc4w6FeYkq1RruK6jW68arxH9OumeV+ZrU0m1PNdREBg5TnQb9gOjMUO66EfnH6oxQzi4S1KXomxNuelY3ffmUj353kqVV4d3PUXHmOjlYEN6FOqak4bqrHF92D4bjBlYrKd+NF6Pf7hWz36yQdVf+sLJ3p1ydOmEvjr/8F7KTtMn1O3NRYf31pgBxXp25ka9+OmmPb6HpEeHbF18eC+dcXB3FcVx4/7+KuS5+ulpQ3TRYb30/JxNen3BVlXV+XIkdS/K0nmH9tDpo7qpMIee7S4/K6TfnTVU60tr9NL8rfp8407VRnwVZWfo6CHFOnl4Z+Vlpd9j25vTtSBTvz55kLbtrNd7K0q1vbpeIddRv+IcHT2wOCWf6JeWf8uO4+jb54/X1792qJ55d74efmW2Vm8qUTgSKCcrQ4eP7KNrzjxMxx0yUF4KJdBE61CQo/t+cb7WbinTlDfn6ZWPlqh0Z61CnqtenQt16UmjdN4xI5SXwynt3R02rIdevu0SzVuxRY+/vUBL1pWoui6sovxsTRjZW5eeMELdipmJ3J3jOJp4/DBdeuxQvb9go96as1alO+sU8hz17JSv8ycM0uCeHZJdpnXysjP0g7NG6tunHaB35m3Usk0Vqq33VZCTofHDu2pU/2LCyl707Jijn5w1TN85dbBmLi9RSVX04N67U65G9yuiZ00Y2DVPN50xRDeePFAbSmtUGw5UkB1Sn445cT+pMx3075yrH54ySD88ZZAifiDXdVr1iPx00bs4Rzcc0y/ZZaScLvmZumD0/nE1QlqGlkYFuVm6+oyxuvqMsZKi3xDNjrZ5fbt10C++fqx+8fVjk11KShk9qJtGD9o/dhztxXUdHXtQLx17UK9kl5JSMkOuvjamt76W7EJSTE6mp2NGpO+lwPHKyfQ0uFt63x8Vr3R9eh8QD7aW3RBYAAAAAPsQWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWCyXrF0ciEUUikWT9+pTS2Cf6FRv6Fjt6Fh/6Fjt6Fh/6Fjt6Fh/6Fjt6Fh/f91u0nGOMMQmuRZI0adIkTZo0Sb7va+nSpZoyZYpyc3Pb41cDAAAAsFB1dbUmTpyo8vJyFRYWNrlcu4WWRhUVFSoqKtKqVav2WRi+EIlENH36dI0fP16hUNJOjqUc+hY7ehYf+hY7ehYf+hY7ehYf+hY7ehaf0tJSDR48uNnQkrSOhkIh/kJjRM/iQ99iR8/iQ99iR8/iQ99iR8/iQ99iR89i43lei5bjRnwAAAAAViO0AAAAALAaoQUAAACA1QgtAAAAAKxGaAEAAABgNUILAAAAAKsRWgAAAABYjdACAAAAwGqEFgAAAABWI7QAAAAAsBqhBQAAAIDVCC0AAAAArEZoAQAAAGA1QgsAAAAAqxFaAAAAAFiN0AIAAADAaoQWAAAAAFYjtAAAAACwGqEFAAAAgNUILQAAAACsRmgBAAAAYDVCCwAAAACrEVoAAAAAWI3QAgAAAMBqhBYAAAAAViO0AAAAALAaoQUAAACA1ULJLiDZIr6vd2Yt14r121VTF1ZhXrYOP7CfDhrUI9mlWe3z5Rv15owlKqusVijkqWfnQp197Ch1Kc5PdmnWKq2o1tNvzdbK9dtVU1evwrwcjTuwv742YaRCnpfs8qwUBIGmzlqitz5epNKKaoU8Vz27dtCFJ4/RkL7dkl2etTZuLdMTr3+ileu3qaY2rIK8bB1+0ACde8IhysnKTHZ5VvL9QG9+vFBvfbxQpZXVCnmeenYp0kUnj9WIQT2TXZ61Nm0v15OvRde16l3H0AE65/hDlJ2VkezyrGSM0ftzlumtGYsb1jVXPbsU6YKTxmhgr87JLs9aFTtr9Ny787Ryw+7jtf464bBh8jzm4JuybO02vT1zafQYGnLVs3ORzjhqpDoU5CS7tJilbWjZUV6lh1+eqYde+kTbyqrkuY4cx1EQGAXGaPSQnvrG2Ufo/OMOYmNoYIzR81M/0z3Pfqg5i9fJcx25riNJivhGv7n7ZZ19zEH6zkVHa9SQXkmu1h7L127VHY+/o6ffmqOIHyjkOTIm+t5dT01Tl+J8XXvOBN1w0THKz81ObrGWCEd83ffs+7rn6Wlat6VUIc+VMZLjSEbSXx95QxMOHqwfXH6iThh3QLLLtcbcJet0+6Nv6LUPPpfjSHIcmcDIdR1Nfvlj/fKOZ/X1M8frR18/WR2L8pJdrhXqwxHd/dRU3fPMe9q8vUIhz1VgjKLtc/TPx9/R2JH99cPLT9LXJhyY7HKt8dnS9br90Tf0ygfzJUmu4ygwRq7j6JEXp+vmO57VlWcdqR9cfqKKC1nXpOgk6UMvTNe/n35PqzfuUGi3sYUx0v899IaOHTNU37/seB0zZmgSK7XLqo07NOnJqfrPG7NVVx9RKOTJNBxE73j8XfXoXKRvnHukrjtvgnKzmZSRouO116cv1r+f+1DTP1st1/livOb7gW6e9KIuPOFgfev8IzW8f+pMAKZlaFm6dpsu/dXD2lK6U0EQXfH9wCg6HIqav3yTvv/XZ/XSB5/rrp9flPYbQjji66a/Pasn3pjzxYofmIa+Rfm+0YvvzdcL783XP392kS444eAkVWuPabOX6opfPaD6iC/fDyRJ4YjZY5ltpTv150fe0AvT5uk/f75e3TsVJqNUa+ysrtPVv3lA02Yv2fWzSEPvdvfxZyv04dzl+t8bztb3LjuhPUu00kvvzdM3f/vIromX6O4suq4FfvSfldV1uufpaXrl/c/0zO3fUf+e6T2rW1lVq8tvvlfTP1u5axC057oW/dmcRWv09V/ep19ce5p+cuUpchwnCdXa45X3P9M3bnlYgTG7jqGB9vxnRVWt7vrPu3rpvc/07O3fVt8enZJWrw2qa+v1jVse0dufLNr1s73t1z74dLnem71Uv7vxHN1w4THtWaKVpn+2Upf/8kHV1Id3O4b6eyyzaXu5/nj/a3ph2md64k/fUOcO6X3Fh+8H+s3dL+uBF2bIaxivBcbsOg5IUn3Y15Nvfaqn3p6rf//yEp0+YUSyyo1J2p1CWLelTOf97AFt3S2w7E3QcAB7a+Yyfeu2/yji+00uu78zxuhn/3hOT745R5L22beIH8j3A91425N65cMF7VWilWYtXKPLfnGf6uoju3a2TQkCo6Vrt+rCm+5WZVVtO1Von0jE17X/+6Dem7NUxmjXGam9aQzMv737Bd3/3PvtVKGd3p25WN/434cU8X35wb7XNT8wWr+1TOf9cJK2l+1spwrtUx+O6Ipf3a8Z878ILE1p3Of96YFXddd/prZDdfZ6b/ZSXfM/D0bXtWb2a35gtG5Lic794STtSON1zfcDffN3j+qdmYtbsF+Lxr7fTPqvHnlxervVaKN5S9frkp/fp+q6+uaPocZowYpNuuhn96qqpr6dKrSPMUb/8+9X9OALMyRpj4nlL/P9QBHf13W3Pq53Zy1rrxJbJe1Cyw9vf07lO2v3+Re5uyAwemfWMj300swEV2avVz9aqCmvzd7njnZvvnPbk2k7AA+CQN+69VH5QbArADfH9wMtX7tV/++h1xNcnb0mv/Kx3v1k8T6D8d7cfMez2rC1NEFV2a0+HNENtz7a7GBod74faOO2cv327hcSW5zFHvrvh/pw7vIWHwsa/e9dL2j1xu0Jqspu4Yivb936SMzr2oZtZfrdv19MbHEWe/y1mXpj+sKY92s/+/sz2ry9PEFV2c0Yo+t/P0XhiN/ivvlBoEWrNuv2yW8luDp7TZuzXPf/92O1dE1r3I6/9YcnVFVrf9iLKbTcdtttOuyww1RQUKCuXbvq3HPP1ZIlS5r/oCWWrN2qj+avbnYm8iuMdN9/P252Nm5/dd/zH+06xdhSRlJtXVhPvfVpYoqy3Dszl2jd5tKYD1J+YDT55Y9VVVOXoMrsZYzRPU+/J8Vx5Y0jpe2s5EvvfaYd5VUtDseN/CDQ02/OVmlFVYIqs5cxRvc8857U4kP7F1zH0cMvfNT2RaWAV96fr+2lO2Nf1/xAT705W2WV1QmqzF6N65obxyWFRtKjL3/c9kWlgA/mrtDKDdtjnlQIAqOHX/xYtfXhBFVmt/ueny7Pje18RGCMKqvr9N+pnyWoqrYT059s2rRpuvHGG/Xxxx/rzTffVCQS0SmnnKKqqtQ46D3y8syY/zKl6I5jzeZSfTBvVdsXZbkV67fpw7krY95xNLr3uQ/TMuzd/9yHca1rklRVU69n306/sDdj/iotXbMl5jN6UjTsPfj8h6oPR9q+MMvd++x7u+4zi1XEDzTl1U/auCL7vT9nmVZv3BHnuhbo4Rc+Um1d+g2K7nnmvZgnsBqFI74eT8N1bfaitVq0anPMQU+KDsAfeP6jr9zDkQ7uf/6juB+CVFFVmxID8La2dnOp3v5kaewT84pOxtz7/HTrx2sxrRGvvfaarr76ao0cOVKjR4/Wgw8+qLVr12r27NmJqq9NvRnnX6YkhTxX785OjWv+2tK02csV7z2nRtLKDTu0aXtFm9ZkO2OMps2Of11zHUfTZi9t46rsN3XWkj2ephOrkooqLVq1qQ0rsl9dfUQzP18d8xm9RsYYTZ2ZOmfL28q02a1b1yqqavX58g1tWJH9whFfM+bHP4EVXdcWt3FV9ps2a2mrnkC6vWynlq7Z0oYVpYZ3Zy5p9j6Wpniuo2lpOF77YO6KOM4dRwXGaNGqLSqpsPtsaKvuaSkvj15r2bFjxyaXqaurU0VFxR6vZCnf2br7K8or0+/+jLKdNXGfMdj9/5FOausjrZoZC4yxfseRCOWV1a1+KlN5ZXqta+VtsG2VlKffDdJlbbCepNulTm2yrqXhfq2ssjquS8P2/H+k134t4vuqacWZTD8wabd9SlJZZW3cZ0IbtcV2nkhxj0aNMfrxj3+so446Sgce2PSz62+77TYVFRXtevXp0yfeX9lqGaHWDb4zM9LvCwCzMkJxXUKx5/8jvfqWGWr9nzc7M/2eRp7ZBn/mzIz06ltb7JPSrWdSdBttbUBOt761xX4t3Xomtc1+LSvNjgee68Z9yasUvccx3XomRY8HrR2v2b6Nxj2K/+53v6vPPvtMjz/++D6Xu/nmm1VeXr7rtW7dunh/Zav17V4c94xHYIx6d+3QtgWlgN7dOsR9mZMkeZ6rrsUFbViR/TzPVZfi+J8TH/Jc9e5W3IYVpYbeXYv3+r0FseiVZttoYV52q75DyvNc9UvD78/o1a24Vfu1xv9HOinIy1Z+blbcn4+ua01flbG/aov9Ws8uRW1UTWpwHEc9Osf/Z/Y8Ny3Ha326dYjr3qlGGSFPnSz/0uG4Qsv3vvc9vfDCC3r33XfVu3fvfS6blZWlwsLCPV7JcvmpY+L+C3UkXXDCqLYtKAWccsQBcR+oPM/VWUcfqIK89PuW9yvOPKJVN0dPPG1cG1dkv/NOPDTuU9ue62j86EHq0z29BkWu62riaYfHfc287weaePrhbVyV/S44aUw8D6mTJLmuo7Ej+mlQ7y5tWpPtHMfR5ae3bl277LT0W9fOOW503PdPea6jow8dop5dOrRtUSng66ePa9Ux9JJTx7ZxRfY7fuwQFRfmxvVZz3N1wQmjlZOV0cZVta2YtiRjjL773e/q2Wef1TvvvKMBAwYkqq6EOO+4g5QXx6yk57o6bfwBaflN5TlZGbri9HFxDSZ9P9C154xPQFX2u/LMI+I6Tes6jg4a0ksHD0veZZTJ0rlDvs494ZC4BkV+YPTN849OQFX2u+bcCXHfsNq/ZycdfeiQNq7Ifj06F+n0o0fFta4FgdF156fnN5VffU5865ojaUCvzjrqkMFtX5TlOhbl6YIT49+vXXfehARUZb+vnz4urokF13E05oC+OnBQzzavyXaZGSFdfWZ8Yc/3A119lv2TCjFtRTfeeKMmT56sKVOmqKCgQJs3b9bmzZtVU2P3jTuNcrMz9f1LYhvYOJJcV/pejJ/bn3zj3PHKycqMaUPwXFeHH9hf40b2S2Bl9urVtViXnjo25ssRA2N005UnJ6gq+914yfFypJieWOd5rgb16aLTjjooYXXZbFj/7vrahAPjmli46apTW31vR6r63mUntvwbEhs0XuJ01rGjE1SV3Yb07aYzjh4V86DISPppGq9r3774OLmOE/N+bVi/bjpl/IjEFWaxbp0K9fUzDo95nQmM0U1XnJSgqux39ZmHqzA3O6bjgec6On7sEB08tFcCK2sbMYWWu+76/+3deXhU5d3/8c+ZJZM9IZCEBEgIgqyyuwRxBWlRqVarWFyL1qq0lfK0P1u1rY9FaWtr+7QiLa31sSqFpy64VAXcAEVljSDKJgjITsi+z8z5/RESCSQkZyCZe5j367rmUibDNV++mTPn/pz73OfMUklJiS688EJlZWU1PubNm9de9Z10P7j2PE28ZGibXlv/JWPp8f/3LQ3uHX2pvUGPzE56+lc3yeNu2+I4t8tSz+w0PfXfN0btTkqSfvujq3XmoJ6Ogsv9371Ul0bp4FuSzujTXX+5/0ZJbdvBu90upSXH69+P3CHvSVgoHKlm3X+D+uVlOdpRTbnuIl339eg7DbHB8P45+tM9325zSHa7XUpOiNW/f3dnVC7ybTDz3us1sFe2o8/aDyeN0bVfO7MdqzLbgF5Zmv3zG2TJatP+wO12qUtqoub99rvyuKP3e+2hKd/Q6KG9HO1DH7zzco05u187VmW2zM5JmvPQTYrxetq0jbpdlvrmZmr2vRM7oLoT5/j0sOYet9xySzuVd/JZlqVH775C0yZdIK/H1ezG0PCL7pQcp2d/dYMuHz2wo8s0zqghvTT/97crPbV+gXlzG0PD9PeoIb302p/uCvncylNFbIxX/37ke5pwQf1aqJZOD7AsS16PW7+derXunjSmI0s00pUXD9MzD9+mOF9MiwPKhstw9+6RrkV//S/lROFi8iMlxcfqlT/9QKOH1Z/q1dJnze2qPxBz722X6oE7vtGRJRpp4tfP0j8e/I58Xq+sFo6EN3zX5WZ11qK/TlOvKFvLcrTEeJ9e/tMPdP6I0yWpxfUaLpcll8vSz2+/XD+//fKOLNFIl58/WM/OuFXxcfWnqDf7WTvcy765mVo46251y4iuiz0cLcbr0ZyHb9U1lwyXpGZvv9Cwj/DFePTHH1+jO74VnaduHml4vx569Q+3N17Aobm+NXzWxpzVVy8/+t2IWXsclYeLXC6XfnLDxZo84WzNXbhG//uf5fpyf8nhn1ka0b+HbvvGOfp6fr+oPnp7tBH9c7Tq2Xu04IPP9MRLH2jZx1sbfxbjdetbY4bplgnnaEgETDF2lDifV3//5U36r5v26MmXlmnuGyuaXH++W3qqbrtqtCaNP0tphl+1oyN9bdRArX/hQf174Ur99fnF+nzngcafWZIuOquvbvvmebr4rH5yneB9hE4VyYlxeu73d2r5J9v0xIvv6aV3C5qsP+iUHK/vXHGubpowKiqvTteSCRcM0XnD+2juG8s1+/kl2rHnUJOfjxraW9+96jyNyx8oD/sDSfVXEvu/R+7Qqk+364kX39OLb69ucoWszikJ+s6Vo3XT5fnKjsKrOLVk7Nn9te65X+qFN1dr9vNLtfGIm0ZalqVLzu6vW795rs4f3ofvtcN8MR79+Z6JunvSRXrqlQ/1zGvLVVFV2/jz3OzO+u43z9W140YoJTEujJWaZeBpWfrwyWl6Z9Vm/eOlD/XOqi2Nd7uPj43RdeOG66bLzlS/nplhrtQZy7ZP9KrOzpSWliolJUU7d+4M65XEjuYPBFRbF1Ccz2vcKU1+v19Lly7VeeedJ4/HnJxZVVOn4rJKeT1upSbFGTeNbWLfqmvrdOBQmapq6pScEKuMtCSjdk4m9sy2bX25r0hFpfWftYy0JHVODf2S0u3BxL6VlFVq1/5iVdXUKjkhTrnZnY26Br+JPQsGg9q2u1DFpRXyeNzK7Jysrp3NutysiX0rLa/S7gPFqqyuU3JirHKzOht1wM/Entm2rb0HS1VUVimvx6WMtGTjBt0m9q3OH9ChkgpVHt6HpiXHGzVmM7FnUv0Yt7SiRl63S4nxPqN6JkmHDh1SXl6eSkpKjpsNzOlomHncbuMG3aaL83kV5zNrh2662Bhv1F2W90RZlqUeXdPom0MpSfFKSYruUzSdcrlchy9lHN2ngDmVnBinZMMG3KazLEtZ6SnKirJ7sJwo7+GDCXDG43Yr7RQ4Zd+cQ7wAAAAA0AxCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIxGaAEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7QAAAAAMBqhBQAAAIDRCC0AAAAAjEZoAQAAAGA0QgsAAAAAoxFaAAAAABiN0AIAAADAaIQWAAAAAEYjtAAAAAAwGqEFAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGM0Trjf2+/3y+/3hevuI0tAn+uUMfXOOnoWGvjlHz0JD35yjZ6Ghb87Rs9AEAoE2vc6ybdtu51okSTNnztTMmTMVCAS0adMmzZkzR/Hx8R3x1gAAAAAMVFlZqUmTJqmkpETJycktvq7DQkuD0tJSpaSkaNu2bcctDF/x+/364IMPlJ+fL48nbJNjEYe+OUfPQkPfnKNnoaFvztGz0NA35+hZaIqKitS7d+9WQ0vYOurxePiFOkTPQkPfnKNnoaFvztGz0NA35+hZaOibc/TMGbfb3abXsRAfAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIxGaAEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7QAAAAAMBqhBQAAAIDRCC0AAAAAjEZoAQAAAGA0QgsAAAAAoxFaAAAAABiN0AIAAADAaIQWAAAAAEYjtAAAAAAwGqEFAAAAgNEILQAAAACMRmgBAAAAYDRPuAsAAACRrbYuoKragBLjPHK7OB7aGtu2tedQhYrKq+V1u5TZKUEpCb5wl2W8Lw+UauvuIlXV1Ck5wadBeRlKiqdvx1NcXq33136horIqeVwuZXVJ1rmDc+VxR952SmhByCprAyqp8svlktLivfJG4AbQ0faWVGvjvgrV1AWVFOvRkB7Jio9xh7sso315qEpvf3pAhyrq5HFZ6tE5TmMHpivBx9dXS3YdqtILy3dp4+4yVdYGlBrv1ai+nTV+aFf61oJ9xVV6/oMden31bh2qqJXbZSmnS4KuGZWjrw3NVizb6TGKymv04rJtmvPuZu0qrJAkWZJGnp6uGy46XRcN6cZ+4ShVNX698uHn+uei9dqw81Dj8y7L0tjhubpx7ACd0z9LlmWFsUqzBIJBvbP6Cz31RoGWrt3R5GexMR5dfUF/3fi1IeqX0yVMFZrpk6179b//WaUXFn+qOn+gyc/SUxN086XDdf24oeqSmhCmCp2L+r3XloOVWrCpUOv2Vqg2EFSSz6PRPVM0pnea0uK94S7POLZt6+NdZXpl3QF9uL1Ytl3/fLzXpa8PSNdlA9OVlcJRj6Ot3l6i/1u5W8u/KG7yvM/j0tcHpevakdnKSokNT3GG+mx3mWa/84U++rxILqt+py5JgaCtP77xuS4bmqnbL+qp5Di20wb7S6r10Isb9P7GQlmWFDy8fVqSFn92UH/4z2ZdN6qH7rykV0QeZWsPtf6AHnruE83/aKekr3omSSUVtfr4iyL9+oX1uuebA3Xl2T3CVKV55i3ZoofmrpY/GGzcD0iSLWnV5oNasemAunaK0+NTzlf/Hp3CVqdJ1m8/qMm/e0OFpdU6OpMEbVtvrdmuhau+0KgB2Zr5w7FKiosJT6EGKa+q1Z2/f1VL1+6Q23VskKuu9WvuW5/o2UXr9NPrR+t73xhB4JP05Ksr9cu/vymX21IgYB/z8wPFFXr0X+/piVdW6plfTtTg3l3DUKVzUbvXqqoL6OG3t+lnb3yuJduKVVhZp7KagHaX1ui5dft154sbNH/9Adn2sb/saOUP2PrDO9t17yubtfyIwCJJlXVBzV+7T7f/6xO9s6kwfEUaaM5Hu/Tj5z7Vyu3Fx/ysxh/Uqx/v03f/uVaf7Crr+OIM9f6mQn3vyQKt2FokqX4g6Q/a8gdt2arv2/xVe3Tr39foQGlNeIs1xJeHqnTDYyv0weZDstV08N3wv9V1QT21eLv+6+m1qgsEw1GmUer8Qd311+V68aOdCtpNeyZ99efyar9+/q+P9dQ7Wzu+SAP9Y+EGPfDsStUFmgaWBsHDTx4oqdak37yp9dsPHfuiKPPp9kJdN/1VFZfXf18117fA4Q/chxv26MZfv6aqGn9Hlmicmjq/bnl4vt7/pP6AQuDoDfSwhud//ex7evzFFR1Wn6meem21fvH3N2VLzQaWBkHbVmlFta69/1lt3HGg4wo8AVEZWmoDQU1/a5s+3lMuqfkdVdCWnl2zVy98Ehm/yPZm27YeW7Jdb22sDyTNbQdBu/75R976Qsu2FXdsgYZ6de0+/f29+unsFr5vFbDrQ/Q9z3+qHYVVHVidmTbtLde9//5UgYDdYs+k+n7uLa7Wj55dp1p/dA/Aq+sCuuuJ1SqurGtxx97AlvT+xkL94T+bO6Y4g/12/not31LY7ACyOb976VMt2xDd+4SPNu7TI88XtOm1gaCt2kBQt/95sSqjeABe6w/otkcXqNYfaHX7lKRg0Nan2wv10JwPO6A6c/3h/z7U6k17FGxDzxo8MneZVm7Y3Y5VmW3Ll4X6+d8Wtvn1gaCt6lq/vjvjhYg4SO84tCxZskQTJkxQdna2LMvS/Pnz26Gs9vX6hkJtPlh13AFRg7kf79Oukur2L8pwG/ZVaOGGQrX1I/3Y4u3yHyfhR4OauoD+unh7m15r21KtP6inPtjZzlWZ7x9Ltito2236rAVsaeuBSr31aXQPJBet3a9dh6rbNCCS6oPLcx/tUmFZ9M5SFZXX6rkPdrQ5sEiSy5L+/uaW9isqAjyxYEOzp+m0JBi0daisRq8ub9t34alo4covtL+4ss3bp1R/FPz5pZtUVB6d44+qmjo9s3Bt46xdW7ldlp58fU07VWW+f76+uvFU6rYKBG1t21Ok99eav406Di0VFRUaMmSIHnvssfaop90FgrZe39j2wbfLkhZsYmr71fUH5HawHRRX+fXhUes3os27mwpVURto/YWHBWxpyaZCHaqobceqzLa/tEbvbSyUkzOXXJb03PJd7VdUBJi7bOcx58i3JmjbemnlnvYpKAK8+NEOR4NIqX52b8WWQm3dV95OVZnty4Plem/9Hsd9syzpmbc3RcSR3Pbw1KL1cjkIeg38waCeX7qpHSoy3yvvb1J5lfN9YSBo6/WPtmh/UUU7VGW2iqpazXtzrePtU5LcLpee/M+qdqjq5HIcWsaPH6/p06frqquuao962t3mg5UqrKxr8+uDtrT48Hn10SoQtLV0S1Gzp4S1xGVJ726J7rD3zoZCxwPJgC29vyV6P29LN7b9VJ0GQVv6bHe59kfp2pYDpTXasLvMcd9sW1q0dl/7FBUB3vx4r+OeSfXfbe+s23vyC4oAi9ftqb+qg0O2LW3eXaK9RZUnvyjDlVXVas2W/Y5OcWpg29Kbq3e0/sJT0LsFXyiEnCepfnbvvXXR17dVG3epsqbt49sjBYJBvbvG/DV77b6mpaamRqWlpU0e4VRU5fy82sq6YFQvWq2sDcgfwhHJQxWhbTynisKKWseDIrfLUnFV9PatuLLO0aknR//daFR8AjNzRVHaM0kqqggt5LpclkqitG/FFTUhb59S/dXYok1J+YkdTCkqi87Tww6VVrbpFP7muCxLJVF4Wt2J/ptr6wKqrjV77Vm7h5YZM2YoJSWl8dGjR3gvGekJ8Qv3RL6oI53XyXlhR4gJ8e+dKkLpm23bIff7VOB1W20+dfNoMVF6CV+vJ/R/dzRvoydy/5AT6Xkki/G4QpqdahCNfYvxntj9fXze6OuZJPliQr8jR9C2T7jvkcjrOfF/s+n3VWr36n72s5+ppKSk8bFzZ3gXGuelxTqa3bYk9ewU63hh06nE53EpI8nZ9eJdlnRaenw7VRQZ+mQmOloHJNXPUOV1id6+5aUnhHQ+bozHpcwovT9QZkqsYkIYDLpdUq/MxHaoKDL06poU0sEof8BWz/TIuRnbyZSbkRTS9ilJHrelrp2i77utU2Ks4kIcgLtdlnp2TTnJFUWGnIyUEzpY3CMj+vqWk3li/+bMtES5oz20+Hw+JScnN3mEU5eEGA3vltTmcyVtSeP7dm7XmkxnWZa+MSjDUdgL2tL4AentVlMkmDA409E6IElKT4zRyNzUdqknEuT3SVOnBGc3i3S7LF06JFNxUXrH8rgYty4b1tXxDj4QlK4+u1s7VWW+a/JzQhqAJ/g8Gjskqx0qMt+FZ2QrJcH5DQ/dLkuXnZmrhNjouxGs1+PSt84/PaQBeCBoa+KF/dqhKvNNvHhQyAE5s1OCzj0j+m4E279nhvr3TA9pLZDLsnTD14ae9JpONrMjVTu5cmB6m6a4XZbUOd6rc3umtntNphvbr7Niva42LSx3WdKZOcnqnhrdd3jvnZGgAdltn22xJF01PCuqT0X0uCx968xsRxcwCARtXTUyOgeRDa45p7ujHbzLkjJSfBp1evQekMnvm67stDhHnzW3y9JV5/SI2oAc43XruvN7O74SViBoa9KFfdqpKvNNuri/8yuuSeqRnqT8/tntU5ThBvXK0KC8DMcDcJdl6aavDZHbFX3DW8uyNPmykaGtBbKk6y4ZctJrOtkc/1bLy8tVUFCggoICSdK2bdtUUFCgHTsi50oN/TISdGd+N1lq+UIoLktKiHHr/jE95YvC83CPlhzr0S/H95bHso77JeKypKxkn/5rTF7HFWewn192upLivK1+8VqSzspL1bdGRPfgW5ImjeqhQd2S27yz+t5FPdU7ik9zkqS+2Um65YLcNr3WUv3Obfq1A6M6ILtclh6YOPi4+4EjuV2WMlNiddslvdu7NKPdPLavMlPj2vzZsSzpG2fn6oyeae1cmbn6dOuka84/3fHVJO+ddE5Il0o+Vfz0+tFycrk6t8tSZlqCJl1yRvsVZbgrzh+gvKxOjr7bLUk3jx+urmlJ7VfYSeJ4NL5y5UoNGzZMw4YNkyRNmzZNw4YN0y9+8YuTXlx7uui0NN03pqf6dIlrfK7hd+y2pHN7puo3l/ZW95Toni040uBuSfrNlX0bZ1AaZhAs1ffOknROz1Q9enU/JceGvojuVJKZ7NPjkwYpt3P95+zoWZeGz9zXBqbrwSv6RvUgsoHP49Kj1w/SiMMznM31xH3483bHxT110+joOw2gOVPGnaYbz8uR1PKFQ1xW/fqfR28crBG9OnVkeUbK75uu390yQm63ddxtz2VJmamx+vuUc5SWGJ1rpxp0SvTpyR9dpPSU2DZ9X100uJt+ddNZsqJ4XagkPXjzaI0ZltPq66zD320P3nKuLhnetgMRp6rRg3P0mzvH1veklY+P22UpNTFWT99/lTolxR3/xaewOJ9Xc/77OnVJTWjT9mlJuuSsPvrF5DHtX9xJ4HhkeeGFF54yN4gakpWkIVlJ2l5UpQ37K1UTCCrJ59aIbskMulvQLzNBsyYO0Gf7KvT2xkIVVtbJ47KU2ylO4/p3cbxgPxp0TYnV328aojU7S/VSwV6t31WmGn9QSbEenX96miYM6apuUX4q3dESfB794YYztGpbsZ5bsVvvbfrq/i3xMW5NGNZV3xyZpZzO0bewtyUul6Wpl/bR6H5dNG/ZTr376YEmpwkk+Ny66qxuuuac7uqWFr079aNdMiRL//rRaP3vO1u1YM1u+YO23C5Ltm0raEupCV5de25P3XhBnlJDWM9xKsrNSNK/7/2aHnvlE83/YJtq6wJyuy3Zdv0gyB+0lZkap1vG9tWNY06PylN1jub1uDTzB2P1l1c+1j8WrFNJRa1cltV4x3e3y1IgaKt3difdM/EsXTiEgzGSdM2FA9UlOV6/+Mc72rm/tLFPDRr+fM7A7vrtnZeoW5fwrps2QfeMFP3nd7foh4++rGWf7DimZ1J9CPS4XJp8+Uj97KYLjV+A34CRuaTcTnHK7cROvK0sy9KAroka0DW6T8lxwrIsDc9J0fCc6LuiSahclqUze3XSmb06qao2oJKq+oCcGu+VJ0K+YMNhZK9OGtmrkw6U1mjr/gpV1waUFOfRgO7Jio3Cy4C2Rb/uKfr1jcP0kysHaMn6fTpUXiuP26XuneN0/oDMqLxUb2u6JMfqgetH6sdXDdEry7/Qlt2lqq71KzHOq/z+XXXewK6ElaO4XS5NuWKYbrt0sBas3KaFq7arsKRKMV6XundJ0tXnn65hp2VE/azU0S4anqfFQ3tq2Sc79czCtdq086Cqav1Kjvdp9OAcXX/JYPXKZub4SJlpiZo3fZI27jigp99Yo9eWbVRpRbXcbpeyOifpukuGaOLFg9UpObLGvoQWAMaLi3FH7eLnUKUn+5SeHN2nMjnVOcmnb57T+ik8+EpinFffviB6F9mHwud16xv5vfWN/OheH+WEy2Vp9OAcjR7M9ulE35x0Tb99nKbfPi7cpZwUHAYBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIxGaAEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7QAAAAAMBqhBQAAAIDRCC0AAAAAjEZoAQAAAGA0QgsAAAAAoxFaAAAAABiN0AIAAADAaIQWAAAAAEYjtAAAAAAwGqEFAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBonnC9sd/vl9/vD9fbR5SGPtEvZ+ibc/QsNPTNOXoWGvrmHD0LDX1zjp6FJhAItOl1lm3bdjvXIkmaOXOmZs6cqUAgoE2bNmnOnDmKj4/viLcGAAAAYKDKykpNmjRJJSUlSk5ObvF1HRZaGpSWliolJUXbtm07bmH4it/v1wcffKD8/Hx5PGGbHIs49M05ehYa+uYcPQsNfXOOnoWGvjlHz0JTVFSk3r17txpawtZRj8fDL9QhehYa+uYcPQsNfXOOnoWGvjlHz0JD35yjZ8643e42vY6F+AAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIxGaAEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7QAAAAAMBqhBQAAAIDRCC0AAAAAjEZoAQAAAGA0QgsAAAAAoxFaAAAAABiN0AIAAADAaIQWAAAAAEYjtAAAAAAwGqEFAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBonnAXAAAAEI0CQVuWJbksK9ylRISquoC2F9eoNmArMcat3FSf3C56dzwHKmr1wY5SHaqsk8uylJPq09k9khXndYe7NMcILQAAHFbtD6qiLijLspQS42JA1IqgbWtPeZ12ldWpLmgr3utSXkqMUmMZXrSkqi6oTw5Uat3+SlXUBSVJmQleDcmMV5+0WAJMM0qr/Xpra7E+3luhgP3V84kxLuX3SNa5Oclsq0cprfHrn6v36pN9FTqyNct2SM9/ckCX9E7T5f07R9TnjW8VAMbzB23VBCVLUqybo5KtCdi2CquCOlQTVCBoy+e2lB7nVnKMJYveNavKH9SmolodrAo0Pue2pO6JXvVK9fKZa0ZJTUBLd5ar0h+UJclW/Ta6uahG2Qlend0tQV4Gkk0UVtbpxY1FqvYHdcTYW/sr6rRwa4k+O1ily3p3ktdN3xocqqzT31btVWVdUEG76c/Ka4Na9Hmxviiu1vWDMwguh5XV+PWbxTtUVFUnScf0rS5o67VNhSqqrtNNw7pGzH4hpDUtjz/+uPLy8hQbG6sRI0Zo6dKlJ7suRAi79ZfgCHVBqaxOKqmTKvxqcsQIx7JtWwerbe2qkg7WSAdqpC8rpbI6GteSuqCtjUV12lURUJXfVm2wvl9bS/3aWeaXbdO7o1X5g1qxt0qFRwQWqX773F5Wp7UHqunbUSrrgnpne5mq/PUzBQ3dafjvnoo6LfuynL4docYf1PxmAov0Vd++LK3VW9tKOro0YwVtW09/vL/ZwHKkzYXVentrcYfVZbo5H+9TUVXdcXsmSR/sKNXyL8s6pqiTwHFomTdvnqZOnar77rtPa9as0Xnnnafx48drx44d7VEfDGVLqpNLtXKrRm75FRkpPZxqAlKpX6q1Jb9tqTpYH14ILi07VCtVNB1Hyj78fKWfxjVnZ5lfNYHmf1ZYY+tQTbBjC4oAnxfXyh9s+SBMYXVQ+ytbaGqU2nioWnVBu8We2ZL2Vfp1oNLfkWUZ7bODVapsJrAcyZa0uahaJdX0TZI+P1Stg5X+VgffkvThl2WqDfD9VlRVp4Ld5W3qmSXp7c8PtXtNJ4vj0PLoo4/q1ltv1W233ab+/fvrj3/8o3r06KFZs2a1R30dxnK55Ha7pQiZIgs3v1yqPyGg/hGQWwGCS4ts+8jBt9X4X1sS+/TmBWxb5cfpTWldx9USKWoDtkpqj7+nOlDFTv1IdQFb+ysDrc4a7yznA9cgYNvaVlzTas8sSZ8X13RESRFh3YHKNr3OkvTpwar2LSZCrN5d3uaBam3A1oYD9G2Fg5kTW9L24hrtK69tv4JOIkehpba2VqtWrdK4ceOaPD9u3DgtW7as2b9TU1Oj0tLSJg8Tud1uudzu+uCCVn0VWBrYh59DcwK2ZB/TM0myxJlOzatt5cA2EwbHqm7DtF01M1RNtHbku0HDgmlItX5bbfkY2ZJKa+lbg7KWpkCbUergtaeyomq/2voJcllSSQ1HAYur/HK6tKe4KjL65ii0HDx4UIFAQJmZmU2ez8zM1N69e5v9OzNmzFBKSkrjo0ePHqFX2wE4/9aJpr0isrSs5Qk8m5sltaC1SU8+b8dqy2eJyeSm2rrA3k3jGrkcfGmxnvwrThaJs6C8npMLOdg226kkuUMYVITyd8IhpDKPvsqAbdstXnngZz/7mUpKShofO3fuDOUt252/rk7+ujoFAxzdaAtP47EPWw3hxd3m4yHRx21JHuurXh3Jx+Res3yu4w94Ern24TESvJY8reyzU2PYqR8p0WvJ18rI2pKUHseG2sDndim1jV9c2Ynedq4mcvRIjmnTwRZbUvekmPYuJyL0TI1t8wEqW1Juqq89y4kIvdLiHK2V9bosdUuOjL45Ci1dunSR2+0+ZlZl//79x8y+NPD5fEpOTm7yMBWzLG3nli2vAnLLlltBxbCipVVJnqMH4bZ8Lik2Qo5wdDTLstS5hf2225JS2Kcfw7IsZSW0PJh0ScqMJ+0dybIs5Sa3PrDunsTg+0inp7U+yLEk5TGIbDQ4I75NpyL63JZ6p8W2ez2RYGS3xDa9zpKUleSNmMF3exrcNVFJbTyo4LKkc3Ii50aTjoZLMTExGjFihBYtWtTk+UWLFmnUqFEntTCYz6X6GRfP4dUaOD6XJaV4pBSPrSSPrVRv/WwBs9kti/NY6horxbvrP29uS0r2SllxnAbQki5xbnVLcB/z5R7jknqnehTb2lRMFOqe6FG3w1N3x646k87o4lOCl6MLR8pNjmkx7DX08OzsBMV56FuDbkkxGtCl9TAyJi9FHk4PkySlxHo0Orf1g92WJY3vk9YBFZnP7bJ0Rf8urb7OkuRxWRrbO3L65viQ27Rp03TjjTdq5MiRys/P1+zZs7Vjxw7dcccd7VEfcEqxLLV6+g6a8rktpUfGQSBjZMS71TnOpdKaoAJ2fQ8TvdxYsiWWZalvpxhlxnv0ZXmdymqDcknqEl8fZhh4H8uyLJ2VlaC02BptPFTd5BLkXeI8GtglVhkJzE4dybIsXdQzRXEet9bsq1DQVuOC6aAtxXlcurhnsnp1YpblSGNPS1XAtrVsR5lcVtMbJVqqH6Rfd0YX5dG3RqN7pqqk2q9XNhQe0zOp/nPndVn6fn53ZSZGzmkLjkPLxIkTVVhYqAcffFB79uzRoEGD9Nprryk3N7c96gMAhMBtWeoUS9prK+twv+hZ21mWpT5pserdyaeSmqD8QVtxXksJEXKqSTi4LEujeiRpeFaCNhZWqbg6IJclZSXGqFcnX5svDBFNXJal8X3SNDwrUSt2lWvDwUrVBWwlxrg1PCtRw7ITFM9n7hiX9eui3p3j9fbnh7R2b0XjqYk+j6XRuam6qFcndYmwAwshndx811136a677jrZtQAAgAhjWZZSCXuOxHpcGpKZEO4yIkpmYowu75umy/tGzulM4dY3PV590+NVXuNXcbVfLstSlwSvYiLlcmFHYUUmAAAAcIpK9HmU6Iv8IX9kRi0AAAAAUYPQAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIxGaAEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7QAAAAAMBqhBQAAAIDRCC0AAAAAjEZoAQAAAGA0QgsAAAAAoxFaAAAAABiN0AIAAADAaIQWAAAAAEYjtAAAAAAwGqEFAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0Twd/Ya2bUuSysrKOvqtI5bf71dlZaVKS0vl8XT4ryxi0Tfn6Flo6Jtz9Cw09M05ehYa+uYcPQtNQyZoyAgt6fCONhQ2YMCAjn5rAAAAAAYqKytTSkpKiz+37NZizUkWDAa1e/duJSUlybKsjnzr4yotLVWPHj20c+dOJScnh7ucY5x55plasWJFuMs4Bn1zjp6Fhr45R89CQ9+co2ehoW/O0bPQmNw327Y1YsQIbdq0SS5XyytXOnymxeVyqXv37h39tm2WnJxs3C9Tktxut5F1NaBvztGz0NA35+hZaOibc/QsNPTNOXoWGlP7FhMTc9zAIrEQP2JMmTIl3CVEJPrmHD0LDX1zjp6Fhr45R89CQ9+co2ehaUvfOvz0MFOVlpYqJSVFJSUlRiZQU9E35+hZaOibc/QsNPTNOXoWGvrmHD0LzanQN2ZaDvP5fPrlL38pn88X7lIiCn1zjp6Fhr45R89CQ9+co2ehoW/O0bPQnAp9Y6YFAAAAgNGYaQEAAABgNEILAAAAAKMRWgAAAAAYjdACAAAAwGiEFgAAAABGI7Qc9vjjjysvL0+xsbEaMWKEli5dGu6SjLZkyRJNmDBB2dnZsixL8+fPD3dJxpsxY4bOPPNMJSUlKSMjQ1deeaU2btwY7rKMNmvWLA0ePLjxDr75+fl6/fXXw11WxJkxY4Ysy9LUqVPDXYqxHnjgAVmW1eTRtWvXcJcVEXbt2qUbbrhBnTt3Vnx8vIYOHapVq1aFuyyj9ezZ85jPm2VZ3JjwOPx+v+6//37l5eUpLi5OvXr10oMPPqhgMBju0oxWVlamqVOnKjc3V3FxcRo1apRWrFgR7rJCQmiRNG/ePE2dOlX33Xef1qxZo/POO0/jx4/Xjh07wl2asSoqKjRkyBA99thj4S4lYixevFhTpkzRhx9+qEWLFsnv92vcuHGqqKgId2nG6t69u379619r5cqVWrlypS6++GJdccUVWr9+fbhLixgrVqzQ7NmzNXjw4HCXYryBAwdqz549jY9169aFuyTjFRUV6dxzz5XX69Xrr7+uTz/9VL///e+Vmpoa7tKMtmLFiiaftUWLFkmSrrnmmjBXZq7f/OY3+stf/qLHHntMn332mX7729/qkUce0Z///Odwl2a02267TYsWLdLTTz+tdevWady4cRo7dqx27doV7tIc4z4tks4++2wNHz5cs2bNanyuf//+uvLKKzVjxowwVhYZLMvSiy++qCuvvDLcpUSUAwcOKCMjQ4sXL9b5558f7nIiRlpamh555BHdeuut4S7FeOXl5Ro+fLgef/xxTZ8+XUOHDtUf//jHcJdlpAceeEDz589XQUFBuEuJKD/96U/1/vvvc3bCCZo6dapeffVVbd68WZZlhbscI11++eXKzMzUE0880fjc1Vdfrfj4eD399NNhrMxcVVVVSkpK0ksvvaTLLrus8fmhQ4fq8ssv1/Tp08NYnXNRP9NSW1urVatWady4cU2eHzdunJYtWxamqhANSkpKJNUPwtG6QCCguXPnqqKiQvn5+eEuJyJMmTJFl112mcaOHRvuUiLC5s2blZ2drby8PF133XXaunVruEsy3ssvv6yRI0fqmmuuUUZGhoYNG6a//e1v4S4rotTW1uqZZ57R5MmTCSzHMXr0aL311lvatGmTJOnjjz/We++9p0svvTTMlZnL7/crEAgoNja2yfNxcXF67733wlRV6DzhLiDcDh48qEAgoMzMzCbPZ2Zmau/evWGqCqc627Y1bdo0jR49WoMGDQp3OUZbt26d8vPzVV1drcTERL344osaMGBAuMsy3ty5c7V69eqIPXe5o5199tn65z//qdNPP1379u3T9OnTNWrUKK1fv16dO3cOd3nG2rp1q2bNmqVp06bp3nvv1fLly/XDH/5QPp9PN910U7jLiwjz589XcXGxbrnllnCXYrR77rlHJSUl6tevn9xutwKBgB566CF9+9vfDndpxkpKSlJ+fr5+9atfqX///srMzNS//vUvffTRR+rTp0+4y3Ms6kNLg6OPbti2zREPtJvvf//7Wrt2bUQe6ehoffv2VUFBgYqLi/X888/r5ptv1uLFiwkux7Fz507dfffdWrhw4TFH2NC88ePHN/7/GWecofz8fJ122ml66qmnNG3atDBWZrZgMKiRI0fq4YcfliQNGzZM69ev16xZswgtbfTEE09o/Pjxys7ODncpRps3b56eeeYZzZkzRwMHDlRBQYGmTp2q7Oxs3XzzzeEuz1hPP/20Jk+erG7dusntdmv48OGaNGmSVq9eHe7SHIv60NKlSxe53e5jZlX2799/zOwLcDL84Ac/0Msvv6wlS5aoe/fu4S7HeDExMerdu7ckaeTIkVqxYoX+53/+R3/961/DXJm5Vq1apf3792vEiBGNzwUCAS1ZskSPPfaYampq5Ha7w1ih+RISEnTGGWdo8+bN4S7FaFlZWcccQOjfv7+ef/75MFUUWbZv364333xTL7zwQrhLMd5PfvIT/fSnP9V1110nqf7gwvbt2zVjxgxCy3GcdtppWrx4sSoqKlRaWqqsrCxNnDhReXl54S7Nsahf0xITE6MRI0Y0XrmjwaJFizRq1KgwVYVTkW3b+v73v68XXnhBb7/9dkR+YZjAtm3V1NSEuwyjjRkzRuvWrVNBQUHjY+TIkbr++utVUFBAYGmDmpoaffbZZ8rKygp3KUY799xzj7l0+6ZNm5SbmxumiiLLk08+qYyMjCaLpNG8yspKuVxNh61ut5tLHrdRQkKCsrKyVFRUpAULFuiKK64Id0mORf1MiyRNmzZNN954o0aOHKn8/HzNnj1bO3bs0B133BHu0oxVXl6uLVu2NP5527ZtKigoUFpamnJycsJYmbmmTJmiOXPm6KWXXlJSUlLj7F5KSori4uLCXJ2Z7r33Xo0fP149evRQWVmZ5s6dq3fffVdvvPFGuEszWlJS0jFrpRISEtS5c2fWULXgxz/+sSZMmKCcnBzt379f06dPV2lpKUdwW/GjH/1Io0aN0sMPP6xrr71Wy5cv1+zZszV79uxwl2a8YDCoJ598UjfffLM8HoZjrZkwYYIeeugh5eTkaODAgVqzZo0effRRTZ48OdylGW3BggWybVt9+/bVli1b9JOf/ER9+/bVd77znXCX5pwN27Zte+bMmXZubq4dExNjDx8+3F68eHG4SzLaO++8Y0s65nHzzTeHuzRjNdcvSfaTTz4Z7tKMNXny5MbtMj093R4zZoy9cOHCcJcVkS644AL77rvvDncZxpo4caKdlZVle71eOzs7277qqqvs9evXh7usiPDKK6/YgwYNsn0+n92vXz979uzZ4S4pIixYsMCWZG/cuDHcpUSE0tJS++6777ZzcnLs2NhYu1evXvZ9991n19TUhLs0o82bN8/u1auXHRMTY3ft2tWeMmWKXVxcHO6yQsJ9WgAAAAAYLerXtAAAAAAwG6EFAAAAgNEILQAAAACMRmgBAAAAYDRCCwAAAACjEVoAAAAAGI3QAgAAAMBohBYAAAAARiO0AAAAADAaoQUAAACA0QgtAAAAAIz2/wG+DFHV1ToQMAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create 100 cells (10*10 grid)\n",
    "N = 100\n",
    "n = int(N**0.5)\n",
    "\n",
    "# create a vector of increasing phases to check the order of the oscillators on the plotted grid\n",
    "l_x0 = np.linspace(0, 2 * np.pi, N)\n",
    "print(l_x0.shape)\n",
    "# reshape l_x0 to match the grid\n",
    "X = np.reshape(l_x0, (n, n))\n",
    "print(X.shape)\n",
    "# plot the oscillators\n",
    "datagrid = np.sin(X)\n",
    "fig, ax = plt.subplots(figsize=(X.shape[0], X.shape[1]))\n",
    "heatmap(datagrid, marker=\".\", ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-hIYlbyKfkWn"
   },
   "source": [
    "#### Question 1: \n",
    "Assume that every neuron is connected to its nearest neighbours in the grid: for example, a neuron in position $(i,j)$ will be connected to the neurons at $(i-1,j)$, $(i+1,j)$, $(i, j-1)$ and $(i, j+1)$. Assume also that the neuron is connected to its second nearest neighbours: $(i-1,j-1)$, $(i-1,j+1)$, $(i+1,j-1)$, $(i+1,j+1)$.\n",
    "\n",
    "Visualize the circadian oscillations of neurons at different times using the provided function heatmap and using subplots.\n",
    "\n",
    "**Hint** : \n",
    "Initialize the parameters of a population of $N$ neurons ($N$~100--1000 cells) with varying intrinsic periods of about 24 hours. For instance, you can modify the period of a neuron $i$ by drawing $f_i$ from a normal distribution with $\\mu =2\\pi/24$ and $\\sigma$ equal to $5\\%$ of the mean (you may have to play with $\\sigma$ to obtain synchronization). Assign coordinates to each cell in a 2D grid, e.g. with x and y coordinates in a rectangle (the real SCN has an ‘egg’ shape). Choose $K = 0.03$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(\n",
    "    l_theta, t, l_f, K, A\n",
    "):  # l_theta: list of initial phase of the N neurons, t = time (needed for odeint),\n",
    "    # l_f: intrinsic frequencies of the neurons, K: coupling strength, A: interaction matrix\n",
    "\n",
    "    # compute interaction matrix\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining some interaction matrices\n",
    "\n",
    "\n",
    "# create a matrix of random connections,\n",
    "# each neuron is connected to another with probability p\n",
    "def A_alltoall(N, p=1):\n",
    "    ### write function\n",
    "\n",
    "\n",
    "def A_NN(N):\n",
    "    # Nearest Neighbours\n",
    "\n",
    "    A = np.zeros((N, N))\n",
    "    n = int(np.sqrt(N))\n",
    "    neuron_index = np.arange(N)\n",
    "    NI = np.reshape(neuron_index, (n, n))\n",
    "\n",
    "    for i in range(n):  # row number\n",
    "        for j in range(n):  # column number\n",
    "            if i < n - 1:\n",
    "                A[NI[i, j]][NI[i + 1, j]] = 1\n",
    "\n",
    "            if i > 0:\n",
    "                A[NI[i, j]][NI[i - 1, j]] = 1\n",
    "\n",
    "            if j < n - 1:\n",
    "                A[NI[i, j]][NI[i, j + 1]] = 1\n",
    "\n",
    "            if j > 0:\n",
    "                A[NI[i, j]][NI[i, j - 1]] = 1\n",
    "\n",
    "            A[NI[i, j]][NI[i, j]] = 1\n",
    "\n",
    "    return A\n",
    "\n",
    "\n",
    "def A_2NN(N):\n",
    "    # First and second Nearest Neighbours\n",
    "\n",
    "    A = A_NN(N)\n",
    "    n = int(np.sqrt(N))\n",
    "    neuron_index = np.arange(N)\n",
    "    NI = np.reshape(neuron_index, (n, n))\n",
    "\n",
    "    for i in range(n):  # row number\n",
    "        for j in range(n):  # column number\n",
    "            # second nearest neightbours\n",
    "            if i < n - 1 and j < n - 1:\n",
    "                A[NI[i, j]][NI[i + 1, j + 1]] = 1\n",
    "\n",
    "            if i < n - 1 and j > 0:\n",
    "                A[NI[i, j]][NI[i + 1, j - 1]] = 1\n",
    "            if i > 0 and j < n - 1:\n",
    "                A[NI[i, j]][NI[i - 1, j + 1]] = 1\n",
    "\n",
    "            if i > 0 and j > 0:\n",
    "                A[NI[i, j]][NI[i - 1, j - 1]] = 1\n",
    "\n",
    "            A[NI[i, j]][NI[i, j]] = 1\n",
    "\n",
    "    return A\n",
    "\n",
    "\n",
    "# plot ieteraction matrices with plt.imshow() to see how they differ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a function for that takes the mean of exp(j*theta) for all theta in l_theta\n",
    "def order(l_angles):\n",
    "    ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a vector of intrinsic frequencies\n",
    "mu = 2 * np.pi / 24\n",
    "sigma = 0.1 * mu\n",
    "l_f = np.random.randn(N) * sigma + mu\n",
    "\n",
    "# define time domain\n",
    "dt = 0.25  # fixed\n",
    "nsteps = 1000\n",
    "T = dt * nsteps\n",
    "\n",
    "tspan = np.linspace(0, T, nsteps)\n",
    "\n",
    "# define coupling strength\n",
    "K = 0.03\n",
    "\n",
    "\n",
    "# define interaction matrix\n",
    "A = A_2NN(N)\n",
    "\n",
    "\n",
    "# # simulate cells, crucial step!\n",
    "X = odeint(model, l_x0, tspan, args=(l_f, K, A))\n",
    "print(X.shape)\n",
    "# reshape X to match the interaction grid\n",
    "X = np.reshape(X, (X.shape[0], n, n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(4, 5, figsize=(20, 16))\n",
    "axs = axs.flatten()\n",
    "\n",
    "# showing only some of the time steps\n",
    "i_show = np.arange(0, nsteps, int(nsteps / 20))\n",
    "\n",
    "# looping over the time steps to plot the oscillators\n",
    "for ax, i in zip(axs, i_show):\n",
    "    # in order to visualize the oscillation on the grid, we need to take the sin of the phases\n",
    "    datagrid = np.sin(X[i, :])\n",
    "    heatmap(datagrid, marker=\".\", ax=ax)\n",
    "    # here it is crucial the part with abs(order(X[i, :]))\n",
    "    # it takes the avrage of the thetas around the circle and than takes the absolute value\n",
    "    ax.set_title(\"Time {0:.2f}, R = {1:.2f}\".format(tspan[i], abs(order(X[i, :]))))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HL7k7OO2fkWo"
   },
   "source": [
    "**Question 2**\n",
    "\n",
    "Do the same plots as for Question 1 but for an all-to-all interaction matrix: $A_{ij}=1$ for all $i,j$. Set $K=0.001$.\n",
    "\n",
    "How does the space-time synchronization dynamics change with respect to the previous case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a vector of intrinsic frequencies\n",
    "mu = 2 * np.pi / 24\n",
    "sigma = 0.1 * mu\n",
    "l_f = np.random.randn(N) * sigma + mu\n",
    "\n",
    "# define time domain\n",
    "dt = 0.25  # fixed\n",
    "nsteps = 1000\n",
    "T = dt * nsteps\n",
    "\n",
    "tspan = np.linspace(0, T, nsteps)\n",
    "\n",
    "# define coupling strength\n",
    "K = 0.001  #\n",
    "\n",
    "\n",
    "# define interaction matrix\n",
    "A = A_alltoall(N)\n",
    "\n",
    "# simulate cells, these is the vector of phases\n",
    "X = odeint(model, l_x0, tspan, args=(l_f, K, A))\n",
    "\n",
    "# reshape X to match the interaction grid\n",
    "X = np.reshape(X, (X.shape[0], n, n))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(4, 5, figsize=(20, 16))\n",
    "axs = axs.flatten()\n",
    "\n",
    "i_show = np.arange(0, nsteps, int(nsteps / 20))\n",
    "\n",
    "for ax, i in zip(axs, i_show):\n",
    "    # X represents the phases therefore we need to take the sin of it to plot the oscillators\n",
    "    datagrid = np.sin(X[i, :])\n",
    "    heatmap(datagrid, marker=\".\", ax=ax)\n",
    "    ax.set_title(\"Time {0:.2f}, R = {1:.2f}\".format(tspan[i], abs(order(X[i, :]))))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AEjMXuDAfkWp"
   },
   "source": [
    "#### Question 3:  \n",
    "Generate a random adjacency matrix where $\\mathbf{A}_{ij} = 1$ with a probability $p$, where a ‘success’ will generate $\\mathbf{A}=\\mathbf{1}$ and $\\mathbf{A}=\\mathbf{0}$ otherwise.\n",
    "Which probability is needed (approximately, based on the final state) in order to retain synchronized oscillations given your chosen $\\sigma$? \n",
    "\n",
    "Use again the function _heatmap_. You should plot only the final state for the different p (both bigger and smaller than the critical probability). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time()\n",
    "\n",
    "fig, axs = plt.subplots(4, 5, figsize=(20, 16))\n",
    "axs = axs.flatten()\n",
    "\n",
    "P = np.linspace(0, 1, 20)\n",
    "R = np.zeros(shape=(len(P),))\n",
    "M = 5\n",
    "\n",
    "for k, p in enumerate(P):\n",
    "    # use function that creates matrix of random connections\n",
    "    # A = A_alltoall(N, p)\n",
    "\n",
    "    K = 0.001\n",
    "\n",
    "    # simulate neurons\n",
    "    X = odeint(model, l_x0, tspan, args=(l_f, K, A))\n",
    "    X = np.reshape(X, (X.shape[0], int(N**0.5), int(N**0.5)))\n",
    "    datagrid = np.sin(X[-1, :])\n",
    "    heatmap(datagrid, marker=\".\", ax=axs[k])\n",
    "\n",
    "    r = 0\n",
    "    for l in np.arange(M):\n",
    "        r += abs(order(X[-(l + 1), :]))\n",
    "    R[k] = r / M\n",
    "\n",
    "    axs[k].set_title(\"p = {0:.2f}, R = {1:.2f}\".format(p, R[k]))\n",
    "\n",
    "end = time()\n",
    "\n",
    "print(\"It took \", end - start, \" s.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VVB7j1FAfkWp"
   },
   "source": [
    "#### Question 4:  \n",
    "With $R$ defined as $R=\\frac 1 N \\sum_j e^{i\\theta_j}$, plot the function $R(p)$ and find the critical probability (the probability of the neuron connections when the system becomes syncronized $\\sim R\\ge 0.5$). Use at least $10$ different values of $p$ for the plot.\n",
    "\n",
    "**Hint:** to have a \"reliable\" estimate of $R$, you could average it on some time (at the end of the simulation).\n",
    "\n",
    "Which is the critical value of $p$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "executionInfo": {
     "elapsed": 25,
     "status": "ok",
     "timestamp": 1638526869096,
     "user": {
      "displayName": "Lorenzo Talamanca",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GhUF0b2-P5yq7rSXU0BUszliuOAWwm2O82n9owJyh4=s64",
      "userId": "06843078895857810674"
     },
     "user_tz": -60
    },
    "id": "WytDzCuhfkWp"
   },
   "source": [
    "$R\\geq0.5$ is attained for $p\\sim 0.5$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B. The Kuramoto branches (Optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import important libraries\n",
    "from scipy.integrate import odeint, quad\n",
    "# from IPython.display import set_matplotlib_formats\n",
    "from matplotlib.markers import MarkerStyle\n",
    "# set_matplotlib_formats(\"png\", \"pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take the self consistent equation:  \n",
    "  \n",
    "$r=\\langle\\cos (\\theta)\\rangle_s=r K \\int_{-\\pi / 2}^{\\pi / 2} \\cos ^2(\\theta) g(K r \\sin (\\theta)) d \\theta$  \n",
    "  \n",
    "And use it to plot the 2 branches stemming from the Kuramoto's model in the plane (K,r)\n",
    "\n",
    "*hint*:\n",
    "\n",
    "the curve is the made by the set of points in plane (K,r) such that:\n",
    "  \n",
    "  $0 =r K \\int_{-\\pi / 2}^{\\pi / 2} \\cos ^2(\\theta) g(K r \\sin (\\theta)) d \\theta - r = K \\int_{-\\pi / 2}^{\\pi / 2} \\cos ^2(\\theta) g(K r \\sin (\\theta)) d \\theta - 1 $  \n",
    "\n",
    "Define function that finds $ F(K,r) = K \\int_{-\\pi / 2}^{\\pi / 2} \\cos ^2(\\theta) g(K r \\sin (\\theta)) d \\theta - 1 $  \n",
    "then find F(K, r) for a mash grid of points, and finally select the curve such that F(K,r)=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=1\n",
    "Kc = 2 / np.pi * np.sqrt(2 * np.pi) * 1\n",
    "\n",
    "\n",
    "# write the function here, use quad to integrate\n",
    "def F(K, r):\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# once define the function F(K,r) we can plot it\n",
    "# we use the contour function have the F=0 slice\n",
    "\n",
    "contr = np.arange(0, 1.1, 0.01)\n",
    "contK = np.arange(0, 5, 0.1)\n",
    "KK, rr = np.meshgrid(contK, contr)\n",
    "V = np.array([0.0])\n",
    "Z = np.zeros_like(KK)\n",
    "\n",
    "for i in range(KK.shape[0]):\n",
    "    for j in range(KK.shape[1]):\n",
    "        Z[i, j] = F(KK[i, j], rr[i, j])\n",
    "plt.contour(KK, rr, Z, V, colors=\"purple\")\n",
    "plt.plot([0, 5], [0, 0], color=\"purple\")\n",
    "plt.scatter(Kc, 0)\n",
    "plt.xlabel(\"K\")\n",
    "plt.ylabel(\"r\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "name": "Solution_11-graded.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "bio_341",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "title": "Exercise 11: Collective synchronization in populations of coupled circadian oscillators"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}