machine-learning/CohortFitDistribution.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "title: Fit Distribution\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dạng hàm phân phối\n",
    "\n",
    "Phân phối bao gồm: Cauchy, Gamma, Exponential, Log-Logistic, Log-Normal Hàm mât đô xác suất (PDF) của phân phối Cauchy:\n",
    "\n",
    "$$\n",
    "f\\left(x ; x_{0} ; \\gamma\\right)=\\frac{1}{\\pi \\gamma\\left[1+\\left(\\frac{x-x_{0}}{\\gamma}\\right)^{2}\\right]}=\\frac{1}{\\pi \\gamma}\\left[\\frac{\\gamma^{2}}{\\left(x-x_{0}\\right)^{2}+\\gamma^{2}}\\right]\n",
    "$$\n",
    "\n",
    "Trong đó:\n",
    "\n",
    "- $x_{0}$ : là thông số vị trí, chỉ định vị trí đỉnh của phân phối\n",
    "- $\\gamma$ : là thông số tỷ lệ chỉ định nửa chiều rộng ở nửa tối đa (HWHM), $2 \\gamma$ là toàn bộ chiều rộng ở mức tối đa một nửa (FWHM). $\\gamma$ là một nửa phạm vi liên phần và đôi khi được gọi là lỗi có thể xảy ra\n",
    "\n",
    "Hàm mât đô xác suất (PDF) của phân phối Gamma:\n",
    "\n",
    "$$\n",
    "f(x ; k ; \\theta)=\\frac{x^{k-1} e^{-\\frac{x}{\\theta}}}{\\theta^{k} \\Gamma(k)}, \\text { for } x>0 \\text { and } k, \\theta>0\n",
    "$$\n",
    "\n",
    "- $k$ : là tham số hình dạng\n",
    "- $\\theta$ : là tham số tỷ lệ\n",
    "- $x$ : là biến ngẫu nhiên\n",
    "- $\\Gamma(k)$ : là hàm gamma được đánh giá tại $\\mathrm{k}$\n",
    "\n",
    "Hàm mật đô xác suất (PDF) của phân phối Exponential:\n",
    "\n",
    "$$\n",
    "f(x ; \\lambda)=\\left\\{\\begin{array}{rr}\n",
    "0, & x<0 \\\\\n",
    "\\lambda e^{-\\lambda x}, & x \\geq 0\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "- $\\lambda$ : là tham số của phân phối, thường được gọi là tham số tỷ lệ.\n",
    "- $x$ : là biến ngẫu nhiên\n",
    "\n",
    "Hàm mật độ xác suất (PDF) của phân phối Log-Logistic:\n",
    "\n",
    "$$\n",
    "f(x ; \\alpha ; \\beta)=\\frac{\\left(\\frac{\\beta}{\\alpha}\\right)\\left(\\frac{x}{\\alpha}\\right)^{\\beta-1}}{\\left(1+\\left(\\frac{x}{\\alpha}\\right)^{\\beta}\\right)^{2}}, \\text { where } x>0, \\alpha>0, \\beta>0\n",
    "$$\n",
    "\n",
    "- $\\alpha$ : là tham số tỷ lệ và là giá trị trung bình của phân phối\n",
    "- $\\beta$ : là tham số hình dạng\n",
    "\n",
    "Hàm mât độ xác suất (PDF) của phân phối Lognormal $\\left(\\mu, \\sigma^{2}\\right)$ :\n",
    "\n",
    "$$\n",
    "f_{x}(x)=\\frac{d}{d x} \\operatorname{Pr}(X \\leq x)=\\frac{1}{x \\sigma \\sqrt{2 \\pi}} e^{\\left(-\\frac{(\\ln x-\\mu)^{2}}{2 \\sigma^{2}}\\right)}, \\text { where } x>0, \\sigma>0\n",
    "$$\n",
    "\n",
    "- $x$ là biến ngẫu nhiên \n",
    "- Một biến ngẫu nhiên $X$ tuân theo phân phối Log-normal $\\left(X \\sim \\operatorname{Lognormal}\\left(\\mu_{x}, \\sigma_{x}^{2}\\right)\\right)$ nếu $\\ln (X)$ tuân theo phân phối chuẩn với giá trị trung bình là $\\mu$ và phương sai $\\sigma^{2}$\n",
    "\n",
    "Để lựa chọn phân phối phù hợp cho từng phân nhóm, Nhóm mô hình lựa chọn mô hình có SSE (Sum of square error - tổng của sai số bình phương) nhỏ nhất với SSE được tính theo công thức sau:\n",
    "\n",
    "$$\n",
    "\\left.S S E=\\sum \\text { (Giá trị thực tế - Giá trị được tính ra từ mô hình phân phối }\\right)^{2}\n",
    "$$\n",
    "\n",
    "Cuối cùng, từ mô hình phân phối vừa được lựa chọn (là mô hình có SSE nhỏ nhất)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAEiCAYAAAAPogpgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACI30lEQVR4nO3deXhTZdoG8PskaZLupXuBUgqyFFktW2EAFSkDDuOoDIyMgAqjTHHBunxWHAV0RFCxjrLIiDJutSq4jYxQRTbBhVpwAZG9LC2lLd3bpEnO90dyTptuNG1O0yT377pyYU9P0rc1b3KePM/7vIIoiiKIiIiIiIiISBEqVw+AiIiIiIiIyJMx8CYiIiIiIiJSEANvIiIiIiIiIgUx8CYiIiIiIiJSEANvIiIiIiIiIgUx8CYiIiIiIiJSEANvIiIiIiIiIgUx8CYiIiIiIiJSEANvIiIiIiIiIgUx8O4kfvzxR9x+++2Ij4+HXq9HQEAArrrqKqxcuRLFxcUuHZsgCLj77rtdOgbyfD/++CPmzZuH3r17w9fXF76+vujTpw/uuusu7N+/39XDI7qsjRs3QhCEZm87duxw9RBdpqqqCkuWLGnybyD93U6dOuXw47bnvqQs6f9NZ3j9XrJkCQRBQGRkJMrLyxt9v2fPnvjDH/7ggpE5x9VXX42rr77a1cOgJnTGeVBYWNihP7dnz5647bbbHLrP3r17sWTJEpSUlDT6njs/3zWuHgAB//73v5GSkoJ+/frhoYcewoABA1BbW4v9+/dj3bp12LdvHz788ENXD5NIMa+88gruvvtu9OvXD/fddx+uvPJKCIKAw4cPIyMjAyNGjMCxY8fQu3dvVw+V6LJef/119O/fv9HxAQMGuGA0nUNVVRWWLl0KAI0umK6//nrs27cPMTExLhgZeZOLFy9i5cqVePLJJ109FCKv8eGHHyIoKMih++zduxdLly7FbbfdhpCQELvvrVmzxomj61gMvF1s3759+Pvf/45Jkybho48+gk6nk783adIkPPDAA/j8889dOEIiZX399ddISUnB9ddfjw8++ABarVb+3rXXXouFCxfi/fffh6+vrwtHSdR6AwcOxPDhw109DLcRERGBiIgIVw+DvMDvf/97vPDCC1i4cCGio6MV+RmiKKKmpobvWUQ2w4YNc+rjufOH2Cw1d7Gnn34agiBg/fr1dkG3RKvV4o9//CMAIDMzE8nJyYiJiYGvry8SEhLwyCOPoLKy0u4+zZVg3HbbbejZs6fdMYPBgGXLliEhIQF6vR5hYWG45pprsHfv3kb3f/PNN5GQkAA/Pz8MGTIE//3vf+Xv7d69G4IgICMjo9H93njjDQiCgO+//741fxLyMk8//TTUajVeeeUVu6C7vj//+c/o2rUrAGD//v34y1/+gp49e8LX1xc9e/bELbfcgtOnT9vdRyrv2r59O/72t78hLCwMQUFBmDNnDiorK5Gfn48ZM2YgJCQEMTExePDBB1FbWyvf/9SpUxAEAc8++yxWrFgh/7yrr74av/32G2pra/HII4+ga9euCA4Oxo033oiCggK7MbR2zpL3ePfddyEIAl5++WW740888QTUajWysrIA1D3/Vq5ciX/+85/o0aMH9Ho9hg8fji+//LLR4+7ZswcTJ05EYGAg/Pz8MGbMGHz22Wd250hz4quvvsLf//53hIeHIywsDDfddBPOnz/f6DEzMzORlJQEf39/BAQEYPLkycjJybE757bbbkNAQACOHTuGqVOnIiAgALGxsXjggQdgMBjk30UKrJcuXSqX3kulh02Vi2dlZeGGG25A9+7dodfrccUVV+Cuu+7q8BJJUl5rnrvSeUlJSdDr9ejWrRv+8Y9/4NVXX3VoqcFTTz0Fk8mEJUuWXPbc4uJipKSkoFu3btBqtejVqxcWL14sP68l0nK8devWISEhATqdDv/5z3/a/R4EWOfLqFGjEBoaiqCgIFx11VXYsGEDRFFs1e9L7qMj58HlfPLJJ0hKSoKfnx8CAwMxadIk7Nu3r9F5H3/8MQYPHgydTodevXrhxRdflMvZ62tYam6xWPDUU0+hX79+8PX1RUhICAYPHowXX3wRgLUk/qGHHgIAxMfHN1qu1VSc40g840rMeLuQ2WzG9u3bkZiYiNjY2Muef/ToUUydOhWLFi2Cv78/fv31V6xYsQLfffcdtm/f7vDPN5lMmDJlCnbv3o1Fixbh2muvhclkwjfffIPc3FyMGTNGPvezzz7D999/j2XLliEgIAArV67EjTfeiCNHjqBXr14YN24chg0bhtWrV+OWW26x+zkvv/wyRowYgREjRjg8RvJsZrMZX331FYYPH97qMtNTp06hX79++Mtf/oLQ0FDk5eVh7dq1GDFiBA4dOoTw8HC78+fPn4+bbroJ7777LnJycvDoo4/CZDLhyJEjuOmmm3DnnXfiiy++wIoVK9C1a1ekpqba3X/16tUYPHgwVq9ejZKSEjzwwAOYNm0aRo0aBR8fH7z22ms4ffo0HnzwQcyfPx+ffPKJfF9nz1lyD2azGSaTye6YIAhQq9X4y1/+gp07d+KBBx7A6NGjMXz4cGzfvh1PPfUUHn30UUyaNMnufi+//DLi4uKQnp4Oi8WClStXYsqUKdi5cyeSkpIAADt37sSkSZMwePBgbNiwATqdDmvWrMG0adOQkZGBmTNn2j3m/Pnzcf311+Odd97BmTNn8NBDD+HWW2+1e04+/fTTeOyxx3D77bfjscceg9FoxLPPPotx48bhu+++s8s41NbW4o9//CPmzZuHBx54ALt27cKTTz6J4OBgPP7444iJicHnn3+O3//+95g3bx7mz58PAC1muY8fP46kpCTMnz8fwcHBOHXqFFatWoXf/e53+Omnn+Dj49O2/znUqbT2ufvjjz9i0qRJ6Nu3L/7zn//Az88P69atw1tvveXQz4uLi0NKSgpeeuklpKamom/fvk2eV1NTg2uuuQbHjx/H0qVLMXjwYOzevRvLly/HgQMHGgVEH330EXbv3o3HH38c0dHRiIyMlJMN7XkPOnXqFO666y706NEDAPDNN9/gnnvuwblz5/D444879LtT59XR86Al77zzDv76178iOTkZGRkZMBgMWLlyJa6++mp8+eWX+N3vfgcA+Pzzz3HTTTdh/PjxyMzMhMlkwnPPPYcLFy5c9mesXLkSS5YswWOPPYbx48ejtrYWv/76q7yee/78+SguLsZLL72EzZs3y9eHzWW6HYlnXE4kl8nPzxcBiH/5y18cvq/FYhFra2vFnTt3igDEgwcPyt+bMGGCOGHChEb3mTt3rhgXFyd//cYbb4gAxH//+98t/iwAYlRUlFhWVmY3dpVKJS5fvlw+9vrrr4sAxJycHPnYd999JwIQ//Of/zj8O5Lna2kOmEwmsba2Vr5ZLJYmH8NkMokVFRWiv7+/+OKLL8rHpefjPffcY3f+n/70JxGAuGrVKrvjQ4cOFa+66ir565MnT4oAxCFDhohms1k+np6eLgIQ//jHP9rdf9GiRSIAsbS0tMlxtjRnyTNIz7mmbmq1Wj6vpqZGHDZsmBgfHy8eOnRIjIqKEidMmCCaTCb5HOn517VrV7G6ulo+XlZWJoaGhorXXXedfGz06NFiZGSkWF5eLh8zmUziwIEDxe7du8tzRxpfSkqK3bhXrlwpAhDz8vJEURTF3NxcUaPRNJo75eXlYnR0tDhjxgz52Ny5c0UA4nvvvWd37tSpU8V+/frJX1+8eFEEID7xxBPN/t1OnjzZ5N9VmjunT58WAYgff/xxq+9LriP9v/n++++bPae1z90///nPor+/v3jx4kX5PLPZLA4YMKBV//+feOIJEYB48eJFsbCwUAwODhZvvvlm+ftxcXHi9ddfL3+9bt26Jp/XK1asEAGI27Ztk48BEIODg8Xi4uImf/+2vgc1ZDabxdraWnHZsmViWFiY3Xtic9d95HqddR40xWw2i127dhUHDRpkd91TXl4uRkZGimPGjJGPjRgxQoyNjRUNBoPdeWFhYWLD8DIuLk6cO3eu/PUf/vAHcejQoS2O9dlnn232d2r4fG9tPNMZsNTcjZw4cQKzZs1CdHQ01Go1fHx8MGHCBADA4cOHHX68//3vf9Dr9bjjjjsue+4111yDwMBA+euoqChERkbalffecsstiIyMxOrVq+VjL730EiIiIhplXIguJzExET4+PvLt+eefBwBUVFTg//7v/3DFFVdAo9FAo9EgICAAlZWVTc6Dhp1qExISAFgbOjU83rBcHQCmTp0KlUpld15z9weA3Nxc+Ziz5yy5hzfeeAPff/+93e3bb7+Vv6/T6fDee++hqKgIV111FURRREZGBtRqdaPHuummm6DX6+WvAwMDMW3aNOzatQtmsxmVlZX49ttvMX36dAQEBMjnqdVqzJ49G2fPnsWRI0fsHlNaviQZPHgwAMjP/61bt8JkMmHOnDkwmUzyTa/XY8KECY06kwuCgGnTpjV6zKbmU2sVFBRgwYIFiI2NhUajgY+PD+Li4gBw7ngKR567O3fuxLXXXmtX0aRSqTBjxgyHf25YWBj+7//+D5s2bbKbl/Vt374d/v7+mD59ut1xqVy24XKPa6+9Fl26dGnysdrzHrR9+3Zcd911CA4Olt9DHn/8cRQVFTVa2kTuyVXzoClHjhzB+fPnMXv2bLvrnoCAANx888345ptvUFVVhcrKSuzfvx9/+tOf7JYIBgQENHovaMrIkSNx8OBBpKSkYOvWrSgrK2vXuB2JZ1yNpeYuFB4eDj8/P5w8efKy51ZUVGDcuHHQ6/V46qmn0LdvX/j5+eHMmTO46aabUF1d7fDPv3jxIrp27Wo3uZoTFhbW6JhOp7P7uTqdDnfddReef/55PPvss6itrcV7772H1NTUJtevE4WHh8PX17fJC/R33nkHVVVVyMvLswsUZs2ahS+//BL/+Mc/MGLECAQFBUEQBEydOrXJeRAaGmr3tfQm0dTxmpqadt0fgPwYSsxZcg8JCQmXba52xRVXYNy4cfjss8/w97//vdmlFk01gIqOjobRaERFRQXKy8shimKT95f6IhQVFdkdb/h6Lr0+S89JqVSwueVBDd8z/Pz87D4ckB6zqfnUGhaLBcnJyTh//jz+8Y9/YNCgQfD394fFYsHo0aM5dzzEpUuXWv3cLSoqQlRUVKPzmjrWGosWLcLLL7+Mhx9+GDt37mz0/aKiIkRHRzdaqxoZGQmNRtNoTrW0VKqt70HfffcdkpOTcfXVV+Pf//43unfvDq1Wi48++gj//Oc/OQ88hCvnQUPSz2luLBaLRR6vKIptHktaWhr8/f3x1ltvYd26dVCr1Rg/fjxWrFjRpsakjsQzrsbA24XUajUmTpyI//3vfzh79iy6d+/e7Lnbt2/H+fPnsWPHDjljBqDJ/e30ej1KS0sbHW/YlCYiIgJ79uyBxWJx2pP173//O5555hm89tprqKmpgclkwoIFC5zy2OR51Go1rr32Wmzbtg15eXl2L/bSWp76zUJKS0vx3//+F0888QQeeeQR+bjBYHD5fvcNOTJnyfu8+uqr+OyzzzBy5Ei8/PLLmDlzJkaNGtXovPz8/CaPabVaBAQEQKPRQKVSIS8vr9F5UsO0hn0PLkc6/4MPPpCzzB3p559/xsGDB7Fx40bMnTtXPn7s2LEOHwspp0uXLq1+7oaFhTW5drSp+dEavr6+WLJkCe68884mG1iFhYXh22+/hSiKdsF3QUEBTCZToznVMEB3hnfffRc+Pj7473//a/fB1kcffeT0n0Wu48p50JD0oWxzY1GpVOjSpYs8L9o6Fo1Gg9TUVKSmpqKkpARffPEFHn30UUyePBlnzpyBn5+fQ+NWIp5RSucenRdIS0uDKIr429/+BqPR2Oj7tbW1+PTTT+UX9YaZ41deeaXRfXr27InffvvNrvNmUVFRo85+U6ZMQU1NDTZu3OiE38QqJiYGf/7zn7FmzRqsW7cO06ZNk5uCEDUlLS0NZrMZCxYsaNTRtSFBECCKYqN58Oqrr8JsNis5TIc5MmfJu/z000+49957MWfOHOzevRuDBw/GzJkzcenSpUbnbt682S4LVl5ejk8//RTjxo2DWq2Gv78/Ro0ahc2bN9tlwCwWC9566y1079692QZSzZk8eTI0Gg2OHz+O4cOHN3lzVMOseks4d7yDI8/dCRMmYPv27XYJBIvFgvfff7/NP/+OO+6Qd5qwWCx235s4cSIqKioaBblvvPGG/H2lCYIAjUZjtwSluroab775puI/mzqOq+dBff369UO3bt3wzjvv2HXOr6ysxKZNm+RO5/7+/hg+fDg++ugju9iloqLCbsej1ggJCcH06dOxcOFCFBcXy8kWR94zlIhnlMKMt4slJSVh7dq1SElJQWJiIv7+97/jyiuvRG1tLXJycrB+/XoMHDgQr776Krp06YIFCxbgiSeegI+PD95++20cPHiw0WPOnj0br7zyCm699Vb87W9/Q1FREVauXNlo8/pbbrkFr7/+OhYsWIAjR47gmmuugcViwbfffouEhAT85S9/adPvdN9998mZm9dff71Nj0HeY+zYsVi9ejXuueceXHXVVbjzzjtx5ZVXyp8Ab9q0CQAQFBSEoKAgjB8/Hs8++yzCw8PRs2dP7Ny5Exs2bEBISIhrf5EGxowZ0+o5S57l559/btTVHAB69+4NPz8/zJgxA/Hx8VizZg20Wi3ee+89XHXVVbj99tsbXeir1WpMmjQJqampsFgsWLFiBcrKyrB06VL5nOXLl2PSpEm45ppr8OCDD0Kr1WLNmjX4+eefkZGR4XA2rmfPnli2bBkWL16MEydO4Pe//z26dOmCCxcu4LvvvoO/v7/dz2+NwMBAxMXF4eOPP8bEiRMRGhoqz+GG+vfvj969e+ORRx6BKIoIDQ3Fp59+Km+1Ru5l+/btTW5zNHXq1FY/dxcvXoxPP/0UEydOxOLFi+Hr64t169bJWzO2JculVqvx9NNP48YbbwRQ1+sAAObMmYPVq1dj7ty5OHXqFAYNGoQ9e/bg6aefxtSpU3Hddde14S/hmOuvvx6rVq3CrFmzcOedd6KoqAjPPfccl+65qc40Dz799FO7vk2S6dOnY+XKlfjrX/+KP/zhD7jrrrtgMBjw7LPPoqSkBM8884x87rJly3D99ddj8uTJuO+++2A2m/Hss88iICDgshWI06ZNw8CBAzF8+HBERETg9OnTSE9PR1xcHPr06QMAGDRoEADgxRdfxNy5c+Hj44N+/fo1OW6l4hlFuKanGzV04MABce7cuWKPHj1ErVYr+vv7i8OGDRMff/xxsaCgQBRFUdy7d6+YlJQk+vn5iREREeL8+fPFH374QQQgvv7663aP95///EdMSEgQ9Xq9OGDAADEzM7NRV3NRFMXq6mrx8ccfF/v06SNqtVoxLCxMvPbaa8W9e/fK5wAQFy5c2GjMDbsU1tezZ08xISGhXX8T8i4HDhwQb7/9djE+Pl7U6XSiXq8Xr7jiCnHOnDnil19+KZ939uxZ8eabbxa7dOkiBgYGir///e/Fn3/+udHzsblOos119Zw7d67o7+8vfy11lX722Wftzvvqq69EAOL7779vd7ypn+fInCX311JXc9g6rt56662in5+f+Msvv9jd9/333xcBiC+88IIoinXPvxUrVohLly4Vu3fvLmq1WnHYsGHi1q1bG/3s3bt3i9dee63o7+8v+vr6iqNHjxY//fTTJsfXcE5Iz+mvvvrK7vhHH30kXnPNNWJQUJCo0+nEuLg4cfr06eIXX3whn9Nw3kikeVbfF198IQ4bNkzU6XQiAHm+NtWZ/NChQ+KkSZPEwMBAsUuXLuKf//xnMTc3t1FndHY177wuNx+k/2etee5K540aNUrU6XRidHS0+NBDD8ldxktKSlocS0vdnMeMGSMCsOtqLoqiWFRUJC5YsECMiYkRNRqNGBcXJ6alpYk1NTV25zV3jdTe9yBRFMXXXntN7Nevn6jT6cRevXqJy5cvFzds2NDoOc+u5p1XZ5wHzd0kH330kThq1ChRr9eL/v7+4sSJE8Wvv/660eN9+OGH4qBBg0StViv26NFDfOaZZ8R7771X7NKli915Da/Pnn/+eXHMmDFieHi4fN958+aJp06dsrtfWlqa2LVrV1GlUtm9RzX1fG9NPNMZCKJYr5aAyAl+/PFHDBkyBKtXr0ZKSoqrh0NE5HZOnTqF+Ph4PPvss3jwwQddPRyiTik5ORmnTp3Cb7/95uqhELlMZ5kHtbW1GDp0KLp164Zt27a5dCydFUvNyWmOHz+O06dP49FHH0VMTIy87QYRERFRe6SmpmLYsGGIjY1FcXEx3n77bWRlZWHDhg2uHhpRh+lM82DevHmYNGkSYmJikJ+fj3Xr1uHw4cN48cUXO3ws7oKBNznNk08+iTfffBMJCQl4//33He5KSERERNQUs9mMxx9/HPn5+RAEAQMGDMCbb76JW2+91dVDI+ownWkelJeX48EHH8TFixfh4+ODq666Clu2bOmQHgjuiqXmRERERERERAridmJERERERERECmLgTURERERERKQgBt5ERERERERECvKY5moWiwXnz59HYGCgvNE8UVuJoojy8nJ07doVKpV7fT7FuUDOxLlAVMdd5wPnAjmbu84FgPOBnMuRueAxgff58+cRGxvr6mGQhzlz5gy6d+/u6mE4hHOBlMC5QFTH3eYD5wIpxd3mAsD5QMpozVzwmMA7MDAQgPWXDgoKcvFoyN2VlZUhNjZWfl65E84FcibOBaI67jofOBfI2dx1LgCcD+RcjswFjwm8pVKRoKAgTiJyGncsQeJcICVwLhDVcbf5wLlASnG3uQBwPpAyWjMX3GtRBhEREREREZGbYeBNREREREREpCAG3kREREREREQKYuBNRESK27VrF6ZNm4auXbtCEAR89NFHl73Pzp07kZiYCL1ej169emHdunXKD5SIiIhIAQy8iYhIcZWVlRgyZAhefvnlVp1/8uRJTJ06FePGjUNOTg4effRR3Hvvvdi0aZPCIyUiIiJyPo/pak5ERJ3XlClTMGXKlFafv27dOvTo0QPp6ekAgISEBOzfvx/PPfccbr75ZoVGSURERKQMj89455VWI23zTyitrnX1UIg6vQqDCZuyz+JYQYWrh0Jebt++fUhOTrY7NnnyZOzfvx+1te17PS+tqsV/9p7C61+fbNfjEHkqk9mCjw+cQ0FZjauHQtRhvjx8Aa/uPoHDeWWuHgp5KI8PvNM2/4SM73KR8na2q4dC1KmJooi73tyPB94/iD+8tBu/5vONh1wnPz8fUVFRdseioqJgMplQWFjY5H0MBgPKysrsbk0pra7FE5/8gme3HnH6uIk8wetfn8J97x7A397ktRMpZ/ny5RgxYgQCAwMRGRmJP/3pTzhy5PKvy0r1/3hv/xk89dlh7D99ySmPR9SQxwfeKVdfAZUAfH2sCGeKq1w9HKJO64fcEnx9rAgAUFNrwQtZv7l4ROTtBEGw+1oUxSaPS5YvX47g4GD5Fhsb2+R5fjo1AKDKaIbFIjpxxESe4b8/5QEADp4pQXGl0cWjIU+1c+dOLFy4EN988w2ysrJgMpmQnJyMysrKZu+jZP8PAU2/txA5i8cH3iPjQ3FVjy4AgL3Hm86SEBGwxXahlRATBADY/msBLvGCi1wkOjoa+fn5dscKCgqg0WgQFhbW5H3S0tJQWloq386cOdPkeQG6uvYm1bVm5w2ayEMUVxrk/z5xkUuPSBmff/45brvtNlx55ZUYMmQIXn/9deTm5iI7u/lKi/r9PxISEjB//nzccccdeO6555w3MJEfyJIyPD7wBoDEOGvgfeg8S2eJmrP/VDEA4O9X90b/6EDUmkXs/O2ii0dF3iopKQlZWVl2x7Zt24bhw4fDx8enyfvodDoEBQXZ3Zo8T6OCypbYqDSanDpuIncniiIulNYF3qeLWC1IHaO0tBQAEBoa2uw5Svb/kIqpGHaTUrwi8O4TFQgA+O0CP7Ulakqt2YLD+eUAgMHdgjG+bwQAVomQ81RUVODAgQM4cOAAAGu54IEDB5CbmwvAmq2eM2eOfP6CBQtw+vRppKam4vDhw3jttdewYcMGPPjgg+0eiyAI8Ndas95VBma8ieorrjTCaLbIX5/mMj3qAKIoIjU1Fb/73e8wcODAZs9Tsv+HHHgz8iaFeEXg3TcqAABwtKDcxSMh6px+u1AOo8mCQL0GcWF+GNPbWsq793iRi0dGnmL//v0YNmwYhg0bBgBITU3FsGHD8PjjjwMA8vLy5CAcAOLj47Flyxbs2LEDQ4cOxZNPPol//etfTttKTFrnzYw3kb2LFQb7r8vZ2ZyUd/fdd+PHH39ERkbGZc9Vqv+HtMZbZORNCvGKfbx7RVgD78IKIyoMJrv1fUQEHLFluwfEBEEQBHl5xtlL1SisMCA8QOfK4ZEHuPrqq1u8mNm4cWOjYxMmTMAPP/ygyHisGW8DKpnxJrJTXmP/YdTFcvb6IGXdc889+OSTT7Br1y507969xXPb2v8jNTVV/rqsrKzp4Ju91UhhXpHxDtBpEKS3Btt5JdUuHg1R55NrKyXsGeYPAAjU+6BXhPW/fzpX6rJxESmFGW+iplUY7OdEYYMMOJGziKKIu+++G5s3b8b27dsRHx9/2fso2f9DHlfrfwUih7Qp8F6zZg3i4+Oh1+uRmJiI3bt3N3vunj17MHbsWISFhcHX1xf9+/fHCy+80Oi8TZs2YcCAAdDpdBgwYAA+/PDDtgytWV1DfAEA5xh4EzUiBd49wvzkY4O7BQMAfj7LwJs8jx/XeBM1qdIWePuorek/Bt6klIULF+Ktt97CO++8g8DAQOTn5yM/Px/V1XXX6h3a/8P2LyvNSSkOB96ZmZlYtGgRFi9ejJycHIwbNw5TpkyxW5tXn7+/P+6++27s2rULhw8fxmOPPYbHHnsM69evl8/Zt28fZs6cidmzZ+PgwYOYPXs2ZsyYgW+//bbtv1kD3WyB9/kSrlUiauhssfVNLja0LvAeaAu8f2TGmzyQtOSIGW8iexW2UnOpAqqwwsA1r6SItWvXorS0FFdffTViYmLkW2ZmpnxOR/b/kNaI89lOSnF4sfOqVaswb948zJ8/HwCQnp6OrVu3Yu3atVi+fHmj8+s30wGAnj17YvPmzdi9ezfuvPNO+TEmTZqEtLQ0ANZPt3bu3In09PRWNVloja5y4M2MN1FDcsa7XuA9oKu1FOu3C2xKSJ7HT2stNa8yMPAmqk8qNY8L88fRggrU1FpQaTSzPw45XWs+0OnI/h91GW+G3qQMhzLeRqMR2dnZjfbPS05Oxt69e1v1GDk5Odi7dy8mTJggH2tuT76WHrO1WwNIYkL0ABh4EzVUa7bggq1rbfcuvvLxPpHWbfhyi6tQU8tyXPIs0nZilUY+t4nqkwLviEAdtGrrZWJJFRuskedrpik6kdM4FHgXFhbCbDY3uX9eww6DDXXv3h06nQ7Dhw/HwoUL5Yw50PyefC09Zmu3BpAfL9AaeBeUc60SUX1FFUaIIqBWCQj108rHwwO0CPHzgSgCxwoqXDhCIueTm6sx401kR5oTgXoNgv2szapKq2tdOSQiIo/QpuZqTe2f19zeeZLdu3dj//79WLduXZMl5I4+ZlpaGkpLS+XbmTNnWvz54YHW7ZDYJITI3kXbh1HhAVqoVHVzThAE9LVlvRl4k6eRMt5VzHgT2amwNRwM0GkQ7GsLvKsYeJPnY3M1UppDC3bCw8OhVqub3D+vYca6IWmLgEGDBuHChQtYsmQJbrnlFgDN78nX0mPqdDrodK3fWzg8wJrJY+BNZE+aE03t1X1FVAC+O1WMowVc502ehRlvoqZJpeb+9QNvZrzJC9Q1V2PkTcpwKOOt1WqRmJjYaP+8rKwsjBkzptWPI4oiDIa6ALi5PfkceczLibAFFcWVRpgtnFBEEinjHRHYOPDuExkAADh6gRlv8ixSoyhmvInsSR9GBejUCGHgTV6EGW9SmsMtKlNTUzF79mwMHz4cSUlJWL9+PXJzc7FgwQIA1hLwc+fO4Y033gAArF69Gj169ED//v0BWPf1fu6553DPPffIj3nfffdh/PjxWLFiBW644QZ8/PHH+OKLL7Bnzx5n/I4AgFB/LQQBsIjW4LupIIPIG11sKeNtC7yPXWTgTZ7FT8vtxIiaIjXT1Puo5Yx3CQNv8gZsrkYKczjwnjlzJoqKirBs2TLk5eVh4MCB2LJlC+Li4gA03m/PYrEgLS0NJ0+ehEajQe/evfHMM8/grrvuks8ZM2YM3n33XTz22GP4xz/+gd69eyMzMxOjRo1ywq9opVGr0MVPi+JKIworDAy8iWxaynhL+7ieKa6C2SJCreK7EnkGf3k7MWa8ieozmCwAAJ1GjSBmvMkLMeFNSmnTpowpKSlISUlp8nsN99u755577LLbzZk+fTqmT5/eluG0WnhAXeBNRFYtZby7hvjCRy2g1izifEk1Yuvt803kzvxspeYVXONNZMdgsn4YpfNRIYRdzcmLCLaUN0vNSSlt6mrurqTAgoE3UZ3CFjLeapUgB9u5xVUdOi4iJckZb5aaE9kx1EoZbxW7mpNXkTZTYnM1Uop3Bt7lRhePhKjzKKq0zodwf22T35fKzU8VVXbYmIiUVrfGm6XmRPXVLzVnV3PyJmyuRkrzqsA71BZYFFcx8CaSSBdU0lq+huLCrBnv00XMeJPnkLuas9ScyI7R1ETGm4E3eQGBbWxIYV4VeHfxswbeJQy8iWTSBVVwc4G3rdT8VCEz3uQ5pH28q2rNsHCLSSKZtMZbX2+Nd0k1r5uIiNrLuwJvf+sbyKVKfnJLBFi3jZGyG8F+zQTe4dZSc67xJk/ibys1F0WgxsRycyKJVGquVasRqLe+L5TXsDKEPF9dczV+GEvK8KrAO8SW8b7EjDcRgLpst0oAArRNb3JQf40334zIU+h9VHJZITubE9WR13j7qBCot74vlNeY+PpPHk9ursanOinEuwJvrlUislNWb323qpk9uruF+EKtElBTa0FBOXcEIM8gCIKc9eZe3kRWJrMFZtvSC51GJWe8zRYRNbZu50Seqq6rOZEyvCrw7sKMN5Gdy63vBgCtRoVuIb4AuM6bPIufbUuxSm4pRgSgLtsNWLua+/mo5WCkvIZJC/J07K5GyvKqwFtqEnKpqpYlU0RoXeANAD1sDdbOXKpWfExEHUXubM4txYgA2AfeWo0KKpUgz5MyrvMmL8EQgZTiVYF3F9t2YkaTBdW1vNAiam3g3b2LNeN99hIbrJHnkDqbV3KNNxGAuq3EfNQC1LblR0G2cnP2QiBPV1dqzsiblOFVgbe/Vg0ftXVWXapiyRTR5fbwltQF3sx4k+fw0zLjTVSftJWYTqOWj0kZb5aak6eTCs2Z8SaleFXgLQhCXWfzSq7zJmp9xttaas6MN3kSf9sab2byiKzkrcQ0dZeH9TubE3kyNlcjpXlV4A0AXWzrvEuY8SZqQ6k5M97kOfxsmTyWmhNZGWydy3VNBN4VDLyJiNrF6wJv7uVNVKes2nohJa3ha46U8c4vrYHJzC1lyDMEsrkauciaNWsQHx8PvV6PxMRE7N69u8XzDQYDFi9ejLi4OOh0OvTu3Ruvvfaa08dVV2ped3kYYHt/KGOpOXk4AdzIm5SlcfUAOlpdxpuBN5GU6QvQt/xSEBmog49aQK1ZxIVyg7y9GJE7q+vWzICCOk5mZiYWLVqENWvWYOzYsXjllVcwZcoUHDp0CD169GjyPjNmzMCFCxewYcMGXHHFFSgoKIDJ5PwMtJGl5uTFWGpOSvO6jHfdXt680KKO4UhmY8+ePRg7dizCwsLg6+uL/v3744UXXlBsbNL+xQE6dYvnqVQCutqC7bPFXOdNniGAJbTkAqtWrcK8efMwf/58JCQkID09HbGxsVi7dm2T53/++efYuXMntmzZguuuuw49e/bEyJEjMWbMGKePzWSxhhwaVROl5lySQR6OzdVIaV4XeLPUnDqSlNlYvHgxcnJyMG7cOEyZMgW5ublNnu/v74+7774bu3btwuHDh/HYY4/hsccew/r16xUZn3QhJXV3bgnXeZOnkTLeDCiooxiNRmRnZyM5OdnueHJyMvbu3dvkfT755BMMHz4cK1euRLdu3dC3b188+OCDqK52/muxWQq81YJ8LJBdzclLCIJw+ZOI2sFrS81LmfGmDlA/swEA6enp2Lp1K9auXYvly5c3On/YsGEYNmyY/HXPnj2xefNm7N69G3feeafTxyeXmutaEXiH+AEoYuBNHkPqbcASWuoohYWFMJvNiIqKsjseFRWF/Pz8Ju9z4sQJ7NmzB3q9Hh9++CEKCwuRkpKC4uLiZtd5GwwGGAwG+euysrJWjU/KeEt7eANAIOcJeRnu401K8bqMdxdmvKmDtCWz0VBOTg727t2LCRMmKDFEVBqsjXT8WxN4yxlvlppT2zjaUOrtt9/GkCFD4Ofnh5iYGNx+++0oKipy2nhYak6u0jCzJopis9k2i8UCQRDw9ttvY+TIkZg6dSpWrVqFjRs3Npv1Xr58OYKDg+VbbGxsq8ZltljXeGvsAm9WhpB3Yak5KcXrAu9gqblaNTPepKy2ZDYk3bt3h06nw/Dhw7Fw4UI5Y94Ug8GAsrIyu1trtXaNNwB0D7UG3udKmPEmxzm67GLPnj2YM2cO5s2bh19++QXvv/8+vv/++xbngqOkSo9yBhTUQcLDw6FWqxu9BxQUFDR6r5DExMSgW7duCA4Olo8lJCRAFEWcPXu2yfukpaWhtLRUvp05c6ZV42sq413XhJDzhDwbm6uR0rwu8A7xZak5dSxHMhuS3bt3Y//+/Vi3bh3S09ORkZHR7LltzWwAdaXmrct4W7cUY6k5tYWjDaW++eYb9OzZE/feey/i4+Pxu9/9DnfddRf279/vtDHVdWvm+wF1DK1Wi8TERGRlZdkdz8rKarZZ2tixY3H+/HlUVFTIx3777TeoVCp07969yfvodDoEBQXZ3VrD3GRzNanUnPOEPJu0nRgz3qSUNgXejpQLbt68GZMmTUJERASCgoKQlJSErVu32p2zceNGCILQ6FZTU9OW4bVIaq7GjDcprS2ZDUl8fDwGDRqEv/3tb7j//vuxZMmSZs9ta2bDYDKj1mx9d3Gkudr5kmr54oyoNdqy7GLMmDE4e/YstmzZAlEUceHCBXzwwQe4/vrrnTYultCSK6SmpuLVV1/Fa6+9hsOHD+P+++9Hbm4uFixYAMD6mj5nzhz5/FmzZiEsLAy33347Dh06hF27duGhhx7CHXfcAV9f527taDI3tcabSzLIO7C3GinN4cDb0XLBXbt2YdKkSdiyZQuys7NxzTXXYNq0acjJybE7LygoCHl5eXY3vV7ftt+qBfX38bYweCAFtSWz0RRRFO2a5DTU1syGtL4bAPy1ly81jwzUw0ctwGQRcaHM+R+Kkedqy7KLMWPG4O2338bMmTOh1WoRHR2NkJAQvPTSS83+HEeXXQTorO8HFTUmiExxUAeZOXMm0tPTsWzZMgwdOhS7du3Cli1bEBcXBwDIy8uzu6YKCAhAVlYWSkpKMHz4cPz1r3/FtGnT8K9//cvpY6vLeNdFIGxCSN6GzdVIKQ53NXe0S3N6errd108//TQ+/vhjfPrpp3bdmwVBQHR0tKPDcViQrdTcIgIVRpP8hkKkhNTUVMyePRvDhw9HUlIS1q9f3yizce7cObzxxhsAgNWrV6NHjx7o378/AOs61+eeew733HOP08cmlZnrfVTQqC//GZxaJSAm2Be5xVU4e6la3tebqLUcWXZx6NAh3HvvvXj88ccxefJk5OXl4aGHHsKCBQuwYcOGJu+zfPlyLF26tNXjkTJ5JouImloLfFvxARSRM6SkpCAlJaXJ723cuLHRsf79+zf6EFcJTa7xts2T6lozas0W+LTi/YLIHcnPesbdpBCHAm+pXPCRRx6xO+5Il2aLxYLy8nKEhobaHa+oqEBcXBzMZjOGDh2KJ5980i4wdxa9jxq+PmpU15pRUlnLwJsUNXPmTBQVFWHZsmXIy8vDwIEDW8xsWCwWpKWl4eTJk9BoNOjduzeeeeYZ3HXXXU4fW4UDW4lJuneRAu8qjIwPvfwdiNC2ZRfLly/H2LFj8dBDDwEABg8eDH9/f4wbNw5PPfUUYmJiGt0nLS0Nqamp8tdlZWUt9jzw06ohCNb1fOWGWgbe5PXkrubqxqXmgPUDW2nJHpGnYXM1UppDgXd7ujRLnn/+eVRWVmLGjBnysf79+2Pjxo0YNGgQysrK8OKLL2Ls2LE4ePAg+vTp0+TjtHWPSgAI8fNBdakZJdVG9IBfq+9H1BaOZDbuueceRbLbTXGksZqkbksxNlij1qu/7OLGG2+Uj2dlZeGGG25o8j5VVVXQaOyfm2q1NTBurixcp9NBp9O1elyCICBAp0F5jQkVNSZEBrb6rkQeqS7jXZfV9lGroPdRoabWgvIaBt7kuaQKLC49IqW0qV6oLV2aASAjIwNLlixBZmYmIiMj5eOjR4/GrbfeiiFDhmDcuHF477330Ldv3xbX8rWnk3Owr7TOmw3WyHtVGq1rvFvTWE3SLcT6QdU5Bt7kIEcbSk2bNg2bN2/G2rVrceLECXz99de49957MXLkSHTt2tVp4+L6VaI6Ta3xBur6IZSxszl5MPZWI6U5lPFuT5fmzMxMzJs3D++//z6uu+66Fs9VqVQYMWIEjh492uw5jpYU1hfCvbyJ5Ix3a/bwlkgZb+7lTY5ydNnFbbfdhvLycrz88st44IEHEBISgmuvvRYrVqxw6rikpRbsbE4EeacLdYPAO0ivQWGFgZ3NySsw4U1KcSjwbku5IGDNdN9xxx3IyMho1VYwoijiwIEDGDRoULPnOFpSWF+Ir7VMqrTK2Kb7E3mCijaUmneTS82rFBkTeTZHG0p1xNKLAHkvbwYURPIa7waBdyDnCXkDrvEmhTnc1dzRLs0ZGRmYM2cOXnzxRYwePVrOlvv6+iI4OBgAsHTpUowePRp9+vRBWVkZ/vWvf+HAgQNYvXq1s35PO138WWpOJK/xdqjUXNrLuwYWiwiVioVZ5N7qAgq+HxA11dUcqPcBlYHzhDyXAGmNt4sHQh7L4cDb0XLBV155BSaTCQsXLsTChQvl43PnzpUzHCUlJbjzzjuRn5+P4OBgDBs2DLt27cLIkSPb+es1LdiW8WapOXmz6lppjXfrS82jg/VQCYDRbEFhhQGRQXqlhkfUIVhqTlSnuTXegfX2vCfyVHVdzRl5kzIcDrwBx8oFd+zYcdnHe+GFF/DCCy+0ZShtIq3xvsRSc/JiNbbmanqf1gfePmoVYoJ9ca6kGmdLqhl4k9uTMt4MKIia7moO1M2TMs4T8mCs4SOltamrubsLsXU1L2WpOXmxGpN1LZ+jexdL5ebcUow8QaDU1ZwZb6K6jLe6mVJzBt7kRLt27cK0adPQtWtXCIKAjz76qMXzd+zYAUEQGt1+/fVXp46LpeakFO8MvNnVnAjVbch4A3UN1rilGHkCqdScAQURYGqmq7n8ARV7IZATVVZWYsiQIXj55Zcdut+RI0eQl5cn3/r06eOU8bRiZ2SidmlTqbm7k9d4s9ScvJi0xtvXwcC7bksxdjYn91cXeDOgIJK6mvs0sZ0YwA+oyLmmTJmCKVOmOHy/yMhIhISEOH08dc3VmPImZXhlxlvqal7KjDd5sbrA27GXAZaakyeR13iz1Jyo2TXeQcx4UycybNgwxMTEYOLEifjqq68ue77BYEBZWZndrSkCtxMjhXll4B0iZ7xr+akWeS1DLUvNidhcjahOc2u8uY83dQYxMTFYv349Nm3ahM2bN6Nfv36YOHEidu3a1eL9li9fjuDgYPkWGxvb5HmsNCeleWWpubTG22QRUWk0y6WGRN5Ezng72Fytexc/AMC5kmqIogiBi6LIjQXopEweAwqi5vbxrlvjzXlCrtOvXz/069dP/jopKQlnzpzBc889h/Hjxzd7v7S0NKSmpspfl5WVNRt8A2yuRsrxyoy33kcNncb6q1+q5Dpv8k5tba4WE2zdQqzKaMYl7gxAbo6l5kR1mt3HW95OjK/51LmMHj0aR48ebfEcnU6HoKAgu1uTbIkE7uNNSvHKwBuoy3pznTd5q+pa23ZiDgbeeh81IgJ1AFhuTu4vgAEFkaz5jDdLzalzysnJQUxMjFMeS3rWM+NNSvHaGusQXy0ulBlQwowdeamaNpaaA9bO5hfLDThXUoVB3YOdPTSiDhPsW1dCa7aIjQIOIm8idTVvnPG2zpMKA+cJOU9FRQWOHTsmf33y5EkcOHAAoaGh6NGjB9LS0nDu3Dm88cYbAID09HT07NkTV155JYxGI9566y1s2rQJmzZtcsp42FyNlOa1gXewvJc3S83JO0mBt17jeODdLcQXObkl7GxObk8KvAFrx+YQP60LR0PkWnX7eNsXREoZb8AafNefN0RttX//flxzzTXy19I67Llz52Ljxo3Iy8tDbm6u/H2j0YgHH3wQ586dg6+vL6688kp89tlnmDp1qlPGI7C9GinMawPvLlLgzYw3eam65mqOrziROpsz8CZ356NWwU+rRpXRjNJqBt7k3Zpb4633UUOrUcFosqC8ppaBNznF1Vdf3eLuQhs3brT7+uGHH8bDDz+s8KhYak7K8d413rYtxbjGm7xVW5urAUB3217e50oYeJP7k4IIvh+Qt2tujTcABHGdN3m4uk1aGHmTMrw38JYz3iw1J+9jsYgwmNrWXA2o21KMGW/yBEF6Bt5EQPP7eAPcUow8H5urkdK8NvCW1nhzOyTyRlLQDbQt4y2Vmp+7VOW0MRG5ipTxLqtmQEHezWRrrtZUxruuszmvm8gzyc3VGHiTQrw28JZKzbnGm7yRtL4baGPgbSs1L6sx8SKM3F4QS82JADS/xhvglmLk+QSBzdVIWd4beMv7eLPUnLyPFHhrNao2bQvjr9PIDQq5zpvcHdd4E1nVrfFufHkYqJNKzTlPyLOJXONNCvHewNuXXc3Je0mN1dqyvlsidzYvZuBN7o2BN5FVazLeZcx4k4djqTkpxXsDb9uWMSW80CIvJO3h3a7Am53NyUMw8CayqtvHu/nmamXMeJOHktd4u3YY5MG8OPC2XWhV1ba4hyCRJ5ICb71P218CpM7mDLzJ3QX72jJ5DLzJy3GNN3kzwdbXnGEBKcXrA2+j2WLXaIrIG1TXtn0Pb4mU8T7Lzubk5oL9mPEmAlruai41IWTgTUTUNl4bePv6qKFVW399bilG3kZe461t/xrvc9zLm9wcS82JrFrex5vbiZFnqys1Z8qblOG1gbcgCHKWo6SKnc3Ju1RzjTeRLEjPwJsIaLmreRBLzcnDyR83Me4mhbQp8F6zZg3i4+Oh1+uRmJiI3bt3N3vu5s2bMWnSJERERCAoKAhJSUnYunVro/M2bdqEAQMGQKfTYcCAAfjwww/bMjSHSJ3NS5nxJi9jqLWWE7an1DzWtsa7sMIoZ9CJ3BEz3kRWLa/x5nZi5NnYXI2U5nDgnZmZiUWLFmHx4sXIycnBuHHjMGXKFOTm5jZ5/q5duzBp0iRs2bIF2dnZuOaaazBt2jTk5OTI5+zbtw8zZ87E7NmzcfDgQcyePRszZszAt99+2/bfrBW6sLM5eSmDyRoo6zRtL3oJ8tUgQGfNgDDrTe5MCrzLamphsfCSi7xXXcabzdXI+who/LwnciaHr7pXrVqFefPmYf78+UhISEB6ejpiY2Oxdu3aJs9PT0/Hww8/jBEjRqBPnz54+umn0adPH3z66ad250yaNAlpaWno378/0tLSMHHiRKSnp7f5F2uNulJzBt7kXQwma8Zb247AWxAElpuTR5CaRokiUG5gUEHey2S2vje0nPHmHCHPxt2OSCkOXXUbjUZkZ2cjOTnZ7nhycjL27t3bqsewWCwoLy9HaGiofGzfvn2NHnPy5MktPqbBYEBZWZndzVFSqXlJNdd4k3eRAu/2ZLwBoHsXdjYn96f3UctzgVuKkTdrTca7wmCSS9KJPAlLzUlpDl11FxYWwmw2Iyoqyu54VFQU8vPzW/UYzz//PCorKzFjxgz5WH5+vsOPuXz5cgQHB8u32NhYB34Tq/p7eRN5E6MTMt4AO5uT5+A6b6L6a7wbvzdIgTdgDb6JPBUT3qSUNl11C4L9J6GiKDY61pSMjAwsWbIEmZmZiIyMbNdjpqWlobS0VL6dOXPGgd/AKsS2xvsSu5qTl6nLeLe9uRrAzubkGEcacwLWyqbFixcjLi4OOp0OvXv3xmuvvabI2OR13gy8yYvJGe8mthPTadTyh7VssEaeSIo7GHeTUjSXP6VOeHg41Gp1o0x0QUFBo4x1Q5mZmZg3bx7ef/99XHfddXbfi46OdvgxdToddDqdI8NvRLrQ4hpv8jZGp5WaWzubM+NNlyM15lyzZg3Gjh2LV155BVOmTMGhQ4fQo0ePJu8zY8YMXLhwARs2bMAVV1yBgoICmEzKZNqY8SZquas5YN16r7DCgLJqE9ClI0dGpDy2ViOlOXTVrdVqkZiYiKysLLvjWVlZGDNmTLP3y8jIwG233YZ33nkH119/faPvJyUlNXrMbdu2tfiYziCVmrOrOXkbqau5s0rNzzLwpstwtDHn559/jp07d2LLli247rrr0LNnT4wcOVKx9wUG3uTtRFGUA++m1ngD9ffy5jwhz8XmaqQUh6+6U1NT8eqrr+K1117D4cOHcf/99yM3NxcLFiwAYC0BnzNnjnx+RkYG5syZg+effx6jR49Gfn4+8vPzUVpaKp9z3333Ydu2bVixYgV+/fVXrFixAl988QUWLVrU/t+wBdJ2YlzjTd7G2aXmF8pr5Cw6UUNtacz5ySefYPjw4Vi5ciW6deuGvn374sEHH0R1tTIf8tQtPeL7AXmn+g3TfJpY4w1wSzHybGyuRkpzqNQcAGbOnImioiIsW7YMeXl5GDhwILZs2YK4uDgAQF5ent2e3q+88gpMJhMWLlyIhQsXysfnzp2LjRs3AgDGjBmDd999F4899hj+8Y9/oHfv3sjMzMSoUaPa+eu1LJhdzclLOau5WniAFjqNCgaTBfmlNegR5ueM4ZGHaUtjzhMnTmDPnj3Q6/X48MMPUVhYiJSUFBQXFze7zttgMMBgMMhfO7LbRai/9f2APT/IW5nqBd5NrfEG6m0pZuAHVOR55Gc9I29SiMOBNwCkpKQgJSWlye9JwbRkx44drXrM6dOnY/r06W0ZTpuFcB9v8lJSqXl713gLgoBuXXxx4mIlzl6qYuBNLXKkiabFYoEgCHj77bcRHBwMwFquPn36dKxevRq+vr6N7rN8+XIsXbq0TWPr4m/NeBdXMvAm71Q/493cGm9mvMmT1TVXY+RNymjfVbebk0oLDSYLamrNLh4NUcdxVsYbqCs3P8vO5tSMtjTmjImJQbdu3eSgGwASEhIgiiLOnj3b5H3as9tFqFRqzsCbvJRdxpuBN3mhVmzQRNQuXh14+2vV8qe6zHqTNzE4qas5AHTnXt50GW1pzDl27FicP38eFRUV8rHffvsNKpUK3bt3b/I+Op0OQUFBdrfWkjPeLDUnL1U/461uJgIJ0nPbPfJ87K1GSvHqwFsQBLncnOv6yJs4azsxoG5LMXY2p5Y42phz1qxZCAsLw+23345Dhw5h165deOihh3DHHXc0WWbeXqH+zHiTdzNZrO8LKgFQNZPxZvd/8mTSs56BNymlTWu8PUmInxaFFUZmvMmrOKurOVBXan6upKrdj0Wey9HGnAEBAcjKysI999yD4cOHIywsDDNmzMBTTz2lyPikXS64xpu8Vd0e3s1/IMveOOTRuMabFMbAW/70lhdb5D2cusZbKjXnGm+6DEcacwJA//79G5WnK0XKeJfVmGAyW6BRe3VBGHkhk7nlPbwBIFjahpUZb/JA0hILM3dHJYV4/ZUFP70lb+SsruZA3RrvvJIauzWCRO4k2NdHbqxTwqCCFLZmzRrEx8dDr9cjMTERu3fvbtX9vv76a2g0GgwdOtTpY6rLeLcQeMvbsHKOkOeRPm+1sNacFOL1gXewr/XTW76JkFIcucDavHkzJk2ahIiICAQFBSEpKQlbt251+picmfGODNRDoxJgsoi4UFbT7scjcgW1SpAroLjOm5SUmZmJRYsWYfHixcjJycG4ceMwZcoUu6UWTSktLcWcOXMwceJERcYldTVvbg9voF6VIPvikAdSyRlvBt6kDK8PvJnxJiU5eoG1a9cuTJo0CVu2bEF2djauueYaTJs2DTk5OU4dlzPXeKtVAmJC9ADYYI3cG/fypo6watUqzJs3D/Pnz0dCQgLS09MRGxuLtWvXtni/u+66C7NmzUJSUpIi43Ik481Sc/JE0jILZrxJKQy8ucabFOToBVZ6ejoefvhhjBgxAn369MHTTz+NPn364NNPP3XquJyZ8QaAWFtn8zPFbLBG7kvey5vZPFKI0WhEdnY2kpOT7Y4nJydj7969zd7v9ddfx/Hjx/HEE0+06ucYDAaUlZXZ3S5H6mre0hpvKVlRaTTL7yNEnkJ67jPjTUph4C1tJ1bJT2/Judp6gVWfxWJBeXk5QkNDmz2nLRdYztzHGwDiwqyB92kG3uTG6jLefD8gZRQWFsJsNiMqKsrueFRUFPLz85u8z9GjR/HII4/g7bffhkbTup64y5cvR3BwsHyLjY297H1a09U8UF/XC4FZb/I0Uqk5M96kFAbeftIab2Y4yLnacoHV0PPPP4/KykrMmDGj2XMcvcASRRFGs3Mz3j1C/QEAp4sqnfJ4RK7AjDd1FEGwzyqLotjoGACYzWbMmjULS5cuRd++fVv9+GlpaSgtLZVvZ86cuex95DXeLWS81SoBgTpr8M/AmzyNXGrOYg5SCLcT4xpvUlhrL7AaysjIwJIlS/Dxxx8jMjKy2fPS0tKQmpoqf11WVtZi8G2oVx7orIx3TynjXcSMN7mvEH/r+wHXeJNSwsPDoVarG334WlBQ0OhDWgAoLy/H/v37kZOTg7vvvhuAtRJKFEVoNBps27YN1157baP76XQ66HQ6h8bWmjXegDVhUVZj4hI98jhyczVmvEkhDLx9meEgZTh6gVVfZmYm5s2bh/fffx/XXXddi+c6eoFlrLdBpdMy3rbAO5el5uTG5Iw3A29SiFarRWJiIrKysnDjjTfKx7OysnDDDTc0Oj8oKAg//fST3bE1a9Zg+/bt+OCDDxAfH++0sbVmH2/AmrDILWbGmzwP13iT0rw+8A4LqOti29pMJFFrOHqBJcnIyMAdd9yBjIwMXH/99U4fl6G2XuCtdtYab2upeXGlEeU1tQjU+zjlcYk6krzGmx/EkoJSU1Mxe/ZsDB8+HElJSVi/fj1yc3OxYMECANYqpnPnzuGNN96ASqXCwIED7e4fGRkJvV7f6Hh7mVtRag7U28ublYLkYbiPNymNgbct8K41iyirNiHYjwEDOY8jF1iANeieM2cOXnzxRYwePVrOlvv6+iI4ONgpY5Iy3jqNymkfNAXoNAjz16Ko0ojTRVUY2M05YyXqSMx4U0eYOXMmioqKsGzZMuTl5WHgwIHYsmUL4uLiAAB5eXmX3dNbCVJXc00L+3gDDLzJc3Efb1Ka1zdX02nUCNJbP3+4WGFw8WjI08ycORPp6elYtmwZhg4dil27drV4gfXKK6/AZDJh4cKFiImJkW/33Xef08ZkqDUDcF6ZuaQH13mTm5My3kUMvElhKSkpOHXqFAwGA7KzszF+/Hj5exs3bsSOHTuave+SJUtw4MABp4+pLuPd8nsD9/ImT8VSc1Ka12e8ASA8QIeyGhMKKwy4IjLA1cMhD5OSkoKUlJQmv7dx40a7r1u62HKWuoy32qmPGxfqh5zcEpwuZmdzck/htgqoogoG3uR9TK1ursbAmzyTmtuJkcK8PuMNWANvgBdb5B2kNd7O6mgukdZ55zLjTW5Kei+orjWj0mBy8WiIOlZr13hLTWlL2AuB2mnXrl2YNm0aunbtCkEQ8NFHH132Pjt37kRiYiL0ej169eqFdevWOW08AkvNSWEMvFG3zruQpebkBaTtxJwfeLPUnNybv04DP621EuRiOd8PyLu0NuPNUnNylsrKSgwZMgQvv/xyq84/efIkpk6dinHjxiEnJwePPvoo7r33XmzatMkp45H38WbcTQphqTnqshwMvMkbGG2Bt7PXeMdxSzHyAOEBOuQWV6GwwoCe4f6uHg5RhzHZliFdtqu5rdS8hIE3tdOUKVMwZcqUVp+/bt069OjRA+np6QCAhIQE7N+/H8899xxuvvnmdo9H6mrOjDcphRlv1M94s2yKPJ/BZG2u5uyMd49Qa5ByvrRa/hlE7iYi0PpBLDPe5G0cznizqzl1sH379iE5Odnu2OTJk7F//37U1rb/+ciu5qS0Nl15r1mzBvHx8dDr9UhMTMTu3bubPTcvLw+zZs1Cv379oFKpsGjRokbnbNy4EYIgNLrV1NS0ZXgOY8abvIlSGe/wAC38tGqIInCmuNqpj03UUcK59Ii8VGu7mrO5GrlKfn4+oqKi7I5FRUXBZDKhsLCw2fsZDAaUlZXZ3ZpSV2rOwJuU4fCVd2ZmJhYtWoTFixcjJycH48aNw5QpU5rdc9JgMCAiIgKLFy/GkCFDmn3coKAg5OXl2d30er2jw2uTuuZqvNAizyd1NXd24C0IAnqESuXm7GxO7okZb/JWre5qLjVXq66FyACFOpjUAE0iPQcbHq9v+fLlCA4Olm+xsbFNnseMNynN4SvvVatWYd68eZg/fz4SEhKQnp6O2NhYrF27tsnze/bsiRdffBFz5sxBcHBws48rCAKio6Ptbh0lnKXm5EVqzdLFlfNXmvS0dTZngzVyV9IHsRf5fkBexiyt8Va3bjsxs0VEWQ27/1PHiY6ORn5+vt2xgoICaDQahIWFNXu/tLQ0lJaWyrczZ840eR4z3qQ0h668jUYjsrOzG62vSE5Oxt69e9s1kIqKCsTFxaF79+74wx/+gJycnHY9niNYak7epNZ2ceWjdn7gzc7m5O6Y8SZv1dqMt95HDX9b9//iSn5ARR0nKSkJWVlZdse2bduG4cOHw8fHp9n76XQ6BAUF2d2aIgXezHiTUhy68i4sLITZbG5yfUXDT6Ac0b9/f2zcuBGffPIJMjIyoNfrMXbsWBw9erTZ+7R2vUZrSM3VqoxmVBn56S15tlq51Lzli6u2iLd1gT5RyFJzck91GW8G3uRdzJbWV0N18bdeNxVXcp5Q21VUVODAgQM4cOAAAOt2YQcOHJCXr6alpWHOnDny+QsWLMDp06eRmpqKw4cP47XXXsOGDRvw4IMPOmU8LDUnpbUp5dXU+oqW1lZczujRo3HrrbdiyJAhGDduHN577z307dsXL730UrP3ae16jdYI0GnkDs9FLC8kDyc1V1Mi490rIgAAcLygwumPTdQRpIx3ITPe5GVam/EGgDA58GaDNWq7/fv3Y9iwYRg2bBgAIDU1FcOGDcPjjz8OwNqguX4Pqfj4eGzZsgU7duzA0KFD8eSTT+Jf//qXU7YSA7iPNynPoX28w8PDoVarm1xf0TAL3h4qlQojRoxoMeOdlpaG1NRU+euysrI2B9+CICA8QIdzJdW4WGFArK1BFJEnktZ4KxN4120pVlNrht5H7fSfQaSkiHoZ7/Z+qEzkTuSu5pdZ4w0Aocx4kxNcffXVLTbo27hxY6NjEyZMwA8//KDIeNTMeJPCHLry1mq1SExMbLS+IisrC2PGjHHaoERRxIEDBxATE9PsOa1dr9FaUoM1ZrzJ05kUXOMd5q9FkF4DUQROstyc3JBUam40WVBu4NIj8h6OZLxD/W27wXCNN3kQaZWFmc3VSCEOZbwBaxnI7NmzMXz4cCQlJWH9+vXIzc3FggULAFgz0efOncMbb7wh30dau1FRUYGLFy/iwIED0Gq1GDBgAABg6dKlGD16NPr06YOysjL861//woEDB7B69Won/IqtI5UXFpR3zN7hRK4ir/FuRVbDUYIgoFdEAA6cKcGJi5VIiGnfB2JEHc1Xq0aAToMKgwkXyw0I0jffsIfIk5gttq7mrSk1tyUripmsIA8il5oz400KcTjwnjlzJoqKirBs2TLk5eVh4MCB2LJlC+Li4gA0Xo8BQF67AQDZ2dl45513EBcXh1OnTgEASkpKcOeddyI/Px/BwcEYNmwYdu3ahZEjR7bjV3NMZJB1z/ALZSybIs9mVLDUHLCWm1sDb67zJvcUEahDhcGEwnIDetv6FhB5Oscy3lKpOQNv8hxyqTkz3qQQhwNvAEhJSUFKSkqT32tqPUZL6zcA4IUXXsALL7zQlqE4TbQUeJcy402eTd5OTKNM4C0FKuxsTu4qPECLk4WV7GxOXsVs+1BW3Yqu5lLgzVJz8iQq24dOotj+xtFETVHmytsNSYF3fhkDb/JsSu7jDQC9bQ3WjjPjTW6KFVDkjdrW1ZyBN3kOdb1Amw3WSAkMvG2igqULLQbe5NmUXOMN1G0pduJi5WWrXYg6I7kCiu8H5EXkruYsNScvpar33Ge5OSmBgbcNL7TIWxhNyq7xjgvzg0qA3JyKyN3E2D6IzePSI/IiXONN3q7+h06Mu0kJDLxtooKsXc0vVdWiptbs4tEQKUfKeGsUCrx1GjW6d/EDABy/yHXeVGfNmjWIj4+HXq9HYmIidu/e3ar7ff3119BoNBg6dKiyA7SJtgXe+aXVHfLziDoDuau5A/t4V9eaUW3kNRN5Bpaak9IYeNsE+/pAZ2s2VcB1feTBlC41B6ydzQHgRCHXeZNVZmYmFi1ahMWLFyMnJwfjxo3DlClTGu2C0VBpaSnmzJmDiRMndtBImfEm7+RIxjtAp4HW9uFtUSWvmcgz1O8ryFJzUgIDbxtBEOqyHCw3Jw+mdHM1AOgVXrfOmwgAVq1ahXnz5mH+/PlISEhAeno6YmNjsXbt2hbvd9ddd2HWrFlISkrqoJECUfWWHnE/V/IWdWu8L//eIAgCy83J49TPePO1n5TAwLueKK7zJi+g9D7eANA70prxPlbAjDcBRqMR2dnZSE5OtjuenJyMvXv3Nnu/119/HcePH8cTTzyh9BDtRAbqIQhArVlEcRWDCvIOjmS8AW4pRp6n/hpvEwNvUkCb9vH2VGywRt7ApPA+3gDQJzIQAANvsiosLITZbEZUVJTd8aioKOTn5zd5n6NHj+KRRx7B7t27odG07q3KYDDAYKgrey0rK2vTeLUaFcIDdLhYbkB+aQ3CA3Rtehwid1K3j3frAu+wAFvgXcHAmzyDIAjwUQuoNYtydSCRMzHjXU9dQx0G3uS5OmKNd98oa6n5uZJqlNfUKvZzyL0Igv1zThTFRscAwGw2Y9asWVi6dCn69u3b6sdfvnw5goOD5VtsbGybx8p13uRtHM14R9g+kOLuFeRJpN4FRhMDb3I+Bt71RAZa30S4xps8WUeUmof4aeWdAo4y6+31wsPDoVarG2W3CwoKGmXBAaC8vBz79+/H3XffDY1GA41Gg2XLluHgwYPQaDTYvn17kz8nLS0NpaWl8u3MmTNtHrNUAcXO5uQt5K7mrQ28ba/xBeW8ZiLPodUw8CblsNS8nq4hvgCY4SDPVmtSvrkaAPSNCsSFMgN+yy/HVT26KPqzqHPTarVITExEVlYWbrzxRvl4VlYWbrjhhkbnBwUF4aeffrI7tmbNGmzfvh0ffPAB4uPjm/w5Op0OOp1zysLZbJO8jZzxbmU1VGSgdY4UMONNHkQKvA0MvEkBDLzr6d7FGnifvVTl4pEQKacjupoD1sB799FCHLlQrujPIfeQmpqK2bNnY/jw4UhKSsL69euRm5uLBQsWALBmq8+dO4c33ngDKpUKAwcOtLt/ZGQk9Hp9o+NKiWapOXkZR7qaA3VVghe5BSt5EOnaiGu8SQkMvOvpZst4F5QbYDCZodOoXTwiIueT13hrlFvjDQD9oqwN1n5j4E0AZs6ciaKiIixbtgx5eXkYOHAgtmzZgri4OABAXl7eZff07kgx7PlBXsbRNd5S4M1Sc/IkLDUnJTHwrifUXwtfHzWqa83IK6lBz3B/Vw+JyOlqO2CNNwD0jbYG3kfyucabrFJSUpCSktLk9zZu3NjifZcsWYIlS5Y4f1DNiA7i0iPyLnUZ71YG3kEsNSfPIzdXY8abFMDmavUIglCv3JwNdcgzGTuo1LxPpLWzeWGFAUUVvDAj9yJVQJ0rqYaF+7mSF2hrxrvKaEaFwaTYuIg6ko4Zb1IQA+8GuM6bPF3dGm9lS839dRrEhlrn028XmPUm9xIToodKsF58XeQHR+QFHO1q7q/TwF9rXZJXwCaE5CG4xpuUxMC7ge5d/AAw402eq6O6mgNc503uy0etQkyw9YOjM8X8IJY8n8ksZbxb/97AcnPyNOxqTkpi4N1Aty515YVEnqij1ngD1s7mAANvck9SxQY/iCVvYHJwjTcARMgN1hh4k2dgczVSEgPvBlhqTp5MFMUOW+MNAP1sDdZ+zWfgTe5HqoBixpu8gdnBfbyBep3NWWpOHsKHzdVIQQy8G2CpOXkyc70mUdoOCLwHxAQBAA7nldn9bCJ3ECsF3vwglryAycE13kBdxvsiM97kIaSMdy0z3qQABt4NSBnv/LIalpmQx5HKzAHAR+F9vAGgV0QA9D4qVBnNOFVUqfjPI3ImlpqTNzGbHetqDgCRgVzjTZ5Fx4w3KYiBdwNh/lr4adUQRWY5yPPUfyPpiFJztUpAgi3r/fO5UsV/HpEzxYYy403eo247sda/N0QFWTPe+dzvnjwE13iTktp05b1mzRrEx8dDr9cjMTERu3fvbvbcvLw8zJo1C/369YNKpcKiRYuaPG/Tpk0YMGAAdDodBgwYgA8//LAtQ2s3QRAQH+4PADh5kRk68iz1t8dwJKvRHgO7BgMADp0v65CfR+QsUgXU+ZIamJj9IA/XljXeXW373eeVsiqEPIPex7pFXnWt2cUjIU/kcOCdmZmJRYsWYfHixcjJycG4ceMwZcoU5ObmNnm+wWBAREQEFi9ejCFDhjR5zr59+zBz5kzMnj0bBw8exOzZszFjxgx8++23jg7PKeTAu5CBN3kWKfDWqlUQhI4JvK/sast4n2fGm9xLVKAeWrUKZouIPGb0yMO1pat5N1vgfb60Bhb28SAP4K+zBt6VBgbe5HwOB96rVq3CvHnzMH/+fCQkJCA9PR2xsbFYu3Ztk+f37NkTL774IubMmYPg4OAmz0lPT8ekSZOQlpaG/v37Iy0tDRMnTkR6erqjw3OKXlLgzTWp5GFqTdJWYh0TdAPAwG7Wef/L+TKIIi/MyH2oVIK8xSTXeZOnkzPeDgTeUUF6CIK1LLeo0qjU0Ig6jJ9WAwCoMppcPBLyRA4F3kajEdnZ2UhOTrY7npycjL1797Z5EPv27Wv0mJMnT27xMQ0GA8rKyuxuzhIfwVJz8kzyVmKajmvv0CcqABqVgJKqWpwrYfBC7kVa551bzPcD8mxt6Wqu1ajkLcXO8/WdPIC/1pbxNjLjTc7n0NV3YWEhzGYzoqKi7I5HRUUhPz+/zYPIz893+DGXL1+O4OBg+RYbG9vmn99QfHgAAJaak+eRSs0daZ7TXjqNGn2jrPt5/8J13uRm4sOsgfcJvh+QhzO3obkaULfOm4E3eQI/nS3jbWDGm5yvTVffDdeGiqLY7vWijj5mWloaSktL5duZM2fa9fPriw+zZrzzy2pQyYlHHqRujXfHlZoDwMBu1nXev7CzObmZXhHWD2JPsAKKPFxb1ngDdeu8WdFEnsDfVmrOjDcpwaHAOzw8HGq1ulEmuqCgoFHG2hHR0dEOP6ZOp0NQUJDdzVmC/XwQ5q8FAO49TB6l1gWl5kDdOu+DZxl4k3vpZVt6dOJihYtHQqQci0WE1ILD0R0v5AZrJWxASO7Pz9ZcjWu8SQkOXX1rtVokJiYiKyvL7nhWVhbGjBnT5kEkJSU1esxt27a16zHbS+psziwHeRKj3FytYwPvobEhAIADZ0rYYI3civRekFtcxS3FyGkc2ZZ18+bNmDRpEiIiIhAUFISkpCRs3brVqeMx1etIrnawIoql5uRJpIx3FbuakwIcvvpOTU3Fq6++itdeew2HDx/G/fffj9zcXCxYsACAtQR8zpw5dvc5cOAADhw4gIqKCly8eBEHDhzAoUOH5O/fd9992LZtG1asWIFff/0VK1aswBdffNHsnt8dQcpyHC1gloPaR4l979tKap7TUXt4S/pHB0GnUaG0upa9E8itdA32hU6jQq1ZZGdzcgpHt2XdtWsXJk2ahC1btiA7OxvXXHMNpk2bhpycHKeNyVwv8Hb0/UEOvLmXN3kAP7m5GjPe5HwOB94zZ85Eeno6li1bhqFDh2LXrl3YsmUL4uLiAFgDh4ZvHsOGDcOwYcOQnZ2Nd955B8OGDcPUqVPl748ZMwbvvvsuXn/9dQwePBgbN25EZmYmRo0a1c5fr+36RVtL14/ksxkUtZ0S+963h8lsvbjSdnCpuVajwiBbuXlObkmH/myi9lCphLoKqEJ+EEvt5+i2rOnp6Xj44YcxYsQI9OnTB08//TT69OmDTz/91Gljkj6UBRxf4901RA+AGW/yDP46ZrxJOW26+k5JScGpU6dgMBiQnZ2N8ePHy9/buHEjduzYYXe+KIqNbqdOnbI7Z/r06fj1119hNBpx+PBh3HTTTW0ZmtMkRFu7MB/JL3fpOMi9KbHvfXtIa7wdvbByhmE9QgAAOWcudfjPJmqPunXerNag9nHGtqwWiwXl5eUIDQ112rjsM96OXRpKa7wLK4yoqWWwQu7Nv17Gm0vjyNk6Nu3lRvrZAu/TxVVssEBtotS+9w05sqe9tI7PpwO3E5MM69EFADPe5H562baY5JZi1F7O2Jb1+eefR2VlJWbMmNHsOY68LwD2a7wd/Vw22NcHgbYs4ZniKsfuTNTJSBlviwhUsbM5ORkD72aEBegQHqCDKAK/XWB5ITlOqX3vG3JkT3vp4krTwduJAXUZ71/zy/lhFrmVumabfC8g52jrtqwZGRlYsmQJMjMzERkZ2ex5jrwvAPX38BYc3h5WEAT0tM0R9vAgd+enVUNra0BbUl3r4tGQp2Hg3YL+crk513lT2ymx7319juxpb3JhqXlMsC+ig/QwW0T8xG3FyI30jrRmvI+x2Sa1U3u2Zc3MzMS8efPw3nvv4brrrmvxXEfeF4C27+EtiQvzAwCcLmLGmxznSBPaHTt2QBCERrdff/3VKWMRBAEhfj4AgEuVRqc8JpGEgXcLpMD7V67zpjZQat/7hhzZ015qrtbR24lJpKz3Dyw3JzfSNyoAgmBdw3qx3ODq4ZAba+u2rBkZGbjtttvwzjvv4Prrr7/sz3HkfQEAzOa6jHdb9AyzZrxPFTHjTY5xtAmt5MiRI8jLy5Nvffr0cdqYQv21AIBLVQy8ybkYeLdAWud9OI8Zb3KcUvvet0d7sxrtlRhnXee9/1SxS34+UVv4aTWIC7Vm9Nhwk9rL0W1ZMzIyMGfOHDz//PMYPXo08vPzkZ+fj9JS51UOSV3NmfGmjuZoE1pJZGQkoqOj5ZtarXbamOSMdxVLzcm5GHi3YKBt+6Ofz5XBYmFnQ3KcEvvet4d0ceXjgjXeADC6VxgA4LuTxXZddIk6u/62LSZ/5dIjaidHt2V95ZVXYDKZsHDhQsTExMi3++67z2ljktd4t7EaKp5rvKkN2tOEdtiwYYiJicHEiRPx1VdfOXVcXfxsGW+WmpOTaVw9gM6sT2QAfH3UqDCYcKKwAldEBrp6SORmZs6ciaKiIixbtgx5eXkYOHBgq/a9l2RnZ+Odd95BXFxcoy342qJWLid0zWduCTFBCNRrUF5jwqHzZRjU3flbphEpoV90ID7/JZ9Lj8gpUlJSkJKS0uT3Nm7caPd1wy1aldD+Nd7WwPt8aTUMJjN0GudlH8lztaUJbUxMDNavX4/ExEQYDAa8+eabmDhxInbs2GG3vXF9BoMBBkPdMqHLdfnvwlJzUggD7xZo1CoM6haM704V48CZUgbe1CaOXGABUHTfSLMt493WdXztpVYJGBUfii8OF+CbE0UMvMltJMRIPT+Y8SbPU7+reVuEB2jhr1Wj0mjGmeJqXGFrSEjUGo40oe3Xrx/69esnf52UlIQzZ87gueeeazbwXr58OZYuXdrq8YTaMt7FzHiTk7HU/DKGxFoDgx/Plrh2IEROIGe8XVRqDtSVm39zoshlYyBylFRqfvRChbw7AJGnaG/Gu/6WYqdYbk6t5KwmtKNHj8bRo0eb/b6jXf6jgnQAgLzSmlaPgag1GHhfxuDuIQCAg2dKXDoOImdo7zo+ZxgVz3Xe5H56hPrB10cNg8mCU2wgRR7GGdVQUuB9opDb7lHrOKsJbU5ODmJiYpr9vqNd/mOCfQEA+Qy8yclYan4ZQ2NDAACH8spQU2uG3ofrlsh9SZk6V5WaA8CArkEI1GlQbuA6b3IfKpWAvtGBOHimBIfzylhKSx5F2mqyPTte9I0MxGfIw28XGHhT66WmpmL27NkYPnw4kpKSsH79+kZNaM+dO4c33ngDAJCeno6ePXviyiuvhNFoxFtvvYVNmzZh06ZNThtTTIgeAJBXWu20xyQCmPG+rO5dfBEeoEWtWcRP55y3dQeRK9RaXNtcDbBe2I2MDwUAfH280GXjoI63Zs0axMfHQ6/XIzExEbt372723M2bN2PSpEmIiIhAUFAQkpKSsHXr1g4cbWODu3HpEXkmkxPeG/pFWz+M+u0CGxBS6zna5d9oNOLBBx/E4MGDMW7cOOzZswefffYZbrrpJqeNScp4F1YYYTCZnfa4RAy8L0MQBLk09pvjXJNK7k0q7XbVdmKScX3CAQA7j1x06Tio42RmZmLRokVYvHgxcnJyMG7cOEyZMqVRV3/Jrl27MGnSJGzZsgXZ2dm45pprMG3aNOTk5HTwyOsMtlVnHDzDD2HJs7R3jTcA9I2yNiD87UI5t2Alh6SkpODUqVMwGAzIzs62a5K2ceNGu87+Dz/8MI4dO4bq6moUFxdj9+7dmDp1qlPH08XPBzqNNURiuTk5EwPvVhjdy5qd++YkA29yb7W2UvP2XFw5w9X9IgEA+08Xo8JgculYqGOsWrUK8+bNw/z585GQkID09HTExsZi7dq1TZ6fnp6Ohx9+GCNGjECfPn3w9NNPo0+fPvj00087eOR1pKVHP50rZYM18ijyGu92fCgbF+YPrUaFmloLzlxiHwRyX4IgoGuINet97hLLzcl5GHi3gtSFOfv0JRhNvNgi92Uyu765GmBtwtMzzA+1ZhFfH2O5uaczGo3Izs5GcnKy3fHk5GTs3bu3VY9hsVhQXl6O0NBQJYbYKr0iAhCg06C61oxjF7mOlTyHM9Z4q1UC+th6H3C/e3J38bZmgcfZpZ+ciIF3K1wRGYAwfy1qai1c20duTSon9HFxxhsAJvSNAADsYLm5xyssLITZbG60PUxUVFSjbWSa8/zzz6OyshIzZsxo9hyDwYCysjK7mzOpVQIGdZPKzUuc+thErtTefbwl/aRycwbe5OZ6R9gC7wJ+yErOw8C7FQRBwChbufkeZufIjUnlsWoXr/EG6srNd/12EaLI9YDeQBDsn3eiKDY61pSMjAwsWbIEmZmZiIyMbPa85cuXIzg4WL7Fxsa2e8wNDbGVmx/gOm/yIM5orgYAfaOtgfevbLBGbk7aueI4q5vIiRh4t5KUnfvq1wIXj4So7eoy3q6f+qN7hUGrUeFcSTW3n/Fw4eHhUKvVjbLbBQUFjbLgDWVmZmLevHl47733cN1117V4blpaGkpLS+XbmTNn2j32hobYGqzl5F5y+mMTuYqc8W7nh7L9bYH34fPOrTYh6mi9I6yB91Fen5ATuf7q201cY8vOHTxbioJydjgk9+SMzrXO4qtVY9wV1u7mn//cunJjck9arRaJiYnIysqyO56VlYUxY8Y0e7+MjAzcdttteOedd3D99ddf9ufodDoEBQXZ3ZwtsWcXAMCRC+UoqTI6/fGJXMFZ7w2Du4cAAE4UVqK0ura9wyJymX7RgRAEIL+sBhfLDa4eDnkIBt6tFBmkl7eS2fEr16SSe5JKzV29nZjk9wOjAQD/+znPxSMhpaWmpuLVV1/Fa6+9hsOHD+P+++9Hbm4uFixYAMCarZ4zZ458fkZGBubMmYPnn38eo0ePRn5+PvLz81Fa6toS78hAPXpH+EMUge9OFrt0LETOInc1b2fgHeqvRWyotRv0T2e5HIPcV6DeB1fYst4H2NODnISBtwOu7W/Nen/56wUXj4SobWo7SVdzyaQBUdCoBPyaX45T7Bzq0WbOnIn09HQsW7YMQ4cOxa5du7BlyxbExcUBAPLy8uz29H7llVdgMpmwcOFCxMTEyLf77rvPVb+CbJRtp4tvGXiTh3BmNdQQW9b7IJvRkpsbKvf04NIico42XX2vWbMG8fHx0Ov1SExMxO7du1s8f+fOnUhMTIRer0evXr2wbt06u+9v3LgRgiA0utXUdK6S7on9rWsRd/1WiCoj9x4m9yNlNTpDqTkAhPhpkdTbGsT8j+XmHi8lJQWnTp2CwWBAdnY2xo8fL39v48aN2LFjh/z1jh07IIpio9vGjRs7fuANSFtMfnOiyMUjIXIOsxKBN7OE5OaG9ggBAOw/xcCbnMPhwDszMxOLFi3C4sWLkZOTg3HjxmHKlCl2mYr6Tp48ialTp2LcuHHIycnBo48+invvvRebNm2yOy8oKAh5eXl2N71e37bfSiEDuwUhNtQX1bVmfHGYTdbI/cjN1TpJqTnAcnNyP6PjrbtcHMor4zpW8ghyNZQTGm9Knf+Z8SZ3N7a3tQ9N9ulLKK/haz21n8OvsKtWrcK8efMwf/58JCQkID09HbGxsVi7dm2T569btw49evRAeno6EhISMH/+fNxxxx147rnn7M4TBAHR0dF2t85GEATcMKQbAODjnHMuHg2R40xOvLhyluQB0VCrBPx4tpTbdpBbiAzSo1e4dZ33t8x6kweQ+n+0t6s5YE1SqATgQpkBeaXV7X48IlfpGe6P+HB/mCwivuZ2wuQEDl19G41GZGdnIzk52e54cnIy9u7d2+R99u3b1+j8yZMnY//+/aitrfv0qKKiAnFxcejevTv+8Ic/ICcnx5GhdZg/DesKANj520VcqmRHW3IvJic10HGmiEAdxvexfqq8+YezLh4NUeuMsz1nvzrCZpvk/py51aSfVoOB3azNaLkcg9zd1f2s2wl/yUpXcgKHXmELCwthNpsb7bsaFRXVaH9WSX5+fpPnm0wmFBZaPz3q378/Nm7ciE8++QQZGRnQ6/UYO3Ysjh492uxYDAYDysrK7G4d4YrIQAyICYLJIuLTH893yM8kcpbO1lxNcnNidwDAhz+cg8V2AUjUmV1ta7a540gBRJHPWXJvUjWU2knLkJJsfRD2HWfgTe4teYC1Avfzn/NRbTS7eDTk7tp09S0I9i/Moig2Ona58+sfHz16NG699VYMGTIE48aNw3vvvYe+ffvipZdeavYxly9fjuDgYPkWGxvbll+lTf483BokvLnvNC+4yK1IDXScUU7oTNclRCFIr8H50hrsY4aE3EBSrzDofVTIK63B4bxyVw+HqF2kaigfJ1VDjbY1zeTrObm7UfGh6N7FF+UGE7YdYhNYah+HAu/w8HCo1epG2e2CgoJGWW1JdHR0k+drNBqEhYU1PSiVCiNGjGgx452WlobS0lL5dubMGUd+lXa5ObE7/LRqHC2owDcnuJ0MuY9ac+crNQcAvY8a04ZYl3G8v7/j5jJRW+l91HLjna+OsASR3Juzq6FG9AyFWiXgTHE1zl6qcspjErmCSiXg5qusCbfXvz7FhBu1i0OvsFqtFomJicjKyrI7npWVhTFjxjR5n6SkpEbnb9u2DcOHD4ePj0+T9xFFEQcOHEBMTEyzY9HpdAgKCrK7dZQgvQ/+NMzaZO3Nb0512M8lai9pHV9naq4mmTnCWrWy5ad8FFYYXDwaosuTys2/OHzBxSMhah9nNlcDgACdBoNs67xZbk7u7tbRcdBpVDhwpoRVHNQuDl99p6am4tVXX8Vrr72Gw4cP4/7770dubi4WLFgAwJqJnjNnjnz+ggULcPr0aaSmpuLw4cN47bXXsGHDBjz44IPyOUuXLsXWrVtx4sQJHDhwAPPmzcOBAwfkx+yM5iTFAbCu+ThZWOni0RC1jrkTbicmGdw9BMN6hMBotiDj26a3JyTqTJIHREEQgJzcEmb1yK3VfSjrvPeGugaErAgh9xYRqJOTA89v+429aKjNHA68Z86cifT0dCxbtgxDhw7Frl27sGXLFsTFWQPRvLw8uz294+PjsWXLFuzYsQNDhw7Fk08+iX/961+4+eab5XNKSkpw5513IiEhAcnJyTh37hx27dqFkSNHOuFXVEb/6CBM7B8Jiwi8vP2Yq4dD1CpSqbm6k5WaS24b0xMA8Na3p+WxEnVWUUF6jLLt6f3Zj9yHntxX3Y4XzquGui7BugRx55GLMJjYlIrc29+v7g0/rRrZpy9hM7cUpjZq0ytsSkoKTp06BYPBgOzsbIwfP17+3saNG7Fjxw678ydMmIAffvgBBoMBJ0+ebJTJfuGFF3D69GkYDAYUFBRg69atSEpKasvQOtQ9E/sAAD46cA6ni5j1ps5P6lzr08m6mkumDIxBRKAOF8oM2PITAxnq/KTeBNzlgtxZ3XuD8z6UHdQtGJGBOlQazeyHQ24vJtgX99qu+5/536+oMJhcPCJyR53z6ttNDI0NwYS+ETBbRKz8/Iirh0N0WVI5YWfNeGs1Ksweba2eWfPVcZZzUac3ZWAMNCoBP58rw/GLFa4eDlGbKLHVpEol4LoB1qz3F4fYB4Hc3x1j4xEX5ofCCgO2sMqJ2oCBdzs9MqU/VALw2U95+IYNF6iTk7eM6YRrvCVzx/REoE6DIxfKuXUHdXqh/lpM6BsBAMj8nh35yT2ZLcrseDHJFnj/7+d8uYEbkbvSalS4Yai1ufLuY4UuHg25Iwbe7ZQQE4RZo3oAAJZ88gvXpVKnJpUTdsau5pJgXx/cPrYnAODFL49x6w7q9KT3gPf3n0FNLdeykvupVaC5GgD87opwhPprUVhhYKBCHkFqGrjn6EVW5ZHDOu/VtxtJndQPIX4++DW/HKu/YqM16rykjHdnLTWX3PG7ePhr1TicV4b/spyLOrmr+0Wia7Ael6pq8b+f+Xwl91O3nZhzLwt91Cr80dYHYVP2Wac+NpErDI0NQYBOg0tVtfjpXKmrh0NuhoG3E4T6a7HshoEArB3OfzrLiUidU2dvriYJ8dPizvG9AVibmDCLSJ2ZWiXglpHWrPfGr0+xSoPcjhLN1STTE7sDALYduoDSqlqnPz5RR/JRq+TlRZ+xCSw5qHNffbuRPw7piusHx8BkEZHyTjZKqoyuHhKRHVEU6/Zq7cRrvCV3ju+FmGA9zpVUY8Oek64eDlGLbhnVA3ofFQ6eLcXXx9jvg9xLXeNN518WXtk1CP2jA2E0WfDu97mXvwNRJ/fHodYqjk8OnGe5OTmEgbcT/fNPAxEb6oszxdW4JyMHZk5G6kTqPx+dvY5PCb5aNR7+fT8AwOqvjuFMcZWLR0TUvPAAHf4ywpr1fvmroy4eDZFjlGy8KQgC7vhdPADgP3tPsRcOub2r+0UgUK9BflkNGyuTQxh4O1GInxbrZw+Hr48au48W4rGPfmLJIXUapvqBdycvNZfcMKQbRvYMRZXRjLTNnE/Uud05vhd81AK+OVGMfcd5MUbuo1bhxps3DO2K8AAdzpfWYAvLc8nN6TRquXfB63tPuXYw5Fbc4+rbjSTEBGHVjCFQCUDGd2fwz88OM1igTsHkZhlvwLoP7Irpg6HTqLDnWCEyvuN2TdR5dQ3xxcwRsQCAf245xBJEcht1zdWUeW/QadS4bUwcACD9i6PMepPbu32stYrji8MXcOJihYtHQ+6CgbcCpgyKwTM3DwYAvLrnJB7/+BeWnZPL1d9D1V0CbwCID/fHQ5OtJefL/vsLfs0vc/GIiJp3/3V9EajT4OdzZdj0A7s4k3swK7SdWH23jY1HmL8WJwsr8f5+zg1yb1dEBuC6hEiIIrAq6zdXD4fcBANvhcwYHosn/zQQggC8+c1ppLydjSqjydXDIi8mlRICnX87sYbuGBuP8X0jUFNrwd/f+gHlNeyMS51TWIAOd197BQDg6S2HcbHc4OIREV2eXGqu4DKkAJ1Gnhurso6wCS25vdRJ/SAIwH9/zMP+U8WuHg65AQbeCpo9Og6rZ10FrUaFrb9cwB9f/hpH8stdPSzyUlJGw0ctQBDcK/BWqQSkzxyKrsF6nCysxL0ZOSxVpE7r9rHxSIgJwqWqWiz+kL0JqPOTm6sp/KHsrFE9cEVkAAorjHjqs8OK/iwipQ3oGoS/2JYX/d+mH1Ft5Nan1DIG3gqbOigG78wfhaggHY4VVOCG1Xvw6u4TdmW/RB1BClTdLdstCfXXYs2tidD7qPDVkYt4ZBMDGuqctBoVnv/zEPioBWw7dAHvfMctlKhzM5is7w86H2UvC3UaNVbcPAiCAHyQfRZfHr6g6M8jUtrDk/sjMlCH4xcrsfTTX3hdQi1i4N0BhvcMxZZ7x8mlsk99dhh/WvM1cnIvuXpo5EWk5mo+CnWt7QhDY0OwetZVUKsEbPrhLB7/+Bc2sKJOaUDXIKROsvYmWPLJLyxDpE7NUGsLvDVqxX9WYlwobh9jbUx1f+YBnC6qVPxnEimli78Wq2YMhSAA735/But2nnD1kKgTc98rcDcTFqDDxttGYPlNgxCktzbeuXHNXtz5xn6Wn1OHMNtKCdUKda3tKBMTovDMTYPk/gkPvH+QZefUKS2Y0AtTB0Wj1izirjezcayAnW+pczKYrCWyOk3HXBY+MqU/hvUIQVmNCfP+sx/FlVzvTe7rd33C8Y/rBwAAVnz+K9bsOMbMNzWJgXcHUqkE3DKyB754YAKmJ3aHSgC2HbqAyem7cMfG77H76EVOVFKM0vu0dqQ/D49F+syhUKsEfJhzDn999Vs2saJORxAEPDt9CAZ2C0JRpRGz/v0NThYyu0edj9HUcRlvwLocY+1fExEdpMexggrMee1blFaxaSa5rzt+F4+7r7E2D1z5+RE8suknNlWmRtz/CtwNRQbq8dyfh2DrovGYMjAaALD91wLM3vAdrlu1E6u/OoYzxVUuHiV5GpO5rrmaJ7hhaDe8Omc4AnQafHeyGH94aTf2HC109bCI7PjrNHjjjlHoFxWIgnIDblrzNcvOqdOR1nhrOyjjDQDRwXq8NX8Uwvy1+PlcGW5et5fXPuTWHpzcD//4wwAIApC5/wymvbQHe4/zuoTqMPB2oT5RgVh7ayK+evBq3DamJ/y1ahy/WIlntx7BuJVf4ea1e7FmxzEczitjJpzaTepa667N1ZpyTf9IfHz3WFwRGYALZQbcuuFb/N8HP3KbGupUQv21eGv+KAzuHoxLVbWY9eq32Pj1SfYnoE7BZLbIPUA6qtRcckVkAN7+2yg58/3Hl/fg85/zOnQMRM4073fxeGuetany8YuVmPXvb/G3N/azrxMBYODdKcSH+2PJH6/EN49OxMrpgzGmdxgEAcg+fQkrPz+CKS/uxphntuPB9w/i3e9ycayggoE4OUxurqbgPq2u0DsiAB8vHIu5SXEArJ8yj1v5FVZ/dYxbe1CnERGow7t3jsbkK6NgNFmw5NNDmPv6dzhxkeu+ybWM9XpkKN3VvCn9o4Pw8d1j5Q+mFrz1A+7JyMG5kuoOHwuRM4y9IhxbF43H3KQ4qFUCsg5dwI1r9mL62r1497tcLqvwYoLoIRFcWVkZgoODUVpaiqCgIFcPp93ySqvxxaEL+OrIRew9XoiaWvvmUV38fDC4ewgSYoKQEBOIATFBiA/3h8bDgipXcefnU3Nj33e8CLf8+xv0iQxAVuoEF45QOd+dLMbjH/+MX20NC7v4+eCWkT0wOykOMcG+Lh6de/LEueBKoijizW9O45+fHYbBZIFGJeDW0XG4c3wvdA3hc7Sz64zPqdZoadyXKo0Y9mQWAODYP6e47DrCaLLghS9+w7qdxyGK1uz7LSN74PaxPREX5u+SMVHz3HUuAB079mMF5Vi74wQ+OXhO7rXjoxaQ1Dsc464Ix9grwtE/OhAqD6pG9DaOPJ8YeLuBmlozvj1ZjO9OFmH/qUs4eLakUSAOWNdmxYf5Iy7MDz3Dbf+G+aNHqB+igvQdunbL3bnz86m5se8+ehGzN3yH/tGB+HzReBeOUFlmi4iPD5zDC1/8hjPF1oyJSgBG9wrDHwZ3xeQroxAWoHPxKN2HJ86FzuD4xQo89d9D+OrIRQCARiVg6qAYTE/sjjG9w/ghaifVmZ9TLWlp3OdKqjH2me3QqlX47Z9TXDTCOj+fK8VTnx3CNyesvRAEARjfJwLThnRF8pVRCNL7uHiEBLjvXABcM/YLZTXY9MNZfHLgvJwckATqNRjYNRiDugdjYLdgXBERgJ7hfvDTajpkbNQ+jjyf2vTOvmbNGsTHx0Ov1yMxMRG7d+9u8fydO3ciMTERer0evXr1wrp16xqds2nTJgwYMAA6nQ4DBgzAhx9+2JaheSS9jxoT+kbgocn9kXlXEn58YjI+WjgWT/1pIP46qgeu6hECf60aRpMFRy6UY9uhC1i/6wQWf/gz/vrqtxi38iv0fex/GP5UFq7/127M2/g9Hv3wJ7z05VFkfp+Lbb/kY/+pYhy/WIFLlUauO3QyJeZLW9Q1V/PsC3q1SsBNV3XHjgevwbpbEzEqPhQWEdh7vAiPfvgTEp/6AlNe3I2n/nsI237Jx/mSai7d6CCdZS50Br0jAvD67SPx5ryRGN0rFCaLiE8Onsec177D6OVf4uEPDuLjA+dQUFbj6qGSAjrTXJB6YgT5do6AdmC3YGT8bTTenDcSE/pGQBSBnb9dxIPvH8RVy7Iwfe1ePL/tCHb9dhGFFdzNwhN0pvmglKggPVKuvgKfLxqPL1LH47HrE3BNvwj4adUorzFh34kirN91Avdm5GDqv3ZjwONbMfrpL3HL+m/w8AcHsWrbEbz97Wl8efgCfj5XiovlBpi4larbcfijlMzMTCxatAhr1qzB2LFj8corr2DKlCk4dOgQevTo0ej8kydPYurUqfjb3/6Gt956C19//TVSUlIQERGBm2++GQCwb98+zJw5E08++SRuvPFGfPjhh5gxYwb27NmDUaNGtf+39DBajQpDY0MwNDZEPmaxiDh7qRoniypxuqgSJwsrcbqoCqeKKnG2uBpGswWFFUYUVhjxy/myFh9fJQAhflp08fNBiJ8WgXoNAnQaBOp95P+2fq2xfe2DAL0Gflo19Bo19FoV9D5q+PqoPT7Iuxwl5ktbSWu8NR7S1fxy1CoBvx8Yjd8PjEZuURX++9N5fPZjHn45X4bDedbbq3tOArCWpF/ZNRg9w/0QF+qPHmF+iO3ih8ggHbr4aT2qIZ2rdKa50JmM6xOBcX0i8NPZUrz7fS62/JSHwgoj3tt/Fu/tPwsA6BbiKy8p6h0ZgO5dfNG9ix8iAnQsT3RDnW0ulFZb15sG+3ae7JogCPLcOHGxAp8cPI///piHYwUV2H/6EvafrmtUFRGoQ//oQMSF+aF7F+trd7cuvggP0CLUX8usYSfX2eZDR7giMhBXRAZi/rheqDVb8NuFcvx8rhQ/nSvFL+fLcLKwEiVVtcgvq0F+WQ32nWj+sQL1GoT4+aCLnxbBvtZ/Q/x8EKDTwF+ngb9WDT+dBv5aDfx0auu/WjX8ddZ/dRoVdBo1tBoVr3U6gMOl5qNGjcJVV12FtWvXyscSEhLwpz/9CcuXL290/v/93//hk08+weHDh+VjCxYswMGDB7Fv3z4AwMyZM1FWVob//e9/8jm///3v0aVLF2RkZLRqXO5c8qI0URRRXGlEflkNLpTVIL/UgPzSauSX1aCwwojiSiMuVVn/La9x7p6DGpUAvY/aGohrVdBr1PCVA3Q1tGoVdBoVfNQCfNQqaDWqev8K0KrV8NEI0KpVdt/3UQu2+1lvGpUAtUqARi1ArbJ+rVELtuP1vi//q7Kda/25DV9snPV8UmK+XE5zY//fT3n4+9s/YETPLnh/wZg2/07u7mK5AftOFGHvsUIcPFuKoxfK5Q8lmqISrF2pwwN0CPXXWj94sn0A5W/7EMpfa32O2z+Hbc9ptQo+Guu/apUAlSBAJVgvLFUCoBIECA3+be4cAfWep/X/s+nDEAT757X99+oft36hVgmNlqR44lzozGrNFuw7XoQ9xwqx93ghfjlfhubepbVqFcIDtOjibw0wpA9M/XUa+PlYX2t9tdYPQf20auh81PBRqeTXSo3ttdD6Olj32qhWW58R0nPO+q9Ng2OCUHeu9dvWbzT3fSXYzQtnPaZgrTZryBnPqc42F6T3hqt6hGBzytg2/U4dQRRFnCmuxr4Thdh3vAgHz5biVFFls/NDovdRIdTPOk8C9Rr4SYFHvUDE1xaA+Kit1wbya7ftesNHo4KPSiXPk/qvzXav3So0eP1u/Dpf93pu08RrMdC61/WmXtPb8z4hNHmuAL2PqtH7Cd8blHWp0oiTRZU4VViJc5eq667hbdfxRZWGyz73HaVRWa+ttbZgXOdjnQc6H1twXu/6XHqv8FE1fd0tv6fUu0b3aeIaXZoT1nkF+b+lOaS2zSv5v6U5pWrmv6VzVPbXUw1/Rv15KNS/DoP1X+l9rP4xQUC7YwaHPgY0Go3Izs7GI488Ync8OTkZe/fubfI++/btQ3Jyst2xyZMnY8OGDaitrYWPjw/27duH+++/v9E56enpjgyPmiEIAsICdAgL0OHKrsEtnms0WVBSbcSlyloUVxpRUmVEucGE8hoTKmpMqDDUosL2dXmNCRUG6bgJ1bVmVBvNqDGZ5RcDk0W0nmNwbkDvTP+8cSD+OirO6Y+r1HxpyGAwwGCoK7crK2u6oqFWynirvLsKISJQhz8O6Yo/DukKwNpD4bcL5TicV4bc4iqcLqrCmeIqnLlUjUtVRlhEyNUinm56Ync89+chTn/czjYXOjMftQrj+0ZgfN8IAEBZTS1+zSvHofOlOJxXbq1isl2EGc0WnC+twflSlqM7W69wf2x/8GqnP25nnAt/f/sHANZS2M5MEAT0CPNDj7AemDnCmgmtMprwa345jl4ox5niapy9VIWzl6pxrqQaRZVGGE0W1NRynrTXgccnIcRP6/TH7YzzobPo4m/9sOiqHl2a/L7JbEFZjQmXqowoqapFie3fS1VGlFZbr9WrDGZUGk2oNJhQaTSjylh3TPq3ft7BZBFhMppRaTQDYOf1prw6ZziuGxDV5vs7FHgXFhbCbDYjKsr+B0ZFRSE/P7/J++Tn5zd5vslkQmFhIWJiYpo9p7nHBNxzErkDrUaFyEA9IgPb/gYsiiIMJgtqas2oqbXYBeQ1RjOq6x03miyoNVtvhnr/XWsWYTRZYDRbUGuqO2Z/jgVG23kWiwiTxboXqckswmwRYbKIMNuO1X1tvdWnUai0Rqn50tDy5cuxdOnSy47HbNvH21tKzVtL76PG4O4hGNw9pNH3TGYLiquMKCw3oqjSIFeFVNo+TKow2N7QDGYYTGYYzaL8fDWaLfWe39bnqVkUIYoiLCJgEUVYLCJE6b9FQITt33rneMLy8842F9xJkN4HI+NDMTI+1O64yWyRq5YuVRlxqbKueqnKaH3Nra41o8poRk29f+u/FtaaLc1+DdvzUbT+J0RRtP3rkj+Dx+iMc6F7F1+UVNVixojYVv4WnYefVoOrenRpMjgRRRGVRjMuVRpRVGmdIxUGE6qM1tfs6lozKg0mVNkCklqz2Oiao/41ifS6XveaLcJisX+9rnv9tn8Nt9Q7Jp0vjVEer93YFf7DdRKdcT64C41ahVBbtVN7mGzX30aT9V+DyVzvv61f232/3vuIyWx/jW0y112LN3xPkc6tu0a3wGS2zgtzvXljttTNE7v/ts03S4P5ZbE089/NnSPNWzSej63R3uqtNi18aVhuIopio2OXO7/hcUcf0xMnkacQhLry8s7IYrFOYOnFwEfhQFSJ+VJfWloaUlNT5a/LysoQG9v4Amra4K6YfGV0q8dN1je29n4Q1V5ivQs3+Zjd9+v9d73vtHTh1tx9VErWA6PzzAVPoFGr0L2LdU2rq0jPzYaBuRSwo97Xyo1BmcdVeCp0qrnwyd2/Q4BO43E7nwiCIPekiQ113TxxJrtAXZpjzX3f7ni9/77M+0RT5/oqfD3XmeaDt9GoVdCoVfDnZi9yEC7/K334XC8xom/n66RDgXd4eDjUanWjT6EKCgoaffokiY6ObvJ8jUaDsLCwFs9p7jEBTiJqO5VKgAoClP5cQKn50pBOp4NOd/lXTOnFldyLvPZIgXWsHaWzzQVyDum5afvKlUNxG51xLrQ3Y0Ydxz5h1eQZHTYWZ+iM84G8l3UdOaDkPHLoKlyr1SIxMRFZWVl2x7OysjBmTNPNmpKSkhqdv23bNgwfPlxeh9HcOc09JmCdREFBQXY3os5EqflC5G44F4isOBeI6nA+kNcRHfTuu++KPj4+4oYNG8RDhw6JixYtEv39/cVTp06JoiiKjzzyiDh79mz5/BMnToh+fn7i/fffLx46dEjcsGGD6OPjI37wwQfyOV9//bWoVqvFZ555Rjx8+LD4zDPPiBqNRvzmm29aPa7S0lIRgFhaWuror0TUiLOeT0rMl44aO5Eoci4Q1eeM5xTnAnkCvjcQWTnyfHI48BZFUVy9erUYFxcnarVa8aqrrhJ37twpf2/u3LnihAkT7M7fsWOHOGzYMFGr1Yo9e/YU165d2+gx33//fbFfv36ij4+P2L9/f3HTpk0OjYmTiJzJmc8nJeZLR42diHOBqI6znlOcC+Tu+N5AZOXI88nhfbw7K3fYk4/chzs/n9x57NT5uPPzyZ3HTp2Tuz6n3HXc1Hm583PKncdOnY8jzyd2WiIiIiIiIiJSEANvIiIiIiIiIgW1aR/vzkiqmC8rK3PxSMgTSM8jd1yJwblAzsS5QFTHXecD5wI5m7vOBYDzgZzLkbngMYF3eXk5AHAvb3Kq8vJyBAcHu3oYDuFcICVwLhDVcbf5wLlASnG3uQBwPpAyWjMXPKa5msViwfnz5xEYGAhBqNv4vKysDLGxsThz5gwbKCjI0/7OoiiivLwcXbt2hUrlXisympsLgOf9f+qsPOnvzLlA7eFpf2d3nQ+cC67naX9nd50LAGMGV/O0v7Mjc8FjMt4qlQrdu3dv9vtBQUEe8T+3s/Okv7O7fYIrudxcADzr/1Nn5il/Z84Fai9P+ju743zgXOg8POnv7I5zAWDM0Fl40t+5tXPBvT6iIiIiIiIiInIzDLyJiIiIiIiIFOTxgbdOp8MTTzwBnU7n6qF4NP6d3QP/P3UM/p07P/4/6hj8O3d+/H/UMfh37vz4/6hjePPf2WOaqxERERERERF1Rh6f8SYiIiIiIiJyJQbeRERERERERApi4E1ERERERESkIAbeRERERERERAry+MB7zZo1iI+Ph16vR2JiInbv3u3qIbm1Xbt2Ydq0aejatSsEQcBHH31k931RFLFkyRJ07doVvr6+uPrqq/HLL7+4ZrBkh3PBuTgX3BfngvNxPrgvzgfn4lxwX5wLzsW50JhHB96ZmZlYtGgRFi9ejJycHIwbNw5TpkxBbm6uq4fmtiorKzFkyBC8/PLLTX5/5cqVWLVqFV5++WV8//33iI6OxqRJk1BeXt7BI6X6OBecj3PBPXEuKIPzwT1xPjgf54J74lxwPs6FJogebOTIkeKCBQvsjvXv31985JFHXDQizwJA/PDDD+WvLRaLGB0dLT7zzDPysZqaGjE4OFhct26dC0ZIEs4FZXEuuA/OBeVxPrgPzgdlcS64D84FZXEuWHlsxttoNCI7OxvJycl2x5OTk7F3714XjcqznTx5Evn5+XZ/c51OhwkTJvBv7kKcCx2Pc6Fz4lxwDc6HzonzoeNxLnROnAsdz1vngscG3oWFhTCbzYiKirI7HhUVhfz8fBeNyrNJf1f+zTsXzoWOx7nQOXEuuAbnQ+fE+dDxOBc6J86Fjuetc8FjA2+JIAh2X4ui2OgYORf/5p0T/790PP7NOyf+f3EN/t07J/5/6Xj8m3dO/P/S8bztb+6xgXd4eDjUanWjT00KCgoafbpCzhEdHQ0A/Jt3MpwLHY9zoXPiXHANzofOifOh43EudE6cCx3PW+eCxwbeWq0WiYmJyMrKsjuelZWFMWPGuGhUni0+Ph7R0dF2f3Oj0YidO3fyb+5CnAsdj3Ohc+JccA3Oh86J86HjcS50TpwLHc9b54LG1QNQUmpqKmbPno3hw4cjKSkJ69evR25uLhYsWODqobmtiooKHDt2TP765MmTOHDgAEJDQ9GjRw8sWrQITz/9NPr06YM+ffrg6aefhp+fH2bNmuXCURPngvNxLrgnzgVlcD64J84H5+NccE+cC87HudAEV7VT7yirV68W4+LiRK1WK1511VXizp07XT0kt/bVV1+JABrd5s6dK4qidXuAJ554QoyOjhZ1Op04fvx48aeffnLtoEkURc4FZ+NccF+cC87H+eC+OB+ci3PBfXEuOBfnQmOCKIpih0X5RERERERERF7GY9d4ExEREREREXUGDLyJiIiIiIiIFMTAm4iIiIiIiEhBDLyJiIiIiIiIFMTAm4iIiIiIiEhBDLyJiIiIiIiIFMTAm4iIiIiIiEhBDLyJiIiIiIiIFMTAm4iIiIiIiEhBDLyJiIiIiIiIFMTAm4iIiIiIiEhBDLyJiIiIiIiIFPT/WahrJSWMGC4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x300 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import cauchy, gamma, expon, lognorm, loglaplace\n",
    "\n",
    "# Define the x range for the plot\n",
    "x = np.linspace(0, 15, 1000)\n",
    "\n",
    "# Create a single row with 5 subplots\n",
    "fig, axs = plt.subplots(1, 5, figsize=(10, 3))\n",
    "\n",
    "# Cauchy Distribution\n",
    "cauchy_pdf = cauchy.pdf(x)\n",
    "axs[0].plot(x, cauchy_pdf)\n",
    "axs[0].set_title(\"Cauchy\")\n",
    "\n",
    "# Gamma Distribution\n",
    "gamma_pdf = gamma.pdf(x, a=2)\n",
    "axs[1].plot(x, gamma_pdf)\n",
    "axs[1].set_title(\"Gamma\")\n",
    "\n",
    "# Exponential Distribution\n",
    "exponential_pdf = expon.pdf(x)\n",
    "axs[2].plot(x, exponential_pdf)\n",
    "axs[2].set_title(\"Exponential\")\n",
    "\n",
    "# Log Normal Distribution\n",
    "log_normal_pdf = lognorm.pdf(x, s=0.7)\n",
    "axs[3].plot(x, log_normal_pdf)\n",
    "axs[3].set_title(\"Log Normal\")\n",
    "\n",
    "# Log Logistic Distribution\n",
    "log_logistic_pdf = loglaplace.pdf(x, c=0.5)\n",
    "axs[4].plot(x, log_logistic_pdf)\n",
    "axs[4].set_title(\"Log Logistic\")\n",
    "\n",
    "# Adjust layout\n",
    "plt.tight_layout()\n",
    "\n",
    "# Display the plot\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy  \n",
    "from scipy import stats \n",
    "import scipy.optimize as optimize"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Các hàm liên quan"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1 Hàm phân phối gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "class opt_gamma:\n",
    "\n",
    "    def __init__(self, actual_pd):\n",
    "        self.actual_pd = actual_pd\n",
    "        self.x_input = range(1, len(actual_pd) + 1)\n",
    "\n",
    "    def func(self, x, scale, a, b):        \n",
    "        predict = stats.gamma.pdf(x, a = a, scale = b) * scale\n",
    "        return predict\n",
    "\n",
    "    def sse(self, params, xobs, yobs):\n",
    "        ynew = self.func(xobs, *params)\n",
    "        mse = np.sum((ynew - yobs) ** 2)        \n",
    "        return mse\n",
    "    \n",
    "    def solver(self):\n",
    "        p0 = [1,1,1]\n",
    "        bounds = [(0.0001, 2), (0.0001, 10), (0.0001, 10)]\n",
    "        res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), bounds= bounds)       \n",
    "        # res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), method='Nelder-Mead')        \n",
    "        return res\n",
    "    \n",
    "    def predict(self, t=30):\n",
    "        res = self.solver()\n",
    "        ypred = self.func(range(1, t+1), *res.x)        \n",
    "        return ypred        "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2 Hàm Phân phối mũ Exponential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class opt_exponential:\n",
    "\n",
    "    def __init__(self, actual_pd):\n",
    "        self.actual_pd = actual_pd\n",
    "        self.x_input = range(1, len(actual_pd) + 1)\n",
    "\n",
    "    def func(self, x, scale, a):        \n",
    "        # change function here        \n",
    "        pred = stats.expon.pdf(x, scale = 1/a) * scale\n",
    "        return pred\n",
    "\n",
    "    def sse(self, params, xobs, yobs):\n",
    "        ynew = self.func(xobs, *params)\n",
    "        mse = np.sum((ynew - yobs) ** 2)        \n",
    "        return mse\n",
    "    \n",
    "    def solver(self):\n",
    "        # change initial guess here\n",
    "        p0 = [1,2]\n",
    "        bounds = [(0.0001, 2), (0.0001, 0.5)]           \n",
    "        res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), bounds = bounds)             \n",
    "        return res\n",
    "    \n",
    "    def predict(self, t=30):\n",
    "        res = self.solver()\n",
    "        ypred = self.func(range(1, t+1), *res.x)        \n",
    "        return ypred       "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3 Hàm phân phối Cauchy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class opt_cauchy:\n",
    "\n",
    "    def __init__(self, actual_pd):\n",
    "        self.actual_pd = actual_pd\n",
    "        self.x_input = range(1, len(actual_pd) + 1)\n",
    "\n",
    "    def func(self, x, scale, a, b):        \n",
    "        # change function here        \n",
    "        predict = stats.cauchy.pdf(x, loc = b, scale = a) * scale\n",
    "        return predict\n",
    "\n",
    "    def sse(self, params, xobs, yobs):\n",
    "        ynew = self.func(xobs, *params)\n",
    "        mse = np.sum((ynew - yobs) ** 2)        \n",
    "        return mse\n",
    "    \n",
    "    def solver(self):\n",
    "        # change initial guess here\n",
    "        p0 = [1,1,1]\n",
    "        bounds = [(0.0001, 2), (0.0001, 15), (-30, 30)]\n",
    "        res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), bounds = bounds)        \n",
    "        return res\n",
    "    \n",
    "    def predict(self, t=30):\n",
    "        res = self.solver()\n",
    "        ypred = self.func(range(1, t+1), *res.x)        \n",
    "        return ypred    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.4 Hàm phân phối log logistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class opt_log_logistic:\n",
    "\n",
    "    def __init__(self, actual_pd):\n",
    "        self.actual_pd = actual_pd\n",
    "        self.x_input = range(1, len(actual_pd) + 1)\n",
    "\n",
    "    def func(self, x, scale, a, b):        \n",
    "        # change function here\n",
    "        predict = scale*(a/b)*(x/b)**(a-1)/(1+(x/b)**a)**2\n",
    "        return predict\n",
    "\n",
    "    def sse(self, params, xobs, yobs):\n",
    "        ynew = self.func(xobs, *params)\n",
    "        mse = np.sum((ynew - yobs) ** 2)        \n",
    "        return mse\n",
    "    \n",
    "    def solver(self):\n",
    "        # change initial guess here\n",
    "        p0 = [2,3,1]\n",
    "        bounds = [(0.0001, 2), (0.0001, 10), (0.0001, 10)]\n",
    "        res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), bounds=bounds)        \n",
    "        return res\n",
    "    \n",
    "    def predict(self, t=30):\n",
    "        res = self.solver()\n",
    "        ypred = self.func(range(1, t+1), *res.x)        \n",
    "        return ypred    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.5 Hàm phân phối log normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class opt_log_normal:\n",
    "\n",
    "    def __init__(self, actual_pd):\n",
    "        self.actual_pd = actual_pd\n",
    "        self.x_input = range(1, len(actual_pd) + 1)\n",
    "\n",
    "    def func(self, x, scale, a, b):        \n",
    "        # change function here\n",
    "        predict = (np.exp(-(np.log(x) - a)**2 / (2 * b**2)) / (x * b * np.sqrt(2 * np.pi)))*scale\n",
    "        return predict\n",
    "\n",
    "    def sse(self, params, xobs, yobs):\n",
    "        ynew = self.func(xobs, *params)\n",
    "        mse = np.sum((ynew - yobs) ** 2)        \n",
    "        return mse\n",
    "    \n",
    "    def solver(self):\n",
    "        # change initial guess here\n",
    "        p0 = [1,1,1]\n",
    "        bounds = [(0.0001, 2), (0.0001, 0.5), (0.0001, 0.5)]\n",
    "        res = scipy.optimize.minimize(self.sse, p0, args=(self.x_input, self.actual_pd), bounds=bounds)        \n",
    "        return res\n",
    "    \n",
    "    def predict(self, t=30):\n",
    "        res = self.solver()\n",
    "        ypred = self.func(range(1, t+1), *res.x)        \n",
    "        return ypred    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Fit models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Đặt giá trị initial\n",
    "\n",
    "- Gamma = 1,1,1\n",
    "\n",
    "- Cauchy = 1,1,1\n",
    "\n",
    "- Expo = 1,2\n",
    "\n",
    "- Log Logistic = 2,3,1\n",
    "\n",
    "- Log Normal = 1,1,1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Đọc dữ liệu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel('results/Cohort Analysis.xlsb', engine='pyxlsb', sheet_name='1.2 Visualize', usecols = 'O:Z')\n",
    "data.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Fit models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fn_export_by_segment(segment, period):\n",
    "    # segment = segment.upper()\n",
    "    actual_pd = data[data.Segment_Group.isin([segment])]\n",
    "    actual_pd = actual_pd.iloc[0, 2:]\n",
    "    actual_pd = actual_pd.values\n",
    "    \n",
    "    res_gamma = opt_gamma(actual_pd)\n",
    "    res_exponential = opt_exponential(actual_pd)\n",
    "    res_cauchy = opt_cauchy(actual_pd)\n",
    "    res_log_logistic = opt_log_logistic(actual_pd)\n",
    "    res_log_normal = opt_log_normal(actual_pd)\n",
    "    \n",
    "    params = pd.DataFrame({\n",
    "        'Segment': segment,\n",
    "        'Distribution': ['Gamma', 'Exponential', 'Cauchy', 'Log-Logistic', 'Log-Normal'],\n",
    "        'SSE': [res_gamma.solver().fun, res_exponential.solver().fun, res_cauchy.solver().fun, res_log_logistic.solver().fun, res_log_normal.solver().fun],\n",
    "        'params': [res_gamma.solver().x, res_exponential.solver().x, res_cauchy.solver().x, res_log_logistic.solver().x, res_log_normal.solver().x]\n",
    "    })    \n",
    "   \n",
    "    params[['scale', 'a', 'b']] = pd.DataFrame(params.params.to_list())    \n",
    "    del params['params']\n",
    "    \n",
    "    dict_predict = {\n",
    "            'Segment': segment,\n",
    "            'Period': range(1, period+1),\n",
    "            'Gamma': res_gamma.predict(period),\n",
    "            'Exponential': res_exponential.predict(period),\n",
    "            'Cauchy': res_cauchy.predict(period),\n",
    "            'Log-Logistic': res_log_logistic.predict(period),\n",
    "            'Log-Normal': res_log_normal.predict(period)\n",
    "        }\n",
    "    \n",
    "    best_distribution = params.sort_values('SSE').head(1)\n",
    "    best_distribution['Tag best distribution'] = 'Best distribution'    \n",
    "\n",
    "    print('Best distribution of ' + segment + ' is '+ best_distribution['Distribution'].item())\n",
    "\n",
    "    return params, best_distribution, pd.DataFrame(dict_predict)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_best_fit = []\n",
    "params_frame = []\n",
    "pred_all_curve = []\n",
    "for seg in data.Segment_Group:\n",
    "    params, df_best_fit, df_all = fn_export_by_segment(seg, 30)\n",
    "    pred_best_fit.append(df_best_fit)\n",
    "    params_frame.append(params)\n",
    "    pred_all_curve.append(df_all)\n",
    "\n",
    "pred_best_fit = pd.concat(pred_best_fit, axis=0)\n",
    "params_frame = pd.concat(params_frame, axis=0)\n",
    "pred_all_curve = pd.concat(pred_all_curve, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "pivot_pred_all_curve = pd.melt(pred_all_curve, \n",
    "        id_vars=['Segment', 'Period'], \n",
    "        value_vars=['Gamma', 'Exponential', 'Cauchy', 'Log-Logistic', 'Log-Normal'],\n",
    "        var_name = 'Distribution').pivot(\n",
    "                index = ['Segment', 'Distribution'], columns='Period'\n",
    "        )\n",
    "pivot_pred_all_curve.reset_index(inplace=True, drop=False)\n",
    "pivot_pred_all_curve.columns = ['Segment', 'Distribution', *range(1,31)]\n",
    "pivot_pred_all_curve = params_frame.merge(pivot_pred_all_curve, how='inner', on=['Segment', 'Distribution']).merge(\n",
    "        pred_best_fit[['Segment', 'Distribution', 'Tag best distribution']], how='left', on=['Segment', 'Distribution']\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xuất OutPut ra Excel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_best_fit.reset_index(inplace=True, drop=True)\n",
    "params_frame.reset_index(inplace=True, drop=True)\n",
    "pivot_pred_all_curve.reset_index(inplace=True, drop=True)\n",
    "with pd.ExcelWriter('cohort_fitted.xlsx') as writer:\n",
    "    pred_best_fit.to_excel(writer, sheet_name='pred_best_fit')\n",
    "    params_frame.to_excel(writer, sheet_name='params_frame')        \n",
    "    pivot_pred_all_curve.to_excel(writer, sheet_name='pivot_pred_all_curve')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rdm",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "b0fc7d2c5e3b32912b68ef575c1bdd2ee7ee3f28f2d636d660b71751dc8fb34a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}