AlexAndorra/pymc_constrained_priors.ipynb

## pymc_constrained_priors.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75ceaa3a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tutorial: Choosing priors under constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "276aeb4b",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "_**NB: you need PyMC >= 4.0 to run this notebook**_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d84adfc1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Problems as a modeler:\n",
    "- You know you want a Gamma prior (positive values) or a Beta ($[0, 1]$ values), but you don't have enough information to parametrize your priors efficiently. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e3c364e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- You need a non-stats person to define your prior. They neither know nor care about Gamma or Beta, but they give you important information:\n",
    "\n",
    "> I bet that we need to spend **between 0.1 and 3 dollars** on Google ads to get one additional sale -- I'm almost sure it's not less, but maybe it's more."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7034f989",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "As a statistician, you can translate that into:\n",
    "\n",
    "> **I need a Gamma that has 95% probability mass between 0.1 and 3**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67be2a18",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "**By hand, it may take you a loooooong time** to find an appropriate combination of parameters that give you a decent Gamma -- the drawback of the flexibility.\n",
    "\n",
    "Now we have two options:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff077aa9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "1. Despair and cry\n",
    "\n",
    "![cry](https://media.giphy.com/media/d2lcHJTG5Tscg/giphy.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4611f74e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "2. Not despair and use the new `pymc.find_constrained_prior` function\n",
    "\n",
    "![BringIt](https://media.giphy.com/media/sG4zmff2zDOp7t2MNA/giphy.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3339912",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "We'll demonstrate how to do option 1 in another video. For now, let's choose option 2!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1438e884",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Under the hood, PyMC is going to ask scipy's optimizer:\n",
    "\n",
    "> **What are the $\\alpha$ and $\\beta$ parameters that fit the user's constraint?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91709c74",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "import numpy as np\n",
    "import pymc as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "08a9b87f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (aesara.tensor.basic_opt): Optimization Warning: The Op gammainc_der does not provide a C implementation. As well as being potentially slow, this also disables loop fusion.\n"
     ]
    }
   ],
   "source": [
    "MASS = 0.95\n",
    "LOWER = 0.1\n",
    "UPPER = 3\n",
    "constrained_priors = pm.find_constrained_prior(\n",
    "    pm.Gamma, lower=LOWER, upper=UPPER, mass=MASS, init_guess=dict(alpha=1, beta=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "742b4369",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "We're getting back the optimized distribution's parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bf6fdbed",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'alpha': 2.2002514675808746, 'beta': 1.7605241820687383}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "constrained_priors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1196353",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**How do we use this in our PyMC model?** \n",
    "\n",
    "Easy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e1de28d5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "with pm.Model() as my_super_model:\n",
    "    my_beautiful_var = pm.Gamma(\"my_beautiful_var\", **constrained_priors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0548e848",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Now you may ask: what does this distribution look like? Go ahead, ask it! \n",
    "\n",
    "My answer: let's draw samples from it and plot them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "98047228",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "draws = pm.draw(my_beautiful_var, draws=10_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f83b193c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABxhklEQVR4nO3dd3gU9drG8e8kmyBhQ0joJRBaKCoikgCCDVHpigcRQcGj2BAUfTl2j+Wgx14AsetRmiBiowhiB4EAgqggSklCQughsLS0ef8YdiEkQDZkM1vuz3Xl2mFnd/aZnQD3/vaZ3ximaZqIiIiIiEi5C7O7ABERERGRYKWwLSIiIiLiIwrbIiIiIiI+orAtIiIiIuIjCtsiIiIiIj6isC0iIiIi4iMK2yIiIiIiPqKwLSIiIiLiIwrbIiIiIiI+orAtIiIiIuIjCtsiIiIiIj6isC3iZ7p27UqLFi3IyMiwuxSfCPb988aff/7JbbfdRnJyMi1btqRFixYsXbrU7rIkQP3111+0atWKf//733aXElRcLhdJSUlcd911dpciAcphdwEi5aFr165kZmYWuS8yMpKaNWvSrl07hgwZQps2bWyqTqS4Xbt2MWTIEHJycqhduzZNmjTBMAyio6PtLs3nsrKymDZtGosXLyY1NRWXy0VUVBTx8fGce+659OjRg/bt29tdZsB54YUXCA8P57bbbvPcZ5omK1as4JtvvmHFihVs3LiRQ4cOUa1aNc4991wGDx5Mx44dS/0a5b29k9m8eTOLFy9m9erVrF69mvXr11NQUMDdd9/N8OHDy7TNXbt28c477/Djjz+SkZFBYWEhtWvXplOnTgwbNoxGjRoVe47T6eSGG27gtddeY8GCBXTr1u10d01CjMK2BJWEhATi4uIAazQiLS2NL7/8kjlz5vD0009z1VVX2VugyBGzZ88mJyeHSy+9lPHjxxMWFvxfNBYUFPD666/zxhtvkJeXR1hYGHXq1CE+Pp69e/eybt06/vjjDyZNmkTbtm358MMPqVSpkt1lB4Tly5fzww8/cPXVV1O/fn3P/UuWLOHGG28EICwsjIYNGxIVFUVqairz589n/vz53HHHHYwaNapUr1Pe2zuZDz/8kA8//PC0t+O2ceNGrr/+enbt2kVERAQNGjQgIiKCtLQ0pk+fzpdffslbb71FcnJysecOHTqU9957j5deeolLL70UwzDKrS4JfgrbElRuu+02rr76as+fc3JyePTRR5k3bx5PPvkkl1xyCTExMTZWKGLZuHEjAF26dAmJoF1YWMgDDzzAF198QUREBCNGjGDAgAHUrl3b85i9e/fy1Vdf8frrr5Oamqqg7YVJkyYBFBtQME2TRo0aceONN9KrVy/Pv3+5ubmMHz+eN998k9dff51zzjmHSy655JSvU97bO5nY2FguueQSzj77bM4++2xmzJjBvHnzyry9J598kl27dtGuXTtefvll6tSpA0B2djYPPfQQ3377LQ8++CALFiwoFqZjYmK45JJLmDNnDkuWLKFTp06ntW8SWoL/X3gJaTExMTz99NNERUWxf/9+Fi1aZHdJIgAcPnwYgDPOOMPmSirGU089xRdffEHVqlX56KOPGDlyZJGgDVC1alUGDBjAl19+yS233GJTpYFn9+7dLFiwgFq1apGUlFRkXZs2bZgzZw6DBg0qMtAQGRnJvffey4UXXgjA9OnTS/Va5b29kxk+fDhvvPEGd955JxdeeCFRUVFl3tbBgwc950M8/vjjnqANVqh/5plnMAyDjIwMzwfh4/Xq1QuAjz/+uMx1SGjSyLYEPafTSUJCAmvWrCnxpLy//vqLr776ikWLFpGZmcmePXuoVq0abdu25aabbqJdu3YlbrdFixYArFu3jh9++IG33nqLNWvWEBYWxrnnnsu9995L69atS3xuZmYmL7/8MosWLeLAgQM0bNiQgQMHMmjQoJPuS3Z2Nm+//TbffPMNWVlZVKpUiZYtW3LNNdfQp0+fYqMxx9b49ddf8+6777Ju3ToqV67MBRdcwOjRo6lZsyYAn3zyCZMnT2bjxo1UrlyZyy67jH/9619l6iEuy/4dW+u8efP48MMPWbduHfv27eObb76hQYMGXh+rwsJCOnbsyL59+1iyZEmRcLBq1SquvfZaAMaMGcM111xT5Lnu8wDcr+3erzfffJNFixaxbds2IiIiiIuLo2XLlvTs2dPzn/HJjBs3jvHjx3v+/OCDD/Lggw8CkJyczMSJE0v9fkBg/E4sX76cSZMmERYWxmuvvcZZZ5110sc7nU6GDRtW7P5NmzYxb948Fi9ezObNm9m5cyfh4eG0aNGCgQMHnrBNLCkpib1797J06VJ+/fVX/ve///H7778TFhZG586defDBBz37/OWXXzJlyhTWrVtHdHQ0l112Gf/3f/9H5cqVfb7Nsu7j119/TV5eHhdeeGGxb0mcTudJ3+vOnTvz448/kpqaetLH+Wp7FSUvL4/CwkIA4uPji62PiYkhJiaGPXv2kJ+fX+I2unTpgsPhYMGCBeTm5hIZGenTmiV4KGxLSDh48CBAif+5Pf300yxevJiqVatSs2ZNatWqxZYtW/j666/59ttvefbZZ+nTp88Jtz116lSeeOIJatSoQePGjdm0aRM//fQTK1asYMaMGTRt2rTI4zds2MCgQYPYs2cPlSpVolmzZmRnZ/Pkk0+yfv36E75OWloaQ4cOJSsri4iICBITE9m7dy8pKSmkpKSwaNEiz+jM8SZOnMiYMWOoU6cODRs2ZOPGjXz22Wf8/vvvzJw5k+eff56JEycSHx9PgwYN2LRpE9OmTWPTpk18+OGHXvUnlnX/3N566y1efPFFatSoQUJCQpETX709VmFhYbRr147vvvuOFStW0LVrV8+6ZcuWFVk+Nmxv2bKFzMxM6tWr5wm1GRkZ9O/fn+zsbCpXrkzjxo0JDw8nKyuLBQsWkJGRUaqwXbduXdq1a0daWhq7du0qcp5BYmKiV+9HIPxOmKbJk08+CcA111xTYj9sab300kvMnz8fp9NJrVq1aN68Odu2bWPlypWsXLmS7du3c+uttxZ5zpYtW9i7dy81a9bkvffe48033yQ+Pp46deqwfv16Zs+ezZYtW5g4cSKjR49m3rx5NGnShNq1a7Nx40YmTpxIQUEBjz32mE+3eTr76P5dLstJ4OX9DYu/fmNTtWpV6tatS1ZWFitXrqRz585F1m/cuJE9e/ZQtWpVEhISStzGGWecQWJiImvWrGH16tU6iVdKzxQJApdccomZmJhofvLJJ8XWbdq0yWzdurWZmJhoLlu2rNj6uXPnmn/++WeR+woLC82vv/7abNu2rdmuXTtz3759xZ6XmJhoJiYmmuecc06R1923b585dOhQMzEx0Rw1alSx7fbr189MTEw0b7rpJjM7O9uzbtasWeaZZ57pqXXz5s1Fnnf11VebiYmJ5vXXX2/u2LHDs+6HH34w27ZtayYmJpqTJ08usca2bduaX375pef+rKws87LLLjMTExPN4cOHm+edd575888/e9b/+eefZnJyspmYmGh+//33xfb9RMq6f8fWeuaZZ5rTpk0zCwsLTdM0zby8PDMvL880zbIdq3feecdMTEw0n3nmmSL333LLLWbLli3Nzp07mxdffHGRdZ9++qmZmJho/utf//Lc9+STT5qJiYnm/fffb7pcriKPX79+vfnRRx+V+n0yTdO8//77T/g7a5qnfj8C5Xdi4cKFZmJiotmiRQszLS2t1M8ryeeff26uW7eu2P0LFiwwW7dubbZr187Mzc0tsu6bb77xvI9dunQxU1JSPOt+/vlns0WLFmZiYqJ58803m7179y6y/Y8//thMTEw027Rp43n/fbXN09nHrl27momJieZvv/12srevmMLCQvOqq64yExMTzSeffNKr51bE9o7n/jvz2muvlen57r/XF154ofnVV1+Zu3fvNvfu3Wv++OOPZvfu3c0WLVqY06ZNO+k2Hn30UTMxMdF88803y1SDhCb1bEvQcrlc/Pzzz9x5553k5+fTrl27Ekciunfv7vlq3c0wDLp168bQoUNxuVx89913J3yd/v37Fzkp0+l0etoCfvrppyKPXbJkCX/88QdnnHEGzz//PNWqVfOs69WrFwMHDizxK8zFixfz+++/ExkZyUsvvUSNGjU86y688ELuvPNOAN555x1M0yyxxt69e3v+XKdOHW6++WYAFixYwIgRI4qc8NOiRQsGDBhQ4j6cTFn371gDBw5kwIABnpFTh8OBw2F9CVeWY+XuYU1JSfHcV1hYyC+//ELLli3p0qULW7ZsKdJi5B4pPLb/1f21+I033kiVKlWKvEbTpk09LSnl7UTvR6D8TsyfPx+Atm3b0rBhw1I/ryR9+/YtcfT/0ksvpUmTJrhcLrZt21Zk3Z9//ulZfv3114sc006dOnm299tvv/H2228X2X7//v2Jiori0KFD7Ny506fbLOs+mqbJ1q1bATxtK6U1ffp01qxZQ0REBEOHDvXquRWxvfJ21VVXMW7cOGJjY7nrrrvo2LEj7du3Z9iwYURERPDWW295fsdPxP0eHz/VrMjJqI1Egsqx/a9uYWFh9OzZk8cff/yEz9uyZQuzZs3ijz/+IDs7m7y8PMA68Qis/1xP1ErSv3//Yve1aNGCSpUqsW/fPrKzs4mNjQWOhpTu3bt7WgeONWjQIE/P7rEWLlzoeV5J/6EOHDiQV199lczMTDZu3FisdaWkGlu1alWq9d5cfKas+3esK6+88qTrvT1WrVu3JioqirVr1+JyuXA6naxdu5Z9+/aRlJREixYt+PTTT1m2bJmnZWT58uVA0bBdt25dAObNm0eLFi0qbOqvE70fgfI78euvvwJw3nnnlbj+k08+4aGHHip2/7PPPlusP/ngwYPMnz+fFStWsGXLFg4cOOD5ILFp0yYAIiIiijzHHYwHDhxYYq+4+4PcHXfcUeSkOTd3D/SxbRG+2GZZ93Hv3r2eD7DHfrg9lT/++IOnnnoKgFGjRp32B6Hy3p4vmKbJ5s2b2bNnD+Hh4UWm/vv777+ZPn06bdq0Oen76D7vIzs7u4KqlmCgsC1Bxd3/apomO3fuZPPmzTgcDs4+++wTTvn36aef8thjj3l6DUuSk5NzwnUn+k8lLi6OrKwsDhw44Anb7tHR44PPsfU7HI5io7/u5zVr1qzE5zmdTurWrUtaWhqpqanFtl9Sje6a4uLiSjzpyR2W9+/fX+JrlqSs+3esEz0XynasHA4H7dq1Y+HChfzyyy9ceOGFnpHr5ORkz0h5SkoK/fr1Y/v27aSmplKzZs0ivZuDBw/ms88+Y8KECXz++ed06dKF9u3b06FDh2KzapSnE70fgfI7kZ6eDnDCPthDhw55Tmw1TZOVK1cCRYM/wA8//MBDDz1U4miwW6VKlahVq1aR+9atWwdAjx49ij3eNE3S0tIwDIPu3bsXW79r1y5cLhexsbFFTgr1xTbLuo/H/l04/oPGiWzevJnbbruNw4cP07t3b883GmVV3tvzlccee4xp06Zx7rnnMmnSJM+H6127dvHwww/z9ddfk56ezqeffkp4eHiJ23B/QDp06FCF1S2BT20kElRuu+02pk6dykcffcSCBQuYMmUKVapU4dlnn+Xzzz8v9vj09HQeffRRDh8+zE033cRnn33GihUr+PPPP1m3bh1jxowBOGk4PNF0VO7Rq2O/wj9w4ABwNNSU9JyS1rmfV9JosVv16tWBkoNQSSeGukdmS1p37HpvlHX/jnWi9/N0jpV7hNodslNSUjAMg/POO89zYpt73bFB/FitWrVi0qRJdOnShW3btjFt2jT+9a9/cdFFF3HzzTezYcOGk+5XWZ3o/QiU3wl3KDnRDCaDBw9m6tSpTJ061TMyGhkZWeTDwfLlyxk+fDjZ2dnccMMNTJkyhaVLl7JmzRrWrVvHq6++CkDz5s2L1Hjw4EHS09OJiIgo8eTBjIwMXC4XDRo0KHEE2j2C3bJlS59u83T28dhR2H379hV7vePt2LGDm266iR07dnDxxRef8ATa0irv7fnKn3/+yfTp04mIiODll1/2BG2w/p688MILxMbGsm7dOubOnXvC7ezZswc48b9xIiVR2Jagdt5553lC2NNPP43L5Sqyfu7cueTl5dGrVy/uv/9+WrVqhdPp9PxnkZWVVa71uIPTib6CLCws9PxjXtLz3K0SJdm1axdAsX7iilTW/SuN0zlW7l79ZcuWeS433bx5c89/mO3bt2fz5s1s3brV09t9/HzFYPUdv/vuu6SkpPDOO+9wyy23UKdOHRYuXMg///lP9u7dW6Z9K4tA+Z1wj5C7PxyczJo1awBrVhZ3KwZYM3Tk5+fz+OOP88gjj3DeeedRrVo1z+jjDz/8ABQfDf/rr78oLCykefPmJY76rl27tsTnubmD8bHrfbHN09nHyMhIz3t8sm/gwAqKN910E+np6SQnJ/Pqq6+WejS8IrbnSytWrMA0TRISEjwtYcdyOp2eD0+///77Cbfjfo9P9iFX5HgK2xL0unXrRtu2bdmzZw/vv/9+kXXuk1zOPffcEp977IlQ5cH9VfqJLpqQlpbm6UEu6XknmjrP5XJ5wuaJvq6vCGXdv9I4nWPVpk0bzjjjDH7//XdWr17Nnj17ioxcu5eXLl1aYr/28apUqeKZk3ru3Lk0bNiQbdu28eOPP3q9X2UVKL8T7jaX1atXn/Kx7rB97Pz0eXl5rFy5EsMwSpxjet++fXz77bdA8SB6olFkN3cwPtX6Y0/K9cU2T2cfj73vZN+u7N+/n1tvvZW//vqLs88+m9dff/20pucr7+35Wmlan9zfQp6sTc39Hp/oGgoiJVHYlpDgvhrdxIkTi/yj674cdEk9khs2bDjpLCRl0aVLFwC++uqrEkd/p0yZUuLzLrjgAs/zduzYUWz9tGnTyM3NpX79+jRp0qQcK/ZOWfevNE7nWEVGRtKmTRvy8vJ48803gaJh2j3yPW/ePDZs2EBcXNwJe6GPV7lyZc/sEdu3by/dzpSDQPmduPzyywHrwi7HzxRyPHcQPTbI5ObmUlhYiGma5ObmFnvO008/7fm25PiA6+6tPtUo84mCU0nP98U2T2cf4ejJpycakc3NzWX48OH8+uuvNG/enHfeeeeUF6c5mfLeXkVwf+BMTU0t8Vswl8vFb7/9VuSxJXE/RnNsizcUtiUkXHrppTRt2pScnBymTp3qud/9n9TUqVM9/9GDddb/qFGjyv0r0U6dOtG6dWsOHjzIfffdV+Rr3zlz5jB16tQiX5+7dezYkbPPPpvc3FzuvfdeT3sAWLNSuK9IeMstt9jaL1nW/SuN0z1W7tFr9wjhsWG7adOm1KhRg2+//RbTNEsc1X7ssceYM2eO5wJJbsuWLWPx4sVAxY52BcrvxIABA2jUqBEul4uhQ4cWuZiQ26ZNm3jppZc8J0ce+z5WqVLF87X/G2+84Rl93L9/P08//TSzZs0CrH7y46eFPNUo9IlaOsAKlBs3biQiIqLIhxVfbPN09hHwXKBlxYoVxdYVFBRwzz33sGTJEho2bMh7773n1awlzz77LF27duXZZ58tl+3NnDmTFi1aFLnAVHk7vmaw3qPY2Fjy8vK45557isyos2vXLkaPHk12djaVKlUq8cRWsL6Z27lzJ02aNCmxFUXkRDQbiYQEwzC46aabePjhh/nf//7HDTfcQKVKlTwtJqtWreIf//gHCQkJhIeH8/fff1OjRg3uuOMOXnnllXKt47nnnuP666/nxx9/5MILL/RcYTEzM5NBgwbxww8/FJvD1TAMXnzxRYYMGUJKSgoXX3wxzZs3x+VykZaWBlhTxA0cOLDcai2Lsu5faZzusXKPRJmmSdOmTT0nDx67/quvvgJKbiFZtWoVH330EQ6Hg0aNGlGlShV27drl2Ze+ffvSsWNHr/errALldyIqKoo333yT2267jU2bNnH99dcTFxdHvXr1yM3NJSsry3NiX6VKlRg0aFCxoHrHHXfw73//m7fffpsvvviCGjVqsHHjRsLCwrjvvvsYM2YMDRs2LDa66h5FLikY5+TksGXLFqpVq1biiYzr168nPz+f1q1bF/kg54ttns4+gvX72qhRI1JSUti5c2eROdfnzp3LggULAOsE5bvvvrvY88GaP3rs2LHF7nf/3XV/U3W62/PGihUrGD58uOfP7r7/t956iw8++MBz/2effVYk/B5fM+A5UX7kyJGsXLmSyy+/nPj4eBwOh6e9zeFw8MQTT5xwdqE5c+YA8I9//OO09ktCj8K2hIy+ffvy6quvsn37dmbMmMHgwYNxOBy8++67vPLKK8ybN4/09HSqV69O//79ueuuuzxzGZen5s2bM2PGDF555RUWLlzI33//TaNGjXj00UcZPHgwl156aYnPa9SoEZ9++ilvv/023377LX///TeRkZEkJSVxzTXX0LdvX7+YBaCs+3cqp3uszj33XCIiIsjLyysxTCcnJ580bD/44IN88803rFixgqysLNLT06lVqxZdunRh8ODBXHLJJWXar9MRKL8TjRs35rPPPmP69Ol8/fXXrF+/nj///JPKlStTu3ZtLr74Ys4//3wuvfTSEqfovPbaazl8+DAffPAB27ZtIzw8nL59+3Lbbbd55vE+Pvxu3rwZl8tF/fr1qVq1arFtnmwE+tj1x27XF9s8nX10MwyDa665hhdeeIE5c+YwZMgQz7pj21JSU1M9U0Yer379+iXef7zT3Z67DexE+3Ks/Pz8Ek+oPnjwYJFvmAoKCk65LYCLLrqIzz//nPfff58lS5awZcsWTNOkVq1atG/fnqFDh3LmmWee8PmzZ88mIiKCfv36ler1RNwMs6RLi4mIiEjAcLlcdOvWjZiYGObOneuZetTf3H777Xz33XdMmjTppCch+5slS5YwdOhQBg0axGOPPWZ3ORJg/PNvo4iIiJSa0+nkjjvuIDU1ldmzZ9tdzgmtWrWKs846K6CCNsBrr71GVFQUd955p92lSABSG4mIiEgQuO6663C5XBQWFtpdSok2btxIdnY2jzzyiN2leMXlctGhQweGDBlSpB9epLTURiIiIiIi4iNqIxERERER8RGFbRERERERH1HYFhERERHxEYVtEREREREfCdjZSI69MlRFiomJKXIJaql4Ogb20vtvPx0D++kY2E/HwH46BhAbG3vKx2hk20v+eqGAUKJjYC+9//bTMbCfjoH9dAzsp2NQOnqXRERERER8RGFbRERERMRHFLZFRERERHxEYVtERERExEcUtkVEREREfERhW0RERETERxS2RURERER8RGFbRERERMRHynQFydWrVzNu3DhWrVpFXl4ezZo1Y+jQofTp08er7bhcLt577z3mz5/P5s2biYiIID4+nksvvZQRI0aUpTQREREREb/hddheunQpN998MxEREfTq1Yvo6Gjmz5/P6NGjyczM5Pbbby/VdrZs2cLQoUPZvHkz559/PhdddBG5ubmkp6czb948hW0RERERCXhehe38/HweeeQRDMNg8uTJtG7dGoA777yTgQMHMm7cOLp3705CQsJJt1NQUMBdd93F9u3b+d///kfHjh2LvY6IiIiISKDzqmd7yZIlpKen07t3b0/QBnA6nQwfPpz8/Hxmzpx5yu3MmzeP3377jZtuuqlY0AZwOMrU3SIiIiIi4le8SrUpKSkAdOnSpdi6zp07F3nMycyZMweA7t27k5WVxffff8++ffuIj4/nwgsvpEqVKt6UJSIiIiLil7wK26mpqQA0atSo2LqYmBhiY2NJS0s75XZ+//13AFasWMF///tfcnNzPevi4uJ45ZVX6NChgzeliYiIiIj4Ha/CtsvlAiA6OrrE9U6nk61bt55yO7t27QJgzJgx3HTTTVx//fVERkYye/Zsnn32We68807mzJlDrVq1TriNmJgYwsLsmbkwNjbWlteVo3QM7KX33346BvbTMbCfjoH9dAxOzZbmaNM0Abj44osZPXq05/4bbriBbdu28fbbbzNjxgyGDx9+wm3k5OT4vM6SxMbGkp2dbctr+5Nt2002b4awMEhoBHFxRoW9to6BvfT+20/HwH46BvbTMbCfjkHpPmx4FbadTicA+/btK3G9y+U64aj38dvJzs6ma9euxdZdcsklvP32255WE/Evv6w0eec9k9W/Hb3PMKD9eSa3DjNo1bLiQreIiIiIv/OqD8M9pV9Jfdk5OTlkZ2eX2M99vMaNGwNQtWrVYuvc9x0+fNib0sTH8vNNXh1XyF33WEE7LAwaNYQGDcA0YdlyuPUOk7ffLaSw0LS7XBERERG/4FXYTkpKAmDhwoXF1i1atAiA5OTkU27HPd3f+vXri61z31e/fn1vShMfys83eehRk48/sf58ZR/4ZJrB5A/D+GhSGNOnGlxxuRW6P5gIT4wxyc9X4BYRERHxKmx36tSJ+Ph4Zs2axdq1az33u1wuJkyYgMPhoF+/fp77d+/ezYYNG9i9e3eR7Vx99dVERkYyadIktm3bVmQ7b775JgA9evQo0w5J+TJNk+dfMvl5MVSqBE89afCv/wujZs2j7SL16ho8+lAYDz9o4HDAN9/C4/8xKShQ4BYREZHQ5lXYdjgcjBkzBtM0GTRoEI8++ijPPvssV155JX///TcjRozwtIgATJ48mZ49ezJ58uQi24mPj+e+++5j165d9O3bl0ceeYQnn3ySvn37snbtWq699lo6depUPnsop2XmpzB7jtU28p/HDS668MQ92T2uMHj6PwYREfD9D/D2uwrbIiIiEtq8no2kY8eOTJkyhbFjxzJ37lzy8vJo1qwZd999N3379i31dm644Qbq16/Pu+++y+zZsykoKKBZs2bcfvvtDBgwwNuyxAf+XGcyboIVmIffbnB+p1Of/Hh+J4MH74MnnzKZNAUaJ5hccblOmhQREZHQZJjuefgCjF1TzYTKNDeHD5vcfJtJaipcdCGMecLAMEofmt98u5CJkyEyAiaMM2hZjrOUhMox8Fd6/+2nY2A/HQP76RjYT8egdFP/2XNVGPF7H06ygnZcLNz3f94FbYBbbjbo0hly86z+7QMHAvIznYiIiMhpUdiWYjK3mEz9yFq+d5RBTIz3o9JhYQYPPWBQqxZkZMIrYxW2RUREJPQobEsx4yeY5OZB+/OsFpKyqhpt8NgjBmFhMOcrWPCNAreIiIiEFoVtKWLZcpOfFkJ4GNw90vv2keOd08Zg6A3W8gsvm2zfrsAtIiIioUNhWzzy801eHW+F4av7QeOE8jmpcegNBq1bgcsFTz9r6gqTIiIiEjIUtsVj3teQmgpVq8JNN5bf7CEOh8EjDxlUqgTLV8DMz8pt0yIiIiJ+TWFbAMjNNXnvf9aI8/WDDKKjy3du7IbxBnfeYW1zwhsmaWka3RYREZHgp7AtAHz+JWzbBjVrwD/6+eY1+l0JyUmQmwv/edokP1+BW0RERIKbwrZw4IDJBxOt4HvjUINKlXxzxUfDMHjofoPoaPhzHZ7XFBEREQlWCtvCZ1/Anj3QoD706uHb16pRw2D0PVaY/3AirFmrwC0iIiLBS2E7xB0+bDJtuhV4h1xv4HD4ZlT7WJd2Neh2KRQUWu0khw4pcIuIiEhwUtgOcV/Nh127oVYtuKxbxb3uvaMMataAzZutEyZFREREgpHCdgjLzzeZMtUKugMHGERE+H5U261qtHU5d7CmAkxZpsAtIiIiwUdhO4T98CNkboGYqtCnV8W/flJ7g/5XW8tPP2uyd68Ct4iIiAQXhe0QZZomEydb4bb/PwwqV664Ue1j3X6rQcN42LkTXnxFYVtERESCi8J2iFqaAus3QOUzfDevdmmccYbBow8bhIfBN9/Cgm8UuEVERCR4KGyHqElTrFDbty9UrWrPqLZbq5YGQ4dYNbzwssmOHQrcIiIiEhwUtkPQX3+brPoVwsPh2v72Bm23IddDq5bgcln924WFCtwiIiIS+BS2Q9CMmVaQvfgiqFXLP8K2w2Hw6EMGlSrBsuXw6Wd2VyQiIiJy+hS2Q8yePSYLFljL/a/2j6Dt1rChwfDbrJomvGmSnq7RbREREQlsCtsh5svZkJsHiYlw1pl2V1Ncv6sgqT0cPmxdXTI/X4FbREREApfCdgjJzzf59PMj0/31MzAM/xrZBggLM3jofgOnE9b+CR9OsrsiERERkbJT2A4hPy+G7duhWgxc2tXuak6sZk2D/7vH+iDwwYcma//U6LaIiIgEJoXtEPLlbCu09uwJlSr536j2sS671ODSrlBQCP95yuTQIQVuERERCTwK2yFi+3aTpSnWcp9e/h203f5vlEGNGpC+Gd54S2FbREREAo/CdoiY8xUUFkLbcyC+QWCE7apVrf5tgBkzYcUvCtwiIiISWBS2Q0BhocnsOVZQ7R0go9puyUkGV11pLb/4sklengK3iIiIBA6F7RDwy0rI2grOKnDxhXZX473bbzGIi7XaST6abnc1IiIiIqWnsB0CZh0Z1e7WDc44I7BGtgGcToM777Dq/t+HJlu2FNhckYiIiEjpKGwHuYMHTRYuspZ79Qi8oO12+WVWv/nhw/DMc/vtLkdERESkVBS2g9zPi+HQIahXD1q2sLuasjMMg3tHGYSHwzff5fHzYvVui4iIiP9T2A5y33xnhdJLL8EvrxjpjSaNDQb0t5ZfHmty+LACt4iIiPg3he0gtn+/yZIl1nLXSwI7aLv9c6hB7VphZGXBpCkK2yIiIuLfFLaD2MJFkJsH8fHQrKnd1ZSPqCiD+++LAmDSFMjMVOAWERER/6WwHcSCqYXkWJd3iySpPeTlweu6sqSIiIj4MYXtILVvn0nKMms5WFpI3AzDYMRwg7Aw+P4H+HW1AreIiIj4J4XtILU0BfLzIaGRdWJhsGnaxKBXT2t5/OsmhYUK3CIiIuJ/FLaD1OIlVvg8v5PNhfjQsH8aVK4Ma9fCgm/trkZERESkOIXtIFRQYLJkqbV8fqfgG9V2q17d4IbB1v698ZamAhQRERH/o7AdhNb+CTl7wemEs860uxrfuvYaqFULtm+HaR/bXY2IiIhIUQrbQWjRkasrdkgChyN4R7YBKlUyuP0Wax8nTzXJydHotoiIiPgPhe0gtHixddspiFtIjtXtUmse8f37YfJHCtsiIiLiPxS2g8z27SbrN4BhQIdku6upGGFhBrfcbH2w+GQm7NylwC0iIiL+QWE7yCw+cmLkma0htlpojGyDNevKWWfC4cPw4USFbREREfEPCttBZtlyK2h27BA6QRusC93cOsza5y9mQVaWAreIiIjYT2E7iJimyapV1vJ57WwtxRbtzjVof551MZ/3P1DYFhEREfspbAeRTamwJwfOOANatrC7Gnu4R7e/mg9paQrcIiIiYi+F7SCycpV1e/ZZEBERWm0kbq1bGVzQGQoL4X/q3RYRERGbKWwHkZWrrHB5btvQDNpu/7zR2v9vvoXNGQrcIiIiYp8yhe3Vq1dzyy23kJSURNu2benfvz9ffvllqZ+/dOlSWrRoccKfVe7GYym1Y/u1z21rZyX2S2xucH4na3R74iSFbREREbGPw9snLF26lJtvvpmIiAh69epFdHQ08+fPZ/To0WRmZnL77beXelvJyckkJxefDLpOnTrelhXy1K9d1I1DDH5ebDJvPgwdYlK/XmiP9ouIiIg9vArb+fn5PPLIIxiGweTJk2ndujUAd955JwMHDmTcuHF0796dhISEUm0vOTmZkSNHel20FKd+7aJatzJITjJJWQaTppjcP1rviYiIiFQ8r9pIlixZQnp6Or179/YEbQCn08nw4cPJz89n5syZ5V6knJr6tYu7cYj1Xsz9CrZuUzuJiIiIVDyvwnZKSgoAXbp0Kbauc+fORR5TGqmpqXz44Ye89dZbzJo1i927d3tTjhxhmiarfrWWQ71f+1htzjZod6417/bkqQrbIiIiUvG8aiNJTU0FoFGjRsXWxcTEEBsbS1paWqm3N2vWLGbNmuX58xlnnMHIkSMZNmyYN2WFvMxM2LMHIiPUr328G4cY/LLSZNZsGHqDSY3qGvkXERGRiuNV2Ha5XABER0eXuN7pdLJ169ZTbicuLo777ruPiy++mHr16rF3716WLl3KCy+8wPPPP4/T6WTgwIEn3UZMTAxhYfbMXBgbG2vL657IwkWHARetWzuoVSvG7nIqRGmPQddLTNq23cuqVfl8OasS946q4uPKQoO//R0IRToG9tMxsJ+Ogf10DE7N69lIykPz5s1p3ry558+VK1emb9++tGzZkquvvppx48YxYMCAk4bpnJyciii1mNjYWLKzs2157RNZtrwQgBaJ+X5Xmy94ewwGXmNNi/jR9ENc84/DOJ0a3T4d/vh3INToGNhPx8B+Ogb20zEo3YcNr4aGnU4nAPv27StxvcvlOuGod2kkJiZyzjnnsHPnTq/aUULdH2us29atFCJLcn4nSEiA/fvhsy/srkZERERCiVdh2z2lX0lBOCcnh+zs7BL7ub3h/oRw6NCh09pOqDh82OTv9dbyma1P/thQFRZmMHig9UFk+scmhw/rZEkRERGpGF6F7aSkJAAWLlxYbN2iRYsASrxITWnl5+ezZs0aDMOgbt26Zd5OKPnrbygogLhYqF3b7mr8V7dLoVYt2J0N8762uxoREREJFV6F7U6dOhEfH8+sWbNYu3at536Xy8WECRNwOBz069fPc//u3bvZsGFDsSn9Vq5ciWkWHV3Mz8/nueeeIzMzky5dulCtWrUy7E7oWXPkMJzZGgxDbSQnEhFhMKC/9f58PMMs9vsnIiIi4gtenSDpcDgYM2YMw4YNY9CgQfTu3Run08n8+fPJyMhg1KhRNG7c2PP4yZMnM378eEaMGFHkSpH/93//B8C5555L7dq12bdvH8uWLWPTpk3Uq1ePJ554opx2L/j9scYKja1bK2ifSu+e8O771qXtl6+ApPZ2VyQiIiLBzuvZSDp27MiUKVMYO3Ysc+fOJS8vj2bNmnH33XfTt2/fUm1j4MCB/PTTT6SkpJCdnY3D4aBhw4bcfvvt3HTTTcTEhMb0deVhzZGTI9WvfWpOp0GvHiYzZlqj20nt9QFFREREfMswA/T7dLummvGnaW527TK58h8mhgHzZhtERYVGeDydY7A5w2TQDSamCVMmGjSMD433rDz509+BUKVjYD8dA/vpGNhPx8AHU/+Jf3H3azdpTMgE7dMV38CgU0dr+ZOZAfk5U0RERAKIwnYAW/eXFRZbtrS5kADjPlFyzlzYt0+BW0RERHxHYTuA/fW3dZvYTKPa3jivnfVtwMFDMHuu3dWIiIhIMFPYDmB/Hwnbx1z5XkrBMAyu+Yf1AWXGJyYFBRrdFhEREd9Q2A5Q2XtMduwEw4BmTe2uJvBcfhlUrQpbt8HSFLurERERkWClsB2g3KPa9evr5MiyqFTJoGd3a/mzzzWyLSIiIr6hsB2gPP3aaiEpsyv7Wh9SFi+FrVsVuEVERKT8KWwHqL/XW+GwuU6OLLP4BgbtzwPThC9mKWyLiIhI+VPYDlB/a2S7XFx1ZHR71mzIz1fgFhERkfKlsB2ADhww2ZxhLTdvZm8tga5LZ6geB7uz4ceFdlcjIiIiwUZhOwBt2Gi1PlSvDnFxaiM5HQ6HQe9e1vLnX2hkW0RERMqXwnYA0smR5atPb4OwMFjxC6SnK3CLiIhI+VHYDkB//+0+OdLmQoJEndoGnTpay5/rREkREREpRwrbAejv9datZiIpP317W+/lvPk6UVJERETKj8J2gCksNElNs5ab6sqR5aZDsnWi5J498PNiu6sRERGRYKGwHWCysuDwYYiMgPr17K4meDgcBldcbi3P+Uoj2yIiIlI+FLYDzMZN1m1CAoSHq42kPPXsceSKkoth924FbhERETl9CtsBxh22GyfYWkZQSmhk0LoVFBTC/AV2VyMiIiLBQGE7wGxKtUZcGzfWqLYv9Oxuva9z5pqYpka3RURE5PQobAeYTUdGtps0treOYHVpV4iMtL5BWLfO7mpEREQk0ClsB5D8fJO0dGtZYds3oqMNLrzAWp6tEyVFRETkNClsB5DNGZCfD5UrQ+3adlcTvHodOVHy6wVw+LACt4iIiJSdwnYA2ZRq3TZOAMNQz7avtDsXatUClwsWLrK7GhEREQlkCtsBZNMma5RVLSS+FR5u0OMKa1lzbouIiMjpUNgOIJ5p/zQTic/1uMJ6j5cthx07FLhFRESkbBS2A4hmIqk4DRoYnNMGCgvhq/l2VyMiIiKBSmE7QBw+bJKRaS03VtiuED2OzLk99yvNuS0iIiJlo7AdINLTrVHW6GioHmd3NaGh68VwxhmQvhnWrLW7GhEREQlECtsBwj2/dkIjzURSUaKiDC46Muf2XJ0oKSIiImWgsB0g0tKtsNeooc2FhBh3K8mCbzXntoiIiHhPYTtAuEe2GzXSqHZFKjLn9s92VyMiIiKBRmE7QHjCtka2K1RYmEH3y63lr+ZpZFtERES8o7AdAAoKTDZvtpYVtite9yNzbqekwM5dCtwiIiJSegrbAWDbdsjNhcgIqFPH7mpCT8N4g7POhIJCmP+13dWIiIhIIFHYDgBpadZtfLx1KXGpeO4TJb+apzm3RUREpPQUtgOAu1+7oVpIbNP1YuubhY2bYN1fdlcjIiIigUJhOwBo2j/7RUcbXHBkzm2dKCkiIiKlpbAdANI17Z9f6HHkRMmvF0BengK3iIiInJrCdgBw92xrZNteSe2henXI2Qs/L7G7GhEREQkECtt+bs8ekz051nJ8A3trCXXh4QZXuOfc1uXbRUREpBQUtv1c+pH5tWvXhsqV1UZiN3cryc9LIHuPAreIiIicnMK2n0vXlSP9SuMEg1YtoaDA6t0WERERORmFbT+XmqaZSPyN+4qSmpVERERETkVh289tzrBu4+PVQuIvunWFiAj4629Yv0GBW0RERE5MYdvPZbjDtk6O9BsxMQadO1nLczW6LSIiIiehsO3HCgpMtmRZywrb/qX7kcu3f/015OcrcIuIiEjJFLb92LbtkJdnXSa8Vi27q5FjdUyG2FjYnQ1LU+yuRkRERPyVwrYfc7eQ1KsPYWHq2fYnDofB5d2sZbWSiIiIyIkobPsxz8mR9e2tQ0rW40gryaKfYe9eBW4REREprkxhe/Xq1dxyyy0kJSXRtm1b+vfvz5dfflnmIvLy8rjyyitp0aIF3bt3L/N2gk1GhhXgGqhf2y81a2rQvJnV6rPgW7urEREREX/kddheunQpgwYNYvny5VxxxRVcd911ZGdnM3r0aN54440yFTFhwgTS3VdvEY+MTOu2QQO1kPgr9xUl1UoiIiIiJfEqbOfn5/PII49gGAaTJ09mzJgx3H///Xz++ec0b96ccePGkZqa6lUBf/zxB2+99Rb33nuvV88LBZs17Z/fu6wbhIfD2rWwcaMCt4iIiBTlVdhesmQJ6enp9O7dm9atW3vudzqdDB8+nPz8fGbOnFnq7eXm5vLAAw9wzjnncP3113tTStDLzzfJOjLtXwP1bPut2FiDCzpby1/MVtgWERGRorwK2ykp1hxnXbp0Kbauc+fORR5TGuPHjyctLY2nnnoKw1CrxLG2boWCAqhUCWrUsLsaOZm+fdyXb4fDhxW4RURE5Civwra7RaRRo0bF1sXExBAbG0taWlqptrV69WreeecdRo4cSePGjb0pIyRsdvdra9o/v9f+PKhbB1wu+O4Hu6sRERERf+Lw5sEulwuA6OjoEtc7nU62bt16yu3k5uby4IMP0qpVK2666SZvSvCIiYkhLMyemQtjY2N9/hq7dx0EDtC4cSSxsSW/36GsIo6BNwZcc4BXxx1kztxwrrs2xu5yfM7f3v9QpGNgPx0D++kY2E/H4NS8Ctvl5ZVXXiEtLY1PPvmE8PDwMm0jJyennKsqndjYWLKzs33+On/9XQhA7Vq5FfJ6gaSijoE3LrnYZPxr8MvKfH5ZuZvGCcH7bYQ/vv+hRsfAfjoG9tMxsJ+OQek+bHg1NOx0OgHYt29fietdLtcJR73d/vjjD/73v/9x++2306JFC29ePqQcnYkkeENbMKlR3eDIaQt8OUt92yIiImLxKmwnJCQAlNiXnZOTQ3Z2don93Mdat24dBQUFjBs3jhYtWhT5Adi0aRMtWrSgffv23pQWdI7OsW1vHVJ6fXu759zWiZIiIiJi8aqNJCkpiTfffJOFCxfSq1evIusWLVoEQHJy8km3kZCQQP/+/UtcN2PGDKKjo7niiiuoXLmyN6UFlbw8E3fru8J24EhqD3Vqw9Zt8MOPcPlldlckIiIidvMqbHfq1In4+HhmzZrFkCFDaNWqFWC1j0yYMAGHw0G/fv08j9+9ezfZ2dnExsYSFxcHQLt27WjXrl2J258xYwY1atTgqaeeKuv+BIWsLCgshMqVoXqc3dVIaYWHG/TuBe+8Z/LFLJPLL1MLkIiISKjzqo3E4XAwZswYTNNk0KBBPProozz77LNceeWV/P3334wYMaLINH6TJ0+mZ8+eTJ48udwLD2bufu0G9dH84wGmVw8ID4NVv0JamlpJREREQp3Xc+d17NiRKVOmcN555zF37lymTJlCtWrVeP7557njjjt8UWPIUb924KpZ0+D8TtayrigpIiIiZZr6r02bNrzzzjunfNzIkSMZOXJkqbe7bt26spQTdDZnWCFNYTsw9e1j8NMik6++gltvNqlUSd9OiIiIhCp7rgojJ5XhnvavvkJaIEpOglq1IGcv/LjQ7mpERETETgrbfsgdtjWyHZjCww369LI+KH3+hVpJREREQpnCtp85fNhk23ZrOV5hO2D16gFhR06UTE9X4BYREQlVCtt+ZksWmCZUqQLVqtldjZRVrVoGnTpay1/oipIiIiIhS2Hbz2RusW7r19O0f4HOc0XJryA3V4FbREQkFCls+5ktR8J2vXr21iGnr0My1Kp55ETJn+yuRkREROygsO1ntmyxRkDrK2wHPIfDoFdPa1mtJCIiIqFJYdvPHB3ZVgtJMOjV0yAsDH5ZCembFbhFRERCjcK2n9mSZd3Wq2tvHVI+6tQ26JhsLc/SFSVFRERCjsK2HyksND0j2/Xr21uLlJ++faxvKebM1YmSIiIioUZh24/s2gW5eRAeZp1YJ8GhYweoWQP25MBPi+yuRkRERCqSwrYfcbeQ1K5jnVwnwcHhMOjdy1r+4kuNbIuIiIQShW0/kplp3WomkuDTq6eBYcCKXyAjQ4FbREQkVChs+5EtWVYI08mRwadObYOOHazlL3WipIiISMhQ2PYjmZr2L6i5ryg55yvIy1PgFhERCQUK235kyzGXapfg06kjVK8O2dmwUCdKioiIhASFbT/imWNbYTsoORwGvXVFSRERkZCisO0nDhwwyc62ltWzHbx6HzlRctlyyNyiwC0iIhLsFLb9hHtUO6YqOJ3q2Q5WdesaJCdZy7qipIiISPBT2PYTWzwnR9pbh/ie+4qSs+dAfr4Ct4iISDBT2PYTnplI1EIS9Dp3gupxsDsbFv1sdzUiIiLiSwrbfsIzx3Z9mwsRn3M4DHoeOVHyc11RUkREJKgpbPsJz7R/ddWvHQr69LKO87LlRz9oiYiISPBR2PYTmerZDin1jpwoaZo6UVJERCSYKWz7gYICk61brWWF7dDhvqLk7Lk6UVJERCRYKWz7gR07IT8fHA6oWcPuaqSidOkMcbGwaxf8vNjuakRERMQXFLb9gLtfu04dCA9Xz3aocDgMevawlnVFSRERkeCksO0HPCdHqoUk5LhPlFyaAlu3KnCLiIgEG4VtP+C+bLf6tUNP/foG7c87cqLkHIVtERGRYKOw7Qc8V4/UtH8hyX1FyVm6oqSIiEjQUdj2A1uyrFuNbIemCzpDtWqwcycsWWp3NSIiIlKeFLb9QJY7bOtS7SEpIuKYEyV1RUkREZGgorBtswMHTPbkWMt1attbi9jHfaLkkhTYuk2BW0REJFgobNts6zbr1umE6Gj1bIeq+AYG57WDwkKYrRMlRUREgobCts08V45UC0nI81xRUidKioiIBA2FbZu5+7Xr1LG3DrHfBV2gWgxs32HNuy0iIiKBT2HbZllHLmRSV2E75EVGGvTobi3ripIiIiLBQWHbZu42kjp11K8t0OdIK8niJbB9uwK3iIhIoFPYtlnWkbCtkW0BaBhvcG5b60TJz75Q2BYREQl0Cts2c49s19UJknJE/6ut0e2Zn8H+/QrcIiIigUxh20YHDpjk7LWWNce2uF3QBRIagctlBW4REREJXArbNnK3kFStClWqqGdbLGFhBtcPtn4fpn1scuiQRrdFREQClcK2jbI8J0faW4f4n25drT7+PXtg1hy7qxEREZGyUti2keeCNgrbchyHw2Dwddbo9uQpGt0WEREJVArbNnLPsa2RbSlJj+5WL/+OnTB9ht3ViIiISFkobNvIffXIuppjW0pQqZLBrbdYvxsTJ5vs3q3RbRERkUCjsG2jrerZllPo1hVatoCDB+G9DxS2RUREAo3Cto10QRs5lbAwgzvvsEa3v/wSNm5S4BYREQkkCts2cblM9u2zljWyLSdzbluDC7pAQSE8/6JJYaECt4iISKBwlOVJq1evZty4caxatYq8vDyaNWvG0KFD6dOnT6mev3TpUqZPn86aNWvYsWMHeXl51KlTh3bt2nHLLbfQpEmTspQVUNyj2tViICpKPdtycqPuMljxi8lvv8Onn8M/+tldkYiIiJSG1yPbS5cuZdCgQSxfvpwrrriC6667juzsbEaPHs0bb7xRqm38/PPPrFixgubNm9OvXz8GDx5M48aN+fzzz7nyyitZsmSJ1zsSaNSvLd6oXcvg9lutD2VvvGWybbtGt0VERAKBYZpmqf/Xzs/Pp0ePHmzdupVp06bRunVrAFwuFwMHDmTTpk3Mnj2bhISEk27n8OHDVKpUqdj9ixcv5sYbb+Sss87ik08+Oek2srOzS1t2uYqNjS2X154+w2TseJOLL4IxT6ibxxvldQwCTWGhyZ13WaPbyUnwwrMGYWEV/61IqL7//kTHwH46BvbTMbCfjoH1HpyKVylvyZIlpKen07t3b0/QBnA6nQwfPpz8/Hxmzpx5yu2UFLQBOnXqRExMDOnp6d6UFZC2HpljWydHSmmFhRnc/y+DyEhIWaa5t0VERAKBV2E7JSUFgC5duhRb17lz5yKPKYuVK1eSk5ND8+bNy7yNQOGZiaSu+rWl9BIaGdx159F2knV/qZ1ERETEn3l1gmRqaioAjRo1KrYuJiaG2NhY0tLSSr29pUuXkpKSQm5uLmlpaXz33XfExsby4IMPelNWQDp6QRt765DAc2VfWLoMfloIj//H5N03dZKtiIiIv/IqbLtcLgCio6NLXO90OtnqPvOvFFJSUhg/frznz40aNeKll17irLPOOuVzY2JiCAuzp9e5NP05p7Jt+27AJDGxKrGxZZoUJqSVxzEIZM88Xcg/rslh8+ZCxr0WwX+fcmIYFRe4Q/399wc6BvbTMbCfjoH9dAxOzdaUN3LkSEaOHMmBAwdYv349EyZM4LrrruPpp58+5TSCOTk5FVRlUeVxMsC+fSb79llf/0dV3kt2tkYlvaETMiyPPmxy1yj4cnYuZ7bOpnevivk90vtvPx0D++kY2E/HwH46Bj44QdLpdAKwz301luO4XK4TjnqfTFRUFG3atGH8+PE0adKEf//73+zevdvr7QQK9+B/tWpQubKCtpTNOW0Mht1s/f68PNZk40b1b4uIiPgbr8K2e0q/kvqyc3JyyM7OLrGfu7QcDgcdOnTgwIED/Pbbb2Xejr87enKkvXVI4Bt8nTUN4OHD8O8nTA4eVOAWERHxJ16F7aSkJAAWLlxYbN2iRYsASE5OPq2Ctm/fDljBO1h5wrZOjpTTFBZm8OjDBjVqQGoavPSKwraIiIg/8Spsd+rUifj4eGbNmsXatWs997tcLiZMmIDD4aBfv6PXkd69ezcbNmwo1hKybNkySrqWzsKFC1mwYAHR0dGce+653u5LwMg6Mse2rh4p5SG2msHjjxqEhcHceTBnrgK3iIiIv/Bq+NjhcDBmzBiGDRvGoEGD6N27N06nk/nz55ORkcGoUaNo3Lix5/GTJ09m/PjxjBgxgpEjR3ruv+OOO4iNjeXss8+mTp06HD58mHXr1rFs2TIiIiIYM2YMUVFR5beXfmarZ2Rb/dpSPtqeY3DzP+Htd01eetWkVStonKDfLxEREbt53avRsWNHpkyZwtixY5k7dy55eXk0a9aMu+++m759+5ZqGyNHjuSnn35ixYoV7N69G8MwqFu3Ltdccw1Dhw4N+ovauNtINLIt5emGwbDqV1i2HB593Jp/u1IlBW4RERE7GWZJ/RwBwK6pZk53mhvTNOne22T/fpj8gUGjRgpD3tJUQyeWnW1y480mu3bDddfCnXeU/1z0ev/tp2NgPx0D++kY2E/HwAdT/8np27cP9u+3ljWyLeUtNtbgvtHWB7hpH8PvfwTkZ2kREZGgobBdwdwtJHGx+opffKPz+Qbdr4DCQnj6GZPDhxW4RURE7KKwXcG2ql9bKsBdIwyqV4f0zfD+BwrbIiIidlHYrmC6oI1UhKrRBqPvsb45mToNNm5S4BYREbGDwnYFy8qyQo8uaCO+dkEXgy6doaAAXnzZLHFuexEREfEthe0KdnTaP/Vri++NGmlwxhnw62qY+5Xd1YiIiIQehe0KtlWXapcKVKeOwU03Wh/sJrxhsm+fRrdFREQqksJ2BTJN82jPtsK2VJAB/SGhEezJgQ8mKmyLiIhUJIXtCrR3Lxw8aC3Xrm1vLRI6HA6DEcOt0e0ZM2FzhgK3iIhIRVHYrkBZWdZt9eqaY1sqVscOBh07QH4+THhdYVtERKSiKGxXILWQiJ1GDDcID4OfFsEvKxW4RUREKoLCdgXK0gVtxEYJjQyu7Gstv/6WpgIUERGpCArbFWjrVs2xLfa6cYhB5TNg7Vr4caHd1YiIiAQ/he0K5Jn2r676tcUecXEGA66xlt9+x6SgQKPbIiIivqSwXYG2qGdb/MB11xpUrQqpaTBvvt3ViIiIBDeF7QpimqZnZFs922Inp9PghsHWtyvv/s/k8GGNbouIiPiKwnYF2ZMDhw5Zy7Vr2VuLyNVXQc0asG0bfP6F3dWIiIgEL4XtCuIe1a5RAyIj1bMt9qpU6ehl3D+cZHLggEa3RUREfEFhu4Jojm3xNz26Q4P61rcun2l0W0RExCcUtiuI++qRCtviLxwOgyHXW6PbU6eZHDqk0W0REZHyprBdQbKOzLGtkyPFn1x+GdSrB9nZGt0WERHxBYXtCuKZY7uO+rXFfzgcBkOOzEwyZapGt0VERMqbwnYFOXpBG3vrEDle9yus9qbd2fDFLLurERERCS4K2xXANE2dICl+y+EwuOFI7/bkqZp3W0REpDwpbFeA7Gw4fBgMA2ppjm3xQz2ugNq1YdcumDXb7mpERESCh8J2BdhyZCaSmjUhIkI92+J/IiIMrh90dGaS/HyNbouIiJQHhe0K4A7b9dSvLX6sZ3eoVg22boPvfrC7GhERkeCgsF0BNMe2BIJKlQz6X32kd3uKiWlqdFtEROR0KWxXgKwsK7TUrasWEvFvV18FZ5wB6zfAsuV2VyMiIhL4FLYrQJam/ZMAUbWqQZ9e1vLkqRrZFhEROV0K2xVAPdsSSK69xiA8DFb8An+uU+AWERE5HQrbPpafb7J9m7Wsnm0JBHXqGFza1Vqe+pHCtoiIyOlQ2PaxHTugoBAiIqBGDburESmd6wZa5xd89wNkblHgFhERKSuFbR9z92vXrg1hYTpBUgJD82YGyUlQWAjTpitsi4iIlJXCto+pX1sC1eDrrA+Hs+dC9h4FbhERkbJQ2PYxz7R/6teWANPuXGiRCIcPw8xPFbZFRETKQmHbxzwXtNEc2xJgDMNg0JHe7U8/g0OHFLhFRES8pbDtY5pjWwLZRRda38rsyYG58+yuRkREJPAobPuYerYlkDkcBtcOsEa3P5pmUlCg0W0RERFvKGz70OHDJrt2Wcvq2ZZA1asHREdD5hZYuMjuakRERAKLwrYPbT3SQlK5MsTE2FuLSFlVrmzQ70precpHJqap0W0REZHSUtj2oS1Hwna9utbJZiKBqv/VBpER8Mca+GVlvt3liIiIBAyFbR86OhOJvXWInK64OIMrrrCW3//goL3FiIiIBBCFbR/SHNsSTAZec+QS7t/nkZ6uVhIREZHSUNj2Ic2xLcGkUSODLp2t5am6hLuIiEipKGz70LE92yLB4LprrQ+O8+bBrl0K3CIiIqeisO1D6tmWYNPmbGh7joPcPJg6TWFbRETkVBS2fcTlMtm3z1pWz7YEC8MwuP3WygB89gVk71HgFhERORmFbR9xj2pXi4GoKPVsS/Do0jmCli3g0CGYpt5tERGRkypT2F69ejW33HILSUlJtG3blv79+/Pll1+W+vnLly/nmWee4eqrr6ZDhw6cffbZdO/eneeff569e/eWpSS/4+7XVguJBBvDMLhxiPUB8pNPISdHgVtEROREvA7bS5cuZdCgQSxfvpwrrriC6667juzsbEaPHs0bb7xRqm3cfffdfPjhh1SpUoUrr7ySQYMGUblyZd555x3+8Y9/sMt9jfMApn5tCWadz4dmTeHgQfjoY4VtERGRE3F48+D8/HweeeQRDMNg8uTJtG7dGoA777yTgQMHMm7cOLp3705CQsJJtzN06FCuuuoqatWq5bnPNE2eeOIJpk6dyvjx43nssce83xs/ojm2JZgZhsE/b4SHHzX5eAZcc7VJXJzapURERI7n1cj2kiVLSE9Pp3fv3p6gDeB0Ohk+fDj5+fnMnDnzlNu59dZbiwRtsP7zHj58OADLli3zpiy/lLnFuq1XTwFEgtOFXaBVK6t3+4OJGt0WEREpiVdhOyUlBYAuXboUW9e5c+cijykLh8MaaA8PDy/zNvxFRoZ1G9/A3jpEfMUwDO641fow+fmXkLlFgVtEROR4XoXt1NRUABo1alRsXUxMDLGxsaSlpZW5mE8++QQ4GtwDVX6+6enZrl/f3lpEfKnduQbJSZCfD+++p7AtIiJyPK96tl0uFwDR0dElrnc6nWzdurVMhaxdu5bXXnuN6tWrM2zYsFM+PiYmhrAwe2YujI2NPen6tPQCCgr3UKkSJDaPJSxMrSTl7VTHQHzr2Pd/9L35DLguh/kL4MahTs5pE2FjZaFDfwfsp2NgPx0D++kYnJpXYdtXNm/ezG233UZBQQEvvfQScXFxp3xOTk5OBVRWXGxsLNnZ2Sd9zB9rrBG++vUgJ2dPBVQVWkpzDMR3jn//69WFnt1hzlfw+H/28vbrBuHh+oDpS/o7YD8dA/vpGNhPx6B0Hza8Ghp2Op0A7HNfGvE4LpfrhKPeJ5KZmcnQoUPZvXs3Y8eOpWPHjl493x+5+7UbqF9bQsQdtxk4q8Bff8EXs+yuRkRExH94FbbdU/qV1Jedk5NDdnZ2if3cJ5KRkcENN9zA9u3beeWVV7jkkku8KcdvZWRYI9sK2xIqYmMNbhlmjWa/+bZJdrb6t0VERMDLsJ2UlATAwoULi61btGgRAMnJyaXaVkZGBkOGDGH79u28/PLLdOvWzZtS/Npm90wk9fVVuoSOq/pCYnNwueDlsQrbIiIi4GXY7tSpE/Hx8cyaNYu1a9d67ne5XEyYMAGHw0G/fv089+/evZsNGzawe/fuIttxB+1t27bx0ksvcdlll53mbviXjEzrViPbEkrCww3uG20QHgbffgfffKvALSIi4tUJkg6HgzFjxjBs2DAGDRpE7969cTqdzJ8/n4yMDEaNGkXjxo09j588eTLjx49nxIgRjBw50nP/kCFDyMzMpG3btqxbt45169YVe61jHx9I8vJM3BOyaI5tCTUtWxgMucHk/Q/gxVdM2p4D1avrGx4REQldXs9G0rFjR6ZMmcLYsWOZO3cueXl5NGvWjLvvvpu+ffuWahuZmdbQ76pVq1i1alWJjwnUsJ21FQoL4YwzoHp1u6sRqXhDrjdY9LPJX3/Ds8+bPPtf6wI4IiIiocgwTTMgv+u1a6qZU01z8/Nik/seNGnaFD541555wIOdphqyV2ne/w0bTW65zSQ3D4bfbjBooMJ2edLfAfvpGNhPx8B+OgY+mPpPTs1zmXZdOVJCWNMmBneNPDI7yVsmv/0ekJ/pRURETpvCdjnbnKlp/0QAruwDl3aFgkJ47EmT7D0K3CIiEnoUtsvZ0Qva6GtzCW2GYXD/aIP4eNi+HR75t0lengK3iIiEFoXtcuae9k8zkYhAVJTBf8cYVKkCv66Gl14xCdDTRERERMpEYbsc5eaabNtmLTdQz7YIAAmNDB5/1MAw4MvZMGOm3RWJiIhUHIXtcpSRYU37V6UKxMXZXY2I/+jU0WD47VZr1bjXTFKWaXRbRERCg8J2OdqUZt0mNNK8wiLHGzgAena3PpA++rjJxk0K3CIiEvwUtstRmjtsJ9hahohfMgyD0fcatDkb9u+Hfz1gsmuXAreIiAQ3he1ylJpmBYdGDTWqLVKSyEjrhMkGDWDbNrjvQZODBxW4RUQkeClsl6PUVOu2cYKdVYj4t5gYgxeeMagWA+v+gifGmBQUKHCLiEhwUtguJ/n5JumbreVGjeytRcTfNWhg8N+nDCIjYOEiGDdBYVtERIKTwnY52ZIF+flQqRLUqW13NSL+7+yzDB552Gq5mvEJTJ+hwC0iIsFHYbucbNxo3TZqBGFh6tkWKY2uFxvccdvRKQF/WqjALSIiwUVhu5ys32CFhObNbC5EJMAMGghX9gHThMf/Y7L2TwVuEREJHgrb5eTv9dZts6Ya1RbxhmEY3HO3QYdkOHzYmqEkK0uBW0REgoPCdjlZv8G61ci2iPccDoP/PG7QrClkZ1tzcO/bp8AtIiKBT2G7HOzdZ7Jtm7XctIm9tYgEqqgog+f+a1CzBqSmwcP/NsnLU+AWEZHAprBdDtYfaSGpUxuio9VGIlJWtWoZPPeMQeXK8MtKeO4FE9NU4BYRkcClsF0O3GG7mVpIRE5b82YGTz5uEB4Gc+fB/z60uyIREZGyU9guB2uOzJ7QsoVGtUXKQ6cOBvfeY/19evd9k6/ma3RbREQCk8J2OfhjjXV7Zmt76xAJJlf2MRh0nbX8zHMmv6xU4BYRkcCjsH2asveYZGVZyy1b2FuLSLC5/RaDSy62rs760KMmqWkK3CIiElgUtk/TmiOj2o0a6uRIkfIWFmbwyIMGZ58FLhf8636T3bsVuEVEJHAobJ+mNWut//hbt7K5EJEgVamSwX/HGNSvB1lb4f6HTQ4dUuAWEZHAoLB9mn5dbd2eeaZGtUV8pVo1g+efNahaFdauhX8/bpKfr8AtIiL+T2H7NBw4YPL7H9Zy0nn21iIS7BrGGzzzlEFkJPy8BP77rElhoQK3iIj4N4Xt07DqV+vErbp1oX59jWyL+Fqbsw3GPGHNwT3vaxg/QRe9ERER/6awfRpSllv/ySe1t7kQkRByfieDhx6wPtxOnwEfTrK5IBERkZNQ2C4j0zRZstRaTm6vUW2RinTF5QZ3jbD+3r39rsknMzW6LSIi/klhu4z+XAcZGRAZqZFtETsM6G9w4xBr+eWxJrPnKnCLiIj/Udguo3lHLh99QReoUkUj2yJ2uPmfBgP6W8vPPm/y7fcK3CIi4l8Utsvg0CGTBd9Yy90vV9AWsYthGIy806BPLygshCf+Y/LzYgVuERHxHwrbZfDBJJM9OVC7tlpIROxmGAaj7zW4rBsUFMAj/zZZ8YsCt4iI+AeFbS8UFJi89c4Bpn5k/fnuEQYOh0a2RewWHm7w8AMGF3SG3Dx44CGT3/9Q4BYREfspbHvhu+/h1XEHyc+HrpdY/doi4h8cDoPH/22Q1B4OHoLR95n89bcCt4iI2Eth2wvntIGrr6rEvx82ePxRA8PQqLaIP6lUyeDp/xi0ORtc++He0SapaQrcIiJiH4VtL9SsafCfJ5xcfplBWJiCtog/qlzZ4Ln/GrRIhD05MOr/TDK3KHCLiIg9FLZFJOg4nQYvPmfQOAF27rQC9/btCtwiIlLxFLZFJChVq2bw8osGDepDVpYVuLOzFbhFRKRiKWyLSNCqUd3glRcNatWC9M1wz2iTvfsUuEVEpOIobItIUKtTxwrccbGwfgP8636TAwcUuEVEpGIobItI0GsYb7WUVK0Kf6yB+x8yOXxYgVtERHxPYVtEQkLTJtZJk1FRsHIVPPKYSV6eAreIiPiWwraIhIxWLa1pAStVgsVL4IkxJvn5CtwiIuI7CtsiElLanmNd+CYiAr7/AZ5/ycQ0FbhFRMQ3FLZFJOR0SLYu7R4WBrPnwDvvKWyLiIhvKGyLSEi66AKD0fdaV4L9YCJ8+rkCt4iIlD+FbREJWX17G9x0oxW4X3rF5IefFLhFRKR8KWyLSEj751Do0xtME5540uTX1QrcIiJSfhxledLq1asZN24cq1atIi8vj2bNmjF06FD69OlTqufv2rWLGTNm8Mcff/D777+TmZkJwLp168pSjohImRmGwf+Nguxsk4WLrDm4J4yDJo0Nu0sTEZEg4PXI9tKlSxk0aBDLly/niiuu4LrrriM7O5vRo0fzxhtvlGob69ev56WXXmL+/PlERERQuXJlrwsXESkvDofB448anHUmuFww+j6T7ds1wi0iIqfPML2Y8yo/P58ePXqwdetWpk2bRuvWrQFwuVwMHDiQTZs2MXv2bBISEk66nZ07d7Jp0yZatWqF0+mke/fubNq0yauR7ezs7FI/tjzFxsba9tpi0TGwVzC//zk5JsNHmqSlQ5PG8NpYg+ho/xvhDuZjECh0DOynY2A/HQPrPTgVr0a2lyxZQnp6Or179/YEbQCn08nw4cPJz89n5syZp9xOjRo1SEpKwul0evPyIiI+FRNjXWWyRg3YuAkefESXdRcRkdPjVdhOSUkBoEuXLsXWde7cuchjREQCUZ06Bi88a1ClCqz6FZ58SleZFBGRsvMqbKempgLQqFGjYutiYmKIjY0lLS2tXAoTEbFLs6YG/x1jXWXyhx/h8f8ocIuISNl4NRuJy+UCIDo6usT1TqeTrVu3nn5VpRATE0NYmD0zF5amP0d8S8fAXqHw/l/aFV5+IZd7Ru/j+x/gyTAHLz4XTWSkf/Rwh8Ix8Hc6BvbTMbCfjsGplWnqP3+Qk5Njy+vqZAD76RjYK5Te/7bnwNNjDB5+xOTb7/L457DdjHnCoFo1ewN3KB0Df6VjYD8dA/vpGPjgBEn3CY379u0rcb3L5TrhqLeISCDq1MHg2f8e7eG+5Q6TjRvVUiIiIqXjVdh2T+lXUl92Tk4O2dnZJfZzi4gEsqT2Bm+8ZlCvHmRlwbDbTaZOMykoUOgWEZGT8ypsJyUlAbBw4cJi6xYtWgRAcnJyOZQlIuJfGicYvDXBIDkJcnPhtddN7hhh8tvvCtwiInJiXoXtTp06ER8fz6xZs1i7dq3nfpfLxYQJE3A4HPTr189z/+7du9mwYQO7d+8uv4pFRGxSrZo1D/f9o622kjVr4Y4RJqPvL2TlKhMvrhEmIiIhwqsTJB0OB2PGjGHYsGEMGjSI3r1743Q6mT9/PhkZGYwaNYrGjRt7Hj958mTGjx/PiBEjGDlyZJFtPfDAA57lHTt2FLvvvvvuIy4urkw7JSLiK4Zh0Kc3dEiG9z8wmTMXliyFJUtN6teDHt3hsm5Qv55/zFoiIiL28no2ko4dOzJlyhTGjh3L3LlzycvLo1mzZtx999307du31Nv59NNPT3rfiBEjFLZFxG/VqmVw/78MBl9nMvkjkwXfQOYWeOc9k3feg1atTC671KDrJVCjuoK3iEioMswA/d7TrqlmNM2N/XQM7KX3v2QHD5p8/wPMX2Cy4hcoLLTuDwuzRsFvGGzQ5uzyCd06BvbTMbCfjoH9dAxKN/VfwM6zLSLiTypXNujRHXp0N9i92+S77+Hrb0x+/wMWL4HFS0zObWtyz10GTZpopFtEJFTYcwlGEZEgFhdn8I+rDd54LYypk6web4cDVq6Cm28z+XCSLv8uIhIqFLZFRHwovoHB/aPDmDbFoPP5kJcHb71jMur/TPbuVeAWEQl2CtsiIhWgdi2DZ54yePSho1ejvO1Ok4wMBW4RkWCmsC0iUkEMw+CKyw1eH2dQuzZs3gy332mycZMCt4hIsFLYFhGpYE2aGLw5wSCxOezJgXv+z2SzRrhFRIKSwraIiA1qVDd45UWDpk1g1264+16TrdsUuEVEgo3CtoiITapWNXj5BYOG8bB9O9z3gMn+/QrcIiLBRGFbRMRGcXEGL79oUL06bNwE/35C0wKKiAQThW0REZvVrmXw7FMGlSrB0hQYO15hW0QkWChsi4j4gZYtDf79sHVlyZmfwYyZCtwiIsFAYVtExE9cdKHBHbdZgXvseJPFSxS4RUQCncK2iIgfGTQQevWEwkKrf/vv9QrcIiKBTGFbRMSPGIbB6HsM2p0LBw/C6PtNtm5V4BYRCVQK2yIifiYiwuCpJw0aJ8CuXVbg3rtXgVtEJBApbIuI+KHoaIMXnjOoWQNS06zA7XIpcIuIBBqFbRERP1W7lhW4q1aFNWvhnn+Z7NunwC0iEkgUtkVE/FjTJgavvmgQUxXWroVR/2eyfXuh3WWJiEgpKWyLiPi55s0NXn3ZoFoMrPsLrrluD6t/0wi3iEggUNgWEQkAzZoavPGaQZPGsHOnychRJh9MNMnNVegWEfFnCtsiIgGiQQMrcF9xeSQFBfD2uyY33myyNMXENBW6RUT8kcK2iEgAiYoyePE5J/9+xCA2FtI3w//dZzLsNpNvvjPJz1foFhHxJwrbIiIBxjAMLu9mMOVDgwH9oVIlq5f7sSdM/jHA5M13CsnKUugWEfEHCtsiIgEqOtrgrhFhfDLN4J9DITYWdu2GiZNgwCCTe/9VyA8/abRbRMRODrsLEBGR01OtmsHN/zQYcr3Jop/h8y9Nli2HlGWQssykehz07GnSt5dB3bqG3eWKiIQUhW0RkSAREWFw8UVw8UUGmVtMvpxlMnvu0dHuSZNNLuhiMuyfBk2aKHSLiFQEtZGIiASh+vUMbr81jJnTDcY8YZDUHkwTfvwJht5s8vh/CtmcofYSERFfU9gWEQli1mi3wcsvhDHpA4NLLrZC94Jv4PohJv99rpCtWxW6RUR8RWFbRCREJDQy+M/jYbz3tsH5HaGgEGbPgYHXm7z0SiE7dyp0i4iUN4VtEZEQk9jc4LlnwnjjNYPz2kF+Psz8zJrB5LXXC9mzR6FbRKS8KGyLiISos840ePWlMMa+bHDWmZCbC1OnwTUDTV54qZC//1boFhE5XZqNREQkxLU71+D18bBkqXUJ+L/+hs++gM++MGna1OTCLtCls0GzphAerllMRES8obAtIiIYhkGnjtCxA/yy0pqr+8efYMMG6+f9D0yiouDM1iZNGkP9+gbxDaBOHYipCk4nhIUpiIuIHE9hW0REPAzD6uM+r51BTo51kZwfF5r8shIOHIBly60fMI97HkRHm0REWMvGkfsKCqGgAAoLobDAWi4otGZEqVzZCulxsVC/PsQ3MGiRCK1aWhfqEREJBgrbIiJSopgYg549oGcPg4ICk40b4Y+1sDnDJDMTMjJh2zY4eNAKz3v3erf9vDzrOVu2wO9/wLEBvn49k7PPhg7JBknnKXyLSOBS2BYRkVMKDzdo3hyaNwdr3PqovDyTffsgZy/k51mR2TStn/BwCA+DsHAIC7OWw8OtTRw8CPv2wc5dkJkJqWkma9dC+mbI3GL9fDXPxDCgZUuTjslW+G7VUr3jIhI4FLZFROS0REQYxMVBXNzpbskK0Hv3WaF7xUqTpSlWz/jatdbP+x+YREdDh2STTh0MOiRr1FtE/JvCtoiI+JWq0VaI7pBsMPw22LHDJGU5LE0xSVlmjYYv+AYWfGONerduZdKpo3WCZ/NmOlFTRPyLwraIiPi1mjUNevWAXj0M8vNN1qyFxUtMFi+B9RvgjzXwxxqTd96D6nHQsaM16p3UHqpUUfAWEXspbIuISMBwOAzanA1tzja47RbYvt1kSYoVvpcvh127rUvQz55jEh4O57QxOb+TwUUXQN26Ct4iUvEUtkVEJGDVqmXQtzf07W2Qm2vy6+ojo95LYfNma87wX1aajJ8ALRJNLr7I4KILoWG8greIVAyFbRERCQqRkVbrSFJ7g7tGQEaGFbp//MkK4ev+gnV/mbz5NjRtYnLRhQYXXwSNE6z5xUVEfEFhW0REglKDBgbXNIBr/mGQnW3y0yL4/geTFb/Aho2wYaPJe/+DhvFw0UUml1xk0LyZgreIlC+FbRERCXqxsUfbTfbuNVm4CL7/0WTZcmte74mTYOIkk7p14eILrXaT1q0UvEXk9Clsi4hISKla9eiVMffvN1m0GH740WTJUsjKgqnTYOo0k1o1oe05JmeeaXDWmdC0iXWCpoiINxS2RUQkZFWpYnB5N7i8m8HBg1bg/v5Hk58Xw/YdMH8BzF9gXUY+IgLiG5gkJEDLFgeIiTGJi7Uu5hMTA5UiweGAyEjrsZrvW0RAYVtERASAypUNLrkYLrnY4PBh66TK3/+A3/8wWbMGXPth4ybr59vvDp5ye+FhJuEO6/L0Dof1E1XZCubVqlk/9eoaNKgP8fEQ3wCiohTQRYKNwraIiMhxKlUySE6C5CQAg8JCk23bIDXN+tm6NZLMLYfZvRt2Z8PeHMjNK7qNgkIoyC16X3Y2ZG459h6zyPq6dUwaN4YmjaFJY4PGjaFRQ2umFREJTArbIiIipxAWZlC3LtStC506Qmysk+zsounaNE3y8yEvD3JzIT/f+ikoOLJcAPv3Q04O7NljXYBnyxaTzRmwOcO6L2ur9fPzYnAH8fAwaNDAHcKNI0Ec6tVTD7lIIChT2F69ejXjxo1j1apV5OXl0axZM4YOHUqfPn1KvY3CwkKmTJnCtGnTSEtLIyoqig4dOnDPPfeQkJBQlrJERERsYxgGERFWv3ZUVKmf5VnKyTHZlOpuVTHZtMmaotDlgrR06+f7H46OhEdGQKNGJo0ToGZNa8aV2GoQGwtVq8IZZ1h95JUqWT8OB5hm8Z/i+2Hdhodbz9OMLCKnx+uwvXTpUm6++WYiIiLo1asX0dHRzJ8/n9GjR5OZmcntt99equ089thjTJ8+nWbNmnH99deza9cu5syZw6JFi/joo49o1qyZ1zsjIiISqGJiDNqeA23PAXcIN02TXbuO9oq7Q/imVDh0CP5eb/1YSkjOpyk8HKpUMakSBVFVoFoM1KgBNWtAzRoGNWpCrZpW2I+L1UmhIiUxTLOkz7Uly8/Pp0ePHmzdupVp06bRunVrAFwuFwMHDmTTpk3Mnj37lCPTS5YsYejQobRv357333+fyMhIABYvXsw///lP2rdvz6RJk066jezs7NKWXa5iY2Nte22x6BjYS++//XQM7Gf3MSgsNMnaCps2WT3ku3ebZO+xesKz91itKocPWz+5uafaWvkIDz8SwmtCrVpWCI+NNYiKsk4MjYqyfsLDiz+3oMBqv8nLh/zjb49ZzsuD/HyTvDwIDz+D/fsPFXmstf7obf5x95lYs8VUijw6a0xkJFQ+w/o2oGpV48gtxBy5rVYNYqtZffy+cviwyYED1km4ZqFV07E/ERH++Q2D3X8P/EFsbOwpH+PVyPaSJUtIT0/n6quv9gRtAKfTyfDhw7nnnnuYOXMm995770m38/HHHwMwatQoT9AG6NSpE126dOGnn35i06ZNNG7c2JvyREREQkJYmEH9elC/HnTpDMe2oxzPNE1yc61AC1abSFjY0eVjM9yxw2/5+VaP+f4DVivL/v1WX/mOnbBzp8mOndb0iDt3WP3nBQWwdZv1c8wWy2eHS3TIB9s8cb2VK5ueWWRiqx2dUaZajMEZZ1htOuHhEGbA4dyjH3YOHjSt9/GY99IdrA8cuS8v74QvCxz9hiGqMlSpcvSDS5UqWN86eP5sfbhxfxNR+Qw8rU2OcHBEQITjyJ+PubVq912YN02TgoKjH6rc5zDkH9lvR4TVFuWuNTzcPz9clJVXYTslJQWALl26FFvXuXPnIo85maVLlxIVFUW7du2KrXOH7WXLlilsi4iInCbDMKhUyfvnRUZaAa5myVst8qf8fJNdu2HHDti23brdvt0kJ8cKlgcOHrk9YM3ScrQ269YRXnIAPPa+CAeEO46GRafzDAoKDuFwGKd8rPt+sGaNyc21fvKO3B44CHv3muTshX17IWcv7N139GTW/Hw4eND6yco6/r0ovw8UUVHWByF3fW4FBbB3r/VzcmWvJSzMLPL+OY6/Pbbn/8hLhYXtIT+/EPNIjfnHBOm8PCg4ZtkbhgEREaYnfEdGHA3kDgdERB5djoyEC7oY9O3tv+Hcq7CdmpoKQKNGjYqti4mJITY2lrS0tJNu48CBA+zYsYPExETCS/guyd2C4n4tERER8W8Oh0HtWlC7Fpx1pvte34af2NgqZGeXZ49MyfWaponLBXuOBO89e6xWHWvZZM+eI7PPHJl1prDw6EmplSpB5crgrGLNoV6lypHlI7dVoqCK07qtXLno6LJpWu0y1ui49YFg//6jH1rcI+UHDsD+AyYH9ruXj94ePFi0DefY22M/9IBV9/Eh/9QKvHqHj2UYVogHq6Zjv1UxTe9q2bTJDJ6w7XK5AIiOji5xvdPpZOvWrSfdxr59+zyPPdE2jn2tE4mJiSHM/T1YBStNf474lo6BvfT+20/HwH46BvarqGMQFwcNK+SVKk5h4dGpKvPyzKK3+cf9+ch9Bkdbj6wfK+AaxpERZwc4Iowj30wYx3xLcWSmHodRYttKkWkz3a+Za7U/5eWZ5B6pxQrgZrHHtWnjIDbWf2ez9t/KTiEnJ8eW19XJAPbTMbCX3n/76RjYT8fAfjoG5csdgitXLv1zvDkGphcj545wcFT2rha7fhVK84HPq6Fh96ize3T6eC6X64Sj3m7u9ScauXbff6KRbxERERGRQOFV2Hb3U5fUl52Tk0N2dnaJ/dzHioqKombNmmRkZFBQULzXx92rrQvbiIiIiEig8ypsJyUlAbBw4cJi6xYtWgRAcnLyKbeTnJzMgQMH+OWXX4qtc2/b/VoiIiIiIoHKq7DdqVMn4uPjmTVrFmvXrvXc73K5mDBhAg6Hg379+nnu3717Nxs2bGD37t1FtjNgwAAAXnnlFXKPad5ZvHgxCxcuJCkpSdP+iYiIiEjA8+oESYfDwZgxYxg2bBiDBg2id+/eOJ1O5s+fT0ZGBqNGjSoSkidPnsz48eMZMWIEI0eO9NzfsWNHrrnmGj7++GP69evHRRdd5Llcu9Pp5PHHHy+3HRQRERERsYvXs5F07NiRKVOmMHbsWObOnUteXh7NmjXj7rvvpm/fvqXezpNPPkmLFi2YNm0aEydOJCoqiksuuYR77rlHo9oiIiIiEhQM0zR9eS1Vn7Fruh9NNWQ/HQN76f23n46B/XQM7KdjYD8dAx9M/SciIiIiIqWnsC0iIiIi4iMK2yIiIiIiPqKwLSIiIiLiIwrbIiIiIiI+orAtIiIiIuIjCtsiIiIiIj6isC0iIiIi4iMK2yIiIiIiPqKwLSIiIiLiIwF7uXYREREREX+nkW0RERERER9R2BYRERER8RGFbRERERERH1HYFhERERHxEYVtEREREREfUdgWEREREfERh90FBILVq1czbtw4Vq1aRV5eHs2aNWPo0KH06dPH7tJCwueff86KFSv4/fff+euvv8jLy+O///0vV199td2lhYRt27Yxd+5cfvzxRzZu3MjOnTuJiYmhXbt2DBs2jHPOOcfuEoPe3r17GTt2LL/99hsZGRnk5OQQGxtL48aNGTx4MJdffjmGYdhdZkh5++23eeGFFwCYNm0abdu2tbegENC1a1cyMzNLXHfttdfy5JNPVnBFoevrr79mypQprFmzhoMHD1KjRg3atm3Lv/71L+rWrWt3eX5HYfsUli5dys0330xERAS9evUiOjqa+fPnM3r0aDIzM7n99tvtLjHovfrqq2RmZhIbG0utWrVO+I+t+MbEiRN5++23adiwIeeffz7Vq1cnLS2NBQsWsGDBAl588UV69uxpd5lBLTs7m08++YRzzjmHSy+9lGrVqrFr1y6+++477rrrLgYMGMB//vMfu8sMGRs2bGDs2LFERUVx4MABu8sJKdHR0QwdOrTY/WeddZYN1YQe0zR57LHHmDZtGg0bNqRnz55UqVKF7du3s2zZMjIzMxW2S6CL2pxEfn4+PXr0YOvWrUybNo3WrVsD4HK5GDhwIJs2bWL27NkkJCTYW2iQ+/nnn2nUqBH169fnrbfe4sUXX9TIdgWaP38+cXFxtG/fvsj9y5cv58Ybb6RKlSr89NNPREZG2lRh8CsoKMA0TRyOouMjLpeLa6+9lvXr1zNr1iyaN29uU4Who6CggGuvvRbDMEhISOCLL77QyHYF6dq1KwDffvutzZWErg8//JCnnnqKwYMH8/DDDxMeHl5kfX5+frF/p0Q92ye1ZMkS0tPT6d27tydoAzidToYPH05+fj4zZ860scLQcP7551O/fn27ywhZl19+ebGgDdC+fXs6dOjAnj17WLdunQ2VhY7w8PAS/wNzOp106dIFgLS0tIouKyS9/fbb/Pnnnzz99NPFgoZIMDt06BCvvfYa8fHxPPTQQyX+/itol0zvykmkpKQAeP4zO1bnzp2LPEYkFLn/YdU/sPY4fPgwS5YswTAMmjVrZnc5Qe+vv/5i/Pjx3HHHHfoWwSa5ubl8+umnbNu2japVq9KuXTtatmxpd1khYdGiRezZs4d+/fpRWFjI/PnzSU1NJTo6mvPPP59GjRrZXaLf0v+QJ5GamgpQ4i9QTEwMsbGxGk2SkLVlyxZ+/vlnatasSWJiot3lhIS9e/fywQcfUFhYyK5du/jxxx/JyspixIgRamfzsfz8fB544AGaNm3Krbfeanc5IWvHjh088MADRe674IILeO6554iLi7OpqtDw+++/A9Y3bX379mXTpk2edWFhYdx4443cf//9dpXn1xS2T8LlcgHWCRklcTqdbN26tSJLEvELeXl53HfffeTm5jJ69Gh9nV5B9u7dy/jx4z1/joiI4L777uOmm26ysarQ8MYbb7Bu3TqmT59ORESE3eWEpKuvvprk5GSaNWtGZGQkGzZsYPz48fz4448MHz6cqVOnalYeH9q1axcA77//Pq1bt+bjjz+madOmrF27lkcffZT33nuP+Ph4Bg0aZHOl/kc92yLilcLCQh566CGWLVvGgAEDuOqqq+wuKWQ0aNCAdevWsWbNGr755hvuuusuXn75ZUaOHEl+fr7d5QWtP//8kzfeeIObbrqJM8880+5yQtaIESNITk4mLi4Op9PJOeecw5tvvsl5553HypUr+eGHH+wuMai559OIiIjgtddeo02bNlSpUoX27dszduxYwsLCeP/9922u0j8pbJ+E0+kEYN++fSWud7lcJxz1FglGpmnyyCOP8MUXX9C3b1+eeOIJu0sKSeHh4TRo0IBbb72VUaNG8fXXXzN9+nS7ywpa999/P/Hx8YwcOdLuUuQ4YWFhnpmpfvnlF5urCW7uTHTWWWdRu3btIuuaN29OfHw86enp7N27147y/JrC9km4eyBL6svOyckhOztbJwRIyHCPaH/yySf07t2bZ555hrAw/RNiN/cJ3DpZ23f+/PNPNm7cyNlnn02LFi08P59++ilgXVClRYsWLFiwwOZKQ1NsbCwABw8etLmS4NakSRPgxK217vsPHTpUYTUFCvVsn0RSUhJvvvkmCxcupFevXkXWLVq0CIDk5GQ7ShOpUIWFhTz88MPMnDmTnj178txzz6lP209s27YNQMfDh/r371/i/cuXLyc1NZWuXbsSFxenKUptsnr1agC9/z7WoUMHADZu3FhsXV5eHunp6URFRelE1RIobJ9Ep06diI+PZ9asWQwZMoRWrVoBVvvIhAkTcDgc9OvXz+YqRXzr2KDdvXt3nn/+eQW7CrZ27VoaNGhQbERpz549vPzyywBceOGFdpQWEp566qkS73/ggQdITU3ltttu00VtfGz9+vXUqlWLqlWrFrl/+fLlvP/++0RGRnL55ZfbVF1oaNiwIV26dGHhwoV8/PHHXHPNNZ51b731Fnv37qVv376aCrYEekdOwuFwMGbMGIYNG8agQYPo3bs3TqeT+fPnk5GRwahRo2jcuLHdZQa9jz/+mBUrVgDWPLfu+9xfm3fr1o1u3brZVl+we+2115g5cyZRUVEkJCTw+uuvF3tMt27dPB9GpfzNnDmTGTNm0KFDB+rVq0flypXZsmUL33//PQcOHOCKK66gT58+dpcp4jNz587lnXfeoVOnTtSvX5/IyEj++usvFi1aRFhYGE888QT16tWzu8yg99hjjzFw4EAeeeQRFixYQJMmTVizZg1Lliyhfv363HfffXaX6JcUtk+hY8eOTJkyhbFjxzJ37lzy8vJo1qwZd999N3379rW7vJCwYsUKT2+k2y+//OI5GaZ+/foK2z6UmZkJwIEDB3jjjTdKfEz9+vUVtn3oiiuuwOVysWrVKpYtW8ahQ4eIiYnhvPPO46qrrqJXr16a8kyCWocOHdiwYQNr1qwhJSWF3NxcqlevTs+ePbnxxhtp06aN3SWGhIYNG/LJJ58wduxYfvrpJxYtWkSNGjUYPHgwd955J9WrV7e7RL9kmO65XEREREREpFxpKgERERERER9R2BYRERER8RGFbRERERERH1HYFhERERHxEYVtEREREREfUdgWEREREfERhW0RERERER9R2BYRERER8RGFbRERERERH1HYFhERERHxEYVtEREREREfUdgWEREREfGR/werSTwPdptLHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.style.use(\"arviz-darkgrid\")\n",
    "ax = az.plot_dist(draws)\n",
    "ax.set_title(\n",
    "    f\"Random draws from $Gamma({round(constrained_priors['alpha'], 1)}, {round(constrained_priors['beta'], 1)})$\"\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4005128a",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "And boom! Indeed, our constraint seems to be respected. To be completely honest, I'm _sure_ it is -- the maths are telling us: PyMC runs some tests in the background, so if any issue appears, it would have warned us that our constraint is not respected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21d9974c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And on that note, PyMCheers, and best Bayesian wishes, wherever you are 🖖\n",
    "\n",
    "![MicDrop](https://media.giphy.com/media/3o7qDSOvfaCO9b3MlO/giphy.gif)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Diaporama",
  "kernelspec": {
   "display_name": "pymc-dev-py39",
   "language": "python",
   "name": "pymc-dev-py39"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
No results found