Simulate_random_variable_with_Inverse_transform_method.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:white;padding:20px;\">\n",
    "\n",
    "<p style=\" color : red; font-size : 40px; font-family : consulas; text-align : left;\">\n",
    "\n",
    "<b>\n",
    "Simulate Random Variables, With Using The Inverse Transform Method\n",
    "</b>\n",
    "</p>\n",
    "\n",
    "<br><hr><hr><br><br><br><br>\n",
    "\n",
    "\n",
    "<p style=\" color : darkblue; font-size : 25px; font-family : courier new; text-align : left;\">\n",
    "<br><br>\n",
    "<b>\n",
    "Example I \n",
    "</b>\n",
    "</p>\n",
    "\n",
    "\n",
    "<hr><br><br>\n",
    "<p style=\" color : darkblue; font-size : 20px; font-family : courier new; text-align : left;\">\n",
    "<b>\n",
    "Using the inverse transform method, simulate a random variable with the following density function.\n",
    "</b>\n",
    "</p>\n",
    "<br>\n",
    "\n",
    "$$\n",
    " f(x) = \\frac{1}{2}x + \\frac{1}{2}, \\quad 0 < x < 1, \n",
    "$$\n",
    "\n",
    "#### Answer:\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& F_X(x) = \\int_0^x f(x)~dx = \\int_0^x \\frac{1}{2}x + \\frac{1}{2} ~dx = \\\\\n",
    "& \\left. \\frac{x^2}{2} + \\frac{x}{2} \\right|_0^x = \\frac{x^2}{2} + \\frac{x}{2}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "based on  <a  href=\"https://math.stackexchange.com/questions/868400/showing-that-y-has-a-uniform-distribution-if-y-fx-where-f-is-the-cdf-of-contin\"> this link </a> $~~F_X(x) \\sim \\text{Unif}(0, ~1) \\implies $\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& u = \\frac{x^2}{2} + \\frac{x}{2}  \\quad \\overset{{\\tiny \\text{Using} ~\\Delta~\\text{Method}}}{\\implies} \\\\\n",
    "& x = \\frac{-\\frac{1}{2} \\pm \\sqrt{\\frac{1}{4} + 2u}}{1} \\implies \\\\\n",
    "& x_1 = -\\frac{1}{2} +  \\sqrt{\\frac{1}{4} + 2u}, \\\\\n",
    "& x_2 = -\\frac{1}{2} -  \\sqrt{\\frac{1}{4} + 2u} \\quad \\overset{{\\tiny 0 < x < 1}}{\\to} \\quad \\text{Unacceptable} \\implies \\\\\n",
    "& x = -\\frac{1}{2} +  \\sqrt{\\frac{1}{4} + 2u} \\implies \\\\\n",
    "& U = [0.1915,~ 0.6221,~ 0.4377,~ 0.7854,~ 0.78] \\implies \\\\\n",
    "& X = [0.2956,~ 0.7224,~ 0.5609,~ 0.8493,~ 0.8453]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## simulate in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy.stats import ContinuousRV, P, E, variance, cdf\n",
    "from sympy import Interval, Symbol, log, exp, sqrt, Rational"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "## define X as a Random Variable\n",
    "\n",
    "x = Symbol('x')\n",
    "pDf =  x + Rational(1, 2)\n",
    "X = ContinuousRV(x, density = pDf, set = Interval.open(0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{cases} \\begin{cases} 0 & \\text{for}\\: x < 0 \\\\\\frac{\\left(\\min\\left(1, x\\right) + 1\\right) \\min\\left(1, x\\right)}{2} & \\text{otherwise} \\end{cases} & \\text{for}\\: x \\geq 0 \\\\0 & \\text{otherwise} \\end{cases}$"
      ],
      "text/plain": [
       "Piecewise((Piecewise((0, x < 0), ((Min(1, x) + 1)*Min(1, x)/2, True)), x >= 0), (0, True))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get cumulative function (CDF) of X \n",
    "\n",
    "cdf(X)(x).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{7}{12}$"
      ],
      "text/plain": [
       "7/12"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get expectation of X \n",
    "\n",
    "E(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{11}{144}$"
      ],
      "text/plain": [
       "11/144"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get Variance of X\n",
    "\n",
    "variance(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1915 0.6221 0.4377 0.7854 0.78  ]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[0.2956, 0.7224, 0.5609, 0.8493, 0.8453]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get above result \n",
    "\n",
    "U = stats.uniform.rvs(0, 1, size = 5, random_state = 1234)\n",
    "print(U.round(4))\n",
    "def f(u):\n",
    "    a = -1/2 + sqrt(1/4 + 2*u)\n",
    "    return a\n",
    "f = np.vectorize(f)\n",
    "xx = f(U)\n",
    "res = list(map(lambda x: x.round(4), xx))\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Mean of Simulation: 0.584672637901349\n",
      "Mean of Theory: 0.5833333333333334\n",
      "Variance of Simulation: 0.0758054328780650\n",
      "Variance of Theory: 0.0763888888888889\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## simulate with transform method and compare variance and expectation of simulation with Text \n",
    "n = 10**5\n",
    "U1 = stats.uniform.rvs(0, 1, size = n, random_state = n)\n",
    "XX = f(U1)\n",
    "sim_var = XX.var()\n",
    "sim_mean = XX.mean()\n",
    "print(\"\"\"\n",
    "Mean of Simulation: {}\n",
    "Mean of Theory: {}\n",
    "Variance of Simulation: {}\n",
    "Variance of Theory: {}\n",
    "\"\"\".format(sim_mean, 7/12, sim_var, 11/144))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAK9CAYAAAA37eRrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABo4UlEQVR4nO3deVhV5eL28XszqwhEKoqiFOVUTjmlpWZiTj/NY4NZIk7lXIo5l6SlKGp6Ss1yHsrxlKeUnDXNLM3h6MnUNGcSJRUUkWmv9w9fOZGgorAXC76f69pXsXjWXvdmQdwtnv0sm2EYhgAAAAALcjI7AAAAAHCvKLMAAACwLMosAAAALIsyCwAAAMuizAIAAMCyKLMAAACwLMosAAAALIsyCwAAAMuizAIAAMCyKLOAxQUGBqpz585mxygwtmzZIpvNpi1bttz12BUrVtzz8Xbt2qX69eurSJEistls2rdvn9577z3ZbLZ7fs6sZOe15aYTJ07IZrNp3rx5uX6s1NRUDR48WAEBAXJyclLbtm1z/Zj3Yt68ebLZbDpx4oTZUYA8hzIL5CE3f2H9/PPPmX7+mWee0eOPP37fx4mKitJ7771338+DG7744gtNmTIlx583JSVFL730ki5evKjJkydr4cKFKleuXKZjx44dq5UrV+Z4hrwit75n58yZowkTJujFF1/U/PnzNWDAgBw/Rnbk9/MI5AoDQJ4xd+5cQ5Kxa9euTD/fqFEj47HHHsuw7fr160ZycnK2jtOnTx+DH/97k5aWZiQmJhppaWnp21q1amWUK1fulrGbN282JBnLly+/p2P9+uuvhiRj5syZGbanpKQYiYmJGbYVKVLECA0Nvafj3HQz7+bNm+/ree6X3W43EhMTjdTU1PRtufU92759e6N06dI5/rz3KqvzmJqaaiQmJhp2u93xoYA8zsXMIg3g/rm7u5sdIdsSEhJUpEgRs2PcEycnJ3l4eDjkWOfPn5ck+fj4ZNju4uIiF5f8+59vm83m0K/x37++eZGzs7OcnZ3NjgHkSUwzACzu73NmU1JSNGrUKD366KPy8PDQgw8+qKefflrr16+XJHXu3FnTpk2TdKM03HzclJCQoIEDByogIEDu7u6qUKGCJk6cKMMwMhw3MTFRb775pooVK6aiRYuqTZs2Onv2rGw2W4Y/B9+c33nw4EG9+uqreuCBB/T0009Lkvbv36/OnTvr4YcfloeHh0qWLKmuXbvqzz//zHCsm89x5MgRdezYUd7e3ipevLjeffddGYah06dP6/nnn5eXl5dKliypSZMm3fHr1q5dOz3xxBMZtrVu3Vo2m01ff/11+raffvpJNptN3377raRb55U+88wzWr16tU6ePJn+tQwMDMzwvHa7XWPGjFGZMmXk4eGhJk2a6OjRo7fN17lzZzVq1EiS9NJLL8lms+mZZ57J8PW4yWazKSEhQfPnz0/PcKd51GfOnFHbtm1VpEgRlShRQgMGDFBSUlKmY3/66Sc1b95c3t7eKly4sBo1aqTt27dnGHMz09GjR9W5c2f5+PjI29tbXbp00bVr1zKMXb9+vZ5++mn5+PjI09NTFSpU0PDhw9M///c5s1l9zxqGocDAQD3//PO3ZL5+/bq8vb3Vo0ePTF/TzWNs3rxZv/zyS/pzbtmyJcu5w5nN5e3cubM8PT119uxZtW3bVp6enipevLjefvttpaWlZdjfbrfrn//8p6pUqSIPDw8VL15czZs3T59WdLvzmNWc2enTp+uxxx6Tu7u7/P391adPH12+fDnDmJvTkw4ePKjGjRurcOHCKl26tCIjIzP92gBWk3//1x6wsLi4OMXGxt6yPSUl5Y77vvfee4qIiFD37t1Vp04dxcfH6+eff9aePXvUtGlT9ejRQ9HR0Vq/fr0WLlyYYV/DMNSmTRtt3rxZ3bp1U/Xq1bV27VoNGjRIZ8+e1eTJk9PHdu7cWcuWLVNISIiefPJJfffdd2rVqlWWuV566SU9+uijGjt2bHoxXr9+vX7//Xd16dJFJUuW1C+//KLPPvtMv/zyi3788cdb3uTUvn17VapUSePGjdPq1av1wQcfyNfXV59++qmeffZZjR8/Xp9//rnefvtt1a5dWw0bNswyT4MGDfTvf/9b8fHx8vLykmEY2r59u5ycnLRt2za1adNGkrRt2zY5OTnpqaeeyvR5RowYobi4OJ05cyb96+Pp6ZlhzLhx4+Tk5KS3335bcXFxioyM1Guvvaaffvopy3w9evRQ6dKlNXbsWL355puqXbu2/Pz8Mh27cOHC9PP9xhtvSJKCgoKyfO7ExEQ1adJEp06d0ptvvil/f38tXLhQmzZtumXspk2b1KJFC9WsWVPh4eFycnLS3Llz9eyzz2rbtm2qU6dOhvEvv/yyHnroIUVERGjPnj2aNWuWSpQoofHjx0uSfvnlF/3f//2fqlatqtGjR8vd3V1Hjx69pRz//WuR2feszWZTx44dFRkZqYsXL8rX1zf9c998843i4+PVsWPHTJ+zePHiWrhwocaMGaOrV68qIiJCklSpUiX9+uuvWWbJTFpampo1a6a6detq4sSJ2rBhgyZNmqSgoCD16tUrfVy3bt00b948tWjRQt27d1dqaqq2bdumH3/8UbVq1cr2eXzvvfc0atQoBQcHq1evXjp8+LA++eQT7dq1S9u3b5erq2v62EuXLql58+Zq166dXn75Za1YsUJDhgxRlSpV1KJFi2y9XiDPMXOOA4CMbs6Zvd3j73Nmy5Url2GOXbVq1YxWrVrd9jhZzT9cuXKlIcn44IMPMmx/8cUXDZvNZhw9etQwDMPYvXu3Icno379/hnGdO3c2JBnh4eHp28LDww1JRocOHW453rVr127ZtnjxYkOSsXXr1lue44033kjflpqaapQpU8aw2WzGuHHj0rdfunTJKFSo0B3nj+7atcuQZERFRRmGYRj79+83JBkvvfSSUbdu3fRxbdq0MWrUqJH+cWbzSu80Z7ZSpUpGUlJS+vZ//vOfhiTjwIEDt82Y1Zzbm1+Pv8rOnNkpU6YYkoxly5alb0tISDAeeeSRDK/Nbrcbjz76qNGsWbMMczWvXbtmPPTQQ0bTpk1vydS1a9cMx/rHP/5hPPjgg+kfT5482ZBkXLhwIct8x48fNyQZc+fOTd+W1ffs4cOHDUnGJ598kmF7mzZtjMDAwDvOMc1sHnpWc4czyxUaGmpIMkaPHp1hbI0aNYyaNWumf7xp0yZDkvHmm2/ekuGvGbM6jzf/23D8+HHDMAzj/Pnzhpubm/Hcc89lmL89depUQ5IxZ86cDK9RkrFgwYL0bUlJSUbJkiWNF1544dYvCmAxTDMA8qBp06Zp/fr1tzyqVq16x319fHz0yy+/6Lfffsv2caOiouTs7Kw333wzw/aBAwfKMIz0P7WvWbNGktS7d+8M4/r165flc/fs2fOWbYUKFUr/9+vXrys2NlZPPvmkJGnPnj23jO/evXv6vzs7O6tWrVoyDEPdunVL3+7j46MKFSro999/zzKLJNWoUUOenp7aunWrpBtXYMuUKaNOnTppz549unbtmgzD0Pfff68GDRrc9rnupEuXLnJzc0v/+Obz3SljbomKilKpUqX04osvpm8rXLhw+tXAm/bt26fffvtNr776qv7880/FxsYqNjZWCQkJatKkibZu3Sq73Z5hn7+f5wYNGujPP/9UfHy8pP/N//33v/99y773onz58qpbt64+//zz9G0XL17Ut99+q9deey1XljDLTGav+6/n91//+pdsNpvCw8Nv2fdeMm7YsEHJycnq37+/nJz+96v89ddfl5eXl1avXp1hvKenZ4ar1G5ubqpTp45p34NATqLMAnlQnTp1FBwcfMvjgQceuOO+o0eP1uXLl1W+fHlVqVJFgwYN0v79++/quCdPnpS/v7+KFi2aYXulSpXSP3/zn05OTnrooYcyjHvkkUeyfO6/j5VulI633npLfn5+KlSokIoXL54+Li4u7pbxZcuWzfCxt7e3PDw8VKxYsVu2X7p0Kcss0o0yXK9ePW3btk3SjTLboEEDPf3000pLS9OPP/6ogwcP6uLFi/ddZv+e++Z5vFPG3HLy5Ek98sgjt5SoChUqZPj45v8QhYaGqnjx4hkes2bNUlJS0i3n6U6vtX379nrqqafUvXt3+fn56ZVXXtGyZcvuq9h26tRJ27dvT//+XL58uVJSUhQSEnLPz5kdN+e//tUDDzyQ4fweO3ZM/v7+GaZC3I+br/Xv58zNzU0PP/xw+udvKlOmzC3n++8ZAauizAL5TMOGDXXs2DHNmTNHjz/+uGbNmqUnnnhCs2bNMjXXX6/C3vTyyy9r5syZ6tmzp7788kutW7cu/apvZuUms3dzZ/UOb+Nvb1jLzNNPP61du3bp+vXr6WXWx8dHjz/+uLZt25ZedO+3zN5PRjPdPAcTJkzI9C8F69evv2V+8J1ea6FChbR161Zt2LBBISEh2r9/v9q3b6+mTZve8oapu/XKK6/I1dU1/ersokWLVKtWrVuK3t3K6kppVvmssMqAVb8HgbtBmQXyIV9fX3Xp0kWLFy/W6dOnVbVq1QwrDGT1y7pcuXKKjo7WlStXMmw/dOhQ+udv/tNut+v48eMZxt3pHfp/denSJW3cuFFDhw7VqFGj9I9//ENNmzbVww8/fNfPcb8aNGig5ORkLV68WGfPnk0vrQ0bNkwvs+XLl8/yjVc3OepP2TmVoVy5cjp27NgtRebw4cMZPr755iMvL69M/1IQHByc4U1Gd8vJyUlNmjTRhx9+qIMHD2rMmDHatGmTNm/enOU+t3t9vr6+atWqlT7//HOdPHlS27dvv6+rsjevJv99VYC/X+3MjqCgIEVHR+vixYu3HXe35/Hmz+Lfz1lycrKOHz+e5c01gPyIMgvkM39f1srT01OPPPJIhmWXbq7x+vdf1i1btlRaWpqmTp2aYfvkyZNls9nS3/XcrFkzSTeWBfqrjz/++K5z3rxS9PdClRt30spK3bp15erqqvHjx8vX11ePPfaYpBsl98cff9R33313V1dlixQpkum0CEcqUqTILeczKy1btlR0dHSG2+xeu3ZNn332WYZxNWvWVFBQkCZOnKirV6/e8jwXLlzIds7Mylz16tUlKculwaSsv2dvCgkJ0cGDBzVo0CA5OzvrlVdeyXa2m8qVKydnZ+f0+dQ3/f37PTteeOEFGYahUaNG3fK5v/4M3O15DA4Olpubmz766KMM+8+ePVtxcXG3XVkEyG9YmgvIZypXrqxnnnlGNWvWlK+vr37++WetWLFCffv2TR9Ts2ZNSdKbb76pZs2apf/yb926tRo3bqwRI0boxIkTqlatmtatW6d///vf6t+/f/qVupo1a+qFF17QlClT9Oeff6YvzXXkyBFJd3d1ycvLSw0bNlRkZKRSUlJUunRprVu37parvbmpcOHCqlmzpn788cf0NWalG1dmExISlJCQcFdltmbNmlq6dKnCwsJUu3ZteXp6qnXr1rkd/5YMGzZs0Icffih/f3899NBDqlu3bqZjX3/9dU2dOlWdOnXS7t27VapUKS1cuFCFCxfOMM7JyUmzZs1SixYt9Nhjj6lLly4qXbq0zp49q82bN8vLy0vffPNNtnKOHj1aW7duVatWrVSuXDmdP39e06dPV5kyZdLXH87q9Um3fs/e1KpVKz344INavny5WrRooRIlSmQr1195e3vrpZde0scffyybzaagoCCtWrUq/SYW96Jx48YKCQnRRx99pN9++03NmzeX3W7Xtm3b1Lhx4/Sfz7s9j8WLF9ewYcM0atQoNW/eXG3atNHhw4c1ffp01a5dO8slyYB8yZxFFABk5l5uZ/v3pbk++OADo06dOoaPj49RqFAho2LFisaYMWMy3PI2NTXV6Nevn1G8eHHDZrNlWPLoypUrxoABAwx/f3/D1dXVePTRR40JEybcssRRQkKC0adPH8PX19fw9PQ02rZtm75M0l+Xyrq5ZFNmSzGdOXPG+Mc//mH4+PgY3t7exksvvWRER0dnubzX358jNDTUKFKkyF19nbIyaNAgQ5Ixfvz4DNtvLlN17NixDNszW7bp6tWrxquvvmr4+PgYktKX6cpqaa3MlnjKTHaW5jp06JDRsGFDo1ChQoakOy7TdfLkSaNNmzZG4cKFjWLFihlvvfWWsWbNmkyXpNq7d6/Rrl0748EHHzTc3d2NcuXKGS+//LKxcePGWzL9/Rz9fUmpjRs3Gs8//7zh7+9vuLm5Gf7+/kaHDh2MI0eO3Pbrc7vv2Zt69+5tSDK++OKL2772v8rqe+XChQvGCy+8YBQuXNh44IEHjB49ehj//e9/M12aK7PvwczOUWpqqjFhwgSjYsWKhpubm1G8eHGjRYsWxu7du9PHZHUe//51vGnq1KlGxYoVDVdXV8PPz8/o1auXcenSpbt6jaGhoZkuKQdYjc0wmP0NIGfs27dPNWrU0KJFi/Taa6+ZHQcFzIABAzR79mydO3fulqvMAPIv5swCuCeJiYm3bJsyZYqcnJxue+ctIDdcv35dixYt0gsvvECRBQoY5swCuCeRkZHavXu3GjduLBcXF3377bf69ttv9cYbbyggIMDseCggzp8/rw0bNmjFihX6888/9dZbb5kdCYCDUWYB3JP69etr/fr1ev/993X16lWVLVtW7733nkaMGGF2NBQgBw8e1GuvvaYSJUroo48+Sl8ZAUDBwZxZAAAAWBZzZgEAAGBZlFkAAABYVoGbM2u32xUdHa2iRYvmiVtQAgAAICPDMHTlyhX5+/vLyen2114LXJmNjo7mndYAAAAWcPr0aZUpU+a2YwpcmS1atKikG18cLy8vk9MAAADg7+Lj4xUQEJDe226nwJXZm1MLvLy8KLMAAAB52N1MCeUNYAAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALitlJQ0HT580ewYmaLMAgAAIEtHjlzU008vVqNGSxQbe83sOLegzAIAAOAWhmFo+vS9ql59gXbuPKekpDQdPPin2bFu4WJ2AAAAAOQt0dFX1bXrGq1de0KSFBxcTnPnNleZMkXNDZYJyiwAAADSLV9+WD17rtfFi9fl4eGiyMiG6tOnhpycbGZHyxRlFgAAALp8+br69t2ozz//VZJUs6afFi5sqUqVHjQ52e1RZgEAAAq4TZtOKTT0W505c0VOTjaNGFFX775bT66uzmZHuyPKLAAAQAGVmJii4cO/15QpuyVJjzziowULWqpePX+Tk909yiwAAEABtHdvjDp2jEpfoaBHj2qaOLGRPD3dTE6WPZRZAACAAiQ11a7IyJ0KD/9Bqal2+fkV1uzZzdSqVZDZ0e4JZRYAAKCAOHbssjp1itIPP0RLktq1e1SfftpUxYoVNjnZvaPMAgAA5HOGYWjWrAMaMGCzEhJSVLSom6ZObaKQkMqy2fLmklt3izILAACQj8XEJKh797Vatep3SVKjRmU0f34LlSvnbXKynEGZBQAAyKdWrvxNr7++TrGxiXJzc9bYsU9rwIBaefYGCPeCMgsAAJDPxMcnqX//zZo797+SpGrVimvhwpaqUqW4yclyHmUWAAAgH9m69bRCQ7/ViRPxstmkwYPraNSo+nJ3z5+1L3++KgAAgAImKSlVI0du14QJu2QYUmCglxYsaKkGDcqYHS1XUWYBAAAsbv/+CwoJidL+/RckSd26VdHkyY1VtKi1boBwLyizAAAAFpWWZteHH/6sd97ZruTkNBUvXkgzZzbT888/YnY0h6HMAgAAWNCJE3EKDf1WW7eekSS1aROkmTOfU4kSRUxO5liUWQAAAAsxDEPz5/+iN9/cpCtXkuXp6aopU55V166PW/4GCPeCMgsAAGARFy5cU48e6/XVV79Jkp56qrQWLGihhx/2MTeYiSizAAAAFrBq1TF1775WMTHX5OrqpNGjn9KgQbXl7OxkdjRTUWYBAABywqHc+RP/1QQ3DYxsrc+WPSlJqhx0TosiF6tG5YHSb7lyyMxVNBx4sLtHmQUAAMijfthbTp2GvqJjp4pJksI6f6cx/dfIwz3V5GR5B2UWAAAgj0lOdtao6U01bmZj2e1OCih1SfMjlqpx3WNmR8tzKLMAAAB5yMGjJdRxcAft/fXGnbtC2uzWRyNWysfrusnJ8ibKLAAAQB5gt9v08aKnNGRSSyUlu8rXO0GfjvqXXmx2wOxoeRplFgAAwGSn//BWl+HttfHHRyVJzRsc0pwPlqlUiSsmJ8v7KLMAAAAmMQxp8erq6j26neKuFFLhQsmaNPgb9Wj/owrg/Q/uCWUWAADABBcvF1KvUe20bE11SVKdqqe0cNxilX8o1txgFkOZBQAAcLC135dX1xEvK/q8t5yd0xTee4OGvbFJLi52s6NZDmUWAADAQa4lumrwxFaa9sVTkqQKD53XosjFqvX4GZOTWRdlFgAAwAF2HSijjoM76MiJEpKkfh2/17iwKBUulGJyMmujzAIAAOSilBQnjf2sid7/pInS0pzlXyJO88YuVdOnHHkv2vyLMgsAAJBLjhwvppChHbRzf1lJ0ist92rau1/J1yfR5GT5B2UWAAAghxmG9Mnienp7wv8p8bqbfLyuafrIr9Sh1T6zo+U7TmYefOvWrWrdurX8/f1ls9m0cuXKu953+/btcnFxUfXq1XMtHwAAQHZFn/dSize6q8/77ZR43U3B9Y7owL8/pMjmElPLbEJCgqpVq6Zp06Zla7/Lly+rU6dOatKkSS4lAwAAyL7la6qqSpswrf2+gjzcU/TRiJVaO2uWypSMMztavmXqNIMWLVqoRYsW2d6vZ8+eevXVV+Xs7Jytq7kAAAC54fLl6+o7uIM+/+YJSdITlc9oUeRiVQo6b3Ky/M/UK7P3Yu7cufr9998VHh5+V+OTkpIUHx+f4QEAAJBTNm06pSpV5uvzb56Qk5Nd7/TcoB2Lp1JkHcRSbwD77bffNHToUG3btk0uLncXPSIiQqNGjcrlZAAAoKBJTEzR8OHfa8qU3ZKkR8pe0IJxS1WvxkmTkxUslrkym5aWpldffVWjRo1S+fLl73q/YcOGKS4uLv1x+vTpXEwJAAAKgr17Y1Sr1qL0ItujRzXt/XIKRdYElrkye+XKFf3888/au3ev+vbtK0my2+0yDEMuLi5at26dnn322Vv2c3d3l7u7u6PjAgCAfCgtza7IyF0KD9+ulBS7/PwKa/bsZmrVKkg6lGx2vALJMmXWy8tLBw4cyLBt+vTp2rRpk1asWKGHHnrIpGQAAKAgOHbssjp1itIPP0RLktq1e1SfftpUxYoVNjlZwWZqmb169aqOHj2a/vHx48e1b98++fr6qmzZsho2bJjOnj2rBQsWyMnJSY8//niG/UuUKCEPD49btgMAAOQUwzA0a9YBDRiwWQkJKSpa1E1TpzZRSEhl2Ww2s+MVeKaW2Z9//lmNGzdO/zgsLEySFBoaqnnz5umPP/7QqVOnzIoHAAAKuJiYBHXvvlarVv0uSWrUqIzmz2+hcuW8TU6Gm2yGYRhmh3Ck+Ph4eXt7Ky4uTl5eXmbHAQAAedTKlb/p9dfXKTY2UW5uzho79mkNGFBLTk5ZXI09lM+v0lZ0XGXMTl+zzJxZAAAAR4iPT1L//ps1d+5/JUnVqhXXwoUtVaVKcZOTITOUWQAAgP9v69bTCg39VidOxMtmkwYPrqNRo+rL3Z3KlFdxZgAAQIGXlJSqkSO3a8KEXTIMKTDQSwsWtFSDBmXMjoY7oMwCAIACbf/+CwoJidL+/RckSd26VdHkyY1VtKibyclwNyizAACgQEpLs+vDD3/WO+9sV3JymooXL6SZM5vp+ecfMTsasoEyCwAACpwTJ+IUGvqttm49I0lq0yZIM2c+pxIlipicDNlFmQUAAAWGYRiaP/8XvfnmJl25kqwiRVz1z38+q65dH+cGCBZFmQUAALkvD6zBeuFiEfUIf0FfbagiSapf44QWjFuioLL9pcPmZsO9o8wCAIB8b9XmSuo+8iXFxBaVq2uqRvdbp0Fdt8jZuUDdOypfoswCAIB862qCmwZGttZny56UJFUOOqdFkYtVo3K0ycmQUyizAAAgX9qxt5xChr6iY6eKSZLCOn+nMf3XyMM91eRkyEmUWQAAkK8kJztr9PRgRcx8Vna7kwJKXdL8iKVqXPeY2dGQCyizAAAg3zh4tIRChnTQnoM37twV0ma3PhqxUj5e101OhtxCmQUAAJZnt9v08aKnNGRSSyUlu8rXO0GfjvqXXmx2wOxoyGWUWQAAYGmn//BWl+HttfHHRyVJzRsc0pwPlqlUiSsmJ4MjUGYBAIAlGYa0eHV19R7dTnFXCqlwoWRNGvyNerT/Udz/oOCgzAIAAMu5eLmQeo1qp2VrqkuS6lQ9pYXjFqv8Q7HmBoPDUWYBAIClrP2+vLqOeFnR573l7Jym8N4bNOyNTXJxsZsdDSagzAIAAEu4luiqwRNbadoXT0mSKjx0XosiF6vW42dMTgYzUWYBAECet+tAGXUc3EFHTpSQJPXr+L3GhUWpcKEUk5PBbJRZAACQZ6WkOGnsZ030/idNlJbmLP8ScZo3dqmaPvWb2dGQR1BmAQBAnnTkeDGFDO2gnfvLSpJeablX0979Sr4+iSYnQ15CmQUAAHmKYUifLK6ntyf8nxKvu8nH65qmj/xKHVrtMzsa8iDKLAAAyDOiz3up64iXtfb7CpKkJk/+pnkRS1WmZJzJyZBXUWYBAECesHxNVfV8r50uxhWRh3uKxg9crb6v/SAnJ8PsaMjDKLMAAMBUl+M91PeDf+jzb56QJD1R+YwWRS5WpaDzJieDFVBmAQCAaTb9GKTQYa/ozDkfOTnZNfyNTXq31wa5uaWZHQ0WQZkFAAAOdz3JRcMnt9Dk+Q0lSY+UvaAF45aqXo2TJieD1VBmAQCAQ+096K+Ogzvo4LGSkqQe7Xdo4qBV8iySbHIyWBFlFgAAOERamk2Rs59R+NTnlJLiIr9iVzT7/WVq9cwhs6PBwiizAAAg1x079aA6DX1FP+wNlCS1a3pAn45aoWIPXDM3GCyPMgsAAHKNYRiaNeuABvQfoIRr7ipa5LqmvrNSIc/vls1mdjrkB5RZAACQK2JiEtS9+1qtWvW7JHc1qn1M8yOWqFzpy2ZHQz5CmQUAADlu5crf9Prr6xQbmyg3N2eN7b9SA0K3cQME5DjKLAAAyDHx8Unq33+z5s79rySpWrXiWriwpaq4hpmcDPkVZRYAAOSIrVtPKzT0W504ES+bTRo8uI5Gjaovd3cXiQULkEsoswAA4L4kJaVq5MjtmjBhlwxDCgz00oIFLdWgQRmzo6EAoMwCAIB7tn//BYWERGn//guSpG7dqmjy5MYqWtTN5GQoKCizAAAg29LS7Prww5/1zjvblZycpuLFC2nmzGZ6/vlHzI6GAoYyCwAAsuXEiTiFhn6rrVvPSJJatw7SzJnPyc+viMnJUBBRZgEAwF0xDEPz5/+iN9/cpCtXklWkiKumTGmsbt2qyMYdEGASyiwAALijCxeuqUeP9frqq98kSfXr+2vBgpYKCvIxNxgKPMosAAC4rVWrjql797WKibkmV1cnjR79lAYNqi1nZyezowGUWQAAkLmrV5M1cOAWffbZfklS5coPatGilqpRw8/cYMBfUGYBAMAtduyIVkhIlI4duyxJCgurqTFjGsjDg+qAvIXvSAAAkC45OU2jR+9QRMRPstsNBQQU1fz5LdS4cVmzowGZoswCAABJ0sGDsQoJ+VZ79sRIkkJCKuujj56Vj4+HycmArFFmAQAo4Ox2Qx9/vEdDhmxVUlKafH099OmnTfXiixXMjgbcEWUWAIAC7PTpeHXpskYbN56SJDVvHqg5c5qrVClPk5MBd4cyCwBAAWQYhhYvPqTevTcoLi5JhQu7aNKkZ9SjRzVugABLocwCAFDAXLyYqF69NmjZssOSpDp1SmrhwpYqX97X5GRA9lFmAQAoQNauPa6uXdcqOvqqnJ1tCg+vr2HD6srFhRsgwJooswAAFADXrqVo8ODvNG3aPklShQq+WrSopWrVKmluMOA+UWYBAMgLDuXePNVdB8qo4+AOOnKihCSpX8fvNS4sSoULpUiHcu2wgENQZgEAyKdSUpw09rMmev+TJkpLc5Z/iTjNG7tUTZ/6zexoQI6hzAIAkA8dOV5MIUM7aOf+G3fuat9in6aP/FK+PokmJwNyFmUWAIB8xDCkTxbX09sT/k+J193k43VN00d+pQ6t9pkdDcgVlFkAAPKJ6PNe6jriZa39/sadu5o8+ZvmRSxVmZJxJicDcg9lFgCAfGD5mqrq+V47XYwrIg/3FI0fuFp9X/tBTk6G2dGAXEWZBQDAwi7He6jfB2216JuakqQnKp/RosjFqhR03uRkgGNQZgEAsKhNPwap8/D2Ov3HA3Jysmv4G5v0bq8NcnNLMzsa4DCUWQAALOZ6kouGT26hyfMbSpIeKXtBC8YtVb0aJ01OBjgeZRYAAAvZe9BfHQd30MFjN+7c1aP9Dk0ctEqeRZJNTgaYgzILAIAFpKXZFDn7GYVPfU4pKS7yK3ZFs99fplbPcAsvFGyUWQAA8rhjpx5Up6Gv6Ie9gZKkdk0P6NNRK1TsgWvmBgPyAMosAAB5lGFIs5bX0YDxbZRwzV1Fi1zX1HdWKuT53bLZzE4H5A2UWQAA8qCYWE91f/clrdpSWZLUqPYxzY9YonKlL5sbDMhjnMw8+NatW9W6dWv5+/vLZrNp5cqVtx3/5ZdfqmnTpipevLi8vLxUr149rV271jFhAQBwkJUbHtPjbQZq1ZbKcnNN1cTB32jTvE8pskAmTC2zCQkJqlatmqZNm3ZX47du3aqmTZsqKipKu3fvVuPGjdW6dWvt3bs3l5MCAJD74q+6q+uIl/SPfp0Ve8lT1SpG6+cV/9TALlu5kxeQBZthGHnip8Nms+mrr75S27Zts7XfY489pvbt22vkyJF3NT4+Pl7e3t6Ki4uTl5fXPSQFACDnbV34sEKHvaITZ31ls9k1uNsWjeq3Tu7cAAF5RUXHVcbs9DVLz5m12+26cuWKfH19sxyTlJSkpKSk9I/j4+MdEQ0AgLuSlJSqkSO3a8KEnjIMJwWWvqgF45aoQa3jZkcDLMHUaQb3a+LEibp69apefvnlLMdERETI29s7/REQEODAhAAAZG3//guqU+dzRUbukmE4qWu7nfrPyg8pskA2WLbMfvHFFxo1apSWLVumEiVKZDlu2LBhiouLS3+cPn3agSkBALhVWppdEybsVO3ai7R//wUVL15IK6fO0+wxy+XlmXTnJwCQzpLTDJYsWaLu3btr+fLlCg4Ovu1Yd3d3ubu7OygZAAC3d+JEnEJDv9XWrWckSa1bB2nmzOfkd6mvyckAa7JcmV28eLG6du2qJUuWqFWrVmbHAQDgrhiGofnzf9Gbb27SlSvJKlLEVVOmNFa3blVks9mkS2YnBKzJ1DJ79epVHT16NP3j48ePa9++ffL19VXZsmU1bNgwnT17VgsWLJB0Y2pBaGio/vnPf6pu3bo6d+6cJKlQoULy9vY25TUAAHAnFy5cU48e6/XVV79JkurX99eCBS0VFORjbjAgHzB1zuzPP/+sGjVqqEaNGpKksLAw1ahRI32ZrT/++EOnTp1KH//ZZ58pNTVVffr0UalSpdIfb731lin5AQC4k9Wrj6lKlXn66qvf5OrqpIiIBtq69RWKLJBD8sw6s47COrMAAEe4ejVZAwdu0Wef7ZckVa78oBYtaqkaNfwy3+GQzXHhgHvBOrMAABQMO3ZEKyQkSseOXZYkhYXV1JgxDeThwa9dIKfxUwUAQA5JTk7T6NE7FBHxk+x2QwEBRTV/fgs1blzW7GhAvkWZBQAgBxw8GKuQkG+1Z0+MJCkkpLI++uhZ+fh4mJwMyN8oswAA3Ae73dDHH+/RkCFblZSUJl9fD336aVO9+GIFs6MBBQJlFgCAe3T6dLy6dFmjjRtvrLzTvHmg5sxprlKlPE1OBhQclFkAALLJMAwtXnxIvXtvUFxckgoXdtGkSc+oR49qN26AAMBhKLMAAGTDxYuJ6tVrg5YtOyxJqlOnpBYubKny5X1NTgYUTJRZAADu0tq1x9W161pFR1+Vs7NN4eH1NWxYXbm4mHoPIqBAo8wCAHAH166laPDg7zRt2j5JUoUKvlq0qKVq1SppbjAAlFkAAG5n164/1LFjlI4cuSRJ6tevhsaNa6jChV1NTgZAoswCAJCplJQ0jR37k95/f4fS0gz5+3tq7tzmeu65QLOjAfgLyiwAAH9z5MhFhYREaefOc5Kk9u0raPr0YPn6FjI5GYC/o8wCAPD/GYahTz7Zp7ff/k6Jiany8XHX9OnB6tChktnRAGSBMgsAsIZDubt+a/R5L3Ud8bLWfn/jzl1NnvxN8yKWqkzJN6VDuXpoAPeBMgsAKPCWr6mqnu+108W4IvJwT9H4gavV97Uf5ORkmB0NwB1QZgEABdbleA/1+6CtFn1TU5L0ROUzWhS5WJWCzpucDMDdoswCAAqkTT8GqfPw9jr9xwNycrJr+Bub9G6vDXJzSzM7GoBsoMwCAAqU60kuGj65hSbPbyhJeqTsBS0Yt1T1apw0ORmAe0GZBQAUGHsP+qvj4A46eOzGnbt6tN+hiYNWybNIssnJANwryiwAIN9LS7MpcvYzCp/6nFJSXORX7Ipmv79MrZ5hmQLA6iizAIB87dipB9Vp6Cv6YW+gJKld0wP6dNQKFXvgmrnBAOQIyiwAIF8yDGnW8joaML6NEq65q2iR65r6zkqFPL9bttxdshaAA1FmAQD5Tkysp7q/+5JWbaksSWpU+5jmRyxRudKXzQ0GIMdRZgEA+crKDY/p9ZEvKvaSp9xcUzV2wLcaELqNGyAA+RRlFgCQL8RfdVf/iDaa+2UdSVK1itFaOH6xqpQ/Z3IyALmJMgsAsLytux5S6LBXdOKsr2w2uwZ326JR/dbJnRsgAPkeZRYAYFlJyc4a+VEzTZjTSIbhpMDSF7Vg3BI1qHXc7GgAHIQyCwCwpP2HSylkyCvaf9hfktS13U5NHva1vDyTTE4GwJEoswAAS0lLs+nDeQ31zj+bKznFRcV9r2rm6BV6vskvZkcDYALKLADAMk6cfUChQ9tr689BkqTWjX/RzNEr5FfsqsnJAJiFMgsAyPMMw9CClTXV74O2upLgoSKFkzRl6Nfq9uJOboAAFHCUWQBAnhYbe009eqzXl1++IkmqX+OEFoxboqCyf5qcDEBeQJkFAORZq1cfU7duaxUTc02urqka3W+dBnXdImdnboAA4AbKLAAgz7l6NVkDB27RZ5/tlyRVrvygFr3/jmpUjjY3GIA8x8nsAAAA/NWOHdGqXn1BepENC6up3btDKLIAMsWVWQBAnpCcnKbRo3coIuIn2e2GAgKKav78FmrcuKzZ0QDkYZRZAIDpDh6MVUjIt9qzJ0aSFBJSWR999Kx8fDxMTgYgr6PMAgBMY7cb+vjjPRoyZKuSktLk6+uhTz9tqhdfrGB2NAAWQZkFAJji9Ol4demyRhs3npIkNW8eqDlzmqtUKU+TkwGwEsosAMChDMPQ4sWH1Lv3BsXFJalwYRdNmvSMevSoJht3QACQTZRZAIDDXLyYqF69NmjZssOSpDp1SmrhwpYqX97X5GQArIoyCwBwiLVrj6tr17WKjr4qZ2ebwsPra9iwunJxYZVIAPeOMgsAyFXXrqVo8ODvNG3aPklShQq+WriwhWrXLmVuMAD5AmUWAJBrdu36Qx07RunIkUuSpL59a2j8+IYqXNjV5GQA8gvKLAAgx6WkpGns2J/0/vs7lJZmyN/fU3PnNtdzzwWaHQ1APkOZBQDkqCNHLiokJEo7d56TJLVvX0HTpwfL17eQyckA5EeUWQBAjjAMQzNm/EcDB25RYmKqfHzcNX16sDp0qGR2NAD5GGUWAHDfoqOvqlu3NVqz5oQkqUmTspo3r4XKlClqbjAA+R5lFgDyg0Pm3Wxg+Zqq6vleO12MKyIP9xSNH7hafV/7QU5XDemQabEAFBCUWQDAPbkc76F+H7TVom9qSpKeqHxGiyIXq1LQeZOTAShIKLMAgGzb9GOQOg9vr9N/PCAnJ7uGv7FJ7/baIDe3NLOjAShgKLMAgLt2PclFwye30OT5DSVJj5S9oAXjlqpejZMmJwNQUFFmAQB3Ze9Bf3Uc3EEHj5WUJPVov0MTB62SZ5Fkk5MBKMgoswCA20pLsyly9jMKn/qcUlJc5Ffsima/v0ytnuHdXQDMR5kFAGTp2KkH1WnoK/phb6AkqV3TA/p01AoVe+CaucEA4P+jzAIAbmEY0qzldTRgfBslXHNX0SLXNfWdlQp5frds5q0CBgC3oMwCADKIifVU93df0qotlSVJjWof0/yIJSpX+rK5wQAgE5RZAEC6lRse0+sjX1TsJU+5uaZq7IBvNSB0m5ycDLOjAUCmKLMAAMVfdVf/iDaa+2UdSVLVCtFaFLlYVcqfMzkZANweZRYACritux5S6LBXdOKsr2w2uwZ326JR/dbJnRsgALAAyiwAFFBJyc4a+VEzTZjTSIbhpMDSF7Vg3BI1qHXc7GgAcNcoswBQAO0/XEohQ17R/sP+kqSu7XZq8rCv5eWZZHIyAMgeyiwAFCBpaTZNnt9QI6Y0V3KKi4r7XtXM0Sv0fJNfzI4GAPeEMgsABcTJsz4KHfaKvtsVJElq3fgXzRy9Qn7FrpqcDADuHWUWAPI5w5AW/Lum+n3QVlcSPFSkcJKmDP1a3V7cyQ0QAFgeZRYA8rHYS4XVI/xFfbm+iiSpfo0TWjBuiYLK/mlyMgDIGZRZAMinVm+pqG7vvqyY2KJydU3V6H7rNKjrFjk7cwMEAPkHZRYA8pmrCW4aGNlany17UpJUOeicFkUuVo3K0SYnA4Cc52Tmwbdu3arWrVvL399fNptNK1euvOM+W7Zs0RNPPCF3d3c98sgjmjdvXq7nBACr2LG3nKq3G5BeZMM6f6fd//onRRZAvmVqmU1ISFC1atU0bdq0uxp//PhxtWrVSo0bN9a+ffvUv39/de/eXWvXrs3lpACQtyUnO+udKc30dMfeOnaqmAJKXdKmeTM0acgqebinmh0PAHKNqdMMWrRooRYtWtz1+BkzZuihhx7SpEmTJEmVKlXS999/r8mTJ6tZs2a5FRMA8rSDB2MV0qGv9hwsI0kKabNbH41YKR+v6yYnA4DcZ6k5szt27FBwcHCGbc2aNVP//v2z3CcpKUlJSf+7o018fHxuxQMAh7LbDX388R4NGbJVSUll5OudoE9H/UsvNjtgdjQAcBhTpxlk17lz5+Tn55dhm5+fn+Lj45WYmJjpPhEREfL29k5/BAQEOCIqAOSq06fj9dxzy9W//2YlJaWpeYND+u/XkyiyAAocS5XZezFs2DDFxcWlP06fPm12JAC4Z4Zh6IsvflWVKvO1ceMpFS7sok8+CVbUp7NVqsQVs+MBgMNZappByZIlFRMTk2FbTEyMvLy8VKhQoUz3cXd3l7u7uyPiAUCuungxUb16bdCyZYclSXXqlNTChS1VvryvdMjkcABgEkuV2Xr16ikqKirDtvXr16tevXomJQIAx1i79ri6dl2r6Oircna2aeTIeho+/Em5uOT7P7ABwG2ZWmavXr2qo0ePpn98/Phx7du3T76+vipbtqyGDRums2fPasGCBZKknj17aurUqRo8eLC6du2qTZs2admyZVq9erVZLwEActW1aykaPPg7TZu2T5JUoYKvFi5sodq1S5kbDADyCFPL7M8//6zGjRunfxwWFiZJCg0N1bx58/THH3/o1KlT6Z9/6KGHtHr1ag0YMED//Oc/VaZMGc2aNYtluQDkS7t2/aGOHaN05MglSVLfvjU0fnxDFS7sanIyAMg7bIZhFKibdMfHx8vb21txcXHy8vIyOw4A3CI11a6xY3/U6NE7lJZmyN/fU3PnNtdzzwVmvdMhm8PyASigKjquMmanr1lqziwA5HdHjlxUSEiUdu48J0lq376Cpk8Plq9v5m9yBYCCjjILAHmAYRiaMeM/GjhwixITU+Xj467p04PVoUMls6MBQJ5GmQUAk0VHX1W3bmu0Zs0JSVKTJmU1b14LlSlT1NxgAGABlFkAMNHy5YfVs+d6Xbx4XR4eLho/vqH69q0hJyfmwALA3aDMAoAJLl++rn79NmnRooOSpCee8NOiRS1VqdKDJicDAGuhzAKAg23adEqdO3+r06evyMnJpuHD6+rdd+vJzc3Z7GgAYDmUWQBwkOvXUzV8+DZNnrxbkvTIIz5asKCl6tXzNzkZAFgXZRYAHGDv3hh17Bilgwf/lCT16FFNEyc2kqenm8nJAMDaKLMAkIvS0uyKjNyl8PDtSkmxy8+vsGbPbqZWrYLMjgYA+QJlFgByybFjl9WpU5R++CFaktSu3aP69NOmKlassMnJACD/oMwCQA4zDEOzZh3QgAGblZCQoqJF3TR1ahOFhFSWzcaSWwCQkyizAJCDYmIS1L37Wq1a9bskqWHDMlqwoIXKlfM2ORkA5E+UWQAFw6HcvyK6csNjen3ki4q95Ck311SN6b9GA0K3yjnRkA7l+uEBoECizALAfYq/6q7+EW0098s6kqSqFaK1KHKxqpQ/Z3IyAMj/KLMAcB+27npIocNe0YmzvrLZ7BrcbYtG9Vsnd7c0s6MBQIFAmQWAe5CU7KyRHzXThDmNZBhOCix9UQvGLVGDWsfNjgYABQplFgCy6cCRkuo4uIP2H75x566u7XZq8rCv5eWZZHIyACh4KLMAcJfS0myaPL+hRkxpruQUFxX3vaqZo1fo+Sa/mB0NAAosyiwA3IWTZ30UOuwVfbfrxp27Wjf+RTNHr5BfsasmJwOAgo0yCwC3YRjSgn/XVL8P2upKgoeKFE7SlKFfq9uLO8X9DwDAfJRZAMhC7KXC6hH+or5cX0WSVL/GCS0Yt0RBZf80ORkA4CbKLABkYvWWiur27suKiS0qV9dUje63ToO6bpGzs2F2NADAX1BmAeAvria4aWBka3227ElJUuWgc1oUuVg1KkebnAwAkBnKLAD8fzv2llPI0Fd07FQxSVJY5+80pv8aebinmpwMAJAVyiyAAi852VmjpwcrYuazstudFFDqkuZHLFXjusfMjgYAuAPKLIAC7eDREgoZ0kF7DpaRJIW02a2PRqyUj9d1k5MBAO4GZRZAgWS32/Txoqc0ZFJLJSW7ytc7QZ+O+pdebHbA7GgAgGygzAIocE7/4a0uw9tr44+PSpKaNzikOR8sU6kSV0xOBgDILsosgALDMKTFq6ur9+h2irtSSIU8kjVp8Cr1fGUHN0AAAIuizAIoEC5eLqReo9pp2ZrqkqQ6VU9p4bjFKv9QrLnBAAD3hTILIN9bu/a4unYaqOjz3nJ2TtPIXhs0vMcmubjYzY4GALhPlFkA+da1aykaMmSrpk7dK8lbFR46r4XjF6t2lTNmRwMA5BDKLIB8adeuPxQS8q0OH74oSer72vcaPzBKhQulmJwMAJCTKLMA8pXUVLvGjv1Ro0fvUFqaIX9/T82d21zPlR1kdjQAQC6gzALIN44cuaiQkCjt3HlOktS+fQVNnx4sX99C0iGTwwEAcgVlFoDlGYahGTP+o4EDtygxMVU+Pu6aPj1YHTpUMjsaACCXUWYBWFp09FV167ZGa9ackCQ1aVJW8+a1UJkyRc0NBgBwCMosAMtavvywevZcr4sXr8vDw0XjxzdU37415OTEHRAAoKCgzAKwnMuXr6tfv01atOigJOmJJ/y0aFFLVar0oMnJAACORpkFYCmbNp1S587f6vTpK3Jysmn48Lp69916cnNzNjsaAMAElFkAlnD9eqqGD9+myZN3S5IeecRHCxa0VL16/iYnAwCYiTILIM/buzdGHTtG6eDBPyVJPXpU08SJjeTp6WZyMgCA2SizAPKstDS7IiN3KTx8u1JS7PLzK6zZs5upVasgs6MBAPIIyiyAPOnYscvq1ClKP/wQLUlq1+5RffppUxUrVtjkZACAvIQyCyBPMQxDs2Yd0IABm5WQkKKiRd00dWoThYRUls3GklsAgIwoswDyjJiYBHXvvlarVv0uSWrYsIzmz2+hwEBvk5MBAPIqyiyAPGHlyt/0+uvrFBubKDc3Z40Z87QGDKgpZ2cns6MBAPIwyiwAU8XHJ6l//82aO/e/kqSqVYtr0aKWqlKluMnJAABWQJkFYJpt286oU6conTgRL5tNGjy4jkaNqi93d/7TBAC4O/zGAOBwSUmpGjlyuyZM2CXDkAIDvbRgQUs1aFDG7GgAAIuhzAJwqAMHLqhjxyjt339BktS16+OaPLmxvLzcTU4GALAiyiyAGw7l7rJXaWk2TZ7fUCOmNFdyiouK+17VzNEr9HyTX6Ro3XgAAJBNlFkAue7kWR+FDntF3+26ceeu1o1/0czRK+RX7KrJyQAAVkeZBZBrDENa8O+a6vdBW11J8FCRwkmaMvRrdXtxp7j/AQAgJ1BmAeSK2EuF1SP8RX25vookqX6NE1owbomCyv5pcjIAQH5CmQWQ41Zvqahu776smNiicnVN1eh+6zSo6xY5OxtmRwMA5DOUWQA55mqCmwZGttZny56UJFUOOqdFkYtVozLv7gIA5A7KLIAcsWNvOYUMfUXHThWTJIV1/k5j+q+Rh3uqyckAAPkZZRbAfUlOdtbo6cGKmPms7HYnBZS6pPkRS9W47jGzowEACgDKLIB7dvBoCYUM6aA9B2/cuSukzW59NGKlfLyum5wMAFBQUGYBZJvdbtPHi57SkEktlZTsKl/vBM1470u91Hy/2dEAAAUMZRZAtpz+w1tdhrfXxh8flSQ1e/qw5oxZJv8S8SYnAwAURJRZAHfFMKTFq6ur9+h2irtSSIU8kjVp8Cr1fGUHN0AAAJiGMgvgji5eLqReo9pp2ZrqkqQ6VU9p4bjFKv9QrLnBAAAFHmUWwG2t215eXYa/rOjz3nJ2TtPIXhs0vMcmubjYzY4GAABlFkDmriW6asiklpr6+dOSpAoPndfC8YtVu8oZk5MBAPA/lFkAt9h1oIxChnTQ4eMlJEl9X/te4wdGqXChFJOTAQCQkZPZAaZNm6bAwEB5eHiobt262rlz523HT5kyRRUqVFChQoUUEBCgAQMG6Pp11rQEckJqqpNGTwtWvQ59dfh4CfmXiNPaWTP18Tv/psgCAPIkU6/MLl26VGFhYZoxY4bq1q2rKVOmqFmzZjp8+LBKlChxy/gvvvhCQ4cO1Zw5c1S/fn0dOXJEnTt3ls1m04cffmjCKwDyjyPHiylkaAft3F9WktS+xT5NH/mlfH0STU4GAEDWbIZhGGYdvG7duqpdu7amTp0qSbLb7QoICFC/fv00dOjQW8b37dtXv/76qzZu3Ji+beDAgfrpp5/0/fff39Ux4+Pj5e3trbi4OHl5eeXMCwEszDAMzZjxHw0Mi1LidTf5eF3T9JFfqUOrfWZHAwDkJRUdVxmz09dMm2aQnJys3bt3Kzg4+H9hnJwUHBysHTt2ZLpP/fr1tXv37vSpCL///ruioqLUsmXLLI+TlJSk+Pj4DA8AN0RHX1XLlv9S794blHjdTU2e/E0H/v0hRRYAYBmmTTOIjY1VWlqa/Pz8Mmz38/PToUOHMt3n1VdfVWxsrJ5++mkZhqHU1FT17NlTw4cPz/I4ERERGjVqVI5mB/KD5csPq2fP9bp48bo8PFw0PmyF+r72g5ycTPtjDQAA2Wb6G8CyY8uWLRo7dqymT5+uPXv26Msvv9Tq1av1/vvvZ7nPsGHDFBcXl/44ffq0AxMDec/ly9cVEhKll1/+RhcvXtcTT/hpz54QvRmynSILALCcbF+ZDQ0NVbdu3dSwYcP7OnCxYsXk7OysmJiYDNtjYmJUsmTJTPd59913FRISou7du0uSqlSpooSEBL3xxhsaMWKEnJxu7ebu7u5yd3e/r6xAfrFp0yl17vytTp++Iicnm4YPr6t3360nNzdnKfM/iAAAkKdl+8psXFycgoOD9eijj2rs2LE6e/bsPR3Yzc1NNWvWzPBmLrvdro0bN6pevXqZ7nPt2rVbCquzs7OkG29iAZC569dTFRa2WU2aLNPp01f0yCM++v77Dnr//advFFkAACwq22V25cqVOnv2rHr16qWlS5cqMDBQLVq00IoVK5SSkr11KMPCwjRz5kzNnz9fv/76q3r16qWEhAR16dJFktSpUycNGzYsfXzr1q31ySefaMmSJTp+/LjWr1+vd999V61bt04vtQAy2rs3RjVrLtTkybslST16VNPevZ1Ur56/yckAALh/9/QGsOLFiyssLExhYWHas2eP5s6dq5CQEHl6eqpjx47q3bu3Hn300Ts+T/v27XXhwgWNHDlS586dU/Xq1bVmzZr0N4WdOnUqw5XYd955RzabTe+8847Onj2r4sWLq3Xr1hozZsy9vAwgX0tLsysycpfCw7crJcUuP7/Cmj27mVq1CjI7GgAAOea+1pn9448/tGDBAs2dO1dnzpzRCy+8oLNnz+q7775TZGSkBgwYkJNZcwTrzKIgOHbssjp1itIPP0RLkv7xj0f16adNVbx44ax3OmRzUDoAgCXl0XVms31lNiUlRV9//bXmzp2rdevWqWrVqurfv79effXV9IN99dVX6tq1a54ss0B+ZhiGZs06oAEDNishIUVFi7rp44+fVadOj8lmo6wCAPKfbJfZUqVKyW63q0OHDtq5c6eqV69+y5jGjRvLx8cnB+IBuFsxMQnq3n2tVq36XZLUsGEZzZ/fQoGB3iYnAwAg92S7zE6ePFkvvfSSPDw8shzj4+Oj48eP31cwAHdv5crf9Prr6xQbmyg3N2eNGfO0BgyoKWdnSy0lDQBAtmW7zIaEhORGDgD3ID4+SQMGbNacOf+VJFWtWlyLFrVUlSrFTU4GAIBjmHY7WwD3Z9u2M+rUKUonTsTLZpMGD66jUaPqy92dH2sAQMHBbz3AYpKSUjVy5HZNmLBLhiEFBnppwYKWatCgjNnRAABwOMosYCEHDlxQx45R2r//giSpa9fHNXlyY3l5cctmAEDBRJkFLCAtza7Jk3drxIjvlZycpuLFC2nmzGZ6/vlHzI4GAICpKLNAHnfyZJxCQ7/Vd9+dkSS1bh2kmTOfk59fEZOTAQBgPsoscDdMuDuWYUgL/l1T/T5oqysJHipSOElThn6tbi/ulO2SpEsOjwQAQJ5DmQXyoNhLhdUj/EV9ub6KJKl+jRNaMG6Jgsr+aXIyAADyFsoskMes3lJR3d59WTGxReXqmqrR/dZpUNctcnZ23D2xAQCwCsoskEdcTXDTwMjW+mzZk5KkykHntChysWpUjjY5GQAAeRdlFsgDduwtp5Chr+jYqWKSpLDO32lM/zXycE81ORkAAHkbZRYwUXKys0ZPD1bEzGdltzspoNQlzRu7VM8+eczsaAAAWAJlFjDJwaMlFDKkg/YcvHHnro6td+vjd1bKx+u6yckAALAOyizgYHa7TR8vekpDJrVUUrKrfL0TNOO9L/VS8/1mRwMAwHIos4ADnf7DW12Gt9fGHx+VJDV7+rDmjFkm/xLxJicDAMCaKLOAg3yxqrp6j26nuCuFVMgjWZMGr1LPV3bI5vj7MQAAkG9QZoFcdvFyIfUe3U5Lv60uSapT9ZQWjlus8g/FmhsMAIB8gDIL5KJ128ury/CXFX3eW87OaRrZa4OG99gkFxe72dEAAMgXKLNALriW6Kohk1pq6udPS5IqPHReC8cvVu0qZ0xOBgBA/kKZBXLYrgNlFDKkgw4fLyFJ6vva9xo/MEqFC6WYnAwAgPyHMgvkkNRUJ4399FmN/iRYaWnO8i8Rp7ljl+m5p46YHQ0AgHyLMgvkgCPHiylkaAft3F9WktS+xT5NH/mlfH0STU4GAED+RpkF7oNhSDOW1NPAyP9T4nU3+Xhd0/SRX6lDq31mRwMAoECgzAL3KPq8l7q985LWbKsoSWry5G+aF7FUZUrGmZwMAICCgzIL3IPla6qq53vtdDGuiDzcUzR+4Gr1fe0HOTkZZkcDAKBAocwC2XA53kP9PmirRd/UlCQ9UfmMFkUuVqWg8yYnAwCgYKLMAndp049B6jy8vU7/8YCcnOwa/sYmvdtrg9zc0syOBgBAgUWZBe7g+vVUDR/XWpPnN5QkPVL2ghaMW6p6NU6anAwAAFBmgdvYuzdGHTtG6eDBG0W2R/sdmjholTyLJJucDAAASJRZIFNpaXZFRu5SePh2paTY5Vfsima/v0ytnjlkdjQAAPAXlFngb44du6xOnaL0ww/RkqR//ONRffr2ayrum2ByMgAA8HeUWeD/MwxDs2cfUP/+m5WQkKKiRd308cfPqlOnx2Q7TJEFACAvoswCkmJiEvT66+v0zTfHJEkNG5bR/PktFBjobXIyAABwO5RZFHgrV/6mN95YpwsXEuXm5qwxY57WgAE15ezsZHY0AABwB5RZFFjx8UkaMGCz5sz5rySpatXiWrSopapUKW5yMgAAcLcosyiQtm07o06donTiRLxsNmnw4DoaNaq+3N35kQAAwEr4zY0CJSkpVSNHbteECbtkGFJgoJcWLGipBg3KmB0NAADcA8osCowDBy6oY8co7d9/QZLUtevjmjy5sby83E1OBgAA7hVlFvleWppdkyfv1ogR3ys5OU3FixfSzJnN9Pzzj5gdDQAA3CfKLPK1kyfjFBr6rb777owkqXXrIM2c+Zz8/IqYnAwAAOQEyizyJcMwtGDBL+rXb5OuXElWkSKumjKlsbp1qyKbzWZ2PAAAkEMos8h3YmOvqUeP9fryy98kSfXr+2vBgpYKCvIxNxgAAMhxlFnkK6tXH1O3bmsVE3NNrq5OGj36KQ0aVJsbIAAAkE9RZpEvXL2arIEDt+izz/ZLkipXflCLFrVUjRp+5gYDAAC5ijILy9uxI1ohIVE6duyyJCksrKbGjGkgDw++vQEAyO/4bQ/LSk5O0+jROxQR8ZPsdkMBAUU1b14LPftsWbOjAQAAB6HMwpIOHoxVSMi32rMnRpLUsWNlffzxs/Lx8TA5GQAAcCTKLCzFbjc0depeDRmyVdevp8rX10MzZjTVSy9VMDsaAAAwAWUWlnHmzBV17vytNm48JUlq1ixQc+Y0l7+/p8nJAACAWSizsITFi39V794bdPlykgoVctGkSc+oZ89q3AABAIACjjKLPO3ixUT17r1BS5celiTVqVNSCxe2VPnyviYnAwAAeQFlFnnWunUn1KXLGkVHX5Wzs00jR9bT8OFPysWFGyAAAIAbKLPIGYdy7s/91xJdNWRSS039/GlJUoWHzmvh+MWqXeVt6WiOHQYAAOQDlFnkKbsOlFHIkA46fLyEJKnva99r/MAoFS6UYnIyAACQF1FmkSekpjpp7KfPavQnwUpLc5Z/iTjNHbtMzz11xOxoAAAgD6PMwnRHjhdTyNAO2rn/xp272rfYp+kjv5SvT6LJyQAAQF5HmYVpDEOasaSeBkb+nxKvu8nH65qmj/xKHVrtMzsaAACwCMosTBF93kvd3nlJa7ZVlCQ1efI3zYtYqjIl40xOBgAArIQyC4dbvqaqer7XThfjisjDPUXjB65W39d+kJOTYXY0AABgMZRZOMzleA/1+6CtFn1TU5L0ROUzWhS5WJWCzpucDAAAWBVlFg6x6ccgdR7eXqf/eEBOTnYNf2OT3u21QW5uaWZHAwAAFkaZRa66nuSi4ZNbaPL8hpKkoLKxWjhuierVOGlyMgAAkB9QZpFr9h70V8fBHXTwWElJUo/2OzRx0Cp5Fkk2ORkAAMgvKLPIcWlpNkXOfkbhU59TSoqL/Ipd0ez3l6nVM4fMjgYAAPIZJ7MDTJs2TYGBgfLw8FDdunW1c+fO246/fPmy+vTpo1KlSsnd3V3ly5dXVFSUg9LiTn4/7auGIb01fHJLpaS46B/BB3Tg35MosgAAIFeYemV26dKlCgsL04wZM1S3bl1NmTJFzZo10+HDh1WiRIlbxicnJ6tp06YqUaKEVqxYodKlS+vkyZPy8fFxfHhkYBjS7BV11H9cGyVcc1fRItf18Tsr1en53bLZzE4HAADyK5thGKYt7lm3bl3Vrl1bU6dOlSTZ7XYFBASoX79+Gjp06C3jZ8yYoQkTJujQoUNydXW9p2PGx8fL29tbcXFx8vLyuq/8uCEmJkGvd+irbzY/JklqWOuY5o9bqsDSl0xOBgAAckxFx1XG7PQ106YZJCcna/fu3QoODv5fGCcnBQcHa8eOHZnu8/XXX6tevXrq06eP/Pz89Pjjj2vs2LFKS8t6eaekpCTFx8dneCDnrFz5m6pUmadvNj8mN9dUTRi0SpvmfUqRBQAADmHaNIPY2FilpaXJz88vw3Y/Pz8dOpT5/Mrff/9dmzZt0muvvaaoqCgdPXpUvXv3VkpKisLDwzPdJyIiQqNGjcrx/AVdfHySBgzYrDlz/itJqlohWosiF6tK+XMmJwMAAAWJpVYzsNvtKlGihD777DM5OzurZs2aOnv2rCZMmJBlmR02bJjCwsLSP46Pj1dAQICjIudL27adUadOUTpxIl42mzR4cB2NevVZuXMDBAAA4GCmldlixYrJ2dlZMTExGbbHxMSoZMmSme5TqlQpubq6ytnZOX1bpUqVdO7cOSUnJ8vNze2Wfdzd3eXu7p6z4QuopKRUjRy5XRMm7JJhSIGBXlqwoKUaNCgjHaLIAgAAxzNtzqybm5tq1qypjRs3pm+z2+3auHGj6tWrl+k+Tz31lI4ePSq73Z6+7ciRIypVqlSmRRY558CBC6pT53NFRt4osl27Pq7//Cf0RpEFAAAwianrzIaFhWnmzJmaP3++fv31V/Xq1UsJCQnq0qWLJKlTp04aNmxY+vhevXrp4sWLeuutt3TkyBGtXr1aY8eOVZ8+fcx6CfleWppdEyfuUq1ai7R//wUVL15IK1e21ezZzeXlxRVvAABgLlPnzLZv314XLlzQyJEjde7cOVWvXl1r1qxJf1PYqVOn5OT0v74dEBCgtWvXasCAAapatapKly6tt956S0OGDDHrJeRrJ0/GKTT0W3333RlJUuvWQZo58zn5+RUxORkAAMANpq4zawbWmb0zwzC0YMEv6tdvk65cSVaRIq6aMqWxunWrIltWd0A4xJ0RAADI1/LoOrOWWs0AuS829pp69FivL7/8TZJUv76/FixoqaAgH3ODAQAAZIIyi3SrVx9Tt25rFRNzTS4uTho9+ikNHlxbzs6mTq0GAADIEmUWuno1WQMHbtFnn+2XJFWu/KAWLWqpGjX8br8jAACAySizBdyOHdEKCYnSsWOXJUkDBtTU2LEN5OHBtwYAAMj7aCwFVHJymkaP3qGIiJ9ktxsKCCiqefNa6Nlny5odDQAA4K5RZgugX3/9Ux07RmnPnht3X+vYsbI+/vhZ+fh4mJwMAAAgeyizBYjdbmjq1L0aMmSrrl9Pla+vh2bMaKqXXqpgdjQAAIB7QpktIM6cuaLOnb/Vxo2nJEnNmgVqzpzm8vf3NDkZAADAvaPMFgCLF/+q3r036PLlJBUq5KJJk55Rz57Vsr4BAgAAgEVQZvOxixcT1bv3Bi1deliSVKdOSS1c2FLly/uanAwAACBnUGbzqXXrTqhLlzWKjr4qZ2ebRo6sp+HDn5SLCzdAAAAA+QdlNp+5di1FQ4Zs1dSpeyVJFSr4auHCFqpdu5TJyQAAAHIeZTYf2bXrD4WEfKvDhy9Kkvr2raHx4xuqcGFXk5MBAADkDspsPpCaatfYsT9q9OgdSksz5O/vqblzm+u55wLNjgYAAJCrKLMWd+TIRYWERGnnznOSpPbtK2j69GD5+hYyORkAAEDuo8xalGEYmjHjPxo4cIsSE1Pl4+Ou6dOD1aFDJbOjAQAAOAxl1oKio6+qW7c1WrPmhCSpSZOymjevhcqUKWpuMAAAAAejzFrM8uWH1bPnel28eF0eHi4aP76h+vatIScnboAAAAAKHsqsRVy+fF39+m3SokUHJUlPPOGnhQtbqHLlYiYnAwAAMA9l1gI2bTqlzp2/1enTV+TkZNOwYXU1cmQ9ubk5mx0NAADAVJRZRzh0b1MArie5aPjkFpo8v6EkKahsrBaOW6J6Nd6Wfs/JgAAAANZEmc2j9h70V8fBHXTwWElJUo/2OzRx0Cp5Fkk2ORkAAEDeQZnNY9LSbIqc/YzCpz6nlBQX+RW7otnvL1OrZw6ZHQ0AACDPoczmIb+f9lWnoa9o+56HJEn/CD6gT0f9S8V9E0xOBgAAkDdRZvMAw5Bmr6ij/uPaKOGau4oWua6P31mpTs/vlo0VtwAAALJEmTVZTKynXh/5or7Z/JgkqWGtY5o/bqkCS18yORkAAEDeR5k10coNj+mN8Bd14aKn3FxTNab/Gg0I3SpnZ8PsaAAAAJZAmTVB/FV3DYhoozlf1pEkVa0QrUWRi1Wl/DmTkwEAAFgLZdbBtv38kDoNfUUnzvrKZrNrcLctGtVvndzd0syOBgAAYDmUWQdJSnbWyI+aacKcRjIMJwWWvqgF45aoQa3jZkcDAACwLMqsAxw4UlIdB3fQ/sP+kqSu7XZq8rCv5eWZZHIyAAAAa6PM5rIPP/xZw4a+peQUFxX3vaqZo1fo+Sa/mB0LAAAgX6DM5rKTJ+OVnOKi1o1/0czRK+RX7KrZkQAAAPINymwuGzeugeoGDFSHVvu4AQIAAEAOo8zmskKFXPXq/+0zOwYAAEC+5GR2AAAAAOBeUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJaVJ8rstGnTFBgYKA8PD9WtW1c7d+68q/2WLFkim82mtm3b5m5AAAAA5Emml9mlS5cqLCxM4eHh2rNnj6pVq6ZmzZrp/Pnzt93vxIkTevvtt9WgQQMHJQUAAEBeY3qZ/fDDD/X666+rS5cuqly5smbMmKHChQtrzpw5We6Tlpam1157TaNGjdLDDz/swLQAAADIS0wts8nJydq9e7eCg4PTtzk5OSk4OFg7duzIcr/Ro0erRIkS6tat2x2PkZSUpPj4+AwPAAAA5A+mltnY2FilpaXJz88vw3Y/Pz+dO3cu032+//57zZ49WzNnzryrY0RERMjb2zv9ERAQcN+5AQAAkDeYPs0gO65cuaKQkBDNnDlTxYoVu6t9hg0bpri4uPTH6dOnczklAAAAHMXFzIMXK1ZMzs7OiomJybA9JiZGJUuWvGX8sWPHdOLECbVu3Tp9m91ulyS5uLjo8OHDCgoKyrCPu7u73N3dcyE9AAAAzGbqlVk3NzfVrFlTGzduTN9mt9u1ceNG1atX75bxFStW1IEDB7Rv3770R5s2bdS4cWPt27ePKQQAAAAFjKlXZiUpLCxMoaGhqlWrlurUqaMpU6YoISFBXbp0kSR16tRJpUuXVkREhDw8PPT4449n2N/Hx0eSbtkOAACA/M/0Mtu+fXtduHBBI0eO1Llz51S9enWtWbMm/U1hp06dkpOTpab2AgAAwEFshmEYZodwpPj4eHl7eysuLk5eXl6OOeghm2OOAwAAkFsqOq4yZqevcckTAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYVp4os9OmTVNgYKA8PDxUt25d7dy5M8uxM2fOVIMGDfTAAw/ogQceUHBw8G3HAwAAIP8yvcwuXbpUYWFhCg8P1549e1StWjU1a9ZM58+fz3T8li1b1KFDB23evFk7duxQQECAnnvuOZ09e9bByQEAAGA2m2EYhpkB6tatq9q1a2vq1KmSJLvdroCAAPXr109Dhw694/5paWl64IEHNHXqVHXq1OmWzyclJSkpKSn94/j4eAUEBCguLk5eXl4590Ju55DNMccBAADILRUdVxnj4+Pl7e19V33N1CuzycnJ2r17t4KDg9O3OTk5KTg4WDt27Lir57h27ZpSUlLk6+ub6ecjIiLk7e2d/ggICMiR7AAAADCfqWU2NjZWaWlp8vPzy7Ddz89P586du6vnGDJkiPz9/TMU4r8aNmyY4uLi0h+nT5++79wAAADIG1zMDnA/xo0bpyVLlmjLli3y8PDIdIy7u7vc3d0dnAwAAACOYGqZLVasmJydnRUTE5Nhe0xMjEqWLHnbfSdOnKhx48Zpw4YNqlq1am7GBAAAQB5l6jQDNzc31axZUxs3bkzfZrfbtXHjRtWrVy/L/SIjI/X+++9rzZo1qlWrliOiAgAAIA8yfZpBWFiYQkNDVatWLdWpU0dTpkxRQkKCunTpIknq1KmTSpcurYiICEnS+PHjNXLkSH3xxRcKDAxMn1vr6ekpT09P014HAAAAHM/0Mtu+fXtduHBBI0eO1Llz51S9enWtWbMm/U1hp06dkpPT/y4gf/LJJ0pOTtaLL76Y4XnCw8P13nvvOTI6AAAATGb6OrOOlp11y3IM68wCAACrY51ZAAAAIGdRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYFmUWAAAAlkWZBQAAgGVRZgEAAGBZlFkAAABYVp4os9OmTVNgYKA8PDxUt25d7dy587bjly9frooVK8rDw0NVqlRRVFSUg5ICAAAgLzG9zC5dulRhYWEKDw/Xnj17VK1aNTVr1kznz5/PdPwPP/ygDh06qFu3btq7d6/atm2rtm3b6r///a+DkwMAAMBsNsMwDDMD1K1bV7Vr19bUqVMlSXa7XQEBAerXr5+GDh16y/j27dsrISFBq1atSt/25JNPqnr16poxY8YdjxcfHy9vb2/FxcXJy8sr517I7RyyOeY4AAAAuaWi4ypjdvqai4MyZSo5OVm7d+/WsGHD0rc5OTkpODhYO3bsyHSfHTt2KCwsLMO2Zs2aaeXKlZmOT0pKUlJSUvrHcXFxkm58kRzmquMOBQAAkCsc2J1u9rS7ueZqapmNjY1VWlqa/Pz8Mmz38/PToUOHMt3n3LlzmY4/d+5cpuMjIiI0atSoW7YHBATcY2oAAICCyNvhR7xy5Yq8vW9/XFPLrCMMGzYsw5Vcu92uixcv6sEHH5TNlvt//o+Pj1dAQIBOnz7tuGkNyFGcQ+vjHFof59DaOH/W5+hzaBiGrly5In9//zuONbXMFitWTM7OzoqJicmwPSYmRiVLlsx0n5IlS2ZrvLu7u9zd3TNs8/HxuffQ98jLy4sfYIvjHFof59D6OIfWxvmzPkeewztdkb3J1NUM3NzcVLNmTW3cuDF9m91u18aNG1WvXr1M96lXr16G8ZK0fv36LMcDAAAg/zJ9mkFYWJhCQ0NVq1Yt1alTR1OmTFFCQoK6dOkiSerUqZNKly6tiIgISdJbb72lRo0aadKkSWrVqpWWLFmin3/+WZ999pmZLwMAAAAmML3Mtm/fXhcuXNDIkSN17tw5Va9eXWvWrEl/k9epU6fk5PS/C8j169fXF198oXfeeUfDhw/Xo48+qpUrV+rxxx836yXclru7u8LDw2+Z6gDr4BxaH+fQ+jiH1sb5s768fA5NX2cWAAAAuFem3wEMAAAAuFeUWQAAAFgWZRYAAACWRZkFAACAZVFmc8C0adMUGBgoDw8P1a1bVzt37rzt+OXLl6tixYry8PBQlSpVFBUV5aCkyEp2zuHMmTPVoEEDPfDAA3rggQcUHBx8x3OO3Jfdn8OblixZIpvNprZt2+ZuQNxRds/h5cuX1adPH5UqVUru7u4qX748/z01UXbP35QpU1ShQgUVKlRIAQEBGjBggK5fv+6gtPi7rVu3qnXr1vL395fNZtPKlSvvuM+WLVv0xBNPyN3dXY888ojmzZuX6zkzZeC+LFmyxHBzczPmzJlj/PLLL8brr79u+Pj4GDExMZmO3759u+Hs7GxERkYaBw8eNN555x3D1dXVOHDggIOT46bsnsNXX33VmDZtmrF3717j119/NTp37mx4e3sbZ86ccXBy3JTdc3jT8ePHjdKlSxsNGjQwnn/+eceERaayew6TkpKMWrVqGS1btjS+//574/jx48aWLVuMffv2OTg5DCP75+/zzz833N3djc8//9w4fvy4sXbtWqNUqVLGgAEDHJwcN0VFRRkjRowwvvzyS0OS8dVXX912/O+//24ULlzYCAsLMw4ePGh8/PHHhrOzs7FmzRrHBP4Lyux9qlOnjtGnT5/0j9PS0gx/f38jIiIi0/Evv/yy0apVqwzb6tata/To0SNXcyJr2T2Hf5eammoULVrUmD9/fm5FxB3cyzlMTU016tevb8yaNcsIDQ2lzJosu+fwk08+MR5++GEjOTnZURFxG9k9f3369DGeffbZDNvCwsKMp556Kldz4u7cTZkdPHiw8dhjj2XY1r59e6NZs2a5mCxzTDO4D8nJydq9e7eCg4PTtzk5OSk4OFg7duzIdJ8dO3ZkGC9JzZo1y3I8cte9nMO/u3btmlJSUuTr65tbMXEb93oOR48erRIlSqhbt26OiInbuJdz+PXXX6tevXrq06eP/Pz89Pjjj2vs2LFKS0tzVGz8f/dy/urXr6/du3enT0X4/fffFRUVpZYtWzokM+5fXuozpt8BzMpiY2OVlpaWfreym/z8/HTo0KFM9zl37lym48+dO5drOZG1ezmHfzdkyBD5+/vf8kMNx7iXc/j9999r9uzZ2rdvnwMS4k7u5Rz+/vvv2rRpk1577TVFRUXp6NGj6t27t1JSUhQeHu6I2Pj/7uX8vfrqq4qNjdXTTz8twzCUmpqqnj17avjw4Y6IjByQVZ+Jj49XYmKiChUq5LAsXJkF7sO4ceO0ZMkSffXVV/Lw8DA7Du7ClStXFBISopkzZ6pYsWJmx8E9stvtKlGihD777DPVrFlT7du314gRIzRjxgyzo+EubNmyRWPHjtX06dO1Z88effnll1q9erXef/99s6PBgrgyex+KFSsmZ2dnxcTEZNgeExOjkiVLZrpPyZIlszUeuetezuFNEydO1Lhx47RhwwZVrVo1N2PiNrJ7Do8dO6YTJ06odevW6dvsdrskycXFRYcPH1ZQUFDuhkYG9/JzWKpUKbm6usrZ2Tl9W6VKlXTu3DklJyfLzc0tVzPjf+7l/L377rsKCQlR9+7dJUlVqlRRQkKC3njjDY0YMUJOTlxry+uy6jNeXl4OvSorcWX2vri5ualmzZrauHFj+ja73a6NGzeqXr16me5Tr169DOMlaf369VmOR+66l3MoSZGRkXr//fe1Zs0a1apVyxFRkYXsnsOKFSvqwIED2rdvX/qjTZs2aty4sfbt26eAgABHxofu7efwqaee0tGjR9P/R0SSjhw5olKlSlFkHexezt+1a9duKaw3/8fEMIzcC4sck6f6jMPfcpbPLFmyxHB3dzfmzZtnHDx40HjjjTcMHx8f49y5c4ZhGEZISIgxdOjQ9PHbt283XFxcjIkTJxq//vqrER4eztJcJsvuORw3bpzh5uZmrFixwvjjjz/SH1euXDHrJRR42T2Hf8dqBubL7jk8deqUUbRoUaNv377G4cOHjVWrVhklSpQwPvjgA7NeQoGW3fMXHh5uFC1a1Fi8eLHx+++/G+vWrTOCgoKMl19+2ayXUOBduXLF2Lt3r7F3715DkvHhhx8ae/fuNU6ePGkYhmEMHTrUCAkJSR9/c2muQYMGGb/++qsxbdo0luayso8//tgoW7as4ebmZtSpU8f48ccf0z/XqFEjIzQ0NMP4ZcuWGeXLlzfc3NyMxx57zFi9erWDE+PvsnMOy5UrZ0i65REeHu744EiX3Z/Dv6LM5g3ZPYc//PCDUbduXcPd3d14+OGHjTFjxhipqakOTo2bsnP+UlJSjPfee88ICgoyPDw8jICAAKN3797GpUuXHB8chmEYxubNmzP93XbzvIWGhhqNGjW6ZZ/q1asbbm5uxsMPP2zMnTvX4bkNwzBshsH1fAAAAFgTc2YBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAAIBlUWYBAABgWZRZAAAAWBZlFgAAAJZFmQUAi7pw4YJKliypsWPHpm/74Ycf5Obmpo0bN5qYDAAcx2YYhmF2CADAvYmKilLbtm31ww8/qEKFCqpevbqef/55ffjhh2ZHAwCHoMwCgMX16dNHGzZsUK1atXTgwAHt2rVL7u7uZscCAIegzAKAxSUmJurxxx/X6dOntXv3blWpUsXsSADgMMyZBQCLO3bsmKKjo2W323XixAmz4wCAQ3FlFgAsLDk5WXXq1FH16tVVoUIFTZkyRQcOHFCJEiXMjgYADkGZBQALGzRokFasWKH//Oc/8vT0VKNGjeTt7a1Vq1aZHQ0AHIJpBgBgUVu2bNGUKVO0cOFCeXl5ycnJSQsXLtS2bdv0ySefmB0PAByCK7MAAACwLK7MAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAsizILAAAAy6LMAgAAwLIoswAAALAsyiwAAAAs6/8B0c247TOa4rcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## define density function \n",
    "def f2(x): \n",
    "    res = x + 1/2\n",
    "    return res\n",
    "\n",
    "f2 = np.vectorize(f2)\n",
    "\n",
    "## create histogram and fit density function \n",
    "\n",
    "fig, ax = plt.subplots(figsize  = (8, 8))\n",
    "XX = np.array(XX, dtype = float)\n",
    "m1 = XX.min(); m2 = XX.max()\n",
    "xs = np.linspace(m1, m2, 100000)\n",
    "ys = f2(xs)\n",
    "ys = np.array(ys, dtype = float)\n",
    "ax.hist(XX, color = \"gold\", density = True)\n",
    "ax.plot(xs, ys, color = \"darkblue\")\n",
    "ax.set_title(\"Histogram with fit density function\")\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"y\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}