adaptive_threshold_selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive Approximate Bayesian Computation Tolerance Selection\n",
    "\n",
    "To use ABC-SMC one has to choose the tolerance threshold or acceptance quantiles beforehand to run the method. However, Simola et al. [1] describe an adaptive strategy for selecting the threshold based on density ratios of rejection ABC posteriors. This tutorial teaches you how to use the adaptive threshold selection method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# from scipy.stats import norm\n",
    "import scipy.stats as ss\n",
    "import elfi\n",
    "import logging\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import gaussian_kde\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Set an arbitrary global seed to keep the randomly generated quantities the same\n",
    "seed = 10\n",
    "np.random.seed(seed)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We reproduce the Example 1 from [1] as a test case for `AdaptiveThresholdSMC`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_mixture(theta, batch_size=1, random_state=None):\n",
    "    sigma1 = 1\n",
    "    sigma2 = np.sqrt(0.01)\n",
    "    sigmas = np.array((sigma1, sigma2))\n",
    "    mixture_prob = 0.5\n",
    "    random_state = random_state or np.random\n",
    "    \n",
    "    scale_array = random_state.choice(sigmas,\n",
    "                                      size=batch_size,\n",
    "                                      replace=True,\n",
    "                                      p=np.array((mixture_prob, 1-mixture_prob)))\n",
    "    observation = ss.norm.rvs(loc=theta,\n",
    "                              scale=scale_array,\n",
    "                              size=batch_size,\n",
    "                              random_state=random_state)\n",
    "    \n",
    "    return observation\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "yobs = 0\n",
    "model = elfi.ElfiModel()\n",
    "elfi.Prior('uniform', -10, 20, name='theta', model=model)\n",
    "elfi.Simulator(gaussian_mixture, model['theta'], observed=yobs, name='GM')\n",
    "elfi.Distance('euclidean', model['GM'], name='d');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABC-SMC Round 1 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 2 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 3 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 4 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 5 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 6 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 7 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 8 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 9 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 10 / 10\n",
      "Progress [==================================================] 100.0% Complete\n"
     ]
    }
   ],
   "source": [
    "smc = elfi.SMC(model['d'], batch_size=500, seed=1)\n",
    "thresholds = [1., 0.5013, 0.2519, 0.1272, 0.0648, 0.0337, 0.0181, 0.0102, 0.0064, 0.0025]\n",
    "smc_samples = smc.sample(1000, thresholds=thresholds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adaptive threshold selection ABC (`elfi.AdaptiveThresholdSMC`) can be used in similar fashion as `elfi.SMC`. One does not need to provide a list of thresholds but user can set densityratio-based termination condition (`q_threshold`) and a limit for the number of iterations (`max_iter`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABC-SMC Round 1 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 2 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 3 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 4 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 5 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 6 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 7 / 10\n",
      "Progress [==================================================] 100.0% Complete\n",
      "ABC-SMC Round 8 / 10\n",
      "Progress [==================================================] 100.0% Complete\n"
     ]
    }
   ],
   "source": [
    "adaptive_smc = elfi.AdaptiveThresholdSMC(model['d'], batch_size=500, seed=2, q_threshold=0.995)\n",
    "adaptive_smc_samples = adaptive_smc.sample(1000, max_iter=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare visually the approximated posterior and the true posterior, which in this case is available. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_mixture_density(theta, sigma_1=1, sigma_2=0.1):\n",
    "    y = 0.5 * ss.norm.pdf(theta, loc=0, scale=sigma_1) + 0.5 * ss.norm.pdf(theta, loc=0, scale=sigma_2)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method: SMC\n",
      "Number of samples: 1000\n",
      "Number of simulations: 1352000\n",
      "Threshold: 0.0025\n",
      "Sample means: theta: 0.0181\n",
      "\n",
      "Method: AdaptiveThresholdSMC\n",
      "Number of samples: 1000\n",
      "Number of simulations: 49500\n",
      "Threshold: 0.236\n",
      "Sample means: theta: -0.042\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(smc_samples)\n",
    "print(adaptive_smc_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute Kernel density estimates of the posteriors based on the approximate posterior samples and visualise them in a density plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAI/CAYAAABtd2SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZhdVZ3v/886c52hKlNVJSRkQBIImUMCCSCDCsSBQVDTaA/QDt0K3f70ts/Vn9qgrc/1Ko19afvXNNo26nXoVsFGWltMt8yEGDBMSUiCBDLVqcpUZ6oz1dm/P3btk6rUXHXqjO/X8+SB2sPaK+aP+OH7XWsZy7IEAAAAAEC5uCo9AQAAAABAYyGIAgAAAADKiiAKAAAAACgrgigAAAAAoKwIogAAAACAsiKIAgAAAADKylOpD8+aNctauHBhpT4PAAAAAJhCzz777FHLslqHulexILpw4UJt3769Up8HAAAAAEwhY8zrw92jNRcAAAAAUFYEUQAAAABAWRFEAQAAAABlVbE1ogAAAADqUy6X08GDB5VOpys9FZRBIBDQvHnz5PV6x/wOQRQAAABASR08eFCRSEQLFy6UMabS08EUsixLx44d08GDB7Vo0aIxv0drLgAAAICSSqfTmjlzJiG0ARhjNHPmzHFXvwmiAAAAAEqOENo4JvJnTRAFAAAA0DDC4XClpwARRAEAAAAAZUYQBQAAAFCXrr/+ep1//vlatmyZ7r333uL1T3ziE1q2bJne+ta3qqurS5K0b98+ve1tb9OqVau0du1avfrqq4PG+/GPf6zly5dr1apVuvTSSyVJ9913n66//npdeeWVWrhwob7xjW/orrvu0po1a7RhwwYdP358zOM3EoIoAAAAgLr07W9/W88++6y2b9+uu+++W8eOHVMymdS6dev08ssv67LLLtMXvvAFSdIHPvAB3XrrrXr++ef11FNPac6cOYPG++IXv6hf/epXev755/Xggw8Wr7/00ku6//779dvf/laf/exnFQwG9bvf/U4bN27Ud7/73TGP30g4vgUAAADAlPnCz1/WzsOxko553hnNuv2aZaM+d/fdd+uBBx6QJB04cEB79+6Vy+XS5s2bJUl/+Id/qBtuuEHxeFyHDh3Su9/9bkn2uZhDufjii3XzzTfrfe97n2644Ybi9SuuuEKRSESRSEQtLS265pprJEkrVqzQCy+8MObxGwkVUQAAAAB155FHHtGWLVv09NNP6/nnn9eaNWuGPGJkpB1fP/vZz2r16tVavXq1JOmee+7Rl770JR04cEDnn3++jh07Jkny+/3Fd1wuV/Fnl8ulfD5fyt9W3aAiCgAAAGDKjKVyORW6u7s1ffp0BYNB7d69W1u3bpUkFQoF/eQnP9Ef/MEf6Ac/+IEuueQSRSIRzZs3Tz/72c90/fXXK5PJqLe3V1/+8pf15S9/uTjmq6++qgsvvFAXXnihfvnLX+rAgQNjmstw4weDwSn5vdcCKqIAAAAA6s6mTZuUz+e1dOlSffrTn9aGDRskSaFQSNu2bdPy5cv13//93/rrv/5rSdL3vvc93X333Vq5cqUuuugidXR0DBrzU5/6lFasWKHly5froosu0qpVq8Y8n7GM30iMZVkV+fC6deus7du3V+TbAAAAAKbOrl27tHTp0kpPA2U01J+5MeZZy7LWDfU8FVEAAAAAQFkRRAEAAAAAZUUQBQAAAACUFUEUAAAAAFBWBFEAAAAAQFkRRAEAAAAAZUUQBQBginRnuvWBX3xA1/7s2uKvGx+8US8fe7nSUwMASLrvvvt02223Tejd/fv36wc/+EHx5+3bt+sv//IvSzW1QR555BG9613vmrLxy40gCgDAFNl9fLde6HpBbU1tWjJ9ic6edrb2nNijbUe2VXpqAIBJOj2Irlu3TnfffXcFZ1RbCKIAAEyRzlSnJOlzGz6nOy+7U3ddfpeCnmDxOgBg6lx//fU6//zztWzZMt17773F6//yL/+iJUuW6IILLtCTTz5ZvP7zn/9cF154odasWaO3ve1tikajkqQ77rhDf/RHf6SNGzdq8eLF+uY3vylJ+vSnP63HH39cq1ev1te//vVixbJQKGjhwoU6efJkcezFixcrGo2qq6tLN954o9avX6/169cP+L5j//79evOb36y1a9dq7dq1euqpp4r3YrGY3vnOd+qcc87Rn//5n6tQKEiS/vM//1Nr167VqlWr9Na3vnXI/z0+/elP67zzztPKlSv1V3/1V5Kkm2++WR/96Ee1YcMGnXXWWXrkkUf0p3/6p1q6dKluvvnm4rtjGX+8PCUZBQAADOIEzrZgW/FaW7CNIAoAZfDtb39bM2bMUE9Pj9avX68bb7xR2WxWt99+u5599lm1tLToiiuu0Jo1ayRJl1xyibZu3SpjjL71rW/pq1/9qv72b/9WkvTCCy9o69atSiaTWrNmjd75znfqK1/5iu6880499NBDkuzWWUlyuVy67rrr9MADD+iWW27RM888owULFqi9vV3vf//79YlPfEKXXHKJ3njjDV199dXatWvXgHm3tbXp17/+tQKBgPbu3aubbrpJ27dvlyRt27ZNO3fu1IIFC7Rp0ybdf//9uuyyy/ThD39Yjz32mBYtWqTjx48P+t/i2LFjeuCBB7R7924ZYwaE5BMnTujpp5/Wgw8+qGuvvVZPPvmkvvWtb2n9+vXasWOH5s6dO+r4E0EQBQBginSmOhX2hhX0BovXCKIAGs4vPy11vFjaMWevkN7+lREfufvuu/XAAw9Ikg4cOKC9e/eqo6NDl19+uVpbWyVJmzdv1p49eyRJBw8e1ObNm3XkyBFls1ktWrSoONZ1112npqYmNTU16YorrtC2bds0bdq0Yb+9efNmffGLX9Qtt9yiH/3oR9q8ebMkacuWLdq5c2fxuVgspkQioXA4XLyWy+V02223aceOHXK73cX5SdIFF1ygs846S5J000036YknnpDf79ell15anO+MGTMGzaelpUWBQEAf/OAH9a53vWvAWtNrrrlGxhitWLFC7e3tWrFihSRp2bJl2r9/vw4cODDq+BNBay4AAFOkM9U5oBoqEUQBoBweeeQRbdmyRU8//bSef/55rVmzRul0esR3/uIv/kK33XabXnzxRf3TP/3TgOeNMQOePf3n023cuFH79u1TV1eXfvazn+mGG26QJBUKBW3dulU7duzQjh07dOjQoQEhVJK+/vWvq729Xc8//7y2b9+ubDY7oXlcffXVWr16tT70oQ/J4/Fo27Ztes973qOHHnpImzZtKj7n9/sl2ZVc59+dn/P5/Ii/z8mgIgoAwBQZKoi2BlvV2dMpy7JG/T8yAFAXRqlcToXu7m5Nnz5dwWBQu3fv1tatWyVJF154oT7+8Y/r2LFjam5u1o9//GOtWrWq+M7cuXMlSd/5zncGjPfv//7v+sxnPqNkMqlHHnlEX/nKV3TkyBHF4/Ehv2+M0bvf/W598pOf1NKlSzVz5kxJ0lVXXaW///u/16c+9SlJ0o4dO7R69epBc583b55cLpe+853vqLe3t3hv27Zteu2117RgwQL967/+qz7ykY9ow4YN+tjHPqbXXnut2Do7Y8YM/epXvyq+l0gklEql9I53vEMXX3xxsao6FsONP1kEUQAApkhnT6cuaLlgwLX2YLvyhbxOZE5oRqA07U0AgIE2bdqke+65R0uXLtU555yjDRs2SJLmzJmjO+64Qxs3btS0adMGhMA77rhD733vezV9+nS95S1v0WuvvVa8t3LlSl1xxRU6evSoPv/5z+uMM85Qa2ur3G63Vq1apZtvvrm41tSxefNmrV+/Xvfdd1/x2t13361bb71VK1euVD6f16WXXqp77rlnwHsf+9jHdOONN+q73/2uNm3apFAoVLy3fv163Xbbbdq3b5+uuOIKvfvd75bL5dK9996rG264QYVCobjGtL94PK7rrrtO6XRalmXprrvuGvP/lq2traOOPxHGsqxJDzIR69ats5xFtwAA1JuCVdD53ztfNy+/WR9f+/Hi9V+//mt98pFP6sfX/Fjnzji3gjMEgKmza9cuLV26tNLTKIk77rhD4XC4uNMshjbUn7kx5lnLstYN9TxrRAEAmALH08eVt/JDrhGVxDpRAEBDozUXAIApEE3Z588NCqJNBFEAqCV33HFHpadQl6iIAgAwBbpSXZJOBU/HrOAsGRmCKACgoRFEAQCYAk7QPL0i6nV5NSMwgyAKAGhoBFEAAKZANBWVy7g0s2nmoHttwbZi6y4AAI2IIAoAwBToTHVqVmCWPK7B2zG0BduKrbsAADQigigAAFOgK9Wl1mDrkPfagm205gJAlXj88ce1bNkyrV69Wj09PZWeTsMgiAIAMAWiqeig9aGOtmCbTmROKNubLfOsAKAxWZalQqEw5L3vf//7+sxnPqMdO3aoqalpUmNh7AiiAABMgc5U54hBVJK6emjPBYCpsn//fp1zzjn64z/+Yy1fvlzf+973tHHjRq1du1bvfe97lUgk9K1vfUv/9m//ps9//vP6wAc+IEn62te+pvXr12vlypW6/fbbhxzrwIEDevjhhweNJ0kLFy7U7bffrrVr12rFihXavXu3JCmRSOiWW27RihUrtHLlSv30pz+VpGHHqXcEUQAASiydTyuWjak92D7kfSeI0p4LAFNr7969+tjHPqZHH31U//zP/6wtW7boueee07p163TXXXfpQx/6kK699lp97Wtf0/e//309/PDD2rt3r7Zt26YdO3bo2Wef1WOPPTZgrJdfflmhUEhf+tKXBo3nmDVrlp577jl99KMf1Z133ilJ+pu/+Ru1tLToxRdf1AsvvKC3vOUtOnr06Ijj1LPBOygAAIBJcTYiGmmNqCR2zgXQEP73tv+t3cd3l3TMc2ecq/95wf8c9bkFCxZow4YNeuihh7Rz505dfPHFkqRsNquNGzcOev7hhx/Www8/rDVr1kiyq5h79+7V/Pnzi2NJ0tatW0cc74YbbpAknX/++br//vslSVu2bNGPfvSj4jPTp08f87zqEUEUAIAScwLmcK25TqW0M0lFFACmUigUkmSv67zyyiv1wx/+cMTnLcvSZz7zGf3Zn/3ZgOv79+8vjjWW8fx+vyTJ7XYrn8+P+L2xzKseEUQBACgxp+V2uNbcZl+zfC4fa0QBNISxVC6n2oYNG3Trrbdq3759Ovvss5VMJnXo0CEtWbJkwHNXX311cb1oOBzWoUOH5PV6Jzxef1deeaX+4R/+QX/3d38nSTpx4sSExqkXrBEFAKDEnCA6XEXUGKO2YButuQBQJq2trbrvvvt00003aeXKldq4cWNxE6H+rrrqKr3//e/Xxo0btWLFCr3nPe9RPB6f8Hj9fe5zn9OJEye0fPlyrVq1Sr/5zW8mNE69MJZlVeTD69ats7Zv316RbwMAMJW++tuv6id7fqJn3v+MjDFDPvMnv/wTGWN036b7yjs5ACiDXbt2aenSpZWeBspoqD9zY8yzlmWtG+p5KqIAAJSYc3TLcCFUstt22TUXANCoCKIAAJTYSGeIOlqDrepKdalSnUkAAFQSQRQAgBIbSxBtC7Yp3WufNwoAQKMhiAIAUEKWZdlBtGnkIFo8woX2XAB1io6PxjGRP2uCKAAAJXQyc1K5Qm5MrbmS1JXiCBcA9ScQCOjYsWOE0QZgWZaOHTumQCAwrvc4RxQAgBIa7egWh3OfI1wA1KN58+bp4MGD6uriP7Y1gkAgoHnz5o3rHYIoAAAl5ATLsQZRWnMB1COv16tFixZVehqoYrTmAgBQQk6r7WhB1O/2a5p/GkEUANCQCKIAAJSQEyxbm1pHfbY12KrOHoIoAKDxEEQBACihaCqqGYEZ8rq9oz7bFmyjIgoAaEgEUQAASqgz1Vk8mmU07cF2gigAoCERRAEAKKHOVGfxaJbRtAXbdKznmHKF3BTPCgCA6jJqEDXGnGmM+Y0xZqcx5mVjzMeHeMYYY+42xuwzxrxgjFk7NdMFAKC6dfV0jbpRkaO1qVWWLB3rOTbFswIAoLqMpSKal/Q/LMs6T9IGSbcaY8477Zm3S1rc9+sjkv6xpLMEAKAGZHuzOp4+PuYg6rTw0p4LAGg0owZRy7KOWJb1XN+/xyXtkjT3tMeuk/Rdy7ZV0jRjzJySzxYAgCrW1WMf3TLWNaKcJQoAaFTjWiNqjFkoaY2kZ067NVfSgX4/H9TgsAoAQF1zzhAdy9EtkoprSaOp6JTNCQCAajTmIGqMCUv6qaT/x7Ks2EQ+Zoz5iDFmuzFme1dX10SGAACgap1In5AkzQjMGNPz0/3TJUndme4pmxMAANVoTEHUGOOVHUK/b1nW/UM8ckjSmf1+ntd3bQDLsu61LGudZVnrWlvH9l+LAQCoFfFcXJIU8UXG9Lzb5VbIG1I8G5/KaQEAUHXGsmuukfTPknZZlnXXMI89KOmP+3bP3SCp27KsIyWcJwAAVc8JlGFfeMzvhL1hgigAoOF4xvDMxZL+SNKLxpgdfdf+X0nzJcmyrHsk/ULSOyTtk5SSdEvppwoAQHVLZBOSpIh3bBVRya6eJnKJqZoSAABVadQgalnWE5LMKM9Ykm4t1aQAAKhFiVxCAXdAXrd3zO9EfJFigAUAoFGMa9dcAAAwvHg2Pq62XMluzY1lJ7QHIAAANYsgCgBAicSzcYW94wyivjCtuQCAhkMQBQCgRBK5hJp9zeN6p9nXTGsuAKDhEEQBACiRibbmxrNx2dstAADQGAiiAACUyERbc/NWXune9BTNCgCA6kMQBQCgRBK5hCK+sR/dIqnYykt7LgCgkRBEAQAokXg2Pu4g6lRQ49n4VEwJAICqRBAFAKAEcr05ZXozE2rNlaR4jiAKAGgcBFEAAErACZK05gIAMDqCKAAAJeAESVpzAQAYHUEUAIAScIIkrbkAAIyOIAoAQAnQmgsAwNgRRAEAKIGJtuY2eZrkNm5acwEADYUgCgBACRRbc33ja801xijkDRFEAQANhSAKAEAJOEFyvBVR551EjtZcAEDjIIgCAFACTpAMeULjfjfii1ARBQA0FIIoAAAlEM/GFfKG5Ha5x/1u2BsmiAIAGgpBFACAEohn4+M+usUR9oVpzQUANBSCKAAAJZDIJSa0PlSyj3Dh+BYAQCMhiAIAUAKJ7MSDKK25AIBGQxAFAKAEYtnYpFtzC1ahxLMCAKA6EUQBACiBybbmWrKUyqVKPCsAAKoTQRQAgBKYbGuuJNpzAQANgyAKAMAkWZY16V1zJSmeI4gCABoDQRQAgElK96aVt/ITrog677FzLgCgURBEAQCYJCdATjiIeu33aM0FADQKgigAAJPkBEhacwEAGBuCKAAAk+QESFpzAQAYG4IoAACTNOnWXB+tuQCAxkIQBQBgkpyK6ERbc/1uv7wuL625AICGQRAFAGCSnErmRCuizru05gIAGgVBFACASZpsa67zLq25AIBGQRAFAGCS4tm43MatJk/ThMcIe8O05gIAGgZBFACASYpn4wp5QzLGTHiMsC9May4AoGEQRAEAmKRELjGptlxJavY105oLAGgYBFEAACYpkZ18EA17qYgCABoHQRQAgEmKZWMTPrrFEfaxRhQA0DgIogAATFIpWnMjvoh68j3KF/IlmhUAANWLIAoAwCSVojU34o0UxwIAoN4RRAEAmKR4Ll6S1lxnLAAA6h1BFACASShYhdJURH1URAEAjYMgCgDAJKRyKVmyStaayxEuAIBGQBAFAGASEjm7gklrLgAAY0cQBQBgEmLZmCTRmgsAwDgQRAEAmAQnODoVzYmiNRcA0EgIogAATILTmusEyYkK+UKSaM0FADQGgigAAJPgtOZOtiLqdXnV5GmiNRcA0BAIogAATIITHCe7RlSyq6q05gIAGgFBFACASSi25pYgiIZ94eJ4AADUM4IoAACTEM/G5XV55Xf7Jz1W2BemIgoAaAgEUQAAJiGejZekGirZVVWCKACgERBEAQCYhEQ2Ubog6o3QmgsAaAgEUQAAJiGeiyvsndyOuQ5acwEAjYIgCgDAJNCaCwDA+BFEAQCYhFK35uYKOWV6MyUZDwCAakUQBQBgEkrdmiuJqigAoO4RRAEAmIRSt+ZKdpUVAIB6RhAFAGCC8oW8evI9xUrmZEW8dhClIgoAqHcEUQAAJiiZS0o6FSAnq9iamyOIAgDqm6fSEwAAoFbFsjFJGrE199DJHu08HNO+zoT2dSZ06GRKn3vneVo+t2XQs7TmAgAaBUEUAIAJcgLjcK25uztievv/eVyWZf/cFvHrZE9O3316v776nlWDnqc1FwDQKAiiAABMUCJnB9HhWnN/8WKHjKQffHiDzjujWS1NXn3y33boP1/q0JeuXyGfZ+AKGSfQOuMCAFCvWCMKAMAEOZXL4Sqiv94Z1fkLpmvjm2aqpckrSbpm5RmKpfN6fG/XoOdD3pCMDBVRAEDdI4gCADBBI1VEDxxPadeRmK48r33A9YvPnqWWJq8eeuHIoHdcxqWQN0RFFABQ9wiiAABMkLNrbtAbHHRvy66oJOnK82YPuO7zuLRp2Ww9/HKH0rneQe8FvcHiuAAA1CuCKAAAE5TKpSTZLbWne/jlqBa3hbVo1uB716w6Q8lsrx55pXPQvZA3RBAFANQ9gigAABOUzCXlMi753f4B10+mstq2//igtlzHhrNmaGbIp58P0Z4b8oSKARcAgHpFEAUAYIJS+ZRCnpCMMQOu/+aVTvUWrGGDqMft0ttXzNZ/7YoqmckPuEdFFADQCAiiAABMUDKXHHJ96K93RtUW8WvVvGnDvvuulWconSvov3YPbM8NeoNK5gmiAID6RhAFAGCCkrnkoPWh6VyvHn2lS29d2i6XywzzprR+4Qy1N/v10POHB1wPeWnNBQDUP4IoAAATlMqnFPQMrIg+/eoxJbO9umrZ0G25DrfL6B0r5uiRV7oUS+eK12nNBQA0AoIoAAATlMqlBlVEH94ZVcjn1kVvmjnq++9aeYayvQVt2RktXuP4FgBAIyCIAgAwQaevES0ULG3ZFdVl57TK73GP+v7a+dMU8Xv0uzdOFq+FPCHlCjnlenMjvAkAQG0jiAIAMEGnrxF96XC3uuKZYXfLPZ0xRktmR/RKNF685gRbqqIAgHpGEAUAYIJOb8198VC3JGndghljHmNJe0SvdMRlWZYkFcdL5dmwCABQvwiiAABM0OmbFe2NJhT0uTV3WtOYxzh3dkTdPTlFYxlJVEQBAI2BIAoAwATkC3llejMD1ojuica1uC084rEtpztndkSStLsjJsleIyoRRAEA9Y0gCgDABDhBsX9r7p5oQovbI+Ma55y+51/piA8Yj7NEAQD1jCAKAMAEOEHRCY4nklkdTWS0pD08rnGmh3xqi/iLGxY54yXzVEQBAPWLIAoAwAQ4mwk5a0T39AXJ8VZEJbs916mIskYUANAICKIAAEyAExSd4LinMyHpVKvteJw7O6K9nQnlewunKqIEUQBAHSOIAgAwAaevEd0bjSvi92hOS2DcY50zu1nZfEH7j6VYIwoAaAgEUQAAJuD0NaJ7onGd3R6WMWPfMdfhVFH3ROPyuXzyGA8VUQBAXSOIAgAwAc5mQs5xK3ujCS1pG39briQtbg/LZaTdHXEZYxT0BgmiAIC6RhAFAGACnIpok7dJxxIZHUtmtXicO+Y6Al63Fs4M6RXnLFFvqLgZEgAA9YggCgDABPRfI7onam9UtGQCGxU5+u+cG/KGqIgCAOoaQRQAgAlI5pJyGZcC7oD2dtoBcrJB9PXjKfVke2nNBQDUPYIoAAATkMqnFPKEZIzRKx1xRQIetTf7JzzeOe0RWZa0tzOukCfErrkAgLpGEAUAYAKSuWTxDNG90YSWtEcmtGOu45zZdjV1d0ec1lwAQN0jiAIAMAGpXEpBb1CWZWlPZ1xLJrhRkWPBzJACXpde6Yjbrbl5gigAoH4RRAEAmIBkPqmQJ6SuREYnUzktnuDRLQ63y2hxW0R7olREAQD1jyAKAMAEpHIphbwh7S3BjrmOJe2RYmtuKpeSZVmTHhMAgGpEEAUAYAKcNaJ7os6OuZNrzZWkc2dH1BXPyFh+9Vq9yvRmJj0mAADViCAKAMAEnAqiCbU0edUamfiOuQ5nw6J4yl38BgAA9YggCgDABPTkexTyhLQ3Gtc5k9wx13FuXxA9kbDH4ggXAEC9IogCADAB/VtzF5egLVeSWiN+NQc8Ot4XRNk5FwBQrwiiAACMU76Qt9dvFvyKpfMl2ahIkowxOnNGUCcT9l/PtOYCAOoVQRQAgHFyAmKix17LWaqKqCSdOT2oo/GB3wEAoN4QRAEAGCdn7Wa8L4ieNat0QXTe9CZ1dvd9J88aUQBAfSKIAgAwTk5ATPS45XO71FaCHXMd86Y3KZPx2t9hsyIAQJ0iiAIAME5Oy2x30qW505vkck1+x1zHvOlBWQX/gO8AAFBvCKIAAIyTExCPxY3mTW8q6dhnzghKBFEAQJ0jiAIAME5Oy+zRmFXyIDp3epMkt9zy0poLAKhbnkpPAACAWuOc73ki4dK86cGSjh32ezQ96JXLBKiIAgDqFhVRAADGqVipLPg1d1ppK6LSqfZcJ/ACAFBvRg2ixphvG2M6jTEvDXP/cmNMtzFmR9+vvy79NAEAqB5OpdIq+EvemivZO+f29vqoiAIA6tZYKqL3Sdo0yjOPW5a1uu/XFyc/LQAAqlcyl5SRS7K8JW/Nleydc3M5n5JZgigAoD6NGkQty3pM0vEyzAUAgJqQyqfkMX55S3yGqOPM6U0q9PoUyyRKPjYAANWgVGtENxpjnjfG/NIYs6xEYwIAUJWSuaSMFdDcaaU9Q9ThnCXanaEiCgCoT6XYNfc5SQssy0oYY94h6WeSFg/1oDHmI5I+Iknz588vwacBACi/VC7Vtz609G25kr1G1Cr4WSMKAKhbk66IWpYVsywr0ffvv5DkNcbMGubZey3LWmdZ1rrW1tbJfhoAgIpI5pPK53xTsmOuZFdEVfAr3cs5ogCA+jTpIGqMmW2MMX3/fkHfmMcmOy4AANUqkU0ql/dOyY65ktTkc6vJHVTOSsuyrCn5BgAAlTRqa64x5oeSLpc0yxhzUK761/QAACAASURBVNLtkrySZFnWPZLeI+mjxpi8pB5Jf2DxtyYAoI51pxOyCj7NmzE1QVSSWgJhHZelnnyPgt6paQEGAKBSRg2ilmXdNMr9b0j6RslmBABAlUtkk1KhfcrWiErSzKaIjhfsjZEIogCAelOqXXMBAGgYqbyzWdHUVURbwy2SpHiWI1wAAPWHIAoAwDhlelNyWQG1RQJT9o3ZfUH0QPeJKfsGAACVQhAFAGAc8oW8epVT2BeUewrOEHXMnTZNkvTGCYIoAKD+EEQBABgH52zP6YHIlH5n/rTpkqSDJwmiAID6QxAFAGAcevI9kqSZwakNogunz5AkHYl3T+l3AACoBIIoAADjcLwnJklq61vDOVWmN9lBtzNBEAUA1B+CKAAA4/B635rNOZGpDaIhb0iSdCwVn9LvAABQCQRRAADG4UBfEJ3bMn1Kv9PksY+GOZHm+BYAQP0hiAIAMA6HYiclSfOnT20Qdbvc8hi/Etmk8r2FKf0WAADlRhAFAGAcon2bB53Zd7zKVPK7mmSZjDpi6Sn/FgAA5UQQBQBgHLqS9mZFzf7wlH8r6A3JuDI6eKJnyr8FAEA5EUQBABiH4z325kFBT3DKvxXxEUQBAPWJIAoAwDicTCckmeJmQlNpWiAi487owPHUlH8LAIByIogCADBGmXyvUvmkvKZJxpgp/17YF5LXk6MiCgCoOwRRAADG6PDJtOTKKuCe+mqoJIU8Ibk9GUXZrAgAUGcIogAAjNHBEykZV0Yhb6gs3wt6g5IroyPdVEQBAPWFIAoAwBgdPNEj48oo4itPEA15Qyooo2gsU5bvAQBQLgRRAADG6Eh3WsaVUUtg6o9ukewg2quMEpms4ulcWb4JAEA5EEQBABijju4eeTw5hctYEZUkubKsEwUA1BWCKAAAY9QRy8jtyZZ3jagk48qoo5v2XABA/SCIAgAwRh3dPZIro6AnWJbvhTx24DVsWAQAqDMEUQAAxqijO62C0mWriJ5qzeUIFwBAfSGIAgAwBqlsXrF0RgXlii2zU835TiSYVwdBFABQRwiiAACMQUd3WnLZ6zSdltmp5lREp4Us+/sAANQJgigAAGPQEUvLuLKSVPbW3OZgLxVRAEBdIYgCADAGHX1niEoqW2uuE0RDTb3smgsAqCsEUQAAxqAj1q81t1zHt/Ttztvky+toIqNsvlCW7wIAMNUIogAAjEG0O61QIC9JZTu+pcnTJJdxye/LSZI647TnAgDqA0EUAIAxONKd1vSwJal8FVFjjIKeoDxee20qR7gAAOoFQRQAgDGIxtKKBO3W2HIFUclej+py20H0CDvnAgDqBEEUAIAxONKdVthpzS3TZkWSHXotYwdQjnABANQLgigAAKPI9xZ0NJFRU5nXiDrfyhXS8ntctOYCAOoGQRQAgFF0JTIqWJLfl5fLuNTkaSrbt0PekFL5lOa0BGjNBQDUDYIoAACjcAKgx5NVyBOSMaZs3w56g0rmkmpvDlARBQDUDYIoAACjiPYFUePOlHV9qGRXRJO5pGZTEQUA1BGCKAAAoygGQJMp6465khTyhJTKpTS7JaDOWEaWZZX1+wAATAWCKAAAo4jG0vJ5XMoVesofRJ2KaHNA2d6CjiezZf0+AABTgSAKAMAoOmJptTf7lcwny96aG/QGlS1k1dbskcRZogCA+kAQBQBgFEe605rT3KRkLqmQp/wVUUlqDhYkiQ2LAAB1gSAKAMAoorG02lsCSuVSFWnNlaRIsFeSXZ0FAKDWEUQBABiBZVl2RbQloGQ+WfYg6rQCB3w5uYzUQWsuAKAOEEQBABjByVRO2XxBbRG/ktnyB1GnFTjd26PWiJ8gCgCoCwRRAABG4LTCzoq4lbfyFWvNTeVSmt0coDUXAFAXCKIAAIzAqUC2hOzNgsq9a64TRJP5pGa3BKiIAgDqAkEUAIAROBVIZ7OgSq0Rdc4SpSIKAKgHBFEAAEZwpDstY6Qmf05S+YNosSKaS2p2S5Pi6bySmXxZ5wAAQKkRRAEAGEG0O61ZYb8yvT2SVLFzRFO5lGa3+CVxhAsAoPYRRAEAGEFHLK3ZzQElc0lJ5V8j6nP55DEeJXNJtTcHJNnhGACAWkYQBQBgBB3dac1uORVEw95wWb9vjFHQG1Qyl9ScliZJdrswAAC1jCAKAMAITq+IlnuNqPPNVN4+vsWZEwAAtYwgCgDAMHqyveruyQ2oiJa7NVeyg2gyl1STz63mgEdRgigAoMYRRAEAGIZTeZzdHFAyX7mKqNOaK0lzWppozQUA1DyCKAAAw+joC3yzWwJK5VLyu/3yuDxln0fIE1Iql5IktTX71UlFFABQ4wiiAAAMoyNmH9nitOZWohoqnWrNlaT25oA645mKzAMAgFIhiAIAMIyObjvwOZsVBT3lXx8q9bXm5p0g6ldnPKNCwarIXAAAKAWCKAAAw4jG0gr7PQr57XM8w77yHt3iCHqCxdbc9uaAeguWjiWzFZkLAAClQBAFAGAY0Vha7c1+SapoRTTktdeIWpaltkigODcAAGoVQRQAgGHYQdQOfpVeI5q38soWssVg3BkniAIAahdBFACAYURjGc3uC6KpfKpiQdQ5uzSZSxaDcTTGhkUAgNpFEAUAYAiWZakznlZblVREnTm0RuyKKK25AIBaRhAFAGAIJ1I55XqtgWtEvZVbIypJqVxKXrdLs8I+KqIAgJpGEAUAYAhOxbG9OaCCVVBPvqdyFVHPqYqoJLVGAuqkIgoAqGEEUQAAhnAqiPqLR6eEvRU6vqXfGlFnTp1xKqIAgNpFEAUAYAidfa2vbZGAErmEJFW8NTeZ7wuikQBrRAEANY0gCgDAEJyg19avIuq0yJZb/zWikl0RPZrIKN9bqMh8AACYLIIoAABDiMbTmh70yu9xF1tiq2HXXElqaw6oYEnHktmKzAcAgMkiiAIAMIRoLFM8s9Npia1Ua+7gNaLOWaK05wIAahNBFACAIXTGBp4hKlWuIup1eeVz+Qa05kriCBcAQM0iiAIAMIRoLKP2iB34imtEKxREnW9TEQUA1AuCKAAAp+ktWOpK9GvNrXBFVLLbc50W4Zkhn1xGnCUKAKhZBFEAAE5zLJFRb8EqtsA6x7dUS0XU43ZpZthPay4AoGYRRAEAOI0T8Jw1oqlcSi7jUsAdqNicQt5QsUVYsteJdsapiAIAahNBFACA0zhrL2f3a80NeUIyxlRsTkFvsFgRlaT2SICKKACgZhFEAQA4TbSv0th/jWiljm5xhDyhAUG0rTlARRQAULMIogAAnCYay8gYaVbYJ0lK5VMVXR8qDd2aezSRVa63UMFZAQAwMQRRAABO0xlLa1bYL4/b/msymUtWRRB1ds2VTlVru+K05wIAag9BFACA00Rj6eKOuVJ1BNGgN6hULiXLsiSpOD/OEgUA1CKCKAAAp4nGMmqPnNohtxqCaMgbkiVLPfkeSVJb3/zYsAgAUIsIogAAnKYzni4e3SJVRxANeoLFuUinWnPZsAgAUIsIogAA9JPrLehoIjuoNdcJgpXiBOFU3t6waGbIJ7fL0JoLAKhJnkpPAACAauJs/uNUHC3LUio3wV1zO16SXvqJ1DxXmrtWal8uefyjvzcE5/gYpyLqchm1hv205gIAahJBFACAfpwKo1MRzRayylv5sQdRy5L2PyE9+XfSvi2SjCR7gyG5vNKcldJbb5fOumxc83K+3/8s0fZmvzrZNRcAUIMIogAA9ONUGJ3NgJzg51QkR3Riv/STD0qHtkuhVuktn5fWf1DKxKXDv5MOPSft+rn0f2+Urv//pJXvG/O8Qp6+1tx+Z4m2NQd04HhquFcAAKhaBFEAAPpxNv9xWnOdIDpqRTR5zA6YyaPSO++SVr9f8jbZ95qmS9PmS+ddJ735k9KPPiDd/2Gp+6B0ySckY0ad13AV0e37j4/3twgAQMWxWREAAP10dKflcRnNDPkknQp+YW94+JeyKemHm+1g+f5/taugTgg9XaBF+sOfSsvfI/3XF6T/+B9SoXfUeRXXiOb7BdFIQCdSOWXyo78PAEA1oSIKAEA/0VhGbRG/XC67Sjlqa25vXvrpB6WD26XN35Pmbxj9Ix6/dMM3pZa50pP/R/JHpCu/MOIrxV1z+7XmFo9wiWV05ozK7uoLAMB4UBEFAKCfoc4QlYZpzbUs6Rd/Jb3yC+kdX5OWXjP2D7lc0pVflNb8kfTU30uHd4z4+OnniEpSW9+GSpwlCgCoNQRRAAD6icbSA84QdSqQzmZBA+z4vvTsv0iXfFK64MMT++BVfyOFZkkP/oVdXR2G2+VWk6fptDWigb45s3MuAKC2EEQBAOgnGssUA540QkU0E5e2fEGad4H01r+e+AebpkvvuFPqeEF6+hsjPhr0BIcJolREAQC1hSAKAECfdK5X3T25IYPooDWiT3xdSnZKm74ypl1vR3TetdK575Ie+V/SsVeHfSzkDQ1YIzo96JXXbaiIAgBqDkEUAIA+ncUzRE+15jq71A6oiJ58Q3rqG9KK90nzzi/Nx99xp+T2ST//uL32dAghb2jArrnGGLVFAqwRBQDUHIIoAAB9oqedISpJyWxSAXdAHle/jea33CEZl/S220v38eY59uZF+x+Xnv/hkI8EvQNbcyV7w6JOKqIAgBpDEAUAoI+z1nJAEM0nB7blvvGM9NJPpYv+QmqZV9oJrP0T6Yw10mN3Dnm26OmtuZJ9lihrRAEAtYYgCgBAH2etZf9dc5O55Km23EJB+tVnpPBs6eKPl34CLpd00V9Kx1+VXvnloNshT2hQRbS92U8QBQDUHIIoAAB9OmNp+TwutTR5i9dSudSpILrzZ9KhZ+1dcv3hqZnE0mulafPts0VPM3RrbkCxdF492cEVVAAAqhVBFACAPh19Z4iafrvgJnNJBT19rblb/1GacZa06qapm4TbI224VTqwVTrw2wG3Qt6QUvnTWnP72ojZsAgAUEsIogAA9OnoTmt2v/WhUr/W3MM7pIPbpPUftltop9KaP5QCLdLTA6uiIW9IPfke9fZbP+q0EXOECwCglhBEAQDoE42lB2xUJEmpfF9r7m+/KXmD0ur3T/1E/GFp3QelXT+Xjv++eNlpEe5fFXXmyzpRAEAtIYgCACDJsixFY5mhK6LGI734E2nl+6SmaeWZ0AUfkYzbbgfu4+ze23+dqHPmKUEUAFBLCKIAAEj2hj+53kEV0WQuqdCJ16V82m7LLZfmOXbw/d3/lVLHJdm75koacIRLS5NXPo9LnXFacwEAtYMgCgCA+p0h2nIqiPYWetWT71HoyIvSgoul2cvLO6mNt0m5lLT925JOteb2r4gaY9Te7FcnFVEAQA0hiAIAIHujIkkDWnOdtZih1Alp/YfKP6n286SFb5Z2fF+yrFOtufnTzhKNBNisCABQUwiiAADIPrpFGhhEncpj0BeRll5TkXlp5WZ7w6JDzw25RlSyNyyKcnwLAKCGEEQBAJCKra1tfcehSFKqa7ckKbToMsntrci8dN61ktsvvfhvQ64Rlew5d1IRBQDUEIIoAACyK6LTgl4FvO7iteTLP5EkhRZvqtS07PNEl1wtvfRThVx2SD49iLY3B5TI5JXI5CsxQwAAxm3UIGqM+bYxptMY89Iw940x5m5jzD5jzAvGmLWlnyYAAFOro/u0o1ssS8nfPypJCjXPq9Cs+qx8n5TsUujQc5KGWCPaV8VlwyIAQK0YS0X0Pkkj/afgt0ta3PfrI5L+cYRnAQCoStFYeuDRLYeeVbKnS9Kp3WorZvFVUqBFTS//TEZm8BrRiD1vNiwCANSKUYOoZVmPSTo+wiPXSfquZdsqaZoxZk6pJggAQDl0xNIDK6Iv3a+kx14XWvEg6vFL510ns/s/FPQ0DbFG1J53JxsWAQBqRCnWiM6VdKDfzwf7rgEAUBPyvQUdTWSKLa4qFKSXH1Cy9VxJKu5WW1ErN0u5pELGPagi6mywFKU1FwBQI8q6WZEx5iPGmO3GmO1dXV3l/DQAAMPqSmRkWVJ7S19F9MAzUvywku3nSaqCiqgkzb9Iap6nYC4zKIhG/B41ed3snAsAqBmlCKKHJJ3Z7+d5fdcGsSzrXsuy1lmWta61tbUEnwYAYPI6uk87Q/Tl+yVPQMlp8+Q2bgXcgRHeLhOXS1pxo8KZhJLpEwNuGWPU3uxXNE4QBQDUhlIE0Qcl/XHf7rkbJHVblnWkBOMCAFAWTktre3NAKvRKO/9dWnyl4oWMIr6IjDEVnmGfFe9TpNCreOzAoFttzQFacwEANcMz2gPGmB9KulzSLGPMQUm3S/JKkmVZ90j6haR3SNonKSXplqmaLAAAU6FYEW0JSK8/JSWi0rIbFD/+jCK+SIVn18/s5Yp4wupIHR10q705oBcPnqzApAAAGL9Rg6hlWTeNct+SdGvJZgQAQJlF4xl53UYzgj67LdcblJZcrfhjWxT2his9vQEi0+Yr3r1PSh2XgjOK19sjfm2JZWRZVvVUcAEAGEZZNysCAKAaRbvTaosE5LJ6pZ0PSks2Sb6Q4tm4mn3NlZ7eAJGZ5yjuMtK+LQOutzcH1JPrVTyTr9DMAAAYO4IoAKDhdcTS9tEt+x+TUkel5TdIkhK5RHW15kqKzDhbGZdL2d3/MeC6c4RLJ+tEAQA1gCAKAGh4HbG0vT50579LvrB09pWSpFg2Vn1B1G9XaOOvPSL15orX2/t2/I1yhAsAoAYQRAEADS/anVZ7xC/t/bX0prdIXjvUxbNxhX3VtUbUWbMazyXsjZX6tEXsiig75wIAagFBFADQ0BKZvJLZXi11H5Bih6TFV0mScoWcevI9VVcRddasxr1+ac9/Fq+3UREFANQQgigAoKE5R7csSzxjX1hst+Ums0lJqr7NivqCcfyMVdIrv5QsS5IU9nsU9nvUGaciCgCofgRRAEBDc1pZzzz2uDRnlRSZLcluy5VUdRXRU0F0pXTiNenonuK9tma/OqmIAgBqAEEUANDQOrrTalZCka7nim25khTLxSSp+s4RdYLorLPtC6/8snivPRJgjSgAoCYQRAEADa0jltalrhdlrMKAIFr1FVG3R5q9YsA60fZmv6K05gIAagBBFADQ0DpjaV3le15qmiHNPb94PZFNSKq+NaJBT1Au47KD8pK3SweekVLHJdlHuERjGVl960YBAKhWBFEAQEOLdqd0qdkhnf02yeUuXq/WiqgxRhFfxJ7fOZskq2AfOyN759xsvqDuntwoowAAUFkEUQBAQwsff0nTrNiAtlxJimXtNaLVFkQlKeKNKJ6LS3PWSOF2aY+9TrS92TlLlA2LAADVjSAKAGho58aeVkEu6ey3Drgez8ZlZBTyhio0s+EVK6Iul3T2ldKrv5EKvWovniXKOlEAQHUjiAIAGlZvwdIF+e06ElkuBWcMuJfIJRT2huUy1fdXZcQXKa5h1ZuukNInpSM71BZxKqIEUQBAdau+v10BACiT49EDWun6vTrbLx10L56NV2VbrmQHUad1WIsus//5+0fUFrErop1xWnMBANWNIAoAaFjpXb+y/7nobYPuxbKxqg6izmZKCrdK7SukV3+jJp9bzQGPOqmIAgCqHEEUANCwvL//L0WtaQrNXz3oXjwbV9gXrsCsRhf2hk8FUUl60+X2MS7ZVPEIFwAAqhlBFADQmAoFTYtu1ROFFZrd0jTodjW35jb7mpXKp5Qv5O0LZ10u9WalN56yg2iciigAoLoRRAEAjanzZQVyJ7TVWqaZYf+g24lsQs2+5gpMbHROQE7mkvaF+RdJbp/06m/U1uxXtJsgCgCobgRRAEBj+v2jkqTdTWvldplBt6u5IurMq7hhkS8onXmh9PtHNacloGg8o96CVcEZAgAwMoIoAKAxvfaYDnvmyT1t3qBbBatQPL6lGjlrVweuE71Cir6oRYGUeguWjiZYJwoAqF4EUQBA4+nNSa8/qe1arrnTBq8PTeQSsmRVbUXUaRkeEETPulySdG7P7yRJh0/2lHlWAACMHUEUANB4Dj0nZRP6r8y5mtMSGHQ7kU1IUtWvEXXmKUmas1oKTNPc41slSR2sEwUAVDGCKACg8bz2mCTp0dxSzRmiIupUGqu1IjpojagkudzSokvVcuQJSZaOEEQBAFWMIAoAaDyvPaqemct0UhHNnTa4IuoEvGoPogNacyXpTVfIFT+scz1RHemmNRcAUL0IogCAxpJNSQeeUXTmhZKkOcOcISqd2hSo2oQ8IUlSPHdaED3rCknS24O7qIgCAKoaQRQA0FgObJV6s9oXOl+SNGeIimgi17dG1Fuda0TdLrfC3vDANaKSNGORNG2BLnK9xBpRAEBVI4gCABrLa49JLo92uJfK53ZpVsg/6JFqXyMq2XMbsEbUcdblWp59UdGTybLPCQCAsSKIAgAay+8fleau0xtxl2a3BORymUGPOAGvWltzJTuIDlojKkkLL1FTIaEZiT3qLVjlnxgAAGNAEAUANI6ek9KRHdJZl+nwyZ4hj26R7Ipok6dJHpenzBMcu7A3PHQQXXCxJGmddupYIlPmWQEAMDYEUQBA49j/hGQVpEWX6kh3WnOHOLpFss/nrOa2XMk+49RZyzpAy1ylwvN1oWuXDrNOFABQpQiiAIDG8dpjkqdJvWesU0csPeRGRZJdEW32VedGRY5hW3MlZeZu1AWu3epgnSgAoEoRRAEAjWP/49L8DerssdRbsIY8ukWyg2i1V0SH3axIkvesN2uaSSp96MUyzwoAgLEhiAIAGkPquNS5U1pwsQ6ftFtWh2vNjWVjCnurd6Miyd5IKZFNqGAVBt0LLblUkhQ8vLXc0wIAYEwIogCAxvBGXyhbsFFHunskDX2GqGSfI1rtFdFmX7MsWUrlUoPumekLdMS0qe34sxWYGQAAoyOIAgAawxtPSW6fNPd8HT7ZF0RrvDVX0rDrRPcEVmpRcodUGFwxBQCg0giiAIDG8PrT0hlrJW+TDp9MK+z3qDkw+HgWy7JqZrMiScOuEz087Xy1WDGpa3c5pwUAwJgQRAEA9S+btM8PXXCRJOlIt32GqDFm0KM9+R71Wr01XxGNtV8oSSrsf6JscwIAYKwIogCA+nfwt1IhXwyih0+mdcYIGxVJ9mZA1SziHTmIBtvO0mFrhrKvPl7OaQEAMCYEUQBA/Xv9Kcm4pDMvkGRXRM8YbqOibEKSaqYimsglhrw/e1pQWwvnyf3GU5JllXNqAACMiiAKAKh/rz8ltS+XAi1K53p1NJEdfqOinF1hbPbW9hrROS0BPVNYKm/6qHR0bzmnBgDAqAiiAID6ls/arbkLLpYkdXTbZ4gO15rrtLpWe0XUaR0erjXXDqLn2j/spz0XAFBdCKIAgPp2ZIeUT0sLNkqSDvedIXpGy9CtubWyRtTr8qrJ0zRsEJ0R8umw+wzFvbOk158s8+wAABgZQRQAUN9ef8r+53w7iB45aVdE5wxTEa2VNaKSvWHRcGtEjTGa09KkVwKrpP1Psk4UAFBVCKIAgPr2+lPSzMVSuE2SdPikXRGdM0xFtFZacyV7jsNVRCVpdnNAz5rzpESHdPz3ZZwZAAAjI4gCAOpXoVd6Y2uxLVeSDnenNTPkU8DrHvKVeDYuv9svv9tfrllOWMQXGXazIskO249mltg/sE4UAFBFCKIAgPrVuVPKdBc3KpLso1vmDHN0i2SvEa2Faqg0ekV0zrQm/TYxU1aozW7PBQCgShBEAQD16/Wn7X/O71cRPdmjM4Y5ukWyz+UMe6t7oyJH2BcurmkdypyWgHK9UmbuBnvDItaJAgCqBEEUAFC/Xn9Sap4nTZtfvHTkZHrYo1skuzW32VfdZ4g6mn3NI1dE+wL3sVnrpdgh6cT+Ms0MAICREUQBAPXJsqQ3nrbXhxojSYqlc4pn8sNuVCTZQbTWWnOtYSqdzu/ztfAa+wLHuAAAqgRBFABQn47/XkpEB7TlOke3jFYRraUgmrfy6sn3DHl/dl8Q3VeYKwVnSvufKOf0AAAYFkEUAFCfnPND+21UdLjbDmxnjLBZUTwbV9hXI2tE+9ayDneW6MyQTz63S0fiGWnBRWxYBACoGgRRAMD/z96dx7d53Xe+/zwASQAkuADc910LqcWR5VWyk9hxbGe3ncRNpmmTdJs76XSSTjvdZtrezDSd6dzcuWmTTpd0ErdpWjuOszuLHe+SvNtaSUrc9xUgCRIECBDP/eMhQckSSQAiJZH8vl8vvEABzzk4VBwBX5xzfmdr6j0GLi8U7ow/tHyG6NaYEV3ay7rSPlHDMCjJdTI8FYKa22CqFyZ7r+QQRURELklBVEREtqaeI9Ys4OL+ULCW5tptBkXZlz4jNLwQZj42v2mKFS0F5tUKFpXkOq0lyUszw5oVFRGRa4CCqIiIbD3Tg1aF2PP2h4K1NLc420Ga/dJvf0uBLjt9c8yILgXR6fnpFa8py3UyND0HRU3g8kCP9omKiMjVpyAqIiJbT3x/6K0XPDw4OUfpGoWKgM2zR3RxnKueJZrnYngqxAIGVN2qgkUiInJNUBAVEZGtp/cYZLihZN8FD/f55qjyZq7YLD4jukX2iAJUeFxEFkxGpkNQc8iaKZ4auEIjFBERuTQFURER2Xp6jkHljWBPiz80H40xNDVHZQJBdNPtEY2sHEQrPdbv2++fg5rD1oM6T1RERK4yBVEREdlagj4YPW0tQz3P4OQcMZMtNSPqsDvIsGWsukd0KXj3+YJQvAccuVqeKyIiV52CqIiIbC19L1n31RcWKur1BYHVg+hSoNssQRSssa62NLcsz4lhQJ8/CDa79feiICoiIleZgqiIiGwtPUfAngHl11/wcCJBdCZiFf1xp2+OYkVgBdHVihU50uwUZzutpblgHePi64DpoSs0QhERkYspiIqIyNbScwzKDkD6hdVx+3xBMtJsK54hCtbS3DQjDVfaypV1rzVrzYiCVbCobzGIU7N4nqj2iYqIfczLawAAIABJREFUyFWkICoiIlvH/CwMvXnRsS1gzYhWeFzYbMaKzafD02RnZGMYK19zrcnOyF51jyhY+0TjM6Il+yEjW0FURESuKgVRERHZOvpfgVh0xSC62rJcAH/Yj8fp2ajRbQiP04Mv5Fv1mkqPi6GpOSILMauScNVN0K0gKiIiV4+CqIiIbB09x8CwWUe3vEVfAkHUF/JtviDq8OAP+Ve9psKTScyEocmQ9UD1IRhvg5mxKzBCERGRiymIiojI1tFzxDqixJl7wcNTwQjToWhCQdTr9G7kCNddviufYDRIKBpa8ZoKr7Xntc+/tE9U54mKiMjVpSAqIiJbQ3Qe+l+1ZvveYqlibuVaS3ND/k0XRD0OawZ3tVnRSo/1e/cvBdGyt0F6poKoiIhcNQqiIiKyNQy9CdG5i84PhcSObonGokyGJzddEF0a72r7REtzndhtBn2+xYJF9nRr+bL2iYqIyFWiICoiIltDz1HrvmrlILrajOhkeBJg8+0RXRzvakE0zW6jNNe5PCMKUH0YRk9DcPVCRyIiIhtBQVRERLaGnqOQ3wjuooue6vUF8WZl4Hakrdh8KchtthnRfGc+sHoQhcWzRJeOcIHzzhM9ulFDExERWZGCqIiIbH6xGPS9eMlluWBVzE1kfyhsviC6NCO6VuXcSk8mfb7zZkTLr4c0p/aJiojIVaEgKiIim9/oGQhNXbJQESR2huhmnRHNSs8iw5aRwIxoJqOBMKHIgvVAmgMqboDuF67AKEVERC6kICoiIpvfKvtDowsxBibnqFo8wmQlS0Fus+0RNQwDj9OzZhCtXPz9BybPW55bfQiGT8Lc5EYOUURE5CIKoiIisvn1HoWcCsiruuipoakQCzEzoRlRm2EjNyN31euuRV6nN4Egav3+FyzPrTkEmND74gaOTkRE5GIKoiIisrmZpjUjWn0LGMZFT/clcYZoniMPu82+IcPcSF6nd809ohUea0a0//yCRRU3gD0DerQ8V0REriwFURER2dx8nTAzcslluZDYGaJgzYhutv2hSxKZES3OdpJuN+g7/wiXdJdVtEjniYqIyBWmICoiIpvb0v7QVQoVpdkMSnNX3yPqD/k33f7QJR6nB3949RlRm82gPM9Fv2/uwieqD8HQcQhNb+AIRURELqQgKiIim1vvMXB5oXDnpZ/2BSn3uLDbLl62e77NPiM6F50jGAmuel2lN5N+/1uuqb0NzAXtExURkStKQVRERDa3niNQfesl94eCtUd0rWW5YAVRj2NzzoguBei1ZkUrPJn0+d8yI1p5k7VPtOvZjRqeiIjIRRRERURk85oeAn/3ivtDwZoRXatQUSQWYXp+Gq9r886IAvjm1jpL1IVvdp7ZcHT5wXSXFUa7ntvIIYqIiFxAQVRERDav3qX9obde8ulAKII/GFlzRnQyZJ2j6XVsziC6tLd1rRnRpUDe/9ZZ0drbrfNEg6sHWRERkfWiICoiIptXz1HIcEPJvks+3bdYmCeRirnApp8RnZibWPW6pSNcLjhLFKDmNsCEbh3jIiIiV4aCqIiIbF49x6DyRrCnXfLpZI5uAbb8HtFKz9KM6FuCaPn1kJ4J3c9vyPhERETeSkFUREQ2p6APRk9D1aWX5cLyzN9ae0Q3+4yoK82F0+5cc49ogTsDZ7rt4oJFaRnWPlvtExURkStEQVRERDanniPWfe1tK17S6wuS60on15W+alf+kDWTuFn3iBqGkdBZooZhWJVz37o0F6x9omOtEBjZoFGKiIgsUxAVEZHNqfsFSHNB2YEVL7Eq5rrW7MoX8mE37OQ4ctZzhFeU1+llIrT6HlGASo/r4mJFYAVR0PJcERG5IhRERURkc+p+AapuspaVrqBnYpZqb9aaXflCPjxODzZj874tep3e+MzuaqyzRC8xI1q6Hxy5Ok9URESuiM37jisiIttX0Acjp6Dm8IqXhCIL9PqC1BcmHkQ3M4/TE9/ruprq/EwCoSi+2fkLn7DZrb9P7RMVEZErQEFUREQ2n6VjRmpW3h/aPTFLzIT6Ivea3flD/njl2c0q35mPP+THNM1Vr6tbDOadYzMXP1l7O/i7YbJ3A0YoIiKyTEFUREQ2n+4XrONGVtkf2j5qBa2GBIKoL+TbtIWKlnicHsILYYLRSyy7PU99ofX30bFSEAXo0j5RERHZWAqiIiKy+XS/AJWr7w9tH53BMJaD12r8If+mPbplydKM7lrLcys8mWTYbXSOzV78ZNFuyCzQ8lwREdlwCqIiIrK5zE5Y54eusj8UrCBa4XHhTLevet38wjyBSACPY/PvEYW1g6jdZlBTkHnpGVHDsI7D6XoO1ljiKyIicjkUREVEZHNZOj90lf2hYAXRhgRmQ5eC22afEc135gPgm1u7YFF9oZuOS82IgrU8NzAIE+3rOTwREZELKIiKiMjm0v384v7Qt614yULMpHN8NqH9oUtHnmyFPaIA/vDaR7jUF7rp9QWZj8YufrLundZ9x1PrOTwREZELKIiKiMjmksD+0H6/FbISLVQEm39GNNGluWBVzl2ImfT6LjEr6q0Fbx20/3y9hygiIhKnICoiIpvH7DiMnrH2Ma4i2Yq5wKbfI+pKc+FKcyUURJcr566wPLf+TmvmORpezyGKiIjEKYiKiMjmkcT+UICGwuw1u9wqM6JgVc5NdEYUVjjCBaDhTogEoffF9RyeiIhInIKoiIhsHvHzQ1feHwpWEC3MdpCbmb5ml/6QnzRbGtnpa4fWa53X6Y3veV1NtjOdomwHHaMrzIjW3Aa2dGh/cp1HKCIiYlEQFRGRzaPreai6GeyrB8z2scQq5oI1I+p1eDEMYz1GeFUlOiMK1vLczvEVZkQdbuvvWQWLRERkgyiIiojI5jAzBmMta54fapqmdXRLAvtDwZoRXSr0s9l5nJ7Eg2hRFh2jM5grnRdafweMnILA8DqOUERExKIgKiIim0PnM9b90vEiKxgLhAmEogkHUV/Ih9e5+feHwvKM6Irh8jx1BW6mQ1EmZucvfUHDu6x7zYqKiMgGUBAVEZHNoeMpcHmgdP+qlyVTMResILpVZkS9Ti/RWJSZyApLbs9Tv/j30zG6wrXFeyCrSMe4iIjIhlAQFRGRa59pQufTUPcOsNlXvbR9LPkgupVmRCGxs0Tr45VzVyhYZLNZy3M7n4ZYbN3GKCIiAgqiIiKyGYy1QmBozWW5YM2IZjvSKMp2rHltKBoiGA1umSC6NLObSOXcslwXznQbnSsd4QLWMS7BCRh6c72GKCIiAiQYRA3DuMcwjDbDMNoNw/j9Szz/ScMwxgzDeHPx9qvrP1QREdm2Op627usTC6L1Re6EquAuBbatEkSXfo+J0MSa19psBrUF7pXPEoXl4N+h5bkiIrK+1gyihmHYga8A9wJNwMcMw2i6xKUPm6Z53eLtq+s8ThER2c46noL8BsirWvPSZCrm+sLWEtattEcUEpsRBWt57opLcwHchdae3HYVLBIRkfWVyIzojUC7aZqdpmnOA/8KfHBjhyUiIrIoGoaeI9Z+xTVMhyKMBsKJB9E5K4hulRnRpUCd6BEudYVu+v1BQpGFlS+qvxP6X4bQ9HoMUUREBEgsiJYDfef9uX/xsbd6wDCME4ZhPGoYRuW6jE5ERKTvJYgEE94fCtBQmOAZouGttTTXYXeQlZ6V1IxozISeieDKFzW8C2JRq2iRiIjIOlmvYkU/AGpM09wHPAE8dKmLDMP4dcMwXjUM49WxsbF1emkREdnSOp4GWxrUHF7z0qSPbtliM6Jg/S6J7BEFqF8M7KsWLKq8yTo2p+3H6zE8ERERILEgOgCcP8NZsfhYnGmaE6Zphhf/+FXg+kt1ZJrm35mmedA0zYOFhYWpjFdERLabjqeg4kZw5qx96egMGWk2Kr2ZCXXtC/tIt6WTlZ51uaO8ZnicnoRnROviR7isEkTtadB4N5z9CSxE12OIIiIiCQXRV4BGwzBqDcPIAH4B+P75FxiGUXreHz8AtKzfEEVEZNuanYCh4wlVywVrRrSuIAu7be2KuWDNiHqd3oQq7G4WXqc34T2imRlplOU6Vy9YBLDrPTDnh74X12GEIiIiCQRR0zSjwG8CP8UKmI+YpnnaMIzPG4bxgcXLfsswjNOGYRwHfgv45EYNWEREtpGuZwAzoUJFAO1j1tEtiRoJjlCUWZTa2K5RRa4iRoIjCV9fX+RefWkuWH//9gxoffwyRyciImJJaI+oaZqPm6a5wzTNetM0/2zxsT82TfP7iz//gWmazaZp7jdN852mabZu5KBFRGSb6HgKnLlQ9rY1L50OReiZCNJUuvYS3iUDMwOUuy9Vf2/zKs8uZyo8xcz8GuFyUV2BdYSLaZorX+TIhtq3Q9vjsNp1IiIiCVqvYkUiIiLryzSh4xkrANnsa15+amAKgD3luQl1vxBbYGh2aOsF0cXfZ2BmYI0rLfVFbmbCUYanQ6tfuOs94O+CMX3XLCIil09BVERErk1jbTDdn/D+0KUgujfBIDoaHCUai1KevbWCaIW7AoD+mf6Erl+aQT49sMY5oTvute5bf5Ty2ERERJYoiIqIyLXp7OJxIY13J3T5yYFpyvNceLMyErp+Kaht2RnRQGIzok1lOdgMOLkY5FeUUwplB6zluSIiIpdJQVRERK5NbT+G0v2Qm1hQPNk/mfBsKCwvXV2aQdwqch25ZKVnJbw0NzMjjfpCN6cH1wiiYC3PHXgNpocuc5QiIrLdKYiKiMi1Z3Yc+l6Gne9J6PKpuQjdE0H2ViQXRA0MSrNK1754EzEMg3J3ecJBFKx9tWvOiALsfK91f/YnKY5ORETEoiAqIiLXnrM/BUzYeW9Cl59Ocn8oWEtXi7OKSbenpzLCa1oqQXRkOsxoYI2CRUW7wVOj5bkiInLZFERFROTa0/Y45JRDyb6ELj+ZShDdgke3LFkKoqseyXKePWUJFiwyDGuWuvNZCCd2PIyIiMilKIiKiMi1JRKyzg/dea8VfBJwYmCKCo8LT4KFisAqVrRVg2hFdgVz0Tl8IV9C1zeX52IkUrAIrCC6EIb2Jy5zlCIisp0piIqIyLWl6zmIBBNelgvW0S3JzIbOL8wzFhzbcoWKliR7lqjbkUZtQVb8CJxVVd8KWUVw6tuXM0QREdnmFERFROTa0vY4ZLih5raELp8KRuiZCLIniSA6ODOIibnlzhBdkmwQBdhTlptYELXZofk+OPszCCVwvYiIyCUoiIqIyLXDNK2KrPV3QJojoSanFo8d2ZdkxVzYemeILkkliO4tz2VwKsTETDiBiz9sLc9t/VGqQxQRkW1OQVRERK4dQ29CYCjhY1tgeV/jnrLEg2h/oB/YukE0Mz0Tr9Mb/z0TsTSjfGpwjYJFABU3QF4VnHw01SGKiMg2pyAqIiLXjrYfg2GDxncn3ORk/xSV3uQKFQ3MDJBuS6cosyiVUW4K5e5y+mcSD6LN5Vbl3ISW5xoG7PkwdD4DM2MpjlBERLYzBVEREbl2tD0OlTdDVn7CTU4mWagIrIq5Ze4ybMbWfRssd5czEEh8aW6OM52a/MzEgihYy3PNBTjz3RRHKCIi29nWfQcWEZHNZbIPhk8mVS13MjhPry+5QkWwtc8QXVLuLmd4dpiF2ELCbZrLcxM7wgWguBkKd2t5roiIpERBVERErg2nH7Pud78v4SanBqz9jPvK85J6qYGZgS17dMuSiuwKomaUkeBIwm32lufS75/DPzufYIMHoO9F60sEERGRJCiIiojIteHko1B+PXjrEm5yYmASgD2L+xsTMTM/w1R4asse3bIk1cq5AKcTKVgEsOcB615nioqISJIUREVE5Ooba4PhE7D3I0k1OzVgFSrKy0yuUBFs3Yq5S5ZmfJOpnNtcZgX6hJfneuusLw9OaXmuiIgkR0FURESuvpOPWtVym+9LrtnAVNLLcpcqyW71pbkl7hJshi2pGdG8zAwqva742awJ2fNha2/v2NkURikiItuVgqiIiFxdpgknvwW1t0N2ScLNhqdC9PnmuK4yyf2hge0xI5puS6c4szipIArWeawJV84F2HO/9SXC8W8mOUIREdnOFERFROTqGnwd/F1JL8s91jkOwC31iR/1AtbS3Kz0LHIdyVXa3YzK3eXJB9HyXHomgkzNRRJrkF0CO+6BN74B0QSLHImIyLanICoiIlfXyUfB7oDd70+q2dH2CXJd6TSVJl6oCJaPbjEMI6l2m1GyZ4kC8Rnm13v8iTe6/lMwOwZtP0rqtUREZPtSEBURkasntmBVXN3xbnAmPkNpmiZHOya4pS4fmy25QLkdzhBdUp5dzujcKOGFcMJtrq/2kJFm40j7eOIv1HAn5FbCa19PfpAiIrItKYiKiMjV0/08zIxYBW+S0OebY2ByjlsbkluWa5rmtgqiSwWZBmcGE27jTLdzQ42HF5IJojY7HPgl6HwGJjqSHKWIiGxHCqIiInL1nPwWZGTDjruTanakwwpJtya5P9QX8jEXnaMie2tXzF2SylmiALfWF9A6HGB8JvGZVN72CTDs8PpDSb2WiIhsTwqiIiJydUTDcOYH1t7QdFdSTY92TFCU7aC+0J1Uu+1yhuiSeBBNcp/o4YYCwPp7TlhOKey8F974ZxUtEhGRNSmIiojI1dH2OISnYO8DSTUzTZNjHePcWp+fdMGh7RZECzMLybBlpFQ5N8eZxpFzSSzPBbj+kxAch9YfJtdORES2HQVRERG5Ol75B8irgrp3JtXs3OgM4zPz3FpfkPRLbrcgajNslLnL6J/pT6qd3WZwS30+L7SPY5pm4g3r74DcKnjta0mOVEREthsFURERufLG2qxCRQc/bRW6ScLR9tTODwXoD/TjdXrJTM9Muu1mVe4upz+QXBAFa3nuwOQcvb5g4o1sdrj+l6DrORUtEhGRVSmIiojIlffKP4A9wypwk6SjHRNUel1UepMPk32Bvngl2e2iIruC/kB/cjObwK2L+0STqp4L1v+mtjR4+e+TayciItuKgqiIiFxZ4Rk4/i/QfB9kJbe8diFm8mLnBIdSWJa7EFvg9MRpdufvTrrtZrbbu5tAJEBvoDepdnUFWZTmOjnankTBIoDsEus4ntcfgqAvubYiIrJtKIiKiMiVdfIRCE/DDb+adNPTg1NMh6IpLcvtnOpkNjLLvsJ9SbfdzJZ+3+Njx5NqZxgGhxoKONIxTiyW3Gwqh/4DRIKaFRURkRUpiIqIyJVjmtay3JK9UHFD0s2XjhNJJYieGDsBwP7C/Um33czq8+pxp7vjv38yDjXkMxmMcGZoOrmGxU3QeDe8/Lcwn8QeUxER2TYUREVE5MrpewlGTlmzoUkevQJWEG0sclOU7Uy67fGx4+Q6cqnKrkq67WZmM2zsKdiT9IwoEF8CfSTZfaIAhz8HwQl44xvJtxURkS1PQVRERK6cl/8eHDmw9yNJNw1HF3ily8ehhuT3h4I1I7qvYF/SZ49uBfsK93HWf5ZgJLnZyaIcJzuK3ckXLAKovgUqb4KjfwULkeTbi4jIlqYgKiIiV8bMKJz5Hlz3ccjISrr582fHmYss8PadhUm3nZ6fpmOqY9sty12yv3A/MTPG6YnTSbe9tb6AV7p9hKMLyb/w4c/BVC+c/k7ybUVEZEtTEBURkSvj5b+DWAQO/kpKzX90cohcVzqHU5gRPTV2CmDbFSpasq8gtYJFYJ0nGorEeKkzhQq4jXdD4S448iVrf7CIiMgiBVEREdl4QR+8+DfQ9EEo3JF081BkgSfOjHBPcwnp9uTfuo6PH8fAYG/B3qTbbgV5zjxqcmpSKlh0uLGAHGca33ljIPkXttmsCrojp+DcE8m3FxGRLUtBVERENt6Lfw3zAXj776XU/NmzY8yEo7xvf2lK7U+MnbCqx2a4U2q/Fewr3MeJsROYSc5MOtPtvG9/GT85NcxMOJr8C+/5MORWwjNf0KyoiIjEKYiKiMjGOn82tLg5pS5+eGIIb1YGt9Qlf2xLzIxxYuzEtt0fumRfwT4mQhMMzCQ/s/nAgXLmIgv8+ORQ8i+clgHv/EMYfEN7RUVEJE5BVERENtaxryzOhv5+Ss3n5hf4ecsI9+wpIS2FZbk90z1Mz09v+yC6v8j6/VNZnnugykNNfiaPvZ7C8lyAfQ9CUTP8/POqoCsiIoCCqIiIbKSgD176W2j6EBQ3pdTF022jBOcXeN/e1JflwvYtVLSkIa8BV5qLE+PJB1HDMLj/QAXHOifo9yd3BAwANju860/B3wWvfT359iIisuUoiIqIyMaJz4amtjcU4EcnhihwZ3BjrTel9sfHjpOdnk1tbm3KY9gK0mxpNOc3c3w0+cq5APe9rRyA76ZStAig8S6oPgzP/g8IB1LrQ0REtgwFURER2RjrMBsanI/y89YR7t1TmtKyXLBmRPcW7sVm6C1vf+F+Wn2thKKhpNtWejO5qdbLt18fSLrgEQCGAXf93zA7Bke/nHx7ERHZUvSuLCIiG+O5/wfmZy5rNvTnLaOEIjHeuy+1ZbmzkVnOTZ7b9styl+wr3EfUjNLia0mp/QPXV9A1PsvrvZOpDaDioFW06uhfwcxoan2IiMiWoCAqIiLrb/gkvPQ3cP0vpzwbCtay3KJsBzfUpLYs9/T4aWJmbNsXKlqyFMhTKVgEcO+eEpzpNh57vT/1QdzxxxANwdN/lnofIiKy6SmIiojI+orF4Ie/DS4P3PknKXcTCEV4um2U9+wtxW4zUupjqTDP3oK9KY9jKylwFVDuLuf4WGr7RLOd6dzTXMIPjg8SiiykOIgGuPHX4bWHoO+V1PoQEZFNT0FURETW1xv/BP0vw7v/K2SmNpMJ8Mir/YSjsXiRnFQcHz1OTU4NuY7clPvYavYV7uP42PHU9nkC9x+oYDoU5UcnUjhTdMkdfwTZpfDDz+o4FxGRbUpBVERE1s/sBDz5J1B9CPZ/LOVuFmImXz/axcFqD/sr81LqYyo8xYtDL3JT6U0pj2MruqnkJkaDoykd4wJwuKGA5rIc/teTZwlHU5wVdWTDe/4njJyCF/86tT5ERGRTUxAVEZH188QfW0dzvPeLVpXUVLs5M0yfb45fvS31I1e+3/F9QgshPrzjwyn3sRXdW3svWelZPNz6cErtbTaD/3TPLvr9c3zzpd7UB7L7fbDzPfD0n4O/J/V+RERkU1IQFRGR9dFzFN78Btzym1C0+7K6+urzXVR6XdzVVJJS+5gZ4+G2h9lfuJ9d3l2XNZatJjM9k/fXvZ+fdP8Ef8ifUh+3NxZwS10+X36qnZlwNPXB3PsXYNjg8d+BFJcKi4jI5qQgKiIily/og8d+HfKq4e3/6bK6erNvkld7/Hzq1tqUixS9NPQSPdM9PLjzwcsay1b14M4HicQifKf9Oym1NwyD3793FxOz8/z9c52pDySv0toveu5ncOa7qfcjIiKbjoKoiIhcnlgMvvt/QWAYPvI1yMi6rO7+4YUush1pfPSGypT7eLjtYTwOD++uefdljWWravA0cLD4II+0PcJCLLV9nvsr83jP3hK++nwnY4Fw6oO58TegZB/86Hd0tqiIyDaiICoiIpfn6F/C2Z/A3V+A8usvq6vByTkePznEgzdU4nakpdTH8OwwT/c9zX2N9+GwOy5rPFvZg7seZGBmgCODR1Lu43fevZNQNMaXnzqX+kDsaXD/38H8DHzvM1qiKyKyTSiIiohI6nqOwc8/D00fght/7bK7e+hYN6Zp8slDNSn38ejZRzFNk4/s+Mhlj2cru7PyTgpcBTzcllrRIoC6QjcfPVjJN1/upWdiNvXBFO2Guz5vLdF95aup9yMiIpuGgqiIiKRmdhwe/RR4quEDf3VZVXIBZsNRvvlSL/fuKaXCk5lSH5FYhG+f+za3VdxGRXbFZY1nq0u3p/NA4wM83/88/YH+lPv57LsacaTZ+ff/8gahSIrHuQDc+OvQ8C742X+G0dbU+xERkU1BQVRERJI3H4R//TdWkaKPPATOnMvu8m+e7SAQivJrt9el3MdTvU8xPjeuIkUJ+vCOD2MYBo+efTTlPopznPyvB6/jRP8Uf/jYScxUl9YaBnzwr609xo/9KkQvY9+piIhc8xRERUQkOdF5eOQT0PcS3Pc3ULrvsrvs8wX52+c6+dB1ZVxXmZdSHzEzxjfOfINydzmHyg5d9pi2g5KsEt5R8Q4eO/cYgflAyv3c1VTMb9+1g8feGOAfXuhKfUDZxfDBr8DwSWvJt4iIbFkKoiIikrjYAnzn16H9SXj/l2DP/evS7Rceb8FuGPzevamf+fml17/Em2Nv8mt7fw27zb4u49oOfmXvrxCYD/C7z/4u0VjqZ4L+5jsbuLu5mC883sIL58ZTH9DOe+GGX4VjX4bTqR0vIyIi1z4FURERSYxpwg8/Z4WDu/4rXP/L69Lt0Y5xfnxqmM+8s57SXFdKfXyv/Xv8n1P/h4/u+Cj3N65PON4u9hXu449u/iOODB7hi69+MeV+bDaDL370OhqK3Pzmv7xO1/hlFC+6+wtQcSN89zMwcib1fkRE5JqlICoiImuLxeAnvw+vPwS3/Uc49Fvr0m10Icbnf3CGCo+LX70ttb2hr4+8zp8e+1NuKrmJ37/p9zEus2jSdvThHR/mE02f4Bst3+CRtkdS7sftSOPvf+kgAPf99RGePzeWWkdpDvjoP4LDDf/6cZjzpzwmERG5NimIiojI6iJz8K1fhpf+Bm7+DNzxX9at6395pY/W4QB/9J7dONOTX047MDPAZ5/+LOXucr74ji+Sbktft7FtN//x+v/IbeW38ecv/TkvDb2Ucj/V+Vl8998dojjbyS//n5f5m2c7UitglFNqhdGpfvj2r1nLwkVEZMtQEBURkZXNjsND74eWH8Ddfw73fOGyj2lZMhmc5//9WRu31OVzz56SpNu/MvwKv/HEbxA1o/zVHX9FriN3Xca1Xdltdv7i9r+gOqeazz3zOX7TAt2IAAAgAElEQVTQ8QNiZiylvmoKsnjs393KvXtL+e8/buUz33ydmXAK+0+rboZ7/we0PwFP/1lKYxERkWuTgqiIiFzaeDt89U6rgulH/xFu+Xfr1rVpmvzOt04wE47yx+9vSmo5bX+gn99+5rf59E8/zfzCPH/5zr+kNrd23ca2nbkz3HzlXV+hMruSP3zhD/nE45/gxNiJlPrKcqTx5Y+9jT+4dxc/OTXM7X/xNH/9THvygfTgp+HAL8HzX4TX/zGlsYiIyLXHSPm8r8t08OBB89VXX70qry0iIqswTTj+r/D471p79T7+MFQcXNeX+PvnOvmzx1v44/c18enDa4fImBmjxdfCT7p+wjdbvondZufTez7NJ5s/iTPNua5jE+vv+/sd3+dLr3+J8blx3lv3Xu5ruI8DRQdItye//PmNXj//35PnePbsGHmZ6fzKoVp+8eZqPFkZiXWwEIFvPgidz8DH/gV23J30GERE5MozDOM10zQv+SFCQVRERJbN+eGHvw2nH4OqW+H+v4W8qnV9idd6/Dz4t8d41+5i/vcvHrjkbKhpmgzNDnF64jQvDLzA8/3PMzY3hoHBe+vey2cPfJbirOJ1HZdcbDYyy1dPfpV/OvNPhBfCuNPd3FJ2C7dX3M6+wn1UZ1cndVTOm32TfPmpczzZMorNgANVHu7YXcQdu4rYWZy9+sx4OABffy+Mn4Nf/iFUXL8Ov6GIiGwkBVEREVlb1/PwnX8LM8Pwjj+Aw5+DdT6P0zc7z3v/8nnS7TZ+8O8P43bYGA4O0zPdQ890D91T3Zz1n6XN30ZgPgCAO93NofJD3F5xO4fLD+N1etd1TLK2YCTIi0Mv8lz/czzX/xxjc1Y1XKfdSUNeAzu9O6nLraMmt4banFpK3aWk2dJW7K9laJofnxzi562jnB6cBiAvM50dxdnsLM5mR0k2NfmZFGU7KXBn4MnMwGYzIDAC/3AXzM/ArzwB+fVX5PcXEZHUKIiKiMjKpgbgiT+GU4+Ctw7u/+q6zjaZpsnY3Bhdk918/qfP0jnZw807Y/jnB+kL9DEfm49f60pz0ehpZKdnJ7u8u9jp3UlTfpOq4V5DTNPk3OQ5WiZaaPO30eZro83fxlR4Kn5Nui2dyuxKanJqqMmtueDe4/Rc0N/IdIinW0c53j/F2ZEAZ4cDBN6yj9RuMyhwZ1CY7WCPY4z/MvJZImlZ/PTGr5FZUEVhtoOibAdlea6Uqi+LiMjGUBAVEZGLRUJw7MtWEZjYAhz6D3D4s5CRlVJ3oWiInukeuqa66JzqpGuqi+7pbnqme5iLzsWvsxtp1ORUU5VTRU1ODVU5VVTnVFOdU02hq1DngG5CpmkyGZ6ke7qb7qluuqa76J7qpnu6m75AH9HYcrDMdeRawfT8kLr430GGPcNalj0Vos8XZHxmnrFAiLGZMGMB6zY+M0/h1Cm+FPlTJswcfmH+PzNMfrz/omwHFR4XVd5MGouz2VWSza7SHMpynfpvS0TkClMQFRGRZZEQvPkNeOFLMNULu98P7/5v4KlZs6lpmvjDfjonO+ma7qJravk2ODOIifWeYmBQ7i6nJreG6pxqWvscPH/G4P49+/nzD9xGmn3lZZuytURjUQZnBume7r7gy4nuqe74El8Am2GjLKvsgnC69HNRZtFFITLW+zLGN+4j4irgxJ3/TG80j37/HP3+IP3+ObrHZxmcCsWvz3GmcaDaww01Xm6o8bKvIlezpyIiG0xBVEREYH4WXv0aHP0rax9oxQ1wx3+Gundc8vKp8BTtk+2c85+L33dMdVywBNNpd8b3Bdbm1VKbW0ttTi3VOdU405yYpslf/LSN//1MB5+4uZrPf7BZs1ISNzM/Y82iT3fFw+mlZtEz0zKpza2lIa+BRk8jjXmNNHgaKBzvwvjGA+AuhE/+CHLKLuh/OhTh7HCA1uEApwenebXbx7nRGQAy0mzcVOvlzl1F3LGrmKr8zCv6u4uIbAcKoiIi29nYWXjt63D8m1ZV3Jrb4PbfhdrbwTAIRUN0TnVeEDjPTZ5jNDga7yI7PZtGTyN1eXXU5dZRm1tLXW4dJVkl2IxLH0ltmiZf/NlZvvx0Ox+/qYr/9sE9VsEZkTXEzBijwdH4Ut/u6W46Jjton2xnfG48fl2uI5cGZzGNAydotGfS+K4/p77qMDkZOSv27Z+d57UeP8c6J3i6bZTOsVkAGorcvGdvKR/YX0ZDkXvDf0cRke1AQVREZLuZn4W2H1sBtPt5sKVh7nwPowc+Tqszk1ZfK23+Ns75z9Eb6CVmxgDIsGVQn1dPo6eRhryG+AxUcWZxUjOZ4egCf/ajFv7xWA+/cEMlX7hvr0KorAt/yH/BFybt/nbafW3MLCzPoJZkldCY1xgveLXLu4vK7MpLfmnSPT7LU62jPNkywoudE8RMaCrN4QPXlXHf28opztE5tSIiqVIQFRHZDqJhaH8STn2baNuP6WGeFk8ZbWXNtDoctE114Q/745dXZVex07szvsyxIa+BquyqpM6FvJR+f5DP/PPrHO+f4tduq+UP7t2tECobyjRNhvuOcO77/5ZzsSDtDbdzNhqgc7KTqGkVSspMy2SHZwe7vLvitwZPAw67I97P6HSIH54Y4vvHB3mzbxK7zeDOXUV87KYqbm8sxK7/jkVEkqIgKiKyVc1OEGz7EWfbvkfr6Ju02mK0OV2cy0gnzPIsZ6PnwtmhHZ4dZKWnVh13NU+1jvC5h48Ti5n8z4/s4549pev+GiIrmh6Cf/4IjJ6BD/wl8/s+SsdkB62+1vitzd/GbMRajms37NTm1saD6U7vTnZ5dpHnzKNrfJZ/faWXR1/tZ2J2nvI8F79wQyUfvaFSs6QiIglSEBUR2SKmZ8doafsup3uf5Yy/jbZogJ70NMzFZbO5aZnsKmhml7cpHjprcms2/BzOsUCY//XkWb75Ui9NpTn89b85QE3B+gddkTWFpuGRX4LOp+HG37AqQqdlxJ+OmTEGAgO0+heDqa+NFl/LBXuiS7JK2OXZxe783ezI28XYRBE/fGOGI+0T2G0Gd+wq4uOaJRURWZOCqIjIJhSYD9AyeoIzXU9yevgVzswO0Gssn8dYZtrYnVnGztIb2FX1DnYXNCW9l/Nyzc0v8A8vdPK/n+kgHI3xiVuq+b17dulYDLm6FiLwxJ/Ai1+B8uvhI1+HvKpVm/hCPtp8bbT52mj1t9Iy0UL3dHd8/3Shq5Ca7B3MzZTS2pPL5GQx5dnF/JubqvnowQry3Y5V+xcR2Y4UREVErnEz8zO0+Fo4M3aS0wNHOeNroScyHX++LBKlyeakObeBporDNO26j7zc1T9Yb6SJmTDffXOQv3+uk+HpEHc3F/N79+yirlDVRuUacub78L3PgGGD+/8OdtydVPNgJEirr5UzE2fit86pzvh5uWlmLqGZUghXcn3pXj518DB3NDboiCIRkUUKoiIi15DZyCwtEy2cnjjNmfHTnBl9g+7gcPz5kmiU5vA8Tel5NBfuY3fdu/E23A2Z3qs4aqsS7jNtY3z7tX6eah0lGjO5vtrD792zixtrr+7YRFY00QHf+mUYPgkHfgnu+jy4PCl3d344PT1xmjdHTtE/2wOL4dQWy6U+ZyfvqDnA/qI9NOU3UZhZuE6/jIjI5qIgKiJylQQjQWumc/FD65nxM3RPd8dnVIoXYjSFQjTPz9PkKKCp/Fby6+60zvp0X70Pr9GFGCOBMGcGp3mtx89rPT6O908xH41R4HZw/4FyHjhQwc6S7Ks2RpGERebg6S/Asa9AZj7c+9+h+X5Yp5nLYCTIGyOn+dbJYxzte5MZerBnjIFh/f+8yFVEU34TTflNNBc005TfRIGrYF1eW0TkWqYgKiJyBQQjQdr8bZwePx0Pnl1TXfHQWUQ6TaE5muZmrBlPVzEF1bdD7duh5jDkpF5hNroQY2gqhG92nqm5CJNzEabmIoQjC0QWTKILMSILMeYXTCKLP0cWYsxHL/zzxOw8Q5MhRgMhYotvD+l2g+ayXK6v9nC4oYDbGgtIs198HqPINW/oBPzgt2DwDWh8tzU7WrR7XV/CNE1e6/HztWOtPNHxBmZ6PyWFY2RkDjIW7l/+92ApnBY00ZyvcCoiW5OCqIjIOpuLztHma7NmOc/bO7ZU2KQgLZNmM53mqXGaZvw0zc9T6C63Zjprb7Pu8ypTeu3ByTle7fFzvG+SzrEZuieC9PmCRGOr/3tuGJBut5Fht5FuN0i32xZv1s9pdhuezHRKc12U5zkpzXNRX+hmX0Wuig/J1hFbgJf/Dp76bzA/C833wdv/07oHUoDxmTAPv9LHN1/qZWByjqJckzv2Raks9dE3e5bT46fpme5ZDqeZVjhdCqYKpyKy2SmIiohchrVCZ35GDs1pedYS24kemgI+ihYWwF2yHDprbwNPbUpLAaeCEZ5qG+GZtjFe7fYzMDkHgDPdRm2Bm9qCTGrys6jOz6TA7SDXlU5eZjo5rnSc6fbF4GnTMRMi5wv64NiX4aW/XQ6kN/46VN28bkt2lyzETJ5uHeWfXuzh2bNj2G0GdzcX84s3V7Ov0nnBntMzE2cuCKfFmcUXhdN8V/66jk9EZKMoiIqIJGjN0On00pRZStOCjaapUZoHWykKBzDACprVt1q3qlvAW5fyB9qxQJgfnRjkZ2dGeKnLx0LMpMDt4KZaLwdrPBys9rK7NFtLZEUuVzyQ/h3MB8BTA/s/Bvt/wfp5nfVMzPLPL/XyyKt9TAYj1Bdm8Ymbq7n/+gpynNZ5v/Eq2udV6+2e7o73UZJVQpP3wj2nXqcKhonItUdBVETkEpb2dL71aIbl0JlPU14jTfYsmoIzNI92UDR4EiMWAQwobl4OnVW3XNYeT7D2eT7TNsbDr/bxVOsoCzGThiI3724q5q6mYvZX5GHTrKbIxgjPQMsP4Pg3oet5wITC3db+7ZrDUH1oXQuIhSIL/PDEEP/0Yg/H+yZxpdv50NvK+cRNVTQVOSAcgPkZMGMQW7DC6XQXZ2Z6OD3dzZnJc/TM9Mf7K80qXS6ItDh76nGmXh1YRGQ9KIiKyLb31tB5evw0XdNdF4TO5vxmmrLKaYos0OQfoKj/dYyxNqsDWzqUH7ACZ/WtUHkTuPLWZWydYzM88mo/3369n7FAmAK3gwcOlPORgxU0FKkqrcgVN9kHp74NXc9B74sQmbUezy4Fbz14ayG/HtzF4MgGR451b08H07TCoxmDaAhC0xBevIXech8OQDjAXMBPMODHNj9DFkEyjIWEhhkwDFqzcjntzuWMI4MzNpMe5uPPl7oKaS7cR1N+czyc5jnX598tEZFEKIiKyLZyqUPozw+dBa4Ca+bAu4smezZNAT9Fgycw+l6C6QGrE0cOVN64PNtZfgDSXes2xtlwlMdPDvHIq3280u3HbjN4585CPnqwknfuKiJdS25Frg0LERh8E3qOwFgb+DrA1wmzY6n1Z0uz/n1xLoZXR+7ifTbhNDctPpNXhqL0ztqxZ2RxY0MRhxqLyM10Wu2jYSvgRkPW3tagzxrL7BjMjDA92UOrLcqZjAxOOzI443TSm7ZcbKzMVUhz4X6aCppp8jaxO3+3Zk5FZMMoiIrIljU+N07LRAtt/rb4fe90b7zQRzx05jfR7G2iyXBR1P86dD4NvS9BeMrqKLt0OXRW3wJFTWBb30qxpmnyRt8kj7zSxw+ODzI7v0BdQRYfOVjJAwfKKcpxruvricgGCk1DcGJ5ZjM0DbEoGLbFmwFpDitoOnOWZ03TXWvuHY/FTI50jPPQ0R5+3jqC3TC4d28pnzpUw4GqNUKjacLMiBWWx8/ByCmmh96gxX+WM/YYpzMyOOPKpM++PIbizGJ2eXex07uT3d7d7PTupMJdgbHORZtEZPtREBWRTS9mxuid7qXV13rBbSI0Eb+m3F0e/xC1y7uLpvwmimwO6HwW2p+Ejqdgqs+6OL/B2ve1FD7zqta9UuaSoak5vvfmIN9+rZ9zozO40u28d18pD95QycFqjz7siciKeieCPHSsm0de6SMQjrK/Mo9P3VrDe/aWkpGWxMqJ2AKMn4XeY9B9hKneI7TM+2jLyKDF7aEtM4vOhSCxxS/x3Onu+L+lS7f63HrS7ekb84uKyJakICoim0ooGqJjsoMWX0s8cJ71n2Uuah1bkmakUZ9XH/9wtNO7k53eneRk5EAsBkNvQPtTVvjsfwXMBcjIhrq3Q8OdUH8neKo39HeYCkb4yekhvvPGAC91+TBNeFtVHg8erOR9+8twO9I29PVFZGuZDUf59uv9fP1oN51jsxRmO/jFm6r5+E1VFGY7ku/QNK1Z046noP3n0PUsoegc7S43reV7ac0rptUMc3aqY/nfXlsaDXkN7PTsZHf+bnZ6rH97szO0l11ELk1BVESuSQuxBfpn+jnnP2fdJq373kBvfD9nVnrWBR96dnl3UZ9XT4Y9Y7mjwPDyh6mOp2DOZz1eeh00vMsKnxU3WIVENlCfL8iTLSM8cWaEl7t8RGMmdQVZfPC6cj54XRk1BVkb+voisvXFYibPt4/ztSNdPNM2Robdxvv2l/KpW2vZW5GbesfRsLUPtvVxq3rwzDDYM1iofTu99Ydp85TTMttPm6+NFl8LvpAv3rQkq4TGvEZ2eHbQ6Gmk0dNIbU6tZk9FREFURK6+ibmJeNBcunWc9027gUFldmX8Q0xjXiO7vbspzy7HZrxl+Vk0bFWy7Pi5NfM5ctJ6PKsI6u+wwmf9OyGrYEN/p/GZMC91+nixc4IXOyc4NzoDQGORm7uairm7uYR9FblaeisiG6JzbIaHjnbz6Gv9zM4vcLDawycP1XB3c8nlFTyLxazVJC3ft26Tvda+1+pD0PRBzF3vZzzNTouvhbP+s/EvErumuojGooC1cqUmt4ZGz2JAzbP+bS/NKtW/iSLbiIKoiFwxs5FZuqa6OOc/Z31AWQyf53977nV642Fz6Rv0utw6MtMzV+54omNxxvPn1hl/kVnrSJWqm5fDZ/EesG1MtVn/7DwtQ9OcGpzi1MA0pwam6By3jnTIzLBzQ42Xww0F3NVUrJlPEbmipkMRvvVqPw8d7abXF6Q018kv3lzNx26swpuVsXYHqzFNGD4BZxZD6fhZwLCOsWr6EDR9ALJLAIgsROie7r7o3/+h2aF4d+50Nw15DRd86VifV6/KvSJblIKoiKwr0zTxhXx0TnXSNdVF51QnnZOddE51MhIciV/ntDsv/MCx+KEj35W/9ouEA1bgbH/SCp/+butxT+3yPs/a26wqlJdpPhpjYjbM6HSYsUCYsRnr5x7fLF3j1m0yGIlfX57norksh+uq8rilLp895bk6bkVErrqFmMkzbaN87Ug3L7SP40iz8aHryvnkoRp2l+asz4uMtsDp78KZ78JYK2BYXwguhdKcsouaBOYDtE+2LwfUxRnUwHwgfo3H4aE2t5a6vDrqcpdvJVklmkEV2cQUREUkJTEzxuDM4AWBc+l+aunYEyAzLdP6AJFbR11eHbW5tTTkNVDhrsCe6BEo0XkYeA26nrWq3Pa/ArEIpGdB7e2L4fMO6xD5BJimydRcxAqW54XLsZlw/LHRQIixQBj/eSHzfKW5TmoLsuK3nSXZNJflXv4Mg4jIBjs3EuDrR7t57PUB5iIL3FTr5VOHannX7iLS1uuLs9FWOPM9K5SOnrEeq7wZmj8Euz8AueUrNjVNk5HgCO2T7fEvMpfeXybDk/HrXGkuanNrl99jFm+VOZWk27QHVeRapyAqIquaCk/RM91Dz3QPfYE+uqe76Zrqonuqm9BCKH6d1+m96MNAXV4dxZnFyX9jHVuA4ZNW8Ox6DnqOWcttMaB0P9S9wwqflTdD2sXBLzgfZXByjoHJEIOTc4s/zzE0GWJwyrqfX4hd1C4jzUZRtoOibAeFSze3k6IcB4Xu5ccK3I7kjkYQEbkGTQUjPPxqLw8d7WFgco7CbAcfuq6MB66vYFfJOs2SAoydtQLp6e/C6GnrsYoboemDsONu68isBN8nfCHfReG0c6qT4dnh+DVpRhqVOZXU5dZRk1NDVU4VVdlVVOVUUegq1CyqyDVCQVRE4mGzN9BL73RvPHT2TPcwPT8dv87AoDSrlNq8twTO3DrynHmpDyASgqE3oe8l6HsZul+A0OK33gU7raNVam+3zvZ0Le8Vmg5FODcyw7mRAGdHZjg3GuDsSICR6fAF3dsMKMlxUpbnoizPRWmuk6Ic52LQtMJlUY6DbEeaPqCIyLYTXYjxVOsoj77Wz1Oto0RjJk2lOXzgujLu3VNCdf467m0fP7cYSr+3XEzOUwuN77ZuNYcg3ZV0t8FI8IJguhRW+wP9RM1o/DpXmiseSquyq6jOqaYyu5LqnGoKXAV6DxC5ghRERbaBmBljYm6CgZkB+mf66ZvuoyfQQ+90L72B3guW0i6FzfO/QV56sy7PLsdhT+FMurcKDFuBcyl4Dr0JC/PWc946qLrVCp81t0FOKYFQhHOjVuA8NzLD2cWfh6aWZ2Sd6TYaitzsKMqmvshNeZ6Lco8VPIuzHeu33ExEZAvzzc7z/TcHeOyNAU70W+8Nu0tzuKe5hDt3F9FUmoPNtk5hzd8N556wbl3PQXQO0lzWF4873g0Nd132uc7RWJSh2aGL3vd6p3sTCqkV2RWUu8spyiwizaYznkXWk4KoyBZwftAcnBlkcHZw+efF23xsPn79FQmbS2YnYOgNGDoOg29aoXOy13rO7oCyt0HVTVB5E7NFBzg36+LsSGB5lnMkwOB5gdORthg4i7NpLHbTWJTNjmI3FZ5M7Ov14UhEROj3B/nJqWF+enqYV3v8mCZ4szK4tT6f2xoLuLkunypv5vrMIkbmoPsInPsZnPvpchG6vCqoPmzNlFYfAk9Nwst415JMSLUbdkqySihzl1GWVUa5u9z62V2moCqSIgVRkU0gshBhbG6MkeAIQzNDawZNsPZslmUtv0me/4ZZkV2xvmETrDL+04NWUYqhNxdD53GY6lu+xlMLpfsJFh2gJ3MPp2LVnPNF4qFzYHIufmlGmo36Qjc7ihdD52L4rPQqcIqIXGmjgRAvnBvnhXPjPN8+zljA2gLhzcpgf0Uu11V62FuRQ0NhNhUe1+XNmpomTLRbx3L1vAA9RyE4YT2XVQjlB6Hieuu+dD9ketfhN7zQUkhdeq8dmBm44Oex4Bgmy5+Tzw+qpVmly++7WWUUZxVTlFmEKy35JcciW5mCqMhVNhuZZWR2hJHgCKPB0eX78x6bCE1c1G61oFmaVbr6uZuXPegJK3COtpx33wLnLfE1vXUE8/cwlLmLNqOeV+crOeUz6BibxTe7HJoz0mzUFWSxo9ia2WwszmZHcTZVCpwiItck0zQ5OzLDqz0+3uyd5M2+SdrHZlj62OhMt1FX4Ka2MIuyXCeluS7K8pwU5zjxZGbgycogx5nEnvxYDMbboOcI9L9mVU6fOLf8fHYZlOyB4marrkBBo1VF3bVx54/OL8wzPDt8QTgdnB1cMagC5GTkxENpcWYxxZnLPxdlFlGSVUJORo72qcq2oSAqsgFM0yQQCTAxN8H43DgToQkm5qzbBWEzOMJsZPai9rmO3BXfqJa+bd3QoAkQ9IGvC/xd4Otc/nmiHWbH4pctOPKYyWlk2FlLt62KU5FyXpgp5vSE7YLKtAXuDOoK3dQXZlFf6I7fyj0uBU4RkU1uOhTh7HCA9tEZ2kdnODc6Q68vyODkHOHoxVXK7TaDHGcarnQ7zgw7rvTFW4YdZ/rynzPSbDjSbDjSbTjSlv/sNmconmkhf+Ysnuk2cqbayJzuwHbeclozswAjrxJyyiG71DrHdOmWXQY5pZCxjoWYznN+UB0Njsbf88//stkX8l0UVh12xwXv+cWZxRS4Csh35Vs3p3Wf58jDZqj2gWxuCqIiCYqZMQLzVriMB8vQYtB8y2MTcxNEYhefP2kzbBS4Ci4MmG/5drQws3Bjl++YJoQDEBiyltIGhiEwCNND1mNT/VbgDE1d0GzOWYzPUcGgrYS2hXJeC5VyLFDEsJkHWEEyw26jwuuirsBNfdH5gTOLvEydrykist2Ypok/GGFwco7RQAj/bAR/cB5/cJ6puQihSIy5yAKh+QXmIou3+QVCiz/PR2OEF28LsdU/l6YTpdIYpc4YonbxVmn3UWr4KWaCbC7+4jdqd7HgyMV05mLP9JCW6cFwecCVB848cOaCww2ObMhYvD//5ww32FPbG7q07WY0OMpwcJjR2dGLvqweDY5e8vOE3bDjcXqskOq8MKR6nd74nwtcBeQ58hI/t1vkClIQlW0pshBhMjwZv02Fp1b9eSo8xVR4igVz4aK+7Ib9gn/0r+ibQSwG4WkrNC7dz01ae2mC49asZnDiwtvsOMzPXPx3kp7DjKMQv72QPoo5O1/A/9/evYZIlt51HP/+z71uXd1z6d3N7uyakEUSJChKVBARNWYRcRM0JEFE8UUIEoxI8BZINBJUBDH4ymgCRoIiRGVfrCRKAioY3RhizG3DZjXJJjOzOzPbt7qey98X53R1dU93z8zuTlXN7u8Dh+d5znm6+mH3THX9znnOU58fneUrk3N8wzeZUAfJc92EC2faPHCmzf1n2lxoyvvPtrmrl71wqymKiIjMKcqKaVkxyQ/KSVHOguokL9mbFOyOC3bGObvjgt1ZWTAZ7mB7l0gGl2iNL7PpV9mwXfoM6FuzMWAjGNK3AS0f3XhQUK/0m/bqwJp0IV07El6bffvhtX0WOufqsn2mDrwnTMd1d3amO8de8D7uYvhxodUwekmP9XSd9XSdfto/VK6n6/SzPhvpxqH9WZQ9n/9dIjd0WhDV0l+ysoqqYG+6x+50l918l93pLnvTPXamO+zlzf5mm2/vTHfYmmwdOx12Xxqmh96cX7n+yll9PV0/ND3mXOsc/bR/8vSYqoJyAsW4/q7M0TegmNRL1Bdz+/P6PwgAAA7pSURBVIvxQbuYa+ej+u7lfsgcHwmdk13g5AtGHneoWmeYJOsMwj5byau5Evf4ZrHOk5M1Hh90eHKyxmXfYDxOYffgruYDm3XAfPNc0Lyw0aaT6q1BREQWLwoDojDghZhg4+5sj3Ku7E25sjfh6t6Ux/cmXN4Zc2l7zLe3RzyztcfezrPE5ZAOI7qM6NqYfjDinqzkrixnM51yNp6yHk7pByM6jMmKAcnk29h0r/kbvlf/3T9JEEHrzOFw2q7r1jlHv3OOfmeTV3Q34e5Xnhhc5x8Lmg+nz06eZWt8cHH9yugKX9v6GluTLYbF8MRhZWE2+zzUTbr0kh69uEcv6dFNuqwla3Tj7rHtXtIjCTUTSp67m7ojamYPAR8AQuAv3P0PjhxPgY8A3wtcBd7s7v932mvqjuiLi7szKScMiyHDfHioHOWjw/tP6TPIB3X4zHcZnfaG3ujEnfrNMW7eHJs3yfV0nfVkjfUwox9mrFvMugWse0jfoVXm9R3D6QCmwyPhcC4gHhsc50Pl6OC7MZ+rIKqvnmb9+mpq1sezNcpkjXHQZRR02aHFtbLFlTzj8jTl2+OE/x21+Nog5ak9mB55NqebRrPv2Dxa3rfe4lw31V1NERERDqYWX9weNQF1zKXtERe3x1zcGnNpZ3zsc7CBwbluyj39epGml61FXOhU3NvKuSceshnucSbYJZ08e3jG0mwmU1M/7mJzmNSrB3fOQ/cu6J6HziZ0N+f2NfXWxqlfdzMtp4dmgh2aETZu6tPt2cX/vby56D/du+751qPSMJ0F007coR23aUftQ2Un7tT1/f3Nsdn+pk8WZSRBooWcXmSe19RcMwuBrwKvA54CHgPe6u5fmuvzy8Br3P3tZvYW4I3u/ubTXldB9PZwd4qqIK/y2TYtp3W9zJlUEybFhHE5ZlJMmJQH27gY1+WRY/PtcTlmWk5n+8bFeBYuK79+oYKTHHqTijLaYUYrTGhbTC9M6QYxPYvpWUTPja47vcrpVSW9IqdbTOnmE8J8OBcoB/UVyf12Obn5/3AWQtyCKIUom9tSPM7wMMPDlCpM8CijDFOqIKUKM8ogoQxTyiClCFLKICEPUoogoSAhDxJGVczQYwZlxG4RMShDdsuI3SJkkNdXbLeGzTaa8uwwvy5c7lvLIjbXMjZ7KXfNlfdutLhvo8V9623WWrewUqGIiIicyt3ZGuZc3B5zaWfEpe0Jl7ZHXNoZc3F7PLvLujMurvvZXhbNwurdaxl39zPOdVM2OglnWyHnoxFn2aJfbRGNrsDe0zB4ui5n9Wfqsrr+9QniJpzuB9XNw8G1u1nffW01z8PG7Zv6ntbKK4b58FAwPTpLbb49zOsbCqPi4ObC/g2Hm2UYWZSRhRlplJKFGVmUkYbpbP+sPVdvRa2DPnPH4zAmCZJZmYTJrB0Hcd0OE+Ig1sJQt8nzDaI/CPyOu7++af8WgLv//lyfjzd9/t3MIuAScN5PefFVD6JPbj1J4QVlVVJ6SVEVVF7N6qWXs2P79cKbPk19dnyu33U/Wx3ePwuQZX4oTB5qHzk2C5pVTnHcG9RzkAUJaRiTBjGZxaRBRGohqYVkFpISkFpAhtEhoEVA26ENtCtoe0W7ctpVQbvIaRdT2vmEdj4my4cE+RjyYb3dQoAFKMMWRdSmiNrkYYs8aDMJWkzCNhNrMbaMkbUYkTEiY0DKgBZ7VcquZ+xWKbtlwnaVslMm7JQxUw/JS6esKorKKSunqJyirLjBugnPWRwarTiknUT0WzH9dsxGO2a9lbDeqcuNdsx6O+ZsN+WuXsbmWkoWazECERGRVTScFlzaru+iXt4PqU17f/8zu5MTP1v0soiznYRuFtFOIjpJ/TmhnYR0YmMjHHHWt+hXz9LOr9HJr9KeXqOdX6U1vUY2uUprepV0cpXQj/9MWFrENF4jj3pM4jWm8RrTqDcrJ9Ea07iHR636In2cQdzG4hYWd7AkI0w6WNomSFrEUUIcGXEYkIQBcRgQh0Yc1e0wgLyaXBdOB/lgNitukA8Yl2NGxWh2w2RcjA9ufDTtSTmp+zQ3S/brt3Iz5CSRRXVgbcLqfkCdD7FREBFaSBRE9WZ1GQbhrD5/LAzCQ/2u+/mj+ywisIAwCAktrOsWXtee73Ohd2Glp0g/32dE7wXmvq2ep4DvP6mPuxdmtg2cBa7c+nBXw8/8wxso7PYt5BQ6BEDoRsBBPQJiN0KHmLqMgMiN2CEBIocII3JIHGJ3YncSd2KsqVckTRl7ReIVqVekVUnbS7KmbDX7M3dSd7LKidlfH/XW5MTkFlM0ZW4xORFjy9gm5TIJQ99gxN0MPWHkCQNPGFR1uVfV+0ak9dbUh2QMPGPY1CtOv2IVBTZb+r0uw7l6QJrWS8N3ooCN5k0yCo0wCIgCIwqNKDi5HQZGHB5u132Cuh4YYfMzUVC/9v5y9e25PyhxqCtvIiIiLybtJOIV57u84nz3xD5l5WwNp1wbTLk6qMv57epgymBSMJgUXB1M+ca1IaNpyWBaMpwW5KUDvWZ74ITf4vQZcM62OW/bnGGHvg1Yo16kqZ8PmvaANbtInyc4Z0P6DIjs1kLd1EPGJORElIQUBAw8pCCctUtCKgspiaisrrcsoG1gZs0srmYzrq/P7uDu10OMNm4dHCjMmJrPbTA1JwcKc3JzSuqywCnmyny/bVDM+owpbHSo79icEqiassQp9+vN61dNvZj1q4/7bZyk9uev/QA/8KofvX2/4DZa6IokZvY24G0A999//yJ/9S37tSsQUhF6/WCsuRECwaxkVgbN/ggIsHqf2axPHTrrtlEHT2vC1EHUtaZt881Ze7b/SLuy+h93SUhpIdWhelC3bb+sj20RUob18cJiCovJacr9enAQKOsymYXL/X4Fh+s298yhNeMzqxceSMKDULYf/OLmqlkU1FfNNgMjCpuwGAUkRwNkHJKE9feMJWFAduT7xpJZqFTAExERkdUUBsbZbsrZbsqDz+Hnp0VFXlZU7vV6ie5N3Zs6VJVj+yGv+Tmz+vOZzX/UtIPPbG6w7Q75gGCyQzkdUk4GVJMRVT6kmgyppkM8H+HTYb1ORl5vlo/wqoAyx6uiqRcEVUlSFVDlUBVYVYIXBFWBe4m7gzvuFc7Bp1z2J1W6Y+7NMW+O1X2Mgz7MHT/4xLycbwY5TgWU0IRdKK3e6sBal4XV/SqjCbxNeUy9al6jAtbt9nxP7iLcTBD9FnBhrn1fs++4Pk81U3P71IsWHeLuHwQ+CPXU3Ocy4EX5+Xd9ftlDEBERERE5JGkuvt8+GfXERpHb62bO4seAB83s5WaWAG8BHjnS5xHgF5r6zwKfPO35UBEREREREXnpuuEd0eaZz3cAH6eeafphd/+imb0P+Iy7PwJ8CPgrM3sCuEYdVkVERERERESuc1PPiLr7o8CjR/a9Z64+Bt70wg5NREREREREXoy0qouIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJSCqIiIiIiIiCyUgqiIiIiIiIgslIKoiIiIiIiILJS5+3J+sdkzwNeX8stv3jngyrIHIStJ54acROeGnEbnh5xE54acROeGnOROODcecPfzxx1YWhC9E5jZZ9z9+5Y9Dlk9OjfkJDo35DQ6P+QkOjfkJDo35CR3+rmhqbkiIiIiIiKyUAqiIiIiIiIislAKoqf74LIHICtL54acROeGnEbnh5xE54acROeGnOSOPjf0jKiIiIiIiIgslO6IioiIiIiIyEIpiN6Amf2emX3ezD5nZp8ws5cte0yyGszsj8zsK8358fdmtr7sMclqMLM3mdkXzawyszt2NTt54ZjZQ2b2uJk9YWa/uezxyOowsw+b2dNm9oVlj0VWi5ldMLNPmdmXmr8p71z2mGQ1mFlmZv9pZv/dnBu/u+wxPReamnsDZrbm7jtN/VeAV7v725c8LFkBZvYTwCfdvTCzPwRw999Y8rBkBZjZq4AK+DPgXe7+mSUPSZbIzELgq8DrgKeAx4C3uvuXljowWQlm9sPAHvARd/+uZY9HVoeZ3QPc4+6fNbMe8F/AG/TeIWZmQMfd98wsBv4NeKe7f3rJQ7sluiN6A/shtNEBlNwFAHf/hLsXTfPTwH3LHI+sDnf/srs/vuxxyMp4LfCEuz/p7lPgb4CHlzwmWRHu/i/AtWWPQ1aPu19098829V3gy8C9yx2VrAKv7TXNuNnuuIyiIHoTzOz9ZvZN4OeA9yx7PLKSfgn4x2UPQkRW0r3AN+faT6EPkyJyC8zsO4DvAf5juSORVWFmoZl9Dnga+Cd3v+PODQVRwMz+2cy+cMz2MIC7v9vdLwAfBd6x3NHKIt3o3Gj6vBsoqM8PeYm4mXNDRETk+TKzLvAx4FePzNSTlzB3L939u6ln5L3WzO64qf3RsgewCtz9x2+y60eBR4H33sbhyAq50blhZr8I/BTwY64Hrl9SbuF9Q+RbwIW59n3NPhGRUzXP/30M+Ki7/92yxyOrx923zOxTwEPAHbXome6I3oCZPTjXfBj4yrLGIqvFzB4Cfh34aXcfLns8IrKyHgMeNLOXm1kCvAV4ZMljEpEV1yxI8yHgy+7+x8sej6wOMzu//20NZtaiXgzvjssoWjX3BszsY8B3Uq+A+XXg7e6uK9mCmT0BpMDVZtentaKyAJjZG4E/Bc4DW8Dn3P31yx2VLJOZ/STwJ0AIfNjd37/kIcmKMLO/Bn4EOAdcBt7r7h9a6qBkJZjZDwH/CvwP9edQgN9290eXNypZBWb2GuAvqf+mBMDfuvv7ljuqW6cgKiIiIiIiIgulqbkiIiIiIiKyUAqiIiIiIiIislAKoiIiIiIiIrJQCqIiIiIiIiKyUAqiIiIiIiIislAKoiIiIiIiIrJQCqIiIiIiIiKyUAqiIiIiIiIislD/DyJCLRSUSVviAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "smc_posteriorpdf = gaussian_kde(smc_samples.samples_array[:,0])\n",
    "adaptive_smc_posteriorpdf = gaussian_kde(adaptive_smc_samples.samples_array[:,0])\n",
    "\n",
    "reference_posteriorpdf = gaussian_mixture_density\n",
    "\n",
    "xs = np.linspace(-3,3,200)\n",
    "smc_posteriorpdf.covariance_factor = lambda : .25\n",
    "smc_posteriorpdf._compute_covariance()\n",
    "adaptive_smc_posteriorpdf.covariance_factor = lambda : .25\n",
    "adaptive_smc_posteriorpdf._compute_covariance()\n",
    "plt.figure(figsize=(16,10))\n",
    "plt.plot(xs,smc_posteriorpdf(xs))\n",
    "plt.plot(xs,adaptive_smc_posteriorpdf(xs))\n",
    "plt.plot(xs,reference_posteriorpdf(xs))\n",
    "plt.legend(('abc-smc', 'adaptive abc-smc', 'reference'));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Simola, U., Cisewski-Kehe, J., Gutmann, M.U. and Corander, J. [*Adaptive Approximate Bayesian Computation Tolerance Selection*](https://projecteuclid.org/journals/bayesian-analysis/advance-publication/Adaptive-Approximate-Bayesian-Computation-Tolerance-Selection/10.1214/20-BA1211.full), Bayesian Analysis 1(1):1-27, 2021\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.9 64-bit ('elfi36': conda)",
   "language": "python",
   "name": "python36964bitelfi36conda8bec451552304612bd355a97e3742bbb"
  },
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}