Created
October 8, 2019 12:02
-
-
Save steven-tey/9829171b9b290d3de090558cf6487ea9 to your computer and use it in GitHub Desktop.
CS146 Session 5.1 Pre-Class Work
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Pre-class work\n", | |
| "\n", | |
| "We consider the eczema medical trial data set again. This time we will compare which of 2 models explain the observed data best.\n", | |
| "\n", | |
| "* Model 1: All studies have the same probability of success.\n", | |
| "* Model 2: A hierarchical model where the probability of success in each study is drawn from a beta prior distribution with unknown $\\alpha$ and $\\beta$ parameters.\n", | |
| "\n", | |
| "\n", | |
| "|Study | Treatment group | Control group |\n", | |
| "|---------------|-----------------|------------------|\n", | |
| "|Di Rienzo 2014 | 20 / 23 | 9 / 15 |\n", | |
| "|Galli 1994 | 10 / 16 | 11 / 18 |\n", | |
| "|Kaufman 1974 | 13 / 16 | 4 / 10 |\n", | |
| "|Qin 2014 | 35 / 45 | 21 / 39 |\n", | |
| "|Sanchez 2012 | 22 / 31 | 12 / 29 |\n", | |
| "|Silny 2006 | 7 / 10 | 0 / 10 |\n", | |
| "|**Totals** | 107 / 141 | 57 / 121 |\n", | |
| "\n", | |
| "\n", | |
| "**Model 1:**\n", | |
| "\n", | |
| "* For each group (treatment and control), all 6 studies have the same fixed, but unknown, probability of success, $\\theta_t,\\theta_c\\in[0,1]$.\n", | |
| "* The data follow a binomial distribution in each study, conditioned on the probability of success — $\\theta_t$ for treatment or $\\theta_c$ for control.\n", | |
| "* The priors over $\\theta_t$ and $\\theta_c$ are uniform.\n", | |
| "\n", | |
| "These assumptions lead to the following model.\n", | |
| "\n", | |
| "* Likelihood: $\\prod_{i=1}^6 \\text{Binomial}(s_i\\,|\\,\\theta, n_i)$, where $s_i$ is the number of successful recoveries, $f_i$ is the number of failures (did not recover), and $n_i=s_i+f_i$ the number of patients.\n", | |
| "\n", | |
| "* Prior: $\\text{Beta}(\\theta\\,|\\,1,1)$ for both $\\theta_t$ and $\\theta_c$.\n", | |
| "\n", | |
| "* Posterior for treatment group: $\\text{Beta}(\\theta_t\\,|\\,108, 35)$.\n", | |
| "\n", | |
| "* Posterior for control group: $\\text{Beta}(\\theta_c\\,|\\,58, 65)$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Since we have closed-form solutions for the posteriors, we can calculate the marginal likelihood by rearranging Bayes' equation: (marginal likelihood) = (likelihood) x (prior) / (posterior).\n", | |
| "\n", | |
| "$$ P(\\text{data}) = \\left[\\prod_{i=1}^6 \\text{Binomial}(s_i\\,|\\,\\theta, n_i) \\right] \\text{Beta}(\\theta\\,|\\,\\alpha_0,\\beta_0) \\,/\\, \\text{Beta}(\\theta\\,|\\,\\alpha_1,\\beta_1)$$\n", | |
| "where $\\alpha_0=1$ and $\\beta_0=1$ are the parameters of the prior, and $\\alpha_1$ and $\\beta_1$ are the parameters of the posterior beta distribution.\n", | |
| "\n", | |
| "Since all factors involving $\\theta$ cancel out, we are just left with the normalization constants of the likelihood, the prior and the posterior:\n", | |
| "\n", | |
| "$$\\begin{aligned}\n", | |
| "P(\\text{data})\n", | |
| "&= \\left[ \\prod_{i=1}^6 \\left(\\begin{array}{c}s_i+f_i \\\\ s_i\\end{array}\\right) \\right] \\frac{\\text{B}(\\alpha_1,\\beta_1)}{\\text{B}(\\alpha_0,\\beta_0)} \\\\\n", | |
| "&= \\left[\\prod_{i=1}^6 \\frac{1}{(s_i+f_i+1)\\text{B}(s_i+1,f_i+1)}\\right] \\frac{\\text{B}(\\alpha_1,\\beta_1)}{\\text{B}(\\alpha_0,\\beta_0)}\n", | |
| "\\end{aligned}$$\n", | |
| "\n", | |
| "We usually compute the log of the marginal likelihood since the results can vary over many orders of magnitude.\n", | |
| "\n", | |
| "**A word on notation** in the derivation above:\n", | |
| "\n", | |
| "* The beta _distribution_ is written as $\\text{Beta}(\\theta \\,|\\, \\alpha, \\beta)$.\n", | |
| "* The beta _function_ is written as $B(\\alpha,\\beta)$. $B$ is the Greek letter _capital beta_.\n", | |
| "* The beta function is part of the normalization constant of the beta distribution.\n", | |
| "\n", | |
| "This is similar to the gamma distribution and the gamma function, where\n", | |
| "\n", | |
| "* the distribution is written as $\\text{Gamma}(x \\,|\\, \\alpha, \\beta)$,\n", | |
| "* the function is written as $\\Gamma(\\alpha)$,\n", | |
| "* the gamma function is part of the normalization constant of the gamma distribution.\n", | |
| "\n", | |
| "**A word on simplifying the expression** in the derivation above:\n", | |
| "\n", | |
| "Just as the gamma function is related to factorials, the beta function is related to combinations:\n", | |
| "\n", | |
| "* $n! = \\Gamma(n+1)$ for integer $n$.\n", | |
| "* $\\binom{n}{k}=\\left((n+1)\\cdot B(n-k+1, k+1)\\right)^{-1}$\n", | |
| "\n", | |
| "The beta function can also be written in terms of the gamma function:\n", | |
| "\n", | |
| "* $B(x,y) = \\Gamma(x)\\ \\Gamma(y)\\ /\\ \\Gamma(x+y)$\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Task 1\n", | |
| "\n", | |
| "1. Take the log of the marginal likelihood above.\n", | |
| "2. Complete the Python function below to calculate the log marginal likelihood.\n", | |
| "3. You can use the built-in function `scipy.special.betaln(a,b)` to compute $\\log \\text{B}(a,b)$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 34, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "//anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:26: RuntimeWarning: divide by zero encountered in double_scalars\n" | |
| ] | |
| }, | |
| { | |
| "ename": "ValueError", | |
| "evalue": "math domain error", | |
| "output_type": "error", | |
| "traceback": [ | |
| "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", | |
| "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", | |
| "\u001b[0;32m<ipython-input-34-e1896a6f6f7d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 27\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[0;31m# Test for treatment group\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 29\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlog_beta_binomial_marginal_likelihood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata_treatment\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", | |
| "\u001b[0;32m<ipython-input-34-e1896a6f6f7d>\u001b[0m in \u001b[0;36mlog_beta_binomial_marginal_likelihood\u001b[0;34m(alpha0, beta0, data)\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[0;31m# Add log\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mbetaln\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbetaln\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m108\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m35\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mbetaln\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malpha0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 27\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[0;31m# Test for treatment group\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", | |
| "\u001b[0;31mValueError\u001b[0m: math domain error" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "from scipy.special import betaln\n", | |
| "import numpy as np\n", | |
| "import math\n", | |
| "\n", | |
| "# data for treatment group\n", | |
| "data_treatment = np.array([[20, 3], [10, 6], [13, 3], [35, 10], [22, 9], [7, 3]])\n", | |
| "\n", | |
| "def log_beta_binomial_marginal_likelihood(alpha0, beta0, data):\n", | |
| " '''\n", | |
| " Compute the log marginal likelihood of the beta-binomial model\n", | |
| " for the eczema data set.\n", | |
| "\n", | |
| " Arguments:\n", | |
| "\n", | |
| " alpha0, beta0: prior beta distribution parameters.\n", | |
| "\n", | |
| " data: 2 x n matrix, where n is the number of samples from\n", | |
| " the binomial distribution. This means the first row contains\n", | |
| " the success counts and the second row the failure counts\n", | |
| " of the samples.\n", | |
| " '''\n", | |
| " # This gave me -inf because the values were too small\n", | |
| " # return np.prod([1/((data[i,0]+ data[i,1]+1)*betaln(data[i,0]+1, data[i,1]+1)) for i in range(6)])*(betaln(108, 35)/betaln(alpha0, beta0))\n", | |
| " \n", | |
| " # Add log\n", | |
| " return math.log(np.prod([1/((data[i,0]+ data[i,1]+1)*betaln(data[i,0]+1, data[i,1]+1)) for i in range(6)])*(betaln(108, 35)/betaln(alpha0, beta0)))\n", | |
| " \n", | |
| " # For some reason, this gave me \"math domain error\" :(((\n", | |
| " \n", | |
| "# Test for treatment group\n", | |
| "print(log_beta_binomial_marginal_likelihood(1, 1, data_treatment))\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Model 2:**\n", | |
| "\n", | |
| "* For each group (intervention and control), each of the 6 studies has a different probability of success.\n", | |
| "* Each probability of success is drawn from a beta prior with unknown parameters $\\alpha$ and $\\beta$.\n", | |
| "* Since $\\alpha$ and $\\beta$ are unknown, we put a broad hyperprior on them — we choose the Gamma(2, 0.5) distribution, which is shown below." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAABaMAAALpCAYAAACuZzrhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XeYVOX5//HPvbssLL0LUkWUolgARUVjNza+Ro0xdmwhRqMpphqjye9r8jXNxN6DvUVj16ixg4UmioKIgNLbSm/L7v374zkzu4wzs4WZPbvL+3Vd5zp7+n3OnBn0Ps+5H3N3AQAAAAAAAACQTwVxBwAAAAAAAAAAaPpIRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAgJwzs9fNzM1sdGM5rpmNjrZ9vab7NbO+0Xyvc9AxyXa+jYWZXR2dw9i4Y0HDYmbFZjbbzFaZWft6PG6j/D1oiMzsLDN7K/oMN5vZQjN70swOjTu2uJlZazNbYWbzzKxF3PEAAFAbJKMBAGikzGxsIvGRMqw2sw/M7M9m1jPuOLdXZnZIlCz9VtyxbI/M7EfR9e/bAGI53szuMrPpZvaVmZVFiaQJZnaTmR1hZoVxx9nEfF/STpJucPeVcQfT2JhZkZkdY2Y3mNnEKgnhRWb2dL5/18zsb5LulXSgpLaSNkjqLukESa+a2Zh8Hj9NPGZm3zOzd8xspZmtMbMpZvYzMyvehv3OzfDveNXh8tTt3H2tpOsl9ZR0yTacGgAA9Y5kNAAAjV+ZpCXRsFRSa0l7Srpc0kdmdmCMsTUmqyR9KunLWmxTFm3zaZplh0i6ShLJ6PxZrnDtF6VZ9iOF69+3PgOqysx2NbNJkp6RdJ6kgQrfz9UKCbbhkn4g6WVJ08xs37hibUrMrLWk30haJ+m6mMNprG6R9LxConOYpBJJGyV1kzRK0r/N7DEza5brA5vZ3pJ+XCWOzu7eTlKHaFqS/mxmbXJ97AzxNJP0rKTbJO2ncC0KJe0l6U+S3o7uuW3xlSr/HU8d1mXY5h+S1kj6lZm13cbjAwBQb0hGAwDQ+I13927RsINCsutsSSsltZf0mJmVxBphI+Du/3b3ge5+di22WRBtMzCfsSE9d78xuv6/ijuWVFFC7T1JQxUSTb+RNNjdm7l7J0nFkvpJGiNpikKi+oCYwm1qzpLURdIT7r4i7mAaqWaSFkr6f5L2ltTc3dtK6iHppmidb0u6Jg/HPiUaT3H3HyQ+w6iF+w8VErRtJO2fh2On87+SjlVIxo+W1FJSK4WkfKmkfRQS1dvipCr/jqcOt6TbILoej0vqKOmcbTw+AAD1hmQ0AABNjLuvd/f7JF0azeomWucC9SZqsfm4wsOgzyTt7e7XuPv0xDoezHH32919qKTTFFp6Y9tdEI0fjjWKxu1mSf3c/bfu/oG7uyS5+0J3v0TS2Gi9i/PwsHNwNH4mdYG7l0uaE012zvFxv8bMukm6LJr8hbvf4+7l0ff3WYU3HiTpNDPbI9/xpPFQND4/hmMDAFAnJKMBAGi6HpVUEf09LDEzteM6MzvDzN6Iath6ai1QM9vZzG6LOgPbGNW8fdPMLqhJnVsz62Bm11XZfr6Z3W5m3TOsX2hmh5rZP8xskpktqdJ51b/N7LCanHwdjlvrDv0sTQeGiXkKJSIk6Zw0NUD7mtnZ0d+LzKwoyzEOjdZbb2btahpbtO2O0TkviK7BbDP7m9WwQzczO9DMHo6u3aboHnnFzE4zM0uz/iFRrHOj6ZFm9qyZLTezDWY21cwuSbdttH5XC7XOp5nZuijmeWY23sx+b2Z9Utb/WgeGiXmSEuu+lnLtX4/WezWa/ks11+CeaL0Ha3LNIhcp1Csul3SKu39R3Qbu/rC735/m+N3N7CIze87MPovug9UW6tX+LtNnmeaz+Gb02ZVaqHn7spntX2X9dmZ2jZnNjD6reWZ2rWVINKbcywPM7IHoXl4fxXZWlXXNQr3diRZq7ZZG91XvDPtubWanRPucFsW7wcxmRffzLpmuo5kNUWiNvlKh/Enq8jFR3BvNbPcM+7g9WufLmn5XciX6zt4cHXtD9Jn/3cy6RMv3iGJ7LZ9xuPv77r4pyypjo3FLSYNyfPhEyYn5GZb3iMb1UQv8ZEnNFco43Z660N2fkjRTkkk6vR7iSfVfSSsk7WnhbQwAABo8ktEAADRRUSIh0dIybT1JM7te0v0KnUSZKpPXieXHS5om6XsKybWNCq8nHyTpDkkvmlmrLGF0kjRBoX5vN0lbFBIJF0qaambpkhiDJL2q0LJ7qKR2kjYrdF71LUn/NbNfZzlmXY+bK+Xaus7nRn29Bmi5pMcUEhzdJB2TZX+JlndPuPuqmgYRneMHCue8o8I16KZQi3WCwqvd2ba/VtJbkk5VuHabFFr6Hi7pQUkPmlnG/5Y0s9GS3lB4vb1IUgtJe0i6QWnq+EaJ5g8Uap3vppAAWh8de39JVyr7dUpYq3CNE/dyai3W0mj+ndH4TMvwMMBCC+dvR5N31+DYCd+Lxs+5+9RabJfODQqtVI+V1F/hc2ilUK/2t5ImWjUdlZrZDyS9IOlQhVq37SQdodAR3IFRovNtSb9WuFcKFDpG+7nCfZrNvgr30+kKtXRbRLHda2Y/jR48PKBQxmAPhd+ZDgr31Vtm1inNPkcrPEw7XeFeKIiGnRXu5ylmdkSGeL4Zjd9397LUhe5+m0L93+aSHrCUzuei37wLJbmk0fXZ+aGFmuEfKTzM6KlQk76/Qsvc/0a/tYn78an6iiuDquVPct35ZuJ35Wufn5kdK6mXwufzQY6Pm86h0fhNd9+YYZ2XonGNHpTmUtRS/L1o8qj6Pj4AAHVBMhoAgCYqatHYJZpMl1AZptA51VWSOrl7R4Uk0fho+50VXnNvoZBUHOju7RVqdY5RSIododCJUiZXRuuPktTa3VsrdOw3J4otXQdYmxUSYKMUkqcl0XY7RPsrl/S/ZjYix8fNCXef5+7dJCVa3D6SpgboPHffoJDUlaRz0+3LQqdUJ0WTNU6GRuf2L4VznS3p4OgatJb0PwrJyN9m2f4yhUTkMoUO9jpE9WJbSfqOQoeB35X0iwy76KKQfLxFUvfovumgkFiVpEvNbLeUba5SeOAwS9I3JBVH92SJpCEKdVsXV3fu7v6X6PrPi2al1mJNXM/HFRLTOygketM5VaHl5xcKLRCrFSWGd44mn6vJNtX4TKHe9G4K34UOCt/JQxSSwDsre73aLgrJ/z8qfM/bKTxYeifaz98k3apQI/gghe9NG4VSF1skHRclADO5XeH3oV/0ObeP9idJv4+GUQp1nFtH+z5I4bPsrfT30AqFe+UASe2je6+FwoOqBxTuwwczPAgbGY0nZYn5fIXOXvdQuK8kSVFSPvGQ4jp3fzXLPnLKzDpI+rfCQ6JXJO0YnffeCvffEIWHHIlk9JP1FVsGB0fjMoWWwXllZp3M7GJVll55xt0X5vu4qiwZ8nGWdT6JxoOihy91cZ2ZLbPwFtBiM3vezE63Grx9JGliND6ojscGAKB+uTsDAwMDAwNDIxwUXpN2Sa9nWH5JtNwlnVxl/ugq8/+QZf93RevMktQyzfLvRcsrJPVPWfZ6lWUHpdl2gEIy2yWdWcvzvjLa7p9pltX5uFWuy9euZ5X9jk6Z3zdxLdNsc3W0bGyWc9k7WmezpC5ZrvFsSVaLa3RWtN0mSQPSLD+oyj3wesqy9pLWKCSZ9s2w//2ia1yqkDROzD+kyn7vyLDth9Hy36bM/ySaf2otzjPjNZY0N1p2SJbt/xGt8+8My8dHy6+uRUxHVrkG+9Xm3q7toJC4XBoda6eUZVU/i3Tfld7RZ5i4//qnWSfxG3B3mmWJfc+UVJSyrEAhiZ5Y5+ws9+jsWp6zKZTfcEnnpFm+IFr23Wr2Mypar1zhYY0UksGu0Dq5eR0/k7S/BzXY7qdVPouuKctOjJbNi8ZTsx27LkMtY21dJZaH83Bfvx7te7TCWz1fpcT7fuo1qrJt3224DnPT7C9x7B9mifeEKvtoU8tznVtl23WSVqfE9LrCA5ls+zgpWndZrj8LBgYGBgaGfAy0jAYAoAmJarP2NbPLJf0pmv2F0nQEpZCE+Vum/SjUypRCC8H1aVa7UyHxY6psrZfqLXd/K3Wmu3+q0HJXWbbNJHEuI7Osk4/j5py7T5E0WaFV6plpVkm0mB7r7l6LXSfO7YnonFOP+5akNzNse7JCsultd38/Q9zvKiTIO6hKPfIUf8wwP1FeILVe7+ponLamd54kWsEeZ2Zdqy4wswEK5UFc0j9rsc+q5U++SreChTrgi9MME2oTvLuXKnqTIYo1k699Fu7+pULCWJIec/dZabZLtAZPW1s58hd335Ky7wqFUjtSqPv7tVrYVfa9UzWlfrYSfQ8SLc63+g2Ifrd2iCazdgbp7s8otOouUCgp8mOFMkCbFR5UZauXnA+J1ufPuvvSlGXPK5SsSZRjydQqOrUcUG2G2rg1imW1pF/WctvaKlB4QFbVLgodJ6Yrr5Mok1SXYVma/SXuzQ1ZYqz672PrLOul86TCb25nd2/loTV8H4U3ayoUWqA/Ws0+Evd653y98QMAQC5l7CwHAAA0GgdblQ70UiyS9C1335xm2Sx3z5Sw6adQykGSXku3grtXWOgM7gyF2s7pvJ5hvhRe7T893bZRiZHvK7Q4G6yQ9Ez975Yds+y7TseNyZ0KNYHPVZVaylHN50QL5LG13Gfi3N7Iss4bCuUwUh0QjUeYWbayGImkay+Fkg9Vlbr77AzbLYjGHVLmPy9phKRrow7q/iXpXQ/lTPLC3T8ys/cV6h6fqa0fziRqdf/Xa9ABYS01V2XStKq0NWmjWsLfV/hseqoyQVZVpu/DRlUmnVMtlbSrQl34dBJJytTPqqqPsuxbkj6JktOZ9i2FZOO6qgujcic/VCgFtLNCeY/UhjSp59xBlfWL0z4ISPEThZrAu6jys7/St73Od10kytake4i2ycymqvKBQ9pktIfyNHllZr9U+M13SRe6+9x8Hs/dV5pZc4V+AHZXeNB2nkKZoSGqLGOUWH+eQnmnnIeSh33K3X+UZt6Xkn5mZnMk3STpSDM7yt1f+toOgqr3eifVoJwRAABxomU0AACNX5kqW3YtlvS5wmvsP5e0m7tn6uQpXSuwhC5V/l6Qca3Q6jF1/aqybZtYttW2ZtZdoWOqvym0CuuiUG5imcI5JhLo2VpT1vq4MXpAoWXdEDOr2so4kQx9JUpO1Ebi3LLVVM10jRItk0sUEqaZhkQLvJZp9rEmy3ETCdfUFnzXSnpaUrFCnepXJa02s/Fm9jMzS20dmSuJ1tHJut1RndazosnadFwoVXaQKGVI4rr7i+5uiUGhw7y0orcc3o3iG6BQO7lqp4yJ65np+7AkS6v68mi8qJrl2VpbVrdt2uUeOl5L2Gr/ZnawpOkKv2GJTkzXqPKcE63oU8+5eZW/0z2AS41hnUI5o4R3VVnrvb4l7pX5GZYnfse/jN6oqHdmNkaVrex/6u7VtdjNCXff7O6L3P1ldx+jyhIrJ2bpyDJXEg9J0v3OKc2ytTk89i0KZTykcM6ZVH2IVZLD4wMAkBckowEAaPzGe2XnbN3dvb+7H+Xuf3b3bK0Dy7Msq6p59avUSaaOnv6u0FpztsLryx3dvbW7d41a/u2Xp+PGwt1XK3TYKEUJ0ej180TZjtomQ2sq03VI/PfhdVUTplmGsbkIxt03ufsJCq0//6SQGPQq0zPNbM9cHCvFQwoJpN3NbHg07xiFpPxKhTrCtTG9yt97bEtgUSeP1yp8VjcqtJ5t7u4dE995VZadaVD3dV1FZQbuVyh38IpC6/0Sd29f5Zx/klg9ZfOqDwJq+vCiauehuyh9i/X6lLZ1vCqT1LnoFLPWzOwshTc4pFBD/bps6+eTu78o6e1o8pA8Hy7xQC/bmziJZWuVw2R09BApUbqnX5ZVqz70WpGr4wMAkC+U6QAAAOlUbTXdR5Wv3KdK1DDN1Mo62//AJ1rgJrc1s2KF0hySdEZUmzhVTZJFtTpuA3CnpHMknW5mP5X0TYVXzUuVuT5sNssUPpuaXIdUifIJg+tw3G0WfebvSlJUS3iUQkK2t8J12ifHx1trZo9IOl8hMTlRla3SH3T3TMnBTPubb2afK5SWOE6hLnFdnazwcOA/7v7DDOvEnTzNtf0V7t1SSSdkqFef9pyjchZrFRLZ2UqLSJLM7AxJ35W0ReGNkgEKD3+OqVvo26RU4Tuf6Y2NRE3pjG+EVFNWJ6tsJT7M7BSFuukFkv7q7r+r63FyKJEkTn2zppcqE7i1Nc/dU39fPlH4LdwtzfoJid/K6bWs7V8TiQcu2fabuNc3Rw83AQBo0GgZDQAA0pmt0CpUCjVVv8bMClTZKm1yhv0cnOUYiWVVt+2sypbYmV5Fr8lr2bU9bj4k6uRW22LV3d+WNEMhqfAtVbbWfLCOHaklzi1dTeiETNcoUf/5YDPrVIdj54y7r3P3hyV9L5o1rBad3dX4+quyVMdpUTLr+Gi6rq3SEwno47axNXfiYU/a70J0Lbb1TYGGJnHOMzMkoqXsvwGJDjt3ynaQ6HO+MZr8vcL3boOko83sBzWMNZc+icZf6xDUzHqrspb7kCz7yFZWp7ohLTMbpVBKqFDSre5+eY3PKL96R+PSlPmFqvs1SPcgINFnwkFm1iJDLEdG4/9mWF4nUYecibc15mZZtW80/lpntQAANEQkowEAwNdErbueiCYvM7N09TIvkNRDocXWv9Isl0JC84DUmVEHdd+OJh+rsmi1KluAfS3pEtWTztRCdFuOmw+JFmo1LReQSIj+RKFFrVT3ZGji3E6Kznkr0bXJlKh+TKFOagtJf852EDOrtvVpTUWt4jNJdGJoCjWla6LG1z9qjT1N4WHAQwo1jKe6+6QaHivVLZLmKCTGHjOzPnXcz6ponCkBeYVCx35NSeKcd0mX/DOzo5ThAVlkXDQenmmFKMk3VuHeeFfSH9x9hqRfRKv82cx2rWXc2+rFaHyqmbVLWXaxKh+qDI46eP2aGpbVSTuk219Uj/kxhe/DPQq13OvT2VHJotS4dlbl5zux6jJ3n7sN16FvmhieUOizoL3Cv3mpsYxSaFHvCr8dNRbdh9mMUWWiOVt5lkRr7q91fgkAQENEMhoAAGTyB4Wk5I6SnjOzAZJkZs3N7EJJ10fr3eXuszLsY7WkJ8zs2MT/eJvZQZJeUGgB/bGkZCdY7r5WUYkGSXeb2V7RNgVmdrikN1Szlq61Om6efByND0yXEE7jXoVO1/ZVSP5M2YaOyh5RaGnZXNLzZnaglLyOxykkWNK+zu3uKyT9Kpo818weNbPdE8vNrIWZHWhmN6ky8ZcL08zsD2a2TyIxbcG+km6I1plQTR30qhLX/7QsLRqrSjwMGBmN61yr293XKJTYWKlQh3iymf3GzAZVXc/MuprZmZIuy7Crl6PxcWb268RDITPrYmZ/VvicmlqN2HEKHXp2knRv9ABKZlZiZudJelzZzzlRSzhbOZcfSzpM4fftrCqdKd6ocM1bSrovXSI0j/6pcL90kvRk1HJbZnZ8FK8U6pE3l3RlvoMxs5EKJYKaS3pY0nm1KUFhZlebmZvZtpStOFTSy2Y2NNpnYfR78LjCb+Ri5bmGtrsvlvSPaPJPZnZW1MGpzOxYhc9Nkh5y9w9TtzezsdF1mJtm99eb2T+i39OSKtv0MrP/U2XL/dfc/YUsYZKMBgA0KiSjAQBAWu7+uaTTFDrUOkTSDDP7StIahTIEzRVeS/5Rlt38P4UOnZ6TtM7M1kh6U6Ge7jJJ33H3spRtfqzQEnaIpClRDdi1Cp2ZdVKo7Vuduhw3115XqEPbUdKnZrbUzOZGQ8/Uld19maSnq8zalmRomaRTFM61v6S3omuwVtKzCp/h77Nsf4NCwsuj/XxkZuvMrFQhgfeWQivJtC0066irQnL1fUnrzWyFQovE9xQ6AlyuNC0Ts7grGp8iaZWZzYuu/cMZ1r8vOp4UHgo8UMv4txI9SNhPoWRKR4V78hMz22xmy6L7ekl03N0VEo1jUvbxkirfULhG0troM1gi6XKFe+TZbYmzoXH3lap8GHKKpIVmtlLh4cldkmZJylaz+HmF+3tA1IJ2K9GDlT9Ekz+t+iAtSraeK+krhYdCv9m2s6k5d1+u8NtWrvB7+4WZrZb0jELi9W5Vdtz4KzNbYWan5jGk/6fK+tRHKHwOizMM+YrDFa7FJDNbp/Db856kPRV+y06vbU33OvqNwn1VovDQcF0Uz3MK/yZNkPT9Ouy3jaRLFX5P15pZqZmtkvSlQiv9QoUHsN/OtIPo35IhCtcmls4tAQCoLZLRAAAgI3d/RuF/dO9QqFnZUqHV4tsKdXy/6e7rsuxihUKrrb8rJNCKFTqeukPSXu7+SeoG7v6eQidmTyokhZopdKB4m6S9JE2tQei1Pm6uRQnhwxWSjQsUSkD0iYZMLS4TicdN2vZk6CcK1+tOSYtU2ZLwOoVrk1prNXX7/1VI+twu6TOFFumton29IOkiSSO2JcYUJ0j6o0LL2IUKndBtlvShpP+TtFu6lodZ4n9V0okKyZwNCiVl+ih0Epdu/dJoXUl6Kmohvk3c/VOFcgKjFFpQfqrw/Wmv8JBnskJJjyMUzu8/aXZzqqRfKiSryxQ+h3GSznH3mjyYaXTc/XpJJ6mylXSRQk31qxRqJ6/Jsu06hTcDpNA5YVLU4v5+hQdpz7n7bWm2X6BQFkOSroha4tYLd39C4bfvCYUHSS0VfgPvlDTG3V9UeCCzUOEBRy4fBqWq+v+JnZW91nK6OBIdpE5Ms6ymrlD4/ZmrcA+4wm/RrQq/469l3jR3ot/yUQoJ53cVfp9d0gcKSeMDo7chautWSX+RNF7hM22hcG/Ok/RvSd+RdFj025TJdxV+Ex6tYwwAANQ7y32HvwAAAKgLM7tDIdn0iLt/t7r1kTtRCYxFktpKOiZK/KERMrPhCq1VP3b33atbP8fHdinUb67P4zY0ZjZDoZbyKHevVet9M3tdoYPVc919bO6jazrMbJKkoZIOcPd3qlsfAICGgJbRAAAADUDUaVkiAX17nLFsp05TSER/IemlmGPBNnD3iQrlS3Yzs6Pjjmd7Y2Y7KCSiJ9c2EY2aM7NDFRLR/yERDQBoTEhGAwAAxCwqH/A3hdIUH0qql9fPEZhZX0lXR5PXu3tFbMEgV34lqUKh1APq1zeicca69MiJ3yiUC/lVdSsCANCQ1GcP0QAAAKjCzL6tUDO0s0I9ZlfoVI06avUg6szwQIX6tgWSZkq6OdagkBPuPs3MLpDUx8zaRx0joh64+2MKdYyRJ2bWWqFT3oejzlIBAGg0SEYDAADEp7VCp3qbJE2R9Ht3fyXekLYr3RQ6NixVaI3+U3ffGG9IyBV3/2fcMQD54O5rJf0u7jgAAKgLOjAEAAAAgBygA0MAAIDsSEYDAAAAAAAAAPKODgwBAAAAAAAAAHlHMhoAAAAAAAAAkHckowEAAAAAAAAAeUcyGgAAAAAAAACQd0VxB9AUmdkcSW0lzY05FAAAAAAAAADYFn0lrXb3nbZ1RySj86NtSUlJx0GDBnWMOxAAAAAAAAAAqKvp06drw4YNOdkXyej8mDto0KCOkyZNijsOAAAAAAAAAKizYcOGafLkyXNzsS9qRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAAAAAAAAA8o5kNAAAAAAAAAAg70hGAwAAAAAAAADyjmQ0AAAAAAAAACDvSEYDAAAAAAAAAPKOZDQAAAAAAAAAIO9IRgMAAAAAAAAA8o5kNAAAAAAAAAAg74riDgAAamtjWbm+WLFec5av1ezl67RyfZkqKlwuyV2qcJckFRcVqEf7EvXqWKJeHVqqZ4eWKikujDd4AAAAAACA7RTJaAANWll5hSbMKdVrny7VjMVrNHvZOi1ctUFRvrnWOrdurj6dWmp4nw7ar18nDe/bQW1aNMtt0AAAAAAAAPgaktEAGpxV68v0+sylemX6Ur3+6VKt2bglZ/tevnaTlq/dpElffKXb3pytApOG9Gin/XbupP36ddL+/TqpRTNaTwMAAAAAAOQayWgADUJFheu/M5Zq7Pg5end2qcorMjd9LjCpV8eW2qlzK+3UuZV2aNtCBSaZTGaSmckkrd+8RfNKN2j+yvWaV7pBC1du0JaU/Va4NHX+Kk2dv0q3vTFbbVsU6bg9dtRJQ3toeJ8OMrM8nzkAAAAAAMD2gWQ0gFiVlVfo6Q8W6tY3PtdnS9emXadH+xIdMairDujfWTt3aa3eHVuquKj2/a9uKa/Q4tUb9cnC1XpvTqne+XyFpi9evVXJj9Ubt+ih97/UQ+9/qd4dW+pbe/fQSXv3UN/Orep6igAAAAAAABDJaAAxWb95ix6ZME93vjVHC1Zu+NryPXu11xEDu+qIwTtoYLc2OWmhXFRYoJ5RR4ZH7dZNkrRy/Wa9P6dU78xeoVemL9G80spYvixdr+v/+5mu/+9nGtm/ky46uL9G9u9Ea2kAAAAAAIA6IBkNoF5VVLgemvCl/vKfT/XV+rKtlrUqLtSZ+/XROQf01Y7tS+olnvYti3XUbt101G7d9NvjB2viF1/pickL9OyHC7eqVT1u1gqNm7VCe/Rsp4sO3llH7dZNhQUkpQEAAAAAAGqKZDSAejNn+Tr98vEP9d7+uFgOAAAgAElEQVSc0q3md2pVrPMO3Elnjuijdi2bxRRdqDW9T9+O2qdvR101arBenbFUT0xeoNc+XZqsYf3h/FW66IHJ6te5lcYc3E/f2ruHmhfR4SEAAAAAAEB1SEYDyLst5RW68+05uu7lmdq0pSI5v0f7Eo05uJ++M7yXWjRrWAndFs0KdeyQ7jp2SHfNK12v29+crUcnzkvGP3v5Ov3i8Y90w6uz9MtjBuq4Id0p3wEAAAAAAJCFedWeu5ATZjZp6NChQydNmhR3KEDsPl64Sr94/ENNW7A6Oa+wwDTmG/106eG7NLgkdDbL1mzS2PFzdO87X2xVwkOShvfpoCuPH6w9e7WPKToAAAAAAIDcGzZsmCZPnjzZ3Ydt675oGQ0gL9xdt7zxuf760sxkiQtJ2m3Htrr25D20e492MUZXN13aNNfPvjlQYw7eWQ+8+6Vuf/PzZN3riV98pRNuGqeT9u6hnx09QN3b1U/NawAAAAAAgMaiIO4AADQ9G8vK9ZNHp+pPL36aTEQXFxXoF0cP1JMXj2yUieiq2rZoposO2Vmv/+xQXXjQTmpWWFme44kpC3ToX17XP175TJurlCQBAAAAAADY3pGMBpBTy9du0hl3vqd/T1mQnDesTwe9eNlBuuiQndWssOn87LQraaYrjhusl398sI4avENy/sayCl33ykyNuuFtTZ23MsYIAQAAAAAAGo6mkxUCELsZi1frhBvHadIXXyXnnbZvLz104X7q16V1jJHlV9/OrXT72cP14IUjNLh72+T8T5es0Yk3j9Mfn5+ujWXlMUYIAAAAAAAQP5LRAHLi1RlLdPLN47Vg5QZJUoFJVx4/WH84cYiKi7aPn5oDdu6sZ354oK4aNVglUceMFS7d9uZsHfOPtzRhbmnMEQIAAAAAAMRn+8gQAcire8bP1QX3TNS6zaH1b+vmRbrznOE6/8CdZGbVbN20FBaYzh25k1768Td0wM6dkvPnLF+n79z2jq5++mNaSQMAAAAAgO0SyWgA2+Ted+bqqqc/VtRPoXq0L9HjFx2gwwbukHW7pq5Xx5Z64IIR+uNJQ9S6eZEkyV0aO36uvnXTOM1auibmCAEAAAAAAOoXyWgAdfboxHn67VMfJ6f37t1eT10yUgO6tYkxqobDzHTavr310o+/oUMHdEnOn7F4jUbdME6PTpgnd48xQgAAAAAAgPpDMhpAnTw9daF++fiHyem9erXXfeePUOfWzWOMqmHasX2J7h69j/5w4hA1j+pnbygr188f/1A/euQDrd20JeYIAQAAAAAA8o9kNIBae+njxfrxIx8kS3MM7t5W95y7b7IcBb7OzHT6iN566pKR6t+1dXL+Ux8s1PHXv6WP5q+KMToAAAAAAID8IxkNoFbemLlMlzw4ReVRJnqXrq113/n7ql3LZjFH1jgM7NZWT18yUqcO75WcN3fFep18y3g9/P6XMUYGAAAAAACQXySjAdTYu7NXaMx9E7W5vEKS1KdT6KSvE6U5aqVlcZGu/fYe+sd390q2Jt9cXqFfPvGRfvvUNJVF1xcAAAAAAKApIRkNoEY+X7ZWF94zURvLQqK0R/sSPXDBCHVt2yLmyBqvE/bqoWd/eKAGdW+bnHfvO1/orLve04q1m2KMDAAAAAAAIPdIRgOo1tpNWzTmvklaE3W017VNcz1wwQj17NAy5sgav76dW+nxi/bXcXt0T857d3ap/ufGcfpk4eoYIwMAAAAAAMitnCWjzaynmd1tZgvNbJOZzTWzv5tZh1rs40gz+6uZ/dfMSs3MzeztLOtfHa2Tbfg8ZZtDqln//7blOgBNjbvrZ49N1aylayVJzYsKdPfofdS3c6uYI2s6WhYX6cbT9tbPvjlAZmHegpUbdPIt4/Xch4viDQ4AAAAAACBHinKxEzPbWdJ4SV0lPSVphqR9JV0m6WgzG+nuK2qwq4slnSBpo6RZkqpLZL+eZdkoSUMlvZBh+RsZts+Y/Aa2R7e88blemLY4Of3Hk4Zo9x7tYoyoaTIzXXxofw3s1kaXPfyB1m7aog1l5br4wcn6fNmu+uFh/WWJTDUAAAAAAEAjlJNktKSbFRLRl7r7DYmZZvY3ST+WdI2k79dgP9dKukIhmd1L0pxsK7v760qTUDazQknnR5O3Z9j8dXe/ugYxAdutN2cu01/+82ly+pz9++ikoT1jjKjpO3zQDnry4gN04b2TNGf5OknS316eqYUrN+j/fWt3NSukuhIAAAAAAGictjmrYWb9JB0laa6km1IWXyVpnaSzzKzad/rd/R13/9jdy7cxrGMl9ZT0rrt/uI37ArZL80rX69KHp6jCw/Q+fTvoN8cPjjeo7UT/rm305MUjdWD/zsl5D0+Ypwvumah1Ud1uAAAAAACAxiYXTewOi8YvuXtF1QXuvkbSOEktJe2Xg2PV1PeicaZW0ZLU38wuMbNfm9l5ZrZLfQQGNAYbNpdrzH2TtHJ9maTQYeFNZwylVW49alfSTHeP3kcnDe2RnPfGzGU69fZ3tHTNxhgjAwAAAAAAqJtcZJYGROOZGZZ/Fo13zcGxqmVmPSQdI2mVpEeyrHqGpBsUSojcJWmmmf2rNh0uAk2Ru+uKf3+kTxatliQ1KzTdcuYwdW3TIubItj/FRQX66yl76tLD+ifnTVuwWifeNF6zlq6JMTIAAAAAAIDay0UyOtGT2aoMyxPz2+fgWDVxgaRCSfe7+/o0y5dJ+qWkIZLaSOqikLyeIulkSc+YWY2ui5lNSjdIGpiLEwHi8PTUhXpiyoLk9FWjdtOwPjyjiYuZ6SdHDdAfTxqiwoLQgeGClRt08i3vaMLc0pijAwAAAAAAqLn6eOfeorHn/UAhiXxeNJm2REdUk/pad5/m7mvdfbm7vyjpEIUOE0dKGpXvWIGGaOmajbrq6Y+T098e1lNnjOgdY0RIOG3f3rrz7OEqaVYoSVq1oUxn3fWe3pi5LObIAAAAAAAAaiYXyehEy+d2GZa3TVkvn46R1Ft16LjQ3VdLejCa/EYNtxmWbpA0o1ZRAw1AKM8xLVknukf7El39P7vJzKrZEvXl0IFd9ciY/dS5dbEkaWNZhS64Z4Je+GhRzJEBAAAAAABULxfJ6E+jcaaa0ImOATPVlM6lRMeFt9Vx+0QTw1Y5iAVoVJ76YKFe/mRJcvpP395DrZsXxRgR0tmjZ3s99v0D1KN9iSSprNx18YOT9a9J82OODAAAAAAAILtcJKNfi8ZHpdZaNrM2CmUvNkh6NwfHysjMdpR0nEIL7EfruJv9ovHsnAQFNBJLV29dnuOMEb01sn/nGCNCNjt1bqVHv7+/+nUOz80qXLr8sakaO25OzJEBAAAAAABkts3JaHf/XNJLkvpKujhl8e8UWhnf6+7rEjPNbKCZ5bqTv/MVOi68L0PHhYljj0zXQaGZnSnpVEmbVfdkNtDouLt+/e9pWrWhsjzHr44dFHNUqE6P9iV6ZMz+GtS9bXLe1c98ohtf/UzueS/RDwAAAAAAUGu5egf/B5LGS7rezA6XNF3SCEmHKpTnuCJl/enReKtitGZ2oKQLosnW0XgXMxubWMfdR6cePEounx9Npu24sIoHJBWY2XhJ8yW1kLSPpH0lbZE0xt3nVrMPoMl46oOFemU65Tkaoy5tmuvhC/fTuWPf1+QvV0qS/vLSTK3dVK5fHD2Aet8AAAAAAKBByUnGyd0/N7Phkn4v6WhJx0paJOl6Sb9z99Ia7qq/pHNS5nVNmTc6zXbflNRHoePCj6o5xi2SjlAoH9JZISG+QNJYSX9396k1jBVo9CjP0fi1a9lM950/QmPum6S3Zy2XJN36xudyuX559EAS0gAAAAAAoMHIWfNHd58n6dwarps2O+LuYxWSwrU99gtKaWWdZd1rJV1b22MATQ3lOZqOVs2LdOc5w3XJg1OSrdxveyOUvichDQAAAAAAGopcdGAIoBF6cdrircpz/JnyHI1ai2aFuvmMoTpy8A7Jebe9MVv/9+IMakgDAAAAAIAGgWQ0sB3aWFaua56fnpw+Y0RvHUB5jkavuKhAN51OQhoAAAAAADRMJKOB7dDd4+Zo/lcbJEntWzbTz785MOaIkCuJhPRRJKQBAAAAAEADQzIa2M4sXbNRN706Kzn9kyN3VbuWzWKMCLlWXFSgG0lIAwAAAACABoZkNLCd+ct/PtW6zeWSpF26ttbp+/aOOSLkQ6aE9PX/nZVlKwAAAAAAgPwhGQ1sR6YtWKXHJs1PTl95/GAVFfIz0FSlS0hf98pM3fHm7BijAgAAAAAA2yuyUMB2wt31+2c/UaJKw2EDu+obu3aJNyjkXXFRgW44fe+tPutrnp+u+9/9IsaoAAAAAADA9ohkNLCdeGHaYr0/p1SSVFRguuK4QTFHhPrSvKhQt505TPvu1DE578qnpunfU+Zn2QoAAAAAACC3SEYD24GNZeX6w/PTk9Nn799XO3dpHWNEqG8lxYW6e/Q+2rNXe0mSu3T5Yx/qxWmLYo4MAAAAAABsL0hGA9uBu96eo/lfbZAkdWjZTJcdvkvMESEOrZsX6Z5z99HAbm0kSeUVrh8+NEWvfbo05sgAAAAAAMD2gGQ00MQtXbNRN782Kzn9kyN3VbuWzWKMCHFq37JY950/Qv26tJIklZW7Lrp/kibOLY05MgAAAAAA0NSRjAaauJtenaV1m8slSbt0ba3T9u0dc0SIW5c2zfXABSPUs0OJJGljWYXOGztBny5eE3NkAAAAAACgKSMZDTRhS1Zv1EMT5iWnf33sIBUV8rWH1L1diR64YIQ6t24uSVq9cYvOvvs9zStdH3NkAAAAAACgqSIrBTRht7z+uTZvqZAk7dmznQ4Z0CXmiNCQ9OnUSmPP3UetmxdJkpas3qSz735fy9duijkyAAAAAADQFJGMBpqopas36qH3v0xOX3bELjKzGCNCQ7R7j3a64+zhKo5azM9Zvk7n/nOC1m7aEnNkAAAAAACgqSEZDTRRt74xW5uiVtF79GynQwd0jTkiNFT779xJ15+2lwqiZxUfLVilMfdN1KYt5fEGBgAAAAAAmhSS0UATtHTNRj3w3hfJ6UsPo1U0sjt69+665sQhyelxs1boJ49MVXmFxxgVAAAAAABoSkhGA03QHW9WtorevUdbHT6IVtGo3mn79tblR+2anH7uo0W65rnpMUYEAAAAAACaEpLRQBOzfO0m3fcuraJRNxcf2l+jD+ibnL573Bzd+dbs+AICAAAAAABNBslooIm5483Z2lgWWkUP7t5WRw7eIeaI0JiYma48frCO2b1bct41z0/X8x8tijEqAAAAAADQFJCMBpqQFWs36d53qrSKPpxW0ai9wgLTdafupWF9OkiS3KUfPfKBJswtjTkyAAAAAADQmJGMBpqQO96aow1l5ZKkgd3a6ChaRaOOWjQr1J1nD1e/zq0kSZu3VOjCeydq1tK1MUcGAAAAAAAaK5LRQBNRum6z7n1nbnL60sN3UUEBraJRdx1aFeue8/ZV59bFkqSV68s0+p/va+majTFHBgAAAAAAGiOS0UAT8c9xc7R+c2gVPWCHNjp6t27VbAFUr1fHlrp79D4qaVYoSZr/1QadP3ai1m3aEnNkAAAAAACgsSEZDTQBG8vK9cB7XyanLzmsP62ikTN79Gyvm87YW4lb6qMFq3TZwx+ovMLjDQwAAAAAADQqJKOBJuDpDxaqdN1mSVKP9iU6ZndaRSO3Dhu4g/73W0OS069MX6I/Pj89xogAAAAAAEBjQzIaaOTcXXePm5OcPnv/Pioq5KuN3Dt9RG+N+Ua/5PSdb8/R/e9+EWNEAAAAAACgMSFjBTRy784u1YzFayRJJc0K9d19esccEZqyXxw9UEcN3iE5fdXTH+uNmctijAgAAAAAADQWJKOBRu6fVVpFnzS0h9q1bBZjNGjqCgpMf//uXhrSo50kqbzCdckDk/Vp9EAEAAAAAAAgE5LRQCP25Yr1enn6kuT0uSP7xhcMthsti4t01znD1b1dC0nSmk1bdN7YCVq2ZlPMkQEAAAAAgIaMZDTQiN3zzly5h78P2qWz+ndtE2s82H50bdtCd4/eR62KCyVJC1Zu0AX3TtTGsvKYIwMAAAAAAA0VyWigkVq7aYsenTAvOX3eyJ1ijAbbo0Hd2+rG04eqwML01Hkr9fN/fShPPCEBAAAAAACogmQ00Eg9Pmm+1mzaIknq17mVDt61S8wRYXt06MCuumrUbsnpp6cu1C1vfB5jRAAAAAAAoKEiGQ00QhUVrrHj5yanR4/sq4JE81Sgnp1zQF+dMaJ3cvrP//lUr3yyJMsWAAAAAABge0QyGmiE3pi5THOWr5MktWlRpJOH9ow5Imzvrhq1m0bs1FGS5C5d9vAUzVyyJuaoAAAAAABAQ0IyGmiE7h43J/n3qcN7qVXzohijAaTiogLdfMZQ9exQIklat7lcF9wzUV+t2xxzZAAAAAAAoKEgGQ00MrOWrtFbny2XJBVYKJEANASdWjfXHWcPV8viQknSl6XrdfGDk1VWXhFzZAAAAAAAoCEgGQ00MlVrRR8xaAf16tgyvmCAFIO6t9XfvrNXcnr85yt0zXPTY4wIAAAAAAA0FCSjgUZkw+ZyPTVlYXJ69Mi+8QUDZHD07t30kyN3TU6PHT9XD73/ZYwRAQAAAACAhoBkNNCIvPjxIq3ZtEWS1LdTS+3fr1PMEQHp/fCw/jpuSPfk9JVPTtN7s1fEGBEAAAAAAIgbyWigEXl0wvzk36cM7yUzizEaIDMz059P2UODu7eVJG2pcF30wGTN/2p9zJEBAAAAAIC4kIwGGokvV6zXO1HL0gKTTh7aM+aIgOxaFhfpjnOGq1OrYklS6brNuuCeiVoXte4HAAAAAADbF5LRQCPxr0nzkn8fvGsXdWvXIsZogJrp0b5Et541TM0KQyv+GYvX6PLHpqqiwmOODAAAAAAA1DeS0UAjUF7h+tekyhIdp+7TK8ZogNrZp29H/e+3dk9OvzBtsa5/9bMYIwIAAAAAAHEgGQ00AuNmLdfCVRslSR1bFeuwgTvEHBFQO6fu01ujD+ibnP77K5/phY8WxRcQAAAAAACodySjgUbg0YmVJTpO3LuHiov46qLx+c1xg3Rg/87J6Z88OlWfLFwdY0QAAAAAAKA+kdECGriV6zfrpY+XJKe/M5wSHWicigoLdOPpe6tPp5aSpA1l5Rpz/0StXL855sgAAAAAAEB9IBkNNHBPTlmgzeUVkqQ9e7bTgG5tYo4IqLv2LYt11znD1bp5kSRpXukGXfrwByqnQ0MAAAAAAJo8ktFAA/foxMqOC0+hVTSagP5d2+gvp+yZnH5z5jL9/ZWZMUYEAAAAAADqA8looAGbtmCVPlkUauo2LyrQqD13jDkiIDeO3r2bLj505+T0Da/O0ksfL44xIgAAAAAAkG8ko4EG7LEqHRceO6S72pU0izEaILd+cuQAHbTL1h0afr5sbYwRAQAAAACAfCIZDTRQG8vK9eQHC5PTpwzvGWM0QO4VFpiu/+7e6tmhRJK0dtMWff++SVq7aUvMkQEAAAAAgHwgGQ00UC9/skSrNpRJknp1LNF+O3WKOSIg9zq0KtatZw5T86Lwz9FnS9fq5/+aKnc6NAQAAAAAoKkhGQ00UI9NqtJx4bBeKiiwGKMB8mf3Hu30x5OGJKef/2ixbn9zdowRAQAAAACAfCAZDTRAy9du0rhZy5PTJw+jRAeatpOG9tQ5+/dJTl/74oytvgMAAAAAAKDxIxkNNEAvTFus8opQpmCfvh3Uo31JzBEB+XfFcYM1vE8HSVKFS5c8OFnzv1ofc1QAAAAAACBXSEYDDdCzUys7Ljx+jx1jjASoP8VFBbr5jKHq0qa5JOmr9WW66P7J2lhWHnNkAAAAAAAgF0hGAw3MktUb9f7cUklSgUnHDOkWc0RA/enatoVuOWOoiqIa6R8tWKXfPDmNDg0BAAAAAGgCSEYDDczzHy1SIu+2X79O6tqmRbwBAfVseN+O+u2owcnpf02arwfe+zLGiAAAAAAAQC6QjAYamGc/XJT8mxId2F6dtV8fnTS0R3L6d898rElffBVjRAAAAAAAYFuRjAYakAUrNyQTboUFpqN3p0QHtk9mpj+cOESDu7eVJJWVu37wwCQtXbMx5sgAAAAAAEBdkYwGGpDnPqzsuHBk/87q2Ko4xmiAeLVoVqjbzhqm9i2bSZKWrN6kSx6coi3lFTFHBgAAAAAA6oJkNNCAbF2io3uMkQANQ6+OLXX9d/eWhf4M9f6cUv35pU/jDQoAAAAAANQJyWiggfhixTp9OH+VJKlZoembu1GiA5Ckb+zaRT89ctfk9G1vzNZLHy+OMSIAAAAAAFAXJKOBBqJqq+iDd+2idiXNYowGaFh+cEh/HTqgS3L6p49N1Rcr1sUYEQAAAAAAqC2S0UAD8czUynrRx++xY4yRAA1PQYHpulP3Uo/2JZKkNRu36KL7J2tjWXnMkQEAAAAAgJrKWTLazHqa2d1mttDMNpnZXDP7u5l1qMU+jjSzv5rZf82s1MzczN6uZhvPMrybZbvjzex1M1tlZmvN7D0zO6c25wzkyqylazVj8RpJUvOiAh0xeIeYIwIanvYti3XzGUNVXBj+6fpk0Wpd/fTHMUcFAAAAAABqqigXOzGznSWNl9RV0lOSZkjaV9Jlko42s5HuvqIGu7pY0gmSNkqaJammiewvJI1NM39+hngvkXSDpBWS7pe0WdK3JY01syHufnkNjwvkxLMfVraKPnRAV7VunpOvJtDk7Nmrva4cNVhXPjlNkvTwhHka1qeDThneK+bIAAAAAABAdXKV8bpZIRF9qbvfkJhpZn+T9GNJ10j6fg32c62kKxSS2b0kzanh8ee6+9U1WdHM+kr6i6RSScPdfW40//eSJkj6qZk97u7v1PDYwDZx961KdIzakxIdQDZnjuitiXNL9dQH4Xvzmyenabcd22nwjm1jjgwAAAAAAGSzzWU6zKyfpKMkzZV0U8riqyStk3SWmbWqbl/u/o67f+zu+SwCep6k5pJuTCSio2N/JekP0WRNEudATsxYvEafLwsdsbUsLtRhA7vGHBHQsJmZ/nDiEO3StbUkadOWCv3ggUlavbEs5sgAAAAAAEA2uagZfVg0fsndK6oucPc1ksZJailpvxwcK5P2Znaemf3azC42s2zHSsT7YpplL6SsA+Rd1RIdhw/aQSXFhTFGAzQOrZoX6ZYzh6pl9H2Zu2K9fv7Yh3L3mCMDAAAAAACZ5CIZPSAaz8yw/LNovGsOjpXJnpLuUigHcqOkd8zsAzMbkmbdjPG6+yKFltw9zaxldQc1s0npBkkD63wm2K64u577cFFy+vg9uscYDdC49O/aRv938h7J6Rc/Xqy73q5pdScAAAAAAFDfcpGMbheNV2VYnpjfPgfHSudvkkZK6iKpjaR9JP1LIUH9qpn1SFm/pvG2y7AcyJlZS9dq7or1kqRWxYU6eNcuMUcENC7/s+eOOmf/PsnpP74wQxPmlsYYEQAAAAAAyCQXyejqWDTOy7vT7v5Tdx/v7svdfa27T3T3UyQ9LqmzpMtrucsax+vuw9INCh0wAtV6efqS5N+HDOiqFs0o0QHU1hXHDdZevcLzzvIK1yUPTtbytZtijgoAAAAAAKTKRTK6upbEbVPWqy+3RuNvpMyvabyrcx4RkOKVTyqT0UcMpuNCoC6Kiwp00xlD1aFlM0nSktWbdOlDU1ReQf1oAAAAAAAaklwkoz+NxplqQu8SjTPVlM6XZdG4Vcr8jPGaWfdo/fnuvj6PsQFatmaTpsxbKUkqLDAdOoBkNFBXPdqX6LpT95JF77aM/3yFrnu5vv/ZAQAAAAAA2eQiGf1aND7KzLban5m1UajnvEHSuzk4Vm3sF41np8x/NRofnWabY1LWAfLm1RlL5FHDzeF9Oqh9y+J4AwIauUMGdNUPD9slOX3ja7P06owlWbYAAAAAAAD1aZuT0e7+uaSXJPWVdHHK4t8ptDS+193XJWaa2UAzG7itxzazoWaW2vJZZraHpGuiyftTFv9T0iZJl5hZ3yrbdJD062jyVgF59vInS5N/Hzl4hxgjAZqOyw7fRQft0jk5/eNHpmr+V7zoAgAAAABAQ1CUo/38QNJ4Sdeb2eGSpksaIelQhfIcV6SsPz0aW9WZZnagpAuiydbReBczG5tYx91HV9nkUkknmdmrkuYpJJkHKrR6LpR0h6SHqh7D3eeY2c8kXS9popk9ImmzpG9L6inpr+7+Ti3OHai1DZvL9fasZcnpwweRjAZyobDA9PdT99LxN7ytRas2atWGMl3y4BQ9OmZ/FRfVR5+9AAAAAAAgk5z8n3nUOnq4pLEKSeifStpZIeG7v7uvqOGu+ks6JxpOjuZ1rTLvnJT1n5T0iqTdo2WXShom6QVJJ7j799z9az1YufsNkv5H0seSzpb0PUmLJY1298trGCtQZ+NmLdfGsgpJUv+urbVT56818AdQR51aN9eNp++twoLwvPODeSv1pxdnxBwVAAAAAADIVctoufs8SefWcF3LMH+sQkK7psd8UiEhXWvu/oykZ+qyLbCtXpleWcf2CFpFAzk3rE9H/fybA/THF0IS+s6352jfnTrqqN26xRwZAAAAAADbL95ZBupZRYXrlelV60V3jTEaoOm68KB+Onxg5ffr8semal4p9aMBAAAAAIgLyWignk2dv1LL126SpP/P3n3HV13dfxx/n5sdEhIg7L2HjDAUAQfOorbVVtwizlpx/bTa2tpWbbVVq9Xi3hPcrW3dWrcsIey9NwQSkhAy773n98e9fHMTCQS4yfeO1/PxyON7z7n3+/1+aA0J73vu56hVs2Tldm7hckVAbPJ4jB44Z4g6ZKVKkkoqvLrutbmq8vpdrgwAAAAAgPhEGA00sdAWHSf2a+P0tQUQfi2aJeuRC4cpMfh9Nn9jke79kP7RAAAAAHNvuOcAACAASURBVAC4gTAaaGKfLalp0XHyAPpFA41teNcW+s24fs74+e/W6uPF21ysCAAAAACA+EQYDTShDQVlWr59tyQpOdGjY3vnuFwREB+uPLa7Tu5P/2gAAAAAANxEGA00odAWHcf0ylF6cqKL1QDxw5hA/+iO2WmSpN0VXl03NY/+0QAAAAAANCHCaKAJhYbRJ/enRQfQlLLTk/XohUNr+kdvKtZfPljqclUAAAAAAMQPwmigiRSXVWvm2kJnfFJIywAATWNolxa67bSa/tEvTlunjxZtdbEiAAAAAADiB2E00ES+XJEvn99KkoZ0ylLb5qkuVwTEpyuO6a5TQjYPvfXtBdpQQP9oAAAAAAAaG2E00EQ+XUKLDiASGGP0wPgh6tSipn/0tVPzVOn1uVwZAAAAAACxjTAaaALVPr++Wr7DGZ88gDAacFNWepIevXCYkhIC/aMXbi7WX96nfzQAAAAAAI2JMBpoAnM3FGl3pVeS1DE7Tf3aZbpcEYDcztn67Wn9nfFL09frg4X0jwYAAAAAoLEQRgNN4JuVNauij+uTI2OMi9UA2OuyMd30oyNqPqnwm7cXaH3BHhcrAgAAAAAgdhFGA03g65U7ncfH9W7tYiUAQhljdP/4IercMtg/utKrSVPyVFFN/2gAAAAAAMKNMBpoZLv2VGnBpiJJksdIo3vmuFwRgFBZaUl6LKR/9OItJbqH/tEAAAAAAIQdYTTQyL5bvVPWBh4P6ZytrPQkdwsC8AODO2Xr9tNr+ke/MmO93luwxcWKAAAAAACIPYTRQCP7ZkVNi45jadEBRKyJo7vptIHtnPFt7yzU2p30jwYAAAAAIFwIo4FGZK2tvXlhb1p0AJHKGKP7xg9Wl5bpkqTSSq+upX80AAAAAABhQxgNNKLVO/ZoS3GFJCkzJVG5nbNdrgjA/jRPDfSPTk4I/HhcsrVEf35victVAQAAAAAQGwijgUYUuip6dK9WSkzgWw6IdIM6Zen3P67pHz1l5gb9Zz79owEAAAAAOFwkY0Aj+mYl/aKBaDTh6K46Y1B7Z/zbdxZozY5SFysCAAAAACD6EUYDjaTS69P01QXO+DjCaCBqGGP017MHqWurQP/oPVU+XTt1Lv2jAQAAAAA4DITRQCOZs36XyoPBVddW6eoSDLUARIe6/aOXbi3R3e/TPxoAAAAAgENFGA00ktotOnJcrATAoRrYMUt/COkf/eqMDXpvAf2jAQAAAAA4FITRQCMJ3byQFh1A9Lr46K46fVA7Z3zbOwu1buceFysCAAAAACA6EUYDjaCgtFKLNpdIkhI8RqN6tnK5IgCHyhije88erC4tA612Siu9uu61PFV66R8NAAAAAMDBIIwGGsG3q2padAzrkq3M1CQXqwFwuJqnJunRC4cqKcFIkhZtLtFf3l/qclUAAAAAAEQXwmigEXy9IrRfNC06gFgwuFO2fnd6Tf/ol6av1wcLt7pYEQAAAAAA0YUwGggza22tftFsXgjEjktHd9OPjmjrjH/z9gJtKChzsSIAAAAAAKIHYTQQZiu2lyp/d6UkKSstSYM7ZbtcEYBwMcbo/vFD1KlFmiRpN/2jAQAAAABoMMJoIMxCV0Uf0ytHCR7jYjUAwi0rLUmPXjjM6R+9YFOx7v1wmctVAQAAAAAQ+QijgTD7emVov2hadACxKLdztn4zrp8zfuG7dfp48TYXKwIAAAAAIPIRRgNhVFHt08w1Bc74GMJoIGZdcUx3ndy/pn/0rW/N18ZC+kcDAAAAAFAfwmggjPI27FKl1y9J6pHTTJ1apLtcEYDGYozRA+cMVsfsQP/okgqvrnttrqqCfwcAAAAAAIDaCKOBMJqxptB5PKpnKxcrAdAUstOT9ciFQ5UY7A0/f2OR7v+I/tEAAAAAAOwLYTQQRqEtOo7uQRgNxINhXVro1+P6OuNnv12rz5Zsd7EiAAAAAAAiE2E0ECYV1T7N3VjkjEf2aOliNQCa0pXH9NCJ/do441+9NV+bi8pdrAgAAAAAgMhDGA2EybyNRU6v2B45zdQmM9XligA0FY/H6MFzhqh9VuD7vri8WtdNzVO1j/7RAAAAAADsRRgNhMnMkH7RI2nRAcSdFs2S9eiFQ5UQ7B89d0ORHvh4uctVAQAAAAAQOQijgTCZuTa0XzQtOoB4NLxrS91yak3/6Ke+XqPPl9E/GgAAAAAAiTAaCItKr09z1u9yxiO7szIaiFdXH9dDY/u2dsa/enO+ttA/GgAAAAAAwmggHBZsKlZlsF90t1bpapdFv2ggXu3tH92ueeDvgV1l1brhtbn0jwYAAAAAxD3CaCAMZq6padHBqmgArTJSNPmCmv7Rs9fv0t8/XeFyVQAAAAAAuIswGgiDGbU2L6RfNADpqO4tdfMpfZzxE1+u1pfL812sCAAAAAAAdxFGA4ep2uev3S+6ByujAQRcc3xPHds7xxnf/OZ8bSuucLEiAAAAAADcQxgNHKYFm4pVXu2TJHVumaaO2WkuVwQgUng8Rg+dl6s2mSmSpMI9Vbrh9bny0j8aAAAAABCHCKOBwzSDftEA9iMn2D862D5as9YW6uHPVrpbFAAAAAAALiCMBg7TzLU1/aKPpkUHgH04ukcr/d/JNf2jH/tylb5ZucPFigAAAAAAaHqE0cBhqPb5NWddyOaF3dm8EMC+XXtCL43pFXjDylrp/16fp/wS+kcDAAAAAOIHYTRwGBZtLtaeqkC/6I7ZaercMt3ligBEqgSP0cPnDVVORqB/dMGeKt34+jz5/NblygAAAAAAaBqE0cBhCG3RMbIHq6IB7F/rzBRNPj9XJtg/evqaAk3+H/2jAQAAAADxgTAaOAwzQzYvPJrNCwE0wOheObrhxN7OePLnK/Xdqp0uVgQAAAAAQNMgjAYOkdfn1/frdjljVkYDaKgbTuqtUT1q+kff+Po85e+mfzQAAAAAILYRRgOHaMnWEpVWeiVJ7bNS1YV+0QAaKMFj9I/zc5WTkSxJ2llaqZveoH80AAAAACC2EUYDh2jmmpB+0d1byuxtAgsADdCmeaoeOq+mf/R3qwr02Ber3C0KAAAAAIBGRBgNHKIZIf2iR/agXzSAg3ds79a67oRezvjhz1Zo+uqC/ZwBAAAAAED0IowGDoHPbzVrXc3K6KMJowEcohtP6q2jugd6zvutdOPrc7WztNLlqgAAAAAACD/CaOAQLN1aot0VgX7RbTJT1K0V/aIBHJrEBI8mnz9ULZsF+kfn7w70j/bTPxoAAAAAEGMIo4FDMGf9LufxUfSLBnCY2mWl6u/nDnHG36zcqSe+Wu1iRQAAAAAAhB9hNHAIQsPoEV1buFgJgFgxtm8bTRrb0xk/+MlyzVpbuJ8zAAAAAACILoTRwCHI21ATRg8jjAYQJjef0kdHdgv8neK30vWv5amA/tEAAAAAgBhBGA0cpPySCm3aVS5JSk3yqH/75i5XBCBWJCZ4NPmCoWqRniRJ2l5SqZvfnE//aAAAAABATCCMBg5S6KrowZ2ylZTAtxGA8Gmflaa/n5vrjL9asUNPfb3GxYoAAAAAAAgPUjTgIOVtKHIeD+tCiw4A4XdCvza6+vgezviBT5Zr9jr6RwMAAAAAohthNHCQQjcvHE6/aACN5JZT+zp/x/j8Vte/Nle79lS5XBUAAAAAAIeOMBo4CJVenxZuLnbGQ7tku1gNgFiWFOwfnZUW6B+9tbhCv3qL/tEAAAAAgOgVtjDaGNPJGPO8MWaLMabSGLPOGPOwMabBS0eNMacYYx40xvzPGFNojLHGmG/38/qOxpjrjTEfBu9XaYwpMMZ8aoz5eT3njA1et76vew/lz4/4sHhLiaq8fklSt1bpyslIcbkiALGsY3aaHjxniDP+fFm+nv2W/tEAAAAAgOiUGI6LGGN6SpomqY2kf0taJukoSTdKGmeMGWOtLWjApa6VdKakCkmrJB0oyL5e0m8krZX0haRtkrpK+rmkk40xD1lrb67n3K8kfbmP+XrDbyAvpEUH/aIBNIWTB7TVlcd017PfrpUk3f/Rco3o1pK/gwAAAAAAUScsYbSkxxUIom+w1j6yd9IY83dJN0m6R9IvG3Cd+yTdrkCY3VmBkHl/Zkkaa639KnTSGNNf0gxJNxljplhr5+zj3C+ttXc2oCbAkbehJoweSr9oAE3k1+P6afb6XZq3sUhev9X1U+fq/RuOUXZ6stulAQAAAADQYIfdpsMY00PSqZLWSXqsztN3SNojaYIxptmBrmWtnW6tXWyt9TXk3tbaf9YNooPzSyW9ERyObci1gIbIW1/kPB7OqkQATSQ50aNHLhiq5qmB95A3F5XrlrcWyFr6RwMAAAAAokc4ekafGDx+Yq31hz5hrd0t6TtJ6ZKODsO9DkZ18Oit5/lexpjrjDG/M8Zcbozp3VSFITptKSrXtpIKSVKz5AT1bZfpckUA4knnlun6W0j/6M+Wbtdz3x7oA0QAAAAAAESOcLTp6Bs8rqjn+ZUKrJzuI+l/YbjfARljmks6W5KV9Ek9L7so+BV63juSrrLW7tr3KT+4z77af0hSvwaWiigyJ6RfdG6XbCV4jIvVAIhHPzqinS4b000vfLdOknTfR8s0oltL5XbOdrcwAAAAAAAaIBwro7OCx+J6nt873yT/UjbGGEnPSmor6Ylgy45QOyTdJmmQpExJrSWdJmmuAgH2f40x4fjfBTEmtF80G4cBcMtvT+uvIZ0CP3qrfVbXTc1TcXn1Ac4CAAAAAMB9TRG67l0+2lSNLR+UdI6kbyTdXPfJYE/q+6y1i6y1pdbandbajxToLb1W0hhJP2nIjay1w/f1pcAGjIgxeesJowG4LznRo0cvHKbMYP/oTbvK9eu359M/GgAAAAAQ8cIRRu9d+ZxVz/PN67yu0Rhj/ibpJklfSzrdWlvZ0HOttSWSpgaHxzVCeYhiFdU+Ld5S4oyHduEj8QDc07lluv42frAz/njxdr04bZ17BQEAAAAA0ADhCKOXB4996nl+78aA9fWUDgtjzEOSbpH0haTTrLWlh3CZHcFjs7AVhpiwYFOxvP7AqsOerZspOz3Z5YoAxLtxA9tr4qiuzvgvHyzVgk1FLlYEAAAAAMD+hSOM/iJ4PLVur2VjTKYCbS/KJc0Iw71+wAQ8Jun/JH0q6QxrbdkhXu7o4HFNWIpDzAjtFz28Ky06AESG353RXwM7Bj6AFOgfPVclFfSPBgAAAABEpsMOo621qyV9IqmbpGvrPH2XAquMX7bW7tk7aYzpZ4zpd7j3Dm5W+LSkSZI+lPRTa235Ac4Zs68NCo0xF0s6T1KVpDcPtzbEFvpFA4hEKYkJeuzCYcpMCfSP3lBYptveWUD/aAAAAABAREoM03UmSZomabIx5iRJSyWNlHSCAu05bq/z+qXBowmdNMYcI+nK4DAjeOxtjHlx72ustZeGnPLH4OvLJc2TdFsgn65lnrX23ZDxFEkeY8w0SZskpUo6UtJRkrySrrbWrjvQHxjxw1pba2X0MFZGA4ggXVs1071nD9a1U/MkSR8s3KaXpq3TpWO6u1wZAAAAAAC1hSWMttauNsaMkPQnSeMknS5pq6TJku6y1hY28FK9JE2sM9emztylIY/3/ks7TdJv67nmS5JCw+gnJJ2sQPuQHAUC8c2SXpT0sLV2fgNrRZzYWFiunaVVkqTM1ET1ap1xgDMAoGmdMbi9pq/poldnbJAk3fPBUg3pnK2hfJIDAAAAABBBwrUyWtbajZIua+Brf7B8OTj/ogKhcEPvealqh9MNOec+SfcdzDmIb3M21LyXMrRLC3k8+/zPFwBc9fszBmjexiIt2lzi9I9+7/pj1KIZG64CAAAAACJDODYwBGJa3voi5/FwVhkCiFCpSQl64qLhap4aeJ95c1G5bnpznvx++kcDAAAAACIDYTRwAHNCNy/smu1iJQCwf51bpuvBc3Od8ZfLd+iJr1a7WBEAAAAAADUIo4H92FPp1bJtJZIkY6TczoTRACLbKQPa6urjezjjBz9Zrmmrd7pYEQAAAAAAAYTRwH7M31SkvZ9w79s2U5mpSe4WBAANcOupfXVUt5aSJL+VbnhtnvJLKlyuCgAAAAAQ7wijgf2Yt7GmX/RQ+kUDiBKJCR49cuFQ5WQENi/cWVqp616bK6/P73JlAAAAAIB4RhgN7MeCjcXO49zOWS5WAgAHp23zVP3j/KHymMB41tpCPfDJCneLAgAAAADENcJoYD8Wbq4Jowd1pF80gOgypleObjq5jzN+8qvV+mzJdhcrAgAAAADEM8JooB47Syu1uahckpSS6FHvthkuVwQAB+/aE3ppbN/WzvjmN+dpY2GZixUBAAAAAOIVYTRQj9BV0Ud0aK6kBL5dAEQfj8fooXNz1SErVZJUUuHVtVPzVOn1uVwZAAAAACDekK4B9QjtFz24Ey06AESvFs2S9ehFw5SUEGggvWBTse5+b6nLVQEAAAAA4g1hNFCPhZuLnMeDOrJ5IYDoNqxLC/3u9P7O+JUZ6/XveZtdrAgAAAAAEG8Io4F6LNgUujKaMBpA9Lt0dDedMai9M/7tPxdqVf5uFysCAAAAAMQTwmhgH7aXVCh/d6UkKT05QT1as3khgOhnjNG9Zw9S95xmkqSyKp+ueTVPZVVelysDAAAAAMQDwmhgH+ZvrGnRMbBjlhI8xsVqACB8MlOT9PhFw5SSGPgVYGV+qW7/1yJZa12uDAAAAAAQ6wijgX1YuDmkRQf9ogHEmP7tm+vuswY643/N3azXZm10sSIAAAAAQDwgjAb2IbRf9CD6RQOIQeeM6KxzR3Ryxnf+Z7EWhbwRBwAAAABAuBFGA3VYa7VgU02bjsGdsl2sBgAaz5/OHKh+7TIlSVU+v66ZMkfFZdUuVwUAAAAAiFWE0UAdm3aVa1cwjMlMTVS3VukuVwQAjSM1KUFPXDxcGSmJkqSNheW65e359I8GAAAAADQKwmigjlr9ojtlyRg2LwQQu7rnNNPfxg92xp8u2a5nvlnjYkUAAAAAgFhFGA3UUatfdEdadACIfacNaq/Lx3R3xvd9tFyz1ha6WBEAAAAAIBYRRgN11O4XzeaFAOLDbaf109AugTfgfH6r61/L087SSperAgAAAADEEsJoIITfb3/QpgMA4kFyokePXThMLdKTJEnbSyp1/dS58vr8LlcGAAAAAIgVhNFAiPWFZdpd4ZUktWyWrI7ZaS5XBABNp0N2mh4+f6j2tsqfvqZAD3yywt2iAAAAAAAxgzAaCBHaomNQRzYvBBB/ju/TWjee1NsZP/nVan20aJuLFQEAAAAAYgVhNBAidPNCWnQAiFc3nNhbJ/Rt7YxveWu+1uwodbEiAAAAAEAsIIwGQiysFUZnu1gJALjH4zF66LxcdWoRaFVUWunVL1+do7Iqr8uVAQAAAACiGWE0EOTzWy3awspoAJCk7PRkPXnxcKUkBn5VWLG9VLe9s1DWWpcrAwAAAABEK8JoIGjNjlKVVfkkSW0yU9S2earLFQGAuwZ2zNKfzxrojP8zf4temrbOvYIAAAAAAFGNMBoImk+/aAD4gXNHdNYFR3Vxxne/v1Sz1xW6WBEAAAAAIFoRRgNBCzcVOY/pFw0ANe74yQDnTTqv3+raqXnK313hclUAAAAAgGhDGA0ELdhcszJ6ECujAcCRmpSgxy8aphbpSZKk7SWVun7qXHl9fpcrAwAAAABEE8JoQFK1z68lW0qc8aCOhNEAEKpTi3RNvmCojAmMZ64t1P0fL3e3KAAAAABAVCGMBiSt2L5bld7ACr+O2WnKyUhxuSIAiDzH9m6tX53Sxxk//fUafbBwq4sVAQAAAACiCWE0IGkhmxcCQINMGttLJ/dv44xvfWu+VuWXulgRAAAAACBaEEYDkhZtqQmjB9KiAwDq5fEYPXhurrq2Spck7any6ZevztGeSq/LlQEAAAAAIh1hNCDV6hd9RIfmLlYCAJEvKy1JT1w0XKlJgV8jVuWX6tfvLJC11uXKAAAAAACRjDAacc/vt1q2bbczHkAYDQAHNKBDc/3lZ4Oc8fsLtuq5b9e6WBEAAAAAINIRRiPurS8sU1mVT5KUk5GsNpmpLlcEANHh58M6acLRXZ3xXz9cpmmrd7pYEQAAAAAgkhFGI+6Ftujo355V0QBwMP7w4wEa2iVbkuTzW10/da42F5W7XBUAAAAAIBIRRiPuLd1aE0bTogMADk5yokdPXDRcORkpkqSCPVX65StzVFHtc7kyAAAAAECkIYxG3FsSGkazMhoADlq7rFQ9cfEwJXqMJGnh5mL9/t1FbGgIAAAAAKiFMBpxL7RNB2E0AByaI7u11B9/MsAZvz1nk16dsd7FigAAAAAAkYYwGnGtcE+VtpVUSJJSEj3qntPM5YoAIHpNOLqrxg/v5Izv+u8Sfb+u0MWKAAAAAACRhDAacS20X3TfdplKTOBbAgAOlTFGd581UIM7ZUmSvH6ra17N07biCpcrAwAAAABEApI3xDVadABAeKUmJeiJi4erZbNkSdLO0kpdM2WOKr1saAgAAAAA8Y4wGnEtdGX0gA6E0QAQDh2z0/TohUOVENzQcO6GIt35nyUuVwUAAAAAcBthNOLakpAwuj8rowEgbEb3zNHvTu/vjF+btUGvzdrgYkUAAAAAALcRRiNuVXp9WpVf6oz7tct0sRoAiD2Xj+mmM3M7OOM7/r1YeRt2uVgRAAAAAMBNhNGIWyu3l8rrt5Kkrq3SlZma5HJFABBbjDG69+eDnU+eVPn8uubVOcrfzYaGAAAAABCPCKMRt2q16GhHiw4AaAxpyQl6esJwZacH3vDbXlKpa6fkqcrrd7kyAAAAAEBTI4xG3Fqyhc0LAaApdG6ZrkcuGKrgfob6ft0u3fM+GxoCAAAAQLwhjEbcWhqyMnoAmxcCQKM6tndr/XpcP2f80vT1emv2RhcrAgAAAAA0NcJoxCVrbe02HayMBoBGd/VxPXTGoPbO+PZ3F2nexiIXKwIAAAAANCXCaMSlTbvKtbvCK0nKSktSh6xUlysCgNhnjNH94werb9tMSVKV16+rX5mt/BI2NAQAAACAeEAYjbhUt0WHMcbFagAgfjRLSdTTlwxXVlrNhoZXvzpHFdU+lysDAAAAADQ2wmjEpVotOugXDQBNqmurZnrswmHOhoZzNxTp9+8ukrXW3cIAAAAAAI2KMBpxqdbKaPpFA0CTO6Z3jm4/Y4AzfnvOJr3w3Tr3CgIAAAAANDrCaMSlJXXadAAAmt7lY7pp/PBOzvieD5bq25U7XawIAAAAANCYCKMRd0oqqrWxsFySlJRg1KtNhssVAUB8Msbo7rMGKrdztiTJ57e6dmqe1hfscbkyAAAAAEBjIIxG3Fm2dbfzuFebTCUn8m0AAG5JTUrQ0xOGq23zFElScXm1rnp5tkorvS5XBgAAAAAIN1I4xJ0lW4qdx7ToAAD3tWmeqqcmjHDeHFyxvVQ3vTFPfj8bGgIAAABALCGMRtwJ7Rfdv32mi5UAAPbK7Zytv/5skDP+dMl2PfzZChcrAgAAAACEG2E04s7SkDYdAzqwMhoAIsXZwzvpymO6O+PJn6/SBwu3ulgRAAAAACCcCKMRV6p9fi3fHhJG06YDACLKbaf107G9c5zxr96cryVbSvZzBgAAAAAgWhBGI66s2bFHVV6/JKlDVqqy05NdrggAECoxwaNHLximbq3SJUnl1T5d9fJsFZRWulwZAAAAAOBwhS2MNsZ0MsY8b4zZYoypNMasM8Y8bIxpcRDXOMUY86Ax5n/GmEJjjDXGfNuA8wYYY940xuQbYyqMMcuNMXcZY9L2c85oY8wHwfuUGWMWGGP+zxiT0NB6EX2WhvSLpkUHAESmrPQkPXPJCGWkJEqSNheVa9KUPFX7/C5XBgAAAAA4HGEJo40xPSXNkXSZpFmSHpK0RtKNkqYbY1o18FLXSrpZ0mhJmxt475GSvpd0lqTPJP1DUomkP0r61BiTso9zzpT0taTjJP1L0mOSkoN1v97AWhGFam9eSBgNAJGqd9tMPXxerowJjGeuLdSf/rvE3aIAAAAAAIclXCujH5fURtIN1tqzrLW3WWtPVCDc7SvpngZe5z5JAyVlSPrJgV4cXMX8gqR0SeOttRdaa38jaaSkdySNkXRTnXOaS3pGkk/SWGvtFdbaWyXlSpouabwx5vwG1osos5QwGgCixskD2uqWU/s641dmrNerM9a7WBEAAAAA4HAcdhhtjOkh6VRJ6xRYYRzqDkl7JE0wxjQ70LWstdOttYuttb4G3v54Sf0lfW2t/U/IdfySfh0c/tKYveuqJEnjJbWW9Lq1dnbIORWSfh8cXtPA+yPKrNxe6jzu2y7TxUoAAA0xaWxPnTGovTO+4z+LNW3VThcrAgAAAAAcqnCsjD4xePwkGAI7rLW7JX2nwMrlo8Nwr/ru/VHdJ6y1ayStkNRVUo+GnKNA644ySaP31d4D0a24vFrbSiokSckJHnVtme5yRQCAAzHG6IFzhmhgx8CnWXx+q2um5GnNjtIDnAkAAAAAiDThCKP3fn52RT3Prwwe+4ThXuG4d73nWGu9ktZKSlTtABsxYOX23c7jnm0ylJgQtv07AQCNKC05Qc9cMkJtMgPvExeXV+vKl2aruKza5coAAAAAAAcjHGlcVvBYXM/ze+ezw3CvcNw7bPUaY+bs60tSvwOdi6a3PCSM7ts2w8VKAAAHq31Wmp65ZIRSEgO/uqzZuUeTps5Rtc9/gDMBAAAAAJGiKZaG7u3XbJvgXuG4t5v1ohGt2FYTRvduS79oAIg2Qzpn68Fzhzjj71YV6E//XeJiRQAAAACAg5EYhmvsXUmcVc/zzeu8LpwO5d5hq9daO3xf88HV0cMOdD6a1orQzQsJowEgKv14cAetyi/Vw58FOnG9MmO9erfN0CWjurlbGAAAAADggMKxMnp58FhfT+jewWN9fZ2b+t71nmOMSZTUjI2jAAAAIABJREFUXZJX0ppwFIjIsSKkTUcfwmgAiFo3ntRbPx7c3hnf9d8l+nrFDhcrAgAAAAA0RDjC6C+Cx1ONMbWuZ4zJlDRGUrmkGWG4V12fB4/j6j5hjOmhQOC8XrWD5XrPkXScpHRJ06y1lWGsEy7bWVqpgj1VkqS0pAR1apHmckUAgENljNED5wzRkE6BDzn5/FbXTs3TqvzSA5wJAAAAAHDTYYfR1trVkj6R1E3StXWevktSM0kvW2v37J00xvQzxoRjk7+vJC2VdJwx5qch1/dIui84fNJaG9r/+W1JOyWdb4wZEXJOqqS7g8MnwlAbIkjtVdEZ8njMfl4NAIh0qUkJevqSEWrXPFWStLvCqyte+l67gm88AgAAAAAiTzh6RkvSJEnTJE02xpykQEA8UtIJCrTIuL3O65cGj7USQWPMMZKuDA4zgsfexpgX977GWntpyGOfMeYyBVY7v22MeVvSBkknSRoh6TtJD4Xew1pbYoy5SoFQ+ktjzOuSCiX9VFLf4PwbB/fHR6Rj80IAiD1tm6fq2YkjNP7Jaaqo9mt9QZmumTJHL18+UsmJTbFHMwAAAADgYITlX2rB1dEjJL2oQAj9K0k9JU2WNMpaW9DAS/WSNDH4dXZwrk3I3MR93HumpCMl/VvSqZJuUmBzwj9JOmVf7Taste9KOl7S18H7XC+pWtLNks6vs5IaMWBFPpsXAkAsGtgxSw+dm+uMZ6wp1B3/WSx+lAMAAABA5AnXymhZazdKuqyBr91njwRr7YsKBNoHe+8lks45yHO+k3T6wd4L0an2yuiM/bwSABBtThvUXrec2kcPfBLYr/i1WRvUp22GLhvT3eXKAAAAAACh+AwrYp61tlbP6L7tWBkNALHm2hN66czcDs74z+8t0ZfL812sCAAAAABQF2E0Yt72kkqVVHglSZmpic5mVwCA2GGM0X1nD1Zu52xJkt9K10+dq5Uhb0YCAAAAANxFGI2YtzwkiOjTNlPG7LNLDAAgyqUmJejpS4arQ1bgTcfdlV5d9uL32rH7B9tHAAAAAABcQBiNmLeyThgNAIhdbTJT9ezEI5WenCBJ2rSrXFe+PFvlVT6XKwMAAAAAEEYj5i3fFhpGs3khAMS6AR2a65ELhsoT/CDM/I1FuvnNefL7rbuFAQAAAECcI4xGzFuRX+o87svKaACICyf1b6s7fnKEM/5w0Tbd9/EyFysCAAAAABBGI6b5/bZ2m452hNEAEC8mju6my8Z0c8ZPfbVGU2ducK8gAAAAAIhzhNGIaZuLylUW7BPaslmycjJSXK4IANCUfn/GAJ3cv40z/sO/F+nrFTtcrAgAAAAA4hdhNGLaiu30iwaAeJbgMfrH+UM1sGNzSZLPbzVpSl6t/QQAAAAAAE2DMBoxbXmtMJoWHQAQj5qlJOq5iUeqQ1aqJKm00qvLX/xe+SUVLlcGAAAAAPGFMBoxbeX2ms0LCaMBIH61bZ6q5y49UhkpiZICbZyufHm2yqq8LlcGAAAAAPGDMBoxLfRj2H3ZvBAA4lr/9s312EXDlOAxkqQFm4p14+vz5PNblysDAAAAgPhAGI2Y5fNbrdoRsjK6DWE0AMS74/u01p/OPMIZf7pku/7ywVIXKwIAAACA+EEYjZi1vmCPqrx+SVLb5inKSk9yuSIAQCS4aGRX/eK4Hs74uW/X6uXp61yrBwAAAADiBWE0YtYKNi8EANTjtnH9NO6Ids74zv8s1hfL8l2sCAAAAABiH2E0YtYKNi8EANTD4zF66LxcDemUJUnyW+m6qXlavKXY5coAAAAAIHYRRiNmLQ9ZGd2XMBoAUEdacoKenXikOmanSZL2VPl0+Yvfa3NRucuVAQAAAEBsIoxGzFoZEkb3bpvhYiUAgEjVOjNFL152pDJTEyVJ20sqdenzs1RcVu1yZQAAAAAQewijEZOqvH6t2bHHGfdmZTQAoB6922bqqQnDlZRgJEkr80t11SuzVVHtc7kyAAAAAIgthNGISWt37pHXbyVJHbPTlJGS6HJFAIBINrpnjh44Z4gznrW2UL96a778wZ8lAAAAAIDDRxiNmLQitF90O1ZFAwAO7Mzcjvrtaf2c8fsLtuqeD5a6WBEAAAAAxBbCaMSk0DC6Dy06AAAN9IvjeujS0d2c8XPfrtWz36xxryAAAAAAiCGE0YhJtcNoNi8EADSMMUZ/+PEAjTuinTN39/tL9d6CLS5WBQAAAACxgTAaMWllfqnzmJXRAICDkeAxevj8XI3o2sKZu/mN+Zq5psDFqgAAAAAg+hFGI+ZU+/zaUFDmjHu0buZiNQCAaJSalKBnLhnh/Ayp8vl11cuza33yBgAAAABwcAijEXPWF5TJ67eSpI7ZaUpPTnS5IgBANGrRLFkvXXaUWmemSJJKKry69PlZ2lZc4XJlAAAAABCdCKMRc1bvqGnRwapoAMDh6NwyXS9ceqSaJSdIkrYUV+jSF2appKLa5coAAAAAIPoQRiPmhIbRPVuzeSEA4PAM7JilJy4erkSPkSQt27Zb17w6R1Vev8uVAQAAAEB0IYxGzFmdv8d53LMNYTQA4PAd16e17j17sDP+blWBfv32fPmDbaEAAAAAAAdGGI2YU3tlNG06AADhMX54J91yah9n/O68Lbr/4+UuVgQAAAAA0YUwGjHFWlsrjO5Fmw4AQBhde0IvXXBUF2f85Fer9ew3a1ysCAAAAACiB2E0YsqO0krtrvBKkjJTEtU6M8XligAAscQYoz+feYRO7t/Gmbv7/aX6Z94mF6sCAAAAgOhAGI2YEtovukebDBljXKwGABCLEhM8euSCYRrRtYUzd+vbC/T5su0uVgUAAAAAkY8wGjGFftEAgKaQlpyg5yYeqX7tMiVJPr/VpCl5mrO+0OXKAAAAACByEUYjptQOo+kXDQBoPFnpSXrp8qPUqUWaJKmi2q/LXvhey7ftdrkyAAAAAIhMhNGIKat31LTpIIwGADS2ts1T9coVI5WTkSxJKqnw6pLnZ2pjYZnLlQEAAABA5CGMRkxZnV+zMrpXG9p0AAAaX/ecZnrxsqOUkZIoSdpeUqlLnp+lnaWVLlcGAAAAAJGFMBoxo7zKp81F5ZKkBI9Rl5aE0QCApjGwY5aevmS4khMCv1qt3blHl73wvUorvS5XBgAAAACRgzAaMWPNzppV0V1bpis5kf+8AQBNZ3TPHE2+IFceExgv3FysX7w8W5Ven7uFAQAAAECEIK1DzAjtF92DftEAABeMG9he9/xskDOetrpAN70xTz6/dbEqAAAAAIgMhNGIGaH9onvSLxoA4JILjuqiW3/U1xl/sHCb/vDvRbKWQBoAAABAfCOMRsxYvSMkjGZlNADARZPG9tTlY7o746kzN+ihT1e4WBEAAAAAuI8wGjEjtE0HYTQAwE3GGP3+jP46K7eDMzf581V68bu1LlYFAAAAAO4ijEZM8Put1tRaGU2bDgCAuzweo7+dM0Rj+7Z25u787xL9M2+Ti1UBAAAAgHsIoxETNheVq9LrlyTlZCQrOz3Z5YoAAJCSEjx6/KJhGtYl25m79e0F+mjRNherAgAAAAB3EEYjJoT2i+5Biw4AQARJT07U85ceqX7tMiVJPr/V9a/l6asVO1yuDAAAAACaFmE0YgL9ogEAkSw7PVmvXDFS3XMCbaSqfVZXvzJbs9YWulwZAAAAADQdwmjEhNX0iwYARLjWmSl69cqR6pidJkmqqPbr8he/14JNRS5XBgAAAABNgzAaMWF1fkgY3YaV0QCAyNQxO02vXjlSORkpkqTSSq8mPj9LK7bvdrkyAAAAAGh8hNGICaFtOnrRpgMAEMG65zTTlCtHKjs9SZK0q6xaFz07U+t27jnAmQAAAAAQ3QijEfWKy6q1s7RSkpSS6FGH4MefAQCIVH3bZeqly45SRkqiJGnH7kpd9OxMbSkqd7kyAAAAAGg8hNGIeqt31rTo6J7TTAke42I1AAA0zJDO2Xpu4gilJAZ+HdtcVK6Ln53pvMEKAAAAALGGMBpRj37RAIBoNbJHKz01YbiSEgJvpK7ZuUcTnpul4rJqlysDAAAAgPAjjEbUC+0X3ZN+0QCAKDO2bxtNPn+o9n6wZ+nWEk18YZZKK73uFgYAAAAAYUYYjai3ekfIyujWzVysBACAQ3PaoPa6f/wQZzxvY5Guemm2Kqp9LlYFAAAAAOFFGI2oVzuMZmU0ACA6jR/eSX868whnPH1NgSZNyVOV1+9iVQAAAAAQPoTRiGrVPr82FJQ54x6sjAYARLFLRnXTr8f1dcafL8vXja/PVbWPQBoAAABA9COMRlRbX1Amr99Kkjpmpyk9OdHligAAODyTxvbSpLE9nfGHi7bppjfmyUsgDQAAACDKEUYjqoW26GBVNAAgVtz6o766fEx3Z/zegq265a358gXfgAUAAACAaEQYjahGv2gAQCwyxugPP+6vS0Z1debenbdFv3lngfwE0gAAAACiFGE0otrq/D3O455tCKMBALHDGKM7f3KELhzZxZl7e84m/e5fCwmkAQAAAEQlwmhEtdoro2nTAQCILR6P0d1nDtR5Izo7c69/v1F/+PciWUsgDQAAACC6EEYjallrtYY2HQCAGOfxGP3154P082EdnbkpMzfozv8sJpAGAAAAEFUIoxG1dpVVq6TCK0lKT05Qm8wUlysCAKBxeDxGfxs/RGfmdnDmXpq+Xne/v5RAGgAAAEDUIIxG1FpXUNMvumurZjLGuFgNAACNK8Fj9OA5Q3TG4PbO3HPfrtW9Hy0jkAYAAAAQFQijEbXW7awJo7vnpLtYCQAATSMxwaOHz8vVaQPbOXNPfbVGD36ygkAaAAAAQMQLWxhtjOlkjHneGLPFGFNpjFlnjHnYGNPiIK/TMnjeuuB1tgSv22kfr73UGGMP8OWrc063A7z+9cP93wJNIzSM7taKzQsBAPEhKcGjf5w/VCf3b+vMPfrFKv3jfytdrAoAAAAADiwxHBcxxvSUNE1SG0n/lrRM0lGSbpQ0zhgzxlpb0IDrtApep4+kzyW9LqmfpMsknWGMGWWtXRNyyjxJd9VzuWMlnSjpw3qeny/p3X3MLzpQnYgMawvKnMeE0QCAeJKc6NFjFw3VNa/m6fNl+ZKkhz9bqUSP0XUn9na5OgAAAADYt7CE0ZIeVyCIvsFa+8jeSWPM3yXdJOkeSb9swHX+okAQ/ZC19uaQ69wg6R/B+4zbO2+tnadAIP0DxpjpwYdP13OvedbaOxtQEyLU+pCe0d1yCKMBAPElJTFBj180TL94ZY6+XrFDkvTAJyvk8RhNGtvL5eoAAAAA4IcOu02HMaaHpFMlrZP0WJ2n75C0R9IEY8x+08Lg8xOCr7+jztOPBq//o+D9DlTTQElHS9os6f0D/iEQday1WhvapoOe0QCAOJSalKCnJwzXMb1ynLn7P1quR2jZAQAAACAChaNn9InB4yfWWn/oE9ba3ZK+k5SuQDi8P6MkpUn6Lnhe6HX8kj4JDk9oQE1XB4/PWWt99bymgzHmamPM74LHwQ24LiJE4Z4q7a7wSpKaJSeodUaKyxUBAOCO1KQEPXPJCI3q0cqZe/DTFXroUzY1BAAAABBZwtGmo2/wuKKe51cqsHK6j6T/HeZ1FLxOvYwxaZIuluSX9Ox+XnpK8Cv03C8lTbTWbtjfPUJeP6eep/o15HwcunUhLTq6tmomY4yL1QAA4K605AQ9f+mRuurl2fp21U5J0j/+t1Jev1+3nNqXn5MAAAAAIkI4VkZnBY/F9Ty/dz67ia5zbvA1H1prN+7j+TJJf5Y0XFKL4Nfxkr6QNFbS/w7UUgTuW7uzZvPC7vSLBgBAackJenbiCB3fp7Uz99gXq3Xvh8tYIQ0AAAAgIoQjjD6QvUtxDvdfQQ29zi+Cx6f29aS1Nt9a+0drbZ61tij49bUCq7dnSuol6cqGFGStHb6vL0nLGnI+Dl3tzQvpFw0AgBRo2fHUhOE6sV8bZ+6pr9fo7veXEkgDAAAAcF04wui9K5az6nm+eZ3XNdp1jDEDJI2WtEnSBwe4Xy3WWq9q2nocdzDnounV2rywFSujAQDYKzUpQU9ePFynDGjrzD337Vrd9d8lBNIAAAAAXBWOMHp58FhfL+fewWN9vaDDeZ2GbFy4PzuCR9LNCBfaM5o2HQAA1Jac6NHjFw3TaQPbOXMvTlun37+7SH4/gTQAAAAAd4QjjP4ieDzVGFPresaYTEljJJVLmnGA68wIvm5M8LzQ63gUaKMRej/VeU2qpAkKbFz43MH8AUIcHTyuOcTz0QSstVoX0jO6KyujAQD4gaQEjyZfMFQ/HtzemZsyc4N+96+FBNIAAAAAXHHYYbS1drWkTyR1k3RtnafvUmCV8cvWWmcpqzGmnzGmX53rlEp6Jfj6O+tc57rg9T+21tYXFJ+jwGaEH9SzceHee480xiTvY/5ESTcFh6/Wdz7cV7CnSqWVXklSRkqicjJ+8H8nAABQIJB++LxcnZXbwZl7/fuNuvXtBfIRSAMAAABoYolhus4kSdMkTTbGnCRpqaSRkk5QoK3G7XVevzR4NHXmfydprKSbjTG5kmZJ6i/pTEn5+mHYHWrvxoVPH6DW+yQdYYz5UoHe0pI0WNKJwcd/sNZOO8A14KJ1O2tvXmhM3f+MAADAXokJHj14bq4SPB69kxf41eedvE3y+f164JwhSkxoiv2sAQAAACBMYbS1drUxZoSkP0kaJ+l0SVslTZZ0l7W2sIHXKTDGjJJ0h6SzJB0rqUDSC5L+aK3dtK/zjDH9JR2jhm1c+Iqkn0k6UtJpkpIkbZf0pqRHrbXfNKRWuIfNCwEAODgJHqO/jR+sRI/RG7MDHyB7d94Wef1WD5+XSyANAAAAoEmEa2W0gq0xLmvga+tdyhoMrm8MfjX03kv1w1XW9b32OR16T2lEgNDNCwmjAQBoGI/H6K8/H6TEBKMpMzdIkt5bsFVen9XkC4YqOZFAGgAAAEDj4l8diDqhmxd2yyGMBgCgoTweo7vPGqiJo7o6cx8t3qarXp6t8iqfi5UBAAAAiAeE0Yg6oSuju+eku1gJAADRxxijO396hK48prsz99WKHZr4/Cztrqh2sTIAAAAAsY4wGlHFWlt7A0PadAAAcNCMMbr9jP668aTeztysdYW66NmZ2rWnysXKAAAAAMQywmhElR2lldoT/BhxZkqiWjZLdrkiAACikzFGN53SR7ef3t+ZW7CpWOc+NV3bSypcrAwAAABArCKMRlSp2y/amAbtWwkAAOpx1XE99JefDdLeH6kr80t1zpPTtbGwbP8nAgAAAMBBIoxGVAntF83mhQAAhMeFI7vo4fNylegJJNIbCst0zpPTtSq/1OXKAAAAAMQSwmhEldB+0d1bsXkhAADhcmZuRz158XAlJwZ+PdxWUqFzn5quRZuLXa4MAAAAQKwgjEZUCV0Z3ZXNCwEACKuTB7TVi5ceqfTkBElS4Z4qXfDMDM1eV+hyZQAAAABiAWE0osraOj2jAQBAeI3ulaNXrxyp5qmJkqTdFV5NeG6Wvlm5w+XKAAAAAEQ7wmhEDWut1oesjO5OGA0AQKMY1qWF3rh6lHIykiVJ5dU+XfHibH28eJvLlQEAAACIZoTRiBo7dleqrMonSWqemqgW6UkuVwQAQOzq37653rx6lDpkpUqSqnx+TZqSp3/N3eRyZQAAAACiFWE0osbakM0Lu+U0kzHGxWoAAIh9PVpn6K1rRqtbcNNgn9/qpjfm6/lv17pcGQAAAIBoRBiNqBG6eWE3Ni8EAKBJdMxO05u/HKV+7TKduT+9t0T3f7RM1loXKwMAAAAQbQijETXYvBAAAHe0yUzV6784WsO7tnDmHv9ytX7zzgJ5fX4XKwMAAAAQTQijETVqb16Y7mIlAADEn+z0ZL16xUid1K+NM/fm7E365atzVB7c0wEAAAAA9ocwGlEjtGd0V9p0AADQ5NKSE/TUhOE6Z3gnZ+6zpfma8NxMFZdVu1gZAAAAgGhAGI2oYK3V+oKaNh3dCaMBAHBFYoJH948frGvG9nTmZq/fpXOemqatxeUuVgYAAAAg0hFGIypsL6lUeXXgI8BZaUlq0SzZ5YoAAIhfxhj9Zlw//eHHA5y5FdtLdfbj07Qqf7eLlQEAAACIZITRiArrQvpFs3khAACR4Ypjuusf5+cq0WMkSVuKKzT+yemau2GXy5UBAAAAiESE0YgK60L6RXdrxeaFAABEijNzO+r5S49UenKCJKmorFoXPjNTXyzPd7kyAAAAAJGGMBpRYW3oymj6RQMAEFGO69NaU686Wi2DbbTKq3266qXZ+mfeJpcrAwAAABBJCKMRFUJXRnenTQcAABEnt3O23vrlKHXMTpMkef1WN785X898vcblygAAAABECsJoRIV1O8ucx/SMBgAgMvVsnaF/Thqtfu0ynbl7Pliqu/67WD6/dbEyAAAAAJGAMBoRz1qr9YX0jAYAIBq0bZ6qN64epaO6tXTmXvhunSZNmaPyKp+LlQEAAABwG2E0It6O3ZWqqPZLkrLSkpSdnuxyRQAAYH+y0pL08hVHadwR7Zy5jxdv1wXPzNDO0koXKwMAAADgJsJoRLwNhTUtOrq0ZFU0AADRIDUpQY9fNExXHtPdmZu3sUg/f3yaVu8odbEyAAAAAG4hjEbEI4wGACA6eTxGv//xAN310yPkMYG5DYVlOvuJaZq1ttDd4gAAAAA0OcJoRLzQMLozYTQAAFFn4uhuemrCCKUmBX71LCqr1sXPztR/529xuTIAAAAATYkwGhGPldEAAES/Uwa01Ru/GKWcjMDeD1U+v65/ba6e+HK1rLUuVwcAAACgKRBGI+JtJIwGACAm/D979x0mV1n3f/xzT9vZXrMl2RTSO5AQAimQhBaaiog8oqj8RJGuqNh7V0SkWR4UH7Aj0hUQkhASQiAhJIQ00rO9952dnZnz+2NmZ2fDJtmQ2T2zu+/Xdc11Zu5zzj33KOLkk+9875NHZ+mxGxZqYn5adOynz+7Q1x/fqkAwZOPKAAAAAAwEwmgkPCqjAQAYOkbnpOjRzy7Q/JNyomN/WX9Q1z60QS0dARtXBgAAAKC/EUYjofk6g6ps6pAkOR1GRVlem1cEAABOVGaKWw996nR94JSR0bFVO6t15W/XqbLJZ+PKAAAAAPQnwmgktJL69ujzokyv3E7+kQUAYChIcjn1yytP0c3LJkbH3i5r0mX3rdWOiiYbVwYAAACgv5DsIaHRLxoAgKHLGKMvnD9FP718lpwOI0kqa/Tp8vtf0YvbK21eHQAAAIB4I4xGQqNfNAAAQ9+V88bowU/OU1qSS5LU6g/q2oc26Her98iyLJtXBwAAACBeCKOR0GLD6NGE0QAADFlnTR6hR69foOLsZEmSZUk/+vcO3f7PLfIHQjavDgAAAEA8EEYjoVEZDQDA8DGlMF2P37hQp43Njo49srFEH3tgvWpbOmxcGQAAAIB4IIxGQqNnNAAAw0teWpL+/On5+tDc4ujYa/vr9P771mpnRbONKwMAAABwogijkbAsy6IyGgCAYSjJ5dTPPzRbX71wqkx4X0OV1Lfr8l+/ohU72NgQAAAAGKwIo5Gwalv9avMHJUnpSS5lpbhtXhEAABgoxhhdd/YE/e7q05TqcUqSWjoC+tT/bdADL+9lY0MAAABgECKMRsKKbdFRnJMi01UaBQAAho3zphfon9cv0Kis7o0Nf/DMdn3l0bfY2BAAAAAYZAijkbB6tuhItnElAADATtOKMvTETQs1N2Zjw79vOKSP/X696lr9Nq4MAAAAwPEgjEbCYvNCAADQJS8tSX++dr4+eOqo6Nhr++r0gfvW6p1KNjYEAAAABgPCaCQsNi8EAACxvG6nfvHhk/Xl5d0bGx6sa9MH7lur596usHdxAAAAAI6JMBoJKzaMHk0YDQAAFN7Y8PolE/Sbj81VSmRjw1Z/UNc9vFF3Pr9ToRAbGwIAAACJijAaCetQXXv0OZXRAAAg1gUzCvXo9Qs0OmZfibtX7Na1D21QY3unjSsDAAAAcCSE0UhI/kBIZY3hMNoYaVQ2GxgCAICephVl6MkbF2nxpLzo2IodVfSRBgAAABIUYTQSUmlDu6zIr2yLMrxKcjntXRAAAEhI2akePfjJebrurPHRsX01rfSRBgAAABIQYTQSEv2iAQBAX7mcDn31omm65yOnKtnds4/0L+gjDQAAACQMwmgkpEOE0QAA4DhdevJI/euGnn2k71mxW5/6v9fpIw0AAAAkAMJoJKTYMJrNCwEAQF9NK8rQUzf17CO9cmc1faQBAACABEAYjYR0kDAaAAC8R1kpHv3xmtN13dnv7iP97Fb6SAMAAAB2IYxGQqJnNAAAOBFOh9FXL5yme6/q2Uf6s3/aqDue26kgfaQBAACAAUcYjYRjWZYO1lIZDQAATtwls8N9pGO/T9y7crc+/of1qmnpsHFlAAAAwPBDGI2E09jeqeaOgCQp2e1UXprH5hUBAIDBbFpRhp68aaHOmjwiOrZ2d60uvvtlvb6/zsaVAQAAAMMLYTQSzuH9oo0xNq4GAAAMBVkpHj34yXm6ZdlEdX21qGzq0P/87lX97+q9sizadgAAAAD9jTAaCYd+0QAAoD84HUa3nT9FD35ynrJT3JKkYMjSD/+9Xdc9vFGN7Z02rxAAAAAY2gijkXAO1bVHn4/OSbZxJQAAYChaMiVfz9yyWKeOyYqOPb+tUu+7d43eLmu0cWUAAADA0EYYjYRzeJsOAACAeBuZlay/f+ZMXbNwXHTsQG2bLrv/Ff3ttYO07QAAAAD6AWE0Es4hwmgAADAAPC6Hvn3pDN131RylJbkkSf5ASF/511v64iNb1O4P2rxCAAAAYGghjEbCoTIaAAAMpItnF+nJmxZqamF6dOzRN0r0gfvWak91i40rAwAAAIYWwmgklEAwpNKG7p7RxdmE0QAAoP+NH5Gmx24mTGCVAAAgAElEQVRYqMvnFEfHdlY26333rNHTW8psXBkAAAAwdBBGI6GUN/oUDIV7NOanJynZ47R5RQAAYLhI9jh1xxWz9dPLZ8njCn9NbvUHddNfNunrj70lXydtOwAAAIATQRiNhEKLDgAAYCdjjK6cN0aP3bBAY3O7v4v8ef1Bvf/etdpV2Wzj6gAAAIDBjTAaCYUwGgAAJIIZIzP11M2LdNGswujYzspmve/eNfrL+oOyLMvG1QEAAACDU9zCaGNMsTHmD8aYMmNMhzFmvzHmLmNM9nHOkxO5b39knrLIvMVHuH6/McY6wqPiKO+zwBjzb2NMnTGmzRizxRjzOWMMfSFsdCgmjC4mjAYAADbK8Lp131Vz9KPLZikp0rbD1xnS1x57Szf9ZZMa2zttXiEAAAAwuLjiMYkxZoKkVyTlS3pC0g5Jp0u6VdJyY8xCy7Jq+zBPbmSeyZJWSPqbpKmSrpF0sTHmTMuy9vZya6Oku3oZ73X7c2PM+yU9Kskn6e+S6iRdKumXkhZKuuJYa0X/oDIaAAAkEmOMrpo/RqeNy9ZNf3lDuyrDXy+featcbx5q0N0fOVVzxx5X7QUAAAAwbMUljJZ0v8JB9C2WZd3TNWiMuVPS5yX9UNJn+zDPjxQOon9pWdZtMfPcIulXkfdZ3st9DZZlfacvCzXGZEj6X0lBSUssy9oQGf+mwgH4h4wx/2NZ1t/6Mh/i6xBhNAAASECTC9L15E2L9P2nt+nP6w9Kkkob2vXh367TbedN1vVnT5DDYWxeJQAAAJDYTrhNhzFmvKTzJe2XdN9hp78tqVXS1caY1GPMkyrp6sj13z7s9L2R+S+IvN+J+JCkEZL+1hVES5JlWT5J34i8vP4E3wPvEZXRAAAgUXndTv3wsln69UfnKMMbrukIhiz9/LmduvoP61XV5LN5hQAAAEBii0fP6GWR4/OWZYViT1iW1SxpraQUSWccY54zJSVLWhu5L3aekKTnIy+X9nJvkjHmY8aYrxljbjXGLD1K7+eu9T7by7nVktokLTDGJB1jvYizZl+n6tvCvRc9Lofy0/mvAAAAJJ4LZxXp37cu7tGeY+3uWl34q5e1cmeVjSsDAAAAEls8wugpkeOuI5x/J3Kc3I/zFEp6WOF2IHcp3G7jHWPM2cfzPpZlBSTtU7h9yYlWYOM4Haprjz4vzk7mp64AACBhFWen6O+fOUM3LZ0oE/nKUtvq1zUPvq4fPL1N/kDo6BMAAAAAw1A8wujMyLHxCOe7xrP6aZ4HJZ2jcCCdKmmWpN9KGifpP8aYk/tpvTLGbOztofCmizhOpQ3dYfSorGQbVwIAAHBsLqdDX7xgiv78qfk9ftH1wJp9+sB9a7WrsvkodwMAAADDTzzC6GPpKm+1+mMey7K+a1nWCsuyKi3LarMsa6tlWZ+VdKfCbT++E4/3Qf8rre/uF12cTb9oAAAwOCyYmKf/3LpYy6bmR8e2lTfpknvW6Pdr9ikU4mslAAAAIMUnjO6qJM48wvmMw67r73m6/CZyPKu/3seyrLm9PSTt6OMaESO2Mro4m8poAAAweOSmJen3nzhN3750ujyu8FdsfyCk7z+9TVf/Yb3KG9uPMQMAAAAw9MUjjN4ZOR6pJ/SkyPFIvaDjPU+Xrt1jUvv6PsYYl6STJAUk7e3j+yBOSupp0wEAAAYvY4yuWXiSnrl5kaYXZUTH1+6u1QW/XK0nN5fZuDoAAADAfvEIo1dGjucbY3rMZ4xJl7RQUrukV48xz6uR6xZG7oudxyHp/MPe71jOjBwPD5VXRI7Le7nnLEkpkl6xLKujj++DOOnRM5rKaAAAMEhNKkjX4zcu1A1LJkQ3N2zyBXTLXzfp1r9tUmNbp70LBAAAAGxywmG0ZVl7JD2v8IaBNx52+rsKVyY/ZFlWa9egMWaqMabHJn+WZbVIejhy/XcOm+emyPzPWZYVDZeNMTOMMTmHr8kYM1bSvZGXfzrs9D8l1Uj6H2PMaTH3eCX9IPLy171/WvSn0nradAAAgKHB43Lo9uVT9Y/rzuzxveaJN8u0/Fer9cruGhtXBwAAANjDFad5bpD0iqS7jTHnSNouab6kpQq31fj6YddvjxzNYeNfk7RE0m3GmFMkvSZpmqT3K9x24/Cw+wpJXzHGrJS0T1KzpAmSLpbklfRvSXfE3mBZVpMx5tMKh9KrjDF/k1Qn6X2SpkTG/358Hx8nqt0fVG2rX5Lkchjlp3ttXhEAAMCJmzcuR/+5dbG+99Q2PbKxRJJU3ujTVQ+s16cWnaQvXTBFXrfT5lUCAAAAAyMebTq6qqNPk/RHhUPoLygcCt8t6UzLsmr7OE+twu017pY0MTLPfEkPSpobeZ9YKyU9pnCf56sk3SbpbElrJH1C0iWWZfl7eZ/HI9etlnS5pJsldUbu/x/LstjyfICVNrRFnxdleeV0HP73FAAAAINTutetn19xsn7zsTnKTnFHx3+/Zp/ed+8abStrsnF1AAAAwMCJV2W0LMs6JOmaPl57xKTRsqw6SbdGHsea5yVJL/V1jYfdu1bSRe/lXsQfmxcCAIChbvnMIs0Zm63b/7lFq3ZWS5J2Vbbo/fet0efPm6zPLB4vlzMutSIAAABAQuLbLhJC7OaFxdkpNq4EAACg/+Sne/XgJ+fpBx+YKa87/FW8M2jpZ8/u1OW/fkXvVDbbvEIAAACg/xBGIyGUUhkNAACGCWOMPnbGWP37lsU6uTgzOr65pFEX371G96/arUAwZOMKAQAAgP5BGI2E0KNNRzZhNAAAGPrGj0jTo9cv0JcumCJPpD2HPxiKVknvokoaAAAAQwxhNBJCjzYdVEYDAIBhwuV06MalE/X0LYs0+7Aq6UvuXqP7VlIlDQAAgKGDMBoJIbZNBz2jAQDAcDO5IF3/un6Bbl/es0r658/t1AepkgYAAMAQQRgN2/kDIVU2+yRJxkiFmV6bVwQAADDwXE6HblgyUc/csqhHL+ktVEkDAABgiCCMhu3KG9tlWeHnBeleeVz8YwkAAIavSQXpevT6Bfry8qnvqpK+7P5XtLOCKmkAAAAMTqR+sF0pmxcCAAD04HI6dP2SCeEq6dFZ0fG3Sht1yT0v654X35E/QJU0AAAABhfCaNiuJHbzQsJoAACAqEkF6Xr0s2f2qJLuDFr6xX936dJ71uiNg/U2rxAAAADoO8Jo2K4ktjI6izAaAAAg1pGqpHdWNuvyX7+ibz2xVc2+ThtXCAAAAPQNYTRsR5sOAACAY5tUkK5/Xb9A37h4mpLdTkmSZUkPrTug8+5crf9uq7R5hQAAAMDREUbDdqUNbdHnVEYDAAAcmdNhdO3i8Xr+82dpyZQR0fGKJp8+/dAGXf+njapq8tm4QgAAAODICKNhu9IePaNTbFwJAADA4DA6J0UPfnKe7v7IqcpL80TH/7O1Qufc+ZL+vP6AQiHLxhUCAAAA70YYDVsFQ5bKG7qrd6iMBgAA6BtjjN538ki9cNvZ+vBpxdHxZl9AX39sq6783Trtrmq2cYUAAABAT4TRsFVlk0+BSNVObqpHyR6nzSsCAAAYXLJSPPrZh07WXz49X+Nyu39l9vr+el30qzW664Vd6ggEbVwhAAAAEEYYDVvFtuhg80IAAID3bsGEPD37ubN049IJcjmMJMkfDOmuF97RRb96Wa/sqbF5hQAAABjuCKNhq9L62H7RhNEAAAAnwut26ksXTNVTNy/SKaOzouN7qlt11f+u1y1/3cQGhwAAALANYTRsVVLfFn1Ov2gAAID4mFaUoUevX6Dvvm+GUmPaoD25uUzLfvGSfr9mnwLBkI0rBAAAwHBEGA1b9WjTQRgNAAAQN06H0ScWjNOKLy7RpSePjI63dAT0/ae36ZJ71mjD/jobVwgAAIDhhjAatiqpj+0ZnXKUKwEAAPBeFGR4dc9HTtWfr52vCSNSo+M7Kpr1od+s0xf+sVk1LR02rhAAAADDBWE0bBVbGU3PaAAAgP6zcGKe/nPrWbp9+RQlu7tbdzz6RomW3bFKD6/br2DIsm+BAAAAGPIIo2Eby7J6bGA4ijAaAACgX3lcDt2wZKJe+MLZWj6jMDre5Avom0+8rQ/ct1ZvHmqwcYUAAAAYygijYZuaFr86AuGNc9K9LmV43TavCAAAYHgYlZWs31w9V3+8Zp7G5na3SnurtFGX3b9WX/3XFtW1+m1cIQAAAIYiwmjYhs0LAQAA7LVkSr6e+9xZuu28yUpyhf9oYFnSX187pCU/X6k/rNmnzmDI5lUCAABgqCCMhm1iW3QUs3khAACALbxup245Z5L++/mzdc7U/Oh4ky+g7z29TcvvWq2VO6tsXCEAAACGCsJo2Kakvi36nM0LAQAA7DUmN0W//+Q8PfDx0zQupnXHnupWXfPg6/rkg69pd1WLjSsEAADAYEcYDdvQpgMAACDxnDu9QM9//mx97aKpSk9yRcdX7azW8rtW63tPbVNje6eNKwQAAMBgRRgN28S26RhFZTQAAEDC8Lgc+sxZE7Tii0v0P/NGy5jweCBk6Q9r92npHav0p1cPKBiy7F0oAAAABhXCaNgmtjKaNh0AAACJZ0R6kn5y+Ww9ddMinX5STnS8rtWvbzy+VRff/bJe2V1j4woBAAAwmBBGwxaWZamknjYdAAAAg8HMUZn6+2fO0P0fndPje9uOimZd9cB6XffwBh2obbVxhQAAABgMCKNhi6b2gFo6ApIkr9uhnFSPzSsCAADA0RhjdNGsIr34hbP1xfMnK8XjjJ577u1KnXvnS/rOk2+rrtVv4yoBAACQyAijYYuShrbo81FZyTJdjQgBAACQ0Lxup25aNkkrvrBEHzx1VHS8M2jpj6/s19k/W6n7V+2WrzNo4yoBAACQiAijYYvYzQuLs1NsXAkAAADei8JMr+688hQ9fuNCzRuXHR1v7gjoZ8/u1NI7VumRDYfY5BAAAABRhNGwRY9+0WxeCAAAMGidMjpL/7juTP3u6rkaPyI1Ol7e6NOX/rlFF9/9slbtrJJlEUoDAAAMd4TRsEVpA5sXAgAADBXGGJ0/o1DPf+4s/fCymcpLS4qe21HRrE8++Lo+9vv12lraaOMqAQAAYDfCaNiiZ5sOwmgAAIChwOV06KPzx+qlLy3R586d1GOTw7W7a3XJPWv0+b+/qZL6tqPMAgAAgKGKMBq2iK2MJowGAAAYWlKTXPrcuZO16otLdNX8MXI6ujerfmxTqZbd8ZJ+8PQ21bX6bVwlAAAABhphNGwRWw0zKosNDAEAAIai/AyvfnTZLD33ucU6d1pBdNwfDOmBNfu0+KcrdOfzO9Xk67RxlQAAABgohNEYcG3+gOrbwn/gcDuN8tOTjnEHAAAABrOJ+el64BOn6e+fOUMnj86Kjrf6g7p7xW4t/ulK/XrVHrX5AzauEgAAAP2NMBoDLrZfdFFmshwxP9sEAADA0DV/fK4ev2GBfvOxuZqUnxYdb2zv1E+f3aGzfrZKf1y7Tx2BoI2rBAAAQH8hjMaAi+0XPSqLftEAAADDiTFGy2cW6tnPnaVfXnmyxuR0t2yraenQd57apmV3vKR/vH5IgWDIxpUCAAAg3gijMeDKG33R50VZXhtXAgAAALs4HUaXnVqsF79wtn502SwVZnR/LyxtaNftj27R+b9crac2lykUsmxcKQAAAOKFMBoDrjymMnpkJpXRAAAAw5nb6dBV88do1ZeW6BsXT1Nuqid6bm9Nq27+6yZdfM8a/XdbpSyLUBoAAGAwI4zGgCujMhoAAACH8bqdunbxeK2+fam+eP5kpXtd0XPby5v06Yc26NJ71+j5tysIpQEAAAYpwmgMuPJGKqMBAADQu9Qkl25aNklrbl+mG5ZMULLbGT23tbRJn3l4oy66e42e3VpO+w4AAIBBhjAaA668gcpoAAAAHF1milu3L5+q1bcv1acWnSSvu/uPLtvLm/TZP72hi+5+Wc9sIZQGAAAYLAijMaAsy1JZTGV0EZXRAAAAOIoR6Un65iXT9fLty/TpxSf1qJTeUdGsG//yhi64a7We3FymIKE0AABAQiOMxoBqaOuUrzMkSUr1OJUR0wsQAAAAOJIR6Un6+sXT9fKXl+q6s8crxdMdSr9T1aJb/rpJ5//yJT2+qZRQGgAAIEERRmNA9aiKzkqWMcbG1QAAAGCwyUtL0lcvnKY1Xw73lE6NCaX3VLfqc39/U+fd+ZL+9UaJAsGQjSsFAADA4QijMaB69IvOpF80AAAA3pucVI9uXz5Va768TDcvm6j0pO5f3O2tadVt/9isJXes0kPr9svXGbRvoQAAAIgijMaAKm8ijAYAAED8ZKd69IXzp2jNl5fplnMmKT2mDVxJfbu+9cTbWviTFbpv5W41tnfauFIAAAAQRmNAlTeweSEAAADiLzPFrdvOm6w1X16m286brOwUd/RcbatfP39upxb+ZIV+/O/tqoopkAAAAMDAIYzGgCpv7P7iPzKLymgAAADEV2ayW7ecM0lrv7JM3750ukbG/BqvpSOg367eq0U/Xamv/ust7a9ptXGlAAAAww9hNAZUGZXRAAAAGAApHpeuWXiSXrp9qe644mRNzE+LnvMHQ/rrawe17BerdONf3tDW0kYbVwoAADB8uI59CRA/VEYDAABgILmdDn1obrE+eOoovbC9Uvev2qM3DzVIkkKW9MyWcj2zpVxnTR6h684arwUTcmWMsXnVAAAAQxNhNAZMKGSpojF2A0MqowEAADAwHA6j82cU6rzpBXp1b51+/dIerd5VHT2/ele1Vu+q1tTCdF27eLwuPblISS6njSsGAAAYemjTgQFT2+qXPxiSJGV4XUpN4u9CAAAAMLCMMTpzQq4e+n+n6+mbF+ni2UVyxBRC76ho1hcf2axFP12pe1e8o/pWv32LBQAAGGIIozFgyhu7+0WPzKIqGgAAAPaaOSpT9101Ryu+sEQfP3Oskt3dldDVzR264/ldOvMnL+obj7+lvdUtNq4UAABgaCCMxoApa4ht0UG/aAAAACSGcXmp+t77Z2rdV5fp9uVTVJCRFD3n6wzpT68e1LJfvKRr/+91rdtTK8uybFwtAADA4EWfBAyY2MroIiqjAQAAkGCyUjy6YclEXbtovJ7eUqYHXt6nbeVN0fMvbK/SC9urNGNkhq5dfJIunjVSHhf1PQAAAH3FNycMmPLYzQszqIwGAABAYvK4HPrgnGI9c8si/eXT83XO1Pwe598ua9Ln/75Zi366Qne9sEtVTb4jzAQAAIBYVEZjwPQIo6mMBgAAQIIzxmjBhDwtmJCnPdUt+sOafXr0jRL5OsObclc1d+iuF97RvSt268JZRfrkgrGaMyZbxphjzAwAADA8URmNAVPeELOBIT2jAQAAMIhMGJGmH142S6985Rx94bzJyk/v7isdCFl6anOZLv/1Ol1yzxr94/VD8nUGbVwtAABAYiKMxoChMhoAAACDXU6qRzefM0lrv7JM9151quaNy+5x/u2yJt3+6Bad8eMX9eP/bNehujabVgoAAJB4aNOBAREMWaqI6aVXRGU0AAAABjG306FLZo/UJbNH6u2yRj287oAef7M02sKjoa1Tv31pr/539V6dM61AnzhznBZOzKWFBwAAGNYIozEgqps7FAxZksLVJF630+YVAQAAAPExY2SmfnL5bH3lwqn6x4ZDemjdAZXUh1vUhSzpv9sq9d9tlRo/IlVXnT5Gl88pVnaqx+ZVAwAADLy4tekwxhQbY/5gjCkzxnQYY/YbY+4yxmQf++4e8+RE7tsfmacsMm9xL9fmGmOuNcY8ZozZbYxpN8Y0GmPWGGM+ZYx51+czxowzxlhHefztRP5zQO/KGrv7RVMVDQAAgKEoK8Wjz5w1QS99aake+PhpWjwpr8f5vdWt+sEz2zX/Ry/q1r9t0vq9tbIsy6bVAgAADLy4VEYbYyZIekVSvqQnJO2QdLqkWyUtN8YstCyrtg/z5EbmmSxphaS/SZoq6RpJFxtjzrQsa2/MLVdI+rWkckkrJR2UVCDpg5IekHShMeYKq/dveJslPd7L+NZjf2Icr/KG2BYd9IsGAADA0OV0GJ07vUDnTi/QnuoWPbzugP65sUQtHQFJkj8Y0hNvlumJN8s0YUSqPkK1NAAAGCbi1abjfoWD6Fssy7qna9AYc6ekz0v6oaTP9mGeHykcRP/SsqzbYua5RdKvIu+zPOb6XZLeJ+kZy7JCMdd/TdJrki5XOJh+tJf3etOyrO/05cPhxJXHVEaPzKIyGgAAAMPDhBFp+s77ZuhLF0zR01vK9Jf1B7W5pDF6fk+kWvpnz+3URTMLddX8sZo3Lpve0gAAYEg64TYdxpjxks6XtF/SfYed/rakVklXG2NSjzFPqqSrI9d/+7DT90bmvyDyfpIky7JWWJb1VGwQHRmvkPSbyMslx/Fx0E/KqIwGAADAMJaa5NKV88boiZsW6embF+mj88coLam7NsgfCOnxN8v04d+u03m/XK3fr9mnhja/jSsGAACIv3j0jF4WOT7fSyjcLGmtpBRJZxxjnjMlJUtaG7kvdp6QpOcjL5f2cV2dkWPgCOdHGmOuM8Z8LXKc3cd58R6U0zMaAAAAkCTNHJWpH142S+u/do5+/MFZml2c2eP87qoWff/pbTr9Ry/q5r9u0upd1dHNwAEAAAazeLTpmBI57jrC+XcUrpyeLOnFE5xHkXmOyhjjkvTxyMtnj3DZeZFH7H2rJH3CsqyDx3qPyPUbj3Bqal/uH07KGmMrowmjAQAAgNQklz5y+hh95PQx2lraqD+vP6gn3yxVqz8oKVwt/dTmMj21uUxFmV5dPqdYH5pbrHF5R/3RKQAAQMKKR2V011/jNx7hfNd41gDNI0k/kTRT0r8ty3rusHNtkr4vaa6k7MjjbIU3QFwi6cVjtRTB8avo0TOaNh0AAABArJmjMvXjD87S+q+fqx9dNkuzRvWsli5v9Onelbu15I5V+vBv1+mRDYfU2nGkH4ECAAAkpnhtYHg0XTtvnOjvyvo0T2Szwy9I2qFwD+oeLMuqkvStw4ZXG2POl7RG0nxJ1yq8YeJRWZY19whr2ChpzrHuHy46gyFVNXdIkoyRCjKojAYAAAB6k5bk0lXzx+iq+WO0vbxJj2wo0eNvlqqutbt/9Gv76vTavjp958m3dfHsIl1x2midNpZNDwEAQOKLR2V0V8Vy5hHOZxx2Xb/NY4y5UeEQeZukpZZl1R3jPaMsywpIeiDy8qy+3odjq2zyyYr8FUJeWpI8rnj8YwcAAAAMbdOKMvStS6fr1a+eo998bI7OmZovR0ze3OoP6h8bSnTFb9Zp2S9e0n0rd6usof3IEwIAANgsHpXROyPHI/VynhQ5HqkXdFzmMcZ8TtIvJW2VdE6kAvp4VUeOtOmIo/KYftEj6RcNAAAAHBePy6HlM4u0fGaRqpp8+temUj2y4ZD2VLdGr9lX06qfP7dTdzy/U/NPytFlp47S8plFykx227hyAACAnuIRRq+MHM83xjgsywp1nTDGpEtaKKld0qvHmOfVyHULjTHplmU1x8zjUHgTxNj3U8z5LyvcJ/pNSedZllXzHj/LGZHj3vd4P3oRW51RlEm/aAAAAOC9ys/w6rNnT9B1Z43XpkMNemTDIT21uVwtkf7RliW9urdOr+6t0zefeFvnTsvXB04ZpSVT8vmFIgAAsN0JfxuxLGuPpOcljZN042Gnv6twlfFDlmVF/9reGDPVGDP1sHlaJD0cuf47h81zU2T+5yzL6hEUG2O+qXAQvVHhiuijBtHGmPnGGE8v48skfT7y8k9HmwPHJ7YyuiiLymgAAADgRBljNGdMtn78wdl6/evn6s4Pn6zFk/J6tPHwB0L691sV+szDG3X6j17QNx5/SxsP1MmyTnQ7HwAAgPcmXhsY3iDpFUl3G2POkbRd4Y0AlyrcVuPrh12/PXI8fIeNr0laIuk2Y8wpkl6TNE3S+yVV6bCw2xjzCUnfkxSU9LKkW3rZtGO/ZVl/jHn9U0kzjDGrJJVExmZLWhZ5/k3Lsl451gdG35XHVEaPpDIaAAAAiKtkj1MfnFOsD84pVmWTT09tLtNjm0r1dllT9JqGtk796dWD+tOrBzUmJ0UfOGWk3n/qKE0YkWbjygEAwHATlzDasqw9xpjTFA6Gl0u6SFK5pLslfbevGwlallVrjDlT0rclfUDSYkm1kh6U9C3LskoOu+WkyNEp6XNHmPYlSX+Mef2wpMskzZN0oSS3pEpJ/5B0r2VZL/dlrei7MiqjAQAAgAFRkOHVtYvH69rF47WrslmPbyrVE2+WqTSmQORgXZvuXrFbd6/YrZmjMnTp7JG6eHaRirNTbFw5AAAYDgw/0Yo/Y8zGOXPmzNm4caPdS0kIl9zzsraWhqsyHr1+geaOzbZ5RQAAAMDwEQpZen1/nR5/s1RPbylXsy/Q63WnjsmKBtMFGRSRAACAsLlz5+qNN954w7KsuSc6V7zadABHVN4QUxmdyZdaAAAAYCA5HEbzx+dq/vhcffvSGVq1s0qPbSrVih1V6gx2FydtOtigTQcb9P1ntmneuBxdevJIXTizUHlpSTauHgAADCWE0ehXvs6galv9kiSHkfLT+SILAAAA2MXrdmr5zCItn1mkxvZOPf92hZ7eUq41u2sUDIWDacuSXttXp9f21enbT2zVggl5uvTkIl0wo1BZKe/aCx4AAKDPCKPRryqbuquiCzK8cjkdNq4GAAAAQJfMZLeuOG20rjhttOpa/Xp2a4We3lKmV/fWKpJLK2RJa3bXaM3uGn39sa1aNClPy2cU6rzpBcqlYhoAABwnwmj0qzJadAAAAAAJLyfVo6vmj9FV88eoqtmn/7wVDqZf318fvSYQsrRqZ7VW7azW1x57S/NPytXymYW6YEahCvmuDwAA+oAwGv2qvLF71+6irGQbVwIAAACgL/LTvfrEgnH6xIJxKm9s1zNbyvXUlnJtPh3rD8cAACAASURBVNQQvSZkSev21mrd3lp9+8m3deqYLC2fUagLZxZpTG6KjasHAACJjDAa/aq8sbsyeiTVEgAAAMCgUpSZrGsXj9e1i8frUF2bnnu7Qs9urdDGg/Wyuvc+jG5++OP/7ND0ogwtn1moC2cWamJ+mowx9n0AAACQUAij0a/KGmIqozOpjAYAAAAGq9E5KdFguqrJp+e2Veq5rRVat7c2uvmhJG0rb9K28ibd+d9dGp+XqvOmF+jc6QWaMyZbTgfBNAAAwxlhNPpVj8roLCqjAQAAgKEgP8Orq88Yq6vPGKv6Vr9e2F6pZ7dW6OV3auQPhqLX7a1p1W9X79VvV+9VTqpHy6bm69xpBVo8KU+pSfxxFACA4Yb/90e/ojIaAAAAGNqyUz264rTRuuK00Wr2dWrlzmo9u7VcK3dUq70zGL2urtWvf24s0T83lsjjcmjhhFydN71Q50zLV0EGhSsAAAwHhNHoV7GV0UVURgMAAABDWrrXrfedPFLvO3mkfJ1BvbKnRv/dVqkXtlepurkjep0/ENLKndVaubNaekw6uThT504Lt/OYWphOn2kAAIYowmj0mzZ/QI3tnZIkt9MoLzXJ5hUBAAAAGChet1PLphZo2dQC/TBkaUtpo17YVqkXtldqR0Vzj2s3lzRqc0mjfvHfXSrK9GrJlHwtnTJCCyfSzgMAgKGE/1dHvylr6K6KLsjwysFmJQAAAMCw5HAYnTI6S6eMztIXL5iig7VtemF7OJhev6+uxwaI5Y0+/fW1g/rrawflcTp0+kk5WjJlhJZMydeEEalUTQMAMIgRRqPflDd294seSb9oAAAAABFjclP0/xadpP+36CQ1tnVq1a4q/XdbpVbvqlaTLxC9zh8Mac3uGq3ZXaMfPLNdo3OStXRKvpZOydcZ43OV7HHa+CkAAMDxIoxGv6FfNAAAAIBjyUxx6/2njNL7TxmlQDCkTYcatHJHlVburNb28qYe1x6qa9dD6w7ooXUHlORy6IzxuTp78ggtnpSniflpVE0DAJDgCKPRbypjwuhCdscGAAAAcAwup0PzxuVo3rgc3b58qioafXppV5VW7qjWmt01aunorpruCIT00q5qvbSrWlL4zxyLJ+Vp8eQRWjQxTzmpHrs+BgAAOALCaPSbyuaePaMBAAAA4HgUZnp15bwxunLeGPkDIW04UKdVO6u1ckeV3qlq6XFtRZNPj2ws0SMbS2SMNHNkphZPytOiSXmaOzZbSS5aegAAYDfCaPSbisaO6PPCTMJoAAAAAO+dx+XQggl5WjAhT1+7aJpK6tv00q5qrXkn3FO6OabXtGVJb5U26q3SRt2/ao+S3U6dMT5HiyeN0KJJeZpESw8AAGxBGI1+U9lEZTQAAACA/lGcnaKPzh+rj84fq0AwpC2ljXp5V41efqdamw41KBiyote2dwa1cme1Vu4Mt/TIS0vSmRNytSDyGJOTQjgNAMAAIIxGv6noEUYn2bgSAAAAAEOZy+nQnDHZmjMmW7eeO0lNvk6t21OrNe+Ew+n9tW09rq9p6dBTm8v01OYySdKorOSYcDqPX3YCANBPCKPRLzqDIdW0dLfpyE/nyxwAAACAgZHhdeuCGYW6YEahJOlgbZte3l2tl3fV6NV9tWpo6+xxfWlDu/65sUT/3FgiSRqflxoJp/N0xvgc5aZRXAMAQDwQRqNf1LR0yIr8Ki4vzSOPy2HvggAAAAAMW2NyU/TR3HBLj1DI0vaKJq3bU6tX9tRq/d5atfqDPa7fW9OqvTWt+vP6g5Kkiflpmn9Sjk4/KUfzT8qlchoAgPeIMBr9oqKRftEAAAAAEo/DYTRjZKZmjMzUtYvHqzMY0luljZFwukYb9terIxDqcc/uqhbtrmqJhtNjc1N0+rgczR+fq/kn5ag4O5me0wAA9AFhNPoFmxcCAAAAGAzcMf2mb1w6Ub7OoDYdbNC6PTV6ZU+tNpc0qDNo9bjnQG2bDtS26ZFIW4+iTG+kcjpXp5+UowkjUgmnAQDoBWE0+gWV0QAAAAAGI6/bqTMn5OrMCbm6TVK7P6hNh+q1fm+dXttXpzcOvrtyurzRp8ffLNPjb4Y3RMxOcWvu2GzNHZuj08Zla9aoTHndThs+DQAAiYUwGv2isrl788JCwmgAAAAAg1Syx6kFE/K0YEKeJKkjENRbJY1avy8cTm/YX/euntP1bZ16YXuVXtheJUlyO41mjsrUaTEBdR6bIgIAhiHCaPSLypjK6MJMvmQBAAAAGBqSXE6dNi5Hp43L0Y1LpUAwpG3lTXptX51e3VunjQfqVN/W2eOezqClTQcbtOlgg/735X2SpHG5KZo7Nkdzx2ZrztgsTcpPl9NBaw8AwNBGGI1+URHTMzqfymgAAAAAQ5TL6dDs4izNLs7StYvHy7Is7alu1cYDddqwv14bD9Rrb03ru+7bX9um/bVtevSNcN/pVI9Ts4uzdOqYLJ0yOkunjMlSfjp/lgIADC2E0egXsWE0bToAAAAADBfGGE3MT9PE/DRdOW+MJKm2pUNvHGzQhgN12ri/XltKGuUP9uw73eoPat3eWq3bWxsdK85O1imjs3TqmGydMjpLM0Zm0HsaADCoEUajX1Q10TMaAAAAACQpNy1J500v0HnTCySF+05vLW2MVk6/eahBVTH77nQpqW9XSX27nt5SLince3p6UYZOHp2lWaMyNas4UxNHpMnldAzo5wEA4L0ijEbctXQE1NIRkCR5XA5lpbhtXhEAAAAAJI4klzPSLzpHkmRZlsobfZG+0uFw+q3SRnUEelZPdwYtbS5p1OaSxuiY1+3QjJGZ4XB6VKZmF2dq/Ig0+k8DABISYTTirqKxZ4sOY/gSBAAAAABHYozRyKxkjcxK1sWziyRJncGQdpQ3a9Oher15sEGbDjVoXy+9p32dIW08EK6w7pLicWrGyAzNGpWl2cXhCuqTclPlIKAGANiMMBpxVxnTL7ogI8nGlQAAAADA4OR2OjQrEiR//MzwWH2rX2+WNGhrSaO2lDbqrZLGHvv1dGnzB/X6/nq9vr87oE5LcmnGyAzNLs7UzFGZml2cpbE5KQTUAIABRRiNuOsZRtMvGgAAAADiITvVo6VT8rV0Sn50rKrZp62ljdpS0qitpeEWHtW99J9u6Qho/b46rd9XFx1LT3JpalG6phVlaHpRhqYVZWhKYTqbJAIA+g1hNOIu9m/m2bwQAAAAAPpPfrpXy6Z6tWxqQXSsssmnt6LV0+H+0zUt/nfd29wReFcFtcNI40ekaVpRhqYVpWt6JKgekZ5EC0YAwAkjjEbcVcb2jM4kjAYAAACAgVSQ4VXBdK/OnR4OqC3LUkWTL1o9vaWkUW+VNqqu9d0BdciSdle1aHdVi57a3D2el+aJBNTdVdTjR6TK7XQM1McCAAwBhNGIu9jK6HwqowEAAADAVsYYFWUmqygzWRfMKJQUDqgrmzq0rbxR28ubta28SdvLmrSvtlWW9e45alr8evmdGr38Tk10zON0aHJhmqYVhtt7TCpI15SCdBVkUEUNAOgdYTTirrKpuz8ZbToAAAAAIPEYY1SY6VVhZs8WH23+gHZWRMLp8iZtL2/W9vImtfmD75rDHwxpa2mTtpY29RhP97o0uSA98kiLPs9L8xBSA8AwRxiNuKukZzQAAAAADEopHpdOHZOtU8dkR8dCIUsH69q0vbypR0hd2tDe6xzNvoA2HqjXxgP1PcazU9y9htTZqZ5+/UwAgMRBGI24CoYsVcXs3JyfkWTjagAAAAAAJ8rhMBqXl6pxeam6cFZRdLyhza/t5c3aUdGkXZUt2lXZrF2VzWr2BXqdp76tU+v31Wn9vroe43lpSZpSmKZJ+emaVJCm8XlpmpCfqhFptPsAgKGGMBpxVdvSoWAo3GAsK8Utr9tp84oAAAAAAP0hK8WjMyfk6swJudGxrl7UOyub9U5ls3ZWNGtXVYt2VzartZdWH5JU09Khmt0dWru7tsd4utel8SPSNGFEqiaMSIs8UjU2N1UeFxsnAsBgRBiNuKqgRQcAAAAADFuxvajPnjwiOh4KWSptaNc7Vc3hKuqKZu2qatY7lS3qCIR6navZF9DmQw3afKihx7jTYTQmJ0Xj81I1IT8cUI+PhNU5tPwAgIRGGI24it28sIAwGgAAAACgcKuP0TkpGp2T0mPDxGDI0qG6tmiLj73VrdpT3aI91a1q6ei93UcwZGlfTav21bTqxR1VPc5lp7g1fkSaxkfaiozLTdXY3BSNy0tVWhIRCADYjX8TI66ojAYAAAAA9JUzph/1+TMKo+OWZam6uUO7I8H03shxT1XLETdOlMJ9qXvbPFEK96Yel5uisbmpOikvfByXm6qxeSnK8Lr75fMBAHoijEZcVTZ2h9EFmYTRAAAAAIDjZ4xRfoZX+RleLZiQ1+Ncuz+ofTVdFdQt0WrqvdWtau/svS+1FOlN3dKhDb0E1bmpnnAFdW44GI8+z01VZgpBNQDEC2E04iq2MrogI8nGlQAAAAAAhqJkj1PTR2Zo+siMHuOhkKXyJp/2VLVof22r9te06UBtq/bVtupQXZs6g9YR56xt9au21a83Dja861yG1xVuMZKdotE5yRqTk6LiyOvi7GR53c64f0YAGKoIoxFXlbTpAAAAAADYwOEwGpWVrFFZyTpLI3qcC4YslTW0h0Pq2jYdqGmNPj9Y2yZ/sPdNFCWpyRfQ22VNerusqdfz+elJGhPphz06OzkaVI/JTVFhhldOh4nr5wSAwYwwGnFV2aMymjAaAAAAAGA/Z8wGiosn9TwXDFkqb2zXgdo27a9t1YHaNu2radWByPOOwJGDakmqau5QVXPv7T/cTqORWcnRKuqRkbC861iY6ZXH5YjnRwWAhEYYjbiqiOkZXUjPaAAAAABAgnM6jIqzU1ScnaKFE3v2p+7aSPFQfZsO1bXrUF2bDtW36WBd+HV5Y7tCR+7+oc6gpQO1bTpQ29breWPCldVd4XRsUN11zEh2yRiqqwEMDYTRiJt2f1BNvoCk8N/+5qR4bF4RAAAAAADvXexGinPHvvt8ZzCk8gZfJKyOhNT14dC6pL5NNS3+o85vWVJlU4cqmzq0qZd+1ZKUluTSyCxvj5C6MMOrokyvCjLDxxQP8Q6AwYF/WyFuYlt05Kd75aAvFgAAAABgCHM7HRqTG+4P3Zs2f0AlkXC6tKFdpQ3tKmvwqbS+TWUNPlU2+2QdpbJaklo6AtpV2aJdlS1HvCbD61JhpleFmckqzEhSYWayijK9KszwqjASWGcmu6mwBmA7wmjETUWPftFJNq4EAAAAAAD7pXhcmlyQrskF6b2e7wyGVNHoi4TU7Sqtb1dZY7tKG3zR1+2dwWO+T5MvoCbf0QPrJJcjHFBHQ+pk5acnKT8jSfnp3uhzqqwB9Cf+DYO4ia2Mpl80AAAAAABH53Y6ohsr9sayLDW0dUarqkvr21XR5FN5o0+VjT6VN7WrsrFD/uDRN1mUpI5ASPtr27T/CP2ru6R6nMrP8GpEelI4oE73RgLrns+ptAbwXhBGI25iNy/MTyeMBgAAAADgRBhjlJ3qUXaqRzNHZfZ6jWVZqmv1q7zRp4pGnyqawsfyRp8qm3wqb2xXRaNPrf5jV1hLUqs/qH01rdpX03rU6zxOh0akJ3WH1odVWOenhwPt7BSPPC7HcX92AEMTYTTiprKpI/qcymgAAAAAAPqfMUa5aUnKTUs6YmAtSc2+zmhY3RVcVzd3qKrZp6rmDlU1dai6uW9V1pLkD4aiFdvHkpnsVm6aR3mpScpN84QfqUnKS/OE154aPualeZThdbMHFTCEEUYjbnq06cggjAYAAAAAIFGke91K97o16Qj9q6VwlXVje2d3ON3iU1VTR/h1c4eqmroC7A61dAT6/N6N7Z1qbO/U3uqjV1tLksthlBMTTncF1T3D7K4A20OPa2CQ4X+xiJueGxgSRgMAAAAAMJgYY5SV4lFWiueImy52afMHosF0OLCOqbBuCQfXNS1+1bV2KGT1fQ2BkBUNv/siyeVQTmp4zdkp7nBbkxS3ciKfI3zOrZxUj7JTwi1PUj1O+l0DNiGMRtzE9owuyEiycSUAAAAAAKA/pXhcGpvr0tjc1KNeFwpZamjvVG1Lh2pa/Kpt7VBtiz/8ujV8rG3xq7bVr5qWDjX7+l5xLYU3ZiyP9MjuK7fThIPplO6gOhxcu6Pj2ZHnWSkeZSa7leF1yeWk9zVwogijEReWZamqOaZNBz2jAQAAAAAY9hyRths5qR5NKjj29R2BoOpa/aptCYfTtTEBdm9htj/Qtx7XsTqDx1d93SUtyaXMZPe7Hym9jEUeWSnh9ihO+mADkgijESd1rX51BsO/u0n3uujZBAAAAAAAjluSy6mizGQVZSYf81rLstTeGQ6vG9o6VdfqV32bX/WtftW3dYaft3VGXnePt3cG39PaWjoCaukI9GnTxsOle13vCqkzk93KSHYrw+tWuteldK8r8rz7dbrXrfQkF5s6YsggMURcVLB5IQAAAAAAGEDGGKV4wgVxxdl9v8/XGVR9mz8aYscG1eExv+oiIXbX5otNvk5Zx9H7+nDNvoCafQGV1B9/kC1J6Ukx4XRMUJ2RHDsWbidyeKid5nUp1eOiOhsJgTAacVHZRIsOAAAAAACQ+LzuvldfdwmFLDX7AtFwuuvR0B4TWMeOt3U/P94+2L1p7giouSMgHUdv7MMlu51KjYTaqUlOpXpcSktyKTUpHFinJYVD69QkZ6/j4TGn0rwuJbmcJ/yZMDwRRiMuKpu6+yzlpxNGAwAAAACAocPhMOHe0Cnu4743GLLU7Os9qO4Kq5t9PY9N0dfh1iDx0N4ZVHtnUDUtx9cruzdup1Hqu0Jqt9IiIXdqUng8JcmpFLdTKR6Xkj1OpSY5lex2KcXjDD+SXEpxO5XscSrJ5ZAxVG8PdYTRiIuKxtjK6CQbVwIAAAAAAJA4nA6jrBSPslI87+n+YMhSS0dATYcH1x2damqPDbDDz5sOC7dbfAG1+t9bn+wj6QxaamgLB+vx4jBSaiS0TvE4lexxKdXjjL4Ot2TpPpficUbOuyJj4SC867nX7VRy5JHkctB3O0EQRiMuKukZDQAAAAAAEHdOh4lufPhehUKW2jqDau0IV1u3doQfLR0BtfoDaukIhkPrrrEe451q7QhGr23tCKgzeAINtI+0RiumHUk/SHI5lOyJCajdTiW7u8eSYsJrr9sRPnq6XsccPU55I3PNHJlJyH2cCKMRF7EbGBYQRgMAAAAAACQMh8MoLdI6oyDjxOfrCHSF18GYQDsSZvu6ngfV1hlQuz+oNn9Qbf5A5BhUuz+oVn/3uXZ/UP5g6MQXdtQ1h9QRCKlB8avm3vOji+I213BBGI24iG3TQRgNAAAAAAAwdCW5nEpKcyo3LX5zdgZD0WA6Nrhu8x8h0I5UekfPdQbV1hE+74v0x27vDD/3dcY/6Pa4HHJSFX3c4hZGG2OKJX1P0nJJuZLKJT0u6buWZdUfxzw5kr4l6QOSiiTVSnpW0rcsyyqJ13sbY6ZL+o6kJZIyJB2Q9DdJP7Esq72v60VYVXN38/vCTMJoAAAAAAAA9J3b6VBmsuOE2pEcSShkqSMQ6hFSt/uD6ggE1e4PHRZch8+1R0Lsnq+7jwTR701cwmhjzARJr0jKl/SEpB2STpd0q6TlxpiFlmXV9mGe3Mg8kyWtUDgcnirpGkkXG2POtCxr74m+tzFmfmR+t6R/SjokaZnCIfg5xpj/396dB1lW1mcc/z6DAWEQXMBQZSxglGVKkqAo2ySsOuIGqFAxRjQkYohYIGqB5YZLTDAVRQeNmqCOWwIKCiZSMCqbYtRSo1kYBRxGhbAIYwBhgAi//HFOa3udZqaXc0/b5/up6jr0ec/yG+qte/t97nvfc2hVzf7RogNx78/vZ91d9wHNOkbbbe0DDCVJkiRJkjQ/LFqUZm3ozTfjEX0XM3CL5ug6f08TBp9YVUdW1Wur6hDgDGA34O2beJ2/pgmiz6iqQ9vrHEkTLD+6vc+s7p1kM+AjwFbAUVX1wqo6FdgHOA9YBpw8nX/80N1yxy9z++233sJPhiRJkiRJkiT9mlmH0UmWAMuBtcD7RppPA+4CjkmyeCPXWQwc0x5/2kjze9vrP72932zufSCwFLiiqj43sbOqHgBOaX89PomJ6ib61YcXOitakiRJkiRJ0q+bi5nRh7TbVW2g+wtVdSdwJc0s5H03cp39gC2BK9vzJl/nAWBV++vBs7z3xDkXjRbQLgFyNbAjsGS0XRt28x0+vFCSJEmSJEnSg5uLMHq3dnv1FO3XtNtdO7jOuM7ZoCTf2tAPzTrXg3HT7b8Mo314oSRJkiRJkqQNmYswett2e/sU7RP7H97BdcZ1jh6EM6MlSZIkSZIkbcxDxnCPibWXq4frdHpOVe21wQs0s6OfNI17/kY7Ys/HsNN2i7n59nvYd8kj+y5HkiRJkiRJ0jw0F2H0xEzibado32bkuLm8zrjO0YPY4zHbssdjpvrfKUmSJEmSJElzs0zH99vtVGss79Jup1qjeTbXGdc5kiRJkiRJkqRZmIsw+tJ2uzzJr1wvycOAZcB64Gsbuc7X2uOWtedNvs4iYPnI/WZ670va7WGjBSRZQhNS/xBYs5F6JUmSJEmSJEmbaNZhdFX9AFgF7AScMNL8FmAx8LGqumtiZ5Ldk+w+cp2fAR9vj3/zyHVe0V7/4qpaM+mcad8buBxYDRyQ5PBJNS0C3tH++oGqmu0a15IkSZIkSZKk1lw9wPDlwFeBFUkOpQl79wEOplnu4vUjx69utxnZ/zrgIOBVSfYEvgEsBY4AbuHXA+dp37uq7k9yLM0M6XOTnAv8CDgUeDJwJXDGNP7tkiRJkiRJkqSNmItlOiZmKD8ZWEkTBL8aeBywAtivqm7bxOvcBuzXnvf49jr7AB8B9mrvM+t7V9XXgacAF9As/3EyzQMN3wo8raru3bR/uSRJkiRJkiRpU8zVzGiq6sfAsZt47OiM6Mlt64CT2p85v/ekc64Cjp7OOZIkSZIkSZKkmZmTmdGSJEmSJEmSJD0Yw2hJkiRJkiRJUucMoyVJkiRJkiRJnTOMliRJkiRJkiR1zjBakiRJkiRJktQ5w2hJkiRJkiRJUucMoyVJkiRJkiRJnTOMliRJkiRJkiR1zjBakiRJkiRJktQ5w2hJkiRJkiRJUucMoyVJkiRJkiRJnTOMliRJkiRJkiR1zjBakiRJkiRJktQ5w2hJkiRJkiRJUucMoyVJkiRJkiRJnTOMliRJkiRJkiR1zjBakiRJkiRJktQ5w2hJkiRJkiRJUucMoyVJkiRJkiRJnTOMliRJkiRJkiR1LlXVdw0LTpLbttxyy0cuXbq071IkSZIkSZIkacZWr17N+vXr11XVo2Z7LcPoDiS5DtgGWNtzKeO0e7v9Xq9VaOjsh5ov7IuaD+yHmg/sh5ov7IuaD+yHmg/sh5qJnYA7qmrn2V7IMFpzIsm3AKpqr75r0XDZDzVf2Bc1H9gPNR/YDzVf2Bc1H9gPNR/YD9U314yWJEmSJEmSJHXOMFqSJEmSJEmS1DnDaEmSJEmSJElS5wyjJUmSJEmSJEmdM4yWJEmSJEmSJHUuVdV3DZIkSZIkSZKkBc6Z0ZIkSZIkSZKkzhlGS5IkSZIkSZI6ZxgtSZIkSZIkSeqcYbQkSZIkSZIkqXOG0ZIkSZIkSZKkzhlGS5IkSZIkSZI6ZxgtSZIkSZIkSeqcYbRmJcnvJPlwkv9Jcm+StUneneQRfdemhS/Jo5K8NMlnk1ybZH2S25N8JcmfJ/E1Tr1JckySan9e2nc9GpYkf5jkvCQ3tu/PNyZZleSZfdemYUjyrLbPXd++P69J8ukk+/VdmxaWJEclOTPJl5Pc0b7vfmIj5+yf5MIk65LcneQ/krwyyWbjqlsLy3T6YZJdkpya5JIkP05yX5Kbk1yQ5OBx166FZSaviSPnf2jSGObxXdaq4XpI3wXoN1eSxwFfBR4NXAB8D9gbOAk4LMmyqrqtxxK18B0NvB+4EbgU+BHw28DzgLOAZyQ5uqqqvxI1REkeC5wJ/AzYuudyNDBJ3gC8DbgV+Fea18jtgCcCBwEX9lacBiHJO4BTgNuA82n64uOBI4DnJ3lxVW3ywFjaiDcAv0/znns9sPuDHZzkCOA84B7gHGAd8BzgDGAZzd+X0nRNpx++Dfgj4Cqa9+R1wG7A4cDhSU6qqhXdlqsFbFqviZMleQ7wZziGUcdiRqOZSnIxsBw4sarOnLT/XcDJwAer6vi+6tPCl+QQYDHw+ap6YNL+HYBvAI8Fjqqq83oqUQOUJMAXgJ2BzwCvAY6rqrN6LUyDkORo4FPAF4HnVdWdI+2/VVX/10txGoT2PfgG4CfA71XVLZPaDgYuAa6rqiU9lagFpu1X1wPXAgfSTFD4ZFW9aAPHbtMety2wrKq+2e5/KE3f3A/446o6e0zla4GYZj/8U+C7VfXvI/sPpPkbsoCdqurGruvWwjOdvjhy3vbAfwKXATu05+5SVdd2WrAGya+wa0aSLKEJotcC7xtpPg24CzgmyeIxl6YBqapLqupfJgfR7f6bgA+0vx409sI0dCcChwDH0rwWSmPRLk30DuBu4IWjQTSAQbTGYEeaMcbXJwfRAFV1KXAnsH0fhWlhqqpLq+qaTfwm3FE0/e/siSC6vcY9NLMJAf6ygzK1wE2nH1bVytEgut1/OU0QuDmw/9xXqSGY5mviZP/Qbk+Y65qkUYbRmqlD2u2qDQSBdwJXAlsB+467MKk1Ebj8vNcqNChJlgKnA++pqiv6rkeDsz/NjPwLgZ+2a/aemuQk1+nVGF0D3AfsnWS7yQ1JDgAeRjNzX+rDxBjmog20XUHzYd7+SbYYX0nSr3AMo7FrZ+sfCRzvUqsaB9eM1kzt1m6vnqL9GpqZ07sCXxpLRVIryUOA7+ZE9wAABd5JREFUF7e/bmiwIc25tt99nGbt8tf1XI6G6Snt9mbg28DvTm5McgXN0kU/GXdhGo6qWpfkVOBdwFVJzqdZO/pxNOuhfgH4ix5L1LBNOYapqp8nuQ54ArAEWD3OwqQkOwKH0nwo4qQGjUXb794DfKKqzu+7Hg2DYbRmatt2e/sU7RP7Hz6GWqRRpwN7ABdW1cV9F6PBeBPNA+L+oKrW912MBunR7fZ44DrgqcDXaZZNeCfwdODTuHyROlZV706yFvgwcNykpmuBlaPLd0hj5BhG81I7G/+TwBbAKVX1055L0gC0S7x9lOaBhSf2XI4GxGU61JW0W5+QqbFKciLwauB7wDE9l6OBSLI3zWzod1bVv/VdjwZrs3YbmhnQX6qqn1XVfwPPpXmYzYEu2aGuJTkFOBdYSTMjejGwF7AG+GSSv+2vOulBOYbR2CXZjObbdcuAc4C/67ciDcjJNA8qPM4PQDROhtGaqYlZA9tO0b7NyHFS55KcQPMVo6uAg6tqXc8laQAmLc9xNfDGnsvRsE0MItZU1XcnN7Sz9Se+KbL3WKvSoCQ5iOZBmp+rqldV1Zqquruqvk3zocgNwKvbh2FL4+YYRvNKG0R/Ajga+BTwohk8eE6atiS7AG8HPlJVF/Zdj4bFMFoz9f12u+sU7bu026nWlJbmVJJXAu8F/osmiL6p55I0HFvTvBYuBe5JUhM/wGntMf/Y7nt3b1VqCCbem/93ivaJsHrLMdSi4Xp2u710tKGq7ga+QTMGeeI4i5JaU45h2g+Xd6Z5cNyacRalYWr73D8DLwD+CXhhVfngQo3LE2iWhTl28vilHcMc2B5zTbvvyP7K1ELkmtGaqYkBxvIki6rqgYmGJA+j+YrReuBrfRSnYWkflHQ68B3gaVV1a88laVjuBT40RduTaAKXr9AMgF3CQ126giZE2SXJ5lV130j7Hu127Vir0tBs0W63n6J9Yv9o/5TG4RLgT4DDaELAyQ4AtgKuqKp7x12YhiXJ5jQzoY8APgYcO3lMLY3BWqYewzwL2IHmWSN34N+OmmOG0ZqRqvpBklXAcuAE4MxJzW+hWRvwg1V1Vx/1aTiSvBF4K/AtYLlLc2jc2uUPXrqhtiRvpgmjP1pVZ42zLg1PVd2a5ByaoOVNwBsm2pI8jeYBhrcDF/VToQbiy8ArgJcl+WBV3TDRkOQZNBMW7gG+2lN9GrZzaZaReUGSM6vqmwBJHgr8VXvM+/sqTsPQPqzwM8AzacLAlxlEa9yq6jtMPYa5jCaMfl1VXTvOujQMhtGajZfTDCRWJDkUWA3sAxxMszzH63usTQOQ5CU0QfT9NIPfE5OMHra2qlaOuTRJ6suraN6LX5/kAJolEXakWav3fpoH1Ey1jIc0F84Fvgg8FVid5LPATTRLGT2b5gFxr62q2/orUQtJ+/Xxia+Q79Bu90uysv3vW6vqNQBVdUeS42j66WVJzgbWAYcDu7X7zxlX7Vo4ptMPgQ/QBNG30qyj/6YNjGEuq6rLOitYC9Y0+6LUC8NozVg7O/rJNGHgYTRvqDcCK4C3OENVY7Bzu90MeOUUx1wOrBxLNZLUs6q6Jck+NLOinwvsC9wJfB74m6py+Sx1qqoeSPJMmm/OvYCmH25FE/hdCKyoqlU9lqiFZ0/gJSP7lrQ/AD8EfhG8VNX5SQ6kmTjzfOChwLU0H+at8OFxmqHp9MOJMcx2NN9kmsplc1WcBmVar4lSH+J7rSRJkiRJkiSpa4v6LkCSJEmSJEmStPAZRkuSJEmSJEmSOmcYLUmSJEmSJEnqnGG0JEmSJEmSJKlzhtGSJEmSJEmSpM4ZRkuSJEmSJEmSOmcYLUmSJEmSJEnqnGG0JEmSJEmSJKlzhtGSJEmSJEmSpM4ZRkuSJEmSJEmSOmcYLUmSJEmSJEnqnGG0JEmSJEmSJKlzhtGSJEmSJEmSpM4ZRkuSJEmSJEmSOmcYLUmSJEmSJEnqnGG0JEmSJEmSJKlzhtGSJEmSJEmSpM79P55e2gx5pMIHAAAAAElFTkSuQmCC" | |
| }, | |
| "execution_count": 2, | |
| "metadata": { | |
| "image/png": { | |
| "height": 372, | |
| "width": 721 | |
| }, | |
| "needs_background": "light" | |
| }, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "import matplotlib.pyplot as plt\n", | |
| "import numpy as np\n", | |
| "import scipy.stats as sts\n", | |
| "\n", | |
| "plot_x = np.linspace(0, 15, 200)\n", | |
| "plot_y = sts.gamma(a=2, scale=1/0.5).pdf(plot_x)\n", | |
| "plt.figure(figsize=(12, 6))\n", | |
| "plt.plot(plot_x, plot_y)\n", | |
| "plt.title(r'Probability density Gamma(x | α=2, β=0.5)')\n", | |
| "plt.show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "These assumptions lead to the following model:\n", | |
| "\n", | |
| "* Likelihood: $\\prod_{i=1}^6 \\text{Binomial}(s_i\\,|\\,\\theta_i, n_i)$, where $s_i$ is the number of successful recoveries, $f_i$ is the number of failures (did not recover), and $n_i=s_i+f_i$ the number of patients. Note that each study has its own $\\theta_i$, whereas Model 1 had the same $\\theta$ for all 6 studies.\n", | |
| "\n", | |
| "* Prior: $\\prod_{i=1}^6 \\text{Beta}(\\theta_i\\,|\\,\\alpha,\\beta)$.\n", | |
| "\n", | |
| "* Hyperprior: $P(\\alpha,\\beta) = \\text{Gamma}(\\alpha\\,|\\,2,0.5)\\,\\text{Gamma}(\\beta\\,|\\,2,0.5)$.\n", | |
| "\n", | |
| "This model has 8 parameters (for each of the treatment and control groups), namely $\\theta_1, \\ldots, \\theta_6$, $\\alpha$, and $\\beta$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Since the posterior does not have a closed-form analytical solution, we have to calculate the marginal likelihood by integrating out all of the parameters in the model.\n", | |
| "\n", | |
| "$$ P(\\text{data}) = \\int_0^{\\infty} \\int_0^{\\infty} \\int_0^1\\cdots\\int_0^1 \\left[\\prod_{i=1}^6 \\text{Binomial}(s_i\\,|\\,\\theta_i,n_i)\\,\\text{Beta}(\\theta_i\\,|\\,\\alpha,\\beta) \\right] P(\\alpha,\\beta)\\ \\text{d}\\theta_6\\cdots\\text{d}\\theta_1\\,\\text{d}\\beta\\,\\text{d}\\alpha$$\n", | |
| "\n", | |
| "This looks like a crazy 8-dimensional integral, but we can actually integrate out the $\\theta_i$ analytically, leaving a 2-dimensional integral over $\\alpha$ and $\\beta$.\n", | |
| "\n", | |
| "First, note that $P(\\alpha,\\beta)$ does not contain $\\theta_i$, so we can move it outside of the $\\theta_i$ integrals.\n", | |
| "\n", | |
| "$$ = \\int_0^{\\infty} \\int_0^{\\infty} P(\\alpha,\\beta) {\\color{blue}{\\int_0^1\\cdots\\int_0^1 \\left[\\prod_{i=1}^6 \\text{Binomial}(s_i\\,|\\,\\theta_i,n_i)\\,\\text{Beta}(\\theta_i\\,|\\,\\alpha,\\beta) \\right] \\ \\text{d}\\theta_6\\cdots\\text{d}\\theta_1}}\\,\\text{d}\\beta\\,\\text{d}\\alpha$$\n", | |
| "\n", | |
| "Next, since there are no factors containing two different $\\theta_i$ variables, we can rearrange the integrals and the products (the blue part) like this:\n", | |
| "\n", | |
| "$$ = \\int_0^{\\infty} \\int_0^{\\infty} P(\\alpha,\\beta) {\\color{blue}{\\left[\\prod_{i=1}^6 \\int_0^1\\text{Binomial}(s_i\\,|\\,\\theta_i,n_i)\\,\\text{Beta}(\\theta_i\\,|\\,\\alpha,\\beta)\\,\\text{d}\\theta_i\\right]}}\\,\\text{d}\\beta\\,\\text{d}\\alpha$$\n", | |
| "\n", | |
| "**Note that we cannot always swap products and integrals.**\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Since the beta distribution is conjugate to the binomial, the blue integrals above can be evaluated analytically (much like we did for Model 1), to get\n", | |
| "\n", | |
| "$$ = \\int_0^{\\infty}\\int_0^{\\infty} P(\\alpha,\\beta) \\left[\\prod_{i=1}^6 \\frac{1}{(s_i+f_i+1)\\,\\text{B}(s_i+1,f_i+1)}\\,\\frac{\\text{B}(\\alpha+s_i, \\beta+f_i)}{\\text{B}(\\alpha,\\beta)}\\right]\\,\\text{d}\\beta\\,\\text{d}\\alpha$$\n", | |
| "\n", | |
| "Finally, move all the factors that do not depend on $\\alpha$ or $\\beta$ out of the integrals.\n", | |
| "\n", | |
| "$$ = \\left[\\prod_{i=1}^6 (s_i+f_i+1)\\,\\text{B}(s_i+1,f_i+1) \\right]^{-1} \\int_0^{\\infty}\\int_0^{\\infty} P(\\alpha,\\beta)\\, \\text{B}(\\alpha,\\beta)^{-6} \\prod_{i=1}^6 \\text{B}(\\alpha+s_i, \\beta+f_i)\\ \\text{d}\\beta\\,\\text{d}\\alpha$$\n", | |
| "\n", | |
| "Unfortunately we cannot evaluate the remaining integrals analytically, so we resort to a numerical calculation." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Task 2\n", | |
| "\n", | |
| "1. Write a Python function to compute the value of the integral above using the `scipy.integrate.dblquad` function.\n", | |
| "2. Write a Python function to compute the marginal likelihood of the hierarchical model for the medical trials data set." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 0, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "from scipy.integrate import dblquad\n", | |
| "\n", | |
| "\n", | |
| "# The integrand for the double integral over a (alpha) and b (beta).\n", | |
| "def integrand(a, b, data):\n", | |
| "\n", | |
| " return dblquad((betaln(a, b)**(-6))*np.prod([a + betaln(data[i,0], b + data[i,1]+1) for i in range(6)]), a, b, 0, inf)\n", | |
| "\n", | |
| "\n", | |
| "def hierarchical_marginal_likelihood(data):\n", | |
| "\n", | |
| "\n", | |
| "\n", | |
| " return ...\n" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.7.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 1 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment