{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Monte Carlo (MC) Methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we will become familiar with the basic notion of Mote-Carlo methods:\n", "\n", "1. Sampling probability distributions\n", "2. MC estimate for $\\pi$\n", "3. Metropolis Hastings algorithm for the Ising model on a square lattice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sampling Custom Probability Distributions\n", "\n", "Our first goal is to learn how to sample probability distributions using numpy. \n", "\n", "First, we define a Gaussian mixture in $1d$, and we visualize it. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "import numpy as np # load numpy library\n", "import matplotlib.pyplot as plt # load plot library\n", "\n", "# set the seed of the random number generator (required for reproducibility)\n", "seed=0\n", "np.random.seed(seed)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2oklEQVR4nO3deXxcZ3Xw8d+Z0b6vlrVbjpd4j2PHCUnIQkhIAolbCDShFEoDaQph6waU96UFutC+bWlKgZBC2MpSoEASEgg0ITtxbMervCqyrdWWZFm7tczMef+YGUdRRtbVcueONOf7+Uykmbn3zokkz5nnOc8iqooxxhgzkc/rAIwxxiQmSxDGGGNisgRhjDEmJksQxhhjYrIEYYwxJqYUrwOYSyUlJbpkyRKvwzDGmHlj586dXapaGuu5BZUglixZwo4dO7wOwxhj5g0ROTHZc9bFZIwxJiZLEMYYY2KyBGGMMSYmSxDGGGNisgRhjDEmJksQxhhjYrIEYYwxJiZLEMaYWQmFlEf2tvNcQxe2fcDCsqAmyhlj4uvlzgE+/uO97DhxBoBLlhTysetXcPkFJR5HZuaCtSCMMTPy9WePcdO9z3C0Y4D/d9t6Prd1DU3dQ7zzP7fxNw/Vex2emQPWgjDGTNuO49187ucHuO7CRfzD29axKDcDgLdvruYzD9fzzeeP8+b15VyypMjjSM1sWAvCGDMtoZDyuUcOUpaXzhffufFccgDISPXzf968mor8DP7vz/YTCIY8jNTMliUIY8y0PLy3jT3NPfzFmy4kK+21nRDZ6Sl8+pY1HDrZzzefPx7/AM2csQRhjHFseCzIP/7iEGsr83jrxspJj3vTmjKuXVnKF359hJO9w3GM0MwlSxDGGMe+/uwx2nqH+dTNq/H5ZNLjRITP3LqWQEj5x18eimOEZi5ZgjDGONI/PMaXf9PADavLeN0FxVMeX1OcxR1banhkbzs9Q6NxiNDMNUsQxhhHfrn/JIOjQe6+5gLH59y2qYrRYIiH97a7GJlxi6sJQkRuFJHDItIgIp+I8fzvi8jeyO15Edkw7rnjIrJPRHaLiG0TZ4zHHtzdRm1xFhurCxyfs6YijwsX5/LjnS3uBWZc41qCEBE/8CXgJmA1cIeIrJ5w2DHgalVdD3wOuH/C89eq6kWqutmtOI0xUzvVN8xzL3exdUMFIpPXHiYSEd52cRV7mnto6BhwMULjBjdbEFuABlVtVNVR4AfA1vEHqOrzqnomcvcFoMrFeIwxM/TwnjZUYet5Ri5NZuvGCvw+4X9eslbEfONmgqgEmsfdb4k8Npk7gV+Mu6/Ar0Rkp4jcNdlJInKXiOwQkR2dnZ2zCtgYE9vPdreyrjKfC0pzpn3uotwMrl5Ryk9eaiEYssX85hM3E0SsdmjMvw4RuZZwgvj4uIevUNWLCXdRfVBErop1rqrer6qbVXVzaWnpbGM2xkzQ0DHA/tY+tl5UMeNrvO3iKk71jfBcQ9ccRmbc5maCaAGqx92vAtomHiQi64GvAVtV9XT0cVVti3ztAH5KuMvKGBNnD+5uxSdw64aZJ4jrVi0iLyPFupnmGTcTxHZguYjUiUgacDvw0PgDRKQG+AnwB6p6ZNzj2SKSG/0euAHY72KsxpgYVJUHd7dx+QUlLMrLmPqESWSk+rl5XTmPH+xgzNZnmjdcSxCqGgDuAR4DDgI/VNV6EblbRO6OHPZpoBj48oThrGXAsyKyB3gReERVf+lWrMaY2Orb+mjqHppV6yHq6hWlDIwE2NXUM/vATFy4uty3qj4KPDrhsfvGff8+4H0xzmsENkx83BgTX88cDdcMrlk5+/re5ctK8PuEp490sqXOlgGfD2wmtTFmUs82dLKyLHdW3UtR+ZmpbKwu4OmjNtpwvrAEYYyJaXgsyPbjZ7hy+dxtH3rVilL2tfZyemBkzq5p3GMJwhgT04vHuhkNhHj9HCcIVXjWhrvOC5YgjDExPXO0kzS/j0vrpl651al1lfkUZKXy1BHrZpoPLEEYY2J65mgXm2oLyUzzz9k1/T7h9ctLeeZoF6o2qzrRWYIwxrxGR/8wh072z2n9Ieqq5SV09o9wsL1/zq9t5pYlCGPMa0SXxLhq+dwvX3PVivA1bTRT4rMEYYx5jWeOdlGYlcqairw5v3ZZXgYXLs7lqcOWIBKdJQhjzKuoKs8e7eLyZSXn3Xd6Nq5cVsLOpjOMBIKuXN/MDUsQxphXaegYoKN/hNcvm/v6Q9TmJUWMBkLsb+117TXM7FmCMMa8yvbj4T28Ll06d8NbJ9q8pBCAHcfPTHGk8ZIlCGPMq+w8cYbi7DSWFGe59holOenUlWSfS0YmMVmCMMa8yktNZ7i4tnBae0/PxKbaQl5qOmPzIRKYJQhjzDmnB0Y41jXIptpC119rc20h3YOjNHYNuv5aZmYsQRhjznkpsldDXBLEkvCS3zutmylhWYIwxpyz40Q3qX5hXWW+6691QWk2hVmpbD/e7fprmZmxBGGMOeelE2dYW5lPRurcrb80GRFhU20hO09YCyJRWYIwxgAwGgixp6WXTTXudy9FbV5SRGPXoO0PkaAsQRhjAKhv62U0EIpL/SFqc+S1dlgrIiFZgjDGAJzr6olnglhbmU+a32fdTAnKEoQxBggniOqizDnZf9qpjFQ/66ry2WGF6oRkCcIYg6qy48SZuNYfojbXFrKvtdcW7ktAliCMMbScOUtn/0hcu5eiNlQXMBZUDtkGQgnHEoQxhl3NPQBs9KAFsb4qPOdib0tP3F/bnJ8lCGMM+1p6SEvxsXJxbtxfu7Igk+LsNPa02NLfiWZGCUJE7prrQIwx3tnT0svq8jxS/fH/zCgirK/KZ58liIQz078Gd5d5NMbETTCk1Lf2sqHK/eU1JrO+qoCjHf0MjQY8i8G81owShKp+da4DMcZ441jXAIOjQdZVFXgWw/qqfEIK+1v7PIvBvNaUCUJEikXkiyLykojsFJF7RcTRVlMicqOIHBaRBhH5RIznf19E9kZuz4vIBqfnGmPmxp7mcNfOeo9bEGCF6kTjpAXxA6ADeBtwG9AJ/PdUJ4mIH/gScBOwGrhDRFZPOOwYcLWqrgc+B9w/jXONMXNgX2svWWl+LijN8SyG0tx0KvIzrFCdYJwkiCJV/ZyqHovc/hYocHDeFqBBVRtVdZRwotk6/gBVfV5Vo3PsXwCqnJ5rjJkbe1t6WFuRj9/nbWlxfVWBtSASjJME8RsRuV1EfJHbO4BHHJxXCTSPu98SeWwydwK/mO65InKXiOwQkR2dnZ0OwjLGRI0FQ9S39bHOw+6lqPXV+Zw4PUTv0JjXoZiISROEiPSLSB/wx8D3gNHI7QfAxxxcO9bHkZibz4rItYQTxMene66q3q+qm1V1c2lpqYOwjDFRR08NMBIIeVp/iNoQrUO09ngah3nFpAlCVXNVNS/y1aeqKZGbT1XzHFy7Baged78KaJt4kIisB74GbFXV09M51xgzO/sib8brPRzBFLW2Mjqj2uoQiSLFyUEicitwVeTuk6r6cwenbQeWi0gd0ArcDrxzwnVrgJ8Af6CqR6ZzrjFm9va09JKbkUJtUZbXoZCfmUpdSTZ7Ist+GO9NmSBE5PPAJcB3Iw99RESuVNXzDj1V1YCI3AM8BviBB1S1XkTujjx/H/BpoBj4sogABCLdRTHPndn/ojFmMvtaellflY/P4wJ11PqqfLY12tLficJJC+Jm4CJVDQGIyLeAXcCUcxNU9VHg0QmP3Tfu+/cB73N6rjFm7owEghw62cedVy71OpRz1lcV8ODuNjr6h1mUG799KUxsTmdSF4z73vtqljFm1g619zMW1IQoUEetrQiXN+vbbEZ1InCSIP4e2CUi34y0HnZGHjPGzGP7WsPF4HWViZMgVkcTRKsVqhPBebuYRMQHhIDLCNchBPi4qp6MQ2zGGBfVt/WSn5lKVWGm16Gck5sRLlTbmkyJ4bwJQlVDInKPqv4QeChOMRlj4qC+rY81FXlEBogkjNUVeTajOkE46WL6tYj8uYhUi0hR9OZ6ZMYY14wFQxxq7z839yCRrK3Ip7n7rM2oTgBORjH9UeTrB8c9pkDiDH0wxkxLQ8cAo8EQayqczHmNr7WV0UJ1L5cvK/E4muQ2ZYJQ1bp4BGKMiZ/9kSLwmorEa0FEY9pvCcJzTibKZQAfAK4k3HJ4BrhPVYddjs0Y45L6tj4yU/3UlWR7HcprFGWnUZGfYYXqBOCki+nbQD/wxcj9O4DvAG93KyhjjLvq23pZXZHn+RLfk1lTmU99mw119ZqTBLFSVTeMu/8bEdnjVkDGGHeFQsqBtj5u21Q19cEeWVuRz/8ePMXgSIDsdEdLxhkXOBnFtEtELoveEZFLgefcC8kY46bjpwcZHA0mZP0ham1lHqpwsN26mbzkJEFcCjwvIsdF5DjwW+BqEdknIntdjc4YM+eiy1isqUy8EUxR5wrVNqPaU07abje6HoUxJm72t/WS6heWL8r1OpRJleWlU5KTZmsyeczJMNcT8QjEGBMfB9r6WLk4l7QUp2t1xp+IsKYin/2WIDyVuH8hxpg5p6rsb+1lTXni1h+i1lbmcfRUP8NjQa9DSVqWIIxJIu29w5wZGjs3WzmRranIJxBSjpzq9zqUpOUoQYhIrYi8MfJ9pogkbuelMWZS0aLv6gQewRQVXQbkgHUzeWbKBCEi7wd+DHw18lAV8DMXYzLGuKS+rQ8RWFWe+J/xqguzyE1PsUK1h5y0ID4IXAH0AajqUWCRm0EZY9xR39bH0pJsstISf/KZzyesKs+zGdUecpIgRlR1NHpHRFIIr8lkjJlnDrT1JvQEuYlWV+RxsL2fYMjecrzgJEE8JSJ/BWSKyPXAj4CH3Q3LGDPXzgyO0tY7nJBLfE9mTUUeZ8eCHOsa9DqUpOQkQXwc6AT2AX8MPAr8HzeDMsbMvQORZSvmUwsiGusBW3LDE072pN6rqmuB/4xPSMYYN0T78lfPoxbEskU5pPl91Lf1cuuGCq/DSTrnbUGoagjYIyI1cYrHzBOq1ic839S39VGen0FRdprXoTiWluJjeVmODXX1iJMupnKgXkQeF5GHoje3AzOJ6/mGLjb97f/yWP1Jr0Mx01Df1jev6g9RayryqG/rsw8lHnAy1u0zrkdh5o1nj3Zx57e2MxII8V8vnOBNaxZ7HZJx4OxokMbOAW5eV+51KNO2piKfH+5o4WTfMOX5mV6Hk1ScLNb3VDwCMYnv6SOdvP/bO6gryWZjTSE/3NFM18AIJTnpXodmpnDoZB8hZd62ICA8o9oSRHw5mUndLyJ9kduwiARFxDoEk0x771ne/+0dLC3N4Xvvv4z3XF5LMKT8Yl+716EZB87tATEPE8SF5XmIYDOqPTBlglDVXFXNi9wygLcB/+Hk4iJyo4gcFpEGEflEjOcvFJHfisiIiPz5hOeORzYl2i0iO5z+Dxl3vNB4mpFAiH9++3qKstNYWZbL8kU5PLzHEsR8UN/WR35mKpUF8+8TeE56CkuKs21GtQemvZqrqv4MeMNUx4mIH/gScBOwGrhDRFZPOKwb+DDwz5Nc5lpVvUhVN083TjO3djX1kJXm58LF4U+gIsItGyp48Xg3bT1nPY7OTOVAWy+ry/MQEa9DmZHVkUK1iS8nXUxvHXe7TUQ+j7OlNrYADaraGFmq4wfA1vEHqGqHqm4HxmYSvImf3c09rK/Kx+975Q3mlsi49Ef2WisikQWCIQ6d7J+X3UtRayryaDlzlt4he6uIJyctiFvG3d4E9DPhjX4SlUDzuPstkcecUuBXIrJTRO6a7CARuUtEdojIjs7Ozmlc3jg1PBbkQFsfG2sKX/V4XUk26yrzeXhvm0eRGSde7hxkJBBK6D2opxKdUV3fbt1M8eQkQXxNVd8bub1fVf8OWO7gvFht2ekMZL5CVS8m3EX1QRG5KtZBqnq/qm5W1c2lpaXTuLxxan9rL4GQsrG64DXP3bKhnL0tvbZWTgI7EHlTnU9LbExke0N4w0mC+KLDxyZqAarH3a8CHH/UVNW2yNcO4KeEu6yMB3Y19QBwUU3Ba557y/pwN9OvbNJcwqpv7SM9xcfSkmyvQ5mxkpx0yvMzzm14ZOJj0nkQIvI64HKgVET+dNxTeYDfwbW3A8tFpA5oBW4H3ukkKBHJBnyq2h/5/gbgs07ONXNvd3MPVYWZLMrNeM1zFQWZlOdncOikbQuZqPa39bKqPI8U//zeYXhNRR77rQURV+ebKJcG5ESOGb/9VB9w21QXVtWAiNwDPEY4oTygqvUicnfk+ftEZDGwg3DSCYnIRwmPeCoBfhoZcZECfE9VfznN/zczR3Y1nWHTkqJJn19elmv7BieoUEipb+1j68b5v9Ddmop8njjUwdBoYF5seLQQTPpTjsygfkpE/ltVD41/TkRKnFxcVR8lvDz4+MfuG/f9ScJdTxP1ARucvIZx18neYdp6h7kzRv0hasWiHL7TeJpgSF81ysl4r/nMEP0jAdbO4/pD1NrKfEIKB9v72VRbOPUJZtactDl/KCKXRe+IyNuA590LySSS3c1nANgYo/4QtaIsl5FAiObuoThFZZza3xrukllbOf8TRLRQbRPm4sdJO+33gQdE5EmgAijGwUQ5szDsauohze877xj65WU5ABw51c+SeVwIXYj2t/WS4pNzv6P5LLpUeX2r1SHixclSG/uAvwPuBq4F7lHVFrcDM4lhV3MPqyvySE+ZfFzC8rJwiepox0C8wjIO7W/tZUVZ7nl/f/OFiEQK1daCiBcnM6m/DnwUWA+8F3hYRD7oclwmAQSCIfa29Jy3ewnCa+VUFmRy2EYyJRRVpb6tj7XzeILcRGsq8jlyqp/RQMjrUJKCkxrEfsJrIh1T1ceAy4CL3Q3LJILDp/oZHgtx0XkK1FHLy3JsJFOCae8dpntwdEHUH6LWVuYxFlT7W4sTJ11MX9BxWzmpaq+q3uluWCYRNES6jFaXT/0JdEVZLo2dgwSC9skuUUQnlc3nGdQTnVtyw7qZ4mLSBCEiP4x83Scie8fd9onI3viFaLxyrGsQEaguypry2OWLchgNhjhhI5kSRn1bHz6BVeW5Ux88T9QWZZGTnnJudJZx1/lGMX0k8vUt8QjEJJ5jXYNU5GeSkTp1gXNFtFB9qp8LSuf/iJmFoL6tlwtKcxbUpDKfT1htheq4mbQFoartka8nVPUEcIbwSq7Rm1ngjncNUudw2OorQ11tJFOi2N/at6DqD1FrKvI42N5HMDSdtT/NTDgZxfTHInIK2AvsjNxsh7cFTlU51jXIkpKpu5cAstJSqC7KtOJhgujsH+Fk3/C83gNiMmsr8hkeC/Fyp30YcZuTtuefA2tUtcvtYEziODM0Rt9wgCXFzie+rViUy1FrQSSEaBF3IRWoo6KtougcD+MeJ8NcXwas8phkovs7LC11niCWl+XS2DXAmI1k8lx0e87VC7AFsWxRDpmpfva2WB3CbU5aEJ8EnheRbcBI9EFV/bBrURnPHY8kiGm1IMpyGAsqx7sGz82uNt7Y29JDbXEW+ZmpXocy5/w+YW1lHvtsbwjXOWlBfBV4AniBV2oQO90MynjvWNcgfp84GuIaFW3uW6Hae/taellfVeB1GK5ZV1lAfVuvzbtxmZMWREBV/3Tqw8xCcuz0IFWFmaROY5OZZYtyEAkv2vdmyl2MzpxPZ/8Ibb3D/FHVwqs/RK2vyueB50Ic7RhglYOJnGZmnPzr/42I3CUi5SJSFL25Hpnx1PGuwWl1LwFkpPqpyM/kxGnbn9pLe1t6ABZ2CyKS/PZZHcJVThLEO4nUIbBhrklBVac1B2K82uIsm03tsb0tvfiEBTnENaquOJvc9BT2tvZ4HcqCNmUXk6rWxSMQkzg6B0YYHA3OKEHUFGXx6wOnXIjKOLW3pYdli3LITl84M6gn8vmEtZX51oJw2fzexdy44nhXuAUwk81/aoqzOD04ysBIYK7DMg6oKnsXeIE6an1VPgfbbelvN1mCMK9xrCs8CqlumjUIgNqi8DlNp62byQttvcOcHhxlwwIuUEetq8pnNBiyfUhcZAnCvMaxriFS/UJFQca0z62JDIttsjqEJ/Y29wCwLhlaEJUFAFaHcJGjTkoRqQRqxx+vqk+7FZTx1vGuQaqLskiZxhDXqJriaIKwkUxe2NvaS6pfFtQS35OpLsqkICs1XIe41OtoFqYpE4SI/CPwe8ABIBh5WAFLEAvU8dODM+peAsjPTCU/M5UT1sXkib0tPVy4+Px7iC8UIsK6ynxbcsNFTloQvwOsVNWRqQ40818opBw/PciVy0pmfI3a4izrYvJAKBQuUN+yocLrUOJmfVU+X32qkeGxoKN9S8z0OOlDaAQW3oIuJqZT/cMMj4VmNIIpqqbIEoQXTnQP0T8cSIoCddS6ygICIeVgu+0w5wYnLYghYLeIPI4t1rfgHesM1w5mMgciqqYoi1/uP0kgGJpRHcPMTDLMoJ5oQ3U4Ge5u7mFjTaHH0Sw8ThLEQ5GbSQLRWdC1xc4X6ZuotjiLQEhp7x2e1mJ/Znb2NPeSkepj+aLk2fK1PD+TxXkZ7Grq4b1XeB3NwuNkJvW34hGISQxN3eEhruX5mTO+Rk1kLsSJ00OWIOJod/MZ1lbkJ12r7eLaAl5qOuN1GAuSky1Hl4vIj0XkgIg0Rm9OLi4iN4rIYRFpEJFPxHj+QhH5rYiMiMifT+dc446m7iGqCrPw+2TG13hlqKvVIeJlJBBkf2sfm2qTr5tlY3UhLWfO0tE/7HUoC46TjxrfAL4CBIBrgW8D35nqJBHxA18CbgJWA3eIyOoJh3UDHwb+eQbnGhc0d8/+U//ivAzS/D5O2FyIuNnf2sdoMJSU/fAX1xYAsKupx9M4FiInCSJTVR8HRFVPqOrfAG9wcN4WoEFVG1V1FPgBsHX8AaraoarbgbHpnmvcceL0EDVFM+9egvCOX1WFmTRbCyJudkW6WKJvlslkTUU+qX6xBOECJ0XqYRHxAUdF5B6gFVjk4LxKoHnc/Racz3eczblmhnqHxug9O3ZuuYzZqCnOsslycfRS0xmqCjNZlDv95VHmu4xUP6sr8q0O4QInLYiPAlmEu4I2AX8AvMfBebE6sdVhXI7PjWxmtENEdnR2djq8vIml+Uz4DT1aZJ6N2qIsmk4Poer0V25m46UTPUlZf4jaWF3A3pYe24J0jk2ZIFR1u6oOqGqLqr5XVd+qqi84uHYLUD3ufhXQ5jAux+eq6v2qullVN5eWljq8vIklWlSeixZEdVEW/SMBeoYm9h6audbWc5aTfcNcnIT1h6iLawsZHgtxyFZ2nVOTdjGJyL+p6kdF5GFifHpX1VunuPZ2YLmI1BHulrqd8O50TszmXDND0QRRPcsaBEBtZC2nE91DFGanzfp6ZnI7T0TqD0mcIDZWFwDhrra1lckzk9xt56tBREcq/fN5jpmUqgYiNYvHAD/wgKrWi8jdkefvE5HFhLcvzQNCIvJRYLWq9sU6dyZxGOdOnB6iKDuN3IzZr6xSO26o60WRf7zGHS81nSEj1ceFSbCC62SqCjMpzU1nV1MP736d19EsHJMmCFXdGfn6VPQxESkEqlV1r5OLq+qjwKMTHrtv3PcnCXcfOTrXuGsuhrhGVRdGEsRpG+rqtpeaelhfVUBqkk2QG09EuLjGJszNNScT5Z4UkTwRKQL2AN8QkX91PzQTb03dQ9TOUYLITPOzKDfdRjK5bHgsyIG23qQuUEdtrCnkxOkhTg/YwtNzxclHjnxV7QPeCnxDVTcBb3Q3LBNvgWCI1p6zc1Kgjqotzjq3tpNxx/7WXsaCmtT1h6joz8DmQ8wdJwkiRUTKgXcAP3c5HuOR9t5hgiGd4wSRzQnrYnJVtEC9sabA20ASwLrK8IS57Se6vQ5lwXCSID5LuFjcoKrbRWQpcNTdsEy8RbuC5nJxvdqiLE71jXB2NDj1wWZGXmo6Q21xFiU56V6H4rnMND/rqwrY1mgJYq44mQfxI1Vdr6ofiNxvVNW3uR+aiadzcyBmscz3RLWRPSVs0T53hELKi8e62bKkyOtQEsaldUXsa+1lcCTgdSgLgpMi9T9FitSpIvK4iHSJyLviEZyJn6buIdL8Phbnzd1SDdGCt3UzueNIRz9nhsa4dGmx16EkjEuXFhMM6bmuNzM7TrqYbogUqd9CeIbzCuAvXI3KxF1z9xBVhZmzWuZ7oiXFr+wLYeZetCvl0jprQURtqi3E7xNePGbdTHPBSYKIzpq6Gfi+qtpPfgE60T0455v75Gelkp+Zast+u+SFxtNUFWbapkzj5KSnsLYyn23HTnsdyoLgJEE8LCKHgM3A4yJSCtjOHAtM0+mhOR3BFLXEVnV1RSikbDvWzWXWvfQal9UVsae5l+ExGxwxW06K1J8AXgdsVtUxYBDbm2FB6R0ao2844EqCqC3O5rjVIObc0Y4BugdHLUHEcOnSIkaDIZtVPQem3A9CRN497vvxT33bjYBM/LkxgimqtjiLn+9tYzQQIi0leZeCmGsvNIa7UKz+8FqblxThk3CN5vILSrwOZ15zsmHQJeO+zwCuA17CEsSCEa0RuNWCCCm09pylrmT2+0yYsBcaT1NZYPWHWPIyUlldkWd1iDkwZYJQ1Q+Nvy8i+TjYk9rMH9EagTsJ4pWhrpYg5ka0/nDtSicbOyanS+uK+a8XTjASCJKe4vc6nHlrJm3+IWD5XAdivPNy5wCL8zLITnfSoJyeVxKEFarnyiv1B+temsyldUWMBELsae71OpR5zUkNYvyGQT5gNfBDN4My8dXYOcjSUnc+3ZfmpJOV5rcEMYeiXSdWoJ7clroiRMJdcVusTjNjTj4yjt8wKACcUNUWl+IxcaaqNHYOcMuGCleuLyLUFGXZbOo5ZPWHqRVkpbGuMp+nj3Ty4eusw2OmnNQgnprqGDN/dQ+O0jccYGlpjmuvUVucRUPHgGvXTyaBYIjnGk5z/eoyr0NJeFevKOXLT75M79kx8jNnv0tiMnKyFtNlIrJdRAZEZFREgiLSF4/gjPuOdYU/2S91sYC8pDib5u6zBEOv2drcTNOelh56z45xzcpSr0NJeFevKCUYUp5v6PI6lHnLSZH6P4A7CC/xnQm8D/iim0GZ+GnsjCQIl2oQEJ5fMRoMcbLPJuDP1pOHO/EJvH6ZJYipXFRdQG5GCk8d6fQ6lHnL0SgmVW0A/KoaVNVvANe6G5aJl8auQVL9QlWhe/3ZryzaZ3WI2frN4Q4urikkP8u6TKaS4vdx5bISnjrSiaq1XmfCSYIYEpE0YHdk6e+PATagfYFo7Bygtjh7TldxnaimyIa6zoWO/mH2t/ZZ99I0XL2ilPbeYauBzZCTBPEHkePuIbwOUzVgGwYtEI1dg67WHwAqCjJJ9YsliFl6+ki4L/0amyDn2FUrwsnUuplmxslifSdUdVhV+1T1M6r6p5EuJzPPBUManuHsYv0BwO8TqguzON5lXUyz8eThDkpy0lldnud1KPNGRUEmK8pyLEHM0KQJQkS2isgHx93fJiKNkdtt8QnPuKnlzBBjQeWCEveGuEZdsCiHhk5r5s9UIBjimaNdXL2iFJ+L3YEL0VXLS9nW2M3QqG1DOl3na0H8JfDQuPvphBfuuwb4ExdjMnHS2OX+CKaolWW5HOsaZCRga/TPRHR467UXWv1huq5eWcpoMHRuBz7j3PkSRJqqNo+7/6yqnlbVJqxIvSBEh7jGYxG95WU5BEN6bt6FmR4b3jpzlywpIiPVx28Od3gdyrxzvgRROP6Oqt4z7q79lS4Ax7oGyM9MpSg7zfXXWrk4F4DDJ/tdf62FyIa3zlxGqp9rVizil/tPErLJmtNyvgSxTUTeP/FBEflj4EX3QjLx0tgZXoJ7wkZQrqgrCQ+lPXLKEsR0NXcPsb+1jzfa8hozdtO6xXT0j9guc9N0vrWYPgb8TETeSXiDIIBNhGsRv+NyXCYOjnUN8roL4rMiaHqKn7qSbI6cskL1dD2yrx2AN68r9ziS+esNFy4iLcXHo/tOsnmJre7q1KQtCFXtUNXLgc8BxyO3z6rq61T1lJOLi8iNInJYRBpE5BMxnhcR+ffI83tF5OJxzx0XkX0isltEdkz3f8yc39BogPbeYdfnQIy3sizXWhAz8MjedjZU5dvqrbOQm5HKVctL+cX+dutmmgYn8yCeUNUvRm5POL2wiPiBLwE3Ed5D4g4RWT3hsJsIbz60HLgL+MqE569V1YtUdbPT1zXOvLIGk/tDXKOWl+XQ1D3E2VEbyeRU0+kh9rX28ub11nqYrZvXLaa9d5g9LT1ehzJvuLmL/BagQVUbVXUU+AGwdcIxW4Fva9gLQIGI2L+EODgWxyGuUSvLclHFlj2Yhmj30s3WvTRr160qI9Uv/GL/Sa9DmTfcTBCVwPhhsi2Rx5weo8CvRGSniNw12YuIyF0iskNEdnR22mxJp452DOCTVxbSi4flZeGRTNbN5Nwj+9q4qLrA1cUUk0V+ZipXLivh0X3ttnifQ24miFhDYyb+Vs53zBWqejHhbqgPishVsV5EVe9X1c2qurm01EbfOrW/tZdli3LISI3fhu5LirNI8/ssQTh0vGuQ/a19VpyeQzetK6flzFn2t9qWNk64mSBaCC/sF1UFtDk9RlWjXzuAnxLusjJzZF9rL2sr8+P6mil+HxcsyrEE4VC0e+mmdYs9jmThuGF1GSk+OfezNefnZoLYDiwXkbrIcuG38+qlO4jcf3dkNNNlQK+qtotItojkAohINnADsN/FWJPKqb5hOvtHWBfnBAGwoizHhro69MjedjbWWPfSXCrISuOqFaX8dFcLgWDI63ASnmsJQlUDhJcIfww4CPxQVetF5G4RuTty2KNAI9AA/CfwgcjjZcCzIrKH8KS8R1T1l27Fmmz2tfQCeJQgcmntOUv/8FjcX3s+2d/ay4H2PrZuqPA6lAXn9y6p5lTfCE8etprlVM43UW7WVPVRwklg/GP3jftegQ/GOK8R2OBmbMlsX2svPoHVFfFfNnpFpFB9tGOAi2sKpzg6eX3vxSbSU3z87sYqr0NZcN5w4SJKc9P5wfYmm50+BTe7mEyC2t/aywWlOWSlufr5IKaV0ZFMtibTpAZGAjy4q5W3rK+wtZdckOr3cdumKp441MHJXtsn/XwsQSShfa29nnQvAVQVZpKZ6rc6xHk8uLuVwdEg77y0xutQFqzf21xNSOHHO5unPjiJWYJIMqf6hunoH4n7CKYon09YXpbD4VM2zDAWVeV725q4cHEuF9cUeB3OgrWkJJvXLS3mv3c029Ib52EJIsmcK1BXeZMgADZUFbC7qcdGkcSwt6WX+rY+fv/SmrisspvMbt9STXP3WZ5/+bTXoSQsSxBJZl9rLyJ4uq/xlroiBkeDHGi3VsRE39vWRGaqn60bJy46YObam9YsJj8zle9uO+F1KAnLEkSSiRaos9PjX6CO2lIXXm75xWO2BeR4ZwZHeWhPG7duqCAvw4rTbstI9fPOS2v4Zf1JWx9sEpYgkoyXBeqosrwMlhRnsc0SxKt87dlGhgNB7nx9ndehJI07r6wjPcXHV5582etQEpIliCTSESlQe50gINyK2H682wqEET1Do3zr+RPcvLb83FwR476SnHTu2FLDz3a30tw95HU4CccSRBLZ1+p9gTpqS10xPUNjHLWmPQBff/YYAyMBPnTdMq9DSTp3XbUUn8B9T1krYiJLEElkT4v3BeqoS8/VIWwESe/QGN987jg3rV3MhYu9/90km/L8TG7bVM2PdrRwqs8mzo1nCSKJPHWkk/WV+Z4WqKOqCjMpz8+wOgTwwHPH6B8J8KE3LPc6lKT1J1dfQFCV+59u9DqUhGIJIkl09A2zp7mH6xNk7RkRYUtdES8e607qzVu6B0d54Llj3LC6zJO1sUxYTXEWv7uxku/89sS53RaNJYik8fihDoCEWpxsS10RHf0jnDidvMXBf3j0IGdHg/zZDSu9DiXp/eWNK0lP8fHpB/cn9YeW8SxBJInHD56iqjDz3GJ5ieDSJJ8Psa3xND/a2cL7Xr+UlYsT5/eSrBblZvCx61fwzNEufmn7VgOWIJLC2dEgzxzt4o2ryhJq+YYLSnMoyk5LyjrEaCDEp362n6rCTD5yndUeEsW7X1fLqvI8PvvzAwyNBrwOx3OWIJLAsw1djARCvHFV4nQvQbgOcdnSIp452pl06zLd//TLNHQM8LnfWUtmWvz2BTfnl+L38bmta2jvHebex496HY7nLEEkgf89cIrc9JRzS1wkkls3VNLRP8JTR5Jnd6/DJ/v59ycaePO6cq5ducjrcMwEm5cU8Xubq/nPpxt5rqHL63A8ZQligQuFlMcPneLqlaWkpSTer/u6VYsoyUnn+y8mx7r8vUNj3PWdHRRkpvLXt672OhwziU/fspoLSnP48Pd30d571utwPJN47xhmTu1u6aFrYDRhhrdOlOr38fbNVfzmcMeCn6QUDCkf/sEu2nrO8pV3bWJRbobXIZlJZKen8JV3bWJ4LMgHvvsSo4Hk6gKNsgSxwD1WfxK/T7hmReJ2Zdx+STXBkPKjHQu7FfEvvzrMU0c6+cyta9lUa/txJ7pli3L4p9s2sKuph8/9/EBSDn21BLGAnRkc5XsvNPHGVYsSem/j2uJsLr9gYe/u9e3fHufLT77MHVtqbCvReeTN68t5/+vr+M4LJ/inxw4nXZKwBLGAfeWplxkYDfCn1yf+JKzbt9TQ3H2W515eeEXB+556mU8/WM8bV5XxN1Z3mHc+edMq3nlpDV958mU+/4tDSZUkvF+Ux7jiZO8w33r+OL97UeW8mIT1pjVlFGal8oMXm3n98lKvw5kTqsoXfn2Ef3+igVs2VPCv79hAqt8+k803Pp/wd7+zFr8IX326kbGg8qk3r8LvS5w5RW6xv9YF6t7HjxJS5WPXr/A6FEfSU/y845JqHt3fzvMLYGhh79AY93xvF//+RAPv2FzFv/3eRZYc5jER4bNb1/CHly/hgeeO8a6vbaNjgQ+qAEsQC9KxrkF+uKOZO7bUUF2U5XU4jn3kuuXUlWTz0f/eTffgqNfhzNhvXz7Njfc+zWP1J/nLG1fy+beuT4pPmwudiPDXt6zmn25bz+7mHm669xmePNzhdViusgSxwAyPBfnUT/eR5vdxzxvm1+YzWWkpfPGOjfQMjfEXP9oz7/p6m7uH+Isf7eGdX3uBzFQ/P/3AFXzgmmX4LDksGCLCOzZX8/CHrqAkJ50//MZ23v/tHTR09HsdmissQSwgQ6MB/uib2/lt42k+u3XNvBxnv6Yin0/efCGPH+rgm88f9zocR5pOD/HpB/fzhn95kgf3tHHnFXX8/MNXJsTOfcYdyxbl8uA9V/Bn16/gty+f5oYvPM3Hf7yXg+19Xoc2p2S+fUo7n82bN+uOHTu8DsMT/cNjvPcb23mp6Qz/8o4N/O7GKq9DmjFV5f3f3sHjhzq484o6/uyGlQm3XtGZwVF+deAk/7OzlRePd5PiE96+uZoPX7eM8vxMr8MzcXR6YIQvPtHA97Y1MRoMsaEqn7dvruaNq8pYnJ/4H9JEZKeqbo75nJsJQkRuBO4F/MDXVPXzE56XyPM3A0PAH6rqS07OjSUZE8TASIAf7WjmgeeO0d4zzL23b+TN68u9DmvWBkcC/P2jB/nutiaWFGfx929dx+uWFnuyGm0opJzoHuLwyT52njjDcw2nOXiyD1VYWprN2y6u4q0XV1piSHJnBkf56a5W/nt7M4dPhbucVpTlcOWyUi6qKWB1eR51JdkJV4/yJEGIiB84AlwPtADbgTtU9cC4Y24GPkQ4QVwK3Kuqlzo5N5aFnCBUlZFAiJ6hMY51DXL4ZB8H2vv4xb6T9I8E2FRbyJ9ev4IrlpV4Heqcer6hi4//ZC/N3Wcpy0vnquWlXLGshOqiTMryMliUmzGtNaZUlUAo/LMcGQsyNBpkYCTAwEiAnqExzgyO0j00ysneYVp7ztJy5izHuwY5OxYEIM3vY2NNAVcuK+GqFaWsr8pPqCXUjfdUlUMn+3nmaCfPHO1i27Huc0t1ZKT6qC3Kpqowk+qiLEpz0ynOTqM4J538zFRy0lPIzUghI9VPRqqPjFS/66PfvEoQrwP+RlXfFLn/SQBV/Ydxx3wVeFJVvx+5fxi4Blgy1bmxzDRB3PLFZxmOvAHEw/if+Pifv0b+E1JFgUBQGQuGGAuGGBwJMjphSezCrFRev7yUP7qyjouqC+IQuTeGRgP8fG87Tx3p5NmjXfSeHXvV8yk+ISPVT1qKj1c+nAmqSkiVkEIgGGIsFP55OvmTz0lPobIgk8rCTGqLs1i1OI8Ly3NZUZZLRmpidXeZxDYaCPFy5wAH2vo42N7Hie4hmruHaD1zlv6RqfecEAmvWZbqE1L8PlJ8gs8n+EXwSbhwXpKTxoP3XDmj+M6XINycKFcJjF9cp4VwK2GqYyodnguAiNwF3AVQUzOzJQwuKM1+zZuv2wQZf+ccn4SfEYEUn49Uv5DiF3LSU8nLTCEvI5XqoixWLc6lNDc9KT69ZqWl8I7N1bxjczWBYIiGzgHae4c51TtMZ/8IZ8eC4RZBIIhqONGqgk/A7wv/PFP8PlL8QqrPR3pK+JNZeqqPzFQ/uRkpZKenkJ+ZSlF2GkXZaWSm+pPiZ2vcl5biY1V5HqvKX7vn+PBYkNODo5weGKF/OBC5jTE8FmR4LMTwWPiD4VhQGQ2ECIZCBFUJhsK3UOQDZU66O2/lbiaIWP+6Jn52m+wYJ+eGH1S9H7gfwi2I6QQY9W+3b5zJacYDKX4fFy7O48LFr/3HZsx8k5HqD7dUCxKzfuVmgmgBqsfdrwLaHB6T5uBcY4wxLnKz+rEdWC4idSKSBtwOPDThmIeAd0vYZUCvqrY7PNcYY4yLXGtBqGpARO4BHiM8VPUBVa0Xkbsjz98HPEp4BFMD4WGu7z3fuW7Faowx5rVsopwxxiSx841isqU2jDHGxGQJwhhjTEyWIIwxxsRkCcIYY0xMC6pILSKdwIkZnl4CJOJWZhbX9Fhc02NxTc9CjKtWVWPu87ugEsRsiMiOySr5XrK4psfimh6La3qSLS7rYjLGGBOTJQhjjDExWYJ4xf1eBzAJi2t6LK7psbimJ6nishqEMcaYmKwFYYwxJiZLEMYYY2KyBDGBiHxIRA6LSL2I/JPX8YwnIn8uIioiCbHxtIj8PxE5JCJ7ReSnIlLgYSw3Rn5vDSLyCa/imEhEqkXkNyJyMPI39RGvY4oSEb+I7BKRn3sdy3giUiAiP478bR2MbF/sORH5WOR3uF9Evi8iGR7F8YCIdIjI/nGPFYnIr0XkaORr4Vy8liWIcUTkWmArsF5V1wD/7HFI54hINXA90OR1LOP8GlirquuBI8AnvQhCRPzAl4CbgNXAHSKy2otYYggAf6aqq4DLgA8mUGwfAQ56HUQM9wK/VNULgQ0kQIwiUgl8GNisqmsJb0Nwu0fhfBO4ccJjnwAeV9XlwOOR+7NmCeLV/gT4vKqOAKhqh8fxjPcF4C+ZZOtVL6jqr1Q1uuv6C4R3/vPCFqBBVRtVdRT4AeFE7zlVbVfVlyLf9xN+s6v0NioQkSrgzcDXvI5lPBHJA64Cvg6gqqOq2uNpUK9IATJFJAXIwqNdLlX1aaB7wsNbgW9Fvv8W8Dtz8VqWIF5tBfB6EdkmIk+JyCVeBwQgIrcCraq6x+tYzuOPgF949NqVQPO4+y0kwJvwRCKyBNgIbPM4FIB/I/yBI+RxHBMtBTqBb0S6v74mItleB6WqrYR7FJqAdsK7X/7K26hepSyyGyeRr4vm4qJu7kmdkETkf4HFMZ76FOGfRyHhroBLgB+KyFKNw1jgKeL6K+AGt2OI5XxxqeqDkWM+Rbgr5bvxjG0cifFYwrS0AEQkB/gf4KOq2udxLG8BOlR1p4hc42UsMaQAFwMfUtVtInIv4e6S/+tlUJE+/a1AHdAD/EhE3qWq/+VlXG5LugShqm+c7DkR+RPgJ5GE8KKIhAgvgtXpVVwiso7wH+UeEYFwN85LIrJFVU96Fde4+N4DvAW4Lh6JdBItQPW4+1V41PyPRURSCSeH76rqT7yOB7gCuFVEbgYygDwR+S9VfZfHcUH4d9miqtFW1o+Zo/70WXojcExVOwFE5CfA5UCiJIhTIlKuqu0iUg7MSfe4dTG92s+ANwCIyAogDY9XblTVfaq6SFWXqOoSwv+ALo5HcpiKiNwIfBy4VVWHPAxlO7BcROpEJI1w8fAhD+M5R8JZ/evAQVX9V6/jAVDVT6pqVeTv6XbgiQRJDkT+rptFZGXkoeuAAx6GFNUEXCYiWZHf6XUkQPF8nIeA90S+fw/w4FxcNOlaEFN4AHggMnxsFHiPh5+K54P/ANKBX0daNy+o6t3xDkJVAyJyD/AY4dElD6hqfbzjmMQVwB8A+0Rkd+Sxv1LVR70LKeF9CPhuJNk3Au/1OB4i3V0/Bl4i3J26C4+W3RCR7wPXACUi0gL8NfB5wl3idxJOZm+fk9ey9z9jjDGxWBeTMcaYmCxBGGOMickShDHGmJgsQRhjjInJEoQxxpiYLEEYY4yJyRKEMcaYmCxBGOMSEbkksldGhohkR/YSWOt1XMY4ZRPljHGRiPwt4fWOMgmvMfQPHodkjGOWIIxxUWS5iO3AMHC5qgY9DskYx6yLyRh3FQE5QC7hloQx84a1IIxxkYg8RHiHuzqgXFXv8TgkYxyz1VyNcYmIvBsIqOr3IvtmPy8ib1DVJ7yOzRgnrAVhjDEmJqtBGGOMickShDHGmJgsQRhjjInJEoQxxpiYLEEYY4yJyRKEMcaYmCxBGGOMien/A+JneM1+rWmCAAAAAElFTkSuQmCC\n", "text/plain": [ "