PyMC3 implementation of Kruschke's BEST

BESTpymc3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unexpected differences in the estimates returned by Kruschke's BESTmcmc and this PyMC3 implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As far as I can tell, this PyMC3 model should be doing the same things as Kruschke's BESTmcmc, but they converge to very different estimates. The difference comes from the sd hyperparameter of the Guassian prior placed on the groups' means. Kruschke uses a $\\tau$ precision paramertiztion that is based on the standard deviation of the data ($\\tau = (5\\sigma)^{-2}$ where $\\sigma$ is the sd of the joint data). I can get the PyMC3 model to come to a similar inference if I use a more diffuse prior, but it doesn't make sense to my why this happens. I'd really like to know why inference in JAGS Gibbs sampling and PyMC3 NUTS yields different estimtes of the groups' means.\n",
    "\n",
    "Code heavily informed by https://docs.pymc.io/notebooks/BEST.html and Meredith, M., & Kruschke, J. (2018). Bayesian Estimation Supersedes the t-Test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "\n",
    "def gammaShRaFromMeanSD(mean, sd):\n",
    "    if mean <= 0:\n",
    "        raise RuntimeError(\"mean must be > 0\")\n",
    "    if sd <= 0:\n",
    "        raise RuntimeError(\"sd must be > 0\")\n",
    "    shape = (mean ** 2) / (sd ** 2)\n",
    "    rate = mean / (sd ** 2)\n",
    "    return shape, rate\n",
    "\n",
    "def gammaShRaFromModeSD(mode, sd):\n",
    "    if mode <= 0:\n",
    "        raise RuntimeError(\"mode must be > 0\")\n",
    "    if sd <= 0:\n",
    "        raise RuntimeError(\"sd must be > 0\")\n",
    "    rate = (mode + np.sqrt(mode ** 2 + 4 * (sd ** 2))) / (2 * (sd ** 2))\n",
    "    shape = 1 + mode * rate\n",
    "    return shape, rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1name = 'controls'\n",
    "y1 = np.array([3.56557903,3.85970963, 3.06846955, 3.85970963, 3.06846955, 4.1256442 , 3.77697393, 4.1256442])\n",
    "y2name = 'patients'\n",
    "y2 = np.array([3.77697393, 2.26464194, 2.07558093, 3.21690876, 3.37651802, 3.85970963, 3.56557903, 4.1256442, 2.88408007, 4.1256442, 2.24897313, 2.80253498, 3.37651802, 4.1256442, 3.77697393, 3.0954668 , 4.1256442])\n",
    "y_joint = np.concatenate((y1, y2), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### the observed means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "empirical y1 mean: 3.6813\n",
      "empirical y2 mean: 3.3425\n",
      "empirical mean difference: 0.3387\n"
     ]
    }
   ],
   "source": [
    "print('empirical y1 mean: {:0.4f}'.format(y1.mean()))\n",
    "print('empirical y2 mean: {:0.4f}'.format(y2.mean()))\n",
    "print('empirical mean difference: {:0.4f}'.format(y1.mean()-y2.mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Set hyperparameters\n",
    "muPriorMean1 = muPriorMean2 = y_joint.mean()\n",
    "muPriorSD1 = muPriorSD2 = y_joint.std() * 5\n",
    "# convert SD to precision\n",
    "tauPriorPrecision1 = ( muPriorSD1 ) ** -2\n",
    "tauPriorPrecision2 = ( muPriorSD2 ) ** -2\n",
    "\n",
    "sigmaPriorMode1 = sigmaPriorMode2 = y_joint.std()\n",
    "sigmaPriorSD1 = sigmaPriorSD2 = y_joint.std() * 5\n",
    "\n",
    "nuPriorMean = 29.\n",
    "nuPriorSD = 29.\n",
    "\n",
    "### Prior for group's SD\n",
    "Sh1, Ra1 = gammaShRaFromModeSD(mode=sigmaPriorMode1, sd=sigmaPriorSD1)\n",
    "Sh2, Ra2 = gammaShRaFromModeSD(mode=sigmaPriorMode2, sd=sigmaPriorSD2)\n",
    "\n",
    "### Prior for normality\n",
    "ShNu, RaNu = gammaShRaFromMeanSD(nuPriorMean, nuPriorSD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "### set some names\n",
    "l_g1_mean = r'$\\mu_{' + y1name + '}$'\n",
    "l_g2_mean = r'$\\mu_{' + y2name + '}$'\n",
    "l_g1_sd = r'$\\sigma_{' + y1name + '}$'\n",
    "l_g2_sd = r'$\\sigma_{' + y2name + '}$'\n",
    "l_g1_stat = f'{y1name}'\n",
    "l_g2_stat = f'{y2name}'\n",
    "l_diff_means = r'$\\Delta  \\mu$'\n",
    "l_diff_sd = r'$\\Delta  \\sigma$'\n",
    "l_effect_size = 'effect size'\n",
    "l_nu = r'$\\log_{10} (\\nu)$'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [nu, $\\sigma_{patients}$, $\\sigma_{controls}$, $\\mu_{patients}$, $\\mu_{controls}$]\n",
      "Sampling 2 chains: 100%|██████████| 41000/41000 [01:16<00:00, 533.17draws/s]\n"
     ]
    }
   ],
   "source": [
    "prior_on_mean_sd = tauPriorPrecision1\n",
    "\n",
    "with pm.Model() as model:\n",
    "    group1_mean = pm.Normal(l_g1_mean, mu=muPriorMean1, sd=prior_on_mean_sd)\n",
    "    group2_mean = pm.Normal(l_g2_mean, mu=muPriorMean2, sd=prior_on_mean_sd)\n",
    "\n",
    "    group1_sd = pm.Gamma(l_g1_sd, alpha=Sh1, beta=Ra1)\n",
    "    group2_sd = pm.Gamma(l_g2_sd, alpha=Sh2, beta=Ra2)\n",
    "\n",
    "    nu = pm.Gamma('nu', alpha=ShNu, beta=RaNu)\n",
    "\n",
    "    lambda1 = group1_sd ** -2\n",
    "    lambda2 = group2_sd ** -2\n",
    "\n",
    "    group1 = pm.StudentT(l_g1_stat, nu=nu, mu=group1_mean, lam=lambda1, observed=y1)\n",
    "    group2 = pm.StudentT(l_g2_stat, nu=nu, mu=group2_mean, lam=lambda2, observed=y2)\n",
    "\n",
    "    diff_of_means = pm.Deterministic(l_diff_means, group1_mean - group2_mean)\n",
    "    diff_of_stds = pm.Deterministic(l_diff_sd, group1_sd - group2_sd)\n",
    "    effect_size = pm.Deterministic(l_effect_size, diff_of_means / np.sqrt((group1_sd ** 2 + group2_sd ** 2) / 2))\n",
    "    nu_transformed = pm.Deterministic(l_nu, np.log10(nu))\n",
    "\n",
    "    trace = pm.sample(20000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<matplotlib.axes._subplots.AxesSubplot object at 0x1c2567e550>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c2856bd68>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28f56b70>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28f7ce48>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28fa7198>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c285597f0>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c2894e978>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28283c50>],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAALICAYAAACThQrSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FNX+x/HP2SxVqkCQEBQBaSEhoYnSgoAg99IUpV0LKFcEVFRAvCqgooBy9YI0Qbii0kWKiCiIKOClBBJAUEEhUqX3Iinz+yNhfwkpbMLsbrJ5v54nT3Z3zsyeM7t7vvudc2bWWJYlAAAAAMCNcfi6AgAAAADgD0iuAAAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkCAAAAABuQXAE3yBjTyRhjGWOq+7ouAACkh1gFeAfJFXDjukmKktTV1xUBACADxCrAC0iugBtgjCkiqZmkx5UUuK4+vtoYUy35diljzE8+qiIAII/LJFbVNsb8YIzZaYxJTB7Zes1nFQX8gNPXFQByuY6SVlqWtc0Yc8EYU8eyrC2SqkjanVwmTNJ2n9UQAJDXpYlVknZKmivpEcuyNhpj3pBUUNJwH9YTyPUYuQJuTDdJ85Jvz5PUzRhzm6SDlmUlJj8eJmmbLyoHAIDSiVWSWkraYlnWxuTHt0m62bIsywf1A/wGI1dANhljSklqIOn+5IfmSvpe0lqlTqbqJi8DAMCrMolVJ5V6VkUdSVu8WzvA/zByBWRfZ0nLLMv6S5Isy9or6U9J4UqaWiFjzB2SOohpgQAA38goVjmVNLNCxpiqSkq+5viqkoC/YOQKyL5uksKMMbEpHisl6YKkQ8aYrUoawfpZ0qOS3vB6DQEAeV1GsWqHpCLJF1w6LqmbZVknfFA/wK8YptYC9jLG/CYpwrKsc76uCwAAALyHaYGAjYwxRSUlklgBAADkPYxcAQAAAIANGLkCAAAAABuQXAEAAACADbJytUDmDwIA7GQ8tF3iFQDATm7HK0auAAAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkCriM+MWtXdc5qeQBA7kNsAJAeY1luf9jpFZBnjYo+7nbZgbVLyelw/+d74hOtLJUH/Ai/c4VcLSuxYUhEaQ/WBICHuR2vsvIjwgDc4HQYAi4AAEAexLRAAAAAALAByRXyHOa9AwAAwBOYFog8h2l7AAAA8ARGrgAAADyMqwsCeQMjVwAAAB7GrAkgb2DkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcAQAAAIANSK4AAAAAwAYkV0AWWZalJW+/pHfa19fYh5rp4M9b0y13cOdW/eehpnqnfX0tefslWVbSZXUP/bpdEx9po3FdIzW+R0tt3LhRkvTOO+8oPDxc4eHhqlWrlgICAnTy5EmvtQsA4HknD/6hCY+01pgODTTrxScUH3clTZkLp0+qefPmKlKkiPr3759q2dy5cxUWFqaQkBANHjzY9fgff/yhFi1aKCwsTJGRkTpw4IDH2wIgLZIrIIt+XbdSJ/bt0cDFG9XplX9r0cjB6ZZbNHKQOr38bw1cvFEn9u3Rrh+/lSR9NfZ1tXhyoJ6Zs1otn3rRFRwHDRqkmJgYxcTEaOTIkWrWrJluvvlmr7ULAOB5y8e9rsY9+mjg4o0qVKyEohbNTFMmX4ECeuONNzRmzJhUj584cUKDBg3St99+qx07dujIkSP69tuk2DJw4EA98sgj2rZtm4YOHaqXXnrJK+0BkBrJFXK82NhYVa9eXU888YRq1aqlHj16aOXKlWrUqJHuuOMObdy4URcuXFCvXr1Uv359RUREaPHixa51mzRpojp16qhOnTr68ccfJUl7otZpSu8Omjmop969/y7NebmPa2Tpen5evVwRf+8iY4xuDauny+fO6OyxP1OVOXvsT/114Zxuq11fxhhF/L2Ldn73lSTJSPrr/DlJ0uXz5xQUFJTmOWbPnq1u3bpld5cBQJ7liZixevVqRUZGqnPnzqpevbp69OjhdsxIybIs/b5prWq1aCdJqvP3Ltr53bI05fIXukmNGzdWwYIFUz2+Z88eVa1aVWXKlJEktWzZUgsWLJAk7dy5Uy1atJAkNW/e3NUmAN7FjwgjV/jtt980f/58TZkyRfXr19esWbO0du1aLVmyRG+99ZZq1qype+65R9OnT9fp06fVoEEDtWzZUoGBgVqxYoUKFiyo3bt3q1u3boqKipKUND3vuflrVbTMLZrc82/6I2aDKkY01NIxr2hP1FrXc88p5NTRS/EKa91JkT2f1Zmjh1Wi7P8nRMUDg3T22J8qVuYW12Nnj/2pYoEpy5TTmaOHJUl/H/impvd/SMv+M1xWYqK2bVqfqq0XL17U8uXLNX78eI/sSwDwd56IGdHR0dqxY4eCgoLUqFEjrVu3TrqpepqYcdXVmJHSxdMnVbBIMQU4k75+FS8blObgXGaqVKmiX375RbGxsQoODtaiRYt05UrStMLatWtrwYIFevbZZ7Vw4UKdO3dOJ06cUKlSpbK7GwFkA8kVcoXbb79doaGhkqSQkBC1aNFCxhiFhoYqNjZWBw4c0JIlS1xTKC5fvqx9+/YpKChI/fv3V0xMjAICArRr1y7XNiuE1FHx5CQpqFotnTq0XxUjGurvA0ekeu4hEaU1Kvp4ikfSHq00Mqnup3dE05ikMus/+6/+/sIbqtWinbZ9s0iPP/64Vq5c6Sr3xRdfqFGjRkwJBIBs8kTMaNCggYKDgyVJ4eHhio2NlUKqp4kZmbHSiR8yJu1jGShZsqQmTZqkLl26yOFw6O6779aePXskSWPGjFH//v310UcfqWnTpipfvrycTr7mAd7Gpw65QoECBVy3HQ6H636ijOLj4xUQEKAFCxaoWrVqqdYbPny4ypYtq61btyoxMTHVFIuAfPldt43DocSEeEm67shV8cAgnT5yyLX8zNFDKlqmbKrnLR5YTmePpixz2DWytWXpXLUb9JYkKbRVB4188/lU686ZM+e6UwLjEy05He4HZADISzKKGQ6HI9sxI+U2AwICFB+ffsy46mrMmN73QZ0/eUy/NWmoyn1H6vL5s0qIj1eA06kzRw6pWOmyadbNTLt27dSuXdK0wilTpiggIECSFBQUpM8//1ySdP78eS1YsEDFixfP0rYB3DiSK+RqTofR8csJqtGgiR4b+rbavzhKxhgd+mWbgqqHaeXuP1W8bJDe3npSUYtnKSEh4brbvN7IVY1mrfW/udNUu3Un7d++WQWLFEs1JVCSipW5RfkLF9G+bVGqEFpX0Uvn6q6uTyQtK32L9m7+UZXqNdLvG9fojjvucG3/8rmz+nrVakUM+s81o2WpDYko7fY+AgCk1rp1a73//vt6//33ZYxRdHS0IiIidObMGQUHB8vhcGjGjBnZihnX6jVxvqT/jyWV6jXST99+odqtO2nL0rmqEXlflup+9OhRBQYG6tSpU5o4caLmzZsnSTp+/LhuvvlmORwOjRw5Ur169crSdgHYg+QKfuGe3i9o6ZhXNLZLM8myVKJcBT02bpYaPtRTMwf20vYVS1SpfiPlL1T4hp+rWuNW+nXtSo3p0ED5ChZS5+HjXMvGdY3UkF9/kiR1/Nc7+mzY04r767Kq3n2PqjVqKUm6/9V39cU7LysxIUHOAgW06KMpWpG8/o7vvtQdDSOVv9BNN1xPAED6Xn31VQ0YMEBhYWGyLEsVK1bU0qVL1bdvXz3wwAOaP3++mjdvrptusr8vvu+ZoZr90j/1zYS3FFQ9VPU79pAk7fx+uQ7ujFGrp4ZIkipWrKizZ8/qypUrWrRokb755hvVrFlTzz77rLZuTfoJkKFDh6pq1aqSki668dJLL8kYo6ZNm2rChAm21x3A9ZksXO0m65fFAbwgsxGe9KQ9hyp3lb+6DuAHPDW3lXgFr/B0bACQY7gdr7gUOwAAAADYgOQKAAAAAGxAcgUAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXAAAAAGADkisAAAAAsAHJFQAAAADYgOQKAAAAAGxAcgUAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXAAAAAGADkisgF4pPtDxaHgAAAFnn9HUFgJTiEy05HcbX1cjxnA6jUdHH3S4/JKK0B2sDAAAAieQKOQxJAwAAAHIrpgUCAAAAgA1IrgAAQJ7GeakA7MK0QAAAkKdldUq6xLR0AOlj5AoAAAAAbEByBQAAAAA2ILkCAAAAABuQXAEAAOQw2bnIBhfmAHyPC1oAAADkMFxkA8idGLkCAAAAABuQXAEAAACADUiuAAAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkC8oD4RMuj5QEAACA5fV0BAJ7ndBiNij7udvkhEaU9WBsAAAD/xMgVAAAAANiA5AoAAAAAbEByBY/i3B0AAADkFZxzBY/iXB8AAADkFYxcAQAAAIANSK4AAAAAwAYkVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcAQAAAIANSK4AAAAAwAYkVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcAQAAAIANSK7gtvhEy9dVgJdk57Xm/QEAAPI6p68rgNzD6TAaFX08S+sMiSjtodrAk3itAQAAso6RKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAA/EpevcBOVtudV/cT4Elc0AIAAPiVrF6Ux18uyJNX2w3kJIxcAQAAAIANSK4AAAAAwAYkVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguYJfWTfrA/3nwSZ6r3NjrZ052fX4yslva2TrUI3rGqnw8HD9snaFJCk2ZoPGPtRM4//RSsf37ZEkXTp3RtP7PijLSv/HFaf07qADO2Nc908d2qf/PNhEkrQnap2KFy+ucd2a693779bKD95xPT68aSWN69Zc/+7UUB883k4///CNR/YBAMB962Z9oFq1amUaN8Z1jUwTN+rXr5+luBEVFeW6f23cuBofUsaN1atXEzeAXIgfEYbf+PO3n7Vp4afq+/HXCsiXX//t30XVm7RS6VsrS5Ia9eijpo/005CI0q4fWVz7yST1GDNdpw7t14bPPtLfnn9dq6b+W5GPD5AxJlv1aNKkiRq/8ZGuXLqgcV2bq0aTeyVJFcMb6rFxsyRJh37drk+ef1T5ChSUIu63ofUAgKy6Gjf2bNus93aczTBupHQ1brQufFavf2xP3LgaH1LFjfIBGcaNKnc2vbGGA/AYRq7gN47t3aUKoXWVv1BhBTidur3u3dqxalmm6zicTsVdvqy4y5cU4HTqxP69Onv0sCrVbXTD9clf6CaVr1FbJw7sTbMsqFqoWvR+Qf+bN+2GnwcAkD1X40bhwlmPGxcvXiRuAEiDkSv4jbKVa+jrCW/pwumTylegoH5du1LBNWu7lv9v7jRFL52nXU3uVIVH/6VCxUoostezWjjiBeUrWFAPvTFRy94bplZ9h1z3uea+3Cdp1ElSQlycjCPtcYoLp09q3/bNuqf3C7pw6kSa5UE1wvTDxxNuoMUAgBtxNW6cOHFCVy5dzDBulK9ZW397/vVUcWNb6aK668WxbseNHj166JzySXIvbkjxaZYTN4Ccj+QqD4tPtOR0ZG8KQ04UWKmqmj32tKb37az8hW5SuaohcgQkvcXvfPCxpGBljC58PlZfvjtUnYePU1C1UPX9eLkkae/mH1W0zC2yLGnWi08owJlPbZ9/TVLpNM/V5c3JCq4ZLilp7vyMZ3u4lq1Zs0ZbujWXMQ5F9nxGZStX156odWkrnMHc/Nwqq+8nf3v/Ach9rsaNVq1a6ZQKZBg3VkwcmSZuDIkorSc/XJJu3ChaKjDNc82cOVMrAypKShs3YmPWa9w1cUNnfkpbYT+LG4A/IrnKw5wO4zr3yB1DItImGTlN/Y7/UP2O/5Akff3+CBUrGyRJqQJd7969Na3lfanWsyxLq6a9p+6jpmrx6CFq2edFnTq0Tz/Oniq1fC9Ldbh6ztX1HPpluwJvr5qlbedk/vh+AuD/6nf8hxa8NkCjoo9nGDca3P9wqmRIyjxutO7/cpbqkPLcqsz4W9wA/BHnXMGvnD95TJJ0+vAB7fjuS4W3SbpYxNljf7rKLFy4MOmoYApbvpij6o1bqlCxEoq7fEnG4ZBxOBR3+ZJH6nl41w6t+vBdNXyol0e2DwBwjztxY8eqZWnixowZM4gbANJg5Ap+ZebAnrp45pQcznxq/+JoFSpWQpL01djXdXjXTzIyql+jsv72wluuda5cuqgtS+eq14T5kqTGPfpo5qCeCnDmU9eRU2yr29VpH3GXL6lIydJqN+gtrvgEAD42c2BPff7XWZ2Kd2QYN0oGVVDHl8e41rly6aJmzJihlqOTRpuIGwCuIrmCX3ly+tJ0H+8yYqLrdspLsUtS/kKF1XvKItf92+vcpQHzfsjwOf45dXGq+yWDbtWA+WskSZXqNdKU3h3STI+rVK+Rhv+wx/2GAAC84snpS9PEBSl13LhW/kKF9d1337nWcSdu1IsorZXJ5a+NG5Xqpb3SYGRkJHEDyIWYFggAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXAAAAAGADkis/Ep/IL7cDAAD3ZPV7A98zgOvjUux+xOkwaS4lm5khEaU9WBsgc/GJlpwO47HyAIDM8b0BsB/JFQCfIKgDAAB/w7RAAAAAALAByZWHWJalZ555RlWqVFFYWJi2bNmSbrmXX35ZFSpUUJEiRdIsmzdvnmrWrKmQkBB179491bKzZ8+qfPny6t+/v0fqD3jTr+u+1b87NdQ77etr9X/Hplm+4bOPFBoaqvDwcDVu3Fg7d+5MtXzfvn0qUqSIxowZ43qsYsWKrnXq1avn8TYA/szdmHblyhX985//VNWqVVW9enUtWLBAkjR58uR0P8NXrlxRz549FRoaqtq1a2v16tXpbpdzfXzrah9dpUoVjRo1Ks3yd999VzVr1lRYWJhatGihP/74w7Vs3759uvfee1WjRg3VrFlTsbGxkqQePXqoWrVqqlWrlnr16qW4uDhvNQfwKKYFeshXX32l3bt3a/fu3dqwYYOeeuopbdiwIU25du3aqX///rrjjjtSPb57926NHDlS69atU8mSJXX06NFUy1999VU1a9bMo20AvCExIUFLRg/R4xPnq1jZIE34x72q0ayNylaq5ipTu80DWvjmQEnSkiVL9Pzzz2v58uWu5c8995zuu+++NNv+7rvvVLo00wmBG+VuTHvzzTcVGBioXbt2KTExUSdPnpQkde/eXX369JGU+jM8depUSdL27dt19OhR3Xfffdq0aZMcjtTHfplG7Dsp++gRrcNUv359tW/fXjVr1nSViYiIUFRUlAoXLqxJkyZp8ODBmjt3riTpkUce0csvv6xWrVrp/Pnzrte2R48e+vTTTyUlvT8+/PBDPfXUU95vIGAzRq48ZPHixXrkkUdkjFHDhg11+vRpHT58OE25hg0bqly5cmkenzp1qvr166eSJUtKkgIDA13LNm/erCNHjujee+/1XAMAL9n/0xaVCq6om4Mrypkvv2q37qifV3+VqkzBIkVdty9cuCBj/v/CFosWLVKlSpUUEhLitToDeY27MW369Ol66aWXJEkOh8N1cKNYsWKuMik/wzt37lSLFi0kJcW5EiVKKCoqytPNQRak7KPz58+vrl27avHixanKNG/eXIULF5aU9L3mwIEDkpJe3/j4eLVq1UqSVKRIEVe5tm3byhgjY4waNGjgWgfI7UiuPOTgwYOqUKGC635wcLAOHjzo9vq7du3Srl271KhRIzVs2NB1lD4xMVEvvPCC3nnnHdvrDPjC2WOHVfyW8q77xQKDdOZo2i9tEyZMUOXKlTV48GCNGzdOUtKXtNGjR2vYsGFpyhtjdO+996pu3bqaMmWK5xoA5AHuxLTTp09LSppZUadOHT344IM6cuSIa3l6n+HatWtr8eLFio+P1969e7V582bt37/fCy2Cu67to6/3fWbatGmumQS7du1SiRIldP/99ysiIkKDBg1SQkJCqvJxcXH65JNP1KZNG880APAykisPsay088NTHm2/nri4eO3evVurV6/W7Nmz9cQTT+j06dOaOHGi2rZtmyrIAbmam5+Vfv366ffff9fo0aM1YsQISdKwYcP03HPPpXvO4rp167RlyxZ99dVXmjBhgn744Qf76w7kEe7EtPj4eB04cECNGjXSli1bdNddd2ngwIGu5el9hnv16qXg4GDVq1dPAwYM0N133y2nkzMWcpQsfJ/59NNPFRUVpUGDBklKek+sWbNGY8aM0aZNm7Rnzx599NFHqdbp27evmjZtqiZNmthedcAX6MFsNGHCBNf88fr166c6+nbgwAEFBQW5va0KFYJ1olyI/v3TGUlFVTCokl5euknrlq1WbPR6jfzPeF25dEEJcVf00wWH2jwzlDnmyJWKBQbpzJ//fxT07NFDKlbmlgzLd+3a1TUvf8OGDfrss880ePBgnT59Wg6HQwULFlT//v1dn7fAwEB16tRJGzduVNOmTT3bGMCPZDWmlSpVSoULF1anTp0kSQ8++KCmTZuWZrspP8NOp1Pvvfeea9ndd9+d5hxk+Na1fXRG32dWrlypN998U99//70KFCggKWmUKyIiQpUqVZIkdezYUevXr9fjjz8uSXrttdd07NgxffDBB15oCeAdjFzZqF+/foqJiVFMTIw6duyojz/+WJZlaf369SpevHi651ZlpGPHjvo9aq0k6cKpEzq+73fdXP42dX1zsoYsi9GLX25R2wHDFfG3h9TmmaGeahLgccEhETq+f69OHvxD8XFXtPXrRarRLPX0kOP7fnfd/vLLL11fvtasWaPY2FjFxsZqwIAB+te//qX+/fvrwoULOnfunKSkqYPffPONatWq5b1GAX4gqzHNGKN27dq5rvj37bffui56sHv3ble5lJ/hixcv6sKFC5KkFStWyOl0prpQAnwvZR995coVzZkzR+3bt09VJjo6Wk8++aSWLFmS6hzx+vXr69SpUzp27JgkadWqVa7X98MPP9TXX3+t2bNnp7mACZCbMXLlIW3bttWyZctUpUoVFS5cWP/9739dy8LDwxUTEyNJGjx4sGbNmqWLFy8qODhYTzzxhIYPH67WrVtrxKeL9d4DjWQCAnTfgOG6qcTNvmoO4DEBTqfavzhS0/s9JCsxUfXad1PZytW1YtIola8ZrprN2uh/c6cp5F+9lC9fPpUsWVIzZszIdJtHjhxRp06dZElKiI9X9+7dmc8P3AB3Y9ro0aP18MMPa8CAASpTpoyr3Pjx47Vy5co0n+GjR4+qdevWcjgcKl++vD755BPvNw6ZStlHL3AmTeUMCQnR0KFDVa9ePbVv316DBg3S+fPn9eCDD0qSbr31Vi1ZskSWcWjMmDFq0aKFLMtS3bp11bt3b0lSnz59dNttt+muu+6SJN1///0aOpSDxcj9SK48xBijCRMmpLvsahCSpLfffltvv/12uuv//YU3pBfeyPA56rbvprrtu914ZQEfq964lao3bpXqsVZPDXHdbjfoLQ2JKO26FPMXV6Qvrrksc8EO/RUvJZcppm4ffctUWcAm7sa02267Ld3zG8eOTfv7dVLS79H9+uuv9lQSHnO1j77aD4+KPq7CnZ7RTkk7o4+r5Ttz1PKadUZFH9eQiNJq1aqVtm3blmab8fHxXqk74G2MwwIAAACADUiuvIRflwcAgHgIwL8xLdBL+HV5wPviEy05He7/BEJWywPIuqzGQ4mYCCD3MOn9dkUGONR0g7KaXLlb/qU6ZbJbJSDHGLnlWKbLs/KZyO46fIHzOk9lssSrHC6zzyUxLee52j9ntU8dWLtUlg9YcZALOZTbb0pGrgAAAGA7RimRF3HOVTYxZxzwP1n9XNMPAHwOYC/6YeR2jFwly+owNOdQAf4nq5/rgbVLZWn7THeBPyIewk70w8jt3D7n6rXXXlsuyZM9YpCkQx7cfm7EPkmN/ZEa+yM19kdaOX2fHB82bJjtv+7sgXiV0/djdvlruyT/bRvtyn38tW15rV3uxyvLsnLE3/Dhwy1f1yGn/bFP2B/sD/YH+8T3f/66H/21Xf7cNtqV+/78tW20K+M/zrkCAAAAABvkpOTqNV9XIAdin6TG/kiN/ZEa+yMt9ok9/HU/+mu7JP9tG+3Kffy1bbQrA1n5nSsAAAAAQAZy0sgVAAAAAORaJFcAAAAAYAOSKwAAAACwgc+SK2PMzcaYFcaY3cn/S2ZStpgx5qAxZrw36+ht7uwTY0y4MeZ/xpgdxphtxpguvqirJxlj2hhjfjXG/GaMGZLO8gLGmLnJyzcYYyp6v5be48b+eN4YszP5/fCtMeY2X9TTW663P1KU62yMsYwx9bxZP29zZ38YYx5Kfo/sMMbM8nYdcxtjzIPJ+yoxs/ePMSbWGLPdGBNjjInyZh2zIwvtcuszlpO4+53CGJOQ/HrFGGOWeLue7vLXOOhGux4zxhxL8Ro94Yt6ZpUxZrox5qgx5qcMlhtjzLjkdm8zxtTxdh2zw412RRpjzqR4vYZ6u47ZYYypYIz5zhjzc3Kf+Gw6ZbL/mvnqOvKS3pY0JPn2EEmjMyk7VtIsSeN9ff17X+8TSVUl3ZF8O0jSYUklfF13G/dBgKTfJVWSlF/SVkk1rynTV9Lk5NtdJc31db19vD+aSyqcfPupvL4/kssVlfSDpPWS6vm63j5+f9whKVpSyeT7gb6ud07/k1RDUjVJqzN7/0iKlVTa1/W1s13ufsZy2p+73ykknfd1Xd1oi1/GQTfb9Vhu/K4nqamkOpJ+ymB5W0lfSTKSGkra4Os629SuSElLfV3PbLSrnKQ6ybeLStqVznsx26+ZL6cFdpA0I/n2DEkd0ytkjKkrqaykb7xUL1+67j6xLGuXZVm7k28fknRUUhmv1dDuo0RtAAAgAElEQVTzGkj6zbKsPZZlXZE0R0n7JaWU++kzSS2MMcaLdfSm6+4Py7K+syzrYvLd9ZKCvVxHb3Ln/SFJbyjpy9Zlb1bOB9zZH70lTbAs65QkWZZ11Mt1zHUsy/rZsqxffV0Pu7nZLnc/YzmNW98pcgl/jYO59b11XZZl/SDpZCZFOkj62EqyXlIJY0w579Qu+9xoV65kWdZhy7K2JN8+J+lnSeWvKZbt18yXyVVZy7IOS0mNlBR4bQFjjEPSvyUN8nLdfOW6+yQlY0wDJR39+d0LdfOW8pL2p7h/QGnf8K4ylmXFSzojqZRXaud97uyPlB5X0pEWf3Xd/WGMiZBUwbKspd6smI+48/6oKqmqMWadMWa9MaaN12rn/yxJ3xhjNhtj/unrytgkq31OTuFu/CxojIlK/izk1ATMX+Ogu++tB5KnYX1mjKngnap5XG79XLnjLmPMVmPMV8aYEF9XJquSp9RGSNpwzaJsv2ZOOyqWEWPMSkm3pLPoZTc30VfSMsuy9uf8AzLusWGfXN1OOUmfSHrUsqxEO+qWQ6T3Ql/7Y2zulPEXbrfVGPMPSfUkNfNojXwr0/2RfEDmPSVNLckL3Hl/OJU0NTBSSaOaa4wxtSzLOu3huuVomfXFlmUtdnMzjSzLOmSMCZS0whjzS/KRXp+xoV05tn+1KX7emvyaVZK0yhiz3bKsnHaA0l/joDt1/kLSbMuy/jLG9FHS6Nw9Hq+Z5+XG18sdWyTdZlnWeWNMW0mLlBRvcgVjTBFJCyQNsCzr7LWL01nFrdfMo8mVZVktM1pmjDlijClnWdbh5EQhvakqd0lqYozpK6mIpPzGmPOWZeWKE2zTY8M+kTGmmKQvJb2SPFTpTw5ISnmkKljSoQzKHDDGOCUVlx8OWydzZ3/IGNNSSV8wmlmW9ZeX6uYL19sfRSXVkrQ6+YDMLZKWGGPaW5aV4y84kA3ufl7WW5YVJ2mvMeZXJQW/Td6pYs6UWV+chW0cSv5/1BizUEnTnnyaXNnQLrf6HF+wI36meM32GGNWK+mIdU5Lrvw1Dl63XZZlnUhxd6qk0V6olzfk2M/VjUiZkFiWtcwYM9EYU9qyrOO+rJc7jDH5lJRYzbQs6/N0imT7NfPltMAlkh5Nvv2opDRH1CzL6mFZ1q2WZVWUNFBJcx9zbWLlhuvuE2NMfkkLlbQv5nuxbt6ySdIdxpjbk9vaVUn7JaWU+6mzpFVW8tmHfui6+yN5GtwHktrngfNpMt0flmWdsSyrtGVZFZP7jfVK2i/+mFhJ7n1eFinpoicyxpRW0jTBPV6tpR8yxtxkjCl69bakeyWle0WtXMad91RO5E78LGmMKZB8u7SkRpJ2eq2G7vPXOOhOPEt5Tkt7JZ0L4w+WSHok+Qp0DSWduTqNNTczxtxy9Vy/5FNVHJJOZL6W7yXXeZqkny3LejeDYtl/zdy98oXdf0qaG/ytpN3J/29OfryepA/TKf+YcuEVZOzeJ5L+ISlOUkyKv3Bf193m/dBWSVdu+V1JU1kk6XUlfUmWpIKS5kv6TdJGSZV8XWcf74+Vko6keD8s8XWdfbk/rim7Wn58tUA33x9G0rtK+hK5XVJXX9c5p/9J6qSko5Z/JX+2vk5+PEhJU9WlpCuebU3+23F13+fkP3faldF7Kqf/uRk/707+DGxN/v+4r+udSXv8Mg660a6RyZ+nrZK+k1Td13V2s12zlXT15rjkz9jjkvpI6pO83EiakNzu7bklLrnRrv4pXq/1ku72dZ3dbFdjJU3x26b//+7U1q7XzCRvAAAAAABwA3w5LRAAAAAA/AbJFQAAAADYgOQKAAAAAGxAcgUAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXAAAAAGADkisAAAAAsAHJFQAAAADYgOQKAAAAAGxAcgUAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXwDWMMauNMdWSb5cyxvzk6zoBAHAtY0xtY8wPxpidxphEY4xljHnN1/UC8jKnrysA5EBVJO1Ovh0maXt2N2SMKWlZ1ilbagUAQDJjTEFJcyU9YlnWRmPMG5IKShqejW0RqwCbMHIFpGCMuU3SQcuyEpMfCpO07QY2+d412zfpPOdHN7B9AEDe1FLSFsuyNibf3ybpZsuyrGxs671rHyBeAdnDyBWQWrhSJ1N1Jc01xjglvS3JkvSHpKmS3pRUQNIpSZMlfSJpiaSGlmV1Mca0kVTdGPOKpLbJy2YYYwalWO8tSReMMRUkDZV0RtJyy7JWerylAIDcrJZSz6yoI2lLVuOVpP8qKVYNlPSppM9FvAKyjeQKSK22kqZVyBhzh6QOkl6R9JSkxZZlfZ+8bLCkWZZlRRljPkteb5FlWWONMTOTt3VcSYHqd0lnLcsal856dSRtkVRd0hVJ4yzL2uetxgIAcq0Tku6RJGNMVUn3S7pbWY9XxyV9alnWeGPMfZLmEK+A7GNaIJBauCSHMWarko7M/SzpUSUFlXUpyoVI2m6MyS/popKC1dfJy65OyQiTtDV5mysyWK++pE2WZa2Q9L6k8caY8h5qGwDAf8yWVCT5oktTJHWzLOuEsh6vrsYqiXgF3DBGroDUwiRFWJZ1LuWDxpgOkj4wxpyUNFLSPCUFs4vJ9wdJ2mWMKS3pz+TVjkt6QtKtkkYnP3btei9IGmeMGS0pQNI+SUc91joAgF+wLOu8pHbpLFqkrMWr45KeMMYcl3SHpF+Tt0O8ArLBZO+8R8D/GGOKStpsWVZVX9cFAAAAuQ/JFQAAAADYgHOuAAAAAMAGJFcAAAAAYAOSKwAAAACwQVauFsjJWQAAOxkPbZd4BQCwk9vxipErAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcAT4Wn5j1C5tlZx0AuR/9BQDkbMay3O506Z0BN8QnWnI6snaF6VHRx7NUfkhE6SyVB3IoLsWeDfQXAOB1bserrPzOFQA3OB0mS19++OIDAADgH5gWCAAAAAA2ILkCAAAAABuQXAEAAACADUiuAAAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkCriM+0fJ1FdLIap1yYhsAAAD8jdPXFQByOqfDaFT0cbfLD4kone7jO1d/pRUTR8k4jBwBTv194AhVjGiYptyU3h107vgR5StQUJLUa+J8Fbm5jPZu/lFL//2K/ty9U1XmzNFvlSNTrXf5/Dm998Ddqtm8rToMGe1WnQDAsiw9++yzWrZsmQoXLqyPPvpIderUSVMuMjJShw8fVqFChSRJ33zzjQIDA/XDDz9owIAB2rZtm+bMmaPOnTu71hk8eLC+/PJLJSYmqlWrVho7dqyMMV5rGwB4G8kV4CWVGzRRjWZtZIzR4V07NHvIE3r+8/+lW7bLm5MVXDM81WMlygWr8/D3teaTiemus2LSSN1e927b6w0gd4tPtOR0ZJzQfPXVV9q9e7d2796tDRs2qM9TT2njhg3plp05c6bq1auX6rFbb71VH330kcaMGZPq8R9//FHr1q3Ttm3bJEmNGzfW999/r8jIyBtrEADkYCRX8FuxsbFq06aNGjdurPXr16t27drq2bOnhg0bpqNHj2rmzJkKCQnR008/re3btys+Pl7Dhw9Xhw4dFBsbq4cfflgXLlyQJN35zAjdVruB9kSt08oP3tZNJW7Wkd9/UVCN2uoyYpJbR2ILFC7iun3l0kVJWTt6WzLoVkmSSedL0sGdW3X+xDFVvfseHdgZk6XtAvCuG+2bPujVVVcuX5QktX9x1HX7puuNvi+cNle3N+2k0TEnpAJVdOb0aR0+fFjlypVzqz0VK1aUJDkcqc80MMbo8uXLunLliizLUlxcnMqWLZu9nQYAuQTJFfzab7/9pvnz52vKlCmqX7++Zs2apbVr12rJkiV66623VLNmTd1zzz2aPn26Tp8+rQYNGqhly5YKDAzUihUrVLBgQe3evVvNOzyo/jNXSpIO/bpdz81fq6JlbtHknn/THzEbVDGioZaOeUV7otZqTiGnjl6Kd9UhrHUnRfZ8VpK0Y9WX+nr8CJ0/eVyPjp2VYb0/G/6MHA6HQlq00z1PPJ9p8paYmKgv3xuqh96YqN83/mDTngPgSTfSN/Wa9JnyFSio4/t+15yXnrxu3/Tcc89p7rIVaepwtW86c/SwSpQNcj0eHBysgwcPpptc9ezZUwEBAXrggQf0yiuvZNo33XXXXWrevLnKlSsny7LUv39/1ahRw4a9BwA5F8kV/Nrtt9+u0NBQSVJISIhatGghY4xCQ0MVGxurAwcOaMmSJa7pLJcvX9a+ffsUFBSk/v37KyYmRgEBATq6d5drmxVC6qh48heRoGq1dOrQflWMaKi/DxwhKen8poyOEofc8zeF3PM37d38o1ZMGqUnJi9IU6bLm5NVPLCc/rpwXp8O6qnoL+epzt+7ZNjG9fOmq1qjlipxS/ns7SQAXncjfdPnbzyvw7t+ksPh0PF9e1zbzKhveu+991T2kczOG017wZv0kqaZM2eqfPnyOnfunB544AF98skneuSRRzLc6m+//aaff/5ZBw4ckCS1atVKP/zwg5o2bXrd/QMAuRXJFfxagQIFXLcdDofrvsPhUHx8vAICArRgwQJVq1Yt1XrDhw9X2bJltXXrViUmJqpAwYKuZQH58rtuG4dDiQlJo1TujFxddXvdu3Vy2NO6cOqEpNQXmygemHS0uMBNRRTe5n7t/2lLpsnVvu1Rio1er/Xz/6srly4oIe6KChS+SW2eGerOLgLgAzfSNxUtVUYPzlktKzFRQ+8Kdi3LqG+63shV8cAgnT5yyPX4gQMHFBQUlKZ8+fJJB3CKFi2q7t27a+PGjZkmVwsXLlTDhg1VpEjSlOj77rtP69evJ7kC4NdIrpCntW7dWu+//77ef/99GWMUHR2tiIgInTlzRsHBwXI4HJoxY4YSExKuu63rjVwd37dHpSrcLmOMDv68VQlxV1S4xM2pyiTEx+vyuTO6qWQpJcTF6Zc136jKnc0yfd6ub0523d68ZLYO7IwhsQJyucz6pqKly8rhcCjqizlu9U3XG7mq0ay1/jd3mmq37qT92zerePHiaaYExsfH6/Tp0ypdurTi4uK0dOlStWzZMtPnvfXWWzV16lS99NJLsixL33//vQYMGODeDgCAXIrkCnnaq6++qgEDBigsLEyWZalixYpaunSp+vbtqwceeEDz589X8+bNlb9Q4Rt+rh2rlmrL0nkKcDrlLFBI3UZNdU29Gdc1Us/MWa2EuL80vd9DSoyPV2Jigqrc2VT1Oz0sSdq/I1qfvvCoLp09oyfXrZCzRBk999naG64XgJwns76paduO2r5iiSrVb2RL31StcSv9unalxnRooHwFC2np7I9dy8LDwxUTE6O//vpLrVu3VlxcnBISEtSyZUv17t1bkrRp0yZ16tRJp06d0hdffKFhw4Zpx44d6ty5s1atWqXQ0FAZY9SmTRu1a9fuhusLADmZsSy3f1yUXyFFnpXV37nyZPnsPgeQA3nqB4/8Ol7RXwCA17kdrxzXLwIAAAAAuB6SK+Q58Yl+fVAbAAAAPsI5V8hzrveDmtdiigwAAADcwcgVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkCAAAAABuQXAEAAACADUiugDwgq5ef53L1AAAAWcel2IE8gMvPAwAAeB4jVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcIVfj95gAAACQU/A7V8jVsvr7TRK/4QQAAADPYOQKAAAAAGxAcgUAAAAANiC5AgAAAAAbkFwBAAAAgA1IrgAAAADABiRXAAD4CD8nAQD+hUuxAwDgI1n9OQl+SgIAcjZGrgCkkZ2j6RyBBwAAeR0jVwDS4MeZgbwrPtGS02E8Vh4A/BnJFQAAcGGqIgBkH9MCAQAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBQAAAAA2ILkCAAAAABuQXAEAAACADUiuAAAAAMAGJFfIUeITLV9XAQAAAMgWp68rAKTkdBiNij7udvkhEaU9WBsAAADAfYxcAQAAAIANSK4AAAAAwAYkVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCgAAAABsQHIFAAAAADYguQIAAAAAG5BcAQAAAIANSK4A2CI+0fJoeQAAgJzO6esKAPAPTofRqOjjbpcfElHag7UBAADwPkauAAAAAMAGJFcAAAAAYAOSKwAAAACwAckVAAAAANiA5AoAAAAAbEByBY/ictsAAADIK7gUOzyKy3MDgH+LT7TkdBiPlQeA3ITkCgAAZBsH0QDg/zEtEAAAAABsQHIFAAAAADYguQIAAAAAG5BcAQAAAIANSK4AAAAAwAYkVwAAAABgA5IrAAAAALAByRUAAAAA2IDkCoBPxCdaHi0PAADgbU5fVwBA3uR0GI2KPu52+SERpT1YGwAAgBvHyBUAAAAA2IDkCgAAAABsQHIFAAAAADYguYLbuKAAAAAAkDEuaAG3ZfUCBFLOvwjBulkfaNPCT2VZlup3+oca9+gjSVo5+W1tWviJbipZSnMKORX++Iuq3riVYmM2aPFbgxWQP7+6vvWBSt9aSZfOndHsF59QzwnzZIxJ8xxTendQ2+deU3DNcEnSqUP7NOPZHhowf432RK1T8chHdFO5WxX/118Ka91JLZ8cpD1R6/Tx8w/r5vK3Ke7yJRW5uYyaPvq0ajS916v7B0DOMnbsWP3n/cmZ9lmSdG//l119VtijL+lkQoDbfVZkZKRq9X4lwz7rat90tc8aMnk0fRYAJCO5Qp71528/a9PCT9X3468VkC+//tu/i6o3aaXSt1aWJDXq0UdNH+mnIRGlXUnl2k8mqceY6Tp1aL82fPaR/vb861o19d+KfHxAul9S3NGkSRM1fuMjXbl0QeO6NleNJklfRiqGN9Rj42ZJkg79ul2fPP+o8hUoqCp3NrWh9QBymz9/+1krp069bp+V0tpPJmnhggUatXKrbX3W1b7pap+1efNDqR6X6LMA5F1MC0SedWzvLlUIrav8hQorwOnU7XXv1o5VyzJdx+F0Ku7yZcVdvqQAp1Mn9u/V2aOHValuoxuuT/5CN6l8jdo6cWBvmmVB1ULVovcL+t+8aTf8PAByp2N7d6lhw4ZZ7rMuXbrk0T7r999/T7OMPgtAXsXIFfKsspVr6OsJb+nC6ZPKV6Cgfl27UsE1a7uW/2/uNEUvnaddTe5UhUf/pULFSiiy17NaOOIF5StYUA+9MVHL3humVn2HXPe55r7cR/kKFJQkJcTFyTjSHte4cPqk9m3frHt6v6ALp06kWR5UI0w/fDzhBloMIDcrW7mGln44WuV6ZN5nla9ZW397/nVXn/XPf/5TxxKcHuuzQkJGaOWPv6VZTp8FIC8iuUKeFVipqpo99rSm9+2s/IVuUrmqIXIEJH0k7nzwMd3T+wXJGF34fKy+fHeoOg8fp6Bqoer78XJJ0t7NP6pomVtkWdKsF59QgDOf2j7/mqS055l1eXNymvMXrlqzZo22dGsuYxyK7PmMylaurj1R69JW2OKCIkBeFlipql588UUNvU6ftWLiyFR91vr16zUq+niGfVbRUoFpniuzPis2Zr3GpeizQkJCpHSSq4z6rPhES05H1qYkZmcdAPAFkivkafU7/kP1O/5DkvT1+yNUrGyQJKX6stG7d29Na3lfqvUsy9Kqae+p+6ipWjx6iFr2eVGnDu3Tj7OnSi3fy1Idrp5zdT2HftmuwNurZmnbAPzL448/rmN1OkjKuM9qcP/DqZIhKfM+q3X/l7NUh5TnVmUmoz7LHy+OBABXcc4V8rTzJ49Jkk4fPqAd332p8Db3S5LOHvvTVWbhwoUqW7l6qvW2fDFH1Ru3VKFiJRR3+ZKMwyHjcCju8iWP1PPwrh1a9eG7avhQL49sH0DucPToUUmZ91k7Vi2jzwIAH2HkCnnazIE9dfHMKTmc+dT+xdEqVKyEJOmrsa/r8K6fZGRUv0Zl/e2Ft1zrXLl0UVuWzlWvCfMlSY179NHMQT0V4MynriOn2Fa3q1Nv4i5fUpGSpdVu0FtcdQvI4x544AHtOng00z6rZFAFdXx5jGudixfpswDAW0iukKc9OX1puo93GTHRdTvlpdglKX+hwuo9ZZHr/u117tKAeT9k+Bz/nLo41f2SQbdqwPw1kqRK9RppSu8OaabIVKrXSMN/2ON+QwDkCWvWrEl3Sl3KPutahQtnrc9avXp1que4ts+qVC/tlQbpswAgCdMCAeQK8YlZv6BHdtYBAADILkauAOQKnAQPAAByOkau8jCO6gMAAAD2YeQqD8vqSACjAAAAAEDGGLkCAAAAABuQXAEAAACADUiuAACwCeeyAkDexjlXAADYhHNZPSM+0ZLTYTxWHgDswshVDnT58mU1aNBAtWvXVkhIiIYNG5amzOTJkxUaGqrw8HA1btxYO3fulCRt3LhR4eHhCg8PV+3atbVw4UJvVx/wuri/LmvCw/dqbJdIvde5sVZMGp1uuXnz5qlmzZoKCQlR9+7dXY+/+OKLqlWrlmrVqqW5c+emWe/pp59WkSJFPFZ/eMby5ctVrVo1ValSRaNGjUqz/N1331XNmjUVFhamFi1a6I8//ki1/OzZsypfvrz69+/vemzu3LkKCwtTSEiIBg8e7PE2IMnVpPV6fz3Hz1WZilVUveod6b7m+/btU/PmzRUREaGwsDAtW7ZMkhQbG6tChQq54mefPn1c68yePVuhoaEKCwtTmzZtdPx41n4SAkDeQnKVAxUoUECrVq3S1q1bFRMTo+XLl2v9+vWpynTv3l3bt29XTEyMBg8erOeff16SVKtWLUVFRbnWe/LJJxUfH++LZgBe48xfQE988Lmenbtaz8z+Trv+t0r7tkWlKrN7926NHDlS69at044dO/Sf//xHkvTll19qy5YtiomJ0YYNG/TOO+/o7NmzrvWioqJ0+vRpr7YHNy4hIUH9+vXTV199pZ07d2r27Nmug1BXRUREKCoqStu2bVPnzp3TJEuvvvqqmjVr5rp/4sQJDRo0SN9++6127NihI0eO6Ntvv/VKe3B9iQkJWjJ6iHq+PyfD13zEiBF66KGHFB0drTlz5qhv376uZZUrV1ZMTIxiYmI0efJkSVJ8fLyeffZZfffdd9q2bZvCwsI0fvx4r7YLQO5CcpUDGWNcR8nj4uIUFxcnY1JPbyhWrJjr9oULF1zLCxcuLKczabbn5cuX06wH+CNjjAoUTvrMJMTHKTE+TjIm1fkvU6dOVb9+/VSyZElJUmBgoCRp586datasmZxOpwoUKqzatWtr+fLlSdtKSNCgQYP09ttve7lFuFEbN25UlSpVVKlSJeXPn19du3bV4sWLU5Vp3ry5ChcuLElq2LChDhw44Fq2efNmHTlyRPfee6/rsT179qhq1aoqU6aMJKlly5ZasGCBF1oDd+z/aYtKBVfUzcEVM3zNjTGugydnzpxRUFBQptu0LEuWZenChQuyLEtnz5697joA8jaSqxwqISFB4eHhCgwMVKtWrXTnnXemKTNhwgRVrlxZgwcP1rhx41yPb9iwQSEhIQoNDdXkyZNdyRbgzxITEjSua6TebFlDVe6M1K2hdVNNJfpy00/6ZM1WVQxvoFtD66nn+LkaFX1cP99UUR9+tkSv/7hPp0+e0Hfffaf9+/dLksaPH6/27durXLlyPm4dsurgwYOqUKGC635wcLAOHjyYYflp06bpvvvukyQlJibqhRde0DvvvJOqTJUqVfTLL78oNjZW8fHxWrRokeu9At87e+ywit9S3nU/vdd8+PDh+vTTTxUcHKy2bdvq/fffdy3bu3evIiIi1KxZM61Zs0aSlC9fPk2aNEmhoaEKCgrSzp079fjjj3unQQByJZKrHCogIEAxMTE6cOCANm7cqJ9++ilNmX79+un333/X6NGjNWLECNdR+jvvvFM7duzQpk2bNHLkSF2+fNnb1Qe8zhEQoGfmrNaQ5dt0YMcW/fnbz6mWJ8TH6/j+Peo9ZbG6jvxAn7/xnC6dO6OqdzVXtUYtNblnW3Xr1k133XWXnE6nDh06pPnz5+vpp5/2UYtwIywr7VX7MhrJ//TTTxUVFaVBgwZJkiZOnKi2bdumSs4kqWTJkpo0aZK6dOmiJk2aqGLFihy8yknceM1nz56txx57TAcOHNCyZcv08MMPKzExUeXKldO+ffsUHR2td999V927d9fZs2cVFxenSZMmKTo6WocOHVJYWJhGjhzprRYByIWICjlciRIlFBkZqeXLl6tWrVrplunataueeuopzZhx7VWqyuiElV8DF6xVcM3wNOtxlSr4o0JFi+v2uo2068dV0v+xd9/xVVT5/8dfnyQUEQSkKB0RpIWQUBRWpYgCiwqrIqK4gIirropd8ecqdlH56oplXSyrIgKKusTGCiIoNlpCRFBAZakKQUE6Kef3x71k0zM33Jq8n49HHrl35szcz5ly5n5mzsy98PS84bWPa0zzTl2Jr1KFY5u0oEGL1mRu+JFmHVPoO/Ym+o69ifEp9bnkkkto06YNaWlprFu3jtatWwOwb98+Wrduzbp16yJVNQlA06ZNC1xV2rRpU7HduebNm8eDDz7IwoULqVatGgBffvkln332Gc8++yx79uzh0KFD1KxZk4kTJ3Luuedy7rnnAjBlyhTi4+PDUyEp0zENG7Pr5/9dqSpunb/44ot53X579uzJgQMHyMzMpGHDhnnrv2vXrpx44omsWbMmL0k/8cQTARg2bFixD8oQETlMV66i0Pbt2/NuoN+/fz/z5s2jXbt2BcqsXbs27/X7779PmzZtAPh183/J8T/A4rctG9m+fh11GxU8+5t8HU0AACAASURBVCpS0ez5LZP9u3cBkHVgPz98vZAGLdsUKNOhzx/5YekiAPb+toPMDT9wbJMW5ObksHfnrwBkZGSQkZFB//79Ofvss/n5559Zv34969evp0aNGkqsYkj37t1Zu3YtP/30E4cOHWLGjBkMHjy4QJm0tDSuvPJKUlNT8+7BA5g2bRobNmxg/fr1TJo0iZEjR+Z9od62bRsAv/32G88++yxjx44NX6WkVE07ppC58Sd+3fxf9h04WOw6b968ed5DSFavXs2BAwdo0KABW3/ZRk5ODuC7t27t2rW0atWKJk2asGrVKrZv3w7A3Llzad++fXgrJiIxRVeuotDWrVsZNWoUOTk55ObmMmzYMM455xzuvvtuunXrxuDBg3n66aeZN28eVapUoW7durzyyisArE/7moUvTyY+IQGLi2PIHY9ydN16Ea6RSGjt3v4Lb064FpeTi3O5dDprCO179efuu+9m/bEn0aH3QE76wxms/WoBT1xwKhYfzx9vuIej6xxL1sEDTLncdyXi3w3r8tprr6mrVwWQkJDA008/zYABA8jJyWHMmDF07NixQDt66623smfPHi688ELA98U7NTW11Plef/31rFixAoC7776bk046KeR1EW/iExIYfPvDvHTNMN5KgNYDL+LdQ8cx7i+30KRDMh16D6Td2L9xz/03csdDkzCDgX97kkfSd9Dmx8+4++67SUhIID4+nueee45jjz0WgAkTJtCrVy+qVKlCixYtePnllyNbURGJavoGEYWSkpJYsmx5kR9AvO+++/JeP/nkk8VO2+WcYXQ5Z1hI4xOJNo1O6si46Z8UGX7ffffldZU1M865+X64+f4CZapUq86Nb30OlN5Vds+ePUGMWMJh0KBBDBo0qMCw/O3ovHnzypzH6NGjGT16dN776dOnBy0+Cb52p51Fu9POYnxK/bx9/6yrx+eNP65VW6761wdFprvgggu44IILip3nVVddVeB3r0RESqPkKkodfspZIHQPlciRyc51RU5qBLO8iIiIVGxKrkRE/AI9qaETGiIiIpKfHmghIiIilVr+HxwP5TQiUvFZcb8FUgK1IkegPN2HytMtsPA0d3RpENA8RGLZw8u3F3hf3D5RmvKUlyMSqj6VYTtelfTbWRIZh9uAUO/7h6cRkUrDc2OvboFhou5GIhWP7tESqby0/4tIcZRclZMaSRHRSRORyivQ/f+WzoH9LIq+Z4jEJiVX5aQvVSISqPJ8WdIXLJGKIdTJGKi9EIkGMZFcBdpYZOU6qgTYuJRnmlhQ+B6U/MLVJz3U97zEekwVoQ7RGlO0Kc9PLAT6BSvQtkxfxkrmZdkUd99ytO0H0bZvhqMOFUE0thehLl+eadSGSbTx/ECLe++9dw4Qi5dfGgNbIh1EkFSUulSUekDFqYvqEX0qSl1Kq0fmhAkTBgb7A2P4eBVKFWV7igQtu/LTsis/LbvyC8Wy8368cs5V6L977rnHRToG1aVi1qMi1UX1iL6/ilKXilKPWP/TetCy07KLrT8tu9hddvqdKxERERERkSCoDMnVvZEOIIgqSl0qSj2g4tRF9Yg+FaUuFaUesU7rofy07MpPy678tOzKL6LLLpAfERYREREREZESVIYrVyIiIiIiIiGn5EpERERERCQIKkRyZWbVzWyxma0ws2/NrEhfSzO7ycxWmVmGmX1sZi0iEWtZPNblKjP7xszSzWyRmXWIRKyl8VKPfGWHmpkzs27hjNELj+tjtJlt96+PdDMbG4lYy+J1nZjZMP++8q2ZvR7uOMvicZ08kW99rDGznZGItSwe69LczD4xszR/+zUoErGWxmM9Wvjb3gwzW2BmTSMRa0VmZgPN7HszW2dm44sZHxNtVSSY2Utmts3MVpYw3sxssn/ZZphZl3DHGK08LLs+ZrYr33Z3d7hjjFZm1szfvq/2t53XF1NG214xPC67yGx7kX5cYjD+AANq+l9XAb4GehQq0xeo4X99NTAz0nEfQV2Oyfd6MDAn0nGXpx7+cbWAT4GvgG6Rjruc62M08HSkYw1SXdoAaUBd//uGkY67vNtWvvLXAS9FOu4jWCdTgKv9rzsA6yMddznr8SYwyv/6DGBqpOOuSH9APPAD0AqoCqwAOhQqExNtVYSWXy+gC7CyhPGDgA/923oP4OtIxxwtfx6WXR/gvUjHGY1/QCOgi/91LWBNMfuttr3yL7uIbHsV4sqV89njf1vF/+cKlfnEObfP//YrICrPmnqsy+/53h5deHw08FIPv/uBR4ED4YotEAHUI+p5rMsVwDPOud/802wLY4ielGOdXAxMD3lg5eCxLg44xv+6NlH4o5Ie69EB+Nj/+hNgSJjCqyxOBtY55350zh0CZqBl7Jlz7lPg11KKDAFe9W/rXwF1zKxReKKLbh6WnZTAObfVObfc/3o3sBpoUqiYtr1ieFx2EVEhkisAM4s3s3RgGzDXOfd1KcUvx3cWICp5qYuZXWNmP+BLTMaFO0YvyqqHmaUAzZxz70UkQI88blsX+C/XzzKzZmEO0TMPdTkJOMnMPjezr8zM26+Rh5nX/d183X9PAOaHM75AeKjLPcClZrYJ+ADflbio46EeK4AL/K/PA2qZWb1wxljBNQE25nu/ieK/aMREWxWFvC5fKV5Pf7fhD82sY6SDiUZm1hJIwXflPz9te2UoZdlBBLa9CpNcOedynHPJ+K5InWxmicWVM7NLgW7AY+GMLxBe6uKce8Y5dyJwO/C3cMfoRWn1MLM44Ang5kjF55WH9fEu0NI5lwTMA14Jd4xeeahLAr6ugX3wXfF5wczqhDfKsnnd34HhwCznXE74oguMh7pcDLzsnGuKr3vIVP/+E1U81OMWoLeZpQG9gc1AdpjDrMismGGFrx7GTFsVhbwsXynecqCFc64z8BTw7wjHE3XMrCbwFnBDod5JoG2vVGUsu4hse1F3gD5SzrmdwAKgyBl3MzsTuBMY7Jw7GObQAlZaXfKZAfwpLAGVUwn1qAUkAgvMbD2+fsSpFoUPtTispPXhnNuRb3t6Huga5tACVsq2tQmY7ZzLcs79BHyPL9mKSh72keFEaZfAwkqpy+XAG/4yXwLVgfphDS4ApewnW5xz5zvnUvC1wzjndoU/wgprE5D/SlRTCnUhjcW2KoqUuXyleM653w93G3bOfQBUMbOobcPCzcyq4EsOpjnn3i6miLa9EpS17CK17VWI5MrMGhw+u25mRwFnAt8VKpMC/BNfYhV195Ec5rEu+b/sng2sDV+E3pRVD+fcLudcfedcS+dcS3z3wQ12zi2NSMAl8Lg+8vd9Hoyv32/U8VIXfGd1+vrL1MfXTfDHcMZZFo/1wMzaAnWBL8MboXce67IB6Ocv0x5fcrU9nHGWxeN+Uj/fFbc7gJfCG2WFtwRoY2YnmFlVfCcWUvMXiJW2KkqlAiP9T27rAexyzm2NdFCxwMyONzPzvz4Z33fPHZGNKjr4l8uLwGrn3OMlFNO2Vwwvyy5S215CqD8gTBoBr5hZPL4F94Zz7j0zuw9Y6pxLxdcNsCbwpn85b3DODY5YxCXzUpdr/VfhsoDfgFGRC7dEXuoRC7zUY5yZDcbXxelXfE/kikZe6vIfoL+ZrQJygFudc9F2EPS6bV0MzHDORXP3CS91uRl43sxuxNcVZHQU1slLPfoAD5uZw/eE0GsiFm0F5JzLNrNr8e3D8fiekPltjLZVYWdm0/Fto/XNd3/jBHwPZsE59xy++x0HAeuAfcBlkYk0+nhYdkOBq80sG9gPDI/CNixSTgX+DHxjvntWAf4f0By07ZXBy7KLyLZn2r5FRERERESOXIXoFigiIiIiIhJpSq5ERERERESCQMmViIiIiIhIECi5EhERERERCQIlVyIiIiIiIkGg5EpERERERCQIlFyJiIiIiIgEgZIrERERERGRIFByJSIiIiIiEgRKrkRERERERIJAyZWIiIiIiEgQKLkSEREREREJAiVXIiIiIiIiQaDkSkREREREJAiUXImIiIiIiASBkisREREREZEgUHIlUoiZLTCztv7X9cxsZaRjEhERyc/MOpvZp2a2ysxyzcyZ2b2RjkukskuIdAAiUag1sNb/Ogn4xuuEZlbXOfdbIB9WnmlERKTyMrPqwExgpHNusZndD1QH7vE4vY5VIiGiK1ci+ZhZC2Czcy7XPygJyAhgFk/453N/oNOIiIh4dCaw3Dm32P8+AzjWOec8Tq9jlUiI6MqVSEHJFEymugIzzawhMAf4D9AWGArcBRwL7HTOTTCzgUA7M/sbkGBmVYB7gRpAVefcX82sKTAVSAV6AP/yT3MLvrOQdwO7gDnOuXmhr66IiMSgRAr2qugCLNexSiTydOVKpKDO+LpWYGZtgCH4DmDdgenOuTuAbUAzoAqwE9+BByATeA1YBqQDfwGO8pepmW/+/3bOPQFkH57GOTcJaAccAibrYCUiIqXYga9nBWZ2EnA+MAMdq0QiTleuRApKBvab2Qp8V7BWA6OAeOALf5nawATgeqABvoMX+A50K4DTgLeBW4FrnHMH882/s38cgMs3Dc65uWa2EXjazK52zm0OSQ1FRCTWTQcG+x+4lAlc7JzbYWbd0bFKJKKUXIkUlASkOOd25x9oZtOBemZ2HvACviTsFqAekOYvlgmMBZoCjwCzgZf9B6H5zrk5+B6WscbM6gM/H57GzDKB0fiSuA34zjiKiIgU4ZzbA5xbzKi26FglElHm/d5HkYrNzGoBy5xzJxUzbqpz7s8RCEtERMQTHatEIk/JlYiIiIiISBDogRYiIiIiIiJBoORKREREREQkCAJ5oIX6D4qISDBZiOar45WIiAST5+OVrlyJiIiIiIgEgZIrERERERGRIFByJSIiIiIiEgRKrkRERERERIJAyZWIiIiIiEgQKLkSKUN2bmAPHgu0vIhIZVKeNlLtqojECnPOc4Ollk0qrYlpmZ7Ljk+pH8JIRCoUPYq9AsjOdSTEBbYqA2lTQe2qiESc50YukN+5EhERESkgIc50AkpExE/dAkVERERERIJAyZWIiIiIiEgQKLkSEREREREJAiVXUunoqVMiIiIiEgp6oIVUOrr5WkRERERCQVeuREREREREgkDJlYiIiIiISBAouRIREREREQkCJVciIiIiIiJBoORKREREREQkCJRciQRZoI9616PhRURERCoGJVcipVi14EOSkpKYPLwPT484k/VpX5Va/tUbLiU5qRMT0zKZmJZJr5HX0vCENjRq05GOZ5zNhE9/YGJaJhc9+ByN2ybSuG0i3bqkEBcXR3p6ephqJSISO5xzjBs3jtatW5OUlMTy5cuLlNm9ezfJycl5f/Xr1+eGG24A4NNPP6VLly4kJCQwa9asAtPddtttdOzYkfbt2zNu3Dic08kuETkySq4kpoX6qs+JJ5/OihUrGDdjARdMeJK377+xxLIrP36PqjWOLjCsdY/eXP/GZ1z/xkLqNz+RBS89CUDKoKGMm7GAcTMWMHXqVFq2bElycnJI6yIiEou+/3wea9euZe3atUyZMoWrr766SJlatWqRnp6e99eiRQvOP/98AJo3b87LL7/MJZdcUmCaL774gs8//5yMjAxWrlzJkiVLWLhwYVjqJCIVl5IriSrr16+nXbt2jB07lsTEREaMGMG8efM49dRTadOmDYsXL2bv3r2MGTOG7t27071rF0Y+MZWJaZnc/v5yTkjpQZP2STRpn8TVL3/AxLRM/vL8bFp1O5VOZ55LwxPaMGLECM9nJ6vVqImZAXBo/z7Aii13cN8eFk37B33H3lRg+Ek9+xKf4Put7uadurJr25Yi006fPp2LL744gKUkIhI6gbbDKSkprFrwIQC/bdnAP8ecw1OXnMFTl5zBf1csBuDHpZ8z5YohTLv1Mtq1a8eMO6/y3A6vXjCHkSNHYmb06NGDnTt3snXr1hLLr127lm3btnH66acD0LJlS5KSkoiLK/iVx8w4cOAAhw4d4uDBg2RlZXHccceVZ5GJiORJiHQAIoWtW7eON998kylTptC9e3def/11Fi1aRGpqKg899BAdOnTgjDPO4KWXXmLnzp207tyV1qf04ui69Rnzj1lUqVadzA0/MOOOK7l22jwAtnz/DTe+uYhaDY7n3WuHUCv9a1qm9OC9SX/jx6WLisSQNOA8+lx2PQDvvPMOj990G3t+zWTUk68XG/PcZydy+qV/pWr1o0qs19LZr5PU/09Fhs+cOZPZs2eXZ1GJiIREKNvhB/sn0jL5FP7rsR3etW0rzZo1yxvetGlTNm/eTKNGjYqNffr06Vx00UV5J8ZK0rNnT/r27UujRo1wznHttdfSvn37I1hqIiJKriQKnXDCCXTq1AmAjh070q9fP8yMTp06sX79ejZt2kRqaiqTJk0CIPvQQXZu3cwxDY5n9iPj2bpmJXFxcWRu+DFvns06dqH2cY0BSE5O5tctG2mZ0oNzbnmgzHjOO+88vm95Oj8t+4K5/5jI2OfeKjB+y/ffsGPjT5xzywP8tmVDsfP45IXHiUtIIHnQ0ALDN3yzjBo1apCYmOh9AYmIhFgo2+G4uDgat03kN8/tsCvSBby0xGnGjBn865VXy6zjunXrWL16NZs2bQLgrLPO4tNPP6VXr15lTisiUhIlVxJ1qlWrlvc6Li4u731cXBzZ2dnEx8fz1ltv0bZtWwAmpmUCMO+5R6lVrwEXzliAy83l7p5N8+YTX6Xq/17Hx5Obkw3g6crVYSd0/QO/TriOvb/t4Oi69fKGb8hYyubVK3jk7C7k5mSz/7dMplwxhL8877satezdGaz+bC5jn3uryBeCjP+8oy6BIhJ1Qt0OW1yc53a4dsPGbN28Ke8zMn74L2/tqMZc//v8tq5Zyba9Bzmle7cy6/jOO+/Qo0cPatasCcAf//hHvvrqKyVXInJElFxJzBkwYABPPfUUTz31FGbGlu8yaNwuiQN7fs87K7r03Rnk5uSUOa+yzphmbvgRl+xLpDavXkFO1iFq1Dm2QJkeF15GjwsvA3z3G7x/+0gu9SdW33/+MZ++/BRXvDCbqkfVKDBdbm4u38xL5dX7i36pEBGJZuFsh9v3HsCrr75KrwfPZOM3y6he8xiOaXB8sWVXzHmbzgPO91SH5s2b8/zzz3PHHXfgnGPhwoV5TxgUESkvPdBCYs5dd91FVlYWSUlJJCYm8tGzEwHoMewylr87k2dHDiRzww9Fkpny+Hb+eyQmJjJ5eB9mTxzPxROfz7v6NHl4nzKnT31kPAf37eGlq4cyeXgf3nnwlrxx65d/Se2GjWnVqtURxykiEk7hbIfbnnYWrVq1YtKQk3n7gZsYcsejeeMKt8MZc1PpPLBgcrVkyRKaNm3Km2++yZVXXknHjh0BGDp0KCeeeCKdOnWic+fOdO7cmXPPPfeI4xWRys0C+E0H/fiDRKWJxXQNKc34lPoBTROO8iKVVOlPHCg/Ha/CLJRtZHmmUbsqIkHm+XilK1ciIiIiIiJBoORKREREREQkCJRciYiIiIiIBIGSKxERERERkSBQciUiIiIiIhIESq5ERERERESCQMmViIiIiIhIECi5EhERERERCQIlVyIiIiIiIkGg5EokwrJzXVimERHxQu2LiEj5JUQ6AJH8snMdCXEW6TDCKiHOmJiWGdA041PqhygaEansAm2T1B6JiPyPkiuJKjqoi4iIiEisUrdAERERERGRIFByJSIiIiIiEgRKrkRERERERIJAyZWIiIiIiEgQKLkSEREREREJAiVXIiIiIiIiQaDkSkREREREJAiUXImIiIiIiASBkisRERGpULJzXUjLi4iUJCHSAYiIiIgEU0KcMTEt03P58Sn1QxiNiFQmunIlIiIiIiISBEquREREREREgkDJlYiIiIiISBAouRIREREREQkCJVciIiIiIiJBoORKJAbpMcMiIiIi0UePYheJQXrMsIiIiEj00ZUrERERERGRIFByJSIiIiIiEgRKrkRERERERIJAyZWElB6kICIiIiKVhR5oISGlBy+IiIiISGWhK1ciIiIiIiJBoORKREREREQkCJRciYiIiIiIBIGSKxERERERkSBQciUiIiIiIhIESq5ERERERESCQMmViIiIiIhIECi5EhERERERCQIlVyIiIiIiIkGg5EpERERERCQIlFyJiIiIiIgEgZIrERERqdSyc11YphGRii8h0gGIiIiIRFJCnDExLTOgacan1A9RNCISy3TlSkREREREJAiUXImIiIiIiASBkiuRSiDQewN0L4GIiIhI4HTPlUglEOj9BLqXQERERCRwunIlIiIiIiISBEquREREREREgkDJlXim+3BEREREREqme67EM/0OiIiIiIhIyXTlSkREREREJAiUXImIiIiIiASBkisREZEKTPfLioiEj+65EhERqcD0O3ciIuGjK1ciIiIiIiJBoORKREREREQkCJRciYiIiIiIBIGSKxERERERkSBQciUiIiIiIhIESq5ERERERESCQMmViIiISIAC/f0w/d6YSOWg37kSERERCZB+P0xEiqMrVyIiIiIiIkGg5EpERERERCQIlFyJiIiIiIgEgZIrESmiPDde62ZtERERqez0QAuJGZ+//k+WvPMazjm6n3cpp424CoB5zz3KknemcnTdegD0v/ZO2p12FuvTv2b2Q7cRX7Uqwx/6J/Wbt2L/7l0MGDCCPg+/hpkV+YwpVwxh0I330rRDMgC/bdlA4p/7cOnUBfy49HNevenPHNukBdkHD5I04DzOvPLWAsOzDuzng+aNOeGCq2jfq3/4Fk6QBXqjNuhmbZHKzmsbnfTEo9Do5BLb6Om3j+WyZ97w3Ea/cv0IbnjzswJt8SuWTZPeg4tto2se24Beo66L6TZaRKKXkqtKLDvXkRBX9OAVjX5et5ol77zGX1/9D/FVqvKvay+i3elngf8L/akjrqLXyGsKTLNo6j8YMeklftuyka9nvczZN93H/Of/jwf+3//jy2IO2l60TO7B6Mmvc2j/XiYP70v70/sXGA4w0DbRd9BgqlSrTutTeh1BrUVEYkNJbXT95icCBdvoQSn1yUjLLLGN7nP5DcUmVl4cbouvO+koWrTvVGwbveX7b5h60yi10SISEkquKrFYeozs9p/W0KxTV6oeVQOAE7r+gW/nfwBDTilxmriEBLIOHCDrwH7iExLYsfEnft+2ld69e/NlgFdlCqt61NE0ad+ZHZt+ouaxDQqMS05Opt8VN/PlGy/qwC0ilUJJbXTv0deVOE1JbXSrrqcecTxHH11yG924bSe10SISMkquJCYcd2J7/vPMQ+zd+StVqlXn+0XzaNqhc974L2e+SNp7b9CkQ2fOvuk+jjqmDn3GXM87D9xMlerVGXb/s3zwxATO+uv4Mj9r5p1XUaVadQBysrJoUKNKkTJ7d/7Khm+WccYVN7P3tx1Fxjdun8Snrz5zBDUWEYkdgbTRV77saxuD1UZbXNHbx3fs2KE2WkQiQsmVxISGrU6i9+jreOmvQ6l61NE0OqkjcfG+zfeUC0dzxhU3gxlzn32Y9x+/m6H3TKZx20789dU5APy07AtqNTge5+Ciiy5ize5cBt10L7XqNSzyWRc9+FyB/vzv3z4yb9z69K+YfHFfzOLoc9k4jjuxHT8u/bxowE4PdxCRyiOQNvrmm2/mpOseLbGNfv32scQnVPHcRr9y/Yi8cYfb6DePrqo2WkQiQk8LlJjR/U+Xct3r87nyxXepcUwd6jVvBUCteg2Ji48nLi6Ok8//M5u+TSswnXOO+S8+Qb8rbubjKY9x7733kjxoKF9Mfz7gGFom92Dc9E+47vWPOWXo6BLLbfnuGxqecFLA8xcRiVVe2+jFixcXmK5wG33mVbcfcRu9bNkytdEiEhFKriRm7Pl1OwA7t27i20/eJ3ng+QD8vv3nvDLfzv+A405sV2C65e/OoN1pZ3LUMXXIOrCfuLg4LC6OrAP7QxJnRkYG8194nB7DxoRk/iIi0chrG52YmFhgusJttIW4jd665lu10SISMuoWKDFj2i2XsW/Xb8QlVGHw7Y9w1DF1APjwyfvYumYlhlG3cTP+dOekvGkO7d/H8vdmMuaZNwE4bcRVXHDBBfyaHcfwh6cELbbDXVGyDuzn/WaNOPfWh3SjtIhUKl7b6I9n/ItX/PlWcW30tFsvIz6hSsja6Jp166uNFpGQUXIlMePKl94rdvhFDzxb4jRVj6rBFVP+nff+hC49+ec335T4lMS/PD+7wPu6jZuzcuVKJqZl0qrbqbTqVvQpVq26nco9n/6Y9358Sv2AfyNKRCTWeW2jGzWqDz/72sji2ugb3vi0xM8oro2+4c3PADy30SIioaRugSIiIiIiIkGg5EpEREQkxLJzA3tCYaDlRSQ6qFugiIiISIglxFlAXcbHp9QPYTQiEiq6ciUiQaGzsiIiIlLZ6cqViASFzsqKiIhIZacrVxWIrgSIiIiIiESOrlxVILpyICIiIiISObpyJSIiIiIiEgRKrkLgwIEDnHzyyXTu3JmOHTsyYcKEImUef/xxOnToQFJSEv369eO///1vgfG///47TZo04dprr80b1qdPH9q2bUtycjLJycls27Yt5HUR8erAgQM88+f+PHlRH54Yehpz//FIkTLLUqfzwBntmDy8D8nJybzwwgt5426//XYSExNJTExk5syZecPnz59Ply5dSExMZNSoUWRnZ4elPiJezZkzh7Zt29K6dWsmTpxYYrlZs2ZhZixdujRvWEZGBj179qRjx4506tSJAwcOAGrvBcaMGUPDhg1JTEwsdrxzjnHjxtG6dWuSkpJYvnx53riBAwdSp04dzjnnnHCFKyJ+Sq5CoFq1asyfP58VK1aQnp7OnDlz+OqrrwqUSUlJYenSpWRkZDB06FBuu+22AuPvuusuevfuXWTe06ZNIz09nfT0dBo2bBjS/TFJOwAAIABJREFUeogEolq1aoz959tcP3MB46Z/wpov57MhY2mRcp36D2HcjAWkp6czduxYAN5//32WL19Oeno6X3/9NY899hi///47ubm5jBo1ihkzZrBy5UpatGjBK6+8Eu6qiZQoJyeHa665hg8//JBVq1Yxffp0Vq1aVaTc7t27mTx5MqecckresOzsbC699FKee+45vv32WxYsWECVKlXyxhfX3uve2spj9OjRzJkzp8TxH374IWvXrmXt2rVMmTKFq6++Om/crbfeytSpU8MRpogUonuuQsDMqFmzJgBZWVlkZWVhZgXK9O3bN+91jx49eO211/LeL1u2jF9++YWBAwcWOMMpEs3MjGo1fNt9TnYWudlZUGi7L8mqVavo3bs3CQkJJCQk0LlzZ+bMmUPfvn2pVq0aJ510EgBnnXUWDz/8MJdffnnI6iESiMWLF9O6dWtatWoFwPDhw5k9ezYdOnQoUO6uu+7itttuY9KkSXnDPvroI5KSkujcuTMA9erVK/PzAr23FnR/baz6w2mns2nDf0scP3v2bEaOHImZ0aNHD3bu3MnGzVto1qQx/fr1Y8GCBeELVkTy6MpViOTk5JCcnEzDhg0566yzCpytLOzFF1/kj3/8IwC5ubncfPPNPPbYY8WWveyyy0hOTub+++/HOZ3BlOiSm5PD5OF9ePDM9rQ+pQ/NO3UtUubb+e/x5LDeDB06lI0bNwLQuXNnPvzwQ/bt20dmZiaffPIJGzdupH79+mRlZeWdZJg1a1beNCLRYPPmzTRr1izvfdOmTdm8eXOBMmlpaWzcuLFIF601a9ZgZgwYMIAuXbrw6KOPFhiv9r5yS4gz/vHtr2QeyGFiWmaRv4WrfmLh/lp573NqH8cvW7dEOmyRSk9XrkIkPj6e9PR0du7cyXnnncfKlSuL7Tf92muvsXTpUhYuXAjAs88+y6BBgwocrA+bNm0aTZo0Yffu3VxwwQVMnTqVkSNHhrwuIl7FxcczbsYC9u/exWs3j+Lndas5vnX7vPHteg2g88DzSahajTpfz2LUqFHMnz+f/v37s2TJEv7whz/QoEEDevbsSUJCAmbGjBkzuPHGGzl48CD9+/cnIUHNlkSP4pKe/D0VcnNzufHGG3n55ZeLlMvOzmbRokUsWbKEGjVq0K9fP7p27Uq/fv3U3osHpW97IhIZunIVYnXq1KFPnz7F9pueN28eDz74IKmpqVSrVg2AL7/8kqeffpqWLVtyyy238OqrrzJ+/HgAmjRpAkCtWrW45JJLWLx4cfgqIhKAo2rV5oSup7Lmi/kFhh9d51gSqvq29SuuuIJly5bljbvzzjtJT09n7ty5OOdo06YNAD179uSzzz5j8eLF9OrVK2+4SDRo2rRpgaupmzZtonHjxnnvd+/ezcqVK+nTpw8tW7bkq6++YvDgwSxdupSmTZvSu3dv6tevT40aNRg0aFDeQwnU3ktZajdszM5f/nelate2LQW2PRGJDCVXIbB9+3Z27twJwP79+5k3bx7t2rUrUCYtLY0rr7yS1NTUAg+mmDZtGhs2bGD9+vVMmjSJkSNHMnHiRLKzs8nM9PWzz8rK4r333ivxCUIikbB9+3b2794FQNaB/fzw9UIatCyYCP2+/ee81+/8ezbt2/uuauXk5LBjxw7A9/S0jIwM+vfvD5D3lLS9+w/wyCOPcNVVV4W8LiJede/enbVr1/LTTz9x6NAhZsyYweDBg/PG165dm8zMTNavX8/69evp0aMHqampdOvWjQEDBpCRkcG+ffvIzs5m4cKFdOjQQe29eNK+9wDS3puJc44NGUupXvMYGjVqFOmwRCo99a8Jga1btzJq1ChycnLIzc1l2LBhnHPOOdx9991069aNwYMHc+utt7Jnzx4uvPBCAJo3b05qamqJ8zx48CADBgwgKyuLnJwczjzzTK644opwVUmkTFu3buX5v4zA5eTiXC6dzhpC+179mfuPiTTpkEyH3gP5YsbzrF74H+LiE2jXpAF/uO0JJqZlknXwAE9f0g+AakfX4k9/e4pJ3/hOUHzwxH1899lH1K1qXH311ZxxxhmRrKZIAQkJCTz99NMMGDCAnJwcxowZQ8eOHQu09yWpW7cuN910E927d8fMGDRoEGeffTZ79+5Vey9cfPHFvD9vPnt3/srDA5M486rbyPX/FMUpQ0fT9rSz+H7RPCYNOZkq1Y9i6D2T86Y9/fTT+e6779izZw9NmzblxRdfZMCAAZGqikilouQqBJKSkkhLSysy/L777st7PW/evFLnkZ3rGD16NKNHjwbg6KOPLtCFSiTaJCUlMW76J0WGn3X1+LzXA6+7i4HX3QX4nmB2+KlnVapV58a3Pi92voNuvIdBN96jJ55J1Bo0aBCDBg0qMCx/e59f4Se4XXrppVx66aUFhqm9F4Dp06eX+mRIM2PIHY8WO+6zzz4LVVgiUgYlV1FKj9sVERGRQGTnOhLivD/UItDyIlI2JVciIiIiFUCgJ2Z1UlYk+PRACxGJCdm5gf/OT3mmERGpLAJtI9WmipTNAvhhQu1RYVba2ac7ujQIYyQi3j28fLuncvnvuQpF+cPTSFQLVX+kCne80u8XVU5e29PDQt2uqk2VSsxzI6wrV2Gisz0iIlKYjg0iIhWL7rkKE/WDFgk/3dwt0a48Dy8SiRS1qSJlU3JVTmowRKJfoF9cb+lcL6D5qx0Qkcok1G0qqF2V2Kfkyi/QnTnSV6KK64cdjr7W5bnnRTEFt3xFiSkaRXq/ltgXqi+Gh9v8itBeVIQ6RGtM0aY8V2Z1kktinecHWtx7771zgGj6JtEY2BLpIMpBcYdPLMYMijvcFHf4FI45c8KECQOD/SFlHK9icbkFQvWLbapfbFP9Yltp9fN+vHLOxeTfPffc4yIdg+KO7r9YjFlxK+6KHHc0xBwNMah+qp/qVzH/VL/Y/gtW/fS0QBERERERkSCI5eTq3kgHUE6KO3xiMWZQ3OGmuMMnGmKOhhhCSfWLbapfbFP9YltQ6hfIjwiLiIiIiIhICWL5ypWIiIiIiEjUUHIlIiIiIiISBFGdXJlZdTNbbGYrzOxbMyvSF9LMepnZcjPLNrOhkYizMI9x32Rmq8wsw8w+NrMWkYg1XzxeYr7KzL4xs3QzW2RmHSIRa6GYyow7X9mhZubMrFs4YywhFi/Le7SZbfcv73QzGxuJWAvF5Gl5m9kw//b9rZm9Hu44i4nHy/J+It+yXmNmOyMRa754vMTc3Mw+MbM0f1syKBKxForJS9wt/O1ehpktMLOmQfjcgWb2vZmtM7PxxYwvsc01s5x86z71SGMJBQ/1K7G9MLNRZrbW/zcqvJF746F+Je6fMbL+XjKzbWa2soTxZmaT/fXPMLMu+cbFwvorq34j/PXKMLMvzKxzvnHr8323WBq+qL3zUL8+ZrYr33Z4d75xpW7b0cBD/W7NV7eV/n3uWP+4qF5/ZtbMf5xc7T8mXV9MmeDuf5F+7GFpf4ABNf2vqwBfAz0KlWkJJAGvAkMjHXMAcfcFavhfXw3MjIGYj8n3ejAwJxaWtX9cLeBT4CugWyzEDYwGno50rOWIuw2QBtT1v28YC3EXKn8d8FK0xwxMAa72v+4ArI+FZQ28CYzyvz4DmHqEnxkP/AC0AqoCK4AOhcqU2OYCeyK93IJQv2LbC+BY4Ef//7r+13UjXadA61eofIH9M9rXnz/GXkAXYGUJ4wcBH/r3nx7A17Gy/jzW7w/5jgl/PFw///v1QP1I1+EI69cHeK+Y4QFt29Fav0JlzwXmx8r6AxoBXfyvawFrimk/g7r/RfWVK+ezx/+2iv/PFSqz3jmXAeSGO76SeIz7E+fcPv/br4AjPnN7JDzG/Hu+t0cXHh8JXuL2ux94FDgQrthKE0DcUcVj3FcAzzjnfvNPsy2MIRarHMv7YmB6yAMrhceYHXCM/3VtouDHHT3G3QH42P/6E2DIEX7sycA659yPzrlDwIzC84y2NjdAZdavFAOAuc65X/375Fwg6D/cfIQCrV/E989AOec+BX4tpcgQ4FX//vMVUMfMGhEb66/M+jnnvjh8TCD29j8v668kR7Lvhk2A9Yup/c85t9U5t9z/ejewGmhSqFhQ97+oTq4AzCzezNKBbfgq+HWkY/IiwLgvx5cxR5SXmM3sGjP7AV+iMi7cMRanrLjNLAVo5px7LyIBlsDjNnKB/xL1LDNrFuYQi+Uh7pOAk8zsczP7ysyi4ouA133S313sBGB+OOMrIZayYr4HuNTMNgEf4DujH3Ee4l4BXOB/fR5Qy8zqHcFHNgE25nu/iaIHz/wKt7nVzWypf3v90xHEESpe61dcexHosokEzzGWsH9G+/rzoqRlEAvrL1CF9z8HfGRmy8zsLxGKKRh6mq879Idm1tE/rEKtPzOrgS+5eCvf4JhZf2bWEkjB16Miv6Duf1GfXDnncpxzyfjOcpxsZomRjskLr3Gb2aVAN+CxcMZXHC8xO+eecc6dCNwO/C3cMRantLjNLA54Arg5UvGVxMPyfhdo6ZxLAuYBr4Q7xuJ4iDsBX9fAPvjOcL1gZnXCG2VRAbQlw4FZzrmc8EVXPA8xXwy87Jxriq9bw1T/Nh9RHuK+BehtZmlAb2AzkH0EH2nFhVFsweLb3ObOuW7AJcDfzezEI4glFLzUr6T2wvOyiaBAYixu/4z29edFScsgFtafZ2bWF19ydXu+wac657rg6y54jZn1ikhwR2Y50MI51xl4Cvi3f3iFWn/4ugR+7pzLf5UrJtafmdXElxTeUKgnFgR5/4v4Qdgr59xOYAFReDm8NKXFbWZnAncCg51zB8McWok8LusZQFSdISwh7lpAIrDAzNbj60ubalHwUIvDSlrezrkd+baL54GuYQ6tVKVsJ5uA2c65LOfcT8D3+JKtqOBh+x5OlHV5KCXmy4E3/GW+BKoD9cMaXClK2ba3OOfOd86l4GsDcc7tOoKP2gTkv7LblGK6SJbU5jrntvj//+iPN+UIYgmFMutXSnvhadlEWCAxFtk/Y2D9eVHSMoiF9eeJmSUBLwBDnHM7Dg/Pt/62Ae/g60oXU5xzvx/uDu2c+wCoYmb1qUDrz6+0/S9q15+ZVcGXWE1zzr1dTJGg7n9RnVyZWYPDZ7zN7CjgTOC7yEZVNi9x+7uq/RPfQT7i96R4jDn/F+SzgbXhi7B4ZcXtnNvlnKvvnGvpnGuJr6/3YOdcRJ9o43F5N8r3djC+fsIR5XGf/De+hwfgP7ichO8m0Ijx2paYWVt8N61+Gd4Ii/IY8wagn79Me3zJ1fZwxlmYx227fr4rbHcALx3hxy4B2pjZCWZWFd8XgAJPjSupzTWzumZW7XBcwKnAqiOMJ9i81K+k9uI/QH9/PesC/f3DokmZ9YPi988YWX9epAIjzacHsMs5t5XYWH9lMrPmwNvAn51za/INP9rMah1+ja9+xT6xLpqZ2fFmZv7XJ+P7fr0Dj9t2LDCz2vh6GszONyzq159/vbwIrHbOPV5CsaDufwlBij1UGgGvmFk8vg31Defce2Z2H7DUOZdqZt3xZcp1gXPN7F7nXMdS5hkOZcaNr0tKTeBN//64wTk3OGIRe4v5Wv+Z3yzgNyAaHgnrJe5o5CXucWY2GF93qV/xPQ0s0rzEfbgxWgXkALfmP0sZIV63k4uBGc65aOi24SXmm4HnzexGfF0VRkdB7F7i7gM8bGYO31M8rzmSD3TOZZvZtfi2vXh8T5L71mOb2x74p5nl+uOd6JyLqi/nHutXbHvhnPvVzO7H9yUP4L5CXXoizmP9oPj9M+rXH4CZTce33dc33z2SE/A97AXn3HP47pkcBKwD9gGX+cdF/foDT/W7G6gHPOvf/7L9XTmPA97xD0sAXnfOzQl7BcrgoX5DgavNLBvYDwz3b6fFbtsRqEKpPNQPfPfHfuSc25tv0lhYf6cCfwa+Md+9wAD/D2gOodn/LPLHYRERERERkdgX1d0CRUREREREYoWSKxERERERkSBQciUiIiIiIhIESq5ERERERESCQMmViIiIiIhIECi5EhERERERCQIlVyIiIiIiIkGg5EpERERERCQIlFyJiIiIiIgEgZIrERERERGRIFByJSIiIiIiEgRKrkRERERERIJAyZWIiIiIiEgQKLkSEREREREJAiVXIvmY2Z4Qz/8oM1toZvGllKlqZp+aWUIoYxERkdhjZhea2Woz+8T/frqZZZjZjQHOp46Z/TXAaV4wsw6BTCNS2ZhzLtIxiEQNM9vjnKsZwvlfAyQ4554so9wEYJ1zblqoYhERkdhjZnOAR5xzn5jZ8cDXzrkW5ZhPS+A951xikEMUqdR05UqkGGZ2k5mt9P/dkG/4XWb2nZnN9Z8tvCXAWY8AZvvnVcfMfs4372VmVtv/9t/+siIiUgmZ2aVmttjM0s3sn2YWb2Z3A6cBz5nZY8BHQEN/mdPN7EQzm+M/nnxmZu388zrOzN4xsxX+vz8AE4ET/dM+Vuizjzaz9/1lV5rZRf7hC8ysm5kN9k+Xbmbfm9lP/vFd/b0zlpnZf8ysUTiXmUg0ULcjkULMrCtwGXAKYMDXZrYQiAcuAFLw7TvLgWUBzLcq0Mo5tx7AObfTfwCr4pzLAlYAScBnwEqge9AqJSIiMcPM2gMXAac657LM7FlghHPuPjM7A7jFObfUzJ7Bd/Up2T/dx8BVzrm1ZnYK8CxwBjAZWOicO8/fLb0mMB5IPDxtIQOBLc65s/3zrZ1/pHMuFUj1j3sDWGhmVYCngCHOue3+hOxBYExQF45IlFNyJVLUacA7zrm9AGb2NnA6viu9s51z+/3D3z08gZm1Au4EajvnhprZ0fgOaoeABf7uffWBnYU+6xfgeGAj0M7/HudcjpkdMrNazrndoauqiIhEoX5AV2CJmQEcBWwrbQIzqwn8AXjTPw1ANf//M4CR4Du+ALvMrG4ps/sGmGRmj+BL3j4r4TNvA/Y7554xs0QgEZjr//x4YGsZ9RSpcJRciRRlAQ7HOfcjcLmZzfIPOh+Y5Zx718xmAtOA/UD1QpNuARr7u2hkOufW5BtXDThQngqIiEhMM+AV59wdAUwTB+ws4UpUQJxza/y9OAYBD5vZR865+woEaNYPuBDolS/mb51zPY/080Vime65EinqU+BPZlbDfwXqPHxd9RYB55pZdf8ZwrNLmUdTfFejAHIAnHO/AfFmlj/B2gL8CV/3jLyuE2ZWD9ju7y4oIiKVy8fAUDNrCGBmx5pZqQ+tcM79DvxkZhf6pzEz65xvflf7h8eb2THAbqBWcfMys8bAPufca8AkoEuh8S3w9c4Ydrg3B/A90MDMevrLVDGzjgHWWyTmKbkSKcQ5txx4GVgMfA284JxLc84twdfHfAXwNrAU2FXCbDbhS7Cg4H72Eb5uh4dtBoYCg51zmfmG9wU+OLKaiIhILHLOrQL+BnxkZhnAXMDLwyFG4OtFsQL4FhjiH3490NfMvsF3r3BH59wO4HP/AyseKzSfTsBiM0vH1+X9gULjRwP1gHf8D7X4wDl3CN/x7BH/56fj66YoUqnoUewiATCzms65PWZWA98Vrr8455b7rzQ9CJwFvIDv5uGn8XXrW3T4kepmlgLc5Jz7cxmf8zZwh3Pu+xBWR0RERESCSMmVSADM7HWgA757p15xzj1cjnmM8U+bU8L4qsBw59yrRxSsiIiIiISVkisREREREZEg0D1XIiIiIiIiQaDkSkREREREJAgC+Z0r9R8UEZFgKvG3446QjlciIhJMno9XunIlIiIiIiISBEquREREREREgkDJlYiIiIiISBAouRIREREREQkCJVciIiIiIiJBoORKREREREQkCJRciVQC2bmBPZk60PIiIlI6tcMilYM553nn1V4uEsMmpmV6Ljs+pX4IIxHJo9+5kkpF7bBIzNLvXImIiIiIiISTkisRERGRKFOeboHqSigSeQmRDkBEAped60iIC1WPqvLNP9QxiYhUJglxFlA3QlBXQpFooORKJAYFetAN9ICrg7qIiIhI4NQtUEREREREJAiUXImIiIiIiASBkisREREREZEgUHIlIiIiIiISBEquREREREREgkDJlYgERaC/r6LfYxEREZGKRo9iF5GgCPXj4UVERESina5ciVQQ235ay7Oj/sjfTmnCp68+U2b51EfGM+HUFnnvv571Mn8f1ovJw/tw2mmn8cuP3wOw9qsFPHVJP/4+rBdPXdKPHxZ/FrI6iIhIQc45Uh+9g8cGd+fJYb3ZvHpFseX+8/SDNGvWjJo1axY7ftasWZgZS5cuzRu2PH0FPXv2pGPHjnTq1IkDBw6EpA4ilYmuXIlUEDVq1+Hc2x5i1ScflFl206p09u/+vcCwzgMv4JShowHosPELbpl4F2OeeYOj6xzLqCencUyD4/l53Wr+dc0w7vjPN6GogoiIFPL95/PYseFHbpm9mI3fLOPfD9/GNa/+p0i59r0G8PoDt9HyxNZFehEc3LuHlx/6P5olduXl73cyLz6TnOxsZl7+Z6ZOnUrnzp3ZsWMHVapUCVe1RCosXbkSKaf169fTrl07xo4dS2JiIiNGjGDevHmceuqptGnThsWLF7N3717GjBlD9+7dSUlJYfbs2XnTnn766XTp0oWULl344osvAFiwYAF9+vRh6NChtGvXjhEjRuCct3uTah7bgGYdU4hPKP3gmJuTw4d/v4c/Xn93geHVa9bKe713717MDIDG7ZI4psHxABx3YjuyDh0k+9BBbwtJRCTGBKtt71JM2z7t1st4/PyezLjzKs9t++oFc0g55yLMjOZJ3Tiwexe/b/+5SLnmSd1o1KhRsfP46NmH6TXqOhKqVcsbtvarT0hKSqJz584A1KtXj/j4+ICWlYgUpStXIkdg3bp1vPnmm0yZMoXu3bvz+uuvs2jRIlJTU3nooYfo0KEDZ5xxBi+99BI7d+7k5JNP5swzz6Rhw4bMnTuX6tWrs3btWvoOuZBrp83jx7W7+GrZcm58cxFdGhzPc5edzdX/ep+WKT14b9Lf+HHpIgAaHpXAtv3ZACQNOI8+l13vOeYvZ75A+14D8xKmguNeZNG056hBNhc9PavI+JUfv0vjtp1IqFqtyDgRkYoiWG37xRdfnNcNLy0tjb/O/Ixa/rb9v+lfF2nbC7jsUhgwll3btlLnuMZ5g2s3bMzv238utg0vzpbvMtj1y2ba9+rPZ1P/12U8878/UNeMAQMGsH37doYPH85tt912ZAtORJRciRyJE044gU6dOgHQsWNH+vXrh5nRqVMn1q9fz6ZNm0hNTWXSpEkAHDhwgA0bNtC4cWOuvfZa0tPTiY+PZ9tPa/Lm2axjF2r7D6SN2yby25aNtEzpwTm3PJBXZnxK/YAeHnHY79t/5pt5qVwxZXax43tedDk9L7qc5qs/YtILjzPsvv8diH/54TvmTL6fMc+8EfDniojEkmC17WvW/K9tP/nkk8ts2/P7Xztf9AqXYZ7qkZuby3v/dxcX3vtU0XE5OSxatIglS5ZQo0YN+vXrR9euXenXr5+neYtI8ZRciRyBavm6WMTFxeW9j4uLIzs7m/j4eN566y3atm1bYLp77rmH4447jhUrVpCbm0u16tXzxsVXqZr32uLiyM3xXaHKf3ZzRjmvXG35LoMdG39i0pCTAcg6sJ/HBnfn1tQlBcoNHz6cMVdelfd+1y9bmHrzKC6872nqNTvB02eJiMQqL237zDdn0bF9uwLTFW7bq+dr2/PPs6S2vQD/lavaDRuz85cteYN3bdtCrQbHearHob17+OWH75hyxZ8A2LNjG6/ecCkj//4atY9rTO/evalf3/fk1kGDBrF8+XIlVyJHSMmVSAgNGDCAp556iqeeegozIy0tjZSUFHbt2kXTpk2Ji4vjlVdeITcnp8x5BePKVbvT+3Pn3FV57yec2iIvscrc8AP1m58IwPvvv0/9Zq0A2L97Fy+Pu4SB1/2NlsmnBPyZIiIVzYABA/jHM0/TZMwEzIwt32XQuF0S89b+TO3jGvPoil9ZOvt1cnJymJiWyY9rd5U4r7KuXLXvPYAvZ75I5wHnsfGbZVSveYznLoHVax3DXfO/z3s/5YohDLrxXpp2SKZe05a8+8Y/2LdvH1WrVmXhwoXceOONgS0IESlCD7QQCaG77rqLrKwskpKSSExM5K677gLgr3/9K6+88go9evRgzZo1VD2qxhF/1u7MX3h4YBKLpv2DT154nIcHJnFgz27Ad0ayuBug8/ty5os8MfQ0Jg/vw+OPP86F9z3tH/4COzb+xPzn/4/Jw/sweXgf9vy6/YjjFRGJVYfb9icv6s3fLzydj56dCECPYZex/N2ZPDtyIJkbfghK2972tLM4tkkLJg05mbcfuIkhdzyaN27y8D55rz/8+700bdqUrAP7eXhgEvOee7SYuf3PUcfU4aabbqJ79+4kJyfTpUsXzj777COOV6SyM69Pq6G4Tr8iEhSBXoUK9MpVqMv///buPDyq6v7j+OckkwSCBlC2BFCMCRASIbIoKMhm2CwBNSL8cEFBqqKIQgVri0C1Yq2iCBbFBWwR0FKKC0tBFBAKiCSimCIoCgGUHUQCWbi/PxLGrGSGuTOTybxfzzNPMnPvnZx7cuac851z7rnn+zcQ9Fy7cMR9tFfwCephIGi43F4xcgUAAAAANiC4AgAAAAAbEFwBAAAAgA0IrgAAAADABgRXgM3yznAtPQAAQDDiPleAzRwhhtWaAAAAghAjVwAAAABgA4IrAAAAALABwRUAv3D32jSuZQMAAJUd11wB8AuuTQMAAFUNI1cAAAAAYAOCKwAAAACwAcEVAAAAANiA4AoAAAAAbEBwBQAAAAA2ILgCAAAAABsQXAEAAACADQiugApw81oAAAC4gpsIAxXgZrcAAABwBSNXAAAAQcjdmRnM5AAqxsgVAABAEGJmBmA/Rq4AAAAAwAYEVwAZMH2MAAAWdklEQVQAAABgA4IrAAAAALABwRUAAAAA2IDgCgAABDVWwQNgF1YLBAAAQc3dVfMkVs4DUDZGrgAEhPP5ZplvowEAgC8xcgUgIPDNMgAAqOwYuQIAAAAAGxBcAQAAAIANCK4AAAAAwAYEVwAAAABgA4IrAAAAALABwRWCDstzAwAAwBtYih1Bx90lvVnOGwAAAK5g5AoAAAAAbEBwBQAAAAA2ILgCAAAAABsQXAEAAACADQiuAAAAAMAGBFcAAAAAYAOCKwAAAACwAcEVAAAAANiA4AoAAAAAbEBwhYCWd8bydxIAAAgK7ra5tNEIRg5/JwDwhCPEaHL6QbeOGXdlHS+lBgCAqsvdNpf2FsGIkSsAAAAAsAHBFQAAAADYgOAKQJXF9QEAAMCXuOYKQJXF9QEAAMCXGLkCAAAAABsQXAEAAACADQiuAAAAAMAGBFcAAKBKYXEaAP7CghYAAKBKYTEbAP7CyBUAAAAA2IDgCgAAAABsQHAFAAAAADYguEKlwkXIAAAACFQsaIFKhYuQAQAAEKgYuQIAAAAAGxBcAQAAAIANCK4AAAAAwAYEVwAAAABgA4IrAAAAALABwRUAAAAA2IDgCgAAALY7n3tXcr9LBDrucwUAhfLOWHKEGK/tDwDBxN17V0rcvxKBj+AKXkXnE4GEm1gDAABPEFzBq+isAgAAIFhwzRVgo7Vvv6KkpCRNSeuoT+fMcL6+YsZf9HTPKzR1YBdNHdhF//t0uSTp+4wNatmypabdlqKDu76TJGX/fExv3H+LLKvseeev3tNPmzZtcj4/sneXXrilkyTpu01rNeG6WE0d1FXP33SNVrzyrCTpk08+cb7+3I3t9crQvspc/R+v5AEABIMXX3xRL9zSyeX6fu3atXpxQGe36/usrzOcz12p74u+Tn0P+B4jV4BNftyRqc8W/kPfbflcU7Ye15sP3KrmnVJU55LLJUnXDr5X190xotgxn/79b1q4YIEmr/hCG/45Szc8MkkrZz6nLkNHyZjzm07ZJLm9hkx9WznZv2jqwK5K6NRDahjqfF2S9m77Un9/5E6FRVRT3NXXeXbiABBkftyRqRUzZ+r+t5YpNCzcpfr+ueee0+C/vqEje3d7rb7//PMBxV6XqO8BX2PkCrDJgZ3fqPEVbRQZGalQh0OXtblGW1cuPucxIQ6HsrOzlXsqW6EOhw7t3qnj+/cpts21HqcnvHoNNUxopUNZO0tti2l2hbrfM1r/fed1j/8OAASbAzu/Ufv27RVe3fX6PiwsTLmnTnm1vv/2229LbaO+B3yLkSvAJvUvT9Cy6X/WoUOHlJN9Uts+XaFGLVo5t/93/utK/+AdNWzRSjc8MknVo2qpy90Pafjw4TqQ79CAP72sxVOeUMr94yr8W4MHD9bPCpMk5efmyoSU/p7kl6OHtevLz9XtntGS8kptj0loqdVvTT//EwaAIFX/8gR98Nozih58WGER1Vyq7x977DGl3j5MYdWquVXfz3/8XoVFVJNUcX2fmPikVqzbUWo79T3gOwRXgE3qxTZV5yEPKiUlRUcUoeimiQoJLfiIXX3LkIIgxxgtf/lpffj8eKVNmKqYZldo/fr1mpx+UDs/X6cL6zaQZUlvjx2mUEeY+jwyURdeXK/U35ozZ45WhDaRVDAHf/ZDg53bvs9Yr6mDusqYEHW5a6TqX95cOvZV6QSXM8cfAHBu9WKbauzYsRp/f5rCq9dwqb5PTk7W/W8tlSS36vtbn5qhRi2SJVVc3ycmJkplBFfU94DvMC0QsFG7/rdp8+bN+u3r7ysyqpYuviRWknThxfUUEhqqkJAQXXXT7craml7sOMuytPL1Kep+z2h99Oqzuv7esUruk6Z1c2e6nYYmye01cu7HevDtj3R12pBy99v7vy9V77Kmbr8/AEAaOnSoHnx7JfU9gGIIrgAbnTh8QJJ0dF+Wtn78oZJ73SRJOn7gR+c+W1cuLhhNKmLz+/PUvOP1qh5VS7mnsmVCQmRCQpR7Ktsr6dz3zVatfO15tR9wt1feHwCquv3790uivgdQHNMCARvNGXOX/nX6uI7khSh17DOqHlVLkrTkxUna981XMjKqHdNY/R//q/OYkydPavMH83X39HclSR0H36s5v7tLoY4wDXz6VdvSdnb6SO6pbF1Qu476/u7PrBwFICBUxhvS33zzzfpmz36FOMJcru9zsqnvgaqO4Aqw0W/f+EDjrqxT6sbJtz75crnHREZG6p5X/+18flnrDhr1zupy9x8+c5HaXllHKwr/Ru2YSzTq3TWSpNi21yq2bemVp7p06aIJq79z61wAoLKojDekX7NmTZlpOld9H17d/fq+KFfq+9i211LfA37EtEAAAAAAsAHBFVyWd4bVhgAAAIDyMC0QLnN3Wobkm6kZAAAAQGXAyBUAnCd3R3MZ/QUAoGpj5AoAzlNlvMgeAAKZuytDVsaVJBHcCK4AAABQKfClFQId0wIBAAAAwAYEVwAAAABgA4IrAAAAALABwVUAWbp0qZo1a6a4uDhNnjy51PaHH35YycnJSk5OVtOmTVWrVi3nttmzZys+Pl7x8fGaPXu2JOnIsePO/ZOTk1WnTh2NGjXKZ+cDBLrdu3dr5vD+ev6mazQlraPWvv1KqX0sy9J7f3lMz6a2U8uWLbV582bntrFjxyopKUlJSUmaP3++8/VOnTo5P5cxMTHq37+/T84HkArK7MiRIxUXF1eqzJ518uRJ3XDDDWrevLkSExM1btw457by2qIffvhBbdq0UXJyshITEzVjxgyfnROqhm1rP9JzN7bXs6nt9MmbL5baPmvWLNWtW9dZ/l577TXntkcffVSJiYlKSEjQyJEjZVms3grvYEGLAJGfn68RI0Zo+fLlatSokdq1a6fU1FS1aNHCuc+UKVOcv7/00ktKT0+XJB0+fFgTJ07Upk2bZIxRmzZtlJqaqtq1a2vgmyt+Peb/uis7qVu5F5Jy0ShQnMPhUJ+HJ6phQiud/uWEXhrcXXHtu6h+bDPnPtvWrtChXd9pzKKN6przre677z5t2LBBH374oTZv3qyMjAydPn1anTt3Vu/evRUVFaU1a9Y4j7/55pvVr18/f5wegtSSJUu0fft2bd++XRs2bHCW2ZLGjBmjrl27KicnR927d9eSJUvUu3fvctui6OhorVu3ThERETpx4oSSkpKUmpqqmJgYn50bAteZ/Hy998w4DX35XUXVj9H023oooXMvqUTf5NZbb9W0adOKvbZu3TqtXbtWW7ZskSR17NhRq1atUpcuXXyVfAQRRq4CxMaNGxUXF6fY2FiFh4dr4MCBWrRoUbn7z507V4MGDZIkLVu2TCkpKbroootUu3ZtpaSkaOnSpcX2P7jrW/1y5KCatO7g1fMAqpLo6Gg1TGglSYqocYHqXdZUx/fvK7ZP5idLdeVvbpUxRu3bt9fRo0e1b98+ff311+rcubMcDodq1KihVq1alfpc/vzzz1q5ciUjV/CpRYsW6Y477ihVZouKjIxU165dJUnh4eFq3bq1srKySr1X0bYoPDxcERERkqTTp0/rzJkzXj4TVCW7v9qsixs10UWNmsgRFq5WPfsr85MlLh1rjNGpU6eUk5Oj06dPKzc3V/Xr1/dyihGsCK4CxJ49e9S4cWPn80aNGmnPnj1l7vvDDz9o586d6tatm8vHfrF0oVr26C9juFcEcD6O7N2lvdu+VOOkNsVeP7Z/n2rV//Wb+bOfv1atWmnJkiU6efKkDh48qI8//li7d+8uduzChQvVvXt3RUVF+eQcAMm99kaSjh49qvfff1/du3cv9nrJtkgqmEp7RcuWaty4scaOHcuoFVx2/MA+1WzQ0Pk8ql6MjpX4MkuSFixYoJYtWyotLc1Zp3bo0EFdu3ZVdHS0oqOj1bNnTyUkJPgs7QguTAsMEGXNDS4vEJo3b57S0tIUGhrq8rFbli3UgD+9bENKgeBz+uQJ/WPMXfrN6CdV7YILS2wt+/PXo0cPffbZZ7rmmmtUt25ddejQQQ5H8Sp57ty5GjZsmBdTDpTmTnuTl5enQYMGaeTIkYqNjS22rWRbJEmNGzfWl1u26Pf/+UqTH7lDe5t11YUX16swTUxLh8opl0VvIty3b18NGjRIERERmjFjhu68806tXLlSO3bsUGZmprKyspR3xlLvnj20evVqXXfddb4+CwQBRq4CRKNGjYp9q52VlVXqG7+8MwUVz7x585zTMFw5dt83Xyk/P08NW7TyVvKBKis/N1dzxtyl5D5pSur+m1Lba9aL0dGf9koq+IwW/fw9/vjjysjI0PLly2VZluLj453HHTp0SBs3blTP3n18cyIIatOnTy+2iEpF7c1Zw4cPV3x8fJmLIZVsi4qKqttA9WOb6fv09facAKq8qHoxOvbjryOox/fvVVTdBs6bDk9OP6iZuyxN+fpnTU4/qMNtb9S6zzZpcvpBjZz2D+Vc2lLTtp9SragL1bt3b61fT9mDdxBcBYh27dpp+/bt2rlzp3JycjRv3jylpqYW28cRYjR64Xp9/9MhraoW76xsMqPb6t0PluqJVTv0xKodeveDpcqMbus87oul/1Krnjf5+pSAgGdZlhZMGqW6lzVVp9vuK3OfhM49lf7BfFmWpU0bN+hkWA3N/jFMf970k/748TeanH5QD81fpZUb05Vet7Xzczvs+Td16TXX64LI6j4+KwSjESNGKCMjQxkZGerfv7/eeustWZal9evXq2bNmoqOji51zB/+8AcdO3ZML7zwQqlt27Zt05EjR9Shw6/X8WZlZSk7O1uSlH38qL7/YqPqXhrnvZNCldIo8Uod3L1Th/f8oLzcHH2x7N8FC1oUcfzAj87fM1ctVb0mTSVJtRo01M7P1yk/L0+5ublatWoV0wLhNUwLDBAOh0PTpk1Tz549lZ+fr7vvvluJiYkaP3682rZt6wy0CgKl4tdORdasrW7DHtG021IkSd3uGa3ImrWd27csf09Dps717QkBVcDatWuV/uE7ahDXQlMHdpEk9Xjgcee3q1enDVGzjina9ukK/bXfVXq71gXq91jBSmr5ebl6dWhfSVJEjQs14MmXFVpkWuAXyxaq85CRvj0hQFKfPn20ePFixcXFKTIyUm+++aZzW3JysjIyMpSVlaWnnnpKzZs3V+vWrSVJDzzwgHMa69y5czVw4MBibVFmZqZGjx4tY4z2n8zTdbePUIP4FgJcEepwKHXs03pjxABZZ86obeog1b+8ucaPH6/vL2qqFp17ad28mcpctUwhoQ5F1qyltIkvSZKSrk/Vt599qhcHXKc51R3q1auX+vbt6+czQlVFcBVA+vTpoz59ik8RmjRpUrHn19/7aJnHtu0/WG37Dy5z26Pvb7IngUCQ6dixo57efOCc+xhj1O+xv0gquG7k7K0OwiKq6eEFa8s9bvjM8lcDBbzJGKPp06eXuS0jI0NSwXTzc90naMKECaVeS0lJcS6FXd4tP4Bzad4xRc07phR7bdKkSc7y1OvBP6rXg38sdVxIaKhu/MNzkrh+D97HtEAAAAAAsAHBVSV1dnEKAAAqM9orBJLzKa+UcbiDaYGV1NnVb9zBUDcAwNfcba9oq+BP9K/gbeZcc6ZLIGz3sfP58LvbwJXc/7HWdd36m0Cwq+iaq6Lc/YxK0phWFzvv4eKKovd8CQDeSijtlY+dT9tDewN3lVXf2tH3ceUYBD2X2ytGrnwkwDo8ACoJRgUAwL/c7cPR5wtuBFc+QgcJABAI6BgCxdGHgzsIrs4TjQ+AyohvWOEpOpKAZ6iHg5vL11xNnDhxqSRv1KAxkvZ64X2rCvLn3MifipFH50b+VMxbeXTwiSee6GX3m3qxvaoqKPOeIf88Q/55jjz0zPnkn+vtlWVZfn1MmDDB8ncaKvOD/CF/yCPyx98P8qhqPfh/kn/kX2A/yMPKnX/c5woAAAAAbFAZgquJ/k5AJUf+nBv5UzHy6NzIn4qRR1UL/0/PkH+eIf88Rx56xqv55859rgAAAAAA5agMI1cAAAAAEPAIrgAAAADABgRXAAAAAGADnwVXxphexphtxpgdxphxZWyPMMbML9y+wRjTxFdpqwxcyJ8hxpgDxpiMwscwf6TTX4wxbxhj9htjvipnuzHGTC3Mvy3GmNa+TqM/uZA/XYwxx4qUn/G+TqM/GWMaG2M+NsZkGmO2GmMeKmOfoC1DLuZPUJehQGaMucgYs9wYs73wZ+1y9ssv8v99z9fprGzot3iGfo1n6Pd4xp/9Ip8EV8aYUEnTJfWW1ELSIGNMixK7DZV0xLKsOElTJD3ji7RVBi7mjyTNtywrufDxmk8T6X+zJJ3r5m29JcUXPoZL+psP0lSZzNK580eS1hQpP5N8kKbKJE/SaMuyEiS1lzSijM9YMJchV/JHCu4yFMjGSfrIsqx4SR8VPi9LdpH/b6rvklf50G/xDP0aW8wS/R5PzJKf+kW+Grm6StIOy7K+sywrR9I8Sf1K7NNP0uzC3/8pqbsxxvgoff7mSv4ENcuyVks6fI5d+kl6yyqwXlItY0y0b1Lnfy7kT1CzLGufZVmbC3//WVKmpIYldgvaMuRi/iBwFW1fZ0vq78e0BAr6LZ6hX+Mh+j2e8We/yFfBVUNJu4s8z1Lphtu5j2VZeZKOSbrYJ6nzP1fyR5JuLhz6/acxprFvkhYwXM3DYNbBGPOFMWaJMSbR34nxl8KpO1dK2lBiE2VI58wfiTIUqOpblrVPKgikJdUrZ79qxphNxpj1xphgD8Dot3iGfo330WZ5zittmsOuN6pAWd/klLzBliv7VFWunPv7kuZalnXaGHOvCr4t6+b1lAWOYC4/rtgs6VLLsk4YY/pI+rcKphIEFWPMBZIWSBplWdbxkpvLOCSoylAF+UMZqsSMMSskNShj0+NuvM0llmXtNcbESlppjPnSsqxv7UlhwKHf4hn6Nd5H+fOM19o0X41cZUkq+o1EI0l7y9vHGOOQVFPBM82pwvyxLOuQZVmnC5/OlNTGR2kLFK6UsaBlWdZxy7JOFP6+WFKYMaaOn5PlU8aYMBUEDnMsy/pXGbsEdRmqKH8oQ5WbZVnXW5aVVMZjkaSfzk4XKvy5v5z32Fv48ztJn6hgBDNY0W/xDP0a7wvqNstT3mzTfBVcfSYp3hhzmTEmXNJASSVXInpP0p2Fv6dJWmlZVrBE4BXmT4l5tKkquCYCv3pP0h2Fq+e0l3Ts7DQYSMaYBmevBTDGXKWCz/4h/6bKdwrP/XVJmZZlPV/ObkFbhlzJn2AvQwGuaPt6p6RFJXcwxtQ2xkQU/l5H0rWSvvZZCisf+i2eoV/jfUHbZtnBm22aT6YFWpaVZ4x5QNIySaGS3rAsa6sxZpKkTZZlvaeChv3vxpgdKvjmZ6Av0lYZuJg/I40xqSpY1euwpCF+S7AfGGPmSuoiqY4xJkvSE5LCJMmyrBmSFkvqI2mHpJOS7vJPSv3DhfxJk3SfMSZPUrakgUHWCbhW0u2SvjTGZBS+9ntJl0iUIbmWP8FehgLZZEnvGGOGStol6RZJMsa0lXSvZVnDJCVIesUYc0YFnYzJlmUFbXBFv8Uz9Gs8R7/HM/7sFxnaRgAAAADwnM9uIgwAAAAAVRnBFQAAAADYgOAKAAAAAGxAcAUAAAAANiC4AgAAAAAbEFwBAAAAgA0IrgAAAADABv8P7w0fF6KakQAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x720 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.plot_posterior(trace, varnames=[l_diff_means, l_diff_sd, l_g1_mean, l_g1_sd, l_g2_mean, l_g2_sd, l_nu, l_effect_size], color='#87ceeb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The small sigma value on the group1/2_mean pulls both estimates towards it reducing the estimated difference from `0.34` to `0.09`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with using a more diffuse sigma: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [nu, $\\sigma_{patients}$, $\\sigma_{controls}$, $\\mu_{patients}$, $\\mu_{controls}$]\n",
      "Sampling 2 chains: 100%|██████████| 41000/41000 [01:04<00:00, 633.23draws/s]\n"
     ]
    }
   ],
   "source": [
    "prior_on_mean_sd = y_joint.std() * 2\n",
    "\n",
    "with pm.Model() as model:\n",
    "    group1_mean = pm.Normal(l_g1_mean, mu=muPriorMean1, sd=prior_on_mean_sd)\n",
    "    group2_mean = pm.Normal(l_g2_mean, mu=muPriorMean2, sd=prior_on_mean_sd)\n",
    "\n",
    "    group1_sd = pm.Gamma(l_g1_sd, alpha=Sh1, beta=Ra1)\n",
    "    group2_sd = pm.Gamma(l_g2_sd, alpha=Sh2, beta=Ra2)\n",
    "\n",
    "    nu = pm.Gamma('nu', alpha=ShNu, beta=RaNu)\n",
    "\n",
    "    lambda1 = group1_sd ** -2\n",
    "    lambda2 = group2_sd ** -2\n",
    "\n",
    "    group1 = pm.StudentT(l_g1_stat, nu=nu, mu=group1_mean, lam=lambda1, observed=y1)\n",
    "    group2 = pm.StudentT(l_g2_stat, nu=nu, mu=group2_mean, lam=lambda2, observed=y2)\n",
    "\n",
    "    diff_of_means = pm.Deterministic(l_diff_means, group1_mean - group2_mean)\n",
    "    diff_of_stds = pm.Deterministic(l_diff_sd, group1_sd - group2_sd)\n",
    "    effect_size = pm.Deterministic(l_effect_size, diff_of_means / np.sqrt((group1_sd ** 2 + group2_sd ** 2) / 2))\n",
    "    nu_transformed = pm.Deterministic(l_nu, np.log10(nu))\n",
    "\n",
    "    trace = pm.sample(20000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<matplotlib.axes._subplots.AxesSubplot object at 0x1c27ac8080>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c29a48940>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x110071ac8>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x110097da0>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c283300f0>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28340668>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28364940>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x1c28387c18>],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAALICAYAAABijlFfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FVXewPHvSUIRERUQNaKLFOkhAVFcZAUVsazYAEHWgm0tWNbKvjZcG772trq6qKg0yyqoWADXF8EVFglSFRAQIu4KSFGKEDLvH4l3yYaSwKR/P8+TJ/fOnJl7ztx7z7m/OWfOhCiKkCRJkiTtvqTSzoAkSZIkVRQGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLCk3RRCOCOEEIUQmpV2XiRJ2hnbLal4GWBJu68PMBXoXdoZkSSpEGy3pGJkgCXthhBCTeAY4CJyG6xfln8cQmia97hOCGFWKWVRkqSEHbRbbUIIE0IIc0IIOXk9XHeWWkalciyltDMglXOnA+OiKJoRQlgXQmgbRdE0oDEwPy9NGjCz1HIoSdJ/FGi3gDnASOC8KIqmhBDuAqoDA0sxn1K5ZQ+WtHv6AK/mPX4V6BNC+BXwbRRFOXnL04AZpZE5SZL+S4F2CzgemBZF0ZS85TOA2lEURaWQP6ncswdL2kUhhDrAEcCZeYtGAv8HTCR/QNUub50kSaVmB+3WD+QfadEWmFayuZMqDnuwpF3XAxgTRdHPAFEULQL+BaSTO7SCEEIT4DQcIihJKn3ba7dSyB1tQQjhMHIDsBGllUmpvLMHS9p1fYC0EMLirZbVAdYBy0IIX5DbkzUXOB+4q8RzKEnSf2yv3ZoN1MybkGkF0CeKopWlkD+pQggOr5XiFUJYAGREUfRjaedFkiRJJcshglKMQgh7ATkGV5IkSZWTPViSJEmSFBN7sCRJkiQpJgZYkiRJkhSToswi6FhCSVKcQgm8hm2XJClOO2277MGSJEmSpJgYYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuKWXZO0WeF3pVtJEmlzzpf0n8LUVToL7m1gVRIgzJXFCn9gIy6xZQTqUzzPlgqc7JzIlKSivbRtM6XKpWdVhBFudGwJElShZaSFIoUMBksSfpvDhGUJEmSpJgYYEmSJElSTAywpDKgqBc8e4G0JElS2eQ1WFIZ4Jh/SZKkisEeLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEk74ZTokiRJKiynaZd2Is4p1KMo4u0H/oevJo6javUa9LjzcQ5q3qZAuuev7MWPK74nZ0s2DTI6cNqA+0lKTmbMIwP58pMPGLbXHjRq1IgXXniBffbZh8WLF9O8eXOaNm0KQIcOHXjmmWeKXlhJ0m754dtvGP7HS9mwZhWpzdLodfefSalSNV+asWPHMmDAADZt2kTVqlV54IEHOPbYY1m/fj09e/bk66+/Jjk5mVNPPZVBgwYltnv11VcZOHAgIQTatGnDsGHDSrp4kgrBHiypBH01aRwrlyzkhlFTOOPWh3jrvpu2me6c+wdzzciPufa1T1i3agUzx40GoHGHY7jm1U+YMWMGhx12GPfdd19im0aNGjF9+nSmT59ucCVJpeT9x//E0X0v44ZRU9ij1j5MfWtogTR169bl7bffZubMmQwZMoRzzz03se6GG27gyy+/JDMzk0mTJvHee+8BMH/+fO677z4mTZrE7NmzefTRR0usTJKKxgBLFdrixYtp1qwZF198Ma1ataJv376MGzeOjh070qRJE6ZMmcK6deu48MILad++PRkZGYwaNSqxbadOnWjbti1PnHMs33wxBYCFUyfx7CWnMfTGfjx85lGMuOUyoqhwwwjnfvw+Gb89mxACh6QdzsYf17B2+b8KpKtecy8AcrKz2bJ5MyFv+WFHdSE5JbfjuUOHDmRlZe3mEZKk8iWuer1t27Z8+umnAHz88cd07tyZHj160KxZsyLV61uLooiv/zmRVsedCkDb357NnL+PKZAuIyOD1NRUAFq2bMnGjRv5+eefqVGjBl26dAGgatWqtG3bNlHPP/fcc1x55ZXsu+++ANSrV6/I+ZNUMhwiqApvwYIFvPbaazz77LO0b9+eYcOGMXHiREaPHs29995LixYtOPbYY3n++edZvXo1RxxxBMcffzz16tVj7NixVK9enRtGTWbEH39P/6HjAFj21Uz+8NpE9trvAJ7pdwrfTJ9Mg4wOvPPgrYyY/Rnfb8jOl4e0bmfQud81rPn+O/bZPzWxfO96qXkBVqsC+X7+ip4snZ1J047H0er47gXXP/88Z599duL5okWLyMjIoFatWtx999106tQppiMoSWVLHPX6/Pnz6dOnD1OnTgUgMzOT2bNnk5qaSoP0I/PV6wunTiyQh1/q9a2tX/0D1WvWSpwI23v/1G2eRNvaG2+8QUZGBtWqVcu3fPXq1bz99ttcc03ua8ybNw+Ajh07smXLFgYOHMiJJ564awdQUrEywFKFd+ihh9K6dWsg90zhcccdRwiB1q1bs3jxYrKyshg9ejQPPvggABs3bmTJkiWkpqbSv39/pk+fzspNESuWLEzs8+CWbdk7L1BKbdqKVcuW0iCjA7+94W4GZNTdwTVbBc+IhkT/VH4X/vk1Nv+8kZG3XMbX//yEJh06J9bdc889pKSk0LdvXwAOPPBAlixZQp06dfj88885/fTTmT17NrVq1Srq4ZKkMi+Oej05OTkRtAAcccQR1K9fHyhYrxdWtI06nrDtOh5g9uzZ3HzzzXz44Yf5lmdnZ9OnTx+uvvpqGjZsmFg2f/58Pv74Y7KysujUqROzZs1in332KXT+JJUMAyxVeFufFUxKSko8T0pKIjs7m+TkZN54443EBBG/GDhwIPvvvz9ffPEF933+PbcfVT+xLnmrC5ZDUhI5W3J7rHbWg7V3vVRW/3tZYvma75ex1377bzfvVapVp/kxJzLn4/cSAdbnb48g6713GD9+PCGv4a5WrVqiXO3ataNRo0bMmzePww8/vNDHSZLKizjq9ZycHKpXr77Nff53vb6jHqxu3boxc/G3HNQinTNve4SNP61lS3Y2ySkprPn3MmrV3XYdn5WVxRlnnMFLL71Eo0aN8q279NJLadKkCddee21iWf369enQoQNVqlTh0EMPpWnTpsyfP5/27dsX9rBJKiEGWKr0unXrxhNPPMETTzxBCIHMzEwyMjJYs2YN9evXJykpicx3XyVny5ad7mtnPVjNj+nGP0YOpk23M1g683Oq16xFrf0OyJfm5/U/8fO6n6i13wFsyc7mq4njaJDRAYCvJo1nwotPMGvyRGrUqJHYZvny5dSuXZvk5GQWLlzI/PnzE2c9JamyKUy9PmTIELYUsl7fkQ8++CBfnd/w8I7MGv82bbqdwbR3RtK880kFtlm9ejWnnHIK9913Hx07dsy37tZbb2XNmjX89a9/zbf89NNPZ/jw4VxwwQWsWLGCefPmWc9LZZSTXKjSu+2229i8eTNpaWm0atWK2267DYArrriCIUOG0KFDB1Ys+Zqqe9TYyZ52runRXal90K948LQj+Nvd13HaH/83se7x3p0B2LRhPS/94Vwe63UMj/fuTM3adTmyxwUAjL5/AD+v/4muXbuSnp7OZZddBsCECRNIS0ujTZs29OjRg2eeeYbatWvvdn4lqTwqTL0+b9489txzz9hf+6Srb2fiK0/zQPf2rF/zA+1Pzx3KPef/3mfs07lTrj/55JMsWLCAu+66i/T0dNLT0/n+++/JysrinnvuYc6cObRt25b09PREoNWtWzfq1KlDixYt6NKlCw888AB16tSJPf+Sdl8owiw53m1VlVZR74NVlPS7ss2O7rUllSPbvzglPrZdKrKSqPMllVs7bbvswZIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsqRzKzinadftFTS9JkqRd432wpHIoJSk466AkSVIZZA+WJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSpBLkzeKlis0bDUuSJJUgbxYvVWz2YEmSJElSTAywJEmSJCkmBliSJEmSFBMDLFU6XiwsSZWHdb6kkuYkF6p0vLhYkioP63xJJc0eLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJVUC2TlRsaaXJElSrpTSzoCk4peSFBiUuaLQ6Qdk1C3G3EiSJFVc9mBJkiRJUkwMsCRJkiQpJgZYkiRJkhQTAyxJkiRJiokBliRJkiTFxABLkiRJkmJigCVJkiRJMTHAkiRJkqSYGGBJkiRJUkwMsCRJkiQpJgZYkiRJkhQTAyyVa9k5UWlnQZIkSUpIKe0MSLsjJSkwKHNFkbYZkFG3mHIjSZKkys4eLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEmSJElSTAywJEmSyrCi3vPRe0RKpcv7YEmSJJVhRb3no/d7lEqXPViSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJamAXZmBylmrJBU36xlJ5YGzCEoqoKgzVoGzVkkqftZNksoDe7AkSZIkKSYGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgqUzxJpKSJEkqz7zRsMqUot5E0htISpIkqSyxB0uSJEmSYmKAJUmSJEkxMcBSpTZp2F94tGcnHulxNBOHPpNYPu6Z/+W+bq15vHdn0tPT+XLiWAAWT5/MY72O4cnfdWXFkoUAbPhxDc9f0ZMo2vb1Y89echpZc6Ynnq9atoRHe3YCYOHUSQz8TUMyMjJ4+MxfM+4vD+Rb/nifLjx0Rgf+ctGpzJ3wYbEcA0mq6ApT1z/eu3Ox1/WP9+nCw2f+mjvvvLPAcut6qeLwGixVWv9aMJd/vvkKV7z0AclVqvJC/7Np1qkrdQ9pBEDHvpfxm/OuZEBG3cR1YRNffpq+Dz7PqmVLmfz6i5xy3Z/46LmH6HzRtYQQdikfDdI7kPnJh/zp0294vHcXmnc6IbH8gseHAbDsq5m8fN35VKlWncZH/iaG0ktS5TBr1qxC1fVbK666/oLHh7FpwzpePv94Tmx8dL7lYF0vVRT2YKnSWr5oHge3bkfVPWqQnJLCoe1+zeyPxuxwm6SUFDZv3MjmjRtITklh5dJFrP3+Oxq267jb+am6x54c1LwNK7MWFViX2rQ1x11yPf94dfBuv44kVSZz584tc3V9u3btrOulCsweLFVa+zdqzgdP3cu61T9QpVp1vpo4jvot2iTW/2PkYDLfeZV5nY7k4PP/hz1q7UPnC6/hzbuvp0r16vS668+MeeQOul4xYKevNfKWy6hSrToAWzZvJiQVPLexbvUPLJn5Ocdecj3rVq0ssD61eRoTXnpqN0osSZVPq1atWDRtwE7r+oNatOGU6/4E1C32uv6zzz7jjB79reulCsoAS5VWvYaHccwFV/H8FT2ouseeHHhYS5KSc78SR/a8gGMvuR5CYN3fHuPdh2+nx8DHSW3amiteeh+ARZ9/yl77HUAUwbCbLyY5pQonX3cnUHDq+LPveYb6LdKB3HH5Q67pm1i3ePpnZGRk8P3GHDr3u5r9GzVj4dRJBTO8nXH/kqTta968eaHq+rF/vo93H76dO48ZVqi6fq869Qq81s7q+sf7dCGEJG4bMIDV1vVShWWApUqt/em/o/3pvwPggyfuptb+qQD5Gs5LLrmEwceflG+7KIr4aPAjnDPoOUbdP4DjL7uZVcuW8Onw5+D4R4qUh1+uwdrZ/b+WfTmTeoceVqR9S5IKV9cfcea5+QIi2HFd363/LUXKw9bXWl221bW9/826Xir/vAZLldpPPywHYPV3Wcz++7ukn3gmAGuX/yuR5s0332T/Rs3ybTft7RE0O/p49qi1D5s3biAkJRGSkti8cUOx5PO7ebP56K8P06HXhcWyf0mqyApT18/+aIx1vaRY2IOlSm3oDf1Yv2YVSSlV6H7z/exRax8A3nvsT3w3bxaBQPvmjTjl+nsT22zasJ5p74zkwqdeA+Dovpcx9MZ+JKdUofd9z8aWt1+Gk2zeuIGa+9bl1BvvLdOzSmXnRKQkFX52raKml6RdVZi6ft/Ugzn9lgcT25Tnun5X6lfrZCk+Bliq1H7//DvbXH723X9OPB7wX0M5qu5Rg0uefSvx/NC2R3HtqxO2+xqXPjcq3/N9Uw/h2tc+AaDh4R1peHjBWakaHt6RgRMWFq4QZURKUtjpMMetDcgoeK2aJBWHwtT1/6081/VFrY/BOlmKk0MEJUmSJCkmBliSJEmSFBMDLBWr7Bynm5UkSVLl4TVYKlZelyNJkqTKxB4sSZJUKhzlIKkisgdLkiSVCkc5SKqI7MGSJEmSpJgYYEmSJElSTAywJEmSJCkmBliVRBRFXH311TRu3Ji0tDSmTZtWIM369es55ZRTaNasGS1btmTAgAGJdRMmTKBt27akpKTw+uuvl2TWVQF9NWk8TZs2pXHjxgwaNKjA+iVLltClSxcyMjJIS0tjzJgxBdbXrFmTBx98sKSyLFUYhWkPADZt2sSll17KYYcdRrNmzXjjjTcA+MMf/kB6ejrp6ekcdthh7LPPPoltbr75Zlq1akWrVq0YOXJkiZRHxev999/fYX398MMP06JFC9LS0jjuuOP45ptvEuuWLFnCCSecQPPmzWnRogWLFy8G4KKLLqJNmzakpaXRo0cPfvrpp5IqjlQiDLAqiffee4/58+czf/58nn32WS6//PJtprvhhhv48ssvyczMZNKkSbz33nsAHHLIIbz44oucc845JZltVUA5W7Yw+v4BvPfee8yZM4fhw4czZ86cfGnuvvtuevXqRWZmJiNGjOCKK67It/4Pf/gDJ510UklmW6owCtse3HPPPdSrV4958+YxZ84cjjnmGAAeeeQRpk+fzvTp07nqqqs488wzAXj33XeZNm0a06dPZ/LkyTzwwAOsXbu2xMql+G3ZsoUrr7xyh/V1RkYGU6dOZcaMGfTo0YObbropse68887jxhtvZO7cuUyZMoV69eoBuZ+hL774ghkzZnDIIYfw5JNPlmi5pOJmgFVJjBo1ivPOO48QAh06dGD16tV89913+dLUqFGDLl26AFC1alXatm1LVlYWAA0aNCAtLY2kJD8y2j1LZ02jTv0GHNLgUKpWrUrv3r0ZNWpUvjQhhMQPszVr1pCampqYzvmtt96iYcOGtGzZssTzLlUEhWkPAJ5//nn++Mc/ApCUlETdugVn8Bs+fDh9+vQBSARhKSkp7LnnnrRp04b333+/eAuj2GxryvwpU6bQuHFjGjZsWKC+/iV9ly5dqFGjBgAdOnRI/G6YM2cO2dnZdO3aFYCaNWsm0tWqVQvI7U3dsGEDIYTiLZxUwvy1XEl8++23HHzwwYnn9evX59tvv91u+tWrV/P2229z3HHHlUT2VImsXf4dex9wUGJ65szsvRk9fQGDMlck/uqcdRWP/HUIe++fyrHdTqLdlXeRkhRYt24d999/P3fccUdpF0MqtwrTHqxevRqA2267jbZt29KzZ0/+/e9/50vzzTffsGjRIo499lgA2rRpw3vvvcf69etZsWIFf//731m6dGkxl0Zx+aVO3vrvyYlz+aHGfonnW9fXKUkFg6LBgwcnRhfMmzePffbZhzPPPJOMjAxuvPFGtmzZkkjbr18/DjjgAL788kuuuuqqEiunVBIMsCqJKCp4Zmp7Z4yys7Pp06cPV199NQ0bNizurKmyKcRn8YsP3qTdqb354/szuODx4bx62xXk5ORwxx138Ic//IGaNWuWVG6lCqcw7UF2djZZWVl07NiRadOmcdRRR3HDDTfkSzNixAh69OhBcnIyACeccAInn3wyv/71r+nTpw9HHXUUKSnebrNcK8Jvh1deeYWpU6dy4403ArmfoU8++YQHH3yQf/7znyxcuJAXX3wxkf6FF15g2bJlNG/e3Ov1VOEYYFVgTz31VOJC5NTU1HxnErOyskhNTd3mdpdeeilNmjTh2muvLamsqhKpVS+VNf/6z9nytd8vo9Z+B+RLM/WtobTuehoAv2rTns2bfmbFihVMnjyZm266iQYNGvDoo49y7733OnZfKoSitgd16tShRo0anHHGGQD07NmzwGQYI0aMSAwP/MUtt9zC9OnTGTt2LFEU0aRJk2IqkUpCYeprgHHjxnHPPfcwevRoqlWrBuT2jGZkZNCwYUNSUlI4/fTTC3yGkpOTOfvssxMTqEgVhQFWBXbllVcmLkQ+/fTTeemll4iiiM8++4y9996bAw88sMA2t956K2vWrOHRRx8tsG5b47OloqrfMoMVSxexaNEisjdv4osP3qL5MSfmS7PPAQfx9ZQJAHy/cB7ZP29kv/3245NPPmHx4sUsXryYa6+9lv/5n/+hf//+pVEMqVwpansQQuDUU0/l448/BmD8+PG0aNEisf6rr75i1apVHHXUUYllW7ZsYeXKlQDMmDGDGTNmcMIJJxR/4VRsfqmvf/j2m+3W15mZmfz+979n9OjRiUksANq3b8+qVatYvnw5AB999BEtWrQgiiIWLFgA5Pamvv322zRr1qzkCiWVAPvuK4mTTz6ZMWPG0LhxY2rUqMELL7yQWJeens706dPJysrinnvuoVmzZrRt2xaA/v37c/HFF/PPf/6TM844g+9XruLVt0Zz5c238ofXJ+70dQdkFLwoWpVbckoK3W++j27durFyw2YO796H/Rs1Y+zTgzioRTotjjmRk6/7E2/e9QcmDv0LIUCPO5/wImgpJoVpDwDuv/9+zj33XK699lr222+/fOmGDx9O7969830vN2/eTKdOnYDcSQxeeeUVhwiWc7/U189f2YsoJydffd3it7+he/fu3Hjjjfz000/07NkTyJ11ePTo0SQnJ/Pggw9y3HHHEUUR7dq145JLLiGKIs4//3zWrl1LFEW0adOGp59+upRLKsXLmq+SCCHw1FNPbXPdL41p/fr1tzk2H3LPRGVlZTEoc0Wx5VGVR7Oju/LiVX3yfZ66Xv6f+67t37Apl70wZlubJgwcOLC4sidVaIVpDwB+9atfMWHChG2m29b3r3r16syYNXubkx+o/Gp2dFeaHd0137Kulw+ge94J1HHjxm13265duzJjxowCyydNmhRvJqUyxgBLUrmQnRMV+Yfbrmwjadf9MhNdYTnKofwqav1qfazKxABLUrlQ1B9u4I83SSouBtPS9jnJRSXmpBWSJElSvML2rrnZBn+NV0BFPfv0S/o/tt2vuLKkSiKKol3+/BVlG5VpJTFeyLarhNlOVCz3TVu+zeVFrZNvaFPHIYWqKHb6wXSIoKQKy2sEpN3jd0JxcUihKhMDrArEhlDKr6gN+g1t6hRp/37nVNH5o1ilxYmNVJ4ZYJVRu1JJ2BBKu6e4AzLwB4BKj589lSe7MrGRJ8lUVhT6Gqw777zzfaAkf5GnAstK8PXKMo/Ff3gs/sNj8R8ei/8oT8dixR133HFicb7AVm1XeTouu8NyVjyVpayWs2KpyOXcedsVRVGZ/Bs4cGBU2nkoK38eC4+Fx8Jj4bHwuFjOylfOylRWy1mx/ipLObf35zTtkiRJkhSTshxg3VnaGShDPBb/4bH4D4/Ff3gs/sNjsW2V5bhYzoqnspTVclYslaWc21SU+2BJkiRJknagLPdgSZIkSVK5YoAlSZIkSTExwJIkSZKkmJSZACuE0DOEMDuEkBNCOHwH6U4MIXwVQlgQQhhQknksKSGE2iGEsSGE+Xn/991Oui0hhOl5f6NLOp/FaWfvcwihWghhZN76ySGEBiWfy5JRiGNxQQhh+VafhYtLI5/FLYTwfAjh+xDCrO2sDyGEx/OO04wQQtuSzmNJKcSx6BxCWLPVZ+L2ks5jaStCm7I4hDAz7zhNLck8xqGytJ0VvV2sTG1eZWjTKkt7ZVu0fWUmwAJmAWcCE7aXIISQDDwFnAS0APqEEFqUTPZK1ABgfBRFTYDxec+3ZUMURel5f91LLnvFq5Dv80XAqiiKGgOPAPeXbC5LRhE+8yO3+iz8tUQzWXJeBHZ0Y7+TgCZ5f5cCT5dAnkrLi+z4WAB8stVn4k8lkKeyZqdtyla65B2n7QYoZVhlaTsrbLtYmdq8StSmvUjlaK9exLZom8pMgBVF0dwoir7aSbIjgAVRFC2MomgTMAI4rfhzV+JOA4bkPR4CnF6KeSkNhXmftz5GrwPHhRBCCeaxpFSWz/xORVE0AfhhB0lOA16Kcn0G7BNCOLBkcleyCnEsKr1CtinlXiVqOytyu1iZ2ryK8FncqcrSXtkWbV+ZCbAK6SBg6VbPs/KWVTT7R1H0HUDe/3rbSVc9hDA1hPBZCKEiNTaFeZ8TaaIoygbWAHVKJHclq7Cf+bPyhhm8HkI4uGSyVuZUlvqhsI4KIXwRQngvhNCytDNThkXAhyGEz0MIl5Z2ZopJRfhuVOR2sTK1ebZpuSrCd7KwKmVblFKSLxZCGAccsI1Vt0RRNKowu9jGsnJ5I68dHYsi7OaQKIqWhRAaAh+FEGZGUfR1PDksVYV5nyvMZ2EnClPOt4HhURT9HEK4jNyznMcWe87KnsrymSiMacCvoij6KYRwMvAWuUNRKpQY2hSAjnn1aD1gbAjhy7yzsmVGZWk7K3G7WJnaPNu0XBXl/dyZStEWbUuJBlhRFB2/m7vIArY+k1EfWLab+ywVOzoWIYR/hxAOjKLou7wu4++3s49lef8XhhA+BjKAst6QFEZh3udf0mSFEFKAvamY3dQ7PRZRFK3c6ulzlNOx+TGoMPXD7oqiaO1Wj8eEEP4cQqgbRdGK0sxX3GJoU7auR78PIbxJ7hCmMhVgVZa2sxK3i5WpzbNNy1UuvpOpYo0LAAAgAElEQVS7q7K0RdtS3oYI/hNoEkI4NIRQFegNlJtZgopgNHB+3uPzgQJnKEMI+4YQquU9rgt0BOaUWA6LV2He562PUQ/goyiKKuLZn50ei/8at90dmFuC+StLRgPn5c3O1AFY88uQosomhHDAL9dnhBCOILeuX7njrSqfEMKeIYS9fnkMnEDupBEVTUVoOytyu1iZ2jzbtFyVor2q1G1RFEVl4g84g9yI/mfg38AHectTgTFbpTsZmEfuGalbSjvfxXQs6pA7S9L8vP+185YfDvw17/GvgZnAF3n/LyrtfMd8DAq8z8CfgO55j6sDrwELgClAw9LOcykei/uA2Xmfhb8DzUo7z8V0HIYD3wGb8+qKi4DLgMvy1gdyZ6f6Ou87cXhp57kUj0X/rT4TnwG/Lu08l8Ix2mmbAjTMO0Zf5B2vctemVJa2s6K3i5WpzasMbVplaa9si7b/F/IOgCRJkiRpN5W3IYKSJEmSVGYZYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJf2XEMLHIYSmeY/rhBBmlXaeJEnanhBCmxDChBDCnBBCTgghCiHcWdr5kiqrlNLOgFQGNQbm5z1OA2bu6o5CCPtGUbQqllxJkvRfQgjVgZHAeVEUTQkh3AVUBwbuwr5ss6QY2IMlbSWE8Cvg2yiKcvIWpQEzdmOXj/zX/sM2XvPF3di/JKlyOx6YFkXRlLznM4DaURRFu7CvR/57ge2WVHT2YEn5pZM/oGoHjAwhpAD/C0TAN8BzwD1ANWAV8AzwMjAa6BBF0dkhhBOBZiGEW4GT89YNCSHcuNV29wLrQggHA7cDa4D3oygaV+wllSRVBK3IP9KiLTCtqO0W8AK5bdYNwCvA37DdknaJAZaUXxtyh1YQQmgCnAbcClwOjIqi6P/y1t0EDIuiaGoI4fW87d6KouixEMLQvH2tILeR+hpYG0XR49vYri0wDWgGbAIej6JoSUkVVpJU7q0EjgUIIRwGnAn8mqK3WyuAV6IoejKEcBIwwnZL2jUOEZTySweSQghfkHtmbi5wPrkNyqSt0rUEZoYQqgLryW2oPshb98uwjDTgi7x9jt3Odu2Bf0ZRNBZ4AngyhHBQMZVNklTxDAdq5k3I9CzQJ4qilRS93fqlzQLbLWm32IMl5ZcGZERR9OPWC0MIpwF/CSH8ANwHvEpuQ7Y+7/mNwLwQQl3gX3mbrQAuBg4B7s9b9t/bXQ88HkK4H0gGlgDfF1vpJEkVShRFPwGnbmPVWxSt3VoBXBxCWAE0Ab7K24/tllREYdeugZQqnhDCXsDnURQdVtp5kSRJUvlkgCVJkiRJMfEaLEmSJEmKiQGWJEmSJMXEAEuSJEmSYlKUWQS9WEuSFKdQAq9h2yVJitNO2y57sCRJkiQpJgZYkiRJkhQTAyxJkiRJiokBliRJkiTFxABLkiRJkmJigCWVAdk5RZvorKjpJVU81huSVDaFKCp0hWvNLBWjQZkrCp12QEbdYsyJVGKcpn03WW9IUolzmnZJkiRJKikGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS4qZM3VJkiRVXimlnQGpoklJCkWa2Quc3UuSJKmisAdLkiRJkmJigCVJkiRJMTHAkiRJkqSYGGBJkiRJUkwMsCRJkiQpJgZYUimb8/F7pKWl8XjvzjzZ93gWZ362zXTZmzfxt7uu48HTj6RZs2a88cYbACxZsoQuXbqQkZFBWloaY8aMAWDx4sXssccepKenk56ezmWXXVZiZZJUPkVRxNVXX03jxo1JS0tj2rRp20zXuXNnmjZtmqhfvv/+ewB+/vlnzj77bBo3bsyRRx7J4sWLARg6dGgibXp6OklJSUyfPr2kiiVJJcpp2qVS1uiITgy59nfcP30l382bzfABF3Pd3/5RIN3f//oINWvX5Ya3JnNTm9r88MMPANx999306tWLyy+/nDlz5nDyyScnftQ0atTIHzGSCu29995j/vz5zJ8/n8mTJ3P55ZczefLkbaYdOnQohx9+eL5lgwcPZt9992XBggWMGDGCm2++mZEjR9K3b1/69u0LwMyZMznttNNIT08v9vJIUmmwB0uVzuLFi2nWrBkXX3wxrVq1om/fvowbN46OHTvSpEkTpkyZwrp167jwwgtp3749GRkZjBo1KrFtp06daNu2LW3btuXTTz8F4OOPP6Zz58706NGDZs2aMeKWy4iiwt1wuFqNmoQQANi0YT0Qtpnu89HD6HzhNQAkJSVRt27uvbNCCKxduxaANWvWkJqausvHRlLJ2t366C8X/pYnzjmWJ845lm++mALAwqmTePaS0xh6Yz8ePvOoItVHo0aN4rzzziOEQIcOHVi9ejXfffddocszatQozj//fAB69OjB+PHjC7z28OHD6dOnT6H3KUnljT1YqpQWLFjAa6+9xrPPPkv79u0ZNmwYEydOZPTo0dx77720aNGCY489lueff57Vq1dzxBFHcPzxx1OvXj3Gjh1L9erVmT9/Pn369GHq1KkAZGZmMnv2bFJTU2mQfiTfTJ9Mg4wOvPPgrSycOrFAHtK6nUHnfrkB05tvvsnD193ETz+s4PzHhhVIu+HHNQB8+OdBLPp8Ep+3asqTTz7J/vvvz8CBAznhhBN44oknWLduHePGjUtst2jRIjIyMqhVqxZ33303nTp1Ko7DKWk37E59dOHTr1OlWnVWLPmaEX/8Pf2H5n7/l301kz+8NpG99juAZ/qdwjfTJ5Odfgo3Xn8df//73wvkoXfv3gwYMIBvv/2Wgw8+OLH8oPr1+fbbbznwwAMLbNOvXz+Sk5M566yzuPXWWwkh5Ns+JSWFvffem5UrVyZOCAGMHDkyESRKUkVkgKVK6dBDD6V169YAtGzZkuOOO44QAq1bt2bx4sVkZWUxevRoHnzwQQA2btzIkiVLSE1NpX///kyfPp3k5GTmzZuX2OcRRxxB/fr1AUht2opVy5bSIKMDv73h7p3m54wzzuCrBp1Y9PmnjH16EBc/80a+9TnZ2az59zIapB/Bb6+/i6p/f4kbbriBl19+meHDh3PBBRdw/fXX849//INzzz2XWbNmceCBB7JkyRLq1KnD559/zumnn87s2bOpVatWXIdRUgx2pz76213X8d28WSQlJbFiycLEPg9u2Za998/tzf6lPkpJCux/3i30Pu+WbeZjUOYKFqz5mVfmrWbiniuA3P70X3rYtzZ06FAOOuggfvzxR8466yxefvllzjvvvG32lG29/eTJk6lRowatWrXatYMlSeWAAZYqpWrVqiUeJyUlJZ4nJSWRnZ1NcnIyb7zxBk2bNs233cCBA9l///354osvyMnJoXr16tvcZ0hKImdLNkCherB+cWi7X/PDHVexbtVK9ty3TmJ5jX1qU6V6DVp0OQWAnj17MnjwYCD3mof3338fgKOOOoqNGzeyYsUK6tWrl8hTu3btaNSoEfPmzStwzYSk0rU79dFedfaj54iPiXJyuP2o+ol1yVWqJh4XpT7au14qq/+9LLE8Kytrm8OODzroIAD22msvzjnnHKZMmcJ5551H/fr1Wbp0KfXr1yc7O5s1a9ZQu3btxHYjRoxweKCkCs8AS9qGbt268cQTT/DEE08QQiAzM5OMjAzWrFlD/fr1SUpKYsiQIWzZsmWn+9pZD9aKJQuJ0nODqW/nfsGWzZuosU/tfGlCCDT/zQksmjqJRkd0Yvz48bRo0QKAQw45hPHjx3PBBRcwd+5cNm7cyH777cfy5cupXbs2ycnJLFy4kPnz59OwYcNdPCKSSsuO6qO96u5PUlISU98eQU4M9VHzY7rxj5GDadPtDJbO/Jy99967wPDA7OxsVq9eTd26ddm8eTPvvPMOxx9/PADdu3dnyJAhHHXUUbz++usce+yxiR6snJwcXnvtNSZMmLCLR0KSygcDLGkbbrvtNq699lrS0tKIoogGDRrwzjvvcMUVV3DWWWfx2muv0aVLF/bcc8/dfq3ZH71Dq//px6rsQEq1Pegz6LnED5LHe3fm6hEfA3Di1bfz6m1X8M6Dt9LikAN44YUXAHjooYe45JJLeOSRRwgh8OKLLxJCYMKECdx+++2kpKSQnJzMM888k+9MsqTyYUf10W9OPp2ZY0fTsH1Hqu5RY7dfq+nRXflq4jgePO0IqlTfg3eGv5RYl56ezvTp0/n555/p1q0bmzdvZsuWLRx//PFccsklAFx00UWce+65NG7cmNq1azNixIjE9hMmTKB+/fqe6JFU4YXCziwEFDqhVNkNylxRpPQDMuoWaZsBGXV3nkgq+7Y9ZWa8KnTbVdR6Y1fqJklSPjttu5ymXZIkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsqRzKzinazNNFTS9JkqRd442GpXIoJSl43yxJkqQyyB4sSZIkSYqJAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsSZIkSYqJAZa0E9k5UWlnQZIkSeVESmlnQCrrUpICgzJXFDr9gIy6xZgbSZIklWX2YEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWJEmSJMXEAEuSJEmSYmKAJUmSJEkxMcCSJEmSpJgYYEmSVAZk50SlnQVJUgxSSjsDkiQJUpICgzJXFDr9gIy6xZibXNk5ESlJodjSS1JFZIAlSZK2qSwGfZJU1jlEUJIkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS6oEsnOiYk0vSZKkXCmlnQFJxS8lKTAoc0Wh0w/IqFuMuZEkSaq47MGSJEmSpJgYYEmSJElSTAywJEmSJCkmBliSJEmSFBMDLEmSJEmKiQGWKh2nIJckSVJxcZp2VTpOWS5JkqTiYg+WJEmSJMXEAEuSJMWiqEOwHbItqSJyiKAkSYqFQ7AlyR4sSZIkSYqNAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLAkSZIkKSYGWJIkSZIUEwMsSZIkSYqJAZYkSZIkxcQAS5IkSZJiYoAlSZIkSTExwJIkSZKkmBhgSZIkSVJMDLCkHZg07C+0atWKR3oczcShzySWj3vmf7mvW2se792Zx3t35suJYwFYPH0yaWlpPPm7rqxYshCADT+u4fkrehJF0TZf49lLTmPq1KmJ56uWLeHRnp0AWDh1EgN/05DH+3Th4TN/zbi/PADAxx9/nFj+0Bkd+MtFpzJ3wofFcgwklV2Thv2FR3t2KlId9VivY4pcR2XNmZ54vqM66s477yyw3DpKUmWTUtoZkMqqfy2Yyz/ffIWFMz7nkdlreaH/2TTr1JW6hzQCoGPfy/jNeVfm22biy0/z5htvMGjcF0x+/UVOue5PfPTcQ3S+6FpCCLuUjwbpHbjg8WFs2rCOx3t3oXmnE+Cg5MRygGVfzeTl686nSrXqND7yN7tXcEnlwqxZs/jnm69wxUsfkFylaqHrqL4PPs+qZUuLpY56+fzjObHx0fmWg3WUpMrFHixpO5YvmsfBrdtRo0YNklNSOLTdr5n90ZgdbpOUksKGDRvYvHEDySkprFy6iLXff0fDdh13Oz9V99iTg5q3YWXWogLrUpu25rhLrucfrw7e7deRVD7MnTuXg1u3o+oeRaujNm/cWGx1VLt27ayjJFV69mBJ27F/o+Z88NS9rFy5kk0b1vPVxHHUb9Emsf4fIweT+c6rHNSiDadc9yf2qLUPnS+8hksvvZTlW1LoddefGfPIHXS9YsBOX6tv3778SBUAtmzeTEgqeO5j3eofWDLzc4695Hogu8D61OZpTHjpqV0vsKRypVWrViyaNoB1q3+gSrXqO62joC6dL7yGN+++nirVqxepjhp5y2VUqVYd2HEd9dlnn3FGj/6sW7WywHrrKEmVhQGWtB31Gh7GMRdcRdeuXVlFNQ48rCVJyblfmSN7XpAb6ITA2D/fx7sP306PgY+T2rQ1n332GYMyV7Do80/Za78DiCIYdvPFJKdU4eTr7mSvOvUKvNbQoUMZl9wAyL2+Ycg1fRPrFk//jMf7dCGEJDr3u5r9GzWDNbMKZng7109IqpiaN2/OMRdcxfNX9KDqHnvutI6685hhpDZtzRUvvQ9QpDrq7HueoX6LdGDHddRtAwawulEzFk6dVDDD26ijsnMiUpIKPzSxqOklqTQYYEk70P703/HGndcyKHMFHzxxN7X2TwXI9wPkiDPPzfdjAyCKIj4a/AjnDHqOUfcP4PjLbmbVsiV8Ovw5uvW/pUh52Po6hh1Z9uVM6h16WJH2Lal8a3/672h/+u8AykQddVlGXQZlrthmum3VUSlJYbvpt2VARt0i5U2SSoPXYEk78NMPywFY/V0Ws//+LuknngnA2uX/SqSZ/dGY3F6lrUx7ewTNjj6ePWrtw+aNGwhJSYSkJDZv3FAs+fxu3mw++uvDdOh1YSz7y84pem/YrmwjafdU1jpKksoye7BUrhX3cJGhN/Tjbz+vZVV2Et1vvp89au0DwHuP/Ynv5s0iENg39WBOv+XBxDbr169n2jsjufCp1wA4uu9lDL2xH8kpVeh937Ox5e2XYTmbN26g5r51OfXGe2ObnauoZ5XBM8tSaRh6Qz/Wr1lFUkqVQtdRmzaU/zpKksqysL37XmyDp6dVJu1KIFDUISnFmb4s5mlXyyAVUUlcTFOu2q6y+L0uS3mynpFUBuy07XKIoCRJkiTFxABLkiRJkmJigCVJkiRJMTHAkiRJkqSYGGBJkiRJUkwMsCRJkiQpJgZYkiRJkhQTAyxJkiRJiokBlgpt48aNHHHEEbRp04aWLVtyxx13FEjzzDPP0Lp1a9LT0zn66KOZM2cOACtXrqRLly7UrFmT/v3759vmlltu4eCDD6ZmzZolUg6VnC1btpCRkcFvf/vbAut+/vlnzj77bBo3bsyRRx7J4sWLAdi8eTPnn38+rVu3pnnz5tx3330ALF26lC5dutC8eXNatmzJY489VpJFUQXw/vvv07RpUxo3bsygQYMKrH/44Ydp0aIFaWlpHHfccXzzzTcAfPPNN7Rr14709HRatmzJM888k9imc+fONG3alPT0dNLT0/n+++9LrDyVUXZO0e8b/e6Y93b4vi9ZsoQuXbqQkZFBWloaY8aMSaybMWMGRx11FC1btqR169Zs3LgRgOHDh9O6dWvS0tI48cQTWbGiaDdwllSxpZR2BlR+VKtWjY8++oiaNWuyefNmjj76aE466SQ6dOiQSHPOOedw2WWXATB69Giuu+463n//fapXr85dd93FrFmzmDVrVr79nnrqqfTv358mTZqUaHlU/B577DGaN2/O2rVrC6wbPHgw++67LwsWLGDEiBHcfPPNjBw5ktdee42ff/6ZmTNnsn79elq0aEGfPn2oVq0aDz30EG3btuXHH3+kXbt2dO3alRYtWpRCyVTebNmyhSuvvJKxY8dSv3592rdvT/fu3fN9fjIyMpg6dSo1atTg6aef5qabbmLkyJEceOCBfPrpp1SrVo2ffvqJVq1a0b17d1JTUwEYOnQohx9+eGkVrVJJSQoMyix8MJOzZQuDr+q/w/f97rvvplevXlx++eXMmTOHk08+mcWLF5Odnc3vfvc7Xn75Zdq0acPKlSupUqUK2dnZXHPNNcyZM4e6dety00038eSTTzJw4MBiKLGk8sgeLBVaCCHRy7R582Y2b95MCCFfmlq1aiUer1u3LrF+zz335Oijj6Z69eoF9tuhQwcOPPDAYsy5SkNWVhbvvvsuF1988TbXjxo1ivPPPx+AHj16MH78eKIoIoTAunXryM7OZsOGDVStWpVatWpx4IEH0rZtWwD22msvmjdvzrffflti5VH5NmXKFBo3bkzDhg2pWrUqvXv3ZtSoUfnSdOnShRo1agC59VJWVhYAVatWpVq1akBuz2tOTs5OX29XeloUv6Wzpu30fQ8hJE4CrVmzJhE4f/jhh6SlpdGmTRsA6tSpQ3JyMlEUEUUR69atI4oi1q5dm9hGksAeLBXRli1baNeuHQsWLODKK6/kyCOPLJDmqaee4uGHH2bTpk189NFHpZBLlQXXXnst//u//8uPP/64zfXffvstBx98MAApKSnsvfferFy5kh49ejBq1CgOPPBA1q9fzyOPPELt2rXzbbt48WIyMzO3+fmTtmXrzxtA/fr1mTx58nbTDx48mJNOOinxfOnSpZxyyiksWLCABx54IN8P6n79+pGcnMxZZ53FrbfeSgihyD0tAAMy6hYpvXZu7fLvdvq+Dxw4kBNOOIEnnniCdevWMW7cOADmzZtHCIFu3bqxfPlyevfuzU033USVKlV4+umnad26NXvuuSdNmjThqaeeKtFySSrb7MFSkSQnJzN9+nSysrKYMmVKgeF+AFdeeSVff/01999/P3fffXcp5FKlYesz9u+88w716tWjXbt2200fRQXP8IcQmDJlCsnJySxbtoxFixbx0EMPsXDhwkSan376ibPOOotHH300X4+ptCPb+7xtyyuvvMLUqVO58cYbE8sOPvhgZsyYwYIFCxgyZAj//ve/gdzhgTNnzuSTTz7hk08+4eWXXy6eAmjXFOJ9Hz58OBdccAFZWVmMGTOGc889l5ycHLKzs5k4cSJDhw5l4sSJvPnmm4wfP57Nmzfz9NNPk5mZybJly0hLS0tcKypJYIClXbTPPvvQuXNn3n///e2m6d27N2+99VaR9uuwmvLrlzP2gzJX8ODfxjL0jbfYN/UQTu3Riw/Hf0TGyT0T6wdlrqB+/fosXboUgOzsbNasWUPt2rUZNmwYJ554IlWqVKFevXp07NiRqVOnArlDU8866yz69u3LmWeeWZrFVTmz9ecNcoewbmtY17hx47jnnnsYPXp0Yljg1lJTU2nZsiWffPIJAAcddBCQO2z1nHPOYcqUKcVUAu2KWvVSd/q+Dx48mF69egFw1FFHsXHjRlasyK2jjjnmGOrWrUuNGjU4+eSTmTZtGtOnTwegUaNGhBDo1asXn376ackVSlKZZ4ClQlu+fDmrV68GYMOGDYwbN45mzZrlSzN//vzE43fffbfIE1ds/SO9MH8qm0686jb++P4Mbn53Gn3ue46Ghx/N2fc8nS9N9+7dGTJkCACvv/46xx57LCEEDjnkED766KPENQ6fffYZzZo1I4oiLrroIpo3b851111XGsVSOda+fXvmz5/PokWL2LRpEyNGjKB79+750mRmZvL73/+e0aNHU69evcTyrKwsNmzYAMCqVauYNGkSTZs2JTs7OzF73ObNm3nnnXdo1apVyRVKO1W/ZcZO3/dDDjmE8ePHAzB37lw2btzIfvvtR7du3ZgxYwbr168nOzub//u//6NFixYcdNBBzJkzh+XLlwMwduxYmjdvXuJlk1R2eQ2WCu27777j/PPPZ8uWLeTk5NCrVy9++9vfcvvtt3P44YfTvXt3nnzyScaNG0eVKlXYd999Ez+gARo0aMDatWvZtGkTb731Fh9++CEtWrTgpptuYtiwYaxfv5769evT7JRzOP6ym0qxpCouY58exEEt0mlxzIlcdNFFnHvuuTRu3JjatWszYsQIIHeIab9+/WjVqhVRFNGvXz/S0tKYOHEiL7/8cuI2AAD33nsvJ598cmkWSeVESkoKTz75JN26dWPLli1ceOGFtGzZMl/9deONN/LTTz/Rs2dPIPeH9+jRo5k7dy7XX///7N15mFTFufjx7zsMi4gLihgRI3pRZN8EMW64XY0a1LhhNMYlet1+7jF6E6PmRsXE664xiXGNgvuFGDVijFFMBBVwj2AUATFBVBQUhGHq90c3k2GYYbqHnpnume/neeaZnnPq9FRVn646b1ed6nOJCFJKnHfeefTv358vvviCffbZh+XLl7NixQr22msvTjzxxGYuqaprk33dh43ci1RZyfajjuT3yzbljJPOq2qLtvv+j7nkf87mwsuvIgLuuOMOIoLOnTtzzjnnMGzYMCKC/fbbj/333x+Aiy++mF133ZW2bduy5ZZbcscddzRvQSUVlahtXnodnLulJpHPyNQFg7s06EbyxvwfLSFPTVUGtXq13wRVWM3Wd5V6O1CMebJtklQE6u27nCIoSZIkSQVigCWpWeS7oIkLoEhqCrZNktaW92CpUVVUJsrLmmIWkEpNvt8T5LQdSU3BtknS2vIeLDW6fDqqC4ds0og5UVO7YupHde7zvgjRAu/Bquu7tdR8UkqNfl+YpFbFe7AkSZIaS0OmCDqtUGrZnCKonDndT82pIeef56ykxpbvlEKA8wZunFd62zKptBhgtWL5NtgN6UScOqFC8SJGUkuRb3tmWyaVFgOsJpJvY7e8MtE2z8Yx32OK8Ubexp4r35BjivV7WoopfUOPaWyNfRHTkPepFz6lKdfXrfp9zaXeDhRjnoqxnWkKjd2WQf7tWbGlB9tXFY+cF7m49NJLnwCKdTiiGzCvuTNRxKyf+llHa2b9rJn1s2Z11c+Ciy++eN/G/MdF3nc1N8/b/Fln+bPO8med5a8p66z+viulVPI/l1xySWruPBTzj/VjHVk/1o/144+vi3VmnZXGj3VW+nXmKoKSJEmSVCAtJcC6tLkzUOSsn/pZR2tm/ayZ9bNm1k9x8nXJn3WWP+ssf9ZZ/oqqzvL5om/JhnsAACAASURBVGFJkiRJ0hq0lBEsSZIkSWp2BliSJEmSVCAlE2BFxBYR8eeIeCsi3oiIM2tJMzIiPouI6dmfnzRHXptDRHSIiCkR8Uq2flabixoR7SPivoh4JyImR0SPps9p88ixfo6NiI+qnT/fb468NqeIaBMR0yLi0Vr2tdrzZ6V66sfzJ2JWRLyWLf9LteyPiLg+ew69GhFDmiOfrUlE7BsRb2fr/IJa9rf687amiLgtIuZHxOt17Pc8riGHOmu112d1yfG61nOtmlKKBUrpi4YrgHNTSlMjYj3g5YiYmFJ6s0a651JKBzRD/prbV8AeKaXFEdEWmBQRj6eUXqiW5gTg05RSz4gYDVwJHNEcmW0GudQPwH0ppdObIX/F4kzgLWD9Wva15vNnpTXVD3j+AOyeUqrrG1G/CWyT/dkB+GX2txpBRLQBbgL2BuYCL0bEhFr6Tc/bVd0B3AjcVcd+z+PV3cGa6wxa7/VZXXK5rvVcW1XJxAIlM4KVUvowpTQ1+3gRmYuczZs3V8UjZSzO/tk2+1NzBZMDgTuzjx8E9oyIVvGV5znWT6sWEd2B/YFb60jSas8fyKl+VL8Dgbuy78cXgA0jYrPmzlQLNhx4J6X0bkppGTCOzGugNUgpPQt8soYknsc15FBnqiHH61rPtWpKKRYomQCruuzUpMHA5Fp275idBvZ4RPRt0ow1s+z0penAfGBiSqlm/WwOzAFIKVUAnwEbN20um08O9QNwSHYY/sGI2KKJs9jcrgXOByrr2N+qzx/qrx9o3ecPZD60eDIiXo6Ik2rZX3UOZc2lSDvHFiLX+m7t522+PI8bptVen9VnDde1nmt1KPZYoOQCrIjoBDwEnJVS+rzG7qnAlimlgcANwP81df6aU0ppRUppENAdGB4R/WokqW20odWM4uRQP78HeqSUBgBP8e/RmhYvIg4A5qeUXl5Tslq2tYrzJ8f6abXnTzU7pZSGkJnWclpE7Fpjf6s9h5pJLvXteZs/z+P8terrszWp57rWc60WpRALlFSAlb135iHgnpTSwzX3p5Q+XzkNLKX0GNA2Iro0cTabXUppIfAMsG+NXXOBLQAiohzYgFY4pF9X/aSUPk4pfZX98zfA0CbOWnPaCRgVEbPITCPaIyJ+VyNNaz5/6q2fVn7+AJBSmpf9PR94hMwUteqqzqGs7sC8psldq1RvfXveNojncZ68Pqtdfde1eK6tplRigZIJsLL3evwWeCuldHUdab628p6QiBhOpnwfN10um09EbBIRG2YfrwPsBfy9RrIJwPeyjw8Fnk6t5Jumc6mfGvOaR5GZ29sqpJQuTCl1Tyn1AEaTOTeOrpGs1Z4/udRPaz5/ACJi3exNx0TEusB/AjVXFJsAHJNdGWsE8FlK6cMmzmpr8iKwTURsFRHtyJy7E6onaO3nbQN5HuepNV+f1SWX61o811ZRSrFAKa0iuBPwXeC17H00AP8NfB0gpXQLmYu+UyKiAlgCjG4tF4DAZsCd2VWjyoD7U0qPRsRPgZdSShPInJR3R8Q7ZEYeRjdfdptcLvVzRkSMIrNKzSfAsc2W2yLh+bNmnj+r2BR4JNuvlQP3ppSeiIiToaqNfgzYD3gH+BI4rpny2iqklCoi4nTgj0Ab4LaU0huet2sWEWOBkUCXiJgLXExmYSTP4zrkUGet+fqsLrlc13qurapkYoHw/JYkSZKkwiiZKYKSJEmSVOwMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsKQaIuKZiOiVfbxxRLze3HmSJKkuETEwIp6NiDcjojIiUkRc2tz5klqr8ubOgFSEegIzs48HAK/lemBEdE4pfZrPP2vIMZIkAUREB+A+4JiU0pSI+B+gA3BJjsfbb0kF5giWVE1EbAl8kFKqzG4aALyax1Nck32e/8n3GEmSGmAvYGpKaUr271eBjVJKKcfj7bekAnMES1rVIFYNqIYC90VEV+AJ4I9AL+BQ4CJgI2BhSuniiNgX2C4ifgyUR0Rb4FKgI9AupXRqRHQH7gYmACOA27PHnEfmE8ifAJ8BT6SUnmr84kqSSlw/Vp1pMQSYar8lNR9HsKRVDSQztYKI2AY4kEzHNQwYm1K6EJgPbAG0BRaS6XAAFgC/A14GpgMnAetk03Sq9vz/l1K6BqhYeUxK6SpgO2AZcL2dlCQpRx+TmW1BRGwLfBsYh/2W1GwcwZJWNQhYEhGvkBnJegv4HtAG+Gs2zQbAxcCZwCZkOi3IdHCvADsDDwM/AE5LKX1V7fkHZvcBpGrHkFKaGBFzgBsj4pSU0geNUkJJUksyFhiVXZBpAXBkSunjiBiG/ZbULAywpFUNAAanlBZV3xgRY4GNI+Jg4FYygdh5wMbAtGyyBcD3ge7AlcB44I5s5/N0SukJMgtozIiILsA/Vx4TEQuAY8kEcrPJfNooSdIapZQWA9+qZVcv7LekZhG53wMptWwRsR7wckpp21r23Z1S+m4zZEuSpLzZb0nNxwBLkiRJkgrERS4kSZIkqUAMsCRJkiSpQPJZ5MK5hJKkQoom+B/2XZKkQqq373IES5IkSZIKxABLkiRJkgrEAEuSJEmSCsQAS5IkSZIKxABLkiRJkgrEAEuqR0VlfouQ5ZtekloK20tJgkgp58bNVlCt1phpC3JOe8HgLo2YE6lFcZn2Fsj2UlIL5zLtkiSpODniJaklyueLhiVJkgqmvCwc8ZLU4jiCJRVYQz5h9VNZSZKklsERLKnA8v1EFvxUVpIkqaVwBEuSJEmSCsQAS5IkSZIKxABLkiRJkgrEAEuSJEmSCsQAS5IkSZIKxABLkiRJkgrEAEuSJEmSCsQAS5IkSZIKxABLWgtvPvM41x2+G9ePHsmNR+3FrGkv1JruttMO57ojRnLNoTvzyGXnUbliRdW+v477Db169aJv376cf/75VduvuOIKevbsSa9evfjjH//Y6GWRpGKVUmLCzy+kZ8+eDBgwgKlTp66WZtGiRQwaNKjqp0uXLpx11lkA3HLLLfTv359Bgwax88478+abbwIwZcqUqvQDBw7kkUceadJySWqZIqWUa9qcE0rFrKIyUV4WeR0zZtqCWrd/9eVi2q2zLhHBhzPeYOwF32f+ezNXS7908SI6dFqPlBL3/OA4+u99IAP3OZh/vDiJP//2Gt549knat2/P/Pnz6dq1K2+++SZHHnkkU6ZMYd68eey1117MmDGDNm3aNLjcUhHK743YMPZdTayu9rI2FwzuklP6v0+ayN/G3crfJz3F5MmTOfPMM5k8efIajxk6dCjXXHMNu+66K59//jnrr78+ABMmTODmm2/miSee4Msvv6Rdu3aUl5fz4YcfMnDgQObNm0d5eXnOZZDU6tTbdzmCpZIya9YstttuO77//e/Tr18/jjrqKJ566il22mknttlmG6ZMmcIXX3zB8ccfz7Bhwxg8eDDjx4+vOnaXXXZh+PZD2bz3AE654zHGTFvASb8Zz9bb70T/vb5F1622YdB+h3LF1I8YM21BvR1/+46diMi8z5Yt+ZK63nMdOq0HQGVFBSuWL69KNfnB2xl53Bm0b98egK5duwIwfvx4Ro8eTfv27dlqq63o2bMnU6ZMWcvak6TcrE1b++m82fzq+AO44Tt7cMN39uD9VzJt17svPc+vTzyQe35wHFd/e0fG/ehkcv2Q961nnmDwAUcQEYwYMYKFCxfy4Ycf1pl+5syZzJ8/n1122QWgKrgC+OKLL6ra7Y4dO1YFU0uXLq3aLklrw49oVHLeeecdHnjgAX79618zbNgw7r33XiZNmsSECRO4/PLL6dOnD3vssQe33XYbCxcuZPjw4ey111507dqViRMn0qFDB84bP5lxF/4Xp9/zFADz3n6Nsx+YxHqbfI1bjtuf96dPpsfgETx61Y8Z98YLzF9SsUoeBuxzMCOPOxOAN57+A3+88Wcs/mQB37vu3jrzfduphzHnjWn02mlP+u01CoAF7/+D96a+wA47/JwOHTpw1VVXMWzYMD744ANGjBhRdWz37t354IMPCl2VklSnhra163buwvG/fJC27TuwYPY/1tjWPv/887Dudjx61Y9596VJq+VhZVv72fwP2XDTblXbV7aJm222Wa15Hzt2LEccccQqAdNNN93E1VdfzbJly3j66aertk+ePJnjjz+e999/n7vvvtvRK0lrzVZEJWerrbaif//+APTt25c999yTiKB///7MmjWLuXPnMmHCBK666iog86nk7Nmz6datG6effjrTp0/n42WJBbPfrXrOLfoOYYNs592tVz8+nTeHHoNHcMB5P6t3CkvfPfan7x77897Lf2XiL8fAsd+sNd3xNz/A8q+Wct+PTuYfLz7HNiNGUrliBUsWLeSFF17gxRdf5PDDD+fdd9+t9VNdP1mV1JQa2tZWVrRn/JUX8OGM1ykrK1tjWztr1izoux0HnPezenKTX5s4btw47r777lW2nXbaaZx22mnce++9/OxnP+POO+8EYIcdduCNN97grbfe4nvf+x7f/OY36dChQ33VI0l1MsBSyVk5nQ6grKys6u+ysjIqKipo06YNDz30EL169VrluEsuuYRNN92UV155hStens9Pduxeta9N23ZVj6OsjMoVmRGrXEawVtpq6Df45OL/x4IFdQdjbdt3oPdu+/LmM4+zzYiRrN91M/rtcQARwfDhwzMXIwsW0L17d+bMmVN13Ny5c+nWrVudzytJhdbQtnbSteez3sabcNi4Z0iVlWtsaysq/t3WrmkEa4Ou3Vj4r3lV29fUJr7yyitUVFQwdOjQWvePHj2aU045ZbXtvXv3Zt111+X1119n++23r/VYScqFAZZanH322YcbbriBG264gYhg2rRpDB48mM8++4zu3btTVlbGtD/cv8pKfnWpbwRrwex32XiLrYgIPnjrFVYsX8bGG28Mcz6uSvPVl4v56ovFrL/J11hRUcHbk56ix+DM9L++u+/HP158Dk48kBkzZrBs2TK6dOnCqFGj+M53vsM555zDvHnzmDlzJsOHDy9MBUlSAdTV1i5d/DkbbNqNsrIyXvr9uJzb2jXpvds+/O2+35IuOJHJkyezwQYbrHF64JFHHrnKtpkzZ7LNNtsA8Ic//KHq8XvvvccWW2xBeXk577//Pm+//TY9evTIofSSVDcDLLU4F110EWeddRYDBgwgpUSPHj149NFHOfXUUznkkEN44IEH6NBnB9qt03Gt/9cbTz/K1Efvp015OeXt1+HIMb+pmrZy/eiRnDHuGZYt+ZK7zv4uK5Yto7JyBf8xbGd2OPRYAIYe+B0euuRM+vXrR7t27bjzzjuJCPr27cvhhx9Onz59KC8v56abbnIFQUlFpa62dsThx3HPecfz2sQJbD1sp4K0tb123pu3Jz1Fz5496dixI7fffnvVvkGDBjF9+vSqv++//34ee+yxVY6/8cYbeeqpp2jbti2dO3eumh44adIkxowZQ9u2bSkrK+Pmm2+mS5cua51fSa2by7SrVWqMZYQbmn7lMVIr5DLtLVBjt6+S1Mxcpl2SJEmSmooBliRJKgkVlfkPSDbkGElaG96DJUmSalVRmSgvK56viCgvC6dgSyp6BlhSEcj3IqbYLnoktUz5BjQGM5JkgCUVBS9iJEmSWgbvwZIkSZKkAjHAkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJkqQCMcCSJEmSpAIxwFJJq6hMzZ0FSZIkqUp5c2dAWhvlZcGYaQvyOuaCwV0aKTeSJElq7RzBkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJkqQCMcCSJEktVr5fSO8X2EtaW37RsCRJarHy/UJ6v4xe0tpyBEuSJEmSCsQAS5IkSZIKxABLkiRJkgrEAEuSJEmSCsQAS5IkSZIKxABLkiRJkgrEAEsqQX6viyRJUnHye7CkEuT3ukiSJBUnR7AkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLBUVCoqU3NnQZIkSWqw8ubOgFRdeVkwZtqCnNNfMLhLI+ZGkiRJyo8jWJIkSZJUIAZYkiS1Ak7BlqSm4RRBSZJagXynYIPTsCWpIRzBkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJysp3MRAXD5FUk4tcSJIkZfl9jJLWliNYkiRJklQgBliSJEmSVCAGWJIkSZJUIAZYUivgTduSJElNw0UupFbAm7YlSZKahiNYkiRJklQgBliSJEmSVCAGWJIkSZJUIAZYkiRJklQgBliSJEmSVCAGWJIkSZJUIAZYkiRJklQgBliSJEmSVCAGWGpUFZWpubMgSZIkNZny5s6AWrbysmDMtAU5p79gcJdGzI0kSZLUuBzBkiRJkqQCMcCSJElqoIZMhXf6vNSyOUVQkiSpgfKdCg9Oh5daOkewJEmSJKlADLDUIj1/76+49rBduObQnZl0zy1V25+65edsvvnmXD96JNePHsnfJ00EYNb0yVx3+G7cePTeLJj9LgBLFn3GbaceRkq1T+X49YkHMvfN6VV/fzpvNtcetgsA7770PJfsujXXH7k7V3/7Gzz1q1+stv1/Dx7Br074Fm89+2Sj1IEkNaXn7/0V/fr1q7XdvWKf/rW2uwMGDMir3R05cmTe7e4zzzxjuyupSTlFUC3OP995ixcf+R2n3vVH2rRtx+2nH8F2u+xNl6//BwBnn302FXseu8oxk+7+JUdddRufzpvD5AfvYP9zfsrTv/lfRp5wFhHRoHz0GDSCY6+/l2VLvuD60bvTe5f/XGU7wLy3X+Puc77Hn/ptChsNbHCZJak5rWx33331Za554/PV2t2djjqZXY85bZVjJt39Sx556CHGPPVK47a7m7eptd1t274DPXfYde0KLkm1cARLLc5H781gi/5DabdOR9qUl7PV0G/wxtOPrfGYsvJyli9dyvKlS2hTXs7Hc97j8/kfsvXQndY6P+3WWZfNew/k47nvrbavW6/+7Hniudx4441r/X8kqbmsbHc7dsyv3V2yZEmztbt/u/+3a/1/JKk2jmCpxdn0P3rzx5su54uFn9C2fQfenvQU3fv8e3Toxhtv5Ktf3c7mfQay/zk/ZZ31N2Tk8WfyyM/OpW2HDhz+Pzfz2DUXs/epF9T7v+770cm0bd8BgBXLlxNlq39m8cXCT5j92svsceK5fPHpx6vt79Z7ABPvv4Ud1qLMktScVra7H3/8McuWfLlau/u3+37LtEfvX63dPemkk/hoRXmjtrtQsdr+br0H8OxdNzW8wJK0BgZYanG6br0tux37/7jt1ENpt866bLZtX8raZE71HQ47lj/efAVXTv+YiTdfwR+u/gmHXnI93Xr159S7ngDgvZf/ynqbfI2U4N4ffp/Xu3Ziy2P/m/U27rra/zrislvo3mcQkLkX4M4zj6raN2v6C1x/5O5ElDHyuDPY9D+2492Xnl89w3Xca9CcKioT5WX5TdFpyDGSWoaV7e7ee+/Np7Rfrd3d48RzIWK1dveFF15gzLQFq7W7bcrbst85lxak3eWz11fPcBG2u5JaDgMstUjDDjqaYQcdDcAfb/gZ62/aDYD1Nu5KmzZtKCsrY/i3v7tKxwyQUuLp317Dd8b8hvFXXsBeJ/+Q/+z4GVeO/Q37nP6jvPJQfc7/msz7+2v07t07r+dubC47LClfww46mocuPYsx0xas1u6ulGu7++m82fy1kdvdrlttm9dzS1KuvAdLLdLiTz4CYOGHc3njz39g0L7fBuDzj/5ZleaNpx/LfLpZzdTfj2O7nfdinfU3ZPnSJURZGWVlZSxfuqRR8vnhjDd4+tarOe200+pPLElFrJDtbjRBuzvi8OMb5fklyREstUj3nHccX372KWXlbRn1wytZZ/0NAXj8up/y+3PfYsHSSjp324KDfnRV1THLlnzJ1Efv4/ibHgBg56NO5p4fHMcz66/DyItuLljeVk5hWb50CZ06d+FbP7icPffckxfzHDGSpGJyz3nH8fBXn/NpRdlq7e6HM14niNXa3S+/rL3dbVPeltFX/Lpgeaut3XUFQUmNxQBLLdJ/3fZorduP+NnNXDC4S63T39qt05ETf/1/VX9vNWRHzrr/2TrTn/Sb8av83bnb1znrgecA2Hr7ndh6+9VXwtp6+5245Nl38yqLJJWC/7rt0VrbyyN+VvcHVB071t7u1uWZZ55Z5flzaXdHjhxpuyupSTlFUJIkSZIKxABLkiRJkgrEAEs5q6h0WVtJKha2yaUr39fO11oqLd6DpZy5dLckFY9822Tb4+Lhaye1bI5gSZIkSVKBGGBJkiRJUoEYYEmSJElSgRhgSZIkSVKBGGCVgDlz5rD77rvTu3dv+vbty3XXXbdamk8//ZSDDz6YAQMGMHz4cF5//fWqfddddx39+vWjb9++XHvttasde9VVVxERLFiQ3wIWUnUVlYmlS5cyfPhwBg4cSN++fbn44otXS3f22WczaNAgBg4axLbbbsuGG24IwPvvv8/QoUMZNGgQffv25ZZbblnt2FGjRtGvX79GL4uUiyeeeIJevXrRs2dPxowZU2e6Bx98kIjgpZdeqtr26quvsuOOO9K3b1/69+/P0qVLAVi2bBknnXQS2267Ldtttx0PPfRQo5dDxe/444+na9eudbZ/KSXOOOMMevbsyYABA5g6dSpQd7u6aNEiBg0aVPXTpUsXzjrrrCYrj9TSuYpgCSgvL+d///d/GTJkCIsWLWLo0KHsvffe9OnTpyrN5ZdfzqBBg3jkkUf4+9//zmmnncaf/vQnXn/9dX7zm98wZcoU2rVrx7777sv+++/PNttsA2SCt4kTJ/L1r3+9uYqnFqK8LPjFm4vY/5r7ad+xEyuWL+eWEw7go6135OsDtq9Kt+kxP2L0MT/igsFduOGGG5g2bRoAm222GX/9619p3749ixcvpl+/fowaNYpu3boB8PDDD9OpU6dmKZtU04oVKzjttNOYOHEi3bt3Z9iwYYwaNWqVdhkyF7LXX389O+ywQ9W2iooKjj76aO6++24GDhzIxx9/TNu2bQG47LLL6Nq1KzNmzKCyspJPPvmkScul4nTsscdy+umnc8wxx9S6//HHH2fmzJnMnDmTyZMnc8oppzB58uQ1tqvTp0+vOn7o0KF8+9vfbqriSC2eI1glYLPNNmPIkCEArLfeevTu3ZsPPvhglTRvvvkme+65JwDbbbcds2bN4l//+hdvvfUWI0aMoGPHjpSXl7PbbrvxyCOPVB139tln8/Of/5yIaLoCqcWKCNp3zARBKyqWU1mxHNZwbo0dO5YjjzwSgHbt2tG+fXsAvvrqKyorK6vSLV68mKuvvpof//jHjZh7KXdTpkyhZ8+ebL311rRr147Ro0czfvz41dJddNFFnH/++XTo0KFq25NPPsmAAQMYOHAgABtvvDFt2rQB4LbbbuPCCy8EoKysjC5dXJ5bsOuuu7LRRhvVuX/8+PEcc8wxRAQjRoxg4cKFfPjhh2tsV1eaOXMm8+fPZ5dddmm0/EutjQFWiZk1axbTpk1b5dNQgIEDB/Lwww8DmY7//fffZ+7cufTr149nn32Wjz/+mC+//JLHHnuMOXPmADBhwgQ233zzqk5eKoTKFSu4fvRILturNz13GMnX+w+tNd3777/Pe++9xx577FG1bc6cOQwYMIAtttiCH/7wh1WjVxdddBHnnnsuHTt2bJIySPX54IMP2GKLLar+7t69+2offE2bNo05c+ZwwAEHrLJ9xowZRAT77LMPQ4YM4ec//zkACxcuBDLn+5AhQzjssMP417/+1cglUSmo74uGa56Pm1c7H+tqV1caO3YsRxxxhB+0SgXkFMESsnjxYg455BCuvfZa1l9//VX2XXDBBZx55pkMGjSI/v37M3jwYMrLy+nduzc//OEP2XvvvenUqRMDBw6kvLycL7/8kssuu4wnn3yymUqjlqqsTRvOGPcMSxZ9xu/O/R7/fOctvtaz92rpxo0bx6GHHlr1yT3AFltswauvvsq8efM46KCDOPTQQ/nwww955513uOaaa5g1a1YTlkSqW0qrX/BWv0CtrKzk7LPP5o477lgtXUVFBZMmTeLFF1+kY8eO7LnnngwdOpSBAwcyd+5cdtppJ66++mquvvpqzjvvPO6+++7GLIpKQHlZ8Ms3PmHB0hW1fkHxO599xe9mLGTSupl9wb/Px9ra1U033bTq2HHjxnmOSQXmCFaJWL58OYcccghHHXVUrfOk119/fW6//XamT5/OXXfdxUcffcRWW20FwAknnMDUqVN59tln2Wijjdhmm234xz/+wXvvvcfAgQPp0aMHc+fOZciQIfzzn/9s6qKphVpnvQ3YauhOzPjr07XuHzduXNX0wJq6detG3759ee655/jb3/7Gyy+/TI8ePdh5552ZMWMGI0eObMScS/Xr3r171WwAgLlz564yMrBo0SJef/11Ro4cSY8ePXjhhRcYNWoUL730Et27d2e33XajS5cudOzYkf3224+pU6ey8cYb07FjRw4++GAADjvssKrFCqQ12aBrNxb+a17V3zXPR1i1XV3plVdeoaKigqFDa59pIKlhDLBKQEqJE044gd69e3POOefUmmbhwoUsW7YMgFtvvZVdd921apRr/vz5AMyePZuHH36YI488kv79+zN//nxOfuQlTn7kJdbr2o3v3TGROz4sZ8y0BbX+SPVZ/OkCPLYgEAAAIABJREFUliz6DIDlS5fwj8l/YZMe26yW7qNZ7/Dpp5+y4447Vm2bO3cuS5YsATKrYj7//PP06tWLU045hXnz5jFr1iwmTZrEtttuyzPPPNMk5ZHqMmzYMGbOnMl7773HsmXLGDduHKNGjarav8EGG7BgwQJmzZrFrFmzGDFiBBMmTGD77bdnn3324dVXX+XLL7+koqKCv/zlL/Tp04eI4Fvf+lbV+f2nP/1ptUUzpNr03m0fpj16HyklZr/6EhtssAGbbbZZne3qStXvg5VUOE4RLAHPP/88d999N/3792fQoEFAZtXA2bNnA3DyySfz1ltvccwxx9CmTRv69OnDb3/726rjDznkkKpVqm666SY6d+7cLOVQy7foo3/xwMWnk1ZUklIl/fc+kN67/icTfzmGzfsMos9u+wLwyhMPM3r06FWmVL311luce+65RAQpJc477zz69+/fXEWR1qi8vJwbb7yRffbZhxUrVnD88cfTt29ffvKTn7D99tuvEmzV1LlzZ8455xyGDRtGRLDffvux//77A3DllVfy3e9+l7POOotNNtmE22+/vamKpCJ25JFH8oennuaLhZ9wxb4D2Ovk86msqABgh0OPpdfOe/P2pKe46sDhtO2wDo+OvQuov129//77eeyxx5qlTFJLZoBVAnbeeeda5/tXt+OOOzJz5sxa91WfDlCXH/7BaShae5tt25czxv55te17n3LBKn/vdfL5XDB41dXR9t57b1599dU1Pn+PHj1W+Y43qTntt99+7Lfffqts++lPf1pr2pqjrkcffTRHH330aum23HJLnn322YLlUS3D2LFj1ziTJCI48MKfV/09aODGQP3t6rvvvlv1uKIyUV7mQhdSIRhgSWoW+Xbmdv6SlJvyssh7an/ND70kNZwBVgviBahKSb4XAHb+kiSpFBhgtSBesEpS6fJDMjUnZxVIhRP13dtTTc4J1XzyDbCqp79wyCaNkSW1QldM/ajeNDXPv1zSq8Vpiquzkuq7ar4nbJe1Ui7t6kr5tq8NOcY2Wa1YvX2Xy7RLKgkVlflfJzfkGEmSpLXhFMEi5dC7tCpv2pak4uGUQqluBlhFyotJae15ASBJjSPf65TzskvH58r2WKXMAKuB8n3jL69MtC3yhqK++d1NNae7MdObp8ZJ39BjGpsXAGouhTo31tQuF+P7utjSt9Y8FaPGbo8h/2st23A1lpwXubj00kufAJpyiKQbMK8J/1+hmf/mVer5h9Ivg/lvXqWQ/wUXX3zxvo35D2r0XaVQJw1l2UpXSy6fZStdLbl8a1u2+vuulFJR/lxyySWpufNg/ps/H601/y2hDObf/BfbT0uuE8tWuj8tuXyWrXR/WnL5mqJsriIoSZIkSQVSzAHWpc2dgbVk/ptXqecfSr8M5r95lXr+G0NLrhPLVrpacvksW+lqyeVr9LLl80XDkiRJkqQ1KOYRLEmSJEkqKQZYkiRJklQgTR5gRcQWEfHniHgrIt6IiDNrSTMyIj6LiOnZn59U27dvRLwdEe9ExAVNm/uc8/+Danl/PSJWRMRG2X2zIuK17L6XmiH/HSJiSkS8ks3/avNQI6J9RNyXrePJEdGj2r4Ls9vfjoh9mjLv2f+fS/7PiYg3I+LViPhTRGxZbd+Kaq/NhKbNfc75PzYiPqqWz+9X2/e9iJiZ/fle0+Y+5/xfUy3vMyJiYbV9zVr/1fLRJiKmRcSjtewr2vO/Wj7WlP+iPf8bS339QinXSQ5lK9r2Ihc5lK/o25O6RMRtETE/Il6vY39ExPXZsr8aEUOq7Svq1y6Hsh2VLdOrEfHXiBhYbV+zXgfVJ4eyFe01ai5yKF/RXsPWJ3K7Rm+a911TL40IbAYMyT5eD5gB9KmRZiTwaC3HtgH+AWwNtANeqXlsMeS/RvpvAU9X+3sW0KW5lqYEAuiUfdwWmAyMqJHmVOCW7OPRwH3Zx32ydd4e2Cr7WrQpwvzvDnTMPj5lZf6zfy9urrrPI//HAjfWcuxGwLvZ352zjzsXW/5rpP9/wG3FUv/V8nEOcG8d7UzRnv855r9oz/9Gqot6+4VSrZMcy1a07UUhylcjfVG2J2vI767AEOD1OvbvBzyebVdHAJNL6LWrr2zfWJln4Jsry5b9exbNeB1UgLKNrKPtbfZr1EKUr0baorqGzSG/ucQYTfK+a/IRrJTShymlqdnHi4C3gM1zPHw48E5K6d2U0jJgHHBg4+S0dg3I/5HA2KbIWy5SxuLsn22zPzVXOjkQuDP7+EFgz4iI7PZxKaWvUkrvAe+QeU2aTC75Tyn9OaX0ZfbPF4DuTZjFNcqx/uuyDzAxpfRJSulTYCLQqF/SWlMD8l9U5z9ARHQH9gdurSNJ0Z7/UH/+i/n8byT19gslXCdr0+c1e3uRg3zLV3TtyZqklJ4FPllDkgOBu7Lt6gvAhhGxGSXw2tVXtpTSX7N5h9J6z+XyutWl2a9Rc5Fn+UrtPZfLNXqTvO+a9R6syEy9GUzmU/CadozMNKTHI6JvdtvmwJxqaeaSe3BWcPXkn4joSObFeaja5gQ8GREvR8RJjZ3HOvLVJiKmA/PJnEw1819VzymlCuAzYGOKpP5zyH91J5D5pGKlDhHxUkS8EBEHNWpG65Bj/g/JDl0/GBFbZLeVVP1HZhrWVsDT1TY3e/0D1wLnA5V17C/q85/6819d0Z3/jSDf16WU6iTXshVte1GPnPNYxO3J2qir/KXw2uWj5nuu2a+DCqDor1HXVrFew+ZqDdfoTfK+K2/ogWsrIjqRedHOSil9XmP3VGDLlNLiiNgP+D9gGzLDeTU1yzrz9eR/pW8Bz6eUqn9SsFNKaV5EdAUmRsTfs58mNJmU0gpgUERsCDwSEf1SStXn4tZVz0VR/znkH4CIOBrYHtit2uavZ+t/a+DpiHgtpfSPpsl5Rg75/z0wNqX0VUScTGY0ZQ9KrP7JTK97MJt+pWat/4g4AJifUno5IkbWlayWbUVx/ueY/5Vpi/L8bwQ5vy4lWCe5lK2o24t65JPHomtPCqBo25pCiYjdyQRYO1fb3OzXQWup6K9RC6Qor2FzUc81epO875plBCsi2pIp+D0ppYdr7k8pfb5yGlJK6TGgbUR0IRNNblEtaXdgXhNkeRX15b+a0dQYWk0pzcv+ng88QjNMMaqWl4XAM6w+BFpVzxFRDmxAZji5KOp/pTXkn4jYC/gRMCql9FW1Y1bW/7vZYwc3RV5rU1f+U0ofV8vzb4Ch2cclU/9Zazr/m6v+dwJGRcQsMtM39oiI39VIU8znfy75L4nzv4Byel1KtE7qLVuptBd1yCePxdierK26yl8Kr129ImIAmanMB6aUPl65vZiugxqi2K9RC6ior2HrksM1etO871LT34AWwF3AtWtI8zX+/SXIw4HZ2ePKydx0thX/voGwb7HlP5tu5UXZutW2rQusV+3xX4F9mzj/mwAbZh+vAzwHHFAjzWmsepP//dnHfVn1Jv93afpFLnLJ/2AyN5puU2N7Z6B99nEXYCZNv0hKLvnfrNrjg4EXso83At7LlqNz9vFGxZb/7L5eZG6GjWKq/xp5HEntNyoX7fmfY/6L9vxvpHqot18o1TrJsWxF214UonzZdEXfnqyhjD2oe7GE/Vn1ZvsppfLa5VC2r5O5T/UbNbY3+3VQAcpWtNeohShfdn9RXsPmUK5cYowmed81xxTBnYDvAq9l7+MA+G8yb0ZSSrcAhwKnREQFsAQYnTKlr4iI04E/klmt5baU0htFmH/IdHRPppS+qHbspmSmVEHmjXhvSumJJsn1v20G3BkRbciMYN6fUno0In4KvJRSmgD8Frg7It4h8wYbDZBSeiMi7gfeBCqA09Kq0zWKJf+/ADoBD2TrenZKaRTQG/hVRFRmjx2TUnqzCPN/RkSMIlPHn5BZJYyU0icR8T/Ai9nn+mladei+WPIPmRtjx2XftysVQ/3XqoTO/1qV0PlfcCmlWvuFllAnOZatmNuLNcqxfFBi7clKETGWzAchXSJiLnAxmYWBVl4rPEZmRbN3gC+B47L7iv61y6FsPyFz7+rN2fdcRUppe4rjOmiNcihbMV+j1iuH8kHxXsPWJ5dr9CZ538Wq7ZUkSZIkqaGadRVBSZIkSWpJDLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsKRqImJxIz//OhHxl4hos4Y07SLi2Ygob8y8SJJKW0QcFhFvRcSfs3+PjYhXI+LsPJ9nw4g4Nc9jbo2IPvkcI7UWkVJq7jxIRSMiFqeUOjXi858GlKeUrqsn3cXAOymlexorL5Kk0hYRTwBXppT+HBFfAyanlLZswPP0AB5NKfUrcBalVskRLKkWEXFORLye/Tmr2vaLIuLvETEx+0nheXk+9VHA+OxzbRgR/6z23C9HxAbZP/8vm1aS1MpFxNERMSUipkfEryKiTUT8BNgZuCUifgE8CXTNptklIv4jIp7I9i3PRcR22efaNCIeiYhXsj/fAMYA/5E99hc1/ve6EfGHbNrXI+KI7PZnImL7iBiVPW56RLwdEe9l9w/Nzth4OSL+GBGbNWWdSc3JKUhSDRExFDgO2AEIYHJE/AVoAxwCDCbz3pkKvJzH87YDtk4pzQJIKS3MdlxtU0rLgVeAAcBzwOvAsIIVSpJUkiKiN3AEsFNKaXlE3AwclVL6aUTsAZyXUnopIm4iMwo1KHvcn4CTU0ozI2IH4GZgD+B64C8ppYOz09U7ARcA/VYeW8O+wLyU0v7Z592g+s6U0gRgQnbf/cBfIqItcANwYErpo2xQdhlwfEErRypSBljS6nYGHkkpfQEQEQ8Du5AZ8R2fUlqS3f77lQdExNbAj4ANUkqHRsS6ZDqzZcAz2al+XYCFNf7Xv4CvAXOA7bJ/k1JaERHLImK9lNKixiuqJKnI7QkMBV6MCIB1gPlrOiAiOgHfAB7IHgPQPvt7D+AYyPQ1wGcR0XkNT/cacFVEXEkmgHuujv95PrAkpXRTRPQD+gETs/+/DfBhPeWUWgwDLGl1ked2UkrvAidExIPZTd8GHkwp/T4i7gPuAZYAHWocOg/olp2isSClNKPavvbA0oYUQJLUYgRwZ0rpwjyOKQMW1jEilZeU0ozszI79gCsi4smU0k9XyWDEnsBhwK7V8vxGSmnHtf3/UinyHixpdc8CB0VEx+xI1MFkpu1NAr4VER2ynw7uv4bn6E5mVApgBUBK6VOgTURUD7LmAQeRmZ5RNXUiIjYGPspOHZQktV5/Ag6NiK4AEbFRRKxxIYuU0ufAexFxWPaYiIiB1Z7vlOz2NhGxPrAIWK+254qIbsCXKaXfAVcBQ2rs35LMjI3DV87wAN4GNomIHbNp2kZE3zzLLZUsAyyphpTSVOAOYAowGbg1pTQtpfQimXnmrwAPAy8Bn9XxNHPJBFmw6vvsSTJTEFf6ADgUGJVSWlBt++7AY2tXEklSqUspvQn8GHgyIl4FJgK5LBhxFJmZFa8AbwAHZrefCeweEa+RuY+4b0rpY+D57CIWv6jxPP2BKRExncxU+J/V2H8ssDHwSHahi8dSSsvI9G1XZv//dDJTFqVWwWXapTxERKeU0uKI6EhmpOuklNLU7IjTZcDewK1kbiK+kcwUv0krl1uPiMHAOSml79bzfx4GLkwpvd2IxZEkSVKBGWBJeYiIe4E+ZO6lujOldEUDnuP47LEr6tjfDhidUrprrTIrSZKkJmeAJUmSJEkF4j1YkiRJklQgBliSJEmSVCD5fA+WcwklSYVU53fLFZB9lySpkOrtuxzBkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAkiRJkqQCMcCSJEmSpAIxwJIkSZKkAjHAklqBisr8VqrON70kqXg0pA233ZcKJ1LK+Q3lO08qYWOmLcg57QWDuzRiTqQqfg+WlIOKykR5WX5vl3zafLDdl/JQ75sxny8aliRJUhMrLws/JJNKiFMEJa3G6SWSJEkN4wiWVIIaMl0kH/l+Wgp+YipJkgQGWFJJcrqIJElScXKKoCRJkiQViAGWJEmSJBWIAZYkSZIkFYgBlqSC8MuMJal02YZLheMiF5IKwoU3JKl02YZLheMIliRJkiQViAGWVAScaiFJktQyOEVQKgKFmpox/72ZPHjJGcz7+6v852n/za7HnLbG55lw5QW8PGEslz7/PgDP/e6XvPTI77i3U3u+7LAhh1x8HZ27bQHAwg/n8tD/nM1n//yAiODYG8bSudvXc86zJCmjsb8sPqXEhJ9fyNuTnqJdh44ceun1bN574Grpfn3igSxa8C/atu/AuHXKOeB/x9Jpo02q+oKyNuWs23njWvuCOxf+k4jgscceo0ePHo1WFqkUGWBJLUjHDTbkW+dfzpt/fqzetHPfnM6SRZ+vsq1br/6c9ruJ/OQbX+eg//4Fj193Kd+58lYA7v/Jaex+wtlsM2IkX325mAgHwCWpIRr7fqfHH3+cj2e/y3njpzDntZf5vyvO57S7/lhr2iMuu4XufQZxweAuVXla2Re0W6cjLzxwe619wW9POZTFixdTVmZfINXku0JaC7NmzWK77bbj+9//Pv369eOoo47iqaeeYqeddmKbbbZhypQpfPHFFxx//PEMGzaMwYMHM378+Kpjd9llF4YMGcKQIUN4/5UpALz70vP8+sQDuecHx3H1t3dk3I9OJqXcphB22mgTtug7mDblbdeYrnLFCh6/9hK+eeZPVtn+H8N2pt06HQH4ev+hfD5/HgD/evdtKldUsM2IkQC079ipKp0ktXYN6QvefOZxAD6dN5tfHX8AN3xnD274zh619gXbbbddXn3B+PHjGXzAEUQEXx+wPUsXfcbnH/0z5/Lk2hd06tSJjh3tC6SaHMGS1tI777zDAw88wK9//WuGDRvGvffey6RJk5gwYQKXX345ffr0YY899uC2225j4cKFDB8+nL322ouuXbsyceJEOnTowMyZM9n9wMM4/Z6nAJj39muc/cAk1tvka9xy3P68P30yPQaP4NGrfsy7L01i3DrlzF9SUZWHAfsczMjjzsw5z3+771Z677ov62/ytTrTvPh/97DtTnsCsOD9f9Ch0wb87txj+WTe+/Qcvhv7nnERZW3aNLDWJKllybcv6DlwKD132JV1O3fh+F8+SNv2HVgw+x+Mu/C/VusLLvvPfvQYtMNqfUFNK/uCDz74gO47fKtq+wZdu/H5R/+stc1/8JIzKCsrY/nRR1C+/8lErDp1sa6+4L5PP2CvvfZizJgxtLEvkFZhgCWtpa222or+/fsD0LdvX/bcc08igv79+zNr1izmzp3LhAkTuOqqqwBYunQps2fPplu3bpx++ulMnz6dNm3aMP+9GVXPuUXfIWywaTcAuvXqx6fz5tBj8AgOOO9nAKtM5cjX5x/9k9eemsCJvx5fZ5rf/e53fPDmK5x0ayZN5YoKZk1/gTPufZoNvtadsRd8n5d/P5ZhBx3doDxIUkuTb19QsewrFn74Aetv8jXGX3kBH854nbKyMhbMfrfqOVf2BWVlZbX2BXWpbaQrWP2eryMuu4UNum7GV18s5rmfnsSGbMSQA46o2j/tDw/U2Rdc/s1BHHHEEdxxxx2ccMIJ+VeY1IIZYElrqX379lWPy8rKqv4uKyujoqKCNm3a8NBDD9GrV69VjrvkkkvYdNNNeeWVV6isrKR9hw5V+9q0bVf1OMrKqFyRGa0qxAjWvL+/ysdz3uOqA4cDsHzpEn4xahg/mPAiAO9M/gt3XncZx9z4MOXtMmXZoGs3uvXqz0bdewDQZ+R+zHntJTgop38pSS1evn3Byg/Jnrrl56y38SYcNu4ZUmUlP9mxe9Xz1NcX1LSyL+jevTsL/zWvavtn8+ex3iabrpZ+g66bZfK+bie+853v8KvHn6sKsN6Z/Bf+/NtrOOnW8bX2BeXl5Rx00EG88MILBlhSDQZYUiPbZ599uOGGG7jhhhuICKZNm8bgwYP57LPP6N69O2VlZdx5551UrlhR73MVYgRru13+kx9NfLPq74t32rIquJr391d55LLzeOHpJ3loceeqNN37DmbJ55+x+NMFdOrchXdffI7N+wxq0P+XpNaoZl8w7++v0m27ASxd/HnVKNVLvx+XV19Ql1GjRnHO5VczcJ+DmfPay3TotP5q0wNXVFSwdNFnrNt5Y1YsX86jjz7K13qNAP7dFxx34zg6bbRJ1THV+wLowtNPP83222+ff2VILZyLXEiN7KKLLmL58uUMGDCAfv36cdFFFwFw6qmncueddzJixAhmzJhRkEUjFi34F1fsO4BJ9/ySP996NVfsO4ClixcBcPv/G13vTc6PXXspy778gsMOO4zrR4/krrMyUwDL2rRhv7Mv4bf/dQjXHr4ricSwb393rfKa73d/+V1hkkpZzb7gyZvHADDi8OOY+vv7uPmYfVkw+x8F6Qv2228/Ntp8S646cDgP/+wcDrzw51X7rh89EoAVy7/ittMO57rDd+P6I3dn8803Z9jBmXZ9ZV9w7/kn1NkX9O/fn5QSJ5544lrnV2ppItcVaQCvbqRGlO+SvcWUvqnypBan8b4I6N/su1SUWkIbLrVS9fZdjmBJkiRJUoEYYEmSJElSgRhgSQXmvUKSJEmtl6sISgVWXhYNmvsuSZKk0ucIliRJkiQViAGWJEmSJBWIAZYkSZIkFYgBliRJkiQViAGWJEmS8pLvirmusKvWxFUEJUmSlJd8V8x1tVy1Jo5gSZIkSVKBGGBJkiRJUoEYYEmSJDWQ9xZJqsl7sCRJkhoo33uRwPuRpJbOESxJkiRJKhADLKkeTv+QJElSrpwiKNXDpWglSZKUK0ewJEmSJKlADLAklYSGTNV0eqckSWpqThGUVBJcqUuSJJUCR7AkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLAkSZIkqUAMsCRJkiSpQAywJEmSJKlADLDU6lRUpubOgiRJrUq+fa99tUpZeXNnQGpq5WXBmGkLck5/weAujZgbSZJaPvtetSaOYEmSJElSgRhgSZIkSVKBGGBJkiRJUoEYYEmSJGW5uIKkteUiF5IkSVkuxiBpbTmCJUmSJEkFYoAlSZIkSQVigCVJkiRJBWKAJUmSJEkFYoAlSZIkSQVigCVJkiRJBWKAJanFyvf7bPz+G0mStLb8HiyVtIrKRHlZNHc2VKT8PhtJktTUDLBU0vK9gAYvoiVJktR4nCIoSZIkSQVigCVJkiRJBWKAJUmSJEkFYoAlSZIkSQVigCVJkqSi0pCvzfCrNlQsXEVQkiRJRcVVglXKHMGSJEmSpAIxwJIkSZKkAjHAkiRJLZb35Uhqat6DJUmSWqx87+XxPh5Ja8sRLEmSJEkqEAMsSZIkSSoQAyxJkiRJKhADLBUVb0aWJElSKXORCxUVb0aWJElSKXMES5IkSZIKxABLkiRJkgrEAEuSJEmSCsQAS5Ky8l1kxUVZJElSTS5yIUlZLrIiSZLWliNYkiRJklQgBliSJEmSVCAGWJIkSZJUIAZYkiRJklQgBliSpP/f3p2HR1Xdfxx/nyzsIMhmWNRgWAIYwiZhqSAhBkMNoKgUUEBRgyhFtIKlKvDQgrXFpeLC5i8gshUxFAgqAqKhYhTCXghCi0AUwipLIIn398ck82QlkzCZmcx8Xs8zD5O5Z2bOOZy533PuPfdcEZEKTyvBiqfQKoJSrrJ+tQjwM+7OhoiIiHg5rQQrnkIDLClXvrazS/rofdo+vJiTl7PoPHAYPYbGAbD+vb+SvHIh1evUBeDupyfRqkcU/03ZSsJfXmBF7Wr0/NMs6t3cjMu/nGPxhFGMnLUMYwoPTmc/3p8+770J/rcCcOb4EeJ/P5Rxy7/i0HdJLBj/MDc2voWsK1cIix5Inyf/wKZNm5h8byw3Nr6FzIzL1LixPncOf4bQO+92Wd2IiPiypI/eJ3nlh1iWRZVn4qDnMKDk+OBfqRKDEpYDtR2KDzHPToH2fQDH4kPu60tDbuPomV8UH0ScQAMsESf56eA+kld+yKGd3/P6nvN88PRDtPpNFPVuvg2A7kPjuPORMfne8/XCdxn6t/lEVzvP1AX/R7/xU9kw5+/0emxckcHTEbeGRzDirY+4evkibw2+i9Df3A2N/e2vAxzfv4uF44cTWLkKIV3uvL6Ci4jINeXGh6cWfIp/YCVWvziM8OCuDsWHM8d/5N1336X+sBfLJz7kvL79q8+YsT1d8UHECXQNloiTnDx8gKa3d6RatWr4BwQQ3LEbezasveZ7/AICyMzI4NKlS/gHBHDqx8OcP5FGs47drzs/lapWp3FoO04dPVxoW6OWtxP5+HP8e9m86/4eERFXqajXzOTGh0pVbfGhZ8+eDseHzIzLBAYGKj6IVCA6gyXiJA1vC+XTWX/h1KlTXL18if1fr6dJ63b27f9eOo/tq5fRuHU7+o2fStVaten16O9ZOe05dtarSdcJb7L29VeIempiid81dOhQfiEQgOzMTIxf4WMlF8+e5siu7+n9+HNAVqHtjULD2LxgVtkLLCLiYqWddg6eMfU8Nz5cPHuawMpVSFy7FnNLG/v2a8WHwCpV+HLlEvo99oxD8WHppDi+fq0GJy5nORQfLp45VWi74oPI9dEAS8RJGjRrQc8RzxAVFcUZKhPUog1+/rafWJcHRtgGOsbw+TvTWTPzZQZNfotGLW/nqQXrmNi+Hk/OXUXN+jdhWfDRhFH4BwQSM34KNes2KPRdixYtYn2Ba7By/TflG9763V0Y40evkWNpeFsrOLe7cIatinkkWESkosmND/OfGkSlqtWJvqMdu3+xbSspPgAcOrRbHgy5AAASI0lEQVTX4fjw0J/f4+2hfZixPd2h+HDou6TCGVZ8ELkuGmCJOFHnAcNYMWUcM7an8+k/plGrYSOAfEHwjvsezhfwACzLYsO81xkyYw4Jr06kT9wEzhw/wpbFc4h+elKp8pD3WqtrOf6fXTQIblGqzxYRkbLpPGAYnQfYFrY4v3wmdevUBhyLD9OmTSNy0juKDyIVhK7BEnGiC6dPAnA27Sh7Nq4hvO99AJw/+ZM9zZ4Na21nlfKIj4+nVY8+VK1Vm8yMyxg/P4yfH5kZl8sln2kH9rBh7kwiHny0XD5fRETyyxsfPv74Y4fjw7Z/LaFfv36KDyIViM5gicN0T6uSLXp+JB9fOc+ZLD9iJ7xK1Vq2I5SJb04l7cBuDIY6jZoyYNLf7O+5evkS8fHx9HnVdlSxx9A4Fv1hJP4BgQyePttpecudGpKZcZkadepx7x/+ohWiRERcZNHzI7l07gx+AYF8+O4skh2MD9tWL+WjLRv5++5z5Rof2rdvb1umXfFB5LppgCUOq6gXF7vSk/NXM7F9vUL19NC0d4p9T6Wq1di4caP9PcEdujJu2eZi0z8xJ4FO7euxPid9nUY3M275VwA069SdZp0KrzDVq1cvJm8+VOryiIiIczw5f7X9eWT7eiTn7MNLig+Pz/6EwEDbokaOxIe8HIkPzTp1Z/LmQ0XGLhEpG00RFBEpo7IsGV1Rl5kWERERx+gMlohIGemsroiIiBSkM1giIiIiIiJOogGWiIiIiPic0k7Z1hRvcZSmCIqIiIiIzyntNG9N8RZH6QyWiIiIuIXOCIiIN9IZLB+m+1qJiIg76QyCiHgjDbB8mAKbiIiIiIhzaYqgiIiIiIiIk2iA5eHWrVtHy5YtCQkJYcaMGUWmWbZsGa1bt6ZNmzYMGTLE/voLL7xAmzZtCA0NZezYsViWba77pEmTaNq0KTVq1HBJGUR8wdmfjjHniQHMvK8brw/qQdJH7xdKY1kWY8eOJSQkhLCwMLZt22bfNmHCBNq2bUvbtm1ZunSp/fUvvviCDh06EB4eTo8ePTh48KBLyiO+5VptM6++ffvSrl072rRpQ1xcHNnZ2QCcPn2aqKgomjdvTlRUFGfOnAFg0aJFhIWFERYWRrdu3dixY4fLyiRSVv+cPJZpkaG88cBvityekJBAWFgY4eHhdOrUia+//tq+zd/fn/DwcMLDw4mNjXVVlsXDaIDlwbKzsxkzZgyJiYns3buXxYsXs3fv3nxpUlNTmT59OklJSezZs4c33ngDgC1btpCUlMTOnTvZvXs3ycnJfPnllwDce++9fPvtty4vj4g38/P3J+bZKYz/eAtPxa/j38vm8/Oh/fnS7E9aT2pqKqmpqcyePZvRo0cDsGbNGrZt20ZKSgpbt27ltdde4/z58wCMHj2aRYsWkZKSwpAhQ5g2bZrLyybeLzExsci2WdCyZcvYsWMHu3fv5uTJkyxfvhyAGTNmEBkZSWpqKpGRkfYDgsHBwXz55Zfs3LmTl156iSeeeMJlZRIpq473Dmbk20uK3R4ZGcmOHTtISUlh/vz5jBo1yr6tatWqpKSkkJKSwqpVq1yRXfFAGmB5sG+//ZaQkBCaNWtGpUqVGDx4MAkJCfnSzJkzhzFjxlCnTh0AGjRoAIAxhoyMDK5evcqVK1fIzMykYcOGAERERBAUFOTawoh4uVr1b6JxaDsAKlevQYPgFpw/kZYvzb5N63jkkUcwxhAREcHZs2dJS0tj79699OzZk4CAAKpXr067du1Yt24dYPst5w62zp07R6NGjVxbMPEJCQkJRbbNgmrVqgVAVlYWV69exRhjf//w4cMBGD58OJ988gkA3bp1s8eniIgIjh496oriiFyX4I7dqHZDnWK316hRw972L168aH8ukksDLA927NgxmjZtav+7SZMmHDt2LF+aAwcOcODAAbp3706XiAh7p6xr167cddddBAUFERQURHR0NKGhoS7Nv4ivOnP8CMf376Jp2475Xj93Iq3I33S7du1ITEzk0qVLpKens3HjRn788UcA5s6dS0xMDE2aNGHhwoVMnDjRpWUR3+BIvMkVHR1NgwYNqFmzJoMGDQLg559/th+4CwoK4sSJE4XeN2/ePO65555yyL2I661cuZJWrVrRr18/5s+fb389IyODTp06ERERYT/QIL5Hqwh6sNxrpvIqeJQkKyuL1NRUNm3axNGjRwmP6M645V9x8cwpVm/dwbNrUgCYN/oB0oO7ENyxm0vyLuKrrly6wIfPj+S3z02jSo2aBbZahe77Y4zh7rvvJjk5mW7dulG/fn26du1KQIBt9zxz5uusXbuWLl268NprrzF+/Hjmzp3rotKIr3Ak3uT69NNPycjIYOjQoWzYsIGoqKgSP3/jxo3Mmzcv37UqIhVN3tvbDBw4kIEDB7J582Zeeukl1q9fD8CRI0do1KgRhw4donfv3oS2aUvL5iHuzLa4gQZYHqxJkyb2o9gAR48eLTQ9qEmTJkRERBAYGEhwcDD1bwkh/cghDn+XRNPbO1G5mm0hi5bdIzmy63sNsETKUXZmJoueH0l4zCDaRv620PYbGjQi7dhR++0Rdv7wP1acqszn29Pxj3mSwTFPArDkj08S6NeAP234Dzt37qBLly4APPTQQ/Tt29d1BRKvNmvWLObMmQNA586dS4w3eVWpUoXY2FgSEhKIioqiYcOGpKWlERQURFpamn26OsDOnTsZNWoUiYmJ1K1bt/wKJFLOiry9Tc3WbNuXyksb9lO9Tl2gEvycDtSiblgEu3akaIDlgzRF0IN17tyZ1NRUDh8+zNWrV1myZEmhFWkGDBjAxo0bAUhPTyf9yA/c2PgWat/UmMPfbyE7K4vszEwOf7+FBsEt3FEMEZ9gWRYrpo6jfnALfjOs6AUCQntGs2DBAizL4sjO76hSoxa16t/Er9nZXDx7GoC0A3v4KXUvzSPuomrN2pw7d44DBw4A8Pnnn2uqrzjNmDFj7BfjDxgwwN42v/nmG2644YZC1+peuHDBfl1WVlYWa9eupVWrVgDExsYSHx9P1q8W8fHx9O/fH7Adzb/vvvtYuHAhLVooBol3SD9yyH7W99i+HWRnXqVa7Ru5fP4sWVevAHDxzCn+l/ItrVu3dmdWxU10BsuDBQQE8PbbbxMdHU12djaPPvoobdq04eWXX6ZTp07ExsYSHR3NZ599RuvWrfH39+eecZOpXvtG2vaJ5Yfkr3nzwTsxxtC8W29Ce0YDkPjGFFLWreDSpUtM7xtG5wHD6BP3gptLK1Kx/S9lK9vXLOOmkNa8NbgXAHc/PYlzP9muY+kyaAQte0RRfX8Sf+t/B4FVqjJo8lsAZGdlMvuxewGoXL0mD057B/+cKYJz5szh/vvvx8/Pjzp16uSb6y/iLDExMaxdu5aQkBCqVavGBx98YN8WHh5OSkoKFy9eJDY2litXrpCdnU3v3r2Ji4sDYOLEiTz44IPMmzcPq04QQ/46jxnb01kx9Y8cO5HOoJG21QP9/AN4etF6+2frBvbiiRa/+ASHv0/i4tnTTO8bRp+4F/g1K4v3ttaALoPYs2E121Yvwz8ggIDKVfndjDkYYzhx+AAr//w8xvhhWb/Sc+RYDbB8lAZYHip3nm9MTAwxMTH5tk2dOtX+3BjDzJkzmTlzJoD91LWfvz8D//T3Ij/7nnGvcM+4V5jYvl7hU90iUia3to9g+raT10xjjGHWrFk0HZX/dxdYuQrPrkgq8j258/xFylNu2yxKSortWt6GDRuSnJxcZJq6devyxRdfAOSLK/e//Ab3v/yGk3MrUr5+N312ka/H5fSbeo4YS88RYwttv6XdHYxbtjnfa3mv23JEadOLZ9IAy0MVOc+3BDoSKCIiIuI5StufU1/OO+gaLBERD1Zw1UFnpxcRERHn0hksEREPpqOf4k6ariQiUnqmqHtfFEOHRa9DWYJUWaYIlrYjNmN7Oi92qF+q7xGRopV0DRaU/XdamvQViCt67opd18mZ7VXxRiqKovbnZbl23cv34b6qxNilM1guoqPQIiIiInItWhTDO2iAVUZq0CIiUpEobol4Ph2Q9w4aYJWRfgAi4onK0olWx7tiKu3/m1anFfE+2ud7JoevwZoyZco6wNl72kbAcSd/prdRHZVMdeQY1VPJVEclc2Ydpb/yyit9nfRZRSqn2FUa3tKmvKUc4D1lUTk8j7eUxVvKAeVTlpJjl2VZbntMnjzZcuf3V4SH6kh1pHpSHXnSQ3Xkm/XlLeXwprKoHJ738JayeEs53FkW3QdLRERERETESdw9wJri5u+vCFRHJVMdOUb1VDLVUclUR6XjLfXlLeUA7ymLyuF5vKUs3lIOcFNZSnMfLBEREREREbkGd5/BEhERERER8RoaYImIiIiIiDiJBlgiIiIiIiJO4pIBljGmrzFmvzHmoDFmYhHbKxtjluZs32qMudUV+fIkDtTRCGPMSWNMSs5jlDvy6U7GmPnGmBPGmN3FbDfGmLdy6nCnMaaDq/Pobg7UUS9jzLk87ehlV+fRnYwxTY0xG40x+4wxe4wxvy8ijdqRY/Xk022pOMaYB3Lq7FdjTKdrpLvmPt/djDE3GmM+N8ak5vxbp5h02XnawCpX57M43tTv8Jb+gbfEcG+Js94SDz02XpX3OvCAP/AD0AyoBOwAWhdI8xTwXs7zwcBSd6+b78qHg3U0Anjb3Xl1cz3dCXQAdhezPQZIBAwQAWx1d549sI56AavdnU831k8Q0CHneU3gQBG/NbUjx+rJp9vSNeouFGgJbAI6FZOmxH2+ux/AX4GJOc8nAq8Wk+6Cu/NalvqtKP0Ob+ofeEsM95Y46y3x0FPjlSvOYN0BHLQs65BlWVeBJUD/Amn6A/E5z/8JRBpjjAvy5ikcqSOfZ1nWZuD0NZL0BxZYNt8AtY0xQa7JnWdwoI58mmVZaZZlbct5/guwD2hcIJnakWP1JEWwLGufZVn7S0hWEfb5eeNyPDDAjXkpLW/qd1SEtuIQb4nh3hJnvSUeemq8csUAqzHwY56/j1K44PY0lmVlAeeAui7Im6dwpI4A7s85RftPY0xT12StQnG0Hn1dV2PMDmNMojGmjbsz4y45U4LaA1sLbFI7yuMa9QRqS2VVEdpYQ8uy0sDWgQEaFJOuijHmO2PMN8YYTxmEeVO/w5f6BxXhd+GoCrVv9JZ46EnxKqC8vwDbacWCCt58y5E03syR8v8LWGxZ1hVjTBy2I2+9yz1nFYuvtyNHbANusSzrgjEmBvgEaO7mPLmcMaYGsAIYZ1nW+YKbi3iLT7ajEurJZ9uSMWY9cFMRmyZZlpXgyEcU8ZrL29i1ylGKj7nZsqzjxphmwAZjzC7Lsn5wTg7LzJv6Hb7UP6go/yclqVD7Rm+Jh54Wr1xxBusokPdoShPgeHFpjDEBwA14wenXUiixjizLOmVZ1pWcP+cAHV2Ut4rEkbbm0yzLOm9Z1oWc52uBQGNMPTdny6WMMYHYdsKLLMv6uIgkakeUXE++3JYsy+pjWVbbIh6ODK7AQ9pYCeX4OXcqUM6/J4r5jOM5/x7Cdt1Zexdl/1q8qd/hS/0Dj/hdXK+KtG/0lnjoifHKFQOsZKC5MSbYGFMJ28WkBVcaWgUMz3k+CNhgWZZHjpDLSYl1VGDOayy2OaaS3yrgkZxVbyKAc7lTXMTGGHNT7nUGxpg7sO0DTrk3V66TU/Z5wD7LsmYWk8zn25Ej9eTrbek6ORIX3S1vXB4OFBo8GmPqGGMq5zyvB3QH9rosh8Xzpn6HL/UPvGLfW1H2jd4SDz01XpX7FEHLsrKMMU8Dn2JbDWe+ZVl7jDFTge8sy1qFrWIWGmMOYjuCNLi88+VJHKyjscaYWCALWx2NcFuG3cQYsxjbSjD1jDFHgVeAQADLst4D1mJb8eYgcAkY6Z6cuo8DdTQIGG2MyQIuA4M9tFNRXroDDwO7jDEpOa/9EbgZ1I7ycKSefL0tFckYMxD4B1AfWGOMSbEsK9oY0wiYa1lWTHH7fDdmuygzgGXGmMeAI8ADAMa29HycZVmjsK2Y+L4x5ldsHZYZlmW5fYDlTf0Ob+ofeEsM96I46y3x0CPjlfHM/3MREREREZGKxyU3GhYREREREfEFGmCJiIiIiIg4iQZYIiIiIiIiTqIBloiIiIiIiJNogCUiIiIiIuIkGmCJiIiIiIg4iQZYIiIiIiIiTvL/ptUIwxLFEKIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x720 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.plot_posterior(trace, varnames=[l_diff_means, l_diff_sd, l_g1_mean, l_g1_sd, l_g2_mean, l_g2_sd, l_nu, l_effect_size], color='#87ceeb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now the estimated difference in means matches the empirical difference (`0.32` vs `0.33`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Same thing in BEST, which uses the $\\tau = (5\\sigma)^{-2}$ hyperparameter on the prior distributions for the groups' means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "*********************************************************************\n",
      "Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:\n",
      "A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.\n",
      "*********************************************************************\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: coda\n",
      "\n",
      "R[write to console]: Linked to JAGS 4.3.0\n",
      "\n",
      "R[write to console]: Loaded modules: basemod,bugs\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <span>ListVector with 2 elements.</span>\n",
       "    <table>\n",
       "      <tbody>\n",
       "      \n",
       "      <tr>\n",
       "      <th>\n",
       "        value\n",
       "      </th>\n",
       "      <td>\n",
       "        function (mu1, sd1, mu2, sd2, nPerGrp, pcntOut = 0, sdOutMult = 2, \n",
       "    rnd.seed = NULL) \n",
       "{\n",
       "    if (!is.null(rnd.seed)) {\n",
       "        set.seed(rnd.seed)\n",
       "    }\n",
       "    nOut = ceiling(nPerGrp * pcntOut/100)\n",
       "    nIn = nPerGrp - nOut\n",
       "    if (pcntOut > 100 | pcntOut < 0) {\n",
       "        stop(\"pcntOut must be between 0 and 100.\")\n",
       "    }\n",
       "    if (pcntOut > 0 & pcntOut < 1) {\n",
       "        warning(\"pcntOut is specified as percentage 0-100, not proportion 0-1.\")\n",
       "    }\n",
       "    if (pcntOut > 50) {\n",
       "        warning(\"pcntOut indicates more than 50% outliers; did you intend this?\")\n",
       "    }\n",
       "    if (nOut < 2 & pcntOut > 0) {\n",
       "        stop(\"Combination of nPerGrp and pcntOut yields too few outliers.\")\n",
       "    }\n",
       "    if (nIn < 2) {\n",
       "        stop(\"Too few non-outliers.\")\n",
       "    }\n",
       "    sdN = function(x) {\n",
       "        sqrt(mean((x - mean(x))^2))\n",
       "    }\n",
       "    genDat = function(mu, sigma, Nin, Nout, sigmaOutMult = sdOutMult) {\n",
       "        y = rnorm(n = Nin)\n",
       "        y = ((y - mean(y))/sdN(y)) * sigma + mu\n",
       "        yOut = NULL\n",
       "        if (Nout > 0) {\n",
       "            yOut = rnorm(n = Nout)\n",
       "            yOut = ((yOut - mean(yOut))/sdN(yOut)) * (sigma * \n",
       "                sdOutMult) + mu\n",
       "        }\n",
       "        y = c(y, yOut)\n",
       "        return(y)\n",
       "    }\n",
       "    y1 = genDat(mu = mu1, sigma = sd1, Nin = nIn, Nout = nOut)\n",
       "    y2 = genDat(mu = mu2, sigma = sd2, Nin = nIn, Nout = nOut)\n",
       "    openGraph(width = 7, height = 7)\n",
       "    layout(matrix(1:2, nrow = 2))\n",
       "    histInfo = hist(y1, main = \"Simulated Data\", col = \"pink2\", \n",
       "        border = \"white\", xlim = range(c(y1, y2)), breaks = 30, \n",
       "        prob = TRUE)\n",
       "    text(max(c(y1, y2)), max(histInfo$density), bquote(N == .(nPerGrp)), \n",
       "        adj = c(1, 1))\n",
       "    xlow = min(histInfo$breaks)\n",
       "    xhi = max(histInfo$breaks)\n",
       "    xcomb = seq(xlow, xhi, length = 1001)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu1, sd = sd1) * nIn/(nIn + \n",
       "        nOut) + dnorm(xcomb, mean = mu1, sd = sd1 * sdOutMult) * \n",
       "        nOut/(nIn + nOut), lwd = 3)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu1, sd = sd1) * nIn/(nIn + \n",
       "        nOut), lty = \"dashed\", col = \"green\", lwd = 3)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu1, sd = sd1 * sdOutMult) * \n",
       "        nOut/(nIn + nOut), lty = \"dashed\", col = \"red\", lwd = 3)\n",
       "    histInfo = hist(y2, main = \"\", col = \"pink2\", border = \"white\", \n",
       "        xlim = range(c(y1, y2)), breaks = 30, prob = TRUE)\n",
       "    text(max(c(y1, y2)), max(histInfo$density), bquote(N == .(nPerGrp)), \n",
       "        adj = c(1, 1))\n",
       "    xlow = min(histInfo$breaks)\n",
       "    xhi = max(histInfo$breaks)\n",
       "    xcomb = seq(xlow, xhi, length = 1001)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu2, sd = sd2) * nIn/(nIn + \n",
       "        nOut) + dnorm(xcomb, mean = mu2, sd = sd2 * sdOutMult) * \n",
       "        nOut/(nIn + nOut), lwd = 3)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu2, sd = sd2) * nIn/(nIn + \n",
       "        nOut), lty = \"dashed\", col = \"green\", lwd = 3)\n",
       "    lines(xcomb, dnorm(xcomb, mean = mu2, sd = sd2 * sdOutMult) * \n",
       "        nOut/(nIn + nOut), lty = \"dashed\", col = \"red\", lwd = 3)\n",
       "    return(list(y1 = y1, y2 = y2))\n",
       "}\n",
       "\n",
       "      </td>\n",
       "      </tr>\n",
       "      \n",
       "      <tr>\n",
       "      <th>\n",
       "        visible\n",
       "      </th>\n",
       "      <td>\n",
       "        \n",
       "    <span>BoolVector with 1 elements.</span>\n",
       "    <table>\n",
       "      <tbody>\n",
       "      <tr>\n",
       "      \n",
       "      <td>\n",
       "               0\n",
       "      </td>\n",
       "      \n",
       "      </tr>\n",
       "      </tbody>\n",
       "    </table>\n",
       "    \n",
       "      </td>\n",
       "      </tr>\n",
       "      \n",
       "      </tbody>\n",
       "    </table>\n",
       "    "
      ],
      "text/plain": [
       "R object with classes: ('list',) mapped to:\n",
       "[SignatureTranslatedFunc..., BoolVector]\n",
       "  value: <class 'rpy2.robjects.functions.SignatureTranslatedFunction'>\n",
       "  R object with classes: ('function',) mapped to:\n",
       "  visible: <class 'rpy2.robjects.vectors.BoolVector'>\n",
       "  R object with classes: ('RTYPES.LGLSXP',) mapped to:\n",
       "[       0]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import rpy2.robjects as robjects\n",
    "r_source = robjects.r['source']\n",
    "r_source(\"/Volumes/cristae/samadhi_daeda/Downloads/BEST/DBDA2E-utilities.R\")\n",
    "r_source(\"/Volumes/cristae/samadhi_daeda/Downloads/BEST/BEST.R\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2.robjects import numpy2ri\n",
    "numpy2ri.activate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "best = robjects.r['BESTmcmc']\n",
    "bestplot = robjects.r['BESTplot']\n",
    "bestshow = robjects.r['show']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"------------------\"\n",
      "[1] 3.450929\n",
      "[1] 0.1005448\n",
      "Calling 3 simulations using the parallel method...\n",
      "Following the progress of chain 1 (the program will wait for all chains\n",
      "to finish before continuing):\n",
      "Welcome to JAGS 4.3.0 on Thu Dec 13 23:52:09 2018\n",
      "JAGS is free software and comes with ABSOLUTELY NO WARRANTY\n",
      "Loading module: basemod: ok\n",
      "Loading module: bugs: ok\n",
      ". . Reading data file data.txt\n",
      ". Compiling model graph\n",
      "   Resolving undeclared variables\n",
      "   Allocating nodes\n",
      "Graph information:\n",
      "   Observed stochastic nodes: 25\n",
      "   Unobserved stochastic nodes: 5\n",
      "   Total graph size: 74\n",
      ". Reading parameter file inits1.txt\n",
      ". Initializing model\n",
      ". Adapting 500\n",
      "-------------------------------------------------| 500\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++ 100%\n",
      "Adaptation successful\n",
      ". Updating 1000\n",
      "-------------------------------------------------| 1000\n",
      "************************************************** 100%\n",
      ". . . . Updating 6667\n",
      "-------------------------------------------------| 6650\n",
      "************************************************** 100%\n",
      "* 100%\n",
      ". . . . Updating 0\n",
      ". Deleting model\n",
      ". \n",
      "All chains have finished\n",
      "Simulation complete.  Reading coda files...\n",
      "Coda files loaded successfully\n",
      "Finished running the simulation\n"
     ]
    }
   ],
   "source": [
    "mcmcChain = best(y1,y2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "postInfo = bestplot(y1,y2,mcmcChain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           SUMMARY.INFO\n",
      "PARAMETER         mean     median       mode     HDIlow    HDIhigh pcgtZero\n",
      "  mu1        3.6912955  3.6938500  3.6864866  3.2836500   4.074370       NA\n",
      "  mu2        3.3570506  3.3578000  3.3324278  2.9861700   3.732540       NA\n",
      "  muDiff     0.3342450  0.3361900  0.3456050 -0.2086800   0.861050 89.76051\n",
      "  sigma1     0.5059501  0.4678200  0.4039889  0.2319230   0.861527       NA\n",
      "  sigma2     0.7323208  0.7126010  0.6690985  0.4721920   1.028840       NA\n",
      "  sigmaDiff -0.2263707 -0.2391440 -0.2306295 -0.6788530   0.253770 13.30433\n",
      "  nu        37.8477728 28.5739000 12.5767150  1.9490500 101.604000       NA\n",
      "  nuLog10    1.4334549  1.4559695  1.4895861  0.7032793   2.138161       NA\n",
      "  effSz      0.5392534  0.5415953  0.6326409 -0.2933158   1.372787 89.76051\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<rpy2.rinterface.NULLType object at 0x12073e088> [RTYPES.NILSXP] [RTYPES.NILSXP]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bestshow( postInfo ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "muDiff is 0.33"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BESTmcmc uses the same $\\tau$ precision prior that the first example does and gets the means (and their difference right). What's going on with my pymc3 implementation such that I don't get the same behavior?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both https://github.com/mikemeredith/BEST and http://www.indiana.edu/~kruschke/BEST/ use this $\\tau$ parameterizion by default and both return the same results as above (correct estimates of $\\Delta \\mu$), whereas this pymc3 implimentation and https://docs.pymc.io/notebooks/BEST.html both return the underestimated mean difference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reference, here's the model code that JAGS runs in `BEST.R`, where \n",
    "\n",
    "`mu1PriorMean = mu2PriorMean = mean(y_joint)`\n",
    "\n",
    "`mu1PriorSD = mu2PriorSD = sd(y_joint)*5`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "  modelString = \"\"\"\n",
    "  model {\n",
    "    for ( i in 1:Ntotal ) {\n",
    "      y[i] ~ dt( mu[x[i]] , 1/sigma[x[i]]^2 , nu )\n",
    "    }\n",
    "    mu[1] ~ dnorm( mu1PriorMean , 1/mu1PriorSD^2 )  # prior for mu[1]\n",
    "    sigma[1] ~ dgamma( Sh1 , Ra1 )     # prior for sigma[1]\n",
    "    mu[2] ~ dnorm( mu2PriorMean , 1/mu2PriorSD^2 )  # prior for mu[2]\n",
    "    sigma[2] ~ dgamma( Sh2 , Ra2 )     # prior for sigma[2]\n",
    "    nu ~ dgamma( ShNu , RaNu ) # prior for nu\n",
    "  }\n",
    "  \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Rpy",
   "language": "python",
   "name": "rpy"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}