15_missing_data_and_other_opportunities-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ksachdeva/rethinking-tensorflow-probability/blob/master/notebooks/15_missing_data_and_other_opportunities.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 15. Missing Data and Other Opportunities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import arviz as az\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Ellipse, transforms\n",
    "from collections import namedtuple\n",
    "\n",
    "import random\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "import tensorflow as tf\n",
    "from tensorflow_probability.substrates import jax as tfp\n",
    "tfd = tfp.distributions\n",
    "tfb = tfp.bijectors\n",
    "# to use 'sample_distributions'\n",
    "import tensorflow_probability as xtfp\n",
    "xtfd = xtfp.distributions\n",
    "xtfb = xtfp.bijectors\n",
    "\n",
    "key = tfp.util.SeedStream(748, salt=\"ch-14\")\n",
    "plt.rcParams['figure.figsize'] = [6, 5]\n",
    "\n",
    "# TFP often requires 64-bit precision (float64) for certain computations \n",
    "# to ensure numerical stability and accuracy, e.g., tfd.Possion.\n",
    "# jax.config.update(\"jax_enable_x64\", True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 15.1 Measurement error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAGvCAYAAAAt2m1cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwGUlEQVR4nO3df3BU9b3/8dcC5geYLAYMScqvqCjG+LOCVVLkOvyIcgF11LFeqOhtb9WAIh0vILVgOxLUllrLGGqv32jLUJw7FQW1AbxALNcqSKSKWMAaflQT8RZMIj9CQs73D7prQrLZs7tn95zdz/MxszPs7uHkfTaw+9rPT59lWZYAAIBRerhdAAAASDwCAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAG6hXJwRUVFaqoqNDevXslSRdddJF+/OMf6/rrr1dLS4t+9KMf6fXXX9cnn3wiv9+vsWPHavHixSooKLD9M9ra2vTZZ58pKytLPp8voosBAMBklmWpqalJBQUF6tGj++/4vkj2AlizZo169uyp8847T5L0wgsv6Mknn9R7772ngQMH6pZbbtH3v/99XXrppTp8+LBmzZql1tZWvfvuu7aL//vf/65BgwbZPh4AAHR04MABDRw4sNtjIgoAXcnJydGTTz6pf//3f+/03NatWzVy5Ejt27dPgwcPtnW+hoYG9e3bVwcOHFB2dnYspQEAYJTGxkYNGjRIX375pfx+f7fHRtQF0N7Jkyf13//93zpy5IiuvvrqLo9paGiQz+dT3759Q56nublZzc3NwftNTU2SpOzsbAIAAABRsNOFHvEgwA8++EBnnnmm0tPTdc8992jVqlUqKirqdNzx48c1d+5c3XHHHd1+kJeXl8vv9wdvNP8DABB/EXcBnDhxQvv379eXX36pP/zhD/qv//ovVVdXdwgBLS0tuvXWW7V//35t2rSp2wBwegtAoPmioaGBFgAAACLQ2Ngov99v6zM05jEAY8eO1bnnnqtf//rXkk59+N9222365JNPtGHDBvXr1y+i80VSPAAA+Fokn6FRjwEIsCwr+A0+8OG/Z88ebdy4MeIPfwAAkBgRBYCHH35Y119/vQYNGqSmpiatXLlSmzZtUlVVlVpbW3XLLbeopqZGr776qk6ePKn6+npJp2YKpKWlxeUCAABA5CIKAJ9//rmmTZumuro6+f1+XXLJJaqqqtK4ceO0d+9erV69WpJ02WWXdfh7Gzdu1JgxY5yqGQAAxCjmMQBOYwwAAADRieQzlL0AAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQBOc/REq4bOfU1D576moyda3S4HAIC4IAAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGiigAVFRU6JJLLlF2drays7N19dVX649//GPwecuytHDhQhUUFCgzM1NjxozRhx9+6HjRAAAgNhEFgIEDB2rx4sV699139e677+q6667TlClTgh/yTzzxhJYsWaKlS5dq69atysvL07hx49TU1BSX4gEAQHQiCgCTJk3SDTfcoPPPP1/nn3++HnvsMZ155pl6++23ZVmWnnrqKc2fP18333yziouL9cILL+jo0aNasWJFvOoHAABRiHoMwMmTJ7Vy5UodOXJEV199tWpra1VfX6/x48cHj0lPT9e1116rt956K+R5mpub1djY2OEGAADiK+IA8MEHH+jMM89Uenq67rnnHq1atUpFRUWqr6+XJA0YMKDD8QMGDAg+15Xy8nL5/f7gbdCgQZGWBKCdoydaNXTuaxo69zUdPdHqdjkAPCriAHDBBRdo+/btevvtt3Xvvffqzjvv1M6dO4PP+3y+DsdbltXpsfbmzZunhoaG4O3AgQORlgQAACLUK9K/kJaWpvPOO0+SdOWVV2rr1q365S9/qTlz5kiS6uvrlZ+fHzz+4MGDnVoF2ktPT1d6enqkZQAAgBjEvA6AZVlqbm5WYWGh8vLytH79+uBzJ06cUHV1ta655ppYfwwAAHBQRC0ADz/8sK6//noNGjRITU1NWrlypTZt2qSqqir5fD7NmjVLixYt0rBhwzRs2DAtWrRIvXv31h133BGv+gEAQBQiCgCff/65pk2bprq6Ovn9fl1yySWqqqrSuHHjJEn/+Z//qWPHjum+++7T4cOHddVVV2ndunXKysqKS/EAACA6EQWA5557rtvnfT6fFi5cqIULF8ZSk6tOtlnBP2+pPaRvDztbPXuEHsQIAEAyYi+Adqp21Gnskurg/emVW1Xy+AZV7ahzsSoAAJxHAPinqh11und5jT5vbO7weH3Dcd27vIYQAABIKQQAnWr2f3TNTlldPBd47NE1Ozt0DwAAkMwIADrV11/XcDzk85akuobj2lJ7KHFFAQAQRwQASQebQn/4R3McAABeRwCQlJuV4ehxAAB4HQFA0sjCHOX7MxRqsp9PUr4/QyMLcxJZFgAAcUMAkNSzh08LJhVJUqcQELi/YFIR6wEAAFIGAeCfSovzVTH1CuVmd9yYKM+foYqpV6i0OD/E3wQAIPlEvBtgKistzteo8/rr4oXrJEnP3zWClQABACmJFoDTtP+wH1mYw4c/ACAlEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQCQEEdPtGro3Nc0dO5rOnqi1e1yAMB4BAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAMBwLNIEmIkAAACAgQgAAAAYiAAQBzSpAgC8jgAAAICBCABIiJNtVvDPW2oPdbgPAEg8AgDirmpHncYuqQ7en165VSWPb1DVjjoXqwIAsxEAEFdVO+p07/Iafd7Y3OHx+objund5DSEAAFxCAEDcnGyz9OianeqqsT/w2KNrdtIdAAAuIADEAf3dp2ypPaS6huMhn7ck1TUc15baQ4krCgAgiQDgOPq7v3awKfSHfzTHAQCcQwBwEP3dHeVmZTh6HADAOQQAh9Df3dnIwhzl+zPkC/G8T1K+P0MjC3MSWRYAQAQAx9Df3VnPHj4tmFQkSZ1CQOD+gklF6tkjVEQAAMQLAcAh9Hd3rbQ4XxVTr1BudnqHx/P8GaqYeoVKi/NdqgxuYalswBsIAA6x24/9wMrtxr3xlRbn643Z1wbvP3/XCG2ecx0f/gDgIgKAQ+z0d+dlmzvYrX0z/8jCHJr944hpqADsIAA4xE5/97wbhie0JpiHaagA7CIAOChcf/e4ogEuVQYTMA0VQCQIAA6jvxtuYBoqgEgRAOKA/m4kGtNQAUSKAACkAKahAogUAQBIASy7DCBSRgQAFh4xj2m/c5ZdBhApIwIAkOpYdhlApCIKAOXl5RoxYoSysrKUm5urG2+8Ubt27epwzFdffaUZM2Zo4MCByszM1IUXXqiKigpHiwbQGcsux5dprUpIfb0iObi6ulplZWUaMWKEWltbNX/+fI0fP147d+5Unz59JEkPPvigNm7cqOXLl2vo0KFat26d7rvvPhUUFGjKlClxuQiv6Z3WS3sXT3S7DBiotDhfo87rr4sXrpN0ahrqt4edzTd/AJ1EFACqqqo63K+srFRubq62bdum0aNHS5L+/Oc/684779SYMWMkSf/xH/+hX//613r33XeNCQCAm5iGCsCOmMYANDQ0SJJycr4eWFRSUqLVq1fr008/lWVZ2rhxo3bv3q0JEyZ0eY7m5mY1NjZ2uAEAgPiKOgBYlqXZs2erpKRExcXFwceffvppFRUVaeDAgUpLS1NpaameeeYZlZSUdHme8vJy+f3+4G3QoEHRlgS4KlQfMX3HALwooi6A9mbMmKH3339fmzdv7vD4008/rbffflurV6/WkCFD9Oabb+q+++5Tfn6+xo4d2+k88+bN0+zZs4P3GxsbCQEpiHERAOAtUQWAmTNnavXq1XrzzTc1cODA4OPHjh3Tww8/rFWrVmnixFNv9pdccom2b9+un/3sZ10GgPT0dKWnp3d6HAAAxE9EAcCyLM2cOVOrVq3Spk2bVFhY2OH5lpYWtbS0qEePjj0LPXv2VFtbW+zVAgAAR0QUAMrKyrRixQq98sorysrKUn19vSTJ7/crMzNT2dnZuvbaa/XQQw8pMzNTQ4YMUXV1tX77299qyZIlcbkAAABidfREq4p+vFaStPMnE9Q7Leoe8qQR0RUGFvQJTPELqKys1PTp0yVJK1eu1Lx58/Rv//ZvOnTokIYMGaLHHntM99xzjyMFAwCA2EXcBRBOXl6eKisroy4oHtrvgb6l9hALowAAjJfybRxVO+q0YPWHwfvTK7cq35+hBZOKulwaldHqAAATpPRmQFU76nTv8hp93tjc4fH6huO6d3mNqnbUuVQZAADuStkAcLLN0qNrdqqrTovAY4+u2dmhewAAAFOkbADYUntIdQ3HQz5vSaprOK4ttYcSVtPpYxEIHwAAt6RsADjYFPrDP5rjYlW1o05jl1QH70+v3KqSxzfQDQEAcEXKBoDcrAxHj4sFYxEAAF6TsgFgZGGO8v0ZCjXZzycp35+hkYU5IY5wBmMRAABelLIBoGcPnxZMKpKkTiEgcH/BpKK4rwfgxbEIAACkbACQpNLifFVMvUK52R03G8rzZ6hi6hVdrgPgNK+NRQAAQDJgIaDS4nyNOq+/Ll64TpL0/F0jEroSoJfGIgAAEJDSLQAB7T/sRxbmJHQZYK+MRUiUoydaNXTuaxo69zUdPdHqdjkAgBCMCABu8spYBAAA2iMAJIAXxiIAANBeyo8B8Aq3xyIkExP35QaARKMFIIHcHIuA+GOpZwDJhAAAOMCppZ4ZRAkgUQgAQIzCLfW8fufnLlUGAKERAIAY2Fnqufz1vyayJACwhQAAxMDOUs/1jazyCMB7CABADFjCGUCyIgAAMWAJZwDJyogJ1r3Temnv4olul2GE06fCpfpaB4Glnusbjnc5DsAnaUB2Bt0AQARYCyQxaAGAY5yaCpdM7Cz1PO+G4QmtCQDsIADAEeGmwqVyCAi31PO4ogEuVQYAoREAEDM7U+EeXbMzpVfGKy3O1xuzrw3ef/6uEdo85zr2eQDgWQQAxMzOVLi6huPaUnvI1vmSdUldlnoGkEwIAIiZ3alwdo7rbhwBy+QCgHMIAIiZ3alw4Y5jSV0ASBwCAGIWmAoXqsHbJynfn6GRhTkhz8GSugCQWAQAxMzOVLgFk4q67RN3ekndZB1H4AZeK8BMBAA4ItxUuHCj4Z1cUtfE9QiixWsFmIsAAMfEMhXOqSV1TV6PIFK8VoDZCABwVLRT4eyMI8jL7j4ksB6BfW6+VnQ5AN5AAIAnBMYRhPoosBR+SV2n1yNwSvspi9FOX3R6CqRbrxXTPAHvIAAgKNnfgJ1cjyDVufFaMc0T8BYCADwh0CQdik/hpwE6tR6BCRL9WjHNE/Ae9lhMILYlDs2JaYB2tubNC7MegSkS/Vo5Pc0TQOxoAYAnONHU3H49gtPZXY/AFE6s3RAJp7tdvN5d5fX6AIkAAI9wqqk5uB5BVnTrEZgk1rUbIpHK3S582CNZ0QUAT7DTJD0gO8NWM3Fpcb4uH9xXVy3aIElaNvUKjSvKM/6b/9ETrSr68VpJ0s6fTFDvtF4qLc7XqPP66+KF6ySdWrvh28POdvy1cvL3C8AZtADAE+w034ebBnj6+QKuHHpW2A+0eH6Laz/P/d29hz037z0R2xjb6XKI5PcLIHYEAHhGoEn67BDN9+OKBrhUWfSqdtRp0q/+N3j/nuU1xi61G67LIRl/v0AyowsAjop1pkNpcb6G52VpzM9OLRbz9O2XauIl31DPHr6k618NzHs//ft+YN67iWMSuutySLbfL5DsaAGA57Rvgr58cPjmey9iWeLQEtHlYCIGIyJStAAADmnf+vHnv/3D9lK7V5/bL0EVAsDXaAEA4oBliQF4HQEAnpOZ1rPLPycTliUG4HUEAHhO77ReXf45ntvIOn1uO9sb50ew1G4ybqGbrH3SyVo3YpOM/8diRQBAUuhuG1kvntvO9sZ2l9qN57UDMPf/GAEAnhduG9lY/pPG89xO8Hp9SDxaKJxl8v8xAgA8LZ7T6RJx7lB8Ns7NVEIgvkz/P0YAQJAX+8C27Ttsezpde3aW37WzRW1X57bDiXNHe+2hePH3C7gpnu8ByYB1ACDpVDPYgtUfBu9Pr9yqfH+GFkwqcnW1ui+amsMfpI7T6ap21OnHr3x9Lfcsr+nyWuI5Vc+Jc0dz7VLXqzF29/sdff7Ztn4OkGpMn65LCwA83Qd2+r4AoQSm0wWu5WBT+GuJ51Q9J84d6bWHEu73u37n57Z+DpBqTJ+uSwAwnNf7wL455Czb0+kivRanp+q1H5xV/I3smM8dybWHYuc1KX/9ryH/Pr7WfsAdg+9Sg9PvAcmGAGA4r/eB2dlGNjCdLtJrieTcsdR9OrvndqI+O69JfWNqNm8mCsEgecXzPSAZEAAMlwx9YOG2kQ3060dzLXbPHVPdIbY3tnPuWOtLpr5LPkjhhni+B3gdgwANlyx9YN1tIxsQ7bXYOXcsdV8+uK+uWrRBkrRs6hUaV5QX0bljqc/t3xuQDOL5HuBltAAYLpn6wMJtIxvLtcRzi9r257pyaHTbG0dbn53XJC+bkACYuE01AcBwqdQH5kS/e6qx85rMu2F44goC4BkEAKRUH5gT/e6pJtxrMq5oQMw/g+VpgeTDGABISq0+MCf63VNNd68JH9iAmWgBQFAq9YE50e+ealLpNWFZYyB2EQWA8vJyjRgxQllZWcrNzdWNN96oXbt2dTruo48+0uTJk+X3+5WVlaVvfetb2r9/v2NFAzDH6d0Lpm7dCjgtogBQXV2tsrIyvf3221q/fr1aW1s1fvx4HTlyJHjM3/72N5WUlGj48OHatGmT/vKXv+iRRx5RRgYjjeFdgfXz9y6eqN5p9Ix51fqdn0e1bDVjFFIbv9/oRPROV1VV1eF+ZWWlcnNztW3bNo0ePVqSNH/+fN1www164okngsedc845DpQKwHSLXv8o5LLGgS2WTR/vAdgV0xiAhoYGSVJOzql51W1tbXrttdd0/vnna8KECcrNzdVVV12ll19+OeQ5mpub1djY2OEGIPnFo5/+9G/+7bm1bLWdrafj9fMY/4BYRB0ALMvS7NmzVVJSouLiYknSwYMH9dVXX2nx4sUqLS3VunXrdNNNN+nmm29WdXV1l+cpLy+X3+8P3gYNGhRtSQA8ws1++kQuf1y1o06TfvW/wfv3LK9RyeMb4rbDIuMf4KSoA8CMGTP0/vvv6/e//33wsba2NknSlClT9OCDD+qyyy7T3Llz9a//+q9atmxZl+eZN2+eGhoagrcDBw5EWxIADwjXTx/v7YcTtfxxd1tPz1q5PW4/j22d4ZSoAsDMmTO1evVqbdy4UQMHDgw+3r9/f/Xq1UtFRR1XHrvwwgtDzgJIT09XdnZ2hxuA5NVdP70UevthO03pA7LT475stZ0m9nDbLLd/PNpm+vYD25qOt7CtMxwXUQCwLEszZszQSy+9pA0bNqiwsLDD82lpaRoxYkSnqYG7d+/WkCFDYq8WgOeF66fvavvhUE3ppzdtP3zDhZLit2y13Sb29/Z/2e02y+050Uy/bd9htnWG4yIKAGVlZVq+fLlWrFihrKws1dfXq76+XseOHQse89BDD+nFF1/Ub37zG3388cdaunSp1qxZo/vuu8/x4oFk0H5aIVMMO+uuKf30pu1xRQPitmx1uCb2DX89GHzs/74KHXK6Em6aYjhfNEX28wA7Ino3qqiokCSNGTOmw+OVlZWaPn26JOmmm27SsmXLVF5ervvvv18XXHCB/vCHP6ikpMSRggETBdYpSDXhmtJ96ty0HY9lq+3U8fN1u4OP9T8zvYsjQ4t1muLZWZH9PMCOiLsAuroFPvwD7r77bu3Zs0fHjh3T9u3bNWXKFCdrBlJWPKd4xbMlon2dZ/U+w/b2w9E2bTu9bPWW2kNh62jfMnD54L7dbrMc6hzRTlP85pCzbGzr/HVIYHog7GAvAMAjknWK1+l1Hz7a0uU36a62H/ZK03akUwcD2yxH8xEbzTTFcNt2W5KOtZwMPpYs/3bgLgIA4AHh+p+9+kYequ6udLX9sFeathM1dTCWnxVq2+6+vc+QJDUc67gErtf/7cB9BADAZeH6n6VTfcdea9Ltru7TPX/XCG2ec12nQXr2mrbj/+E8sjAnbB0D2n3wBq49Ek5MUywtztcbs68N3v9/d16p9F5dv417+d8OvIEAgJSUTCPv7fQ/O7HErdPjC8LV3V6ofvr2Tdun66rLIF7CNbFL0oPjhgUfW7nlgO1rb3+OWKcpSh3HP/To4VO9B5dHRnLw9jsjEsrrI829Xl+07PYJx7LEbdWOOi1Y/WHw/vTKrcr3Z2jBpCKNPv/sqM7p1JK7gabtH7/yYYepgHkx1hdtHQtWf9ihSyPPn6HJl+Zrybo9wcd+8caerk4RUuBaYpmm2BW7YygSuTwykgcBAHCZ3T7haPuOA/30p3/fD/QRP3X7ZVGd18l+89LifF0+uK+uWrRBkrRs6hXB6XKJ3N61qymGR5tPqmxF59fPLiemKYZidwxFIsc4IHnQBQC4zE7/c7R9x3bGF0S7hKyduiPR/gPyyqFnubalb/uf+80hZ+mnr9kb59Be+8qdmKYYip0xFE4sj4zURABASgp0F+xdPNHzYwDs9D9H23dsZ3xBtEvI2qk72YVbp6Aribz2eP7bQeojAAAeEGqKV6xL3Ma777e7uqPtWvCSaNYpSPS1x+vfDlKft78aAQaJxxK3iej7DVV3c+vJMH/T+yJdpyAwdiHR1x6PfztIfbQAAA5xYpqd00vc2umnd2KefVd1tx+8Z2cgn53tgBPNTh97+/UBvDJ2IZ7jDtzUfovkRA4OTVUEAMABXl3G104fcSLm2YdjdzvgRLOzTsEPx5+fuIKQMrwQZggAQIzCLePbfjtbN4TrI26/NK8bItkO2A3B1y+r69fvuuG5LlUGxIYxAEAMotnO1g3d9RG72ZSaTK9fqHUK/u8rFtlBcqIFAJ4Tzy1xnRbPaXZO82IfcbTbAbvBK+sUoLNkes/wEloA4CndLVnrxelMLLEaG69sB4zkFeo9Y+71w/XAyu2SpJ0/meD59UDcQAsAPCMZt8RlidXYeGU7YCSn7t4zZv3zwx+hEQDgCcm6JW6iptmlKq9sB4zkY+c9A90jAMATErUlrtOSZZqdG+zM6/fKdsBeFKpfO5mWuY6WnT59O+NH0D0CADwhEVvixovXp9m5IZJ5/eGm2Zn6+nlxXYlEsHvtjB+JHQEAnhDvLXHjrbQ4X2/MvjZ4//m7RmjznOs8OXAx3tbv/Dzief2lxflaM3NU8P6yqVcY+/ol41gYp0SypgbjR2JHAIAnxHNL3ETx4jQ7Nyx6/aOoth9mml3yjoVxQqRbV9sZP4LuEQDgCam8rakJfbaBa/z997/V6dtbe16a1+9FyToWxgmRrqlhwnbU8UYAgGewrWnysztG45e3XxZTIErVUJXMY2FiFc01pfp21PGWOv9zkBLY1jS5JftYDreZ/PpFe02pvB11vNECAM+hLz15eXUshxNLxSZiu2Kvvn6JEMuaGrxnRIcAAMAxXhzL4cSUukRtV+zF1y9RWFMj8QgAABzlpbEcTkypC7dd8Ya/HnS0Zi+9fonGmhqJxRgAAI7zwlgOO1sNP7pmZ3Bb34DAAEO751iybk+Hv+sEL7x+bvHq1tWpiAAAGK79B56T3O6XjWRK3dXn9ov6HPGa1uj26+cmk689kQgAAFKSE1PqIpnWOOWyb9g6FvAKxgAASElOTKkzeVoeUh8BAEBKcmJKncnT8pD6CACAIVJ19bxQnJhSZ2e74lSdludFTqznEIpp/z8kAgCAFObElLpw2xWn8rQ8LzF5i+R4MSPmADCWE1PqSovzdfngvrpq0QZJp7YrPn36IOInsBbD6d/3A2sxsO5/dAgAAOIiXtMLo+HEtDK2K3aHnbUYQm0xje7RBQAA8Kxt+w67thZDqqMFAAAc5qXWj2T3xWlLMMM5BAAASAGpGjrOPm3wJZxDFwAAwLO+OeSsqLcJRvcIAAAAz2Kb4PghAADwhPaLr5iyEAvsYZvg+OB/GQDA89gm2Hm0AAAAkgLbBDuLAAAAgIEIAAAAGIgAAABIOfHcOTBVMAgQAJDUTl8EqWpHnRas/jB4f3rlVuX7M7RgUhG7N7ZDAADgCam6kh0SK9zOgWzh/DUCAOAQPsAAd9nZOfDRNTvZyvmfCAAA4EEEyshtqT0UdufAuobj2lJ7SFef2y9xhXXh9DEKgTUNEolBgACAlHCwyd62wHaPi5eqHXUau6Q6eH965VaVPL5BVTvqEloHAQAAkBJys+xtCmT3uHgIjFH4vLHjNseBMQqJDAEEAABAShhZmBN258B8f4ZGFuYksqygcGMUpFNjFBI1ZZEAAMDTmM8Nu+zsHLhgUpFrAwAjGaOQCAQAAJ7llb5SJI9wOwe6OQXQa2MUCAAAPMlLfaVILqXF+Xpj9rXB+8/fNUKb51zn+vx/r41RIADAUUdPtGro3Nc0dO5rbNGJqHmtrxTJx4s7B3ptjAIBAIAntA+Pf9rzhaf6SgEneG2MAgEAgOd80dQc/iC5P58biJSXxiiwEiAAzzk7Kz38QXJ3PjcQrdLifI06r78uXrhO0qkxCqwECACSvjnkLE/1lQJO88IYBQIAAM/xWl8pkIoiCgDl5eUaMWKEsrKylJubqxtvvFG7du0KefwPfvAD+Xw+PfXUU7HWCcAwXuorBVJRRGMAqqurVVZWphEjRqi1tVXz58/X+PHjtXPnTvXp06fDsS+//LLeeecdFRQUOFowAHN4pa8U3sAOic6KKABUVVV1uF9ZWanc3Fxt27ZNo0ePDj7+6aefasaMGVq7dq0mTuz+l9Xc3Kzm5q9H/DY2NkZSEoAU54W+UiAVxTQGoKGhQZKUk/P1QJy2tjZNmzZNDz30kC666KKw5ygvL5ff7w/eBg0aFEtJAADAhqgDgGVZmj17tkpKSlRcXBx8/PHHH1evXr10//332zrPvHnz1NDQELwdOHAg2pKApBdo4ty7eKJ6pzFLF0D8RP0OM2PGDL3//vvavHlz8LFt27bpl7/8pWpqauTz2WumS09PV3q6vTm/AADAGVG1AMycOVOrV6/Wxo0bNXDgwODjf/rTn3Tw4EENHjxYvXr1Uq9evbRv3z798Ic/1NChQ52qGQAAxCiiFgDLsjRz5kytWrVKmzZtUmFhYYfnp02bprFjx3Z4bMKECZo2bZruuuuu2KsFgCgwehzoLKIAUFZWphUrVuiVV15RVlaW6uvrJUl+v1+ZmZnq16+f+vXr1+HvnHHGGcrLy9MFF1zgXNVIabxZA0D8RdQFUFFRoYaGBo0ZM0b5+fnB24svvhiv+gAAQBxE3AUQqb1790b8dwAAQHyxFwAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAgA3t92ZgnwakAgIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAA8ISTbV/vNrql9lCH+wCcRwAA4LqqHXUau6Q6eH965VaVPL5BVTvqXKwKSG0sZwXAVVU76nTv8hqd/n2/vuG47l1eo6duv8yNspDkeqf10t7FE90uw9NoAQDgmpNtlh5ds7PTh7+k4GPlr/81kSUBxiAAwFH04yISW2oPqa7heMjnLUn1jaGfBxA9AgAcQz8uInWwiQ93wC0EADgi0I/7eWNzh8cD/biEAHQlNyvD7RIAYxEAEDM7/biPrtlJdwA6GVmYo3x/hnwhnvdJyssmJADxQABAzOz049Y1HNeW2kOJKwpJoWcPnxZMKpKkTiEgcH/eDcMTWhNgCgIAYma3H5f+XnSltDhfFVOvUG52eofH8/wZqph6hcYVDXCpMiC1sQ4AYma3H5f+XoRSWpyvUef118UL10mSnr9rhL497Gz17OHT0ROtLlcHpCYCAGIW6Metbzje5TgAn059mxtZmJPo0pBEevb4uhNgZGFOh/swC4v4JAZdAIiZnX7cBZOKeEMHAA+hBQCOCPTjLlj9YYepgHn+DC2YVKTS4nwXqwNix7dSOMkL/54IAHBMd/24AABvoQsAjqIfFwCSAwEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwEAEAAAADEQAAADAQAQAAAAMRAAAAMBABAAAAAxEAAAAwUC+3CwAASeqd1kt7F090uwzAGLQAAABgIFoAAHgaLQNAfNACAACAgQgAAAAYiAAAAICBCAAAABiIAAAAgIGYBQBHMWIbAJIDLQAAABgoogBQXl6uESNGKCsrS7m5ubrxxhu1a9eu4PMtLS2aM2eOLr74YvXp00cFBQX67ne/q88++8zxwgEAQPQiCgDV1dUqKyvT22+/rfXr16u1tVXjx4/XkSNHJElHjx5VTU2NHnnkEdXU1Oill17S7t27NXny5LgUDwAAouOzLMuK9i9/8cUXys3NVXV1tUaPHt3lMVu3btXIkSO1b98+DR48uNPzzc3Nam5uDt5vbGzUoEGD1NDQoOzs7GhLAwDAOI2NjfL7/bY+Q2MaA9DQ0CBJysnJ6fYYn8+nvn37dvl8eXm5/H5/8DZo0KBYSgIAADZE3QJgWZamTJmiw4cP609/+lOXxxw/flwlJSUaPny4li9f3uUxtAAAAOCMSFoAop4GOGPGDL3//vvavHlzl8+3tLTo9ttvV1tbm5555pmQ50lPT1d6enq0ZQAAgChEFQBmzpyp1atX680339TAgQM7Pd/S0qLbbrtNtbW12rBhA9/kAQDwmIgCgGVZmjlzplatWqVNmzapsLCw0zGBD/89e/Zo48aN6tevn2PFAgAAZ0QUAMrKyrRixQq98sorysrKUn19vSTJ7/crMzNTra2tuuWWW1RTU6NXX31VJ0+eDB6Tk5OjtLQ0568AAABELKJBgD6fr8vHKysrNX36dO3du7fLVgFJ2rhxo8aMGRP2Z0QygAEAAHwtboMAw2WFoUOHhj0GAAC4j70AAAAwEAEAAAADEQAAADBQ1AsBxUtgDEFjY6PLlQAAkFwCn512xuN5LgA0NTVJEnsCAAAQpaamJvn9/m6PiWk3wHhoa2vTZ599pqysrJDTDrsT2EvgwIEDKTuNkGtMDVxj6jDhOrnG5GBZlpqamlRQUKAePbrv5fdcC0CPHj26XF44UtnZ2Un7C7SLa0wNXGPqMOE6uUbvC/fNP4BBgAAAGIgAAACAgVIuAKSnp2vBggUpvcUw15gauMbUYcJ1co2px3ODAAEAQPylXAsAAAAIjwAAAICBCAAAABiIAAAAgIEIAAAAGChlAsCnn36qqVOnql+/furdu7cuu+wybdu2ze2yHDV06FD5fL5Ot7KyMrdLc0xra6t+9KMfqbCwUJmZmTrnnHP0k5/8RG1tbW6X5qimpibNmjVLQ4YMUWZmpq655hpt3brV7bKi9uabb2rSpEkqKCiQz+fTyy+/3OF5y7K0cOFCFRQUKDMzU2PGjNGHH37oTrFRCneNL730kiZMmKD+/fvL5/Np+/btrtQZi+6usaWlRXPmzNHFF1+sPn36qKCgQN/97nf12WefuVdwFML9HhcuXKjhw4erT58+OuusszR27Fi988477hQbZykRAA4fPqxRo0bpjDPO0B//+Eft3LlTP//5z9W3b1+3S3PU1q1bVVdXF7ytX79eknTrrbe6XJlzHn/8cS1btkxLly7VRx99pCeeeEJPPvmkfvWrX7ldmqO+973vaf369frd736nDz74QOPHj9fYsWP16aeful1aVI4cOaJLL71US5cu7fL5J554QkuWLNHSpUu1detW5eXlady4ccHNv5JBuGs8cuSIRo0apcWLFye4Mud0d41Hjx5VTU2NHnnkEdXU1Oill17S7t27NXnyZBcqjV643+P555+vpUuX6oMPPtDmzZs1dOhQjR8/Xl988UWCK00AKwXMmTPHKikpcbuMhHvggQesc88912pra3O7FMdMnDjRuvvuuzs8dvPNN1tTp051qSLnHT161OrZs6f16quvdnj80ksvtebPn+9SVc6RZK1atSp4v62tzcrLy7MWL14cfOz48eOW3++3li1b5kKFsTv9Gturra21JFnvvfdeQmtyWnfXGLBlyxZLkrVv377EFOUwO9fY0NBgSbLeeOONxBSVQCnRArB69WpdeeWVuvXWW5Wbm6vLL79cv/nNb9wuK65OnDih5cuX6+67745q10SvKikp0f/8z/9o9+7dkqS//OUv2rx5s2644QaXK3NOa2urTp48qYyMjA6PZ2ZmavPmzS5VFT+1tbWqr6/X+PHjg4+lp6fr2muv1VtvveViZYhVQ0ODfD5fyrW2Bpw4cULPPvus/H6/Lr30UrfLcVxKBIBPPvlEFRUVGjZsmNauXat77rlH999/v37729+6XVrcvPzyy/ryyy81ffp0t0tx1Jw5c/Sd73xHw4cP1xlnnKHLL79cs2bN0ne+8x23S3NMVlaWrr76av30pz/VZ599ppMnT2r58uV65513VFdX53Z5jquvr5ckDRgwoMPjAwYMCD6H5HP8+HHNnTtXd9xxR1LvnNeVV199VWeeeaYyMjL0i1/8QuvXr1f//v3dLstxKREA2tradMUVV2jRokW6/PLL9YMf/EDf//73VVFR4XZpcfPcc8/p+uuvV0FBgdulOOrFF1/U8uXLtWLFCtXU1OiFF17Qz372M73wwgtul+ao3/3ud7IsS9/4xjeUnp6up59+WnfccYd69uzpdmlxc3pLlWVZKdV6ZZKWlhbdfvvtamtr0zPPPON2OY77l3/5F23fvl1vvfWWSktLddttt+ngwYNul+W4lAgA+fn5Kioq6vDYhRdeqP3797tUUXzt27dPb7zxhr73ve+5XYrjHnroIc2dO1e33367Lr74Yk2bNk0PPvigysvL3S7NUeeee66qq6v11Vdf6cCBA9qyZYtaWlpUWFjodmmOy8vLk6RO3/YPHjzYqVUA3tfS0qLbbrtNtbW1Wr9+fcp9+5ekPn366LzzztO3vvUtPffcc+rVq5eee+45t8tyXEoEgFGjRmnXrl0dHtu9e7eGDBniUkXxVVlZqdzcXE2cONHtUhx39OhR9ejR8Z9lz549U24aYECfPn2Un5+vw4cPa+3atZoyZYrbJTmusLBQeXl5wVkr0qm+1erqal1zzTUuVoZIBT789+zZozfeeEP9+vVzu6SEsCxLzc3NbpfhuF5uF+CEBx98UNdcc40WLVqk2267TVu2bNGzzz6rZ5991u3SHNfW1qbKykrdeeed6tUrJX59HUyaNEmPPfaYBg8erIsuukjvvfeelixZorvvvtvt0hy1du1aWZalCy64QB9//LEeeughXXDBBbrrrrvcLi0qX331lT7++OPg/draWm3fvl05OTkaPHiwZs2apUWLFmnYsGEaNmyYFi1apN69e+uOO+5wserIhLvGQ4cOaf/+/cF58YEvJXl5ecFWEK/r7hoLCgp0yy23qKamRq+++qpOnjwZbNXJyclRWlqaW2VHpLtr7Nevnx577DFNnjxZ+fn5+sc//qFnnnlGf//731NqunWQu5MQnLNmzRqruLjYSk9Pt4YPH249++yzbpcUF2vXrrUkWbt27XK7lLhobGy0HnjgAWvw4MFWRkaGdc4551jz58+3mpub3S7NUS+++KJ1zjnnWGlpaVZeXp5VVlZmffnll26XFbWNGzdakjrd7rzzTsuyTk0FXLBggZWXl2elp6dbo0ePtj744AN3i45QuGusrKzs8vkFCxa4WnckurvGwPTGrm4bN250u3TburvGY8eOWTfddJNVUFBgpaWlWfn5+dbkyZOtLVu2uF12XPgsy7LiHzMAAICXpMQYAAAAEBkCAAAABiIAAABgIAIAAAAGIgAAAGAgAgAAAAYiAAAAYCACAAAABiIAAABgIAIAAAAGIgAAAGCg/w+iLwu3IW6oYgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = pd.read_csv(\"./data/WaffleDivorce.csv\", sep=\";\")\n",
    "\n",
    "plt.errorbar(x=d['Divorce'], y=d['MedianAgeMarriage'], yerr=d['Divorce SE'], fmt='o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above plot, the lenght of the vertical lines show how uncertain the observed divorce rate is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 15.1.1 Error on the outcome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Code 15.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "dlist = dict(\n",
    "  D_obs=d.Divorce.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(np.float32),\n",
    "  D_sd=(d[\"Divorce SE\"].values / d.Divorce.std()).astype(np.float32),\n",
    "  M=d.Marriage.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(np.float32),\n",
    "  A=d.MedianAgeMarriage.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(np.float32)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_15_1(A, M, D_sd):\n",
    "\tdef _generator():\n",
    "\t\ta = yield tfd.Normal(loc=0.0, scale=0.2, name=\"a\")\n",
    "\t\tbA = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bA\")\n",
    "\t\tbM = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bM\")\n",
    "\t\tsigma = yield tfd.Exponential(rate=1.0, name=\"sigma\")\n",
    "\t\tmu = a + bA * A + bM * M\n",
    "\t\tD_true = yield tfd.Normal(loc=mu, scale=sigma, name='D_true')\n",
    "\t\tyield tfd.Normal(loc=D_true, scale=D_sd, name='D_obs')\n",
    "\n",
    "\treturn tfd.JointDistributionCoroutineAutoBatched(_generator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "jd_15_1 = model_15_1(dlist[\"A\"], dlist[\"M\"], dlist[\"D_sd\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "jd_15_1_pin = jd_15_1.experimental_pin(D_obs=dlist[\"D_obs\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jax.jit\n",
    "def run_chain(init_state, seed, step_size=1e-3, num_steps=500, burnin=1000):\n",
    "\n",
    "  def target_log_prob(*states):\n",
    "    return jd_15_1_pin.unnormalized_log_prob(states)\n",
    "  \n",
    "  bijectors = [tfb.Identity(), tfb.Identity(), tfb.Identity(),\n",
    "               tfb.Exp(), tfb.Identity()]\n",
    "\n",
    "  kernel = tfp.mcmc.NoUTurnSampler(\n",
    "      target_log_prob_fn=target_log_prob,\n",
    "      step_size=step_size)\n",
    "  \n",
    "  kernel = tfp.mcmc.TransformedTransitionKernel(\n",
    "    inner_kernel=kernel,\n",
    "    bijector=bijectors)\n",
    "\n",
    "  kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(\n",
    "    inner_kernel=kernel,\n",
    "    num_adaptation_steps=burnin)\n",
    "\n",
    "  chain_state = tfp.mcmc.sample_chain(\n",
    "      num_results=num_steps,\n",
    "      num_burnin_steps=burnin,\n",
    "      current_state=init_state,\n",
    "      kernel=kernel,\n",
    "      seed=seed,\n",
    "      trace_fn=None)\n",
    "  \n",
    "  return chain_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_chains = 4\n",
    "init_state = jd_15_1_pin.sample_unpinned(num_chains, seed=key())\n",
    "post_1 = run_chain(init_state=init_state, seed=key())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "post_1_dict = {k: jnp.moveaxis(v, 0, 1) for k, v in post_1._asdict().items()}\n",
    "dict_1 = dict(posterior=post_1_dict)\n",
    "idata_1 = az.from_dict(**dict_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_5.5%</th>\n",
       "      <th>hdi_94.5%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>D_true[0]</th>\n",
       "      <td>1.159</td>\n",
       "      <td>0.369</td>\n",
       "      <td>0.604</td>\n",
       "      <td>1.763</td>\n",
       "      <td>0.011</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1141.0</td>\n",
       "      <td>496.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[1]</th>\n",
       "      <td>0.690</td>\n",
       "      <td>0.532</td>\n",
       "      <td>-0.107</td>\n",
       "      <td>1.597</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.013</td>\n",
       "      <td>885.0</td>\n",
       "      <td>1107.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[2]</th>\n",
       "      <td>0.433</td>\n",
       "      <td>0.344</td>\n",
       "      <td>-0.140</td>\n",
       "      <td>0.936</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.006</td>\n",
       "      <td>2302.0</td>\n",
       "      <td>1486.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[3]</th>\n",
       "      <td>1.404</td>\n",
       "      <td>0.475</td>\n",
       "      <td>0.659</td>\n",
       "      <td>2.174</td>\n",
       "      <td>0.012</td>\n",
       "      <td>0.009</td>\n",
       "      <td>1480.0</td>\n",
       "      <td>1006.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[4]</th>\n",
       "      <td>-0.899</td>\n",
       "      <td>0.135</td>\n",
       "      <td>-1.120</td>\n",
       "      <td>-0.703</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2010.0</td>\n",
       "      <td>1117.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[5]</th>\n",
       "      <td>0.662</td>\n",
       "      <td>0.408</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1.289</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1728.0</td>\n",
       "      <td>1169.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[6]</th>\n",
       "      <td>-1.365</td>\n",
       "      <td>0.337</td>\n",
       "      <td>-1.917</td>\n",
       "      <td>-0.875</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2438.0</td>\n",
       "      <td>1446.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[7]</th>\n",
       "      <td>-0.336</td>\n",
       "      <td>0.479</td>\n",
       "      <td>-1.130</td>\n",
       "      <td>0.389</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.010</td>\n",
       "      <td>1335.0</td>\n",
       "      <td>1110.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[8]</th>\n",
       "      <td>-1.899</td>\n",
       "      <td>0.593</td>\n",
       "      <td>-2.742</td>\n",
       "      <td>-0.868</td>\n",
       "      <td>0.019</td>\n",
       "      <td>0.013</td>\n",
       "      <td>1019.0</td>\n",
       "      <td>1305.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[9]</th>\n",
       "      <td>-0.618</td>\n",
       "      <td>0.171</td>\n",
       "      <td>-0.895</td>\n",
       "      <td>-0.358</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>2454.0</td>\n",
       "      <td>1146.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[10]</th>\n",
       "      <td>0.766</td>\n",
       "      <td>0.291</td>\n",
       "      <td>0.282</td>\n",
       "      <td>1.195</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2642.0</td>\n",
       "      <td>1383.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[11]</th>\n",
       "      <td>-0.553</td>\n",
       "      <td>0.472</td>\n",
       "      <td>-1.316</td>\n",
       "      <td>0.193</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.010</td>\n",
       "      <td>1313.0</td>\n",
       "      <td>1407.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[12]</th>\n",
       "      <td>0.185</td>\n",
       "      <td>0.490</td>\n",
       "      <td>-0.583</td>\n",
       "      <td>0.978</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.013</td>\n",
       "      <td>731.0</td>\n",
       "      <td>834.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[13]</th>\n",
       "      <td>-0.875</td>\n",
       "      <td>0.241</td>\n",
       "      <td>-1.205</td>\n",
       "      <td>-0.459</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3598.0</td>\n",
       "      <td>1239.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[14]</th>\n",
       "      <td>0.554</td>\n",
       "      <td>0.306</td>\n",
       "      <td>0.062</td>\n",
       "      <td>1.029</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2818.0</td>\n",
       "      <td>1316.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[15]</th>\n",
       "      <td>0.282</td>\n",
       "      <td>0.369</td>\n",
       "      <td>-0.285</td>\n",
       "      <td>0.886</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.007</td>\n",
       "      <td>2083.0</td>\n",
       "      <td>1523.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[16]</th>\n",
       "      <td>0.509</td>\n",
       "      <td>0.414</td>\n",
       "      <td>-0.161</td>\n",
       "      <td>1.144</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1825.0</td>\n",
       "      <td>1201.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[17]</th>\n",
       "      <td>1.248</td>\n",
       "      <td>0.344</td>\n",
       "      <td>0.681</td>\n",
       "      <td>1.778</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.006</td>\n",
       "      <td>1485.0</td>\n",
       "      <td>834.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[18]</th>\n",
       "      <td>0.424</td>\n",
       "      <td>0.377</td>\n",
       "      <td>-0.134</td>\n",
       "      <td>1.055</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.007</td>\n",
       "      <td>1810.0</td>\n",
       "      <td>1285.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[19]</th>\n",
       "      <td>0.388</td>\n",
       "      <td>0.540</td>\n",
       "      <td>-0.497</td>\n",
       "      <td>1.208</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.013</td>\n",
       "      <td>870.0</td>\n",
       "      <td>1224.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[20]</th>\n",
       "      <td>-0.565</td>\n",
       "      <td>0.322</td>\n",
       "      <td>-1.070</td>\n",
       "      <td>-0.047</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.006</td>\n",
       "      <td>1853.0</td>\n",
       "      <td>1011.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[21]</th>\n",
       "      <td>-1.097</td>\n",
       "      <td>0.265</td>\n",
       "      <td>-1.555</td>\n",
       "      <td>-0.698</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>2924.0</td>\n",
       "      <td>1406.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[22]</th>\n",
       "      <td>-0.274</td>\n",
       "      <td>0.258</td>\n",
       "      <td>-0.677</td>\n",
       "      <td>0.131</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2380.0</td>\n",
       "      <td>1176.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[23]</th>\n",
       "      <td>-1.002</td>\n",
       "      <td>0.305</td>\n",
       "      <td>-1.509</td>\n",
       "      <td>-0.542</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2681.0</td>\n",
       "      <td>1147.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[24]</th>\n",
       "      <td>0.429</td>\n",
       "      <td>0.401</td>\n",
       "      <td>-0.190</td>\n",
       "      <td>1.090</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.007</td>\n",
       "      <td>1730.0</td>\n",
       "      <td>1251.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[25]</th>\n",
       "      <td>-0.024</td>\n",
       "      <td>0.327</td>\n",
       "      <td>-0.552</td>\n",
       "      <td>0.449</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.010</td>\n",
       "      <td>2236.0</td>\n",
       "      <td>1430.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[26]</th>\n",
       "      <td>-0.018</td>\n",
       "      <td>0.503</td>\n",
       "      <td>-0.743</td>\n",
       "      <td>0.854</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.012</td>\n",
       "      <td>1402.0</td>\n",
       "      <td>1306.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[27]</th>\n",
       "      <td>-0.157</td>\n",
       "      <td>0.396</td>\n",
       "      <td>-0.824</td>\n",
       "      <td>0.427</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.009</td>\n",
       "      <td>1439.0</td>\n",
       "      <td>1404.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[28]</th>\n",
       "      <td>-0.252</td>\n",
       "      <td>0.504</td>\n",
       "      <td>-0.998</td>\n",
       "      <td>0.591</td>\n",
       "      <td>0.014</td>\n",
       "      <td>0.011</td>\n",
       "      <td>1375.0</td>\n",
       "      <td>1155.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[29]</th>\n",
       "      <td>-1.792</td>\n",
       "      <td>0.237</td>\n",
       "      <td>-2.147</td>\n",
       "      <td>-1.400</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>2397.0</td>\n",
       "      <td>1183.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[30]</th>\n",
       "      <td>0.168</td>\n",
       "      <td>0.423</td>\n",
       "      <td>-0.517</td>\n",
       "      <td>0.818</td>\n",
       "      <td>0.011</td>\n",
       "      <td>0.010</td>\n",
       "      <td>1385.0</td>\n",
       "      <td>1093.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[31]</th>\n",
       "      <td>-1.661</td>\n",
       "      <td>0.158</td>\n",
       "      <td>-1.899</td>\n",
       "      <td>-1.406</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2271.0</td>\n",
       "      <td>1119.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[32]</th>\n",
       "      <td>0.108</td>\n",
       "      <td>0.248</td>\n",
       "      <td>-0.275</td>\n",
       "      <td>0.511</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.005</td>\n",
       "      <td>3136.0</td>\n",
       "      <td>1252.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[33]</th>\n",
       "      <td>-0.052</td>\n",
       "      <td>0.515</td>\n",
       "      <td>-0.874</td>\n",
       "      <td>0.778</td>\n",
       "      <td>0.016</td>\n",
       "      <td>0.013</td>\n",
       "      <td>1080.0</td>\n",
       "      <td>1228.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[34]</th>\n",
       "      <td>-0.126</td>\n",
       "      <td>0.231</td>\n",
       "      <td>-0.468</td>\n",
       "      <td>0.245</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.006</td>\n",
       "      <td>3413.0</td>\n",
       "      <td>1289.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[35]</th>\n",
       "      <td>1.266</td>\n",
       "      <td>0.417</td>\n",
       "      <td>0.606</td>\n",
       "      <td>1.923</td>\n",
       "      <td>0.011</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1554.0</td>\n",
       "      <td>1343.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[36]</th>\n",
       "      <td>0.218</td>\n",
       "      <td>0.355</td>\n",
       "      <td>-0.314</td>\n",
       "      <td>0.815</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.006</td>\n",
       "      <td>1883.0</td>\n",
       "      <td>1109.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[37]</th>\n",
       "      <td>-1.028</td>\n",
       "      <td>0.220</td>\n",
       "      <td>-1.374</td>\n",
       "      <td>-0.676</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>2828.0</td>\n",
       "      <td>994.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[38]</th>\n",
       "      <td>-0.933</td>\n",
       "      <td>0.557</td>\n",
       "      <td>-1.807</td>\n",
       "      <td>-0.089</td>\n",
       "      <td>0.016</td>\n",
       "      <td>0.012</td>\n",
       "      <td>1167.0</td>\n",
       "      <td>1124.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[39]</th>\n",
       "      <td>-0.679</td>\n",
       "      <td>0.318</td>\n",
       "      <td>-1.169</td>\n",
       "      <td>-0.171</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2258.0</td>\n",
       "      <td>1281.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[40]</th>\n",
       "      <td>0.219</td>\n",
       "      <td>0.545</td>\n",
       "      <td>-0.650</td>\n",
       "      <td>1.066</td>\n",
       "      <td>0.015</td>\n",
       "      <td>0.011</td>\n",
       "      <td>1393.0</td>\n",
       "      <td>1202.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[41]</th>\n",
       "      <td>0.749</td>\n",
       "      <td>0.342</td>\n",
       "      <td>0.179</td>\n",
       "      <td>1.261</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.006</td>\n",
       "      <td>2178.0</td>\n",
       "      <td>1365.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[42]</th>\n",
       "      <td>0.201</td>\n",
       "      <td>0.188</td>\n",
       "      <td>-0.087</td>\n",
       "      <td>0.498</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2264.0</td>\n",
       "      <td>1244.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[43]</th>\n",
       "      <td>0.807</td>\n",
       "      <td>0.414</td>\n",
       "      <td>0.143</td>\n",
       "      <td>1.436</td>\n",
       "      <td>0.012</td>\n",
       "      <td>0.009</td>\n",
       "      <td>1127.0</td>\n",
       "      <td>1535.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[44]</th>\n",
       "      <td>-0.423</td>\n",
       "      <td>0.519</td>\n",
       "      <td>-1.244</td>\n",
       "      <td>0.406</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.011</td>\n",
       "      <td>1471.0</td>\n",
       "      <td>1066.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[45]</th>\n",
       "      <td>-0.386</td>\n",
       "      <td>0.259</td>\n",
       "      <td>-0.793</td>\n",
       "      <td>0.023</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2805.0</td>\n",
       "      <td>1113.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[46]</th>\n",
       "      <td>0.132</td>\n",
       "      <td>0.296</td>\n",
       "      <td>-0.324</td>\n",
       "      <td>0.593</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.006</td>\n",
       "      <td>2521.0</td>\n",
       "      <td>1081.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[47]</th>\n",
       "      <td>0.555</td>\n",
       "      <td>0.495</td>\n",
       "      <td>-0.219</td>\n",
       "      <td>1.367</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.009</td>\n",
       "      <td>1425.0</td>\n",
       "      <td>1277.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[48]</th>\n",
       "      <td>-0.638</td>\n",
       "      <td>0.284</td>\n",
       "      <td>-1.084</td>\n",
       "      <td>-0.190</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2448.0</td>\n",
       "      <td>1006.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[49]</th>\n",
       "      <td>0.855</td>\n",
       "      <td>0.586</td>\n",
       "      <td>-0.050</td>\n",
       "      <td>1.800</td>\n",
       "      <td>0.020</td>\n",
       "      <td>0.015</td>\n",
       "      <td>869.0</td>\n",
       "      <td>1121.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>-0.055</td>\n",
       "      <td>0.100</td>\n",
       "      <td>-0.216</td>\n",
       "      <td>0.097</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1653.0</td>\n",
       "      <td>1639.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bA</th>\n",
       "      <td>-0.615</td>\n",
       "      <td>0.156</td>\n",
       "      <td>-0.869</td>\n",
       "      <td>-0.368</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1313.0</td>\n",
       "      <td>1604.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bM</th>\n",
       "      <td>0.054</td>\n",
       "      <td>0.160</td>\n",
       "      <td>-0.198</td>\n",
       "      <td>0.314</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1179.0</td>\n",
       "      <td>1620.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.585</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.418</td>\n",
       "      <td>0.765</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>447.0</td>\n",
       "      <td>413.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             mean     sd  hdi_5.5%  hdi_94.5%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "D_true[0]   1.159  0.369     0.604      1.763      0.011    0.008    1141.0   \n",
       "D_true[1]   0.690  0.532    -0.107      1.597      0.018    0.013     885.0   \n",
       "D_true[2]   0.433  0.344    -0.140      0.936      0.007    0.006    2302.0   \n",
       "D_true[3]   1.404  0.475     0.659      2.174      0.012    0.009    1480.0   \n",
       "D_true[4]  -0.899  0.135    -1.120     -0.703      0.003    0.002    2010.0   \n",
       "D_true[5]   0.662  0.408     0.008      1.289      0.010    0.008    1728.0   \n",
       "D_true[6]  -1.365  0.337    -1.917     -0.875      0.007    0.005    2438.0   \n",
       "D_true[7]  -0.336  0.479    -1.130      0.389      0.013    0.010    1335.0   \n",
       "D_true[8]  -1.899  0.593    -2.742     -0.868      0.019    0.013    1019.0   \n",
       "D_true[9]  -0.618  0.171    -0.895     -0.358      0.003    0.003    2454.0   \n",
       "D_true[10]  0.766  0.291     0.282      1.195      0.006    0.004    2642.0   \n",
       "D_true[11] -0.553  0.472    -1.316      0.193      0.013    0.010    1313.0   \n",
       "D_true[12]  0.185  0.490    -0.583      0.978      0.018    0.013     731.0   \n",
       "D_true[13] -0.875  0.241    -1.205     -0.459      0.004    0.003    3598.0   \n",
       "D_true[14]  0.554  0.306     0.062      1.029      0.006    0.005    2818.0   \n",
       "D_true[15]  0.282  0.369    -0.285      0.886      0.008    0.007    2083.0   \n",
       "D_true[16]  0.509  0.414    -0.161      1.144      0.010    0.008    1825.0   \n",
       "D_true[17]  1.248  0.344     0.681      1.778      0.009    0.006    1485.0   \n",
       "D_true[18]  0.424  0.377    -0.134      1.055      0.009    0.007    1810.0   \n",
       "D_true[19]  0.388  0.540    -0.497      1.208      0.018    0.013     870.0   \n",
       "D_true[20] -0.565  0.322    -1.070     -0.047      0.007    0.006    1853.0   \n",
       "D_true[21] -1.097  0.265    -1.555     -0.698      0.005    0.003    2924.0   \n",
       "D_true[22] -0.274  0.258    -0.677      0.131      0.005    0.005    2380.0   \n",
       "D_true[23] -1.002  0.305    -1.509     -0.542      0.006    0.004    2681.0   \n",
       "D_true[24]  0.429  0.401    -0.190      1.090      0.010    0.007    1730.0   \n",
       "D_true[25] -0.024  0.327    -0.552      0.449      0.007    0.010    2236.0   \n",
       "D_true[26] -0.018  0.503    -0.743      0.854      0.013    0.012    1402.0   \n",
       "D_true[27] -0.157  0.396    -0.824      0.427      0.010    0.009    1439.0   \n",
       "D_true[28] -0.252  0.504    -0.998      0.591      0.014    0.011    1375.0   \n",
       "D_true[29] -1.792  0.237    -2.147     -1.400      0.005    0.003    2397.0   \n",
       "D_true[30]  0.168  0.423    -0.517      0.818      0.011    0.010    1385.0   \n",
       "D_true[31] -1.661  0.158    -1.899     -1.406      0.003    0.002    2271.0   \n",
       "D_true[32]  0.108  0.248    -0.275      0.511      0.004    0.005    3136.0   \n",
       "D_true[33] -0.052  0.515    -0.874      0.778      0.016    0.013    1080.0   \n",
       "D_true[34] -0.126  0.231    -0.468      0.245      0.004    0.006    3413.0   \n",
       "D_true[35]  1.266  0.417     0.606      1.923      0.011    0.008    1554.0   \n",
       "D_true[36]  0.218  0.355    -0.314      0.815      0.008    0.006    1883.0   \n",
       "D_true[37] -1.028  0.220    -1.374     -0.676      0.004    0.003    2828.0   \n",
       "D_true[38] -0.933  0.557    -1.807     -0.089      0.016    0.012    1167.0   \n",
       "D_true[39] -0.679  0.318    -1.169     -0.171      0.007    0.005    2258.0   \n",
       "D_true[40]  0.219  0.545    -0.650      1.066      0.015    0.011    1393.0   \n",
       "D_true[41]  0.749  0.342     0.179      1.261      0.007    0.006    2178.0   \n",
       "D_true[42]  0.201  0.188    -0.087      0.498      0.004    0.004    2264.0   \n",
       "D_true[43]  0.807  0.414     0.143      1.436      0.012    0.009    1127.0   \n",
       "D_true[44] -0.423  0.519    -1.244      0.406      0.013    0.011    1471.0   \n",
       "D_true[45] -0.386  0.259    -0.793      0.023      0.005    0.004    2805.0   \n",
       "D_true[46]  0.132  0.296    -0.324      0.593      0.006    0.006    2521.0   \n",
       "D_true[47]  0.555  0.495    -0.219      1.367      0.013    0.009    1425.0   \n",
       "D_true[48] -0.638  0.284    -1.084     -0.190      0.006    0.005    2448.0   \n",
       "D_true[49]  0.855  0.586    -0.050      1.800      0.020    0.015     869.0   \n",
       "a          -0.055  0.100    -0.216      0.097      0.002    0.002    1653.0   \n",
       "bA         -0.615  0.156    -0.869     -0.368      0.004    0.003    1313.0   \n",
       "bM          0.054  0.160    -0.198      0.314      0.005    0.003    1179.0   \n",
       "sigma       0.585  0.110     0.418      0.765      0.005    0.004     447.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "D_true[0]      496.0   1.00  \n",
       "D_true[1]     1107.0   1.00  \n",
       "D_true[2]     1486.0   1.00  \n",
       "D_true[3]     1006.0   1.00  \n",
       "D_true[4]     1117.0   1.00  \n",
       "D_true[5]     1169.0   1.01  \n",
       "D_true[6]     1446.0   1.00  \n",
       "D_true[7]     1110.0   1.00  \n",
       "D_true[8]     1305.0   1.00  \n",
       "D_true[9]     1146.0   1.00  \n",
       "D_true[10]    1383.0   1.00  \n",
       "D_true[11]    1407.0   1.00  \n",
       "D_true[12]     834.0   1.01  \n",
       "D_true[13]    1239.0   1.00  \n",
       "D_true[14]    1316.0   1.00  \n",
       "D_true[15]    1523.0   1.00  \n",
       "D_true[16]    1201.0   1.00  \n",
       "D_true[17]     834.0   1.01  \n",
       "D_true[18]    1285.0   1.00  \n",
       "D_true[19]    1224.0   1.00  \n",
       "D_true[20]    1011.0   1.00  \n",
       "D_true[21]    1406.0   1.00  \n",
       "D_true[22]    1176.0   1.00  \n",
       "D_true[23]    1147.0   1.00  \n",
       "D_true[24]    1251.0   1.00  \n",
       "D_true[25]    1430.0   1.00  \n",
       "D_true[26]    1306.0   1.00  \n",
       "D_true[27]    1404.0   1.00  \n",
       "D_true[28]    1155.0   1.00  \n",
       "D_true[29]    1183.0   1.01  \n",
       "D_true[30]    1093.0   1.00  \n",
       "D_true[31]    1119.0   1.00  \n",
       "D_true[32]    1252.0   1.00  \n",
       "D_true[33]    1228.0   1.00  \n",
       "D_true[34]    1289.0   1.00  \n",
       "D_true[35]    1343.0   1.01  \n",
       "D_true[36]    1109.0   1.00  \n",
       "D_true[37]     994.0   1.00  \n",
       "D_true[38]    1124.0   1.00  \n",
       "D_true[39]    1281.0   1.00  \n",
       "D_true[40]    1202.0   1.00  \n",
       "D_true[41]    1365.0   1.00  \n",
       "D_true[42]    1244.0   1.00  \n",
       "D_true[43]    1535.0   1.00  \n",
       "D_true[44]    1066.0   1.00  \n",
       "D_true[45]    1113.0   1.00  \n",
       "D_true[46]    1081.0   1.00  \n",
       "D_true[47]    1277.0   1.00  \n",
       "D_true[48]    1006.0   1.01  \n",
       "D_true[49]    1121.0   1.01  \n",
       "a             1639.0   1.00  \n",
       "bA            1604.0   1.00  \n",
       "bM            1620.0   1.00  \n",
       "sigma          413.0   1.01  "
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idata_1, hdi_prob=0.89)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Code 15.5\n",
    "\n",
    "What happens when there is a measurement error on predictor variables as well ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "dlist = dict(\n",
    "  D_obs = d.Divorce.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(jnp.float32),\n",
    "\tD_sd = (d[\"Divorce SE\"].values / d.Divorce.std()).astype(jnp.float32),\n",
    "\tM_obs = d.Marriage.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(jnp.float32),\n",
    "\tM_sd = (d[\"Marriage SE\"].values / d.Marriage.std()).astype(jnp.float32),\n",
    "\tA = d.MedianAgeMarriage.pipe(lambda x: (x - x.mean()) / x.std()).values.astype(jnp.float32),\n",
    "\tN = d.shape[0]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_15_2(A, M_sd, D_sd, N):\n",
    "\tdef _generator():\n",
    "\t\ta = yield tfd.Normal(loc=0.0, scale=0.2, name=\"a\")\n",
    "\t\tbA = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bA\")\n",
    "\t\tbM = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bM\")\n",
    "\t\tsigma = yield tfd.Exponential(rate=1.0, name=\"sigma\")\n",
    "\t\tM_true = yield tfd.Sample(tfd.Normal(loc=0.0, scale=1.0), sample_shape=N, name=\"M_true\")\n",
    "\t\tyield tfd.Normal(loc=M_true, scale=M_sd, name=\"M_obs\")\n",
    "\t\tmu = a + bA * A + bM * M_true\n",
    "\t\tD_true = yield tfd.Normal(loc=mu, scale=sigma, name='D_true')\n",
    "\t\tyield tfd.Normal(loc=D_true, scale=D_sd, name='D_obs')\n",
    "\n",
    "\treturn tfd.JointDistributionCoroutineAutoBatched(_generator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "jd_15_2 = model_15_2(dlist[\"A\"], dlist[\"M_sd\"], dlist[\"D_sd\"], dlist['N'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "jd_15_2_pin = jd_15_2.experimental_pin(M_obs=dlist[\"M_obs\"], D_obs=dlist[\"D_obs\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jax.jit\n",
    "def run_chain(init_state, seed, step_size=1e-3, num_steps=500, burnin=1000):\n",
    "\n",
    "  def target_log_prob(*states):\n",
    "    return jd_15_2_pin.unnormalized_log_prob(states)\n",
    "  \n",
    "  bijectors = [tfb.Identity(), tfb.Identity(), tfb.Identity(),\n",
    "               tfb.Exp(), tfb.Identity(), tfb.Identity(),]\n",
    "\n",
    "  kernel = tfp.mcmc.NoUTurnSampler(\n",
    "      target_log_prob_fn=target_log_prob,\n",
    "      step_size=step_size)\n",
    "  \n",
    "  kernel = tfp.mcmc.TransformedTransitionKernel(\n",
    "    inner_kernel=kernel,\n",
    "    bijector=bijectors)\n",
    "\n",
    "  kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(\n",
    "    inner_kernel=kernel,\n",
    "    num_adaptation_steps=burnin)\n",
    "\n",
    "  chain_state = tfp.mcmc.sample_chain(\n",
    "      num_results=num_steps,\n",
    "      num_burnin_steps=burnin,\n",
    "      current_state=init_state,\n",
    "      kernel=kernel,\n",
    "      seed=seed,\n",
    "      trace_fn=None)\n",
    "  \n",
    "  return chain_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_chains = 4\n",
    "init_state = jd_15_2_pin.sample_unpinned(num_chains, seed=key())\n",
    "post_2 = run_chain(init_state=init_state, seed=key())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "post_2_dict = {k: jnp.moveaxis(v, 0, 1) for k, v in post_2._asdict().items()}\n",
    "dict_2 = dict(posterior=post_2_dict)\n",
    "idata_2 = az.from_dict(**dict_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_5.5%</th>\n",
       "      <th>hdi_94.5%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>D_true[0]</th>\n",
       "      <td>1.126</td>\n",
       "      <td>0.350</td>\n",
       "      <td>0.577</td>\n",
       "      <td>1.663</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2396.0</td>\n",
       "      <td>1264.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[1]</th>\n",
       "      <td>0.765</td>\n",
       "      <td>0.564</td>\n",
       "      <td>-0.093</td>\n",
       "      <td>1.697</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.011</td>\n",
       "      <td>1765.0</td>\n",
       "      <td>1350.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[2]</th>\n",
       "      <td>0.423</td>\n",
       "      <td>0.336</td>\n",
       "      <td>-0.111</td>\n",
       "      <td>0.923</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.006</td>\n",
       "      <td>3288.0</td>\n",
       "      <td>1540.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[3]</th>\n",
       "      <td>1.460</td>\n",
       "      <td>0.456</td>\n",
       "      <td>0.761</td>\n",
       "      <td>2.241</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>2007.0</td>\n",
       "      <td>1264.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D_true[4]</th>\n",
       "      <td>-0.895</td>\n",
       "      <td>0.132</td>\n",
       "      <td>-1.102</td>\n",
       "      <td>-0.687</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2031.0</td>\n",
       "      <td>1189.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M_true[49]</th>\n",
       "      <td>1.273</td>\n",
       "      <td>0.705</td>\n",
       "      <td>0.249</td>\n",
       "      <td>2.455</td>\n",
       "      <td>0.022</td>\n",
       "      <td>0.015</td>\n",
       "      <td>1046.0</td>\n",
       "      <td>1187.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>-0.044</td>\n",
       "      <td>0.099</td>\n",
       "      <td>-0.193</td>\n",
       "      <td>0.123</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1911.0</td>\n",
       "      <td>1675.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bA</th>\n",
       "      <td>-0.546</td>\n",
       "      <td>0.159</td>\n",
       "      <td>-0.807</td>\n",
       "      <td>-0.305</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1162.0</td>\n",
       "      <td>1565.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bM</th>\n",
       "      <td>0.201</td>\n",
       "      <td>0.209</td>\n",
       "      <td>-0.134</td>\n",
       "      <td>0.529</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>818.0</td>\n",
       "      <td>1264.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.558</td>\n",
       "      <td>0.108</td>\n",
       "      <td>0.385</td>\n",
       "      <td>0.723</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>605.0</td>\n",
       "      <td>746.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>104 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             mean     sd  hdi_5.5%  hdi_94.5%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "D_true[0]   1.126  0.350     0.577      1.663      0.007    0.005    2396.0   \n",
       "D_true[1]   0.765  0.564    -0.093      1.697      0.013    0.011    1765.0   \n",
       "D_true[2]   0.423  0.336    -0.111      0.923      0.006    0.006    3288.0   \n",
       "D_true[3]   1.460  0.456     0.761      2.241      0.010    0.008    2007.0   \n",
       "D_true[4]  -0.895  0.132    -1.102     -0.687      0.003    0.002    2031.0   \n",
       "...           ...    ...       ...        ...        ...      ...       ...   \n",
       "M_true[49]  1.273  0.705     0.249      2.455      0.022    0.015    1046.0   \n",
       "a          -0.044  0.099    -0.193      0.123      0.002    0.002    1911.0   \n",
       "bA         -0.546  0.159    -0.807     -0.305      0.005    0.003    1162.0   \n",
       "bM          0.201  0.209    -0.134      0.529      0.007    0.005     818.0   \n",
       "sigma       0.558  0.108     0.385      0.723      0.004    0.003     605.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "D_true[0]     1264.0    1.0  \n",
       "D_true[1]     1350.0    1.0  \n",
       "D_true[2]     1540.0    1.0  \n",
       "D_true[3]     1264.0    1.0  \n",
       "D_true[4]     1189.0    1.0  \n",
       "...              ...    ...  \n",
       "M_true[49]    1187.0    1.0  \n",
       "a             1675.0    1.0  \n",
       "bA            1565.0    1.0  \n",
       "bM            1264.0    1.0  \n",
       "sigma          746.0    1.0  \n",
       "\n",
       "[104 rows x 9 columns]"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idata_2, hdi_prob=0.89)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Code 15.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 500, 50)"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "post_2_dict[\"D_true\"].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAHACAYAAAB5+ch6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACzxklEQVR4nOzde1xM6R8H8M/MdNPVpYtSmQiVogkrrSj3S5Fkkc1l2WVZ1srPyu66LWJXbj+77C6WRbEq2Y11L6LBWkUp91LoolCj6DLz/P44v5mtFDXNzJmp5/169Wo758w537FT853nfJ/vwyGEEFAURVEURakJLtsBUBRFURRFVUWTE4qiKIqi1ApNTiiKoiiKUis0OaEoiqIoSq3Q5ISiKIqiKLVCkxOKoiiKotQKTU4oiqIoilIrNDmhKIqiKEqtaLEdgLqTSCR48uQJjIyMwOFw2A6HoiiKojQGIQQikQhWVlbgcus/HkKTk3d48uQJbGxs2A6DoiiKojRWdnY2rK2t6308TU7ewcjICADzD2tsbMxyNBRFURSlOYqLi2FjYyN7L60vmpy8g/RWjrGxMU1OKIqiKEoODS2LoAWxFEVRFEWpFZqcUBRFURSlVmhyQlEURVGUWqHJCUVRFEVRaoUmJxRFURRFqRWanFAURVEUpVZockJRFEVRlFqhyQlFURRFUWqFJicURVEURakV2iGWoiiKUhqJBMjKAkQiwMgIsLUFGrD+G9VM0eSEoiiKUor0dODwYeDWLeD1a0BPD3BwAMaMARwd2Y6OUmc0OaEoiqIULj0d2LIFKCgAbGwAAwOgpARISgKys4F582iCQtWNDq5RFEVRCiWRMCMmBQWAkxNgbAzweMx3Jydme0wMcxxF1YYmJxRFUZRCZWUxt3JsbAAOB3jy5Cr+/nsbxOIKcDiAtTUzspKVxXaklLqiyQlFURSlUCIRU2NiYMD8nJCwBpcvbwaXywPAbH/9mjmOompDkxOKoihKoYyMmOLXkhLg2bN7uHUrBn36LACHw7zllJQw+42MWA6UUls0OaEoiqIUytaWmZWTnQ1curQJ+vqm6NYtCABACPDoEVMMa2vLcqCU2qLJCUVRFKVQXC4zXdjQ8BmuXfsVLi6zweG0QFERkJYGmJoCfn603wlVN415aYSGhqJXr14wMjKCubk5/Pz8cPv27Xc+7ty5c+jRowf09PTQoUMHbN++XQXRUhRFNW+OjoCJyXZwOGLY2MzGnTtAYSHg5kanEVPvpjF9Ts6dO4c5c+agV69eqKysxFdffYUhQ4YgLS0NBtKqqxoyMjIwYsQIfPzxx9i3bx8uXryI2bNnw8zMDGPHjlXxM6Aoimo+ysrKcPDgfzF16mR89ZU57RBLNQiHEELYDkIeT58+hbm5Oc6dO4d+/frVesyXX36JP/74A+np6bJts2bNwvXr1yEUCut1neLiYpiYmKCoqAjGxsYKiZ2iKKqp2717N6ZNm4b09HQ4ODiwHQ7FEnnfQzU2fy0qKgIAtG7dus5jhEIhhgwZUm3b0KFDcfXqVVRUVNT6mLKyMhQXF1f7oiiKouqPEIINGzZg5MiRNDGh5KKRyQkhBAsWLEDfvn3h7Oxc53G5ubmwsLCots3CwgKVlZUoKCio9TGhoaEwMTGRfdnY2Cg0doqiqKbu1KlTSElJQXBwMNuhUBpKI5OTzz77DDdu3EBERMQ7j+VwONV+lt7FqrldKiQkBEVFRbKv7OzsxgdMURTVjISFhUEgEMDLy4vtUCgNpTEFsVJz587FH3/8gfPnz8Pa2vqtx7Zt2xa5ubnVtuXn50NLSwtt2rSp9TG6urrQ1dVVWLwURVHNSUpKCk6ePIl9+/bV+SGQot5FY0ZOCCH47LPPEB0djbNnz8LOzu6dj+nTpw9OnTpVbdvJkyfRs2dPaGtrKytUiqKoZmvDhg2wtrbGBx98wHYolAbTmORkzpw52LdvH8LDw2FkZITc3Fzk5ubi1atXsmNCQkIwefJk2c+zZs3Cw4cPsWDBAqSnp2PXrl3YuXMnFi5cyMZToCiKatJycnKwf/9+zJs3j34ApBpFY5KTbdu2oaioCF5eXrC0tJR9HTx4UHZMTk4Osqosc2lnZ4djx44hPj4erq6u+Pbbb7Flyxba44SiKEoJtm7dCl1dXXz88cdsh0JpOI3tc6IqtM8JRVHUu5WUlMDGxgaTJ0/Gpk2b2A6HUhPNrs8JRVEUpT52796NoqIifP7552yHQjUBNDmhKIqiGkUsFmPjxo0YO3ZsvSYrUNS7aNxUYoqiKEq9/PHHH7h//z7279/PdihUE0FHTiiKoqhG2bBhA95//3307t2b7VCoJoKOnKiYRAJkZYGu0ElRVJNw5coVXLhwAdHR0WyHQjUhNDlRofR04PBh4NYt4PVrQE8PcHAAxowBHB3Zjo6iKKrhwsLC0LFjR4waNYrtUKgmhCYnKpKeDmzZAhQUADY2gIEBUFICJCUB2dnAvHk0QaEoSrNkZmYiMjISW7ZsAY/HYzscqgmhNxRUQCJhRkwKCgAnJ8DYGODxmO9OTsz2mBjmOIqiKE2xefNmmJiYYOrUqWyHQjUxNDlRgaws5laOjQ3A4QCXLm1CePhIAMzP1tbMyEqV5rYURVFq7cWLF9ixYwc+/fRTGBgYsB0O1cTQ5EQFRCKmxkT6+8vh8PDgwWmIxRUAmO2vXzPHURRFaYJffvkF5eXl+Oyzz9gOhWqCaHKiAkZGTPFrSQnzc9u2rhCLy1FQkA6A2a6nxxxHURSl7ioqKrBlyxYEBgbC0tKS7XCoJogmJypga8vMysnOBggB2rbtDgDIyUkCIcCjR0wxrK0ty4FSFEXVw++//45Hjx5hwYIFbIdCNVE0OVEBLpeZLmxqCqSlAa9fG6NVK3s8fJiEtDRmu58f7XdCUZT6I4QgLCwMQ4YMgYuLC9vhUE0UnUqsIo6OzHRhaZ8TAwMBcnKS8cEHTGJCpxFTVNPU1BovxsfHIykpCSdOnGA7FKoJ4xBCCNtBqDN5l3uui/QP1caNofj113V4/vw5eDyOAiKlKErdNMXGiz4+PsjKysL169fB4dC/XdTbyfseqsH5u2bicgE+Hxg2zBUiUREePsxgOySKopRA2ngxKYm5ddulC/M9KYnZnp7OdoQNl56ejqNHj2LBggU0MaGUiiYnLBEIBACApKQkliOhKErRmmrjxY0bN6Jt27aYOHEi26FQTRxNTljStm1btG3bliYnFNUE1Wy8+OzZfaSmHkBlZZnGNl7Mz8/Hb7/9hrlz50JXV5ftcKgmjiYnLBIIBEhOTmY7DKqZk0iAzEwgJYX5rmmf5tVRzcaLd+7EIiZmKjgc5k+uJjZe/PHHH8Hj8TBr1iy2Q6GaATpbh0UCgQC7d+9mOwyqGWuKBZvqoGrjRWNjIDc3CRYW3cDjaQPQvMaLr169wg8//IBp06ahdevWbIdDNQN05IRFrq6uePLkCfLz89kOhWqGmmLBprqo2XgxJ+caLC3dAEAjGy/u3bsXhYWFmD9/PtuhUM0ETU5YRItiKbY01YJNdVG18WJKyis8fZoGc3M3FBVB4xovSiQSbNiwAX5+frC3t2c7HKqZ0IBfjaarQ4cOMDIyonUnlMrVLNg8cWIBbt48BICulK0o0saLFhYpIESM8nIBCgsBNzdmu6bcNjt27Bhu376N4OBgtkOhmhFac8IiLpcLV1dXOnJCqVzNgs0HD05BLC5H167jADDbHz/WnIJNde3C6ugIuLhcA4/Hw9q1LjA1VZ/Y6issLAy9e/eGh4cH26FQzQhNTlgmEAhw/PhxtsOgmpmaBZsmJrYoLs6W7dekgk11L+q9fj0JXbt2Rc+eemyH0mDXrl1DfHw8fv/9d9p0jVIpDcrfmyZXV1fcvXsXL1++ZDsUqhmpWbBpbGyDoiImOdGkgk1NKOq9du0a3Nzc2A5DLhs2bACfz8eYMWPYDoVqZmhywjKBQABCCK5fv852KFQzUnOlbF1dGxQXZ2tUwaYmFPVWVFTgxo0bGpmcPHr0CAcPHsT8+fOhpUUH2SnVUuM/Pc2Dk5MTtLW1aVEspXLSgk1m0pgNSksLkJf3SmMKNmsW9UZFBeLIkekA1KeoNy0tDeXl5bKZeZpky5YtMDAwwEcffcR2KFQzRNNhluno6MDZ2ZkWxVKscHRkboXY2dkgMRH4+ONH6N+/k1qPmEjVLOp9+TIHhoaWsv3qUNR77do1cDgcdO/enb0g5CASifDzzz9j5syZMKqj8Ehdi5CppoEmJ2qAztih2MTlAr162QAACMkGl9uJ5Yjqp2ZRb0VFKbS1DWT71aGoNykpCZ07d67zDV5d7dy5EyUlJZg7d26t+9W9CJnSfDTPVQMCgQCpqamoqKhgOxSqmbK2tgYAZGdnv+NI9VGzqLe8vAQ6Okxyoi5FvZpYDFtZWYlNmzZh/PjxstdFVZpQhExpPpqcqAGBQIDy8nKkpaWxHQrVTOnp6cHMzEyjkpOaRb1lZSXg8QzUpqhXLBYjOTlZ4+pNoqOj8fDhw1qbrtUsQjY0FCM392+1KkKmmgaNSk7Onz8PX19fWFlZgcPhICYm5q3Hx8fHg8PhvPF169Yt1QRcT927dweHw6FFsRSrbGxsNCo5AaoX9ZaVlaC4WF9turDevXsXJSUlGjVyQghBWFgYvL29a02qahYhJyaux86dffD8eYbaFCFTTYNG1ZyUlJSge/fumDZtGsaOHVvvx92+fRvGxsayn83MzJQRntyMjIxgb2+PpKQkTJkyhe1wqGZKE5MT4N+i3lWrSjB8uAFmz1aP4kxpHZkmjZxcvHgRV65cQWxsbK37qxYhP3nyD+LivoaHx0K0amUHQD2KkKmmQaOSk+HDh2P48OENfpy5uTlatmyp+IAUiBbFUmyzsbFBfHw822HISYJXr0rRqZMB+Hy2Y2Fcu3YNfD4frVu3ZjuUegsLC4ODg0Odf2elRcgvXpQgOjoQFhbd4e29UrZfHYqQqaZBo27ryEsgEMDS0hIDBw5EXFzcW48tKytDcXFxtS9VxZicnAwJvVlLsURTR04A4NWrVwAAAwODdxypOppSDCsWixEfH48NGzbgyJEjmD9/Prh1DDtJi5CPHv0CxcWP4O+/HzyeDgD1KUKmmoYmnZxYWlri559/RlRUFKKjo9GlSxcMHDgQ58+fr/MxoaGhMDExkX3Z2NioJFaBQIDi4mJkZGSo5HoUVZONjQ2Kioog0sAx+dLSUgDqk5wQQnDt2jW1v6UTHR0Ne3t7eHt7Izg4GIQQrFmzBtHR0bUez+UCRkaHkZHxC1xcNkFbuwsqK6E2RchU09GkX0JdunTBxx9/DDc3N/Tp0wc//vgjRo4cifXr19f5mJCQEBQVFcm+VPVJUvpHjBbFUmyRJuKaOHpSUlICQH2Sk8zMTLx48UKtR06io6MREBAAFxcXHD9+HC1atMCMGTPQvXt3BAQE1JqgPHnyBMuXz8DAgX7w95+BwkLgzh2oTREy1XRoVM2JIri7u2Pfvn117tfV1YWurq4KI2JYWFjA0tISSUlJDSr2pShFqZqcODk5sRxNw6hbciKtH1PX5EQsFiM4OBg+Pj6IiYlBaGiobNSkTZs28PPzw8KFCzF69GjweDwAgEQiwZQpU6Crq4sDB35B69Yc2iGWUppml5wkJSXB0tLy3QeygBbFUmySTtHX5JETfX19liNhXLt2DZaWlmjbti3bodQqISEBmZmZiIiIQEVFBbZu3YrJkyfLZjKGhITAw8MDCQkJ8PLyAgBs3LgRp0+fxqlTp2BqagoAalN8TDU9GpWcvHz5Evfu3ZP9nJGRgeTkZLRu3Rq2trYICQnB48eP8dtvvwEANm3aBD6fj65du6K8vBz79u1DVFQUoqKi2HoKbyUQCPDrr7+yHQbVTGlra8PS0lKjkxN1GTlR92LYnJwcAICzszMqKiowc+ZMTJw4Ubbf2dm52nHJyckICQlBcHAwBg0apPqAqWZHowbhrl69CoFAIKvPWLBgAQQCAZYuXQqA+UXKqtL9p7y8HAsXLkS3bt3g6emJCxcu4OjRo/D392cl/ncRCATIyclBXl4e26FQzZSmzthRp+SEEIJ//vlHrYthpaPHqampMDQ0xPLly9GlSxfZ/tTUVNlxpaWlCAwMRNeuXbF69WpW4qWaH40aOfHy8gIhpM79u3fvrvbzokWLsGjRIiVHpThVi2KHDh3KcjRUc0STk8bLyclBfn6+Wo+ceHp6gs/nY82aNYiJiak2dVgikSA0NBR2dnbw9PTEvHnzkJGRgWvXrrFSj0c1Txo1ctLU2dnZwdjYmNadUKzR1OREnaYSq3sxLADweDyEhYUhNjYWfn5+EAqFEIlEEAqF8PPzQ2xsLNavX49jx47hxx9/xIYNG+BIp+FQKqRRIydNHZfLRffu3WlyQrFGmpwQQsDhcNgOp95KSkqgq6srm1nCpmvXrsnq4NSZv78/IiMjERwcDA8PD9l2Ozs7REZGwsPDAy4uLvDx8cGsWbNYjJRqjmhyomYEAgH++usvtsOgmikbGxuUlpbi+fPnGtV2vaSkRC1GTYB/i2E1Ibnz9/fH6NGjkZCQgJycHFhaWsLT0xMcDgcjRowAj8fDzp07NeK5UE0LTU7UjEAgwH//+1+IRCIY0QUqKBWr2utE05ITdZpGPH78eLbDqDcejyebLiy1efNmnDhxAn/99RfMzc3ZCYxq1mjNiZoRCAQghODGjRtsh0I1Q5raJVZdRk4KCwuRlZWl1vUm73Ljxg0sWrQIn3/+OYYNG8Z2OFQzRZMTNePo6AgdHR1ad0KxwsLCAlpaWjQ5kZMmFMO+zatXrzBp0iR07twZa9euZTscqhmjt3XUjI6ODrp27UqTE4oVPB4P7dq1o8mJnK5duwZDQ0PY29uzHYpcFi9ejLt37+Lq1avQ09NjOxyqGaPJiRoSCAQ0OaFYo4nTiUtLS1lNTsRiMRISEnD48GHY2dm9tR+Tuvrrr7+wZcsWbNmyRdYhlqLYQm/rqCGBQICbN2+ivLyc7VCoZkgTkxM2R06io6Nhb28Pb29vXLp0CSkpKbC3t691VV91lZ+fj6lTp2LYsGH47LPP2A6Homhyoo4EAgHKy8uRnp7OdihUM0STk/qLjo5GQEAAXFxccPr0aQDA119/DRcXFwQEBGhEgkIIwUcffQRCCHbv3k2nDVNqgSYnaqhbt27gcDj01g7FChsbGzx69AgSiYTtUOqNjeRELBYjODgYPj4+iImJgY6ODgDggw8+QExMDHx8fLBw4UKIxWKVxtVQP/74I44ePYpff/0VFhYWbIdDUQBocqKWjIyMYG9vT5MTihU2NjYoLy9Hfn4+26HUi1gsRkFBAbKyshAfH6+yZCAhIQGZmZlYsmQJuFwurl27Bj09PTg6OoLL5SIkJAQZGRlISEhQSTzyuHnzJhYuXIg5c+Zg5MiRbIdDUTI0OVFTtCiWYosm9TqR1ntkZ2fj6NGj8Pb2Vlm9R05ODgDIikczMzPRrVs3aGlpVdsuPU7dlJWVITAwEHZ2dvj+++/ZDoeiqqHJiZoSCARITk7WqKF1qmnQlOSkar2HkZERVq5cCaFQqLJ6D0tLSwBAamoqAGDjxo04d+6cbL90u/Q4dbNkyRLcunULERERaNGiBdvhUFQ1NDlRUwKBACKRCBkZGWyHQjUzpqam0NPTU+vkpGa9x+vXr9G6dWu4u7urrN7D09MTfD4fa9askX2IkPYGkUgkCA0NhZ2dHTw9PZUWg7xOnjyJDRs2YO3atejevTvb4VDUG2hyoqZcXV0BgN7aoVSOw+HA2tparZOTqvUeYrEYFRUVsoJYVdV78Hg8hIWFITY2Fn5+fhAKhRCJRBAKhfDz80NsbCzWr1+vFislV1VQUIApU6Zg8ODB+Pzzz9kOh6JqRZMTNWVhYQFLS0uanFCsUPfpxFXrPQghWL58OQQCgWy/quo9/P39ERkZiZSUFHh4eMDY2BgeHh5ITU1FZGQk/P39lXr9hiKEYMaMGaioqMDu3bvB5dK3AEo90Q6xakxad0JRqmZjY4O7d++yHUadqtZ7uLu7Y9myZdX2q7Lew9/fH6NHj0ZCQgJycnJgaWkJT09PtRsxAYCff/4ZR44cQUxMDKysrNgOh6LqRNNmNUZn7FBsUfeRk9rqPaTYqPfg8Xjw8vLCxIkT4eXlpZaJya1bt/DFF19g5syZGD16NNvhUNRb0eREjQkEAuTk5CAvL4/tUKhmxsbGBk+ePEFlZSXbodRKU+s92FJeXo7AwEDY2toiLCyM7XAo6p1ocqLGaFEsxRYbGxtIJBK17dEBaF69B5u++eYbpKamIjw8XC1Wb6aod6HJiRqzs7ODsbExTU4olZPWI+zatUulXVcbyt/fH/fu3UNcXBzCw8MRFxeHu3fv0sSkirNnz+L777/H6tWr4ebmxnY4FFUvHKKJa3urUHFxMUxMTFBUVARjY2OVX79///5o27YtDh48qPJrU81TdHQ05s+fX63mhM/nIywsjL7pa5jCwkJ0794dXbp0walTp+jsHErl5H0Ppa9UNUeLYilVknZd7d69O/T19bFq1SqVdl2lFIcQgpkzZ6K0tBR79uyhiQmlUeirVc0JBALcvXsXIpGI7VCoJq5q19UjR46gffv2yM/PV2nXVUpxdu3ahaioKPzyyy+wtrZmOxyKahCanKg5aVHs9evX2Q2EavJqrrJbdTqxpqyySzHu3LmDefPmYfr06Rg7dizb4VBUg9EmbGrOyckJOjo6SEpKQt++fdkOh2rCaq6y+91330FXV1e2X91X2aUYFRUVmDRpEqysrLBp0ya2w6EoudDkRM1pa2vD2dmZdoqllK5m19WaC8Kp+yq7bJJIgKwsQCQCjIwAW1uArRKP5cuXIzk5GYmJiTA0NGQnCIpqJJqcaACBQIBr166xHQbVxFXtuhoTE1OtgFLdV9llU3o6cPgwcOsW8Po1oKcHODgAY8YAjo6qjeXcuXMIDQ3F6tWr0atXL9VenKIUiNacaABXV1ekpqaivLyc7VCoJox2XW249HRgyxYgKQkwNQW6dGG+JyUx29PTVRfL8+fPERQUBE9PTyxatEh1F6YoJaDJiQYQCASoqKhAWloa26FQTRztulp/EgkzYlJQADg5AcbGAI/HfHdyYrbHxDDHKRshBLNmzYJIJMLevXtpAklpPHpbRwN0794dHA4HSUlJstk7FKUsmrTKLpuysphbOTY2AIcDXLiwDgYGZhAIPgKHA1hbMyMnWVkAn6/cWPbu3Yvff/8dBw8ehK2trXIvRlEqoFEjJ+fPn4evry+srKzA4XAQExPzzsecO3cOPXr0gJ6eHjp06IDt27crP1AFMzQ0RKdOnWhRLKUymrDKLsCMSmRmAikpzHdVjFJIiURMjYl0qZp//tmO/PxU2X4DA2a/slsU3b9/H3PmzMGUKVPwwQcfKPdiFKUiGjVyUlJSgu7du2PatGn1mrufkZGBESNG4OOPP8a+fftw8eJFzJ49G2ZmZho39592iqWo6tguRDUyYq5ZUgJwODl48SITNjYesv0lJcx+IyPlxSCdNmxubo4tW7Yo70IUpWIalZwMHz4cw4cPr/fx27dvh62trWyuv6OjI65evYr169drXHLi6uqKY8eOQSKR0DbUVLMnLUQtKGBuqxgYMMlAUhKQnQ3Mm6f8BMXWlkmGkpIADkcIALC27gMAIAR49Ahwc2OOU5ZVq1bh6tWruHDhAitrf1GUsjTpdzmhUIghQ4ZU2zZ06FBcvXoVFRUVtT6mrKwMxcXF1b7UgUAggEgkwoMHD9gOhaJYpS6FqFwuM0pjagokJyfCyMgW+vrtUFQEpKUx2/38lNfv5MKFC1i1ahWWLVsGd3d35VyEoljSpJOT3NxcWFhYVNtmYWGByspKFBQU1PqY0NBQmJiYyL5sbGxUEeo7CQQCAKB1J1SzV7UQ9dmzu9i40RZPnvwDAG8UoiqboyMzSlNWJkTr1n1w5w5QWMiMmChz9KaoqAgffvgh+vTpg5CQEOVchKJYpFG3deTB4XCq/UwIqXW7VEhICBYsWCD7ubi4WC0SFHNzc1hZWSEpKQkBAQFsh0OxSJ26kbKhaiFqXl4xiouzweH8+w9gYAA8fqz8QlSpDh3KkJt7FYsXj0dAgGr+n8yZMwfPnz9HfHw8tLSa/J9xqhlq0q/qtm3bIjc3t9q2/Px8aGlpoU2bNrU+RldXt9p6IuqEFsVSbBeBqoOqhagVFSUAAG1tfdl+VRSiVpWUlITy8nL4+PSBi4vyr7d//37ZF1/Zc5QpiiVN+vNWnz59cOrUqWrbTp48iZ49e0JbW5ulqOTn6upKk5NmTJ26kbJJWoianQ2UlTHJiY4OM59XWojq6KjcQtSqEhMT0aJFC5X0IMrIyMDs2bMxadIkBAYGKv16FMUWjUpOXr58ieTkZFndRUZGBpKTk5H1/5vLISEhmDx5suz4WbNm4eHDh1iwYAHS09Oxa9cu7Ny5EwsXLmQj/EYTCATIzc19YzSIavpqFoHq6LzAixd3WelGyraqhagZGUxywuEYqKwQtSahUKiSDzyVlZUICgpCq1at8MMPPyj1WhTFNo1KTq5evQqBQCArDl2wYAEEAgGWLl0KgFnKPatKFZydnR2OHTuG+Ph4uLq64ttvv8WWLVs0bhqxFC2Kbb5qdiNNSFiD/fuHAVB9Eag6kBaiWluXAgAePjRQSSFqTYQQJCYmwsPD490HN1JoaCiEQiH2798PExMTpV+PotikUTUnXl5esoLW2uzevfuNbf37928yK/ra2dnBxMQESUlJGDZsGNvhUCpUsxupiUl7FBVlgxAJOByuyotA1YGjI+DlVYKYGC0sW6bDSnFwdnY2njx5gj59+ij1OkKhECtWrMDXX3+N999//53Hi8ViuvwApdE0Kjlp7jgcDq07aaaqFoEaGwMmJraQSCrw8mUujIysVF4Eqi5evSqBgYGBSgpRa5OYmAgASk1OiouLMWnSJPTq1QvffPPNO4+Pjo5GcHAwMjMzZdv4fD7CwsKa1cKNNEHTbBp1W4eiRbHNVdUiUEKAli3bAwBevHjIShGouigpKYG+vv67D1QSoVCIjh07wtzcXGnXmDdvHp4+fYp9+/a9c9pwdHQ0AgIC4OLiAqFQCJFIBKFQCBcXFwQEBCA6OlppcaqT6Oho2Nvbw9vbG4GBgfD29oa9vX2zef5NAU1ONIxAIMC9e/fUpnMtpRpVi0DT0gAOh0lOcnKyWCkCVRclJczICVuUXW9y8OBB7NmzB1u3bkXHjh3feqxYLEZwcDB8fHwQExMDd3d3GBoawt3dHTExMfDx8cHChQshFouVFq86oAla09DM/pRpPmlR7I0bN1iOhFI1aRGoQAC8fGkCLS1j5OU9VHkRqDopLS1lLTkpLS1FcnKy0m7pZGVlYebMmRg/fny1WYh1SUhIQGZmJpYsWQIOh4MdO3bgyZMnAAAul4uQkBBkZGQgISFBKfGqg5oJmrGxMQwMDJpdgtYU0OREwzg6OkJXV5fe2mmmHB2BxYuBlSsBO7v26No1C19+2TwTE4DdkZOrV6+isrJSKSMnYrEYQUFBMDY2xrZt2+rsaF1VTk4OAMDZ2Rn379/H7Nmz0a5dO7z//vvYuHGjbIaP9LimqGqC9ujRI7i5uWHRokUghDSbBK2poMmJhtHW1oazszNNTpoxLhfg84HOnW3x/PnDZncrpyo2kxOhUAhDQ0M4Ozsr/NzfffcdEhISsHfvXrRq1apej7G0tAQApKamwt7eHnl5edizZw9MTU0REhKCrl27AgDi4+Nx//59hcesDqomaLa2tli3bh3Wr18vazch/X/VlBO0pqIZ/1nTXLQolgKA9u3bV+vr0xyxWRCbmJiI3r17K3wGyN9//42lS5ciJCQE/fv3r/fjPD09wefzsWbNGkgkErRq1QqTJ0/GkSNHkJubCzc3N+jr62Pv3r2wt7eHm5sb1qxZg9u3bys0fjZVTdAA4PPPP8e6deuwatUqrFq1SrZdehylvmhyooEEAgFu3ryJ8vJytkOhWGRra4uHDx+yHQar2Bo5IYRAKBQqvN7k5cuXCAwMhKurK5YvX96gx/J4PISFhSE2NhZ+fn7VikEnT56MpKQk7N27F0+fPsWhQ4fQuXNnrFmzBg4ODnBxccGKFStw8+bNt/aSUnc1EzQAWLRoEVauXIlvvvkGM2bMgJ2dHTw9PVmOlHoXmpxoIIFAgIqKCqSlpbEdCqVCEgmQmQmkpDDfbWzao6ioCEVFRWyHxgqxWIy8vDzk5OQgPj5epUWO9+/fx9OnTxVebzJ//nw8efIE4eHhcrXD9/f3R2RkJFJSUuDh4QFjY2N4eHggNTUVkZGR8Pf3h4GBAQICAnDgwAE8ffoUMTExcHV1xYYNG+Ds7AxHR0d8/fXXSE5O1rhEpa4EbeDAgbC3t0daWhoGDBhA+51oAkK9VVFREQFAioqK2A5FRiQSEQ6HQ3bt2sV2KJSKpKURsno1IUFBhIwbx3z/5JOLBAC5ceMG2+GpXFRUFOHz+QSA7IvP55OoqCiVXP+3334jAMizZ88adZ7KykoSFxdHwsPDyfLlywkAsmPHjkbHV/W8cXFxpLKy8p2Pef36NYmNjSVTp04lrVq1IgBIx44dyaJFi8iVK1eIRCJpdFyqUtfrY9SoUQQA+emnn9gOsdmQ9z2UdojVQIaGhujUqROSkpIwbdo0tsOhlEy6GnFBAbO2joEB0yn2zh2m18nFi1lwYatFKgukfSx8fHzw+vVrBAUFwd/fH2vWrEFAQIBshECZEhMT4ejoWO9i1drU1slVX18fLVu2bHR8PB4PXl5eDXqMrq4uRo4ciZEjR6KiogJxcXGIjIzErl278N1338HW1hYBAQEYO3Ys3N3dwVXjSmx/f3+MHj36jQ6xXC4Xc+fOxaxZs6Crq4spU6awHSpVB/V9dVFvJRAIaFFsM1BzNWIjIwIOpxLGxoBA0BYcjhb++uths1iNGHizj0VlZSVatmyp8j4Wja03qdoo7OLFi+jXrx/MzMzg6emJcePGsd4oTFtbG0OGDMHPP/+MnJwcnD17Fj4+PggPD8f7778PGxsbzJs3D+fPn1fbniHSBG3ixInw8vICj8cDh8PBli1bMH36dHz00UeIiIhgO0yqDjQ50VACgQDXr1+XFX1RTVPV1YjF4jJ8/70ZbtzYD4D542tkZIP797NqXY24Zo1KY18qij6fPKr2seByudUKYlXVx0IkEslqOuRRM8G6ePEiEhIScODAARw7dkztGoVpaWnB29sbP/zwAx49eoTz58/LEqj+/fujXbt2+PTTT3HmzBlUVlayHe47cblc/PTTT5g0aRKCgoJYTwSp2tHbOhpKIBBAJBLhwYMHsLe3ZzscSkmqrkbM4+miRYvWyM1NBsAMR7dsaYP8/L9x8GAEevf+d3Gz9HRmxOXWLebxenrM2jxjxsjXsE3R55NX1T4WhBDY2dmhbdu2sv2q6GNx5coVSCQSuUdOpAlWREQEbt26ha+++gr/+c9/MGDAAABASEgIPDw8kJCQ0OBbM8rG4/Hg6ekJT09PbNiwAVeuXEFkZCQiIyOxfft2tGnTBmPGjMHYsWMxYMAA6OjosB1yrbhcLnbt2oXy8nJMmDAB0dHR8PHxYTssqgqanGgoV1dXAEBSUhJNTpqwmqsRt23ritxc5nZeeno0cnL+QUVFCRYvPguAWX32iy/CkJ7u/0aNSlISs3BgQ1vd11XzIu/5GqNqHwt3d3fcvHmz2n5V9LFITExEy5Yt4eDgINfjqyZYLVq0wLZt2xAUFCTbrymNwrhcLtzd3eHu7o7vv/8e165dkyUqO3bsQMuWLTF69GgEBARg8ODB0NXVZTvkarS0tLB3716Ul5dj7Nix+PPPPzFkyBC2w6L+j97W0VDm5uawsrKidSdNXM3ViJnkJBlpaVH4/fcA6OpawcDADEVFTD8LZ2cXzJ8fgOvXo+HkxCQ0PB7z3cmJSTBiYup/S6ZmzUtjz9dYtfWx+DdWCUJDQ5XWx0IsFiM+Ph5RUVHo3Lmz3NNsqyZYPB4P06dPrzbCoImNwjgcDnr06IHQ0FDcuXMHycnJmDt3Li5fvgxfX1+YmZlh0qRJOHz4MF69esV2uDLa2to4cOAABg8ejNGjRyMuLo7tkKj/kys5KSsrk7VW/umnnxAdHY2MjAxFx0a9Ay2KbfpqrkZsZOSKsrIi/PXX57C09IGb20KUlhagRQtduLu7Y/PmGLRr54ObNxeCEDEKC+8gP595s+NwAGtrZiSkvo1lq9a8ECJGdnYiSkry5T5fY72t0Zifnx9iY2Oxfv16hfexiI6Ohr29Pby9vXH9+nVcuXIF9vb2ctUrsJlgqQKHw0H37t2xcuVKpKen4+bNm1i4cCFSU1Ph7+8PMzMzfPDBB/j999/x8uVLtsOFjo4OIiMj0bdvX/j4+ODChQtsh0ShgclJYmIiJk6ciJYtW8LLywvz58/Ht99+iw8//BD29vbo1KkTvv/+e4hEImXFS1UhEAiQnJzMdhiUklVdjVhLi1mV+uXLxxg9egmCgvgghODRo0cAgJISLjp3DkFxcQayshJw8mQwTp9eLDuXgQFTM1JUVL/i1qo1L2JxOXbteh/3759843yq/JWvq9FYcnIqwsIi4ebmr9CRnKoza6SzOzZv3gwXFxcEBAQ0OEFhK8Fii5OTE5YuXYrr16/j9u3b+Oqrr3D//n2MHz8eZmZm8Pf3x/79+1ltJqinp4cjR46gV69eGDFiBK5cucJaLBSDQ+o5Njl69Gj8/fffCAwMxKhRo9CzZ89qa1o8ePAACQkJiIiIwPXr1/Hbb79h8ODBSgtcVYqLi2FiYoKioiIYGxuzHU410dHRGDt2LHJycqoVBVJNk0QCPHxI0K1bS7x8WYyiIhFyc5/Ay8sLsbGxcHNzQ2YmEBIiwoEDxvD3D8fDh+fx+PElzJzJjLAVFQEPHgCdOwN5ee8ubs3MBJYuZUZuDA3F+PZbLYwatRMCwUey8xUWMqsk8/kq/eeAWCxGQkICrl7Nwa1blnj92hPl5TyFFuuKxWLY29vDxcUFMTEx+PXXX/HJJ5/g+fPnMDQ0hJ+fH1JTU3H37t0GJxO19Tmxs7PD+vXrld6nRR1kZGQgKioKkZGRuHz5MnR0dDBkyBAEBARg1KhRjeohI6+XL19i6NChSEtLw5kzZ+Dm5qbyGJoaed9D6z1yMmTIEGRmZmL9+vXo16/fG4ttdejQAVOmTMHx48dx+vTp+kdOyU0gYD5F01s7zQOXC9jZceDg0BkAkJaWis6dO+PJkyeyP6K2toCJCXMbx9DQEsbG7VBczIyqEMLcGsrNBTIymISjSxfme1ISU/Sanl79mlVrXjgcHjgcHiory2Tne/SISQBsbVX0j1AFj8eDhYUX7t+fiKIiL5ib8975fBqq5tRlpq7HGcbGxo2euuzv74979+4hLi4O4eHhiIuLw927d5tFYgIwidjChQtx6dIlPHz4EOvWrcOLFy8wbdo0mJubY9iwYdixYwcKCgpUFpOhoSGOHTsGe3t7DBkyBCkpKSq7NlVdvZOTOXPm1HtaWNeuXZvEqIm64/P5MDExoclJMyNtKFVbzQIgwZ07oTAyskNJiSe0tKxRWlqAwsLXuHkTePECaNkS6Nq1fsWtNWteeDxdVFSUo6iI+dnUFPDzY45TNVUU61adWQMA2dnZ1fqbNHZmTW2NwpojW1tbzJ8/HwkJCXj06BE2bdqEsrIyzJw5E23btsWgQYOwbds25ObmKj0WExMTnDhxAjY2Nhg0aBBu3bql9GtSb6KzdTQYh8OBq6srTU6aGTc3N4jFYvz555+11izEx8di1ar1cHPjQSJpBwB49OgJOnQA2rZlEhMOB3j6NA1xcctQUfHqrcWtVWteuFwd5OaWobAQcHNT7TTimqoW63I4QHr6Ydy6FQNAccW6VWfWAMCJEyewZcsW2X5NnFmj7qysrDBnzhzExcUhJycHP/74o6ztvJWVFfr374///ve/ePz4sdJiaN26NU6ePAkzMzMMGDAA9+7dU9q1qNrVOzlp1aoVWrduXa8vSnVoUWzzI72dt3z58jpXn503zx+LFwNffskkJxMnPsakSYCuLlPECgBPn6bj/PmVqKxkpnbq6zP1I//882aRrKMjsHgxYGysi/79y7FyJfDll+wlJsC/xbr6+gQXL36P338fK0tOAMUU69Y2s0a6WnBTmFmj7szNzfHJJ5/g5MmTyMvLw86dO2FoaIjg4GBYW1vDw8MDGzZswMOHDxV+bTMzM5w+fRpGRkYYMGBAtdogSvnq3YRt06ZNsv8uLCzEqlWrMHToUFmXRKFQiBMnTuCbb75ReJBU3QQCATZt2oTi4mK1K9illKNTp05o0aIFjI2Nce/evTcWN5PeGuBygd69rQEAYvEjmJhUb+jG4XAAAIQQPH0KJCcDOTnArl3AkSNvFpVyuUCLFjpo2bJM5cWvtTEyAnR0KnDkyGykpu5A375LMGDAt7L9JSXM8zUykv8a0pk1AQEB8PPzQ0hICJydnZGamorQ0FDExsYiMjKy2d6OUaU2bdpg2rRpmDZtGl68eIE///wTkZGRWLJkCYKDg9GrVy+MHTsWY8eOVVhjyrZt2+Ls2bPo168fBgwYgPPnz8Pa2loh56bert7JSdXVG8eOHYuVK1fis88+k22bN28etm7ditOnT+OLL75QbJRUnaSfoq9fv04/vTUTPB4P3bp1Q1JS0jtXnzU2NoahoSEeP34sK25NSmJqMjgcZuC0oIDg+nVm9g6fz9y+KS2tvQOsrq4uysvLlf8ka5DOzKmahBkbi3DpUgAePDiPUaN+hUAwVXa8tFjXza3xxbrSqcvBwcHV6k3s7OxUsgIy9aaWLVsiKCgIQUFBEIlEOHr0KCIjI7FixQosXrwYrq6ushWU5e3kK9WuXbtqCcq5c+fobTwVkKvm5MSJExg2bNgb24cOHUpn6qiYg4MDdHV1ad1JM+Pq6lrv23nW1tZ4/PjxG8WtpaXMyElqKkFeHlOP0r07oKVVd1Gprq4uysrKlPOk6lC1AVpgYCC8vb3Rvn17dOvmjKdPr8HH5yR0dKaiqAiorIRSinWb+8wadWZkZIQJEyYgMjIST58+RWRkJBwcHLB27Vo4OjrC2dkZy5cvR2pqqtxdfdu3b4+zZ8+ipKQEAwcORH5+voKfBVWTXL+2bdq0weHDh9/YHhMTgzZt2jQ6KKr+tLW14ezsTJOTZsbV1RXp6el4/fr1O49t166drElb1eJWkYhJTvLyJODzgd69gYKCw9i9uz+A2otKdXR0VJqcVG2AJi38/emnn/D06VM8fvwYq1evxLp1XhAImHqZO3egtGJdOrNG/RkYGGDs2LGIiIjA06dPceTIEbi5uWHjxo1wcXGBg4MDvvrqKyQlJTU4UenYsSPOnj2LZ8+eYfDgwXj27JmSngUFyLnw34oVKzB9+nTEx8fLak4uXbqE48ePY8eOHQoNkHo3gUCAq1evsh0GpUICgQCVlZVIS0t7Z6Moa2tr3LlzR/azoyPT38TcnIP4eKBTJwJ3d2bEJC+vDA8fnsfr10XQ0zOBgQHw+PG/RaWqvK0jFosRHBwMHx8fxMTEgMvlIiIiAvPmzUPv3r2hr6+PDRs24NNPP8XixTxkZTFxGhkxt3LYmN5MqQ89PT2MGjUKo0aNQllZGc6cOYPIyEhs27YNa9asQYcOHWS3fnr16iWrwXqbLl264MyZM/Dy8sKQIUNw+vRptGzZUvlPphmS69d36tSpspU5o6OjERUVBRMTE1y8eBFTp05VcIjUuwgEAty8eZOVWgCKHS4uLuByufW6tVN15ESKywXatmX+GLdsSVBaymw3NrYBABQXZwN4s6hUlSMnVRugcTgcfPvttwgMDMT48eNx6tQpLFu2TNYAjctl6mVcXJjvNDGhqtLV1cWIESOwa9cu5OXl4eTJkxg0aBB+/fVX9O7dG3w+HwsWLMDFixdr6R1UXdeuXXHq1Cncv38fw4cPp8u1KIncv8K9e/fG/v37ce3aNSQlJWH//v3o3bu3ImOj6snV1RUVFRVvLB9PNV36+vro3LlzvW7nWVlZ4cmTJ9i/fz/i4+MhFosB/Dtbx96eyFY9NjFhqkeLirJr7QCrypGTqg3Qnjx5go0bN2LVqlXYvXs3dHV1G90AjWqetLW1MXjwYPz000948uQJ4uLi4Ovri4iICPTt2xc2NjaYO3cuzp07J/tdqcnV1RUnT55EWloaRo4ciZKSEhU/i6ZPruSEx+PVWhBUWFhI78OyoFu3buBwOLTupJmpT1FsdHQ0Vq5cCbFYjA8//BDe3t6y1XS5/x9eGDaMyIpkJRJLcDg85OZm1VpUqsqRk6oN0Nq1a4c7d+7gq6++kiVVtAEa1VhaWlrw8vLC1q1b8fjxYyQkJGDcuHGIiYmBl5cXrKysMGvWLJw+fRqVlZXVHturVy/89ddfuHbtGkaPHo1Xr16x9CyaJrmSk7oKicrKyurd4p5SHENDQ3Tu3Jk2Y2tmXF1dcf369TqHoaXFpF26dAEAnDt3DkKhULaabmJiIgBm5ERaJPv8uRZatGiH/PysWotKVTlyUrMBmqmpqWxf1QZoNjae71xdmaLehcvlom/fvti0aRMePnyIS5cuYfLkyThx4gQGDx6Mtm3bYsaMGfjrr79kvwMeHh44evQoEhMTMXbsWJXPZGvKGlQQK23bzOFwsGPHDhgaGsr2icVinD9/vtFzyt/lxx9/xPfff4+cnBx07doVmzZtqrO/R3x8PLy9vd/Ynp6ervQ4VU0gENCRk2ZGIBBAJBIhIyMDHTt2rLavajHpTz/9BCsrKzx79gx+fn6IiYmBn58ffvnlFwDMhw1pkWxWFpCaaov27bPw5Zdv1m7o6uri+fPnKnl+VRugjR49GsOGDUOLFi3w6tUrHD9+HEePHsXEiZFYsYL3ztWVKaohuFwuevfujd69e+O7775DUlISIiMjcejQIezcuRMmJiYYPXo0AgICMHjwYBw5cgS+vr4YP348Dh06JOsiTMmvQcnJxo0bATB/zLZv317tFo6Ojg74fD62b9+u2AirOHjwIObPn48ff/wR77//Pn766ScMHz4caWlpsH1Lp6Xbt29X655qZmamtBjZIhAIEBsbC4lEIhuup5q27t27AwCSk5PfSE6kxaQRERGwsLCAtra2bC0S6Wq60oZi0pEXaVGpg4MNHj/OqrWoVEdHR6WF1/7+/li4cCE2btyI2NhY2XYeTwvdui1Eebk/bGyYVvUlJbU3jqOoxuBwOHBzc4ObmxtWr16NlJQUREVFITIyEr/99huMjIzg4+OD4OBgfPfdd5g0aRLCw8OhpSXXZFjq/xr0r5eRkQEA8Pb2RnR0NFq1aqWUoOqyYcMGTJ8+HTNmzADAtNQ/ceIEtm3bhtDQ0DofZ25u3uSne7m6uuLly5e4f/8+OnXqxHY4lApYWFjA0tISSUlJGDt2bLV9VYtJuVwufv31V/To0UO2X1pMCrx5m9bW1hZCobDWa6q6CVt0dDTWr1+PkSNHYvjw4WjRogVKSl5h69a/cP36enTu7A5jY6YRmrRxXFoa0ziuSxc6a4dSLA6Hg27duqFbt25YsWIF0tLSEBUVhaioKEREREBXVxeRkZF48uQJjh49ChMTE7ZD1lhypXZxcXHVfhaLxUhJSUH79u2VlrCUl5fjn3/+weLFi6ttHzJkiOzeeV0EAgFev34NJycnfP3117Xe6pEqKyur9se3uLi4cYGriLSNfVJSEk1OmpG6imKrFpO6u7tj0qRJ1fZLi0mB2pOTx48fQywWv1HgrsqRk9r6nABMbcnly7NQVjYSx4/PBUDQocMg6OmZvNE4Th3WAKKaLicnJzg5OeGbb77B3bt3ERUVhV9++QUXL15EmzZt4OPjg3HjxsHHx4cmKg0k1+eK+fPnY+fOnQCYPyD9+vWDm5sbbGxsEB8fr8j4ZAoKCiAWi2FhYVFtu4WFBXJzc2t9jKWlJX7++WdERUUhOjoaXbp0wcCBA3H+/Pk6rxMaGgoTExPZl42NjUKfh7KYmZmhXbt2tCi2makrOaltNV0paTFp27ZtAdSenFRUVCAvL++N86py5KRqnxMul4u//voL48ePh49PL/z+uzkyM49DJHqCQ4cCkJd3XfY4RaxGTFEN1alTJyxevBj3799HWFgYxGIxLl26hA8//BDm5ubw8fHB7t27aWfZepIrOTl06JDsfveff/6JzMxM3Lp1C/Pnz8dXX32l0ABrqtnFjxBSZ2e/Ll264OOPP4abmxv69OmDH3/8ESNHjsT69evrPH9ISAiKiopkX9nZ2QqNX5loUWzz4+rqisePH+Pp06fVtkuLSWNjY+Hn5ydr/S4UCuHn54fY2FjMmjULQO3JCQBkSXvWV6HKqcRVb00BQFFREQoKCtC1a3c4OX2BwYOZD0gjRmyDre2/RfGKWI2YohpjwYIF+Pnnn5GXl4fp06dj3bp1KC4uxkcffQQLCwsMHToUv/zyyxu/t9S/5EpOCgsLZZ+6jh07hnHjxqFz586YPn06UlJSFBqglKmpKXg83hujJPn5+W+MpryNu7s77t69W+d+XV1dGBsbV/vSFDQ5aX6qrkpdk3Q13ZSUFHh4eMDY2BgeHh5ITU1FZGQk+vdn1tCpmZxYWVkBACIiIqo1bROLxcjJyUFRUVG17cpS9dYUAEyYMAFnzpxBRMQOjBv3FSormYpXMzMH2QeU2hrHURQbPv74Y2zduhU7d+5EXl4ezp07h8ePH2Pz5s2oqKjArFmz0LZtWwwcOBDbtm2r8w5AcyVXcmJhYYG0tDSIxWIcP34cgwYNAgCUlpYqrQmbjo4OevTogVOnTlXbfurUqWrLmL9LUlJSk23a5Orqiry8PNoxsxnp2LEjDAwM6ryd97bVdP99Q/83OYmOjkbPnj0BMK0DpE3bFi1aBHt7e+zevRvPnz+v1sxNWeq6NcXlAqNHS/DgQSgMDOxgYuKptNWIKaox5syZg7CwMKxduxYrV66EpaUlZs+ejbNnzyI3Nxfbt2+HlpYW5s6dCysrK/Tr1w9btmx5Y7mJZonIYdmyZcTExIQ4ODgQW1tb8vr1a0IIITt37iTu7u7ynLJeDhw4QLS1tcnOnTtJWloamT9/PjEwMCCZmZmEEEIWL15MgoKCZMdv3LiRHD58mNy5c4ekpqaSxYsXEwAkKiqq3tcsKioiAEhRUZHCn4+iPXjwgAAgR48eZTsUSoX69OlDBg4cSMLDw0lcXByprKys1+POnTtHAJDbt28TQgiJiooiHA6H+Pr6Ejs7OzJr1iwiFApJz549CQDSs2dPMm/ePGJkZESEQiHx9fUlHA6nQb9PDVU1psTERFJcXEwSExNl1w4MjCJBQYSMG0dIUBAha9YQkpamtHAoSi6rV68mAEhoaGit+wsLC8mvv/5KRo4cSbS1tQkA4u7uTtavX08yMjJUG6yCyfseKldyQgghhw4dIhs2bCDZ2dmybbt37yYxMTHynrJefvjhB9K+fXuio6ND3NzcyLlz52T7pkyZQvr37y/7ed26daRjx45ET0+PtGrVivTt27fBb9yalJxIJBLSsmVLsnr1arZDoRpALCYkI4OQGzeY72Jx/R8bFRVFjIyMCADZF5/Pr1fCcP78eQKApKenk8rKSsLn84mvry8Ri8Vk2LBhZPTo0bLtFhYWhM/nky1bthAdHZ3/xy2WJTL1TYjkERUVRfh8frXnaGdnR6Kiohr1b0dRqrR06VICgGzcuPGtxz1//pzs3buXjB49mujq6so+GISGhpK7d++qJlgFUnly0lxoUnJCCCFeXl4kICCA7TCoekpLI2T1alLt0//q1fX79C8dVejWrRvhcrkkPz+/QSMaCQkJBABJS0sjcXFxBAARCoWEEEI+/vhj4ubmJtv+888/EwAkODiYACASiYQQQkhiYiIBQOLi4hr7T/FWlZWVJC4ursGjQxSlLiQSCVm0aBEBQH788cd6Paa4uJgcOHCABAQEEH19fQKAdO/enXz77bckTUOGCOV9D633XdkDBw7U+1ZRdnY2Ll68WO/jKcWhRbGaIz0d2LKF6Wpqaso0DTM1ZX7esoXZX5eqPUB++eUXSCQSZGZmwt3dHTExMfDx8cHChQvfWrQq7RtCCHljZkyPHj3QsWNH2XYfHx8ATLHs3r17ZXUqqloZmMfjwcvLCxMnToSXlxddYJTSOBwOB2vXrsXnn3+O2bNnY9euXe98jJGRkawl/tOnTxEVFQUnJyesW7cOTk5O6Nq1K5YtW4aUlJQ617zTVPVOTrZt2wYHBwesW7cO6bX81SwqKsKxY8cQGBiIHj160LncLHF1dcX9+/dRVFTEdijUW0gkwOHDQEEB09XU2Bjg8f7tclpQwHQ5rWshu6o9QFxcXMDj8WRFsdL29BkZGUhISKgzhqoFsTVnxsycORO///67bLu0dbybmxs+/PBDWWJDVwamqPrjcDjYuHEjZs2ahRkzZmD//v31fqy+vj78/f0RHh6Op0+f4siRI+jZsyc2b96Mbt26wcHBAUuWLMG1a9eaRKJS7+Tk3LlzWL9+Pc6ePQtnZ2cYGxujU6dOcHFxgbW1Ndq0aYPp06eDz+cjNTUVvr6+yoybqsPbppZS6iMrC7h1C7CxATgcICLCF+fOrQSAN7qc1qbqSEeLFi0wbtw4GBgYyPbXZ0SjanJS18wY6fZvvvkGfD6/2iKbVVcGrmvxTYqiquNwOPjhhx8wZcoUTJ48GYcOHWrwOfT09DBq1Cjs2bMH+fn5OHbsGPr27YuffvpJNuq5aNEiXL58WWMTlQa1r/fx8YGPjw8KCwtx4cIFZGZm4tWrVzA1NYVAIIBAIKCLzrHMwcEBurq6SE5ORr9+/dgOh6qDSMR0MZXmEzyeLh48OI3+/ZcCYLY/flx3l9Oa7ekjIiKq7a/PiEbV5KTqCsB+fn4ICQmBs7MzUlNTYWpqiqtXr6Jnz564cuWKbHtoaChiY2MRGRlJb7NQVANwuVzs2LED5eXlCAwMhI6ODkaPHi3XuXR0dDB8+HAMHz4c27dvx7lz5xAZGYndu3fj+++/h42NDcaOHYuxY8fCw8OjzvdoiYT5MCQSMQ0MbW3ZnY7PIZqaVqlIcXExTExMUFRUpDEN2Xr16gVnZ2f8+uuvbIdC1SEzE1i6lKkxMTYGrlzZipMng/Hll8+hra2PoiKgsBBYubL29WHEYjHs7e3h4uJSbd0ZgBnR8PPzQ2pqKu7evVtn4nDlyhX07t0b169fR7du3QAwfU6Cg4ORmZkpO87Ozg4BAQE4dOjQG9vXr18Pf39/BfyLUFTzU1lZiQkTJuCPP/7AkSNHMHz4cIWdWywW48KFC4iMjERUVBRycnJgaWkJf39/BAQEwNPTU/a3IT2duc186xbzoUlPD3BwAMaMafzq3vK+h9Lk5B00MTn55JNPcOXKFbrOjhqTSIC1a5niVycn4OnTVGzb5oLJk8+Azx+AtDTAzQ348su6P71ER0cjICAAPj4+1UY6qo5ovC1x+Pvvv/Hee+8hKSkJrq6usu1isRgJCQmyP2bSP2J1bacoSn7l5eUICAjAyZMncfToUQwcOFDh15BIJBAKhbJEJTs7G+bm5hgzZgx69RqLy5e98Py5NmxsmFHbkhIgO5v58DRvXuMSFJqcKIkmJic//vgjPv/8c7x8+RK6urpsh0PVQTpbp6AAaNdOgl9+sYCLy6ewsVlZ7z8KdY101GdE459//kHPnj1x7do1Wa0SRVGqV1ZWhtGjR+P8+fM4fvy4Um/JE0Lw999/IzIyEpGRkcjIyICOTms4OfnByWksOnQYBB5PB4SgXh+S3oUmJ0qiicmJUCiEh4cHfdPRAFWHU8+eDUBFRT7mzz8PP7/6f1qRd0Tj2rVr6NGjB/755x+4ublV26fs+89VY7awsET79p4oLeWpxb1uimLDq1evMHLkSPz99984efIk+vTpo/Tfw4wMgrlzk/H0aSTu3TuEZ8/uwtV1GkaPZqY5v+v2cn3I+x7aoIJYSjN069YNHA4HSUlJNDlRc46OTH+TrCzA2toLYWHBmDu3FIaG+vU+h7QHSEPVtrYOoNz7z0Dtoz0GBnx06xYGe3t/hV6LojRFixYt8Oeff2LYsGEYNmwYduw4g7t3eyrt9xAAXr7kQF9fgCFDBBg2bBXy81PB4fyb/byrMF+ZGpWDlZeX4/bt26isrFRUPJQCGBgYoEuXLrQZm4bgcplPJYGBXigvL8eVK5dUct3akpPGNIarD2mdjIuLCyIihJg+XQRvbyEsLFwgFAbg+fNohV2LojSNgYEBjh49Cj7fEUFBQxAXd10pv4dSRkZM0lNSwvw9sLBwgbl5V9n+khJmv5GRYq7XEHIlJ6WlpZg+fTr09fXRtWtXZP2/GcO8efOwdu1ahQZIycfV1ZUmJxrGyckJpqamiI+PV8n1aiYnjW0M9y5Vu9pGR8fg1q1uKCoyhKenOz78MAadO/vgwoWFcHAQN/paFKWpDA2NMWrUcbRoYYfExEEoK0tT6O9hVba2zGhMdjZQs8CDEODRI2aUxta28ddqKLmSk5CQEFy/fh3x8fHQ09OTbR80aBAOHjyosOAo+XXv3h3Xrl3D/v37ER8f/9Y25pR64HK56N+/P2vJSc3GcMnJe5CXl/L/Y9/dGO5dqna1vXmzCCtXGuHMGQdER09EYuJ6ODl9gBcvMpCdndDoa1GUpsrKAh4+bIkJE07C0NASv/02EIWFdwAo5vewKi6XuU1kasoUvxYVAZWVzPe0NGa7nx87NWByXTImJgZbt25F3759ZX/gAOaT3/379xUWHCWf6OhobNq0Ca9evcKHH34Ib29v2NvbIzo6mu3QGkUiYfqDpKQw35vip2ovLy9cvnwZpaWlSr+WtBPsiRMnEB8fjxcvxNUaw/3111w8eHBKdryBAXPvW977z1W72r56pYUePXagY8fBKCrKxrlzK1FUxPy1FYlyGn0tTSQWixEfH4+IiAj6gaIZkzZoNDVtg8mTT0NPryWuXv1Jtl/RvxuOjszMQIGAKX69c4f57ubW+GnEjSFXQezTp09hbm7+xvaSkpJqyQqletJ7+oMHD8bJkyexZ88edO7cGWvWrEFAQMA7e1+oK2UXaaoLLy+m7uTSpUsYMGCA0q4THR2NuXPnAgCWLmW60lpb89GxYxisrf1hbAxoaemisrJM9pjG3n+u2tW2bVt3ODpOg6npNBgbA4RIkJl5HgBgZGTJ6r1uNtRWJMzn8xEWFqaRv6+U/KrWgRgbm2PatAvQ02sp26+M342qhfnq0iFWrkv36tULR48elf0sTUh++eUX9OnTRzGRUQ1W9Z7+X3/9BWtra6SnpzdopVp1pOwiTXUirTvZvXu30j5BSxPYzp07AwBOnToFoVAIV1cXnD8fgEuXokEI01JfLGaSE0Xcf666fo+1teSNe92XLm1Ay5Z2sLHxZPVet6pVLRIWCoUQiUQQCoVwcXFBQECAxo94Ug1Tsw5EX78NuFymNYAy60CkhfkuLsx3tqfzyzVyEhoaimHDhiEtLQ2VlZXYvHkzbt68CaFQiHPnzik6RqqepPf0IyIiwOVyqxXFSleq9fDwQEJCglxTT9lQs0hTOjAnLQ5LS2OKw7p0Yf+XSRFiYmJQWlqKvXv3Yu/evQAU+wm6agIbGhoqWzjQ3d0dR47EYNAgP1y9uhBWVqPB5eqioqIMRUXMH8TG3n+uun6Pv78fxo9fhAsXihEZ+TeePz+BnBwhfHyicOsWj9V73apU9f9HTEyM7IOe9AOFn58fFi5ciNGjR9NuvM2EtA4kO5v5+2Zt/W/XVkX8HmoKuZ6eh4cHLl68iNLSUnTs2BEnT56EhYUFhEIhevTooegYqXqqek8fAL799lts2LBBtr8+K9Wqm6pFmoSIERv7KdLTmU+Sii4OY5v0E3SHDh2gra2NvLw8hX+CrlqUKn2zkxbEcrlcrF4dApEoAy1bJgDQxdOnZQq9/+zv74/IyEhcunQJH37oib/+Gom0tOXIyRFCV9ccIhH797pVqer/j9zcXLi4uODixYsA/v1AkZGRgYSEBJYjpVRJXetAVEnuJmwuLi7Ys2ePImOhGqnmSrVV10uRbq96nCaounovl8vD06c38fz5fTg6MqMIbDYJUqSqn6BXrVqF7t27IzU1FQMGDFDoJ+iqCezjx48BVO9zIk1g+/d/jNOny9GmTSpGjIjH2LGe0NZW3Cf3goICuLl5wN19KNzceuH1a2Ps3r0O588HYM6cSNn/36au6v8PfX19iEQihIeHo2/fvrLtVY+jmg91rANRJbmeJo/HQ35+/hvbCwsL6dAji6re05fUmMoikUgQGhoKOzs7eHp6shRhw1UtDgMAV9dpePDgNIqKsgGw2yRIkap+gnZ2dq7W70SRn6CrJrB8Ph8pKSnVWtdLE9j//CcYDx/ex6VLpzBxojc6d1bMbC+xWIx584LRpYsPnJwSkJY2C5s3v48dO95H69YxaNfOBzNnLkRqqmbVRcmr6v8PLpeLCRMm4Pfff0dFRYVse9XjqOZF3epAVEmup1rXcjxlZWXQ0dFpVECU/KT39GNjY+Hn51etuM7Pzw+xsbFYv369RiWQNYvDnJwCoK3dAtev/8Z6kyBFqvoJmsvlYvv27Rg3bpxsv6I+QVdNYLW1teHs7AyD/88dlkgkmDdvHgCm6N3V1RUTJ05U6K2l336Lx+PHmeByeyIvLx6XLw/D/fsf4OXLChQUcOHiEoLnzzOwZElCkyp0rkvNDxSBgYEoLCzE6dOnNfYDBUUpQoNu62zZsgUAMztnx44dMDQ0lO0Ti8U4f/48HBwcFBsh1SDSe/rBwcHw8PCQbbezs9PIacRvFocZwcEhAElJu9G69RKYmXGaRHFYzVtyY8eOrbZfUZ+gqxal+vn5ISQkBM7OzkhNTcWaNWtw9epV9OzZE0eOHMHQoUMhFosVVpwZGRmNuXNnAQDS0pYhLQ3gcs0gkdzAy5ezoaX1M0pKmCQsLy+nSRU616Xm/4/FixfDwcEBmzdvxrZt2xAbG4vIyEiN+kBBUQpBGoDP5xM+n084HA6xsbGR/czn80nnzp3JkCFDyKVLlxpySrVXVFREAJCioiK2Q2mQyspKEhcXR8LDw0lcXByprKxkO6RGSUsjZPVqQoKCCOnf/ywBQGbOvEDS0tiOTDEqKysJn88nvr6+RCwWV9snFouJr68vsbOzU9j/x6ioKMLn8wkA2ZelpSUBQIRCISkvryQCgTtxcupBIiLiSHl5JUlMTCQASFxcnFzX43A4xNTUgwAgY8acIW3afEE4HD3Z9W1tw0jHjsw1xo6NI0FBhGRkKOTpqr3a/n+0b9+eREVFsR0aRTWKvO+hDUpOpLy8vMizZ8/keajG0dTkpCkSi5k3q+RkMbG25pOPPprOdkgKJX0D9/X1JYmJiaS4uJgkJiYSX19fwuFwFP5GVTOB3bdvHwFA1q3bR1q1qv5G2aoVn6xbx+wPDw9v8HX4fD7p39+XjB1bToyN2xMTk04EAGnXLpi0bu1DtLTaEDe3NGJo6EuMjOxISEglGTeOkBs3FPqU1Zr0/8fGjRsJABIREcF2SBTVaCpNTpoTmpyop2XLlhEjIyPy8uVLtkNRqNo+QdvZ2ankE3RcXNz/r8khlpa+ZMIEIVm0SEQmTBASS0tfAnDkGjmRnjc6WkiCggjp0GEUAUB0dBxJ166JxMXlFAFADAz6EIBDRo6MIvPnk2Y1clJT7969yejRo9kOg6IaTd73UA4hdVS3vsOjR4/wxx9/ICsrC+Xl5dX2Ve2toemKi4thYmKCoqIiGBsbsx0O9X8ZGRno0KED9u7diw8//JDtcBRKLBYjISEBOTk5sLS0hKenp0pqDl6/LoeBgQG0tdvgP/95BB7v35I0sbgS339vjYqKQpSUlEBPr/6F7xEREQgMDERRkQhbtxri+PHTKC2Nxs2bf+H160zZcTyeGRwdt2PMGH+kpzM9Hb78smnXnNRly5YtWLhwIfLy8tCqVSu2w6Eoucn7HipXn5MzZ85g1KhRsLOzw+3bt+Hs7IzMzEwQQqpNS6QoZbGzs4OXlxd2797d5JITHo/HSgffmJhESCSVKCvLw759Q/D8+QOMHXsAAMGFC6EoK8uTHTdhQv3jkxbxpqWlYswYd2RnD0JBwSA4OIhx+XICCgvP4fnz5bCzi4Cb20CkpzefLph1+eCDD/DFF18gKioKM2bMYDscilI5uX71Q0JCEBwcjNTUVOjp6SEqKgrZ2dno379/temPFKVMU6dOxdmzZ/Hw4UO2Q2kSHj5kpin7+e3D8+f3UVT0ELt29cGuXR7Iz0/FqFH7qh1XX1Wny3bpIpF1vuRyeejQoR/E4n+go2OHDh28QEjz6oJZl7Zt22LAgAGIiIhgOxSKYoVcyUl6ejqmTJkCANDS0sKrV69gaGiIlStXYt26dQoNkKLqMnbsWOjr6+O3335jO5QmoX17ZoRDT68j5s17gJEjt2P06N2YMiUOc+fehZ5eh2rH1VfN/jsvXggxZ44Ifn5ClJX5QSSKxaZN67F+PQ8rVzK3cppzYiIVGBiIuLg4PHnyhO1QKErl5EpODAwMUFbGrFZqZWWF+/fvy/YVFBQoJjKKegdDQ0N88MEH2L17d52NAan6GzvWE61a8XHu3BoAHPTsOROurlPA53sB4CAhIRStWtlh7NiGNwST9t9JSUmBh4cHWrY0xtixHsjKSkVkZCQ+/dS/WXbBfJsxY8ZAW1sbBw8eZDsUilI5uf4MuLu7yxanGjlyJIKDg7F69Wp89NFHcHd3V2iAFPU2U6dOxYMHD+jCaAqgrc3D8uVhyMmJxY4dfkhPF6KkRIT0dCF27PBDTk4sli9fL/caO/7+/rh37x7i4uIQHh6OuLg43L17V+MaA6pKy5YtMXLkSHprh2qW5Jqt8+DBA7x8+RLdunVDaWkpFi5ciAsXLsDe3h4bN25E+/btlRErK+hsHfUmkUjQqVMn9O/fH7t27WI7nCZhy5ZoLF8ejOfPM2XbWrWyw/Ll6zFvHk0kVCkyMhLjxo3DnTt30KlTJ7bDoagGk/c9tMHJiVgsxoULF9CtW7dmMcWNJifqb+XKlfjuu++Qm5tbbUkFSn4VFWJERSXg4cMctG9vqfBVian6efXqFSwsLLBw4UIsXbqU7XAoqsHkfQ9t8G0dHo+HoUOH4sWLFw19qEL8+OOPsLOzg56eHnr06PHO4fxz586hR48e0NPTQ4cOHbB9+3YVRUqpypQpU1BSUoKoqCi2Q2kytLV5mDDBC19+ORETJnjRxKSBxGIx4uPjERERgfj4eIjF8q2y3KJFC4wZMwbh4eG0ropqVuSqOXFxccGDBw8UHcs7HTx4EPPnz8dXX32FpKQkeHp6Yvjw4cjKyqr1+IyMDIwYMQKenp5ISkrCkiVLMG/ePPom1sS0b98eAwYMwO7du9kOhaIQHR0Ne3t7eHt7IzAwEN7e3rC3t5d7RefAwEDcvn0bycnJig2UotSZPO1oT5w4QVxdXcmff/5Jnjx5QoqKiqp9Kct7771HZs2aVW2bg4MDWbx4ca3HL1q0iDg4OFTbNnPmTOLu7l7va9L29Zrht99+IwDIgwcP2A6Fasaqro8kFAqJSCQiQqGwUesjVVRUEDMzM7Jw4UIlRExRyqXS9vXcKnP9OBxO1UQHHA5H7iHMtykvL4e+vj4OHTqEMWPGyLZ//vnnSE5Oxrlz5954TL9+/SAQCLB582bZtsOHD+ODDz5AaWkptLW133ldWnOiGUpKSmBpaYkFCxZg+fLlbIfTYGy1rKcURywWw97eHi4uLoiJian2d1IikcDPzw+pqam4e/dug//ffvbZZzhy5AgePnxY7bwUpe5U2r4+Li5Onoc1SkFBAcRiMSwsLKptt7CwQG5ubq2Pyc3NrfX4yspKFBQUyNpqV1VWVibr4QIw/7CUakgkQFYWIBIBRkaArW39e14YGBjggw8+wJ49e7B06VK1/ANeVwISHR2N4OBgZGZmyo7l8/kICwuj02w1SEJCAjIzMxEREQEul4v09HRUVFSgW7du4HK5CAkJgYeHBxISEhq8PEFgYCB++OEHXLhwAf369VPOE6AoNSJXctK/f39Fx1FvVUdqgH9HaxpyfG3bpUJDQ7FixYpGRkk1VHo6cPgwcOsW8Po1oKcHODgAY8bUv1vo1KlTsXPnTpw/f56VtWnepq4EZNy4cVi/fj18fHwQEREBZ2dnpKamYs2aNQgICEBkZCRNUDRETg7T1t/Z2RkA8MUXXyA5ORnnz59H586dZdulxzVEnz59wOfzER4eTpMTqllQv4+XdTA1NQWPx3tjlCQ/P/+N0RGptm3b1nq8lpYW2rRpU+tjQkJCUFRUJPvKzs5WzBOg6pSeDmzZAiQlMQu+denCfE9KYranp9fvPO+//z7s7e3VrjA2OjoaAQEBcHFxgVAohEgkglAohLOzM77//nv06NEDMTEx6NGjBwwMDODu7o6YmBj4+Phg4cKFSrlNSimedCQ2NTUVALB37160atUKgwYNQlZWlmx7bSO278LhcDBhwgQcOnTojVXgKaop0pjkREdHBz169MCpU6eqbT916hQ8PDxqfUyfPn3eOP7kyZPo2bNnnfUmurq6MDY2rvZFKY9EwoyYFBQATk6AsTHw6lUejI2ZnwsKgJgY5rh34XA4mDp1Kg4dOgSRSKT02OtDLBYjODgYPj4+iImJQZs2bVBWVgZ3d3d88cUXAIDHjx+DEILQ0FB4eXlBIpHIbgNkZGTQ7rcaouoChxKJBGZmZjh9+jS0tLQwcOBALFu2DHZ2dvD0bHj7f4C5tfPs2bM3/qZRVFOkMckJACxYsAA7duzArl27kJ6eji+++AJZWVmYNWsWAGbUY/LkybLjZ82ahYcPH2LBggVIT0/Hrl27sHPnTixcuJCtp0DVkJXF3MqxsQE4HODevRPYvNkOT578Aw4HsLZmRk7qmC0OiQTIzARSUpjvkyYF4dWrV4iMjFTl06iTtA5h8eIlyMri4qOPPoev7zhIJEBeXh4AZpj/zJkz2Lp1K5ydnWX1Mo25DUCpXs0FDoVCIYyNjfHdd98hOzsbp06dwrJly+QudHZxcYGzszPCw8MVHDlFqR+5ak7YMn78eBQWFmLlypXIycmBs7Mzjh07JmuXn5OTU63niZ2dHY4dO4YvvvgCP/zwA6ysrLBlyxaMHTuWradA1SASMTUmBgbMz3y+F8zNuyIy8gN88sk1GBiY4PFj5riaaq9TsYW7+0Ds3r0b06ZNU+2TqYU0sThxwhl37rxEYuIZuLmtw9q1QLt2/w7vHzp0CAUFBbLRFACNug1AsUO6wGFwcHC1Ed127dqhuLgY27Ztg7+/P4yMjOQ6f2BgIFatWoWSkhIYSH9pKKoJkmsqMQBUVlYiPj4e9+/fR2BgIIyMjPDkyRMYGxs3qRbidCqxcmVmAkuXMjUm0n/e588z8NNPAnToMBBDhkTi2TMOVq5kVqyVktapFBQwoy4GBkBJCZCdDTx/Ho4zZybh3r176NixIxtPS2bPnnhMneoNb28hDAweIzY2ANOm3YdI1AGtW4tx+LANnj7NgY2NDXr27Clr1NXYqacUu2qbmZWcnIwBAwagR48eOHr0KFq0aNHg82ZkZKBDhw6IiIjAhAkTlBA5RSmWytrXA8DDhw/h4uKC0aNHY86cOXj69CkA4LvvvqO3TKgGsbVlZuVkZwPSNLlVKzuMHv0r0tOjER//Xzg6MsdJ1VanwuNBVqdiZOQHXV1j7N69h50nVSXO7GxP6Ovz8fjxGmRlxcDc3Bm2th3g5AQUFnKgo9MOAJCdnY1hw4bJimX9/PwQGxuL9evX08REA/F4PHh5eWHixInw8vICj8eTJSWXLl3CuHHjUFFR0eDz2tnZoU+fPvTWDtXkyZWcfP755+jZsyeeP39eLfsfM2YMzpw5o7DgqKaPy2WmC5uaAmlpQFERUFkJWFmNgb39fKSkLESHDleq9TuprU6lrKwYEokYDx/G4/XrIzA19cSuXXsgqU8lrZJkZQF37vAwYEAY7tz5E6mpB2Fp2QNlZSI8eiREUpIfHj/+B5aWNtDV1cXMmTNhbGwMDw8PpKam0mnETVDfvn1x+PBhnDx5EkFBQXLNxAoMDMRff/2FwsJCJURIUepBrts6pqamuHjxIrp06QIjIyNcv34dHTp0QGZmJpycnFBaWqqMWFnRHG7rNKb5maLUVj/SqVM5Dh3yRHFxHpKSkmSrYKekAN9+y0w5Li9/js2b7WBm5giRKAdFRQ+rnXfZsmWsdYytGqdQuBZnzoRU29+ypR1sbGYhJeVL/P777zAzM6MdYpuJ6OhojBs3Dh999BF+/vnnt/ZqqikvLw9WVlbYtm0bPvnkEyVGSVGNp9IOsRKJpNaM/9GjR3IXelHsUETzM0VwdGTexKsnSTqYPPl3CAQCTJs2DYcPHwaHw4GRERNnSQlgbNwKPXp8gsTE72Fm1hXTpwuhq+uM+/dTkJAwGCtWrEC3bt1YGYGoGmffvovh4OAPkegRXr7Mg5GRJUxMPHHkyCTY2naAv78/TUaaEX9/f/z666+YMmUKjIyMEBYWVu8ExcLCAoMGDUJERARNTqgmS67Px4MHD8amTZtkP3M4HLx8+RLLli3DiBEjFBUbpWSKan6mKFwuU/Tq4sJ853KZFYf37NmDI0eOYOPGjQCq16mIxWKkpR1Cq1YdUFh4G2KxGPn5hujfvw+WLFkCLpeLBQsWsNLIrGY9jalpZ9jZDYCLy0S0b++Fu3ezkZ0dif/8Z4HciYlYLEZ8fDwiIiIQHx9PG7ZpkMmTJ2Pr1q3YuHEjVq5c2aDHBgYG4ty5c3j06JGSoqModsmVnGzcuBHnzp2Dk5MTXr9+jcDAQPD5fDx+/Bjr1q1TdIyUEryrqLQhzc+UzdfXFwsXLsSXX36JS5cuVatTOX8+AS9eZMLXdw/atnXHgQMfwMAgH35+wJQpTM+bhw8fstLIrK56mqIi5udHjzbDyMgE06ZNlev80dHRsLe3h7e3NwIDA+Ht7Q17e3vZjB9K/c2ZMwdr1qzB8uXLsWHDhno/bsyYMdDR0cHBgweVGB1FsUeu5MTKygrJycn4z3/+g5kzZ0IgEGDt2rVISkqCubm5omOklKBmUWlKSjj++ednAKhX8zNVW7NmDd577z188MEHKCwshKMjMG8e0LYt00dEJHJFnz4HYWzcEkOH3oejI2BtbQ1vb28A7DUyk8YpEACFhcCdO8x3R8cXuH9/Bz777FO5+lXU1RLfxcUFAQEBNEHRICEhIVi8eDGCg4Pxyy+/1OsxxsbGsvWYKKpJItRbFRUVEQCkqKiI7VAU6sYNQsaNI+TrrwlZtoyQ996bS/T0WpIvv3xOli0j5KuvmP03brAd6b+ysrJImzZtyMiRI4lYLCaEEHLmTBwBQPbuFZKMDEIqKsTVHrNixYr/79/LQsT/EosJychg/j0zMghZu3Yd0dHRITk5OQ0+V2VlJeHz+cTX11f27/DvdcTE19eX2NnZkcrKSsUETymdRCIhs2fPJhwOh0RERNTrMVFRUQQAuXXrlpKjoyj5yfseKtfISWhoKHbt2vXG9l27dtHbOhqiarEmAPTtG4LKytcQCpm6jpISZr861Tfb2Nhg7969OHr0KNavXw8A6N+fWc/k99/XwNZWAi2tf1/SEokEly9fBpfLxc2bN9kKG0D1ehorq3Js2bIZQUFBaNu2bYPPJW2JL62pmTRpkmykhK7Jo5k4HA7++9//4sMPP0RQUBBiY2Pf+ZgRI0bA2NiYjp5QTZJcyclPP/0EBweHN7Z37doV27dvb3RQlPLVLNY0MrJEr15zcOnSRpSUFOLRI7zR/EwdDB8+HIsXL8aSJUtw4cKFWtczqdrI7K+//sKgQYOwb98+tSkWPXDgAJ48eYIFCxbI9XjpLSrp2jsJCQlISkqS7adr8mgmLpeLXbt2wdfXFwEBAYiLi3vr8Xp6evD390d4eDiIfI2+KUptyZWc5Obm1rreh7RPA6X+aivW7N37SxAiwZ9/fg9TU8DPT/X9Turj22+/hYeHByZMmICnT5/K1jNJSUmBh4fHG43MVq5ciUePHuHs2bNshw5CCNavX48RI0bAyclJrnNIf/eka++Ym5sjPz9ftp+uyaO5tLS0EBERgf79+8PX1xeXLl166/GBgYG4e/curl27pqIIKUo15HrrsbGxwcWLF9/YfvHiRVhZWTU6KEo1ahZrPn5shi5d5uHBg/8iMDBPpX1OGkL6B7y8vBxBQUGQSCTw9/fHvXv3EBcXh/DwcMTFxeHu3bvw9/fHe++9BwcHB+zevZvt0HH69GmkpKQ0apkHT0/mVtaaNWsgkUiqJScSiQShoaGws7ODp6enosKmVEhXVxfR0dFwdXXF8OHDcePGjTqP9fb2hoWFBW1nTzU98hS4rF27lrRp04bs2rWLZGZmkszMTLJz507Spk0bsmbNGnlOqbaaakFsVVWLNZOSComxsTH54osv1CamjAzm55pOnDhBOBwOWb169TvPt27dOqKnp0eeP3+u6FAbZMiQIUQgEBCJRNKo80RFRREOh0N8fX3JiBEjyHvvvUcSExOJr68v4XA4JCoqSkERq4fKykoSFxdHwsPDSVxcXLMo9n3x4gVxc3Mj5ubm5Pbt23UeN3fuXGJlZdUs/k0ozSPve6hcyYlEIiGLFi0ienp6hMvlEi6XS/T19cmKFSvkOZ1aaw7JSU3Lli0jenp65PHjx6xcPy2NkNWrCQkKYmYMBQUxP6elvXns119/TbhcLomPj3/rOR8/fky4XC756aeflBT1u12/fp0AIPv371fI+aKiogifzycAZF92dnZNLjGp7Xny+fwm9zxrk5+fTxwdHYmNjQ15+PBhrccIhUICgMTFxak2OIqqB5UlJ5WVlSQ+Pp4UFhYSkUhErly5QlJSUsjr168beiqN0ByTkxcvXpBWrVqROXPmqPzaaWmEzJpFSEAAIV98wUx1/vzzStKvXxwZODCc7N5d/VNzZWUl8fLyIpaWliQ3N/et5x4+fDhxd3dn7RP4lClTiI2NDSkvL1fYOSsrK8mnn35KWrRo0SRHFKqOEAmFQiISiYhQKGyyI0S1efToEbGzsyP29va1Tj2XSCTEzs6OfPzxxyxER1Fvp9KRE11dXfLgwQN5HqpxmmNyQgghq1evJjo6OnV+WlMGsZgZIQkIIGTpUkK++CKb+PvvJy1bvv1T85MnT4iFhQUZNGjQW9+cFyxYUO08qvwE/ujRI6KtrU3Wr1+v8HPv3buXACClpaUKPzebaD+Xf92/f59YWVkRFxcXUlhY+Mb+JUuWkFatWpGysjIWoqOouqm0z4mLiwsePHggb5kLpQHmzZsHY2NjrFq1SmXXrNq1FpBg1673ER09CaamTpg+XYjZs0UYNkyIjh2rd0G1tLTE/v37cebMGaxevbrWc0dHR2PDhg3Q1tZGUFCQSjqqSiRAZiazOvGqVf9FixYt8PHHHyv8OtKuzFVn7DQFNfu5REdHY8aMGbh161az6+fSoUMHnDp1Ck+ePMHw4cMhEomq7Q8MDMTz589x4sQJliKkKMWSKzlZvXo1Fi5ciNjYWOTk5KC4uLjaF6X5DA0N8eWXX+LXX39VWSIqEjErIxsYMFNuxeJycDg8EFKJtm1d0aqVIYyM3BEWFgMfHx8sXLhQ1rtk4MCBWLZsGZYvX/7GlGGxWIzg4GD4+vpixowZOHPmDFq0aAF3d3fExLx5LkVITwfWrgWWLgWWLhVh166f0K3bJ3j8uP5LhtdXU01OavZzef36NY4dOwZHR0f4+fnh5cuX1Y5r6pycnHDixAncunULo0aNwqtXr2T7unbtim7dutFZO1STIVdyMmzYMFy/fh2jRo2CtbU1WrVqhVatWqFly5Zo1aqVomOkWDJ79myYmpo2eMVUeVXtWpuVlYCXL3MxbNhGPHx4HocOjUNxcTn09AATk9o/NX/99dcYOHAgAgMDkZubK9te9RP4tGnT8OTJE5w6dQqAcjqq1lztubR0FyorX8LYeJ5SVnuWJid5eXmKPTHLavZzCQwMREZGBnbt2oU7d+5gyJAhAID79+9Dog4rVKpAjx49cPToUVy+fBnjxo1DRUWFbF9gYCCOHDkiS9ooSpPJlZzExcXJvs6ePSv7kv5MNQ36+voICQnB3r17cfv2baVfr2rXWpGI+TTs6joN48cfxv37p/DPPydlXWtr64LK4/Gwb98+cDgcBAYGykZCqn4C79mzJ7p27Vqt+6YiO6rWXO3Z0LASly9vRNeu49Gzp41SVns2MzMD0PRGTmr2cwGYHiDTpk3DjRs30KtXL+jq6uKbb76Bk5MTdu7cibKyMpajVr6+ffsiJiYGp06dQlBQkOx1PmHCBLx69QpHjhxhOUKKajy5kpP+/fu/9YtqOj755BNYWVlhxYoVSr9W1a61BQXMp+acnFSYmQ3D0KF30a2bj6xrbV1dUC0sLHDgwAGcO3dOFnPVT+AcDgcJCQnV1oBSZEfVmqs9p6VFoajoIfr0CVbaas/a2tpo3bp1k0tO3rY0gb+/P65evYrw8HBcvHgRjo6O+Pjjj8Hn87Fu3Tq8ePGC7fCVasiQIThw4AAiIyMxc+ZMEELQvn17vP/++3StHappkLcC9/nz52T9+vVk+vTpZMaMGWTDhg3kxYsX8p5ObTXX2TpVbdu2jXA4HJKSkqKS66WlEfLtt5XEwIBPLC19yYcfismaNf/2OanPTI1Vq1YRDodDTpw4odJZHzVXe3Z2nkDs7AaSZcuIUld7dnBwIPPnz1foOdWl8Vl9+7mkp6eTGTNmEB0dHWJkZESCg4NJdnY2KzGryp49ewgAMn/+fCKRSMgPP/xAtLS0yNOnT9kOjaIIISqeSvz333+T1q1bk3bt2pExY8YQPz8/Ym1tTdq0aUP++ecfeU6ptmhyQkhZWRnh8/lk7NixKrumWEzItm1Mj4uBA33JhQuJpLi4uN5dUMViMRkyZAgxMzMjjx8/rtYvIzGxYedqiIwMpmncF18wycjSpWLyn/88lSUn8+cz+zMyFHI5mf79+5PAwECFnU/dGp81JFF68uQJCQkJISYmJkRLS4tMmTJFZYk1G7Zu3UoAkGXLlpH8/HzC4/HItm3b2A6LogghKk5O+vbtS6ZOnUoqKipk2yoqKsiUKVOIp6enPKdUWzQ5YezcuZMAIElJSSq9bmO6oObn5xMrKyvi6elJKioqVNJRtWavFmlSwiQqzPY1a2pvx98Y48aNI4MGDVLIuZpK47Pi4mISFhZGrK2tCQAycuRIcu7cuUYvHaCOQkNDCQASFhZGhg0bRvr168d2SBRFCFFxcqKnp0fS09Pf2H7z5k3SokULeU6ptmhywqioqCD29vbE19dX5dduzO2F8+fPEx6PR0JCQhp9rvqq2uV2/nzmVs78+czPs2bV3oa/sebMmUO6devW6PPUvAVWWFhIbvz/HpSmNj4rKysje/bsIV27diUASO/evUlUVJRGPYf6CAkJIQDItGnTCACSlZXFdkgUpdrkxNzcnJw4ceKN7cePHyfm5ubynFJt0eTkX9JOpJcvX2Y7lAaRfqo8duyYyq5Z2/pAVetmFG3lypXEwsKi0eeJi4sjAIhQKCSEMIvKde7cWbY/MTFRY9dxkUgk5OjRo6R///4EAOnUqRP56aefyKtXr9gOTSEkEgmZM2cO4XA4RFtbm3z33Xdsh0RRqk1O5s6dS6ytrcmBAwdIVlYWyc7OJhEREcTa2pp8/vnn8pxSbdHk5F+VlZXE0dGRDB06lO1QGkQsFpMRI0aQNm3aqPTTZH1WVlaU7du3Ey6X+0bBb0OFh4cTAEQkEhFCCNm1azcBQC5efEEyMgh58aKYACDh4eEKiJo9ly9fJgEBAYTD4RBzc3OyatWqWtvCq/L/oSKIxWIyefJkwuFwSIcOHdgOh6JUm5yUlZWRefPmER0dHdmqxLq6umT+/PlNbgFAmpxU9/vvvxMA5MKFC2yH0iBPnz4l1tbWxMPDg7x69UotZqEoUnR0NAHQ6FkaVUdO0tIImTcvlQAg/fvHkaAgQmbO1NyRk9rcuXOHzJo1i+jq6hIDAwMyf/582XpSVUe/AgIqyeDBceSDD95cfFLdVFRUkN69exMA5Ndff2U7HKqZk/c9lEMIIfJOQy4tLcX9+/dBCIG9vT309fXlPZXaKi4uhomJCYqKimBsrPjW45pGIpFAIBCgTZs2GtdwLzExEf369YOBgUG1ZRb4fD7CwsLg7+/PYnSNk5iYiPfffx+pqano2rWr3OcRi8Wwt7eHnZ0LOneOQUEBwZEjxujffwVcXBbg4EE/vHyZisTEu3B25inwGbArLy8PW7duxQ8//IDi4mKMGDERenr/ASHd8Pp1NBISglFUlCk7vl07PrZsUd/XTHFxMdq0aQMOh4Pz58/D3d2d7ZCoZkru91B5MqHdu3eTly9fyvNQjUNHTt4UExNDAJCzZ8+yHUqDREVFyWbpfP/99xo7C6U2d+/eVdj/k0OHmNk6lpa+ZNq0RGJt7U7s7AaRzp19CcAhffpEKWXGkToQiURk48ZNpGVLWwKAtG0r+H99ig+ZPl1IFi8WEW9vIXFwUP/XTFBQENHT0yMtW7Yk169fZzscqplS6ciJmZkZSktL4evriw8//BDDhg2DlpZWQ0+jEejIyZsIIXjvvfegq6uLhIQEcDgctkN6J+mIgLRV/cWLF5GUlIT27dtDIpHAz88PqampuHv3Lng8zRsRkL5ODxw4gPHjxzfqXJmZwNSp0UhKCkZxcaZse8uWdhgyZD2srPxRWAisXAnw+Y26lFrKzAS+/roCBQUHcOrUdEgkFeDzB2D48C0wN++KoiKgoECCnBw/PHigvq+ZM2fOYNCgQXBwcMCzZ8+QkJCAzp07sx0W1czI+x4qV/v6nJwcHDx4EDweDxMmTIClpSVmz56NxMREeU5HaRgOh4OVK1fi4sWLOHnyJNvh1It08b+vvvoKe/bsgbGxMcaPH4/y8nKFLP4nkTBvaikpzHdVr0NnYGAEXV09JCXlNfr6IhFgbu6Pzz67hylT4tCv31L06jUHn312B46O/jAwYFaPFokUFr5aEYmA8nJtmJnZQCKpwNChm8Dne2H79m6IjJwAIB9lZVxMnarYBSMVzcvLC23btkX//v3Rpk0bDBo0CFmKXDeBopRIruRES0sLPj4+2L9/P/Lz87Fp0yY8fPgQ3t7e6Nixo6JjBAA8f/4cQUFBMDExgYmJCYKCgt65fsbUqVPB4XCqfdF7r4oxbNgw9OnTB9988w3kGHxTuaqL/7Vu3RoHDx5ERkYG0v+/RHBjFv9LTwfWrgWWLgW+/Zb5vnat4lcfftv1163jQEvLHH/8kd/o60tXh371igc+3wve3iswYsRW8HjM6GhJCbPfyEiBT0KNSJ9/YSHzWnBzm45+/b6Gn99vePDgNLZvd8Tjx7+he3emtkcRC0YqA4/Hw/jx43HkyBEcP34cWlpaGDhwYLUVuylKXcmVnFSlr6+PoUOHYvjw4ejUqRMyMzMVENabAgMDkZycjOPHj+P48eNITk5GUFDQOx83bNgw5OTkyL6OHTumlPiaGw6Hg2+//RZ///03YmNj2Q7nnaou/gcAvXv3RkZGBrp3715te0MX/0tPB7ZsAZKSmAULu3RhviclMduVnaBUvb6hoTn09PIbff2qq0PXzDsJAR49gmx16KZI+vxLSpjXQn4+s2Bkt26TMHt2OszMhiExcQpmzBgCQDELRipLYGAgcnNzce/ePZw+fRqlpaUYMmQInj17xnZoFPVWcicnpaWl2L9/P0aMGAErKyts3LhRdt9e0dLT03H8+HHs2LEDffr0QZ8+ffDLL78gNjYWt2/ffutjdXV10bZtW9lX69atFR5fczVgwAD0798fS5culS1pr648PT3B5/OxZs0aWazS2WUSiQShoaGws7ODp6dnvc8pkQCHDwMFBYCTE6CnV4KCghQYGzM/FxQAMTHKu8VT8/rGxhYoLc1r9PWrrg6dlgYUFQGVlcz3tDRmu3R16KZI+vwdHDyhr89HXNwaVFRIUFQEPHxohpEj9+OHH/7AjRvXweFwkJSUBLFYzHbYterVqxc6duyI8PBwdOjQAadOnUJOTg6GDx8OUVO9L0c1CXL9eZk4cSLMzc3xxRdfwM7ODvHx8bh//z5WrVoFR0dHRccIoVAIExMT9O7dW7bN3d0dJiYm76xziY+Ph7m5OTp37oyPP/64yS0rzybp6ElycjIOHz7MdjhvxePxEBYWhtjYWPj5+UEoFEIkEkEoFMLPzw+xsbFYv359gwobs7KAW7cAGxuAwwFOnFiAffuGoKQkHxwOYG3NjFwo6zZ/zesbGJijpIR5fTf2+o6OwLx5gEAAFBYCd+4w393cmO1K+DVXK46OwPz5PPj5heHBg1j8+qsfbt0SomtXEfr1E+L48V/w+vVrDB06FAsWLEDfvn2RlpbGdthv4HA4mDhxIiIjI1FWVgYnJyecOHECt27dwqhRo/Dq1Su2Q6SoWsmVnHA4HBw8eBBPnjzBDz/8AA8PD0XHVU1ubi7Mzc3f2G5ubv7W+6fDhw/H/v37cfbsWYSFheHvv//GgAEDUFZWVudjysrKUFxcXO2LqpunpycGDx6MpUuXqu2nRyl/f39ERkYiJSUFHh4eMDY2hoeHB1JTUxEZGdngnhUiEVMYamDA/OzltRyESHD4cBDE4goUFMTjzp0InD0br5R/m5rX9/RcgtGjf5Xtb2zhqqMjsHgxMyvnm2+Y719+2fQTEylHR2DvXn9s2xYJDicFx497YMUKYwQGMq+ZqKgo/PXXX0hISMDz58/h6uqKFStWoLy8nO3QqwkMDERRURGOHz8OAHBzc8PRo0dx+fJljBs3Tu3ipSgAaFQTtsZavnw5VqxY8dZj/v77b5w8eRJ79ux54xZOp06dMH36dCxevLhe18vJyUH79u1x4MCBOt+I6oqJTiWu2+XLl+Hu7o79+/cjMDCQ7XDeSSwWIyEhATk5ObC0tISnp6dcU0EzM5niV1NTQPrSuH//FPbtGwo9PRO8fv1CdqwyGr3Vdv2qiorQpKf8qtK7XjOvX7/GqlWrsG7dOnTp0gU7d+6sNtLLNoFAgM6dO+PgwYOybSdPnoSvry/GjBmD/fv3q+V0aErzyTuVuN7NSbZs2YJPPvkEenp62LJly1uPnTdvXr3O+dlnn2HChAlvPYbP5+PGjRvIy8t7Y9/Tp09hYWFRr2sBTOFa+/btcffu3TqPCQkJwYIFC2Q/FxcXw8bGpt7XaI569+6NkSNHYvny5fjggw/UvucNj8eDl5dXo88jLZxMSmJqPDgcoLxcBIDg9esiCATbMWrUJAwZkoq1a9cgICBArhGahlxfSlq46ubWdAtXVeldrxk9PT2sWrUK48aNw/Tp09GnTx98/vnnWLVqFQykQ1ssCgwMxNKlSyESiWD0/2lWQ4YMwYEDBzBu3DgYGhril19+0YieRVQzUd9ubXw+nxQUFMj+u64vOzu7BnWBq4+0tLQ3VsO9dOkSAUBu3bpV7/MUFBQQXV1dsmfPnno/hnaIrZ9//vmnWa7lkZZGyKxZhAQEEDJvXiUxMeETO7uRpE2b94mhoQ1JTGR+Z8RiMfH19SV2dnYKXZel6vXnzyfkq6+Y7wEBzHZlrYJM1a2iooJ8//33pEWLFoTP59e6gruqPXz4kAAgv/322xv7fvvtNwKAzJ8/n0gkEhaio5oylS78x4Zhw4aRbt26EaFQSIRCIXFxcSE+Pj7VjunSpQuJjo4mhDBtqIODg0liYiLJyMggcXFxpE+fPqRdu3akuLi43telyUn9jRkzhtjZ2ZHy8nK2Q1Ep6QJxgwczi+YNGyYkixZlEROT1mTu3Lmy4xITlbNoXtUF6saNY76vWUMTE7bdu3ePeHt7EwBkypQpta56rEqenp5k+PDhte774YcfCACybNky1QZFNXlNPjkpLCwkkyZNIkZGRsTIyIhMmjSJPH/+vNoxVT+5l5aWkiFDhhAzMzOira1NbG1tyZQpU0hWVlaDrkuTk/q7ceMG4XA45KeffmI7FJUTiwnZvDmcACCpqSIiFhNy9erVamtQFRcXEwAkPDxcKdfPyCDkxg3me1Nc90YTSSQSsmPHDmJiYkLMzc3J77//ztroxLZt2wiPxyP5+fm17l+7di0BQMLCwlQcGdWUKX1tnap1GO+yYcOGeh+r7ujaOg0zYcIEJCYm4u7du9DV1VXptRVV6Cqv+Ph4eHt7QygU1tqJWCgUwsPDA3FxcQqpeaE0R05ODj777DNER0dj9OjR+OGHH9CuXTuVxlBQUABLS0ts3rwZs2fPrvWYJUuWIDQ0FD///DM+/vhjlcZHNU1KX5XYy8ur2peRkRHR19cnAoGACAQCYmBgQIyNjYm3t3cD8yr1RkdOGiY9PZ1wuVzy3//+V6XXjYqKInw+X7bqMADC5/NVumpsZWUl4fP5xNfXl4hrDF0oq+aE0iyRkZHEwsKCGBsbk59++umN14myjRgxgvTt27fO/RKJhMyZM4dwOBwSERGhwsiopkqlt3XCwsKIr68vefbsmWzbs2fPyOjRo8n69evlOaXaoslJwwUFBRFLS0tSWlqqkutFRUURDodDfH19iVAoJCKRiAiFQuLrq/pl7avGkpiYSIqLi0liYiIrsVDq6dmzZ+Sjjz4iAEj//v3JnTt3VHbtffv2EQAkMzOzzmPEYjGZPHky0dLSIn/88YfKYqOaJpUmJ1ZWViQ1NfWN7SkpKcTS0lKeU6otmpw03L179wiPx1PJvWvpaIWPjy+5f19creZCOlrB5/PJ6dOnSXh4OImLi1P6yEVtozh2dnY0MaGqOX36NOnQoQPR09Mja9euJRUVFUq/pkgkIi1atCBr165963EVFRXE39+f6OrqkjNnzig9LqrpUmlyYmhoWOsL9syZM8TQ0FCeU6otmpzIZ8aMGcTMzIyIRCKlXicujpkhM3OmkAQFEeLrm0dsbf3JwoUZJC2NkDVr1lRLElR1u6eyspLExcWpLCGiNNPLly/JggULCJfLJQKBgFy7dk3p1xw/fjzp3r37O497/fo1GTp0KDEwMCBCoVDpcVFNk0qTk6CgIGJra0sOHTpEsrOzSXZ2Njl06BDh8/lk8uTJ8pxSbdHkRD6ZmZlEW1ubhIaGKvU633/PzJDx8xORL74gxMNjMdHWNiSjRhWSIUOYWywAyIoVK1i93UNRb3PlyhXi4uJCeDweWbx4sVJviR45cuT/s8reHP2uqaSkhPTt25e0bNmSXL9+XWkxUU2XSpOTkpIS8umnnxJdXV3C5XIJl8slOjo65NNPP602dbIpoMmJ/GbPnk1at26ttH87sZiQ6dOZkZOPPhKSRYueER0dI+LhsYh8/XUl0dfnk/bt+xIAJCoqirx48eL/j1NccSodIaEUpaysjHz77bdER0eHdOrUiZw7d05p12nVqhX56quv6nX8ixcviJubGzE3Nye3b99WSkxU08VKn5OXL1+S69evk+Tk5CaXlEjR5ER+jx49Irq6umTFihVKOX9GBiGTJlUSY2M+6dzZl3h5rSBaWnokODiXTJnCJC1mZn2JjY0dmT17NnFwcJA9VhEN0dRhhhDV9KSlpREPD4//366cKUuqFWnGjBmkQ4cO9e658vTpU+Lk5ERsbGzeWkxLUTXJ+x4q16rEUgYGBujWrRu6d++uFutHUOqlXbt2+PTTT7FhwwY8f/5c4ecXiYDych6GDAnDnTuxyMiIg4fHImhr6yMz8xwA4OnTC/jii/W4evUqevXqJXuss7MzAKb/hDyio6MREBAAFxcXCIVCiEQiCIVCuLi4ICAgANHR0Y1/glSz5OjoiISEBPz3v//F/v370bVrV/zxxx8KvUZgYCAePHiAK1eu1Ot4U1NTnDp1Ctra2hg0aNBbV4OnKIVQUrLUZNCRk8bJzc0l+vr69R5CboiMDKZV+xdfEPLBB1GkZUv+G8WvAkEouXXrNdHR0anWe6UxIye0nwmlKg8fPiTDhw8nAMj48eNJXl6eQs5bWVlJrKysyLx58xr0uAcPHhArKyvi7OzMejt+SjOwMnJCUe9iYWGBzz77DJs3b0ZBQYFCzy1dlTc7G3Bw8MfcufcwZUoc/P3D8eGHp6Gvz8erV4l49iwZ5eXleO+99wAAEokEoaGhsLOzg6enZ4Ovm5CQgMzMTCxZsgRcLherV6/GsmXLAABcLhchISHIyMhAQkKCQp8v1fzY2tri6NGj2LdvH06fPg1HR0f89ttvIPVr7F0nHo+H8ePH4+DBg6isrKz34+zs7HD69Gnk5uZi+PDhEIlEjYqDoupCkxNK6f7zn/8AAL777juFnpfLBcaMAUxNgbQ0QCTiwdraC7a2E/H69UD07RuG27dj8cknM6ClpYUOHTpAKBTCz88PsbGxWL9+vVzt7aW3gqS3hv78809kZGTI9jf2lhFFVcXhcDBp0iSkp6dj2LBhmDJlCoYPH47MzMxGnTcwMBB5eXmIj49v0OMcHR1x4sQJ3Lp1C6NGjcKrV68aFQdF1YYmJ5TSmZqaYv78+di6davC71U7OgLz5gECAVBYCNy5w3x3cwM2bPDHjz9G4v79DFRWVsLMzAweHh5ITU1FZGQk/P395bqmpaUlACA1NRUSiQQpKSno1q2bbH9qamq14yhKEczMzLB//37Exsbi5s2bcHZ2xubNmyEWi2XHSCRAZiaQksJ8l0jqPl+PHj3QqVMnhIeHNzgWNzc3HDt2DFeuXMG4ceNQXl7e8CdEUW9R74X/miu68J9iPH/+HHZ2dpg6dSo2bdqk8PNLJEBWFlMka2QElJQAR44At24BMTGOMDHpAg+P8RgxwhIffti4BQHFYjHs7e3h4uKC9evXo0uXLjhx4gSGDBkCiUQCPz8/pKam4u7duypdeJBqPoqLixESEoIff/wR7u7u2LFjB7jcrjh8mHnNv34N6Okxtz3HjGGS+NosX74cGzduRF5eHvT09Bocx6lTp+Dj4wM/Pz+Eh4fT1zv1BnnfQ+nICaUSrVq1QnBwMLZv345Hjx4p/PxcLsDnAy4uwKtXwNatQFISYGRUBJHoFpycxkAimYhLl7xw507j/oDyeDyEhYUhNjYWkyZNAgCF3TKiqPowNjbGDz/8gISEBDx//hyurgJMnrwC//xTDlNToEsX5nZnUhKwZQuQnl77eSZOnIji4mL89ddfcsUxePBgHDhwAFFRUZg5c2aja2EoSoomJ5TKfP755zAwMMCaNWuUdg2JBDh8GCgoAJycAJHoKgDA3v49ODkx22Ni3j7cXR/+/v6IjIzEvXv3AACdOnVSyC0jimqIvn374tq1ZPTtuwj//LMKCQluKC6+DB4PMDbGO1/zXbp0gZubm1y3dqTGjBmDX3/9FTt37sSCBQtogkIpBE1OKJUxNjbGokWLsGPHDjx8+FAp18jKYoa1bWwADgd4/PgKdHSMYGraBRwOYG3NfIrMymr8tfz9/eHt7S374x4XF4e7d+/SxIRSqfx8PdjYrMKkSVehrd0CO3f2wfHj81FeXlKv13xgYCD+/PNPFBcXyx1DUFAQfvjhB2zatAkrVqyQ+zwUJUWTE0qlPvvsM7Rq1QrffvutUs4vEjH326U9ASsqSmBn5w0Oh3mpGxgw+xU1AzIlJQX9+vXDxIkT4eXlRW/lUConfc3z+d0xfboQgwd/jzt3YiGRVAB492t+/PjxKC8vx+HDhxsVx+zZs7F27VqsWLECYWFhjToXRdHkhFIpAwMDLF68GLt375bdElEkIyOmELCkhPl5wIBVmDDhiGx/SQmz38io8dcqKSnB/fv3q83UoShVq/qa53K14OERjDlz0qGn1xLAu1/z1tbW6NevHyIiIhody5dffoklS5Zg4cKF+OWXXxp9Pqr5oskJpXKzZs2ChYUFVq5cqfBzV23MVvPWNyHAo0fMzAVb28Zf6+bNmyCE0OSEYlVtr3keTxtA/V/zgYGBOH36NPLy8hodz6pVq/DZZ59h5syZOHDgQKPPRzVPNDmhVK5FixZYsmQJ9u/fj/S6phHIqWZjtqIioLKS+Z6Wxmz382OOa6wbN26Ay+XCycmp8SejKDkp4jU/duxYcDgcHDp0qNHxcDgcbN68GZMnT0ZQUBD+/PPPRp+Tan5on5N3oH1OlKOsrAydOnVCnz59cPDgQYWfPz0db/R8cHRk/kjX1fOhLmKxGAkJCcjJyYGlpSU8PZk+KfPmzcOpU6cUnmBRlDwa+5r39fXFs2fPcPHiRYXEU1lZifHjx+Po0aM4duwYBgwYoJDzUppF3vdQmpy8A01OlGfHjh34+OOPcf36daXcGqnZmM3WtuEjJtHR0QgODq7WKpzP5+PLL8Owc+cWmJub488/f1fISAxFNVZjXvMRERGy1Yrt7OwUEk9ZWRlGjx6NCxcu4PTp03B3d1fIeSnNQZuwURpnypQp6NChg2zRPEWr2piNz5cvMQkICICLiwuEQiFEIhEiIoTQ03PBp5+ORVLSP3j8uBvWrq27yRVFqVJjXvOjRo2Cvr6+QutEdHV1ER0dDTc3NwwfPhzXr19X2Lmppo0mJxRrtLW1sWzZMsTExOCff/5hO5xqxGIxgoOD4ePjg5iYGLi7uyM72xDnzrmja9cY2NoOglj8Era2zu/swinv9ePj4xEREYH4+Phq66dQlDIYGBhg9OjRjWrIVht9fX38+eef6NixI4YMGYI7d+4o9PxU00STE4pVkyZNQpcuXbB06VK2Q6kmISEBmZmZWLJkCbhcLjIyHmLp0t+Qm1uKrl256NJlOABAR+e1QjvPAsyIjb29Pby9vREYGAhvb2/Y29sjOjq68SenqLcIDAxEamoqUlJSFHpeExMTHD9+HKamphg0aJDSmjBSTQdNTihW8Xg8LF++HMeOHYNQKGQ7HJmcnBwAgJOTMzIzgY0bTyAq6iO0a8d0nq2oeAkAkEgkCu08K72V5OzsguhoIS5dEiE6WghnZxcEBATQBIVSqiFDhqB169YK6XlSk6mpKU6dOgVtbW0MGjRItkJ5Q1ZSppoPWhD7DrQgVvkkEgm6d+8OCwsLfP3112/MimFDfHw8vL29MXOmEKWl7jh7dgby86/i/feT4eQE5OcfQmTkB5gyJQ58vhcqK4E7d4BvvmHu98tDutqxnZ0L+vePRHJyNjictjA0NEDnzhKcPeuHzEy62jGlXDNnzsTJkyfx4MEDcDgchZ8/IyMDffv2RevWrfHzz+cQF9e6QSspU5qFFsRSGovL5WL48OE4c+aM2tzKMDX1hJERH3/8sQZt2kjw+vXfMDDohUePgEuXJLh6dS9atrSDra0nAMV0npXeSmrZcgkuX36CmBh7tGghhKkpcP06FyYmIcjIyEBCQoKCniVFvSkwMBCZmZm4dOmSUs5vZ2eH06dP4/HjXPj5DceVK6IGraRMNQ80OaFYFx0dje+//x7GxsZwdXVFcXExhEIhXFzYuZUhkQB//MGDs3MYcnJiERvrg2fPbsLMzAUVFULcueOHzMxYDB68HlwuT2GdZx8/Zm4lEeKM7t1NAQBlZQWy1WUJca52HEUpg6enJ9q1a6fwwtiqunRxRGDgSTx/fhvJyaPQosWreq+kTDUPNDmhWCWdFePr64u9e/ciOTkZly9fhru7O2JiYuDj44OFCxeqdLaKdGVjd3d/fPBBJHJzk0CIGHfvfo6HDz3w+nUq2rWLhLm5v0I7z/J4lgAAPb1U6OgYgMfTRUnJUwBMnYuubmq14yhKGbhcLiZMmIDff/8dlZWVSrlGVhZQXCzAmDFH8eTJFRw6NA5icTkAKHz1cEoz0eSEYlXVWTG+vr7o3bs3vvnmGxBCwOVyERKi+lsZVVc2dnT0h6fnV+DxdODn9xtGj47D++/fBY/nj7t3gcJCwM0NmDev8ffIHRw8oa/Pxz//rAFAoK/fBqWlTHJCiARJSaEwMLCDg4Nn458kRb1FYGAg8vPzcfbsWaWcX/o71rnz+xg/PgYPHpzG/fsnZfsVvXo4pXm02A6Aat6ks2KcnZ3B4XCwYcMGPHv2TLbf2dm52nGqUHWVV2Nj4L33PkPXruNhYGAGAHjxAmjdGvjkE6Z4T57Os7Vp2ZKHXr3CcO5cAA4c8MOoUb+iXbv3kJ0txIULobh7Nxb9+0eiZUtaDEspl0AgQJcuXRAeHo4hQ4Yo/PxVf8c6dhyMuXPvwMTk33uiilw9nNJMGjNysnr1anh4eEBfXx8tW7as12MIIVi+fDmsrKzQokULeHl54ebNm8oNlGoQS0vmFkVqKnPLwsPDAz4+PrJZAtLt0uNUobZVXqWJCSHA48dAr17AkCHydZ5923WHDPFHnz6RyM9Pwf79Q/Hdd/9r797DoizzPoB/h1GH04DmyAAJzCTJIbBES6EQ8JSWLkhwedj10vQts7RQ3Hyx1WzbFV2lNN0yd683t6uVLAQqrFXb0GjFAoMURVBhHFRIQeUUB5l53j9mmRg5CCo8M8z3c11z2XPPc/jNxDA/7ud33/cQ/N//heDKlQIEB6fgySej78mKykRdkUgkmDdvHlJTU9HQ0HDPz3/rZ6xtYnKvVw8ny2QxyUlzczNiY2OxdOnSbh/zl7/8BW+99RZ27NiBnJwcuLq6YsqUKahlX6HZCA0NhUqlwoYNG6C/pfpNr9cjMTERarUaoaF9dyujL1c27ui6Dz8cjUmTzuGZZzIRFbUHzzyTiUmTzuLhh6N75bpEHZk7dy5qa2uxf//+e35usT5jZDksbp6T3bt3Iy4uDjdu3OhyP0EQ4O7ujri4OKxevRqAYREqpVKJTZs2YcmSJd26Huc5+VVnq/PerdaJx2bMmIGEhAQEBASgoKAAiYmJyMjIQEpKCqKjo+/BK+iZe7mysSVcl+hWjz76KDw9PbFv375eOT9/1vu/O/0O7bc1J6WlpaioqDC5XyqTyRAWFoajR492Ozkhg45W5/Xy8sKiRYvw4IMPmiQrPV0ZNTo6GikpKYiPj0dISIixXa1Wi5aYAIZfjj4+d7+ysaVcl+hW8+bNQ0JCAm7cuNHt2+k9wZ916ky/TU5ap0ZWKpUm7Uqlsst1HZqamtDU1GTcrqmp6Z0ALUjbno3k5GQEBARgx44d+OMf/2iyorBKpcKKFUmoq4vu8YyP0dHRiIyM7JWembvRusqrtVyXqK3Zs2cjPj4eaWlpePbZZ3vlGvxZp46Imp+uX78eEomky0dubu5dXePW6ZcFQehySubExEQ4OzsbHx4eHnd1fUt36+q8jz32GA4cOIA1a9Zg0qRJeOKJJ+Dl5YXvvvsOanUgXnklBhkZqXc046NUKkV4eDjmzp2L8PBw0RMTImvn7u6O8PDwXllrh6grovacLFu2DHPmzOlyH9UdptSurq4ADD0obUd6XLlypV1vSlsJCQlYuXKlcbumpsaqE5TWeUiSk5NhY2OD7OxsxMbGIiwsDJ999hm+//57hISEoKnpJiZOTMeZM1EoLFyFiIjpsLGxhZOTBP7+hiK39HRDwsIuWyLLMW/ePCxZsgQVFRXG36tEvU3UrwmFQgFfX98uH7a2tnd0brVaDVdXVxw6dMjY1tzcjCNHjpjUNdxKJpPBycnJ5GHN2s5DAgBffPEFBEHAH/7wB9jY2MDmv5lGQUE5iottMGFCAm7cKMX27Q+iqOhzAJzxkciSPfPMM5BKpfjkk0/EDoWsiMX8DavVapGfnw+tVgudTof8/Hzk5+ejrq7OuI+vry/S0tIAGG7nxMXFYcOGDUhLS0NBQQEWLlwIe3t7zJs3T6yXYXFunYek8L/3Zlxdx+HkSWDp0jgAgKOjGxobAU9PQxJjb6/AgQMr0NLSCIAzPhJZqiFDhmDatOnYvTsZJ08CGg3XvKHeZzEFsevWrcM//vEP4/bo0aMBAJmZmQgPDwcAFBUVobq62rjPq6++ioaGBrz44ou4fv06xo0bh4MHD0LOaQe7re08JOnp6bh48ToAYNWqAgjCYOTlHYO9vQL33x8KW1tAqzUkMcHBK/H554tx9OgWeHo+gZ9/LodO5wZ7+1AArCUhshSFhYC9/Tzk5c3B6tUlUCge6FaRO9HdsLh5Tvoa5zn5dbROePgM5OaeQkPDFbi4PIL6+ouortZg3LhkPPLIHMhkenz6aRQaGgqwfPlZ7N07C8XFGQB+/RFTqVRISkoSbXgwEXVfYaGhmL2i4hfs3++Cxx9PwOjRr6GszFDsfi/WlKL+7U6/Qy3mtg6JJzo6Gp98koK8vBOorS1BS0sdLl/+DtXVGnh5hSM8/GkUF2cjIyMK5eUZ8PPbgry8z1Bc/AWk0oFQKqdi8eJaJCdnIzAwEDExMUhNTRX7ZRFRF/R6wwRplZXAqFH28PWNwqlTeyCXC/D3N7Snp/MWD/UOJifULWPHRiMs7GsAQEjIajz44NMAgAsXDmPTJidkZhrWf1m3LgVPPRWJzMx4uLnNxKOP/hU//3wQERE/Yc6c8UhPT8eMGTOwatUq6HQ6MV8SEXVBqzXM3OrhYShqDwych6tXT+PKlZMscqdeZzE1JySu2lrg2rVzAICxY5fA0dEVWu13kEoHora2HPb2bmhsDEVMjBRXrx5Gfb0GO3cmIyTkMcye/T6SkpZjzpwcSKVSJCQkICQkBFlZWcZ6ISIyL7W1hiJ2BwfD9gMPTMH8+YcwbJg/AEP7pUsscqfewZ4T6ha5HGhoKIJUKoOzsycGDrTDiBFToFKFIzBwLhSKcNjZSSGXAz//bBh+HBUVgAcesMH27dshk8lQWVkJ4Ndhya3DlInI/Mjlhhme6+sN21LpQDzwwGTY2Bj+pq2vNzzP8QXUG5icULcY1rsogoPDg5BITEfb3LrE+a3Dj8ePH4+jR48aJ79rbW87OR4RmRdPT8PSE2Vlhs94W7d+5onuNSYn1C02NoBEUgSFwve2S5y3HX6s/2+1XOuSAXq9HomJiVCr1QgNDRXxFRFRV2xsDMOFFQrc9jNPdK/xx4q67eLFIkyZ4oPRo4GqKqC42PBvUJDpkEKpVIqkpCRkZGQgKioK2dnZqK2tRXZ2NqKiopCRkYEtW7Zw7RwiM+fnZ/hs3+4zT3SvsSCWuqWurg6XLl3C44/74Le/vf0S59HR0UhJSUF8fLzJcgFqtRopKSmc54TIQvj5GdbEut1nnuheYnJC3VJcXAwA8PHx6fYS59HR0YiMjERWVhbKy8vh5uaG0NBQ9pgQWZjufuaJ7hUmJ9QtRUVFAAzJSU9IpVIOFyYioh5hxxx1y5kzZ6BUKuHs7Cx2KERE1M8xOaEu6fWGVUh/+KEIXl4+nKqaiIh6HW/rUKcKCw1ra5w5Axw7VgSF4lFs3MjVSImIqHex54Q61LoaaV4eMHSoHvX1xXB390VenqG9sFDsCImIqL9ickLttF2N1N8fAC7h5s1fcP/9PlyNlIiIeh2TE2rn1tVIq6oMI3UUCh+uRkpERL2OyQm1c+tqpEOGjMCUKZsxeLAKgKG9sZGrkRIRUe9gckLt3Loa6ZAhaoSErOJqpERE1CeYnFA7XI2UiIjExKHE1E7raqRlZYbVR4cPN9zKqa83JCaWtBqpTqfj9PlERBbGAr5eSAz9YTXS1NRUeHt7IyIiAvPmzUNERAS8vb2RmpoqdmhERNQF9pxQpyx5NdLU1FTExMRgxowZSE5ORkBAAAoKCrBhwwbExMRwZWQiIjMmEYRbqwqorZqaGjg7O6O6uhpOTk5ih0PdoNPp4O3tjcDAQKSnp8OmTTal1+sRFRWFgoICnD17lrd4iIh60Z1+h1rA38BEPZOVlQWNRoM1a9bAxsYGf/vb35CXlwcAsLGxQUJCAkpLS5GVlSVypERE1BEmJ9TvlJeXAwACAgIgCAJefvllHDlyxPh8QECAyX5ERGRemJxQv+Pm5gYAKCgowM8//4zGxkao1Wrj8wUFBSb7ERGReWFyQv1OaGgoVCoVNmzYgJKSEgCASqUCYKg5SUxMhFqtRmhoqIhREhFRZzhah/odqVSKpKQkxMTE4NKlSwCAoUOHIjs7G4mJicjIyEBKSgqLYYmIzBSTE+qXoqOjkZKSgsWLFwMAPDw8AABqtZrDiImIzByTE+q3oqOj8dVXX+HIkSN44403OEMsEZGFYHJC/ZpGo8GoUaMwd+5csUMhIqJuYkEs9WsajcZYDEtERJbBYpKTP//5zwgJCYG9vT0GDx7crWMWLlwIiURi8hg/fnzvBkpmQ6fT4cKFCybDiImIyPxZTHLS3NyM2NhYLF26tEfHTZs2DeXl5cbHl19+2UsRkrkpLy/HzZs32XNCRGRhLKbm5I033gAA7N69u0fHyWQyuLq69kJEZO5KS0sBgD0nREQWxmJ6Tu7U4cOH4eLigpEjR+K5557DlStXxA6J+ohGowEAeHl5iRsIERH1iMX0nNyJ6dOnIzY2Fl5eXigtLcXatWsxceJEHD9+HDKZrMNjmpqa0NTUZNyuqanpq3DpHistLYWLiwscHBzEDoWIiHpA1J6T9evXtytYvfWRm5t7x+efPXs2nn76aQQEBGDmzJn46quvUFxcjP3793d6TGJiIpydnY2P1sm7yPJoNBre0iEiskCi9pwsW7YMc+bM6XKfe1nM6ObmBi8vL5w9e7bTfRISErBy5Urjdk1NDRMUC1VaWspiWCIiCyRqcqJQKKBQKPrselVVVSgrK+tyNVqZTNbpLR+yLBqNhkPHiYgskMUUxGq1WuTn50Or1UKn0yE/Px/5+fmoq6sz7uPr64u0tDQAQF1dHVatWoXs7GxoNBocPnwYM2fOhEKhwKxZs8R6GdRHWlpaUFZWxp4TIiILZDEFsevWrcM//vEP4/bo0aMBAJmZmQgPDwcAFBUVobq6GoBhZdqTJ0/iww8/xI0bN+Dm5oaIiAjs3bsXcrm8z+OnvnXx4kXodDrWnBARWSCJIAiC2EGYs5qaGjg7O6O6uhpOTk5ih0PdlJmZiYkTJ6KoqAgjR44UOxwiIqt0p9+hFnNbh6gnOMcJEZHlYnJC/Y5eD+TllcLFxR3l5TLo9WJHREREPWExNSdEbel0OmRlZaG8vBxubm4IDQ2FVCpFYSGQlgbs36+BTqfGunWAry8waxbg5yd21ERE1B1MTsjipKamIj4+3njrBjDMh7NiRRIKC6NRWQk0NpbCzU0NhQLIywPKyoCXX2aCQkRkCXhbhyxKamoqYmJiEBgYiOzsbNTW1iI7OxsBAYGIi4vBTz+lwt8fqKvTYMgQFZycAH9/oLISSE8Hb/EQEVkAJidkMXQ6HeLj4zFjxgykp6dj/PjxcHR0xPjx47FtWzruv38GTp1ahZs3f0FNzSUMGWIYRiyRAMOHA4WFgFYr8osgIqLbYnJCFiMrKwsajQb/+79roNXa4KOPsrFixRtobGxGfb0NRo5MQE1NKS5cOIzHH1+N++9/zHisgwPQ2AjU1or4AoiIqFtYc0IWo7y8HABw4EAASkuB7OxklJd/DoViHcaNA1xcAgAAjY3VmDw50eTY+nrA1hbg/HtEROaPPSdkMZqbDWsiZWUVYOhQAZWV+zFixFPIz5cgJQXQ6wsAAI6OpmsnCQJw8aKhGNbTs8/DJiKiHmJyQhZBrwfKykJhb6/CpUsbcPNmEW7cKMFDDz3134JXPXJzEyGXq1FfH4rqaqClBaiuBk6fBhQKICoKsOFPPBGR2eNtHbIIWi1QXCzFxIlJyMiIwfXrJbCxGQg3t7G4eDEbeXmJKCnJwLp1KZDJpDhzBrh0yXArJyjIkJj09TDizuZiISKirjE5IYtQW2soaH3kkWjY26fgX/+Kg15/E2+9ZbiFM3iwGsHBKYiJicZDDxmSmdpaQ42Jp2ff95h0NhdLUlISoqOj+zYYIiILw+SELIJcbugFqa8H/Pyi4eMTCa02C7W15ZDL3eDsHIrr16WQyw2JiEolXqytc7HMmDEDycnJCAgIQEFBATZs2ICYmBikpKQwQSEi6gJXJb4NrkpsHvR6YONGw2yv/v6GuUtaCYKhriQoCFi9Wty6Ep1OB29vbwQGBiI9PR02bYLR6/WIiopCQUEBzp49y1s8RNTvcVVi6tdsbAzr4ygUhkTEXAteW+diWbNmDRobG7F161Zcv379v6/BBgkJCSgtLUVWVpa4gRIRmTHe1rFillaw6ednWB8nLQ1mUfDakda5WL7++mtERkaiqqoKKpUKUVFRAICAgACT/YiIqD0mJ1bKUgs2/fwAHx/xC147Ultbi8zMTADA+vXrsWjRIiQkJECtVhv3KSgwzMXi5ubW4TmIiIi3daxSZ4vnBQYGIiYmBqmpqWKH2KXWgtfAQMO/Yicm1dXV+NOf/gSVSoXdu3dDLpcjPDwcO3fuNElM9Ho9EhMToVarERoaKmLERETmjQWxt9HfCmJvLdiUSCSQ/Le6lAWbPXP9+nW888472Lp1KxoaGvDcc8/h1VdfRU5OjnG0TkJCgnG0TmJiIjIyMjhah4isBgtiqVvaFmxKJBJ4eXnhgw8+AMCCze6qqqrC2rVroVKpsHHjRixYsAAlJSXYvn07PDw8EB0djZSUFJw8eRIhISFwcnJCSEgICgoKmJgQEXUDa06sTGshZkBAAC5evIiysjIMGzbM+Lw1F2zq9V3XslRWVuKtt97C9u3bodfrsXTpUqxatQqurq7tzhUdHY3IyEiLKjgmIjIXTE6sTGshZkFBASorKwEAo0aNMj5vrQWbhYW/jgJqbDSMAvL1NQxfHjr0CrZs2YJ3330XALBs2TKsXLkSLi4uXZ5TKpUiPDy8D6InIupfWHNyG/255mTcuHHYsmULrl27BolEYrU1J4WFwDvvAJWVgIcH4OBgmIn27NkKaLWbUVz8HgYOHIDly5djxYoVUCgUYodMRGQR7vQ7lD0nVkYqlSIpKQkxMTE4fvw4VCoV6urq2hVsWktiotcbekwqK3+deba29jL+859N+PHHXZBIZAgN/T327n0FCsV9YodLRGQVmJxYodaCzTlz5uDy5cvGbFatVltdwaZWa7iV4+FhSEyKir7Ap5/GYuBAOzzxRAJ8fV9Gff1g1NUZZqElIqLex+TESk2fPh06nQ7x8fEYM2aM1RZstq527OBg2PbwCEFY2Ot47LGXIJM5oaUFqKoy7EdERH2DyYmVOn36NPR6PWJjYzFu3DixwxFN29WOnZwAe/uhCA1NMD5fX294Xi4XMUgiIivDeU6s1IkTJyCRSIxDh62Vp6dhVE5ZmWF147YEAbh40TBlvqenOPEREVkjJidW6sSJE/D29oZD6/0MK2Upqx0TEVkT/sq1Mno9oNEAR4+ewIgRo6DXix2R+FpXOx492lBfUlxs+DcoyNBuDqsdExFZE9acWJHWicYKCwXk5/8EP7+XsXGjoefA2r+AzXm1YyIia8PkxEq0nWhsyJByNDdXwdNzFPLyDPUW7CH4dbVjIiISF/8utAK3TjT2yy8nAABq9cPw9ze0p6eDt3iIiMgsWERyotFosHjxYqjVatjZ2WHEiBF4/fXX0dzc3OVxgiBg/fr1cHd3h52dHcLDw3Hq1Kk+itp83DrR2M8/n8CgQXIMHuwFiQQYPtzQs6LVih0pERGRhSQnZ86cgV6vx/vvv49Tp07h7bffxs6dO7FmzZouj/vLX/6Ct956Czt27EBOTg5cXV0xZcoU1FrZjFq3TjTW0HAd7u5jIJEY/vc7OBiet7K3hYiIzJTFLvy3efNmvPfeeygpKenweUEQ4O7ujri4OKxevRoA0NTUBKVSiU2bNmHJkiXduk5/WPhPowHWrTMMi219CYIgQCKRADAMm62qAv74R9ZcEBHRvXOn36EW0XPSkerqatx3X+cLsZWWlqKiogJTp041tslkMoSFheHo0aOdHtfU1ISamhqTh6XraKKx1sSEE40REZG5scjk5Pz589i+fTteeOGFTvepqKgAACiVSpN2pVJpfK4jiYmJcHZ2Nj48PDzuTdAi4kRjRERkSUT9Olq/fj0kEkmXj9zcXJNjLl++jGnTpiE2Nhb/8z//c9trtPYQtGp7O6MjCQkJqK6uNj7Kysru7MWZGU40RkRElkLUeU6WLVuGOXPmdLmPqk0RxOXLlxEREYHg4GDs2rWry+NcXV0BGHpQ3NzcjO1Xrlxp15vSlkwmg0wm60b0locTjRERkSUQNTlRKBRQKBTd2vfSpUuIiIjAmDFj8MEHH8DmNt+oarUarq6uOHToEEaPHg0AaG5uxpEjR7Bp06a7jt1ScaIxIiIydxbxN/Ply5cRHh4ODw8PbNmyBVevXkVFRUW72hFfX1+kpaUBMNzOiYuLw4YNG5CWloaCggIsXLgQ9vb2mDdvnhgvg4iIiLrBIqavP3jwIM6dO4dz585h+PDhJs+1HQldVFSE6upq4/arr76KhoYGvPjii7h+/TrGjRuHgwcPQi6X91nsRERE1DMWO89JX+kP85wQERGJwermOSEiIqL+ickJERERmRUmJ0RERGRWmJwQERGRWWFyQkRERGaFyQkRERGZFSYnREREZFYsYhI2MbVOA1NTUyNyJERERJal9buzp1OqMTm5jdraWgCAh4eHyJEQERFZptraWjg7O3d7f84Qext6vR6XL1+GXC6HRCIRO5x+q6amBh4eHigrK+NMvH2I73vf43ve9/iei6P1fT99+jR8fHxuu2BvW+w5uQ0bG5t26/lQ73FycuIvDxHwfe97fM/7Ht9zcdx///09SkwAFsQSERGRmWFyQkRERGaFyQmZBZlMhtdffx0ymUzsUKwK3/e+x/e87/E9F8fdvO8siCUiIiKzwp4TIiIiMitMToiIiMisMDkhIiIis8LkhIiIiMwKkxMyS3/+858REhICe3t7DB48WOxw+qV3330XarUatra2GDNmDLKyssQOqV/79ttvMXPmTLi7u0MikSA9PV3skPq9xMREPProo5DL5XBxcUFUVBSKiorEDqtfe++99zBq1CjjhHfBwcH46quvenweJidklpqbmxEbG4ulS5eKHUq/tHfvXsTFxeG1115DXl4eQkNDMX36dGi1WrFD67fq6+vx8MMPY8eOHWKHYjWOHDmCl156CceOHcOhQ4fQ0tKCqVOnor6+XuzQ+q3hw4dj48aNyM3NRW5uLiZOnIjIyEicOnWqR+fhUGIya7t370ZcXBxu3Lghdij9yrhx4xAUFIT33nvP2Obn54eoqCgkJiaKGJl1kEgkSEtLQ1RUlNihWJWrV6/CxcUFR44cwYQJE8QOx2rcd9992Lx5MxYvXtztY9hzQmRlmpubcfz4cUydOtWkferUqTh69KhIURH1vurqagCGL0vqfTqdDh9//DHq6+sRHBzco2O58B+RlamsrIROp4NSqTRpVyqVqKioECkqot4lCAJWrlyJJ554AgEBAWKH06+dPHkSwcHBaGxshKOjI9LS0uDv79+jc7DnhPrM+vXrIZFIunzk5uaKHabVkEgkJtuCILRrI+ovli1bhhMnTiA5OVnsUPo9Hx8f5Ofn49ixY1i6dCkWLFiA06dP9+gc7DmhPrNs2TLMmTOny31UKlXfBGPFFAoFpFJpu16SK1eutOtNIeoPli9fjs8//xzffvsthg8fLnY4/d6gQYPg7e0NABg7dixycnKwbds2vP/++90+B5MT6jMKhQIKhULsMKzeoEGDMGbMGBw6dAizZs0yth86dAiRkZEiRkZ0bwmCgOXLlyMtLQ2HDx+GWq0WOySrJAgCmpqaenQMkxMyS1qtFteuXYNWq4VOp0N+fj4AwNvbG46OjuIG1w+sXLkS8+fPx9ixYxEcHIxdu3ZBq9XihRdeEDu0fquurg7nzp0zbpeWliI/Px/33XcfPD09RYys/3rppZewZ88efPbZZ5DL5cbeQmdnZ9jZ2YkcXf+0Zs0aTJ8+HR4eHqitrcXHH3+Mw4cP41//+lfPTiQQmaEFCxYIANo9MjMzxQ6t3/jrX/8qeHl5CYMGDRKCgoKEI0eOiB1Sv5aZmdnhz/SCBQvEDq3f6uj9BiB88MEHYofWby1atMj4e2XYsGHCpEmThIMHD/b4PJznhIiIiMwKR+sQERGRWWFyQkRERGaFyQkRERGZFSYnREREZFaYnBAREZFZYXJCREREZoXJCREREZkVJidEdFd2796NwYMHix2GWZk/fz42bNhwT895+PBhSCQS3LhxAwCQkZGB0aNHQ6/X39PrEJkDJidEdFdmz56N4uJiscO4K+Hh4YiLi7sn5zpx4gT279+P5cuX9+q1Z8yYAYlEgj179vT4WCJzx+SEiG6rubm5w/abN2/Czs4OLi4ufRxR99y8ebPPr7ljxw7ExsZCLpf3+rWeffZZbN++vdevQ9TXmJwQWbDw8HAsX74ccXFxGDJkCJRKJXbt2oX6+no8++yzkMvlGDFiBL766ivjMTqdDosXL4ZarYadnR18fHywbds2k/MuXLgQUVFRSExMhLu7O0aOHAmNRgOJRIJPPvkE4eHhsLW1xUcffdTuts758+cRGRkJpVIJR0dHPProo/j6669Nzl9eXo6nn34adnZ2UKvV2LNnD1QqFbZu3Wrcp7q6Gs8//zxcXFzg5OSEiRMn4qeffur0vegsvqqqKsydOxfDhw+Hvb09AgMDkZycbPJajxw5gm3btkEikUAikUCj0QAATp8+jaeeegqOjo5QKpWYP38+KisrO41Br9fj008/xW9+8xuT9nfffRcPPvggbG1toVQqERMTc9trf/nllxg5ciTs7OwQERFhbG/rN7/5DX744QeUlJR0GhORRbrnq/4QUZ8JCwsT5HK58OabbwrFxcXCm2++KdjY2AjTp08Xdu3aJRQXFwtLly4Vhg4dKtTX1wuCIAjNzc3CunXrhB9++EEoKSkRPvroI8He3l7Yu3ev8bwLFiwQHB0dhfnz5wsFBQXCyZMnhdLSUgGAoFKphH379gklJSXCpUuXhA8++EBwdnY2Hpufny/s3LlTOHHihFBcXCy89tprgq2trXDhwgXjPpMnTxYeeeQR4dixY8Lx48eFsLAwwc7OTnj77bcFQRAEvV4vPP7448LMmTOFnJwcobi4WIiPjxeGDh0qVFVVdfhedBbfxYsXhc2bNwt5eXnC+fPnhXfeeUeQSqXCsWPHBEEQhBs3bgjBwcHCc889J5SXlwvl5eVCS0uLcPnyZUGhUAgJCQlCYWGh8OOPPwpTpkwRIiIiOv3/kZeXJwAQKioqjG05OTmCVCoV9uzZI2g0GuHHH38Utm3b1uW1tVqtIJPJhFdeeUU4c+aM8NFHHwlKpVIAIFy/ft3kmi4uLsLu3btv/8NCZEGYnBBZsLCwMOGJJ54wbre0tAgODg7C/PnzjW3l5eUCACE7O7vT87z44ovCM888Y9xesGCBoFQqhaamJmNb65f/1q1bTY69NTnpiL+/v7B9+3ZBEAShsLBQACDk5OQYnz979qwAwJic/Pvf/xacnJyExsZGk/OMGDFCeP/99zu8RmfxdeSpp54S4uPjjdthYWHCK6+8YrLP2rVrhalTp5q0lZWVCQCEoqKiDs+blpYmSKVSQa/XG9v27dsnODk5CTU1NR0e09G1ExISBD8/P5PzrF69usPkZPTo0cL69es7e6lEFmmASB02RHSPjBo1yvjfUqkUQ4cORWBgoLFNqVQCAK5cuWJs27lzJ/7+97/jwoULaGhoQHNzMx555BGT8wYGBmLQoEHtrjd27Ngu46mvr8cbb7yBjIwMXL58GS0tLWhoaIBWqwUAFBUVYcCAAQgKCjIe4+3tjSFDhhi3jx8/jrq6OgwdOtTk3A0NDTh//nyX1781Pp1Oh40bN2Lv3r24dOkSmpqa0NTUBAcHhy7Pc/z4cWRmZsLR0bHdc+fPn8fIkSPbtTc0NEAmk0EikRjbpkyZAi8vLzzwwAOYNm0apk2bhlmzZsHe3r7TaxcWFmL8+PEm5wkODu5wXzs7O/zyyy9dvhYiS8PkhMjCDRw40GRbIpGYtLV+wbUOOf3kk0+wYsUKJCUlITg4GHK5HJs3b8b3339vcp7Ovrxv96X++9//HgcOHMCWLVvg7e0NOzs7xMTEGItqBUHo8Li27Xq9Hm5ubjh8+HC7/W43bPnW+JKSkvD2229j69atCAwMhIODA+Li4jot8m0bw8yZM7Fp06Z2z7m5uXV4jEKhwC+//ILm5mZjYieXy/Hjjz/i8OHDOHjwINatW4f169cjJyen09fS2XvUkWvXrmHYsGHd3p/IEjA5IbIyWVlZCAkJwYsvvmhsu11vRE/Pv3DhQsyaNQsAUFdXZ1LM6evri5aWFuTl5WHMmDEAgHPnzhnn7wCAoKAgVFRUYMCAAVCpVHcdT2RkJH73u98BMCQdZ8+ehZ+fn3GfQYMGQafTmRwXFBSEffv2QaVSYcCA7v2qbO19On36tElP1IABAzB58mRMnjwZr7/+OgYPHoxvvvkG0dHRHV7b398f6enpJm3Hjh1rd73GxkacP38eo0eP7lZ8RJaCo3WIrIy3tzdyc3Nx4MABFBcXY+3atcjJybmn509NTUV+fj5++uknzJs3z2SiMF9fX0yePBnPP/88fvjhB+Tl5eH555+HnZ2dsZdn8uTJCA4ORlRUFA4cOACNRoOjR4/iD3/4A3Jzc3scz6FDh3D06FEUFhZiyZIlqKioMNlHpVLh+++/h0ajQWVlJfR6PV566SVcu3YNc+fONY6IOXjwIBYtWtQumWg1bNgwBAUF4bvvvjO2ZWRk4J133kF+fj4uXLiADz/8EHq9Hj4+Pp1e+4UXXsD58+excuVKFBUVYc+ePdi9e3e76x07dgwymazTWz5ElorJCZGVeeGFFxAdHY3Zs2dj3LhxqKqqMulFuVtvv/02hgwZgpCQEMycORNPPvmkSX0JAHz44YdQKpWYMGECZs2aheeeew5yuRy2trYADLeivvzyS0yYMAGLFi3CyJEjMWfOHGg0GmMNTXetXbsWQUFBePLJJxEeHg5XV1dERUWZ7LNq1SpIpVL4+/tj2LBh0Gq1cHd3x3/+8x/odDo8+eSTCAgIwCuvvAJnZ2fY2HT+q/P555/HP//5T+P24MGDkZqaiokTJ8LPzw87d+5EcnIyHnrooU6v7enpiX379uGLL77Aww8/jJ07d3Y442xycjJ++9vfdlm/QmSJJEJPbm4SEfWCixcvwsPDA19//TUmTZokdjh3pbGxET4+Pvj44497tUfj6tWr8PX1RW5uLtRqda9dh0gMrDkhoj73zTffoK6uDoGBgSgvL8err74KlUqFCRMmiB3aXbO1tcWHH37Y5WRt90JpaSneffddJibUL7HnhIj63IEDBxAfH4+SkhLI5XKEhIRg69at8PLyEjs0IjIDTE6IiIjIrLAgloiIiMwKkxMiIiIyK0xOiIiIyKwwOSEiIiKzwuSEiIiIzAqTEyIiIjIrTE6IiIjIrDA5ISIiIrPC5ISIiIjMyv8DSYmUmUF8oVMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "post_D_true = post_2_dict[\"D_true\"]\n",
    "post_M_true = post_2_dict[\"M_true\"]\n",
    "\n",
    "D_est = jnp.mean(post_D_true, axis=(0, 1))\n",
    "M_est = jnp.mean(post_M_true, axis=(0, 1))\n",
    "plt.plot(dlist[\"M_obs\"], dlist[\"D_obs\"], \"bo\", alpha=0.5)\n",
    "plt.gca().set(xlabel=\"marriage rate (std)\", ylabel=\"divorce rate (std)\")\n",
    "plt.plot(M_est, D_est, \"ko\", mfc=\"none\")\n",
    "for i in range(dlist[\"N\"]):\n",
    "    plt.plot([dlist[\"M_obs\"][i], M_est[i]], [dlist[\"D_obs\"][i], D_est[i]], \"k-\", lw=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above figure demonstrates shrinkage of both divorce rate and marriage rate. Solid points are the observed values. Open points are posterior means. Lines connect pairs of points for the same state. Both variables are shrunk towards the inferred regression relationship."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 15.2 Missing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 15.2.2 Imputing primates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Code 15.16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = pd.read_csv(\"./data/milk.csv\", sep=\";\")\n",
    "d[\"neocortex.prop\"] = d[\"neocortex.perc\"] / 100\n",
    "d[\"logmass\"] = d.mass.apply(np.log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat_list = dict(\n",
    "\tK=jnp.array(d[\"kcal.per.g\"].pipe(lambda x: (x - x.mean()) / x.std()).values),\n",
    "\tB=jnp.array(d[\"neocortex.prop\"].pipe(lambda x: (x - x.mean()) / x.std()).values),\n",
    "\tM=jnp.array(d.logmass.pipe(lambda x: (x - x.mean()) / x.std()).values)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat_list[\"B_impute_idx\"] = jnp.isnan(dat_list['B'])\n",
    "dat_list[\"B_impute_N\"] = dat_list[\"B_impute_idx\"].sum()\n",
    "dat_list[\"N\"] = dat_list['M'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([-0.9239385, -0.9389385], dtype=float32)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def gen_B_likelihood(nu, sigma_B):\n",
    "\t\tdef _generator():\n",
    "\t\t\tyield tfd.Normal(loc=nu, scale=sigma_B, name='B_merged')\n",
    "\t\treturn tfd.JointDistributionCoroutineAutoBatched(_generator)\n",
    "\n",
    "bb = gen_B_likelihood(0., 1.)\n",
    "bbb = bb.experimental_pin(B_merged=jnp.array([.1, .2]))\n",
    "bbb.log_prob()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO\n",
    "TFP should include the fuctionality equivalent to numpyro's (1) numpyro.deterministic() and (2) **obs** argument, such as numpyro.sample(\"x\", dist.Normal(mu, sigma), obs=x)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def model_15_5(M, B_obs, B_impute_idx, B_impute_N):\n",
    "\t\n",
    "# \tdef _generator():\n",
    "# \t\ta = yield tfd.Normal(loc=0.0, scale=0.5, name=\"a\")\n",
    "# \t\tnu = yield tfd.Normal(loc=0.0, scale=0.5, name=\"nu\")\n",
    "# \t\tbB = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bB\")\n",
    "# \t\tbM = yield tfd.Normal(loc=0.0, scale=0.5, name=\"bM\")\n",
    "# \t\tsigma_B = yield tfd.Exponential(rate=1.0, name=\"sigma_B\")\n",
    "# \t\tsigma = yield tfd.Exponential(rate=1.0, name=\"sigma\")\n",
    "# \t\tB_impute = yield tfd.Sample(\n",
    "# \t\t\ttfd.Normal(loc=0.0, scale=1.0), sample_shape=B_impute_N, name=\"B_impute\")\n",
    "# \t\tB_merged = B_obs.at[B_impute_idx].set(B_impute)\n",
    "# \t\tB_ll = tfd.Normal(loc=nu, scale=sigma_B).log_prob(B_merged).sum()\n",
    "# \t\tyield tfd.Deterministic(B_ll, name=\"B_ll\")\n",
    "# \t\tmu = a + bB * B_merged + bM * M\n",
    "# \t\tyield tfd.Normal(loc=mu, scale=sigma, name='K')\n",
    "\n",
    "# \treturn tfd.JointDistributionCoroutineAutoBatched(_generator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# jd_15_5 = model_15_5(dat_list[\"M\"], dat_list[\"B\"], dat_list[\"B_impute_idx\"], \n",
    "#                      dat_list['B_impute_N'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# jd_15_5.sample(seed=key())._fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# @jax.jit\n",
    "# def run_chain(init_state, seed, step_size=1e-3, num_steps=500, burnin=500):\n",
    "\n",
    "#   def target_log_prob(a, nu, bB, bM, sigma_B, sigma, B_impute, B_ll, K):\n",
    "#     # a, nu, bB, bM, sigma_B, sigma, B_impute, B, K = states\n",
    "#     # B_ll = tfd.Sample(tfd.Normal(loc=nu, scale=sigma_B), sample_shape=29).log_prob(B)\n",
    "#     return  B_ll + \\\n",
    "#       jd_15_5.log_prob(\n",
    "#         a=a, nu=nu, bB=bB, bM=bM, sigma_B=sigma_B, sigma=sigma, \n",
    "#         B_impute=B_impute, B_ll=B_ll, K=dat_list[\"K\"])\n",
    "  \n",
    "#   bijectors = [tfb.Identity(), tfb.Identity(), tfb.Identity(), tfb.Identity(),\n",
    "#                tfb.Exp(), tfb.Exp(), tfb.Identity(), tfb.Identity(), tfb.Identity()]\n",
    "\n",
    "#   kernel = tfp.mcmc.NoUTurnSampler(\n",
    "#       target_log_prob_fn=target_log_prob,\n",
    "#       step_size=step_size)\n",
    "  \n",
    "#   kernel = tfp.mcmc.TransformedTransitionKernel(\n",
    "#     inner_kernel=kernel,\n",
    "#     bijector=bijectors)\n",
    "\n",
    "#   kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(\n",
    "#     inner_kernel=kernel,\n",
    "#     num_adaptation_steps=burnin)\n",
    "\n",
    "#   chain_state = tfp.mcmc.sample_chain(\n",
    "#       num_results=num_steps,\n",
    "#       num_burnin_steps=burnin,\n",
    "#       current_state=init_state,\n",
    "#       kernel=kernel,\n",
    "#       seed=seed,\n",
    "#       trace_fn=None)\n",
    "  \n",
    "#   return chain_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# num_chains = 2\n",
    "# init_state = jd_15_5.sample(num_chains, seed=key())\n",
    "# post_5 = run_chain(init_state=init_state, seed=key())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# post_5_dict = {k: jnp.moveaxis(v, 0, 1) for k, v in post_5._asdict().items()}\n",
    "# dict_5 = dict(posterior=post_5_dict)\n",
    "# idata_5 = az.from_dict(**dict_5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summary_5 = az.summary(idata_5, hdi_prob=0.89)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summary_5[summary_5.index.str.startswith(('B'))]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}