{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "Stochastic timeseries can be characterized by the noise color, which is defined by the slope of the power spectral density. The more negative this slope is, the more structure there is in the timeseries. This also results in a more slowly decreasing auto-correlation function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-19T15:19:47.209082Z",
     "start_time": "2020-02-19T15:19:46.783523Z"
    }
   },
   "outputs": [],
   "source": [
    "# Data manipulation\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# Options for pandas\n",
    "pd.options.display.max_columns = 50\n",
    "pd.options.display.max_rows = 30\n",
    "\n",
    "from IPython import get_ipython\n",
    "ipython = get_ipython()\n",
    "\n",
    "# autoreload extension\n",
    "if 'autoreload' not in ipython.extension_manager.loaded:\n",
    "    %load_ext autoreload\n",
    "\n",
    "%autoreload 2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "%matplotlib inline\n",
    "\n",
    "import time\n",
    "np.random.seed(int(time.time()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specific imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-19T15:26:42.782519Z",
     "start_time": "2020-02-19T15:26:42.747678Z"
    }
   },
   "outputs": [],
   "source": [
    "from generate_timeseries import Timeseries\n",
    "from noise_properties_plotting import PlotTimeseriesComparison, example_noise_fit, noise_cmap_ww, noise_lim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Settings figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-19T15:20:00.980685Z",
     "start_time": "2020-02-19T15:20:00.648003Z"
    }
   },
   "outputs": [],
   "source": [
    "from elife_settings import set_elife_settings, ELIFE\n",
    "\n",
    "set_elife_settings()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-19T15:25:11.739430Z",
     "start_time": "2020-02-19T15:25:10.950249Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAADtCAYAAADOSW5rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxcd3nv8c9zzuyjfbNky7udfSdAHCCx2RsSStmhpaWldAFabimFUm7L7aUbpSsNSwvthbK0lNwWCBdoSsFZIITYSZzNiRPHji1btrVLM5r1nOf+cY40sixpRrYWW37er5dfHs3ZfvOd35znLDPniKpijDHGGHCWuwHGGGPM2cKKojHGGBOyomiMMcaErCgaY4wxISuKxhhjTMiKojHGGBOKLHcDZtPW1qYbNmxYkHlls1nS6fSCzGsls5xqZ1nVxnKqjeVUu4XKavfu3f2q2j79+bO2KG7YsIFdu3YtyLx27tzJ9u3bF2ReK5nlVDvLqjaWU20sp9otVFYi8uxMz6/4w6e9o3kypZMvUKCqFMv+MrVo+ZU8n3P1og2qWlPbcyWPsnf+vsdQe1bzmZ/nn5v9ZoJ/jrf/dC1GX1iOLJdivXXW7ikulL3HM7z3R+N85NE7EQSAoucTc4PtAUWZyNl1BF91cjwAXxVfwXVAEBTF86Eu7pIreThSmcZXJeY6FDwfAWKuQ9HzcUSmNwvXEVQ1HFeIOCePo4TLnTKtSNCeYlmJuDI5TAmWX/L9yflMtBVgaj8qeUoi6qBaee0x16Hk+4yNZWl48K6T2jCRxcRjRYk4TrCC1JPbWPaV9roY/dkiAI4E7RCptGFimRNNkvB1TYw30V5fFUcqr8HXyriqQX5T36PJjML2xiIOZc/HC4dNvD8T03vh/CfaONGuibZMtHVi2Y5MtCsYJ5vJ0vTQXZR8ReCU93jq/KYuP3jfg/mVfZ2xb8zEU6UlGWU4XzrlvZ2a09QcZspw4r2ayGLi+XjYV6e+/qmvZfqyGuJRMsXyZD+CyutxpsxjLJOl7sE7J9s88T4G8w0eT7TclZPfh4n2TjznSNBHyr5OZjhh6uuc6F/OtNcytT9PlYq5ZAplomEGUz9zvgZznfpZizrBZ3ziOU+VSPi+Tje17YpOri8ijhNsnAIRRxgdy1L/4F2TbZzaLSbWL1PfY1dksg9PX2fN1IayX2njxOuZ2rapy5po60QfmnifJnjh8qKuTH725lr+xDRlL+j/EUdOeS9EKp/76R+Jif6gCu/fvpnVcy7pzK34ovjirW187oYUO3bcuKDz7csUaEnFUFUibmWHO1/ySERdAAplj3jEnXH68WIZ1xHiERdVRWpcOU5sKU0f3/eDzlTrfKZOV/aVWMQJD0vcUHWaTKFMPOIQdU8+0KCqHBnJs6YxMWMbq73OYGsWnHCDYeLxYpgtx1p9/wc/4IUveiGxyPwOtni+4kjwv+vIvN73gfEiLclY1UzKno+InLThUG3e40WPdLy21YGqMpovUxePTC5jtvc26FM3Tu5VzNT2kuef0pdmM9HPofLeTe8rvq+TfahavqpKtuiRjrmThaPaNLmSRzJa+VzPp69OvO8lzz+pkN511501ffYmluf5wXqn7PknrX9mMvUzXota5lltebOti+azrpvNzp1Pn9H01az4oginv+KbS3tdfGLuJz2fmPJhma0gAqRilejn077Zxj3d4uE4Qmye09bNsvIUEbqbkrNOV+11ishJK7xFeNtqbks1jsi8CyJU9nAj7vyWLyK0pePVR4R5r9BEpOaCODF+YzJ6ynNzmat/1loQZ5vP9L4yMU4t77GITPbnaI3vydSCONPy5zLx/k+85tPphiIy2X9qea/n+xk/k4I4sbzZLMa6eKGt+HOKxhhjTK2sKBpjjDEhK4rGGGNMyIqiMcYYE7KiaIwxxoSsKBpjjDEhK4rGGGNMyIqiMcYYE7KiaIwxxoSsKBpjjDGh8+IybwDqlfBLWbQ8jjd+HC1lcZJt+PlBSkP7iDSsx022gxNBnCg4EbSURSIpIvXdIO6p1xstjuHlB3FiDTixevziGE6sgfLIfsqZo7jpVUQaNgCCuHH8XD9erg+JJNBSFjfViZcfRNxoMI9EC17mCF5+EC97jGjjRqLNFwDg5QcBcCJpyuPHQD38/DASSeBle4k0bkLcOOJEkFgDWsrilzLhlbh9xInh1q3Bz/XhxJsRN7hM18Q1QL3MESL5HrJPfwNxIjjJdkTCbSaR4DVEEkSbtuAXRpBoGvXy+IVhInVrUb+IFjP45SzixMCJ4GWO4pcyOLEGYu1XhFf2dfAyPTjxJiSSwsv14WWOEm25EMRBC6NIrI7SwN5g2Y7LxKX0vEwPbqqTSMM6QHCSbZNt94ujRBs34+X6KA09iRNvxk214yTbAcUb68Evj+NEUkQaN1Ia2kd59BDF/odJb30diEukYR3e+PHgNdWvp3DsPrQ0TqRhHZH69eDGKA3uDfLzxoM7BeQHkGgdfmEIN96CapnCsfvBLxPruBov14ebbEciSfziKCJO0A+LY0E2iRZiLRdNeS/Ci3cXhhE3gRNNTQ4rDT2JlnNEm7ZS7H8YcIi2XISbaMIvjlEeOYCTaMZNdVbeX7+MXxjBiTfh5/qReEPw/qiPuFG83ABe7kQwTSSJP34Ct24NiINfGMYbP4ETq8NNdeIXR9FSlkh9N172OOXs0aB/ikNp8An8wjCJ1S9E3Ch+MUOx/2GiuYOo76FeHi1mgmlaLkJL4+FnyqGcPYq4CbScpTx2hFjbpeF7sI7C8QeItV+JN34MJ5pGvSJOohl8DyfZil8YpjS4D4kkiNStRiJJSiPP4KY6gtc00/U3vRJaHseJN+Jlj1Po24MTawAU9Us4sXq88T7cZCtOvJGJy5ZrOYdE64Is8wN44ydw012oVwS/hFu3GnETlfWHV8DLHsNNd+IXhnAiKSSaxi8MB+9xKYs4UcpjPeCXUFXKw0+DOLjJtsnXEm3chLgx/HIOvzASfCTdOKgPKE68icLx3cRaLkKiaRAXP9ePE29E3Bhefhgn3ohfGMLLHsOJ1gXrHXFw61bjFzO4iWZwY5MZiROlNPxUJbNyAb8wRHzVc3FidfilcfBLQRuHnyZSv3byyt2lwSdxUx1Bdk4keJ8n/hcJ18Vj4XzzQb+MNeCmO9HiGBJJUux/FPWLxNuvBHEm18fqFeZe0S8AOVtvIXTttdfqQtxPsXBsF09/+300xXM4iRYkksRNtIATRUsZnHgT0ZaLKY/sxy+MoH4J/DLql5FIEi3nKGeOgHrBDFVRLQdXso+mcVMd+MUMfnEUJ5IMPszNW4nUr8PLHqM8ejCYxivgxBtx69agXh5x45RHDuKmO0F9/PwQ3ngvkaatuIkW3HQXpcEnKI08gzhRJFaH4OAXhog0bgQnghNrwC8M46a78DJHJ9s+uUJNNAdtdly0nMfL9uIkWoIPpV8OOpv6gOCmO+kfHmfNhS8CFG+8H9CwqIaFMz9IeeQZ3GQHfnEUUNy61XhjRyASx4mkkVgaLedBfdx0V/gBHKDU/0iYeZZI40a88eOgipvuxE22U+zbA24MN2xftO0yBAfFB98LlpXuxM8NUBraBwh+YSi4a0ddN06sYXJlGG2+ED8/gJ8fojx2GJwIkbo1wcZIOU957BBuuoto42airZcw/vTXQSTIJ9mOm2ylPHKQaPMFROq7KY8dojx6CL+UJdZ2OX5+gMHjz9CQdHGSbcEGVrw52NgqZ0mufwUgFI79hEjDerz8QDBOrBFQJJrCidbjxOop9j+KX86CVwrfq2A71YnV4+eHwo0CQP1gw0ccytleYm1XgHphMRoBkaBthZGgHV4xvOWHO9lPnGgdfimDlseRSHCN2mCjb224oTiOk2zFCzNz4k24yXa8XF/Qp6JpxE0EGwJugkjDesqZnqBtdd24yTbyR+4ONpjKeeJd13H84B6aEoVggy2SJFK3htLQvuCz5ZfQUiYokuVc8HlKtFDsfxQ33UWx/2FibZfjZY4GxaecBTeOnx8GvxR85uLNRJsvCAtQL35xjGjTFrzxE3i5E0xem1jD24iph6K48Wa8/CBuahXxjqvwS1kQB5EIfnEEN9WBl+tHS1lQP7irQ9h/vfwgWhwj0rQZ/HJYiJzgM+gVw89hCXETuKkOypke3GR7UADKWZx4U/hWJ1G/jJtaxeAzd5FOp4g2bQHHxc8N4JcyQcajzwb5RFK48eZg49EvhbeO8PFz/cTar6I48FjwmfbLOMlWtDiGlvNIrC7of4nW8HX14cTqcaJ1lEYP4kRSaDmHqjeZl3oFIg3rEHFBFYkkw+J7f9iWJOJE8QsjRFsupJztRcRF/TKRdCd+aRwtj6N+OVh3+h6q5cn1jZNoCftfAjfRHOxgZI8hsXrwirj13biJFgp9e8LPf6D5+j/kvv0s1P0Ud6vqtdOfX7Y9RRG5AXiXqr55MZfjxBsY7volLv2pty/mYlaEJ3bu5LJt25e7GTVTvwwI4sx+4fVaNFz+jnlPs3fnTi6e4YOpfjnYKq6RqgYrVjdafeQFoupXjgIsAr84FhRQcXh0504utZvnVvXYzp3ceOONs14weyHuLrFi7N+5qLNftKIoIpuADwONqvp6EUkDnwKKwE5V/bKIXL9Yy58Qbb6AUvLoYi/GLIP5FJ+lMt82iQgsYUEMlrm4XyVwYvWLOv+Vaq6iZwVx6Sz64VMRuS0sim8DhlX1dhH5qqq+SUR+V1X/bKbpOjs7tbGxcfLvm2++mVtuueW02pDJZKirqzutac8nllPtLKvaWE61sZxqt1BZ7dixY9kPn3YDj4SPPRG5EniRiDykqt89ZeTubhbinCJM3Oh0+4LMayWznGpnWdXGcqqN5VS7xc6qalEUkYiqlhdgWT0EhfEhwFHVPcCrFmC+xhhjzIKY8+SCiHwU+Hz4+K/nM2MRaRWRzwBXi8iHgH8HXicinwZuP73mGmOMMYun2p5iPbAvfFyaz4xVdQD4tWlP/+J85mGMMcYspWpfQ1OgS0RuBjqXoD3GGGPMsqlWFD8C7AHWAu9e/OYYY4wxy6fa4dPXqOpnAETkLcC/LH6TjDHGmOVRbU/x8imPr1jMhhhjjDHLrdqeYoOIvIPg3GLLErTHGGOMWTbV9hTfBRwFesPHxiwbVcWfuLCzAcDzfX7Sd4jBwvjk30VvIX5WfCpfffrzWbKlAo8O9fKpvT/ksaFji7Ksc4nn+xS8MsOFHP35LLv7exgr5VFVerLDeP7S9Nmz9eYO55pqe4o7gJ8G4sAbgF9a9BYtsK/sf4CPHrmbm34ySsqNUvDLDBfzpNwoP7flOTgIPzxxgNsPPU7P+DCrkw2ko3HWpptYm26iO93E5c2ddCbrccShOZbk6bF+DowNciAzyIGxAZ4eHSBTKrChvoXLmjppiCV4ZKiXl3RtJRWJ0pGsoyc7wsODR4k5ERpiCe46tp8DY4NkygVU4bLmThpjSVriSdKRGD/pP8wTw8dJR2LUReP05bNc0rSKQ5khkpEoG+paeO2Gy9nWvgERaIwlyZdLPJMZoDmW4v7+w+zuP8y2jg1sqm+lNZ7CEeH/PvsIDsLx3BiHs8OUfI8dXVs4nhvj2yce4Ed7yhzKDOOpT8/4CEOFcS5pWkVnsoGRUo6nRwdojCW4uLGDzmQ927u2sLm+lZgb4ZmxASLi8GxmiDuP7efo+CgXNLbzwlUbiTsRnhg5QX8hy4+OH+SCxjYOZ4P59+Uz1EcTXNTYQW9ulEeHjpGOxMiWi0SdynZbeGcaIuLwyu6LeF77On6q++JwmOKpz57BXlYl69g7fIKOZB3j5SLfO7qP7nQTOzq3sK6uid7xMR4bPsaTIyfoz2f5ds9efFXSkRhtiTryXok3bLySA2ODtMRTXNy0iqQb4Yrm1TTGEif1r/FykS88vYuHB4+yd/g4FzR2sDbdyH19h+hKNnBVy2qubu3muo51HMuN8dVnHqK/kKXs+zy3fS3dqSauaOkiFYmRL5foGR9hV/9h7jl+gE31LUHuxRyvXncZq5J1ZMtFEm6U2w8/xpeefoD9Y/1c09rNo0O9pCIxRoo5BOHS5k66043cvPYSViXrua/vEN87uo9HBnu5rLmL61dtwBWH57Wtpaw+T470sXf4OC9ZvRVVZf/YAAk3Sn00zn88+wj7R/sZKxVYlayn4JdJR2Lcsu5S/uCB73JkfIQbOzezsb6FjfWt5Molns0MUvZ9WoqjbPPK3N9/mEeGesmXS/TmRtnZu58tDW0kI1EOjg0yEPazumicN2y4khs7NxFzI5R9j4OZIR4Y6OF4boyo43J//2Ei4pApF3lwoId0JIavSsEr0xBLEHVc2uJp9o4cZ0tDG6uTjdRH4zwzNkBTPElEHJpiSU7kM7xizYVsqGuhO93IkfERdvUd5r+OPkXCjfDyNRcyUszRnW7CFYeh4jgdiTq+eegxhos5Lmzs4FB2iCdHTuD5SjoaIyIOG+tb2Tt8nEypQGeqgWPjo1zRshpF2T/az8b6VrpSDUTEIeI4PDF8ggMnetl27wCN0STHcqNsrG/h0qYuFOV4bozVqUYeHeol7kb49uG9bGlooymWZLAY3MLsseFjjBTzvGDVRo5kRxgoZFGF+mict2y+hk31LSTdKFe2rGawME7McelMNeCrz8ODvVzQ2E4qEsNXn8eGjtORrCPuRDicHaYlnmJ1qoH+Qpb+fJae7AiuCDu6tsx4TdaxUp7d/T3E3QirUw2UfJ9UJMqBseD2d187sIeRUp4nR07gq1L0PZpjSW5ZdykPDPRwYGyQbLkIBPfsiDku77nkRaxNN+IvQeGf89qnInIr8LuEv1FU1cW/mVVooW4dBXDHD/4bf+saFCXuRKiPxjmUHea7PU9QVp/r2tdzU/fFrEk3MFzMky0VOZwd5nB2mJ7sMPecOMBgYRxflf58litauthc38aG+hY21bewqb6VlBvlcHaYhwaPMlDIcnVrN//Z8wQiwpHsCGvTTVzVuoaS7zFYyPKiVZu4qKmDdCS4h9mu/p5ga7OYY7RU4JrWNVzU2EG2XKToeyTdKLsHeri0aRWjpTyHsyN8/dlH2NV/mPFyibgbIVMqcHHTKkaKeS5r7uSq1jXsGTjC/rEBTuQzwcp+w5Uk3CirkvV0pxtRhbuPP0NrIo1/4BhNF6xnbbqJuBOhO91IXTTO3uETHBkfoSWe4sLGdkaKeR4Z6mXv8HEeHuzliZETxByXVcl6Io5DYzTBLesuZU2qkQcHenho8Cg5r8RFjR00RBNs69jAgcwA3akmOpJ1NEYTZMoFnhrtZ1WinkuaVlHwy0TEIeaeut02WBjn+71P8W8H9jBUGGe0FLxnrjhc0tzJWCnPxY2r6MtniDguL19zAUfHR7njyJMcz43RmWxgQ30zz2tbR3M8xfUdG2hLpMmWCvTls2TLRb51+HEuauxguJjj8eHjjJeLPDzUy0gxj6c+xfEcDfX1+Kq8aeNV3NC5ma0Nbewb7ePo+CjPae2mL5/hocGj3Nd3iN39h+lI1vOadZeytq6ZiDjs7u/hcHaYR4d6yXtlYq5Ld7qJCxra2dG1hd7cKAfHBmmMJbj98OP05TIk3Ahl9XnZ6gt525bncHHTKiDcQ/Q94q6LKjw5eoLD2RG+9PRuin6Z57at5cVdW7m8pYu9w8f5ryPBz48fHuol4UZoT9SxJtXAkyN9eOpzUWMHJfUZKea4rmMDN6zaRNRxiMxwR5Khwjj39x/m6dF+9o30kY7GWJduRkT46p4fccQtcmPnZi5v7qIxlqQrVc9zWrs5MDaIiLAu3URLPMWjQ8cYLI7zpad3s3ugh20d63l06BhNsSRXtaxmdaqR8XKR7V1bKPplXHF4Tms3nvq44iAijJeLeOrz1Eg/V7eu4YmRE4yVCoyVCqxONTBeLqEEn+PGWILvHd3HsdwYhzLDtCfSPL99PS9ZvZVcucQdR56kKZ7kRC5DWT3qInEGCuO8Ys2FdKUa2Dt8nO50E13JejpTDbOuf0q+x97h4zgibKlv40BmkL58hpLvU/I9NtW38tTuhyht7sJXZVN9C8+MDfLoUC8Rx6EtUcdTI308p62bglfm+o4NDBVzZEoFmmLBLcAuaVqF6wR9al1dE6tTwXWjj2RH+NrBPfRkh8l5JR4Z7CUZiVH2PYaKORyETQ2tHMoMBf3a87i0uZP+fBZPfdamm3g2M8SJfIZVyTpa42k21bdO9m1VxQ0vMO+jOAipSJSrW7vJeyX681lcxyFbKrKxPjgDt6NrC6uS9VzR3EXcjRB3I9zfd4if9B3ixq7NbKxrIR2N46uPIw7Hxkf5xOP3sGfwKJ++/nU8s2vPot46qlpR/BDwOSALoKrjZ9ySGi1kUVzIa+UVvfKMK+rlNlrMUxeN4ZzBHRBONyfP98l5Jeqi8dNe9unqHR+lPhqf3LhY7LsJjJeD4nvv3ffYtSprcLp9aiCf5emxfppiSS5s7Fj4hp1lzvZrn55Nt65aqKxO936KFwAfY+LW0+fg4dOFdjYWRICGaYf1lpLrONQ5S18QAbrm2EJfDKlIrPpI5oy1JtK0JtLL3QwTOlsK4lKYcw2vqr8oIu1AionbrxtjjDEr1JxFUUQ+DlwHPAVsBV60FI0yxhhjlkO1E1AOsFNVfwn4+hK0xxhjjFk21U6QPQO4IvKPQHIJ2mOMMcYsm2rnFD8JICLNwPCStMgYY4xZJrMWRRH5GrCG4DeKRaAReN4StcsYY4xZcrOeU1TVNwDfV9UbVfVlwFeWrlnGGGPM0qt2TnGLiNxIsKd48RK0xxhjjFk21YribwBvCh9/eJHbYowxxiyraj/JWAXUAa3YXTKMMcascNX2FN8H/BXhBcGNMcaYlaxaUXxUVR9dkpYYY4wxy6zq/RRFZDtQAFRV37j4TTLGGGOWR7Uf79+yVA0xxhhjllu1C4L/C8HdMeqAdap61ZK0yhhjjFkG1fYU3zLxWET+x+I3xxhjjFk+1fYUJ36GEQVOuUPxmRKRG4B3qeqbF3rexhhjzHxV+6LNYwSHTwvAP8w1oohsIviBf6Oqvl5E0sCnCK6Gs1NVvzx9GlW9S0SuP62WG2OMMQusWlGMAL8C+MDngP+ebURVfQZ4h4jcFj71WuA2Vb1dRL4qIlngNeGwO1TVrqVqjDHmrCKqOvtAkf8DvAMQ4LPhzYbnnqHIbeGe4oeA76jqQyLyFVV96wzjXgn8CfB3qvrdqcM6Ozu1sbFx8u+bb76ZW245vS/DZjIZ6urqTmva84nlVDvLqjaWU20sp9otVFY7duzYraqnnBac69ZRKaCH4PZRChyd5zJ7gG7gIWa5nJyq7gFeNdOw7u5udu3aNc9Fzmznzp1s3759Qea1kllOtbOsamM51cZyqt1iZzXX4dNPEhTDPwRigDvXjESkFfhj4OpwL/ETwK0i8irg9oVprjHGGLN45iqKvw68HrgFuAB421wzUtUB4NemPf2LZ9Q6Y4wxZgnNdZeMJ8L/3wr82K6BaowxZqWbqyj+DMFvE/8PsFlE4kvTJGOMMWZ5zHr4VFUfBB4Mi+HrgS+H/xtjjDErUrWbDKOqBVX9sqpaQTTGGLOiVS2KxhhjzPnCiqIxxhgTsqJojDHGhKwoGmOMMSErisYYY0zIiqIxxhgTsqJojDHGhKrdT9HUKPfsgxQO7cGJp0ld8mIida1Vp/GLeSQaR0SWoIXGGGOqWfFFsTTYQ+zITygeX0u0YxPF3idwEvVEW7pPa35+IYs3PsL4vrsZ/N6taCELCE6qkfqrbqE0cIiB7/4l8bVX0vKSdyNuBNwo+QO7UK/I8A//GT87BKpINI6WS4Dil/I0PPf1pC99GZG6VtyGDoon9lMe7iV94Q0Ujz9NpHn1vNvt5zOUho8Sa9+EeiVGfvQlCkceJ33Zy6m77GWgCo6DOHPeBKVm4/vvY/TeL+Mk6omvvRI31URyy3W4yYYZx1fVyY0CLzeKOC5OPL0gbTHGmPla8UWxeGI/0aMPcvxf76XUf5Bo+0a83AjiRklf8hIijV14Y304yQbcVCPDP/wi5cEekpufT+7ALty6Vpx4CnyP0sAhcKO46RYSay9n7bu/RqSp85Rltv307zN2/230f/tjCIJfyJDY8BwkmmT1L36WaPvGoCg6laPXWi4ycu9XGLv/NspjfZRHjxPr2IKbamTg2x8nvvpi8ocfpuF5b6TlZb+BEz31UrSqSubh7yIiaLnA4Pc/g5cZINq2nmJvcH33huveQt2VNzH24Dfp+/ffB8dFS3lSW7YRp4vxfRFiXReB+ogbw0034RdzIA7F409T6nuGfM8jlAcOUzj2JJGmLmKrtlIe7KHY9wxuuoWmF709KMZ9z5DPDHD8q79D/VW3kD/0EOXhXnAcom0b8EaOB++FOCACbhQt5UluvJb0pS+j7sqbZi2m5xotF/HGR9ByEdRDfY/y0BFK/c8isSSoT7HvAJH6drzcKPHVFxNffQnRljXzW47vn9Sv5jXtlA2UxXQmbTTLQ1XJH9yNuDES664443kBs/Y1LzdG/sD9RFrWknv6XrzsIKmLbqR4dC9OsgFY3Jsxy0QDzzbXXnutLuZNhgu9T5I/+ADlkV6ceB1+fgy/mKPu8lcQ67qQ8X33kL54B15mEC0XwHGJtnQv2B7V6fCLeQa+8xeM7rqNlpe8m4bnv5nRH/8L5eFeJBKjNHCI8kgvsVVbQRyabvgl4l0XAsGKCL+MRGIzzDdH7pmf8Nh3Ps/a5gTFvmfA9yhn+hFxkGgC9UpE2zYQ69hMtG0Die7Lia++mNJgD8Xj+4i2bSDS0EG0dd0p8/dyY2Qe/g6JdVcS77oQ9T2KvU/i1rcRaegI2lAq4ETjeNlhisefIrPn22Qe+S5aLiLRBNG29fi5ESQSp/7qV9PwvDcGGyuhqR+07BN3cuJrvwcixDo242WH0FKO8lgf4kapv/rVJNZeSfbJuyideBq/MI6X6UciMSSeJtF9BbHOrfiFLJmHvoWTakLcKLFVWxHHJfPoHYyNDtO0ejOtL38vqYt3UDjyGJH6diQa5BdtXkPx2D7G9vw/Mnu+jUTjROo7kGgcfA9EiLauJ9q2Pj5oIXMAABx6SURBVCiUCLH2jZSGjuAk6ime2E/+4G783AjJrS8g0tSFiIOWCjRe/3NE2zdOFha/mCO3/z6Of+W38MsFOl73x6S2Xk+kcdUp70N5pBcvO8TgHX9LoefR4CiBuKh6UC5BJAqqxNo3IZEYbl0r9de+jvyB+4mvvoT0pS89KffpVJWh7386PALh8tS+J9nU1Up89UU4qWZGfvRFxp+4k/qrb6H1pg9MHvkoj/ZRPPYk0baN894QWAmW6ibD+WcfIvvETuLdlxNtXUup7yDj++4m2rKW2KotSCRG9ok70WIO9cskN15L9ok7GX/yLhJrr6A8FnxOoi1rcRJ1pC95CfGui/ByIxQOP4xfzNFw7Wtx4nXkDu5i/Mm78POZ4HM4Pkz+4APkDz2IOBHU93DTzcG6t5Alsf5qUKV4/ClSF95I8cR+Emsvx21YReHwHpxEA+2v+QPuvn/PgmQlIrtV9dpTnj9fi+K5zMsOM/CdvyDz6H9Sf81rSHRfjnpFAOqvfvWMha8WZ2tOXnaI8sgx3FQzXn6M0R//C6O7/i/RlrV4mf6wqACqqO8RW7WFVW/9a9x0C6UT+3HSzTixFG59G1rIMvbgNykceZzExmtJdF8eFN2WNfjFHH4hS/7ZByke2xcU0GtfixbH0XKRwtG9+MVx6q9+NXf96D6u29TC8N3/xPi+e4h1XkB59DhazBNffTHl0RO4yQYatr2V1NYX4KYaT+u1T2yweGN9+MU8qMfY7v+gNNgDfrDHKW6USFMXXW//e1CPwf/6O3IHduGNngiKMAACjkOsfRM4Li0v+w2Sm687aWt9cl2gPqWBQ6hXpnhsH9lH/5PEuqsp9D7B+JN3oqUCkZa1uHWtiBslf+hB1CsR69hMqf9ZUluvJ9q+CYkm2PfUU1xyzTYKRx7HGx8mfcmLqbvs5WQe/g593/wjnHgdkfp2iv0HSG6+jsLRx9FijtiqrTixJOnLXk79VTdP9unZ9jJLQ0cpDx0hsfHak15TefgY5bE+4t2XnbJn4hdzwUau4866seuXCpRO7Kc81o9b10p8zSWU+g6ACF5mEHFcou0bcWLJySMd3lg/TjyNk2zAyw6SP/ww5cGe4OhIdgg/P0r20f8i3n0ZbrqFuitv4kf3P8ANL70JcYNiMbU96nt4mQEKRx4n3n0Zkfq2qv1Gy0Wye39A7uBuxh74Bm6yEVUPLYxTd/kr8bKDlIaOINE4jc9/M6XBHgqH9wBC6sIbkGgC1CN3YBfpS19K+sIbkUgUgHzPo3iZAbzsEOP77qF4bB9Osp7E2itBhLHd/4GWCiQ2XEP6ou3gRin1PUOkYRXxtVeQWH814jj4pQJ+bhQn2YBEYhQOPwzizPheTbVQ6ykrimfhyv5scy7lFBx+PIrb0DHjoeTFdi5ktViHKScO/XrZIfxiLljJhSu1WMdmnETl8Fa1nMqZAcqDR4h3XzpZCLzxEUr9z+IXMow98A0yD38HicRQvxzs1frlySMYTTe8A3FcBv7rb4l3Xkix/yCxjs2TpwucVBNuupny0FGad/wq6UtejJcZZPieL5B99A4knkbLBcRxkXgaN9mIE27AlAd78LJDYSFqp3jiaUoDh4l1bAbHxUnUA0qp7wB+cRy/kMVJ1BNp7MTPjeKPD+M2dASHwdvWU+o/iJtuwUnUBRsARx6jPHKM3IFdDB4/Qn2kHITi+5N7UDgR8MtBQV59MeNP/5hYxybKI8eItqwj2roWv5Ald3D35EYSXjCf1AUvJLnpedRdeRNaKoDjTB6ZOZctdlFc8ecUzcokjku0de1yN+Ostljn7YLs151yqDyx7sp5zytS13rKN7XdVCNueN4qtfV6Ot7wp+Fyg9ejquCVKI+eYPQn/4ZfyrP+d+4g0tCBlx2iNHQ0KM6xxOQ8y6MnGLzjE4w98A0kmqDhua9n1Vv+YrIQqypaHMcbH8HPjQJKpGn1vPbw53tONrX1+snHz+zcyVVTVvTlcG9TSwW0XJg8FK6+z/i+u4PDkZkBSgOHcJKNtP307+PEUqhXWjHn4ZfLeVEUb7/99rN+q/5sYDnVzrKqzULkNL24iwhEYkRbuml95ftOGuamm3HTzafMI9LQQcfr/2j2ZYgg8XTwzefm1afXzjP4ktL0nCYPkcaSJy/DcUhfdCMA0ZbuGTZEkqx0i/3ZOy++Avatb31ruZtwTrCcamdZ1cZyqo3lVLvFzuq8KIrGGGNMLc7aL9qISB/w7ALNrg3oX6B5rWSWU+0sq9pYTrWxnGq3UFmtV9X26U+etUXRGGOMWWp2+NQYY4wJWVE0xhhjQiv6JxkikgY+BRSBnar65WVu0llHRDYBHwYaVfX1IvJWYAcQB35dVbPL2sCzhIi8BngV0AF8Ergc2AhEgV9TOw8BgIhcDLyX4LzPfwMjWH+aUbh+ugv4CHAh1p9mJCLbgY8CjwH/CjyHRcxqpe8pvha4TVXfCbx6uRtzNlLVZ1T1HVOe+pkwr38jyM8Aqvr1MJe3A28BrlHV9wCPAC9czradTVR1r6r+GvBG4FqsP83lgwS5OFh/mosCGSABHGWRs1rpRbEbOBw+9pazIeeQia2uZwnyMyf7n8DngL7wb8tpGhF5NXAPwZ6i9acZiMhLgceB40Aj1p/mcreq/hTBRsSnWeSsVnpR7KES2kp/rQttHUF+BpDAx4DvAPcTHB4Ey+kUqvpNVb0e+NkpT1tOJ9sBXAe8Nfw3cVFSy2kaVfXDh0MEh+MX9bO3on+SER6zvxXIA/fYOcVTiUgr8MfAywj2gJ4FXkRwvah32zmggIj8JvALBAXxISAFrKdyrmzlfpDmITz/81qCXB4mWJFZf5qFiLyd4Dd3F2D9aUYi8lrgFUATwZ7iNSxiViu6KBpjjDHzYYcUjTHGmJAVRWOMMSZkRdEYY4wJWVE0xhhjQiv6ijbGGGNWFhFJAn8NrAaagb3Ak6r6lwsyf/v2qTHGmHNN+POfy4BvAe8h+Pnd54GfEPwEqBd4PvC7wDjw24AA+1X1b2abr+0pGmOMWSn2qeoHROQ/CK7pvBN4JdAJ5MJ/l881AyuKxhhjVorR8P+Cqo6KSJHgR/4O8EVVfbjaDKwoGmOMWeluBf5ERHqBMVX9w9lGtHOKxhhjTMh+kmGMMcaErCgaY4wxISuKxhhjTMiKojHGGBOyomiMMcaErCgaY4wxISuKxhhjTMiKojHGGBOyomiMMcaErCgaY4wxISuKxhhjTMiKojHGGBOyomiMMcaErCgaY4wxISuKxhhjTMiKojHGGBOyomiMMcaErCgaY4wxoWUtiiLyGhH5rIh8Q0RevpxtMcYYY0RVl7sNiEgz8Beq+o7lbosxxpjzV2S5GxD6n8Anpz6RSqXU9/3Jv5uammhsbDxpIt/3cRw7AgyWxQTLocKyCFgOFZZFxb59+/pVtX3688taFEVEgD8DvqOqD0wddskll7Br1645p9+5cyfbt29fvAaeQyyLgOVQYVkELIcKy6JCRJ6d6fnl3lP8DeClQKOIbFHVzyxze4wxxpzHlrUoquongE8sZxuMMcaYCXZw2RhjjAlZUTTGGGNCVhSNMcaYkBVFY4wxJmRF0RhjjAlZUTTGGGNCVhSNMcaYkBVFY4wxJnRGRVFE/llEfmqhGmOMMcYspzPdU/xloF1E/lVE3isi6YVolDHGGLMczrQotgKbgFHgGPCP85lYRDaJyD+KyG1n2A5jjDHmjJ3ptU/fD3xSVZ8BEJHD85k4nO4dVhSNMcacDc7oJsMi8i5V/VT4+H+o6t+c5nxuU9XXT32us7NTp94/8eabb+aWW245abpMJkNdXd3pLHLFsSwClkOFZRGwHCosi4odO3bsVtVrpz9/2nuKIvJxYJuIbAAE2AKcVlGcSXd3t91PcR4si4DlUGFZBCyHCsuiujM5fHor8BBwN6AE5xTnRURagT8GrhaRD6nqn55Be4wxxpgzciZF8WeANcCVBHuKCnxgPjNQ1QHg186gDcYYY8yCqVoURUR05hOPXwViC98kY4wxZnnU8pOMD870pKr2AtsJDpv+GfCLC9csY4wxZunVcvj0hSLyIWAEYOLbpqHLgBcD/xz+b4wxxpyzaimKH5/yePph1Abg7cAvADcvUJuMMcaYZVHL4dNu4GZVvRO4fOoAVf1VVX2TquaB31yMBhpjjDFLpZY9xW3AifDxhqkDROSvgFUExVWBty5k44wxxpilVEtRLAOISCPQOW1Yr6q+b8FbZYwxxiyDWori54H3AZ8B/nzasBeKCEAWTvkSjjHGGHNOqaUo3gB8XFUfmWHYX015fPoXUTXGGLOo7j04yLGxAvceHGTbhpblbs5Zq5Yv2vwIeLmIfGWGu1nM+iUcY4wxZ4d7Dw7y0r+/l6MjeV769/dy78HB5W7SWauWovhSgnOJB4EvThu2DRgLH2+Y78JFJC0iXxCRz4rIz853emOMMdXt3D9Asexzbz8Uyz479w8sd5POWlVvHSUi/5Pg94hZ4F5VvWPKsL8BBoG/BW5V1bfNa+EibwOGVfV2Efmqqr5pYpjdOmp+LIuA5VBhWQQsB8gWPfb1Zdg7Ahc3wgXtdaRj7nI3a1mdya2j/p3gcm4vAK4B7pgy7PNUvoTzsdNoVzcwca7SO2nAmjXs+vF9c06885572P6CF57GYlceyyJgOVRYFgHLIdDy7BBdex9g/cXXsG1983I3Z3lFZy99tRTFm4AfAJ+eemFwEbkpfPjV8P/VBLeSmo8egsL4ENMP5SpQ9maYZNpIVcc5X1gWAcuhwrIIWA4A29Y0UDgcY9uaBsvjDIvi3cB7ARWRT6vqj8Pn24GfBu4FisCLgG/Ps2n/DtwqIq8Cbj95kILvzz21Un2c84VlEbAcKiyLgOVQYVlUVUtR/HXglwnumfgp4McAqvoFEblYVT8OICLzPkCtqllmu7tGLXuKaluAkyyLgOVQYVkELIcKy6KqWm8yPHHYVKY/LyL/RLCn2LtgrZpYZC1bNLbVU2FZBCyHCssiYDlUWBZzqqUofhr4J4KCeOvUAar6IRFpCB+PLmjLFPCqHT7V6uOcLyyLgOVQYVkELIcKy6KqWoriIeAxgjL17PSBC14MKzMGr4bd/FrGOV9YFgHLocKyCFgOFZbFnGopiv+L4FxiDPhr4M2L2aCT2J5i7SyLgOVQYVkELIcKy6KqWoriY6q6B0BEnpp4UkQ+zsnnGlVVP7BgLbM9xfmzLAKWQ4VlEbAcKiyLOdVSFN8oIjeH47oi8m+q+kamnV9kMS4IXvUnGTV+Ged8YFkELIcKyyIwjxzu7c2w88gY29fUs61rBV4Fx/pEVVWLoqrOeCkIVX1WRJ4PvA1IhU//0oK1rJavDtf0A//zhGURsBwqLItAjTnceyzLS29/mnxZSUSE792yhW2d6cVv31KyPlFVrT/JmM07gWME91n8hTNvzjRVt2hsq6fCsghYDhWWRaC2HHYeGaPoKQoUPWXnkTG2dSQXv3lLyvpENWdaFI8DCcAHOs68OVOoQrmGK9pUG+d8YVkELIcKyyJQYw7bO1LEHKHoKzFH2N6RWnn5WZ+oqmpRFJH3Apep6jtF5PdV9aNTBn8ZKAAfAL4/nwWLyPOA3wYOq+r7TxlBqeGEcI1fxjkvWBYBy6HCsgjUlsO2tjjfe8U6dh4bZ3tnim1t8RWYn/WJamrZU9wMHA4f108b9mpV/TPgN+e7YFX9iYh8EHjPLGNUf/NqKpznCcsiYDlUWBaBeeSwrTXOttZ48MdKzM76RFW1FEUFkiJyGcGdMKZ6oYh8CBgBUNVPzTQDEbkc+NNpT8/5pZye/j4ufH3l/ok3v/KV3PJTN500TqZcZOeRAzW8hJXPsghYDhWWRcByqLAsQkOzX5W0lqL4l8C7CL5l+nvThn18yuNZf5Khqo8AN09/XkQ2zDbNmuZWfvSPX5yzYT8aPs62xvY5xzlfWBYBy6HCsghYDhWWRSD+gufMOqyWotgANBEUvYZpw16gqn8CICJ/CNxVa6NE5ALgI8ClIvIrqvoPU4er+pSK+Tnnob5WHed8YVkELIcKyyJgOVRYFoH4HMNqvXXUB8P5/A3wVgAR+RpwiYhcFY5XnE+jVHUf8LNzDMcrFKrMw686zvnCsghYDhWWRcByqLAsqqulKD4NjAHjwORl3lT1DSLyAlX94WI0TH2fctU3T2sY53xhWQQshwrLImA5VFgW1dRSFF8MvBRwgNKUy7wBvBz4IYCIfFhV/3ihGqaqlPK5uceJOZSKc49zvrAsApZDhWURsBwqLIvqarnM2y0icmn4+LFpg5unPG5dyIbVsqeokbht9YQsi4DlUGFZBCyHCsuiulp+vP+3QH/4+FdVdepvEg+JyBcIrmjz+EI2TFUpVSuKqVjVcc4XlkXAcqiwLAKWQ4VlUV0th0/LE1exEZE/nzbsr4A1wBDz/KJNVb5PuVDl26earjrO+cKyCFgOFQuZxYP5GD/Jx3heosjViYX9qC826xMVlkV1tRTFiIh8JHw8/eq4fwOkVfUdIvL3wK8uVMNU/ernFP3q5x3PF5ZFwHKoWKgsHiom+OXBVooIMZTPtRzlqti5s2K1PlFhWVRXyznF985xTtEDng0fjyxkw9Sv/i0p1Vq+oXp+sCwClkPFQmVxX66eIoKPUALuG49ymS7ox31RWZ+osCyqm7Mohr9FVEDCv3XKN08huBj4xSLyHk7+0s0Zq21Psfo457qHvTS7/Xqe44xxhZuddbzzIYtaWA4VC5XFVf4gUdoo4xDB5yp/8JzK2PpEhWVR3ZxFUVXfMPFYwoo4bZQPAy8jKJqfXsiGqa94xSpbNFrDOOewR7Se9/gXUMIhis+tzmNcLmMzj7zCs6iZ5VCxQFlcSoFbncd4QBu4Rka51Bs7t64pbX2iwrKoqtqe4vtV9S9E5J3Aq0XkcVX94JRR3kpwTdMIwU8yvrRQDVP1q3/71K/+DdVz2f3SSckRfBFKKtxfSnOR9s847krPolaWQ8VCZnERBS4KvoROaUHmuHSsT1RYFtVVO6e4Lvz/eeHvFT8xbfgmVX0TgIh8knkURRH5ZeB5QBfw+6r60Ekj2GXeuNI5QTSxgZL6RFGuLJ7A82d+vSs9i1pZDhWWRcByqLAsqqtWFLtE5BXA8VnGbxCR6yeeF5FLVLWm3yuq6ueAz4nI1cAtwElF0fd9ijWcU6w2zrlsKzk+VrybPZE2riz3s7U8OOvvXlZ6FrWyHCosi4DlUGFZVCenniacMlDkEoJLuX1OVTMi8tOq+o0pwyd+qjHxZRxV1f89w3xmu5/iIPBZ4A9U9fDUgR2trdrQULkpxytuvJFXbN9+0gy0vh4Zm+Uc23nGsghYDhWWRcByqLAsAg0bNrBjx47dqnrt9GFzFsVqJvYMReTngadU9d55TBsFPgl8UlX3TB++pbVVP37TTadOOHUeL34x+v3vz7fZK5JlEbAcKiyLgOVQYVkEfuaLX0REZiyKtfx4fy4/LyJfATYC24GaiyLwMeAC4NdF5L9V9WtTB/pArspX3BJA/pz6GtzisSwClkOFZRGwHCosi+rOtCiuJ7iKzZ8CvzWfCVX1fXMN91XJ+/6c84jXMM75wrIIWA4VlkXAcqiwLKo706L4h0CjqvaIyO0L0aAJPpCr8ubV1zDO+cKyCFgOFZZFwHKosCyqO6OiqKpPTHm884xbc/K8KVR583yoOs75wrIIWA4VlkXAcqiwLKo70z3FReND1d18tUMBkyyLgOVQYVkELIcKy6K6s7Yo1vLm1VI4zxeWRcByqLAsApZDhWVR3VlbFH2gUOXnIlrDOOcLyyJgOVRYFgHLocKyqO6sLYpK9Z9k+KpVxzlfWBYBy6HCsghYDhWWRXVnbVGs5ScZdiigwrIIWA4VlkXAcqiwLKo7a4uiAsVqh09Vq45zvrAsApZDhWURsBwqLIvqztqiWMtufi1XvTlfWBYBy6HCsghYDhWWRXVnb1Gk+m7+zjvvtEMBIcsiYDlUWBYBy6HCsqjOWa4Fi8iLROTTIvJNEXnN9OET35IqqHK4r2/Gx3fdfffk44X4N3XeCzH+XMNnGjb9ufn8vZxZ1DLubOPU+vxcf58tOSxFn5jrta+ULM63PjHfLOa77jiXsliqdeZclm1PUVXvBu4WkWbgD4CvTx3uq05eeeFYXx+NLS2nPFYW9uoMU+e9EOPPNXymYdOfm8/fy5lFLePONk6tz8/199mSQy3jn2mfmP7c2donahnf+kT1cU6nT0z/+1zKYinXmbM5o1tH1WqO+ym+CngP8H5V/cG0acao7Mm6wDGgH2gERsLn28LnFsrUeS/E+HMNn2nY9Ofm8/dyZlHLuLONU+vzc/19tuRQy/hn2iemP3e29olaxrc+UX2c0+kT0/8+l7JYynXmelVtn76QJSmKcwnvq/g1VT3lEKoxxhizlJbt8KmIvBbYAaSALy1XO4wxxpgJy76naIwxxpwtlu3bp8YYY8zZZsUVRRF5voh8QETevtxtWW4i8ioR+YiI/Nxyt2U5icgNIvKvy92O5RK+/g+KyNuWuy3L7XzvCxNs3TC7s7ooisgmEflHEbkt/DstIl8Qkc+KyM/ONI2q3gfU/h3hc8RpZvH/gL8EVi9lWxfTaeZwF/DQkjZ0CdWQyXWq+jFWUD+YTbUsVnpfmFBDDitu3bBQzuqiqKrPqOo7pjz1WuA2VX0n8GoReY2IfD7899Yp0/0uUL/U7V1Mp5OFiDjA7wD/sAxNXhSn2ydWsmqZLFOzloVlEajhc7Li1g0L5ay9zNssuoFHwseeqn6daT/6F5FXA9cAvUvctqVWNQvgw0AzcD3w7SVs21KqpU9cCbxIRB5S1e8udQOXwUmZAD8Wkd8Fjixfk5bNSVmch31hwvQ+cT6sG07LuVYUewje3IeYZS9XVb8JfHMpG7VMasnio0vaouVRSw57CC4Ucb44KZPwkOFdy9ukZTM9i/OtL0yYnsP5sG44LWf14VMRaRWRzwBXi8iHgH8HXicinwZuX97WLS3LImA5nMoyqbAsApbD6bPfKRpjjDGhs3pP0RhjjFlKVhSNMcaYkBVFY4wxJmRF0RhjjAlZUTTGGGNCVhSNMcaYkBVFY4wxJmRF0RhjznMi8nYReSh8HBWRZ0Xk5kVYzm3LOX0tzrXLvBljjFkcT4rI9cAq4F4AEVkP/DYgwH6CC4j/HtAE7FHVz4rIPQRXzLkWeL+qHg2nTU0fF+gQkY8CFwK/BbyT4ELlj4rIv6rqm0XkEeDzwHOAdwDPBX4FeGq2+U5vAxAF/gg4AfwHkAXeTlDz7lXVL80Wgu0pGmOMAbgNeB3wCuCO8Ll3ATlgALgcUILCMgi8MRxnTFX/CvgX4MYp85tpXE9Vfx/4OPDzs7SjR1X/kqAwXwX8OvDLwGfmmO/0Nrwb+N+q+tuqeg/wPmAI6AOunisE21M0xhgDQfEDOAb44WMH+KKqPgwgIq8DHlfVfxaRH4TjZMP/S0B8yvxummFcnfK/AgUgIiICJGeZnw+Uw3Fnm+/0aWTKawCIAX+rqkPVQrCiaIwxZsIHCIrVxF7crcCfiEgvMAZ8EfgzEekC3CrzenCGcSMi8kfAVoK9t3bgN4DHCQrZTP4e+HOCYj3bfKf7FPC/wnZ/E/gY8Hcichw4qKp/N1uj7YLgxhhjTMjOKRpjjDEhK4rGGGNMyIqiMcYYE7KiaIwxxoSsKBpjjDEhK4rGGGNMyIqiMcYYE/r/cM8uUaIWHnMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 510.236x283.465 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "new = False #True\n",
    "\n",
    "N = 4\n",
    "\n",
    "params = {}\n",
    "\n",
    "steadystate = np.array([0.01, 2, 15, 100]).reshape([N,1]) #np.logspace(-1,2,N).reshape([N,1])\n",
    "\n",
    "# no interaction\n",
    "omega = np.zeros([N, N]); np.fill_diagonal(omega, -1)\n",
    "\n",
    "params['interaction_matrix'] = omega\n",
    "\n",
    "# no immigration\n",
    "params['immigration_rate'] = np.zeros([N, 1])\n",
    "\n",
    "# different growthrates determined by the steady state\n",
    "params['growthrate'] = - (omega).dot(steadystate)\n",
    "\n",
    "\n",
    "params['initial_condition'] = np.copy(steadystate) * np.random.normal(1,0.1,steadystate.shape)\n",
    "\n",
    "params['noise'] = 1e-1\n",
    "\n",
    "params['noise_linear'] = 1e-1\n",
    "params['noise_sqrt'] = 1e-2\n",
    "\n",
    "files_ts_impl = ['results/without_interactions/timeseries_autocorrelation1.csv']\n",
    "\n",
    "implementations = []\n",
    "\n",
    "def create_new_files():\n",
    "    for f in files_ts_impl:\n",
    "        if os.path.exists(f):\n",
    "            os.remove(f)\n",
    "        \n",
    "        Timeseries(params, f = f, noise_implementation = NOISE.LANGEVIN_LINEAR, \n",
    "                            dt = 0.01, tskip=4, T=50.0, seed=int(time.time()))\n",
    "\n",
    "if new:\n",
    "    create_new_files()\n",
    "    \n",
    "fig = plt.figure(figsize=(18/2.54,10/2.54))\n",
    "\n",
    "PlotTimeseriesComparison(files_ts_impl, composition = ['ts', 'nc'], fig=fig)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-19T15:26:46.070244Z",
     "start_time": "2020-02-19T15:26:44.965605Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAC1CAYAAAC6Xo9GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXwbxdnHf7MrrVa3LMv3mfskdwKBAIEESCCEu9wtR4FytKUlQFugXKWl0NKWAm+hQGmBclMg4SYQEkIuEiCQkDu+D9mWrVu70u68f6wky4pkS7ZsJ2G/n48TaTU7M/to9cwzzzzzLKGUQkVFRUXl8IQZ7g6oqKioqAweqpJXUVFROYxRlbyKiorKYYyq5FVUVFQOY1Qlr6KionIYoyp5FRUVlcMYVcmrqKioHMaoSl5FRUXlMGbYlTwhRDPcfVBRUVE5XBlWBUsIuRfACACXEEL+Qin9RX/qcTgctLq6Oqd9OxzZvHlzO6W0INVnqgwzozcZAqocM0GVYW7oS44xhtuKNgPYFX0d7m8l1dXV+OKLL3LTo8MYQkhtus9UGWZGbzIEVDlmQqYy3NzQhWaPgCUTi4aqa4cUfckxxnC7ayiAEkLIEgDFvRUkhMwnhKwhhPyDEDJ/SHqnoqIybNR1BvFlo3u4u3HIM9xK/k4AXwOoAHB9H2UpAB8AHkDDIPdLRUVlmNEwBJKsJlAcKMPtrjmTUvoPACCEXAjghV7KrqGUfkoIKQLwEICLYx+0tbVh1qxZ8YJXX301rr766kHq8uGJKsPcoMpx4MRk6LaPht9SgdLaMaoMB8BwK/kjEl5PQS9KnlIqR192AtAlflZQUKD6QQeIKsPccDjJMSBGwGtYMAwZ0nZjMnx/pxOr9nbg6lMnDGn7hxvD7a6xEEKuJIRcAcDeW0FCyNmEkMcBPAvgkSHpnYrK95i8O97Ds5uHzzOqumtyw3Bb8tcBODnhdVoopa8DeH3Qe6SikkAoLIHXssPdjWEhLFE4fcKwtc8Sgoiq5AfMcFvyJwA4A8APAPxzqBq98LnNOO/fh8eU+vuIJFMExMigt+MNRWD49TsDrmfNvo5BjRLxCxFsa/EOSt0DVbLMsuXwhPoXHa1hFUt+f0cATe7QgPoxmNS4Anh6Y138fUSSeyk99Ay3kj8TwDIAPwFw7VA12uwJ4bVvmhFO+DLUxyAeOvxtzT5U/e6jQW/HHVVO7mC/t3BAlimufPlr3PjGt7nq1gH8c0MdjvjTqkGpW4r+LgaiuELh/p0bs+QXP7keFz+/JaNzZJli3P0f96u9TBEjMtoSZjivbW3Gj1/+Ov6eu/VtXPrfzPo7FAy3km8EoAfARv8GHU8ojLBEcewIO5o9inXgFyJgb14Bzc3Lh6ILhxy723zD3YUeLFu+HR2BME56fN2gtuMJKbOFgbgsJv9pFfa0++EVBm/mEYpIg1Z3RFKUPHfr2/CGsruGDr8IAAjL/VPyGoaBJFOYdZr4gNsXEqXY3e7HjIc+xb4Of4/PfEIEy7e19Dh2/KNr8e53rVn168+f7kXx3R/E38spDMTntzRmXF+zJwRm2eDpnuFW8mMB/BHAoxiixdSjHv4M62o7MbrAiI92t8MbisB827sAAJkqPtjekGWKPe3+XsvUuALY2uTJWZ+HEzEiY9wfP8n6PGbZcjy9sS6nA0SHX0RDVzD+fuXu9rQLcx/tasPVr3yd8rNM2ecKYP6ofFTm6ftdxw6ncv3l1uzqWFfjwu3v7oi/FyISXvwyteIYW2AC0HM2+nWTG1saugAAN721LeV5dZ0BOH77Xvx9fVcQ37X2dPtIlEKIDiLZDlQFd74PYACWPKO0zxKCTL1Gsdn5V00ejP5Dt0UvyxQrtrfijH9t6lF+zX4XtjZn5+q67d0dSNTrUh9egD+s3H3AMVmm8e/rw11tAJR7FgBavbldBxlWJU8pvRzArQDuiv4NOjoNA+nBJVg6sRh1nUEsW7EN1Xl6eO9bjAumlWK/K9Dr+av2dmBsH9PBm97ahmkPfZrLbg8pK3e3QY7+qlwBxRoTI33/UDc3dOHt7d1W0Y9f/jrtABGRZLwVtapW7WnvU6braztRcOf7GJNUrtEdTFl+5Z52PLmh7gCllQ1Ln96IVXs7oNMMfJKZzUAhRmR8tt+F30eVg0+I4O4PduGiNC6LmLX9eU0n7nx/JwJiBNe8shWXv/QVwpKMv6ze16O8LFNsb/Gi+r6VcAXC8e+66ncf4bhH1+KlrxqxIvo97nL6oP+Vsi7xwa421HcF44qUWba8x/f9j89rUvYvlMG9kwoNw0CSKBJV6NdN6dc2OvwiRCm1wtXcsiIuv4AYwZvftuC/W5TIoZiSnvTAJ1nPVgCkHYDa/QJe+LIRtyUM1gBwzN8/w+lPb8TP31QG38te/AoAcPIT67H0qY0oSZgl5IJhVfKEkAehRMzcCeD5XNT59vZWdAbE+I2bSIdfRJ5eC0IIyqw8nt/SgH+ur0OFTQ+jToMxBSa8ua0FZzy9Ma1y+HhPO4DeffiaaFzxW0lTw0OFkx5fH1cMHQFlmvxY0g94zl9X45Wvm9CV4K+e/dc1OP3pjXgnafqbyp/b6A7hzH9tAqUUrmAYe9r9cAXEtHJduVuxcoSowth283ycNNYR718ysR/rpAdXHfAZs2w5PtzV1qef+dTxhXjqB1N7LZMpibHmbT4hfp2STOOznqI738d3rV7wv3o7HtETECP41dvf4f6P98TPv3XFdty8fHv8fcxd0+wJ4d4Pd8H0m3dBos2t2tsBAPHfw9vbW6G5ZQUmJ/jwN9Z3wXa7MpvtCkVw/evfYOnTGwEAQoKMnD4BVb/7CM9uboj3/6WvGvHfLQ145LP9uO71b1Jee3/dSSxDIFGKTfVdIAQ49pHPMP2h1WCWLcdfVu9FozuIWlcA9V1B/GHlbhTc+T5m/mV1jzo21nX28J8DwMQHVuGsZzbhkv9+CQDxGc93Th+2tXrxwU4nKKVxF8qja/fj2te24q1tLXGrOwalFDud3bPVc//dPVMovPOD+FpC4n29rrYTmxu68F2rF79MmmWtyNJ1lAnDHULJAFhFKb2DEHLTQCsLSzJufPNb7O0I4HeLx+M3C8aAUgoSveOPengNjh/lAABU5emxt0Ox2t+96kgAQJmVx09e3QoAWL69Fb7fL4aB6xYRs2w5OJZBVZ4er3/TjHOmlKbsxytbmwEAv/9oN5ZOKsaja/djxfZW5Om1+O8lMwd6mYPOPYvGYV2tC8yy7fj0uqMBdP8QYnzR4Mb5z24GAMh/Oh3BBDfXkqc29ii7q82PicXmHsc6ojME9uYV8WOO374PXsMgcP9pPcq+8W0zXv66qcexCUVmTCq2pI3c0Gl6t19OeWI9PrpmLk4c4+hxfF+HHyPzjQAUC+/UCblJjiVGZDR7QvjPFw349TvfYcPPjkW1XY/aTmUmcvcHO9HmF+OD0ldRi9UVCKPF2x1ZcvYzm/DGt4rx8ODpEwF0D3zr6zrj5TbUdWFysRnXRF1WrT4BJRYepz/d87sBlJlCbP3BkBQu6vSJ8dds9HfkFyXURftNoViivUXh9Mc6BhRjKVbv10nuz5ve2o6b3tp+wDnJM/GjHv4Mt544usexuq6es7/Xv2nBqf9cDwD47Xs78NHudnxy7VwAQPXvPoqXf2tbC5o93QPGnnZ/jxnotD9/iq3Nqd20oiRDp2GxMfodhSWKhq4gNjeknpkwy5ZDenAJCCEIhSV4hQgKTLqUZftiuH3y+wC0EUKeAjBg7adlGbx31VEAEPdnnvbkBpz1r42498Nd2NsRwG0LxgAACkw67PrVifjfZbPjinx6qbVHfTsSRuiYW0CUZEwpseC8/2w+wOrc1uKFX4hgVL4BH10zF2ZeA2bZcvz0f9/i/Z1tePGrnorqYOXIyjy8/o2iSJqii9PPbWmMr1ckWi6AYqX89H89o0euOqoy/vqov68BoKx3bKhVbvJZf12Tsu1QREZt9If6xPpa7HT6cNs7O/BNgt90560nAAAsOk1cOcX6EbPOLbwGP5lbBQA9BiBngr9z4ePr4BciyItasU6v0MOP2+YTUWDkUvYzW55YX4sznt6IX7/zHQDg+S8bUHTXB5jzN0UO9V09QwSboyGDgbCEYIJPO6bgE7kmapg89GlPt0xnMIyaqDJeubs97eJeYpSZgWN7+NDX13YPHDevUJSqGJEx4vcrlevY0ohk9b65oauHYj/pifUp2+0LNkebobZnEF763s6oX3y3MlM/4f+URf3EAaErKcrqn+t7JoFMp+ABxF1eRz38GQBFj+xs8x9QZyKX/vdLMMuW476Vu1F0V/9dOMPtk3+UUvowlDDKi/sqnwmjHEaEopbgPz6vwbraTry5rRV3vr8TV8ypxIh8Q7zsaIcRZ0zuTn45u9IG+U+no/b2hbhgWikaE2JzP9nTgXsWjcO4AiPeuHw2xheacHX0xwUoo/oRf1oF823vYrTDiGNH2rEyesPE2gLQ75jh3hAjMmr6WEvIhnIrH3+9qb4Ly+aPAgD88RPFZdCQ5Ad3hyI94oQBoMLW7YMujFog7+xwYu7fP+szXHVLNKb8J69uxaOf14DX9rxNx0QXGi28Bp0JP5I1+1zgbn0bgOIOuny2MtAYE2LdL32h26+9eHwhFj+5Ae6oQnIFFav1z6v2xsvkckv/FwlW28Nr9vdaNuaGOv6xz/HuDmfKMje++S12ONMrsMT7d3cvwQKnJcy8Wr0CAtFBcUz0nk3m8STllqyIZ/91DZ75oj7t55miYUiPQaa/CP1cE0gmmLSA/GDCfZIJc/7a7Uryi327sP4bXWi/76MDF26zYdiUPCHkFULI54SQTwG8DGBDrurmNAxmV9hw3evfxH/Ap44vxK+Tpm3pqLDpMX+0I/4j2d7ixV0f7MSschu+u/VEEEJw4fQyPLWhLq6wEuN47100HlqWAccyeOPy2XjhkhnY9asTAQAb67oObLAfSDJFiyeEf2+qh+k372Dk71f2UE4DYYS9eyBcu9+FeSOUjBN3f6Ck/l+1twMnju52c6zZp/h9X/vRrLjrq9Kmx5VHVuLuU8bhlHGFALoXcZ/Z1K0AAOD0iUXYfsv8+HuZ0rgPeXN9F75sVCwk/x9ORetdJ8fLlVn5+KIVADRGZx1Or4AnN9RhtKP7OnZFo3yml9kAAFZeA68QiYfR7m33wx1U7pWYxTpQ1u539fvcTfVdWDDG0WukxcNr9mPiA6syqu/eD3f1XSiJdONbb9FlMYWeeOqSfrq8WIYc4FrpDx8k+dGHiy/SuGYGm2FT8pTS8wB8TCk9nlJ6EoD/5rL+64+pBgCcO6UE0oNLsOLHR2JUGsskFRVWHtuii6//2lSPfR0BTErwK99ygmLdsjevwClPrMOm+i68dcUcSA8uwawKRZEE7z8VSycV4/xpZQCAB5dM7LFpYiDc//FulN7zIS5/SfGHHlOdh5tXbM/Jpi5ey6L29oX4xXEjsaGuC+MLTT0+/7LBjYfPnAyTTvHf/i/qQtAwBKeMK0TTb0/CxTPK8c/zpuIHU0vjrq6rX1FmPhpWUQEVNh5HlJjx0qUzkadX3CI3HT8KXkFCm1/ElBJL3Ge95vpjoNeyPfyS508rw4wyxcXmCoh4L2rxzv27MiXOM3S7Wh5dWwMAeOCTPThupB2Nvz0JfjGCeSPsOH1iEXa2+fDgqu7FzWBYii+g95djH10LALh30bh+nT+xyJzy+MKkdYS+uOn4Uf1qnyHZX38s3PJnCZu/+rvwmih/vTa3qqrCxvdd6CDjng+yH6iB4ffJjyaEHE8ImQug11RzhBAjIeTfhJB/EkL6dO38cFYF3L9bjJd/OCu+8JoNM8tt2N2mWCzBsITHz53Sw42h07D4QzQ73oe72lGdp8eCMY4ebSW3e+3RVajrCuL6179Bu7/bQnt07X78a2MdmGXL8d4OJz7e3Y6j//4Z6joDeOPbZnzd5AazbHmPrd13vLcTAHDFnEo033ky1twwD5U2Pa5PE+GQLRU2fTzsb3S+EV/ceCwAJZb6nR1OjC80IT+qRJ/ZVI/ZFTYcOzIfAFBs4cFGf6B2gxZvbWvF7jZf3L+9ep9i4d63eAKWzR8FXsvCbtDinCNKMLnYjK3NHqze14FppRY0eUK4Zm4VjhmRPn9drSsAx2/fx7ObG2DQsj0W3845ogQlFh3+/lm3eyQUlmHgNNje6kOjOwSHicNFz2+Jr0MAiuVv4fsfl5Do575t4VjU3r4QAFBkPnDx7PxpqRfwbXptj/cFRg6FptRrBKN7MWDqOrvlMbO8e92pzJpe0U0uNh/gaweAUQnuzmSmlVpg+s27PY7xGqbf7hI24fdzwzEj4q+LU8gwW+aP6h4oqwawDyIdG352bM7rvOuDnf06b7iV/E+hpBueCeC2PsqeDeBVSulVAJZmUrl5AD/SfCOHzQ1deGJ9LR77vAbzR+UfoLRvPXE09v1mAervWIh9ty2Evo9EVgZOgw+vOQr/93kNCu/sXkh5bG0Nroxa+Ne8+jUWPr4O62s7UX3fSpz9zBeY/pDiy4tFT2yq60KFjYf8p9Px5A+mxhXHsxdNx8d72uMRQgOlzKIoAYYhmFGuzE6qfvcRfjizHAxDsOrao+Nl/3fZ7AOUEgDkRY9d+9o3aIvugHxqQx1umDcCl8wsx6UzKwAoi+av/GgWKIC/rt6H85/dDJtBObc3i1LLEnQlrHMw0Tv6jEmKi+CVH83C8dHBJ5av6IojFV+9EJHx0e52/Pn0SXBEB6DwA0vw4yMrMf2h1X1+n72R7HONrVEcUXygdR6bfLEMQb5BG48MSox0qb19IXb96kRwbPqfbF6C/I1c97lcQqRRzM/+tzMn472rjorLBgDyDdq4HEIROeWic28K+6sUGwAfWDIRa/a7+rVmxCZY8jG36wuXzOixkN4bPz92RNrPvklYJF06qXtd7kezyvusN91Am8iEIhPuOGlsys+mlVr6PD/micgFw63kiwCYAOSjjyyUAMoBxJy5Pb7l2EMGYn9PPPHEgDvGMgRLJxfHFWa5LfVoX203oCyL3YwLxhSg895FGJlviCfZyjdyePnSmai7fSHqu0L42bEjELz/VADAw2dOhk7D4F/nT8Nd7+8Es2w5jnx4TdwFlMjR1Xbs7QjgifW1fe7cTSaVDM+ZUoK2u085oOx/oulnq+wGyH86HWdOLkaJJbV1pYkqpdj+gnEFipJ5Pk0K20SlI0Rk1N6+EA8uST/Jy9Nrcds73ZtNfIJy3dfMrY4fe/r8aTBoWbz2jRLaelZ0sf2SGYoMbXot9kXDaVmG4MkNPReRsyEmx+MWnAyNcOCiaGx3KgDccoKyRiTJFFfMqVSScd22EH88TbnexO3yFTY9rHotHjtnCu5bPAH/vXhG/LM/L1VCKaeWWmDWKYaNQcticnRAMSQo/NiA+dN5IzCp2IxPoiGyv1s8Hm33LIorxltPGI2fHzcSp4wrQIWNR81tC7Diyjn44ayKA67pvsXj08ojpqx2ZrHzOSbDBQsXxI/dvnAMXr50Js6fVhZX+DGOrs4DACTbAmcfUZK2DSEi486TFSWcqLSPS7j/EiPEYtx64mg033lyj2N/P2sy5D+d3sPoMek0uPuUnm66Zy+ajmXzR+GzG445oN7kmViqEOBfL8hsTTGZ4Y6T/yWUpzxlEnLSAEXRf4WkwWmwHtRQnWdIGbc9UKx6LY4dYcem+i6c95/NaPeLWH298sXv/fUCFFt00GlY7P/NApRaedwwbwQkmeLyl7oXGX8670ArhWUIxD+ehiP/tiaePTEWa9sX6WSYn0EI4euXze718+lllvji6cafHwfr7e/iJ0dXpSw7It8AC6+ERuq1bI8onVTY9Fq8+FUTzppcjNcumx0PE1w0vjBehteyuGhGWVx5W3nF4u1ME772i+NG4i+r9x2w8SUTYnLc1ebD+c9u7hHfHfjDqXj7u9b4xrJbTxiFBz7ZA4nSuIIy6TTQMOltr9hDrWdX2jC7wobLXvwS186txvs72lBs1uGhpZPQ5hehZQjmj3YoIXiLx+Of65VrT3crxPzfsR2jV0ZnO4+vqwFLCCrzDKjMMxyQ2uCoqjxMKDKh0qY/YJE0dP9p8Xvvw11t8QX4vki8Fx9fV4N3dzhRbtPj3Oi98PKlMzG5xBxfdI4NXMnLUfNG2BF5YAk0t6yAw8ih3d8d839kZR60UQMkUd6lVh5/OWMSfvHmNkwvteI3C8bg9yt3o+OeU9AZDKPQpAMhBNV5etR0BvHrBaNx3dHVAIDjRuUjFZ9cOxdalsHR1XZcPKOnKw8AOu9dBFcgjFF/WBk/xiftspb/dHoGkkvNcCv5bymlmabnex3AI4SQ0wAMSSaxG46pxqUz+56+9YcFYwpwwv+tw5QSS4+FscQQz6qEKBeWIXj83CmYW5WHySXpp3uEEPztzMk45pG18fe54oczy/GfzQ348Jqjsjrv72cdgXmPrMVXvzw+7kK76+T0i5EGLQsrr4nvaegNa9RFMTO62N3425N6/JhjzBthx5Mb6tBwx0lx90ViSJz4x9Pgi7pYLp1Zjr+s3odroz/e/uAXJYwvNGF6WbcPnNeyOGdKKX51ohv3f7wn7t6SKcUVsyswJ3oNsa8sprwSI4oSGeUwYs0N8wAAb10xB4QgrrhiPHPBNNj1HMY4jNjd7u8RPJBI7C6ZXGzGonEF8eM3zx8NT4JiP29KKU65txA3Ld+Gf2+qx+c/Vdo/c3LJAXH4iW6iZKWVKdfMre4xKwOAc6f2XMNg09zjhBAQoszY7jplXHwPxLVHV+PRs49AQIzgwmlleGWrsn8lpkhPGVeIX7y5DSxD4hFGeQaux0L+vtsWQozI0DCkx2/Mc99iWG7rXpdIpZyTF/Stei20bM9jF88oi6e1iO0L6S/DreRPIITMByAAoJTSH6QrSCn1A7h8qDoGHPjF5pKLZ5ThtAmFWdV/1VGprd9k5lbbBzTyp+OZC6fjP5sbDoi26YuYbzkW0fDOj488QBklQgFMKbFkNIt4MxrZc0pUMZVYeJRYDlxQtEflXJqw2Hj/aRPiOVo0LAObXunTtKhiTpVdMFMCogSHkcPfzzrigM/uWzwe08usIITg7SvnwG7gcGRVHuZWK4vLFTY97AYtOI3yw89kpyOXZodvzL3yxY3HQadhoGUJrj86va/6vKmlOC9BiS4cW9Djc4YhsOm1oHLm0Tf/u2x2Rn7o/jC3Kg+Tis0IhiX8/awjcP3rW/FFgxuPnd0t9/9cNKNHRsqHz5wMQFkjG5GvwZVzKnFUZd4BdXMs0+s1ppJ54lpIOhIHhdjv1MBp4Lp3ETQMAaU967EPUAcNq5KnlOZeEx0iEEIGbQAZTN68fDZKzNmFnx1RYsGzF02PX2+iKyUVP503AqPzMwt3rbTpEYrImBldGE6HNcUi/KwKWzzcNZmfHTsCC8cUpPwsE/xiBEYu9c+LEBJXpItTxJAvmViE9nsW4Z3vWg9IM9BfEoMQUgUksFmGi44vMvXwXwPA5l8cd0DuGAA9Nhzmmnkj7Lg/uoZBCMHGG48Ds2z5AW6pWHqEkfmGA64138ildLUYOTateysdhBC8f3XfM91URliqwAUg/X6FTBlWJU8IeQGK4WYCUEkpnTac/VHpm9MnZf+DZRmCi2dk7vb6TQZumhgf/WRuRvHslX349pP56xmTsyqfjE+UMrLqeuPUCUXw/n7xgOrIhBllVsytOtCS7Y1bThgdXziOMb3Mij2/PhFmnaZHjpfBRKapXZKxfRcxYpu0tt08P6N6r5hTiTmVNnznzD5V9klj+28cpGKgmR2G25K/MPaaEHLjcPZF5dDEpMvsFo5FAg0VfjES3yw2EHK5ppKOL35xXM7qiiV3628yrWxJ5VJz/27xAbKPuQczTRv9ZDT7aA6zWmTNR9fMRUdATGvhZ8pwW/KxsEktgFnD2RcVlVziF6S07hqV3PDY2UekdLelckeNdhjhundR1m30Z9dvrkjOkNpfhvsu3AbFXSMAGHhwu4rKQYJflFCcZu+ASm74SZbRT/2xiMcXmjB1kBaNh4rh3gylAXA9gBsBHN1H2YzIxUao73tbQ93eULU1lNfkEyNYu2pl3wVzxOH4fR0MbZ11RAm+/OXxQ9LWYDHcSv4SABdG/89JquHhvikOh7aGur3DUclTCqx8Z+geDH84fl9qW7lhOFMNG6DsYi0DUALg0HiihopKBtx1yjiYu2qGuxsqKiC5SE3br4YJ+RcQT3THAWATo22yrKsNQOxJBg4A7b0UzyWHWltVlNKU8V1JMsxVe5kyVG0NqgwB9V7MEFWGQ3AvxhhOJc8DOBfA6QDGArg0ixQHKioqKioZMJw++VjqwIsArFcVvIqKikruGU4lfxaU2Ph/ARhFCFHjzVRUVFRyzLC5a+IdUJT7uQDOopSeO8C6jAAeAyACWEUpfT4HXQQh5EwApwEoBPAolAedjICyiesnUBaOH4SS5/5flNJPBtieEcBqAHcCGDeYbaVpW5XhwNvOuQyjdQ+ZHFUZHvr3IgCAUnrY/AG4FMDp0dcvDUL9eVBmHs9H398A4FgAd0C5URgA/81BO/cAuBXKE7AGtS1VhoeeDIdKjqoMD/17kVI67HHyuSbt06NyxO0AngQQe5pEbbTNcgD1lNL+PcwyAULIQgDbAbQCsA5mW2lQZThwBluGwCDLUZXhYXMvHnZKPvb0KCCH10YU/gjgXQCboIQ/AUBltM0GAOWEkFy0eQKAo6AsSF8EZUo5WG2lQpXhwBkUGQJDKkdVhofHvTj8PvlcEvV9PQIgBOAzmjt/8s8A/AjKDfEVAAOAKgA6ANdC8avdDyAC4DlK6cc5aPMyKHG0Ywe7raR2VRkOvN1BkWG07iGVoyrDQ/teBA4zJa+ioqKi0pOcZqEkhGgATAdgAbCLUlrfxykqKioqKoNIzix5QsjdUHxo+wB0AhgFIB/A3yilX+akERUVFRWVrMilki+ilLamOF5AKW1LdY6KioqKyuCSsxVdSmkrIeQ/hJDFScdVBa+ioqIyTOR04ZUQwgG4AMCpANYBeJJS6s9ZAyoqKiqHOIQQPYC/ACiFsiHrOwA7KVnsA5cAACAASURBVKV/Hoz2cv34v3wAIwF4ALQAeAqK0ldRUVFRAUApDQL4CSFkPoDJAFYAuIEQUg3gGQAbAegBNAM4EsCvAAQA3ASAANhLKf1rpu3lWskvA/AopXQfABBC1OgaFRUVlczZRSm9hRDyPwC3AVgFYBGAYgDB6N8R2VSYayW/N0HB35jNaKOioqKiAk/0f4FS6iGEiFA2TjEAnqWUbs22wpwpeULIgwDmRqccBMBoAKqSV1FRURk4jwD4PSGkGYCXUnp3pifmMoSyCsA8AGugPNavhVIazknlKioqKir9ok8lTwghAH5KKX24j3I3QnkoN4ViyVNK6S256qiKioqKSvZkZMkTQp4F8A4ANwBQSt9JUaYEygO541BKa5PL9dFOj2T+lNIPosefgZLAJwLg55RSIZt6VVRUVL6vZLoZ6iMoTzIpQHdazh5QSpsBzIcSOnk/gMuz7Qyl9A1K6VUALgNwfsJHQSgzhC4AqgtIRUVFJUMyXXi1AZhMKb2KEHJHL+UmAzgRwH+i//eX26E8livG9ZRSOZoidAmAtxILG7Vawc7zQQoCiTCEpTIlULNrJtPg84UppQWpPjPpdEK+wRgc6j4datR1daaVIQBYTSah2O5Q5dgLu+pre5eh1SoUFxerMuyDXbt29SrHGJkq+VHofkqLuZdyFihW+I+gKOOsiPr/7wfwLqV0S+x4wlNTnABMyedxLCsBgGSyk1DldO0CqaPm2sqivdm2f7iz8H+vlqT7LCbDGKeOm1D749lH1g1+rw4tZj3yl7QyBACe0/WQ4yWnnFp7x2VXqXJMgMyb0bsMeb6nDC+5pPaOO+5QZZhE1EXeJ5kqeQpATwiZDGUrbupClF6T0IGfZVh3Ij8FsBCAlRAyGsAxlNJLCSF/hrIDLA/Aj5NPMmm14jMnLVojgiGtWpOuKMwJwKA9TeuwxKzjxdcvuWzNcPfjUMdmMos7X/ifKscBYLPZxJ07d6oyzBGZKvk/A7gOysN1f5OuECHkIQBFUHz9FMrjrjImGsGTGMXzj+jxmzI5n4NMK8KeUDZtAkD34OATOMiqnyeHiCCkieh1DhoS2wnPldKgwEF9Uo2KylCR6cLrVAB/oJTeSintbdrUTCm9mFJ6IaU0KwU/VIhgSL3WwotgSOxYq9ak+9RQVdiqNekyPSfVMZUDaSJ63UriKNxKbOaVxFFYSwx8DTHwIkhKuYkgpLfPVdITkijZ4Q3zIYmqslOJk6mS1wG4O5pK+E+9lJtHCLmZEHIdIeS6HPQv56RS6EVhn3B8oNZZFPb1CM2MKfIGrYVPPqevgUFFoZQGhQW03TmFdnkX0HYnQLGSOAqbiD6l3GKDQmww8IFhVKWfGTWBiO7VxkBhTSASl62q+FUyVfIcAAFAK7oXYFPxEJQMatsAfDuwrg0OqRR6zM2T7KqJKXIASD4n3cDwfSAba5sDpdU0EDJBlqtpIFRFg6EFtN3poCExVR2lNCgcRzva6qHjX2VKyrYQm6W3QUGlm2qDRji3zOCsNmji92RM8e/yhfn+KvvDbaAIRSSyw+njQxHpkLmegfQ5UyV/QrTsbgDv9lKuHMASSumnyDJT2lCRTqGnIqbIy8OeUOI5ffnwD3dXTsza7o/ijSn9dsJzsTqkqIs+5r8HKDYTmz0AVuOAKC6g7c5SGhRUV043iYo39hoAxpu1IQCIfRZT/AAQs/JDEiVb3aJ+q1vUJ57vFmUmWZmHJEpWtYUsL9Z3zxAOdaVf4wrqXv26qbDGFRyw4ZCJ8u2Pgg5FJLK1yaP/or7TsLXZrd/V5uP72+dMF15/CeBYAFcAuB6Kjz4Vc6GEOQJAdbadOdhIXsiNKfcwGPK5oaLg+ECtM3mh1w8Ns9pU6WhiLfopYmvX1GCr51BbzI0p23SLpKU0KCyAonhTnQMoA0Fvi6ylNCgch462MEB+2YxprRIxTNCzLq3ZLJ+j89afidYmUIoSGhKaiV63hxj07dBx24jZMhVe93jq9u0gVtMU2uU1QT7oQqlCEiU1gYiu2qAReJZk9P1nc07MQo8p8Ls37BlViJBwQkW+K89iDX/SRfOOytO55xfwnvFmbSgkUbK0BG2+iMy80xKwvdUcKrRqmciV1campqDEfdoeslo0jOyNUHaunXNXGTVClV4jrHMJ5s86BOu8fJ07NkPY5QvzT9X4ShcW8h0nFerdPEtof653uKi264Vzp5Y6q+36jGfhoYhEalxBXbVdL/AaNn59sQHj3KmlzvGFplCqcjWuoO7FLxsLj6rKc8+ttnl3OwN8szfIHVmV523xiBwIRZ5eG35nuzN/XJHRz7EM/bKhy/R5TWdeSJLZQoMuNK3M6ls0vqA9mz7HyFTJ3wPgMwDXUEo7eikXAQBCiBVK/uPDipj75uhAfVs6V812fYF5lbG6dFqopf07rsBaHPEL/Yn4GU5ilvoCtDuraeCAvses8XTnAEBv58fq0ILSlcRReH4Zv39ThM+Ht41u62jLv0WU5lRx8Byph3OsWe9dwxUV+IhGo4cUGSv7vVuI2boRVvsu1my9QCI1lQiG0g0ofQ1Yg0WiEo5Z14mkUorJ58TKlOhYsVmQuMSyya6Zs6vzWr5tbDW9tK2mZJfLY3EGRP4pjZZOc1g6ZxfbuqYU5XntVlv4DadU6BIlDQAc7+Bce3xh/kNnyM4zRP68PWgfY9b632kJFXAs5CkWzrvHH+bLdawYiFBGkCiJte8Oy5qXGwJFYZmSSoNGqAtEdCudgn1R8dAq/nTKt6+y4wtNff4mE8vvcvr5ZzbVlV42u7JpSqklvlErNmCUWDjxi7ouw4a6DvN3Tr/56qOqG6eUWoKhiER8oQhTYuFCn+xps62vcVnX1bryfIKkXTKpuKXBHdR5AmENw0L+ZHd7yYh8o9dh5MRaV8DAamhEDENjLmWFf6zbX3WGp6R5TIGxsa/rTCZTJf8JgEsAnEcIeYFS+maacs9Asfr/AeCP2XTkYKE3V0yiH56DTGUKPBK2jZnFhtpnMaFODQEmBtu85wL7xwQ7fF4tr7WHA2K91sIfSuGZqSz1bM/J5PzYOQ4aEks1vFBqiwicje6XKLAlCOtqPwreaQ5W+ul+7Wie6RxtZL3H6MOutWyR/V22oPRIqbPdgZD4LnEUlkIMHUvbOjxEp42FazpoSNxKbOZvidl6MtrSDjiDQSr/eCKpBoHkc2JlZtk49xddojWxLM8Smjh4nDe6uOO80cVxAywkUbK21WPe7nQZ29we7sVtNaW7XR5zRyjMGTltpNpm9n/SZclvgEEf5oww6jhxlEnj3+uLmM4t1TdrWCLv8oX1ekLoWy2hIp8ka9sEqX6uQ+cxMCRi1hCxxi+ZXmsMFgVlmZVkCn+EsrV+Qb/BJZgXFug7XWFJ+2m7aL+y2tg0xcoNyg7WvizpxPfJZTOpO2aBO4zasCRT7Gn38iAUVXl6odkjctV2vVBt1wur9nRYnt9SX7y1yWMbmW/07Wn38j4xzHzZ0GXaUNeVZ+Q0YaNOI7+xvdkREiTNgrGFrTtbvbwgRcie9oCpwxfS6bWasF2vEUbmGzwOAxdaV9tu39UatDV3+viQTLQvfd1YMbPC5j19UnFXNjLKVMkvoZT+AAAIIf8H4AAlTwg5Nfrypej/pQC+yqYzhBAjgMcAiABWUUqfjx6fDODX0WJ/oJQO2qJuzFpP5YpJdt/IAEYxomdFxFT5MLVPLWNlz3zibTwh0NJmIFS2hUWpXmvhPzZUF44V2z32SEjUQqZFYZ/g0hq4g1Xxp7LUsz0nk/MTz+FoUIhb3AR0tgHuqQbiaSL6BlMkEPncH8lf640UvNhOxuSxzaEKg9tTbKTBck4IFTKW0FtsUUWbxHEyIcxkeN3fEpM1H2FhDzGYZlO3KwwQEYQMlTWfrISTSTUIJJ8TK1OiY8Vqo0aoNmiEROsYUAaCVJYyzxJaZjaI6zywnTuizMkxpDV2zsY2r6m1s5P7qqXT3NXSzDV0eQ3toTDPc1pphM3sXR+w5IU5I6rttuCWAMxNwQh/jEPX9r4zVLCxSzB7wrJ2fYdYMMGsdc8p0Ha93xrIbxVkg4YQeWdINqxz+QvfbQl5izgE3RFoK/VMoEqvEawck3O3WrLrJVmRJ77vzU2TakZQbdcLR1XluT+rcVmPrLB5Fo51dHy4uz3/070uTCm1euvdQf6CaWVOAPhsv8s60m4MhCMyGVNgDLzxbUuROxRmdzp9llIrHzhuhMPvEyPM0VX57fWdAVOhmQu9ur+9vLkzaGQYIjME0NGIFBAl8uwXDSMNLJUYQiSWQHaGYAAokxeWBDPPRrKVUaZKXk8IqYy+NqYpUwDgDCgP8Bah+PAPyFbZB2cDeJVSupwQ8hKA56PHfw5lLYACeADANWnOHzDZRM1oCLBYE2idrNe43zJWlrWLov4Lj1z4imAaaydycA4Tap0nB50TNW3ujXyZw89qNQDFxFB7Zxtr5CeLTvfUYKsHABq0Fj4IlvFoeO2UYKvHiMhB52ceTFK5iOLHNO3ORZag02HNC40iZs+EUKtnoz/seLyDnfRwGGyJrt0XNhPsNBosR8HXXkG9gW9gsn7G2AvNiIQ1VKafEEfBZPjcs2inhwOlqTZpxdocCtdOX4NAchkrx8QXVBN98cmzgcRBoNqgEZaW6NtqAxHd5s6w5YIK5ZzVXTTv3PJS5xHlpX59Y0BeWqJvAwCnL6Dd2tJpeq2mrdTnbNR+sG1nmRAOszzHSetbTKST0WutJrPgpAa+xMQFp1g1nY1BSdcZIXyzQA0VOnj0BIJbgrbRL+k7BWhaBZi2e7x5TQFJ99uJtrpcK3pew9JEqzxZkSe+Ty6bSI0rqLv71bUjP9+4qZDXMBIARGRKWAIqSJR9MkJZI8eEo8eZFyOU1WuI9MLLTJUkg4RlmQlFKKtlifxRRGYiMmUYEBqKyJpGhli2fkYKtAwjB8MSK1MwnxKUBMOyhiGgBIQSAsiUkiYKi0RBmFjyLeVfAgAuR5He5Z+szVZGmSr5u6CkHACAe1MVoJT+mxAygVL6IAAQQthsOwMlOueb6OvE/BVWSmlXtN4Dcuf4wmHusg/fOzb2fmFFVe0l4yf0K9dFNrtmY64dezggLvbvawKAco0nxGlkulPWGldJhuJ7xbyZVJQxWlvvyjeagm3GIv1+zmYeIXb6NvJlDktECLs0PLfGUFXkYzQal8ZoAOjuucGmrKZkA8UrhLizn3smLsOhzl2TykWUeKyJ6HXfEpN1MvW6Z3GC50gO7rl2S9unKHPUB0WT7Ovitrg6Sr8mtPhNvd7nM/MarZ6TzpRdte3Qamuojt/HGExGWYwYQOV66PhNTJ59hBzw7SYGUxHCwmy4OjcQe/5oBLwARTFCwmgaCGaj8Lt8Xm7chWfF5Zjr3DXJM4B0IZMxxc8xhG7uEi1H2XXuEh0r1gYjuqUl+rbk83mWUFi5YKHJEN5KzeamoKQLBSLQMYQucjDN+zq6DHs6PJa9zc02WQiQCCLcP7Zr7JLOSMrMRs/sQlt7kdkaXNEmlQsAy8qQgwL4IICgDN3LTYHKk4r1nYuK9e4+ZdjVxY0bN65bhlnkrklW5L0p9kTrvdquF354wrTGkSMqffOq7e4Sq05c/m1rwTlTS51VebwQt/JZloYkiexy+nkhIpEWT4jb2OC2zCm3eSrtvFBh48XP9nWa/725vqzVHeJnVOS5XH5BN6vC6l69z2Xf0eqzcCwjLxyT3/zR7vZimVIyKt/g2dUesGoIlZq6ggYNy1BXUNbpWFCHkQ2OLTR31rr8RgqNlmVo1oNkRkqeUroXwM2ZFCWEPA3Fkm/OtjMAGqAo+q/QM7zTHV3MpQC8ySfFctf0o70BkejaGRnu6uFzHMeE/SMY776TDKTJJRPtprCmYJVfqoq4a9lxbKS9XCsFZJsN3+nyzdt5hy1ENJrJQadrH2+PSGBIFzj2G32BhYBgerDFbUREHsz0C8OduyaViyjxWCkNCidHFX5M6c6inZ4IIeRjQ75Gb+R9F8nuvZsivH1LiBZqXPVwCyHja1rtKJe5hA+ZrIyV1YglCPr2MwazB1qdFRGxjuH5Ouitfkar2yvxTeOoz/sKU1zZSnhjOQ14psjerrNpU7MdkpSq38kMdu6a5BlA8mwgeRCoNmiEC8qNzmqDRqgJRHRvNQcLYrOAVK6esSZt6JxSQ+vbLcH804p4Z4VBK8zL13m2uK2mj9uC/paQ5NrmkfJ+Mcq4d4PTa3l9b3sF5AC21DY6PP6dXDAc0fJarRjmjNoQZ2JZnVEGb5QmmPSdk8xafybXGMtdk6iEcyO9niS7dk4Y4/AQQrC+rtN6trWk7bxpZfEZwKTibjmzlKLVJ3LLt7c4fEJEc3R1XudXLV7L+GKLs8BkkE4cq/HQaJhvsUUnvvp1c3G5zSSYeV/k6qOraopMvLjD6TfoeY00odDsCYoyO66I9UYkSmu7BOukEouzxSMa3UFRd+zoorYFYx2uhz7dO3JPe8C0cldH/hmTy/ocKBPJSMkTQm4CMCtafmPMWk+GUvprQogl+tqTqkwfvA7gEULIaQCWE0KepZReCuBvUHLaECjumoOCvlw7DVoL/7plfKWeRiKLvXualsDVyIvByEewlHwRIoV7WlrytTpjxCqHgiU85y2QAoIhJHV+wxVaNhjKHLs4ex4IIRR0V3XEGwyDIWsMlQUTxTZ3utDMwzUPT7p1ghKEhB/QxnotpbQZvK5RZzMt0XTVnWcUmtrAcw+JBWNcwTDHNW1n3LKsf0NvHktMnCwZzJhBva1+ouHbGJ2RAljL5lfUUKO7nXB6GSAijNadrMUhSQy5nNY1HAo5d5IHgdj7kESJKFMSs+J7i/6pMmqE60aZm8aatKFYlEypnhXPLDV0bHIJZiDCjrfogiPN2tBqNy3Y5bNbrSYI9hIEPDJ0BWw4sK/DZ5UFP9G4Wwht87O7pCrbNq/VWGHUZKygsl0ozZZk1w6vYen80fmeartB6C1ap8YV1K2v7bROKbb4tjm9pplleT5Ow3hFSSKhiESaPSK3vdVnjtVtN+gafUKYKTTpIieNLe5q8oS4DV805B83ssDlFyIsoRJdMLagbUNdl3mUw+hp94WNFXl6/+LxRa3+sMS+uKWp1G7QCEcUWzpPn1Tclu11ZuquYSilFwJAH2kN+qvcY+f60fNhI89Hj38LJX3xQUW6OPpEBaunkciMYHNHeXQzVT1n4SVDKflJoHZniEbIc5x9VMjXpt3WSYu/0+qKinmtt4KPeLbyZYVVQqdLAqtpY026nVyBuUDyiQZJCL9vGFHSyJr5k3z72pJ9970tHB9uNBG9bjXJL1hAFT9+CUICJFLjQEjkQGkBBPFYXcSp09uozW4INUa0hq8FlGh87Qi3NWg3sWx5Ac97HGaHX9YZ0cXwxk5odAEwOhmgHqLhLJDCtYze1CTrdUMZnZNrEq14niU0cVF3hzfMxyz65HKJ5xZwTOi1pmDZ+WWGBo4FFcOE2HVMyBGWdBYNI3YIMj/KwHovLTfVfWzVWba68xwyBXGHoZtewLkyteRj9CeePRtSuXJ6c+8k9uuC6WXOEgsn1rpCHhAKgOKtb1sLzp3KOqvtemHppOI2UVI2PykhmDJ/2ZyKprEFxlCxhRMvmF7edPI4h6szENHGonVmltl8a0uM/k/2dDiMnCYyp9LmXlfrstR3BQ3FVl3wR7NL6+ZU2bKSIZC5kp9ECLkoWr6EEHJqqkcAft9JVrDlYU9oqXdXY0zpi2BIGAw5OlDfZg8HxFatSfejcP3eIsYnNJgs/FYQ2+5gwLrJpykTXR62UW+05RtMwf2c2RhgeO5rfVGhXhbEZo3FXMPl2WSAFEgBoSjiF8rC3lCr1qQLgmVGiS6PORwKH2qhm9mS7Mc3QZYrEQytJI5CA6gTALqIRrdUbm0AKD5h8/IDRrOr0JjvK5H1/k8lS4UrEDDKrhZWEoIkwpuJ12jjYbRFIhpOA7BaDQ0LxXLQb6HCIf1EsmQ3TszCT1zIHW/WhlJF/cSOuUSJHdOp9ZXoWfH1pkBBiY4NjTdz/nNLta0jDJrghi7BcoJD3xWWKXmqPlDlkaEzEUQqjazfGab6bd5wVpZ8Jgp3OEjsF6cR6atfNxUunVTclrjAy2kYqsxCWCcAvLWtpWDppOK2GldQJ0oSafaEeL8oa6aUWoLukMisq+kyO4zacGdQ5n44q7KxzS9oP9jVVvBti9ds1LESC4Z85/SZmz2i28pzWcmkTyVPCOGhxMlroLhLPgDgIITwlNJQQrkHgfjjmL6XD/JOdt8kW/qtWpMutlPWpTVwnxsqCo4O1Ldt1xeYN/JlDoZQTNd1drSbRhk40Su6AgFjS0eLeQeYQmq005HaYPtYGurymHU6lsryan15cbvWbKoKd3XNDja17eIdNgDQy+HIbi7P0qQ1G85276gbG+4MDJdMck3y5qZk69pBQ+IE+NwOqljzJycMAiUQhDBtaQaAT5HvEDkjO42VmqrNVt9ulFh3isQR9rs1YmMLz1AK2WijPqOVW63PL5sm+zzH0Y4hXQzPJemiedIp/+Rzqw0aQZQp//PR5roqvUbQEkLfaw3ln1qs7zgun/escwlmAoYxaYn82B5feZ1fNh9foGk6p8zkLNKx4r5ARD/DyvmG6nqHiuTonVTHAeDcqaVOUZJI8oAAAOtqusx/W71vxHXHVNdcML3MKUoS2drsNZ82ochZlccHvmz0Wq28JnzKuMHb8Xo/lJDw7VAe5F0NYCSAvQDWJpR7JOm8w9J67I2+InOSB4HjA7XOMBiynSuwzgg1dRRFAoI9HBBtUkgUoGE2msoKOiz5Or/GzDC+TlLj73LsFkLFVGiV7XrezxjzJA+j0e7XWq0dDKcjFCiQfD4f4bSfG0tLOMh0uy6/szrsDh4u1ny63bgx5R8GyHfEZC2AICqfKJfNgdIx1B/0gWHeIUUFOxiTpYwGfUZGljvB6yfLvk6jzhDmOKPss5tYMULZ3aFwvtvTYZA0DljZQ9uST0cmoZxAT3ePjiU0TCnxS1RTotOItcGI7t3WoGNhId/RHJK4/YGw4exyfd0vR5tr64KS3qhh5J2+iLlFlLyFevagf6xfNrto0802ko/HNmrF3DmxekMRiTgMXPi6Y6pqjhtl91h5Tg5FJLJkQlH75sZOU6svzOk0jHzOlNLWk8YVurPd7QpkoOQppTcSQvIATANgB/AepfSbFOVqCSFHQnmwiCF6+IpsO3Q4kzwIFIV9QoPWwh8bqGuL+exFMKQoEhDCYMgZ7p31dRqznppHAjoiu4zlhiBYrSbQRd1BHx9w79VzLBvuMBfqOo15BsJy2MPk54PKlCVAUSTo38fZLa1ap/tw8c+n240bU/7H0Y62BbTdGQbIclJcCkJwttzcGBsQthKbeQVbXBkGw54mNdX6oWHricHsZnS6EjksfsNa7PlSyL+Hs9kt2nAwbCnVFlJv0CAffPlxhpJEi78mENF92i7YWQJwbE9j7kibznvrONv+uXadd3VHyPLn3d5RZ5XwjV1hWSNKOCQSmg3Wgm+qAaHGFdS9t9PpOHdqqdPKc3KsXHtA1L6+taV03sj8dl7L0tEOU6g/Ch7IPISyE4rLpi+uAtACJQLmoFsoPVhIlegMAOq1Fj4MhrxvGlUKAizx7m48Pljvqox4g+sNpfmfswZDAQ0FxsLfsdJRXU00OpEJusOcz6URm1q1BJTC5ADMBZA4XupiOJ1eEsP2cEAc7mvOFemibBKVf2yz0+nRJGeJA8IU2uW9TGL2xB4P+TrrqDpBam8uoIKwluTlG2lEDDEsW0iD3gmy29XCGAOL5NamKho8qAfJwc4Tk2jxVxs0wpXVpiZACbkEgMVF+vb1LtE6xqQNLSrWu0MSJWGJkiq9xlfMs2JtUI4kDwgHK+kWfLOx8Afa1txqmxcYtX9GucXnCkS0A1l8znThFYSQiwGURd82xlIOJNEKgIfyCyrsb6cOd1IlOks8dopvbxOgWPocZGpARJbAkmOD9c0drEHHIwIrRMEUCYfDWp1GyiuVmLxKKUAoB287gXMvg4jIBox27XZzQb5La2gxHiaWfDoSlb8IQmqJgQcoqmgwFAt9FEFIO+G5ebSjiwOlPjAML2H/eOr2dRCeO5W2tmynZvNGkufII1K4k+H5ibK3axbtch/s4ZN9JUQbKMmDSHIumvkFvKfaqBFi0TqiTMlWb9h8YYWhZYaV81m1rFSlT53HZzhJpbjTuWAGw8JP15aV5+RF4wvdAFBoQkZ7NNKRaT55ACijlD5AKX0AyoalVDwP4P8A3ALg42w6Qgj5MSHkCULIckLItKTP3iWE/KOv8M1DhVR56hOPjQx3BWObq+q1Ft4eDognBfY7z3LvbDw62NAmg8XZnp0117dt2jY32NhQJnZ2mRAKcZSJcJbCcFlheTtTPimk1VsiMwNN37sHmzQRvW45KSpdzpSUJua8ryUG/nVSXKYMAMrAUIqQ0EJ43XKmuLQDHLeDMVtNRAofJbvaJMpoNrD5RTuI1TR8V5MZfSVEGyipnjqVSGxxdp1LML9YHygEgAvKjc4ZVs63vCWYv7ZDtDULEjcYfRsI2eSW729I53A/pCRjSx7AOkLIzVBWstamKbOUUno/gJ9l2xFK6ZMAniSETAdwOnomNwtAGZBas633YCTVAm2qY6li3qcGWz3FEb9QFPYJX+uLLA2cxVoa8fvrmXy7DYJQJnjcNZwtz0HDoZDBxJnFcPhwWXTNlFIaFFK5agAKEIKYUV5L9PxypqR0mtzlUo4x0NNIZAwNek+gzo7J8HrbJY6bQrsO2GV9sJHpAmp/yWQQqQlEdOs7BetRdp07tpHqvZag9eXGQOkPygxNgzUADYRsFHd/QzoHe1NXXxDaj1koIeRkSukHKY6vgDIAuAGAUvpYmvOPAPCHpMNXAHABfcyxywAAHclJREFU+CeA31JK6xPKM5RSmRDyEIBnKKVbE0/M4/mASauN+50HkrvmYEIEQxq0Fh4AylM8zcoPDbNdX2AuC7oDm0yl+QYpHKYA+cxUWTQ62OH+zlicf4Hrm71Tw20+AFj4v1dLKKWzUrVlNxgCZh0fl+FQ564ZCpLDL3cTo345KSpdRJ3NBsiyg4bELSTPspHY8s+gLU1VNBBKdvvMeuQvaWUIAEX2/IDNZI7LMde5aw5mUq0LuEWZWecSzHPtOm8sORmZN6N3GRYVBWw2W7cMs8hdczAyGL58ACCE9CrHGNn45F8A8CWUGPh5UOLlk0lMd5D2YqLROUuS6tdCSTP810QFHy0fi2xwAjhg6jxcuWsGGw4y1UKmMWs+5ruP+eqNiMhTg62eBq2FP0JodwPAKkNVwZxgU7tJEsNe0aAzI5yRP2+4c9cMBcmLtlU0EDobLY2J+XBKEBQYYgMojbt9kiN0emOwc9cczKSaTVg5Rs4kKVkisdw1ue3d8DHcm7qycdc8Sin9DFB85GnKHEMp/X20zN0AVmdR/x8BjAVwLSFkJaX0lVjuGkLIv6G4bDQ4iHLXDAWJsfWp3DetWpMuFo1zindv0xGi0/05X+bo0PP8jKCz/fvmj8+GVJE6VTQYOltuboy5eVK7fVRUDh0yTVD2LgANURKr+QD8G8DWpDKvAJiYsGiaVdgepfSXKY5dGv3/exuOmeirT5UQrSjsE07x7W0KRx8aPjHY5m1ljbpd/Gj7Ll6yzRRa3Id7ZE0uSVb8Y6j/oN+8o6LSG5la8hug5JGnAH4LYAGANxILUErPI4QcQylNtyirMkDS7ajVRn31sZj7E3w17flSQMyPBEXVkldR+X6TaQjlOAAlUB7pNxZAuod5nxx7QQi5bWBdOzh4bsd3lX2XGr72Yi4cQEmTUBT2CUZE5LnBpq6x4c7AwRBZ8+SmDUMqw/62J4KQGmLgxWgu8IOJe5/555DKcKjbGyruvffeoZXjELeXikyV/N0Aboz+3QMlv3sq8hJe5w+gXwcNH9XXVh3M7aWKuR+svvWXd3Z+N6Qy7G97sdQIibH1BwvPvf/OkMpwqNsbKp577rmhleMQt5eKTNMa7ACwLIOiddFF0lhCM5VBJpvHFar0Trq8OCoqhzIZxckTQpYCuASK8n6BUvpmmnIMlNQHnQBESumQ5EwhhLQBqB2k6h0A2gep7qFur4pSWpDqA1WGGZNWhsCgylGVYW743sgxRqZK/glK6dXR1/9HKb02TbmHARgppVcSQh6nlF6Tba9VVFRUVHJHptE1ekJIbAHB2Es5Cd0jcFYbIFRUVFRUck+fC6+EkFMBrAfwWvRvYy/FBQATCCE3oOci7CFNb8nTBqm9OYSQlw6XhGyAKsNcoMpw4HwfZdinu4YQkrwRiVJK/5OmLAvgJEQfE0gpHVCKzIONaPK0JZTSe4egrWoAN1BKM1nwPmRQZThwVBkOnO+TDDN5MtS/s6jvIig5aTRQQiif62e/ho0+kqf9DMpmsKFo75BFleHAUWU4cFQZKmSTuyYTRlJKzwcAQsijOASVfLbJ0wajvWib1blsZyhRZThwVBkOHFWGCtk8NCQTLISQowkhxwEAIWRijusfLhKTp5032I0RQsYCuA/AQkLI1YPd3hChynDgqDIcON87GfYrn3zaygi5M/qSQvHLU0rpPTlrQEVFRUUlK3Kt5CdSSrcTQn4IYDeldF3OKldRUVFRyZpcu2t+SAiZAmAEgKtyXLeKioqKSpbkWslXAbgGwFNQUhuoqKioqAwjuXbXjAdgpZRuIITMp5SuylnlKioqKipZk1Mlr6KioqJycJHrOPmsIIRMhhJXygGAGomjoqKiklty7ZPPll8CWAHgpeifioqKSr8hhFxGCDlgQ1I/6vkbIWQSIeStxHj6aDr1YYEQMpMQknVm32G15AF8Syn9dpj7oKKicnhRTQh5DkAIwHJK6ZuEkAkA7gSwE/j/9s48Oqoq28PfrqoMZIIkBBKmMIPMQxQFAaFRsBX7aYPiAA3PGbVbxddo+/rpa22Xs9iC3YL9nBUVpQUbZxBEBiGMhnkIYwYSMs9Vtd8ftwIBQ8aCm4TzrVWrqu4959xf7ara995z9tmHIao67kyVRSQKa67PNCwfeUREpgJjgPUishCYgTUXaC/wJfA4sAu4CJgPZKjq5yIyX1UniUj8aXWygZHAPsCrqn/1lXkSSAcWAr8FnsfKR/+6qk4WkQeA12pjDLud/CgRuQzry0BVr7dXjsFgaAJMBSaq6n4R+Rj4DLgNmAkcAb7yzcqf7nPAI4BLgKOq+g4wAGtlu11Asqqu8s1c/UJV3xORZ4Ai36Mv0L5i22fQNP20Oj8AX6rqhyLyga/MPcBfVHU3gIikAndgnQjm+8oUi0hrVU2rqTFsdfKqOl5EevteJ9mpxWAwNBkE60qcCs/l2xRAVVeIyFDfvotV9RkRmel7H4V1pX065WtkOIB3VHULgIi8WLFtrJTr5b419Ax1pgIFFbSVP3vLD6aqe0SkDXABMMm3OQuIABqHkxeRl/EtjSUid6rq7+3UYzAYmgRvAU+ISCFQfpU8DytvzS4gv5r6u4BxWN0mlTEbeEpEUoC8Cm3vBcqA5cCzItIJaHGGOpUtb/gq8LivzCJVXQV8B8SparnzbwscrEb/KdgaQikiL6jqDN/rZ1X1j7aJMRgMTRZfP/v9WCnQv8NyyE8BrwCFwFDgiKq+IyIC/F1V76rDcRao6gQ/aR4G/AG4VVXzRCQMeL62uux28i9j5XYGaKmq99kmxmAwGHyIyHAgUVUL7dZSjm9cwKOqe2tVz+7JUDXtkxeRzsCjWDNq/XKmNBgMhqaObU5eRJ4FOpa/xUpLXG10TWW3Qy1bttSOHa2mcovdRATbHTTUMElMTMxQ1ZjK9hkb1oyqbAgn7VhY5iHA6SDAIWcqet5SUxuWerx4vRAcYPd0noZJdXYsx7Z/sqr+UUTaquoRABFpV9e2QkOtAWwVBzl9xhPYvR+L/vBrOkeHVlPz/EJEKhvsAU7aEOBYws1MuzCKR+659ZzoakxUZUM4aces7lcQnr2feyaO4Y47msp6G/6hpjYsbN2LotAYHht7gbFhJVRnx3Lsvly7FShPZfB7oE4DrzExMaxfv/7E+1XJx5n0TiI3DGjLH4Z3wuU0VwLVUdGGD/97OwO7RNusqHFSbsfHv9rJoHbNuaZ3rN2SGh3lNlyclEri4RzuGNvDbkmNGjun6H4M3CAiH4nIR5yMGT1T+WgR+QcwUEQeqars0I5R/HDvMPJL3Vz291VsOpJTVXHDaVzQKozt6Xl2y2jUBAc4KC7zVl/QcEZcDsHtNQkU64ud3TUTK3bX1KB8JlDj0KEgl5PHrujBxH5tmP7pFi6Oj2LmqC5EhgTWWfP5wgWtw5i3tlahuIbTCHY5KfF47JbRqHE5HcbJ+wG7+zEmi8haEVkqIlvOxgF6xYaz9K6hdI4OYdy8tVz/9noWJaVS6jZXWWeiV+twfk7JtVtGoybIZa7k60uAQyjzGBvWF7udfCxW/obRWLPUzgoOh3DHxfGs/cNwHh3TjRX7Mrnw5RXc++lWvt11zDj80wgLcuF0CNlFZXZLabQEuxwUm99VvXA5TXeNP7B74DUDcIjIn4F+5+KA/ds0p3+b5rg9XpbtzeSzpFQeWryNHjGh9G/bHIdAXombvBIPxwtKSc4qJCzQxdCOUVzVqxUD2zTHcR6ExQ3vHM0P+zIZbwYO60RwgJPMQnOSrA8uh+mu8Qd2Jyh7EkBEBgKzzuWxXU4Hl3eP4fLuMagqSal5bE/PRxXCg5yEB7uIahZIfGQz8krcLNubycsr9rE5JZfuLcMY1bUlo7tF071lWJN0+qO6RPPlzmPGydeRIJeDYrfpk68P1sCruRuqL7Y5eV90zekZ4mxJNSwi9ImLoE9cRKX7Q4Nc3DiwLTcObIuqsvNYPkt3Z/LfX+zgQFYRXlUEQVGCXE7iI5sxsnM0Y7rH0CU6BCsVRuNiWMcoHvtqp90yGi3Bpk++3pjoGv9ga3QNgIgEqGqjua8VEXq2Cqdnq3CmD+v4i/0lbg97MgpZuieDhxYnkXy8iMHtm/Orri35VbcYWocHnXvRdSA0yEWQy8Gx/BJiwhqH5oZEsMtJiemTrxcBTgcej3Hy9cXuVMNPAJ2AW0TkJVV9wE49/iDI5aR3bDi9Y8O579JOuD1eEg/n8O3uY9z0biL5pR6u7RvL1IT2xEYE16hNj1dJyS0mo6CUjIJS0vNLCAl00jkqlJ6twgh0nZ3x8/G9YvnXz6ncfnH8WWm/KWO6a+qPuZL3D3YPvIZj5W4GKw9zk8PldDAkPpIh8ZE8OqY7ecVuPtx8hIlvryc00MV1/eIY1SWari1DT3TrHM0p5utdx/h6Zzo70vNxiNAmIpiWYYFEhwTSOjyIglI3C7emsi0tj/YtmvHomG4MbteiGjW1Y0K/OO76ZItx8nUgJNBJQalx8vXBZUIo/YLdTl6BON/Cu+fFCF94sIvbhsRz25B4ko8X8unWFB5Zsp29mYU4RfCo0ioskF91i+Hh0d3oGxdebZ/++kPZ/M+XO4kIdvHEuJ50bemfnD3xUSEcyy+luMxDcIDTL22eL7QIDiDHhKDWC3Ml7x/sdvIvAb8GOgD32qzlnNMxKoQHR3bhwZFd6tVOQvsW/Pu2ISzdncHU+RsZ0KY5M0d3pX2LZvXWOKxTFCv3H2dM92qT3RkqEBkSYOYZ1BMTJ+8f7Hbyk1T1WZs1NBlGd2vJqK7D+GRLCre8t4Ho0EBuH9KBMd1jCKhjkrbf9I5l/qYjxsnXkoggF7klbrtlNGpMnLx/sNvJXyMiXfAtkGuW/6s/IsKE/m2Y0L8Nm47k8Ma6Qzy8ZDsXx0fWqb3LukQzY3ESBSVuQoPs/rk0HhwOweb1eBo9AU4TJ+8P7P7XTgNKbdbQZBnQtjkvt22Ox6us2JfJ3Dq04XAI1/WN4/2NR87KAGxRmYcgp6PBTShLSs2jqMxDQvv6DWZ7vIqzgX22xoLLIbhNCGW9sTt3zbWqeqD8YbOWJovTIYzq2rLO9e+7tBNzfkwmr9g/3Q9er7JkexrTP9nCJX9bSe/nlvFu4mHsXooS4Fh+CTe+m8jdn2xhxqIkHv9qJ+46Rnj0aBXG9rTzL2Xz/sxCHA8tJquwftdvZuDVP9jt5KeIyFci8rEvp/wZEZFQEXlLROaJyM3nSqABWjQL4J5hHXni213VF66GnKIyrn1zHQs2pzC0YxQbHxzBsruH+iaPbfOD2rrz3obDXDlvLWN7tGLFPcNYevdQiso8jJ27hrUHsmrdXkL75qw7lF3j8pWdRHek5/FzSi7vbThc6+PXlIonsR3peXUaMFbVEyfpfyWlENksgB+Ta2+zitSkT/5oTjHL9mTU6zhVoaoczSmu8gIku6iMY/klp2w7lF1Uadnk44Ws3J/pV43VYXfumj61KH4dsEBVF4vIh8B7Z0mWoRL+86IOjP/nWj7dmsJ1fePq1EZSah5TPtjAI6O7MaF/mxPbYyOC+ef1/Zn8/kZeXL633tFGlbFsTwa7juVz+5D4SruGfjqYxawV+1g+feiJsQenQ3jm6l6sPZDFvQu38tcre9bqmEM6RPL62oNc2imKral59IwJIyYssNIZxDvS8+j17PdsmTGSJ7/dzf+O7UFEsIuJbyWSlJaH0yH0iQ0nwOGgbfNgxJdILzY8mNdWH2B0t2iiQwIp8yiRIQEMeGE5Y7rHMOe6vpVqS80tJjYimOOFpYydu4aL4yN5eHRXLpz1A1f2bMVHUxJOKb98bwaT39/IwT9f/ou2Nh7J4cp5a3A5HAzrGElqXgn/d8MAnvpuN9EhAVzSMapWdiunPE7+821pXNopilK3lx+Tj7Mno4Al29O5dUgHdh7LZ+HWFFStuQkfTU5g1YHj/KZ3LGFBLvJL3CQfL6RPXAT3LdxK37gIvF4lq6iM3BI3F7ZvQe/W4fywPxNVaB0eRHp+Ca/+mMz43rGk5ZUwd80B7r20EzcPbMvD/97OsulDAeuO9HBOETe9twGnCN9PH8qWlFwEYcCLy0n6r8uIiwjm9o8288DIzuSXuFmVnMW+zAKGxkdx3VvrGNE5mq0pudw9tCNpeSU4HcKDi5JYctsQokMCKXF7iQkLRETYcDibQe1asO5gNg8s+rnGdrRtIW8AEbkduBwrXn6Zqv6jirKPAF+o6iYReV9Vbyrfl5CQoBWX/zNUjogkqmpCZftqYsOcojLGzVvDlIT23Hlx5c7yTCSl5nHzexv44JZBXNA6vNIyZR4v0+Zv4mB2Eb9LaI/LIbQMDeRX3VrWK05/1op9LEpKpW9cBEdzinlufC86RoUAUFDi5sHFSexKL2DuxH50iwmrtI2U3GKmzd/EV3deckYbwql29HqV3761jlXJWeSVuCl2e7lpYFtmX9eXOT/uZ2SXaMKDXKxOzmL6p1vp3yaCzUdziWwWQH6pm6iQQN6cNIAeMWFsScnlzgVbSMs79YpxyuB2vJ148irf6RA+nDyYb3YdY+2BLIbER+JVZWTnaN5ef5ihHaMQgaeX7mZKQnu+3nmMKQntyC4q48NNR/n98E589nMq913aiYn925BVWMacVck88c0uglwOHh7dlQCnle7imat6Uez2MPLVVSyYkkBKXjGP/Hs729LyyXxiHNM/2cI/Vh+g5JmrTkR3FZV5CAl01diGjocWA9bdZGV3GE6H4HIIl8RH8v3eU6+QJw9uR+LhbLal5fP4FT14/OvKczG1bR5MfombnAp3UnERQaTkllRaPjY8iNQK30Og00Gpx8sHtwzixnc3nOljncLTV13Aa6sPsP944SnbJw1ow/xNR89Y75/X9+eF5XvZlpaPvnBNlXYsx24n/4qq3ud7PUtV76+i7GQgS1U/F5H5qjqpfF98fLzGxJwM8bvjjjvMwr+VUJWTr6kNC0vd/OWb3azYl8mtF3VgdNeWdIoOOaWMqrJkezqfbEnhYHYRLoeQklfMh5MH07NV5Q6+Iul5JcxZlQxAVlEZK/dl0iU6lMfH9qB3bPX1K+q4/7MkMgtKeWPSAAKcDhZsPsqsH/bRvkUzerUOZ8n2NH7brw0PjOhc7QBpqdtLUICzyj9WZXacNGUa7284QkGphwVbjtInLoKSMg/vbji5KNpjV3RnbI9WDH1lJd7nx5NZUEpqXskpnzcpNY9OUdbch0PZxXSKCmHcvDW4HMKs3/ShzOvlo81HWbA5hYVTL6RLyxBueCeRyYPbMfHtRIZ0aMGgdi34+6pkerUOo09sBNf0bs1Ng9oB1knJ4RCO5hQz5rXVdI0O5Yf9mQS5HDx5ZU/+seoAezMLyCl2M/XC9nSJDqF3bDjrD+XwpO8uJz2vhGK3hw6R1m8iaObnlHmUT36XwE+HsvlgwxEO/PnyGtvQK05uu+1WljcbxMdbUnj3poGICKsPZDF75X48z12Nx6u4nA6+3pnOuHlr6dU6jG1p+QBMvbA9wztF8WNyFlMGt+PFFXt5YERnLmzfgrA/fUH2k+PILXYTFRJA2J++oFfrMGaM7ML2dOsO4bUJ/RnWKZJlezLp3yaCFs0CCH1kCSM6R7E7o4CMglJGdo6mZ+twZq/cz82D2vLJlhTenDSQSe8mnvK5Pp2awDvrDzO2ZyvuWrCFDQ+MYNBLK/i/GwZQUOrmcE4xzyzdw0OXdeG11Qd49bd9eXPdIb7bncG0i9qzOjmLHen5tIkI5rp+cbxybd8aOfkTfWl2PIC3gZHACODtasqGAm8Afwdurrhv8ODBaqgeYL2ewb61tWFyZoE+v2yPJry0XAe+8L32ePo7ffKbnfrENzv14pdX6LT5G/WrHWmaVViqx/KLtbjMXS/ty/dk6EWzVugbPx3UjPwSnbs6WR/7cof+5eudumJvhpaUeU4p7/F49U9LtumdH29Wr9d7yj6v16uJh7L0/Q2HdX9mQa10VGVDrYEdX19zQDs++Y2WuT3q9nh197F8LXWf1F5bOxWXubWgpKzaciv3ZWpOUamqqsqMRfrmTwerLL/xcLa+vGKvJh7K0iPZRaqqmldcpm6PVwtL3VpY6tauT32r/Z5fphsPZ5+xna5PfatzVyfrNf9cqyNmr9Rp8zfWyYY/p+Tqwq1HT7z3er26oYrjJry0/Bff++mk5haf8v6Jb3aq22PV2ZGWpwezCiuttyb5+IlyXq/3xHH2ZuTrgeMnf0+FpW71eLx6/7+2qsxYdMpvtNymFb+7NcnHVWYs0l3peaf8JsopdXtO+UzV2bH8YfeVfBRwo+/tB6p6vC7tmO6amlHf7pozoWr1cS7Znk6g08HwzlHE1TD5Wm3IL3Hzpy92sO5gFlf3ak3fuAhyisv4fm8mPx3I5s5L4rm2bywRQQH8bv5GuseE8dcre/o1hLEqG0LN7KiqtqafLihxExLorLeG/i98z5AOkcyd2L9W9fxhw8ZEcZmHYreXFs0C/NpudXYsx+44+atVdQ6AiNwIfGCzHkMdEBGiQgK5ZXC7s3qcsCAXf/uPX47VTx7cnvwSN3N+TOaGdxLJLCjlsSt6MGlg27Oqp67Yvb6Avya1rbznUoID7A7Qa/gEBzhtzf1kt5OvOPTfD+PkDXUkLMjFzNFdmTm6q91SzhvCg+12H4aaYPe3FCEit2JF19QtzspgMBgMZ8Tue63pwFEgBbjHZi0Gg8HQ5LDbyf8aK3/N74Cr/NHg3Ll1ydDifxq7joagvyFoAGNDf9CYbQgNQ0ddNdjt5Mer6vVqxbyP80eDDeHLgMavoyHobwgawNjQHzRmG0LD0FFXDXb3yTcTkQ6+1/5ZzshgMBgMJ7A7Tr4LcJfv7WuquqeO7RwDyrNYtgTOXsaimtMQdcSraqWrf5xmw9Pr2UVD0AA1tCE0yN9iQ9AAjduG0DB0nK6hSjuWY7eTn66qr/pe36+qs2wTYzAYDE0Q25y8iDwHDAV+BAToqqrX2iLGYDAYmih2Ovl44FLgB8ALpKmqWfnYYDAY/Ijd3TVvYE2EEkBV9T/r2V4o8CrWkoLfq+pZzzkvIp2BR4HmqjpBRG4CRgFBwN2+YmdVk4j8B1YIaitgDtZM4k5AANaYRxzwHOAB3lDVZVW0dV7a0KfDL3Y0NmycNvQd13Y7+vP/DNiehTLI94gCZvqhvclYYZkAH57jz7LA9/yx7/lqn55zpgmIxMrU+Z7v/b3AcODPvh+KA3jf2PDs2tHYsHHbsKHY0R//Z1W1PU7e6Xu4gTbVlK0J7YBDvtceP7RXF8pvjQ5g6TmXmv4beB04VpkGVa3JYqXnuw2h/nY0NmwaNoTG/3+23cnPAWYDC4GL/NDeYSwDgP2frQOWnrOuSSyeAb4A1mGFWv1Cg4jU5PjnpQ3Br3Y0NmxaNoTG+3+2deA1CJgIjAe6A1NUdWs92wzFOmkUAyv13PSFRgN/xVrG8HWss+1woBkn8/GcVU0i8nus1BDrgE1ACBDPyX7EOOBprDumd1V1aRVtnZc29Onwix2NDRunDX3Htd2O/vw/g71OPhnrduQDYLaq3l11DYPBYDDUFjtvg64FErAGFrr4ruwNBoPB4EdsDaGEE902E4BrVXWCrWIMBoOhiWG7kzcYDAbD2aMhjFobDAaD4Sxhd6phg8FgaLCISDPgJax5PJHAdmCnqr5Qz3anAhmq+vlp2y8D+qjq7Pq0XxHj5A0Gg+EMqGoRcFe58wU+B+4VkY7Am8BPWOGVKcAQ4GGgEJiBla5lr1bIrisiT/u29wBeF5GBwFQsX7waKwYeEWmDFbIZDXzp2z5KVZ8TkZeAF1W1fFJWlZjuGoPBYKgbu1T1j1iTo2YDz2CtcDcdKAIysdIPACAizYFYVZ2J5dABHgSysGa1DqzQthsrLj4NuFlV1wP9RKQFEFFTBw/mSt5gMBjqSq7vuURVc0WkFMsxO4B3VHVLJXVKy+v4ngOBl1U1C05014CVI2cRsBb4zLftU2A+8FRtRBonbzAYDP5lNvCUiKQAear6vwCqmiMiKSIyAxgG7MG6+n9FRNKAZKB81v8qrIyTwzh5YlgMPKSqK2ojxoRQGgwGQwNHRIKBvwFLVPVftaprnLzBYDA0XczAq8FgMDRhjJM3GAyGJoxx8gaDwdCEMU7eYDAYmjDGyRsMBkMTxjh5g8FgaML8Pywn4Mf3a0GWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 396.85x180 with 12 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tss = pd.read_csv('results/without_interactions/timeseries_autocorrelation1.csv')\n",
    "\n",
    "fig = plt.figure(figsize=(ELIFE.TEXTWIDTH,2.5))\n",
    "gs = gridspec.GridSpec(3,N, hspace=0.7, wspace=0.1, top=0.95, \n",
    "                       bottom=0.15, left=0.15, right=0.95)\n",
    "\n",
    "autocorrelation = {}\n",
    "\n",
    "for species in ['species_%d' % i for i in range(1,N+1)]:\n",
    "    R = np.zeros(len(tss))\n",
    "    \n",
    "    ts = tss[species].values.flatten()\n",
    "    mean = np.mean(ts)\n",
    "    std = np.std(ts)\n",
    "    \n",
    "    for i in range(len(ts)):\n",
    "        R[i] = np.mean((ts[0:len(ts)-i] - mean)*(ts[i:]-mean))/std**2\n",
    "    \n",
    "    autocorrelation[species] = R\n",
    "    \n",
    "for i in range(N):\n",
    "    ts = tss['species_%d' % (i+1)][:500]\n",
    "    ts = ts - np.mean(ts)\n",
    "    ts = ts/np.std(ts)\n",
    "    ts += 1\n",
    "    \n",
    "    ax_ts = fig.add_subplot(gs[i], sharey=ax_ts if i > 0 else None)\n",
    "    \n",
    "    ax_ts.plot(range(len(ts)), ts)\n",
    "    if i == N-1:\n",
    "        ax_ts.set_xlabel('Time', ha='right', x=1)\n",
    "    ax_ts.set_ylabel('Abundance')\n",
    "    \n",
    "    ax_psd = fig.add_subplot(gs[N+i], sharey=ax_psd if i > 0 else None)\n",
    "    \n",
    "    slope = example_noise_fit(ax_psd, ts)\n",
    "    \n",
    "    if slope > 0:\n",
    "        slope = 0\n",
    "        \n",
    "    c = noise_cmap_ww((slope - noise_lim[0])/(noise_lim[1]-noise_lim[0]))\n",
    "    #ax_psd.set_facecolor(c) #, alpha=0.4)\n",
    "    \n",
    "    if i < N-1:\n",
    "        ax_psd.set_xlabel('')\n",
    "    elif i == N-1:\n",
    "        ax_psd.set_xlabel('log$_{10}$(frequency)', x=1, ha='right')\n",
    "        \n",
    "    ax_psd.patch.set_facecolor(c)\n",
    "    ax_psd.patch.set_alpha(0.7)\n",
    "\n",
    "    ax_corr = fig.add_subplot(gs[2*N+i], sharey=ax_corr if i > 0 else None)\n",
    "\n",
    "    ax_corr.plot(np.arange(len(ts))[:500], autocorrelation['species_%d' % (i+1)][:500])\n",
    "    ax_corr.set_xlim([-10,210])\n",
    "    \n",
    "    if i == N-1:\n",
    "        ax_corr.set_xlabel('Time delay', ha='right', x=1)\n",
    "    if i == 0:\n",
    "        ax_corr.set_ylabel('Autocorrelation')\n",
    "        ax_psd.set_ylabel('log$_{10}$(power \\n spectral density)')\n",
    "    \n",
    "    if (i > 0):\n",
    "        for ax in [ax_ts, ax_psd, ax_corr]:\n",
    "            ax.set_ylabel('')\n",
    "            ax.tick_params(axis='both', left=True, labelleft=False)\n",
    "    \n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}