{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Characteristics #\n",
    "\n",
    "Defining the various characteristics of the CarBen and Phat distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density ##\n",
    "\n",
    "### Phat ###\n",
    "The Phat distribution density is simply s weighted average of the constituent left and right CarBens.\n",
    "$$\n",
    "f(x) = \\sum\\limits_{i=1}^nw_if_i(x)\n",
    "\\\\f(x) = 0.5*f_{\\textit{left}}(x) + 0.5*f_{\\textit{right}}(x)\n",
    "$$\n",
    "\n",
    "### Right CarBen ###\n",
    "The density of the CarBen is a piecewise function as detailed by [Carreau and Bengio](#references.ipynb).\n",
    "\n",
    "$$\n",
    "f_{\\textit{right}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma} f_{\\mu,\\sigma}(x) & \\text{if } x\\leq a \\\\\n",
    "      \\frac{1}{\\gamma} g_{\\xi,a,b}(x) & \\text{if } x > a\\\\\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "### Left CarBen ###\n",
    "When reflected to the left-side, we simply flip the valence of the $x$ and Pareto location, $a$.\n",
    "\n",
    "$$\n",
    "f_{\\textit{left}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma} g_{\\xi,-a,b}(-x) & \\text{if } x < a \\\\\n",
    "       \\frac{1}{\\gamma} f_{\\mu,\\sigma}(x) & \\text{if } x \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative Distribution Function ##\n",
    "\n",
    "### Phat ###\n",
    "Like the PDF, the Cumulative Distribution Function (CDF) of the Phat distribution is just the weighted average of the component CDFs.\n",
    "\n",
    "$$\n",
    "F(x) = \\sum\\limits_{i=1}^nw_iF_i(x)\n",
    "\\\\f(x) = 0.5*F_{\\textit{left}}(x) + 0.5*F_{\\textit{right}}(x)\n",
    "$$\n",
    "\n",
    "### Right CarBen ###\n",
    "The CDF of the CarBen is also piece-wise, expressed as follows for the right.\n",
    "\n",
    "$$\n",
    "F_{\\textit{right}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma} F_{\\mu,\\sigma}(x) & \\text{if } x\\leq a \\\\\n",
    "      \\frac{1}{\\gamma} \\left(F_{\\mu,\\sigma}(a) + G_{\\xi,a,b}(x)\\right) & \\text{if } x > a\\\\\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "Thus, the CDF is simply the sum of the Gaussian and Pareto components.\n",
    "\n",
    "### Left CarBen ###\n",
    "The formula for the left CarBen is only slight more complicated. As with the other functions, the CDF of the left CarBen is a summation  of its components. The CDF is scaled left-to-right, so in the left instance the first component incorporated is the generalized Pareto. The Pareto tail is a reflected version of its right-tailed cousin, so the CDF of the left is equivalent to the survival function of the right. And we know the reflection is achieved by reversing the valence of $a$ and $x$.\n",
    "\n",
    "$$\n",
    "G_{\\textit{left}}(x) = \\frac{1}{\\gamma}\\left(1 - G_{\\xi,-a,b}(-x)\\right)\n",
    "$$\n",
    "\n",
    "And this handles the CDF up to the junction $a$. Equal to and beyond $a$, we must determine how to sum the CDF of Pareto and the Gaussian. When $x=a$, \n",
    "\n",
    "$$ 1 - G_{\\xi,-a,b}(-a) = 1$$\n",
    "and\n",
    "$$G_{\\textit{left}}(a) = \\frac{1}{\\gamma}$$\n",
    "\n",
    "For the Gaussian, we must exclude the left-tail portion and so calculate the CDF on the interval $[a, x)$, which is equivalent to:\n",
    "\n",
    "$$\n",
    "F_{\\mu,\\sigma}(x) - F_{\\mu,\\sigma}(a)\n",
    "$$\n",
    "\n",
    "This results in:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\gamma}\\left(1+ F_{\\mu,\\sigma}(x) - F_{\\mu,\\sigma}(a)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGACAYAAACeIXc9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACaP0lEQVR4nOzdd3gUxRvA8e+kAikkhCRAQm+hdxClqkiVLlKligVFVGzYEPwpIoIooKJgoSugdKT3TkgghCSUIL2T3pP9/XHJkXJJLskll/J+nuceuNvZ2Xcme3dzs7MzStM0hBBCCCGEaVmYOwAhhBBCiOJIGllCCCGEEPlAGllCCCGEEPlAGllCCCGEEPlAGllCCCGEEPlAGllCCCGEEPlAGllCCCGEEPlAGllCCCGEEPlAGlmiUFBKOSilZiulLiql4pRSmlLqsrnjykxRizcnlFJ7ksszKifb8jkmLflRrSCPK4oOpdTl5HOkk7ljESKFNLLMQClVRin1ilJqg1LqilIqSikVqZQKVkqtVkoNV0qVTrfPb6m+aFIe8Uqp+0qpC0qpf5RSU5RS1Y04vqG8DD0m5VslZLQWeBOoAUQDt4G7xuyYqjx78i+8DLKMVyk1Sik1VSnVNDeZG/n3MfTYY4KyFUu5ed8VUFyZvR8TlVIPlFIHlFJvmSM2IUTeWJk7gJJGKfUssBCokOrlSCAJqJb8GAB8pZQaoWnarnRZxAMPUrIDHIFyQE2gD/C5Umo18KqmafeyCSd1XoZEZlceU1BKNQCeTo6ng6ZpRwriuLllZLyjgI7AZcAnF4e5ncnr5QBrIAYINbA9q7+nsa4AgZnkXySZ4H1XENK/H0sBzsATyY+xSqlOmqYZ9eNDCGF+0pNVgJIvsfyD7oM+EBgBlNc0zV7TNEfACRgI7AEqAR0MZHNI07QKyQ93TdNKo/sg7g6sAjTgOcBHKeWZTUip8zL0+DmvZTZSg+R/Txf2BlayfI83s78JcCg5yapM0vQ3wbFf0DTNS9O0v/OaV2FgovddQUj/fnRKjm0yusZgfWCGmWITQuSCNLIKiFKqCfAjujrfDDTTNG2ppmn3U9JomhaqadoaTdM6A4OBcGPy1jQtRNO0rZqmDQZ6ouvl8ABWm7oc+STlMkiEWaMwXlGLt8TKz/ddQUiO7RtgUfJLz5ozHlE0KaXaKKXqmjuOkkgaWQXnc8AWuA4M1TQtOqvEmqatAmbn9CCapm1F98sXoE3yZZICo5Tqr5TaqpS6q5SKVUpdU0otU0o1N5B2qlJKA35LfqljujEpnQog3oZKqcXJ43JilFIhSqmDSqmXlVLWOY03eSyWhu5SIcCv6dJczu8yJcdaXin1qlJqnVIqQCkVnjz+yF/pBuxXymLfXA9uz0l9ptvPQin1ulLKVykVnXz+bFBKtc1pDOnk+X2X27pUqQZiK6U8lFILlFKXkt8XPjksx+nkf+2ySpTT+k8XY7nk8gQnx3hdKfWzUqpiDmPNltLdMKIppXoa2PZ9qvdLGwPbVyRvm5pF/rkqS27O3/yqQ6XUIPVo7G3NTNL8kZwmWCnlnkV23YGA5LKMU0o55jQeI2MurZTqo5T6QSl1Uil1Kzn+20o3bviJ/DhuoaZpmjzy+YGuVykJ3aW8d3OZx2/J++8xIq0NujE9GrAiL3nlID4L4PfkfDUgAXiY6nki8Eq6fSYDt9CN/dGAuOTnKY/HTV036fZ7LTmulBjDk+NOeb4bKJOTeIHnk/8fl5wmNF2a43ms5z3J+f6WTbpZqcoRD9xPV7Y7QONsjjEqh9tyVJ+p9rNCdzkvdbwPU/2/f6pt1QryfZeXukQ3Hk8DxqO7KUJDNw4sAvDJyfkL/JScxtdU53O6GIen+n8kut7wlP2CAWdTfVYkH3dxct5fGdh2OtWx3zGw/Ubyts6mLEsezt98qUN0Y259kvdfbGD7tORt9wGvbPIax6PPJA2IApYATwLKhH/X0amOkZAcW+p6SASeNeW5VNgfZg+gJDyAYalOsizfDFnkke0Hcbr0y5PTX8trXkYe7/3kPJOAjwCH5Nc9gD9TvcE6GNh3VF7iyU15gL7J+4QB76AbowO6BmpXICh5+0+5iZcsGiN5rOeUfH/LJt1E4AOgEWCV/Jol0ALYmpyHn6EP2Kxiz2xbHuvzw1Tnx2SSv8iA6sAWICTV+6daDuoqz++7vNQlj75ww9E1HB5Pta2WMecvuhtbJvHoy9/g+ZTb+k8V40PgFNA2+XUroDePGrszTXwep7yHjqR73QXdZ0hY8vaN6bbXTn49FihtqrLk8fzNtzoEevGocV891espjZkYoJ2RebkCbwDePHpfaMnxT02dfx7+rq+j+/xvCdiker0xsDP5eIGmPJcK+8PsAZSEB7pLFilviFz9aiDnjawPUr2JrDPJK31PTOrHrzmIzZ5HvTtfGthuCexP3r7PwPZROSmbCerGMtUHY9dM0tRE92s0HqiY03gxcyMrmzxsgbPJ+XTMSeyGtuWlPtFd/kr5Qp2aTaw5bWTl+X2Xl7ok7ZevuxHnb/r3Y0iqcnsDI0x9Pqfa7xbgYmC/t5O3XzJxvVXnUePBPtXrfZNf/wFdL0gIYJFq+7jk7fsN5JmrsuSl/gqiDoHDyfv/kvy8S3IcScBzucyzEfBNcswp51gSut66FwC7fHivpDSgNaCcqfMvrA8Zk1UwXJL/fagln20F4GGq/5fLJI014J7JwzkHx+qC7hd3HDAz/UZN0xKB6clP2yulKqRPU8A6AVUBP03T/jWUQNO0i8ARdL9GOxVYZAVA07RYYHvyU1OMkehE7uvzGcABXc/EnExinZXLuPL9fWdkXf6haVpmU3Kklv79WDbVtnKAm1JKGdivE3k/nxdqqW4GSOWf5H+rK6WyHA+WE5qmBQPXkuN5PNWmjsn/7gYOoKuDpga2780i+5yWpROm+TzIrzr8MPnfF5RSvdHd0GQFTNY07a9c5IemaWc0TXsb8ETXW7Ya3ed3J3TDPm4qpX5RSrXLTf6ZSJkyBXQ9hCWCzJNVsu3VNK2TCfJJGdTuq2naw0zS7EN3ycMyOf1mExw3t1I+1GsrpW5lkS7lS65yPseTL5RSXujGmXRANw+UPbpxHqllOgA+B/JSnynnjo+maZnNy5XVF2qByGNdHjbyMGnej0opS6AK0A3dj5RZ6KZxGJtuP1Ocz8cz2ed6qv87Ydq58/aiu6TbEdiW/FrqRpQnusttHdH15KXfnpmclsVUnwf5Uoeapu1SSu0GOgPrkl/+TtO0HN8YZSDvBGATsEkp5Yzu7tpRQGt059lYpdQUTdO+NCY/pZQrMAbdOVsP3Y/19A2qeCC7ORyLDWlkFYyUXzfOSilVQL1ZqXuiTDFBZVZck/+9nlkCTdNilFL30P1Cd80sXQFJudPHFl082Slj6gCUUscx/GE9S9O03PbcpM5/MPAHut4R0P2CDEXXYwS6RoId2dytZqS81GfKuXAji/SZnlfZMMn7zgR1mavJQ5N7gIOBH5RSl9CN/xqjlPpV07QDqZKa4nw2OG1F8vs25Wmmd4fmUupGFkqpskATIEDTtNtKqZSGVEdgjtKtZlEZ3YDqQwbyS5HTspjq8yA/63AeukYWwHp0q02YlKZpD5VSa9DVQyV0jVySn2dLKdUP3Q0NTqleDuPRVZUy6Hqtzyc37koEuVxYMM4l/2sLFNRcJY2S/72maVp8AR2zVAEdJ69Szvt1mqYpIx5T8yEGVwxfprXPa8bJvyZ/RveBvgrdINRSmqY5a48mNU25NGfo8lNOFYb6NCTP7zsT1WVibo6dWvJlrJRelkHpNhfW+s/OvuR/WyndkkHt0ZUlpXHlg+5Lun3yZdKUXqyTmqaZsketUNefUsoFSN2TZKlpWlJm6XORv61S6jml1AZ0P2jmoGtg3UE3nckfRuTRFvgLXQNrJbq/lZ2maWVTvU9S8vExVexFgfRkFYy96Ab7KXTd3wH5eTCllA3wVPLT/fl5rGQpv9SrZJZAKVWKR2NkzL0sSMr4mEzjzW+aplXLx+y7o2us+aObG8rQB7Ixv9iNlZf6TDkXsrpsmdtLmqZ43xV0XWblCrpZ62uke93s53NuaJoWqJS6ja7+2vKoEbUneXuiUuoA0APdj0ZjLhXmRqGtv+TPzXVAHXR3LjYBeiql2mqaZuxl6MzyfgwYiW7amZQrH/HJx/sV2JSDHqdP0A0F+UXTtBcNHEsB/ZKf+uQh7CJHerIKgKZp13g0Bul1YyeCy2SQqzFeBNyS/78sl3nkRMp4idpKKY9M0nTgUaPeO5M0BSXlw6lxFvHmVcqXsSl6inIqpZv/tKFGQfJ59aQJj5eX+kw5F5pm8b7omMnrWTLR+66g6zIrKXWbvme6IM7n/JLSm9URw42o7LabQqGsv+Rzawm6GyoC0a2X+mfy5v/lMs8qSqkPlVKB6Mr9MroGlh+6uyA9NU3rq2nauhxe0ku56WNlJtsH8OjH0qlchF5kSSOr4HyEbgyHJ7A8+RdKppRSg4C3cnoQpVRX4Ovkp4c1TduU0zxyYRu6bn1rdHPMpI/JEvg4+el+TdOyGlxaEHYCV9H98vo6q4TJg0FzIyz5X6dc7p8XKQPIG2bSUH8R3S3pppKX+kw5d2zRzeGTPr0Nug//3Mrr+66g6zKzuJ7gUSMr/Y+Ugjif80tKg6kXupsggjRNu2lg+wh00z4korvr0JQKa/19g25NzTtAD03THqCbgDQJ6KyUMrpxr5R6Qim1C910E5+j6xl7CCwAWmma1kjTtNmapt3JZawpS425pN+glKqCbjLdFD65PEaRJI2sAqJpmg8wAd3li57AKaXUcKWUfnoFpVRZpVuWZje68R8OxuSdvF9XpdQKdL/cS6P70Bho4mIYlDw+4ovkpxOTfynZJ8fmAawA2vFootL8Yq10S6Bk9bBOHqP2Grq/xZDk5R6apmSilLJWSrVUSs1EN/A4N84m/9s/eUBvQdqBrmwNge+UUk4ASilHpdQ7wHweDQrPs7zUZ/K5kzLtx6dKqbeSx+eglKoG/E0e7u40wfuuQOsyPaVbpqQvuvcQ6GbqXpw6TQGdz5nF10nlbRmslJ6q5ugaOel7qU6guxuvVfJzH03TwjAhc9ZfZpRSE9ENbo9GN0P6peRYz5G73qwu6AbOa+h+2AwBKmmaNkHTtBMmCNk3+d+pSql6oB/r9Txp7669rmlaibmzEJDJSAv6gW6yvZQlb1Ie4TyakDHlcZlUs6NjeMLC2+g+dFPvl4Tui6J8FjGk5LXHhOWyJOOyOg94NPlcIvBqJvuOyks8qcpjzKNTqv1Go+vlSNkWRcYlU7TcxAt4pco7Ht2A0svAgTzW8x6MmIwU3YDV1OV+yKNZw7fyaKLODPmQixnf81if+bKsjined3mpSx5NUtkpm9hSzt/0k5HeTXfcCKBXFvnkuP6NiTGrukc3r1KG91YO/i4K3e38KXkMNZBmW6rt32SRV17LktvzN0/HNZC2X/L5lQj0NbC9fqrzL9PzId0+Y9DNt+WZm/ePEfn35NFnvYauBzhlGZ9N6BqEGulm8C8JD+nJKmCapv2DbuDqBHS9TikT8lmhe7OuBoYCdTVN22cgi9QTFrqgm836Errbej8Eamia9rxWwL8WNE1L1DRtJLres23oZmq2B26i+xXeWtO0BQUZU3Y0TfsV3V1n36LreUpEN6nqfXSNiU/J5V1pmqYFoPv1uBXdB04FdBMeema1n6lomvYWujXzTqH74rBM/v8kdB+IJr+FOrf1qenGfgxAt3zN6eTYEtF9OHfUNG2tCWL7h1y+7wqwLtNPRloeXcPqNLpLRw00TduY2c75eT5nIWX6gyh0NwfkiKb7hk59c46h8VaGxmiZnJnqL43kwejL0F1leiv5vE0fpz+6O/kAphszdlfTtMWapv1P041TNDlNNyylF3AM3XskAd1l3WGapvXk0c0aPvlx/MJMJbdChRBCiBxRSv0IvISuh2myueMRorCRRpYQQohcUUqdQ9dDW10zbukgIUoUuVwohBAix5RuolYv4CdpYAlhmFGNLKVUN6VUoFLqglLqfQPb31JK+SulTiuldiqlqqbalqiU8kl+rDdl8EIIIcxD07S7mm4GdJMv8SJEcZHt5cLkOY6C0A3ivYZuEcwhyYPvUtJ0Bo5qmhallHoF3V0Wzydvi9A0Lc9LhQghhBBCFCXG9GS1Bi5omnZJ07Q4dDO69kmdQNO03ZqmRSU/PUIB3UElhBBCCFFYGbN2oQe6iS1TXAPaZJF+LLAl1fNSSqkT6G7pnGHollSl1Hh0t0hjZ2fXwsvLy4iwhBBCCCHM6+TJk/c0TXM1tM2kC0QrpYajW6W+Y6qXq2qadl0pVQPYpZQ6o2naxdT7aZq2EFgI0LJlS+3ECVNMQCuEEEIIkb+UUv9lts2Yy4XXSbushWfya+kP8jS6yTB7a5oWm/K6pmnXk/+9hG5Ct2ZGRS2EEEIIUYQZ08g6DtRWSlVPXqx1MLrZxfWUUs3QLQDZW0u1wKRSylkpZZv8//LoVurO8azAQgghhBBFTbaXCzVNS1BKvQb8i245icWapp1VSk0DTmiath7dyuX2wF/JM/xf0TStN1AP+EkplYSuQTcj9V2JQgghhBDFVaGb8V3GZAkhhBCiqFBKndQ0raWhbTLjuxBCCCFEPpBGlhBCCCFEPpBGlhBCCCFEPpBGlhBCCCFEPpBGlhBCCCFEPpBGlhBCCCFEPpBGlhBCCCFEPjDp2oVCiOIvLiGJO+ExPIiMIyQiBqWglI01TmVs8HQuQ2kbS3OHKIQQhYI0soQQWboXEcu2U8Hs8btGwJ1IrkdbkJhFJ7iTZQJ1XWxoU7M8TzevSSNPZ5JXghBCiBKlSDay4uLiuHnzJrGxsRS2GevzSimFjY0NFSpUoFSpUuYOR5RQUXEJrDkUxIrDlzgXChoKB6Kpa3GLNpZ3qWr1kNqOkThbxZCkKWI1S+7ElyY41okLseU4c8eDo3eS+O7wHdysE+hWrzyjuzSmuquDuYtmdvv376dHjx7UrFmTqKgoKlasyKpVq6hQoQIDBgzg/fffp1WrVuzdu5dXXnmFs2fPopTizp07PP7448yaNYu+ffuaNKbY2FieeOIJdu7cSdmyZfOUl7HlAwqsjAVVPiBDGU2pqJ4TJVmRa2TFxcURFBSEk5MT7u7uWFgUr2FlmqYRERFBUFAQderUkYaWKFD3wmP4dv1J1vjdJ1qzpBIhjC/jT8+yF2hoexsLBbfDdO85d/ukLPO6k2DHv+E1+CukPkt9LVhyei+PuVvyRq/mPFbbvSCKUyh5e3vTp08fli5diqZp9OrVi/nz59OrVy8iIiL0X84dO3bE2dmZv//+m+7du/Pss88yYcIEo79M79+/j4uLi1FpbW1tGT58OLNnz+azzz7LbdEA48sHeStjfpWvWrVqXL58Ocflmz59OkePHs1QxtyIj4/nwIEDdOrUKU0vcFE9J0qyItfIunLlCi4uLlSsWNHcoeQbR0dHrK2tuXz5Ml5eXuYOR5QAEbEJzFx7jJW+94jXLOhofZ6XXbxpU/oaub3S52YVyQjnMzxjeZZ7SXasiW3KqtvNGLzoBK1dLfhs8GPU83A2bUGKAG9vbxo1agToeq6rVKlCQkICCxcuZOjQoWnSTpkyhalTp7J8+XJatWrFm2++mWXeUVFRrF+/nmXLluHv78/FixfTbB8zZgxVqlRh6tSpnD9/np49e7Jy5UqaN2/O4MGDadmypUkaWcaWL6dlzK58AEeOHGHy5MlERUURHh7Om2++yauvvprv5QMylDGzWAy5evUqW7ZsYcuWLRw5coRWrVrRsWPHDJfai+I5UZIVuQWi/fz8qFWrVrHv4UlISMDPz4+mTZuaOxRRjGmaxl+HzvPFZn9CEq152uoc77sdpJbtg0z30fdkOWbdk2UofXSSFT/eb8HCiMeIwYb+Xo58NuRx7G2L3O+9XGvUqBFfffUVPXr0IDAwkE6dOvHPP/8wdOhQNm3alOGHVfXq1alduzZbtmzB0jLjTQUJCQns2LGDZcuWcfDgQbp168aQIUNo165dhi/o69ev06JFC7Zs2cKwYcNYuHAh7dq102+vUaMGmzdvThPDwIEDuXDhgsGyHD58mNKlS+epfNmVMSflA3j48CEODg5YWVkRGhpKtWrVePjwYablSy+7nqzMytemTRtq1qyZpoxZxZJi2rRp/PnnnyQkJNC9e3e6d+9Ox44dsbW1zTQGc58TIq2sFogucp9sSUlJxe4SoSEWFhYkJRn3JSZEbtwOieK1RXs5fjeJKkkP+Lrsdhrb3IJYuB2b+XvsdnjO7h5Mmz6JobbH6Wl9hrmh7Vh7rgl7pm9ixnPN6NKkSi5LUnTExMQQEBDAu+++y8cff4yTkxOLFi2iTZs2XLt2DXf3tJdRFyxYQExMDLGxsQa/TAFatGjB3bt3+f777/n111+xssr8Y93Dw4ORI0fSvn171qxZk+bLFKBChQpcu3YtzRfq6tWr8618xpQxJ+UD2LRpEz///DMhISEkJiam+b4wVD6A3r17c+XKFQBu3Lih/3FrZWVF6h/9WZUPyFDGrGJJERUVRVRUFG5ubjg5OeHk5IS1tXWm5SsM54QwXpFrZAkh8m798Yt88LcfcYkwSu2nj60vlWwSCuz4ZS1i+MR5B73jzvJZWHdeXHGG533/Y/rQJ7CxKr4/onx9fXFzc8PPzy/DttKlSxMTE6N/vmHDBubMmYO3tzcdOnTg0KFDPP744xn2W7RoEUuWLOHdd99lzZo1DBkyhK5du2JjY5Mh7Z07d9i8eTP29vZUqZKxURsTE5OhZyonPVk5KZ+xZcxJ+TZu3MgPP/zAmjVrqFChAn/99Rc//fRTluUDWL9+vf7/1apVw8fHx2B5sypf+jJmF0uKGTNmMGPGDPz9/dmyZQtTpkwhICCAjh07smTJkjQNs8JyTgjjFd9PswK2detW6tatS61atZgxY4ZR+4wZMwY3NzcaNmyYz9EJoZOQmMSUZQeZuCaAitp9Nnv8xisuJ6nkmIC7Y5JxD4dE3cME6buWv87u6r8yqNRJVvmH0WPGJq7cizB3NeUbb29vWrY0eFWBRo0aERgYCMCJEyd49dVXWb9+PRUrVuS9997jf//7n8H9WrZsydy5cwkKCmLkyJGsXr2a2rVrM27cuDTpQkJC6N69O5999hlTp07l3XffTbM9MTGRS5cuZfg8Wr16NT4+PgYf6b98jS1fTspobPlSjt+iRQsqVKjA9evX+eijj/TxZFa+nMiqfOnLmFUshtSvX5+3336bHTt2EBgYyPPPP5/m0l5hOieE8aSRZQKJiYlMmDCBLVu24O/vz4oVK/D39892v1GjRrF169YCiFAIuB8Ry8DZW1l+JoQBtqfYXHUJtbMYe1VQbFUiMyvu5DuXv7keEU+P2Ts5GHDD3GHli6y+pPv378+///7L5cuXGTBgAEuWLKFevXoAvPDCC5w5cwZfX99M87a0tKRr1678/vvvBAQE0KVLF/22qKgoevXqxauvvkr//v0ZN24cQUFB7N69W5/m4MGDtGnTJk+36xtTPiBXZcyqfClGjBjBnj17aNKkCZMnT8bT05MWLVoUSPnSlzGrWFJr1aoVFSpUSPOoWbMm48eP1w8ZKcrnRElX5Aa+nz59Gi8vrwzdnp06deLHH3/Ey8uL+/fv07Fjx0y7dH19fXn99de5d+8eAQEBaJrGxx9/zLRp03IV8+HDh5k6dar+zfXll18C8MEHH2S77+XLl+nVq1eGWJOSkvDx8aF58+a5ikmI1P67G87z83ZxL0Yx1flfhjs/Ot/yMpDd1Okvxjkz7NpA7qqyTO9Ri6Ed6hl1jOIgLCyMdu3acfToUbNcnhk6dChjx47lqaeeypf8i3v5wPxlNLWCqLPiIKuB78WmJ+vChQvUqVMH0DXEUm6xTS8mJobnn3+eWbNm4e/vz4cffsjkyZMN3qLavn17mjZtmuGxY8eONOmuX79O5cqV9c89PT25fv26CUuXd/b29ibN74svvkjz3NC4gNwYNWoU1atXp2nTpjRv3pzDhw9nmd6Y43777bdERUWZJL6iyPe/e/T+dhfhsQksrbAyTQOrsKlp85Al5ZZSj2tM2XyJOetPmjukAuPo6Mg333xDcHBwgR87NjaWDh065OuXaXEvH5i3jKZWUHVW3BWLge///fcfHh4e+gGCp0+fpnHjxgbT7tixg+bNm9O6dWsAGjduzNatWw3eCrx///78C7qI++KLL5gyZYr++aFDh0yW99dff83AgQPZtm0bL730EqdPn840rTHH/fbbbxk+fDhlypQxWYxFxX7/67y45CRltBj+rLSKurb3zR1SthwtYvmp3Go+Ce/B3EMQFn2ETwa1KRFL8xi6BFYQbG1tefnll/P9OMW9fGC+MppaQdZZcVYserJ8fX3TNKpOnjyZaSPLz88vTS+Xt7d3ppfkjO3J8vDw4OrVq/rn165dw8PDIy9Fyhd79uyhV69e+uevvfYav/32G6C7o+bTTz+lefPmNGrUiICAAAAiIiIYPXo0jRo1onHjxqxZs4b333+f6OhomjZtyrBhw4BHPWWapvHOO+/QsGFDGjVqxKpVq/TH7tSpEwMHDsTLy4thw4ZluyRShw4d9Hc1zZ49m4YNG9KwYUO+/fZbfZqU42aW/3fffceNGzfo3LkznTt3znslFiH7/a8z5o+TlNdC2ej5R5FoYKWwUkks9NhAb9vT/HrqPu8tOVjsltASQhR/xaIny8fHR3/b7Pnz51m3bh2ff/45AE899RR//PGHvtHj4uLCrl27AAgKCmLt2rWZ9oYY25PVqlUrzp8/T3BwMB4eHqxcuZLly5frt6ePobAqX7483t7eLFiwgFmzZvHLL78wffp0ypYty5kzZwDd5HoDBgxg3rx5Bm9zXrt2LT4+Pvj6+nLv3j1atWpFhw4dADh16hRnz56lUqVKPPHEExw8eDDDnCypbdiwgUaNGnHy5El+/fVXjh49iqZptGnTho4dO9KsWbM06Q3lP3HiRGbPns3u3bspX7686SqrkDt47gZj/ziJuwplredy3KyK3uVSS6Uxt+JW7G7FscK/JWrpQWYMf6JE9GgJIYqHYtOTlZSURJMmTZg2bRr169fn999/JykpiQsXLlCuXDl92iFDhhAREUHDhg0ZP348K1asMHotp8xYWVkxb948unbtSr169Rg0aBANGjQAMBhD6ljatm1LYGAgnp6eLFq0KE9x5FX//v0B3UR2KTMe79ixgwkTJujTODtnvQzKgQMHGDJkCJaWlri7u9OxY0eOHz8OQOvWrfH09MTCwoKmTZtmOqvyO++8Q9OmTVm4cCGLFi3iwIED9OvXDzs7O+zt7enfv7/BBrCx+Rd3R4NuMfr3E7gSylqPFUWygZVCKfiiwi6eL3WSVWdD+fyvo+YOSQghjFYserJOnz6Nt7c3Dg4OaV738/NjwIABae7ysLe3Z8OGDSaPoUePHvTo0SPD6/7+/hliSLFixQqTx5EVKyurNLPIp58YMGUZB0tLS/1aXKaUepmIrI6RMiYrxc6dO02af3F27vpDRv96FGctgrWVV+BmHWnukPJMKfiywk4ib9qwyLsRDmW8mdRL7roVQhR+Rb4nKzw8HKVUhgYWQMOGDZk9e7YZoipcMaSoWrUq/v7+xMbGEhISYlTjpUuXLsyfP1//PGXdLWtra+Lj4zOkb9++PatWrSIxMZG7d++yb98+/U0GudW+fXv++ecfoqKiiIyM5O+//6Z9+/ZG7+/g4EB4eHieYigKbjyMZNgP+7BKiuUHpz9xLwYNrBQWCuZW3Epn60C+PXCTZXvPmTskIYTIVpHvyXJwcCAoKMjcYRRqCQkJ2NraUrlyZQYNGkTDhg2pXr16hjFNhnz00UdMmDCBhg0bYmlpyaeffkr//v0ZP348jRs3pnnz5ixbtkyfvl+/fhw+fJgmTZqglGLmzJlUqFBBP5A+N5o3b86oUaP0jbVx48YZFXuK8ePH061bNypVqpRmor3iJDQqjue/20lUfBJfWP+NVXSkfn6q7ORtLcKCTT/dYRMvPyjNJ1uSqOzqQIf6njnKW+SvBw8e0KVLFwIDA4mIKL4z9wthrGIzGWlxY8rJSH19fXnxxRc5duyYCSIThU1CYhIDZ2/lzL0EvrFfQ4143Rxt7g6JRu2f0qgpKunDk2x54cEwHlo4sm5iR+pUdDIqX5H/4uPjCQ8PZ9CgQRnuwhaiuMpqMtIi35Mlsvbjjz/y3XffpZn2QBQvHyw9iM99jc+d/6Wv8+Ucz8ieoqikdyealWX+otf1EQz/YS+b3+lKeYdSOcpb5A9ra2uDN/kIUVIV+TFZImsvv/wy/v7+PPPMM+YOReSD33ed5a9zYQwpfaJQz+RualVsQlnkvoYHsYpxP+4iITFnDThz2r9/Pw4ODjRt2pQ6derQsWNHbt26BcCAAQP0d+Pu3buX+vXr6+cHu3PnDrVq1eKff/4xeUyxsbG0bNmS0NBQk+edWuryZVUPec0bCq7+TFl3pq6TrKSvL1MzV/0X1LlsLGlkCVFEHTt/m2nbLtHM8grT3YvnWLOstCpzk6nO/+JzX+PTlVkvv1SYeHt706dPH3x8fAgMDMTe3p758+dz9OhRIiIiaNWqFQAdO3bE2dmZv//+m+joaJ599lkmTJhA3759jTrO/fvGTz5ra2vL8OHD8/UmnfTly6weTJE35K3+8qvuqlWrluV2U9ZJVgzVV27Ex8eze/dugxMFm6v+C+JczglpZAlRBN0Ji2H8b0coTziLPf7BShWusZUFZbizHz0tfVh2JoS1R86bOxyjeHt761edUEpRpUoVEhISWLhwIUOHDk2TdsqUKXz55ZeMGDGCVq1a8eabb2aZd1RUFCtXruTZZ581eFfvmDFjmDp1KqCbuLlOnTp4e3sDMHjw4Hydqy99+TKrB1PknSIn9Zdd3QEcOXKEdu3a0bx5c2rXrs2CBQsA09WdKeskK+nrK7NyGXL16lUWLlxIv379qFKlCnPmzMl0NQZT139m52/6+s/vczknpJGVia1bt1K3bl1q1arFjBkzMmy/evUqnTt3pn79+jRo0IC5c+fqt1WrVo1GjRrRtGlTWrY0OBZOiFxLStJ46addRCbArxXX4GwZk/1OxdgHZXdRm5t88M85Am6EmDucbKX+Ig0MDOSff/6hb9++7NmzhzZt2qRJ27NnT+7du0dYWFiaz5jUEhIS2Lp1KyNGjKBhw4bs27ePd999V78kVWrTp0/nxx9/5NSpU/Tp04fFixfrb66pUKECNjY2Ge4EHjhwoMHlxZo2bUp0dHSGYzz99NOcOnWKp59+Gj+/R5ew05cvs3oAWLduncG8U0ybNo3Y2FiGDh1KdHS0wbozpv5yUncAdevWZc+ePXh7e3PixAk+/PDDLOsup7Kqk8GDB+fojs1nn302023p6yuzcqU2bdo0GjZsSJcuXTh37hyvvPIKly9fZv369fp1g9Mzdf1ndv6mr39T/T1MQQa+G5CYmMiECRPYvn07np6etGrVit69e1O/fn19GisrK7755huaN29OeHg4LVq0oEuXLvo0hXEZl86dOzNlyhS6dOnCRx99RGhoKN9//32O81m9ejWzZs0iOjoaBwcH/v77b1xdXfMhYmHInPUnOHVf41OnndQvdc/c4ZidtUpittM/DAsZyYs/72P7lJ6Uss7ZtBEFJSYmhoCAAN59910+/vhjnJycWLRoEW3atOHatWu4u7unSb9gwQJiYmKIjY3F0tJwmVq0aMHdu3f5/vvv+fXXX7Gyyvxj3cPDg5EjR9K+fXvWrFmTYVmrChUqcO3aNby8vPSvrV69OkdlzOyuwtTly6oeABYvXkzv3r0zPcYnn3wC6Ho/SpcubbDuIPv6y0ndAWzatImff/6ZkJAQEhMT0zQuDNUdQO/evbly5QoAN27coGnTpoDuOyT1nfTZ1Ul4eLh+rdbsJCUlZXkHfvr6yqpcKaKiooiKisLNzQ0nJyecnJywtrbOMg5T139W52/6+s/s71HQik1P1tatW/W/rtq0aZNmZvOcOnbsGLVq1aJGjRrY2NgwePBg1q1blyZNxYoV9b8AHRwcqFevHtevX89TGfLbZ599xv/+9z+WLVvGqVOnMtxxaOyC2J07d+bIkSP4+vrSpUsX/vzzzwIsRcl24sJN5h+5RUfrIEY5+5o7nELDzTKSb8pv4mq0Je8vMW7NUXPw9fXFzc0NPz8/Tp48yc6dO/UrRZQuXTrNKgwbNmxgzpw5eHt7c+PGjUzXWF20aBHPPfcc7777Li+88AIbNmwgLi7OYNo7d+6wefNm7O3tqVKlSobtMTExGVanyElPllIq00fq8mVVD3/99ReHDh2ic+fOnD17lldffZUnnniCdu3aER0dTVJSEgMGDCAxMVH/pZy+7oytv5zU3caNG/nhhx9YtWoVvr6+fPrpp2nm6zNUdwDr16/Hx8cHHx8fKlWqpP9/+qmKsqqTxMREbt++zZAhQ2jSpIl+bsiPP/6Ytm3b0q1bN6Kiopg6dSrt2rXj7bffpkaNGqxcuVJ/2WzdunXMmzcvQ31lV64UM2bM4NKlSyxevBh7e3umTJlClSpVGDZsmMHvW1PXP2R9/qav/8z+HgWt2DSyXn/9dbZs2YKPjw9Hjx7N0BI3tgEBcP36dSpXrqx/7unpmWUD6vLly5w6dUr/i0MpxTPPPEOLFi1YuHChiUqYdx06dEDTNGbPns3KlSsz/LLYv3+//gMg9ePpp59Ok+63336jdevWNGnShAULFlCqlNw+XxDCY+KZ8PsRnIlkbsWtyDrJaXVxCGZY6eP8ExTJP0cNX24wN29v70yHEDRq1IjAwEAATpw4wauvvsr69eupWLEi7733Hv/73/8M7teyZUvmzp1LUFAQI0eOZPXq1dSuXZtx48alSRcSEkL37t357LPPmDp1Ku+++26a7YmJiVy6dImGDRumeX316tUGPxd8fHwyfIlt3ryZsWPHomkaN27cYPjw4WiahqZpacqXVT307NmTbt26sWfPHho0aMDXX3/NwYMH9ZcfL126RPXq1bl48SI1a9bMUHc5qT9j6y4l5hYtWlChQgWuX7/ORx99pC9DZnWXE1nVycWLF7G0tOSPP/5g1qxZrFixgu3btxMbG8vhw4fp3bs333zzDRERERw4cAAnJyfq169Ps2bNOH36NKBr0Lz44osZ6iurchlSv3593n77bXbs2EFgYCDPP/98hkXb86P+szp/09e/Kf4eplJsLhf26NGDxo0bM2zYMINzQhlaUNgUIiIiGDBgAN9++y2Ojo6AbpFkDw8P7ty5Q5cuXfDy8qJDhw75cvycOHPmDDdv3sTFxcXgMkTt27c3uPzMrFmz9A2tP/74g2PHjrFr1y7s7e3p0KGDfjFskb/eWrSLO3FWLK+wGqcSPg4rM5+67+XYf1WY8k8CLWtXxLOcnblDSiOrL9L+/fvz77//Uq1aNQYMGMCSJUuoV68eAC+88ALTpk3D19eXJk2aGNzf0tKSrl270rVrV6Kjo1m/fr1+W1RUFL169eLVV1+lf//+JCQkMGfOHHbv3k3nzp0BOHjwIG3atKFs2bK5Ll/37t31PyxnzZqV5oswpXydOnXKsh7Onj2r/0y5cuUKkydP5s6dO1y6dIkXX3yRkydP0qBBA/z9/fX1kzrvy5cv57j+sqq7FCNGjKBPnz40adKE+vXr4+npSYsWLUxWd1nVib+/Py+88ALW1taULl2acuXKsWnTJiZMmACAhYUF69evZ+nSpQBomka9evWoXbs2ly5dYsOGDXTp0kW/vmvq+sqqXKm1atWKq1evGozv+vXr+h/t+VH/2Z2/6evfFH8PUykWjaxDhw6haRo3b97M9JquMQ2IFB4eHmlOpmvXruHh4ZFh3/j4eAYMGMCwYcPo379/mv0B3Nzc6NevH8eOHTN7I+vmzZsMGzaMdevWMXHiRLZu3Uq3bt3SpDGmIXrmzBkef/xx7O3tWbNmDYcOHdIP1BT5Z/WxS2y/msBL9kdoW+aaucMptGxUEj9XWke366MYv3APG97tgaVF4eny+/nnnzPdNnbsWNq1a8cnn3zCf//9l2abjY2NflyPMUqXLs3zzz+vf16mTBkOHDigf25lZZWm5wd0Exen793KDTs7O4KDg7l7926az4bU5cuqHs6dO0ft2rUB3RCHjz/+mJo1a/LUU0/h4eHBH3/8wZNPPsnOnTv1n92p865WrVqe6i993aWoXr26vlcoPWPr7vLly5luy6pO/P39cXFxAWDVqlVMnDiRuXPnkpCQQFxcHGvWrKFhw4bExcURGRnJ8uXLeeutt7CwsEApxS+//MKqVav0+aWur6zKlZqxc2rlR/1nd/6mr39TncumUCwaWX/99Rd16tTBysoKTdMIDw/X9yqlyElPVqtWrTh//jzBwcF4eHiwcuVKli9fniaNpmmMHTuWevXq8dZbb+lfj4yMJCkpCQcHByIjI9m2bZt+kKa5REVF0b9/f7755hvq1avHxx9/zHvvvZehkWWMUaNG0b9/f5YtW8YzzzxDjRo1sLMrXL0FxU3wnTCmrDlN5YT7vGB7ONs1CYvSWoT5kb40YbxVagdfhPTg+82nmNQr70tTFQRHR0e++eYbgoOD09xkUxBiY2Pp0KEDTz31VJ7zqlatGq+99lqGu7KNLZ+Xlxfjxo3jzp07PP7444waNYrOnTvrP9PPnTvHa6+9xvz583n99ddzlHd+MGXdZebcuXO4ubnRt29fnnzySWrXrs2LL77I6NGjAXjvvfcoX748gwcPplq1ajg6Ourrq0qVKtSpUyfNsA5z1peppa//gvh75ESxWLvw2LFjjB07Vj+4csGCBQa7O3Ni8+bNTJo0icTERMaMGaO/pbVHjx788ssvXLp0ifbt29OoUSP9+K8vvvgCLy8v+vXrB+huTx06dKjB22GzY8q1C0XRcOHCBdauXcuNGze4du0afn5+JGkaVUbM4kpYAi+H/c2wmoa761Mz99qChSG9psHEB33xpjobJ7bHq5KzUXmJvFu0aBFbtmzJ8V2JIn/079+f5cuXy9jZfFTs1y5s3bo1Z86cMWmePXr00N/ZkdrmzZsBqFSpUqYTsPn6yl1fIue8vb15//3305xX9dt05kKEBY/FX6IcUTla/6+orEWYX+nn2W3lqatjmbD4AP9+0BMry2Jzn0+hppRi2LBh5g6jxLtz5w5jxoxhzJgx0sAyI/nUEaKQaNWqVZrFdZ1dyhP+xKtUUuFUVMZPQih03KyimF7uXy5GWDBr3Ulzh1NinDt3Tj/gWZiPm5sbGzduTDNeWBS8YtGTJUqOEydOULZsWf3A2OLA39+fGTNmsGnTJvr168eyZctQSuH8zGuUtrbA08EKYs0dZdHUt2wQG8LPsfBYEj1b3adRZZcCO/bevXsN3mxT3HXs2JELFy5kOmu3EIWFg4MDHTt2zNdjSCNLFCmBgYFMmjSJZ599Vn8nUVF17NgxvvzySw4dOsQbb7zBd999h5OTE87OzqzaeQKtcnMqli+DZfxDc4dapM2qsJ1OV6sy6fdDbJvSq8DuNgwPDy8U8/QIIQxLvexTfpHLhaJIGTZsGEFBQXh4eNCiRQvGjx/PqVOnzB2W0TRN0996PnDgQJ588kmCg4OZMmUKTk5OANTo8gKqy9t42kGp5HltRO45W8XwSbkdXIywYMEWH3OHI4QoQaSRJYocZ2dnpk+fTmBgIJ6envTt25eWLVvy008/ERYWZu7wDAoODmbmzJk0b96cCRMmMHz4cC5cuMDrr79OmTJl9OlOX3nArF2XcbWMxTG50SXyrr9jAK0tLzFv/1Wu3pfxbUKIgiGNLCNt3bqVunXrUqtWrQzzv6SXmJhIs2bN6NWrl/61atWq0ahRI5o2bZrlkgXCeOXLl+eTTz7h0qVLfP7552zbto0qVaowZswYDh8+nOndnwXl6tWrzJ49mzZt2tC6dWsuXrzIrFmzOHv2LKNGjcqwgGtkbAIvLT6IjZaAq2v5DEtViNxTCr6puA2NJN76bZ/Zzw0hRMkgY7KMkJiYyIQJE9i+fTuenp60atWK3r17ZzqJ29y5c6lXr16GXpXdu3dTvnz5ggi5RLG0tKRbt25069aN27dv89tvv+mXoHjmmWfo0KED7du3x9XVNV/jSExMJCAggJ07d7Jq1SoCAgLo27cv06dPp3PnztmuWP/OkoPcjFbUKmeHlVXOJuQU2atsHcYbDvv5+u6TrD1ygQFti8/NE0KIwqnY9GStXr2axx57jCZNmtCuXTvu3r1rsryPHTtGrVq1qFGjBjY2NgwePJh169YZTHvt2jU2bdpkcIFLkf/c3d157733CAoKYtGiRbi7u/PLL79Qq1Yt6tWrx0svvcTSpUv577//8tSbkZCQgJ+fH7///jsTJ07kiSeeoGzZsvTp0wdvb28+/PBDbt68yaJFi3jmmWeybWD9efg8my9E4FEqgdKpLh8K03rJ5STVuMNnG/0JjYo3dzginZdeeokmTZrwyiuvmDsUIUyi2PRkde7cmYEDBwK69a7+/PNP/eKZKXKyfmFq169fp3Llyvrnnp6eHD161GDaSZMmMXPmzAzHUUrxzDPPoJTipZdeYvz48UaXTeScUoq2bdvStm1bPvjgAxITEzl9+jT79u3j77//5u233yY0NBR3d3eDDzs7O8LCwggLCyM0NDTNvw8ePODcuXNUqlSJFi1a0KJFC/r27Uvz5s31g9dz4vLdcD5eH4CzisPZxc30lSH0rJTGpw5bGRM2nE9XHuLbMfl7+7bImdGjR/Pcc8+xZs0ac4cihEkUm0bWb7/9xqpVq4iNjeXWrVt88cUXGdLkZP3C3Ni4cSNubm60aNGCPXv2pNl24MABPDw8uHPnDl26dMHLy8vsi0aXJJaWljRr1oxmzZrxxhtvALo1HW/fvp3hERQURFRUFI6OjpQtWxYPDw/q16+vf162bFnq1atnkhXe4xKSGLdwH1pSItUdNGwjb2dI4xBxC5vYCJKStGzXLYTCv7agudO7xN6ni+bLuqCmjLp8j6bV5BJ+YfHYY49x5MgRc4chhMkUi0bWH3/8wbFjx9i1axf29vZ06NCBBg0aZEhnbE/W/Pnz9Suib968GQ8PD65efbRm3LVr1/Dw8MiQz8GDB1m/fj2bN28mJiaGsLAwhg8fztKlS/Xp3dzc6NevH8eOHZNGlpmVKVOG6tWrU716dbPFMHXVES6EQ71SsdhayjisgjLK9jCH47x4b/lhtrzfC4sCmjvLWKGhoTz11FOsWbOGqlWrmjucNF599VWaNWvGiy++aO5QirWcngPF5e+Sm3O/MJe9WDSyzpw5w+OPP469vT1r1qzh0KFDNGrUKEM6Y3uyJkyYkOZSo5ubG+fPnyc4OBgPDw9WrlzJ8uXLM+z35Zdf8uWXXwKwZ88eZs2axdKlS4mMjCQpKQkHBwciIyPZtm0bn3zySS5LK4qL9SeCWX7mIRWs47Aq70FWEwvYRd/HIlHJ2oUmSu9OFO8l7uGjkO6s2B/AsI4FuwzMsGHDOHToEKDrZfXw8OCll15i6NChgO6HXufOnfVfMpqm0a1bN9q3b89HH32kz+f48eOMGDGCTz75RL9vfps4cSKDBw/m+eefx9HRsUCOmV5SUhJnz57Vf84fPXqUn3/+GT8/P27fvs3XX3+tHz5iaqn/dqmlXK0wlfTnQHby++9SWMsNmZf99OnTNG7c2GSx5UaxGPg+atQoFixYQOvWrTl16hQ1atTAzs7OZPlbWVkxb948unbtSr169Rg0aJC+p6xHjx7cuHEjy/1v375Nu3btaNKkCa1bt6Znz55069bNZPGJoue/exG8u+YMZYnBxVUuV5nDUKcz1FU3+WpbEGExBTsI3s/Pj0mTJnHs2DF2795Nx44d+eijjzh79izR0dGsWrWKQYMG6dMrpZgwYQLLly/nwYMHAFy8eJHx48czduzYXDWwwsPDczWvnJeXF1WqVOGff/7J8b6msnTp0jRfwlFRUdStW5dPPvkkx4shT548mW+//dbo9H5+fkyePJljx46leZiyoWHoHMhOfv9dCmu5IfOy3717l8DAQJPFlxtGNbKUUt2UUoFKqQtKqfcNbH9LKeWvlDqtlNqplKqaattIpdT55MdIUwafokGDBgQGBnLs2DE+//xzgoKCTH6MHj16EBQUxMWLF/nwww/1r2/evJlKlSplSN+pUyc2btwIQI0aNfD19cXX15ezZ8+m2V+UPDHxiYz6cS9JSRoV3FywkPmwzMJCwQz37YQnWPK/PwtuHNB///1HWFgYLVq0wNXVlcqVK/P666+jaRp+fn7s3r0bpVSG+fR69eqFu7s7ixYt4u7du4wePZpOnTrxzjvvGH3sxMRE9u7dyxtvvEHr1q3x9/fPkGbz5s3UrVuXa9eu6V/77LPP6Nixo/6u7aeffpoNGzbksgbyJjg4mLi4uDQ9Fp07d+add96hR48eWFjkX99Byt+uTZs2uLq6pnmYkqFzwJx/F3OWO8X3339Pt27daNCgAS1atGDy5MnExMTotxsq+1NPPcXSpUtNGmNOZXs2KqUsgflAd6A+MEQplX6CqFNAS03TGgOrgZnJ+5YDPgXaAK2BT5VSzqYLX4ii54NlhwiOAE8nW2ysi8UV+yKrWelbPGt7mr/8Qzh3PaRAjnnmzBkA6tV7dIny1q1bgG6C3ePHj9OwYcMMk9FaWFgwYcIElixZwpgxY6hcuTJfffWVUccMCgriyy+/5PHHH+f111+ndOnS/Pbbb7Rp0yZD2u7du1O3bl3mzZsHwMKFC9mwYQO//fab/ku1SZMm+Pr6pvmSS23+/Pk0aNAgy8exY8cy7Dd8+HAmTJjAnj17aNu2Ld7e3hnS/Pjjj/Ts2dOocpvamTNnsLS0NDjm15QMnQOm+LvkljnLnSIxMZHPP/+cf//9l7lz53LgwAEWL16s355Z2Z2dnfHx8cnXuLNizCd8a+CCpmmXAJRSK4E+gP4nkKZpu1OlPwIMT/5/V2C7pmkPkvfdDnQDVuQ9dCGKnlWHzvN3QBhVVRh29pWz30Hku6kV9rHzihdTlh9i7eTu+T7Tvp+fH+7u7ri4uAC6XoJp06ZRtWpV2rVrx6pVq3B3dze477PPPsvUqVOJiIhg2bJlGVYNSO3hw4f8888/rF27lsDAQDp06MAnn3zC008/jW0Wa2IqpZg8eTJjx46latWqLFiwgKVLl6a5QcTNzY34+Hhu375tcOzMsGHDsm0IVahQIcNr2fU6pIzFqlixYpbp8oufnx+JiYlpelo8PDzYtm2bSY9z/fr1DOeAKf4uuWXOcqeYNGmS/v+enp507tyZixcv6l/LrOxNmjRh+/btNG3a1KSxGsuYRpYHcDXV82voeqYyMxbYksW+GS7gKqXGA+MBqlSpYkRIQhQ9gTdD+XhDAC5aFBUdrYk0d0ACgHKW0UxwOMDM+0+x3fcKzzTN37v5/Pz8uHv3Lg0aNCAxMRFN0+jatStz587F1taWmJiYTFeG+PTTT0lISODhw4fZXhb7/fffmTt3Li1atGDXrl14enoaHWOHDh1o3Lgx33zzDb/88gtNmjRJsz1l3FNmPSZOTk65mjMuO0FBQfrGaW7Nnz+fBQsW6J/HxcWhlNLfUQ7w66+/0rp16wz7+vn50b17d9599139azkdA2aMzM6BvP5dZs2axfz587M89ooVK3jsscfSvGbuct+4cYOFCxdy5MgRbt26RXx8PHFxcWnmm8ys7HXr1uWHH34weazGMum1CqXUcKAlkKMZ/jRNWwgsBGjZsqUsKiaKnai4BMYs3IdFUiLVHCyRYViFyziXUywJb8Fn//jwZKPKWFnm37ies2fP6gerlypVCjc3tzQNpnLlyhEaGpphv7lz57J161bWrl3LuHHj+OOPP3j11VczPc6QIUOwsrJi7dq1dO3alWeeeYZ+/frxxBNPYJnNdCGHDh3i3LlzaJpm8EsvJCREH6sh6RsyhqRvyGQ3lUpwcDC3b9/O8/x06XvZvvrqK9zd3Rk1apT+NUO9bKD7202cOJFq1arlKYbsZHYO5PXvMmbMGPr165flsQ2NMTZnuR8+fEifPn1o3bo1H3zwARUrVsTS0pI+ffqkWdous7I7Ojpy+3bG+QcLijGfJNeB1Nc1PJNfS0Mp9TTwIdBb07TYnOwrRHGmaRpvLN7H9SiFh3MZbGQ6rELHRiXxkcsursdY8evOs/l2nKtXrxISEsITTzxBtWrVqFChQoYeqfr163PhwoU0r61atYoffviBhQsXUr9+fcaPH8/ixYuJjo7O9Fju7u689tpr7Nq1i6VLl2JnZ8fEiRNp27Ytn3/+ucFB7wD+/v68/PLLTJ06lWeeeYaZM2dmSBMUFESFChUyHfg8bNgwNm3alOUj/a3127Zt06+mERwczKlTp/Dy8iIgIIDg4GBA90VaJo/LTjk5OVGtWjX9w87OLsNrhnppUv52WY1L6t+/v34Qup+fX5bLA/n4+NC/f3+efPJJ/Xx9KTdLGToHTPF3KVeuHDVr1szyUbp06UJV7l27dhEbG8u8efPo2LEjderUwcfHh8jIyDSNrMzKbmdnZ7DBWlCMaWQdB2orpaorpWyAwcD61AmUUs2An9A1sO6k2vQv8IxSyjl5wPszya+Z3J49exgxYoRJ89y6dSt169alVq1azJgxI8P2mJgYWrduTZMmTWjQoAGffvqpfltISAgDBw7Ey8uLevXqcfjwYZPGJoqOhdv92H45Gs/SSZSxk3UJC6seDhdoZHGV7/ZcIjyfpnRIGfSe1dw9HTp04MKFCzx8+BDQ3XH1ySef8M033+h7foYMGYKFhYXB+foMadasGZ9//jlHjx5l6tSpBAcH06dPnwyDz69du8bo0aMZN24cgwYNYtKkSRw4cCDDLOzHjx/PcjLl9I0WQ4/0DZnatWvTu3dv/Wfl999/z/jx49OMH7O3tzd4KSwyMhJ/f3/8/f1JSkrixo0b+Pv7c/266X7TnzlzBqVUmi/21DRN4/79+/ov+cDAQLy8vAymjY6O5s033+SLL75g165dDBgwgOnTp+t72NKfA6b6u+SGOcsNuoHrUVFRbNu2jcuXL/Pbb7/x1VdfYW9vn6ZnLbOyR0dHY29vn9vi51m2jSxN0xKA19A1js4Bf2qadlYpNU0p1Ts52deAPfCXUspHKbU+ed8HwHR0DbXjwLSUQfCm5uvrS7NmzUyWX2JiIhMmTGDLli34+/uzYsWKDL/8bG1t2bVrF76+vvj4+LB161b9Sf/GG2/QrVs3AgIC8PX1TXMnkSg5DgXd5Ktd/1HeIgancnJjbWGmFHzutouIREtmrs1455sp+Pn5UblyZZydMz8XvLy8aNKkCRs2bMDX15fXXnuNd999N80lLltbW1588UV+/vlnYmNjM80rPVtbW3r06MGiRYs4fPgwtWvX1m8LCQlh1KhRPPXUU0ycOBHQjWfp0aNHml6T2NhYtm3bxuDBg3NSdKPUqVOH4OBgrly5wqlTp+jbt2+a7a6urgZX7Thz5gw9e/akZ8+exMTEMGfOHHr27MmcOXNMFpufnx/VqlXL9Av7ypUrada4DQgIoG7dugbT7tu3j9atW+sbI7Vr1yYsLEx/00Xqc8DcfxdzlTtF586dGTJkCJMnT2bgwIFcvnyZPn364OXlpd8vq7KHh4ebfKqJnDBqTJamaZuBzele+yTV/zNdXVnTtMXA4sy2m4qPjw/lypWjTZs23L17l8WLF9OpU6dc53fs2DFq1apFjRo1ABg8eDDr1q1L05pXSulPvPj4eOLj41FKERoayr59+/jtt98AsLGxyfIuIFE83QyJ5uXfj2OnxVPdLgmr5HUJHSJu5SgfWbuw4NJX4C6dLPxZcboe47tGUNnFtL+A33333TSDhzPzxhtv8Nlnn7F9+3bOnjV8+fLFF1/M0zIi6cf0ODk5sWPHjgzpUqYMSLFq1SqaNGli0h+1KapVq8auXbuYOXMmkydPznCnZ926dfXTXaT22GOP6S8p5tSsWbOMSpfd3y5948Lf3z/TSWLPnz9PnTp19M/PnTuXYaxUyjkwbNgws/5dzFluS0tLlFJMnz6d6dOnZxpDVmW/dOkSzZs3z3Tf/FYsZnwHXU+Wg4MDR48e5ccff+Tjjz/OkKZ9+/Y0bdo0w8PQCXz9+vU0rXNPT0+DXc+JiYk0bdoUNzc3unTpQps2bQgODsbV1ZXRo0fTrFkzxo0bR2Sk3EtWksQmJDLqh91EJWjUtEvCqpCtjScy92bZ/SgtiU9WmO8Sf8eOHRkxYgQ3b940WwyZsbKyYurUqfmSd/Xq1Tlw4ABRUVE8/vjjGbbb2tpSuXJls46xyUxAQIC+VzHl6kbKVALDhg1L0zgsW7YsAQEB+rQXLlzIUN6cngP5+XfJirnLDVmX/cyZMzz55JM5KZJJFYuZEOPj47l37x5TpkwBoGnTpty7dy9DOmPXLswJS0tLfHx8CAkJoV+/fvj5+ZGQkIC3tzfff/89bdq04Y033mDGjBlZtsRF8fL+0sMEhmrUKGtDooOLwXUJI+wN38FkiKxdWHDp3QnjhbjjLLrWFu/gezSvbp5lj0aPHm2W42YnP9dITJnF/b333ss0zbhx49i0aVOBrdVorICAAGxtbenevTstW7bEzc2NdevW0bdvXy5fvpxmSos+ffqwY8cOunXrhru7O/Pnz8fa2jpDnjk5B8xVH+YuN2Rd9uDgYF566aUc5WdKxaKRFRAQQK1atfSX5Ly9vTPMHwK6nixD1/NnzZrF00+nveLp4eHB1auPpvi6du1alms0OTk50blzZ7Zu3crw4cPx9PTUz6Y8cOBAgwPnRfH0+55z/B0QiodNHHYO5hsLIHLvTdej/PlfUz776yjr3jXP7OIlUXx8PI6OjpmO6QFo2bIlR44cIT4+3uAXtLmcP3+eTZs2ZZjo9cKFC3Tv3j3NQH9HR0d+//33gg4xXxTmcp84cYLevXvn+wTDWSkWjSwfHx+Cg4OJjY0lPj6ezz77zOCAx5z0ZLVq1Yrz588THByMh4cHK1euzHAnz927d7G2tsbJyYno6Gi2b9/Oe++9R4UKFahcuTKBgYHUrVuXnTt3ZnpnhiheTl66y7StFyin4nAu72bucEQu2VvE8arjIb568BS7/a7SuaHMzl8QLl68SM2aNbNNN2bMGE6dOmVwwlBziIqKwtLS0uBM+rVq1eKjjz4yQ1T5r7CXOzIyko4dczRtp8kVi0aWr68v/fv35/HHHyc6OpqPP/44w4y1OWVlZcW8efPo2rUriYmJjBkzRj9PSI8ePfjll1+4d+8eI0eOJDExkaSkJAYNGkSvXr0A3e3Hw4YNIy4ujho1avDrr7/muZyicLsbHsO4xYcpRSLu7uXN+utJ5N2Ycj4sDmvF9L9P0amBp/w9C4CXlxe//PJLtunKlClTaBpYoIvn33/zZXaiQq2wl9vcDSwoJo0sY+8OyakePXrQo0ePDK9v3qy70bJSpUqcOnXK4L5NmzblxIkT+RKXKHziEpIYuWAXYfFQw9URq2xm1BaFn61K5C3nA3zwsAd/H71A/8dqZ7+TEEKkUmzuLhTCXDRN4+0/DuD/UKOKoxWlslh8VxQtg8qexVO7z8zN/sQn5mxgvRBCSCNLiDz6cbsfG4LC8bCNxz757ihRPFgqjdfs93Irzoolu/NvuR0hRPEkjSwh8mCX3zVm7voPdy0cZxcXc4cj8kEn20vU4Qbf7b5EdFyiucMRQhQh0sgSIpcu3A5lwjJvHLVoqpa1koHRxZRSMMl+LyGJVvyw1cfc4QghihBpZAmRC6FR8Qz/YR9JSUnUdABLaV8Vay1sr9Pa8hK/HLlOaFT+LB4thCh+isXdhfkpMTGRli1b4uHhwcaNGzNs37p1K2+88QaJiYmMGzeO999/36htouhKSExi9A+7uBMNDcrEUj76bo72l7ULi2b618rs54XwGszdeJJPBuVtihghRMkgjaxszJ07l3r16hEWFpZhW2JiIhMmTGD79u14enrSqlUrevfuTf369bPcJoq2KcsP4303gTo2kdjbWEGcuSMSBaGuzV3aqiCWn6rJaz1iKWef9V2kDg4O+Pn5FVB0QoiccnBwyPdjFJtG1tmzZ3njjTe4cuUKI0aM4M6dO7zwwgu0atUq13leu3aNTZs28eGHHzJ79uwM248dO0atWrWoUaMGAIMHD2bdunXUr18/y22i6Pp5hx9/ng2hok08tq6eadYkzMlahDlNL2sXFo70n9rup/v12sxZf4LpQ5/IMn1hmAhRCGFexaKRFRMTw3PPPcdff/1FjRo18PLyokWLFhkaWDlZuxBg0qRJzJw50+A+ANevX6dy5UfLbXh6enL06NFst4miaZvvVb7YcRkXi1hcZMmcEsnL9j5PWgfw5+k6THw2BleHUtnvJIQosYpFI2vHjh00a9ZMv+xNXFwcb7/9doZ0OVm7cOPGjbi5udGiRQv27NljqlBFEeV39QGvr/TBkTjc3d3kTsIS7AO3g3S5XpfZ/xznyxHtzR2OEKIQKxaNLB8fH5o1awbAjRs3sLe354knMnbl56Qn6+DBg6xfv57NmzcTExNDWFgYw4cPZ+nSpfo0Hh4eXL16Vf/82rVreHh4ZLtNFC03QqIY8dMBLJMSqVihPJYW0sAqyWrZPuAZm3OsPuvFm2ExuDlKb5YQwrBi0ciysbHh+vXrAHzwwQfExRkeiZyTnqwvv/ySL7/8EoA9e/Ywa9asNA0sgFatWnH+/HmCg4Px8PBg5cqVLF++PNttouiIiE1g6LzdRMRD9fKOWFvJmoQC3nc7yPZr9fj672N8PbKDucMRQhRSxWKerKFDh7Jv3z7q1q1LkyZNaNu2LZMmTcq34/Xo0YMbN25gZWXFvHnz6Nq1K/Xq1WPQoEH6S5ZZbRNFQ0JiEmN+2M1/4RpVnGywLSVrEgqd6jYhdLc9y9/nQrn5MMrc4QghCqli0ZPl6enJyZMn8y3/Tp060alTJ/3zzZs36//fo0cPevToYXC/rLaJwk3TNN5bdphjt+KoaqdRxt7e3CGJQuY910NsudaAmX8fY86YTuYORwhRCBWLniwhTG3eltOs8Q+hkk08js7O5g5HFEJVbELpZXuGDUHhXHsQae5whBCFkDSyhEjnn2MX+WbfVVwtYylXXhZ9Fpl71+0QaBpfrT1m7lCEEIWQNLKESOVAwE3eXuuPsxaFm7urTNUgsuRpHc6ztmfYfCFCerOEEBkUizFZQpjC2WsPGff7ceySYmhqeZeYSON+g+RmLcKcppe1Cwtv+pGlj7I+thGz/jnOtzI2SwiRivRkCQFcvR/J0B/3Y5GUQFPLe1hbaOYOSRQRFa3CecryLBuDwrkZEm3ucIQQhYj0ZKUzZswY/WzvmS3uGhISwrhx4/Dz80MpxeLFi2nbti2BgYE8//zz+nSXLl1i2rRp+TqdhMi7B5FxPD9vN9HxGtVdnYiPTyKe/F2LMKfpZe3Cwp3+w1KH2XGtIbPXHZd5s4QQetKTlc6oUaPYunVrlmneeOMNunXrRkBAAL6+vtSrVw+AunXr4uPjg4+PDydPnqRMmTL069evIMIWuRQVl8Dg73ZwO0qjSrky2NjamDskUQRVtQmlq60//5wL4U5YjLnDEUIUEsWikXX+/HmqVavGhQsXAIiPj6dp06ZplrUxVocOHShXrlym20NDQ9m3bx9jx44FdLPNOzk5ZUi3c+dOatasSdWqVXMcgygY8YlJjF6wm/MhSVQta0XpMqXNHZIowt5xPUyCppiz7oS5QxFCFBLFopFVu3Ztxo8fz7///gvAvHnz6N27N5UrVwZ0axY2bdo0w2PHjh05PlZwcDCurq6MHj2aZs2aMW7cOCIjM95VtHLlSoYMGZK3gol8o2kab/y6n6O34qhiD3YODuYOSRRxNWwe8pRNIGvOPuB+RKy5wxFCFALFopEF0LBhQwIDA3nw4AGLFi3ivffe02/bv3+//jJe6kf6RaGNkZCQgLe3N6+88gqnTp3Czs6OGTNmpEkTFxfH+vXree655/JcLmF6mqbxyaqjbL4QgYdtAo4GeiKFyI13XQ8Rr1nw3Yb8W4FCCFF0FJtGVp06dQgMDGTq1KlMnjwZOzs7/TZT9mR5enri6elJmzZtABg4cCDe3t5p0mzZsoXmzZvj7u6et0KJfPHN+lMs8blPBet4nF0yvzQsRE7Vsb1PW4vzrDp9j5AowwvVCyFKjmJzd2HNmjXx9vYmNDSUb7/9Ns22/fv3m+w4FSpUoHLlygQGBlK3bl127txJ/fr106RZsWKFXCospH7a7se8wzdxs4ylvKtMNipM7xWHQ4wIq8O8jd58NOgxc4cjhDCjYtOTZW1tjaOjIzNmzMDCIvfFGjJkiH46Bk9PTxYtWgToFnu+ceMGAN9//z3Dhg2jcePG+Pj4MGXKFP3+kZGRbN++nf79++etQMLkVhwI4sud/1HeIkZmcxf5po71PVpzgeWnbhMWE2/ucIQQZlRserJAd1dhx44d85THihUrDL6+efNm/f+bNm3KiROG7yCys7Pj/v37eYpBmN6GE5eZsjGIcioW9wpu0sAS+eoVh8OMDq/FD5t9eK9/K3OHI4Qwk2LTk3X58mWqVq0qX54ig91+15i0+gxlVSwVKrhhIeeIyGf1bW7T2vISf5y4SURsgrnDEUKYSbHpyapWrZpJx16J4uHI+VuMX+qNQ2I0tZ0ssIq6ne0+BbEWYU7Ty9qFRS/9i2UO8WJ4DX7edoY3n22WozyEEMVDsenJEiI93//uM2rxMUolxdLE6j5WFtKDJQpOY5tbNOQKvx25Qkx8ornDEUKYQZHrybK2tiYuLg4bm+K9/ElcXBzW1tbmDqPI8r/2gKE/HcRCS6S2oyLO0q1QrUWY0/SydmHRTP+u5SFeuDOYpXvPMe7phjnaXwhR9BW5niwXFxcuX75MfHzxvWsnISGB4OBgXFxczB1KkRR4M5RBPxwgKTERTzdnbHN2tUcIk2lf5gp11C1+2neRhMScNdKEEEVfkevJcnNzIykpibNnz1KqVKliN9Bd0zRiYmIoX748FStWNHc4Rc75W6E8N28vCQlJVHVzwsbaCmSFE2EmSsEkl0O8eq8/qw+fZ3C7uuYOSQhRgIpcIwt0E4K6uLgQGxuLpmnmDseklFLY2NgU+8uh+eHSnXAGzttLXEISVdzKYmMjl1uF+XWzv0Dl+/eZvzORQY/XwULGBgpRYhTJRhboxmbJmCWR4vLdcAZ8t5uY+CSqujpiKw0sUUhYKHjd+TDvPujFZu/L9GpZ3dwhCSEKSJEbkyVEelfuRdD/u91Exmu6Bpat9AKKwqW/4zncCGPuVr9i1/suhMicNLJEkXbtQST9v9tFRJxG1fL20sAShZKV0ni57BHOR1iw1/+GucMRQhQQaWSJIuvq/Qj6ztlJaKyugVWqlK25QxIiU0OdzuBEJHM2+Zg7FCFEAZFGliiSgu+G0+fbXYTFaVQrbycNLFHolbJIZIzDcXwfwIlLd8wdjhCiAEgjSxQ552+F0m/ubiLjkqjmak+pUqXMHZIQRhlTzgc7Yvhmnbe5QxFCFIAie3ehKJkCboQwaP4+4uIT8LJPpEz8Q8hmXtrCuBZhTtPL2oXFJX0Cfa1Osex2W85dD6Geh1OO8hZCFC3SkyWKDL+rDxg4bx/x8fG0sLxNGWuZb0gUPSMdTmKjxfPNP8fNHYoQIp9JT5YoEnz+u8/Qnw5CUiJeDkkkWbkWqrUF8zu9rF1YfNK7E8WAaF9WXW3OlXsRVClvn6N8hRBFh/RkiULv+MU7DP5R18Cq4uZMKSvpwRJF28Tyx7BAY/a6E+YORQiRj6SRJQq1vf7XGf7LESy1RCq7lcPaWjpfRdFX0TqC7jZn2Xg+nNuh0eYORwiRT6SRJQqtTd6XGfuHNzZaApXdy2FjnbMByEIUZm+5HiFRg+82yp2GQhRXRjWylFLdlFKBSqkLSqn3DWzvoJTyVkolKKUGptuWqJTySX6sN1XgonhbeSCI1//0w544PCu4YG0lDSxRvFS3CaGzdSCr/e4TGpXNLbJCiCIp20aWUsoSmA90B+oDQ5RS9dMluwKMApYbyCJa07SmyY/eeYxXlAA/bjvD+xvP46RiqVTRDStLaWCJ4ult1yPEapYs2HLK3KEIIfKBMT1ZrYELmqZd0jQtDlgJ9EmdQNO0y5qmnQZydvuNEKlomsaMf04yY9cVylvEULGiG5YWMshdFF8NSt2lOZdYdvImUXEJ5g5HCGFixjSyPICrqZ5fS37NWKWUUieUUkeUUn0NJVBKjU9Oc+Lu3bs5yFoUF0lJGh8sP8KPR27hZhlLhQpuWChpYIni7yWHI0QkWbF451lzhyKEMLGCGPheVdO0lsBQ4FulVM30CTRNW6hpWktN01q6uroWQEiiMElITOK1xftYeeYBFa3jcHN3RUkDS5QQTW1uUI9rLD74H3EJcjFAiOLEmEbWdaByqueeya8ZRdO068n/XgL2AM1yEJ8o5qLjEhkxbyebL0TgWSoBF9fy0sASJc44uyM8SLBk5cEgc4cihDAhYyYdOg7UVkpVR9e4GoyuVypbSilnIErTtFilVHngCWBmboMVxcuDyDiGfLeDwJAkGmi3cLR1gsjb2e5XGNcWzO/0snZh8U5fS7tKVe0OP+xKZFj7ujIWUYhiIttPbE3TEoDXgH+Bc8CfmqadVUpNU0r1BlBKtVJKXQOeA35SSqUMLqgHnFBK+QK7gRmapvnnR0FE0XLlXgQ9v/6XCyEJNOM6HtYx5g5JCLNRCkaXOcLNWEs2nLxs7nCEECZi1PTZmqZtBjane+2TVP8/ju4yYvr9DgGN8hijKGZOX7nP8IUHiYnXqF6uDKWSXAincK0VWNjSy9qFxT/9CIdAfrocwtx/z9KnZTW5bC5EMSAzvosCtfvsdZ774RAJ8YlUc3WgdJnS5g5JiELBUmm8XPYwwRGKXWdvmDscIYQJSCNLFJiVB4MYu+QUtlocVSqUw9bWxtwhCVGoDHE+SzkimLPJ19yhCCFMQBpZIt9pmsacDd68v0E3i7tnBVdZJkcIA2xUEmPsj+L3UOPoxTvmDkcIkUfSyBL5Kj4xiUm/7WfuwZu4WsbqZnG3lNNOiMyMcTmNPTHM3uBj7lCEEHkk33Yi34RGxzPo2+2sCwzHwzYed3dXmcVdiGyUsYhnaOkTHL0Vj//1EHOHI4TIA2lkiXyRMkWD7504qtlBufIyyagQxnq1/Elsieeb9d7mDkUIkQfSyBImd+zCHXrM2cWdyERqlLPFwamsuUMSokhxsoqlv60Pu/6LIvhuhLnDEULkkjSyhEmtOXKBob8chYR4qrk5UrpMGXOHJESRNMn1GJZaEt9uPGXuUIQQuSSNLGESmqYx85+TvP1PII5aNF4OSdjaWJs7LCGKLHfrSLrb+LEpMITbYbIighBFkVEzvguRlZj4RCYu3se24CgqEk4Dy3tYREO4pXFjsArjWoGFLb2sXVgy048sc4yNIY34frMPnw9+LEf5CiHMT3qyRJ7cDImm19db2RYcRXXLcKo6WiFr2wphGpWtQnlcBfGX711CouLMHY4QIoekJ0vk2vFLdxm7+DDR8VDDyQY7e08i0S12C4Vr7b+inl7WLiy56adYH6LnLS9+2u7He32a5yhPIYR5SU+WyJXl+4MYvPAIiQkJVHe1x87e3twhCVEsNSh9j9YWl1hy7DqRsQnmDkcIkQPSyBI5kpCYxPvLDjNl03nKqjiqVHCRNQiFyGdvlz9ERKIFv+8NMHcoQogckEaWMFpIVByDvt3OyjMPqGgdR6WKblhZyhqEQuS3NnY3aKCu88v+YGITEs0djhDCSNLIEkbxv/aAZ776F5/kGdzLu7nKDO5CFKBJ5Q7xIN6CPw9fNHcoQggjSSNLZOvPQ+fpM+8gYbGJ1HApJTO4C2EGTzsEU527LNgVSGKSZu5whBBGkEaWyFRcQhLvLj3Eu+uDcNSiqObuTOnSpc0dlhAlklLwmvMhbkZbsMH7P3OHI4QwgjSyhEE3Q6LpPetf/vR7SBUVRq2yllhbyfgrIcypb9kg3LWHfPfvWTRNerOEKOykkSUyOBBwk2e+3s6Fh/FUd1B4OFrLBKNCFAKWSmOE7VEuhcMu/5vmDkcIkQ1pZAk9TdOYu8mHF347CYkJ1HC1x97R0dxhCSFS6Wvvj7MWwZxNp80dihAiGzLjuwAgPCaeCb/sZd+1WCoQQRVHKyzjH0B84VzLr6Sll7ULJX2KBxGKXpxgyYNOHL14jzY1y+dofyFEwZGeLIHvlfs89cUW9l+NpqZlGNUcrTBybWchhBk8a3MGey2a2Ru8zR2KECIL0pNVgmmaxk/b/fh612VsSKSmiz2lSzsTmUn6wrSWX0lLL2sXSvq0EhkadZKFt9px9noIDTyccnQsIUTBkJ6sEio0Kp4X5u9ixq4rOFvEUa1COZmeQYgi5NXyJ7HV4pm9XnqzhCispJFVAp24dJenZmzl4NUoqpROpEIFWR5HiKLGySqWAaVOseu/KILvRpg7HCGEAdLIKkGSkjS+2+zLoIVHiI6No6ZLacqWKyfL4whRRE1yPY6llsQcGZslRKEkjawS4n5ELEO+287sfddw1SKp75hEqdKlzB2WECIP3Kwj6WFzhs1BodwOizF3OEKIdKSRVQLs9rvGkzP+5fjNWKqU0ahe1gormV1UiGLhLbejJGkwd6P0ZglR2EgjqxiLiU/kg+WHGb3Ul6SEBGq6lqGssxNydVCI4qOaTSidrQNZc+YeDyPjzB2OECIVaWQVUwE3Q+g6YysrTj+gonUsVSq6UsrW1txhCSHywduuh4nVLPnxX19zhyKESEUaWcWMpmn8uN2PXnP3czsinpplLSnv5oaFXB4UotiqX+oebSwvsvTETSJjE8wdjhAimTSyipE7YTEM+nY7M3b+R1mLOKpVcKaMvb25wxJCFIDJroeITLLklx1nzB2KECKZzPheTGz2vsy7q08TnQi1rSNxKWONirlrMG1hXJtP0medXtYulPTZpa/CbRpyhcWHEhjfpTGlbWTuOyHMTXqyirjQ6Hhe/nkPr/55FsvEWBrZxVDezloGtwtRAr1if5DQRCt+3X3W3KEIIZCerCJtl9913l7lTUgcVLOMwK6CJ4lKYezcz4VpbT5JnzVZu1DSG5O+J1f5MeoaC/cnMubJBpSylt4sIcxJerKKoIjYBN78bT9jlvqQEB9PI7toKjpYYyHdV0KUeO+4HCAkwZI/9p4zdyhClHjSyCpiDgTcpNP/tvD3uVAqWcdRpZIrdtbyZxRC6LS3u0I9dYMf914iNiHR3OEIUaLJt3MRER2XyHtLDzH8N29i4uKpVc4WFzdX6b0SQqShFLxTfj8P4i1Ztj/I3OEIUaJJI6sIOBx0i85fbGaV30MqWsdStaIrpcuUNndYQohCqrPdf9RRt1iw+zxxCTkb6yWEMB1pZBVi4THxvPX7AYYsPklYTAI1na1lYlEhRLaUgrfL7edenCUrD0pvlhDmIo2sQupf3yt0+N9W1vqHUNE6jmoVy1OmTBlzhyWEKCKecQimlrrNvJ1BxCdKb5YQ5iCNrELmfkQsY3/cxUsrzpAQH08tl1KUd3PFUnqvhBA5oBS8Ve4Ad+Is+fPwBXOHI0SJJI2sQkLTNFYeOk+HL/5ld3Ak1S1CqVLJldKlS5k7NCFEEdXd4SJVtbt8vz2ABOnNEqLASSOrELh6P5Ln5mzj/fVBWCfF09guhgoONnLnoBAiT5SCF8sc5FasJauPXjR3OEKUODLjuxnFJyaxYOtp5u+/QpIGtawjKV/GGsfI2xBrfD6Fca09SW/a9LJ2oaTPbfqGCcF4JN3ju22JPPdYLRl6IEQBkp4sMzl6/g5P/m8zc/Zfx0mLorFDHK6y5qAQwsQsFAy2PMKNGEv+Pia9WUIUJOnJKmAPI+P4eOURNp4Ppwzx1ChbGjt7T+KB+HRpC9PaeZLevOll7UJJn5f0A7UAVj18nG//9adf65rSmyVEAZGerAKSlKSxdH8g7f63lc1BoXgkT8tgZ29v7tCEEMWchYI3nQ5wLdqStTI2S4gCI42sAnDuRgi9Zm7ho00XsNHiqelqRzmZlkEIUYD6lQ2girrHnH/95U5DIQqINLLyUXhMPFOWH6Lnd/u5+DCOanbgUdGNUrY25g5NCFHCWCiYXG4/N2IsWSXzZglRIKSRlQ+SkjSW7Q/kielbWO77gIpE0MghHgensigZ2S6EMJNnHc5TQ91h7vYAmQVeiAIgjSwT8758j64ztvDhpgtYJMZTu3xpqpS1xtpSGldCCPNSCt4pt487sZYs3R9o7nCEKPaMamQppboppQKVUheUUu8b2N5BKeWtlEpQSg1Mt22kUup88mOkqQIvbO6ExfDKz3vo/+NRroXFUd0ePCu5UaqUrblDE0IIvW4Ol6itbjFv53liExLNHY4QxVq2jSyllCUwH+gO1AeGKKXqp0t2BRgFLE+3bzngU6AN0Br4VCnlnPewC4+4hCS+2+xL+y+38++FcDxsYqle0QX7snJpUAhR+CgFU8rv5X68Jb/vCTB3OEIUa8b0ZLUGLmiadknTtDhgJdAndQJN0y5rmnYaSH+RvyuwXdO0B5qmPQS2A91MELfZaZrGNt+rdPx8E7P3XcOBWGq7O1DO1U3uGhRCFGqd7P6jgbrOgj0XiYmX3iwh8osxjSwP4Gqq59eSXzNGXvYttM5ee0i/b7YyfsVpwmISqOFkRcVKFbCxlrldhRCFX0pvVkiCJb/sOGvucIQotgpFq0ApNR4YD1ClShUzR5O5W6HRTPvzKJsvRmCrxVPLOgaXMtZYaOEQEZ7pfoVxLTxJX7TSy9qFkt7U6Wtxg4Zc4acDiYx+sj52toXi60CIYsWYnqzrQOVUzz2TXzOGUftqmrZQ07SWmqa1dHV1NTLrghMZm8D/1hynw4wdbL0QTq3EOzxucQNXO2vkyqAQoqh63X4/4YmW/PjvaXOHIkSxZMxPl+NAbaVUdXQNpMHAUCPz/xf4ItVg92eAD3IcpZkkJi+FM3t7EKEJlrhZxOLsXp5ysbHE4Fio1raT9MU7vaxdKOnzI31XrtMq6hKLjyQy/plGOJSyzlGeQoisZduTpWlaAvAaugbTOeBPTdPOKqWmKaV6AyilWimlrgHPAT8ppc4m7/sAmI6uoXYcmJb8WqGmaRrbT1+l0+eb+HTLRSwS4qlVzhb3ijLuSghRvHzotp/IJEu+33TK3KEIUewY1WLQNG0zsDnda5+k+v9xdJcCDe27GFichxgL1LGLd5m2+jh+DzXsiKe6Y2ns7N1kOgYhRLHUtPRt2lsG8fuJRF7qGouLvcztJ4SpyIzvyc7dCGHI3G0M+vkYFx7GUcVOo1ql8tg72EsDSwhRrH3ivo+4JMVXfx83dyhCFCsl/trXf/cj+fyvo2y/HIW1lkjl0uDo7IKFjGgXQpQQtW0f0MPGj7VnG/L6g0gql7Mzd0hCFAsltifrTlgMb/9+gCe/3sWu4AiqqjBqVXTGyUUaWEKIkucj9wMoLYnP/zpq7lCEKDZKXE9WTHwiM9edYsnJmyQkKSpYx1GxlIatpQ0RliW2zSmEKOEqWkcwqLQ3y4Nb4X8jhPqVnMwdkhBFXolrVVhZKDb5XqO8FkFtd3vKu7tjm7O5/YQQolh6x/UwpbQ4pv15zNyhCFEslLxGlqUFb1e6RS27JGxlOgYhhNBzsopliM0xjtyK58j52+YOR4gir8Q1sgBsLTRzhyCEEIXSSIeTlNUi+Wz1CTRNPiuFyIsS2ZWTlJhI6egHJFjp5oMpjGvVSXpJnzq9rF0o6QsqfWiExiB1mJ9Dn2arz1W6Nyu868kKUdiVyJ4sIYQQmeth7Ye7FsIX631ITJLeLCFyq0T2ZFlYWhJdOuPag4VprTpJL+lTk7ULJX1Bp5+i7eGNB31Zui+AkZ3q5WhfIYSO9GQJIYTIoLdjEF7qBrN3nCcyNsHc4QhRJEkjSwghRAZKwTS3nYQmWDJnoyweLURuSCNLCCGEQa3L3KSTVQC/n7jFjZAoc4cjRJEjjSwhhBCZmua+By1JY9qfstyOEDkljSwhhBCZqmITxvOlT/LvpUh8rzwwdzhCFCnSyBJCCJGld90OYUcMH688KhOUCpED0sgSQgiRpbKWcbzueIDTD5LY4nPV3OEIUWRII0sIIUS2xpbzoRIPmb7Oh/jEnM25JURJJY0sIYQQ2bJSGh+77OJmjCW/7PAzdzhCFAklcsZ3WbtQ0he19LJ2oaQvDOmbacE05Arf703g+XZelLOzyVHeQpQ00pMlhBDCKErBB447iUm04LNVR8wdjhCFXonsyZK1CyV9UUsvaxdK+sKS3p27DIg9xeqg5oy98oDGVcrlKF8hShLpyRJCCJEjH7kfwJ4Y3l9xRKZ0ECIL0sgSQgiRI2UtY5lcdg/+DzVWHb5g7nCEKLSkkSWEECLHRjifoZa6zYzN54iITTB3OEIUStLIEkIIkWMWCr5020ZIgiVf/X3C3OEIUSiVyIHvQoi8cR48F5zcjEpbPmX8tJE/6SR9AaePuA3rJxqXWTqtytzkSeXHcp/6jOwcSi33srnKR4jiShpZQogcs3R0g5gQo9ImJSrdPpbGDZCW9AWc3t7dqHwy83bZfRx6WIf3lh5m9VtdUUrlKT8hihNpZOXBl/0b4eZgm206pSUCoCnjJgWU9IU3/Z3wWD5Ye8aofIQoCcpbRjHW9gDz7z7J38cu0b9NTXOHJEShIY2sPHBzsCU0Oj7bdJZJujSJFtZG5SvpC296YxrVQpQ0w+xPseVhQ6ZtSOTpJlVwLGXce02I4k4GvgshhMgTS6Uxy20rofGKaX8eNXc4QhQaJbIny1RrFyotUd/rkRWLpJzd3mwo/ddfTKO8qxujX3w5X/I3Nn2/Hk/z1Zx51KnrlS/5F/b0SrPAPt35UhLXLnTWlH6sT3biE1P+Z/r0H/7vG8qXL89r40bmS/7mTv949+f4+dsvaFC3dr7FY5EE91KdV7k9fzwc7tDT0pc1/k0ZevkezauVz1E+QhRH0pOVD0JDQqhZyYVGtarQoIYnbVo0ZvEvP+U6v/v37/H36lUMGT5S/9p/l4NpUMOTO7cffWGvW/sXbZvV58b163mKPyvjXn6Nb7/+Mt/yF0XXw5BQbCp64VyzOWWrN6Vq0/Z8t/D3fDve3XsPWLZ6HWOHPZ9vx0hx6fIVylZvys3bd/SvLV+zgapN23P1+s18O+6br4zhs5nf5Vv+pvZm2f04alG8vfQICYk5W8ZHiOKoRPZkmWrtQk1ZGhy34+d/jnIu5Tly9hIAuzav4+UxI+jS+zkqVPIwOv+UvFf/+ScdnuqKtZ0jKT9KPWvUofMz3Vn8y8+8/+HHeJ84ztQP32fxyr9xr1xNny4rt+8/pLyrcbfhp8TTuduzfPz+29y69wBXt6zvSjJ2DFRRSq8py0zPk5K0dqGl0jLcrXb23Dlcy5fjjv9BANZt2Unfka/Rv1dXPCtVwNrIu+FSelyyS7/sr7X0eLoDjvaljEqf0/xTp69bswrPdu3E/F9+56tPJnP4+CkmfTidbX/+Qo0qFQAtTXpD+d++cw93N0O9O5nH0797Z157byr3792hgrtrruPPMr2F4b99bs4fd2L4lB289bAPP233Y0K3xjnKQ4jiRnqy8sG5s6dp1LS5/nnT5i0AiIuPy1V++3Ztp3XbJzK8Pv61Saxc8iuBAed4ddxIps38lsbNWmSZV1hoCMt+X0y/Hk/z3sSMlx4B7t+9y8svPE/bhjVpVrMSL74wlPDwMABsS5WiQeNm7N+9I1dlEcWXj985WjVtqH/eprnuCzYuPvtL6rmxZdd+OrZtlea1O3fv03vEq7jXb4dD9RY8O/wVwsIjsswnJDSMv9ZvZcwbH/Lc2EmZpnvv9XH89Mef+J0Lov/oifw0ayqtmjXKNu8ffl1B666DGDVxisE070//mgGjJuifvzP1a54aMJq4uDhKlbKlReP6/Lv7YJbHKUz6lQ2klWUw3+29zPWHUeYORwizkkZWPvA/c1rf2AkLDeGbGf+jYeMmVK5SLVf5BZ47S/WatTO83qBxUxo3a86Ans8w7IXR9OjT3+D+SUlJHNizkzdfHk2nlg3Zv3c3r058ix//WGUwfUREGCPGvsTek/7sPuHHw4f3WbHk0WWfmrXrEODvl6uyiOLr1JlztG6ma1iFhIbx4Zff0qJJA6pX8cyX4505F0TdWtXTvBYWHsHrY4dx5dQu/vPeyb37D/np94znue/ZAL6cu5AOvYdTtflT/L7qH1o2acCsqe9kerzmjRvQulkj2nQfzCujBjOoT3eD6ZKSkti+5wBDXnqbqs2fYtueg3w46SXWL5lvMP07E8ax99BRTp3x58ffVrJ1937W/vodNjY2ANSrUxPfswHGVovZKQXfVNhGUpLGO38ckAWkRYlWIi8X5reAs2fYvnkDf/zyA45lnWj7eDt+/n05SimOHNhHlWrVqeRZ2eC+Sxb9RMeOHaleoyZffPI+fZ8bQnhYKHb2DhnSJiUlYWFhiYWFYvwEwzM2L1n0E78smItzORf6Pz+Uj/83C1dnRyDzy2dVq9ekanXdXDc2trY80b4TYaEh+u129g7cvZ2zwdui+PPxC2Dtph3M/XkJzmUd6dyuNRuWLkApxfvTZvJ8n660bm785aPp3yzgsRZN6NIpYy8uQEhoOA72dmleq1WjKrVqVAXA1taGLh0f52FoWJo0/Ua+wqZtu2ndvDGfTn6VJ9s9RqlS2U/NkZSUhKWlJRZK8d7rLxpMM2/RMmZ+vwiXcs6MHtKX77/4iPIuzlnm61LOmYnjRzLytQ8IDQ/nwIZllHV89H53sLfj5u272cZXmFSxCeUVh4N8d7Mj/xy/RL/WMneWKJmkJ8vE4mJjuXg+kK0HTnA84Ao7j57my2/m4uauG3+zeuWSLGdEHjH2JarX0H0gXTwfSK06XjiWdSIyIjxD2hlTpxAeFkq16jVZv3a1wfyuXfmPsJAQ6jVohFf9RjiXK5dtGbas/5vBz3bhsQY1aFGnMj/Nn0u1Go8+JCMjwnEoK8tniEdiY+M4d/4SAYc28SDoCBePb+OXOZ9T0V035u9swHka1K2Vozz9gy7SuH7dTLc7OzkSHhGZ5rW/1m/liZ5Dcav/BE61WjPj+1+oU7NamjQfvzWB9994mfj4eEa+/gFjJn3Ikj/Xcefu/SzjefvTrwgJDaN2jaosW7PBYJrg/67xMDSUJg29aNKgLi7lnIwqa9OG9TlzLogvP3yLyh4V02wLj4jEyTHjj6zCbmL5Y9RUd/hk3VkeROZuqIQQRZ00skwsKMCf0mXsDA5w3/nvZvZs/5d3XhvPP3+t4NuvpvPCgF48++TjnDx2BIDh/Xro08fGxGJja0vd+g25fOlCmrxW/rGY7Zs38sNvKxk/YSI///C9wW75Dz77gh1HfKntVZ/pH77Dk60aMXvmFwRfumgw/sMH9vL1558yZdoMDvgGcdQ/GJfy5anf4NHYk4vng/Cq39Dg/qJk8gs4j12Z0nhWMjyY/8btO7z5yQyaP9WfBYuXAzB/0TIe6/48bbo9z9GTvgD88OsK2vV8njET3+POvQe4u5Vn8469TJ05D4AHD0PoOfQlABrXq0PQxcv6Y+zaf4T3pn3Dt59/wI3Te7kXcAi38uVo2iDtdCPNmzRk2vuT8N65ljN71tGl4+Ns2LaHeu168sKE9wzG//MfK/l78w7W/TGf914fx9fzFxt8v30z7T0CjmynQd3avP7B/6je8mk+njGX85cuZ8w02Zlzgbz+/lRGPt+XxcvXZNh+LugiTRpkPmVKYWWlNL6rsJGoBHh3yQFzhyOEWUgjy8T8/U5TO5M5pDp36UaDxk1Z+vdm+j43hPGvvcUfazbyv2++Z9um9Ty4fx+X8ro7iO7du4tL8p1/HZ96hmOHHw18PbhvN7O//IyflvyJi6sr3Xv1Jj4+nh1bNxk8rourK2Nefo0Nuw/z/aIlhIWG8tyz3fhg0qsZ0gacPUPFSh54NWhEWEgIH0x6lfv37lGrjq5HITYmhrOnT/FExyfzVE+ieDl1xj/Tnqp79x8SEhLGVx+/zaFNK/h56WqOeZ/mwDFvDm9eyZ8/z+br+Ys54x/E3sPH2b9xJeNfGIxFco9vyyYNOeV3DoAv5/7M+xN1l+p6PN2RvYeP64/jezaAyh4VaNKgLg9DwhjzxofcufeA+nUzv1RVwd2V0UP68+cvc7jjf5AP38x4M8jOfYf4eMYcNi79ATdXFwY+25W4uHjWbdlpME+38i5Menk0p/euY83i7wgJDadtjyGMeePDDGmv37xN/xdeYd5Xn7Hgq084c+48ew4e02+PiYnl5Gl/unR6PNMyFGYNSt1jjN0RdlyOZvOpK+YOR4gCJ40sEzvnd5raXvUNbvsv+CLVksc63b97l0/em8SI/j35aPJE3CtWIvCcH3Xq6fYNPOdP3eT/931uCHt3biMmOpqL54N48+XRzPx+oT6tpaUlY8a/ws/z5mQbX8MmzZj6v684dOosg18Yk2F77/7Pk5AQT2uvqowf8RzVatSkVu26+kG4u7ZtofXj7XGvUDHDvqLk8vELoKFXxpszQNdTM7h/L5ydylKqlC1lSpfiny07eXX0EJRSWFhYULqULf9s2cH4EYNQSqGUolG9OgC4ubpw/0EIV6/f5PLV67R/rCUALwzqw+Yd+4iOjgFg2IBniY9PoFydtvQa/gq1a1Slfp2a+nM3xdP9R2BT0QvlVi/Nw6piQ8ZO+ihN2oDzlxj+8lv8+v1MGibHY2lpyVsvj+Sr73/Jtl5aNGnA919+xI3Te3l5ZNr5vMLCI+gx5CXeeGkUz3Z9kjJlSvPOhDF8+MW3+jQbtu2m0+OtqFTB+KlWCpt33A5TRd1nyhofo5YhE6I4kYHvJvbpl99kuu18YADVa+l+7c+f8xW9+w+iw5NdeP+NV6jXsBGB/n7UqdcAgKCAc/r/l3Nxoe9zQ1i5ZDGjxk/gmP/lDHmPGD2OoWNfMTpOGxsbmjRvmeF1F1dXVm1MOz3DxElv6f+/6Ifv+GK24bukRMk1b8bHmW474x+IpaVuVvA1G7bxbNfO3H8QQlyc7gt33qJl9On+FAeOniQ+PiH5tSV06/xowLu7qwuvf/A/pr//6AaP8i7OvDCoDz8vWcXE8SNxc3Xh0OYVaY798dsZe2t3rF0CGDfPlFftGtz0P5L87FH6CWOHMWHssGz3T2FjY5Nh0L+jgz2+e/4hPtXM+ZMnjGHyhEc/fmYt+JVFc6YbfZzCyEYlMdd9IwNujmDKsoPMH9fJ3CEJUWCkkVWAatSszdyvPuf61Ss0b/0YX332EZvXreV84Dm86jdkw5o/ebKrbkxW4Dl/OnXrpd/37SmfmivsNFZv2W3uEEQR4xcQRJnSpRg07k3s7crw49efcuXaTUa+/gEfz/iODm1bMvDZrlSv4smICe/hWfE3bt+7z7sTRuvzaNGkAYEXgmlUv06avL/48M00jZTi5uhWw9OsFDXNSt9meJnj/HGhDf1OX+XpxobvrhaiuCmRjSxzrV1Yt3ZNtu5+NAC0T+8+adJ99c23+vQzZ80mycIKCmhtRElvXHpZuzDnaxfO//p/AFinOkzVqtXYs/5Rr1N8IjRu2BDfvZv0a+1ZWz5ad8/nbBBffPi2wQZVUVqLsDCmN9XahdkZV/oQO6Jq886f3uyu5Y5TGZvsdxKiiJMxWUKIQuvWnbv0G/kK7dq0oGpl45ekEoWPrUrkc8dNhMZb8PbvcrehKBlKZE9Wfq9dmJmisDafpM86vaxdqGNo7cLM5X6tvcoVy7Nx6YLk1zPb38Rr+ZW09CZcuzDbNNxkePQh/vivHauPXGDgYzmbO02IokZ6soQQQhSYlxyPUlO7xSfr/bkRImsbiuJNGllCCCEKjJXS+NJpIwmJSUz4ZR9JSbK2oSi+pJElhBCiQFW1CuHdsrs5dS+RH7edMXc4QuQbaWQJIYQocGPL+dDa8hJz9vyH/7WH5g5HiHwhjaxCZtb/pvLrzz+aO4xCadb/pvLbQtNOhDqgWyfOB5xL81q/Hk8TFBhg0uMIIdJSCuZV2kwZYhi/6ABRcTmbakWIokAaWYXIg3v3+OevFQwZPlL/2pJFP9H/mY40qFKe9yZmXFctL2Kio5k78390aduUpjUq0qlFAyaOewGfk8ey37mApdTN4BEZlwLKi7GvTGRu8jxOKca9/Brffv2lSY8jhMjIzSqKb8uv53qU4t0lh8wdjhAmJ42sQmTtqmV0fOoZSpUurX/NrUJFXnnzHQYOHmHSY0VFRTLq+T6cPe3DD7+txOfSTTbtO8rjHTtzYM8ukx7LFAzVjSk81bUHRw/u4+6d2/rXnn6mG0cOHUjzmiiaPvh8Nt/+9LtRaVt3HcTZgPP5HJFpfPi/b4wu1+Pdn+NsYOEtV2eHK4wsc4SN58P583DhjVOI3DCqkaWU6qaUClRKXVBKvW9gu61SalXy9qNKqWrJr1dTSkUrpXySH3IdLAv7dm2nddsn0rzWtWdvunTvhVM5Z5Mea+b/poGmseC3FdSq6wWAnZ09g0eM5rW3M/yJzc5Q3Rjr3r27jB85jLYNa9KsZiVeGjGIiPAwAGxLlaJB42bs3/1ovUbbUqVo2KhJmtdE/nkYEopyq4d9tRbYV2uBc83mVGvWAU3L211nd+894I8/1/HSC89nnxiY/OpoPvnq+zwdM7X05bKv1gKPxh1NUq5lq40v15uvjOGzmd/l6Zj57SP3A9RTN/h4/TmC74abOxwhTCbbRpZSyhKYD3QH6gNDlFL10yUbCzzUNK0WMAf4KtW2i5qmNU1+mPZ6VyE3fvhztKhTmaZe1WnqVZ0WdSrrH+OHP5chfeC5s1SvWTvf4wp5+JCVS3/n7Q8/w8qqaMxHm5e6iQgP54Ux49h70p/dJ/x4+OA+K//4Vb+9Zu06BPj7pdnH0Gsif/j4BVDOuSwRl08ScfkkDy96c/nUPpTK25qEv638mx5Pd6B06VJGpe/d9Ul2HzzGrdt383TcFOnLFXH5JNdP781zuf7482+6PdnR6HI9+8yT7D1kunLlByul8bPHOiy1BMb9tJfYhMTsdxKiCDDmG7Y1cEHTtEsASqmVQB/AP1WaPsDU5P+vBuapvH6S5KOCWrtw0R/LgUdr4SVZpKvudPuGh4XiWKaUwbXzlJaEIsng8XK6Nt+R/Xtwcnbmsdat9DH06/E0wZcuEhcXx2/L/6L1Y48DsHjhD8yY/imVK1ehz4BBTHz73Wzzz8tagV9O+5Qez/amSbMWadKk1I1lUnyO869RtQo1qlYhycKC0mUdaNe+I+Eh9/V1aW9Xhru3b+ufWyQlYG9Xhjt37maob1m7MOdrF2a3dt7J0wHUq11Lvy6hqdbm27xzPyMHD0iz3uH7078m6OJlVi2aDyjenzaTU2f82bBsITY2pWjeuAGbdh3khUH9jI4/s3jSlyun8Wdmy859vPB82nLduXefl97+iOPep4mKjqbD4635Y/4sHB3ssbS2pWmjjOVKUVBrF2aX3ppIPiq9mSkR/fhw+WFmvdAuR/kKURgZ08jyAK6men4NaJNZGk3TEpRSoYBL8rbqSqlTQBjwkaZp+/MWcvHlWNaJyMiIXO07eGA/jh4xPHC0Ras2/Llus/753bt3cHNzT5Pm7807uHfvLm0ae1Gnbj396/8FX+Kvtetp1qJlxkZiPggKPMcbkzM25AzVzdABvTl6+KDBfFKXedPG9fz6y0KCg4OJj48jJjqaz2fO1qeNjIjAsWzZNPsbek0Yp++Ilzh4zNvgtidaN+efJT+lec3H7xz16tQ0eRx+54KoW6t6mtfemTCOuo91wcfPn1M+p/l39372rFuOjY1usWKv2jU4c9Y0d5bmV7n8A4KoUzNtucLCI5gwZjjtF7YiMiqKPiNe4uclq3j71bEAeNUyXbny01NlLtIv7gSr/VvS8mAQg5+oY+6QhMiT/P7WvAlU0TTtvlKqBfCPUqqBpmlhqRMppcYD4wGqVKmSzyEV3NqFY4f05+TRwwa3tWjTlkUr1qZ5rW79hlwMvkzjps2BtGvnacoCDYtMj7d87YYsY0nd+V7Rswo3b9wgTrPA0vLRr0r/cwFU9PDEwcVNn97P7wyTp3xMkoVVntf++/ar6XgfO8rDhw+YOmM2LVo/liH97du3mT71E/x8T/Hc0BcYNvpF4FHdNGj+qH2/5O8tWR8fOHxgL1998Tnf/fgLdRvresc6t2pI3UbN9PFduHCB3gMGpYn3woULPDtwcIYyyNqFOlmtXbhpedqhlyk9Lo/Wzku73+mz57hw+Qp/rdf9PT0rVeDU7g15XpsvJCwcZ8cyaV6v4OrExPEjeenN9wkLD+fAhmWUd7bXx1TWwY6bt++myyt3awWmL1dljwqc2bve6PgzExIWjpOjXZr09WpVoV4t3WenfZmyPNPxccLCQpPTKBwd7LhzN325khXg2oXGpJ/psIegKxX5eEMSDauUp2HlcjnKV4jCxJiB79eByqmeeya/ZjCNUsoKKAvc1zQtVtO0+wCapp0ELgIZfppomrZQ07SWmqa1dHV1zXkpCqlFK9bic+kmZy5c4cyFK/hcuql/pG9gAXR86hmOpeuZSUhIIDYmhsTEJBKTEomNiSEhIW/zyTzRviN29vb87+P3CHn4AE3TuPrfZVavWIJX/YZp0jZp3pLHmzVk+idT8nRMgPGvvcUfazbyv2++Z9umjF82D+7fJyw0hHc++oxVG7bz57JHd08ZqhtjBJw9Q8VKHnjVb0hYSAgfTHqVB/fuUquObrB/bEwMZ0+f4omOT+r3iY2Jwe+Mb5rXRP6IjY3j3PlL7Fu3hJALxwi5cIxTuzeYJG9nJ0fCIyIzvN60YX38zgXx5YdvUdmjYppt4RGRODk6ZNjn6f4jUG71DD7a9RpmVLkMN7ByUa6yjkREpi3XX+u38kTPobjVfwKnWq2Z8f0v1KlZTb89ItJwuQoja5XErx5/Y0c0Y385QGhU5kMyhCjsjGlkHQdqK6WqK6VsgMFA+k+L9UDK5E4DgV2apmlKKdfkgfMopWoAtYFLpgm9+On73BD27txGTHS0/rUFc2bSqJobC7+fzfrVq2hUzY0Fc2bm6Tily5Th95VruHXzBt3ataRNg+pMGDOMylWr8+mX3+jT/bnsd0qVLsMJvyA+nvZFpvndvXOb4f166B9DB/Rm6IDeaaZAuH/3Lp+8N4kR/Xvy0eSJuFeslCGfwHN+PNt/EGWdnLEtVYrSqaZrMFQ3xujd/3ni4+NpUb8W40c8R7UaNalZx0t/eWjXti20frw97hUefdnu3P4vbdo+keY1Ybzug8enuaPOuWZznGs2x75aC7oPHp8mrV/AeSwsFI3rZ7wslJSUxFMDRuufd+r7QqbH/GbBIh7vMQT3+u0oU7UZP/y6gsb16hB08XKadGf8g3j9/amMGNSXxcvXZMjnXNBFmjTwyvD6jrVL0O6cM/g4sHFZhvRnA4KMLtfT/Q1Pz/L1vIxlAmhUv26acu3af4T3pn3Dt59/wI3Te7kXcAi38uVomqocAecvGSxXYeVmFcVPbn9zLwbGL9wt6xuKIivby4XJY6xeA/4FLIHFmqadVUpNA05omrYeWAQsUUpdAB6ga4gBdACmKaXigSTgZU3THuRHQYqDci4u9H1uCCuW/s7oF3U3Yk58ZwoT38l7L1J6VatVZ8Gvy7NMs3fnNt756DNsSz26i+nWjetUqOSRJp2rmztL/3405itlsHjqS23z53xF7/6D6PBkF95/4xXqNWykz8+jghsAgf5+WFrq2v3/blxH52e66/dPqZuVSxYzdlzaL+qsuLi6snrD1jTxTHjrPf32RT98xxez084i/8uP8/jym8J9y3thtmXlwjTPM14ufOTUGX8aetXWN3pTu3T5KjWr6TrRHzwMwbmso8HjHT91mqMnfTi4aTlnA84zffYPvDJ6CNExsew9fJxhA58F4PrN2zw74hXmffUZT7VvS502T7Pn4DE6PdEagJiYWE6e9uf3eTNyX/hkPn7njC6Xk4Gxf8e8T3PkpG+GMgF0e7ID+w8fY+SgXgD4ng2gskcFmjSoy8OQMN7+9Cvu3HtA/bo19eXyPn2WP+YVrQl229jd4L2yu/jiVhdm/H2CKQNamTskIXLMqDFZmqZtBjane+2TVP+PATLMSaBp2hog489Fkam3p3ya5R2LBen54aN465WxODjY8/KEN2jV/knefGUMK9b9m+O8mrd+jK8++4jN69ZyPvAcXvUbkpCQwJuvjOHPvzcCEHTOn1KlS/PGiyMpY2fHZzO/TZPH21M+1f3HhPWzesvuDK+t3bQdSDuOTeQPH78AWjZpaHDbmXNBNK5fF4DT/kE0qmd4EPSm7XsYM/Q5lFJYWlpibaVrTL8wqA9Nn+xHdHQM8QkJ9BjyEm+9PIpnu+ouA78zYQwffvEtBzfpfmxs2LabTo+3olJyoz8vfM+eM7pcDQ2Ua+O2PYwbPjBDmQCGP9eXVl36Eh0dQ+nSpRg24Fn+Wv8v5eq0pYFXLXp16Uj9OjX1DbyN23fToW1rk5SroL1Y7hRHwiux8Dg0rHKJ3q1qmDskIXKkaEySJMyiw5Nd6PBkF877+bDkt0WUKVuOvs8NyVVevfoOpFffgWle8/U+kSa/L+aYdl1CUfjNm/FxptvOnAvSjytau2k7Hdq2BODajVt4Vnp000BMbCxx8bqG94+/r2JAry4AlHdx5oVBffjpj1VMemkkvnv+AR5NmTB5whgmT3i0TNOsBb+yaM50k5Rr7hefZDqQPX25Hm+TsVwxsbHExWUsU0q5hg18VC43VxcObV6R5hgfv/2q/v9zfljMT7M/N0m5CppS8Knzdv574MLkNYlUcy9L4you2e8oRCEhy+qIbC34bg7PDR5Kk+YteX74KJPla+r8RPFy5tx5/t68g/6jXue/azc45n2ahIQEhrw0OU26UYMH8MWcBXTsM4JyTmXp1/NRg+SLD99k0ksj02dt0NGtqwz2Kpla+nIdP3UmQ7nGDB3A9Nk/GCwTwOdT3jK6XAc3/0lDr6I7FUIplcB3TmuxJ5pRPx/kVmjOxmUKYU7SkyWy9d2PvwBy+UwUrAchoexc82ua1455n+aFQb3TvOZVuwZH/l2TgykWzCt9ueITFcdPpS2XV+0anNi+2hzhFUrlLaP4vcJfDLw5nOHzdrHhnW6UtsnZJKhCmIP0ZAkhCp2YmFhsbTLOy9a6eWNeHDHIDBGZRmblatWsaJerIDQqfZfZLhu4GKbx8s975I5DUSRII0sIUeiUKmXL5hULs09YxBTXchWUnmUvMMlhD3uvxjDtr6PmDkeIbJXIy4UFtXZhirys5SfpC1d6WbvQ9GsXSnrzpi8saxcam/55m5P4W7jw26nGVHQ+w0vPNMpRfkIUpBLZyBJCCFE0KQWfOG3n/gM7ZuzScHOyo19rmdpBFE4lspFVUGsXppfXtf8kvfnTy9qFOlmtXZhR7tb+k/QFlL6QrV1obPqlDuvod2Uw76zVKO9YivZeGVeREMLcSmQjy1TuhMfi5mCbbTql6briNWVcN7mkL7zp74THGpWHECJ/lbJIYLnnap69OpwXfz/BX68+QaPKMoeWKFykkZUHH6w9Y1S6lDE8xvZcSPqilV4IYR7OVjH86bmSXtdeYPiPB/hn0pNUdy0aC2GLkkEaWUKIHEsMu4Olk3HLtFikXOUx8l5mSV/A6SNuZ7KhaKhkHcHyiqsYeHMYA//f3r1HR13eeRx/f2dymSSTALmHkAQIkYKiKEEQLW7rDS9HxBti1+qurT271LWr53jWY3fXdc962uOdyta6K7VVK3W9gW0txBu4VBbQgiaogHJLIFwSICRhcn32j4wVFSSJM/zm8nmdM4ffzHzz48vDJR+e+f2eZ97rvPSjcynLy/K6LRFAIUtEBmHfwlv6fU3Np3euqT4+6uPRNwJNPF28kDmNc7ji4dd46dZzGT400+u2RLROloiIxL9TMnbz66LfcqCjhysefJVd2n5HYoBCloiIJITqzEYezH6OplAvVzxUw56DIa9bkiSnkCUiIgljUvoO7gu+wK52x2X3L2Xnfs1oiXcUskREJKFMCWznicJn2Xuol8serGF7U6vXLUmSUsgSEZGEc2awnqeKF3Iw1M2seW+yeW+b1y1JEkrKuwsjtXeh6lV/vOpjbe9C1as+HurL2cW/+/+Hu7rncPUv3uapG6cwtljraMnxo5ksERFJWKNsD3eW7cBncOWjf2LlJ01etyRJJClnsiK1d6HqVX+86mNt70LVqz5e6ts6jbLMXp6fPY3rF6ziu4+v4qFrJnLRhJIBnU9kMDSTJSIiCW/EsEye/7tpTBgxhLm/eZcnVmz2uiVJAgpZIiKSFIZmpvH096Zw3rgi7np5PXe/vJ7unsRdCV+8p5AlIiJJI5Dq5+d/PYkbpo1kwYrN3PirNbSEurxuSxKUQpaIiCQVv8+469ITuWfWBFZs2sus+Su0xINEhUKWiIgkpWunlPPU96bQ3NbJZfNX8NbGPV63JAlGIUtERJLW1NF5LJp7FsU5Ab67YBUP1mygp9d53ZYkCIUsERFJauV5mbw4dxqzTi3l4dc2cv2CVext7fC6LUkAClkiIpL0MtNSuP+qU/jpFRNYvaWZi+e9xarNzV63JXFOIUtERAQwM2ZPLufFvz+TjFQ/c/5rJQ/UbKBLyzzIICXliu/au1D18VavvQtVr/rB1e9v99Hb0zOgc4wfnsPim8/irsV1zHttI29+tJsHrp7ImMLggM4jopksERGRL8gJpPLA1RP5z++cxrbmdi6e9xa/+tMWenVRvAxAUs5kae9C1cdbvfYuVL3qB7934SH/wGbDDnfRhBKqK4Zx+/Pv8a+L61hS18h/zJrAqPysQZ9TkodmskRERL5CYU6AX94wmXtmTeD9+gNc8NBy5r+xic5uXaslX00hS0RE5BjMjGunlPPabWdz3rgi7l3yEZf87C3e2ao7EOXoFLJERET6qTAnwPzvnMbj11fTGurmip+/zW3PrmNXS8jr1iQGKWSJiIgM0Dnjiqi59Wx+cPZoXl63g2/d9yaPvL6RUNfA7mSUxKaQJSIiMghZ6SncceE4am6dzvSqAu5buoFz7l/GorUNugtRAIUsERGRr6UiL4tHr5vEM9+fypCMVG5ZuJaL5r3FkrpGnFPYSmYKWSIiIhFwRmUev7v5LObNOZXO7l5+8OQ7zJy/gjc/2q2wlaQUskRERCLE5zMuPWU4S/9xOvdeeTLNbZ3c8MvVzJy/gt+/t5MefYyYVJJyMVIREZFoSvH7uKq6jJkTS3nunXoeW/4xc3/zLhV5mXz/m6O5ctIIAqmDXyRV4kNShiztXaj6eKvX3oWqV/3g6gezd2EkpaX4uHZKObMnl7G0rpFHl33Mj1+q5YGaDcyeXMa1p5dTlpvpWX8SXUkZskRERI4nv8+4cEIJM04q5u1PmnhixRZ+sexjHl32Md8aW8h1UyuYfkIBfp953apEUFKGLO1dqPp4q9fehapXvTd7F0aamTGtMp9plfns2H+Ihau28czq7fzNE6spHZrB5aeVMuvUUkYXBL1uVSIgKUOWiIiI14YPzeDW88dy8zlVLK3bxcLV25j/xiZ+9vomJpYN5fLTSrnk5OHkZqV53aoMkkKWiIiIh1L9Pi4+uYSLTy5hV0uIRWsbeOHdBv5lUR13v7yeMyrzuODEYs4fX0RhTsDrdmUAFLJERERiRFFOgJumV3LT9Eo+2NnC4nU7WFLbyI9fquWfF9VyWvkwLjixiL8aW0hVYRAzXcMVyxSyREREYtC4khzGleRw+wVj2bS7lT/WNvLHukbu+cOH3POHDykZEuCbVflMP6GAs8bkMzRTHyvGGoUsERGRGGZmVBVlU1WUzc3nVNGw/xDLN+xh+YY9vFLbyLNr6vEZnFQ6hMkjc5k8MpfqkcPID6Z73XrSU8gSERGJI6VDM5hzejlzTi+nu6eXdfX7WbZhLys/aeLJlVt5/H83AzC6IIvJFblMqhjGSaVDqCoKkurXRi/Hk0KWiIhInErx+5hUkcukilwAOrp7qG04wKrN+1izpZlXanfy2zXbAUhP8TGuJIcJpUOYUDqE8cNzGFMY1MrzUaSQJSIikiDSU/yHha5Kensdm5vaqG04wPv1B3i/4QAv/rmBJ1duBcBnUJ6byZjCbKqKgpxQFKSqMJtR+VlkpSsifF0aQRERkQTl8xmVBUEqC4LMnFgKQG+vY0tTG+t3trBxVyubdreycfdBlm3YTVfPZxtY5wfTKMvNpDw3k4rcTMpyM6nIy2L40ACF2QHSUvTR47EkZcjS3oWqj7d67V2oetUPrt7rvQtjkc9njC4IfmlV+a6eXrY2tbFxVytbmtrZ1tzGtuZ23tm6j5fX7aDXff48+cE0inICFOcEKBrS92NxToCCnHTystIYlplGXjCNjFR/0i410a+QZWYzgIcBP/DfzrmffOH9dODXwCSgCZjtnNsSfu8O4EagB/gH59ySiHUvIiIiEZHq9zGmMJsxhdlfeq+rp5eGfYfY2tzOzv2HaGwJsaslROOBEA37D/Hutn3sa+864nnTU3x9oSsrjdzwY1hmGtmBFLIDKQTTU/9ynB34/HFWWnwHtGOGLDPzA/OB84B6YLWZLXbOrT+s7EZgn3NujJldA/wUmG1m44FrgBOB4cCrZnaCc87T/1Zo70LVx1u99i5UveoTY+/CeJXq9zEyP4uR+VlHrQl19bC7pYM9rR00t3Wyr62TprZO9rV30tz22WNrUzv72ztp7ej+0uzYF5lBIMVPRpqfjFQ/gVTfYcd9j4xPH2l+0lN9pKf4SU/xkZ7iIz+YzmWnlkZ4NPqvPzNZpwObnHOfAJjZQmAmcHjImgncFT5+DnjE+qLnTGChc64D2Gxmm8Lnezsy7Q+O6+2hZM/7ZIf2ApAWagFgWGtDv75e9ao/3vVDWupxPvhzw7G31DgY6vtxx8F+nV71qk/o+h0HfOSUDiygyeAEUv2U52VSnpfZr3rnHG2dPRwMddEa6qYl1N133NHNwU+PQ92Euns51NnDoa6+Ryh83NrRzZ6DHYTCrx/q7CHU1Utnz2e/31WFwZgPWaXA9sOe1wNTjlbjnOs2swNAXvj1lV/42i/9as3sJuAmgPLy8v72PmhTz7+AMd/YEvWfRyRSnJtCT3cXKanHXtE5d4DnVr3qE7m+AigcVTnAs8jxYGYE01MIpqfAkMidt7fX0dnTF7Z6jzVVFmUxceG7c+4x4DGA6urqqI9IcWUVxZVV0f5pRERE5Djz+YyAzx8T63/15/7LBqDssOcjwq8dscbMUujLpE39/FoRERGRhNOfkLUaqDKzUWaWRt+F7Iu/ULMYuD58fCXwunPOhV+/xszSzWwUUAWsikzrIiIiIrHrmB8Xhq+x+iGwhL4lHBY45+rM7G5gjXNuMfA48GT4wvZm+oIY4bpn6btIvhuY6/WdhSIiIiLHg/VNOMWO6upqt2bNGq/bEBERETkmM3vHOVd9pPe0Jr6IiIhIFChkiYiIiESBQpaIiIhIFChkiYiIiESBQpaIiIhIFChkiYiIiESBQpaIiIhIFChkiYiIiESBQpaIiIhIFMTciu9mtgfY6nUfMSQf2Ot1EwlGYxpZGs/I05hGlsYz8jSmn6lwzhUc6Y2YC1nyeWa25mjL9cvgaEwjS+MZeRrTyNJ4Rp7GtH/0caGIiIhIFChkiYiIiESBQlbse8zrBhKQxjSyNJ6RpzGNLI1n5GlM+0HXZImIiIhEgWayRERERKJAISsOmNm9Zvahmb1nZi+a2VCve4pnZnaVmdWZWa+Z6e6Yr8HMZpjZR2a2ycz+yet+4p2ZLTCz3WZW63UvicDMyszsDTNbH/47f4vXPcUzMwuY2SozWxcez3/zuqdYp5AVH2qAk5xzJwMbgDs87ife1QKXA8u9biSemZkfmA9cCIwH5pjZeG+7intPADO8biKBdAO3OefGA1OBufoz+rV0AN92zp0CTARmmNlUb1uKbQpZccA5t9Q51x1+uhIY4WU/8c4594Fz7iOv+0gApwObnHOfOOc6gYXATI97imvOueVAs9d9JArn3E7n3Lvh44PAB0Cpt13FL9enNfw0NfzQhd1fQSEr/vwt8IrXTYjQ981q+2HP69E3MIlRZjYSOBX4P49biWtm5jeztcBuoMY5p/H8CileNyB9zOxVoPgIb93pnFsUrrmTvunvp49nb/GoP+MpIsnBzILA88CPnHMtXvcTz5xzPcDE8LXBL5rZSc45XUN4FApZMcI5d+5XvW9mNwCXAOc4rbtxTMcaT4mIBqDssOcjwq+JxAwzS6UvYD3tnHvB634ShXNuv5m9Qd81hApZR6GPC+OAmc0Abgcudc61e92PSNhqoMrMRplZGnANsNjjnkT+wswMeBz4wDn3gNf9xDszK/j07nYzywDOAz70tKkYp5AVHx4BsoEaM1trZo963VA8M7NZZlYPnAH83syWeN1TPArfjPFDYAl9FxQ/65yr87ar+GZmzwBvA2PNrN7MbvS6pzh3JnAd8O3wv51rzewir5uKYyXAG2b2Hn3/yapxzv3O455imlZ8FxEREYkCzWSJiIiIRIFCloiIiEgUKGSJiIiIRIFCloiIiEgUKGSJiIiIRIFCloiIiEgUKGSJiIiIRIFCloiIiEgU/D8Nthhc3yemvQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import phat as ph\n",
    "\n",
    "shape, mean, sig = 1/4, 0, 1\n",
    "x = np.linspace(-2.5, 3.5, 1000)\n",
    "\n",
    "dist = ph.dists.CarBenHybrid(-shape, mean, sig)\n",
    "\n",
    "fig, ax1 = plt.subplots(1,1,figsize=(10, 6))\n",
    "\n",
    "tailprops = dict(ls='--', lw=2, alpha=.5)\n",
    "ax1.plot(x, dist.pdf(x), label='CarBen')\n",
    "\n",
    "paramtxt = r'$\\mu, \\sigma$ = ' + f'{dist.mu:.0f}, {dist.sig:.0f}'\n",
    "paramtxt += '\\n'\n",
    "paramtxt += r'$\\epsilon$ = ' + f'{dist.xi:.2f}'\n",
    "paramtxt += '\\n'\n",
    "paramtxt += r'$a$ = ' + f'{dist.a:.2f}'\n",
    "paramtxt += '\\n'\n",
    "paramtxt += r'$b$ = ' + f'{1 / dist.b:.2f}'\n",
    "paramtxt += '\\n'\n",
    "paramtxt += r'$\\alpha$ = ' + f'{1 / dist.xi:.2f}'\n",
    "paramtxt += '\\n'\n",
    "paramtxt += r'$\\gamma = $' + f'{dist.gamma:.2f}'\n",
    "ax1.text(\n",
    "    .02,.69, paramtxt,\n",
    "    transform=ax1.transAxes,\n",
    "    bbox=dict(boxstyle='round', ec='.8', fc='w')\n",
    ")\n",
    "\n",
    "txt = 'Junction Point\\n' + r'$x = a$' \n",
    "arrowprops = dict(\n",
    "    arrowstyle=\"-|>\", connectionstyle=\"arc3,rad=.4\", fc='black', ec='black'\n",
    ")\n",
    "ax1.annotate(\n",
    "    txt, xy=(dist.a, dist.pdf(-dist.a)), xytext=(dist.a*4, dist.pdf(-dist.a)),\n",
    "    arrowprops=arrowprops\n",
    ")\n",
    "plt.rcParams['hatch.linewidth'] = .5\n",
    "txt = r'$P_{tail}(X<a)$'\n",
    "txt += '\\n'\n",
    "txt += r'$= (1 - G_{\\xi,-a,b}(-a))$'\n",
    "ax1.text(\n",
    "    dist.a*5, dist.pdf(dist.a)/3, txt, fontsize=12, ha='left', va='top',\n",
    "    bbox=dict(boxstyle='square', ec='w', fc='w', alpha=.8)\n",
    ")\n",
    "fill1 = ax1.fill_between(\n",
    "    x[x<dist.a], dist.pdf(x[x<dist.a]), \n",
    "    alpha=.8, hatch='+', ec='C5'\n",
    ")\n",
    "\n",
    "x_max = 2\n",
    "txt = r'$P_{body}(a\\leq X<x)$'\n",
    "txt += '\\n'\n",
    "txt += r'$= F_{\\mu,\\sigma}(x) - F_{\\mu,\\sigma}(a)$'\n",
    "ax1.text(\n",
    "    0, dist.pdf(dist.a)/3, txt, fontsize=12, ha='left', va='top',\n",
    "    bbox=dict(boxstyle='square', ec='w', fc='w', alpha=.8)\n",
    ")\n",
    "fill2 = ax1.fill_between(\n",
    "    x[(x>=dist.a) & (x<=x_max)], dist.pdf(x[(x>=dist.a) & (x<=x_max)]), \n",
    "    alpha=.8, hatch='+', ec='C5'\n",
    ")\n",
    "txt = r'$P(X<x) = P(X<a) + P(a\\leq X<x)$'\n",
    "txt += '\\n'\n",
    "txt += r'$P(X<x) = \\frac{1}{\\gamma}\\left(P_{tail}(X<a) + P_{body}(a\\leq X<x)\\right)$'\n",
    "ax1.text(\n",
    "    .875, dist.pdf(dist.a), txt, fontsize=11,\n",
    "    bbox=dict(boxstyle='square', ec='w', fc='w', alpha=0)\n",
    ")\n",
    "txt = r'$P(X<x) = \\frac{1}{\\gamma}\\left(1 + F_{\\mu,\\sigma}(x) - F_{\\mu,\\sigma}(a)\\right)$'\n",
    "ax1.text(\n",
    "    .875, dist.pdf(dist.a)/1.1, txt, fontsize=14,\n",
    "    bbox=dict(boxstyle='square', ec='black', fc='grey', alpha=.25)\n",
    ")\n",
    "ax1.grid(False)\n",
    "ax1.set_title('CDF of Left-Tailed CarBen' + r', when $x>a$', size=24)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting the pieces together results in:\n",
    "\n",
    "$$\n",
    "F_{\\textit{left}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma}\\left(1 - G_{\\xi,-a,b}(-x)\\right) & \\text{if } x < a \\\\\n",
    "      \\frac{1}{\\gamma}\\left(1+ F_{\\mu,\\sigma}(x) - F_{\\mu,\\sigma}(a)\\right) & \\text{if } x \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantile Function ##\n",
    "\n",
    "### Phat ###\n",
    "Once again, the Quantile Function $Q(.)$, known in `scipy` as the percent-point function, `ppf`, for the Phat distribution is a weighted average of the CarBens.\n",
    "\n",
    "$$\n",
    "Q(x) = \\sum\\limits_{i=1}^nw_iQ_i(x)\n",
    "\\\\Q(x) = 0.5*Q_{\\textit{left}}(x) + 0.5*Q_{\\textit{right}}(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Right CarBen ###\n",
    "Before we derive $Q(x)$ for each CarBen, it is helpful to know the quantile of the junction point between the tail and body distributions. We know this point to be the location parameter, $a$, of the generalized Pareto. We also know the relationship between the junction point and $\\gamma$. For the right CarBen,\n",
    "\n",
    "$$\n",
    "P(X>a) = \\frac{1}{\\gamma}\n",
    "$$\n",
    "\n",
    "Quantiles are, of course, merely the cumulative distribtion up to the given point. So the junction quantile, available via the `qjunc` method, is merely:\n",
    "\n",
    "$$\n",
    "P(X\\leq a) = 1 - P(X>a) = 1 - \\frac{1}{\\gamma}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, for the right-tailed CarBen, first remember the definition of the Quantile function:\n",
    "\n",
    "$$ x = Q(P(X<x)) = Q(F(x))$$\n",
    "\n",
    "where, in this framing, we refer to $P(X\\leq x)$ or $F(x)$ as the *quantile*. One important difference between the Q derivation and the PDF/CDF derivations above is that the Q function is *not* a summation of the body and tail Q functions. Q(.) identifies a *value derived specifically and solely from the distribution from which it originates*. That means, for a right-tailed CarBen, \n",
    "\n",
    "$$\n",
    "x = \\left\\{ \\begin{array}{ll}\n",
    "    Q_{\\textit{body}}(P_{\\textit{body}}(X<x)) & \\text{if } x < a \\\\\n",
    "    Q_{\\textit{tail}}(P_{\\textit{tail}}(X<x)) & \\text{if } x \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "Remember, the input to the Quantile function is a cumulative probability for the entire piecewise, hybrid CarBen. So we must conver that probability, $P(X<x)$, into $P_{\\textit{body}}(X<x)$ or $P_{\\textit{tail}}(X<x)$ to find the correct value.\n",
    "\n",
    "So, for all $x$,\n",
    "\n",
    "$$\n",
    "P(X<x) = P(X\\leq a) + P(a<X<x)\n",
    "$$\n",
    "\n",
    "Recall:\n",
    "$$\n",
    "\\text{quantile, q} = P(X<x)\n",
    "\\\\\\text{qjunc}, q_a = P(X\\leq a)\n",
    "\\\\P_{tail}(X<a) = 0\n",
    "$$\n",
    "\n",
    "so, where $x\\leq a$,\n",
    "\n",
    "$$\n",
    "P(X<x) = \\frac{P_{\\textit{body}}(X<x)}{\\gamma}\n",
    "\\\\P_{\\textit{body}}(X<x) = \\gamma*P(X<x) = \\gamma q\n",
    "$$\n",
    "\n",
    "and where $x>a$,\n",
    "\n",
    "$$\n",
    "P(X<x) = q_a + P_{\\textit{tail}}(X<x) / \\gamma\n",
    "\\\\P_{\\textit{tail}}(X<x) = \\gamma(P(X<x) - q_a) = \\gamma(q - q_a)\n",
    "$$\n",
    "\n",
    "So, for the right-tail CarBen, the Quantile function for a given $q = P(X<x)$:\n",
    "\n",
    "$$\n",
    "x = \\left\\{ \\begin{array}{ll}\n",
    "    Q_{\\textit{body}}(\\gamma q) & \\text{if } x \\leq a \\\\\n",
    "    Q_{\\textit{tail}}(\\gamma q - \\gamma q_a) & \\text{if } x > a\\\\\n",
    "\\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Left CarBen ###\n",
    "\n",
    "As we've seen, the left-tailed CarBen requires slight modifications to the formulas. For the `qjunc`, remember that that CDF of the left-tailed Pareto is equivalent to the Survival Function of the right-tailed, with location, $a$, reflected so:\n",
    "\n",
    "$$P_{\\textit{right}}(X>a) = \\frac{1}{\\gamma}$$\n",
    "\n",
    "$$\n",
    "\\text{qjunc } q_a = P_{\\textit{left}}(X<a) = P_{\\textit{right}}(X>-a) = \\frac{1}{\\gamma}\n",
    "$$\n",
    "\n",
    "again available via the `qjunc` method.\n",
    "\n",
    "For the Quantile function, we rearrange the piecewise components as follows:\n",
    "\n",
    "$$\n",
    "x = \\left\\{ \\begin{array}{ll}\n",
    "    Q_{\\textit{tail}}(P_{\\textit{tail}}(X<x)) & \\text{if } x < a \\\\\n",
    "    Q_{\\textit{body}}(P_{\\textit{body}}(X<x)) & \\text{if } x \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to express the target left-tailed `x` value and the Quantile function in terms of the right-tailed `x` and right Quantile function. This means:\n",
    "\n",
    "$$\n",
    "P_{\\textit{left}}(X<x_{\\textit{left}}) = (1 - P_{\\textit{right}}(X<x_{\\textit{right}}))\n",
    "\\\\x_{\\textit{left}} = -x_{\\textit{right}}\n",
    "$$\n",
    "\n",
    "Beginning again with, for all $x$,\n",
    "\n",
    "$$\n",
    "P(X<x) = P(X<a) + P(a\\leq X<x)\n",
    "$$\n",
    "\n",
    "so, where $x<a$,\n",
    "\n",
    "$$\n",
    "P_{\\textit{left}}(X<x_{\\textit{left}}) = \\frac{P_{\\textit{left, tail}}(X<x_{\\textit{left}})}\n",
    "{\\gamma}\n",
    "\\\\q = \\frac{(1 - P_{\\textit{right, tail}}(X<x_{\\textit{right, tail}}))}\n",
    "{\\gamma}\n",
    "\\\\P_{\\textit{right, tail}}(X<x_{\\textit{right, tail}})) = 1 - \\gamma q\n",
    "\\\\x_{\\textit{right}} = Q_{\\textit{right}}(.)\n",
    "\\\\Q_{\\textit{tail}}(.) = x_{\\textit{left}} = -x_{\\textit{right}} = -Q_{\\textit{right}}(1 - \\gamma q)\n",
    "$$\n",
    "\n",
    "And where $x\\geq a$, first, recall in the body, \n",
    "\n",
    "$$P_{\\textit{body}}(a\\leq X<x) = P_{\\textit{body}}(X<x) - P_{\\textit{body}}(X\\leq a) $$\n",
    "\n",
    "We also know that,\n",
    "\n",
    "$$\n",
    "P_{\\textit{body}}(X<a) = 2\\gamma q_a - \\gamma = \\gamma(2q_a-1)\n",
    "$$\n",
    "\n",
    "So,\n",
    "\n",
    "$$\n",
    "P_{\\textit{left}}(X<x_{\\textit{left}}) = q_a + \\frac{P_{\\textit{body}}(X<x_{\\textit{left}})}{\\gamma} - \\frac{P_{\\textit{body}}(X<a)}{\\gamma}\n",
    "\\\\\\gamma(q-q_a) = P_{\\textit{body}}(X<x_{\\textit{left}}) - P_{\\textit{body}}(X<a))\n",
    "\\\\P_{\\textit{body}}(X<x_{\\textit{right}}) = \\gamma(q-q_a) + P_{\\textit{body}}(X<a)\n",
    "\\\\P_{\\textit{body}}(X<x_{\\textit{right}}) = \\gamma(q-q_a) + \\gamma(2q_a-1)\n",
    "\\\\P_{\\textit{body}}(X<x_{\\textit{right}}) = \\gamma(q-q_a + 2q_a-1)\n",
    "\\\\P_{\\textit{body}}(X<x_{\\textit{right}}) = \\gamma(q + q_a - 1)\n",
    "$$\n",
    "In the body, realize that $x_{\\textit{left}} = x_{\\textit{right}}$\n",
    "$$\n",
    "\\\\x_{\\textit{right}} = Q_{\\textit{right}}(.)\n",
    "\\\\Q_{\\textit{body}}(.) = x_{\\textit{left}} = x_{\\textit{right}} = Q_{\\textit{body}}(\\gamma(q + q_a - 1))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "$$\n",
    "-\\gamma(1-q-q_a) = \\frac{-\\gamma(1-q-q_a)(q-q_a)}{q - q_a} =  -\\gamma + \\gamma q + \\gamma q_a = \n",
    "\\gamma(q-q_a) + P_{\\textit{body}}(X<a)\n",
    "\\\\P_{\\textit{body}}(X<a) = -\\gamma + \\gamma q + \\gamma q_a - \\gamma q + \\gamma q_a\n",
    "\\\\P_{\\textit{body}}(X<a) = 2\\gamma q_a - \\gamma = \\gamma(2q_a-1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The end results is:\n",
    "\n",
    "$$\n",
    "x = \\left\\{ \\begin{array}{ll}\n",
    "    -Q_{\\textit{right, tail}}(1 - \\gamma q) & \\text{if } x < a \\\\\n",
    "    Q_{\\textit{body}}(\\gamma q + \\gamma q_a - \\gamma) & \\text{if } x \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moments ##\n",
    "\n",
    "Here we show how to find the moments of the Phat distribution. As we know, the Phat distribution is a mixture model of two Carben distributions, each weighted 50%. Thus, the PDF of the Phat distribution is simply:\n",
    "\n",
    "$$\n",
    "f(x) = \\sum\\limits_{i=1}^nw_if_i(x)\n",
    "\\\\f(x) = 0.5*f_{\\textit{left}}(x) + 0.5*f_{\\textit{right}}(x)\n",
    "$$\n",
    "\n",
    "The Carben distribution has a piecewise density function, defined for a right-tail as:\n",
    "\n",
    "$$\n",
    "f_{\\textit{right}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma} f_{\\mu,\\sigma}(x) & \\text{if } y\\leq a \\\\\n",
    "      \\frac{1}{\\gamma} g_{\\xi,a,b}(x) & \\text{if } y > a\\\\\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "And for the left\n",
    "\n",
    "$$\n",
    "f_{\\textit{left}}(x) = \\left\\{ \\begin{array}{ll}\n",
    "      \\frac{1}{\\gamma} g_{\\xi,-a,b}(-x) & \\text{if } y < a \\\\\n",
    "       \\frac{1}{\\gamma} f_{\\mu,\\sigma}(x) & \\text{if } y \\geq a\\\\\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "The MGF for a mixture model is simply the weighted average of MGF's of its components, shown as follows:\n",
    "\n",
    "$$\n",
    "M_x(n) = \\int_{-\\infty}^{+\\infty}e^{nx}f_X(x)\\delta x\n",
    "\\\\ = \\int_{-\\infty}^{+\\infty}e^{nx}\\left(.5*f_{\\textit{left}}(x) + .5*f_{\\textit{right}}(x)\\right)\\delta x\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\\\ = \\int_{-\\infty}^{+\\infty}.5*e^{nx}f_{\\textit{left}}(x) + .5*e^{nx}f_{\\textit{right}}(x)\\delta x\n",
    "\\\\ = .5\\int_{-\\infty}^{+\\infty}e^{nx}f_{\\textit{left}}(x)\\delta x + .5\\int_{-\\infty}^{+\\infty}e^{nx}f_{\\textit{right}}(x)\\delta x\n",
    "\\\\ = .5*M_{left} + .5*M_{right}\n",
    "$$\n",
    "\n",
    "Thus, to find the moments of the Phat, we must find the MGFs of the component Carbens."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Carben Right ###\n",
    "\n",
    "The MGF of a piecewise distribution is simply the sum of the integrals along the bounded ranges.\n",
    "\n",
    "$$\n",
    "M_{right}(n) = \\frac{1}{\\gamma}\\int_{-\\infty}^{a}e^{nx}f_{\\mu,\\sigma}(x)\\delta x + \\frac{1}{\\gamma}\\int_{a}^{+\\infty}e^{nx}g_{\\xi,a,b}(x)\\delta x\n",
    "$$\n",
    "\n",
    "The body term is a truncated normal distribution. The moments of such functions are known and available in `scipy`. We know that the generalized Pareto is always restricted the interval $(a,\\infty)$, so the tail term is merely the MGF of generalized Pareto scaled by $\\frac{1}{\\gamma}$. Thus,\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\gamma}M_{\\mu,\\sigma,u}(n) + \\frac{1}{\\gamma}M_{\\xi,a,b}(n)\n",
    "$$\n",
    "where: $u$ is the upper bound of the body"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Carben Left ###\n",
    "\n",
    "Now for the left:\n",
    "\n",
    "$$\n",
    "M_{left}(n) = \\frac{1}{\\gamma}\\int_{-\\infty}^{a}e^{nx}g_{\\xi,-a,b}(-x)\\delta x\n",
    " + \\frac{1}{\\gamma}\\int_{a}^{+\\infty}e^{nx}f_{\\mu,\\sigma}(x)\\delta x\n",
    "\\\\=\\frac{1}{\\gamma}\\int_{-a}^{+\\infty}e^{nx}g_{\\xi,-a,b}(x)\\delta x\n",
    " + \\frac{1}{\\gamma}\\int_{a}^{+\\infty}e^{nx}f_{\\mu,\\sigma}(x)\\delta x\n",
    "$$\n",
    "\n",
    "Once again, the tail portion is the MGF of the generalized Pareto. And once again, the body portion is a truncated normal, this time with a lower bound, so the MGF is found as:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\gamma}M_{\\xi,-a,b}(n) + \\frac{1}{\\gamma}M_{\\mu,\\sigma,l}(n)\n",
    "$$\n",
    "where: $l$ is the lower bound of the body"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean ###\n",
    "\n",
    "The first moment is the first derivative of $M$, so:\n",
    "\n",
    "$$\n",
    "\\mu_{\\textit{right}} = \\frac{d}{dx}\\frac{1}{\\gamma}M_{\\mu,\\sigma,u}(1) + \\frac{d}{dx}\\frac{1}{\\gamma}M_{\\xi,a,b}(1)\n",
    "\\\\\\mu_{\\textit{left}} = \\frac{d}{dx}\\frac{1}{\\gamma}M_{\\xi,-a,b}(1) + \\frac{d}{dx}\\frac{1}{\\gamma}M_{\\mu,\\sigma,l}(1)\n",
    "$$\n",
    "\n",
    "In both instances, the constant multiplicative survives the differentiations. For the left tail, in order to reflect its true location, $a$, we must invert the sign:\n",
    "$$\n",
    "\\mu_{\\textit{phat}} = .5*\\mu_{\\textit{right}} + .5*\\mu_{\\textit{left}}\n",
    "\\\\\\mu_{\\textit{right}} = \\frac{1}{\\gamma}\\text{Mean}_{\\mu,\\sigma,u} + \\frac{1}{\\gamma}\\text{Mean}_{\\xi,a,b}\n",
    "\\\\\\mu_{\\textit{left}} = -\\frac{1}{\\gamma}\\text{Mean}_{\\xi,-a,b} + \\frac{1}{\\gamma}\\text{Mean}_{\\mu,\\sigma,l}\n",
    "$$\n",
    "\n",
    "Finding mean of the Phat distribution programmatically is thus fairly trivial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as scist\n",
    "\n",
    "import phat as ph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "phat = ph.Phat(.1, .2, .2, .25)\n",
    "\n",
    "bmu = scist.truncnorm(\n",
    "    -np.inf, \n",
    "    phat.right.a, \n",
    "    *phat.right.body.args\n",
    ").mean()\n",
    "tmu = phat.right.tail.mean()\n",
    "rmu = (bmu + tmu) / phat.right.gamma\n",
    "\n",
    "bmu = scist.truncnorm(\n",
    "    phat.left.a,\n",
    "    np.inf, \n",
    "    *phat.left.body.args\n",
    ").mean()\n",
    "tmu = -phat.left.tail.mean()\n",
    "lmu = (bmu + tmu) / phat.left.gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the Phat distribution is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1426574558489754"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((lmu, rmu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the mean of the Phat in the above example is close to, but slightly higher than the mean of the Gaussian body, $\\mu$. This results from the greater tail index in the right tail used in the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1426574558489754"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phat.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variance ###\n",
    "\n",
    "The second moment is the second derivative of $M$, so:\n",
    "\n",
    "$$\n",
    "\\text{Var}_{\\textit{right}} = \\frac{d^2}{dx}\\frac{1}{\\gamma}M_{\\mu,\\sigma,u}(2) + \\frac{d^2}{dx}\\frac{1}{\\gamma}M_{\\xi,a,b}(2)\n",
    "\\\\\\text{Var}_{\\textit{left}} = \\frac{d^2}{dx}\\frac{1}{\\gamma}M_{\\xi,-a,b}(2) + \\frac{d^2}{dx}\\frac{1}{\\gamma}M_{\\mu,\\sigma,l}(2)\n",
    "\\\\\\text{Var}_{\\textit{right}} = \\frac{1}{\\gamma}\\text{Var}_{\\mu,\\sigma,u} + \\frac{1}{\\gamma}\\text{Var}_{\\xi,a,b}\n",
    "\\\\\\text{Var}_{\\textit{left}} = \\frac{1}{\\gamma}\\text{Var}_{\\xi,-a,b} + \\frac{1}{\\gamma}\\text{Var}_{\\mu,\\sigma,l}\n",
    "\\\\\\text{Var}_{\\textit{phat}} = .5*\\text{Var}_{\\textit{right}} + .5*\\text{Var}_{\\textit{left}}\n",
    "$$\n",
    "\n",
    "Note that no negative is applied to the second moment as it eliminated bystanders the square.\n",
    "\n",
    "And so the variance of the Phat is found as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "bvar = scist.truncnorm(\n",
    "    -np.inf, \n",
    "    phat.right.a, \n",
    "    *phat.right.body.args\n",
    ").var()\n",
    "tvar = phat.right.tail.var()\n",
    "rvar = (bvar + tvar) / phat.right.gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "bvar = scist.truncnorm(\n",
    "    phat.left.a,\n",
    "    np.inf, \n",
    "    *phat.left.body.args\n",
    ").var()\n",
    "tvar = phat.left.tail.var()\n",
    "lvar = (bvar + tvar) / phat.left.gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4828404012432034"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lvar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variance for the Phat distribution is then: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5732923559721449"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((lvar, rvar))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is available as `var` method of the Phat distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5732923559721449"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phat.var()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}