{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SMC samplers\n",
    "\n",
    "This tutorial gives a basic introduction to SMC samplers, and explains how to run the SMC samplers already implemented in ``particles``. For a more advanced tutorial on how to design new SMC samplers, see the next tutorial. For more background on SMC samplers, check either Chapter 17 of the book or [Dau & Chopin (2022)](https://doi.org/10.1111/rssb.12475) for waste-free SMC. Arxiv version is [here](http://arxiv.org/abs/2011.02328). \n",
    "\n",
    "## SMC samplers: what for? \n",
    "\n",
    "A SMC sampler is a SMC algorithm that samples from a sequence of probability distributions $\\pi_t$, $t=0,\\ldots,T$ (and compute their normalising constants). Sometimes one is genuinely interested in each $\\pi_t$; more often one is interested only in the final distribution $\\pi_T$. In the latter case, the sequence is purely instrumental.\n",
    "\n",
    "Examples of SMC sequences are: \n",
    "\n",
    "1. $\\pi_t(\\theta) = p(\\theta|y_{0:t})$, the Bayesian posterior distribution of parameter $\\theta$ given data $y_{0:t}$, for a certain model. \n",
    "\n",
    "2. A tempering sequence, $\\pi_t(\\theta) \\propto \\nu(\\theta) L(\\theta)^{\\gamma_t}$ ,where the $\\gamma_t$'s form an increasing sequence of exponents: $0=\\gamma_0 < \\ldots < \\gamma_T=1$. You can think of $\\nu$ being the prior, $L$ the likelihood function, and $\\pi_T$ the posterior. However, more generally, tempering is a way to interpolate between any two distributions, $\\nu$ and $\\pi$, with $\\pi(\\theta) \\propto \\nu(\\theta) L(\\theta)$. \n",
    "\n",
    "We discuss first how to specify a sequence of the first type.\n",
    "\n",
    "## Defining a Bayesian model\n",
    "\n",
    "To define a particular Bayesian model, you must subclass `StaticModel`, and define method `logpyt`, which evaluates the log-likelihood of datapoint $Y_t$ given parameter $\\theta$ and past datapoints $Y_{0:t-1}$. Here is a simple example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sb\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "import particles\n",
    "from particles import smc_samplers as ssp\n",
    "from particles import distributions as dists\n",
    "\n",
    "class ToyModel(ssp.StaticModel):\n",
    "    def logpyt(self, theta, t):  # density of Y_t given theta and Y_{0:t-1}\n",
    "        return stats.norm.logpdf(self.data[t], loc=theta['mu'],\n",
    "                                 scale = theta['sigma'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In words, we are considering a model where the observations are $Y_t\\sim N(\\mu, \\sigma^2)$ (independently). The parameter is $\\theta=(\\mu, \\sigma)$. Note the fields notation; more about this later. \n",
    "\n",
    "Class `ToyModel` implicitly defines the  likelihood of the considered model for any sample size (since the likelihood at time $t$ is $p^\\theta(y_{0:t})=\\prod_{s=0}^t p^\\theta(y_s|y_{0:s-1})$, and method `logpyt` defines each factor in this product; note that $y_s$ does not depend on the past values in our particular example). We now define the data and the prior: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 30\n",
    "my_data = stats.norm.rvs(loc=3.14, size=T)  # simulated data\n",
    "my_prior = dists.StructDist({'mu': dists.Normal(scale=10.),\n",
    "                             'sigma': dists.Gamma()})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more details about to define prior distributions, see the documentation of module `distributions`, or the previous [tutorial on Bayesian estimation of state-space models](Bayes_estimation_ssm.ipynb). Now that we have everything, let's specify our static model: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_static_model = ToyModel(data=my_data, prior=my_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, object `my_static_model` entirely defines the posterior. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.37897769e+02, -8.61877866e+01, -4.94969347e+02, -2.32511323e+05,\n",
       "       -2.63076061e+02])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thetas = my_prior.rvs(size=5) \n",
    "my_static_model.logpost(thetas, t=2)  \n",
    "# if t is omitted, gives the full posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input of `logpost` and output of `myprior.rvs()` are [structured arrays](https://docs.scipy.org/doc/numpy/user/basics.rec.html), that is, arrays with fields:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.241375242422195"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thetas['mu'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Typically, you won't need to call `logpost` yourself, this will be done by the SMC sampler for you. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IBIS\n",
    "\n",
    "IBIS (iterated batch importance sampling) is the standard name for a SMC sampler that tracks a sequence of partial posterior distributions; i.e. $\\pi_t$ is $p(\\theta|y_{0:t})$, for $t=0,1,\\ldots$. \n",
    "\n",
    "Module `smc_samplers` defines `IBIS` as a subclass of `FeynmanKac`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=0, ESS=30.38\n",
      "t=1, Metropolis acc. rate (over 49 steps): 0.179, ESS=320.47\n",
      "t=2, Metropolis acc. rate (over 49 steps): 0.265, ESS=739.71\n",
      "t=3, ESS=495.11\n",
      "t=4, Metropolis acc. rate (over 49 steps): 0.239, ESS=699.01\n",
      "t=5, ESS=362.60\n",
      "t=6, Metropolis acc. rate (over 49 steps): 0.357, ESS=855.81\n",
      "t=7, ESS=563.24\n",
      "t=8, ESS=429.52\n",
      "t=9, Metropolis acc. rate (over 49 steps): 0.348, ESS=944.20\n",
      "t=10, ESS=821.65\n",
      "t=11, ESS=672.69\n",
      "t=12, ESS=600.56\n",
      "t=13, ESS=512.21\n",
      "t=14, ESS=812.56\n",
      "t=15, ESS=738.51\n",
      "t=16, ESS=675.90\n",
      "t=17, ESS=618.24\n",
      "t=18, ESS=558.24\n",
      "t=19, ESS=482.73\n",
      "t=20, Metropolis acc. rate (over 49 steps): 0.349, ESS=972.66\n",
      "t=21, ESS=925.96\n",
      "t=22, ESS=889.08\n",
      "t=23, ESS=970.02\n",
      "t=24, ESS=950.08\n",
      "t=25, ESS=732.06\n",
      "t=26, ESS=714.42\n",
      "t=27, ESS=453.15\n",
      "t=28, Metropolis acc. rate (over 49 steps): 0.316, ESS=978.24\n",
      "t=29, ESS=951.31\n"
     ]
    }
   ],
   "source": [
    "my_ibis = ssp.IBIS(my_static_model, len_chain=50)\n",
    "my_alg = particles.SMC(fk=my_ibis, N=20, \n",
    "                       store_history=True, verbose=True)\n",
    "my_alg.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: we use option `verbose=True` in `SMC` in order to print some information on the intermediate distributions. \n",
    "\n",
    "**Note**: Since we set `store_history` to `True`, the particles and their weights have been saved at every time (in attribute `hist`, see previous tutorials on smoothing). Let's plot the posterior distributions of $\\mu$ and $\\sigma$ at various times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG0CAYAAADgoSfXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA58UlEQVR4nO3df1QU973/8Re4SwT5GYELCEYwgj8iFduaxLT1V6L9gT/IiTk2TVvi9SatxCQ3MWm85oeYHhPaRHMa7c9g1FjitRgatb2RGm2j10Zv0kaaEE0i0EYLylZWICLuwn7/8MuSlUV2F5bZhefjHM9hZj+z89436/Ca2dmZEIfD4RAAAIBBQo0uAAAADG6EEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKJPRBQAY+MrKynTkyBGdOnVKYWFhyszM1J133qmUlJRul3n//fdVWFjYZf66des0YsQIf5YLoJ8RRgD4XWVlpebMmaPRo0erra1N27Zt0w9/+EOtXbtWQ4cOveKyzz//vCIiIpzT0dHR/i4XQD8LqjDS0NAgu93u0diEhATV19f7uaKBh775ZqD3zWQyKS4uzuflV65c6TK9dOlSLVmyRFVVVRo/fvwVl42JidGwYcN8XndP242B/rvzBr3oRC869aYXnm47giqM2O122Wy2HseFhIQ4x3NTYs/RN9/QN++dP39ekhQZGdnj2EceeUQ2m02pqam69dZbdd1117kdZ7PZXLYPISEhCg8Pl91u7zaMdPzu2traBv3vjl50ohed+qsXQRVGAAQ/h8OhzZs3a+zYsRo5cmS34+Li4nT33XcrIyNDdrtdb775pp566ik9+eSTbo+mlJWVqbS01Dmdnp6uoqIiJSQk9FhTUlKSby9mAKIXnehFJ3/3gjACoF8VFxfrH//4h1avXn3FcSkpKS4nuGZmZspisWjXrl1uw0heXp5yc3Od0x17dPX19Vc8MpKUlKS6ujr2gOmFE73o1NtemEwmj3YICCMA+s3GjRv1zjvvqLCwUMOHD/d6+czMTB04cMDtY2azWWaz2e1jPW1EHQ7HoP+j04FedKIXnfzdC64zAsDvHA6HiouLdfjwYT3xxBNKTEz06Xmqq6sVGxvbt8UBMBxHRgD4XXFxsQ4ePKhHHnlE4eHhslqtkqSIiAiFhYVJkkpKSnT27Fnde++9kqTf/e53SkhIUFpamux2uw4cOKDDhw/roYceMuplYACx2+3OE6ndaWlp0cWLF/uxosDVUy8iIiJkMvUuThBGAPhdeXm5JGnVqlUu85cuXarp06dLuvQVXIvF4nzMbrfr5Zdf1tmzZxUWFqa0tDQ9+uijmjx5cn+VjQHKbrfr008/VVRUlEJD3X9AYDabPfr25mBwpV60t7erqalJw4YN61UgIYwA8Lvt27f3OKagoMBlev78+Zo/f76/SsIgdv78+SsGEXguNDRUUVFRam5u7tUFCflNAAAGHYJI3+mLXvLbAAAAhiKMAAAAQxFGAACAoTiBFQAw6LXvLHGZtoUOUXt7m9/WFzrvDq+Xue222zR+/Pger17cndOnT2v16tWqqKhQdXW1Fi9e7PNz9TWOjAAAMAhcvHhRw4cP13333dfj3bL7G2EEAIAA98ADD+jPf/6ziouLNWLECI0YMUKffPKJV8+Rlpam1atXa+HChb36Gq4/8DENAAABbvXq1aqqqtLYsWO1fPlySdLw4cM1ZsyYKy53/fXXa+vWrf1RYq8QRjBgvFJRL0n6ZnbPd4gEAlnHe7kD72lER0crLCxMQ4cOdbm3U8fVjbszdOhQf5fWJwgjAAAEqfT0dKNL6BOEEQAAgtSg/ZimsrJSO3fuVHV1tRoaGrR8+XJNmTLlisvYbDaVlpbqwIEDslqtGj58uPLy8jRz5kyfCwcAYDAxm81qb293mTdoP6ZpbW3VqFGjNGPGDD333HMeLbNu3TqdO3dO3/ve95SUlKTGxka1tfnv+9sAAAw0aWlp+utf/6pPPvlEw4YNU2xsrNcf07z33nuSpE8//VRnz57Ve++9p7CwMGVmZvqjZI95HUZycnKUk5Pj8fh3331XlZWVWr9+vSIjIyXJ5eQboLcuP9kPAAaie+65Rw888ICmT5+uCxcu6K233lJaWppXzzFnzhznzxUVFSorK1NqaqoOHz7c1+V6xe/njLz99tsaPXq0XnvtNb355psaOnSoPv/5z2vRokUKCwtzu4zNZpPNZnNOh4SEKDw83PlzTzrGeDIWnYK3b6719nf9wds3AB0uvyKq2Wx2+TsUCEaPHq1du3b16jlOnTrVR9X0Lb+HkdOnT+vYsWMym816+OGH1djYqOLiYjU3N2vp0qVulykrK1NpaalzOj09XUVFRUpI8O7rbUlJSb2qfbAKtr5FVbW4TCcnJxtSR7D1DQAChd/DiMPhkCTdd999ioiIkHTpyMfatWu1ZMkSt0dH8vLylJub65zu2OOsr6+X3W7vcZ0hISFKSkpSXV2dc/3oWbD2ramp2WW6tra2X9cfrH3zhslk8npnAAA85fcwEhsbq6uvvtoZRCRpxIgRcjgc+te//uV2L9ZsNstsNrt9Pm829g6HY8D+cfCn4Ouba61G1R58fQOAwOD3e9OMHTtWDQ0NunDhgnNebW2tQkJCNHz4cH+vHgAABDivw8iFCxdUU1OjmpoaSdKZM2dUU1Mji8UiSSopKdH69eud47/0pS8pKipKP/3pT3Xy5ElVVlZq69atmjFjRrcnsAIAgMHD649pTpw4ocLCQuf0li1bJEnTpk1TQUGBGhoanMFEunTBlccee0wbN27Uo48+qqioKN14441atGhRH5QPAACCnddhZMKECdq+fXu3jxcUFHSZN2LECD3++OPergoAAAwCfj9nBAAA4EoIIwAAwFCEEQAAYCi/X2cEAIBAd/k9rkJDh6i93X83dP1mtvcXEbzttts0fvx4rV692qd1/v73v9eWLVv0/vvv6+LFi8rMzNRDDz2k6dOnO8fYbDatX79ev/nNb1RXV6eMjAw9+eST+vKXv+zTOj3FkREAAAaBt956S1/5ylf08ssv63/+5380depU5efnO+/kK0k/+tGPtHXrVj311FPav3+/vv3tb3cZ4w+EEQAAAtwDDzygP//5zyouLtaIESM0YsQIffLJJ149x+rVq7V06VJNmjRJGRkZWrFihdLT0/WHP/zBOWbHjh1atmyZZs2apWuuuUbf/e53NX36dP3iF7/o65fkgo9pAAAIcKtXr1ZVVZXGjh2r5cuXS5KGDx+uMWPGXHG566+/Xlu3bnX7WHt7u5qbmxUbG+uc19raqquuusplXHh4uI4cOdK7F9ADwggAAAEuOjpaYWFhGjp0qBITE53zy8vLr7jc0KFDu33sF7/4hc6fP6+5c+c6502fPl2//OUvdf3112vUqFE6ePCgXn/9dbW1+e/8GYkwAgBA0EpPT/dpud/+9rd67rnntHHjRsXHxzvnr169Wg8//LCmTZumkJAQXXPNNVq0aJFeeeWVvirZLcIIAABBypePaV577TU99NBD+sUvfqGvfOUrLo8NHz5cGzdu1IULF9TQ0KCkpCQ988wzGjlyZJ/X/lmEEQAAgoDZbFZ7e7vLPG8/pvntb3+rhx56SBs2bNDNN998xeWSk5Nls9m0e/du5ebm+l64BwgjAAAEgbS0NP31r3/VJ598omHDhik2Ntarj2l++9vf6v7771dhYaEmT56sM2fOSLoUPKKjoyVJf/nLX1RXV6cJEyaorq5Ozz33nNrb27V06VK/vKYOfLUXAIAgcM899yg0NFTTp0/XxIkTderUKa+W37p1q+x2u1auXKmcnBznvyeeeMI5prW1VT/60Y80Y8YM/fu//7uSkpK0e/duxcTE9PXLccGREQDAoHf5FVHNZrNsNptB1bg3evRo7dq1y+flS0tLexxz44036o9//KPLvP7oBUdGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACABh0Lr94GHzXF70kjAAABpWIiAg1NTURSPpAe3u7mpqaFBER0avn4TojABDgXqmod/58+fUw4D2TyaRhw4apubm52zFhYWG6ePFiP1YVuHrqxbBhw2Qy9S5OEEYAAIOOyWRyXgL9ciEhIUpOTlZtba0cDkc/VxZY+qsXfEyDAeeVinqXPUkAQGAjjAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhvI6jFRWVuqZZ57RPffco9tvv11HjhzxeNljx45p0aJFevjhh71dLQAAGKC8DiOtra0aNWqUFi9e7NVy58+f14YNGzRx4kRvVwkAAAYwr+9Nk5OTo5ycHK9X9Mtf/lI33XSTQkND9X//939eLw8AAAamfrlR3v79+3X69GktW7ZMO3bs6HG8zWaTzWZzToeEhCg8PNz5c086xngyFp2Ct2/u6+2v1xG8fQOAwOD3MFJbW6uSkhIVFhZqyJAhHi1TVlam0tJS53R6erqKioqUkODdrbOTkpK8Go9Lgq1vUVUtbucnJyf3ax3B1jcACBR+DSPt7e36yU9+ooULFyolJcXj5fLy8pSbm+uc7tjjrK+vl91u73H5kJAQJSUlqa6ubtDf/tkbwdq3pqZmt/Nra2v7Zf3B2jdvmEwmr3cGAMBTfg0jLS0tOnHihKqrq7Vx40ZJksPhkMPh0KJFi/TYY4/puuuu67Kc2WyW2Wx2+5zebOw71gXvBF/f3Nfa368h+PoGAIHBr2EkPDxczz77rMu88vJyvffee3rwwQeVmJjoz9UDAIAg4HUYuXDhgurq6pzTZ86cUU1NjSIjIxUfH6+SkhKdPXtW9957r0JDQzVy5EiX5aOjo2U2m7vMBzBwlZWV6ciRIzp16pTCwsKUmZmpO++8s8ePbysrK7V582adPHlScXFxmjdvnmbPnt1PVQPoL16HkRMnTqiwsNA5vWXLFknStGnTVFBQoIaGBlkslr6rEEDQq6ys1Jw5czR69Gi1tbVp27Zt+uEPf6i1a9dq6NChbpc5c+aMnn76ac2aNUvLli3T8ePH9eKLLyo6Olo33HBDP78CAP7kdRiZMGGCtm/f3u3jBQUFV1z+9ttv1+233+7tagEEsZUrV7pML126VEuWLFFVVZXGjx/vdpny8nLFx8crPz9fkpSamqoTJ05o165dhBFggOmX64wAwGedP39ekhQZGdntmI8++kjZ2dku8yZNmqT9+/fLbrfLZHLdfPlyfaLAvUZM9/X4q9bA7UX/oxed+qsXhBEA/crhcGjz5s0aO3bsFc8ds1qtiomJcZkXExOjtrY2NTU1KS4uzuWx3lyfKNCuEdPdtXMk/18/J9B6YSR60cnfvSCMAOhXxcXF+sc//qHVq1f3OPbyvbGOr06720vz5fpEgXqNmO6unSP57/o5gdoLI9CLTr3thafXKCKMAOg3Gzdu1DvvvKPCwkINHz78imNjY2NltVpd5jU2NmrIkCFuP97pzfWJAu8aMd3X4u86A68XxqEXnfzdC6/v2gsA3nI4HCouLtbhw4f1xBNPeHSNoTFjxqiiosJl3tGjR5WRkdHlfBEAwY0wAsDviouLdeDAAd1///0KDw+X1WqV1WrVxYsXnWNKSkq0fv165/Ts2bNlsVic1xnZt2+f9u3bp7lz5xrxEgD4EbsXAPyuvLxckrRq1SqX+UuXLtX06dMlqcs1ihITE7VixQpt3rxZe/bsUVxcnO666y6+1gsMQIQRAH53pWsTdXB3jaLx48erqKjIHyUBCCCEEQStVyrqjS4BANAHOGcEAAAYijACAAAMRRgBAACGIowAAABDEUYAAICh+DYNABiMb4ZhsOPICAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADCUydsFKisrtXPnTlVXV6uhoUHLly/XlClTuh1/+PBhlZeXq6amRna7XampqVq4cKEmTZrUm7oBAMAA4fWRkdbWVo0aNUqLFy/2aPwHH3yg7OxsrVixQs8884wmTJigoqIiVVdXe10sAAAYeLw+MpKTk6OcnByPx+fn57tM33HHHXr77bf1zjvvKD093dvVAwCAAcbrMNJb7e3tamlpUWRkZLdjbDabbDabczokJETh4eHOn3vSMcaTsegUfH27cp399TqCr28AEFj6PYzs3r1bra2tuvHGG7sdU1ZWptLSUud0enq6ioqKlJCQ4NW6kpKSfK5zMAuWvkVVtVzx8eTk5H6q5JJg6RsABJp+DSMHDx7Ub37zGz388MOKiYnpdlxeXp5yc3Od0x17nPX19bLb7T2uJyQkRElJSaqrq5PD4eh94YNEsPWtqan5io/X1tb2Sx3B1jdfmEwmr3cGAMBT/RZGDh06pJ///Od68MEHlZ2dfcWxZrNZZrPZ7WPebOwdDseA/ePgT8HTtyvX2N+vIXj6BgCBpV/CyMGDB/Wzn/1M999/vyZPntwfqwQAAEHC6zBy4cIF1dXVOafPnDmjmpoaRUZGKj4+XiUlJTp79qzuvfdeSZeCyIYNG5Sfn6/MzExZrVZJUlhYmCIiIvrmVQAAgKDldRg5ceKECgsLndNbtmyRJE2bNk0FBQVqaGiQxWJxPr537161tbWpuLhYxcXFzvkd4wEAwODmdRiZMGGCtm/f3u3jlweMVatWeV0UAAAYPLg3DQAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhjIZXQCAga+yslI7d+5UdXW1GhoatHz5ck2ZMqXb8e+//74KCwu7zF+3bp1GjBjhz1IBGIAwAsDvWltbNWrUKM2YMUPPPfecx8s9//zzioiIcE5HR0f7ozwABiOMAPC7nJwc5eTkeL1cTEyMhg0b5oeKAAQSwgiAgPXII4/IZrMpNTVVt956q6677rpux9psNtlsNud0SEiIwsPDnT+70zG/u8f7j+fr91etgdML49GLTv3VC8IIgIATFxenu+++WxkZGbLb7XrzzTf11FNP6cknn9T48ePdLlNWVqbS0lLndHp6uoqKipSQkNDj+pKSkvqsdl9EVbV4PDY5OdmPlRjfi0BCLzr5uxdehxFvT0TrWGbz5s06efKk4uLiNG/ePM2ePdvnogEMbCkpKUpJSXFOZ2ZmymKxaNeuXd2Gkby8POXm5jqnO/bk6uvrZbfb3S4TEhKipKQk1dXVyeFw9OEr8E5TU7PHY2tra/1SQ6D0IhDQi0697YXJZPJoh8DrMOLtiWhnzpzR008/rVmzZmnZsmU6fvy4XnzxRUVHR+uGG27wdvUABqnMzEwdOHCg28fNZrPMZrPbx3raiDocDoP/6Hi+bn/XaXwvAge96OTvXngdRrw9Ea28vFzx8fHKz8+XJKWmpurEiRPatWsXYQSAx6qrqxUbG2t0GQD8wO/njHz00UfKzs52mTdp0iTt379fdrtdJlPXEnw5Ee2zOPnIN8HXtyvX2V+vI/j61v8uXLiguro65/SZM2dUU1OjyMhIxcfHq6SkRGfPntW9994rSfrd736nhIQEpaWlyW6368CBAzp8+LAeeugho14CAD/yexixWq2KiYlxmRcTE6O2tjY1NTUpLi6uyzK9ORHtszj5yDfB0reeTvrz94l+lwuWvhnhxIkTLhcx27JliyRp2rRpKigoUENDgywWi/Nxu92ul19+WWfPnlVYWJjS0tL06KOPavLkyf1eOwD/65dv01y+x9jxuVN3e5K+nIh2+fo4+ch7wda3nk7689eJfpcLtr75wtOT0LozYcIEbd++vdvHCwoKXKbnz5+v+fPn+7w+AMHF72EkNjZWVqvVZV5jY6OGDBmiyMhIt8v05kS0y8cO1D8O/hQ8fev5pMT+FDx9A4DA4vcb5Y0ZM0YVFRUu844ePaqMjAy354sAAIDBxeswcuHCBdXU1KimpkZS54loHZ/3lpSUaP369c7xs2fPlsVicV5nZN++fdq3b5/mzp3bN68AAAAENa8PTXh7IlpiYqJWrFihzZs3a8+ePYqLi9Ndd93F13oBAIAkH8KItyeiSdL48eNVVFTk7aoAAMAg4PdzRgAAAK6EMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxlMroAwFuvVNQbXQIAoA9xZAQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKF8+jbNnj17tHPnTlmtVqWmpio/P1/jxo3rdvyBAwe0c+dO1dbWKiIiQpMmTdK3v/1tRUVF+Vw4AAAYGLw+MnLo0CFt2rRJt956q4qKijRu3DitWbNGFovF7fhjx45p/fr1mjFjhtauXasHH3xQJ06c0M9//vNeFw8Ag80rFfUu/4CBwOswsnv3bs2cOVOzZs1yHhWJj49XeXm52/EffvihEhMT9fWvf12JiYkaO3asbr75ZlVVVfW6eAAAEPy8+pjGbrerqqpKCxYscJmfnZ2t48ePu10mKytL27Zt01/+8hfl5OTo3Llzeuutt5STk9Ptemw2m2w2m3M6JCRE4eHhzp970jHGk7HoFDx986y+/nodwdM3AAhMXoWRxsZGtbe3KyYmxmV+TEyMrFar22WysrJ033336fnnn5fNZlNbW5u+8IUvaPHixd2up6ysTKWlpc7p9PR0FRUVKSEhwZtylZSU5NV4XBLofYuqavFoXHJysp8rcRXofQOAQOXTCazu9gC72ys8efKkXnrpJd1222363Oc+p4aGBm3dulW/+tWv9P3vf9/tMnl5ecrNze3y3PX19bLb7R7Vl5SUpLq6OjkcDk9eEhQ8fWtqavZoXG1trZ8ruSRY+tYbJpPJ650BAPCUV2EkOjpaoaGhXY6CnDt3rsvRkg5lZWXKysrSvHnzJEnXXHONhg4dqieeeEKLFi1SXFxcl2XMZrPMZrPb5/NmY+9wOAbsHwd/Cvy+eVZbf7+GwO8bAAQmr05gNZlMysjIUEVFhcv8iooKZWVluV2mtbW1y1GT0NBLq2XDDQAAvP42TW5urt544w3t27dPJ0+e1KZNm2SxWHTLLbdIkkpKSrR+/Xrn+C984Qs6cuSIysvLdfr0aR07dkwvvfSSrr32Wl199dV990oAAEBQ8vqckalTp6qpqUk7duxQQ0OD0tLStGLFCufnyQ0NDS7XHJk+fbpaWlr0+uuva8uWLRo2bJgmTJigO++8s+9eBQAACFo+ncA6Z84czZkzx+1jBQUFXeZ97Wtf09e+9jVfVgUAAAY47k0DAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCiT0QUAGPgqKyu1c+dOVVdXq6GhQcuXL9eUKVN6XGbz5s06efKk4uLiNG/ePM2ePbufKgbQnzgyAsDvWltbNWrUKC1evNij8WfOnNHTTz+tcePGqaioSHl5eXrppZf01ltv+blSAEbgyAgAv8vJyVFOTo7H48vLyxUfH6/8/HxJUmpqqk6cOKFdu3bphhtu8FOVAIxCGAEQcD766CNlZ2e7zJs0aZL2798vu90uk6nrpstms8lmszmnQ0JCFB4e7vzZnY753T3ef3xff1/VHji9MB696NRfvSCMAAg4VqtVMTExLvNiYmLU1tampqYmxcXFdVmmrKxMpaWlzun09HQVFRUpISGhx/UlJSX1vuheiKpq8XnZ5OTkPqzE+F4EEnrRyd+9IIwACEiX74k5HA638zvk5eUpNze3y/L19fWy2+3driMpKUl1dXXO5zdCU1Ozz8vW1tb2SQ2B0otAQC869bYXJpPJox0CwgiAgBMbGyur1eoyr7GxUUOGDFFkZKTbZcxms8xms9vHetqIOhwOg//o+L7uvq7b+F4EDnrRyd+94Ns0AALOmDFjVFFR4TLv6NGjysjIcHu+CIDgRhgB4HcXLlxQTU2NampqJF366m5NTY0sFoskqaSkROvXr3eOnz17tiwWi/M6I/v27dO+ffs0d+5cI8oH4GfsYgDwuxMnTqiwsNA5vWXLFknStGnTVFBQoIaGBmcwkaTExEStWLFCmzdv1p49exQXF6e77rqLr/UCAxRhBIDfTZgwQdu3b+/28YKCgi7zxo8fr6KiIn+WBSBA8DENAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQPl1nZM+ePdq5c6esVqtSU1OVn5+vcePGdTveZrOptLRUBw4ckNVq1fDhw5WXl6eZM2f6XDgAABgYvA4jhw4d0qZNm7RkyRJlZWVp7969WrNmjdatW6f4+Hi3y6xbt07nzp3T9773PSUlJamxsVFtbW29Lh4AAAQ/r8PI7t27NXPmTM2aNUuSlJ+fr6NHj6q8vFx33HFHl/HvvvuuKisrtX79eufdNhMTE3tZNgAAGCi8CiN2u11VVVVasGCBy/zs7GwdP37c7TJvv/22Ro8erddee01vvvmmhg4dqs9//vNatGiRwsLC3C5js9lks9mc0yEhIQoPD3f+3JOOMZ6MRafg6Ztn9fXX6wievgFAYPIqjDQ2Nqq9vV0xMTEu82NiYmS1Wt0uc/r0aR07dkxms1kPP/ywGhsbVVxcrObmZi1dutTtMmVlZSotLXVOp6enq6ioSAkJCd6Uq6SkJK/G45JA71tUVYtH45KTk/1ciatA7xsABCqfTmB1twfY3V6hw+GQJN13332KiIiQdOnIx9q1a7VkyRK3R0fy8vKUm5vb5bnr6+tlt9s9qi8pKUl1dXXO9aNnwdK3pqZmj8bV1tb6uZJLgqVvvWEymbzeGQAAT3kVRqKjoxUaGtrlKMi5c+e6HC3pEBsbq6uvvtoZRCRpxIgRcjgc+te//uV279VsNstsNrt9Pm829g6HY8D+cfCnwO+bZ7X192sI/L4BQGDy6jojJpNJGRkZqqiocJlfUVGhrKwst8uMHTtWDQ0NunDhgnNebW2tQkJCNHz4cB9KBgAAA4nXFz3Lzc3VG2+8oX379unkyZPatGmTLBaLbrnlFklSSUmJ1q9f7xz/pS99SVFRUfrpT3+qkydPqrKyUlu3btWMGTO6PYEVAAAMHl6fMzJ16lQ1NTVpx44damhoUFpamlasWOH8PLmhoUEWi8U5fujQoXrssce0ceNGPfroo4qKitKNN96oRYsW9d2rAAAAQcunE1jnzJmjOXPmuH2soKCgy7wRI0bo8ccf92VVAABggOPeNAAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGMhldAAAMRq9U1BtdAhAwCCMIWo7j713x8faaWq+fM3TeHb6WAwDwER/TAAAAQxFGAACAoQgjAADAUIQRAABgKMIIBqxtrcna1ppsdBkAgB4QRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKJMvC+3Zs0c7d+6U1WpVamqq8vPzNW7cuB6XO3bsmFatWqW0tDT9+Mc/9mXVAABggPH6yMihQ4e0adMm3XrrrSoqKtK4ceO0Zs0aWSyWKy53/vx5bdiwQRMnTvS5WAAAMPB4HUZ2796tmTNnatasWc6jIvHx8SovL7/icr/85S910003acyYMT4XCwAABh6vPqax2+2qqqrSggULXOZnZ2fr+PHj3S63f/9+nT59WsuWLdOOHTt6XI/NZpPNZnNOh4SEKDw83PlzTzrGeDIWnYKnb97W5/l4X1578PQNAAKTV2GksbFR7e3tiomJcZkfExMjq9Xqdpna2lqVlJSosLBQQ4YM8Wg9ZWVlKi0tdU6np6erqKhICQkJ3pSrpKQkr8bjkkDvW1RViyTpQliYZ+OjIj1+7pjkZJ9qkgK/bwAQqHw6gdXdHqC7ee3t7frJT36ihQsXKiUlxePnz8vLU25ubpfnrq+vl91u96i+pKQk1dXVyeFweLzewS5Y+tbU1CxJcly86NV4T5yvrfW6nmDpW2+YTCavdwYAwFNehZHo6GiFhoZ2OQpy7ty5LkdLJKmlpUUnTpxQdXW1Nm7cKElyOBxyOBxatGiRHnvsMV133XVdljObzTKbzW5r8GZj37EueCfw++Ztbd69Z3wV+H0DgMDkVRgxmUzKyMhQRUWFpkyZ4pxfUVGhL37xi13Gh4eH69lnn3WZV15ervfee08PPvigEhMTfSwbAAAMFF5/TJObm6sXXnhBGRkZyszM1N69e2WxWHTLLbdIkkpKSnT27Fnde++9Cg0N1ciRI12Wj46Oltls7jIfAAAMTl6HkalTp6qpqUk7duxQQ0OD0tLStGLFCufnyQ0NDT1ecwQAAKCDTyewzpkzR3PmzHH7WEFBwRWXvf3223X77bf7sloAADAAcW8aAABgKJ+OjACAt7y5p9X777+vwsLCLvPXrVunESNG+LtUAP2MMALA7zruabVkyRJlZWVp7969WrNmjdatW6f4+Phul3v++ecVERHhnI6Oju6PcgOW4/h7Xea113h/bZwOofPu6E05QJ/hYxoAfufrPa1iYmIUGxvr/BcayiYLGIg4MgLAr3y9p5UkPfLII7LZbEpNTdWtt97q9iKJHXy5p5Wx9xXqq3X6/jyffd3cY6kTvejUX70gjADwK1/uaRUXF6e7775bGRkZstvtevPNN/XUU0/pySef1Pjx490u05t7WhlxX6GOeyx5w939mLy599Ll3N2LiXssdaIXnfzdC8IIgH7h6T2tJCklJcXlflaZmZmyWCzatWtXt2HEl3taGXlfIW/umdTB3f2YfHmeDp+9F9NguMeSp+hFp972wtP7WhFGAPiVt/e06k5mZqYOHDjQ7eO9uaeVMfcV6qv19e5+Su7mDfY/wB3oRSd/94KzwQD41WfvafVZFRUVysrK8vh5qqurFRsb28fVAQgEHBkB4Hfe3NNKkn73u98pISFBaWlpstvtOnDggA4fPqyHHnrIyJcBwE8IIwD8ztt7Wtntdr388ss6e/aswsLClJaWpkcffVSTJ0826iUErG2triehLrrK9+uOAEYhjADoF97c02r+/PmaP39+f5Tl1L6zpE+fjwuKAZ7jnBEAAGAojowgaLxSUW90CQAAP+DICAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChTEYXAAADUfvOkis+7mhN7qdKgMDHkREAAGAowggAADAUYQQAABjKp3NG9uzZo507d8pqtSo1NVX5+fkaN26c27GHDx9WeXm5ampqZLfblZqaqoULF2rSpEm9qRsAAAwQXh8ZOXTokDZt2qRbb71VRUVFGjdunNasWSOLxeJ2/AcffKDs7GytWLFCzzzzjCZMmKCioiJVV1f3ungAABD8vA4ju3fv1syZMzVr1iznUZH4+HiVl5e7HZ+fn6/58+fr2muvVXJysu644w4lJyfrnXfe6XXxAAAg+Hn1MY3dbldVVZUWLFjgMj87O1vHjx/36Dna29vV0tKiyMjIbsfYbDbZbDbndEhIiMLDw50/96RjjCdj0Snw++ZrXZ4v58trD/y+AUBg8yqMNDY2qr29XTExMS7zY2JiZLVaPXqO3bt3q7W1VTfeeGO3Y8rKylRaWuqcTk9PV1FRkRISErwpV0lJSV6NxyWB2reoqhaX6QthYZ4tF9V98L1cTLLv134I1L4BQKDz6QRWd3uAnuwVHjx4UL/5zW/08MMPdwk0n5WXl6fc3Nwuz11fXy+73e5RfUlJSaqrq5PD4ehxPC4J9L41NTW7TDsuXvRpuSs5X1vrVU1S4PetL5hMJq93BgDAU16FkejoaIWGhnY5CnLu3Lkrhgvp0omvP//5z/Xggw8qOzv7imPNZrPMZrPbx7zZ2DscjgH7x8GfArdvvtbk3XvGV4HbNwAIbF6dwGoymZSRkaGKigqX+RUVFcrKyup2uYMHD2rDhg267777NHnyZN8qBQAAA5LX36bJzc3VG2+8oX379unkyZPatGmTLBaLbrnlFklSSUmJ1q9f7xzfEUS+853vKDMzU1arVVarVefPn++7VwEAAIKW1+eMTJ06VU1NTdqxY4caGhqUlpamFStWOD9PbmhocLnmyN69e9XW1qbi4mIVFxc750+bNk0FBQV98BIAAEAw8+kE1jlz5mjOnDluH7s8YKxatcqXVQAAgEGCe9MAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAzl03VGgIGqfWeJD0uF6FxUpNqamnX5fXBC593RJ3UBwEDGkREAAGAojowAwCDleiSw+yN83uBoIHzBkREAAGAowggAADAUYQQD3rbWZG1rTTa6DABANzhnBAAGkM8G70VX1RpYCeA5jowAAABDcWQEAAaoyz+e5EgJAhVHRgAAgKE4MoJ+49vVTTs5OAkVAAYkjowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAU36ZBwONS7gAwsHFkBAAAGIowAgAADMXHNAAwSHATPQQqwggAoM/09krLlwudd0efPh8CEx/TAAAAQxFGMGhsa03mmzkAEID4mAZu9fWhVgAAukMYAQAELH/sGHEeSuDxKYzs2bNHO3fulNVqVWpqqvLz8zVu3Lhux1dWVmrz5s06efKk4uLiNG/ePM2ePdvnogEEH7YbALrjdRg5dOiQNm3apCVLligrK0t79+7VmjVrtG7dOsXHx3cZf+bMGT399NOaNWuWli1bpuPHj+vFF19UdHS0brjhhj55EQACG9sNBJKej7aE6FxUpNqamiU5enw+jrT0ntdhZPfu3Zo5c6ZmzZolScrPz9fRo0dVXl6uO+7o+gspLy9XfHy88vPzJUmpqak6ceKEdu3axUYFA14wnHvTHxtSthuB5/KTuRddVWdQJYCXYcRut6uqqkoLFixwmZ+dna3jx4+7Xeajjz5Sdna2y7xJkyZp//79stvtMpm6lmCz2WSz2ZzTISEhCg8PdzvWnZCQEEmS2WyWw9Fzqh0I2v/0eq+fIyREagqP0JCW8wqN/7c+qKpvJF6M69PnM4e19+nzhYRIoeERChsaoWB8u4WazT2O8fT/njuBvN347Laiv9/zff2+7q2wq9qD+n3cl7z+P33oDb/XZJTP/l24vBeh077a4/Kebju82sI0Njaqvb1dMTExLvNjYmJktVrdLmO1Wt2Ob2trU1NTk+Liuv6HLCsrU2lpqXP6pptu0v333+927JW4O/w7YN327T57qqg+e6a+8R2jC/BQoPUtUATDdiM+Pr5P/w95IlDf17yPO9GLTv7uhU/XGenYm+hpXnePdRyt6G6ZvLw8bdq0yfnvP/7jP1z2eHrS0tKiH/zgB2ppafF4GdA3X9E3zwTidoPfXSd60YledOqvXnh1ZCQ6OlqhoaFd9mbOnTvXZS+mQ2xsbJfxjY2NGjJkiCIjI90uYzabZfbg0HF3HA6HqqurB81HNH2FvvmGvl1ZIG83+N11ohed6EWn/uqFV0dGTCaTMjIyVFFR4TK/oqJCWVlZbpcZM2ZMl/FHjx5VRkZGrz6HBhAc2G4A6InXH9Pk5ubqjTfe0L59+3Ty5Elt2rRJFotFt9xyiySppKRE69evd46fPXu2LBaL83oB+/bt0759+zR37ty+exUAAhrbDQBX4vUuxtSpU9XU1KQdO3aooaFBaWlpWrFihRISEiRJDQ0NslgszvGJiYlasWKFNm/erD179iguLk533XWXX7+eZzabddttt/Xqo57BiL75hr71LFC3G/zuOtGLTvSiU3/1IsTBh2IAAMBA3LUXAAAYijACAAAMRRgBAACGIowAAABDDfgv7BcUFKi+vt5l3vz58/Wtb33LoIoCk7e3d4e0fft2l8uPS5cuWf6rX/3KoIrgDd7zly6hf+TIEZ06dUphYWHKzMzUnXfeqZSUFKNLM1xZWZleeeUVff3rX3fesHGwOXv2rLZu3ap3331XFy9eVHJysr7//e8rIyOjz9c14MOIJN1+++26+eabndNDhw41sJrA4+3t3dEpLS1Njz/+uHM6NJSDjcGA9/wllZWVmjNnjkaPHq22tjZt27ZNP/zhD7V27dpBvZ38+OOPtXfvXl1zzTVGl2KY5uZmPf7445owYYL+67/+S9HR0Tp9+rQiIiL8sr5BseUMDw9XbGys899g/k/mzmdv796xhxgfH6/y8nKjSwt4oaGhLu+t6Ohoo0uCB3jPX7Jy5UpNnz5daWlpGjVqlJYuXSqLxaKqqiqjSzPMhQsX9MILL+iee+7RsGHDjC7HMK+99pqGDx+upUuX6tprr1ViYqImTpyopKQkv6xvUBwZee2117Rjxw4NHz5cN954o+bNm8clpf8/X27vjk51dXW65557ZDKZNGbMGH3zm9/Uv/1b/96KHt7hPd+98+fPS1K39/8ZDF588UXl5OQoOztbr776qtHlGObtt9/W5z73Oa1du1aVlZW6+uqrNXv2bJdPGfrSgP+L/LWvfU0ZGRkaNmyYPv74Y5WUlOjMmTP63ve+Z3RpAcGX27vjkjFjxqigoEApKSmyWq169dVX9dhjj2nt2rWKiuLm44GK97x7DodDmzdv1tixYzVy5EijyzHE//7v/6q6ulpPP/200aUY7syZM/rDH/6gb3zjG8rLy9PHH3+sl156SWazWdOmTevz9QVlGHF34uDlnn76aY0ePVq5ubnOeddcc42GDRumtWvX6lvf+hZ/MD7D29u7Q8rJyXH+PHLkSGVmZmrZsmX605/+5PK+Q2DiPe+quLhY//jHP7R69WqjSzGExWLRpk2btHLlSoWFhRldjuHa29s1evRo3XHHHZKk9PR0ffLJJyovLyeMdPjqV7+qm2666YpjOu55cbnMzExJlw6vE0Z8u7073Bs6dKhGjhyp2tpao0vBFfCe72rjxo165513VFhYqOHDhxtdjiGqqqp07tw5Pfroo8557e3t+uCDD/T666+rpKRkUJ2gHhcXp9TUVJd5qampOnz4sF/WF5RhJDo62ucTBaurqyVdajRcb+8+ZcoU5/yKigp98YtfNLCy4GOz2XTq1KlB9/XQYMN7vpPD4dDGjRt15MgRrVq1SomJiUaXZJiJEyfq2WefdZn3s5/9TCkpKZo/f/6gCiKSlJWVpX/+858u8/75z392u6PfW0EZRjz14Ycf6sMPP9R1112niIgIffzxx9q8ebO+8IUvDKqv7/UkNzdXL7zwgjIyMpSZmam9e/e63N4d7m3ZssX5Xjp37px27NihlpYWvxzCRN/iPX9JcXGxDh48qEceeUTh4eHOo0URERGD7qOK8PDwLufKXHXVVYqKihqU59B84xvf0OOPP65XX31VU6dO1ccff6w33nhDd999t1/WN6Dv2ltVVaXi4mKdOnVKNptNCQkJmjp1qubPn6+rrrrK6PICSscFoDpu7/7d735X48ePN7qsgPb888/rgw8+UGNjo6KjozVmzBgtWrSoy6FNBCbe85euweTO0qVLNX369P4tJgCtWrVKo0aNGrQXPXvnnXdUUlKiuro6JSYm6hvf+Ibfvk0zoMMIAAAIfIPrQzAAABBwCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAfrVhwwYVFBQYXQYCGFdgBQD4VV1dnVpaWpSenm50KQhQhBEAAGCoAX3XXvjX9u3bVVpaqh//+MfasWOHjh49qtDQUE2fPl133nmnTp8+rZdeeknHjx9XVFSUZs+erfnz50uS/vjHP+qnP/2p1q9f73Lb8vfff1+FhYV68sknNWHCBKNeGgAvNDY26pVXXtG7776rc+fOKTw8XCkpKVq4cKGys7O1YcMGVVZWasOGDc5lPv30U23ZskVHjhyR3W7X+PHjtXjxYi1btky33Xab8yZ+vdnOSNLFixe1bds2/e1vf9OZM2cUGhqqlJQULViwQF/84hf7vVdwjzCCXlu3bp2+/OUv6+abb1ZFRYV27typtrY2/e1vf9Ps2bM1d+5cHTx4UL/+9a+VlJSk66+/3uiSAfShF154QdXV1Vq0aJFSUlL06aefqrq6Ws3NzW7Ht7e3q6ioSCdOnNDChQuVkZGhDz/8UGvWrOl2Hb5uZ+x2u5qbmzV37lxdffXVstvt+tvf/qZnn31WS5cu1bRp0/zSE3iHMIJeu/nmm5WbmytJys7OVkVFhV5//XUtX75cU6ZMkSRNmDBBf/nLX3TgwAHCCDDAHD9+XDNnznS5vfyVjjq8++67OnbsmJYsWaLZs2dLurTtMJlMKikpcbuMr9uZiIgILV261Pk87e3tmjhxoj799FP9/ve/J4wECMIIem3y5Mku0yNGjNDf//53TZo0yTlvyJAhSkpKksVi6efqAPjbtddeqz/96U+KiorSxIkTlZGRIZOp+z8vlZWVkqSpU6e6zL/pppu6DSO92c78+c9/1u9//3vV1NSotbXVOd9sNnv0+uB/hBH0WmRkpMu0yWRSWFiYwsLCusxvaWnpz9IA9IMHHnhAr776qvbt26f//u//1tChQzVlyhTdeeedio2N7TK+ublZQ4YM6bLtiImJ6XYdvm5nDh8+rHXr1umGG27Q3LlzFRsbqyFDhqi8vFz79+/34dXCHwgjMETHHondbneZ39TUZEQ5AHohOjpa+fn5ys/Pl8Vi0dtvv61f//rXOnfunFauXNllfGRkpNra2tTc3OwSMqxWa5/XduDAASUmJuo///M/FRIS4pxvs9n6fF3wHRc9gyESEhIkSX//+99d5r/99ttGlAOgj8THx+urX/2qsrOzVV1d7XbM+PHjJUmHDh1ymX/5dF8xmUwuQcRqtbKtCTAcGYEhrr32WqWkpOjll19WW1ubIiMjdeTIER07dszo0gB44fz58yosLNRNN92kESNGKDw8XB9//LHefffdbk9WnzRpkrKysrRlyxadP3/e+W2aN998U5IUGtp3+8mf//zndeTIEb344ou64YYbZLFYtGPHDsXFxam2trbP1oPeIYzAEKGhofrBD36gjRs36le/+pXMZrOmTp2qxYsX65lnnjG6PAAeMpvNuvbaa3XgwAGdOXNGbW1tio+P1/z5812u9/FZHf//t2zZotdee012u11ZWVlatmyZVq5cqYiIiD6rb8aMGTp37pz+8Ic/aP/+/UpMTNSCBQv0r3/9S6WlpX22HvQOV2AFAASEgwcP6ic/+YmeeuopZWVlGV0O+hFHRgAA/e7gwYM6e/asRo4cqdDQUH344YfatWuXxo0bRxAZhAgjAIB+Fx4erkOHDunVV19Va2urYmNjNW3aNC1atMjo0mAAPqYBAACG4qu9AADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAICh/h/hMWF4LKPlgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use('ggplot')\n",
    "for i, p in enumerate(['mu', 'sigma']):\n",
    "    plt.subplot(1, 2, i + 1)\n",
    "    for t in [1, 29]:\n",
    "        plt.hist(my_alg.hist.X[t].theta[p], weights=my_alg.hist.wgts[t].W, label=\"t=%i\" % t, \n",
    "                 alpha=0.5, density=True)\n",
    "    plt.xlabel(p)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the posterior distribution concentrates progressively around the true values. \n",
    "\n",
    "As always, once the algorithm is run, `my_smc.X` contains the final particles. However, object `my_smc.X` is no longer a simple numpy array. It is a `ThetaParticles` object, with attributes:\n",
    "\n",
    "* `theta`: a structured array (an array with fields); i.e. `my_smc.X.theta['mu']` is a (N,) array that contains the the $\\mu-$component of the $N$ particles; \n",
    "* `lpost`: a 1D numpy array that contains the target (posterior) log-density of each of the particles;\n",
    "* `shared`: a dictionary that contains \"meta-data\" on the particles; for instance `shared['acc_rates']` is a list of the acceptance rates of the successive Metropolis steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['18%', '26%', '24%', '36%', '35%', '35%', '32%']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGdCAYAAADXIOPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAq50lEQVR4nO3df3DU9YH/8dcmu4GEhGzIDxIsaRIMCFFIPARE7pBfajVfaVRisVQtxLtrKAOd46SKWmCggFNFe0O9q+IBuYKESBQrM8cULBJRqcoctLmBKQETlL0EyRp+heyS/f7BZOuakGTDbvad3edjhpl89vNj3y8+CXnx3s9+1uLxeDwCAAAwTFSoBwAAANARSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMJI11AO4Ho2NjXK73aEeRlCkpqaqoaEh1MPoVWSODGSODGSODP5mtlqtSkpK6v72PRmUKdxut1wuV6iHEXAWi0XS1XyR8tFKZCZzuCIzmcNVb2Tm5R4AAGAkv2ZSKisrdfDgQX3xxReKiYnR8OHDNWfOHA0ZMsS7jcfj0fbt27Vnzx6dP39eubm5mjdvnoYOHerdxuVyqaysTB988IFaWlp08803q6SkRMnJyYFLBgAA+jS/ZlKqq6t19913a9WqVXrmmWfU2tqqlStXqrm52bvN22+/rXfffVdz587V6tWrZbfbtXLlSl26dMm7zcaNG3Xw4EEtXLhQK1asUHNzs9asWaPW1tbAJQMAAH2aXyVl6dKluvPOOzV06FBlZWWptLRUZ86cUU1NjaSrsyi7du1SUVGRxo8fr8zMTM2fP1+XL19WVVWVJOnixYvau3evHn30UY0ePVrZ2dlasGCBamtrdfjw4cAnBAAAfdJ1XTh78eJFSVJ8fLwkqb6+Xk6nU2PGjPFuY7PZNGrUKB09elQzZsxQTU2Nrly5otGjR3u3GTRokDIzM3Xs2DHl5+e3ex6Xy+VzgazFYlFsbKz363DTlikcs10LmSMDmSMDmSNDb2TucUnxeDzatGmTbrrpJmVmZkqSnE6nJCkxMdFn28TERJ05c8a7jdVq9Rabb27Ttv+3VVZWqqKiwrucnZ2ttWvXKjU1tafD7xPS09NDPYReR+bIQObIQObIEMzMPS4pGzZsUG1trVasWNFu3bdbVXfemtTZNkVFRSosLGx3/IaGhrC8T4rFYlF6erocDkdEvZWNzOGPzGQOV2TuXmar1erXBEOPSsrrr7+uTz/9VMuXL/d5R47dbpd0dbbkmzdraWpq8s6u2O12ud1unT9/3mc2pampSSNGjOjw+Ww2m2w2W4frwvmbwePxhHW+jpA5MpA5MpA5MgQzs18Xzno8Hm3YsEEff/yxnnvuOaWlpfmsT0tLk91u97kA1u12q7q62ltAcnJyFB0d7bNNY2OjamtrNXz48OvJAgAAwohfMykbNmxQVVWVnnzyScXGxnqvIYmLi1NMTIwsFovuvfdeVVZWKiMjQ+np6aqsrFS/fv00adIk77ZTp05VWVmZEhISFB8fr7KyMmVmZvpcTAsAACKbXyVl9+7dkqRly5b5PF5aWqo777xTkjRz5ky1tLTotdde04ULF3TjjTdq6dKl3nfjSNJjjz2m6OhorVu3znsztyVLligqihvgAgCAq/wqKeXl5V1uY7FYVFxcrOLi4mtuExMTo7lz52ru3Ln+PD0AAIggTF0AAAAjUVIAAICRruuOswAAILiuPHF/j/eNfnVnAEfS+5hJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYyervDtXV1dq5c6dOnDihxsZGLV68WOPGjfOuLy4u7nC/OXPm6P7775ckLVu2TNXV1T7rJ06cqEWLFvk7HAAAEKb8LimXL19WVlaWpkyZohdeeKHd+t/+9rc+y4cOHdK///u/a/z48T6PT5s2TQ8//LB3OSYmxt+hAACAMOZ3SSkoKFBBQcE119vtdp/lP/3pT8rLy9PgwYN9Hu/Xr1+7bQEAANr4XVL84XQ6dejQIc2fP7/duv3792v//v1KTExUfn6+Zs2apdjY2A6P43K55HK5vMsWi8W7rcViCc7gQ6gtUzhmuxYyRwYyRwYymyOY4+mNzEEtKfv27VP//v19rlmRpEmTJiktLU12u111dXXasmWLPv/8cz377LMdHqeyslIVFRXe5ezsbK1du1apqanBHH7Ipaenh3oIvY7MkYHMkYHMgVF3HftmZGQEbBzXEszzHNSS8t577+nv//7v211vMn36dO/XmZmZysjI0M9//nPV1NQoJyen3XGKiopUWFjoXW5rbQ0NDXK73UEafehYLBalp6fL4XDI4/GEeji9gsxkDldkJnMonT59OmjH7klmq9Xq1wRD0ErK//7v/+rLL7/s1jt2srOzFR0dLYfD0WFJsdlsstlsHe5r0jdDoHk8nrDO1xEyRwYyRwYyh15vjCWYmYN2n5S9e/cqJydHWVlZXW5bV1enK1eucCEtAADw8nsmpbm5WQ6Hw7tcX1+vkydPKj4+XikpKZKkixcv6qOPPtKPfvSjdvs7HA5VVVWpoKBACQkJOnXqlMrKypSdna2bbrrpOqIAAIBw4ndJOX78uJYvX+5d3rx5syRp8uTJ3nfxHDhwQB6PR5MmTWr/hFarjhw5ol27dqm5uVnJycm69dZbNWvWLEVFcQNcAABwld8lJS8vT+Xl5Z1uM336dJ+LY78pJSXFp+QAAAB0hKkLAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADCS1d8dqqurtXPnTp04cUKNjY1avHixxo0b512/fv167du3z2ef3NxcrVq1yrvscrlUVlamDz74QC0tLbr55ptVUlKi5OTk64gCAADCid8l5fLly8rKytKUKVP0wgsvdLhNfn6+SktL//YkVt+n2bhxoz799FMtXLhQCQkJ2rx5s9asWaO1a9cqKorJHQAA0IOSUlBQoIKCgs4ParXKbrd3uO7ixYvau3evFixYoNGjR0uSFixYoJ/85Cc6fPiw8vPz/R0SAAAIQ36XlO6orq5WSUmJBgwYoJEjR2r27NlKTEyUJNXU1OjKlSvegiJJgwYNUmZmpo4dO0ZJAQAAkoJQUgoKCnT77bcrJSVF9fX12rZtm1asWKE1a9bIZrPJ6XTKarUqPj7eZ7/ExEQ5nc4Oj+lyueRyubzLFotFsbGx3q/DTVumcMx2LWSODGSODGQ2RzDH0xuZA15SJk6c6P06MzNTw4YNU2lpqT777DONHz/+mvt5PJ5rrqusrFRFRYV3OTs7W2vXrlVqampgBm2o9PT0UA+h15E5MpA5MpA5MOquY9+MjIyAjeNagnmeg/JyzzclJSUpNTVVp0+fliTZ7Xa53W6dP3/eZzalqalJI0aM6PAYRUVFKiws9C63tbaGhga53e4gjj40LBaL0tPT5XA4Oi1v4YTMZA5XZCZzKNXdN7bH+1pfe6fT9T3JbLVa/ZpgCHpJOXfunL766islJSVJknJychQdHa3Dhw97Z10aGxtVW1urH/7whx0ew2azyWazdbjOpG+GQPN4PGGdryNkjgxkjgxk7tu6myOYmf0uKc3NzXI4HN7l+vp6nTx5UvHx8YqPj1d5ebkmTJggu92uhoYGbd26VQkJCd57qcTFxWnq1KkqKytTQkKC4uPjVVZWpszMTJ+LaQEAQGTzu6QcP35cy5cv9y5v3rxZkjR58mQ98cQTqqur0/vvv68LFy4oKSlJeXl5WrRokfdCV0l67LHHFB0drXXr1nlv5rZkyRLukQIAALz8Lil5eXkqLy+/5vqlS5d2eYyYmBjNnTtXc+fO9ffpAQBAhGDqAgAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMZPV3h+rqau3cuVMnTpxQY2OjFi9erHHjxkmS3G633njjDR06dEj19fWKi4vTLbfcokceeUSDBg3yHmPZsmWqrq72Oe7EiRO1aNGi60sDAADCht8l5fLly8rKytKUKVP0wgsv+KxraWnRiRMn9OCDDyorK0vnz5/Xpk2b9Pzzz2vNmjU+206bNk0PP/ywdzkmJqaHEQAAQDjyu6QUFBSooKCgw3VxcXF69tlnfR778Y9/rKefflpnzpxRSkqK9/F+/frJbrf7+/QAACBC+F1S/HXx4kVZLBbFxcX5PL5//37t379fiYmJys/P16xZsxQbGxvs4QAAgD4iqCWlpaVFW7Zs0R133OFTUiZNmqS0tDTZ7XbV1dVpy5Yt+vzzz9vNwrRxuVxyuVzeZYvF4i00FoslmBFCoi1TOGa7FjJHBjJHBjKHh66y9EbmoJUUt9utl156SR6PRyUlJT7rpk+f7v06MzNTGRkZ+vnPf66amhrl5OS0O1ZlZaUqKiq8y9nZ2Vq7dq1SU1ODNXwjpKenh3oIvY7MkYHMkYHMgVEX8CN2T0ZGRre2C+Z5DkpJcbvdWrdunRoaGvTcc8+1e6nn27KzsxUdHS2Hw9FhSSkqKlJhYaF3ua21NTQ0yO12B3bwBrBYLEpPT5fD4ZDH4wn1cHoFmckcrshM5r7q9OnTna7vSWar1erXBEPAS0pbQXE4HPrFL36hhISELvepq6vTlStXrnkhrc1mk81m63BduHwzdMTj8YR1vo6QOTKQOTKQuW/rbo5gZva7pDQ3N8vhcHiX6+vrdfLkScXHxyspKUkvvviiTpw4oSVLlqi1tVVOp1OSFB8fL6vVKofDoaqqKhUUFCghIUGnTp1SWVmZsrOzddNNNwUsGAAA6Nv8LinHjx/X8uXLvcubN2+WJE2ePFmzZs3SJ598Ikl68sknffb7xS9+oby8PFmtVh05ckS7du1Sc3OzkpOTdeutt2rWrFmKiuIGuAAA4Cq/S0peXp7Ky8uvub6zdZKUkpLiU3IAAAh3V564P9RD6JOYugAAAEaipAAAACNRUgAAgJGCflt8AADCQVfXlYTqpmvhjJkUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJz+4BAESMrj5/B2ZhJgUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGsvq7Q3V1tXbu3KkTJ06osbFRixcv1rhx47zrPR6Ptm/frj179uj8+fPKzc3VvHnzNHToUO82LpdLZWVl+uCDD9TS0qKbb75ZJSUlSk5ODkwqAADQ5/k9k3L58mVlZWVp7ty5Ha5/++239e6772ru3LlavXq17Ha7Vq5cqUuXLnm32bhxow4ePKiFCxdqxYoVam5u1po1a9Ta2trzJAAAIKz4XVIKCgr0gx/8QOPHj2+3zuPxaNeuXSoqKtL48eOVmZmp+fPn6/Lly6qqqpIkXbx4UXv37tWjjz6q0aNHKzs7WwsWLFBtba0OHz58/YkAAEBY8Pvlns7U19fL6XRqzJgx3sdsNptGjRqlo0ePasaMGaqpqdGVK1c0evRo7zaDBg1SZmamjh07pvz8/HbHdblccrlc3mWLxaLY2Fjv1+GmLVM4ZrsWMkcGMkeGSMwcjro6f71xngNaUpxOpyQpMTHR5/HExESdOXPGu43ValV8fHy7bdr2/7bKykpVVFR4l7Ozs7V27VqlpqYGbvAGSk9PD/UQeh2ZIwOZI4OJmetCPYA+JCMjo1vbBfM8B7SktPl2q/J4PF3u09k2RUVFKiwsbHf8hoYGud3uHo7SXBaLRenp6XI4HN36uwsHZCZzuCJzZGQOR6dPn+50fU/Os9Vq9WuCIaAlxW63S7o6W5KUlOR9vKmpyTu7Yrfb5Xa7df78eZ/ZlKamJo0YMaLD49psNtlstg7XhfMPgMfjCet8HSFzZCBzZIjEzOGku+cumOc5oPdJSUtLk91u97kA1u12q7q62ltAcnJyFB0d7bNNY2OjamtrNXz48EAOBwAA9GF+z6Q0NzfL4XB4l+vr63Xy5EnFx8crJSVF9957ryorK5WRkaH09HRVVlaqX79+mjRpkiQpLi5OU6dOVVlZmRISEhQfH6+ysjJlZmb6XEwLAAAim98l5fjx41q+fLl3efPmzZKkyZMna/78+Zo5c6ZaWlr02muv6cKFC7rxxhu1dOlS77txJOmxxx5TdHS01q1b572Z25IlSxQVxQ1wAQDAVX6XlLy8PJWXl19zvcViUXFxsYqLi6+5TUxMjObOnXvNG8IBAAAwdQEAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjBTQT0EGACDYrjxxf6iHgF7CTAoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAka6APOH/+fDU0NLR7/K677lJJSYnWr1+vffv2+azLzc3VqlWrAj0UAADQhwW8pKxevVqtra3e5draWq1cuVK3336797H8/HyVlpb+bRDWgA8DAAD0cQFvBwMHDvRZfuuttzR48GCNGjXqb09qtcputwf6qQEAQBgJ6hSG2+3W/v37dd9998lisXgfr66uVklJiQYMGKCRI0dq9uzZSkxMvOZxXC6XXC6Xd9lisSg2Ntb7dbhpyxSO2a6FzJGBzJEhEjOHo67OX2+cZ4vH4/EE6+AHDhzQr3/9a/3mN7/RoEGDvI/1799fKSkpqq+v17Zt29Ta2qo1a9bIZrN1eJzy8nJVVFR4l7Ozs7V27dpgDRsAYLC6+8aGeggRYei7n4R6CMGdSXnvvfeUn5/vLSiSNHHiRO/XmZmZGjZsmEpLS/XZZ59p/PjxHR6nqKhIhYWF3uW21tbQ0CC32x2k0YeOxWJRenq6HA6HgtghjUJmMocrMkdG5nB0+vTpTtf35DxbrValpqZ2ewxBKykNDQ06fPiwFi9e3Ol2SUlJSk1N7fQvw2azXXOWJZx/ADweT1jn6wiZIwOZI0MkZg4n3T13wTzPQbtPynvvvafExETdeuutnW537tw5ffXVV0pKSgrWUAAAQB8UlJmU1tZW/fGPf9TkyZMVHR3tfby5uVnl5eWaMGGC7Ha7GhoatHXrViUkJGjcuHHBGAoAAOijglJSjhw5ojNnzmjKlCk+j0dFRamurk7vv/++Lly4oKSkJOXl5WnRokXed+sAAABIQSopY8aMUXl5ebvHY2JitHTp0mA8JQAACDN8dg8AADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGMka6AOWl5eroqLC57HExES9+uqrkiSPx6Pt27drz549On/+vHJzczVv3jwNHTo00EMBAAB9WMBLiiQNHTpUzz77rHc5KupvEzZvv/223n33XZWWliojI0M7duzQypUr9dJLLyk2NjYYwwEAAH1QUF7uiYqKkt1u9/4ZOHCgpKuzKLt27VJRUZHGjx+vzMxMzZ8/X5cvX1ZVVVUwhgIAAPqooMykOBwO/dM//ZOsVqtyc3M1e/ZsDR48WPX19XI6nRozZox3W5vNplGjRuno0aOaMWNGh8dzuVxyuVzeZYvF4p11sVgswYgQUm2ZwjHbtZA5MpA5MkRi5nDU1fnrjfMc8JKSm5ur+fPna8iQIXI6ndqxY4eeeeYZvfjii3I6nZKuXqPyTYmJiTpz5sw1j1lZWelznUt2drbWrl2r1NTUQA/fKOnp6aEeQq8jc2Qgc2QIVua6oBwV35aRkdGt7YL5vR3wklJQUOD9OjMzU8OHD9eCBQu0b98+5ebmSmrfujweT6fHLCoqUmFhoXe5bf+Ghga53e5ADd0YFotF6enpcjgcXf7dhAsykzlckTkyMoej06dPd7q+J+fZarX6NcEQlJd7vql///7KzMzU6dOnddttt0mSnE6nkpKSvNs0NTW1m135JpvNJpvN1uG6cP4B8Hg8YZ2vI2SODGSODJGYOZx099wF8zwH/T4pLpdLX3zxhZKSkpSWlia73a7Dhw9717vdblVXV2vEiBHBHgoAAOhDAj6TsnnzZo0dO1YpKSn6+uuv9eabb+rSpUuaPHmyLBaL7r33XlVWViojI0Pp6emqrKxUv379NGnSpEAPBQAA9GEBLylnz57Vyy+/rKamJg0cOFC5ublatWqV9zWomTNnqqWlRa+99pouXLigG2+8UUuXLuUeKQAAwEfAS8qiRYs6XW+xWFRcXKzi4uJAPzUAAAgjfHYPAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIwU8M/uAQCgK1eeuD/UQ0AfwEwKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjUVIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjERJAQAARqKkAAAAI1FSAACAkayBPmBlZaUOHjyoL774QjExMRo+fLjmzJmjIUOGeLdZv3699u3b57Nfbm6uVq1aFejhAACAPirgJaW6ulp33323hg0bpitXruiNN97QypUr9eKLL6p///7e7fLz81VaWvq3gVgDPhQAANCHBbwZLF261Ge5tLRUJSUlqqmp0ahRo/72xFar7HZ7oJ8eACLKlSfuv679o1/dGaCRAIEX9OmLixcvSpLi4+N9Hq+urlZJSYkGDBigkSNHavbs2UpMTOzwGC6XSy6Xy7tssVgUGxvr/TrctGUKx2zXQubIQGbzBGNcpmdG93R1/nrjPFs8Ho8nWAf3eDx6/vnndeHCBa1YscL7+IEDB9S/f3+lpKSovr5e27ZtU2trq9asWSObzdbuOOXl5aqoqPAuZ2dna+3atcEaNgD0GXX3jb2u/Ye++0mARuKf6x03gi9U3xvfFNSZlA0bNqi2ttanoEjSxIkTvV9nZmZq2LBhKi0t1Weffabx48e3O05RUZEKCwu9y22traGhQW63O0ijDx2LxaL09HQ5HA4FsUMahcxkDlemZz59+nTAj2l6ZnRPV98bPTnPVqtVqamp3R5D0ErK66+/rk8//VTLly9XcnJyp9smJSUpNTX1mn8hNputwxkWSWH9A+DxeMI6X0fIHBnIbA53yf/r8b5dXc9iamZ0T3fPXTDPc8Dvk+LxeLRhwwZ9/PHHeu6555SWltblPufOndNXX32lpKSkQA8HAAD0UQGfSdmwYYOqqqr05JNPKjY2Vk6nU5IUFxenmJgYNTc3q7y8XBMmTJDdbldDQ4O2bt2qhIQEjRs3LtDDAQAESWfvLKrrxXEgfAW8pOzevVuStGzZMp/HS0tLdeeddyoqKkp1dXV6//33deHCBSUlJSkvL0+LFi3yvmMHAAAg4CWlvLy80/UxMTHt7qUCAADwbXx2DwAAMBIlBQAAGImSAgAAjERJAQAARuKjhwEgAK7n7bh8yB/QMWZSAACAkSgpAADASJQUAABgJEoKAAAwEiUFAAAYiZICAACMREkBAABGoqQAAAAjcTM3AAixzm4EB0QyZlIAAICRKCkAAMBIlBQAAGAkSgoAADASJQUAABiJkgIAAIzEW5ABGOd63pIb/erOkDwvgMBjJgUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJF4dw+AoLjWO2XqenkcAPouSgrQSzr6pd3dX9ihelvt9TwvAFwvXu4BAABGYiYFIfPt/+H78zIAMwsAEP6YSQEAAEaipAAAACNRUgAAgJG4JgXoA/jgu+7j7woIHyEtKf/93/+tnTt3yul06jvf+Y4ef/xxjRw5MpRDAgAAhghZSTlw4IA2btyokpISjRgxQn/4wx/0y1/+UuvWrVNKSkqohtVnRdo7VkL1v+VI+196pOUFYJaQlZTf//73mjp1qqZNmyZJevzxx/U///M/2r17tx555JFQDUtS5P3CBwDARCEpKW63WzU1Nfr+97/v8/jo0aN19OjRdtu7XC65XC7vssViUWxsrKzW4Aw/atiIHu8bbbNd9/NbLBZJks1mk8fj6dY+oR5zT1zPmAEAwdXV74ae/K7y9/d2SEpKU1OTWltblZiY6PN4YmKinE5nu+0rKytVUVHhXb7jjju0cOFCJSUlBWeAv/5dcI7rJ79e9jJkzH7pi2MGAPgI5iUaIX0LclsL6+qxoqIibdy40fvniSee8JlZCTeXLl3SkiVLdOnSpVAPpdeQOTKQOTKQOTL0RuaQzKQMHDhQUVFR7WZNvv7663azK9LVqSRbiF6SCAWPx6MTJ050e/osHJA5MpA5MpA5MvRG5pDMpFitVuXk5Ojw4cM+jx8+fFgjRnCdAgAACOG7ewoLC/Vv//ZvysnJ0fDhw/WHP/xBZ86c0YwZM0I1JAAAYJCQlZSJEyfq3LlzevPNN9XY2KihQ4fqqaeeUmpqaqiGZAybzaaHHnoool7iInNkIHNkIHNk6I3MFk8kvYAGAAD6DD5gEAAAGImSAgAAjERJAQAARqKkAAAAI4Xs3T3o2Pz589XQ0ODz2MyZM/XDH/7Qu3zkyBFt27ZNtbW16t+/v/7hH/5Bs2fPVnR0dG8PNyC6k/mvf/2rtmzZopqaGlksFg0bNkxz5sxRVlZWL482MLrK/Mc//lG/+c1vOtz31Vdf7fCmh6brznmWrmb//e9/r9OnTysuLk4TJkzQvHnzenOoAdOdzMXFxe32Kykp0V133RX08QVDd8+zJJ07d07/+q//qrNnz+o///M/NWDAgN4aZkB1lfncuXP69a9/rdraWp07d06JiYkaO3asZs+erbi4uFAM+bp1lfnkyZN66623dPToUTU1NSktLU0zZszQvffe69fzUFIMVFxcrOnTp3uX+/fv7/36888/1+rVq/XAAw/opz/9qc6ePatXX31Vra2tevTRR0Mx3IDoLPOlS5e0atUq3XbbbSopKdGVK1dUXl6uVatW6ZVXXgnaB00GW2eZJ06cqPz8fJ/t169fL5fL1ScLSpvOMktXPx39nXfe0Y9+9CPdeOONcrlc+r//+7/eHmZAdZVZkkpLS33Od1/9xdWmO5kl6ZVXXtF3v/tdnT17treGFjSdZbZYLLrtttv0gx/8QAMHDpTD4dCGDRt0/vx5LVy4MBTDDYjOMtfU1GjgwIFasGCBkpOTdfToUf32t79VVFSU7rnnnm4/R9/81z3MxcbGym63d7jugw8+0He/+1099NBDkqT09HTNnj1bL7/8smbNmqXY2NheHGngdJb5yy+/1IULF1RcXOz9IKtZs2Zp8eLFOnPmjNLT03txpIHTWeaYmBjFxMR4l5uamvTnP/9ZP/nJT3ppdMHRWebz58/rjTfe0JIlS3TLLbd4Hx86dGgvjS44OsvcJi4urstt+pLuZN69e7cuXryohx56SIcOHeqdgQVRZ5nj4+N9ZsZSU1N111136Z133uml0QVHZ5mnTp3qszx48GAdO3ZMH3/8MSWlr3v77bf15ptvKjk5Wbfffrvuv/9+72yB2+1ud+OcmJgYuVwu1dTUKC8vLxRDvm6dZR4yZIgSEhK0d+9ePfDAA2ptbdXevXs1dOjQPn3zv84yf9u+ffvUr18/TZgwoZdHGVidZT58+LA8Ho/Onj2rn/3sZ7p06ZKGDx+uRx99NKifshps3TnPr7/+uv7jP/5DaWlpmjJliqZPn66oqL57yWBXmU+dOqWKigr98pe/7PMzZW38+Xk+e/asDh48qJEjR/byKAPLn8ySdPHiRcXHx/v1HJQUw3zve99TTk6OBgwY4L0Oo76+Xv/8z/8sSRozZozeffddVVVVaeLEiXI6ndqxY4ckqbGxMZRD77GuMsfGxmrZsmV6/vnn9eabb0q6WlyWLl3aZ6/D6Srzt7333nuaNGmSz+xKX9NV5vr6erW2tqqyslKPP/644uLitG3bNq1cuVK/+tWv+uTLet05zw8//LBuueUWxcTE6MiRIyorK9O5c+f04IMPhnDkPddVZpfLpZdffllz5sxRSkpKWJSU7v48v/TSS/rkk0/U0tKiv/u7v7vmz3tf4O+/YceOHdOHH36op556yq/n4Y6zvaC8vFwVFRWdbrN69WoNGzas3eMfffSRXnzxRW3YsEEJCQmSrr5uv337dl2+fFk2m00PPvigtmzZokWLFmnixIlByeCvQGZuaWnRsmXLNGTIEN1zzz1qbW3VO++8oy+//FKrV6825hd3oM9zm2PHjumZZ57RmjVrlJOTE9AxX69AZt6xY4feeOMNLV26VGPGjJF09WWuJ554Qk899VS7a3RCJVjnuc0777yjiooKbdq0KSDjDYRAZt60aZMaGxu1aNEiSdJf/vIXLV++3LgLZ4Nxnp1Opy5cuKAvv/xSW7du1ahRo1RSUhLwsfdUsL636+rqtHz5cn3ve9/zu3z3vf+a9EH33HOP7rjjjk63udbLFsOHD5ckORwO74kvLCzUfffdp8bGRsXHx6u+vl5btmxRWlpaYAd+HQKZuaqqSg0NDVq5cqV3CnzhwoX68Y9/rD/96U9dPk9vCfR5brNnzx5lZWUZV1CkwGZOSkqSJH3nO9/xbjNw4EANHDhQZ86cCdCIr1+wznOb3NxcXbp0SU6n05jrVAKZ+c9//rNqa2v10UcfSZLa/p88b948PfDAAx2+2ykUgnGe7Xa77Ha7brjhBiUkJOi5557Tgw8+6P3eD7VgZD516pRWrFihadOm9Wh2kJLSC9r+oe2JEydOSFK7b2KLxaJBgwZJunoxbXJyslG/xAKZ+fLly7JYLLJYLN5t2r42aSIwGOe5ublZH374oR555JHrHl8wBDLziBEjJF29UDo5OVnS1Ytpm5qajLr2KBjn+ZtOnjwpm81m1KxCIDP/y7/8i1paWrzrjx8/rldeeUUrVqzQ4MGDr3+wARLs89z2b5fL5erRcwRDoDPX1dVpxYoVmjx5smbPnt2j41JSDHLs2DEdO3ZMN998s+Li4vTXv/5VmzZt0tixY30uHNy5c6fy8/NlsVj08ccf66233tLPfvazPnmhXXcyjx49Wv/1X/+lDRs26J577pHH49Fbb72l6OjoPnmhcHfPsyQdOHBAV65c0aRJk0I02sDoTuYhQ4Zo7Nix2rhxo/7xH/9RsbGx2rJli2644YawPc+ffPKJnE6nhg8frpiYGP3lL3/R1q1bNX369D75abrdyfztd+OdO3dOknTDDTcYVcy6qzuZP/vsM3399dcaNmyY+vfvr1OnTul3v/udRowYYdQMeHd1J3PbSzyjR49WYWGhnE6nJCkqKsqvIkRJMYjVatWHH36oiooKuVwupaamatq0aZo5c6bPdocOHdKOHTvkcrmUlZWlJ598UgUFBSEa9fXpTuYbbrhBS5Ys0fbt2/XMM8/IYrEoOztbTz/9tDHTpP7o7nmWpL1792r8+PF+XxFvmu5m/ulPf6pNmzZpzZo1slgsGjVqlJ5++uk+edFsdzJbrVbt3r1bmzdvlsfjUVpamh5++GHdfffdIRx5z/nzvR0uupM5JiZGe/bs0aZNm+RyuZSSkqJx48bp+9//fugGfh26k/nDDz9UU1OTqqqqVFVV5X08NTVV69ev7/ZzceEsAAAwUt97fQAAAEQESgoAADASJQUAABiJkgIAAIxESQEAAEaipAAAACNRUgAAgJEoKQAAwEiUFAAAYCRKCgAAMBIlBQAAGImSAgAAjPT/AZsejjcZfqKmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print([\"%2.f%%\" % (100 * np.mean(r)) for r in my_alg.X.shared['acc_rates']])\n",
    "plt.hist(my_alg.X.lpost, 30);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You do not need to know much more about class `ThetaParticles` in practice (if you're curious, however, see the next tutorial on SMC samplers or the documentation of module `smc_samplers`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Waste-free versus standard SMC samplers\n",
    "\n",
    "The library now implements by default waste-free SMC ([Dau & Chopin, 2020](https://arxiv.org/abs/2011.02328)), a variant of SMC samplers that keeps all the intermediate Markov steps (rather than \"wasting\" them). In practice, this means that, in the piece of code above:\n",
    "\n",
    "* at each time $t$, $N=20$ particles are resampled, and used as a starting points of the MCMC chains; \n",
    "\n",
    "* the MCMC chains are run for 49 iterations, hence the chain length is 50 (parameter ``len_chain=50``)\n",
    "\n",
    "* and since we keep all the intermediate steps, we get 50*20 = 1000 particles at each iteration. In particular, we do O(1000) operations at each step. (At time 0, we also generate 1000 particles.) \n",
    "\n",
    "Thus, the number of particles is actually `N * len_chain`; given this number of particles, the performance typically does not depend too much on `N` and `len_chain`, provided the latter is \"big enough\" (relative to the mixing of the MCMC kernels).  \n",
    "\n",
    "See Dau & Chopin (2020) for more details on waste-free SMC. If you wish to run a standard SMC sampler instead, you may set `wastefree=False`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_ibis = ssp.IBIS(my_static_model, wastefree=False, len_chain=11)\n",
    "my_alg = particles.SMC(fk=my_ibis, N=100, store_history=True)\n",
    "my_alg.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This runs a standard SMC sampler which tracks $N=100$ particles; these particles are resampled from time to time, and then moved through 10 MCMC steps. (As explained in Dau & Chopin, 2020, you typically get a better performance vs CPU time trade-off with wastefree SMC.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single-run variance estimates\n",
    "\n",
    "An advantage of waste-free SMC is that it gives you the possibility to estimate the (asymptotic) variance of a given estimate from a single run. Here is a quick example, but check also the documentation of `smc_samplers.var_wf` and the collectors `smc_samplers.Var_logLt`, `smc_samplers.Var_phi` for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimate of E(mu]: 2.9956255516306944\n",
      "asymptotic variance of the above estimate: 0.03999515708204\n",
      "95% confidence interval for E[mu] 2.9956255516306944 +-0.012395377987232371\n"
     ]
    }
   ],
   "source": [
    "phi = lambda x : x.theta['mu']  # scalar function\n",
    "est = np.average(phi(my_alg.X), weights=my_alg.W)\n",
    "var_est = ssp.var_wf(my_alg, phi)\n",
    "print(f'estimate of E(mu]: {est}')\n",
    "print(f'asymptotic variance of the above estimate: {var_est}')\n",
    "print(f'95% confidence interval for E[mu] {est} +-{1.96 * np.sqrt(var_est / 1000)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The confidence interval accounts for the *Monte Carlo error* (not the statistical error, which you would measure through the posterior variance). Classically, you must divide the asymptotic variance by the number of particles to get an estimate of the actual variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regarding the MCMC steps\n",
    "\n",
    "The default MCMC kernel used to move the particles is a Gaussian random walk Metropolis kernel, whose covariance matrix is calibrated automatically to $\\gamma$ times of the empirical covariance matrix of the particle sample, where $\\gamma=2.38 / \\sqrt{d}$ (standard choice in the literature). \n",
    "\n",
    "It is possible to specify a different value for $\\gamma$, or more generally other types of MCMC moves; for instance the following uses Metropolis kernels based on independent Gaussian proposals: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=0, ESS=58.32\n",
      "t=1, Metropolis acc. rate (over 9 steps): 0.383, ESS=423.19\n",
      "t=2, ESS=165.63\n",
      "t=3, Metropolis acc. rate (over 9 steps): 0.526, ESS=626.73\n",
      "t=4, ESS=705.49\n",
      "t=5, ESS=493.02\n",
      "t=6, ESS=368.40\n",
      "t=7, ESS=243.54\n",
      "t=8, ESS=196.84\n",
      "t=9, Metropolis acc. rate (over 9 steps): 0.748, ESS=942.05\n",
      "t=10, ESS=801.91\n",
      "t=11, ESS=622.29\n",
      "t=12, ESS=531.48\n",
      "t=13, ESS=426.85\n",
      "t=14, ESS=807.19\n",
      "t=15, ESS=727.41\n",
      "t=16, ESS=655.57\n",
      "t=17, ESS=586.12\n",
      "t=18, ESS=514.90\n",
      "t=19, ESS=439.62\n",
      "t=20, ESS=376.15\n",
      "t=21, ESS=319.90\n",
      "t=22, ESS=322.03\n",
      "t=23, ESS=403.06\n",
      "t=24, ESS=351.35\n",
      "t=25, ESS=486.05\n",
      "t=26, ESS=443.87\n",
      "t=27, ESS=442.08\n",
      "t=28, ESS=445.23\n",
      "t=29, ESS=412.93\n"
     ]
    }
   ],
   "source": [
    "mcmc = ssp.ArrayIndependentMetropolis(scale=1.1)\n",
    "# Independent Gaussian proposal, with mean and variance determined by \n",
    "# the particle sample (variance inflated by factor scale=1.1) \n",
    "alt_move = ssp.MCMCSequenceWF(mcmc=mcmc)\n",
    "# This object represents a particular way to apply several MCMC steps \n",
    "# in a row. WF = WasteFree\n",
    "alt_ibis = ssp.IBIS(my_static_model, move=alt_move)\n",
    "alt_alg = particles.SMC(fk=alt_ibis, N=100,ESSrmin=0.2,\n",
    "                        verbose=True)\n",
    "alt_alg.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the future, the package may also implement other type of MCMC kernels such as MALA. It is also possible to define your own MCMC kernels, as explained in the next tutorial. \n",
    "\n",
    "For now, note the following practical detail: the algorithm resamples whenever the ESS gets below a certain threshold $\\alpha * N$; the default value $\\alpha=0.5$, but here we changed it (to $\\alpha=0.2$) by setting `ESSrmin=0.2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGxCAYAAABiPLw8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiiElEQVR4nO3deXxU9b3/8dd3spGQkBAWEyBBdkVWEXHBuhVtkRaxFpd6xaJd1Frv0lt7a3vbem3Vtle9P69dbrW1xWpVFHGhimKLglbcAUFQFsMqIAlZyD7f3x9nzmQhgSScmXNm5v18PHwwMzlzzpfjkHzy/X6+n4+x1lpEREREUkDI7wGIiIiIxIsCHxEREUkZCnxEREQkZSjwERERkZShwEdERERShgIfERERSRkKfERERCRlKPARERGRlKHAR0RERFJGut8DCJry8nKampo8PeeAAQPYu3evp+dMBbpvPaP71jO6b92ne9Yzum8909l9S09Pp2/fvl0+jwKfdpqammhsbPTsfMaY6HnVHaTrdN96RvetZ3Tfuk/3rGd033rGy/umpS4RERFJGQp8REREJGUo8BEREZGUocBHREREUoYCHxEREUkZCnxEREQkZSjwERERkZShwEdERERShgIfERERSRkKfERERCRlKPARERGRlKHAR0RERFKGAh+RgLHNzX4PQUQkaSnwEQmQ8IuLCV//ZeyGtX4PRUQkKSnwEQkQ+8EaaG7CvvK830MREUlKCnxEgqT2IAD2/bexYS15iYh4TYGPSJDUOYEP1VWw5UN/xyIikoQU+IgESWTGB8CuftPHgYiIJCcFPiJBUlcbfWjXvOHjQEREkpMCH5EgaTXjw7Yt2PJP/RuLiEgSUuAjEhC2sRGaGp0nxwx2Xlv7lo8jEhFJPgp8RIKi1TKXmXoGoDwfERGvKfARCQp3R1dmFmbSyc7j9e85M0EiIuIJBT4iQeHm92TnQMlwyO8L9bXw4fv+jktEJIko8BEJCnfGp1cOJhTCjJsCgF2j5S4REa8o8BEJitpIjk+vbADM+JMA5fmIiHhJgY9IQNi6VktdAGMnQVo67NmJ/WSnb+MSEUkmCnxEgqK2ZakLwGTnwKixgJa7RES8osBHJCgiMz4mOzv6khmvPB8RES8p8BEJimiOT070JTN+qvNg41psqzo/IiLSMwp8RIKifY4PQNFgGFAETU3wwXv+jEtEJIko8BEJinY5PgDGGO3uEhHxkAIfkYBo2dWV3eb1ljyft7DWxntYIiJJRYGPSFDUHZrjA8CY8ZCZBRWfwvatcR+WiEgyUeAjEhS17q6utoGPyciE4ycCYFe/EfdhiYgkEwU+IkFRd2iOjyvavmLtW/EckYhI0lHgIxIU7nb2djk+0NK+gk0bsNWVcRyUiEhyUeAjEhS1Nc6fHc349BsAg4eCDWPffyfOAxMRSR4KfEQCwDY3Q0O98yT70MAHWs36aFu7iEiPKfARCYLWVZk7mPGBVt3a338bG26Ox6hERJKOAh+RIHATm9MzMBkZHR8z4jjIyYWaKti8MX5jExFJIgp8RIKgtoN2Fe2YtDTMCZMBNS0VEekpBT4iQRDdyn7ojq423OUuBT4iIj2iwEckCKJb2Tuf8QEw404EY2DbFmz5p3EYmIhIclHgIxIA9jDFC1szefkwbLTzHs36iIh0mwIfkSDoQo6Pq3XTUhER6R4FPiJBEJnxMUfK8QHM+KnOg/XvYhsbYzkqEZGko8BHJAi6mOMDQOlwyC+E+jr4cG1sxyUikmQU+IgEQRdzfACMMU6SM1ruEhHpLgU+IkFQ28Xt7BFmQmRbu9pXiIh0iwIfkQCI7urqylIXwNhJkJYOe3ZiP9kZs3GJiCQbBT4iQeDm+HRhqQvA9MqB0ScAYNe8EatRiYgkHQU+IkHg7urq6owPYMZpW7uISHcp8BEJgm7m+EBLng8b12Jbd3cXEZFOKfARCYK6bmxndx0zGAYUQVMTrH8vNuMSEUkyCnxEgqAb29ldxhiMmpaKiHSLAh8Rn9lwuGczPtAq8HkLa63XQxMRSToKfET81lAHbtDSjRkfAMaMg8wsqPgUtm3xfmwiIklGgY+I3w5GlrlCIcjM7NZbTUYmHD8R0HKXiEhXKPAR8Vu0eGFvjDHdfrvyfEREuk6Bj4jferCVvTUz3qnnw+aN2OpKjwYlIpKcFPiI+K2Hic0uUzgAhhwLNox9/23vxiUikoQU+Ij4rQdb2dtzZ33UtFRE5PAU+Ij4zNZ2s0FpB8z4qc651r6NbW72YlgiIkkp3e8BNDc389hjj/HKK69QUVFB3759Oeuss7jooosIhZy4zFrLY489xrJly6iurmbUqFFcffXVlJSURM/T2NjIggULWLlyJQ0NDYwbN45rrrmGfv36+fVXE+kat09XD3N8ABg+BnJyoaaKhg1roe9AjwaXGMLPPQ5lmzHzvo3JyvJ7OCISYL7P+CxevJgXXniBq6++mrvuuosrrriCp556iueee67NMc8++yzz58/ntttuo6CggFtvvZXa2pb+RA888ACrVq3ixhtv5JZbbqGuro7bb7+dcDjsx19LpOtqjy7HB8CkpWFOmOyc7o0VXowqYVhrsc88gn3jFezfl/g9HBEJON8Dn40bN3LSSSdx4oknMnDgQE455RQmTJjApk2bAOeb2pIlS5gzZw7Tpk2jtLSU66+/nvr6elascL7BHzx4kJdeeokrr7ySCRMmMGzYMG644QbKyspYvXq1n389kSPzIMcHgEjT0ro3Vh7lgBJMbQ3U1wFgn1uIde+niEgHfF/qOu6443jhhRfYuXMngwYNYuvWrWzYsIF58+YBsGfPHioqKpg4cWL0PRkZGYwdO5YNGzYwY8YMNm/eTHNzMxMmTIgeU1hYSGlpKRs3bmTSpEmHXLexsZHGxsboc2MM2dnZ0cdecc/l5TlTQUrdt8iuLpOdc1R/39C4KTQDjVs2kl5Themd59EAA65if8vj6ipY9gxm1iXdOkVKfd48onvWM7pvPePlffM98Jk9ezYHDx7kX/7lXwiFQoTDYS699FKmT58OQEVFBQD5+flt3pefn8++ffuix6Snp5Obm3vIMe7721u0aBELFy6MPh82bBh33HEHAwYM8Ohv1lZRUVFMzpvsUuG+7SNMLdDnmCLyiot7fqLiYnYNOZam7Vsp2P8J2SNHezbGIKvduZV9AKE0CDdjX1jMMZddTSi3+4FfKnzevKZ71jO6bz3jxX3zPfB59dVXeeWVV/j2t79NSUkJW7du5YEHHogmObvaR3ldach4uGPmzJnDrFmzDjn/3r17aWpq6ubfonPGGIqKiti9e7eaSHZDKt235v3OjEVlQxPVu3Yd3bmOHQXbt7L/jVcJlYz0YniBF978IQDm+AnY8n3YndvY+eBvSZv9lS6fI5U+b17RPesZ3beeOdx9S09P79akhe+Bz4MPPsjs2bM5/fTTASgtLWXv3r08+eSTnHXWWRQUFABEd3y5Kisro7NABQUFNDU1UV1d3WbWp7KykjFjxnR43YyMDDIyMjr8Wiw+jNZafch7IBXuWzQnpVf2Uf9dzYjjsCteIPzRekyS3zeX3e/M/FI4gNAZ5xP+ze3YFxYTPmcWJrdP986VAp83r+me9YzuW894cd98T26ur6+Pblt3hUKh6F9s4MCBFBQUtElSbmpqYt26ddGgZvjw4aSlpbU5pry8nLKyMkaPTo3pfklg0To+R7GdPcKMPN55sHUj1sOZy0Cr+NT5s6AQJp8CJcOgrhb7/CJ/xyUigeT7jM+UKVN44okn6N+/P0OGDGHr1q0888wznH322YAzvTVz5kwWLVpEcXExRUVFLFq0iKysrGgeUE5ODueccw4LFiwgLy+P3NxcFixYQGlpaZuEZ5FAcltWHO2uLoBjBhPKyydcdQC2bYFho47+nAFnyyOBT9/+mFCI0OyvEP7fW7EvPYOd8UVMn76HP4GIpBTfA5/58+fzyCOPcN9993HgwAEKCwuZMWMGF198cfSY2bNn09DQwH333UdNTQ0jR47k5ptvju7CApg3bx5paWncdddd0QKGN9100yGzSSKB40HlZpcJhcg4bjx1b6zAblqPSYHAh0jgYwoixUonTIVho2HLRuxfn8BccrWPgxORoPE98MnOzuaqq67iqquu6vQYYwxz585l7ty5nR6TmZnJ/PnzmT9/fgxGKRIb1lrv6vhEZB0/gbo3VsCmD+CzX/TknIFWEcnx6esEPsYYZ9bn7h9h/74Ee96FmL6q4C4iDk2HiPipoQHc6uIe5PgAZI51al7Zj9YnffKkbah3avdANPABYOwkGDkWmhqxSx7zZWwiEkwKfET85M72GAOZvTw5ZeaoEyAUcpJ+3R1PycotXpiZ6fQqizDGELrQ2c5uX1mK/XSPH6MTkQBS4CPip9qWrezGo3y0UK9eUDocALtpvSfnDCw3sbmg3yG1vsyY8XDcBGhuwj77qA+DE5EgUuAj4ieP83tcZuRY58GmDzw9b9DYipYdXR0JRYoY2pUvYvfsjNewRCTAFPiI+KnVjI+XzIjjACfPJ6mVO0t5nSUvm5HHw7gpEA5jn34kniMTkYBS4CPiJzfwyent6WnNiEghw+1bsG6doGTk5vgUdL5rKzT7cgDs68uxu7bHY1QiEmAKfER81LpdhZdMYX8o7O/sGNv6oafnDhJb3nYre0fMsaNg0jSwYezTD8dpZCISVAp8RPxU68zGGI9zfKBl1iepl7vaFy/sRHTW541XsNu3xHxYIhJcCnxE/FTnXdXmQ7iBTzInOJcfPrnZZYYMw5zktLgJL9asj0gqU+Aj4qfa2OzqAjAjnQRnNn+AdYskJhHb3AyV5c6TvoVHPN588TIwIXj3H9iPP4rx6EQkqBT4iPipzrvO7IcYMgwys+BgDSRjUm9lhZPDFApBn4IjHm6KSzDTPgNAePFDsR2biASWAh8RP9V62Jm9HZOW5jTrJEkLGbo1fPILMaG0Lr3FfOFSJ1Ba82ZyLwGKSKcU+Ij4yMYyx4dW29qT8Yd8F3Z0tWcGDsKcdi4A4cV/jsWoRCTgFPiI+ClGlZtdZmTyJjjb8iPX8OmImXUJpKXD+vewG9bGYGQiEmQKfET8FEluNrHI8QEYPsb585Md2KoDsbmGX45Qtbkzpt9AzBkzAAgvfjDpO9iLSFsKfET8VBe7HB8A0zsXikucJ8k26xPt09W9wAfAzJwL6Rnw4TpY/6634xKRQFPgI+Kn2tjm+EDyLnfZVp3Zu8v07Yc56/MAhJ/8s2Z9RFKIAh8RP8U4xwdoKWSYbBWcIzM+3V3qcpnPf8nZ7r9lI3b1G16OTEQCTIGPiE9sYyM0NTlPYpXjQ0undrZ+iG1qjNl14sla2+WqzZ0xffpizr4AcHZ4adZHJDUo8BHxizvbA543KW3jmEGQ2weaGqFsc+yuE081VdDY4DwuOHLV5s6Y8y+CrGwo20ztq3/zaHAiEmQKfET84ub3ZPXqcgG+njDGQGTWJ2nyfNzE5tw+mIzMHp/G5PXBfPYLABz4071YdwZORJKWAh8Rv8Qjvyci6Tq1H0Vic3vmvDmQl0/T9o+xf19y1OcTkWBT4CPiF7ddRQzze1zRPJ9N65MilyW6o6uHic2tmZzehC78CgDhpx7G1lQd9TlFJLgU+Ij4JY4zPhw70qlWfKAcPt0T++vFWrm7o6tnic3tmennkXHsSDhYjX36L56cU0SCSYGPiE9sHGr4uExmFpQOd66bDMtd0eKFPU9sbs2kpVHwtX8FwP7tWeyubZ6cV0SCR4GPiF9i3KC0vWRqWGqjDUq9mfEB6DXpZMykaRAOE37sD56dV0SCRYGPiF8iOT4mHktdgBnp7uxKhhkfp0Gp8SC5ubXQxV+FtDRY8yZ27duenltEgkGBj8SMra7EhsN+DyO4amucP+M04+NuaWf7x9jWNYQSUXTGx9vAxxQNxpw9C4Dwo/djm5s9Pb+I+E+Bj8SEffcfhP/1SuyTD/o9lOCKJjfHflcXRGZH+g0EG4bNG+NyzViw9XVwMBI0ejzjA2BmXQK5ebBrG/bl5z0/v4j4S4GPeM42NxN+7AGwYfVAOpzodvY4zfjQqp5PIuf5uFvZs7Jjcu9M71zMFy8HwD71Z2xNtefXEBH/KPARz9nXXoI9O50nu7ZhG+r9HVBA2XhuZ3clQ55Pqx1dxpiYXMJ85nNQXALVVdhnH4nJNUTEHwp8xFO2sbFtHZRwGHaU+TegIKuN71IXtNrZtXkDNpyY+Sv2KJuTdoVJSyM092rnei89g929I2bXEpH4UuAjnrKvPA/79zqNI0ef4LxWtsnnUQVUXWRXVxyXuhg81Fkiqj0IOxO0Vk1kxsccRXPSrjDjToRxU6C5mfBCbW8XSRYKfMQztr4eu+QxAMwFc1vaJCRLR3Cv1cZ/qcukpcHw0UAC5/nEoIZPZ0Jz50MoBO+twq5/L+bXE5HYU+AjnrF/e8ZpidBvIGb6DCgZ4by+TYFPh+JcwNAVXe5K0ArOttyp4ROLHV3tmeISzFkzAQg/cl/CLg+KSAsFPuIJe7AG+9wTAJgvXoZJz8BEWiSwfavqoXQkztvZXe5MXMImOEdmfIzHNXw6Y75wKeTkwo6PsSteiMs1RSR2FPiIJ+yLi6GmCoqGYE45y3lxQJHzQ72xAXZv93V8QWObmqChwXkS5xkfho8GY2DvbmxleXyv7YVI1Wavixd2xuT2cYIfwD75Z6xbQ0hEEpICHzlqtroS+8JiAEKzL8eE0gAwoRCUDHOOUZ5PW/W1LY/jPeOTkwuDSp0nmzbE9dpHyzY1gRusxSnwATBnzYSiwVB1IJrHJiKJSYGPHDX73OPODqWSYXDiaW2+ZkqdPB8lOLfjJjZnZGLSM+J++Wghw0TL8zlQDtZCWjrk5sftsiY9ndCX5wNglz2F3bMrbtcWEW8p8JGjYiv2Y//2LAChC69wZnlaK3HyfLSlvR2f8nuiEjXPxy1eWFB46Gct1safBGMnQVMT4ccfiO+1RcQzCnzkqNgljzq5KiOOc34wtGOGRhKct23BWhvn0QWYD+0qWnM7tfPxR9jGRl/G0CPuVvYY1/DpiDHGKWpoQvD2a9gNa+I+BhE5egp8pMfsvk+wLy8FIDT7Kx23DygqgfQMpxP5vk/iPMIA86NdRWsDiiEvH5qa4OOP/BlDD1i3eGEcavh0xAweijnzfCDSvV3b20USjgIf6TH7zF+guQmOm4A5fmKHx5j0dKdaMICWu6JsrT81fFzGGEjEhqVuu4o41PDpjPni5ZDdG8o2Y199ybdxiEjPKPCRHrG7t2Nf/Rvg5PYcjlvPRzu7WvE7x4eW5a6EyvOJ9unyMfDJy8fMmguAXbSgpdmsiCQEBT7SI/aph8GGYcLUltYUnVHgc6haH/p0tdO6gnOi5F/Fo0FpV5hzZsHAYqiswC5Z6OtYRKR7FPhIt9ltW7BvvAI4uT1HYkrcBGcFPlF+5/gADB0B6elQdQD27vZvHN0RzfGJf3JzayY9g9CXvwqAfWExVvlrIglDgY90W3jxnwEwJ01vaUtxOEOGOTthDpRj3aq7qc7nHB8Ak5EJQ0cCiZHnY61t2c7u84wPABOnwZjx0NSIffyPfo9GRLpIgY90i928Ad5bBSbkJHl2gcnKcqregmZ9XD41KG0vukyZCHk+1ZXOLjSA/L7+joXI9vZLrgFjsG+uwG583+8hiUgXKPCRbgk/+SAA5tSzMcVDuvw+JTi3Fd3V5edSFwlWwdmt4dOnwJdq1x0xJcMwZ5wHQPjh/9P2dpEEoMBHusxuWAPr34O09GjTxi4rVQXnNqJLXf7t6gKiFZzZWRb85pvlkWVSH7eyd8RceAXk9IbtW7AvP+/3cETkCBT4SJdYawkvWgCAOeM8TP9juvX+aIKzZnwcdZFdXX7P+OT3hQFFTv+rLRt9HcuRWHfGx8et7B0xefmYSJK/ffLP2OpKn0ckIoejwEe6Zu1bsOkDp6nmBV/u/vvdZqX7PsEerPZ2bIkoAMnNroRZ7oru6ApW4ANgzvy8U6izpgobWQ4WkWBS4CNHZMPhltyesy/A9GCpwfTOhX4DnSfbtng5vMQUhO3srkRpWBqAqs2dMWlphC77BgD25ee1pCsSYAp85Mjeec1ZouqVjfncl3p+HjfP52P9UGhpUupzjg+tGpZu3hjo5Fxb4X/V5sMxY8Zhpp4B1jqJzglSFFIk1SjwkcOy4WbCix8CwHx2NiavT4/PFa35k+Jb2m24GeojgU8QZnwGlTpLbvW1sKPM79F0rtzfBqVdYS6+CjKznGrYry/3ezgi0gEFPnJY9h/LYdc2yMnFzJh9VOcykTyflN/SXlfX8jgIOT6hNBg2Bgh4nk9FcJe6XKZwAGamkwNnFz6gPl4iAaTARzplmxqxTz8MgPnclzA5vY/uhO6Mz67t2Pr6oxxdAnN/GKalQ1Dq0QS8kKGtPdiSEO5zu4ojMedd6OyUO7Af++xjfg9HRNpR4COdsitehH2fOAXjzrng6E+YXwh5+U5z0x1bj/58iapVfo8xxt+xRLR0ag9o6wp3tic7x/cSAEdiMjKdis5E+njt3uHziESkNQU+0im7dBEAZuZcTFavoz6fMcZpjEmKL3cFaUeXa9gYp5/avk+C2U8twDu6OjRhKoybAs1NhB+5T4nOIgGiwEc6ZJsaox27zclneHZedWqnZckmQIGPyc5x6tCAU68pYGx5sHd0tRft45WW7tTAWv2G30MSkQgFPtKxqkj12VAIeud5dlqjLe2tGpT6v5W9NTMyUsjw3dd9HkkHAly8sDOmaDDms18EcGZ9Ght8HpGIgAIf6UxlhfNnXgEm5OHHxE1w3vEx1u20nWKC0qC0PXPaOQDYN15pqZkTFG67ikRZ6oows+Y6uW17d2OXPun3cEQEBT7SmWjgk+/tefsXOVu4mxph93Zvz50o3D5dAdjK3poZNhpGjoXmJuxLz/g9nDaieUcBruHTEdMrx6ntA9glj2H37/V3QCKiwEc6ZqsqnAd9Cjw9rwmFoGSYc41UTXAO6IwPQOi8CwGwy5/DRgK0QIjM+CTSUpfLTDsTRh4PDfXYhQ/4PRyRlKfARzoWmfExHgc+0LpTe4rm+QQ0xweAiVNh4CA4WINd+aLfo2mRaLu6WjHGELrs62CMs4y4Ya3fQxJJael+DwBg//79PPjgg7z77rs0NDRQXFzMtddey/DhkURYa3nsscdYtmwZ1dXVjBo1iquvvpqSkpLoORobG1mwYAErV66koaGBcePGcc0119CvX+J9owwEd6mrj8dLXdDSsytVd3bVBahdRTsmlIaZ8UXsn3+DffEp7NkzncrOPrKNjVB1wHmSYEtdLlM6AvOZ87HLnyP88G8J/fBuTJq/91UkVfk+41NdXc0Pf/hD0tPT+f73v8+dd97JlVdeSU5Oyw+FxYsX8+yzzzJ//nxuu+02CgoKuPXWW6mtbZmKf+CBB1i1ahU33ngjt9xyC3V1ddx+++2Ew2E//lqJz/1BE4sZn0jrCrZtwabi/x93qSv7KCthx4g59VzIzXOKV77zD7+HAwci+T3pGc64EpS58Apnh+SOj7HL/+r3cERSlu+Bz+LFi+nXrx/XXXcdI0eOZODAgYwfP56ioiLAme1ZsmQJc+bMYdq0aZSWlnL99ddTX1/PihUrADh48CAvvfQSV155JRMmTGDYsGHccMMNlJWVsXr1aj//egnLttrV5bmiIc4PsdqDsG+39+cPOBvkpS7AZGVhzpoJQPj5Rf4X32tVwycola57wuT2wVz4FQDs4j9j3ZIRIhJXvi91vfnmm0ycOJE777yTdevWUVhYyHnnncdnP/tZAPbs2UNFRQUTJ06MvicjI4OxY8eyYcMGZsyYwebNm2lubmbChAnRYwoLCyktLWXjxo1MmjTpkOs2NjbS2NgYfW6MITvyg8jLb67uuRLuG7ab45Nf4PnYTUYG4SFDYetHsG0L5pjBhx6TqPetK2rdXV29vb+3Ht230NkX0PzcE7BlI2bTB5hRY70YXo9Ed3QVxC7widfnLXTm52h++XlntvPJBYSu/FZMrxdLSf1vNIZ033rGy/vme+CzZ88eXnjhBS644ALmzJnDRx99xB/+8AcyMjI488wzqaioACA/v22uSX5+Pvv2OTs9KioqSE9PJzc395Bj3Pe3t2jRIhYuXBh9PmzYMO644w4GDBjg3V+uFXcGK1HsqK4iDAwYPpLM4mLPz79/zHhqtn5E7/17KDjM+RPtvnXF7qYGGoF+g4fQKwb3Fjy4b8XF7D9nJjVLF5P58nP0/8y53gysByqbGzgA5AwaQr8Y3S9XPD5v9d/6Pntu+hr2laX0u+gKMkcdH/NrxlIy/huNB923nvHivvke+ITDYUaMGMHll18OOAHItm3bWLp0KWeeeWb0uPZRXlem3w93zJw5c5g1a9Yh59+7dy9NHhbWM8ZQVFTE7t27/V8y6CIbDhM+UA7AvoYmzK5dnl8jPMD5AVa97j1qOzh/It63rmqK5E/tP1jr+b318r7Z6efB0sXU/uPv7Hz3LcwxgzwaZfc0l20BoLZXDrti8FmEOH/e+hVhpp2JfX05n9zzU9JuusPbIqFxksz/RmNJ961nDnff0tPTuzVp4Xvg07dvX4YMGdLmtSFDhvD6607Z/IKCAsCZ1enbt2/0mMrKyugsUEFBAU1NTVRXV7eZ9amsrGTMmDEdXjcjI4OMjIwOvxaLD6O1NmE+5La60umgDtjefSAW43Zr+Xy8iXA43On0ZSLdty6LLHXZXtmxubd4dN+KS2D8SbDmTcIvLCb0lW96M7ju2t9StTnWn4V4fd7Ml65yWoNs+oDwP/5O6NSzY37NWEnKf6NxoPvWM17cN99/zRgzZgw7d+5s89rOnTuj0dvAgQMpKChok6Tc1NTEunXrokHN8OHDSUtLa3NMeXk5ZWVljB49Og5/iyTjJjb3zsOkxyg2Hnys0w286kDLrp0UYK0NZnf2ToRmzAbAvvqiExD7wCZgn64jMX37YS6YC4B9/IGWNiYiEnO+Bz4XXHABH374IU888QS7d+9mxYoVLFu2jPPPPx9wprdmzpzJokWLWLVqFWVlZdx7771kZWUxffp0AHJycjjnnHNYsGABa9asYcuWLdxzzz2Ulpa2SXiWLorW8CmI2SVMVhYUR2b6UqmCc31dyyxPwFpWdOi4Cc7sXEMDdvlz/owhgYsXHo757GynWOSBcuwzj/g9HJGU4ftS18iRI/nOd77DQw89xOOPP87AgQOZN28eZ5xxRvSY2bNn09DQwH333UdNTQ0jR47k5ptvju7CApg3bx5paWncdddd0QKGN910E6EEXDv3m41D4ANOp3a7swxbtgkzYWpMrxUY7myPCUFmlr9j6QJjDOa8Odj778S+9Az2vAsxGZlxu74Nh1tmBJNoxgec3Y2hS68h/P9uwS57GvuZ833LoxJJJb4HPgBTpkxhypQpnX7dGMPcuXOZO3dup8dkZmYyf/585s+fH4shppZIny7jdYPS9kqGwz/+nlo9u9yim9nZCbOd1Zw0HfvEn6B8H/b15ZjpM+J38aoD0NzsBIp9+h75+ARjxp8E46bA2rcIP/Z70r71A7+HJJL0NB0ih6qMXdXm1kyp27MrhQKfBMrvcZn0dMy5zg5I+8Li+CZkRpqT0qcgdvlmPgvNvRrS0uC9Vdh17/g9HJGkp8BHDhWnpS63Zxef7sHWVMf2WkERbVeROIEPgDnjPMjKhp1l8P7b8btwRUvV5mRliocQrZT9yP3Y5mZ/BySS5BT4yCHiluOTkwv9j3GepEqn9uiMTzDbVXTG5OQ6wQ8QXvpk3K5rkzSxuT3zhcucPmQ7y7Av+5RELpIiFPjIodx2FbHO8YGU69Ruozk+iTXjA2A++wUIhWD9e/HLyypPvq3sHTG9czFfdPt4PYStqfJ5RCLJS4GPHCqGndnbMyUplucTmfExCZTj4zL9BmKmnA44uT5xUZ78S10u85nzYfBQqKnCPv0Xv4cjkrQU+Egb1tr45fgAZugI57qpEvgkaI6Py5x3IQD2jZexbkXlGLIpkOPjMmlpTqIzYP/2LHbXNp9HJJKcFPhIW7UHoSnStT6vIPbXc2d8du/A1tfH/np+S9AcH5c5dhSMPgGam7EvPRP7C0aXuvrH/loBYMZOgoknQzhM+NH7/R6OSFJS4CNtuctcWdlOdeUYMwWFzsySDcP2LTG/nu/cHJ8EXOpyhWZcCIB9+XlsXexaLVhrW3Z1JXlyc2uhL8+HtHRY+zZ2zZt+D0ck6Sjwkbaiy1xxSGx2lUaWu1IhwbkusZe6AJgwFY4ZDLU12BUvxu46tTVOiw9IqcDHHDMIc+4XAAg/ej+2qcnnEYkkFwU+0lYc83tcqVTIMNqMMkGXugBMKIRxm5e++FTs6s6UR1pV5OTGZfYxSMwFcyEv31kC/vuzfg9HJKko8JE2bKRdRVzyeyLcwCclEpzdXV2JPOMDmFPPhtw+TvHJt1+LzUXcqs0pkNjcnsnpjbnwCgDs03/BVlX6PCKR5KHAR9pya/jEccYnmuC8Y2vyT+tH6/j09nccR8lkZuFWG7ZLF8WkjUUq7ejqiJn+WRgyDA7WYJ/6s9/DEUkaCnykLR+WuhhQ5AQCTU2Q7Ft4E3xXV2vm7JmQngFbP4QP13l/gRTb0dWeCaURuvRrANjlz2O3b/V3QCJJQoGPtGF9SG42xkDJMOf6yb7clQzJzRGmT4Gz5AWEX3jS+wtEd3QVen/uBGHGjIMTTwPrbG+Pa4NYkSSlwEfaimxnj+tSF60SnJN4Z5e1Nim2s7dmIlvbeW8V9pOdnp472qcrRWd8XKGLr3Jm1ta/B++97vdwRBKeAh9py53xiWNyMxDN87HJ3Ky0sQGaIzlMSTDjA05ncSZMBWuxL3rcxsJd6kqhrewdMQOKorvowo/+HtvY6POIRBKbAh9py486PrS0rqBsCzYcjuu146Z1sb+sXv6Nw2Mht43FymXe7j6qSN1dXe2ZmRdDfl/Yuxv70tN+D0ckoSnwkSjbUA91kaWYOC91UTQEMjKhvhb27o7vteMlusyVjQkl0T+90eOcIpSNDdjlSzw5pW1sgOpIh3IFPpheOZg5VwJgn3kEW1nu84hEElcSffeVo+a2q0hPj/t2a5OW5nSmJokTnKM7upJjmctljGlpXvrSs94sxbj5PZmZkJN79OdLAubUs2HoSKirxT6p7e0iPaXAR1q02spujIn75VsSnJM0zyfBO7MfjplyutNWouoA9l0PEnDLW3p0+fFZDCITChG69BoA7IoXkjsfTiSGuhX41NXVsW/fvkNe37VrF3fffTf/9m//xk9/+lPWrl3r2QAljiojMz7xTmx2uT27Pk72GZ/Er+HTnklPx5x+LgB2xdKjPl9L8cLU3tHVnhk5FjP1DLCW8CP3aXu7SA90K/B56KGH+K//+q82r1VWVvKDH/yA1157jf3797NmzRp+9rOf8dFHH3k6UIm9aN5AvPN7IlpvaU/Gb+g2WrU5+WZ8AMzpn3UerHsXe7R5WpF2FSaFa/h0xnzpKmcJcOP78Parfg9HJOF0K/DZsGEDp59+epvXlixZQnV1NfPmzeMPf/gDv/rVrxgwYABPPfWUpwOVOIi2q4jvjq6owUMhFHJyjSr2+zOGWErSHB+XGVAEx08EwL667OhO5v7/14zPIUy/AZjzLgIg/NgfnERwEemybgU++/btY+jQoW1ee+edd+jXrx8zZzp9ewoLC7ngggvYuHGjd6OU+Kjyd6nLZGZBcQkA9uMkzF+odRuUJt9Sl8uccR4AdsWL2HDPu7Zbt0Fpitfw6Yz53EXOvfl0D/YFj+sniSS5buf45OXlRZ83NDRQVlbG2LFj2xw3ePBgKivVTTjh+NGnqx1TksQJzkk+4wNgJp0CvfOcdhPvv9PzE0X7dCnw6YjJ6oX50jwA7JLHsMk4QyoSI90KfAoLC9m7d2/0+UcffUQ4HGb48OFtjguHw2RlZXkzQokbG4DAh1K3gnMSJjgneY4PgMnIaOnf9cpRJDmrXcURmWlnwvAxUF+Hfeohv4cjkjC6FfiMHj2a5557jvr6egBefPFFACZNmtTmuG3btlFYqKTEhOPm+OT5lONDS4JzUgY+KTDjA2Cmz3AerH4De6D7hfZsczO4ifZ99X2kM8YYQl+eD0SWFndt83lEIomhW4HPRRddxPbt2/n617/ON7/5TVauXMnUqVMZNGhQm+Nef/11Ro0a5elAJQ7cHB8/Z3wiXdr5dA/N7gxUkrDROj7Jm+MDYAYPdWYimpuxr73U/RNUVkA47CS6+/lZTABm5PEw6RSne/sTf/J7OCIJoVuBz+DBg7nllluYOnUqxx57LJdccgn//M//3OaYiooKcnJyOO2007wcp8SYbW6GmkiLAD9zfHJyYUARAI2bkyxBvi65OrMfjjvrY195ofulCdwaPvmFmFCaxyNLPqGLrnSCxHdfx360zu/hiAReenffMGzYML71rW91+vWCggJuuummoxqU+KDqAFgLJgS5eUc+PpZKhsPe3TRs2gDHlPg7Fi9Fd3WlQOAz9QzsI/fDnp3w4ftOP6+uKldz0u4wxUMw02dgX36e8MIHCN10h6pdixyGZy0rDh48yKZNm9i/X7sLEpK7zJWb5/tv2W6eT+OmD3wdh+dSJMcHwPTKxpx8BuDM+nSHLY98D9FW9i4zX7gMMrNg0wfwzj/8Ho5IoHUr8Fm/fj1PPPHEIa8//fTTfO1rX+P73/8+1157Lb/61a+SsvJuUgvCjq4IE2ld0bBpg88j8VgS9+rqSLSmz1srsTXVXX+jW7VZMz5dZgoKMTNmAxB+4k/YpiafRyQSXN0KfJ577jk++KDtb+EffPABDz74ILm5ucycOZOJEyeyfPlyXnihe7/lib8CsZXdFZnxadrxMba+zufBeKgu+bezt3HsKKcad2MDdtXyrr8v2qdLgU93mPMvgtw+8MkO7Ap9/xXpTLcCn82bNzNlypQ2ry1btoxQKMQPf/hD5s2bx/e//31OPvlkXn75ZU8HKjEW3cpe4OswAEx+X8gvBGuxWz70eziesE2N4LYWSIGlLnC2W0dnfV5e2uVZYNuqM7t0ncnOwcy6FAD79MNYN9AWkTa6FfhUVlZSXFzc5rXVq1czfPhwhgwZEn1t+vTpbN++3ZsRSnxUVTh/BmHGBzAnTALAvvOavwPxSusfQknYnb0z5pSzID0Dtm+Bsi5W465Q1eaeMmee7+yKrKxQKwuRThxVcnNFRQUVFRWH1OzJz8+noUGN8xJKkJa6ADPFaYZr334VGw77PBoPuPk9mVmYtNTZom1652FOPBUA24VKztZaVW0+CiY9AzPnSgDs84talrAlMMIrXmD//7sVe7DG76GkrG4FPgMHDuSjjz6KPl+7di3gVHRuraqqqk1PLwm+QOX4AGbsZEx2b+eH4JYkqOeTavk9rURr+qx6+cg5WwerW5YEC1S1uSfMSac7+VX1tdhn/uL3cKQV21BP+KHfUvP8kzT/v59oOdIn3Qp8Tj/9dJ566imWL1/Oe++9x8KFC+nVqxeTJ09uc9yGDRsoKirydKASY5Ht7KaPf+0qWjMZGWRP+wzg7ApKeLWR3+5SJL+njTHjneWX2oNH/n/p1vDJ7YPJyIz92JKQMYbQxVcBYF9+HvvJTn8HJC0+WA0NTssnPlpP+H9vxUZaQEn8dCvw+fznP09JSQm/+tWv+NnPfsaePXv46le/SnarEvxNTU2sWLGCE044wfPBSgwFbMYHIHv6uQDYt15N/PIIboPSFMrvcZlQCHP6Z4Eu1PRRDR9PmDHjYfxJTtuQRQv8Ho5E2PdWAZA1carzvWDDGsK/+hm2Uakh8dStys1ZWVn85Cc/Yf369VRVVTFixAgGDBjQ5pi6ujquuuqqQ5a/JLhsONxSwDAAu7pcvU48BbJ6wf69sPUjGJa4/d9sXWrV8GnPnH4udvFD8NE67K7tmOIhHR5nVbXZM6GLriS89i2njtKWjZhh+p7sJxsOY1e/AUDel66kqeYgzXf/CNa9Q/g3dxC69nuY9AyfR5kaujXjs27dOhoaGjjhhBM45ZRTDgl6AHJzcxk7dizvvvuuV2OUWDtYDc3NzmMfO7O3F8rqhZkwFUiC5a7a1Kna3BFT0A8mnARw+Boz5drR5RUz5FjMqecAEF74QOLPmia6sk1QsR+ysuk1YQpm1FhC3/oBZGTC6jcI/+6XTs9EibluBT4/+clP2mxTD4fDXHbZZWzZsqXNcZ988gm//e1vvRmhxJ4725PTG5MRrN842uzuSuRv3HVun67UW+pyhdwk59decuoadUTFCz1lZl/u/GDduBbWvOn3cFKau8xlxp0YzV8zx00gdN33IT0d3n4N+/u7sGEFP7F21L26wuFwYv9AkkDm97jM+CmQmQl7d8O2zX4Pp+dqU6cze6fGn+QUpqw6AO+90eEh0aUu5fh4whQOwJwzC4Dw43/UD1Uf2Xcjgc+kk9u8bsadSOib/wFpac7Oxz/+b3KU8Agwz5qUSuKKbmUP0DKXy2T1gnGRJZK3XvV5NEchxXN8AExaGua0yNLLik5q+lQ4yc1GNXw8Yz5/MeTkws4y7Ksv+T2clGQ/3esU8TQhzPiTDvm6mTiV0Nf+HUIh7KvLsH/+jSYUYkiBj0BlZKkrgDM+AGbKaQDYN1cm7jeDFM/xcZnpzu4u3n/H+WHQnpKbPWd652Iu+DIAdvFD2j7tA7vame1h5HGY3D4dHmOmnIaZ/y9gDPbl57CP3Je43+8CToGPtPTpCmrgM+Ekp+3Bnp2w42O/h9MjLbu6UjfHB8AMHOTU9bEWu/LFNl+z9XXgVrPVUpenzNkXQOEAqPgU+9LTfg8n5USXuSZOO+xxoWlnYuZ923nPsqexjyspPRa6Hfjs3LmTzZs3R//r6LUdO3Z4PlCJoYD16WrP9MqBcScCCbzcVaccH1e0cenKF9vmnLitKrKyU3pJMBZMRibmwisAsH99HFtd6fOIUoetPQgb1gBgJp58hKMhdPq5mCuuc977/CLsUw/HdHypqFt1fADuvffeQ1675557PBmM+KMlx6fAz2EclplyGvbd151t7bMv93s43Vfr7urSD3Rz4qnYnFynPtO696JBbcuOrkKMMf4NMEmZaWdilz4J27dgn30Mc8nVfg8pNax7B5qb4JjBmKLBXXpL6MzPEW5qxP7ld9hn/kI4PZ3QBXNjPNDU0a3A59prr43VOMRPAV/qAjATTsampcOubdidZZhBpX4PqXvqlOPjMhmZmFPPxi57mvCKpaS5s3lqThpTJhQi9KV5hP/nx9i/P4s9dxam/zF+DyvptSxzHXm2p7XQuV9wgp+FD2CffJBwRiah8y6MwQhTT7cCn7POOitGwxBfBXg7u8vk9Iaxk2DNm9i3X028wMfdzp7iOT4uM/2z2GVPw7ursJUVTtAdmfExak4aOydMhuMnwvr3sE8+iLnm3/weUVKzzc3YSP2k7gY+AKHzLyLc2OAkpT/2e8IZGYTOvsDrYaYcJTenOGttS45PALeztxYtZpiIeT6a8WnDDBkGw0ZDcxP2tb85L0Z3dGnGJ1aMMYS+NA8A+/pybNkmn0eU5DZ9ADVV0DsPRhzXo1OYCy7BzIzsynvot4Rf6aQURIxYa7G1B7F7d2M/3oRtSPxdgd3O8ZEkU18HDZEGeQGe8QGn8JdNS4PtW7G7d3R5vdxvNtzs3GdQ0m4rZvoM7JaN2BVLseddiFWD0rgwQ0diTj4Tu2o54cf/SNq/3OL3kJJWtFrz+JMwaWk9OocxBi68ApoasUufxC64l+b3VmGye0NWltPP0P0vsxdkZTn1z6LPe0WOy3b+NCEnGKuuhOpKJ9G9uhKqq1o9r2rzGs1NLQPq29/pK5bAvd8U+KQ6d5krMwsT8M7hpnceHDfBqQHz9qvR34ICz93RBZrxacWcfAb20fth9w74aH10xkd9umLPXPgV7NsrYd272HXvYMZO9ntISSka+Ezq/jJXa8YYuPir0NiA/dsSeG8Vh9vkHpMN8JlZEApB+T7CP/8e5vJvEors0Ew0CnxSXQLk97RmppyOff8dZ7krUQIfN78nPSNwvdD8ZHrlYE6a7mxrf2VptGqzihfGnhlQhDlrJvbFpwg/9gdCN4/HpOvHgZfs7u3wyQ6nD9cJRx9YGmPg0q87Gz327Yb6eqivdf5sqIP6Oqc4ZZvX6p3Z5vo657nbBDUjE3L7QG4e5PZxiipGHrv/mejjPOjdB5OVha09SPj3d8G7r2P/9L+Et36IufTrCfd9TZ/0VOc2KA14fo/LTDoF++CvoGwTdu9uzIAiv4d0ZGpX0SlzxnlO4PPWCmiMNC5V4BMXZuZcJ79q+1bsC4sxn/+S30NKKtbtRzd6vFOLzAMmFIJxJ9LTYg+2qRGaw5isrJ5dPzuH0LX/gf3rQuziP2Nffh67bQuhb34PU5g4uXlKbk5xNtFmfPL6OJV/cTq2J4Rou4pgLyX6YvgYKC5x8syshbR0yE2MIDzRmbw+mLnzAbBPP4z9ZKfPI0ou9r3XgaNf5vKSSc/ocdATPUcoROiCuYS+/SOnB9yWjYRv/RfshrUejTL2FPikugSo4dOeOTHSuytRdndpxqdTxphoJWcACgqd32olLsyp58DYydDYQPhP6gruFVtVCR99ADg1yJKRGXcioR/cCUOGQdUBwnf+gPCLixOixYa+w6S66Fb2Aj9H0S3mxFPAGNiyEfvpHr+Hc0S2Vu0qDseccrYz0wOgGj5xZYwh9E/XOYmrG9diV8R3q3SysmveBBuGIcMw/Qb4PZyYMQOKCH3v55hpZ0I4jH3kfux9dzp99wJMgU+KS7SlLgDTpy+MOgEA+/ZrPo+mC2ojjTc149Mhk9cHM/kU57Fq+MSd6X8MZk6kj9fCB1oqaEuPud3Yg7TMFSsmKwtz9b9iLv06pKU5ZRJu/y52zy6/h9YpBT4BF17wK8J//nXspg8TMPABp3cX4PTuCrrIUlfQywX4ycy+HI6fiDl7pt9DSUnmnFlOQcnag4Qf+k1CLFcElW1shLXvAD2r1pyIjDGEzp1F6F//y/lZsn0r4Z/+K3bNW34PrUMKfALMVpZjX34O+/e/woH9sblIpbOry/RJrIRSc+KpzoNNHwT/N9RouwrN+HTGFA0h7V//CzN6nN9DSUkmlEZo3g3OkuO7r0Mi/EIRVBvWOFvKCwqhdITfo4krM3ocoR/c5WxaOFhD+J5bCD/zl8DljinwCbL9+1oe79oem2u4OT6JNuNT0A9GHg8kwHKX2lVIAjCDh2I+fzEA4Yd+i62p8nlEiSm6zDVhakom6pu+/Qh952eYsz4P1mIXP0T4Vz/DHqzxe2hRqfd/JZG0msmwu70PfGxjI7gfxgQLfKDVctfbAf/tVNvZJUGYmV92ygtUHcA++nu/h5NwrLUt1ZpTZJmrIyYjg9BXrsVc9W1Iz4D3VhH+6b9hd5T5PTRAgU+g2TYzPtu8v4A725OW5tRjSDBmshP48OE67IFyfwdzGFbb2SVBmIwMZ8nLGOyry7Dr3vF7SIll2xZnpj4zy2mvk+JCp3+W0PfugMIBsGcn4du+E4i8TAU+QVa+N/rQ7t7h/flbVW02pqe1QP1j+g1wEjKtxb4T4OUubWeXBGJGHIc5y0kyD//p3sBvTQ4Sd7aHsZMwmUdXKDBZmKEjnXo/x0902mq4TbF9FKiWFYsWLeLhhx9m5syZXHXVVYAzdfjYY4+xbNkyqqurGTVqFFdffTUlJSXR9zU2NrJgwQJWrlxJQ0MD48aN45prrqFfvwQvfd86aTcWMz4JuqOrNTPldKfD91uvwlkB3RHk7urSjI8kCHPRPzmVhz/dg33yz5hLrvZ7SAlBy1wdM3n5hG78Max5AzPpFL+HE5wZn48++ogXX3yRoUOHtnl98eLFPPvss8yfP5/bbruNgoICbr31VmprWzpeP/DAA6xatYobb7yRW265hbq6Om6//XbCAcsk7642S10V+z1PDkvEGj7tRXd3bViLdWewgkY5PpJgTK8cQldcD4Bd9jR2y0afR+QtW3fQ8+RtW/4pfPwRGIOZcJKn504GJi0tEEEPBCTwqaur45577uEb3/gGvXv3jr5urWXJkiXMmTOHadOmUVpayvXXX099fT0rVqwA4ODBg7z00ktceeWVTJgwgWHDhnHDDTdQVlbG6tWr/foreaN8X9vnn3i83OVuZU+QBqUdMQOKYOhIsGHsO//wezgdq9N2dkk8ZvwUpyKvDRP+4z1Og8skYBvqCd/yz4T/4+ueJtva1ZGmpMPHOEVWJbACsdR13333MXnyZCZMmMATTzwRfX3Pnj1UVFQwceLE6GsZGRmMHTuWDRs2MGPGDDZv3kxzczMTJrQkkhUWFlJaWsrGjRuZNGlSh9dsbGyksbHlH7Ixhuzs7Ohjr7jn6u45bTgMFZGlroHFsGcX7NqOGT7Gs7G1bGXvG7gcn+7ct9CU0wh//BG8/SrmzM/FemjdF13q6h3z+9zTz1uq033rWOjSr9H8/juw42N47gnMFy6Nfi1R75l9ZSns3Q1A+Nc/I+3mOzE5vY/wri6ILHOFJp582HuSqPfNb17eN98Dn5UrV7JlyxZuu+22Q75WUVEBQH5+2xmJ/Px89u3bFz0mPT2d3NzcQ45x39+RRYsWsXDhwujzYcOGcccddzBgQGz6qhQVFXXr+Ob9+9jZ3AyhEL1PPIWa5xbRu6qcguJiz8b0aWM9B4E+Q0rp4+F5vdSV+9b4uQvZ/cSfsOtXM7B3NmkBWrqz4TDbIzM+xww9lrTC+LRk6O7nTRy6b+0UF1Nz7XfZ/4sfEH72UQZ+7kIySoe1OSSR7lm4vo5dz0d+uc7IhE92kvHgvfT/4X8fVc2dcF0tOz5wVhgGfvYCMrrw/TSR7luQeHHffA189u3bxwMPPMDNN99MZmZmp8e1j/C6Uk79SMfMmTOHWbNmHXKNvXv30tTUdMTzd5UxhqKiInbv3t2tMvB2c2RNPb8vtX2dYKz6ow+o3eVd/5PmT5xzVRGixsPzeqFb9y2UASXDYNsWdj//NKHpn43PILvA1h2EyPg/qazC1Md2uaCnn7dUp/vWOTt6Amb8Sdg1b7L7l/9J2k23Y0KhhLxn4aVPEi7/FPoNJO3r/07zL75P3apX2Hnf/xBqNZvV7fO+8xo0NsCAIvZmZGMO8/00Ee9bEBzuvqWnp3dr0sLXwGfz5s0cOHCA733ve9HXwuEw69ev57nnnuPuu+8GnFmdvn1b1kwrKyujs0AFBQU0NTVRXV3dZtansrKSMWM6XxbKyMggIyOjw6/F4sNore1e4LM/spW9b38oGuK8tmu7p2OLJjfn5Qf2H2BX75s58TTsti2E31qJOf3cOIysa+zBSGJzKITNyIwGQTG/bjc/b+LQfeuY+cq12B99CzatJ/y3Zwmd0/JLY6LcM1tfR/ivziy/uWCuk4tzxXXYB/6H8FMPwdARmPE9S0q2777unDeym6urv5wnwn0LGi/um6/JzePHj+eXv/wlP//5z6P/jRgxgunTp/Pzn/+cY445hoKCgjZJyk1NTaxbty4a1AwfPpy0tLQ2x5SXl1NWVsbo0aPj/nfyinUTm/v2cyqpAuzd5W2CYRLs6nKZKac7D9a9iz1Y7e9gWosWL4x9fo9IrJh+AzBfuhIA+8QC7Kd7j/CO4LF//6tTu6z/MZhTzwEgdPq5uK0Vwvf9N3bPzu6fN9yMXf0m4LSpkODzNfDJzs6mtLS0zX9ZWVnk5eVRWlqKMYaZM2eyaNEiVq1aRVlZGffeey9ZWVlMnz4dgJycHM455xwWLFjAmjVr2LJlC/fccw+lpaVtEp4TTiTwMX0HOM3uemVDOOwkOXvAhpuhOrKdMxkCn+IhMKgUmpuw773h93BaaCu7JAlz5ued/nj1tYT//OuEmq2w9XXYSG6PmXUJJr1lscNccg2MOM5pqvmr27pfsHHLh05Ald0bRp3g5bAlRgKxnf1wZs+ezcyZM7nvvvv4j//4D/bv38/NN98c3YEFMG/ePKZOncpdd93FD3/4QzIzM7npppsIJXKDOLeGT2F/Z6YgstzlWbPS6iqwYTAGcvt4c06fRXt3BaAkepTaVUiSMKEQoStvgPR0WPMmdtVyv4fUZfZvzzrByYAizClnt/maSc8g9M2bnF8Ad3yM/eM93UtLeC+yzDV+SpuASoIrcP+XfvzjH7d5boxh7ty5zJ07t9P3ZGZmMn/+fObPnx/j0cWPjc74ONWnTfEQ7NYPsbu24cmCibvM1TsPk5bmxRl9Z6acjn36L/D+O9jag8GolKx2FZJETPEQzAWXYBf/mfDDv6P57M/7PaQjsnW12OcXAZHZng6+35mCfoS++T3C/30z9o1XYNhozIzZXTv/u5E2FVrmShgJPCWS5KI5PpHtz+6Mj1dd2pMovydqUCkUDYamxpZiYj5Tg1JJNuZzF8HgoVBdScX//bffwzki+7dnoboSBg7CTDur0+PMqLGYuU5rDrvwD9gPjlwA1+7Z5bQTSkvDjJvi1ZAlxhT4BJANN0PFfudJobNFz0QSnK1HS122VYPSZGGMwZzoJDnbt1/1eTQRkRwfoxwfSRIm3e3gHuLg3/9K2J3xCCBbe/CIsz2tmbMvcJbCwmHC//eLlt21nZ1/deTvPuoETO/cwx4rwaHAJ4gOVDiJzKEQ5Bc4rxW3zPhYL3qQRWZ8TDLN+NCS58Oat7B1tYc/OB404yNJyAwbjZnxRQDCD/wP1v1FLWDsS89ATRUUDcac/JkjHm+MwfzTdU5dsKoDhH99O7ax827i7jKXmpImFgU+QeT+llFQiAlFfkPpXwRpadBQ37Zre08l41IXON+wBhQ5xcTWvuX3aJTjI0krNOefyBg2GqorCf/+Lm9+IfOQPViDXfokAGbWpV3OZTSZWYSu+z70zoOtH2If/r+Oz19TDR++77xHgU9CUeATRO3ze8DZLTBwkPNk17ajv0aSBj7GmGhNH/tmAHZ31dY4f2ZrqUuSi8nIpN9NP4PMLFj/XnRJKSjsS0/DwWooLsFMnd6t95r+xxD6+nfAhLCvLCX88nOHnn/tW87M/KBSp1myJAwFPgFkIzM6prBdCe6iwc7XPUhwTsYcH1d0W/vat7D19f4Opk4zPpK8MkqOJXTZ1wGwix/Ebt7g84gc9mA19oXFQCS3J9T9natm7GTMnCuc8z30f9hNH7Q94D0tcyUqBT5BtL9V1eZW3ARnT2r5RHN8+h7+uEQ0dCT0Gwj1dfD+274OxdYqx0eSm5k+AzP1DGhuJvy7X2IP1vg9JOyLT8HBGme256TTe3we87kvwYmnQnMT4d/cga0sd87f1Ihd63xvUeCTeBT4BJAtb9Wnq7VIgrPdraWuw3GWu9xihj7v7qpzd3Up8JHkZIzBXHGt88vGvk+wD/7K16rOtqbaCXyA0Bcv69Fsj8sYQ+irNzrlRCo+Jfzbn2ObmuDDdc4ydl4+DBvl1dAlThT4BFFkxscUtg18jEfVm621UFXhPOmTfEtd4DQtBWe76eF2ZcRcdMZHOT6SvExOLqGvfcdpxvvGK9hXX/JtLPbFxc6/u8FDIfJ94GiYXjmErv++03Zm4/tOjR93mWvC1KMKrMQfCnyCyN211X7Gxw18qg5ga6p6fv7aGmhqch4n4YwPAMNGQ0E/J8dm3bv+jUM5PpIizIjjMF+8HAD78G89yUXsLltT1TLb84XLMB61LTJFQwjN/xfnGsuexq54wXl9kpa5EpECn4Cxzc1wwFlHbh/4mF7Z4M4CHc3OLneZKzsHk5HZ8/MEmAmFgtG7Szk+kkLM578EY8ZDfZ2T79PYGNfr26WLnV82hhwLk0/x9Nxm8imYmZHWSfV1kJ4Bx0/y9BoSHwp8gqZiv9M8NC2949mYyKzPUVVwrnR3dHVw/iQSXe56bxW2Kb7fgCGypOgWMNSMj6QAE0ojdPW/Qm4elG3GPvGnuF3bVlVilz0NQOiLl3s229OamX0ZnDDZeTJ2Eiarl+fXkNhT4BM0bg2fgsIO/+FGd3YdzTRykuf3RI08DvL7Ors7utB3x3MNDU6dD1COj6QM07cfoatuBJx8G7vmzbhc176wCOproXQ4TJoWk2uYUBqhb9yE+fJXCV36tZhcQ2JPgU/A2A6KF7bhwYyPTeIdXa2ZUBpm8qmAT7u73NkeYyBTvxlK6jATT8acMwuA8B9i39LCVh3AvvQsEMntMSZm1zLZOYTOm6OihQlMgU/QdLKjy2WKPejSnqR9ujpiTowEPu/8w9mGGk9ufk+v7JhMu4sEmbn4KifXpupAzFta2OcXOXk3Q0eC6urIEei7cdAcacbHDXz2fYJt6GFV4srkrdp8iNHjILeP06hw49r4Xlv5PZLCTEYmoa//O2RmxrSlha2swP4tPrM9khwU+ATMEZe68gogJxeshU929uwaKbLUBWDS0jCR3R1xX+5qNeMjkopMcQnm0ti2tLDPP+E0bz52FEw4yfPzS/JR4BM0R1rqMqZVBeceLndFkptTYakLaGla+s5r2HBz/C7s1vDJ6R2/a4oEjJk+A3PS9Ji0tLAHyrF/XwJEdnJptke6QIFP0LjFCzsJfKB1Bece1vJxZ3ySfDt71JjxzixZ1QH4cH3cLms14yPitLT4p+taWlr8+deetbSwzz3h7J4cNhrGnejJOSX5KfAJENvUCJVu8cJ+nR9YfJStK9wcn1SZ8UlPx0x2trfGtZih+nSJAJGWFtf8m9PSYtXL2NeOvqWFrdiPXf5XQLM90j0KfIKkYr+Tu5OeDrmdJx6bIqeWT0+Wumx9vVPrAlIm8IFWxQzffi2mu0vaUNVmkSgz8viWlhYP/Ra7e8dRnc8+9zg0NsCI41qKCop0gQKfINnfkth82O3P0S3tO7qfs+IWL8zITK0lmOMnOQHIgf2w+YP4XFN9ukTaaNvS4hdHbGlhmxqx1ZXYvbux27dgP1qHXfMW4X/8Hbv8OUCzPdJ96X4PQFq07Og6zDIXQP+BTp+Ypkb4dC90p5BWNL8nP6W+WZiMDMzEk7H/+Dv2rVcxI8fG/qJ16swu0prb0iL8k29D2WbC//NjTJ8CJx+urtb5N1NX2/L4SLW3Ro6F4yfGZ/CSNBT4BEkk8DGdbWWPMKE0OGYQ7PjYSXDuSeCTQstcLjPlNCfweftV7NyrYx/41aqOj0h7TkuLbxO+96ewYQ1dSnPOzISsbGeWOjvH+bN3HqHZX0mpX+DEGwp8gsRd6jrMji6XKS7B7vgYu2s7ZsLULl8ilWr4HGLsZOeb5/59sPVDZydIDFl3qUs5PiJtmEnTCF37PezObdFAxvTKdn5J6BUJcFo9Nmlpfg9ZkogCnwCx7lb2vgOOfHBRD1tXVDk7ukwqVG1ux2RmYSachH3jFexbKzExDnw04yPSOXPiaRjtQBcfKLk5SKJLXUfI8YGWIobdreWTyjM+OMtd4FRx9qqWSKdqnUJtRjk+IiKBocAnSPbvdf7s4lIXALu2d+8HeIoHPoyb4uQL7PsEyjbH9lra1SUiEjgKfALCNjZGl6G6tNR1zCAwBg5Wt7yvK9epSq3ihe2ZrF4wzunnE/NihqrjIyISOAp8gqIikt+TkQm5eUc83GRmOSXgoXsVnCMzPqmY4+OK23KXurOLiASOAp+g2N9Sw6fL2zMjy13dyvOJLnX17fp7koyZcJJTB2nPTqckQAzYxsaWGiTK8RERCQwFPgHRUrzwyPk9LlM02HnQxZ1dtqkJaqqcJym61AWR3lmRhob2rVdjcxF3tgdSq0K2iEjAKfAJii4WL2wjOuPTxaWu6kh+TygEvXO7M7qkE+3dFas8Hze/J6uXU3BSREQCQYFPUHSjeKHLRGv5dHGpq3W7isP1AksBZuJUSEuHXdu6XxKgK5TfIyISSKn90y9AerLUFW1Wun9fS5Xgw4kGPgXdGVpSMjm5MHYSEKPlrlq3arOWuUREgkSBT1C4S13dmfHJ7QPu7qxPdhzxeFvpbmVP3R1drbXe3eU5zfiIiASSAp+g2N+DGR+ASIJzl5ZrqioAMCmc2NyamTQN0tJg+xbsJzs9PbdVDR8RkUBS4BMAtqEeqiudJ92Y8YHWFZyPPOOT8lWb2zG982DMBADs2x7P+tQp8BERCSIFPkHgFi/MzIKcbu62cnt2dSXBuTK1qzZ3JGbLXZEcH6OlLhGRQFHgEwStlrm6XLwwIrqzqwtb2m2rXV3iMJNPAROCjz/C7t3t3Yk14yMiEkgKfALA9mAre5S71LVnF7a5+fDHuu0qNOMTZfLyYfQJANh3XvPuxG6Oj4oXiogEigKfIOhJ8UJX3/7OEllzE+zddfhjU7xBaWfMlNMBj5e7NOMjIhJICnyCoLylT1d3mVAourPrcK0rbDgc3dWlOj5tOctdBjZvwO7f68k5rVvHRzk+IiKBosAnAI5qqQswRV1oXVFTDeGw81g5Pm2YgkIYeTwA9m2PlrvqtNQlIhJECnyCoNzZ1WX6DujZ+4sjMz6Hq+XjJjb3zsOkp/fsOknM8+WuSI6P0VKXiEigKPAJgqNY6oKWWj5292Fq+bjLXMrv6ZCZfKrzYNN6bMX+oz9hnZa6RESCSIGPz2x9PdRUOU96uNRFkVvEcBvW2o6vo63sh2UK+8PwMWAt9p1/HP0JldwsIhJICnz8Vh5Jps3KhuzePTvHwGKnFk1dLXQ2W6Gt7EfUUsxw5dGfrLbG+VM5PiIigaLAx2+R/B4Ku1+80GUyMmBAkfOks51d2sp+ROZEJ/Bh4/stM2Q9YJuaoKHBeaIZHxGRQFHg47Pojq4e5vdEua0rOktw1lLXEZn+x8DQkWDD2HePYrmrvrblsWZ8REQCRYGP3yJLXT0qXtjKkVpXWDUo7RJPdne5VZszMjHpGR6MSkREvKLAx2+tlrqOSnRn1+GXupTjc3hmSmR31wersdWVPTuJaviIiASWAh+f2VYNSo+GKT5Cs1ItdXWJGTgISoZBOIx98NdOxevucqs2K79HRCRwFPj47Wj6dLXmLnUd2I89WNPmS9balsBHMz5HFLrka5CWjn1rJfbxP3b/BNEZHwU+IiJBo8DHb+VH167CZXJ6Q36h86T9clddLTRGdhn16XtU10kFZsw4zFdvBMAuXUT4pWe69X5bqxo+IiJBpcDHR7auFtzZmaOd8YFWO7vaBT5u1easbExW1tFfJwWEpp2JmfNPANi/3Ne9XV7K8RERCSwFPn5yZ3uyczzp6RTd2dV+xie6zKX8nu4wn78Y85nzwYYJ/+6X2M0buvbGSI6P+nSJiASPAh8/lXuT2BzVWS0f5ff0iDEGc/k3YfxJ0NBA+H9vxe7ZdeQ3KsdHRCSwFPj4yLPihRGd1fKxlZGqzXkFnlwnlZi0NEJf/3coHQFVBwj/z0+wVUfY5q4cHxGRwFLg46dI4GMKB3hzvkgtH/btxjY2trwe7dOlpa6eML2yCX37P6HfQNizk/C9t2Ib6jt/gxqUiogElgIfP3m91FVQ6CTUhsPQeknGTW7WUlePmfy+TvCT0xs2fUD4/rs6rfFj3To+WuoSEQkcBT4+sh5tZXcZY1rq+exuyfOJLnUp8DkqZlApoetuhvR0ePtV7MI/dHxgdMZHu7pERIJGgY+f3KUuj3J8oKWCc5st7dGlrgLPrpOqzJhxmKsiNX5eWEx42dOHHhTJ8TGa8RERCRwFPn6KLnV5lOMDLXk+HQQ+alfhjdC0MzEXzQPAPnIf9p12NX7q1LJCRCSoFPj4xNYebPkB6eWMT2Spy7Za6lKOj/fM5y7CnPk5sNap8bPpg5Yv1mo7u4hIUKX7PYBFixaxatUqduzYQWZmJqNHj+aKK65g0KBB0WOstTz22GMsW7aM6upqRo0axdVXX01JSUn0mMbGRhYsWMDKlStpaGhg3LhxXHPNNfTr511Q4Sl3K3tOb4yXFX7dZqW7dzjJt81NLT+IFfh4xhgDl33DKUmw5k3C/3srof/4udPkVDk+IiKB5fuMz7p16zj//PP56U9/yg9+8APC4TC33nordXV10WMWL17Ms88+y/z587ntttsoKCjg1ltvpdbdPQM88MADrFq1ihtvvJFbbrmFuro6br/9dsI96a4dD+V7nT+92tHl6l8EaenQUO8spbmJzenpkN3b22uluGiNn6EjobrSqfFTWdEyk6cZHxGRwPE98Ln55ps566yzKCkp4dhjj+W6665j3759bN68GXBme5YsWcKcOXOYNm0apaWlXH/99dTX17NixQoADh48yEsvvcSVV17JhAkTGDZsGDfccANlZWWsXr3az79ep2z5p84Dr2r4RJj0dBhY7DzZta1Vfk+BM0shnjK9sgnd8MNIjZ9dhP/nxy1fVI6PiEjg+L7U1d7Bg84yQW5uLgB79uyhoqKCiRMnRo/JyMhg7NixbNiwgRkzZrB582aam5uZMGFC9JjCwkJKS0vZuHEjkyZNOuQ6jY2NNLYq8meMITuyNOFlgOCe65Bzlrfs6PI6IDHFQ5y2Fbt3wDHWebFPYgU+nd63ADIFhZgbf0zz7d+FMidgJy0dk5EZ9/En0n0LEt237tM96xndt57x8r4FKvCx1vLHP/6R4447jtLSUgAqKioAyM9vuyMpPz+fffv2RY9JT0+PBkutj3Hf396iRYtYuHBh9PmwYcO44447GDDA2xkYV1FRUZvn++tqqAHyhg4jv7jY02sdGHU8lW+/RvaBT8ksKqYc6DVgIAM8vk48tL9vgVVcTN2P7mTvzddDUyOh3r3b5KnFW8Lct4DRfes+3bOe0X3rGS/uW6ACn/vvv5+ysjJuueWWQ77WPsqz1h7xfIc7Zs6cOcyaNeuQ8+/du5empqauDvmIjDEUFRWxe/fuNuNp3uHsuqpO78XBXV1ofNkN4dwCAGo2b6Q2Jw+A+sxsdnl8nVjq7L4FWr9iQvP/mfDvfkl4wCBf7ndC3rcA0H3rPt2zntF965nD3bf09PRuTVoEJvD5/e9/z1tvvcVPfvKTNjuxCgoKAGdWp2/fvtHXKysro7NABQUFNDU1UV1d3WbWp7KykjFjxnR4vYyMDDIyMjr8Wiw+jNbaNudt3aDU8+sVtzQrtaUjnMd9ChLyH1n7+xZ0ZuoZhEqGQ16+r+NOtPsWFLpv3ad71jO6bz3jxX3zPbnZWsv999/P66+/zn/+538ycODANl8fOHAgBQUFbZKUm5qaWLduXTSoGT58OGlpaW2OKS8vp6ysjNGjR8fnL9IN1lrv+3S1dsxg58+qA9idkXo+2soeN6ZoMKZ37pEPFBGRuPN9xuf+++9nxYoVfPe73yU7Ozuak5OTk0NmppMcOnPmTBYtWkRxcTFFRUUsWrSIrKwspk+fHj32nHPOYcGCBeTl5ZGbm8uCBQsoLS1tk/AcGAdroD6yXT8GgY/ple30/9q/Dzatd15U1WYRERH/A5+lS5cC8OMf/7jN69dddx1nnXUWALNnz6ahoYH77ruPmpoaRo4cyc033xzdhQUwb9480tLSuOuuu6IFDG+66SZCId8ntQ7lzvbk5mGysmJzjaISJ/CJBFjq0yUiIhKAwOfRRx894jHGGObOncvcuXM7PSYzM5P58+czf/58L4cXG27gUxCDZa4IUzwEu+6dlhcU+IiIiPif45OKoonNhbELfIj07Irqo6UuERERBT5+cIsXxjDwMcUlrZ6EILdPzK4lIiKSKBT4+MGd8SmIYQPV4sEtj3PzMKG02F1LREQkQSjw8YF1c3w87tPVRl4B5ES2VCu/R0REBFDg449Ig9KYLnUZ01LIUIGPiIgIoMAn7pzihXudJ31juNQFmEiCs8kriOl1REREEoUCn3irqYKGBudxLKo2t3bCZOfPER237RAREUk1vtfxSTmRZS7y8jEZmTG9VGjqGdixkzC982J6HRERkUShGZ942x/DHl0dUNAjIiLSQoFPnNk45feIiIjIoRT4xNv+2BcvFBERkY4p8Ik3N8enbwxr+IiIiEiHFPjEWbR4oZa6RERE4k6BT7ztd3J8tNQlIiISfwp84sgpXugudSnwERERiTcFPvFUXQlNjc7jWDYoFRERkQ4p8Iknt4ZPnwJMRoa/YxEREUlBCnziqTy+xQtFRESkLQU+cRTd0aXEZhEREV8o8Iknt3ihZnxERER8ocAnnjTjIyIi4isFPnFkleMjIiLiKwU+8aSlLhEREV8p8IkTGw5DRaR4oZa6REREfKHAJ16qDkBTExgD+YV+j0ZERCQlKfCJFze/p09fTHq6v2MRERFJUQp84sTu144uERERvynwiRft6BIREfGdAp84sdEdXWpOKiIi4hcFPvGi4oUiIiK+U+ATJy3FCwf4OxAREZEUpsAnXtylLs34iIiI+EaBTxzY5uaW4oXK8REREfGNAp84CB8oh+ZmMCEVLxQREfGRAp84aNq723lQUIhJS/N3MCIiIilMgU8cNO/7xHmgZS4RERFfKfCJg+Z9ewB1ZRcREfGbAp84aN4bmfHRji4RERFfKfCJg6boUpcCHxERET8p8IkDd8ZHNXxERET8pcAnDpo/1YyPiIhIECjwiTEbbqb5U/XpEhERCQIFPrFWUQ7hZkhLgz4Ffo9GREQkpSnwiTW3OWl+ISak4oUiIiJ+UuATY9Gu7FrmEhER8Z0Cn1hzu7IrsVlERMR3CnxiTDM+IiIiwaHAJ9aaGiE9QzM+IiIiAaDAJ8bSvnItQxatxJz1eb+HIiIikvIU+MSBCYUw6Rl+D0NERCTlKfARERGRlKHAR0RERFKGAh8RERFJGQp8REREJGUo8BEREZGUocBHREREUoYCHxEREUkZCnxEREQkZSjwERERkZShwEdERERShgIfERERSRkKfERERCRlKPARERGRlJHu9wCCJj09NrckVudNdrpvPaP71jO6b92ne9Yzum8909F96+69NNZa69WARERERIJMS10xVltby0033URtba3fQ0koum89o/vWM7pv3ad71jO6bz3j5X1T4BNj1lq2bNmCJta6R/etZ3Tfekb3rft0z3pG961nvLxvCnxEREQkZSjwERERkZShwCfGMjIyuPjii8nIyPB7KAlF961ndN96Rvet+3TPekb3rWe8vG/a1SUiIiIpQzM+IiIikjIU+IiIiEjKUOAjIiIiKUOBj4iIiKQMNQuJseeff56nnnqKiooKhgwZwlVXXcXxxx/v97AC69FHH2XhwoVtXsvPz+d3v/udTyMKpnXr1vHUU0+xZcsWysvL+c53vsPJJ58c/bq1lscee4xly5ZRXV3NqFGjuPrqqykpKfFx1P460j279957Wb58eZv3jBo1ip/+9KfxHmqgLFq0iFWrVrFjxw4yMzMZPXo0V1xxBYMGDYoeo89bW125Z/q8HWrp0qUsXbqUvXv3AjBkyBAuvvhiJk+eDHj3OVPgE0OvvvoqDzzwANdccw1jxozhxRdf5Gc/+xl33XUX/fv393t4gVVSUsIPf/jD6PNQSBOT7dXX13Psscdy9tln89///d+HfH3x4sU8++yzXHfddRQXF/PEE09w6623cvfdd5Odne3DiP13pHsGMGnSJK677rroczWSdALG888/nxEjRtDc3Mxf/vIXbr31Vu6880569eoF6PPWXlfuGejz1l5hYSGXX345RUVFACxfvpyf//zn/PznP6ekpMSzz5l+osTQM888wznnnMO5554bne3p378/S5cu9XtogRYKhSgoKIj+16dPH7+HFDiTJ0/m0ksvZdq0aYd8zVrLkiVLmDNnDtOmTaO0tJTrr7+e+vp6VqxY4cNog+Fw98yVnp7e5rOXm5sbxxEG080338xZZ51FSUkJxx57LNdddx379u1j8+bNgD5vHTnSPXPp89bWSSedxIknnsigQYMYNGgQl112Gb169eLDDz/09HOmwCdGmpqa2Lx5MxMnTmzz+oQJE9iwYYNPo0oMu3fv5hvf+AbXX389d999N5988onfQ0ooe/bsoaKios1nLyMjg7Fjx+qzdwTr1q3jmmuu4cYbb+Q3v/kNBw4c8HtIgXPw4EGA6A9pfd6OrP09c+nz1rlwOMzKlSupr69n9OjRnn7OUnteLYYqKysJh8Pk5+e3eT0/P5+Kigp/BpUARo0axfXXX8+gQYOoqKjgiSee4Ac/+AF33nkneXl5fg8vIbifr44+e/v27fNhRIlh8uTJnHrqqfTv3589e/bwyCOPcMstt3D77berym6EtZY//vGPHHfccZSWlgL6vB1JR/cM9HnrTFlZGTfffDONjY306tWL73znOwwZMiQa3HjxOVPgE2PGmC69Jg43iQ2gtLSU0aNHc8MNN7B8+XJmzZrl48gST/vPmYq0H95pp50WfVxaWsqIESO47rrrePvttw+7PJZK7r//fsrKyrjlllsO+Zo+bx3r7J7p89axQYMG8Ytf/IKamhpef/117r33Xn7yk59Ev+7F50xLXTHSp08fQqHQIbM7Bw4cOCRilc716tWL0tJSdu3a5fdQEkZBQQHAIZ+9yspKffa6oW/fvgwYMECfvYjf//73vPXWW/zoRz+iX79+0df1eetcZ/esI/q8OdLT0ykqKmLEiBFcfvnlHHvssSxZssTTz5kCnxhJT09n+PDhrF69us3rq1evZsyYMT6NKvE0NjayY8cO+vbt6/dQEsbAgQMpKCho89lrampi3bp1+ux1Q1VVFZ9++mnKf/astdx///28/vrr/Od//icDBw5s83V93g51pHvWEX3eOmatpbGx0dPPmZa6YmjWrFncc889DB8+nNGjR/Piiy+yb98+ZsyY4ffQAutPf/oTJ510Ev379+fAgQM8/vjj1NbWcuaZZ/o9tECpq6tj9+7d0ed79uxh69at5Obm0r9/f2bOnMmiRYsoLi6mqKiIRYsWkZWVxfTp030ctb8Od89yc3N59NFHOeWUUygoKGDv3r08/PDD5OXltan1k4ruv/9+VqxYwXe/+12ys7Ojv3Hn5OSQmZmJMUaft3aOdM/q6ur0eevAQw89xOTJk+nXrx91dXWsXLmS999/n5tvvtnTz5m6s8eYW8CwvLyckpIS5s2bx9ixY/0eVmDdfffdrF+/nsrKSvr06cOoUaO49NJLGTJkiN9DC5T333+/zbq368wzz+T666+PFvp68cUXqampYeTIkVx99dVtkitTzeHu2de+9jV+8YtfsGXLFmpqaujbty8nnHACl1xyScrX3Jo7d26Hr1933XWcddZZAPq8tXOke9bQ0KDPWwd+/etfs3btWsrLy8nJyWHo0KHMnj2bCRMmAN59zhT4iIiISMpQjo+IiIikDAU+IiIikjIU+IiIiEjKUOAjIiIiKUOBj4iIiKQMBT4iIiKSMhT4iIiISMpQ4CMiIiIpQ4GPiCS1DRs28Oijj1JTU+P3UEQkABT4iEhS27BhAwsXLlTgIyKAAh8RERFJIerVJSJJ69FHH2XhwoWHvP6jH/2IE044wYcRiYjf0v0egIhIrJx77rlUV1fz3HPP8Z3vfIeCggIAhgwZ4u/ARMQ3CnxEJGn169eP/v37A3DssccycOBAn0ckIn5Tjo+IiIikDAU+IiIikjIU+IiIiEjKUOAjIkktIyMDgIaGBp9HIiJBoMBHRJJaaWkpAEuWLGHjxo1s2rSJ2tpan0clIn5RHR8RSXoPPfQQy5cvp6KiAmut6viIpDAFPiIiIpIytNQlIiIiKUOBj4iIiKQMBT4iIiKSMhT4iIiISMpQ4CMiIiIpQ4GPiIiIpAwFPiIiIpIyFPiIiIhIylDgIyIiIilDgY+IiIikDAU+IiIikjL+P041MEOOpkhEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(alt_alg.summaries.ESSs)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('ESS');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the algorithm waits until the ESS is below 200 to trigger a resample-move step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SMC tempering\n",
    "\n",
    "SMC tempering is a SMC sampler that samples iteratively from the following sequence of distributions:\n",
    "\n",
    "\\begin{equation}\n",
    "\\pi_t(\\theta) \\propto \\pi(\\theta) L(\\theta)^\\gamma_t\n",
    "\\end{equation}\n",
    "\n",
    "with $0=\\gamma_0 < \\ldots < \\gamma_T = 1$. In words, this sequence is a **geometric bridge**, which interpolates between the prior and the posterior. \n",
    "\n",
    "SMC tempering implemented in the same was as IBIS: as a sub-class of `FeynmanKac`, whose `__init__` function takes as argument a `StaticModel` object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=0, ESS=5000.00, tempering exponent=0.000798\n",
      "t=1, Metropolis acc. rate (over 9 steps): 0.270, ESS=5000.00, tempering exponent=0.0111\n",
      "t=2, Metropolis acc. rate (over 9 steps): 0.245, ESS=5000.00, tempering exponent=0.062\n",
      "t=3, Metropolis acc. rate (over 9 steps): 0.253, ESS=5000.00, tempering exponent=0.227\n",
      "t=4, Metropolis acc. rate (over 9 steps): 0.291, ESS=5000.00, tempering exponent=0.788\n",
      "t=5, Metropolis acc. rate (over 9 steps): 0.333, ESS=9573.62, tempering exponent=1\n"
     ]
    }
   ],
   "source": [
    "fk_tempering = ssp.AdaptiveTempering(my_static_model)\n",
    "my_temp_alg = particles.SMC(fk=fk_tempering, N=1000, ESSrmin=1., \n",
    "                            verbose=True)\n",
    "my_temp_alg.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: Recall that `SMC` resamples every time the ESS drops below value N times option `ESSrmin`; here we set it to to 1, since we want to resample at every time. This makes sense: Adaptive SMC chooses adaptively the successive values of $\\gamma_t$ so that the ESS equals a certain value ($N/2$ by default). \n",
    "\n",
    "We have not saved the intermediate results this time (option `store_history` was not set) since they are not particularly interesting. Let's look at the final results: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGxCAYAAABlfmIpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOvklEQVR4nO3de3gUVZ4//nf1TRJi0p0LhJAwSZCOXIzAjBHx63JTcCWKYXUY1J1fVMbZSXRkZvDC+DgY8MHJuBpcgZ11RQk4ERSEBYImYkDhQcPqPBKWOOCERAUNSSfdIVzT1V2/PyA13enupLvT93q/nidPuqpOVZ9TXX36U6dOnRIkSZJAREREpDCqcGeAiIiIKBwYBBEREZEiMQgiIiIiRWIQRERERIrEIIiIiIgUiUEQERERKRKDICIiIlIkBkFERESkSAyCiIiISJE04c5ApDObzRBF0au0aWlpaG9vD3KOwoNli07RXjaNRgODwRDubPjFl7oj1KL9uAg07g9n0b4/fKk3GAQNQBRFWK3WAdMJgiCnj7UnkbBs0SmWyxYNvK07Qo3HhTPuD2dK2x+8HEZERESKxCCIiIiIFIlBEBERESkSgyAiIiJSJAZBREREpEi8O4yIYtK2bdtw6NAhnDp1CjqdDkajEQ888AAyMjI8rnP06FGUlZW5zK+oqMDIkSODmV0iCgMGQUQUkxobGzFnzhyMHj0aNpsNmzZtwvPPP4+XX34ZQ4YM6XfdVatWIT4+Xp5OTEwMdnaJKAwYBBFRTHrmmWecpktKSrBo0SKcOHEC48aN63fdpKQkDB06NJjZI6IIwCCIiBTh/PnzAICEhIQB0z755JOwWq3IzMzE/PnzMWHCBI9prVar06CIgiAgLi5Ofh1pevMUiXkLB+4PZ0rbHwyCiCjmSZKEyspKXHvttRg1apTHdAaDAY888ghyc3MhiiI++eQTrFixAsuWLfPYerRt2zZs2bJFns7JyUF5eTnS0tICXo5ASk9PD3cWIgr3hzOl7A8GQUQU89atW4dvv/0Wy5cv7zddRkaGU8dpo9EIk8mEnTt3egyCioqKUFhYKE/3nkG3t7dH5LPDBEFAeno6WltbFfFYhIFwfziLhf2h0Wi8PglhEEREMe2NN97AF198gbKyMqSkpPi8vtFoxP79+z0u12q10Gq1bpdF8o+IJEkRnb9Q4/5wppT9wXGCiCgmSZKEdevWob6+Hn/4wx8wbNgwv7bT3NwMvV4f2MwRUURgSxARxaR169bhwIEDePLJJxEXFweLxQIAiI+Ph06nAwBUVVWhs7MTjz76KACguroaaWlpyMrKgiiK2L9/P+rr6/G73/0uXMUgoiBiEETkgePdEUpoFo41tbW1AIDnnnvOaX5JSQmmT58OADCbzTCZTPIyURSxceNGdHZ2QqfTISsrC08//TQmT54cqmxHPH4vKJYwCCJyQxAEoHoz7B3tUKWkQZi7gBV+lHnnnXcGTFNaWuo0PW/ePMybNy9YWYp6giBAem8D7JYOqPQpEOb/nN8LimoMgog8sHe0w376ewDsPEfUy27pgNTRDjsAdbgzQzRIrNuJiIhIkRgEERERkSIxCCIiIiJFYhBEREREisQgiIiIiBSJQRAREREpEoMgIiIiUiQGQURERKRIDIKIiIhIkRgEERERkSIxCCIiIiJFCvuzwxobG7Fjxw40NzfDbDZjyZIlKCgo8Jh+zZo1+Pjjj13mZ2Zm4uWXXwYA7Nu3D2vXrnVJ89Zbb0Gn0wUu80RERBS1wh4EXbp0CdnZ2ZgxYwZeeumlAdM/+OCDuP/+++Vpm82GJ554AlOmTHFKFxcXh1deecVpHgMgIiIi6hX2IGjSpEmYNGmS1+nj4+MRHx8vTx86dAjnzp3DjBkznNIJggC9Xu/1dq1WK6xWq9P6cXFx8uuB9KbxJm20UWLZBEGAAPzjLwrLHsufGxFRIIQ9CBqsuro6XHfddUhLS3Oaf/HiRZSUlMButyM7OxsLFixATk6Ox+1s27YNW7ZskadzcnJQXl7ust2BpKen+1aAKKKksomiiE6NFtDpAI0Wyamp0Gii8+sSy58bEdFgRGetfoXZbMaXX36JX//6107zMzIyUFJSglGjRuHChQvYvXs3nn32Wbz44osYMWKE220VFRWhsLBQnu49e25vb4coigPmRRAEpKeno7W1FZIkDaJUkUeJZRMEAXbRCntPD1SiFSaTKerKHgufm0aj8flEhCKXY6tktB6TFFuiOgjat28fhg4d6tKR2mg0wmg0ytN5eXl46qmn8P777+Ohhx5yuy2tVgutVut2mS9fVkmSYvbLrbSySY5/UVz2aM47xQ5BECC9twF2SwdU+hQI83/O45LCLmqDIEmSsHfvXtxyyy0DXqZQqVQYPXo0WltbQ5Q7IiLqy27pgNTRDjsAdbgzQ4QoHieosbERra2tmDlz5oBpJUnCN99841NHaSIiIoptYW8JunjxolMLTVtbG1paWpCQkIDU1FRUVVWhs7MTjz76qNN6dXV1GDNmDEaNGuWyzXfffRdjxozBiBEj5D5BLS0tePjhh4NeHiIiIooOYQ+CmpqaUFZWJk9v2LABADBt2jSUlpbCbDbDZDI5rXP+/HnU19ejuLjY7TbPnTuH1157DRaLBfHx8cjJyUFZWRmuueaaoJWDiIiIokvYg6Dx48fjnXfe8bi8tLTUZV58fDzeeustj+sUFxd7DJCIiCjweOcXRaOwB0FERBTdeOcXRSsGQURENGi884uiUdTeHUZEREQ0GAyCiIiISJEYBBEREZEiMQgiIiIiRWIQRERERIrEIIiIiIgUiUEQERERKRKDICIiIlIkBkFERESkSAyCiIiISJEYBBEREZEiMQgiIiIiRWIQRERERIrEIIiIiIgUSRPuDBARUewSBMHpP1EkYRBERERBIQgCpPc2wG7pgDorFwIAKdyZInLAy2FERBQ0dksHpI522Lst4c4KkQsGQURERKRIvBxGiuLYL0GS2DBPRKRkDIJIMQRBAKo3w97RDlVKGoS5C8KdJSJlEgSekFBEYBBEimLvaIf99PcAeC2YKFyEJAOkrZWwWTqg0qdAmP9zBkIUFgyCiIgo5OQO0wDU4c4MKRaDICIiChyHS10cG4giHYMgIiIKGMdLXV6NDcT+QRRGDIKIiCig5Etd+uQB07J/EIUTgyCiKxzPRtmMTxQ67B9E4cIgiAhXhvfftQn2jnYAgDrXCIZBRESxjUEQ0RWOt88LyWlhzg0REQUbh0ohIiIiRWIQRERERIrEy2EUkXjLLBERBRuDIIo47p7xxUCIiIgCjUEQRSQ+44uIiIKNvy9ERESkSAyCiIiISJHCfjmssbERO3bsQHNzM8xmM5YsWYKCggKP6Y8ePYqysjKX+RUVFRg5cqQ8/dlnn2Hz5s04ffo0hg8fjoULF/a7XSIiIlKWsAdBly5dQnZ2NmbMmIGXXnrJ6/VWrVqF+Ph4eToxMVF+ffz4caxatQoLFixAQUEBDh06hIqKCixfvhxjxowJaP6JKDJt27YNhw4dwqlTp6DT6WA0GvHAAw8gIyOj3/UaGxtRWVmJkydPwmAw4K677sLs2bNDlGsiCqWwB0GTJk3CpEmTfF4vKSkJQ4cOdbusuroa+fn5KCoqAgAUFRWhsbER1dXVWLx4sdt1rFYrrFarPC0IAuLi4uTXA+lNE4vPnAp12QRBgAD84y9A79vfdh2X4cp/KQh5CKVYPia90djYiDlz5mD06NGw2WzYtGkTnn/+ebz88ssYMmSI23Xa2trwwgsvYNasWXjsscdw7NgxvP7660hMTMSUKVNCXAIiCrawB0H+evLJJ2G1WpGZmYn58+djwoQJ8rLjx49j7ty5Tumvv/567N692+P2tm3bhi1btsjTOTk5KC8vR1qab49PSE9P9yl9NAlV2URRRKdGC+h0gEaL5NRUaDSDP1T7225KSso/lgGARnP5L8B5CIdYPib788wzzzhNl5SUYNGiRThx4gTGjRvndp3a2lqkpqaiuLgYAJCZmYmmpibs3LnTYxA02BOoUBtMcNx7sgA4nxz0ne94AuHr61DvM6WfLPSltP0RdbW6wWDAI488gtzcXIiiiE8++QQrVqzAsmXL5IrNYrFAr9c7rafX62GxWDxut6ioCIWFhfJ07wHQ3t4OURQHzJcgCEhPT0dra2vMjWkT6rIJggC7aIW9pwcq0QqTyRSQ93W3XeBykNDR0QHblWUAoBZFSKIY8DyEUiwckxqNxucTEU/Onz8PAEhISPCY5uuvv0Z+fr7TvIkTJ2Lv3r0QRdFtIByoE6hQ8yc4ttls6NBeOVnQapGSlga1Wu08X6OBpNVA8Oe1wzZDTaknC54oZX9EXRCUkZHhdE3faDTCZDJh586dHs/ugMujDvcX2Wq1Wmi1Wo/rekuSpKj9wRlIKMsmOf4F8H37bleeL0nyfPT5H+g8hFo05z1QJElCZWUlrr32WowaNcpjOovFgqSkJKd5SUlJsNls6O7uhsFgcFlnsCdQoTaY4FgQBNisVkg9PRCsVrS3t8t1a+98lShCsop+vXbcZqjEwslCIMXC/vDl5CnqgiB3jEYj9u/fL0+7a/Xp6upyqdyISBnWrVuHb7/9FsuXLx8wbd+Tpd4fAk8nUYE6gQo1f4Njp5MFh21IfZf58Xow+Rosniw4U8r+iIlxgpqbm50ufxmNRhw5csQpTUNDA4xGY4hzRkTh9sYbb+CLL77AsmXLkJKS0m9adydQZ86cgVqt7vcyGhFFp7AHQRcvXkRLSwtaWloAXL47o6WlRe6vUVVVhdWrV8vpq6urcejQIfzwww/47rvvUFVVhfr6etx+++1ymjvuuAOHDx/G9u3bcerUKWzfvh1Hjhxx6SxNsU8QBKc/Ug5JkrBu3TrU19fjD3/4A4YNGzbgOmPGjEFDQ4PTvMOHDyM3NzdqO8YTkWdh/1Y3NTU5DX64YcMGAMC0adNQWloKs9ksB0TA5Tt8Nm7ciM7OTuh0OmRlZeHpp5/G5MmT5TR5eXlYvHgxNm3ahM2bNyM9PR2LFy/mGEEK4/ggVgBQ5xrBMEg51q1bhwMHDuDJJ59EXFyc3MITHx8P3ZW7AKuqqtDZ2YlHH30UADB79mzU1NSgsrISs2bNwvHjx1FXV4fHH388XMUgoiAKexA0fvx4vPPOOx6Xl5aWOk3PmzcP8+bNG3C7U6ZM4bge5PQgViE5su/WocCqra0FADz33HNO80tKSjB9+nQAcDnJGjZsGJYuXYrKykrU1NTAYDDgwQcfZF1CFKPCHgQREQVDfydXvfqeZAHAuHHjUF5eHowsEVGECXufICIiIqJwYEsQERH5zuFmA950QNGKQRDRQASVx7FjiJRKSDJA2loJm6UD6qxc+fEXRNGEQRDRAARDCqRdm+S7zFQpaRDmLmAgRIpnt3RA6miHXZ8c7qwQ+YVBEJEXJIe7zAB2piMiigWsy4mIiEiRGAQRERGRIjEIIiIiIkViEERERESKxCCIiIiIFIlBEBERESkSgyAiIiJSJAZBREREpEgMgoiIiEiRGAQRERGRIvGxGUQBxoetEvnJ4cn0AL87FHwMgogCSBAEoHozH7ZK5AfHJ9Or9CkQ5v+c3x0KKgZBRAFm58NWifwmP5kegDrcmaGYx/qZiIiIFIlBEBERESkSgyAiIiJSJPYJIiKiiMe7xigYGAQRBUBvBd339ngiGjxBECC9twF23jVGAcYgiJRJUMkBiyiKgwpeHG+LV+cawTCIKPB41xgFA4MgUiTBkAJp1yZIHe3o1GghjcoZVPDSe1u8kJwWsDwSEVFwMQiiqBOoEZml3vF8dDpIifoA5IyIiKIJgyCKKhyRmYiIAoVBEEUdjshMpGB8vhgFEIMgIiKKGny+GAUSgyAiIooqvFOMAoVXEoiIiEiRGAQRERGRIjEIIiIiIkViEERERESKxCCIiIiIFIlBEBERESlS2G+Rb2xsxI4dO9Dc3Ayz2YwlS5agoKDAY/r6+nrU1taipaUFoigiMzMT9957LyZOnCin2bdvH9auXeuy7ltvvQWdTheMYhAREVGUCXsQdOnSJWRnZ2PGjBl46aWXBkz/1VdfIT8/HwsXLsTQoUOxd+9elJeXY+XKlcjJyZHTxcXF4ZVXXnFalwEQERER9Qp7EDRp0iRMmjTJ6/TFxcVO0/fddx8+//xzfPHFF05BkCAI0Ov1Xm/XarXCarU6rR8XFye/HkhvGm/SRptQl00QBAjAP/4c3tdxGdws97Sd3rSSw3bdjTHbX1pP79k3v/2lDaVYPiaJiAIh7EHQYNntdly4cAEJCQlO8y9evIiSkhLY7XZkZ2djwYIFTkFSX9u2bcOWLVvk6ZycHJSXlyMtLc2n/KSnp/tWgCgSqrKJoohOjRbQ6QCNFsmpqdBoNK7LAJflHrcDABrN5T+dzvk1AHWf6f7S9punAdKGQywfkxQ8fD4XKUHUB0G7du3CpUuXcNNNN8nzMjIyUFJSglGjRuHChQvYvXs3nn32Wbz44osYMWKE2+0UFRWhsLBQnu6tANrb2yGK4oD5EAQB6enpaG1tjbkKI9RlEwQBdtEKe08PVKIVJpNJfl/HZQBclnvaDgCoRRGSKMLe0+P0WqfTweYw3V9ad+/p+D4DpQ2lWDgmNRqNzyciNHiCIEB6bwPsfD4XxbioDoIOHDiAd999F0888QSSkpLk+UajEUajUZ7Oy8vDU089hffffx8PPfSQ221ptVpotVq3y3z58kuSFLOVRSjLJjn+9Xnf3vnwsNzddtB3nX7e11Pa/t5T8iFtqIX7/Sny2Gw2+UTP07HB53OREkRtEHTw4EH8+c9/xm9/+1vk5+f3m1alUmH06NFobW0NUe6IiCKTIAjoXLcKtvbTENjKQwoXleMEHThwAGvWrMGvf/1rTJ48ecD0kiThm2++8amjNBFRrJIsnZdbeSwd4c4KUViFvSXo4sWLTi00bW1taGlpQUJCAlJTU1FVVYXOzk48+uijAP4RABUXF8NoNMJisQC4fPt7fHw8AODdd9/FmDFjMGLECLlPUEtLCx5++OGQl4+IiIgiU9iDoKamJpSVlcnTGzZsAABMmzYNpaWlMJvNMJlM8vI9e/bAZrNh3bp1WLdunTy/Nz0AnDt3Dq+99hosFgvi4+ORk5ODsrIyXHPNNSEqFRERDYogcJgHCrqwB0Hjx4/HO++843F5b2DT67nnnhtwm8XFxS7jCRERUfQQkgyQtlbCZumAOivX49heRIMRlX2CiIgo9sl3qHVbwp0VilFhbwkiijqCymXEaCIKA4dLZgAHdSTfMQgi8pFgSIG0axPsHe0AAHWuEQyDiELP8ZIZB3UkfzAIIvKD1NEO++nvAQBCMkc0JgoXDupIg8E+QURERKRIDIKIiIhIkRgEERERkSKxTxAREXnGQQsphjEIIiIijzhoIcUyBkEU2TgmD1HYyXdg6ZPDnRWigGIQRBGNY/IQEVGwMAiiiMcxeYiIKBh4dxgREREpEluCKLr16TPEIfPJUWNjI3bs2IHm5maYzWYsWbIEBQUFHtMfPXoUZWVlLvMrKiowcuTIYGaViMKAQRBFNcc+Q6qUNAhzFzAQItmlS5eQnZ2NGTNm4KWXXvJ6vVWrViE+Pl6eTkxMDEb2iCjMGARR1HPsM8Tru+Ro0qRJmDRpks/rJSUlYejQoUHIERFFEgZBFDt4Oz0FyJNPPgmr1YrMzEzMnz8fEyZM8JjWarXCarXK04IgIC4uTn4daZy+I3CfR0EQ5Lswe8cFEqLgtT/7mwNBOlPa/mAQRDGDt9PTYBkMBjzyyCPIzc2FKIr45JNPsGLFCixbtgzjxo1zu862bduwZcsWeTonJwfl5eVIS4vMOxltNhs6AOh0OkCrRUpaGtRqtWsarRbQ6QCNBpJWAyHSX3soi7fS09MDsHdjh1L2B4Mgiim8nZ4GIyMjAxkZGfK00WiEyWTCzp07PQZBRUVFKCwslKd7z6Db29shimJwM+wHlUoFAUBPTw8EqxXt7e0u/egEQYDNaoXU0wOVKEKyihH/2lNZBiIIAtLT09Ha2sr+hIiN/aHRaLw+CfE7CBJFERoNYygiCqxIq1uMRiP279/vcblWq4VWq3W7LBJ/RCRJkltIpSvT7vIpOfwhCl73ls3ffT6YdWORUvaH3/1If/nLX6KqqgomkymQ+SEihVu2bBl27twJs9kc7qwAAJqbm6HX68OdDSIKAr9Pt3784x/j/fffx86dOzF58mTcfvvtuO666wKZNyJSoPHjx2P//v3Yu3cvrrvuOtx1111+1y0XL15Ea2urPN3W1oaWlhYkJCQgNTUVVVVV6OzsxKOPPgoAqK6uRlpaGrKysiCKIvbv34/6+nr87ne/C0jZiCiy+B0ElZSU4Oc//zk++ugjfPjhh3j++eeRkZGBOXPmYNq0afLdEUREvrjvvvtw991349NPP8XBgwcHVbc0NTU5DX64YcMGAMC0adNQWloKs9ns1JotiiI2btyIzs5O6HQ6ZGVl4emnn8bkyZMDV0AiihiDuvCekJCAefPm4a677sIXX3yBDz74AG+++Sbefvtt/NM//RNuv/12jrJK1AdHuB5YfHw8Zs2ahTlz5uDbb7/1u24ZP3483nnnHY/LS0tLnabnzZuHefPmDTr/RBQdAtL7UBAE/OQnP0FKSgo2bNiAxsZG1NbWora2FjfccAN+8YtfICkpKRBvRRTVBEEAqjdzhGsvsW4homAadBBks9nw6aefoqamBsePH0dqairuv/9+TJ06Ff/7v/+LrVu3YvXq1XjmmWcCkV+iqGfnCNdesdlsOHz4MA4ePMi6hYiCwu8gqLOzEx9++CE++ugjdHV14dprr8VvfvMbFBQUQKW6XLX/8z//M5KTk/Hqq68GLMNEMaPPCNcAL48BgMViwcGDB/HZZ5+hu7ubdQsRBY3fQVBpaSlUKhVuvvlm3HHHHcjOznabbvjw4WyuJnKj7wjXvDx22fLly6FWqzFp0iTMnDnTY6dk1i1ENFh+B0H33HMPbrvttgGfrpydnY01a9b4+zZEMc1xhGuAl8cA4Pbbb8fUqVORkJDgcRBCgHULEQ2e33VuamqqxwesnT17Fh9//LHfmSIi5dLr9axbiCgk/A6C1q5di9OnT7td1tbWhrVr1/qdKSJSrrffftvjSPSsW4gokILS+t7T0yN3YCQiChTWLUQUSD71CTKZTGhra5Onm5ub0dPT45Smp6cHe/bsQWpqamBySEQxz2w2o7OzU54+efIkRFGEWq1Ge/vljuOsW4go0HwKgvbu3YstW7bI06+//rrHtMXFxX5nioiUpb6+HjU1NfK0Yz3TF+sWIgoUn4Kgm266CVlZWQCAiooKLFy4EOnp6U5ptFotsrKyMGzYsMDlkohi2sSJEzFixAhIkoTKykrMnTsXaWlpUKvV8h2orFuIKNB8CoIyMzORmZkJAPjVr36FH//4x7j66quDkjEiUo709HT5hKqnpwfjx4/H0KFDodVqkZaWFubcEVGs8nucoOnTpwckA42NjdixYweam5thNpuxZMkSFBQUDLhOZWUlTp48CYPBgLvuuguzZ892SvPZZ59h8+bNOH36NIYPH46FCxcOuF0iCj9+T4koVHwKgrZs2YKZM2ciOTm532v2ve65554B01y6dAnZ2dmYMWMGXnrppQHTt7W14YUXXsCsWbPw2GOP4dixY3j99deRmJiIKVOmAACOHz+OVatWYcGCBSgoKMChQ4dQUVGB5cuXY8yYMQMXlIhC6oMPPsBNN92EpKQkfPDBB/J8tVqNoUOHuqT3pm4hLwiC05hMSh+tnJTHpyDo3XffxcSJE5GcnIx33313wPTeVFSTJk3CpEmTvM5DbW0tUlNT5c6RmZmZaGpqws6dO+UgqLq6Gvn5+SgqKgIAFBUVobGxEdXV1Vi8eLHX70VEoVFTU4OxY8ciKSnJqYO0JwyCAkNIMkDaWgmbpQMqfQqE+T9nIESK4lMQtHnzZrevQ+nrr79Gfn6+07yJEydi7969EEURGo0Gx48fx9y5c53SXH/99di9e7fH7VqtVlitVnlaEATExcXJrwfSm8abtNEm1GUTBAECIP9JV/5jgGl/0rq8d4C267jMcb8NVLZA7uNoOiZXrVrl9rVGo2GfoCCzWzouP74FgDrcmSEKMb/7BIWLxWJxeWhiUlISbDYburu7YTAYYLFYoNfrndLo9XpYLBaP2922bZvTJb6cnByUl5f7XAH3vVsuloSqbKIoolOjBXQ6QKO5/KfTXV7Y3/Qg0qqDtF1otEhOTYVGo/GibM5pAyWWj0kiT3iZj7zhd23b09MDURQRHx8vzzt48CCam5tx3XXXubTWBFLfM9veA7y/M15JkvpdXlRUhMLCQpf3aG9vhyiKXuUpPT0dra2tMfeFC3XZBEGAXbTC3tMDtShCEkXYrwzK2d+0v2l1Oh1sQdguAKhEK0wmk9Mx6qlsfdMGYj9G4zFptVphs9kwZMgQuSUoVHWLojn0D4qG1sP+CIIA6b0NsPMyHw3A7yBo9erVuOqqq1BaWgoA2L17NyorKwEAO3bswFNPPYXJkycHJpcO3LXonDlzBmq1GgkJCR7TdHV1ubQgOdJqtR6fWO3Ll0eSpJj9soWybJLDH7z8709ad+8biO3K/wXXRzz0V7Zg7ONoOyY3btwInU6H+++/H0Do6halc+wfpM7K9XjJOFrwMh95w++H8Pz973/HxIkT5en3338ft9xyC958803ceOON2LlzZyDy52LMmDFoaGhwmnf48GHk5ubKlxGMRiOOHDnilKahoQFGozEoeSLyRDCkQNq1CfYNq2HfsBrYX4PoPscOvm+//RZjx46Vp0NVt5BD4NBtCXdWiELC7yDozJkzSE5OBnD5tvW2tjbcfvvtiI+Px8yZM/Hdd995tZ2LFy+ipaUFLS0t8rZaWlrkp0hXVVVh9erVcvrZs2fDZDLJ4wTV1dWhrq4Od955p5zmjjvuwOHDh7F9+3acOnUK27dvx5EjR1w6SxOFgtTRDvvp7y//Wczhzk7EO3v2rNxq2/u8Qn/qFiKigfh9Oeyqq67C+fPnAQBfffUVhgwZgtGjRwO4fGnp4sWLXm2nqakJZWVl8vSGDRsAANOmTUNpaSnMZrMcEAHAsGHDsHTpUlRWVqKmpgYGgwEPPvigfHs8AOTl5WHx4sXYtGkTNm/ejPT0dCxevJhjBBFFAZ1OJ9cfTU1NftctREQD8TsIGjVqFGpqapCWloba2lqMHz9e7kxnMplc7s7yZPz48XjnnXc8Lu/tc+Ro3LhxKC8v73e7U6ZMcQqMiCg6jBgxAvv374fBYMAnn3zid91CRDQQvy+H/cu//AsaGxvxxBNPoKWlBXfddZe87K9//StycnICkkEiUpbZs2ejqakJL774Ik6dOsW6hYiCxu+WoAkTJqCiogInTpxAdnY2hg8f7rQsOzs7EPkjIoUxGo1YunQpvvvuO2RnZ+Paa6+Vl7FuIaJAGtSobGlpaW4HE7ztttsGs1kiUrjk5GQkJye7DFvBuoWIAmnQQ9N2dXWhvb0dPVcGfHM0bty4wW6eiBSqu7sb3d3daG9vd1nGuoWIAsHvIMhsNmP16tX4v//7P49pwvV8MSKKXl1dXfjLX/6Cr7/+2mMa1i1EFAh+B0Hr1q1Dc3Mz7r//fvzoRz/yONoyEZEvtm7dipMnT+LOO+/EqFGjkJqaGu4sEVGM8jsI+uqrr/Cv//qvmDFjRiDzQ0QK19TUhHnz5uHGG2+EVqvlU+SJKGj8vkUeAFJSUgKVDyIiGccCIqJQ8DsIuummm/DXv/41kHkhIsLEiRNx9OjRcGeDiBTA78thN910E/7rv/4LdrsdP/nJT+QnuDvKzc0dVOaISHkmTpyIzZs3Q5IkXH/99eju7nZJw7qFiALB7yBo+fLlAICamhrU1NS4TcM7OIjIV2vXrgUAHDhwAAcOHHCbhnULEQWC30HQr371q0Dmg4gIALBw4UL5tVqtxtVXXx3G3BBRLPM7CJo+fXoAs0FEdFlBQYH8mneHEVEwDerusF7ff/89/va3v+HixYuB2BwREQDg9OnTrFuIKGgG9diMjz/+GG+//TbMZjMA4IUXXkBubi5efvll5Ofn49Zbbw1IJolIWQ4dOoTq6mqcOXMGAOsWIgoOv1uCPv30U6xduxY5OTl4+OGHnZbl5ubi008/HXTmiEh5vvzyS7z99tvIzMzET3/6U6dlrFuIKJD8DoK2b9+O6dOn46mnnnI5Kxs5ciROnjw56MwRkfLs2bMHBQUF+MUvfoH/9//+n9My1i1EFEh+B0EnT57EzTff7HZZQkICzp4963emiEi5Tp8+jcmTJ7tdxrqFiALJ7yDoqquuwvnz590u6+zsxNChQ/3OFBEpl1arxYULF9wuY91CRIHkdxCUl5eHDz74AJIkuSzbt28fxo0bN6iMEZEy5eTkYP/+/axbiCjo/A6C7rnnHnz99df4/e9/j927dwO4fEfHH//4R3z11VeYP39+wDJJRMoxZ84cfPPNN3j55Zexb98+AKxbiCg4/A6CRo8ejaVLl+LixYvYuHEjAGDbtm344YcfsHTpUowaNSpgmSQi5Rg1ahR++ctf4tKlS3jvvfcAsG4houAY1DhBEyZMQEVFBVpbW9HV1YWrr74aGRkZgcobESnUmDFj8Pvf/x4Wi0V+dAbrFiIKNL+CoDNnzuDDDz/EV199JQ+UaDAYMH78eNx666181g8R+eXs2bM4ePAgmpqa0NXVBUEQkJaWxrqFiILC5yDoyJEjeOmll3DhwgWoVCpcffXVkCQJ33//PY4cOYKdO3diyZIl7LxIRD45fvw43njjDVy6dAmCICAhIQGSJKGtrY11CxEFhU9B0JkzZ7Bq1SrEx8fj3/7t3zBp0iRcddVVAIBLly7hiy++wMaNG/Hyyy+joqKCZ21E5JWzZ8+isrIScXFx+NnPfoZx48ZBp9NBq9UiMTGRdQsRBYVPHaPr6upgt9uxYsUKTJkyRQ6AgMvjBk2dOhUrVqyAzWZDXV1dwDNLsUsQBKc/UpbPPvsMdrsdjz/+OCZOnAidTicvY91CRMHiUxB0+PBhzJgxAykpKR7TpKamYvr06fjyyy8HmzdSCEEQgOrNsG9YDfuG1cD+GjAMUpZjx47hxhtvhF6v95iGdQsRBZpPQdCpU6dw7bXXDphu7NixOHXqlN+ZIuWxd7TDfvr7y38Wc7izQyF2+vRp5ObmDpiOdQsRBZJPQdC5c+eQmJg4YLrExEScO3fO70wRkbJcuHABCQkJA6Zj3UJEgeRTECSKIjSagftSq9VqiKLod6aISFlEUYRarR4wHesWIgokn2+R//7776FS9R87sbmaiHzV1tbmUrdoNBp0d3fL06xbyCOHmyp4cwV5y+cgaM2aNcHIBxEpXFVVVbizQFFMSDJA2loJm6UD6qxcCABcH8FL5MynIOhXv/pVsPJBRAq2cOFCt/N7H5lB5A27pQNSRzvs+uRwZ4WihE9B0PTp04OUDSJSsoKCArfztVot0tLSQpwbIlIKv58iT0RERBTNGAQRERGRIvn1FPlAq6mpwY4dO2CxWJCZmYni4mKMHTvWbdo1a9bg448/dpmfmZmJl19+GQCwb98+rF271iXNW2+95TQcPxERESlX2IOggwcPYv369Vi0aBHy8vKwZ88erFy5EhUVFUhNTXVJ/+CDD+L++++Xp202G5544glMmTLFKV1cXBxeeeUVp3kMgIiUpbGxETt27EBzczPMZjOWLFnisf+R4zqVlZU4efIkDAYD7rrrLsyePTtEOSaiUAr75bBdu3Zh5syZmDVrltwKlJqaitraWrfp4+Pjodfr5b+mpiacO3cOM2bMcEonCIJTuv6eSUREsenSpUvIzs7GQw895FX6trY2vPDCCxg7dizKy8tRVFSEN998E5999lmQc0pE4RDWliBRFHHixAncfffdTvPz8/Nx7Ngxr7ZRV1eH6667zuUOkosXL6KkpAR2ux3Z2dlYsGABcnJyPG7HarXCarXK04IgIC4uTn49kFgepCvYZRMEAQIgPzS1d3wPoc/rvssCkdYlL2HIg+OyQO7jWD4mvTVp0iRMmjTJ6/S1tbVITU1FcXExgMuX2ZuamrBz506X1mYiin5hDYLOnDkDu92OpKQkp/lJSUmwWCwDrm82m/Hll1/i17/+tdP8jIwMlJSUYNSoUbhw4QJ2796NZ599Fi+++CJGjBjhdlvbtm3Dli1b5OmcnByUl5f7fHtuenq6T+mjSbDKJooiOjVaoPdypUZz+U+nc37dd1kA06qDtF3f0mqRnJrq1aNpfBHLx2Sgff3118jPz3eaN3HiROzdu9fjY4MGewIVao558hSUx9rr/j4Hniw4U9r+CHufIMD9zvbmA9i3bx+GDh3qco3faDTCaDTK03l5eXjqqafw/vvve2wWLyoqQmFhocv7t7e3e/WsIkEQkJ6ejtbWVkhSbI1TGuyyCYIAu2iFvacHAKAWRUiiCHtPj9PrvssClVan08EWhO36mlYlWmEymQK2j2PhmNRoNCEdJ8hisbg9KbPZbOju7obBYHBZJ1AnUKFis9nQgSt9JDUaSFoNhFh+rdUiJS1twGfT8WTBmVL2R1iDoMTERKhUKpdWn66uLpeKqC9JkrB3717ccsstA545q1QqjB49Gq2trR7TaLVaaLVaj+/lLUmSovYHZyDBLJuEf1yekgaY5+m/P2nd5SPUeXCaF4R9HMvHZDD0PQHr3XeeTswGewIVaiqVCgKAnp4eqEQRklWEFMOvBasV7e3tHr8DsXCyEEixsD98OXkKaxCk0WiQm5uLhoYGp9achoYG3HDDDf2u29jYiNbWVsycOXPA95EkCd988w2ysrIGnWciil16vd7lpOzMmTNQq9VISEhwu06gTqBCRZIkuR+ap6A8ll73lnmgz4InC86Usj/CfjmssLAQr776KnJzc2E0GrFnzx6YTCbcdtttAC4/VLGzsxOPPvqo03p1dXUYM2YMRo0a5bLNd999F2PGjMGIESPkPkEtLS14+OGHQ1ImIopOY8aMwRdffOE07/Dhw8jNzQ14Xy0iCr+wf6unTp2K7u5ubN26FWazGVlZWVi6dKnclGU2m2EymZzWOX/+POrr6+U7OPo6d+4cXnvtNVgsFsTHxyMnJwdlZWW45pprgl0cIoogFy9edLoM3tbWhpaWFiQkJCA1NdXlJGv27NmoqalBZWUlZs2ahePHj6Ourg6PP/54uIpAREEU9iAIAObMmYM5c+a4XVZaWuoyLz4+Hm+99ZbH7RUXF3sMkIhIOZqamlBWViZPb9iwAQAwbdo0lJaWupxkDRs2DEuXLkVlZSVqampgMBjw4IMP8vZ4ohgVEUEQEVEwjB8/Hu+8847H5e5OssaNG4fy8vJgZouIIkTYR4wmIiIiCgcGQURERKRIDIKIiIhIkRgEERERkSIxCCIiIiJFYhBEREREisQgiIiIiBSJ4wQREVHsEgSPD8UlYhBEREQxS0gyQNpaCZulAwCg0qdAmP9zBkIEgEEQERHFOLulA1JH++XXANThzQ5FEPYJIiIiIkViEERERESKxCCIiIiIFIl9gogihaDiXSxERCHEIIgoQgiGFEi7NsF+pQOnKiUNwtwFDISIiIKEQRBRBJE62mE//b08zevVRETBwzqWiIiIFIlBEBERESkSL4cREZFyuHmMBikXgyAiIlIMx8doqPQpUP3L/xfuLFEYMQgiIiJF6X2Mhh3sE6J0/PyJiIhIkRgEERERkSLxchgRkQL0dgZmp2Cif2AQRGHhWBGzUiYKLkEQIL23AXZLB4SsXAAchZwIYBBEYSAIAlC9WX48hDrXCIZBRMEldwbWJ0Md7swQRQgGQRQWdofHQwjJaWHODRERKRE7RhMREZEisSWIgqZvXx8+DZ2IiCIJgyAKir79flQpaRDmLmAgREREEYNBEAWNY78fCCqoeYuubwSV075iAElEFFgMgigkBEMKpF2bYO9o591gXnLcZ2xJIyIKPAZBFDLSlZYh3g3mPcmhNY13MRARBRbrVSIiIlIkBkFERESkSAyCiIiISJEiok9QTU0NduzYAYvFgszMTBQXF2Ps2LFu0x49ehRlZWUu8ysqKjBy5Eh5+rPPPsPmzZtx+vRpDB8+HAsXLkRBQUHQykBERETRJexB0MGDB7F+/XosWrQIeXl52LNnD1auXImKigqkpqZ6XG/VqlWIj4+XpxMTE+XXx48fx6pVq7BgwQIUFBTg0KFDqKiowPLlyzFmzJigloeIiIiiQ9gvh+3atQszZ87ErFmz5Fag1NRU1NbW9rteUlIS9Hq9/KdS/aMo1dXVyM/PR1FREUaOHImioiJMmDAB1dXVwS4OUXBcGTPI8Y+IiAYnrC1BoijixIkTuPvuu53m5+fn49ixY/2u++STT8JqtSIzMxPz58/HhAkT5GXHjx/H3LlzndJff/312L17t8ftWa1WWK1WeVoQBMTFxcmvByLE8ECA/pRNEAQIgDwekABAuvLf8XXfZaFO65LvCM2v6sqYQdKVEbiFlDSoCn/W77hBsXxMEhEFQliDoDNnzsButyMpKclpflJSEiwWi9t1DAYDHnnkEeTm5kIURXzyySdYsWIFli1bhnHjxgEALBYL9Hq903p6vd7jNgFg27Zt2LJlizydk5OD8vJypKX5NqZNenq6T+mjiS9lE0URnRotoNNdnqHRXP7T6Zxf910WprTqCMjDgGk72wGz6coyLZJTU6HRDPwVjuVjkohoMMLeJwhwf6bq6ew1IyMDGRkZ8rTRaITJZMLOnTvlIMgdSZL6PSMuKipCYWGhy/u3t7dDFEWvypCeno7W1taYG9XXn7IJggC7aIW9pwcAoBZFSKIIe0+P0+u+y8KRVqfTwRbmPPiaViVaYTKZBmwJivZjUqPR+HwiQkTkrbAGQYmJiVCpVC4tNF1dXS6tQ/0xGo3Yv3+/PO2u1WegbWq1Wmi1WrfLfPkBkSQpan9wBuJr2ST845KTNMA8T/+DndZdniM5v07zvPw8YvmYJCIajLB2jNZoNMjNzUVDQ4PT/IaGBuTl5Xm9nebmZqfLX0ajEUeOHHHZptFoHFR+iYiIKHaE/e6wwsJCfPTRR6irq8PJkyexfv16mEwm3HbbbQCAqqoqrF69Wk5fXV2NQ4cO4YcffsB3332Hqqoq1NfX4/bbb5fT3HHHHTh8+DC2b9+OU6dOYfv27Thy5IhLZ2kiIiJSrrD3CZo6dSq6u7uxdetWmM1mZGVlYenSpXI/ALPZDJPJJKcXRREbN25EZ2cndDodsrKy8PTTT2Py5Mlymry8PCxevBibNm3C5s2bkZ6ejsWLF3OMICIiIpKFPQgCgDlz5mDOnDlul5WWljpNz5s3D/PmzRtwm1OmTMGUKVMCkj8iIiKKPWG/HEZEREQUDgyCiIiISJEYBBERkTJdeQSNzWbjyOoKFRF9goiIiEJNSDLAvmU9Os6dgTQ0EcL8n3NMLYVhEERERIplt3RA3d0FuyhC7dAaxGBIGRgEERGR4glJBkhbK2GzdEClT2GrkEIwCCIiIsLlViGpox12AOpwZ4ZCgkEQEVGM6u3sy06/RO4xCCIiikGCIEB6b8PlPi9ZuRDg+cHBRErFW+SJiGKUfHmn2xLurBBFJAZBREREpEgMgoiIiEiRGAQRERGRIjEIIiIiIkXi3WFEFNNqamqwY8cOWCwWZGZmori4GGPHjnWb9ujRoygrK3OZX1FRgZEjRwY7q0QUYgyCiChmHTx4EOvXr8eiRYuQl5eHPXv2YOXKlaioqEBqaqrH9VatWoX4+Hh5OjExMRTZJaIQYxBERDFr165dmDlzJmbNmgUAKC4uxuHDh1FbW4v77rvP43pJSUkYOnSoV+9htVphtVrlaUEQEBcXJ78OF0EQ0PvuvWME9c2N43wlvPZU5r77RIByB5hU2gCbDIKIKCaJoogTJ07g7rvvdpqfn5+PY8eO9bvuk08+CavViszMTMyfPx8TJkzwmHbbtm3YsmWLPJ2Tk4Py8nKkpaUNKv+DZbPZ0KHVAjodoNFA0mog9L4GoHM3P9ZfA26XXZ7tMF+rRUpaGtRq5T48Iz09PdxZCAkGQUQUk86cOQO73Y6kpCSn+UlJSbBYLG7XMRgMeOSRR5CbmwtRFPHJJ59gxYoVWLZsGcaNG+d2naKiIhQWFsrTvWfQ7e3tEEUxMIXxgyAIsFmtkHp6oBJFSFZRfq0G0ONmfqy/BuB2mRaXg+be+YIowmQyyQ9QVdKDVAVBQHp6OlpbW6O23BqNxuuTEAZBRBTT3DXre2rqz8jIQEZGhjxtNBphMpmwc+dOj0GQVquFVqt1uyzcPyKSwx/geulH6pMu1l97KrMjCZefKG/fsl7RT5SXJEkRZeYt8kTRSFBd7vPh8EfOEhMToVKpXFp9urq6XFqH+mM0GtHa2hrg3AUPj4nAkB85YukId1YoiNgSRBSFBEMKpF2bYO9oBwCoUodBNXeBIpvvPdFoNMjNzUVDQwMKCgrk+Q0NDbjhhhu83k5zczP0en0Qchh4fGgqkW8YBBFFKamjHfbT3wMAhOQ04EpQpEpJgzB3QZhzFxkKCwvx6quvIjc3F0ajEXv27IHJZMJtt90GAKiqqkJnZyceffRRAEB1dTXS0tKQlZUFURSxf/9+1NfX43e/+104i+ETuQVDnxzurBBFPAZBRDHCMSjide7Lpk6diu7ubmzduhVmsxlZWVlYunSp3GnSbDbDZDLJ6UVRxMaNG9HZ2QmdToesrCw8/fTTmDx5criKQERBxCCIiGLanDlzMGfOHLfLSktLnabnzZuHefPmhSJbRBQBeMJIREREisQgiIiIiBSJQRAREREpEvsE0aA4jkXC27KJiCiaMAgivwmCAFRvdrotm4EQEcUynvjFFgZBNCh23pZNRArhOBilUh+nEWsYBBEREXlJHowSgHKfMR87ePJOREREisSWIAqMKw/0lCf58EYiIopwDIIoIPo+0FOdawTDICIiimQMgihgXB7oSUREFMHYJ4iIiIgUKSJagmpqarBjxw5YLBZkZmaiuLgYY8eOdZu2vr4etbW1aGlpgSiKyMzMxL333ouJEyfKafbt24e1a9e6rPvWW29Bp9MFqxhEREQURcIeBB08eBDr16/HokWLkJeXhz179mDlypWoqKhAamqqS/qvvvoK+fn5WLhwIYYOHYq9e/eivLwcK1euRE5OjpwuLi4Or7zyitO6DICIKBb13ojAGxKCQBC4f2NY2IOgXbt2YebMmZg1axYAoLi4GIcPH0ZtbS3uu+8+l/TFxcVO0/fddx8+//xzfPHFF05BkCAI0Ov1wcw6EVHYOQ7gp87KhQCAw/cFjpBkgLS1Ejbu35gU1iBIFEWcOHECd999t9P8/Px8HDt2zKtt2O12XLhwAQkJCU7zL168iJKSEtjtdmRnZ2PBggVOQVJfVqsVVqtVnhYEAXFxcfLrgcTymYKnsgmCAAGQ/6Qr/zHAdCSldSlrhOfX27QcrkBZ5AH89MnhzkpM4v6NXWENgs6cOQO73Y6kpCSn+UlJSbBYLF5tY9euXbh06RJuuukmeV5GRgZKSkowatQoXLhwAbt378azzz6LF198ESNGjHC7nW3btmHLli3ydE5ODsrLy5GW5ttdTunp6T6ljyZ9yyaKIjo1WkCnAzSay3+9lxz7m47AtOoIyEPg0mqRnJoKjUbj9nMjIqLLwn45DHB/purN2euBAwfw7rvv4oknnnAKpIxGI4xGozydl5eHp556Cu+//z4eeught9sqKipCYWGhy/u3t7dDFEWvypCeno7W1taYe5aMp7IJggC7aIW9pwdqUYQkirD39ABAv9ORllan08EWRfkdKK1KtMJkMgFA1B+TGo3G5xMRIiJvhTUISkxMhEqlcmn16erqcmkd6uvgwYP485//jN/+9rfIz8/vN61KpcLo0aPR2trqMY1Wq4VWq3W7zJcfEEmSovYHZyDuyiU5/MHL/5GU1l15Ijm/3qZ1/Kxi+ZgkIhqMsAZBGo0Gubm5aGhoQEFBgTy/oaEBN9xwg8f1Dhw4gP/8z//E448/jsmTJw/4PpIk4ZtvvkFWVlZA8q1UgiBwVOho4PAIE1EUL39uDIKIiFyE/XJYYWEhXn31VeTm5sJoNGLPnj0wmUy47bbbAABVVVXo7OzEo48+CuByALRmzRoUFxfDaDTKrUg6nQ7x8fEAgHfffRdjxozBiBEj5D5BLS0tePjhh8NSxlhi56jQEa/3ESZSRzs6NVpISXoIcxcwECIi6iPsQdDUqVPR3d2NrVu3wmw2IysrC0uXLpX7AZjNZrl/AwDs2bMHNpsN69atw7p16+T506ZNQ2lpKQDg3LlzeO2112CxWBAfH4+cnByUlZXhmmuuCW3hiMJEfoSJTgdJtLLFjojIjbAHQQAwZ84czJkzx+2y3sCm13PPPTfg9oqLi13GEyJSLIfLY73YKkREFCFBEBEFT+/lsd6+XKqUNF4eIyICgyAiRZAc+nIBfHIy0aA5PE4DYOtqtGIQRC74xSYi6p/j4zRU+hQI83/O+jIKMQgiJ4IgANWbYe9oly+bEBGRK/lxGoIADU8eoxKDIHLheBs8L5sQRaZYfl5htGGrUPRiEEREFGX45PjII7cKAVCHOzPkNZ7oExFFIflHt9sS7qwQRS0GQURERKRIvBxGREQUKH1unXfEfkKRh0EQERFRgDh2klZn5QLdXewwHcEYBBEREQWQ3F9LnwxYOtlhOoKxTxAREREpEoMgIiIiUiReDiPPHJ4+LooiB2UjIqKYwiCIPOp9+rjU0Y5OjRbSqBwwDCIioljBIIj6JT99XKeDlKgPd3YoBPq2+PFuFqLA44OqIwODICKlcbjMCThXwI4P0AUgP0SXlTRR4Dg+9oS3zocXgyAihem9zGnvaHcb5Dg+QBeCCmqesRIFHJ81FhkYBBEpkOQQ6PR3i+hAARMRUTRjEEROl0Z4Bxj15W3AREQUbRgEKVzfPiDqXCPvACMiIkVgEEROfUCE5LQw54aIiCg0GAQRkXf63FUGsKN0KPBWaqLgYRBERF5x7CQN8Pb5UOCt1ETBxSBIgdgRmmR9WncGOh4cO0kD7CgdCvKt1IIAzZXPh9/bKCQI8ufW3+fHlr/QYhCkMOwITY76tu7weIhcQpIB0tZK2CwdUGflQgDAn8jo4c3nx5a/0GMQpBCOZyA2doQmBxKPh6ghtwrpk8OdFfKD28+vTwuRjYMohhSDIAVwbP3hmT4RUeRgC1948ZK+QvTeBm+3mMOdFSIiciC3EHVbwp0VxWFLUIxi52cioijmcJkMYCfpYGEQFCP6Bj3s7EpBx3GDiILG8TIZO0kHD4OgKOLpB8fdHV/s7ErBxnGDiILL3fAIAE82AolBUJToG+j0/cHhoy8oHDhuEFHwsVUoeBgERRE7f3CIFMGbQfVIWewebp1nv6HBYRBERBQBHAMf+9ZK2HnLNA2AgysOHoMgIgoKdpr2nuOPmTorFxIHRSQvBarfkGMQrqTvKoOgaOVwZw6bzCnSDNSHjVxxNGgajMH0G3IKwvUpsP/bEt/eO4ovyUVEEFRTU4MdO3bAYrEgMzMTxcXFGDt2rMf0jY2NqKysxMmTJ2EwGHDXXXdh9uzZTmk+++wzbN68GadPn8bw4cOxcOFCFBQUBLsog+btweR4Zw5vgaeI0Ccwt0VIH7Zg1C9EYdPPg1g9tQo5cvxd8fS4Du+y4Xr5NhovyYU9CDp48CDWr1+PRYsWIS8vD3v27MHKlStRUVGB1NRUl/RtbW144YUXMGvWLDz22GM4duwYXn/9dSQmJmLKlCkAgOPHj2PVqlVYsGABCgoKcOjQIVRUVGD58uUYM2ZMUMohCAJEURxUU6Lj2bMqdRhUDmfO7lp7eu/M4d1gFAn6Dcz7jCkUqkoyGPVLILE1l3zl1YNY+6RBd5dLC1HfS7DydgQBNptNPibdfVc9Xr5F9D3vLOxB0K5duzBz5kzMmjULAFBcXIzDhw+jtrYW9913n0v62tpapKamori4GACQmZmJpqYm7Ny5U66kqqurkZ+fj6KiIgBAUVERGhsbUV1djcWLFwe8DL2DE3Z2WSAl6X1q9u87yKHNMbDhgIcUZTwF5o4BUigvjQWjfgkUjz9CRAPw5tKpUxpLp0sLkVPrj8N2hCQDLOtfha39NARDKlQOLTuO39nBXL719YQomCdQYQ2CRFHEiRMncPfddzvNz8/Px7Fjx9yu8/XXXyM/P99p3sSJE7F3716IogiNRoPjx49j7ty5Tmmuv/567N6922NerFYrrFarPC0IAuLi4qDRDLyLBEGAlJIG9VVxEBISoNLpvP9gD9dDOtsNAFClpsNmHAcMz4AwchSks2cAjfZyWn0yJOM4aIdnXJ4eOQpSahq0vWmvvO67LBBpMTwDarUGqvSMgG43EtKybCHMb+/xnKSHoNV69R3x5vvnSbDql778rTsEQYCUOhzSkCFQJadCumoIpPQuqIaPhJR6NqivVZcuQEgZHpL3iqTXAFyXpZ2F+sI5qFKGhz1/oXiNvx6EdP4sVCnDoL5mnMfjA2qNU3ohPgHCjdPkViS7m3WFhKQBf/8EQYBU/7HLNr1Nj4J/GrDu8KXeCGsQdObMGdjtdiQlJTnNT0pKgsVicbuOxWJxm95ms6G7uxsGgwEWiwV6vd4pjV6v97hNANi2bRu2bNkiT9988814/PHHYTAYvCtM0f3epetr5h3+rUdE/QpW/dLXoOqOBcUDpyGKRIM5du+YH9z0PoiI8fbcXQ/v7xp5v4+P8KA3evWkqKgI69evl/9+8YtfOJ3dDeTChQt46qmncOHCBa/XiRYsW3SK5bL5Itj1y2DrjlDjceGM+8OZ0vZHWFuCEhMToVKpXM7Kurq6XM7Gerlr0Tlz5gzUajUSEhI8pulvmwCg1Wqh1Wp9LkMvSZLQ3NwcVb3ivcWyRadYLps3glW/9DXYuiPUlH5c9MX94Uxp+yOsLUEajQa5ubloaGhwmt/Q0IC8vDy364wZM8Yl/eHDh5GbmytfBzQajThy5IjLNo1GYwBzT0SRLFj1CxHFjrBfDissLMRHH32Euro6nDx5EuvXr4fJZMJtt90GAKiqqsLq1avl9LNnz4bJZJLH8airq0NdXR3uvPNOOc0dd9yBw4cPY/v27Th16hS2b9+OI0eOuHSWJqLYFoz6hYhiR9hPbaZOnYru7m5s3boVZrMZWVlZWLp0KdLSLt9iazabYTKZ5PTDhg3D0qVLUVlZiZqaGhgMBjz44INOt6/m5eVh8eLF2LRpEzZv3oz09HQsXrw4aGMEAZebxO+5556oahb3FssWnWK5bN4KRv0S7XhcOOP+cKa0/SFISrnwR0REROQg7JfDiIiIiMKBQRAREREpEoMgIiIiUiQGQURERKRIYb87LBps27YNhw4dwqlTp6DT6WA0GvHAAw8gIyOj3/UaGxvlW20NBgPuuusuzJ49O0S59o4/ZTt69CjKyspc5ldUVGDkyJHBzK5PamtrUVtbi/b2yw+hzczMxD333INJkyZ5XCcaPjPA97JFy2dGg9fY2IgdO3agubkZZrMZS5YsQUFBwYDrRMNx7w9f90csf1di+bfMXwyCvNDY2Ig5c+Zg9OjRsNls2LRpE55//nm8/PLLGDJkiNt12tra8MILL2DWrFl47LHHcOzYMbz++utITEyMqNtt/Slbr1WrViE+Pl6eTkxMDHZ2fZKcnIz77rsP6enpAICPP/4Yf/rTn/CnP/0JWVlZLumj5TMDfC9br0j/zGjwLl26hOzsbMyYMQMvvfTSgOmj6bj3h6/7o1csfldi+bfMXwyCvPDMM884TZeUlGDRokU4ceIExo0b53ad2tpapKamori4GMDlM/Wmpibs3Lkzog4cf8rWKykpCUOHDg1m9gblJz/5idP0woULUVtbi6+//tptoBAtnxnge9l6RfpnRoM3adKkfls7+4qm494fvu6PXrH4XYnl3zJ/MQjyw/nz5wHA47OEAODrr79Gfn6+07yJEydi7969EEUxYofg96ZsvZ588klYrVZkZmZi/vz5mDBhQrCz5ze73Y5PP/0Uly5d8vj4lGj9zLwpW69o+swoNKL1uA82JXxXYvm3zFvRnfswkCQJlZWVuPbaazFq1CiP6SwWi8tDGpOSkmCz2dDd3Q2DwRDsrPrM27IZDAY88sgjyM3NhSiK+OSTT7BixQosW7ZswNajUPv222/xzDPPwGq1YsiQIViyZAkyMzPdpo22z8yXskXTZ0ahFW3HfbAp5bsSy79lvmAQ5KN169bh22+/xfLlywdMKwiC03Tv4Nx950cKb8uWkZHh1JHOaDTCZDJh586dEVdJZGRk4MUXX8S5c+dQX1+PNWvWoKyszGOwEE2fmS9li6bPjEIvmo77YFPKdyWWf8t8wVvkffDGG2/giy++wLJly5CSktJvWr1eD4vF4jTvzJkzUKvVXl1qCjVfyuaO0WhEa2trEHI2OBqNBunp6Rg9ejTuu+8+ZGdnY/fu3W7TRttn5kvZ3InUz4xCK9qO+3CIte9KLP+W+YpBkBckScK6detQX1+PP/zhDxg2bNiA64wZMwYNDQ1O8w4fPozc3NyIuobqT9ncaW5uhl6vD2zmgkCSJFitVrfLouUz86S/srkTLZ8ZBVe0H/ehECvflVj+LfMXgyAvrFu3Dvv378fjjz+OuLg4WCwWWCwW9PT0yGmqqqqwevVqeXr27NkwmUzy2Ap1dXWoq6vDnXfeGY4ieORP2aqrq3Ho0CH88MMP+O6771BVVYX6+nrcfvvt4SiCR1VVVfjqq6/Q1taGb7/9Fm+//TaOHj2KW265RV4ejZ8Z4HvZouUzo8G7ePEiWlpa0NLSAuDyLc4tLS0wmUwAovu494ev+yOWvyux/Fvmr+gP40KgtrYWAPDcc885zS8pKcH06dMBAGazWf5SAcCwYcOwdOlSVFZWoqamBgaDAQ8++GDE3VLoT9lEUcTGjRvR2dkJnU6HrKwsPP3005g8eXKosu2Vrq4urF69GmazGfHx8fjRj36EZ555Rr7TIVo/M8D3skXLZ0aD19TU5DTY34YNGwAA06ZNQ2lpaVQf9/7wdX/E8nclln/L/CVIvT2ciIiIiBSEl8OIiIhIkRgEERERkSIxCCIiIiJFYhBEREREisQgiIiIiBSJQRAREREpEoMgIiIiUiQGQURERKRIDIKIiCjmrFmzBqWlpeHOBkU4jhhNREQxp7W1FRcuXEBOTk64s0IRjEEQERERKRIfoEpR5Z133sGWLVvw4osvYuvWrTh8+DBUKhWmT5+OBx54AKdPn8abb76JY8eO4eqrr8bs2bMxb948AMC+ffuwdu1arF69GsOGDZO3efToUZSVlWHZsmUYP358uIpGRD44c+YM3n77bXz55Zfo6upCXFwcMjIycO+99yI/Px9r1qxBY2Mj1qxZI69z7tw5bNiwAYcOHYIoihg3bhweeughPPbYY7jnnnvw05/+FMDg6hkA6OnpwaZNm3DkyBG0tbVBpVIhIyMDd999N2644YaQ7yvyjEEQRaWKigrccsstuPXWW9HQ0IAdO3bAZrPhyJEjmD17Nu68804cOHAAf/nLX5Ceno4bb7wx3FkmogB69dVX0dzcjJ/97GfIyMjAuXPn0NzcjLNnz7pNb7fbUV5ejqamJtx7773Izc3F8ePHsXLlSo/v4W89I4oizp49izvvvBPJyckQRRFHjhzBv//7v6OkpATTpk0Lyj4h3zEIoqh06623orCwEACQn5+PhoYGfPDBB1iyZAkKCgoAAOPHj8df//pX7N+/n0EQUYw5duwYZs6ciVtvvVWe118ry5dffom//e1vWLRoEWbPng3gct2h0WhQVVXldh1/65n4+HiUlJTI27Hb7bjuuutw7tw57N69m0FQBGEQRFFp8uTJTtMjR47EN998g4kTJ8rz1Go10tPTYTKZQpw7Igq2a665Bh9//DGuvvpqXHfddcjNzYVG4/knrbGxEQAwdepUp/k333yzxyBoMPXMp59+it27d6OlpQWXLl2S52u1Wq/KR6HBIIiiUkJCgtO0RqOBTqeDTqdzmX/hwoVQZo2IQmDx4sV47733UFdXh82bN2PIkCEoKCjAAw88AL1e75L+7NmzUKvVLnVHUlKSx/fwt56pr69HRUUFpkyZgjvvvBN6vR5qtRq1tbXYu3evH6WlYGEQRIrRewYmiqLT/O7u7nBkh4gGITExEcXFxSguLobJZMLnn3+Ov/zlL+jq6sIzzzzjkj4hIQE2mw1nz551Cm4sFkvA87Z//34MGzYMv/nNbyAIgjzfarUG/L1ocDhYIilGWloaAOCbb75xmv/555+HIztEFCCpqam4/fbbkZ+fj+bmZrdpxo0bBwA4ePCg0/y+04Gi0WicAiCLxcK6JgKxJYgU45prrkFGRgY2btwIm82GhIQEHDp0CH/729/CnTUi8sH58+dRVlaGm2++GSNHjkRcXBz+/ve/48svv/R4E8TEiRORl5eHDRs24Pz58/LdYZ988gkAQKUKXJvAj3/8Yxw6dAivv/46pkyZApPJhK1bt8JgMOCHH34I2PvQ4DEIIsVQqVR46qmn8MYbb+C///u/odVqMXXqVDz00EP44x//GO7sEZGXtFotrrnmGuzfvx9tbW2w2WxITU3FvHnznMbrcdT7/d+wYQP+53/+B6IoIi8vD4899hieeeYZxMfHByx/M2bMQFdXFz788EPs3bsXw4YNw913342Ojg5s2bIlYO9Dg8cRo4mISLEOHDiA//iP/8CKFSuQl5cX7uxQiLEliIiIFOHAgQPo7OzEqFGjoFKpcPz4cezcuRNjx45lAKRQDIKIiEgR4uLicPDgQbz33nu4dOkS9Ho9pk2bhp/97GfhzhqFCS+HERERkSLxFnkiIiJSJAZBREREpEgMgoiIiEiRGAQRERGRIjEIIiIiIkViEERERESKxCCIiIiIFIlBEBERESnS/w9eW6ags6KPtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i, p in enumerate(['mu', 'sigma']):\n",
    "    plt.subplot(1, 2, i + 1)\n",
    "    sb.histplot(my_temp_alg.X.theta[p], stat='density')\n",
    "    plt.xlabel(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks reasonable!\n",
    "You can see from the output that the algorithm automatically chooses the tempering exponents $\\gamma_1, \\gamma_2,\\ldots$. In fact, at iteration $t$, the next value for $\\gamma$ is set that the ESS drops at most to $N/2$. You can change this particular threshold by passing argument ESSrmin to TemperingSMC. (Warning: do not mistake this with the `ESSrmin` argument of class `SMC`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=0, ESS=1000.00, tempering exponent=0.037\n",
      "t=1, Metropolis acc. rate (over 9 steps): 0.246, ESS=1000.00, tempering exponent=0.79\n",
      "t=2, Metropolis acc. rate (over 9 steps): 0.335, ESS=9537.70, tempering exponent=1\n"
     ]
    }
   ],
   "source": [
    "lazy_tempering = ssp.AdaptiveTempering(my_static_model, ESSrmin = 0.1)\n",
    "lazy_alg = particles.SMC(fk=lazy_tempering, N=1000, verbose=True)\n",
    "lazy_alg.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm progresses faster this time, but the ESS drops more between each step.\n",
    "Another optional argument for Class `TemperingSMC` is `options_mh`, which works exactly as for `IBIS`, see above. That is, by default, the particles are moved according to a certain (adaptative) number of random walk steps, with a variance calibrated to the particle variance.\n",
    "\n",
    "Again, all the indications above should be sufficient if you simply want to run a SMC sampler with random-walk Metropolis steps. If you are interested in more exotic SMC samplers (based on different MCMC kernels, or such that the state-space is not $\\mathbb{R}^d$), see the next tutorial on SMC samplers. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}