analytical_hmc.py

from numpy import sin, cos, arctan, sqrt, exp, random, pi, linspace
import matplotlib.pyplot as plt

def draw_sample(xold, sigma):

    t = 3.0
    vold = random.normal()
    phi = arctan(-vold / xold * sigma)
    A = vold * sigma * sqrt(xold ** 2 / sigma ** 2 / vold ** 2 + 1)
    xnew = A * cos(t / sigma + phi)
    vnew = -A / sigma * sin(t / sigma + phi)

    E = lambda x: 0.5 * x ** 2 / sigma ** 2
    K = lambda v: 0.5 * v ** 2
    H = lambda x, v: E(x) + K(v)

    p_acc = min(1, exp(-(H(xnew, vnew) - H(xold, vold))))

    if random.random() < p_acc:
        return xnew, True
    else:
        return xold, False


sigma = 2.0
samples = [2.0]
accepted = 0
n_samples = 100000
for _ in range(n_samples):
    new_state, acc = draw_sample(samples[-1], sigma)
    samples.append(new_state)
    accepted += acc


fig, ax = plt.subplots()
ax.hist(samples, bins=40, density=True)
gaussian = lambda x: exp(-0.5 * x ** 2 / sigma ** 2) / sqrt(2 * pi * sigma ** 2)
xspace = linspace(-5, 5, 300)
ax.plot(xspace, list(map(gaussian, xspace)))
plt.show()

print("Acceptante rate:", accepted / n_samples)

imgs.sh

#!/bin/bash

img_list=$(ls -v output*.png)
b=$(<$2)

while read strA <&3 && read strB <&4; do
    rstring="..\/..\/img\/posts\/${strB}"
    echo $rstring
    sed -i "s/${strA}/${rstring}/g" $1
    mv $strA $strB
    # cp $strB ~/projects/tweag/www/app/assets/img/posts/
done 3<<<"$img_list" 4<<<"$b"

# cp $1 ~/projects/tweag/www/app/views/posts/

mcmc_introduction_gibbs.ipynb

mcmc_introduction_HMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to MCMC, part III: Hamiltonian Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is part 3 of a a series of blog posts about MCMC techniques:\n",
    "- [the basics](https://www.tweag.io/posts/2019-10-25-mcmc-intro1.html)\n",
    "- [Gibbs sampling](https://www.tweag.io/posts/2020-01-09-mcmc-intro2.html)\n",
    "- Hamiltonian Monte Carlo\n",
    "\n",
    "So far, we discussed two MCMC algorithms: the Metropolis-Hastings algorithm and the Gibbs sampler.\n",
    "Both algorithms can produce highly correlated samples—in the case of Metropolis-Hastings because of the pronounced random walk behaviour, while the Gibbs sampler can easily get trapped when variables are highly correlated.\n",
    "In this blog post, I'd like to introduce you to Hamiltonian Monte Carlo: [originally proposed to simulate quantum field theories on lattices](https://doi.org/10.1016%2F0370-2693%2887%2991197-X), it is a MCMC algorithm which is able to generate less correlated proposal states with high acceptance probability and thus leads to faster mixing of the Markov chain.  \n",
    "\n",
    "As always, you can [download the corresponding IPython notebook](https://github.com/tweag/blog-resources/blob/master/mcmc-intro/mcmc_introduction.ipynb) and play around with the code to your heart's desire!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonian Monte Carlo\n",
    "As mentioned above, when doing Metropolis-Hastings sampling with a naive, uninformed proposal distribution, we are essentially performing a random walk without taking into account any additional information about the distribution we want to sample from. \n",
    "We can do better!<br>\n",
    "When considering not discrete probabilities, but differentiable densities, one such additional piece of information is the local shape of the density.\n",
    "The simplest measure of that is the density's derivative.\n",
    "It essentially tells us, given a certain value $x$ for the variables, how to change $x$ in order to increase the value $p(x)$ of the density.\n",
    "That means we should somehow be able to use the derivative of $p(x)$ to find high-probability proposal states&mdash;which is exactly the key idea of Hamiltonian Monte Carlo (HMC).<br>\n",
    "As we have done in previous blog posts, we will rather consider the log-probability of the density $p(x)$.\n",
    "If you were to just turn your log-density upside-down, meaning, multiplicating it by $-1$, you can imagine the variable $x$ as the position of a particle (for example, a little marble) in a mountainous landscape.\n",
    "You place the marble somewhere and it will roll downhill towards a region with higher probability (remember: a valley means high probability, because everything is upside-down).\n",
    "Your upside-down log-density (let's call it $E(x)$) thus gives rise to a gravitation-like force on the marble.\n",
    "Here's how all that looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAEGCAYAAADGwUaDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3QUVRvH8e+kNwiQhBp6CS1AQu8CUixUEUREQEFBAbu8FhTsil0EbKgoShCkS2/SS0IaJUBCDYSEBNLr7rx/TECUlrK7s5t9PufkGLKTmd9u4rM3d25RVFVFCCGEZTjoHUAIIeyJFF0hhLAgKbpCCGFBUnSFEMKCpOgKIYQFOd3hcRnaIIQQxafc6gFp6QohhAVJ0RVCCAuSoiuEEBYkRVcIISxIiq4QQliQFF0hhLAgKbpCCGFBUnSFEMKCzFN0k47B2legIM8spxdCCLPa+iGc2mmWU5un6F4+BXtmw7G1Zjm9EEKYTcpJ2PoenLalotugF5SrDmHzzXJ6IYQwm4O/guIArR42y+nNU3QdHCFoJMRugtRzZrmEsC3Z2dl0794dg8Fwy2Py8vLo1q0bBQUFFkwmxHUMBRD+G9TvBd7+ZrmE+W6ktRoJqlF7AsLuzZs3jyFDhuDo6HjLY1xcXOjVqxchISEWTCbEdWI3Qfp5CH7UbJcwX9GtVBfqdoeDv4DRaLbLCOsyYsQIhg8fTrt27ahduzarV68GYMGCBQwcOPDacT169GDDhg0AvP7660yePBmAQYMGsWDBAssHFwK0LlEPX2jUz2yXuNPSjqUT/CgseRxOboP6Pcx6KfGPGSsPcfh8mknP2bR6ed7s3+yOx0VERDBw4EBCQkLYsWMHzz//PL179yYuLo46der8k3HGDN544w0SExM5ePAgK1asAKB58+bs37/fpNmFKJL0i9rN/w4TwcnFbJcx7zjdxveDWwW5oWYncnJySEpK4s033wSgadOmXL58mUuXLlGhQoV/HdutWzdUVeXTTz9l4cKF17odHB0dcXFxIT093eL5hZ2L+B2MBRBkvq4FMHdL19kNWgyH0B8hKwU8Kpn1ckJTlBapOURHR9OwYUPc3NwACAsLo2XLlri7u5OTk/OvY6Oiorhw4QI+Pj6UK1fuX4/l5uZeO4cQFqGqWldozQ7g18islzL/jLTgR8GQBxELzX4poa+IiAjOnDlDTk4OmZmZvPnmmzz33HNUrFgRg8FwrfBeuHCBkSNHsnz5cry8vFi79p/x3MnJyfj6+uLs7KzX0xD26PQuSD4BwaPMfinzF92qzaFGGwj9SXs3EWVWREQEQ4YMoX379rRt25aJEyfSuXNnAPr06cOOHTvIyspiyJAhfPLJJzRp0oRp06YxY8aMa+fYsmUL9913n15PQdir0J/A1RuaDTH/tVRVvd2HaYT+rKpvllfV07tNdkphfbp166YePXr0po+FhoaqjzzyyB3PMXjwYDUmJsbU0YS4tcxkVX3LT1VXPW/Ks96yrlpmwZtmQ8ClnPZuIsqs2NhYGjZseNPHgoOD6dGjxx0nRwwaNIhGjczbpybEv0SGgCEXWo+xyOUU9fZ/8puuP2DVc9pEiReOgntFk51WCCFKTFVhdkdw8YDxm015ZivYDTh4NBTkQOQfFrukEELc1tl9kHREq08WYrmiW70VVGslN9SEENYj7Gdw8YLmD1jskpZdxLz1GEg8BOcOWPSyQghxg+wrEP0nBD4Irl4Wu6xli27gUO1dRW6oCSH0FrkICrKhteW6FsDSRde1nPauEr0Esi9b9NJCCHGNqsKBeVA9GKoHWfTSlt8jrc1Y7d0lQpbvE0Lo5Mwe7QZam7EWv7Tli261ltoMtQPz5IaaEEIfB+aBa3mL3kC7Sp/dgNuMhUsx2nxnIYSwpMxkOLwMWj4ELp4Wv7w+RbfZEG2e84F5ulxemM8333xDtWrVaNWq1bWPqKgoq92ux8vL/HetrfW5263wBdoiXK0t37UAehVdFw9oNQKOrIDMS7pEEOYRFRXFO++8Q3h4+LWPwMBAu96ux56fu9UxGrXRUzU7QJWmukTQp+iC9i5jyNPWsBRlRmRkJK1atbrh67awXc+nn35K8+bNad68OZ9//vm1r7/99tsEBATQpUsXRowYwccff3zT75etimzAyW2QEqvLDbSrzLuI+e1Ubgy1u8CBH6HTFG0HYWEaa/4HCVGmPWfVQLjngzsedujQIcaOHYuDg/Z+/tRTTzFmzJhSbdfTtWvXm+4k8fHHH3P33XeX4kn9IzQ0lB9//JG9e/eiqirt27ene/fuFBQUsGTJEiIiIsjPzyc4OJjWrVvf9ByyVZEN2P89uFeCpoN0i6Bf0QVo+zgsHgsnNkGjPrpGEaV39uxZ/Pz8iIyM/NfXz58/f9vterZu3XrT7Xqu7iixffv2Yme5++67SUhIuOHr77777r9anVft2LGDwYMH4+mp3VgZMmQI27dvx2g0MnDgQNzc3HBzc6N///43vV5Jtyq603MXJpQaDzF/QafJ2q42OtG36Da+H7yqaO8+UnRNpwgtUnOIioqiSZMmN3y9tNv1lKSlu3HjxpI8hRKTrYpswNV1X9o8pmsM/fp0QdtxM3g0HF8Pl0/pGkWUXmRkJI0bN77h66Xdrmf79u3/ujF39cNUXQugFfZly5aRlZVFZmYmS5cupWvXrnTu3JmVK1eSk5NDRkYGq1atuun3y1ZFVq4gT1vcpmEfqFhH1yj6Fl3QFsFRHLS+XWHToqKi+OWXX64NFQsKCiIjIwOw/u16goODGTNmDO3ataN9+/aMGzeOoKAg2rZty4ABA2jRogX33HMPgYGBeHt73/D9slWRlTu6CjIuQttxeiex0HY9d/L7w6r6YV1Vzcu22CWFZdnydj3p6emqqqpqZmam2rp1azU0NPSGY2SrIis3715V/SxQVQ0Flrqiztv13EnbcZCVDIeX651EmIktb9fzxBNP0KpVK4KDg3nggQcIDg6+4RjZqsiKJR6B0zu0vlwrGCVlue16bsdohK/bgZs3jN9kkUsKIezEqufh4K/w/BHw9LHUVa1gu57bcXCAduMh/gDEh+qdRghRVuSkQsRCbS1vyxXc27KOogvQcoS2wPm+7/ROIoQoK8J/g/xMaPeE3kmusZ6i61ZeK7zRS2Q9BiFE6RmNWiPOv522R6OVsJ6iC1oXgyFPtvMRQpRe7GZtnQUrauWCtRVdvwCod5e25KNBlrgTQpTCvm/BszI0vXHat56sq+iC9q6UFg8xq/VOIoSwVSlx2kzXNmO1ma9WxPqKbqN+UKEW7P1G7yRCCFu17zttTK5OC5XfjvUVXQdHaDseTu+EC5F3Pl4IIa6Xm66Ny206CMpX0zvNDayv6AIEjwJnD2ntCiGKL/x3yE2DDhP1TnJT1ll03Stqw8ei/pDhY0KIojMaYe9cbcdx/zZ6p7kp6yy6AO2fBEOurD4mhCi6Exu1YWLtJ+id5Jast+j6BUD9ntoC5wV5eqcRQtiCvXPAq6rVDRO7nvUWXYD2EyEjQVYfE0LcWVKMNiGi7TirGyZ2Pesuug3uBp8GsGe2ts2GEELcyt654OiqbYxgxay76Do4aH0z58Pg7F690wghrFVWijZqocUw8PLTO81tWXfRBWj1MLhVgN2z9E4ihLBWB+ZBQTZ0fFrvJHdk/UXXxVObynd0NaSc1DuNEMLaFORpM9Dq94TKN+5GbW2sv+iCth6D4iCTJYQQNzr0p3bD3QZauWArRbd8dWg2BA7+oq0EL4QQoN1g3/01+DWG+r30TlMktlF0ATo+BXkZEDZf7yRCCGtxagckRGpTfpVbbktmVWyn6FYPgtqdtS4GQ77eaYQQ1mD3LPDwgRbD9U5SZLZTdAE6TYbUs3Bomd5JhBB6S4qBY2u1ez7O7nqnKTLbKroN+4JvI9j1pUyWEMLe7Z4FTm7aDDQbYltF18EBOk7S+nBO/q13GiGEXtIvalurtxoJnr56pykW2yq6oPXdePrBrq/0TiKE0Mu+b7V7OzYyTOx6tld0nd2g3ZNwYgNcPKx3GiGEpeVlaqsPNr4PfOrrnabYbK/oArR9HJzcZWqwEPbo4ALIuQKdpuidpERss+h6VNK29IlcBKnxeqcRQliKoQB2fwU120Ot9nqnKRHbLLqg3VBTjdqyj0II+3BoKVw5A52f1TtJidlu0a1YG5oPgdCfIPuy3mmEEOamqrDzC/ANgEb99E5TYrZbdAE6P6NNDd7/g95JhBDmdmITXIyCzlO04aM2ynaTA1QN1HaX2DsX8rP1TiOEMKedn0O56hA4TO8kpWLbRRe0vp3MJAj/Te8kQghzORcKp7ZrC19Z8f5nRWH7RbdOF6jRWpsabCjQO40Qwhx2fgau3la//1lR2H7RVRTo8jxcPgWHZSEcIcqcpBg4shLaPwGu5fROU2q2X3QBAu7VFjHe/qkshCNEWbPjc3D2gPYT9U5iEmWj6Do4aK3dxENwbJ3eaYQQpnL5NESGaN0Knj56pzGJslF0AZo/ABVqwfaPpbUrRFmx6yttf8SOk/ROYjJlp+g6Omnjds/t17bwEELYtvSL2vZcLR8C7xp6pzGZslN0AVo9Ap6VtdauEMK27ZkNxnzo8pzeSUyqbBVdZzdtS5+4rXDugN5phBAllZWiLd/YdJBNLt94O2Wr6AK0eQzcK8G2j/ROIoQoqT1ztCn+3V7UO4nJlb2i6+qlzVo5vg7Oh+udRghRXDmp2q7fTfpDlWZ6pzG5sld0Qdsd1M0b/p6pdxIhRHHt/RZyU6HbS3onMYuyWXTdvKH9BDi6ChKi9U4jhCiq3HTY87W2dGO1lnqnMYuyWXRBK7ouXjKSQQhbsv8HbX3sbi/rncRsym7R9agE7cbDoWWQeFTvNEKIO8nL1CZD1O8J/q31TmM2ZbfoAnScrM3Z3vah3kmEEHey/3vIugTd/6d3ErMq20XX00dbmejQUkg8oncaIcSt5GVqW/HU72mzG04WVdkuuqC1dl08ZdyuENZs33eQlQx3vaJ3ErMr+0XX00cbQiatXSGsU26GtglB/V5Qs53eacyu7Bdd0KYGu3hK364Q1mi//bRywV6KrkclaP+kNpLh4iG90wghrspNh51XW7lt9U5jEfZRdEFbj9O1HGx5T+8kQoir9syF7BTo+ZreSSzGfoquRyXo+LQ2S+38Qb3TCCGyL2vjcgPu1TaXtRP2U3QBOkwE94qw+V29kwghds3S1ljo8areSSzKvoqum7e2u8SJDXBmr95phLBfmZe05RubDYaqgXqnsSj7KrqgDR/z9IMt7+idRAj7tfNzKMi2mxEL17O/ouviCV1fgJN/aztMCCEsK+28NhkicBj4BeidxuLsr+gCtB4L5f1h4wzZOVgIS9v2ERgN0MP+Wrlgr0XX2U37gZ8PgyMr9U4jhP1IjtV2+G0zFirW0TuNLuyz6AK0eAh8G8Hmd8BQoHcaIezDlnfBybXM7gpRFPZbdB2doOfrcCkGIhfqnUaIsu9CBEQv0YZuelXWO41u7LfoAjQZANWDYMv7kJ+jdxohyrZNb4NbBeg0Re8kurLvoqsocPd0SDunLaAshDCPk9u18fFdngP3Cnqn0ZV9F12AendpCyf/PVObliiEMC2jETZMg/I1tIWn7JyT3gGswt0z4JtusOMz6P2W3mnsnqqqXEzLJf5KFgmpuVxMy+FKVh7puQVk5haQV2C8dqyjgwNero54uTnh7e5MlfJuVCnvRnVvd2pUdMfRQdHxmQgADi/V1jsZNAec3fVOoztFvf04VfsZxPrnk9pC51PCwNtf7zR2I7fAwJEL6USeu0LkuVSOXUwnNjGDzDzDv45TFPByccLT1QkXJweUwlpaYFDJyC0gI7cAg/Hfv64uTg7U8/WkfmUvAmt408Lfm8Aa3pRzc7bU0xMFefB1W3D2hAnbwcFR70SWcst3eym6V105A1+1gcChMGi23mnKLINRJfzsFXaduMSu2GRCz1y+1nL19XKhSbXy1Pfzor6fJ/6VPKha3o2q5d3wdnfG4TatVlVVSc8tIDEth4tpuZy7nEVcUiaxSRkcTUjn3OVsABwUaF7Dm471fehU35f2dSvh5mw3hcDy9n4Da16GkYuhYW+901iSFN0iWf+6tvLRhB1QtbneacqMnHwDfx9LYsPhi2w+mkhyZh4ATauVp2N9H9rUrkiLmhWo7u2GopinOyA5I5eo+FTCzlxhT2wyB89eJt+g4uHiSNeGvvRuWpXeTavg7S6tYJPJSYUvg6ByUxi9Esz0s7VSUnSLJPuy9ktSrSWMWmZvvyQmZTSq7DmZzLKD8ayJTiA9p4Bybk70CKjM3U2r0LWBLxU9XXTLl5VXwN6TKWw6cpGNhxNJSMvBxdGBHo39GNiqBr2aVMbVSVrApbLhTW1hmye2akMz7YsU3SLbMwfW/g8e/gMa9dE7jc1JTM/hjwPnWLj/DGdTsvFydaJvs6oMaFWdTvV9cHa0vgEzqqoScS6VFeHnWRl5nqT0XCp6OPNAsD8PtatFg8peeke0PZdPw6y22tKNQ77RO40epOgWWUEezO4ADk4wcZc2c03cUeS5K/yw4ySrIy9QYFTpUK8SI9rVom+zqjbVZ2owquw8cYmF+8+w/tBFCowqnRv4MK5LPbo38rttv7K4zuLH4OhfMPmAvd6YlqJbLEdWQchIuO8TaDtO7zRWS1VVNh9NZO62WPafuoyXqxPD2tTkkQ61qOdn+63DpPRcFh04yy+7T5OQlkM9P0+e6FqPIcH+uDhZX4vdapw7AN/3gm4v29XeZ/8hRbdYVBV+ug+SYrQhZG7eeieyKkajytpDCczafILDF9KoUcGdx7rUZVgb/zI5HCvfYOSvqAt8v/0kUfGpVPd244lu9XioXS2basVbhKrCvL5w+RRMDgNX23/zLSEpusV2/iB8e5c2T7zP23qnsQqqqrLpSCIfr4/haEI69Xw9eapHAwa2qm6VfbWmpqoqfx+/xKzNx9l/6jJVyrsypVdDhrWpaRfPv0iil2hdC/2/hNaj9U6jJym6JbLsaYgMgaf3gk99vdPoam9cMh+sPcrBM1eo4+PBc70bcX+L6nY742t3bDKfrI/hwOnL1KrkwQt9GtG/RXX77vPNy4Kv22mbvz6x1Z4mQtyMFN0SSb8IXwVD3W4w4ne90+ji5KVMPlhzhHWHLlLN241nejXkgdb+0rJDa/lujUli5roYDl9Io6W/N6/f35S2dSrpHU0fWz+Ere/B2DVQu5PeafQmRbfEdnwGG6fDqKXawjh2Ij0nny82HuenXadwdXLgqR4NeLxLXenDvAmjUWXpwXhmroshIS2HewOr8tp9TalRwY7WGUg9p83oDOgHD/6kdxprIEW3xApy4ev22mr3E3aW+SFkqqoVkPfXHOVSRi7D29Tk+T6NqFzOTe9oVi87z8C3f8cxZ9sJACb1aMC4rvXs441q8eNwdBVM2g8VaumdxhpI0S2Vo6th4cPQ70PoMEHvNGZzIjGdV5dGs+9kCi1rVuCtAc1oWdO+1z4tifgr2byz6jBrohOo4+PBu4MD6dzAV+9Y5nN6N/zYD7pPhR6v6p3GWkjRLRVVhV8GQ3wYTA4FLz+9E5lUTr6B2VtOMGdbLB4uTrxyT2OGtalp3zeFTGD78SSmLYvmVHIWQ4Jq8Np9TfDxctU7lmkZCuDb7pB9BSbtAxdPvRNZCym6pXbpOMzuCC2Gw6Cv9U5jMvtPpTB1cSRxlzIZXFgYfMtaYdBRTr6Br7ecYO62WLxcnZg+oBkDWlY328I+Frf3W1jzEgybD00H6p3GmkjRNYkNb8DOL+DxjVCzrd5pSiUzt4CP1h5l/p7T1KjgzvtDAunasGy14K3J8YvpvLwkkoNnrtCrcWXeHRxIVW8b7yfPSIJZrbXFbGSBqP+SomsSuRkwqw14VYHxm212HOLeuGReXBzBucvZjO5Yh5f6BuDpWrZvEFoDg1Hlx50n+Xh9DM6ODkzv34whwTVst9W7fBJE/K6tUeIXoHcaayNF12SiFsOSx+G+T6Ht43qnKZacfAMz18Uwb+dJalb04OMHW9Kurp2OKdXRqUuZvLQ4gv2nLtOnaRXeGxJoe106V9dXkBmbtyJF12RUFeYPgAsRMMl2bqpFx6fybEg4JxIzGNWhNv+7p7G0bnVkMKrM23GSmetj8HJ14sMHWtC7aRW9YxWNoQC+uwsyk7WbZ67l9E5kjW5ZdGVaUXEpCtz7iTblcf3reqe5I4NRZe62WAbP3kl6Tj7zH2vH24OaS8HVmaODwvhu9Vg9uQvVvN0YP/8Ar/wZRVZegd7R7mzft5AQBfd8IAW3BKSlW1Kb3obtH2vbkNTtpneam7qQms1zIeHsiUvhnuZVeW9woK67NYibyysw8umGY3zzdyx1fDz58qEgAv2tdGW7tPPa4uS1OsLIP+Tm2a1J94LJ5Wdri507umgz1Zysq5itO5TA1CWR5BUYmTGgGUNb+9vuDRs7sScumedDwknKyOXlvo15vEtd6xsrvWg0HFsLT+2BSnX1TmPNpHvB5Jzd4d6P4dIx2PWF3mmuyck3MG1ZNE/+EkrNih6sntKVB9vUlIJrAzrU8+GvZ7rSq3EV3v3rCGN+2k9Seq7esf5xfCMcXgbdXpSCWwrS0i2tRaMhZo02bMa3ga5RYpMyeHpBGEcT0hnftS4v9W0sOxzYIFVV+W3fGd5aeZjy7s58PryV/tOI8zK1v+yc3GHCdm0tEnE70r1gNukJMKsdVGuh6zbTf4ad4/Vl0bg5O/LJsJb0CKisSw5hOjEJ6Tz9WxixSRlM7tGAKb0a4qTXkprrXoPds2DsWqjdUZ8MtkW6F8ymXFXo8xac2g4Hf7X45bPzDLy8OILnF0XQvIY3f03pKgW3jAioWo4VkzozNNifLzef4JEf9pKYlmP5IOcPwp7Z0HqsFFwTkJauKRiN2p5qiYe1pe28LFP0ru9OmNSjAc/erWNLSJjV4tBzvL4sCi9XZ758qBWdLNXdYCiA73pARqK2g4q7rDpXRNLSNSsHB+j/BeRnwZqXLXLJ1ZEXGPDVDi6m5fDj2La82DdACm4ZNrS1P8uf7oK3uxOP/LCXWZuPYzRaoE20exYkRMK9H0nBNRFp6ZrS3zNh8zsw/Fdo0t8sl8grMPL+miP8uPMUQbUq8PXDwVS3px0K7FxmbgGv/BnFiojz9Ajw47PhrajgYabhipeOw5zO0LC39jstI2CKQ26kWYQh/58/xZ7aAx6mXdfgQmo2Ty8II+zMFcZ2rsMr9zSR0Ql2SFVVft1zmrdWHaZyOTfmPBJMC38Tt0KNBvjxHkiK0boVylU17fnLPulesAhHZxg4G7KSYZ1pV9DfdeIS93+5g5iEdGY9HMSb/ZtJwbVTiqIwqmMd/pigbf44dM5uft93hjs0oIpn33dwdi/0+0AKrolJS9ccNr+jdTU8/Ac06lOqU6mqytxtccxcd5R6fl7MfSSYBpVlvrvQpGTm8czCg2w/fokHW/vz9qDmpd+TLeUkzOkEtTvLVN+Sk+4FiyrIhW+6QU4aPLUL3CuW6DTpOfm8+EcE6w5d5L7Aanw4tAVeslCN+A+DUeXzjcf4avMJmlUvz9xHWlOzkkfJTmY0ws/9tVX0nt4D3v6mDWs/pHvBopxcYdAcyLgIa6aW6BTHL6YzcNZONh5J5PX7mjDr4SApuOKmHB0UXugTwA+j23AmJYv+s3aw7VhSyU62dy6c3qGtICYF1yyk6JpLjWBtjnpkCBxeUaxvXR15gYFf7yQtJ58F49ozrms9WTtB3FGvJlVYOakLVcu7MebHfXy5qZjDypKOwaYZ0KgftBppvqB2TroXzKkgT1tdP+28NprhDgueFxiMfLQuhm//jiO4VgVmj2xt+/toCYvLzjPwyp+RLAs/z91NKvPp8FaUd3O+/TcZCmBeH0iJ035X5eZZaUn3gi6cXGDwN5CbBque1XaduIXkjFxG/bCPb/+OY1SH2ix8oqMUXFEi7i6OfDa8FdP7N2VrTBIDZ+0kJiH99t+04zOID9W2oZKCa1bS0rWEnV9oOwkPmAXBo254OPzsFSb+GkpKZh7vDg5kaGvpSxOmsf9UCk8tCNN2fx7agvtbVL/xoPhQ+KGPtoX60HmWD1k2yegFXRkNMH8gxIfBxB1Qqd61h37fd4Y3lx+icnlX5j7SmuY1rHTHAGGzLqbl8NSCMEJPX2Z817pM7df4nynjeZkwt6s24mbijhKPtBE3kKKru9Rz2thH30Ywdi05RoXpKw6xcP9Zujb05cuHgmQrHWE2eQVG3ll9mPm7T9OhXiVmPRys7UC88hkI/blw26muescsS6ToWoXoJbD4MdLav8AjsT2JPJfKU3fV54U+ATha27YsokxaEnqOV5dGUdHDhQXdkqm/cTx0fgZ6v6V3tLJGiq61SPx5ND4nV/CY+iYPD3uIvs3kpoWwrOj4VKb9sp4fsp+F8v5UnLINxVlu2pqYjF7Qm6qqzNkaS6+j/UlwqMIP5ebSt650JwjLWXYwns4fbGbgV38zLfcz3JUChl4ax0vLYsjJN+gdz25I0bWAtJx8Jvwayodrj9I9sB6VHv0Vp6xLsHzSbYeRCWEqyw7G88qfUcRfyeZpx2UEq4eYYRhLw6ZBLA49xwNzdnE2JUvvmHZBiq6ZxSSkM6hwOu+0+5vy1Ygg3Ou0gd4zIGY17P9e74jCDsxcF0N2voE2HGGK42L+NHRhYX4XouPT+GF0G86mZHH/VzvYcjRR76hlnhRdM1p2MJ5BX+8kPbeA38a15/Eudf+ZztvhKWjYV1sC8vxBfYOKMu/8lWwqqqk4rp3OV5HuTMsfCyicv5KtTR+e3IXqFdwZ+9N+Pl0fg8ESu1LYKSm6ZpBbYOCN5dE8GxJOYA1vVk/uQvt6Pv8+SFG0RXE8/bRt3LMv6xNW2AV/bxd6Rb/B4Qs5rG0wlUy03Uau7jpS28eTpU91YmhrbRPMMT/uIyUzT8/IZZYUXRM7dzmLYXN3M3/3acZ3rcuC8e2pXP4Wd4Y9feDBnyAtHpY9Lf27wmxecAjh561xdBg0ihinRgC4OzvyUt+Aa8e4OTsyc2gLPhgSyN6TKdz35XZCT9Eseq8AABnkSURBVEtjwNRkyJgJbTmayLMh4RiNKjMfbEm/5kUcDrZ7Nqx7BXq/DZ2nmDeksDvpEatp22sAEwZ14I/6b3E+NYfqFdx5qW8Ag4Jq3PR7ouNTmbgglAtXcnjl3iY81rmOrHRXPDJO15wKDEY+2XCMOVtjaVKtPHNGBlPH17PoJ1BVWPQoHF0Njy6XmUHCZNTLZxjVqymubu78sPUUuBT99zI1K58X/ohg45GL3NO8Kh8ObXHn1crEVVJ0zeViWg6Tfz/IvpMpjGhXkzf7NyvZdik5adoykFkp8OQ2WUBalF5+Nj883pLP1p1k3/5QPGq1KPYpVFXlu+1xfLg2hhoV3Jk9MljWBykaKbrmsO1YEs+HhJOVZ+C9Ic0ZHFTKQpl0DL7rCb4NYewakFlCoqRUlegvH6LHq3+wLWQ2Te+fUKrTHTiVwqTfDpKSmcfr9zdhVIfa0t1we1J0TSnfYOST9ceYuy2WgCrl+HpkkOk2izyyCkJGQtAj2lKQ8ostSiBj61e0ffBZ/vf4YEZ/sNgk50zJzOOFReFsiUmiXzOtu8HbXbobbkGKrqmcu5zFMwvDCT19mRHtavFm/6al3331v67uJnzPR9D+SdOeW5R9J7czevDdKOWr8dPWOHAw3SAlo1Hl+x1xfLQ2hqrebnw5IojgWrIc5E1I0TWFv6IuMHVJJKoK7w5uzsBWN7/zW2pGI4Q8AsfWwMjF0KCXea4jyp6Uk/z0ZFs+2pHJ/uhYPH1usmi5CYSduczk3w6SkJbD870bMbF7fRxkpbzrSdEtjay8At5edYTf952hZc0KfPVQELV8SrjFdVHlpsMPfbV1eMdv0vp5hbidnDQOv9OF7p8fYsvaVTTvco9ZL5eanc+rS6NYHXmBTvV9+HRYK9li6h9SdEsqOj6VKQsPcvJSJk90q8eLfQJwdrTQnJLLp7Uba+4V4PEN4FHJMtcVtsdoIPPHB2n3ygpeeOEFHpv6oUUuq6oqiw6c5c0Vh3BzduTDB1rIcqUaKbrFZTBqQ2U+WR+Dj6crnw5rSacGvpYPcno3zB8A/u1g1J/g5Gr5DMK6qSr89RKPvf4l+dVaM3/1LouPLDiRmMGzIQeJjk9jRLuaTLu/KR4uThbNYGWk6BbH2ZQsXvgjgn0nU7ineVXeGxyo71Y6kX/An+OgxXBtd2EZ0SCut/trfvnoRd494MmBI6fw8vLSJUZegZFPNxzjm79jqV3Jg0+Ht7Lnm2xSdItCVVWWhMUzfcUhAKYPaMYDwTWsYzzi3zO1UQ3dp0KPV/VOI6zFkZUc/fphus7PZ/P2vQS2bKl3IvbEJfPCoggupGYzqUcDJvdqaLkuOeshRfdOEtNzeG1pNBsOX6Rd3Up88mBLalYy882y4lBVWDEJDv4KA76C4Ef1TiT0dmYvWd/3p/28LKa8/iHjJzytd6Jr0nLymbHiMEvCztG8Rnk+ebAVAVVNNJbdNkjRvZ1VkeeZtiyazDwDL/UJ4LEuda1zo0hDPvw+AmI3wfBfofF9eicSekk8AvP6MX5FJpn+3VkQstg6/iL7j7XRCby2NIr0nAKe7d2QJ7rW+2f797JNiu7NJKXn8sbyaNZEJ9DS35tPhrU03cwyc8nLhJ/7w8VDMGop1O6kdyJhaVfOwg99+C3sCtP3uhN6MIJy5az39zY5I5fXl/3z/9lHQ1vaQ6tXiu71VFVlWXg8M1YeJivPwDO9GvJkNxt6B85Mhnl9ISMRxv4FVZvrnUhYSuYl+PEeYuLO0uXnPDZu3kJLK+jHvRNVVVkZeYHpKw6RnpPP5J4NmXhX/bLc1ytF96qzKVlMWx7N1pgkgmtV4KOhLay/dXszV87AvH5gyNMWx5HJE2Vf9hX4uT/ZF2LosKg8E595kQkTSreQjaUlZ+QyfeVhVkacp3HVcrw/JJCgsjnCQYquwajy486TfLL+GIoCL/YJYHSnOtbZd1tUl47Dj/eAo4tWeCvW1juRMJfcDPh1CMSHMeFIB66oXvz+++9W2Y9bFBsOX2TasmgupucwumMdXuwbgJdrmRrXa99FN/zsFV5bGsWh82n0bFyZtwc1p0bh3lA2LyEafrpPm7U2dg2UN89ce6Gj/Gz4bTic2kGI1+O8PncZoaGhlC9fXu9kpZKek8/MdTH8suc0Vcq58Wb/pvRrXtVm30j+wz6LbmpWPjPXH2XB3jP4ebnyRv+m3BdYraz8UP9xLhTmDwSvyjBmlRTesiQ/WxuxEreV48Fv0unx91i/fj1BQUF6JzOZsDOXeW1pNEcupNEjwI8ZA5qbf20T87OvomswavPBZ66L4UpWHmM61eW53g0pV5a3Gjm7D34ZIoW3LLmu4Ob0+4yOE79k3LhxPP209YzHNZUCg5Gfd5/m0/Ux5BtVJnSrx8S7GuDuYuJlUy3Hfopu2JnLTF9xiMhzqbStU5EZA5rTtLpt/xlWZNcKrx+MXgXeZlp6UphfXhYsHAFx22DQbJ7+bieJiYksWrSo7P2ldp2E1BzeX3OE5eHnqVHBnVfvbcK9gTbZ5VD2i278lWw+XHOUFRHnqVzOlVfvbcLAVtVt8YdVOmf3azdc3Ctom1xWqqd3IlFcOWlaH+6Z3TBoNotPuDB16lTCwsLw9raP/cn2xiXz5opDHE1Ip12dSky7vymB/jb13Mtu0U3PyeebbXF8tz0OgPFd6zHxrvp4lq07ocVz/qDW4nV00Qpv5cZ6JxJFlZWivWkmRMGQb4l1b0XHjh3566+/aNOmjd7pLKrAYCTkwFk+WX+My1l5DAny54U+jahuGzfBy17RzSsw8vu+M3yx6TgpmXn0b1mdqf0C8K9o8x3wppF4BOYP0sbxPrIEagTrnUjcSdoF+GUwpMTBsPnk1ulBp06dGDNmDJMnT9Y7nW7ScvL5evMJftx5CkWBx7rUZUL3+ta+P1vZKboGo8qKiHg+33ic08lZdKhXiVfvbUIL/wp6R7M+KXHaqIbMZBg2HxrerXcicStJMfDrA1pLd8TvUK87U6ZMIT4+nsWLrXNdBUs7m5LFpxuOsfRgPBU8nJnQvT6jO9ax1ptttl90VVVl3aGLfLbhGDEX02lSrTwv9w3grgA/+YW8nfQEWDBUa/kO+ApaPax3IvFfZ/ZofbiOLjByEVQP4s8//+TFF18kLCyMChWkQXG96PhUPl4fw9aYJPzKuTK5ZwOGt62Jq5NVFV/bLbpGo8r6wxf5YtNxjlxIo56vJ8/3acS9zavJRnhFlZMGi0ZB3Fbo8Rp0e0kWQrcWh5fDn09A+RpaN1ClusTFxdGhQwdWrVpFu3bt9E5otfafSmHm2hj2nUqhmrcbT91Vn2HWU3xtr+gajCqroy4we8sJjiakU9fXk8k9GzCgZXXbWZjGmhTkwYrJELkQAh+EAbPAWTYR1I2qwvaPtYXp/dvBiIXg6UNeXh5dunTh4Ycf5tlnn9U7pdVTVZVdscl8tuEYB05fpmp5N8Z3q8eIdjX13i7IdopuTr6BpQfj+WZbLKeSs6jv58nTPaTYmoSqwvZPYPPb4N8WHvpNm0whLCs/B1ZOgcgQCBymdfsUvgE+99xznDx5kqVLl0q3WTFcLb5fbDrOvpMpVPRwZmznuozqUFuvrbasv+imZObx657TzN99iksZeQTW8ObpHg3o07SKdCOY2uHl8OeT2u7Cw+aDv30NRdJV6jkIGQXnw6Dn69D1xWtdPcuXL+eZZ54hLCyMSpVk5+eSCj2dwuwtsWw6moi7syPD2vjzWJe61PbxtGQM6y26RxPS+HnXKZYejCcn30iPAD/Gd61Hx/o+8k5vThciIOQR7UbbvTOh9Ri9E5V9J/+GP8ZCQS4MngNN+l976NSpU7Rv357ly5fToUMHHUOWHTEJ6Xy/PY5l4fEUGFXublKFsZ3qWKq2WFfRzSswsv5wAr/sPs3ekym4OTswOKgGYzvXpVEVG1zb1lZlpcCSxyF2M7R6BO79CFws2hqwD0Yj7PoCNr0FPg1g+ALwa3Tt4by8PLp168aDDz7ICy+8oGPQsikxLYf5u0/z274zpGTm0aiKF490qM3goBrmXI/FOoru6eRMQvafZdGBc1zKyKVGBXdGdazN8DY19d3i3J4ZDbDlPa2v17cRPPgjVGmmd6qyIyMRlj6pvbE1HQQDZ4HrvxsWL774IjExMaxYsUL+ujOjnHwDKyPO8/PuU0THp+Hh4sjAVtUZ3rYWLf29Tf3a61d0s/IKWBudQMj+s+w9mYKDAj0bV2Fkh1p0a+hn24uIlyWxW7ShS7lp0PddaPO4DCsrrRObYOkE7TXt94HWhfOf13TlypVMmjSJsLAwfHx89MlpZ1RVJfJcKgv2nmZFxHly8o0EVCnHg238GRRUA18vV1NcxrJF12BU2XniEksPxrPuUAJZeQZq+3gwrE1NHgj2p6q3DFWyShmJWpGI3QT1e2mtMlkisvjyMmH9NDjwA/g1hqHzbvrXw5kzZ2jbti1Lly6lUyfZYFQPaTn5rIq4QMiBs0ScvYKjg0K3hr4MDvanT9MquDmXeMyvZYtuyP4zTF0SRTk3J+5vUY3BQf60rVNR/nSyBUajViw2vAGOznDPTGgxTFq9RXV6NyybCJdPQcentREKzjcu0JKfn0/37t0ZNGgQL7/8suVzihscu5jOn2HxLA+P50JqDqM61ObtQSXe9NWyRTc1K5+dsZfo2bhyad4phJ6SY7XicXav1uq97xOoVFfvVNYr+zJsnA6hP0GFWjBoDtTpcsvDp06dSlRUFKtWrcLBQcafWxOjUWXPyWQql3Mtzaa11nEjTdgYowH2f6/ddTca4K6p0OFpcJKbnteoKkQvgbWvQNYl6PAU9Hj12iiQyMhI3vvgI2Z9+Tm+vr4ArF69mokTJxIWFnbta6LMkaIrSiH1HPz1MsSshkr1tRttjfpJl8P5g1qxPbMbqrWC/l9A9Vb/OmTmzJlMe/9Tyrk6sPSPEOrUqUObNm1YvHgxXbrcuiUsbJ4UXWECxzdoRSb5ONTrAb1nQLWWeqeyvCtnYesHEL4APHyg1zQIGgUON3alPTD8Ybam+eHoVZGsjbPw8S7HE+PH8dprr+kQXFjQLYuudCaJomvYG57aDX3f11p533SDRaPh0nG9k1lGRpL2pvNVMEQt0m6UTQnThoLdpOAChEdE4lK5Dh7121JhxMck5jmzfvNWkpKSin35L7/8kiZNmjBy5MhSPhGhJ2npipLJSYVds2D311CQDc0GQ+dnoVoLvZOZXuo57bmG/QwFOdqaxN3/BxVq3vbb8vPz8fQqR9VJC3AoXNBGNRSQues3OPE3f/4RQrdu3Yoco3HjxmzcuBF/f/87HltQUICTkx1vWaU/6V4QZpKRBLu/gv3zIC8dGvSGDhO17gdbvyt/Phz2fQuRi0A1akPnujz/rym8txMdHU3X3vfjPfrrf33dkJ1GyoqP4NIJLicn4+Jy5xuTEyZMYN68eQQEBDBmzBi2b99OXFwcHh4efPvtt7Ro0YLp06cTGxtLXFwctWrV4tdff2Xq1KmsXbsWBwcHxo8fz+TJkwkNDeX5558nIyMDX19ffvrpJ6pVq1ail0jc0i2LrrwVitLx8oPeb0GX57SRDnu/0TZWrFQf2o6Dlg9pq5nZirwsOLoK9n0H5/aBswe0GQudJmtDwYohMjIS58p1rv3bkJ1OVugysiPX8uDQobz15rIiFVyAuXPnsnbtWrZs2cKMGTMICgpi2bJlbN68mUcffZTw8HAADh8+zI4dO3B3d2fOnDmcOnWK8PBwnJycSElJIT8/n8mTJ7N8+XL8/PwICQnhtddeY968ecV6bqLkpOgK03CvqO1I0WmKtnTkvm9h3SvaJItGfaHFcGjYxzoXTjca4PQubYH3Q8u1FnulelrfdauHte3sSyDsYDj53jUxZKeTHbacrIg1PDBkCG//Gk6dOnVKHHfHjh0sWbIEgJ49e5KcnExaWhoAAwYMwN1dm4yxceNGJkyYcK2boVKlSkRHRxMdHU3v3r0BMBgM0sq1MCm6wrScXLU/w1sM07YRj1io/Xl+dBU4e2qbYza+X5tw4anjWgN5mXBqBxxZCTFrtDG2Ll7QdKD2BlGna6m7R9Zt2kp2ci654SsZPHgw78w/SN265p1g4ul5+1XiVFWlWbNm7N6926w5xK1J0RXmUzVQ+7h7BpzcphXeo6u1lvDVx+vdBbU6QvVgKG/GFlf2ZW3Exdn9Wpaz+8CYDy7loFEfaHyfNvbYhEtbnj93hi5BQXw7Zyn16tUz2Xm7du3KggULmDZtGlu3bsXX15fy5cvfcFzv3r355ptv6NGjx7XuhYCAAJKSkti9ezcdO3YkPz+fY8eO0ayZrCxnKVJ0hfk5OkGDXtrHvZ9ouybEbYG4bVof8K6vtOPKVdcWhvELAN+GUKG2tuBOuWrgdmNRuUF+NqRfgLQLkHoWLh3TPi4ehpTYwoMUbWxxx6egbndtqq6TSVaVukFyYoJZzjt9+nQee+wxWrRogYeHBz///PNNjxs3bhzHjh2jRYsWODs7M378eCZNmsTixYuZMmUKqampFBQU8Oyzz0rRtSAZvSD0lZ+jdUPEH4D4UEg8qk2+KMj593GOLto6tC5e4OT2z2w4Qx7kZkBuujZ07XoOTlrfrG8jqBEMNVpD9SBw87bMcxP2TIaMCRtiNGot1dSzWqs1/by2y0VuOuRl/LsgOziDq5dWkN28tdZy+WpQ3l9boMfRbDsDCHE7UnSFEMKCZBqwEEJYAym6QghhQVJ0hRDCgqToCiGEBUnRFUIIC5KiK4QQFiRFVwghLEiKrhBCWNBtJ0coirIWKOl2pb7ApRJ+rzlZay6w3mySq3gkV/GUxVyXVFXtd7MH7jQjrcQURTmgqmobs5y8FKw1F1hvNslVPJKreOwtl3QvCCGEBUnRFUIICzJn0f3WjOcuDWvNBdabTXIVj+QqHrvKZbY+XSGEEDeS7gUhhLAgKbpCCGFBJiu6iqLMVBTlqKIokYqiLFUU5ab7ViuK0k9RlBhFUU4oivI/U13/NrkeVBTlkKIoRkVRbjn8Q1GUU4qiRCmKEq4oygErymXR16vwmpUURdmgKMrxwv9WvMVxhsLXK1xRlBVmynLb568oiquiKCGFj+9VFKWOOXKUINcYRVGSrnt9xlko1zxFURIVRYm+xeOKoihfFuaOVBQl2Epy3aUoSup1r9cbFspVU1GULYqiHC78//GZmxxj2tdMVVWTfAB9AKfCzz8EPrzJMY5ALFAPcAEigKamynCLXE2AAGAr0OY2x50CfM2Zpbi59Hi9Cq/7EfC/ws//d7OfZeFjGWbOccfnDzwFzC38/CEgxAKvT1FyjQFmWer36brrdgOCgehbPH4vsAZtZ4MOwF4ryXUXsEqH16saEFz4eTng2E1+liZ9zUzW0lVVdb2qqgWF/9wD+N/ksHbACVVV41RVzQMWAgNNleEWuY6oqhpjzmuURBFzWfz1KjQQuLrF7M/AIAtc82aK8vyvz7oY6KUoyi23SrFgLl2oqvo3kHKbQwYC81XNHqCCoijVrCCXLlRVvaCqaljh5+nAEaDGfw4z6Wtmrj7dx9DeGf6rBnD2un+f48YnqBcVWK8oSqiiKE/oHaaQXq9XFVVVLxR+ngBUucVxboqiHFAUZY+iKOYozEV5/teOKXzTTwV8zJCluLkAHij8c3Sxoig1zZypqKz5/8GOiqJEKIqyRlEUi+8JX9g1FQTs/c9DJn3NnIoZaiNQ9SYPvaaq6vLCY14DCoAFJQ1VXEXJVQRdVFWNVxSlMrBBUZSjhe/Oeucyi9tlu/4fqqqqiqLcalxh7cLXrB6wWVGUKFVVY02d1UatBH5XVTVXUZQn0VrjPXXOZM3C0H6fMhRFuRdYBjS01MUVRfEClgDPqqqaZs5rFavoqqp69+0eVxRlDHA/0Est7Az5j3jg+nd8/8KvlcqdchXxHPGF/01UFGUp2p+QpSq6JshlltcLbp9NUZSLiqJUU1X1QuGfUYm3OMfV1yxOUZStaK0EUxbdojz/q8ecUxTFCfAGkk2YoUS5VFW9PsP3aP3k1sBsv1OlcX2hU1X1L0VRZiuK4quqqtkXwlEUxRmt4C5QVfXPmxxi0tfMlKMX+gEvAwNUVc26xWH7gYaKotRVFMUF7caHWe56F4eiKJ6KopS7+jnaTcGb3mW1ML1erxXA6MLPRwM3tMoVRamoKIpr4ee+QGfgsIlzFOX5X591KLD5Fm/4Fs31nz6/AWh9hdZgBfBo4R35DkDqdV1JulEUperVvnhFUdqh1SZzv3lSeM0fgCOqqn56i8NM+5qZ8C7gCbR+j/DCj6t3lKsDf/3nTuAxtBbRa+a4I/mfXIPR+mBygYvAuv/mQrsLHVH4cchacunxehVe0wfYBBwHNgKVCr/eBvi+8PNOQFThaxYFPG6mLDc8f+AttDd3ADfgj8Lfv31APQu9RnfK9X7h71IEsAVobKFcvwMXgPzC36/HgQnAhMLHFeDrwtxR3GZEj4VzTbru9doDdLJQri5o93Mir6td95rzNZNpwEIIYUEyI00IISxIiq4QQliQFF0hhLAgKbpCCGFBUnSFEMKCpOgKIYQFSdEVQggLkqIrbIqiKG0LF5FxK5xJeEhRlOZ65xKiqGRyhLA5iqK8gzYTzR04p6rq+zpHEqLIpOgKm1O43sF+IAdtuqhB50hCFJl0Lwhb5AN4oa3076ZzFiGKRVq6wuYo2n5sC4G6QDVVVSfpHEmIIivWerpC6E1RlEeBfFVVf1MUxRHYpShKT1VVN+udTYiikJauEEJYkPTpCiGEBUnRFUIIC5KiK4QQFiRFVwghLEiKrhBCWJAUXSGEsCApukIIYUH/B2MrQRCj8CEfAAAAAElFTkSuQmCC\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",
    "np.random.seed(42)\n",
    "\n",
    "xspace = np.linspace(-2, 2, 100)\n",
    "unnormalized_probs = np.exp(-0.5 * xspace ** 2)\n",
    "energies = 0.5 * xspace ** 2\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(xspace, unnormalized_probs, label=r\"$p(x)$\")\n",
    "ax.plot(xspace, energies, label=r\"$E(x)=-\\log\\ p(x)$\")\n",
    "x_index = 75\n",
    "ax.scatter((xspace[x_index],), (energies[x_index],))\n",
    "prop = dict(arrowstyle=\"-|>,head_width=0.4,head_length=0.8\",\n",
    "            shrinkA=0,shrinkB=0)\n",
    "a_start = np.array((xspace[x_index], energies[x_index]))\n",
    "a_end = np.array((xspace[x_index] - 0.5, energies[x_index] -0.5 * xspace[x_index]))\n",
    "ax.annotate(\"\",a_end, a_start, arrowprops=prop)\n",
    "text_pos = (a_start[0] + 0.5 * (a_end[0] - a_start[0]), a_end[1])\n",
    "ax.text(*text_pos, r\"force\")\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_yticks(())\n",
    "for spine in ('top', 'right', 'left'):\n",
    "    ax.spines[spine].set_visible(False)\n",
    "ax.legend(frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you now find a way to predict where the marble will end up given an initial position and velocity, you can use the result of that prediction as a proposal state for some kind of fancy Metropolis-Hastings move.<br>\n",
    "This kind of prediction problem is actually very common and well-studied in physics&mdash;think of simulating the movement of planets around the sun, which is determined by gravitational forces.\n",
    "We can thus borrow well-known theory and tools from classical (\"Hamiltonian\") mechanics.<br>\n",
    "While HMC is quite easy to implement, the theory behind it is actually far from trivial.\n",
    "We will skip a lot of detail, but provide additional information on physics and numerics in detailed footnotes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In simple cases like ours, the force acting on a particle can be calculated as the derivative (or, in higher dimensions, gradient) of a _potential energy_ $E(x)$. The negative log-density plays the role of that potential energy: \n",
    "$$\n",
    "E(x)=-\\log p(x)\n",
    "$$\n",
    "Our marble does not only have a potential energy, but also a *kinetic* energy, which is the energy stored in the movement of the particle.\n",
    "It depends on the particle's momentum&mdash;the product of its mass and velocity.\n",
    "In HMC, the mass is often set to one, which motivates denoting the momentum with the letter $v$ from now on.\n",
    "The kinetic energy is then given by $K(v)=\\frac{|v|^2}{2}$[^1].\n",
    "Just as the position $x$, the momentum $v$ generally is a vector, whose squared norm (\"length\") is denoted by $|v|^2$.\n",
    "We can then define the total energy of the particle as\n",
    "$$\n",
    "H(x,v)=K(v)+E(x) \\mbox .\n",
    "$$\n",
    "HMC introduces the momentum of the particle as an additional variable and samples the joint distribution of the particle's position and momentum.\n",
    "Under some assumptions[^2], this joint distribution is given by\n",
    "$$\n",
    "p(x,v) \\propto \\exp\\{-H(x,v)\\} \\mbox .\n",
    "$$\n",
    "The sampling itself is being done efficiently by means of a fancy Metropolis-Hastings method:\n",
    "starting from a current state $x$, to obtain a proposal state $x^*$, a random initial momentum $v$ is chosen[^3] and the movement of the particle is numerically simulated for a certain time chosen by the user."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the simulation stops, the fictive particle has some position $x^*$ and some momentum $v^*$, which serve as our proposal states.\n",
    "If we had a perfectly accurate simulation, the total energy $H(x,v)$ would remain the same and the final position would be a perfectly fine sample from the joint distribution.\n",
    "Unfortunately, all numerical simulation is only approximative.\n",
    "As we have done before in the original Metropolis-Hastings algorithm, we can correct for the bias with a Metropolis criterion.[^4]\n",
    "In our case, the acceptance probability is given by\n",
    "$$\n",
    "p_\\mathrm{acc}^\\mathrm{HMC} = \\min\\left\\{1, \\exp\\{-[H(x^*, v^*) - H(x, v)]\\} \\right\\} \\mbox .\n",
    "$$\n",
    "We are not interested in the momentum $v^*$ and store, if the move is accepted, only the position $x^*$ as the next state of the Markov chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## copied and pasted from Gibbs sampling blog post\n",
    "def plot_samples_2D(chain, path_length, burnin, ax, lims=(-5, 5)):\n",
    "    chain = np.array(chain)\n",
    "    bins = np.linspace(lims[0], lims[1], 100)\n",
    "    ax.hist2d(*chain[burnin:].T, bins=bins)\n",
    "    ax.plot(*chain[:path_length].T, marker='o', c='w', lw=0.4, \n",
    "            markersize=1, alpha=0.75)\n",
    "    ax.set_xlabel('x')\n",
    "    ax.set_ylabel('y')\n",
    "    ax.set_xlim(*lims)\n",
    "    ax.set_ylim(*lims)\n",
    "    \n",
    "## copied and pasted from MH blog post\n",
    "def p_acc_MH(x_new, x_old, log_prob):\n",
    "    return min(1, np.exp(log_prob(x_new) - log_prob(x_old)))\n",
    "\n",
    "def sample_MH(x_old, log_prob, stepsize):\n",
    "    x_new = proposal(x_old, stepsize)\n",
    "    # here we determine whether we accept the new state or not:\n",
    "    # we draw a random number uniformly from [0,1] and compare\n",
    "    # it with the acceptance probability\n",
    "    accept = np.random.random() < p_acc_MH(x_new, x_old, log_prob)\n",
    "    if accept:\n",
    "        return accept, x_new\n",
    "    else:\n",
    "        return accept, x_old\n",
    "    \n",
    "def build_MH_chain(init, stepsize, n_total, log_prob):\n",
    "\n",
    "    n_accepted = 0\n",
    "    chain = [init]\n",
    "\n",
    "    for _ in range(n_total):\n",
    "        accept, state = sample_MH(chain[-1], log_prob, stepsize)\n",
    "        chain.append(state)\n",
    "        n_accepted += accept\n",
    "    \n",
    "    acceptance_rate = n_accepted / float(n_total)\n",
    "    \n",
    "    return chain, acceptance_rate\n",
    "\n",
    "def proposal(x, stepsize):\n",
    "    return np.random.uniform(low=x - 0.5 * stepsize, \n",
    "                             high=x + 0.5 * stepsize, \n",
    "                             size=x.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "In terms of implementation, the only not immediately obvious part of HMC is the simulation of the particle's movement.\n",
    "All simulation of the movement of physical objects works by discretizing time in discrete steps of a certain size.\n",
    "Position and momentum of a particle are then updated in an alternating fashion.\n",
    "In HMC, this is usually done as follows[^5]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def leapfrog(x, v, gradient, timestep, trajectory_length):\n",
    "    v -= 0.5 * timestep * gradient(x)\n",
    "    for _ in range(trajectory_length - 1):\n",
    "        x += timestep * v\n",
    "        v -= timestep * gradient(x)\n",
    "    x += timestep * v\n",
    "    v -= 0.5 * timestep * gradient(x)\n",
    "\n",
    "    return x, v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we already have all we need to proceed to the implementation of HMC!\n",
    "Here it is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_HMC(x_old, log_prob, log_prob_gradient, timestep, trajectory_length):\n",
    "    # switch to physics mode!\n",
    "    def E(x): return -log_prob(x)\n",
    "    def gradient(x): return -log_prob_gradient(x)\n",
    "    def K(v): return 0.5 * np.sum(v ** 2)\n",
    "    def H(x, v): return K(v) + E(x)\n",
    "\n",
    "    # Metropolis acceptance probability, implemented in \"logarithmic space\"\n",
    "    # for numerical stability:\n",
    "    def log_p_acc(x_new, v_new, x_old, v_old):\n",
    "        return min(0, -(H(x_new, v_new) - H(x_old, v_old)))\n",
    "\n",
    "    # give a random kick to particle by drawing its momentum from p(v)\n",
    "    v_old = np.random.normal(size=x_old.shape)\n",
    "\n",
    "    # approximately calculate position x_new and momentum v_new after\n",
    "    # time trajectory_length  * timestep\n",
    "    x_new, v_new = leapfrog(x_old.copy(), v_old.copy(), gradient, \n",
    "                            timestep, trajectory_length)\n",
    "\n",
    "    # accept / reject based on Metropolis criterion\n",
    "    accept = np.log(np.random.random()) < log_p_acc(x_new, v_new, x_old, v_old)\n",
    "\n",
    "    # we consider only the position x (meaning, we marginalize out v)\n",
    "    if accept:\n",
    "        return accept, x_new\n",
    "    else:\n",
    "        return accept, x_old"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, analogously to the other MCMC algorithms we discussed before, here's a function to build up a Markov chain using HMC transitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_HMC_chain(init, timestep, trajectory_length, n_total, log_prob, gradient):\n",
    "    n_accepted = 0\n",
    "    chain = [init]\n",
    "\n",
    "    for _ in range(n_total):\n",
    "        accept, state = sample_HMC(chain[-1].copy(), log_prob, gradient,\n",
    "                                   timestep, trajectory_length)\n",
    "        chain.append(state)\n",
    "        n_accepted += accept\n",
    "\n",
    "    acceptance_rate = n_accepted / float(n_total)\n",
    "\n",
    "    return chain, acceptance_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, let's consider a two-dimensional Gaussian distribution.\n",
    "Here is its log-probability, neglecting the normalization constant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_prob(x): return -0.5 * np.sum(x ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to calculate the gradient of the log-probability. \n",
    "Its $i$-th component is given by $\\frac{\\partial}{\\partial x^i} E(x)$, where $x^i$ is the $i$-th component of the variable vector $x$.\n",
    "Doing the math, you'll end up with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_prob_gradient(x): return -x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to go, that is, to test the HMC sampler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance rate: 0.622\n"
     ]
    }
   ],
   "source": [
    "chain, acceptance_rate = build_HMC_chain(np.array([5.0, 1.0]), 1.5, 10, 10000,\n",
    "                                         log_prob, log_prob_gradient)\n",
    "print(\"Acceptance rate: {:.3f}\".format(acceptance_rate))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To display the result, we plot a two-dimensional histogram of the sampled states and the first 200 states of the chain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3RcaXnwf3d606j33iWrWLbl3u1dby+BBTaQhNBDgJCQL3zAF04gpBCSUFIgLARCW1jYXWd3Wbz22usu23KRLKv33stIo+nlfn+MrZmRLFnr9a6L3t85Oudqbnvnzsz71Pd5JFmWEQgEAsHKQ3G7ByAQCASC24MQAAKBQLBCEQJAIBAIVihCAAgEAsEKRQgAgUAgWKGobvcA3gwaSSvrMN7uYQgEAsFdhZWpcVmW4+e/flcJAB1GNkp7b/cwBAKB4K7isPx8z/VeFy4ggUAgWKEIASAQCAQrFCEABAKBYIUiBIBAIBCsUIQAEAgEghWKEAACgUCwQhECQCAQCFYoQgAIBALBCkUIAIFAIFihCAEgEAgEKxQhAAQCgWCFIgSAQCAQrFCEABAIBIIVihAAAoFAsEIRAkAgEAhWKEIACAQCwQpFCACBQCBYoQgBIBAIBCsUIQAEAoFghSIEgEAgEKxQhAAQCASCFYoQAAKBQLBCEQJAIBAIVihCAAgEAsEKRQgAgUAgWKHcdgEgSZJSkqQaSZJ+e7vHIhAIBCuJ2y4AgM8CTbd7EAKBQLDSuK0CQJKkNOAR4Ie3cxwCgUCwErndFsC3gc8D/sUOkCTp45IkXZAk6YIH1zs3MoFAILjHuW0CQJKkR4FRWZYvLnWcLMvPyLJcKctypRrtOzQ6gUAguPe5nRbAVuBxSZK6gV8BeyRJ+vltHI9AIBCsKG6bAJBl+YuyLKfJspwFPA28IcvyH9yu8QgEAsFKQ3W7ByAQCO48FHp92P9+h+M2jUTwdnJHCABZlo8Bx27zMAQCgWBFcUcIAIFAcHsI1fRDtXyh8a8MbncaqEAgEAhuE8ICENxzLKbVChYins/KRlgAAoFAsEIRFoDgnmO5Wu1yM11ERozgXkVYAAKBQLBCEQJAIBAIVijCBSS4o3gn3S3LvbZw+QjuVYQFIBAIBCsUYQGsYO7EdMk7ZRwCwUpAWAACgUCwQhEWwApGaNs3hzImem5bdjjD9olnKribEBaAQCAQrFCEBSC4a7hTYhahWv+tGMed8r4EKw9hAQgEAsEKRVgAghtyp5RCuJn7hvrrYWntfbma+M2MY/4zfKvXu1O4U74bgptDWAACgUCwQhECQCAQCFYowgUkuCF3s1nvm5xadN9S7oul3EG3wp2z2L1u9no3y1sNQN/N3w2BsAAEAoFgxSIsAMFdw52SLinpdXPb8y2Mt9pj9+14X0s9t+VaPYJ7E2EBCAQCwQpFWACCu5Llph+qUpLDj5uy3PAcAIqyg/dq7grbNb/8Qyg6tcSGhyqoPlCLk8U1alVOVnDf0Mii1wsdv3dwaNHjlqvlL4XQ+lcewgIQCG4hH/jrd/On3/4wGx6quN1DEQhuiLAABG+at9NXvJRm/2buZTAb5jRx+4x90esvhr+mcW57/mKyUOZrUF3dkxz89TkGxl1h451/DW9n9/X3zdfeQyyW5WYL3YrFWWKB18pAWACCe5IHP7SLP/vux/j9LzxJRIzpHbmnKcrAxPAMP/mn35KcFUdUQuQ7cl+B4GYRFoDgTfNWNco3c/7NnKfQ62m80M3P//F/OfKTo2SVpBMZF4Hb5aW1ro+ZidkbXkNZXDC37WtqXfy4EO09f0MBzSfr8c86OfXsCXY8vmbOApF0ukXPWyqmsBS3IoPnrWYtCe5uhAAQ3JPYZhxcOFTH9NgMV8ZmAFBrVBRvLyEyzozH5aHpVCPT4zO37J56kw7HbHAyP/Gbs+x5/zZOPH8W/y27i0Bw6xACQHBPIvv9KNThX2+P20t9VUCbV2tU5JWmUrqtCK/HS0t1O5axmxcGkiQhI4e95vf7Of7rKna+bwsnjnXg9wkxILizEAJA8I5ws6mIiwU+5wdVQxdkSXodskaDUqdekAYqOwMauhdoON0MBIRBQWUu5Y9txuvx0Xq5l8m2vhuOAYLum7SCZPrre8MWieFw4HF7Ob2/mm2PV/LGs6fmdoWmgYa9L324q2i57qE34x67XYR+ZkuV6BC8c9y2ILAkSemSJB2VJKlRkqQGSZI+e7vGIrj38Pv8KJTL+3p73F4aqlqoeq2OSyeaSc6MY8tj69j4YAXRicsL5KblJ9PfNnzdfY5ZJxcP1bHjqU3LHv98DGY9O5/aiMF8+yZwwb3H7bQAvMBfyrJ8SZKkCOCiJEmvy7LceKMTBfcWSy2KUkRHBfeFaLlLlWCQM5PxJ5oh0oD/YkfYcUulZqLT4gWaGgbB6UKlVpJfnkNRhRevx0vbxU4mhy0o1qwKnnNtkZgEsiyjiAoRGCFjtNo8tNb1s+nJTVQfvAxOV3C8S2j518b7nv/3JI9+ch//+flfcnz/+cDORbT+N1Nc7p20CITWf+dx2wSALMtDwNDVbaskSU1AKiAEgOAt4/fLKBTSW7qG1+Oj6WIX3sEhlColBetyKN5UgD8jhbamQSbHA9lEepMOl911g6vBaN8EWoOW8u1F1Ndf31q4Hpsfr+TYc6ex2j3MTN44g0kgWC53RAxAkqQsYA1w7vaORHArWSzFcL7WGbooaimNdNkdu2oa8SVFQWrMkhrvfM07VEMNjR0o9HpkoKWuH+r6UaoayavIpCg/BrksFZ1Rw6XjzQFrxbm0IOhrGURn0FKwIY/Wur4F7yX0vpIksWPXBuqrWpgcnqbnW6+ydm8psdF6xgenFn1W0rw4Qtj7X+LY0H03mwa63M9oqeNErOCd47YvBJMkyQS8APy5LMsL0jAkSfq4JEkXJEm64OHGWpZAAFctgGXGAN4sPq+PlgudVL1ykfNHG8kqSqFsUx4b7y8lLjnqhue31XSh1qjIyE9c9BilSsGuJ9dRc7SByeHpudcvHamnbHsRas0dobsJ7nJu67dIkiQ1gcn/F7Isv3i9Y2RZfgZ4BsAsxcjXO0bwzrJczX5Jf3vo9RY5580gZaXPbaumrShjTKgS4oD28HuFxBTQacP2edNjguOobrr+OYA/Mfhe3DWNNB6p4czLF1AoFRTtW0/B6gx8Pj/tFwOaOoAUFR5MvvLKWSrvL8Oh9jEWEgORLdNo9Rq2PVTBqZeqsQ+OB0/atBqAU13TbPzUI1T9z/Hg2EMKxSlDngWA3N3HYrydGvbNlvUQWv87x20TAJIkScB/A02yLH/zdo1DcG/i9/vfcgxgOSRmxjPaE5ik/T4/rZd7ab3ci0KpICc/noKKTPyyTGe3hbEhS9i5F16/wtbHK3F3TDI9aQPAaNaz8YFyTrx0AY/be917ulxeWpuHKN2YS/25juseIxAsh9tpAWwF/hC4IklS7dXXviTL8u9u45hWLG+m+NfNaOmh/valfNTz8/YXK4G8wP89Mhbc1unw+2UkhbTwuBCtX7ZMh+1ShWQghY7Rm5EQfo2zl+c2c//wPi4ca5ob9zXrwA+01fbQVtuDQqkg6+FKcrYW4JdlOtpGmLh676rX69n14b1UH29Bo1WxujCeYy9U4/fLC56H3No7tz3W2kvyhmziTEpG+ybChjdf41ckB11NC95zyL7QDKzl+u/nI0pI3F3cziygU8Dbr6IJViTy2xgDCEWpUtxwha/f56etZZi2lmEUComcvEQKH67AL8t01fdz4rV6nv7YTgb7Jjj2oyPLvvfl403sfGojU6MzuMTEK7gJRCRJcE/i9/tRSG+vfqHWqhd10yyG3y/T3jpMV2svCoVEdkkau57YSHJ6DJIkoVDeWKCEUvXKJbY8tpajPzv2JkcvEAgBIFgGyzX5l7vIaKmFT/NdPoulBErzA526kK9yzxB+yY/k9S4Yk6TVBLfnVegkMT64HeJSUo2Gu0246pYpWJNJR8NA2C5HWrD0tKEn+LrC5QkfR1QkMiCZI2is6eHXPzzBqjUZfPF/Psmxl2vpbB5kuG8SOXS8LMQLNF/uY+1TO6iragu8Nu8ZKkJSU+e/58W6kS33M59f0sIfktJ7JyL6HIQjBIDgnsTv8yPdZBA4LtHMU5/dw/P/cWgui+d6xCSYaanpWXT/jShZl4XH4+PKhW4AGi71MNAySNnGXHQGLVvuL8WvUtLdPsJQ7+Si1xntnyQxN5mkjFiGeycW7FcoFag1KjRaFWqtCo1WffV/NUpXLBpdwJJJSI/l/KE67DMre1JcSQgBsMJYrA79m9GEFjt2vna1WLBXzgwP9Eo9QY11vtc+uSyLj37hcX76rQP0XGibe92Zbg47TnuyIfhPVjp+lQJFbDT+nWvDjvNog3dQR4aPd6rYAMCH3r+Lx3dXIClV/NffvLBQS77aL9gfZVqQSmq8GBKADdG2xyvC00D3KbKwTNoY6J0kQqecm6BVJj0eWWLrIxVcPteBVimxY28xKZlx+I1aurrGGBwICCVViHCSZZk/+MuHqT5cj9NeghTi/vJMz+BxevC4PHhQ4HF5cNvs2Cc8eNxeouIi2P5kJQ9+YBsAx58/t+zvg3eexn+nL+Ja6Rr/fIQAENzRPPnB7Wy+fxWy7OfwLwIT9FDnKK0SyEusCpHfQimIn712kXyVgf6OEaITzKgiVGh0atRaFWqtGnV+Bslp0SSnxYDVET4Od8DVEzoBy7LMdHnC3Os7KnIZqepFoRgnIlKPxxqoGup2e3H7/LQ3DmCzOjBG6Lj0RgMXT7QELp0eQ3Z2PCmp0SDL9HeMMNAVdFXVnGxh68OrqT4SIgwJbyt5PaG8ZtcqXvzuYXrrezl/qO6mnpng7kSSl/oV3WGYpRh5o7T3dg9jRRCafrhYKuY7cb2kLav56BcfZf+PTtB8vB6AlJwE0tcXIkkSLoebtqZBZkNSJOXMZAwmLZ/6/MNUH7yCyxHifw+ZoGWjntDvv98Q8Lf7/TJS3whPf+4ROuv7eHFqBMusE7fXi9vrQ93rYm15BtWXupAnwlena/qCrppQ//21e27bV0LDxR6mBoKuGnfIAjSPOaiTFeYloht10No4CIBiOtjbGCAj0UBqbiKyDL1XuulvGyYuJZrMygJqTrYExxEZMbc9v7tZ+Z88xuDAFONjVpShz3CJOE2oEFlKy79Zf/vb2XN6pXJYfv6iLMuV818XFoDgjmZ8eJqvf/YXbHuonOarrw12jjJoC0zcOoOG3KJkInIDroeRnnG6HTLrt+RTsT6HC4fqqT4aXNUbVoUzZGIE8IW4hJStvazeXsSTn9iL9VQdz/z2LLNONwAmhxdkcHt8qJf5PpQqBTseKOP8qVZmpx3Lyn9uaR9hQ0IsmTnx9HSOLdjf0zJET0tAmGZkxbLlsXUAmFOjSc2OD7MOrofepMMYoWN8zLrMdyG41xAC4B7gzWhay9Wulqulh15vvnshVDu81ogFlm7mMv+9XNOixyftxK4vZHw0UC7Kaw7cywPU9k+gcPkASEqJZV1WHLJGxaEDV2hqGg7XxOdp5aEop4PPw741n1/VdJHTPMh0xzTrdbHoozT4/H6c8iyP7VxFV3UvA+GhCKZ3Bi0djykwzevUKu5bnc/+unYc642AkZQXuuaOU18J2Q7NRAKuHL5C5b5yHLYZxiy2sH2hGT1dGiNdvYGJPNOj4P1feoILVW20Ng7SFVJBVDMdHN/GRyo489OjKD2BZ7fcmFBonGapz/JmudUF6gSLIwSA4K6gqbaXbY+Wc/Lw0tXChwctDA8GfN5ntC0UpEaTXZAEwOiQhe62kWXn2U9N27lQ2019Qz8J8WaOn2pBqZD4xleeYk1lFl/8qo745Eie/+0lBq6WeXCZg1OWTydh0mvYWZLDwdpWKrKTcXp8ON0eIvMdOB1uXA43rmkZl8Oz6LguHKpj2xOVuC/1MD1+Y229r3ucf/nKfnbcXwLA1q35IEn0908y1BYQ7Ol5iQz1jOO9OvkLViYiBiC4ITdbnmGxEgTzc9FDr7egVPSG4rntCpWPwc6RsOqYc9cM0RpDNdnQMSSkRJG6NR+lUoHb7aXrUi/TU0Gt2rIqmKkTfSngo9+4qwidXk1bwyCxCRFcudDNA+9aR3ZhEht3FJGUHk171ygf+YufAmBsGJ27hjk5moqNOZw8WI+lOBq1SolWo0KnVaE2aYLbBhVatQqlUoFuInxhmXo64LJSKCR23VfCuWMt2G0uJAm8Hj8upxuXy4urfwS304PP62PVzlLOH29Bq1eTU5wyF0ROzY4nLUaDUqUkpyyDn37thWUX7AtlOQ1sYOmssKWsvuVcW/DmEDEAwV1Pw5k2tj1ZyclrHbHeJKODFvrPBoqnaTQqCtOiWVWeBsDYyAw1Sie+62jhTbV9xMRH4PfLpGbGMmOxc2j/JTqbh3j603v41+8dWnBOQkoUeWsyOX7gytxrHq8Pj9fHrN2Fzx4M6PpCvFKmgWDA2mjQsLUwjfNVbdhnXVQfbmL3w+VcqmrH4/aiVCnQatVo9WqSzXq2PrIGWZZ55MO7ATj+6mVmU+xk5CXQ2z7KQNcYfVUjVO4rp/q1WjY/uhY8Hgbah+muDwaABSsHIQAEdw2yLDM9biUqwYxldEHriDeF2+2l6Ur/3P9xCWY27s1BpVTg8foY7XUzOW5FkgKuo8KyNE4fbuCBd1eiN2hw2Fwc2n+JM8PjmE3hFk1aTjxJaTGcvoG7ailioo28512VPLKvnO//62scfvUyfp+fEwfr2fFgKUdfrcPn9WOI1VK8OoPJ1j6e/ddXMUToGRqZ5fzxgNbfcrmPLftKGR+exj7rIjIuAr/PT/P5gCD0Oxyk5iWx5Yn1KCNMDHWP0dU4sNTQBPcQQgAIbkhoABcWN9kXmOiLVJecf/5ii9MANG0h7qboKBobR9jy0FpOvnA27LjQevu+goy5bUd8+EIt7VgwC0gxEnRDTI5M0ZQV8N+rVUpKdmdQFmlk7dpcumMUnJ6cIfcPKrDG6NmwJo9fnmnElxfFlMvGjtIkXBMBl1HJtkIMRi1na3pQhpST0MxbuKaeCWr6rqtjjDTrWV2aglqtZHLCRsvFXnbvKCK9PIWc6Vk6T3XgBqpOt/Guj2ynzeBhxGLlhc4hou0+KMtGHWfCF2nAG2/Gawhc9+yrNex8fC1H919g7VNbOPHq5TnXmGSZZnjSxfDJdmSHk5ScBDbtWYUyNprJcSuJqdGcP9bMbNvyVjwv5c5ZLEAsXDu3DyEABHcVfr+MdcqGOcb0tvXH9Xh9NHQEevb6ZZmxqVkKMuP5vT3lmPRavv3d19m2KZ/jVfNy6vNTUPbYqXsT5SEiInRUrM1Co1ExPWPnUlUXHo8PtVrJv/3gQyQmRFK6Kh2LxcETu0s5sP8i8Ulm+nvG8RVFUNcZFJCrV2fw2ONrKCsLuLWOXbVA/D4/l04084HPPcSx39WxVNxvfHCKqHgz0TFR3Pd7a9l8XynIMkeXKQAEbw8Gs4END1VQfaAW+4z9xicsEyEAVjDLTQmdr5UvFjhUXC2RMEdzML0xNJAcujIVwoPC8+/lKE2d29a3BxqvXLncz+b7Kzh9MOhf90cagtta5dy24XSwfASEl6EYeyArbF98dXBcrvjA9TTTXqZGZqmzOEmOjGB1USrWDVoMGVH88d49nNnfgnPGzYOVhXT0jtNgnYK0wHM1ErRE3CELvIx6DZXFsWi1aqxWJ+dbB3B7AsFfjV4BegVbthfy9e8c4BNPb+GZb/wOl8OL/NmtfOlb76O6pY+fHb6I7ltHeagohSsnm9nyB7vpahni5PMezrzewIWqtrACdHarM1ALKEIflhbrK8ggNS2ajKw4FD1DuJ0eOup6aKrrp/ZgDXVHGzh/uJ7Q4nvS9OKZSLdi0eBi3KuF3CRJwhRlxBxrwhwbQUSMaUEp89yKLB79xD4Ajj1XdcvuLQSA4K7D7/Njn3ViitQzO/32TwIqpYJdGwtoaBtCpVTQZ5+m/vwIlampxKgUPLijmNgoI88fqMFmdzEysXCCNOjUlOSnoNepsdldXDnWicsVmPTdaeETW3FhMl1dY/T2T/LvX32ZwrI0RocslGYlERdlwueXuW9tAd5piezSdPY+vYVvf/VlSipzuFzVzqxSueD+lbuLefbbB1nzYAUOu5uk1Gj0Ri0+g4bB/knOVbVDTXD1sKTXYbc6OH4t4G4wLbimYGkkScIUbcQcGxGc2K9TnkSWZWYtdmbGZxjqHKHlQseClOC6E00MtA1RfaB2wflvaYwiDfTeZrla080uJgtlqU5foaWWl+pYNb/wWmgaqGomJL1TqWDzjkJOX13lG1qGOdTPb5s3ufq0wR+gMyb8x6gMLPTFpNPwaGImp061sroiAznDwNGadjatyuRUXSePryrA4/NTfamTP/v8g3zv9bNkxkVxsqmb1cSQGBdYYTwV6UOlUqJWKnDaPdR1D2G/qpVHdgfTPUPHG6OQKCxLp6djlILSVLp0fhrbA+6o9+woo3JtJv/23cOMjc/yyXdvQqVWMjNlQ69VoVIr+fG/HMAdG77COTE7lo0bcxkbnUEzOs26XcX899+9hG3GEV7MLnSV9Ly4z2Kpn0t95ktZjktxp5eCUCgUIRO7CVP09Sd2v19mdmoW6+Qs0+NWZi22N9Xr4VYi0kAF9xR+nx+n043RpMM2u3he+s2wpTSLP3nPHuITzaSnxfAPr5zAdbXxi88vU1mRxfo1WRw7nco/7D/K1sIs7C43MjLdAxMYDRoijDpiIwy4vT48Ph8alZJoo35OAFwPSYJH3reR8ZFpcorKefnZswxnGgFQSBIutwelQoHBoOX+PRn878+rcDk9VG4rICUlErVPxYbdxchRRqYmbZgitKjUSnbcX8J//+A4IyPTaNqGqKtqo3RjLuder7+lz+1uR6FUEBFtmnPFmKKNYUX9rnFtYp8et9LXMojNYsfvvz0T+1tFCIB7kKU0qJsp17vchToLGqyEIIdkxMzXGuf3qQ0lNAsotDevK17L2aExNj9YzMmzbShc1/8BGvvnvf8Ql1H7H8fOOzbwYz/Z2s3Xf/kGJdlJRKVHMpnoQ6nQMmH2MJsucdTST447kW83X8DVaadmvIsvfWQfpZHx9O6vpfn5Gjotdmzrgn7z2Rw1uUmx7CrNAcBS4KKxfwSnxwv5KtJiIvnkvk08P9zIpq1pPFRYysVVLrpfbAdgXX4aj99fTnFpGgqfzN984Te4jCrQaznU0MOHdu2mKDcZS4KO8dphsnMTsNmcZOcm8OxXXmSqawzN1Wc4BQz7JTLvK6PnVHgg+xrzP8vQLKtQS2Epn/9yvzfzWSxj7Ga/ywqlAnNMYFI3x0VgjDRcf2L3+bFOzjIzYaW3aYBZi23JgPm9gBAAgrsWr9eP2+1Fr1fjcrlufMIymXW6mbLaef7oZdYXp/OxreuZtNt5ua4Zk1ZDdlw0XRNT7C7IwaG04/X6+Olvq9mzvoALoYXnQvD5ZVoHx2kdDASyIyJ1rMlKIT8ljqKUBIYtVnrGpkCGswP9NE2Mcby3i7XZycRFGSnPTqa5cRCjScd3vnEA2S+TkRZDelqglaTJoOUr/3WATz+9g8HzA7z60iX0Bg0KhYLkzDiyipIZ7B6n3RawZDraRti8LZ/JKANWy63LKnknUaqUREQbiIwxYVAmLjqx+7y+qxP7LN31fdim7ff8xL5cRAzgDuZm/ffLOWf+ecv13wNILvfc9lLa+2Ja2YKSAyHXn9/oRV9//UVJ17KD1GolG9Zlc+GVYB37axk8EJ59M5/R9eGTRcL5wG+hND+ZkTQfLSPj7C7IofG3XXzzS+/iYkMf0WY9VZphHsrJx+H18tVvvjZ3/p+/azvPNl+heyKQTaQIPiYefl/4uoWJv9lCXlEyg32TVFV4+cdH97EtJ5Oj+2v46b8dpmRtFgqFRDU2vD4ff/XBvVwaHybWZKBhcASXx4vlzDh9g5PERptI1htobgqUjX7fZ3YyPmVDIcHLJxqIPBV4hsmZsWQXJuPz+miq6cFqcbD9j7fxxhuN+P3ykqWsQyulcvby3KZizaqw40LXVixoTblM375KrSIiJuiKMV7N8JI04WPy2h3MTAQ09pkJK7bpu1OQvROIGIDgnsTj8eHz+9Hq1Lici/vXl0teZjxIEi0jAU1doZDIzYgDYF1JOrN2F9UtF7C6XWSao1iTl0pNe2CCPV7XycaCdCZsDqzO61skcZpkUvWZnDFoOPF6oHHLrqc28u8nzpAbF0PL5X6yC5I4f7Il0Dzmz3eyKjcJq81JfISR315uxuZy0zNhIW4goM0X5SVx/o1gBs//vHSOT753G263j4rCVLquCoChngmGeiZQqZUUVWRgjjYyO+ti2/ZCThxvXjjYW4xaoyIixkhEjIkIkwZDxPUVF6/HNzepd1zumct7v1fTQG8nwgIQAEv/uOYXgwvN41/SwgjRDkPbPob5kwm3Iua3i5zfBOUaI/enzG1r1Ur2KeM4e7XomTU3aEXox9xh5w1tDlo6cQ3hhdcM22NJjY7kTEcvW0fNpKREsXFjLo1D4+h1Gn7x8nlSE6N4//0VqFRKfvV8NVv+soTv1J/AJ8volWo2u/KINuh5tb6Z9OeCcQ/50/kUJcbTPTlF08gYpqtVKHaU57BKG8n6ddnUXO7lYskBMgypaBRqKqJLePHHQ0xbnahVSmYTZFqHx1mXlUpVey8JlwLvbfvWAi78JCQ3PNrIfY+sxufzU3exm7LSdK5c6GJ8ZCZckweiDRoefHItBoOGgz88OtdTeEEWUMjnEvZZ6nSBiT3agDnGhMHrmJvYF2jsHl9AY5+0Yp2cxW4N3ENM5G8/wgIQ3LO4PD5khfyWrIC4GBPbywroHJtkb3EuUz0DnDrZSnZ2AqmJUXz3FycAGBix8MyPjvMPX303WzflYfO62ZqYzYnhThw+DxqVkjdaOrivKI8WxkjLjiM7P4lzKi8HQrpx5abEsr4wnbzkOLLNZi7W9JCcHIVJZaTfMUR+RDY/7nqOlvosdqzNxefzc2GoF4/Pj1YV/NmmpkTT3x/eMH7tplzOHG/B5/ezYWs+R1+qpaZ7YNwAACAASURBVHx9NjlFyZy70offF1T6pqdsPPfjk2zcXkBKdjw5Jak47W4aTjXh9/kxRxuJiDYSsSoTnS7Q/kYqCrrsJLUGj9vLzJSNmSkbA42dOK5mZd1sEFjwziEEgOCeoPZ8FxXrszl38vpZLYuRnhpNaVEq6yuy+OaFc3SPB3zYyYNOSq7GGerbhvD5g5Omw+nhv39ykrTUaEqjslkTl8qF8T7sXg9+Wcbu9uD2+fjT//cYB35TzclD9XTGxaFTqdiYmUZJciIqi5/CtARMeg0HDzTw2utX2LwxF8N79ZhUBk6Ph1c8VSoVeK6TQ56dFcepqjauTbWRMcZAuYwZx9xz2XrfKk4fbsQYoWPHvhL6eieYnrITGWUgUq9Gp9cgSbD3fZs48fIlIiINfPjL70KhlGi/3EvN8SYG6wdwOgIWx3wLIBTvLU7JFby9CBfQO8Q7ubjlVt9rKfdQaEB3/mQQ6ioKXdAVVuCN8PRO1ei8oPIivnR3friraGizjj3FOZxq6yH2ZLDG/2xaeDE4U7+L5JQocrLjUU/YGBuykJ6XwNGXauj8QPCaj08mYJmyEZ8Qwce+cIxpx+/m9n3mM39GdKSBpMRIFI2j/MmXH2dmysY3/++vKXrfOorykjh6ugVVioHJGRvREQb2PViCX5ap7umnfmiESnsMZpMevU7Nr393Ca1GxRf/ZB/fOVCF1e5kS2kWVfXdZD3dSmlkIR6/h9c/F5hct2zKo+psO55IFdsr8zhxvh1rhgKdWsXvb62g8Y0uIkw6dFoVsgzJiZGkp0RTXdONdshGZl4CeoOGwy/VMFbdiOvqxG7Kz6BiSz5Vh+rn3HJJmXHklKTh1+tprulhetKGrST4eWlfqQ57vjfz3Xs7FisKwhEuIME9z9mOPjblpNN2cmFAMyk2goLMBPRZXoYGpzh9ug39kIVdj67h+Ku1YWmBuwtzGHmxn8GBKdZtyMbpbV9wvalpO8UFyVxuGeK5773B+z65h6888yFermqms2ecCKOO7esLyUmN5Ux9N39/6BhObyDmUJwYT7wugsgIPT/df4705GiKcxPpHpgkNyWWXatz2Vmewz9zDI3RjlGpp9c1RX5eKuYIPevXZqFWK0nLi2fW7mLruhwc8QryEmPpm7AwZbHR0zcxt3gNICs9FnOEjtpjzZw71oxGq2L99gL0M9O01QRqNtlnXfS2j1BUkUHTsUCdpeGecYZ7xlHGRlO8JpOSSiMziUYiTFrOXOjEu+DJCO4mhAB4m7id2sqtWHK/VIqoKjrq+vdNDE/vDC1nFZZiOC/AGJpW6N20OmxXaCmH0IDufM1eNwV+PBiTlIzuMeLx+Ug0m9iqikehUDA0McPxhk7Gy/0QD9JqDR9ybONwVRvOzMC4nTkuHs4s5MJoH9mfqycGiE2r4CPVTzLsCNb38T0YSB8tyJLwrSrGEhVLd7yHBL2RfQmB8R+qb+Mb9oPYGt08llmKV2lBkmUi1Fruy0snfTqW5v4x/v6vn0COkGidHKc4LoHBKjtag4qXLzbh1cDult04HW7qX7AynO6mZWSGvlN2/H4Zf6vMydMBl1dEXQ9JD5TT8Mol2v9v0dWRBmoC5X+nm7GLfURXZJAbpaS9thsXcKqli/SyLHa9fzsXjtRjd7rob+xn7c4iYkpymAxpP+luauVyWydavYZPfv9P2XRfKZqOCY7dgp7AC1KQF/n+Co3/1iMEgOCGGCL0rL+vlPOH67Fb7+wfYefwJF94eAdnOnsZsc5yurorzH9/jQfzCqj+cSfOq6UZlAqJJ7KLOT7YicXl5FpdU4/fg1a58GeyPjmVj5Sv48RgF90zFvyynymXA9eUF5NGw7vWleAbddBnsxCjNfLjbX9In22KdGNgwlSolKhVSn527BI1xlEsTic7M7JIlo383YtvcH95Pqdbuom0WvD5ZLpaR5gyBVYuD09YWVecjtca7OdbuTfw+SxFa20vZWXJpBck09cacMP1tw8z2DVK5Z4SnF6ZK2fauXS8mV3v38bJQ/VhPYONkQY2PbqOn37jt1ypauP80ZtveCO4MxAxgDuYW7EQbD7L7QEbWsjr4U8/zEf/3xNcOdfBkeercQ4F+t7KsozTYGZ6yobVYmeqpims2FWYFRFSTtiTEF5ZMtQ6cOTFhe27VgIaCPM9z6YEJ+Uok56CNUmolUrGrDbKjo9z+uWL+Lw+ev5qTdj1DCOwvTSblv4xHBofH9mznmdP1bImO5XXBtrx+v0YNBr0eX5yzDHsSMqlb3yWfusMaxKTWZ2QhNkUfIb20cAiJavDRYRei8IDibERpCZE0jkxSe3wML+6XIfPJfPu1SVM2AKrUHuaxjl5pROfX0bhllEqJD7/gT388t+Os3lLHh/7+G5+8MxRLKs0nOvpo3vSQuarwdjGR/90Dy/+5BSWSRsJSZHE5MRTXx/ILQ0tmgcgOQOOmoSUKOI/sgq/X6aub5ihaSuZvx6cOy420Uz5xlwun+3ANjDK+r2lnHz5IgARG4upWJfJiTeakJqDvQGW6g+8oORHyLFLNQhShFiYb2d56ZWEiAEIbpoTL1/CYXXSdLGLzMJkVGolgx0jdNf3YV5diDnKSHJGLDlx5ShVyuByfLV67hoOYwQzFjtWi4MJlYTX+9aKZ5kNWkqzk9GpVUzNOjjd0YvHF9BWzU0jfO3Fv+RX//wyESmJGDRqDBo1EhLbd6QzPGlFrVbyRw9UkpUQzZbCTDpHp0jLiEKvUuPyeZk12dmeFKjbU2gKVCt1er38rL6W3KxBVJIKo8pA6+lIHG4PCkliTU4qboubobFp9m4s4Atnj+D0etmRncW4xU5ubAy5sTH805ET9FweDXs/960v5MDZZpRKBePjVg4dusLUlA2nR6J7Mrx/AkB0tBHLZEAglFRk8HpV24JjQilfn82OB8tY/1AJ3zhwnLL0JBye8JTZiZEZjr5cQ/mmXLSZUXQ1DrBqQy4TQxayy9LmGswsLDYtuFu5rRaAJEkPAt8h8J36oSzLX1/q+JVmAbwdLDdL40ZFuJKz4slelYrf6ablYieTwwsnqVB0W8owR+oDf2Y9KlVgGpEkUAyOI8sykiRhi45kZtrOzLSDmWkH7pBAptkvU7Q6HZ1OzXC0ioa24blGKqqQYnCfe89WHn//ZjqaBzl0sJ7e7nFUaiUlZWlMpCpoG5ng4fJCWuuH2LE2j5+8Uo1apSR2h5mWqTGUCgXqcTWrU5K42D/IjpJMRuw2Jux2SuITaG+bYMw6i1+GyBwdv2y5jAxIwLulVSgkCbfXy+uDXfhlGZ1axfc+9AQAR/o6eK61jld2B7/q3zr8fWq7BzHrdXhVMrkJMfhlUEigsUkcuxzo35v4ekBb1+rUfPSLj/IPL51hzao0+oYsjE8Fu6PZksObidyvSqSgIIkjhxtYnZc012R+90NlHG/ow3E1C8gRHxTYGZdGWLezmPIteYDE9z/9zNyK3KUsx9AYwHKzx5a6huDWcMdZAJIkKYH/BO4H+oHzkiS9LMuycCzeBQx1jzHUPYbkdlOwNpviDXk4bE4azrTisrsXHG+3ubDbXAwPLhQUytbeuW316lzMkQYiow1kZMURFW2kuCQVvUGDfXKWliv9qFRKstOiyE4L+MStNhfuWRfXdJmOpiFqzrTT3TrM4MAUySnRxMaZaG8d5jdt/TxcXsh/HD5DqSKGv/vhQbKSY5CRaZkaI8FgotUyjnJWRbzJyNCMFavbzYfL1/Gfl87xlVNvoO4J1hDalJTBNRVKBmbsTiasdlZnJ5MTF0P72ATff/+T9M5a+N+ORuINRn4vdxUBnceHTrud1qFxJmcdqJVKzIbgpBltNNA3PLHgeRWvyaTuXAeREXqMem3Y5D+f7auysZwbp69vkvHxWY53BuMExw81sP096zh6ohnPVV9/QrSJ9ORoElQm0vMTiUuOpnhdFnUH13HkFycXvY/g7uSGFoAkSZ8Bfi7L8i0Vy5IkbQa+IsvyA1f//yKALMv/uNg5wgK4PSwVYwht5qLzOlm1IQ+dXsPowCRNfsXcpLxUoTFPQqChhsmkY7Xay8e+/gGazrbh8/gYLo+idmgIm9tNTKMfWZaxzrqwqn3IcqBOvsmgxewPOiamcwLbD6zK5+ixFox6DZFGHYMTMzz97nX4/TIjtllOHmunb2ya7uFJ/FcHWpyRyMaiDMrTk/jBi1V09k/wUEYGSoVE9ZkO7DYXmq8H3Te5uo0cGwnW4ZHPxbC7ODfwnnUq3rO2lKbhUf6o8DvgPgsoQfcY/9XzLg72tOH2+Wh1BTT7KI2eQlUqSkkiwxxFlE7HkS+fmZucr5WxuK8kj87RSXavy+UH5y/g9fvJ3h8smdH1e4G4xL78POpHRijWxnGxewC720PGM4FqpZIkEZ8aTcbD69j3QBlHjjQwnaRi1DJL7+gUG+oteN1e+tqG2PLoWhRKBYd+dmMBsFRRwaWsg1BEts+t561YAIkEtPNLwI+Ag/Kt8RulAqGtofqBjfMPkiTp48DHAXQY5u8W3EE4Zl1cfCNQ4Cw+NYYtj1WiUCoYHp7Gp1YSEWkgIsoAGnXYeb5IHX5Zxmp18uFP7yElNwmVWsVnNn2JS58LLiBLOR/MSHGZwz3R2pngvlFv4Gt9eWCIh+JzKMlKYmTKyt7KAjzITDmdODwedlfkMTA+TWv/GIXpCZh0Gs409rD/9BWGMqdIjI3AZNDinfFy8lTAx24watkSv5YLk43YfQsDoD6/jN3lxu728NktW2kaHqNuYBiyr1UD9YHzf5lxP8IXKnfywUO/gatzpsPnoSw5kZfbmrk/O4/qwf65yf8aZp2WGYcTvUaN2+/De51GJBLwcFEhZ3v7sLpcaE0qYk0GKmIiSXtQA1Kg9v3YwBSXLnVTV9fHhg05nLjSCcDGogxsVYN0NwSCyoefPU1SVjxl2wq5cqplwf0Edy83FACyLP+1JElfBvYBHwL+Q5KkXwP/Lctyx9s9QFmWnwGegYAF8HbfT7A0BrOe9fvKuXikHqVKQWROPJHRRoxmPdjDC7f5/DI+n49Vq1Ix5SUyPmSh+mgTs+7wSSs0K+ifPvgf/MX3P8G/f+qHJOckcH9+HrIs0zg6CiwdZ7hGotlEYWIccUYDDxUV87UfH6IkJ4mW3lF6kuxcHgm0WIxqUrCvsoB9lYU0945ytLYdj9dHVmIMG0uzOFXTQV56PPeX5lC6Op01ldnMWp3k5Mbza80hjo5eYNrjIEZjZNJtCxvDn+7ZBIDF4eD1pnbYFrJTkUjz5Bj1EyN8vGwDB4cb6bRO4PJ5yYyMYshmJUKjZdaz0JW2OiOZi90D7CzKpvs6vnKtVsWfbd1Cx8QEa1KS2ZWTTXVjH1qVigvdA/S/Fp4q6inLxuXyUlvby+59ufhkmc7BCVxXJ/9rDHePkZAeS2JGHCO94wjuDZYVA5BlWZYkaRgYBrxANPC8JEmvy7L8+Zu89wCQHvJ/2tXXBG8jN9MRLJT1+8r59Lf+iLO/vUhf8wA+pZrpiVlsVgdySJVPWZZRDU7h8/nobx/C0TeEpFCwc2s2mtgYhvsmqD/XjsPmgpCEmJqSeP7oH16E1BgAIr52EYVSQVl+IvHFWcjI9HaP0+IKugl0WhUleYlEmgNuhdrcWc7M9LI3KpcvPnuAjzyygZquQfots9i/0U4+UFiayif+uR6r+2Vc3gHWxivpfOOLxMRGoHaq+F1jK+d6eqEHepJtfPmJ3SQYjYzabPxj7XEaJ6dQKwvIsCSxLb6cmpEhdGoVuz+ZQb/NgjoaVM12dhsSqf7dNLsaP4nX5yfCqKWyLJOjkQFXjEmtIVEZTUZsPEf7Oojoh4zXZDK3mZh6bpipoqBLxR0to4xSUl6awrGhHh6fSGZHt4mk5Ch8T2rQqFVsX5/LD2pqaB0eR6mQUFqh4dnAyugowFmZG/Z5ak8GLLZZSWLV7mIsE7M0/a6H0LRdAGnaSn39MDsfX8tUbR/uqw3tQ+v/Q3j5D8UiCwZh6TTQ5fJ2lp1YKdxQAEiS9Fngj4Bx4IfAX8my7JEkSQG0ATcrAM4D+ZIkZROY+J8G3n+T1xK8Q5w/VMe/f+qHVB+oxT5jD8/mmN84ZtqKUqVAo1WjtEyi1mmYmbBiyEghKSOW3//sg+gMGvoGLAz1TSLLYEs3EdrUSZ8StCp8Pj+SQmLbziI+vT1v7vWLDX00Xuylq2cct9vLSL7ME7nF/Ka1nkezCshMjOYXJ2roHbOwqzSN+CQzLfUDjNtfCxmtj8bGgP6xY2cRJr2GveV5DExMsyczhw+9+gL/sudBXD4vHdPjNFgCzeub6qbYXZRD19gkRcnxfK32AB/M28Trg82sajOQlB5D4eoMEtNNGA1adm0s4PVTTexND0zEPtnP4zlF7G9v5BNlG0jVRpGWHJg4p6btYAz8RLVqJduzUtmenYXX50ejVFKRnkZH+wjnqzuRUw3sWJ/H33/3IGNZAUN5S142Zzp6CW+xsxCFUsHOpzZyYP9FYuIiKKnIoPFs53WPrXqtji1Pb+bYQdFP+F5gOUHgrwI/kmW55zr7imVZvn4PvOXcXJIeBr5NICXiR7Is//1Sx4sg8FvnrVoAAKqcrLlt/9BIcHueNhV6XCjehPB+AHnRBlIyY/G4vdT3TWG1Bn3r13r4mqMMZOzMQ6dV45dlRi/0MjJgISMnnoKyNFzlUUzbnfSMTvFEWTFn6rrYsTaPTL2B5st9rFqTwdjQNFf6xhgdmwFgrPJqHOFqDqfSEUifvK8gh5MHW3G5vWypyKZoUwJHh1qJ05nosI6xeiwdu9PN91+oYmSLn69tv4/9rY1cHh3ig3GrUauU9E9MI6sldGoVkgQHr7TxyOoiXrvSisfnC1uApd+ew55tRYxNWHnovcGFaz1TFlRukGVwe73kSCaGBqZ48blqpqftZDyRx/nWfgxaNXtS0nnjZDM+v4yh3YJWp2b1xhyqj7fQ+0RwcV1Ue7j7La59hu37Sjl9pBFnV0AAFq7NxjU+SUdt99xxoZpzfFosqXmJ1B5rXPCZe/YF44zXrIu5a4RYBKGWwnyW2xNYsHxuOggsy/LfLLHvpif/q+f/DvjdDQ8U3NN0twzR3TKEWqMkf3cJZrMOm82Ny+khVqdGlmVmLHbqGwbmSjfo+gKTeEfzEB3NQwxNxRFvNvKZR7dSf2WAylUZTFhs1Dxfi0ajwjpt58WfnA5rOdkbG/71V1kVKCQJCahv7icqQs/opJXegQlmPE6SDWaezFhN8+AoJoOW737xPbgjZWY9bmL0er6592FifTpcHh9DUzOoNSres76MiVk7791QzosX6tmSl4Hb6yN6nQa3y4tKrSSiOJW4GBNbNuQybLcTYwgkO3zz6Cnsg0Fh+LcbtnLqeAt6g4bE5EjS0+LZUJRBrNnI6cPNbFwXWLimTJphVUUGJw5cwRixeFaOQa9h+75STrx2JSzY3HKpi7LVqaQXptDXMrjgvLH+CRIzYknJSaC/YYFeKLiLEKUgVjBLpncusRw/9LwwrW5eMTh/TXBJR2h3sPl+42vXT8iII+/+9SiVCgwmHVq9GsvwFP3to/S2DdP/dPbcOSnfqwm/xoZiNm3Ko719hPT7M3nfzgqmrHaqZgfomJpgzGbDoNEw/b3ghBXqXweYTYP8+FhcXi+d3ikeyS/kUEcb7+1JQaNRceJoE5u25nN6eox3713NmqI0jC6Z4b5JLpxuQ6GQ0G9J4fkTddhdHmz3W7kvpZD7kotoap3ixbpGsmOiSY00k9AnodGo8Hp8NCfbsTicfGHXDrD5UV4tXuf2+nCM2GloGSQ7I47KtVn86tWLc+NdpdSydlsB+390gvFNsciyjMvjJbFHZsPGHGpre9Hp1OgmrXhcXlxOD65II7IMkVEGKtZl8cqrtUxN2bBandjGZ+fSdj1mNZs35NLeOcrYuHVO4AL4rja22fGezVyo6pwrJw3gK8gIfs6Xw3sz3IoihYKb445bCCYQ6I1a8srT0XvTABjpGePiyRZ8oWUinC7SchPY8lA5I6viudI7xOTswkmhYk0mAwOTvPd9G8koS+Jvf/46A+PTWIr9GNVqVicm896SUs5sUNLcPMTMzPUnlozoSI60drI6L4nGsVE8fj/JqdHU1fZgMGhwOj2sK07n7JVufvrb85RNSnzmb57g0ac34vf5+ZvfncDu8gTqE8VEsjupgBRDJA3yFE+tLuWl+iYuDw6TVB1c4Ty4LeCK+p+Ll/jIqrXUtA/Q2DPCiSudFMkmCvKSKMxL5K/+aT/Oq8HXzJQYtDM+etpGOPdGE/3mFBSShEal5D1R2fzmN9VISGg0SowjFjQ6NRqtGkWEgYSkSIqKUzh9vJn0jBhKSlOJiNBj0mlQKAMBGL864A77g9/fhFqt4sDPzjDQO4HH5cWhdeN2uOmo7ebhD+7g4C9O43F58LgXFoe+ljV2/lAds2Iyv+MQAuAe5GY0qPnHLdUTeNECXfP8umHNYqbtKBQSabkJpKxOQ5IC6wbaG/qxytdy+jVIs9awa7jzk+kEOtvHUDePsrY0laioWEb+Oo0rbYO4PT5K85LJMcfwxGNrmLU6+cR/vRCoAKqBpFMAHjrp5acJ0yRr9aTnx2M267E8YKdpeogxV+Ce7n4zXpMfbZKSvEETJy6NEI9E+mfS+H5XAxUVmeQXZ/K9lotYI11QDoMpMfx4son3lZYB8MVP7aV9eoIDPS3EXNLw0olGHt26iiOHmlhfmI6n1YFJhqn84E8ve7+NqGgj6pYJrNku4qOMKPsCE3H/kAWFQsHjD6wm+oF40iIjGbFaUSgUnPl+B5u35OFOj0E7JQMyep2Smc5Oxo9cP0ibsr0UyePl1V8G1iWEWmOhjXmuae9/t/8vqby/nJn7S/i7v/glGo0KXWoKUTFGoqKNaNUKvvyjj9NR14MhQo/s8+GYdTI9Ns10SSzRiVE89skHkN1ujj1XxXIQWv87hxAAgreVqPgIcssy0MZE4vf76esY4+yhK+EHRUYs61o+n5/6y4G1g/6yaNaXZPD7D65DIUnUHGulq2OUn//oJL71xuuePzg6TVlFPEdPBNIiZ9Y6KI5MoSw6Db8s0+H0YnW72Z2Rw+mfBxZ+xUYaGZ+xUZaZxPaSbF671EJJZgImrYbkiAhyEmLmmsn87fGjFGXGYFCpeSCjAL1SyUunrvDsoUt8YO9afnHkEpWF6VQ39y4YW3ZOPM1Ng0TotXz9V2/wmSe3cfRqDaCH7yvje/9zjIefXsWTJav46cUa/vPsOTIB5nlw1xdlcPkfD173/eeUpmFKNFN9qg2FUoFer8GoiMJg0qE3adFkJ6HVqpBlGUWyASQ48+olIuPMVB1ppHx9NpIk4ZXAbncxOjRNz+sddDf0M9I7RueVvgWTt8FsYLBjmOoDtcv6jAXvLEIA3IPcrNYfSpgVMU+zD7UIQq0BhV6PWqsmpyydmKQoFCYjlvFZGi714tSG+NvnTfgzJTFz2xEd8wKzIaWNr5U1liSJ/LhY8mNjOPR6Pemp0ayvyOB/f3iMRLeLmcnginGfVgq73hv5oxhXJdIwMUr8fjMtzNLCLJIEH03NIrMomZ6aIdp2xWBzu/njzZvI1kYxM+MkRq9H61fSeWk0ULqiMIODVYPcv76QF0/VMnZpAmuOgyfWrAKvxBFHDyM5XvJj4igtTsHbVoeUrMLjVJLwWnDhmC1Nj5RuwjqlxaeHzthZrjhGKdqRimyQaY+bxfRwHI0n+3B2OnjlyGViHaCZdWHwy6SadKRKJuLjI1ifnIRx89UU2atvXWEyUbQ2E4D64RnWP1CKz+fH4XDjruvCYbVhGZrA2Tk+l98fbMXZzYEX5sVb5n3mw91jbP+99Qx1jeEMKQUhO5w4PTInXg6cH2oRhpaGnl8+YrnZaSJW8NYRAkDwlknOSSRzVRpKnQ6P20PnlV5aLnSGLwTSLp6NslwUConSdVmYow10zDrpH5oiJsqI1erkC+/9dxQKifzVGaRW5uFwerjcNjDnM79G98wUj2QX0jARXH2mkCSSY80UrcvCNuNguHeCHTsLWJOaQpROR2PzCD8/cpGtpdk09Y6g0CvZtyafU43dfO6hbfzHC6coyw0IRZfXx4TNzqh1lg73JMVx8bRMjPFacytPlBYzYp3l/+zeznNHD+NyeZEkMOjUxEWb2LO5AIBN6elkx8SQFxtLRVoyFqeDPusMg6MTdPSPs6YoEDPRR88SEWkgvziF5vFpUpOj+fkvz+B5LVzbXvPIek7+tpbOxoEFvZTV9cEVv1JUeHrum+H0yxfZ+dRGjr1cc+ODBXcMQgAI3jRGs5680jQMETq8E6kMd41SfaAGNNobn3wTKBQSZWszidBrqb/Yjdfjo+ypCmKjjPj8fp576QJGwO+XaanpYTJKwqBTs7ogFb1WzejkLE1dw8gyKCUJl8/HJ8o2MG6zIEkSPr+f4UkrOr2Gf/rUT5BlmZh9kSSbI/iLV17lg6nlbC3NZvfqXKZnnaSnRnHwUiuffmQLvzl6GcusA7vTTWlOMhs2Z3O8uROdWs0fbqhgyuEgNzqGTUmBlbU/rr6EWafl//zVI3g8XgYGpmgcm8Sg15CZGsP5/gGuDA9j0mp4vLiIDHMkPdMWvnW+Cs/B8PhIVOM0Or2GwtJUpqftDA5bmLW5CP0UNjywmqEhCwNdY2/LZ3MNv89P7bFGKveWcmGR+IPgzkOkgd5hvJNm7ZL3CunNqzh/hcxVaSRl///27jw66vM89Pj3nZFm0WhHO0IgITaxCbGDWYwB28Rr4qZJmu04bW5u09O6y01vmvaeni6n7U1bp7dJmrpNrm/S1Enj2DHewQazg5BAIJCEWCQhgUAbGm2zSfPeP0ZIvxkYMRaSBpjnc46Ptcz85vlJYt79eQKLhP3OAS6cbGSgZ2DM2sHhqkUZM4gC6G7nyMe+xaNb1lCt2wAAIABJREFUPU0mRdnMTBKT7ZypbKRxlg2rJY4tK+eQkZpI34Cb1/cG1hMGsoOnegDizWaWtFpZWlrA/AWBXTL1753k1P5alm9dRNr2En78XgVtN/r41KYlNKT3UdfRzjPzSuiodTIvP4uK+mY2rCrmo7pLPLygiJzkRHRDHynpgXWGt1+vxO3ykjM9lc9+ZRN//s4eTrVco8/tpV/5mJ+dSWaSg9pDLfz915+ioq6ZX+47xaa0XPbvP0dubipbP7+c1o4e5s3MovLlY5w+epFnn9/I4pWFpCzO5Gtv7aTReWN4QXvUzRrJa9cWQ1I8h45fZGjIT8KFwJTd+m0LOVfdQnv36HSTyRmcryn0UN5NxhTdocaavimePQ13v5ummpbQp4X9e5P0DJNPtoGKj2VaRiJFxdnEW8wMTTPTVNPCkZ0VU/LaJpNiaWkBiYk2aj+sxXljuBBJkZ0tK+cwIzuN1g4nb+4fPWlqiTNTkJFKUVY6iTYrDms8lrg4fGeddHb2cehAPcqkmJuZzKd/fwelm0vIyJ/G8nkzOHDqEkuK89jZeo4VudN541wt66x5/HTPCbp6BhhMUuw718C64gL+YuceZu3v5enfWMvL/2c3c5YVYLPH88iOUlLTHDy3cjG7awKLt9oCtdfa6egbYFPpbBpaO/m1LUu50uGko7qTrKxkWlu7qb7YStWFKywpzmPlxvkkJtlZviEwHfQvFcdodI49J56UbOf6wABDw+U4TWYTGx9dxMkjFwI/O9vU/TM/d/wCDz27io4rXfSHNDbi3iMNwD1msns/4RbiHEUFFC3MJzUjkJmz48IVTv/f0/i8gxHHpGy2sJ9rt6EerCF9BIweJjOZTSydl0vatEQu1LXS7hsipTCT3IUWEhxWHn1kAanTEnEPeDG1u/mzT23GZAqc3PV19dLa0MHlk5c5vzotcBDL46Uj2cOAycvQ8Ei3+OVAgZW8Vyr48p89QeXpy2xaNxd8mnldSZw6cJnUa14uf7qH+qFADYOkajc7/NkMVN/AvLuTz/3tJ/nTf3sHz4IEmns6oAcu/mwPf/zlrRw43cDmglnUNl/HUxPoeQ/Ry77E86wrmcVP367gaGUD/S4vWzfM5dKxehwtHsr8ybz50iGGpllZu2Y2Qw4L7R29mP62nNn7AgfqTMtKMJlNpKYlkJrmIK04nXhLHJtWF7Oz4RL0JuAZHOLJOSUc3V3NQJ+HW8ZFI4u7ATdTbYQaq9avUeiCrXnBXI6d72HzVz/B7m//4ra/Zwj+Ow/9+5LF3akjDUAMm16cQ8G8PFDgM1u4eLaF2soGYOxcLZGIizdjd1ixJ9pISLRii1fYHTYstnj8IW9C8WnJLF5TjCPJzqmaq3S295KYZGNgwIvWGqstnseeWEpKgoWu6z3s+mUFl/s8tF93jpxcjbs8uqjbkjs6jeSy3n6K82pTJ3/6t2+wae1cLjW1k57qoLa+ldY2J2vKimiKd7FtxhyG/H7K4vPp7Ozn+PEGfvvrW/nuLw/gCTn0dK6xnef//BXc6Yp4s5mSGVnkpeYy5PdT39DGtf5+/uEne3n+6dXkZSZz/nIHF1s6mTszkyG3n/wZ6bzz5kkeWrKU0sUz8PmG8Ps1s5fO4vHnHyHOGkeby0/dmRbaWp103+in5dx1bPZ45i3O54XPbCTObEb7NQf/ft/obp4o8A/5qTx0ntWfKOPY2yeiFoe4M2kAYkhSWiLFy/KxJ9rQWnP1wjWOvHYE/22KipjT0zDHmbEnWgP7xIffzO2JNiwmRva++z2ekSLw5pzAGoHWMDg4RF99I66+Xnqvt9ERZ8Xd78HnHRwZGZjMJkrXzyUe+OmL76H9mrwVs0mblojVGkdqugNXWw+FGYm4u/p5q+4ieyrPQzxk7K/DWPlWGzKRBg5FBWRVBjc2bRtHE6MteaiA339mG98/eIw3y8tJffUahcXZLMrP4GzTJaqb24g3m/nSxlLsmQ7WbJ5Lb7+HhG4/SV1+hvwaa3vw9e3tgWI3l2uaOLrBglkpFqzKYtmcaWg0znQPT81bAq2Kg7WNfGFzGQdOXaJ0bj4dmSaeXDIfvH4q9tejlKLmoQyqhrrIcjh4YfNqLtS1cr42kJ/Hb49nxeb5/OJXlbQMekgE3v3wDKas4PyfceWjKbtU6PqLYXrIeChMjZHK2TiaM4dWALseWGx2XofOmakULi6goTqwnhCuUxFp7WCpFTzxpAF4QCUkJ7DmiTLamjvJmZWJIyWBId8QV8634upzY7VbyC3KJrcoe6Qgu5FKdOAf0rj63Az0uXH1uenu6GWgz423e7QGbdCJ4ZDsn4OXGkc+Nv5DTkpL4NO/vY04q5mL1S30OV0sWjWbPucALdednDs7WhZidWkBcxbnc3jXGfZ4bq2Pa5SRncxzz2/gpaoztI1RJ/emv/7EVtIS7Hxt/Sp+XFFFKmBPsPDDf/6AmvUW1s+eya+VLaKu8goFeem4PT7+4YcfUpiZSumiAuLjTMT1D9F+3UlLcxdaa/zJFkwmhdlkIiMvHpMyYVYKp9dNnMnEnJQMNuQUcdneTWF2OkN+zbwZWXR09/GFbctxtw9wubGD9Rvn8b0X32dff2AePdFqIeFgN8cPnx+JPyU1gb5eN1ZLPH39HvYdqQ+UybzjnU+NCycbWPvUCjpauuiN4Pchpp40AA8gk93O6qdX8d///gu896M9HHu7MvBG3uPC1evC1efG4/KGXQ8IFek8rHE3DxC0k0gZpmgefbqUrc+t5OW/eJVdP/ooaG5Yu90j89YLVhZR+sll7D19mf1X2jAbykAak44BxLU5+dLXt/DIs8uZvmw6h6oaSLBbSBs0Y7dbGBz04/H48F5sxT3gwd3v4d3KOj65ejFnm66zPbeIjOWdPPT4Uva/VcXMLgebSmfy7vePsCDLQW5+JpeON7DZauNSaz/nr/TiH9L05lvJnZZMycLAoq3topNrLV001l/BaU7D79cM+TWJP3ejteZfrG0cX1HEibdPkZmTwte++QQ5Dgf29FRmZadxpK+X17oaeLLAxqGCATJO3RydudnVcwnsZrAHDrptmJdDw7lrzCvM5L0zFyEp8PNxVDYH/Wy08XBWSD1mv3W0PGfQiMod/PdgPCOgP0Yq5/IPatjymbXs+dkRwu04HGsNYKw6wrI+cPekAXhAHd91mu//0U849kY5Az331m6Mt364l+vNnRzfdTrsY2YtmM6yTQs4cfoy+4+eD/s4ox+/+D4DvW7+/WQtbV2BHmdiS2CKxmw2YbPF46htxOawYkuw0DvbSlOHk5zUwMJ34fw8GuquYoozsWnbIsqPXMBsNrF0w3wAqg7UkVeYxdaNsxjo91JV2Uibx01NwzVO1LXg9vpIO9xKVl4qBcXZZMwOTEu13eijvf4KPt8QbrePDz88S4LLR1JKAm2t3WTNyyJzmoPrHT0cOt/EU8sWUJw1jezkRDz03OZOIS01gbSMJPp73Bz7qA7KpkX2w59iWmvK3zvFmh3LOCLrAfccOQcQY8azwyKo6hfhe2VjnQMIOhVsCz4wZtyLbna6cCTZ+PIL2zh3qoXdxrxBdQ0jH7q2Lg66huPs6Aija+3oade0E8HTRsppOExls7Jy83xOHKxnaNDPnL/bxHt159lcXEhzeRvt3X08uX4h+c19mMwmfv6ddzGZFEMPL8CRYGXZ4gJmJDoYGvJz8fx1lEnhyQi+N6UgMz2RwlmZ2Cxx5GUG7vXG2atcqrnK0OAQT/y3LUzPS2NwcIh//NfdbN24gLLFBbz29km++70PRq4V3+Mb+fg3f3sLRw+dH8mNZEy97F86NygG457+sU77GqfsQhMAhiv1qEPSetxMFQ3BfxuFi2ZgjjNxoerW+gFT2ZOP1VGEnAMQ941tzyzj0WeX03T++p0ffJfMcWaGBv3kFkyjudvJ8vw8mrudtHf3UZCdyuLZufzHX/6QjU8vBwKnjV1uHy63j937arC2e3AkWilZlI8l3syVrj7qG9uCXiMvK4WOfhcazfmmdry+IZZ1DbHpqTJSMhLpdvbjcnk5caqJjs4+2jp6aGjupP7iNdauGq3hax4I7OzZvG0h06YlUl11mTXr5zA4OITXAe4BD54BLwNZybjdPtxuLz7vEPeChjPNrH6slLZpnfR0ynrAvUIaAHHPyS/MJC4+jrmL8uFXkzdtEG8xMzi8nXP2wukcHhrCYjZzsaOLaXFmXvj0Jr7x/TeZaTGTmZ9OQpKdgd5be4z9fR6OHw0c/kopy2HjymKSE22YzSa6ewa4et3JvsoLQc85sb+R+ctncfn8NWoHPJQumUHduVYs8XF0O9288lo5uz6qweQZ3aEV3+Nj5ZrZXLncxf/+izdGtsCazSYSLrZgtVuwOawkJdvIzErGaovHYjFjmpUyssivHAmE45ufjrvfjWfAi8+RhMflxT3gxePy4upWYefwI3XsvSq2fGYdH/3iKP6hW3eeiaknDUCMGc+Qd6wF4tAF3LDXMEwV6JCTqcbpC5Wbzc//eTdDLi8tDe0kpybQ0x1YwzCmkIjvCd7n7izLGfnY2hO+12vMeT83O4lLda2kZSZhNptYVjlIxcFqioDf+5Md/Ph3XyHtUhuLv7qF7V/aRMPAIHv31mJvGe3Btq8avX9bfByz5qXitZm4ONBDUq2blNQEsrHQc6oT98Bo5Sz1zEJy1hfx3uE60tu8vPK/foXWGmdZDvZUG8XF2az1+bB2eunv89BwsY25z86ns2eAI5db6Ztup7t4dNk2NXfOyMdJvwzJGWEQOk0XNCV0vQ1rghVrggWr3YMtwUpyihVrbiL2jQtGGhHdP5paAo8v6Hr+ojIAvMONx80Gxd3vxj3gpXxnOZueXsbQoJ/yd6vuuD410YfCYmXKJ1LSAIh7Tsc1Jz/4q50AbPi1VVSVX8LZPfEL2enZKdRXt7D9uZV4XF4qDgYaom3PlFFzvIErlwJTOeXll4L+b2RSipKZ2WSlJuL2+jh5uZU+d+CNPqM60CDa7RaWls3ClmCh45qTpvprrFtRTFePi8YrXfhb+4N6185eNydqmjlUeQlHi4ukJBtf+q1N1HV1k5Zsp9/lxWG30M3EHvby+wPbfl19tzbkptzRdR/jbq/QzsHNN1iLzYI1wYLNYcOWYGFaXvrI50VLC9ny2YcAIi4SIyaHNAAxLNLe1Vjb9MJWBwt9rcbRrYmh1zOeHwjdprj/33az6dmVVB1ow2lIKW3aFzw1lGS8l5BF0HD65qWRkrKQbb+5gW/981u4ZyRTVJjJrFUz+acf7YXlgQye1vYB9r0XyHBpIpBQbfqsDGbOyaF7QSq1TdepqQ/8HDL3tZIc8jpDQHl34J4zMpP4/X/6PF39bnD7+eLDpbz8W/+Kf7hEpePVRmwuL0opHG9VYnpkBSs3zuWHLx/A29DJpscX87v/Ywf2lj4+en/0gJdxgTz0/n3Jo1s9rQfOBn0v3KJwaMI+Y/oO4ygidETB8O/W6/bidXvp7bp1vr/mSD3nKy9GVCRmPHWEP87zYp00AOKet+/142x8ZgWnznXg7Oq/8xMikJrmoKtngM88VsaCwmw+t2M5P7myn888t4oXv7v7ts9JS3cwb0Eelj43LY3tHN59hnZv7m0fG05RcRbf/84uFq4r4g9+dztoOLt9CftePXbLY612C+sfLuHA/jo8nkHigeMHAqOUm/+/Hw30DEjP/x4h20BjzEQfrQ/teY3HWL01Y7wbHl/CqQO1ONt7UbNmBD3OWFUs+WzX6DeuB+fBv5luesXKQmo/qCU5NYHnnt/Aay8f5NnfWMPuX1bQUNc6kjTNlmChZEURNoeFG229nKtqClrADN0GadxmalwT6dkyj3mFWWgN9Y1tOOwWNq8OzNuXv17NQH/g9fxWM5s2zWdpaQEXLlzno+/tZtA3vKYRsn3WO2P0nuPbDD3tkHsOOsQVclgv3N9A6O/V2NM3PifietFjXD/S37+kghg/2QYq7nv7Xytnw7MrqT50DuedHz4miyUOr2eQjus9/OBv3mb7s8u5ePYqDXWtmMwm5i0vJC0zCfeAl5qKS0ELuOORkZZIWkoCR6saAeh3eXn7o8B0jKM/OJ9Q2fJZbN22iOpvvz365i/EJJAG4AE0Vu9qMntRkab1Dc0ZpEJTSBiEHjo78PpxNjy7krMnmrjRNnpKNsmws8h4GMm8IHg+3LTvBCazCZ3ISE95wcoictNt7P9VNevWF+HzDXKu3cXZk4F1C09RZtA17GdGcxW5ZwTP+NsMve+bSe+stnhWzshhz+4abv5ETJ7RN3ZjJtMMyxDW7h52fvddjv/D6wwadsmE9sothqyqznUzRz4OXYPAM9p43ZK+OUwP+5bfZZjnhCZ4m+j0ItLrn1z3St4oISJ24PXjLFozh/Ts8dWwnVmSP1KxasacHF548Ys4O3vxenwcevsk5buq6Z6gtQaAdY8v5cD+c3d8XMGcHPLn5bHvvw5Tvb/mnkvhIR48MgJ4AEW1lGQYxp6hMeUABI8IQr8Xrkd5YGclDz1ZRk35RbquO/GfrBm9nmFeOnSFKy4vl+nL53N0VzVkp7L6U2tJzU7hhkdz5aKhULwhNbI9pLKVcW6/b3p80PdszaOjBb8tjtUb5lJx5gpUXSD4kYbrZWdSXJKHPcFC1aF6Vm6aT1y8mbiiWbcUzwl+sdE1gZQT10avF7KTylgEJtI1m9CRU7i1jVt+/xGWehwrDtnBM3VkBCDuWwffPMGClbPHPRKIt8TR0z3A9//6zUnZVTO3JI+2VifOG2OPJhaWzcRkUlRXNALgdnux2sM1F0JMHBkBiPvaobdOsP4Ty/B3O+lsvfN8cWKKnf7hPfdrtyzgo3dOT0r1rIycFBKT7ZwYThERTtnmEjp7XDRdGB19eFw+rHbLGM8SYmJIAyA+tkiH6EGPC1PlCW6d9jGKJKPo4UOXWPWVR6mrvUpHRx+cH92KGJp2org4l3NVTSxYXsjFqiZ8PQMj9QeMVcUwTHm4ijOCrmFrHn1zztwXvO3RWZaDzRpH0fLZHPmvKhjOt28mmLLZWLN9MU3nWmk+Mprx1GS340qLI74wDX/rdQZXLRj5Xnx1Q9A1jFs6jVs9gzKeErxQG7p9NnTLaLiv+yOsETzew4UiOqIyBaSU+rZSqk4pdVop9bpSKnz9OSEicPBAPfMX5JGRkTTm42wOKxZbPI5kG1ebOiYllk2r5vDRsfA1DJRSbHhiGfWnmmi9TQyeAS/mOJmdFZMvWiOA3cA3tdaDSqm/A74J/HGUYhF3EOliXqQHekK3BwYt2oYklAvaBmpMTxBaWL66gWPVDax9vJQ65aHjyvBhsNCDSSXZlK2bzZ5fHLslDYI2pKsw1sS1Xwh+kw5aZA2J49H0adT+pBxHV39QnYOb92GOM7Px02s50tlDT2YSZCYRb4jRlJaKz2JD2e2Y0lKDUl7osdIdRLpdsjb8WkekC8RjLeZKz/7+EpVuhtZ6l9b65sTrUSA/GnGIB8+Rd6uYu2I2mfm3VsjKn5vH9OJsyndX3+aZd29+2SyuNoVPV2GxWXj4M+s5/EYFPb3hp1Q8Lt8tNZqFmAz3whrA88DPox1ELJrK+Vpjrz80gZgxfUBoj3Ks7aNGxucdfuM4a59aAUrR1T964Grp48uJs1nQ8RaUzR/xFstbGObHb24zzZ4xDasZTvf6YDhFg7H3npidzpodyzjwWjk+vwpOyrasZPQeT9bQZzbh71kY+LkY6ipzKrj3Hq4aW+jhqfGk64j0ANbH+duY6NTO4u5N2ghAKfWBUurMbf572vCYbwGDwE/HuM5XlVIVSqkKH55wDxMiyJGdFRQvm0VW/mi+nLj4OJ744gZWPlwyxjM/PnuilXnLC6k6UHfb76dmJrPq0aV89Iuj+Lx33nHkH/KPFKoRYjJN2ghAa711rO8rpb4MPAE8osfISKe1fgl4CQLJ4CYyxlg3lal2w+0IulNM4XoooTuJggxf48jOCtZ9ah3a7aGtuZMPf7KP3hv9HP/gDNrtHjt1xRijA+NuJKUU655ewb43KlE2W9BOHW23k5k/jXnLC9nz473h4zUcYrvJPRDo7BiL5YSmaDYyxivz8iJSUZkCUko9BnwD2KS1lvPuYtIcfecka3YsQynF9csd7Hv9+IRef+2jizn+4VmGBm8tcZg/J4fcoiwOvlHxsa/rdd1d8jkhIhGtNYDvAlZg9/Bi11Gt9deiFIu4C2PtCBlrXtrIuAsoNJ1wuOvfTLQ28jjDXvfQXUXH9tay+tElKLuVq1XBtXmNgs4ZhLm28fOSNXNoqjyPs+HWFMhFSwpw2OM49qtjt8Qeel+h95+QnMCsRQUkJJ/E7Rsd9IaOSiLdjWX8XqSFXsJV+oLxp2iWkci9JyoNgNa6OBqvK2LXsfdPs2r7YvzOHq41tN35CXeQW5hFXLyZ5vpb3/xLVhfj92tOfXRmXNde9Xgp27+0mYbqJvbvPHm3oQoR1r2wC0iIKVG+q5oVawMFYe6mEUhIslFcOpMDt5lOKt1cQk9nL5eqm2/zzAjjHC6VGEnJRCHuhjQAMSzSaYNxpX6AMRd7g54XMsVyt3EYp5FCaw8ce/sEq3aUoZSirdsb9nnGaQ7j1IjJpFi9bSF7f3YQ/5A/KL7Vj5Vy1enjSocXU252UJ2DsaZKjIffzOlpeIADu2shzoq/5+NvxxxzEXiMaaSxfl/Ga0qO/geHnDcXMaf8nRPMmJdHbmHmnR8cYu3jpRx9qzKoLCTAQ8+s5OLpJq40Tk56CSEmg4wAYthEH/6ajO2H4a4R+nVjj91sPGjmvv3ZkfJ3T7Lq6dXg8XL1YqBXbBwt3G4b6KJ1c7l08hIurx65V3NiAhufWcGJvTX03HChXYbUFYZF5dBtq8ZF1vH2qCfiYFW0krdN9KEw2fo6PjICEDGrYtdp8mZnkzc7/P76m24+5mZjAYF6Aps/uZJj75+ip6sv3FOFuGfJCEBMmMnodY21RdQo0l50aE+xYtdpVmxbjFKKZsMhLuPjHCkJFC2aMbKf3+9yYU+0sf7JZez/+WF8Ht/IY4O2Uroj6+WPtzqW8bXMEW65vVd6xhMdx71yX/cbGQGImFexu5qcWZnkz8275Xsms4k1O0o5/GblyNeS0hNZ8+QK9rxyMOjNX4j7jYwARNSN1Xsbq9cfiUhTWVd+UM2yh+agFDSfuzryuIeeW8OhV48w2B84sD4tL52F6+ay95WDgWuEvJ7xgNpYSe7GE3/ozykoVfY4ry9im4wAhBhWufs0WQUZzJgXGAks2VTC+RMNDPQE3vynF+cwp6yQ/a8ejWaYQkwYGQGI+8Z4zi1EOocOgMtF5e7TlG1dQvGqueTPyeX1771PVnEeX/nLX+foOyfZ9+qxMeMI93qhrxWu9z7WfY1lrPscz46b0F1Lsvf/wSQjACFCnPjgNCkZSex4/mFWbl/Cr//BE2x+bg0lqyWDiXiwyAhAiNv44D8P4ezo5fiu05w9Uo/f7+cX33kn2mEJMaHUGKn47znJKl2vVo9EOwxhcK9XeRrvAaGxMl6O557HG8e9/vMV94cP9KuVWusVoV+XKSAhhIhRMgUk7sq93isdb3wTfZgq0kXa8b7WRFzjQSWjqPBkBCCEEDFKRgDivjGVc+/jiWms60e6rfLjxBfpzyPWe8CxeM+RkhGAEELEKBkBiPvGRM+9T4RIrz8ZCdoirdMrPWARjowAhBAiRskIQMSsiVgfuFcKkUiqBjEeMgIQQogYJQ2AEELEKJkCEjFlordEygKruJ/JCEAIIWKUjABETJnsHnusH7oS9xcZAQghRIySEYAQE+hut5LKqEFMJRkBCCFEjJIRgBAfk+wkEg8KGQEIIUSMimoDoJT6Q6WUVkplRDMOIT4Ov8s18p8Q97OoNQBKqRnAduBytGIQQohYFs0RwIvAN4D7pyq9EEI8QKLSACilngauaK1PRfDYryqlKpRSFT48UxCdEELEhknbBaSU+gDIuc23vgX8CYHpnzvSWr8EvASQrNJltCCEEBNk0hoArfXW231dKbUYKAROKaUA8oETSqlVWutrkxWPEEKIYFN+DkBrXQ1k3fxcKdUIrNBad0x1LEIIEcvkHIAQQsSoqJ8E1lrPinYMQggRi2QEIIQQMUoaACGEiFHSAAghRIySBkAIIWKUNABCCBGjpAEQQogYJQ2AEELEKGkAhBAiRkkDIIQQMUoaACGEiFHSAAghRIySBkAIIWKUNABCCBGjpAEQQogYJQ2AEELEKGkAhBAiRkkDIIQQMUoaACGEiFHSAAghRIySBkAIIWKUNABCCBGjpAEQQogYJQ2AEELEKGkAhBAiRimtdbRjiJhSqh1oinYc45ABdEQ7iCkUa/cLcs+x4n6955la68zQL95XDcD9SilVobVeEe04pkqs3S/IPceKB+2eZQpICCFilDQAQggRo6QBmBovRTuAKRZr9wtyz7HigbpnWQMQQogYJSMAIYSIUdIACCFEjJIGYAoppf5QKaWVUhnRjmWyKaW+rZSqU0qdVkq9rpRKjXZMk0Up9ZhS6pxS6oJS6n9GO57JppSaoZTaq5SqUUqdVUr9XrRjmgpKKbNS6qRS6q1oxzJRpAGYIkqpGcB24HK0Y5kiu4FFWuslQD3wzSjHMymUUmbge8DjQAnwWaVUSXSjmnSDwB9qrUuANcDXY+CeAX4PqI12EBNJGoCp8yLwDSAmVt211ru01oPDnx4F8qMZzyRaBVzQWl/SWnuBnwFPRzmmSaW1btVanxj+uJfAm+L06EY1uZRS+cAngH+PdiwTSRqAKaCUehq4orU+Fe1YouR54N1oBzFJpgPNhs9beMDfDI2UUrOAZcCx6EYy6b5DoAPnj3YgEyku2gE8KJRSHwA5t/nWt4A/ITD980AZ65611m8MP+ZbBKYMfjqVsYnJp5RKBH4JvKC17ol2PJOtPfgnAAAB5klEQVRFKfUE0Ka1rlRKbY52PBNJGoAJorXeeruvK6UWA4XAKaUUBKZCTiilVmmtr01hiBMu3D3fpJT6MvAE8Ih+cA+cXAFmGD7PH/7aA00pFU/gzf+nWuvXoh3PJFsPPKWU2gHYgGSl1H9orT8f5bjumhwEm2JKqUZghdb6fswoGDGl1GPAPwKbtNbt0Y5nsiil4ggscj9C4I3/OPA5rfXZqAY2iVSgJ/P/gC6t9QvRjmcqDY8A/khr/US0Y5kIsgYgJst3gSRgt1KqSin1g2gHNBmGF7p/B3ifwGLofz3Ib/7D1gNfALYM/26rhnvH4j4jIwAhhIhRMgIQQogYJQ2AEELEKGkAhBAiRkkDIIQQMUoaACGEiFHSAAghRIySBkAIIWKUNABC3AWl1Mrhmgc2pZRjOD/+omjHJUQk5CCYEHdJKfVXBHLE2IEWrfXfRDkkISIiDYAQd0kpZSGQA8gNrNNaD0U5JCEiIlNAQty9aUAigdxHtijHIkTEZAQgxF1SSu0kUAmsEMjVWv9OlEMSIiJSD0CIu6CU+iLg01r/53B94MNKqS1a6z3Rjk2IO5ERgBBCxChZAxBCiBglDYAQQsQoaQCEECJGSQMghBAxShoAIYSIUdIACCFEjJIGQAghYtT/B5M2NWQTadwFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plot_samples_2D(chain, 100, 200, ax, lims=(-5.5, 5.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the HMC states are indeed quite far away from each other&mdash;the Markov chain jumps relatively large distances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does Metropolis-Hastings do on the same distribution?\n",
    "Let's have a look and run a Metropolis-Hastings sampler with the same initial state and a stepsize chosen such that a similar acceptance rate is achieved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance rate: 0.623\n"
     ]
    }
   ],
   "source": [
    "chain, acceptance_rate = build_MH_chain(np.array([5.0, 1.0]), 2.6, 10000, log_prob)\n",
    "print(\"Acceptance rate: {:.3f}\".format(acceptance_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9Z3gc53Wwfc/uYtF77yAKQbA3sIlFpCRKVLGqbVnujqPYTvLaif06caqTLz2Onbxxj60kbpJldatQJMVeAXYQRCXRewcWwGIL5vuxJHZmgV0uQYAAuee+Ll0aYp6ZeRZlTj9HUVUVQRAEIfAwzPUGBEEQhLlBBIAgCEKAIgJAEAQhQBEBIAiCEKCIABAEQQhQTHO9gZvBrASrIYTP9TYEQRDuKIbo61ZVNdHz63eUAAghnPXKfXO9DUEQhDuKfeorDVN9XVxAgiAIAYoIAEEQhABFBIAgCEKAIgJAEAQhQBEBIAiCEKCIABAEQQhQRAAIgiAEKCIABEEQAhQRAIIgCAGKCABBEIQARQSAIAhCgCICQBAEIUARASAIghCgiAAQBEEIUEQACIIgBCgiAARBEAIUEQCCIAgBiggAQRCEAEUEgCAIQoAiAkAQBCFAEQEgCIIQoIgAEARBCFBEAAiCIAQoIgAEQRACFBEAgiAIAcqcCwBFUYyKopxTFOXtud6LIAhCIDHnAgD4MlAx15sQBEEINOZUACiKkgE8AvxkLvchCIIQiMy1BfDvwNeBcW8LFEV5XlGU04qinLYzdvt2JgiCcJczZwJAUZRHgU5VVc/4Wqeq6o9VVV2rquraIIJv0+4EQRDufubSArgH+JCiKPXAS8AORVF+MYf7EQRBCCjmTACoqvoNVVUzVFXNAZ4F9quq+om52o8gCEKgMdcxAEEQBGGOMM31BgBUVT0IHJzjbQiCIAQUYgEIgiAEKCIABEEQAhQRAIIgCAGKCABBEIQARQSAIAhCgCICQBAEIUARASAIghCgzIs6AEG4WzAlJ00cOzo653AnU2OMiND922mxzNFOhPmAWACCIAgBilgAgjAFyvrlE8fqqYt+XzcftX4tnhq/1iIQayDwEAtAEAQhQBELQAgoPH3gWrQasFbrny9+c899KOFhE8fTtTxE6w9sxAIQBEEIUMQCEAKK6Wi8t1NL9mVtTNqHaO/CLSIWgCAIQoAiFoAA6PPXYf5ns8w2I09vmDgOe/XkLd/P3+/vbFgb2merwyOz/jzhzkEsAEEQhABFBIAgCEKAIi4gAZjskpiPBUK+XBlaprtf7WeeCbePFl8uNV/f6+m0lvAMJGuv83RFaQPJxmWL3Psoq/TrWcKdjVgAgiAIAYpYAMKUzBetX8tsB6ads6gNm3JzdP9WO7unfK4n/n5mfy0FX+dup9Y/Hy3MQEQsAEEQhABFLADhlphtTU57//EluRPHvhq0abVtx9X6aT13utqwt1YT/u5j/N7Vun8bDp7167o7LW1XtP75gVgAgiAIAYpYAHcBc9mszN9nrX60mD/410/wwt++SuXhS4yN2rBb7disdsbHxyfW+WrW5m+DNl/ati/rYDrWjK8Gbb4ylbzhqfH727xuuogvPrARASDcFr74T8+RsTCVz/3V0/zDc02YQ4IIDjUTFGLGYFAm1hlCQrzeY9xq1a2z2xzYrHZsY3asvQPYRm3YrHZUVWXp5kWUvHeekcGbfwkLQqAgAuAu4E7Q3L79pZ/yx9//Hfa9eIzac3W3fD9jRARBZhNBIUGYg4MITk/CnBxEeHAQD354HZsfWoYx6S0++M83ddf5sg602ru/jdYmDVjx837eNG9fGr8ntl3FE8fm90r9usbflhS+1s1EjEWYH4gAEG4L1Wfr+cKGvyQhPY6NH1rLibdO3/I97TYHdpuDEUZRFPPE13/+7++jqiqlh6pu+RmCcDcjQWDhttLd0ktdWSPrdq2atWdYBkd543+OkF2QPGvPEIS7AbEAhJvGW9GRv73s2y0WjCYjax5Yzpm9Fxn8xMaJc7FvlHm9h69WEFpXhCk5ifbKEXIeXE5ESgKjFitTsiBD90+1rnnKz2Jfu1C3LqjH/Wyls1d/D82+fBZd+eu20+5Rsz+AsLMN7n/4WQjmKzDtq5hMd87je63lTktHDXTmzAJQFCVTUZQDiqJcVhSlXFGUL8/VXoTbT0tNG13Nvay4d8msPaP0g0sU71x+44WCEKDMpQXgAL6qqupZRVEigTOKouxVVfXyHO5J8ANvzcX8bSh3/estTf2YoyLYVN/HpaPXCq98BEG9BSIBTEkJE8fX2yw4gCtWyHtqEzVV7QAYyq+69+RR7KUNqppGnRPHQaerdet0wV0PtOdMPlJCvd1jUr9+zR617SkAVC/rfOH53MEnlk0cx+294vU6X5aeN+bLLGXBO3NmAaiq2qaq6tlrx0NABZA+V/sR5oa68mbGRmwUbVh448XToKmhm+TUGIKDxdspCJ4oqqreeNVsb0JRcoDDwFJVVQe9rYtS4tT1yn23a1vCDHKjZmWL1uUz7hyn+szVSeemuocnWs3Zs6lbkNnE+q2FHN1XrvOje2rDalKc+x5azdtHsZcvn7fWogg9UqE7563x3HiYWbfO2OX+c9A2kAO8xgdmoqW0L6Zj9Qlzyz71lTOqqq71/PqcZwEpihIBvAp8ZaqXv6IozyuKclpRlNN2xm7/BoXbQmVJLUHBQeStyJnxe9ttDtqae8nO8y5ABCEQmVMLQFGUIOBt4H1VVb99o/ViAcws/s6KnaR5azVgbUaIh1/e2znP4iHt/ZdtzGegZ5jG6jbXWj99z972q71+20c2cbJmAIfD5d+flMGjsQC055w5Kbp12jjCJCtC+330yDLS4ZHRM/EsH9p764fzdeeSvnt8ynWe+9W20PBZ4OWnpSBa/p3HvLMAFEVRgJ8CFf68/IXAoOxELXHJUWTkz3wO/6l3zlK8ZXZiDYJwJzKXkbF7gE8CZYqinL/2tT9TVfXdOdzTnHM7tSt//cG+1mnbFzt8tC7Wav2emrxO8x62cv5sE2u2FuIMMtOq0aj9bbSmsxo0/nU7YG3p4snHV7DnpRMMe2bmaD6mdk/DGXotP+yU5ufi2eJBm6njI94w+ODSieOQHtvEsblRb5WMFCZOHKce0p9TvGnsvuoANJ9r0lof2U26506jZcZsIJbIrTNnAkBV1aOAcsOFQkBy5nAV67YX4ejoprOpZ8buGx0fycf+aBdGk4FXvvXbGbuvINyJzHkQWBC8UXKggvwV2cSnxkx53mA0kFGQwran1xEW5Z/2Wrq/nO9/49dUn29gw66VuDyRghCYSHL0PGOmTdnZKMbRBY81RVKKR4DRW1qlZ3BUG3D1dO2cKIHNu1bT1tBNcIiZqLhw1/3aOnA6nBTvWsWWpzdAUjyHdrvaSGhTKZUufWJZd2Ey7zW50imTuh1s//2HOXO8hoG+ERh2t4zwTMf09vk93SZOzXUGzed0etwv8v1LE8dDGndQ0Gl9qmeYj9YSynp3lfPgA3kTx5MKurRBex/Fb/jZUXS+IG6fW0cEgDCviIgNJzUnkbjkaBSDASUmCrvNwb2PrWL3r09x+Ww94I4p1J6vZ3hghNKj1d5v6oWhgVH2v3OBdVsWMtA/QvXJ2hn8JIIw/5kXhWD+Immg8wNv6YKeKYa+CquCw8ykZCWQmBmPwWGfiAYNZKXR3t5PT48FVdW3YbjnL57i5OmrWIbHdP3vN35oLSV1U2uD9ni9hq5r5KbR+DPzkshcu4ATh6twOMYZKAj3+vkdIW63UeIpfWBWazlog8dRR/UFbrYit3WgnQI2KdXVVyqpBu3n9JwqNvL0honjyOp+/YVeAtX+poEqmhYcnkyauKYJkPvbumKmi9gCFW9poGIBCLOKyWQkMTWaxKRlhIQFA2AIDcU6bKW9oZtzBy5j6x+YWO85FF3L/iOV7NhaxLFTtTi9rpoeTVc6abc72bJjMZcuNDIw408QhPmHCADhhvgqHorKz6J4xxJOH7hM0MgQKQuSiIyLABUMKUk47U462/qpGoSxzuta+rX/K8GQmYopOHjifq1L9SMhUzVzY4IPXebE0Uru+/B6juWkMzZqB6Cjz05aQjhtDdeyhTR+c/OwvhW0LcttlQR5FGNZE6M4sOcSK9fkEFoYQmmN63xE27hundYC8IVpxC1Emj6pL+JKPu3el/b7O7g5V7dOmyLq2ZRudEvRxHFYVdfEsXO9vgOqNt7gq9DMX3RFdx4xEJ8tsL1o/b5+v0Trn11EAAi3xAMf2cDHv/owu391nP0v7KX+UhOW/mFAX/2rhnuf9XszjDvHOfRqCU9++RG6Wvso3V9Oc207xY+spq2hh7wlaTz3xR0ceLWU7rZ+7EaTa3LYmGt6mCUqFLvdic3m8PqM82fqMcSnsmtNIQfLrgA2r2sF4U5GBIBwQ7wNB4mIDmVo1MH3vvk6pYeqGFaDITMTY+a1hdpWCxoB4NmCQXv/iDa9BqwrOgoPIzkjlkUrs7EFm/n83zyDGh3Bod1l9Bcn0BqSxp8+/yHWFy2g0+nkpz8/Sti4QlBEMGazkbAgE9FBRsxmE0FBRsY+VAiaNNDBHLemn7lnGNOFAf5u+xbqa9pprHNr2F0bolz7djpRFmZgczix2R3YHE6GQ53YnE5MBgObE1M5WlaHZdRGdJ3eitDGIno1GTyeA3G0MYDxJfrvTXCvpjeWxurp2ZGmWxcT6q5+9owPaGMRnkVoWrwWXflZPOaJL99+QkY8H/uTJ3jxn9+gu3nmakCEyYgAEKaF0WSg+P6lHHjrHBOJBDOk5U9FeFQoq3Yup6u1j0Nvnyc0OZr+3uGJ7B+bw0mQycgH52vputrHS6+WMDhkZWTUuy+/v1+/3x7F/ZJedMgVLD32wWUKl2UQHRNGyRHXsxpiXW4kk8FAdI8Rs8mIOciI2WTCHGVifUEW6/MyKV6Qwbiqsqf05jOUApXIuAj+6EfPs3bnShKzEvjbZ/4Nh927tSbcGpIFJAC+6wU8m7wxPMK2J9Zw8v0yhus1YxQ9BpZo0WapaLVfgL4VsRPHccfb9Psas7J6WxGqqnL2dCPjTtdLutNDy41NCicpJoKEqHCqv/bGxNfbPqP3hyeXuj+XYUTv2ula744PRLTYJ45Dj1QQnRDJmvuWUrLnIsNBoRPnqr+YCriSmNakp5MYHkZZewf9Visf7s+gpOQqw8NjNN2nFzZ533WnnGq1cEeoUbeuZ0nQxHHqt47rzk36uVy/X5a+3YM2duCzeZ0PZjvn3mAwsGbnclQVGi8388xXH+WtH+whv3ghrVc7uXKxUfL+bwHJAhJmjLU7FnPpZC1jo7PrGy9alU18dAhnD1UwYrH6dDd0DQ6zJDuZ2dJnBrqHOPibk6x7cAU9VpXa8hYAjIrC+swMIoNDONvayumWlolrDhyo8HY7QcOCZVlkFWVwdu9FhvpcL/nvf+V/AGhrGSSrMI1tT6+j5M1TE/ElYWYQASDcFAuXZdDTPkBP+8CNF0+T9JwE8pekU3GugYqj/r9EQ4KC6B3yT6OdDuPjKiffO8+CzUvZ9uhKRofHyMlLoKS5mX6rl8Hzglci4yJY88By6subOfTyca/rGqtaaa5tZ+XGApyOcc7uu+h1rXBziAAQgMkm/lS995MzYgmLCOH8Sf982tq0P23TKU+vfCwu11FkdCgrn91Aa0sf+0vqXPvQaP3azpieJFyyEplo47nHivneyss01XRMua76E+77LfqGvjgrIiNy4tiS7na9hL0dNXEcbDCT9L+xGIwGElNjeHekme4QC1zz8MQkD02svZqscXs16Pdx+e+yJ45TDrn/DPsK9SmmqcfdrihtQRfAwAL3dzW+3L2u/in9s4p63C4m1SP1tU8zEzjqFycmjsOfv5/nHivmV78tpavHQtS5dvc9OrvJLEzl+X94lh985QWaKtxWjz8dOrXuntN7LrB25wo6G7sZGZwsvK/fwwmU7j5PTGIU25+9h8qSWtquTv0zFvxHmsEJfhEeGULBskzOn5j5dglBQUY2bFtEweI0jhyooKbSFQcICw9m2yMrCIvwL7i888FlLF+RxfPffNrnuuWJyXxi6QrCokJ9rtMSbgxlU/xylkfnc/54DQfePMs7vzpBYWwCxcn+VeveSXzlM9t5ZtdKPvXE+inPf/FfPs7a+5fxhW996qbuu2BZFlue2UDlqVpOv3+ep778CF/89mdYt2ulX9f3dw1y4KVjxCRGsfnJdZhDvPdsupsIiwrj3o9u8rvpob+IBRBg+Ftar9XeDKMjFH9sM/t/8Dbj4/p0Rt39PAp9ej+/aeJY2zJB24Jh2bpcIhUnZ35zBOvwGEMPLkVRYElBGg9uKWL7hoWM/Wgf+w9XEnbWrUY3POxRMNVn5t9/dZg/NBn48c+OYV6XzwP3L+XFpmosVnesIv0gfPi5Ah5btYw9347gUlkzzS19NDZ04/iKO+VwU/y1IjBTFOrBhxmx2nivogmncxzTYnf7g/rjLWTGRvNM+hI+qLxC+G53EDdGUwgW2qxPsRzVWBsND7sDF0q0fuxpf4db+CVc0ruZLKnuc9pCsKI/0jeU0wZ+RzTFY6BvHKfNtTny3T2MPdRBzaHLPGQOpmtpBrXlLaiqyuiWIg5XtWJIiOLQyyfwhvZ3IzI2nBUrUie5eyx9Fl743gHONA67kwh8zDe+TsWpGswhZtbtWklPax819XqX5N1WQPbRrz/Ozk/fC8DBX3t3l90sIgCEG7L1mQ0ce6Nk0sv/VsgqSCanIIVLp+voueRyxYSEB7N01QLMQUbKa9t58e3T9HUOceq090HxWhpaevnaP75GSI+NP/3TR1mzdgFNb9nYfVbvsnp970Vsdidv/vwkw8NjpKfHsnFjAfbEYfptw1QNtRAbFE92eB4WxyAvnall3Ed0ualvgI4hC/cvyqO5q53m9n6va+8UznxwiZHBUc4dvAxA0op81t+3GIOiYMuIYsWSTP7uW78lu7GTpKwEOhu7p7yPwaCw+t4iVBUOv3x40u9QYlYib/7PqWnt0Wa1cfT1ElJzk7n3ybWcO1zJQM/dmSlUV9bAj772v5S8d/7Gi28CEQBzwGy0aNaibROsnQc7aR9+NOfa8dxmVj9STGVVN/Y214vNW+qg5+eKqXVrrNeLv2ISIllanEbD5Wb2v+/6w0/MiKewOA/riI3L//H+RHZRSEIkF0YVHM1DBAOdj7nbKSSf0r+Qw5vd2SHrPrGaS53dHHqxkdKLdYSMutc6QhVaLENUtnbRkqoAIfQzSnlnE2O/glXZafyfdVvpax+ltLaJzj4Lv/+7S/nvQ2foGLjmj9a0goivHCcoyITZZKCqsZE1S7NYlJuMCuzpa5iwPox5+tRMS7Z7T3mvurX+bs9WGP/j/vlVfEevvUdWuY91LS46p34Zg0fxGPqfpfbnZx0eIyTc3aKj80ItxsEBcpdm8aVnHyIpO4H+c1f40eul3PvMeo69dQa7zaErVsse7id7cQZnd5+eyO7RYjAYCI0ImVwYeJN/D21XO2j74busum8ZQdnRnN5z4aauvxPoauqhsbJlyjjJrSACQPDJyu1LeeCjG0CFf//qL6d9n+DQIFZvK2JkyMqhl4+hqioL1+aRlJVAd3MPR18vAfQvIYNBwTl+c3mdazfkER8TxusHyhjz0e7BOT6OUXP//PQEMvLiaOju4y9/sxejFUKCTHz7dx5lU1EO6/OyeKWkjJa+QZzu+DDhoU7s9mv/OZzUN/ewakkmH398LaPvHea9C3duEZg6rhIZF05RcR5B5iBa6zrpbO7hXz77XbY+s4GXr01UO/HOOTY8soojr7s6tEZGhbJmXS51bx7zmd2TlJ1Ae93MBXLPfVBGREw4W5/ZwNWLjTRWNN/4ojuEjoYuUnKSGOqdYWVRCsEEmDrrByAhPZYP/8mTXDpTh83q4PThKmw1U7tkJjUWCw9DURSWr8/FuTaVU2UuH/49LYOEhAdTfeYqnY3dOktE2045Lj6CuPZ2as66nqftFOrZGM2+diHLlmdiMhroWBXEqQbXH3/uy3rt8upHXJpycmQEW7d3oigGEsxxNI20cPLvoyfWWWNd+RHJsRH830e38MGhClrb+0lNjmYs3kRrzyC1Ld303z+ku3/m902Ehwezbl0upw9VMzLs0rjrH9V/fxPXubNqeo+kTByHteuW6awoz2IyZ4jGinjFbQGNxQXr1mm1fs9zoUfcabbXYwVGk4H8p9axeHE6YzYH77x9HuViAwaDwqaHV3Lk18fwJC0vmdikKMISolFVOHeogvFQfZDdsz302p0rCI0M5cirJyfdzxthUWGs27WSkvfO+9SG81bmkJ6fQunu84xa7vwUXUVRWPfwKk69433uti+kEEyYFt0tffzw794CXJlAmx9cRo3BSlNV6w2vXbAolczcJC6crGUwVmHjihzsdiflv3nNrz9KxaD4HXdYuDCFMaudhPRY3mqo831fIDMmml2py/l10xtcsVxfHz1prd3h5PS5egaHRomNDuPoyVos6UGkxEWyaekC7GljDNitlPW2Mup0pWIOD49x4EAFxtE7q4VBRn4yWQtTcNqdlDX1UlbWzKZNBdjtTszAyi2FnD9aNeW15uAgdn5iC6/95CCN1e1TrvEkKCQI+5j9xgs1rH94NX/4vd/hg18c5u0f7aWpqnWiOlzLlfP1NJQ3U/zQSoYHRrh4+PJNPcdf/BVIt4qqqhgMM5+0KQIggNG2EvDUznRc88kOAPsvVLLosXu49741lB6sxGp0Z71c10fjkyJZujqHioEhrqgOFj6xkqH95Vz82SlUVcW+drHu9jZN+wOTpnePGmrCmJ6G6dpUxxHNOocmmyUnK57hZeGMORw0WvuIvTS1VWswGthpyiQiNJgLpS38MjSYD6qTiA7JJCkygvyPhhNs0rdiWJGSwrMrlvNyxSW6hof5/U8sYU9zNSowqDpYn9yI2RDMRwtycYzbKR2Jo6S6ia7BYWzaUcaFetO9q8St9WcedQtDS4ZeQ7fGuy2i2Cr954q90DdxrG2nMZyiryWIL3cf9+cF6c5FpBSzpCAVk9FAz9sXOf6uy3+uJsURBBiHxgjqGSE0KRZTZDgWu74pXVRrB2t2LKGhopXv/P4LbP/jx2lwKow7Vd/xp4gIDCGhqDeZWHDq3bOoXxpnsGeITY8XMzI0OhGA7mzoprGiGbvNgTEiAhUoOVBJbLSZ7R/bTNnhy3S3eG94Nx2e+vLDPPnlh9n9wgEuHiqnubqNtisdM5owMZuIABBumtryFhpqOii+dxHdg1aqylzulpBQM2s25TNisTI0MEJ2ZjxX6jo5dKya0JKbrx9wWQC+XZSpKdEkJkSyt7GV+5cW8P7FajxneQWHBrH9oeUkp8dS5RzBMjrGktwUkuLi+PyGOM63tNExZKG0pYUxh15rP9XUTOVANwcb6xiy2UgIDWPjggzeqavEqaqEJ1zApJgINYYTZgqna2A9H9m8nIVpiazIT2Xv5VoOVdfhTHW95K1OO4N2K46kYBYmJrC/9ipwe10UQSYjKxekEhUagrVxmFMX6nE4xomaQnO3DI4SERXKyuUZnNrv1qINBoU163KhIYwjr5dO/JxOHK1m4+ZCjh3yPfErKNiEfcyOKcjoc50nI4MjE2mQZ/eVsXBNLikLkrh6sRGb1cbKHUsxmU0YQkLpaeujsbKV7pZeDrx4lGVbiihYk0vpe+dnpMHcpseLKXn3LM3VrZS8dx7rsJWMglSKd63EYHRp6z0tvTRcbp71tinTRWIAwpRo/fKKx1AVrbWQuSKPwjULGHeOExwdgc1qxzoyxqVTVxh2ajRRH8PNtYPJQ5vdPvXktBhMqVHU17s0vP58tw88osVOXGw4SxalceREDU98J4qusW7MBjNjr7qyhYKDTWz/yDI6LRZev1zB1b4+Pnnvb1BxaWevf+kBNm0q4PjxGgCuPO3WttUQtyUS7JGbn7Enlq1r8thfUkP6717SnTtxphCAdz/ySYoSEmm42kV1ZRuDuUZK2lrYfbWGEKOJP1+2lU0Ls/mX3x5kzyl3PGOgUK85Rle5zX5tIzvQj5zUjnrUWgMAEc2u/efmJpKyIAG73cmlihaGLFbd91sbf7neKC8k1MzWB5fScu4ql0tdsZjs7DiyFqVx9kA5A02Tg7g5S7MIMpu4Wu09Fz8zJw7rsJWcJZmUNuj7+/g7LlJL7vJsMham0nC5mYbLLoUkPi2OnCUZmENdn2uga5D2uk6Wbl5ES207V87X3/RzAIJDzWx+aj0l7527YVA2IT2O7MXuPQwPjFB/qYnBniGf103Fkk2FNFa2TCsQLDEAYVYIjwolqzCNrEVpNNR28qtvv4vTce0lNs1e8ddRFAWtehIZFkxidASJMeFkFIZQvHoB+49U8tD9S8mPMNI62kantZuaa2nlO7YuYkVKCn+zfz+nml0vhesvf/czbn5fw1Yb+05Vc/+6hTQZrzDiHNWd35Gdy18d/oA/SSvmh/9vD031PVR/Iox1qen8yYYt5MbE4hwe50hlHWfrbxxLuRXio8NZk5WFwahQd7WLw8dvLivJOmpjzaZ89v54HxExYazeuoiGC/Ucfq3E6zX1lxrZ+NhaurqHGeqdunlbcnYCJe+eI2dJ5pTnb5arFxu4erGB7MUZbH1mw8QLvqfV7fKJSYxiwbIs7DYHxQ+tZPOT62mv7+TIq6f89t/Hp8WxfNtiDv/mBHYfWWbX6W7p1bmdwqPDyFmSSVSCqxDQYXPQVNVKe92NC9fa62c+E0gEgDAtkrISeOT5B7COOfngpeM017QTm53MlsdWc7n0Kp3NN+9rzc6M44//6mnefrkE25iDtKx4DNEhpKW6nOlt0Spd/RYqGztJS8/kH779DuPjKtu3LOLFxjdR0Vuzp05f5Z8OH2HY5j3Q2NU1RGJiJF1d/mlkEaZgPrJzFe8cKef9E5V87L6lVA5W02d3VaKuTk6laWiAmt4e/vxbL+quLWlroX/MSoQ5mJVKEn+wcyMnahrY0zGzqaLBZhPLC9IIDzXT0z/MycO1blda6M25XAoWp1FzuZXi+5YwPq5y9O1zOIdu3JHz5OwqUVYAACAASURBVNtnuP9z93Pg5ZNM5WUwGAwoBgWnY2ZnL1+3ANLzU9j6zAY6GrqpKnW5H/u7Bjm332WxnXjrNA99bgef/buPMTYy5ld1bd7KHOJSYjjw4tFp7294YITy4+5AepDZREZhGhseXYNiUEB1pXw2VbZMEjC9bX3kr8qZ9rOnQgTAXYg/DbluhDczPCo+ks/87I+oq2zjjdcuMFjdOHGur6GDAw0dLNlQwIK8Qk69fmLC16rdk+fcgOuByX/4wn2sKMrBEBHMl//0RVLHrKhhJhqvaXHWWAMmo4EHiwt5paYKR5GJvOR4LioD9Px/7qla1xu5WVA58MtL3Ls2n+QyV87/i8vd7qbmjxupU7rZlpFLbdMIwe1ud0veFvfnunIkZ+L4d7Zu5OMPraQ+cZjf1lbx4/+ysK1oM4NdvahAdLaJC63tGFEYKHBHI9QQBylhEaTFRXCg+SpN7b0MOKwcrq+nf4VbQJn69X+Sg+vc7reES/pz11NVAYbvjaMoI5nU2EiMjVYulTYzcs3v3PCkO/Cb96reF60twnJq5hIMZ4QRFR7Mk8/voKG1F8sb7ZRdSxn1ljKsRVVVjlrGWPFHj3Ci1NVuwvxe6cR5Q2goESnxjDmVSS7GmaC93UL77kskZ8Wz9ZkNdLf0cvmEXtAefuUk1mGrX9W1ax5YzlCvhdLdM1uJa7c5qCtrpK7M/fuWnJ04EcsAGOyxUH+pkeGBkRnPBBIBINwUq+9fxqYHllBeepWh/qnN5vKTNYRGhLDp8bU0VbXqfrl98a+vH+JPR8Y5e6GRYLMJg6Lg0GiPEaFmvvT4Pfz03VM4DOMoCixMiee9C9XE+7jvucpmVhZmcKaiadI5p6pi9NMPlBQRTnVXN7+8dJ6Dje5U00MVV3lgWQGL0hL5l/OTc+QBIoLMrE3O4O06l2C1jNl499Kta/4pcZEUZSXhNENFUyeXmzqIuXJzqZXe2LGpkFVLMjh4qobIuBu/9D0ZGBylt2+YnKx46hvdfZbCo8MYHhghNCIU6/CYjzvcOh2NPbRebiAhPY6tz2ygr2OAS0crUVVVF1D2htFkZMvT6yk/Xk1Xk/cK6xndc0MXHQ3u3k6RcRHkrcghPCaMtLxk/v6db/DDr/1M14V1uogAuAuZzclJLm3pBUreO49jcMRrczmLxcKxqgQWFCzg3p3rOHG5mdFrGqlnfnzqoWsaPr38XkcfJoOBh54ooqqlC/uQHavVpfX858cfYtnidMZ7xzjV3s6y3BROXW4k2RGKYrVODINJ+oT7D+PKkRz6GGNhQQgDqgFFU3T12Y8cBCArrI+CjCZ+Wl44ca510N0CeizF7pr4lZdGSXsTA44rROW0EgW0W5MxG40o8Ubeb7/CuowMSq7FG97+p3+7docg/r7mu7x1xV105dRkeypWt1smokEvjEat7sC3I9T1Yg8NNbO0KA3b4lDaBoZ4t+0KQX14peAPvRda2TTFddqU0/DmEY69fA570xAXT9Sywt27DiXJ3QzP06Gk/d0Lq+qiqaqLex5cxlDzINedbFmL0mmsaiU6PpKRoVHfKcjTxPNvoLull8OvnCQ2OYZ7nlyHpW+Yi4cvT1lDcJ2ImHA2PLqGo6+XYJ0FK8VfhnotE3UMf//ONyh+aBUAf/7IP97yvUUACDeFP1qTlrqaDprqulm1aylDQ6OUld24PN8xPs4HF2r4nQfW8fZhVxL7thW5vP72WS6UNfLb3RcYSjERFhxEe+8QGYnRpK+Lw2BwvTxNcS57QAWiM6LpGR6hqrObtZnpnKmZ/Py20Vqywl21CSkhkXy+cBO/aDhEh3VwYs32zDwONl8lMzKGPpt+UPmDBQW8X1OD1eEgLy6OXQUFGA0GFCUSVR0iOPQRdtfX4LyFjDujorC4MJX4+AhGR2yUlbfQOI0A9s0wMjzGwX2u7/+Q3U5kbMSUPX1uxIl95Wx7ZAV7y1y+75jkaKrL20jJSWSg/uazYW6Fvo5+jr52iqj4SDZ9aC2jFivnD5RPikWk56ewYFkWH/zyyJQxjLnih1/7me7/t4qkgc4Snm0RZqI9rb+tnH3hb3xAu85zjqwW7T4m+YYXuH3KSmcvSRlxLC7O5cTAMH0a95G2zXPvA25f/qrCDNZEx3G1poO21j7qRkbYtD6f46dq2XbPQo6dqsVxLeNIOwTls+vdbphRZwqRpkgigqLZGL+N137Si6qqhASbSMyP46Uj57FYbTy0eiGt4aN8fu0aNmdn87Oy8/zN0QMA3Fc8SlRQOJVD9ayNLeKDfxud6CG0pSiHqqsd9F77PF0rDXyqeCVfvGc9v2z9JbZxGzVDVzEb9Dn277/nzsjTFnh1r9B/C9MzI/hk0UoudrfT/Z2r9GiC1a1bNHEVjYIa0abXai2pbr+xZ0tp3XzmOu/COSQ8mKLiPFd3UO3P1UeKsJb41FhylmZyZu9FNn5oLSfeOs2q+5Zx+XjVtHPkZyLWFR4dxsrtS7FZbZzffwm7zcGSTYUYjIaJmMfdgLc0UBkII9w2Opt7Ofj6aXJzElm/ZsENUzAHh60YTQZyC5JpbXb5OBQF0lJi6Owemnj5+2LEaaFjrI0rlkreavk1QxYrx85cBRQ+sX0VX35sM89tW8V9KwroGh7mz/bs5efnz/Pjc66ApdlgpCAyg8qhegCMimHi5b80L5X2rsGJl/91Xrt4mbfKKhh1jtIy2sagY/pa7vqUDD68cBmxwaG6l//txjo8RnRCJNueXkdYRPCNL/Cgp60P64iNtDy3Cy441DznBVLDAyMce6OES0creejzO/jhuX8lPDrsrnr5+2JOXUCKojwE/AcuV+JPVFX9p7ncz0wyGwMptPf0u42DB1pNSZeN40P788zn1z5Puw/Vow2xNsPEpskwKfvHV4iMi2DnzhVUn75Ca477pXC9aAkgLMjG+IpgSs7Xk39/AXUvniEoLYHVsZEMf+Yl3FfBqScWThz/6rP3ThyPpegDos8lhZNqjSUuOYLzlgaCc8bZN1bKbzttfH7bCPXDl6jnHf7pXtcLICFsJ898kI3V6XpaujOW9m39ZEfGYAmHkr42cLvRSdoTDIzRndhPgymbA1VOIJ6VD+tfKEZN7HMwxy0JHTH6/R4fOcs/llsYMo5S+LGlHL9UP3Eu7YhmaM+I+0XqWQgW0u+2MDzjL9psL1/twU3JSeQX53PvE2vpf//YxIyFpO/67w6svNDMk7+/k+qzvns1+ctMxrpGLVY2PLKG3OXZPP4HD1Hy3rkZu/d8Zs4EgKIoRuB7wANAM1CqKMpbqqrOTtcmYV4x1GvhwEvHWLQun3u2LaLkeA12u94Pm5OTwOMPLKextY+2zgFyi9KIiY+gu32AIC/39UZ8SBjrUzLZFbuYhzfD9359hP+O+gAFhYfTl7O7tYwQg76JRIS5iFF7PVZnuu7r0eYQiuIS2d1Qg2EKZTglPpKWrgEMabf+52VxjLG7tQyAhPp0dq1bRGlVE90DN87Fn2mqzzVQe7GJo5HTb3rW1dLLfc9uovbU1E3l5pKZ9q/fCdwwBqAoyh8Cv1BV1UeewTQerCgbgW+qqvrgtX9/A0BVVa+h7TspBjBf8ebb97RYfFkY3mIRnjEAbbYIXobIAHR9bhGbC3No7RukY487gyc8zMym5TmUlFxleHiMxR8tIjclnrqOXmr+RV+MU/9c2sSxeZ3L8gg1BrF1fwE5CxJJSozi4sVG3kxqZ0tONu9V1zCe7FLDw0xB3JuWh6FC4XhFA47xccJ2WdiQmMf+9gpMb7mHrexYnY89z8ib1RWogBKtd2GkJ46zISGPvW3lrI9YwYGWKzjUcaJK9K2cz/3ZDyaO8/Z9duI4Kkb/fXIedj87+bTL3756dQ7jqsprq9yxk4Lvea9K1baAtsbrhZJDM9wmosVtfWhbRABQ18yyewq5crGBoTa3pefzZ+5BYmoMT3xmM9s/tJp3fnaY2guNnNp76a4b3zgfuZUYQDIu7fxlRVEeUpTpFM9PSTqgTcxuvvY1HYqiPK8oymlFUU7bmd2cYWFusNod7LtUy6jNzo6ti4i4NolqeMTGgQMVDF/LFTcaDLT3DZEY5dnuzY3RoLAiNpOdqUv5nfxtJCRGUnaxiR/84AOOHavBYrNxvq2dlampE9eMOOxc7GkjMjSYwoxEALYmF3K4Q6+lKgrsWJXP7qs1TKU2mQ1GtiYVsq/NZcSWdjazJmnmB8afPVtPc1MPj+csIcbs/2D7W2VkaHTaQ8lTs+NZsCiVX/7nPn70d2/SVNPB+p3LCIsMufHFwqxxQwGgqupfAAXAT4HPADWKovyDoih5Pi+cIVRV/bGqqmtVVV0bxM0Hn4Q7h7quPg4drWLp4nRWLc/SnUtIiGBw1Mrus1Wszk8nIkr/4stPjOeBonyeXLmEhVEp2MYd/KruBO/vLqPao8tl29AQqZF6zbV5eIDmngHWL8xiVW4aZf3NOFR9kPmRDYv54GwtVsfU2vaurEXsayufaEnRbxsldpZe0N3dFt6qL2dlQhprEifpTbPC8OAo4dN4YWfmJZGaFU/JgQpGLFYOvXMBBdj44AqKdyyZ+Y0KfuOXk1JVVVVRlHagHXAAscAriqLsVVX169N8dgug7QSVce1rwm3Cl+mtdft4mvlqktstof0FGlmdrVundSM4vQSOQR/MdIaaOH24hoSECB7ctIgDYx109ltYuzqXg/urCAF+/sopnv70Jt7/1XFyitIxmgwMro9m1DnGlbEeSi62T+Tcx+XrX1g577jcKH1nmthWkEWFZiTh+j/7CU9k/hvbN5v55C/NWGyuYKpxp4XlcWkMhg9zsbeF6PfcetNAoat46sHcfA5WNjESYwDcBVWNln4yI2Jo36oPkC/66Rcnjjfc5z3jJOFZ9/S1y6eX6s4pHcEc7mgmMzqa4j9YzeHKeobHbAwV6gPJ+b90x1ZMo/rSLW3Q3exjhKLDYmGwFTIL0/xOEXYmRpGbn0xkVAilZxtQNb8DJ17qxDkwMONDzoWb44YCQFGULwOfArqBnwD/V1VVu6IoBqAGmK4AKAUKFEVZgOvF/yzw3DTvJdxldHdb+OCDcnKeyOeBFa72zuGhZoZHbXT3DxOzPJJPf+Nxjvz2DKqqcqqrnrZRV+GWU9X/WgeZjESHhxAdHkJ6mIHway6me3Yu4YMSdzuGrPBiBuwtLI55hF0LD/CbS64iqMzwGEKNQfSODdNltaCgFyqbMrK43N1F/5hV8+p3cbmvgwczF9LeNXttBJoGBuirHWLrogV0Dg5zktmZhTsyOEp4tP8uoIJFqQSZjVzQ1Hm473VzBYXC7OCPBRAHPKWqqu6nqKrquKIoj073waqqOhRF+QPgfVxpoC+oqlp+g8uEW2Q6qXOeWp5DkyKondMb5vGHrktb1RbGeQSEh1e5EzrDm/Xn6t+oJXeHmed/bweGXidZaXEMj45xz33LOfzWWbILUynZV87C1gKWmly/zkaNAmx5xoFt3En/mJXWsVEuW0cZsrsCt8crh4hIMFPe6drnY85XMTgAtvLsyuXsbqplXFVZaM7g3dpq7s/KZ6DXiNYL8tjGcVTq6bM0UwTYvun+LM//+DUAUsI2cfxkvs51lHrcvcnzY+7pZrZofXTBEeO+JsfDm6SdCTwaqfJ+y1Vy42P5nHM5qqpy+OwVhkdtGEfdWr4lVS+8oo5OLSy0absA/Y+5hPDI8lxocRd/TYqFXEsnLlqXj/3yFcpLb34Q0EwwE0WTgcANBYCqqn/t49wtVUuoqvou8O6t3EO4+zl16gqqqrJ0fQ4fun855ytaeP9Xx7l4oobxcZWFK7P537oWxq6V82urYvuN3hujNfT3s2thwYQAAOgefp9x1cpfHFzNF9esIyU8gq8dc/2KembMZcZFE2Ea5Xy/76Zu3aPnKM7axpGrkzXhmcYyZqNocTL3ry/E4Rhn7xykWy7fXIhlYISaOXr5C/4jvYDuUGaiDN6XL1erNfnSoAwHz04cqz7aBKsard8zVVCr9Ru7BnXnSIzCOurg8LtlvLukmb5qO/9VdYrQV93xhYwgE0u2pVJa79I+Ldnu4O2/bPuN7nZfP/ThiePQDoX68R5WmhKpauvm63/1hYlzjpBuYpOCeHf/JXblFNLUP4C530RwexB/+fyvMBuiSQxdze/+IBXImbguKsMtJP5f3Y6J400DQUT68f5POelhAYS6/0Q71+gT8IxW99qFliiW5CTTPTDMS998hyuba7h8tJpYyxhDC93DiU0ePc20rTeifnFi4tjTlZV0TdWLseotCHu8/vdm/YIl9LT20VA5edCNsn75xLGvecEzgWj9/iGtIIQ7hraRIf729D7aRvQ56s1NvQyP2ViSlnzT96zv6iM7IXbS1zctzeGF90qICg/hWF0jIzY7jywuJDEiHKMSTGrYRlqGD/j9nJrmbvLTvefIT5cFsbE8vHAhsRGh7DtTw/naVkYsYxzaXcaI5famTW/YkE97Q/eUL39hfiLN4IRbwrPpnRanpsWDodydzTL0oD6bJeqcO01zpDBRd06bSdS13p19lLRf/5K58k9RFCdk0TM2TGWte0SjuV+v49hi3NZBdLZrjm5uRCIKCh32a/7rqAws78dQ09yNQVG4f/si3j1fydZFCxix2fl0Si7/8dVfMmKx0vTJfN39R5Pdf0/mAbfGHlWvsn1tAQdOu+YP9xW6z2nbQowm61NPs9/VzCbudS9ckJdE0EfTqOvro7Krm5gq9/1GUnS3IHOvdwtR20JC27pDm+kF0LbN9e9tS3Op/tob7muuWY6bdq2gtqyJjip3TMGXZTodH71nNpqurYmPc4I0gxMCgNLuRtYnZvOrxz9MfkzcjS+4xlVLF7mRLsGTFhpHmCmYmmZX1s64qnLw8lUeXVWEZcxGdGgI//XNV1m1bRGLi3Nvan9DI1Yiw26tlmVhUSpbdhThcDh5r7qGylnMLvKXLY+tpuJM3bTGgApzi1gAwi2hy+n3yO7xzCS5jrnR40Whvc4jFtH6cBpT4enLvs7P/+w5irKTOFpVzxdfeIPoB/RFYM2N7tlh2kEshXEJhMVCWlgU7zbqG6FFlYRwT1EOT21cyrdeP0RXm8sqSYmLpGhNGucaW2kdcH3N6eX9nnrcTpDJSPHqHI6XXKFnibub0fZnT08c11r0bqLen7hqK5YXpBGdHsHlxg5ael1xEu3kr7Aq9wQpW5Ze+Dk0c4C1zwWPhnIaK82bBr3h0TWcfPuMa73BwLaPbuLs3osMdA9OuV6LVku3r3U37/P8fZiNATGBjlgAwl1NeIiZXesX8a1fH+BoVT3/+vahm7p+cMzK/1l2D+d7pvZf56bEUd3ajcno/pNp7x3inbIqkqIieGBxPqFBvnMq7A4nRpPhhm2wr2PAwJqiTHYUF9DWPcje8zUTL/+5Iio+ki/9+2dIykpg+8fuoeTds369/IX5iWQBCXc8OSlxLEiNY09pFc5xlS++8MaNL/Lg+VXFbE3N5fNF6/jb0/smvq4A9y7N5dHiIl47XsbKBWkMdo8ybHX7zs83tRFkNLI5PxuLaudkw+TZw9epqGpjcWEah+nyusakGFkek0+IIZijzd30D12LacTM8ggwP9j2kY2s2bmCzEXp/O0z/8aoZe5GJQq3jriA5oD5ErDy1cvf19QnrcvGV3qntw6gnkG/wU9snDgO6dH32dEWl3U+5g64Jp5yuQ3WbMrHMmTlaJj7WTG17v3WP6L3yWTsd99f6w5Jjo7gE0+s4kcXS2gbsbCiNZE1+RkUZSaRlRfHupQMRu12vnHsfeKqTLx/opJxVcXq8VIO2mBmQ2oml3s66SzRtMLQZE8+sDift3rctQMLf+Hae5DZSNGDizCZjJwrb8IyMoY11m1xhPTpA8QhPW4h1L3U/YDrM5avo50PEHe8TXdODXdfp/05T9UBdsn6PEJCzTz5pQf43l+/Rltjr093jb+/57q5FOiD0ZLOOTN4cwGJBSDckQSHBrFpexEXSq7S222B9f4HfaeiY8DC9y+WsDktm6XxyQxeGeVCXRvvnK7AsNbEH6zcSGVvF/ubrpBXGc5XntvGxZpW+ox2ugeH6Ryw0GcZpWt0mN9erWRpfDIPLl/I8eoGhqz6dMyuoWGSwsLpHHH19A8JDWLlmgWAyvGKZkat3ovX5orMgmRCw4OJiovgm8//Nw6P2Q3CnYlYAHcoM13qfjNWibciNM8mb94mh3laBoOb3dk02pRQ0Guolc+7Cpqyo2PYeT6cEwcqcVyr/tVqrxV/5JEHqSH+glujNllVIsKCWVOUSVpiFP1L+rjQ18RVSxcDDe7iKTVE/7L7XPcivvT57Xz/Jwd4Ob+MhJAIkkMjiQkOI6g6dMLHH5obxIcLlvFS9UVOvFDN+LW/NYOisG5HHiVXmlmXl4Gp2c7Z8iZsdie9RfqwnLaq2bPJmzaInX7Q/XXTiH6/jjD3uq6V+vtn/6W7H4+3lN7o+AiK71uK0+Hk0JtnsGe612kDxzC91M/5YhHfzYgFINwVbMrIwu50cmTvrbWNigg186mHV3O1uYezlc0cOlOLNcq/mUenTl+d+L89b5y20cGJRnTmCvfL7KvFW1mTnM6o08GDX8ijobWXk2X19A2Nsn1JHnbnOIcq6gi/4n2Yy1xjDgli+1PF9HUNcvjNs5PaYQh3NiIA7lBm2jfqqXX5sjC8aWg+J4dprQEPTVOn9XtYB6PXCsPMQUY+ayng0ustdHYNYfJoQaBNbjRa3X75hT/Q+7yvF5rt2JrHo/cu49tvHKLVMAoxim7yV7rG326N1adODuBkz7WCrsg9eu1VW4T149LT4FD48blShi6PkJccz7biHJZlpVJ2vpmC4Ggq+4z0pLqfFaY3gAj5sPsL4T/XVzpri8mijrr77nim32pTRMNerded01lzUxTuPfb5XfS09nDoHXfrBm0bh5txBHmbae35Mzf5mFQnzCwiAIR5T3JiFEuK0jh6ogbbDPmeT52+ykiykaOX62fkflPRNmzhb4662kVEYOBKRw8LU+L5+ObV9LYOsX1dAT39I/y2ynczubnioU9vZbDHwoGXT/rs+y/cuUgMIMCYjh/WtqtYd06bmeOt3YMnvsr2fc2RXZwXi9Fk5OIRfXFW17PLdf/WzrPVzr31zHrRapsVf6Gv5F30jUsTx9pCpdbN3qdgmfv1/45oc2vzkdXuk1c/4rIuIoLNbMvLofxcK6ty0zh6uZ6Q0+48+n6PATZai8KzyEzbQiLtqDtY4FlYpc3wGl+i/8xabV6rlS/69HqiIkPZd8z1fddaadpCM20zQJiZJoW3s2lcoCCFYMJtJSE9lj/8j0+TkD650Zo/RESH8Xt/8SFGLdZJL/87GcuYjXcuV9PRb2H32WosVtuNL7rNrNm6kK89fz/N7f7FRIQ7F3EBBRjefKqe2ppWkzO/V6q/h+ZY9wukcRN8+KuPsetTW8lemMLbP9zDJQt0tg8AoHi2fNZo5dfbE6/dUMj2Z9bScq6GmhMVk/bkmROv83Nr7qetHQBX5s91kk/pt6HV+rVE1eutZNOod6tZm7ePpg3zWIo+gyf3ZU1H0zp3A7WRzXrLJq7C/Tm1ef8Algy3SaBt92DwiMVoNerhDL0rJ6orZ+K4d1MqAJ/868dJS4rmG8/ewxdXXxv4p7EcDQf199cyExk8ovXfPkQACLPCK9/fx/jICC9/67fEJUfz8O8/hqqq1Fa10Xq+gcZa38G9YxfqUP93P+qVdpbcs4jyY3ePFTDf+ccX9vAXT97Dd37vR3O9FWGWEQEgzArdbf18/yv/4zpu7qHmRwdISYuhYFEqy9blkZmXzLhznO72Aa6crmXcqdfoh0dt7D1VRdS750gvSGXzU+s59UHFpHVaVt1TQFRMGKXvnWNkSFoUTJeKuk635i/c1UgQWLhpvBaCeRYSaTNHNG6ZMMcoy7YU4bA7aK1tJ3PZAgxGA4O9FkoywnWZPtdbLUSGBvO4M5ZTH1xmeMg6KXUw5am1bFiby1MfWs33f3KA/YddFkPzDr2Oo00Lve7yuI42YNz0lLsLafJpvTDxNyiccMl9ndZdA9Cx3v13F9rhdhtlv6kP4GqneXm6gHTP9eywqkEbtA3q8WjPoXE/+QrG61qFaNqEOMu8W2a+CgOF24sUggnzhsGeIY69UUJwqJnl2xZjCjJRfrIGdVxlzQP3Yg4yMWK1cflKO1ZcvvOh0TEO7j7HhvuX0FjbQXOZOxMpNDyY9LQY3nz3HO2dAxOFWoIg+EYsgABjOml6ns26fGl908Gcv4Ci1TnEJkZSs6eUlpo2QiNCWLQun5CNS7DbnVRXtzP4xlEAFm9cSHBRHhdL6wC4d9dyjuy9hNPhcg9pm4mNrM7WPUvbFkEXsEXfRE6rsQ/m6Bu+aSdsec4w9rQqrhO394ru39r02bE497NCj1To1lV8p2jiuPAn+p5CTs28YO1+tbN94Saa/nV6tyK8JQ9IG4c7A7EAhHnL+LhK+WnXyzwzPIitz2ygs7Gbc/svMT5uJijIyMKFKSx9vBh1XOXKhXpM0XFsfXAplsFRys81TLz8BUHwH7EAAhitZng7/bPa9s+g11iva5SJGXEsXL2AoY5eyg5X4LA7MOXmYDAaWLAoleT4MJIy4njo4/fwz1/6b+out0zcQ6tde6Y9aouzPDXe+t91p4xqUz8jmvWat7bFdMHP9Lnydk2LCq1f3nNKl9YXr/Xze6ItwHImRunONT3g1r5z/svdCsJzqppW6/es6NXO/tVadr6K9cSXf+chhWDCHUVXcy/H3jpD9ekrFO9aycbH1hIWEcy4c5wr5S0cf+8CKVnxZOQn86k/eZScoqlHRwqC4B1xAQUw09XkptNOQnsu9o0y/TrNsWerghHgZIeNoCAjKxclExoRQuWZOrrD4/jNy6dRw0J45X+OEpcSz+rCDM4cr8VY0WZ0ywAAC81JREFU79aao+r1e9e1nvYYdDOSrZmx2+5uAOeZwZNwwX18vcXDdTI/cN9Tl31zWt/vR6uJh4eZJ44NIx6ZPppsJ6NH5lP2m5pn++jVo8vU8lg3adjPFNcAmDQWgPbnr3rsSWIAdxYiAIQ7ArvdydlDlSiKwqI1ORQVZtB4tYsf/vO7AHR3DJKQFMW9Dy3j+EvHsNvmb4tlQZgvSAxAuC1oLQBfoyM9WxlrNWetdeAMNZGdHU92TiI93UNUfv9tVFXFHGJm22cf4PzRKvo6Bydrxl7GWXriaYnoPosm88dXLYHWZ6/N9AF9w7rYC+44gqdGrr2HZ8aRFq2ff5L2ron1XG+HfR1t1pHPhn1+tmie6UFFwswgWUDCXUdDQw8NDT3Ex0dwzxPFWIfHuHCwnAOvlbJ2x2L6u4eovdIz19sUhHmLCADhjqenx8LRg2cJjQhhzc4VBMXHUnaihqSMeNZtL6LkQMWNbyIIAYgIAGHW8FZ0ZttSpFundUMYR/W+e51bQjNvQJ2ie6kNKD1YhbEgi2VblhARFUp7n4VtH9/EsUOVjI05MGqucfpyZWgCyZ5uJK3bx7Poyrle383zOtpupQDaO2rvF9Kj//zawrVIj3kLWjfVmMa14/ksbbDf6uGyCtO64zTfU083nb8JA+L2ubOYkzRQRVH+VVGUSkVRLiqK8rqiKN4ToQXhJnE6xjlfcpWj+8pRrk1o//jntpKYHHWDKwUhsJiTILCiKDuB/aqqOhRF+WcAVVX/5EbXSRDYP6YTiJvUyE3DTGh1050U5e26mw1SpuensOWZDTiCgtnz4glXt1APzV4bINVOPfPVTkI3zxi9pqxtoeFZdKYtwNLi2WZD+zlHPSwnLaHN7vkCnoFk3Z78bN3g+fswV5q9BJVnhnlVCKaq6h5VVa/buieBDF/rBeFWaaltp72uk6e/cD/FO5bM9XYEYV4wH2IAnwN+PdebuJuYjqbk65qZaAvgS3ufTvGQryHlnve/zpkjNfzXv++htKzVpYFrWiEDmDWfRdXcw9Onrk0l9fy+aadvOTWTrXzNVfaF9nvjy7evLabz1N51hVueloeXxn6eKbIzMet3OojWP7vMmgBQFGUfkDLFqT9XVfXNa2v+HNeEwV/6uM/zwPMAIXj/oxeEGzEyOMqh3WU3XigIAcKcFYIpivIZ4PeA+1RV9V6Ro0FiADPLXGl1N4M3bd7fLJXpWhu+2h3oWODde+lv22xfRXLe4hKgb3pnKPc+A2E6sRNf67SxE2kMd2cwrwrBFEV5CPg6sM3fl78gCIIws8xVDOC7QDCw91qa3klVVb8wR3sJWOaj1q945NFr/ei6jBAfmqfPASgaPAfdaGMCWl+54uGG1loEnlq+dv+6kYgeVoRW29b9HDx+Jt7iEuDRGkJjOXi2jTZp1k36vnn5HbCvXaj7t+HgWff9fGSMCXcWcyIAVFXNv/EqQRAEYTaReQCCIAgBynxIAxUCHF0wWuPy8cRXSqDOdaRxeXgGVbWtIMY1ffhBrw1pC7duJhVR18/fR7roTOA18N2pdxU5fLj6vCUCaF0+k54rqZl3DWIBCIIgBChiAQjzCl9pm75mGKsay0H1N73Vw9rQFlPh5bmez550zkvqp7/T0nzu1zPl1MsM35sJ7s/HRADh9iEWgCAIQoAiFoAw5/irhfqb+umtIRv4X5zl73M9z3nbh+dnnI7G7mvvipe2zoLgC7EABEEQAhSxAIS7Am9aui+tebotj321KJ5Oc7yZ2JO0ZBCmg1gAgiAIAYpYAMKM4StbZr40nvOlvfu7R3+1cm+N7DzRNlrz2XhOEGYYsQAEQRACFBEAgiAIAYq4gIQZw1cg0l+3j98zaz3cTdqunzPhvvFVdDYdfH7+abrEZF6ucKuIBSAIghCgiAUgzCsmzfr1oh179vmfjoXheY2vvvzTwdeeZiIo7m2Cl7/PvZVnC3cHYgEIgiAEKGIBCDfNbKZ0+uvLnu5zZ2K/M+F711k6s2xtTGedEBiIBSAIghCgiAUg3DTzQYscv3e17t/aASazXXQ2Exk3s5m1I35+wV/EAhAEQQhQxAIQbjszoaH6Glk42xrvfM+/F41f8BexAARBEAIUEQCCIAgBiriAhNvOfHFRTNcVNdNun/nSKVUIPMQCEARBCFDEAhCEOUa0fmGuEAtAEAQhQBELQAhY5ovmLTEAYa4QC0AQBCFAEQtACCjmo7Y9X/YhBB5iAQiCIAQocyoAFEX5qqIoqqIoCXO5DyFwcFosE/8JQqAzZwJAUZRMYCfQOFd7EARBCGTm0gL4DvB1QJ3DPQiCIAQscyIAFEV5HGhRVfWCH2ufVxTltKIop+2M3YbdCYIgBAazlgWkKMo+IGWKU38O/Bku988NUVX1x8CPAaKUOLEWBEEQZohZEwCqqt4/1dcVRVkGLAAuKIoCkAGcVRRlnaqq7bO1H0EQBEHPba8DUFW1DJiYqKEoSj2wVlXV7tu9F0EQhEBG6gAEQRAClDmvBFZVNWeu9yAIghCIiAUgCIIQoIgAEARBCFBEAAiCIAQoIgAEQRACFBEAgiAIAYoIAEEQhABFBIAgCEKAIgJAEAQhQBEBIAiCEKCIABAEQQhQRAAIgiAEKCIABEEQAhQRAIIgCAGKCABBEIQARQSAIAhCgCICQBAEIUARASAIghCgiAAQBEEIUEQACIIgBCgiAARBEAIUEQCCIAgBiggAQRCEAEUEgCAIQoAiAkAQBCFAUVRVnes9+I2iKF1Aw1zvYxokAN1zvYnbSKB9XpDPHCjcqZ85W1XVRM8v3lEC4E5FUZTTqqqunet93C4C7fOCfOZA4W77zOICEgRBCFBEAAiCIAQoIgBuDz+e6w3cZgLt84J85kDhrvrMEgMQBEEIUMQCEARBCFBEAAiCIAQoIgBuI4qifFVRFFVRlIS53stsoyjKvyqKUqkoykVFUV5XFCVmrvc0WyiK8pCiKFWKotQqivKnc72f2UZRlExFUQ4oinJZUZRyRVG+PNd7uh0oimJUFOWcoihvz/VeZgoRALcJRVEygZ1A41zv5TaxF/7/9u4nxKoqgOP490djTKRt2hSO0GyHUQhsiFwkTgvLQdeKhrR1wIFCqNkG7apFggs3ghMhGNQm/INLKcRhZiGCRATNYCS4qI3Y1M/FvcEspnrw3n0H3/l9Vu++eYvfZeD87rmHey7TtvcA94APC+fphKRngLPA28AUcFTSVNlUndsA3rc9BbwOnKrgnAFOA3dLhxikFMDwfAacAapYdbd91fZGe/g9MFEyT4dmgB9t/2T7MfAVcKRwpk7Zvm97uf38B82guLNsqm5JmgAOAedLZxmkFMAQSDoCrNteLZ2lkPeA70qH6MhO4JdNx2uM+GC4maRXgFeBH8om6dznNBdwf5cOMkhjpQOMCknXgZe2+NMi8BHN7Z+R8l/nbPub9jeLNLcMloaZLbonaTtwGViw/XvpPF2RNAf8Zvu2pP2l8wxSCmBAbL+11feSdgOTwKokaG6FLEuasf3rECMO3L+d8z8knQTmgFmP7gMn68CuTccT7XcjTdI2msF/yfbXpfN0bB9wWNI7wDjwgqSLto8XztW3PAg2ZJJ+Bvbafhp3FOyZpIPAp8Cbth+UztMVSWM0i9yzNAP/LeCY7TtFg3VIzZXMBeCh7YXSeYapnQF8YHuudJZByBpAdOULYAdwTdKKpHOlA3WhXeieB67QLIZeGuXBv7UPOAEcaP+3K+3VcTxlMgOIiKhUZgAREZVKAUREVCoFEBFRqRRARESlUgAREZVKAUREVCoFEBFRqRRARB8kvda+82Bc0vPt/vjTpXNF9CIPgkX0SdLHNHvEPAes2f6kcKSInqQAIvok6VmaPYAeAW/Y/qtwpIie5BZQRP9eBLbT7H00XjhLRM8yA4jok6Rvad4ENgm8bHu+cKSInuR9ABF9kPQu8KftL9v3A9+UdMD2jdLZIv5PZgAREZXKGkBERKVSABERlUoBRERUKgUQEVGpFEBERKVSABERlUoBRERU6gl/Z+LtCBjq6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plot_samples_2D(chain, 100, 200, ax, lims=(-5.5, 5.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What you see is that Metropolis-Hastings takes a much longer time to actually find the relevant region of high probability and then explore it.\n",
    "This means that HMC produces much less correlated samples and thus will need fewer steps to achieve the same sampling quality.\n",
    "\n",
    "\n",
    "While this advantage is even more pronounced in higher dimensions, it doesn't come for free:\n",
    "numerically approximating the equations of motions of our ficticious particle takes up quite some computing time, especially if evaluating log-probability gradient is expensive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choosing the parameters\n",
    "Two parameters effectively determine the distance of the proposal state to the current state.\n",
    "These are the integration time step and the number of steps to perform.<br>\n",
    "Given a constant number of steps, the larger the integration time step, the further away from the current state the proposal will be.\n",
    "But at the same time, larger time steps incur less accurate results, which means lower acceptance probability.<br>\n",
    "The number of steps has a different impact:\n",
    "given a fixed step size, the more steps performed, the less correlated the proposal will be to the current state.\n",
    "There is a risk of the trajectory doubling back and wasting precious computation time, but\n",
    "[NUTS](http://jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf), a neat variant of HMC which optimizes the number of integration steps, addresses exactly this problem.<br>\n",
    "For a basic implementation of HMC, one can use a fixed number of integration steps and automatically adapt the timestep, targeting a reasonable acceptance rate (e.g. 50%) during a burn-in period.\n",
    "The adapted timestep can then be used for the final sampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drawbacks\n",
    "Unfortunately, there's no such thing as a free lunch.\n",
    "While the big advantage of HMC, namely, less correlated samples with high acceptance rates, are even more pronounced real-world applications than in our toy example, there are also some disadvantages.<br>\n",
    "First of all, you can only easily use HMC when you have continuous variables since the gradient of the log-probability for non-continuous variables is not defined.\n",
    "You might still be able to use HMC if you have a mix of discrete and continuous variables by embedding HMC in a [Gibbs sampler](https://www.tweag.io/posts/2020-01-09-mcmc-intro2.html).\n",
    "Another issue is that the gradient can be very tedious to calculate and its implementation error-prone.\n",
    "Luckily, popular probabilistic programming packages (see below), which implement HMC as their main sampling workhorse, take advantage of [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation), which relieves the user from worrying about gradients.<br>  \n",
    "Finally, as we just discussed, there are quite a few free parameters, the suboptimal tuning of which can severely decrease performance.\n",
    "The algorithm's popularity also means that implementations come with all sorts of useful heuristics to set these parameters.  \n",
    "As you can see, with HMC's popularity comes user-friendly software, which should make it easy for you to profit from HMC in spite of these drawbacks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I hope that at this point, you have a general understanding of the idea behind HMC and both its advantages and disadvantages.\n",
    "HMC is still subject to active research and if you like what you just read (and have the required math / stats / physics skills), there's many more interesting things to learn[^6].<br>\n",
    "To test your understanding of HMC, here's little brain teaser, which should also invite you to play around with the code:\n",
    "what would happen if you made a mistake in the gradient calculation or implementation?\n",
    "The answer is: you still get a perfectly valid MCMC algorithm, but acceptance rate will decrease.\n",
    "But why?\n",
    "\n",
    "I would also like to point out that there are several great, accessible HMC resources out there. For example, a classic introduction is [MCMC using Hamiltionian dynamics](https://arxiv.org/pdf/1206.1901.pdf) and [these amazing visualizations](https://chi-feng.github.io/mcmc-demo/app.html) explain HMC and the influence of the parameters better than any words.\n",
    "\n",
    "Whether you now decide to take a deep dive into HMC or you're just content to have finally learned what kind of magic is working under the hood of probabilistic programming packages such as [PyMC3](https://github.com/pymc-devs/pymc3), [Stan](https://mc-stan.org), or [Edward](http://edwardlib.org/)&mdash;stay tuned in order not to miss the last blog post in this series, which will be about Replica Exchange, a cool strategy to fight issues arising from multimodality!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Although hidden here, in general, the kinetic energy of a physical object involves the mass of the object. It turns out that HMC even permits a different \"mass\" (measure of inertia) for each dimension, which can even be position-dependent. This can be exploited to tune the algorithm.\n",
    "\n",
    "2. These assumptions are defined by the [canonical ensemble](https://en.wikipedia.org/wiki/Canonical_ensemble). In statistical mechanics, this defines a large, constant number of particles confined in a constant volume. The whole system is kept at a constant temperature.\n",
    "\n",
    "3. The temperature of a system is closely related to the average kinetic energy of its constituent particles. Drawing a new momentum for every HMC step can be seen as a kind of thermostat, keeping the average kinetic energy and thus the temperature constant. It should be noted that there are also [HMC variants with only partial momentum refreshment](https://arxiv.org/abs/1205.1939), meaning that the previous momentum is retained, but with some additional randomness on top.\n",
    "\n",
    "4. HMC combines two very different techniques, both very popular in computational physics: Metropolis-Hastings, which, in a sense, circumvents a real simulation of the physical system, and molecular dynamics, which explicitly solves equations of motion and thus simulates the system. Because of this, HMC is also known as (and was originally named) Hybrid Monte Carlo.\n",
    "\n",
    "5. The movement of a classical particle follows a set of differential equations known as [Hamilton's equations of motion](https://en.wikipedia.org/wiki/Hamiltonian_mechanics). In the context of HMC, numerical methods to solve them have to maintain crucial properties of Hamiltonian mechanics: _reversibility_ assures that a trajectory can be run \"backwards\" by flipping the direction of the final momentum and is important for detailed balance. [_Volume conservation_](https://en.wikipedia.org/wiki/Liouville%27s_theorem_(Hamiltonian)) essentially keeps the acceptance criterion simple. In addition, [_symplecticity_](https://en.wikipedia.org/wiki/Symplectic_integrator) is highly advantageous: it implies volume conservation, but also guarantees that the energy error stays small and thus leads to high acceptance rates. We use [leapfrog integration](https://en.wikipedia.org/wiki/Leapfrog_integration), which  checks all these boxes and is the standard method used in HMC.\n",
    "\n",
    "6. How about [using differential geometry to continuously optimize the mass](https://rss.onlinelibrary.wiley.com/doi/full/10.1111/j.1467-9868.2010.00765.x) or [viewing HMC in terms of operators acting on a discrete state space ladder](https://arxiv.org/abs/1409.5191)?\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "venv"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

mcmc_introduction_MH.ipynb

mcmc_introduction_MH.zip

mcmc_introduction_RE.ipynb

replace_inlinemath.py

#!/usr/bin/env python3

import sys
from itertools import cycle
import re

with open(sys.argv[1]) as ipf:
    lines = ipf.readlines()
# ## replace \{ with \\{ and \} with \\}
lines = [l.replace('\\{', r'\\{') for l in lines]
lines = [l.replace('\\}', r'\\}') for l in lines]

## replace \\ with \\\\
lines = [l.replace(r' \\', r' \\\\') for l in lines]

## replace ^* with ^\*
lines = [l.replace(r'^*', r'^\*') for l in lines]

## alternatingly replace $ with \\( and \\)
## if it's not part of $$
lines2 = []
for line in lines:
    if '$$' in line:
        lines2.append(line)
        continue
    else:
        cycler = cycle((True, False))
        matches = re.finditer('\$', line)
        offset = 0
        for match in matches:
            replacement = '\\\(' if next(cycler) else '\\\)'
            line = line[:match.start()+offset] + replacement + line[match.start()+1+offset:]
            offset += 2
        lines2.append(line)

with open(sys.argv[2]) as ipf:
    header = ipf.readlines()

with open(sys.argv[3], 'w') as opf:
    for line in header + lines2[2:]:
        opf.write(line)

Reads nicely @simeoncarstens !

Here are a few thoughts that I had while reading the beginning (will continue another day). Some of them are resolved by reading it several times, but maybe it is a good hint for where a reader might hang a bit:

A sufficient (but not neccessary), convenient condition for a Markov chain to have a distribution $p(x)$ as its invariant distribution is called detailed balance: $p(x)T(y|x) = p(y)T(x|y)$.

what is y in the equation?

A distribution $\pi$ of the state space of a Markov Chain

Calling the stationary distribution eigenvector \pi confuses me because I associate it so much with a number. But maybe it is the standard. In that case never mind ...

So all that is fun, but back to sampling arbitrary probability distributions. If we could design our transition kernel in such a way that the next state is already drawn from $\pi$, we would be done, as our Markov Chain would... well... immediately sample from $\pi$. Unfortunately, to do this, we need to be able to sample from $\pi$, which we can't. Otherwise you wouldn't be reading this, right?

I feel this is a nice paragraph but it could be easier to understand without prior knowledge: You don't really introduce what "design our transition kernel" means. Also, how can I "draw from \pi" I thought I only draw a new state based on the row in T of the previous step?

n $q(x_{i+1}^*|x_i)$, from

the ^* is difficult to see and read. Probably better if it has different support. Otherwise, can you use a different symbol?

probability $p_\mathrm{acc}(x_{i+1}=x_{i+1}^|x_{i+1}^, x_i)$

For this transition kernel to obey detailed balance, we need to choose $p_\mathrm{acc}$ correctly. One possibility is the Metropolis acceptance criterion:

how can we know when it is set correctly? Is this difficult to grasp intuitively?

I don't think the method you use for the mixture of Gaussians can be called Gibbs sampling. As you say, in Gibbs sampling you'd have to draw from p(k|x), but it does not equal p(k) because k and x are not independent. You'd have to use the Bayes theorem to calculate p(k|x).

Your method is to draw hierarchically: first draw k from the marginal distribution, then draw x from the conditional distribution, whereas in Gibbs sampling you draw from both conditional distributions, each time using the value from the previous step.

Your method works, but it works for a different reason than the reason Gibbs sampling works. I don't think it has a name.

Edit: moreover, I think your method works well precisely because you are not using GIbbs sampling. If you did, you'd have the same problem of being stuck in one of the mixture components, because p(k|x) would very likely give you the same k as the one that generated that x.

Finished reading — very nice! I corrected a few typos here: https://gist.github.com/feuerbach/5ccaeee166b45be13a8e375f980f405a.

Maybe one day I'll make a similar post except everything is done in Stan :)

Thanks for your feedback, @MMesch and @feuerbach! Highly appreciated!
For now I mostly worked on the MH part and put it into a separate notebook, where I addressed some of @MMesch's points.

@MMesch: I fear the \pi is indeed somewhat standard - as p often denotes other probabilities such as p_acc.

@feuerbach: This is very interesting - and you're right. The background to this example is that, a few years ago, I did something like that to sample from a different kind of mixture model and back then we used to call it Gibbs sampling. Looking this up, it seems like a very common method to sample from Gaussian (and other) mixture models (where you can easily marginalize out x to obtain p(k)). It seems to be called "Collapsed Gibbs sampling". I will rewrite this part to feature an actual Gibbs sampler. I just need to think of a nice, concise example where you can sample from the conditional distributions without resorting to MCMC. And as for being stuck in one of the mixture components: highly correlated samples often are a problem when using Gibbs sampling and I should discuss this.

I am also thinking of discussing how to assess convergence for MCMC methods, but this is a very difficult topic. It would make the series much more complete, but is also a really ugly rock to look under...

@feuerbach: Just FYI, I implemented the actual Gibbs sampler for this problem and it turns out that

you'd have the same problem of being stuck in one of the mixture components, because p(k|x) would very likely give you the same k as the one that generated that x.

is a slight understatement. If the sampler starts in the mode on the right, there's (on average) a ~1/1e6 chance for it to jump to the mode on the left. That was fun and very instructive and might even serve in a next version of that blog post / notebook as a negative example.