doc/examples/early_stopping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Early stopping of model simulations\n",
    "==================="
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "This notebook can be downloaded here:\n",
    ":download:`Early Stopping <early_stopping.ipynb>`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "For certain distance functions and certain models it is possible to calculate the\n",
    "distance on-the-fly while the model is running. This is e.g. possible if the distance is calculated as a cumulative sum and the model is a stochastic process. For example, Markov Jump Processes belong to this class. However, we want to keep things simple here and only demonstrate how to use the pyABC interface in such cases. So don't expect a sophisticated (or even useful) model implementation here.\n",
    "\n",
    "In this example we'll use in particular the following classes for integrated simulation and accepting/rejecting a parameter:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "* :class:`IntegratedModel <pyabc.model.IntegratedModel>`\n",
    "* :class:`ModelResult <pyabc.model.ModelResult>`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with the necessary imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# install if not done yet\n",
    "!pip install pyabc --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import tempfile\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import pyabc\n",
    "from pyabc import (\n",
    "    ABCSMC,\n",
    "    RV,\n",
    "    Distribution,\n",
    "    IntegratedModel,\n",
    "    LocalTransition,\n",
    "    MedianEpsilon,\n",
    "    ModelResult,\n",
    "    NoDistance,\n",
    ")\n",
    "from pyabc.sampler import SingleCoreSampler\n",
    "\n",
    "pyabc.settings.set_figure_params('pyabc')  # for beautified plots\n",
    "\n",
    "db_path = \"sqlite:///\" + os.path.join(tempfile.gettempdir(), \"test.db\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We define here a (very) simple stochastic process, purely for demonstrative reasons.\n",
    "First, we fix the number of steps *n_steps* to 30."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_steps = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We then define our process as follows:\n",
    "\n",
    "$$\n",
    "   x(t+1) = x(t) + s \\xi,\n",
    "$$\n",
    "\n",
    "in which $\\xi \\sim U(0, 1)$ denotes a uniformly in $[0, 1]$ distributed\n",
    "random variable, and $s$ is the step size, $s = $ step_size.\n",
    "\n",
    "The function `simulate` implements this stochastic process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate(step_size):\n",
    "    trajectory = np.zeros(n_steps)\n",
    "    for t in range(1, n_steps):\n",
    "        xi = np.random.uniform()\n",
    "        trajectory[t] = trajectory[t - 1] + xi * step_size\n",
    "    return trajectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We take as distance function between two such generated trajectories\n",
    "the sum of the absolute values of the pointwise differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance(trajectory_1, trajectory_2):\n",
    "    return np.absolute(trajectory_1 - trajectory_2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Let's run the simulation and plot the trajectories to get a better\n",
    "idea of the so generated data.\n",
    "We set the ground truth step size *gt_step_size* to "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_step_size = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "This will be used to generate the data which will be subject to inference later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4+UlEQVR4nO3deXxU1d3H8c8vCwSSQDa2EELY17CGvSgCKm5AUVFUFhWQWq1LW6WLT9VHW2p9rFSpSkUBRcBCEbTgjgrKGghrgLAkZA8J2SH7ef64kxAgQAhJZsnv/XrNayb3zkx+NwPfnJx77jlijEEppZTzcrN3AUoppa6NBrlSSjk5DXKllHJyGuRKKeXkNMiVUsrJaZArpZST0yBXdUZE3haR5+xdh1KuToNc1ZiIxIrIWRHJFZEsEflJROaIiBuAMWaOMeZ/q/k+Y+u+4msnIiNFJO+CmxGRO2377xWRwyKSLSJpIrJERJpVen2AiKwRkXwRiROR+y7zvZ4XkeILvlfHSvtHi8guEckRkeMiMrtuj145Kg1yda3uMMb4Au2BecCzwCL7llR3jDGbjDE+5TfgdiAP+Nz2lB+BEcaY5kBHwAN4qdJbLACKgFbA/cBbItLrMt9yZeXvZ4w5DiAinsAa4B2gOXAP8JqI9K21g1VOQ4Nc1QpjTLYxZh1WoEwXkd4islhEXgIQkSAR+czWcj8tIptExE1EPgBCgU9tLc5nbM//t4ik2Fq2P1QOO9v7LhCR/9r+GtgmIp0q7e8lIl/Zvk+qiPzett1NROaKyDERyRCRj0Uk4BoPfTqwyhiTb/s5xBtj0ivtLwU6276/N3An8JwxJs8YsxlYB0ytwfcNAJoBHxjLDiAa6FnzQ1HOSoNc1SpjzHYgARh5wa5f27a3wGqN/t56upkKnMRq2fsYY16xPX8D0AVoCewCll3wfvcCLwD+wFHgZQAR8QW+xmohB2OF6De21zwOTASut+3LxGohY3tt1mVucy88Vlsw3wUsuWD7z0QkG8jFCu7Xbbu6AiXGmCOVnr4HuFyL/A7bL6QDIvKL8o3GmFRgOfCgiLiLyDCsv4o2X+a9lIvysHcByiUlYbUYKysG2gDtjTFHgU2XewNjzHvlj0XkeSBTRJobY7Jtm9fYfmkgIsuA12zbbwdSjDH/Z/u6ANhmezwHeMwYk1DpfU+KyFRjTIkxxu8qj3MSkA58f0Htm4HmItIWmAXE2nb5ADkXvEc24HuJ9/8YWAikAkOA1SKSZYxZbtu/HHgXmG/7+hfGmPirPAblArRFrupCW+D0Bdv+htVy/tJ2Yu6iFm45Wwtznq0LJIdzQRhU6WkplR6fwQpJgHbAsUu8dXtgTXkrG6srohTrL4SamA4sNZeYec4Yk4j1l8EK26Y8rO6Qypphtdyrev1BY0ySMabUGPMTVmDfBSAi3W3vOw1ohNWqf0ZEbqvhsSgnpkGuapWIDMIK8vP+xDfG5Bpjfm2M6QiMB54WkTHluy94m/uACcBYrBN5YeVvX40S4rFOMl5q3y3GGL9KNy9b4FLFaJTKt99fcJztgFHA0ivU4wGU998fATxEpEul/X2BA9U4LrB+TuU/g97AEWPMF8aYMmPMYeC/wC3VfC/lQjTIVa0QkWYicjtWK/FDY8y+C/bfLiKdRUSwuhNKgTLb7lTOD19foBDIAJoCf76KUj4D2ojIkyLSWER8RWSIbd/bwMsi0t5WUwsRmVD+wgtGh1x4u7CGqcBPxpjzWv8icr+IhNoet8fqu//G9v75wH+AF0XEW0RGYP3C+qCqAxGRCSLiL5bBwK+Atbbdu4EutiGIYjvZezuw9yp+VspFaJCra/WpiORitXb/gNVX/WAVz+uCdRIyD9gC/NMYs9G27y/AH21dHr/BauXGAYnAQWBrdYsxxuQCNwJ3YHW/xAA32HbPxxol8qWt5q1Yfc81MY0LTnLa9AR+EpF8rKGIh7H6ycs9CjQB0rD6uH9hjDkA58aoV3ruvVjdUblYP5O/GmOW2I7zGPAQ8A+sfvfvgdVYfeaqgRFdWEIppZybtsiVUsrJaZArpZST0yBXSiknp0GulFJOrl6v7Bw3bpz5/PPPr/xEpZRSlV32Gop6bZGnp6df+UlKKaWuinatKKWUk9MgV0opJ6dBrpRSTs7u09gWFxeTkJBAQUGBvUtRqka8vLwICQnB09PT3qWoBsruQZ6QkICvry9hYWFY8ykp5TyMMWRkZJCQkECHDh3sXY5qoOzetVJQUEBgYKCGuHJKIkJgYKD+Ransyu5BDmiIK6em/36VvTlEkCullCsqLCllU8wp/vezgxSVlF35BTWkQQ68/PLL9OrViz59+tCvXz+2bbOWeHz99dc5c+ZMnX3fpKQk7rrrrjp7/6p89913NG/enH79+tGvXz9efPHFKp9njGH06NHk5FhLTKampnLffffRsWNHBg4cyLBhw1izZk19lk5sbCy9e/c+b9u+ffsqjiUgIIAOHTrQr18/xo4dW+33/Oijjyq+Xrx4MY899liVzx07diyZmZk1PwDVIKTlFvDxjnge+WAnA178iqmLtvPh1jiOpFa5ol+tsPvJTnvbsmULn332Gbt27aJx48akp6dTVFQEWEH+wAMP0LRp0zr53sHBwaxatapO3vtyRo4cyWeffXbZ56xfv56+ffvSrFkzjDFMnDiR6dOnV4ReXFwc69atu+h1JSUleHjU3z+r8PBwoqKiAJgxYwa33377Rb8cL1dTeZDfd999V/xeU6dO5Z///Cd/+MMfrrlu5TqMMRxIyuGb6DS+PZTKngRrffA2zb2Y2L8tY3q0ZFjHIJo0cq+zGhp8izw5OZmgoCAaN24MQFBQEMHBwfzjH/8gKSmJG264gRtusBaY+fLLLxk2bBgDBgzg7rvvJi/PWswlLCyMZ555hvDwcAYPHszRo0cv+j7ff/99Rcuxf//+5ObmntfCnDlzZsX+Fi1a8MILLwDwt7/9jUGDBtGnTx/+9Kc/1cePBIBly5YxYYK1Ctq3335Lo0aNmDNnTsX+9u3b8/jjjwNWK3b8+PGMHj2aMWPGcPr0aSZOnEifPn0YOnQoe/daq489//zzvPrqqxXv0bt3b2JjY4mNjaVHjx7MmjWLXr16cdNNN3H27FkAIiMj6du3L3379mXBggXVrn/UqFE8+eSTREREMH/+fGbMmHHeL00fH2ut5rlz57Jp0yb69evH3//+d8D6S2ncuHF06dKFZ555puI148ePZ/ny5SgFEBl3mt/9Zy9D//INt7+xmde/OYK7m/Dbm7ux4YmR/DR3NC//PJzR3VvVaYiDg7XIX/j0AAeTcmr1PXsGN+NPd/S65P6bbrqJF198ka5duzJ27Fjuuecerr/+en71q1/x2muvsXHjRoKCgkhPT+ell17i66+/xtvbm7/+9a+89tpr/M///A8AzZs3Z9++fSxdupQnn3zyohbvq6++yoIFCxgxYgR5eXl4eXmdt//dd60VuuLi4hg3bhwzZszgyy+/JCYmhu3bt2OMYfz48fzwww9cd9115732nnvu4fDhwxcd29NPP820adMu2r5lyxb69u1LcHAwr776Kr16Xfzz+fHHH3nnnXcAOHDgAAMGDLjkzxBg165d7N27l4CAAB5//HH69+/PJ598wrfffsu0adMqWs2XEhMTw/Lly/nXv/7F5MmTWb16NQ888AAPPvggb775Jtdddx2//e1vL/seFyoqKmLnzp2A1Vqvyrx583j11VcrPq/FixcTFRXF7t27ady4Md26dePxxx+nXbt2+Pv7U1hYSEZGBoGBgVdVi3ItG/Yl8/jy3TTxdOe6ri0Y3b0lo7q1INCnsV3qcaggtwcfHx8iIyPZtGkTGzdu5J577mHevHkX/cffunUrBw8eZMSIEYAVEsOGDavYP2XKlIr7p5566qLvM2LECJ5++mnuv/9+Jk2aREhIyEXPKSgo4O677+aNN96gffv2vPHGG3z55Zf0798fgLy8PGJiYi4K8pUrV1b7eAcMGEBcXBw+Pj6sX7+eiRMnEhMTc9HzTp8+ja+vb5Xv8ctf/pLNmzfTqFEjduzYAcCNN95IQEAAAJs3b2b16tUAjB49moyMjIq+9ksp79sGGDhwILGxsWRlZZGVlVVxvFOnTmXDhg3VPtZ77rmn2s+tbMyYMTRv3hyAnj17EhcXR7t27QBo2bIlSUlJGuQN2IZ9yTy2fDf92vmx+MFB+HrZ/0Iwhwryy7Wc65K7uzujRo1i1KhRhIeHs2TJkouC3BjDjTfeeMk/rSsPQatqONrcuXO57bbbWL9+PSNGjOCLL764qFU+Z84cJk2aVHGizhjD7373Ox555JHL1n81LfJmzZpVPL711lt59NFHSU9PJygo6LzneXh4UFZWhpubG7169aoIZoAFCxaQnp5ORERExTZvb+/L1lj5PctVHntd3rUF1udR3rVyLSrXVPl7l5WVVZwHqcqFtZSUlJxXc5MmTa65NuWc1tta4v3a+bHkocH4NHaMCG3wfeSHDx8+r0UaFRVF+/btAfD19SU31zrTPHToUH788ceK/u/8/HyOHDlS8bryVvHKlSvPa6mXO3bsGOHh4Tz77LMMGjSIQ4cOnbd/wYIF5ObmMnfu3IptN998M++9915FX3xiYiJpaWkXvffKlSuJioq66FZVt0pKSgrlC25v376dsrKyKluX3bp14/jx44DVqi4oKOCtt96q2H+50TwjR45k2bJlgDVKJigoiGbNmhEWFsauXbsAqyvmxIkTl3wPAD8/P/z8/Ni8eTNAxXvWRFhYGJGRkQCsW7eO4uJi4PzP+EqMMaSkpBAWFlbjOpTz+u9eK8T7O1iIg4O1yO0hLy+Pxx9/nKysLDw8POjcuTMLFy4EYPbs2YwbN47g4GA2btzI4sWLmTJlCoWFhQC89NJLdO3aFYDMzEz69OlD48aNq2y1v/7662zcuLGihXvLLbeQnJxcsf/VV1/F09Ozonthzpw5zJkzh+jo6IpfDD4+Pnz44Ye0bNmyxse7atUq3nrrLTw8PGjSpAkrVqyo8i+I2267je+++47OnTsjInzyySc89dRTvPLKK7Ro0aLiPEFVnn/+eR566CH69OlD06ZNWbJkCQB33nknS5cupVevXgwZMqTiZ3c577//Pg899BAiwk033VTj4541axYTJkygb9++jBs3rqK13qdPH9zd3enbty8zZszA39//ku8RGRnJ0KFD63VUjnIMn+1N4okVUQwI9eP9Bx0rxAGkvHVWHyIiIkz5yady0dHR9OjRo95qqAthYWHs3Lnzou4JZ5acnMy0adP46quv7F2Kw3jiiScYP348Y8aMuWifK/w7VlX7dE8ST660e4g7zgpBynm0adOGWbNmXfEkZUPSu3fvKkNcua7yEB8Y6s9iB2yJl3PMqpxMbGysvUuoE5MnT7Z3CQ5l1qxZ9i5B1aN1e5J4csVuIsICeH/GILwdNMRBg1wppS6yNiqRp1ZGOUWIgwa5UkqdpzzEB4UF8P6Dg2jayPFj8op95CLSTUSiKt1yRORJEQkQka9EJMZ2f+nT/Uop5QS+OJDCUyujGNzBeUIcqhHkxpjDxph+xph+wEDgDLAGmAt8Y4zpAnxj+1oppZzS/sRsnlwRRZ8QP96b4TwhDlc/amUMcMwYEwdMAJbYti8BJtZiXfWqIU1ju2zZMvr06UN4eDjDhw9nz5499fr9lXJEqTkFPLxkB/5NPVk4baBThThcfR/5vUD51S6tjDHlV7SkAK2qeoGIzAZmA4SGhtakxjrV0Kax7dChA99//z3+/v5s2LCB2bNnV/ziUqohOltUyswlO8krKGHVL4bT0tfryi9yMNVukYtII2A88O8L9xnrqqIqrywyxiw0xkQYYyJatGhR40LrSkObxnb48OEVVy8OHTqUhISEa35PpZxVWZnh6Y+j2J+Uzfx7+9OjTbMrv8gBXU2L/BZglzEm1fZ1qoi0McYki0gb4OJJQK7WhrmQsu+a3+Y8rcPhlnmX3N0Qp7Ett2jRIm655ZbL//yUcmGvfXWEDftT+MOtPRjbs8pOBadwNUE+hXPdKgDrgOnAPNv92lqsq940tGlsy23cuJFFixZVTEilVEOzZncCb248yr2D2jFzZAd7l3NNqhXkIuIN3AhUnk91HvCxiDwMxAHXfhngZVrOdakhTWMLsHfvXmbOnMmGDRt0Xm3VIO2MPc2zq/YxtGMAL07oXeX/WWdSrSA3xuQDgRdsy8AaxeLUDh8+jJubG126dAGqnsY2KCiIoUOH8stf/pKjR4/SuXNn8vPzSUxMrJjBb+XKlcydO/eK09iGh4ezY8cODh06VDHTIVx6GtvnnnuO+++/Hx8fHxITE/H09Lxo9sOraZGfPHmSSZMm8cEHH1Rr9kGlXE386TM88kEkbf2b8PYDA2nk4fxTTjnXGJs60NCmsX3xxRfJyMjg0UcfBawFFy6ckVIpV5VbUMzDS3ZQXFrGoukR+DVtZO+SaoVOY1sLXHEaW3V1XOHfsasrKS1j5tKdbI5JZ8lDgxnR2an+v16276fBt8iVUg3DS/+N5rvDp/jzz8OdLcSvSIO8FrjqNLZKObvsM8UcTs1l4+E0Fv8Uy0MjOnDfEMe7MPFaOUSQG2Oc/qyxarjqs3tSVe1sUSlH0/I4nJrL4ZQcDqfmcTglh9Scworn3NyrFX+4zTW7v+we5F5eXmRkZBAYGKhhrpyOMYaMjIyLhpKquhednMMb38YQnZxLbEY+5b9PG3m40aWlDyM6B9GtlS9dW/vSvbUvrZt5uWzG2D3IQ0JCSEhI4NSpU/YuRaka8fLyqvICL1V3Vkcm8IdP9tG0kQdDOgQwoV9wRWiHBXrj7uaagX0pdg9yT09POnRw7quqlFL1o7CklBc/PciybScZ2jGAf0zp75STXNU2uwe5UkpVR0LmGR5dtou9CdnMub4Tv7mpKx7uzn8xT23QIFdKObzvDqfx5MooSksN70wdyM29Wtu7JIeiQa6UclhlZYZ/fBvD/G9i6NbKl7ceGEiHIG97l+VwNMiVUg4pM7+IJ1dG8f2RU0wa0JaXJ4bTpJG7vctySBrkSimHsyc+i0eX7eJUbiF//nk4Uwa3c9mhg7VBg1wp5TByC4r5aNtJ/u/LI7TwbcyqXwyjT4ifvctyeBrkSim7i07O4cOtcazZnciZolJu6NaC1yb3w9/bNWYnrGsa5EopuygqKWPD/mQ+3BrHjthMGnu4cUffYKYObU/fdn72Ls+paJArpepVYtZZPtoWx8od8aTnFdE+sCl/uLUHdw0M0RZ4DVV3qTc/4F2gN2CAh4DDwEogDIgFJhtjMuuiSKWUcysrM2w6ms4HW+L49pC1fvvo7q2YOqw9IzsH4dbALqmvbdVtkc8HPjfG3CUijYCmwO+Bb4wx80RkLjAXeLaO6lRKOaFTuYX8OzKeFdvjOXn6DIHejfjFqE5MGRxKiH9Te5fnMq4Y5CLSHLgOmAFgjCkCikRkAjDK9rQlwHdokCvV4Blj2HIsg2XbT/LlgRSKSw1DOgTw65u6Mq53axp76Fjw2ladFnkH4BTwvoj0BSKBJ4BWxpjyRSdTgFZVvVhEZgOzAUJDXW9Cd6WU5XR+Easi41m+PZ4T6fk0b+LJtGFhTBkcSueWPvYuz6Vdcc1OEYkAtgIjjDHbRGQ+kAM8bozxq/S8TGOM/+Xeq6o1O5VSzssYw7YTp/lo20k+359CUWkZEe39uW9IKLeGt8HLU1vfteSa1+xMABKMMdtsX6/C6g9PFZE2xphkEWkDpF1bnUopZ5KZX8TsD3ayIzYTXy8P7hsSyn1DQunaytfepTU4VwxyY0yKiMSLSDdjzGFgDHDQdpsOzLPdr63TSpVSDiMx6yzTFm0jPvMsL03szZ0DQnQeFDuq7qiVx4FlthErx4EHATfgYxF5GIgDJtdNiUopR3I0LZepi7aTV1DC0ocGM7RjoL1LavCqFeTGmCggoopdY2q1GqWUQ9t9MpMHF+/Aw82NFY8MpVdwc3uXpNArO5VS1fT9kVPM+SCSFr6N+eDhwbQP1HnBHYUGuVLqitZGJfKbf++hc0tfljw4iJbNdJ1MR6JBrpS6rMU/nuCFzw4yKCyAf02LoHkTT3uXpC6gQa6UqpIxhr9/dYR/fHuUG3u24o0p/XVcuIPSIFdKXaS0zPDc2v18tO0kkyNC+PPPw3XFegemQa6UOk9hSSlPrYxi/b4UfjGqE8/c3E2XWXNwGuRKqQrGGOau3sf6fSn88bYezBzZ0d4lqWrQv5WUUhXmfxPDmt2J/PrGrhriTkSDXCkFwJrdCbz+dQx3DgjhsdGd7V2Ougoa5Eopth3P4NlV+xjaMYC/TArXPnEno0GuVAN3Ij2fRz6MJCSgCe88EEEjD40FZ6OfmFINWGZ+EQ++vx03Ed6fMYjmTfViH2eko1aUaqAKS0qZ/cFOkrILWD5riM6d4sS0Ra5UA2SM4ZlVe9kRm8mrd/dlYPsAe5ekroEGuVIN0Otfx7A2Konf3tyN8X2D7V2OukYa5Eo1MP/ZlcD8b2K4a2AIj47qZO9yVC3QIFeqAdl2PINnV+9lWMdA/vxzHWboKqp1slNEYoFcoBQoMcZEiEgAsBIIA2KBycaYzLopUyl1rY6dymP2B5GEBjTl7QcG6jBDF3I1n+QNxph+xpjyJd/mAt8YY7oA39i+Vko5oOwzxTy8eAcebsL7MwbrMEMXcy2/kicAS2yPlwATr7kapVStKykt47Hlu0jMOss7UwcSGtjU3iWpWlbdIDfAlyISKSKzbdtaGWOSbY9TgFZVvVBEZovIThHZeerUqWssVyl1tf6y4RCbYtJ5eWI4EWE6zNAVVfeCoJ8ZYxJFpCXwlYgcqrzTGGNExFT1QmPMQmAhQERERJXPUUrVjX/vjGfR5hPMGB7G5EHt7F2OqiPVapEbYxJt92nAGmAwkCoibQBs92l1VaRS6upFxmXyhzX7+VnnIP54Ww97l6Pq0BWDXES8RcS3/DFwE7AfWAdMtz1tOrC2ropUSl2d5OyzPPJBJG38vHjzvv66TJuLq07XSitgjW28qQfwkTHmcxHZAXwsIg8DccDkuitTKVVdBcWlzF4aSUFxKR/NGoJf00b2LknVsSsGuTHmONC3iu0ZwJi6KEopVTPlc6jsT8rmX1Mj6NrK194lqXqgf28p5ULe/v446/Yk8ZubujG2Z5UDyZQL0iBXykV8E53KK18c4o6+wTqHSgOjQa6UCzialssTK6LoFdyMV+7so3OoNDAa5Eo5uewzxcxcshMvT3cWTo2gSSN3e5ek6pkGuVJOrKik8uX3Awj2a2LvkpQd6FJvSjmpTTGneH7dAY6dyueVO/voKj8NmAa5Uk4m/vQZXvrvQb44kEpYYFPemxHB6O46QqUh0yBXykmcLSrlre+P8c73x3AT4bc3d2PmyA409tA+8YZOg1wpB2eMYcP+FF7+bzSJWWcZ3zeY393anTbNtT9cWTTIlXJgR1JzeeHTA/x4NIPurX1ZOXsoQzoG2rss5WA0yJVyQDkFxbz+VQxLtsTi09iDFyf04r7BoTr5laqSBrlSDmR/YjbLt59kbVQS+UUl3DsolN/e3I0Ab534Sl2aBrlSdnamqIRP9yTx0baT7EnIprGHG7f3CebBEWH0btvc3uUpJ6BBrpSdHEzK4aPtcXyyO4m8whK6tPThT3f0ZFL/EF0cWV0VDXKl6tGZohI+25vMR9tOEhWfRSMPN24Pb8OUIaFEtPfXOVJUjWiQK1UPjDEs3x7PXzZEk1tQQqcW3jx3e0/uHNBWF35Q10yDXKk6lp5XyNzVe/k6Oo3hnQJ5cmxXBoVp61vVnmoHuYi4AzuBRGPM7SLSAVgBBAKRwFRjTFHdlKmUc/r2UCrPrNpLTkEJz93ekweHh+HmpgGuatfVDEp9Aoiu9PVfgb8bYzoDmcDDtVmYUs7sbFEpf/xkHw8t3kmQT2PWPTaCh3/WQUNc1YlqBbmIhAC3Ae/avhZgNLDK9pQlwMQ6qE8pp7M3IYvb3tjEh1tPMvu6jqx9bATdWzezd1nKhVW3a+V14BmgfCXXQCDLGFNi+zoBaFvVC0VkNjAbIDQ0tMaFKuXoSssMb313lNe/jqGFb2M+mjmE4Z2D7F2WagCuGOQicjuQZoyJFJFRV/sNjDELgYUAERER5mpfr5QziD99hqdWRrEzLpM7+gbz0oTeOhZc1ZvqtMhHAONF5FbAC2gGzAf8RMTD1ioPARLrrkylHEv22WLiMvKJyzhDTFoe720+gQDz7+3HhH5V/nGqVJ25YpAbY34H/A7A1iL/jTHmfhH5N3AX1siV6cDauitTqfpljOFUXiFxGWeITc/n5OkzxGWcscL79BmyzhSf9/zhnQJ55a4+hPg3tVPFqiG7lnHkzwIrROQlYDewqHZKUqp+FRSXcjQtj+jkHKKTczmUkkN0cg6ZlcLaTaCtfxPCAr25LbwN7QObEhrgTVhQU0IDmtK0kV6Soeznqv71GWO+A76zPT4ODK79kpSqO9lni9l9MpPo5Fyik3M4lJLDsVP5lJZZp2+8PN3o1roZN/dqTffWvoQFeRMW6E1b/yZ46hSyykFpM0I1CNlni3l303He23yC/KJSANr6NaF7a19u6tmaHm2a0b2NL2GB3rjrWG/lZDTIlUvLLSjm/R9j+dem4+QWlHBbeBvuHxpKrzbNdVSJchka5MolnSkqYclPcbzzwzGyzhRzY89WPDW2Kz2D9cIc5Xo0yJVLKSgu5cOtcbz9/THS84oY1a0FT9/YlT4hfvYuTak6o0GuXEJhSSkf74jnzY1HSc0pZHinQN6Z2pWB7QPsXZpSdU6DXDm9zTHpPLt6L4lZZxkU5s/r9/RnWCddaV41HBrkymmVlJbx+tcxLPjuKJ1a+LD0ocGM7BKk83yrBkeDXDmlpKyzPLFiNztiM5kcEcLz43vpRTmqwdJ/+crpfHUwld+u2kNxSZnObaIUGuTKiRSWlDJvwyHe/zGWXsHNePO+AXQI8rZ3WUrZnQa5cgpxGfk89tFu9iVmM2N4GL+7tTuNPdztXZZSDkGDXDm8dXuS+P1/9uHuJrwzdSA392pt75KUciga5MphnS0q5cXPDrB8ezwD2/vzjyn9aevXxN5lKeVwNMiVQ9qbkMXTH+/haFoej47qxFM3dtXZB5W6BA1y5VCKSsp489sYFnx3jBY+jfng4cGM7NLC3mUp5dA0yJXDOJSSw9Mr93AwOYdJA9rypzt60byJzlCo1JVokCu7Kykt450fjvP610do3sSThVMHcpOe0FSq2q4Y5CLiBfwANLY9f5Ux5k8i0gFrvc5AIBKYaowpqstiles5diqPX3+8h6j4LG4Nb81LE8MJ8G5k77KUcirVaZEXAqONMXki4glsFpENwNPA340xK0TkbeBh4K06rFW5kLIyw+KfYvnr54fw8nTnH1P6c0efNjpPilI1cMUgN8YYIM/2paftZoDRwH227UuA59EgV9UQf/oMv/n3HradOM3o7i2ZNymcls287F2WUk6rWn3kIuKO1X3SGVgAHAOyjDEltqckAFVOeCEis4HZAKGhoddar3JCZWWGQym5bDmewdbjGWyOScfdTXjlrj7cPTBEW+FKXaNqBbkxphToJyJ+wBqge3W/gTFmIbAQICIiwtSgRuVkysoMh1Nz2Xo8gy3HMth24jTZZ4sBaB/YlIn9g/nlDZ0J8W9q50qVcg1XNWrFGJMlIhuBYYCfiHjYWuUhQGJdFKicQ3L2Wb7Yn8LW46fZdiKDzDNWcIcGNOXmXq0Y2jGQoR0DCdYrM5WqddUZtdICKLaFeBPgRuCvwEbgLqyRK9OBtXVZqHJc+xKymfbeNjLPFBPi34QxPVoxrGMgQzoGaKtbqXpQnRZ5G2CJrZ/cDfjYGPOZiBwEVojIS8BuYFEd1qkc1I7Y0zz0/g6aNfFkxexhdGvta++SlGpwqjNqZS/Qv4rtx4HBdVGUcg6bYk4xa+lOgv2asGzmENo0124TpexBr+xUNfLlgRQe+2g3HVt48+HMIQT5NLZ3SUo1WBrk6qqtjUrk6Y/3EN62OYsfHIRfU70SUyl70iBXV+WjbSf5wyf7GNIhgHenD8Knsf4TUg2YMVCUB/npcCbDuuWnw5n0c9vKv77nQ2gWXCdl6P9CVW3vbjrOS/+NZlS3Frz9wEC8PHWpNdUAGAO5KZARA+kxkHHUdh8DOclQWlj169wbQdMg8A607kvrbioqDXJ1RcYY5n8Tw+tfx3BreGtev6c/jTx0kQflYgpz4fRxyDh2LqjTY6yvi3LPPc+jCQR2huD+0OMOW1gHVboPtO4b+UA9XbWsQa4uyxjDn9dH869NJ7hzQAh/vTMcD12pRzmrghw4fcwK7NPHIcN2f/oY5J86/7nN21mB3W8KBHaBoM7WfbO24OZY/wc0yNUllZUZnlu7n2XbTjJtWHuev6MXbm46L4pyEmezIGEHnNwK8dsgLdrqq67Mtw0EdIKu4yCwEwR0tN06QSPnuZhNg1xd0svro1m27SRzru/Es+O66eRWynEZA1lxcHIbxG+17tMOAgbEHVqHQ/dbrYAO6GiFtn8YNPK2d+W1QoNcVWnR5hMs2nyCGcPDNMSV4yk6A6n7ITHyXIs7N9na18gX2g2CnhMgdAi0jYDGPvatt45pkKuLrN+XzEv/PcjNvVrx3O09NcSVfRUXQNoBSNptu0VZ3SSm1NrfvB20HwGhQ6HdEGjVC9wa1ogqDXJ1nh2xp3lyZRT92/kx/97+uGufuKpvp45A3I/ngjvtIJTZlj5oEmCNFuk6zroP7g/Nq1wKoUHRIFcVjqblMXPJTkL8mvDu9EE6TlzVn5wk2LcK9n0MKfusbV7NraAe/nil0G5Xb0P6nIkGuQIgLbeA6e9tx9NdWPzgYF0AWdW9gmw4uM4K7xObAAPBA2DcPOh6M/h30NCuJg1yRX5hCQ8t3sHp/CJWPjKU0EDnGXalnExJIcR8CXs/hiNfWFdFBnSE65+F8LutsdrqqmmQN3DFpWU8umwXB5NyeHd6BH1C/OxdknI1pcUQuwkOrIGDa62WuHcLiHgQwidD2wHa8r5GGuQNmDGGP67Zz/dHTvHnn4czunsre5ekXEVxARzfaHWdHF4PBVnWJevdb4c+d0OHUeCu8VNb9CfZgL3x7VFW7ozn8dGduW9IqL3LUc6uMA+OfgXRn1rdJkV51gnLbrdCj/HQ6Qbw1MVH6kJ11uxsBywFWgEGWGiMmS8iAcBKIAyIBSYbYzLrrlRVm/69M57XvjrCpAFtefrGrvYuRzmrs5lw5EuIXgdHv4aSAmvyqN53Qs/xEHYdeOiJ87pWnRZ5CfBrY8wuEfEFIkXkK2AG8I0xZp6IzAXmAs/WXamqtnx3OI3f/WcfI7sEMW9SH73gR1WPMdYEU/HbbFdTbodT0dY+3zYwYJrV8g4dpt0m9aw6a3YmA8m2x7kiEg20BSYAo2xPWwJ8hwa5w/vPrgSeXb2XLq18+ef9A3Q6WnVpJYXWVZTxttCO33ZuhsDGza3L4HvfCR2vty6Dd7AZARuSq/q1KSJhWAsxbwNa2UIeIAWr66Wq18wGZgOEhmo/rL0YY3j96xjmfxPDsI6BvP3AQHy9PO1dlnI0xQXWuO7dyyBp17nFEAI6Quex1iXw7YZAi+4a3A6k2kEuIj7AauBJY0xO5T/HjTFGRExVrzPGLAQWAkRERFT5HFW3CktKmbt6H2t2J3LXwBD+/PNwbYmr8+Wnw45FsONfVqu7ZU8Y8gi0GwrtBoNPS3tXqC6jWkEuIp5YIb7MGPMf2+ZUEWljjEkWkTZAWl0VqWouM7+IRz6IZHvsaX5zU1d+eUNn7RNX55w6AlsXwJ4V1onKLjfBsMegw3U6ttuJVGfUigCLgGhjzGuVdq0DpgPzbPdr66RCVWOx6fk8uHgHiZlnmX9vPyb008mFFNZJy9hNsGUBHPkc3BtD33th6KPQsru9q1M1UJ0W+QhgKrBPRKJs236PFeAfi8jDQBwwuU4qVDWyI/Y0s5fuBGDZrCEMCguwc0XK7kqLYf9/YMubkLLXGiY46ncQ8TD4tLB3deoaVGfUymbgUn9jjandclRtWBuVyG//vZe2/k14f8YgwoJcYxUUVUNns2DXEtj6NuQmQVA3uOMf0GeyXqDjInSwpwsxxvDmt0f5v6+OMLhDAO88MBB/ncWw4cqMg21vw66l1lWWHa6DO+Zbo090xIlL0SB3EWeKSnjukwOs3pXApP5t+cud4TT20PnEG6TESPjpTWuCKhHoNQmGPwZt+tq7MlVHNMidnDGGdXuSmLfhEMnZBTw5tgtPjOmiI1MamrIy68Tlljet1XUaN4Nhj8KQOdA8xN7VqTqmQe7E9idm88KnB9gRm0nvts14Y0p/IvSkZsNSkA37V1sjUDKOWivo3Pxn6D8VvJrZuzpVTzTInVBGXiGvfnmYFTviCWjaiHmTwrk7op2ur9lQFGTD4Q1w4BM49o119WWbfnDnIug5Uec5aYD0E3cixaVlLPkplvnfxHC2qJSHRnTgV2O60LyJXmrv8qoK72YhMGgW9Po5hEToBTwNmAa5k/j+yCle/PQAx07lc13XFvzP7T3o3NLX3mWpunSl8G47UEefKECD3OHFnz7DC58e4OvoNNoHNuXdaRGM6dFST2a6qsI8K7z3r64ivCfqLIOqShrkDiwtt4C7395CbkExz47rzkM/C9Mhha6o+Ky1IPH+1dYiDSVnwTcYBs20tbw1vNXlaZA7qMKSUn7x4S6yzxaz+hfD6RmsIxBcSkkRHPvWCu/D660LdrxbQP8HoPcka9ZBDW9VTRrkDsgYw5/WHiAyLpMF9w3QEHcVBTmQsN1aTT76U6sP3MvPCu7ed0L7n+mIE1Uj+q/GAX24NY4VO+J57IbO3Nanjb3LUVfLGMhNgZR9kLLHdr/PWiYNoJEvdL/NCvCON+ialuqaaZA7mK3HM3jh04OM6d5SF0V2BGVl1gnH0iJr9sAqHxdB1klrRsHy0C5fEg3APwxa94G+91mXyXcYqZNVqVqlQe5AEjLP8OiyXbQPbMrf7+2Hm17gYx+Ju2D7QmvYX8nZ6r/OzdOaz7vLzdA63HbrDV7N66xUpUCD3GGcLSpl9tJIikvL+Ne0CJrpepr1q7TYmmRq2ztWP3YjH2ua12Ztwd0T3BvZblU99gSfVtY6ltpNouxAg9wBGGP47ao9RKfk8N70QXRs4WPvkhqOvFMQuRh2LoLcZGuR4XF/hX736VwlymlokDuAt74/xmd7k3l2XHdu6K6L3NaLpN2wbSHsX2X1cXcaYy22oHN1KydUnTU73wNuB9KMMb1t2wKAlUAYEAtMNsZk1l2ZrmvjoTT+9sVh7ugbzJzrO9q7HNeVnwGnDkHaQdi3CuK3gqc3DJgOg2dDCz2xrJxXdVrki4E3gaWVts0FvjHGzBORubavn6398lzbsVN5/Gr5bnq2acYrd/bRy+5rQ346pEVboX3qsO3+0AWjSDrAzX+B/vfriUjlEqqzZucPIhJ2weYJwCjb4yXAd2iQX5WcgmJmLd2Jp4cb70wdSJNGeul9tZSVQV4qZMVZS5mV32eesAL7TMa55zZuBi26Qddx1onIFt2tUSXN2upMgcql1LSPvJUxJtn2OAVodaknishsYDZAaGhoDb+dayktMzy5IoqTGWf4cOYQQvyb2rsk+ysugIIs62rHgmxrweCCbOsEZFYcZMbagvsklBae/1qf1tZY7e63nQvsFt2hWbAGtmoQrvlkpzHGiIi5zP6FwEKAiIiISz6vocgtKOaJFVF8eyiN/53Qi6EdA+1dUt0ryLF1c0RD2iHIiIGzmecH9oXhXJlXc/Brb7Wmu95shbZ/mLXNr51eXKMavJoGeaqItDHGJItIGyCtNotyVfGnzzBzyU6Onsrjfyf2ZurQ9vYuqXYV5Z/rl06LPtdXnR1/7jkeTSCwM3gHWl0cXs2tWxM/22M/28223acFNPG30wEp5RxqGuTrgOnAPNv92lqryEVFxp3mkQ8iKSwpY/GDgxjZpYW9S7p6hbmQFQ/ZCVY4Z5c/TrC25yQCtj+63BtBUFdoNwQGzoCWPazuDv8wcNPzAUrVpuoMP1yOdWIzSEQSgD9hBfjHIvIwEAdMrssind0nuxN5ZtVe2vh5sWL2IDq3dPALfkpLrMmeTmyC+O1WH3V2vNUFUpmbh9Wqbt4Own4GAR1sgd3DurBGZ/JTql5UZ9TKlEvsGlPLtbicsjLD378+whvfHmVIhwDefmAg/t4OeAl3Wak14dOJTRC7GU5ugcIca19gF6srJHQYNA+xbn6h1r1PK21dK+UAtMlUR84WlfLrf0exfl8KkyNCeGliOI08HOSKwbIySN1nhfaJTRD3ExTaWtuBXSD8LquFHTYSfPRKU6UcnQZ5HUjNKWDW0p3sS8zm97d2Z9bIjo5xsU9OMuz+ACKXQE6CtS2gk7UWZIfroP0IaKbznyvlbDTIa9n+xGxmLtlJTkExC6dGcGPPSw6xrx9lZXDie2tSqEPrwZRaixmMec5qcTdva9/6lFLXTIO8FhhjOJGez8bDp3j1i8P4N/Vk1Rw7r7OZnwFRyyDyfWtlmiYBMOyX1giSwE72q0spVes0yGsoNaeAH4+m8+PRDLYcSycpuwCAge39eeuBAbT09ar/ooyB+G2w8z1rUYTSQusk5ajfQ487wNMONSml6pwGeTVlnylmy/EMfjqWzo9H0zl2Kh8Av6aeDOsYyC9uCGJEp0A6BHnXX394fgakHbAuvEk9YIX4qUPWHCMDp8PAB6FVz/qpRSllNxrkl1FYUsraqCSWbTvJvoQsygw08XRncIcA7hnUjuGdgujZplndL8lWdMZ2eXs0pB48F955qeee0yTAWlZs6KPWiuyNHXysulKq1miQVyGnoJjl207y3o8nSM0ppHtrXx4f3YURnYPo186vbocRlhRB6n5IjDx3S4+h4opJjybWnCOdx0LLnlaLu2VPa0y3I4yMUUrVOw3ySlJzCnjvxxN8tPUkuYUljOgcyN/u6svILkF1011ijDX9akIkJO60Qjt577kJpLxbQkiE1cJu2RNa9dJL3JVSF9EgB2JSc1n4w3E+iUqktMxwa3gbHrmuE+EhdbDoQPFZOLgO9q+GhB1w9rS13bMptOkHQ2ZD24HQNsK6elJb2UqpK2iwQW6MYWdcJu98f4yvo9Pw8nRjyuBQZv6sI6GBdTA/eFIU7FpqLTNWmG1Nwdr9Niu0QyKs+Ul0bhKlVA00uOQoLCnlv3uTWbIljj3xWfg39eTJsV2YNiyMgNqeB+VsphXcu5ZAyj7w8IIe42HANOsqSl3kVylVCxpMkCdlneWjbSdZvv0kGflFdAzy5sUJvbh7YLvaXWbNGGsOk11LIXodlBRA63C49VVrDhOdW1spVctcOsiNMWw9fpqlW2L58mAqZcYwpnsrpg9vz4hOQVceNlh0xpoJMGGH1bdtys7dykptj0srbSuDuB+tE5iNm0O/+63Wd3C/ejlepVTD5JJBnl9YwprdiSzdEsuR1Dz8mnoyc2QHHhjSnnYBl+n/Li2BpN1w4js4/r11gU1pESDg7gniDuJmjRoRsR6ft83Nmod71FyrC6WRrsWplKp7LhPkxhgOJOWwelcCq3YmkFtYQu+2zXjlrj6M7xuMl2cV3SfGQPoROP6ddYvdfG4e7tZ9YMgj0HGUdZl7I+96PBqllKo+pw/yE+n5rItK4rOoOHwz9jHQ4xjzWnkyILgprb0FSSuCzwutC21KC6Gk0GpllxRa60vmpVhv5B8GvSdBh+utKV29g+x6XEopVV3XFOQiMg6YD7gD7xpj5tVKVVeQmlPAZ7vjOLTrB1qk72Co+0EecY/Bq7E1cRXpQIYbuDcGj0a2+8bWOpKV79sPt1rcHa+3glwppZxQjYNcRNyBBcCNQAKwQ0TWGWMO1lZxlWXnnmH7T1+Tsf9bgrN3cq8cwVsKwROKg3rg2XG6tapN6FBr3hEdk62UaiCuJe0GA0eNMccBRGQFMAGo9SDf9sY0eqd/zo1iXbp+yqczJZ3uhx43QPsReHoH1va3VEopp3EtQd4WiK/0dQIw5MInichsYDZAaGhojb5RWbN2HHQfT6s+Y2nXfywttP9aKaUq1Hn/gzFmIbAQICIiwtTkPYZNf7lWa1JKKVdyLdeIJwLtKn0dYtumlFKqHl1LkO8AuohIBxFpBNwLrKudspRSSlVXjbtWjDElIvIY8AXW8MP3jDEHaq0ypZRS1XJNfeTGmPXA+lqqRSmlVA3oPKpKKeXkNMiVUsrJaZArpZST0yBXSiknJ8bU6Bqdmn0zkVNAXA1fHoQ1HZYrcbVj0uNxfK52TK52PFD1MaUbY8Zd6gX1GuTXQkR2GmMi7F1HbXK1Y9LjcXyudkyudjxQs2PSrhWllHJyGuRKKeXknCnIF9q7gDrgasekx+P4XO2YXO14oAbH5DR95EopparmTC1ypZRSVdAgV0opJ+cUQS4i40TksIgcFZG59q7nWolIrIjsE5EoEdlp73pqQkTeE5E0EdlfaVuAiHwlIjG2e3971ng1LnE8z4tIou1zihKRW+1Z49UQkXYislFEDorIARF5wrbdmT+jSx2TU35OIuIlIttFZI/teF6wbe8gIttsebfSNk345d/L0fvIbYs8H6HSIs/AlLpa5Lk+iEgsEGGMcdoLGUTkOiAPWGqM6W3b9gpw2hgzz/YL198Y86w966yuSxzP80CeMeZVe9ZWEyLSBmhjjNklIr5AJDARmIHzfkaXOqbJOOHnJCICeBtj8kTEE9gMPAE8DfzHGLNCRN4G9hhj3rrcezlDi7xikWdjTBFQvsizsiNjzA/A6Qs2TwCW2B4vwfpP5hQucTxOyxiTbIzZZXucC0RjrbPrzJ/RpY7JKRlLnu1LT9vNAKOBVbbt1fqMnCHIq1rk2Wk/PBsDfCkikbbFqV1FK2NMsu1xCtDKnsXUksdEZK+t68VpuiEqE5EwoD+wDRf5jC44JnDSz0lE3EUkCkgDvgKOAVnGmBLbU6qVd84Q5K7oZ8aYAcAtwC9tf9a7FGP12Tl2v92VvQV0AvoBycD/2bWaGhARH2A18KQxJqfyPmf9jKo4Jqf9nIwxpcaYflhrHg8GutfkfZwhyF1ukWdjTKLtPg1Yg/UBuoJUWz9meX9mmp3ruSbGmFTbf7Qy4F842edk63ddDSwzxvzHttmpP6OqjsnZPycAY0wWsBEYBviJSPnqbdXKO2cIcpda5FlEvG0nahARb+AmYP/lX+U01gHTbY+nA2vtWMs1Kw88m5/jRJ+T7UTaIiDaGPNapV1O+xld6pic9XMSkRYi4md73ARrQEc0VqDfZXtatT4jhx+1AmAbTvQ65xZ5ftm+FdWciHTEaoWDtWbqR854PCKyHBiFNeVmKvAn4BPgYyAUa7riycYYpziBeInjGYX157oBYoFHKvUvOzQR+RmwCdgHlNk2/x6rT9lZP6NLHdMUnPBzEpE+WCcz3bEa1R8bY160ZcQKIADYDTxgjCm87Hs5Q5ArpZS6NGfoWlFKKXUZGuRKKeXkNMiVUsrJaZArpZST0yBXSiknp0GulFJOToNcKaWc3P8Dk+kagLHWKwMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gt_trajectory = simulate(gt_step_size)\n",
    "trajectoy_2 = simulate(2)\n",
    "\n",
    "dist_1_2 = distance(gt_trajectory, trajectoy_2)\n",
    "\n",
    "plt.plot(gt_trajectory, label=f\"Step size = {gt_step_size} (Ground Truth)\")\n",
    "plt.plot(trajectoy_2, label=\"Step size = 2\")\n",
    "plt.legend()\n",
    "plt.title(f\"Distance={dist_1_2:.2f}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "As you might have noted already we could calculate the distance on the fly.\n",
    "After each step in the stochastic process, we could increment the cumulative sum.\n",
    "This will supposedly save time in the ABC-SMC run later on."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "To implement this in pyABC we use the :class:`IntegratedModel <pyabc.model.IntegratedModel>`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with the code first and explain it afterwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyStochasticProcess(IntegratedModel):\n",
    "    def __init__(self, *args, **kwargs):\n",
    "        super().__init__(*args, **kwargs)\n",
    "        self.n_early_stopped = 0\n",
    "\n",
    "    def integrated_simulate(self, pars, eps):\n",
    "        cumsum = 0\n",
    "        trajectory = np.zeros(n_steps)\n",
    "        for t in range(1, n_steps):\n",
    "            xi = np.random.uniform()\n",
    "            next_val = trajectory[t - 1] + xi * pars[\"step_size\"]\n",
    "            cumsum += abs(next_val - gt_trajectory[t])\n",
    "            trajectory[t] = next_val\n",
    "            if cumsum > eps:\n",
    "                self.n_early_stopped += 1\n",
    "                return ModelResult(accepted=False)\n",
    "\n",
    "        return ModelResult(\n",
    "            accepted=True, distance=cumsum, sum_stat={\"trajectory\": trajectory}\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Our `MyStochasticProcess` class is a subclass of `IntegratedModel <pyabc.model.IntegratedModel>`.\n",
    "\n",
    "The `__init__` method is not really necessary. Here, we just want to keep\n",
    "track of how often early stopping has actually happened.\n",
    "\n",
    "More interesting is the `integrated_simulate` method. This is where the real thing\n",
    "happens.\n",
    "As already said, we calculate the cumulative sum on the fly.\n",
    "In each simulation step, we update the cumulative sum.\n",
    "Note that *gt_trajectory* is actually a global variable here.\n",
    "If *cumsum > eps* at some step of the simulation, we return immediately,\n",
    "indicating that the parameter was not accepted\n",
    "by returning `ModelResult(accepted=False)`.\n",
    "If the *cumsum* never passed *eps*, the parameter got accepted. In this case\n",
    "we return an accepted result together with the calculated distance and the trajectory.\n",
    "Note that, while it is mandatory to return the distance, returning the trajectory is optional. If it is returned, it is stored in the database.\n",
    "\n",
    "We define a uniform prior over the interval $[0, 10]$ over the step size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = Distribution(step_size=RV(\"uniform\", 0, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "and create and instance of our integrated model MyStochasticProcess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = MyStochasticProcess()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We then configure the ABC-SMC run.\n",
    "As the distance function is calculated within `MyStochasticProcess`, we just pass\n",
    "`None` to the `distance_function` parameter. \n",
    "As sampler, we use the `SingleCoreSampler` here. We do so to correctly keep track of `MyStochasticProcess.n_early_stopped`. Otherwise, the counter gets incremented in subprocesses and we don't see anything here.\n",
    "Of course, you could also use the `MyStochasticProcess` model in a multi-core or\n",
    "distributed setting.\n",
    "\n",
    "Importantly, we pre-specify the initial acceptance threshold to a given value, here to 300. Otherwise, pyABC will try to automatically determine it by drawing samples from the prior and evaluating the distance function.\n",
    "However, we do not have a distance function here, so this approach would break down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "abc = ABCSMC(\n",
    "    models=model,\n",
    "    parameter_priors=prior,\n",
    "    distance_function=NoDistance(),\n",
    "    sampler=SingleCoreSampler(),\n",
    "    population_size=30,\n",
    "    transitions=LocalTransition(k_fraction=0.2),\n",
    "    eps=MedianEpsilon(300, median_multiplier=0.7),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We then indicate that we want to start a new ABC-SMC run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:History:Start <ABCSMC id=1, start_time=2021-02-21 17:22:41.018237>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<pyabc.storage.history.History at 0x7f16e89e95e0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abc.new(db_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We do not need to pass any data here. However, we could still pass additionally\n",
    "a dictionary `{\"trajectory\": gt_trajectory}` only for storage purposes\n",
    "to the `new` method. The data will however be ignored during the ABC-SMC run.\n",
    "\n",
    "Then, let's start the sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:ABC:t: 0, eps: 300.\n",
      "INFO:ABC:Acceptance rate: 30 / 111 = 2.7027e-01, ESS=3.0000e+01.\n",
      "INFO:ABC:t: 1, eps: 140.6687680285705.\n",
      "INFO:ABC:Acceptance rate: 30 / 140 = 2.1429e-01, ESS=2.1150e+01.\n",
      "INFO:ABC:t: 2, eps: 72.25140375797258.\n",
      "INFO:ABC:Acceptance rate: 30 / 433 = 6.9284e-02, ESS=1.8744e+01.\n",
      "INFO:pyabc.util:Stopping: maximum number of generations.\n",
      "INFO:History:Done <ABCSMC id=1, duration=0:00:01.067589, end_time=2021-02-21 17:22:42.085826>\n"
     ]
    }
   ],
   "source": [
    "h = abc.run(minimum_epsilon=40, max_nr_populations=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "and check how often the early stopping was used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "594"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.n_early_stopped"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Quite a lot actually."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Lastly we estimate KDEs of the different populations to inspect our results\n",
    "and plot everything (the vertical dashed line is the ground truth step size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEHCAYAAAC0pdErAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKK0lEQVR4nO3deXxU5dnw8d+VfYNASEjIAoQlrGENOygKIiDivrVqqa27bbX6PFWrbV+1ta3ax9oqaq1ibaviDoiKIIissoVAWENYkrBkIQkhezL3+8cZIMAkmUlmMgm5vp9PzMw5933OlUhy5dyrGGNQSimlzuXj7QCUUkq1TpoglFJKOaQJQimllEOaIJRSSjmkCUIppZRDft4OwFXTp083X375pbfDUOoskydPBmDFihVejUOpBoirFdrcE0R+fr63Q1BKqXahzSUIpZRSLaPNNTEp1Ro98cQT3g5BKbfTBKGUG0ydOtXbISjldh5tYhKR6SKyW0QyROTResrcKCI7RCRdRP7ryXiU8pTU1FRSU1O9HYZSbuWxJwgR8QVeBi4DsoENIrLAGLOjTpm+wGPABGNMoYh09VQ8SnnSgw8+COgoJnVh8eQTxGggwxiTaYypAt4DrjqnzJ3Ay8aYQgBjTK4H41FKKeUCTyaIOCCrzvts+7G6koAkEVktIutEZLoH41FKKeUCb3dS+wF9gclAPLBSRJKNMUV1C4nIXcBdAN27d2/hEJWq37a8bWQWZ1JZW0mgb6C3w1HKrTyZIHKAhDrv4+3H6soG1htjqoH9IrIHK2FsqFvIGPM68DpASkqKbmChWoU3tr3BXzf/FYCD+QdJDE/0ckRKuZcnm5g2AH1FJFFEAoCbgQXnlPkU6+kBEYnEanLK9GBMSrnF8kPL+evmvzIzcSbzZ81nzE/GUHV5FVklWY1XVqqN8FiCMMbUAA8AXwE7gfnGmHQReUpEZtuLfQUUiMgOYDnwP8aYAk/FpJQ7VNdW8/zG5+kV3otnJj7DgC4D+Nc9/6JTv048vfZpb4enlNt4tA/CGLMYWHzOsd/UeW2AX9o/lGoTFmUu4lDJIV6e8jL+Pv4A7E/bz6W2S1lwZAHpBekM6jLIy1Eq1Xy6FpNSLvpw74f0Cu/FpLhJp489/vjjrPzHSkL9Q3l7+9tejE4p99EEoZQL9hbuJS0vjWv7XovI2asn+4ov1/a9lq8Pfk1RRZF3AlTKjTRBKOWCxfsX4yu+XNn7Sofnr+x1JTWmhiUHl7RwZEq5nyYIpVywImsFI6NHEhEU4fB8/4j+JIYnsnj/YofnlWpLNEEo5aSskiwyijKYnDC53jIiwozEGWw6tom8sryWC04pD9AEoZSTvs36FoDJ8ZPPO/fiiy/y4osvAnBJwiUArMpZ1VKhKeURmiCUctK6I+vo0bEHCR0Tzjs3bNgwhg0bBkC/zv3oGtKVldkrWzhCpdxLE4RSTqi11bLp2CZGxYxyeH7p0qUsXboUsJqZLoq/iLVH1lJdW92SYSrlVpoglHLCrsJdnKw+yahoxwnimWee4Zlnnjn9fmLcREqrS0nLT2upEJVyO00QSjlh49GNAKTEpDhVfmTXkWfVU6ot0gShlBM2HttIj4496Bri3KaHnYI60bdzXzYe0wSh2i5NEEo1whhDWl4aQ6OGulRvVPQotuZtpdqm/RCqbdIEoVQjck7mcLziuMsJIiUmhfKactLz0z0UmVKe5e0d5ZRq9bblbwMgOTK53jKvvfbaecdGRtv7IY5tZFjXYR6JTSlP0icIpRqRlpdGkG8QfTv3rbdMv3796Nev31nHIoIi6B3eW/shVJulCUKpRqTlpzGwy0D8fOp/4F64cCELFy4873hKTApbjm2hxlbjyRCV8ghNEEo1oNZWy97CvQzsMrDBci+88AIvvPDCecdTolMoqylj1/FdngpRKY/RBKFUAw6VHKK8ppz+Ef2bVP9UP8TmY5vdGZZSLUIThFIN2H18N0CTE0RUSBQxoTGnO7qVaks0QSjVgF3Hd+Hn40ev8F5NvkZyZLImCNUmaYJQqgG7CnfRO7w3/r7+Tb7GkMgh5JzMoaC8wI2RKeV5Og9CqQbsOb6HcbHjGi33zjvv1HsuOcqaP7Etf1uDmw0p1droE4RS9cgvzyevPM+p/oeEhAQSEs7fJwJgYJeB+IovaXm6sqtqWzRBKFWPPcf3AM51UL///vu8//77Ds8F+wXTt3Nf7YdQbY4mCKXqsavQmruQ1Dmp0bJz585l7ty59Z5Pjkxme/52bMbmtviU8jSPJggRmS4iu0UkQ0QedXB+jojkiUiq/eOnnoxHKVfsOr6LbqHdCA8Mb/a1kiOTOVl9kgPFB5ofmFItxGMJQkR8gZeBGcBA4BYRcTQd9X1jzDD7xxueikcpV+05vod+Ef0aL+iEIVFDANiat9Ut11OqJXjyCWI0kGGMyTTGVAHvAVd58H5KuU1lbSX7T+ynX2cnEkRNVaNFEsMTCfMP034I1aZ4MkHEAVl13mfbj53rOhFJE5EPRcThMBARuUtENorIxry8PE/EqtRZDhQfwGZs9Oncp/5ChQfh3R/As/FwcDUc2w6Htzgs6iM+DIocpAlCtSne7qReCPQ0xgwBvgbedlTIGPO6MSbFGJMSFRXVogGq9mlf0T4Aeof3dlwgZzO8OgkOfAcpP+bDPz3AhzeGwT+nQdoHDqsMiRzC3sK9lNeUeypspdzKkxPlcoC6TwTx9mOnGWPqTi19A/izB+NRymkZRRn4ii89O/Y8/2RxDvz3JggOhx8thM49iQQoexLevw0+uQuCwiFp2lnVkiOTqTW17CjYcXoRP6VaM08+QWwA+opIoogEADcDC+oWEJFudd7OBnZ6MB6lnJZZnEn3jt3PX2LDGPj8l1B1En74IXTuCcC8efOYN38B/HA+RA+GD++A4/vPqnp6RnWeNjOptsFjCcIYUwM8AHyF9Yt/vjEmXUSeEpHZ9mI/F5F0EdkK/ByY46l4lHLFvqJ9jpuXdn0Oe76ES34NUWc6sOfNm8e8efMgIBRu/i+IwKf3gq32dJnI4EhiQ2NJy9cZ1apt8GgfhDFmsTEmyRjT2xjze/ux3xhjFthfP2aMGWSMGWqMucQYo7uqKK+rqq0iqySLXp3OWcHVZoNvnoHIJBhzT/0X6JQA0/8Ih9bC1nfPOjUkaoh2VKs2w9ud1Eq1OgdPHKTW1J7/BLFzAeTthIt/Bb6NdN8N+wHEj4JlT0PlydOHkyOTOVp6lNyyXA9ErpR7aYJQ6hz7iu0jmDqdkyDWvQIRvWDQNY1fRAQufxZOHoXVfz19+NSEOe2HUG2BJgilzpFZlImP+NCjY48zB4+kQdZ6GPVT8PF17kIJo2DwdbDmb1CcDcCALgPw8/HTfgjVJmiCUOoc+4r2ER8WT5Bf0JmDG98Ev2Cr6ciBxYsXs3jx4vNPTPktmFr47gUAAn0D6d+5vy79rdoETRBKnWNf0b6zO6irK2D7xzBwNgR3dlgnJCSEkJCQ80907gHDfghb/g0nDgNWM1N6QTo1thpPhK+U22iCUKqOals1B08cpE+nOkts7F0ClcUw5KZ6673yyiu88sorjk9OfMga7mrvi0iOSqa8pvz0bG2lWitNEErVkXUiixpTQ6/wOk8Qae9DWDQkXlxvvfnz5zN//nzHJzv3gKE3w6Z5UHKMoZFDrctqP4Rq5TRBKFXHeSOYqkohYykMvKrxoa0NmfQw1FbB2r8R3yGezoGdtR9CtXqaIJSqY1/RPgQhMTzROpCxFGoqYMDshis2pktvGHQtbHwLqSgiOSpZh7qqVk8ThFJ1ZBZlEhsWS7BfsHVg50IIjoDu45p/8YkPWms4bXiD5MhkMoszKakqaf51lfIQTRBK1bGveN+Z5qXaGtizBPrNbF7z0ikxydB3Gqyby5BO/TAYtudvb/51lfIQTRBK2dXYajhQfODMEhs5G63RS30va7TuihUrWLFiReM3mfgQlBWQnLMdQbQfQrVqmiCUsssuyabKVnVmDkTGUhAf6FX/6CWX9RgPCWPpsP41EsN76sJ9qlXTBKGU3ekRTKeeIDKWWQvu1TM5rq7nn3+e559/3rkbTXwIirMY4tuBtLw0jDFNDVkpj9IEoZRdZlEmgPUEUVpg7S/de4pTdRctWsSiRYucu1HfadB1IEMO76KwspBDJYeaGrJSHqUJQim7zOJMYkJjCPUPhczlgIE+ziUIl/j4wMSHGJl/AIDNxza7/x5KuYEmCKXsztpFLmOZ1bQUO9wzNxt0LYmhsXQ2wsZjGz1zD6WaSROEUoDN2NhfvN+aIGcM7FsGvS5xfmlvV/n6IeN/zsiyUjblrPHMPZRqJk0QSgFHSo9QUVthzYHI3QEnj0HvS52uHxwcTHBwsGs3HX4rI21+5FTkc7T0qIsRK+V5miCUgtMrq/bu1BsOrLIOJl7kdP0vvviCL774wrWb+gczsv91AGza9bFrdZVqAZoglKLOCKbwXlaCCE+wVmH1sKTxjxBmM2xKf9fj91LKVZoglMKaA9ElqAvhAR3h4BroMcGl+k8//TRPP/20y/f1DYlgeFgCmyrz4NB6l+sr5UmaIJTCGuLau1NvyNsNZfnQ07UEsWzZMpYtW9ake4/sO5vMAH8Kvvldk+or5SmaIFS7Z4whsyjTal46aO9/cPEJojlGxlorxW7K3QKZK1rsvko1xqMJQkSmi8huEckQkUcbKHediBgRSfFkPEo5kluWy8nqk/YO6tXQIRYiejVe0U0GRQ4izD+UNeGRsOQJa3tSpVoBjyUIEfEFXgZmAAOBW0RkoINyHYBfANoAq7wis/hUB3UiHFxtNS+JtNj9/X38GdNtLGs6dsIc3Qap/22xeyvVEE8+QYwGMowxmcaYKuA94CoH5Z4G/gRUeDAWpep1OkHYfK35D01oXurSpQtdunRpcgzjY8dzpKqY/fEjYNlTUKkbCSnv82SCiAOy6rzPth87TURGAAnGmM8bupCI3CUiG0VkY15envsjVe3avqJ9hAeG0+WIffOenhNdvsZHH33ERx991OQYJsRZSWnNwMugNBe+e6HJ11LKXbzWSS0iPsBfgIcbK2uMed0Yk2KMSYmKivJ8cKpdySzOpHd4byR7A4R0gS59WjyGuLA4enbsyeqyLBh6C6z5uzWiSikv8mSCyAES6ryPtx87pQMwGFghIgeAscAC7ahWLS2zKNNagylrPcSPblL/w2OPPcZjjz3WrDjGxY5j49GNVF76JASEwsIHwWZr1jWVag5PJogNQF8RSRSRAOBmYMGpk8aYYmNMpDGmpzGmJ7AOmG2M0aUtVYs5XnGcwspCeod0g4K9kDC6SddZu3Yta9eubVYsE2InUFFbwebSQzDtaTi0BlL/06xrKtUcHksQxpga4AHgK2AnMN8Yky4iT4nIbE/dVylXnF6DqarSOtDEBOEOo2JG4e/jz+qc1TDsVug+3hr2WqIL+Snv8GgfhDFmsTEmyRjT2xjze/ux3xhjFjgoO1mfHlRL21+8H4BehYdBfCF2hNdiCfEPYXS30XyT9Q1GBGa/BDUVsPAX1hLkSrUwnUmt2qXyqlpW7sljwY7N+EkwFenryA/rx4bDFVTXeq/df2r3qWSVZLGncA9E9oUpv4U9X2pTk/IKTRCqXcnIPckjH2xlxNNfc/ub37PpyG6qy7vQtXgHCwsTuOHVtaQ8s5RnF+/kaLHzU3Pi4+OJj49vdnyXJFyCICw9tNQ6MOYe6DERvnwMirObfX2lXCGmjT26pqSkmI0btSVKuaa8qpYXluzmzdX7CfTz5erhsUwf3I0nN1/PJV2S+X+r36F45qusDZnMwrQjfLHtCL4+wh0TEnnosiSC/D20s5wDc76cQ3FlMZ9c9Yl14Ph+mDvB6h+57ZMWneWtLigu/8PRJwh1wTtSXM6Nr63ljVX7uXl0d7771SU8e+0QBiUIhZXH6Vtr/ZEUnjSB6YO78fIPRvDt/1zCNcPjeG1lJle89B1bs4paLN6p3aeSUZTBgeID1oGIRLj8GchcDhv/2WJxKKUJQl3QtmUXM/vvq8nMO8kbt6fwh2uSiQwLBLDa+YGk4lzo0M3aJMguISKEP18/lHd+MpqKahs3vLaWxduO1HufBx98kAcffNAtMU/pPgWAZYfqLB8+8sfWFqhLnoSCfW65j1KN0QShLlh7j5Vw25vrCfD14ZP7JzB1YPTZ5wv3AtD36C6IH+Ww6WZS3ygW/mwig2M7cv9/N/PGd5kO75Wamkpqaqpb4u4W1o3BXQaz5OCSMwdFYPbfwccfPr1PV3xVLUIThLog5RSVc/ub3+Pv68O7d44lKbrDeWX2FO6ha1AXOhcegoQx9V4rIjSA/945lhmDY3jm853MW73fk6EDMD1xOjsKdpwehgtAeBzMfA6y1sGav3k8BqU0QagLTnlVLXe8tYGTlTX8647RdO8S4rDc3sK99A2IsN40kCAAgvx9eenm4UwbGM3vFu7gky2eHVE0M3EmPuLDwn0Lzz4x5EYYcCUs/z0cS/doDEppglAXnKc/38HuYyW8/IMRDOjW0WGZGlsNGUUZJNXawDcAug1p9Lp+vj68dMtwxvfuwiMfpLEmI9/doZ8WFRLF2G5j+Tzzc2ymzrwMEZj1IgSFwyd3Q02Vx2JQqtEEISK+IrKrJYJRqrm+2HaE/64/xN0X9+KipPpX/j144iDVtmr6nsiH2OHgF+jU9YP8fXn99hQSI0N54N0tHC4qByApKYmkpCS3fA2nzOo1i8Olh9l8bPPZJ0Ij4cq/wtFtsPLPbr2nUnU1miCMMbXAbhHp3gLxKNVkuSUV/OqjNIYmdOKRaf0aLHt6BFPuPpfXXwoL9OPVW0dSVWPj3v9sprKmltdff53XX3+9ybE7MqX7FIL9glmUuej8k/2vgCE3w6r/g9ydbr2vUqc428TUGUgXkWUisuDUhycDU8pVf/h8JxXVNv7vxqH4+zb8T3tv4V78xJdeFWXWEt8u6tM1jOdvGMLWrCL+9IVn9m0I8Q9havepLDmwhIoaB7O6L/8DBHaARQ/psuDKI5xNEE8Cs4CngBfqfCjVKqzLLODT1MPcfXEvekWFNVp+T+EeevqH4w9NXsF1+uBu3Da2B2+u3s9VN9/OXXfd1aTrNGR2n9mUVJfwzaFvzj8Z2gUuexoOrdW1mpRHOJUgjDHfAruwNvnpAOy0H1PK66prbfzms+3Edw7mvsnO7Qa3p3APSTUGOnWHDjFNvvdjM/vTs0sIK75PY+cu9z9JjI4ZTWxoLJ9mfOq4wHD7suBfPwnlhW6/v2rfnEoQInIj8D1wA3AjsF5ErvdkYEo5693vD7Hn2El+e+UgggMaXzOpuLKYI6VHSCrJb3R4a2NCAvx44cZhVNXUcvB4WbOu5YiP+HBVn6tYd2QdR046mMktYs2NKC+Clc+7/f6qfXO2ienXwChjzI+MMbcDo7GanZTyqvKqWv72TQajEyOYOqCrU3XS8635A4NPND9BAIzs0ZmY8GByT1Sw6aD7/4qf3Xs2BsOCffV0+8UMhmE/hO9fh8IDbr+/ar+cTRA+xpjcOu8LXKirlMe8s+4AeSWVPDKtH+LkKqfpBVaCGFBVZS2x4QbxnYMJ8PPhiU+3U+Pm/STiO8QzOmY0n2Z8evaciLou/bW14dGyp9x6b9W+OftL/ksR+UpE5ojIHOBzYLHnwlKqcScra5i7Yh8XJUUxOjHC6XrpBen08A2ho28wRA92SywjRwzn0gmj2XnkBPPWHHDLNeu6us/VZJ/MZtOxTY4LdIyF8Q/A9o/g8Ba331+1T852Uv8P8DowxP7xujHmV54MTKnGvL3mAIVl1Tx8mWsT1NIL0hlYVQNxI8HXzy2xvPjii3zyr9e5pF8ULy7dS15JpVuue8rUHlMJ8w+rv7MaYPzPrRnW2heh3MTpZiJjzEfGmF/aPz7xZFBKNaaiupa3Vu9ncr8ohiZ0crpefnk+R0uPMrg4z23NS6eICE/OGkhFdS3/t3SPW68d7BfM5T0v5+uDX1NaXeq4UFBHGHMv7Fqk6zQpt2gwQYjIKvvnEhE5UeejREROtEyISp3vky055J+s4q6LerlUb0fBDgAGVVRA97Fui+fWW2/l1ltvpVdUGLeO7cF73x9i99ESt10f4Jq+11BeU85XB76qv9CYuyEgDL7TaUqq+RpMEMaYifbPHYwxHet8dDDGOF4FTSkPs9kM//guk8FxHRnXq4tLddPz0xHc20ENkJ2dTXa2tcLrL6b0JSzQj98vdu8SGEMih5AYnthwM1NIBIz6KWz/GPIz3Hp/1f7oYn2qzVm2K5fMvFLuuqi30yOXTkkvSKcXgYRE9LF+mXpA59AAfj6lLyv35LHajSu+ighX97maLblbzmxH6si4B6zFB9fqnhGqeXSxPtXm/OO7TOI6BTNzsGszoI0xpBekM6j8pFvmPzTktnE9iA0P4rmvdmOMcdt1r+x1Jb7iy2f7Pqu/UFgUJF8PafOtCXRKNZFHF+sTkekisltEMkTkUQfn7xGRbSKSKiKrRGSgq1+Aal/2HCvh+/3HuW1cD/waWZDvXMfKjpFfns+g0hJIcG8H9bkC/Xz52ZS+pGYV8c2u3MYrOCkqJIqJcRNZkLGA2oa2HR11J1SXwdZ33XZv1f54bLE+EfEFXgZmAAOBWxwkgP8aY5KNMcOAPwN/cT501R79Z91BAnx9uGFkvMt1t+ZtBWBwVZXbnyDGjRvHuHHjzjp2/ch4enQJ4fkle7DZ3PcUcXWfq8ktz2XN4TX1F4odZq1S+/0/dKVX1WSuLNZ3APC3v94AbG6wkrUcR4YxJtMYUwW8B1x1znXrjoQKBdz3U6QuOKWVNXy8OYeZyTF0CXNug5+6Nh3bRDC+DCAYIhveL8JVzz77LM8+++xZx/x9fXhoahI7j5xg8XYH6yg10cXxF9M5sHP9S2+cMvouOL4PMpe77d6qfXF2sb47gQ+B1+yH4oBPG6kWB2TVeZ9tP3bute8XkX1YTxA/r+f+d4nIRhHZmJeX50zI6gK0cOthSipruHVsjybV33RsE8NqwT8+BXxaZqWYK4fG0rdrGH/5eo/bluDw9/VnWs9pfJv9LeU15fUXHDgbQqNgwxtuua9qf5z9KbkfmACcADDG7AWcWxmtEcaYl40xvYFfAU/UU+Z1Y0yKMSYlKqr+bSTVhe3f6w/SL7oDI3t0drlucWUxewv3MvLEcY90UF933XVcd9115x339REenpZEZl4pn6Yedtv9pvWYRnlNOd9lf1d/Ib9AGPYD2LsESj23f7a6cDmbICrtzUQAiIgfjTcH5QAJdd7H24/V5z3gaifjUe1M+uFituec4Adjurs8tBVg87HNGAwjKyo80kFdUFBAQUGBw3OXD4phYLeOvLw8g1o39UWMjB5JRFBEw5PmwNqW1FZjrdGklIucTRDfisjjQLCIXAZ8ACxspM4GoK+IJIpIAHAzcFajqYj0rfP2CmCvk/GodubjzTn4+wqzh8Y2qf6mY5sIwIfkyiqIS3FzdA0TEX52aR/255eyKM09TxG+Pr5c1uMyvsv5jrLqBvahiB4IMUN0NJNqEmcTxKNAHrANuBtYbIz5dUMVjDE1wAPAV8BOYL4xJl1EnhKR2fZiD4hIuoikAr8EftSEr0Fd4KprbXyWmsOU/tF0Dg1o0jU2HttIMgEERg201ixqYZcPiiEpOoyXl2e4bUTT5T0vt5qZchpoZgIYeou1wmuuzndVrnE2QfzMGPMPY8wNxpjrjTH/EJFfNFbJGLPYGJNkjOltjPm9/dhvjDEL7K9/YYwZZIwZZoy5xBijK4yp83y3N4/8k1VcO+K8MQ5OKa0uZefxnYwsKWry/tPN5eMj3H9JH/YcO8mSHUfdcs0RXUfQJahL481Myddbe0WkveeW+6r2w9kE4egv+zlujEOpen20KYeI0AAm92vauIjU3FRsxsbIkyc8NoN6ypQpTJkypcEys4bEkhgZyt++yXDL7GpfH1+m9pjKqpxVVNRU1F8wrCv0mWrNrNY5EcoFja3meouILAQS686gFpEVwPEWiVC1a8Vl1Xy98xizh8YS4Ne0oanrj67HDx+GVVZ67AniySef5MknG96F19dHuG9yb9IPn3Db7OrJCZMprynn+6PfN1ww+Xo4kQPZG9xyX9U+NPYTtwZrxvQuzp5B/Uvgcs+GphQs2naYqhob141wfeb0KatzVjPctwMhwREQ4dry4O529fA44jsHu+0pYlTMKIL9glmZvbLhgkmXg48/7Gx0hRylTmtsue+DxpgVwFTgO/ss6iNYQ1ZdH2uolIs+3pxD365hDI5rWsdyblkuewr3MKH0pLX0RBOGyDpjxowZzJgxo9Fy/r4+3De5D6lZRaxyw0qvgb6BjOs2jm+zv2044QSFQ+9LYccCcOPigerC5uwz+0ogSETigCXAbcA8TwWlFMD+/FI2HSzkupHxTZr7ANbTA8DE/GyPdlCXl5dTXt7ArOY6rhsZR7fwIP62zD37NUxOmMzR0qPsKWxkF7uBs6H4EBxJdct91YXP2QQhxpgy4FrgFWPMDcAgz4WlFCxIPYwIXDWsaXMfAFZmr6RrQDhJ1dUeX+LbWYF+vtx9US++P3Cc9ZmOJ9e5YlL8JABWZK1ouGC/mdZoph3azKSc43SCEJFxwA+Bz+3HfD0TklKWxduOkNKjM93Cg5tUv6KmgtWHV3OJfyTi4wexw90cYdPdPLo7kWGB/O2b5j9FRAZHkhyZ3Hg/REgEJF4EOz7TZiblFGcTxIPAY8An9sluvQBdIlJ5TEbuSXYfK2FmcrcmX2PN4TWU15Rz6Yki6DYMAkLcFl9zBflbTxGrMvLZdLCw2de7KP4ituVvI7+8kX6NgbOtFV5z3bsdqrowOb3ctzFmNvCyiITZl/B2uPKqUu6weJu1PPaMwU1PEMsOLaODfxijctKhx7jGKzTDrFmzmDVrlkt1fji2OxGhAby0rPkrzEyKm4TBsO7IuoYL9ptpfd7zZbPvqS58zi73nSwiW4B0YIeIbBIR7YNQHnOqeSkmPKhJ9StrK1l+aDmXRAzGv7YKuo93c4Rne+SRR3jkkUdcqhMS4Medk3rx7Z48UrOKmnX//hH9CQ8MZ+3htQ0X7BAD3YZaK7wq1Qhnm5heA35pjOlhjOkOPAz8w3NhqfZsX95Jdh1tXvPSquxVlFSXMJMw60D3sW6Kzr1uG9eDTiH+/K2ZTxG+Pr6MiRnDusPrGp9fkTQdstZDmc51VQ1zNkGEGmNO9znY50aEeiQi1e4tTrM3LyXHNP0a+xcTERTBmLyDEDXA6qD1oMmTJzN58mSX64UFWk8Ry3blsi27uFkxjI8dT255LpnFmQ0X7Hs5GBtkLGvW/dSFz9kEkSkiT4pIT/vHE0Aj/wqVaprPtx1hZDNGLxVWFLI8azkzel6OX/aGVvv0cMrt43rQMciPl75p3lPEuFirn6XRZqbY4dZOc3sbWeRPtXvOJog7gCjgY+AjINJ+TCm3ynRD89JnGZ9Rbavm+oihUHkCeni2/6G5OgT585OJvfh6xzHSDzf9KSI2LJYeHXuw9kgjCcLHB/pOg71fQ21Nk++nLnyNLdYXJCIPAk9jdVCPMcaMNMY8aIxp/tg8pc5xavTSzCY2L9mMjQ/2fMCIriPoczzbOtjdsyOY3GHOhJ50CPJr9uzqcd3GseHoBqprqxsu2HcaVBTp4n2qQY09QbwNpGBtFDQDeM7jEal27fNtRxnRvVOTm5e+P/o9h0oOcUO/G+DgGghPgE4JjVf0svBgf348IZEv04+y88iJJl9nXOw4ymvKSc1Lbbhg70vBx09HM6kGNZYgBhpjbjXGvAZcD1zUAjGpdmp/fik7j5xoVvPS/N3z6RTYicu6T4VDa1vs6eHGG2/kxhtvbNY17pjQk7BAP15c2siaSg0YFTMKH/Fhw9FGngyCOkL8KMjU+a6qfo0liNPPqfYtRJXymDPNS01LEEdLj7L80HJm955NYHEOnDzm8Qlyp9x3333cd999zbpGp5AA7pzUi6/SjzV5XkSHgA4MiBjQeIIA6HUJHE7V4a6qXo0liKEicsL+UQIMOfVaRJr+HKyUA5+nHWF4907Edmpa89Lb6W8D8MMBP7SeHsDjE+ROKSsro6ysrNnX+cmkRLqEBvCnL3Y1eb+IlOgU0vLSqKytbLhg70sAA/u/bdJ91IWvsf0gfI0xHe0fHYwxfnVet/zO7+qCdSC/lB1HTnBFE58eCsoL+HDPh1zR6wpiw2Lh4FoIjoCofm6O1LGZM2cyc+bMZl8nLNCPBy7tw9rMgibvFzEqZhRVtirS8tIaLhg7AgLDYZ82MynHmraHo1Ju9vmptZeamCDe2fEOlbWV/CT5J9aBg6ut/gcPbRDkST8Y0524TsH8+cvdTXqKGB49HEHYeHRjwwV9/SBxktUPoau7Kgc0QahW4YvtRxiW0Im4JjQvFVcW897u97i85+UkhidCcTYU7oeeEz0QqecF+vnyy8uS2JZTzBfbj7pcv2NAR/pH9GfDMWf6ISZD0SE4rvNe1fk0QSivO1RQxvacpjcvvbvrXUqrS/lp8k+tA/u/sz4nTnJThC3v6uFxJEWH8dxXu6mutblcf1TMKCf7IS61PutoJuWAJgjldaeal6YPdn1yXFl1Gf/e+W8mx0+mX4S9v+HAd1b/Q9e2u+Cwr4/w6Iz+7M8v5V9rD7pcPyU6hcraSrblbWu4YEQvCO+u/RDKIY8mCBGZLiK7RSRDRB51cP6XIrJDRNJEZJmI9PBkPKp1WrztCEPjw0mIcH1Dn/m751NcWcydQ+48c3D/d9BzgrWkRAuZM2cOc+bMces1L+nXlYuSonhx6R4KTjbyJHCOEdEjEKTxZiYR6D0Z9q/UZTfUeTz2EyQivsDLWDOwBwK3iMjAc4ptAVKMMUOAD4E/eyoe1TodKihjW05xk+Y+VNZW8vaOtxnTbQxDooZYBwsPQPEh6Nmyczo9kSBEhN/MGkBZVS1/+dq1yXPhgeH0i+jHpqObGi/c6xJrzarDm5sYqbpQefJPrNFAhn33uSrgPeCqugWMMcuNMacGj68D4j0Yj2qFFm9v+uS4T/Z+Qn55Pncl33XmoJf6H/Lz88nPb9qw1Ib06dqB28f14N3vD7E9x7WF/FKiU0jNS6Wqtqrhgr0mA6LNTOo8nkwQcUBWnffZ9mP1+QnwhQfjUa3Q4m1HGNKE5qVqWzVvbn+TYVHDGBUz6syJA99ZS1lH9XdzpA27/vrruf766z1y7QenJhERGshjH2+jxoUO61Exo6x+iPxG+iFCIiB2mHZUq/O0ik5qEbkVa1FAh4sBishdIrJRRDbm5eW1bHDKY7KOl5GW3bTmpUX7FnGk9Ah3DrkTOTXXwRirLb3nxDY5/6E+4cH+/PbKgWzLKeZtFzqsR0aPdG4+BFhPEdkboPJk0wNVFxxPJogcoO4ymvH2Y2cRkanAr4HZxhiHPXHGmNeNMSnGmJSoqCiPBKta3qm1l1wd3mozNt7c/ib9I/ozKa5OU1LBPig5Aj3b7vDW+swa0o3J/aJ4YclucorKnaoTHhhOUuck5+ZDJF4MthprBVyl7DyZIDYAfUUkUUQCgJuBBXULiMhwrP2uZxtjcj0Yi2qFFm87QnKc681LK7JWcODEAe4YfMeZpweAAyutz4kX3qLDIsLTVw0G4H8/3IrN5tzM55SYFLbmbm18f4juY8E3UNdlUmfxWIKwr/76APAVsBOYb4xJF5GnRGS2vdhzQBjwgYikisiCei6nLjBZx8vY2sTmpbe2v0VcWByX9bjs7BP7v4MO3aBLHzdF2bokRITwxBUDWZ1RwNtrDzhVZ2T0SCpqK0gvSG+4oH8wdB8DmSuaHae6cPh58uLGmMXA4nOO/abO66mevL9qvb7Y3rTmpc3HNpOal8pjox/Dz6fOP19j4MAqqy3dC/0P9957b4vc55bRCSzdeYw/frGLSX0j6dO1Q4PlR0aPBGDTsU0M6zqs4YsnXgzfPA0n8yBMm3JVK+mkVu3P59uOMjiuI927uNa89Nb2t+gU2Ilr+l5z9om8XVCa67XlNW666SZuuukmj99HRPjjdcmEBvrxwH+3UF5V22D5iKAIeoX3YuMxZzqqL7E+azOTstMEoVpcdmEZW7OKXG5eyjqRxbfZ33JTv5sI9jtnUb9931ife012T5AuysrKIisrq/GCbtC1QxB/uXEou4+V8OtPtzW64uvI6JFsyd1Cra3hZELsMGv5b00Qyk4ThGpxX2yzViidOdi1BPH+7vfxFV9u7Odga89930CXvtCpuztCdNltt93Gbbfd1mL3m9yvK7+Y0pePN+fwn/WHGiybEp1CaXUpuwp3NXxRH1/78t+aIJRFE4RqcYvSDjM4riM9I0OdrlNeU87HGR8zpccUuoZ0PftkdQUcWH1mZdJ24ueX9mVyvyh+tyCd1Q1sLnS6H8KZZTcSL4aig3B8v7vCVG2YJgjVog4VWKOXrhwS61K9L/Z/QUlVCTf3u9nBRddCTTn0meKmKNsGHx/hpVuG0zsqjHve2cSuo453AY4OjSahQ4KT/RAXW5+1mUmhCUK1sIVphwG4YojzzUvGGN7d9S59O/c9/dfwWfYtAx//NrtBUHN0DPLnrR+PIiTQlx+/tYGs4473xR4ZPZLNuZuxmUaW6ohMsoYK63BXhSYI1cIWpR1hRPdOxHd2fvRSekE6u47v4uZ+N589Me6UfcutiV4BzjdZXUhiOwXz1pzRlFbWcMs/1jmcaZ0SnUJxZTH7ivY1fDERq5lp/0qwub5RkbqwaIJQLSYj9yQ7j5xglovNSwv2LSDAJ4DpidPPP1lyFI5t93rz0sMPP8zDDz/stfsPjO3Iv386huLyam55fR3ZhWc/SZx68nKumWkylBVAbiOT69QFTxOEajGL0g4j4lrzUrWtmi/3f8nkhMl0DOh4foFTw1u93EF95ZVXcuWVV3o1hiHxnXjnJ2MoLKviurlrzuqTiAuLIyY0hk3HnNkfwt4Poc1M7Z4mCNUijDEs3HqY0T0jiO4Y5HS9NTlrKKws5Mre9fzy3feNtbx3dLKbIm2a3bt3s3v3bq/GADAsoRMf3DMOQbhh7trTo5tEhJHRI9l4dGOj8yboGGsNGdbhru2eJgjVInYdLWFfXimzhrrWvLQwcyGdAjsxIW7C+SdtNitB9L60RbcXdeTuu+/m7rvv9moMp/SP6chH940nJjyI2/65ntdX7sMYQ0p0CgUVBRw84cSS4b0mWyu71jSy2ZC6oGmCUC1iUdphfH2EGYNjnK5TUlXC8kPLmd5zOv4+/ucXOJpmtZW3s/kPzojrFMwn909g+uAY/rB4F/f/dzMDOg8DnO2HuBiqSyHHibLqgqUJQnmc1bx0hPG9uxAZFuh0va8Pfk2VraqB5qVl1udTawips4QF+vHyD0bw+Mz+fLn9KA/8K4vwgAjn+iF6TgTx0Wamdk4ThPK4jQcLOXS8jGuGN7Tj7PkW7ltIj449SI6sp39hzxLoNhQ6RLshyguTiHDXRb35z0/HUlFloyA/juUH11Hd2NalwZ2h2zDtqG7nNEEoj/t4czYhAb5cPsj55qUjJ4+w8dhGZvWa5XjuQ2kBZH8PSTPcGOmFa1zvLnz54EUMihhOaW0+V726kIzcRrYX7XWx1cSk25C2W5oglEdVVNeyaOsRpg+OITTQ+e1HPt//OQCzes1yXGDvEjA26OdgboQXPPHEEzzxxBPeDqNB4cH+PDvzagCyK9K44qXv+Oeq/fXvTtdrsm5D2s55dMMgpb7ecYySyhquHxHvdB1jDAv3LWRE1xHEd6in3p4vICwGYoa6KdLmmTq1bex91adTH7oEdWHYsGJOZEXy9KIdfLX9KH++fsj5iycmjLG2Ic1cAUnTvBKv8i59glAe9fHmbGLDgxjbq4vTdXYc30FmcSZX9LrCcYGaKsj4BpIu9/rw1lNSU1NJTU31dhiNEhFGdxvN1oJNvHH7SJ6/YSg7j55gxl+/Y97qc54m/IOtJUx04b52q3X8dKkLUm5JBSv35nP18Dh8fJzfBnTRvkX4+/hzec/LHRc4uBqqSqBf6+l/ePDBB3nwwQe9HYZTxnYbS355PpnFmVw/Mp4lD13E6MQIfrdwB7f8Yx2HCuos09HrYmspk5N53gtYeY0mCOUxC1IPU2szXOtC81KNrYbF+xdzcfzFhAeGOy6050vwC7IWlVMuG9NtDADrj64HoFt4MPN+PIo/XzeEHYdPMP2vK3ln3UFrxnXiZKuSPkW0S5oglMd8tDmHoQmd6NM1zOk6aw+v5XjFcWb1rqdz2hjYvdjqQA1wbT9rZYkLiyM+LJ51R9adPiYi3Dgqga8euoiRPTrz5KfbufNfGykMHwhB4WfWvFLtiiYI5RHph4vZeeQE17o69yFzIeGB4VwUd5HjAkdSoegQDJjd/CDbsTHdxrDx6EZqbDVnHY/tFMy/7hjNb68cyMo9+cz8+xoKYibB3q91+e92SBOE8oj/rj9EoJ8PVw9zPkGUVpeeWVrD18HSGgA7PgMfv1bV/9AWje02lpPVJ9lRsOO8cyLCjyck8tG94wnw8+H3exOgNJfaw1u9EKnyJk0Qyu1KK2v4LPUws4bEEh5Szy96B74++DUVtRX1z30wBnYsgJ6TICTCTdG6xx/+8Af+8Ic/eDsMp43uNhqwmvTqkxwfzqKfTSSw/+XYjPDp/DcpLqtuqRBVK6AJQrndgq2HOVlZww/GdHep3qLMRSR0SGBoVD1zG3J3wPF9MPAqN0TpXuPHj2f8+PHeDsNpEUERDOoyiFU5qxos1yHInz/cOpnjnQbTu2gNV728iozckhaKUnmbRxOEiEwXkd0ikiEijzo4f5GIbBaRGhG53pOxqJbzn/UH6R/TgRHdOzld52jpUb4/8n39S2uA1bwkPtC/nicML1qzZg1r1rStGceT4ieRlp9GUUVRg+VEhMjhVzLUZx9+FQVc/fIalu441jJBKq/yWIIQEV/gZWAGMBC4RUQGnlPsEDAH+K+n4lAtKy27iO05J/jBmO71/6J3YPH+xRhM/c1LYDUvdR8PYVFuiNS9Hn/8cR5//HFvh+GSSXGTsBkbaw47kdiSpiEYPphSSmJkKHe+s5GXl2c0vvmQatM8+QQxGsgwxmQaY6qA94Cz2gaMMQeMMWmADo+4QMxbc4CQAF+udmH00qmlNYZGDaV7x3qapY6lQ97OVtm81FYN6jKIzoGd+S7nu8YLxwyFsGg65yzng3vGMXtoLM99tZv/+TCNqhr98b1QeTJBxAFZdd5n24+5TETuEpGNIrIxL09ndLZWuScqWLj1MDemJNAxyPnO6d2Fu8koymj46SFtPogvDL7WDZEqAF8fXybETWB1zmpqbbUNF/bxsZY22buUIKnhxZuG8dDUJD7clM0d8zZwokI7ry9EbaKT2hjzujEmxRiTEhXV+poXlOWddQepsRl+PKGnS/UW7VuEn48f03vWszKrzQbbPoA+UyE0svmBqtMmxU2isLKQ9IL0xgsPmG0tcbJvOSLCL6b25fkbhrIus4AbX13L4aJyzwesWpQnE0QOkFDnfbz9mLoAVVTX8p/1h5g6IJoeXUIbr2BXa6tl8f7FTIqbRKegTo4LHVwNJ3JgyI3uCVadNj52PD7i41wzU+LFEBgOOxecPnT9yHjm/Xg0OYXlXPPKatIPF3swWtXSPJkgNgB9RSRRRAKAm4EFjdRRbdQnW3I4XlrFHRMSXaq35vAa8srz6t9WFCDtfQgIg34zmxml57z44ou8+OKL3g7DZZ2COjE0aijfHHJiKQ2/AGv/jV2fQ+2ZJqWJfSP54N5x+Ihw46tr+XaPNgNfKDyWIIwxNcADwFfATmC+MSZdRJ4SkdkAIjJKRLKBG4DXRMSJ51zV2tTaDK+vzGRQbEfG9nJtAtunGZ/SObAzk+MnOy5QXWENbx1wZatee2nYsGEMGzbM22E0yWU9LmNP4R4OFB9ovPCA2VBRBAfOfuLoH9ORT+6bQPcuodwxbwPvbzjkkVhVy/JoH4QxZrExJskY09sY83v7sd8YYxbYX28wxsQbY0KNMV2MMYM8GY/yjEVph9mfX8rPLu3j0tDWoooilmct54peV9S/tMauRVB5otU3Ly1dupSlS5d6O4wmuazHZYA1k71RfaaAf6g15PgcMeFBzL97LON7d+FXH23jL0t26zDYNq5NdFKr1stmM/z9mwz6RXdg2kDn95wGa1vRals1V/e5uv5Cm+ZBp+5nlp1upZ555hmeeeYZb4fRJDGhMQyJGsKSg0saL+wfDH0vsxJ3bc15pzsE+fPmnFHcmBLPS99k8PD8rToMtg3TBKGa5cv0o+zNPcn9l/ZxaVMggM8yPmNAxAD6RfRzXCBvj9WUMfLHrWbnuAvVtB7T2HV8F4dOONE0NPg6KM2ztiJ1wN/Xhz9dN4RfXpbEx1tymPPW9zoMto3SnzrVZLU2w0vL9tIrKpQrkru5VHf38d3sPL6Tq/o0MPFt0zzw8YfhtzYvUNWoaT2sPaedeopIuhyCO0Pqf+otIiL8fIo1DPb7/ce5Ya4Og22LNEGoJvtkSw67jpbwiyl98XXx6eH93e8T6BvIFYn17DtdXW79AhowC8K6uiFa1ZBuYd1IjkxmyQEnEoRfICTfYI1mKi9qsOj1I+N5+47RHC7SYbBtkSYI1SQV1bW8sGQ3Q+LDuXJIrEt1T1SdYFHmImYmzqx/7kP6p9ZomZQ7mhuqctKMxBnsPL6T3cd3N1542A+gthLSP2606IQ+Ogy2rdIEoZrkn6v2c6S4gsdnDmhS30N5TTk397/ZcQFjYM3fIKq/tfdDG/Daa6/x2muveTuMZpnVaxb+Pv58kvFJ44W7DYOoAZDq3Dqbp4bBJkSEcMe8Dbz7vQ6DbQs0QSiX5ZVU8uqKfUwd0JWxvbq4VNdmbLy36z2GRQ1jYJdzF/e12/s15KbDhF+AC8Nmvalfv37061dPZ3sb0TmoM1O7T2XhvoVU1lY2XFjEeorI3gB5TjxxYA2D/eCecUzoE8ljH2/jyU+3U12rI5xaM00QymVPL9pBZY2Nx2YOcLnu6pzVHCo5xC39b2mg0IvQMR4Gt50tQhYuXMjChQu9HUazXZt0LSeqTrD0oBNzOobeAr6BsN75J6cOQf68+aMU7pyUyDvrDvLDN9aTf7KRZKS8RhOEcsmK3bks2HqYeyf3pndUmMv1305/m8jgyNOTs86T9b219tK4+62lHdqIF154gRdeeMHbYTTb6JjRxIfF89HejxovHBZlTWBM/S+UHXf6Hn6+Pvz6ioG8eNMwtmYVMftvq9ieo53XrZEmCOW08qpanvxsO72iQrnvkt4u10/NTWX90fXMGTSn/pnTK5+3hlCOuL2Z0aqm8BEfru17LRuObiCzKLPxCuPuh5py2PhPl+919fA4PrzH2qb1urlrePf7QzrzupXRBKGc9uwXO8k6Xs4frkkm0M/X5fqvpr1K58DO3JB0g+MCB1bB3q9g/M8g0PWnE+Ue1yVdR5BvEP/c7sQv/a4DoPcU+P4fUON6U1FyfDgLfjaRUT0jeOzjbfzs3S2U6KS6VkMThHLKV+lH+dfag/x0YqLLHdMA2/O3szpnNbcPup0QfweL7tlssOQJ6BgHY+9zQ8SqqSKCIrg+6Xo+z/yc7JLsxiuMux9OHrM2dWqCyLBA/nXHaP7n8n58sf0oV7y0irTsoiZdS7mXJgjVqMNF5fzvh2kkx4Xzv9P7N+kac7fOJTwwvP7O6e0fweEtcOmT1no/yqvmDJqDj/jw5vY3Gy/c+1LoNhRW/NFafbcJfHyE+y/pw/t3jaWm1sZ1c9fw6rf7qLVpk5M3aYJQDaqoruXe/2ymutbGS7cMJ8DP9X8yq3JWsTJ7JT8e9GNC/R1sJlR5EpY9BTHJMOQmN0Td8t555x3eeecdb4fhNtGh0VzT5xo+zfiUY6XHGi4sApc9DSeyYf2rzbpvSs8IFv9iElP6R/PHL3Zx7dw17D1W0qxrqqbTBKHqZbMZfjk/lbTsIv7vpmEkRjq/U9wp1bXV/On7P9GjYw9uG3ib40JLfwfFWTDjuTa7KF9CQgIJCQmNF2xD7ki+A2MMc7fObbxwr4uh7zT47i8ujWhypFNIAHNvHcFLtwznUEEpV7y0ipeXZ+icCS9omz+NyuOMMfzpq10s3naUx2cM4PJBri3lfco7O9/hwIkDPDr6UQJ8HQxbzfwWNvwDxt4LPcY1M2rvef/993n//fe9HYZbxYXF8YMBP+DjvR+zo2BH4xUue8ras/rbPzX73iLC7KGxfP3Li7lsYDTPfbWbq19ezZZDhc2+tnKeJgh1HmMMLyzZw2vfZnLr2O78dJJr24iekl2SzatbX+WShEuYGDfx/AIVJ2DBAxDR2+p7aMPmzp3L3LlO/KXdxtw99G46B3XmqbVPUWM7f/+Hs3QdYC3N/v3rkLXBLfePDAvk5R+OYO4PR5BXUsk1r6zhkQ+2kleik+tagiYIdRZjDH/6cjd/X57BLaMTeGr2YJd2iTulxlbD46sex0/8eGz0Y+cXsNXCRz+B4hy4em6r3k60PesY0JHHxzxOekE689LnNV5h6u+gQyx8eq/Vt+QmM5K78c0jk7nn4t58lprDpc+v4OXlGZRVNZK0VLNoglCnlVfV8vP3Unn12338cEx3fn91sssL8Z3y0uaX2JK7hV+P/TXdwhzsFbHkCdi7BK54HrqPaWbkypOm9ZjGtB7T+PuWv5Oam9pw4aCOcM1cOL4PFj1kLbzoJmGBfjw6oz9fPXgRY3pF8NxXu7nozyt4e80BKmtq3XYfdYYmCAXAoYIybnhtDYvSDvOr6f155urBTU4On2Z8ylvpb3FTv5u4opeD/R7W/A3WvQJj7tHlvNsAEeF3439Ht9BuPLTiIY6cPNJwhcSLYPJjsG0+rPqL2+PpFRXGGz8axUf3jqN3VCi/XZDO5OdW8MZ3mZys1CcKd9IE0c7ZbIZ/rT3A9L+u5GB+GW/cnsK9k3s3qVkJYOnBpfxuze8Y220svxr1q7NPGgPLnraeHgZeBdN+74avQLWEDgEd+Nulf6OipoJ7lt5Dfnl+wxUu+h9rU6FlT8GGNzwS08geEbx311j+dcdoukeE8MznOxn37DL+9OUuso6XeeSe7Y20tbVPUlJSzMaNG70dxgVhfWYBz36xi9SsIi5KiuKP1yYT26npk9Q+2vMRT697msGRg3ntstfOnvNQeRK++F9rl7gRt8OsF8HH9eU6Wqv8fOsXZmRkpJcj8awNRzdw/7L7iQ6J5pUpr5DQsYGhvTVVMP922PMFTHoYLvm1R/+fp2YV8frKfXy5/SgGmNgnkptHdWfqwK5NWhrmAuTyX32aINoZYwzf7c3nn6v28+2ePGI6BvHwtCSuHxnf5KeGsuoyntv4HB/u+ZDxseP5v8n/d/ZyGlnfwyd3w/H9cPH/Ws0PbWSfB3W+zcc28/PlPwfg8dGPMyNxRv3/dmoqYfEjsPlf0GcqXPMahHo2ieYUlfPBxizmb8jicHEFHQL9mDowmhmDY7goKYog/3abLDRBKMf255eyaOthPk3NYV9eKZFhAfx4QiJ3TEgkOKBpPzA1thoW7FvA3K1zOVZ6jDmD5/Dz4T/Hz8fPKnBsB6z4A+xcCOHd4ZpXoecEN35Vrce8efMAmDNnjlfjaCmHThzi0e8eZVv+Ni6Ov5iHRj5E704NrPC78S1Y/D/WMiqTHrb6n/yDPBpjrc2wKiOfz9MOs2THMYrKqgny92F0Yhcm9unC+N6R9I/pgJ9vu2lpb10JQkSmA38FfIE3jDF/POd8IPAvYCRQANxkjDnQ0DU1QTSuptbGvrxStucUs/lQIasz8jlQYLXJju4ZwU2jEpg1tFuTH7uzS7JZvH8xn2V8xqGSQyRHJvNIyiOMiB5hDVvd+xVsfR+y1kFgR2sxt7H3WSNcLlCTJ08GYMWKFV6NoyXV2mr5z87/8PfUv1NeU87EuIlc2/daxseOd7ykSt4e+PpJ2PMlhHSBITfD0JsgOtnjM+ira22s3VfAN7tyWZ2Rz95cawhusL8vyXHhDIkPZ0C3jiRGhdI7MozwkHqWo2/bWk+CEBFfYA9wGZANbABuMcbsqFPmPmCIMeYeEbkZuMYY0+BiPO01QRhjqKq1UVFt42RlDYWlVRSVVVNYVsXx0iqyC8vILiwnq7CMvcdOUlljLUsQGuDLuN5dmNgnkssHx9At3Lk+hmpbNSerTpJblkt2STbZJ7PZU7iHLblbyCrJAmBkeBK3RY7k0lp/JH835GyC4/Y9BKL6W+sqjZwDIRGe+Ja0Ku0xQZxyvOI483fP571d71FQUYCfjx9Do4YyIGIASZ2TiAmNITo0muiQaEL8QpADq6zZ87sWg60agiOsJ8vowRDVz9pNMKyr9eGhhRuPnahg7b4CUrOK2JpdRPrhE1TVnFnKo0toAD0jQ4npGERUh0C6dgwkukMQnUP9CQv0p0OQH2GBfnQM8icsyA/fJo74a2GtKkGMA35njLnc/v4xAGPMs3XKfGUvs1ZE/ICjQJRpIKiIHiFm2q/7cW4Bc96RU8cbft9Y/frrnF2joe+io3JyXh2DOed/n7H/59z6Uvf8KfaDIuCDgICPWEMUz7+/Oee9xQaUAifFUOngn1JErY1hFZWMqKjgsrIyYuuOPe8YD92GQM+JkHgxRA9qV/0M7TlBnFJtqyY1N5WV2SvZfGwzewr3UFF79uquPuJDsF+w9eETQGBtNT5VpfhWlSG1Vfhg8DHWP2dfDIIPItYHPr4gp540xP5vXs68hyb/mzNYg+xqjcFmM9gM2IzBGOt1fb+R5JwXde8u55x0OrImfAkOq5wT83t3b3H5yn6uh+K0OCCrzvts4NwZUafLGGNqRKQY6AKcNYZORO4C7gLo1D2Y2jpf+Vn/Q8669Jn/Ked+Vxx/l6TB/41n3cd+eyNnajmsaerEUedkvf+XzPnXE5HTX8OpjkAfAR8RROT0a5865c674bnxn/e1nfnKQ8WXDvgSKr6Eig9R4k+8TyjxgZ0ID+xkzXgOCLP/hRcDHaKhU48LuvlIOcffx59RMaMYFTMKsJqgDp88zNGyo+SW5ZJblktJVQnlNeWnPyprK7EZm/WL2FaNreoktupybLWVmJoq65ipxdhsYOwfp37+T//WrvOb0NExZwn4CnWeBs7+SbEShTl9G5v9hal7N3POH2HmzJ+SrkbU7D/d3fD3mScThNsYY14HXgeriemDu9pfE5NSbY2vjy8JHRMaHgqrWjVPJogcoO6/jHj7MUdlsu1NTOFYndVKtSmLFy/2dghKuZ0nhw5sAPqKSKKIBAA3AwvOKbMA+JH99fXANw31PyjVWoWEhBASogsOqguLx54g7H0KDwBfYQ1zfdMYky4iTwEbjTELgH8C74hIBnAcK4ko1ea88sorANx3n+6nrS4cOlFOKTfQUUyqDXC527rdTCFUSinlGk0QSimlHNIEoZRSyiFNEEoppRxqc53UIlIC7PZ2HK1EJOfMOm/H9Htxhn4vztDvxRlBxpjBrlRoEzOpz7HbGJPi7SBaAxHZqN8Li34vztDvxRn6vThDRFwe/qlNTEoppRzSBKGUUsqhtpggXvd2AK2Ifi/O0O/FGfq9OEO/F2e4/L1oc53USimlWkZbfIJQSinVAjRBKKWUcqhNJQgRmS4iu0UkQ0Qe9XY83iIiCSKyXER2iEi6iPzC2zF5k4j4isgWEVnk7Vi8TUQ6iciHIrJLRHbat/5td0TkIfvPxnYReVdEgrwdU0sSkTdFJFdEttc5FiEiX4vIXvvnzo1dp80kCBHxBV4GZgADgVtEZKB3o/KaGuBhY8xAYCxwfzv+XgD8Atjp7SBaib8CXxpj+gNDaYffFxGJA34OpNgnhvnS/rYSmAdMP+fYo8AyY0xfYJn9fYPaTIIARgMZxphMY0wV8B5wlZdj8gpjzBFjzGb76xKsXwJx3o3KO0QkHrgCeMPbsXibiIQDF2Hts4IxpsoYU+TVoLzHDwi271QZAhz2cjwtyhizEmuPnbquAt62v34buLqx67SlBBEHZNV5n007/aVYl4j0BIYD670cire8CPwv9j3k27lEIA94y97k9oaIhHo7qJZmjMkBngcOAUeAYmPMEu9G1SpEG2OO2F8fBaIbq9CWEoQ6h4iEAR8BDxpjTng7npYmIrOAXGPMJm/H0kr4ASOAucaY4UApTjQjXGjsbetXYSXMWCBURG71blSti31r50bnOLSlBJEDJNR5H28/1i6JiD9WcviPMeZjb8fjJROA2SJyAKvJ8VIR+bd3Q/KqbCDbGHPqafJDrITR3kwF9htj8owx1cDHwHgvx9QaHBORbgD2z7mNVWhLCWID0FdEEkUkAKvTaYGXY/IKERGsduadxpi/eDsebzHGPGaMiTfG9MT69/CNMabd/qVojDkKZIlIP/uhKcAOL4bkLYeAsSISYv9ZmUI77Kx3YAHwI/vrHwGfNVahzazmaoypEZEHgK+wRiW8aYxJ93JY3jIBuA3YJiKp9mOPG2MWey8k1Ur8DPiP/Y+oTODHXo6nxRlj1ovIh8BmrBF/W2hnS26IyLvAZCBSRLKB3wJ/BOaLyE+Ag8CNjV5Hl9pQSinlSFtqYlJKKdWCNEEopZRySBOEUkophzRBKKWUckgThFJKKYc0QSillHJIE4RSgIg8KCIhLXCfFBF5ydP3UcoddB6EUoB9uY4UY0y+t2NRqrXQJwjV7ohIqIh8LiJb7RvK/BZrUbflIrLcXmaaiKwVkc0i8oF9YURE5ICI/FlEtonI9yLSp4H73GC//lYRWWk/NvnUxkYislhEUu0fxSLyI/vmR8+JyAYRSRORuz3/HVHKMU0Qqj2aDhw2xgy1byjzItZ+AZcYYy4RkUjgCWCqMWYEsBH4ZZ36xcaYZODv9rr1+Q1wuTFmKDD73JPGmJnGmGHAqaUPPrW/LjbGjAJGAXeKSGIzvlalmkwThGqPtgGXicifRGSSMab4nPNjsXYtXG1f6+pHQI8659+t87mhLT1XA/NE5E6s9cPOY09G7wA/sMcxDbjdft/1QBegrwtfm1Ju02YW61PKXYwxe0RkBDATeEZElp1TRICvjTG31HeJel6fe597RGQM1o53m0Rk5Fk3sbbRfQ94yhhzau9gAX5mjPnK+a9IKc/QJwjV7ohILFBmjPk38BzWngklQAd7kXXAhFP9C/Y+i6Q6l7ipzue1DdyntzFmvTHmN1g7vSWcU+SPQJox5r06x74C7rXv94GIJLXHXeFU66BPEKo9SgaeExEbUA3ci9VU9KWIHLb3Q8wB3hWRQHudJ4A99tedRSQNqATqe8rAfo++WE8Fy4CtwMV1zj8CpNdZsv03WHtr9wQ22/cyyMOJvYOV8gQd5qqUC3Q4rGpPtIlJKaWUQ/oEoVQzicivgRvOOfyBMeb33ohHKXfRBKGUUsohbWJSSinlkCYIpZRSDmmCUEop5ZAmCKWUUg79f0aNdfeYjeo7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pyabc.visualization import plot_kde_1d\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for t in range(h.max_t + 1):\n",
    "    particles = h.get_distribution(m=0, t=t)\n",
    "    plot_kde_1d(\n",
    "        *particles,\n",
    "        \"step_size\",\n",
    "        label=f\"t={t}\",\n",
    "        ax=ax,\n",
    "        xmin=0,\n",
    "        xmax=10,\n",
    "        numx=300,\n",
    "    )\n",
    "ax.axvline(gt_step_size, color=\"k\", linestyle=\"dashed\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "That's it. You should be able to see how the distribution\n",
    "contracts around the true parameter."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}