examples/weibull_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# survival analysis\n",
    "\n",
    "A set of statistical approaches used to investigate the \"time to event\",\"churn prediction\"\n",
    "- censored data (special missing data type), so can't apply the regression directly\n",
    "- survival function: $ S(t) = Pr(T > t) $\n",
    "- hazard finction: $ \\lambda(t) =\\lim_{\\delta t\\to0}{\\frac{Pr(t\\leq T\\leq t+\\delta t\\mid T>t)}{\\delta t}}= -\\frac{S^\\prime(t)}{S(t)}=\\frac{f(t)}{S(t)} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.134783454153106 0 0.9980276199467129\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VOX9/vH3J3vYEkJCCCSQAIEQQFYBcUVRNhGtiqC2uFTUqq3ab61Wq5bW2qoV677WfQWVRRFEwBVRwg5hD1sSICHskD3P748Z/I0xkAmZmWeWz+u6cjnLmZmbw/Hm5DxnniPGGJRSSgWXMNsBlFJKeZ6Wu1JKBSEtd6WUCkJa7kopFYS03JVSKghpuSulVBDScldKqSCk5a6UUkFIy10ppYJQhK0PTkxMNOnp6bY+XimlAtKSJUv2GGOS6lvOWrmnp6eTk5Nj6+OVUiogicg2d5bTwzJKKRWEtNyVUioIabkrpVQQ0nJXSqkgpOWulFJBSMtdKaWCkJa7UkoFIS13pZQKQlruSikVhKx9Q1UFrslzNzTq9Xec38VDSZRSx6N77kopFYS03JVSKghpuSulVBDScldKqSCk5a58xxhSDq6AhU/Blm+gqtx2IqWClp4to7zLGBKPbiSreA5d9swlrnwnrHI+FxEL7QdBx7Mh4yxI6Q1h4VbjKhUstNwDUCCcihhXuoOsPXPoWvw5rUq3UEM421oO5Pv2NzF89FgoXA5bvoYtX8EXDzpeFBMH6WdCxtmOwk/sAiJez6pUMNJyVx7TtLyYLnvmkrVnNm0OrwUgv0Uf5qXczYbE8yiLjAdgeIu20KItZI10vPBwkaPo8750/HfdJ47Hm7Vx7NF3PNtR+PFpFv5USgUmLXfVKDGV+8ksmU/X4s9JPbgUwbC7aRZfp/+e9Ynnczi6Tf1v0qw19LzM8QOwbyvkfeUs/AWw6gPH4y0zHEXf6TzIuhDCdMhIqePRclcNFll9lI57vyareA4d9n9PuKlmb2wHFqXdwPrE89nXJL1xH9AyHfqlQ78JYAwUrXUcvsn7ClZ/BEtegy7D4ZLnIbZl4/9ASgUhLXfVMKumcuOPtxBZU8ahqNYsSxnPuqThFDf10vFxEUjOdvwMuhmqqyDnfzDnL/DC2XDFm5DSy/Ofq1SA03JX7ls5BT6eSFHzU/iu/e8oaNELxMeHRsIjYOBEaNsHpkyAl8+HUf+Bvr/2bQ6l/JwetFTucRY77QfzUfaTFMT18X2xu0o7FW78GjqcBjNuhRm3QWWZvTxK+Rndcw9BDT2VMqvoM4ZtfJCCFn2YlvJPqsJjvZSsgZomwtUfwYJ/wjePwc4VMPYNxzF7pUKc7rmrEzpW7PlxfZmWPdl/iv2YsHA4768w/n3HWTYvnAUb5thOpZR1Wu7quLKKZv1U7NO7Pe5/xe6q63CY+BXEt4d3xsL8h6Cm2nYqpazRcld1chT735zF7od77HVJyIDr50Lvq+HrR+CtS+FIie1USlmh5a5+IatoFsN/2mOfTFV4jO1I7ouMhYufgdFPwraFjsM0+Tm2Uynlc1ru6me6OYt9R1y/wCt2V/0mwPVzHN9i/d9w+PElxxeilAoRWu7qJ92KPmXYxgfZEdc/sIv9mLZ9HMfhOw2BWf8HH98IFUdsp1LKJ7TcFXCs2P/G9rhTnYOnAV7sxzRJcJxJM+ReWPkBvDwU9myynUopr9NyV2QXffJTsc/o9p/gKfZjwsLg7Lvg6g/h0C54aQgULLGdSimv0nIPcdlFn3DBxknOPfYgLHZXnc+DG7+C2Hh45wrYt812IqW8Rss9hGXvnvmzYq8O5mI/Jr49XDUVqivg7cuhdL/tREp5hVvlLiLDRWS9iGwSkbvreL69iCwQkWUislJERno+qvKk7N0zuGDT39kePyB0iv2YpK5wxduwNw/evxqqKmwnUsrj6i13EQkHngFGANnAeBHJrrXYfcAHxpg+wDjgWU8HVZ7jKPZ/sC1+INOzHgutYj8m40wY8zRs/QZm/kFPk1RBx5099wHAJmNMnjGmAngPGFNrGQO0cN6OAwo9F1F5Usbeb38q9hlZj4ZmsR/Taxyccw+seAe+ftR2GqU8yp1ZIdsBO1zu5wMDay3zIPC5iNwGNAWG1vVGIjIRmAjQvn37hmYNGo29wPXJiqw+yrmb/0VJk45Wi92vLvB99p8dE44teAjiO0CvKzz33kpZ5M6ee12X16n9O+x44DVjTCowEnhT5JeTfRtjXjTG9DfG9E9KSmp4WtUog7a/RIuK3czrdE9o77G7EnFMVZB+Jky/BbZ+azuRUh7hTrnnA66XnU/ll4ddrgc+ADDGfA/EAImeCKg8I/HIRvoWvsuq5DEUttDL0v1MRJTjcn0JGfDelVBs5zcrpTzJnXJfDGSKSIaIROEYMJ1Ra5ntwHkAItINR7kXezKoagRTw3mbH6YsojnfdLjNdhr/FNsSrpoC4VHw9mVwWDdfFdjqLXdjTBVwKzAHWIvjrJg1IjJJRC5yLvZH4AYRWQG8C1xjjJ5+4C967J5O20Or+Drjdsoj42zH8V8t0x1TFRwugnfHQWWp7URKnTS3LrNnjJkFzKr12P0ut3OB0z0bTXlCbMVeztz2NDta9GVtkn79oF6p/eDSl+D9X8NHE+Hy1x3TFygVYHSrDXJnbf0vkdVHmd/pbsfgoapft9Ew7CFYOwO+uL/+5ZXyQ3qB7CCWuj+H7OJZ/JB6HXubZNiO4zfcOhXTDOWclOX0WfgU83Y1YWXKZT895dFTMU+SX51OqvySlnuQCq+p4Ly8f7E/ph0/pF5rO07gEeGrjDuJKytkSN6jHIxOYWuCHnlUgUMPywSpfgVvklC6jfkd79Jz2k+SkXBmdX2I4qaZjFp/D0mH19uOpJTbtNyDUFzpDgbu+B8bWg1lW8vBtuMEtMrwJkzvNpmyiBZcvPYOmpXvth1JKbdouQcbYzg37xGqwyL5MuNO22mCwpHoJKZlP0Fk9VEuzr0Dyg7ajqRUvbTcg0yXki9I37+Ihe1v5ki0TvHgKSVNO/NJ1r9JKM2DKddAdZXtSEqdkJZ7EImqOszZeY+zu2k3Vric3aE8Y3v8QOZ1ugc2z4Ov/mU7jlInpOUeRAZvf56mlSXM63Q3RsJtxwlKa5LHQO+r4evHdJIx5df0VMggkXwol947P2B5yuXsbl77WirBxdaUyT8Z8W/Y/r3jG6w3f+eYl0YpP6N77kFATDXnbX6Yo5EJLGx/s+04wS+6GVz6MhzerVdxUn5Lyz0I9No5leQj6/gy404qIprZjhMa2vWFc/8KudNh2Zu20yj1C1ruAa5peTGDtz/H1vhBbEg833ac0DL495BxFnz2Z9iz0XYapX5Gyz3Anb3lccJrKpnf8S6dGMzXwsLgkhcgIgamXgdV5bYTKfUTLfcA1mHfQrqWfMGPaddyIDat/hcoz2vRFsY8DbtWwvy/206j1E+03ANUeHUZ5+Y9wt7YDuS0+43tOKEtaxT0vx4WPgWb59tOoxSg5R6wBuS/SnxZAfM63k11WJTtOOqCf0BSFnx8ExzZYzuNUlrugajl0a2cWvAGuUkjyY/vbzuOAohqApe+AqX7Yfotenqksk7LPdAYw3mbH6YyLJav0/9gO41y1aYHnD8JNsyGxS/bTqNCnJZ7gOlS8gVpB5fybfqtlEYl2I6jaht4I3Q+H+bcC7tzbadRIUzLPZAYw6n5r1ISm86q5Ittp1F1EYGLn4OYOPjweqgstZ1IhSgt9wDSYf/3tD6y0XF2jOhfnd9qluQo+KJcmKsX2FZ2aEMEkFPzX+dQVGvWJQ23HUXVJ3MoDLoFfnwR1s+2nUaFIC33ANHm0CrSDi5lSdurqAmLtB1HuWPoA5DcE6b/Dg7tsp1GhRgt9wBxav4blEW0YHUbPdYeMCKi4bJXoOKo4/z3mhrbiVQI0XIPAAlHt9B575csb3M5leFNbMdRDZHUFYY/DHkL4PunbadRIUTLPQD0K3iTyrBolre9wnYUdTL6XQNZF8K8SVC4zHYaFSK03P1cs/LddCv+jNXJYyiN1Cv+BCQRuOgpaJoEU6+H8sO2E6kQoOXu5/oWvoMYw9K2V9mOohqjSQL86gXYmwez77adRoUALXc/Fl15gJ67prE+6XwOxrS1HUc1VsZZcMYdjis3rZ1pO40KclrufqzXrqlE1RxlcbsJtqMoTxnyF2hzCsy8XWePVF6l5e6nIqrL6FP4Hnktz6CkaWfbcZSnhEfCJc9D2QH49I+206ggpuXup7oXzaBJ1X69EEcwSu4O59wNudNg9Ue206ggpeXuh8RU0a/gbQqbn0JBi9624yhvOP12aNvXsfd+uMh2GhWEtNz9UJc9XxBXXsjidr/Ri14Hq/AIx+GZiiOO4+96cQ/lYW6Vu4gMF5H1IrJJROo8j0tExopIroisEZF3PBszhBjDqflvUBKbQV7CmbbTKG9K6grn3gfrP4WVH9hOo4JMveUuIuHAM8AIIBsYLyLZtZbJBO4BTjfGdAdu90LWkJC+fyFJRzeSk6rT+oaE026BtIHw2Z/g4E7baVQQcac9BgCbjDF5xpgK4D1gTK1lbgCeMcbsAzDG6EHEk3Rq/uscjEpmXeIw21GUL4SFw5hnoaoCZv5eD88oj3Gn3NsBO1zu5zsfc9UF6CIi34nIIhHRCcdPQsrBlaQeXMbSdlfqtL6hJLGzY3rgjZ/D8rdtp1FBwp1yr2tEr/buRQSQCZwDjAdeFpH4X7yRyEQRyRGRnOLi4oZmDXr9C96gNCKO1XoJvdAz4EbocDrMvgcO5NtOo4KAO+WeD6S53E8FCutYZroxptIYswVYj6Psf8YY86Ixpr8xpn9SUtLJZg5Kjml9v2J5ylid1jcUhYXBmGegphqm36qHZ1SjuVPui4FMEckQkShgHDCj1jLTgCEAIpKI4zBNnieDBrv+BW84pvVNGWs7irIlIQMumOSY+33Jq7bTqABXb7kbY6qAW4E5wFrgA2PMGhGZJCIXORebA5SISC6wAPiTMabEW6GDTbPyXWQVz2Z18sWURf7iaJYKJf2ug4yzYc59sG+r7TQqgLl1rp0xZpYxposxppMx5iHnY/cbY2Y4bxtjzJ3GmGxjTE9jzHveDB1s+ha+ixjDEp3WV4WFwZinHafBTr9VL82nTpqeSG2ZY1rfj1mXNIxDMSm24yh/EN8ehj0EW7+BxS/bTqMClJa7Zb13TSGqppScdr+2HUX5k76/gc5D4YsHoGSz7TQqAGm5WxRRXUbvwvd1Wl/1SyIw+kkIi4Rpv3OcRaNUA2i5W3RsWt/FqdfYjqL8UVw7GPFv2LEIFj1nO40KMFruloTVVNGv4C0KmveisEUv23GUv+o1DrqMgPl/h+INttOoAKLlbkmXPXOJK9/J4lS9hJ46AREY/V+IjIVpN0N1le1EKkBoudtgDP0L3mBPk45saXm67TTK3zVPhpGPQUEOfP+U7TQqQGi5W5C+byFJRzeR026CTuur3NPjUuh2ESz4J+zOtZ1GBQBtFgtOLXidg9FtWJ94ge0oKlCIwKjHIbo5TLuJsBo9PKNOTMvd13b8SOrBZSxpexU1YRG206hA0iwJLpwMO1cwIF/nnlEnpuXuaz++RFl4M1Yn177eiVJuyB4DPccyIP8Vkg/p4Rl1fFruvnR0L+ROZ23rkVSFx9pOowLVyEc4GtmKYRsfJLy6zHYa5ae03H1pxbtQXc6q5EtsJ1GBLLYln2feT6vSLZy+Xb/cpOqm5e4rxkDOq5A6QKcaUI22PX4gy9tcTr/Cd0jdn2M7jvJDWu6+sm0hlGyE/tfaTqKCxDfpt7Evpj3DNv2NqKrDtuMoP6Ona/jKklchOg6yL4av9BqZgWzyXP+YBqAqPJbZXR7kipW/5ewtjzM3837bkZQf0T13XzhSArnTHfOEROn1UZXn7Grek8WpE+hRNJOOJV/ZjqP8iJa7L6x4F6oroN81tpOoILQo7QaKmmYydPM/ia3cZzuO8hNa7t5mDCx5DdIGQnK27TQqCNWERTI7cxLRVYc4b9M/HducCnla7t627TvHQGo/HUhV3lPStDML299M5t4v6Vb8me04yg9ouXtbzqsQEwfdL7adRAW5pe2upKBFb4bkPUKz8l224yjLtNy96UgJrJ0BvcY75uNWyouMhDMn8wHE1DBs4yQwNbYjKYu03L1pxTs6kKp86kBMKl9l3EH7A4vpvXOK7TjKIi13b/lpIHUQtO5mO40KIauTLyav5emcue0pWh7dajuOskTL3Vu2fgslm3SvXfmeCF90vo/KsBiGb3wAMTr3eyjScveWJa9CTLwOpCorjkQlMr/Tn2lzOJcB+a/ZjqMs0HL3hiN7IFcHUpVdGxLPZ13iMAbueJnWh9fZjqN8TMvdG5a/AzWVekhGWTe/458ojUxg+Ib7Ca8ptx1H+ZCWu6cdG0htfxq0zrKdRoW48sg4Pu/8V1qVbmHwNp37PZRouXva1m9g72bda1d+Y1vL01jR5lL6Fb5DuwNLbMdRPqLl7mk5zoHUbL1GqvIfX6f/gQMx7Ri2cZLO/R4itNw96XAxrJ2pA6nK71SFxzI780Gal+/irC1P2I6jfEDL3ZNWOAdS9WpLyg/tbNGLnNTf0LNoumMnRAU1LXdPqalxDqQOhqSuttMoVafv0yayq1k3mH4L7NtmO47yIi13T9n6DezN04FU5ddqwiKZ1dU55/vUa6GqwnYk5SVa7p6yRAdSVWA4EJMKY56GgiXwxYO24ygvcavcRWS4iKwXkU0icvcJlrtMRIyI9PdcxABwuBjWfgK9r4TIGNtplKpf9hgYMBEWPQPrPrWdRnlBveUuIuHAM8AIIBsYLyK/uF6ciDQHfg/84OmQfm/52/qNVBV4LvgHpPSCaTfD/u220ygPc2fPfQCwyRiTZ4ypAN4D6jr28HfgEaDMg/n8nw6kqkAVEQ2Xv+Y4/j5Fj78HG3fKvR2ww+V+vvOxn4hIHyDNGPOJB7MFhq1fw74tevqjCkwJHeGiJ6EgB+b9zXYa5UHulLvU8dhPl1cXkTBgMvDHet9IZKKI5IhITnFxsfsp/VnOqxDbErpdZDuJUien+yVw6m/h+6dhvV5cO1i4U+75QJrL/VSg0OV+c6AH8KWIbAUGATPqGlQ1xrxojOlvjOmflJR08qn9xeEiWPcJ9NKBVBXgLngI2pwCH98E+3fUv7zye+6U+2IgU0QyRCQKGAfMOPakMeaAMSbRGJNujEkHFgEXGWNyvJLYnyx/G2qqdCBVBb7IGMfx95pqmHodVFfaTqQaqd5yN8ZUAbcCc4C1wAfGmDUiMklEQvdYxLGB1A6nQ1IX22mUarxWneCi/0L+jzBvku00qpEi3FnIGDMLmFXrsfuPs+w5jY8VALZ8Bfu2wpD7bCdRynN6XOq4/u/CJyH9DOgyzHYidZL0G6ona8mrEJsA3UbbTqKUZw17GJJ7wsc3woF822nUSdJyPxmHixzf6tNvpKpgFBkDY193HHfX4+8BS8v9ZCx7yzGQ2neC7SRKeUerTjD6v7DjB5j/D9tp1EnQcm+omhpY+jp0OEMHUlVw63mZ40yw756AjXNtp1ENpOXeUFu+dAyk6umPKhQM/xck94CPJsKBAttpVAO4dbaMcrH0Tec3UnUgVQWuyXM3uL1sy7YPcGXxbyh+5Uqm9HgOI42vjTvO1996vU333Bvi6F7HN1J7jtWBVBUy9jVJZ16ne2h3cDmDt79gO45yk5Z7Q6yaCtUV0Odq20mU8ql1rUewKnkMA/Jfo8O+723HUW7Qcm+IZW865t9IOcV2EqV8bkHG/1HcpDPDN9xPs/JdtuOoemi5u2vnSti1Evr82nYSpayoDo/h064PE24qGbP2TiKrjtiOpE5Ay91dy9+G8CjH6WFKhah9TdL5tOvDJB7JY9SGvyCmynYkdRxa7u6oKoeV70PWhdAkwXYapaza1vI05ne6i4x9CxmS95jjSk7K7+ipkO5YPwtK9+lAqlJOq9r8iriyAk4teIP9MWksbXeV7UiqFi13dyx7C1qkQsdzbCdRym982+EW4soKOGvrfzkQ05bNrYbYjqRc6GGZ+hzIh03zoPd4CAu3nUYp/yFhzM58kJ3NezBiw19JPrTGdiLlQsu9PiveBYxjBkil1M9Uh8cwI+sxjka2YszaO2lRVlj/i5RPaLmfSE2N45BM+pmOq8QrpX6hNCqBadlPEF5TycW5txNddch2JIWW+4ltX+iYJEwHUpU6ob1NMpjZ7VHiy3Zw4bo/E1ajc8DbpuV+Isvegqjm0C10LxWrlLvy4/rxRed7aX9gMedtflhPkbRMz5Y5nrKDsGYa9LoCoprYTqNUQMhtfSFxZQUM2vEy+2PSWJx2re1IISsky92d6U577JrG+VWlvFt5NrtqLa/TlSp1fN+nTSSurIAztj/LwZi2rE/Si2zbEJLl7o7uRTMoic1gV7PutqMoFVhEmNv5PpqX7+KCjX/jUHQyhS16204VcvSYex0Sjm6h7aFVrE6+CERsx1Eq4FSHRTEj61EORqdw0dr/I650h+1IIUfLvQ7di2ZSLeGsSxphO4pSAas8Mo5p2U9gEC7JvZ2Yyv22I4UULfdawmqq6Fb0KVtansnRqFa24ygV0A7EpjGz22M0L9/F6HV3EV5TYTtSyNByryV933c0rdzLmtYX2o6iVFAobNGLOZkPkHpwGedvnKSnSPqIDqjW0r1oJkciE9ja8nTbUZQKGhuSLiCuvIAztj3LgZhU4HHbkYKelruLJhV76Lj3W5a2u5KasOOvmoZcOV4p5bC43TXElRYwKP8VFrwez/K24076vfR05PppubvoVvwZYVSzpvVo21GUCj4izO90NzHVBxmy5T9E1pSxOPUa26mClh5zP8YYuu+eSWHznuxtkmE7jVJBqSYsgk+7/pO1ScM5Y9sznLbteT0G7yW65+7U5vBqWpVuYW6ne21HUSqoGYlgTuaDVIVFMyj/FSJrSvk6/Xb9TomHabk7dd89k8qwGDYkDrUdRamgZyScLzrdS1VYDP0K3yGippz5He8C0YMJnqLlDkRUl9F1z+dsSDyPiohmtuMoFRpE+DLjj1SGxzIg/zUiasqY2/k+jGgteYKuRSCzZD7R1UdY01qn9lXKp0T4rsMtVIbFcvr254ioLmd2l7+f8Gw15R5dg0D33TPYF5NGQYs+tqMoFZJ+TLuOyrAYztk6mYh15Xya9TDVYdHHXb6xpyOHwqmUbh3gEpHhIrJeRDaJyN11PH+niOSKyEoRmSciHTwf1TviSvNJO7iE3NYX6oCOUhYta3cl8zreTad93zBm7R+JqC61HSmg1VvuIhIOPAOMALKB8SKSXWuxZUB/Y8wpwFTgEU8H9ZbsopkYhNzWo2xHUSrkrUy5lNmZD5K2fzG/yv09UVWHbUcKWO7suQ8ANhlj8owxFcB7wBjXBYwxC4wxR513FwGpno3pHWKq6V70CVvjB3E4Otl2HKUUsLb1KGZ1fYg2h1Zx6ZpbiK48YDtSQHKn3NsBrpMx5zsfO57rgc8aE8pX2u//keYVRaxJ1oFUpfzJxsShfJL1CIlHNnL56puJrdhrO1LAcafc6zoQXedXykTkaqA/8Ohxnp8oIjkiklNcXOx+Si/pvnsGpRFx5CWcZTuKUqqWvISzmJ49mfiy7Vy++kaalhfZjhRQ3Cn3fCDN5X4qUFh7IREZCtwLXGSMKa/rjYwxLxpj+htj+iclJZ1MXo+JrjxAp71fsS5pONVhUVazKKXqtj1+IB9nP0WzimLGrp5Ii7JfVI86DnfKfTGQKSIZIhIFjANmuC4gIn2AF3AUe0D885q1Zw4RptJxKT2llN8qiOvDh92fIabyIJevmkh86XbbkQJCveVujKkCbgXmAGuBD4wxa0Rkkogca8ZHgWbAFBFZLiIzjvN2fqP77hnsbprFnqbBf76rUoFud/PuTOn5PBGmgrGrbiD5UK7tSH7PrfPcjTGzjDFdjDGdjDEPOR+73xgzw3l7qDEm2RjT2/nj17vDSYfXk3xkPWuSdWpfpQLFnqZdmNLjBaolkrGrfkuvnR/ojJInEJKz9HQvmkmVRLEucZjtKEqpBtjbJIO3e7/F9vgBnJv3KCPX/0XPhT+O0Cv3yjK6FX/G5lZnUx4ZZzuNUqqByiLjmd7tcb7tcAuZJQu4csUEEo/o1dFqC71yX/8pMVUHWa2ThCkVuCSMxanXMLXHs0RWH2X8yuvosWuaHqZxEXrl/sOLHIhuy474U20nUUo1UkFcX97q/TYFzXtx/uaHGLbxQZ2Txim0yr1wGexYxLKUKzASbjuNUsoDSqMS+Lj7k3yfNpFuxZ9x5YoJJBzNsx3LutAq90XPQ1QznW5AqSBjJJxF7W/go+5PEVN1gCtXTCCraJbtWFaFTrkf2g2rP4TeV+nVlpQKUtvjB/J2r7fY3awbIzY+wNBN/yC8usx2LCtCp9xz/gc1VTDwRttJlFJedCQ6iak9nuXH1GvouXs641deR3zpNtuxfC40yr2qHHJegcwLoFUn22mUUl5mJILvOtzCx92eoFlFEVeumECXPXNtx/Kp0Cj31R/CkWIYdJPtJEopH9qacDpv936LkiYdGbX+LwzZ/AjhNRW2Y/lE8Je7MbDoOUjKgo5DbKdRSvnYoeg2TOnxIjltr6L3rilcsfK3ULTWdiyvC/5y3/497FrpONau10hVKiTVhEXwTcbtzMh6lLjyAnjudPj0j3CkxHY0rwn+cl/0HMTEwynjbCdRSlm2udU5vNr3Q+h/HeS8Ck/1cXREdaXtaB4X3OW+bxus+wT6XQNRTWynUUr5gbLIeBj1GNz8HbTrB7PvhmdPgw1zgmr6guAu98UvAQIDbrCdRCnlb1p3g6s/gis/AAy8MxbeuhSK1tlO5hHBW+4VR2DpG9BtNMSl2k6jlPJHItBlGNz8PQx7GApy4LnBMOtPcDSwL8odvOW+4l0oOwCDbradRCnl7yKi4LTfwW3LoP+1sPhleLKPY8qSAD0eH5zlXlMDP7wAbftA2kDbaZRSgaJpKxj1H7jpO2jbG2b/2Xk8/nPbyRoswnYAr8ibD3s2wCUv6umPSqmGS86GX0+DDbPpWhn1AAAJVElEQVRhzr3wzuXQeShc8BC0zmLy3MZdHOSO871/7ebg3HNf9Bw0S4bul9hOopQKVCLQdQT8bpGj1Hcsdh6Pv4um5cW209Ur+Mq9eANs+gL6X+84jqaUUo0REQWDb4XfL4V+E2DxS/w2ZzSj1t1N6oElfnv6ZPAdlvnxBQiPcgyKKKWUpzRNhAsnw+DbWDr1P3QvmkmXknmUxGawIuUy1iaN9KvpxINrz710Pyx/F3pcBs1a206jlApGCR35JuMPvNT/E+Z0vp/K8FjOzXuUGxaP5NzN/6LVkU22EwLBtue+7E2oPKKzPyqlvK46PIbc5NHkJo8m+dAaeu2aSveiT+i160PyW/RhRZvL2NRqCDVhkVbyBU+5V1fBDy9Ch9MhpZftNEqpELK7eXc+b96dr9P/QLaz4EdtuJcjkQmsSr6YVW1+xeHoZJ9mCp5yXz8LDmyHYQ/ZTqKUClFlkfEsbXc1S9teSYf9i+i1cyoD819lQP7rbE44k5Upl7E9boBPsgRkudd1junlqx6neXQKr+ZnYgoadw6qUko1ioSxreVgtrUcTIuyQnru+ogeu6eTufdL9sZ2gA6PQNfhXo0QFAOqSYfXk3pwGStSLsdIuO04Sin1k4Mxbfku/VZePvUTPsv8G2URLSDc+8fhA3LPvbY+O9+jIiyW1cljbEdRSqk6VYdFs671SNa1HskdnTK9/nkBv+ceW7GXrsVzyG09ivKIFrbjKKVU/XwwLUrAl/spuz4kwlSyPOUK21GUUspvBHS5h9VU0mvXh2yJP419TdJtx1FKKb8R0OXeZc9cmlaWsKytXh9VKaVcBe6AqjH02fkee2M7sC1+kO00SqkA0tgpewNBwO65pxxaSZvDa1mWMg4kYP8YSinlFQHbin0L36MsvDm5rUfZjqKUUn7HrXIXkeEisl5ENonI3XU8Hy0i7zuf/0FE0j0d1FWz8l10LlnA6uQxVIXHevOjlFIqINVb7iISDjwDjACygfEikl1rseuBfcaYzsBk4N+eDuqq984pgGF5ylhvfoxSSgUsd/bcBwCbjDF5xpgK4D2g9ldBxwCvO29PBc4T8dJZ+hVH6bl7Gptbnc2hmBSvfIRSSgU6d8q9HbDD5X6+87E6lzHGVAEHgFaeCPgLK98npuqgYyBVKaVUndw5FbKuPfDaFw10ZxlEZCIw0Xn3sIisd+Pz65II4/ec5Gu9KRHQXO7TXA3nr9k0VwPc2bhcHdxZyJ1yzwfSXO6nAoXHWSZfRCKAOGBv7TcyxrwIvOhOsBMRkRxjTP/Gvo+naa6G0VwN56/ZNFfD+CKXO4dlFgOZIpIhIlHAOGBGrWVmABOcty8D5hvjp5cEV0qpEFDvnrsxpkpEbgXmAOHA/4wxa0RkEpBjjJkBvAK8KSKbcOyx6wFxpZSyyK3pB4wxs4BZtR673+V2GXC5Z6OdUKMP7XiJ5moYzdVw/ppNczWM13OJHj1RSqngE7DTDyillDo+vyv3xkx1ICL3OB9fLyLDfJzrThHJFZGVIjJPRDq4PFctIsudP7UHo72d6xoRKXb5/N+6PDdBRDY6fybUfq2Xc012ybRBRPa7POfN9fU/ESkSkdXHeV5E5Eln7pUi0tflOa+sLzcyXeXMslJEFopIL5fntorIKue6yvFUpgZkO0dEDrj8fd3v8twJtwEv5/qTS6bVzm0qwfmcV9aZiKSJyAIRWSsia0TkD3Us47vtyxjjNz84Bmw3Ax2BKGAFkF1rmd8BzztvjwPed97Odi4fDWQ43yfch7mGAE2ct28+lst5/7DF9XUN8HQdr00A8pz/bem83dJXuWotfxuOgXqvri/ne58F9AVWH+f5kcBnOL67MQj4wQfrq75Mg499Fo5pQH5weW4rkGhxfZ0DfNLYbcDTuWotOxrHGXxeXWdACtDXebs5sKGO/x99tn352557Y6Y6GAO8Z4wpN8ZsATY5388nuYwxC4wxR513F+H4PoC3ubO+jmcYMNcYs9cYsw+YCwy3lGs88K6HPvuEjDFfU8d3MFyMAd4wDouAeBFJwYvrq75MxpiFzs8E321bxz67vvV1PI3ZNj2dyyfblzFmpzFmqfP2IWAtv/w2v8+2L38r98ZMdeDOa72Zy9X1OP51PiZGRHJEZJGIXOyhTA3JdanzV8CpInLsC2l+sb6ch68ygPkuD3trfbnjeNm9ub4aova2ZYDPRWSJOL4BbsNpIrJCRD4Tke7Ox/xifYlIExwl+aHLw15fZ+I4XNwH+KHWUz7bvvztSkyNmerArSkQTpLb7y0iVwP9gbNdHm5vjCkUkY7AfBFZZYzZ7KNcM4F3jTHlInITjt96znXztd7Mdcw4YKoxptrlMW+tL3fY2L7cIiJDcJT7GS4Pn+5cV62BuSKyzrlX6ytLgQ7GmMMiMhKYBmTiB+vLaTTwnTHGdS/fq+tMRJrh+MfkdmPMwdpP1/ESr2xf/rbn3pCpDpCfT3Xgzmu9mQsRGQrcC1xkjCk/9rgxptD53zzgSxz/ovsklzGmxCXLS0A/d1/rzVwuxlHrV2Yvri93HC+7N9dXvUTkFOBlYIwxpuTY4y7rqgj4GM8dinSLMeagMeaw8/YsIFJEErG8vlycaPvy+DoTkUgcxf62MeajOhbx3fbl6UGFRg5IROAYSMjg/w/CdK+1zC38fED1A+ft7vx8QDUPzw2oupOrD44BpMxaj7cEop23E4GNeGhgyc1cKS63LwEWmf8/gLPFma+l83aCr3I5l+uKY3BLfLG+XD4jneMPEI7i5wNeP3p7fbmRqT2OMaTBtR5vCjR3ub0QGO7JdeVGtjbH/v5wlOR257pzaxvwVi7n88d2/Jr6Yp05/9xvAE+cYBmfbV8e3Qg8tIJG4hhl3gzc63xsEo69YYAYYIpzY/8R6Ojy2nudr1sPjPBxri+A3cBy588M5+ODgVXOjXsVcL2Pcz0MrHF+/gIgy+W11znX4ybgWl/mct5/EPhXrdd5e329C+wEKnHsLV0P3ATc5HxecFycZrPz8/t7e325kellYJ/LtpXjfLyjcz2tcP4d3+vJdeVmtltdtq9FuPwDVNc24KtczmWuwXGShevrvLbOcBwuM8BKl7+rkba2L/2GqlJKBSF/O+aulFLKA7TclVIqCGm5K6VUENJyV0qpIKTlrpRSQUjLXSmlgpCWu1JKBSEtd6WUCkL/D0FtbXqoIamqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Weibull disribution\n",
    "%matplotlib inline\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sample=stats.weibull_min.rvs(c=2,scale=1,size=300,loc=0) # c:shape parameter, scale: scale parameter\n",
    "#print(sample)\n",
    "\n",
    "#stats.exponweib.fit(sample,floc=0,f0=1)\n",
    "shape,loc,scale=stats.weibull_min.fit(sample,floc=0)\n",
    "print(shape,loc,scale)\n",
    "plt.hist(sample,bins=np.linspace(0,int(max(sample)),20),density=True,alpha=0.5)\n",
    "x=np.linspace(0,int(max(sample)),20)\n",
    "y=[(shape / scale) * (u / scale)**(shape-1) * np.exp(-(u/scale)**shape) for u in x]\n",
    "\n",
    "plt.plot(x,y,label='weibull')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kaplan-Meier estimate\n",
    "- a non-parametric statistic used to estimate the survival function from lifetime data\n",
    "- \n",
    "$ S(t)=\\prod_{i:t_i\\leq t}(1-\\frac{d_i}{n_i}) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   ID Retail date Production date  Group\n",
      "0   1  2014-12-09      2014-10-23      1\n",
      "1   2  2014-12-09      2014-10-17      1\n",
      "2   3  2014-12-14      2014-10-14      1\n",
      "       ID  Mileage Repair date  Censored\n",
      "0  1689.0    82441  2017-08-31         0\n",
      "1  2736.0    66951  2017-08-30         0\n",
      "2  1515.0    75500  2017-07-24         0\n",
      "   ID Retail date Production date  Group  Mileage Repair date  Censored  \\\n",
      "0   1  2014-12-09      2014-10-23      1      NaN  2018-12-19       1.0   \n",
      "1   2  2014-12-09      2014-10-17      1      NaN  2018-12-19       1.0   \n",
      "2   3  2014-12-14      2014-10-14      1  69879.0  2018-05-28       0.0   \n",
      "3   4  2014-12-15      2014-11-11      1  31563.0  2018-08-22       0.0   \n",
      "4   5  2014-12-16      2014-10-20      1      NaN  2018-12-19       1.0   \n",
      "\n",
      "   Operating_days  Operating_time  \n",
      "0            1471              48  \n",
      "1            1471              48  \n",
      "2            1261              41  \n",
      "3            1346              44  \n",
      "4            1464              48  \n",
      "    time  n_censored  n_event  n_risk  hazard_rate  survival_rate  survival  \\\n",
      "6      1         0.0     12.0   10091     0.001189       0.998811  0.998811   \n",
      "7      2         0.0      5.0   10079     0.000496       0.999504  0.998315   \n",
      "8      3         0.0     12.0   10074     0.001191       0.998809  0.997126   \n",
      "9      4         0.0     15.0   10062     0.001491       0.998509  0.995640   \n",
      "10     5         0.0     14.0   10047     0.001393       0.998607  0.994252   \n",
      "\n",
      "    failure_probability  std_error  lower_95_ci  upper_95_ci  \n",
      "6              0.001189   0.000343     0.998138     0.999483  \n",
      "7              0.001685   0.000408     0.997515     0.999115  \n",
      "8              0.002874   0.000533     0.996082     0.998171  \n",
      "9              0.004360   0.000656     0.994354     0.996925  \n",
      "10             0.005748   0.000753     0.992777     0.995727  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4FNX6wPHv2c2md1KANBIIJbQAAaSKgDQF7IANsWDv5erVq6g/rl7btXKFa7mKgoIFFFABqdJD7yEJpAKpJIT0zfn9MSsGCBAgyaa8n+fZh52ZM7PvGc27Z8/MnKO01gghhGg6TPYOQAghRN2SxC+EEE2MJH4hhGhiJPELIUQTI4lfCCGaGEn8QgjRxEjiF0KIJkYSvxBCNDGS+IUQoolxsHcAp/Pz89OtWrWydxhCCNGgbN68OUtr7V+dsvUu8bdq1YrY2Fh7hyGEEA2KUiqpumWlq0cIIZoYSfxCCNHESOIXQogmRhK/EEI0MedN/Eqpz5RSGUqpXWfZrpRS7yul4pVSO5RS3Sttm6iUOmB7TazJwIUQQlyc6rT4/weMOMf2kUCk7TUZ+A+AUsoXeAnoDfQCXlJK+VxKsEIIIS7deRO/1noVkHOOImOBL7VhPeCtlGoBDAeWaK1ztNa5wBLO/QUihBCiDtTEffxBQEql5VTburOtrxVaa+JyElmdvhxnszMuDi64OLjg7OBMj8AeeDl5UVZRhgkTZpO5tsIQQoh6ryYSv6pinT7H+jMPoNRkjG4iQkNDLyqIvKIyJk3/mONBv56xLVr/jWDPLhQ4rmVV9sf4OvvR3C2Qlu7NCXQL5M5Od+Ln4kdReREWkwUHU717rk0IIWpMTWS4VCCk0nIwkG5bP+i09SuqOoDWegYwAyAmJuaiZn93MJt4KTyEK/YnU2QyUaQUhSZFkVJMczzOktSjhFvjGOTpRIpDGanmdOIsyZRZCklK7E3n5sEcLJ/Hb2kzCXYPJsQzhFCPUII9ghnXbhyOZke01ihV1feZEEI0HDWR+H8CHlJKfYNxITdPa31YKfUb8M9KF3SHAc/VwOdVyd3JgeHjH4bSO3EsysWrKBeKcqEoh48irgAXb4q3ZKNWrcSSsxOTtgJgxcRYy2EWbi2gtWsOPm5RpB7XHD6WxDpzLGClOUOJDPDki7i3WZW2klCPUCK8Imjt3ZpIn0h6Nu9ZW9USQogad97Er5SajdFy91NKpWLcqWMB0Fp/DCwCRgHxQCEwybYtRyn1KrDJdqhXtNbnukh86ZQCJ3fj5R1yxmbn7uOh+3goL4WcBMjchzkrngUDx5FXVE7Z99/hl7AQjsNhcxA7KkJZX96Se/ZuMfb3NuPpG05cYTY7MhZSqk/Qwi2IxTcY3Uvvb3mfEmsJ3QK6ER0QjZ+LX61WVwghLobS+qJ6VmpNTEyMttsgbQWZcHib8UrfBkd2YjU7s23MbyRkFNB57SOYCzPYbQ3hj6IQNptakGJyw9exNR1bepLpMoP00q2U61IAwjzDGNN6DJO7TAaQriIhRK1RSm3WWsdUp6xcxazM3R8irzReNubyUno4ONIjzAfyu8OhP2h79A+utRwHINlvAO8FTGV3eh6+8e1JsQ6myCWPoBZHsJakkJh1nKJSKw4OFYz4bgRtfNrQ2a8zXfy70NmvMz7O8miDEKJuSYv/YlRUQFYcpMWCszd0uBpKC9GvBQOaw+4dWaWjmZ3bnh3WMCxmB7qGWajwWcgJEjlcdIgKXQHAs72e5ZYOt5Bfms/mI5uJ9ImkpXtLTEpG0xBCVJ+0+GubyQQB7Y3Xn8wW1B0LIHElLeOXMD7tK8ZbNPG9/s4cyxg2x6eStKUXWXooZnMpEcG5BPgdpfh4GBn5xcQf38Ujyx8BwNXBlTY+bYjwiuDuzncT5hnGgdwDrExdiZPZCSezEz7OPgS7B9PauzWOZkc7nQghREMkLf7aciIL4n+H0N7g0wp2/whz76DAuz37XLrxe3F7vs0MIafcGYBmHorwFnn4+GRjcjrC8YoUjhSm8N4V79HZvzPz4ufxjzX/OONjvhv9He1827EseRm/HvqVUI9Quvh3oat/V7ycvOq40kIIe7mQFr8k/rqSnWAk/4MrIXkDWEvQyszOG1YRm+vG4UN72JdRxNosF6xGLxAuFjPuzg44W0w4OZhwsmgcHaw4Wipwdy3EwyOf/kH9aRfQjB15v/HVvs85fOIwVtutqm282/DVqK9ws7hRVF6Es9lZLi4L0UhJV0991Kw1DHzKeJUVQ8oGVOomukR1pItSkPEa5H+L9vHnuG8XDjm3Yxsd2OMcTXGZFdeiw5SWl1NQbiGv2IF9mR6k51uYp/faPsANP/cnifR3IKxlNk7uSZSZMnCzuAHwjzX/YNORTUQ1izr56tisI83dmtvvnAgh7EJa/PXFkV2QvA7StkDaZuPicUAHeGCdsf2ToZC66a/yFjeskcOJH/g+B7NOUBy/ioQCZ9Ye82RbeiHWCo2j2UR0qDd9IpqhPLZwpHQH+3P3kpiXSIWuIKpZFN9e/S0Ac+Pm4u3kTXvf9gS7B8svAyEaGOnqaQyK8yEvFQKjjOX4pXD8CJQVQVkh5KeDqx9c/rSx/c1IOJEBZkfKg3uT5N2bpRU9WJDuya70PLQGi1kR1dKLzsEuBPplE+7vxFWRAwDo/01/8kvzAfBw9KC9b3tGR4zm2shrAXkGQYj6ThJ/U5S+FTLj4MgOSFwBR3dBn4dg+FTyCgpJWTWTFWUdWH3Ewo7UPIrKjOsAfu5OdAv1plOwK/4+uVgdUziUf4C9OXu5IuQK7ulyD8dLjzP6x9F09utMZ3/jGYROzTrh7uhu3zoLIU6SxC+MXwe6AjxbQtI6+Nw2FYJ3GBWhfTns1Y21Dj1Zd9TEtuRjJGadAIxRL1r7u9M12JvoEC+6h/ng61nMh9veZ2fWTg7mHTTKoZjafyqjW4/meOlxMgozCPcKl+cPhLATSfziVBUVxi+BpDWQtNa4llCYDXcsglb9IG0LRQlr2evYkbUFzdmaWsC2lGNknzCGnvBwciCmlQ89w33pGGRBOyWzN2cXI8JHEO4VzqLERfxt9d/wsHjQ2b8znf0608mvE72a98LV4mrnygvRNEjiF+emtXHx2KcVODjBqjdh2f8Z2yxuEByDDulNWsd72ZRezMaDuWw8mE1CpvGrwMnBRHSIN73DfekV3owQ/1K2ZW1ie+Z2tmduJ/5YPBW6gsXXL6aFewtWpa4iLjeO7gHd6RbQTa4VCFELJPGLC5eXBinrIXm98YsgLw2eTjCeUl73ERTnkde8L+vLItiQVMDGQ9nsSc+nQoODSdEpyMv2ReBLVJAzR4sP0tW/K0op3tj0BjP3zASgnU87JnacyIhWI7CYLXautBCNhyR+cenKS4xfAwBzJsKe+YAGiyuE9YWosRyPmsDmpFw2Hcph48EctqfkUWp7+qxNgDsxYT50D/MhJswHX49ylqcs54vdX5CQl0Dv5r35ZPgn9qufEI2MJH5R84py4dAfkLjSePo4qAdc+7HRbfTzo9AympKQAWwt8GFz8jE2J+WyOSmXvKIyAHxcLfRt7cfdA8IpMO1CoxkYPJDCskI+3v4xw1oNo4NvB5kPWYiLJIlf1D5rGZgtUJABMwZBfpqx3isUwgdAzJ1UtOxBYlYBm5NyiT2Uy5K9RzlWWMblbf15ZEgkPcJ8WJu2lgd/f5ByXY6HxYOY5jH0btGbkeEj8XX2tWsVhWhIJPGLuqU1ZMcbzw8cXAmH1sDodyFqrPFE8vppENaPE8H9+XJ3Of9dnUjOiVL6t/Hj0aGRtG4OGw5vYMPhDWw8spGU4yksvHYhoZ6hbD66mbSCNPq27CszmglxDpL4hX1VVBjPEJgdYN9CmP8QFOWAMkG3Wynq9zQzd5cxY1UiWQWl9IloxvheIQyI9MfXzZHDBYdp7tYcpRRT1k7h+wPfA9DBtwP9g/rTL6gf3QO6y91BQlQiiV/ULxUVkLkPtnwBmz4FJw94Yg9F2pFZG5OZvjKBjOMlKAVdgr0Z1Nafy9v50zXYG6U0+3L2sSZtDX+k/cH2zO2EeYYx/5r5ABSWFcqzAkIgiV/UZ7mHjPmMO15jdBHt+BZr+zHszChl5f5MVsRlsD3lGBUavFwsDGzrz9AOAQxqF4CXi4XjpcdJL0innW87CssKGfnDSHo3782kTpPo0KyDvWsnhN1I4hcNQ8om+HQoeLSEHndAu5HQvDO5hWX8EZ/FyrhMVuzPIKugFAeToneEL0M7BDK0QyAhvq7kl+bzyc5PmLt/LgVlBfRp0YdJnSZxWYvLpBtINDmS+EXDcegPWP6aMZwEGrxC4Pb5xvwFgLVCsy3lGEv2HGXp3qPEZxQA0L65B7dcFsZNMcGUVhQyN24uX+35isyiTOZcPUda/6LJkcQvGp6CDIj7DRKWwXUzjFtFV74BGXuNu4Mih4GjKwezTrB0z1EW7Ehne2oeLbyceeCKNtwUE4xSVv5I+4PBoYMB+HDrh3g4enB95PUykqho9CTxi8Zhxeuw8b9QmGU8MRx5JXQZB+2vQmvNH/FZvLv0AJuTck/5AnByMFOhK3hg6QOsSV+Du8WdG9veyM0dbpYZx0SjJYlfNB7WckheC7vnwd6fjYfDbvjM2Ja5H+3XljXx2fx7adxfXwCDWnNjTAjOFjO7s3bzxe4vWJy0GIXi5X4vM6b1GEqsJZRYS/B09LRv/YSoIZL4ReNUYYXiPHD1NSaemTHI6AIa8hI6sCNr4rN5d2kcsUm5BHo6MXlga27uFYqLo5m0gjS+i/uOEa1G0M63HcuTl/Po8keJ9ImkW0A3+rTow4DgATiaHe1dSyEuiiR+0fiVFcGG6fDHO8Y0lV3GwRV/R3uHsi4hm/eXHWB9Yg7N3By5e0AEt/UJw93J4eTuSflJ/HLwF7Yc3cL2zO0Ulhfi6ejJ3NFzaene0o4VE+LiSOIXTUdRLvzxb+NLwOIKT+wFizMAmw7l8MGyeFbFZeLlYuHOfuGMjW5JWDPXU273LK8oZ+PhjaxOW80zPZ9BKcX07dPRaK6OuJpgj2B71U6IaqvxxK+UGgG8B5iBT7TWr5+2PQz4DPAHcoBbtdaptm1WYKetaLLWesy5PksSv7goeWnGLGPtRhoPhsV+Cp1vBGcvtqUc48Nl8SzdexSAFl7O9IloxmWtm9Enohkhvmc++fv48sdZmrwUgO4B3RndejTDWg2TawKi3qrRxK+UMgNxwJVAKrAJmKC13lOpzFxggdb6C6XUYGCS1vo227YCrXW176WTxC8uWWosfDIEXHyg32PQazI4unIo6wSr47NYn5DN+sTsk1NLBvu4cHlbf8b3DKVzsNfJw6QXpLMgcQE/J/zMofxDjGs3jhcuewGtNeW6HItJJpIR9UdNJ/4+wBSt9XDb8nMAWuvXKpXZDQzXWqcq4zd0ntba07ZNEr+oe+lbjekk45eCeyAMeMp4OtjBuHirtSbuaAHrErJYl5jNyrhMissq6NjSkwm9Qhkb3RIPZ8vJsruzd+Pp6EmoZyjbMrbxyLJHuCriKm5seyMR3hF2rKgQhppO/DcAI7TWd9uWbwN6a60fqlRmFrBBa/2eUuo64HvAT2udrZQqB7YB5cDrWut55/o8SfyiRiWthd9fhbwUeHizMatY5n7wCT/5JQCQX1zG/K1pzNqYwt7D+bhYzIzu2oLxvULpFuJ9yjWB/Tn7mbFjBstSllFeUU5MYAzj2o1jSNgQ+RUg7KamE/+NGK35yom/l9b64UplWgIfAuHAKuB6oKPWOk8p1VJrna6UigCWAUO01gmnfcZkYDJAaGhoj6SkpGpWVYhq0Np4Mtgj0Lgl9M3WUF4KrfpD6yug9WDwawtKobVmR2oe32xKZv62dApLrbQL9GB8rxCu7RaEt+tfXxbZRdnMi5/H3DhjrKDfb/wdJ7MTxeXFODs427HCoimq866e08q7A/u01mfcCqGU+h/GtYDvzvZ50uIXtcpaDnG/QuJySFgOObY2SN9HYNirpxQtKCnnp23pfLMpmR2peTg6mBjVqTnje4XSO9z35K+ACl1BUn4S4V7hVOgKRv84mmJrMRFeEYR7hRPhFUFX/64yfpCoVTWd+B0wLu4OAdIwLu7erLXeXamMH5Cjta5QSk0FrFrrF5VSPkCh1rrEVmYdMLbyheHTSeIXdSo3yfgSCO0D/u0gbQts+dK4IBwYdbLY7vQ8vt2Uwo9b0zheXE64nxs39wrlpp4heLn81b1TXF7M7H2zOZB7gIN5BzmYf5ATZSeY0H4Cf+/9dwDKKsqkS0jUuNq4nXMU8C7G7Zyfaa2nKqVeAWK11j/ZrgO8BmiMrp4Hbcm+LzAdqABMwLta60/P9VmS+IVdbfkSFj0N5cXQagD0vhfajjRmEwOKSq0s2nmY2RuTiU3KxdXRzE0xIUzq14qwZm5nHE5rTUZhBhpNc7fmfLPvG35O+JmPr/wYD0ePuq6daMTkAS4hLkVhjjFb2MZPID8VAjvDfavhtDH+d6Xl8dkfB/l5RzrlFZqhHQK5q3/4Kd1Ap1uevJwnVj5BlG+UJH9RoyTxC1ETrOWwfxEUZkPMJGMKyd9fho7XQsvok8WO5hczc10SX29IIrewjKgWnvRp3Yx2gR60be5BZIA7bpWGi1iWvIwnVzxJlF8U04dOlyGjRY2QxC9EbciMgxmXQ1khhPSG3vcZcwWYzIDRDfTj1jTmxKaw70g+xWUVJ3cN9nGhXaAHo7u25JpuQfye/DtPrXiKjn4d+XzE59LnLy6ZJH4hakvRMdg2Czb9F3ISwbc13DL35Ixhf7JWaFJyCtl/9DgHjh5n/9ECdqYe41B2IW9c34Wbeobwe9LvHD5xmFujbrVTZURjIolfiNpWUQH7FsD2b+CmL4wZw9K3Gs8DOJ55kRegtLyCu77YxJr4LP5zaw+Gd/xrUpi92Xvxd/XHz8WvrmogGhlJ/ELUtfJSeLczWEuNO4FaD4aAKHA6tf++sLScWz7ZwO60fP53Z0/6tvYjvzSfq364ilJrKXd3vpvbom6TB8DEBbuQxG+q7WCEaBIcHOGmLyGkF6x4DT69El4LhjXvG9tLC2H/L7haC/j8jp608nPlni9i2ZF6DE9HT74a9RWXtbiM97e+z5h5Y/jl4C/Ut0aZaDykxS9ETctLg8Pb4chOiLgcQi+D5A3w2TDwCoG7l3KkwpsbPl5LYamVOff2oU2A8ctg4+GNvBn7Jvty9jFz5EyiA6LP82FCGKSrR4j6pqwIDq6GuXeAf1u4YyEH8+HGj9fiaDbx3f19aentAoC1wsqGwxvoG9QXgAWJC+gW0I0g9yA7VkDUd9LVI0R9Y3GBtsOMieIPb4fv7yHc15n/TerF8eJybvt0A4eyTgBgNplPJv0TZSf45/p/MvrH0bwT+w75pfn2rIVoJCTxC1GX2o2AEa/DodWQHU+nIC8+mRjD0fwShr27ineWxFFcZj1Z3M3ixg9jf2Bk+Eg+3/05V/1wFbP2zqKsosyOlRANnXT1CGEPx4+Ax1+3cx7NL2bqwr38tD2dUF9XXh7TkSvaB5yyy57sPbwV+xZbj25l3jXzCPMMq+uoRT0mffxCNARaw/pp0KwNtB0OwNr4LP4xfxcJmScYFhXIi6OjCPZxrbSLJuFYAm182gAwe99shoUNo5lLM7tUQdQf0scvRENgLYUd38LcSUa/P9C3jR+/PDqQZ0a0Y/WBLIa+s5KPlsdTUm50/yilTib9lOMpvLnpTa6dfy2/HfrNbtUQDY8kfiHsxcEJbp5jTAo/a9zJ5O/oYOKBQW1Y+uTlXN7Wnzd/28+Id1ezMi7zlN1DPEKYO3ouQe5BPLXyKZ5c8SQ5xTn2qIloYCTxC2FPHs3hljmgK2DGFbD67ZObgrxdmH5bDF/c2QuAiZ9tZPKXsaTkFJ4s09q7NTNHzeTR7o+yLGUZE3+ZiLXCesbHCFGZ9PELUR8U5cLiF4zJX7qOP2NzSbmVT1Yf5MNl8VRozYNXtGHywAicLeaTZeJy4zhy4ggDgwdirbCSXZxNgGvAGccSjZNc3BWiodswAzL3wdAp4Ox5cnXasSKmLtzDop1HCGvmyls3dqVnK98zdp+zfw5vxb7F/V3v59YOt2Ixy7DPjZ1c3BWioTuRAZs/h2mXQdxfF26DvF2YdksPvrqrN1rDuOnr+PeSOMqtFafs3rdlX3q36M07m9/hhp9vYMPhDXVdA1GPSYtfiPoqNRbmPwSZe6HjdTDyX+D+V9fN8eIyXvppNz9sSaNHmA/vjosmxNf1lEOsTFnJaxtfI60gjTs73cnjPR6v61qIOiJdPUI0FuWlsOZdWPUW3Po9hA84o8j8bWm88OMuAP7v2k6MjT51TJ/i8mI+3/U50QHR9GnZhwpdgUnJj/3GRhK/EI1NQcZfrf3N/zMuAlea9Sslp5DHvt3G5qRcrusWxMtjO+LhXHW//jux75BdnM3zvZ/H1eJaZRnR8EgfvxCNzZ9JvzgPlr4M0/oYvwKs5QCE+Lry7eTLeHRIJPO2pTHi3dUs359R5aFcHFz4OeFnxi0Yx76cfXVVA1GPSOIXoiFx9oIH1hmDvS17Fb4ca4z7AziYTTx+ZVvm3tcXF0czkz7fxGPfbCW7oOSUQ9wffT+fDv+UwrJCbll4C7P3zZZJX5oY6eoRoqHa/g0seBxcfOHhWGPoZ5uScivTlicwbUU87k4OvDg6imuig1BKnSyTW5zLC2teYH36euZdM48QjxB71ELUEOnjF6KpOLoHju6CLjcZy1pDpeS+/8hx/vb9DralHOPytv5MvbbTGYO+xeXG0c63HVprvtzzJcPChtHCvUVd10RcIkn8QjRFu+fBtllw7cfg+tdDXdYKzcx1h3jjt/0AvH59F8Z0bXnG7gfzDnLd/OvQaEaEj2BSx0m0821XV9GLSyQXd4VoiorzIGEZTL/ceAbAxmxS3NEvnMWPDySqhSePzN7Kqwv2UHbaQ1/hXuEsum4RN3e4meXJy7nh5xu4d8m9HDlxpK5rImqZtPiFaEzSNsOcOyA/DQY8AQOfAQfHk5tLyyv456K9/G/tIXqF+/Lhzd0I8HA+4zB5JXnMjZvLLwd/4etRX+Ps4Mx3cd+RV5JH7xa9ae/bHgeTQx1WTJxPjXf1KKVGAO8BZuATrfXrp20PAz4D/IEc4Fatdapt20TgBVvR/9Naf3Guz5LEL8QlKjoGvz4H22fBhG+NO4BO8+PWVJ77YSdeLham3dKdHmFnjvcDxjWAPy8IP7niSRYnLQbA3eJOTGAMg0MHc23ktbVXF1FtNZr4lVJmIA64EkgFNgETtNZ7KpWZCyzQWn+hlBoMTNJa36aU8gVigRhAA5uBHlrr3LN9niR+IWpI6mYI7mG8T9sCLbqC6a/RPPek53PfV5tJP1bEP66O4vY+Yafc9VOVrKIsYo/EsuHIBjYc3kAH3w68PcgYSvpfG/9FO992XNbiMpq7NT/ncUTNq+nE3weYorUeblt+DkBr/VqlMruB4VrrVGX8n5OntfZUSk0ABmmt77WVmw6s0FrPPtvnSeIXooblpcIHPYzEf81/TnniN6+wjMfnbGPZvgyuaOfPY0Pb0jXEu9qHLrGW4GR2Iq8kjzHzxpycCKaDbwemDZ2Gn4tfjVdHVK2mL+4GASmVllNt6yrbDlxve38t4KGUalbNfYUQtckzCMZ8YAzz/HF/4+4fGy9XC5/cHsPfR7VnS/Ixxn60hts/20jsoerN5OVkdjKO4+TFiptW8P2Y73kq5ikO5R/i0eWPUmItOc8RhD1UJ/FX9dvv9J8JTwGXK6W2ApcDaUB5NfdFKTVZKRWrlIrNzMysYhchxEVTyrjP//51ENgJ5k6ElW8a9/wDJpNi8sDWrHl2MH8b0Z5daXnc8PE6bv7vetYlZFf7qV6lFG192jKx40Sm9p9KC7cWVOiK8+8o6lyNdPWcVt4d2Ke1DpauHiHqmbJi+PlRcPeHYf9XZZHC0nJmbUhm+qpEMo+X0D3Um+EdmzOwrT/tm3uc9zrA6awVVsyVri2I2lHTffwOGBd3h2C05DcBN2utd1cq4wfkaK0rlFJTAavW+kXbxd3NQHdb0S0YF3fP+jtSEr8QtUxr42UywZGd4OZvzP17muIyK99uSuHrDUnEHS0AwN/DiQGRfgyM9Kd/pB9+7k7n/KgjJ47w0O8P8Wj3RxkQfOaQ0qLmXEjiP++NuFrrcqXUQ8BvGLdzfqa13q2UegWI1Vr/BAwCXlNKaWAV8KBt3xyl1KsYXxYAr5wr6Qsh6oBSxstaDnMmQnkxjJ8FLaNPKeZsMTOxbysm9m3FkbxiVh3IZPWBLJbvy+CHLWkAdAry5Ip2AQxq5090iA9m06m/BjwdPVFK8cyqZ/j6qq+J8Iqos2qKs5MHuIRoyo7shFnjoTAbRr4O0beC+dztQWuFZnd6HqviMlmxP5MtyblUaPB2tTAg0p8r2vlzeVt/mtl+DRwuOMz4hePxcPTg61Ff4+XkVRc1a3JkrB4hRPUVZBgt/+S10CwS7vwN3JpVe/djhaWsPpDFiv2ZrIzLIKugFDdHM/Mf6k+bAHcAtmZs5c7f7qRnYE+mDZ0mT/3WAhmrRwhRfe4BMGkRjPsKWvX7a4C3zLiTd/6ci7erI6O7tuTtm7qy8e9D+fGBvlgcTDz93XasFcb+3QK68eJlL5Jbkkt+aX5t1kZUg7T4hRBnyk+H96KNh76ueA7CLz/lqd/zmb8tjUe/2cbfR7Vn8sC/Hhgrs5ZhMVc9JaS4NNLiF0JcGjd/GPWG8dTvzGvhrUj4YTJkJ1Rr9zFdW3JlVCBvLY4jPqPg5HqL2UJeSR5T108lrySvtqIX5yGJXwhxJrMFetwBj2yF6z+FNkMhfqmxHuDAUvjj35Cxr8ruIKUUU6/thKuj+ZQuH4DU46l8F/cdUzdMraPKiNNJ4hdCnJ3FGTrfANfNgKfiwTvUWJ+4HJZOgWm94YPu8NvzkLT2lC+BAA9nXh7Tka3Jx/j0j8ST6zv6deT+6Pv55eAvLExcWMcVEiCJXwhRXaZK6WL4VHh8D1z1NviEw4bpsPDJv6Z9zIwDzt7lc1enu+gW0I3/W/9/pBek12UtBJJ87pfYAAAaVUlEQVT4hRAXyysIet4Nt/0AzyQaXUIApSfgv4Mhef1Zu3zMJjP/7P9PNJp/bfyXHSvRNEniF0JcOmdPCIwy3ldYwdHN6P7R+qxdPsEewbxz+Ts8f9nzdgq66ZLEL4SoWc6eMOQfkBYLu74Hzt7l0zeoLwGuAVToCjILZWTeuiKJXwhR87pOgOadYenLUFZ8ssvHxWLmH/N2nTHU8/N/PM/EXyey8fBGOwXctEjiF0LUPJMZhk2FvGTYb9y5E+DhzFPD2rIuMZtfdh05pfiNbW+kxFrCXYvvYvLiyezK2mWPqJsMSfxCiNoRcTlMXgmdrj+5akKvUNo392Dqwr0Ul1lPru8e2J1F1y3i6Zin2ZezjwkLJzA3bq49om4SJPELIWrPn0M9F+UC4GA28dLojqQdK2L6ysRTijqZnbi94+38cv0vPBj9IIOCBwEQnxtPUn5SXUbd6EniF0LUrr0L4J0o4ylfoE/rZlzVuQX/WRlP2rGiM4q7Wdy4r+t9+Lv6A/DGpjcY/eNonljxhHQB1RBJ/EKI2hV6GZgssOQfJ1c9N6o9AP9ctPe8u/9zwD+5u/PdrD+8ngkLJ3Dnb3ey6cim8+4nzk4SvxCidrn5wcAn4cBiSFgGQLCPK/dd3pqFOw6zLiH7nLv7ufjxSPdHWHLDEp6KeYqk/CR2Z+0+5z7i3GRYZiFE7Ssrho96gaM73LcaTGaKy6wMeXslHs4OLHi4Pw7m6rVDy6xlWLUVZwdnFh9aTEFZAddFXlfLFaj/ZFhmIUT9YnGGoVMgcx+kGg07Z4uZ56/qwL4jx5m9KaX6hzJbcHZwRmvNTwk/MWXtFH479FvtxN1ISeIXQtSNjtfCw7EQ2vvkqpGdmtMnohlvL97PscLSCzqcUoo3L3+T6IBonl39LGvT1tZ0xI2WJH4hRN1QCnwjjPfFebZVipfGRJFfVMbbi+Mu+JAuDi58OORDWnu15rEVj7EtY1tNRtxoSeIXQtStpVPg4wFgLQOgfXNPbu/Tipnrk3ho1hYyjhdf0OE8HT35+MqP8XPxY036mloIuPGRxC+EqFuhfeFYEmyffXLV30d14Ikr27J491GGvL2SrzckUVFR/RtP/Fz8+Obqb3ig6wMAZ4wFJE4liV8IUbcir4SW3WDVmydb/Y4OJh4ZEsmvjw2gU0svnv9xFzdOX8f+I8erfVhPR0+UUuzP2c/EXyeSVZRVWzVo8CTxCyHqllIw6Dk4lnxKqx8gwt+dWff05q0bu5KYWcBV76/mjV/3UVRqPcvBzlRqLWVv9l5eWvuStPzPQhK/EKLuRQ4zWv0bZ5wxWbtSiht6BPP7k4MYGx3EtBUJDHl7BfO3pVUrkXf278wj3R9hVeoqFh1cVFs1aNDkAS4hhH1k7gc3f3D1PWexDYnZvLJgD7vT8+ke6s0/ro6iW6jPOfexVli5/ZfbSTmewrxr5uHrfO7PaAzkAS4hRP3n385I+lob0zWeRe+IZvz0UH/euKELKblFXDttLY99s5X0KgZ4+5PZZGZK3ykcLzvOt/u/rY3oG7RqJX6l1Ail1H6lVLxS6tkqtocqpZYrpbYqpXYopUbZ1rdSShUppbbZXh/XdAWEEA3YiSyYPvCMvv7TmU2Km2JCWP7UIB66og2Ldh1h8NsrmLYi/qzdP5E+kcwcOZN7u9xbG5E3aOdN/EopM/ARMBKIAiYopaJOK/YCMEdr3Q0YD0yrtC1Bax1te91XQ3ELIRoD12agTKfc4XMu7k4OPDW8HcuevJxBbQN449f9zN2cetbynfw6YVImsoqyOFF2oiYjb9Cq0+LvBcRrrRO11qXAN8DY08powNP23gtIr7kQhRCNllIw6FnIPQQ75lR7t2AfVz66pTt9Wzfjxfm7OHD07Ld95pfmc93863h387s1EHDjUJ3EHwRUHkEp1bausinArUqpVGAR8HClbeG2LqCVSqkBlxKsEKIRajsCWnSFFa8bXwDVZDYp3h0XjbuTAw/O2nLWWz49HT0ZFTGKb/d/y9aMrTUUdMNWncSvqlh3eqfaBOB/WutgYBQwUyllAg4DobYuoCeAWUopz9P2RSk1WSkVq5SKzczMvLAaCCEaNqVgxL+g+Bisv7DLgAGezvx7XDQHMgqY8tPZx+h/pNsjtHBrwUtrX6LEWnKpETd41Un8qUBIpeVgzuzKuQuYA6C1Xgc4A35a6xKtdbZt/WYgAWh7+gdorWdorWO01jH+/v4XXgshRMMW1gfu+wOGvGgsZydAafX65AdE+vPAoNZ8G5vCvK1pVZZxtbjyYp8XOZh3kOnbp9dU1A1WdRL/JiBSKRWulHLEuHj702llkoEhAEqpDhiJP1Mp5W+7OIxSKgKIBBIRQojT+YSBoytYy2HWOJgxCI7srNaujw9tS89WPjz/404SMwuqLNMvqB9jWo/haOHRJv9E73kTv9a6HHgI+A3Yi3H3zm6l1CtKqTG2Yk8C9yiltgOzgTu0cWYHAjts678D7tNa59RGRYQQjYTZAa56C4rz4b9DYMOZT/eezsFs4v0J3XB0MPHQrK0Ul1Xd3z+lzxSm9p+KUlX1YDcd8uSuEKJ+OpEF8+435uptcyVc/wm4eJ9zl2X7jnLn/2K57bIwXr2m01nLJR5LZHXaaiZ2nFjTUduNPLkrhGj43Pzg5jkw8g2wloKTx3l3Gdw+kHsGhDNzfRJL9xw9a7l58fN4K/YtVqeursmIGwxJ/EKI+ksp6H0v3D4fTGbbr4AHoSDjrLs8Pbw9bQLc+dev+7CeZUz/B7s9SBvvNry49kVyi3NrK/p6SxK/EKL++7NPPjUWds6Fj3rD7h+rLOroYOKxoZEcyChgwY6qnyV1Mjvx+oDXOVZyjFfWvdLkLvZK4hdCNBztRsC9q8CnFcy9A368r8oLv6M6taB9cw/eW3qAcmtF1YfybcfD3R5mafJSFictrt246xlJ/EKIhiWgPdy1BPo/bgzuVsUAbyaT4rGhbUnMOsG8bWcfQWZi1ESe7fUsg0IG1WLA9Y+DvQMQQogLZnaAwS+Cs5cx5EMVhncMpFOQJ+//foCx0S2xmM9s55pNZm7pcAsAJ8pO4Gx2xmwy12ro9YG0+IUQDZPJZLT6/xzT/zRKKZ64si3JOYV8f44RPAFyinO4/qfr+ffmfzeJ/n5J/EKIhi33EMy4HA6tOWPTFe0CiA7x5oNl8ZSUn32yFx8nHwYEDeCLPV/wyc5PajHY+kESvxCiYXMLgKJcWPA4lJeesunPVn/asSLmxJ691a+U4rnez3F1xNW8v/V9vt3XuGftksQvhGjYHF1h1NuQtR/Wvn/G5gGRfvRs5cNHy+LPOpQDgEmZeKXfKwwKHsTUDVNZlrysNqO2K0n8QoiGr+0wiBprzOSVc+o4kEopHr+yLUfyi5m1Ifmch7GYLLw16C1u6XALPQJ71GbEdiWJXwjROIx4HUwWWPvBGZv6tvajT0Qzpq1IOOuELX9yMjvxt15/w8vJi1JrKftz9tdWxHYjiV8I0Th4toSJ841JXarwxLC2ZBWUMHP9oWof8o1Nb3D7L7ezO+vsk7w0RJL4hRCNR1APcHA0hnQuOXUe3p6tfBkQ6cd/ViRwvPj8E7sD3NP5HnycfbhnyT1sPrq5NiK2C0n8QojGpfQETOsDa8680Pv08HbkFpbx39UHq3WoQLdAPhv+Gc2cmzF58WSWJC2p6WjtQhK/EKJxcXSD5p1g8+dQfur8ul2Cvbmqcws+WZ1I5vHqzb3b0r0lM0fOpEOzDkxdP5XCssLaiLpOSeIXQjQ+ve6BE5mw5/RZYuHJYW0pKa/gw2UHqn04b2dv/jvsv3w6/FNcLa5orRv0E76S+IUQjU/EYGjWBjbOOHOTvzvjeoYwa2MyydnVb727OLjQ2rs1AB9s/YAX1rxAWUX1rhXUN5L4hRCNj8kEPe+B1I1wZNcZmx8dEonZpHhnycXdquloduSnhJ94cOmDfB/3Pd/HfU/K8RQAjp44yvdx37Pl6JZLqkJtktE5hRCNU/QEaN4ZAjuesSnQ05lJ/cL5eGUCkwe2Jqql5wUd+r6u9xHoGsgr615h3eF1ALw58E1CPEI4mH+QKeumADAweCBP9niSCO+IS65OTZLJ1oUQTVJeURkD31hO91BvPp/U66KOkV+af/Jir5eTFy4OLpRYS8gtzuWXg78wY8cMisqLuLHtjTzT6xksJktNVuEUMtm6EEIAVFhh0dOwYfoZm7xcLNw/qDXL92eyITH7og7v6ehJc7fmNHdrjouDC2A8+dvcrTmTOk1i4XULuaHtDRwtPHoy6deHxra0+IUQjdsXoyHnIDy63ZiwvZLiMiuD3lxBC29nfri/L+rPuX1rWIWuwKRMpOSncP/v9xPoGnjK9se6P0Zn/86X9BnS4hdCiD/1mgx5KRD36xmbnC1mHhsaydbkYyzZc7TWQjApI9Xml+UT7BFMeUX5KS9N3TbApcUvhGjcrOXwXlfwawO3zz9jc7m1gmHvrsKsFL8+NhCzqXZa/bVNWvxCCPEnswP0vBMSV0DmmbdvOphNPD2sHQcyCpi1Ianu47MDSfxCiMav+0TocQc4OFW5eUSn5gyI9OOfi/ZxKOtE3cZmB5L4hRCNn5sfjH4PfFpVuVkpxZs3dMViVjwxZxvl1oq6ja+OVSvxK6VGKKX2K6XilVLPVrE9VCm1XCm1VSm1Qyk1qtK252z77VdKDa/J4IUQ4oKkxkJC1VMqNvdy5tVrOrEl+RjTVyVWWaaxOG/iV0qZgY+AkUAUMEEpFXVasReAOVrrbsB4YJpt3yjbckdgBDDNdjwhhKh7vzwDi56Biqpb9GO6tuSqLi3495I4dqXl1XFwdac6Lf5eQLzWOlFrXQp8A4w9rYwG/nzm2QtIt70fC3yjtS7RWh8E4m3HE0KIutfrXsg+AIlVt/qVUky9phO+bo48/u22c07O3pBVJ/EHASmVllNt6yqbAtyqlEoFFgEPX8C+QghRNzpeC+6BsP7jsxbxdnXkjRu6cCCjgLcXN775dqF6ib+qm1pPv/l/AvA/rXUwMAqYqZQyVXNflFKTlVKxSqnYzMzMaoQkhBAXwcERYu6C+CWQFX/WYoPaBXDrZaF88sdB1l/kcA71WXUSfyoQUmk5mL+6cv50FzAHQGu9DnAG/Kq5L1rrGVrrGK11jL+/f/WjF0KICxUzyWj1Z527Nf/3UR1o1cyNJ+dsr/YcvQ1FdRL/JiBSKRWulHLEuFh7+rQ2ycAQAKVUB4zEn2krN14p5aSUCgcigY01FbwQQlww9wB4fA+0v+qcxVwdHXj7pq4czitiyk976ii4unHexK+1LgceAn4D9mLcvbNbKfWKUmqMrdiTwD1Kqe3AbOAObdiN8UtgD/Ar8KDWunFeLRFCNBxmB+POnmMp5yzWPdSHhwZH8v2WVObGnrtsQyJj9QghmqYfJkPyenhk6xmjdlZmrdDc9ukGNifl8uMD/S540pa6ImP1CCHE+bQbCceSIO63cxYzmxTvje+Gl4uFB77eTH4j6O+XxC+EaJrajwbPYNjwn/MW9fdw4qNbupOSW8TfvttRLyZTuRSS+IUQTZPZAXrdDQdXwdHd5y3es5Uvz45ozy+7jvDZmkO1H18tksQvhGi6uk8EBxfYNqtaxe8eEM6wqEBeW7SXzUk5tRxc7ZHEL4Roulx9YdIiGDqlWsWVUrx5Y1eCfFx48OutZBWU1Gp4tcXB3gEIIYRdBXW/oOJeLham3dKda6et5bFvtvHsyPZVlgtr5oqHs6UmIqxxkviFEGLnd7D+P3Dnr2A+f7Lu2NKLV8d25G/f7+TqD/6oskwLL2cWPjIAXzfHmo72kkniF0IIRzdIi4Xd86DLjdXa5aaYENoGepBVUHrGtvyiMp77YSfPfLed/94eg1L1ax5fSfxCCBE5DAI7wfL/g6ixxmBu56GUoluoz1m35xWV8cqCPXy5LomJfVvVYLCXTi7uCiGEyQxXvgy5hyD2sxo55KR+rRjcPoCpi/ayJz2/Ro5ZUyTxCyEEQOshEDEIVv4LSi99wnVjHt8ueLlYeHj2FgpLyy/5mDVFEr8QQgAoBcNfg5u+MPr8a0AzdyfeHRdNYtYJXl1Qf0b4lMQvhBB/CoyC8IHG+xoalqFfGz/uu7w1szemsHDH4Ro55qWSxC+EEKdb8hL8/GiNHe6JK9vSNcSbZ3/YQWpuYY0d92LJXT1CCHE6bYUtX0Kve6B550s+nMVs4oPx3bjq/dU8MnsrDw+OrLKcp4sDPcJ8L/nzzkfG4xdCiNMV5cJ70RDUA277ocYO+9P2dB6ZvfWs26NDvJn3YL+LOvaFjMcvLX4hhDidiw8MfBoWPw8Jy6D14Bo57JiuLekS5MWxoqrH9Hd1PPuEMDVJEr8QQlSl1z2wcTr8/gpEXGHc9VMDWvnVzB1Dl0ISvxBCVMXBCa75GNwDayzp1xeS+IUQ4mxaVepvr6gAU+O4EbJx1EIIIWqLtQxm3wwrX7d3JDVGEr8QQpyL2WJ0+6x5D44l2zuaGiGJXwghzmfYq4CCxS/YO5IaIYlfCCHOxysY+j8Oe+Ybk7M3cJL4hRCiOvo9Al6hsPINe0dyyeSuHiGEqA6LC4ybCd6h9o7kkkniF0KI6moZbfxbYYXyEnB0tW88F0m6eoQQ4kJYy+CTIbDkRXtHctGqlfiVUiOUUvuVUvFKqWer2P5vpdQ22ytOKXWs0jZrpW0/1WTwQghR58wWCIqB2E/h6G57R3NRzpv4lVJm4CNgJBAFTFBKRVUuo7V+XGsdrbWOBj4AKg9nV/TnNq31mBqMXQgh7OOKv4OzF/zytxqbsKUuVafF3wuI11onaq1LgW+AsecoPwGYXRPBCSFEveTqC1c8D4dWw8veYLXNp7vgCZjiZby+GFNvvxSqc3E3CEiptJwK9K6qoFIqDAgHllVa7ayUigXKgde11vMuMlYhhKg/Yu4EkwMcPwLK1oZuOxzc/CEnAXbOhaQ10Kq/feOsQnUSf1XD0p3ta2w88J3W2lppXajWOl0pFQEsU0rt1FonnPIBSk0GJgOEhjb8W6WEEE2AyQwxk05d13a48SorgoTlsOb9epn4q9PVkwqEVFoOBtLPUnY8p3XzaK3Tbf8mAiuAbqfvpLWeobWO0VrH+Pv7VyMkIYSoxywuMPJfcNl99o6kStVJ/JuASKVUuFLKESO5n3F3jlKqHeADrKu0zkcp5WR77wf0A/bUROBCCFGvdb6hxmbuqmnnTfxa63LgIeA3YC8wR2u9Wyn1ilKq8l06E4Bv9KmT+HYAYpVS24HlGH38kviFEE1DQYYxsFv+2TpJ7KNaT+5qrRcBi05b9+Jpy1Oq2G8tcOlT1AshRENUegLWfWRc/L3yFXtHc5I8uSuEELXFNxw6jIHYz6E4397RnCSJXwghalO/R6AkH7Z8ae9ITpLEL4QQtSmoB4T1h/X/Mcb5qQck8QshRG3r/ziE9YGS4/aOBJBhmYUQovZFDjVe9YS0+IUQoq4c2QWHd9g7Ckn8QghRJ6zl8PWN9WLCdkn8QghRF8wO0PteOLgSPuwFH/WG7+/5a/usccasXnVA+viFEKKu9LwLsg/8dZHXJ+yvbb4RVD0mZs2TxC+EEHXFyQPGflT1thGv1VkY0tUjhBBNjCR+IYRoYiTxCyFEEyOJXwghmhhJ/EII0cRI4hdCiCZGEr8QQjQxkviFEKKJUadOkWt/SqlMIOk8xfyArDoIpz6TcyDnoKnXH+QcwF/nIExr7V+dHepd4q8OpVSs1jrG3nHYk5wDOQdNvf4g5wAu7hxIV48QQjQxkviFEKKJaaiJf4a9A6gH5BzIOWjq9Qc5B3AR56BB9vELIYS4eA21xS+EEOIiNbjEr5QaoZTar5SKV0o9a+946oJS6jOlVIZSaleldb5KqSVKqQO2f33sGWNtUkqFKKWWK6X2KqV2K6Ueta1vSufAWSm1USm13XYOXratD1dKbbCdg2+VUo72jrU2KaXMSqmtSqkFtuWmVv9DSqmdSqltSqlY27oL/jtoUIlfKWUGPgJGAlHABKVUlH2jqhP/A0actu5Z4HetdSTwu225sSoHntRadwAuAx60/XdvSuegBBiste4KRAMjlFKXAf8C/m07B7nAXXaMsS48CuyttNzU6g9whdY6utItnBf8d9CgEj/QC4jXWidqrUuBb4Cxdo6p1mmtVwE5p60eC3xhe/8FcE2dBlWHtNaHtdZbbO+PY/zhB9G0zoHWWhfYFi22lwYGA9/Z1jfqc6CUCgauAj6xLSuaUP3P4YL/Dhpa4g8CUiotp9rWNUWBWuvDYCRGIMDO8dQJpVQroBuwgSZ2DmzdHNuADGAJkAAc01qX24o09r+Hd4FngArbcjOaVv3B+LJfrJTarJSabFt3wX8HDW3O3apmIpbbkpoIpZQ78D3wmNY632jwNR1aaysQrZTyBn4EOlRVrG6jqhtKqauBDK31ZqXUoD9XV1G0Uda/kn5a63SlVACwRCm172IO0tBa/KlASKXlYCDdTrHY21GlVAsA278Zdo6nVimlLBhJ/2ut9Q+21U3qHPxJa30MWIFxvcNbKfVnA64x/z30A8YopQ5hdPEOxvgF0FTqD4DWOt32bwbGl38vLuLvoKEl/k1ApO1KviMwHvjJzjHZy0/ARNv7icB8O8ZSq2x9uZ8Ce7XW71Ta1JTOgb+tpY9SygUYinGtYzlwg61Yoz0HWuvntNbBWutWGH/3y7TWt9BE6g+glHJTSnn8+R4YBuziIv4OGtwDXEqpURjf9GbgM631VDuHVOuUUrOBQRij8B0FXgLmAXOAUCAZuFFrffoF4EZBKdUfWA3s5K/+3b9j9PM3lXPQBePCnRmjwTZHa/2KUioCowXsC2wFbtVal9gv0tpn6+p5Smt9dVOqv62uP9oWHYBZWuupSqlmXODfQYNL/EIIIS5NQ+vqEUIIcYkk8QshRBMjiV8IIZoYSfxCCNHESOIXQogmRhK/EEI0MZL4hRCiiZHEL4QQTcz/A8cy0ueTMiDvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  Kaplan-Meier Estimate, create the life table\n",
    "\n",
    "production=pd.read_excel(\"./data/production.xlsx\",\n",
    "                         parse_dates=['Retail date','Production date'],\n",
    "                         date_parser=lambda x: pd.datetime.strptime(str(x), \"%Y-%m-%d %H:%M:%S\"))\n",
    "failure=pd.read_excel(\"./data/failure.xlsx\",\n",
    "                     parse_dates=['Repair date'],\n",
    "                     date_parser=lambda x:pd.datetime.strptime(str(x),\"%Y-%m-%d %H:%M:%S\"))\n",
    "failure['Censored']=0\n",
    "print(production.head(3))\n",
    "print(failure.head(3))\n",
    "sample_list=production.merge(failure,on='ID',how='left')\n",
    "sample_list['Repair date'].fillna(value=np.max(sample_list['Repair date']),inplace=True)\n",
    "sample_list['Censored'].fillna(value=1,inplace=True)\n",
    "sample_list['Operating_days']=list(map(lambda x: x.days,sample_list['Repair date']-sample_list['Retail date']))\n",
    "sample_list['Operating_time']=sample_list['Operating_days']//30.5\n",
    "print(sample_list.head(5))\n",
    "\n",
    "life_table=sample_list.loc[:,['Operating_time','Censored']].groupby(['Operating_time']).agg({\"Censored\":[\"size\",\"sum\"]}).reset_index()\n",
    "life_table.columns = ['_'.join(col) for col in life_table.columns.values]\n",
    "\n",
    "life_table['n_event']=life_table['Censored_size']-life_table['Censored_sum']\n",
    "life_table=life_table.loc[life_table['Operating_time_']>0,:]\n",
    "life_table.sort_values(by=['Operating_time_'],ascending=False,inplace=True)\n",
    "life_table['n_risk']=life_table['Censored_size'].cumsum()\n",
    "life_table.sort_values(by=['Operating_time_'],ascending=True,inplace=True)\n",
    "life_table.rename(columns={'Operating_time_':'time','Censored_sum':'n_censored'},inplace=True)\n",
    "life_table.drop(columns=['Censored_size'],inplace=True)\n",
    "\n",
    "life_table['hazard_rate']=life_table['n_event']/life_table['n_risk']\n",
    "life_table['survival_rate']=1-life_table['hazard_rate']\n",
    "life_table['survival']=life_table['survival_rate'].cumprod()\n",
    "life_table['failure_probability']=1-life_table['survival']\n",
    "\n",
    "# standard deviation according to Greenwood formula\n",
    "d_i=life_table['n_event']\n",
    "n_i=life_table['n_risk']\n",
    "life_table['std_error']=life_table['survival']*np.sqrt(np.cumsum(d_i/(n_i*(n_i-d_i))))\n",
    "life_table['lower_95_ci']=life_table['survival']-stats.norm.ppf(1-.05/2)*life_table['std_error']\n",
    "life_table['upper_95_ci']=life_table['survival']+stats.norm.ppf(1-.05/2)*life_table['std_error']\n",
    "print(life_table.head(5))\n",
    "plt.plot(life_table['time'],life_table['survival'])\n",
    "plt.plot(life_table['time'],life_table['lower_95_ci'], linestyle='dashed')\n",
    "plt.plot(life_table['time'],life_table['upper_95_ci'], linestyle='dashed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# weibull regression\n",
    "- Parametric regression model for survival data\n",
    "- univariate regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weibull model train finished with score: 97.11571195761121%\n",
      "Weight: 1.551033921803207, Bias: -7.779230574903931\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAE4CAYAAAB47tUIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlA1HX++PHnzMAAw40IiAioeN+ipQJmpna5HmXp12y3ds12t9qtvm25323L1m9b3992bWatqbXVtpZlUml50IYHCHngfSAqh+KogDicMwMzvz9GRkbOgRmYkdfjL/nM5zOfF59BXrzPl8JsNpsRQgghHETZ2QEIIYS4sUhiEUII4VCSWIQQQjiUJBYhhBAOJYlFCCGEQ0liEUII4VCSWIQQQjiUJBYhhBAOJYlFCCGEQ3k09cKqVaucckONRsP8+fOd8t5CCCE6X5OJZevWrU65YVBQkCQWIYS4gTWZWADUajXx8fEOu9muXbsc9l5CCCFcU7OJRaPR8OSTTzrsZpJYhBDixieD90IIIRyqyRbLbbfdhq+vr0Nv5oz3FEII4VoUUo9FCCGEI0lXmBBCCIeSxCKEEMKhmp0V1lpFRUXs2bOHixcvAhAWFkZ8fDzdu3d3xNsLIYRwI+0eY1m7di3r16/HZDLZHFcqlcyaNYu5c+e2K0AhhBDupV0tli1btrBu3To0Gg0JCQlERERgMBg4efIkWVlZfPXVVwQHBzNt2jRHxSuEEMLFtSuxbN68maCgIF555RVCQkJsXtu5cyfLli1j8+bNkliEEKILaXbwfsuWLc1erNVqGTFiRIOkApCYmIharUar1bYvQiGEEG6l2cSyevVq/vznP3P27NlGX/fz86OgoKDR1woLCzEYDPj5+bU/SiGEEG5DtWTJkiVNvajX68nMzCQlJQWj0ciAAQNQqVTW1y9dusSePXs4c+YMgYGBKBQKSktLycrK4v3336eiooKJEycyevTojvhehBBCuIAWZ4Xl5+ezYsUKcnJyiIiIYOHChQwbNgyAyspK/vKXv3DmzJlGr42OjmbJkiWyjYsQQnQhrZpubDab2bx5M5999hlVVVUkJSXxi1/8An9/f2pqati6dSu7d+/mwoULgGUdy9ixY5k6dSqenp5O/yaEEEK4DrvWsVy+fJnVq1eze/du/Pz8ePDBB5k0aZITwxNCCOFu2rRAcs+ePaxevZqSkhKGDBnCokWLiIiIcEZ8Qggh3EybV95XV1fz2WefsXnzZlQqFffccw8zZ860GdwXQgjR9bQ6sZhMJiorK9FoNCiV12Ypnz59mhUrVpCbm0tUVBSLFi1iwIABTgtYCCGEa2vVrLBPP/2Uw4cPU1NTg4eHB0OHDuW//uu/iI2NBSxJ57vvvmPt2rUYDAYmT57MggUL0Gg0bQrqp59+Yu3atRQWFhIcHMydd97J9OnTbc4xm82sX7+erVu3otPpiIuL4+GHH7bGJIQQonM0u44lPz/fukDSZDLh4+ODwWBAq9WyY8cORo8eTVBQEAqFgv79+5OUlERhYSHp6els27aNbt260atXL7sCOn78OH/961+Jj49n3rx5BAQE8Pnnn+Pj40P//v2t5yUnJ7Nu3Truv/9+7rzzTs6cOUNycjK33HIL3t7ebX4gQggh2qfZxLJixQoKCgoYO3YsL730EnPnzuX222/nwoUL5OfnU1JSQkJCgvV8jUZDYmIiUVFR7N+/n23btpGTk0NSUlKrA1qxYgXBwcH84Q9/ICIigsGDB1NRUcF3333H9OnTUSqVGAwG/va3vzFz5kxmzJhBeHg4Y8eO5bvvvsNsNjN06NB2PRQhhBBt1+yWLseOHUOpVPL4448TEBAAQEBAAL/97W9RKBQcO3as0evGjx/Pm2++yZQpUzhw4IBdAeXm5loXYNYZMWIEFRUVZGdnA5CdnU1VVRXjx4+3nuPt7U18fDxZWVl23U8IIYRjNZtYDAYDHh4eDbqWvL298fT0xGAwNHmtRqPhkUce4S9/+YtdAdXds766r+v2LDt37hxKpZIePXrYnBcVFUVhYaFd9xNCCOFYzSaWqKgoDAYDP/74o83xH3/8EYPBQFRUVIs3qD8u0hoRERHk5OTYHKv7ury8HICKigq8vb1tZqcB+Pr6otfrqampafC+KSkpLF68mMWLF9sVjxBCCPs0W4/l7rvvZvny5fzjH/8gNTWV8PBwLly4wPHjxwG48847HR7Q1KlTWbVqFSkpKYwbN46cnBw2bNgAYJNIFApFg2ubm+A2ZcoUpkyZYv1aWjaiLUJDQykqKursMNyGuz2vzoy3I+7d3ntERka26rxmE8vEiRMpKSnhyy+/5Pjx49aE4uHhwZw5c5yyncvkyZPJy8tj1apVvP/++3h5efHAAw/wwQcfEBQUBFhaJlVVVZhMJptkU1lZiZeXV4OuNCGEEB2nxd/As2bNYsqUKRw/fhydToe/vz+DBg1yWp0VpVLJr371K+bOnUtJSQlhYWGcO3cOgH79+gHQs2dPTCYTWq3WJoOeO3eu1RlVCCGEc7TqT3s/Pz/GjBnj7Fga3LMueW3evJkBAwbQs2dPwDJu4+Pjw65du7j33nsBS+2YvXv32nR3CSGE6Hgu12eUnZ3N8ePHiY2NpaqqirS0NA4cOGAzu0ytVjNr1izWrVuHr68vPXv2ZMOGDZjNZu64445OjF4IIYTLJRYPDw927drFF198gVKpZODAgSxdupTo6Gib82bNmoXZbCY5OZmysjL69u3L888/bx2HEUII0Tma3Cvs+eefx9/fn+eee85hN3PGe7aVzAoTbeFus5w6m7s9L5kV1rx2zwo7efKkw//6d8Z7CiGEcC3NLpAUQggh7NXsGItOp+Opp57qqFiEEELcAJpNLCaTqcPHIpYsWcLRo0cbfe1///d/6d+/v9RiEUIIF9ZkYlm4cKFTbujl5dXs6wsXLqSystLm2Nq1azlz5gx9+/YFrtViefDBB4mMjGTjxo0sXbqU119/XcZwhBCikzWZWKZOndqRcVhdv7FlTU0Np06dYsKECahUKgwGA8nJycyePdu6ZqV///489thjbNq0iXnz5nVG2EIIIa5y+cH7/fv3U1FRYS0oJrVYhBDCtbl8YklLSyMkJIRBgwYBUotFCCFcncutvK+v/v5fddvkt6YWS2O7G6ekpJCSkgLAq6++SmhoqPO/AXHD8fDwkJ8dO7jb8+rMeDvi3h31/bl0Ytm7dy/V1dXWbrA69tZigYb1WNxpNbBwHe62kryzudvzkpX3zWvtynuX7gpLS0sjIiLCOhsMbGux1Ce1WIQQwjW4bGKprKxk//79DVor9Wux1Ce1WIQQwjW4bGL56aefMBqNDRJL/VosderGYkaNGtXRYQohhLiOy/YbpaWlERMT02Bdi9RiEUII1+aSiUWn03H48GHmzp3b6OtSi0UIIVxXk/VYGrNq1SqmTJlyQ+zJJWteRFu42yynzuZuz0tmhTWv3fVYGrN161a2bt1KXFwct912GwkJCS3u/WWv1NRU3n333QbHFy5cyLRp0wBkE0ohhHBhdiWWPn36cPr0aXJycsjJyeHjjz8mMTHRKa2YF154AbVabf06PDzc+m/ZhFIIIVyXXYnllVdeITc3l5SUFNLS0qisrLRpxUyZMoWEhASbhNBWcXFxeHt7Nzgum1AKIYRrs3u6cWxsLAsXLmTFihX85je/oX///gDk5OTwj3/8g0WLFrF69Wpyc3MdHSsgm1AKIYSra/OsMLVazaRJk5g0aRJnz54lJSWFHTt2UF5ezpYtW9iyZQt9+/Zl6tSpbWrFPPHEE5SVlREeHs706dOt2/g3twll/bUtQgghOodDphtHRUXx0EMPsWDBAjIyMti4cSOnT5/m1KlTnDp1io8++oiJEydy1113ERER0ex7BQUFMXfuXOLi4jCZTKSlpbFy5Ur0ej3Tp09v8yaUQgghOoZDfwOfPn2agwcPcvbsWZvjVVVVbN68ma1btzJ9+nTmz5/f6EaSACNHjmTkyJHWr0eNGoXRaOSrr77irrvuAtq2CaXsbiwcwd126+1s7va8ZHdjB92nvW9QUVHBtm3b+OGHH2wSSq9evZg6dSrjx4/nwIEDbNq0iZycHL755hv8/PyYOXNmq+8xbtw4du3axaVLl2w2oazfamlpE0rZ3Vg4gruty+hs7va8ZB1L85yyjqW+Y8eOkZKSQmZmJkaj0fJmHh6MGzeOqVOnMnDgQOu5SUlJJCUl8d133/HRRx/xn//8x67EUkehUNhsQln/m5RNKIUQwjXYlVjKy8tJTU3lhx9+sFm5Hh4ezpQpU7j11lvx9/dv8vq77rqLL774gkuXLtkVZGZmJv7+/oSGhhIUFGTdhPLee+8FbAuCCSGE6Fx2JZZHH32UmpoaAJRKJfHx8UybNo3hw4e3+j00Gg2VlZVNvv7aa68RFxdHTEwMJpOJ9PR00tPTefjhh1EqlbIJpRBCuDi7EktNTQ0hISHcdttt3HbbbQQHB9t9w8ceewyDwdDk65GRkfz4448UFxdjNpuJiori8ccfZ+LEidZzZBNKIYRwXXZtQrl7927i4+MbTPV1R7IJpWgLdxuM7mzu9rxk8L55TilN3L9/fyoqKlp9fllZGVeuXLHnFkIIIdycXYll0aJFPPPMM60+f/HixTz66KN2ByWEEMJ9Ob1Py46eNiGEEDcAp+59UlNTY/d4TEZGBhs2bKCwsBC9Xk9oaCgTJ05k5syZ1sWPUo9FCCFcl9MSS2lpKVeuXCEgIMCu68rKyhgyZAgzZsxAo9GQk5PDF198QWlpKb/61a8AqccihBCurNnEcuLECY4fP25zrLq6mq+//rrJa8xmM5WVlezduxez2UxcXJxdAdXtYlxn6NCh1r3GfvnLX2I0GqUeixBCuLBmE8vBgwf58ssvbY5VV1fz73//u3Vv7uHRpq1brufv729dmNlSPRZJLEII0bmaTSzdunWzaXHk5OTg4eHR7FiGUqnEx8eHXr16ceuttxIVFdWmwEwmE0ajkTNnzvD9998zbdo0FAqF1GMRQggX12ximTx5MpMnT7Z+PXfuXPz8/Hj55ZedHtiDDz5o3dxy4sSJLFiwAEDqsQghhIuz6zfwwoUL8fLyclYsNpYuXYrBYCAnJ4cvv/ySDz74gIULFwJSj0V0HnerL9LZ3O15ST0WB93HnpOvH1h3pj59+gAwcOBA/P39Wb58OdOnT5d6LKJTudsWJZ3N3Z6XbOnSPKds6dJZevfuDcDFixdt6rHUJ/VYhBDCNTTZYlm1ahVgqUE/Z84cm2P2quvCaqsTJ04AEBYWRkhIiNRjEUIIF9ZkYtm6dStgafrUJZa6Y/ayJ7G8/PLLDBs2jF69eqFUKjlx4gTffvstEyZMICIiAkDqsQghhAtrMrHUrRMJCQlpcMyZ+vbty7Zt27h48SIqlYrw8HDmz59vM74j9ViEEMKWKSMV8/pPoOQSKJVgMkFIdxSzH0Q5blKHxmJXPZYbidRjEW3hboPRnc3dnpe7DN5fSyJF4OsHRgMY9I2frPZC8eBjKMdNksF7IYQQDZkyUjF/stzSMsEMFWVNJxUAg96ShDqQrCQUQgg3YcpIxfzhW5ZuLnuUdGwrTFosQgjhBqwtFXuTCkBIxy76bHG6sSO0d7qxEEJ0VTaD8m2h9kIx+0HHBtWCFqcbO4IkFiGEsJ+1ldLcGEpjFAowmzttVliL042FEEJ0PLvGU3z9oaIcQkI7JZFcr8nE8uSTT3ZkHEIIIa5q9XhKvanErkRmhQkhhIsxr/+k5e6vTurmag1JLEII4SJaNVDvoq2U+iSxCCGEC6jatrnlgXql0uWTCrjJ7sZCCHGjqmul6FqaTuwGLZU6Lre7sRBCdAWmjFTMn620bMnSEhceT2mMy+1uLIQQN6o2LXYM6Y7q/1Y7LygnsGu6sUxBFkKItmnTYsdOWDXvCDJ4L4QQHaBVU4jrc7Pur/oksQghhBPZ3f3lRoP0TWlzYjGZTBw5coRTp06h0+kACAgIoG/fvgwZMgSlUjZOFkJ0bXZ3f/n6o5j3iFsnFWhjYvnhhx9Yu3YtpaWljb4eFBTEfffdx5QpU9oUlFar5ZtvvuHkyZPk5+czaNAglixZYnOO2Wxm/fr1bN26FZ1OR1xcHA8//DCxsbFtuqcQQjhaq7u/QroT8PPfUjEk3vlBdQC7E8uqVatsph37+/tbZ45dvnwZnU5HaWkpK1eu5MyZMzzyyCN2B1VQUEBWVhb9+vWjpqam0XOSk5NZt24dDz74IJGRkWzcuJGlS5fy+uuvExQUZPc9hRDCUVrV/XXdGIpPaCgVblTGuTl2JZZdu3ZZk8rNN9/MnDlziI6OtjmnoKCAdevWsWvXLlJSUhg6dKjd05Tj4+MZO3YsAK+//jplZbbzvA0GA8nJycyePZs77rgDgP79+/PYY4+xadMm5s2bZ9f9hBCiLRrUnodWr0txtynE9rArsWzatAmA2267jUWLFjV6Tq9evXjyySfx9fUlJSWFTZs22Z1YWhqfyc7OpqqqyuZ9vb29iY+PJysrSxKLEMIpGiSS6iqovdqr0pqEAm47hdgedo2w5+bmolAoWvWLe+7cudZrHO3cuXMolUp69OhhczwqKorCwkKH308IIawD8SWXALMlkdQ23lXfpJDubj/jqzXsHmPx9fUlICCgxfMCAgLw9fWltra2TYE1p6KiAm9v7wYtG19fX/R6PTU1NXh42H5rKSkppKSkAPDqq68SGtqxNaDFjcHDw0N+duzgbs+rsXirtm2m/NN/YL50oV3vreweTvf319t1b0frqM/DrsTSo0cP8vLy0Ov1eHl5NXtudXU1VVVVTpulpVAoGhwzm81Nnj9lyhSbWWpFN8ggmehYoaGh8rNjB3d7XtfH2+bSwNdTe2Ge8UCzz6IjnlV77xEZGdmq8+zqCps0aRImk4ktW7a0eO7WrVsxmUzccsst9tyiVXx9famqqsJ0XXW1yspKvLy8GrRWhBCitUwZqdQ+9ytqH5lpKQ3c3qTSRbq/6rPrN/Dtt9/O4cOH+fe//01tbS133XUXarXa5hyj0ch3333HZ599xtixY62zthypZ8+emEwmtFqtTQY9d+5cqzOqEEJcr0FNFJOZfSEDOBrYmwVnNjW8QKUCb42l3rx1Vpjr1J7vLC3WY7leYGAg3t7erFmzhuTkZPr160dISAgKhYLi4mJOnjxJVVUVGo2GwMBAVq1a5fBt8/v374+Pjw+7du3i3nvvBUCv17N37942L8oUQnRd9WuiVKq8qFb7E2KwzPK64ulHcvQtzDi7nQBj5bWL3HgvL2drsR5Lc6qqqjh48GCjr1VWVloHy+1NLHq9nqysLABKSkqoqqoiIyMDgFGjRuHl5cWsWbNYt24dvr6+9OzZkw0bNmA2m53SQhJC3LhMGalUf/o+e/17s2PIHezrNpDJ5/fw6EnLQPtNRUfgOHjVGi0X3AB7eTlbi/VYOsOVK1d44403bI7Vff3OO+8QFhbGrFmzMJvNJCcnU1ZWRt++fXn++edl1b0QokWmjFQMyZ9ygBB2ho3kpzHPUe1hmZCkMJu4ova1nutbW82tl7LAZJZWSispzM1NpbqByXoX0RbuNsups7nK82psYeO/oqfwVcxk6zn9dPkkXtzPhIsH6WbQXbu4g1ooN9KsMJk+JYS4odXsSuXE19+wM2Q8/TwLuOWCpZt9XNFh9nYbROLF/SRcPEBEdUnDi6WF0iaSWIQQNxyz2cyZy3p25OnYftyHomGWLagGl562Jpa+ZWd5c8+bjb+BjKO0iyQWIcQNJfXMFdYeLuaczmA5oA6kW3UpiRcPkHhxv/W8BkuslUowmyG4a08VdoQ2JZaysjK2bdvG8ePHKS4uRq/XN7nqXaFQNBiIF0IIR7lQbkCBgjA/TwBqTGbO6QwE1FSQcOEAiZcOMKA0FyXNDCdfbaGETZ/jEmNC7s7uxLJv3z6WLVtGZWVlyyc7UWpqKu+++26D4wsXLmTatGmdEJEQoqOUVNWQlqdjR56OE0XV3BFQwaKd70JJETcFdOMFz1CGlWSjMpsaf4P6Cxu7+GJGZ7ArsRQWFvLGG29gNBoZOnQoo0eP5uOPP0aj0TBv3jxKS0s5dOgQJ0+exN/fn3vuuQdPT09nxQ7ACy+8YLP6Pzw83Kn3E0J0Dp2+loyCMrbn6jh8odLa/vBSmFCcOGQtquWrK2IkjbQ6lMqrU4YlkTibXYllw4YNGI1Gxo8fz5NPPgnAxx9/jFqt5vbbbwcs2+X/9NNPvP3222RkZDQoKexocXFxeHt7O/UeQojOt/ZwEd8evwyAB2ZGX8kh8exPjLl8HO+aVuznZTKjWvm1k6MUYGdiOXLkCAD33HNPs+fddNNNLFiwgA8//JDvvvuO6dOntz1CIUSXoq8xsaewnB25ZYzt6cutl/ZjXv8JiUZvCvpPJ9Gvmpt/+grfqiv2vXGI+2zf7+7sSiwlJSWoVCp69eplc9xoNDY4d9KkSXz00UekpaU5NbE88cQTlJWVER4ezvTp05k6darT7iWEcA5jrZkD2gp25OrIOFtOdY1lbKTsUhGT/mPZFLIf8MLe99p2gy5QtdGV2JVYlEolHh4eNrVQvL29qaysxGQy2RTe8vb2xtvbG61W67ho6wkKCmLu3LnExcVhMplIS0tj5cqV6PV6aSEJ4Ua+PV7C54eKKDNcG2iPqzxPYuFeJhQdbNu29TI436nsSiwhISEUFhZiMBisA+ahoaGcPXuWvLw8evfubT23srKSyspKpw3ejxw5kpEjR1q/HjVqFEajka+++oq77rqrQXVJqSApHMHdKiJ2tuufl9ls5oi2jCAfT6KCfADoFlRDmeEivbtpmORZypjvl9Pjynn7b3Z1HYoyNAy/B36Nzy23tzvejtRlK0hGRkZSWFiIVqslOjoagAEDBnD27Fk2btzI448/bj137dq1AA3q0jvTuHHj2LVrF5cuXWowO0wqSApHcJW9r9xFaGgoly5dsq6C35mn42JFDXcPCGbRmHBMGakM/fpz3qyqJcarFvTVllry9rpupXwFUNGGz6kzP98uu1fY6NGj2bNnDz/99JM1sUydOpX//Oc/7Nixg4KCAmJjY8nPz+f06dMATqkg2ZLGyhYLITrWOZ2Bb3Ly2XxMy9m6VfBANx8Pgr1V1rK/vgY9vmDJBq2h9oLxk+HQHsumktLV5XLsSizx8fGMHz+e2tpa67HevXvzwAMP8Omnn5Kbm0tubq71tbFjx3LXXXc5LNiWZGZm4u/vL10VQnQSs9ls/cNuY/ZlNp6wTA8O8FIxwauMxH3JDDx7EGVmN8z66taPn8gaFLdiV2IJCgqyrl+p72c/+xnDhw9n165dFBUVodFoGDlyJKNHj3ZYoNd77bXXiIuLIyYmBpPJRHp6Ounp6Tz88MMNxleEEM5TWlVDWr5l4eKUvoFMjbPURJoUG4BZ6clNPbwYmrsb5b/qlfy9upixVWRDSLfjsE0oY2JiiImJcdTbtSgyMpIff/yR4uJizGYzUVFRPP7440ycOLHDYhCiqyrX17KroIzteZZV8Kary+A1nkprYukf6sOEgb0oKiqi9q1PWt868fUHL2/p5nJjbru78fz585k/f35nhyFEl/Ov/ZdYf6yYq0tN8FBCfKQfE2MDGNvTD7hWWOvC5SIIDm19C0XthWLeI5JI3Fy7EotWq+X06dNcuWJZARsYGEifPn2IiIhwSHBCiM5lqDWxt7CCXgFqogItpXu7aTwwmWFEhIaJsQGMi/LHz0tlvaZuUL5V3V7SOrkhtSmxHDx4kDVr1lhnfl2vT58+zJs3jxEjRrQrOCFEx6sxmTmorWB7ro6MgnKqakz8bGAwC+MtU/hv6R3A+F7+BPl4WJLI+59QWy8xmNe3sttLWic3LLsTS3JyMmvWrLE5ptFoAKxb6Z8+fZq//vWvzJs3j9mzZzsgTCGEsx2/VMV/Tl8hvaCMMv21mZ99Q7yIvtpaAdB4qtB4Nt4ysfm6MSHdpXXSBdiVWA4dOmRNKv369WP27NkMGTLEurtwdXU1R44cITk5mezsbD777DPi4uIYNmyY4yMXQrSL2WzGZAaV0jI9OPXMFTbnlAIQFaBmYmwAiTEB9AxQY8pIpXblJzZJodGWiUF/dWpwI3VQQrqj+r/Vzv62hAuwK7F8++23gGV9ytNPP91gWq+3tzfx8fGMGjWKN954g927d7NhwwZJLEK4CLPZTF6pnh15ZezI03H/0G5M6WuZxXVb30B81SqSYvyJCfKyrkexu2ViMlkWMdZ/XTaB7FLsSiw5OTkAPPTQQ82uFVEqlTz00EPs3r2b7Ozs9kUohGi3Qp2BHXk6tufqbFbBZ52vsCaWft186NfNp8G1bWmZWFs0V2eFSbdX12JXYjEajfj6+rZqZXtoaCi+vr6NbqkvhOg472ZqrV1cAP5eKib08icp1p/B3TU259ZNE67f5UVJE3tLNdMyUY6bBOMmyd5qXZRdiSUsLAytVktNTQ0eHs1fWlNTQ3V1dYduQilEV1e3Cn5YuIboIMuAe2ywFz4eSsb18iMpJoARPXzxUDbcT6/JLi9fv8Y3hqzfMpEBeVGPXYllwoQJrF27lp07dzJp0qRmz925cye1tbUkJCS0Jz4hRAvK9bVknC1jR66Og1dXwc8aFMLDo8MAuK1PILf1CcTLw9J9bcpIpbaRZNBkl5enusWWiRD12ZVYZs6cyb59+1i9ejVqtZoJEyY0el56ejqrV6+mX79+zJgxwyGBCiFsZZ4tI+XUFfYVll+3Ct6XwWHXxkrqEgo03SoxQdNdXhXlKH71lLRMRKs1mVi+/vrrRo+PHDmSwsJC/v73v/P5558zZMgQQkJCUCgUFBcXc/ToUc6fP2/diHLjxo3MnDnTad+AEF2FodaEAvBUWRLFvsIKfjpbjlIBwyM0JMVYFi76e1m2pK9d1vpWiXn9J5aa8I2tkg8JlZaJsEuTieXf//53ixdrtdomSw9XVlbyxRdfAEhiEaKN6lbB78izrIJ/ZEw4k/sEAjAtLohegWoSogMI9rn2X7lNrZKSIkur5PppxDJNWLRBk4klLi6uI+MlDndbAAAgAElEQVQQQlxlMps5drGKHXk60vLL0NVbBZ9TXGVNLL2zM4i92j1V66BWiQmky0u0W5OJ5eWXX+7IOIQQV726/RyZZ8utX/cMUDMxJoDEWH+iAiwzvZzVKpEuL+EIbrttvhA3grxSPdtzdUyMDSDm6vTg4REaTmuvkFi4h8S8XcR6GVF2fxBlwCTrddIqEa5MEosQHex8mYEduTp25OnIv2JZBW8ym/nFKMv04KklB7lj+3IUdYmjAmtrxJoApFUiXFibE4vJZOLIkSOcOnUKnU4HQEBAAHFxcQwePFjKAwtxnS05pWw+WUpOSbX1mL9ayYToAG6O8rce80hupjVSlxCkVSJcWJsSy7Zt21izZg2XL19u9PXg4GDmz58vZYJFl3alugZvD6V1HUl2URU5JdV4eygZ511BwqGNDM/fh2dwsKU10X2S5cJmWiN1FLMflFaJcFl2J5bPP/+cr776yvq1n58f3bp1A6C4uJjy8nIuX77M8uXL0Wq13H///Y6LVggXV26oJbOgjB15ZRzQVvD78T2Y1Nsyi+vuAcGMjvRl1Nl9qD9tfOBdOW5Ss62ROtIqEa7MrsRy4sQJa1IZNWoU8+bNIzY21uac3Nxc1q5dy969e1m3bh0jRoxgwIABDgtYCFdTXWNi99lyduTp2FtYQY3JDIBKAdqya5uw9g72pnewN7XLm+/qaqk1UkdaJcJV2ZVYvv/+ewCSkpJ4/PHHGz0nNjaWZ599luXLl7N9+3a+//57SSzihrb0xwIOX6wCQAEMC7+6Cj7an4B6teCtWujqktaIcHd2t1gAHnjggRbPnT9/Ptu3b+f48eNti0wIF1NrMvNT3mU2HjrPjIEh1unB43r5YzRBUow/CTEBhPi08N+qlV1d0hoR7squxKLT6fD19SU4OLjFc4ODg/H19aWsrJHttoVwEyazmeOXrq2Cv1JtWQUf5O3BgyO7AzB9QDA/GxhiOb+JnYPra21XlxDuyq7E4u3tTVVVFUajEU9Pz2bPNRqNVFVVodFomj1PCFf1+aEituSUUlRZYz3WK8ibCb18mRgbYD3WUglfm/UnSFeXuPHZlViio6M5evQo27ZtY8qUKc2eu23bNkwmE9HR0e0KUIiOkn9FT5ivJ95XpwdfKDdSVFlDqMaDpJgAJsYGMLZfT4qLixu9vtnV8NclDenqEjcyuxJLQkICR48e5Z///CdeXl4kJSU1et727dv55z//CUBiYmK7gxTC0epK8GorTaTFjCctZjy5Bk/+OyHS2hqZNTiEqX0DGdDdBzK3Yf6/T7jYXA33Vqw/EaIrsCuxTJ48mW3btpGdnc0777zDF198wdChQ23qsRw+fJgLFy4A0L9/f2699VanBC5EWxXt3M7O1H3sjL2fkwExloMG8FWaKKu3k3B0YMsbPtokl1YMygvRFdiVWJRKJYsXL2bZsmVkZWVx4cIFaxK53ujRo3nsscectrWLVqvlm2++4eTJk+Tn5zNo0CCWLFnilHuJG8srx2rJ6X0XAN61esYWHSXx4n5GUoL3f61scH5ru7hkUF4IC7tX3vv6+rJ48WKOHj1Keno6p0+f5sqVKwAEBgbSp08fEhISGDRokMODra+goICsrCz69etHTU1NyxeILqfCUEvm2XJ25Op4eHQY0VenB086l0m3oL4kXjzAmOJjeJnqFjEqGn+jVnZxyaC8EBZt3oRy8ODBDB482JGx2CU+Pp6xY8cC8Prrr8u0ZgGAvsbEnnPlbM/TsfdcBcarq+D75upYcHV68F1VJ7nrXHrDi5vqsrKji0sG5YWwM7GsWrUKgBkzZhAWFuaUgFpLdk8WdeoG4t/vlkBqxBiqVWrA0v4YGq4hKcafCb2u7R5sb5eVdHEJYR+7EktKSgoqlYpf/vKXzopHiFapNZk5fLGS/qd3Wzd01IepqVap6VdWQFJsIIkTR9FN03C9lb1dVjbnNzcrTAgB2JlYAgICqKmpkdaC6BQms5kT9VbBl1bX8kxeJhOutiTm5P3AfbkpRFSXwJnuqO5Y3eR72dtlVXd+aGgoRUUyfViI5tiVWPr06UNWVhaXL19u1bYuriQlJYWUlBQAXn31VUJDZQqouzhxsZyU7Ev8kF3EhbJr3VFRgd6YKq7Vhu9RVW/h4uUip3zGHh4e8rNjB3d7Xp0Zb0fcu6O+P7sSy5133klWVhZffvkljzzyiLNicoopU6bY7BYgf3W6j//ddIbTly0JpdvVVfBJMQH0DfHCtKOw8YuCndOykBaLfdzteXVmvB1x7/beIzIyslXn2ZVYRowYwQMPPMCaNWvQ6/XMmjWLqKioNgUoxPUulBvYkVfGzjwd/50QSa+rCxTv6BdMbmk1STEBDOzug1JxbVqwDKwL4XrsSixPPfWU5SIPD3bs2MGOHTvQaDQEBAQ0Oe6iUCh444032h+puCEVVxpJzy9jR56OE0XXasHvXPUJcxPiUI6bxO39gpq8XtaOCOF67EoshYUNux0qKyuprKx0WECtpdfrycrKAqCkpISqqioyMjIAS3VLLy+vDo9JtJ7ZbOav28+x+2w55qvHvGoN3FR0hISLBxhVcgLzaY+G26Y0QtaOCOFa7EosCxcudFYcdrty5UqDllDd1++8806nr7MRtiqNtWQWlDMh2h8vDyUKhQKNhxKVUkF8pC8Jaf9mTF4m3qZrpXwx1Da6M7AQwrXZlVimTp3qrDjsFhYWxtq1azs7DNEMfY2JPYXl7MgtY29hOYZaM14eCiZEW3YPXjCyO4vGhuOrVlH7rzSwtl3qkZ2BhXA7bd7SRYjGmM1m9hZWsCNXR8bZcqprTIBlFfyQMB98PK/VgO/uW2/xouwMLMQNQxKLaDeT2WypV3J1AP3jm/+bPB9LV2S/bt4kxQSQEONPaCOr4OvI7C4hbhxtTiz79u0jLS2t0d2NJ0yYQHx8vMOCFK7HbDaTXVzN9lwdaTlFLN39GZE6S4vjZ3k/UuLTjcSkEfRsohjc9WR2lxA3DrsTy+XLl3nrrbc4fvx4g9cqKiooLCxk586dDBgwgCeffJKQkBCHBCo6n9lsJrdUz/ZcHTvzyrhYUTfQrmJPQF9m6CyzBidr91oOX8mCViYWkNldQtwo7EosVVVVLFmyBK1WC0Dv3r2tFSTBMu33yJEjnD59mhMnTvDSSy/x6quv4uPj4/jIRbvU7Qjc2taB2WzmmU155JRcW2vSzceDxBh/Jnz6F+LKChpeJAPvQnRJdiWW9evXo9Vq0Wg0PPHEE4wePbrR8/bv38/f//53tFot69evZ/78+Q4JVjhGa0rtXiw3kpav4+4BwahVlunBkf5qLlYYSYj2JykmgEFhllXwtWurG7+RDLwL0SXZlVgyMzMBy3qWppIKwMiRI1m4cCFvv/02mZmZklhcTFOldku+/YqM4OHsyC3jeFEVAD381Yy7Wstk4ZgwfNUqPJS2lRZl4F0IUZ9diaWoqAgPDw/Gjx/f4rnjxo3j3XffdasN6LqMel1UJhT8J2IMO8NHcjioL6Y9FwFQqxTcFOVHsM+1H5FA78Z/XGTgXQhRn12JxdfXF6PR2Kp6LCqVCrVajadn01NMRePsHf+wl75bBF7F5wFQYubbXkkU+EbgYaplTC8/kmICGNvTDx/P1tfdkYF3IUQduxJLv3792LNnDxcuXCA8PLzZc7VaLZWVlda69KJ1WjP+0RaGWhN7z1WwI0/H7uFP8taeN4kos0zCuC/vB/RqDeOmTCAgYUj7vwkhRJdmVynImTNnolQqWb16NTU1NU2eV1tbywcffIBSqWTWrFntDrIraWr8w7z+E7vfq8ZkZu+5ct5KL+TnX+bw6o5zpOWXYTArOHrHryCkO6AgsaaQKbePJyDhFsd8E0KILs2uFkv//v353e9+xz/+8Q/++Mc/MnPmTIYOHUpQkGVb89LSUg4fPsw333zDhQsX+P3vf09cXJxTAr9hNTVF186pu7UmM7/++hSXKq/9AdA3xJukGH8SYwIs26lMu7k9kQohRKPsSiwLFiwALC2S/Px8li1bBmAdczGZTNZzlUol77zzDu+8806D91EoFHzyif1/gXcJbdgzq24V/K78Mh4YEYqnyrJrcP9QH7xK9UyMDSAxJoCeAWonBi6EEBZ2JRaj0djo8foJpf6xxo6L5rV26q7ZbCavVM+OPEuRrAvlls9mcJgPN0VZpgf/fnwP1CoFCoXt9GAhhHAmuxLLH//4R2fFIa5qaepujcnMuiPFbM/VcVZnsF4XfHUVfGS9VomXh11DaEII4RB2JZaRI0c6Kw5Rz/VTd0uraqgrzqtSYE0q/l6qa6vgu/ugUkrLRAjR+WTb/CY4ey1JS0qraki7Wgv++KUqVszsQ7ifGoVCwc9HdsdTpWB4hG+DVfBCCNHZJLE0wllrSVpSrq8l42wZ23N1HLpQielqQUW1SsHpy3rC/SzdXDdf3WJFCCFckSSWRjS7lsRJicVYa+KRr09RabRMePBQwugeviTFBnBTlB+aepUXhRDClUliaYyD1pI0xVBrYl9hBZlny/ntTRF4qhR4qpSM6uFLmaGWiTEBjOvlj7+XJBMhhPuRxNIYJ9RfrzGZOai1bKmSUVBubZkkRPszpqcfAP+dECkD8EIItyeJpRGO3AZeX2Pig30XSc8vQ6evtR7vE+xFUmwAvYO9rMckqQghbgSSWBrRnm3gzWYzBVcMRAdZEoZapeCgthKdvpaoADVJsQEkxvgTFeDVwjsJIYR7ksTSBHu3gc+z1oLXoS03snJmX8L8PFEoFDw6NpxAbxWxQV6yCl4IccOTxNIO58sM7MjVsSNPR/6VeqvgvVWcLzcQ5mepRTOyh29nhSiEEB2u3YmluroavV5PYGCgI+JxG1VGE09sOIPx6mITf7WSCdGWbq4hYRoZLxFCdFltSiz5+fmsX7+eAwcOUFFRgUKh4LPPPrO+Xl5ezhdffIFCoWDBggV4eLh3w+hKdQ3p+WXsLSxn8cQoPJQKfDyVJMb4YwaSYgIYEeGLp0qSiRBC2P0bPz09neXLl9sU+jKbzTbn+Pn5kZ+fz9GjRxkxYgSjRo1qf6QdrNxQS2ZBGTvyyjigrbCugj+orWB0pGV68JMTIjsxQiGEcE12JZbz58/z7rvvUlNTw2233cYtt9zC3/72N8rKyhqcO2nSJI4ePcru3bvdKrFUGmt5K/08ewsrqLmaTVQKiI/0JSkmgIHdfTo5QiGEcG12JZZvv/0Wo9HIbbfdxqJFi4BrRb6uN2zYMACys7PbGaJzGWtNZBdVMyRcA4CPh5K8Uj21JjPDwjUkxQQwPtqfAFkFL4QQrWJXYjl8+DAA99xzT4vnhoSEoFarKSpyzDYojrb/fAXbc3VkFJRRaTSxanZfQjWW6cFPTYiku68H3TSenR2mEEK4HbsSS0lJCd7e3oSGtm5rE7VaTWVlZZsCc7YX/1Ng/XfvYC8uV9UQejWRSHeXEEK0nV2JRaVS2QzaN8dgMFBZWYlGo2lTYM4W6e9JUmwASTEB9AqUVfBCCOEodiWWsLAw8vPz0Wq1RERENHvu/v37MZlM9OzZs10BOsu7P+sjq+CFEMIJ7CqKPnz4cAA2bdrU7HkVFRV8+umnAC47I0ySihBCOIddieXuu+9GrVazadMmkpOTMRgMNq/X1NSQmZnJH//4R7RaLX5+fkybNs2hAQshhHBtCvP1qxtbkJmZyVtvvYXJZEKtVlNTU4PJZCI2Npbz58+j11u2mlepVDz77LOMHDnSKYG3V2FhYWeHINxQaGioy850dEXu9rw6M96OuHd77xEZ2bpF4Xa1WABuvvlmXnzxRaKjozEYDJhMloJVubm51qQSFRXFCy+84LJJRQghhPO0aROvgQMH8re//Y2cnByOHz9OSUkJJpOJoKAgBg4cyIABA2QMQwghuqh27Q4ZFxdHXFyco2IRQghxA7CrK+ypp57i6aefRqvVOiseIYQQbs6uFsvFixdRqVQtrmERQgjRddnVYgkKCpKxEyGEEM2yK7EMHTqU6upq8vLynBWPEEIIN2dXYpk5cyZqtZoPP/wQo9HorJiEEEK4MbvGWHx9ffnNb37DihUreOaZZ7j77rvp378/gYGBTdZlAQgMDGx3oEIIIdyDXYmlrrgXgFarZfXq1S1eo1Ao+Oyzz+yPTAghhFuye+W9vezcMUYIIYSbs6vF8uabbzorDiGEEDcIuxJLazcgE0II0XU5vStMCCFE1yKJRQghhEPZ1RWWnp7epptMmDChTdcJIYRwP3Yllr///e9230ChUEhiEUKILsSuxBIQENDsXmFVVVXWcsVqtRofH5/2RSeEEMLt2JVYVq5c2eI5+fn5rFu3jn379rFw4UJuuummNgcnhBDC/Th88D46OpqnnnqKm2++mbfffpv8/HxH30IIIYQLc9qssLlz52I0Glm/fr2zbiGEEMIFOS2xdO/eHY1Gw5EjR5x1CyGEEC6oXTXvm2MwGKiurrYO5gshhOganNZi2bFjByaTieDgYGfdQgghhAuyq8Vy5cqVZl83Go0UFRWRmZnJ1q1bARg7dmzboxNCCOF22lyPpTUiIyO599577bpGCCGEe3PKGEtISAgJCQncc889aDQaZ9xCCCGEi3JoPRalUomfnx9+fn7tCkoIIYT7knosQgghHEq2zRdCCOFQkliEEEI4VJNdYW2tvdIY2TZfCCG6jiYTS1tqrzRG6rEIIUTX0mRiaan2ihBCCNGYJhNLa2qvCCGEENeTwXshhBAOJYlFCCGEQ0liEUII4VAtTjfWaDSMHDnS5pi9ZFaYEEJ0HS1ON46MjLQmlrZMQZbpxkII0bW0ON04ICCgwTEhhBCiKQqz2Wzu7CCEEELcOGTwXgg7LF68uLNDcCvu9rw6M96OuHdHfX+SWIQQQjiUJBYhhBAO1ebSxDqdjuzsbIqLi9Hr9TQ3VDNz5sy23kYIlzJlypTODsGtuNvz6sx4O+LeHfX92T14X15ezgcffMCuXbswmUytuubzzz9vU3BCCCHcj10tFoPBwEsvvUR+fj5KpZLo6Gjy8/Px8PAgOjqa0tJSSkpKAPD19aVHjx5OCVoIIYTrsiuxbNmyhfz8fCIiIvjzn/9MaGgoc+fOxc/Pj1deeQWAwsJCPv30U/bt28eECRO4++67nRK4EM6g1Wr55ptvOHnyJPn5+QwaNIglS5bYnGM2m1m/fj1bt25Fp9MRFxfHww8/TGxsbKfE3Fl27drF9u3bOX36NJWVlURGRvKzn/2MxMREm/NSUlL45ptvKC4uJioqigULFjBs2LBOiTkjI4MNGzZQWFiIXq8nNDSUiRMnMnPmTDw8LL8OO+rzLSkp4fe//z16vZ6PP/4Yb2/vdt0/NTWVd999t8HxhQsXMm3atA793uxKLJmZmQDMnz+f0NDQRs+JjIzkD3/4A2+88QaffPIJvXv3ZvDgwe2PVIgOUFBQQFZWFv369aOmpqbRc5KTk1m3bh0PPvggkZGRbNy4kaVLl/L6668TFBTUwRF3ng0bNhAWFsYvfvELAgIC2LdvH2+//TZlZWXceeedAKSlpbFy5Uruu+8+Bg4cSGpqKq+++iqvvPIK0dHRHR5zWVkZQ4YMYcaMGWg0GnJycvjiiy8oLS3lV7/6FdBxn+8nn3yCt7c3er3e5nh77//CCy+gVqutX4eHhzvsvVvLrsRy9uxZAEaNGmVzvLH/gPPnzyczM5PvvvtOEotwG/Hx8YwdOxaA119/nbKyMpvXDQYDycnJzJ49mzvuuAOA/v3789hjj7Fp0ybmzZvX4TF3lueee85mZ46hQ4dy+fJlNmzYYE0sa9eu5ZZbbmHOnDkADB48mDNnzpCcnMzvfve7Do956tSpNl8PHTqUqqoqNm/ezC9/+UuMRmOHfL7Hjh1j//79zJ49m3/961/W4474+YqLi7O2furryJ9du6YbGwwGNBqNTTb09PRskHEBIiIi0Gg0nDx5sv1RCtFBlMrm/0tkZ2dTVVXF+PHjrce8vb2Jj48nKyvL2eG5lPpJpU7v3r3R6XQAXLhwgfPnz9vsFahUKhk/fjz79+/vsDhb4u/vb/3juCM+X5PJxAcffMCcOXMaPENn3r8jf3btSiyBgYFUV1fbTC329/fHaDRaB+3rmEwm9Ho95eXljolUCBdw7tw5lEplg4kpUVFRFBYWdlJUruPEiRNERUUBlmcFlu7x+nr27El5ebk1AXWGut9Px48f5/vvv2fatGkoFIoO+Xy3bNmC0Wjk9ttvb/CaI+7/xBNPMG/ePH7/+9+zdetWh753a9nVFRYaGkpxcTGXL18mJCQEgJiYGEpKStizZ491gAhg37591NbWdqk+Z3Hjq6iowNvbu0HLxtfXF71eT01NjXUQuKs5dOgQe/bs4Te/+Q1geVZgeTb11X1dXl7eaKunIzz44IMYjUYAJk6cyIIFCwDnf75lZWV8/vnnPPHEE42+T3vuHxQUxNy5c4mLi8NkMlnHt/R6PdOnT+/Qn1273mXo0KGcOHGCI0eOkJSUBMD48ePJysri008/pba2ltjYWPLy8li7di0AI0aMcEigQriKxnb47up7uV68eJG3336bMWPGMGnSJJvXmtoRvTN3Sl+6dCkGg4GcnBy+/PJLPvjgAxYuXNhkXI76fNesWUO/fv0YPXp0k+e09f4jR460ljgBy1i40Wjkq6++4q677mrXe9vLrsQyZswY1q1bx86dO62JJSkpiZSUFLKzs/nnP/9pc76/vz/333+/w4IVorP5+vpSVVWFyWSy+cuvsrISLy+vLtlaKS8v55VXXiE0NJQnnnjCeryuZVJRUYFGo7Eeb6ol05H69OkDwMCBA/H392f58uVMnz7dqZ9vQUEBP/74Iy+99JL1GdSNT1dWVqJUKh1+/3HjxrFr1y4uXbrUoT+7dr1Tnz59+Ne//mWT9ZRKJX/6059Yu3Yt6enpXL58GW9vb0aMGNHstGQh3FHPnj0xmUxotVqbsYNz5841GEvoCvR6Pa+++io1NTUsXrzYZjZSz549Acuz6d69u/X4uXPn8PPz67RusOv17t0bsLS6nPn5nj9/ntraWp5//vkGr/36179m8uTJJCYmOuX+CoWiQ3927U5Rnp6eDY55e3vz85//nJ///OcNsqEQN5L+/fvj4+PDrl27uPfeewHLL9e9e/e63b5Y7VVbW8sbb7zB+fPnWbp0KYGBgTavh4eH06NHDzIyMqxdNCaTyeZrV3DixAkAwsLCCAkJcdrnO3DgQF588UWbY/v37+frr7/mj3/8I+Hh4YSGhjr0/pmZmfj7+xMaGkpQUFCH/ew6vN0uSUW4M71eb516WVJSQlVVFRkZGYClz9rLy4tZs2axbt06fH196dmzJxs2bMBsNlvXBnQVq1atIisri4ceeojy8nKys7Otr/Xu3RtPT0/uu+8+li1bRvfu3RkwYADbtm3j/PnznbKGBeDll19m2LBh9OrVC6VSyYkTJ/j222+ZMGECERERAE77fAMCAhgyZIjNsUuXLgEwaNAga2uvrfd/7bXXiIuLIyYmBpPJRHp6Ounp6Tz88MMolUrUanWH/ex2vQ5hIZpx5coV3njjDZtjdV+/8847hIWFMWvWLMxmM8nJyZSVldG3b1+ef/75LjcD8uDBgwANxlbh2rNKTEykurqar7/+mnXr1tGrVy8WL17cKavuAfr27cu2bdu4ePEiKpWK8PBw5s+fb7NwsrM/37bePzIykh9//JHi4mLMZjNRUVE8/vjjTJw4sd3vba9mdzeeO3cuQUFBrFixosFrZ8+epba2lpiYGIcGJIQQwr21ucXyl7/8BZ1Ox2effebIeIQQQri5dg2IdPW5+0IIIRqSkXYhhBAOJYlFCCGEQ0liEUII4VCSWIQQQjiUJBYhWunIkSPcf//9sv+dm0pNTeX+++/nscce6+xQbngtTjcuLS1l7ty5Tb7e3Gtg2aNGpiR3HoPBwLZt29i7dy95eXnodDo8PDwICQlh4MCBJCQkMHTo0M4Os1NVVFSwceNGAO6+++5O3RxR2E8+P9fj9JX3MiW58xw8eJD33nuP4uJi6zEfHx9qamo4d+4c586d44cffmDUqFE8/vjj+Pv7d2K0naeiooIvv/wSgEmTJjX5i8nLy6tLbjTp6lr7+YmO02xiqatTLdxPeno6y5Yto7a2lpCQEO6//35uuukm/Pz8AMuOplu3bmXz5s1kZWXxpz/9qdGNBMU1cXFxvPXWW50dhhAur9nEct9993VUHMKBzp07x3vvvUdtbS3R0dG88MILDbYo79mzJw899BDDhw/ntddeQ6vV8vbbb/PnP/+5k6IWQtwoZBPKG9CaNWvQ6/V4enry9NNPN1v3YvTo0dxzzz2sXbuWQ4cOsW/fPpvqdhcvXuTxxx8HLBsL1tbW8tVXX3Ho0CF0Oh2BgYGMGjWKOXPmWMtVN+Wnn34iNTWVU6dOodPp8Pb2Jjo6moSEBCZPntxooaElS5Zw9OhR5syZwz333MP3339PWloaWq2WyspKXnzxRYYMGYLJZCI7O5u9e/dy9OhRiouLuXLlCj4+PvTq1avJe9S9f52677XO4MGDWbJkCWAZvH/ppZcArBVS66SmpvLuu+/SvXt3li9fzunTp0lOTub48eOUl5cTEhLC2LFjuffee62txsYcPXqUb775hpMnT1JdXU1oaCjjxo1j9uzZZGRk2NzDHsuXL2fbtm3ccsstPPbYY6SmprJ161bOnj2LUqmkT58+3HvvvQwePBiwbIm/ZcsWUlNT0Wq1gGXb97lz51qLZDWmsrKS7777jt27d6PVaqmpqaFbt24MGzaMGTNmEB4e3uh1dRMiXnzxRfr06cPXX39NRkYGly5dwsvLi/79+3PvvffSr18/m+vs+fyu15bP6OTJk3z//fecOHGC0tJSlEol/v7+dO/enWHDhnHrrbfSrVu3Jp9PV7Xv9GMAAA2USURBVCGJ5QZz+fJldu/eDUBCQkKrxgSmT5/Ot99+S1VVFZs3b26ybGpOTg4rVqygqqrKWju7uLiYlJQUMjIyeP755xv9pVNdXc1bb73Fvn37rMd8fHyorKzk2LFjHDt2jO3bt7N48eIm/0MbjUZeeuklTpw4gUqlsikoBVBUVMQLL7xg/VqlUqFWqykvL7feIy0tjT/96U+o1WrreX5+fvj7+1NWVgZYqp7WL/3QXBJoys6dO1m+fDm1tbVoNBpqa2u5ePEiGzdu5ODBg7z88ssN4gf4/vvv+ec//2kdl9RoNFy6dIn169fz008/OaxmRl2SqXtGFRUVHDp0iKNHj/LMM88wfPhw/t//+38cOHAADw8PVCqVtZzA0aNHeemllxr9nAsKCvjrX/9qHdPz9PTEw8MDrVaLVqslNTWVJ554gnHjxjUZ2+XLl3nuuefQarV4enqiVCopLy9n3759HDhwgOeee86mlktbP7+2fEapqam899571s+nLr6ioiKKioo4duwYoaGhDUozd0WSWG4wR44csf7g33zzza26xtvbm+HDh5OZmcmxY8eora1FpVI1OO/9998nLCyMRx99lLi4OMxmMwcPHmTFihUUFRXx2muv8frrr+Pj42Nz3bJly9i3bx8RERHcf//9xMfH4+Pjg8Fg4ODBg3z00UdkZ2fz3nvv8Yc//KHRGDdv3gzAb3/7WyZMmIBaraasrMxazVSlUjFmzBgSExMZOHAgQUFBKJVKqqurycjIYM2aNRw7dow1a9bwi1/8wvq+zzzzjE2r7JVXXiEsLKxVz60xOp2O9957j1tuuYU5c+YQGhqKXq/nxx9/5KOPPqKgoICvv/66wWzKEydOWJPK8OHD+eUvf0lkZCS1tbXs3r2b999/3zpA3R579uzBaDSyaNEiJk6ciFqtprCwkLfffpvTp0/zwQcfEB8fz6lTp3jqqacYO3YsKpWKM2fO8Oabb3LhwgU+/PBDli5davO+VVVV/N///R/FxcWEhISwaNEiRo4ciVKpJDc3l5UrV3Ly5EmWLVtGREQEsbGxjca3evVqgoODeeGFFxg8eDAKhYJTp07xzjvvUFhYyMqVK1m2bJk1ebTl82vLZ6TX6/nwww8xm80kJSVx3333Weu3VFdXU1hYSHp6ustUxexsso7lBnP27Fnrv+tKrrZG3X/06upqa/Gh66lUKp5//nni4uIAy1TyESNG8D//8z94eHhQVFTE1q1bba7Zt28fu3fvJigoiCVLlpCYmGhNPGq1mjFjxrBkyRK8vLzYvXs3ubm5jd67urqa3/3ud0yaNMna4vD397f+RdqtWzeeffZZJkyYQEhIiPUXj7e3N5MmTeLZZ58FICUlBYPB0OrnYi+9Xs+ECRP49a9/bS3L7eXlxR133MGdd94JQFpaWoPr1q5da62h8dxzz1lbmiqVinHjxvH0009b66S3R0VFBY8++ihTpkyxPsfIyEieeuopFAoFly5dYtOmTTz77LOMHz8eDw8PFAoFffr04dFHHwUsSbD+TEOwJP66Gif/8z//w+jRo62fQWxsLM8//zzdu3fHaDQ2u/xApVLx4osvMnToUJRKJQqFgri4OJ5++mnAUhirfkGxtmjLZ1RQUEBVVRVeXl789re/tSYVsPyM9enThwULFjTZ2u9qJLHcYOq6BMC+bpz6U43Ly8sbPWfq1KmNzhqLioqydm9c/x/yhx9+AGDixIlNjsF069bNWllv//79jZ7Tq1cvxowZ08J30bS+ffsSGBiIXq9vMnk5Sl3Z1+vVxa/VatHr9dbj5eXlHD58GIAZM2Y0Wv576NChDBo0qN2xhYaGkpiY2OB4eHi4dfxj0KBBDBw4sME5gwcPtsaWl5dn89quXbsA+P/tnW1IU18cx7/q0k1nNsus1JVmz2WTfKKCMJoEJdWbCoqIIElKNIiiIoJKRpES9SKjCBaxGIHREwQrLbIHozJzlKtQmqklNXNNHe7O/V+MHTa7d865dPr/fV6NnXvuzr2/u/t7POfk5OTwbuIlkUiwYcMGAEBdXR16enp4x7dmzRreZ0wulzNPxGg0Cl6frwxVRpGRkQAAjuM8/mMEPxQKI3zG20TKxYsXo6amBkajERzHsSR5Y2MjAKen8OTJE8H+rhfNz58/edvnzZs36Pg4jkNVVRVevXqFlpYWWCwW2Gy2v44zmUyDnstfpFKphzXrjrti7e7uRkREBACgubmZhS9dyXM+Fi5ciI8fPw5rfCkpKSx8OJCYmBh8//4ds2fP5m13JapNJpOH98RxHFM0S5YsEfzttLQ0AM65bU1NTbzPk8sb5kMmk6Gjo0PQ8PEVf2Q0bdo0JCQkoLW1FUePHoVSqYRCoYBcLqft2HkgxTLOGOh5DFap5cIXT8fbuVxtdrsdFosFkyZN8rDuhCzUgbhbie4MFrvu6urCyZMnPazZCRMmeCRzzWYzHA4HrFarT2Pxh4H5JXfcX0Acx7HPZrOZfZbJZIL9fZWlN7yNz5VX4yssGHiM3W5n31ksFvT39w86Rvc292se6vjc750/+COj0NBQFBcX4+zZs+jo6IBGo4FGo2EVa9nZ2Vi1ahVTRP93SLGMMxITE9nnpqYmn19Gzc3NAJwvlbi4uICMxfWyAYCSkhIsX77c73MNZhWq1WoYjUZER0dj+/btSE9P/2sf78LCwr9yA8GA++oUQt7EwOOCFW/j99Y2Fpg1axbOnTuHN2/eoL6+HgaDAd++fUNDQwMaGhpw69YtHD58mDcU+H+DfLhxxqJFi9gfuLa21qc+VqsVDQ0NAJzxdb6KMMB7CMnVFhYWxjye8PBwFpsORFxcCI7j2LXu2rULubm5fymV/v5+QSt5tHHPKXi7x52dnSMxnCEjlUqZ4vemuN3bxmr1lEgkQnZ2NgoKClBWVoYrV65g9+7dkEql+PXr15DnF41XSLGMM2QyGTIzMwE4l3Vpa2sbtM+9e/fQ29sLAMjLyxM8zpVg9tYml8s9JiG6ciMvXrzw8GACidlsZrkUoUq4xsZG3nwLMLg39K9JTk5mxoD7ZL+BeGsbTUQiEWbOnAkAzEDh4/379wDAqswCxWjKLzo6GkqlEtu2bQPg9PwpuU+KZVyyZcsWhIeHw2azoby83KulXldXh8rKSgBOb8dbuaROp+M9V1tbG/MYBoa7XJP62tvbcefOHa/jtlqtfsXPIyMj2YuZr+LLbrfjxo0bgv3dY+6BKOkdKlKplFXF3b17l/cefPjwYdiJ+3+JS+4vX77k9U6tViuTf3p6OvNkA8FIyE/IKHHhPul2tA2VYIDuwDgkKSkJe/bsQWhoKIxGIw4dOoSqqiqPP11bWxvUajXOnDkDjuMQHx+P4uJir3Fwu92OU6dO4cuXLwDAJkiWlpbCZrNh8uTJUCqVHn0yMzORlZUFANBoNLh8+bKHF8VxHD5//ozr169j79696OrqGvL1isVi5hldu3YNer2eeUdGoxEqlQpNTU2CidWoqCiWi6qurvZITI8UmzdvRkhICFpaWnD69Gm0t7cDcN7z2tpalJWVBfWqvXl5eZg6dSrsdjtUKhXq6uo8ZFBaWoqOjg6IRCJs3bo1oL89EvJ79uwZjh07Bp1Ohx8/frDv+/v78e7dO2g0GgDA3Llzg1pOIwUl78cpK1euhFQqZcvmV1RUoKKiApGRkbDZbB4W2NKlS1FUVDRo3LugoACXLl3CkSNHIBaL4XA4WBVXVFQUDhw4wGuJFhUV4eLFi3j+/Dl0Oh10Oh0iIiIgEonQ09Pjc/LaGzt37sTx48dhMplw4sQJtpxIb28vwsLCUFhYCK1WKzj5U6lUQqvV4sGDB3j06BEmTpyI0NBQzJkzByUlJX6NaSjMnz8fO3bsgFqtRn19PYqLixEVFYW+vj7YbDYkJSVh9erVUKvVvPNcRhuJRIKDBw+yJV1UKpWHDABnlV5RUZHgrPvhMBLyMxgMMBgMAJzXIhaLYbFY2PMrk8lQWFgYkN8a65BiGccoFApcuHABjx8/Zht9/fnzByKRCFOmTMGCBQuwYsUKr3MP3ElNTYVKpUJlZSX0ej3MZjNiY2PZIpRCi+9FRESgpKQESqUS1dXVMBgM6OzshNVqRUxMDBITE6FQKJCVleV3SW1KSgpUKhVu3rwJvV6Pnp4eiMViKBQK5OfnIzU1FVqtVrD/pk2bIJFI8PTpU7S2tsJkMsHhcASsQs4X1q1bh+TkZNy+fRufPn1CX18f4uLikJOTg40bN+Lhw4cAELQWsVwuR3l5Oe7fv++xCGV8fDzS0tKQn58vOH9kuPxr+WVkZGDfvn3Q6/Vobm7G79+/YbFYIJFIMGPGDCxbtgxr164NWtmMNCGOsVDDSIwaA1c3Hs46WsTwOH/+PGpqapCbm0uWMRHUUI6FIMYA7gUS7qv7EkQwQqEwgggStFotYmJikJGRwRbStFqtePv2LdRqNWw2GxISElg5OUEEK6RYCCJI+Pr1K16/fo2rV68iLCwMEokE3d3dLDkcGxuL/fv3826IRhDBBD2hBBEkrF+/HrGxsWx3QldyePr06Sw57M/GYwQx0lDyniAIgggolLwnCIIgAgopFoIgCCKgkGIhCIIgAgopFoIgCCKgkGIhCIIgAgopFoIgCCKg/AeWDEuXjWo2nAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "import logging\n",
    "\n",
    "class Input_builder(object):\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def __call__(self, model,x,y=None,train_window=20,train_window_2=None):\n",
    "        if model=='weibull':\n",
    "            return self.create_weibull_input(x,y,train_window)\n",
    "        elif model=='svm' or model=='lstm':\n",
    "            return self.create_RNN_input(x,train_window=20)\n",
    "        elif model=='seq2seq':\n",
    "            return self.create_seq2seq_basic_input(x,train_window,train_window_2)\n",
    "        elif str(model)=='arima':\n",
    "            return x.iloc[:,-1].values\n",
    "        elif str(model)=='xgb':\n",
    "            return self.create_xgb_input(x)\n",
    "\n",
    "    def create_weibull_input(self,x,y,train_windows=20):\n",
    "        index_end=len(y)-1\n",
    "        y=list(y)\n",
    "        for yy in y[::-1]:\n",
    "            if yy!=y[-1]:\n",
    "                index_end=y.index(yy)\n",
    "                break\n",
    "        index_begin=index_end-train_windows if (index_end-train_windows>0) else 1\n",
    "        x,y=x[index_begin:index_end],y[index_begin:index_end]\n",
    "        logging.info(\"Weibull train data {}\".format(len(x)))\n",
    "        return np.array(x),np.array(y)\n",
    "        \n",
    "        \n",
    "alpha=10e-7\n",
    "class Weibull_model(object):\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def train(self,x,y,train_window):\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        train_x, train_y = Input_builder().create_weibull_input(x, y, train_windows=train_window)\n",
    "        self.x_weibull = np.log(train_x)\n",
    "        self.y_weibull = np.log(-np.log(1 - train_y) + alpha)\n",
    "        self.x_weibull=self.x_weibull.reshape(-1,1)\n",
    "        self.model=linear_model.LinearRegression()\n",
    "        self.model.fit(self.x_weibull,self.y_weibull)\n",
    "        print(\"Weibull model train finished with score: {}%\".format(100*self.model.score(self.x_weibull,self.y_weibull)))\n",
    "        self.weight_weibull=self.model.coef_[0]\n",
    "        self.bias_weibull=self.model.intercept_\n",
    "        print(\"Weight: %s, Bias: %s\" % (self.weight_weibull,self.bias_weibull))\n",
    "        return self.weight_weibull,self.bias_weibull\n",
    "\n",
    "    def predict_by_interval(self,predicted_interval,return_full=True):\n",
    "        x_max=max(self.x)\n",
    "        #assert x_max< (len(self.x)+5)\n",
    "        x_future=[i for i in np.arange(x_max+1,x_max+predicted_interval+1)]\n",
    "        x_future_weibull=np.log(x_future).reshape(-1,1)\n",
    "        y_predict_weibull=self.weight_weibull*x_future_weibull+self.bias_weibull\n",
    "        y_predict=1.0-1.0/(np.exp(np.exp(y_predict_weibull-alpha)))\n",
    "        y_predict=y_predict.reshape(y_predict.shape[0],)\n",
    "\n",
    "        if return_full:\n",
    "            self.x_future = list(self.x) + list(x_future)\n",
    "            self.y_future = list(self.y) + list(y_predict)\n",
    "            return self.x_future, self.y_future\n",
    "        else:\n",
    "            return list(x_future),list(y_predict)\n",
    "\n",
    "\n",
    "    def predict_by_calendar(self,interval_samples,calendar_failures,predicted_interval,output_file,return_full=True,include_future_samples=False):\n",
    "        \"\"\"\n",
    "        interval_samples: [sample0,sample1,sample2] cumulative value for each operating time\n",
    "        calendar_failures: {Date1: failure1,Date2: failure2}\n",
    "        \"\"\"\n",
    "        x_future,y_future=self.predict_by_interval(predicted_interval+1)\n",
    "        failure_rate_interval =[y_future[i+1]-y_future[i] for i in range(len(y_future)-1)]\n",
    "        failure_rate_interval=[max(alpha,x) for x in failure_rate_interval]\n",
    "        logging.info(\"Use {} data to further predict {} future\".format(len(interval_samples),predicted_interval))\n",
    "        assert len(interval_samples)+predicted_interval==len(failure_rate_interval)\n",
    "\n",
    "        if return_full:\n",
    "            pass\n",
    "        else:\n",
    "            calendar_failures=pd.DataFrame()\n",
    "\n",
    "        if include_future_samples:\n",
    "            #Todo : check if it's better or not\n",
    "            samples_future=interval_samples+[np.mean(interval_samples)]*predicted_interval\n",
    "\n",
    "        else:\n",
    "            for i in range(1,predicted_interval+1):\n",
    "                samples_failure_rate_interval=failure_rate_interval[i:i+len(interval_samples)]\n",
    "                failure_interval=sum(interval_samples*samples_failure_rate_interval)\n",
    "                calendar_failures=calendar_failures.append({\"Failures\":failure_interval},ignore_index=True)\n",
    "        calendar_failures.to_csv(output_file,index=False)\n",
    "        return calendar_failures\n",
    "\n",
    "    def plot(self):\n",
    "        plt.style.use('ggplot')\n",
    "        fig,ax=plt.subplots()\n",
    "        x, y = np.array(self.x), np.array(self.y)\n",
    "        x_weibull,y_weibull=np.log(x),np.log(-np.log(1 - y) + alpha)\n",
    "        ax.plot(x_weibull,y_weibull,marker='o',linestyle='')\n",
    "        x_plot = np.arange(np.log(1), np.log(len(x)), np.log(2))\n",
    "        ax.plot(x_plot,self.weight_weibull*x_plot+self.bias_weibull,'--',linewidth=2)\n",
    "\n",
    "        ax.set_yticks(list(map(lambda y: np.log(-np.log(1 - y) + alpha),np.array([0.01,0.05]+[i/100 for i in range(10,100,20)]))))\n",
    "        ax.set_yticklabels(np.array([1, 5] + [i for i in range(10, 100, 20)]), fontsize=15)\n",
    "        ax.set_xticks(list(map(lambda x: np.log(x),[i*10 for i in range(1, 10)])))\n",
    "        ax.set_xticklabels([i*10 for i in range(1, 10)])\n",
    "        ax.set_xlim([1, np.log(len(self.x)+10)])\n",
    "        ax.tick_params(axis='both', which='major', labelsize=15)\n",
    "        ax.set_xlabel('Operating months', fontsize=25)\n",
    "        ax.set_ylabel('Failure probability [%]', fontsize=25)\n",
    "        plt.show()\n",
    "\n",
    "        \n",
    "x,y=life_table['time'].values,life_table['failure_probability'].values\n",
    "weibull_model = Weibull_model()\n",
    "weibull_model.train(x,y,train_window=50)\n",
    "x_full_pred,y_full_pred=weibull_model.predict_by_interval(predicted_interval=len(x),return_full=False)\n",
    "weibull_model.plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# log-rank test\n",
    "- To compare the survival times of two or more groups\n",
    "- univariate analysis for categorical groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# cox proportional hazards regression\n",
    "- To describe the effect of categorical or quantitative variables on survival\n",
    "- assume that the effects of the predictor variables upon survival are constant over time and are additive in one scale, the hazard for any individual is a fixed proportion of the hazard for any other individual\n",
    "- A semi-parametric model, multivariate regression\n",
    "- $ h(t|x)=h_0(t) \\times exp(b_1x_1+b_2x_2+ \\ldots +b_px_p) $ \n",
    "- Cumulative hazard at a time t is the risk of dying between time 0 and time t, and the survivor function at time t is the probability of surviving to time t\n",
    "- cox regression aims to estimate the hazard ratio, while logistic regression aims to estimate the odds ratio\n",
    "- [reference](http://courses.washington.edu/b515/l17.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tree based survival model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# mixed weibull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian restoration maximization\n",
    "- [Bayesian estimation of Weibull mixture in heavily censored data setting](https://hal.archives-ouvertes.fr/hal-01645618/document)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}