daily_lab_notebooks/day13-GibbsSampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# day13 of CS 136 : Gibbs sampling for a 2D multivariate normal\n",
    "\n",
    "## Outline\n",
    "\n",
    "* **Part 1: Defining the conditional distributions**\n",
    "* **Part 2: Writing a Gibbs sampler**\n",
    "* **Part 3: Sanity check results against NumPy's built-in sampler**\n",
    "\n",
    "## Takeaways\n",
    "\n",
    "* Gibbs samplers are only possible if we can exactly define the conditional distributions (for some models this is easy, for many others this is quite challenging)\n",
    "* When possible, Gibbs samplers are almost always a better choice than Random Walk or other proposal-based MCMC methods\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"notebook\", font_scale=1.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Background\n",
    "\n",
    "## Random variables\n",
    "\n",
    "(We'll use 0-indexing so we have close alignment between math and python code)\n",
    "\n",
    "We are interested in a 2-dimensional random variable $z = [z_0, z_1]$\n",
    "\n",
    "Each entry $z_d$ is a real scalar: $z_d \\in \\mathbb{R}$\n",
    "\n",
    "## Target distribution\n",
    "\n",
    "Our target distribution is a multivariate Gaussian with known zero mean and a covariance matrix that establishes positive correlation between $z_0$ and $z_1$.\n",
    "\n",
    "\\begin{align}\n",
    "p^*(z_0, z_1) =  \\mathcal{N}\\left(\n",
    "    \\left[ \\begin{array}{c}\n",
    "    0 \\\\ 0\n",
    "    \\end{array} \\right],\n",
    "    \\left[\n",
    "    \\begin{array}{c c}\n",
    "    1 & 0.8 \\\\\n",
    "    0.8 & 2\n",
    "    \\end{array} \\right] \\right)\n",
    "\\end{align}\n",
    "\n",
    "Our goal is to draw SAMPLES from this distribution, using the Gibbs sampling method.\n",
    "\n",
    "Naturally, you could draw samples from this distribution just by calling `np.random.multivariate_normal` with the right mean and covariance parameters. This is simpler (of course), and we'll use this fact to *check our Gibbs sampler*, since both procedures should generate samples from the same distribution.\n",
    "\n",
    "One way to think about what we're doing is to say, what if you were given a software library where it was easy to sample from a \"univariate\" normal, but not a *multivariate* normal. The Gibbs sampler here would be one way we could still draw *joint* samples.\n",
    "\n",
    "## Gibbs Sampling algorithm sketch\n",
    "\n",
    "To begin, we initialize the vector $z^0$ to *any valid value in the sample space*.\n",
    "\n",
    "Then, at each iteration $t \\in 0,1 , 2, \\ldots$, we iteratively sample each entry of $z$ as follows:\n",
    "\\begin{align}\n",
    "    z^{t+1}_0 &\\sim p^* (z_0 | z_1 = z^t_1) \\\\\n",
    "    z^{t+1}_1 &\\sim p^* (z_1 | z_0 = z^{t+1}_0)\n",
    "\\end{align}\n",
    "\n",
    "## Things to remember\n",
    "\n",
    "This is a simple example to illustrate the idea of how Gibbs sampling works.\n",
    "\n",
    "There are other \"better\" ways of sampling from a 2d normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Prepare for Gibbs sampling\n",
    "\n",
    "## Exercise 1a: Using the code below, complete the method to sample from p(z0 | z1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we used Bishop textbook formulas to compute the conditional mean/var\n",
    "# See Eq 2.96 and 2.97\n",
    "# https://www.microsoft.com/en-us/research/uploads/prod/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf#page=110\n",
    "\n",
    "# Start by defining parameters of the given target distribution\n",
    "mean_2 = np.zeros(2)\n",
    "cov_22 = np.asarray([[1.0, 0.8], [0.8, 2.0]])\n",
    "lam_22 = np.linalg.inv(cov_22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6799999999999999\n"
     ]
    }
   ],
   "source": [
    "# Implement Eq 2.96: Solve for the conditional variance\n",
    "# Note that this is a CONSTANT no matter what value of z1 we're conditioning on\n",
    "\n",
    "cond_v_of_z0 = 1.0 / lam_22[0,0] # grab diag entry of lam_22 relevant to 0\n",
    "print(cond_v_of_z0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement Eq 2.97: Solve for the conditional mean\n",
    "# Note that this is a FUNCTION of the value of z1 we're conditioning on\n",
    "def calc_cond_mean_of_z0(z1):\n",
    "    return mean_2[0] - cond_v_of_z0 * lam_22[0,1] * z1 + cond_v_of_z0 * lam_22[0,1] * mean_2[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.4       , -0.31111111, -0.22222222, -0.13333333, -0.04444444,\n",
       "        0.04444444,  0.13333333,  0.22222222,  0.31111111,  0.4       ])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_cond_mean_of_z0(np.linspace(-1, 1, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_z0_given_z1(z1, random_state):\n",
    "    ''' Sample a value of z[0] from its conditional given z[1] \n",
    "    \n",
    "    Args\n",
    "    ----\n",
    "    z1 : float\n",
    "        Value of random variable z[1]\n",
    "    random_state : numpy RandomState object\n",
    "        Random number generator\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    z0_samp : float\n",
    "        Sample from the conditional p(z[0] | z[1])\n",
    "    '''\n",
    "\n",
    "    cond_mean_of_z0 = 0.0 # TODO compute the conditional mean of z[0] given z[1]\n",
    "    cond_var_of_z0  = 1.0 # TODO compute the conditional variance of z[0] given z[1]\n",
    "        \n",
    "    # Draw a value u from the Standard Normal (zero mean, unit variance)\n",
    "    u_samp = random_state.randn()\n",
    "    \n",
    "    # Transform this value into a sample from Normal(mean_01, var_01)\n",
    "    # Remember, if u ~ Normal(0, 1), a \"standard\" normal with mean 0 variance 1,\n",
    "    # then using transform: x <- T(u), with T(u) = \\mu + \\sigma * u\n",
    "    # we can say x ~ Normal(\\mu, \\sigma^2)\n",
    "    z0_samp = u_samp # TODO update this line to perform the transformation\n",
    "    return z0_samp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1b: Using the code below, complete the method to sample from p(z1 | z0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the same process as above! Just flip the indices so you get the mean and variance for z1 instead of z0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_z1_given_z0(z0, random_state):\n",
    "    ''' Sample a value of z[1] from its conditional given z[0] \n",
    "    \n",
    "    Args\n",
    "    ----\n",
    "    z0 : float\n",
    "        Value of random variable z[0]\n",
    "    random_state : numpy RandomState object\n",
    "        Random number generator\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    z1_samp : float\n",
    "        Sample from the conditional p(z[1] | z[0])\n",
    "    '''\n",
    "    ## First, use Bishop textbook formulas to compute conditional mean/var\n",
    "    cond_mean_of_z1 = 0.0 # TODO compute the conditional mean of z[1] given z[0]\n",
    "    cond_var_of_z1  = 1.0 # TODO compute the conditional variance of z[1] given z[0]\n",
    "    \n",
    "    \n",
    "    # Draw a value u from the Standard Normal (zero mean, unit variance)\n",
    "    u_samp = random_state.randn()\n",
    "    \n",
    "    # Transform this value into a sample from Normal(mean_01, var_01)\n",
    "    # Remember, if u ~ Normal(0, 1), a \"standard\" normal with mean 0 variance 1,\n",
    "    # then using transform: x <- T(u), with T(u) = \\mu + \\sigma * u\n",
    "    # we can say x ~ Normal(\\mu, \\sigma^2)\n",
    "    z1_samp = u_samp # TODO update this line to perform the transformation\n",
    "    return z1_samp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Execute the Gibbs sampling algorithm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs sampler algorithm\n",
    "\n",
    "Perform 6000 iterations.\n",
    "\n",
    "Discard the first 1000 as \"not yet burned in\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = 6000\n",
    "sample_list = list()\n",
    "z_D = np.zeros(2)\n",
    "\n",
    "random_state = np.random.RandomState(0) # reproducible random seeds\n",
    "\n",
    "for t in range(S):\n",
    "    z_D[0] = draw_z0_given_z1(z_D[1], random_state)\n",
    "    z_D[1] = draw_z1_given_z0(z_D[0], random_state)\n",
    "    \n",
    "    if t > 1000:\n",
    "        sample_list.append(z_D.copy()) # save copies so we get different vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert list of samples to a 2d numpy array\n",
    "z_samples_SD = np.vstack(sample_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the sampled values for z0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEgCAYAAACq+TSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3zN1//HXzebSIwUEStGboIQyY3YI6gRNcLXLKoitKqoUbOlqNJfi1YoLR1aFK1Ro5TaUiMRI0ESIUQkkUX2Pr8/kntzx+d+7h6J9/Px8HDzGee8z+dzPud9zvu8z/sIGGMMBEEQBKEEC1MLQBAEQZg3pCgIgiAIXkhREARBELyQoiAIgiB4IUVBEARB8EKKgiAIguCFFIWWHD58GP/73//QsWNHeHt7Y9y4cThx4oSpxcKkSZOwbNkyre9PTk6Gu7s7rl27pkepKtFVPkPz6tUrvPvuu2jfvj1Gjx5tMjk+++wzeHt7w9fXF5mZmSaTAwBevnyJDRs2ICAgAB07doRIJMLEiRMV6vuCBQsQFBQEAAgNDYW7uztSU1OVpturVy98//33epX13LlzePvtt+Ht7Y2OHTti5MiR2Lt3r17z0AbpZ1MVsTK1AFWRffv2Yf369Vi+fDlEIhGKi4tx5swZzJs3D4WFhQgMDDS1iISWHDt2DDdu3MCePXvQsGFDk8jw4MED7NmzB2vWrEHXrl1Rt25dk8gBAImJiZg0aRJq166NuXPnom3btigsLMTZs2exZMkSPH36FO+99x4AYMWKFSgrKzOZrOfPn8eHH36Ijz/+GKtXr4ZAIEBoaCjWrl2Lly9f4v333zeZbFUdUhRasG/fPowZMwYjR46UHGvdujUeP36MXbt2kaKowmRlZaF+/fro0KGDSWUAgB49eqBRo0YmkwMAli1bBmtra+zduxd2dnaS461atUKdOnWwevVqvP3223BwcICDg4MJJS3/Lvv374/JkydLjrVo0QIvXrzArl27SFHoAJmetMDCwgI3b95Edna2zPFFixZh8+bNkr8fPHiA4OBg+Pr6wtPTEwMHDsThw4cl5ydNmoRvv/0WixYtQseOHdGjRw/s378fYWFhGDZsGLy8vDB+/Hg8ffoUAPDs2TO4u7vj6NGjGDx4MLy8vDBp0iRER0crlTUsLAzjxo1Dhw4d0K9fP3z99dcoLCyUnE9MTMT06dPh7e2Nfv364dKlS0rTSkhIgIeHh4JZKjg4GAsWLFCrzNIcPHgQbdu25T326tUrLFmyBJ07d4afnx+Cg4Px6NEjyflHjx5h6tSp8PHxgUgkwsyZM/Hs2TOlZcjMzMSnn36Knj17wsvLC++88w7u3bsHAFi8eDE2bdqE58+fw93dHQcPHuSUeeDAgVi5ciVEIhE+/vhjmfOMMfTt21emHgDADz/8gD59+qjscR84cACTJk0CAPTp0wfLli1DaGgovL29sXnzZnTq1AlTpkwBUP6sp02bBl9fX3Tu3BmLFy+WmKlKSkrg7u6OEydOYNy4cWjfvj3eeust3LlzB7/99ht69eoFkUiEBQsWoKioiFOWuLg4/Pfff5gzZ46MkhAzfPhw/P333xIFwWVeOXnyJPz9/eHl5YX3338faWlpMueTkpLwzjvvoH379hg4cCAuXLggOffixQt88MEH8PPzg7e3N9599108ePBA6bOztLTEvXv3FMxdQUFBMuanhIQEzJ49G507d0a7du3Qt29f/PTTT5LzCxYswIoVK7B27VqIRCJ07twZ27ZtQ3R0NMaOHYsOHTpgxIgRiIqKknnW+/fvx8iRI9GhQweMGjUKYWFhSmV98OAB3n33XXh5eaF3795YsWIFcnJyJOd79eqFL7/8EgMGDEDXrl15y20MSFFoQVBQEO7cuYOePXvivffew86dO3H//n3Uq1cPTZo0AQDk5eVh6tSpaNCgAfbv348jR46gU6dOWL58uczH8sMPP0ga/379+mHVqlX47LPPsHz5cvz2229ISUnBhg0bZPJft24d5s6diz/++AMODg549913FZQWANy/fx9BQUF48803cfToUaxZswbnzp3DypUrAQDFxcWYNm0a8vPzsXfvXnz++ee8NuOmTZtCJBLh+PHjkmMZGRkIDQ1FYGCg2mVWF8YYpk+fjhcvXmDHjh3Ys2cPXFxcMGHCBEmDuGDBAri4uODQoUPYvXs3MjMzsXTpUs70SktLMXXqVNy9exebNm3C/v37UbduXUycOBHPnj3DsmXLEBwcDGdnZ1y+fBkBAQGc6cTHxyMnJweHDx/GjBkzZM4JBAIMHz4cx44dkzl+7NgxDB8+HBYW/J/c0KFDJUrm0KFDWLJkCYDy+hQREYE///wTS5cuxdOnTzFhwgQ4OTlh7969CAkJQVRUFIKCgmSU0dq1a/Hee+/hyJEjsLW1xbRp03Dx4kXs3LkTq1evxsmTJ5UqcnFD16VLF87zNjY2aNy4MW95du3ahXXr1mHv3r3IzMxEcHAwpKMG7du3D/7+/vjrr78wdOhQzJgxQ6K4V6xYAcYY9u7diz/++AN2dnaYM2eO0rymTJmC58+fw9/fH0FBQdi2bRtu3boFBwcHuLq6AqisUyUlJfj1119x4sQJDB06FOvWrZPpcP3555+wtrbG4cOHMWnSJGzcuBGzZ8/GzJkzsX//flhYWGD16tUy+X/55ZeYMGECDh06BA8PDwQFBSExMVFBzqSkJEyaNAnt2rXD4cOHsWnTJkRHR2P27Nky1+3duxdr167Ftm3b4O7uzvucDQ4jtCI8PJzNnj2biUQiJhQKmVAoZCNGjGAxMTGMMcbS0tLY9u3bWV5enuSep0+fMqFQyG7cuMEYY2zixIls7NixkvMxMTFMKBSyP/74Q3Lsyy+/ZAEBAYwxxhISEphQKGS7du2SnM/KymIdO3Zke/fulaS5dOlSxhhj8+fPZx9++KGM3GFhYUwoFLKUlBR2/vx55u7uzhITEyXnz58/z4RCIbt69Spnuffv38/8/PxYUVERY4yx3377jfXs2ZOVlpaqXWaxfH/++Sdr06aNTPrSx65cucLatGnDsrOzZa4ZMGAA27ZtG2OMMR8fH/bVV1+x4uJiyTOKiIjglF1ctkePHkmOFRYWst69e7N169YxxhjbsmUL8/f357xfLJ9QKGRxcXFKr3ny5AkTCoUsMjKSMVb5Xvnukebq1atMKBSypKQkxlj5cxAKhSw0NFRyzdq1a1nfvn0l5WaMsQcPHjChUMguXrzIiouLmVAoZN98843k/M8//8yEQiF7/vy55FhgYCBbvXo1pxwhISFMKBSy0tJSybHExETWsWNHmX/Hjx9njJXXt6lTp8rIfOXKFYXnIq5bPXv2ZPPnz5fJc/z48WzRokWMMcYCAgLY4sWLWUFBAWOMsdTUVHb16lVWVlam9Nk9ePCALVq0iHXu3FnyXQ4YMICFhYUxxhjLzc1lO3bsYCkpKZJ7ioqKmFAoZEePHpWUw9/fX5JPVlYWEwqFLCQkRHLPrl27mI+PD2OMSZ71l19+KTlfXFzM+vTpwzZu3KjwbP7v//6PjRkzRkbuxMREJhQK2e3btyXPZsGCBUrLaWxojkJLfHx84OPjg9LSUkRFReHs2bP47bffEBwcjH/++QdOTk6YMGECDh8+jPv37yM+Pl4yfCwtLZWk07x5c8nvGjVqAACaNWsmOWZnZ6dgGujUqZPkt4ODA1q1aoWYmBgFGe/fv48nT57A29tbcoxV9Obi4uIQGxuLunXrwsXFRXLey8uLt9yDBw/GmjVrcOXKFfTp0wdHjx7FsGHDYGFhoXaZ1eXevXsoLS1Fz549ZY4XFhYiLi4OADBnzhysX78ee/bsQZcuXdCnTx8MHTqUM72YmBjUrVsXLVq0kByzsbFBhw4dEBsbq7ZcAoFAMnLkolmzZvDx8cGxY8fQrl07HD16FF5eXmjZsqXaeXAhnWdsbCy8vLxgZVX5Cbu7u8PR0RGxsbHo2rWrRBYxNWvWhLW1tcy8h62trVLTU506dQCUm//EE+oNGjSQGYEMGDAAxcXFSmWWrnvNmjVD7dq1ERsbi86dOyucB4D27dvjxo0bAIAPPvgAixcvxt9//41OnTqhd+/eGD58OAQCgdL83N3dsW7dOjDGEB0djYsXL2LXrl2YPn06zpw5IxlBnjhxAnfu3MGTJ08462jTpk0l+aj7Xfr5+Ul+W1lZoV27dkq/y8jISIWyA+WmVPH8WNOmTZWW09iQotCQpKQkbN++HR988AHq168PS0tLdOjQAR06dICvry+CgoIQHR2Nhg0bYuzYsWjYsCH8/f3Rp08fNGjQAKNGjZJJT/pDF8P3IQCAtbW1zN9lZWWcJg1ra2uMGDECwcHBCufq16+Pe/fuyZgBuNKWp1atWujfvz+OHTuGVq1aISIiAmvWrAFQblNWp8x8SH+s1tbWqFOnDvbv369wXc2aNQEAkydPRkBAAM6dO4fQ0FB88cUX2LNnD/bt2wcbGxuZe2xtbTnzLCsr43wPyrCwsFBIW57AwEBs2bIFCxcuxLFjx/TiGik9T6Asf/myyL9PVXVLGpFIBAAIDw9H//79AZTXV+nOjSosLS0V5JOWne98QEAAunfvjgsXLuDKlSvYvHkzfv75Z+zfvx/16tWTuS8nJwcbN27EhAkT0KpVKwgEAnh4eMDDwwN9+/bFkCFDEB4ejq5du2LcuHEAgIEDB6Jr165o3749+vTpI5MeV31QZTaUv6e0tFShfED5O+nZs6fErCiNdLm45oVMBc1RaIitrS3++OMPBRs0ADg6OkIgEMDJyQmnT59Gbm4udu/ejRkzZqBv374Su7p846wpkZGRkt+vXr3C48eP0aZNG4XrWrdujbi4ODRv3lzyLyMjA+vXr0dubi7atGmDzMxMxMfHc6atjMDAQJw7dw5Hjx6Fp6cnWrduDQAal9na2hqlpaXIz8+XHJOWxc3NDS9fvgQAifxNmjTBpk2bcOPGDWRmZmL16tUoKSnB6NGjsXHjRvz888+4d+8e5+Sfm5sbMjMzZSbDi4qKcPfuXUkZ9MXgwYPx8uVL7N69Gy9evMCQIUP0mr6bmxtu376NkpISybEHDx4gJydHb2Xx8PCAn58fNm/ejIKCAoXzycnJKtMQzzcA5aOg7OxsGfmkzwNAREQE3NzcUFRUhC+++AJJSUkYNmwY1q9fj7/++gsJCQmck8R2dnY4evQoZ6dCPNnu5OSE//77D7GxsdizZw9mzZqFAQMGSCaRdf0uxZPbQHm9unfvHu932bhxY0m9Bsrnk1JSUnSSwVCQotCQevXqISgoCF9//TU2b96M6OhoPHnyBKdPn8aSJUsQGBgIFxcX1K1bFzk5OTh16hQSExPx77//YsWKFQCgdKivLhs2bMDly5cRExODjz/+GHXr1sXgwYMVrgsODsadO3fwxRdfIC4uDtevX8eiRYuQnZ2N+vXrS7w+Fi5ciLt37+LmzZuS0QEf3bp1g729PXbs2IERI0ZIjmta5o4dO0IgEODbb7/Fs2fPcOLECRw6dEhyvmvXrujYsSPmzp2LsLAwPH78GMuXL8e5c+cgFApRu3ZtXLx4EZ9++ikePHiAJ0+e4ODBg3B0dJQxL4np0qULvL29sWDBAoSHhyMmJgZLlixBVlYWxo4dq9azV0Zqaipyc3Mlfzs4OKBfv37YuHEj/P39JWYcACgoKEBqaqpOaw4mTZqEzMxMLFu2DA8fPkRYWBg+/vhjtGvXTsYEoivr169Hfn4+Ro8ejWPHjuHp06eIjY3Fjh07MHz4cDg5OcHNzU3p/cuXL0dYWBju3LmDjz/+GH5+fvDx8ZGcP3jwIPbs2YNHjx7hyy+/xIMHDxAUFAQbGxtERkZi5cqVuHPnDhISErB//35YW1tzNr5WVlb46KOP8Msvv2Dt2rWIiopCQkICLly4gNmzZ6Nbt27w9vZGvXr1wBjDX3/9hcTERISGhmLevHkAwGtCU4cdO3bg1KlTiIuLw6effoq8vDzORZuTJk1Ceno6li5dipiYGNy5cwfz5s1DQkKCRqM1Y0KKQgs++ugjrFq1CleuXMGECRMwZMgQbNy4ESNGjMCqVasAlPco33nnHaxZswZDhgzBN998g5kzZ6J58+a4e/euTvmPGTMGq1atwpgxY8AYwy+//CIxxUjj7u6O7du34+bNmxgxYgTmzp2LTp06ISQkBED5sP+HH35Ao0aNMHnyZMyZM0fiesmHhYUFhg0bhqKiIpmesqZlbtq0KVauXIlTp05h8ODB2L9/v4y7qUAgwJYtW9C6dWvMnDkTgYGBiI+Px44dO9C6dWtYWFhg+/btAMo/vmHDhuHhw4fYuXMnp0+/QCBASEgIWrRogRkzZmDs2LHIzMzEnj17dLYH9+jRAz/++KPMsREjRiA3NxfDhw+XOX706FH06NEDL1680Dq/Bg0a4KeffkJiYiJGjhyJWbNmoX379ti5c6dGZjRVuLi44ODBgxg8eDB++OEHDB8+HGPGjMHff/+NqVOn4uTJkwouztIEBQXho48+wuTJk9G8eXN88803MueDg4Px119/YdiwYbh8+TK2b98u8VDasGEDGjZsiOnTpyMgIAAXLlzAd999p/RdjR8/HiEhIRLX08GDB2P16tXo0qULtmzZAqB8bnH+/PnYvn07AgICsHr1agQGBqJTp056+S5DQkIQGBiIxMRE/Prrr2jQoIHCdQ0bNsTPP/+MpKQkjB49GsHBwWjatCl+/PFHlaZfUyFguo63CKPx7Nkz9OvXD7t374avr6+pxSEIAuXrKNq1a4cNGzbo3cRoLtCIgiAIguCFFAVBEATBC5meCIIgCF5oREEQBEHwUq0W3JWVlSE3NxfW1tYaLSwiCIJ4nWGMobi4GPb29pwLC81WUZSWlmLKlClo3Lgx1q1bp9Y9ubm5nEvmCYIgCNUIhUJO13KzVRQhISEICwtTGZ1SGrEPslAoVBligYvIyEh4enpqfF9Vhsr8ekBlfj3QtsxFRUWIiYlRuo7DLBXFf//9h5MnT2LQoEEa3Sc2N9nY2CiN66MKbe+rylCZXw+ozK8HupRZmcne7Caz09PTsWzZMnz99decq40JgiAI42JWiqKsrAwLFy7ElClTeMMCEARBEMbDrExP27dvh42Njcyet9qgTgRUZYSHh+uUd1WEyvx6QGV+PTBEmc1KURw5cgQvXryQxDEShzaOiorC0aNH1U7H09NTKztdeHi4JAb/6wKV+fWAyvx6oG2ZCwsLeTvYZqUoTp48KfP3smXLUFpaqrZ7LEEQBKF/zGqOgiAIgjA/zGpEIc/nn39uahGMzqYzMdh0Jhbx66pnuGKCIKoeNKIwMzadiTW1CARBEDKQoiAIgiB4IUVBEARB8EKKgiAIguCFFAVBEATBCykKgiAIghdSFARBEAQvpCgIgiAIXkhREARBaMDiP+9g8DeXTC2GUTHrldkEQRDmxu83EkwtgtGhEQVBEATBCykKgiAIghdSFGbEvH23TC0CYULO3EvBqahkU4tBEArQHIUZcTAi0dQiECZk2q4wAKDIwYTZQSMKgjAzIhNf4crDNFOLQRASaERBEGbGW5svA6CRBWE+0IiCIAiC4IUUBUFUcYpKynDkViIYY6YWhaimkOmJIKo4m8/GYvPZh7C1ssAgz0amFoeohtCIgiDMgILiUq3vTckqAAC8yi9WeS1jDH+EP9Mpv9eZp+l5phbBJJCiIAgzwOOTkzqnoY7l6fLDNCw4cBufH7+vc36vIwUlr6eCJdOTHCWlZQAAK0vSoYTmfH78HhIy8tHOxREf9nMztTgK5BSUAABSswtNLAlRlSBFIYfv52dQWspw97OBphaFqIL8cOkxAOBkVLLRFYVAYNTsXksepeaaWgSTQN1mOV7mFSO7sMTUYgAAcgtL0GP9WdyIzzC1KISZ8Cg1B7cSXnKeI6cnw/Peb+GmFsEkkKIwY6KeZ+FZZj6+PPnA1KKYNSlZBRi06SKSXxWYWhSD0/frCxix5YqpxSBeM0hRmCkRTzPx1T/RphajSrD3+lM8SM7GnutPTS2KSSHTE2EoaI7CTAncGir5LQC1AK8jv159ApfadujXpqGpRSFMyP2kLKRmF6KXsL7JZKARBVHl0cY2n5JVgEV/3EFRSZkW+TFsPB2DF1n8pi5dV0p/cjgSQb+E6ZRGVSHgm0sYs+0/U4thlgz+5hIm/3jdpDKQoqgKmNmAgjGGR6k5phZDAU0e04ojUdgXloAz91M0zici4SW++TcWH+2n/UP0xb2kLFw3oNPGs8w8xKe9nh5L+oAUBaExv99IQN+vL+Dao3RTiwIA0KXfro0OLisrz7GgWPPRCBcZuUV6Sac6kVVQrFfnhB7rz6HPV+f1lt7rBimKasqRW4nILqgM6fBPVDKeZeon/MDtCvfMRxw9tNCHacgrMo17sbqTuVvDXuFkxU5yuigZfQXhe+vbS3pJx1CcuZeCrl/8i0Ijrkru+9V5dPniX6PlpwslpWVIy6neCxhJUVRDIp5mYs7vt/DxH3ckx6b/Go4h3142aL4JGXmYsOMaFkrla0zUbbf/fZxvWEE05Lmees6GWkex8mgUkl4V4EWW8RrDtJzKUVbiS/XfV15RiVbzTrqw8mgUfNecQa6ZrL8yBKQoqgCamEcycoskHlPJcpOt6gSN04XcipHEwxTzm78wBAIVQxhzWgAXnlQI18XHq1zP95+oZHRfdxZnH6g3l9T201NGX2dyMrJ8dJpXVH3jQJGiqGa8yK4ai862nHuo98lFY68jMLf9H/jK//fD8mcdmZhlJGn0w93EVwA0k/teUtUqY1WAFIUGrPv7AQZ/Y3x7clVbSMVUWP5jU7Lxf6ei8faOa3rKUHl+jDG0WHIcOy8/xuXYNBy781yjpB+l5mD7hTil50tKy8xGYWgihqp3VBXZd8OwCy7vPnvFe746PlMxtOBOA7bxNBjmgnRjEfH0JRhjKk0k+kLdhYFvbrwIAMjX854IyvJnDFh97J5WaY79/ipSswvxdpfmqGUr+7kUFJeh9bK/8WHf1pg/wF2r9NXlSXouckxkAzcTPchLWk4hFv1516B56Lu+ViVoRFHNURZAThkTfriKXf/FG0QWXSkoLkXwrjCFzWN081ziP5/H0TiL9a644d5rhNAhvf/vvEpnhKo28tQnxph7uRSbKvO32E36dYAURTVDvuErLtWsMofGpePTI1Fa5WVozken4vS9FKw5rv7oQBsZy8qYWTUC0cnZSs+l5xQipsJ5wFDvQ14BRSa+guvi4wgzo6jGhq6LT9PzsPnsQ5lj31UBC4O+IEVRBTCHWE8vsgoQp+fV2Pqy7evydLjsyj5rTqPburNqpyHtyimbtn4Y973y0BYDNl7UeNSoK5di0wAAp7VY1V5V4TL7/d+paORLeTqJv9OkV/nVbqtZUhTVDPnen74aY7+1/6Lf1xd48zIVupQwI7cIK/+KQnFppe/9y7xiiWuxumm/zDPc6upSntFNutSqbnXehzYTruY+R3HibhJm7blpkLRLyxiGh1zGuegXnOd/vPJY8lv8bLt+cRYzfq1e+1bQZHYVQJMG2dw/akMiEJRH2lx/8gG2TxLB1spSZbO4+tg9FJcy+LrWRT+Phui/4YKKO2TzE5NdUII6NW20E1xPGPrdP0nPk20wjVjX+Mo2c7dhlARQ3gG4/ewVbivxeJJW4tIyXohJ5bi66kIjCj3y9T/RiEzkd6EzNCVlxl2VykVMSg5vL1iMtt5Y8imHxlXGnFp66C7OR6eq7XcvnsNhDLiX9EqjVcCaNswv84qMvmpYn0zceQ2fHpEdfekb18XH8ePlx6ovNBMYqzQ9bjwdI3Nuz7Wn1WaTKVIUeqK0jGHz2YcmrxjiPZtNgXS7n5VfrJayUIeNp2Pguvi4UiUY/iQTAFAm03KLFYB6MjAAo75TPhegSqWpk03HVadf2600NWGVlCuzOGyIuZg5+ZCPRLz00F2jzx8ZCrNTFHfu3MHEiRPh6+uLnj17Ys2aNcjPN6/YPFyIG6QyA47/XRcfxxd/3+e9Rp+28l9C49Fx1T9a3eu9+jTm/B7Be426jbi4h5mvIkSCLo9+9l5+Wbl4msEfZJGrfGcfcNu6NUXV4i99It9Ib5DrOasiLacQmVpGyN0XlqDVfdLo8k1oMuqtzmZfs1IU2dnZCA4ORv/+/XH16lUcOHAAERER+Oqrr0wtmtoYYnGbdJLbLzzSe/rKWPFXFF7maR8f6tidJLWuKy1jCl4id5+9UtjzQhxssAp0LmWIeJqJLeceqr5QA4aGaBbgUeyRo2ljllVQjBfZuq1R8F1zBt6rT+uUhi6N8IIDt7W6b/O/sRqtkVEmYviTTL17DALlHccPDDSJL49ZKYrExET4+vpiypQpsLKygrOzM4YPH44bN26YWjSCA331oGb8GgaPT07KHBsachl9xV5WOmgGXUXURxkDt4bi/05pv/+5Pkx4Yo8cTfsxPdadVWtepaikDCv/ipKMHG7EZ6DFEu4ghAXFpfjufBxKDDjXIU26lqOZr0/HqHxv0l5kykbIo74LVfAY1BfH1eyM6YpZKQoPDw9s2bJF8jdjDKdPn4anp6cJpVIPGes4Y3qN/3PlYbrSODZv77iKgRUhMdSh/cpT+hILTzPy4Lr4OGJSlC8I81l9Gr+r6JWdua+ZOaYaj/A5ydUhKmlMSjYOcJhvTkWlqGWSySpQHjZE+j38HZmEn0PjseZ4uWn0h4uPwBgQFp+pcN/Wcw+x/uQDjcxKugzUtfkUDRGdwHXxcY3NduaC2brHlpSUYOXKlUhISNDY9BQZGalz/uHhyicdZ/94Dr4udiguZWhgbwkXBysUl1VOni757SJ+j8rB7sAGsLPi18VHY3JxMi4PWwbzb5wuHcdGWrYrD9NljmVlyXr7RMfEwCqz0m0zu6BEcu39tCII61lzllX6mPzvtLRy+/jh8HgAQMjxMLzj5YinrxTNVBm5RVh26C7cLBXdBUtKSpTmI32stFS2oXz58iXntUlJz5Fb4YHy4EE0kGYjeS/aEB4ejrKKCfSIWxGoUfEuo9MVG9i7kXeRWkv2c7p58yYsLRRbOL66pal80lpM5zQAACAASURBVDx58gThcs951IHyENg+zrYAgLS0Sg+xGTsv4uNudbXOPzk5GeHh5fM0cU/K5xFT09IQHh6Oly8rNrd6VLl6WSxvXEJ5HY199AThVmlq5ZX4/DnCw9WPCivOa/0fl3ArIUvhOB9ljOHTI+otJkx6XtmjL5Wrz1xsOxeL3vWUd6wA4NLTfBSVMthZCeDbyE5yXFna4eHhuPG8AK3qWqNeDUu91S9pzFJRpKenY968ecjIyMDu3bvRsGFDje739PSEra2txvlKP+DaTYVo3cBB9oIDxwEAf8Xk4a+YyonM/5b0rRiep0AgEOBcQnmD2cqjPZxr24GPURVpikQiSfqqEIlECjKJjzneugakVH587kIhRC2dZNIWiUS4+TQTyw+EYlQbe3z9Th/u9Cp++/j4AAdOSI6/EX8HeJxQ/oxz8tDI2RkiURs4pGQD/3CMbgQCTpmtrKxk8lFWLsujp4Diyp5tnTp1OK9t1MgFsTmpQMZLeHi4Q9S8Xvl7+fNvCASa9yxFIhEKK9Ju064DnGpV1KknGcBZWQ8pT09PNHeyl3nOPj4+sLK0UHivMrJLya8p8nWmefPmEImacaYdkVxuAnJycgLinwEABLb2irIASH5VgHtJr9Bb2EDy3rlwrnjvAPDU4hlw/Tbq1qsHkcgbdaLCgOcpaNmyFRAaXikvgGPPo4DYeDRt2hQiUQsZOZXh0sgFIpEb90mOe0UiEcLDw/FdWJbCcVX8EhoPQD1F0cilEXAvFgBgKVefuRBYWKiUYZTU/eP9HCW/ldUbkUiEUQeOo1m9mtjYz1GtMspTWFjI28E2O0URHR2N6dOnw8vLC1u3boW9vb1J5Oi/obzBi183ROW1Xb+QDfdgSu8HdfMWux0mvFIdkVRlmgKZ//SOPtLV9Z3sD3uG9/u0UnqeK8yKOZnINJElcOsVJL0qwOy+rdW+R3mYGcWc9fF9JL8qwBu1bMoVsZ7RZC2NNOoMXjUt+vOXqveXEW89XO6F58h/sZaY1RxFWloapk6dikGDBuHbb781mZLQBenPpcsX/yKrwLC7yolRtk/18sO6m+F0Rdl8TaaWHlWMAbuvPVGIv8OVi772CHiZXyQ1MculFPQ7L2UMlDXuSRVbs16RWsiojMjEV5i156Zkwl2dR/AkPbcif/WRnqNIyylEly/+xbq/H2iQAj/quJ5zYehXrk7ybT/V37yjMsxKUfz+++9IS0vDvn374O3tLfk3ZIjqXr0hufIwTeswxs+17J1oSn8lXhWxL/jd8pRVxKRX+SqvMVW7eO1xOpYdisQKNaPcArov2Np+4RF815xBcWkZPj3CrXz/uq3epkjd153FOT2tp9CFlCzdd0Ocufsmjt1JkvTCo56/wk6ZldWKD/5ctObhLaTrmngSXln8JW3R1fX8VX6xymCAuqzMLywpRWkZQ3ZBsdbrUrTFrExPs2bNwqxZs0wthgJv77iGlvXNe3Tz/JVmH7244bzxnFsB5hZWVnhlPWVxb/3s/RfYcekx/ny/m0YyKOP0PVn7sPzaFPHexH/efIavx3hxppGSJVsufSm17y8+QtRzxUnVyMQs3OM4zkXiy3x8djQKvq51YWdtqR/BOHihQhE80nEr2u8vVjas4jcUl5qL1cfuYUBb8byi8gdvytXWZWUM476/iuvxGfB3l3Uk0UUseTdvfeK+/CQGtXPGhZhUo2+iZFYjCnPmUap+93c2FOqYWgpLSpGsQrHM3Ve5Ull6T4vT91Lw+41yt0Zx4xv7ojy2k7I9sMtYebA+QL3FT8G7wnjPSysu+UV5YmbuvonwJxl6d0dU5lf/wZ6bGhm5GID2K/9RWVZtiXiaCb+1/6q8znXxcZSWMVx7lK7XOGXKlMCknXra/hb8ZpmVf0XhWqLyOp5TVILrFftpKIxwzHhF58moZJPstGdWI4rqSGkZw5jt/2GWf2v0EvK7wOoDdXrOs/dG4FQUv1eHdFC9Np9W9pIW/XlHaV58SmrwN5dwZl5v/BH+TLWAUqjaQIhvYya+2E2GQH7kxfcunlTs0ndeCzOMOvCtbZGnsKQUY7+/CkA95w11EMffkmbtifuSvSyA8nAmL7ILsVCNbWRTsqUbfdUt+c+h8UrPnbibpLeIs1VrVkp7SFEYmFf5xbj+OAOzUyJw69MBphYHAFQqCT50WSWsSQhvMX5rzyjusaHsYhNPJpt6LvtekvrrDKRRx7NGU17lKzoqSJuqgPINkC7FpmFgO2eV6YknwHXFdfFx2NsYztynby6aSbhyMj3pGYFAriGr+EO6rSsrYzgX/UKmB6oPj5k7z17KhNzmk1Gevl+fR5//O6fyXq4GwJBw7R4n/aj05dWkD0wtya7/nmh1nyajD2Vs+jdW63uN7S2maqW7OewoaW6QopDCEBVWnKL0hOzOy4/x7k83cDIyWXJs9zX1g48pY1iI9iHOH6XmIj6dPxqqqZD/bJXZv03dUJt6RKEtXKvHNUV+pKlJYxu4NVTlNdWt8S6V2pf9fPQLlZGRTQ0pCim2hmk3dOeDyY0oGGP45165gkiVcrlV17USAELj0tB57Rnkcuzja0y0XZikKzIjiiraOBsKbcKPv1GLO4pBYYn5NF6cnQNzePdaVsBWS09g1LZQPHyRjSk/3cCyQ3dV36QGuUWGCbRIikKKs/GGa/jEFf2v289xoyJQmrZ9pAk/XENKVqHaJoOqtKvaP1HJCsf4FubJ77amz36npg2lvBnM0GYxrtDVqsKPc0lUt6Y157Xm2ouvCpsYqUPE05eSoItxOroqi9kTqf9w5gApCoOQyhO//0FyZeMusz7AgG3KQxWL7vSFPnr30zXclD6bJ7qprmg6ca/L3h3acP1xhl7Ska6HrouPSx3XS/J6x5CjyDP3UjQq97dnddtnxEKcmZ4KZaiN00hRGJjKXqVi7dP1Q9R2kyRDff/a1NFtFyqji96I17zhkw/joU807VEfikhUODZ6m2r7u7YsOchtrthxyXibW6ni+F3FEaI2XIpNw0f7buGLv+9r5T2nLtMMtK5Fmv87VRl6RFzDbj97hf162M3PUJCiUIGmE9zyjYu4U8rVcdB1aK+vBl9fH542eks6Xg/X3gWqmCE1AmEMuPnUvPYovqFFmXRFvCcEF4ZatyFPUYVJ8KgGc2+qOBSRaNQdHg3FlnOVnSMLqY/m4z/ucF1uFpCiMDBiRcPVhurB2UQ9GeTsWvIWFX2ZpnQd9a4/qb8gbwQ3XDHL/L86z3ntHSPuy60NZYwZJJaaMS1u5mrek4cUhZ4pkptc/e1quW87V4XQ3fSk2/36Js8EoQUMia7Pt8yEPgTSG11VV+LT89Bt3VnVF2qIsSI+VyVIURgY+W0+pRufPB19p83NNXT1sXsmzT/knG4Ti/omIdM816UQ/Px2Vfc1Teqi7Wp6Y0OKQgX6aozF8xHS6X12tLJh1caV0lAeDtpSldxwCcIc0Pe8hA4RdnghRWEkxCOJmxzB0rRFNu6/csxMn1QZTnGs6dAEeu6EsbmuZNsAXSFFYSTEFqdcJTvR5RRqboZ6lqneRN5bm/kXYRHczPn9lqlFIAiNKOGJpqwLpChUoK/H/vxVAcrKmNIJUmVbmfJxK8G8XEEJgtANnaPFGsjBhcKMq6DV0hN6S6vdilNKNx0hM0X1w5wi2xJVg8k/XtctAZqjqPqYYmcqwnSQ8ieqC6QoCIIgqgsGMj2RojATnmaQzz1BEOYJKQqCIAiCF1IUBEEQBC+kKAjCQNBkNlFdIEVBEAaC3GMJY2OoOKGkKAiCIAheSFEQBEEQvJCiqCA2JVv1RQRBEGYMmZ4MzId7I0wtAlHNoMlsorpAiqKCB8k0oiAIompjqL4JKQqCIAiCF1IUBGEgZu6+aWoRiNcMmqMgiCoGxe8iqgukKAiCIAheSFEQBEEQvJCiIAiCqC7QfhQEQRAEL7QVqmHp417f1CIQBEHoBK2jMDCGcisjCIIwFuQea2A6Nq1rahEIgiDMElIUFdSpaW1qEQiCIHTCUKYnK00uzsvLww8//IB//vkHycnJAABnZ2e8+eabmDZtGmrVqmUQIQmCIAjTodGIYsGCBSgqKkJISAhCQ0MRGhqKkJAQFBcXY8GCBXoTKi0tDe+99x5EIhG6d++OkJAQvaWtDEahPgmCIDjRaEQRFxeHrVu3yhxr0aIFFi5ciIEDB+pNqI8++ghNmjTBlStXkJCQgOnTp8PV1RVvvfWW3vIgCIIg1EOjEYW9vT0uXbqkcPzixYuoWbOmXgR68uQJrl+/jvnz58POzg5ubm545513sHv3br2kr4ynGfkGTZ8gCKKqonJEkZubC3t7ewDA+vXrsWLFCixfvhzOzs4QCAR4/vw5mjRpgnXr1ulFoNjYWNSpUwdvvPGG5FjLli0RHR2tl/SVcTDimUHTJwiCMDQmm8z29fVF79698c0338DNzQ179uxBeno6kpOTwRiDs7OzTKOuK7m5uahRo4bMsRo1aqCgoEDtNCIjIzXOt6y0VON7CIIgzAkBgPDwcL2nq1JRMMYQERGB4OBgbNu2DTVr1oSTkxOcnJwQExODTz/9VGHeQhdq1qypoBTy8/M18qjy9PSEra2tRvnWOHUGWYWFGt1DEARhTjRxtIJIJNL4vsLCQt4OtlpzFF9++SVyc3MxdepUZGdXbhman5+Pc+fOaSwUH25ubsjMzERGRobkWFxcHFq3bq3XfOSxsqAlJQRBVG1a1DHMejC1Wse6devil19+gUAgwDvvvIOXL18aRBgAcHV1hbe3N7788kvk5eXh4cOH2LVrFwIDAw2WJ0EQBKEctbvRtWrVwo8//ghHR0dMnDgRaWlpBhPq22+/RU5ODvz9/TFlyhSMHj0ao0ePNlh+BEEQhHI0WkdRo0YNfP/995g1axbefvttzJ071yBCNWjQwCiL7KQRUFRAgiAITlSOKDZs2AAXFxfJ3zY2Nti6dSvc3d2xZMkSgwpHEARBmB6ViiIgIEDB/dXKygqbNm3C0KFD9bbQztTQiIIgCIIbrV19LCwssHr1aoP47JoCAe1IQRBEFcdQHV7yCSUIgiB4IUVRAZmeCIIguCFFQRAEQfBCioIgCILghRRFBWR5IgiC4IYURQUCmqQgCILghBRFBaQmCIKo6hiqHSNFQRAEQfBCikIMDSkIgiA4IUVBEARB8EKKogIaUBAEQXBDiqICr6Z1TC0CQRCEWUKKooLlQ9qaWgSCIAizhBRFBVaWZHwiCILgghQFQRAEwQspCoIgiGrC4NaG2UiOFAVBEEQ1wd7aMCZ0UhQV0AyFebNpbEdTi0AQZk9tO0uDpEuKgiAIguCFFAVBEATBCykKgjAy15b2M7UIBKERpCgIwsjUs7cxtQgEoRGkKAjCyFjQJllVmnurBppaBKNDiqIC2uHOfPiwb2tTi6AzV5coNy9ZUFWr0tS0sTK1CEaHFAWhEU3q1jB4HpZVpCXlU2jOte2UnjO3Tsnkrs1NLQJh5pCiIDTCwc5a8tuzsaMJJTE98we4m1oEvdC6QS29p1nL9vXrdVdnSFEQGvE/URODpn/r0zcNmj4AzO1c2+B5EFWLbRN9TC2CWUOKwkSMNnCDayjcGzpIfgvUXM/ex72+2unXqWmjVrrj/ZqpnaY5Yk7WJzMSxSQ0q1cTgzwbmVoMLBrkoXDsz/e7mUASRUhR6Mj6Ue1NLQInsZ8Pxnu9W+klrdo1Ks1NPdze0Oje3sL6mNlHs8npoJ4tVF7zxUjzfO7qYq6N85RurqYWQcKeaZ1NLYJRqe9gq3CsQxPzGP2SojACY32bKhwb0M7ZoHlaW1qglq1+4r6cmNNT63tbvGGv8T2q7NtVxf797Xhvpee6t1Zf4Q5pb7ze7vReLdW67tO3DL/RV9N6homEKg8DM0o+muLTrA6sLZU30csC2sj87cLjQKErpCgqMGQPz0LqKdewtkT8uiEamWNMTeM62ns6GXpOw5wZ5uWCwx90BwDYyH3wK4e1Uzsda0NvqqWFHcyjkQPveVOOmLZNFOEDf/2MpsWIOyf7Z3TVa7p8qPKOkx+BGNKbjhSFkWhmpN6RuWFpIQBjuvfYPBubxxBcUzo2rYNHawPQU85kJ684TAGXDOq8qe/e9oF307r6F0hPDPJ0hq2V+qNpdebE3qhVvpperLT/fF//CsMQ3mf6wvS1tZohbMj1sgU4OquH0WUxB399PegIALIfkelLpRkWFgKFBliT9Sj6NozImyg0fZ7uzvyjCU2Y298N73Z31egefY/G1TE9yX9Loub19CoDUN6pkKZrSye956EtpCh0RL438s9HvRWvEQDWVsZv3vTRk9cH5iFFJaZ4LPLvwiyUuNSbOfxBd1z62N/oMjg72uGTIYaf7zAF9extMN5PcX6Si/5tGigc++hNIee1ykajhvzeSVHoiPTH1s6FewGapUCgVePk56r/Xos+YWAKJhV5BAKgTk1r3mteBzQxhWjK8iFtVF8khbySEgjKe7PqTh6rUnJfjfbSaJhiYSHAkYq5HPnjxkCdb1ObRtjKQgBHO/Xqvq21Yv0wpwgFpCh0QD6Eg7LvR/qFq9uRtLGyQMgE5V4zfNhZ6/+1Hp/dAxcXKvY4fw3id2FkDPBwdsRuHlfHHhp4AHHhZsa2XTHNncob4RlqehVJI11lhnd0UTg/0MAedGJcndRTJJo4P/A1hm/UsuHsiAgAdG9teLOMvY1i463OSNDWqvL701S9qKMcTOGlRYqiAmXv36+F5r3668tU7zdgqaLCbR7vjQaOmru7bRjjhTPzys1f2po3bKwUq0U7l9poJtdQqLvgDuB3Bw3qoXrdBB8NtXhOpqJOTc1DjDOUr2/4qL8Q34xT7DzYcryvllq4JQPKe84z+7SS2OXtbS05v5caHL1iPqb1aIHhHRvzXjOjF7f30veTfCUeZfpCvGZhGk99VDWyWD28HX6f3kWvcpkDpChUoMkHJ244GzjINlxcYaUtLASIXzdEN+GkGOzpjNsrBmCkTxM0qatez0/ZiMWYfvv/EzWBv4eifdYYwQf1wZl5vXBopvqrZ9XV3fKN7sph7TCnvxvntZwdCj1aLeLXDcHHgzzweaAn/p7TEw0c7PQyz7P8rbacnRJ1sLe1QgcpT7jTH/WSOa+NfC61y+tcu4oYZkKpSXt1O12BPk3QqDZ/3b240F9BXk0wxRwbKQodkHdnUzYkNIYnpIVAILOCWh28mtThPL5iqPEmF5VV+oPvd0OXlkpGc3LfrKpvmC/kt660buAA72bqu4qq85HvmdYZiwbpFnDQXo1Q2GJZ1B0Z2llbok0j5YEgpb2hlKW4NEAxTAUAtKzP3yGTNjWJG23pOQy3itAyXItbVSH/TuysLLEnuDN+nuKncVrS5Zaul9Ihe5o51ZTIK8/Adg0ByE5Yj+ukXploHYUJ4fuwh3m5qPXhWwgU3SO5aK6mDVieCZ2bYcUwzRv3mhw2WEB/k4gD2zWEG4e7sDrujQ0c7bQKJ+HM0buuKbdCvaa1CSYJebLcNtEHnwd6Sv72blaX96NfP6o9Lizsw5vdlG6uWDzYA1Yc71KcNFedVLe3Ki+eg62VZHTA114NVhJTycHOGvHrhsisN9LEtAkAj9YGYJ2eQup0a/UGaks5YWgymS3dYdS0pn02rLweSNfZdaM6yFwT3KslatewVgin81p5Pd25cwcTJ06Er68vevbsiTVr1iA/P9+oMqwaLrtqdvVw7lW0Ct4jSqqF9HV8Fefs/D7o36ahzDHf5uW91cWDPXBiNncojbWB7RXMXergVEsxtoy+cHa0w/ZJvpwhCH5+t7KnxteoaBOo7SrPftQ2lhb4arQXRI0MV26lcHzDOyb7YllAGwzybIS3O6u/J8TYTs3Q3Im/B25tZYH3ereClRqruvXREW3r4qjWzK0qk+LuaZ3xgX8r7J7WWalZSpm4FhYCtXrVFgJgREcXhThtmq6n+GlKJ47zsiM1TZtusWLna/PbuTji9ooBeMOA3688ZqUosrOzERwcjP79++Pq1as4cOAAIiIi8NVXXxlVjk5SbqkMDM1UfJSq8NBigZK4Ov42rTPCl/fHe71blX+MRsBOypWzoaP6lTH288Gc/uCG4v/+10H1RagszzvdmuN/oiYGGaJvHOul8aR8/7YNEczhBSUQAA04AsRpgjq9y0EV3lI+HKYzVZ5LfDGIBNDeDNK0Xk0sHOghcX7gary12QdEWukIBAJsGueNbq1087ZzqsXvmCD9CHStcm+2rexAvvZzFImJifD19cWUKVNgZWUFZ2dnDB8+HDdu3DB43uoOc7WZZB3kqb77YtdW5bZYcY/RztqSs+fvr+PqVC5b7jAvF8SvGwIbKwvUrRh2z/JXHvlV/iO2trTA6hGeSq7mhssFUZr/iZpglA93vCh1vZ1srCwQs2YwlgZott5AEwK9m/C7j2rYUAzydFbLqUBVA8TVqIjv6e1eH/Hrhmi10tpSzhlDIDCs26b096lJhynQu9yrysm+slEXp2QGax6VIi/bd2/76PzN64LRw3AWFhYiNTWV81yLFi2wZcsWyd+MMZw+fRqenpo1Pobk8qK+cF183CBpCwBM7e6KwZ7OcFHRo9P1kxzTiT9Yn2MNa2TmFaOHm3aVU91G42OOGPzSfDXaS6v85dHWu8ZUCAQCvNm2IY7fTVJ6zZl5veFYo/wTvr1iALw++0fhGmN1Pg3Vy7WvCMbXwNEWXk1q4/azVxrdz+VxqKqzZwoFIq3IAMXnaWVpgRm9W+FcdKqMxaNDk9q4o+Ez0QajK4pbt25h8uTJnOd27dqFzp3LF2aVlJRg5cqVSEhI0Nj0FBkZqbFchSWVb+bevXuS32lp6Sh8USj5Ozw8XOb3kyd5kr/z8nIl53cMrY9pR8sV4q2ICJSUladfWlYmk4Y04TfDJRVbefNQTtarLE6ZpHmemCP5LUBlo7GoWx0g7THC0x7LXJ+RkSFJq7CwvMyRkZHIcOCuJvl5eQrPIz2vFABQXFSsvJzh4RjUqiZOxuXhQeRtyfFWda2U3iNNWWkpwsPDEZNcLmNWFv+zUCdNbZBP9+nTPKXnUpKzAQCJic8QHv6SN92IiAjYWArw+En53Jz0e5HnFYAEALnFZTLHHz9+jPDSZLCyMoV7CguLAABP4uMRLnghK2duCQCguKgI7/s64lGm8vcoTU5ONiqqOKJjYlBaUsJ5nXx9UYcP/Wqjc+NiWLjaocDbRq13/Px5ed3PyEgHABQVFWGkhz0OPsiFg2UJwsPD8SK3vK4WFhUiPDwcL1+Wv5eHcY8QXvhcJr2Ciu/hwYP7KH1R3qg/zChWkONWxC1kFZVJ8kxOTgZQ3o4AQAN7SwVZ3xc5oqdzAcLDw5FdcW9pSYnCddYA/hztjMSHUUisOJabmys5X1RUxPks9IHRFUXnzp0RHR3Ne016ejrmzZuHjIwM7N69Gw0bNuS9Xh5PT0/Y2mpm480vKgUOnQQAtGnTFjh9CQDg5OSEAT28gGPlowiRSAQcqPwdW/oUCL8LAKhZ0778fAWWx0+gtIzB28cbRSVlwOF/YGlhIXONOC0AEPmIVHscVVxfu7YjkJxaKRMHV7MeApHlz7pzy3q4+igDAPD+MKmFSlL516tXDyJR+doKz8gwJN9PQSdvLzjLx7mvuEdSXqnnkfQqHzh+FtY21opySV0nf+pssxzUd7CV2ZNbWdktLS0hEomQE5MKXLoOR0dHBLk5wMpSAJGojUyZuJ6POh9S03o1kJDB70Qhn+69onggIorznEvrfMTtCsecYZ04N6gBIJHb29sbdtaWSLBIBK7fQt26le9FGVkFxcDhyhFFixYtIOrYGIKDfwOQVRa2tjZAbj6au7pCJJI1QSZk5AEnzsHG1gaL/qd6H5Kv8QzzD9xGrVoOKGMMSMuEh7s7LK+FAcWKDalIJAL+OFH5Ww14L5OqU9JcyowB7sXCyckJiH8GGxsbjOzuiYMPrsHR0QEikQjPMivKam0DkUiEOvfCgMQUtGrZEiKx2a8ifTtbWyAnDx4ebSTB+6yfvQT+vSKTr7e3NzLyioDjZ2FrYwNnZ2cg5hHqv/EGEJ+AWjXsKmWtSHvR6Mrn/DKvCDhyGpZWVmo9n1r/XQYyy0cUNjY2nM9CHQoLC3k72Ga3A0x0dDSmT58OLy8vbN26Ffb2uk0ka4qmw06+ETfXhCLfJJ8meWs6WTjKpwkWD26Df6KS1bp+07iOuPX0paKSkJFB8Zg4BIGmazpa1tctDMcnHBvp2FhaoKhUsUetDhcX+uP3GwlYcvCuxvdO7KK4TWuj2jVw9EP1IgiLn6s+TCC9hG/gzP0XaFXXCnGZ5b18vvm4xnVqYELnZpjcVT0vrEZ17CSyMqlHbS4BKflQ9g2p+9y5nmO515Ni2WtUzMWpW8/VfX7GespmpSjS0tIwdepUvPXWW1iyZImpxQGgnxehqT+4vnmvdyuMrpi8lg9lrIxatlYqtz3lqssNHOzw2bB2Ml4apuLCx33w/KV2rtUCnkCOjevUwBIlC8f0iSZtrbIaFjLBB8mvCpD8+D7G/ZkCAJjdzw17rz/ljBFlYSHA2kDN1yJIy8rX0Joyaq6yeTN1HjPnmhOOo9IL5aTL2qi2HXZN9YNPc/Pdx4MPs5rh+/3335GWloZ9+/bB29tb8m/IEP2FutAUjT5WLb4BVQundMFY3+TAdg1ltv18p5urysl4Xdg4tqNa1zWqXUPlvgGB3o3xj4bhFP4naoK3OigG59NX706+Y6HLe7SztoTrG/awljJpejg74M/3u2k86tOED/tyhxvRN+90ba6wJagqxM9Xk8eq7rOysODuYDAAvYT11d7G1xzC0EtjViOKWbNmYdasWSbJW/q96NvNj++dq1o4pQvGGv1vn+RrnIwq6NdGf6MVeaWz8x1fieOBtvXA1CNIVRhit8By99hKgnu1xOcn7us9lXgn7QAAFj9JREFUH3k+G669R6T4uxSvCVHne5EJ0WHm71mfmJWiMAd0efVN5YLxMSW/leZtgl7E2r71cDbJGmfup6h9T6v69ohLzcUYNWPQmCPWlgIE91Rc7KZPJaQt8tVAG4VfBaYIjIL0oxQrBHGv3tnRDrP7uUnWWqiD9GNtoMFiVK4v+9DMbkpHGNrM8RjylZOi0BNO9jZYr2SlsEDJb11QZwirjt5xd7JBzfrOOHM/RW0Th0udGohLzTX6PuDXl/ZDVgG326WmxH4eoPa16nhAGQIzsz5wIjbJtHjDHjEpOSquNh5c7WznFvWweLCHJMieQCDAPCW7yMnD9SoaOtohfHl/iNac0UpGrmCS5jpKMas5CnNAvlevyvwwzMsFA9o2xIk5PdW2P+rC4sEe+O5tH0kFl1+oow3ila59OcJ9c9GrYhFeUyOHAm/gaGeSDeh7tJZddKisRgzydEZDR1u8o0UwQ2n00VRos4+KNrRzqY2f3+2EFUMr46GZZ1NX/m2/17uVGnuCKJZA2Tt3qmWLP9/vhpb17RX3IjfXB6EFNKJQQTsXfnuuva0Vvp/Mb6MXCASoZWOFfh4NdN6k573e5Ru5pGQVANDPdonChg54sHoQ7NTceGZazxYY7u2iVSDC3kL9hSHo2tIJI70bK91bWFc0Hf03cLDDtaX9TZY/UB4SJWrVIL3JoA593Ms7GPLmkhrWlsgvLjWqLGK0aaTbudTGP/dSFFzC3+3uivAnmeXpctwnal4XZ+f3kfwt/RjE2wA7GtBxwBiQouBA+kVP7e6ql3QsLATYyRFt0hioMymrrpIAyhWfNkri5idvwt5Wf3tH21hZYIOaHlD6RGTGLo7m4C0jFuGtDo1wIPyZSWTQRsnO6tsa/h710UFqnxZxPKthIZc1Tk8gAIJ7locEH6PFPhka52fAtElR8NCmkaNePjxDvkC+78Hc7J319GAmMzW3Pn1Tq+1MNUFc5/Td5jd0tEVKVqHqC7XEHOfPNXmGlhYCGSWhD6wtLTQKIW+ukKKQQ5/fpnRseqLqIt3YGFpJGJIz83ojr8gYpiDDKDptqFvxvnQ1+WqCIaPo8udrOGgymwdd67kxlASfjKaqsNUFUz09Q7WvDnbWaodm14aeFU4O4r00Gtcp94ob6qW4ONFY2FlbIn7dEEzjcIXWBLErtasG6550GdGb25dLIwo5DNELMoeeFaE99PrUY24/N4z3a4pGtWW94VpoucWvOTHUy0VthVcdrQg0ouBBuoFfM8ITO1R4N8kztbvhhrvqrSKlJk4XRvk0xlsdGmFuf8N4VSlDYcGdUXPXHgsLgYKSAKqO/PpGlw6iuX25NKLgQboxnthF8wmpT95qg+VD2piFJwrBz7aJIoXtP2vaWCFkgo/RZPBuVgcRT/n3qlCGuTfGCwYI8dU/MaYWwyjo412Y2/skRVGBYUxOAoOZnUj3KOeLke1xOTZNo3s02a7WUPwy1Q/xabkKHYuq/qoFAGb1dcMsIwUKBMyjodXqvenwsh359nLRETI9ySFtrlHWGOtjkZtRMYevxoiM92uGLW8bbySgLxztrLV2zzTHGmkOzhSm6FA5VzgMGNNkOdK7MX5613DrtGhEUYG0mUmV/f/Kor5IzTacP7o6qDVHYY6tRxXG2tJ8H6jpm2RCTA0bS8lCPa3R8IWKQ/ur2kJZW0hRyCPVFihraJ1r2/Hu/GZMSBkYhxOze8KplmnWUFRVJWBKZ4oaFZEGNIk4QCiHTE8ciCNieqqI86Qvfg3yw0gNQh1LwzeyqI5ueqairYujQdcgVEeUmZ70EchSFVN7uGJ8u1p4V4cQPCZFTR07700h7G0s0crAwTJpRMFBM6eaODSzG9o0cjRKfj3d6ksWKxGELpjlAFNu2PvXhz0QmfjKoFnaWlnif21rwdaqeo8o+rg3MEoQSFIUSuCKFW+O8O9PbDw5CMNR3V5j4zo1FFyRCTnMzBpAikKO6vZRElUfddqMWjZW8GlWB7P6tja4POrybrcWuJ+UhXd13J/jdcJcO3ekKCogez5RlbGwEODgzO6mFkOG2jWtjb6fOmEYaDJbDnPV6PLUs7eBg50Vlg9pq/Ja0oEEUTWwqdjX29fVvEzfNKKowBwWB2mCjZUF7q4cyHtNFdF5BEFUYGdtib/n9ERzMwukSIpCDgqkRxCEKTGWt6UmkKKooIa1JUSNbDE3wMvUohAEAMCrIpzH0A6NTCwJ8bpDiqICgUCApT3qQkTrGQgzwfUNe91DQRCEHqDJbIIgCIIXUhQEQRAEL6QoXgMYLRIhCEIHSFFUY6rKmhCCIMwbUhQEQRAEL6QoCIIgCF5IURAEQRC8kKIgCIIgeCFFQRAEQfBCioIgCILghRTFawAtoyAIQhdIUVRjKBIuQRD6gBQFQRAEwQspCoIgCIIXUhTVGO9m5fsZdG3lZGJJCIKoytB+FNUYX9d6uLtyABzsrE0tCkEQVRgaUVRzSEkQBKErZqsoSktLMWnSJCxevNjUohAEQbzWmK2iCAkJQVhYmKnFIAiCeO0xS0Xx33//4eTJkxg0aJCpRSEIgnjtMfpkdmFhIVJTUznP1a9fHzk5OVi2bBlCQkKwe/dulJaWGllCgiAIQhqjK4pbt25h8uTJnOd+/vln/PDDD5gyZQratm2rdR6RkZFa3xseHq71vVUVKvPrAZX59cAQZTa6oujcuTOio6M5z3333XewsbFRqkjUxdPTE7a2thrfFx4eDpFIpFPeVQ0q8+sBlfn1QNsyFxYW8nawzWodxZEjR/DixQv4+voCAAoKCgAAUVFROHr0qMr7WUX0u6KiIq1lKCws1PreqgqV+fWAyvx6oE2ZxW0mUxJBVMCUnTEDli1bhtLSUqxbt06t67OzsxETE2NgqQiCIKonQqEQDg4OCsfNakShK/b29hAKhbC2toZAQJFTCYIg1IExhuLiYtjb23OeN+sRBUEQBGF6zHIdBUEQBGE+kKIgCIIgeCFFQRAEQfBCioIgCILghRQFQRAEwQspCoIgCIIXUhQEQRAEL6QoAKSlpeG9996DSCRC9+7dERISYmqRtCY+Ph5+fn549uyZ5Ni5c+cQEBAALy8vDB8+XGafj5KSEqxduxZdunSBj48P5s+fj5ycHMn5R48eYdKkSfD29oa/vz9+//13o5aHjzt37mDixInw9fVFz549sWbNGuTn5wOovmU+c+YMhg0bBm9vb/Tp0wfffvstysrKAFTfMovh2sysupb5zJkzaNOmDby9vSX/Fi5cCMBEZWYEmzhxIlu8eDHLz89nMTExrE+fPuzo0aOmFktjzp49y7p27cqEQiFLSEhgjDH2+PFj1qFDB3bu3DlWVFTE9uzZw/z8/Fhubi5jjLFvv/2WjRgxgqWkpLDMzEw2depU9sknnzDGGCsuLmZvvvkm27RpEyssLGQ3btxgvr6+7ObNmyYro5isrCzm5+fHfvrpJ1ZcXMySkpLYyJEj2apVq6ptmePj41m7du3YuXPnJH93796d/fHHH9W2zNJs2rSJeXh4sEWLFjHGqm/dZoyxjRs3so8++kjhuKnK/Norivj4eCYUCllqaqrk2E8//cTGjRtnQqk0Z+vWrSwgIIDt379fRlFs2LCBBQUFyVw7dOhQtn//fsYYY7169WLHjx+XnIuKimKenp4sJyeHXblyhXl7e7OioiLJ+dWrV7P58+cboUT83L9/n82cOVPm2C+//MKGDh1abcvMGGPZ2dmMMcbKysrYzZs3mZ+fHzt79my1LjNjjIWGhrJBgwaxuXPnShRFdS5zUFAQ27lzp8JxU5X5tTc9xcbGok6dOnjjjTckx1q2bKk0FLq5MnLkSBw7dgxdu3aVOf7w4UO4ubnJHGvVqhWio6ORnZ2N5ORktG7dWuZcUVER4uPjERsbC1dXV1hbWyvca2o8PDywZcsWyd+MMZw+fRqenp7VtswAUKtWLZSUlKB9+/YYN24cOnfujJ49e1brMqenp2PZsmX4+uuvUbNmTcnx6lzmqKgoXLlyBX379kXPnj2xfPlyvHr1ymRlfu0VRW5uLmrUqCFzrEaNGpIQ51WFhg0bcgZCzM3NhZ2dncwxOzs75OfnIzc3FwBkym9rawuBQCA5XxWeTUlJCT755BMkJCRgzpw51b7MlpaWuHnzJk6dOoXHjx9j1apV1bbMZWVlWLhwIedmZtW1zNnZ2WjdujX69++PY8eO4cCBA0hISMDChQtNVubXXlHUrFlT4UHl5+ejVq1aJpJIv3BVhIKCAtSqVUtSaaTPFxYWgjGGWrVqVYlnk56ejqCgINy+fRu7d+9Gw4YNq32ZBQIBbGxs4OrqipkzZ+L48ePVtszbt29XuplZdS2zg4MDfv31V4wfPx41a9aEs7MzFixYgIsXL4IxZpIyv/aKws3NDZmZmcjIyJAci4uLkxm+VWWEQiEeP34sc+zhw4do3bo1ateujQYNGuDRo0cy56ytreHq6go3Nzc8efJEZt9yc3o20dHRGDlyJGrXro3ff/8djRs3BlB9y3z27FmMGjVK5lhhYSFq165dbct85MgRXL9+Hb6+vvD19cWRI0dw7NgxDB06tNqW+dGjR1i/fr2MbIWFhbCwsECHDh1MU2aNZ1mqIWPHjmWLFi1iubm5LDY2lvXp00cyOVTVSEhIkJnMjo2NZR06dGCnT5+WeEmIRCKWmZnJGCufHBs+fDhLSkqSeEksWLCAMcZYUVER8/f3Zxs2bGCFhYUsLCyM+fr6stDQUJOVT0xqairr1q0bW7t2rcK56lrmtLQ05ufnx77//ntWUlLCYmNj2YABA9j27durbZnlWbp0qWQyu7qWOT09nYlEIrZ582ZWVFTEnj9/zsaMGcOWLl1qsjKTomCMpaSksA8++ID5+fmx7t27sy1btphaJK2RVxSMMXb+/Hn21ltvsY4dO7LAwEAWFhYmOVdYWMjWrVvHunfvznx9fdm8efMknjWMMRYXF8cmT57MfHx8mL+/P9u3b59Ry6OMzZs3M6FQyLy8vFjHjh0l/wICAhhj1bPMjDF2+/ZtNm7cOObj48PefPNN9uOPP7LS0lLGWPUtszTSioKx6lvmiIgINn78eObj48M6d+7MVq1axQoKChhjpikzbVxEEARB8PLaz1EQBEEQ/JCiIAiCIHghRUEQBEHwQoqCIAiC4IUUBUEQBMELKQqCIAiCF1IUBKEm165dg7u7O5KTk7W6/+DBgzLxih4+fIjz58/rSTrlnD9/Hg8fPgQAPHv2DO7u7jJ7GBCEKkhREISRCAgIwMWLFyV/z5w5E3fv3jVonikpKZgxYwbS09MBAI0aNcLly5fh5eVl0HyJ6oWVqQUgiNcFOzs7mcifxljrKp+HpaUl6tevb/B8ieoFjSgIs+XgwYMYPHgwPD094e/vL7PtZ1lZGbZu3YoBAwbA09MTvr7/3979hzS1v3EAf+92FVErk4KwyMKaxJjaapoyk/wjy0poQv3jwiGmsKaM/aFYRg0ZsyKt+SNWYFmUBCFR5I8/ImdFImv203DtNGc/9A+xJArMfO4f93poaXNJ8N393ucFgjue83zOeQ6eZ/vscJ5N0Ov14sMde3p6IJfLYbfbkZWVhYSEBOTn52N4eBgmkwkbN25EWloabDabOF55eTnKyspQWVmJDRs2QKVSoa6u7qcX9ImJCVgsFqhUKigUCuTl5aGvr8/v8UxPPWk0Gni9XtTV1SEzMzOgeFarFRqNBiUlJVAoFKipqZkzDxkZGQCA/fv3o7y8fMbU0+TkJM6dO4dt27ZBLpdj9+7duH37ts+YBQUFqK+vh0qlglKpRHFxMUZGRn7tZLJ/t9/2cBLGfqP+/n6SyWTU0dFBb9++pc7OTkpISKDW1lYiIjp//jylpKRQV1cXvXnzhu7evUupqalUVVVFREQPHz6k+Ph4ys3NpSdPntCjR49IqVSSUqmkEydOkCAIVFtbS1KplFwuFxERlZWVkUwmI4PBQC6Xi27evElJSUl09uxZMaZUKqX3798TEVFpaSmp1Wrq7e0lQRDIarWSXC4nQRBmPabr16/T+vXriYhobGyMtm7dShaLhUZHRwOKd+bMGZJKpXT8+HHyeDzk9XrnzMPz589JKpVSR0cHjY+Pi88C6+3tJaK/O5ylpKRQW1sbCYJAjY2NFB8fT+3t7eKYMpmM9Ho9DQwMUHd3NyUnJ1NFRcXvPeEsqPHUEwtKQ0NDkEgkiImJEX+ampqwfPlyAMCaNWtQXV2NLVu2AABWrFiB9PR0DAwMiDGICAaDAXK5HACwefNmPHv2DEajERKJBEVFRWhoaIDL5RIftbxkyRJYLBaEhoZi7dq1cLvduHTpEg4cOOCzf4ODg2hra8OtW7fEjmMHDx6Ew+FAU1MTTCaT3+OLiorCggULEB4ejujo6IDjSSQS6PV6cQprrjxER0cDABYvXoyFCxfi48eP4j58+vQJV69exZEjR7B9+3YAQHFxMV6+fAmbzYasrCwxj2azGZGRkVi3bh1ycnLw4MGDAM8k+3/AhYIFpfT0dCQmJiI3NxexsbFQqVTIzs5GTEwMACAzMxNOpxM1NTV4/fo1BEGA2+3Gpk2bfOKsWrVK/D08PBwrV64UOwFOX2wnJibEdRITExEaGiq+TkpKQkNDA8bGxnzivnjxAgCwd+9en+UTExM+8QIVaLxly5b5fM8RaB5mIwgCJicnoVAofJYrlUrcuXNHfL106VKf5jaLFi3C169ff+0A2b8aFwoWlMLCwnD58mU8ffoUdrsd3d3duHLlCoxGIwoLC9HY2AibzQa1Wo309HQUFRWhubkZ796984nzfX9gAPjjD/9fy/35p++/xHSTlx+3m47b0tIyozXl94UmUIHG+/FvgeZhNj/bz2/fvvnkYbb1iB86/Z/ChYIFpfv376Ovrw86nQ5yuRw6nQ5Hjx5Fa2srCgsLcfHiRZSUlECr1YrbDA4OzrjQ/6r+/n5MTU2JheHx48eIiYlBVFSUz3rT00Ojo6NIS0sTlx87dgxxcXHIy8ubc6zve5zPN95ceZitj/q01atXIyQkBA6HQxwfABwOR1B0emPBg+96YkEpJCQE9fX1aG5uxtDQEJxOJ3p6esT7/6Ojo3Hv3j243W64XC6YTCY4nc55Tft8z+PxwGw2QxAE3LhxA83NzSgoKJixXmxsLLKzs1FZWYmuri54vV7U1NSgpaUFcXFxAY0VEREBj8eDkZGRecebKw8REREA/m4b++P0WVhYGLRaLWpra9He3g6PxwObzYbOzk6fwsMYFwoWlJKTk2E2m3Ht2jXs3LkTOp0OSqUShw4dAgBUV1djfHwce/bsgVarxYcPH2A0GvHq1St8+fJl3uMqFAp8/vwZarUap0+fhsFg+Om7+aqqKmRkZKCiogK7du2C3W6H1WpFampqQGPl5+fDbrcjJycHU1NT84o3Vx4iIyOh0Whw8uRJHD58eMb2paWl2LdvH8xms3hr7KlTp7Bjx47AEsb+E7jDHWP/KC8vx/DwMC5cuPC/3hXGggp/omCMMeYXFwrGGGN+8dQTY4wxv/gTBWOMMb+4UDDGGPOLCwVjjDG/uFAwxhjziwsFY4wxv7hQMMYY8+svAWqPkXeIYT0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(z_samples_SD[:,0])\n",
    "plt.xlabel('sample iteration');\n",
    "plt.ylabel('$z_0$')\n",
    "plt.title('Sampled values of r.v. from Gibbs Sampler');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2a: Compute the empirical MEAN of your samples\n",
    "\n",
    "You should estimate the mean of your samples, which is effectively:\n",
    "\n",
    "$$\n",
    "\\frac{1}{S} \\sum_{s=1}^S z^s \n",
    "$$\n",
    "\n",
    "If your sampler is correct, and you have enough samples, you should get *very close* to the true mean of your target distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO compute the mean of z_samples_SD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2b: Compute the empirical covariance of your samples\n",
    "\n",
    "You should estimate the covariance of your samples, which is effectively:\n",
    "\n",
    "$$\n",
    "\\frac{1}{S} \\sum_{s=1}^S (z^s - \\bar{z} )  (z^s - \\bar{z} )^T\n",
    "$$\n",
    "\n",
    "If your sampler is correct, and you have enough samples, you should get *very close* to the true covariance of your target distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO compute the covariance of z_samples_SD. Hint: use np.cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Compare to samples from built-in routines for 2D MVNormal sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw 1000 samples from target distribution using `np.random.multivariate_normal`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_mu_2 = np.zeros(2)\n",
    "true_cov_22 = np.asarray([[1.0, 0.8], [0.8, 2.0]])\n",
    "true_samples_SD = random_state.multivariate_normal(true_mu_2, true_cov_22, size=S-1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize results of your Gibbs-sampler and Numpy's built-in sampler side-by-side\n",
    "\n",
    "Do the results visually match?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAEiCAYAAAD3Z6a+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOyde1xUdf7/X+ecmUEtzaSS3Fb7dcG8gKDFNlsqhlmaFmVtFwtDVwTxgrah5teW1l0xsmTFGyIQs2G1pWneyjTJ0ilXxKTUTFO8FF2mSC1hLufz+4POp3NmzsAgw2WY9/Px4KFzzpkznzkz5z2vz/vzvgiMMQaCIAiCIIggR2zpARAEQRAEQbQGSBQRBEEQBEGARBFBEARBEAQAEkUEQRAEQRAASBQRBEEQBEEAIFFEEARBEAQBgERRm0GWZbz66qt4+OGHccsttyAyMhKjRo3CihUrUFNTw49bu3YtevfuzR/37NkT69ev93reWbNm4cknn2zKobc4d9xxB5YtW9bSwyCINsucOXPwxBNPtPQwGsXp06fRs2dP7N27FwBQVVWFN998k+9valuZk5ODO++80+fjv/nmG2zatKnJxtOaaei1UmPw81iIFsDpdGLixIk4ePAgUlNTYTabERISgrKyMmRnZ+Pjjz9GYWEhBEHAiBEjMGjQoJYeMkEQRECzcOFCVFRU4MEHHwRQK/xkWW6y1xs3bhzGjBnj8/HPPPMMunbtinvuuafJxtQWIVHUBigoKMAnn3yCNWvWoGfPnnz7Nddcg379+mH48OH44IMPEBsbi3bt2qFdu3YtOFqCIIjAx73ucceOHZv09S655BJccsklPh9PdZkvDlo+C3AYY1i9ejXi4+M1gkihe/fu2Lx5MwYPHgzAc/kMAI4ePYqHHnoIffv2xX333Yf9+/dr9judTjz77LOIjo7GbbfdhhUrVvAbzuVy4fnnn8fAgQPRt29fjBo1Clu2bPE63q+++grjxo1D//79MWDAAEyaNAmnT5/m+z/55BM8/vjjiI6O5uPZuXMn33/HHXeguLgYEydORGRkJOLi4vD+++9j69atGDZsGKKjozFhwgT8+OOP/Hy9e/fGO++8gzvuuAPR0dGYOHEivvnmG69j3LZtG+69915ERETg7rvvRn5+Pp8BKi70FStWwGw2Y/jw4bDb7V7PRRD+omfPnnjzzTcxZswYREZG4u6778brr7/O98+aNQszZ87E3LlzER0djdtvvx1Llizx+uPo7bvsyz1YUFCA5ORk9OvXD3FxcViyZInm3MXFxbjjjjvQr18//O1vf0N1dbVm/9dff43p06fDbDYjOjoakyZNwqlTpzSv0ZD73Nt7Kykp4ffygw8+iOPHjyMnJwe33norYmJi8M9//pM/R2/JxdsyTE5ODt58803s2bMHPXv2xOnTp/nyGWMMd9xxB3JycjTPycvLQ2xsLGRZRlVVFWbPno3bb78dffr0we23347nn3+e25mcnBw88cQTmDp1Kvr3749FixZ5jKWuz2nWrFmwWq146623+O+CLMtYsWIFhgwZgqioKIwePRoffPCB7vXzRn3fQQB48803MXLkSERGRuLOO+/EK6+8wvetXbsWd911FzIyMjBgwACkp6dj7dq1uPvuu7Fu3TrExsaiX79+mDZtGr799ls89dRTiIqKwuDBg/HWW2/x89R3/RoFIwKakydPsvDwcLZhwwafjl+zZg3r1asXfxweHs4iIyPZm2++yY4ePcqee+45FhERwb755hvGGGMzZ85k4eHh7Omnn2ZHjx5lmzdvZv3792erVq1ijDFmsVhYXFwcKysrY6dOnWJLlixhvXr1YidPntR9/fvvv5/NmTOHnThxgh06dIg98sgj7IknnmCMMfb111+ziIgI9tJLL7GTJ0+yQ4cOsYkTJ7Jbb72V1dTUMMYYGzJkCIuKimJvvfUWq6ioYCkpKSw6Opo99NBD7MCBA8xqtbJbbrmFZWVlMcYY+/jjj1l4eDgbOnQo27VrF/v888/Zww8/zEaMGMEcDgc/59KlSxljjJWUlLB+/fqxN954g1VUVLBt27axQYMGsZycHMYYY6dOnWLh4eFs5MiR7OjRo+yzzz5r0OdFEBdLeHg4+9Of/sQ2bdrEjh49yubMmcNuuukmdvr0acZY7b3ap08fNn36dPbll1+yDRs2sKioKLZixQrd8+l9l329ByMjI9lrr73GvvrqK7Zo0SIWHh7O/ve//zHGGHvrrbdYnz592Ouvv86OHTvGXnjhBRYeHs4ef/xxxhhj586dY4MGDWKJiYnss88+Y+Xl5ezxxx9ngwYNYmfPnuWv0ZD73Nt7GzZsGPvkk0/YoUOHWFxcHLvllltYeno6O3bsGFu9ejULDw9nJSUljDHGFi9ezIYOHao5j3qbcs7//e9/7Pz582zGjBns4YcfZt999x1zOp1s5syZbOzYsYwxxrKzs9mwYcM057r33nvZSy+9xBhjLCkpib+XkydPsldeeYX17NmTvffee/x1w8PDWVZWFjtx4gQ7efKkZiz1fU5nz55ljz32GJs2bRr77rvvGGOMZWVlsTvvvJPt3LmTnThxgv3nP/9hffv2ZR9//HG93z2F+r6DBQUFLDIykv33v/9lx48fZ6+++iqLiIhg+fn5jLHa35/w8HD21FNPsZMnT7KjR4+yNWvWsD59+rDExET2xRdfsJKSEta7d28WExPDXn75ZXbs2DE2Z84c1rdvX1ZVVeXz9XP/LH2FRFGAs2/fPhYeHs4++ugjzfZRo0axqKgo/jd37lzGmL4oeuGFF/hjl8vFhg4dym/emTNnskGDBnGDyBhjOTk5bODAgYwxxubNm8dGjhzJbzxZltnOnTu5cXOnf//+bOHChVyQnDp1ipWVlTHGGKuoqGCrVq1isizz461WKwsPD2dff/01Y6zWWM6YMYPv37FjBwsPD2dWq5VvmzZtGpswYQJj7HdRtH37dr6/oqKChYeHs507d/JzKqLokUceYQsWLNCMef369SwyMpK5XC5uGF999VXd90cQTYXyI6nw888/s/DwcLZ582bGWO29evvtt2vu1ezsbHbbbbdp7ikFve+yr/fgtGnTNOe6+eabWV5eHmOMsQcffJDNnj1bs3/06NFcFL3yyissKiqK/fTTT3y/zWZj/fr1Y6+88gp/jYbc597e2+uvv863LViwgPXt25dVV1fzbWazmY+7IaKIMcaeeeYZ/p4YYxpRpNgYZdJ05MgRFh4ezo4dO8YYY+w///kPO3LkiOa1YmNj2ZIlS/jr9uzZk124cEF3LL58TmPHjmUzZ85kjDF2/vx51rdvX7Zjxw7Na86ZM4eNGzdO9xrqUdd3UJZl9uc//5m9+OKLmudkZWWxW2+9lcmyzEWRch0Y+10oHT9+nG974IEH2JgxY/jjo0ePsvDwcPbpp58yxny7fhcriiimKMDp3LkzAODnn3/WbF+xYgUcDgcAYObMmXUu8URHR/P/i6KI3r1748svv+TbIiIiYDKZNI9zcnJw9uxZPPbYY3jvvfcwaNAg9O3bFwMHDsR9993ndX192rRpeP7557F69WrceuutiI2NxahRowDULvXFx8ejqKgIX3zxBSoqKnDo0CEAtct0Cj169OD/b9++PX+uQrt27VBVVaV53ZiYGP7/7t27o0uXLjhy5AgGDhyoOe7QoUMoLy/Ha6+9xrfJsozq6mqcOXMGgiAAAP74xz/qvj+CaEquvfZa/v9OnToBAL/PAaBfv36aezUqKgrLli3DTz/9hC5duuieU/1d9vUeVI8DqI2nUcbx5ZdfIj4+XrM/KioKX3zxBd9/3XXXcdsFAF26dMH111+PI0eO8G0Xc5+7oz6+Q4cOuOqqqxASEqI5R1Msf3fv3h39+/fHxo0b0adPH2zYsAH9+vXDddddBwB49NFHsX37drzxxhs4ceIEvvjiC1RWVmqWf6688kqv8Z++fk4Kx44dg91ux7Rp0yCKv0fNOBwOXHHFFQ16b96+gz/++CN++OEH9O/fX3P8LbfcglWrVsFmswEABEHANddco/ueFDp06KA5RvnMlM/Kl+t3sZAoCnC6d++OK664Anv37sWIESP49m7duvH/1xdYLUmS5jFjTGNY1TcRAP7FMxqNuO6667Bt2zZYrVbs2rULmzZtQkFBAXJzc/GnP/3J47USEhIwYsQI7NixA7t370ZmZiZWr16N119/HRUVFXjsscfQr18/mM1mjBgxAk6nE8nJyZpzGAyeX1v3Mbrj/hxZlnWfYzQa8de//pULNTVdu3bFd999BwAaw0oQzYX6vlRgqpgh9++58gNZ1/2h/i5/+eWXPt2DdY1DmTioMRqNuq+nRpZlzXEXc5+7436Ohj7f6XQ26Hg1999/P5YuXYqnn34aGzduxPjx4wHUXqekpCQcP34co0aNwn333YfIyEiMHTtW8/y67Lavn5OC8nnl5ORoxCbQ8Gvi7bP39rkq30HlsxBF0eMckiR5jMPbuHy9fhcLBVoHOJIkYcyYMVi7di2OHTvmsd9ut3sNRlQ4ePAg/7/D4UB5eTluuOEGvu3w4cMaw7tv3z5cc801aN++PYqLi7F161YMGjQIs2fPxpYtW3DNNdfoBlv/9NNPmDdvHpxOJx566CEsWrQIL7/8Mg4ePIjDhw9j7dq1uPrqq7Fq1SqMHz8eAwcOxLfffgug8ZkUn332Gf//8ePHUVVVhV69enkcd8MNN+DEiRPo0aMH/zty5AgWLVrUqNcniObg0KFDmtnyp59+im7dumm8MnXhj3vwpptuwr59+zTb1PffDTfcgK+++krj5fnxxx9x/PhxXH/99T69RlNgNBrxyy+/aLZVVFR4PV5P/KkZPnw4qqqqUFxcjO+++46nxh89ehQfffQRcnJyMH36dNxzzz24/PLL8f333/t8jX35nNTj69GjB4xGI7799luNbduwYQPWrl3r02vWx6WXXoqwsDCUlpZqtpeWluLKK6/EZZdd5pfX8cf1qwsSRW2ApKQkmM1mPProoygsLMSXX36JU6dOYcOGDRg9ejS++uorDBgwwOvzV61ahY0bN+LYsWOYO3cuzp8/j8cee4zvP3XqFP7+97/j6NGjWLduHf7zn//wGYkidHbs2IEzZ85g+/btOH36NPr16+fxOpdddhl27tyJZ599FocPH0ZFRQXWrl2LTp064f/9v/+HLl264MyZM9i1axfOnDmD9evXczHSWBf3c889h3379qG8vBzp6emIiIjQLKkppKSkYNOmTVi5ciVOnDiBkpISPPvss2jXrp3uDIkgWhMnTpzA/Pnz8dVXX2H9+vWwWCzcQwHUio9z5855fb4/7sHx48djy5YtKCoqwvHjx7Fs2TLND+W9996LLl26YMaMGTh48CA+//xzzJgxA506dWpUTZ363lt9REVFwWaz4eWXX8bp06exevVqTdadO5dccgm+/fZbnDp1Stej1LFjR8TFxWHRokUYMmQIF6adOnWCwWDAli1bcPr0aZSVlWHSpEmw2+0+X2NfPqdLLrkEp0+fxpkzZ9C+fXs8+eSTePHFF7F582acOnUKFosFS5cu9WsoQEpKCiwWC9544w1UVFTgv//9L1555RU8+eST9YpIX/HH9asLWj5rAxgMBixbtgzr16/H2rVrsWLFCvz666/o1q0bbr/9duTk5HjEAKiZNGkS8vLycOzYMfTp0wf5+fma+IM777wTNTU1GD16NC677DJMnToVDz30EAAgOTkZ1dXVeO655/DDDz/g6quvxpQpU3D//fd7vI4oisjNzcWCBQvwxBNPwG63IyIiAvn5+ejYsSMSEhJw7NgxTJ8+HS6XC9dffz2ee+45zJ49G+Xl5Y2aRcbHxyMtLQ3nz5/HkCFDMGfOHF337KBBg5CVlYWVK1di8eLF6NKlC+Lj4zF9+vSLfm2CaC769++PX3/9FQ888AC6dOmC6dOn4/HHH+f7H3zwQcTExGDBggW6z/fHPTh06FBkZmZi2bJlWLhwIf785z/jL3/5C/dkh4SEID8/HwsWLMCYMWMgSRLMZjOKi4t5jMrFUN97q49bb70VU6ZMQV5eHhYtWoRBgwZh6tSpKC4u1j3+gQcewLZt2zBixAivx8THx2PTpk247777+LauXbti/vz5yMnJQVFREbp27Yrhw4eja9euKC8v92msvnxOY8aMwd/+9jeMGDEC27ZtQ1paGoxGI7KysvDDDz/gj3/8I/7xj3/ggQceaPjF8sIjjzyC6upq5Obm4rnnnsMf//hHzJo1SzPJbiz+uH51ITB/+JsIopXyySefICEhAR988AHCwsJaejgE0WTMmjULlZWVePnll1t6KAQRsNDyGUEQBEEQBFrh8llVVRWef/557Nq1CzU1NRg8eDD+/ve/N6i8OUEQBEEQvvGPf/xDUzFaj7fffjsoSpG0uuWzhx9+GN26dcO8efPgcDgwefJk3HjjjcjIyGjpoREEQRBEm8OXIPVu3bppSia0VVqVKDpw4ACeeOIJWK1WdOjQAQDw/fff45dffqkzUJggCIIgCKKxtKrlMyVq/pVXXsHrr78Ou92Ou+66C0899ZRPz5dlGb/88guMRqPf0v8IgvgdxhgcDgcuueSSBhd9I+qG7BdBND312bBWJYp+/vlnHD58GP3798eGDRtQVVWFKVOm4Pnnn/dp+eyXX37RlIknCKJpCA8P99rKhbg4yH4RRPPhzYa1KlEUEhICURQxc+ZMGI1GdOjQASkpKfi///s/n0SRst4ZHh7e6EJ7n332Gfr27duoczQ3gThmIDDHHaxjttvtOHLkSFDEFjQ3/rRfQPB+R5ubYB7zsmXL8NJLL8HlckGSJMyYMQOTJk3ywwj1aQ4b1qpE0Q033ADGGOx2Ox9wQ3rPKC5nk8nkl95UgdjfKhDHDATmuIN5zLS843/8bb+A4P6ONifBOuZbb70VNpsNdrsdJpMJt956a5Nfi6a2Ya0qKMBsNuPqq6/Gs88+iwsXLuDbb79Fbm4u7r333pYeGkEQBEEQKsxmM7Zv34558+Zh+/btMJvNLT2kRtOqPEUmkwmvvPIK5s+fj7i4OLhcLtxzzz0+B1oTBEEQBNF8mM3mNiGGFFqVKAKAsLAwLF68uKWHQRAEQRBEkNGqls8IgiAIgiBaChJFBEEQBEEQIFFEEARBEAQBgEQRQRAEQRAEABJFBEEQBEEQAEgUEQRBEARBACBRRBAEQRAEAYBEEUEQBEEQBAASRQRBEARBEABIFBEEQRAEQQAgUUQQBEEQBAGARBFBEARBEAQAEkUEQRAEQRAASBQRBEEQBEEAIFFEEARBEAQBgEQRQRAEQRAEABJFBEEQBNGmsFqtyMzMhNVqbemhBByGlh4AQRAEQRD+wWq1Ii4uDna7HSaTCdu3b4fZbG7pYQUM5CkiCIIgiDZCSUkJ7HY7XC4X7HY7SkpKWnpIAQWJIoIgCIJoI8TGxsJkMkGSJJhMJsTGxrb0kAIKWj4jCIIgiDaC2WzG9u3bUVJSgtjYWFo6ayAkigiCIAgiwLBarV6Fj9lsJjF0kZAoIgiCIIgAgoKpmw6KKSIIgiCIAEIdTF1TU4OMjAxKv/cTJIoIgiAIIoBQgqlFUYQsy9i2bRvi4uJIGPkBEkUEQRAEEUAowdRDhw7lwojS7/1DqxVFLpcLTzzxBGbNmtXSQyEIgiCIBtOUlaXNZjMyMjIQEhJC6fd+pNUGWi9ZsgR79+7FH/7wh5YeCkEQBEE0iMYGQ9eVXaZA6ff+p1WKIqvVinfeeQd33313Sw+FIAiCIDi+iBVAv7K0r6KlIYKK0u/9S6sTRTabDXPmzMGSJUtQXFwMl8vV0kMiCIIgCF2xYjKZdI9VgqGVY0NDQ5GZmelVTKnFVmMEFdE4WpUokmUZTz/9NJ588kn07t37os/z2Wef+WU8paWlfjlPcxKIYwYCc9w0ZqIp8Jf9AgLz827NYy4uLkZNTQ0PbC4uLkZiYqLumE0mE5YuXYrS0lJcdtllmDp1KhwOB4xGI5YvX47IyEh+7IEDB5CSkgK73Q5JkvDYY4/BYKj9eTYYDAgLC/N6XQ4cOIDS0lIMGDBAc876aM3XuS6aetytShTl5ubCZDIhISGhUefp27cvQkJCGnUO5UsWSATimIHAHHewjrmmpsavP9qEJ/6wX0DwfkebErvdjoKCAu79GTNmDAB4HfOAAQOQmJiIzMxMOJ1OyLIMp9OJyspKJCYm8uO2bt0Kh8MBxhicTidWr16NpUuXwmaz1blMZ7VakZqaipqaGkiShCVLliApKane99Har7M3msOGtSpRtH79enz33Xe4+eabAQDV1dUAgM8//xwbNmxoyaERBEEQQY5eYLMvngv3pTT3LLHY2FhIkgRZlgEATqcT+fn5yM7OBgDdZTer1YqMjAxUV1eDMQZZlpGamoqIiIgGL7X5GicVDLQqUfTOO+9oHs+ZMwculwsLFixooRERBEEQxO9cTGBzfVliZrMZS5YswaRJk3gc7Z49ezB48GAIggCXy6UJuFZimxRBpOByuRocf0QtQ7S0KlFEEARBEG0RPTGl9tAkJSWhrKwMK1as4PsdDgcEQQBjTBNwrQRiqwURAEiS5LVWkfq11MHhFNStpVWLon/9618tPQSCIAiC8ImGLENZrVbExsby4OuSkhIkJCTwmCUAMBqNGk+RksEWGhoKk8nEg76BWkG0dOlSr5ltam/Q0qVLeWxOfUt7wUarFkUEQRAEEQg0dBnKYrFw8WO322GxWLB8+XKUlJTAYrEAAE86slgsqKysRGpqKlwuF4xGI3JyclBWVoZVq1bB5XJBkiRERETovpa7N6i0tJQHelMBSC0kigiCIAiikXhbhmpoELP7MpvVakVRUZEmfshut6OsrAzdu3cHYwyMsTrjidy9Qe4ZXFQA8ndIFBGtBsqAIAgiUNFbhtJbIlNsW0JCAvLz8/k+b6VovMUPKa+pZK3pxROpbaraG+St4CRBoohoBRw4cACrVq1CYWEhnE4nZUAQBBFw6C1DKQUZgd+XyHyJNVKfQxExak9RSEgIF1GCIGj+VZ/HfTlv9uzZAAK3cGNzQKKIaFGsVis3HGrXcLBnQBAEEXi4L0NVVlZ6PdZiscDhcACozTJT4ojchUx5eTm6du2KEydO8OdeffXVKC8vh81mg9Pp5EUf1XaTssouDrGlB0AENyUlJbySK1A726EMCIIgAh2r1YotW7bwx5Ik8e3eUAuZ6upqzJo1CxMnTtQIIgA4ceIEJk6ciKqqKphMJoiiCFEUERoayo+JjY2FwWCAIAgwGAxkU32ERBHRosTGxsJoNEKSJJhMJkycOJGWzgiCCFhWrlyJu+66C1lZWXA6nQB+X9rKy8tDXFwcrFYrEhISYDKZIAgCJElCdHQ0jxECAMYYPvroozpfa//+/cjOzoYoinC5XEhLS9OILmWyqRePROhDy2eEX7jYIGmz2Yzly5ejsrKSAqwJgghoVq5ciYkTJ/LHRqMRALhoURrJlpSUYPbs2XjooYdQXFwMl8uFqVOnYseOHRg3bhxyc3N56w41HTt2xLlz5/jjqKgo2Gw2yLIMWZZRU1MDi8WCkpISnDx5Ei6Xq97MNEILiSKi0TS2THxkZKSmOSJBEEQgsmbNGs3j6OhoxMfHIzQ0FGlpaaipqYEgCKiqqsL999+PdevW8WMVQdOpUyddz86wYcPw/vvv88eCIKBz584IDQ3l4kmWZeTl5UGWZRgMBu51opAE36HlM6LR6AX0EQRBBBujR4/WPFaESEREBKZMmQKgtj9ZVlaWRhApVFZWYtGiRbrn/uijj/hyHFDrfYqNjYXNZtNknineIYfDgREjRmDevHkUktAAyFNENJqmLBNPtYsIgmgq3O1LQ+2N1WpFcXEx7HY7zGYzkpKSANR6jKKiopCdnc0zzIC6Y3sUr45a+Kj59ddfNY8FQcC6deuwceNGr+cNCwvjafiEb5AoIhpNY8rEuxsV933UvZkgiKbA3b5kZ2cjLS3NZ3ujPL+mpgYFBQXYvn07AMBms2H06NHIz8/nNYr0EEURkiThhhtuQM+ePTF8+HBMnTrV56Bop9OJrKysOs9fWVkJq9VKdrMBkCgivNKQWdPFlIlXG5W8vDyMGzcOCQkJVGeDIIgmx92+rFmzpkFtOpTnK8HTFosFRUVFmiateoiiiEcffRTff/89Ro8ejaSkJFitVmRkZGi8SvUhCIKugFKW0mRZxrp167B582aynQ2ARBGhiy9emsYubbkbldzcXBQVFfHXou7NBEE0Fe72ZfTo0fjwww892nR4s4PuzwfA7VlddOvWDa+99hoYY/jwww8BgAdh1/dcwLsYUnDf52slbaIWEkWELvV5afyxtOVevp4xpnktf3VvprgkgiDc0bMvERERmseZmZle7aDy/OLiYowZMwYAsGrVqnpf9/Tp0/z/NTU1WLNmTZ2CSC2ClD5nQMNqD+Xn5wOAxhNP6EOiiNClPi+Nt4yzhnaD3r59O1588UVs3LiR9z1Tv1ZjuzdTXBJBEN5wty/uj701eVXbOaXr/MyZM70GSddHXR4itfiJi4vDtm3bGlyM0eFwaDzximeLJoyekCgidKnPS+NuLEJDQy9KfJjNZsyePRtPPfVUnTfnxd68FJdEEATg3YbUZVvc7SDg2ZvMZDLBarVi4cKFuq9b13KXLMvYunWrz+9h27ZtHgJKFEV+LvftkiTx3mhqT/ywYcNowugFEkUEAH3DoOelUR+nNha+LLfVJWrq8gi537xTpkzB/v37eZBiXVBcEkEQ3gRAfcLA3W7pLacNGzYMJSUlHsInJiYGe/bs8WuLDbXwufbaazF79myUlZXh4MGD2LlzJ983aNAgLFiwAACQlZWFDRs2gDGmsYEWi4WHLvgyYQwWrxKJIsLnGYP6OEmSNNli5eXlEEXR48ZryPndX0tPcNXU1PA0VGWGVZcw8ldcEkEQgYu3SVtdkzk9u6WeZBkMBpw8eRIHDhxAbGws2rVrh+rqagiCgEcffRSffPJJk76nEydOYOfOnXjzzTdRU1PDt4uiiLvvvhtArfDZsmULGGOQJAnZ2dkwm80oLCxEYX5HFyUAACAASURBVGGhRrCpm8m6E0xeJapoTXiND3JHmVkox+Xm5iIuLg4rV65EWloaXC4XRFHkN15Dz6+g3IBz587FkCFDsGfPHk3JejXuZfWV52dmZvLGiMoSna83sfvzCYIIbBQxozSeViZt3rYDnh3rlQyu7du3Y8KECWCMIS8vDykpKQCA7du341//+heWL1+OtWvX4ujRo03+voqLizWCCABCQkJ4OENubi4P4pZlGTabDQBQWlqqiX+SZdmjmaya+mx4W7KZ5CkifFpislqtKCgo0MwsFLerUt9DlmUIgsBvvIacX/06GRkZ3K3rcrmwbt06GAwGjBo1Cj/++KPGTVxdXY2UlBTusVq5ciUmT54Ml8uFkJCQBs9ogmlGRBDBgjePcV2e5NjfOtYrbTPy8vKwb98+jB8/Ht27d4fT6YQsy2CMwWKx8Oft2bMH1dXVzf4eFe666y6UlZWhpqaG22tBEDS2d8CAAbqZv0ozWb1r4c2GtzmbydoQ1dXVbO/evay6urrR59q7d68fRtS81DXm3bt3s/nz57Pc3Fw2f/58tnv3bt397tuVfcOGDWOiKDIADACTJIlJksTat2/PcnNzWfv27fljb+dwP7+yraCggO3evZslJyezkJAQJggCfx31n/KaRqOR9erVi0mSxPeFhISw9PR0zRhFUWTz589v0DWcP38+P68kSV6f39a+H77iz3uM0OLvaxus39HGorZVycnJuvYoJiZG81htd1rDn8Fg0DyOj4/X2N69e/fy92cymZgkScxkMrGQkBCvdtzbb4SvNtMfNIcNI09REKCuHC3LMkRR9PCieAt0dn+uIAiQJAkzZsxA586dvdb3cA/Kcz+/1WrFkCFDeHySKIpwOBx1BiW6XC7+/z/+8Y84fPgwf2y327Fw4UJNIKIkSQ0OrKbAbIIIPhR7pXSzV7f+aNeuHS5cuKA5fs+ePZrHvhRd9BeCIKBjx444e/as12PcG8fGxMToZtaZzWYkJCSgpKQEJ0+eRF5entf4Kl8zkQPdZpIoCgKU9WDlxlUqSCvrwu5fdr0gZ+W5giBAlmXk5OR4FVVWqxWxsbFwOBwwGo26WQ0Wi4Wvhfta20OJKZIkCR06dIDBYOBl8UVR1BgmQRCwZMmSBrtxW0NgdrBkeRBEa0C9/KPYN8VG2mw2ZGdnIzs7G4cOHWrpoQKoDVuoSxCpEQQBBoPBq1BR2xoAKCoq8hA39S2PtQab6U9IFAUBipJXe4q81RYC4NEk0WAwcC+NIjy8pXBarVY+01KOy8rKQlhYGADfK6oaDAb07dsXn376KRhjCAkJweLFi1FWVobCwkJs2LABkiQhPj4eYWFh6NSpExYuXMg9TaIoIiIiQvfcjSkP0JDzXAxtbn2eIFo5JSUlGk+4wWDg/+7ZswebNm2C0+mst71GU9DY12SMQRAElJeXc0+YzWZDWFgY7Ha7h63REzd1ZegpNLbIbmuCRFEQoFbyyk3hrbYQAM02m82GxMRE5Obm6gbtqd3OZWVlKCgo8Ghq+Pbbb3MxVVhYiB07diA6OrrOMTudTuzfv58/njZtGpKSkpCZmQmHw8EDHGNiYhAbG4vBgwdrPEUul0u334+/RMeBAweQmprqd/HiiwEiCOJ3Gjs5qaqq0rTOmD59Os6ePYuCggKsX7++TlFiMBguuoq1L5hMJvTq1UtjCxuK3W5Hamoq94CJogij0YjExEQPW6OXpdvWlsfqo9WJogMHDiArKwuHDx9G+/btcdddd+Gpp55C+/btW3poAY03Ja/3ZdfbprhVRVFEdHQ0YmNjYbFYkJ+fzyum6s1q3Je11OKrIbOgtWvXIj4+XmPAZFlGaGgosrKy6uwurTaaWVlZPD6gMaKjtLS0ScRLsBkggmgMjZ3kWK1WLFq0iD8WBAGdO3dG586dedaZmhtvvBFffvklf9yUggiobc9x5syZRp2D/ZbFq7wXWZb5uJVeanXFX7a15bH6aFWi6Ny5c5gwYQJSUlLw8ssv44cffkBqaioWLlyIuXPntvTw2hx6X3ar1YqxY8eisrKSL3kB4Ns2b96MvXv3egQaAtoePUohR/cARGV9u7y8vEFu4aNHj3p4g0RRRFlZGTZs2OBxvNFoREJCgiagWxRFTbC2KIoIDQ1FZmZmg292JaXV3+Il2AwQQTSGxnpWS0pKPGyCenKobtRqNBrRp08fjSjSw5/LbLIs4/vvv2/0eZQQCGWJUBAEdOrUCYIgAAD/1xttaXmsPlqVKDpz5gxuvvlmPPnkkwCAsLAw3HfffXjzzTdbdmBtDL3MMGW7e5aaUvVUKcyo3Fh6qI2Bt2PGjx/PK8mqvUiCIHiIFnfcvUHKjaw2QIIg4L777kN6ejrMZjPuv/9+HtDtfu4bb7xRk2nSkErbYWFhTSZegskAEURjaKxnNTY2FkajkdsIURRRXl7OA6y3bNnCl9CcTifWrVtX5/kEQdDNVmtJ2rdvj+zsbNhsNlRVVeGll16C0+nkMZjKe6Ol+lpalSi66aabsHTpUv6YMYb33nsPffv2bcFRtS3qcjd7y1JTCx1JkviNpMZoNCI6Ohp79+6tUzQpsUTuxuiSSy7Br7/+Wu/41UJKlmUcPHgQRqMRDocDoihi6dKlvO2H1WrV9SIphIeH44svvvB5lqm+dgaDATt27MDs2bPrHTNBEE2De7yksjTvSyq5ctzw4cO52HG5XJg8eTK3dUpoAACfvD+MsRYVRNdccw1Onz7NH8fHx/MJIgCkpKR4JM0oiTe+CMpgyIxtVaJIjdPpREZGBk6dOuW1+zDRcOpyN8fGxnq4UdXeG8YYZsyYgbNnzyI/P58LkXvvvRfp6ekoLy/Hvn37vIoixhgmT56M/Px83HjjjTxDDQDOnz/vdcyKB8pgMGDkyJHcgDHGsHPnTgiCgF69evFgbPV7dTdkkiShZ8+eCA8Px/Dhw7F58+Z619T1rp3yuK0aBoIIFJR7sL5MWj1PsNVqxaZNmzTbFG+4kswRSKgFkSAICAsL0whE935ngiBg6NChyMjI8MlLHgyZsa1SFNlsNsyYMQM//vgjiouL0bVr1wY9/7PPPvPLOEpLS/1ynuakvjGHhYXBYKj92A0GA8LCwvhzDh8+7CFo1EtOgiDg2LFjCAsLw6OPPorVq1fD6XRi/fr1qKmpwfvvvw+Xy8VL4+vhcDiwZ88e3ZgkPURR257v4MGDHscwxnDw4EFMnjwZRqMRkZGR/L0ajUaNt4sxhiNHjuDQoUN4++23NUtwhw8fhslk8jqWuq5doBBo4w1G/GW/gMD8vC9mzEoPMMW7XVxcDAAe20wmEw4cOIDS0lJcdtllWL9+vWZZXpZlGI1GzXJ+cxZm9Cfst9YkXbp0wQMPPIDi4mLNe1WyiB955BGYTKZ6r7veNa7LXjYVTf6dbnTNbD9z+PBhNmjQIDZlyhR2/vz5Bj2X2nz4Nua6yrWry9W7l7cXBIEZjUYmCIKmxYb7n7c2HU39JwgCb2UybNgwlpubq2kfIkkSMxgMuuPztTy9ujVJoEFtPlo31Obj4se8e/duj1ZDdW2rqy1HcnIytyPx8fEt3rKjsX8Gg8HjehiNRhYfH8+Sk5N12zL5eo2bm6Br8/HDDz9g3LhxGDlyJMVq+Ih6jddX1a4E8iqdjZX14djYWISEhKC6uhqCIOCqq65CZWUlf95ll12GqqoqAJ5By2pYC7mcDQYD3nnnHd4wduvWrcjNzcXy5ct5KfvQ0FCkpqZqUmndmyUCnmvn6sezZ88OyBk4QbRVvGVtum/LzMzUxE26YzAYkJCQgPLycqSmpvIEk0D1FgG1tlqpQaRcj/Pnz+PFF1+E3W7ntePqWwoLmszYRssuP5KTk8PCw8NZv379WFRUFP8bMWKET89v654idw+Pu3JviPdC/VyTycRnDLm5ucxoNDJRFJnRaNTMOKKiolp81qP3JwgCEwRBd/YXExPjcf3S09P5e5QkyaNZovt11Wt42xq/H/VBnqLWDXmKmn7M9XmKQkJCuA1sabvmz7/c3FzNdYiNjdXsT05ObrLr7a3R+MUQdJ6iyZMnY/LkyS09jFaJXpCbe9B0aWkpEhMTdZ9rsVgA/N5mQ/1cl8uFFStWoLCwEImJiTzIENBmex04cIBnerUmWB3ZIaWlpUhJSUF0dDRPv5ckCWazGbt27QJjDO+++y6GDx/OZ0Du1zU/Px/V1dVgjKGmpgYlJSUYNmxYUGRiEEQg4u3eVHs7qqqqsGjRIk2GWU1NDf7v//6vyYsyNgV11Udas2YNgNp43dDQUHz44Ycex/jbngVqYHarEkWEd/SyxmJjYzUVSQcMGODxPKvVisGDB3MhU1BQwJeRFMGjNgiVlZW87gegXSaTZRm33347qqurfQ6UbmlcLhdyc3M1xctcLhdfYgNq33dqaiqvyh0eHs4DvCVJwr59+zRlCUJDQ5uszQdBEI3DlwamypL42bNncfDgQY098KVYYkhICC8n0hqQJAlhYWFeq1+/99572Lp1Kw8gVy8HSpKE6OhovwuYxhbWbCnE+g8hWgNKzJAkSZr4l/oqklosFo1nx263w2KxIC0tTTcuSClKOGHCBN1zVlRUoH///j6NuSUyE/RgvxUnUyq5qhEEgQtLl8sFp9OJgwcPwuFwYNSoURg3bpxHBp7NZtNt8wGAx2lZrVbN/wmCaB70fozdUYRTbm6uRhD5SmsSRIB+5WuDwYCYmBjeXQCAphCv0gNt2bJlsNls9V6zhuLtN6u1Q56iAEEvyC0zM5O7fp1Op9flMzWiKGLfvn2artDKTaO0xlCW1/SCCysrK9GpUyefxqyuQ9TSKEZB7WIWRRE9e/bEqFGjkJOT41F07ddff0VCQgIKCwu5ETQajYiNjeXp++pKuuoZqiRJEAQBTqeTPEkE0YzUVeVaWSI6efIkampqAq4OUV1ERUVpPPgPP/wwUlNTERcXx5f/FQYOHIhhw4ZpEkn83bYoUAOzSRQFEO7tH9xvfr3lM+VH3W63c6+Iuuo0+60o4vjx47kgslqteOedd3RFUU1NTcAW03T3jMmyjEOHDuHIkSMYNWoUfvzxR82scfTo0TCbzdixY4dHTJY6rss9s0Vd24QxFlCuY4IIdLz9GKsnLWrvSVuAMabJiBVFEampqTCbzcjOzkZ+fj727dsHl8sFk8mEhIQEzQRauWaKnfMXgdiyiERRAON+8+stVyk/6srsKC8vT7dAY/fu3bkgUscg6RHI6al6uFwurF+/Hu3atUN6ejr279+PqKgo2Gw2WK1Wrzd2XSLV3VMUKK5jggg0Vq5ciTVr1iAqKgqdO3f26Omo4J5c0tZwf0/KEpg6weSvf/0rEhISvIY2FBUVwW63o6ioKGi92ySKAhz1za/MFNRZBAC4+o+OjobJZPJwpQqCoDm2tWWXNQeKR6dz587IyMjQbRmg9hZ5E6BqkQog4FzHBBFIrFy5EhMnTgQAHkhsMpmQmJjIvboKyqTF3f61RZSEELUQlGUZ27ZtQ3R0NIxGI7Zu3aqxTYEaGO1vSBS1MdzjWhhjXOSEhIRg8eLF+Pe//61pl3HbbbdxL9GqVas8zhkTExMw2Wa+4p6BIYoiDAYDTp48CYvFwo1DdXU1Zs2ahV27dvGZWGFhIZYvX667XOk+Qw1Go0IQzYWSaq6glM3Izc1FUVER7w6v/PgrS0RK78a2iiiK/H2rheDRo0cxceJESJIEAJpJ38mTJ3kbo2D2blP2WRvDXe27Z57ZbDYMGjRI85wuXbogMzMTs2bN0q3PsXfv3iYfd3MiiiLPRJMkCenp6UhKSuK9ggoKCnhKPvut6azaNa3UhCIIomUZPXq07nZFHE2aNAnPPPMMBg4ciJUrVwIADxVo6ygT2e3bt+P666/X7FOWEJVs5Li4OOTl5YExhgkTJgTt0hlAnqI2h1K7SG/N3GAwcPVfUFAAh8MBQRCwadMmvP32215jhdpaDBEAXp/J5XLh7Nmz6N69uybWYMCAAfjf//6n62b3FtROEETzkpSUBAA8pujs2bPIz8/nkzvlfna5XEhJSYHRaNQUa2yryLKMdevWYfPmzSgpKcHTTz/NlxkVlKVGAHwiDQSPaPQGeYraGGazGePGjeP1eBSvCFB7o5SXl8NsNiMtLQ2CIECWZTgcjjYpfLzRtWtXjVHMy8tDaGiopqbG+PHj0a5dO+4xUrxK8fHx2LFjBwBQDSKC8CNWqxWFhYVe7ylvdb+SkpLw7rvv4vnnn0dCQoLXmm1Kd3d1dmhbR/EEJSUlYcyYMZp99913H7Zv385jJJXfitDQ0BYabeuAPEUBjLey7AkJCTyLQBRFzaxp0qRJAICXXnrpog2DIAia5rCBhvss0eVyYcuWLZqU1IiICB40HRoaqolLsFqtSElJoRpEBOEnlFjImpoaFBQUeNxT9VWpVmzhnj174HA4eHV6deyg0WgEgDYdS6RHYWEhEhIS0KdPH16nTRAExMTE8GuYnZ2NyZMnw+VyIS0tDREREbrXNxiSRkgUNTP++nLpGQnFFarOggoNDcWkSZM0buT58+c3KiWVMRawgggAvv32W49tb7/9NsLDw7Fq1Sq4XC6sWrUKS5cuxezZszXHWa1WpKWl8WKOykwsWAwGQTQFSiyk4s1xz3yqKzNKLajUEz1BEPDUU0/h7NmzqKysBABs3Lixed9YM9GrVy+MGjWKLx+qhZ/D4eC/BerK1lVVVcjMzERsbCxsNhvveent+gZLSyMSRc2IP79cekZi2LBhfL86C2rLli1Yt24d31dRUdG4NxLg6MUTyLKMF154ge9zOp2YPHmyZsZktVoxZMgQTYl/xhjy8/Mhy3JQGAyCaArqqkJd3/6srCyPavRA7QRw0aJFmD59OoqKitp0Kv7hw4dx4sQJZGdn45577sH69es9+jXabDZNRf9FixZxu5Wdne31+gZbqj7FFDUjvvTk8ZWG9JVJT0/nKZiEd9wNptPp1FR4VT4/NUpMlj97BhFEsKF4t5OTk3UnFsr+efPm8f1WqxX333+/ZsLnjsPhwMKFC9tcSw931Nl269at83ivNptN4ykCau2bYrdsNpvH9VXQ+61py30dyVPUjNQ3G2oIeoUClXYeas+Gsn/ZsmWaZTRCi9Ig0eVyaVzMynq82Wzmn59eM0glkyNYa3sQRGNR2uf4ktmpeN31PETuyLLM729JktpsTFFdMaLvvPMOevfuzWOs1I2wleDquir3u//WqFc83GtBBTysDVFdXc327t3LqqurG32uvXv3+mFEnuzevZvNnz+f7d6922/n2b17N2vfvj0TBIEZDAaWm5vLt0mSxNq3b892797NkpOTGQD68/JnNBpZcnIyi4+PZ4IgMABMkiQ2f/58zXVPTk5msbGxLCQkhEmSxEwmE0tOTm70Z9rU+OM77c97jNDi72vbVDasKfE2Znd7lpyczCRJatD9LYoi69evX4vbmZb6M5lM3GaFhISw9PR0ZjQamSiK/DfCF+bPn8+vvSiKzGAwaH5nmpLmsGHkKWpm/NEgzz02aezYsdw9rMTC3HPPPXwNXVnWSUhIQF5eHnmLvKDMINPT0/Huu+/y6xsaGoqUlBQAta1SunfvjpiYGMyfP58CrAnCjxw4cMCj/QTgGXoAgHtt3T0k7tXqFWRZxqefftq0b6AV43K5MGHCBHTv3h1hYWGorKz0GlxdF+oVD8W73tBztGZIFLVi3DPVlMcnT570MBBqQ+B0OrFhwwa+DKQUbTSbzRg1alSda/DBTn5+PiorK3HXXXchLCwM0dHRmDJliiaWSBRFGI1G7NixwyM7jSCIi6OuUhfuoQcJCQlISEhARkYGtm3bphFBypIQa8MxRPWRnp4OAFi9ejW++eYbAODXzWw2o7Cw8KLberhnNysNZ9tK+ACJolaKuzcoOzubf/kMBoOmd01CQgI6deqErKwsAOCVmoHaWJfExEQuqoi6cTgcXDQajUZER0d7xCDIsgyn0+mRtqonYMmLRBC+UVJSwgvJunsd3ONalHvsuuuug9FohN1uD2oRBNRO1nr27MnrDKl7YI4bN44LIrX4lCQJEyZM8GieWx/qFY+IiIg2ZetIFLVS3N3Fa9as0ZRiV9ygyhexpKREd3YkiiKio6O5yKqurm6JtxOQOBwO3Ua4SvNYZVZUl4ClNH2C8I3Y2FjehkPP66D+IXZvfH3LLbd43KveltHaKrIs48SJE1ykeGvdoRaf7vsuBn+EhLQmSBS1UtzdxaNHj8aHH36ocR/rpU26z5iUCqXquCOi4fTq1Qs9e/YEAISFhaFLly48/b4uAdtW1tkJoqkxm81Yvnw5KisrNV4HPa+r+49+//798emnn/LMUKPRiLS0NE3l/mAQSIq9CQ0NhSiKYIx5CMz6xGewQ6KolaLnLq7PTTly5Ei4XC5s3LhR0+1eMRRKCibRcK688koefG0wGHiAol7hM3cBS0aHIHwjMjISiYmJ/LHVakVsbCwcDgeMRiOfYHiLMVLqikVHR8Nms2Hp0qUoKytDXl5eS72lZkUURezZswebN2+Gy+WCKIp44IEHkJaWhnbt2qF3795ISEjQFZ/EbzQ6v60VEQgp+U2Bkq6qpFaqU8rxW1r57t27WXp6eounhQbqnyiKTBRFBoAJgsCvryAIPB1fXWrBX6UX/Aml5LduKCXfc8zuZUSSk5P5Pm/3WHp6OhNFkQmCwO1hS9uP5vxT237FZqn/QkJCWEFBQbN8nv6GUvIJn1BcyUqAYlhYmKbIIGNMU5mZaDiyLEOSJAiCwD1FSqHHgoICJCQkaDLR2to6O0G0JrwlNlRVVfGEEwCorq7G22+/3YIjbX7YbyES3jLw7HY7SktLNR454ndIFLVifM1gcnclR0dH44MPPsChQ4cA1P6gr1ixApIkBV3wYVMgCAIGDhyIDz74gGf6UdwQQfifhIQE3uDUaDQiISHBI7FhypQpmtghd4LR3hmNRowfPx4HDx7Ezp07Nft8rRoerJAoaqU0pHmsEn9UXFyMyMhIj7o6ClS0sXHIssyLYVZUVECSJLhcLkiSpIkbonR8gvAfgiDwP0AbZF1TU4OFCxd6FT533nkntm7d2pzDbXEEQcD48eOxfPlyrFy5Eh999BFvdXLvvfciPT0dJpOppYfZavFbQ9i5c+f661T44YcfkJycjAEDBuC2227DkiVL/HbuQME9o6m+RqNmsxmJiYkoKyvTFURE3fTo0QNRUVF1HiOKv98ux48fh9PpBGOMG2vgdzE7d+5cxMXFUW0oguNPGxkslJSU8GVqxSOrblCqZFipMRgM6NGjB3r37o0OHTpo7s9gQAk8t1qtmDp1Kl/6X758Od566y2aqNWD30TRhx9+6K9TYfr06bj88suxa9cuvPzyy1izZg02btzot/MHAnqdib2hdCw+cOCAx76YmBiMGTOmCUfaNqioqMD+/fu97hdF0avLWSnkCDRczBLBgz9tZLCgZwfNZjOys7MRFxeHGTNmoF27dprnuFwuVFRU4ODBg1i3bh0vdBssJCYmory8HH/5y194XKnL5UJZWVkLjywwaNDyWa9evXS3u8+WG0NFRQX27NmDRYsWoV27drjxxhsxduxYFBcXY+TIkX55jUBALyVfD/Uym8FgwOLFixESEqIpJKgOPCQajiAICAkJwfjx47F//36Plh9q0aoY8ZqaGt59mggemsNGBjoNXV4eO3YsgNr4IgBISUlBQUEBnE4n3n//fTz88MN47bXXeHiAu+foiiuuQGVlpZ/fRevEYDDg4MGDWLFiRUsPJWBpkCi68sorsW7dOnTp0sVj3+DBg/0yoC+//BKdO3fGFVdcwbddd911+OKLL/xy/kDClwwm9yJmNpsNO3bs4AW8SkpK6vSAEPUjCAKys7ORlJSEY8eOobi4GFdddRWSk5M9Zl/KLHby5Mm8cGZERAS5rIOE5rCRgYy3ukN6x1ksFhQWFvIig9HR0Zg6dSr3fgC1cX6vvvpqna/Z1gVRVFQUvv/+e1xxxRW6gdXA70tqRP3UK4p++eUXXHLJJQCAO+64A1999ZXuDT9w4EC/DOiXX35B+/btNdvat28f1O0p6ppZqTPP1I1fASAuLk63izTRcGw2G1auXMm9bmfOnMGxY8dQVFQEu92OoqIiHgxvs9kuqvs0EZg0t40MZCwWC/e02u12WCwWj3tD3ZJI8frY7Xbk5+drBJGCEkSsBGMzxoKqcv/nn38OWZbxzTff6Nr6Hj164JlnniEb5CP1iqKbb74ZgwcPxr///W9kZGR4Pe6f//ynXwbUoUMHDwF04cIFXHrppT6f47PPPvPLWEpLS/1ynsZw4MABpKSk8JnV8uXLERkZyfebTCYsXboUpaWlGDBgAEwmE0pLS1FcXKwxKgpXX301Kisrg8poNBZBEFBaWorDhw9rtr/66qtcdNrtdhQXF8NkMiEsLIx3nzYYDAgLC2sV3yWgdXyn2xr+tpH+sl9A6/u8v//+e4/H7mMsLi7WtCRSaoN16NDB63kVMRBsdk0QBJ7w4Y2KigpMnjwZRqNR89vR2r4bvtLk466v+mPPnj1ZTEwMe+KJJ9ivv/6q2ffFF1+wlJSUhpeUrIPjx4+z8PBwZrPZ+LaCggL26KOP1vvctljRev78+UySJAbUVqZWKrjqVXJVjzk3N1e3kmlubi7r0aNHi1ddDaQ/pYK1wWDQbB8zZgwzGo28krj682iNFa6ponXT4C8bGQwVrXfv3s1CQkKYIAgsJCRE14YpFfolSWImk4klJyez3NxclpyczG0h/dX+GY1GZjQafbJh8+fP11znQKTVVLTOysrC4sWLkZiYiLy8PHTs2BFArQdnx44dvpzCZ6699lpER0cjKysLzz77LL7++mtYLBZMmjTJr68TKLgXZgwNDfWpfpFepoHT6URqaqqmLxpRP+y3WZjT6cSgv8N5+wAAIABJREFUQYPQrl07hIWF4Y033uD9hbKzszWfg7eO3vXVnCICk+a0kYGM2WzmMY/eAq3dk0wA76EASj9HFmQeIoV77rkHW7Zs4cuGSjC/+/UwGo3Ug9FHfErJv/zyy1FUVARBEDB27FhUVVU16aAWL16M8+fPY8iQIXjyySfx0EMP4aGHHmrS12ytKAZi3rx52L59O2w2m0fKt3tKvtVqRWFhoce5XC4XCaJG0qVLF7z77rvo2LEjb63CGIPNZvP6HErTb/s0t40MZMxmM2bPnl3nxEB9jLqNkRpBEHgNo2BEKTVgt9vBGIMoirjzzjs9YteuvfZaimtsAD5nn1166aUoKChASkoKHn/8cbz88stNNqirrroqKAs2esM9C82b58hgMOCmm27SBDMS/iUsLAwAePyW+nPIzMysNxi+vppTRODSnDYymFCXuVALo2AVQ0CtIHzkkUfwxhtv8OtgNBqRkZEBi8WiyUC7++67SRA1gAal5Ldv3x4rV67E5MmTMWbMGKSlpTXVuAgvKCnfa9aswejRozWeI6A2uyM/Pz+oDUZTIQiCJq1VqZ8SHR2NtLQ0r8tjvtacIgIfspH+QZ1xC9Tea5WVlfjxxx+xa9cuMMaCPqv2v//9L7f7giAgMTGR25aCggJNvzjCd+oVRS+99BK6devGH5tMJixbtgwzZszQdAUnmger1cp/gD/88ENkZ2drUvIB0BJZI/DWWRqonZlmZWVh+PDhmDx5Mq+fAsBjecxd+PhSc4oITMhG+hf3grSMMTidTi6CRFFE9+7dceLEiXrP1b59e1y4cKGJR9z8sN/aniilCNR1iJQlR5qEXRz1iqIRI0Z4PslgQHZ2Nv7+979j8+bNTTIwQh/3+BSbzca9EGFhYbjppptQWFioW8+DqJtevXph8ODB+OCDD3Do0CHdY9atW4cNGzbwGVpNTQ0qKyt9Wh5zn/2S0WobkI30L2obp5dqL8uyT4IIABwOR1MMsVUgSRKWLFkCm81GdsSPNGj5TI0oipg3bx7mzZvnz/EQ9aAXn6J4IZRaRTt27MCsWbOwb98+nD9/3uu5Lr300jr3BxunTp1CXl5evW55RRAphIWF1bs8pp79SpLE64tQNlrbhWxkw7BarSguLkZkZKTG+92YBJG26jUXRRFLlixBUlISn2wBtV4iynZtHBctioiWwVt8imJQ7HY7ysvLNYF2vXv3xsiRI/HCCy9oZlzBLIg6deqEs2fParbVdT06d+6Ms2fP8iwP9lvVXGXNXhGmSiaguzjyNvulitcE8fukoaamBiEhIcjOzuZlRaxWKz799NMWHmHzYjQaIQgCHA6Hx3K+KIr429/+BpvNhpkzZ2LRokVwuVwICQnhvw31LecT3iFRFIC4x6eoDUpBQQEiIiI0x19zzTWIj4/HCy+80NxDbbW4CyJ3JEmCwWCA0+mEJEm4cOECF0KyLMNgMGD8+PFcEAF193VSe/jcPUWUjUa0ZXxpAKtOu7fb7SgrK0NRUVHQtinq06cPrr32Wqxfv95jX/fu3ZGTk+Nxbaqrq2GxWJCQkMCz9URRpKbUDYREURvA3aCogz6B2oaBJSUllJHmI6IoYtmyZYiIiEBJSQlOnjyJvLw8fv0UYdS9e3eNka+rr5NeQTqKKSLaOr4u5bj3cNy3b1+DBVFdSRKBxqeffuq1kXeHDh106zYxxpCfn4+EhARkZ2cjNTWVmlJfBD4VbyRaN4pBkSQJJpMJ6enpSE9PhyjWfrwLFy7EO++8A6PR2MIjDQySkpKQlJTEC8gpMy/leoqiyJvvNgR1QTpfCtgRRKDja+FSZdIQHx8Pxhj27t0LWZZ5hWb1v0rRQncuv/zyJnkPLYE3cWc0GjFt2jSNPVLjcDhgsVhgs9n45I0KxjYM8hS1ARSDUlxcjDFjxsBsNmPWrFl8JiHLMnbu3AlJkiCKYlC6oxvCuXPnPFz+ipcnNDQUNpuNF3FUxw8lJCSgsLCQz4qpPggR7DSkcKnZbEZYWBhPvxdFEUOHDkVUVBT279+PqKgonD17Fnl5ebrP/+mnn5roXbQ88fHxiImJ4bZG8WKHhoYiPz8fe/bs0RxPBWMvHhJFAYjeGr3ZbIbJZMKAAQMwc+ZMTaC1gnvWFKHPq6++yvuaqV3+aq/OnDlzMGHCBH6MEj9UX18ngggm9BJD6ooxuuyyyzSTuerqavz73/+Gw+HA+++/j5EjR3qd1LWVpTN3BEFAenq6196KADSiKDo6mgrGNgISRQGG+xq9OksjJiYGpaWlFFDdSBSXMwDd7A2r1Yrnn3+ei0x1/BAVaSQILQ1pjvzzzz9rvNnqyZ0sy1i/fn2bih3yBaVorLswUrDZbPyaiaLI+zCSLbo4SBQFGOo1+pqaGqSkpHADoriVg8lgNDWSJHm4nktKSsjrRhAXQX3p4gMGDEBISAiqq6t17ZiSARpsrFu3Dps3b9ZNr4+NjUVISAgtlfkJEkUBRmhoKK+To3SJVvD2Q200GnUruwbbjOtiGDdunK4RUl9T6i9EEL5RX6xLZGQktm/fDovFgsLCQjidTgiCAFmWgz4WUvFIA9rMVVoq8y8kigIIq9WKqVOnwuFwQBRFXHPNNaioqOD7BUGAKIoacRQTE+M1tZMEUd0YDAZdsWM2m5Gbm8vX8aOjo2GxWHiNEDJKBOGJEkuUnZ1dZ2sK5Yc+OjqaN74GgNTU1DZbodpX8vLyuFh0j3cku+MfSBQFEBaLhfc0k2VZI4gkScLMmTPRo0cPpKamQpZlhISEoH///igtLW2pIbc4HTp0QHV1tccsUxRFREZGeghGQRB42u+SJUs8YomU2VhkZCQSExN5wUYlBqmwsBA7duwgA0UQKnytV3TgwAFs3boVoaGhmDJlChwOB29uunPnTlgsFqxcuTJovUZKRfz6quH7UjCT0IdEUYCjpK1mZGTw7DMlXVNxTRcUFPgcA9OrVy+vzVADkV9//RWiKKJXr144fPiwxjtWXl6uOVaSJCxbtowHrqsrg7tXq16xYgUGDBiAkpISzdIkldUnCE98aT1htVqRkpLCvUHuiQzLly/npS/++te/4uDBg83+PpoD97AG9WOj0QhRFOushk+9zxoHiaIAQFH90dHRMJlMvB+Oslw2evRo3hAW8Mw6GDduHHJzc+tdLhNFEWlpaW3OTS3LMg4fPsyvnWJk1EJREASMGjUKAFBUVAS73Y6ioiLu6t+zZ48mI23Tpk1ITEzkMV7KuSjQkSA88RZLpPZoKBMMPS9QZWUl/7/ZbMaqVaswePBg3VjJQEfpr6hkkynbFO+1etKrJ3ao91njIFHUynFX/Tk5ObDZbKiqquKNAJUy7iaTSfccCQkJPvURYoyhrKwMS5cu1WS1tQUYYzCZTLjyyitxxRVXeDSYZIxh/fr12LRpE5xOJxhjuHDhApKTk/mSmjtWqxVpaWncYI0aNcpr2ixBBDPe6hW5lxdREhh8sT1ms1m3HlsgoYge9ftVe4bU2xljsNls9cYPUeHGxkGiqJXjrvptNhtmz56NzMxMnpGhzAaGDRsGwHM9WW2QqqqqkJWVpXkN5SZkjKGgoADjxo1ribfa5Jw7dw7nzp3D6dOndfczxrggUm9Te+VkWYbRaMQ999zDPxtZliFJEmJiYkgQEYQX3H/M9Wzb8uXL8cILL3gs4e/fvx/3338/wsLC0KlTJ7z44osBXxZDaVeiCB+1HVY/VuNLc1fKRmscJIpaOd5Uf13uaL31ZLVBKikp0VRAvfTSS3H+/Hm+pFRZWekxc1HSYts6kiTpLh2KoshLISizO19mZBTwSBD66N0/Gzdu1I1pPHHiBE6cONH8g2xC3Jfv//CHP2gmbL169cLRo0d5uIQsy5g8eTKA2v6MdUHZaBcPNYRt5Siqf968eZqAOW/b9daT3Rk/frzm8blz5/iPvbcluGAQRADw8MMPe1wDZWlMyfpwOp0oLS31+hkoKAJ17ty5iIuLg9VqhdVqRWZmJqxWa3O+LYJodejdP++//35LD6tFYIzhm2++0Ww7evQocnJycOedd/Lle4fDgdTUVLIfTQh5igIAb6pfb3t93gur1QqbzYb09HQUFxfjzJkzALRZbLNmzdI8J5jqGX3//ffIyclBfn4+2rVrh969e/NaRe+++y6/rgMGDABQ94zMXaBaLBYexE1ZIQTx+/2jTBbCw8Px8ccf8/2BWmDWYDDg0ksvRVVVVZ3HeYsfAmo9STabDRkZGdixYwcPKpdlmYKnmxASRW2MutaT1UtrkiRpjI0oirjuuuuwbt26gA9ebAxRUVFIS0vjwmXBggX8Gqqvq9qb5L5EpjwODQ3VCFQAlBVCEG6o7ZIoihg0aBCqq6vRrVs3vP32215FUceOHXHu3LlmHq1vOJ3OegURoBVF7llnyqTWbDZjyZIlmDRpElwuFwwGAwVPNyEkitog3rwXas+F+6xElmWsXLnSwwApP+ptmTFjxmDfvn0QBAFHjhzxKlzU11Upf6CXQaMWVerqvQA0niIybAShtUsulws7d+706f44f/588wywCXGfmA4dOhSjR4/2qPgdEREBg8HAr5FSY43iFf0PiaIgIjY2FpIk8WwpQRA0Llk92mIdEHdOnTrFgzsPHjwIo9EIwLeaQ+5LZGvWrNHNFlSgrBAi2Kgv2UDxvF64cIFvs9vteOGFF7wWMfz/7d17cFTl+Qfw7zlnL1mQyhjbbtGhgIFRNCQIJk0ZJAoNQnBAaOsoNAoIRPASRCOobVMvoIjADKQYbinROGqJI0RE0Nag1SAXgwHKVUXkEoFY5Zrdze77+4PfOe7ZPXsJ2Wvy/cxk2uyePfvmuPvynPd93ucFknNa/7LLLsOQIUOwZs0a3UozWZZhtVpRUlIS8IZW7YvdbjemTp0KRVHgdrs5FR9hTLRuZ9SEPUVRkJ+fb1h/R5IkpKWlYdSoUYbPJ6IOHToEfV5RFFgsFiiKAqvVqq0gA+A3Xdi3b9+AydO+1A5dURSYTCZ06NABJpNJez/foConJwezZs1iB0btgtFiA/VxdcGBOuXfrVs33Wt9g55kDIJ8nT17Fh07doTJ9NN4hCRJGDJkSND+Jjc3V9dneTweuFwuwwU1XMzROhwpakdqamq0OjzNzc2w2+2GS9BTUlK03Zg3bNigu4NLFOrcu6qpqSngser2HWol2N27d6OysjLg8RMnTgy55FWldugVFRVYuXIlqquroSgKJk2axM1hqd0LtBrWt2wIABw+fDiOLY2dtWvXYvr06Zg/f762R2WgESJVTk4OSktL8cADD8DtdsNsNmslVMIpyULhS7igqL6+HnPnzsXevXths9kwdOhQzJgxAzabLd5NS3q+K9P69u2LESNGoLq6GkIImEwmTJgwQfvHvLa2FkOHDsX+/ftx4MABbUl6vJfnZ2RkYOjQoXjppZfg8Xj8aihJkgSLxYJhw4bBbrfr/p7Dhw/jtdde8zunJEm45ppr8NhjjwUNiNSpALvdrluBVlNTo833A0DXrl3ZGVG7593nmEwmHD58GBUVFYaBUlsYCQrHmTNnsGDBAowYMULXP6kCTTdOnjzZb19L3+O4xUcEiARy+vRpkZWVJcrLy4XL5RLHjx8Xo0ePFk8//XRYr29qahLbtm0TTU1NrW7Ltm3bWn2OWAunzZ9++qmYPXu2KCsrEzabTSiKIiwWiygsLBSffvqp7nmz2SwACADCbDaLwsJCcfPNN2uPxfLnuuuu0/2uKIqQZVmYzWa/NnXr1k37e9S/ubCwUFitViFJkt+5JUkSNptNOz7YtVOvmdVq1R3v/Zz39UwkkfhMR/I7RnqRvraJ0oep3z+LxaJ9P6xWq1AURfveffrpp7r+pj38GPU73v1IOH2S0bUO5/WJ8tloqVj0YQmVU3T06FH0798f9957L0wmE+x2O0aOHImtW7fGu2lthprT0tjYqFvxAQAVFRW45ZZb8Oc//xlTp07VJVmrU2yfffZZi97PbDbDbrdr8+GSJCEzM7NF55BlGd99953uMXUFndou7/n2Q4cOYenSpRg8eLD2v2VlZXA4HH6rPcaOHYvnnnsurGFm77uw5uZm3Ty+Oo02adIkSJKEZcuW6XIoiNqrnJwcdO3aVetr3G43xo8fr8vby8nJ8Ssq29YJIeBwOHT9SDjFd4MJVVCWQov59JnD4cDJkycNn+vevTtKS0u134UQeP/993HDDTfEqnnthvdKNEmSUF5eDqfTqdt3x5v6uxqESJKEkSNH4p133jHcFkPlcrnQ0NAAk8mkbcr6m9/8Bl988UXI4XJ1Ww2TyYT+/ftj48aNfscIIQzrKql7wlVVVfkFQ8BPeUZGU2WBhq99pwKMkqjVvC0OXxP9xHfq3ijfrqCgACtWrLjkFa/qtLnD4YhEk2PC4/Ho9jOLxGau3OKjdWIeFO3YsUOrEOyroqIC2dnZAC7+41tSUoJvv/0W8+bNa9F77Nq1q9XtBH6qRZNMwm3z3r17dZVU3W63LiBSa2J4b1Z4xRVXaMGNLMvo3bs3FEVBVVVVyPdTR6Oam5uxd+9eyLIcckNH9b2FEH61kkJVulX3azt9+nTAHKj6+nqUl5dj+/bt6NevH/r06YP6+nrcf//9cLlcMJvNWLJkCfr06QPg4hL90tJS7XiLxeJ3ve12u7ayRB3tTKTPUSK1hYxFqv8C4vffu76+Xve9Cue7s3fv3oDfVfUGyeg7772RajIERN6LRGRZ1q4VEF4fEynJ2hdEvd2tnqCLglOnTomCggIxYsQIceTIkbBfx5yi8Ns8e/ZsoSiKAKDl5iiKIkwmk8jKyhJlZWWisLBQy8FRFEXMnj1bFBcXC1mWtTyj4uJi7TzBfrxzeSRJMsztCfQjy7L2nqGOy8rKEmPHjhUmk8nwPdS/02az6fKq1Pl37+ui/s0tvdbqeZhTRC3RVnKKwslrMfqOzJ492/A7a7VaRXFxsTCZTH7PZWRkiG7dusU9P6glP5Ik6foh71zOWPUZyfjvmxCx6cMSbvXZvn37MHnyZGRkZODvf/87OnbsGO8mtUm+w7QLFy7E+vXrUV1dja1bt+Lzzz/HI488gpSUFO2Y1NRU/OUvf9HuclwuF+bPn48ZM2bgpZde0kZ+vEdxZFnGXXfdhbfeegtNTU0B7/YCkWUZiqKENaTu8Xiwbds21NXV6Ua+VCaTCYsXL9aqxRrN33P4mqh1Qq2ACrRsPDc3FyaTSfuuq1P0w4YNQ1VVlW7UWu1jvvjii7j8jb7UQorhkGVZ1w8BP5UoMJlMGD9+PMt5xFFCBUWnTp3ChAkTMGLECF0VYIosNWfGdwsKtQYGcHGaa8GCBX5BhNH2IKdPn8aAAQO03B7fYKSurg7du3fHlVdeiU8++STszgOAVvb+wQcfDGu7EY/H49dGdTpw8eLFfjlEvgFQsL3jiCi0UDcWgYKmnJwcDBkyBOvXrwdwsR/p1asXioqK4HA4tD3B1NpqLbm5aom0tDQcPHiwRa8ZPHgwevTogYaGBqxbty7oTdztt9+u64fmzJmjW/RSVlaGVatWMVE6ThIqKHr99ddx6tQpvPHGG3jzzTe1x7t06YJ169bFsWVtR6C7tDlz5vgFK83Nzairq0NBQYG2wanVatVGfCRJgizLWLZsmd9r1dVgHo8H//3vf7XHO3fuHNZGiQB0Rc3S09Mxd+5crTx+S/zud7/TzuO9WWtjY6MuMDTa44yIWibUjUWgoKm2thYbNmzQHfvWW2/B6XRqAVH//v2RkpIS1U2rv/766xa/ZuPGjSguLsaSJUtQW1uLoqIibNmyRXtezSMym80oLi7WvVa9Ht4j6VykEUetnqBLIMwpCt3mQDkzah6Ab+6O2WzW6ouoeThqHaNRo0YFzA0qLi4WeXl5Yc2xd+/eXWRmZuoe6927t+H8ellZWVg5TN4/ZWVl2mvNZrPWZlmWL7mmUFv9fITCnKLoaSs5ReHwzqFR6xhlZWUZ9iNqfTG1L2pJPmIsf9LS0nR/n3e71RxHi8USMMfKu5bTpdQoaolE/mwE0+7qFFHkGO1/o1Z0VhTFb28u9e7u2Wef1e151tzcrNtjRx1Vqaurw9q1a/1GbSRJQnFxMUaNGoUePXpAURS/tkmShM6dO6NTp06QJAmHDh3Cnj17tGPNZjOWL19ueJeUnp6ubWYbDkmS0NjYiNraWkybNg0ul0u36s7pdKKsrIw1hYhiSK2XBlwcKXn55Ze1kRV1BLq4uBjXXHONNlXmcrl0ZUN8/eIXv4hZ+4106dJF14cMGzYMN910E/Lz87URILfb7Vd7SB29VkfkWWMovhJq+owiw2iKDNAn8xntzaVOG6nD2Oqx6pdZTbYePHiwNtTrS5ZlnD59WvdeAwYMwKlTp7StQhRFwblz53Tz7i6XS5dAOXfuXMMS+Op2GuL/ywIIr8RtRVFw++2362onqYGfUT6USvgMV4fa1ZuIIsN793eVut1Oeno6Bg0aFHYO4okTJ6LRxLB9/PHHyM3NxaJFi/DQQw9p5QFMJpNWpsM3x8qor2Y+bXwxKGqDjBIZDx8+rAtkfPfmWrp0KaqqqjBmzBhMnjxZlxOgnjM1NRUrVqwIuEGsmgT5+eefa4mRAHDbbbdh1qxZ2nt06NABa9as8Xu9GuA0Nzfj7bffBgCUl5fjww8/1NrqnY+gKAqGDx8OAH57nKkb2noHVVarFQ6HA7Is484778SBAwdQV1cHj8ejdVbcUJEodnxXnAHAV199haKiIgwdOjRoYdhEo95cPfPMM7pFIW63G/fddx+6du3qd6PFvcoSD4OiNsg3kTE1NRV/+9vfdCMqqampmDNnDnJzc7Fz505MmTIFALSq0erqCDUwSk1NxdSpU3V3bWazGfn5+bDb7fjZz36Gmpoa1NXVYdu2bVpipHewUVRUBKfTCVmW/Yov+i7bVzkcDhQVFWHhwoXaSJa6K315eTmqq6u1Zawqo0Rp7+TP1NRUrS2+O9p7rwRhJ0UUfb5T4R6PBw6HQ9uouiXMZjNycnLwySefQAgRdjmPSDpy5Ijud0VRAi6xj0QJEIosBkVtkO/qD3XrCeBiBzR8+HAtKLBYLEhPT9e9vqqqCunp6dqIiboE1nf6aeLEidpqC98pNVmWteX06kiV715rKlmWcf3112PatGmYO3euNkqk2rJlC2655RZtxMh3O41wl7Gqr/UOfAD9qBk7KaLYqaioMKxWD0DX34SqYK/Kz8/Hhg0btNWx+fn5OH/+vOEWQcF07NgR586da9FrfKltNsqr9J6iZwmQxMKgqI3yHS3x/ofebrfrRkO6dOmie60ayKjHeDwevw5J3b8I+GkIWD1GkiRYrVaMGTNGNyJjMpn8ziVJEsxms9YhZGVlYe3atX4BmO+ojdEy1qamprBGdoIFPqxTRBQbtbW1KC8v93vce9sOlclkQn5+ftAaQN59m9p/+N5ghatnz57YvXu37r0URYHdbsfRo0cNX9O7d2/s379fuwFV2+90OlFRUaH1JcwjSmxcfdYO+O6cXFBQAIvFoq1AKy4uRllZGfLy8lBWVobJkydrgYOiKLod6AHg17/+NSZMmKD97n2sxWLBlClT8K9//QuNjY260aHx48djypQp2rFWqxVTpkzBkiVLdMGO1Wr1e0/feiZq8cmBAwdqxwghwqqBFGonaXVlDAMioujxHsEGLgZDiqIYVr33eDyw2+3Iz8/H1VdfbXg+SZLQt29fWCyWsFenBrJz504sXrwYhYWFuPnmm3HdddcBAI4dO+bXXkVRYLPZsHz5ctx3331+7y2EwMqVK7WVaUZ5RJRAWr3oP4GwTlH4bQ5nrx31mLKyMq0+iMlkElar1a+WhtH5Au2BVFZWJvLy8rT6Qb7tVmt2qDWFFEXRjvU9p29tk7y8vLCvQWu09c9HIKxTFD3tqU6REPrvslovrKCgwK9OmqIowmq1Gu595v0jy7LWBwWroRbOj/e5bDab37kkSdLqnI0ZM0bXD6p/k6IofntH+h4T7XpEgST6ZyOQdrn3GcVGOFWbvY9JT0/XcoPUCta+JfqDJTerU1HeCdcff/wx0tPTYbFYtNeoo0AAdFNtjY2NAPzvsoym/ogo8fhWk1fzabxXilZWVmrVn2VZRn5+Po4dO4YzZ85gz549Qc8vy7I2mqzmFV0qj8eD1NRUv9QAlfCqO2S323XV8I0WdBjVhOMUfWJiUERh8a5htGrVqrATkX2DJaOh47y8PAAXO83c3Fy4XC4tBwnQT5355gMVFxdrG0aq5QSIKP68k4mBi3XSvPcwM5vNGDZsGNavX4/m5masWrUK06dP10pnSJKkqznmKzMzEzt27NB+f+SRR3QLKVqrsbERqampfo+bzWZIkqTVbuvXr5/ueaObSfUaqCt+uZVQ4mJQRC0Szl1OsOKHwZKcvVeiNDc3Y8SIEcjKyvLbl8y3hlJjYyNKSkoA6DsdIooP32Tie+65R5cArS67906Edjqd+PHHH7Fw4UI88MADITd97datG+rr67Ugq3PnzgCA1NRUrSq2yWTC8OHDYbfb0dDQEHbitaIoePvtt3X7lymKgn79+mHixIkAoN2I9enTJ+B5vG8mWf8sOTAoaqOiVZU51HlDffmNgqrt27cbvpfdbjdclWHU0ahbfzQ3N7PTIYoz3xFh4OKIr3dRV2+SJGmjLg0NDYYrXhVF0R6XZRl2ux1Wq1VXj+3+++9HeXm5Vjl/0aJF2uhxbW1twNVrsizj0Ucfxf79+1FdXQ23260LiICLRRi3bt2KL774Qps6+/jjj1FaWuo3WhTqerD+WeJiUNQGReuuJJzzhvPlDzR0XFBQgPLycu386pL/QHzLBgD+W3YQUez5jggXFBRoe3ulpqairq4OK1euhNvthqIomDBhgrYq9tprr9VeK0kSevXqhV69emHYsGG6Eh8AsHDhQm2aq6ioSFcrTQhdDiHiAAASSklEQVSByspKrFixAl26dEGvXr0Cbhmi7seYlZWF6urqgH+XEEKb2hP/XwZk3bp1uuKx4VwP1j9LXAyK2qBo3ZWEc97WfPlzcnLw4Ycfhj3C5bvlh/dIETsdovgJNM3u/Z1WgyTv57dv3669du7cuaiursaePXtw4MABDBs2TFfNftmyZdrNWUVFhV81fCEEPvroo5BtVUep1D7DYrEE3MoIuDiqJMuyNr23du1a1NbWBu2vmFydPBgUtUHRuisJ57yt/fK3JAHRKL+InQ5RYgj1XQ71/DvvvKON7LhcLkybNg0fffQRunbtqlWzv3DhAmbOnIna2tpLXm0myzKGDh2KnTt3orGxURt92r17Nz777DNkZ2dj9erV2gavsixjxIgRWLNmDYQQ8Hg8Yd14Mrk6OTAoaoOidVcS7nlj+eX3fS92OkSJy2hZvtF3tqamxi/3SA0+cnNzIcuyFjAZjQapy/rD4Xa7sWbNGrz99tuQZRlWq9WvynSnTp1QVlamLcW32+1ISUmB0+mEyWTiyHQbwqCojYpWYMK7HSK6FGpOoveyfDUAMZqGt1qt2pSYJEkwmUw4fPgwACAtLS1g3SKLxYKioiK8+OKLfqNH6tSXx+OBx+Px21LE4/EYpgYUFBToSpF450h51ymi5MegiAKK1go2Imp/1JxE72X5wRZjeBdBrKur0/KIysvL8atf/crwPbxXnG3evFk3itStWze89tpr2LlzJ6qqqpCZmamtNlOnwWRZNkwNCJYjFWj1LCUnBkVkiHU1iCiS1JxE75EiWZYNCyQC+lHpO+64Q8vpcbvd+OabbwxfI4TQqt8///zzGDRoEFwuF8xmM1577TUA0Fawbdq0SRslUpfkd+7cOeBNoO8ouXrTaLfbQy7Jp+TBoIgMsa4GEUWS92jLDz/8gPnz58PtdqOoqAjp6ekB+5elS5f6FV1Up9SuuuoqHDlyRHtcURTddhqbNm3Sje6o1a59y3ioS/LD3a3e+6bRZDLh2muvZf/YRsihD6H2SL2rUxSFS9yJKCJycnIwa9YsdO7cWZuyCrVTfFVVleHjQgjcfffdsNlsWvXqxYsX+y28mDVrlvaYb79mNpsvqY/zvmlsbm4O2n5KLhwpIkOt3c6DiCiQlpQNGTNmDDZu3Oj3uLq1h3fuUWNjY9CaQWq/pm5C27dv36Cr4ADjfs67/Vx91rYwKKKAgq00Y84REV2qlpQNUbfpWLFiBVJSUvDZZ5/pirR6Bypq/lCo6X7vlWTB+q5A/Zx3+7n6rG1hUESXhDlHRNQaLSnvoQZGVVVVePjhh/0Sor03k3Y6naioqPA7tzris2XLFm2pf6i+K1g/p7afq8/aFgZFdEm4lw8RxUJtbS3mzp2rJVtv3LgRZWVlQQOqhoYGv3N410hSeSdmG/GdJjt8+HDILT0ouTHRmi6JOnz8zDPPcOqMiKJCDWZ8V5/5Jl8XFBTAbDZrv69fvx61tbXa7zU1NX4BkSRJmDBhQlh7lk2aNAlCCCxbtgyDBw/WnZvaloQNitxuN/70pz9h5syZ8W4KBeC7soOIKFJqa2tRUlKi1SfyNmbMGN3vOTk5mDhxIiRJAgC/FWGpqam6gEiWZaSkpKCgoCBkO3JyctC1a1e43W7dNBq1TQkbFC1evBjbtm2LdzOIiCjG1BGiDz74QCv0aDKZkJWVhbKyMi3HyFtBQQFSUlIMl9g3NjZCli/+cydJEoYMGdKiEW6WKGk/EjKnqLa2Fu+99x5uu+22eDeFiIhizHu6S5ZlDBkyBCUlJWFNdRmtaFP3UnM6nVAUBT169GhRe6K1yTYlnpgHRQ6HAydPnjR87uc//znOnj2LJ598EosXL0ZlZaW2EzIREbUP3tNdHo8HY8aMQU5OTsjaaIFWtHnXJ1L3UFu1alWLRou4GXb7EPOgaMeOHQHncf/xj39g2bJluPfee9G7d+9Lfo9du3Zd8mu9JeNSy2RsM5Cc7WabKRoi1X8Byfnfe/v27aivr9d2sJdlGfX19SgvL8f999+v1SJasmQJ+vTpE/Z51ekvl8ulVdKurKyExWKJSJuTTTK2GYh+u2MeFGVnZ2Pfvn2Gzy1ZsgQWiyWs5LdgbrjhBlit1ladY/v27Um3yV8ythlIzna31zY7HI6I/qNN/iLRfwGJ+xkNNtqjttnpdGLlypVayY+xY8eipqYGzc3N8Hg8aG5uRkNDA8aPH9+i9zY6b2uvUaJe52CSsc1AbPqwhMopWrNmDU6cOIH+/fsDAJqamgAAu3fvRnV1dTybRkRErRRuJfxAOTytrRnE3CAKJaGCovfee0/3+5NPPgm3243nn38+Ti0iIqKWCjQa1JJK+L45PN55QStXrrykvCCj8xJ5S6igiIiIkluw0aBwK+EHCqpycnJQU1PjVzMokkEON7pu3xI6KHruuefi3QQiImqBUPuFhZq+CjXFlpubC0VR4PF4Qm7T0VLc6JoSOigiIqLkEmo0KNj0VX19PV5//XWtRlGgkSC1crX6v5HCja6JQREREUXMpSYz19bWakvu1aKNRkGVugpNCKFt59GawMV7uowbXRODIiIiiqhLSWauqanRBUSBqli3NHAJliNkNF3G1WntG4MiIiKKu9zcXJjNZjQ3N8NisQTc1qMlI1GhcoSMpsu4yXX7xqCIiIjiLicnB0uWLEFDQ0PIYCfckahQOUKcLiNfDIqIiCgqWrq8vU+fPi2uUh1MOEnfnC4jbwyKiIgo4hJheXs4QQ+LOZI3BkVERBRxibK8nUEPtYQc7wYQEVHbo05dKYrCfB1KGhwpIiKiiGO+DiUjBkVERBQVnLqiZMPpMyIiIiIwKCIiIiICwKCIiIiICACDIiIiIiIADIqIiIiIADAoIiIiIgLAoIiIiIgIAIMiIiIiIgAMioiIiIgAMCgiIiIiAsCgiIiIiAgAgyIiIiIiAAyKiIiIiAAwKCIiIiICwKCIiIiICEACBkU//PADZs2ahZtvvhnZ2dkoLi7GuXPn4t0sIiIiauMSLiiaMmUKmpqa8O677+K9997D0aNH8eKLL8a7WURERNTGmeLdAG/19fXYu3cvysvL0aFDBwDAwoULOVJEREREURfzoMjhcODkyZOGz+3cuRPXXHMNXn31VbzxxhtwOp0YOnQoZsyYEeNWEhERUXsT86Box44dKCgoMHzu4Ycfxt69e3HjjTeiuroaP/zwAx588EG88MILKCkpCXluIQQAwOl0RqStDocjIueJpWRsM5Cc7W6PbVa/W+p3jSIn0v0X0D4/o/HANsdOtPswSSRQ77ZixQosWLAAdXV1MJvNAIAPPvgATz31FDZv3hzy9WfOnMH+/fuj3Uyidq9Xr17o1KlTvJvRprD/IoqdQH1YQuUUpaWlQQgBp9OpBUXNzc1hv75jx47o1asXzGYzJEmKVjOJ2i0hBFwuFzp27BjvprQ57L+Ioi9UH5ZQI0VOpxPDhw9HRkYGnn32WZw+fRqFhYW46aab8MQTT8S7eURERNSGJVRQBAANDQ2YPXs2tm3bBrfbjfz8fDz++OOwWq3xbhoRERG1YQkXFBERERHFQ8IVbyQiIiKKBwZFRERERGBQRERERASAQRERERERAAZFAC4WiLzuuuvQt29f7eexxx4LePzq1asxePBgZGZm4u6778bBgwdj2Nqf1NfXY9y4cejfvz8GDhyIZ599FhcuXAh4/KhRo9CnTx/d3/nll19GvZ2nTp1CYWEh+vXrhwEDBmDx4sUBj03Gaxuv6+qrJZ/jRLnOFBnJ2Iex/4qeZOy/gATpwwSJBQsWiOnTp4d17ObNm0W/fv3Ejh07hMPhEPPnzxe33Xab8Hg8UW6l3unTp0VWVpYoLy8XLpdLHD9+XIwePVo8/fTThsc7HA5x/fXXiy+//DKm7RRCiHHjxomZM2eKCxcuiP3794vc3FxRXV3td1wyXtt4Xldf4X6OE+U6U+QkWx/G/it6krX/EiIx+jCOFAHYtWsXbrjhhrCO/ec//4n8/HxkZGTAYrHgoYceQmNjI2pra6PcSr2jR4+if//+uPfee2EymWC32zFy5Ehs3brV8Pg9e/bAarWie/fuMW3nN998gy1btmDGjBlISUlBz549cc8996CystLv2GS8tvG6rkbC/RwnynWmyEm2Poz9V/Qka/8FJEYf1i6CIofDgSNHjhj+OBwO7N69G5988gluvfVWDBw4EE899RR+/PFHw3MdPHgQPXv21H5XFAXdunXDvn37Ytru7t27o7S0VDtWCIH3338/4Adq165d6NChAwoKCpCdnY3Ro0fj3//+d8Tb7OvAgQPo3LkzrrzySu2xHj16GF6vWF7bYK699tqwr228rquRcD/HiXKdKXzJ2Iex/2L/1VKJ0Icl1N5n0bJjxw4UFBQYPldaWoq0tDQMGTIEI0eOxOnTp/H444/jsccew9KlS/2OP3fuHFJSUnSPpaSkBJ0Lj0a7KyoqkJ2dDeDi/nAlJSX49ttvMW/ePMPjhRBIT0/Ho48+iquvvhrvv/8+Hn74YbzyyivIzMyMeNtV586dg81m0z1ms9nQ1NRkeGysrm24Ql3beF1XX2fOnAn7c5yI15mCS8Y+jP1X/L9XydJ/AYnTh7WLoCg7OztoBDlkyBDt/3fo0AGPPvoo/vCHP+Ds2bO47LLLdMfabDY4HA7dY01NTX7HRUKodgNAY2MjHnnkEXz//feorKzEL3/5S8Pjxo0bh3Hjxmm/5+fnY+3atdiwYUNUP/wdOnTw60AuXLhgeL1ieW3DEc61jdd19dWpUye88sor2u/BPseJdp0ptGTsw9h/sf9qiUTpw9rF9FkwX331FV544QW43W7tMYfDAVmWYbFY/I7v2bMnvvrqK+13t9uNQ4cO6YbyYmXfvn0YPXo0Lr/8crz++uu46qqrAh5bWVmJTZs26R5zOBxR31OuZ8+e+N///ofvv/9ee+zLL79EWlqa4bHJdm3jdV19teRznEjXmVovWfsw9l/Rk2z9F5BAfVirU7WTXGNjo+jXr59YtGiRcDqd4tixY+KPf/yjeOKJJwyP/89//iP69+8vtm3bpmW9Dx48WDidzpi2++TJk+K3v/2tmD17dljHz58/X9x6663i66+/Fi6XS6xevVr06dNHHDp0KMotFeLOO+8Ujz/+uDh37pw4cOCAyM3NFW+++abfccl4beN5Xb215HOcKNeZIiMZ+zD2X9GTjP2XEInTh7X7oEgIIerq6sRdd90lbrzxRpGdnS2efvpp0dTUJIQQ4ujRoyIzM1Ns3bpVO3716tUiLy9PZGZmirFjx4qDBw/GvM2LFi0SvXr1EhkZGSIzM1P7GT58uGG7nU6nmDt3rhg4cKDIyMgQo0ePFps3b45JW7/77jsxbdo0kZWVJQYMGCBKS0sN2yhE4l/bRLquvgJ9jhP1OlPkJFsfxv4repK1/xIiMfowSQghWj/eRERERJTc2n1OERERERHAoIiIiIgIAIMiIiIiIgAMioiIiIgAMCgiIiIiAsCgiIiIiAgAgyIiIiIiAAyKiIiIiAC0kw1hKTG8++67mD59uuFz27dv54akRJSw2H+1D6xoTTFz4sQJHDt2TPv9xx9/RHFxMTIyMrB06dI4toyIKDj2X+0DgyKKi7Nnz+Kee+6BzWbDsmXLYLPZ4t0kIqKwsP9qu5hTRDF3/vx5TJo0CbIs4+WXX4bNZsP+/ftxxx13IC8vD4WFhTh79my8m0lE5If9V9vGoIhiyuFwYOrUqTh//jyWL1+uzcP/9a9/RVFRETZu3IgePXpg+fLlcW4pEZEe+6+2j0ERxYzL5cKDDz6IEydOoLy8HJdffjkA4NSpUzhy5AgGDRoEAPj973+PjRs3xrOpREQ67L/aB64+o5hwu92YMWMGvv76a7z66qu44oortOcaGhpgt9u137t06YLjx4/Ho5lERH7Yf7UfDIooJkpKSrBp0ybMmTMHx48f1zoNm80GIQQkSYpzC4mIjLH/aj8YFFHUCSGwbt06NDU1+dX5uOWWW/DMM8/o7qyOHTumu/MiIooX9l/tC4MiijpJkvD5558HPeaqq67Cpk2bMGjQIKxevRp5eXkxah0RUWDsv9oX1imihLB3717MnDkT58+fR/fu3TFv3jx06tQp3s0iIgqJ/VfbwaCIiIiICFyST0RERASAQRERERERAAZFRERERAAYFBEREREBYFBEREREBIBBEREREREABkVEREREABgUEREREQFgUEREREQEAPg/anwHyLfQPDYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax_grid = plt.subplots(nrows=1, ncols=2, sharex=True, sharey=True, figsize=(10,4))\n",
    "\n",
    "ax_grid[0].plot(z_samples_SD[:,0], z_samples_SD[:,1], 'k.')\n",
    "ax_grid[0].set_title('Gibbs sampler')\n",
    "ax_grid[0].set_aspect('equal', 'box');\n",
    "ax_grid[0].set_xlabel('$z_0$');\n",
    "ax_grid[0].set_ylabel('$z_1$');\n",
    "\n",
    "ax_grid[1].plot(true_samples_SD[:,0], true_samples_SD[:,1], 'k.')\n",
    "ax_grid[1].set_title('np.random.multivariate_normal')\n",
    "ax_grid[1].set_aspect('equal', 'box');\n",
    "ax_grid[1].set_xlim([-6, 6]);\n",
    "ax_grid[1].set_ylim([-6, 6]);\n",
    "ax_grid[1].set_xlabel('$z_0$');\n",
    "ax_grid[1].set_ylabel('$z_1$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}