{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains all the code to interpret and plot the data obtained from the experiment. The cells depend on eachother so it's best to run them in order to load all the variables. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matchgate import MatchgateSequence\n",
    "from orth_matchgate import OrthMatchgateSequence\n",
    "from esprit import esprit\n",
    "import numpy as np\n",
    "import pickle as pk\n",
    "from scipy.optimize import curve_fit\n",
    "from toolkit import *\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import f\n",
    "from sklearn.utils import resample\n",
    "from pandas import DataFrame\n",
    "from scipy.special import binom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data from experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"qpu_q31q32_oct24_400shots_ORTHO_4pi.pkl\"\n",
    "data = DataFrame()\n",
    "with open(data_path, \"rb\") as fl:\n",
    "    data = pk.load(fl)\n",
    "\n",
    "gate_path = \"qpu_q31q32_oct20_randomsequences_ORTHO.pkl\"\n",
    "\n",
    "with open(gate_path, \"rb\") as fl:\n",
    "    sequences = pk.load(fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4  5  6  7  8  9 10 12 14 16 18 20 22 24]\n",
      "[64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64]\n"
     ]
    }
   ],
   "source": [
    "## Check what the different sequence lengths are\n",
    "\n",
    "#We discard the first pair of sequence lengths. \n",
    "#This is good practice in RB experiments to avoid any gate-dependent noise related outliers.\n",
    "\n",
    "sequence_lengths = np.array(list(sequences.keys()))[2:]\n",
    "print(sequence_lengths)\n",
    "\n",
    "##check how many ramdom sequences per sequence length were sampled\n",
    "print([len(sequences[x]) for x in sequence_lengths])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Recompute correlators\n",
    "#Optional! (This takes a while)\n",
    "# for m in sequences.keys():\n",
    "#     for seq in sequences[m]:\n",
    "#         seq.compute_correlator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Normalization factors for each decay\n",
    "normalisation = [1, 4, 4 / 3, 4, 1]\n",
    "\n",
    "corr_data = {}\n",
    "average = {}\n",
    "variance = {}\n",
    "\n",
    "\n",
    "for k in range(5):\n",
    "    if k%2 == 0:\n",
    "        basis = \"Z\"\n",
    "    else:\n",
    "        basis = \"X\"\n",
    "\n",
    "    corr_data[k] = np.real(\n",
    "        [\n",
    "            normalisation[k]\n",
    "            * np.array(\n",
    "                [\n",
    "                    sum(sequence.correlators[k] * data[m][sequence.sequence_id][basis])\n",
    "                    for sequence in sequences[m]\n",
    "                ]\n",
    "            )\n",
    "            for m in sequence_lengths\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    average[k] = np.array(\n",
    "        [np.mean(corr_data[k][m]) for m in range(len(sequence_lengths))]\n",
    "    )\n",
    "\n",
    "    variance[k] = np.array(\n",
    "        [np.var(corr_data[k][m]) for m in range(len(sequence_lengths))]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit data to exponential\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Stores the extracted Majorana fidelities\n",
    "poles = {}\n",
    "\n",
    "#Stores the associated amplitudes\n",
    "amplitudes = {}\n",
    "\n",
    "#Starting points for curve_fit for each curve\n",
    "inits = [[1, 1], [0.9, 1], [0.9, 1], [0.9, 1],[0.9, 1]]\n",
    "\n",
    "#Fit an exponential to each curve, extracting pole and amplitude.\n",
    "for k in range(5):\n",
    "    for i in [0]:\n",
    "        popt, pcov = curve_fit(\n",
    "            exp, sequence_lengths, average[k], p0=inits1[k], maxfev=10000\n",
    "        )\n",
    "        poles[(i, k)] = popt[: i + 1]\n",
    "        amplitudes[(i, k)] = popt[i + 1 :]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEOCAYAAACXX1DeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABmPElEQVR4nO2deXhdVbm432/vM5/MTZqm6TzRgbYUwjxJQWSSyZ8KgopyRRyu81UcQBS9KiJyvUwioIgiDuAVZZ4pk3SgtNDS0rlJ0zRpmvnM+/v9sU/aJM10kpMmTdf7POc55+y99lrfPsP+9lrfJKqKwWAwGAw9YQ23AAaDwWAY2RhFYTAYDIZeMYrCYDAYDL1iFIXBYDAYesUoCoPBYDD0ime4Bcg2xcXFOmXKlOEWw2AwGA4qli9fXqeqJd3tG3WKYsqUKSxbtmy4xTAYDIaDChHZ2tM+s/RkMBgMhl4xisJgMBgMvWIUhcFgMBh6xSgKg8FgMPSKURQGg8Fg6BWjKAwGg8HQK0ZRGAwGg6FXjKIwGAwGQ6+MOkURc8DU2DAYDIbsMeoUxYZW2BoZbikMBoNh9DDqFAXAy7uHWwKDwWAYPYw6ReG3jKIwGAyGbDLqkgLmeGBpA0RSStCW4RYnYzTZRqpqCU5rNVa4DLv8ZMQTylr7AzHGSJRpIJgxRkb/ZozhR4bL8Csi9wLnAbtU9fBu9gvwP8A5QBtwhaqu6Kvf2UdUqPfmpdy2AE4pPrgURWr3O8T+/QNQB1IxsP0gFv5jv489Zt6g2x+IMUaiTAPBjDEy+jdjHDhEZLmqVnS7bxgVxSlAC/D7HhTFOcB/4iqKY4H/UdVj++r3qIoK/fnjSzmuEEKeg0dRaLKNyFNXQLIbS7wnSPDM+xBPcMDtD8QYI1GmgWDGGBn9mzEOLCNSUQCIyBTgXz0oil8DL6jqn9Lv1wHvU9Xq3vo8anaZvnLPFZ222eNPwjv1XDQZdTV6FzwTT8cz6Qw01khs2U/33z/lHDzlJ+NEaomvuHn//dMvwjPuGJyWSuJv3bbffu+sj2KXHIHTuIn427/Zf/+cT6DN24iv/jU4if1PyvLim381EhpLYv2fAdBoPdpaDXT3/QkSHo8ECgHwHfk1rGAJsZW3ktr2VPfH2H7s8SejbTv3burvGIETfwJAbOlPSFW/1nP/JUegidZOmzXZhrZUuXdY3YxhFcwicMpNAMTf/g2pmmU9yyQWEirbe94AVv40fId/xpVv+S/QaF2nQ6zC2fjmfjIt/3+j8eaMPlu79Gi8My4GIPrKt/c/7R5+e/0dYzC/PY3Wo5Ha7n9T6TF8i76MXTSHVP1aEmt/v18r3+GfwcqfRqp25d7fXib9+0+4AStYQrJqCcktj+3Xyl9xDeLPJ7ntGZLbn91vvz3+RBJrftfj78O38It4Jp9JYsPDpGqWdjnYR+A49zNPrHuQVN1bnY/25eI/+jsktz5JfNUdoKn9h7D9+A6/Cqd5K07jps7H55TjX/hFAGJv3er+jjvQ8bcXfenrOA3v0dN/Q3Im7Kco7OKFeA+7pJvzHhp6UxQj2ZhdDmzv8L4yvW0/ROQqEVkmIsviiQR1msM9ydPYoQUHQs6s4LRW9/CHA5yEu78DmorT/UUGQNFu/lga29PzMakYGmsY1Bgab+69/3jz/tuT0R4uAukxutyB9SqTOt2ed6YM5LMdsWP09JvKwhhD3T+A07qz199H1//FwMao7l5JAKRiWRnD/R33/N8gGR30GEPJSJ5RPAr8RFVfTr9/Fvimqi7vrc+Kigr9+5KlnPMaXDMTLpt4cCw/Jbc+6c42uvtTpO9qPJPPHHD7AzHGSJRpIJgxRkb/ZowDy8E6o6gEJnZ4PwHY0Z8DJwaFKSFYchC5ydrlJ4P08HWI5e4fRPsDMcZIlGkgmDFGRv9mjJHDSFYUjwCfEJfjgMa+7BMdOWnMPjfZgwHxhPAf+33wBF1vCHCfPUH8x35/v/XLTNsfiDFGokwDwYwxMvo3Y4wchtPr6U/A+4BioAb4PuAFUNU70+6xtwJn4brHfkpVl/XVb0VFhS5btozX6pWrVsL/LoD3HURuspqMdONn3fOPKNP2B2KMkSjTQDBjjIz+zRgHhhHr9TQUtCuKhKOc8Qp8ZgpcfpDYKQwGg2G46E1RjLrI7Ha8lvDMiYrXMkrCYDAYBsNItlEMmnYlMdpmTQaDwXAgGdWKIukoly9Tbts83JIYDAbDwUvGikJE/CIyVUTmikjJUAiVLTyW4LHg+bq+2xoMBoOhe/qlKEQkV0Q+JyIvAY3ABuBtYKeIbBeR34jI0UMp6EA5rRjWt8COiFl+MhgMhoHQp6IQka8CW4BPA08DFwBHALOA43HdWj3A0yLyhIjMHCphB8Kpxe7zCwdR8J3BYDCMJPrj9XQCcKqqvt3D/jeAe0XkauBK4FTgvSzJN2imhIQpIeWFWvjYhOGWxmAwGA4++lQUqvrh/nSkqjHg9kFLNFhSrpeTG6/n8ulJ4AyjSAaDwXAwM+riKJwGh9TWJJ4p3r3bLhpvYikMBoNhoGSkKERkAvA1oAzYDLwJvKmqG4ZAtoFhQ2x5DHuiB+lQCnV3XFnXAicUGaVhMBgMmZCpe+xDuMbsCHAM7lLTehFpTHtEDTsSFrRFSaztnCf/tk3wldUQPUiSBBoMBsNIIdOlp8OB41R1dfuG9CxjEbAgm4INFPEKdrlNfFUM7wwPEnB14Rkl8Ncd8Fo9nDaioz8MBoNhZJHpjGIZkNNxg6pWquo/VfXH2RNrcPgr/KCQqt1nwj66EHI98EztMApmMBgMByGZKoqvAzeISMEQyJI1rAKb8Idz8EzcN2HyWsJpxW6UdsIxy08Gg8HQXzJVFM1AGFgnIneIyCdEZL6I2EMg26AQn6CqpHal9iYFPL0EmpPwTjelmw0Gg8HQPZkqij8DBennMuAG4C2gWUTeyK5ogye1NUnk8TZSVW7h9BOL4KkT4Ih84/lkMBgM/SVTY/ZM4BhVfad9g4gUAkfipvUYUdgTPUieRWxpFLssjN8Wykbc3MdgMBhGNpnOKF4HijpuUNU9qvqsqv4ie2JlB7EF/zF+tGmfu2xVRPnSKmVlo7FTGAwGQ3/IVFH8GviBiBQPhTBDgafcgz3BdZd1Ig4FXni1Hh6vGW7JDAaD4eAg06WnB9PP60XkX7gzjDeBt1S1LauSZRF/RYDI021ok0O41MMpY5SndsE3Zyq2GHuFwXCwkEgkqKysJBqNDrcow4pt2xQUFFBcXIxlDX39uUwVxURcW8TC9POXgemAish7qjo3q9JlCSvfInRxGEmXRv1AKTxdC8sb4JjC4ZXNYDD0n8rKSnJzc5kyZUqnxJ+HEqpKIpGgpqaGyspKJk2aNORjZqQoVLUKqAIebd8mIiFcxTEiIrN7QixBHSW5KcnJk22CtvB4jVEUBsPBRDQaPaSVBICI4PP5KC8vZ926dQdkzD4VhYhMVdUeq06nl5xeA14T99uboKrbsyhj1kjtSBF7JYov4efyCV6KfMMtkcFgyJRDWUl05EAsOe0dqx9tXhORe0Tk+J4aiEihiHwOWIObNHBEYpfb2GU28TdjfHG8cvlE84MzGAyGvujP0tNs4LvAoyKSApYD1UAUKATmAnNwK919RVWfHCJZB42I4D82QNsjrcSXxZATAmxohXl5RmEYDAZDT/Q5o1DVBlX9L6Ac+BzwLm509lQgCdwHLFLVE0eykmjHyrfwHu4juSnJH5cm+cxKiJvcTwaDYZBMmTKFZ555ZrjFGBL6vcilqhFV/ZuqfkVVL1LVs1T1clX9RS/1tEckvvk+7IkejigWmpPw8u7hlshgMBzq1NfXc9FFFxEOh5k8eTIPPPDAcIu0l0FZQ0TkC9kS5EAiHiG4OMj8aTaFXnjMBN8ZDKMeTbaR3Pok8TW/I7n1STQ5skK/vvCFL+Dz+aipqeGPf/wjn/vc53jnnXf6PvAAMFiz+WIR+UD7GxHJEZGH+3uwiJwlIutEZIOIXNPN/nwR+aeIvCUi74jIpwYpbyc8SfhxW4z3qh2ak2b5yWAYraR2v0PkqSuIv/0bkhseIv72b4g8dQWp3UNzIX733XeZOnUqDz74YN+NgdbWVh566CFuuOEGcnJyOOmkkzj//PO5//77h0S+TBmsovgEcK2IzBOROcAS4O/9OTCdmvw24Gxcg/ilItI1YO8LwBpVXQi8D/iFiGTNqVWTyoK6BF/bHeOlWqMoDIbRiCbbiP37B5CMQCrmbkzFIBkh9u8foMlIVsdbsWIFZ555Jv/7v//LJZdcwnnnnUdBQUG3j/POOw+A9evXY9s2s2bN2tvPwoULR8yMItPIbABE5FfAyvTjSuBP6V2f6FgmtQ+OATao6qZ0nw/iutau6dBGgdx0fEYOUI9rQM8KVsgiWOHnqNdj+JqTUGYCKwyG0Uaqagmo0/1OdUhVLcEz+cysjLVkyRLuuece7r//fk477TQA/vWvf/V5XEtLC/n5+Z225efn09w8Morn9DmjSBcoCnXZ/CRQCnwT+CswBdgJnCsi5/Vz7HKgY2BeZXpbR27Fdb3dAawGvqy6/zcuIleJyDIRWVZbm1mtU+8sL1apTXx5DKethx+TwWA4aHFaq/fNJLqSirn7s8Sdd97JCSecsFdJ9JecnByampo6bWtqaiI3Nzdrsg2G/iw9XQnkAojIQyKSo6qPqupPVPUSVV2AW8ToWmAXcHo/x+4ueKHr+s8HcGct43FzS90qInn7HaR6l6pWqGpFSUlJP4dPCyGC93g/yQSsf6mHH5PBYDhoscJlYPu732n73f1Z4s4772Tbtm189atf3bvt7LPPJicnp9vH2WefDcCsWbNIJpO89957e4976623mDdvXtZkGwz9URTbgWPTry8Cus4uUNWYqi5X1XtV9atd9/dAJW6SwXYm4M4cOvIp4GF12QBsxg0AzCrefJv7Jvi5OeTbWzbVYDCMDuzyk0F6uNSJ5e7PErm5uTzxxBO89NJLXHON65/z+OOP09LS0u3j8ccfByAcDnPxxRdz3XXX0drayiuvvMI//vEPPv7xj2dNtsHQH0XxU+AhEVmBe8f/aRE5pbs7+wxZCswUkalpA/UlwCNd2mwjPUMRkVLgMGDTIMftlvGHeXktabGmSdGUURYGw2hBPCH8x34fPMF9MwvbD54g/mO/j3iCWR2voKCAp59+mscff5xrr72238fdfvvtRCIRxo4dy6WXXsodd9wxYmYU0p87aBE5DNfQ/FNgPTANsIGt7DNqr1TVrhf6vvo9B7gl3de9qvpjEbkaQFXvFJHxwO9wl7YE+Kmq/qG3PisqKnTZsmWZiAFAY0I5bYlyb2uU2YVC4KTs/ngMBsPgWbt2LXPmzBnQsZqMkKpagtNajRUuwy4/OetK4kAzmM+jKyKyXFUrutvXL68nVV0H3CgiVwInAc3A4bh2gyOA9wNfw03t0W9U9THgsS7b7uzwegeQHXeEPsj3CqeNhTc2WszYmCA5OYln4oCcwgwGwwhEPMGseTcdamRaj+KwDm+Xpx+jho+Mh2dtH7o+Rey1KPbYMOI3CQMNBsOhzYFLaH4QcGyR8J05FqGTAmhUiS09tMstGgwGAxhFsR+qyirbIjnHS6oqhUZNbIXBYDi0MYqiC9VR+MQK+GuRj9AFISRgPiKDwXBoY66CXRgfFI4thIdrBPW7dbYTWxImvsJgMByyGEXRDReWQWUUljVAclOS2ItRkpuylmLKYDAYDioyVhQicraI/EtE1ojIxPS2/xCR/qbuGPGcUQK5HnhoB3imebDG2sT+HcVpMfYKg8Fw6JGRohCRy4C/AO/hlkL1pnfZuAkCRwUBWzi3ND2jAAInBwCILomipmyqwWDoBlMKdR/fBD6TzufUcS3mddzAu1HDF6fB48eD1xKsHAv/sQGcXSkSb8eHWzSDwTAKufXWW6moqMDv93PFFVcMtzidyDT0eCbwWjfbW4DB5n4aUeR73UA7J23E9kzz4DT4sMebaG2D4WDiuBeV1lTP+8M2vH7q8AfWjh8/nu9973s8+eSTRCLZLaY0WDKdUewAZnWz/RRg4+DFGVlsalXOeQ1er3fTkfuP8mMX2wBmCcpgOEjoTUn0Z/9AyLQUKsDFF1/MhRdeyJgxY7Iv0CDJVFHcBfxKRE5Mv58oIp8EbgTuyKpkI4AJQYik4C8dkp+rKrE3osReiRqXWYPBsB8DKYU60sk019ONIpIPPA0EgOeBGHCTqt42BPINKz5LuLBMuW871MSUUr8gIohfSKxNYJcl8c7w9t2RwWA4JBhoKdSRTsbusar6XaAYt+b1cUCJqvY/6foQE1Nnr10hG/y/cnAUHu4wq/DO92GPS7vMNhqXWYPB4DLQUqgjnQEF3Klqm6ouU9U3VLUl20INhm3xKI831WWtv4lB4YQi+FsVJNJ2CbEE/8kBsIXoixE0aZagDAbDwEuhjnQyWnoSkXt72KVAFNgA/DldR2JYCFk29+7ewXHhAsZ4srMs9LmpsDsOVgfHCCtkETgpQPTFCE69gz3WzspYBoPh4KW9FOrpp5/ONddcw09/+tO95U77IplMkkwmSaVSpFIpotEoHo8Hj2f4PS0znVGUABcDFwIz0o8L09sOw42zWCciR2RNwgwp8fiIq8NddZVZ63NhvrC4RLClswudZ4KH8IfCRkkYDIa9DLQU6o9+9COCwSA//elP+cMf/kAwGORHP/rREErafzJVVa/gxkxcqaptACISAn4DvAWcA/we+AXpWtcHGp8IlxSO4/76as5oLeLocH5W+m1OKr/bBouLYV7ePoUhAQtVJfleAnusjVUweKXhRBySm5N4pnqwgiYdl8EwGMJ27y6w4Szd523ZsmXv66KiIt56662Mjr/++uu5/vrrsyNMlslUUXwZWNyuJMC1V4jIj4Fn015RPwOGNY79w4WlvNbaQF0ykbU+LeCB7bA9Ajd2rXceh/ibccQvBM8NId6BBe9oSom/Eyexyo3+jq+I4V3gwzfPh9jDHxBkMByMjIRguoOdTG9Xc4CybraPS+8DaCJzBZRVfGLxPxNmc3Z+cdb6DHuEi8bD07tcV9mOiF/wnxLAaXKIvTqw+ApNKW1/byWxOg4p9j4Sq+K0/b0VTRmDucFgGB4yVRR/B+4RkQ+LyBQRmSwiHwbuAR5OtzkGWJ9NIQeCLYKq8lLLHtZGW7PS56UTIKXw527MH54yD75FPpJbkiTWDmAmkwKNaOcMWh23D0H0qMFgMPSHTBXF1cCTwB9wU3ZsSr9+Avh8us1a4DPZEnAwxFW5u66Km2u2EHcGH+8wMSi8rxj+tgMi3dzhew/3YU/0EF8WMynJDQbDqCEjRZGOn7gaKAIWAUcCRar6OVVtTbdZqaorsy7pAPBbFl8eO4ntiRgP7KnOSp+fnAQVBdDSTR0jESFwUoDA4iBWjjFCGwyG0cGAbAlppbAqy7IMCUeF8nh/bhF/2VPDieFCZgZCg+uvQDiqoOf94hM8E9yPNVWbwiqy+meItkGCgka7LDN5XBsIxgPXYDAMExkrChHx4NohJgG+jvtU9fdZkiurXFU8gRVtzfx81xZumzgbrwz+bn9zq9KSgvl53SsBp8Uh8kQbnike/CcFEOldWYgthC4Kd/J6AvAt8OGda7yeDAbD8JFpZPZs4J+41e0E997XAyRwkwOOSEWRa3v4eulkmlPJrCgJVeVLqyFkw4MV2q0SsHIsfAt8xFfGsYoS+Ob5uumpM2IL/gV+vDO9Jo7CYDCMGDK9Ct0CLAfygTZgDlABrAQ+lE3Bss2RoTxOzS0CIDXIpIEiwicmwppmWNrQczvvAh/2ZA/x5TGSVd0YNXrAClr45vqMkjAYDiJMKdR9HA38KG2jcACPqq7ATd3xi0wHF5GzRGSdiGwQkWt6aPM+EVkpIu+IyIuZjtGVp5t284Xta4k6g/M3PX8cFHnh3q09txERAicGsAosNyeU8YQyGAzdEIvFuPLKK5k8eTK5ubksWrSo3zmiDgSZKgrBnUkA1ALl6deVuHmf+t+RiA3cBpwNzAUuFZG5XdoUALcD56vqPODDGcq7H6VeH1vjUe6qqxpUP35buHwivFIP7zb3PEMRrxBYHMS3wIeEjZ3BYBhOHFWeatqd1VIE2SCZTDJx4kRefPFFGhsbueGGG/jIRz7SKS3IcJKpongbWJh+/QbwLRE5FfgBbubYTDgG2KCqm1Q1DjwIXNClzceAh1V1G4Cq7spwjP1YEMzl4oKxPNZUx2utDYPq66PlUOyDtc29t7NyLHyH+xERnBbHRFkbDMPEK60N3LxrK68O8r/fF5mWQg2Hw1x//fVMmTIFy7I477zzmDp1KsuXLx9SOftLporixx1efw+YiFvl7kzgSxn2VQ5s7/C+kn0zlHZmAYUi8oKILBeRT3TXkYhcJSLLRGRZbW1tnwN/csx4pvuC3FKzjfpB5IPK8wpPHA8Xje/fTEFjSuTRNmKvmTKqBsOBxkkH4ALcXVc1ZLOKbJRCrampYf369cyb1zWx3PDQb68nEfECPwQ+AaCqm4C5IlIE7NHMr3zdXV279uEBjsLNRBsEXhOR11W1U4oQVb0Lt543FRUVfcrhE4tvjZvCf25fx7K2Js7MG3gxc7/tpgrZ1AbT+1haEr/gOcxL4q04Vm4c30L/gMc1GAyZ8UprA40p16mkIZXk1dYGTsopzOoY2SiFmkgkuOyyy/jkJz/J7NmzsyrfQOn3jEJVE7husdple/0AlAS4M4iJHd5PALoWPKoEnlDVVlWtA15i39LXoJjkC/LbyfMGpSTa+XMVXPRvN7aiL3wLfXime4ivjJPYmL3stgaDoWfaZxNRdR1KouoMyaxisKVQHcfh4x//OD6fj1tvvTWrsg2GTJee7iN7eZyWAjNFZKqI+IBLgEe6tPkHcLKIeNJ1L47FzSWVFYrSFfBWR5rZEGvro3XPvH8s+C24uxcPqHZEBP/xAbfm9itRUrtMtj+DYajpOJtop31WkU0GUwpVVbnyyiupqanhoYcewuvNToXObJBpZHYYuExE3o8bT9EpLauq9ttOoapJEfkibpJBG7hXVd8RkavT++9U1bUi8gRuuhAHuFtV385Q5l5JqMPPa7biEeHWibMJWZnnyhjjEz5crjxQCZ+dokwK9R2FHTgtSPytGFaRiZUwGIaSrrOJdtpnFSeEC7D6yJzQXwZTCvVzn/sca9eu5ZlnniEYDGZFnmyR6VVqDrAC2ANMA+Z3eBye6eCq+piqzlLV6ar64/S2O1X1zg5tfq6qc1X1cFW9JdMx+sIrFt8sncLORIz/3bVtwEbmT00Cr8Adm/vXXnyC/+gA4hE0rjitJsbCYBgKXm1toDYZ73bfrmQ867OKgZRC3bp1K7/+9a9ZuXIl48aN2zvj+OMf/5hV2QZKRjMKVR3YwtsI5/BgDpcVlXF/fTULQ7mclZd5waMSv3DpBOWfO6E1qYQ9/fSEUiX6fASnzSF0dggJmBmGwZBNSj1+Lswf2+v+bDCYUqiTJ08e0Z6Qw1qJbiRxSeE4VkdauL12O/MDOZT7Ahn3cdUUuHoKhPqpJMC1WfiO8BF5KkLk2QjBMwdeStVgMOzPzEBo0FmjD3Uyvn0VkbNF5FERWSsiE9Pb/kNETs++eAcOW4RrSqfw6THllHkHdoeR4xFCHiHpKLvj/b87sEs9BE4N4Ox2iD4fMQF5BoNhRJGRohCRy4C/4JY6nQK0m+Vt3HxPBzUFHi8XFozFEqE+mRhY7WtVrlgB312T2XGeSV78JwRIVaeIL41lPK7BYDAMFZnOKL4JfEZVv0rn6s6vA0dkS6jB0FoTo7VmcBfanYkYn922hocaMs8YIiIsLnFzQC1vyEzReGd48Z8YwHt43ynJDQaD4UCRqaKYCbzWzfYWIG/w4gyeaH2Sl7+/EXUGvnxT6vGxMJjLvbureLOtKePjL50AY31w8wYynpV4Z3ixcizUURKbBzarMRgMhmySqaLYgZt/qSunABsHL87gCY/zUfVaA2sf3DngPkSEr5dOZqIvwE92bmZnIrMZStAWvjANVjXBU32nnuqW5JYksZeixFfEjbIwGAzDSqaK4i7gVyJyYvr9RBH5JHAjcEdWJRsggUIvE04uYOktW2nYNPBo66Blc924aaSAG6o3EXUyi3O4oAxmhOGxAeorz1QP3sO8JN6OE3+rex9wg8FgOBBkpChU9UbgYeBp3Cjt54E7gTtV9bbsizcwTrp+Bt6Qxas/3jSofsp9Aa4pncIkX4D98xX2ji3CHQvh5vkDG1tE8B3rxzPDTSIYX2UM3AaDYXjIOI5CVb8rIj/GLTZkAWtUtSXrkg2CULGP9/3sMHLKBm8UPjqcz9HhfACSqngyCPUfF3DbNiUUS1z32Uxw80L5wVHiq+J4prn2C4PBMPKYMmUKd999N2ecccZwi5J1MnWP/buIXAwkVXWZqr4x0pREO+OPzSdvUhBVpaV68Hfju5Nxvrh9LS+37MnouKaEct7rcGc/U3t0RSzBf2KA4DkhoyQMhlHM5ZdfTllZGXl5ecyaNYu77757uEXaS6ZXngjwe6BGRH4jIqcMgUxZZfmvtvHIpW/RVju4df5cy0NAbH5eszWjTLN5XuG0YvhjJWzqRxry7hBLsIvcZIWJ9WYZymAYKE7EIb4mjhMZebnVvv3tb7Nlyxaampp45JFH+N73vndwVrhT1Y8BY4H/xK1G94yIbBWRn4jIyCjF1IUZHywhEXF48bvv4Qwi4tlnWXy/bBp5ts31OzZS10OSse740nQIWHDjez27yx73ojL/uZ4fx72oqCqpXSnib8aJvRkz3lAGQz/RlBJbFaPtoVbiK9zn2KrYkGVByLQUKsC8efPw+92sECKCiLBx44hwJs08hYeqtqnqH1T1HFxl8XPgg0D/M2AdQAqmhTj+W1Op/ncjb/2mclB9FXq8/KBsOq1Oiut2bKTN6V8tiTE+4XNT3SC8F3d336a1j65aU2mbxQkBPDO9JFbFiS83ysJg6AtNKW1/byWxOg4p9j4Sq+K0/b0168piMKVQP//5zxMKhZg9ezZlZWWcc845WZVtoAx40VtEAsBi4AO4sRXbez9i+Jh50Vimn1fCm3dup+q1hkH1Nc0f4rvjpuERychl9tIJMC0Er9cPanjXZnG833WdfSdB/N9GWRgMvZICjWjnXBIdt2exdtiSJUs4//zzue+++/YqgX/96180NDR0++haJvX222+nubmZJUuWcPHFF++dYQw3mRqzLRE5U0TuA2pwYyeqgTNUdepQCJgNRIQTvjuNcUfldV+pO0MqwnncMuEwijxeUqr9ulB7LeH+o+CaWYMXoN111jvPi4TdKarBYBh+BlsKFcC2bU466SQqKyu5444REZ42oMjsfwC5wKeAcap6laq+lHXJsow3ZHP23fMoP64gK/1ZIsTV4YfVm3hgT/+i6vLS6cM3tyrV0cHNAkQE31F+fPPdO47UnpTJOmswDDODKYXalWQyedDaKK4DylT1YlV9WFUPqpBhEUFVeevuSl7/2QD9VTvgRci1be6vr+axxrp+HRNJKZ9YATesyzwPVFfaZxIaUyJPtBF5OoJmkN7cYDgksEGC4ua47oinh+2DoL0U6ksvvcQ111wDwOOPP05LS0u3j/Yyqbt27eLBBx+kpaWFVCrFk08+yZ/+9CcWL16cPeEGQaZeT3epasMQyXJAEBFiTUnWPFDN+r/XDLqvr4ydzNGhPG6t3cYrLQ19HhO0hc9MhiW7B54Haj85/IL/2ADOrhSRJ9pw2kae65/BMFyILYQuCuNd4HOVQvrhW+AjdFEYsbO7dDuQUqgiwh133MGECRMoLCzkG9/4BrfccgsXXHBBVmUbKJLpXa2IeIBjgElAp9BnVf199kQbGBUVFbps2bJe2zhJ5akvrGHn8ibOufdwxi7IHdSYUSfFNVXvsTEe4UdlM1gY6r2/pKNcugx2x+Efx0GuR5j/XN/fw+rFvf+gk1VJoi9EkIAQfH8IK88E6BlGF2vXrmXOnDkDPt6JOCQ3J/FM9WAFD/7/x2A/j46IyHJVrehuX6bG7NnAWuAl4I/A3cDvgN8Atw5OzAOH5RHe97NZhMb6eO5r6wYdjBewbH44fgaz/WFCVt/zWI8lfH821MXhV+klyHAfh/W1H8BT7iH4gRCagPiag2pV0GA4IFhBC99c36hQEgeSTHM93QIsxy1StDP9nI/r/fS9LMo15AQKvJxxy2we/dTb1KxsYur7iwfVX57t4cbymXvtBk2pJHl2zx/v4XnCZRMVj7i2itdPzc701y62CZ0bQkJp+0VKsz61NhgMhxaZKoqjgVNVtVVEHMCjqitE5JvA/wILsi7hEFI0K8xHHjsSf76378b9oF1JPLSnhr811HBT+SzKfYEe239zBkPi2mrlundL7UZuz0wvvrmmap7BYBgYmc6/BGhPdFSLG5kNUAnMyJZQB5J2JbHthXreumdwkdvtHB3Ow1H41o73qO6l6FG7kli6R7lv2xB4K9kguRbxpTFir0cHVfXPYDAcumSqKN4GFqZfvwF8S0ROBX4AbMimYAearc/Xs/xX29j46OBdkSb5gvykfAYxx+FbVe/1WSHv8Rq3bOpbjdm9kItHCLwvgHeel8S6BNHnjPuswWDInEwVxY/ZF9v8PWAibvGiM4EvZVGuA84J35vGuIo8lly3geo3Ggfd3zR/iJ+Uz6TNSfHtHRuI95Lu42szYKwfrl0L0SwHzYkl+CsC+I/3k9qRIvbvaFb7NxgMo59M4yieVNWH0683qepcoBgoVdUXhkC+A4bttTj9l7PJmxTg2a+9y573Wgfd5wx/iP8eP4PLi8bhs3r+qHM8wg/nwOY2+OVGcFR5qmk3ThZzOHln+QicGcR3lBvJbfJDGQyG/jJoHzFVrdcBXnVE5CwRWSciG0Tkml7aHS0iKRH5fwOXtG/8eR7OvG0OdsBi89M9pHnNkFmBMKfnjgFgRVsTO+Ld39EfXyRcPgEeqIT7dzZz866tvNrakBUZ2vGM82CFLNRRos9EiDwfIdWWxYxoBoNhVDJszsQiYgO3AWfjllW9VETm9tDuZ8CTB0KunPEBLnhwIYs+NzGr/cYdh5t3beW/qt5jew/K4svT4fuHwbPRbQDcXVeV1VkF7MvLn9qRIrUtSeRvbURXRk2eKINhkEyZMoVnnnlmuMUYEoYz6uQYYEN6CSsOPAh0F6/+n8BDwK4DJVio2IeI0LCpjef+ax2JLNx1+yyLH5XNIIXyX1Xr2dhNlbyALZTmN9DkJEmlbPYkk1mdVbTn5U++k+iwEZJvJWh9qMUoC8Oop7kqyp/PWk5z1ci11b333nsEAgEuv/zy4RZlL8OpKMrpXMOikn3utgCISDlwEXBnbx2JyFUiskxEltXWZimBEtC4JcLWZ3bz3DfWkUoMPn/SFH+Qn5fPwiPCN6veY02kc7lxR5W766pojnupqZ5LfXNhdmcVPeXlB4iAJo2iMIxuXvnhRlprYrxyw8jIytodX/jCFzj66KOHW4xODKei6C7SrOuV6hbgW6ra6y19OllhhapWlJSUZEs+Ji8ew4nXTqfqlQZeGmQp1XYm+gL8onwWBbaHJV2SCL7S2kBjKonHE8Pna6NhzwRqIp6s2yq6RdLZdR018RaGUcn2JXvYtbIZHNj1ZjPbX94zZGMNpBQqwIMPPkhBQQGnn376EEk2MPqMzBaRe/vbmap+OoOxK3Hda9uZgFvvoiMVwIPpwLRi4BwRSarq/2UwzqCYdXEpsaYkS3+5Fdu/gZN/MAOxBhdNXer1c8uEwwin80JFnRQ+sbi7roqoOohA4Zit7No5m521k/m1fzMnhAuweojizkqis3TX8RUxnHoH/8kBkw/HMGpIxhxevn4Dyai7MpCMuu8//OhRePzZ/Z2vWLGCCy+8kNtvv53zzjuP8847j5dffrnbtieddNLeKndNTU1cd911PPvss9xzzz1ZlWmw9CeFR9db9FMAB1idfn847swk0+JFS4GZIjIVqAIuAT7WsUHHqnki8jvgXwdSSbQz/4pyklGH6qWNpGIOnuDgE9jnpvNA7U7G+WrlehYGc6hN7kvkZ9tJisZsoW7XDNbvKuWVkgZOzi3s1IemlPg7cRKr3OPiK2J4F/jwzfN1n98pnZdfo13KP3rcVOXYYOVbJN5NEPlnG/6TA3jKMs3yYjCMPFbdU0m8pfPCRKI5xap7Kjny85OyNs6SJUu45557uP/++/dWueta7rQnrr32Wq688komTsyuI0026PMqoKofbH8tIt8GIsCnVLU1vS0M3MM+xdEvVDUpIl/E9WaygXtV9R0RuTq9v1e7xIHmiM9OYMGV5dhei2Q0he23spKnKdfyMMsf4unmemb6QhweCO/rNx/eoJn6WIDiLhf+dsO0xjpf9BOr4iTXJ7rNs9+el7+jcgE3L793rqtcvDN9WGNsoi9FiT4VcRXPQt+gZ1EGw3Cy5oFqUtHOdsZk1GHNn6qzqijuvPNOTj311IxLoa5cuZJnnnmGN998M2uyZJNMbxe/BJzeriQA0gkCbwCexY3c7jeq+hjwWJdt3SoIVb0iQ1mziohge4VkNMUTn11D8Zwwx35r6qCVhc+y+Pa4qdxZV8k/G2sZ5/PxjbFT8KcD9P6jWBG6SR7YbpjuamPvWDC+m4mP2IJ/gR/vTG+Py1V2kZuBNvbvKIm343ineZF8oygMBy9zP1bG6vt2dFIWnoDF3EvLsjrOnXfeyc9+9jO++tWv8stf/hJwS6EuWbKk2/Ynn3wyjz/+OC+88AJbtmxh0iRXabVXuluzZg0rVqzIqowDIdPFuRxgfDfby4DQ4MUZ+dh+i5L5Oaz5005e++9NWTH82iJ8vngCnxlTzsstDdxXv89UY4kgImxtU65eqdRlKVdTX3n5xSsETgoSuiCMle+2SdWa4DzDwcmCKyfgy+l85+TNtVlw5YSsjjPQUqhXXXUVGzduZOXKlaxcuZKrr76ac889lyefPCDhY32S6YziIeC3IvJfwOvpbcfhBsQ9nE3BRioiwjFfn4JlC6t/twNNuXmiBrs0IyJ8qLCUib4AcwNhwE2z0T6TiKZgeQN84234zRGa8Rc3UNqr5CWrkkSfieCZ5sF/bADxmRmG4eDB47c46QczeP7r60hGHTwBi5Oun5F1QzbsK4V62mmn4fV6ueGGG/o8JhQKEQrtu9fOyckhEAiQTS/OwZDp9eZzwC9wq9q1F3FI4toovpE9sUY2IkLFVyYjHmHV3VV4c2yO+dqUrPR9TDgfcCO5f7BzExfkl3BMOJ/DcoXvz1a+vQZ+9h58Z3rfhulsYpfZ+Bb6iK+Kk6ppxX+iMXQbDi4mnlTI2CNy2fFGI2MX5TLxpMK+D8qALVu27H1dVFTEW2+9NeC+rr/++sELlEUyTQoYUdXPA2OARcCRQJGqfl5V9w81HsWICEd9cRLH/tcUDvtQadb7b3VSNCQTXF+9kYf21KCqnDdO+NQk+HMV/KmaA1owXmMKPsG/OAA2RJ+KEFvZe/p0g2GkceJ10wmX+jnx2unDLcpBxYBuCdPG7FVZluWgQ0SYd7lrslFV1jywk1kXjsXbnwLXfVDo8fKLCbP4ec1WfrO7ik3xCF8qmcRXpgtb2uCxGvhoOX0apgdLVxdcAM/hXoiy13aRLUZb4XvDyCO3PMBHnzhquMU46MhYUYjI2cAXganAB1R1u4j8B7BZVZ/NtoAHC7vXtvLGLzaz4Z+7OPO2OQTHDL70aMCy+e64qfxpz07+UF9NUpVvj5vKT+e5nlCetF2k3TCdbXpywU2+nUACgv9oN2V5fG0cbXbwLfIj3sxnMhnHgxgMhgNKRrdtInIZ8BdgPa6iaLdT2MA3syvawUXx3BzOuGU2jZsj/OsTq2naFslKv5YIlxWVcX3ZdC4vcl35ghYEbaE1qVy7VqmJDVHKjZ5yQ3V0wQW01SGxNkHbI60kd3SXSKpn2pVRYnXc7S/9SKyKu0rKJCo0dMHUUnE5kJ9DpvP7bwKfUdWv0vny8TpwRLaEOliZeEoRZ/1mHvGWFP/6xGpqVzdnre9jw/lM9AVQVX6xayv3797B9ojy1C74/FvQPIwJ/fwVAYJnBcGC6NMRoi9HXCN7f+inMjIYAGzbJpFI9N3wECASieD1evtumAUyVRQzgde62d4C5A1enIOfsQtyOe+++fjyPMSbs3+Va+/xj3t28vumjfxoXoJNrfCVVRAfxmR+dqmH0PlhvPN9JDcncRrNFd6QfQoKCqipqcHppbTwaEdVaWtro6qqirFjxx6QMTO1UewAZgFbu2w/BRi5eXsPMPlTglz80BFYXlcP177dTPG8nKyk/PCI8PWxk5kXyOH2uu1sia/jqpmTuX19Lt9eAzfOU+wsjAP0KzdUR8QW/Ef68c72YoXcc0+8G8ceZ2MVZNlf13BIUlxcTGVlJevWrRtuUYYVr9dLaWkpeXkH5v48U0VxF/CrtPEaYKKInAzcCFyfTcEOdtqVxK5VzfzrE6uZef5YTvjeNGzf4L15RISz84s5LBDiv3du5vn4Zr48fR5/3G5TE4XxwUEP4Y7Tj9xQ3dGuJDSuxFbGIa545/nwLfDtb+zOUBkZDm0sy9qb5sJw4JBMDSIi8mPgq0AgvSkG3KSq12ZZtgFRUVGhy5YtG24x9qKO8uad21n560rGLsxl8c2HESrOnodS1EmxNR7lsECYPXGHpCQp9rjrll1nMMe9qLT2siIUtuH1U7ObylyjDrHlMZIbkkhI8FX48UzxdJKtOxdc38LelZHBYMguIrJcVSu63ZeJohCRSbh1JAK4da4tYA3QCkxU1W2DF3dwjDRF0c7mp+pYcu0GfHkeTv/lYZQcnpv1Mf7RsIvf7d7BVGcCk6WIL02XThfk+c/1/V2vXtyDolDlmeZ6zsgt6rEuRm+kdqWI/TuK0+gQuiiMFd5f2Zg4CoNh+OhNUWT6b9wMFKtqm6ouU9U3VLUFKErvM/TA1DOLOfe++VgeoXZVS98HDIDjwwXM8IdYI9t4uG0zv9iUyJoL3SutDdy8a2u/q+05EYf4mjhOxDU62mNtgueGCJ4dwgpbqLqziPb90HeiQoPBMDxk+o8U9i9XCm5W2ZFbrXyEMGZ2mAv/upA5l44DXCN3Mrr/WlCmBeAdVZ5q2k2xx8tPy2fy6aLxhIJNPJl8l2s3Ng5aWbTX8gb6rOGtKSW2KkbbQ63EV7jPsVUxNKWIJdhjXKODNjrEl8doe7iV+OqYqddtMIxg+mXMFpFfpV8q8BMR6ZjXyQaOAVZmV7TRiS/H/chjTUmevHoNOeMDnPbzWeRP3meB7lgA/qw75/XZZ/vdfsiyOCmnkI8UjePocB7f2LaNx2ssxlnwxWkDl7m9ljdAQyrJq60NnJSzf0K1TIopWQU2oQvCxJbHiK+Ik1iXwLfIj2eqxxRJMhhGGP2dUcxPPwSY0+H9fGAGsAK4YgjkG7X48zyc8uOZtO6M8Y9L3mLj47VA5gXge7rbn+oP8Zfps7ioOJfZOfCn+mqCoXq6nxD23X9U3SWiqDo9zyoyDJ6z8i2Ci4MEPhBEAkJ8WcwE2BkMI5B+zShU9TQAEfkt8GVVbRpSqQ4RJp1axAV/WcgL31rPi9e8R9VrjVS+XJ9RAfje7vZty+J7h0FSlW9WNVFU3Eo0soeG+omkUv3zvOrYfzu9zSoGgmecB/tcG21RxCuoo8ReieKd5cUuNanMDYbhJtM0458ySiK75Izzc87d81hwZTn161pJtHaOOG0vAN8d/b3b94jwpcKZNOwpx+dvYWzZWnJya+hrdtG1/3Z6nVUMEBHByk3HXzQ5pKpTRJ6IEHmmjVRd99OMrgZzg8EwNGTsXiIiHxWRu0Tk/0TkkY6PoRDwUMDyWlR8aTLNlZEeC8B3R293+12ZErZItIxlV/UcYtEc8gqqsT3xTm26Zkd/tbWB2mTnNu3sSsb3HycdPLdfkJynh+09YBXYhC4O4zvSR6ouReTRNiLPRVzbB70bzA0GQ/bJaF4vIj8HvgI8j5vOw/wzs8i8y8bvVwDe9km3BeD7uts/IVywX7zD8tOEp3b5+PaaaeSGY/zraD8TgsL/Nezi1JxCCj2dE4yVevxcmN9zLplSj7/T+4FGcneHeATffD/ew3wk1sZJ7kiBz1USrQ+3QJw+DeYGgyE7ZLoA/AngUlX921AIc6iz4MoJrHuohkgHRZGKK96wjTrayRuoP3f73dkQzhwrFPvgP1cFeGoXfKAsyt11Vfy+fgeXFZZxfkEJXnEnmjMDIWYGQvv10RtiS1aLKYlP8C30413g1g93Ig50V0uxo8HcpP0wGLJKporCwrjBDhldC8DbfqFwZpilN2/Fl+PpVHK1093+ziT212tI/aIUxnn27u+JIwuEh49VxvpAJMCNZbN5oLGK3+yu4rGmOj5TXM6xofxBJTHMdjGlvbIIPUfzGAyGISHTW727gMuHQhCDS3sBeCwoPTKPD/5hPot/cRgzzi8BoHFLhFTCYWYgxFUlE7iqZAKTbm1FalNMuq1177a+ZgKlfje9x46I8tkVAUoi0/nBuGkI8MuabfstaY0UxCuuojAYDAeMPmcUHYLtwFUsl4nI+3FrZneqIKKqX8queIcmJ143nceufIcTr52OiDDljDEAJCMpHv/MO/hybE743jTGHZXfbdzFxJP677Za4ofTS+DurbCuJZ+fzcmjUaMELZuUKvft3sG5+cWUenueoRxQ+sg2G18dQwIW3plexGc0isGQDfpMCigiz/ezL1XVxYMXaXCM1KSA2WLbi/W8/tPNtOyI8U5FCZPX7iGndZ/nU0uel3u/dySBgNVjJtiuqCp/qYKfvgfj/PDL+TA7V1gXbeUbVetB4Zz8Yi4tHEeB58BU1OpV3h6yzXrmeIk9FyVVnQIveGd68c7xYeWY3FEGQ19kLXvswcBoVxTgziy++N3tVDy3A9HOKzFxr8XyxeN5/ayJPWaChe6zwa5sVL6+Go4rgh/PdbfVJuL8cU81TzXtxicWFxaM5WOF4/BZw3/x7SnbbKouRWJNnOQWV4H6j/fjnZk9e4nBMBrJZvbYrCIiZ4nIOhHZICLXdLP/MhFZlX68KiILh0POkYYnaPPKeZNJ+Kz9lut9CYdFL3Ufd9GR7rLBHpEv/PkY+O4s931VRAmJl6+MncyvJ83lmHAer7U27K2g11vAXXuiwmwG5XWlp2yzdrFN4JQgoQ+F8R7u2xvdndqVIv5uHI2Prpsjg2GoyTSO4t4edilu9tgNwJ9VdUc/+rKB24D349a4WCoij6jqmg7NNgOnquoeETkb15h+bCYyj2ZWnDqeiud34E3sMzw7AmuPLO71uK75oTrGXBT79imBr6yGtpRbXnVeXoDvjJtG1HGwRWh1Unxx+7sszinkwoKx5Nqdf0pdExUOB1bYwn/kPttKcmuCxJoE8eUxPNO8eA/zYhf17ktramQYDJnPKEqAi4ELcZMBzki/vhg4DPgmsE5EjuhHX8cAG1R1k6rGgQeBCzo2UNVXVbU9K97rwIQM5R3VLD2jnFig84VOFBa+VsNpf9tEW133cRbd5YfqiiXCNbMg5sDly+HXm5WkowTSS05tToqpvgB/3LOTT2x5m3vqqqhPur4NmaQlP5D4KvwEzwnhmewhuTFB5J9tRJ6PdNvWRH8bDPvIVFG8AjwOTFDVU1T1FNyL92PAU8Bk4FHgF/3oqxzY3uF9ZXpbT1yZHns/ROQqEVkmIstqa2v7MfToIOW1ePrS6STSdbgTXovHPj6D1ceXsuC1Gh46/02iDZ0c0zLKBntUgfDQMfD+Erh1M3xyBdSll21KPD6uK5vOHRPncEw4n4caavjk1repScT6pYiGAxHBLrEJnBQk/OEcfMf48ZSn62M4SuyNKKldKZykQ9vfW0msjrueVelHYlXcTaOeZWVhclYZRjqZBtx9GVisqntjY1W1LV1H+1lVvVFEfgY804++urO0dvsPFJHTcBXFSd3tV9W7cJelqKioOKRu+bbMKWTHlFwmvtdI1bRc1h9ZwvojS1hxahl30UigwPVS2vrcbsYfV8DrTlNG2WDzvcKNh8P7dip/3QH5XX4xU/1Bvj1uKp+Ml3Hp6gbO2OqldHwVHs8+RXT99ipqqvMJ29JvT6yhRvyCb84+A7fT4JB4L0FibQLJFbRN9/81Zjn6u6v3VnxFDO8CH755pla4YWSR6YwiB9g/8RCMS+8DaKJ/CqgSmNjh/QTc/FGdEJEFwN3ABaq6OyNpDxGe+eg0Wgp8PPuRfdWJGkqCzPmo+1U1VUZ59mvr+Ou5K3j8nk0koilCuxzOvTpKaJfTr2yw54wT7l0EXktoTipfWa1sbN3XfrwvQH1jKYFgA5bVeRZj2QkCwQZaR3CtCbvIJvyRHPwnBJCADHnkd3uRpwM1azEYBkOmiuLvwD0i8mERmSIik0Xkw8A9wMPpNscA6/vR11JgpohMFREfcAnQKQOtiExK9/txVe1Pn4cMHTO9NhUFuOe6o2gqCnS7P29CgHN/ezgy1cuMe2Oce3WMU38QJ1gLFXe6F/Vus8F2oT2NxqZWWN4A/+8N+NVGJbL3oqbkF+7Asjpf5CxLKSregj8wsjPUi1fwzvQSPCPUc/S3kp2LeIZFngaLWd4yDIZMl56uBm4G/tDh2CRwL/CN9Pu1wGf66khVkyLyReBJ3In8var6johcnd5/J3AdMAa4PX2RSvbk53uokekSTumiPCpun8Xzr+7Af9se/BsSCFD6Lnzo3Rz02FCv+aE6sjBf+Mexyi82wG+2wmM18J1ZSiDYiG13b0AHsO0YABtjbUQch3mB8KDySQ0ZNkiom+hvAIW2/2vFM9mLd5oHe9zILqxklrcM2WBAAXciEgam4953bVDV1mwLNlAOhYC7wZCMOfz1nOVE6vYtD1leYeaFY5l/RTl5EwK9HL0/S/coN6yDCUF4vamNYLi+x7aR1iJWnBzmpzs380LLHib7ApydV8zpuUX7udcON91Ff3sXeLEKbFJbkyS3J/FM9hI4MYCqktqZwi61+13vW+NK659boLsbfAvCH80ZdAqSTjXMO85cbJCAmJTshk70FnA3oH9nWjGsGpRUhmFh1T2VxFs63yaro6x/qIb1D9Uw6bQi5l0+ntJFuZ3u9purojx25Tucc888csv3KZOjC4WHjlGak3DqyyESDX2nJf/K2EksDObyRFMdd9ZVcs/uKs7PL+HKMeX7RYsPF72lS/dO8qJJRRPuTZZT7xB9KgI+8Ez04JnkwR7vQTy9nEMfOauykiq9fRmrqzIyKdkNGdKfpICPAJeralNfVexU9fysSWYYEtY8UL1fFT1NgTfHYs4lZaz7aw1bn63ntJ/P4lL/mL0G6Ivv2MjEnTF+/tWN/P3quYBrB3n9VMFrCUU+6K8FOGDZnJ1fzNn5xWyKtfF4024Kbe/eIL0VbU18oqiM8b7MZjdDQU/p0sUjexWBVWAROC1AcmuS5LYkyY1J8EDwzBB2iY2q7rfEls0iTwbDUNOfGcVu9l0BjNfRQc7cj5XtV0XPE7CYd9l4jvz8JBZeOYGNj9Yy8eRCWl+DuW/sYuo7exi/uQlLoXxTM1PW7mHLnMJevZjseIqUz0ZSDmr37DMxzR/iCyUhHFU+tfUdAF5o2cMLLXuYEwhzem4RJ+cUkj/ClqY6IrbgmeTFM8mLppRUTYrUtiRWQTq+5e04ycoUnok2ngkeJN9CRLJe5MlgGCr6/Pep6qe6e204OOmuip4312bBlW7QuzdkM/vD49J7lJLKVmauqt/rBORNOLz/Txu499qjSHk7X9TCNnuVxwX3rMMReOFDU2koCQIQstyAv47LSse9qLSmIBDcQ+GYJJYFjiNE2gpYFY+wNrqdO2q3c9vEOUzxB4fkM8kmYgue8R484/f9tcQvkFTiy+PEl8eRHMEzyYP/aHfGlO0iT3s5EMtbhkMCc/tyiNFeRc8TcL96T8DipOtn4PF3/1OIBT0ku6y1h1qSHP1s1X5tXz9VWL1YeMzbwPRtzUxd18h/3LyKv3sbWL1YuHUhXPRv+OdONx0ItCuWdrdaV3lZluL3t7Jr52HUVM/mY0VlTEovQ/2mrpIfVm/k+eZ6Wp0RHJjRAe8sH6EPhgn9vzD+4/1YBRbaIQYl+nqU+NsxUvUpspnNuX15y7vA5yqF9MO3wGcM2YaMyHg+n07O90VgKvABVd0uIv8BbFbVZ7MtoCH7tFfR2/FGI2MX5fZa6GjRkmq8yS5xEQqLXqpm9fGl/OPSd5l+TgnTzykmOMZHMubw8vUbSKZnLKmow9IfbGDKo0cBgi3wnTVwx2a4YpLbrxuk1zmgwLKTBIJNRCMFXFa0z0AesmzejbbxamsjXhEmegNcUljKKblFWfp0hg4rbGHN8uGdtW/2oCkltStFcp0Dy+NIQLDLbDwzvXjKBr/cZpa3DNkgo1+MiFwG/AU3oG4q0F7FxsZNCGg4SDjxuumES/2ceO30Xtu9eXIZiS5LTHGvxZunlBFqTiACb9y0hQfPWMaTV6/hhW+uI9bc+aKfaE6x6p5Kji4U/nYM/M98NxXIDesAlLyCfbOJdizLIb+giq4G8suKyvjDlMO5qXwWi4K5bIpHeKypDnALMD3aWMuORGxAn8lwILawYnGcz1RsYttRKezxNqnqFNrofh5Oq0P0tSiJzQmctoEHy/WUkt1g6A+Z3rJ8E/iMqj6YnkW08zrww+yJZRhqcssDfPSJo/pst/SMcua/VtMplXk8aLP09HJSXovzH1jIno1tbHyslk2P1dGyY/+LdDLqsOZP1Rz5+UlYIiwugdOKlWUN8Nl32vB44qhCU0M5ofBuvL4oALYnTiDYCHSe8VgizA2E2Rp32+2Ix3BU2RaP8r+1bp7JiV4/x4bzOTqUz7xgDp6RGNjHviSNe/wpbvRUcu+8ea49KP1xO00Oyc0JkuvduBfJs7BLbXwL+q7c127/6Yl2rzWDoS8yvb2YCbzWzfYWIG/w4hhGGt1lqH36kumdDNmF00NU/OdkPvzYkRz24VLsQOeflVgw7qg8Yo37gvxEhKMLBVRoaS6hcc94WpqL2bVzDjur5lJfN4nmphJSye6NvB0z1DY5KV5tbWCKP8hvJ8/j6uIJlHh8/F9DLd/a8R5L2xoB2JNMUDPCZhvdZdpt94gC8JR5CF+SQ/C8EL4KP1aekNyS2JtiJLExQfSlCIl3466Nw9k3A+srt9ZIzr1lGFlkOqPYAcwCtnbZfgqwMSsSGUYM7V5MXTPUbplTuHd/R0SEY/9rKtuer+/kVQWw7fk9PPDSUmZ/eBzHf3tf8sJEYl+QXu6eKLOX1bFs8XgiqTE9ytVTqvQTwgWUef1cWDCWCwvG0uakWNnWzBHBXAAeb6rj9/XVTPD6OTKUx6JgLguCuYTt4XH/6e08OnqGiSXYY2zsMTbM83WKy9CYGxWe3Jxe7vPiplI/3fUQs1VJjdDZlOHgIVNFcRfwqw7LThNF5GTgRuD6bApmGH46Lks0Hzadx658h//65XR+WN7zhafdq+r5r68jGXXwBCzed9MsgoVetj5bT+5E13sp0ZbisU+9zcll+WyaW0D11Fze/+AmJr7XSPnGJv4vHdTXzp8qFVvgA2NhVayhX6nSQ5bNCTkFe9+/L7eIoGWzvK2JJ5t280hjLQGx+Nu0hXhEqIpHKfJ4OW2JdUCWbDrOJno7j650DN7zzfXhneNFW1yjuFObQmOaTiWi/LI2SmnK4R2fzRqfxVqfzXs+i7hRHoYMyEhRpOtN5ANPAwHgeSAG3KSqtw2BfIYRQn9tGrC/V9Wkk12PpJLDc/e2idYn8OV7WPRSNRXP7yDhtbCTDpbChA5Bfe2zlqd3wdIG+Ol6CActrFAuhc0NnP79GM//0Efb2P3LunZlvNfPRQVjuahgLAl1eDfaSnUivtd+8bOaLWyMtREqDmHHcohHc4jFw6jT+W+SjSWbrrOJvZ9LD7OK3hARJFewci2Y7u2076WgTUUMKqIpzm5zldLrAZuvpmNbEpsSWHkWVqFl3GUNPdKvpIAi8j7g1XTJUkQkBMzFtXGsUdWWIZQxI0xSwJFBT7mhupJoTbF9yR6WfH9Dp2hxX77N5NPGUH58AWVH5xEo8rK2BX69LcoLtTaO42XGuzWce9cmahYIL10bwBblO+OmDrhG95ttTayMNPP76hZ8/jZElEhbPvV17lJZMLSHRDxIMuln9eLBeQ+93LKH/965uaecgIM6j3bmP9f5v12SdJgbd4hY8EbAQ8BRnt/R6jqWiZuKxCqy8Ez14ikfuZHwhqEhG0kBnwOiIvJa+vVzwBuqasxhhm7p7wzEG7Zp2NS23/Zkq8PGx2p57/92AVAwLci4o/N588QSwrl1FK4XKh6vw1JwmsLs3DoPyx/la9U2zx2v5PSWkK8HFoXyWBTK4+a3FMTB52tFVQiFdxON5FJUvAWAVMrm2h1hZgfCnBguGFDEeKnHz4X5Y3vdn21qPRYvevYpuKhA6OIwzu4Uqd0OTn2K1I4UVpEN5eC0OESeaMMqtLAKbfe5wMLKt/qdJdcwOuivopgFnAa8D/g8cAPQKiIv4yqN54Hlms2wUsMhQ3eJCp2k4s21OPfew6le2kj10iaq/93InmOnYNeF+MTNS/HF3WPCDXFmLa/j3YoSkiKcvAQqCpQfzIbxwQFe0NQiHsulpLWGi/5nE3//8jRqqmfj87Xi87dSE2xjWX0TJR4fU/xBtsej/KG+mln+EDPSj96M5DMDIWYG+s60O6SIYOUIVo6FZ/K+zXs9pxywx9o4exwSVfG9IS2BUwN4pnhxGlIkt6ZzWgWEVG0K73TvkMVqOBHHBA0OE/1SFKq6AdgA/AZARGbjKo5Tga8DPwUagZEfHmsYcfSYqPBj4ymZn0vJ/FwWfNoNqLv+eTj6mUrs5L62BXvinPXARo5YspMHv7aAj0+EpTtSFPncC/WfKpXNbXBiEVQUQLjfsw3l9D9vJVgHp/9lKw9+6miSiSBtrcXcdYTQmkrRfmNdl4yzJtrCiy179h5d5vVz7bipTPOHaEglSKoyxvbul0m2v8t0mdIx91ZP+7ujfbZg5VkETnFnS5pSnEYHZ4+DVeoemNrtEF/ZuVBVYnkczxwv/qP8aJuiMcXKswZVW8MUXxp+BlqP4l0RqQfqcRXEJeyrmW0wZERfiQrbcS+wyqIl1dhdFvcFKK52l7C+Og3+8MnlPJLvoXheDrtLwryRn8PDE8KkQh4W5CmnlcAVk3q/yMzaUEnJ+hSWQsm6FLM2VrF++j6ZOs4YFoXyuH/KfBqSCTbE2ngv1sbGWIQi2zUuP9m0m9/u3kGeZTPNH+LftQHaYgHaWot6TeE+GLIZTCe2YBfZ2EX7ztkzxUNsBa47SweFlFyXILXNvfNPvJ0OFAwKkmdh5Vn4j/EjHkGjDnRI194dnYovdRgjsSpOcn3C5Kw6QPRbUYjIGNylp9OAxcA0YDnwIvAR4OUhkM9wCNCdS21viQrfPLmMiud3dI4W91osXzwegFRCOeIzE6h9u4W6d1oIVe3mfKDk4+VsumgSb+xIUXXvDraeEqZwVpgbG32UBoUj82FhvtufnUix+I9VeNLxeZ4YnPaHSjZ+p4yUt+clpQKPlwpPPhXh/E7bjwvnExSLTfEIm2MRvKHdFISVkmUW47c0s/b/2dQdFmFiaiO78wpIJAI0pwIjrvJfJ1JAlP0LIzluYSTPNC92sY3T5LiPRodUdRJs1/4SWxYjuTGJhATJtbByBSvfxne4G2SpCUUdNcWXRgD9+hWKyFu4dopluIrhy8Arqrq/FdJgGACZJCrsLa0IuIpn/hXle/dFGxLsXttKTpmfD04RaiMR/vl/lTz7d3f/pIBNbWmQGz44mR3T8/DFUsxaVY23Td1leQFR8EWU41/cxJsfmJnx+U32BZns22f0nv+cg9eJccWf3sYbd5CUTawAKG+k0NcEwOe3e7l/ynwAHmnYRUyVCV4/431+yjx+fFb/1+kd1QNePdAKW0hhz1dxz3QvVq6F0+zgNCupqhSpWmevoog+FyFVl+q+XCyAusuRgplRDDX9vV2ZCewBNuFGYG8wSsKQbU68zg3q6ytRYXtakfN+tx5v3Ok2rUhHAgVeyo8v2Pu+5PBcPv7qsezZ0Eb9+lb2vNfG7vfaOH6exboy2PbUHmY+4OaMSnqE9UeMIem1aM31kVcf4YGK/etqZI5w1JO1+KPuesrch1PMfThFzGfxxvmlvH16Id+Zv8835PmWPayNtnY4Gk4MF/C9Mtd198mmOsKWzTivq0S6GtLbqweGLGvQbrfZwlPmgS4ZcjumIPHM8CJ51t48V/uh0PbXViQsrlIKCfY4D94Z7nKf0+ggQQEv+9mFDJnRX0WRDxyLu+z0ceA2EakFXmh/qOrmoRDQcOiQSVBfT2lF+os3ZDN2QS5jF+R22r4YeO5Uh9deFUpXuRetGavr8cUcfnvNETSUBtn182re/8ZOxk0LwLgAu8YEmDA1wKJT88kP9n8dZNGS6k6zIgB/3OHoR3ex7MQpnJK77+L2ywmH0ZxKsiMRoyoRoyoeZYwnvUSjyp21lUQ6BO/lWDbn5hfzqTHlOKr8atc2AO6o3c78QA55tmfwF88hKIzU0e3WO92LZ4qHVFWy2zGw3FmJtina6uDscMASvDO8qCpt/2x1j/GAhAQrZOGZ6sE7y02DktqaXvYKWq4NZQAu1YcK/fV6SuDaIF4GbhARP3A8rs3iCuB2EalR1SlDI6bB4NLRk+eZj07jw7e+w7MfmdZp/2CZOCWXDf89Abm8Ck+9A+l6HB+/ezW1d02meI7N2OYwbduj1C9tgqjDRoGv3HgsxWHllGcqWbClkTET/Njj/ITGBxg3OcC4IzvnzezJ1vLmKWXdypVrezjM9nBYINxpu4jwhynzqU7G2JmIsTMRZ2cyxkSv60H1bPNumtNFnnanknx0y2qCYnF5URkfKiwl4qT4254aij0+ij1exni8FHt85Fp2r8qkvTBS7J0YrW/F9gbuhecH8GfJIynTMfZ66Cv4Twy4SqTNQVsVp81B2ycnMYi+GO08mBe883yIV7DLbZLrEq4CCaQVSdCNfh+MBxcMzzLgYBmopczp8Eh/dUzMllAGQ0909uQJwv+r4CtZHmNmIMTM8hDbb8jpZGB//w9mMXF6IUwHPloKgOM4bN2RYP2WGF+cYLGxFZo9FlZKqXqtkbZa16WzLc/La7+sYHIIpt67kY++10ZLvg/H2vcHgrStJW2Uz4SwbTPDduM3OuKocv/uao77eYxQHVRXWFR+NMDpuUVMTlcNrE3GeWDPzi6VP+CLJRM5L7+EHfEov63fQZHtpcjjpcj2Umh7mRkIkW97WDotwt3eSo6qDbOspJWrJk7gJDt7AYNiS7/HaFdsYgneqd799u/FB8HzQ2lF4iqRVFXSdcEVXDcd2M8+4j/ej3eWj1R9iugLESSQViQBQfyCd7oHq8BGY4rT5CB+QQKdl79G4jJgX/TXmO0BjsFdejoNdzYRALbhBtvdk342GEYN/TGwW5bF1Al+pk7ocNGaVw64xvTVtSne3BijpS5J0Iale6DZ42OMJ0pJVRt2UhGgblyQV86dRF1ZiI/c+g65DXH+cZcP/xgvOSU+xi7IZdbFrnLa/W4r3hybYJEXb6j3KdQrrQ0ElsYZ/6biiUHB9hStM5IcfnbOXs+sSb4gj0w/gj3JJHXJOHXJBLtTCQ4PuB7vzU6KzbEIy1NNtDn7rpzXl03jmFA+t+3aToOdYIkTo1WFm2q2sjrSwoUFYynz+tmdjFOdiJNveyiwPYQtu9930m5NDaV0fBV5e2LID5tp+aGP66mipjqfsC0DcgMWS7ALbSjs4oLbUTHYQAiCZ4Yg7npZWUVW+niwi200qmiLg1OnaFSxS22sAkjVJIk+32HGIu5ynO+0APdv3s6nbrGIfzRCtDyI5Rd3luIXPOUexCdoIi2LlxERBd/fGUUDEASqcRXCfwLPG7uEYbTTXwN7T8wvsZlf0vkuP3H8RE5ekl5CU2XO0lrEUWrHh2guCrBuUTGluyIc5o2zdmsc+61W6rcm2TJlLKV+mPm5d6DRTfBnBS38hV6mnVXMsV92w6tX3L4Nb9jGX+jl4bYqjvjfWCc334W3xbhnYSUnzNyXeNArFmO9PsZ696//cVggzN2T5wEQdVLUp5LsSSaY5AvwSmsDbU6KcW86WCnXO6x5osPjTXWckTsGvPDv1iZ+Vbttb38WkGt5uHnCLMp9Af7d2sjLLXvItT3kWR5ybZscy+a4cAGtKSEYqseyEhx1Z4JgHVTcmeCF73kIBBtpjRQM6HvpRLurrQOtjXGe/8smTvvINML5Poi6n7Hkd75YWwX23mDEdjomprBKbAKLg2hMOz1WaDPTbm2jIJXLnJoA8ZY4Vmpf3/aFIcRnk1ifIL4s/aV5XWVlFw+fH3B/kwJ+FnhOVd8bepEGh0kKaDjY6BiZnTPeTyQFobRh9R/VynstUBNVauJCTQyO2tTAp3ITRHbH+cPqBNqYYMeUXDadOo4Sj8MFX3oD6VLnvL4kQGuej2BznKLaKIkcyJ8cpKQwgC/Pw9QPjGHSqUUk2lJseWY3/jwPvlwbX64HX45NoMiLJ9D5QnXciw45pWuY8FaEE36RwBODpB9e/YaXyvlBWnbN5fVTLXYnE2yORWh0EjSnUjSmkjSmknxqzHhybQ//bKzlz3t20pxKEetgkP/r1AWc8KJF2YRVbqlcR/FGwNsG8TBEvT5qqudxy9F7WBdtI2TZhC2bkGWRa3k4JdedAdYm46RUCVk2QcvCK5294zSutP65BRx48W+bqdneQumkHE790FSwIPzRnEHbJcBdBvz6QyuZc2Nk72e19ltBbrpgIZIAjaaj2G0hVZ8itTMFcUXjineeDys8tGlLeksK2C9FcTBhFIXhUOK1eqUmBnUxqI3D7jjMzVEuK3J4t6aZb75o4Y+kaCr0sXt8GF8kyQV3v0v5tmYKjshhS43ijyTZdfZ4ImeNo6g2QulXVu43zqSvTmXGR8bhbG3jja+sxRe2eVeBnCjFaxVPBw/WRBBWXephj47lJwvz8IZsig4LESjwkow5JCMpvCEbyyv7GcvjjkOzk6Q5lWKSL8CxrzcwJn8T8x5KoraQCEE8LMTyoWaeRV3TND4zs5WnmnbT5qT2rhzlWDZ/m7YQgB/v3MSSloa9Y3gRyn1+7pzkRsG/f0UlR6RasJpSNG5pw44ouXVwJRMYMzWXj8xKcs08h4BlERCbgGWRb3uYmLbxtKZSeETwyf7n05EXd+/mnYvXEdgnCtFC4fCHZnHKmJ4LdR0ospE91mAwjECOL+ruwiSAxcK8Qm54ehtL/1FLq8cmEvbQluPFb8GiKycw57MT+eG70JCExgQ0NUOz+Pn4nUdwYb7Dzt1Jrl2WxB9JUu3Ppf51yK+zuHJePiWkSGyP46sP0pKb4q0TS6kvC1G2qYnjntnBUfcmgR08yw4Acn44m7zjCnFeq2fn99e7Utpg+W08QYuTbj6M8gW57HqtkdX37cATsNgcsDi10aF4p0NhVQpvB2WUsmH7ggBrF7VxVjzEp48pww5Y7KmL0VAfI2VBa00M22fxQe8YKkryiKhDxEnR5jj4OlzQoyjNdooGK0p8jkXKDzk7lTd+XEnFVbPRvJ3ctCvS6ROeGwhz84TDAPhq5Tq2JVx7hF8Ev1gcFcrjW+OmAvDD6o20pFJs3dRM+CovnigUv+sw46kUnjbl9y9tpurUOD6x8FkWPhEmegPMDbo2ojfbmpgVCBO2RvjS05ANLnIW8D+4ZqO7VfWnXfZLev85QBtwhaqu6K1PM6MwGPaRjDn89ZzlROr2XWWDJV4+/OhRPaZIaSfhuMkUm5PQktz3fEQ+zM4V5j+nzHqzjqTHIhGwifttEj7hpH9up3xTI3/66gJ8sRTeWIrdZSGiYS/5tRG+3tDANI/D5t0pntmewht3WHp6OY3FAWas2cNHXq0imHRobk1R3+LgjaV49v9NI5Lj5bA365j/+q79ZG369ZF4x/rx/7USKx0s2ZHyh4/Gn++h8bfbaP7HTmyfBV5BPRaVKeGdo0tY+PouAk0J/NEUAjgW1JcG2TklzNd/NYmYOlQ+vpum9RG8XouxQT+WV9joiZI8P4eY4xBd1kayLsEYn4+KvHwsj/DXxC6emuwwpjJCaJdryylar8z+p+uy/MQvfTRNcr8Lf6MiDizOL+Kz4yaiFly87S3+Z9qcIc82PCJnFCJiA7cB7wcqgaUi8oiqrunQ7GzcqPCZuAF/d6SfDQZDP8g0j1ZHvJYwq49Un0U7Iz3m3WoYG2TpqRBJuY+2FESdIBODQfK9QmlUkT0QTcERjvscnVzISZ8uZGJIeLNB+cQKKNgVIbchTsojvHNMCeUbm8jZE+Ol8yezY0Y+OY7DzmYvsTYomlTEz38UpAiHJTsdnqtysJMOb260cDwwzZ/Dl88uwZ9S1tQ7vNfgYCeVFe8bz4rTyql4tooTH92GKETCXuykUraplbNf9mMLnPFYlMNX7sZJKFXpU47ke/nXzOl4LDjxD2sZu6qBFuAFagAIlgbYdc0CcutiHHXfesbUdJ6dnPATi999Zx5njLWZ+aO38W5qA6r5I9UAnD03jOee/Z0MDiTDufR0DG4qkE0AIvIgcAHQUVFcAPw+XefidREpEJEyVa0+8OIaDAcnmeTRypS+8m4FbCHgeqHux7iAcEH3sYUALCpwswU3F/r58K3vkNO8b1bUkudlzXGlpLwWLyx2l5FUlYSG8EgIS4RxSeWMBCQdSCgkHEgcVcS83CK8ljC+VZnWBl9ZDTPe2s3kdxsQQEVAlZryMCtOK2P7rAKumgKOQuqbM/jkDDfX1yNVDit2K8mkcoQfUgotX5jB5yc4aFL561aHt/coji2gFv58H++cUMopf98CwKrjxtJQGsQqDuA4XtY0CdY5E/h0fhLHUR7erlS1Km35Ps5MDW/mw+FUFOVAxzliJfvPFrprUw50UhQichVwFcCkSZOyLqjBcLAzWDffnsg071amhG1opecxOkbiiwgdnZNyPEJOL1e4qWFhahhA2Ty3kMUPbyTctC9PSGl1CzumutH0/zltf1vQ+eUW55d33brvzv8/puzbevwzDrfVRtkNvOoRUkll0bJaTvjgZMbk+zhX4bHjBY4v3nvM3J5FP+AMZ5mo7qxwXQ0m/WmDqt6lqhWqWlFSUpIV4QyG0UR7Hq1sFkZqpz3vliMMKO9Wb7x+qrB6sfDPLxQx+chcsGByRS7//EIRqxcPLNiuO1Je4Y3P+0im4yaTfnjjai8pb3aD3cqm5lI8PgwCxeVhyqbm9n3QCGA4FUUlndN+TIC0i0RmbQwGwzDQ8W7+mY9Oo6XAl/W8Wx058brphEv9WZ8VAQSCDdQcJdQdJjgCdbOFmgohEGzMSv9Rgd22EBGoeH85oRwvFWeUE0lvjw5/8HWvDOfS01JgpohMBapwq+R9rEubR4Avpu0XxwKNxj5hMIwMDkTerY5kkl04E8K2klO4A8tyWPY5L6ddF2fZ1V4syyG/oAo7nk/3ixv9JynCR8tCXNIc51PAGZ+ZjQrcm+fjz7lekiM8OeCwKQpVTYrIF4Encd1j71XVd0Tk6vT+O4HHcF1jN+C6x35quOQ1GAyjk5uOauS/d8ZxgLaxFo/euW95zueNc9NRjXRvjs+MhAj35/l5NOzljLYkT4c87LGHc1Gn/wxrwJ2qPoarDDpuu7PDawW+cKDlMhgMhw6lHj8X5o/tdf9g6Zgev962+Euub7/9IxkTmW0wGA5pZgZCQx7Mli2j+3BxcMx7DAaDwTBsGEVhMBgMhl4xisJgMBgMvWIUhcFgMBh6xSgKg8FgMPSKURQGg8Fg6BWjKAwGg8HQK0ZRGAwGg6FXjKIwGAwGQ68YRWEwGAyGXjGKwmAwGAy9Im7evdGDiDQD64ZbjmGgGKgbbiGGAXPehxbmvIeOyarabeW30ZgUcJ2qVgy3EAcaEVlmzvvQwZz3ocVwn7dZejIYDAZDrxhFYTAYDIZeGY2K4q7hFmCYMOd9aGHO+9BiWM971BmzDQaDwZBdRuOMwmAwGAxZxCgKg8FgMPTKqFIUInKWiKwTkQ0ics1wy3OgEJEtIrJaRFaKyLLhlmeoEJF7RWSXiLzdYVuRiDwtIu+lnwuHU8ahoIfzvl5EqtLf+UoROWc4ZRwKRGSiiDwvImtF5B0R+XJ6+6j+zns572H7zkeNjUJEbGA98H6gElgKXKqqa4ZVsAOAiGwBKlR1VAciicgpQAvwe1U9PL3tRqBeVX+avjkoVNVvDaec2aaH874eaFHVm4ZTtqFERMqAMlVdISK5wHLgQuAKRvF33st5f4Rh+s5H04ziGGCDqm5S1TjwIHDBMMtkyCKq+hJQ32XzBcB96df34f6hRhU9nPeoR1WrVXVF+nUzsBYoZ5R/572c97AxmhRFObC9w/tKhvnDPYAo8JSILBeRq4ZbmANMqapWg/sHA8YOszwHki+KyKr00tSoWn7piohMARYB/+YQ+s67nDcM03c+mhSFdLNtdKyr9c2JqnokcDbwhfRShWF0cwcwHTgCqAZ+MazSDCEikgM8BHxFVZuGW54DRTfnPWzf+WhSFJXAxA7vJwA7hkmWA4qq7kg/7wL+jrsMd6hQk17TbV/b3TXM8hwQVLVGVVOq6gC/YZR+5yLixb1Y/lFVH05vHvXfeXfnPZzf+WhSFEuBmSIyVUR8wCXAI8Ms05AjIuG0wQsRCQNnAm/3ftSo4hHgk+nXnwT+MYyyHDDaL5RpLmIUfuciIsA9wFpVvbnDrlH9nfd03sP5nY8aryeAtLvYLYAN3KuqPx5eiYYeEZmGO4sANxvwA6P1vEXkT8D7cFMu1wDfB/4P+AswCdgGfFhVR5Xht4fzfh/uEoQCW4DPtq/bjxZE5CRgCbAacNKbv4O7Xj9qv/NezvtShuk7H1WKwmAwGAzZZzQtPRkMBoNhCDCKwmAwGAy9YhSFwWAwGHrFKAqDwWAw9IpRFAaDwWDoFaMoDAaDwdArRlEYDAaDoVeMojAYDjAi8jsR+ddwy9GOiLwgIrcO4LhCEakRkemDHP9vIvK1wfRhGFqMojBkFREpEZHb08WUYukLybMi8v7hlu1QZ6AKoRe+AzymqhsH2c8PgO+JSH4WZDIMAZ7hFsAw6ngICAFXAhtwU0CfCowZTqEM2UVEQsB/AB8cbF+qulpENgGXA7cNtj9D9jEzCkPWEJEC4GTgGlV9VlW3qupSVb1JVR9MtxER+aaIbBSRSLqE6+Vd+gmll2da0jOS74jIv0Tkdx3a7Hd33HVJp6+x0n3cLiL/LSJ14pYbvUlErC59fD1ddjMmIpUi8pP+jtHPz60/n0mvsqaTQ/6+w2f27Y6fWfr5VNw09Jp+TOkwhNXb59AN5+DmIXqli5wT0n1/VESeE5E2EXlLRGaLSIWIvJTe9oaITOpw6CO4uYwMIxFVNQ/zyMoDd4baDPwKCPTQ5sfAOuAsYCrwMaAVOLdDm9uBKuADwOHAX4Em4Hcd2rwA3Nql798B/+rvWOk+GoEfArNwS00mcUvotvfxE6AB+DQwAzge+Hwm59PNZ5CRnP2RFbgT2IpbCngeboXHxvbPDMgHXgXuBcalH3Z/P4duzuF/gKe62X4ebtK6F3AV03zgXeB14DngJGAhsAm4pcNxZwFxIDjcv2Pz6Ob7Hm4BzGN0PYAP4ZbtjAKvATcBx6b3hYEIcHKXY27BXesGyAFiwGUd9uekL9a/67DtBXpRFP0c6wXgtS77nwbu7jBuFLi6h3Ptc4wejstIzr5kTcsZBy7pItuevj6z/nwOPZzD/wH3dbP9u+nvqrTDtv8FaoExHbb9Fvhzh/cLcBXM9OH+DZvH/g9jozBkFVV9SEQexV2COh73TvHrIvJd3ItPAHhCRDqmLfbipk0Gt4KXD1fJtPfZIiKrMxRlbj/GAljV5bgd7CutORfwA88OcoxsyNmbrNPT7d9o36GqrSKSSb2C3j6H7gjipjzvyhG4SrDjvknAw6q6u8u2pR3eRzr0axhhGEVhyDqqGsVVCk8DPxSRu4Hr2XfB/SBuHYGOJNLP3ZW07Q6nm7beDq/b19d7G6vra3DvatuP7UuW/o6RrT56klU6vB8ovX0O3VEHdFezeSHu0mNHFgHXddPurg7vi9LPtb2LaRgOjKIwHAjW4P7W3sVdVpqsqs/10HYD7kXrONx17PbKfYcDHd0wa4GyLscuZN9d+Jp+jNUfuWPA6cB7vezPxhiD6aP9MzsG2Ax7vZK6fmZx3KJe2eBN4IqOG9Lf0/T0vvZtRbglijtum4jrBfdmh8MPB3Z0mYkYRghGURiyhoiMwTU834u7lNEMVADfBJ5V1UYRuQm4SUQEeAl3ff04wFHVu9LLTPcAPxORWtwlkOvY/wL3HHCLiJyPawj+LO4FaQuAqjb3NVZf55Pu43+An4hILN3HGOAoVb0ji2MMto8WEbkX9zOrA6qB7+HOCDrOMrYAx6S9nVqAenXrLw+EJ9PjjemwpLQg/fxWh3aLcBXhmg7bjkiPv6HDtpOBJwYoi2GIMYrCkE1acL1bvozrIeTH9V56APhRus21uGvb3wDuwPVmWgnc2KGfb+AaY/8OtOEaQ8NdxroX98J0b/r97en2xR3a9Gesvvg2rlH4WmBCur/fZ3mMbPTR/pk9gvs9/BIoxTXGt3MTcB/uRTuI62G1JYMx9qJu7MMbuLXp22MfFgLvqWpLh6aLgLdVtePS1kLgrXYlJSIB3BrQHxiILIahx5RCNRwUpOMj6lT1iuGW5WBARPy47rI/V9VfDNEYZ+G6yc5V1dQg+vkCcIGqnpk14QxZxcwoDIZRgIgsAubgej7lAt9KP/95qMZU1SdE5DbcmdbWQXSVAP4zO1IZhgKjKAyG0cPXgMNwg+VWAqeoauVQDqiqXT2cBtJHn3YYw/Bilp4MBoPB0Csm15PBYDAYesUoCoPBYDD0ilEUBoPBYOgVoygMBoPB0CtGURgMBoOhV4yiMBgMBkOvGEVhMBgMhl75/+B6x4lWLYekAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "CB91_Blue = '#2CBDFE'\n",
    "CB91_Green = '#47DBCD'\n",
    "CB91_Pink = '#F3A0F2'\n",
    "CB91_Purple = '#9D2EC5'\n",
    "CB91_Violet = '#661D98'\n",
    "CB91_Amber = '#F5B14C'\n",
    "\n",
    "\n",
    "ax = plt.subplot()\n",
    "colors = ['#F5B14C','#2CBDFE', '#47DBCD', '#F3A0F2','#9D2EC5','#661D98']\n",
    "markers = [\"o\", \"s\", \"^\", \"h\",\"d\" ]\n",
    "lines = [\"--\", \"--\"]\n",
    "offset = np.array([0,0.5,0.2,-0.2,-0.3])+0.25\n",
    "x_points = np.linspace(0, sequence_lengths[-1] + 1, 100)\n",
    "for k in [0, 1, 2,3,4]:\n",
    "\n",
    "    for fu in [0]:\n",
    "        ax.scatter(\n",
    "            sequence_lengths,\n",
    "            average[k],\n",
    "            marker=markers[k],\n",
    "            color=colors[k],\n",
    "            s=50,\n",
    "            label=f'k={k}'\n",
    "            \n",
    "        )\n",
    "        ax.plot(\n",
    "            x_points,\n",
    "            [\n",
    "                np.real(exp(m, *poles[(fu, k)], *amplitudes[(fu, k)]))\n",
    "                for m in x_points\n",
    "            ],\n",
    "            ls=lines[fu],\n",
    "            color=colors[k] \n",
    "        )\n",
    "\n",
    "ax.legend(loc=(0.74, 0.35), fontsize=12)\n",
    "ax.set_xlabel(\"Sequence length ($m$)\", fontsize=14)\n",
    "ax.set_ylabel(\"Weighted average $f_k(m)$\", fontsize=14)\n",
    "ax.set_ylim([-0.1,1.1])\n",
    "ax.set_xlim([0,27])\n",
    "\n",
    "plt.savefig(\"experiment-final.png\", dpi=1200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error bars by bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run bootstrap resampling on corr_data to obtain error bars on the decay rats of each fit. (this takes a while to run)\n",
    "\n",
    "pole_dist = {}\n",
    "NUM_RUNS = 10000\n",
    "for k in [0,1,2,3,4]:\n",
    "    pole_dist[k] = np.array(\n",
    "        [\n",
    "            (\n",
    "                curve_fit(\n",
    "                    exp,\n",
    "                    sequence_lengths,\n",
    "                    np.array(\n",
    "                        [\n",
    "                            np.real(\n",
    "                                np.mean(\n",
    "                                    resample(\n",
    "                                        corr_data[k][m], replace=True, n_samples=200\n",
    "                                    )\n",
    "                                )\n",
    "                            )\n",
    "                            for m in range(len(sequence_lengths))\n",
    "                        ]\n",
    "                    ),\n",
    "                    p0=inits[k],\n",
    "                    maxfev=1000000,\n",
    "                )\n",
    "            )[0][:1]\n",
    "            for i in range(NUM_RUNS - 1)\n",
    "        ]\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD6CAYAAABUHLtmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAswUlEQVR4nO3de3TU5b3v8feXICiXCEECCSFGSyqBgBGiwNksDy07lqNuKEJTkNYgKLvUtlTr6YkNq61/KLGnbvG23Zs2bdHuAwe6t8JGRGO0y3WomI0SMKA0KrQkhpuBcr8kPOeP+SVMkpkQksnMZObzWisrM8/8Ls88mcz393uu5pxDREQkkB6RzoCIiEQvBQkREQlKQUJERIJSkBARkaAUJEREJCgFCRERCUpBIkzM7EozKzez7Wa208we9dKTzKzUzKq83wP99nnEzD4xs91m9rXI5V5E4pV113ES11xzjcvIyIh0NtrNOceFCxdISEjAOcfHH3/M8OHDOXr0KD179mTo0KHs37+f+vp60tLSOH36NHv27OH06dOHgVuAN4EvO+cagp2ju5VJR73//vuHnXOD27t9PJSLyiSwyykXlUlgPbsyM10pIyODrVu3RjobHXLq1CkmT57MCy+8wD333MMf//hHUlJSqK2tZcqUKWzdupVly5YB8JOf/OQvzrk9ZvYJvmDxbrDjducyuRxm9pfL2T4eykVlEtjllIvKJDBVN4VRQ0MDOTk5JCcnk5eXx4QJEzhw4AApKSkApKSkcPDgQQBqamoYPny4/+7VwLCwZ1rC7qmnnmL06NFkZ2czd+5czpw5Q11dHXl5eWRmZpKXl8eRI0eatvcuKLJVLSldQUEijBISEqioqKC6upry8nIqKyuDbhukGrBVopktMrOtZrb10KFDIcytREJNTQ3PPPMMW7dupbKykoaGBlavXk1xcTFTp06lqqqKqVOnUlxcDMCuXbtYvXo1wE5gGvDPZpYQwbcgMUZBIgIGDBjAlClT2LRpE0OGDKG2thaA2tpakpOTAUhLS2Pfvn3+u6UBn7c8lnNuhXMu1zmXO3hwu6sZJYrV19dz+vRp6uvrOXXqFKmpqaxbt46CggIACgoKeOWVVwBYt24dc+bMAXDOuT1AY7WkSEgoSITJoUOHOHr0KACnT5/mzTffZOTIkUyfPp2VK1cCsHLlSmbMmAHA9OnTG68QzcyuAzKB8kjkXcJn2LBhPPzww6Snp5OSksLVV1/NbbfdFpJqSd11SkcoSIRJbW0tX/nKVxg7diw333wzeXl53HnnnRQWFlJaWkpmZialpaUUFhYCMHr0aPLz8wFGA5uAB9rq2SSx4ciRI6xbt449e/bw+eefc/LkSX7/+98H3b691ZLetrrrlMt2yd5NZvYb4E7goHMu20tLAv4vkAHsBfKdc0e81x4BFgINwA+cc6976eOB3wFXARuBJc45Z2a9gReB8cAXwDedc3tD9g6jxNixY9m2bVur9EGDBlFWVhZwn6KiIpYuXVrpnMvt6vxJdHjzzTe57rrraPwSv+uuu/jTn/7UVC3Z2AvucqslRTqqPXcSv8PXIOavEChzzmUCZd5zzGwUMAff1W/LRrQXgEX4qk0y/Y65EDjinBsBPAU80dE3I9Ldpaens2XLFk6dOoVzjrKyMrKyslQtKRFzySDhnHsHqGuRPANY6T1eCXzdL321c+6sfyOamaUAic65d53v/vjFFvs0HusPwFQzs469HZHubcKECcyePZtx48YxZswYLly4wKJFi1QtKRHT0cF0Q5xztQDOuVozS/bShwFb/LZrbEQ77z1umd64zz7vWPVm9jdgEHC4g3kT6dYeffRRHn300WZpvXv3VrWkRESoG64D3QG4NtLb2qf1wdU7I6pkFL4a6SzEFJVn+6mswqejQeKAV4WE9/ugl14N+PfHa2xEq/Yet0xvto+Z9QSupnX1FqDeGSIi4dbRILEeKPAeFwDr/NLnmFlv/0Y0r2rquJlN9Nob7mmxT+OxZgNvue4666CISIxpTxfYVcAU4BozqwZ+BhQDa8xsIfBX4BsAzrmdZrYG2AXU07wRbTEXu8C+5v0AlAAveRPY1eHrHSUiIlHgkkHCOTc3yEtTg2z/GPBYgPStQHaA9DN4QUYkXqmOXaKVRlzLZdGXmUh8UZAQEZGgFCRERCQoBQkREQlKQUJERIJSkBARkaAUJEREJCgFCRERCUpBQiLq6NGjzJ49m5EjR5KVlcW7775LXV0deXl5ZGZmkpeXx5EjR5q2X7ZsGUC2me02s69FLONdRONQJNooSEhELVmyhGnTpvHxxx+zfft2srKyKC4uZurUqVRVVTF16lSKi4sB2LVrV+MCOztpvaiVxIAFCxaQnJxMdvbFyRkuddEwYsQIbrjhBl5//fWm9Pfff58xY8aA74LiGa1R03EKEhIxx44d45133mHhwoUA9OrViwEDBrBu3ToKCnxzPhYUFPDKK68AsG7dOubMmQPg/Be1ikTepWvMnz+fTZs2NUu71EXDzp072bRpE9/97ndpaPBNFbd48WJWrFgBUEnzlTDlMilISMR89tlnDB48mHvvvZebbrqJ++67j5MnT3LgwAFSUlIASElJ4eBB30z0NTU1DB/uPxN9s8WrJAbceuutJCUlNUu71EVD7969ue666xgxYgTl5eXU1tZy7NgxJk2a1HgI/5Uw5TIpSEjE1NfX88EHH7B48WK2bdtG3759m64SAwkyg3zMLVC1e/ducnJymn4SExNZvnx5m9UuwFAz+yQW22qCXTQs+8PmZhcNaWlp1NTUUFNTQ1qa//I1upjoDAUJiZi0tDTS0tKYMGECALNnz+aDDz5gyJAh1NbWAlBbW0tycnLT9vv27Wt2CC4uXtVMd1qgqmVj9Q033EBFRQUVFRW8//779OnTh5kzZ7ZZ7QIk4VvnOubbao6dPg9AoGsGM4ubi4lwUZCQDglFL5yhQ4cyfPhwdu/eDUBZWRmjRo1i+vTprFy5EoCVK1cyY8YMAKZPn97YcG3+i1p1OiNRrKysjC996Utce+21bVa7AHXOubOx2FbT8qKhR98BAPTsP6jZRUN1dTWpqamkpaVRXV3tf4iYuJiIFAUJiahnn32WefPmMXbsWCoqKvjJT35CYWEhpaWlZGZmUlpaSmFhIQCjR48mPz8ffFfMm2i+qFVMWr16NXPn+pZ0aautBjjnt1tMVa+0vGjoM8J353nViAmsXr2as2fPsmfPHqqqqrjllltISUmhf//+bNmypfEQ/ithymW65KJDIl0pJyeHrVu3tkovKysLuH1RURFLly6tdM7ldnXeIu3cuXOsX7++cWxIUO2tXjGzRcAigPT09FBkMeTmzp3LH//4Rw4fPkxaWhqPPvoohYWF5OfnU1JSQnp6OokT7weg1+Bryc/PZ9SoUfTs2ZPnn3+ehARfLdsLL7zA/PnzwbfQ2a+5uBKmXCYFCZEo9dprrzFu3DiGDBkCXKx2SUlJadVWA/Ty2zVg9YpzbgWwAiA3Nzcq15FftWpVwHT/iwb/qs6ioiKKiopabZ+bm0tlZSVmVumc+17ocxo/VN0kIaGRwqG3atWqpqomaF3t4t9WAySZWe94aauR8FGQCJN9+/bxla98haysLEaPHs3TTz8NwM9//nOGDRvW1N1x48aNTfvE+hQUEtypU6coLS3lrrvuakprq60GqAN2ESdtNfEiGi6+VN0UJj179uTJJ59k3LhxHD9+nPHjx5OXlwfAgw8+yMMPP9xs+xZTUMwG3jSzL+ufPz706dOHL774olnaoEGDgrbVAPvjoZ1Gwk93EmGSkpLCuHHjAOjfvz9ZWVmNvVIC0hQUIhINFCQiYO/evWzbtq1pENlzzz3H2LFjWbBgQdMoWk1BISLRQEEizE6cOMGsWbNYvnw5iYmJLF68mE8//ZSKigpSUlL40Y9+BFxet0aNGBWJPdHQHgEKEmF1/vx5Zs2axbx585oaJIcMGUJCQgI9evTg/vvvp7zc1ymlvVNQRGLEaLR8eEWk6ylIhIlzjoULF5KVlcVDDz3UlN443QDAyy+/3DSPfjxOQSEi0Ue9m8Jk8+bNvPTSS4wZM4acnBwAHn/8cVatWkVFRQVmRkZGBv/6r/8KXJyCYseOHXEzBYWIRB8FiTCZPHlywHaG22+/Peg+8TQFhYhEJ1U3iYhIUAoSIiISlIKEdIp6OonENgUJEREJqlNBwsweNLOdZlZpZqvM7EozSzKzUjOr8n4P9Nv+kUDr8JrZeDP70HvtGTOzzuRLRERCo8NBwsyGAT8Acp1z2UACMAcoBMqcc5lAmfccMxvlvR5oHd4X8C2Gkun9TOtovkREJHQ6W93UE7jKzHoCffCNCJ4BrPReXwl83Xs8A1jdch1eM0sBEp1z7zpfH9EX/fYREZEI6nCQcM7VAL8E/grUAn9zzr0BDHHO1Xrb1ALJ3i7DAP95JhonrBvmPW6ZLhLz1PAv0a4z1U0D8d0dXAekAn3N7Ftt7RIgzbWRHuicmswuAvRFJhK/OlPd9PfAHufcIefceeA/gP8GHPCqkPB+H/S2rwb8575unLCu2nvcMr2VSExmJyISzzoTJP4KTDSzPl5vpKnAR8B6oMDbpgBY5z1eD8xpuQ6vVyV13Mwmese5x28fERGJoM60SbwH/AH4APjQO9YKoBjIM7MqIM97jnNuJ7CGwOvwLgZ+ja8x+1PgtY7mSyQWHD16lNmzZzNy5EiysrJ49913qaurIy8vj8zMTPLy8poWqPIMDdS9PNao6jP8OtW7yTn3M+fcSOdctnPu217PpS+cc1Odc5ne7zq/7R9zzn3JOXeDc+41v/St3jG+5Jz7nguy4o7EnoyMjKaZcXNzffMYtvVluGzZMoDsWP8yXLJkCdOmTePjjz9m+/btZGVlUVxczNSpU6mqqmLq1KkUFxcDvvXQgSQCdy8X6RSNuJZ2C3YV19mru7fffpuKigq2bt0K0OaXobfGxk5i+Msw/cG1vPPOOyxcuBCAXr16MWDAANatW0dBga8mt6CggFdeeQXwrYcO1LXsXh6JvEvsUZCQqNPWl+GcOXMAXCx/GdYf3c/gwYO59957uemmm7jvvvs4efIkBw4cICUlBYCUlBQOHvT1CampqQE453cIdSOXkFGQkDZ1dR2wmXHbbbcxfvx4VqxYAdDml+Hw4f4d5GLzy9BdaOCDDz5g8eLFbNu2jb59+zbdTQXcXuuhSxfSokMSUZs3byY1NZWDBw+Sl5fHyJEjg27b3i9D8H0h4pvqhfT09FBkNWx69r+GtLQ0JkyYAMDs2bMpLi5myJAh1NbWkpKSQm1tLcnJvnGqaWlpAL38DhF0PXR8nUvIzc1Vu5+0i+4kpF266o4iNTUVgOTkZGbOnEl5eXnTlyHQ6stw3z7/QfuxOaYmod9Ahg8fzu7duwEoKytj1KhRTJ8+nZUrfTPerFy5khkzZgC+9dCBpJbdyyOSeYk5upOQiDl58iQXLlygf//+nDx5kjfeeIOf/vSnTV+GhYWFrb4M7777bgCL9S/DZ599lnnz5nHu3Dmuv/56fvvb33LhwgXy8/MpKSkhPT2dtWvXAr710IE6fN3L69F66BJCChISMQcOHGDmzJkA1NfXc/fddzNt2jRuvvnmoF+G+fn57NixYzStx9rElJycnKbeXv7KysqC7bI/1tdCf+qpp/j818vB4IrBGVxz+w+pq6vjm9/8Jnv37iUjI4M1a9YwcKBvdQL/7tLAD5xzr0cw+92WqpskYq6//nq2b9/O9u3b2blzJ0VFRQAMGjSIsrIyqqqqKCsrIykpqWkfb5vKlmNtJLbV1NTwzDPPMLTgKVIX/jNcuMDJj96J++7S4aAgISLdQn19Pa7+HO5CA67+LAn9kniq5P/EdXfpcFCQEJGoN2zYMB5++GFqXriX6ue+jfXuw1XXjaPh5NG47i4dDgoSIhL1jhw5wrp16xj2nRLSHngRd/4sJ3a+HXR7jR0JHQUJEYl6b775Jtdddx0Jfa7GEnrS58uTOFvzEQl9B3Squ3R37iodLgoSIhL10tPT2bJlCxfOn8E5x5m/bOeKQcPpM2JC0LEjXsN1zHeX7moKEhJSmspZusKECROYPXs2tb/7IbW/eQCco/+N00icOJvS0lIyMzMpLS2lsLAQuNhdGt/MuDHdXbqraZxEmOzbt4977rmH/fv306NHDxYtWsSSJUvUz1uknR599FFWnm3eQSmh5xVBx44UFRWxdOnSylgfP9LVdCcRJj179uTJJ5/ko48+YsuWLTz//PPs2rVL/bxFJKopSIRJSkoK48aNA6B///5kZWVRU1MT99Nii0h0U5CIgL1797Jt2zYmTJjQ6Wmx1YVPJPZEU9uegkSYnThxglmzZrF8+XISExODbtfeft7qwiciXUlBIozOnz/PrFmzmDdvHnfddRdASKbFFhHpKgoSYeKcY+HChWRlZfHQQw81pbe1RoD6eYtIpClIhMnmzZt56aWXeOutt8jJySEnJ4eNGzdSWFioft4iErU0TiJMJk+eHKydIWb6eUdTY5uIhIbuJEREJCgFCQk53VF0XkZGBmPGjCEnJ4fcXN+NZF1dHXl5eWRmZpKXl8eRI0f8dxlqZp+Y2W4z+1pEMi0xSUFCgtKXfWS9/fbbVFRUNC1j2tbofCAJX/uVRudLSClIiHQTbY3OB+qcc2c1Ol9CTUFCJAqZGbfddhvjx49nxYoVAG2OzgfO+e2uVdgkZNS7SSQKbd68mdTUVA4ePEheXh4jR44Muu3lrMIGLALf+gzSPWQUvsre4jsidn7dSYhEodTUVACSk5OZOXMm5eXlbY7OB3r57a5V2CRkFCREosyFc2c4fvw4ACdPnuSNN94gOzu7zdH5QJKZ9dbofAk1VTdJxDU0NJCbm8uwYcPYsGFD3C/E1HDqKJMnTwagvr6eu+++m2nTpnHzzTeTn59PSUkJ6enprF27FvCNzgfqgF1APRqdLyGkICER9/TTT5OVlcWxY8eAi109CwsLKS4upri4mCeeeKLlQkyzgTfN7Mux9oV4xYChbN++vVX6oEGDgo7OB/Z3l5H50r10qrrJzAaY2R/M7GMz+8jMJplZkpmVmlmV93ug3/aPBBrwY2bjzexD77VnzMw6ky/pPqqrq3n11Ve57777mtK0EJNI9Ohsm8TTwCbn3EjgRuAjoBAoc85lAmXec8xsFDCHwAN+XsDX6yLT+5nWyXxJN/HDH/6QX/ziF/TocfGj2NmFmECLMYmESoeDhJklArcCJQDOuXPOuaPADGClt9lK4Ove4xnA6pYDfswsBUh0zr3rfH35XvTbR2LYhg0bSE5OZvz48e3avr1dPb1t1ZNHJAQ60yZxPXAI+K2Z3Qi8DywBhjjnagGcc7VmluxtPwzY4rd/41Xgee9xy3SJcZs3b2b9+vVs3LiRM2fOcOzYMb71rW81dfVMSUnRQkwiEdaZ6qaewDjgBefcTcBJvKqlIAK1M7g20lsfQFUIMWXZsmVUV1ezd+9eVq9ezVe/+lV+//vfayEmkSjSmSBRDVQ7597znv8BX9A44FUh4f0+6Le9f4Vy41Vgtfe4ZXorqkIIn0hO7qeFmESiR4eDhHNuP7DPzG7wkqbi66e9Hijw0gqAdd7j9cCclgN+vKqp42Y20evVdI/fPhInpkyZwoYNG4CLXT2rqqooKysjKSmpabuioiKASufcDc651yKTW5H40dlxEt8H/s3MegGfAffiCzxrzGwh8FfgGwDOuZ1mtobAA34WA78DrgJe835EYpqmYpfuoFNBwjlXAQQawDM1yPaPAY8FSN8KZHcmLyIiEnqau0lERIJSkBARkaA0d5O0orpyiUZHjx7l0MuPc+7wXwG45vYl9ExKIy8vL24ngwwH3UmISLewZMkSrrx+PMPu/xdSFzzLFYOGc2zL2qDrfvtNBql1vztBQUJEot6xY8d455136Df2NgAs4Qp6XNmPU5+8p8kgu5iqm0Qk6n322WcMHjyYXRuXc+7gHnoPHcHAqYtoOHk06GSQEydO9D9EwOl+tKTrpelOQkSiXn19PR988AH9b7qd1Hufwa7ozbEta4Nu397JIDWLw6UpSIhI1EtLSyMtLY3eqb4JHvrc8HecO/ApCX0HBF33W5NBhoaCRJgsWLCA5ORksrMvjhn8+c9/zrBhw8jJySEnJ4eNGzc2vbZs2TJGjBgBvt4ZX2t9RJH4MXToUIYPH875L3wTRp/5y3auuCadPiMmaDLILqYgESbz589n06ZNrdIffPBBKioqqKio4Pbbbwcu9szYuXMnwJ9Rz4y409DQwE033cSdd94JQF1dHXl5eWRmZpKXl8eRI0eatvXv6hnLFxTPPvsshzf8ks9/8z3OHdxD4qR8EifO1mSQXUxBIkxuvfXWZhPVtaWxZ0bv3r0BzqGeGXGncd3vRo3rfsdzV8+cnBxSCpaTuuA5ku9aSsKV/Ui4KlGTQXYxBYkIe+655xg7diwLFixoujqMhWU6NSCv4+qPHda639JMJP+fFCQiaPHixXz66adUVFSQkpLCj370I0DLdMa7I2UrumTdb5GOUJCIoCFDhpCQkECPHj24//77KS/3taupZ0b8OvVJOT36DuiSdb+j9a5TopuCRAQ1dt0DePnll5t6PjX2zDh79ixAL9QzI26crdnF6ar3yMjIYM6cObz11lvN1v2Gjnf11F2ndISCRJjMnTuXSZMmsXv3btLS0igpKeHHP/4xY8aMYezYsbz99ts89dRTwMWeGaNGjQL4MuqZETcG/vf5pD2wUut+S9TQtBxhsmrVqlZpCxcuDLp9UVERRUVFmFmlembEn4zCV/ndtL5NzwsLC8nPz6ekpIT09HTWrvWNNm68oNixY4e6ekqXUJAQiVJTpkxhypQpwMV1vwMpKipi6dKllc65QKtEinSKqptEpFtTd+uupSAhIiJBKUhIxJw5c4ZbbrmFG2+8kdGjR/Ozn/0M0BQUEpjuGCJDQUIipnfv3rz11lts376diooKNm3axJYtWzQFhUgUUZCQiDEz+vXrB8D58+c5f/48ZqYpKESiiIKERFRDQwM5OTkkJyeTl5fHhAkTNAWFSBRRkJBmwl3vm5CQQEVFBdXV1ZSXl1NZWRl0W01BIRJ+ChISFQYMGMCUKVPYtGmTpqAQiSIKEhIxhw4d4ujRowCcPn2aN998k5EjR2oKCpEoohHXEjG1tbUUFBTQ0NDAhQsXyM/P584772TSpEmagkIkSihISMSMHTuWbdu2tUrXFBQi0UPVTSIiEpSChIiIBKUgISJRL56m5Ii299rpIGFmCWa2zcw2eM+TzKzUzKq83wP9tn3EzD5pOe+OmY03sw+9154xM+tsvkREpPNCcSexBPjI73khUOacywTKvOeY2ShgDjCa1vPuvAAswtelMdN7XSRmRdvVokgwnQoSZpYG3AH82i95BrDSe7wS+Lpf+mrn3Fn/eXfMLAVIdM6963xDal/020dERCKos3cSy4EfAxf80oY452oBvN/JXvowwH+4bOO8O8O8xy3TRUQkwjocJMzsTuCgc+799u4SIM21kR7onJqPR2Kaqz9H7YsP8vlvvqc1NiQqdOZO4u+A6Wa2F1gNfNXMfg8c8KqQ8H4f9LavBvyn8Gycd6fae9wyvRXNxyMxL+EKhsx5nNQFz2mNDYkKHQ4SzrlHnHNpzrkMfA3SbznnvgWsBwq8zQqAdd7j9cAcM+vtP++OVyV13Mwmer2a7vHbRySumBk9el0FaI0NiQ5dMU6iGMgzsyogz3uOc24nsAbYRet5dxbja/z+BPgUeK0L8iXSLbgLDXz+2++HfI0NVddKR4QkSDjn/uicu9N7/IVzbqpzLtP7Xee33WPOuS85525wzr3ml77VOZftvfY9F2ThAJF4YD0SSL332ZCvsdHdq2t9wfMHHPzDowA0nD7OgdVLqVlxf8B2mhEjRoCvrUbtNJ2gEdfSZTQWoHNCvcZGd3d863quGHTxrunYlrVcmXEjwxb9KmA7zc6dOwH+jNppOkVBQiSKNJz6GxfOnAC0xoa/+mOHOf3Zf9Hvxtua0k598h59s6cCgdtpevfuDXAOtdN0iqYKF4kiDSfqOPzqU+AucPOGvlpjw3OkbAUDpizAnTvVlNZw8ig9+yUBrdtpJk6c6L97m+00+GZ7ID09vUvy3t3pTiJMFixYQHJyMtnZ2U1pl+r7Hgt1qqpyujy9kq8j9d5nSF3wHJWVlfz0pz8FLq6xUVVVRVlZGUlJSU37FBUVAVS2bOuLFRs2bKBH3wH0HjqiXdvHUztNOChIhMn8+fPZtGlTs7RL9X1XnaoIbN68mdNV71H9wgIOrf8FZ/6yg8P/+UsS+g6g/oSvX0y8ttOEg4JEmNx6663Nrv6AS/Z9D3edqq76JRotW7aMtAdWkrb4Nwye/mOuvHYs1/zDw/QZMYGTlb4VDAO105w9exagFzHaThMuapOIoLb6vqtOVaRtiRNnc3hdMSd2vEHpTVmt2mlGjRoF8GVgZiy204SLgkQUutw6VWAFQG5ursaXSEy7Mn0sV6aPBSDhqkSGzHkcgLLiO5ptV1RURFFREWZWGYvtNOGk6qYIUt/3+KRqPelOFCQi6FJ931WnKiKRpiARJnPnzmXSpEns3r2btLQ0SkpKKCwspLS0lMzMTEpLSyksLAQC1ql2ed93Xd2KSCBqkwiTVatWBUwvKysLmB4Pdar79u3jnnvuYf/+/fTo0YNFixaxZMkS6urq+OY3v8nevXvJyMhgzZo1DBzoWyrdf+0E4AfOudcj+R5EYp3uJCRievbsyZNPPslHH33Eli1beP7559m1a5fWThCJIgoS0uWCVWWlpKQwbtw4APr3709WVhY1NTVaO0EkiihISFTYu3cv27ZtC9naCSISGgoSEnEnTpxg1qxZLF++nMTExKDbXc74ES2wIxIaChISUefPn2fWrFnMmzePu+66CwjN+BFN3CYSGgoSErHur845Fi5cSFZWFg899FBTeryvndBI3ZIlGihISMRs3ryZl156ibfeeoucnBxycnLYuHHjJcePADG9doJIIJG6aNA4CYmYyZMnB2tnaHP8yNKlSyudc7ldmTcR8dGdhIiIBKUgIRJl6o8dYv+qR6j51XcYPXo0Tz/9NND2SobAUDP7xMx2d+eVDCX6KEiIRJseCQz8ykKG3f8v7R6JDiTha6vRSHQJKQWJOKceNNGnZ7+kpvWc2zsSHahzzp3VSHQJNQUJkSjW3pHo+Ja5baSR6BIyChIiUerCudMhHYmuUejSEQoSIlHINdRz6OXH2z0SHd/iVI0CjkTvrqPQVSUaWQoSIlHGOccXrz3NFYOGt3skOpBkZr3jYSR6IAokXUeD6USizNmaXZzc+TZXDM4gJycHgMcff5zCwkLy8/MpKSkhPT2dtWvXAr6R6EAdsAuoRyPRJYQUJESizJVpo7n2f20AoKL4jmavBRuJDuzXKHTpCqpuEhGRoBQkREQkKAUJEREJqsNBwsyGm9nbZvaRme00syVeepKZlZpZlfd7oN8+jwSaX8bMxpvZh95rz5iZde5tSXuoR4iIXEpn7iTqgR8557KAicADZjYKKATKnHOZQJn3HO+1OQSeX+YFYBG+rnuZ3usiIhJhHQ4Szrla59wH3uPjwEf4pgKYAaz0NlsJfN17PANY3XJ+GTNLARKdc+8639DRF/32ERGRCApJF1gzywBuAt4DhjjnasEXSMws2dtsGLDFb7fG+WXOe49bpgc6zyJ8dxykp6eHIusi0g3UHzvE4Vf/iYYTRzDrQb+cr5GYO4OG08c5vO4J6o8dIO/90axZs4aBA3013MuWLQPINrPdwA+cc69H8j10V51uuDazfsC/Az90zh1ra9MAaa6N9NaJ3XRagWik9gjpVvymTx/67V9y/INXOXf4rxzbspYrM25k2KJftZo+3VsPfSfdaPr0aPy/7FSQMLMr8AWIf3PO/YeXfMCrQsL7fdBLrwaG++3eOL9Mtfe4ZXrcyMjIYMyYMeTk5JCb6xsP1bjADL4roWYdALqjaPzwR4LKoWP8p0/v0bsPVwwaTsPxLzj1yXv0zZ4KtJ4+fc6cOQBO06d3Tmd6NxlQAnzknPsnv5fWAwXe4wJgnV/6nJbzy3hVU8fNbKJ3zHv89okbb7/9NhUVFWzduhWgaYEZoBK/DgAi8a7+bwc4d+AzeqfeQMPJo/TslwS0nj59+HD/a1JNn95RnbmT+Dvg28BXzazC+7kdKAbyzKwKyPOe45zbCazBN7/MJprPL7MY+DW+aP8p8Fon8hUT/BeYoXkHAJG4deHcaQ69/DhJU++nR+8+QbfT9Omh0+GGa+fc/yNwewLA1CD7PAY8FiB9K5Dd0bx0d2bGbbfdhpnxj//4jyxatKjZAjMtOgB0mqo8pDtqnD6976gp9LnhvwGQ0HcA9Sfq6NkvqdX06fv27fPfPej06cAKgNzc3ICRJd5pgr8osHnzZlJTUzl48CB5eXmMHDmy3fuqx5fEA//p0xNvmdmU3mfEBE5WlnH1xG+0mj797rvvBl/NeFxOnx4qmpYjCqSmpgKQnJzMzJkzKS8vb7bATIsOAM109x5fCxYsIDk5mezsizeSjY32mZmZ5OXlceTIkabX/Ls1+o/al9jWOH36mb/u4PPffp/Pf/t9Tn/6XyROnM2ZvduoWXE/paWlFBb6mu5Gjx5Nfn4++AbvtqzelsugIBFhJ0+e5Pjx402P33jjDbKzs5stMEPzDgAxZf78+WzatKlZWmOjfVVVVUx0a5SOa6wabZw+PXXBc6Te+yyp9z7LVV+6mYSrEhky53GGLfoVZWVlJCUlNe1bVFQEUOmcu8E5F/ftnB2lIBFhBw4cYPLkydx4443ccsst3HHHHUybNo3CwkJKS0vB11bT1AGgOwvUFnLrrbc2+8eG5o326tYoEllqk4iw66+/nu3bt7dKHzRoEGVlZZhZpXMuYEeAWOXfaN+yW+PEiRP9N9XofJEupjsJ6Tba263R27Zbt9WIRAsFCYk6/o32HenW2N0d3ricfc/O4/OS7zalqTFfIkVBQqKOf6N9y26NXsN1THdr7Dfm70n+xqPAxXYcNeZLpChISETNnTuXSZMmsXv3btLS0igpKWlqtM/MzIzLbo1XDs8m4ar+zdLUmC+RooZriahVq1YFTC8rKwuYXlRUxNKlSyudc7ldma9oo8Z8iRTdSYh0Y2rMl66mICFhpXmjOibeG/MlchQkRLqBeG/Ml8hRm4RIlDm0/hec/euHNJw+RvXzBZRkPkFhYSH5+fmUlJSQnp7O2rVrgYuN+Tt27IjpxnyJHAWJOKPqnug3ePqPmz1fuPAOQI35EhmqbhIRkaAUJOJARuGrTT8isSoePt+ReI8KEiJhEg9fYtJx0fr5UJAQkagUrV+a8UZBQkREglKQEBGRoBQkREQkKAWJGKd6XRHpDAUJkTBQsJbuSkEihumLSUQ6S0FCwi7egldn32+8lZdEFwUJEREJSkFCRESCUpAQEYmwaK5SVJCQiIjmfwrpvvS5Cj0FiRilfxbpzvT5jR4KEiIi3Ui4A2jUBAkzm2Zmu83sEzMrjHR+okhirJZLJz7sMVsmwbSjrOKuTC5l06ZNANkqk86JiiBhZgnA88D/AEYBc81sVGRzFXkNDQ0A6ahcmqhMWou1MgnFlXJDQwMPPPAAwJ+JgTKJpKgIEsAtwCfOuc+cc+eA1cCMCOcp4srLywHOxnK5XO6Ked2pTEJdLRDseN2pTC4lVGVWXl7OiBEjAM519zKJtGgJEsOAfX7Pq720uFZTUwNwzi8pZsulvV8O0VomLfPf+DwcgSJay6S9uqKsampqGD58uH9SVJZJR5cVDme7RM+wnaltFiDNtdrIbBGwyHt6wsx2e4+vAQ53Ud78hfs8A4HUFq81K5dYKhN7ol3nuWSZQNBy6dIy8ct/uM7TeK4GOl4mjceI6GflEn/7y2JPNJVJInCt30vdqkwu5TLLzP8817a1YSvOuYj/AJOA1/2ePwI8chn7bw1TPsN6ns6Ui8ok8nkN17n0/6My6crzREt1038BmWZ2nZn1AuYA6yOcp2igcmlNZdKayqQ1lUmIREV1k3Ou3sy+B7wOJAC/cc7tjHC2Ik7l0prKpDWVSWsqk9CJiiAB4JzbCGzs4O4rQpmXaDpPJ8pFZdLGMbpYuM7TdC79/7Q+j8okNOcxr75KRESklWhpkxARkSgU1UHiUlN1mNlAM3vZzHaYWbmZZfu9ttfMPjSzCjPbeonz/MbMDppZZZDXzcye8fKxw8zGtTePITxPwPfTjjL6n94+FWZWaWYNZpbUVj6D5O1S57nazP7TzLab2U4zu/dyz9HO8wT9m1/GMVQmrV/vVmXSznO1WS7hKpN2nitqPiuthKP7VQe7bCUAnwLXA72A7cCoFtv8b+Bn3uORQJnfa3uBa9p5rluBcUBlkNdvB17DN55jIvBee/MYivMEez8dOP8/AG910d/iJ8AT3uPBQB3QK5x/c5VJfJRJKMolXGXS3T4rgX6i+U6iPVN1jALKAJxzHwMZZjbkck/knHsH3x8lmBnAi85nCzDAzFLamcdQnCeYy53OZC6wqo3XO3MeB/Q3MwP64Xuf9V1wnkv9zVUmsV8m7T1XW+USrjJp77mi5bPSSjQHifZM1bEduAvAzG7BN5IwzXvNAW+Y2fvmG1XZFXkJ9XQibR0v0Ptp9/nNrA8wDfj3EOer0XNAFvA58CGwxDl3oQvO09bfvL3HwNtfZdJCNymT9p6rrXIJV5m091zR8llpJZqDRHum6igGBppZBfB9YBsXo+/fOefG4ZsZ8wEzu7UL8tKu6URCcB4I/H4u5/z/AGx2zrV1J9ORfDX6GlCBb3qIHOA5M0vsgvO09Tdv7zEaqUxa6w5l0t5ztVUu4SqT9p4rWj4rrUTNOIkAqgH/GbrS8EXZJs65Y8C94Gv0BfZ4PzjnPvd+HzSzl/Hdir0T4rz0ulQeQ3SeYO9n82Wcfw4dv12+5N8C39+h2PkqOz8xsz346jzLQ3metv7ml5HXRiqT1rpDmbTrXJcol3CVSXvPFS2fldY60hATjh98Aewz4DouNsKMbrHNALzGHeB+fPX5AH2B/n6P/wRMu8T5MgjeoHwHzRuUy9ubxxCdJ+D7ae/5gavx1XH27cK/xQvAz73HQ4Aa2tlxIBR/88s5hsqke5dJKMolXGXS3T4rAY/b0Tcejh98vX3+jK/FvshL+w7wHe/xJKAK+Bj4D2Cgl369V0DbgZ2N+7ZxnlVALXAeXzRe2OI8hm9RpE/x1RfmtpXHUJ+nrfdzqTLyns8HVnfx3yIVeMPLdyXwrXD+zS/nGCqT7l8moSiXcJVJd/ustPzRiGsREQkqmhuuRUQkwhQkREQkKAUJEREJSkFCRESCUpAQEZGgFCRERCQoBQkREQlKQUJERIL6/xBZHR9pcC5BAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Histogram of results to get a sense of what's going on\n",
    "\n",
    "fig, ax = plt.subplots(1, 5)\n",
    "ax[0].hist(np.real(pole_dist[0][:, 0]), bins=100, range=[0.95, 1.05])\n",
    "\n",
    "for k in range(1,5):\n",
    "    ax[k].hist(np.real(pole_dist[k][:, 0]), bins=100, range=[0.7, 0.9])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% interval for pole 1 in  sector 0: [1.00000,1.00000]\n",
      "95% interval for pole 1 in  sector 1: [0.71689,0.84102]\n",
      "95% interval for pole 1 in  sector 2: [0.82377,0.87234]\n",
      "95% interval for pole 1 in  sector 3: [0.85138,0.89265]\n",
      "95% interval for pole 1 in  sector 4: [0.81450,0.84550]\n"
     ]
    }
   ],
   "source": [
    "# Bootstrapped confidence intervals.\n",
    "alpha = 0.05\n",
    "cut_off_max = int(NUM_RUNS * (1 - (alpha / 2)))\n",
    "cut_off_min = int(NUM_RUNS * ((alpha / 2)))\n",
    "\n",
    "for k in range(5):\n",
    "\n",
    "    top = np.sort(pole_dist[k][:, 0])[cut_off_max]\n",
    "    bottom = np.sort(pole_dist[k][:, 0])[cut_off_min]\n",
    "\n",
    "    print(\n",
    "        f\"{int(100*(1-alpha))}% interval for pole {i+1} in  sector {k}: [{bottom:0.5f},{top:0.5f}]\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Fidelity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.3825\n",
      "0.8765000000000001\n",
      "0.9675824445967607\n"
     ]
    }
   ],
   "source": [
    "##Compute average fidelity from Majorana fidelities (extracted from fitting)\n",
    "\n",
    "s = (1/4) *(binom(4,0)* 1+ binom(4,1)*0.78 + binom(4,2)*0.85 + binom(4,3)*0.87 + binom(4,4)*0.83)\n",
    "\n",
    "F = (s+ 1)/5\n",
    "print(s)\n",
    "print(F)\n",
    "print(F**(1/4))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}