My personal solution for probabilistic programming: binf

binf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistic programming with the Next Big Thing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... or rather YAPPP (yet another probabilistic programming package), this time my own: [https://bitbucket.org/simeon_carstens/binf/src/default/](`binf`), which stands for... you guessed it... Bayesian inference.  \n",
    "So what I do in `binf` is:\n",
    "- providing interfaces for components of posterior distributions: a posterior object, a likelihood, which includes a forward- and an error model, and priors. Varaibles / parameters of subordinate objects (like the error model) are propagated to the likelihood object, which then passes it on to the posterior object etc. You can then fix (\"condition\" on one or more variables) \n",
    "- providing samplers: there's a Gibbs, Metropolis-Hastings (in an example script) and HMC sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, let's jump right in with a basic example: fitting a polynomial. First, some data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEEVJREFUeJzt3X+MHOddx/HPh7PTbNPCpfWVxJcEJ1J1bcFCF1ZVGkNVJQVXaZUYA1L+KKRQZEX8ShE6sBVBJSTkFqOqRYCQSYpaEaUVqWtCaXEKToT4w6brnJNLenGbhLjN2TTXoktacWod98sfO+vYl93bGd/Ozj4375dkeW92due7s+ePn3lmnmccEQIApONHqi4AAFAMwQ0AiSG4ASAxBDcAJIbgBoDEENwAkBiCGwASQ3ADQGIIbgBIzIYy3nTTpk2xZcuWMt4aANalY8eOfTsiJvKsW0pwb9myRa1Wq4y3BoB1yfbJvOvSVQIAiSG4ASAxBDcAJIbgBoDEENwAkBiCGwASU8rlgABQFwdnF7Tv0AmdWlrW5vGGZrZPacf0ZKnbJLgB4CIdnF3QngNzWj5zVpK0sLSsPQfmJKnU8KarBAAu0r5DJ86FdsfymbPad+hEqdsluAHgIp1aWi60fFAIbgC4SJvHG4WWDwrBDQAXaWb7lBobxy5Y1tg4ppntU6Vul5OTAHCROicguaoEABKyY3qy9KBeia4SAEgMwQ0AiSG4ASAxBDcAJIbgBoDEENwAkBiCGwASQ3ADQGIIbgBIDMENAIkhuAEgMQQ3ACSG4AaAxBDcAJAYghsAEkNwA0BiCG4ASAzBDQCJyRXctn/f9pO2n7B9v+1Lyy4MANBd3+C2PSnp9yQ1I+KnJI1Jur3swgAA3eXtKtkgqWF7g6TXSjpVXkkAgNX0De6IWJD0F5K+Iem0pBcj4qGyCwMAdJenq+RySbdJulbSZkmX2X5/l/V22W7Zbi0uLg6+UgCApHxdJe+W9N8RsRgRZyQdkHTjypUiYn9ENCOiOTExMeg6AQCZPMH9DUk32H6tbUu6WdJ8uWUBAHrJ08d9VNIDkh6VNJe9Zn/JdQEAetiQZ6WI+LCkD5dcCwAgB0ZOAkBiCG4ASAzBDQCJIbgBIDEENwAkhuAGgMQQ3ACQGIIbABJDcANAYghuAEgMwQ0AiSG4ASAxuSaZAoD17ODsgvYdOqFTS8vaPN7QzPYp7ZierLqsnghuALV2cHZBew7MafnMWUnSwtKy9hyYk6SRDW+6SgDU2r5DJ86FdsfymbPad+hERRX1R3ADqLVTS8uFlo8CghtArW0ebxRaPgoIbgC1NrN9So2NYxcsa2wc08z2qYoq6o+TkwBqrXMCkqtKACAhO6YnRzqoVyK4R0Rq15ECqA7BPQJSvI4UQHU4OTkCUryOFEB1CO4RkOJ1pACqQ3CPgBSvIwVQHYJ7BKR4HSmA6nBycgSkeB0pgOoQ3CMitetIAVSHrhIASAzBDQCJIbgBIDEENwAkhpOT6wRznQD1QXAPSJXByVwnQL3QVTIAneBcWFpW6JXgPDi7MJTtM9cJUC+5gtv2uO0HbD9le972O8ouLCVVBydznQD1krfF/QlJ/xoRb5H005LmyyspPVUHJ3OdAPXSN7ht/6ikd0q6V5Ii4gcRsVR2YSmpOjiZ6wSolzwt7uskLUr6e9uztu+xfVnJdSWl6uDcMT2pvTu3anK8IUuaHG9o786tnJgE1ilHxOor2E1JRyRti4ijtj8h6aWI+OMV6+2StEuSrrnmmp85efJkSSWPJi7HA7AWto9FRDPXujmC+wpJRyJiS/bzz0naHRHv7fWaZrMZrVYrf8UAUHNFgrtvV0lE/I+kb9ruHPffLOmra6gPALAGeQfg/K6k+2xfIulZSb9eXkkAgNXkCu6IOC4pVxMeAFAuRk4CQGIIbgBIDMENAIkhuAEgMQQ3ACSG4AaAxBDcAJAYghsAEkNwA0BiCG4ASAzBDQCJIbgBIDEENwAkJu+0rljnuIMPkA6CGzo4u6A9B+a0fOasJGlhaVl7DsxJEuENjCC6SqB9h06cC+2O5TNnte/QiYoqArAaghs6tbRcaDmAahHc0ObxRqHlAKpFcEMz26fU2Dh2wbLGxjHNbJ/q8QoAVeLkJM6dgOSqEiANBDcktcOboAbSQFcJACSG4AaAxBDcAJAY+rgBJK9uUzYQ3ACSVscpG+gqAZC0Ok7ZQHADSFodp2wguAEkrY5TNhDcAJJWxykbODkJIGl1nLKB4AaQvLpN2UBXCQAkhuAGgMQQ3ACQmNzBbXvM9qztL5RZEABgdUVa3HdJmi+rEABAPrmC2/ZVkt4r6Z5yywEA9JO3xf1xSX8o6Ycl1gIAyKFvcNt+n6QXIuJYn/V22W7Zbi0uLg6sQADAhfK0uLdJutX2c5I+I+km2/+wcqWI2B8RzYhoTkxMDLhMAEBH3+COiD0RcVVEbJF0u6TDEfH+0isDAHTFddwAkJhCc5VExCOSHimlEgBALrS4ASAxzA6YqdvNRgGki+BWPW82CiBdBLdWv9kowT0cHPEA+RHcqufNRkcJRzxAMevm5OTB2QVt+8hhXbv7X7TtI4d1cHYh92vreLPRUbLaEQ+AV1sXwd1psS0sLSv0Sostb3jX8Wajo4QjHqCYdRHca22x7Zie1N6dWzU53pAlTY43tHfnVg7TC+CIBxieddHHPYgWW91uNjpIa+2jntk+dcHrJY54gNWsixY3LbZqccQDDNe6aHHTYqsWRzzAcK2LFjcttmpxxAMM17pocUu02KrEEQ8wXOsmuFGdzn+YjHwEhoPgxkBwxAMMz7ro4waAOiG4ASAxBDcAJIbgBoDEENwAkBiCGwASQ3ADQGIIbgBIDMENAIlh5CTWBW42nDa+v2IIbiSPmw2nje+vOLpKkDxuNpw2vr/iCG4kj5sNp43vrziCG8njRg5p4/srjuBG8ma2T6mxceyCZdzIIR18f8VxchLJ40YOaeP7K84RMfA3bTab0Wq1Bv6+ALBe2T4WEc0869JVAgCJIbgBIDEj08fNyCkAyKdvi9v21bYftj1v+0nbdw26iM7IqYWlZYVeGTl1cHZh0JsCgOTl6Sp5WdIfRMRbJd0g6bdtv22QRTByCgDy6xvcEXE6Ih7NHn9X0rykgfZhMHIKAPIrdHLS9hZJ05KODrIIRk4BQH65g9v26yR9TtKHIuKlLs/vst2y3VpcXCxUBCOnACC/XMFte6PaoX1fRBzotk5E7I+IZkQ0JyYmChWxY3pSe3du1eR4Q5Y0Od7Q3p1buaoEALroezmgbUu6V9J8RHysrEJ2TE8S1ACQQ54W9zZJvyrpJtvHsz+3lFwXAKCHvi3uiPhPSR5CLQCAHBjyDgCJIbgBIDEENwAkZmQmmQKqxCRnSAnBjdrrTHLWmS+nM8mZJMIbI4muEtQek5whNQQ3ao9JzpAaghu1xyRnSA3BjdpjkjOkhpOTqL3OCcgqryrhqhYUQXADqnaSM65qQVF0lQAV46oWFEVwAxXjqhYURXADFeOqFhRFcAMV46oWFMXJSaBio3BVC9JCcAMjgFv3oQi6SgAgMQQ3ACSGrhJgAKoe+Vj37dcNwQ2sUdUjH+u+/TqiqwRYo6pHPtZ9+3VEcANrVPXIx7pvv44IbmCNqh75WPft1xHBDaxR1SMf6779OuLkJLBGVY98rPv268gRMfA3bTab0Wq1Bv6+ALBe2T4WEc0869JVAgCJIbgBIDEENwAkhpOTABiynhiCG6g5hqynh64SoOYYsp4eghuoOYasp4fgBmqOIevpIbiBmmPIenpyBbft99g+Yftp27vLLgrA8OyYntTenVs1Od6QJU2ON7R351ZOTI6wvleV2B6T9NeSfl7S85K+YvvBiPhq2cUBGA5uVpyWPC3ut0t6OiKejYgfSPqMpNvKLQsA0Eue4J6U9M3zfn4+W3YB27tst2y3FhcXB1UfAGCFPMHtLsteNaVgROyPiGZENCcmJtZeGQCgqzzB/bykq8/7+SpJp8opBwDQT57g/oqkN9u+1vYlkm6X9GC5ZQEAesl1IwXbt0j6uKQxSZ+MiD/rs/6ipJMXWdMmSd++yNeWibqKoa5iqKuYUaxrrTX9RETk6mcu5Q44a2G7lfcuEMNEXcVQVzHUVcwo1jXMmhg5CQCJIbgBIDGjGNz7qy6gB+oqhrqKoa5iRrGuodU0cn3cAIDVjWKLGwCwisqD2/Y+20/Zftz2522P91hvqDMU2v4V20/a/qHtnmeKbT9ne872cdutEapr2PvrDba/bPvr2d+X91jvbLavjtsubTxAv89v+zW2P5s9f9T2lrJqKVjXB2wvnrePfnMINX3S9gu2n+jxvG3/ZVbz47avL7umnHW9y/aL5+2rPxlCTVfbftj2fPbv8K4u65S/vyKi0j+SfkHShuzxRyV9tMs6Y5KekXSdpEskPSbpbSXX9VZJU5IekdRcZb3nJG0a4v7qW1dF++vPJe3OHu/u9j1mz31vCPuo7+eX9FuS/jZ7fLukz45IXR+Q9FfD+n3KtvlOSddLeqLH87dI+pLa01/cIOnoiNT1LklfGPK+ulLS9dnj10v6WpfvsPT9VXmLOyIeioiXsx+PqD2kfqWhz1AYEfMRMXI33ctZVxUzOt4m6VPZ409J2lHy9laT5/OfX+8Dkm623W1enmHXNXQR8R+S/neVVW6T9OloOyJp3PaVI1DX0EXE6Yh4NHv8XUnzevWke6Xvr8qDe4XfUPt/qpVyzVBYkZD0kO1jtndVXUymiv314xFxWmr/ckt6U4/1Ls1mkTxiu6xwz/P5z62TNRxelPTGkuopUpck/VJ2iP2A7au7PD9so/zv7x22H7P9Jds/OcwNZ91r05KOrniq9P3V90YKg2D73yRd0eWpuyPin7J17pb0sqT7ur1Fl2VrvhwmT105bIuIU7bfJOnLtp/KWgpV1jX0/VXgba7J9td1kg7bnouIZ9Za2wp5Pn8p+6iPPNv8Z0n3R8T3bd+p9lHBTSXX1U8V+yqPR9UeJv69bFqOg5LePIwN236dpM9J+lBEvLTy6S4vGej+GkpwR8S7V3ve9h2S3ifp5sg6iVYoZYbCfnXlfI9T2d8v2P682ofDawruAdQ19P1l+1u2r4yI09lh4Qs93qOzv561/YjaLZZBB3eez99Z53nbGyT9mMo/LO9bV0R857wf/07t8z5VG8kZQs8PzIj4ou2/sb0pIkqdw8T2RrVD+76IONBlldL3V+VdJbbfI+mPJN0aEf/XY7WRnKHQ9mW2X995rPaJ1q5nwIesiv31oKQ7ssd3SHrVkYHty22/Jnu8SdI2SWXcAi/P5z+/3l+WdLhHo2Goda3oC71V7T7Uqj0o6deyqyVukPRip1usSrav6JyXsP12tfPsO6u/as3btKR7Jc1HxMd6rFb+/hrmGdkeZ2mfVrs/6Hj2p3Omf7OkL644U/s1tVtndw+hrl9U+3/O70v6lqRDK+tS++qAx7I/T45KXRXtrzdK+ndJX8/+fkO2vCnpnuzxjZLmsv01J+mDJdbzqs8v6U/VbiBI0qWS/jH7/fsvSdeVvY9y1rU3+116TNLDkt4yhJrul3Ra0pnsd+uDku6UdGf2vNW+7+wz2ffW8yqrIdf1O+ftqyOSbhxCTT+rdrfH4+dl1i3D3l+MnASAxFTeVQIAKIbgBoDEENwAkBiCGwASQ3ADQGIIbgBIDMENAIkhuAEgMf8PBc+UXga5mPYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_data_points = 20\n",
    "real_coeffs = np.array([2.0, -4.0, 1.0, 1.5])\n",
    "real_precision = 2.5\n",
    "polynomial = np.polynomial.polynomial.polyval\n",
    "xses = np.linspace(-2, 2, n_data_points)\n",
    "ys = np.random.normal(loc=polynomial(xses, real_coeffs), \n",
    "                      scale=1.0 / np.sqrt(real_precision))\n",
    "\n",
    "plt.scatter(xses, ys)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we build up our likelihood from two components, a forward- and an error model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from csb.statistics.pdf.parameterized import Parameter as ScalarParameter\n",
    "\n",
    "from binf import ArrayParameter\n",
    "from binf.model.forwardmodels import AbstractForwardModel\n",
    "from binf.model.errormodels import AbstractErrorModel\n",
    "from binf.pdf import ParameterNotFoundError\n",
    "\n",
    "class ForwardModel(AbstractForwardModel):\n",
    "\n",
    "    def __init__(self, xses, polynomial):\n",
    "\n",
    "        super(ForwardModel, self).__init__('polynomial')\n",
    "        \n",
    "        self.xses = xses\n",
    "        self.polynomial = polynomial\n",
    "        \n",
    "        self._register_variable('coefficients', differentiable=True)\n",
    "        self.update_var_param_types(coefficients=ArrayParameter)\n",
    "        self._set_original_variables()\n",
    "\n",
    "    def _evaluate(self, coefficients):\n",
    "\n",
    "        return self.polynomial(self.xses, coefficients)\n",
    "\n",
    "    def _evaluate_jacobi_matrix(self, coefficients):\n",
    "\n",
    "        return np.vstack([self.xses ** i for i in range(len(coefficients))])\n",
    "\n",
    "    def clone(self):\n",
    "\n",
    "        copy = self.__class__(self.xses, self.polynomial)\n",
    "        self._set_parameters(copy)\n",
    "\n",
    "        return copy\n",
    "\n",
    "\n",
    "class GaussianErrorModel(AbstractErrorModel):\n",
    "\n",
    "    def __init__(self, ys):\n",
    "\n",
    "        super(GaussianErrorModel, self).__init__('error_model')\n",
    "\n",
    "        self.ys = ys\n",
    "\n",
    "        self._register_variable('mock_data')\n",
    "        self._register_variable('precision')\n",
    "        self.update_var_param_types(mock_data=ArrayParameter,\n",
    "                                    precision=ScalarParameter)\n",
    "        self._set_original_variables()\n",
    "\n",
    "    def _evaluate_log_prob(self, mock_data, precision):\n",
    "\n",
    "        logZ = len(self.ys) * 0.5 * np.log(precision)\n",
    "        return -0.5 * np.sum((mock_data - self.ys) ** 2) * precision + logZ\n",
    "\n",
    "    def _evaluate_gradient(self, mock_data, precision):\n",
    "\n",
    "        return (mock_data - self.ys) * precision\n",
    "\n",
    "    def clone(self):\n",
    "\n",
    "        copy = self.__class__(self.ys)\n",
    "        copy.set_fixed_variables_from_pdf(self)\n",
    "\n",
    "        return copy\n",
    "\n",
    "def make_likelihood(xses, ys, polynomial):\n",
    "    \n",
    "    from binf.pdf.likelihoods import Likelihood\n",
    "    from binf.example.likelihood import ForwardModel, GaussianErrorModel\n",
    "\n",
    "    FWM = ForwardModel(xses, polynomial)\n",
    "    EM = GaussianErrorModel(ys)\n",
    "    L = Likelihood('points', FWM, EM)\n",
    "\n",
    "    return L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make a likelihood from those:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from binf.pdf.likelihoods import Likelihood\n",
    "likelihood = Likelihood('points', ForwardModel(xses, polynomial), GaussianErrorModel(ys))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need priors for all variables, that is, for `coefficients` and for `precision`. We choose a Gaussian distribution centered around 0 for the former and a Gamma distribution (the conjugated prior for the precision of a Gaussian likelihood) for the latter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from csb.statistics.pdf.parameterized import Parameter as ScalarParameter\n",
    "from binf.pdf.priors import AbstractPrior\n",
    "\n",
    "class GammaPrior(AbstractPrior):\n",
    "\n",
    "    def __init__(self, shape, rate):\n",
    "\n",
    "        super(GammaPrior, self).__init__('precision_prior')\n",
    "\n",
    "        self.shape = shape\n",
    "        self.rate = rate\n",
    "        \n",
    "        self._register_variable('precision')\n",
    "        self.update_var_param_types(precision=ScalarParameter)\n",
    "        self._set_original_variables()\n",
    "\n",
    "    def _evaluate_log_prob(self, precision):\n",
    "\n",
    "        return (self.shape - 1.0) * np.log(precision) - precision * self.rate\n",
    "\n",
    "    def clone(self):\n",
    "\n",
    "        copy = self.__class__(self.shape, self.shape)\n",
    "        copy.set_fixed_variables_from_pdf(self)\n",
    "\n",
    "        return copy\n",
    "            \n",
    "\n",
    "class GaussianPrior(AbstractPrior):\n",
    "\n",
    "    def __init__(self, means, variances):\n",
    "\n",
    "        super(GaussianPrior, self).__init__('coefficients_prior')\n",
    "        \n",
    "        self._register('means')\n",
    "        self._register('variances')\n",
    "        self['means'] = ArrayParameter(means, 'means')\n",
    "        self['variances'] = ArrayParameter(variances, 'variances')\n",
    "        self._register_variable('coefficients')\n",
    "        self.update_var_param_types(coefficients=ArrayParameter)\n",
    "        self._set_original_variables()\n",
    "\n",
    "    def _evaluate_log_prob(self, coefficients):\n",
    "\n",
    "        means = self['means'].value\n",
    "        variances = self['variances'].value\n",
    "\n",
    "        return -0.5 * np.sum((coefficients - means) ** 2 / variances)\n",
    "    \n",
    "    def _evaluate_gradient(self, **variables):\n",
    "\n",
    "        coeff = variables[self.variable_name]\n",
    "        \n",
    "        return (coeff - self['mu'].value) / self['sigma'].value ** 2\n",
    "\n",
    "    def clone(self):\n",
    "\n",
    "        return self.__class__(self['means'].value, self['variances'].value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's instantiate those and store them in a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "PP = GammaPrior(1.0, 0.2)\n",
    "CP = GaussianPrior(means=np.array([0.0, 0.0, 0.0, 0.0]), variances=np.ones(4) * 5)\n",
    "priors = {PP.name: PP, CP.name: CP}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have everything we need to create a posterior object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from binf.pdf.posteriors import Posterior\n",
    "posterior = Posterior(likelihoods={likelihood.name: likelihood},\n",
    "                      priors=priors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What can we do with the posterior object? We can, for example, look at all variables and parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set(['precision', 'coefficients'])\n",
      "('means', 'variances')\n"
     ]
    }
   ],
   "source": [
    "print(posterior.variables)\n",
    "print(posterior.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All variables and parameters got inherited from posterior components, in our case, the likelihood.  \n",
    "Now let's make a posterior distribution which is conditioned on a certain value for `precision`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set(['coefficients'])\n",
      "('precision', 'means', 'variances')\n"
     ]
    }
   ],
   "source": [
    "conditioned_posterior = posterior.conditional_factory(precision=1.0)\n",
    "print(conditioned_posterior.variables)\n",
    "print(conditioned_posterior.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conditioned posterior has only the `coefficients` left as variables, while `precision` has been turned into a parameter. This functionality is of paramount importance for Gibbs sampling.  \n",
    "We can also clone a posterior (or any PDF in `binf`, for that matter), which will preserve the conditioning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set(['coefficients'])\n",
      "('precision', 'means', 'variances')\n"
     ]
    }
   ],
   "source": [
    "clone = conditioned_posterior.clone()\n",
    "print(clone.variables)\n",
    "print(clone.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To sample from this posterior, we will use the in-built implemented Gibbs sampler, which draws samples by cycling over the conditional distributions for `coefficients` and `precision`. This means we need samplers to draw samples from these conditional distributions. We use a Metropolis-Hastings sampler for the coefficients and a wrapper around the `numpy` built-in sampler for Gamma-distributed random variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "\n",
    "RWMCSampleStats = namedtuple('RWMCSampleStats', 'acceptance_rate')\n",
    "\n",
    "\n",
    "class GammaSampler(object):\n",
    "    \n",
    "    def __init__(self, pdf, state):\n",
    "\n",
    "        self.pdf = pdf\n",
    "        self.state = state\n",
    "\n",
    "    def _get_prior(self):\n",
    "\n",
    "        prior = filter(lambda p: 'precision' in p.variables,\n",
    "                       self.pdf.priors.values())[0]\n",
    "        \n",
    "        if not isinstance(prior, GammaPrior):\n",
    "            msg = 'Prior for precision is not a Gamma distribution'\n",
    "            raise NotImplementedError(msg)\n",
    "        \n",
    "        return prior\n",
    "\n",
    "    def _calculate_shape(self):\n",
    "\n",
    "        prior = self._get_prior()\n",
    "        n_data_points = len(self.pdf.likelihoods['points'].error_model.ys)\n",
    "        \n",
    "        return 0.5 * n_data_points + prior.shape - 1\n",
    "        \n",
    "    def _calculate_rate(self):\n",
    "\n",
    "        args = dict(coefficients=self.pdf['coefficients'].value,\n",
    "                    precision=1.0)\n",
    "        r1 = -self.pdf.likelihoods['points'].log_prob(**args)\n",
    "        r2 = self._get_prior().rate\n",
    "        \n",
    "        return r1 + r2\n",
    "    \n",
    "    def sample(self, state=42):\n",
    "\n",
    "        rate = self._calculate_rate()        \n",
    "        shape = self._calculate_shape()\n",
    "        sample = np.random.gamma(shape) / rate\n",
    "\n",
    "        self.state = sample\n",
    "        \n",
    "        return self.state\n",
    "\n",
    "\n",
    "class RWMCSampler(object):\n",
    "\n",
    "    def __init__(self, pdf, state, stepsize):\n",
    "\n",
    "        self.pdf = pdf\n",
    "        self.state = state \n",
    "        self.stepsize = stepsize\n",
    "        self._n_moves = 0\n",
    "        self._n_accepted_moves = 0\n",
    "\n",
    "    @property\n",
    "    def last_draw_stats(self):\n",
    "        \n",
    "        return {'coefficients': RWMCSampleStats(self.acceptance_rate)}\n",
    "\n",
    "    @property\n",
    "    def acceptance_rate(self):\n",
    "        if self._n_moves > 0:\n",
    "            return self._n_accepted_moves / float(self._n_moves) \n",
    "        else:\n",
    "            return 0.0\n",
    "\n",
    "    def sample(self):\n",
    "\n",
    "        E_old = -self.pdf.log_prob(coefficients=self.state)\n",
    "        change = np.random.uniform(low=-self.stepsize, high=self.stepsize,\n",
    "                                   size=len(self.state))\n",
    "        proposal = self.state + change\n",
    "        E_new = -self.pdf.log_prob(coefficients=proposal)\n",
    "\n",
    "        accepted = np.random.random() < np.exp(-(E_new - E_old))\n",
    "\n",
    "        if accepted:\n",
    "            self.state = proposal\n",
    "            self._n_accepted_moves += 1\n",
    "\n",
    "        self._n_moves += 1\n",
    "\n",
    "        return self.state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can instantiate a Gibbs sampler, we need an initial state for the Markov chain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from binf.samplers import BinfState\n",
    "start = BinfState(dict(coefficients=np.ones(4), precision=1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set up the samplers for the conditional distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mh_stepsize = .1\n",
    "coeffs = start.variables['coefficients']\n",
    "precision = start.variables['precision']\n",
    "coefficients_pdf = posterior.conditional_factory(precision=precision)\n",
    "coefficients_sampler = RWMCSampler(coefficients_pdf,\n",
    "                                   coeffs,\n",
    "                                   mh_stepsize)\n",
    "\n",
    "precision_pdf = posterior.conditional_factory(coefficients=coeffs)\n",
    "precision_sampler = GammaSampler(precision_pdf, precision)\n",
    "\n",
    "subsamplers = {'coefficients': coefficients_sampler,\n",
    "               'precision': precision_sampler}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, here's the Gibbs sampler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from binf.samplers.gibbs import GibbsSampler\n",
    "# the naming here is funny only for German-speaking people\n",
    "gips = GibbsSampler(posterior, start, subsamplers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that thing now can actually sample from the posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'precision': 0.024801358393265305, 'coefficients': array([0.97042051, 1.04955161, 1.0045029 , 0.90398355])}\n"
     ]
    }
   ],
   "source": [
    "sample = gips.sample()\n",
    "print sample.variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get a few samples..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#### Gibbs sampling step 5000 ####\n",
      "RWMC acceptance rate: 0.513076520722\n",
      "#### Gibbs sampling step 10000 ####\n",
      "RWMC acceptance rate: 0.512471072255\n",
      "#### Gibbs sampling step 15000 ####\n",
      "RWMC acceptance rate: 0.511839684275\n",
      "#### Gibbs sampling step 20000 ####\n",
      "RWMC acceptance rate: 0.511512212195\n",
      "#### Gibbs sampling step 25000 ####\n",
      "RWMC acceptance rate: 0.510870332463\n"
     ]
    }
   ],
   "source": [
    "from copy import deepcopy\n",
    "samples = []\n",
    "for i in range(30000):\n",
    "    samples.append(deepcopy(gips.sample()))\n",
    "    if i % 5000 == 0 and i > 0:\n",
    "        print \"#### Gibbs sampling step {} ####\".format(i) \n",
    "        print 'RWMC acceptance rate: {}'.format(gips.last_draw_stats['coefficients'].acceptance_rate)\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and plot the data again together with a bunch of samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvWmUZGl61/e7N+LGjRs39oxcIveq6qruanXPounRMJpFmpHMCGkQAwIMWICNdcbLMcjLGYw+yDoHY+PDGLDkgz8MyBI2WIYjBgkLpNGCpGFAGmbVLF1d3V1ZVblF5BL7dvfXH7LfdzJreq3KyvX9nZOnqzKzMm52Zvzjuf/3ef6PIYRAo9FoNOcf87QvQKPRaDTHgxZ0jUajuSBoQddoNJoLghZ0jUajuSBoQddoNJoLghZ0jUajuSBoQddoNJoLghZ0jUajuSBoQddoNJoLQvokH6xWq4nV1dWTfEiNRqM593zpS1/aF0JMv9Hnnaigr66u8sUvfvEkH1Kj0WjOPYZh3H8zn6ctF41Go7kgaEHXaDSaC4IWdI1Go7kgaEHXaDSaC4IWdI1Go7kgaEHXaDSaC4IWdI1Go7kgaEHXaDSax4jv+8RxfCKPpQVdo9FoHhNCCC3oGo1GcxGI45jhcKgFXaPRaM47QRAghMCyrBN5PC3oGo1G85gIgoDRaEQURSfyeFrQNRqN5jHheR5hGGrLRaPRaM4zSZIwHo9JkoR0+mSCbbWgazQazWMgiiI8zyOOY4QQJ/KYWtA1Go3mMRAEAZ7nkUqltKBrNBrNeWYymRCGIYC2XDQajea8IoRQ/vnOzg7dbvdEHlcLukaj0RwzSZIwmUzwPI+dnR3G4/GJPO4b3gcYhvF/Ah8FdoUQz7zyvirwT4BV4B7wp4UQncd3mRqNRnN+CIIA3/f5ja9t8Au/9WUmv+tz7eoqn/jIk3zsnQuP7XHfTIX+88APPPC+vwb8lhDiOvBbr/xdo9FoNBwI+q99dZ2f/e3n6foxKSfHVnfCT3z66/zSV7Ye2+O+oaALIT4LtB949x8D/uErf/6HwMeO+bo0Go3m3DIcDvm/f+8e3mQC6SymZQMwCWM++Znbj+1xH9ZDnxVCNABe+e/M8V2SRqPRnF/kQFGjO0BEPknkI4KJ+vh2d/I6//rReOyHooZhfNwwjC8ahvHFvb29x/1wGo1Gc6pEUcRkMqFmxSSxR+yPiMdD9fH5svPYHvthBX3HMIw6wCv/3X2tTxRCfEoI8ZwQ4rnp6emHfDiNRqM5H/i+j+/7fPSZGVJxhGm7kCsA4FgpPvGRJx/bYz+soP8L4C++8ue/CPzy8VyORqPRnG+GwyFJkvDMTIYfeKbOVN4hZaRYKDv8zT/x7GPtcnkzbYu/AHwvUDMMYxP4KeB/Af6pYRj/KbAO/KnHdoUajUZzTpADRUEQMB6PeXahxB95zyJ/9I9+H1NTU4/98d9Q0IUQf/Y1PvR9x3wtGo1Gc66J45jxeKyGigzDoFAokM1mT+Tx9aSoRqPRHBNBEKilFsPhENu2qVarWtA1Go3mvDEcDomiSO0RTafT2LZNEAQn8vha0DUajeaYGAwGhGHIYDDANE2y2ayyYE4CLegajUZzDMhArvF4jOd5qjq/ffs229vbJ3INWtA1Go3mGJD955PJhH6/TzqdJo5jPM/Dtu0TuQYt6BqNRnMMSP98MBgQRRGZTEZ1u1QqlRO5Bi3oGo1GcwwMh8Mj/nkul2M8HmMYxoldgxZ0jUajeUTkQJH00NPpNJZlMRqNsCyLKIpO5Dq0oGs0Gs0jEoYhnucxHA4ZDAZkMhmCICCKIhzHIZVKnch1aEHXaDSaR0T2ncv/2raN53mMx2Oq1SqO8/gSFg+jBV2j0WgekX6/rw5ETdPEcRzV3WLbNkKIE7kOLegajUbzCAghGA6HatzfNE0sy1Kj/47jnNjBqBZ0jUajeQRk//loNGI0GpHJZPA8jyAIyOfzqh/9JNCCrtFoNI+AzD+X/ee2beP7PmEYUigUcF0Xy7JO5Fq0oGs0Gs0jMBgM8H2ffr+PYRiq/zydTuM4Dq7r6rRFjUajOQ+MRiPVsphKpdSSaMdxKBQKpFIpwjA8kWvRgq7RaDQPied5+L7PYDBgPB5j2zZRFOH7PpVKBdd1KRaLpNNvuEvoWNCCrtFoNA+JzG/p9/skSYJt24zHYwAKhQL5fJ58Pq8HizQajeas0+/38X2fTqeDEIJMJoPv+2SzWRzHIZvNksvlTux6tKBrNBrNQyDzW+REqNxM5HkeuVwO13VV2+JJoQVdo9FoHoIgCPB9n16vh+/7WJZFHMf4vq/sFnkoelJoQddoNJqHoN/vE8cxvV6POI6xLAvf9wHI5/PkcjkKhQKmeXIyqwVdo9FoHoJ+v6+2EwEYhsFkMlF2i+M4J7apSKIFXaPRaN4istd8MpmozPMwDAnDkFKppOyWk/TP4REF3TCM/8YwjG8ahvENwzB+wTCMkxmH0mg0mlNkMpkQBAGdTocgCMhkMoRhSBzH5PN5XNelXC6fqH8OjyDohmEsAH8FeE4I8QyQAv7McV2YRqPRnFX6/T5hGNLtdomiSAVwmaZJPp/Htu0TbVeUPKrlkgYcwzDSQA7YfvRL0mg0mrPNYDBQdotpmiRJotIV5XRoJpM58et6aEEXQmwB/yuwDjSAnhDi14/rwjQajeYsEkURnuep/HO5MzQMQ4rFIoVCgVKpdOL+OTya5VIB/hhwBZgHXMMwfvRVPu/jhmF80TCML+7t7T38lWo0Gs0ZYDweEwQB7XabOI5JpVJqf2g+nyebzVIoFE5sqcVhHsVy+X7grhBiTwgRAp8GvvvBTxJCfEoI8ZwQ4rnp6elHeDiNRqM5ffr9Pp7n0e12laAnSaLaFXO53Im3K0oeRdDXgT9kGEbOOHgp+j7g1vFclkaj0Zw9kiRR6+Ymk4k6DPV9n1wuR7lcplQqndhCiwd5FA/988AvAl8Gvv7K1/rUMV2XRqPRnDmkfy7jcmW7YpIkqjovlUon3q4oeSTXXgjxU8BPHdO1aDQazZlmNBoRhiF7e3skSYIQgiRJSKVS1Go1JeqnhZ4U1Wg0mjeBEIJ+v6+WQcv3RVFENptV6Yqn5Z+DFnSNRqN5U8RxzGQyYTgc4nke6XQa3/eJ41gtgz7pdMUH0YKu0Wg0bwLf91V3SxiGGIaBaZrEcUy1WlUDRafRrig5+c53jUajOYcMBgM8z6PVahHHMUmSKLulUCiQy+XIZk83zkpX6BqNRvMGJEnCaDSi1+vheR5wYMGEYYjjOBQKBYrF4qm1K0q0oGs0Gs0bEIahyj4PguDIxyqVivLPT3KZxauhBV2j0WjegMlkgud57O3tEUURSZKoKdFyuUwulyOfz5/2ZWpB12g0mtdDCMFwOKTX6zEej0mSRA0TOY5DqVSiWCyearuiRAu6RqPRvA5RFDEejxkMBkRRhBBCtSbm83mKxSKu65663QJa0DUajeZ1ke2KrVZLVebSbqnVajiOcybsFtBti2cWIYT682n2tWo0lxkhBOPxmF6vx2AwUIssAEzTpFqtks/nT3Xc/zBa0E8RIQRxHKueVvn2ahiGoQYZ5FsqlToTt3kazUUljmNGo5GKzI2iSFXnchl0Pp8/lWUWr8bZuIpLhDxQkSflEinS6XRaVeSGYahKXYYAyWEGiWEYpNNp0uk0qVRKV/MazTEShiG+77O/v6/6zuVzbGpqCtd1cV33lK/yW2hBPwFkgE8QBErE0+k0lmWRSqUeKvtB+nhxHKv1V4e/7lmpGDSa84zneQyHQ2W3hGGIaZpYlqX6z7WgXxKEEARBQBiG6mQ8m80eqcLlL4mswIUQR/xzOKjepeUibRb5JifTDgt7FEXqY5Zl6apdo3kIZBjX/v4+k8mEIAiI4xjLslS74lkY9z+MFvTHRBAEBEGAEIJ0Ok0mkyGVSqlqXXpxD4o3oERdVuHy70KIIz56JpM5YrfYto1t2+puwPd9giDAsiwymYwWdo3mLRCGIaPRiE6nQxzHRyZEp6amlH9+ltCCfszEcYzneSr0PpvNqkS2yWSi/G9ZcR8mDENlo8Rx/G1fWx6iykoeUL57JpNRb+l0WlUN8oUlCAL1cS3sGs0b4/s+4/FYLYP2fV81I8juFi3oFxRprwRBgGmaOI5DOp1WGRBJkqjqWn6+POCUAi0tFfkicLizRf6bw2/ykEb2yXqeh2maZLNZ5c9LcZcVRhiGStg1Gs2rkyQJvu+zu7ur7JYoipTFMjMzg+M4Z6ZdUaIF/RhIkkSJthRLmc52WMiliAPf5pkftlg8z3tVK0ZaLVL4U6mUOgAVQij/3Pd9oigik8kQx7F6fMuyVKURhiHZbPZUw/g1mrNKGIaMx2O63S5CCHzfxzAMkiRhenqaXC535sQctKA/MmEY4nkehmGQy+UwTVP1qz4o5FKMpbct7ZfDh5yyS+WwJfNg+6J8k4ephwVb2juyGpe+urxW+fWSJGE8HmNZFrZtaxtGozmEtFtku6Lv+8BBF1mlUqFQKJw5uwW0oD8SUphTqRSO4xBFEZPJBOCIkKfTaVWx+76v3ic7VGRlffig9EEbRgq2PFw9PFwkP0++aKRSKTKZDL7vM5lMGI/HajzZNE31QiArDxnSr1sdNRrUNGi73WYymajnpyx+ZmdnyWQyWtAvEvKAU1osh6tyQB2KAmqpbBRFpFIpwjBkOBwqET9czR9Gvig82OUivXH52NI6kWIuhdlxHBzHwfM8tdw2n89TKpW+zaKRdpGu1jWXnSiKVFSubGYIw5B0Ok2tVsN1XXK53Jl8nmhBf4sIIZhMJsRxrIRURmqapqk8c7lAttvtMh6P1b+R/rgUcemDS/GW3S1SvA3DUB9Pp9Mq7W08HjMajVSlns1m1dthf92yLOX3yQjQwWBApVKhVCoRhuGRThh58KMjBTSXFWm37O7uqsPRdDqNaZqqXbFQKJz2Zb4qWtDfAlJIZQ4ywHg8/jbLQwjBzs4OOzs7asOJFGdptTiOo27tZIUsBxdka6KsAGT1LTtWbNtWh6GHN6nI6ly2U8lqWz5esVikWCzSarXY29uj1+sxPT2N67rqBWg0GhGGIfl8/tTXaWk0J81hu2U4HKoOMtu2cRyHubk5bNs+k3YLPKKgG4ZRBv4B8AwggL8khPi947iws4YUcyEEjuMghFC7BeXHhRC0Wi3u37/P/v4+ALZt8+XGmF/491u0PMF0Pst/9F113rtaVlEA0lOXlbFsKZSCftgeCcOQMAyVnyer8nQ6rcaUd3Z22N/fJ5PJUK1WKRQKTCYTbNvGdV1qtRphGLK/v8/W1hbFYpHp6WnVwz4ej+l0OhQKBbLZ7Jm8tdRoHgdyKG9vbw/DMPjynR0++7W7jI0M07N1ik+OuH797HaHPWqF/tPArwkh/qRhGBng7PXxHAPSZkmShFwud+TUG1CeeKPRYH19HYBCoUC1WuV3X9rnU/+uiReEEMc0dgf8zL/apfueBT74VB3LslTfuOxmkQIqbQ9ZZctOF+nrycpc2irFYlH544PBgHa7zf3798lms5TLZSqVihL2YrHIwsIC/X6fbrfLaDRiZmaGfD5PJpNhMBjQ7XbJ5XJnYleiRnMSBEHAeDxmb2+PL9zZ4ze+uckkTrBSEV3T4e/9m03ml67wo0tLp32pr8pDC7phGEXgg8B/DCCECIDg9f7NeeSwZy47WeRIfxzHDIdDWq0WW1tbTCYT8vk8tVqNTCbDZDLh5z57m8lwSBKHiDjEIEVgWfzLF4f8yAemlU1y2FaRHrvMZzkcryv7ynO5nOqG8X2fTqdDr9dTOw4rlQq1Wo3RaMT+/j7NZpP9/X31Md/3lVi7rsve3h7NZpNCocDU1BSVSoXhcKgsmHK5rLtgNBcaabe0Wi16vR6//fwWYeAjEoO042LlioSk+Huf2+JHv/eZ077cV+VRnqFXgT3g5wzDeDvwJeDHhRCjY7myM4LneUrMZX93HMdqJLjdbtPv90mn08zOzlIulxkOh0wmE3zfp7mzj0hiDNMkZRcwMjZmKsN+fHBi7vs+g8FAHZzKx5PV+quFcR1uVTw8XCQPNre2ttjZ2aFUKlGr1bhy5crBtbwi6p1Oh1KpRLlcZjweUywWmZmZYTgc0u128X2fqakpCoUCmUyGbrerXgzOUhCRRnOcyGKt2WwihKA/8kiSGCvtYNh5rEwBI5WmOYze+IudEo8i6GngO4G/LIT4vGEYPw38NeAnD3+SYRgfBz4OsLy8/AgPd/LIHu1fu7XP3/n1F9na7zGbN/lL31XnHdOmmriUG0tSqRT7+/tqHH80GjHlWrTjPJZTQJgpEpEQT4ZUIp/f//3fJwgCJpOJGjqS7YhSvOWdgOyKORywL3vSZctjJpPBcRwl7Pv7++zt7SmPfGVlhSAI2N7ept1uMxgMKBQK6lS/WCwyNzdHu91md3eXcrlMsVikVqvR7XZpt9sUCoUze8Kv0TwK8rm4t7dHGIa4ZsQkMUmRYObLYGUxrRzzZee0L/U1eRRB3wQ2hRCff+Xvv8iBoB9BCPEp4FMAzz333LfPs59RZDvfrz2/x0/+8vOMxkNiz2Oz2+d/Xr/Hf/b+Vd57raJOu4MgOFKZZzIZ6vU6H//BMj/zm7cZDzvEky7Cn5AyYt7/TJ1Op6PsE1l5S3tFHn4KIfjCWptffX6b/VFEzc3yx942z3tuHFTUhmGow0z5QiDzJvL5PKPRiMFgQK/Xo9FoMD8/z/LyMp7nsbOzQ6fTYTgcKhumUChQqVTU2q0wDNV5QL/fZzAYEIYhlUpFH5ZqLgzSbpHWZRAEfOdCjt8ajTEzDpZbIpWxyLk5PvGRJ0/7cl+ThxZ0IUTTMIwNwzCeFELcBr4PeP74Lu30kP3i6XSav/ubdxgNB4STIcIfkkQ+ESa/8OUG/8Hbl9SruhwcAtRqqiAIWM0M+KGFCf/qa016k4By3uWPP3eVD71tBdu2j/SuywodvtU187u3mvzzF3pMogymCe2Rx//z+bv4scd7rs6qHBlpwaTTaWzbplAoYFkWrusqG6jf7/Piiy9SKBSo1+ssLS0xHA5pNpuqjbFararOmGKxyGg0Io5j1XubTqfp9/vs7+9TrVbP7Gm/RvNWkAVcs9lUz//lqsP3P13nmz2bMFukXsnzk3/6XXzsnQunfbmvyaOecv1l4B+/0uGyBvwnj35Jp4s8BJWWx0Zzh3A0IAlGB4ea6TSQYq/vMRgMGI1GeJ6nvG7btonjmK2tLTqdDqlUivder/P973iCYrGokhTlYahsP5RDQIcPSAH+2r/uYkw/gYsBCKIwQEQ+v7uT8Ce/Z4nRaKSmVqW3P5lM2NnZIZvNKq/ctm3m5ubwPI9Op8NLL71EtVplbm6OlZUVut0u3W6Xvb09BoMBMzMzlMtlZeEMh0PiOMZ1Xaampmi320rUdb+65rwjbdJms6kmRQ3D4MZsgb/0I+/nmWeeYXZ2lqmpqdO+1NflkQRdCPFV4LljupYzgcxiyWQytFotqsaYRjhCiBjDymIYKUQUUE75tFot1SMux+bluqpUKqWmyorF4pGRfNmqeHhBxeEArsPLond6Y17pRgfDIG1nMbI5+ga87W1vYzKZqBbD8XisBo3G4zH9fp/NzU22t7fJ5/PMzMxQrVZZWFhgPB7TarXodDrMz89TqVSwbZt2u81oNGJjY4PBYMDs7Cy2bZMkicp5z+Vy1Go1Op0O+/v7VCoVfViqObfIZgdZ1MhhIsuyKBQKzMzMkMlkztSquddC96Edwvd94jgmnU7TarVoNBr82XfV+Xu/PcDDAsMgDiekJ30++u55NVgkK9ThcEiSJExNTR2pxuXWIjlhZtu2EnHZqigrc1mdS396ab7OVncMCBCCJI5BxMwXs8p/z+VyTE9PE4YhrVaLwWCgxvhlW2W73abVauG6LrOzs9Trdebn5+n1emxubtLr9VhcXGRmZkZ55a1Wi/F4zOzsLMViEUB9zVwuR7lcptfrqa6ZsxgnqtG8EbIibzabRwb4TNNkbm5O3YXK1NKzjBb0V5AtS4A6QAzDkA9enyKdfpqf//11mnstinGfH3rHHE9OHSQoypyVIAiUX53NZtWhZi6XU5664zjKbpH2yuH9oq/GX/3Bp/mJT3+dSXiQ8ZJKg2Ol+Ks//DSFQuFIVIAcMArDUA0WZTIZyuUyy8vL7O7u0mw2eemll9jY2FAHpPl8nna7ze3bt6nX61QqFdWuKMefp6enmZqaUo8zHo/VgJLMiJFeuz4s1ZwnZFdYo9FQtqVsNpifn1eV+nn4vdaCDmqMX/rQjUZDiXsmk+Ejb1vkuQWHtbU14ngaQImxrMDn5ubUIeJwOFSdIdLKeFDIX++X43C64kefnSWOIv72r7/Idm/CfDnHJ37gKf74dy4e+TeHkxvl1qNqtaoq7cFgwPLyMvPz8zQaDTY3N3n55ZfZ2triypUrzM/PMxqN2NraYjAYsLi4yNzcHL/0+Rf5R//mK+wPAmZnpvkrf+Sd/NB3HhzoymlZ2bI5HA4RQpybX36NRtotnU5H2S2+7yOEoFwuMzMzo7KQzgNa0EGtbxNCsL29rdbIwYGgb29vs76+fmQvp/TDZZtfHMd0Oh1s22ZpaYnp6Wl1iyZzzF9tfF5Of8pe8zAMVRa69NE/sOLw/h9725GFF/v7++pFRb7JbBf5dQA1GToYDNjb22M4HLKyssL8/Dybm5usr6/zta99jc3NTa5fv069XlfV+kvDNP/gCy3GiY0wQ3Z2dvkf/+nn8IKAP/5d17BtW7VWyvyZ8XgMoLLXNZqzjDwMbTQaKt5DpisuLCxQKpVURPV54NILutzsLcOq5MIKGY+7tbXF2tqaSjqUgzxyUWypVFJbwefm5lhcXFReshTZB4VNCq4cTJJLLQ5PhB6O1T3M4QiAIAjUQeXhqVHbttUhrfTSZQLjaDRiZ2eH8XjM6uoqc3NzrK2tsbOzwxe+8AUWFha4evUqcRzzs7/4WYZRGrM0jWVYhJM+k0Gf/+NXfp8PXjt4oZBxu/L7dRxHhZjl83nd1qg5s8iGBt/3VSE3GAxUmmq9fpC1dJ5sxEst6EmSMBgM8H1fxdzCwW2YbdtsbGywtrZ2ZAIzjmNKpRLVahXTNNnZ2SGXy3H9+nVmZmaAg7hb27aPiJl85fc8T4m4rKxzudyRdsU3W9keFnd5BiC7XGRPukxjlEusXdfl6tWrqkUxlUrx9NNPMzc3x/3799nc3KTVarG6usp+YGLGPsn+FpSmsQo14lSP/e7BQWoQBNRqNZIkUQNWckG2PDDWoq45q8jlLt1uV50BxXGMaZoqC8m27XN12H+pBb3f76u+7cPpiel0mvX1dTY2NpTYjsdjMpkMtVqNcrmM53kEQcDU1BTXr19XHpusjiVySEG2FB62RqSIPyyH95DKaFwZ7CW9QNlGKavnTCajwraKxaLqfqnVahQKBba2ttje3ubWrVsUgg6DdB4DA9FpIgplUm6F2XKRJEnY3NwkDENqtRrAEVG3bVvZWPl8Xgd7ac4c8nmyubmpzn+kDhy2W85Dd4vk0j7LZP+253lqGhIObsPu37/Pzs4OrusSxzGDwUC1BhYKBUajEUIIlpeXWV1dVXaHbFME1Gm5XIghWxZl/srjQkbtyltF+YIiv1/DMNSLjmmaTE9Pqy6XbrfL6uoqpVKJzc1Nvnt+xK+/tEOULWGaGRj3sJOI//xPfYBr16ZYX1+n0Wjg+z5zc3NKvKVVlMlk1FCSFnXNWULesQZBwM7ODlEUMR6PCYKASqXC8vKyWhZzXuwWuKSCHkWRyjABVNthGIYqkbBYLDIejxkMBlSrVbXZZzAYkM1mWVxcZGFhQbU3yVfxw+vhpHUjRe40fjH+v681+eRnbrPdnTCXT/NffXCJD1+v4HmeOtzNZrPU63Vc16XdbpNKpXBdF8dxsHNr/MbX7jNKuRRyVX7k7dMsmh3iuMyNGzfUi18QBNTrdbUARJ43yCfOaDRSmTUazWkjrclOp0On01HVuWEYagBP3kmfJy6doCdJQrfbpd/vq/RC6W3v7e3R7XYpFotq+4/MB7dtW6UTrq6uKg9dWidw0M8qvXjbtimVSqd6u/ZLX9k60sPeGEb8jd9Yx3EL/OEnD/LO5Xo8eb2O49But+n1ejzxxBO4rst3PTHDYDDAtm3K5Qzj8Zg7d+4wNzfHlStXsG1b5cGvrq5SqVQQQqh2Rnl+AGhR15wJZHbL1tYW/X5fTUJblsXi4iKu6547uwUumaALIRgOh7TbbeBbB5WTyYRms8lwOFRBVq1WS43LZ7NZfN+nXC6zsrJCuVwmlUqpQSG5IUgubZaj8Kd9q/bJz9xWYi6ZhDF/+zde4k+868Oq60WmRE4mE7LZrPqe5cKMfD7P+vo6o9GIdrtNNptVA0y+77OyctCXvra2xksvvcSVK1eo1WpK1C3LUk8g0KKuOV3k8z6KIra3t9XvpjwPWlxcVKmlp/0cfqtcKkH3PI9Wq0UURZimqXwzKebT09MMBgMajQa5XE7lmIRhyNTUFPPz82pzj+M4GIZBGIZ0u13VvneW1rVtdydv+H7XdcnlcoxGI0ajkZqSk36/9Nsty6LRaNDtdtUuVHl3E4YhKysr3Lx5k5deeonbt28TRREzMzNHfHU5PSvzcrSoa04D+TvYarVotVoYhqGKjcO5RudlmOgwl0bQpR0yGo1UL3m322Vra0uNtnuex+bmJpZlUa/XyWazJEmiUglLpdKRqTGZGW4YxpkMqJovO2y9iqg/GNAvBdx1XYbDoRrll0umZReNbds0Gg12dnbUWYFs5QyCgCtXrvDMM89w69YtXnzxRcIwpF6vAwcvHOl0+kilDlrUNSePrMY3NjZUd8t4PCabzbK0tHRwdvRAt9p54VIIuqzEW62WSgvc29tTwwTVahUhBI1GgyRJVGUOUKvVmJ6eVgE9juMcrKd65cXhwFcun8le60985MkjHjoc5MC8VkC/YRgUCgUcx1EWUr/fVxaJnHjNZrNsbm6qASm5EGA8HvPEE0/w7LPP8sILL3D//n2CIGBlZQVAtYDKSl3yeqL+S1/ZUoe682WHT3zkyTNtlPk+AAAgAElEQVSdR6052xy+q9ze3iaOY2Udzs/Ps7i4qOzU82a3wCUQdJnTsr+/r8KjBoMBm5ubjMdjqtUq2WyW9fV1ut0uc3NzyiKQyx5kWJWs2NvtNkEQUCgUznRbkxS+tyqI6XSaSqVCLpej3+8zHA5JpVIUCgU1jZpKpVR2umEYato2CAKeeuopnnzySdbW1tje3iZJElZXVxFC4LoulmWpKv3BxR6HefBQd6s74Sc+/fUj35tG81aQdovsZjNNk/F4TDqdpl6vq0aG8zRMdJgLL+ie5/FP/u0L/O+/+gfsDXxqxSw/uGpyvZBQLpfJ5/O0Wi02NzepVCpUq1Ucx1FCfrh9KY5j9vf3SZKDf3sefugfe+fCQ4ufbdtMTU3hOA6dTgfP89Rhr/TVm82mquLlXdA3vvENbty4wZUrV/h3ax3++s/+Dt0ozcLiKj/+g2/nR/7QE6pSlziO822i/lqHup/8zG0t6Jq3jDwMjeOYjY0NJpMJqVRKrV68du2aauM9r1bghRb0IAj4xc+v8Tf++VeY+D5kHDa2tvn7t3f5C+9/ij88n8P3fe7du6cWUkgxL5fLStwdxyGKIrXQYmpq6lz6aw+DaZoqw73X66nOl1wux/LyMrZts729Tb/fx3VdRqMRnU6Hb3zjG2xT5ue+3GVkFoiDFpsbd/nr/+xgevRH/tATpFIplTs9mUxUm6PkzRzqajRvljAM1WS4tFflGdDs7Cxzc3Nq+cx55Wy0YzwG5Kvx3/1XX2USeJDOEI86iP4esWHxL291yWQy3L9/n+FwyNTUFLlcTgVuFYtFtWxZBnddNjE/jLSgZmZm1HYleZu6vLx8pJUznU7T7Xb5+X/5bxm2G1hukXShRhKFDFtb/PSvfEnd5hqGoaqmyWSiJnbh2w9v3+j9Gs3rIb3zra0tNUA3Go2wLIuVlRVc11XZR+eVCyvok8mE4XBIoz0AUog4JOzuIhBgObQCaDQaNBoNFbYlNw25rqsOB+UWIMMwqNVq5/ZW7DgwTVMtmJZ2k2maVKtVnnjiCbXcI5vN4jgOrVaPsLeD32pipG1S+TJJLGhsbdJoNBgOh+r/52FRl7HBn/jIkzjW0cPm1zvU1WheC3kYGkUR9+/fVzlDURRRqVRYWloCUFbieeVCWi5y8tPzPGZyaXbHMd7+JsKfYObypKwsJTthfX2dVCrF9PQ0MzMzKl1N+uNRFNFutzFNk1qtdmb6y0+bTCbDzMyMWp4h2xqvXr3KvXv36PV6OI5DuVyk2+8SxDGWSEgXKhhOkWr6IK40iiI1lSfzZmT8bi6Xe+hDXY1GIrukNnY7TFsBf/YdVfxmE8MwGA6HRFHEwsKCWmRxHs7FXo8LJ+hxHOP7vhqS+bEPXuNv/YsvwaQPVhYjY5HJWnzffES/36derzM7O6vaEqvVKq7rkiSJqsynpqa0mD9AKpVSffnSjkqShGvXrrG+vk673eaPvmuFf/J7Ef5kRDTaI45D8pUZ/tz7rmKaJru7u8RxzOrqKo7jMJlM8DxP/VmKuhZwzcMgu6TGQYSIQxqjCX/rn/4u7y92uVEvMplMKBQKLC8vK+/8PNstcAEFXaYnynVy779WZu2mwy+PHXqJRa1Q4KPXXZzRNtmsq1atyf2YhUIBgFarBcDU1NSZ7DE/CxiGoQaG2u02hmHgeZ4anX6vaZIyTH75y3fptHuUrDH/4dMF3nftYFhJrv6K45irV68qIQdU18x5HL/WnA1Ul1QckUQBSRIy6ezwu80uK1VbnYktLy+rHQbnvXC7UILu+77a0CMjbhuNBm+vu3zoL3w3uVwOIQRf+MIX6L9yq7W0tKSyxEulEul0WrUm1mo1Hfn6JrBtW92y9no9hsMh8/PzALw3leI9T0wzHA7pdDrk8wdnEsViUbWIDQYDXnrpJW7cuIFt22pBx+FKXaN5q8huqCTyQST4vX3i8ZBxcjCbYhgGKysrFAoFMpnMhfg9O98vR4eQXS2j0UgdsEmPV7bDZTIZ5fGWy2VWV1cpFAqk02lKpRLZbJZOp0MURcqC0bw5ZNtnrVZTVky9XqdWq6lKXmbl7O/vq4GlIAjI5/N4nscLL7xAFEVks1llv8iDUo3mrTJfdhBJjIhD4igm6eyQxD6OZTCZTKhWq6yurqrq/CJ0rz2yoBuGkTIM4yuGYfzKcVzQwyIP1eAgY8X3fbViLZPJkM/n6XQ63Lt3D8uyuHnzJtPT0wBUKhUKhQL9fh/P8ygWi+cuNvMsYBgGpVKJmZkZlW0ju4fkqr3p6YNqXe5vHQwGDAYDisUiYRjy/PPPE0URlmUxGo3UHZdcaafRvFk+8ZEnsQkRcUTsD4jGHdKG4D1LJUzTZHZ2ltnZ2SNhe+ed4/ATfhy4BRSP4Ws9FDIbJAgCJpMJYRiqyUY5xmsYBnfu3GE8HnPjxg1WVlZUKFWxeHBAMhqNVEiV5uGRd0RyRyocRDB0u13y+TymabK3t4cQgpmZGXq9nsrQabfb3Lp1i5s3b2JZlhJ7QOW2azRvhh9+e51u6xo/85svcG9zh5zp854bsywWD9Y2rqyskMlk1M7gi8AjCbphGIvADwH/E/DfHssVvUXkppHRaAQcVOfD4VDtAJU+7ebmJs1mk2q1ytNPP61uscrlMnCwXzSbzapDUc2jYds209PTStjhQJA7nQ6u65JKpVSeRq1Wo9vtEscxy8vL7O3tKVE3TZN+v69E3TRNbYVp3hSe5/H9T9X47uXv5Fd/dYe7dyeqC+7wmrnDS2rOO49qufxvwF8Fktf6BMMwPm4YxhcNw/ji3t7eIz7ctxMEAb7vkyQJ/X4f3/cZDAZqIbMcR3/55ZcJw5CrV68qq0Vu6Ol0OqRSKcrl8oW47TorpNNppqamVMhZtVqlWCwSxzGO4zAzM0O/32d/fx/bttnf3+fevXvMzMwQRRG3bt0CUOmWcj9qFEWn/J1pzjoyEtcwDO7fv0+73SZJEuI4JpVKqd25MsbiovDQgm4YxkeBXSHEl17v84QQnxJCPCeEeE4K6XFxeOOQTFGTUa6maapFxffu3aPT6bCwsMCTTz6JEIJCoUCpVFK3+5VK5dy3LJ1F5GFptVpVnmWpVFKHofV6neFwyO7uLq7r0mq1uHv3LjMzM4RhyIsvvohlWQgh1JPywYgAjeZBZMeb7/usr6+rM5goiigWi6ysrKjq/CLZeI+iYO8DftgwjHvA/wt82DCMf3QsV/UmkROhgLJZJpMJQggymQylUol2u839+/dxHIenn34a13XV/kw5TVosFvVt/GNELgCpVCoqBKlarTKZTHBdl7m5OSaTCd1uF9d12dvb486dO6pSv337trpz6nQ6StRlRIBG8yBycUWz2aTRaCiBT6fTzM3NHTkMvUiF3EN/J0KInxBCLAohVoE/A/xrIcSPHtuVvQFRFKlXYNniNhgMlJjLjfNra2uMRiNWV1dZWVlBCEGxWCSTyTAcDsnlcvoQ9ISQmTn1ep35+Xmq1SrD4RDHcZidnVW96nIBydraGlNTU8RxzN27d1X2hvTb5Yu3RnOYIAjUmsT79++rrH7DMMhms1y7dg3btpWgXyTO7UmAFHJAibn0VoUQVKtVNjY22NraolKpcPPmTQC191PGwMrDNs3JkM/n1QpAWXXv7e3hui7z8/M0X8nZcF1XLSBYXFxUS0nq9TpBEKhZApn/os8+NJLhcAgc/F5tbW0xHo+Vdz49PU29Xlfe+UU5DJUcy3cjhPgd4HeO42u9GWSbYhRFamO9/KGl02kymQyj0Yh79+4BcPXqVSqVihogkoNH1Wr1Qt1unRccx1EDR7IDZnd3l3w+z+LiIltbW8DBi2+z2QRQtkyz2WR6ehohBL1ej2KxqKZKNRo5syCEYH19nX6/TxiGatPW8vKyKiouYqzEuXt5km2K8iBULn+Ooki92ubzeW7fvs3Ozg6zs7Ncv35d7cuUvc1yaYPmdLBtm1qtBqCeVM1mE8dxWFhYoNlsMhqNyGazyO6o6elpgiBgf3+fqakpFREgRf28BytpHh1ZnbfbbRXRLLuqisUiS0tLymq5CJOhD3LuBF2KuOxqkQehpmmqqdBer8f6+jqZTIZr166Rz+fJZrPk83nVn65989NHxvDKjiQ4yKjPZrPKfhmNRuRyOTqdDgDVahXf92m322qGQAhBqVRSXU2ay0kcx4zHYwA2Njbo9/uMx2Nl8S0vL1OpVFR1fhFD986VoAshCIIAz/OOiHoURdi2TZIk5HI5Xn75ZQaDAYuLiywsLGBZFoVCgTiOdYviGSOdTn9b1ry0XJaWltje3mY8HuM4Dr1eD4ByuUwYhnS7XTUIJj15wzB0x9IlZTAYAAedUM1mk3a7TRzHKvpjeXlZTYZe1Lu5cyXoUszDMFS2y+GtN7Zts7u7y+bmJq7rqkjWXC5HJpNhMpmQz+f1E/6MIQeQpPUihGB7e5skSZSnPhqN1PILuTlJCMFwOCRJkiNeqLxb01weZHVumibb29t0u116vZ66a5ubm2NqagoA13UvrAacG0GX3rl8C8NQ+WPZbJYoijBNk/X1dXzf58qVK0xPT6t1aGEYYlmWtlrOKHIA6bAwr6+vI4RgeXmZzc1NRqOREus4jikUCpimqXqO5b81DINcLqfvwi4Rw+EQwzDo9Xo0Gg329vaUNti2zbVr19SfL9Jk6IOcG0GXvaXyzff9I/6Y3JzTbDbV7ZWsztPpNGEYaqvljCP3k8qfUZIkbG5uArC4uMj6+jrdbhc46GmXcwe2bR8Rdfk5ruteuC4Gzbcjq3PDMNje3qbf76sl0DLGWW4du+h36OdC0A9X51LQ5T5AubQCvuW9LiwsUCwWyefzqjrP5XK6q+UcYJomlUpF/V2KumEYrK6ucvfu3SMHpDKULZvNqj+Xy2V6vZ6q1LWoX2xkZ0u/32dnZ4dms0kYhhSLRXK5HCsrKyqA66LfuZ0LQZf9yocrdM/zVFdLHMe0Wi329vYoFovMzMzguq7qTb5oATwXHRkVcNhT39zcxDRNnnjiCe7cuUOn08G2bfL5PP1+HyEEruuqW28ZC2Capu5Rv8BEUaQ6W7a3t+n1eirsLZ1OU61WmZ6eVm3LF72oOxeCLpPT5EDReDwmCAJyuRxJkhAEAdvb28RxzOzsrAreymQyJEmi1p1pzg+GYaj0SyEEQgg2NjYAeOKJJ3jxxRfZ2dkBDuYORqMRSZLgOI7qdpBfR/eoX1zkz7rf79NqtdjZ2VF35LlcjitXruA4DqlUSsU2X2TOhcpJy0VOh8pXZMdx8H2ffr/P3t6e2o6Tz+dVP/JFblG66MgNSFLU4zhmY2MDwzB48sknuXXrFjs7O+oObDQaqalR6aN/5vkmP//5z7EfZliYKvKJjzzJx965cMrfmeY4kHfqSZKws7NDq9Wi2WySSqXIZrPUajVmZ2fVIpvLoAPnQtDhQJgPty1K73wymaiDs3q9rpIU5cGH67oX2jO76MhJ0JWVFZVn3Wg0SKVSXL9+nRdeeIHt7W3VY3w4bfNXvvAS/+Dze4RmGjPrsrEPP/HprwNoUT/nCCEYDAYkSaJC3ba3t/F9n1wuh+M4LC8vY9s2qVRKTYlfdM6F0pmmiWEYqvdc5pnLsX85NVgsFtWyZ1mp6cnB84/0P69evcp3fMd3MDMzw/b2NoPBgKeeeopUKsXm5iZhGCp7LggCfuH31piM+8RhQDIZEHsDRp7PJz9z+7S/Jc0jIhskkiRR3W1yh7DjOCqEyzRNstnspTlDOxeCHscxu7u7JEnCeDxWU6Hj8ZiNjQ1SqRT1eh3LsiiXy6RSqQu1+FWDum2+du0azz77LLVajY2NDXzf58aNGxiGQaPRUIehk8mE3cEEEXrgH3RExeMh8aTPZnt4yt+N5lGQVblsV+x0OjQaDYIgUHMnckLcNM0jd+wXnXMh6LJTQaaoVSoVfN+n2+3S7/eZmprCdV0KhQLZbFa9Sl+WH+JlQYr69evXefbZZykWiypR89q1a4RhyM7ODqZpHvyeZASxgCScYPgDEhERT/pMW6FejnGOkbZrFEW02222trbY3d1Vh98zMzMsLBxYajLD6bIUdudC0IUQpNNpxuMxlmWRTqcZjUbs7u5imiYzMzOqOrcsSx+EXmBkVvpTTz3FO9/5TjKZDGtra1iWxcrKCuPxmEajgWVZfOxtc9gIYkxi34dxn4xp8uffWVHr7DTnC9kUIc/TWq0WGxsbxHGMbdsUCgWWlpZUMVepVC6V7XouBN00TdW2WKlU8DyPXq9Hp9OhXC6Ty+XI5/NqaMBxHN2meIGRon7z5k2ee+45LMvizp07lEolFhcXGQ6HNJtNPvD0In/uuTlqNghTUM4k/NhzVb73+hSNRoN2u603Hp0j5PmI7GyR6yXl4aicCq3X6yRJogT+MjVFnAvVkz88Kdbb29vqAGRubk6N9Nq2feGWvmpeHWm/PPvss8RxzJe//GXW1ta4du0aQRCwu7tLOp3me59Z4t1XquqAPAxDWq0WlUqFRqMB8G0ZMpqzyeGk1SAI2NnZodFoEIYhmUyGSqXC0tKSWj9ZqVQunRacG0GXG4Z836fT6dDpdCgUCqpFqVgsKu/8og8PaL5FLpfjHe94B2EY8gd/8Aesra1x9epVTNOk2WxiWRaVSkXFAliWpaZJi8UijUZDxQ1oUT+7yMXgcqlJp9Phzp07avYgk8mwuLioFotnMhnK5fKlqs7hnAh6KpViaWmJXq/H9vY2rVaLJEnU0EAul1OtSZftFVlzIOrvfve7iaKIr3/962xubrK4uEgURWrwqFgsqh20pmnS6/VIkoRyucz29jbApRSA84JcBC9bUzc3N9nb2yNJEnWOtri4CHxrp/Bl1IJzIehCCJIkIYoiWq0W3W6XYrFIsVgkk8lQLBaxLItsNqufkJeUXC7He9/7XoIg4IUXXqDZbKpOh8OCLTdWASrAy3VdLepnGGmxjMdj0uk0W1tbrK2tEQSBSttcWVlR1XkqlaJarV7KO/VzIegyVnV9fV3lHMslw7lcjkKhgOu6l/IVWfMtcrkcH/zgB4njmJdffhkhhDogazQapNNptYYwm82SJAmDwUCFeUlPXYv62SFJEjzPUwefURTx0ksv0el01AayhYUF5ufn1efXarVL2+V2LgRdCMFoNKLZbDIcDsnlctRqNZIkoVAokMlksG1bPwk1uK7L93zP9xBFERsbG1iWpSr1zc1NlpaWyOfzKq0zjmO12u5wSmOpVLqUFd5Zw/M8PM/D932y2Sy3bt1ic3NTibnruqysrFCtVgnDENM0mZqaurQ/u3OjgJ1Oh3a7TRAEqitBtiXp6lxzmEKhwIc//GHm5+fpdrvs7e2xsLBAqVRia2tLFQUSKerD4ZDBYKBWmMVxfIrfhUaO9w8GAzV7cuvWLZV/DrC8vMzi4iJJkiCEoFwuX+qtZOeiQk+SRD3JbNumUqkQxzGVSkWJua7ONYcpFot86EMf4rd/+7dpNpsAzM7Oqp2TcGCtyJ7mOI7V1iP5JheK65mGkyeOYyXmct3gV77yFdrtNnBwSFosFrlx4wau66qMp1qtdmmrc3gEQTcMYwn4v4A5IAE+JYT46eO6sMMcblmamprCcRxM06RYLOrqXPOaVKtVPvCBD/DZz36W/f19AKanp0mlUjSbTdUN4XkecRwThiGTyYQkSY5MkcoJZM3JIIRQVotc7P7yyy9z584dJpOJ2h987do1pqeniaJI5efn8/nTvvxT5VFKjwj474QQXzYMowB8yTCM3xBCPH9M16aQi4Ft26ZWqxFFEfl8nnK5rDtbNK/LzMwM73vf+/jc5z5Hq9UCDoTeMAwV+DY3N6fa4oIgAFD7SaVXWy6XdeFwQsjdB/1+H8uyCMOQ559/Xs0ShGHI9PQ0Tz31FJZlqdTFy16dwyN46EKIhhDiy6/8eQDcAh5LyLRhGERRhOu6FItFhBCUSqUjiyw0mtdibm6O97znPdRqNTzPU3HLpVKJdrtNs9nEcRxV3UlBGQwGDIdDFRMg+9g1j48wDFUstnwh/drXvsbu7i5RFBFFEel0mps3b1IulwnDUH3eZa/O4ZgORQ3DWAXeCXz+VT72ccMwvmgYxhf39vYe6utLP61UKql842q1qqdCNW8KwzBYWFjgXe96F9PT04RhqCaNy+Uy7XabjY0NstkspVIJQNkwvV6PwWDAzs4O+/v7qkrUHD9xHCubxfM8KpWKslp831dnHXNzc9y4cUMdhALUajV91sExHIoahpEH/hnwXwsh+g9+XAjxKeBTAM8999xDJSH98le3+alf36DRajHt3OcvfuhZ3v3ukq7ONW8aOW0M8NWvfpW9vT36/b7aW9rv97l7965qgWu326pfvdvtKk9dWjCXKZL1JJDbx+RBqOM4jEYjvvnNb6p9sePxGMdxeNe73kUmk1EvuuVymUKhcNrfwpngkSp0wzAsDsT8HwshPn08l3SUX/rKFj/5L27RMYoYhkV7FPL3P7/Db77Y0a/ImrdEOp2mXq/ztre9jfn5eVWpu65LpVJhMplw//59DMOgVqthmia+7yOEYDgc0uv16Ha77O7u0ul0dFLjMSLFvN/vq0ymb3zjG+zu7qqR/3Q6zcrKCsvLy6rnPJVK6er8EA8t6MZBefKzwC0hxN85vks6yic/c5uRHyIiH0NEkLYJMy4/8zv3HtdDai4w2WyWubk5bt68ycLCAlEU0e12yWQyyo65e/cucRwzMzNDKpViMpkQxzH9fp9ut8twOKTVarG7u6t71Y8BKdiyEi+Xy9y5c4c7d+6oA09Zib/jHe8giiL1vmKxqKvzQzxKhf4+4M8DHzYM46uvvP3gMV2XYrt7kK5mptIgBIbtkLJzNAfhcT+U5pKQy+VUl8T8/Dy+79Pv9zEMg+npaUzT5P79+3iep1YbSlHvdru02221XKHZbKrOGM1bR+a0DIdDdU7WbDa5desWnucRRZGKx7158yazs7MEQYBhGKTTad3Z8gAPfZ8ihPgc8NhNxPmyw3prSBL6GGYaK1fETGWYLzuP+6E1FxSZpZ4kCTdv3iSdTrOxscFgMCCXy1GpVOh0Oqyvr7O0tES9Xmd3dxff90mn08pumZ6eptPpEIYhc3NzlzY/5GGJogjP8xiNRoRhqHJ2XnzxRdrttuo2MgyD+fl5nnnmGRWXC6gpcc23OPMN3J/4yJM4VgoEGGkbM1fCydp84iNPnvalac4xpmmqLpcbN26wvLxMFEWqw6JWq2HbNhsbG3Q6Her1Oo7jkCQJhmHQ7XbVsFK/32dzc1NlwmjeGLngWW4is22bKIq4e/cuGxsbSuTltO53fMd3qOwd6Z1Xq1XtnT/Amf+/8bF3LiCE4H/4+W12+7A8XeG//+izfOydj6XlXXOJSKVSFAoFhBDcvHkTOAjw8n2fKIqoVqv0ej0ajQZJkqiKXI6Zt1otleg4mUxoNptMJhNmZmb0sNvrEMcxk8lEibkU5Wazyd27d9W0bhRF2LbN6uoqKysrKu5Y/tx0df7tnHlBB/jht9e5+V9+ANd1WVhY0E8WzbFhWRb5fB4hBE8//TSGYdBoNJRfXqlUEELQbDaJ41gFw8mWOTltWq/XCYKA/f19giBQ3vsvfWWLT37mNtvdCfNlh0985MlLXYzI9kN5AJpKpUilUuzt7bG2tka73SaKIoIgwDRN6vU6Tz/9NL7vK9/cNE3K5bKuzl+Fc/F/JJVKUa/X1Q9TozlObNtW+S3PPPOMEvUoitQictu22d/fV5V6Op1Wor6zs6NyuQHlq//77Yi//ut3mYQHnTBb3Qk/8emvA1xKUU+ShNFoxHA4JEkSJcidTofNzU22traAg0ldOMjduXHjBvl8XrUzplIpSqUSjqPP0F6Nc6GOQgiVea7RPA7k6L9pmjzzzDPMz8+TyWTUlizXdSmXy3S7Xba3t8lkMjiOQ6lUIpvNsrOzw/r6uop1Hg6H/K1//u8YtPeOBH1NwphPfub2KX6np0Mcx0fE3LZtdRaxu7vL+vo6vu8ru6VYLLK8vMzKyoraAWvbNplMhkKhoKvz1+BcCLphGDiOo3+ImseK4zi4rntE1F3XVX65bduUSiWCIKDRaCgLoFQq4bouu7u73LlzhyRJyOVyNLsjokmfqL9HHPrqcba7lysT5kExdxznSAvo/fv36fV6apArk8kwPz/P9evXj6yZM00T13V1N9HrcC4EXaM5CeRKQ8dxMAyDmzdvUq/Xyefzyn5Jp9MUi0WCIFD+uYxunZqaotVqcfv2bTzPY65cBCGIvRFRb5d43EcIcalabqMoYjgcqrCtfD5PGIYqI2dzc5OdnR0ymYzaGbqwsMDKygrZbBbP8wDIZDJkMhlc19WF3eugBV2jOYRpmuRyOWXvPfHEEywuLlKtVkmShG63ixCCqakpVblHUYQQglwux+zsLP1+n29+85v8mXfVyDk5DMtChD7RsE163OHHP7R6qt/jSSFH+aX/XSwW8TyPXq+H53msr6+zsbGB4zgMBgOiKGJ2dpa5uTnm5+fxPE/dncs4AO2dvz5a0DWaB0ilUkrUU6mUWnM2PT2tFkt7nke5XEYIQa/XU9Oi2WyW+fl5giBgkQ4/9u4K9VIB084xk8/wX3z3NG8vh7RarQsdG+B5nopJsCyLUqnEaDSi2+0SBAGbm5usr6+rylxGLczOznLlyhWVoWPbNtlsFsuycF1XT4W+AfreRaN5FSzLOlINLiwsYFkW6XSa3d1dhsMhYRhSqVRUFSrtGsuymJ+fZ29vj2VrxN/+wXlWVlaURzyZTNRu06mpKXK53IXp3jrcyRKGIblcjkKhwP7+PsPhkDiO2djY4M6dOyoxcTgcUqlUmJ+fZ2lpiTiO1RpAKeaWZWnv/E2gBV2jeQ0ymYzK3DZNk5mZGSXqjUZDRb2Wy2X6/b4aSMrn81iWxezsLK1Wi83NTUajEXtsA84AABh5SURBVDdu3CCbzZJOpxmPx7TbbUajEdVqVbVGnudIXrkURI7nF4tF1QEkvfD19XXW1tawLIsoiuj1epTLZer1OrOzsxQKBRXIlcvlcF1X/T+9KC96jxMt6BrN65DNZhFCqByXSqWCZVmkUil2dnbUgIw8KJVBU3LX7czMDJ1Oh263yze/+U2uXr2qRtblyrvt7W0GgwFTU1O4rksmkzlXwi6EYDweKytKtoDKgawwPAjSW1tb4969e2qeZH9/n1wux/z8PPV6Xe13lR+XU7y2beuW5TeJFnSN5g14UNTz+TzXrl1Tlfp4PKbf76vKXB78lUoHS1impqZIpVJ0u11u377N/9/eucdIdl91/nPurbr1ru56dHveyThOZnAQEWYShyTEyBAwZpfA8hAggVGCrAixj3+AoEhI+9CyLPvQLgRQSCJAywaySQhelJDHxjGbRTZ2gpNJ4ok9djKe6Znp7umpd9163Krf/lF1fr7d7p7pyfRren4fqVTVXXeqTv2q53vPPb/zOHbsGHfccYcVby2D73a7FItFSqWSDd3sda+03+/TarXodrsAdnB7t9ulVqshIozHY86cOcPi4iLJZBIR4fLlyySTSY4dO8b8/DyVSoV+v082m2U0Gtnxkh9/6lu874klLrcjV2m7CZygOxzXQTMtjDHW24RJBkwqlWJhYYFGo0EURTbMoOGHbDZLOp2mVCqRSCRoNBp885vfpNVqcfToUfL5PIlEguFwSLfbpdFo0G637czcTCZDEAR7bjNwOBzSbrdtP5ZsNmv3Aur1Op1OB9/3abVanDlzhnq9bk9QFy9eJJlMcuTIEebn56lWq6sadKVSKWZnZ/noE8/z259+gYE/iZ3f7pW2m8EJusNxDeK9WA4WU/zz+47xI6+dxxjDaDTi+PHjBEHA+fPnaTabNBoNW0E6GAxsBkculyOfz1tPfWlpiTAMOXz4sJ2Pm0gkbNimVqvRbDZtEyr12NXD3S2iKLIeufYpL5fLNltFKz2TySQLCwucO3eObrdr9wcuXryIiNisoVKpxGg0wvd9ezUyPz/PcDjkDx97gb4kV/Xo1kpbJ+jr4wTd4diAj//jAr/5sdO2F8vFZp9//cmzADz4nXcAE4F7xSteQRAEvPjiizQaDfr9PisrK5RKJQqFAp1Oh1qtZjNakskktVqNTqfDiy++SLvdtgVMQRAwGAwYDod2RF6n0yGVSpHL5chmsyQSCbs5u1PiHoYhnU7HThcKgoBSqUQQBPT7fdsj3vd9O/VJe7PoFYsOiY975prNksvlCMOQubk5fN+n3W5zOQQv+fIrk9ut0vZGcILucGzA737qG1bMlV5k+L3HzvFPX3cIeGm48dGjR23/9KWlJaIoYnFxkVKpxMzMDGEYcuXKFRtjLpVKNr1veXmZMAwpl8vMz8/bMMtgMLD9ZMIwpN/v2ysAPSaRSJBIJGzjqq0SeGMMvV6PXq9Hv99nMBgwGo1Ip9MUCgWSySTD4dAWWqXTaaIoYmFhgUuXLrGysmKvJkajEVevXmU8HnPs2DHm5uaoVCo2g6hardJut8nn85TLZdtm4Uh1louN3stsu50qbW8UJ+gOxwZs5Aleag7IZDI2lNDr9Wi1Whw8eNDmTS8tLeH7PlevXmU0GlEoFMhkMtTrdfr9PqVSaTJa0fPodDq2UrLValGtVimXy2SzWdtKNplMYoyxWTRhGOJ5ns2XV1GPi7vneYjIqlsc9Y41fDQajWzIZzgcEkURn/zKAn/8hRdZ7I44XJ3hX9x/Fw/cnaLT6QDY915eXub8+fMsLy8zHA5tJtAXnr3MR77wdVrdIfMHDvKL+aOcOFHGGGPnh3Y6HTKZDAcPHqTdbmOMoVgs8usPnFx1hQSQSfpuuM01cILucGzAodkMC+uI+qHZl0rRwzC099o/XfuOLC8v4/s+nU6HMAypVqvMzc3ZaUezs7MUi0USiQRhGNq5per5lkol5ubmKBQKDIdD+v2+3RyNosjea2Mrz/NsOAawIRkVchHBGAO8XMx18LIe5/s+jz5X5z//3SXCoYAIF5Zr/NZHnqL/o3fzk/e+ilQqxcrKCufPn+fKlSv0ej3bSdH3fT7zpef5n3//LNF4hBTK1L0sH/iHJQqlMt/3qsmVy2g0IggCDhw4YAdf6BWAxsldP/nN4wTd4diAX/vhE9f0EBOJxCpR7/V61Go1yuUyJ0+eJAgCrly5YvOzL126RKVSoVQq2YwWrS5Vb3o0Gtmy9zAMaTQaFItFqtUq+XweeKlfuAp1XKTH47HtgRKG4Srv3Pd9K+L67/UqIZ1Or9qYNMbw+48+Qbfbgenx+EkGfoo/enyRB7/rEGfPnmV5edkOcwZsV8R6vc5fP/k8w5GQKMyTyBXxkil6Bj742HP8yOvux/M8xuMx8/Pztt+8lvgrP/7dh52A3wBO0B2ODdiMh5hIJGxXwGw2aytAS6USJ06cIJPJsLS0RKvVsvHyTqfD/Pw86XSaTqezqrpUxR8m7Qc+e/oCH37yHMvhiMNzVX7lh7+Ln37TCQqFAv1+f5WYxrsQxr1ymAh03FuPp0HGQy6A9dQvN7qIH4Dvg+cj4xFRr8m5s+c5fXpS7Toej+0mqc78XFlZmWT8DCAozeNl8vh+EuMnIBqwHHp2b6BarVIoFKjX63Z+6K1UVLXXcILucFyDzXiIGuJQUQ/D0Ir68ePHSaVSLC4ukk6nSaVStm3s3Nwc+Xyefr9v48gzMzMkk0m63S6PffUC/+PJi/RJgAiXLi3ybz/0GI2rV/jR75lUnOZyObtB2e/3GY1GdrMxLuB6BRAnLvrqwesAZt/3OVQpcbHWZjzoYYYDorDBeNChnJxsmKqnXygUbM65jpCrVCocOHSYOgEeCYznQzRE/ICD5SKj0YhSqcTs7Cy9Xs+GWoIg2Oqv8LbCCbrDsQWsFXURsePrjhw5QjabtQU16XSaWq3G5cuXyefzVKtVYBIP932ffD5PKpXiY6dP0+t1GXsJ/ESScZBmgPAnjz3DPYcmnn8ul6NQKNixbNlsFsB63NraV1lvk1TvdTCz3r/j9RX+/SMLDDptRsMOZhgR+MJPfM9hoiiyfWl0tqqGgnRT9+G3lfm9z5+lZzxkFCGJgEwmzcP33cns7CwzMzOMx2MajQapVMqGlBzfPk7QHY4tIi7q2hmwXq/b4dKpVIrLly9Tq9VIpVI0Gg3q9Trnzp1jfn6eXC7HcDjEGEMmk6Fp0nhZD4Yhpt/BDAaMg4ArozT5fN628m02mywsLNhc9Xw+b68GNBSiG6LqvavI6896Auj3+/b22pmIX3htlg8/uchSP6SaS/Ezr38l93/npPOkhle63a6Nw8/MzNixfD9YTpDMBLz/sedZaRsOlIq86/7X8M/ufZUNMa2srAAwOzvrQi1bwE0Juog8APw3wAfeb4z5D1tilcNxixIXdZ181Gq1Vs3JzGQyXL161XqlS0tLrKys0Gg0KJfL+L5Pv99nfjbHUiuJl0gyCvKMe228XodiMMlxLxaLtlBJY9n1ep2VlRWMMasKkDTzJL45GhfzKIqskMfj8vcen+Etr3n9qtfSQqIwnGQAJRIJW/qvWTvpdBoR4f4T89x3V8U+rzn06XSaZrPJYDCgVCrtudYGtyrftqCLiA+8F3gbcAF4UkQeMcZ8fauMczhuRbSgJgxDK6SdTofRaGRbxWrYJZlMUigUWF5eti0BtIf4z99ziD/8v99kkEjhJcaIJ6STs/zc6+9gOByysLBg+8xonnuxWLSx8CiK7K3RaFgRVzTermP0fN9nPB7bUIqeFLT6U4uber0eQRDYk5e+rxY5FYtF28wMsJ8nnU7b4dq9Xo9Op2PbGji2hpvx0N8AnDXGvAAgIn8BvB1wgu647YmnNCYSCQqFAu12226WVioV0um0LTY6cOAAlUqFxcVFO6Lt7kqad7zxIB99+gpLnTEHKiUeuvcIb75rMg7P8zza7TaNRoPFxcVVAqxesPZQ0asFDbeodw6sym6JZ8pEUUSv12M0GtlUSA3hRFFEJpOxrW0TiQS5XM6GjTRlUqtiVewzmQxRFNlmXcVicVe+n/3KzQj6YeB87OcLwL1rDxKRh4GHAY4dO3YTb+dw3FpoKELFbWZmhmazaYuKtH2uhh86nQ5Hjx6lXC5z9epV2u023zEr/LsHjpLJZKhUKja8ohWduVyOmZkZ65G3Wi2Gw6Ft4RtFEaPRCBGxYh1vE6AFSRry0FF68Rh7EAQUCgUb1hER6217nmdtA+h0OgyHQ9vPvFKp2Awb3bCt1WqAi5tvBzcj6Ot9E+ZlvzDmfcD7AE6dOvWy5x2O/YzOJ9Wc7ZmZGVqtFrVajUKhQDabtdWlmUyGTqdDEARks1k7g7PdbtuipVKpRLFYtGEN7beiYZNCoWDDJXrTjJb4hujazVFgVS66DtnwPI/BYGCrWPXEEASB9b7T6bS10RhDEAQ2i0UFW68QtBVCpVJZdTXg2BpuZkUvAEdjPx8BLt6cOQ7H/sPzPDv0QXum61AMba1bKBRsVkoYhjZTJpvNMhgMqNfr9kRQr9dJp9MUi0UKhQLlchmYdEQcDAY2H93zPCvea1MUlXi1qFZ5ak8XjZlrrF3TIrPZrG35q7nnOqYvlUpRLpcpFot2YzWTmbRKqNVqto+NyzffHm5G0J8EXi0ix4EF4GeBn98SqxyOfYaGHMIwZDgc2qZWKvLaOldni6qoa9fFdDpNtVql2+3SbrdptVosLi6ysrJCNpu1sez4wGkR4dEzy3zgCy+w1OgyX0jx0Pce462vmVvV12U8HtuURb2S0FCNdlfU+Z75fB4RYTAYcOXKFbvJqk3CqtUq2WzWVruqPc1m0xYPuU3Q7ePbFnRjTCQivwp8ikna4geNMV/bMsscjn2Girr2FNdxdmEY0mw27XQjja1rJ0cV9H6/TxAE5PN55ubm7Ni6drtNp9OxG6IaF//7s8t84P99i340RsTjYs/jv/5tnV7/NXzfXVWbrqjirf3MU6mUHZ+nsXJt9KUVqfBSLN73fdtvRrN7NLbueZ5te5DP5ykUCrv8LexvbiqIZYz5BPCJLbLF4bgtSKfTeJ5Hv9/H8zzy+Ty9Xo92u20rMFVEdYDEYDAglUoxHA4ZDAY2JFMsFq13rW1vYRIP//AXF+gPxxhGEA2BMWE45k8++zQnZ15nNyu1elVb/+otfnJQwfe8SR8WDekEQWDno6oNmmWjOfg6is9ltGw/blfC4dgFNP1PNxs191tj1jpjM5VKkUwmV/Up19i6Cuh4PLbNuvTfjsdj6qMUiXwWg+AlPGAi4G3g5MmTwEtNu+K9XjSPPpFI2JDKcDi0r63ph5VKxW7OhmFIFEX25ADQbDZpt9tkMhlmZ2d3Y5lvO5ygOxy7hO/7dvSa9nHJZDKrNjaTySSpVMp67SrqehKIosimKapnrWGNg3NlFhtdzNggIwMyAs/jQCGzahiGFhVpXF03QfVko2i4RzdjNVVSs1vUwwdoNBp0Oh2y2awT8x3ECbrDsYtoXF29be2JorFtY4z1fFXYU6mUFfZkMmlTCHW60WAwIIoiHn7rnfynzzxHbxAxKeweE/ge77zvLkTExs5hdU8XjaVrCb9eKeRyORsuinvlvu+v+n29XicMQyv+jp3DCbrDsQfQFgHaR0U3IbU0XwVcS+5VZNVDV2FNpVLAJJTyc/eVyeaLvPdzz3Kp1uGOQoZffssredvdd6zKP9c4ucbEE4mELTbSOH4ymbTZM3oFoXZrCqLODh0Oh8zMzKwaVOHYGST+xW43p06dMk899dSOvZ/DcauhxUJRFNnYdnxDUv+/qtBqhWe8yZaGTJT1hl3ovZ404s/rKDsVdD1WN2S1xYDmrQO2YZcxhnK5bE8sjq1BRL5ojDl1veOch+5w7CG070q8j4o2zVJR11F12rZWuyDqDVgVC1/bu0XR11KPXL3yePGRniS0re/aE8l4PKbZbNLtdkkkEpTLZVcBuou4lXc49iDa7Krf79v+KfHKT83/jvcwj3vWce/6RtCY/drB0WuFHCZeebPZtIVRhULB9WbZZZygOxx7lHgTLPXW431ZRqOR3bgE7KaoEo+Nx9MS41Wi8d4u+pqKbnZqPF8ZDoe0Wi36/T6+77sQyx7CCbrDscfxPI9sNmuHUGhvFo2vR1FkQy+pVGrVBCIt49/Me+imaLwbo6Ix9DAMbfdIbQmgcXTH7uME3eG4RdAYeVzY4SWPW0MvGnaJb1rGuy8qa4dEr0XDL1EU2cZfGuPP5XIuVr4Hcd+Iw3GLocKuXRHXeuAqwv1+3xYNaRm/evawOuslvnmqnn28iEmHV2sIxrE3cd+Mw3GLopWlGg6JFwrFvXH93bVCIxpD11u8b4u2H3BzP/c+TtAdjlscHQ0XBIH1sDW2Ht/8jKcuxgU/3s9FvX/NknHe+K2F+7Ycjn2Ebm5qT5V4Bst6cXRg1Si6eDaM49bDCbrDsY+Jx9Ad+x+Xb+RwOBz7BCfoDofDsU9wgu5wOBz7BCfoDofDsU9wgu5wOBz7BCfoDofDsU9wgu5wOBz7BCfoDofDsU/Y0RF0IrIMnLuJl6gCV7bInK1kL9q1F20CZ9eN4uy6MfarXa8wxsxd76AdFfSbRUSe2sxcvZ1mL9q1F20CZ9eN4uy6MW53u1zIxeFwOPYJTtAdDodjn3CrCfr7dtuADdiLdu1Fm8DZdaM4u26M29quWyqG7nA4HI6NudU8dIfD4XBswJ4WdBH5XRE5IyJfEZG/EpHZDY57QES+ISJnReTdO2DXT4vI10RkLCIb7lyLyLdE5LSIPC0iT+0Rm3Z6rcoi8hkReW56X9rguNF0nZ4WkUe20Z5rfn4RSYnIX06ff0JEXrldttygXb8kIsuxNfrlHbDpgyKyJCJf3eB5EZH/PrX5KyJyz3bbtEm7vl9EGrG1+q0dsuuoiDwqIs9M/y/+y3WO2d41i08x2Ws34IeAxPTx7wC/s84xPvA8cCcQAF8G7t5mu74DOAF8Hjh1jeO+BVR3aK2ua9MurdV/BN49ffzu9b7D6XPtHVij635+4FeAP5o+/lngL/eIXb8E/P5O/C3F3vOtwD3AVzd4/kHgk4AAbwSe2CN2fT/wNzu5VtP3PQjcM31cAJ5d53vc1jXb0x66MebTxhgdaf44cGSdw94AnDXGvGCMGQB/Abx9m+16xhjzje18jxtlkzbt+FpNX/9Pp4//FPjxbX6/a7GZzx+39yPAD8j2z2Tbje/luhhj/g64eo1D3g78mZnwODArIgf3gF27gjHmkjHmS9PHLeAZ4PCaw7Z1zfa0oK/hHUzObGs5DJyP/XyBly/ibmGAT4vIF0Xk4d02ht1ZqzuMMZdg8gcPzG9wXFpEnhKRx0Vku0R/M5/fHjN1JhpAZZvsuRG7AH5yepn+ERE5us02bYa9/H/ve0XkyyLySRF57U6/+TRU993AE2ue2tY12/WZoiLyWeDAOk+9xxjz19Nj3gNEwJ+v9xLr/O6mU3c2Y9cmeLMx5qKIzAOfEZEzU+9it2za8bW6gZc5Nl2rO4HPichpY8zzN2vbGjbz+bdlja7DZt7zfwMfMsb0ReRdTK4i7t9mu67HbqzVZvgSk1L5tog8CHwcePVOvbmI5IGPAv/KGNNc+/Q6/2TL1mzXBd0Y84PXel5EHgL+CfADZhqEWsMFIO6tHAEubrddm3yNi9P7JRH5KyaX1t+2oG+BTTu+ViKyKCIHjTGXppeWSxu8hq7VCyLyeSbezVYL+mY+vx5zQUQSwAzbf3l/XbuMMSuxH/+YyZ7SbrMtf083S1xEjTGfEJE/EJGqMWbbe7yISJKJmP+5MeZj6xyyrWu2p0MuIvIA8BvAjxljuhsc9iTwahE5LiIBk42sbcuS2CwikhORgj5mssG77q78DrIba/UI8ND08UPAy64kRKQkIqnp4yrwZuDr22DLZj5/3N6fAj63gSOxo3atibP+GJP47G7zCPCL08yNNwINDa/tJiJyQPc9ROQNTHRu5dr/akveV4APAM8YY/7LBodt75rt9E7wDe4an2USb3p6etPsg0PAJ9bsHD/LxKN7zw7Y9RNMzrR9YBH41Fq7mGQsfHl6+9p227UZm3ZprSrA/wGem96Xp78/Bbx/+vhNwOnpWp0G3rmN9rzs8wP/honTAJAG/tf0b+8fgDu3e402addvT/+Ovgw8CpzcAZs+BFwChtO/rXcC7wLeNX1egPdObT7NNTK+dtiuX42t1ePAm3bIrrcwCZ98JaZZD+7kmrlKUYfD4dgn7OmQi8PhcDg2jxN0h8Ph2Cc4QXc4HI59ghN0h8Ph2Cc4QXc4HI59ghN0h8Ph2Cc4QXc4HI59ghN0h8Ph2Cf8f/c8YbvVsrBvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(xses, ys)\n",
    "xspace = np.linspace(xses[0], xses[-1], 100)\n",
    "for sample in samples[5000::-100]:\n",
    "    plt.plot(xspace, polynomial(xspace, sample.variables['coefficients']), alpha=0.05, c='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Done!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}