{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" }, "colab": { "name": "MCMC.ipynb", "provenance": [], "collapsed_sections": [] } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "IzanQ5HN0xCX", "colab_type": "text" }, "source": [ " $\\quad \\quad \\quad$ *Credit: [Japan Times](https://www.japantimes.co.jp/news/2017/11/10/national/japans-casinos-work-keep-yakuza-deal-problem-drinking-experts/#.W-RJPJMzaUk).*\n" ] }, { "cell_type": "markdown", "metadata": { "id": "25VtwX9J0xCb", "colab_type": "text" }, "source": [ "# **Practical session: MCMC algorithms** (Geneviève Robin, Inria Paris)\n", "\n", "This practical session will consist of three parts. \n", "- In the first part, we experiment two Monte-Carlo Markov-Chain algorithms to sample from probability distributions with unknown normalizing constants: Randow-walk Metropolis-Hastings (RWMH) and the Metropolis adjusted Langevin algorithm (MALA). We will also study the discretization of the overdamped Langevin dynamics.\n", "- In the second part, we experiment how these algorithms can be used to compute empirical averages (empirical mean, etc.), and a Central Limit Theorem (CLT) for such averages. \n", "- Finally in the third part, we experiment the Hamiltonian dynamics and an MCMC algorithm based on it: the Generalized Hamiltonian Monte Carlo (GHMC) algorithm." ] }, { "cell_type": "markdown", "metadata": { "id": "kzjN3Ukr4I9r", "colab_type": "text" }, "source": [ "# 1. MCMC Sampling algorithms" ] }, { "cell_type": "markdown", "metadata": { "id": "-Wj-x3QLDVBT", "colab_type": "text" }, "source": [ "## 1.1. Random walk Metropolis-Hastings algorithm (RWMH)\n" ] }, { "cell_type": "markdown", "metadata": { "id": "A-PRnq88Dtyf", "colab_type": "text" }, "source": [ "The Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. For instance, when the normalizing constant of the probability distribution is intractable.\n", "\n", "The goal is to generate samples according to a probability measure $\\nu(q)\\propto \\exp(-\\beta V(q))$, where $q$ denotes the position, $\\beta = (k_B T)^{-1}$, and $V$ is a potential function. Recall that the Metropolis-Hastings algorithm produces a Markov Chain $q^1,\\ldots q^n, \\ldots$ iteratively as follows:\n", "\n", "\n", "* Given $q^n$, propose a position $\\tilde{q}^{n+1}$ according to transition probability $T(q^n, \\tilde{q}^{n+1})$\n", "* Accept the proposed position with probability $\\min(1,r(q^n, \\tilde{q}^{n+1}))$, where\n", "$$r(q^n, \\tilde{q}^{n+1}) = \\frac{\\nu(\\tilde{q}^{n+1})T(\\tilde{q}^{n+1}, q^n)}{\\nu({q}^{n})T(q^n,\\tilde{q}^{n+1})}.$$\n", "\n", "In the simplest setting, the proposals are generated through a random walk with Gaussian steps:\n", "$$\\tilde{q}^{n+1} = q^n + \\sigma G_n,$$\n", "where $G_n\\sim\\mathcal{N}(0, Id)$, and $\\sigma$ is the standard deviation of the Gaussian displacements. In this case, $T(\\tilde{q}^{n+1},q^n)=T(q^n,\\tilde{q}^{n+1})$, and thus we simply have $r(q^n, \\tilde{q}^{n+1}) = \\nu(\\tilde{q}^{n+1})/\\nu(q^n)$.\n", "\n", "\n", "Consider the double-well potential function defined by\n", "$$V(q) = (q^2-1)^2 + \\frac{1}{\\sqrt{2\\pi w}}\\exp\\left(-\\frac{(q-c)^2}{2w}\\right),$$\n", "$a\\geq 0$. The energy landscape associated with this potential is represented using the code below.\n" ] }, { "cell_type": "code", "metadata": { "id": "7sgqnIH_GlAd", "colab_type": "code", "outputId": "c2b60368-a332-4ca1-f36d-96711b2b526c", "colab": { "base_uri": "https://localhost:8080/", "height": 283 } }, "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import math\n", "\n", "center=0.2\n", "width=0.1\n", "magn=1\n", "\n", "def V(q):\n", " return 0.5*q**2 + magn*math.exp(-(q-center)**2/(2*width))/math.sqrt(2*math.pi*width)\n", "\n", "q=np.linspace(-2, 2,100)\n", "plt.plot(q, list(map(V, q)))\n", "\n" ], "execution_count": 42, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "[]" ] }, "metadata": { "tags": [] }, "execution_count": 42 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhUVZ7/8fc3+x6WLGQl7PsWAggi\niguiiKi44K6tg3a33T09PdOj0/OzZ+zp1ell7La1UXFXXFBABQFZBGVNMOxbCAlkIQkEQhayn98f\nKfqpjglUQlXdWr6v56knlXtv1f1UUXxz69xzzxFjDEoppXxXgNUBlFJKuZYWeqWU8nFa6JVSysdp\noVdKKR+nhV4ppXxckNUBOhIXF2cyMjKsjqGUUl4jJyfnpDEmvqN1HlnoMzIyyM7OtjqGUkp5DREp\n7GydNt0opZSP00KvlFI+Tgu9Ukr5OC30Sinl47TQK6WUj7tooReRNBFZJyL7RGSviPyog21ERJ4T\nkTwR2SUimXbrHhSRw7bbg85+AUoppS7Mke6VzcBPjDE7RCQayBGR1caYfXbb3AAMst0mAS8Ak0Sk\nF/BzIAswtscuM8acduqrUEop1amLHtEbY0qNMTts96uB/UBKu83mAG+YNluAHiKSBFwPrDbGVNqK\n+2pgplNfgU19UwsLNhzh67yTrnh6pZRyqbUHylj41VGaWlqd/txdaqMXkQxgHLC13aoU4Ljd70W2\nZZ0t7+i554tItohkV1RUdCUWAMGBASzYcJS3tnR6zYBSSnmsV78u4PXNBQQFiNOf2+FCLyJRwGLg\nn40xZ50dxBizwBiTZYzJio/v8CreCwoMEG4ancSaA+VU1zc5O55SSrlMRXUDX+edZPboZEQsKvQi\nEkxbkX/bGPNRB5sUA2l2v6falnW23CVmj0mmsbmV1fvKXLULpZRyuhV7Smk1cPPYZJc8vyO9bgR4\nBdhvjPlDJ5stAx6w9b65DKgyxpQCK4EZItJTRHoCM2zLXCIzvQcpPcL5ZGeJq3ahlFJO98nOEoYk\nRjM4Mdolz+9Ir5vLgfuB3SKSa1v2H0A6gDHmRWA5cCOQB9QBD9vWVYrIL4Dttsc9Y4ypdF78fyQi\n3DQmiVc2HuV0bSM9I0NctSullHKKkjPn2F5wmn+dMdhl+7hooTfGfAVcsNHItM0w/v1O1i0EFnYr\nXTfMHp3M377M5/O9J7h7Yrq7dquUUt3y2a5SAG4a7ZpmG/DBK2NHJMfQPz5Sm2+UUl5h2c4SxqTG\nkhEX6bJ9+FyhFxFmj05mc/4pys/WWx1HKaU6dfRkLbuLq5g9xnVH8+CDhR5g9pgkjIFPbV+JlFLK\nE31qa3mYNTrJpfvxyUI/MCGa4UkxfLJLm2+UUp7JGMPSnSVMzOhFUmy4S/flk4UeYM7YZL45doZj\np+qsjqKUUt+yv7SavPIa5oxzbbMN+HChP9/mtWyny67PUkqpbluaW0xQgHDjSNc224APF/rkHuFM\n7NeLJbkltPX+VEopz9Daali2s4QrB8e75Xofny300NZ8k1dew75Spw/No5RS3ba9oJLSqnqXDXnQ\nnk8X+htHJhEUICzL1ZOySinPsXRnCeHBgVw3PNEt+/PpQt8zMoQrB8ezbGcJra3afKOUsl5jcyvL\nd5cyY0QiESGOjEJz6Xy60EPbaHClVfVsL3DZEDtKKeWwjYcrOFPXxBw3NduAHxT664YnEh4cyBJt\nvlFKeYCluSX0jAjmikFdn3eju3y+0EeEBDFjRCLLd5fS0NxidRyllB+raWhm1b4T3DgqieBA95Vf\nny/0ALeOS6HqXBPrD3Z9ikKllHKWlXtOUN/Uym2ZHc6o6jJ+UeinDowjLiqEJd/oxVNKKessyS0m\nrVc4mek93bpfvyj0QYEBzB6TzJr95VSd0/lklVLuV3a2nq/zTnLr2BSXzAt7IX5R6KGt+aaxpZUV\nu3VES6WU+32ys4RWA3PGubfZBvyo0I9KiaV/fCQfa/ONUsoCS3KLGZMay4D4KLfv228KvYhw69gU\nth6tpPjMOavjKKX8yOGyavYUn+UWC47mwYFCLyILRaRcRPZ0sv7fRCTXdtsjIi0i0su2rkBEdtvW\nZTs7fFfNGdv2Ji/N1aN6pZT7LMktJjBAXDov7IU4ckT/GjCzs5XGmGeNMWONMWOBp4AvjTH2l6FO\nt63PurSoly69dwRZfXvy0Y5iHdFSKeUWra2GJd+UMHVgHPHRoZZkuGihN8ZsABwdP+Bu4N1LSuRi\nt2Wmkldew55iHdFSKeV655uL545PtSyD09roRSSCtiP/xXaLDbBKRHJEZP5FHj9fRLJFJLuiwnUX\nNs0alURIUACLdxS5bB9KKXXeRzuKiA4NYoabRqrsiDNPxs4Gvm7XbDPVGJMJ3AB8X0SmdfZgY8wC\nY0yWMSYrPt51Y0DERgRz7bAElu0soaml1WX7UUqpc40tLN9dyg2j+hAWHGhZDmcW+nm0a7YxxhTb\nfpYDHwMTnbi/brttXCqVtY18qUMiKKVcaNW+E9Q2tnBbpnXNNuCkQi8iscCVwFK7ZZEiEn3+PjAD\n6LDnjrtdOSSeXpEh2qdeKeVSi3cUk9IjnIkZvSzNcdFR70XkXeAqIE5EioCfA8EAxpgXbZvdCqwy\nxtTaPTQR+Nh2qW8Q8I4x5nPnRe++4MAAbh6TzDvbjlFV10RsRLDVkZRSPqbsbD1fHa7g+9MHEhDg\n3iEP2rtooTfG3O3ANq/R1g3Tflk+MKa7wVxtbmYqr20q4LPdpdwzKd3qOEopH7M0t5hW0zb8itX8\n5srY9kamxDAoIUp73yilnM4Yw+KcYsal96C/BUMetOe3hV5EmDs+lZzC0xw9WXvxByillIP2lpzl\nYFm15Sdhz/PbQg9tX6kCBBbn6FG9Usp5PswpIiQogJstGvKgPb8u9IkxYUwbHM/iHUW0tOqQCEqp\nS9fQ3MKS3GJmDE/0mI4efl3oAW4fn0ppVT2bjpy0OopSyges3V/OmbombrdwyIP2/L7QXzsskZiw\nID7U5hullBN8mFNEYkwoVwxy3RX+XeX3hT4sOJA5Y1P4fM8JztbrNINKqe4rr65n/aEKbstMJdDi\nvvP2/L7QQ1vzTUNzK5/t0mkGlVLdt/SbElpaDXM9pLfNeVrogdGpsQxKiOKD7ONWR1FKeSljDB/m\nFDEuvQcDE6zvO29PCz1tfervyEplx7Ez5JXXWB1HKeWFdhVVcbCs2qNOwp6nhd7m1nFtbWof5OhR\nvVKq697PPk5YcACzx3hG33l7Wuht4qNDuXpoAotzinWceqVUl5xrbGFZbgk3jkwiJswz+s7b00Jv\n586sNE7WNOg49UqpLvl8bynVDc3ckZVmdZQOaaG3c9WQeOKiQnlPT8oqpbrg/e1FpPeKYFI/a8ed\n74wWejvBgQHMzUxh7YFyyqvrrY6jlPICx07VsTn/FHdmpVo+7nxntNC3c0dWGi2thiU6+5RSygEf\n5hxHBOZ6YG+b87TQtzMwIYrxfXvy3vbjGKMDnSmlOtfS2tZ3ftqgeJJiw62O0ykt9B24MyuVIxW1\n5BSetjqKUsqDbThcQUlVPXd66EnY87TQd+Cm0clEhgSyaLuelFVKdW7RtmP0jgzhuuGJVke5oIsW\nehFZKCLlIrKnk/VXiUiViOTabk/brZspIgdFJE9EnnRmcFeKDA3i5rHJfLqrRAc6U0p1qLy6njX7\ny5k7PpWQIM8+ZnYk3WvAzItss9EYM9Z2ewZARAKB54EbgOHA3SIy/FLCutO8CenUN7WyLLfE6ihK\nKQ+0OKeY5lbj8c024EChN8ZsACq78dwTgTxjTL4xphFYBMzpxvNYYnRqLMOSYli0/ZjVUZRSHsYY\nw3vbjzExo5fHDWDWEWd935gsIjtFZIWIjLAtSwHsG7mLbMs6JCLzRSRbRLIrKqy/MlVEmDchjT3F\nZ9lTXGV1HKWUB9mSX0nBqTrmTfT8o3lwTqHfAfQ1xowB/gws6c6TGGMWGGOyjDFZ8fGeMTPLLWNT\nCA0K0KN6pdQ/WLT9GNFhQdwwMsnqKA655EJvjDlrjKmx3V8OBItIHFAM2P+5S7Ut8xqxEcHMGpXE\n0m9KqGtstjqOUsoDnKlrZMWeE9w6LoXwkECr4zjkkgu9iPQREbHdn2h7zlPAdmCQiPQTkRBgHrDs\nUvfnbvMmplPd0MynO3X2KaUULN5RTGNzK/MmpFsdxWFBF9tARN4FrgLiRKQI+DkQDGCMeRG4Hfiu\niDQD54B5pu2S0mYReQJYCQQCC40xe13yKlxoQkZPBiZE8fa2Y9w5wTva45RSrmGM4Z2thYxN68Hw\n5Bir4zjsooXeGHP3Rdb/BfhLJ+uWA8u7F80ziAj3TEznmU/3sbekihHJsVZHUkpZZNvRSo5U1PK7\n20dbHaVLPLuXv4eYm5lKaFAA72zVk7JK+bN3trWdhJ092vNmkboQLfQOiI0IZtboJJbmllDboCdl\nlfJHlbWNrNh9gtu86CTseVroHXTvpHRqGpr5ZKdeKauUP1qcU0RjSyt3T/Kek7DnaaF3UGZ6T4Yk\nRvPONm2+UcrfGGN4d9sxMtN7MLSP95yEPU8LvYNEhHsmpbOrqIrdRXqlrFL+ZHP+KfJP1nLPpL5W\nR+kWLfRdcGtmCuHBgby1pdDqKEopN3p7yzFiw4O5abR3XAnbnhb6LogJC+aWccks3VlMVZ0OX6yU\nPyg/W8/KvSe4MyuVsGDvOgl7nhb6Lrrvsr7UN7WyeEeR1VGUUm6waPtxmluN1zbbgBb6LhuRHMu4\n9B68taVQ55RVysc1t7TyztZjXDEojn5xkVbH6TYt9N1w36S+5J+sZdORU1ZHUUq50JoD5Zw4W899\nl3nv0Txooe+WWaOT6BERrCdllfJxb20pJCk2jGuGJlgd5ZJooe+GsOBA7sxKY9W+Mk5U1VsdRynl\nAkdP1rLx8EnumZhOUKB3l0rvTm+heyel02qMXkCllI96c3MhQQHCXT4waq0W+m7q2zuS6UMSeGfr\nMRqbW62Oo5RyotqGZj7IOc6No5JIiAmzOs4l00J/CR6Y3JeTNQ2s2KOTkijlS5bkFlNd38yDU7z7\nJOx5WugvwbRB8WT0juCNzXpSVilfYYzhjU2FjEiOITO9p9VxnEIL/SUICBDun5xBTuFp9hTr+DdK\n+YKtRys5WFbNg1MysM2S6vW00F+i28enEh4cyBubC6yOojxYQ3MLW/JP8X9fHOb97cc519hidSTV\nidc3FdAjIpibx3jX5CIXctGpBNWFxYYHc2tmCotzinjqhmH0jAyxOpLyIDUNzTy5eBer95XRYHfS\n/pfL93PXhDQevjyDpNhwCxMqe6VV51i1r4xHr+jntePadOSiR/QislBEykVkTyfr7xWRXSKyW0Q2\nicgYu3UFtuW5IpLtzOCe5MHJGTQ0t7Jo+3GroygPcqaukXtf3sqKPSeYNyGNlx/IYufPZ/De/MuY\nOjCOV746yi3Pf83xyjqroyqbNze3DW1ynxePa9MRR5puXgNmXmD9UeBKY8wo4BfAgnbrpxtjxhpj\nsroX0fMN6RPNlAG9eXNzAc0t2tVStY14eNfftrC/9Cwv3jee/54zkmuHJxIbHsyk/r15/t5MPv3B\nVOqbWrn35a2UndUL76xW39TCu9uOcd3wRNJ6RVgdx6kuWuiNMRuAygus32SMOW37dQuQ6qRsXuXh\ny/tRUlXPyr1lVkdRFqttaGbegi0cP13Hqw9N4LrhiR1uNywphte/M5FTNQ3c9/JWKmsb3ZxU2Vvy\nTTGn65p4+PJ+VkdxOmefjH0EWGH3uwFWiUiOiMy/0ANFZL6IZItIdkVFhZNjud7VQxNI6xXOq18f\ntTqKstjvVx0i/2QtLz+YxeUD4y647di0Hrz84AQKK+t4/K0cWlt1RFQrGGN49esChiXFMKlfL6vj\nOJ3TCr2ITKet0P+73eKpxphM4Abg+yIyrbPHG2MWGGOyjDFZ8fHxzorlNoEBwoOTM8guPK1TDfqx\nXUVneG3TUe67LJ0pAy5c5M+bPKA3/zNnJNuOVvJhjs5zYIXNR05xsKyahy/3nS6V9pxS6EVkNPAy\nMMcY8/exe40xxbaf5cDHwERn7M9T3TkhjciQQF7dpEf1/qi5pZUnF+8mLiqUn84c2qXH3j4+lQkZ\nPfn1iv2c1iYct1v4dQG9I0N8qkulvUsu9CKSDnwE3G+MOWS3PFJEos/fB2YAHfbc8RUxYcHcPj6V\nT3eWUlHdYHUc5Wavfl3AvtKz/PfNI4gJC+7SYwMChF/cMpKz9c389vMDLkqoOnLsVB1rDpRxz6R0\nn+pSac+R7pXvApuBISJSJCKPiMjjIvK4bZOngd7AX9t1o0wEvhKRncA24DNjzOcueA0e5cEpGTS2\ntPKmjlXvVyqqG/jD6kNcOyyBmSP7dOs5hvaJ4ZGp/Vi0/Tg5hZ32f1BO9tqmAgJFvH5ykQu56AVT\nxpi7L7L+UeDRDpbnA2O+/Qjf1j8+imuHJfD2lkK+d9UAnz1CUP/otU1HqW9u4Wezhl9SG++PrhnE\nJztLeHrpXj79wVSfbC/2JGfrm3hv+zFuHpNMog+MUtkZHQLBBb4ztR+nahtZ8k2x1VGUG9Q0NPPm\n5kJuGNnnkucVjQwN4iczhrC35CzrD3pf7zNv896249Q2tvCdqb7XpdKeFnoXmNy/N8OTYnj5q6M6\ngbgfWLTtGGfrm3ls2gCnPN+cscmk9Ajn+XV5Tnk+1bHmllZe/fool/XvxciUWKvjuJQWehcQER69\noh955TV8eUiPynxZY3Mrr3zVVizGpPVwynMGBwYwf1p/sgtPs+2ottW7yoo9JyipqufRqf2tjuJy\nWuhd5KbRySREh/LKV9rV0pd9srOE0qp6HrvSOUfz592ZlUbvyBD+ul6P6l3BGMPLG/PpFxfJ1V4+\n8bcjtNC7SEhQAA9OyWDj4ZMcOHHW6jjKBYwx/G3DEYYkRnPVYOde5BceEsh3pvZj/cEKnevABXIK\nT7OzqIrvTO1HQIDvn/DWQu9C90xMJyw4gJc36lG9L9p85BSHymqYP62/S3rH3D+5L9GhQbzw5RGn\nP7e/W7Ahnx4RwczNTLE6iltooXehnpEh3JmVxtLcYk5U6eiEvubDnCKiw4KYNTrJJc8fExbMvZf1\nZcXuUkqrzrlkH/7oSEUNq/eXcf9lfYkI8Y8pObTQu9ijU/vT0mp0WAQfU13fxPI9pcwek+zSayXu\nnZSOARZt07kOnOXljfkEBwbwwOQMq6O4jRZ6F0vvHcENI5N4Z8sxquubrI6jnGT57lLqm1q5Y7xr\nR+VO6xXBlYPjWbT9mM514AQV1Q0s3lHM3MxU4qNDrY7jNlro3WD+tP5UNzTzns5A5TM+yC5iQHwk\nY53UpfJC7p3Ul7KzDaw5UO7yffm61zcV0NTSyj9d4dsXSLWnhd4NxqT1YFK/Xiz86ihNelTm9Y6e\nrCW78DR3ZKW5ZYiC6UPiSYoN4+2tx1y+L19W29DMm1sKmTE8kf7xUVbHcSst9G7y2JX9Kamq59Nd\nJVZHUZfow5zjBAjcOs49PTaCAgOYNyGdDYcqOHZK55ftrve2H6fqXBPznXQFszfRQu8mVw1OYHBi\nFC+uz9dhEbxYS6vhox3FXDk43q2DYN01IY3AAOGdbXpU3x2Nza28vDGfCRk9Gd+3p9Vx3E4LvZsE\nBAiPXzmAg2XVrDuoba3eakv+KUqr6pnr4pOw7fWJDePaYQl8kH2cxmZt/uuqpbnFlFTV872rBlod\nxRJa6N1o9pi2wapeWK8XwHirz3aXEhESyLXDOp7w25XumpDGqdpGPVDootZWw4tfHmFon2iuGuJ9\n05Q6gxZ6NwoODOCfrujH9oLTbC/Qwaq8TXNLKyv3nODqoQmWzDMwbVA8cVGhLNZ5Zbtk1b4yjlTU\n8t2rBvjt+P5a6N3srgnp9IoM0aN6L7TtaCWnahuZNco1V8JeTFBgALeMTWbdwXIqdV5ZhxhjeOHL\nI6T3irDs380TaKF3s/CQQB6aksHaA+XsL9XBzrzJZ7tLCQ8O5Koh1o12OHd8Kk0thmW5OqmNIzYf\nOcXO42eYP60/QYH+W+7895Vb6IHJfYkMCdSjei/S0mpYubet2SY8xLrpIYclxTA8KYbFO7TQO+L5\n9XnERYVyu5tPnnsahwq9iCwUkXIR2dPJehGR50QkT0R2iUim3boHReSw7fags4J7sx4RIdx7WV8+\n3VVCwclaq+MoB2w7WsnJmkZu9ICv/3PHp7K7uIpDZdVWR/FoOYWn+TrvFPOn9fP7uZsdPaJ/DZh5\ngfU3AINst/nACwAi0gv4OTAJmAj8XET8rxNrBx69oh9BgQF6VO8llu8uJSw4gOlDre+1MWdsMkEB\nwuIdelL2Qv6y9jA9I4K5d1Jfq6NYzqFCb4zZAFyom8gc4A3TZgvQQ0SSgOuB1caYSmPMaWA1F/6D\n4TcSosO4e0Iai3cUUXxGh6D1ZC2thhV7TjB9SIJHDGsbFxXKVUPiWfJNMS2tevFdR/YUV7HuYAWP\nTO1HZKj1/2ZWc1YbfQpgP2JXkW1ZZ8u/RUTmi0i2iGRXVPjHPKvzbdPPLdCJJTxadkElJ2saPKLZ\n5rzbMlMpO9vApiMnrY7ikf6yNo/osCAemJJhdRSP4DEnY40xC4wxWcaYrPh4678eu0NKj3DmZqby\n7vbjlFfrxCSeavW+MkICA5juQXOLXj00gejQIJZ8o2MntXfwRDWf7z3Bw1MyiAkLtjqOR3BWoS8G\n0ux+T7Ut62y5svnuVQNobmnV6QY9lDGG1fvLmDKwN1Ee1AQQFhzIjaOS+HxPKecaW6yO41H+si6P\niJBAHr7cv4YivhBnFfplwAO23jeXAVXGmFJgJTBDRHraTsLOsC1TNhlxkcwek8ybmws5VdNgdRzV\nTl55DYWn6iwZ8uBi5oxLpraxhdX7y6yO4jEOl1Xz6a4SHpicQc/IEKvjeAxHu1e+C2wGhohIkYg8\nIiKPi8jjtk2WA/lAHvAS8D0AY0wl8Atgu+32jG2ZsvODqwdR39zCS3pU73FW7Wsrop5Y6C/r15uk\n2DCWfqNfks97bm0e4cGBzJ/W3+ooHsWh76LGmLsvst4A3+9k3UJgYdej+Y+BCVHMHp3MG5sLmD+t\nP730SMRjfLG/jNGpsfSJdd+QxI4KCBBuHpvMKxuPUlnb6Pefm0O2o/nHrxzg9+9Fex5zMtbf/fCa\ngZxrauGljflWR1E25dX15B4/45FH8+fdMjaF5lbDZzqhDc+tOUxEcCDzr9Cj+fa00HuIgQnRzBqV\nxBubCjitA1Z5hLX7yzEGrhvuuYV+WFIMQ/tE87GfN98cKqvms92lPHS5ts13RAu9B/nhNYOoa2rh\n5a/0qN4TfLG/jJQe4QztE211lAu6ZVwKO46dofCU/w6n8X9fHCYyJIhHp+rRfEe00HuQwYnR3Dgq\nide+LtBhaC1W19jMxsMnuW54osePYX7zmGRE8Ns+9XtLqvhsdykP69F8p7TQe5gfXzuIc00t/E2v\nlrXUV4dP0tDc6tHNNucl9wjnsn69+fibIr+cj/j3qw4RGx7Mo9o23ykt9B5mYEI0c8am8PrmAr1a\n1kJrD5QTHRrEhIxeVkdxyK2ZKRScqiP3+Bmro7hVTmElaw+U89iV/YkN16tgO6OF3gP96JpBNLUY\nHdnSIq2thrUHypk2OJ6QIO/4L3LDyD6EBgX41UlZYwzPrjxIXFQoD+mYNhfkHZ9iP5MRF8nczBTe\n3nqM0iod2dLd9pacpby6gas9aGybi4kOC+a64Yl8srOEppZWq+O4xdd5p9iSX8n3pw/wiFFFPZkW\neg/1g6sHYYzh+XV5VkfxO2sOlCECVw3xrsH1bh2Xwum6Jr486PujvxpjeHbVQZJjw7hnUrrVcTye\nFnoPldYrgrsmpLFo23GOnaqzOo5fWXugnHFpPegdFWp1lC6ZNjieXpEhfOwH88mu2HOCncfP8M/X\nDSY0yL9nj3KEFnoP9oOrBxEUKPzxi0NWR/Eb5dX17Cqq4hoPvhq2M8GBAcwencTqfWWcrW+yOo7L\nNLW08uzKgwxJjGZupn/PBesoLfQeLDEmjIem9GNJbjEHTpy1Oo5fWH+grdlj+hDvaZ+3d2tmKo3N\nrXy2q9TqKC6zaNsxjp6s5d9vGEJggGdf4+AptNB7uO9eOYCo0CD+d+VBq6P4hTUHykiKDWNYkmdf\nDduZMamxDEqI4oPs4xff2AvVNDTzf2sOM6lfL6/9Y2wFLfQeLjYimMevHMAX+8vJKdQRnl2pobmF\njYdPcvXQBI+/GrYzIsIdWansOHaGvPJqq+M43Usb8jlZ08hTNw7z2n8jK2ih9wIPX55BXFQov/38\noF9e+eguW/MrqWts4Zph3n2keMu4FAIDhA9yiqyO4lTlZ+t5aWM+s0YlMTath9VxvIoWei8QERLE\nj64ZyLajlaw7WG51HJ+19kA5oUEBTO4fZ3WUS5IQHcb0IfF8tKOYZh/qU//syoM0txj+feZQq6N4\nHS30XmLexHT6xUXymxUHaGnVo3pnM8aw5kAZUwfGER7i/d317shKo6K6gQ2HfaNP/e6iKj7cUcTD\nUzNI7x1hdRyvo4XeSwQHBvDT64dwqKyGxT72ldwT5JXXcLzyHFd7ebPNeVcPTaB3ZAjvb/f+z4ox\nhl98uo9eESE8MX2g1XG8kqNzxs4UkYMikiciT3aw/o8ikmu7HRKRM3brWuzWLXNmeH8zc2QfxqX3\n4PerD3KuscXqOD5lzYG2JjFvGvbgQoIDA7hlXAprDpR5/ZDXn+85wbaCSn4yYwjRYTpwWXdctNCL\nSCDwPHADMBy4W0SG229jjPmxMWasMWYs8GfgI7vV586vM8bc7MTsfkdEeOqGYZSdbWDh1zqRuDOt\n3V/O8KQYkmLDrY7iNHdmpdHUYvgwx3u7WtY3tfCrFfsZ2ieauyakWR3HazlyRD8RyDPG5BtjGoFF\nwJwLbH838K4zwqlvm9ivF9cOS+SF9Uc4WdNgdRyfcKaukezCSq/vbdPekD7RTMjoydtbj9Hqped1\n/vZlPscrz/H0TcP14qhL4EihTwHsDwmKbMu+RUT6Av2AtXaLw0QkW0S2iMgt3U6q/u6pG4dS39TC\nn3RoBKdYf7CCVoNXDntwMbTJuwsAABWzSURBVPdPzqDwVB1feuFJ2WOn6vjr+jxmj0lmykDv7gll\nNWefjJ0HfGiMsW9A7muMyQLuAf4kIgM6eqCIzLf9QciuqPC+D6U7DYiP4r7L+vLO1mMcKvO9i2Lc\nbc2BcuKiQhidEmt1FKebOaIPcVGhvLW50OooXfbMp3sJChB+duMwq6N4PUcKfTFg3ziWalvWkXm0\na7YxxhTbfuYD64FxHT3QGLPAGJNljMmKj/eu4WGt8KNrBhEVGsSvlu+3OopXa2pp5cuD5UwfkkCA\nDzYNhAQFcPfENNYeLOd4pfeMgvrFvjK+2F/Oj64dRJ/YMKvjeD1HCv12YJCI9BORENqK+bd6z4jI\nUKAnsNluWU8RCbXdjwMuB/Y5I7i/6xkZwg+vGcT6gxV8eUi/AXVXTuFpztY3+1z7vL17JqUTIMJb\nW73jqP5cYwv//eleBiVE8fDl/ayO4xMuWuiNMc3AE8BKYD/wvjFmr4g8IyL2vWjmAYvMP16jPwzI\nFpGdwDrgN8YYLfROcv/kvvTtHcEvP9vnU1dAutOa/WUEBwpTB/nut8ik2HCuG5bI+9uPU9/k+d1y\n//TFIY5XnuOZOSMJDtRLfZzBoXfRGLPcGDPYGDPAGPNL27KnjTHL7Lb5L2PMk+0et8kYM8oYM8b2\n8xXnxvdvoUGBPHXDMA6V1fDOtmNWx/E6xhhW7ytj8oA4okJ9eyq6+yf35XRdE8tyS6yOckF7iqt4\naWM+8yakMXlAb6vj+Az9c+nlrh+RyOUDe/P7VYc47eUXxrhbXnkNBafquG647/W2aW/KgN4MT4rh\nxQ1HPHYIjaaWVn764S56R4XylJ6AdSot9F5ORPj57BHUNDTz+9U6Zn1XrNpXBsB1Ptitsj0R4fvT\nB5JfUcvKvSesjtOhlzceZV/pWX4xZwSx4XoFrDNpofcBgxOjud/W3XJfic5E5ajV+8oYnRrrN706\nZo7sQ/+4SJ5fl+dxw13nldfwpy8OMXNEH2aOTLI6js/RQu8jfnztYGLDg/mvT/Z63H9iT1R+tp7c\n42f84mj+vMAA4fGrBrC35KxH9dRqamnlx+/lEhESyDNzRlgdxydpofcRsRHB/Nv1Q9l2tJJlOz37\nhJsn+GJ/2yBm143wn0IPcMvYFJJjw3h+XZ7VUf7uuTWH2V1cxa9vG0VCjH98u3I3LfQ+5K4JaYxO\njeWXn+2nur7J6jgebfW+E6T1CmdIonfODdtdIUEBzJ/Wn+0Fp9maf8rqOOQUVvL8ujxuH5+qTTYu\npIXehwQGCL+YM5KKmgb+uPqw1XE8Vk1DM1/nnWLG8D5+Oe/oXRPSSYgO5dcrDlg62FlNQzM/fm8n\nyT3C+fns4Rd/gOo2LfQ+ZkxaD+6emM7rmwvYX6onZjuy4VAFjS2tftGtsiPhIYH82/VDyD1+hqU7\nOxvNxLWMMfz7h7soOl3HH+4cq+PMu5gWeh/00+uHEBsezNNL9+iJ2Q6s2nuCHhHBZPXtaXUUy8zN\nTGV0aiy/XXGQusZmt+//la+O8tnuUn46cygT+/Vy+/79jRZ6H9QjIoQnZw5le8FpPsj2/qnknKm+\nqYUv9pczY3giQX58eX1AgPD0TcM5cbaeF7/Md+u+t+af4tcrDnD9iEQem9bfrfv2V/77Sfdxt49P\nZWJGL365fL9OUGJn4+GT1DQ0c+MoPfGXldGL2WOS+duXRyg+c84t+zxRVc8T735Deq8Inr1jjF+e\nI7GCFnofFRAg/Oq2kdQ1NvM/n+o4cuct311KbHgwl+tEFgA8ecNQROCpj3a7/MRs1bkmHnp1G3UN\nzbxwXyYx2i7vNlrofdjAhGi+e+UAluSWsNELZxhytobmFr7YV8aM4Yk6KqJNSo9w/nPWcDYcquDl\nr1zXhNPQ3MJjb2aTV17Di/ePZ2ifGJftS32bftp93PemD6RfXCQ/+3gP5xo9f4haV9p46CTVDc3c\nOFqbbezdOymdG0b24XefHyT3+BmnP39rq+Ff3tvJlvxK/veOMVzhw0NCeyot9D4uLDiQX946kmOV\ndX4/x+zy3aXEhAVx+QBttrEnIvzmttEkxoTxg3d3cNaJF9s1tbTykw928tnuUv7jxqHcMq7D6aaV\ni2mh9wNTBsQxb0IaL23MZ6cLjti8QUNzC6v3lzFjRB9CgvRj315sRDDP3T2OkjP1fP/tHU759neu\nsYX5b2Tz8TfF/OS6wfzTFdrDxir6ifcT/zFrGAnRYfz0w100NvvfbFRf552kur6ZWdrbplPj+/bk\nN7eN4qu8k3znte3UNnS/f/3p2kbue2Ur6w9V8MtbR/KDawZpDxsLaaH3EzFhwfzy1pEcLKv2qAGt\n3OXTXbZmG+1tc0F3ZKXxxzvHsvXoKR5cuK1bYyatO1jO9X/awO6iKp6/J5N7J/V1QVLVFVro/cg1\nwxK5ZWwyz6/L86vhEeoam1m55wQ3jEzSZhsH3DIuhT/fnUnu8TPc9Oev+GJfmUNXWFeda+Kpj3bz\n8Kvb6RkRwkffm6LXK3gIhz71IjJTRA6KSJ6IPNnB+odEpEJEcm23R+3WPSgih223B50ZXnXd07NH\n0CMimH95f6ffNOGs3HuC2sYW5o5PtTqK15g1Ook3H5lEUIDw6BvZPPzadnIKT3f4mdlfepanPtrN\nZb9aw6Ltx3jsyv4s+8HljEyJtSC56ohc7C+1iAQCh4DrgCJgO3C3MWaf3TYPAVnGmCfaPbYXkA1k\nAQbIAcYbY05faJ9ZWVkmOzu7yy9GOWb1vjL+6Y1snpg+kH+9fojVcVzu/le2UnCqli//dToBAdpO\n3BVNLa28vqmA//viMNUNzYQGBTAqJZb46FBKquopOXOOiuoGQoMCuHlMMg9dnsGIZC3wVhCRHGNM\nVkfrghx4/EQgzxiTb3uyRcAcwJHLLa8HVhtjKm2PXQ3MBN51JLhyjeuGJ3L7+FT+uj6Pa4YlMC7d\ndwf3Kq06x1d5J/nh1YO0yHdDcGAAj17Rn7mZqWw6coodx06z49hpDpVVk9wjnKFDEhiaFM2t41Lo\nERFidVzVCUcKfQpw3O73ImBSB9vNFZFptB39/9gYc7yTx3bYkVZE5gPzAdLT0x2IpS7F07OHs/nI\nKX7ywU6W//AKwoIDrY7kEku+KcEYuC1T+29fip6RIcwancQsvdjMKznrzNQnQIYxZjSwGni9q09g\njFlgjMkyxmTFx+uVc64WExbMs7ePJr+ill8t3291HJcwxvDRjiKy+vakb+9Iq+MoZRlHCn0xkGb3\ne6pt2d8ZY04ZY84PkfgyMN7RxyrrTBkYxyNT+/HG5kLW7C+zOo7T7S6u4nB5Dbdl6klY5d8cKfTb\ngUEi0k9EQoB5wDL7DUTE/vvczcD5Q8SVwAwR6SkiPYEZtmXKQ/x05hCGJcXwbx/uory63uo4TvXR\njmJCggK0uUH5vYsWemNMM/AEbQV6P/C+MWaviDwjIjfbNvuhiOwVkZ3AD4GHbI+tBH5B2x+L7cAz\n50/MKs8QGhTIc/PGUtvQzE/e32npHKLOVNfYzEc7irh+RB9iw3U4XOXfHDkZizFmObC83bKn7e4/\nBTzVyWMXAgsvIaNysUGJ0fy/m4bzn0v28MpXR/knH5j1Z2luCWfrm3lgsl6VqZReJqiAtqFqrx+R\nyG8/P0BO4QUvc/B4xhhe31TAsKQYv54XVqnztNAroG2o2t/dPoakHmH84J0dnK5ttDpSt20vOM2B\nE9U8OLmvDqSlFFrolZ3Y8GD+es94TtY08i/v53pte/3rmwuIDQ9mzljtO68UaKFX7YxKjeU/bxrG\nuoMVvPDlEavjdNmJqno+33OCuyakER7imxeBKdVVWujVt9x/WV9mj0nmf1cdZP3BcqvjdMk7Wwtp\nNYb7dGhcpf5OC736FhHht3NHMSQxmh+++w0FJ2utjuSQ2oZm3tp6jKuHJJDeO8LqOEp5DC30qkMR\nIUG89EAWAQHC/DezL2m2IXd5Y3MhlbWNfP/qgVZHUcqjaKFXnUrrFcGf7x5HXnkNP34vlxYPPjlb\n09DMgg1HuHJwPJk+PBqnUt2hhV5d0BWD4vnZrOGs2lfGb1Z47uBnr28q4HRdEz++brDVUZTyOA5d\nGav823cuz6DwVC0vbTxKeu9I7r/Ms050Vtc38dLGfKYPiWdsWg+r4yjlcbTQq4sSEZ6+aThFp8/x\n86V7SO0ZzvQhCVbH+rvXNxVwpq6Jf75Wj+aV6og23SiHBAUG8Oe7xzG0Twzfe2sHOYWeMTZdeXU9\nf9uQzzVDExijR/NKdUgLvXJYZGgQr39nIokxoTz06nb2lZy1OhL/8+l+Gppa+Y9Zw6yOopTH0kKv\nuiQ+OpS3Hp1EVGgQDyzcylEL+9h/eaiCZTtL+N70AQyIj7Ish1KeTgu96rLUnhG8+cgkjIG7F2zh\nSEWN2zOca2zhP5fspn98JN+9aoDb96+UN9FCr7plYEIUbz06iebWVu7622YOnHBvM85zaw9zvPIc\nv7p1FKFBOqaNUheihV5127CkGBbNn0xggDBvwRZ2FZ1xy343Hq5gwYZ87sxK5bL+vd2yT6W8mRZ6\ndUkGJkTxwWNTiAoNYt6CLazae8Kl+8uvqOH7b+9gUEIUT88e4dJ9KeUrtNCrS5beO4LF353CoIQo\nHnsrhxfWH8EY5w+XUFXXxKOvZxMUGMBLD2QRFaqXgSjlCIcKvYjMFJGDIpInIk92sP5fRGSfiOwS\nkTUi0tduXYuI5Npuy5wZXnmOxJgw3ntsMjeNTua3nx/gx+/lUl3f5LTnr29q4Yl3d3D8dB0v3JtJ\nWi8dnVIpR130kEhEAoHngeuAImC7iCwzxuyz2+wbIMsYUyci3wV+B9xlW3fOGDPWybmVBwoLDuS5\neWMZkhjFH1YfYnvBaZ69YzRTBsRd0vOWV9cz/40cdhad4bdzRzNJ2+WV6hJHjugnAnnGmHxjTCOw\nCJhjv4ExZp0xps726xYg1bkxlbcQEZ64ehAffncKIUEB3PPSVp5euodTNQ3der79pWe55S9fc/BE\nNS/cO547s9KcnFgp3+dIoU8Bjtv9XmRb1plHgBV2v4eJSLaIbBGRWzp7kIjMt22XXVFR4UAs5cky\n03uy/IdX8NCUDN7cUsgVv1vHsysPcKbOsUnHT9c28rvPD3DbXzfRauCDxyczc2QfF6dWyjc59WyW\niNwHZAFX2i3ua4wpFpH+wFoR2W2M+dZkpMaYBcACgKysLM8d+Fw5LDwkkP+6eQT3XZbOH784zPPr\njrDwqwKmDY7juuF9uGpIPL0jQxARAKrONbG3pIqNh0/y5uZCahubmTUqif9303ASY8IsfjVKeS9H\nCn0xYP99OdW27B+IyLXAz4ArjTF//55ujCm2/cwXkfXAOMD7Zp1W3TYwIZrn78nkielneXtrIV/s\nK2fl3jIAggKEHhHBhAQGUFJV//fHzBqVxI+uHcTgxGirYivlM+Ri3eBEJAg4BFxDW4HfDtxjjNlr\nt8044ENgpjHmsN3ynkCdMaZBROKAzcCcdidyvyUrK8tkZ2d38yUpT2eMYXdxFduOVlJZ28jpuibO\nNTYzKDGakSmxjEyOoXdUqNUxlfIqIpJjjMnqaN1Fj+iNMc0i8gSwEggEFhpj9orIM0C2MWYZ8CwQ\nBXxg+xp+zBhzMzAM+JuItNJ2PuA3FyvyyveJCKNTezA6VYcVVsodLnpEbwU9oldKqa650BG9Xhmr\nlFI+Tgu9Ukr5OC30Sinl47TQK6WUj9NCr5RSPk4LvVJK+Tgt9Eop5eM8sh+9iFQAhd18eBxw0olx\nnEVzdY3m6hrN1TW+mKuvMSa+oxUeWegvhYhkd3bRgJU0V9dorq7RXF3jb7m06UYppXycFnqllPJx\nvljoF1gdoBOaq2s0V9dorq7xq1w+10avlFLqH/niEb1SSik7WuiVUsrHeX2hF5FnReSAiOwSkY9F\npMPZLERkpogcFJE8EXnSDbnuEJG9ItIqIp12lxKRAhHZLSK5IuLyQfi7kMvd71cvEVktIodtP3t2\nsl2L7b3KFZFlLsxzwdcvIqEi8p5t/VYRyXBVli7mekhEKuzeo0fdkGmhiJSLyJ5O1ouIPGfLvEtE\nMl2dycFcV4lIld179bSbcqWJyDoR2Wf7v/ijDrZx7ntmjPHqGzADCLLd/y3w2w62CaRtntr+QAiw\nExju4lzDgCHAeiDrAtsVAHFufL8umsui9+t3wJO2+0929O9oW1fjhvfooq8f+B7wou3+POA9D8n1\nEPAXd32ebPucBmQCezpZfyOwAhDgMmCrh+S6CvjUne+Vbb9JQKbtfjRtU7W2/3d06nvm9Uf0xphV\nxphm269baJu8vL2JQJ4xJt8Y0wgsAua4ONd+Y8xBV+6jOxzM5fb3y/b8r9vuvw7c4uL9XYgjr98+\n74fANWKbR9PiXG5njNkAVF5gkznAG6bNFqCHiCR5QC5LGGNKjTE7bPergf1ASrvNnPqeeX2hb+c7\ntP0VbC8FOG73exHffmOtYoBVIpIjIvOtDmNjxfuVaIwptd0/ASR2sl2YiGSLyBYRcdUfA0de/9+3\nsR1oVAG9XZSnK7kA5tq+7n8oImkuzuQIT/7/N1lEdorIChEZ4e6d25r8xgFb261y6nt20cnBPYGI\nfAH06WDVz4wxS23b/AxoBt72pFwOmGqMKRaRBGC1iBywHYlYncvpLpTL/hdjjBGRzvr99rW9X/2B\ntSKy2xhzxNlZvdgnwLvGmAYReYy2bx1XW5zJU+2g7fNUIyI3AkuAQe7auYhEAYuBfzbGnHXlvryi\n0Btjrr3QehF5CLgJuMbYGrjaKQbsj2xSbctcmsvB5yi2/SwXkY9p+3p+SYXeCbnc/n6JSJmIJBlj\nSm1fUcs7eY7z71e+iKyn7WjI2YXekdd/fpsiEQkCYoFTTs7R5VzGGPsML9N27sNqLvk8XSr74mqM\nWS4ifxWROGOMywc7E5Fg2or828aYjzrYxKnvmdc33YjITOCnwM3GmLpONtsODBKRfiISQtvJM5f1\n2HCUiESKSPT5+7SdWO6wh4CbWfF+LQMetN1/EPjWNw8R6Skiobb7ccDlwD4XZHHk9dvnvR1Y28lB\nhltztWvHvZm29l+rLQMesPUkuQyosmums4yI9Dl/XkVEJtJWD139xxrbPl8B9htj/tDJZs59z9x9\nxtnZNyCPtrasXNvtfE+IZGC53XY30nZ2+whtTRiuznUrbe1qDUAZsLJ9Ltp6T+y03fZ6Si6L3q/e\nwBrgMPAF0Mu2PAt42XZ/CrDb9n7tBh5xYZ5vvX7gGdoOKADCgA9sn79tQH9Xv0cO5vq17bO0E1gH\nDHVDpneBUqDJ9tl6BHgceNy2XoDnbZl3c4FeaG7O9YTde7UFmOKmXFNpOze3y65u3ejK90yHQFBK\nKR/n9U03SimlLkwLvVJK+Tgt9Eop5eO00CullI/TQq+UUj5OC71SSvk4LfRKKeXj/j/9vmb9hp6T\nJgAAAABJRU5ErkJggg==\n", "text/plain": [ "