diff --git a/02_kaplan_meier.ipynb b/02_kaplan_meier.ipynb new file mode 100644 index 0000000..f479662 --- /dev/null +++ b/02_kaplan_meier.ipynb @@ -0,0 +1,3007 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "view-in-github", + "colab_type": "text" + }, + "source": [ + "\"Open" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "q403LSMFTTlt" + }, + "source": [ + "# Kaplan-Meier Estimation" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_hSUxOBlTTlu" + }, + "source": [ + "[Run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/SurvivalAnalysisPython/blob/master/02_kaplan_meier.ipynb)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "IXpLRANsTTlv" + }, + "source": [ + "This notebook introduces Kaplan-Meier estimation, a way to estimate a hazard function when the dataset includes both complete and incomplete cases.\n", + "To demonstrate, I'll use a small set of hypothetical data. " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "jVxKfh1hTTlv" + }, + "source": [ + "## Dog adoption data\n", + "\n", + "Suppose you are investigating the time it takes for dogs to get adopted from a shelter. You visit a shelter every week for 10 weeks, and record the arrival time for each dog and the adoption time for each dog that was adopted.\n", + "\n", + "Here's what the data might look like. " + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "id": "DB1SvRfYTTlw" + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 269 + }, + "id": "rr3TgmIQTTlx", + "outputId": "60d67ffd-8571-4d60-b51a-66fbc945d484" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " start end status\n", + "0 0 5 1\n", + "1 1 2 1\n", + "2 2 6 1\n", + "3 2 9 0\n", + "4 4 9 0\n", + "5 6 8 1\n", + "6 7 9 0" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
startendstatus
0051
1121
2261
3290
4490
5681
6790
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "obs", + "summary": "{\n \"name\": \"obs\",\n \"rows\": 7,\n \"fields\": [\n {\n \"column\": \"start\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 0,\n \"max\": 7,\n \"num_unique_values\": 6,\n \"samples\": [\n 0,\n 1,\n 7\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"end\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 2,\n \"max\": 9,\n \"num_unique_values\": 5,\n \"samples\": [\n 2,\n 8,\n 6\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"status\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 24 + } + ], + "source": [ + "obs = pd.DataFrame()\n", + "\n", + "obs['start'] = 0,1,2,2,4,6,7\n", + "obs['end'] = 5,2,6,9,9,8,9\n", + "obs['status'] = 1,1,1,0,0,1,0\n", + "\n", + "obs" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "N8FlCd7yTTlx" + }, + "source": [ + "This `DataFrame` contains one row for each dog and three columns:\n", + "\n", + "* `start`: arrival time, in weeks since the beginning of the study\n", + "\n", + "* `end`: adoption date, for dogs that were adopted, or `9` for dogs that had not been adopted at the end of the study\n", + "\n", + "* `status`: `1` for dogs that were adopted; `0` for dogs that were not." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "oSpxyFXUTTlx" + }, + "source": [ + "## Plotting lifelines\n", + "\n", + "The following function visualizes the data." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "id": "GxBWpgBtTTly" + }, + "outputs": [], + "source": [ + "def plot_lifelines(obs):\n", + " \"\"\"Plot a line for each observation.\n", + "\n", + " obs: DataFrame\n", + " \"\"\"\n", + " for y, row in obs.iterrows():\n", + " start = row['start']\n", + " end = row['end']\n", + " status = row['status']\n", + "\n", + " if status == 0:\n", + " # ongoing\n", + " plt.hlines(y, start, end, color='C0')\n", + " else:\n", + " # complete\n", + " plt.hlines(y, start, end, color='C1')\n", + " plt.plot(end, y, marker='o', color='C1')\n", + "\n", + " plt.xlabel('Time (weeks)')\n", + " plt.ylabel('Dog index')\n", + " plt.gca().invert_yaxis()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "1sPBWtrSTTlz" + }, + "source": [ + "Here are the results:" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 449 + }, + "id": "e_wzTSo3TTlz", + "outputId": "c5b62064-1f3d-4b48-d86d-b806b57985b5" + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJ/1JREFUeJzt3Xl0FGWi/vGnk5BFkm6RJSSXBllGISwCEjGArJHlIIgzwpVRtkHPHQyBiHoHmKvATyXgHR3mgIPAeFhmwBVRLiOCooR9gCAIRnBAWS6rDJJOgjaSrt8fkb6GhKVDJ/U2/f2c04fUW9VdD7TSD1VvVTssy7IEAABgoAi7AwAAAFwORQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFhRdge4Hj6fT8eOHVNCQoIcDofdcQAAwDWwLEsFBQVKTk5WRMSVj5mEdFE5duyY3G633TEAAEAFHDlyRPXq1bviNiFdVBISEiSV/EadTqfNaQAAwLXweDxyu93+z/ErCemicvF0j9PppKgAABBirmXaBpNpAQCAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxQvrOtABQYb5i6dAmqfCkFJ8oNeggRUTanQrAJYw4ovLKK6/o1ltvVWxsrNq3b6+tW7faHQnAjSxvuTSjhbTwPmnpyJJfZ7QoGQdgFNuLyptvvqlx48Zp0qRJ2rFjh+644w716tVLp06dsjsagBtR3nLpraGS51jpcc/xknHKCmAUh2VZlp0B2rdvr9TUVM2aNUuS5PP55Ha7lZmZqfHjx1/xuR6PRy6XS/n5+ZXzpYTni4L/mgDs4yuWXrlLKjh+mQ0ckjNZytrNaSCgEgXy+W3rHJXz588rNzdXEyZM8I9FREQoPT1dmzdvLrO91+uV1+v1L3s8nsoNODW5cl8fgGEsyXO0ZO5Kw3vsDgNANp/6OX36tIqLi5WYmFhqPDExUSdOnCizfXZ2tlwul//hdrurKiqAcFJ40u4EAH4SUlf9TJgwQePGjfMvezyeyi0rE49dfRsAoePQJmnxg1ffLj7x6tsAqBK2FpVatWopMjJSJ0+W/tfLyZMnVbdu3TLbx8TEKCYmpqriSdHVq25fACpf4+4lc1A8xyWVNz3vpzkqDTpUdTIAl2HrqZ/o6GjdeeedWrNmjX/M5/NpzZo1SktLszEZgBtSRKTUe/pPC45LVv603HsaE2kBg9h+efK4ceM0b948LVy4UF9++aVGjRqloqIijRgxwu5oAG5EKf2lQYskZ1LpcWdyyXhKf3tyASiX7XNU/v3f/13ffvutnn32WZ04cUKtW7fWhx9+WGaCLQAETUp/qWlf7kwLhADb76NyPSr9PioAACDoAvn8tv3UDwAAwOVQVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGCsKLsDAFflK5YObZIKT0rxiVKDDlJEpN2pAABVwNYjKuvWrVO/fv2UnJwsh8Oh9957z844MFHecmlGC2nhfdLSkSW/zmhRMg4AuOHZWlSKiop0xx136JVXXrEzBkyVt1x6a6jkOVZ63HO8ZJyyAgA3PFtP/fTp00d9+vSxM0J4OF9kd4LA+Yqllf8pySpnpSXJIX04Xmral9NAAHADC6k5Kl6vV16v17/s8XhsTBNCpibbnaASWJLnaMnclYb32B0GAFBJQuqqn+zsbLlcLv/D7XbbHQl2KzxpdwIAQCUKqSMqEyZM0Lhx4/zLHo+HsnItJh67+jamObRJWvzg1beLT6z8LAAA24RUUYmJiVFMTIzdMUJPdHW7EwSucXfJmVwycbbceSqOkvUNOlR1MgBAFQqpUz8IIxGRUu/pPy04Lln503LvaUykBYAbnK1FpbCwUDt37tTOnTslSd9884127typw4cP2xkLpkjpLw1aJDmTSo87k0vGU/rbkwsAUGUclmWVd1y9Sqxdu1bdunUrMz5s2DAtWLDgqs/3eDxyuVzKz8+X0+mshIQwAnemBYAbSiCf37bOUenatats7EkIFRGRXIIMAGGKOSoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYK8ruAACA6+Arlg5tkgpPSvGJUoMOUkSk3amAoLH1iEp2drZSU1OVkJCgOnXqaMCAAdq3b5+dkQAgdOQtl2a0kBbeJy0dWfLrjBYl48ANwtaikpOTo4yMDG3ZskUfffSRfvzxR/Xs2VNFRUV2xgIA8+Utl94aKnmOlR73HC8Zp6zgBuGwLMuyO8RF3377rerUqaOcnBx17tz5qtt7PB65XC7l5+fL6XRWQULY6jwFFpBUcrrnlbukguOX2cAhOZOlrN2cBoKRAvn8NmqOSn5+viTplltuKXe91+uV1+v1L3s8nirJBUNMTbY7ARAiLMlztGTuSsN77A4DXBdjrvrx+XzKyspSx44d1aJFi3K3yc7Olsvl8j/cbncVpwSAEFJ40u4EwHUz5tTPqFGjtHLlSm3YsEH16tUrd5vyjqi43W5O/YQLTv0AJQ5tkhY/ePXthq3giAqMFHKnfkaPHq0VK1Zo3bp1ly0pkhQTE6OYmJgqTAajRFe3OwFghsbdS+ageI5LKu/fmj/NUWnQoaqTAUFn66kfy7I0evRoLVu2TJ988okaNmxoZxwACA0RkVLv6T8tOC5Z+dNy72lMpMUNwdaikpGRob/97W9asmSJEhISdOLECZ04cULff/+9nbEAwHwp/aVBiyRnUulxZ3LJeEp/e3IBQWbrHBWH49J/CZSYP3++hg8fftXnc3kygLDHnWkRgkJmjooh83gBIHRFRDJhFjc0Yy5PBgAAuBRFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGCsgIvK5MmT5fP5yozn5+dr8ODBQQkFAAAgVaCovPbaa+rUqZO+/vpr/9jatWvVsmVLHThwIKjhAABAeAu4qHz++eeqV6+eWrdurXnz5unpp59Wz549NWTIEG3atKkyMgIAgDAVFegTatSoobfeeksTJ07Uf/zHfygqKkorV65Ujx49KiMfAAAIYxWaTDtz5kz96U9/0uDBg9WoUSONGTNGu3btCnY2AAAQ5gIuKr1799aUKVO0cOFCLV68WJ999pk6d+6su+++Wy+++GJlZAQAAGEq4KJSXFyszz//XA8++KAkKS4uTrNnz9Y777yjP/7xj0EPCAAAwpfDsiwrWC92+vRp1apVK1gvd1Uej0cul0v5+flyOp1Vtl8AAFBxgXx+V2iOyvr16/XII48oLS1NR48elST99a9/1d69eyvycgAAAOUKuKgsXbpUvXr1UlxcnD777DN5vV5JJTd8mzp1atADAgCA8BVwUXn++ef16quvat68eapWrZp/vGPHjtqxY0dQwwEAgPAWcFHZt2+fOnfuXGbc5XLp7NmzwcgEAAAgqQJFpW7dutq/f3+Z8Q0bNqhRo0ZBCQUAACBVoKg89thjGjt2rP7xj3/I4XDo2LFjWrx4sZ566imNGjWqMjICAIAwFfAt9MePHy+fz6cePXro3Llz6ty5s2JiYvTUU08pMzOzMjICAIAwVeH7qJw/f1779+9XYWGhUlJSFB8fH+xsV8V9VAAACD2BfH4HfETloujoaKWkpFT06QAAAFd1TUXll7/85TW/4LvvvlvhMAAAAD93TZNpXS6X/+F0OrVmzRpt377dvz43N1dr1qyRy+WqtKAAACD8XNMRlfnz5/t//t3vfqdBgwbp1VdfVWRkpKSSLyp8/PHHmScCAACCKuDJtLVr19aGDRt0++23lxrft2+fOnTooH/9619BDXglTKYFACD0VOqXEl64cKHcLx/cu3evfD5foC8HAABwWQFf9TNixAiNHDlSBw4c0F133SVJ+sc//qFp06ZpxIgRQQ8IAADCV8BF5Q9/+IPq1q2rl156ScePH5ckJSUl6emnn9aTTz4Z9IAAACB8VfiGb1LJOSZJts0PYY4KAAChp0pu+CbZV1AAAEB4CHgy7cmTJzVkyBAlJycrKipKkZGRpR4AAADBEvARleHDh+vw4cN65plnlJSUJIfDURm5AAAAAi8qGzZs0Pr169W6detKiAMAAPB/Aj7143a7dR3zb0uZPXu2WrVqJafTKafTqbS0NK1cuTIorw0AAEJfwEVlxowZGj9+vA4ePHjdO69Xr56mTZum3Nxcbd++Xd27d9f999+vL7744rpfGwAAhL6AL0+uUaOGzp07pwsXLuimm25StWrVSq0/c+bMdQW65ZZb9N///d8aOXLkVbfl8uTwcu78BbsjAEDYuSn6ui4QLlelXp48Y8aMiua6ouLiYr399tsqKipSWlpaudt4vV55vV7/8sX7uCA8pDy7yu4IABB2Dk7ra+v+Ay4qw4YNC2qA3bt3Ky0tTT/88IPi4+O1bNkypaSklLttdna2pkyZEtT9AwAAc13TqR+Px+M/NHO1oxiBnoI5f/68Dh8+rPz8fL3zzjv6y1/+opycnHLLSnlHVNxuN6d+wgSnfgCg6tl96ueaikpkZKSOHz+uOnXqKCIiotx7p1iWJYfDoeLi4oonl5Senq7GjRtrzpw5V92WOSoAAISeoM9R+eSTT3TLLbdIkj799NPrT3gFPp+v1FETAAAQvq6pqHTp0qXcn6/XhAkT1KdPH9WvX18FBQVasmSJ1q5dq1WrmDQJAACu80sJr9epU6c0dOhQHT9+XC6XS61atdKqVat077332hkLAAAYwtai8tprr9m5ewAAYLiA70wLAABQVSgqAADAWBQVAABgrIDnqLRp06bc+6g4HA7FxsaqSZMmGj58uLp16xaUgAAAIHwFfESld+/e+vrrr1W9enV169ZN3bp1U3x8vA4cOKDU1FQdP35c6enpev/99ysjLwAACCMBH1E5ffq0nnzyST3zzDOlxp9//nkdOnRIq1ev1qRJk/Tcc8/p/vvvD1pQAAAQfq7pFvo/53K5lJubqyZNmpQa379/v+68807l5+dr7969Sk1NVUFBQVDDXopb6AMAEHoC+fwO+NRPbGysNm3aVGZ806ZNio2NlVRyG/yLPwMAAFRUwKd+MjMz9dvf/la5ublKTU2VJG3btk1/+ctfNHHiREnSqlWr1Lp166AGBQAA4SfgUz+StHjxYs2aNUv79u2TJN1+++3KzMzUr3/9a0nS999/778KqDJx6gcAgNATyOd3hYqKKSgqAACEnkA+vyv8XT+5ubn68ssvJUnNmzdXmzZtKvpSAAAA5Qq4qJw6dUoPPfSQ1q5dq5tvvlmSdPbsWXXr1k1vvPGGateuHeyMAAAgTAV81U9mZqYKCgr0xRdf6MyZMzpz5oz27Nkjj8ejMWPGVEZGAAAQpip0H5WPP/7Yf8XPRVu3blXPnj119uzZYOa7IuaoAAAQeir1Pio+n0/VqlUrM16tWjX5fL5AXw4AAOCyAi4q3bt319ixY3Xs2DH/2NGjR/XEE0+oR48eQQ0HAADCW8BFZdasWfJ4PLr11lvVuHFjNW7cWA0bNpTH49HMmTMrIyMAAAhTAV/143a7tWPHDn388cfau3evJKlZs2ZKT08PejgAABDeuOEbAACoUpV2wzefz6cFCxbo3Xff1cGDB+VwONSwYUM9+OCDGjJkiBwOx3UFBwAA+LlrnqNiWZb69++vRx99VEePHlXLli3VvHlzHTp0SMOHD9cDDzxQmTkBAEAYuuYjKgsWLNC6deu0Zs0adevWrdS6Tz75RAMGDNCiRYs0dOjQoIcEAADh6ZqPqLz++uuaOHFimZIilVyyPH78eC1evDio4QAAQHi75qLy+eefq3fv3pdd36dPH+3atSsooQAAAKQAisqZM2eUmJh42fWJiYn67rvvghIKAABACqCoFBcXKyrq8lNaIiMjdeHChaCEAgAAkAKYTGtZloYPH66YmJhy13u93qCFAgAAkAIoKsOGDbvqNlzxAwAAgumai8r8+fMrMwcAAEAZAX8pIQAAQFWhqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYy5iiMm3aNDkcDmVlZdkdBQAAGMKIorJt2zbNmTNHrVq1sjsKAAAwSJTdAQoLC/Xwww9r3rx5ev755+2OA+AanTt/we4IAKrATdH2VgXbi0pGRob69u2r9PT0qxYVr9crr9frX/Z4PJUdD8BlpDy7yu4IAKrAwWl9bd2/rUXljTfe0I4dO7Rt27Zr2j47O1tTpkyp5FQAAMAUthWVI0eOaOzYsfroo48UGxt7Tc+ZMGGCxo0b51/2eDxyu92VFRHAFeT9v152RwAQBhyWZVl27Pi9997TAw88oMjISP9YcXGxHA6HIiIi5PV6S60rj8fjkcvlUn5+vpxOZ2VHBgAAQRDI57dtR1R69Oih3bt3lxobMWKEmjZtqt/97ndXLSkAAODGZ1tRSUhIUIsWLUqNVa9eXTVr1iwzDgAAwpMR91EBAAAoj+2XJ//c2rVr7Y4AAAAMwhEVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjBVldwAAAG4ovmLp0Cap8KQUnyg16CBFRNqdKmTZekRl8uTJcjgcpR5Nmza1MxIAABWXt1ya0UJaeJ+0dGTJrzNalIyjQmw/otK8eXN9/PHH/uWoKNsjAQAQuLzl0ltDJVmlxz3HS8YHLZJS+tsSLZTZ3gqioqJUt25du2MAQGg7X2R3gvDmK5ZW/qfKlBTppzGH9OF4qWlfTgMFyPai8s9//lPJycmKjY1VWlqasrOzVb9+/XK39Xq98nq9/mWPx1NVMQHAbFOT7U6AK7Ikz9GSuSsN77E7TEixdY5K+/bttWDBAn344YeaPXu2vvnmG91zzz0qKCgod/vs7Gy5XC7/w+12V3FiAACuQ+FJuxOEHIdlWeUdp7LF2bNn1aBBA7388ssaOXJkmfXlHVFxu93Kz8+X0+msyqgAYBZO/djr0CZp8YNX327YCo6oqOTz2+VyXdPnt+2nfn7u5ptv1m233ab9+/eXuz4mJkYxMTFVnAoAQkB0dbsThLfG3SVncsnE2XLnqThK1jfoUNXJQp5RN3wrLCzUgQMHlJSUZHcUAACuXUSk1Hv6TwuOS1b+tNx7GhNpK8DWovLUU08pJydHBw8e1KZNm/TAAw8oMjJSgwcPtjMWAACBS+lfcgmy85J/bDuTuTT5Oth66ud///d/NXjwYP3rX/9S7dq11alTJ23ZskW1a9e2MxYAABWT0r/kEmTuTBs0Rk2mDVQgk3EAAIAZAvn8NmqOCgAAwM9RVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADCW7UXl6NGjeuSRR1SzZk3FxcWpZcuW2r59u92xAACAAaLs3Pl3332njh07qlu3blq5cqVq166tf/7zn6pRo4adsQAAgCFsLSrTp0+X2+3W/Pnz/WMNGza0MREAwE7nzl+wOwIucVO0rVXB3qKyfPly9erVSwMHDlROTo7+7d/+TY8//rgee+yxcrf3er3yer3+ZY/HU1VRAQBVIOXZVXZHwCUOTutr6/5tnaPy9ddfa/bs2frFL36hVatWadSoURozZowWLlxY7vbZ2dlyuVz+h9vtruLEAACgKjksy7Ls2nl0dLTatWunTZs2+cfGjBmjbdu2afPmzWW2L++IitvtVn5+vpxOZ5VkBgBUHk79mKcyTv14PB65XK5r+vy29dRPUlKSUlJSSo01a9ZMS5cuLXf7mJgYxcTEVEU0AIAN7J4PAfPYeuqnY8eO2rdvX6mxr776Sg0aNLApEQAAMImtReWJJ57Qli1bNHXqVO3fv19LlizR3LlzlZGRYWcsAABgCFuLSmpqqpYtW6bXX39dLVq00HPPPacZM2bo4YcftjMWAAAwhK2Taa9XIJNxAACAGQL5/Lb9FvoAAACXQ1EBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIwV0l9TefGmuh6Px+YkAADgWl383L6Wm+OHdFEpKCiQJLndbpuTAACAQBUUFMjlcl1xm5D+rh+fz6djx44pISFBDocjqK/t8Xjkdrt15MgRvkfIALwfZuH9MAvvh3l4T67MsiwVFBQoOTlZERFXnoUS0kdUIiIiVK9evUrdh9Pp5D8yg/B+mIX3wyy8H+bhPbm8qx1JuYjJtAAAwFgUFQAAYCyKymXExMRo0qRJiomJsTsKxPthGt4Ps/B+mIf3JHhCejItAAC4sXFEBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUyvHKK6/o1ltvVWxsrNq3b6+tW7faHSlsZWdnKzU1VQkJCapTp44GDBigffv22R0LkqZNmyaHw6GsrCy7o4S1o0eP6pFHHlHNmjUVFxenli1bavv27XbHCkvFxcV65pln1LBhQ8XFxalx48Z67rnnrun7bHB5FJVLvPnmmxo3bpwmTZqkHTt26I477lCvXr106tQpu6OFpZycHGVkZGjLli366KOP9OOPP6pnz54qKiqyO1pY27Ztm+bMmaNWrVrZHSWsfffdd+rYsaOqVaumlStXKi8vTy+99JJq1Khhd7SwNH36dM2ePVuzZs3Sl19+qenTp+vFF1/UzJkz7Y4W0rg8+RLt27dXamqqZs2aJank+4TcbrcyMzM1fvx4m9Ph22+/VZ06dZSTk6POnTvbHScsFRYWqm3btvrzn/+s559/Xq1bt9aMGTPsjhWWxo8fr40bN2r9+vV2R4Gk++67T4mJiXrttdf8Y7/61a8UFxenv/3tbzYmC20cUfmZ8+fPKzc3V+np6f6xiIgIpaena/PmzTYmw0X5+fmSpFtuucXmJOErIyNDffv2LfX/CeyxfPlytWvXTgMHDlSdOnXUpk0bzZs3z+5YYatDhw5as2aNvvrqK0nSrl27tGHDBvXp08fmZKEtpL+UMNhOnz6t4uJiJSYmlhpPTEzU3r17bUqFi3w+n7KystSxY0e1aNHC7jhh6Y033tCOHTu0bds2u6NA0tdff63Zs2dr3LhxmjhxorZt26YxY8YoOjpaw4YNszte2Bk/frw8Ho+aNm2qyMhIFRcX64UXXtDDDz9sd7SQRlFByMjIyNCePXu0YcMGu6OEpSNHjmjs2LH66KOPFBsba3ccqKS8t2vXTlOnTpUktWnTRnv27NGrr75KUbHBW2+9pcWLF2vJkiVq3ry5du7cqaysLCUnJ/N+XAeKys/UqlVLkZGROnnyZKnxkydPqm7dujalgiSNHj1aK1as0Lp161SvXj2744Sl3NxcnTp1Sm3btvWPFRcXa926dZo1a5a8Xq8iIyNtTBh+kpKSlJKSUmqsWbNmWrp0qU2JwtvTTz+t8ePH66GHHpIktWzZUocOHVJ2djZF5TowR+VnoqOjdeedd2rNmjX+MZ/PpzVr1igtLc3GZOHLsiyNHj1ay5Yt0yeffKKGDRvaHSls9ejRQ7t379bOnTv9j3bt2unhhx/Wzp07KSk26NixY5nL9b/66is1aNDApkTh7dy5c4qIKP2xGhkZKZ/PZ1OiGwNHVC4xbtw4DRs2TO3atdNdd92lGTNmqKioSCNGjLA7WljKyMjQkiVL9P777yshIUEnTpyQJLlcLsXFxdmcLrwkJCSUmRtUvXp11axZkzlDNnniiSfUoUMHTZ06VYMGDdLWrVs1d+5czZ071+5oYalfv3564YUXVL9+fTVv3lyfffaZXn75Zf3mN7+xO1pos1DGzJkzrfr161vR0dHWXXfdZW3ZssXuSGFLUrmP+fPn2x0NlmV16dLFGjt2rN0xwtr//M//WC1atLBiYmKspk2bWnPnzrU7UtjyeDzW2LFjrfr161uxsbFWo0aNrN///veW1+u1O1pI4z4qAADAWMxRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBIEkaPny4BgwYYNv+hwwZ4v8WYDtMnjxZrVu3rtBz8/LyVK9ePRUVFQU3FACKChAOHA7HFR+TJ0/Wn/70Jy1YsMCWfLt27dIHH3ygMWPG2LL/65WSkqK7775bL7/8st1RgBsOX0oIhIHjx4/7f37zzTf17LPPlvrW3fj4eMXHx9sRTZI0c+ZMDRw40NYM12vEiBF67LHHNGHCBEVF8VcrECwcUQHCQN26df0Pl8slh8NRaiw+Pr7MqZ+uXbsqMzNTWVlZqlGjhhITEzVv3jz/t4knJCSoSZMmWrlyZal97dmzR3369FF8fLwSExM1ZMgQnT59+rLZiouL9c4776hfv37+sVmzZpX6Rub33ntPDodDr776qn8sPT1d//Vf/+Vffv/999W2bVvFxsaqUaNGmjJlii5cuOBff/bsWT366KOqXbu2nE6nunfvrl27dl0214EDB9SoUSONHj1almXp0KFD6tevn2rUqKHq1aurefPm+uCDD/zb33vvvTpz5oxycnIu+5oAAkdRAXBZCxcuVK1atbR161ZlZmZq1KhRGjhwoDp06KAdO3aoZ8+eGjJkiM6dOyeppAx0795dbdq00fbt2/Xhhx/q5MmTGjRo0GX38fnnnys/P1/t2rXzj3Xp0kV5eXn69ttvJUk5OTmqVauW1q5dK0n68ccftXnzZnXt2lWStH79eg0dOlRjx45VXl6e5syZowULFuiFF17wv+bAgQN16tQprVy5Urm5uWrbtq169OihM2fOlJupU6dO+vWvf61Zs2bJ4XAoIyNDXq9X69at0+7duzV9+vRSR4Cio6PVunVrrV+/vsJ/3gDKYfO3NwOoYvPnz7dcLleZ8WHDhln333+/f7lLly5Wp06d/MsXLlywqlevbg0ZMsQ/dvz4cUuStXnzZsuyLOu5556zevbsWep1jxw5Ykmy9u3bV26eZcuWWZGRkZbP5/OP+Xw+q2bNmtbbb79tWZZltW7d2srOzrbq1q1rWZZlbdiwwapWrZpVVFRkWZZl9ejRw5o6dWqp1/3rX/9qJSUlWZZlWevXr7ecTqf1ww8/lNqmcePG1pw5cyzLsqxJkyZZd9xxh7Vx40arRo0a1h/+8IdS27Zs2dKaPHlyub+Hix544AFr+PDhV9wGQGA4kQrgslq1auX/OTIyUjVr1lTLli39Y4mJiZKkU6dOSSqZFPvpp5+WO9fkwIEDuu2228qMf//994qJiZHD4fCPORwOde7cWWvXrlV6erry8vL0+OOP68UXX9TevXuVk5Oj1NRU3XTTTf79bty4sdQRlOLiYv3www86d+6cdu3apcLCQtWsWbPMvg8cOOBfPnz4sO6991698MILysrKKrXtmDFjNGrUKK1evVrp6en61a9+VerPR5Li4uL8R5cABAdFBcBlVatWrdSyw+EoNXaxXPh8PklSYWGh+vXrp+nTp5d5raSkpHL3UatWLZ07d07nz59XdHS0f7xr166aO3eu1q9frzZt2sjpdPrLS05Ojrp06eLftrCwUFOmTNEvf/nLMq8fGxurwsJCJSUl+U8d/dzNN9/s/7l27dpKTk7W66+/rt/85jdyOp3+dY8++qh69eqlv//971q9erWys7P10ksvKTMz07/NmTNn1Lhx43J/nwAqhjkqAIKmbdu2+uKLL3TrrbeqSZMmpR7Vq1cv9zkX712Sl5dXavziPJW3337bPxela9eu+vjjj7Vx40b/2MX97tu3r8w+mzRpooiICLVt21YnTpxQVFRUmfW1atXyv05cXJxWrFih2NhY9erVSwUFBaUyud1u/fa3v9W7776rJ598UvPmzSu1fs+ePWrTpk0F//QAlIeiAiBoMjIydObMGQ0ePFjbtm3TgQMHtGrVKo0YMULFxcXlPqd27dpq27atNmzYUGq8VatWqlGjhpYsWVKqqLz33nvyer3q2LGjf9tnn31WixYt0pQpU/TFF1/oyy+/1BtvvOG/Kig9PV1paWkaMGCAVq9erYMHD2rTpk36/e9/r+3bt5fab/Xq1fX3v/9dUVFR6tOnjwoLCyVJWVlZWrVqlb755hvt2LFDn376qZo1a+Z/3sGDB3X06FGlp6df958jgP9DUQEQNMnJydq4caOKi4vVs2dPtWzZUllZWbr55psVEXH5v24effRRLV68uNSYw+HQPffcI4fDoU6dOkkqKS9Op1Pt2rUrdYSmV69eWrFihVavXq3U1FTdfffd+uMf/6gGDRr4X+uDDz5Q586dNWLECN1222166KGHdOjQIf88m5+Lj4/XypUrZVmW+vbtq6KiIhUXFysjI0PNmjVT7969ddttt+nPf/6z/zmvv/66evbs6d8ngOBwWJZl2R0CQHj7/vvvdfvtt+vNN99UWlqa3XECdv78ef3iF7/QkiVLSh3pAXD9OKICwHZxcXFatGjRFW8MZ7LDhw9r4sSJlBSgEnBEBQAAGIsjKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWP8fzQp4rkpqe3cAAAAASUVORK5CYII=\n" + }, + "metadata": {} + } + ], + "source": [ + "plot_lifelines(obs)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9v-FAT6MTTlz" + }, + "source": [ + "Each line represents the time a dog spends at the shelter. Each dot represents an adoption.\n", + "We can see, for example:\n", + "\n", + "* The dog with index 0 arrived during week 0, and was adopted during week 5.\n", + "\n", + "* The dog with index 3 arrived during week 2, and had not been adopted at the end of week 9.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "gP-olgxWTTlz" + }, + "source": [ + "## Estimating survival\n", + "\n", + "Now suppose we want to know the distribution of \"survival time\" from arrival to adoption.\n", + "For the dogs that were adopted, we have all the data we need. \n", + "For the others, we have only partial information: if a dog hasn't been adopted yet, we don't know when it will be, but we can put a lower bound on it.\n", + "\n", + "When we have a mixture of complete and incomplete observations -- adopted and unadopted dogs -- we can't compute the Survival function directly.\n", + "Instead, we have to work backwards: we estimate the hazard function first, then use it to compute the survival function, CDF, and PMF.\n", + "\n", + "Specifically, we'll use Kaplan-Meier estimation, which is based on two key ideas.\n", + "\n", + "The first idea is that we can ignore the arrival time in the observed data, and consider only the durations. In effect, we can take the actual lifelines and shift them so they all start at 0, like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "tags": [], + "id": "X4RK8YAHTTlz" + }, + "outputs": [], + "source": [ + "duration = obs['end'] - obs['start']" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "tags": [], + "id": "pFB7dzT2TTl0" + }, + "outputs": [], + "source": [ + "shifted = obs.copy()\n", + "shifted['start'] = 0\n", + "shifted['end'] = duration" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 449 + }, + "id": "hUSE7kJYTTl0", + "outputId": "61b1c0af-b72a-4001-ae1d-8903025c010b" + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAK6tJREFUeJzt3Xl0FGWixuG3k5BOhHSzh0RCWEUChEWWCShrFBABmYMoIgQGPY6DEEG8AkdFHCU4d7wyowwCOoAKoiPreG/YSYABlX0RjIIsDgTQEdIJOg2m6/4R7bElLIGk64v9e87pA/V1ddXbLdIvVV9VHJZlWQIAADBQmN0BAAAALoWiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgrAi7A1wPn8+nEydOKCYmRg6Hw+44AADgKliWpfz8fMXHxyss7PLHTMp1UTlx4oQSEhLsjgEAAK7Bl19+qdq1a192nXJdVGJiYiQVvVGXy2VzGgAAcDU8Ho8SEhL83+OXU66Lyo+ne1wuF0UFAIBy5mqmbTCZFgAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYq1zfmRYAUEZ8hdLRzVLBKalSrJTYQQoLtzsVQpARR1SmT5+uunXrKioqSu3bt9fHH39sdyQACF37l0vTmknz7pIWjSj6dVqzonEgyGwvKu+++67Gjh2rSZMmaceOHWrRooV69Oih06dP2x0NAELP/uXSe0Mlz4nAcU9u0ThlBUHmsCzLsjNA+/bt1bZtW7366quSJJ/Pp4SEBI0aNUrjx4+/7Gs9Ho/cbrfy8vLK5ocSnj9X+tsEAFP5CqXp7aT83Eus4JBc8dJjezkNhOtSku9vW+eonD9/Xtu3b9eECRP8Y2FhYUpNTdWWLVsuWt/r9crr9fqXPR5P2QacEl+22weAcsWSPMeL5q7Uu83uMAgRtp76+frrr1VYWKjY2NiA8djYWJ08efKi9TMyMuR2u/2PhISEYEUFAPyo4JTdCRBCytVVPxMmTNDYsWP9yx6Pp2zLysQTV14HAH4pjm6W5g+48nqVYq+8DlBKbC0q1atXV3h4uE6dCmznp06dUq1atS5a3+l0yul0BiueFFkxePsCALs16FY0B8WTK6m46Ys/zFFJ7BDsZAhhtp76iYyM1C233KK1a9f6x3w+n9auXauUlBQbkwFACAoLl3q++MOC42dP/rDccyoTaRFUtl+ePHbsWM2ePVvz5s3TgQMH9Mgjj+jcuXMaPny43dEAIPQk9ZUGvim54gLHXfFF40l97cmFkGX7HJV7771XX331lZ555hmdPHlSLVu21IoVKy6aYAsACJKkvtLNvbkzLYxg+31UrkeZ30cFAACUupJ8f9t+6gcAAOBSKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwVoTdAWAgX6F0dLNUcEqqFCsldpDCwu1OBQAIQbYeUdmwYYP69Omj+Ph4ORwOLV261M44kKT9y6VpzaR5d0mLRhT9Oq1Z0TgAAEFma1E5d+6cWrRooenTp9sZAz/av1x6b6jkORE47sktGqesAACCzNZTP7169VKvXr3sjHB558/ZnSB4fIVS5n9Jsop50pLkkFaMl27uzWkgAEDQlKs5Kl6vV16v17/s8XjKdodT4st2++WKJXmOF81dqXeb3WEAACGiXF31k5GRIbfb7X8kJCTYHSn0FJyyOwEAIISUqyMqEyZM0NixY/3LHo+nbMvKxBNXXueX4uhmaf6AK69XKbbsswAA8INyVVScTqecTmfwdhhZMXj7sluDbpIrvmjibLHzVBxFzyd2CHYyAEAIK1enflCGwsKlni/+sOD42ZM/LPecykRaAEBQ2VpUCgoKtGvXLu3atUuSdPjwYe3atUvHjh2zM1boSuorDXxTcsUFjrvii8aT+tqTCwAQshyWZRV3nD8osrKy1LVr14vG09LSNHfu3Cu+3uPxyO12Ky8vTy6XqwwShijuTAsAKEMl+f62dY5Kly5dZGNPwqWEhXMJMgDACMxRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgRdgcAAOP4CqWjm6WCU1KlWCmxgxQWbncqICTZekQlIyNDbdu2VUxMjGrWrKm7775bOTk5dkYCEOr2L5emNZPm3SUtGlH067RmReMAgs7WopKdna2RI0fqww8/1OrVq3XhwgXdcccdOnfunJ2xAISq/cul94ZKnhOB457conHKChB0DsuyLLtD/Oirr75SzZo1lZ2drU6dOl1xfY/HI7fbrby8PLlcrtIPdJ7CBIQMX6E0vZ2Un3uJFRySK156bC+ngYDrVJLvb6PmqOTl5UmSqlatWuzzXq9XXq/Xv+zxeMo20JT4st0+gHLEkjzHi+au1LvN7jBAyDDmqh+fz6fHHntMHTt2VLNmzYpdJyMjQ2632/9ISEgIckoAIa/glN0JgJBizKmfRx55RJmZmdq0aZNq165d7DrFHVFJSEjg1A+A63d0szR/wJXXS/uAIyrAdSp3p34effRRffDBB9qwYcMlS4okOZ1OOZ3O4AWLrBi8fQGwV4NuRXNQPLmSivv32w9zVBI7BDsZENJsPfVjWZYeffRRLVmyROvWrVO9evXsjAMglIWFSz1f/GHB8bMnf1juOZWJtECQ2VpURo4cqbffflsLFixQTEyMTp48qZMnT+q7776zMxaAUJXUVxr4puSKCxx3xReNJ/W1JxcQwmydo+Jw/PxfLUXmzJmjYcOGXfH1ZX55MoDQxJ1pgTJVbuaoGDKPFwAChYUzYRYwhDGXJwMAAPwcRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgrBIXlWeffVY+n++i8by8PA0aNKhUQgEAAEjXUFTeeOMN3Xrrrfriiy/8Y1lZWWrevLkOHTpUquEAAEBoK3FR2bNnj2rXrq2WLVtq9uzZeuKJJ3THHXdoyJAh2rx5c1lkBAAAISqipC+oUqWK3nvvPU2cOFEPP/ywIiIilJmZqe7du5dFPgAAEMKuaTLtK6+8oj/96U8aNGiQ6tevr9GjR2v37t2lnQ0AAIS4EheVnj17avLkyZo3b57mz5+vnTt3qlOnTvrVr36lP/zhD2WREQAAhKgSF5XCwkLt2bNHAwYMkCRFR0drxowZev/99/Xyyy+XekAAABC6HJZlWaW1sa+//lrVq1cvrc1dkcfjkdvtVl5enlwuV9D2CwAArl1Jvr+vaY7Kxo0b9cADDyglJUXHjx+XJL311lv69NNPr2VzAAAAxSpxUVm0aJF69Oih6Oho7dy5U16vV1LRDd+mTJlS6gEBAEDoKnFRef755/Xaa69p9uzZqlChgn+8Y8eO2rFjR6mGAwAAoa3ERSUnJ0edOnW6aNztduvs2bOlkQkAAEDSNRSVWrVq6eDBgxeNb9q0SfXr1y+VUAAAANI1FJWHHnpI6enp+uijj+RwOHTixAnNnz9f48aN0yOPPFIWGQEAQIgq8S30x48fL5/Pp+7du+vbb79Vp06d5HQ6NW7cOI0aNaosMgIAgBB1zfdROX/+vA4ePKiCggIlJSWpUqVKpZ3tiriPCgAA5U9Jvr9LfETlR5GRkUpKSrrWlwMAAFzRVRWVX//611e9wcWLF19zGAAAgJ+6qsm0brfb/3C5XFq7dq22bdvmf3779u1au3at3G53mQUFAACh56qOqMyZM8f/+yeffFIDBw7Ua6+9pvDwcElFP6jwd7/7HfNEAABAqSrxZNoaNWpo06ZNaty4ccB4Tk6OOnTooH/961+lGvBymEwLAED5U6Y/lPD7778v9ocPfvrpp/L5fCXdHAAAwCWV+Kqf4cOHa8SIETp06JDatWsnSfroo480depUDR8+vNQDAgCA0FXiovLHP/5RtWrV0ksvvaTc3FxJUlxcnJ544gk9/vjjpR4QAACErmu+4ZtUdI5Jkm3zQ5ijAgBA+ROUG75J9hUUAAAQGko8mfbUqVMaMmSI4uPjFRERofDw8IAHAABAaSnxEZVhw4bp2LFjevrppxUXFyeHw1EWuQAAAEpeVDZt2qSNGzeqZcuWZRAHAADgP0p86ichIUHXMf82wIwZM5ScnCyXyyWXy6WUlBRlZmaWyrYBAED5V+KiMm3aNI0fP15Hjhy57p3Xrl1bU6dO1fbt27Vt2zZ169ZN/fr10yeffHLd2wYAAOVfiS9PrlKlir799lt9//33uuGGG1ShQoWA57/55pvrClS1alX993//t0aMGHHFdcv68uRvz39f6tsEAKA8uSHyui4QLlaZXp48bdq0a811WYWFhfrb3/6mc+fOKSUlpdh1vF6vvF6vf/nH+7iUlaRnVpbp9gEAMN2Rqb1t3X+Ji0paWlqpBti7d69SUlL073//W5UqVdKSJUuUlJRU7LoZGRmaPHlyqe4fAACY66pO/Xg8Hv+hmSsdxSjpKZjz58/r2LFjysvL0/vvv6/XX39d2dnZxZaV4o6oJCQkcOoHAIAyYvepn6sqKuHh4crNzVXNmjUVFhZW7L1TLMuSw+FQYWHhtSeXlJqaqgYNGmjmzJlXXJdb6AMAUP6U+hyVdevWqWrVqpKk9evXX3/Cy/D5fAFHTQAAQOi6qqLSuXPnYn9/vSZMmKBevXqpTp06ys/P14IFC5SVlaWVK5nECgAArvOHEl6v06dPa+jQocrNzZXb7VZycrJWrlyp22+/3c5YAADAELYWlTfeeMPO3QMAAMOV+M60AAAAwUJRAQAAxqKoAAAAY5V4jkqrVq2KvY+Kw+FQVFSUGjZsqGHDhqlr166lEhAAAISuEh9R6dmzp7744gtVrFhRXbt2VdeuXVWpUiUdOnRIbdu2VW5urlJTU7Vs2bKyyAsAAEJIiY+ofP3113r88cf19NNPB4w///zzOnr0qFatWqVJkybp97//vfr161dqQQEAQOi5qlvo/5Tb7db27dvVsGHDgPGDBw/qlltuUV5enj799FO1bdtW+fn5pRr257iFPgAA5U9Jvr9LfOonKipKmzdvvmh88+bNioqKklR0G/wffw8AAHCtSnzqZ9SoUfrtb3+r7du3q23btpKkrVu36vXXX9fEiRMlSStXrlTLli1LNSgAAAg9JT71I0nz58/Xq6++qpycHElS48aNNWrUKN1///2SpO+++85/FVBZ4tQPAADlT0m+v6+pqJiCogIAQPlTku/va/5ZP9u3b9eBAwckSU2bNlWrVq2udVMAAADFKnFROX36tO677z5lZWWpcuXKkqSzZ8+qa9euWrhwoWrUqFHaGQEAQIgq8VU/o0aNUn5+vj755BN98803+uabb7Rv3z55PB6NHj26LDICAIAQdU33UVmzZo3/ip8fffzxx7rjjjt09uzZ0sx3WcxRAQCg/CnT+6j4fD5VqFDhovEKFSrI5/OVdHMAAACXVOKi0q1bN6Wnp+vEiRP+sePHj2vMmDHq3r17qYYDAAChrcRF5dVXX5XH41HdunXVoEEDNWjQQPXq1ZPH49Err7xSFhkBAECIKvFVPwkJCdqxY4fWrFmjTz/9VJLUpEkTpaamlno4AAAQ2rjhGwAACKoyu+Gbz+fT3LlztXjxYh05ckQOh0P16tXTgAEDNGTIEDkcjusKDgAA8FNXPUfFsiz17dtXDz74oI4fP67mzZuradOmOnr0qIYNG6b+/fuXZU4AABCCrvqIyty5c7VhwwatXbtWXbt2DXhu3bp1uvvuu/Xmm29q6NChpR4SAACEpqs+ovLOO+9o4sSJF5UUqeiS5fHjx2v+/PmlGg4AAIS2qy4qe/bsUc+ePS/5fK9evbR79+5SCQUAACCVoKh88803io2NveTzsbGxOnPmTKmEAgAAkEpQVAoLCxURcekpLeHh4fr+++9LJRQAAIBUgsm0lmVp2LBhcjqdxT7v9XpLLRQAAIBUgqKSlpZ2xXW44gcAAJSmqy4qc+bMKcscAAAAFynxDyUEAAAIFooKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLGMKSpTp06Vw+HQY489ZncUAABgCCOKytatWzVz5kwlJyfbHQUAABgkwu4ABQUFGjx4sGbPnq3nn3/e7jgBvj3/vd0RAMBWN0Ta/jWBEGf7n8CRI0eqd+/eSk1NvWJR8Xq98nq9/mWPx1Om2ZKeWVmm2wcA0x2Z2tvuCAhxthaVhQsXaseOHdq6detVrZ+RkaHJkyeXcSoAAGAK24rKl19+qfT0dK1evVpRUVFX9ZoJEyZo7Nix/mWPx6OEhISyiqj9z/Uos20DAIArc1iWZdmx46VLl6p///4KDw/3jxUWFsrhcCgsLExerzfgueJ4PB653W7l5eXJ5XKVdWQAAFAKSvL9bdsRle7du2vv3r0BY8OHD9fNN9+sJ5988oolBQAA/PLZVlRiYmLUrFmzgLGKFSuqWrVqF40DAIDQZMR9VAAAAIpj++XJP5WVlWV3BAAAYBCOqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGCsCLsDAEbyFUpHN0sFp6RKsVJiByks3O5UABBybD2i8uyzz8rhcAQ8br75ZjsjAdL+5dK0ZtK8u6RFI4p+ndasaBwAEFS2H1Fp2rSp1qxZ41+OiLA9EkLZ/uXSe0MlWYHjntyi8YFvSkl9bYkGAKHI9lYQERGhWrVq2R2jeOfP2Z0AweQrlDL/SxeVFOmHMYe0Yrx0c29OAwFAkNheVD7//HPFx8crKipKKSkpysjIUJ06dYpd1+v1yuv1+pc9Hk/ZhpsSX7bbRzljSZ7jRXNX6t1mdxgACAm2zlFp37695s6dqxUrVmjGjBk6fPiwbrvtNuXn5xe7fkZGhtxut/+RkJAQ5MSAiibYAgCCwmFZVnHHuW1x9uxZJSYm6n/+5380YsSIi54v7ohKQkKC8vLy5HK5Sj8Qp35Cy9HN0vwBV14v7QOOqADAdfB4PHK73Vf1/W37qZ+fqly5sm666SYdPHiw2OedTqecTmfwAkVWDN6+YL8G3SRXfNHE2WLnqTiKnk/sEOxkABCyjLrhW0FBgQ4dOqS4uDi7oyAUhYVLPV/8YcHxsyd/WO45lYm0ABBEthaVcePGKTs7W0eOHNHmzZvVv39/hYeHa9CgQXbGQihL6lt0CbLrZ2XZFc+lyQBgA1tP/fzzn//UoEGD9K9//Us1atTQrbfeqg8//FA1atSwMxZCXVLfokuQuTMtANjOqMm0JVWSyTgAAMAMJfn+NmqOCgAAwE9RVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADAWRQUAABiLogIAAIxFUQEAAMaiqAAAAGNRVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBYFBUAAGAsigoAADCW7UXl+PHjeuCBB1StWjVFR0erefPm2rZtm92xAACAASLs3PmZM2fUsWNHde3aVZmZmapRo4Y+//xzValSxc5YAADAELYWlRdffFEJCQmaM2eOf6xevXo2Jgr07fnv7Y4Am90Qaev/IgAQ8mz9W3j58uXq0aOH7rnnHmVnZ+vGG2/U7373Oz300EPFru/1euX1ev3LHo+nTPMlPbOyTLcP8x2Z2tvuCAAQ0mydo/LFF19oxowZatSokVauXKlHHnlEo0eP1rx584pdPyMjQ2632/9ISEgIcmIAABBMDsuyLLt2HhkZqTZt2mjz5s3+sdGjR2vr1q3asmXLResXd0QlISFBeXl5crlcpZ6PUz/g1A8AlD6PxyO3231V39+2/i0cFxenpKSkgLEmTZpo0aJFxa7vdDrldDqDEU0SX1IAANjN1lM/HTt2VE5OTsDYZ599psTERJsSAQAAk9haVMaMGaMPP/xQU6ZM0cGDB7VgwQLNmjVLI0eOtDMWAAAwhK1FpW3btlqyZIneeecdNWvWTL///e81bdo0DR482M5YAADAELZOpr1eJZmMAwAAzFCS72/bb6EPAABwKRQVAABgLIoKAAAwFkUFAAAYi6ICAACMRVEBAADGoqgAAABjUVQAAICxKCoAAMBY5frHA/94U12Px2NzEgAAcLV+/N6+mpvjl+uikp+fL0lKSEiwOQkAACip/Px8ud3uy65Trn/Wj8/n04kTJxQTEyOHw1Gq2/Z4PEpISNCXX34Zkj9HKNTfv8RnwPsP7fcv8RmE+vuXyu4zsCxL+fn5io+PV1jY5WehlOsjKmFhYapdu3aZ7sPlcoXsH1CJ9y/xGfD+Q/v9S3wGof7+pbL5DK50JOVHTKYFAADGoqgAAABjUVQuwel0atKkSXI6nXZHsUWov3+Jz4D3H9rvX+IzCPX3L5nxGZTrybQAAOCXjSMqAADAWBQVAABgLIoKAAAwFkUFAAAYi6JSjOnTp6tu3bqKiopS+/bt9fHHH9sdKWg2bNigPn36KD4+Xg6HQ0uXLrU7UlBlZGSobdu2iomJUc2aNXX33XcrJyfH7lhBNWPGDCUnJ/tv8JSSkqLMzEy7Y9lm6tSpcjgceuyxx+yOEhTPPvusHA5HwOPmm2+2O1bQHT9+XA888ICqVaum6OhoNW/eXNu2bbM7VlDUrVv3oj8DDodDI0eOtCUPReVn3n33XY0dO1aTJk3Sjh071KJFC/Xo0UOnT5+2O1pQnDt3Ti1atND06dPtjmKL7OxsjRw5Uh9++KFWr16tCxcu6I477tC5c+fsjhY0tWvX1tSpU7V9+3Zt27ZN3bp1U79+/fTJJ5/YHS3otm7dqpkzZyo5OdnuKEHVtGlT5ebm+h+bNm2yO1JQnTlzRh07dlSFChWUmZmp/fv366WXXlKVKlXsjhYUW7duDfjvv3r1aknSPffcY08gCwHatWtnjRw50r9cWFhoxcfHWxkZGTamsocka8mSJXbHsNXp06ctSVZ2drbdUWxVpUoV6/XXX7c7RlDl5+dbjRo1slavXm117tzZSk9PtztSUEyaNMlq0aKF3TFs9eSTT1q33nqr3TGMkZ6ebjVo0MDy+Xy27J8jKj9x/vx5bd++Xampqf6xsLAwpaamasuWLTYmg13y8vIkSVWrVrU5iT0KCwu1cOFCnTt3TikpKXbHCaqRI0eqd+/eAX8fhIrPP/9c8fHxql+/vgYPHqxjx47ZHSmoli9frjZt2uiee+5RzZo11apVK82ePdvuWLY4f/683n77bf3mN78p9R/+e7UoKj/x9ddfq7CwULGxsQHjsbGxOnnypE2pYBefz6fHHntMHTt2VLNmzeyOE1R79+5VpUqV5HQ69dvf/lZLlixRUlKS3bGCZuHChdqxY4cyMjLsjhJ07du319y5c7VixQrNmDFDhw8f1m233ab8/Hy7owXNF198oRkzZqhRo0ZauXKlHnnkEY0ePVrz5s2zO1rQLV26VGfPntWwYcNsy1Cuf3oyUJZGjhypffv2hdz5eUlq3Lixdu3apby8PL3//vtKS0tTdnZ2SJSVL7/8Uunp6Vq9erWioqLsjhN0vXr18v8+OTlZ7du3V2Jiot577z2NGDHCxmTB4/P51KZNG02ZMkWS1KpVK+3bt0+vvfaa0tLSbE4XXG+88YZ69eql+Ph42zJwROUnqlevrvDwcJ06dSpg/NSpU6pVq5ZNqWCHRx99VB988IHWr1+v2rVr2x0n6CIjI9WwYUPdcsstysjIUIsWLfSnP/3J7lhBsX37dp0+fVqtW7dWRESEIiIilJ2drT//+c+KiIhQYWGh3RGDqnLlyrrpppt08OBBu6METVxc3EWlvEmTJiF3Cuzo0aNas2aNHnzwQVtzUFR+IjIyUrfccovWrl3rH/P5fFq7dm3InZ8PVZZl6dFHH9WSJUu0bt061atXz+5IRvD5fPJ6vXbHCIru3btr79692rVrl//Rpk0bDR48WLt27VJ4eLjdEYOqoKBAhw4dUlxcnN1RgqZjx44X3Zbgs88+U2Jiok2J7DFnzhzVrFlTvXv3tjUHp35+ZuzYsUpLS1ObNm3Url07TZs2TefOndPw4cPtjhYUBQUFAf9yOnz4sHbt2qWqVauqTp06NiYLjpEjR2rBggVatmyZYmJi/HOT3G63oqOjbU4XHBMmTFCvXr1Up04d5efna8GCBcrKytLKlSvtjhYUMTExF81JqlixoqpVqxYSc5XGjRunPn36KDExUSdOnNCkSZMUHh6uQYMG2R0taMaMGaMOHTpoypQpGjhwoD7++GPNmjVLs2bNsjta0Ph8Ps2ZM0dpaWmKiLC5KthyrZHhXnnlFatOnTpWZGSk1a5dO+vDDz+0O1LQrF+/3pJ00SMtLc3uaEFR3HuXZM2ZM8fuaEHzm9/8xkpMTLQiIyOtGjVqWN27d7dWrVpldyxbhdLlyffee68VFxdnRUZGWjfeeKN17733WgcPHrQ7VtD9/e9/t5o1a2Y5nU7r5ptvtmbNmmV3pKBauXKlJcnKycmxO4rlsCzLsqciAQAAXB5zVAAAgLEoKgAAwFgUFQAAYCyKCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAJSauXPnqnLlykHZV05OjmrVqqX8/Pyg7K84DodDS5cuvabX3nfffXrppZdKNxDwC0RRAcqZYcOGyeFwyOFwqEKFCoqNjdXtt9+uv/71r/L5fEHLUbduXU2bNi1g7N5779Vnn30WlP1PmDBBo0aNUkxMTFD2V9qeeuopvfDCC8rLy7M7CmA0igpQDvXs2VO5ubk6cuSIMjMz1bVrV6Wnp+uuu+7S999/f83btSzrul4fHR2tmjVrXvPrr9axY8f0wQcfaNiwYWW+r7LSrFkzNWjQQG+//bbdUQCjUVSAcsjpdKpWrVq68cYb1bp1a02cOFHLli1TZmam5s6dK0k6cuSIHA6Hdu3a5X/d2bNn5XA4lJWVJUnKysqSw+FQZmambrnlFjmdTm3atEmHDh1Sv379FBsbq0qVKqlt27Zas2aNfztdunTR0aNHNWbMGP/RHan4Uz8zZsxQgwYNFBkZqcaNG+utt94KeN7hcOj1119X//79dcMNN6hRo0Zavnz5Zd//e++9pxYtWujGG2+UVFSwatSooffff9+/TsuWLRUXF+df3rRpk5xOp7799lv/Z/Hggw+qRo0acrlc6tatm3bv3h2wn2XLlql169aKiopS/fr1NXny5MsWuUmTJikuLk579uyRJP3lL39Ro0aNFBUVpdjYWA0YMCBg/T59+mjhwoWXfa9AqKOoAL8Q3bp1U4sWLbR48eISv3b8+PGaOnWqDhw4oOTkZBUUFOjOO+/U2rVrtXPnTvXs2VN9+vTRsWPHJEmLFy9W7dq19dxzzyk3N1e5ubnFbnfJkiVKT0/X448/rn379unhhx/W8OHDtX79+oD1Jk+erIEDB2rPnj268847NXjwYH3zzTeXzLtx40a1adPGv+xwONSpUyd/ATtz5owOHDig7777Tp9++qkkKTs7W23bttUNN9wgSbrnnnt0+vRpZWZmavv27WrdurW6d+/u3+/GjRs1dOhQpaena//+/Zo5c6bmzp2rF1544aI8lmVp1KhRevPNN7Vx40YlJydr27ZtGj16tJ577jnl5ORoxYoV6tSpU8Dr2rVrp48//lher/dy/3mA0Gbrz24GUGJpaWlWv379in3u3nvvtZo0aWJZlmUdPnzYkmTt3LnT//yZM2csSdb69esty7Ks9evXW5KspUuXXnG/TZs2tV555RX/cmJiovXyyy8HrDNnzhzL7Xb7lzt06GA99NBDAevcc8891p133ulflmQ99dRT/uWCggJLkpWZmXnJLC1atLCee+65gLE///nPVtOmTS3LsqylS5da7du3t/r162fNmDHDsizLSk1NtSZOnGhZlmVt3LjRcrlc1r///e+AbTRo0MCaOXOmZVmW1b17d2vKlCkBz7/11ltWXFxcQPa//e1v1v333281adLE+uc//+l/btGiRZbL5bI8Hs8l38fu3bstSdaRI0cuuQ4Q6jiiAvyCWJblPw1TEj89OiFJBQUFGjdunJo0aaLKlSurUqVKOnDggP+IytU6cOCAOnbsGDDWsWNHHThwIGAsOTnZ//uKFSvK5XLp9OnTl9zud999p6ioqICxzp07a//+/frqq6+UnZ2tLl26qEuXLsrKytKFCxe0efNmdenSRZK0e/duFRQUqFq1aqpUqZL/cfjwYR06dMi/znPPPRfw/EMPPaTc3Fz/6SNJGjNmjD766CNt2LDBfypKkm6//XYlJiaqfv36GjJkiObPnx/wOqloTo+ki8YB/EeE3QEAlJ4DBw6oXr16kqSwsKJ/h1iW5X/+woULxb6uYsWKAcvjxo3T6tWr9cc//lENGzZUdHS0BgwYoPPnz5dJ7goVKgQsOxyOy17BVL16dZ05cyZgrHnz5qpataqys7OVnZ2tF154QbVq1dKLL76orVu36sKFC+rQoYOkoiIWFxfnP1X0Uz/OsSkoKNDkyZP161//+qJ1flqSbr/9dr3zzjtauXKlBg8e7B+PiYnRjh07lJWVpVWrVumZZ57Rs88+q61bt/r38eNppho1alz6wwFCHEUF+IVYt26d9u7dqzFjxkj6z5dfbm6uWrVqJUkBE2sv5x//+IeGDRum/v37Syr60j5y5EjAOpGRkSosLLzsdpo0aaJ//OMfSktLC9h2UlLSVeW4lFatWmn//v0BYw6HQ7fddpuWLVumTz75RLfeeqtuuOEGeb1ezZw5U23atPEXstatW+vkyZOKiIhQ3bp1i91H69atlZOTo4YNG142S9++fdWnTx/df//9Cg8P13333ed/LiIiQqmpqUpNTdWkSZNUuXJlrVu3zl9+9u3bp9q1a6t69erX8WkAv2wUFaAc8nq9OnnypAoLC3Xq1CmtWLFCGRkZuuuuuzR06FBJRacVfvWrX2nq1KmqV6+eTp8+raeeeuqqtt+oUSMtXrxYffr0kcPh0NNPP33REY66detqw4YNuu++++R0Oov9sn3iiSc0cOBAtWrVSqmpqfr73/+uxYsXB1xBdC169OihBx98UIWFhQoPD/ePd+nSRY8//rjatGmjSpUqSZI6deqk+fPn64knnvCvl5qaqpSUFN199936wx/+oJtuukknTpzQ//7v/6p///5q06aNnnnmGd11112qU6eOBgwYoLCwMO3evVv79u3T888/H5Cnf//+euuttzRkyBBFRERowIAB+uCDD/TFF1+oU6dOqlKliv7v//5PPp9PjRs39r9u48aNuuOOO67rswB+8eyeJAOgZNLS0ixJliQrIiLCqlGjhpWammr99a9/tQoLCwPW3b9/v5WSkmJFR0dbLVu2tFatWlXsZNozZ84EvO7w4cNW165drejoaCshIcF69dVXrc6dO1vp6en+dbZs2WIlJydbTqfT+vGvkp9PprUsy/rLX/5i1a9f36pQoYJ10003WW+++WbA85KsJUuWBIy53W5rzpw5l/wMLly4YMXHx1srVqwIGN+5c6clyXryySf9Yy+//LIl6aJ1PR6PNWrUKCs+Pt6qUKGClZCQYA0ePNg6duyYf50VK1ZYHTp0sKKjoy2Xy2W1a9fOmjVr1iWzv/vuu1ZUVJS1aNEia+PGjVbnzp2tKlWqWNHR0VZycrL17rvv+tf97rvvLLfbbW3ZsuWS7xOAZTks6ycnsAGgnJg+fbqWL1+ulStX2h3lmsyYMUNLlizRqlWr7I4CGI1TPwDKpYcfflhnz55Vfn5+ubyNfoUKFfTKK6/YHQMwHkdUAACAsbiPCgAAMBZFBQAAGIuiAgAAjEVRAQAAxqKoAAAAY1FUAACAsSgqAADAWBQVAABgLIoKAAAw1v8Dfsl01SDWp4kAAAAASUVORK5CYII=\n" + }, + "metadata": {} + } + ], + "source": [ + "plot_lifelines(shifted)\n", + "plt.xlabel('Duration (weeks)');" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "wYv8QX6_TTl0" + }, + "source": [ + "Notice that the x-axis in this figure is duration, not time.\n", + "\n", + "The second key idea is that we can estimate the hazard function by considering:\n", + "\n", + "* The number of dogs adopted at each duration, divided by\n", + "\n", + "* The number of dogs \"at risk\" at each duration, where \"at risk\" means that they *could* be adopted.\n", + "\n", + "For example:\n", + "\n", + "* At duration 1, there is 1 adoption out of 7 dogs at risk, so the hazard rate is `1/7`.\n", + "\n", + "* At duration 2, there is 1 adoption out of 6 dogs at risk, so the hazard rate is `1/6`.\n", + "\n", + "* At duration 4, there is 1 adoption out of 4 dogs at risk, so the hazard rate is `1/4`.\n", + "\n", + "And so on. Now let's see how that works computationally." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RoNynN_uTTl0" + }, + "source": [ + "## Computing \"at risk\"\n", + "\n", + "For each observed duration, we would like to compute the number of dogs that were at risk.\n", + "Here are the unique durations, in order:" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "9ElKOYRnTTl0", + "outputId": "0ffb2ed7-46ee-479d-c03d-5b36ac693ac5" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "array([1, 2, 4, 5, 7])" + ] + }, + "metadata": {}, + "execution_count": 30 + } + ], + "source": [ + "ts = duration.unique()\n", + "ts.sort()\n", + "ts" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "OM-V9gg5TTl0" + }, + "source": [ + "To compute the number of dogs at risk, we can loop through `ts` and count the number of dogs where `t` is less than or equal to `end`." + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "scrolled": true, + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "3e3xfHhLTTl0", + "outputId": "7b1b6c08-c968-4c43-ae60-dbc6e6f38568" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "1 7\n", + "2 6\n", + "4 4\n", + "5 3\n", + "7 1\n", + "dtype: int64" + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
0
17
26
44
53
71
\n", + "

" + ] + }, + "metadata": {}, + "execution_count": 31 + } + ], + "source": [ + "at_risk = pd.Series(0, index=ts)\n", + "\n", + "for t in ts:\n", + " k = (t <= shifted['end'])\n", + " at_risk[t] = k.sum()\n", + "\n", + "at_risk" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "LRF-nVAyTTl0" + }, + "source": [ + "If you don't like mixing for loops with array operations, we can do the same computation using mesh grids." + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "IPbcUPxETTl1", + "outputId": "cf15369c-d2f0-4c70-e00f-02c744aab826" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "(5, 7)" + ] + }, + "metadata": {}, + "execution_count": 32 + } + ], + "source": [ + "E, T = np.meshgrid(shifted['end'], ts)\n", + "T.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "B-7MVJfmTTl1" + }, + "source": [ + "The results are arrays with one row for each value of `t` and one column for each dog.\n", + "Now we can use comparison operators to compare all values of `t` to all values of `end` at the same time." + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "DaILrZrXTTl1", + "outputId": "d42b5a4d-129c-4a50-aadd-a85204101369" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "array([7, 6, 4, 3, 1])" + ] + }, + "metadata": {}, + "execution_count": 33 + } + ], + "source": [ + "at_risk = (T <= E).sum(axis=1)\n", + "at_risk" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "LyWqWdv0TTl1" + }, + "source": [ + "The result is an array with the number of dogs at risk for each value of `t`." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "lg5zAb0hTTl1" + }, + "source": [ + "## Estimating the hazard function\n", + "\n", + "Now, to compute the hazard function, we need to know the number of dogs adopted at each value of `t`." + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "KZUyAvh4TTl1", + "outputId": "b26580cb-829d-43f7-fc8d-0bef22044315" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "1 1\n", + "2 1\n", + "4 1\n", + "5 1\n", + "7 0\n", + "dtype: int64" + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
0
11
21
41
51
70
\n", + "

" + ] + }, + "metadata": {}, + "execution_count": 34 + } + ], + "source": [ + "adopted = pd.Series(0, index=ts)\n", + "\n", + "for t in ts:\n", + " k = (shifted['status'] == 1) & (t == shifted['end'])\n", + " adopted[t] = k.sum()\n", + "\n", + "adopted" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "hhQQRVHZTTl1" + }, + "source": [ + "Or here's the same computation with array operations:" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "NIq0xjVkTTl1", + "outputId": "1abfed8f-6051-4c8d-9c23-56638f4f7eb7" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "array([ 5., 1., 4., nan, nan, 2., nan])" + ] + }, + "metadata": {}, + "execution_count": 35 + } + ], + "source": [ + "adopt_times = np.where(shifted['status'], shifted['end'], np.nan)\n", + "adopt_times" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "FVasqpkBTTl2", + "outputId": "f05116e3-66c4-43d2-c9ed-de254e76e32f" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "(5, 7)" + ] + }, + "metadata": {}, + "execution_count": 36 + } + ], + "source": [ + "A, T = np.meshgrid(adopt_times, ts)\n", + "T.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "xZq57sPrTTl2", + "outputId": "fd561f0d-7e5a-4420-8e84-de8728f008bd" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "array([1, 1, 1, 1, 0])" + ] + }, + "metadata": {}, + "execution_count": 37 + } + ], + "source": [ + "adopted = (T == A).sum(axis=1)\n", + "adopted" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "TIzyAaluTTl2" + }, + "source": [ + "For the next step, it will be easier to see what we're doing if we put the results in a table." + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "0Lq98277TTl2", + "outputId": "69815111-ec30-4a71-dd8b-2a14bc41e00a" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " adopted at_risk\n", + "1 1 7\n", + "2 1 6\n", + "4 1 4\n", + "5 1 3\n", + "7 0 1" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
adoptedat_risk
117
216
414
513
701
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 5,\n \"fields\": [\n {\n \"column\": \"adopted\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"at_risk\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 1,\n \"max\": 7,\n \"num_unique_values\": 5,\n \"samples\": [\n 6,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 38 + } + ], + "source": [ + "d = dict(adopted=adopted,\n", + " at_risk=at_risk)\n", + "df = pd.DataFrame(d, index=ts)\n", + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "5BRwXmNzTTl2" + }, + "source": [ + "Finally, the hazard function is the ratio of `adopted` and `at_risk`:" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "wxrpXNhmTTl2", + "outputId": "74ad305c-e067-4962-f304-6f55222cfef7" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " adopted at_risk hazard\n", + "1 1 7 0.142857\n", + "2 1 6 0.166667\n", + "4 1 4 0.250000\n", + "5 1 3 0.333333\n", + "7 0 1 0.000000" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
adoptedat_riskhazard
1170.142857
2160.166667
4140.250000
5130.333333
7010.000000
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 5,\n \"fields\": [\n {\n \"column\": \"adopted\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"at_risk\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 1,\n \"max\": 7,\n \"num_unique_values\": 5,\n \"samples\": [\n 6,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"hazard\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.12485819621073233,\n \"min\": 0.0,\n \"max\": 0.3333333333333333,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.16666666666666666,\n 0.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 40 + } + ], + "source": [ + "df['hazard'] = df['adopted'] / df['at_risk']\n", + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "D3hAv2CHTTl2" + }, + "source": [ + "## Working backwards\n", + "\n", + "Given the hazard function, we can work backwards to compute the survival curve.\n", + "\n", + "The hazard function is the probability of being adopted at each duration, so its complement is the probability of *not* being adopted.\n", + "\n", + "In order to survive past `t`, a dog has to *not* be adopted at all durations up to and including `t`.\n", + "\n", + "So the survival function is the cumulative product of the complement of the hazard function." + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "8kzq6cNuTTl3", + "outputId": "23638e12-e03e-401c-e0fc-249fce481d09" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " adopted at_risk hazard surv\n", + "1 1 7 0.142857 0.857143\n", + "2 1 6 0.166667 0.714286\n", + "4 1 4 0.250000 0.535714\n", + "5 1 3 0.333333 0.357143\n", + "7 0 1 0.000000 0.357143" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
adoptedat_riskhazardsurv
1170.1428570.857143
2160.1666670.714286
4140.2500000.535714
5130.3333330.357143
7010.0000000.357143
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 5,\n \"fields\": [\n {\n \"column\": \"adopted\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"at_risk\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 1,\n \"max\": 7,\n \"num_unique_values\": 5,\n \"samples\": [\n 6,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"hazard\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.12485819621073233,\n \"min\": 0.0,\n \"max\": 0.3333333333333333,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.16666666666666666,\n 0.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"surv\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.2207362448623206,\n \"min\": 0.35714285714285726,\n \"max\": 0.8571428571428572,\n \"num_unique_values\": 4,\n \"samples\": [\n 0.7142857142857144,\n 0.35714285714285726\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 41 + } + ], + "source": [ + "df['surv'] = (1 - df['hazard']).cumprod()\n", + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "I9BsBjMqTTl3" + }, + "source": [ + "The CDF is the complement of the survival function." + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "95J9o2rSTTl4", + "outputId": "cc019da4-c3e2-4303-c443-8a4a1b856bce" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " adopted at_risk hazard surv cdf\n", + "1 1 7 0.142857 0.857143 0.142857\n", + "2 1 6 0.166667 0.714286 0.285714\n", + "4 1 4 0.250000 0.535714 0.464286\n", + "5 1 3 0.333333 0.357143 0.642857\n", + "7 0 1 0.000000 0.357143 0.642857" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
adoptedat_riskhazardsurvcdf
1170.1428570.8571430.142857
2160.1666670.7142860.285714
4140.2500000.5357140.464286
5130.3333330.3571430.642857
7010.0000000.3571430.642857
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 5,\n \"fields\": [\n {\n \"column\": \"adopted\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"at_risk\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 1,\n \"max\": 7,\n \"num_unique_values\": 5,\n \"samples\": [\n 6,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"hazard\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.12485819621073233,\n \"min\": 0.0,\n \"max\": 0.3333333333333333,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.16666666666666666,\n 0.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"surv\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.2207362448623206,\n \"min\": 0.35714285714285726,\n \"max\": 0.8571428571428572,\n \"num_unique_values\": 4,\n \"samples\": [\n 0.7142857142857144,\n 0.35714285714285726\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"cdf\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.22073624486232063,\n \"min\": 0.1428571428571428,\n \"max\": 0.6428571428571428,\n \"num_unique_values\": 4,\n \"samples\": [\n 0.2857142857142856,\n 0.6428571428571428\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 42 + } + ], + "source": [ + "df['cdf'] = 1 - df['surv']\n", + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "hHsDwfegTTl4" + }, + "source": [ + "And the PMF is the difference between adjacent elements of the CDF." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "tags": [], + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "FuCLNtTKTTl4", + "outputId": "ea68d925-375f-4389-d342-fcaec922b389" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " adopted at_risk hazard surv cdf pmf\n", + "1 1 7 0.142857 0.857143 0.142857 0.142857\n", + "2 1 6 0.166667 0.714286 0.285714 0.142857\n", + "4 1 4 0.250000 0.535714 0.464286 0.178571\n", + "5 1 3 0.333333 0.357143 0.642857 0.178571\n", + "7 0 1 0.000000 0.357143 0.642857 0.000000" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
adoptedat_riskhazardsurvcdfpmf
1170.1428570.8571430.1428570.142857
2160.1666670.7142860.2857140.142857
4140.2500000.5357140.4642860.178571
5130.3333330.3571430.6428570.178571
7010.0000000.3571430.6428570.000000
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 5,\n \"fields\": [\n {\n \"column\": \"adopted\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0,\n \"min\": 0,\n \"max\": 1,\n \"num_unique_values\": 2,\n \"samples\": [\n 0,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"at_risk\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2,\n \"min\": 1,\n \"max\": 7,\n \"num_unique_values\": 5,\n \"samples\": [\n 6,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"hazard\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.12485819621073233,\n \"min\": 0.0,\n \"max\": 0.3333333333333333,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.16666666666666666,\n 0.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"surv\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.2207362448623206,\n \"min\": 0.35714285714285726,\n \"max\": 0.8571428571428572,\n \"num_unique_values\": 4,\n \"samples\": [\n 0.7142857142857144,\n 0.35714285714285726\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"cdf\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.22073624486232063,\n \"min\": 0.1428571428571428,\n \"max\": 0.6428571428571428,\n \"num_unique_values\": 4,\n \"samples\": [\n 0.2857142857142856,\n 0.6428571428571428\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"pmf\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.07405871911902757,\n \"min\": 0.0,\n \"max\": 0.1785714285714286,\n \"num_unique_values\": 3,\n \"samples\": [\n 0.1428571428571428,\n 0.1785714285714286\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 43 + } + ], + "source": [ + "df['pmf'] = np.diff(df['cdf'], prepend=0)\n", + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "U-tre_xiTTl4" + }, + "source": [ + "## lifelines\n", + "\n", + "Kaplan-Meier estimation is available in a library called `lifelines`.\n", + "First I'll import it and create a `KaplanMeierFitter`." + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "r4gAQpn2TTl4", + "outputId": "d05f84c7-abfa-47e7-b3f5-f394a2d43756" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Collecting lifelines\n", + " Downloading lifelines-0.30.3-py3-none-any.whl.metadata (3.5 kB)\n", + "Requirement already satisfied: numpy>=1.14.0 in /usr/local/lib/python3.12/dist-packages (from lifelines) (2.0.2)\n", + "Requirement already satisfied: scipy>=1.7.0 in /usr/local/lib/python3.12/dist-packages (from lifelines) (1.16.3)\n", + "Requirement already satisfied: pandas<3.0,>=2.1 in /usr/local/lib/python3.12/dist-packages (from lifelines) (2.2.2)\n", + "Requirement already satisfied: matplotlib>=3.0 in /usr/local/lib/python3.12/dist-packages (from lifelines) (3.10.0)\n", + "Requirement already satisfied: autograd>=1.5 in /usr/local/lib/python3.12/dist-packages (from lifelines) (1.8.0)\n", + "Collecting autograd-gamma>=0.3 (from lifelines)\n", + " Downloading autograd-gamma-0.5.0.tar.gz (4.0 kB)\n", + " Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n", + "Collecting formulaic>=0.2.2 (from lifelines)\n", + " Downloading formulaic-1.2.2-py3-none-any.whl.metadata (7.0 kB)\n", + "Collecting interface-meta>=1.2.0 (from formulaic>=0.2.2->lifelines)\n", + " Downloading interface_meta-2.0.1-py3-none-any.whl.metadata (6.4 kB)\n", + "Requirement already satisfied: narwhals>=1.17 in /usr/local/lib/python3.12/dist-packages (from formulaic>=0.2.2->lifelines) (2.22.1)\n", + "Requirement already satisfied: typing-extensions>=4.2.0 in /usr/local/lib/python3.12/dist-packages (from formulaic>=0.2.2->lifelines) (4.15.0)\n", + "Requirement already satisfied: wrapt>=1.0 in /usr/local/lib/python3.12/dist-packages (from formulaic>=0.2.2->lifelines) (2.2.1)\n", + "Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (1.3.3)\n", + "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (0.12.1)\n", + "Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (4.63.0)\n", + "Requirement already satisfied: kiwisolver>=1.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (1.5.0)\n", + "Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (26.2)\n", + "Requirement already satisfied: pillow>=8 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (11.3.0)\n", + "Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (3.3.2)\n", + "Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.0->lifelines) (2.9.0.post0)\n", + "Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.12/dist-packages (from pandas<3.0,>=2.1->lifelines) (2025.2)\n", + "Requirement already satisfied: tzdata>=2022.7 in /usr/local/lib/python3.12/dist-packages (from pandas<3.0,>=2.1->lifelines) (2026.2)\n", + "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.12/dist-packages (from python-dateutil>=2.7->matplotlib>=3.0->lifelines) (1.17.0)\n", + "Downloading lifelines-0.30.3-py3-none-any.whl (409 kB)\n", + "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m409.1/409.1 kB\u001b[0m \u001b[31m13.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", + "\u001b[?25hDownloading formulaic-1.2.2-py3-none-any.whl (118 kB)\n", + "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m118.9/118.9 kB\u001b[0m \u001b[31m6.3 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", + "\u001b[?25hDownloading interface_meta-2.0.1-py3-none-any.whl (15 kB)\n", + "Building wheels for collected packages: autograd-gamma\n", + " Building wheel for autograd-gamma (setup.py) ... \u001b[?25l\u001b[?25hdone\n", + " Created wheel for autograd-gamma: filename=autograd_gamma-0.5.0-py3-none-any.whl size=4030 sha256=3bdfcf25da65de9ea3c0a5bb4050a130af0e3d05dd642afbcbf7c6045ef9b765\n", + " Stored in directory: /root/.cache/pip/wheels/50/37/21/0a719b9d89c635e89ff24bd93b862882ad675279552013b2fb\n", + "Successfully built autograd-gamma\n", + "Installing collected packages: interface-meta, autograd-gamma, formulaic, lifelines\n", + "Successfully installed autograd-gamma-0.5.0 formulaic-1.2.2 interface-meta-2.0.1 lifelines-0.30.3\n" + ] + } + ], + "source": [ + "# If we're running in Colab, install lifelines\n", + "\n", + "import sys\n", + "IN_COLAB = 'google.colab' in sys.modules\n", + "\n", + "if IN_COLAB:\n", + " !pip install lifelines" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": { + "id": "_reKAgIJTTl4" + }, + "outputs": [], + "source": [ + "from lifelines import KaplanMeierFitter\n", + "kmf = KaplanMeierFitter()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Jbqs07QYTTl5" + }, + "source": [ + "Now we need two sequences, the durations, including complete and ongoing cases." + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 304 + }, + "id": "aB8mrOgiTTl5", + "outputId": "e704c78b-81b0-44e4-cb66-6c073af3a201" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "0 5\n", + "1 1\n", + "2 4\n", + "3 7\n", + "4 5\n", + "5 2\n", + "6 2\n", + "Name: end, dtype: int64" + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
end
05
11
24
37
45
52
62
\n", + "

" + ] + }, + "metadata": {}, + "execution_count": 46 + } + ], + "source": [ + "T = shifted['end']\n", + "T" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "LJVpjzlOTTl5" + }, + "source": [ + "And an event flag that indicates whether a case is complete." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": { + "id": "sXDZIyiwTTl5" + }, + "outputs": [], + "source": [ + "E = shifted['status']" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pM0MTgO9TTl5" + }, + "source": [ + "The `fit` method does the Kaplan-Meier estimation." + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "FhEdJUFFTTl5", + "outputId": "3202c08c-027b-481c-d87b-ddc9303383c0" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 48 + } + ], + "source": [ + "kmf.fit(T, E)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "jpqJ3EPJTTl6" + }, + "source": [ + "Now the `Fitter` object contains the estimated survival function." + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 269 + }, + "id": "_CWAN-YKTTl6", + "outputId": "5c21ee25-28cb-4e69-d189-b4f9258bdc8a" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " KM_estimate\n", + "timeline \n", + "0.0 1.000000\n", + "1.0 0.857143\n", + "2.0 0.714286\n", + "4.0 0.535714\n", + "5.0 0.357143\n", + "7.0 0.357143" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
KM_estimate
timeline
0.01.000000
1.00.857143
2.00.714286
4.00.535714
5.00.357143
7.00.357143
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "summary": "{\n \"name\": \"kmf\",\n \"rows\": 6,\n \"fields\": [\n {\n \"column\": \"timeline\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 2.6394443859772205,\n \"min\": 0.0,\n \"max\": 7.0,\n \"num_unique_values\": 6,\n \"samples\": [\n 0.0,\n 1.0,\n 7.0\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"KM_estimate\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.2657456458708585,\n \"min\": 0.3571428571428571,\n \"max\": 1.0,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.8571428571428572,\n 0.3571428571428571,\n 0.7142857142857143\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 49 + } + ], + "source": [ + "kmf.survival_function_" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "g3hQpEZOTTl6" + }, + "source": [ + "`timelines` includes an element at `t=0`, but other than that it is identical to what we computed (except for floating-point error)." + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "Zd3UwtacTTl6", + "outputId": "1897f2d1-247e-4cc4-ccef-b9104f7149e7" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "1.6653345369377348e-16" + ] + }, + "metadata": {}, + "execution_count": 50 + } + ], + "source": [ + "max(abs(kmf.survival_function_['KM_estimate'] - df['surv']).dropna())" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "n-hL5zkiTTl6" + }, + "source": [ + "`lifelines` also computes a confidence interval for the survival function." + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 238 + }, + "id": "xaUmvSyQTTl6", + "outputId": "ec3308ab-973f-4566-85ee-6517ba45aa2c" + }, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " KM_estimate_lower_0.95 KM_estimate_upper_0.95\n", + "0.0 1.000000 1.000000\n", + "1.0 0.334054 0.978561\n", + "2.0 0.258154 0.919797\n", + "4.0 0.131988 0.824997\n", + "5.0 0.051977 0.698713\n", + "7.0 0.051977 0.698713" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
KM_estimate_lower_0.95KM_estimate_upper_0.95
0.01.0000001.000000
1.00.3340540.978561
2.00.2581540.919797
4.00.1319880.824997
5.00.0519770.698713
7.00.0519770.698713
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "ci", + "summary": "{\n \"name\": \"ci\",\n \"rows\": 6,\n \"fields\": [\n {\n \"column\": \"KM_estimate_lower_0.95\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.35889775427078546,\n \"min\": 0.051976519457499495,\n \"max\": 1.0,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.3340538792922218,\n 0.051976519457499495,\n 0.2581536654587949\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"KM_estimate_upper_0.95\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 0.13433413016538617,\n \"min\": 0.6987130319850487,\n \"max\": 1.0,\n \"num_unique_values\": 5,\n \"samples\": [\n 0.9785610585261756,\n 0.6987130319850487,\n 0.9197974560448583\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 51 + } + ], + "source": [ + "ci = kmf.confidence_interval_survival_function_\n", + "ci" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "c1EeXiHFTTl6" + }, + "source": [ + "With such a small dataset, the CI is pretty wide." + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 449 + }, + "id": "xb4ohYm3TTl9", + "outputId": "45b0b386-38b6-4a87-960a-77f521af653b" + }, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYSNJREFUeJzt3Xd4lGW+PvD7nd5beqNIkV6kCYgiIliWI1aOoiKWPavoqiwWLKDoApbdVVeOBde2yg9cj2UVAREBWcUFQZpAIBAIhPSpmWRmkpn39wfObCIQM5CZd8r9ua65JG+mfCcmM/c87/M8X0EURRFEREREKUImdQFEREREHYnhhoiIiFIKww0RERGlFIYbIiIiSikMN0RERJRSGG6IiIgopTDcEBERUUpRSF1AvIVCIRw7dgxGoxGCIEhdDhEREbWDKIrweDzIz8+HTNb22EzahZtjx46hqKhI6jKIiIjoNBw5cgSFhYVtXiftwo3RaARw/IdjMpkkroaIiIjaw+12o6ioKPI+3pa0CzfhU1Emk4nhhoiIKMm0Z0oJJxQTERFRSmG4ISIiopTCcENEREQpJe3m3BARkXSCwSCampqkLoMSlEql+tVl3u3BcENERDEniiIqKyvhdDqlLoUSmEwmQ9euXaFSqc7ofhhuiIgo5sLBJjs7Gzqdjpuo0gnCm+xWVFSgU6dOZ/Q7wnBDREQxFQwGI8EmIyND6nIogWVlZeHYsWNobm6GUqk87fvhhGIiIoqp8BwbnU4ncSWU6MKno4LB4BndD8MNERHFBU9F0a/pqN8RhhsiIiJKKZKGm2+++QaTJk1Cfn4+BEHAJ5988qu3WbduHc455xyo1Wp0794db7/9dszrJCIiouQhabjxer0YOHAgFi1a1K7rl5aW4vLLL8eFF16Ibdu24b777sPtt9+OVatWxbhSIiKi5HPLLbdg8uTJUpcRd5KGm0svvRRPP/00rrzyynZd/9VXX0XXrl3xpz/9Cb1798bdd9+Na665Bn/5y19iXGn7fFdSC6+/WeoyiIiog5wsHHz44YfQaDT405/+hFtuuQWCIOB3v/vdCbedMWMGBEHALbfcEvM6Dx06BEEQsG3btlbHX3zxxbic4Ui0EJVUc242btyI8ePHtzo2ceJEbNy48ZS38fv9cLvdrS6xsP2IE7e8vRlXvPwv7KtwQBTFmDwOERFJ54033sDUqVPxyiuv4A9/+AMAoKioCEuXLkVjY2Pkej6fD0uWLEGnTp2kKhUAYDabYbFYJK1BCkkVbiorK5GTk9PqWE5ODtxud6tfqpYWLFgAs9kcuRQVFcWktuaQCJNGgZIaLyYv2ojFX3yPPXv24ODBgzh27Biqq6vhcDhQX18Pn893xsvciIiSmSiKaAg0S3I53Q+fzz77LO655x4sXboU06dPjxw/55xzUFRUhI8++ihy7KOPPkKnTp0wePDgdt9/KBTCggUL0LVrV2i1WgwcOBAffvhh5PsOhwNTp05FVlYWtFotevTogbfeegsA0LVrVwDA4MGDIQgCxo4dC+DEEZWxY8finnvuwX333Qer1YqcnBwsXrwYXq8X06dPh9FoRPfu3bFixYrIbYLBIG677bZIXWeffTZefPHFyPefeOIJvPPOO/j0008hCAIEQcC6desAAEeOHMF1110Hi8UCm82GK664AocOHWr3z+R0pfwmfrNnz8bMmTMjX7vd7pgEnCGdrfjnjFG47W/fYU9tAPM32HF1jQ9Xn60BWvwhyWQyKBQKyOVyKJVKqNVqaDQaKJVKKJVKKBSKVv/m0kkiSkWNTUH0mSPNfMnd8yZCp4ru7e+hhx7C//7v/+Lzzz/HRRdddML3b731Vrz11luYOnUqAODNN9/E9OnTI2/y7bFgwQK89957ePXVV9GjRw988803uPHGG5GVlYULLrgAjz/+OHbv3o0VK1YgMzMTJSUlkQ/2mzZtwvDhw/HVV1+hb9++bbYveOedd/Dggw9i06ZNWLZsGe688058/PHHuPLKK/HII4/gL3/5C2666SaUlZVBp9MhFAqhsLAQ//jHP5CRkYHvvvsOv/3tb5GXl4frrrsOs2bNwp49e+B2uyNhy2azoampCRMnTsTIkSOxYcMGKBQKPP3007jkkkuwY8eOM26x0JakCje5ubmoqqpqdayqqgomkwlarfakt1Gr1VCr1fEoD/lWPV64oiv+svYwVh704f/2NuCQB3j4ghyY1HIAxxNw+OL3+9HQ0BAZxRFFEYIgQC6XQy6XQ6FQQKVSQaPRQK1Wtwo/4f/K5fK4PDcionS1YsUKfPrpp1izZg3GjRt30uvceOONmD17Ng4fPgwA+Pbbb7F06dJ2hxu/34/58+fjq6++wsiRIwEAZ511Fv71r3/htddewwUXXICysjIMHjwYQ4cOBQB06dIlcvusrCwAQEZGBnJzc9t8rIEDB+Kxxx4DcHwAYOHChcjMzMQdd9wBAJgzZw5eeeUV7NixA+eeey6USiWefPLJyO27du2KjRs34oMPPsB1110Hg8EArVYLv9/f6rHfe+89hEIhvPHGG5EP6m+99RYsFgvWrVuHCRMmtOtnczqSKtyMHDkSX3zxRatjq1evjvwiJAKr2YSb+mowsNCMF76txpbyBtz9zyOYMy4P3TPUkeByKqFQCKFQCM3NzQgGg6ivr4fL5UIoFIpc59dGf34Zgjj6Q0SJRquUY/e8iZI9djQGDBiA2tpazJ07F8OHD4fBYDjhOllZWbj88svx9ttvQxRFXH755cjMzGz3Y5SUlKChoQEXX3xxq+OBQCByauvOO+/E1Vdfja1bt2LChAmYPHkyRo0aFdVzCT+fMLlcjoyMDPTv3z9yLDz9o7q6OnJs0aJFePPNN1FWVobGxkYEAgEMGjSozcfZvn07SkpKYDQaWx33+Xw4cOBA1HVHQ9JwU19fj5KSksjXpaWl2LZtG2w2Gzp16oTZs2ejvLwc7777LgDgd7/7HV5++WU8+OCDuPXWW/H111/jgw8+wPLly6V6CifQ6/VQKpUY00mDLtZCzPu6AhWeZty//Ch+PzILF/cwtXl7mUwWCS8nI4oiQqHQCaM/zc3/WaUlCAIUCkXkfk41+hP+N0d/iCjeBEGI+tSQVAoKCvDhhx/iwgsvxCWXXIIVK1ac8IYNHD81dffddwNAu7c4CauvrwcALF++HAUFBa2+Fz77cOmll+Lw4cP44osvsHr1alx00UWYMWMGnn/++age65c9mwRBaHUs/IE4/KF66dKlmDVrFv70pz9h5MiRMBqNeO655/Dvf//7V5/TkCFD8P7775/wvfBIU6xI+pv1ww8/4MILL4x8HZ4bM23aNLz99tuoqKhAWVlZ5Ptdu3bF8uXLcf/99+PFF19EYWEh3njjDUycKE36PxmtVguNRgO/34+zbAb8dVIRnv2mCpuONuD5f1WjuNaP/xmeCaX89EZTWp62OpWW4ae9oz/h8KNSqU449cXRHyJKd507d8b69esjAWflypUnBJxLLrkEgUAAgiBE/b7Up08fqNVqlJWV4YILLjjl9bKysjBt2jRMmzYNY8aMwQMPPIDnn3++w3oyncy3336LUaNG4a677ooc++XIi0qlOuGxzznnHCxbtgzZ2dkwmdr+YN/RJA03Y8eObXPW+snW5o8dOxY//vhjDKs6MzKZDBaLBeXl5QAAo1qOJ8fnYck2B97bZsdne10oqfPjsQtzkamPzY8/PPpzqo6qvxz98fl88Hq9CIVCkXk/AE4Y/dFqtVCr1ScNPxz9IaJUV1RUhHXr1uHCCy/ExIkTsXLlylbfl8vl2LNnT+Tf0TAajZg1axbuv/9+hEIhnHfeeXC5XPj2229hMpkwbdo0zJkzB0OGDEHfvn3h9/vx+eefo3fv3gCA7OxsaLVarFy5EoWFhdBoNDCbzR3yvHv06IF3330Xq1atQteuXfH3v/8dmzdvjqzQAo7P/1m1ahWKi4uRkZEBs9mMqVOn4rnnnsMVV1yBefPmobCwEIcPH8ZHH32EBx98EIWFhR1S38kkx5hgkjEYDBBFMRIUZIKAGwfb0CNTjWe/qcKeGh9m/PMIHrswF/1zTz4ROpZiPfqj1+thMBhiOhOeiEgKhYWFrQJOXl5eq++fyQjFU089haysLCxYsAAHDx6ExWLBOeecg0ceeQTA8dGR2bNn49ChQ9BqtRgzZgyWLl0K4PiH0Zdeegnz5s3DnDlzMGbMmKhWarXlf/7nf/Djjz9iypQpEAQB119/Pe66665Wy8XvuOMOrFu3DkOHDkV9fT3Wrl2LsWPH4ptvvsFDDz2Eq666Ch6PBwUFBbjoootiPpIjiGm225zb7YbZbIbL5YrZD9fn82HXrl1QqVQnrNQ65m7CvK8rUOoIQCYAdwzLxJV9zEl32ueXoz/hCdDhYcnwpwaLxQKDwRC3FWtElHh8Ph9KS0vRtWtXaDQaqcuhBNbW70o0798cuYkBtVoNnU4Hr9d7wpt6vkmJF35TiBe/rcbXB+vx2qZaFNf4cP/obGiUybOnYlujP6FQCH6/HzU1NaiqqoJarYbJZILFYoHRaOSLGxERxRTDTQwIggCr1Qqn03nS72sUMjx4fg56ZWnw2qZarCutxyFHAHPG5aLAnPyncmQyGbRaLbRaLURRhN/vR11dHWpqaqBWq2EwGGC1WiNBJ9lGrYiIolVWVoY+ffqc8vu7d++WvFVDKmG4iRGdTgeZTIZgMHjS0Q1BEHBFHwu6Zajxx7WVOOQM4O7PjuKh83Nwbie9BBXHhiAI0Gg00Gg0EEURgUAATqcTtbW1UKlUMBqNkaCj1WoZdIgoJeXn55/Q1PKX36eOw3ATI3q9Hmq1GoFA4JS7JwNAvxwtXv6vIvxxbSV+qvZh7poK3DDQihsH2SCXpdYbvSAIkR2jw0HH7Xajrq4OKpUKer0eNpsNBoMBer2eQYeIUoZCoUD37t2lLiNtMNzEiEKhgMlkQk1NTZvhBgAydAo8c0kBFm+uxad7XFiy3YF9tX481KJtQ6ppGXSA47twer1eOBwOKJVK6PX6yIiOXq+HTJY885GI6ORarrYkOpmOWuPEcBNDRqMRlZWV7bquUi7grnOzcHaWBi9+W40fyhtwz89tG7plpP5KI5VKFVk63tTUhIaGBrhcLsjlcuh0OthstkjQ4Z46RMlFpVJBJpPh2LFjyMrKgkql4sgsnUAURdTU1JywY/Lp4FLwGPJ6vfjpp5+g0+mi+h91wO7HUz+3bVDJBfx+VBYu7h7f3R0TRXNzc6SPiUwmg06ng9VqhclkgsFgYNAhShKBQAAVFRVoaGiQuhRKYIIgoLCw8KT9u6J5/2a4iaFQKIRdu3ahqanppP+j2uLxByNtGwBgUi/zGbVtSAXNzc3w+Xzw+/2RFVktg86p+nERUWIQRTGyJxbRySiVylN+aGW4aUM8ww1wfPlfeXk5bDZb1LcNiSLe32bHe9scAIDeWRo8Pi4XGTq+iYfbRvj9fgDHNw20Wq0wm80wGAxnPKRJRESJheGmDfEON3a7HXv37oXNZjvtc8zfH/Hi2W+q4A2EYNXK8ehYado2JKpQKITGxkb4/X6IogitVguz2Qyz2Qyj0cg2EEREKYDhpg3xDjdttWKIRrk7gKe+rkz6tg2xFgqF4PP54PP5IIpipA1EOOiwDQQRUXJiuGlDvMONKIrYs2cPvF7vGT+erymEF76rxtqD9QCAsWcZcP+o5GrbEE/hNhA+nw+hUIhtIIiIkhh7SyWQX2vFEA2NUoaHzs/B2ZkaLN5ci3UH63HIHsCci3JRYOKpl19qqw1Ey92RDQYDd0cmIkohDDdx8GutGKIhCAKu7GtB9ww1/rjueNuGez47igfPz8G5RanTtqGjnaoNRF1dHZRKZaugo9PpGHSIiJIYT0vFQXNzM3bu3AlBEH51t+Jo1DU04+m1ldhd7QOAlG3bEEuiKKKpqQmNjY1obm6GUqmETqdDRkYG20AQESUQzrlpgxThBgAOHDiAmpoaWK3WDr3fpqCI1zfX4p97XACAYQU6PJjCbRtiLRAIwOfzIRAIRIJOy92R2QaCiEgaDDdtkCrcVFdXo6SkBBkZGTG5/69K3HjxuxoEgiJyDYq0adsQS+ERnaamJsjlcmi1WthsNphMJraBICKKM04oTkB6vR4KhQJNTU0x2WBufHcTulrVmPd1BSrrm3Hf8qO4d1QWxqdp24aOoFQqI/+vwm0gysrK2AaCiCjBceQmTs6kFUM03P4gnllfhR/K2bYhVoLBYGTTQEEQIkHHaDTCaDSyDQQRUQzwtFQbpAo3wPFWDEePHo3ZqamwYEjEe9vsWLL9eNuGPtkaPHYh2zbEQrgNhM/ni6zICu+jo9PpoNFoOE+HiKgDMNy0Qcpw0xGtGKLxfZkXz2443rbBppXjEbZtiKnw7sh+vx+hUAhKpRJqtRoWiwV6vT4Sdrj6iogoegw3bZAy3Ph8Pvz000+RN714KHcHMG/N8f1w5D+3bZjMtg0xF+5+7Pf7EQgEIIoiVCpVpO9VOOywHQQRUftwQnGCUqvV0Gq18Hq9cXtTKzCp8OJvCvGX76qx7mA9Xt1Ui+JaH+5j24aYEgSh1YTklvvpuFzHl+2r1WrodLpWYYfdzImIzhzDTRyFWzE4HI64Pq5GKcPD5+egV6YGr2+uxdqD9TjkCODxcWzbEC+CIEClUkU6lIdCIQQCAdTX18PhcEAQBKjVahgMBphMJuh0Ouh0Ok5OJiI6DXzljDOdTge5XN4hrRii8cu2DaUOtm2Qkkwmi7SDAP7T5NPhcKCmpgZyuTzS6DM8OVmr1XLJORFRO3DOTZyFWzEAx4OOFOoamvHU15XYU3O8bcONg6yYOsgGGefhJIxgMAi/3x+ZnCyXy6HRaFqFHa7EIqJ0wgnFbZA63ACxa8UQjaagiNc21eKzvcfnfwwv1OHB83NgZNuGhNRycjJXYhFROmK4aUMihJtYt2KIxuoSN176uW1DnlGBx8floZuNK3gSXVNTU6uVWOE+WC0nJ6tUKoYdIkoZXC2V4GLdiiEaF7do21Dhacb9nx/FvaOzcVE3o6R1UdvasxIrvOzcYDBwJRYRpRWO3EggXq0YovHLtg1X9DbjjmFs25CMQqFQZGSnqamJK7GIKCXwtFQbEiHcAPFrxRCNX7Zt6JutwaNs25D0wiux/H4/mpubuRKLiJISw00bEiXcxLsVQzR+2bbh0Qtz0S+HbRtSRVsrscKnsLRaLVdiEVFCYbhpQ6KEG7/fj127dsW1FUM0yl0BzPv6P20b/md4Jv6rN9s2pKLwSiy/3x+ZnMyVWESUaBhu2pAo4UYURezduxcejwdms1myOtriawrhL99WY11pPQBg3FkG3Ds6GxoFP9GnspOtxNJqta3CDldiEVG8cbVUEhAEARaLBXa7XepSTkmjlOHhC3JwdpYGizfX4uuD9Sh1BDBnXB7yTVx5k6pOthLL5/OhrKwMwIkrsbRabaStBBFRImC4kZBUrRiiIQgCrgq3bVgbbttwBA+dn4PhbNuQ8k7WE6upqQlerxdOpzOyEkuv18NsNnMlFhElBJ6WklAitGKIRq23GU+vrcCeGj8EADcOsuGGQVa2bUhj4ZVYgUAATU1NXIlFRDHDOTdtSKRwAyRGK4ZosG0DtYUrsYgoVhhu2pBo4SaRWjFE45dtG+aMy8NZbNtAv9ByJVYoFIJKpYJarW61czJXYhFRe3BCcRJJpFYM0fhl24b7Pj+K+0ZnYxzbNlALCoUCCoUCev3x+VnhlVjHjh3jSiwiihmO3EgsEVsxRMPtC2LhN1XY0qJtw2+HZ0Ih45sTtS28Eis8ZwdovRKrZdghIuJpqTYkWrgBErMVQzROaNuQo8GjY9m2gaIjiiICgcAJPbG4EouIAJ6WSjrhERtRFJNyOF4uEzDtnAz0zNTg2W+q8FOVD3f/8wgeuzAXfdm2gdopHGbCO3aHQiEEAgG4XC7U1tZGVmIZjcZIA1CuxCKik+HITQJI9FYM0Tj6c9uGw2zbQB2MK7GI0htPS7UhEcNNMrRiiEbjz20b1v/ctuGibkb8flQW2zZQh+JKLKL0wtNSSSbcisHhcEhdSofQKmWYfUEOev3ctmHNAQ9KHX7MGZeHPGPyrAijxMaVWER0Kgw3CUKv10MmkyV0K4ZohNs2dLOpMH9dFQ7aA7j7n2zbQLETTU8srsQiSm08LZUgkq0VQzRqvM34I9s2kIS4Eoso+XHOTRsSNdwAydeKIRqBoIjXNtXg871uAMCIIh0eHJMDA9s2kATCK7H8fj+am5shk8mgVqthMBgiYYcrsYgSC8NNGxI53NTU1GD//v1Ju99Ne3y5342XNtagKSgiz6jEnHG5bNtAkuNKLKLExwnFSSo8LJ5srRiiMaGHCV2tKjy1thIVnia2baCEIJfLI6emgP+sxKqsrEQoFIJSqYRGo+FKLKIkwZGbBJLsrRii4fYFsXB9JbYcawQATO5jxh3D2LaBElPLNhEtV2K1DDtciUUUWzwt1YZEDjdA8rdiiEYwJOLvP9rx/3YcXwLf7+e2DTa2baAE1p6eWAaDgZOTiToYT0slsWRvxRANuUzALUMy0DNTjec2VGNXlQ8z2LaBEpwgCFCpVJFl5OGVWF6vF06nE4IgwGw2o7CwMCE/QBGlA86OSzB6vR5qtTryiTAdjOpswF8nFaKTRQV7YxAPrCjHP/c4kWaDipSkwsvKTSYTMjIyYDab4Xa7UVxcjGPHjiEYDEpdIlHaYbhJMCqVCjqdDj6fT+pS4qrQrMJLvynE+V0MCIrAou9r8dyGaviaQ1KXRhQVuVwOq9UKhUKB0tJSlJSUwOv1Sl0WUVphuEkw4VYMzc3NUpcSd1qlDI+MzcEdwzIgE4A1Bzy4f/lRVHiapC6NKGo6nQ5WqxV2ux3FxcWorq5GKMSwThQPDDcJqGUrhnQjCAKu6WfFwon5MGvkkbYNm4/yky8ln/AojiiKKCkpwcGDB9NuVJZICpKHm0WLFqFLly7QaDQYMWIENm3a1Ob1X3jhBZx99tnQarUoKirC/fffn3IvFjqdDmq1Gn6/X+pSJDMwT4dF/1WEXllq1AdCeHx1Bd7fZkeI83AoyQiCAIPBAJPJhKqqKhQXF6Ouro5zyohiSNJws2zZMsycORNz587F1q1bMXDgQEycOBHV1dUnvf6SJUvw8MMPY+7cudizZw/+9re/YdmyZXjkkUfiXHlsKRQKmEymtA43AJClV+C5Swtx2dkmiADe/dGOJ9dUoN6ffiNalPyUSiUyMjIQCASwf/9+HDp0KK0WDhDFk6Th5s9//jPuuOMOTJ8+HX369MGrr74KnU6HN99886TX/+677zB69GjccMMN6NKlCyZMmIDrr7/+V0d7kpHJZOL5eQAquYB7R2Vj5nnZUMoFfH+kAfd8dhSl9vQOfpScBEGAyWSCXq/HsWPHsG/fPjidTqnLIko5koWbQCCALVu2YPz48f8pRibD+PHjsXHjxpPeZtSoUdiyZUskzBw8eBBffPEFLrvsslM+jt/vh9vtbnVJBi1bMRAwsYcJf7msADkGBY55mnDv8qNYd9AjdVlEp0WlUsFms8Hr9WLfvn04cuRIWi4iIIoVycJNbW0tgsEgcnJyWh3PyclBZWXlSW9zww03YN68eTjvvPOgVCrRrVs3jB07ts3TUgsWLIDZbI5cioqKOvR5xIpWq4VGo0n7U1Mt9cjU4OVJRTgnXwt/s4gF66vw6r9r0Bzi3AVKPjKZDBaLBWq1GmVlZdi/fz88HgZ2oo4g+YTiaKxbtw7z58/H//7v/2Lr1q346KOPsHz5cjz11FOnvM3s2bPhcrkilyNHjsSx4tMXfuFjuGnNpJHj6Yvz8d8DrACAj3e78PDKctgb+KmXkpNGo4HVaoXT6cS+fftQUVHBU9JEZ0iy9guZmZmQy+Woqqpqdbyqqgq5ubknvc3jjz+Om266CbfffjsAoH///vB6vfjtb3+LRx99FDLZiVlNrVZDrVZ3/BOIg3RqxRANuUzA9CEZODtTjec2VGFnlQ93f3a8bUOfbLZtoOQjl8sjp6lKS0vh8XhQWFgY6VJORNGRbORGpVJhyJAhWLNmTeRYKBTCmjVrMHLkyJPepqGh4YQAI5fLASAll1WmYyuGaIzqbMBLk4rQyaJCXcPxtg2f7XGl5O8CpQe9Xg+z2Yza2loUFxejpqaGv89Ep0HS01IzZ87E4sWL8c4772DPnj2488474fV6MX36dADAzTffjNmzZ0euP2nSJLzyyitYunQpSktLsXr1ajz++OOYNGlSJOSkEpVKBb1en3L7+HSkop/bNozpYkBzCHj5+xo8v6EafrZtoCSlUChgs9kQDAZRUlKC0tJSnp4mipKkXcGnTJmCmpoazJkzB5WVlRg0aBBWrlwZmWRcVlbWaqTmsccegyAIeOyxx1BeXo6srCxMmjQJf/zjH6V6CjEV7i5st9ulLiWhaZUyPDo2B//3kxp/+6EOXx3w4JAjgMfH5SLXqJS6PKKoCYIAo9GIpqYmVFRUoL6+HoWFhbBarTxFTdQOgphmY55utxtmsxkulwsmk0nqcn6V2+3G7t27YTKZUnJ0qqNtq2jA/HVVcPmCMKhkmD02B0ML9FKXRXTaRFGE2+1GKBRCXl4e8vPzoVQytFP6ieb9O6lWS6UjtmKIzqA8HV6eVIizM4+3bXjsywosYdsGSmLhEVydToejR49i3759SbNfF5FUGG4SHFsxRC/boMTzl/2nbcM7P9rx5JpKeANs20DJS61Ww2azwePxoLi4GOXl5WnZXJeoPRhukgBbMUTvxLYNXtzz2VEccjAkUvKSyWSwWq1QKpU4dOgQ9u/fD6/XK3VZRAmH4SYJsBXD6Qu3bcjWK1DubsLvP2fbBkp+Wq0WVqsVdrsde/fuRVVVFT8AEbXAcJME2IrhzPTI1ODl/yrC4BZtG17bxLYNlNzkcjkyMjIAAAcOHMDBgwfR2NgocVVEiYHhJgmwFcOZM2vk+OPF+Zjyc9uGj35yYfaqcjga2baBkpvBYIDZbEZ1dTWKi4tRW1vLjf8o7THcJAmDwQBBEPiidQbkMgG3DsnAnHG50CkF7Kj0YcY/j2BPNTdJpOQW3vivqakJ+/fvx6FDh7izOaU1hpskodfroVKp+ILVAUZ3NuDF3xShk1mJuoYgZq04is/2sm0DJTdBEGAymWAwGHDs2DEUFxfD6XRKXRaRJBhukoRarWYrhg7UyaLCi5OKMKaL/njbho01+NO/2LaBkp9KpYLNZkNDQwP27duHsrIyNDfz9CulF4abJGI2m/ki1YF0ShkeHZuL24dmQCYAq0s8mLm8HJUerkqj5Baep6dWq3HkyBHs27cPHg9XCVL6YLhJInq9HjKZjBt3dSBBEHBtfysWTMiHWS1Did2Puz87gi3lDVKXRnTGNBoNrFYr3G43iouLcezYMb5+UFpguEkibMUQO4PydXj5v4pwdqYaHn8Ij355DEu2s20DJT+5XA6r1Qq5XI5Dhw6hpKQEDQ0M75TaGG6SCFsxxFakbUPPn9s2bLVj3tds20CpQafTwWw2o66uDnv37kV1dTUn0VPKYrhJMmzFEFsquYB7R2fj/tHH2zZsLGPbBkod4SXjoihGNv7jhyVKRQw3SUan00GpVLIVQ4xd0tOEP7do23Dv50exvpQTMin5CYIAg8EAo9GIyspKFBcXw263cxSHUgrDTZLRarWcdxMnPVu0bfA1i5i/rgqvbapFkG0bKAUolUpkZGTA5/NFlozzQxOlCoabJMNWDPEVadvQ3wIA+OgnJx5edQxOtm2gFCAIAsxmM3Q6HY4ePYri4mK4XC6pyyI6Yww3SchgMAAAh5HjRC4TcOvQTDx+YS60CgE7KhvZtoFSilqths1mg9frxb59+3D06FHuqUVJjeEmCen1eqjVarZiiLPzuhjw0qQiFJmVqP25bcPnbNtAKSI8KqxUKnH48GGUlJSgvr5e6rKITgvDTRJiKwbpdLKo8NKkIpzX+Xjbhr9urMGfv2XbBkodWq0WVqsVDocDxcXFqKio4ApNSjoMN0mKrRiko1PK8NiFubjt57YNX+73YOYXbNtAqUMul8Nms0EQBJSWlqKkpASNjY1Sl0XUbgw3SYqtGKQlCAKu62/F/HDbhjq2baDUo9frYTabUVtbi71796K2tpanYSkpMNwkKbZiSAyDf27b0LNF24ZXvq9BXQNH1Sg1hDf+CwaD2L9/P0pLS/m6QwmP4SZJsRVD4sg2KPGnSwtw6c9tGz7Z48K0Dw/j5Y01qK7nqSpKfoIgwGg0wmAwoKKiAsXFxXA4HBzFoYQliGn22+l2u2E2m+FyuWAymaQu54zU1NRg3759yMzMlLoU+tmW8ga8v92On6qOT/ZWyIAJPUyY0t+KXKNS4uqIzlwoFILH44EoisjLy0NeXh6USv5uU+xF8/7NcJPEvF4vdu/eDa1WyxeXBCKKInZUNuL9bQ5srzw+CVMuABd3N2HKACvyTfx/RcnP7/fD4/HAYrGgqKgo6V9PKfEx3LQhlcKNKIrYuXMnmpqaIhv7UWLZWdmIJdvt2HrseMiRCcC4s4y4fqAVhWaVxNURnZlQKASXywW5XI78/Hzk5uZCLpdLXRalKIabNqRSuAGAI0eO4MiRI8jIyJC6FGrDnmoflmy3Y9PR46upZAJwflcDbhhoQ2cLQw4lt8bGRni9XmRmZqKwsBB6vV7qkigFMdy0IdXCjd1ux969eyN7UlBi21frw/vbHPj+iBcAIOD4zsc3DLTiLJta2uKIzkAwGITL5YJKpUJBQQFPlac5lUoFo9HYoffJcNOGVAs3fr8fu3btgkKhgEajkbocaqeSOj+WbLfj28PeyLFRnfS4YaAVPTL5/5GSkyiK8Hq98Pl8kMm4GDddhUIh2Gw29O7du0PvN5r3b0WHPjLFXbgVg9vtZrhJIt0z1JgzLg+ldj/+3w4Hvimtx3dlXnxX5sWIIh1uGGhDryz+/6TkIggCDAYD5wCmOY/HI3nLDkbrFMBWDMmrq02NR8bm4vUrO2HcWQbIBODfRxpw7+dH8eiXx/BTFbe8JyKKFkduUkDLVgxcqZCcOllUeOiCXNw4KID/t8OBNQc8+KG8AT+UN2BQnhY3DrKhf65W6jKJiJLCaYebQCCA6urqE4aeOnXqdMZFUXRatmLQ6XRSl0NnoMCswqwxOZg6yIalOxxYvd+NbRWN2FZRjgG5GkwdaMPAPC0njxMRtSHqcLN//37ceuut+O6771odF0URgiCwkaMEwq0YampqGG5SRJ5RiftHZ+OGgVZ8sNOBVfvc2FHpw47KY+iTrcHUQVYMydcx5BARnUTU4eaWW26BQqHA559/jry8PL64JgiTyYTKykqpy6AOlmNQ4p6R2fjvATb8Y6cDK/a5sbvah0e/rEDPTDWmDrJhRCFDDhFRS1EvBdfr9diyZQt69eoVq5piKtWWgoexFUN6qGtoxoe7nFi+1wV/8PifbnebGlMHWXFuJz1kDDlEJDGPxwOtVou+fft26P1G8/4d9WqpPn36oLa29rSLo9gIz7vx+XxSl0IxlKFT4H+GZ+Ldazvj2n4WaBQCSux+PPl1Je769Ai+Ka1HKL22riIiOkHU4eaZZ57Bgw8+iHXr1qGurg5ut7vVhaQhCAKsVisCgYDUpVAcWLQK3D4sE+9e2wXXD7BCpxRQ6gjgj+sq8T+fHMHagx4EQww5RJSeoj4tFd518pfn+JNlQnGqnpYCAIfDgT179sBqtXJ30DTj8QfxyW4nPt7tgjdwfAVjoUmJ6wdaceFZRshlPF1FRPGRCKelop5QvHbt2tMujGIrfGoqEAhwt+I0Y1TLcdPgDFzV14JP97jw0U9OHHU34bkN1XhvmwPXD7Diou5GKBhyiCgNsLdUitm7d2/kOVL6amgK4bM9LvzfLgdc/uMjOTkGBab0t+LiHiao5Aw5RBQbSTlyAwBOpxN/+9vfsGfPHgBA3759ceutt/INNQGYzWbY7XapyyCJ6ZQyTBlgxRW9zVhe7MI/djlRVd+MlzbWYMl2B6YMsOCSHiaoFDx9SUSpJ+qRmx9++AETJ06EVqvF8OHDAQCbN29GY2MjvvzyS5xzzjkxKbSjpPrIjdvtxu7du2EymdiKgSL8zSGs2OfGBzsdqGs4Pi/OppXj2v5WXHa2CRqGHCLqIIkwchN1uBkzZgy6d++OxYsXQ6E4PvDT3NyM22+/HQcPHsQ333xz+pXHQaqHm+bmZuzcuRMAuFsxnSDQHMKqEg+W7XCgxnu82apFI8c1/Sz4TS8ztEqGHCI6M0kZbrRaLX788ccTNvHbvXs3hg4dioaGhugrjqNUDzcAcODAAdTU1MBqtUpdCiWopqCI1SVuLN3hQFX98ZBjVstwdT8rJvU2Q8eQQ0SnKRHCTdSvYCaTCWVlZSccP3LkCIxGY7R3RzFgMpkSfkk+SUspF3DZ2Wa8eXVn/OG8bOQblXD5Q3hzSx1u+uAQ3t9mR72fv0NElJyiDjdTpkzBbbfdhmXLluHIkSM4cuQIli5dittvvx3XX399LGqkKOn1eqhUKjQ1NUldCiU4hUzAhB4mvHFVJzx4fg4KzUrUB0J490c7bv7wMN79sQ5uhhwiSjJRr5Z6/vnnIQgCbr75ZjQ3Hx/OViqVuPPOO7Fw4cIOL5Cip9VqodFo4PP52GeK2kUuE3BRNyPGdjVgw6F6LNnuwGFnAO9vc+Djn5yY1NuCq/taYNZwkjoRJb7T3uemoaEBBw4cAAB069YtaSavpsOcG+D4acKjR4/CZrNJXQoloZAo4tvDXry/zY5Sx/GWHhqFgN/0MuOafhZYtae1iwQRpYFEmHNz2q9QOp0O/fv3P92bU4wZDAaIoohQKMRWDBQ1mSBgTBcDRnfW4/syL97f7kBJnR8f7nLisz0uXHa2Cdf2tyJDx5BDRImnXa9MV111Fd5++22YTCZcddVVbV73o48+6pDC6MywFQN1BJkgYFRnA0Z20mPz0Qa8v92OvTV+fLzbhc+L3bikhwnX9bcg28DTn0SUONoVbsxmc6RRpslkOqFpJiUetVoNvV4Pt9vNcENnTBAEDC/SY1ihDluPNeL9bXb8VO3DZ3tdWLHPhQk9TJjS34pcI0MOEUmPvaVSWGVlJQ4ePIiMjAypS6EUI4oitlc24v1tDuyobAQAyAXg4u4mTBlgRb6JIYcoXSXCnJuoJ2OMGzcOTqfzpA86bty4aO+OYkin00Emk3HPG+pwgiBgUJ4Oz11agOcvLcA5+VoERWDlfjdu++gwnvumCkddAanLJKI0FfVswHXr1iEQOPFFy+fzYcOGDR1SFHWM8Lwbv9+fNKvZKPn0z9ViQW4Bdlc3Ysk2BzaXN+CrAx58fdCD87sacMNAGzpbVFKXSURppN3hZseOHZF/7969G5WVlZGvg8EgVq5ciYKCgo6tjs6IQqGA2WxGdXU1ww3FXJ9sLZ6eoEVxjQ9Ltjvw/REv1h2sx/qD9TiviwE3DLTiLJta6jKJKA20O9wMGjQIgiBAEISTnn7SarX461//2qHF0ZkzGo2oqKiQugxKI2dnafDk+DyU1PmxZLsd3x72YsOhemw4VI9RnfS4YaAVPTI5yZ2IYqfd4aa0tBSiKOKss87Cpk2bkJWVFfmeSqVCdnY25HLuXppoWrZi4G7FFE/dM9SYMy4PpXY/lmx3YMOhenxX5sV3ZV6MKNLhhoE29MpiyCGijtfucNO5c2cAQCgUilkx1PHYioGk1tWmxqMX5qLMGcD/227HutJ6/PtIA/59pAFDC3S4YaAVfXO0UpdJRCkk6tVSCxYswJtvvnnC8TfffBPPPPNM1AUsWrQIXbp0gUajwYgRI7Bp06Y2r+90OjFjxgzk5eVBrVajZ8+e+OKLL6J+3HQhCAIsFgubaJLkOllUeOiCXCy+shMu7m6ETAB+KG/AzC/K8dDKcuz8eUk5EdGZijrcvPbaa+jVq9cJx/v27YtXX301qvtatmwZZs6ciblz52Lr1q0YOHAgJk6ciOrq6pNePxAI4OKLL8ahQ4fw4Ycfori4GIsXL+ZE5l/RshUDkdQKzSrMGpODN6/ujEt6miAXgG0VjZi1ohwPrDiKbccakGbbbxFRB4t6Ez+NRoM9e/aga9eurY4fPHgQffr0gc/na/d9jRgxAsOGDcPLL78M4Pgpr6KiItxzzz14+OGHT7j+q6++iueeew579+497VMs6bSJX5jf78euXbugUCi4WzElnKr6Jnyw04FV+9xo+jl/98nWYOogK4bk67gjOlGSScpN/IqKivDtt9+ecPzbb79Ffn5+u+8nEAhgy5YtGD9+/H+Kkckwfvx4bNy48aS3+ec//4mRI0dixowZyMnJQb9+/TB//vw2N6nz+/1wu92tLukm3IrB7/dLXQrRCXIMStwzMhtvXdMFV/Q2QykXsLvah0e/rMC9nx/F90e8HMkhoqhEvYnfHXfcgfvuuw9NTU2RJeFr1qzBgw8+iD/84Q/tvp/a2loEg0Hk5OS0Op6Tk4O9e/ee9DYHDx7E119/jalTp+KLL75ASUkJ7rrrLjQ1NWHu3Lknvc2CBQvw5JNPtruuVGWxWGC326Uug+iUsvQK3HVuFqYMsOLDXQ4s3+tGca0fc7+qQHebGlMHWXFuJz1kHMkhol8Rdbh54IEHUFdXh7vuuiuyU7FGo8FDDz2E2bNnd3iBLYVCIWRnZ+P111+HXC7HkCFDUF5ejueee+6U4Wb27NmYOXNm5Gu3242ioqKY1pmIWrZi4JJ9SmQZOgX+Z3gWpvS34sNdTny214USux9Pfl2JrlYVbhhow3ldGHKI6NSiDjeCIOCZZ57B448/jj179kCr1aJHjx5Qq6PbeTQzMxNyuRxVVVWtjldVVSE3N/ekt8nLy4NSqWz15ty7d29UVlYiEAhApTpxi3e1Wh11bamIrRgo2Vi0Ctw+LBPX9rfio5+c+OceJ0odAfxxXSU6WVS4YaAV53cxQC5jyCGi1qKecxNmMBgwbNgw9OvX77TCg0qlwpAhQ7BmzZrIsVAohDVr1mDkyJEnvc3o0aNRUlLSatXPvn37kJeXd9JgQ/8RbsXAeTeUbMwaOaYPycA713bBjYOs0KtkKHMGsHB9FX77cRm+KnEjGOKcHCL6j6jDjdfrxeOPP45Ro0ahe/fuOOuss1pdojFz5kwsXrwY77zzDvbs2YM777wTXq8X06dPBwDcfPPNrU513XnnnbDb7bj33nuxb98+LF++HPPnz8eMGTOifRppyWg0IhgMcnImJSWTWo6bBmfg79d2xrTBNhhUMhx1N+G5DdW47aMyrNrnRjNDDhHhNE5L3X777Vi/fj1uuukm5OXlndEyzSlTpqCmpgZz5sxBZWUlBg0ahJUrV0YmGZeVlUEm+0/+KioqwqpVq3D//fdjwIABKCgowL333ouHHnrotGtIJ+FWDM3NzdytmJKWXiXHDYNsmNzXgs/2uPB/uxyo8DThz99W4/3tdvz3ACvGdzdBJefpKqJ0FfU+NxaLBcuXL8fo0aNjVVNMpeM+N2GiKGLXrl3w+/0wGo1Sl0PUIXxNIXxe7MKHu5xwNB7fFiJTp8CUARZc0sMEleK0z74T0WlIyn1urFYrbDbbaRdH0mErBkpFGqUM1/Sz4p1rOuPOEZnI0MlR29CMRd/XYtqHh/HRT074mrk7N1E6iTrcPPXUU5gzZw4aGhpiUQ/FGFsxUKpSK2SY3MeCt6/ujLvPzUKWXgF7YxCvbarFtH8cxj92OtDYxN97onQQ9WmpwYMH48CBAxBFEV26dDlh7sbWrVs7tMCOls6npQC2YqD00RQUsbrEjaU7HKiqbwYAmNUyXN3Pikm9zdApebqKKBYS4bRU1BOKJ0+efLp1UQIIt2Jwu90MN5TSlHIBl51txoQeJqw54MH/23584vGbW+rwwU4HruprwRW9zTCouaklUaqJeuQm2aX7yA0AVFZW4uDBg8jIyJC6FKK4CYZErCutx5Ltdhx1HZ93plfJMLmPGZP7WGBiyCHqEEk5ckPJj60YKB3JZQIu6mbE2K4GbDhUj/e3O1DmDOD9bQ58/JMTk3pbcHVfC8wa/k0QJbuow41MJmtzb5u2OnRTYtDr9WzFQGlLLhMw9iwjzu9qwLeHvXh/mx2ljgCW7XDg091O/KaXGdf0s8Cq5Wc/omQV9V/vxx9/3OrrpqYm/Pjjj3jnnXfYfTtJyOVymM1mVFVVMdxQ2pIJAsZ0MWB0Zz2+L/Pi/e0OlNT5jzfr3OPCZWebcG1/KzJ0DDlEyabD5twsWbIEy5Ytw6efftoRdxcznHNzXE1NDfbt24eMjIwz2mWaKFWIoojNRxvw/nY79tYc78GmlAu4pIcJ1/W3INvAXb2J2iMR5tx02FrIc889t1UTTEpsLVsxENHxTS6HF+nxwuWFmD8hH32zNWgKivhsrwvT/+8wXvyuGpUeboBJlAw6ZLy1sbERL730EgoKCjri7igOtFottFotfD4f+0wRtSAIAoYU6HBOvhbbKxvx/jYHdlQ24otiN1btc+Pi7iZMGWBFvol/N0SJKupwY7VaW53GEEURHo8HOp0O7733XocWR7EjCALMZjM8Ho/UpRAlJEEQMChPh0F5OuysbMT72+348VgjVu5348sSN8adZcT1A60oNKukLpWIfiHqcPPCCy+0+lomkyErKwsjRoyA1WrtqLooDlq2YmjZfZ2IWuufq8XC3ALsrm7Ekm0ObC5vwFcHPPj6oAcXdDXg+oE2dLYw5BAlinaFm6uuugpvv/02TCYTBEHAlClToFarY10bxZhOp4NarUYgEOBuxUTt0Cdbi6cnaFFc48OS7XZ8f6QBaw/WY93BepzXxYAbBlpxlo2vjURSa9dqKZVKhcOHDyMvLw9yuRwVFRXIzs6OR30djqulWtu7d2/kZ0JE0Smp82PJNju+LfNGjo3upMfvRmRydRWlrURYLdWukZtevXph9uzZuPDCCyGKIj744INT3vHNN98cfcUkGYvFArvdLnUZREmpe4Yacy7KQ6ndjyXbHdhwqB7flnlRYvfj2UsKkGtkwCGSQrtGbr777jvMnDkTBw4cgN1uh9FoPOneKIIgJPwbJUduWnO73di9ezdMJhNbMRCdoUMOP+Z9XYlydxOy9Ao8e0kBV1VR2kmEkZuoN/GTyWSorKzkaakUEQwGsWPHDgDgbsVEHaCuoRkPrizHUVcTMnVyPHtJAQq4oorSSCKEm6iXyJSWliIrK+u0i6PEEm7F4PP5pC6FKCVk6BR47pICdLKoUNsQxKwV5ShzBqQuiyitRB1uOnfuzO36U4zRaEQoFEIHdeIgSns2nQLPXZKPLhYV7I1BPLiyHIccfqnLIkob3NyE2IqBKAYsWgWevbQAZ9lUcDQG8eDKYyi1M+AQxQPDDbVqxUBEHceskeOZSwrQPUMNl+/4CM6BOgYcolhjuKFIK4ZAgPMCiDqaSS3Hwon56JmphtsfwoMry7G/lh8kiGKJ4YYAHG/FAAChUEjiSohSj/HngNM7S436QAgPrzqG4hoGHKJYadcmfoMHD273JOKtW7eeUUEkDb1ez1YMRDGkV8nxxwkFeHz1MfxU7cPDq45h/oR89M7m3xtRR2tXuJk8eXKMyyCpqVQq6PV6uN1uhhuiGNGrZPjjhHw8vvoYdlb58MiX5Xj64nz0zdFKXRpRSol6E79kx038Tq2yshIHDhxAZmam1KUQpTRfUwhzvqrA9spGaBQCnr44H/1zGXAoNSTlJn6UunQ6HeRyOYLBoNSlEKU0jVKGeRfn4Zx8LXzNIh5dfQzbjjVIXRZRyog63ASDQTz//PMYPnw4cnNzYbPZWl0oeYXn3XBJOFHsaRQyPHFRHoYW6OBvFvH4VxXYUs6AQ9QRog43Tz75JP785z9jypQpcLlcmDlzJq666irIZDI88cQTMSiR4iXcisHv5z4cRPGgVsgwd1wuhhfqEAiKmLumApuPeqUuiyjpRR1u3n//fSxevBh/+MMfoFAocP311+ONN97AnDlz8P3338eiRoojk8kEURTZioEoTlQKGR4fl4eRnfRoCop4ck0Fvj/CgEN0JqION5WVlejfvz+A43ujuFwuAMBvfvMbLF++vGOro7jT6XRQKpVsxUAURyq5gEfH5mJ0Zz2aQsBTX1fgu8P1UpdFlLSiDjeFhYWoqKgAAHTr1g1ffvklAGDz5s1Qq9UdWx3FHVsxEElDKRfwyNhcnN/FgOYQ8PTaSmw4xIBDdDqiDjdXXnkl1qxZAwC455578Pjjj6NHjx64+eabceutt3Z4gRRfbMVAJB2FTMDDF+TgwrMMCIrA/HWVWF/qkbosoqTTrk38Wlq4cGHk31OmTEHnzp3x3XffoUePHpg0aVKHFkfSaNmKQSbjbgFE8SSXCXhgTA7kgoCvDniwcH0VgiFgXDej1KURJY2ow43P52u1g+25556Lc889t0OLImmxFQORtOQyATPPy4ZcBqza78FzG6oQFEVc3J0bjxK1R9Qfy7OzszFt2jSsXr2aTRZTVLgVA5eEE0lHLhNw3+hsXNbThJAI/GlDNVbtc0tdFlFSiDrcvPPOO2hoaMAVV1yBgoIC3Hffffjhhx9iURtJyGKxoKmpSeoyiNKaTBBwz6gsTOplhgjgz99W44til9RlESW805pQ/I9//ANVVVWYP38+du/ejXPPPRc9e/bEvHnzYlEjSUCv17MVA1ECkAkCZpybicm9zQCAF7+rwWd7GHCI2nLas0WNRiOmT5+OL7/8Ejt27IBer8eTTz7ZkbWRhHQ6HVsxECUIQRDwuxGZuKqvBQDw8vc1+GS3U9KaiBLZaYcbn8+HDz74AJMnT8Y555wDu92OBx54oCNrIwnJ5XJYLBbOuyFKEIIg4LfDMnBdfwsA4JV/1+L/djmkLYooQUW9WmrVqlVYsmQJPvnkEygUClxzzTX48ssvcf7558eiPpKQ0WhERUUFRFGEIAhSl0OU9gRBwK1DMqCQCViy3YHXN9ehOQRMGWCVujSihBJ1uLnyyivxm9/8Bu+++y4uu+wyKJXKWNRFCSDciqGpqQkqlUrqcogIxwPOtHMyIBcE/H2bHW9uqUNQFHHDQJvUpREljKjDTVVVFYxGbiaVDlq2YmC4IUosNw62QSYD3tlqxztb7QiFgKmDrBxlJUI7w43b7YbJdHzzKFEU4Xafeq+F8PUo+QmCAIvFgrKyMqlLIaKTuGGgDXJBwJtb6vD3bXYERRE3D7Yx4FDaa1e4sVqtqKioQHZ2NiwWy0n/cMLzMrh0OLXo9XoAbMVAlKimDLBCIQNe31yHJdsdaA6JuHVIBgMOpbV2hZuvv/4aNpst8m/+0aQPtmIgSnxX97NCLhPwyr9r8cFOJ5pDwG+HMeBQ+mpXuLngggsi/x47dmysaqEEpFKpYDAY4HK5GG6IEtjkPhbIBQEvf1+Dj35yIhQS8bsRmQw4lJaiPs/Qo0cPPPHEE9i/f38s6qEEZDab2YqBKAlM6m3GvaOyAACf7HFh0fe1CImixFURxV/U4eauu+7C8uXL0atXLwwbNgwvvvgiKisrY1EbJQi2YiBKHpedbcbM0dkQAHy214W/flfDgENpJ+pwc//992Pz5s3Ys2cPLrvsMixatAhFRUWYMGEC3n333VjUSBJjKwai5DKxpwl/GJMNmQB8sc+NF76tRjDEgEPp47SXv/Ts2RNPPvkk9u3bhw0bNqCmpgbTp0/vyNooQbAVA1Hyubi7CQ+MyYFMAFbt9+DP/2LAofQR9SZ+LW3atAlLlizBsmXL4Ha7ce2113ZUXZRg2IqBKPmM62aEXAYsXF+Frw54EBRFPDAmB3IZ/4YptUU9crNv3z7MnTsXPXv2xOjRo7Fnzx4888wzqKqqwtKlS2NRIyWAlq0YiCh5XNDViEfG5kIuAGsP1mPh+io0cwSHUlzUIzfhicQzZszAf//3fyMnJycWdVGCYSsGouQ1posB8gtz8cd1lfjmUD1CooiHL8iFUs4RHEpNUY3cBINBvPbaa1i5ciXuvfdeBps0Em7FEAgEpC6FiE7DqM4GPD4uD0oZ8K/DXvxxXSUCQY7gUGqKKtzI5XLcc889cDqdMSqHElnLVgxElHzOLdJj7kV5UMoFbCzz4qmvKxBo5t8zpZ6o59z069cPBw8ejEUtlOBatmIgouQ0rFCPJy/Kg0ouYNPRBjz5dSX8DDiUYqION08//TRmzZqFzz//HBUVFXC73a0ulLrCrRi4JJwouQ0p0OGp8XlQKwT8UN6AJ9ZUwMeAQylEEMXotq5s2Rm65ZLgZOkK7na7YTab4XK5YDKZpC4n6VRWVuLAgQPIzMyUuhQiOkM7Kxvx2Opj8DWLGJirxbzxedAoT3v7MyIAgMfjgVarRd++fTv0fqN5/476t3jt2rWRy9dffx25hL8+HYsWLUKXLl2g0WgwYsQIbNq0qV23W7p0KQRBwOTJk0/rcSl6bMVAlDr652oxf0I+tAoB238OOo1NHMGh5Bf1UvCWHcI7wrJlyzBz5ky8+uqrGDFiBF544QVMnDgRxcXFyM7OPuXtDh06hFmzZmHMmDEdWg+1rWUrhvAEYyJKXn1ztFgwsQCPfHkMO6t8ePTLY3jq4nzoVRzBoeQV9Wmpb775ps3vn3/++VEVMGLECAwbNgwvv/wygOMrcYqKinDPPffg4YcfPultgsEgzj//fNx6663YsGEDnE4nPvnkk3Y9Hk9LnbnS0lJUVlbCZrNJXQoRdZDiGh8e+fIY6gMh9M5S448T8qFXyaUui5JQIpyWinrkZuzYsSccazn3JprTFYFAAFu2bMHs2bMjx2QyGcaPH4+NGzee8nbz5s1DdnY2brvtNmzYsKHNx/D7/a0mwHLS85ljKwai1HN2lgYLJ+bj4VXHsKfGj4dXHcP8CfkwqhlwKPlEPe7ocDhaXaqrq7Fy5UoMGzYMX375ZVT3VVtbi2AweMJmgDk5OaisrDzpbf71r3/hb3/7GxYvXtyux1iwYAHMZnPkUlRUFFWNdCK2YiBKTT0yNXj2kgKY1DLsqz0ecNw+zq+j5BN1uGkZFMxmMzIzM3HxxRfjmWeewYMPPhiLGiM8Hg9uuukmLF68uN2rdWbPng2XyxW5HDlyJKY1poNwKwYuCSdKPd0y1Hj2kgKYNXKU1Pnx0KpyuBhwKMmcUVfwlnJyclBcXBzVbTIzMyGXy1FVVdXqeFVVFXJzc0+4/oEDB3Do0CFMmjQpciy8W65CoUBxcTG6devW6jZqtRpqtTqquqht4VYMZWVlUpdCRDHQ1abGc5cW4KGV5ThoD+ChleVYODEfFm2HvWUQxVTUv6k7duxo9bUoiqioqMDChQsxaNCgqO5LpVJhyJAhWLNmTWQ5dygUwpo1a3D33XefcP1evXph586drY499thj8Hg8ePHFF3nKKY5atmJoufcREaWGzhYVnr3keMApdQTwwMpjeGZiPmw6BhxKfFH/lg4aNAiCIOCXi6zOPfdcvPnmm1EXMHPmTEybNg1Dhw7F8OHD8cILL8Dr9WL69OkAgJtvvhkFBQVYsGABNBoN+vXr1+r2FosFAE44TrHVshWDRqORuhwiioFOFhWev7QAD64sR5kzgAdWluPZSwqQwYBDCS7q39DS0tJWX8tkMmRlZZ32G9yUKVNQU1ODOXPmoLKyEoMGDcLKlSsjk4zLyso4MpCAwq0YXC4Xww1RCiswq/DcpYV4cGU5jrqa8MCKcjxzSQGy9Aw4lLii3ucm2XGfm47DVgxE6aPS04QHV5ajqr4ZeUYFnr2kANkGpdRlUQJKhH1u2j0ksnHjRnz++eetjr377rvo2rUrsrOz8dvf/parZ9IMWzEQpY9coxLPXVqAPKMCFZ5mzFpRjkoPt4OgxNTucDNv3jz89NNPka937tyJ2267DePHj8fDDz+Mzz77DAsWLIhJkZSYWrZiIKLUl2NQ4rlLC5FvVKKqvhkPrCjHMTcDDiWedoebbdu24aKLLop8vXTpUowYMQKLFy/GzJkz8dJLL+GDDz6ISZGUmORyOSwWCwKBgNSlEFGcZOkVeP6yAhSalKj2NuOBFUdR7uJrACWWdocbh8PRaifh9evX49JLL418PWzYMG6Ql4aMRiNCodAJq+eIKHVl6BR47tICdDIrUdsQxKwV5TjCgEMJpN3hJicnJ7JSKhAIYOvWrTj33HMj3/d4PFAqObks3bAVA1F6sv0ccLpYVLA3BvHAinIcdjLgUGJod7i57LLL8PDDD2PDhg2YPXs2dDodxowZE/n+jh07TtgdmFIfWzEQpS+LVoFnLy3AWTYVHD8HnEMOvhaQ9Nq9UcFTTz2Fq666ChdccAEMBgPeeecdqFSqyPfffPNNTJgwISZFUuJiKwai9GbWyPHMxALMXnUMJXY/HlhRjntGZkOnEqQujSTS2OCH1SBDxy4Ej07U+9y4XC4YDAbI5fJWx+12OwwGQ6vAk4i4z03Hczgc2LNnD6xWKzdcJEpTHn8Qj3x5DPtqOXJDQK8sNVb+YXyH3mc0799RbzFpNptPetxms0V7V5Qi2IqBiIxqORZOzMei72tx2MG5N+ksGAqiwCjtDtbcP5vOGFsxEBEA6FVyPHh+zq9fkVJaeIdiKfEcAnUIs9mM5uZmqcsgIiJiuKGOodfrIZPJ2IqBiIgkx3BDHUKn00Gj0bAVAxERSY7hhjqEXC6H2WxmKwYiIpIcww11GLZiICKiRMBwQx1Gr9ezFQMREUmO4YY6jEajgdFohMfjgdPpZMghIiJJcJ8b6jCCIKB79+5wOp2oqamBy+VCKBSKTDYWBG7HTkREscdwQx1KoVAgMzMTGRkZ8Hg8sNvtqKurQ11dHdRqNXQ63QmtO4iIiDoSww3FhCAIMJlMMJlMyM3NbTWaIwgCdDod1Gq11GUSEVEKYrihmNNoNMjNzUVWVhbcbjdqa2vhdDojW3RrtVo23CQiog7DcENxI5fLYbVaYbFY0NDQAIfDgZqaGjgcDigUCuj1eigU/JUkIqIzw3cSijtBEKDX66HX65GTkxM5ZeV2uyGKIrRaLScgExHRaWO4IUkplUpkZWUhMzMTHo8HdXV1kUnInIBMRESng+GGEkLLCch5eXlwOp2orq6OTEDW6/VQqVRSl0lEREmA4YYSTssJyC6XC7W1tXC5XHC73ZyATEREv4rhhhKWXC6HzWaD1WpFQ0MD7HY7amtrYbfboVQqOQGZiIhOiu8MlPBaTkBuuWeO2+1GKBSCXq+HWq3mBGQiIgLAcENJJjwBueUOyHa7HfX19ZyATEREABhuKEnJZDKYzWaYzWbk5eVFTlk5nU7IZDJOQCYiSmMMN5T0NBoN8vPzkZOT02oCssfjgUaj4QRkIqI0w3BDKaPlBGSv1xvZL4cTkImI0gtf6SnlCIIAg8EAg8HQagKyx+OBKIqRpp2cgExElJoYbiilqVQqZGdnt9oB2eFwRCYg6/V6nrIiIkoxDDeUFlpOQG5sbITD4UBtbS0cDgcnIBMRpRiGG0o74V2Os7OzT5iAHP4eT1kRESUvhhtKWwqFAhkZGbDZbK0mINfV1XECMhFREuMrN6U9TkAmIkotDDdELbScgOx2u1tNQNZoNNDpdJyATESU4BhuiE5CJpPBYrHAYrG0moDsdDojIz1KpVLqMomI6CQYboh+xakmIDc3N3MCMhFRAmK4IWqnlhOQ6+vrI6M5dXV1UKlU0Ol0nIBMRJQA+EpMFCVBEGA0GmE0GjkBmYgoATHcEJ2Bk01AdjqdqK+vh0KhgFarhUqlYtAhIoojhhuiDtByArLP54Pb7YbD4YDb7YbH44FSqYROp+MkZCKiOGC4IepgGo0GGo0GWVlZaGxsjPS08nq9cLlcnJ9DRBRjfHUlihFBEKDT6aDT6ZCdnY2GhobIqav6+no0NzdDrVZDq9Uy6BARdSC+ohLFgSAI0Ov10Ov1yMnJgdfrjQQdj8eDUCgUCTpyuVzqcomIkhrDDVGcyWSyyGqrvLw81NfXw+VywW63w+VyQRRFaDQaaLVa7oZMRHQaGG6IJCSTyWAymWAymZCfnw+PxxMJOg6HA4IgQKvVQq1WM+gQEbUTww1RgpDL5ZEVVwUFBXC73XA6nXA6nXA4HJDJZJGgw6XlRESnxnBDlIAUCgVsNhtsNhsCgQA8Hg8cDgdcLhe8Xi/kcjn30CEiOgWGG6IEp1KpkJGRgYyMDPj9frjdbtjtdng8nsgeOuGgQ0REDDdESUWtViMrKwuZmZnw+XzweDyRoON2u6FSqaDVarlZIBGlNYYboiQUnmis1WqRlZWFhoaGVpsFNjc3R4IO99AhonTDVz2iJHeyPXTCK664hw4RpSOGG6IUIggCDAYDDAZDZA+d8Bwdt9uNUCjEPXSIKOUx3BClqJZ76LTcLDDcuRz4Tx8sBh0iSiUMN0RpQC6Xw2w2w2w2RzYL/OUeOuGgw6XlRJTsGG6I0oxCoYDVaoXVakVTU1OrzQLtdjtkMhl0Oh330CGipJUQY9GLFi1Cly5doNFoMGLECGzatOmU1128eDHGjBkTeXEeP358m9cnolNTKpXIyMhAt27d0K9fP/Ts2TOycWD49FUgEJC6TCKiqEgebpYtW4aZM2di7ty52Lp1KwYOHIiJEyeiurr6pNdft24drr/+eqxduxYbN25EUVERJkyYgPLy8jhXTpRa1Go1MjMz0bNnT/Tr1w89evSA2WyGz+dDXV0dXC4XmpqapC6TiOhXCaIoilIWMGLECAwbNgwvv/wyACAUCqGoqAj33HMPHn744V+9fTAYhNVqxcsvv4ybb775hO/7/X74/f7I1263G0VFRXC5XDCZTB33RIhSkCiKaGxsjOyhU19fj6ampsjScu6hQ0S/5PF4oNVq0bdv3w69X7fbDbPZ3K73b0lHbgKBALZs2YLx48dHjslkMowfPx4bN25s1300NDSgqakJNpvtpN9fsGBBZCKl2WxGUVFRh9ROlA4EQYBOp0NOTg569+6NPn36oGvXrtBoNKivr48EnmAwKHWpREQRkoab2tpaBINB5OTktDqek5ODysrKdt3HQw89hPz8/FYBqaXZs2fD5XJFLkeOHDnjuonSUXgPnfz8fPTp0wd9+vRBp06doFQq4Xa7I7sjM+gQkdSSekx54cKFWLp0KdatWweNRnPS66jVaqjV6jhXRpTaZDIZjEYjjEZjqz10wp3LRVGEVqvlHjpEJAlJw01mZibkcjmqqqpaHa+qqkJubm6bt33++eexcOFCfPXVVxgwYEAsyySiNvxyD51w0LHb7XA4HJE+WNxDh4jiRdKPVCqVCkOGDMGaNWsix0KhENasWYORI0ee8nbPPvssnnrqKaxcuRJDhw6NR6lE1A4KhQIWiwWdO3dGv3790KtXL+Tk5CAUCkXCjs/ng8TrGIgoxUl+WmrmzJmYNm0ahg4diuHDh+OFF16A1+vF9OnTAQA333wzCgoKsGDBAgDAM888gzlz5mDJkiXo0qVLZG5OuJ8OESUGpVIJm80W2TfH7XbD4XDA7XbD6/VCLpdDq9Vys0Ai6nCSh5spU6agpqYGc+bMQWVlJQYNGoSVK1dGJhmXlZW1Omf/yiuvIBAI4Jprrml1P3PnzsUTTzwRz9KJqJ1UKhUyMzORmZkJn88Hj8cT6Vru8XigVCojQYeI6ExJvs9NvEWzTp6IYkcURfh8vlYrrQKBAFQqFbRaLZRKpdQlEtFpSIR9biQfuSGi9BSeaKzVapGdnY2GhoZWmwUGg8FI0OFmgUQUDb5iEJHkBEGAXq+HXq9HTk4OvF4vXC4X6urq4PF4EAqFIrsiy+VyqcslogTHcENECSW8WaDBYGi1h47dbofL5QKASNDhHjpEdDIMN0SUsGQyGUwmE0wmE/Lz8+HxeCJBx+l0AgC0Wi3UajWDDhFFMNwQUVKQy+WwWCywWCwoKCiAx+OB0+mE0+mEw+GATCaLBB0uLSdKbww3RJR0FAoFrFYrrFYrAoFAq6BTX18PuVwOnU7HPXSI0hTDDRElNZVKhYyMDGRkZMDv90f20HG73aivr4dCoeAeOkRphuGGiFJGuFFuZmYmGhsbW20W6Ha7uYcOUZpguCGilBTeQycrKwuNjY2tNgtsamqKrLjiHjpEqYd/1USU0gRBgE6ng06ni+yh43a7IyM63EOHKPUw3BBR2mi5h05ubi7q6+sjQcftdiMUCkGj0XAPHaIkx3BDRGmp5R46LTcLdDgckT10NBoNNBoNgw5RkmG4IaK0J5fLYTabYTabI5sFhpeWO51OCIIQCTpcWk6U+BhuiIhaaLmHTlNTEzweT2Q0x263c7NAoiTAcENEdApKpRI2mw02mw2BQAButxsOh+Oke+gw6BAlDoYbIqJ2UKlUyMzMRGZmJnw+X6vNAj0eD5RKJXQ6HffQIUoADDdERFEKz79puVlgeA8dl8sFtVoNvV7PichEEmG4ISI6TS330MnOzkZDQwPcbjdqamrgdDohk8mg1+s5mkMUZww3REQdQBAE6PV66PV6ZGVlweVyRUJOKBSCXq+HRqORukyitMBwQ0TUwRQKBTIyMmCz2SJtH+x2O+rr6yNtIXjKiih2GG6IiGJEEITI/jm5ubmw2+2oqamB3W6HSqWCXq9nyweiGGC4ISKKg5Zzc5xOJ6qrq+FyuTgvhygGGG6IiOJIpVIhOzsbGRkZkXk5LpcLzc3N0Ol03AWZqAMw3BARSUAul8Nms8FqtaK+vh61tbWw2+2oq6uDRqOBTqfjvByi08RwQ0QkIUEQYDQaYTQakZubC4fDgZqaGjgcDigUCuj1eigUfKkmigb/YoiIEkR4JVVWVhacTidqamrgdrsBHJ+zo1arJa6QKDkw3BARJRilUomsrCxkZGTA7XajtrYWDoej1VJyzsshOjWGGyKiBCWTyWCxWGA2m+H1emG321FbW4u6ujq2eCBqA8MNEVGCEwQBBoMBBoMBOTk5cDgcqK6uhsPhgFwu51Jyol9guCEiSiJqtRq5ubnIzMxs1eJBFEXo9XrOyyECww0RUVL6ZYuHlvNywkvJOS+H0hXDDRFREmvZ4qGhoSHS4qGuro4tHihtMdwQEaWIli0ewvvluFyuyJwdzsuhdMFwQ0SUYlQqFXJyclrNywm3eAjPy+EpK0plDDdERCnqVC0eWs7L4VJySkUMN0REKe6XLR7C++WwxQOlKv42ExGlEa1Wi4KCAmRnZ8PpdKK6uhputxuCIECv10OlUkldItEZY7ghIkpDLVs8uFwu1NbWwul0wuPxQKfTQaPRcF4OJS2GGyKiNCaTyWC1WmGxWNjigVIGww0REf1qiweDwcB5OZQ0+JtKRESttGzx4HQ6I0vJ2eKBkgXDDRERnZRCoUBmZiYyMjIiLR5+uZSc83IoETHcEBFRm1q2eMjNzY3sfswWD5SoGG6IiKjd9Ho99Ho9WzxQQmO4ISKiqJ2sxYPT6UQwGOTpqjTX3NwsdQkMN0REdPpatnjweDyoq6uLhBxKT0qlUvJJ5ww3RER0xgRBgMlkgslkQnNzM0RRlLokkpDUc7AYboiIqENxPxySGredJCIiopTCcENEREQpheGGiIiIUgrDDREREaUUhhsiIiJKKQw3RERElFIYboiIiCilMNwQERFRSmG4ISIiopTCcENEREQpheGGiIiIUgrDDREREaUUhhsiIiJKKWnXulUURQCA2+2WuBIiIiJqr/D7dvh9vC1pF248Hg8AoKioSOJKiIiIKFoejwdms7nN6whieyJQCgmFQjh27BiMRiMEQejQ+3a73SgqKsKRI0dgMpk69L6TQbo/f4A/Az7/9H7+AH8G6f78gdj9DERRhMfjQX5+PmSytmfVpN3IjUwmQ2FhYUwfw2Qype0vNcDnD/BnwOef3s8f4M8g3Z8/EJufwa+N2IRxQjERERGlFIYbIiIiSikMNx1IrVZj7ty5UKvVUpciiXR//gB/Bnz+6f38Af4M0v35A4nxM0i7CcVERESU2jhyQ0RERCmF4YaIiIhSCsMNERERpRSGGyIiIkopDDcdZNGiRejSpQs0Gg1GjBiBTZs2SV1S3HzzzTeYNGkS8vPzIQgCPvnkE6lLiqsFCxZg2LBhMBqNyM7OxuTJk1FcXCx1WXH1yiuvYMCAAZFNu0aOHIkVK1ZIXZZkFi5cCEEQcN9990ldStw88cQTEASh1aVXr15SlxVX5eXluPHGG5GRkQGtVov+/fvjhx9+kLqsuOjSpcsJ//8FQcCMGTMkqYfhpgMsW7YMM2fOxNy5c7F161YMHDgQEydORHV1tdSlxYXX68XAgQOxaNEiqUuRxPr16zFjxgx8//33WL16NZqamjBhwgR4vV6pS4ubwsJCLFy4EFu2bMEPP/yAcePG4YorrsBPP/0kdWlxt3nzZrz22msYMGCA1KXEXd++fVFRURG5/Otf/5K6pLhxOBwYPXo0lEolVqxYgd27d+NPf/oTrFar1KXFxebNm1v9v1+9ejUA4Nprr5WmIJHO2PDhw8UZM2ZEvg4Gg2J+fr64YMECCauSBgDx448/lroMSVVXV4sAxPXr10tdiqSsVqv4xhtvSF1GXHk8HrFHjx7i6tWrxQsuuEC89957pS4pbubOnSsOHDhQ6jIk89BDD4nnnXee1GUkjHvvvVfs1q2bGAqFJHl8jtycoUAggC1btmD8+PGRYzKZDOPHj8fGjRslrIyk4nK5AAA2m03iSqQRDAaxdOlSeL1ejBw5Uupy4mrGjBm4/PLLW70epJP9+/cjPz8fZ511FqZOnYqysjKpS4qbf/7znxg6dCiuvfZaZGdnY/DgwVi8eLHUZUkiEAjgvffew6233trhDarbi+HmDNXW1iIYDCInJ6fV8ZycHFRWVkpUFUklFArhvvvuw+jRo9GvXz+py4mrnTt3wmAwQK1W43e/+x0+/vhj9OnTR+qy4mbp0qXYunUrFixYIHUpkhgxYgTefvttrFy5Eq+88gpKS0sxZswYeDweqUuLi4MHD+KVV15Bjx49sGrVKtx55534/e9/j3feeUfq0uLuk08+gdPpxC233CJZDWnXFZwolmbMmIFdu3al1VyDsLPPPhvbtm2Dy+XChx9+iGnTpmH9+vVpEXCOHDmCe++9F6tXr4ZGo5G6HElceumlkX8PGDAAI0aMQOfOnfHBBx/gtttuk7Cy+AiFQhg6dCjmz58PABg8eDB27dqFV199FdOmTZO4uvj629/+hksvvRT5+fmS1cCRmzOUmZkJuVyOqqqqVserqqqQm5srUVUkhbvvvhuff/451q5di8LCQqnLiTuVSoXu3btjyJAhWLBgAQYOHIgXX3xR6rLiYsuWLaiursY555wDhUIBhUKB9evX46WXXoJCoUAwGJS6xLizWCzo2bMnSkpKpC4lLvLy8k4I8r17906rU3MAcPjwYXz11Ve4/fbbJa2D4eYMqVQqDBkyBGvWrIkcC4VCWLNmTdrNN0hXoiji7rvvxscff4yvv/4aXbt2lbqkhBAKheD3+6UuIy4uuugi7Ny5E9u2bYtchg4diqlTp2Lbtm2Qy+VSlxh39fX1OHDgAPLy8qQuJS5Gjx59whYQ+/btQ+fOnSWqSBpvvfUWsrOzcfnll0taB09LdYCZM2di2rRpGDp0KIYPH44XXngBXq8X06dPl7q0uKivr2/16ay0tBTbtm2DzWZDp06dJKwsPmbMmIElS5bg008/hdFojMy1MpvN0Gq1ElcXH7Nnz8all16KTp06wePxYMmSJVi3bh1WrVoldWlxYTQaT5hjpdfrkZGRkTZzr2bNmoVJkyahc+fOOHbsGObOnQu5XI7rr79e6tLi4v7778eoUaMwf/58XHfdddi0aRNef/11vP7661KXFjehUAhvvfUWpk2bBoVC4nghyRqtFPTXv/5V7NSpk6hSqcThw4eL33//vdQlxc3atWtFACdcpk2bJnVpcXGy5w5AfOutt6QuLW5uvfVWsXPnzqJKpRKzsrLEiy66SPzyyy+lLktS6bYUfMqUKWJeXp6oUqnEgoICccqUKWJJSYnUZcXVZ599Jvbr109Uq9Vir169xNdff13qkuJq1apVIgCxuLhY6lJEQRRFUZpYRURERNTxOOeGiIiIUgrDDREREaUUhhsiIiJKKQw3RERElFIYboiIiCilMNwQERFRSmG4ISIiopTCcENEREQpheGGiOJu3bp1EAQBTqczpo/z9ttvw2KxRL5+4oknMGjQoJg+JhFJj+GGiGJu7NixuO+++yJfjxo1ChUVFTCbzXGtY9asWa2a3BJRamLjTCKKO5VKhdzc3Lg/rsFggMFgiPvjElF8ceSGiGLqlltuwfr16/Hiiy9CEAQIgoC333671Wmp8Omjzz//HGeffTZ0Oh2uueYaNDQ04J133kGXLl1gtVrx+9//HsFgMHLffr8fs2bNQkFBAfR6PUaMGIF169adspZfnpa65ZZbMHnyZDz//PPIy8tDRkYGZsyYgaamptN+DCKSHkduiCimXnzxRezbtw/9+vXDvHnzAAA//fTTCddraGjASy+9hKVLl8Lj8eCqq67ClVdeCYvFgi+++AIHDx7E1VdfjdGjR2PKlCkAgLvvvhu7d+/G0qVLkZ+fj48//hiXXHIJdu7ciR49erSrvrVr1yIvLw9r165FSUkJpkyZgkGDBuGOO+7osMcgovhiuCGimDKbzVCpVNDpdJFTUXv37j3hek1NTXjllVfQrVs3AMA111yDv//976iqqoLBYECfPn1w4YUXYu3atZgyZQrKysrw1ltvoaysDPn5+QCOz6lZuXIl3nrrLcyfP79d9VmtVrz88suQy+Xo1asXLr/8cqxZswZ33HFHhz0GEcUXww0RJQSdThcJNgCQk5ODLl26tJojk5OTg+rqagDAzp07EQwG0bNnz1b34/f7kZGR0e7H7du3L+RyeeTrvLw87Ny5s0Mfg4jii+GGiBKCUqls9bUgCCc9FgqFAAD19fWQy+XYsmVLq3ACIKpJw/F4DCKKL4YbIoo5lUrVaiJwRxg8eDCCwSCqq6sxZsyYDr3veD4GEXU8rpYiopjr0qUL/v3vf+PQoUOora2NjIyciZ49e2Lq1Km4+eab8dFHH6G0tBSbNm3CggULsHz58g6oOj6PQUQdj+GGiGJu1qxZkMvl6NOnD7KyslBWVtYh9/vWW2/h5ptvxh/+8AecffbZmDx5MjZv3oxOnTp1yP3H6zGIqGMJoiiKUhdBRERE1FE4ckNEREQpheGGiIiIUgrDDREREaUUhhsiIiJKKQw3RERElFIYboiIiCilMNwQERFRSmG4ISIiopTCcENEREQpheGGiIiIUgrDDREREaWU/w+uGVmjQoC83wAAAABJRU5ErkJggg==\n" + }, + "metadata": {} + } + ], + "source": [ + "ts = ci.index\n", + "low, high = np.transpose(ci.values)\n", + "\n", + "plt.fill_between(ts, low, high, color='gray', alpha=0.3)\n", + "kmf.survival_function_.plot(ax=plt.gca())\n", + "plt.ylabel('Survival function');" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "FhDF_F2YTTl9" + }, + "source": [ + "Part of [Survival Analysis in Python](https://allendowney.github.io/SurvivalAnalysisPython/)\n", + "\n", + "Allen B. Downey\n", + "\n", + "[Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "YG7xWmS_TTl-" + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "celltoolbar": "Tags", + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11" + }, + "colab": { + "provenance": [], + "include_colab_link": true + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/CHAMOIS_ONTOLOGY_processing_script.ipynb b/CHAMOIS_ONTOLOGY_processing_script.ipynb new file mode 100644 index 0000000..ad7080f --- /dev/null +++ b/CHAMOIS_ONTOLOGY_processing_script.ipynb @@ -0,0 +1,631 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "provenance": [], + "authorship_tag": "ABX9TyP0nAQijjnyl13zl5mbjomA", + "include_colab_link": true + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + }, + "language_info": { + "name": "python" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "view-in-github", + "colab_type": "text" + }, + "source": [ + "\"Open" + ] + }, + { + "cell_type": "code", + "source": [ + "from google.colab import files\n", + "uploaded = files.upload()" + ], + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 73 + }, + "id": "2jHM3dUBlw8L", + "outputId": "e197276a-9a2c-4f3b-fa65-7f9be2102724" + }, + "execution_count": 11, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "" + ], + "text/html": [ + "\n", + " \n", + " \n", + " Upload widget is only available when the cell has been executed in the\n", + " current browser session. Please rerun this cell to enable.\n", + " \n", + " " + ] + }, + "metadata": {} + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Saving chamois-appendix-f-hazard-and-ontology lists.csv to chamois-appendix-f-hazard-and-ontology lists.csv\n" + ] + } + ] + }, + { + "cell_type": "code", + "source": [ + "print (uploaded['chamois-appendix-f-hazard-and-ontology lists.csv'][:200].decode('utf-8') + '...')" + ], + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "9lzlYcL-l_So", + "outputId": "75482dc0-4b22-491e-8be9-2d153512fc31" + }, + "execution_count": 12, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "id,value\r\n", + "Infrastructure.Signalling.Interlocking/ProtectionAWS/TPWS,1\r\n", + "Infrastructure.Signalling.TrainDetection,1\r\n", + "Infrastructure.Signalling.Linesidesignals,1\r\n", + "Infrastructure.Signalling.Runningsignals...\n" + ] + } + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "jqfwEU3mf5-x", + "outputId": "a38bc4d8-12e3-4ad5-f258-60ee46ffccec" + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Drive already mounted at /content/gdrive; to attempt to forcibly remount, call drive.mount(\"/content/gdrive\", force_remount=True).\n" + ] + } + ], + "source": [ + "import csv\n", + "import json\n", + "import io\n", + "import pandas as pd\n", + "\n", + "from google.colab import drive\n", + "drive.mount('/content/gdrive')\n", + "\n", + "\n", + "def insert_path(root, path, value):\n", + " parts = path.split(\".\")\n", + " node = root\n", + "\n", + " for part in parts[:-1]:\n", + " found = None\n", + "\n", + " if \"children\" not in node:\n", + " node[\"children\"] = []\n", + "\n", + " for child in node[\"children\"]:\n", + " if child[\"name\"] == part:\n", + " found = child\n", + " break\n", + "\n", + " if not found:\n", + " found = {\"name\": part, \"children\": []}\n", + " node[\"children\"].append(found)\n", + "\n", + " node = found\n", + "\n", + " # leaf node\n", + " leaf = {\n", + " \"name\": parts[-1],\n", + " \"value\": int(value)\n", + " }\n", + "\n", + " if \"children\" not in node:\n", + " node[\"children\"] = []\n", + "\n", + " node[\"children\"].append(leaf)\n", + "\n", + "\n", + "def csv_to_flare(csv_file, json_file):\n", + " root = {\"name\": \"Train Accident\", \"children\": []}\n", + "\n", + " with open(csv_file, newline=\"\") as f:\n", + " reader = csv.DictReader(f)\n", + "\n", + " for row in reader:\n", + " insert_path(root, row[\"id\"], row[\"value\"])\n", + "\n", + " with open(json_file, \"w\") as f:\n", + " json.dump(root, f, indent=4)\n", + "\n", + "\n", + "# run conversion\n", + "csv_to_flare('chamois-appendix-f-hazard-and-ontology lists.csv', 'gdrive/My Drive/flare.json')" + ] + }, + { + "metadata": { + "id": "CJ9ijZC3Q1Xl", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 423 + }, + "outputId": "2d365dbe-2a17-42e5-b4ca-8ec04786b044" + }, + "source": [ + "df = pd.read_csv(io.StringIO(uploaded['chamois-appendix-f-hazard-and-ontology lists.csv'].decode('utf-8')))\n", + "df" + ], + "cell_type": "code", + "execution_count": 16, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id value\n", + "0 Infrastructure.Signalling.Interlocking/Protect... 1\n", + "1 Infrastructure.Signalling.TrainDetection 1\n", + "2 Infrastructure.Signalling.Linesidesignals 1\n", + "3 Infrastructure.Signalling.Runningsignalsandsub... 1\n", + "4 Infrastructure.Signalling.Colourlightsignals 1\n", + ".. ... ...\n", + "86 People.Signallingstaff.Signaltechnician 5\n", + "87 People.Signallingstaff.Pilotman 5\n", + "88 People.Stationstaff.Dispatchstaff 5\n", + "89 People.Depotstaff(includingstablingpoints/yard... 14\n", + "90 Organisation.Securingcooperationandcompetence.... 2\n", + "\n", + "[91 rows x 2 columns]" + ], + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idvalue
0Infrastructure.Signalling.Interlocking/Protect...1
1Infrastructure.Signalling.TrainDetection1
2Infrastructure.Signalling.Linesidesignals1
3Infrastructure.Signalling.Runningsignalsandsub...1
4Infrastructure.Signalling.Colourlightsignals1
.........
86People.Signallingstaff.Signaltechnician5
87People.Signallingstaff.Pilotman5
88People.Stationstaff.Dispatchstaff5
89People.Depotstaff(includingstablingpoints/yard...14
90Organisation.Securingcooperationandcompetence....2
\n", + "

91 rows × 2 columns

\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + "\n", + "
\n", + "
\n" + ], + "application/vnd.google.colaboratory.intrinsic+json": { + "type": "dataframe", + "variable_name": "df", + "summary": "{\n \"name\": \"df\",\n \"rows\": 91,\n \"fields\": [\n {\n \"column\": \"id\",\n \"properties\": {\n \"dtype\": \"string\",\n \"num_unique_values\": 91,\n \"samples\": [\n \"RailwayVehicles/RollingStock.Brakingsystem.Emergencybrakeequipment\",\n \"Infrastructure.LevelCrossings.Signagefortraindrivers:Levelcrossingsightingboard\",\n \"Operations.Operatinginservice.Inplatformdutiesforpassengertrains(includingcheckstoppingposition/doorrelease/trainsafetycheck/traindispatch)\"\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n },\n {\n \"column\": \"value\",\n \"properties\": {\n \"dtype\": \"number\",\n \"std\": 4,\n \"min\": 1,\n \"max\": 19,\n \"num_unique_values\": 13,\n \"samples\": [\n 15,\n 10,\n 1\n ],\n \"semantic_type\": \"\",\n \"description\": \"\"\n }\n }\n ]\n}" + } + }, + "metadata": {}, + "execution_count": 16 + } + ] + }, + { + "cell_type": "markdown", + "source": [], + "metadata": { + "id": "aj2OooXQlv7T" + } + } + ] +} \ No newline at end of file