Tutorial1_fittingAutocorrelations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 1: How to use abcTau package to fit autocorrelations or PSDs\n",
    "\n",
    "Details of the method are explained in:   \n",
    "Zeraati, R., Engel, T. A., & Levina, A. (2022). A flexible Bayesian framework for unbiased estimation of timescales. Nature Computational Science, 2(3), 193-204. https://www.nature.com/articles/s43588-022-00214-3\n",
    "\n",
    "To start you need to have:\n",
    "Python >= 3.7.1,\n",
    "Numpy >= 1.15.4 ,\n",
    "Scipy >= 1.1.0. \n",
    "\n",
    "and for visualizations:\n",
    "Matplotlib >= 3.0.2,\n",
    "Seaborn >= 0.9.0 \n",
    "\n",
    "You can install seaborn from https://seaborn.pydata.org/installing.html "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns # comment this line if you don't want to use seaborn for plots\n",
    "\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# add the path to the abcTau package\n",
    "import sys\n",
    "sys.path.append('./abcTau')\n",
    "import abcTau\n",
    "\n",
    "\n",
    "# graphical properties for figures\n",
    "sns.set_context('talk',font_scale= 1.5)\n",
    "plt.rcParams[\"axes.edgecolor\"] = \"k\"\n",
    "plt.rcParams[\"axes.facecolor\"] = \"w\"\n",
    "plt.rcParams[\"axes.linewidth\"] = \"0.8\" \n",
    "plt.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: Extracting required statistics from real data\n",
    "\n",
    "First load your data and use \"extract_stats\" functions to compute statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#---------- load real data  (OU process with one timescale)\n",
    "dataload_path = 'example_data/'\n",
    "filename = 'OU_tau100_T1000_fs1000_trials500'\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "#---------- extract statistics from real data\n",
    "# select summary statistics metric\n",
    "summStat_metric = 'comp_ac_fft'\n",
    "ifNorm = True # if normalize the autocorrelation or PSD\n",
    "deltaT = 1 # temporal resolution of data.\n",
    "binSize = 1 #  bin-size for binning data and computing the autocorrelation.\n",
    "disp = None # put the dispersion parameter if computed with grid-search\n",
    "maxTimeLag = 500 # only used when using autocorrelation for summary statistics\n",
    "lm = round(maxTimeLag/binSize) # maximum bin for autocorrelation computation\n",
    "data_ac, data_mean, data_var, T, numTrials = abcTau.preprocessing.extract_stats(data_load, deltaT, binSize,\\\n",
    "                                                                                  summStat_metric, ifNorm, maxTimeLag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2 (optional): Check the bias in timescales estimated from exponential fits (preprocessing module using parametric bootstrapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45% bootstrap-error\n",
      "The true error from the ground truth can be larger\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3cAAAFgCAYAAAAYf4AcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABfuUlEQVR4nO3dd7RcVfnG8e+TBEjCDQkl9BKQ3kx+IE2QUKSIVAULCAEsgAjYQAQBAREVFRREFCGggIjSREFASIxUUSIiVSAUkRIg5aZASN7fH3tP7mQyM3fO3Jqb57PWrDPnnL332TNn7l3zzm6KCMzMzMzMzGzR1q+nK2BmZmZmZmYd5+DOzMzMzMysD3BwZ2ZmZmZm1gc4uDMzMzMzM+sDHNyZmZmZmZn1AQ7uzMzMzMzM+oABPV2BxcEee+wRt912W09Xw8y62KRZs1j7gQcAeG7rrRkxaFAP18jMuoh6ugJmZtW45a4bTJ48uaerYGZmZmZmfZyDOzMzMzMzsz7AwZ2ZmZmZmVkf4ODOzMzMzMysD3BwZ2ZmZmZm1gc4uDMzMzMzM+sDHNyZmZmZmZn1AQ7uzMzMzMzM+gAHd2ZmZmZmZn2AgzszMzMzM7M+YEBPV8DMzMwWTU+f8DStE1urnmsZ2cJ656/XzTUyM1u8FQruJH0EuCki3u2i+lgvNU7jqh4fHaO7tR5mZtZ7tE5spXViKy0jWxY6bmZm3a9oy911wGuSfgn8IiKe6II6mZmZ2SKiZWQLo8aNWuDYw6Mf7qHamJkt3poZc7ci8CXg35L+KukwSYM6uV5mZmZmZmZWQNHg7v3AZcAMQMB2ef8VST+V9L5Orp+ZmZmZmZk1oFBwFxH3RcSngZWBTwP3kIK8IcBngPsl/VPSsZKGdWZFJQ2VdLakf0uaIektSfdKOkpSp836KWntXH7kx+jOKtvMzMzMzKyrNBUURcTMiLgsInYANgTOA14jBXqbARcAL0v6laSdOlpJSesA/wJOATYGAlga2Ba4GLhD0sCOXie7BBjcSWWZmZmZmZl1iw63eEXEUxFxIrA6sB9wCzAXGAh8ArhT0n8knSxplaLlSxoA/B5YA3gJGB0RLaTg7rPA28DOpICyQySNAT4I3NfRsszMzMzMzLpTp3VnjIi5EXFzROxDCsSuIbXkCVgbOBt4XtJvJW1doOgjaGutOyAixufrzYmInwOn5nSflrRhs/WXtBLwfWAKacIYMzMzMzOzRUanBXcAkoZI+ixwM/BxUkAGqSXvbdLSC/sD90q6RFL/Boo9NG/vioi/VTn/E2A66bUc3IHq/xhYDvga8EoHyjEzMzMzM+t2nRLcSdpB0ljgf6QxcO8jtdg9B3yd1JK3CnAc8HQ+92ngy+2UO5g0rg7gtmppImImMCHv7tpk/fcBDgTuBX7WTBlmZmZmZmY9qengTtLKkk6S9CQwDvgUaSKSOaTFzj8YEetGxLkR8WpETI2IC4FNgJ+TArzD27nMxmV1fLROutK5jZt4HcuQWv/mAJ+LiGgni5mZmZmZWa8zoEjivOTAh4EjgT2B/qQgDeAp4FJgbERMrlVGRMyVdCpp6YS127nkqmXPX66TrnRuGUktEdHaTrnlvgusBnw7IuoFkFVJmtJemi222KJosWZmZmZmZoUUCu5Is1WulJ+LNI7ueuBnpYlOGhERr0sCWKKdpC1lz2fWSVd+bgjQUHAn6QOkGTefAc5qJI+ZmZmZmVlvVDS4WzlvHyd1rbwiIt5q8tpn0jbhSrfL6+KVuoceHRGzmiknIoa1l2bLLbd0V08zMzMzM+tSRYO7X5Ja6e7p6IUj4owGkpW3wNVbWLz83PQGq3A6sD5wVUTc0WAeMzMzMzOzXqlQcBcRh3VVRWooH2e3KvBIjXSlsXnTGhlvJ2ld4CvANOA0SS0VScqDxUH5/Lw8M6eZmZmZmVmvU3RClbtIXSk/FBFvN5C+H3AnEBGxSxP1ewyYR5oxcxNqLIeQz5XSN2J10mtfhjTerp4/5u3zwIgGyzczMzMzM+tWRZdCGJ0fjSw+Dmk8WylPYbml7L68u0fVC0iDgB3y7p3NXMfMzMzMzGxR1ymLmHexK/N2F0nV1hQ4mtQCNw+4qpECI2JcRKjWgwWXaNgpHx/RkRdhZmZmZmbWlbo6uBuSt7M7UMZlpO6WAq7PyxcgaYCkI4FzcrpLI+KJ8oySxkoKSZM6cH0zMzMzM7Ner+hsmSWNTu1/YN7+t8nrEBHvStoHuBtYExgvaQZpjbwlc7K7gOObvYaZmZmZmdmirm5wlydQqeY2SXPrZO1PmsFyHVIgeHdz1Usi4hlJmwEnAvuTJjaZBfyD1G3zkoiY15FrmJmZmZmZLcraa7kbTQrOVHZMtE1g0ojXaOs62bSImAqckh+N5hkDjGniWpNY8DWbmZmZmZn1au0Fd39hwS6YO+b9e4B6LXdzgMnAg8AvI+KNjlTSzMzMzMzM6qsb3EXE6PJ9SaWuj3t4QW8zMzPrq8444wy++c1vsuOOOzJu3Lierk6PkLQj8DXgfcCypIn4roiIMZLOAE4Hxld+XzSznlN0QpUzSS1373RBXczMzMyaNmbMGK644ooFjvXr148hQ4YwdOhQ1l13XUaNGsWHP/xhdtxxR6S+NwIjB10AY/Mwk2bL2Q74M2kehQDeIPXamtpA3v2AkcDEiLix2TqYWXGFgruIOKOL6mFmZmbWKfr168fw4cPn77e2tvLCCy/wwgsvcNddd/H973+f9ddfn4svvpidd965ahkrrLACG2ywAWuuuWZ3VbuznJ6344BJHSjneFJgdw+wT0S8WXF+MvAk8EKVvPsBhwFXADd2oA5mVtCisIi5mZmZWcPWWGMNXnnllfmP1tZWZs+ezX333cdJJ53EsGHDeOqpp9h111256KKLqpZx7LHH8sQTT3DllVd2c+17jc3y9tdVAjsi4sKI2DAiDu3meplZHTVb7kqLhQNExF8qjxVVKsPMzMysuy211FJss802bLPNNhx77LHsvffeTJw4keOOO45NN92UHXfcsaer2NsMytvWHq2FmRVSr+VuHGl9uruqHCv6qLVenpmZmVm3Wn311bnxxhtpaWlh3rx5nHLKwqssnXHGGUhi9OjRC52TNFZS5K0kfVbSvZLezMf3q0j/PklXSHpO0mxJ0yQ9JOlkSUPq1VXScEln5vRv5fzPSLpZ0mGSBpbXqSzr3bkupcekRt6bUnrSmsIAl5eXU5bujHxsXNmx0TnNYfnQYRV1CEkLv6Fm1mna65apKmnU5MPMzMysV1hrrbUYM2YMAPfccw/PPPNMM8UIuA64BNg6H4sFEkjfIi0NdSgpYHqX1Cq2BWkd4L9LWrtq4dKewNPAN3L6pYHZwDrA3sBYYMOcfCrwaln2t/J+6fF6g6+plL40Q/q0inLqeSenmZ33Z1fkfRVPymfWpeoFd2uXPaodK/JYp7MrbmZmZtYRe+211/zn48ePb6aIA4B9gC8CwyJiOdKSAfcCSPoy8HVgCvAVYMWIaAEGAx8A/gGsB1wvaYHvZJK2AG4AhgIPA7sCgyNiGLAMsAPwM3KwFBHHR8TK5XWLiJXLHu9r5AWV0gMv5kPHl5fTTt57c5pr86FrK+qwckTc20g9zKw5NcfcRcTzjRwzMzMzq+bdefN46e23O6281ZdaigH9Om8uuM0333z+82effbaZIlqAz0fET0oHImIqpO6UwFmk5QP2iYgJZWnmABMk7Qb8i7RswD4sOLPkBcBSwKPADhExoyz/dOCv+WFmNl/Rde7MzMzMGvLS22+z9gMPdFp5z229NSMGDWo/YYOWW265+c/ffHOhCSEb8Sbw8xrnDiZ1v5xQHtiVi4g3JN0KHAHsTg7uJK0PvD8nO6k8sDMzq8fBnZmZmVlz/pZb4aopBWdbSnqlThktebtW2bFt8/Zt4M4O1M/MFjOFgrs8I9P6QETEv6qc35e0eOZGwBvAVcDpETG7Mq2ZmZn1basvtRTPbb11+wkLlNeZylvrll9++WaKqDdJyap5O4i2ZQXqGVz2vDS27ZWI8AQkZtawoi13B5MG7/4DWGBgrqTdgetLu6R/al8BNgX2wszMzBYrA/r169RulJ3tkUcemf/8Pe95TzNFzK1zrn/e/jgijmumcCpm3jQza0/RUcm75e1VVc6dSwrqpgNXkgb5Ctijcr0XMzMzs572hz/8Yf7zLljEvNQVc9MO5F1V0pKdVB8zWwwUDe5K/6AeLD8oaWPgvaRfmHaJiDER8QHS2i8CDuloRc3MzMw6ywsvvMDYsWMB2GGHHVh77apLzXXEPXn7fkmrFMxbWi5gSeCDBfOWWvt6co3h0hp5XufYrJsVDe5WyNuXKo7vnLf3RcTfy47/OG+3KFoxMzMzs67w0ksvse+++9La2kr//v0555xzuuIyvyIt4r0kcIGkmoGOpCUklSZWISKepm2Zg3MlLV3gutPydlix6naq3lAHs8VS0eBu2bydXnF8B9IvRbdXHP9P3q5U8DpmZmZmnebtt9/mwQcf5OSTT2bzzTdn4sSJSOLCCy9k++237/TrRcT/SAuYAxwI/FHSVqUgT1J/SZtLOpX0fWlkRRHHk2bL3BT4i6SdJfXPeYdI2lHSr3LvqXKP5u3BkgbTM0p12EHShj1UB7PFUtEJVWYAy5BmcXqr7PgH8vbeivSllUvrDTg2MzMz6zQvvvgiK6+88vz9mTNn0traSkTb/CQbbLABF198MTvttFOX1SMifphnGj8b2CM/3pbUCgxlwe9hUZH3H5IOAH4N/B/wZ+AdSTNZsEXsvIrL/pS0DMNHgH0kvQa8C7wUEZ0fxVb3O+AcYDjwuKTJpO+QAB+PiPu7qR5mi52iwd3TpC6WewKPA0j6AKll7l2g8o+11GL3agfqaGZmZtawefPm8eqr6atHv379aGlpYfXVV2fddddliy22YK+99mL06NHdUpeI+LakG4DPk4axrEkK7KYAT5HWsbshIiZWyfvHvKD5CaTvXusASwHPklrHriN/HyvL86vcOPg5YDNgFYr31OqQiHgrfz88ndS7a0XahvYM7M66mC1uVP4rVruJpTOA00i/vpwMvAx8G1gXuCMi9qhI/0lSn/PxEdF1P431cltuuWU89NBDPV2NDhmncVWPj47R3VoPs95s0qxZrP3AAwA8t/XWvXoKeLPO8PDohwEYNW5UQ8f7EE8UYma9UtGWuwtIvwStlJ9D+gc3j9TloNKHSd0M/lrlnJmZmZmZmXWSQs30EfEWsBNpel/lx/PARyNigQBO0vLAfnn3zx2uqZmZmZmZmdVUtOWOiHiCNPvREKB/REypk3zPvL2nThozMzMzMzProMLBXUlEVC6HUHn+DWB8s+WbmZmZmZlZ47p19iQzMzMzMzPrGg7uzMzMzMzM+oCmumVK2hk4EtiKtKD54HayREQ03QXUzMzMzMzM6isccEn6CWk5BPA6L2ZmZmZmZr1CoeBO0uHAUXn3VeA64HFgZifXy8zMzMzMzAoo2nL3mby9B9gzIlo7uT5mZma2iJv1n1nMbZ3Lw6MfXuhcy8gW1jt/vR6olZlZ31c0uNsECOBEB3YGME7jap4bHaO7rR5mZtZ7zG2dW/V460R/dTAz60pFg7vS7JqPd3ZFzMzMrO9oGdnCqHGjFjhWrSXPzMw6T9GlEJ7N22U7uyJmZmZmZmbWvKLB3bV5+6HOroiZmZmZmZk1r2hw9wNSl8wzJG3QBfUxMzMzMzOzJhQdc7cicARwFfA3ST8CbgNeBt6tlzEiXmiqhmZmZmZmZtauosHdJNJsmZAWMD85P9oTTVzLzMzMzBYzksYBOwLfjIgzuuF6Y4DLgecjYkRXX8+sKzUTcKnGczMzM7OFVC6bU28ZnVoaWV5nzJgxXHHFFQsdHzx4MKussgrbbbcdRx11FNttt13h61ttkk4AhgE3RsTEHq3MYkTSSGA/YEpEnN/JZe8HjAQmRsSNnVm2da2iwd3hXVILMzMzs07Sr18/hg8fPn//jTfe4JlnnuGZZ57hV7/6FaeeeipnnnlmD9awzzkBWIvUw2tiJ5T3AvAkMLkTyurLRgKnA88D53dy2fsBhwFXADd2ctnWhQoFdxGx8M9hZmZmZr3IGmuswaRJk+bvz5kzhwkTJnD00Ufz1FNPcdZZZ/G+972Pvffeu+cqaTVFxKE9XQezRVXR2TLNzMzMFilLLLEEO++8MzfddBNLLrkkABdddFEP18rMrPM5uLMuM07jqj7MzMx6woYbbsiWW24JwEMPPTT/+AsvvMBFF13EXnvtxfrrr8/SSy9NS0sLG2+8MSeccAIvvFB7wm9J4ySFpDMkLSnpq5L+IWlqPj4yp1tC0gcl/UjSQ5L+J+kdSa9J+pOkT0iqOpeBpNG5rMj775V0raRXJM2U9Iikz1bk2U/SXZLekDRD0j2Sdq/3/ig5UNLvc9nvSJos6Q5Jh0jqV5H+jFyntfKhy0v1LK9vTjui7PgISetKulTSJElvS5pY7T2tU9dNJf1U0pP59U2X9Jiky9t7ne2RtFu+J5Pz+/tPSV+SVLfHm6StJF0t6cX8mt7Mr+Uzkvq3k3cDST+X9Kyk2ZKmSHpA0pclDaySPkiTwACsVfm+V753kvaUdJOkl/N9nSrpP5JukXRs6RqlzxqpSybAYVXKHl1W7qR8bIykIZLOlvRvSa35+LCcbpCkffJrnCjp9fwevSzpRkl71nlvxuSyJnXk/iwumn4TJC0JfBDYChgOLBURR5ad7wcMBOZFxOyOVlTSUOCrwP7ACOAd0pp7VwI/i4h5BctbgjSGcEtgFLAKsAJpSYcXgPHARRHxaEfrbmZmZr3DaqutBsC0adPmHzv00EMZP378/P2hQ4cyffp0Hn/8cR5//HHGjh3LLbfcwvbbb1+v6IHAOGBb0neJ1orz7wduL9t/Oz+GA7vlx/6SPl7vO42kDwHXA0sBU4FBwGbAJZLWiYivSTqHNJv5XGAGMBjYDvijpAMi4qYq5Q4BfgPsUXZ4GrAcsGt+fFzSRyLi7Xy+FXg1v4Z+Of2sWnUvsx1wCdCS61d3Oa0qdT0FOIu2if1m5ecb5cf+pAleCpP0eeDHubwpwJLA5sD3gb0k7VXte62krwLfKavTFGAIadbPHYFDJO0TEVOr5P0kMBZYIh+aRrqvW+XHEZJ2i4j/lmV7NadZBpgHvF5R7PzPXw70Ti87NxPoD7wnP/YCbiGNmXwnlz2U9JmeTfqclXun8jUAywN/B9bL5ys/Bx+jLRgln3+X9P17X2BfSd+PiK9UKXu+Zu/P4qSplrv8IXwOuBk4FTgKGFORbG3Sh3NKDsyaJmkd4F/AKcDGpKUVlib9A70YuKParxrtGEr6x/IZUoC3KunDPoj0j+EoYKKk4ztSdzMzM+s9nn/+eQCWW265+cc23XRTzj33XB577DFmzpzJlClTePvtt7n//vvZY489mDp1Kh/72MeYNatu3PJ5YBPgUGBIRCwLrAw8m8/PAq4mfZFeGRgUEUNIX4qPJ31nOhA4tp2XcBXwW2D1iBgGLAv8Ip/7qqSvAScCJwHLRsRQYF3gQdL3vgtrtCJdRQrsnsz1GJLztgCHAK/lun+3lCEizouIlYEX86HjI2Ll8keN13AJ8AiwRUS0RMTS+ZrtknQscDbpy/11wGYRMTgiWoCVSIHCrY2UVcVw4IfAtaT3d1nS98UvkwLlnUkBXGWd9iW9LyIFyGvmvEOAo0lB/AeAS6vk3ZK2wO5uYKP8vi9NCojeIn33/W35fcvvbek76ouV73tEnJfLXwv4Rk73PWDliFg6v1+lwH0sOWCLiHtz2dfmPNdWKfveKu/dGbnO+wJL58/m2qTgHVIg9jNgJ2CFfM+WJn3/PgOYA3xZ0j5Vyi5p6v4sbgoHd5K+AvySFGkLeKNauoh4Bvgz6cP64WYrmJtYfw+sAbwEjM4fyKWBz5L+YHYGLihY9NvAj4CDSN0JloyI5Ui/AGwF/In0q8YPJdX9qc7MzMx6v4ceemh+d8xtttlm/vELL7yQk046iY022ohBgwYBMGDAALbeemtuueUWNt98c15++WV+97vf1Su+BfhERPyy1HIQEa9GxLT8/IGIODgi/piPRz7+ZkT8CCj1fjqunZfx74g4pNSKExFTSN+HniF9r/s2cG5EfDcipuc0zwAfJ/04vjqp5Ww+SXsDe5OCtB0j4rcR0ZrzzoyIq4ADcv6jJa3YTh3bMxnYLSL+UToQEU+3l0nSssA5efeyiDiovIdVRLwWETdHxCearNdgYAJwcNn7OyMifkBqzAA4RtJqFfnOzds/kz4DL+a8syPip8Ax+fxHJb2vIu/ZpO/KjwEfiognct53I+I3wEdzum1I96CorUmfiycj4sSIeLV0IiLeiog/R8ThEfFyE2WXGwjsmd//d3P5kyJiTn5+Y0R8LiLGRcT82CEi/hcR3wS+ng/V+/w3e38WK4WCO0mb0xYRX0vqHrlOnSw3kALAjkxHdQRtrXUHRMR4gIiYExE/p+1mflrSho0WGhHTI+L4iLguIl4o+yDOjYi/kX55eC7Xf0wH6m9mZmY96JVXXuGaa65h3333Zd681OPx+OMb65jTv39/9tgj9VT861//Wi/poxHxxw5U8w95+x5Jq9RJd3blgdyN8668O4cqrRcR8Rzwn7y7WcXpUmA5tvzLf0X+e0jfi5Ygtb50xIURMaP9ZAs5kNQaNpvUOtkVzqnRLfYCYDppSNNHSgclvRcoff88s0bey0ldHgHmB545WN0t736nWnfCiLiLNFRogbwFlLpULiNp6SbyN+q2iHikA/lLn/9t2xmfWOj+LI6Kjrk7nhTs/LH0q0g7H5S/5e17m6hbSWk63Lty0FXpJ8BppD/2g2lreu6QiCgN7l0bWKx/ATAzM1uUPP/886j63CQMGDCA8847j512WjA+mTBhAr/4xS+4//77eemll5gxY+HY46WXXqp32Xvaq1ce13YUqUfTRqRxYUtUSboa8L8axdT6Al0KyiaVWuxqpFmP1JWz3Pvz9jhVTMxSodSXda06aRrR7ntVw7Z5e395608neheoGsFHxCxJ95Pmm9ii7NSWeTubGq8rIkLSn0lBdHne/6NtjN6ddep1J2nc3hZ10tTyAKmX3SrA/ZJ+CtwREU81UVY9jXz+VyK1Yu4GrE/qUlkZyA0mfT6rrXHYzP1Z7BQN7kaTWtAa7c9a6oPdVHAkaTBtf8i3VUsTETMlTQA+ROo33CnBnaRBpD86aOsvb2ZmZr1c+SLmkhg4cCCrrLIK2267LUceeSQbb7zxAulPOukkvvvd+UPJ6N+/P8suu+z8ZRNaW1uZMWNG1YCvTOWEFguQtD6p297qZYdnksYilVoiVsrbmj+c1+k+Nzdv63WvK6WZH1DmCeZWyLtD86M9gxtIU0/d96qO0hi+5zt4/VomR9tkMdWUJjQp75Y6PG9fiYi51Fb6ZaBa3nnUDuZr5W1IREzJc2VcDWwKXAgg6U3SBEBXAzfUaA0ror3P/7bAH1lwoptW0t9AkIK80udwaaoHd83cn8VO0eBu1bxtdAbJ0g0oOtlJyca0dR2td81HScHdxnXStEvpZ74VgJGkyVvWIv1K4MVwzMzMFhGVi5jXc8cdd8wP7I455hiOPvpoNtpoI/r3b2tQ+MY3vsHZZ59NHiZXS70v9pC65q1O6p73VVKPpDdLJ3NXtNKskdWbHbtGecvJRyLi+m64ZnvvVXvq3oguLLfefWm0Th2pe1N5I+J2SSNI3RV3JbXUrk0aw3cAcI+kD0ZEIzOd1lLznub5M64hBXYTSePr/lrewizpPbR1G671Pnfk/iw2ik6oUgrWqnUhqKb0C9RC0742aNWy5/V+iSqdW0ZSS9GL5DU5gvTLyWukqYp3JP0CsGdEPFYn75T2HkXrY2ZmZt3j17/+NQC77747F110EZtuuukCgR2kMXsdIWkN2iYx+USesOTNimS1ZpbsUnmcV+l72qY9UYcCSjdiRBeVPzwv9VVL6Xvpa2XHSs9XaWedtVKLbXkLVylvPxb8zttI3kIiojUiroiIT0XEOqTg7izSd9/3k4Y4dZVtSQ0mc4EPR8StVboON/L5b+b+LHaKBnelZuHKgbi1jM7bdmdAqqE8UJtZJ135uSFNXGc6qR/667T9KvAScAJpWlozMzPrg158MY0gGTVqVNXzEcFdd91V9VwBa5Q9f7hGml07epEOKI2X+mjdVLWVuvR1dctJaQr+bSQt3wXlD6Bt/OEC8pJbW+fdv5edeihvB1IxC2lZXtE2Cc1DZaf+Qdv3zl3q1Kt07qGK402/73kmy9NoW3tu584qu4rS5//1WHCtvnKNfP6buT+LnaLB3d2km3xUewnzeLkvkz609QaJ9riI+E6kdTtWJK1ztx1pStrrgNslLVcn77D2Ht3zKszMzKyooUPTELN//vOfVc//9Kc/5dlnOzz0vrwH00KTzOWJVk6tPN6NSuuvbZbXkaupxnei0orwwzqzUlVcR/pBfiBl6+11sq+r+mw8XyAtGD6XtIg8AHmGyMfz7mk18h5K2+zyvy7LO4W2he1PqrZms6QdaWssuabidLvvezstXdDWQFJ57c68p6XP/0p5UpUFSFqd9pcAKSl0fxZHRYO7C0lv2gGSvl4rUb5xN5FmZJpDWqyyGa1lz+sN3i0/V2uGqIZExNsRcR9pIc+bSb9kXNaRMs3MzKx3Ki1zcOutt3LWWWfNnzRlypQpnHPOOXzhC19g+eU73Ej0GPBCfn6ZpPmz+eWJJsax8AyW3elG0vc2gB9J+qGk+UtdSRok6QOSfkxaT69SaV6Ej9b7QbyjcjB0ct49QtK1kjYpq+dwSQdJuqHJS8wkDcv5paRVc5mDJZ1A2/p6F0dE5bSppTrtAlyTu+EiaWCeffSn+fxvq8z8fgpprOXGwB9Ky3pJGiDpo7QFKveTlhgrV3rfl5H0sRqv6WuS/ijpk6XXVPa6Dgc+kw9VLuNRKnuHIkuN1fBX0mLmAn6TJxdCUn9Ju5M+/42MJ2z2/ixWCgV3EfEk8E3SzTlL0uOUzZwp6URJ15L+8Hcm3aiT6zTBtqd8nF29vsilc9NKi252VKRR0+fn3X3bWXPGzMzMFkGHHnooO+ywAwCnnXYaQ4YMYbnllmP55ZfnlFNOYY899uDoo4/u0DXyd4rPk77EbwI8JGmGpBmkroYbArW+nHe5XL+Dgd+SvuOdADwjaXqeVXEGaa21Y1lwyEzJz0jf+bYDXpP0sqRJkiZ1QV0vAk7P1zsIeFRSq6RW0lira2l+Hb7XgS8BnwReyq99KvBDUpfA8cBJVep0Uz4epPv4Qs47jdTAMZC0+Panq+T9O2k95Tmk786P5/kaWkktlcuRfhz4aOVsnBHxH9IMrAC/ljSt9L7ngAfSd/09gauA/+bPXemeXpbrdi/wrYqq/S6/H8vmOr1eVvY29d7EKq9xKvCVvPsB4ElJ0/NrvI00Q+vhDRTV1P1Z3BRtuSMizgbOIH2ANwCOpi3a/japv3apJe30iPhhB+r3GG19fjepk650rubEJ00qD0pHdHLZZmZm1sOWWGIJbr/9dk4//XTWX399llhiCSKCrbbaiosvvpibb755oQlWmhERt5C+2P6BtPzBANJ075cD/xcRf66du+tFxIyIOJC0TtjVpOUGBpCmpX+Z9CX8C6SJOCrz/gXYizQMZyppQr216Ph6eLXqeiZpuarLSQurl2YafZwUsDQ7dpCIuJAUDN1J+n77LqkV6yvAByOi6hwQEfFdYBtS18n/koLgGcBfgM8CO+cgp1req4DNgV+QZlMdCLxDWi/6q8CWdRpKPkoKbp4iTXhYet+H5fM/Iw2nupb0PXk2aX6KyaTA8DPAjpWNIxHxFunz+uv8eoaWlV14FvyI+CnpMzKOFNQNyOX+mNRV+V8NltPU/VmcqJ1pfWtnlDYFvkj6J1C+Zstk4E/AeRFRvQN7sev8lTR48s6I+GCV84NIsyctA5wdEZ2yzl0uew/g1ry7cUQ8Xi99LVtuuWU89FDlGNhFyziN67SyRsfoTivLrDeZNGsWaz/wAADPbb01IwYN6uEamXWth0enuUFGjVtwMpIJwybQMrKl4eP1yuqlPOW6WTeRNIYUyD8fESN6tja9X9F17uaLiEeBI2H+5CnDgOlVpjbtqCtJwd0ukrbIzdfljiYFdvNITc4NkTQgIt6td570awmkmTSfLFRrMzMzMzOzblS4W2Y1ETEzIl7ugsAOUvP6Y6Rfya6X9AGYP9D0SNoGUF4aEU+UZ5Q0VlLU6PP9Y0kXSxpdvjaepKUk7Uxq7i1NC3t6RMyrUoaZmZmZmVmv0HTLXXeJiHcl7UNahmFNYHwegLwEUJre9S7g+IJFDwIOI/VDjjywcw6pBbLUuX4OKbBrdrZPMzMzMzOzbtHrgzuAiHhG0mbAicD+pMlNZpEWf7wSuKSJlrVzgX+TplTdAFiRNFh0KmnR9XHALyKi2QXYzczMzMzMuk3N4E5Sh1fsLBMR8Z4OFjCVtBbIKQXyjCFNL1vt3BPAE8D3OlIvMzMzMzPrGhExFhjbw9VYZNRruRvRQP6g+oxRlcebm5LTzMzMzMzMGlIvuPtmnXOfAtYhjUkbT5rwpJW0psfGpHUxliItZv6rTqmpmZmZmZmZ1VQzuIuIqsGdpF+QArtfAl+KiDeqpFmetKDiIcAaEfHpzqmumZmZmZmZVVNoQhVJHwMOB34bEYfVSpcDvkPzAuOHS/pzRFzTsaqamZmZmZlZLUXXufscafxco5OQfJc09u4zBa9jZmZmZmZmBRQN7jbL20aXB/hPRT4zMzMzMzPrAkWDu6XzdpUG06+ct4MLXsfMzMzMzMwKKBrcPZe3RzSYvjSRyqSC1zEzMzMzM7MCigZ315PG0H1R0vH1Eko6ATieNEbvd03VzszMzMzMzBpSaLZM0gQphwBrAT+QdDRwAwuvc7c/sB4pEJyU85mZmZmZmVkXKRTcRcR0SbsAtwAbAusDJ1ZJqrx9HPhwRLR2qJZmZmZmZmZWV9FumUTEs8BI4EvAv/JhlT3Ix78IjIqI5yrLMDMzMzMzs85VtFsmABHxDnA+cL6kocAIUpfMVuD5iJjSSfUzMzOzRdjTJzxN68TUgWdu61xaJ7by8OiHF0pXSlPtHEDLyBbWO3+9rquo9QmSJpGGDx0eEWObyB/56U4RMa7zama2IEljgcOAKyJiTGeVW7jlrlJETI2If0bEPXk7pRPqZWZmZn1A68TW+YHb0O2H0jKypWq6lpEtNc+Vl1HPmDFjkFT10dLSwiabbMLRRx/No48+2vwL6iKSRks6Q9KYLih7TC57dGeXbbaokbRf/nvYr6fr0hWaarkzMzMza1TLyBZGjRvVdP5arXm19OvXj+HDh8/fnzt3Lm+88QaPPfYYjz32GD//+c+55JJLOPLII5uuUxcYDZwOjAfGdnLZY4Ad8/NxnVx2X/Nk3s7s0VpYV9qP3GIG3NijNekCHW65MzMzM+tN1lhjDV555ZX5j9dff523336bW2+9lXXWWYe5c+dyzDHHMGnSpJ6uqvUyEbFhfjzY03Uxa4aDOzMzM+vzllhiCfbYYw+uuOIKAN555x1uvfXWHq6VmVnncnBnZmZmi42RI0fOf97aWnsc3+23387+++/PKquswpJLLsnw4cPZc889ue6669q9hqStJF0t6UVJb0t6U9I4SZ+R1L8i7Yg8icfp+dCOkqLiMaYsvSR9QtLtkl6TNCeX/6Sk6yQdWZZ2TC671CXz9CpljyhLXzo2WtJKki6Q9LSkWZKmlKUbKunjkq6S9K98/dmSns+ve5s6780Z+Rrj8v7HJU2Q9JakVkkPSDpckmqV0ShJy0j6TtlreE3SbyRtWifP/PegyrkNJH1V0p2SnsllTpP0sKSzJa1Qp9wBko6W9BdJb+T7NlnS45KulPSRDr7WTSX9NH8OZkiaLukxSZdL2r1GnmUlnSXpnzn9zFyfH0hatUaeEeWfHUnrSrpU0qT8WZ+Y040upcv7W0i6RtJ/Jb0r6cYqdTld0kOSpuTP03O5/jXvV87bT9InJd0i6X+5Hq9Iuj+XuXZ5nUhdMgEOq/L3MLpK+e+TdEWuz+x8zx+SdLKkIe3UbT+lv/2pOd/fJB3VGZ/vWjzmzszMzBYb//znP+c/32CDDRY6HxEcd9xxXHjhhQBIYtiwYbz55pvcdttt3HbbbQC/AQ7Js4cvQNJXge/QtjzUFGAIKcDaEThE0j4RMTWfnwu8Spp1fGlgDvBmRbGzyp5fTtuXU4BpwCDS2sPrAx8FflGW71VgOWAJYAZpZvNycxd6E2Bd4FpgRWB2rlO5L9IWjJLLFLBmfnxc0gkR8aMqZc8n6bvAV4EgvU+DgK3yY2dJh0XEvHpl1LEs8CCwAfA28A4wHDgQ2FfSRyLiloJl/ok0Eye5zlOBoaQlwkYCYyTtEhFPlmdSCuhvA3YpOzwFWAZYnrR29AeA3xWsT6n8U4CzaPvMzcrPN8qP/YFhFXk2z3VapSzP3FyXDYEj8+d0fJ1LbwdcQvrszgDerVG/jwDXkD6D0yrTSdqONPatNFD2XdLnbgRpvOjBkj4XEZdXKXs4cAPw/nyodF9WAFYCtiZ9Fk4gfQZeJd2zgfkaUxcskQX+piV9C/h62aEZpM/pFvlxuKTdqy39Juk7tK0HXqrXKOBi0hjb2ZV5OoNb7szMzKzPmzNnDrfffjtjxowBYOONN+ZDH/rQQul+9KMfzQ/sjjvuOF599VXefPNN3nrrLc466yzyD+4HAt+qzCtpX+C7pC/WvwHWjIhlScHd0aQg4wPApaU8EfFiRKwMnJcP3RsRK1c8rs3l70AK7OaRgqJlI2JoRAwiBWL7ANeVlX1tLvvefOi8KmW/WOXt+iHwBukL6OCIWAZ4X9n5V3KabXIdhpC+LL8HKAV0P5BUbxadkfk1nA8Mj4jlSF/Iv5vPHwIcXyd/e07P5X0EWDq/hpHA34AlgaslrVmwzPuBL5CC34H53g4EdiUFkqsBV1fJ90lSYDcLOCLXZ1lgKVJwdTBwe8G6ACDpWOBs0mfuOmCziBgcES2k4GZf4NaKPMsAN+drPw/slus0hHRP/0UKPG+QtHqdy18CPAJsEREtEbE06W+j0thch/UjYigpOPpKrss6wB9Jgd01wHtJ7+0QYA3gQlJQ+DNJW1W8jiWAm0iB3UxSALdC2X1Zn7Qu9/MAEXFv/nu4NhdxbZW/h3vLyv8yKbCbkuu7Yn5fB5P+jv8BrAdcL2mBmErSgbQFdlcAq+Z6LQucAhxEujedzi13ZmZm1qe8+OKLrLzyyvP3S7NlRgTDhw/nc5/7HOeccw4DBiz4NWjmzJl885vfBODwww/nggsumH9umWWW4dRTT+Wdd97hrLPOAjhe0g8j4uWyIs7N2z8Dnyi1OkXEbOCnkt4htap9VNL7IuJvBV/atnl7R0ScV34iIl4Hfp8fHTUX2CUi/ldW/tNlz39amSEiAniW9L70Bz6fH5+ucY2hwNiI+GJZGW8BJymtofw54BuSLs7vX1FDqVirLiL+KWk34HFgZeBkUtDdkIj4eJVj7wB/lrQL8B/g/yRtHxF/LUtWum9Xlrc+5ffsFVJAWC0orEvSssA5efeyiFhg+teIeI0UxN1ckfUYUgvkbGC3iHiqLM8D+bU8TmpVPJl0H6uZnPPPKMv/dJV0jwIHRMTcnCZI7xXA90j36ucR8dmK+r8EfEHSUsBnSIHWfmVJxpDe23nAhyPi7rK87wJPk36EKCy3CJ5F+lvYJyImlJU9B5iQP0v/Iv1osA8Lzrx5dt7eSlpzMXLe6cA5kgYBpzZTt/a45c7MzMz6lHnz5vHqq6/Of0yePJn83YrW1lZaW1uZPn36QvnuuOMO3nrrLSRx+umnL3Qe4MQTT4TUNWsJUhdIACS9l9SdDeDMGt0JLwcm5eefaOKllbqQraiKsXud7JflgV0T/pC327eT7uw6x4PUyrFbk3X4a7VFyPN6zBfm3YM6a+xTRLSSlrGAhV936b6tQuc6kNQqPJu2VqJGfCxvrykP7EryDwWl9+jjdd6jC8sDuzrOKwV25SQtR1uwdm7l+TK/zNtdKz73h+ft78oDu05yMKmF8d7ywK5cRLxBW6vo/HGNkkaSWg0BvlUK7Cp8jwW7W3caB3dmZmbWp6y11lpExAKPadOmcffdd7PVVltx1VVXsc022/DYY48tkO+hhx4CYN1112WttdaqVjQtLS2QuudBGnNTsmXezgbuqZY3f8n7c5W8jbqT1LVzFPAXpQlTinYtbETV+peTtI6k8yT9PU+AMbdsAo0/5mT1uvS9GBHPVDuRW2xKLTvNvE8A9b7sl84tRxrX1TBJH5Z0raRnlSYuibLXfVBOVvm6S+/HPpL+KOljklam40otgvfnQKOR+i8JbJZ376yTtHSu3nvU7ueknXTb0haL3JsnQVnoAVyf0yxNak1E0gDa/ub+QOcrjeHbsla9ct1KgXL5P4zSZ3Y2bf8rFhAR04C/d0G9iwV3kn6cf5kyMzMzW2QMGTKE0aNHc8sttzBixAheeeUVjjrqqAXSvP766wCsttpq7RX3Ut6uWHasNBnEK9VaKdrJ25AcDH2W1HK4Hakl8HlJL0v6paQPFi2zhtfrnZS0P/AY8GXg/0jd6lqB10gTVryVky5dp5j/tlOH0vnC71MD5Zefa6h8pRkZryZ1ez0IWJs0du8t0mt+lbYJMhZ43bnl52TSxDR7Ar8G/pdnX7xE0tZVrve3GgHFBWXJSgHi8428hmw5oNT6Ve89eqnsea33qO7npIF05TNyrlTnUT4L6eC8XZ7Ueg7FXn+jSnUb1E7dSvd6cFneRv8XtPc30JSiLXefB/6Rp//8XHvTf5qZmZn1Ji0tLXzsY+nH9gkTJvC//y3c+7BAT71q3a2qHWs0b/uZIq4ktaQcQ5pA42VSd79DgNsl/bZycocm1PxCKml50gQZSwF30TbpytCIWClPWFFtUo1K7b3+jnaXrFd+M2UfSepKOxc4kzSRxlIRsVxpMg7gt7XKj4hzgXVIE3P8njRhzQhSsH6/pB9UZBlO9WBiaJW6NfVZKpCvVrp6gUtb5toBTinInB4RavAxqUD9OqJUtx83WK/R3VSvdhX9459G+sD+H/AT0q8Ol0tqr0+1mZmZWa+w5pptPRknTZo0//mKK6YGihdfrDaB5AJK3e7KWyRey9tVcpexInkLiYjJEXFxRBwUEasBG5O+l0GaHfIzzZbdgA+RZlJ8C9g7IsZHROXYoUa6HLbXPFpqOXmtbqrmyi9vMWq0/NJkKpdGxOkR8Z8q4yrrvu6IeCkivh8R+5CCt/eRZlUF+KLK1qOLiBE1gogxZUW+krcjGnwNkJbZKAVba9RJV961tOnPajtK9R/SRPfiN2hbomNEp9WoTaluddfYq6H0fq3cztjYdrsINKNocLcycCgwLu8PzvvjJT0h6St5dhkzMzOzXumll9p6nC29dFsPui23TEN4nnnmGZ5/vnpPrxkzZkBaOwvgobJTpecDSV0mF5InptipSl5IM/5BE61KEfF4RHyetvF8O3dW2VWUAoInI2JmjTS7NlDOmnka/IVIWo203AA0Py5pdAPn3qRtgpv2lF73w9VOSmqh7XPRrkgeIrUGlmaYrLxv7SlN279NblFt5LrvkJYvgAXX3atUOvdGVFnDrZPcR1vrViOtvfPl2TBLs81+uOB1G/l7KI0TfL+kohPhlD6zpXUbF5J7PzY7nrSuQsFdRMyOiF9FxM6kP7pzSN0BRJoV5jvAS7lLwJ6dNQORmZmZWWd45513uOGGG4AU2JUvZP7BD36Q5ZZbjoiYvyRCpfPOOw/Sos1zKFt0OiIeIU0fD3Baje9Ah5K65kEad1VuWt4Oq1X3PBlGPaVga2DRsgsozfy4vqTK65RmCvxkg2V9vc5xkVoH7yhawWyHvC7gAvIyC6Wp/a+rMZNhNaXXXWvuiW+QZq5cSL37llv/SmP1Fno/23EdMD3n+247acuVWgs/IWn9ypOSVqDtPar8nHaaiHiVtqU7TpZUfRajtnotV3Hosrw9QNLoApdu5O/hV6T7siRwQb2YRtISObgHICIm0hawn1Ij25dJwV+na7pPdkQ8FxGnkmaH2Yu0Ovy7pMGN+wO3kAb5ntHezTIzMzPrak888QSf+MQneOKJJwD4/Oc/z1JLLTX//KBBgzjttNMAuPzyyzn++OOZPHkyANOnT+ecc87hzDPPLCW/IBZc4w7SpBmQWj2ukbQGgKSBkj4LlNaH+20svMbdo3m7iaT3U91FeabG/fMXcHL5QyV9lfR9DNpmZ6ws+0O5Vawjbie1fCwHXFUqT9KSkg7K5xdeZ2JhU4EjJX2/9KVd0jBJ55DGE0KaRr7Z6eKnkRaX3q/UNU7S5sCfSN0yW4FvFyjvtrz9jKTPlgI2SStL+iFpKYJaM1beKOkXkvbIwSU570qSvkfb7JWV962uvKxD6TN3RP5sbFJW/nBJB0m6oSLrT0iTkAwE/iRp11LwIul9pJkyh5OC63pLFHSGL+XrLA/cJ+lT5YGSpFUlHSzpdlIjUrkrSa1//YBbJB2ntPYfkgZIWk/SNyR9pSJf6e9hB0kbUkVeCqT048OBwB8lbVX2PvWXtLmkU0kzu46sKOIbebuXpMuUZ0eV1CLppHx+SjvvTVM6vIh5/sXhVuDW/I/mUOAIUv/v1UmVP1XSn4GfAzfmplQzMzNbDLRObOXh0VV7szWcv2VkS/sJs8pFzCGtb5e7VALw0Y9+tLQY+QKOO+44nn76aS666CJ+9KMfceGFFzJs2DCmTp3K3Lnz54X4LVV+kY+Im/IXt3NJU6R/TNJbpJa+0sx+E6i+sPc44ElgA+CvkqbQ1lr0lYj4bS7joPxA0nRSoFU+ycYNwKUVZV9BailYF3hB0uu0tRZtn5ceaEhEPJ0DkpOAA0itJlNJQ3WWAJ4jLc58VTtFTQQeBL4KnJDLGEpbw8OvaXIB6uybpIXQbwBmKy0gv0w+Nwc4JCKKzLL4fdK6hhsClwAXS5qW66x8bCBwWJW8g0nfjY8AyK+1P+lzUfKjiPhTgfoAEBEX5S6ZZ5A/G5JKH/RSn+OpFXmmSdqHFOiOILWOzpQ0r6xO00gLjzf82WhGRDyjtBj49aSur1cC8/LnfyALzkJ5aUXeOZL2A24CtgEuAM7PeYfQFueUzzAKqcX9HFIA+7ikyaQZaAE+HhH35/J/mFunzwb2yI+3JbWS7nt5HLVAC3BEXCtpK1LwejgwJtdrGdK9v5b0N1jt89IhnbrOXR7g+wPgA8DVFdfZlfRCXpT01Qa6FpiZmdkirmVky/zAbOpfp9I6sbVqutaJrTXPlZfRiMpFzF999VXmzp3LiBEjOOigg7jlllu47rrrWHLJhb+KSOLCCy/ktttuY5999mGFFVZg+vTpDBs2jN12241rr72WiDgwj11aSER8l/RF8xrSVOctpC+OfyHNjLhzREytku9dUovfpaRxYINJvaPWou0L91nACaQvs0+SvlAuTZr84Q/AQRFxQOXshBHxNGms382kyR6WLyu78A/9EfE10o/5D5IWYl6C1HpxDmkNvsoWzVrlnEjqwllqfZlNGkf1aeCTVSYsKeIt0nin75Gm9V8KmEwKzLeMiJuKFJZbybYDzifdn7mkHmvjgE9ExFE1sgJ8AfgaqfXvP6T3fEngxVyf3SLi+CL1qajbmaTJDi8nBdf9c90eJ3Vd/GiVPI8Am5ACl3+RPksDgKfya9woqiwC3xXy2MONSJ/tcaSxkMvkOj1Oel37kt7HyryvkRaNH0MKUt8g/b1MBh4ATiO9nvI8b5FilV+T/kaH0vb3MLAi7bdJ79OFpOU/5uT0U0hjHs8ERkXEQmv5RcSXSRMcTSD9D1iC9KPGMaSxll1CjXc1bqAwaRfSVLH7kf6ISv1THyW9qNJA08jHdomIyZ1WgV5qyy23jNLCqIuqcRrXaWWNrjpbrNmib9KsWaz9wAMAPLf11owY1CXd6c16jVJr3KhxoxY4PmHYBFpGtjR8vF5ZvZTnFGiSpDOA04HxNaaPN7MO6HDLnaTVc3/WZ0n9rD9GinpnkdZB2S4iNo+IbUl9iseSgrtNSb9AmZmZmZmZWQc1NeZOaf2W/UitdLuSgsTSr1gTSWPrfhURCwyojYh/kwZ83k8aVPyhpmptZmZm3eLpE56u25WySHdJMzPrWoWCuzwDz5HAIaT+2pCCulZSv9Wf5X6z7bmKFNx1yeJ91pzO7HppZmZ9Q2ksXLUgbm7r3Co5zMyspxRtuSsNuCy10j1EaqW7OiJm1MxVISJmlGYSLXh9MzMz62a1xslNGDahB2pjZma1NNMtczqp5e3neZG+Zu3UgbxmZmZmtoiJiDNI0/abWRcoGtwdCfy6AwtKzhcR4ztahpmZmZmZmSWFgruIuLyrKmJmZmZmZmbNK7QUgqR5kt6VNLj91CCpfylPc9UzMzMzMzOzRjSzzl0zk6B44hQzMzMzM7Mu1OFFzNtRCuqii69jZmZmZma2WOvq4G7FvG14mQQzMzMzMzMrrtngrm5LnJLVgNPzoWeavE55mUMlnS3p35JmSHpL0r2SjpJU+HXk8YC7Svq+pPslvSlpjqQ3JP1F0hclLd3RepuZmZmZmXWHurNlSppb7TDQmhchb0QANxesV2U91gHGAWvkQzOApYFt8+NASXtFxOwCxf4U+HTZ/lygFVgO2CE/jpW0Z0Q81ZH6m5mZmZmZdbX2WrxU8ah1vN7jbuA7zVZQ0gDg96TA7iVgdES0kIK7zwJvAzsDFxQsegng1Vy3rYGBETEMWBb4KjATWAf4o6SBzdbfzMzMzMysO7S3zt03K/ZPJ7XEnQu8UyffHGAy8GBETGy6dskRwMb5ugdExN8AImIO8HNJQ4HvAZ+W9MOIeKLBci8Gjq5ckD0ipgDnSXoJuAZ4D3Ag8MsOvg4zMzMzM7MuUze4i4gFgjtJpTF034qImV1WqwUdmrd3lQK7Cj8BTgOGAAcD32ik0Ih4oJ0kvwF+lsvdEgd3ZmZmZmbWixWdiGQnUhfIWe0l7Ax5sfRt8+5t1dLkIHNC3t21s64dEfNILZAA/TurXDMzMzMzs65QKLiLiPH50V3r1m1MWx0frZOudG7jzrqwpM1Ik6sA/KtOuintPTqrTmZmZmZmZrV09Tp3HbVq2fOX66QrnVtGUksnXfvbeTsVuK6TyjQzMzMzM+sSNcfcSTqt9Dwizqw8VlSpjILKA7V6Y/zKzw0hLWnQNEknAHvl3a9HxJu10uYZNuvacsstu6ul08zMzMzMFlP1JlQ5g7bFys+scqyoZoK7bidpH+C8vHtNRPykJ+tjZmZmZmbWiHrB3QssHMhVO9aVylvgBtdJV35uerMXk7QbaZbM/sAtwGHNlmVmZmZmZtadagZ3ETGikWNdrHyc3arAIzXSlcbmTYuIprpkStoFuBFYCvgT8NG8lp6ZmZmZmVmv19snVHkMmJefb1InXencY81cRNJOwM3AIOAOYL+IeLuZsszMzMzMzHpCrw7u8hp29+XdPaqlkTQI2CHv3ln0GpJ2JHXBHJzz7xsRs4vX1szMzMzMrOf06uAuuzJvd5G0RZXzRwPLkFr4ripSsKQdgD/QFtjtExHdskC7mZmZmZlZZ6q3FMKanXmhiHihyayXAceTFii/XtKnIuIvkgaQJjw5J6e7NCKeKM8oaWxO83zleEFJ25ACu6VxYGdmZmZmZou4erNlPteJ14l2rlU7Y8S7eXmCu4E1gfGSZgBLAEvmZHeRAsAiziGtiQfwf8BzkmqlvTciDihYvpmZmZmZWbepF3DVjHS6W0Q8I2kz4ERgf2AEMAv4B6nb5iURMa92CVWVd0ldrp207Z03MzMzMzPrUfWCu526rRYNiIipwCn50WieMcCYGudGd0a9zMzMzMzMeoN669yN786KmJmZmZmZWfMWhdkyzczMzMzMrB0O7szMzMzMzPoAB3dmZmZmZmZ9QL117u7KTyMidqk4VtT8MszMzMzMzKzz1Zstc3TeRsWxoPgyCdF+EjMzMzMzM2tWveDuiirHrsSBmpmZmZmZWa9TbymEw6scG9OltTEzMzMzM7OmeEIVMzMzMzOzPsDBnZmZmZmZWR9Qb8xdQyRtCowAhgDTgUnAvyPCY/PMzMzMzMy6SVPBnaRhwKnAGGDZKkneknQ58K2ImNJs5czMzMzMzKwxhbtlShoJ/Bv4IrAcaVmEysdywJeARyW9t7Mqa2ZmZmZmZtUVarmTNBy4A1g+H7oL+B3wGNAKtAAbAx8BdgZWBW6XtGlEvN5ZlTYzMzMzM7MFFe2W+XVSYDcbOCQirq+SZjxwsaQDgKuAFYCTSS15ZmZm1gs9fcLTtE5sXeh468RWWka2dPn1Z/1nFnNb5/Lw6IcXOtcysoX1zl+vy+tgZraoKxrc7U1axPzsGoHdfBFxvaSzgbNyPgd3vcQ4jeuV1x8do7u1HmZm1qZ1YmvVQG5u69xuuX6t61QLOM3MrLqiwd1qeXttg+l/TQruVmsvoZmZmfWslpEtjBo3aoFjE4ZN6NHrV2vJMzOz6opOqDIlb99qMH0p3bSC1zEzMzMzM7MCigZ3pZ/PNmkw/aZ5+4+C1zEzMzMzM7MCigZ3PyYtdXC6pLp58/nTSWP0LmyuemZmZmZmZtaIQsFdRNwKfJO0zMFtkjaqlk7SxsCfgNHAmRHxxw7W08zMzMzMzOqoOaGKpMvq5JsE7EpapPwpFl7nbv2ydGtJ+kVEHNkZFTYzMzMzM7OF1ZstcwypS2UtQeqiuQFtwRz5WOn8iPwAcHBnZmZmZmbWReoFdy9QP7gzMzMzMzOzXqJmcBcRI7qxHmZmZmZmZtYBRWfLNDMzMzMzs17IwZ2ZmZmZmVkf4ODOzMzMzMysD6g3oUpDJA0GhrVXVkS80NFrmZmZmZmZWXVNBXeSNgS+AuwOrNpAlmj2WmZmZmZmZta+wgGXpAOBK4ClaFvTzszMzMzMzHpQoTF3kt5DCuwGAr8F9sqnIj8/APgm8J98/FFgD2DnzqismZmZmZmZVVe05e4LpMDu7oj4GIA0v/FufETMBG6UdDZwLvBl4MsRsXsn1dfMzMzMzMyqKDpb5s6kVroL6yWKiLkR8VXgemBXSUc2WT8zMzMzMzNrQNHgbq28/VeVc0tVOfZD0ri8Qwpex8zMzMzMzAooGtwNytvJZcdm5O2wKukfz9uNCl7HzMzMzMzMCiga3L2Vt8uXHXs1b9evkr6UbljB65iZmZmZmVkBRYO7J/N2pbJjD+ft3lXS75O3rxe8jpmZmZmZmRVQNLi7O2+3LDv2W9K4us9KOkbSMpKWzZOonEGagOWODtfUzMzMzMzMaioa3P2eFMjtV3bsN8CDpGUVfkzqujkZ+BmwNGlM3jkdraiZmZmZmZnVVii4i4iHSIHd98uOBfAh4I+kwK/88TSwe0T8Z6HCzMzMzMzMrNMUXcSciLi5yrE3gQ9LGgGMIi2L8CzwUETM62glzczMzMzMrL7CwV09ETEJmNSZZZqZmZmZmVn7io656zGShko6W9K/Jc2Q9JakeyUdJamp1yFpQ0lHSLpQ0n2SZkoKSdHZ9TczMzMzM+tKHWq5k7Q8sC0wAhgCTCe13N0XEW90tHJl11kHGAeskQ/NIE3Wsm1+HChpr4iYXbDonwI7dlY9zczMzMzMekpTwV0Otr5LWseuf5UkcyXdBJwUEc92oH5IGkCapXMN4CXgkIgYL2kJYAxphs6dgQuAzxUsfi7wOPCP/BgOfK0j9TUzMzMzM+sJhYM7SbuT1rYbTJoRs1a5BwB7SDogIjqyzt0RwMak9fIOiIi/AUTEHODnkoYC3wM+LemHEfFEgbJ3i4i5pR1JH+9APc3MzMzMzHpMobFqktYGrid1iZwLXAbsBaxNavVaO+9fns8vDdyQ8zXr0Ly9qxTYVfgJqTtoP+DgIgWXB3ZmZmZmZmaLsqITkXwdGARMAbaPiE9HxK0R8XxEvJG3t0bEkcAOOd0g4ORmKidpMGlMHcBt1dJExExgQt7dtZnrmJmZmZmZLeqKBne7k7pHnhoRD9ZLGBEPAN8gdd3cvbnqsTFtdXy0TrrSuY2bvE7TJE1p79HddTIzMzMzs8VP0eBuxbz9fYPpSwuer1g3VW2rlj1/uU660rllJLU0eS0zMzMzM7NFVtHgbnLezmowfSndmwWvU1IeqM2sk6783JAmr9WUiBjW3qM762NmZmZmZounosHdvXm7RYPp31eRz8zMzMzMzLpA0eDuh8A84BxJS9dLmM+fA7yb8zWjtez54Drpys9Nb/JaZmZmZmZmi6xCwV1E3EdaKHwT4D5Je0paoAxJ/SR9CLgf2BA4KiKabbkrH2e3as1UbeemRURrnXRmZmZmZmZ9Us1FzCXdVSff68CmwC3ATElPkVrZWoD1aWtJ+y9wiKSDI2KXJur3GKmlsB8poKy6HEI+V0pvZmZmZma22KkZ3AGjScseqOJ4lD0XaaHyUTXKWA1YvSJPwyJipqT7gPcDewDfr0wjaRBpTT2AO5u5jpmZmZmZ2aKuXnD3F5oMyjrZlaTgbhdJW0TE3yvOHw0sQ2rhu6q7K2dmZmZmZtYb1AzuImJ0N9ajnsuA40kLlF8v6VMR8RdJA4DDSJO2AFwaEU+UZ5Q0Nqd5PiJGVBYsaSkWXDphSNm5FcrTRsRkzMzMzMzMeql6LXe9QkS8K2kf4G5gTWC8pBnAEsCSOdldpACwqE8Al9c493rFfmX3VDMzMzMzs16j6FIIPSIingE2I7XSPU6q9yzSjJzHAB+MiNk9V0MzMzMzM7Oe1eGWO0mrAyNIXRqnA5Mi4qWOllspIqYCp+RHo3nGAGPqnB8LjO1YzczMzMzMzHpeU8FdHqt2AmnNu7WqnJ8E/BS4ICLe6UD9zMzMzMzMrAGFu2VKWg14kNRFcgRpLFrlY23gXOBvkuotPm5mZmZmZmadoFDLXZ6h8hbS+DeAfwPXAo+SumQOIS1uflDebgrcIul9ETG3syptZmZmZmZmCyraLfMI4L2kNeW+EhHnV0lzI3C2pOOBH+T0hwOXNl9NMzMzMzMzq6dot8yDSAub/7hGYDdfRFwA/JjUTfNjTdXOzMzMzMzMGlI0uCt1x/xFg+l/nrebF7yOmZmZmZmZFVC0W+awvG10qYOX83ZoweuYmZmZ1fX0CU/TOrG16rmWkS2sd/563VwjM7OeVTS4ewsYTlr+YEoD6dfI20bSmtU0TuOqHh8do7u1HmZm1nu0TmyldWIrLSNbFjpuZrY4KhrcTQQ+CHwW+HwD6T9Xls/MzMysU7WMbGHUuFELHHt49MM9VBszs55VdMzdtaQJUo6S9DVJqpZIycnA0aQJWK7pWDXNzMzMzMysnqItd1eQArYtgW8BR0r6HfAYbevcbQIcAKyT8/wNuLJTamtmZmZmZmZVFQruImKepA8DfwC2AN4DfLVK0lKL3t+AfSIiOlRLMzMzMzMzq6tot0wi4jVgW+ALwCP5sMoe5OPHAu/P6c3MzMzMzKwLFe2WCUBEvAtcBFwkaShp9swhpK6Zz0fE1M6ropmZmZmZmbWnUHAn6a789JyIuBMgB3KP1M5lZmZmZmZmXa1oy90OpK6ch3dBXczMzMzMzKxJRcfcvZq30zu7ImZmZmZmZta8osHdQ3m7UWdXxMzMzMzMzJpXNLi7mDQj5oldUBczMzMzMzNrUqHgLiL+BHwb2FvSVZJW7ppqmZmZmZmZWRFFZ8u8LD/9L/Bx4EBJE4FngZl1skZEHNlUDc3MzMzMzKxdRWfLHANEfq6cf4v8qEU5j4M7MzMzMzOzLlI0uHuBtuDOzMzMzMzMeolCwV1EjOiiepiZmZmZmVkHFJ0t08zMzMzMzHqhhlvuJK0ErE4aQ/ffiPhfl9XKzMzMzMzMCmm35U7SRyT9E3gZeBB4AHhJ0mOSDu7qCpqZmZmZmVn76gZ3kk4CfgNsSmqxK39sCFwp6VtdXUkzMzMzMzOrr2a3TEmbAWeTArl5wJ3ARNJsmZsBuwH9ga9J+n1E3N/ltTUzM7PFyqz/zGJu61weHv3wQudaJ7bSMrKlB2plZtY71Rtz9zlS8DYF2CMiHiw/KWkkcDuwPHAU4ODOOmScxvV0FczMrJeZ2zq3qXNmZoujesHdDqRWutMqAzuAiJgo6WvApTmtmZmZWadrGdnCqHGjFjo+YdiEHqiNmVnvVW/M3Vp5+4c6aX6ft2t0TnXMzMzMzMysGfWCuyF5+1qtBBHxen7aX9LATquVmZmZmZmZFVIvuFPeRieUZWZmZmZmZl3IAZmZmZmZmVkfUG9ClZLtJc3uaLqI+Evj1TIzMzMzM7MiGgnubm3nfKnbZr100eC1zMzMzMzMrAntBVxq57yZmZmZmZn1AvWCu8O7rRZmZmZmZmbWITWDu4i4ojsrYmZmZp3j6ROepnVi60LHZ/1nFgCD1h200LnWia20jGzp8rqZmVnX8Ti4Rdw4jevpKpiZWS/TOrG1arD2zivv0L+lf9U8c1vndkfVzMysCzm4MzMz64NaRrYwatyoBY5NGDah6vHSOTMzW7QtMuvcSRoq6WxJ/5Y0Q9Jbku6VdJSkDr0OSbtJulnSq5JmS3pe0i8kbdhZ9TczMzMzM+tKi0RwJ2kd4F/AKcDGpKUVlga2BS4G7pA0sMmyzwH+BOwNDAfeBtYEjgAelrRPh1+AmZmZmZlZF+v1wZ2kAcDvgTWAl4DREdFCCu4+SwrGdgYuaKLsTwIn590fA8tFxFBgXWAcMBD4taT3dPBlmJmZmZmZdaleH9yRWtBKrXUHRMR4gIiYExE/B07N6T5dpBtlDhrPzbu3RMRxETEll/0MsC/wMjAIOLMzXoiZmZmZmVlXWRSCu0Pz9q6I+FuV8z8BppNey8EFyh1Nag0E+E7lyYiYRuryCXCAJM8PbWZmZmZmvVavDu4kDSaNqwO4rVqaiJgJlKb42rVA8aW004D7aqT5U94OBN5foGwzMzMzM7Nu1auDO1J3zFIdH62TrnRu4wJlb5q3j0dErcV9yq+5SYGyzczMzMzMulVvX+du1bLnL9dJVzq3jKSWiGgtUHbNciNilqQpwLCKusyXz9e1xRZbNFAdMzOzztM6sZWHRz+8wLG5rXOrHq93rjPzdNf1qy3gbma2OFBE9HQdasqzWV6Vd9eLiP/USPcZ4Gd5d9WI+F8DZT8FrAdcFRGH1En3X1Jg97OI+FyV81PauxYwB3i+gXQGI/N2Yg/WwZo3Mm8n9mAdrGNG5u3EHqyDNW9k3k7swTosDiZHxB49XQkzs0q9veWu14uIYT1dh76kFCxHxJY9XBVrgu/fos/3cNHm+2dmtnjr7WPuyrtXDq6Trvzc9IJl1yu3/Hyj5ZqZmZmZmXW73h7clY+HqzrmreLctAbH25WXXbNcSYNI4+0q62JmZmZmZtar9Pbg7jFgXn5eb7bK0rnHCpRdmglzI0m13ofya/67QNlmZmZmZmbdqlcHd3kNu9IadFUHLufWtR3y7p0Fii+lXYa2tfQq7Z63s4F7CpRtZmZmZmbWrXp1cJddmbe7SKq2psDRpABtHm0zazZiHPBifn5i5UlJQ3LZANcX6O5pZmZmZmbW7RaF4O4yUndLAddL+gCApAGSjgTOyekujYgnyjNKGispJE2qLDQi3gW+lnf3kXS+pKE53zrATcBqwCzgtM5/WWZmZmZmZp2n1y+FEBHvStoHuBtYExgvaQawBLBkTnYXcHwTZV8taVPg5Jz/C5KmA0NzktnAxyPimQ6+DDMzMzMzsy61KLTckYOrzUitdI+T6j0LuB84BvhgRMxusuyvk8bW3QK8AQwCXgAuB0ZFxM0dfgHWsIgY5rUDF12+f4s+38NFm++fmdniTRHR03UwMzMzMzOzDlokWu7MzMzMzMysPgd3ZmZmZmZmfYCDOzMzMzMzsz7AwZ11O0kX5yUqQtK4dtJuIOlSSc9Lmi3pVUm/l7R7vXzWOSSNKbtX9R671inD97CXkLSspFMlPSjpjXw/XpD0J0lfk7REjXy+hz2kwb+/0uPyOuX4HpqZLQY8oYp1K0nbA38hrVsIMD4iRtdIuw9wLTAwH5oGDCnL++0826l1EUljSDPHzgNer5P0wIiYUCW/72EvIemDwNXACvnQO8BMYFhZsmUjYkpFPt/DHiTplXaSLEXbPTw2Ii6qUobvoZnZYsItd9ZtJC0FXArMBR5qJ+26wK9JX0bGAetGxFBgOeDHOdnJkj7RZRW2ci9GxMp1HtUCO9/DXkLSB4DfkwK7W4FtgIERsSzQArwf+AEwpyKf72EPa+fvbmXS/1RI67JeXZnf99DMbPHiljvrNpK+BXwd+A6wMnAYNVruJF0NfAJ4GdgoIqZVnP898GHSmoTrRMTcrq394qms5e75iBhRMK/vYS8gaWngUWAEcBXwqWjwH7/vYe8maQDwErAScHVEHFwlje+hmdlixC131i0kbQ6cCDwLfLOdtC3A/nn34sovI9l38nZNYHQnVdM6ie9hr3IoKbCbDny+QGDne9j7fZgU2AFcVnnS99DMbPHj4M66nKT+wC+AAcAxETGrnSzb0zY25LYaae4jjRsBqDmZh/UY38Pe41N5+7uImFogn+9h73dE3k4C7qpy3vfQzGwx4+DOusMJwJbANRHxpwbSb1r2/N/VEuSuQ0/k3U06VDtrxHBJf5fUKmmWpGcl/UrS6BrpfQ97gTzOdYu8Oz7PmHiFpJclvS3pv5KuyxMdVfI97MUkrQzsmXcvr9Ei63toZraYcXBnXUrSOsCZwBTgiw1mWzVv32qnle/livTWdQYD/0eaYbEfsDZwMHC3pMvy2J9yvoe9wwhgyfx8XWAiqZvmcsAs0vv+UWCCpG9U5PU97N0OI/WGmAeMrZHG99DMbDHj4M662s9IgcGJEfFqg3la8nZmO+lK54c0UzFryMukMZLvJc2uuBzpfr4fuDOnORz4YUU+38PeYdmy518HWkljsFoiYhiwAW338UxJB5Sl9z3s3Q7P2zsi4oUaaXwPzcwWMw7urMtIOgLYBbiHtum6bRESEbdHxBkR8UhEvJ2PzY2Ie4HdgZty0mMkrddjFbVayv/HCzg8Im6MiHcBIuIpYD/aWm5O697qWTNyN9oN8u5CE6mYmdniy8GddQlJKwHnkdbN+myjM/RlrXk7uJ10pfPTC1bPOkFEzAO+knf7AXuXnfY97B3K39enIuKWygQRMQP4Sd59bx7LBb6HvVlpIpU3aPuBpRrfQzOzxYyDO+sq55K6hP0YeEFSS/mDNFYEoH/Z8dKxUivCspIG1blGaXzIy3XSWBeKiP8Ak/PuOmWnfA97h/+WPX+iZip4vOz5mnnre9gL5f+fB+XdX5Va1GvwPTQzW8w4uLOusnbefon0a3Dlo7TY7vZlxw7Jxx4tK2fjaoVL6gdsmHerzgJnPcr3sBeIiDdp7Au7yrPlre9h7/QxYOn8vL0umb6HZmaLGQd31hv9FZidn+9RI822wDL5+Z010lgXk/QeYIW8+1zZKd/D3uP2vN2oTprSuSCtmQa+h71VqUvmQxHxSDtpfQ/NzBYzDu6sS0TE6IhQrQdwRU46vuz42Jy3Fbghnz9KUrUZ3E7M2xeAcV34UhZbktTA+e/k3XnA/PFcvoe9yti8XU/S3pUnJQ0Gjs67D0bE6+B72BtJ2gDYLu/+or30vodmZosfB3fWW51GWodrdeCmvF4ekoZKOh/YJ6c7OS/Ca51vLUkPSvqcpHVKwZ6kfpK2Av4AfCSnvSQinqzI73vYC0TEeNq+4P9C0r6S+gNIWp80IceqpAD9lIrsvoe9y5F5Owu4psE8vodmZosRFZvE0KxzSBpLWoR3fESMrpFmH+BaYGA+NJW0DlPpR4lvR8TXu7amiy9JI1iwq+XbpLGRQ4Clyo5fTpoR9d0qZfge9gK5xeY22lp9ZufHsLz/LvCFiPhplby+h71AnnDqJWAl0kQqnyqQ1/fQzGwx4ZY767Ui4mZgJGnSgBeAQaSpv28B9vCXkS73KvAF4GrgMWAaKRiYQ5p58TJg+4g4olpgB76HvUVETAc+ABxDWndyFmn6+xdIXaS3qBbY5by+h73DXqTADgqubed7aGa2+HDLnZmZmZmZWR/gljszMzMzM7M+wMGdmZmZmZlZH+DgzszMzMzMrA9wcGdmZmZmZtYHOLgzMzMzMzPrAxzcmZmZmZmZ9QEO7szMzMzMzPoAB3dmZmZmZmZ9gIM7MzMzMzOzPsDBnZk1TdIISdHBxxllz8f09GtaHPn9NzMz6xsc3JmZmZmZmfUBDu7MrCOeB4bUeGxSlu7bddKd0431NTMzM+uzBvR0Bcxs0RURAbRWOydpZtnuOxFRNV0peadWzMzMzGwx5JY7MzMzMzOzPsDBnZn1uHoTekgam89NyvvrSvqZpEmSZkl6TtJFklYpyyNJh0n6q6Q3JM2Q9KCkQxqszwclXVV2jemSHpF0rqSV2sm7jaQrJT2T886S9KKkv0n6oaSd6uSVpI9I+o2k53PetyQ9msvcW9ISFXkGSNpR0vfya3xT0pyc7++Szil/b5qVr3OkpFsl/U/SO/m9nSDpOElLdfQaZmZm1jFKvarMzDqXpBHAc3n3mxFxRp20pX9Eh0fE2IpzY4HDSOP7jgRuII3Vq/QcsAMwGfg1sF+Ny50WEWfVqMfSwC+B/WvVFZgKfCQi/lwl/5eA86jfzfSfETGySt5VgN8C29XJCzAqIiaW5TseOL+dPG8B+0fE+Gon673/+fy6wM3ARnWu8W/gQxHxQjt1MTMzsy7iljszW1QMA34DPAXsCawIrA2cmc+vTZq45VzgQ8AZwIbA8sDOwOM53emSNqwsXFI/4EZSYPcO8H1gK2AFYFXgIOAJYChwY2UZktYHvksK7P4BHACsk6+/Sa7TT4DXq1x7GeBu2gK7XwI7ASsDw4Etga8A/6zyvrwN3AocBbwfeE+u8yakoPgfwLLA7yStWCV/XZJWBv5CCuxeA04ANgaWy9c6gRTwbgLcImlQ0WuYmZlZ53DLnZl1iS5ouYMUqGwfEbMq0lwJfAqYS/rR6qCI+G1FmvVJAV4/4NsR8fWK88cAFwHvArtFxN1V6jkU+BuwHnBTROxXdu4LwI9yHVaMiDdrvd4q5f4YODbvfjYifl4n7YCIeLdA2UsA40iB4+kRcWaVNPXe/98ABwIvAttGxH+r5B8J3A8sBXwxIs5vtH5mZmbWedxyZ2aLkpMqA7vs2rztD9xTGdgBRMRTwMN5d+sqZXwpb6+oFtjlMqYC38q7e0saVna6NPvwTGBKrRdQSdIQ4NN595Z6gV2uQ8OBXU4/B7gq7+5eJK+ktYCP5N2vVwvs8jUmAtfk3U8VuYaZmZl1Hi+FYGaLireBqmPGgGfLnv+pThnPAFsAC0wwIuk9pC6GAH+W1FKnjH/nbb9cVmnsXSlwHAL8QtI3IuKlOuWUbA8MzM8vbyD9QvJkJocC+wKbk7plVuseuUHBonel7UfA8e28L4/k7UhJS0bEOwWvZWZmZh3k4M7MFhWv51aoaspb8/5Xp4xSusrAp3z83NUF6jS89CQixkm6kTSRyxjgMEmPABPy484aXTXXLXv+cJXzdUlaB7iN1FW0PcMKFl/+vjQ6UUo/0ni8Vwpey8zMzDrI3TLNbFExtxPTVc5mOaxYVeYbWLF/EGnik2fyNd5LGkt3LfCKpF9VWZZgmbLn04tcXFJ/0uyh65EC1+8Cu5Aml1me1Io4BDgmZ+lfpHw6730xMzOzbuCWOzOzBYOqbSPi/mYKyS2L3we+n1vUtiUtz/BhYDXgYGB7Se/N4/cAppUVMYS0lEOjPkDqhgnwsYj4fbVEkpoNtkrvy2xgcHgGLjMzs17NLXdmZqmlrWSLzigwIp6NiKsi4ihgLeDsfGot4PCypE+XPR9Z8DKj8nZKrcAue2/BcktK78tA0vIHZmZm1os5uDOzxV5E/BsoTX5yhKR6i5A3U/5c0np8pZkuyxcD/yttYwHHFCx6qbztX6vOefmG/QqWW3J72fNP10xlZmZmvYKDOzOz5Ly8/T/gnPYSV1nEfL28EHota9HWFf6N0sGIaAVKyx/sI+mIdq5bPm6uNEvoENLMlpVp+wGXkBZeLywiniYt7A5wrKQPt1O3AXnmUTMzM+sBDu7MzJILgTvy869JukvSAZLWlDRU0uqSdpT0dUkTgV9V5D8FeEbSuZJ2y/mGSVpH0sHArTndXOA3FXlPBZ7Mzy+VNDZfa0VJy0saJel4Sf8ANivLdytta+pdJekzktaWNFzSbsCdwMdoW76hGccA/yUFpjdJulTSTpJWlrRsfn17Sfoh8DxwdAeuZWZmZh3gCVXMzEhdJyXtR2rpOgTYKT9qubfKsRHASflRzbvAsXnR7/JrT5e0M2nmy62Aw/KjvTpPk/Rp4NekZRl+ViXZr4C7gMvaK6/GNf4naQfgOtJ4xCPzo5a3m7mOmZmZdZxb7szMsoiYGRGfIgVYlwCPkWaznAtMBf6Zj3+IhQO/k0hB4eWk9epeIQVzrcC/gB8Bm0bEJTWu/TJpds1PADcDLwPvAG8CjwJXkGbdfLQi3+9yvt8BrwNz8rX/RJpB81NAh2a5jIjnSO/JR0mtjs+TxgnOyde8l9SV9X0RcUpHrmVmZmbNk2e2NjMzMzMzW/S55c7MzMzMzKwPcHBnZmZmZmbWBzi4MzMzMzMz6wMc3JmZmZmZmfUBDu7MzMzMzMz6AAd3ZmZmZmZmfYCDOzMzMzMzsz7AwZ2ZmZmZmVkf4ODOzMzMzMysD3BwZ2ZmZmZm1gc4uDMzMzMzM+sD/h+Nl3SHLwOeGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gt_tau = 100 # ground-truth timescale\n",
    "\n",
    "# fit \n",
    "popt, poptcov = abcTau.preprocessing.fit_oneTauExponential(data_ac, binSize, maxTimeLag)\n",
    "tau = popt[1]\n",
    "\n",
    "# check if estimated timescales with exponential fit are biased or not\n",
    "theta = np.array([tau])\n",
    "numTimescales = 1\n",
    "taus_bs, taus_bs_corr, err = abcTau.preprocessing.check_expEstimates(theta, deltaT, binSize, T, numTrials,\\\n",
    "                                                                     data_mean, data_var, maxTimeLag, numTimescales,\\\n",
    "                                                                     numIter = 100, plot_it = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result inidicates that the direct exponential fitting gives biased estimates with the current amount of data. We can use the aABC algorithm to get a more accurate estimation of timescales. aABC additionally provides the posterior distribution of timescales that we can use for model comparison (see Tutorial 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Estimation of timescales using the aABC algorithm\n",
    "\n",
    "In the following, we explain how to customize your own Python script step-by-step such that it best fits your data and your desired generative model. We provided 3 different examples:\n",
    "1) Fitting one timescale in time domain\n",
    "2) Fitting one timescale in frequency domain\n",
    "3) Fitting two timescales in time domain\n",
    "\n",
    "The duration for fitting the generative model with the aABC algorithm depends on the size of your data, selected summary statistics (e.g., fastest is with the PSD), number of accepted samples in posteriors (min_samples), final selected accepatance rate (minAccRate), and computational resources (e.g., number of cores in parallel processing). Depending on these parameters the fitting duration is in between ~5-20 min (usually for fitting PSD as summary statistic and parallel processing, see Supplementary Table 1 in the paper for more details) up to ~4-8 hours (wihtout any parallel processing).\n",
    "\n",
    "\n",
    "\n",
    "### [1]  Fitting a general type of generative model in the time domain:\n",
    "This example takes about 2 hours to run on a normal computer without any parallel processing, using parallel computing (parallel = True) can speed up the fitting at least to half of this duration. You can check the provided example Python scripts for parallel processing.\n",
    "\n",
    "The following is the example in Fig. 3a.\n",
    "\n",
    "1) import the package and set the parameters for parallel processing (if applicable):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add the path to the abcTau package\n",
    "import sys\n",
    "sys.path.append('./abcTau')\n",
    "# import the package\n",
    "import abcTau \n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# stetting the number of cores for each numpy computation in multiprocessing\n",
    "# import os\n",
    "# os.environ[\"OMP_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"OPENBLAS_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"MKL_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"NUMEXPR_NUM_THREADS\"] = \"2\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Define the directories and filenames for loading the data and saving the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# path for loading and saving data\n",
    "datasave_path = 'example_abc_results/'\n",
    "dataload_path = 'example_data/'\n",
    "\n",
    "# path and filename to save the intermediate results after running each step\n",
    "inter_save_direc = 'example_abc_results/'\n",
    "inter_filename = 'abc_intermediate_results'\n",
    "\n",
    "# define filename for loading and saving the results\n",
    "filename = 'OU_tau20_mean0_var1_rawData'\n",
    "filenameSave = filename"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Load the real data and extract required statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data time-series as a numpy array (numTrials * time-points)\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "# select summary statistics metric\n",
    "summStat_metric = 'comp_ac_fft'\n",
    "ifNorm = True # if normalize the autocorrelation or PSD\n",
    "\n",
    "# extract statistics from real data\n",
    "deltaT = 1 # temporal resolution of data.\n",
    "binSize = 1 #  bin-size for binning data and computing the autocorrelation.\n",
    "disp = None # put the dispersion parameter if computed with grid-search\n",
    "maxTimeLag = 50 # only used when using autocorrelation for summary statistics\n",
    "data_sumStat, data_mean, data_var, T, numTrials =  abcTau.preprocessing.extract_stats(data_load, deltaT, binSize,\\\n",
    "                                                                                  summStat_metric, ifNorm, maxTimeLag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Select the desired generative model from the list of 'generative_models.py' and the distance function from 'diatance_functions.py'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select generative model and distance function\n",
    "generativeModel = 'oneTauOU'\n",
    "distFunc = 'linear_distance'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Define the prior distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a uniform prior distribution over the given range\n",
    "# for a uniform prior: stats.uniform(loc=x_min,scale=x_max-x_min)\n",
    "t_min = 0.0 # first timescale\n",
    "t_max = 100.0\n",
    "priorDist = [stats.uniform(loc= t_min, scale = t_max - t_min)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6) Set the aABC fitting parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set fitting params\n",
    "epsilon_0 = 1  # initial error threshold\n",
    "min_samples = 100 # min samples from the posterior\n",
    "steps = 60 # max number of iterations\n",
    "minAccRate = 0.01 # minimum acceptance rate to stop the iterations\n",
    "parallel = False # if parallel processing\n",
    "n_procs = 1 # number of processor for parallel processing (set to 1 if there is no parallel processing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7) Create the model object: Just copy paste the following (this a general definition of the model object, but all parts including generative models, summary statistics computation or distance function can be replaced by your own functions. You can add your handmade functions inside respected modules: \"generative_models.py\", \"distance_functions.py\", \"summary_stats.py\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating model object\n",
    "class MyModel(abcTau.Model):\n",
    "\n",
    "    #This method initializes the model object.  \n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    # draw samples from the prior. \n",
    "    def draw_theta(self):\n",
    "        theta = []\n",
    "        for p in self.prior:\n",
    "            theta.append(p.rvs())\n",
    "        return theta\n",
    "\n",
    "    # Choose the generative model (from generative_models)\n",
    "    # Choose autocorrelation computation method (from basic_functions)\n",
    "    def generate_data(self, theta):\n",
    "        # generate synthetic data\n",
    "        if disp == None:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var)')\n",
    "        else:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var, disp)')\n",
    "               \n",
    "        # compute the summary statistics\n",
    "        syn_sumStat = abcTau.summary_stats.comp_sumStat(syn_data, summStat_metric, ifNorm, deltaT, binSize, T,\\\n",
    "                                          numBinData, maxTimeLag)   \n",
    "        return syn_sumStat\n",
    "\n",
    "    # Computes the summary statistics\n",
    "    def summary_stats(self, data):\n",
    "        sum_stat = data\n",
    "        return sum_stat\n",
    "\n",
    "    # Choose the method for computing distance (from basic_functions)\n",
    "    def distance_function(self, data, synth_data):\n",
    "        if np.nansum(synth_data) <= 0: # in case of all nans return large d to reject the sample\n",
    "            d = 10**4\n",
    "        else:\n",
    "            d = eval('abcTau.distance_functions.' +distFunc + '(data, synth_data)')        \n",
    "        return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "8) Run the aABC algorithm and save the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# fit with aABC algorithm for any generative model\n",
    "abc_results, final_step = abcTau.fit.fit_withABC(MyModel, data_sumStat, priorDist, inter_save_direc, inter_filename,\\\n",
    "                                                 datasave_path,filenameSave, epsilon_0, min_samples, \\\n",
    "                                                 steps, minAccRate, parallel, n_procs, disp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Done :) \n",
    "\n",
    "Just copy paste the content of all the above cells in a text file and rename it to \"scriptName.py\"\n",
    "\n",
    "Then, run it from terminal as: python scriptName.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [2]  Fitting a general type of generative model in the frequency domain:\n",
    "\n",
    "Depending on your application you can also use the aABC method to extract timescales and other parameters from the frequency domain by fitting the sample power spectrum density (PSD).\n",
    "\n",
    "For this purpose, we use exactly the same script as for the time domain, and only replace the autocorrelation computation with PSD computation.\n",
    "\n",
    "### Attention! You should use exactly the same method for computing the PSD of real data and synthetic data!\n",
    "\n",
    "We have implemented a simple PSD estimation function in the abcTau package but depending on your application, you can replace it with any other method such as methods avaiable in Scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html\n",
    "\n",
    "However, the result of aABC algorithm generally does not depend on the choise of PSD computation method, as long as the same method is used for both real and synthetic data.\n",
    "\n",
    "This example takes about 1 hour and 20 min to run on a normal computer without any parallel processing, using parallel computing can speed up the fitting to a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# add the path to the abcTau package\n",
    "import sys\n",
    "sys.path.append('./abcTau')\n",
    "# import the package\n",
    "import abcTau \n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# stetting the number of cores for each numpy computation in multiprocessing\n",
    "# import os\n",
    "# os.environ[\"OMP_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"OPENBLAS_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"MKL_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"NUMEXPR_NUM_THREADS\"] = \"2\" \n",
    "    \n",
    "    \n",
    "# path for loading and saving data\n",
    "datasave_path = 'example_abc_results/'\n",
    "dataload_path = 'example_data/'\n",
    "\n",
    "# path and filename to save the intermediate results after running each step\n",
    "inter_save_direc = 'example_abc_results/'\n",
    "inter_filename = 'abc_intermediate_results_psd'\n",
    "\n",
    "# load real data and define filenameSave\n",
    "filename = 'OU_tau20_mean0_var1_rawData'\n",
    "filenameSave = filename\n",
    "print(filename)\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "# load data time-series as a numpy array (numTrials * time-points)\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "# select summary statistics metric\n",
    "summStat_metric = 'comp_psd'\n",
    "ifNorm = True # if normalize the autocorrelation or PSD\n",
    "\n",
    "# extract statistics from real data\n",
    "deltaT = 1 # temporal resolution of data.\n",
    "binSize = 1 #  bin-size for binning data and computing the autocorrelation.\n",
    "disp = None # put the dispersion parameter if computed with grid-search\n",
    "maxTimeLag = None # only used when suing autocorrelation for summary statistics\n",
    "data_sumStat, data_mean, data_var, T, numTrials =  abcTau.preprocessing.extract_stats(data_load, deltaT, binSize,\\\n",
    "                                                                                  summStat_metric, ifNorm, maxTimeLag)\n",
    "\n",
    "\n",
    "\n",
    "# Define the prior distribution\n",
    "# for a uniform prior: stats.uniform(loc=x_min,scale=x_max-x_min)\n",
    "t_min = 0.0 # first timescale\n",
    "t_max = 100.0\n",
    "priorDist = [stats.uniform(loc= t_min, scale = t_max - t_min)]\n",
    "\n",
    "\n",
    "# select generative model and distance function\n",
    "generativeModel = 'oneTauOU'\n",
    "distFunc = 'logarithmic_distance'\n",
    "\n",
    "\n",
    "# set fitting params\n",
    "epsilon_0 = 1  # initial error threshold\n",
    "min_samples = 100 # min samples from the posterior\n",
    "steps = 60 # max number of iterations\n",
    "minAccRate = 0.01 # minimum acceptance rate to stop the iterations\n",
    "parallel = False # if parallel processing\n",
    "n_procs = 1 # number of processor for parallel processing (set to 1 if there is no parallel processing)\n",
    "\n",
    "\n",
    "# creating model object\n",
    "class MyModel(abcTau.Model):\n",
    "\n",
    "    #This method initializes the model object.  \n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    # draw samples from the prior. \n",
    "    def draw_theta(self):\n",
    "        theta = []\n",
    "        for p in self.prior:\n",
    "            theta.append(p.rvs())\n",
    "        return theta\n",
    "\n",
    "    # Choose the generative model (from generative_models)\n",
    "    # Choose autocorrelation computation method (from basic_functions)\n",
    "    def generate_data(self, theta):\n",
    "        # generate synthetic data\n",
    "        if disp == None:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var)')\n",
    "        else:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var, disp)')\n",
    "               \n",
    "        # compute the summary statistics\n",
    "        syn_sumStat = abcTau.summary_stats.comp_sumStat(syn_data, summStat_metric, ifNorm, deltaT, binSize, T,\\\n",
    "                                          numBinData, maxTimeLag)   \n",
    "        return syn_sumStat\n",
    "\n",
    "    # Computes the summary statistics\n",
    "    def summary_stats(self, data):\n",
    "        sum_stat = data\n",
    "        return sum_stat\n",
    "\n",
    "    # Choose the method for computing distance (from basic_functions)\n",
    "    def distance_function(self, data, synth_data):\n",
    "        if np.nansum(synth_data) <= 0: # in case of all nans return large d to reject the sample\n",
    "            d = 10**4\n",
    "        else:\n",
    "            d = eval('abcTau.distance_functions.' +distFunc + '(data, synth_data)')        \n",
    "        return d\n",
    "\n",
    "\n",
    "# fit with aABC algorithm for any generative model\n",
    "abc_results, final_step = abcTau.fit.fit_withABC(MyModel, data_sumStat, priorDist, inter_save_direc, inter_filename,\\\n",
    "                                                 datasave_path,filenameSave, epsilon_0, min_samples, \\\n",
    "                                                 steps, minAccRate, parallel, n_procs, disp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [3] Fitting the specific case of a two-timescale generative model in the time domain:\n",
    "\n",
    "We developed a bit different function for fitting the two-timescale generative models to help them converge faster. While you always can use the above script for all models, the following options might help to speed up the fitting when your model have two timescales.\n",
    "\n",
    "The following is the example in Fig. 3c.\n",
    "\n",
    "1) define directories and prepare the model object same as above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inhomPois_tau5_80_coeff04_T1000_trials500_deltaT1_data_mean1_data_var1.25\n"
     ]
    }
   ],
   "source": [
    "# add the path to the abcTau package\n",
    "import sys\n",
    "sys.path.append('./abcTau')\n",
    "# import the package\n",
    "import abcTau \n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# stetting the number of cores for each numpy computation in multiprocessing\n",
    "# import os\n",
    "# os.environ[\"OMP_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"OPENBLAS_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"MKL_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"NUMEXPR_NUM_THREADS\"] = \"2\" \n",
    "    \n",
    "    \n",
    "# path for loading and saving data\n",
    "datasave_path = 'example_abc_results/'\n",
    "dataload_path = 'example_data/'\n",
    "\n",
    "# path and filename to save the intermediate results after running each step\n",
    "inter_save_direc = 'example_abc_results/'\n",
    "inter_filename = 'abc_intermediate_results'\n",
    "\n",
    "\n",
    "# Define the prior distribution\n",
    "# for a uniform prior: stats.uniform(loc=x_min,scale=x_max-x_min)\n",
    "t1_min = 0.0 # first timescale\n",
    "t1_max = 60.0\n",
    "t2_min = 20.0 # second timescale\n",
    "t2_max = 140.0\n",
    "coef_min = 0.0 # coeffiecient or weight of the first timescale\n",
    "coef_max = 1.0\n",
    "priorDist = [stats.uniform(loc= t1_min, scale = t1_max - t1_min),\\\n",
    "             stats.uniform(loc= t2_min, scale = t2_max - t2_min),\\\n",
    "             stats.uniform(loc= coef_min, scale = coef_max - coef_min)] \n",
    "\n",
    "\n",
    "# select generative model and distance function\n",
    "generativeModel = 'twoTauOU_poissonSpikes'\n",
    "distFunc = 'linear_distance'\n",
    "\n",
    "\n",
    "# load real data and define filenameSave\n",
    "filename = 'inhomPois_tau5_80_coeff04_T1000_trials500_deltaT1_data_mean1_data_var1.25'\n",
    "filenameSave = filename\n",
    "print(filename)\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "# select summary statistics metric\n",
    "summStat_metric = 'comp_ac_fft'\n",
    "ifNorm = True # if normalize the autocorrelation or PSD\n",
    "\n",
    "# extract statistics from real data\n",
    "deltaT = 1 # temporal resolution of data.\n",
    "binSize = 1 #  bin-size for binning data and computing the autocorrelation.\n",
    "disp = None # put the dispersion parameter if computed with grid-search\n",
    "maxTimeLag = 110 # only used when suing autocorrelation for summary statistics\n",
    "data_sumStat, data_mean, data_var, T, numTrials =  abcTau.preprocessing.extract_stats(data_load, deltaT, binSize,\\\n",
    "                                                                                  summStat_metric, ifNorm, maxTimeLag)\n",
    "\n",
    "\n",
    "\n",
    "# set fitting params\n",
    "epsilon_0 = 1  # initial error threshold\n",
    "min_samples = 100 # min samples from the posterior\n",
    "steps = 60 # max number of iterations\n",
    "minAccRate = 0.01 # minimum acceptance rate to stop the iterations\n",
    "parallel = False # if parallel processing\n",
    "n_procs = 1 # number of processor for parallel processing (set to 1 if there is no parallel processing)\n",
    "\n",
    "\n",
    "# creating model object\n",
    "class MyModel(abcTau.Model):\n",
    "\n",
    "    #This method initializes the model object.  \n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    # draw samples from the prior. \n",
    "    def draw_theta(self):\n",
    "        theta = []\n",
    "        for p in self.prior:\n",
    "            theta.append(p.rvs())\n",
    "        return theta\n",
    "\n",
    "    # Choose the generative model (from generative_models)\n",
    "    # Choose autocorrelation computation method (from basic_functions)\n",
    "    def generate_data(self, theta):\n",
    "        # generate synthetic data\n",
    "        if disp == None:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var)')\n",
    "        else:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var, disp)')\n",
    "               \n",
    "        # compute the summary statistics\n",
    "        syn_sumStat = abcTau.summary_stats.comp_sumStat(syn_data, summStat_metric, ifNorm, deltaT, binSize, T,\\\n",
    "                                          numBinData, maxTimeLag)   \n",
    "        return syn_sumStat\n",
    "\n",
    "    # Computes the summary statistics\n",
    "    def summary_stats(self, data):\n",
    "        sum_stat = data\n",
    "        return sum_stat\n",
    "\n",
    "    # Choose the method for computing distance (from basic_functions)\n",
    "    def distance_function(self, data, synth_data):\n",
    "        if np.nansum(synth_data) <= 0: # in case of all nans return large d to reject the sample\n",
    "            d = 10**4\n",
    "        else:\n",
    "            d = eval('abcTau.distance_functions.' +distFunc + '(data, synth_data)')        \n",
    "        return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) But instead use the following function for fitting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# fit with aABC algorithm for the two-timescale generative model\n",
    "abc_results, final_step = abcTau.fit.fit_withABC_2Tau(MyModel, data_sumStat, priorDist, inter_save_direc,\\\n",
    "                                                      inter_filename,\\\n",
    "                                                 datasave_path,filenameSave, epsilon_0, min_samples, \\\n",
    "                                                 steps, minAccRate, parallel, n_procs, disp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resuming a previous aABC fit and let it run for more iterations \n",
    "\n",
    "It's possible to resume a previous fit with the aABC method and allow it to run from the point it stopped. This functionality is useful when the initial stopping criterion (e.g., minimum acceptance rate) was not sufficient to obtain narrow posteriors and you need better statistics for the fits without running the aABC fit from the beginning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OU_tau20_T1000_trials500_deltaT1_data_mean0_data_var1\n"
     ]
    }
   ],
   "source": [
    "# add the path to the abcTau package\n",
    "import sys\n",
    "sys.path.append('./abcTau')\n",
    "# import the package\n",
    "import abcTau \n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "# stetting the number of cores for each numpy computation in multiprocessing\n",
    "# import os\n",
    "# os.environ[\"OMP_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"OPENBLAS_NUM_THREADS\"] = \"2\"\n",
    "# os.environ[\"MKL_NUM_THREADS\"] = \"2\" \n",
    "# os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"2\"\n",
    "# os.environ[\"NUMEXPR_NUM_THREADS\"] = \"2\"\n",
    "    \n",
    "    \n",
    "# path for loading and saving data\n",
    "datasave_path = 'example_abc_results/'\n",
    "dataload_path = 'example_data/'\n",
    "\n",
    "# path and filename to save the intermediate results after running each step\n",
    "inter_save_direc = 'example_abc_results/'\n",
    "inter_filename = 'abc_intermediate_results'\n",
    "\n",
    "# load real data and define filenameSave\n",
    "filename = 'OU_tau20_T1000_trials500_deltaT1_data_mean0_data_var1'\n",
    "filenameSave = filename\n",
    "print(filename)\n",
    "data_load = np.load(dataload_path + filename + '.npy')\n",
    "\n",
    "# select summary statistics metric\n",
    "summStat_metric = 'comp_ac_fft'\n",
    "ifNorm = True # if normalize the autocorrelation or PSD\n",
    "\n",
    "# extract statistics from real data\n",
    "deltaT = 1 # temporal resolution of data.\n",
    "binSize = 1 #  bin-size for binning data and computing the autocorrelation.\n",
    "disp = None # put the dispersion parameter if computed with grid-search\n",
    "maxTimeLag = 50 # only used when suing autocorrelation for summary statistics\n",
    "data_sumStat, data_mean, data_var, T, numTrials =  abcTau.preprocessing.extract_stats(data_load, deltaT, binSize,\\\n",
    "                                                                                  summStat_metric, ifNorm, maxTimeLag)\n",
    "\n",
    "\n",
    "# Define the prior distribution\n",
    "# for a uniform prior: stats.uniform(loc=x_min,scale=x_max-x_min)\n",
    "t_min = 0.0 # first timescale\n",
    "t_max = 100.0\n",
    "priorDist = [stats.uniform(loc= t_min, scale = t_max - t_min)]\n",
    "\n",
    "\n",
    "# select generative model and distance function\n",
    "generativeModel = 'oneTauOU'\n",
    "distFunc = 'linear_distance'\n",
    "\n",
    "\n",
    "# creating model object\n",
    "class MyModel(abcTau.Model):\n",
    "\n",
    "    #This method initializes the model object.  \n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    # draw samples from the prior. \n",
    "    def draw_theta(self):\n",
    "        theta = []\n",
    "        for p in self.prior:\n",
    "            theta.append(p.rvs())\n",
    "        return theta\n",
    "\n",
    "    # Choose the generative model (from generative_models)\n",
    "    # Choose autocorrelation computation method (from basic_functions)\n",
    "    def generate_data(self, theta):\n",
    "        # generate synthetic data\n",
    "        if disp == None:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var)')\n",
    "        else:\n",
    "            syn_data, numBinData =  eval('abcTau.generative_models.' + generativeModel + \\\n",
    "                                         '(theta, deltaT, binSize, T, numTrials, data_mean, data_var, disp)')\n",
    "               \n",
    "        # compute the summary statistics\n",
    "        syn_sumStat = abcTau.summary_stats.comp_sumStat(syn_data, summStat_metric, ifNorm, deltaT, binSize, T,\\\n",
    "                                          numBinData, maxTimeLag)   \n",
    "        return syn_sumStat\n",
    "\n",
    "    # Computes the summary statistics\n",
    "    def summary_stats(self, data):\n",
    "        sum_stat = data\n",
    "        return sum_stat\n",
    "\n",
    "    # Choose the method for computing distance (from basic_functions)\n",
    "    def distance_function(self, data, synth_data):\n",
    "        if np.nansum(synth_data) <= 0: # in case of all nans return large d to reject the sample\n",
    "            d = 10**4\n",
    "        else:\n",
    "            d = eval('abcTau.distance_functions.' +distFunc + '(data, synth_data)')        \n",
    "        return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the previous aABC fit\n",
    "data_abc_path = 'example_abc_results/'\n",
    "filename = '1on_OU_tau20_T1000_lag50_steps43'\n",
    "abc_results = np.load(data_abc_path + filename + '.npy', allow_pickle=True)\n",
    "ind = filename.find('steps') \n",
    "final_step = int(filename[ind+5] + filename[ind+6])\n",
    "abc_results_old = abc_results[:final_step]\n",
    "\n",
    "# set new fitting params\n",
    "epsilon_0 = 1  # initial error threshold\n",
    "min_samples = 100 # min samples from the posterior\n",
    "steps = 60 # max number of iterations\n",
    "minAccRate = 0.005 # new minimum acceptance rate to stop the iterations\n",
    "parallel = False # if parallel processing\n",
    "n_procs = 1 # number of processor for parallel processing (set to 1 if there is no parallel processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# resume the aABC fitting\n",
    "abc_results, final_step = abcTau.fit.fit_withABC(MyModel, data_sumStat, priorDist, inter_save_direc,\\\n",
    "                                                             inter_filename,\\\n",
    "                                                 datasave_path,filenameSave, epsilon_0, min_samples, \\\n",
    "                                                 steps, minAccRate, parallel, n_procs, disp, resume = abc_results_old)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}