QingFranke2014.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats, special, integrate\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Qing & Franke (2014)\n",
    "\n",
    "## General stuff about model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $\\phi (h)$ is the distribution of $h$ in the comparison class (prior)\n",
    "- $\\Theta$ is the threshold for \"tall\"\n",
    "- $\\phi (h \\mid u_0, \\Theta) = \\phi(h)$\n",
    "- $\\phi (h \\mid u_1, \\Theta) = \\phi(h \\mid h \\geq \\Theta) $, i.e. the survival function of $\\phi$ at $\\Theta$, if $h\\geq$ and 0 otherwise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected communicative success:\n",
    "\\begin{align}\n",
    "ES(\\Theta) & = \\int_{-\\infty}^{\\Theta} \\phi (h) \\phi (h \\mid u_0, \\Theta) dh + \\int_{\\Theta}^{\\infty} \\phi (h) \\phi (h \\mid u_1, \\Theta) dh \\\\\n",
    "& = \\int_{-\\infty}^{\\Theta} \\phi (h)^2 dh + \\int_{\\Theta}^{\\infty} \\phi (h) \\frac{\\phi(h)}{ \\int_{\\Theta}^{\\infty} \\phi (h) dh} dh \\\\\n",
    "& = \\int_{-\\infty}^{\\Theta} \\phi (h)^2 dh + \\frac{1}{ \\int_{\\Theta}^{\\infty} \\phi (h) dh} \\int_{\\Theta}^{\\infty} \\phi(h)^2 dh\n",
    "\\end{align}\n",
    "\n",
    "Calculate utility:\n",
    "\n",
    "\\begin{align}\n",
    "U(\\Theta) & = ES(\\Theta) - c \\int_{\\Theta}^{\\infty} \\phi(h) dh \\\\\n",
    "& = \\int_{-\\infty}^{\\Theta} \\phi (h)^2 dh + \\frac{1}{ \\int_{\\Theta}^{\\infty} \\phi (h) dh} \\int_{\\Theta}^{\\infty} \\phi(h)^2 dh  - c \\int_{\\Theta}^{\\infty} \\phi(h) dh\n",
    "\\end{align}\n",
    "\n",
    "Now the probability of production:\n",
    "$$\n",
    "Pr(\\Theta) \\propto e^{\\lambda U(\\Theta)} \n",
    "$$\n",
    "\n",
    "Probability of using positive form given an observed height $h_0$:\n",
    "$$\n",
    "\\sigma (u_1 \\mid h_0, Pr) = \\int_{-\\infty}^{h_0} Pr(\\Theta) d\\Theta\n",
    "$$\n",
    "\n",
    "Understanding, i.e. distribution over degrees given signal:\n",
    "$$\n",
    "\\rho (h \\mid u_1) \\propto \\phi(h) \\sigma ( u_1 \\mid h, Pr )\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifically for Beta distribution. \n",
    "\n",
    "Start by considering the Chebyshev integral: http://mathworld.wolfram.com/ChebyshevIntegral.html:\n",
    "$$\n",
    "\\int x^{p} (1-x)^q dx = B(x; 1+p, 1+q)\n",
    "$$\n",
    "\n",
    "I am interested in the integral of the square of the beta density:\n",
    "$$\n",
    "\\int_0^\\Theta \\left( \\frac{1}{B(\\alpha, \\beta)} x^{\\alpha-1} (1-x)^{\\beta-1} \\right)^2 dx \\\\\n",
    "= \\frac{1}{B(\\alpha, \\beta)^2} \\int_0^\\Theta x^{2\\alpha-2} (1-x)^{2\\beta-2} dx \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting the two things together:\n",
    "$$\n",
    "\\frac{1}{B(\\alpha, \\beta)^2} \\int_0^\\Theta x^{2\\alpha-2} (1-x)^{2\\beta-2} dx = \\frac{B(x; 2\\alpha-1, 2\\beta-1)}{B(\\alpha, \\beta)^2}\n",
    "$$\n",
    "Check numerically that they are actually the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.713713169744543\n",
      "4.713713169744545\n"
     ]
    }
   ],
   "source": [
    "theta = 0.1\n",
    "alpha, beta = 2,21\n",
    "\n",
    "func = lambda x: (x**(2 * alpha - 2)) * ((1-x)**(2*beta - 2))\n",
    "\n",
    "print(quad(func, 0, theta)[0] / (special.beta(alpha, beta)**2))\n",
    "\n",
    "print((special.betainc(2*alpha-1, 2*beta-1, theta) * special.beta(2*alpha-1, 2*beta-1))/ (special.beta(alpha, beta)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, I need it from $\\Theta$ to 1:\n",
    "$$\n",
    "\\frac{1}{B(\\alpha, \\beta)^2} \\int_\\Theta^1 x^{2\\alpha-2} (1-x)^{2\\beta-2} dx \\\\\n",
    "= \\frac{1}{B(\\alpha, \\beta)^2} \\left( \\int_0^1 x^{2\\alpha-2} (1-x)^{2\\beta-2} dx - \\int_0^\\Theta x^{2\\alpha-2} (1-x)^{2\\beta-2} dx \\right) \\\\\n",
    "= \\frac{B(2 \\alpha-1, 2\\beta-1) - B(x; 2\\alpha-1, 2\\beta-1)}{B(\\alpha, \\beta)^2} \n",
    "$$\n",
    "\n",
    "In sum, when the property has a beta distribution:\n",
    "$$\n",
    "U(\\Theta) = \\frac{B(\\Theta; 2\\alpha-1, 2\\beta-1)}{B(\\alpha, \\beta)^2} + \\frac{B(2 \\alpha-1, 2\\beta-1) - B(\\Theta; 2\\alpha-1, 2\\beta-1)}{(1-I_\\Theta (\\alpha, \\beta)) B(\\alpha, \\beta)^2} - c (1-I_\\Theta (\\alpha, \\beta))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [],
   "source": [
    "def U(theta, alpha, beta, cost):\n",
    "    incomplete_beta_for_cdf = special.betainc(2*alpha-1, 2*beta-1, theta) * special.beta(2*alpha-1, 2*beta-1)\n",
    "    beta_for_cdf = special.beta(2*alpha-1, 2*beta-1)\n",
    "    beta_squared = special.beta(alpha, beta)**2\n",
    "    survival = 1 - special.betainc(alpha, beta, theta)\n",
    "    \n",
    "    first_term = incomplete_beta_for_cdf / beta_squared\n",
    "    second_term = (beta_for_cdf - incomplete_beta_for_cdf)/(survival * beta_squared)\n",
    "    third_term = cost * survival\n",
    "    return first_term + second_term - third_term\n",
    "\n",
    "def unnorm_pr(theta, alpha, beta, cost, lamb):\n",
    "    return np.exp(lamb * U(theta, alpha, beta, cost))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [],
   "source": [
    "def move_legend(leg, ax, x_offset):\n",
    "    # Get the bounding box of the original legend\n",
    "    bb = leg.get_bbox_to_anchor().inverse_transformed(ax.transAxes)\n",
    "\n",
    "    # Change to location of the legend. \n",
    "    bb.x0 += x_offset\n",
    "    bb.x1 += x_offset\n",
    "    leg.set_bbox_to_anchor(bb, transform = ax.transAxes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:24: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:24: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \"\"\"\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:3: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAADnCAYAAABSbO4uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydeVxU5ffH35d9FQQFFETcRVFxS80yza3MzDWz1bRc0spss/Vn9W3PNUtt0azcNS0zyyWzNJdEUdx3xQVF2WSHmef3x8MAw8zAADPMDN7363VfA3d57hm9fObMeZ5zjiKEQEVFRUXF+jjZ2gAVFRWVWwVVcFVUVFSqCFVwVVRUVKoIVXBVVFRUqghVcFVUVFSqCFVwVVRUVKoIFyuNq641U6ksig3uafq5PXcOGjSAhQth5MgqM0jFdnTuDDk5sH9/uS4r9blVPVwVFXNwKvhT0Wpta4dKlXDuHOzeDcOHW3ZcVXBVVMxBJ7gajW3tUKkSVqyQrw8+aNlxVcFVUTEHZ2f5qnq4twTLl0PHjtCwoWXHVQVXRcUcVA/3luHUKdi3z/LhBFAFV0XFPFQP95Zh+XL5aulwAqiCq6JiHqqHe8uwfDncfjvUq2f5sVXBVVExB9XDvSU4ehTi4qwTTgBVcFVUzEP1cG8JVqwARYGhQ60zviq4KirmoHq41R4hZDihWzeoW9c691AFV0XFHFQPt9pz6JAMKVgrnACq4NqE+fPnExISQnR0NA0bNuS7774zel5WVhZ33XUXmoI/8vj4eAYMGECTJk1o1KgRzz//PLm5ueW690MPPVRZ840SHx9Pjx49iIyMpGXLlsyaNavwWG5uLt26dSM/P98q964SLOjhvv3227Rq1YqmTZvy1VdfVWiMUaNGERQURFRUVKXtUZEsWyY/V4cMsd49VMG1AQcPHmTq1KnExsayatUqXnzxRYNzNBoNCxYsYPDgwTg7OyOEYNCgQQwaNIiTJ09y4sQJ0tPTeeONN8p173PnzlnoXejj4uLCtGnTOHr0KLt27eKLL77gyJEjALi5udGzZ0+W69bbOCIW8nD/+OMP9u/fT2xsLKtXr2bt2rUVGmfkyJH8/vvvlbJFpQitFhYvhj59ICjIevdRBdcGxMXFERkZCUBYWFihBzts2DAmT55Mjx49+PDDD1m8eDEPPPAAAFu2bMHLy4snn3wSAGdnZ2bMmMGCBQvIzMystE2XL19myJAhtG3blubNm7Nnz55yXV+nTh3atWsHgK+vL5GRkVy6dKnw+MCBA1m8eHGl7bQZFqql8MsvvzBy5Ejy8vKYM2cOQyroTnXr1o2AgACzzz98+DC9evWiadOmvPfeezz77LP8999/Fbp3deTff+H8eXj0Uevex1rVwuyeSZMgNtayY0ZHw8yZZZ8XFxdH8+bNEUIwe/Zs+vfvX7g/MjKSrVu3kpuby5w5c4iIiADgyJEjtG/fnn79+nH58mUABgwYQHh4OKdOnaJ169YVtjs/P597772X999/n/79+5OZmVn4IQBw5513cvPmTYPrPvvsM3r16mWw/9y5c+zfv59OnToV7ouKinLsP3BFkVslPdyYmBg6duxIYGAgERERzJgxQ+94ef+tzSE7O5thw4axcuVKGjZsSPPmzWnfvj0dO3a06n0diR9/BC8vKPBvrMYtK7i2Ij4+nvT0dPr27Yurqyu33XYbX3zxBdnZ2SQlJfH2228DcP36dfz9/QuvE0KgKAq//fab3ni//PILilJUES4uLg4/Pz/Cw8PNtmnt2rVERkYWCr+Xl5fe8X/++cfssdLT0xkyZAgzZ86kRo0ahfudnZ1xc3Pj5s2b+Pr6mj2eXeHsXCkPV6vVcvHiRUaOHMnw4cMZO3Ys06dP58033yw8pzz/1uayefNm2rZtS8uWLQEZUy8ZxrLGfR2FnBy5HGzQIPDxse69blnBNccTtQYHDx6kZ8+eBvG3mJgYOnXqhIuL/C/x9PQkOzu78HhUVBSrV6/WuyYtLY34+HgaNWqkN05ERES5BDc2NpbOnTubPG6u95OXl8eQIUN45JFHGDx4sMH5OTk5eHh4mG2X3eHkVCkP9/jx4zRp0gSQ/79du3YlISFB7xxreJr79+8vDPdcvnwZHx8funbtavX7OgobNkBysvXDCXALC66tiIuLo02bNkb3Fw8L1KxZE41GQ3Z2Nh4eHtx9991MmTKF77//nscffxyNRsMLL7zAqFGj8PLy4siRI8yePZvjx4/j6+vLsmXLeOeddwgODi7TppCQEA4cOFD4e2JiIrVr1y783RzvRwjB6NGjiYyMZPLkyQbHb9y4Qe3atXF1dS1zLLulkh7u/v37ycnJQaPRkJ+fz5IlS5g9e7beOZbwNHv27Mn3339PaGgoAO7u7ly8eBGA1157zejKllvZw128WE6UVcXnijppVsWUFNbS9vfp04ft27cDoCgKa9asYdWqVTRp0oQmTZrg7e3N+++/D0CLFi2YN28eTzzxBJMnT2bevHlmiS3IGe+rV6/SsmVLoqOj2blzZ7nf144dO/jhhx/4888/iY6OJjo6Wi/8sXXrVvr161fuce2KSnq4sbGxZGVl0ahRI7p27coTTzxh9MPXHEaMGEGXLl04fvw4YWFhfPvtt4AMW5w6dUpvQu3hhx/m77//plmzZrRp04YuXbowadKkCr+P6kRKCqxbBw89BC5V4X4KIayxqViAffv2iUcffbRc1yxcuFBs3brV5PFOnTpV0qqKMWjQIHHs2LHyXGKtZ7Piz62vrxCTJpXnPejRq1cvERcXV+HrzSEuLk688MILVr1HdeKbb4QAIfbssdiQpT5jDh1SyNfmo9FqcHdxt7UpVqFt27b06NEDjUaDs27hfRmMtMN+W7m5uQwcOJBmzZrZ2pTKUcmQwrFjx2jevLkFDTIkKiqK6dOnW/Ue1YnFi6FpU+jQoWrupwhhlX6PVm8iqRVa7vruLk7eOMn+sfup41vH2rdUqVrsq4kkQGAgjBgBc+ZUkTkq1iQ+HurXh6lToWBxkCWonk0kd8bvZPuF7VzNuMqSuCW2NkflVqCSHq6KfbF0qSxY88gjVXdPhxXcLWe3oKAQ4hPCxjMbbW2Oyq1AJSfNVOwHIeD776FLFyi2qtLqOKzgxlyJoXmt5gxsNpDdF3djpdCIikoRqodbbdi7Fw4fhqqe8nBYwT2SeISWQS1pGdSS1JxUEtITyr5IRaUyqB5utWHhQvD0tG4pRmM4pOBm52dzJvkMLWq1ILKWLAJzJPGIja1SqfaoHm61IDtbxm8HDwY/v6q9t0MK7tnks2iFlma1mtG8llxmc+LGCRtbpVLtcXICR67pqwLA2rUy4aGg8F6V4pCCeyH1AgDhfuGE+ITg4uRCfFq8ja1Sqfa4ujpsSOH333+nWbNmNG7cmI8++sjg+LFjx+jSpQvu7u589tlnNrCw6li4UC4H69Gj6u/tkIKrE9d6Nerh7ORMXd+6XEy7aGOrVKo9Li4O6eFqNBomTJjAhg0bOHLkCEuXLi0sDq8jICCA2bNn89JLL9nIyqohPh42bYInnigqcVyVOKbgpsajoFDXV3Z6q1ejnkMJrtpix0FxUMHds2cPjRs3pmHDhri5ufHQQw/x888/650TFBREx44dHbu4kBn88INcEmarhEyHTO2NT4unjm8dXJ3lwxFWI4yYKzE2tsp8dC12xo0bx759++jdu7dBSq6pFjsTJkzgl19+QaPRMGbMGN544w0+/fRTs+9t7RY77dq14+bNm7Rv357evXvTokULvRY7j1TlKnNLU1nBtVHV+0uXLlGvXr3C38PCwti9e7dl7XAAhJDhhO7doUED29jgmB5uWjz1ahR7gGqEcTHtosOsxVVb7DgoDurhGvu7KF60/lZhxw44dco2k2U6HNLDvXzzMs0Ciwqh1KtRj+z8bJKykgj0CjRrjEm/TyI2wbLeRnRINDPvKbuyudpix0GprODaqOp9WFgY8fFFk8oXL16kbt26NrHFlixcKDs6WLMrb1k4pOBey7hGt/Buhb+H1QgD4GLaRbMF11ZYu8VORVBb7JiJg3q4HTt25OTJk5w9e5bQ0FCWLVvGkiW3Vv2RmzdlG50HHwRvb9vZ4XCCq9FquJF5gyDvol7Gup8TMxPNHsccT9QaWLvFTkVQW+yYiYuLbIDlYLi4uDBnzhz69u2LRqNh1KhRtGzZknnz5gEwbtw4EhIS6NChA2lpaTg5OTFz5kyOHDmi96HpyCxdCunpMGaMbe1wOMG9kXUDgdAT3FpetQBIzDBfcG2FtVrsxMXF8dprr+mNuWDBAoKCgkreygC1xY6ZuLpCRoatragQ/fr1M+i4MW7cuMKfQ0JCCtvwVEe++gpat4bbbrOtHQ43aXYt4xqAnuDW9pbiUB4P11ZYq8VOq1at+PXXX/U2c8QW1BY7ZuOgIYVbnZgYuY0dKzvd2xKH83CNCW6AZwBOipNDeLimZuqnTZtmsG/ixIlMnz698Gt7WFgYv/zyi8Vt8vHxqfS4d9xxR6mrRJYsWcKHH35YqXvYHFVwHZKvvpKFauxhRWK18HCdFCcCPQO5nnndVmZZheItdhyZatNiRxVch+PmTViyRDaJrOpCNcaoFh4uyDiuI4QUysuoUaMsOt6uXbssOp45uLm58fjjj1f5fS2OKrgOh26ybOxYW1sicUgP11lxpqZnTb39tb1rV0vBVbEjVMF1OOxlskyHwwluYkYigV6BOCn6ptf2ql3tQgoqdoYquA6FbrJszBjbT5bpcDjBTcpOIsAzwGB/La9aDjFppuLAuLhAXp6tragQZZVn/Pnnn2ndujXR0dF06NChcHWMI6ObLHv0UVtbUoTDxXCTsowLbm2v2tzIuoFWaA28XxUVi+Dq6pAerq4846ZNmwgLC6Njx44MGDCAFi1aFJ7Ts2dPBgwYgKIoHDx4kAcffJBjx47Z0OrKYW+TZTocTpmSs5KNC653bbRCS1JWkg2sUrklcNCQgjnlGX18fApTxDMyMhy+uM3339vXZJkOh/Rwo4KiDPbrss2uZ14v/FlFxaJUUnAnTZpErIXLM0ZHRzPTQuUZ16xZw2uvvca1a9dYv369Re2sSoSAOXOgY0coVj/JLnA4D7e0kAI4RnqvioPioB6uueUZBw0axLFjx1i7di1vvfVWVZhmFf78E44dg4kTbW2JIQ7l4eZp8riZe9Oo4OqWiaVkp1S1WSq3CpUU3LI8UWtR3vKM3bp14/Tp01y/fp1atRzv2+Lnn0OtWrIymL3hUB6uTkxretQ0OKbb5wiCWx1b7IDpmXC1xY5tKV6eMTc3l2XLljFgwAC9c06dOlXoCe/bt4/c3FwCA+271Kkxzp2DdevkUjB7LEznUIKrmxAz5uH6e8jascnZyVVqU0XQtdiJjY1l1apVvPjiiwbnmGqxM2jQIE6ePMmJEydIT0/njTfeKNe9rdVip7RGhcVb7Dg0Li6g1crNgShenjEyMpIHH3ywsDyjrkTj6tWriYqKIjo6mgkTJrB8+XKHnDgreDsUK4RmXwghrLFZhZ3xOwVTEetPrDc4lqfJE0xFvLN1qhDXdgiRddVaZlSaO++8U/z1119CCCGuXr0q/Pz8hBBCDB06VLzwwguie/fu4r333hNdunQRZ8+eFUIIsWnTJnHnnXfqjZOamioCAgJERkaG2ffu1KmT0f2XLl0SgwcPFtHR0aJZs2Zi9+7d5XpP//77r+jTp0/h7x988IH44IMPCn+PjY0V9957b3mGtNazWfHn9n//EwKEyMkpz/tQqSIyM4UIDBRi8GCbmlHqM+ZQMdzSPFwXJxd83HxokvovbJoKnqHQ/xi4+hgfzEYN/aB6ttgpaya82rTYARlWcHOzrS0qBixbBjdu2OdkmY5qI7gANT386Z65CxRnyLoEZxZCs2er0sQysXaLnbi4OPz8/AgPDzfbJku02BFlzIRXixY7uuLpDhjHre4IISfLWraUXXntFYcS3OQsGZ81NmkGEO3lSR1xETp+CSe+hEu/mBZcG80YW7vFTkxMDBEREeUSXEu02DFnJrxatNgBVXDtkJ07Yf9+mDvXfuomGMOhBFfn4ZasFKaju0fBZEadeyDtJJz8EvKzwMWzqkwsE2u12Dly5AizZ8/m+PHj+Pr6smzZMt555x2Cg4PLtMkSLXbKalRYLVrsqIJrt8yYAf7+9lU3wRgOt0qhhnsNXJyMf060ds0lSesC3hEQcjdocyDJvuKG1mqx06JFC+bNm8cTTzzB5MmTmTdvnlliC5ZpsWNqJlxHtWmxA6rg2hlnz8JPP8k0Xh8TUzb2gkN5uMnZySbDCQDNlJscyHOhh6JAzbYFFx2EoG4mr6lqqmuLHTDeqFBHtWmxA6rg2hmzZoGTEzxrX9M1RnEoDzclO6Vwva0B2nzqiDQO5BSEFTzrgnsgpBwwfr4DUJEWOyNHjqS7nc0aVKsWO+CQgltWeUYd//33H87OzqxataoKras4KSnw7beyKlhoqK2tKRuHEtzUnFT8PEzUWss4hwta4rJy0QqtjJz7t4ZkxxVckC12nJ2dLTae2mKnEugE18Fq4paWlFLyvFdffZW+ffvawMqK8fXXsirY5Mm2tsQ8HEpwU7JT8HM3IbhpxwE4lgtpOWlyn38bSD0EwrEyg1TsFN2En4MJrjnlGQE+//xzhgwZQlBQkJFR7I+8PJg9G3r0gLZtbW2NeThUDDc1OxX/YBMhhQLBPZ4rl4/5e/hDjaagyYKsy+AVVoWWqlRL3N3la05OhS6f9PskYhMsXJ4xJJqZ91S+POOlS5dYs2YNf/75p8MkqKxYARcvFqXzOgIO5eGm5qSa9nDTz5Dn7MUNbbECNj4F61Nvnq4aA1WqN7rssnIWDLI1ZSWlgKzV+/HHH1s0fGVNhIDp06F5c7j3XltbYz4O4+FqhVZ6uKYmzTLOk+MeApwxFNz00xB8V5XYqVKNqaSHW5Ynai3MSUrZu3dvYSW569ev89tvv+Hi4sLAgQOr1FZz2bYN9u2D+fPlCgVHwWEENz03HYEwPWmWeR6NVxh6gusdLtN801UPV8UCVFJwbUVZSSkAZ8+eLfx55MiR9O/f327FFuDTT2XN28ces7Ul5cNhPhtSs1MBTIcUMi6geEcAxUo0OrnKJAg1pKBiCRxUcM0pz+hIxMbCb7/J+lOe9pNEahYO4+HqvFajIYXcVMhLxdW3sd65gAwrqB6uiiVwUMEF40kp40wUjTVVEN9e+Ogj8PWFCRNsbUn5cRwPN6fAwzUWUsi8AIC7XxMUFH3B9VUFV8VCOOikWXXi1ClYuRLGj5e1ExwNhxFcnYgaDSlkSMF18qqPn4dfCQ+3IeQmSy/YTlBb7DgoDuzhVhc++UQuh37hBVtbUjEcRnB1MVyjIYUsWZAbr1BqetTUb7PjGap/jh2gtthxUFTBtSmXLsGiRTBqFISE2NqaiuE4gltaSCHrinz1CMHfw1/fw/UsWP5iR4Kr6+wAcsmOzoMdNmwYkydPpkePHnz44YcsXryYBx54AIAtW7bg5eXFk08+CciC3jNmzGDBggVkZmZW2qbLly8zZMgQ2rZtS/PmzdmzZ0+5ri8rm2ngwIEmC/c4DKrg2pTp00GjgZdftrUlFad6TJplJ8hCNc5uhoLrVeDhZl7SvyZmEiRbuMVOzWhor7bYgWraYkcVXJtx44Zcc/vQQ9Cgga2tqTgOI7ip2am4Obvh4WKkY0DWFfCoA0hBPpl0suiYZ52Cc+zDw7V2i52KoLbYMRPdpJkquFXOnDmQkQFTptjaksrhOIJbWlpv1hXwlEGdmh41C1vxAODiDa5+hoJrhidqDazdYqciqC12zMTFRaY1OeAqhd9//53nn38ejUbDU089xZQSyvXpp58Whnzy8/M5evQoiYmJBAQY7x9YlaSmypq3AwZAVJStrakcDiO4pdbCzboCvk0BDEMKIMMKWZeMXFj1WKvFTlxcHK+99premAsWLDCr8pPaYqccuLs7nIerm9DctGkTYWFhdOzYkQEDBtCiRYvCc15++WVeLgiOrlu3jhkzZtiF2IKsCJacDP/3f7a2pPI41KSZ0QkzIWQM17MopJCRl0GeplgJPc+6kGkfIQVrtdhp1aoVv/76q95mbpk9tcVOOXBAwTW3PKOOpUuXMmLEiCq00DSpqXKybMAAaNfO1tZUHrM8XEVRVgMLgA1C2Ka4bGq2iZBCbjJocwsFV9dgMjUnlVpeteQ5nnUh7a8qsrR01BY7Dk4lBHfSJJmWakmio8tuQG1OeUYdmZmZ/P7778yZM8eSZlaYWbNkV4epU21tiWUw18OdCzwMnFQU5SNFUZpb0SajmAwp6JaEFfNwAf04rmeojOE6WCHyirTYsUeqTYsdkBNnDubhmlOeUce6devo2rWrXYQTUlKkd/vAA45TYLwszPJwhRCbgc2KovgBI4BNiqLEA18DPwohrF4C3+SkWXaCfPWQk2Y6wTVYiyvyIec6eDhGNXsdo0aNsuh4aoudSuLuXuFJs7I8UWthzoSmjmXLltlNOGHWLBlSqA6xWx1mx3AVRQkERgJPAfuBWUA7YJNVLCuBaQ+3QHBLeLj6a3HtL/lBxUFxwBhu8QnN3Nxcli1bxoABAwzOS01NZdu2bYXJNrYkJQVmzICBA6uPdwvmx3B/ApoDPwD3CyEKvsezXFGUvdYyTkeeJo/MvEzjk2Y51+Rrgeeqa6Ou7+EWS36oGW1NU1WqO56ekJVlayvKRfEJTY1Gw6hRowrLM0JR1bA1a9bQp08fvL29bWkuIL8N2Jt3e/DqQb7d9y1X0q+wYtiKCo1h7rKwb4QQeivuFUVxF0LkCCE6VOjO5UDXFNJ4SOGarHvrKo8VxnD16imoHq6KhfDyAgukUlc15pRnHDlyJCNHjqxCq4yTlCS920GD5KSgLUnNTmXZoWV8s/8b9l7ei5uzG4MjB5OvzcfFqfyras294n/AbyX27USGFKxOqWm9OYngXku2RcdUDDcEUAzTe1VUyou3NyQm2tqKas1HH8HNm/Duu7a5vxCC7Re2883+b1h5eCVZ+VlEBUUxo+8MHm39aNHqpwpQquAqihIChAKeiqK0BXRTmzUAL5MXWphSC9dkJ4J70SJ9L1cvXJxc9AXXyVWGHFQPV6WyOKiH6yhcvAiffy5b51R1VllCegLfH/ieb/d/y4kbJ/B18+Wx1o8xut1oOtbtWOkUeijbw+2LnCgLA6YX238TeL3SdzeTUtvr5OgLrqIo1PSoaZht5llXFVyVyuPtLZP6VazCu+/KimDvvFM198vX5vP7qd/5Zt83/HriVzRCwx3hd/DaHa8xrMUwvN0sG88uVXCFEIuARYqiDBFCrC7tXGtSeqWwRAjUDyP7e/jrx3ChQHBtH1JwdnamVatW5OfnExkZyaJFiwyKxYCMuS1ZskSvgI2KHaAKrtU4fhwWLJCtcwqK5FmNU0mnWLB/Ad/FfseV9CsEeQcxuctkRrUdRfNa1kszKCuk8KgQ4kcgQlGUySWPCyGmG7nM4pQaUijh4UKB4GaVEFyvUEiyfXlAT09PYgvSjR555BHmzZvH5MlF/7RCCIQQBlXBSkN3jZMj9Yt2VNSQgtV48025CKScNfXNJleTy9pja/kq5iu2nN2Ck+JEvyb9GN12NPc1uQ9XZ+vX+SjrL1TnT/sAvka2KsFkSEGTC3mpRgVXJ9KFeNaVKxq0Vs/RMJs777yTU6dOce7cOSIjI3nmmWdo164d8fHxREREcP36dQCmT59OVFQUUVFRzCxYPW/sGpUqwNtbCq7WsbIW7Z3//oNVq+DFF8HMEiBmczrpNFM2T6HejHoMXzWcU0mneK/He1yYdIF1I9YxsPnAKhFbKDukML/gtYoiKsYp7GdW0sPNkYKEh77g1vSsSXxaCQEqXBp2BbzDmTRpUqGnaSmio6MLBbEs8vPz2bBhA/fccw8Ax48fZ+HChXz55Zd658XExLBw4UJ2796NEIJOnTpx1113UbNmTZPXqFgR3RrV7Gzp7apYhNdfh1q1YLLB9+iKkafJ4+fjPzM/Zj6bz2zGWXHm/mb3M7b9WHo37I2zk7NlblROzPoOqijKJ4qi1FAUxVVRlC2KolxXFOVRaxunIzUnFW9Xb8N1bzkFy3NKerjuRkIKdtLbLCsri+joaDp06EB4eDijR48GoH79+kZr0m7fvp1Bgwbh7e2Nj48PgwcPLiyXaOoaFSuiE1k1jmsxNm6EzZul6NaoUbmxziSf4fUtr1NvRj2GrRzGiRsneLf7u5yfdJ41w9dwT+N7bCa2YP463D5CiFcURRkEXASGAVuBH61mWTFSslNMx2/BwMM1XhNXP/nBXE/U0hSP4RbHVHaPscIjZV2jYkV0/+YZGVC7dunnqpRJfr4MIzRsCM88U7Ex8jR5rDuxjvkx89l4eiNOihP9m/ZnbPux9G3U16YCWxJzBVcX4OgHLBVCJFliTZq5pOakFqbs6pFtwsP18CdHk0N2fnZRSx5PE73N7Jxu3boxcuRIpkyZghCCNWvW8MMPP9jarFsXneCqE2cWYcECOHRIxm91LePM5VzKOb7Z9w3f7v+WhPQE6tWoxzvd32FU21GE1QizjsGVxFzBXacoyjEgC3hGUZTaQHYZ11iM5Kxk01lmYFRwQXrGIT4F/ZTdA2UChIOtxW3Xrh0jR47ktttuA+Cpp56ibdu2Vmt3rlIGakjBYqSlwVtvwZ13wuDB5l2jFVo2nt7IF/99wfoT61EUhfua3MeY9mO4t/G9duXNGsPc8oxTFEX5GEgTQmgURckAqqykUEp2CnV9jZSTy04ExQnc9Wt36oqQJ2clFwmu4iQbTdpYcNPT0w32RUREcOjQIb19xQV18uTJekvHTF2jUgWoHq7F+OADuHYN1q8vzMw3SVJWEgv3L2Tu3rmcTj5NsHcwb3Z7k6fbPU09v3qlX2xHlKf6QiRyPW7xa763sD1GSclOoUXtFoYHchLBLVCKaTGM1lMAuRbXwUIKKnZG8RiuSoU5e1YWqHn8cehQSvmrfVf28cWeL1h6aClZ+VncEX4H/7v7fwyOHIybs1vVGWwhzC3P+APQCIgFdO0HBFUouCZDCh6GExcmBSY8QqoAACAASURBVNezLqQesYaJKrcKakjBIkyZAs7OUNCST4+c/BxWHlnJl/99yc6LO/Fy9eKx1o/xTMdnaBNi2IDVkTDXw+0AtBClTZlbCa3Qlj5p5l5OwU2oknrpKtUV3bolI23jVcxjxw5YsULWug0rNrd1IfUC8/bO45t935CYmUjTwKbM7DuTJ6KfMN2x28EwV3APASHAlbJOtDTpuelohda0h+tnWFLIaE1ckCGFvDTISwdXn9JvnJgoK2ikpsqP42IdaFVuYXS1LVJSSj9PxSgaDTz7LISGwssvy2WPW89tZfbu2aw7sQ6A+5vez4SOE+jZsCdOSvVKVzdXcGsBRxRF2QMU9hcRQhj26bAwZdbCLW9IAWS2mWsT0zfNyoI+feDwYfDwgF9/hf37rV9RQ8X+8fEBJydVcCvI/PnyT+mHJbksP/EjM3fNJO5aHLW8avFq11cZ234s9f3r29pMq2Gu4E61phGlYVJwtRrISTIaUvBw8cDDxaMUwb0ENUoR3JkzZT/rdesgMhLatZMljNavr8xbUXFAkpOTeeeddxg+fDhdunSRYuvnpwpuBUhMhNff0BLR9hwvxHfm+olEWge3ZsGABYxoNaJozXw1xix/XQixDTgHuBb8/B+wz4p2FWJScHNvAMKo4OrONxRcM9J7MzLgs8+gXz/o3x8aNYLXXoPffoM9eyr4LlQclaysLGbNmsXBgweLdvr7Q3Ky6YtUDIi5HEOXEdtITdNw7vb+3B7ehT8f/5PYsbE82fbJW0JswfxVCk8DY4AA5GqFUGAe0NN6pklMCm5hlpnxdhdGa+Ka07137VrZVOmVV4r2TZggFw3OnQsFCQgqtwbuBelPOcU79fr7qx6uGeRr81l7bC2zds9i+4582LKT9kO3sPzddTQKaGRr88xGo9GQkJDAhQsXuHDhArm5uTz22GMVGsvckMIE4DZgN4AQ4qSiKBYuomYcXREaA8E1USlMh9GuDy6+4OJd+lrc77+H+vVl+osOX18YPhyWLpX9P3zKmHBTqTa4ucm1nrm5uUU7VcEtldTsVL7e9zWf7/mcC6kXaODXmLB/dqAN1bJ1QU98q6ywq3mkpaURHx9fKKglt4sXL5Kfn194fp06dawuuDlCiFxd/YSC5IcqWSJm0sM1kdarw9/Dn8TMEs3+FEWGFUx5uDduyLJFU6bIWF1xRo6Eb76RHvCjVVYoTcXGmPRwT560kUX2y8W0i8zaNYv5MfO5mXuT7hHdmX3PbOK33M+zJ5xYtowqF1utVktCQgLnzp3j3LlzRgU1NVW/draLiwthYWGEh4dzxx13EB4errfVq1fxzDZzBXeboiivI5tJ9gaeAdZV+K7lwHQtXOOVwnT4e/hzMsnIH0Vpvc02bZKFpe+/3/BYly4QEiIn0lTBvWVwdZV1m1QP1zRxV+P4bOdnLIlbghCCYS2H8fLtL9OuTjsuXYLHX4eePeHBBy1/75KCWnI7f/68/v8dEBAQQHh4OA0aNOCuu+4yENSQkBCcna1Tk8FcwZ0CjAbigLHIlunfWMWiEqRkp+Dr5mtYC1cXw3ULNHqd0TY7IAX3+k7jN/v9dwgIgI4dDY85OcF998HKlZCbC26Ol1aoUn4URcHNzU2N4ZZACMFf5/7i038/ZcOpDXi5evFMh2d4ocsLRPhHFJ737LPyz2XevLLrJRhDq9Vy5coVk4Kqi6kWJygoiIiICNq2bcugQYOIiIggIiKC+vXrEx4ejo8NQ4LmFq/RKoqyFlgrhEgs8wILkpJjKq33Orj6gYl8at0qBSGEfntjr4KQghD6T4AQ8Mcfcv2tqU+3/v3h22/hn3/kR7bKLYG7u7uh4Kan35IfvPnafFYfWc2n/35KzJUYgryDeK/He4zvMJ5AL33nZ80auX34ITRubHrMpKQkzpw5w+nTpzlz5gxnzpwpVVCDg4OJiIigffv2DB48uFBQdaJqrCmrvVBWE0kF+D9gIqAU7NIAnwsh3q0C+0qvo2Aifgty0kwjNGTkZeDjVuwTzbMuaHMgN0mWbNRx7BgkJECvXqaN6dULXFxknFcV3FsGNzc3/T/64GD5mpgoU6ZuAbLzs1mwfwGf/fsZZ1PO0jSwKfP7z+fxNo8bXdKVmioX97RpA5MmaTh3Ll5PVE+fPl34c0qJbwtBQUE0aNCA9u3bM2TIED1BDQ8Pt2tBLYuyPNxJQFegoxDiLICiKA2BuYqivCCEmGFtA00KbnaiySVhoJ9tZiC4IL3c4oL777/y9Y47TBvj4yOXhf31l5nWq1QHDDxcneBevVrtBTc9N535e+fz2c7PSEhPoFNoJ6b1mcaAZgP0as+mp6dz9uzZQiFduPA2rly5HReXB6lRYx15eUXNW11dXYmIiKBhw4Z07tyZRo0a0bBhQxo1akSDBg1s+pXf2pQluI8DvYUQ13U7hBBnCvqZbQSqRHDr1TAyK5iTCN6mUwAL6ylkJetXf/cq1vnBv1XR/n//hcBAaNq0dIO6d4ePP5bFS+xtfYuKVTCI4RYX3GpKSnYKn+/+nJm7Z5KUlcTdDe5mYf+FhOWFcfLYSaatm8aJEyc4efIkJ06cICEhodjVXYEXCQpaQufOzjRq9KKeqIaFhVltUsreKUtwXYuLrQ4hRKKiKFXSVzglO4VWQa0MD+RchwDThTRN18QtEN/Mi/r7d+yA228vO7LfvbtMgtixAwo67qpUb7y8vMgoXo5RJ7h6IlM9uJJ6hXc3vcuio4vI0mZRP6c+jY434vTC0/S70E+vx15wcDBNmjShX79+NG7cmEaNGhEa2ognn4wmNxcOHXoYH5+Hbfhu7I+yBDe3gscshtGQghAmC9foMCm4HnXka1ax5Ifr1+H4cbnWtixuvx1cXWVYQRXcW4KgoCCuXbtWtKMaeLjXr1/n6NGjHD16lOPHj3Pw3EFi3GNIbpgsVeEI8A8kZyZTu2ltunbtypNPPknTpk1p0qQJTZo0wc/PsLHryy/LJcp//KHmBxmjLMFtoyhKmpH9CmD15Get0JKabaQWbl4aaPNKjeHq2uwYCK6zG3gE6Web7SxYJnb77WUb5e0tl439/bc5b0GlGhAcHExMTEzRDm9vudm54AohuHjxYqGwFt8SEwsWG/mAc3dntG214ARR2ihG1B1Bt17daDqnKbVr18bchrHbt8O0aTB2rFzso2JIqYIrhLBpoCUtJw2BKHeWGZRSExcKss2KCe6uXXIpmLH1t8a4/XaYPRtycsrfalTF4QgODtb3cOVOuxHc/Px8zpw5Y1RYi/fQq1mzJpGRkQwYMIB6zetx0PcgvyX+Rp42jyfbPMkb3d6gYc2GFbIhI0N+QYyIgE8/tcz7qo6Up6dZlWM6rbcgrGyG4CZlJRke9AyFzPii32NjoUUL8PQ0z7DOnWVFsdhY6NTJvGtUHJbg4GDS0tLIzMwsWpIUEmITwU1MTCQuLo4DBw5w8OBBDh48yOHDh/Um9erWrUtkZCQjR44kMjKSFi1aEBkZSVBQEElZSUzbOY1Pd39K1tUsHmn1CG91e4smgaWUKzWDKVPg9GkZaVPnkk3jmIKbXXpaL4CLkwt+7n7GBdcrFG7sKvo9NrZ862q7dJGvO3eqgnsL0Lx5cwDi4uLopPv/rlMHrNg1OTc3l2PHjhWKqm67cqWo6UpwcDCtW7dm4sSJtGzZkhYtWtC8eXOjsdWU7BSm/jWVGbtmkJ6bzvCo4bzd7W0ia0dW2tY//4Q5c2DSJLjrrkoPV61xTMHNKb00o44AzwBuZN0wPOAZKr1kTQ7cSIXLlyE62nzD6taFevVkKEKl2tOhoK3srl27igQ3IkIWpC+ZsVhOhBBcuXLFQFiPHj1aWKHKzc2Nli1b0qdPH1q3bk3r1q1p1aoVwbrJu1LIzs9mzp45fPDPByRnJzMkcghTu08lKsiwNVVFSE2FUaPkasoPPrDIkNUaBxXcskMKAIFegdzINCK4XsUKkR84JX8uj+CCDCvsNFGTQaVaUa9ePaKioli6dCnPPfecnERq0ACys2VYISTErHEyMzM5cuSIgbjeuFH0jNarV4/WrVvTv3//QnFt0qRJYREdc9FoNfx48Efe2voW8Wnx3Nv4Xt6/+33a1mlbrnFKQwgYPx4uXpQTZuZG5G5lHFRwE8HZQ9a2LYVAz0DTMVyQKxViY+XPbcrZfrlLF1nI5vJl6fGqVFsURWHs2LE8++yzfPfddzz55JNF/e3OnjUQXI1Gw+nTp4mLiyMuLo5Dhw4RFxfHqVOn0Gq1gFzb26pVKwYPHqzntdasaaQ7dTkQQrDh1AambJ5C3LU4OtbtyKKBi+jRoEelxjXGDz/IEtH/+5/0P1TKxjEFV9cevYyvcoFegcZLNBZ6uAWCW6+ezDIrD7onbNcuGDy4fNeqOBzjxo3jp59+YtSoUaxevZq7W7SgNpD33Xek7d5NfHw8Z8+e5cyZMxw/fpzs7GxAinXjxo2JiopixIgRheLasGFDnErWXK4key7t4ZVNr7Dt/DYaBzRmxdAVDG0x1OxlXeXh1ClZK6FbNzlhpmIedi24NzJv4KQ4Ga+FW0b8FiDAI6D0kILOwy1vOAFkY0k3NxlWUAW32uPi4sL69ev56KOPWLhwIet1DUW/+goAT09PGjRoQIMGDejZsyetWrUiKiqKFi1aWL3YysW0i0zZPIXFcYsJ8g7ii35f8HS7p3F1tk4yaG4ujBgh839+/NF0cT0VQ+xbcLNuUNOjpmFv+pzrZcZvQXq4qTmp5Gvz9evpuvqDsyfcPC8zzIYMKb9x7u5SqHfvLv+1Kg6Jp6cn77zzDlOnTiUpKYmkpk1x69UL37lzqVmzplU8ydLIysvis38/46MdH6HRanjjzjd4teur+Lpbd13W22/D3r2wapX8cqhiPnYtuNczr1PLy4gnm50IvmWvGwz0lGGC5KxkansXE2hdq52rR0GjqZiHCzKs8M03kJ9PVp4L338vv2rddZesVV7Ff38qVYSiKAQGBhLYujVcuCCL1lchQghWHVnFy5te5nzqeYa2GMonvT6hQc0GVr/3pk3wySfw9NMV81NudSwbRLIwN7JuGBQ1BsqshatDd63RpWFeoZByVv5cUcHt1AkyM0n+5xBdusC4cTBzpuzQ89BD8quXSjUmKgoOH5bT9VXEwasH6b6oOw+uehB/D3/+euIvVg5bWSVie/EiPPywzBGaYfU6gdUT+xbczBuFXmohmmzITzcrhqu71ngcNwzyrkKNGkUzzuWlYE3md+N3c+wY/PILZGXJCvcrVkgvQKUa07KlLNMZH1/2uZUkPTedlze+TLv57Th87TDz7ptHzJgY7oqomkyD3FzZkyw7G1avlqUkVMqPXYcUbmTdoF2ddvo7y2iPXpwAz4DCcQzwrg+u6RDd1bBDr7k0bEi2by38j+9i9vyxhb0np0yRwvvuu9C7t9pzstoSVZA8EBcH4eFWu80vx39h4m8TiU+L56m2T/Fx748Ln+2q4pVX5Pzw8uXQrFmV3rpaYdce7vXM64YebnbZhWt06EIKxtfi1gNnAR1LabZUBnn5CjvyOnGXx26eekr/2FtvyRo3kyZBspH6OSrVgFatZKB+3z6rDB+fGs/AZQN5YNkD+Hn4sf3J7Xw94OsqF9uVK2HWLHjuOet03r2VsFvBzczLJDs/2zCGa2ZaL5QRUkgraP7X0rwsIWOsXg1/ZXeiQc4xnG6W7G0PX3whxfaddyp8CxV7xs8PIiMtnnGoFVrm7Z1Hiy9bsPH0Rj7u9TH7xuyja3hXi97HHI4fl6m7nTurVcAsgd0Krk4kDVYplMPDreFeA2fF2XhI4WxBmd+IigejZs2Ci3U7oQgB//1ncDw6GkaPhi+/lJPZKtWQLl1k8ktBBlllOZdyjt4/9Gb8+vF0DuvMkQlHeKXrK1ZbU1saqakwaJBcAblixS3XoNgq2K3gXs+UsVqDkEI5YriKohDgGWA8pHDwsnwNqNgM84kT8u+s/fjb5A4ThWzefFNOYqveQTWlSxf5NebEiUoNI4Rg3t55tJrbij2X9jC//3w2PrqRCP8Iy9hZTjQaeOQR2b1h5Up1va2lsFvB1XmlRkMKihO4mZdzHugVaNzD3X8EMp0h55LhMTNYulSG7wY96Q/Nm5tMgAgPhyeegK+/hmKV9VSqC10LvuZv21bhIS7fvEzfH/sWerWHxh9iTPsxVZ5IUZw33pDF0GbPhh6WL8Nwy2K/glsQUjD0cBPBLVCKrhkEepqoGBYbC3l+kHG+3LYJIQW3W7eCLtmdOknBNbEe87XXIC9PXbtYLWnWDOrXhw0bKnT5ryd+pc28NuyI38Hc++ay8dGN1Pc33Y26Kli8WDamHjdOVgNTsRx2K7i6kILRGK4Z4QQdtbxqFY6Vn5/PlStXuBwXR+6lS+Bet0KCe+SInEwYPrxgR+fOkJgI584ZPb9RI1lu4euvZSsSlWqEokC/frB5s2y5ZCbZ+dk8v+F57l96P6G+ocSMiWFch3E29WpBTkWMHi2zJWfNsqkp1RK7FVxdGMBgCYyZdRR01HKrxflt57n77rvx9vambt26hLZujQ/Q9dPLTF9xihvXDTrBl8pvv8lX3brbwq4PpRQkf/55SEmRJe1Uqhn9+slP0j//NOv000mn6fxNZ2bvmc3znZ5n11O7aF6ruZWNLJv4eBg4UDazWLlSnSSzBvYruJk3qOFew3B21sy0XoAtW7bw03M/kbY8jcuXLzNhwgS+/PJL5g4ezCQgR/HgxR/yqB9Rn/fee4/MzEyzxt2wQa55Dwsr2NGqlay+XEohm65dZYGx2bOrNBNUpSro3Rtq1pSls8pgw8kNdPi6A/Fp8awbsY6Z98zEw8XqDbDLJCVFfm6kp8uMydrm+zQq5UEIYY2t0jy8+mHRcFZDwwMr/IT4b2Kp12q1WvH+++8LQNQOqy14GHH15tWiEx55RIiwMCEu/CQOfIgY3L+HAETz5s1FTExMqWOnpQnh6irEK6+UOHDHHUJ07lzqtYsWCQFC/PFHqaepSKz1bFrnuR03TghPT/mAGEGj1Yh3/3pXKFMVET0vWpxJOlOp21mSnBwh7r5bCBcXITZvtrU1Dk+pz5hde7gGE2b5WZCXCp51TF4nhOC5557jjTfeYMSIEUxbMw2aQmJmYtFJ+/dD27bgHUHrcFj9xTNs2rSJmzdv0rlzZxYsWGBy/C1b5ATYvfeWONC5sxy3lDje8OGyu7YaG6uGPPGEzOc2EjPKyM1g6IqhvP3X2zzS+hF2jNpRJcVmzEEIeOopGQ359tvy9VJVKT92K7hGSzNmF6yr8jAtuG+88QZz5sxh8uTJLF68mHq15ALCaxnX5AmZmXDsmMxK8G0k96WfplevXhw8eJAePXowevRoXn/99cJ2KMXZsEG2ge5aMumnUycptgcOmLTN3V3O/P72W6WXbarYG506yW3aNCho/ghwNf0qPRb14OfjPzOj7wy+H/g9Xq7WLUheHt5+W35GvPsuPP64ra25BSjLBa7gVmlCp4WKJ9c+qb/z2nYhFiPEpQ1Gr1m0aJEAxJgxY4RWqxVCCHH42mHBVMTSuKXypF275Pf6n36Sv68OEmLXU4Vj5ObmirFjxwpAjB8/Xmg0Gr17NGkixP33G7n5hQty3NmzS31fV67IkMTE0qMiKo4WUhBCPlMgxOLFQgghjiYeFQ1mNhCe//MUPx/7udLDW5q5c6W5o0cLUfDnolJ5HC+koBVarmZcJcSnRJ2DrAIP10hI4eDBg4wbN47u3bvzxRdfFC6vCfIOAop5uLqmkW0Lupf6NIabpwrHcXV1Ze7cubz66qvMnTuXiRMnFnq6V67IzJu7jFXECwuTzSTL6AAREiJr5X73nUydVKlGPPAAtG4Nr79OzNkddF3QlYy8DLaN3MaAZgNsbZ0eS5bAM8/IQvlz56rF8qsKuxTc5Kxk8rX5BHsH6x/QCa6HvhDn5OTw8MMP4+/vz7Jly3BxKao6GeAZgLPiXCS4+/fLGeX6BYvLfRtD+im98RRF4cMPPywU3ZdffhmAf/6Rx7t1M2K0osivlKUsDdPx/PNyNnjhwjJPVXEknJxg5kx2as5z96Ke1HCvwc7RO+kY2tHWlumxbp0MH9x1l1z+Vc4O7CqVwC7r4SakJwAY93AVZ4PEhw8++IDDhw+zfv16goP1RdpJcaK2d219wY2OLvpI92kMZ7+XE3IunoXX6UQ3MzOT6dOnU69ePU6enIS3d5FzbECnTrBmDVy/DrVMVzNr316Wbvz8c3j2WbUJX3Xi74bO9BvlQp3kHLb0mEV4zYa2NkmPrVth2DC5RPGXX+RqRpWqwy49XJOCm30FPIL10noPHTrEBx98wKOPPkq/fv2MjhfkHcTVjKtyMuPgQf2WOr4F9XDTzxhcpygKM2bMYNCgQUyePJl161bRtassvWgUXQLEnj1lvsfnn4czZ2S+ukr1YO/lvdy35D7qBTZg2x91CR/1gvzwtRP27IEBA6Bx46LJX5WqxS4F92rGVQCCfUqGFBL04rdCCCZNmkSNGjWYOXOmyfGCvYOlh3vihOwRUtxFLRTcU0avdXZ2ZvHixXTs2IX4+Edp1MiwDGMhHTrIr5VmhBUGDZJh39mzyzxVxQE4fv049y6+l0DPQDaP3Erdhavg0iWZjmhmQo01+e8/6NsXgoJg40YINNIqUMX62KXglhpSKLYk7LfffmPLli1MnTqVwFKeoGCfYK7cvCLDCWBccG8aF1yQ7bEnTvwZqMOqVYNISEgwfqKPj/SedcHeUnB1hQkT5LreQ4fKPF3FjklIT6DPj31QUNj02CZCa4TKso2LF8tJ1Pvug7Q0m9m3Z09RMtzWrXJuV8U22K3guju74+fup38g+0qhh5uXl8dLL71E06ZNGTduXKnjhfqGcvnmZcS+fXIxbPNieetuNcEtAG6WvjB2//5auLquJT09iaFDh5JrqiVvjx6yA0B2dpnv8+mnZQxN9XIdl1xNLkNXDCUxI5ENj2ygSWCTooNDhsh03+3boXt3k8WNrMnu3VJsAwPhr7+s2npNxQzsVnBDfEL0Kydp82WlsALB/fHHHzl27Bgff/wxrmVMs4b6hpKnzSPvv11y2U7JIKxfJKQeLXWMv/+G229vw8KFC9mxYwfPP/+88RO7d5cJEGa0XQkMlA0mf/gBbhipIKli/zy34Tl2xO9g4QMLaV+3veEJDz8MP/8sA/Zt28JPP1WZbbt2QZ8+si6CKrb2gV0K7tWMq4bx2+yrgADPOuTn5/P+++/Tvn17HnjggTLHC60RipMWnPfFwm23GZ5QowWkHjZZVSYtTUYjunWD4cOH8+qrrzJv3jy+/fZbw5PvvFPGcf/6q+w3ipw8y85W030dkWWHljE/Zj6v3P4Kw6OGmz6xXz/ZaLJRI+n1DhwI58tfFrQ8/POPvtiqHRvsA7sU3Cs3rxjGbzMLOjN41mXp0qWcPn2at99+26z6oaG+oUQmgnNGZtFKguL4tYTcJMi+ZvT6f/+VLat062/ff/99evXqxcSJEzlQMpXXz0+u+9q6tUy7AFq2lLVyZ89WEyEciYtpFxm/fjydQjvxfs/3y76gYUP5refjj2HTJmjaVAbx4+Mtbtuvv0qxDQ2VYltY1U7F5til4F5Mu0iYb4mnJFN6BBqPMP73v//Rpk0b7i8sSFs6oTVC6aTrpGNMcP1byte0I0av37ZNRiG6dJG/61YuBAQEMHToUNJKToh07y6DZ2bOTr/xhhTbzz8363QVGyOEYPQvo8nV5PLDoB9wcTJzOburK7zyiqzl8eSTsiJ9w4ZyYezmzRZpRPnjj9KBjoqSXq4qtvaF3QnuzZybJGcnE+5XIuCUIdvertkcx4kTJ3jrrbfMro4f4hNC50uQ5eMBTZoYnlCjhXxNOWz0+r//lk6rd7EGv0FBQSxfvpyzZ88yevRoRPFwRI8ekJtrdvvsdu3kRPaMGXDzplmXqNiQn47+xMbTG/mw54f6k2TmUq8ezJsHp07JmNLWrXJmKyICJk2SSqnRlHvY2bPhscdkBtmff5aae6NiK8oqtlDBrcLois0sObhE/8B/zwmx3Fd0795dREREiPz8/HKNe6iui4iLDjV+UKsVYoW/EHvGGxzKzJTFZl5+2filn3zyiQDErFmzinaaLJprGl1NnQ8+MPuS6o5dFq/JyM0Q4TPCRasvW4k8TZ4l3qcQWVmy4M399wvh7i4fhIAAIYYMEeKLL4Q4erTU6jJarRBvvikvGzRIDqdiM0p9xuxOcH8/+btgKmL7+e36B7Y9II7MbSwA8eGHH5Zv0PR0ka8gfhxopKC5jo1dhdh4p8HurVvlv9K6dcYv02q14v777xeurq5i165dRQfuvluIqKhymXn//ULUqCHEtWvluqy6YpeCO+3faYKpiK1nt1rgLRohLU2IpUuFGDlSiHr15MMHQvj7C9GzpxCvvirEypVCnDghRF6eyM6W9fRBiKeeEiLPQp8BKhXGsQT3q71fCaYizqec1z/wW1vx3OAI4ebmJq5evWr8YlNs2yYEiPHjw02fs+cZIVbUEEKrX47xnXeEUBQhkpJMX5qUlCQiIiJEeHi4uH79utz52Wfyn/f8edMXluDIESGcndXSjQXYneBm5WWJkM9CRM9FPS31HktHqxXi1CkhvvlGiLFjhWjfXn5zKhDhRJcQcYfnXvnNqPsfQrvwO+khnDkjRG5u1dioUpJSnzG7K15zIfUCTooTdX3102Eybpxj0cYMhg4dSlBQUPkG/ftvtAqsD0xCCGE89hvQHk5+CTdPQ40mxS+ldWuZpWOKmjVrsnLlSrp27crjjz/OunXrcLrvPnjpJVltvIzEDB2RkTIZYt48WdSmadPyvU0V6/L9ge9JSE9gyeAlVXNDRZFLyRo1kq10Qa7xjovj5Obz9PvkLuJTa7CszgsM/+dz+Eujf23durIqXni4rAsaHGy4LmyMMAAAGGFJREFUBQWp3SKrELsT3Pi0eEJ9Q/VnfvPSWfZ3MqnpMH78+PIPum0bNxrV4YLLFZKykgj0MpIGHNBOvibFFApuXp6c99I966XRoUMHZsyYwYQJE/joo494/bXXoEEDWZ3GTMEFmDpVZoROnAh//KHWKbUnvor5ijbBbege0d12Rri7s/VmB4Z+2gEnZ/jzH7j99hmQ+zFcuCDX954/r//zf//B1auyJqgx/P2lR6HbSv6u2/z8ZPp68c3XV6ZLOtnd/LtdYneCeyH1guEKhcx45m6GqGZhdDXobVMGeXnw77+kDe4BrOdsylnjguvXEpzcIDkGIh4C5Fr1zEyZy2AO48ePZ/v27bz11lt07NiR3vfdJxtFZWWZXQcvOBg++kgu0Vy8WGaiqdieAwkHiLkSw+x7Zpu9OsbSCCETZF56SX77WbdOOr+A9FIbN5abKTIzpfDqtoQE+ZqYCMnJRdulS0U/m0phL46iyCU8JcXY21s+9x4eRa/l+dndXb4vV1fTr66uDiX2dim4JQs2/7djIzFn4YuPHiz/wx4TA5mZON/VHS6t52zyWTrU7WB4npMr+LeWHm4Bf/8tX80VXEVR+Prrr4mLi2PEiBHsmzmT8DlzpKs6cKDZJo8bJ9dTvvAC3HOPurzHHlh1ZBVOihMjWo2wyf0zM2HMGPkhPGgQLFpUgfKKXl7yW1cDMxtYCiGdBZ343rwpveSSm6n9SUkyjVK3ZWUV/Vys71ulcXYuW5iLv7q4yGtcXIq2kr+Xdk5AgPzjrAB2Jbh5mjzOpZxjRJT+Q/3lt8vwdodHn3i6/INu2wZAYN9BsOBlziQb1r0tJKA9nF8GQguKE//8Iz2JkBDTl5TE29ub1atX07FjR4bOmsU/gYG4L19eLsF1cpJr4tu2lTHdn35SQwu2ZsOpDXQO62zY2LQKOHdOiuyBA/D++zBlShU5dYoiRdrLS6atWZL8fONCXPLnvDy55eZa7jU3V95fo5GvxbeS+4z9Hh5eYcG1q1UKx68fF0xFLIpdVLjvxo0bwsPdWYzt5VKxTnf33CNEZKQQQojAjwPF2HVjTZ97+jvZpDL5oNBo5Eqcp54yfXpp/PTTTwIQYyMjhfDyEiI9vdxjTJ8uJ6TnzKmYDQ6O3axSuJp+VTAV8d629yz9Hsvk55/lklw/PyHWr6/y26sYo0Rj2RKU+ozZVfDjxA1ZIrFpYNH0/KJFi8jO0TB+YKPyu3lZWTKZvE8fABrUbMDZlLOmzw8qiB1c+4cDByAlxUT/MjMYNGgQU6ZMYf7RoyzMzKxQa4dJk2QG2osvwt69FbNDpfLsvigbg1blZFlOjvz/f+ABmYC2d6+sgaNiB1Ti64VdCe7JGyeBIsEVQjBv3jxub+5BmzZtyj/gX3/Jryb33gtA44DGhfcwincD8KwLif+waZPc1bNn+W+r47333qPn3XczHtj35Zflvl5RZHffOnVk44ALFypui0rFiU2IRUGhTXAFnsEKcOqU7Hk3axY895wsnlTaXJiK42BXgnvixgkCPAMI8AwA4M8//+TEiROMvzsHfCuwKHXDBjnbWdDXPLJWJOdSzpGZZ6KojKJAUDe49g+bNwtatqxcdXwXFxeWLltGkK8vD2zbxhVdx4lyUKuWdI4zM6F/fzkPoVK1HLh6gMYBjfF1t24TMCHkpFi7dnD2LKxdK0XX3d2qt1WpQuxLcJNO6IUT5s6dS2CAP0M7CqhRQcG9+265zAQpuALB8evHTV9T+07IusTFY2fo3bv8tzQYrnZtfvnxR5KBAfffT2YF+lu1aAGrV8Px47LGiSq6VcuF1As0tHL33cREWTTs0Udlok1srAwnqFQv7EZwhRDEXY0jslYkAJcvX2bt2rWMerAHHm6Ab7PyDXj4sPxudt99hbsia8uxj14vpbtDiFTZu5tvoFev8t3SFNEDBrCkdWtiLl3iiccfR1uBMny9ekmP59AhGea4dKnsa1Qsw+Wbl6njW6fsEyvIL7/Icorr1sk12Nu2qd0Zqit2I7gJ6QkkZibSNkQ2ePz666/RaDSMHVDgWZTXw12+XAa3Bw8u3NUkoAnOijNHE0sR3BpNSMxpQv9263WRCIsw4M03+RRYtXo1b731VoXGuPde+cd56hR07KhOpFUFWqHlasZV6vhYXnBTU2VZ3AcekKGrvXvh1Vflck+V6ondCG5sQiwA0SHR5Ofn8/XXX9O3b18a+cbLySw3f/MHEwJWrJCFwIOLWvW4u7jTKKARhxON173V8fuB+7i7xVZ8PDIq8laMM3Agk+vX56natfnggw+YO3duhYbp21dOori5QdeuMH26RepWq5jgeuZ18rX5FhfcNWtkqOiHH+DNN2W9+latLHoLFTvEbgR3f4KcUGod3Jp169Zx6dIlnnnmGUjeDzWjyzfYwYMy4PnggwaHWge3LhR3Y5w+Dd9v6YebSw4kbC7ffUvD1RXl1Vf5MjGR/p07M2HCBJYtW1ahoVq1kunx994rl4z16CHfsorlcVacef2O1+lSr4tFxrt4UebADB4s+43t3AnvvafWj7llKGuhbgW3cjNsxTDRcJasV9u7d29Rr149kZ+VLMRiRYiD75ZvsOefl2XsEhMNDn28/WPBVMT1jOtGL502TQhX5xyRvzxQiO0Plft9lEpWlhAhISKza1dx5513ChcXF7Fhw4YKD6fVCvHtt3JhvJOTEE8/LcTp0xa017bYTeKDJcjPF2LWLCF8fITw9BTi00/V2rXVFPtPfBBCsCN+Bx3rduTkyZNs2rSJMWPG4HzzMCDK5+FmZclE8yFDjBYh0NVRiLkSY3AM5MRUZEs3nBsMh4trIS/N6HkVwsMD3n4bzx07WDdmDFFRUQwePJi/zOzwWxJFgVGj4ORJWV1s0SLZQWj4cNmnsAJdWlSswO7d0Lmz7KZzxx1yPvell2RavsotRlmKXMGtXJy6cUowFfHlni/Fiy++KFxcXMSVK1eEODZbptqmXzB/sO++k/mwW7caPZyclSyYinj/7/cNjl29Kj3Ft94SQlzbIe99emF5307p5OUJ0bKlEA0biqsXLojIyEjh6ekpNm7cWOmhL12SDQH8/OQ/QZ06spj5unVC3LxpAdurFof3cC9dEuKxx4r+L5YsqVh2uopDUeozZhcP7rf7vhVMRcScjxEBAQFi2LBh8sA/DwrxU6j5T6lGI8WsZctSr2n6eVPRf0l/g/2zZsl/kbg4Ia//pakQGzpY/q9k0yZ5oylTxNWrV0Xr1q2Fu7u7+PXXXy0yfFaW7MIyYID8+goywtKhg2wcMH++EH//LUR8fFlp4TbFYQU3K0uIDz8UwttbCDc3IV57zSE/8FQqRqnPmF18qdl6biu1vGqxb+M+kpKSmDhxolxpcG0bBPc0v4bCzz/L72tLlpR6Tff63Vl2eBn52ny9QueLFsksn6goAAWaT4L/noHEHRB0R+XeZHF69YKnnoJPPiFowAC2bt1Knz59GDRoEAsWLODRShbB9fCAoUPllp0NO3bIEMN//8nVcvPnF53r7i6byNauLSMwtWpBYKAsZ1q8lKnuZ2dn+U/r5GR6EwU9YLRaw59L7uvdW96rOiCEXLY3eTKcOSOXe02bVqxmrYpKWYpcwc1scvNzhf9H/uKxnx4T0dHRolWrVkKr1QqRekx+pT8537yB8vJk08bGjeUMRSksP7RcMBWxM35n4b64OCkDM2cWHzNdiJUBQmw19IYrTWqqEPXry0aBCQkiJSVF9OjRQwBi6tSp8t/ACmi1suXVH38IMXeu7EY8fLjsTxgdLURYWJFXXBXbiRMmTXUoD3fvXiG6d5fvKTJSCAtEiFQcE/v2cP869xcp2Sm0yGrBD7E/MH/+fFlkPGGLPCGou3kDzZ8v07BWry5z5XiPiB4AbD6zmc5h/9/e2QdXWd15/PO7uTcJCSSB8Opqedug0sUyUalYVqRayoiFrtoWu9UdZQaq9W3ozNYOu5AK3cFpYTpYa9FFl65Q0bZLESHo4tJCuwIS3jG2SQGDwqZALhBIyMv97R/nud4Q8nITcu99En6fmTPPk/OcPOcl3/zuuefld24BnP/ZYBDub+qKN5gN138X9syFyi0xb2JdQU6Oc3Q7YQLccw+5mzZRXFzMrFmzKCoqoqysjGXLlpGVldV1eeJ6p/H4oI5EYi5Jo6G21sVHQ7SX2jQ0NsZ6wCIX37cU1913VB05AnPnOh8I/fvDc8/B7NnO17VhXEJ7FrmTIW4eWvOQZv8wW++9717Ny8vT6qjf2E13qr55bXzjp0eOuJmiL34x7vHWcS+N08JlhaqqWlXlxtseeKCFhPXnVP/rajeW25iAdTyvv+66RV/+smpNjUYiEV2wYIECOmbMGC0tLe36PLsHvu7hhsNugjIjQzUz043ThsOXU12jh9CmxlIq3L+e+6tmLMjQb738LQ0Ggzpnzhz3oPak6qo01V1Pt/+S+nrV225zCxw7sAh18R8XK0Xon0/+WRctci1RUtJK4sOvueGN/f8W9/s7xPLl7iz2SZNUvWPWN2zYoPn5+Zqdna0rVqxI2BCDj/Glwa2rU126VDU/32nmgQfc571hePjX4M57d55ShH5z5jc1FArpRx95y7/+9IIzcCd2tP2CSMSt9gfVX/wi3mxVVbXidIVShH5vw3zt10918uR28tnyddVfhlT/73cdyiduXn3VTWkPG6a6fbsrY0WFTpgwQQGdOnWqVlRUJCZvf+I7g7t5s+ro0U5ukyap7tzZVVU1ehD+NLiHqw5rr4W99CsvfkUzMzN15syZ7kEkorrus6rrC9seHmhsVH3qKVeFuXPjbYyLmLpyqmbNH6gEa3TXrnYSXzjlhjje6Kdata9T+bXLtm1u1ioQcHULh7WhoUGXLFmivXr10pycHF28eLHW1tYmJn9/4RuDW12t+vDDTmrDhrljb668LxxGnPjP4F5ouKC3Lr9Vs36YpbMfn62BQED/FJ2uPvpm+xsOTpxQ/epXXfGfeKLT6n+heJNShN7y+PPx/cKZMtXfDHFGt3Jrp/Jsl3BY9ZFHXN1yc1Xnz1c9flzLy8t1ypQpCuiIESN09erV2ujjRbRdgC8Mbk2N6s03uxGfp59WPXcuEVU1ehD+MriV1ZU6+T8nK0Xocxuf04yMDH3wwQfdw4YLrhe5tsDdN6e+3o13DhniVvIvWdJpY3vqlOqoayOaPmui9l3Ur1XfCpdwtlz1tyPdGPO+hS2XsysoKYl9qASDbhfDypW68fXXdcyYMQroqFGj9KWXXuqpPV5fGNziYtf8v/pVIqpo9ED8YXDPXjirc4rnaN6iPA0+E9TlJct1xowZmpmZGRub3P4d17s92uR40khEtbRU9Zln3Pc5UB03ro0ZrvY5e1Z14kRns19+a6+GngnpXSvv0obGttfvfsqFKjemuxJnfMtfUa0/3+nytMnBg26x7ODBru4i2lBYqKsnT9bCoUMV0P75+frkk0/q7t27E1OG1OALg6uqeuhQF9fM6Mm0qTFR1YSsNmse0RhpZNRPRzF28FgWTFpA+XvlTJs2jXnz5vGDef8KO74Lf1kKufdDzT3wwQdw8CBs2RI73uCOO5yXlunTO36Cr0dZmXPusns3vPqqW3f7wo4XeHT9o9z/d/fzyvRXyAjGcYiUKhzbCLu/B+G9kN4PrvkHuGoqDLod0vt2qnyt0tgIO3dCcTG8+y6UlKBnz7IJeBH4LVAHXN+nD3ePHMndY8cy/sYbCQ0c6LZyRUN2dmwbWSjkFh+HQrEtZP4hFYVJyD+DcUXRpm6TZnAB9G9HIufOc+jCBT4fDjNEYMfNAdLvboQRwDvAiia/PXw43HSTM7RTpsDQoZ0u0EcfwfPPw9KlzvfoqlWx03dUlWf/8Czf3/R9Rg8YzcJJC5l27TTSAnG43leFys1Q9iJ8sgHqT7v43iOgbyH0KYDsoV74DKTnO2fqaZd5MmAk4tyE7d0L5eWcPHCA17ZtY83HH/O78+epB7KAm4Hx3vV6YCTQquvVtLSLjXD0OOioIY7n2tozcG3V9AruA3VEi+eFmcE1uiM+MbiqsHA87x+v4uuvHSJc28Af5wvXXR2Buhyovhv6fAkGDEQHD0ELRqFZ2WgL+/Jbi6uthXDYhaoq50y8tBS2boU9e5z9+MY34Mc/bvk03uKyYh5961EOhQ/RP6s/dwy/gxsG3UBBvwIG9R5EbkYuuZm55GTkkCYxYyxRoxKpJ3hyO2mnthOo2k1aeC9y/giil/pJ1EAmpOeioTwI9kYD6c4IBzIgEIK0TAiku5CW4a4EgJhR+zRfYlcR4Ux1He/87yF+v+Mw7+07xq6yEzQ0uj9JME0YMSCLoX0zuSo3naty0xnSJ0S/XmnkpAfIyQi4a0jITBNC0SCxayAgF/+FtflVL/65WRljF4GHV8DVLZ5XZwbX6I74w+DW1dVxzYAMKs/AwJwARfeOpb7xNt7ZP5m3932Jhsbgp4azK+nTB2680R1NM2MGDBvWdvqGSANrP1zLmtI1bD68mYozFZeVfwC4KghDg3BNEPqmQV7Aheh97wCkC2R4oaX7dLwOpPfeeK5Re1xbBx98DKUfw4fH4MNP4OgpOBZ2oaGDfnMD0qTzG82vece2lfjmz9a9/R9M/Pt/aikbM7hGdyT5BldEioFLvX87+gMnujzT7oe1Q4yW2uKEqk5JZiFMt+1ibeBoqx3a1G2ieritIiLvq+pNSc3Uh1g7xOgObdEdyphorA0cl9MOvjhixzAM40rADK5hGEaSSIXBfTEFefoRa4cY3aEtukMZE421gaPT7ZD0MVzDMIwrFRtSMAzDSBIJMbgiMkVEPhSRMhF5uoXnIiJLved7RaQwEeXwA3G0xe0iclpEdnthXirKmWhE5GURqRSR/a08T7kmTLcO02wC9dqes4WOBiANKMdt1k0H9gCjm6W5C9iAWyR8C7Ctq8vhhxBnW9wOrEt1WZPQFrcBhcD+Vp6nVBOm2w61Q4/XbKL0moge7jigTFX/oqp1wGvA9GZppgPRIxreA/JEZEgCypJq4mmLKwJV/T1wqo0kqdaE6dZhmiVxek2Ewf0boOl+2KNeXEfT9ATired4EdkjIhtE5LPJKZrvSLUmTLcO02x8dEoLiTgmvaW9xM2XQsSTpicQTz1LgKGqWi0idwFrgIKEl8x/pFoTpluHaTY+OqWFRPRwjwLXNPn5auCTTqTpCbRbT1U9o6rV3v16ICQire3n78mkWhOmW4dpNj46pYVEGNwdQIGIDBeRdGAGsLZZmrXAg95M3y3AaVU9loCypJp220JEBovnZ1FExuH+JieTXtLUk2pNmG4dptn46JQWunxIQVUbROQxYCNuxvNlVT0gIt/2nv8cWI+b5SsDzgMPdXU5/ECcbXEf8IiINAA1wAz1pkF7EiLyS9zsdn8ROQrMB0LgD02Ybh2mWUei9Go7zQzDMJKE7TQzDMNIEmZwDcMwkoQZXMMwjCRhBtcwDCNJmME1DMNIEonYadYjEZFGYB9uaUgDsAL4iapGUloww2gF06z/MIMbPzWqOhZARAYCq4Bc3Pq8y0JE0lS1g4eVG0a7mGZ9hg0pdAJVrQRmAY95O03SRORHIrLD8405G0BEAiLyMxE5ICLrRGS9iNznPTssIvNEZCvwNREZKSLFIrJTRLaIyHVeugEi8mvv3TtE5Ate/MQm/kh3iUifFDWH0Q0wzfqEVPud7C4BqG4hrgoYhBPyv3hxGcD7wHDcjpz1uA+2wV76+7x0h4F/bvKuTUCBd/954F3vfhUwwbv/DPCBd/8m8AXvvjcQTHUbWfBXMM36L9iQwuUR9Rg0Gbgh2hPAfW0rACYAb6gbMzsuIv/T7PdXA4hIb+BW4A1vizq4fwKAO4HRTeJzvJ7BH4AlIrIS+I2qHu3Smhk9FdNsCjGD20lEZATQCFTiRPy4qm5slmZqO685510DQFi98bZmBIDxqlrTLH6RiLyF28/9nojcqaqlHa2HceVgmk09NobbCURkAPBz4Kfqvh9txDnzCHnPR4lINrAVuNcbFxuEc4ZxCap6BjgkIl/zfl9E5HPe47eBx5rkHZ0EGamq+1T1WdzXwesSUFWjh2Ca9QdmcOOnlzfYfwD4b5yofuA9+3fgIFAi7tC5ZbhvD7/G+c2Mxm0DTrfy/n8EZorIHuAAsWNNngBu8iY2DgLf9uKfEpH9Xvoa3PlKhtEU06zPMG9hCUZEeqvzjJ8PbMdNGhxPdbkMozVMs4nDxnATzzoRycOdgLrAhGt0A0yzCcJ6uIZhGEnCxnANwzCShBlcwzCMJGEG1zAMI0mYwTUMw0gSZnANwzCShBlcwzCMJPH/uAsdROJ1/aEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 410.976x246.586 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2)\n",
    "\n",
    "# left plot\n",
    "\n",
    "# cost, alpha, beta, lamb, polarity\n",
    "par_sets = np.array([\n",
    "    [2, 3, 7, 4, 1],\n",
    "    [2, 7, 3, 4, 0],\n",
    "    [0, 3, 7, 4, 1],\n",
    "    [0, 7, 3, 4, 0]\n",
    "])\n",
    "\n",
    "colors = [\"blue\", \"green\", \"red\", \"orange\"]\n",
    "labels = [\"$Pr(\\Theta_{+} \\mid c=2)$\", \n",
    "          \"$Pr(\\Theta_{-} \\mid c=2)$\",\n",
    "          \"$Pr(\\Theta_{+} \\mid c=0)$\",\n",
    "          \"$Pr(\\Theta_{-} \\mid c=0)$\"]\n",
    "\n",
    "for index, parameter_set in enumerate(par_sets):\n",
    "    cost, alpha, beta, lamb, polarity = parameter_set\n",
    "    func = lambda theta: unnorm_pr(theta, alpha, beta, cost, lamb)\n",
    "    # xs = np.linspace(0,1,1000)[1:-1]\n",
    "    # ys_grid = [func(x) for x in xs]\n",
    "    norm_const = integrate.quad(func, 0.0, 0.996, limit=1000, points=(0,1))[0]\n",
    "\n",
    "    xs = np.linspace(0,0.99,1000)\n",
    "    ys = np.array([func(x) for x in xs]) / norm_const\n",
    "    if polarity:\n",
    "        ax1.plot(xs, ys, color=colors[index], label=labels[index])\n",
    "    else:\n",
    "        ax1.plot(np.flip(xs), ys, color=colors[index], label=labels[index])\n",
    "\n",
    "ax1.plot(xs, stats.beta.pdf(xs, 3, 7), color=\"black\", label=\"Prior\")\n",
    "ax1.set_ylim(0,10)\n",
    "\n",
    "ax1.legend(frameon=False)\n",
    "\n",
    "# right plot\n",
    "\n",
    "cost, lamb = 0.2, 4\n",
    "alphas, colors = [0.1, 0.3, 0.4, 0.7], [\"red\", \"black\", \"green\", \"blue\"]\n",
    "for index in range(len(alphas)):\n",
    "    func = lambda theta: unnorm_pr_fixed_beta(theta, alphas[index], cost, lamb)\n",
    "    xs = np.linspace(0.001,1,1000)[1:-1]\n",
    "    ys_grid = [func(x) for x in xs]\n",
    "#     print(ys_grid)\n",
    "    norm_const = integrate.simps(ys_grid, xs)\n",
    "    \n",
    "    xs = np.linspace(0,1,1000)\n",
    "    ys = np.array([func(x) for x in xs]) / norm_const\n",
    "    ax2.plot(xs, ys, color=colors[index], label=str(alphas[index]))\n",
    "\n",
    "leg2 = ax2.legend(title=r\"$\\beta=1, \\alpha=$\", frameon=False)\n",
    "move_legend(leg2, ax2, 0.1)\n",
    "ax2.set_ylim(0,2.5)\n",
    "\n",
    "ax1.set_ylabel(\"Density\")\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    ax.set_xlabel(\"Degrees\")\n",
    "\n",
    "    for text in [\"top\", \"left\", \"right\"]:\n",
    "        ax.spines[text].set_visible(False)\n",
    "    ax.set_yticks([])\n",
    "\n",
    "w = 5.708\n",
    "fig.set_size_inches(w, w*0.6)\n",
    "plt.savefig(\"qingFranke2014.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beta distribution with $\\beta=1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Various useful bits:\n",
    "\n",
    "$$\n",
    "\\int_{0}^{\\Theta} \\phi (h)^2 dh = \\int_0^\\Theta \\alpha^2 x^{2\\alpha-2} dx = \\alpha^2 \\frac{\\Theta^{2\\alpha-1}}{2\\alpha-1} \\\\\n",
    "\\int_{\\Theta}^{1} \\phi (h)^2 dh = \\int_{\\Theta}^{1} \\alpha^2 x^{2\\alpha-2} dx = \\frac{\\alpha^2}{2\\alpha-1} - \\alpha^2 \\frac{\\Theta^{2\\alpha-1}}{2\\alpha-1} = \\frac{\\alpha^2 \\left( 1 - \\Theta^{2\\alpha-1} \\right)}{2\\alpha-1}\\\\\n",
    "\\int_{0}^{\\Theta} \\phi (h) dh = \\int_0^\\Theta \\alpha x^{\\alpha-1} dx = \\Theta^\\alpha \\\\\n",
    "\\int_{\\Theta}^{1} \\phi (h) dh = 1-\\int_0^\\Theta \\alpha x^{\\alpha-1} dx = 1-\\Theta^\\alpha \n",
    "$$\n",
    "\n",
    "Finally, ES:\n",
    "$$\n",
    "ES(\\Theta) = \\alpha^2 \\frac{\\Theta^{2\\alpha-1}}{2\\alpha-1} + \\frac{1}{1-\\Theta^\\alpha}\\frac{\\alpha^2 \\left( 1 - \\Theta^{2\\alpha-1} \\right)}{2\\alpha-1} - c (1 - \\Theta^\\alpha)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [],
   "source": [
    "def U_fixed_beta(theta, alpha, cost):\n",
    "    first_term = np.square(alpha)*(theta**(2*alpha-1)/(2*alpha-1))\n",
    "    second_term = (1/(1-theta**alpha))*((np.square(alpha)*(1-theta**(2*alpha-1)))/(2*alpha-1))\n",
    "    third_term = cost*(1-theta**alpha)\n",
    "    return first_term + second_term - third_term\n",
    "\n",
    "def unnorm_pr_fixed_beta(theta, alpha, cost, lamb):\n",
    "    return np.exp(lamb * U_fixed_beta(theta, alpha, cost))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \"\"\"\n",
      "C:\\Users\\s1569804\\AppData\\Local\\Continuum\\anaconda3\\envs\\computational_models\\lib\\site-packages\\ipykernel_launcher.py:3: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAADYCAYAAACnQFtWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3iUxdbAf5OEEkBaCKGEFqpIlSZFpIOgiCKoXFBE8VqvXL1eRUVFLCh8iiiC9YqooHQQBIEgCIL0bgIJBAghBAgkEDBlc74/ZgMJpGw22zO/55nn3X3LzHnfnTl73jMzZ5SIYDAYDAbn4+duAQwGg6G4YBSuwWAwuAijcA0Gg8FFGIVrMBgMLsIoXIPBYHARRuEaDAaDizAK12AwGFyEUbjXoJQappTappS6qJQ6qZT6RSnVpQj5xSiletl5bSul1Hal1CXrtlU+505WSh1SSl1QSkUopR60V2aD7+Ap9Vkp1UgptVgpdVoplaiUWqmUamyvHN6KUbjZUEo9B0wB3gFCgNrAp8BdbpClJLAY+A6oBMwEFlv350YKcCdQAXgI+Egp1ckVsho8E0+qz0BFYAnQ2CrLFnT9Ll6IiEl6tl0F4CIwJI/jpdCVN86apgClrMeqAD8D54FE4Hf0n9ksIBO4bM37v4WQpw9wAlDZ9h0D+tl4/RLgeXc/V5PckzytPudSfmVAgCB3PytXJmPhXqUjUBpYmMfxV4BbgFZAS6A98Kr12PNALBCM/vd+GRARGYFWkneKSDkReR9AKXU+n/SSNc+bgD1irZ1W9lj354tSKhBoB+y39eYNPoen1edr6QrEi8jZot+q9xDgbgE8iCDgjIhk5HH8H8AzIpIAoJQaD3wGjAPSgepAHRGJQlsEeSIiFW2QpxyQdM2+JOAGG66dAewGVtpwrsE38bT6fAWlVCgwDXiuMNf5AsbCvcpZoIpSKq8/oRrA0Wzfj1r3AUwCooBflVKH8/lXLwwXgfLX7CsPXMjvIqXUJKAZMPQa69hQvPC0+gyAUioY+BX4VERmOypfb8Eo3KtsAv4GBuVxPA6ok+17bes+ROSCiDwvImHojqvnlFI9reddp/SsPcZ5pZetp+0HWiilVLZLW5CPm8BqpdwO9BGR5ALv2ODLeFp9RilVCa1sl4jI20W+Qy/EuBSsiEiSUuo1YJpSKgNdMdKBXkB3YDbwqlJqK7rSvYYeQYBS6g4gAogGkgGLNQGcAsKuKaucDSL9Zs3jX0qpGcBo6/7w3E5WSo0FhgFdi5tfzHA9nlaflVLl0S6ujSLiMIvZ63B3r52nJbRvaxt6mFU8sAzohO6AmAqctKapQGnrNf8GYqzXxALjsuV3F7qj4Tzwn0LK0hrYju4V3gG0vkbO/dm+C5CKdkVkpZfd/TxNcm/ylPqMHqoo1jyz19Ha7n5GrkzK+jAMBoPB4GSMD9dgMBhchFG4BoPB4CKMwjUYDAYXYRSuwWAwuAijcA0Gg8FFOEvhSp4pNlZQSvjyy7zPMckk95C7LHffLbRs6e7nYZJ3pHxxvYXrZy0yM9PlRRsMBoM7MQrXYLAFM17d4ACMwjUYCiJHOAuDwX6MwjUYbMFYuAYHYBSuwVAQxsI1OAijcA0Gg8FFGIVrMNiCcSkYHIBRuAZDQRiXgsHKSy9Bjx72X28UrsFgC8bCNQB798L58/ZfbxSuwVAQxsI1WDlyBOrVs/96r1e4r732Gs2bN6dRo0Z8/vnnduUxatQoqlatSrNmzRwik8Hgbky7cDwiEBMDdevan4dXK9yVK1eyc+dOdu3axfz581m0aJFd+YwcOZIVK1YUWR6DD+NFLgXTLpxDfDxcvgxhYQWfmxderXCXLFnCyJEjSU9P55NPPmHw4MF25dO1a1cqV65s8/n79++nV69eNGrUiAkTJvDMM8+wdetWu8o2eAFe5lIw7cI5REfrbf369ufh+lV7syqvAxTu9u3badeuHUFBQdStW5cPP/wwx/Fbb72VCxcuXHfd5MmT6dWrl11l/v333wwZMoS5c+cSFhZGkyZNaNOmDe3atXNqub5GRAS89hpMmgR16rhbGhvwIgvXtAvncPiw3hZF4Tprdcr8UUrk1VcLPC0/LBaL1KxZU0RELl26JCNGjJAJEybYnd+RI0fkpptuKvC8pUuXyrBhw658r1atmmzYsMHucosrjzwiUrq0SEJCnqe4Y1XV3Bk6VKRJE3tv1aWYduE8XntNxM9PJDU139PyrWOut3AB/P2LbOFGRkbSsGFDAAIDA+ncuTPx8fE5znHGP+rOnTu5+eabAYiLi6NcuXJ07tzZ6eX6EnFxMGsWPPooBAe7Wxob8RIL17QL5xEdDbVqQcmS9ufhPoWbkVGkLHbu3ElqaioWi4WMjAx++OEHpk6dmuOc33//vUhlAPTs2ZNvv/2WmjVrAlCqVCliY2MBGDt2LGlpaddd44hyfZmpU/XP/9xz7pbERrzIh2vahfOIji6iOwF3LbETEAAWS5Gy2LVrF5cvX6Z+/fp07tyZhx56iJYtW9qV1wMPPEDHjh2JjIwkNDSUr776CoDMzEyioqJydBwMGzaM9evX07hxY1q2bEnHjh0ZM2ZMke6lOJGcDDNmwL33Fr3yGq7HtAvn4QiF69UW7qxZsxwyRnD27Nm57j9w4ACDBw8mMDDwyr7Q0FC2b99e5DKLK198AUlJ8MIL7pakkHiJS8G0C+dw4QKcPl10havEORUp/0yDgmDYMPj4Y7sLqFWrFkeOHCEgwD3/GYbCk5amxzA2agTh4QWe7o73+Nzr7QMPwI4dEBnpYnEKj2kXzmHXLmjdGubO1W9n+ZBvvfVaC/f48eMOEsbgKmbPhhMn4Msv3S2JHXiJhWvahXNwxBhc8GIfrsG7ENFjbps3h7593S1NIfGiTjODc8gag1uUWWbgxRauwbtYtgz279fDwYz+Mngb0dHaE1qhQtHyMRauwemIwFtv6aAf993nbmnsxEtcCgbn4IgRCuAuhWss3GLFmjXw5586eHOJEu6Wxg6MSV7s8W6FayzcYsVbb0GNGjBypLslKQLGwi22pKfDsWOOUbjGh2twKr//DuvWwZQpUKqUu6WxE2PhFmuOHtX2ofdauP7+XmvhrlixgsaNG9OgQQMmTpx43fGIiAg6duxIqVKlmDx5shsk9CzeekvHSxg92t2SFBFj4RaIr7aNrOHXjRoVPa8CLVylVGfgDaCO9XwFiIjYP0AiIMArLVyLxcJTTz3FqlWrCA0NpV27dgwcOJCmTZteOady5cpMnTrV7qDPvsSWLfDrrzBxIpQp425pioCxcAvEl9tGlsJt3Ljoedli4X4FfAB0AdoBba1b+/FSC3fLli00aNCAsLAwSpYsyf3338/ixYtznFO1alXatWtHCa/sHXIsb78NlSrBk0+6WxKDs/HlthEZqYeEBQUVPS9bfLhJIvJL0YvKXqoDOs3GjNHz7RxJq1ba2ZgHJ06coFatWle+h4aG8ueffzpWBh9h925YsgTGj4cbbnC3NA7AW1wKbmgX4Ntt4+BBx7gTwDaFu1YpNQlYAKRm7RSRHXaX6qWdZrnFnVDmdTNX3npLK9pnnnG3JA7A/MYF4sttIzLScbMjbVG4Hazbttn2CdDD/lIdYOEW8I/rDEJDQ3PMVY+NjaVGjRoul8PT2bUL5s2DV1/VLgWfwFssXDe0C/DdtpGcDCdPOsZ/CzYoXBHp7piisuHvD6mpBZ/nYbRr145Dhw5x5MgRatasyZw5c/jhhx/cLZbH8frrULEiPP+8uyVxED5iqTkTX20bBw/qrcsUrlKqAvA60NW6ax3wpogk2V9qAFy6ZPfl7iIgIIBPPvmEvn37YrFYGDVqFDfddBMzZswA4PHHHyc+Pp62bduSnJyMn58fU6ZM4cCBA5QvX97N0ruGLVu073bCBK10DcUDX20bWQrXUT7cAuPhKqXmA/uAmdZdI4CWInJPPpfln2n//jqarw8toWzQ9OsH27bBkSNF7izznHi4Dz2kZ2/ExLhWGoPbef113R9x6ZLNE3eKHA+3vohkX9h+vFKqaN2gZmqvT7JxI6xcCe+95yMjE7IwLoViS2SkDrrkqFmStozDvayU6pL1xToR4nKRSvXSUQqG/Bk3DkJC4Kmn3C2JE/CWTjODQ4mMdJz/FmyzcJ8AZlp9uQpIBEYWrVRj4foa4eGwdq3uJC9b1t3SOBhj4RZLMjO1D7dbN8flacsohV1AS6VUeev35CKXaixcn0JEW7c1a8I//+luaZyEsXCLHSdOaN+tozrMIB+Fq5QaLiLfKaWeu2Y/ACLygf2lGgvXl/j5Z/jjD5g+HUqXdrc0TsBYuMWS/fv19qabHJdnfhZu1othbt0fRfu7Nxauz5CRoQOLN2oEjzzibmkMBsfhDIWbZ6eZiHxm/bhaRMZnT8CaIpXqpcFroOAQdIsXL6ZFixa0atWKtm3bsmHDBjdI6TpmzoQDB+Cdd7x0NQdbMS6FAvG1trF/v+4EdkTQmiuISL4J2GHLvmtS/oweLVKtWoGneRoZGRkSFhYm0dHRkpqaKi1atJD9+/fnOOfChQuSmZkpIiK7d++Wxo0bu0NUl5CSIlKjhsgtt4hYb9mRFFg3nZByZ9QokdBQB9+eb+GLbaN9e5EePQp9Wb51LD8fbkegExB8jR+3POBfJC1fooRXuhSyh6ADroSgyx7zs1y5clc+p6Sk+EwAj9z46COIi4M5c4qBm9NYuPnia21DRL+5PfywY/PNz4dbEihnPSe7HzcZuLdIpZYsCWlpRcpizJgx7HJwGLpWrVoxxQHhGRcuXMjYsWNJSEhg2bJlDpXRUzhzRgcWv/NOuPVWd0vjZDxYMVyLO9oF+F7bOHYMLl50rP8W8vfhrhPtr71FcvpwPxCRQ0Uq1QEK1x2IjSHo7r77biIiIli0aBHjxo1zhWgu5+23dYXMxVVnKIb4WttwRocZ5D8sbIqIjAE+UUpd9zRFZKDdpTpA4Rb0j+sMChuCrmvXrkRHR3PmzBmqVKniChFdwpEjMG2aft3K9sbo23iJS8Ed7QJ8r224XOECs6xbx6/2VrKknsZhsegRC16CLSHooqKiqF+/PkopduzYQVpaGkEO7eZ0Py++qIdSv/GGuyVxEV7kUnAXvtY29u+H6tUdH885T4UrItut23VZ+5RSlYBaIrKnSKWWLKm3aWkQGFikrFyJLSHo5s+fz7fffkuJEiUIDAzkxx9/9OjOgcKyfj3MnauVbWiou6VxIV5i4boLX2sb+/c73roF28Iz/gYMRCvnXcBpYJ2IPJfPZfln+uGH8NxzcP48VKhQKIEN7sNigXbtdIdZRITTV+L1nPCMjz2mp9PFxblYHIM7yMzU0e5Gj7ZrAY18660t0cIqiI6fcA/wPxFpA/QqtBjZyW7hGryGb76BnTt1+EWvXvbcHoyFW2yIidExFJxh4dqicAOUUtWBocDPDinVKFyvIzkZXn4ZOnWC++93tzQuxkNfew3OIWtUXatWjs/blvCMbwIrgQ0islUpFQYUfVgYGIXrRbz9NiQk6Ddro38Mvszu3eDnB82aOT5vW8IzzgXmZvt+GBic9xU2YBSuVxEdrX1ZI0dqH26xxLgUig27dumg487oz7dlEclgYDRQN/v5IjLK7lKzopwYhesVPP+8/o985x13S+ImjElfrNi1Czp3dk7etrgUFgO/A6sBx4T4Mhau1/Dzz7B4se4oq17d3dK4EWPhFgvOndPTep980jn529JpVkZEXhSRn0RkflYqUqlerHALCkGXxdatW/H392fevHkulM6xXLoEzzyjZ5ONGeNuadyIsXBtwhfaxu7deuuMDjOwTeH+rJTq79BSvVThWiwWnnrqKX755RcOHDjA7NmzOXDgQK7nvfjii/Tt29cNUjqOd9/VQ2Q+/fTqT2Yw5IavtA1njlAA2xTus2ile1kplayUuqCUKtq6ZlmtNz29SNm4muwh6EqWLHklBN21fPzxxwwePJiqVau6QUrHcPAgvP8+DB8Ot93mbmk8AONSyBdfaRu7d+ug4yEhzsnfllEKuS2xUzQcYOGOWTGGXfEODkNXrRVT+hUtPOOJEydYuHAh4eHhbN261aHyuQoRvdR5YCBMmuRuaTwAL3IpuKNdgO+0jV27nGfdgm2dZlkxFBoCV5YIFJH1dpfqpS4FW0LQjRkzhvfeew9/LwrKcy1z58Lq1fDJJ1Ctmrul8RCMhZsvvtA20tJ0DAVnejtsGRb2KNqtEIqOpXALsAnoYXepDlC4Bf3jOgNbQtBt27aN+61Tsc6cOcPy5csJCAhg0KBBLpXVXpKS4N//hptvhscfd7c0HoIXWbjuaBfgG21j717t5Wzd2omFFLQGD7AXbdnusn5vAvxYwHX5c+iQCIjMmlWYtYLcTnp6utSrV08OHz58Zd2mffv25Xn+Qw89JHPnznWhhEXn8cdF/PxEtmxxtyQetKbZk0+KBAU5+PZ8C19oG9Ona7UUHV2kbOxb0ywbf4vI30oplFKlRCRCKdW4SFreS10KtoSg82bWr4cZM3Qgt2I7oyw3lDIuhQLwhbaxdateobdePeeVYUt4xoXAw8AYtBvhHFBCRPIbKpZ/pvHxehT99OnmvdVD+PtvaNlSv1Lt3Qtly7pbIg8Kz/ivf8F330FioovFMbiSli21WlqxokjZ5FtvbRmlcLf14xtKqbVABaBoImVZuKmpRcrG4DjefFMPBVu1yiOUrWehlA6SavBZLl3SHWYD7V84zCbyVbhKKT9gj4g0g5yrPxSJrKgQly87JDtD0di1S4+5HTkSehUt0rFv4udnXAo+zs6dVwPsO5N8Jz6ISCawWylV26GllraOLjMK1+1kZMAjj0CVKvB//+duaTwUY+H6PFnDgp2tcG3pNKsO7FdKbQFSsnZKUVbtVUovGXDpkt1ZGBzD5MmwYwf89BNUruxuaTwUY+H6PFu3Qs2azg/QZIvCHe+UkgMDjcJ1M3v2wGuvweDBcO+97pbGgzEWrs+zdSu0bev8cmyJpdBfRNZlT0DRg9kYC9etpKbCiBHaqp0xw6vG9rseY+H6NImJcOiQa4ZC2qJwe+ey7/Yil1ymjFf6cAsKQTdp0iRatWpFq1ataNasGf7+/iR64HCi8eO1hfvFF9p/a8gHY+HahLe2jU2b9LZTJxcUlteMCOAJ9CyzFGBPtnQE+K6AGRUF07q1yJ13FmlKh6vJyMiQsLAwiY6OvjKbZv/+/Xmev2TJEunevbsLJbSNP/7Qs8lGjXK3JPniOTPNxo4VKVHCwbfnW3hz23j5ZRF/f5GLFx2SXb51LD8L9wfgTmCJdZuV2ojI8CJrei90Kdgagi6L2bNn88ADD7hQwoJJSYEHH4RateDDD90tjZdgLNwC8ea2sXGjjp9gy/jz2XtnM/638aRZ7Jslm2enmYgkAUmAc55KYKBu/XYyZszVYMGOolUrvVhiXtgSgi6LS5cusWLFCj755BPHCllEXnwRoqIgPBzKl3e3NF6CF/lw3dEuwHvbRno6bNkCjz1W8LlpljTGrhlLcNlgXrvtNbvKsyk8o1MoUwbOnHFb8fYgNoSgy2Lp0qV07tyZyh401mrpUpg2TTfK7t3dLY0XYSzcAvHWtrFzp+5KssV/+9WOrziadJTP7vgsz3sriDwVrjVQjfPm3hbRpVDQP64zsCUEXRZz5szxmFcmgLg4ePhhPV/83XfdLY2X4Wf1vIl4/HAOd7QL8N628ccfelvQKr2X0y/z1u9v0blWZ/rU72N3efn5cDcBKKVm2Z17fnihD7ddu3YcOnSII0eOkJaWxpw5cxiYy+TrpKQk1q1bx1133eUGKa/HYtFL5Vy+DHPmXJ3oZ7CRLCXrJW4Fd+CtbWPjRqhTR096yI/Ptn9G3IU43urxlt3WLeTvUiiplHoI6KSUuufagyKywO5SQftwvWxYmK0h6BYuXEifPn0o6yFRYN57D9auha++giZN3C2NF5LdwjXkije2DRGtcAtyr6WkpfDuhnfpWa8n3ep2K1KZeYZnVEp1Af4BDEWPVMgpq8iofPItuGa++CJ89JGOC2hwGn/8AV27wpAh8MMP7n8jTs1IZVPsJlZFryI8JpyVw1dSvlSuvXeeE57xrbdg3Dgdv7lECReLZHAWUVHQsKFelfqJJ/I+b+KGiYxdM5Y/Rv1Bx1odC8rWvvCMIrIB2KCU2iYiXxVUSqGpUEFPd0pNhVKlHJ69Ac6fh2HDoHZt980mExH2Jexj1eFVrDq8ivVH13Mp/RL+yp8OoR04dfFUXgrXczAWrk+ydq3e5mfhnr10lokbJjKg4QBblG2B2DJKYZZS6l9AV+v3dcAMESnaGucVKuhtUhJ46JLJ3oyIDrd44gRs2HD1cbuCE8knWH14NasOr2L14dWcSjkFQOOgxoxqNYre9XvTrW43z1e0WWT9U5mRCj5FeLgOVtM4n/VrJqyfwIW0C7zX6z2HlGmLwv0UKGHdAowApgOPFqlko3Cdyvvvw+LFute6QwfnlnUp/RLrYtaxMnolqw6v4sDpAwAElwmmV1gveof1pldYL2pVqFVATh6KsXB9DhFt4fbsmfebX3RiNJ9u/ZRRrUZxU9WbHFKuLQq3nYi0zPY9XCm1u8glV6yot+fPFzkrQ07WroWXX4ahQ/XqMI4my02wMnolK6NX8vvR30m1pFI6oDRd63Tl4VYP0zusN81DmuOnbAnX4eEYC9fniIiAU6fydyeMXTOWEv4leLP7mw4r1xaFa1FK1ReRaAClVBhgKXLJ2S1cg8M4cQLuvx8aNYIvv3Sc3/bspbOsOryKldEr+TX6V+IuxAFwU/BNPNXuKfo26MuttW8lsESgYwr0JIyF63MU5L/dHLuZuQfm8lrX16h+g+OC5NqicF8A1iqlDqN74OqgF5UsGlkK11i4DiM9XVu1KSnw229www3255WRmcHm2M2sjNJW7La4bQhCpdKV6BXWi771+9K3QV9Cy4c6TH6PxVi4Pkd4uI4nEhZ2/TER4YVVLxBSNoQXOr/g0HJtWURyjVKqIdAYrXAjHDIDLculYCxch/HCC3oY2Jw5cOONhb8+5nzMFQW75sgaklOT8VN+dKjZgddve52+DfrSrkY7/P38HS+8J2MsXJ8iM1MbJHfckfsb4E/7f2LDsQ3MGDCDciXLObRsm2IpWBXsHoeWbFwKDuWbb/Sw5mefhfvus+2a1IxU1h9dz/JDy/kl6hciz0YCUKt8LYY2HUrfBn3pWa8nlQIrOU9wb8BYuD7Fjh1w9mzuC6ampKXwn1X/oVW1Vjx6c9HGBeSG+4LX3HCDrsjGpVBkNm6Ef/5TV6DJk/M/NzY5ll8O/cKyQ8tYfXg1KekplPIvRbe63Xi87eP0rd+XJlWaFGn6os9hLFyf4pdftOrp2/f6Y+/8/g6xybHMGTzHKW9y7lO4fn46PqCxcIvEsWNwzz16PvhPP0HANb9oRmYGf8b+ybJDy1h+aDm7T+kBJrUr1ObBlg/Sv2F/etTrQZkSZdwgvZdgLFyfYsUKaNMGgoNz7o9KjGLypsmMaDGCzrULiGZjJwUqXKXUfOBr4BfRy6Y7jooV4dw5h2ZZnEhJgYED9WS9JUugkvXN/8ylM6yIWsGyQ8tYGbWSc3+fw1/506V2F97r9R4DGg6gaXBTY8XairFwfYbERNi8GV555fpjY1aMoZR/KYdNcsgNWyzc6ehRCVOVUnOBb0QkwiGlV62qB8MZCk1mJjz0EOzdCz//LFyusIsJ635medRy/oz9E0GoWrYqAxsPZEDDAfSu35uKpSu6W2zvxFi4PsPq1fpn7Ncv5/6lkUtZdmgZk3tPzncYWFRUFNHR0fTNzR9hA7aMUlgNrFZKVUCv/rBKKXUc+AK9tpnNU3x//PFHKlWqRJ8+1niSISGQLYamwXbGvmxh/nx/Oj0yn8cixhC7JRaAdjXa8fptr9O/YX/a1GjjGxMP3I2xcH2GX37Rb4Lt21/ddyH1Ak8tf4qmwU35V4e8ZwodPXqUnj17kp6eTnR0NIGBhR9zbpMPVykVBAxHT+vdCXwPdAEeArrZWtibb75J06ZNcyrcbdsKJ3Ex5tzlcyw/tJwPp11i++ejoc1n7Kz3HH1r9OHNbm/Sv2F/QsqFuFtM38NYuD5BZqb23/bunbOv49XwV4lNjmXjqI2U8M89GtzJkyfp1asXSUlJhIeH26VswTYf7gKgCTALuFNETloP/aiUKpS2DAwM5HL2GLghIXD6tH4SfsYSy40j546wJHIJiyMXs/7oeiwR/WDOYmq13cPHM0Pp0/CMR83uysjI4MyZMyQkJHD69Ok8U9bxHTt2ULt2bXeLnT/GwvUJtm6F+Hg9/jaLP2P/5OMtH/NkuyfzjAZ2+vRpevXqxcmTJ1m1ahU333yz3TLYYuF+KSLLs+/IWn5HRNoWprAyZcpwKfsqDyEhejmCs2ev7zIspogIW+O2siRyCUsil7A3YS+gp9AOD/qIHxc9TrM2fqxd24Jy5Vo4XZ60tLR8Fea16VwenaB+fn4EBQURHBxMcHAwzZo1Izg4mBLeEF/WWLg+wcKF2rLNUrhpljQeXfooNcvX5J2e7+R6zfnz5+nbty+HDx9m+fLldOxYtBCNtijct4Dl1+zbBBRazQcGBpKUfRhYiPX199SpYq1wLZkWNh7fyPwD81kQsYDY5Fj8lT+31rmVD/p8wMDGA1Hn69OxI1SvBj//DOXsnAAjIiQlJXHq1KkcKT4+nlOnTl1RpFnb5OTkXPPx9/enSpUqVxRo69atr3zOLVWuXBl/fy+doWYsXK9HRCvcbt2ujuaZtHES+xL2seT+JbmGCk1KSuL2229n3759LFmyhO4OWHk1v0UkqwE1gUClVGuuRjIvD9g1aLNMmTLEx8df3VGtmt7Gx0OzZvZk6bWkW9L5LeY35v81n0URiziVcopS/qXo26Avb/d4mwENBxBUJgjQ/0e39oOMDO30D7nGTSsiJCYmXqdE80ppaWnXyePn50fVqlWvpPbt21+nNKtWrXrlc8WKFfErLm4gY+F6PX/9BQcP6pmYAPsS9vHm+jcZ0nQIdza+87rzsyzbHTt28NNPP9Hv2mENdpKfhdsXGAmEAh9k23VI5v0AACAASURBVH8BeNmewgIDA3O6FLJ8d0eP2pOd15Gakcqqw6uY/9d8lkQuIfFyImVLlKV/w/4MvnEw/Rv254ZSVyPOpKSkEBFxivvvr8rx46X55z/n87//7bxOgSYkJJCefv1gkYCAAKpWrUpISAghISHcdNNNVz5npWrVqhESEkJQUFDxUaCFxVi4Xs/ChXp7113alTB8wXAqlq7ItP7Trjs3MTGRPn36sGfPHubNm+fQBS/zW2JnJjBTKTVYROY7orDrfLi1aoG/Pxw54ojsPZLL6Zf5JeoX5h2Yx88Hf+ZC2gUqlKrAbdVvo13ZdtROq82Z42fYtmUbS+KWcPLkSeLi4jh58iTJyRbgV/SLxu1MnbqaEiVK5FCWLVu2zFWBhoSEUKlSJaNEHYGxcL2ehQt1IP6aNeGVNePZfWo3i+9fTHDZnK7Ms2fP0rt3b/bv38+CBQu4I3sPmwPIz6UwXES+A+oqpZ679riIfJDLZflyww03kJSUhIjoWU4BAdrKPXy4sFl5JKmpqcTFxRFzLIYVh1aw5tQa9qbvJU2lEZAWQJljZSi9uzRJB5JYYlnCkmxrc5YuXZoaNWpQvXp1WrRoQc+e/fn112eIjq7DuHF7GTp0CtWrV6dSpUpmhpirMRauVxMTA9u3w8SJsOn4JiZunMioVqMY2DjnMu5ZoxEiIyNZtGgRt99+u8Nlyc+lkLWOscPik9WuXZuUlBQSExMJCtL+SerV8woLNzU1lRMnTnD8+HFiY2OJjY3N+fnEcRJKJ0AzoCnay/03+B/0p3pCdcJUGDWr16R6j+rUGK4Va5aCrVGjBhUqVLiiSDMydFzbqCgdBeyhh1rmI5nB6RgL16v54Qe9vfPuS9y16EFqla/Fh/0+zHFOXFwcffr0ITo6miVLllydK+Bg8nMpfGbdjndUYfXq1QPgyJEjORXu0qWOKsIuLBaLtkxjYjh69GiuSvX06dPXXVehYgWCWgSR2TaTi/0vgj+UVCXpVLkTg5sMZmjroQRXDi6URZqRAQ8+qF+Bpk7V03cNbsZYuF6LCHz/PXTuDB9E/IvoxGjCHwrPMSohKiqK3r17c+bMGZYtW0aPHj2cJo8tEx/eRw8NuwysAFoCY6zuhkLRoEEDAPbv30/bttYhvGFhkJAAFy/aP9apADIyMjhx4sQVhRoTE3MlHT16lGPHjpGRkZHjmsqVKxMaGkpoaCjt27e/8jk0NBRVSRF+NpwfI34k+lw0Jf1LcnuD23mg2QPc0egOypYsm4ckBckJI0boAOITJ8Izzzji7g1Fxli4XsuePXDgAIx8ZQtf7fyKV259hW51u105vmvXLvr27YvFYiE8PJx27do5VR5bxuH2EZH/KqXuBmKBIcBaoNAKt2nTplStWpWlS5fyUJbplrU0wYEDOSc4F4L09HRiY2PzVKjHjx/HYsm5DFuNGjWoU6cOHTp04L777qNu3brUrVuXOnXqUKtWLcqUyTnyLTk1mbn75/Lm7jf5fePvKBTd63XnlVtf4e4b7y5yYJiMDBg+HH78Ua+4+4JjV/YwFAVj4Xot338PAQHCTzKELrW78Ea3N64cW79+PXfeeSfly5dn1apVNGnSxOny2KJws6YC9Qdmi0iivZ02fn5+jBw5kkmTJrFmzRp69uwJLa3+yd2781S4aWlpHD9+PIcSzf45NjaWzGzWh1KKGjVqULduXTp37nxFkWYp1Vq1alG6dOkC5bVkWlh9eDUzd89kYcRC/s74m0ZBjXi7x9sMbzGc2hUcMyU1PR3+8Q+YOxcmTYL//Mch2RochbFwvZLMTPhhtlDmxvWUqJDC7MGzCfDTKm/p0qUMHTqUOnXq8Ouvv7pserktCnepUioC7VJ4UikVDPxtb4Fjx45l2bJl9OnTh4EDB9KhfXuqli5N5uzZpKWnc/bs2SvjS0+ePElMTAxxcXFINuvCz8+PmjVrUqdOHbp27XpFkWa3UEuVKmWviEQnRvPlji+ZuXsmJy+epFLpSoxqNYoHWz5I+5rtHTpKIC1NK9t58/RqDc8/77CsDY7CWLheydq1cCJWQYfp/Dxo5pUFT6dPn87TTz/NzTffzPLlywl24SxXJTZUIqVUJSBZRCxKqTJAeRGJz+eSfDM9d+4cb7/9Nj/99BPHcwnPWKlSpStjSrMs0+wWamhoqMPn4KdZ0lgUsYgvdnzB6sOr8Vf+9G/Yn5GtRjKg4QBKBdivwPPi0iW49149e+z//g+eu27wXbHGHWPfcq+3CxfqZTV27oRWrVwsksFe2vc9zNb1Ffn3vMl8MOAdMjMz+e9//8v//d//MWDAAObMmUM5x/cb5V9vRaTABHQChgEPZqUCrrGZixcvSszIkXKsdGk5efSopKamFubyInPwzEF54dcXJPj9YOENpM6HdWTCugkSmxTr1HLPnRPp0kVEKZHPP3dqUd6KTXXTwSl3liwRAZGtWx18iwZnsXzXVsE/VWr3nS/plnRJSUmRe+65RwB56qmnJD093VlF51vHbBmlMAuoD+wCsnqeBPjW7v+AbJQtW5ayAwfqAaexsVen+zoRS6aFpQeX8vGWjwk/Eo6/8mdg44E81uYxeof1dvoy4KdO6QXsDhzQnWRDhji1OENRyXqbumYki8EzOXnhJPe/uhwsbZn9Tg/Onj7LwIED2bp1Kx9++CHPPvus2yYP2eLDbQs0FXGiA+vWW/X2t9+gUyenFXPu8jm+2vkV07ZOI+Z8DLXK1+Kt7m8xqvWofJfVcCQxMToAclycHn5s50odBleSFa06l3gVBs8iNSOVwT/dy4VN39C6Qwpl/GLo0GEQCQkJLFiwgEGDBrlVPlsU7j6gGnCyoBPtpkoVaN4cwsPhZbvi4uTLgdMH+PjPj/l2z7dcSr9E1zpdmdx7Mnc1uetKr6Ur2L5dx+L8+29Ytcqp/y0GR2IsXK9ARBi9dDSb1peCsw25ZehmOnXqQeXKlVm/fv3Vsf9uxBZtUwU4oJTaAqRm7RSRgXlfYgf9++ueo7NnIWsWWhHZcGwDEzdMZNmhZZTyL8Ww5sN4pv0ztK7e2iH5F4aff4b77tP/LatWFbtolN6NsXC9gtd/e51Ze2bROCqCY4EXmT69O126tGXevHmEXBvT1E3YonDfcLYQgHZkvvceLFoEjzxidzYiwvJDy3l3w7tsPL6RoMAgxncbzxNtn7guMpCrmDYN/vUvaN1auxGqu8Z7YXAUxsL1eL7e+TUT1k9gYMX/sOSPxsB4Hn98JB999BElS5Z0t3hXKahXzeq6rQP0sn4uA9xQwDWFJzNTJCxMpHt3uy63ZFpkzt450vzT5sIbSO0Pa8vUzVPlYupFu/JzBBaLyHPP6Q7uO+8Uueg+UbwRzxmlsH27/hEXLXLwLRocwYpDK8R/vL90nNZRypefJXBZJk361l3i5K9LCzwBRgNbgWjr94bAmgKus49339Ui7d1r8yWZmZmy6K9FVxTtjZ/cKDN3zZS0jDS7xXAE58+L3HGHvp1nnhHJyHCrON6I5yjcPXv0DzlvnoNv0VBUNh/fLOXeKSc1J9SUgLIhotQlueOOeHeKlG8dsyU69VNAZyDZahEfAqo62NDWjB4NgYEwZUqBp4oIK6NW0v7L9gz6cRB/Z/zND/f8wL4n9/FgywfzXO7YFURE6GDHK1bAJ5/oqF/eupyXAePD9VB2x++m73d9kYvCifdPUK/GZEQCefddz/DX5oYtCjdVRK4sgqWUCqCAmWR2ExSk/bfffAORkXmetj9hP32/60u/7/uRkJLAVwO/4sBTB3ig+QP4KfeucLB0qQ4JkZgIq1fDU0+5VRyDIzA+XI8j8kwk3f/XnYuJF7k84zJvPD+FhIR/cPfdnt0hbYt2WqeUehm9mGRvYC7gvAC248ZBmTLw739fN3f97KWzPL38aVrOaMnWuK182PdDDj59kFGtR7l0eFduZGbCW2/pNZMaNoRt2+C229wqksFRGAvXo4g6E0WHTztw7tw5qiyvwvol6xF5lqQkxWuvuVu6AijI54BWyqPRinae9bMq4Lqi8eGH2mf2xRciov20X27/UipNrCR+4/3kyZ+flNMpp4tcjKM4fVrk9tu1yMOHi1y65G6JfALP8eHGxuof97PPHHyLhsKyds9aKfVSKeFFpMcDPeT06dOSmChSvrzI3Xe7WzoRKWqnmWilGwwE23KuOELhWiwiPXuKBAbKobXzpfs33YU3kNv+d5vsPWV7h5orWL9epGZNkZIlRaZN04MtDA7BcxRufLxuKtOmOfgWDbaSmZkp73zxjqj/KFH/VTL+s/GSaW1sr7yif55du9wspMY+hYuOevMGcAY4CyQCp4HXCsrUEVJnxsXJlNsrS+lXkfJvl5PPt30ulkyLI7J2CBaLyDvviPj7izRoILJjh7sl8jk8R+GePaubypQpDr5Fgy2cOnVKeg3vJbyAlBhbQlbsWHHl2LFjIqVLi9x/vxsFzEm+dSw/H+4Y9OiEdiISJCKVgQ5AZ6XUvx3v3LjK6ZTT3Ll2NGM6JNIztiQHPgtgdFozt3eIZXHypJ4Y9/LLOrzi9u16UoPBR8nqNEtLy/88g0MREWbNmkXDng1ZXXM1FcpWYPeY3fRtfTUAycsv666eiRPdKGhhyEsTAzuBKrnsDwZ2FqDJ7WbjsY1SfXJ1KTmhpEzdPFUyDx0SqV9fpEQJkUmTtGnpRubOFQkK0v+q06cbF4IT8RwLNyNDW7hvvOHgWzTkRUxMjPTr10+4EVHjlNSZXEcOJx7Occ6ff+qfZexYNwmZO3a7FPbZc0yKoHBn7Z4lJSeUlPof1ZedJ3dePXD2rPaIg0iHDiKbNtlbhN2cO6c7xECkbVuRv/5yuQjFDc9RuCIipUqJ/Pe/Drw9Q25kZGTI1KlTpWzZslLy1pKi3lByy5e3XNdJnpEh0r69SNWqIsnJLhQwKUkkKiq/M+xWuDvsOSZ2KNzMzEx587c3hTeQbt90kzMpZ3I7SWTmTJHq1bXYAwaIhIe7xMRcs0YkNFT7a19/XSTNvZPYiguepXArVxZ56ikH3p7hWjZt2iRt2rQR/JA6/6wjvIEMmjNILqVdP+xn6lStBr7/3sFCZGaKJCSI/PGHyKxZusEPHy5yyy0iwcG60Hr18svBboVrQc8uuzZdANILyLgQ95cpr6x5RXgDeXDhg5KaUcCKDxcuiEyYcPXmmzYVefttkcOH87/ODs6dExk9WhfTqJF+hTG4DM9SuLVqiYwc6cDbM2Rx6tQpGTVqlABSrV41af6+nqb/zPJnJMNy/Zz4Y8dEypUT6dfPTnvLYhE5flxk7Vo99PSll0TuvVekdWuRG27QDT4rKaV/++7dRR59VGTiRJEFC/LLPd86ZtOaZva4hm098a31bzFu7ThG3zyaGXfMsL1j7PJlvQbyzJmwYYPe16wZ9Omjo3p36gR2rlckAvPnwzPPQEKCXmts/Hg9H8PgMjxnTTOAJk2gRQv46ScXiuPbZGRkMGPGDMaNG8fFixcZ8fwIfqv2GycunuDT/p/yyM3XRw0UgQEDYN062LcP6tXLI/P0dDh6FKKjISoq5zY6GlJTr54bEKAzql8fGjTIua1XD2xY4Tsb+dZbtyrc2XtnM2zBMEa0GME3g76xfxRCTIxeY3zlSvj9d92b7OcHTZvqebbt22tlfOONULlyvlmdOKGn4y5erEcefPEFtGljn1iGIuFZCrdNG6hWDZYtc6E4vomIsHz5cl588UX2799Pr1696PdCP17f9jrlS5Vn/tD5dKzVMddrP/lEG0JTp8Izj16Gw4evV6hRUVrZWixXLwwMzF2hNmgAtWpdnU1YdDxT4e5L2Efbz9vSIbQDvw7/1XGr4qakaKW7eTNs2aLT2bNXj1etqhVv48ZQp45eQ61OHdKq1ebjRaGMn+BPejq8+aaeXey438FQSDxL4Xbtqv/Ef/vNddL4INu2beOFF17gt99+o2HDhox/dzyrS6zm651f0zG0I/OGzqPGDTWuXpCcrBVoVBT7NibRdtpIelbYxs+lh6DiTuTMvGLF3BVq/fo6CLVr1jHzPIWbZkmj/RftibsQx94n9hJSzonRfUT0v92BA/DXX1fToUNw+jQAv9KbZ/mICG7k9lJr+DhsCvVrpWnlHBKit1WrQnCw/lGzpzJlXPVDFjc8S+HefjucOQNbt7pQHN/hyJEjvPLKK8yePZvg4GBef/11Og7qyPBFw4k4E8HYhqN4I6AnJaKO6LZ56BAcPHiljaZQhlvYzGm/EPa0GUXVplWuV6wFvL26CM9TuFM2T+HfK//NwvsWMqiJ+xZ1O3LgMs89m8Gi1TfQoGoSU7ovYUC5dXpZ3YQEnU6d0v7ivAgIyKmAK1TQvuOyZXNPZcpcv690aZ1KlcqZSpcuzia2Zynce+/Vf9oHDrhQHO/n6NGjvPvuu3z99dcE+Pvz3IABPN/8Jr4+v4JXbthKpb8Vs+Zn0utwtotq1NARoKxJ6jfgvs97MH91BX75RdGnj9tuxxY8S+EmpyZTd0pd2tZoy68jfnVG2QVy7pyemfLRR1qfvfqqdh+Uys2rIaLdFAkJ+t82KQnOn78+Zd9/8aK+JnvKzLRPWD+/nAo4N6VcsqS+kRIldMr6nNe2MOf4++vk53f1syP3BQXlFSzYsxTu44/DggW6HhhyJzVV+1QPHeLYli28s2ABX0dEoER4FHgZuFAFHr4LNteCO09X4ovLfQip3+Kqgm3Q4LrO7nfegVdegUmT4D//ccudFYZ8663LzadZu2dx7u9zTOg+wdVFc/mydrq/847Wj8OH68+hoflcpJSuAOXKQViYfQWL6Mp4rRLOSqmpeinf1NSr6drvBe1LSdE9sxkZuW9z2+cJHDyoG5qnU62adimkp1+d6lscSU+HI9le+7OnY8c4kpnJ+8BXaM0zOiSElzp3pnrL5nxY5RDjzs6jbMmyfN//Ex5o9gCqAHfcTz9pg+gf/4Dnn3fFDToXlyvcz3d8TpvqbegQ2sFlZVos8O238NprEBur4yC8+64e5eMSlLrqNnDQisRFRkRb3QUpZYtFn2ex5ExF2Zd9f7B7FvYsNNWq6Wd2+rR+5fVlRLQlHxmply+JiNCfDx7UyjZ773+FCtCwIdsbN2ZS2bLMPXAAf39/Hh0xgrHjx1OrVi3WH13PgOVPszdhL4OaDGL6gOlUK1etQDF+/VUbRZ0769FCvtBV4lKFeyzpGHtO7eGDPh+4pLyMDPjxRx0YPCJCjw6bNQu6dXNJ8Z6NUldf7w0FU82qIOLjfUfhpqfrYVTZlWrW5/Pnr55XurQe1dO6NQwdCo0aad9qgwas2LaNSZMns3blSsqXL8/z//kPzz77LDVr1iTuQhzDFwzn+73fU6dCHRYMXcCgJoMKtGoB/vgD7r5bDyhaulSP6vIFXKpw1xxeA0Cf+s71emdk6DkRb7+t33SaN4d58+Cee3zjX9LgBurU0dvoaLj5ZvfKUljOnr3eWo2I0PeS3VqtXl1P8HjgAa1gmzTRqVYt7XO3cuHCBb777jumPfYY+/fvp0aNGrz//vs89thjVKhQgQupF3hz3ZtM+mMS6ZZ0xnUdx0tdXqJMCdtmDoWHw8CB+n9t5UrdF+0ruFTh7kvYR2BAIDcG3+iU/NPS4LvvtKI9fBhatdL9HHfdlaO+GAyFp0kTXYn274chQ9wtzfVkZOgJQLlZq2fOXD2vZEntM2/eXI+8yFKqjRtD+fL5FrFv3z6mT5/Ot99+y8WLF2ndujXffPMNDzzwACVLliTNksbHf37MhPUTOH3pNPfceA/v93qf+pXr23wbS5fqx9uwoXYpVCvY8+BVuFThRp2Lon7l+g6Pa3vuHMyYAR9/rGPVtmmjZ4rdeaexaA0OIjBQv0pv2eJeOZKScrdWDx3K2REaHKwV6d1357RW69YtlBvp0qVLLFy4kM8//5z169dTqlQp7rvvPp588knat2+PUorUjFS+2P4F7254lyPnj9C9bncm9ppI+5rtbS5HRI8aev55/QKxYoXndHc4Etcq3MQoGgc1dlh+hw/rFdW//lp30vfuDf/7nw6nYBStweH06QOffw6XLjk3sEZmJhw7lru1Gh9/9byAAD3ov0kTuOOOnNZqESYBiAgbNmxg5syZ/PTTT1y4cIF69erx/vvv8/DDD1OlShUALqZd5IvtXzB502TiLsTRtkZbpg+YTp/6fWzy02aRmgpPPqnb8T336PAodoZB8XhcOg53VfQqAksE0qV2F7szzszUPp4ZM2DhQv1nPWyYDjDjslEHBlfgWeNwQUdM6dYNpk3TGqKoXLige/4jI3Mq1YMH9XC/LCpVyqlMsz6HhTl0iNqhQ4eYPXs23377LdHR0ZQtW5YhQ4bw0EMP0bVrV/ysfrnjSceZsW0Gn23/jLOXz9K9bnfGdhlLr7BehVK0oOeRDBsGu3frUUSvv+717j/PmvhgL/Hx2nr98ktt2VaqBP/8pw5k4SudxoYceJ7CFdEKd+dO/c7bqVPBOaana99q1rCq7NuTJ6+e5+enI1Ndq1QbN9buASe9sh06dIi5c+cyd+5cdu3aBUD37t0ZOXIk99xzD+WspqaIEH4knGlbp7E4cjEiwp2N7+Slzi/lGWgmPywW+PRT+O9/4YYbtHV7xx0OvTV34b0K12LRjvMvvtDO9IwMuO02eOwx/epRuKhpBi/D8xQu6IHcXbtqJdqjhx5rGBSkX+8vXNA+1hMntEvg2DH9Ofssw8qVtRJt1CjntkEDl1RoEeHAgQMsXrw4h5K95ZZbGDJkCPfeey+1a9e+cv6xpGN8v+d7Zu6eSeTZSIICg3j05kd5vO3j1K1Y1y4Ztm+HJ57QYSn69dOGlA91jnmfwj1xQv/jffmlrrNVqsDDD8Ojj+r6aSgWeKbCBa1UP/hAD4E5cCCnQi1dWr9y1a6tU61aWplmKVY39ARdvnyZtWvXsmzZMpYtW8bRo0cBrWSHDh3K4MGDcyjZ5NRkFvy1gG93f8tvMb8hCF1qd+GR1o9wf7P7KR1g3x/DsWM6Ct///qeN9g8+0CPQfKy/xTsUbkYG/PKLtmaXLdN1uFcvbc3edZcezWIoVniuws2OxaJ7bDMy9LuxB0z7FRH++usvwsPDWblyJWvWrOHy5cuUKVOGXr16MWDAAPr3709otjntpy6eYnHkYhZGLGTN4TWkZ6bToHIDRrQYwfAWwwmrZOe0dvRLwaRJut8FtHX7xhu+Nb42G56tcI8eha++0hbtiRP61eLhh+GRR3QHrKHY4h0K1wMQEQ4fPkx4eDjh4eGsXbuWU6dOAVCvXj0GDBjAHXfcwW233UZpq9vCkmlh+8ntrIpexfKo5Ww6vglBqF+pPnc3uZvBTQfToWaHQneCXZUJ/vxTD/WaO1fve/hhGDdOG/4+jOcpXItF9zl8+qm2akH7ckaP1o5zDzASDO7HKNw8SE1NZceOHWzatInNmzezadMmYmNjAahevTo9evSgR48edO/enXrWNWhEhMizkayLWceqw6sIPxLOub/PAdC6WmsGNRnE3U3uplnVZnYrWdDW7Pff6wlI+/bpUAuPPgpPP62HABcDPEfhpqdra/a993SfQ7VqWsk+8sjVmZMGgxWjcNHrfkVGRrJr1y62b9/Opk2b2LFjB2lpaQDUqVOHjh070qVLF3r27Enjxo1RSnEx7SJbTmxh0/FN/BH7B5tjN5N4ORGA0PKh9A7rTe+w3vQM60nVslXtlk8E9u7VbsBly3QMBBHo2BFGjNDJV8fU5oFnKNzkZD3za/16/WP8+98waJCxZg15UuwU7tmzZ/nrr7/YtWsXu3fvZteuXezdu5dU64KHpUuXpm3btnTs2JGOHTtyyy23UK1aNY4nH2d3/G52n7Km+N1EJUYh1tu5scqNdKrViY6hHelSuwuNghrZbcVaLNpy3bhRp3XrtCsQ9Ayxu+7SoRSLsTvQMxRuaqqeZThsmP5BfKxn0uB4fFLhpqWlcfToUSIiIoiMjMyxPZMt5kFQUBCtWrWiVatWtGzZkrCbwqAyHEk+wqGzhziYeFBvzx7kQtqFK9eFVQqjZUhLWoa0pH3N9twSeguVAisVWk4RSEzUgzD27tVp3z49QeGCtbjq1aFLF71I9u23m/HwVjxD4YL+EY2iNdiI1ylcESE5OZn4+HiOHTtGTEwMR48eJSYm5srnEydOkL3NBdUMot5N9ajRoAZValchMCSQgEoBnJfzxCbHcjz5OLHJsVxKv3TlGj/lR50KdWgU1IiGlRtyY/CNtAxpSfOQ5pQvlX8AGtBWamJizpWk4uO1my8mRoe8jYnRb6VZVKig4920aKHne3TqpH2ypj1fh+coXIOhEHiEwr2iRE/Fczz+OMdPHScuIY64M3HEn4kn4VwCp8+f4WxyIudTkskQC5RQEOCntyX8KVuxHIE3lKV0uUBUoD+WEpn8TRrJ6RfIsFj0rYofiN4q5U/VwGpUv6EGIWVqULVUKFVK1aB6YF1CStcmqFR1JKPklUU/0tKuLgCSkqIt0OTk3LeJiTp4WG4rPpUpoye7ZaW6dfVkt+bN9aooRrnahFG4Bq/EIxTujFl7eeKRsCvK8FrlCJ438T8gQEdaLF9eDw3Ovq1YMedC1NlTUJBRqg7A9QpXKbUCqJLH4SrAmTyOFSfMc9Dk9RzOiEg/Vwpi6m2BmGegye855FtvnWXh5olSapuItHVpoR6IeQ4ab3kO3iKnMzHPQFOU5+B570MGg8HgoxiFazAYDC7CHQr3czeU6YmY56DxlufgLXI6E/MMNHY/B5f7cA0Gg6G4YlwKBoPB4CKconCVUv2UUpFKqSil1Eu5HFdKqanW43uUUjc7Qw53Y8Nz6KaUSlJK7bKm19whp7NRSn2tlEpQSu3L47hH1AdTbzWm3jqxzoqIQxPgD0QDYUBJYDfQ9Jpz+gO/oAcJ3wL86Wg53J1sfA7dgJ/dLasLdC73pwAABDBJREFUnkVX4GZgXx7H3V4fTL0t1HPw+XrrrDrrDAu3PRAlIodFJA2YA9x1zTl3Ad+KZjNQUSlV3QmyuBNbnkOxQETWA4n5nOIJ9cHUW42ptzivzjpD4dYEjmf7HmvdV9hzvB1b77GjUmq3UuoXpdRNrhHN4/CE+mDqrcbUW9uwqy4EOEGQ3OYSXzsUwpZzvB1b7nEHUEdELiql+gOLgIZOl8zz8IT6YOqtxtRb27CrLjjDwo0FamX7HgrE2XGOt1PgPYpIsohctH5eDpRQSuU1l9+X8YT6YOqtxtRb27CrLjhD4W4FGiql6imlSgL3A0uuOWcJ8KC1p+8WIElETjpBFndS4HNQSlVT1tD7Sqn26N/jrMsldT+eUB9MvdWYemsbdtUFh7sURCRDKfU0sBLd4/m1iOxXSj1uPT4DWI7u5YsCLgEPO1oOd2Pjc7gXeEIplQFcBu4XaxeoL6GUmo3u2a6ilIoFXgdKgOfUB1NvNabeapxVZ81MM4PBYHARZqaZwWAwuAijcA0Gg8FFGIVrMBgMLsIoXIPBYHARRuEaDAaDi3DGTDOfRCllAfaih4ZkADOBKSKSy4LTBoP7MXXW8zAK13Yui0grAKVUVeAHoAJ6fF6RUEr5i4ilqPkYDNdg6qyHYVwKdiAiCcBjwNPWmSb+SqlJSqmt1tiY/wRQSvkppT5VSu1XSv2slFqulLrXeixGKfWaUmoDMEQpVV8ptUIptV0p9btSqon1vGCl1Hxr3luVUp2t+2/LFo90p1LqBjc9DoMXYOqsh+DuuJPekoCLuew7B4SgK/Kr1n2lgG1APfSMnOXoP7Zq1vPvtZ4XA/w3W15rgIbWzx2AcOvnH4Au1s+1gb+sn5cCna2fywEB7n5GJnlWMnXW85JxKRSNrIhBfYAWWZYA+rWtIdAFmCvaZxavlFp7zfU/AiilygGdgLnWKeqgGwFAL6Bptv3lrZbBRuADpdT3wAIRiXXonRl8FVNn3YhRuHailAoDLEACuhI/IyIrrzlnQAHZpFi3fsB5sfrbrsEP6Cgil6/ZP1EptQw9n3uzUqqXiEQU9j4MxQdTZ92P8eHagVIqGJgBfCL6/WglOphHCevxRkqpssAGYLDVLxaCDoZxHSKSDBxRSg2xXq+UUi2th38Fns5WdlYnSH0R2Ssi76FfB5s44VYNPoKps56BUbi2E2h19u8HVqMr1XjrsS+BA8AOpRed+wz99jAfHTcza9+fQFIe+f8DeEQptRvYz9VlTf4FtLV2bBwAHrfuH6OU2mc9/zJ6fSWDITumznoYJlqYk1FKlRMdGT8I2ILuNIh3t1wGQ16YOus8jA/X+fyslKqIXgF1gqm4Bi/A1FknYSxcg8FgcBHGh2swGAwuwihcg8FgcBFG4RoMBoOLMArXYDAYXIRRuAaDweAijMI1GAwGF/H/8VsOgNIGJRQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 410.976x205.488 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "cost, lamb = 0.2, 4\n",
    "alphas, colors = [0.1, 0.3, 0.4, 0.7], [\"red\", \"black\", \"green\", \"blue\"]\n",
    "for index in range(len(alphas)):\n",
    "    func = lambda theta: unnorm_pr_fixed_beta(theta, alphas[index], cost, lamb)\n",
    "    xs = np.linspace(0.001,1,1000)[1:-1]\n",
    "    ys_grid = [func(x) for x in xs]\n",
    "#     print(ys_grid)\n",
    "    norm_const = integrate.simps(ys_grid, xs)\n",
    "    \n",
    "    xs = np.linspace(0,1,1000)\n",
    "    ys = np.array([func(x) for x in xs]) / norm_const\n",
    "    ax1.plot(xs, ys, color=colors[index], label=str(alphas[index]))\n",
    "\n",
    "legend = ax1.legend(title=r\"$\\beta=1, \\alpha=$\", frameon=False)\n",
    "move_legend(legend, ax1, 0.1)\n",
    "ax1.set_ylim(0,3)\n",
    "ax1.set_title(\"Cost=0.2\")\n",
    "ax1.set_xlabel(\"Degrees\")\n",
    "ax1.set_ylabel(\"Density of transition\")\n",
    "\n",
    "cost, lamb = 2, 4\n",
    "alphas, colors = [0.1, 0.3, 0.4, 0.7], [\"red\", \"black\", \"green\", \"blue\"]\n",
    "for index in range(len(alphas)):\n",
    "    func = lambda theta: unnorm_pr_fixed_beta(theta, alphas[index], cost, lamb)\n",
    "    xs = np.linspace(0.001,1,1000)[1:-1]\n",
    "    ys_grid = [func(x) for x in xs]\n",
    "#     print(ys_grid)\n",
    "    norm_const = integrate.simps(ys_grid, xs)\n",
    "    \n",
    "    xs = np.linspace(0,1,1000)\n",
    "    ys = np.array([func(x) for x in xs]) / norm_const\n",
    "    ax2.plot(xs, ys, color=colors[index], label=str(alphas[index]))\n",
    "\n",
    "ax2.legend(title=r\"$\\beta=1, \\alpha=$\", frameon=False)\n",
    "ax2.set_ylim(0,8)\n",
    "ax2.set_title(\"Cost=2\")\n",
    "ax2.set_xlabel(\"Degrees\")\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    for text in [\"top\", \"left\", \"right\"]:\n",
    "        ax.spines[text].set_visible(False)\n",
    "    ax.set_yticks([])\n",
    "\n",
    "w = 5.708\n",
    "fig.set_size_inches(w, w*0.5)\n",
    "plt.savefig(\"qingFranke2014Cost.png\", dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:computational_models] *",
   "language": "python",
   "name": "conda-env-computational_models-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}