{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script>\n",
    "        var code_show=false; //true -> hide code at first\n",
    "    \n",
    "        function code_toggle() {\n",
    "            $('div.prompt').hide(); // always hide prompt\n",
    "            $('div[style*=\\\"border-style: dashed\"\\]').hide(); // hide extra comments, usually log from treasuredata\n",
    "            $('div.output_stderr').parent().hide();\n",
    "    \n",
    "            if (code_show){\n",
    "                $('div.input').hide();\n",
    "            } else {\n",
    "                $('div.input').show();\n",
    "            }\n",
    "            code_show = !code_show\n",
    "        }\n",
    "        $( document ).ready(code_toggle);\n",
    "    </script>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "plt.style.use(\"bmh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "## Background\n",
    "\n",
    "On weekdays, I commute using either train or TransJakarta. Since those are public transportations, we can't expect the vehicle has enough space for all passengers (yes, I don't even mention *comfortable*). Thus, sometimes I have to wait for the next train or bus, wishing that there is more space to get in. The question is, will the train/bus be spacious enough so that I could get in?\n",
    "\n",
    "## Research Questions\n",
    "\n",
    "How likely will I get in a bus after waiting for certain minutes?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assumptions and Known Limitations\n",
    "\n",
    "- This analysis focuses on bus, not train, assuming that I don't have any other commuting options if I am already waiting at the train station. Meanwhile, if the bus doesn't arrive, I can use ride-hailing service.\n",
    "- The survival analysis is based on Kaplan-Meier estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Preparation\n",
    "\n",
    "The dataset contains 1,000 observations which simulate bus waiting time, which ranges from 1 to 45 minutes. `observed` attribute describes whether I get in the bus or not after waiting for `wait_time` minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import randint\n",
    "import datetime\n",
    "\n",
    "n = 1000 # number of observations\n",
    "start_date = datetime.datetime(2018,1,1,8,0,0)\n",
    "wait_time = []\n",
    "end_date = []\n",
    "\n",
    "for i in range(n):\n",
    "    date = start_date + datetime.timedelta(minutes=randint(1,45))\n",
    "    wait_time.append((date.hour-8)*60 + date.minute)\n",
    "    end_date.append(date.strftime(format='%H:%M'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(\n",
    "    data = {\n",
    "        'start_time': np.repeat(start_date.strftime(format='%H:%M'), n),\n",
    "        'end_time': end_date,\n",
    "        'wait_time': np.array(wait_time).T,\n",
    "        'observed': np.random.binomial(n=1, p=.5, size=n)\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset preview:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>start_time</th>\n",
       "      <th>end_time</th>\n",
       "      <th>wait_time</th>\n",
       "      <th>observed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>995</th>\n",
       "      <td>08:00</td>\n",
       "      <td>08:20</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>08:00</td>\n",
       "      <td>08:15</td>\n",
       "      <td>15</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>08:00</td>\n",
       "      <td>08:42</td>\n",
       "      <td>42</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>08:00</td>\n",
       "      <td>08:11</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>999</th>\n",
       "      <td>08:00</td>\n",
       "      <td>08:18</td>\n",
       "      <td>18</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    start_time end_time  wait_time  observed\n",
       "995      08:00    08:20         20         1\n",
       "996      08:00    08:15         15         1\n",
       "997      08:00    08:42         42         1\n",
       "998      08:00    08:11         11         1\n",
       "999      08:00    08:18         18         0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"Dataset preview:\")\n",
    "df.tail(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis\n",
    "\n",
    "We use Kaplan-Meier estimator to calculate probability of getting a bus after waiting for *x* minutes. Wait, how does it work?\n",
    "\n",
    "We aggregate the dataset based on the time attribute (`wait_time`), then we calculate number of `observed` events (get in bus). There are cases when we haven't observed the final result due to **censoring**: at the time of data collection, the `observed` event doesn't occur yet and we lost track of the final outcome. For example, I haven't got in the bus after waiting for 30 minutes, and we don't know whether I continue waiting until I get in the bus or decide to use other means of transportation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from lifelines.plotting import plot_lifetimes\n",
    "\n",
    "# plt.title(\"Observed vs censored events\")\n",
    "# plot_lifetimes(df['wait_time'].head(10), df['observed'].head(10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct survival table from the aggregated data, then we use survival function $S(t)$ to calculate the survival probability.\n",
    "\n",
    "<br/>\n",
    "<center> $S(t) = Pr(T > t)$ </center>\n",
    "\n",
    "> In plain English: the survival function defines the probability the death event has not occured yet at time $t$, or equivalently, the probability of surviving past time $t$. - [Lifelines Documentation](https://lifelines.readthedocs.io/en/latest/Survival%20Analysis%20intro.html#survival-function)\n",
    "\n",
    "The `observed` event on our case is \"getting in the bus\", then the survival function calculates probability of **not** getting in the bus at time $t$; it is displayed on `KM_estimate` field below. To make it more intuitive, we can use $1 - \\text{KM_estimate}$ to resemble probability of getting in the bus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import KaplanMeierFitter\n",
    "kmf = KaplanMeierFitter()\n",
    "kmf.fit(durations=df['wait_time'], \n",
    "        event_observed=df['observed']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Survival table\n",
      "==============\n",
      "\n",
      "          removed  observed  censored  entrance  at_risk  KM_estimate\n",
      "event_at                                                             \n",
      "0.0             0         0         0      1000     1000     1.000000\n",
      "1.0            33        23        10         0     1000     0.977000\n",
      "2.0            22        10        12         0      967     0.966897\n",
      "3.0            26        13        13         0      945     0.953595\n",
      "4.0            19         7        12         0      919     0.946332\n",
      "5.0            19        10         9         0      900     0.935817\n",
      "6.0            22         9        13         0      881     0.926257\n",
      "7.0            21         9        12         0      859     0.916552\n",
      "8.0            20        11         9         0      838     0.904521\n",
      "9.0            18         9         9         0      818     0.894569\n"
     ]
    }
   ],
   "source": [
    "from lifelines.utils import survival_table_from_events\n",
    "\n",
    "table = survival_table_from_events(df['wait_time'], df['observed'])\n",
    "table = table.merge(\n",
    "    kmf.survival_function_,\n",
    "    how='left',\n",
    "    left_index=True,\n",
    "    right_index=True\n",
    ")\n",
    "print(\"Survival table\")\n",
    "print(\"==============\\n\")\n",
    "print(table.head(10))\n",
    "\n",
    "table['prob_get_in_bus'] = 1-table['KM_estimate'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median: 34 minutes\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAFzCAYAAACZwbV4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4XGXZ+PHvk0km+540XdOW0oUWqEixdKEsLRTEooALvMoiKP60oIi8ihsCboAKKooLigjqC4qAFSgFpGVrWQsttJTubdKQPZ0kk2WSyfP7Y9IY0kk66cyZM/fk/lwXl+lkOrmHr0OfnjnzHGOtRSmllFJKxU+K2wMopZRSSo00ugBTSimllIozXYAppZRSSsWZLsCUUkoppeJMF2BKKaWUUnGmCzCllFJKqThLdXuA4VqzZo1NT0939GdYazHGOPozlLO0oWzaTz5pDU2gHgDrLXF5ksQhrWGiamtrq1+8eHHpwNvFLcDS09OZMWOGoz+jvr6ekhJ9EUqmDWXTfvJJa1i4sgiAprMaXZ4kcUhrmKjWr1+/J9zt+hZkGH6/3+0RVJS0oWzaTz5tKJ82dJYuwMIYPXq02yOoKGlD2bSffNpQPm3oLF2AhVFdXe32CCpK2lA27SefNpRPGzpL3Dlg4VhraW1tJVbXtUxPT6e5uTkmj6WGzxhDTk5OVCd/pqWlxXAiFW/aTz5tKJ82dFZSLMBaW1tJT0/H6/XG5PGys7PxeDwxeSw1fIFAgNbWVnJzcw/7MfLz82M4kYo37SefNpRPGzorKd6CtNbGbPEF0N3dHbPHUsPn9XqjPppZX18fo2mUG7SffNpQPm3orKQ4AhZrevRLPv2bm2zaTz5pDXX7iYNJayhNUhwBi7VYnUum3BMIBNweQUVB+8mnDeXThs7SBVgYPT09cf15y5Yt44033nDs8W+77bZD3mfp0qWH9dizZ8+moaHhsH6vk9rb290eQUVB+8mnDeXThs7SBVgYTnzyIxgMxvwxI3X77bcf8j6rVq2KwyTxo/vXyKb95JPWMPfFU8l98VS3x0go0hpKowuwMLq6uoZ1/7179zJ37lyuuOIK5s6dyyWXXEJbWxuzZ8/mhhtu4JRTTuGRRx7hrbfe4vTTT2fhwoVcdNFF7N+/v+8xHnjgARYtWsT8+fN5/fXXB/1Z9fX1nHvuucybN48vf/nLHHvssX1HoP7+97+zZMkSFi1axFe/+lWCwSA33ngj7e3tLFq0iCuuuGLQx50wYQIAL7zwAsuWLeOSSy7pe06Hekv2l7/8JQsWLGDJkiXs3LkTgOXLl/Ovf/3roMevrq7m7LPP7nuu69atO8S/3cOj+9fIpv3kk9YwtXkDqc0b3B4joUhrKE3SnYR/xh+ceSvvyc8dN+T3t23bxi9+8QtOPPFErrzySv74xz8CUFRUxJo1awBYuHAht9xyCwsWLOBHP/oRt9xyCz/+8Y+B0KHe5557jrVr13LVVVexdu3asD/n1ltv7VtgPf300/zlL38B4N133+Xhhx9m5cqVpKWlce211/KPf/yD733ve/zhD3/gueeei/i5bty4kbVr1zJmzBjOPPNMXn75ZU488cRB75+Xl8eLL77I/fffz7e+9S3uv//+Qe/74IMPctppp/G1r32NYDBIW1tbxHMNRyw/FaviT/vJpw3lS+aGdf4A6Z4U8jLcWwbpEbAYGTduXN8i5ZOf/CQvv/wyAOeeey4Azc3N+Hw+FixYAMCFF174vqM/559/PgDz58+npaUFn88X9ue89NJLnHfeeQAsWbKEgoICAJ577jk2bNjA4sWLWbRoEc899xy7d+8+rOfywQ9+kHHjxpGSksIxxxzD3r17h7z/gdnPP/98Xn311UM+9t/+9jduvvlmNm/eHNVeX0Nx6nFVfGg/+bShfMnc8DfrKrn4gU28UhH+z9p4SLojYIc6UhWJzs5O0tPTh/V7Bu7afuDXWVlZUf3+SFlrueCCC7j++uuH9fvC6f/cPR7PIfdF6z/rga89Hk/fhxl6enr6Pk0zf/58Hn30UZ588kmWL1/Ol770JS644IKoZx6ooaGBnJycmD+uig/tJ582lC9ZG26tb+OF3T68HsOU4sj+jHaCHgELIzV1+OvSyspKXnnlFSD0NtvcuXPf9/28vDwKCgr6jno98MADzJ8/v+/7Dz/8MBA6wpWXl0deXl7YnzN37lweeeQRAJ555pm+88gWLVrEihUrqKurA6CpqYmKioq+5zPc89qG48DsDz/8MCeccAIA5eXlbNgQOp9i5cqVfT+/oqKCUaNGcckll3DRRRf13SfWCgsLHXlcFR/aTz5tKF+yNvzza+8BcM7MUoqz3LvcUtIdAYuFnp6eYW/GOnXqVP74xz9y1VVXMX36dC677DLuuuuu993nzjvv5JprrqG9vZ1Jkybxq1/9qu976enpnHzyyXR1dXHHHXcM+nO+/vWv8/nPf54HHniAE044gbKyMnJyciguLuZb3/oW559/Pj09PaSlpXHrrbcyYcIELrnkEhYuXMjs2bP5/e9/P7x/GRHYv38/CxcuJD09ve85X3zxxXz605/mpJNOYvHixWRnZwOhk/zvuOMO0tLSyM7O5je/+U3M54HQOXWDLWJV4tN+8mlD+ZKx4abqVl6tbCYrLYVPzS5zdRYjbdPRdevW2RkzZrzvtubm5pj+n2S4b0Hu3buXCy64YNAT52Ops7MTj8dDamoqr7zyCtdee+2wTrCXItqme/bsYeLEiTGcSMWT9pNPWsOst68GoO3on7s8SeKQ1vBQrLX872Pb2VjdymeOG83Fx4+Jy89dv37964sXL54z8HY9AhZGIl8BvrKykssuu4yenh68Xi8//7n+xyIc3b9GNu0nn7SGuvA6mLSGh/JGVQsbq1vJTfdw/jGj3B5HF2DhdHV1DesIWHl5ecyPfv31r3/ld7/73ftumzt3Lj/5yU949tlnD+sxGxsb+djHPnbQ7Y888ghFRUVD/t6LLrqIPXv2vO+2733veyxevPiwZnFadXV1Uv3NbaTRfvJpQ/mSqaG1lj/1nvv1iWNHke11/5rPugALIyXF/c8mfPrTn+bTn/50TB+zqKjosN+uvO+++2I6i9MyMjLcHkFFQfvJJ62hx/cmAMH8D7g8SeKQ1nAoL+1t5t26NgoyUvnozFK3xwF0ARZWIizAVHQyMzPdHkFFQfvJJ61h3trTAGg6q9HlSRKHtIaD6bGWP79eBcCFHygjM839o1+g21CEdah9r1Tia2pqcnsEFQXtJ582lC9ZGj63cz87GzsoyU7j7Bklbo/TRxdgYRzOPmAqsRQXF7s9goqC9pNPG8qXDA2DPZZ714fO/frMcaPxpibOsidxJkkgwWDQtZ+9bNky3ngjdD3LT37yk4NekkgNraWlxe0RVBS0n3zaUL5kaPj09kYqfZ2MzfNyxrTEWlDqoZ4wEmVvtL///e9ujyDWgUsfKZm0n3zaUD7pDbuCPfxlfTUAnzluDKkpw7vEn9P0CFgYw90HbO/evcydO5fly5dzwgkncMUVV7BmzRrOPPNM5syZw+uvv47f7+fKK69kyZIlnHzyyTz++ONAaKfhyy+/nLlz53LRRRfR3t7e97izZ8+moaEBgM985jOceuqpzJs3j3vuuafvPhMmTOAHP/gBJ510Eqeffjq1tbXR/wtIAsm2f81Io/3k04bySW+48t0GaloDTCzI4NQpiXdZJV2AhXE4103cuXMny5cv5+WXX2bbtm08+OCDrFy5kptuuonbb7+d2267jUWLFvH000+zYsUKvve97+H3+7n77rvJzMzk5Zdf5rrrrhv02oh33HEHq1ev5plnnuH3v/89jY2hT+r4/X7mzJnD888/z/z587n33nujeu7Jorq62u0RVBS0n3zaUD7JDTu6e/jbm6H5Lzl+DJ4EO/oFSfoWZOHKwTcV9c+6jUD5pQB4995D9qZrBr3vcD6OPHHiRGbOnAnAjBkzOPnkkzHGMHPmTPbu3UtVVRUrV67su/5jR0cHlZWVrFu3jiuuuAKAWbNmMWvWrLCP/7vf/Y7HHnsMgH379rFjxw6Kiorwer0sXboUCB0xW7NmTcQzJ7Nk+fj0SKX95JPWsHn+M26PkHCkNezv0c11NLZ1c2RxJgsm5bs9TlhJuQBzg9fr7fs6JSWl79cpKSl0d3fj8Xj485//zNSpU4f92C+88ALPPvssq1atIisri2XLltHZ2QmE3i41JrSy93g8rn6AIJH076Hk0X7ySWuoG7AeTFrDA9oCQR7YGDod59I5Y/r+jEw0SbkAi/TIVaD80r6jYf0N92LckTjttNO46667uOWWWzDGsHHjRo499ljmzZvHgw8+yKJFi9i8eTObNm066Pc2NzdTUFBAVlYWW7du5bXXXovpbMnI5/NRUFDg9hjqMGk/+bShfFIbPrSpDl9HN7PKsjlhfJ7b4wxKzwELw4l9wK699lq6urpYuHAh8+bN40c/+hEAl112GX6/n7lz53LzzTcze/bsg37v4sWL6e7uZu7cudx0003MmXPQRdXVACUlibPZnho+7SeftIZZb19N1ttXuz1GQpHWEKC5o5sHN9YAcOnxiXv0C8AkypYLkVq3bp2dMWPG+25rbm4mLy92q9xAICD20GuyiLZpVVUVY8eOjeFEKp60n3zSGh44d1gvRfRf0hoC3P1qFfdvqOG4sbnc8uEj3R4HgPXr17++ePHig46c6BGwMKQtStXBDueTrCpxaD/5tKF80ho2tXXx8KY6IHTuV6LTBVgYw90HTCUe6fvXjHTaTz5tKJ+0hvdvqKGzu4cTy/M4alS22+Mcki7AwpC26lcHk7x/jdJ+yUAbyiepYaWvg0ffqQdC+35JoAuwMFJS9F+LdNnZif+3HzU47SefNpRPSsMea7nt+b109ViWTitiSnGW2yNFJClWGsaYmF6zKpE/NTESBAKBqBt4PJ4YTaPcoP3k04bySWn4+JYG3q72U5iZyhVzx7k9TsSSYh+wnJwcWltb6ejoiMnj+f1+MSv/ZGSMIScnJ6rHaG5uprAw8a79pSKj/eST1rA77+AtgEY6CQ3r/AH+8Mo+AJbPH09uupxljZxJh2CMITc3N2aPl5qaSlaWjEOYKrzS0lK3R1BR0H7ySWvYsmC12yMknERvaK3lly9U0NbVw/yJ+Zw0SdamsUnxFmSsHbjQtZJLG8qm/eTThvIlesM1O/fzckUz2V4PV82fIO70IV2AhaH7gMmnDWXTfvJpQ/kSuWFzRzd3rqsE4IoPjaU4W972UXFbgBljzjTGvGuM2W6MuS7M98uNMauNMW8YYzYaYz4cr9kGSvTDrurQtKFs2k8+aQ0LVxb17YavQhK54W9fqsTX0c3sMTmcOb3Y7XEOS1wWYMYYD/Br4CxgJnChMWbmgLt9B/i7tfY44ALgznjMFk5NTY1bP1rFiDaUTfvJpw3lS9SGr1Y08/T2Jrwew9ULy8W99XhAvI6AfQjYbq3daa0NAPcDHx1wHwscuPhfPlAVp9kOEu0n8JT7tKFs2k8+bShfIjZsCwT5xYt7Abj4+DGMy093eaLDF69PQY4DKvr9uhKYO+A+NwBPGmOuArKBJfEZTSmllFIS/Om196ht7WJqSSbnHz3K7XGikkjbUFwI3GOt/ZkxZh5wnzHmaGttT/871dbWcvnll5OamkowGOS8885j+fLlVFdXk52djcfjobm5mdLSUhobG7HWUlpaSk1NTd9qvrW1lbKyMurq6jDGUFRURF1dHXl5eQSDQaqqqsjJyaG6upq0tDTy8/Opr68nPz+fQCBAe3s7o0ePprq6Gq/XS25uLg0NDRQWFtLe3k5HR0ff9zMyMsjMzKSpqYni4mJaWloIBAJ938/MzMTr9eLz+SgpKcHn89HV1dX3/Vg9J7/f3/eYI+E51dTU0NramlTPKRk7Dfac9u3bR1paWlI9p2TsNNRz8vl8tLa2inlOB3a72rNnz4jqNNRz8vv9eDyehHlOL2ypYMXmFjwGLjwyDX9ri4hOgzHx+JRD74LqBmvt0t5ffxPAWvvjfvfZBJxpra3o/fVO4ERrbW3/x1q3bp2dMWOGo/N2dHSQkZHh6M9QztKGsmk/+aQ1PHACftNZib31QjwlUsNAsIcvPfwue/d3cOHsMj57wli3R4rY+vXrX1+8ePGcgbfH6xywV4GpxpjJxhgvoZPsVwy4z15gMYAx5iggA6iL03zvU1fnyo9VMaQNZdN+8mlD+RKp4f+9WcPe/R2Mz0/n08eNdnucmIjLW5DW2m5jzJXAKsAD3G2t3WSMuQl4zVq7AvgacJcx5quETsi/1Lq0CYnUT1So/9KGsmk/+aQ19M+6ze0REk6iNNzV2M79b4beyvvqSeV4U5NjC9O4nQNmrX0ceHzAbdf3+3ozsCBe8wylqEj3gpFOG8qm/eST1jBQfqnbIyScRGgY7LHc9vxeghY+clQJx4xOvE9mHq7kWEbGWCIddlWHRxvKpv3k04byJULDhzfV8W5dGyXZaVwu6LyvSOgCLIy8vLxD30klNG0om/aTT1pD79578O69x+0xEorbDXc0tPHn10Jbgn5lwQSyvR5X54m1RNqGImEEg0G3R1BR0oayaT/5pDXM3nQNoG9F9udmwzeqWrjxqZ10Bi2nTSlkbnm+a7M4RY+AheH3+90eQUVJG8qm/eTThvK51XD1jia+/cQO2rp6OHlyAdcsKndlDqfpEbAwRo9Ojo+4jmTaUDbtJ582lM+Nhg++VcvvX94HwLlHl/KFueNISZBPY8aaHgELY6ida5UM2lA27SefNpQvng17rOV3L1X2Lb6u+NBYvnji+KRdfIEeAQsrLS3N7RFUlLShbNpPPm0oX7waBoI9/Oy5vaze0URqiuFri8pZfKT7W2A4TRdgYeTnJ9/JfiONNpRN+8mnDeWLR0N/IMiNT+/kzapWstJSuH7JZD44TtYnaA+XvgUZRn19vdsjqChpQ9m0n3zaUD6nGzb4u/jao1t5s6qVosxUfvaRqSNm8QV6BCws/ZubfNpQNu0nn7SGehHugznZcO/+Dr79xA5qWgOMz0/nh2dOYUxuumM/LxHpAiyMQCDg9ggqStpQNu0nnzaUz6mGm2pauf7JnbR0BjlqVBY3nTGF/IyRtxwZec84Au3t7W6PoKKkDWXTfvJpQ/mcaPhqRTM3Pr2TQNAyd0Ie3148mYwkubj2cOkCLAzdv0Y+bSib9pNPWsPcF08FoGXBapcnSRyxbvhOrZ+behdfZ00v5ssLJuBJSd5tJg5lZC47D0H3r5FPG8qm/eST1jC1eQOpzRvcHiOhxLJhpa+D767aQWfQsnRaEVcvHNmLL9AFWFher9ftEVSUtKFs2k8+bShfrBo2tnXxrSd20NwZ5ITxeXxlYTkmiTdYjZQuwMLIzc11ewQVJW0om/aTTxvKF4uGbYEg31m1g+qWANNLs/jO4kmkjvAjXwfoAiyMhoYGt0dQUdKGsmk/+bShfNE27Ar2cNN/drG9oZ2xeencdMYRZKZ5YjSdfLoAC6OwsNDtEVSUtKFs2k8+bShfNA2ttdz2/F7W72shPyOVH505hcJMvTxVf7oAC0M/Pi2fNpRN+8mnDeWLpuHdr1bxn+1NZKSm8MOlUxibN7I2WY2EbkMRRkdHh9sjqChpQ9m0n3zSGnZOuNjtERLO4TZ8ZFMdD2ysJcXAdxdPZlppVownSw66AAtD2v416mDaUDbtJ5+0hm1H/9ztERLO4TR8blcTv1lXCcA1J5VzwoSRc23H4dK3IMOQtn+NOpg2lE37yacN5Rtuw43vtXLLmj1Y4LNzxnDGtGJnBksSugALIyMjw+0RVJS0oWzaTz5pDT2+N/H43nR7jIQynIa7Gtv53lM76Qpalh1VwgWzyxycLDnoW5BhZGZmuj2CipI2lE37ySetYd7a0wBoOqvR5UkSR6QNa1sDfPuJHfgDQRZMzOdL88brRqsR0CNgYTQ1Nbk9goqSNpRN+8mnDeWLpGFta4DrVm6nvq2LWWXZXHfqpBF/iaFI6RGwMIqL9X1r6bShbNpPPm0o36Ea7m3q4LontlPv72JyYQY3nn4E6al6XCdS+m8qjJaWFrdHUFHShrJpP/m0oXxDNXyn1s9XH91KvT905OunH5lKXoYe0xkO/bcVRiAQcHsEFSVtKJv2k08byjdYw9cqm7nx6V10dvcwd0Ie3148mQw98jVsugALQ9r+Nepg2lA27SefNpQvXMPVOxr5ybN76e6xLJlaxDUnlevFtQ+TLlnD0P1r5NOGsmk/+bShfAMb/mtTHTev3kN3j+Xjx4zi2kW6+IqGHgELQ9rHp9XBtKFs2k8+aQ2b5z/j9ggJ50BDay33ra/mL2+EFmSfO2Esn9R9vqKmC7AwvF6v2yOoKGlD2bSffNIaBvM/4PYICcfr9RLssfx6XSWPvlNPioGrF5Zz5nT9hGss6FuQYfh8PrdHUFHShrJpP/m0oXz1Tfv58erdPPpOPWkew/VLJuviK4b0CFgYJSUlbo+goqQNZdN+8klrmPX21YBelPuAtkCQ373dyYaaNrLSUrjpjCM4dkyu22MlFT0CFob+zU0+bSib9pNPWsP0intJr7jX7TESgj8Q5Bsrt7Ohpo3CzFR+9pGpuvhygB4BC6Orq8vtEVSUtKFs2k8+bShTINjDDU/t5N26NkoyU/jpsmmMzUt3e6ykpAuwMHT/Gvm0oWzaTz5tKE+Ptdy6Zg8b3mulKCuVm5dO0sWXg/QtyDB0/xr5tKFs2k8+bSiLtZbfrNvHc7v2k5WWwg+XTsG2Nro9VlLTI2BhZGdnuz2CipI2lE37yacNZXlgYw3/2lxHWorhhtOPYEpxFvVWGzpJj4CF4fF43B5BRUkbyqb95NOGcjy5tYG7X30PA3zjlIl8YGzohHtt6CxdgIXR3Nzs9ggqStpQNu0nn7SG3Xmz6c6b7fYYcfdKhY/bnt8LwBfnjWfREYV935PWUJqI3oI0xlwIvGmtfccYMx24CwgCX7TWbnFyQDeUlpa6PYKKkjaUTfvJJ61hy4LVbo8Qd+/U+vn+f3bTY+GC2WV8bNb7m0lrKE2kR8B+ABw4G++nwCvAs8CdTgzltsZGPfFQOm0om/aTTxsmtor9HXx31Q46u3s4Y2oRn50z5qD7aENnRXoSfqm1tsYYkwEsBD4OdAH1jk3mImut2yOoKGlD2bSffNowcTX4u/jWEzto7gzyoQl5XH1SOcaYg+6nDZ0V6QKszhhzJHAM8Kq1ttMYkwUcXCwJ6GFX+bShbNpPPmkNC1cWAdB0VnIf9fEHgnx71Q5qWgPMKM3i26dNIjUl/B/l0hpKE+lbkN8HXgf+CPyk97YlwAYnhnJbTU2N2yOoKGlD2bSffNow8RzY5X5nYzvj89P5/tIpZKYN/klHbeisiI6AWWvvMcb8vffrtt6bXwIucGowN+Xk5Lg9goqSNpRN+8mnDRPLwF3uf3TmFPIzhl4CaENnRfopyBSgo9/XkKTnfymllFLJJNhjue35ve/b5X50rl5iyG2RvgXZTeik+3D/JJ3W1la3R1BR0oayaT/5tGFi6O6x3LxmN09tayQ9NYWbzgjtch8JbeisSE/Cnzzg12OA64B/x3acxFBWVub2CCpK2lA27SefNnRfINjDD5/Zzbo9vr4jX7NGR/62ojZ0VkRHwKy1ewb88xJwCfANZ8dzR11dndsjqChpQ9m0n3za0F0d3aET7tft8ZGb7uGWDx85rMUXaEOnRXMx7jwgKT+jGm4/FCWLNpRN+8knraF/1m1ujxAz7V1Brn9yJxveayU/I5Wbz5oS8duO/UlrKE2kJ+HfB/TfkS0LWAT8xYmh3FZUVOT2CCpK2lA27SeftIaB8kvdHiEm/IEg335iB5tr/RRlpXLrWVMpL8w4rMeS1lCaSE/C3w7s6PfPS8D/WGuvcmowN+lhV/m0oWzaTz5tGH/NHd18/fFtbK71MyonjZ+dPe2wF1+gDZ0W6T5gNzo9SCLJy8tzewQVJW0om/aTT1pD7957ALlHwprau7ju8e3saupgTK6XWz88lbJcb1SPKa2hNBGfA2aMuQy4EBgLVAH3A3fbJLxYVDAYdHsEFSVtKJv2k09aw+xN1wAyF2D1/gDfeHw7Fb5OJuSnc8uHj6QkO7rFF8hrKE1Eb0EaY24l9InHh4D/Bf4JXAvc4txo7vH7/W6PoKKkDWXTfvJpw/ioaQnwtUe3UeHr5IiiDH76kakxWXyBNnRapOeAXQosttb+xlr7uLX2t8AZwGcj/UHGmDONMe8aY7YbY64b5D6fNMZsNsZsMsb8LdLHjrXRo0e79aNVjGhD2bSffNrQeft8HXztsa281xJgWkkWt354KoWZaTF7fG3orEgXYC29/wy8rTmS32yM8QC/Bs4CZgIXGmNmDrjPVOCbwAJr7Szg6ghni7nq6mq3frSKEW0om/aTTxs667ldTSx/5F1qW7uYOSqbWz58JHmHuLbjcGlDZw1ayxhzRL9f/hx4yBhzM1AJTCD0VuTtEf6cDwHbrbU7ex/7fuCjwOZ+9/k88GtrbROAtbY20icRa2lpsfsbhHKHNpRN+8mnDZ0R6O7hty/v49F3QpdjXjAxn6+fMpHMNE/Mf5Y2dNZQy+XthPb+6r8T26kD7nMa8KsIfs44oKLfryuBuQPuMw3AGPMi4AFusNY+EcFjx1x+fr4bP1bFkDaUTfvJpw1jr2J/Bz98Zjc7G9tJSzF84cRxLDuqxLENU7WhswZdgFlrI317MlZSganAKcB44DljzDHW2v3971RbW8vll19OamoqwWCQ8847j+XLl1NdXU12djYej4fm5mZKS0tpbGzEWktpaSk1NTXk5IQuw9Da2kpZWRl1dXUYYygqKqKuro68vDyCwSBVVVVMnz6d6upq0tLSyM/Pp76+nvz8fAKBAO3t7YwePZrq6mq8Xi+5ubk0NDRQWFhIe3s7HR0dfd/PyMggMzOTpqYmiouLaWlpIRAI9H0/MzMTr9eLz+ejpKQEn89HV1dX3/dj9Zz8fn/fY46E51RVVUV2dnZSPadk7DTYc9q5cyeTJk1KqueUjJ2Gek4+n4/09HQxz6mw98+YPXv2JGQrROA0AAAgAElEQVSnDftTuGt9A51By5hcL5+bmcGUgi7a2toc+/+e3+9n7NixCdVJ4utpMCYeu0gYY+YROqK1tPfX3wSw1v64331+C7xsrf1T76//A1xnrX21/2OtW7fOzpgxw9F59+/fT0FBgaM/QzlLG8qm/eTThrHR3hXkznWVrNraCMApRxTwlYXlZHtj/5bjQNowNtavX//64sWL5wy8PbZn7A3uVWCqMWYysA+4APifAfd5hNA+Y38yxpQQektyZ5zme59AIODGj1UxpA1l037yacPo7W5q54f/2c2e/R14PYbl88Zz5vTiuF2jURs6Ky4LMGtttzHmSmAVofO77rbWbjLG3AS8Zq1d0fu9M4wxm4Eg8L/W2oZ4zDdQe3u7Gz9WxZA2lE37yacND5+1lie2NnLn2go6g5byggy+fdokJhdlxnUObeiseB0Bw1r7OPD4gNuu7/e1Ba7p/cdVuveJfNpQNu0nn7SGuS+GPmPWsmC1q3O0BYL84sUKVu9oAmDptCK+NG+8I59yPBRpDaUZ9gLMGPO+k/OttT2xGycxVFdXM3HiRLfHUFHQhrJpP/mkNUxt3hD3nxno7mF3UwfbG9rY0dAe+qexnc7uHjJSU/jyggksmVoU97kOkNZQmogWYMaYDxLaSPVY4MCl1Q2hbSrivyx3mNcbm8s4KPdoQ9m0n3za8P38gSA7GtrY3tDO9oZ2dtS3sXd/B8Ewn4ObXprF/548kfKCjIO/GUfa0FmRHgH7M/Bv4DKgzblxEkNubq7bI6goaUPZtJ982jDkveZOvv+fXWxvOPh8qhQD5QUZTCnOZEpxJkcWZzKlOIv8GO9of7i0obMirTwR+LaNx54VCaChoaFvrw8lkzaUTfvJpw2huaObb6/aQaWvkzSPYXLhfxdaR5ZkMakww5VzuyKlDZ0V6QLsYUIX317l4CwJo7Cw8NB3UglNG8qm/eQb6Q0DwR5uenoXlb5OJhdmcNuyaXHZuyuWRnpDp0W6AMsAHjbGvAC8b1tXa+3FMZ/KZe3t7eTl5bk9hoqCNpRN+8k3khtaa7n9+b1srG6lOCuN7y+dIm7xBSO7YTxEugDbzPsvnJ3UOjo63B5BRUkbyqb95JPWsHNC7I4l3Le+mv9sbyIjNYXvn3EEo3JknswuraE0ES3ArLU3Oj1IItG9T+TThrJpP/mkNWw7+ucxeZwntzbwlzeqSTHw7dMmcWRJVkwe1w3SGkoz6AW3jTGL+n192mD/xGfM+Brq4plKBm0om/aTbyQ2fKOqhduf3wvAl+aNZ255vssTRWckNoynoY6A3Qkc3fv1Hwe5jwWOiOlECSAjw929V1T0tKFs2k8+aQ09vjcBCOZ/4LB+/56mdm56ehdBC+cfXco5M0tjOZ4rpDWUZtAFmLX26H5fT47POIkhMzO+19tSsacNZdN+8klrmLc29IZO01mNw/69TW1dfGfVTvyBIAsm5vP5ueNiPZ4rpDWUZtC3IEeypqYmt0dQUdKGsmk/+UZKw47uHq5/aic1rQGml2bxjVMnkWKM22PFxEhp6BZdgIVRXFzs9ggqStpQNu0n30hoGOyx3Lx6N+/WtVGW4+Wm048gIzV5/lgdCQ3dlDz/T4mhlpYWt0dQUdKGsmk/+UZCw7te2cfaPT5yvB5+uHQKhVlpbo8UUyOhoZt0ARZGIBBwewQVJW0om/aTL9kb/mtTHQ+9XUdqiuH6JZMpL0y+E9aTvaHbIlqAGWPeGOT212I7TmLQvU/k04ayaT/5krnh6h1N/OalSgCuXjiBD4xNzotWJ3PDRBDpEbAjB95gjDEk4RYUoHufJANtKJv2ky8ZG1preWBDDT9evZseC585bjRnTEve86SSsWEiGXInfGPMvb1fevt9fcAkYJMTQ7lNP3ornzaUTfvJJ61h8/xnhvx+sMdyx9oKHt/SAMDnPjSWTxwzKh6juUZaQ2kOdSmiHYN8bYEXgX/EfKIE4PXKvG6X+i9tKJv2k09aw6E2YG0LBPnBM7t4rbIFr8fw9VMmsmhyYRync4e0htIMuQA7cA1IY8xL1tpV8RnJfT6fj4KCArfHUFHQhrJpP/mSpWGdP8B3V+1kZ2M7+Rmp3Hj6Ecwsy3Z7rLhIloaJKtKLca8yxpwOXACMstYuM8bMAfKstUMftxWopKTE7RFUlLShbNpPPmkNs96+Gnj/Rbl3NLTxnVU7aWjrYnx+Oj9YOoWxeelujRh30hpKE+mnIK8CfgNsAw5cpLsd+IFDc7nK5/O5PYKKkjaUTfvJJ61hesW9pFf891TnVyuauebRbTS0dXF0WTY/XzZtRC2+QF5DaSI6AgZcDSy21u42xnyj97YtwHRnxnJXV1eX2yOoKGlD2bSffJIbPralnjterKDHwqlTCvnaSeV4k2iH+0hJbihBpAuwXKCi92vb+79pQFLu0qZ7n8inDWXTfvJJbfjHV/bxwMZaAC78QBmXHD8maa7tOFxSG0oR6ZL+OeC6Abd9GVgd23ESg+59Ip82lE37ySe14QMba0kx8NWTyvnsnLEjdvEFchtKEekRsKuAfxtjPg/kGmPeBVqAjzg2mYuys0fGJ1ySmTaUTfvJJ6lha2c3BzaVyEpL4buLJ3P8+DxXZ0oEkhpKFOmnIN8zxpwAfAgoJ/R25CvW2h4nh3OLx+NxewQVJW0om/aTT0rDprYuvrVqB/f3Xsrx9mXTmFykG5CCnIZSRXxWoQ152Vr7D2vtS8m6+AJobm52ewQVJW0om/aTT0LDmpYA1zy6jR0N7ewITqMj+1hdfPUjoaFkER0BM8ZU8N+T7/vrBCqBh4DfWGu7Yziba0pLS90eQUVJG8qm/eRL9IZ7mzq47ont1Pu7mFKciTntedqz0tweK6EkekPpIj0C9kugCbgR+BxwE9AA/Al4gNAJ+T9yYkA3NDY2uj2CipI2lE37yZfIDbfWtXHNo1up94f2+Prp2VMp1MXXQRK5YTKI9CT8S4HTrbVVB24wxqwEnrTWzjLGrAaeBr4e+xHjz9pwB/uUJNpQNu0nX6I2fLOqhe89tZP2rh5OGJ/Hd5dMJmME7vEViURtmCwiXYCNAVoH3OYHxvZ+vRVImgtG6WFX+bShbNpPvkRsuHbPfn74zG66gpZTpxRy7aJy0jyhxVfhyiIAms7Soz4HJGLDZBLpsv/fwL+MMUuMMTOMMUuAf/beDjAP2O3AfK6oqalxewQVJW0om/aTL9EaPrWtgZue3kVX0PKRo0r4xikT+xZfKrxEa5hsIj0C9gXgBuB3hI56VQH/IHQuGMBO4OxYD+eWnJwct0dQUdKGsmk/+RKp4cNv1/Kbl/YBod3tLz1+DGYEb7AaqURqmIwi3Qesg9BO+AN3wz/wfd0uVymlVEKx1nLf+mr+8kboj6grPjSWjx9b5vJUSoXo8dcwWlsHnu6mpNGGsmk/+dxuWOcP8PMXKvjLG9WkGLjmpHJdfA2T2w2TXaRvQY4oZWX6IpVOG8qm/eRzo6G1lreq/fxrcx0v7t5Pj4W0FMM3T53EwslJ8zmxuNHXobN0ARZGXV0dEyZMcHsMFQVtKJv2ky+eDdu7gjyzo4kVm+rY1dQBQIqBkycX8InZZUwryYrLHMlGX4fOGnQBZox5wFr7qd6vP2ut/VP8xnKXnpwpnzaUTfvJF4+G+3yd/PudOlZtbcQfCAJQmJnK2TNK+PCMYkqyvRE/ln/WbU6NKZa+Dp011BGwpcYYY0M7sf2C0K73I0JRUZHbI6goaUPZtJ98TjXssZbXKptZsbmeVyua+66RN3NUNufMLOGkyQWHtb1EoPzSmM6ZDPR16KyhFmDPA+uMMVuBDGPMveHuZK292JHJXFRXV8fEiRPdHkNFQRvKpv3kc6Lhuj0+fv/yPvY1dwKQ5jGcNqWQZTNL9W1GB+jr0FlDLcA+AXwcmEjoQtw74jJRAsjLy3N7BBUlbSib9pMvlg1rWgLcua6SdXt9AIzKSWPZUaWcOb2Y/IzYnMrs3XsPoEfC+tPXobMG/X9u795ffwEwxqRZa2+M21QuCwaDbo+goqQNZdN+8sWiYVewh4feruMvb1TT2d1DVloKlxw/hnNmluJJie35SdmbrgF0Adafvg6dFelGrDcYY6YCFwLjgH3A/1lrtzk5nFv8fj8lJSVuj6GioA1l037yRdtw43st3PFiJXv2hz7VePIRBfy/ueMpzk6L1YjqEPR16KyIFmDGmGXAX4FHgT3AdOA1Y8xF1toVDs7nitGjR7s9goqSNpRN+8l3uA2b2ru465Uqnt4Wuij22Lx0rpo/nuPH69th8aavQ2dF+ub5j4CPWmtXH7jBGHMK8Csg6RZg1dXVeuKhcNpQNu0n33Ab9ljL41sa+NNrVbR0BknzGC6YXcanji3Dm6oXbXGDvg6dFekCbDyhT0X290Lv7UknLU0PcUunDWXTfvINp+GOhjZ+8UIFW+raADh+XC5Xzh/PuPwMp8ZTEdDXobMiXYC9CXwNuKXfbdf03p508vPz3R5BRUkbyqb95IukYSDYw32vv8c/3qqlx0JRVipfPHE8iyYX6CagCUBfh86KdAH2ReDfxpivABXABKANWObUYG6qr68nOzvb7TFUFLShbNpPvkM13NXYzi1r9rCzsZ0UA+fOKuXi48eQ7fXEcUo1FH0dOivST0FuMcYcBZwIjAWqgJettV1ODucWXfXLpw1l037yDdawx1oefruOu1+roitoGZPr5eunTGRWWU6cJ3y/prMaXf35iUhfh86KeAc7a203ofO+kl4gEHB7BBUlbSib9pMvXMPa1gA/eXYPG95rBeCs6cX8vxPHkZmmR70Skb4OnRWbLYSTTHt7u9sjqChpQ9m0n3z9G1preWZHE79aW4k/EKQgI5WvnlTOvIl6hCWR6evQWboAC0P3PpFPG8qm/eQ70LC5o5s7Xqzg2V37AZhXns/VJ02gMDOxPmGX++KpALQsWH2Ie44c+jp0lm6uEkZ1dbXbI6goaUPZtJ981dXVvF7ZzBce2sKzu/aTmZbCNSeVc8PpkxNu8QWQ2ryB1OYNbo+RUPR16KxId8K/HfiztTYpt50YyOv1uj2CipI2lE37ydbR3cPft3Xw1O4dAMwclc03TpnImLx0lydTw6GvQ2dF+hakB1hljKkD7gP+aq2tdG4sd+Xm5ro9goqSNpRN+8nU3NHN4+/Ws2JTPfVtXXgMXHz8GD55bFnML56tnKevQ2dFug3Fl40xXwXOAj4NfMcY8zJwL/CQtbbVwRnjrqGhgZwcdz8SraKjDWXTfrLsaWrnkU11PL2tkc6gBWB8jodvLjmSqSVZLk+nDpe+Dp01nG0ogoQuxv2oMWYW8DfgHuBOY8z9wPestfsG+/3GmDOBXxA6mvYHa+3Ng9zvfOBB4ARr7WuRzhdLhYWFbvxYFUPaUDbtl/h6rOW1ymYe2VTHa5UtfbefMD6Pc48uZWquJT9fF1+S6evQWREvwIwxecAngM8AxwL/BL4E7CV0maKVvbeH+70e4NfA6UAl8KoxZoW1dvOA++UCXwFeHvYziaH29nby8vLcHEFFSRvKpv0SV3tXkKe3NfLwpjoqfZ0ApKemcPrUIj42q5TygtD1G2tqanQjT+H0deisSE/CfxBYCjwH/BZ4xFrb2e/71wC+IR7iQ8B2a+3O3vvfD3wU2Dzgft8ndL3J/430CTiho6PDzR+vYkAbyqb9Ek9ta4AVm+t4fEsDrYEgACXZaXxsZilnTi8mL+P9f5xIa9g54WK3R0g40hpKE+kRsJeAK621YT+Taq3tMcaUDfH7xxG6huQBlcDc/ncwxnwQmGCtfcwY4+oCTPc+kU8byqb9Eke9P8Df3qzhiXcb6O4Jnd81c1Q25x5dyoJJBaQOcnK9tIZtR//c7RESjrSG0kS6ADvJWvvTgTcaYx6y1p4HYK1tO9whjDEpwG3ApYe6b21tLZdffjmpqakEg0HOO+88li9fTnV1NdnZ2Xg8HpqbmyktLaWxsRFrLaWlpdTU1PSdTNja2kpZWRl1dXUYYygqKqKuro68vDyCwSBVVVVMnz6d6upq0tLSyM/Pp76+nvz8fAKBAO3t7YwePZrq6mq8Xi+5ubk0NDRQWFhIe3s7HR0dfd/PyMggMzOTpqYmiouLaWlpIRAI9H0/MzMTr9eLz+ejpKQEn89HV1dX3/dj9Zz8fn/fY46E51RVVUV2dnZSPadk7DTYc9q5cyeTJk1KquckrZNNz+Zvr1exprKTrh6LARZNymPhKJhemkVubhr7KvYO+px8Ph/p6ekJ9ZySsZOTz8nv9zN27Nikek5udBp07WOtjWSB1GytPeiNYGNMo7W2KILfPw+4wVq7tPfX3wSw1v6499f5wA7gwKcpRwONwDkDT8Rft26dnTFjxiFnjkZNTQ1lZUMd0FOJThvKpv3c09LZzYMba3l4Ux0d3T0ALJpcwMUfHEN5YUbEjyOtoccX2uYymP8BlydJHNIaJqr169e/vnjx4jkDbx/yCJgx5qbeL739vj7gCGBPhD//VWCqMWYysA+4APifA9+01vqAkn4/dw1wrVufgszMzHTjx6oY0oayab/4awsEeXhTHQ++VYu/9xyvuRPyuHTOGKYUD//TjNIa5q09DYCmsxpdniRxSGsozaHegpzQ+78p/b4GsITO6bohkh9ire02xlwJrCK0DcXd1tpNvYu616y1K4Y1tcOampr0kx/CaUPZtF/8dHb38O/NdTywsRZfRzcAx43N5dI5YzhqVPZhP642lE8bOmvIBZi19rMAxpi11tq7ovlB1trHgccH3Hb9IPc9JZqfFa3i4mI3f7yKAW0om/ZzXneP5fEt9fztzWoa20ILr5mjsvnsnDHMHhv9DujaUD5t6KxBF2DGmEnW2t29v/yPMeaIcPc7sLVEMmlpadHdf4XThrJpP2e9WdXCr9dWsmd/aJuBI4szuXTOGE4Yn4cxsblkkDaUTxs6a6gjYG8BB/4atJ3Q244DX5mW0FuKSSUQCLg9goqSNpRN+zmjzh/g9y/t49ld+wEYm5fO504Yy4JJ+TFbeB2gDeXThs4adAFmrc3t93VKfMZJDLr3iXzaUDbtF1tdwR4eeruOv75RTUd3D+kew4UfGM3Hjx2F1+PMf961oXza0FkjamEVqaH27VAyaEPZtF/svF7ZzBce2sIfX62io7uHhZMK+MPHZ/I/x412bPEF2jAZaENnDXUO2POE3mIckrV2UUwnSgD60Vv5tKFs2i96ta0BfvvSPl7YHXq7cXx+Ol+aN5454+PzqTZpDZvnP+P2CAlHWkNphjoH7A9xmyLBeL1et0dQUdKGsmm/wxcI9vDgxlr+781qOoOW9NQUPn1cGecd7dzbjeFIa6gbsB5MWkNphjoH7M/xHCSR+Hw+CgoK3B5DRUEbyqb9Ds9rlc38am0lVc2dAJw8uYDPzx3HqJz4/0GqDeXThs4a6i3Ii6y19/V+fdlg97PW3u3EYG4qKSk59J1UQtOGsmm/4Wlo6+K3L1Xy7M7Q243lBRksnzee48ZFv5/X4ZLWMOvtqwG9KHd/0hpKM9RbkBcC9/V+fdEg97FA0i3AfD4f2dmHvwO0cp82lE37RabHWh57p54/vlpFW1fo042f+eAYzju6lLQ4vt0YjrSG6RX3AroA609aQ2mGegvyw/2+PjU+4ySGrq4ut0dQUdKGsmm/Q9vR0MYvXqhgS10bAB+akMeV88czOjfd5clCtKF82tBZh7oWZB9jTAFwNjAWqAIes9bud2owN+neJ/JpQ9m03+Dau4Lct76ah96upcdCcVYaX5w3jpMmFcR8M9VoaEP5tKGzIjpGbYw5DdgNfBk4AbgK2G2MWezcaO7RvU/k04ayab/w1u3x8fl/vsODb9ViLXx0Zil/+PhRLJpcmFCLL9CGyUAbOivSI2C/Aq6w1v79wA3GmE8AvwZmODGYm/Q9b/m0oWza7/3q/AHuXFvJi3t8QOjajVcvLGdaaZbLkw1OG8qnDZ0V6QJsLPDPAbc9DNwV23ESg8eTdJe3HHG0oWzaL6THWlZsrudPr1XR3tVDZloKlx4/hnNmluJJSawjXgNpQ/m0obMi/ZjMfcDyAbd9Ebg3tuMkhubmZrdHUFHShrJpv9BRr2+u3M6d6ypp7+ph4aR8/vDxozj36FEJv/gCeQ2782bTnTfb7TESirSG0kR6KaIU4P8ZY74O7APGAWXAS45P6ILS0lK3R1BR0oayjfR+a3Y08csXK2gNBMnPSOUrCyewcJKsDTGlNWxZsNrtERKOtIbSDOdSREn5dmM4jY2NZGUl7rkV6tC0oWwjtV9rZze/WlvJMzuaAJg7IY9rTiqnMCvN5cmGb6Q2TCba0Fl6KaIwrD3kNchVgtOGso3EfhuqWrj12T3U+btIT03hC3PHcfaM4oT7dGOkRmLDZKMNnTWcfcDKgA8BJUDffxGS8VJEethVPm0o20jqFwj2cM9r7/HPt2qxwPTSLL5xykTG52e4PVpUpDUsXFkEQNNZjS5PkjikNZQm0n3APgbsAG4CfkdoH7DfMfglikSrqalxewQVJW0o20jpt6uxnS//610efKsWY+Azx43m9mXTxC++YOQ0TGba0FmRHgH7AfBZa+0/jDFN1trjjDGfBWY5OJtrcnJy3B5BRUkbypbs/Xqs5aG36/jTq1V09VjG5qXzjVMmctSo5Nl3KdkbjgTa0FmRLsDKrbX/GHDbn4Fq4NrYjqSUUsnJWsub77Xy1/XVbKxuBeDDM4r5wtxxZKbpnktKjSSRLsBqjTFl1toaQpcgmgfUA0n5X4zW1laKi4vdHkNFQRvKlmz92gJBnt7eyIrN9ezd3wFAQUYqXz2pnHkT812ezhnJ1nAk0obOinQBdhewkNBu+LcDq4Ee4GcOzeWqsrIyt0dQUdKGsiVLv4r9HazYXM9T2xpo6+oBQhfPPvuoEpYdVUJ+RsSfgxInWRqOZNrQWRG9+q21t/T7+l5jzBog21r7jlODuamuro4JEya4PYaKgjaUTXK/YI/llYpm/rW5jvX7WvpuP2Z0Dh+dWcL8SQWkCtjJPlqSG6oQbeis4WxD4QFOJHRdyCqSdBd8QOy+O+q/tKFsEvs1d3TzxLsN/PudempaAwCkewynHVnEOTNLmFI8sja0lNbQP+s2t0dIONIaShPRAswYcyzwCJABVALjgQ5jzLnW2g0OzueKoqIit0dQUdKGsknpF+yxbHivhVVbG3lx934CwdDGlWNyvSybWcrSaUXkpifv24xDkdLwgED5pW6PkHCkNZQm0v8y3A38GrjNWmtNaFn81d7bj3dqOLfU1dUxceJEt8dQUdCGsiV6v32+Dp7c1sjT2xqp83f13T5nfC4fnVnKCRPySBnhRw8SvaE6NG3orEgXYNOAn9ve6xL0LsJ+Adzg1GBuysvLc3sEFSVtKFsi9vMHgjy3az9PbW3g7Rp/3+1lOV7OmFbEkqlFjMlNd3HCxJKIDYfi3XsPoEfC+pPWUJpIF2CPA+cAD/e7bRnwWMwnSgDBYNDtEVSUtKFsidKvx1o2vNfKU1sbeH63j87u0CcZ01NTWDS5gDOmFnHMmJwRf7QrnERpGKnsTdcAugDrT1pDaQZdgBlj7gMOXInTA9xvjHkdqAAmEHrr8V+OT+gCv99PSUmJ22OoKGhD2dzuV+8P8MS7Daza2th3Qj3AsaNzOGNaEQsnFZDlTcptEGPG7YYqetrQWUMdAds+4Ndv9/t6M7Aq9uMkhtGjR7s9goqSNpTNjX491rJ+XwuPvVPPur0+enr/+lmW4+X0qaG3GMfm6VuMkdLXoHza0FmDLsCstTfGc5BEUl1drSceCqcNZYtnv/3tXTy5tZHHttTzXkvoaJfHwKLJBZw9o4TZY/UtxsOhr0H5tKGzhrMP2CnAxcA4YB9wn7V2tUNzuSotLc3tEVSUtKFsTvez1vJWtZ/HttTzwq79dPUe7irL8fLhGcUsnVZMUZb+fyga+hqUTxs6K9J9wD4H/Aj4A/AyUA78nzHmu9bauxyczxX5+cl5bbaRRBvK5lS/1s5untrWyGNbGvquyZhi4MTyPD5yVAnHj8vDMwJ2qY8HfQ3Kpw2dFekRsK8Dp/ffdNUY8wCha0Mm3QKsvr6e7Oxst8dQUdCGssW6X70/wD/fquWxLQ109H6SsSgrlbOml3DW9GJG5Xhj9rNUiL4G5dOGzop0AVZM6MT7/t4FknKbXF31y6cNZYtVv4r9HfxjYy1Pb2+ku/dtxg+MzWHZUaXMm5g/Iq7J6BZpr8GmsxrdHiHhSGsoTaQLsBeA24wx37DWthljsoEfA2udG809gUDg0HdSCU0byhZtv611bdy/oYYXd+/HEnqb8eTJBXxqdhlHloysazK6RV+D8mlDZ0W6APt/wP2AzxjTSOjI11rgQqcGc1N7e7vbI6goaUPZDqeftZY3q1q5f0MNb1S1AJCWYjh9WhGfOGYU4/IzYj2mGoK+BuXThs465AKs97qPmcBiYDQwFqiy1lY6PJtrdO8T+bShbMPp12Mta3f7eGBjDe/WtQGQmZbCR2aUcN4xoyjWTzO6QtprMPfFUwFoWZCUH+4/LNIaSnPIBVjvdR/fAnJ7F11Ju/A6QPc+kU8byhZJv7ZAkKe3N/LIpjoqfZ0A5Gekct7RpXzkqBJy0yPeZUc5QNprMLV5w6HvNMJIayhNpP+FeoPQBbm3ODhLwvB69RNR0mlD2YbqV7G/gxWb63lqWwNtXaFPNJblePnEsaNYOq2Y9NSUeI2phqCvQfm0obMiXYCtAZ4wxtxD6FqQB64RibX27tiP5a7c3Fy3R1BR0oayDewX7LG8UtHMis11vL6vpe/2Y0bn8NGZJcyfVKCfaEww+hqUTxs6K9IF2AJgF3DygNstkHQLsIaGBnJyctweQ0VBG8p2oF9zRzertjbw73fqqe69TFC6x3DakUWcM7OEKcX6icZEpfTOrfEAACAASURBVK9B+bShsyJagFlrT3V6kERSWFjo9ggqStpQNp/J4vbn9/LM9kY6g6ED7qNzvZxzVAlnTCsmL0PP70p0+hqUTxs6a8j/ihljsoDvAEcD64EfW2s74zGYm9rb28nLy3N7DBUFbShPTUuAZ3c1sWZHE9sb/vvx9znjczlnZiknjNfLBEmir0H5tKGzDvXXyF8Dc4CVwMcJ7Yh/ldNDua2jo8PtEVSUtKEMDW1dPLeziWd37mdzrb/v9sxUw9LpJSw7qoQJBbp/l0TSXoOdEy52e4SEI62hNIdagJ0JfNBa+54x5g7gOUbAAkz3PpFPGyYuX0c3L+zez5odTWx8r7XvEz3pHsOJE/M55YhCjh2VTm5WpqtzquhIew22Hf1zt0dIONIaSnOoBVi2tfY9AGtthTFmRFwYSvc+kU8bJhZ/IMiLu/ezZmcTb+xrofe0LtJSDHMm5HHKEYWcWJ5HZpoHgD179pCr/UTT16B82tBZh1qApRpjTgXMIL/GWvuMU8O5JSND3/KQThu6r8da3qxq4cmtjby4e3/fyfQpJnRe1ylHFDJ/Yj45YTZM1X7ySWvo8b0JQDD/Ay5PkjikNZTmUAuwWt6/zUTDgF9b4IhYD+W2zEx960M6beiequZOntzawFPbGqnzd/XdfszoHE6dUsjCSfkUZA59eSDtJ5+0hnlrTwOg6axGlydJHNIaSjPkAsxaOylOcySUpqYm/eSHcNowvtoCQZ7fvZ9VWxt4u/q/J9OX5Xg5fWoRp08rYkxuesSPp/3k04byaUNn6WY6YRQXF7s9goqSNnRej7W89V4rT25r5Pld++noDl0WKD01hZMmF3DG1CKOHZNDihn+1hHaTz5tKJ82dJYuwMJoaWnR3X+F04bOaWrv4smtjTy+pZ73enenBzh6dDZnTC1m0eQCsryeqH6G9pNPG8qnDZ2lC7AwAoHAoe+kEpo2jC1rLRvfa+XRLfW8uNtHd0/ohPrS7LTQW4xTixmXH/lbjIei/eTThvJpQ2fpAiwM3ftEPm0YG80d3Ty1rZHHttRT6QtdBCPFwLzyfM4+qpjjxzmzO732k08byqcNnaULsDB07xP5tOHhs9byTm0bj26p57mdTQR6t48ozkrjrOnFnDm9mFE5Xkdn0H7yaUP5tKGz4rYAM8acCfwC8AB/sNbePOD71wCfA7qBOuAya+2eeM3Xn370Vj5tOHztXUGe7j3atbPxv5cgmTM+l7NnlHBieX7crsWo/eST1rB5ftJtaRk1aQ2licsCzBjjIXRdydOBSuBVY8wKa+3mfnd7A5hjrW0zxnwRuBX4VDzmG8jrdfZv98p52jBy1S2drNhcz8p3G/AHggDkZ6Ry5rQizppRwti82J3bFSntJ5+0hroB68GkNZQmXkfAPgRst9buBDDG3A98FOhbgFlrV/e7/0vAZ+I020F8Ph8FBQVu/XgVA9pwaNZa3q7x8/Dbtazd46P3nHpmlWVzzsxSFkzKx+tJcW0+7SefNpRPGzorXguwcUBFv19XAnOHuP/lwEpHJxpCSUmJWz9axYg2DC8Q7OHZnU08/HYd2xvaAUhNMZw6pYBzZ41iWmmWyxOGaD/5pDXMevtqQC/K3Z+0htIk3En4xpjPAHOAk8N9v7a2lssvv5zU1FSCwSDnnXcey5cvp7q6muzsbDweD83NzZSWltLY2Ii1ltLSUmpqavr2M2ltbaWsrIy6ujqMMRQVFVFXV0deXh7BYJDq6mqmTp1KdXU1aWlp5OfnU19fT35+PoFAgPb2dkaPHk11dTVer5fc3FwaGhooLCykvb2djo6Ovu9nZGSQmZlJU1MTxcXFtLS0EAgE+r6fmZmJ1+vF5/NRUlKCz+ejq6ur7/uxek5+v7/vMUfCc6qtrSU9PT2pnlM0nTzZ+fx7cx3/2e2nORA63JXrTeGMKXmcWApHjivF59vPnj11CfGcdu/eTXl5+YjrlEzPye/34/F4xDynD1TcC8A7uV8dUZ2Gek6dnZ2MGjUqqZ6TG50GXe9Ya4e7Rho2Y8w84AZr7dLeX38TwFr74wH3WwLcAZxsra0N91jr1q2zM2bMcHTePXv26Cc/hNOGIbsa2/nnW7Ws3tFEV+/7jEcUZXDu0aM49YhCvKnuvc04FO0nn7SGhSuLAL0WZH/SGiaq9evXv7548eI5A2+P1xGwV4GpxpjJwD7gAuB/+t/BGHMc8DvgzMEWX/Gie5/IN5IbWmt5q7qVBzbU8mplMwCG0N5d5x5dyuwxOZjDuDxQPI3kfslCG8qnDZ0VlwWYtbbbGHMlsIrQNhR3W2s3GWNuAl6z1q4AfgLkAP/o/cNhr7X2nHjMN5DufSLfSGwY7LGs2+PjgY01vFvXBkC6x3Dm9GLOPXqUK59mPFwjsV+y0YbyaUNnxe0cMGvt48DjA267vt/XS+I1y6FkZ2e7PYKK0khqGOju4entjTz4Vm3fbvV56R4+OquUc2aWkp+RcKd6HtJI6pestKF82tBZ8v7LHAceT3QXElbuGwkNWzu7eXRLPY+8XUdjezcAZTlePn7MKM6YVkRmmtx/ByOhX7LThvJpQ2fpAiyM5uZmCgsL3R5DRSGZG9b7Azz0dh2Pb6mnrasHgCOKMvnU7FEsmlwYt93qnZTM/UYKaQ2782a7PULCkdZQmv/f3p0Hx3Gedx7/PriIc4YYAAQPAKQoUjzE6CRFyjpikbZDJlKUKI7ilOPYKTlVW2VXOau4nKOSbDa1qbU3G8euWqeSiq84ceI4lu3IjmVbB63DoiiKlBUdFAWSEgiQAjDAgDM4ievdP2YAjqgmKRLo6XmB36eKpeme5vQz+rGJh91vv60GLEBTU1PUJcgcLbQMJ6cdz3am+eGRfp7tzMxOnHrdylruvaaZG1fVFf3A+kux0PJbjHzLcPCWvRffaJHxLUPfqAELkEqlqK4ujgkp5fIslAxPpsf44ZF+Hm5PzV5mLCsxblsT5/3XLGND08Ico7FQ8lvMlKH/lGG41IAFKMTcaBIunzMcm5zmqddP89CRfl7sHppd37a0kt1XJdi1PkF9VXmEFYbP5/wkSxn6TxmGSw1YAJ129Z+PGbb3jfDQkX72HhuYfSj2krIS3r12Kbs3NLB5Wc2Cusx4IT7mJ2/lW4aaiPXtfMvQN2rAAvT09GjuE8/5kuHoxBR7jw3w/cN9s89mBNjQVM2eDQ38/Np6aioW351IvuQn56cM/acMw6UGLMDMs57EX8WeYVd6jO+90seP21OzZ7vqlpTynnUJdm9o4IpEVcQVRqvY85OLU4b+U4bhUgMmUiBT045nTqT53uE+Dp0cnF2/eVkNd25q5PYrlhbtsxlFRGR+qQELMDQ0RENDQ9RlyBwUU4YDIxM8dKSf/3y1j+TwBJB9RNDOdQnu3NTI+kbdZXSuYspPLo8y9J8yDJcasADNzc1RlyBzFHWGzjle6RnmwcN9PPn6aSZzE3etjC3hrk2NvO+qBHVLdPidT9T5ydwpQ/8pw3DpJ0CAZDJJa2tr1GXIHESV4eS044njAzzwUi/tfdlB9SUGN6+Oc9emRm5YVUfJIrmTcS50DPpPGfpPGYZLDViAxXKr/0JW6AwHz0zyg1f7+Y+Xk/SNZC8zxivL2LOhgV/a2EhzXUVB6/GdjkH/+Zbh8NWfjbqEouNbhr5RAxYgkUhEXYLMUaEyPJke4zsvJ/nRaynOTGafy9i2tJJf29LEznUJlmhQ/WXRMeg/3zIcb/tI1CUUHd8y9I0asADJZFJzn3guzAydc7zYPcwDL/XyTEeambmib1xVxz1blrG1ZWE9lzEKOgb9pwz9pwzDpQYsQCwWi7oEmaMwMpyYmuaJ10/zwIu9s5OmlpcYO9fVc8+WZYt+7q75pGPQf75lWHHiq4DOhOXzLUPfqAELMDU1FXUJMkfzlWFmbJIDXRme6UhzoCvDyET2MmO8soy7NjVy16ZG6qsX9nMZo6Bj0H++ZVjz8v2AGrB8vmXoGzVgAYaHh2lsbIy6DJmDuWR4Mj3GvhPZpuulniGm855HuzZRyd2bNb4rbDoG/acM/acMw6UGLMDy5cujLkHm6FIynJp2HO4d5pkTafZ1pOlMn5l9r9Tg+pW17GiLs6MtzorYkjDKlXPoGPSfMvSfMgyXGrAA3d3dGnjouYtlOO0cL745xCNHUzxzIkN6bHL2vdqKUra1xtjRFmdbSx21mjC14HQM+k8Z+k8Zhks/WQKUl2tMj+/Ol2Hn6TEeOZri0aMpeocmZtevqKtgx+o4N7fF2bK8lrIS3cUYJR2D/lOG/lOG4VIDFiAej0ddgsxRfoaZsUl+cnyAR9pTvJocmV2/rLacXesS3HFlPauXVmrqiCKiY9B/ytB/yjBcasAC9PX1UVNTE3UZMgfdvUm6bYJH2lPs78zMPouxqryE269YynvWJfi5FbV6LFCR0jHoP2XoP2UYLjVgAdT1+2dq2tGVHqO9b5RXeoZ5/HiawfEBIPssxq0tdbxnXYJ3rVlKpe5eLHo6Bv3nW4YDe1JRl1B0fMvQN2rAAoyPj0ddglzA5LSjY2CUo/2jtPeNcLRvlGOp0dlHAc1YU1/Je9cn2HllgoYajWXwiY5B/ylD/ynDcKkBCzA6Ohp1CZLnVOYMz58a5GjfCEf7RzmeGmViyr1tu+baCtY3VrGuoZrWsmFu3bJW47o8pWPQf8rQf8owXGrAAmjuk2hNTE3zUvcw+zvTPNuZoStvXq4ZK2NLWN9QxbrG6tmmK1Z59o/zmTNn1Hx5TMeg/3zLsO6ndwAweMveiCspHr5l6Bs1YAE090nh9Q9P8Gyu4Tp0apDRibOXE2srSrlhVR0bm6pZ31jNlQ1VF52bSxn6Tfn5z7cMyzIvRF1C0fEtQ9+oAQtQUVERdQkL3tS040hyhP2daQ50ZmYfbj1jTX0l21tjbGuNc3VzDaWXOC+XMvSb8vOfMvSfMgyXGrAAdXV1UZewIDnneK1vhMeODfD4sQFSo2dnn19Saly3so7tbXFuao2xrHZuB74y9Jvy858y9J8yDJcasAD9/f3U1tZGXcaC0ZUe47GjA+w9NsDJzNnxXM21Fexoi3FTa5xrV9RSMY/TQyhDvyk//ylD/ynDcKkBC1BfXx91Cd7rH57gJ8cHeOxYiva+s5cX66vKePfaeu64sp4NTdWhDZRXhn5Tfv5Thv5ThuFSAxZgdHSUWCwWdRneGTozyVNvpHnsWIoXTg0xM1FEdXkJt65Zyh1X1nPdyrpLHs91OZSh35Sf/5Sh/5RhuNSABRgbG4u6BC845zieGuW5rkEOdGZ4uWeImem5ykuMm1pj3LGunu2tcZYUePZ5Zeg35ec/3zI80/rbUZdQdHzL0DdqwAJo7pPzGzozyaFT2Ybrua5B+kcmZt8rMbh2RS071yW4dU2cuotMFREmZeg35ec/3zIc2fK5qEsoOr5l6Bs1YAE098lZzjmO9Y9yoCvDgc4Mr/QOM503CX2iuoxtLTG2tcS4flVdpE1XPmXoN+XnP2XoP2UYruL4aVlkKisroy4hUpmxSQ6eHOS5rgwHuzJvmS6ixODnlteyrbWObS0x1iaqinLG+cWeoe+Un/98y7A0/TMApuLXRVxJ8fAtQ9+oAQtQVVUVdQkFNTXteDU5zMGuQQ50ZXgtOUL+kxYbq8vZ2hJjW2uMG1bVUVNRGlmt79Riy3ChUX7+8y3D2NM7ARjYk4q4kuLhW4a+UQMWYGBgYMHf+ZEcHue5rkEOdmU4dHKQofGp2ffKS4wty2vY2hJja0uMNfWVRXmW60IWQ4YLmfLznzL0nzIMlxqwAA0NDVGXMO+mnePV3hH2dZxmf2eGNwbeenfLqtgStrbUsbUlxjUraqkqL/6zXBeyEDNcTJSf/5Sh/5RhuNSABRgcHFwQs/+emZzm0MlB9nWkeeZEmtNjZ8dyVZaVcP3KOm7MNV0rY0sirHT+LZQMFyvl5z9l6D9lGC41YAHGx8ejLuGypccm2X8izdMdaQ6eHOTM5PTse821Fdy8Os7NbXG2LK+hvLSwc3MVks8ZivJbCJSh/5RhuNSABfBp7pOpaUdneoznOjM8fSLNKz1vnSZifWMVN69eyrva4lyR8G8s1+XyKUN5O+XnP2XoP2UYLjVgAYp57pP+4QkOJ4c5khzh1d5h2vtGGJk4e5arrMS4fmVt9kzX6jhNNRURVhudYs5QLk75+U8Z+k8ZhksNWIBiufV2ZHyK9r4RXk2OcCQ5zKu9I/TlzTw/Y1ltOVc313JzW5xtrTEvpokIW7FkKJdH+fnPtwwz73os6hKKjm8Z+kYNWICKiujOGqXHJnn0aIpH2lMcT42+5XIiQE1FKVc1VrNxWTUbm2rY0FRNoro8mmKLWJQZytwpP//5lqEmYH073zL0jRqwAOl0mqVLlxZsf9PO8fzJQX54pJ+nO9JM5LqushJjXUMVG5qyDdeGphpa4ksoWSTjuOai0BnK/FJ+/lOG/lOG4VIDFqCxsbEg++kdGufHr/Xzo9dS9Axl7zYpMdjeGuMXNjRwU0uMirKFe6dimAqVoYRD+fnPtwyrX/o9QA/lzudbhr5RAxYgnU5TU1MTymdPTE2z/0SGh47081xXZvaRP821Feze0MD7rkos2oHz8ynMDCV8ys9/vmW4pPNrgBqwfL5l6Bs1YAEmJt4+0H0upp3jSHKEJ18/zcPtKdK5CVHLS4xb1sTZs6GRa1fW6tLiPJrvDKWwlJ//lKH/lGG41IAFmI+5T4bHpzh0cpD9J9I825l5yyz0a+or2bOhgV3rEsQqFUEYNH+N35Sf/5Sh/5RhuPTTP8Dlzn1yMn2G/Z1p9p9I82L3MJN5tzA211awoy3GznUJNjZVL5oJUaOi+Wv8pvz8pwz9pwzDpQYswDu95j017Xixe4j9J9Ls78zQlT4z+16JwZblNexojXNTW4zVSxfPLPTFQOMW/Kb8/KcM/acMw6UGLEBp6cUnMn2tb4S/eryDjoGx2XW1FaVsa42xvTXG1paYLi9G6J1kKMVL+flPGfpPGYZLHUKATCZDfX194HsTU9P8y896+NefdTPtsrPQ//wV9Wxvi3N1cw2lJTrLVQwulKEUP+XnP98ynIxdG3UJRce3DH2jBixAU1NT4Ppj/SP83ydOcKx/FIBfvbqJ39m2kkrN1VV0zpeh+EH5+c+3DAdv2Rt1CUXHtwx9owYsQCqVorq6enZ5ctrxby/08PXnu5mcdiyvq+CTt7dxzYq6CKuUCzk3Q/GL8vOfMvSfMgyXGrAAzp29e/GNgVH+6vEO2vuyZ73u2tTIR29aSVW5ro0Xs/wMxT/Kz3/K0H/KMFwFa8DMbDfweaAU+KJz7tPnvL8E+BpwI9AP/IZz7o1C1ZevqamJqWnHt17s5WsH32Ri2rGstpz7b2vjhlWxKEqSS6RT535Tfv7zLcP6hxIADOxJRVxJ8fAtQ98UZPCSmZUCXwD2AJuB3zSzzedsdh8w4JxbB/wN8JlC1BbkheOnuP/7r/GlA6eYmHbs2dDA39+zSc2XR3p6eqIuQeZA+flPGfpPGYarUGfAbgKOOueOA5jZN4C7gVfytrkb+PPc628B/8/MzBXwHOi0c3znpSRfPpBhYtrRWF3Of7+tjW2tarx8U1tbG3UJMgfKz3/K0H/KMFyFasBWAZ15y13A9vNt45ybNLM00AD05W/U29vLfffdR1lZGVNTU9xzzz187GMfo7u7m5qaGkpLS8lkMjQ1NZFKpXDO0dTURE9Pz+wfpqGhIZqbm0kmk5gZiUSCZDJJ7+QS/n7/mwDsXBvn7jYjXjrE8HApfX19xONxxsfHGR0dZfny5XR3d1NRUUFdXR39/f3U19czOjrK2NjY7PuVlZVUVVUxMDBAQ0MDg4ODjI+Pz75fVVVFRUUF6XSaxsZG0uk0ExMTs+/P9TvFYjGmpqYYHh6e/czy8nLi8fiC/k7pdJqhoaEF9Z0WYk7n+07JZJLy8vIF9Z0WYk4X+k7j4+MMDQ15851mJlvo6OhYVDld6DtNTU1RWlq6oL5TFDmdjxXiBJOZvR/Y7Zz7aG75Q8B259zH87Z5KbdNV275WG6btzRg+/btcxs3bgyt1q8dfJP49BB3b1sf2j4kfB0dHXqEhseUn/98y1BjwN7OtwyL1aFDhw7u2rVr67nrC3UG7CTQmrfcklsXtE2XmZUBcbKD8Qvqt29cwdjY2MU3lKLW3NwcdQkyB8rPf8rQf8owXIWaQfQAsN7MrjCzCuADwIPnbPMg8OHc6/cDjxVy/Fe+ZDIZxW5lHilDvyk//ylD/ynDcBXkDFhuTNfHgR+RnYbiy865l83sL4DnnHMPAl8C/snMjgIpsk1aJPTQbP8pQ78pP//5luHw1Z+NuoSi41uGvinYPGDOuR8APzhn3Z/lvR4Dfr1Q9VxIIpGIugSZI2XoN+XnP98yHG/7SNQlFB3fMvSNHmIYQKdd/acM/ab8/KcM/acMw6UGLEAspnm/fKcM/ab8/OdbhhUnvkrFia9GXUZR8S1D3+hZkAGmpqaiLkHmSBn6Tfn5z7cMa16+H9ClyHy+ZegbnQELMDw8HHUJMkfK0G/Kz3/K0H/KMFxqwAIsX7486hJkjpSh35Sf/5Sh/5RhuNSABbjQowPED8rQb8rPf8rQf8owXGrAAnz3u9+NugSZI2XoN+XnP2XoP2UYLjVgAb797W9HXYLMkTL0m/LznzL0nzIMlxqwAJOTk1GXIHOkDP2m/PynDP2nDMNlET1u8bI9+uijSaAjzH2kUqnGRCLRF+Y+JFzK0G/Kz3/K0H/KcN6s3rVrV9O5K71rwERERER8p0uQIiIiIgWmBkxERESkwNSA5TGz3WZ2xMyOmtkfRl2PXJyZfdnMes3spbx1CTN72Mzac/+tj7JGuTAzazWzvWb2ipm9bGafyK1Xjp4ws0oze9bMXshl+D9z668ws/25v1P/zcwqoq5Vzs/MSs3seTP7fm5Z+YVIDViOmZUCXwD2AJuB3zSzzdFWJe/AV4Hd56z7Q+BR59x64NHcshSvSeD3nXObgR3Ax3LHnnL0xxlgp3PuWuA6YLeZ7QA+A/yNc24dMADcF2GNcnGfAA7nLSu/EKkBO+sm4Khz7rhzbhz4BnB3xDXJRTjnngBS56y+G/jH3Ot/BH6loEXJJXHOvemcO5R7PUj2B8AqlKM3XNZQbrE898sBO4Fv5dYrwyJmZi3ALwFfzC0byi9UasDOWgV05i135daJf5qdc2/mXncDzVEWI++cma0Brgf2oxy9krt89TOgF3gYOAacds7NTCalv1OL2+eATwHTueUGlF+o1IDJguay86xorhUPmFkt8ADwe865TP57yrH4OeemnHPXAS1kryhsjLgkeYfM7E6g1zl3MOpaFpOyqAsoIieB1rzlltw68U+Pma1wzr1pZivI/otcipiZlZNtvr7unJt5/oly9JBz7rSZ7QVuBpaaWVnuLIr+Ti1etwC/bGa/CFQCMeDzKL9Q6QzYWQeA9bm7PiqADwAPRlyTXJ4HgQ/nXn8Y+I8Ia5GLyI01+RJw2Dn32by3lKMnzKzJzJbmXlcB7yU7lm8v8P7cZsqwSDnn/sg51+KcW0P2Z99jzrkPovxCpZnw8+S6/88BpcCXnXN/GXFJchFm9q/Au4FGoAf4H8B3gW8CbWQfW3Wvc+7cgfpSJMzsVuBJ4EXOjj/5Y7LjwJSjB8zsGrKDtEvJ/sP+m865vzCztWRvaEoAzwO/5Zw7E12lcjFm9m7gk865O5VfuNSAiYiIiBSYLkGKiIiIFJgaMBEREZECUwMmIiIiUmBqwEREREQKTA2YiIiISIGpARMRb5nZH5vZFy/w/gfN7McFqqXNzIbMrLQQ+xMRv2kaChEpGDP7I+B259yevHXtwNGAdX/qnPvGJXz2GuB1oDzv+XWhMbM3gI865x4Je18isvDoDJiIFNITwLtmzhLlHjFUDlx/zrp1uW1FRBYkNWAiUkgHyDZc1+WWbyP7uJMj56w75pw7BWBmnzezTjPLmNlBM7tt5sPM7M/N7J9zizMN2+ncpcCbzewjZvZU3vbOzP6bmbWb2Wkz+0LuUUiYWamZ/bWZ9ZnZ62b28dz2b3tmrpn9E9kZ+r+X29enzGxN/vZm9hMz+19m9nRum++ZWYOZfT33XQ7kztrNfOZGM3vYzFJmdsTM7p3b/2oRKWZqwESkYJxz42QfMXR7btXtZB9D9NQ56/LPfh0g25wlgH8B/t3MKgM+fub3L3XO1Trn9p2njDuBbcA1wL3AL+TW/y6wJ7evG4BfucD3+BBwArgrt6//c55NPwB8CFgFXAnsA76S+y6HyT46CzOrAR7Ofb9lud/3t2a2+Xw1iIjf1ICJSKE9ztlm6TayDdiT56x7fGZj59w/O+f6nXOTzrm/BpYAG+aw/0875047506QPfs2c+btXuDzzrku59wA8Ok57GPGV5xzx5xzaeAhsmf2HsmNUft34PrcdncCbzjnvpL7ns8DDwC/Pg81iEgRUgMmIoX2BHCrmSWAJudcO/A02bFhCWALeWfAzOyTZnbYzNJmdhqIk334+uXqzns9AtTmXq8EOvPey399uXryXo8GLM/sezWwPXdZ9HTue34QWD4PNYhIEXrb2AYRkZDtI9tE/S7wUwDnXMbMTuXWnXLOvQ6QG+/1KWAX8LJzbtrMBgAL+Ny53tL9JtCSt9x6ke3n8xbyTuBx59x75/EzRaSI6QyYiBSUc24UeA64n+ylxxlP5dblj/+qAyaBJFBmZn8GxM7z0UlgGlh7maV9E/iEma0ys6XAH1xk+5457Otc3weuMrMPmVl57tc2M9s0T58vIkVGDZiIROFxsoPNn8pb92RuXX4D9iPgh8BrQAcwxnkuDTrnRoC/BH6au4y34xJr+gfgx8B/Ac8DPyDb/E2dZ/v/DfxJbl+fvMR9vYVzbhB4H9nB96fIXib9DNnxbiKy0tzFEAAAAIZJREFUAGkiVhGRAGa2B/g759zqqGsRkYVHZ8BERAAzqzKzXzSzMjNbRXaKiO9EXZeILEw6AyYiAphZNdlLoxvJ3qH4n8AnnHOZSAsTkQVJDZiIiIhIgekSpIiIiEiBqQETERERKTA1YCIiIiIFpgZMREREpMDUgImIiIgUmBowERERkQL7/6zRpkkV4QRjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "table['prob_get_in_bus'].plot(ax=ax)\n",
    "ax.axvline(kmf.median_, linestyle='--', c='orange', label='median')\n",
    "\n",
    "plt.xlabel(\"Waiting time\")\n",
    "plt.ylabel(\"Probability of get in bus\")\n",
    "plt.legend();\n",
    "\n",
    "print('Median: {:.0f} minutes'.format(kmf.median_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the survival probability, I am more likely to get in the bus after waiting for 34 minutes. It surely is crowded at the city, eh?\n",
    "\n",
    "Given this kind of analysis, we could suggest to add more buses so that passengers don't have to wait for such a long time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "- [Lifelines documentation](https://lifelines.readthedocs.io/en/latest/Quickstart.html#kaplan-meier-and-nelson-aalen), accessed at November 15, 2018."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Extras: hazard function\n",
    "# # to estimate probability of observing death event\n",
    "# from lifelines import NelsonAalenFitter\n",
    "# naf = NelsonAalenFitter()\n",
    "\n",
    "# naf.fit(df['wait_time'], df['observed']);\n",
    "# naf.plot(title=\"Hazard Function using Nelson Aalen Fitter\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Extras: weibull fitter\n",
    "# from lifelines import WeibullFitter\n",
    "\n",
    "# wf = WeibullFitter()\n",
    "# wf.fit(df['wait_time'], df['observed'])\n",
    "# # print(wf.lambda_, wf.rho_)\n",
    "# # wf.print_summary()\n",
    "\n",
    "# wf.cumulative_hazard_.plot(title=\"Hazard Function using Weibull Estimate\");"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}