{
"cells": [
{
"cell_type": "markdown",
"id": "34ec9328-355b-46dc-a240-ebace87ac8ca",
"metadata": {},
"source": [
"#### Survivalpredict general walkthrough"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "77602431-2f2e-480e-8d0d-621d3971e941",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('ggplot')\n",
"\n",
"%matplotlib inline\n",
"\n",
"from survivalpredict.estimators import CoxProportionalHazard, ParametricDiscreteTimePH, KaplanMeierSurvivalEstimator, KNeighborsSurvival, CoxNeuralNetPH\n",
"from survivalpredict.strata_preprocessing import StrataBuilderDiscretizer, StrataBuilderEncoder\n",
"from survivalpredict.metrics import brier_scores_administrative, integrated_brier_score_administrative\n",
"from survivalpredict.validation import sur_cross_val_score\n",
"from survivalpredict.model_selection import Sur_GridSearchCV\n",
"from survivalpredict.datasets import load_iranian_telecom_churn\n",
"\n",
"from sklearn.preprocessing import StandardScaler"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ba38571c-4299-4314-bdb6-babc035cfeb9",
"metadata": {},
"outputs": [],
"source": [
"#loading some stock data\n",
"iranian_telecom_churn = load_iranian_telecom_churn()\n",
"\n",
"#X is our design matrix/features\n",
"X_raw = iranian_telecom_churn['X']\n",
"ss = StandardScaler()\n",
"X = ss.fit_transform(X_raw)\n",
"\n",
"times = iranian_telecom_churn['times'].astype(np.int64)\n",
"events = iranian_telecom_churn['events'].astype(np.bool_)\n",
"column_names = iranian_telecom_churn[\"column_names\"]"
]
},
{
"cell_type": "markdown",
"id": "e670f71f-3c0a-4016-8ab7-8136ecf9426d",
"metadata": {},
"source": [
"**A quick note on time**\n",
"\n",
"The times array should be the last known interval of survival, regardless of whether the individual experienced the event (i.e., death, churn, conversion) or is still 'alive' or in an unknown state. The times array is assumed to possess the type of integer for 'survivalpredict'. It is up to the user to encode the time array. It is recommended, if possible, to maximize the times array to a few thousand for significant datasets. Presence of large times can trigger a lot of computation for various estimators. It is common to map time to the age of an entity, like the length of time of an individual as a customer."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "ec5e2f4c-66dc-4248-9854-c44b44f7f59c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([38, 39, 37, ..., 18, 11, 11], shape=(3150,))"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"times"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a4d10825-9b19-4fe6-9c55-b11a11b447c2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.int64(47)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"times.max()"
]
},
{
"cell_type": "markdown",
"id": "58401410-72c9-44de-84b6-b554a2a2646f",
"metadata": {},
"source": [
"**Events**\n",
"\n",
"The ‘events’ array should indicate if an individual has experienced the event of interest. Meaning that if said individual has churned, failed, converted or whatever else, at the time interval; the row is coded as ‘True’, and ‘False’ if otherwise. It is assumed to be boolean, or castable to boolean type."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "5cb51c27-f9b1-4251-beb0-7449aec12982",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([False, False, False, ..., False, False, True], shape=(3150,))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"events"
]
},
{
"cell_type": "markdown",
"id": "302d3288-3e93-4bce-bb2f-acbe22f2cfa1",
"metadata": {},
"source": [
"The estimators of survivalpredict are very much like scikit-learn's. But on fit,instead of 'estimator.fit(X,y)' we call 'estimator.fit(X,times,events)'"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "3a45a721-77e2-43aa-ba52-e222f43c73df",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
CoxProportionalHazard()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
\n",
"
\n",
" \n",
" Parameters\n",
"
\n",
" \n",
" \n",
"
\n",
"
\n",
"
alpha
\n",
"
0.0
\n",
"
\n",
" \n",
"\n",
"
\n",
"
\n",
"
max_iter
\n",
"
100
\n",
"
\n",
" \n",
"\n",
"
\n",
"
\n",
"
ties
\n",
"
'breslow'
\n",
"
\n",
" \n",
"\n",
"
\n",
"
\n",
"
tol
\n",
"
1e-09
\n",
"
\n",
" \n",
" \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"CoxProportionalHazard()"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cox = CoxProportionalHazard()\n",
"\n",
"cox.fit(X,times,events)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a317e4f7-3deb-417b-8026-ed113863ab17",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'call_failure': 0.4656662503827164,\n",
" 'complains': 1.8829341579938057,\n",
" 'charge_amount': -0.40937529187142296,\n",
" 'seconds_of_use': 0.29837689062637685,\n",
" 'frequency_of_use': -2.173154255802621,\n",
" 'frequency_of_sms': -3.0279468588329532,\n",
" 'distinct_called_numbers': -0.29707817534951325,\n",
" 'age_group': -0.3465183450012194,\n",
" 'tariff_plan': 0.134561966604143,\n",
" 'status': -0.06233669084308485,\n",
" 'age': 0.268829547578736,\n",
" 'customer_value': 1.87481967428526}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#getting coef for each feature\n",
"dict(zip(map(str,column_names), map(float,cox.coef_)))"
]
},
{
"cell_type": "markdown",
"id": "a67cc107-2db8-49f4-84a3-a5a18ce1e41e",
"metadata": {},
"source": [
"When we call '.predict' on an estimator, we will get the predicted survival curve for each individual. Each row in the predicted array corresponds to the row we ran predict on. Each column represents a point in time, starting with the ‘1’ interval of the ‘times’ array we trained on. By default, the left-most column goes till the max time seen in the training data, but we can set the max time by the max_time key word argument on predict."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "9e1ccdb5-554d-40de-a4cf-42f925581a3c",
"metadata": {},
"outputs": [],
"source": [
"max_time=times.max()\n",
"preds = cox.predict(X,max_time=max_time)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "e9a9ae9a-bb87-4a3d-b86d-e7a2bea98e8c",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n_failed_to_show = 5\n",
"n_censored_to_show = 5\n",
"n_rand_rows_to_show = n_failed_to_show + n_censored_to_show\n",
"\n",
"random_rows_failed = np.random.choice( np.arange(0,preds.shape[0] )[events.astype(np.bool)] ,size=n_failed_to_show, replace=False)\n",
"random_rows_censored = np.random.choice( np.arange(0,preds.shape[0])[~events.astype(np.bool)] ,size=n_censored_to_show, replace=False)\n",
"n_rand_rows_to_show = np.concat((random_rows_failed,random_rows_censored))\n",
"np.random.shuffle(n_rand_rows_to_show)\n",
"\n",
"\n",
"fig, axs = plt.subplots(len(n_rand_rows_to_show))\n",
"fig.set_figheight(13)\n",
"fig.set_figwidth(7)\n",
"\n",
"for a,i in enumerate(n_rand_rows_to_show):\n",
" axs[a].set_ylim(0.0,1.05)\n",
" axs[a].plot(preds[i],color='b')\n",
"\n",
" if events[i]:\n",
" axs[a].axvline(x = times[i], color = 'r', label = 'event')\n",
" else:\n",
" axs[a].axvline(x = times[i], color = 'g', label = 'event',linestyle='--')\n",
"\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "31eea7cb-a0a5-447a-9a75-d03950786c5c",
"metadata": {},
"source": [
"Within survival analysis literature, it is not uncommon to use some form of the ‘brier score’ to assess the performance of a model’s predicted survival curves. Lower the values the better."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "89e9d9f0-313b-406b-84eb-5907d565f02c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG0CAYAAADU2ObLAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAR31JREFUeJzt3Xl8VPW9//H3mS17CGELiyTEEBdIRCxocQHER0Ubq1a0im21KO2t3tr2tqU/0Wv1KrZY24dWsLe9rtSVBxUUBCoiLoBLBQUiCmpAgRBJIJOQbTLL+f0xziSBBLJMZk5mXs/HI4/MnDkz+Wa+kXn7/X6+52uYpmkKAAAgTtli3QAAAIDeRNgBAABxjbADAADiGmEHAADENcIOAACIa4QdAAAQ1wg7AAAgrhF2AABAXCPsAACAuOaIdQOsorq6Wj6fr9PnDxo0SJWVlb3YInQG/WAN9IM10A/WQD9Eh8PhUP/+/Tt3bi+3pc/w+Xzyer2dOtcwjPBz2G0jdugHa6AfrIF+sAb6wZqYxgIAAHGNsAMAAOIaYQcAAMQ1wg4AAIhrhB0AABDXCDsAACCuEXYAAEBcI+wAAIC4RtgBAABxjbADAADimiW3i1i9erWWL18ut9ut3NxczZo1SwUFBe2e+/rrr+vhhx9uc8zpdOrpp5+ORlMBAIDFWS7sbNy4UYsWLdLs2bM1evRovfzyy5o3b54eeOAB9evXr93npKSk6MEHH4xySwEAQF9gubCzYsUKTZs2TVOnTpUkzZ49W5s3b9a6det02WWXtfscwzCUlZXVqdf3er1tNvw0DEMpKSnh250ROq+z56N30A/WQD9YA/1gDbHqh+ZmyeWK6o/sUywVdnw+n8rKytqEGpvNpqKiIu3cubPD5zU1Nemmm26SaZoaNWqUrrnmGp1wwgntnrt06VItWbIkfH/UqFGaP3++Bg0a1OX25uTkdPk5iDz6wRroB2ugH6whmv3wzDPSD38oLVokzZwZtR/bp1gq7NTW1ioQCBw1SpOVlaXy8vJ2nzNs2DD99Kc/VW5urhoaGvTSSy/p9ttv15///GcNGDDgqPMvv/xylZSUhO+H0ndlZaV8Pl+n2mkYhnJyclRRUSHTNDv52yHS6AdroB+sgX6whmj3g98v3XbbIPn9Dq1eXa+pU2t7/WdahcPh6PRAhaXCTncUFhaqsLCwzf1f/vKXWrNmja6++uqjznc6nXI6ne2+Vlf/ME3T5B8VC6AfrIF+sAb6wRqi1Q+vvpqk3buDH+WHDxv0fQcstfQ8MzNTNptNbre7zXG3293pmhyHw6FRo0apoqIi8g0EAMBCHn00PXz78GFLfaRbiqXeGYfDofz8fJWWloaPBQIBlZaWthm9OZZAIKAvv/xS/fv3761mAgAQcx9/7ND69Unh+/X1FKd3xHLTWCUlJVq4cKHy8/NVUFCglStXyuPxaMqUKZKkBQsWKDs7WzO/rsJasmSJRo8erZycHNXX1+ull15SZWWlpk2bFsPfAgCA3vXYY2mSpIED/aqqsjOycwyWCzuTJk1SbW2tFi9eLLfbrby8PM2dOzc8jVVVVdVmSV9dXZ3+9re/ye12Ky0tTfn5+brnnns0YsSIGP0GAAD0rkOHbHrhhVRJ0k9+Uq958zJVV8fITkcMk2omScHVWK2vv3MshmFo6NCh2r9/P8VgMUQ/WAP9YA30gzVEqx8eeihdf/hDpsaObdaDD7o1bdpgZWf7tW3bV732M63G6XR2ejUWY14AAPQhXq/0xBPBKawbbqhXRkYwVNXX85HeEd4ZAAD6kJUrk1VRYdfAgX5demmj0tMDkiSPx5DHE+PGWRRhBwCAPiS03PyHP2xQUpKUltYyXcboTvt4VwAA6CM+/NCpTZtccjpN/eAH9ZIkh0NKSQmO7hw+TJFyewg7AAD0EY8+GqzV+c53GjV4cCB8PFS3w4qs9hF2AADoAyoqbFq+PEWSdOON9W0eS08PhR0+1tvDuwIAQB+waFGavF5DEyZ4VFzc9lIpoSJlprHaR9gBAMDimpqkp54KXkTwhhvqj3o8NLLDlhHtI+wAAGBxL76YooMH7Ro2zKeLLmo66vGMjNDIDh/r7eFdAQDAwkxTeuSR4HLz669vkKOdjZ5Cy88pUG4fYQcAAAt75x2Xtm93Kjk5oJkzj57CklpWYzGy0z7eFQAALCy03PyKKxrVv3/7+22FprEY2WkfYQcAAIvas8euf/0rWdLRy81ba1l6TthpD2EHAACLevzxNAUChs47r0mFhb4Oz2tZes7Hent4VwAAsKD6ekPPPtvxcvPWWHp+bIQdAAAs6OWXk1Vba1Nenk/nn3/s7cwpUD423hUAACxo9epgrc53v9so23E+rdPSKFA+FsIOAAAW09Bg6I03gmHnwgsbj3t+y8gOYac9hB0AACzmjTeS1NRk6IQTfBozpuPC5JBQgXJ9PR/r7eFdAQDAYlatCo7qTJ/eJKMTgzWtR3bM9i/Fk9AIOwAAWIjXK61d2xJ2OiO0GisQMNTYyFTWkQg7AABYyDvvuOR225Sd7deECc2dek5qqinD4MKCHSHsAABgIatXp0iSvvWtJtntnXuOYVCkfCyEHQAALMI0W5acd3YKK6Rl53M+2o/EOwIAgEVs2eJURYVdqakBnXvusS8keKTQZqCM7ByNsAMAgEWERnWmTvUoOblrz23ZMoKP9iPxjgAAYBGhsHPRRV2bwpIY2TkWwg4AABbw2Wd2ffqpUw6HqfPP73rYaanZIewcibADAIAF/OtfwVVYZ5/tUb9+Xb8yIJuBdox3BAAAC2h91eTuCG0ZwcjO0Qg7AADEWEWFTR984JIUvL5Od4RGdlh6fjTeEQAAYuxf/wqO6owf36ycnEC3XiM0skOB8tEIOwAAxFgo7HR3CktqvfScsHMkwg4AADFUU2Now4YkSdL06Y3dfh0KlDvGOwIAQAy99lqyfD5DhYVenXiiv9uvk5ZGgXJHCDsAAMRQaBXWhRd2fwpLYmTnWHhHAACIkaYmad264BRWd66a3FqoQJmanaMRdgAAiJG33kpSQ4NNQ4f6VVzs7dFrtYzsEHaORNgBACBGQnthTZ/eKKOHGSW0GquhwSZ/90t/4hJhBwCAGPD7pVde6fmS85DQNJbEVNaRCDsAAMTAv//t0qFDdmVlBXTmmc09fr2kJMnlYiqrPYQdAABiIDSFdcEFTXI6I/OaLcvP+XhvjXcDAIAoM83W9To9n8IKoUi5fYQdAACibPt2h/bscSg5OaApUzwRe92WLSP4eG+NdwMAgChbvTpFkjR5skcpKWbEXjcjg81A20PYAQAgynpjCkuS0tKCwYktI9oi7AAAEEV79ti1fbtTdrupCy6IbNgJjexQoNwW7wYAAFG0dm1we4gJE5qVnR25KSyppWaHaay2CDsAAETR2rXBKaxp0yJXmBwSCjuM7LTFuwEAQJQ0NhrauDE4sjNtWmSnsCQKlDtC2AEAIErWr3epqcnQiBE+FRb6Iv76LUvPCTutEXYAAIiS1lNYPd34sz0tIzt8vLfGuwEAQBSYZktxcm9MYUksPe8IYQcAgCj45BOHysuDV02eNCnyxclS6+0i+HhvjXcDAIAoCE1hnX12s1JSeudnpKcHp7Go2WmLsAMAQBT09hSWxHV2OuKIdQPas3r1ai1fvlxut1u5ubmaNWuWCgoKjvu8DRs26MEHH9Q3vvENzZkzJwotBQDg+KqrDb3/vkuSdMEFvTOFJbVMY3GdnbYs925s3LhRixYt0owZMzR//nzl5uZq3rx5qqmpOebzDhw4oH/84x865ZRTotRSAAA65403khUIGDr5ZK+GD/f32s8JTWM1Nxvy9F6m6nMsF3ZWrFihadOmaerUqRoxYoRmz54tl8uldevWdficQCCghx56SFdddZUGDx4cxdYCAHB80ZjCklqmsSSpvt5yH/ExY6lpLJ/Pp7KyMl122WXhYzabTUVFRdq5c2eHz1uyZIkyMzN1/vnn6+OPPz7mz/B6vfJ6veH7hmEo5etKMaOTFz0IndfZ89E76AdroB+sgX6whvb6we+X1q0LFidfcIGnV/vI4ZBSUgJqbLSprs6mAQMiu/dWX2WpsFNbW6tAIKCsrKw2x7OyslReXt7ucz755BO99tpruu+++zr1M5YuXaolS5aE748aNUrz58/XoEGDutzenJycLj8HkUc/WAP9YA30gzW07oeNG6Xqaql/f6mkZKAcvfzJ26+f1NgoJScP1tChvfuz+gpLhZ2uamxs1EMPPaSf/OQnyszM7NRzLr/8cpWUlITvhxJ2ZWWlfL7OXbrbMAzl5OSooqJCpklqjhX6wRroB2ugH6yhvX54/vkMSek677xGVVa6e70NqamDJDm0e3eVhgzxHvf8vsrhcHR6oMJSYSczM1M2m01ut7vNcbfbfdRojyR99dVXqqys1Pz588PHQn9cV199tR544IGj/i/H6XTK6XS2+/O7+g+EaZr8o2IB9IM10A/WQD9YQ+t+ePXVlnqdaPRNqEi5ttbgb+Frlgo7DodD+fn5Ki0t1cSJEyUFi49LS0s1ffr0o84fNmyY7r///jbHnnvuOTU1Nen666/XwIEDo9JuAADaU15u0/btThmGqalTo7M8KlSkzPLzFpYKO5JUUlKihQsXKj8/XwUFBVq5cqU8Ho+mTJkiSVqwYIGys7M1c+ZMuVwujRw5ss3z09LSJOmo4wAARNtrrwULk8eP9yo7OxCVnxnaDJT9sVpYLuxMmjRJtbW1Wrx4sdxut/Ly8jR37tzwNFZVVRWrDQAAfUK0lpy31jKyw2dliOXCjiRNnz693WkrSbrzzjuP+dybb765F1oEAEDXNDVJb70Vu7DDZqAteCcAAOgF77yTpMZGm3Jy/BozpnOrfSMhNI3F/lgtCDsAAPSC1lNY0ay+CI3ssPN5C8IOAAARZprS2rXB4uRp06K7SVVo6TnTWC14JwAAiLDPPrPriy8ccrlMnXNOtMMOBcpHIuwAABBhoVGdb37To7S06F7YLyOD6+wciXcCAIAIa6nXie6ojtQyjcXITgvCDgAAEVRTI737rktSdJech7QsPSfshBB2AACIoDVrJJ/P0IknepWX54/6z28Z2eEjPoR3AgCACHr55eD3WExhSa1rdgyxD2gQYQcAgAgJBKSVK4O3YzGFJbVMYwUChhobmcqSCDsAAETM1q1OHTgQnEqaOLE5Jm1ITTVlGNTttEbYAQAgQl59NbgKa/Jkj1yu2LTBMFqmsgg7QYQdAAAiJJZLzlsLFSnX1/MxLxF2AACIiLo6Q9u2OSVJU6bEOuwwstMaYQcAgAjYssWpQMDQyJFSTk4gpm1p2TKCj3mJsAMAQER88EGwSOfMM2PcEEkZGVxFuTXCDgAAEfDBB8EpLCuEHTYDbYuwAwBAD5mmtUZ2Wmp2+JiXCDsAAPRYeblNX31ll91uavz4WLeGzUCPRNgBAKCHQqM6p5ziU2pqjBuj1ltG8DEvEXYAAOixUNgZPz42V00+Umhkh6XnQYQdAAB6aPPmYHHy6ad7Y9ySIAqU2yLsAADQA15vcE8syTojO0xjtcW7AABAD+zY4VBTk02ZmQGdeKI/1s2RRIHykQg7AAD0wObNwXqdceOaZbPIpypLz9viXQAAoAdCxclWqdeRGNk5EmEHAIAeCF05+fTTrVGvI7Wu2SHsSIQdAAC6rabG0KefhoqTrTSyEww7DQ02+a1RRhRThB0AALppy5bgFNbIkT4NGBDbnc5bC01jSYzuSIQdAAC6LXR9HassOQ9JSpJcLqayQgg7AAB0kxWLk0NaipT5qOcdAACgG4I7nVuvODmkZfk5IzuEHQAAumHPHrsOHrTL6TQ1ZowVR3a4inII7wAAAN0QGtUZM8ar5OQYN6YdGRlcayeEsAMAQDeErpxsxSksSUpLo0A5hLADAEA3tIQd601hSS0jO2wZQdgBAKDLmpuljz6ybnGy1Lpmh5Edwg4AAF20fbtTHo+hrKyARo2y5iWKW7aM4KOedwAAgC4KFSePH98sw6IDJ6Hr7LD0nLADAECXWb04WWIaqzXCDgAAXWTlKyeHtCw956OedwAAgC44dMjQrl0OSdK4cdYd2WHpeQvCDgAAXfDhh8FRnVGjfOrf34xxazoWKlBm6TlhBwCALmmZwrLuqI7UeiNQRnYIOwAAdEFoJdYZZ1g77LQsPSfsEHYAAOik4E7n1i9OltgItDXeAQAAOmnXLrvcbpuSkkydcorVw05wGqu52ZDHE+PGxBhhBwCATgqN6owd65XLFePGHEdoZEeS6usT++M+sX97AAC6IFSvY/XiZEmy26WUFK6iLBF2AADotNCVk8ePt37YkVovPyfsAACA42hqCm4AKlm/ODmEIuWgxP7tAQDopNJSp7xeQwMG+HXCCdbc6fxILVtGMLIDAACOI1ScPH6817I7nR+pZcuIxP64T+zfHgCATupLxckhoZEdanYAAMBx9ZVtIloL1ezU1xN2AADAMVRV2fTllw4Zhqlx4/pGcbLUEnYSfTNQR6wb0J7Vq1dr+fLlcrvdys3N1axZs1RQUNDuue+++66WLl2qiooK+f1+5eTk6JJLLtF5550X5VYDAOLV5s3BKayCAp8yM6270/mRmMYKslzY2bhxoxYtWqTZs2dr9OjRevnllzVv3jw98MAD6tev31Hnp6en67vf/a6GDRsmh8OhzZs36+GHH1ZmZqbGjRsX/V8AABB3+sp+WEdi6XmQ5X77FStWaNq0aZo6dapGjBih2bNny+Vyad26de2eP2bMGE2cOFEjRoxQTk6OLr74YuXm5uqTTz6JcssBAPHq/ff7Xr2OxNLzEEuN7Ph8PpWVlemyyy4LH7PZbCoqKtLOnTuP+3zTNFVaWqry8nJde+217Z7j9Xrl9bYkc8MwlJKSEr7dGaHzOns+egf9YA30gzXQD72nvt4Ih51zzmk+5ntstX5ITw9+r6uzWaZNsWCpsFNbW6tAIKCsrKw2x7OyslReXt7h8xoaGvSTn/xEPp9PNptNN9xwg4qLi9s9d+nSpVqyZEn4/qhRozR//nwNGjSoy+3Nycnp8nMQefSDNdAP1kA/RN6KFVJzs5SXJ5199uBOXWPHKv0wcmTwu8eTpKFDh8a2MTFkqbDTXcnJyfrjH/+opqYmbdu2TYsWLdKQIUM0ZsyYo869/PLLVVJSEr4fSrqVlZXy+Xyd+nmGYSgnJ0cVFRUyzb5TqBZv6AdroB+sgX7oPS+8kCkpTeedV6+Kitpjnmu1fvB6XZIGqLraq/37q2LdnIhyOBydHqiwVNjJzMyUzWaT2+1uc9ztdh812tOazWYLp+i8vDzt27dPy5YtazfsOJ1OOZ3Odl+nq3+Ypmla4o850dEP1kA/WAP9EHmvv54kSZo82dPp99Yq/ZCW1rIaywrtiRVLFSg7HA7l5+ertLQ0fCwQCKi0tFSFhYWdfp1AINCmLgcAgO744gu7du1yyOEwdfbZnlg3p8vS00MFypb6uI86S43sSFJJSYkWLlyo/Px8FRQUaOXKlfJ4PJoyZYokacGCBcrOztbMmTMlBWtwTjzxRA0ZMkRer1cffPCB3nrrLd14440x/C0AAPEgNKrzjW80KyOj742MhNpcV2fINNVn9vSKNMuFnUmTJqm2tlaLFy+W2+1WXl6e5s6dG57GqqqqalNR7vF49Mgjj+jgwYNyuVwaPny4fvazn2nSpEkx+g0AAPGi9RRWXxQKO4GAocZGQ6mpfS+wRYJhJvIkXiuVlZWdnvoyDENDhw7V/v37E3oONNboB2ugH6yBfoi85mZp7Ngc1dfbtHp1pYqKjv8ZYbV+ME1p5MihCgQMbd5coSFDArFuUsQ4nc5OFygn9iQeAAAd2LTJpfp6mwYM8GvMmL5ZB2oYrffHStA5LBF2AABoV+spLFsf/rSkSJmwAwBAu0JhZ8qUvlmvE9K6SDlREXYAADhCZaVNpaXBLSL6anFySFoam4F2ezVWY2Oj6uvrNXDgwPCxQ4cOac2aNfJ6vTrrrLNUUFAQkUYCABBNb7wRHNUpKmrWwIF9u6g3tBloItfsdDvs/O1vf1NlZaXmzZsnKbg/1W233aZDhw7JMAytWrVKc+fObfcqxgAAWFko7PT1UR2ppUC5vj5xw063x7R27Nih8ePHh++/9dZbqq6u1t13363HH39cI0eO1AsvvBCRRgIAEC2BQEu9ztSp8RN2Dh9O3Gmsbv/mtbW1ys7ODt9///33dfLJJ6uwsFApKSmaPHmydu/eHYk2AgAQNaWlTh06ZFd6ekBnnNEc6+b0WMtqLEZ2uiwtLS28YWdzc7M++eQTFRcXt7ywzabm5r7/RwIASCyhUZ2zz/aog32j+5TQaqxEHtnpds1OYWGhXnnlFQ0fPlwffvihmpubNWHChPDj+/fvbzPyAwBAXxAvS85DQgXKjOx0w/e//33Z7Xb96U9/0tq1a1VSUqITTjhBUnDX8XfeeUennHJKxBoKAEBvq601tGlTcMl5vISdlqXniRt2uj2yk5OTowceeEB79+5VamqqBg8eHH7M4/Fo1qxZys3NjUgjAQCIhg0bkuTzGcrP92nkSH+smxMRTGP1cNdzh8OhvLy8o46npKS0mdICAKAvaJnCaopxSyInVKCcyEvPexR2Ghoa9Morr+ijjz5STU2NfvzjH6ugoEB1dXV6/fXX9Y1vfEM5OTmRaisAAL3GNOPr+johLD3vQdg5ePCg7rzzTlVVVWno0KHat2+fmpqCSTg9PV1r1qxRZWWlfvSjH0WssQAA9JbPP7drzx6HXC5TkybFz2pilp73IOz84x//UGNjo/74xz8qMzNTs2fPbvP4hAkTtHnz5h43EACAaHjjjWRJ0sSJzUpNNWPcmshpqdlJ3LDT7TGtrVu36qKLLtKIESNkGEe/gUOGDNHBgwd71DgAAKKl5arJ8VOvI7WEncZGm/zxUXPdZd0OO83NzcrMzOzw8cbGxu6+NAAAUdXUJG3cGB+7nB8pLa1lI9NEncrqdtgZMWKEPv744w4f//e//93uSi0AAKzmvfdcamqyKSfHr5NP9sW6ORGVlCS5XKFr7SRmkXK3f+uLL75YGzZs0LJly9TQ0CApeDHBiooKPfTQQ9q5c6e+/e1vR6yhAAD0ltdfD9brTJ7sUTuVGX1eohcpd7tA+bzzzlNVVZWef/55Pffcc5Kke++9V6Zpymaz6ZprrtHEiRMj1lAAAHpLy5Lz+KrXCUlPN3XoUOIWKffoOjvf/e53dd555+mdd95RRUWFTNPUkCFDdOaZZ2rIkCGRaiMAAL2mvNymTz5xymYzde658VWvExK61k6iTmN1K+x4PB7dcccdmjZtmr71rW+ppKQk0u0CACAq3nwzOKpz2mleZWfHz5Lz1kKbgSbqyE63Il5SUpIOHDjQ7pJzAAD6klC9Trxs/Nme0MhOom4Z0e3xrHHjxmnLli2RbAsAAFHl90tvvRV/+2EdKVSgnKhbRnT7t77iiiu0f/9+PfTQQ/rkk0906NAh1dXVHfUFAIBVffihU263Tf36BTRunDfWzek1LTU7iTmy0+0C5V/96leSpL1792r9+vUdnvf8889390cAANCrQquwzj3XI0ePluxYW8uWEYk5stPtrr3iiiuo2QEA9Gmh/bDi7arJR+I6O9101VVXRbIdAABEVW2toQ8+cEqSzjsv3sNOYk9jRWw8q7m5Wc3NzZF6OQAAetXGjUny+w3l5/s0YkR875DZsvScaawuq6qq0uLFi/XBBx+otrZWkpSZmanTTz9dV155pQYNGhSRRgIAEGnxftXk1hJ96Xm3w86+fft0xx13qL6+XsXFxRo+fLgkqby8XG+++aY2bdqku+++W8OGDYtYYwEAiJTQxQTjfQpLagk7jOx00dNPPy3DMHTfffdp5MiRbR778ssvdffdd+vpp5/Wb37zmx43EgCASPryS7t273bI4TD1zW/GfwlGohcodzviffzxx7rooouOCjqSNHLkSF144YXavn17jxoHAEBvCI3qjB/fHF6WHc9alp4TdrrE5/PJ5XJ1+HhSUpJ8Pl93Xx4AgF4TqtdJhCksqXXNTmJOY3X7tx41apRee+01NTQ0HPVYQ0ODXnvtNeXn5/eocQAARJrfH1yJJSVS2AlOYzU3G/Ikxq/cRo+uszNv3jz94he/0JQpU8KFyOXl5XrjjTd0+PBh3XDDDRFrKAAAkbB1a8sWEaedFr9bRLQWGtmRpLo6m5KSAjFsTfR1O+yMHTtWt956q5566im9+OKLbR7Ly8vTf/7nf2rs2LE9biAAAJEUmsI6++z43iKiNbtdSk0NqKHBpro6QwMGxLpF0dWjbi4uLtZ9990nt9utyspKSdKgQYOUlZUVibYBABBxoV3OE2UKKyQ93VRDQ2IWKUck02ZlZRFwAACWV1dn6P33g4trEjHsHDgQnMZKNN3+jVeuXKl58+Z1+Pi9996rV155pbsvDwBAxL39tks+n6G8PJ9yc+N7i4gjtWwZkXgjO90OO+vWrQtfNbk9I0aM0KuvvtrdlwcAIOJC19c599zEGtWREvsqyt3+jSsqKjRixIgOHx82bJi++uqr7r48AAARFwo7kycnXtg54YTgte8+/TRBqrJb6XbYcTgccrvdHT7udrtlGIk3VAYAsKZ9+2z67DOnbDZTkyYlXtgpKgous9+2zRnjlkRft8NOYWGhXn/9dTU2Nh71WENDg9atW6fRo0f3qHEAAERKaBXWuHFe9esX/1tEHCl0TaEtW5wyE+zX73bYmTFjhqqrqzVnzhytWrVKpaWlKi0t1cqVKzVnzhy53W5deeWVkWwrAADd9sYbyZIScwpLkk45xSuHw9ShQ3aVl9tj3Zyo6vbE3ejRo/Xb3/5Wf//73/XEE0+0eWzw4MGaM2eOCgsLe9o+AAB6LBCQ3norMZechyQnSyed5NNHHzm1ZYtTw4cnzmq0Hl9U8C9/+Yt2796tiooKSVJOTg57YgEALKW01KnqarvS0wM6/fTmWDcnZoqLm/XRR05t3erUxRc3xbo5UdPtaazdu3dr/fr1stlsys/P16RJk5Samqonn3xSc+fO1cqVKyPZTgAAui20CmvSJI+ciVefG1ZcHKzb2bo1sd6Eboedp556Shs3bgzfP3DggO6//34dOHBAkvTkk09ynR0AgCWE9sNK1HqdkNZhJ5GKlLsddr744gudfPLJ4ftvvPGGbDab5s+fr3vvvVdnnXWW1qxZE5FGAgDQXQ0NLVtEJOLFBFs75RSvnE5T1dV27d2bOEXK3Q47DQ0NysjICN//4IMPVFxcrMzMTEnBep5QHQ8AALHyzjsuNTcbGjHCp/z8xCnKbU9SknTSSYk3ldXtsJOVlaV9+/ZJkqqrq1VWVqbi4uLw401NTVxUEAAQc6F6nfPO84iPpZbr7SRS2On2aqwJEyZo1apVam5u1meffSan06mJEyeGH//iiy80ZMiQiDQSAIDuah120HIlZcJOJ1x99dWqra3VW2+9pdTUVN10003KysqSFJzieuedd3ThhRdGqp0AAHRZRYVNO3Y4ZRimzj6bsCO1HtlxyTSVEKNd3Q47ycnJuuWWWzp87H//93/lcrm63TAAAHoqNKpz2mleZWcn0PKjYzjppGCRsttt0549do0cGf91TL2yz7vNZlNqaqocjsTbWRUAYB2h/bASfRVWa0lJwVVZUnCfrERgyTSyevVqLV++XG63W7m5uZo1a5YKCgraPffVV1/Vm2++qT179kiS8vPzdc0113R4PgAgMQQCLSM7iX59nSMVFXm1datL27Y5dckl8X8l5V4Z2emJjRs3atGiRZoxY4bmz5+v3NxczZs3TzU1Ne2ev337dp199tn63e9+p3vuuUcDBgzQPffco0OHDkW55QAAK/n4Y4eqquxKTQ3ojDMSd4uI9rSu20kElhvZWbFihaZNm6apU6dKkmbPnq3Nmzdr3bp1uuyyy446/8i6of/4j//Qu+++q23btmny5MlHne/1euX1esP3DcNQSkpK+HZnhM5jaX1s0Q/WQD9YA/1wtDffDO5yPmlSs5KSovO+9JV+CF1Jeds2pyQj7ouULRV2fD6fysrK2oQam82moqIi7dy5s1Ov4fF45PP5lJ6e3u7jS5cu1ZIlS8L3R40apfnz52vQoEFdbm9OTk6Xn4PIox+sgX6wBvqhxbvvBr+XlCRr6NChUf3ZVu+HAQMkl0tyu21qahqqeN+/21Jhp7a2VoFAILyEPSQrK0vl5eWdeo2nn35a2dnZKioqavfxyy+/XCUlJeH7ofRdWVkpn8/XqZ9hGIZycnJUUVEhM5E2F7EY+sEa6AdroB/aamyU3nwzR5Kh00+v1P79nfv3vaf6Uj+ccsoAbdni0po11frOd/pe3Y7D4ej0QIWlwk5PLVu2TBs2bNCdd97Z4bJ3p9MpZwdb3nb1D9M0Tcv/MScC+sEa6AdroB+CNm5MksdjKCfHrxNP9EZ908u+0A/FxV5t2eLS1q0OXXKJtdvaU5YqUM7MzJTNZpPb7W5z3O12HzXac6SXXnpJy5Yt0+23367c3NzeayQAwPKeeipVknTRRY1xX4/SXaG6nS1b4r9I2VJhx+FwKD8/X6WlpeFjgUBApaWlKiws7PB5L774ov75z39q7ty5OvHEE6PRVACARe3ZY9eaNcHi5Ouvb4hxa6yruDi4Qm3bNmfUR76izVJhR5JKSkq0du1avf7669q7d68eeeQReTweTZkyRZK0YMECPfPMM+Hzly1bpueff14//elPNXjwYLndbrndbjU19b35RwBAz/3jH6kKBAyde65HBQXRqdXpi046yaekJFO1tTbt3m2PdXN6leVqdiZNmqTa2lotXrxYbrdbeXl5mjt3bngaq6qqqs2SvjVr1sjn8+nPf/5zm9eZMWOGrrrqqmg2HQAQY42N0jPPBKewfvSj+hi3xtqcTunUU7364AOXtm51atSo+N02wnJhR5KmT5+u6dOnt/vYnXfe2eb+woULo9AiAEBf8NJLKaqutmvECJ8uuIAR/uMpKgqFHZcuvTR+3y/LTWMBANAdpik99liaJOm66xpkj++ZmYg47bRg3c7WrfG9RxZhBwAQFzZtcqq01KXkZFNXX80UVmcUFbVcSTkQiHFjehFhBwAQF554Ijiqc+mljcrOjvPlRRFSWBgsUj58OL6LlAk7AIA+78ABm1asCO5zSGFy54WKlKX43hSUsAMA6POefjpVXq+hM85oDk/NoHNCFxeM57odwg4AoE/zeqWnngpOYTGq03WJUKRM2AEA9GmrViWrosKuQYP8+va3G2PdnD4nEYqUCTsAgD4tVJh87bUN6mAPaBxDYaFPycmm6upsKiuLzyJlwg4AoM/66COH3n03SQ6HqR/8gCms7nA4WoqUt22Lz7RI2AEA9FlPPhkc1bnooibl5MTpHEwUhOp2tmyJz7odwg4AoE9yuw39858sN4+E1nU78YiwAwDok557LlVNTTadcopXEyc2x7o5fdppp8V3kTJhBwDQ5/j90qJFLcvNDSPGDerjCgp8Sk4OqL7eprIyS+4R3iOEHQBAn7NuXZK++MKhfv0Cuvxylpv3lMMhjR3rkxSfdTuEHQBAnxNabv697zUoNZV9sCKhuDh+Ly5I2AEA9CllZXatW5cswzD1wx9SmBwp8bxtBGEHANCnhJabT53q0ahR/hi3Jn6Ewk5pqVP+OHtbCTsAgD6jvt7Q4sWpkqRZsxjViaSCAp9SUgJqaIi/ImXCDgCgz1iyJEW1tTbl5fk0ebIn1s2JK3a7NHZscHQn3oqU4yu6AQDilmm2FCZff329bPzvesQVF3v1738naetWp2bMaH+VW1OTtH27U1u3OuV222SzSTabZLebMoxgaAoeM8OPDR4c0IUXNkX5t2lB2AEA9Anr17u0c6dTqakBfe97DbFuTlw6skjZ75c++8yhDz906sMPXfrwQ6c+/tgpr7drFzYaP76ZsAMAwPGERnVmzGhUZibLzXtDy5WUXZoxY4C2bnWqvv7oIbQBA/waN86rnBy/TFPy+w0FAsFwZJr6+nbwWCAgnXiiL9q/ShuEHQCA5e3da9crryRLYh+s3pSf71NGRkCHD9v09ttJkqTU1ICKi70aN86rceOaNW6cVyNG+PvUVasJOwAAy3vyyVQFAobOOcejwsLYjhLEM7tdeuABt9avd2ns2GDAGT3aJ7s91i3rGcIOAMDSGhulZ54JTmGx3Lz3TZ/epOnTY1df0xuoZQcAWNqyZalyu2064QSfLrggvj6EER2EHQCAZZmm9NhjwVGd666r7/PTKYgNwg4AwLLee8+l7dudSk4O6OqrWW6O7iHsAAAsKzSq893vNqp/f5abo3sIOwAASyovt2nVKpabo+cIOwAAS/rHP9Lk9xs66yyPTj2V5eboPsIOAMByPB7p6aeDu5szqoOeIuwAACxn+fIUHTxo19Ch/ri75guij7ADALCcxx8PFib/8If1cnD5W/QQYQcAYCmbNwd32E5KMnXttSw3R88RdgAAlhIa1fnOdxo1YEAgxq1BPCDsAAAs48ABm5YvT5HEPliIHMIOAMAynn46VV6voTPOaFZxsTfWzUGcIOwAACyhuTl4bR2JUR1EFmEHAGAJq1Yl66uv7Bo82K+LL26MdXMQRwg7AICYM03p0UfTJUk/+EG9XK4YNwhxhbADAIi5NWuStGmTS8nJLDdH5BF2AAAx5fVKd9/dT5J04411GjKE5eaILMIOACCmnnoqVWVlDg0c6Nd//mddrJuDOETYAQDETE2NoT/9KUOS9KtfHVZGhhnjFiEeEXYAADHz0EMZqq62q7DQq5kzqdVB7yDsAABi4ssv7Xr00eB1dW6/vZYNP9FrCDsAgJj4wx8y1Nxs6JxzPDr/fE+sm4M4RtgBAETd5s1OvfhiqgzD1B131MgwYt0ixDPCDgAgqkxTuuuu4FLzq65q1Jgxvhi3CPGOsAMAiKqXX07W+++7lJIS0Jw5tbFuDhIAYQcAEDUej3TvvZmSpJ/+tF45OVxAEL2PsAMAiJonnkjTF184NGSIXz/9KRcQRHQQdgAAUXHokKEHHwxeQHDOnFqlpnIBQUQHYQcAEBUPPJChmhqbTjnFqyuvbIx1c5BACDsAgF5XVmbXk08GLyB4xx21sttj3CAkFMtdr3L16tVavny53G63cnNzNWvWLBUUFLR77p49e/T8889r165dqqys1HXXXadvf/vbUW4xAOB47r03Uz6fofPPb9J553EBQUSXpUZ2Nm7cqEWLFmnGjBmaP3++cnNzNW/ePNXU1LR7vsfj0ZAhQzRz5kxlZWVFt7EAgE555x2XVq1Kkc1m6vbbWWqO6LNU2FmxYoWmTZumqVOnasSIEZo9e7ZcLpfWrVvX7vkFBQX6wQ9+oLPPPltOpzPKrQUAHE9Tk3TbbcELCM6c2aCTTuICgog+y0xj+Xw+lZWV6bLLLgsfs9lsKioq0s6dOyP2c7xer7xeb/i+YRhKSUkJ3+6M0HmdPR+9g36wBvrBGqzaD3/4Q4Y++cSpAQP8mjOnznLtizSr9kOis0zYqa2tVSAQOGo6KisrS+Xl5RH7OUuXLtWSJUvC90eNGqX58+dr0KBBXX6tnJyciLUL3Uc/WAP9YA1W6odXXpH+7/+Ct594wq6ioiGxbVAUWakfYKGwEy2XX365SkpKwvdD6buyslI+X+eGVw3DUE5OjioqKmSaXCciVugHa6AfrMFq/XDwoKEf/GCQJLuuu65eZ5xRq/37Y92q3me1fohnDoej0wMVlgk7mZmZstlscrvdbY673e6IFh87nc4O63u6+odpmiZ/zBZAP1gD/WANVugH05R+85ssHThgV0GBV//937Uxb1O0WaEf0MIyBcoOh0P5+fkqLS0NHwsEAiotLVVhYWEMWwYA6Ipnn03V6tUpcjpNLVxYrZQUPvQRW5YZ2ZGkkpISLVy4UPn5+SooKNDKlSvl8Xg0ZcoUSdKCBQuUnZ2tmTNnSgoWNe/duzd8+9ChQ9q9e7eSk5OZLwWAGPj8c7vuuCO40edvf1ursWNZfYXYs1TYmTRpkmpra7V48WK53W7l5eVp7ty54WmsqqqqNhXuhw4d0pw5c8L3ly9fruXLl+vUU0/VnXfeGeXWA0Bi83qlW27pr8ZGmyZN8ugnP6mPdZMASZJhMqkoKVig3HpJ+rEYhqGhQ4dq//79zMnGEP1gDfSDNVihH+bPz9Bf/pKhfv0CWrPmgIYPD8SkHbFkhX5IFE6ns9MFypap2QEA9F3vvefSggXpkqT5890JGXRgXYQdAECP1NYa+tnPshQIGLryygZdcklTrJsEtEHYAQD0yG239dPevQ6NHOnT3Xe3v5chEEuEHQBAty1blqIXXkiV3W7qoYeqlZFBnQqsh7ADAOiWvXvtuvXW4CafP/95nb7xjc4t8gCijbADAOiyhgZDN97YX7W1No0f36yf//xwrJsEdIiwAwDokkBA+vnPs7Rtm0vZ2X49/HC1HJa6ahvQFmEHANAl99+foZUrU+RymXr00WqdcII/1k0CjomwAwDotKVLU/TggxmSgtfTmTixOcYtAo6PsAMA6JRNm5z61a+yJEk333xYV13VGNsGAZ1E2AEAHNe+fXbdcEO2PB5DF17YqP/3/yhIRt9B2AEAHFN9vaHrr89WZaVdp57q1UMPuWXj0wN9CH+uAIAOBQLSz36Wpe3bnRo40K8nnjiktDQuHIi+hbADAOjQH/6QoX/9K0VJSaYeffSQhg9n5RX6HsIOAKBdixenaOHC4Mqr++93c4Vk9FmEHQDAUd57z6U5c7IkSbfccljf/S4rr9B3EXYAAG188olDN9zQX16voYsvbtRvfsPKK/RthB0AQNirrybpO98ZqEOH7Bo7tlkPPsjKK/R9/AkDAGSa0t/+lqbrr89Wfb1N3/ymR88+e0ipqay8Qt/H1m0AkOCam6Vbb+2n555LkyRde2297rmnRi5XjBsGRAhhBwAS2KFDNs2e3V/vvJMkm83UHXfU6sYb62UYsW4ZEDmEHQBIUDt2OHT99dn68kuHMjICevjhap1/vifWzQIijrADAAlo7dok3XRTf9XV2ZSb69MTTxxSYaEv1s0CegUFygCQQExT+vvfg4XIdXXBQuQVK6oIOohrjOwAQIJoaDD03/+dGS5EnjmzXvPmUYiM+EfYAYAE8N57Lv3yl1navdtBITISDmEHAOJYY6N0332Z+r//S5NpGho61K8//7la553XHOumAVFD2AGAOLVpk1O//GWWPv/cKUm6+up6/e53tcrM5EKBSCyEHQCIM01N0p//nKG//jVdgYChnBy/7rvPrWnTWFaOxETYAYA4smWLUz//eT/t3BkczZkxo0F33VWjrCxGc5C4CDsAEAeam6X//m/p978fIL/f0KBBfs2fX6MLL2yKddOAmCPsAEAf9tlnDv3znyl64YUU7d0rSYYuu6xBd99do+xsRnMAibADAH3OgQM2vfhiMOBs3dpykZxBg6R7763WxRc3xrB1gPUQdgCgD6ivN7R6dbJeeCFFb76ZpEAgeIEch8PUlCkeXXFFo667rr/c7iaZDOgAbRB2AMCi/H5p/fokLVmSolWrktXY2LLDz/jxzbriigZdckmTBgwIyDAMpaRIbnfs2gtYFWEHACxm9267Fi9O1eLFqdq/3x4+npfn0xVXNOjyyxs1apQ/hi0E+hbCDgBYQEODoRUrkrV4carefjspfLxfv4Auu6xRM2Y06PTTvWzvAHQDYQcAYsQ0pfffd+n551P00kspqq8PTlMZhqnJkz266qoGXXhhk5KTY9xQoI8j7ABAFB04YNN777n03nsurVuXrLKyln+G8/J8uuqqBs2Y0aDhwwMxbCUQXwg7ANBLTFMqK7Pr3/926d13k/Teey7t3t32n93U1IBKSpr0ve816Mwzm5mmAnoBYQcAIsDrlb780q5duxz67DOHNm0Kjt5UVdnbnGcYpk45xaeJE5s1caJH06Z5lJ7OWnGgNxF2AKCTfD5p375goCkrc2jXruDtXbsc2rPHLr//6GEZl8vUuHHNX4ebZn3jG83q149wA0QTYQdAQmtslGpqbDp40KaqKrsqK21ff9nD36uqgscOHrTJNDueZ0pJCSgvz69Ro3w67TSvJk5sVnFxMwXGQIwRdgBYmmkGp4h8PkPNzS3fvV5DTU1HfzU2GmpqUvh+fb1NNTWG3G6bamqOvG2Tx9O1IpmkJFN5eT6NGuXTqFH+r78Hv3JyAtTcABZE2AH6KNMMfqA3NLT98vkM+XzBKRe/32jz/chjfr/a3A7eD90OnhcIKPyaoeOh80LPdbmkurqsNo+HzgkE2t5u+d7Rz20bbHy+3k8PNpup/v0DGjQo9OUPfx84MKDBgwMaODB4bODAgGy2478mAOsg7AC9zDSDF4yrqwt+1dfbdPiwofp6Q3V1NtXXtwSV1uGlsdFQQ4PtiPttHz/WlEr0pUTtJzmdphwOU8nJppKTpZSU0O2jv1JTTWVlBdSvX+h78Kv1sfR0kxEZII4RdoDj8Hql6upgvcbBgzZVV9t0+HAwsNTVBb8Hv2yqqzPCjx0+HAwy9fVGeNPG3pKcHFBqqqmUFFMul2S3m3I4JLs9uFFk2+/Bx493LPjV/uvYbMHnBI8Z6t8/U/X1NbLZgufbbC2vGzpmt7e+3f7joeNOpymns/3vDocIJgC6hLCDhOTxSJWVdlVU2PTVV/avv4Jh5tAhmw4etIdv19REZs7CZjOVlmYqPd1UenpwNCEtzVRaWjCohMJKSorZ5n7L7ZbzWn+lpAQDQqwYhqGhQzO1f3+DTLbbBmBBhB3EJY9HKitzaOdOhz791Kl9+4Jh5sCBYMCpru5aOjCMYE3HgAEBZWcHlJFhKjMzGFgyMoL3MzJC91tuh0JNenowlDAiAQDRR9hBn9bUJH30UTDU7Njh0KefBm/v3u1o95onrblcpgYP9mvIkICGDPGHi1AHDGgJNaHbWVmBmI6eAAC6j7ADyzJNqbraUHm5vcOvffukQGBQu8/PzAxo9GifCgu9ys31a/Bgv3JyAl8HHL/692ekBQASAWEHllBZadPWrU5t2+ZUaalTO3c6VF5uV2Pj8etlMjMDKiz06aSTvBo92qeTTvJp9Ggv1zwBAEgi7CDKTFOqqLBp2zantm1zff3dqYqKjueIBgzwa9iwlq/hw0O3A5owYaCkryRRGAsAaB9hBxFjmlJdnaGKCrv27w+ucqqoCH599ZVNFRV27dljP2pjRClYAHziiT4VFXlVVOTVKad4NWKEX0OH+pXSweVbgquApP37gz8bAID2EHbQLbW1hkpLneGpp48+Cq54amg4/rSTzWaqsLAl2BQVeTVmjFdpaSQWAEDkEXZwXDU1Rni6aetWl7ZudWr37o7/dDIzg6ubcnICysnxt/oKaOhQvwoLfUpJIdgAAKLDkmFn9erVWr58udxut3JzczVr1iwVFBR0eP7bb7+t559/XpWVlcrJydG1116r8ePHR7HFfY/fH7wqcGWlTVVVwd2eg99tX+/wHLyo3oEDNpWXt/9nMny4T8XFLaMzubnBjRAZoQEAWInlws7GjRu1aNEizZ49W6NHj9bLL7+sefPm6YEHHlC/fv2OOn/Hjh168MEHNXPmTI0fP17r16/XH//4R82fP18jR46MwW/Qe0wzeLG80H5JoT2VWvZWsn2935Kh2lqbamsN1dTYwrdra4M7PtfW2lRX17WrAo8Y0RJsiouDX9nZgV76TQEAiBzLhZ0VK1Zo2rRpmjp1qiRp9uzZ2rx5s9atW6fLLrvsqPNXrlypcePG6Tvf+Y4k6eqrr9a2bdu0evVq/fjHP45m09sIbUfQGX6/5HYHR1kOHgyOsoRuV1YGR1xCWxdEcgfo0FWBBw5s+Ro0KHhRveDuzsHbo0b5lJ3NaA0AoG+yVNjx+XwqKytrE2psNpuKioq0c+fOdp+zc+dOlZSUtDl22mmn6d///ne753u9Xnm93vB9wzCU8vVyH6OTF2UJnXes8z/6yKlLLhnYqdfrjqSk4L5IoX2VgnssBbckSEszv97ZObilQb9+AWVmhm6bX98P3nZ0+i/Aehes6Uw/oPfRD9ZAP1gD/WBNlgo7tbW1CgQCysrKanM8KytL5eXl7T7H7XYfNb3Vr18/ud3uds9funSplixZEr4/atQozZ8/X4MGtX8V3mPJycnp8LFBg6Tk5M69jmFI2dnS4MHSkCHB761vh74PHChlZkppaZLDYSgYQCKzSWVfdqx+QPTQD9ZAP1gD/WAtlgo70XD55Ze3GQkKpe/Kykr5fL5OvYZhGMrJyVFFRUWHuzyPHCmVlfW8vUdqaAh+oXP9gN5HP1gD/WAN9EP0OByOTg9UWCrsZGZmymazHTUq43a7jxrtCcnKylJNTU2bYzU1NR2e73Q65XQ6232sq3+Ypmnyx2wB9IM10A/WQD9YA/1gLZaaA3E4HMrPz1dpaWn4WCAQUGlpqQoLC9t9TmFhobZt29bm2NatWzV69OhebSsAAOgbLBV2JKmkpERr167V66+/rr179+qRRx6Rx+PRlClTJEkLFizQM888Ez7/4osv1pYtW7R8+XLt27dPixcv1ueff67p06fH6DcAAABWYqlpLEmaNGmSamtrtXjxYrndbuXl5Wnu3Lnhaamqqqo2Ve4nnXSSbrnlFj333HN69tlnNXToUP3mN7+Ju2vsAACA7jFMJhUlBQuUWy9JP5bgBpRDtX//fuZkY4h+sAb6wRroB2ugH6LH6XR2ukDZctNYAAAAkUTYAQAAcY2wAwAA4hphBwAAxDXCDgAAiGuEHQAAENcIOwAAIK4RdgAAQFwj7AAAgLhmue0iYsXh6Ppb0Z3nIPLoB2ugH6yBfrAG+qH3deU9ZrsIAAAQ15jG6obGxkb99re/VWNjY6ybktDoB2ugH6yBfrAG+sGaCDvdYJqmdu3axSZvMUY/WAP9YA30gzXQD9ZE2AEAAHGNsAMAAOIaYacbnE6nZsyYIafTGeumJDT6wRroB2ugH6yBfrAmVmMBAIC4xsgOAACIa4QdAAAQ1wg7AAAgrhF2AABAXGPzjm5YvXq1li9fLrfbrdzcXM2aNUsFBQWxblbc2r59u1566SXt2rVL1dXV+vWvf62JEyeGHzdNU4sXL9batWtVX1+vk08+WTfeeKOGDh0aw1bHl6VLl+q9997Tvn375HK5VFhYqO9///saNmxY+Jzm5mYtWrRIGzdulNfr1WmnnaYbb7xRWVlZsWt4nHnllVf0yiuvqLKyUpI0YsQIzZgxQ6effrok+iBWli1bpmeeeUYXX3yxrr/+ekn0hdUwstNFGzdu1KJFizRjxgzNnz9fubm5mjdvnmpqamLdtLjl8XiUl5enG264od3HX3zxRa1atUqzZ8/Wvffeq6SkJM2bN0/Nzc1Rbmn82r59uy688ELNmzdPt99+u/x+v+655x41NTWFz3nyySe1adMm/dd//ZfuuusuVVdX609/+lMMWx1/srOzNXPmTP3hD3/Q73//e40dO1b33Xef9uzZI4k+iIXPPvtMa9asUW5ubpvj9IXFmOiSW2+91XzkkUfC9/1+v/njH//YXLp0aewalUCuvPJK89133w3fDwQC5uzZs80XX3wxfKy+vt6cOXOmuX79+lg0MSHU1NSYV155pfnRRx+Zphl8z6+++mrz7bffDp+zd+9e88orrzR37NgRq2YmhOuvv95cu3YtfRADjY2N5i233GJu2bLF/N3vfmc+/vjjpmny34MVMbLTBT6fT2VlZSoqKgofs9lsKioq0s6dO2PYssR14MABud1uFRcXh4+lpqaqoKCAPulFDQ0NkqT09HRJUllZmfx+f5v/NoYPH66BAwfSD70kEAhow4YN8ng8KiwspA9i4JFHHtHpp5/e5t8fif8erIianS6ora1VIBA4as41KytL5eXlsWlUgnO73ZKkfv36tTner1+/8GOIrEAgoCeeeEInnXSSRo4cKSnYDw6HQ2lpaW3OpR8i78svv9Rtt90mr9er5ORk/frXv9aIESO0e/du+iCKNmzYoF27dun3v//9UY/x34P1MLIDoEseffRR7dmzR7/4xS9i3ZSENGzYMP3xj3/Uvffeq29961tauHCh9u7dG+tmJZSqqio98cQTuuWWW+RyuWLdHHQCIztdkJmZKZvNdlQyd7vdVNjHSOh9r6mpUf/+/cPHa2pqlJeXF5tGxbFHH31Umzdv1l133aUBAwaEj2dlZcnn86m+vr7N/83W1NTw30aEORwO5eTkSJLy8/P1+eefa+XKlZo0aRJ9ECVlZWWqqanRb3/72/CxQCCgjz/+WKtXr9Ztt91GX1gMYacLHA6H8vPzVVpaGl76HAgEVFpaqunTp8e4dYlp8ODBysrK0rZt28LhpqGhQZ999pm+9a1vxbZxccQ0TT322GN67733dOedd2rw4MFtHs/Pz5fdbte2bdt01llnSZLKy8tVVVWlwsLCWDQ5YQQCAXm9XvogioqKinT//fe3OfbXv/5Vw4YN06WXXqqBAwfSFxZD2OmikpISLVy4UPn5+SooKNDKlSvl8Xg0ZcqUWDctbjU1NamioiJ8/8CBA9q9e7fS09M1cOBAXXzxxXrhhRc0dOhQDR48WM8995z69++vCRMmxLDV8eXRRx/V+vXrNWfOHKWkpIRHN1NTU+VyuZSamqrzzz9fixYtUnp6ulJTU/XYY4+psLCQf9wj6JlnntG4ceM0cOBANTU1af369dq+fbtuu+02+iCKUlJSwvVqIUlJScrIyAgfpy+shV3Pu2H16tV66aWX5Ha7lZeXpx/96EcaPXp0rJsVtz766CPdddddRx2fPHmybr755vBFBV999VU1NDTo5JNP1g033NDmgnfomauuuqrd4zfddFM46IcuorZhwwb5fD4uotYL/vrXv6q0tFTV1dVKTU1Vbm6uLr300vBqIPogdu68807l5eUddVFB+sIaCDsAACCusRoLAADENcIOAACIa4QdAAAQ1wg7AAAgrhF2AABAXCPsAACAuEbYAQAAcY2wAwAA4hphB4Dl3XzzzVq4cGGsmwGgjyLsALCMHTt2aPHixaqvr491UwDEETYCBWAZO3bs0JIlSzRlyhSlpaWFjz/wwAMyDCOGLQPQlxF2AFie0+mMdRMA9GFsBArAEhYvXqwlS5YcdXzBggW66667dOqpp+rmm2+WJL3++ut6+OGH9T//8z/auHGjNmzYIL/fr0mTJmnWrFnyeDx6/PHHtWnTJknStGnTdO2117YZHQoEAlq1apXWrl2rr776SqmpqZowYYJmzpyp9PT08Hmff/65nnvuOZWVlampqUlZWVkaM2aMbrrppl5+RwBECiM7ACzhzDPP1P79+7VhwwZdd911ysjIkCRlZmZ2+JzHHntMWVlZuuqqq/Tpp5/q1VdfVWpqqnbu3KmBAwfqmmuu0ebNm/XSSy/phBNO0OTJk8PP/fvf/6433nhDU6ZM0UUXXaQDBw5o9erV2rVrl+6++245HA7V1NTonnvuUWZmpi699FKlpaWpsrJS7777bq+/HwAih7ADwBJyc3M1atQobdiwQRMmTNDgwYOP+5x+/frp1ltvlWEYuvDCC1VRUaHly5frggsu0OzZsyVJF1xwgW6++WatW7cuHHY++eQTvfbaa7rlllt0zjnnhF9vzJgxuvfee/XOO+/onHPO0Y4dO1RfX6/bb79dJ554Yvi8q6++OsK/PYDexGosAH3W+eef32ZqqqCgQKZp6vzzzw8fs9lsys/P11dffRU+9vbbbys1NVXFxcWqra0Nf+Xn5ys5OVmlpaWSFC6S3rRpk3w+X5R+KwCRxsgOgD5r4MCBbe6npqZKkgYMGHDU8dbL2SsqKtTQ0KAbb7yx3detra2VJJ166qk688wztWTJEr388ssaM2aMJkyYoHPOOYeiaaAPIewA6LNstvYHp9s73notRiAQUL9+/fSzn/2s3eeH6oQMw9CvfvUr7dy5U5s2bdKWLVv017/+VStWrNC8efOUnJwcgd8CQG8j7ACwjGhdS2fIkCHatm2bTj75ZLlcruOeX1hYqMLCQl1zzTVav369/vKXv2jDhg2aNm1aFFoLoKeo2QFgGUlJSZKkhoaGXv05kyZNUiAQaHepu9/vD0951dXV6circ+Tl5UmSvF5vr7YRQOQwsgPAMvLz8yVJzz77rM4++2zZ7XadccYZEf85p556qi644AItW7ZMX3zxhYqLi2W321VRUaG3335bP/rRj3TWWWfpjTfe0CuvvKIJEyYoJydHjY2NWrt2rVJSUjR+/PiItwtA7yDsALCMgoICfe9739OaNWv04YcfyjRNLViwoFd+1o9//GPl5+fr1Vdf1bPPPiu73a5Bgwbp3HPP1UknnSQpGIo+++wzbdy4UTU1NUpNTdWJJ56oW265pVNL4wFYA1dQBgAAcY2aHQAAENcIOwAAIK4RdgAAQFwj7AAAgLhG2AEAAHGNsAMAAOIaYQcAAMQ1wg4AAIhrhB0AABDXCDsAACCuEXYAAEBcI+wAAIC49v8BXPGT3kYP/4UAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"scores = brier_scores_administrative(preds,times,events)\n",
"\n",
"plt.xlabel('times')\n",
"plt.ylabel('scores') \n",
"plt.plot(times_of_survival_curve,scores,color='b')"
]
},
{
"cell_type": "markdown",
"id": "8811b363-4349-4207-9969-3db01ef20384",
"metadata": {},
"source": [
"The integral of the ‘brier scores’ are used as a singular metric."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "8ea2d3ad-948b-4bfa-aef5-0d8d47ec7a24",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.float64(4.225046295847599)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrated_brier_score_administrative(preds,times,events)"
]
},
{
"cell_type": "markdown",
"id": "c4a1b1aa-466f-48cd-97df-ac559b0307dc",
"metadata": {},
"source": [
"Within statistical learning methodologies, cross-validation is considered an important cornerstone of prediction evaluation. 'survivalpredict' comes with some tooling to evaluate models with cross-validation."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "1db8f896-b320-4590-8b6a-7179033e3633",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.float64(3.9861410461140165)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sur_cross_val_score(CoxProportionalHazard(),\n",
" X,\n",
" times,\n",
" events,\n",
" brier_score_max_time=max_time,\n",
" cv=10,\n",
" scoring='integrated_brier_score_administrative').mean()"
]
},
{
"cell_type": "markdown",
"id": "75181cb4-a190-429e-825f-0ac98fcac490",
"metadata": {},
"source": [
"It is a good idea to compare the performance of your models against the Kaplan Meier univariate survival curve. This would be analogous to sklearn’s dummy estimators."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "727fbee3-4dd2-42d8-bd80-ff0bfd83207c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.float64(4.815747737919164)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sur_cross_val_score(KaplanMeierSurvivalEstimator(),\n",
" X,\n",
" times,\n",
" events,\n",
" brier_score_max_time=max_time,\n",
" cv=10,\n",
" scoring='integrated_brier_score_administrative').mean()\n",
"\n",
"#the lower the better, it seems that the cox model has a lower integrated brier score than the dummy KaplanMeier model."
]
},
{
"cell_type": "markdown",
"id": "7297b7b7-ed3f-4fb5-a279-3ef6f533d585",
"metadata": {},
"source": [
"If you wish to train a model with strata, fitting becomes 'estimator.fit(X,times,events, strata=strata)’.\n",
"and predict becomes 'estimator.predict(X, strata=strata)'. The strata array should be an array with integer values. We also have tools for building strata. It is common practice to remove the column/columns used to build the strata."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "e0cf1f73-e6f9-4e69-8cd2-e995d2d6dc25",
"metadata": {},
"outputs": [],
"source": [
"#getting the position of the age column of the numpy array\n",
"position_of_age_col = int(\n",
" np.argwhere(column_names == \"age\")[0][0]\n",
")\n",
"age = X[:, position_of_age_col]\n",
"\n",
"\n",
"\n",
"sbd = StrataBuilderDiscretizer(n_bins=3,strategy='uniform')\n",
"\n",
"strata = sbd.fit_transform(age)\n",
"X_without_strata = X[:,~np.isin(range(X.shape[1]), position_of_age_col)]"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "f37569de-b3db-48a9-af62-92a9cb02843d",
"metadata": {},
"outputs": [],
"source": [
"cox_with_strata = CoxProportionalHazard()\n",
"\n",
"cox_with_strata.fit(X_without_strata,times,events,strata=strata)\n",
"_ = cox_with_strata.predict(X_without_strata,strata=strata)\n"
]
},
{
"cell_type": "markdown",
"id": "51dff29c-9ded-4ba9-9c06-df498a2ba26f",
"metadata": {},
"source": [
"We can see a significant improvement in the cross-validation score after adding the strata to our cox model."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "523c9b71-be52-483e-81be-02af228d195a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"np.float64(2.10385484255526)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sur_cross_val_score(CoxProportionalHazard(),\n",
" X_without_strata,\n",
" times,\n",
" events,\n",
" brier_score_max_time=max_time,\n",
" cv=10,\n",
" strata = strata,\n",
" scoring='integrated_brier_score_administrative'\n",
" ).mean()"
]
},
{
"cell_type": "markdown",
"id": "bfff1a73-8891-414a-aac7-b5ded4d63dc9",
"metadata": {},
"source": [
"we can also make our strata from categorical data."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "7349ba53-209c-4f18-9681-b04ca03c7fd5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 1, 2])"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"StrataBuilderEncoder().fit_transform(['a','b','b','c'])"
]
},
{
"cell_type": "markdown",
"id": "b864511f-07d2-48ea-8e8a-5554709b2ea4",
"metadata": {},
"source": [
"**On Base Hazards**"
]
},
{
"cell_type": "markdown",
"id": "4717c0ed-549c-4d28-bad6-a5890c556edd",
"metadata": {},
"source": [
"The Cox Proportional Hazards model famously avoids estimating the base hazard. Rather, it uses Breslow’s non-parametric estimator for relative risk. The Breslow estimator is simply a function of the sum of failures and the total relative risk at each interval of time. This can lead to stepwise and sporadic base hazard. See here."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "638626ff-19f5-45fb-9af1-27d0277f69b6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.00000000e+00, 0.00000000e+00, 4.22362623e-06, 2.11185727e-06,\n",
" 3.23049888e-06, 4.36391890e-06, 1.00938202e-05, 5.73938244e-06,\n",
" 1.06450882e-05, 5.99679788e-06, 1.13141900e-05, 5.06075672e-06,\n",
" 3.81145763e-06, 2.55188701e-06, 2.57106780e-06, 1.28562813e-06,\n",
" 1.28623916e-06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
" 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.29690188e-06,\n",
" 2.62215977e-06, 2.62261275e-06, 6.70512152e-06, 1.10353441e-05,\n",
" 1.39089749e-05, 2.81904663e-05, 3.69440115e-05, 4.53947482e-05,\n",
" 7.88527319e-05, 1.29201916e-04, 1.57957811e-04, 2.07469988e-04,\n",
" 2.38680946e-04, 3.23737833e-04, 3.37293381e-04, 3.19367161e-04,\n",
" 3.58466069e-04, 4.16940175e-04, 8.15280093e-04, 9.12927357e-02,\n",
" 5.87992106e-02, 0.00000000e+00, 0.00000000e+00])"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#notice the 0's at different points in time\n",
"cox._breslow_base_hazard"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "740537ca-4ad9-4ac2-ba1b-87ace4cb6457",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 0.0001)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAGiCAYAAAAC4AllAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAS7VJREFUeJzt3Xt8VNW9///XnkxCCCEJCDkJBhIjhGgNIBqwYCsXaxGpKNAWb21NSbXWLz3f/mrLQc73tB4Bqe051dJTc0o4lrYIkYIIBqogtNVg5VIuUSCGCCRcJAiTkARCJrN+f0yydQ5BMrnNJe/n4+Ejmdlrdj7Jysg7a6+1tmWMMYiIiIgIjkAXICIiIhIsFIxEREREmigYiYiIiDRRMBIRERFpomAkIiIi0kTBSERERKSJgpGIiIhIEwUjERERkSYKRiIiIiJNFIxEREREmjjb8qKNGzeybt06XC4Xqamp5OTkMHjw4Mu237ZtGytXrqSyspKkpCQeeOABRo4caR83xlBQUMDmzZupra0lMzOTWbNmkZycbLepqalh6dKl7Ny5E8uyGD16NA8//DDR0dEAXLx4kd/+9reUlZVx7NgxRo4cyY9+9KNLannvvfdYtmwZ5eXlXHXVVUyfPp1x48a15ccgIiIiYcbvEaOioiKWLVvGjBkzWLRoEampqcyfP5+qqqoW2x88eJDnnnuOCRMmsGjRIrKzs3n22Wc5evSo3Wbt2rVs2LCB3NxcFixYQI8ePZg/fz4XL1602zz//POUl5czb9485syZw/79+8nLy7OPezweoqKiuPPOO8nKymqxllOnTvHMM8/wuc99jp/97GfcddddvPDCC+zevdvfH4OIiIiEIb+D0fr165k4cSLjx48nJSWF3NxcoqKi2LJlS4vtCwsLGTFiBHfffTcpKSnMnDmT9PR0Nm7cCHhHiwoLC5k2bRrZ2dmkpqby+OOPc/bsWbZv3w5ARUUFu3fv5tFHH2XIkCFkZmaSk5NDUVERZ86cASA6Oprc3Fxuv/12EhISWqzl9ddfJzExkW984xukpKQwadIkbrnlFl577TV/fwwiIiIShvwKRm63m7KyMp8RGYfDQVZWFiUlJS2+pqSk5JIRnOHDh/PBBx8A3lEcl8vFsGHD7OMxMTEMHjzYPmdJSQm9evXi2muvtdtkZWVhWRalpaWtrv+DDz5osZbL1Q7Q0NBAXV2dz38NDQ2t/poiIiISOvyaY1RdXY3H47lkRCYhIYHjx4+3+BqXy0V8fLzPc/Hx8bhcLvt483Of1SYuLs7neEREBLGxsXab1rhcLefPn+fixYtERUVd8po1a9awatUq+/HYsWP5/ve/3+qvKSIiIqGjTZOvu5N7772XKVOm2I8tywLg7NmzuN3uDvs6lmXRr18/Tp8+jTGmw84r/lE/BAf1Q3Do1H6oq6HxF/MAiJj3H2BpkfTlhMP7wVOQjzm4D8ddX8ca+fmA1OB0OunTp8+V2/lz0ri4OBwOxyWjNC6X67LzehISEi6ZmF1VVWW3b/5YVVXlU3BVVRVpaWl2m+rqap9zNDY2UlNTc9mv608tPXv2bHG0CCAyMpLIyMhLnne73R16Sa05cDU0NITsL344UD8EB/VDcOjMfjCus3gOHYQe0XjcjUBjh54/nITD+6HxeAUcOoh19jSOIJ+O4ldEdzqdpKenU1xcbD/n8XgoLi4mIyOjxddkZGSwb98+n+f27t3LkCFDAEhMTCQhIcGnTV1dHaWlpfY5MzIyqK2tpayszG5TXFyMMeYztwn434YMGdJiLZerXUREOsn5Ou/H6J6BrUPkf/F77HLKlCls3ryZrVu3UlFRwZIlS6ivr7f3Alq8eDHLly+320+ePJk9e/awbt06jh07RkFBAYcOHWLSpEmANwlPnjyZ1atXs2PHDo4ePcrixYvp06cP2dnZAKSkpDBixAjy8vIoLS3lwIEDLF26lDFjxtC3b1/7a1VUVHD48GFqamo4f/48hw8f5vDhw/bxO+64g1OnTvGHP/yBY8eO8ec//5lt27Zx1113teVnJyIibXXhvPdjdExg6xD5X/yeYzRmzBiqq6spKCjA5XKRlpbG3Llz7Utap0+ftof9AIYOHcrs2bNZsWIFL730EsnJyTzxxBMMGjTIbjN16lTq6+vJy8ujrq6OzMxM5s6d63N5a/bs2eTn5/PUU0/ZGzzm5OT41LZw4UIqKyvtx80bPBYUFADe0ak5c+bwu9/9jsLCQq666ioeffRRRowY4e+PQURE2uOCRowkOFkmVC9YBlhlZWWHzzFKTk7mxIkTIXsNORyoH4KD+iE4dGY/eN79K+a3P4ehWUT8cH6HnjvchMP7ofHXC2D3O1gPPYbji5MCUkNkZCT9+/e/YjstAxARka6nESMJUgpGIiLS9ZrmGFk9NcdIgouCkYiIdL3zzZOvNWIkwUXBSEREup5WpUmQUjASEZGupzlGEqQUjEREpOtpxEiClIKRiIh0OaMRIwlSCkYiItL17FVpCkYSXBSMRESk6+leaRKkFIxERKTraY6RBCkFIxER6XoXtI+RBCcFIxER6VLGGI0YSdBSMBIRka7lboBGt/dzjRhJkFEwEhGRrtU8WgQQHR24OkRaoGAkIiJdq3lFWo9oLEdEYGsR+V8UjEREpGtp4rUEMQUjERHpWvau15p4LcFHwUhERLqWRowkiCkYiYhIlzLa9VqCmIKRiIh0rXqNGEnwUjASEZGudb75BrKaYyTBR8FIRES6luYYSRBTMBIRka6lVWkSxBSMRESka2nESIKYgpGIiHSt8xoxkuClYCQiIl3KaFWaBDEFIxER6VpNI0ZWTwUjCT4KRiIi0rU0x0iCmIKRiIh0LTsYaY6RBB8FIxER6VoaMZIgpmAkIiJdxhijfYwkqCkYiYhI13E3QGOj93ONGEkQUjASEZGu07yHEUB0dODqELkMBSMREek6zfOLekRjOSICW4tICxSMRESk62h+kQQ5BSMREek6WpEmQU7BSEREus55BSMJbgpGIiLSZYx9KU3BSIKTgpGIiHQd3UBWgpyCkYiIdJ2mS2lWT02+luCkYCQiIl1Hk68lyCkYiYhI19FyfQlyCkYiItJ1NGIkQU7BSEREus55jRhJcFMwEhGRLmO0Kk2CnIKRiIh0naYRI6ungpEEJwUjERHpOppjJEFOwUhERLqOHYw0x0iCk4KRiIh0HY0YSZBTMBIRkS5hjNE+RhL0FIxERKRruBugsdH7uUaMJEgpGImISNdo3sMIIDo6cHWIfAYFIxER6RrN84t6RGM5IgJbi8hlKBiJiEjX0PwiCQEKRiIi0jW0Ik1CgIKRiIh0jfMKRhL8FIxERKRLGPtSmoKRBC8FIxER6RrNl9J6ao6RBC8FIxER6RpNwcjSiJEEMQUjERHpGlqVJiFAwUhERLqGVqVJCFAwEhGRrnFek68l+CkYiYhI17BHjHQpTYKXsy0v2rhxI+vWrcPlcpGamkpOTg6DBw++bPtt27axcuVKKisrSUpK4oEHHmDkyJH2cWMMBQUFbN68mdraWjIzM5k1axbJycl2m5qaGpYuXcrOnTuxLIvRo0fz8MMPE/2p++0cOXKE/Px8Dh06RFxcHJMmTWLq1Kk+tbz22mu8/vrrnD59mri4OEaPHs39999PVFRUW34UIiLSSlquL6HA7xGjoqIili1bxowZM1i0aBGpqanMnz+fqqqqFtsfPHiQ5557jgkTJrBo0SKys7N59tlnOXr0qN1m7dq1bNiwgdzcXBYsWECPHj2YP38+Fy9etNs8//zzlJeXM2/ePObMmcP+/fvJy8uzj9fV1fH000/Tr18/nnnmGR588EFefvllNm3aZLd56623WL58OV/96lf5z//8Tx599FG2bdvGSy+95O+PQURE/NW8Kq2ngpEEL7+D0fr165k4cSLjx48nJSWF3NxcoqKi2LJlS4vtCwsLGTFiBHfffTcpKSnMnDmT9PR0Nm7cCHhHiwoLC5k2bRrZ2dmkpqby+OOPc/bsWbZv3w5ARUUFu3fv5tFHH2XIkCFkZmaSk5NDUVERZ86cAbyhx+1289hjjzFw4EDGjh3LnXfeyfr16+1aDh48yNChQ7n11ltJTExk+PDhjB07ltLSUr9/cCIi4idNvpYQ4NelNLfbTVlZGffcc4/9nMPhICsri5KSkhZfU1JSwpQpU3yeGz58uB16Tp06hcvlYtiwYfbxmJgYBg8eTElJCWPHjqWkpIRevXpx7bXX2m2ysrKwLIvS0lJGjRpFSUkJ1113HU6n0+frrF27lpqaGmJjYxk6dCh/+9vfKC0tZfDgwXz00Uf84x//4Atf+MJlv+eGhgYaGhrsx5Zl0bPprx3LslrxU2ud5nN15DnFf+qH4KB+CA4d3g/2iFEv9a0fwuH9YFlgAIvg/z78CkbV1dV4PB4SEhJ8nk9ISOD48eMtvsblchEfH+/zXHx8PC6Xyz7e/NxntYmLi/M5HhERQWxsrE+bxMTES+pqPhYbG8utt95KdXU1//qv/wpAY2MjX/rSl5g2bdplv+c1a9awatUq+/E111zDokWL6N+//2Vf0x5JSUmdcl7xj/ohOKgfgkNH9cOx+gt4gP4DBxH5qTmk0jqh/H44HR3Nebz/tscGed+3afJ1qHrvvfdYs2YNs2bNYsiQIZw8eZL/+Z//YdWqVcyYMaPF19x7770+I17NSbeyshK3291htVmWRVJSEidPnsQY02HnFf+oH4KD+iE4dGQ/GGPwnK8FoPJcLdaJEx1RYrcQDu+HxgsXAKiqquJcgPre6XS2alDDr2AUFxeHw+GwR2mauVyuS0aRmiUkJFwyMbuqqspu3/yxqqqKPn36+LRJS0uz21RXV/uco7GxkZqaGp/ztFTXp7/GypUr+eIXv8jEiRMBGDRoEBcuXOC///u/mTZtGg7HpVOuIiMjiYyMbPF764xfUGNMyP7ihxP1Q3BQPwSHjugHc7EeGhu9n/eIBvWr30L5/dBctjGd829nR/Jr8rXT6SQ9PZ3i4mL7OY/HQ3FxMRkZGS2+JiMjg3379vk8t3fvXoYMGQJAYmIiCQkJPm3q6uooLS21z5mRkUFtbS1lZWV2m+LiYowx9jYBGRkZ7N+/32cUZ+/evQwYMIDY2FgA6uvrL7m22VIYEhGRDtY88RrgU9usiAQbv1PBlClT2Lx5M1u3bqWiooIlS5ZQX1/PuHHjAFi8eDHLly+320+ePJk9e/awbt06jh07RkFBAYcOHWLSpEmAd4hw8uTJrF69mh07dnD06FEWL15Mnz59yM7OBiAlJYURI0aQl5dHaWkpBw4cYOnSpYwZM4a+ffsCcOutt+J0OnnhhRcoLy+nqKiIDRs2+FwGu+mmm3jjjTd4++23OXXqFHv37mXlypXcdNNNCkgiIp2ptsb7sWcvLEdEYGsR+Qx+zzEaM2YM1dXVFBQU4HK5SEtLY+7cufblqtOnT/uMygwdOpTZs2ezYsUKXnrpJZKTk3niiScYNGiQ3Wbq1KnU19eTl5dHXV0dmZmZzJ0712fTxdmzZ5Ofn89TTz1lb/CYk5NjH4+JiWHevHnk5+czZ84cevfuzfTp07n99tvtNtOnT8eyLFasWMGZM2eIi4vjpptu4r777vP3xyAiIv441zSlonfcZ7cTCTDLBPvFviBVWVnps4y/vSzLIjk5mRMnTgT99ddwpn4IDuqH4NCR/WB2bcPzm4WQPpSIf3m2gyrsHsLh/dD46wWw+x2sBx/DcdukgNQQGRnZqsnXun4kIiKdztQ0LaCJ1YiRBDcFIxER6XxNl9IsXUqTIKdgJCIina/mnPejRowkyCkYiYhI59OlNAkRCkYiItLpNMdIQoWCkYiIdD57jlH8FRqKBJaCkYiIdD6NGEmIUDASEZHOp8nXEiIUjEREpFOZhotQ33SvNC3XlyCnYCQiIp3rXNNlNIcDevYKbC0iV6BgJCIinetT84s+fS9NkWCkYCQiIp2rORhpRZqEAAUjERHpVKZpqb4mXndjITRQqGAkIiKdy16R1juwdYi0goKRiIh0rqZLaZZGjCQEKBiJiEjnqmm6lKY5RhICFIxERKRzndOu1xI6FIxERKRT6QayEkoUjEREpHM1zzHSrtcSAhSMRESkc2nESEKIgpGIiHQaY4yCkYQUBSMREek85+ugsdH7uYKRhAAFIxER6TzNS/V7RGNF9QhsLSKtoGAkIiKdR0v1JcQoGImISOexbweiYCShQcFIREQ6jb2HkZbqS4hQMBIRkc7TNMdI90mTUKFgJCIinUdzjCTEKBiJiEjn0R5GEmIUjEREpNNojpGEGgUjERHpPM33SdOIkYQIBSMREek8mmMkIUbBSEREOo99KS0+sHWItJKCkYiIdArjdkNdjfeBRowkRCgYiYhI56hr2vXasqBXbGBrEWklBSMREekc55qCUa9YLEdEYGsRaSUFIxER6RxNu17rMpqEEgUjERHpHNrcUUKQgpGIiHQKo6X6EoIUjEREpHM0b+6opfoSQhSMRESkc9iX0noHtg4RPygYiYhI59ClNAlBCkYiItIp7BvIxupSmoQOBSMREekcTcv1rd4aMZLQoWAkIiKdQ8v1JQQpGImISOdQMJIQpGAkIiIdztTXw8WL3ge6lCYhRMFIREQ6XvPtQJxO6NEzsLWI+EHBSEREOt6nLqNZlhXYWkT8oGAkIiId75yW6ktoUjASEZEOZ5ovpWl+kYQYBSMREel4zfdJ04o0CTEKRiIi0vHOnfN+VDCSEKNgJCIiHU97GEmIUjASEZEOpzlGEqoUjEREpOPpBrISohSMRESk451rnnzdO8CFiPhHwUhERDpe84iRLqVJiFEwEhGRDmU8HqjRqjQJTQpGIiLSsc7XgvF4P1cwkhCjYCQiIh3rXNOKtJ4xWM7IwNYi4icFIxER6Vjaw0hCmIKRiIh0LAUjCWHOtrxo48aNrFu3DpfLRWpqKjk5OQwePPiy7bdt28bKlSuprKwkKSmJBx54gJEjR9rHjTEUFBSwefNmamtryczMZNasWSQnJ9ttampqWLp0KTt37sSyLEaPHs3DDz9MdHS03ebIkSPk5+dz6NAh4uLimDRpElOnTvWppba2lpdeeol3332Xmpoa+vfvzze/+U2fekREpO3MOQUjCV1+jxgVFRWxbNkyZsyYwaJFi0hNTWX+/PlUVVW12P7gwYM899xzTJgwgUWLFpGdnc2zzz7L0aNH7TZr165lw4YN5ObmsmDBAnr06MH8+fO5ePGi3eb555+nvLycefPmMWfOHPbv309eXp59vK6ujqeffpp+/frxzDPP8OCDD/Lyyy+zadMmu43b7ebpp5+msrKSH/zgB/zyl7/kkUceoW/fvv7+GERE5HJ0A1kJYX4Ho/Xr1zNx4kTGjx9PSkoKubm5REVFsWXLlhbbFxYWMmLECO6++25SUlKYOXMm6enpbNy4EfCOFhUWFjJt2jSys7NJTU3l8ccf5+zZs2zfvh2AiooKdu/ezaOPPsqQIUPIzMwkJyeHoqIizpw5A8Bbb72F2+3mscceY+DAgYwdO5Y777yT9evX27W8+eab1NTU8MQTT5CZmUliYiLXX389aWlp/v4YRETkcrSHkYQwvy6lud1uysrKuOeee+znHA4HWVlZlJSUtPiakpISpkyZ4vPc8OHD7dBz6tQpXC4Xw4YNs4/HxMQwePBgSkpKGDt2LCUlJfTq1Ytrr73WbpOVlYVlWZSWljJq1ChKSkq47rrrcDqdPl9n7dq11NTUEBsby86dOxkyZAj5+fns2LGDuLg4xo4dyz333IPD0XJGbGhooKGhwX5sWRY9e/a0P+8ozefqyHOK/9QPwUH9EBza3A/NI0a949WHHSAc3g8WFgawCP7vw69gVF1djcfjISEhwef5hIQEjh8/3uJrXC4X8fG+98qJj4/H5XLZx5uf+6w2cXG+f3lEREQQGxvr0yYxMfGSupqPxcbG8tFHH1FZWcmtt97Kv/zLv3Dy5EmWLFlCY2MjX/3qV1usf82aNaxatcp+fM0117Bo0SL69+/fYvv2SkpK6pTzin/UD8FB/RAc/O2Hyov1XADirx5I7Kfmikr7hPL74XR0NOeB+Pi4oP+daNPk61BljCEuLo5HHnkEh8NBeno6Z86c4dVXX71sMLr33nt9Rryak25lZSVut7vDarMsi6SkJE6ePIkxpsPOK/5RPwQH9UNwaGs/uD+uBKCq0XDuxInOKq/bCIf3Q+OFCwBUVVUH7HfC6XS2alDDr2AUFxeHw+GwR2mauVyuS0aRmiUkJFwyMbuqqspu3/yxqqqKPn36+LRpnvuTkJBAdXW1zzkaGxupqanxOU9LdX36ayQkJOB0On0um1199dW4XC7cbrfPZbhmkZGRREa2vEFZZ/yCGmNC9hc/nKgfgoP6ITj43Q+fWq6v/us4ofx+MHjrDoXvwa/J106nk/T0dIqLi+3nPB4PxcXFZGRktPiajIwM9u3b5/Pc3r17GTJkCACJiYkkJCT4tKmrq6O0tNQ+Z0ZGBrW1tZSVldltiouLMcbY2wRkZGSwf/9+n1GcvXv3MmDAAGJjYwEYOnQoJ0+exOPx2G1OnDhBnz59WgxFIiLSBtrHSEKY36vSpkyZwubNm9m6dSsVFRUsWbKE+vp6xo0bB8DixYtZvny53X7y5Mns2bOHdevWcezYMQoKCjh06BCTJk0CvEOEkydPZvXq1ezYsYOjR4+yePFi+vTpQ3Z2NgApKSmMGDGCvLw8SktLOXDgAEuXLmXMmDH2Uvtbb70Vp9PJCy+8QHl5OUVFRWzYsMHnMtgdd9xBTU0NL774IsePH2fXrl2sWbOGL3/5y23+AYqIyCdMQwOcr/M+0Ko0CUF+D5OMGTOG6upqCgoKcLlcpKWlMXfuXPty1enTp31mnA8dOpTZs2ezYsUKXnrpJZKTk3niiScYNGiQ3Wbq1KnU19eTl5dHXV0dmZmZzJ07l6ioKLvN7Nmzyc/P56mnnrI3eMzJybGPx8TEMG/ePPLz85kzZw69e/dm+vTp3H777Xabfv368eSTT/K73/2OJ554gr59+3LnnXf6rLITEZF2qG0aLXI4oGevwNYi0gaWCfaLfUGqsrLSZxl/e1mWRXJyMidOnAj666/hTP0QHNQPwaEt/WAqPsTz0+9D73gi/uP3nVxh9xAO74fG/1oA/3gH64Hv4hh3Z0BqiIyMbNXka90rTUREOo5uByIhTsFIREQ6jNGu1xLiFIxERKTj2CvS4j+7nUiQUjASEZGOc043kJXQpmAkIiIdp6ZpQ18FIwlRCkYiItJxas55P2qOkYQoBSMREekwRrteS4hTMBIRkY5zznspTXOMJFQpGImISMfRcn0JcQpGIiLSIYwxWq4vIU/BSEREOkb9eXC7vZ/rUpqEKAUjERHpGM23A4mKwurRI7C1iLSRgpGIiHQMXUaTMKBgJCIiHUNL9SUMKBiJiEiHMOe067WEPgUjERHpGE0jRpaW6ksIUzASEZGOoUtpEgYUjEREpGM03ydNwUhCmIKRiIh0CM0xknCgYCQiIh3DnmOk5foSuhSMRESkY2iOkYQBBSMREekY5xSMJPQpGImISLuZxkaoq/E+0HJ9CWEKRiIi0n51NWCM9/NevQNbi0g7KBiJiEj7Nc8vionFiogIbC0i7aBgJCIi7ael+hImFIxERKT9mkeMNL9IQpyCkYiItJvRUn0JEwpGIiLSfk1L9S0FIwlxCkYiItJ+9qU07XotoU3BSERE2k+X0uSzWFagK2g1BSMREWk3zTGScKFgJCIi7dc8x0ir0iTEKRiJiEj7acRIwoSCkYiItJ+CkYQJBSMREWkXc7Ee6i94HygYSYhTMBIRkfZpHi2KcELPmMDWItJOCkYiItI+n7qMZoXQsmyRligYiYhI+9jBqHdg6xDpAApGIiLSLuacdr2W8KFgJCIi7VOj+6RJ+FAwEhGR9tFSfQkjCkYiItI+56q8H7XrtYQBBSMREWkX3SdNwomCkYiItE/NOe9HBSMJAwpGIiLSPpp8LWFEwUhERNrMuN3w8Snvg7iEgNYi0hEUjEREpO1KiuHCee8eRgMGBroakXZTMBIRkTYzu98BwBoxGssREeBqRNpPwUhERNrEeDyYf/wd8AYjkXCgYCQiIm1z5BC4PoYePeG64YGuRqRDKBiJiEib2JfRbhiJFRkV4GpEOoaCkYiItIn5hzcYceMtgS1EpAMpGImIiN/MyQo4UQ4RTqysmwNdjkiHUTASERG/NU+6ZmgWVkyvwBYj0oEUjERExG/2/KIbtRpNwouCkYiI+MW4zkDZQUDL9CX8KBiJiIhfzO6my2jXZGAlXBXYYkQ6mIKRiIj45ZPLaFqNJuFHwUhERFrN1NXCgX2AgpGEJwUjERFpNbNvBzS6ISkFKykl0OWIdDgFIxERab2m+UVajSbhSsFIRERaxTQ0YPbtBMC68fMBrkakczjb8qKNGzeybt06XC4Xqamp5OTkMHjw4Mu237ZtGytXrqSyspKkpCQeeOABRo4caR83xlBQUMDmzZupra0lMzOTWbNmkZycbLepqalh6dKl7Ny5E8uyGD16NA8//DDR0dF2myNHjpCfn8+hQ4eIi4tj0qRJTJ06tcWa3n77bZ577jluvvlmfvSjH7XlxyAi0r0c2AP15yGhL6Re/v/5IqHM7xGjoqIili1bxowZM1i0aBGpqanMnz+fqqqqFtsfPHiQ5557jgkTJrBo0SKys7N59tlnOXr0qN1m7dq1bNiwgdzcXBYsWECPHj2YP38+Fy9etNs8//zzlJeXM2/ePObMmcP+/fvJy8uzj9fV1fH000/Tr18/nnnmGR588EFefvllNm3adElNp06d4ve//z3XXXedv9++iEi31XxvNGvEaCyHLjhIePJ7xGj9+vVMnDiR8ePHA5Cbm8uuXbvYsmUL99xzzyXtCwsLGTFiBHfffTcAM2fOZN++fWzcuJHvfOc7GGMoLCxk2rRpZGdnA/D444+Tm5vL9u3bGTt2LBUVFezevZuFCxdy7bXXApCTk8PChQt56KGH6Nu3L2+99RZut5vHHnsMp9PJwIEDOXz4MOvXr+f222+36/F4PPzqV7/ia1/7Gvv376e2tvYzv9+GhgYaGhrsx5Zl0bNnT/vzjtJ8ro48p/hP/RAc1A/B4dP9YDyNmN3veh/f+Hn1TRcKh/eDBZimj8H+ffgVjNxuN2VlZT4ByOFwkJWVRUlJSYuvKSkpYcqUKT7PDR8+nO3btwPe0RuXy8WwYcPs4zExMQwePJiSkhLGjh1LSUkJvXr1skMRQFZWFpZlUVpayqhRoygpKeG6667D6XT6fJ21a9dSU1NDbGwsAKtWrSIuLo4JEyawf//+K37Pa9asYdWqVfbja665hkWLFtG/f/8rvrYtkpKSOuW84h/1Q3BQPwSHpKQk6t/fzalzLqxesQy47UtYkZGBLqvbCeX3w+noaM4D8fHxxH5qmkww8isYVVdX4/F4SEhI8Hk+ISGB48ePt/gal8tFfHy8z3Px8fG4XC77ePNzn9UmLi7O53hERASxsbE+bRITEy+pq/lYbGwsBw4c4M033+RnP/tZK75br3vvvdcn2DUn3crKStxud6vPcyWWZZGUlMTJkycxxnTYecU/6ofgoH4IDp/uB/em17xP3nAzJ0+fDmxh3Uw4vB8aL1wAoKqqinMnTgSkBqfT2apBjTZNvg5F58+f51e/+hWPPPLIJSHrs0RGRhJ5mb+MOuMX1BgTsr/44UT9EBzUD8HB4/Fgdm0DvMv01SeBEcrvB2N/DP7vwa9gFBcXh8PhsEdpmrlcrktGkZolJCRcMjG7qqrKbt/8saqqij59+vi0SUtLs9tUV1f7nKOxsZGamhqf87RUV/Oxjz76iMrKShYtWmQfb+6cmTNn8stf/jKkhylFRDrN8aNQeRKckfC5kVduLxLC/ApGTqeT9PR0iouLGTVqFOD9S6K4uJhJkya1+JqMjAz27dvHXXfdZT+3d+9ehgwZAkBiYiIJCQns27fPDkJ1dXWUlpZyxx132Oeora2lrKyM9PR0AIqLizHG2NsEZGRk8NJLL+F2u+15Rnv37mXAgAHExsYSFRXFz3/+c5/aVqxYwYULF/jWt75Fv379/PlRiIh0G82r0bh+BFZ0z8AWI9LJ/F5vOWXKFDZv3szWrVupqKhgyZIl1NfXM27cOAAWL17M8uXL7faTJ09mz549rFu3jmPHjlFQUMChQ4fsIGVZFpMnT2b16tXs2LGDo0ePsnjxYvr06WOvUktJSWHEiBHk5eVRWlrKgQMHWLp0KWPGjKFv374A3HrrrTidTl544QXKy8spKipiw4YN9vygqKgoBg0a5PNfr169iI6OZtCgQT6TtkVE5BOe5stoI7TbtYQ/v9PAmDFjqK6upqCgAJfLRVpaGnPnzrUvaZ0+fdpnKd7QoUOZPXs2K1as4KWXXiI5OZknnniCQYMG2W2mTp1KfX09eXl51NXVkZmZydy5c4mKirLbzJ49m/z8fJ566il7g8ecnBz7eExMDPPmzSM/P585c+bQu3dvpk+f7rNUX0RE/OM+dRKOHgLLgTV8VKDLEel0lgn2WVBBqrKy0md/o/ayLIvk5GROnDgR9BPTwpn6ITioH4KDZVnEbv8Lrryfw5DrifjRM4EuqVsKh/dD428Wwq5tWA88imPc5IDUEBkZ2apVadq6VERELuv8tq0AWCNuCWwhIl1EwUhERFpkaqqpL/4HANaNCkbSPSgYiYhIi8zeHeBphJQ0rP7azkS6BwUjERFpkX3TWI0WSTeiYCQiIpcwDRcx7+0CwDFcy/Sl+1AwEhGRSx3YCxfribgqEVKvvXJ7kTChYCQiIpcwe94FIHr0F3z2phMJdwpGIiLiwxiD2bMdgJ6jvxjgakS6loKRiIj4OloGro+hRzTRw24OdDUiXUrBSEREfJg9fwfAuv5GrKgeAa5GpGspGImIiI/my2i6N5p0RwpGIiJiM2dON9001sLSZTTphhSMRETEZvZ6V6ORPhQrLiGgtYgEgoKRiIjYdBlNujsFIxERAcBcOA8H9gAKRtJ9KRiJiIjX+7vB7Yb+SZA8MNDViASEgpGIiACf7HZtDR+l3a6l21IwEhERjKcRs28HANaw7ABXIxI4CkYiIgJlJXCuCnr2giGfC3Q1IgGjYCQiIp9cRsu6CcvpDHA1IoGjYCQiInYwQpfRpJtTMBIR6ebMqRNwohwiIrBuuCnQ5YgElIKRiEg3Z+92Pfh6rF6xgS1GJMAUjEREujmzu2l+0Qht6iiiYCQi0o2Z2hr44D0ArGEKRiIKRiIi3Zgp3gkeDyQPxEpMDnQ5IgGnYCQi0p3t0WU0kU9TMBIR6aaM240p3gXoMppIMwUjEZHu6oP34Hwt9I6H9IxAVyMSFBSMRES6KbN3OwBW1s1YjogAVyMSHBSMRES6IWPMJ7cBGa7LaCLNFIxERLqjE+VQeRKcTrh+RKCrEQkaCkYiIt2QfW+0zOFY0T0DW4xIEFEwEhHphnQZTaRlCkYiIt2MqXZB2UEArGHZgS1GJMgoGImIdDNm304wBgZdi9W3X6DLEQkqCkYiIt3Nof0AWJ8bEdg6RIKQgpGISDdjyj8EwBp0bYArEQk+CkYiIt2IaWyEY0e8D1KuCWwxIkFIwUhEpDs5dRwaLkKPaEhMCnQ1IkFHwUhEpBsxR8u8n1ydqtuAiLRAwUhEpDupOAyANVCX0URaomAkItKNmPKmEaOB6YEtRCRIKRiJiHQnzSNGKWkBLUMkWCkYiYh0E6b6LFSdBcsCBSORFikYiYh0F+WHvR8TB2D1iA5oKSLBSsFIRKSbaJ5fpInXIpenYCQi0l00jxjpMprIZSkYiYh0E/aI0SCtSJOuZgW6gFZTMBIR6QbMxXr46Jj3gW4FInJZCkYiIt3B8aPg8UBsHCT0DXQ1IkFLwUhEpBsw5R96Pxl4DZYVOpc1RLqagpGISHfQFIy0Ik3ksykYiYh0A58eMRKRy1MwEhEJc8bjgYqmESNNvBb5TApGIiLh7uNTcOE8OJ2QlBLoakSCmoKRiEi4a76MNmAQltMZ2FpEgpyCkYhImDOaeC3SagpGIiJhzjTNL9LGjiJXpmAkIhLu7BEj3QpE5EoUjEREwpipq/FOvgbdPFakFRSMRETCWflh78erErF6xQa0FJFQoGAkIhLGPplflBbQOkRCRZvWbW7cuJF169bhcrlITU0lJyeHwYMHX7b9tm3bWLlyJZWVlSQlJfHAAw8wcuRI+7gxhoKCAjZv3kxtbS2ZmZnMmjWL5ORku01NTQ1Lly5l586dWJbF6NGjefjhh4mOjrbbHDlyhPz8fA4dOkRcXByTJk1i6tSp9vFNmzbx17/+lfLycgDS09O57777PrN2EZGQVl4GaH6RSGv5PWJUVFTEsmXLmDFjBosWLSI1NZX58+dTVVXVYvuDBw/y3HPPMWHCBBYtWkR2djbPPvssR48etdusXbuWDRs2kJuby4IFC+jRowfz58/n4sWLdpvnn3+e8vJy5s2bx5w5c9i/fz95eXn28bq6Op5++mn69evHM888w4MPPsjLL7/Mpk2b7Dbvv/8+Y8eO5d/+7d94+umnueqqq3j66ac5c+aMvz8GEZGQYJoupWmpvkjr+D1itH79eiZOnMj48eMByM3NZdeuXWzZsoV77rnnkvaFhYWMGDGCu+++G4CZM2eyb98+Nm7cyHe+8x2MMRQWFjJt2jSys7MBePzxx8nNzWX79u2MHTuWiooKdu/ezcKFC7n22msByMnJYeHChTz00EP07duXt956C7fbzWOPPYbT6WTgwIEcPnyY9evXc/vttwMwe/Zsn9oeffRR/v73v7Nv3z5uu+22Fr/fhoYGGhoa7MeWZdGzZ0/7847SfC7d9Tqw1A/BQf3QMYzbDcePAGANSvf756l+CA7h0A+WBQawCP7vw69g5Ha7KSsr8wlADoeDrKwsSkpKWnxNSUkJU6ZM8Xlu+PDhbN++HYBTp07hcrkYNmyYfTwmJobBgwdTUlLC2LFjKSkpoVevXnYoAsjKysKyLEpLSxk1ahQlJSVcd911OD+1q+vw4cNZu3YtNTU1xMZeOumwvr4et9vd4rFma9asYdWqVfbja665hkWLFtG/f//LvqY9kpKSOuW84h/1Q3BQP7TPxcOlfOR2Y8X0IvmG4ViOtk0rVT8Eh1Duh9PR0ZwH4uLi6f2paTLByK9gVF1djcfjISEhwef5hIQEjh8/3uJrXC4X8fHxPs/Fx8fjcrns483PfVabuLg4n+MRERHExsb6tElMTLykruZjLYWfP/7xj/Tt25esrKwWawe49957fYJdc9KtrKzE7XZf9nX+siyLpKQkTp48iTGmw84r/lE/BAf1Q8fw7HoXAHN1Kic/+sjv16sfgkM49EPjhQsAVFdXUXPiREBqcDqdrRrU6LY3zXnllVd4++23+clPfkJUVNRl20VGRhIZGdnisc74BTXGhOwvfjhRPwQH9UP7mOaJ1ynXtOvnqH4IDqHcD81lG9M5/3Z2JL/GVePi4nA4HPYoTTOXy3XJKFKzhISESyZmV1VV2e2bP16pTXV1tc/xxsZGampqfNq0VNenv0azV199lVdeeYV58+aRmpraYt0iIqGu+R5paOK1SKv5FYycTifp6ekUFxfbz3k8HoqLi8nIyGjxNRkZGezbt8/nub179zJkyBAAEhMTSUhI8GlTV1dHaWmpfc6MjAxqa2spKyuz2xQXF2OMsZfaZ2RksH//fp/LW3v37mXAgAE+l9HWrl3Ln/70J+bOneszZ0lEJJwYYz51KxAFI5HW8nsm3pQpU9i8eTNbt26loqKCJUuWUF9fz7hx4wBYvHgxy5cvt9tPnjyZPXv2sG7dOo4dO0ZBQQGHDh1i0qRJgPfa6eTJk1m9ejU7duzg6NGjLF68mD59+tir1FJSUhgxYgR5eXmUlpZy4MABli5dypgxY+jbty8At956K06nkxdeeIHy8nKKiorYsGGDz/ygV155hZUrV/Ld736XxMREXC4XLpeLC03XPkVEwkbVGaipBssBAwYFuhqRkOH3HKMxY8ZQXV1NQUEBLpeLtLQ05s6da1+uOn36tM9SvKFDhzJ79mxWrFjBSy+9RHJyMk888QSDBn3yRp06dSr19fXk5eVRV1dHZmYmc+fO9Zn7M3v2bPLz83nqqafsDR5zcnLs4zExMcybN4/8/HzmzJlD7969mT59ur1UH+CNN97A7XbzH//xHz7f04wZM/ja177m749CRCR4NV9GS7oaK6pHYGsRCSGWCfZZUEGqsrLSZ3+j9rIsi+TkZE6cOBH0E9PCmfohOKgf2s9T+DJmze+xRt2GI/f/a9M51A/BIRz6ofE3z8CuIqz7H8UxfnJAaoiMjGzVqjTdK01EJBzZE6/TAlqGSKhRMBIRCUPNN4/VPdJE/KNgJCISZkz9BfioadNdjRiJ+EXBSEQk3Bw74t1JL74PVlyfQFcjElIUjEREwoy9sWNKWkDrEAlFCkYiIuFG84tE2kzBSEQkzGjESKTtFIxERMKI8Xig4jAA1iCNGIn4S8FIRCScVJ6E+gsQGQWJAwJdjUjIUTASEQkn5U032746FSsiIrC1iIQgBSMRkTBiyg8DYA28JrCFiIQoBSMRkTBimkeMFIxE2kTBSEQknDRPvFYwEmkTBSMRkTBhaqrh7GnvAy3VF2kTBSMRkTBhdhV5P+mfhBUdE9hiREKUgpGISBgw+3ZglucBYN0yLrDFiIQwBSMRkRBnPngfzwvPQGMj1ujbsKbMDHRJIiFLwUhEJISZig/x/Orf4eJFyLoZ61vfx3Lof+0ibaV3j4hIiDKnTuD55U/gfC0Mvg7HIz/GcjoDXZZISFMwEhEJQcZ1Bs8v/w2qzkJKGo7/869YPXoEuiyRkKdgJCISYkxtjTcUVZ6E/kk4/vmnWDGxgS5LJCwoGImIhBBTX49n8b/DsSMQ3wfH/30KK75PoMsSCRsKRiIiIcK43d7VZ6X7IaaXd6Sof1KgyxIJKwpGIp9ijpZRvep3GLc70KWI+DAeD+Z/noPinRAVheP//D8s7W4t0uG0fEGkiXE34Pn1fKo+PoWj7jzW7XcHuiQRmynIx7z7F4iIwPHdf8EafF2gSxIJSxoxkpBnGi5ijGn/eYrehI9PAeDZ8hrG42n3OUU6gtnzLmbzOrAsrIf/GeuGmwJdkkjYUjCSkGZKivHMnolZ/kL7zuNuwBS+/MkTp07Ae/9oZ3Ui7WfqL+B56b8BsO64B8fo2wJckUh4UzCSLmeM6ZDRGOPx4Fm5BNxuzNYNmLKDbT/Xti3e0aL4PvSaPB0Az5vr212jSHuZ1wq8v5t9+2N95b5AlyMS9hSMpEsZjwfPz/4Fz9zvYKpd7TvX9r/B0TL7sWflkjZdUjPuBu8/PoDjy9Pofe+DYFlQvBPz0fF21SjSHuZEOeb1VwBw3JeL1SM6sAWJdAMKRtK19rwLpe/Dx6cwf/pdm09j3A2YV/4AgDV+MvSIhrKDmHf/6v+5mkeL4hKwbptE5ICB9hwOs7WwzTWKtIcxBs8fX4BGNwwfhTXilkCXJNItKBhJlzHG4Nn4p08eF23GlLzXtnP9ZSOc/gji+2BN/xbWnTO8z6/+Haa+vvXncbvt0SJr0nT7L3Jr4hTv8bc3YS6cb1ONIu1h/r4VDu7zLs2fmRvockS6DQUj6TofvA9lB8EZCSM/D4Bn+Qt+7xlkztdh1q8EwPrKfVg9orG+NBX69oczpzFvrGn9uba9+clo0Rcn2c9b198IiQPgfB3mna1+1SfSXqa2BlOwFADrrq9j9funAFck0n0oGEmXaR4tssZMxPHQ9yC2Nxw7gnlznV/nMX9eDTXV8E9XY936Je85o3pgzfiW9/iGP2HOfnzl83x6tOjL03xuwGk5HN5LdIDZ8lqHbAcg0lrmld/DuSpIHoh1xz2BLkekW1Ewki5hjh2BfTu8+7DccQ9WbBzW9G95j736EubM6dadx3UG88ZaABzTHsKKiLCPWTffCoOvg4v1mDXLrnyud5rmFvWOx7rtzkuOW2MmeucuHT/qvaQh0gXMhyXeS8WA44FHsZyRAa5IpHtRMJIuYf7cdHlr5Oex/mkA0BQ8rs2E+gt4Cpa07jzrV8DFekgfCjd+3ueYZVk4vj7L227bFsyHH1z+PD5zi3xHi+zzxfTC+vx4QEv3pWsYTyOeP/wGjMG6ZTzW0KxAlyTS7SgYSaczZyq9tzIAHF+ebj9vORw4HvwuOBywswhTvPOzz3OyAvO3173nmf4tLMu6pI2VNgTr8xMA8Kz87WUvgZl3tngnb19mtMg+3/i7vJ/sfhfTtCu2SGcxWzfA0UMQ0wvrq98KdDki3ZKCkXQ6s+lVaGyEoVlY1wzxOWalXIM14SsAeJbnYRouXvY8njV/AI/Hu3Q543OXbWdNewiiesChA5gdb11azyWjRZffG8YaMAgyh4HxYP6y4TO/T5H2MK4zn2xBce9DWHF9AlyRSPekYCSdytTWYP7aNMozaVqLbay774OEvlB5ErPhTy22MYcOwK4isBw47n3oM7+mlXDVJ8v3V72Iuei7fN/8fWurRouaOSY0Ld3/2+uXnEuko5iXl8L5OkgbgvXFLwe6HJFuS8FIOpXZWgj15yElDT43ssU2Vs8YrK81zQ3asApzyne3aWMMntXezSCtMeOxrk694te17rgH+vaDM5X2zsHQ0kq0VuwkPDzbuxVAzTnM9ktHoETay+zf492c1PJeXrYcEVd+kYh0CgUj6TTmYr33juA0hZAW5gQ1s24eC9ePAHeD95Lap+cG7dsBJe9BZBTW3fe36mtbUT0+WfW2YRXG5V2+b/6+FSpPekeLxl15tAjAckRgjWtauv/mei3dlw5lGhq8O1wD1rg7sVIHB7YgkW5OwUg6jdm2xbsXS9/+3qX0n8GyLBz3PQJOp/eu9ru2ec/hacSz2rv03powBatv/1Z/fSv7C95VbxfrMat/37bRouZzfeFLEBnlnRh76ECrXydyJebPq+GjY95d3O95MNDliHR7CkbSKYynEfO6d4m+dcc9WE7nFV9jJV2NNanpzvYrfou5UIfZthWOHYGYWHveUGv5Lt9/E1OQ7/dokX2u2DisUV/wnmvLa369VuRyzMmKT8L6V3OwYnoFuCIRUTCSzvGPd+DUCejV296dujWsO2dAv38C18feUZ5X/+h9fvJXsXrF+l2GdU0G1i3evYiaA4315XvbdJdyq3kS9s63Ma4zfr9e5NOMx4Pnd4vB3QCfuxFr1BcDXZKIoGAkncAYg6dpdZk1/i7/LllF9cBx/yPe82x5Dc6chr79sCbc1eZ6rGnf8C7fB4iNo3m+kN/nGXStd2ftxkbMX//c5npEoOlGyKXvQ4+eOB763mfOwRORrqNgJB3v4D44UgpRUW0KNFbWzXDjLZ88nvoAVmRUm8ux+lyFNfWBpnPd36bRIvtcTRs+mr9uxLgb2nwe6d7Mx6cwf2paaTn9G1hXJQa4IhFpduWJHyJ+sm8WO/Z2rN7xbTqH4+u5eA4dgMRkrFvGtbsmxx33YMZMwIqNa9d5rJGfx8T3haozmJ1FWKNva3dt0r0YY/D84b+821gMvq5Ve2mJSNfRiJF0KFP+oXdVmeXA+tI9bT6PdVV/HIvycTyxoMP2dGlvKAKwnJH25ntm8zpMg0aNxD/mna1QvAuckTi++X+wHPrfsEgw0TtSOpTZuBrw7ktk9U9q17ksZ2RQbnRn3TYJIpzwYQmeubl4Xn8Fc+F8oMuSEGCqz2JWem+YbN19H1ZSSoArEpH/TcFIOow5/RFmx98A7z3IwpUV3wfr2//XexsT1xnMy0vx/PjbeNYux5yrDnR5EsTM8v+G2nMwKL1dI6oi0nk0xyhMmSNNd+hu56hNq76WMXDoAJ5Xmm7yev0I7wquMObI/gJmxC2Yd7Z4R8lOHcesX4F5fQ3WF+7w7t3kx2aUEv7Mrm2YnW+Dw+G9hNaKvb1EpOvpnRmGzP49eP7z/0GEE+uhx3CMmdg5X6fhImb73zCb13t3hAaIiMBx19c75esFGysyEusLd2DGToR/vIOncBUcPeSde7R1A9Ytt2F9eTpWsi6XdHemtgbP8qbbfkyaHvZ/OIiEMgWjMGMaLuL5w2/AGHA3YP7nOTyHS7G+9u0O+wvVnDmN+ctGzF83Qk3TpaPIKKxRX8S6/StYKdd0yNcJFZYjAm4ai2PkGHh/N54Nq+DgPszbmzFFb8KNt+C4cwZW2pBAlyoBYl5eClVnIelqrCnd4w8HkVClYBRmzGsFcOo4xPfFGjsRU/gyZstrmIoPcTz6Y6y4Pm07rzHwwfvem6j+Y5v3khl4N18cdxfWF77UIau+QpllWfC5G4n43I2YQwe82xbs/jvs2oZn1za4fgSOO2fA0Cxt5teNmPd3Y97eBJblvYTWjj25RKTzKRiFEXP8qL0qzHHfd7BuGoO5JgPP0v+ED97H8+8/wPHdOVjpQ1t/zov1mHf/inlzPZR/+MmBjBtwTJwCw0djRQTfyrFAs67NJOJ7T2KOHcVs/BPm3b94R5Pe3w3pQ70BaVi2lmqHOXPhPJ5li4GmXeAHXx/gikTkShSMwoTxePD8/r+g0Q3DR8HIzwNgjRiNY+7P8fx6AZyswPPsv2Dd/yiOL9zx2eerPIn5ywbMW5u8q2jAu5P16HFYE+7qdpfL2sq6ehDWt/8vZur9mD+vwbz1BpQdxPPr+XB1qne+SfYXFC7DlHnlD/DxKejbH+vehwJdjoi0goJRmDBvvdF036VoHPc94nOpxkpK8Yajpb+E3e9gli3Gc6QUa2YuljPyk3N4PLB/D54tr8He7d55SgBXJWKNuxPrVl0uayur3z9hPfAoZsrXMZtexWwthGNHMPn/gVn7R6wvT8NKHQwY72VKY5r+84Ch6aPx7p+UnuHTb9K1jOtjzM5t4L742Q0vXPCOtAKObzyOFd2zC6oTkfZSMAoDpuos5k8vAk33Fbvq0mXiVs8YHN+dg9mwCrP2j97J0xWHcTz6Y4iKxmx7E7OlED469smLrr8Rx4S7IOumoNxoMRRZ8X2wpn8Tc+d0zJZCzKZX4fRHmD/+BtPak/Tt5x1puvVLmq/ShYzrY8zG1d6bv/pxnzxrzESsz93YiZWJSEdSMAoDZuUSqKuFQddiTZhy2XaWw4F119cwg9Lx/PYX3r2Hfvp9aLgI9Re8jXrGeP9HPu5O7crbiayYWG9f3D4V89brmL+9DufrwLI+9Z8DHBbQ9NjhgKozcOY0Znke5rUCrDvuxbptUrtujCufrcVAlD4U65+uvvKLe8VifeW+zi1QRDqUglGIM8U7Mdv/BpbDO1zfirkqVtbNOOb9As9/LYRjR7xPJg/0zh26ZRxWdEwnVy3NrB49sCZ+BSZ+pVXtzcV6zNubMBv/5A1ILy/FbHgZ6/ap3sm9Mb06ueLuw5z92Dtx/q9//iQQDb4ex933QeYwrSwUCVMKRiHM1F/w7lkEWBO/gpXa+k3jrMQBOOb8DPP2JqyrU7WEPERYUT2wxt+F+cIdmG1bMBtWQeVJzCt/wPx5DdbEKd7fBc0FazMFIpHuTcEohJl1K5pWvPTDmnq/36+3ont6Rysk5FjOpl23x0z07j5e+DKcKMesX4l5Yy3WLeOgV+8rn6hPP6zPj9elOJo2Lv3zagUikW5OwShEmfIPMW+8AoDj/ke14qWbsiIisG4Zhxn1Re9tSV5bCeUfeufDtJJZ+0esLzVdiuvZ/S6jmtMfYTb8CVO0Cdxu75MKRCLdloJRCDKeRjy//7V3WffIMVjDRwW6JAkwy+GAm8bgGPl52LsDc2DPJ9stXI7Hgyne6b0Ut+b3mD+vxprwFe9tXVoz2hTizKnjmMJVmHe2QGOj98mMG3BM+boCkUg31qZgtHHjRtatW4fL5SI1NZWcnBwGDx582fbbtm1j5cqVVFZWkpSUxAMPPMDIkSPt48YYCgoK2Lx5M7W1tWRmZjJr1iySk5PtNjU1NSxdupSdO3diWRajR4/m4YcfJjr6k0sAR44cIT8/n0OHDhEXF8ekSZOYOnWqX7WEArN1A3xYAj1jcNyXG+hyJIhYlgXDs7GGZ7eqvWlsxGz/K+a1l+FkBWb9Cu+luHF3Yt0xtc23kAlm5kQFprAA8/e/eveHArhuOI4pX8fKuCGwxYlIwPl9P4KioiKWLVvGjBkzWLRoEampqcyfP5+qqqoW2x88eJDnnnuOCRMmsGjRIrKzs3n22Wc5evSo3Wbt2rVs2LCB3NxcFixYQI8ePZg/fz4XL36ygdrzzz9PeXk58+bNY86cOezfv5+8vDz7eF1dHU8//TT9+vXjmWee4cEHH+Tll19m06ZNftUS7MzZjzFrfg+Ade83sBKuCnBFEsqsiAgct4zH8dPF3j2tUq6B+vOYP6/GMycXz4rfYs6cDnSZHcJUHMaT9zM8//Y9zDtbvaEo62Ycc35GxA/+XaFIRIA2jBitX7+eiRMnMn78eAByc3PZtWsXW7Zs4Z577rmkfWFhISNGjODuu+8GYObMmezbt4+NGzfyne98B2MMhYWFTJs2jexs71+5jz/+OLm5uWzfvp2xY8dSUVHB7t27WbhwIdde6115lZOTw8KFC3nooYfo27cvb731Fm63m8ceewyn08nAgQM5fPgw69ev5/bbb29VLS1paGigoeGTzdwsy6Jnz544O+hO9c1MYQEfnz2No8GN+ayt/s5VEXH1ILg6jYgv3e3d30Y6TPPlk8jISO+Nc7uTW8bBLeMwH7yH529vwLHDcLgElvwCBgyELvxVs7D4ONJ55fdDazW64Xi59/P0DKyhWd7b4iQPbP+5w1i3fj8EkXDoh4gBKZiqoVCyD44fvnL7u74G/ZOv2M4frf13269/3d1uN2VlZT4ByOFwkJWVRUlJSYuvKSkpYcoU300Hhw8fzvbt2wE4deoULpeLYcOG2cdjYmIYPHgwJSUljB07lpKSEnr16mWHIoCsLO/y8tLSUkaNGkVJSQnXXXedzzc+fPhw1q5dS01NDbGxsVespSVr1qxh1apV9uOxY8fy/e9/nz59OvgSwze/17Hnk3bp169foEsInP7jYMy4QFchQaRbvx+CSEj3w3d/FOgKWs2vS2nV1dV4PB4SEhJ8nk9ISMDlcrX4GpfLRXx8vM9z8fHxdvvmj1dqExfnuy9LREQEsbGxPm1aquvTX+NKtbTk3nvv5cUXX7T/y83N9RlB6ijnz5/nxz/+MefPn+/wc0vrqR+Cg/ohOKgfgoP6oWtpVdoVREZGEhnZ+TfsNMbw4YcfhuwwabhQPwQH9UNwUD8EB/VD1/JrxCguLg6Hw3HJCEtLozXNEhISLpmYXVVVZbdv/nilNtXV1T7HGxsbqamp8WnTUl2f/hpXqkVERES6N7+CkdPpJD09neLiYvs5j8dDcXExGRkZLb4mIyODffv2+Ty3d+9ehgwZAkBiYiIJCQk+berq6igtLbXPmZGRQW1tLWVlZXab4uJijDH2NgEZGRns378fd/MGbU1fZ8CAAcTGxraqFhEREene/F6uP2XKFDZv3szWrVupqKhgyZIl1NfXM27cOAAWL17M8uXL7faTJ09mz549rFu3jmPHjlFQUMChQ4eYNGkS4J1tP3nyZFavXs2OHTs4evQoixcvpk+fPvYqtZSUFEaMGEFeXh6lpaUcOHCApUuXMmbMGPr27QvArbfeitPp5IUXXqC8vJyioiI2bNjgM9n6SrUEUmRkJDNmzOiSy3ZyeeqH4KB+CA7qh+CgfuhalmnDRcuNGzfy6quv4nK5SEtL4+GHH7ZHXX7yk5/Qv39/vve9T1ZZbdu2jRUrVlBZWUlycvJlN3jctGkTdXV1ZGZm8u1vf5sBAwbYbWpqasjPz/fZ4DEnJ+eyGzz27t2bSZMmXbKFwJVqERERke6rTcFIREREJBz5fSlNREREJFwpGImIiIg0UTASERERaaJgJCIiItJEO18HiY0bN7Ju3TpcLhepqank5OTYezRJx3v//fd59dVX+fDDDzl79iw//OEPGTVqlH28eaXk5s2bqa2tJTMzk1mzZpGc3LE3Nezu1qxZw7vvvsuxY8eIiooiIyODBx980GdF6sWLF1m2bBlFRUU0NDQwfPhwZs2apY1ZO9Drr7/O66+/TmVlJeDdImXGjBnceOONgPogEF555RWWL1/O5MmT+da3vgWoH7qKRoyCQFFREcuWLWPGjBksWrSI1NRU5s+ff8ku3dJx6uvrSUtL49vf/naLx9euXcuGDRvIzc1lwYIF9OjRg/nz53Px4sUurjS8vf/++3z5y19m/vz5zJs3j8bGRp5++mkuXLhgt/nd737Hzp07+cEPfsBPf/pTzp49yy9+8YsAVh1++vbty/33388zzzzDwoULueGGG/jZz35GeXk5oD7oaqWlpbzxxhukpqb6PK9+6BoKRkFg/fr1TJw4kfHjx5OSkkJubi5RUVFs2bIl0KWFrRtvvJGZM2f6jBI1M8ZQWFjItGnTyM7OJjU1lccff5yzZ8+yffv2AFQbvp588knGjRvHwIEDSUtL43vf+x6nT5+2d7mvq6vjzTff5Jvf/CY33HAD6enpPPbYYxw8eJCSkpIAVx8+br75ZkaOHElycjIDBgzgvvvuIzo6mg8++EB90MUuXLjAr371Kx555BF69eplP69+6DoKRgHmdrspKysjKyvLfs7hcJCVlaVf9gA5deoULpeLYcOG2c/FxMQwePBg9Uknq6urA7Bv41NWVkZjY6PP++Pqq6+mX79+6otO4vF4ePvtt6mvrycjI0N90MWWLFnCjTfe6PP/H9B7oStpjlGAVVdX4/F4LrlGnJCQwPHjxwNTVDfXfPPh+Ph4n+fj4+MvuVGxdByPx8OLL77I0KFDGTRoEODtC6fT6fOXM6gvOsPRo0d58sknaWhoIDo6mh/+8IekpKRw+PBh9UEXefvtt/nwww9ZuHDhJcf0Xug6GjESkaCQn59PeXk5//zP/xzoUrqlAQMG8Oyzz7JgwQLuuOMOfv3rX1NRURHosrqN06dP8+KLLzJ79myioqICXU63phGjAIuLi8PhcFyS+F0ul1YaBEjzz72qqoo+ffrYz1dVVZGWlhaYosJcfn4+u3bt4qc//SlXXXWV/XxCQgJut5va2lqfv5Srqqr0/uhgTqeTpKQkANLT0zl06BCFhYWMGTNGfdAFysrKqKqq4sc//rH9nMfjYf/+/WzcuJEnn3xS/dBFFIwCzOl0kp6eTnFxsT0R2OPxUFxczKRJkwJcXfeUmJhIQkIC+/bts4NQXV0dpaWl3HHHHYEtLswYY1i6dCnvvvsuP/nJT0hMTPQ5np6eTkREBPv27eOWW24B4Pjx45w+fZqMjIxAlNxteDweGhoa1AddJCsri5///Oc+z/3mN79hwIABTJ06lX79+qkfuoiCURCYMmUKv/71r0lPT2fw4MEUFhZSX1/PuHHjAl1a2Lpw4QInT560H586dYrDhw8TGxtLv379mDx5MqtXryY5OZnExERWrFhBnz59yM7ODmDV4Sc/P5+33nqLH/3oR/Ts2dMeOY2JiSEqKoqYmBgmTJjAsmXLiI2NJSYmhqVLl5KRkaF/DDrQ8uXLGTFiBP369ePChQu89dZbvP/++zz55JPqgy7Ss2dPe25dsx49etC7d2/7efVD17CMMSbQRYh3g8dXX30Vl8tFWloaDz/8MEOGDAl0WWHrvffe46c//eklz992221873vfszd43LRpE3V1dWRmZvLtb3/bZ+NBab+vfe1rLT7/2GOP2X8YNG9q9/bbb+N2u7WpXSf4zW9+Q3FxMWfPniUmJobU1FSmTp1qr4xSHwTGT37yE9LS0i7Z4FH90LkUjERERESaaFWaiIiISBMFIxEREZEmCkYiIiIiTRSMRERERJooGImIiIg0UTASERERaaJgJCIiItJEwUhERESkiYKRiIiISBMFIxEREZEmCkYiIiIiTf5/vimWSvee2y0AAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(cox._breslow_base_hazard,)\n",
"#zooming in\n",
"plt.ylim(0.,.0001)"
]
},
{
"cell_type": "markdown",
"id": "6ca020f0-19d0-4d30-a844-2b7d459d09c4",
"metadata": {},
"source": [
"The ParametricDiscreteTimePH estimator allows us to estimate a proportional hazards model with its coefficients and base-hazard concurrently. ParametricDiscreteTimePH has multiple options for the base hazard distributions. The typical hazard distributions are implemented: namely, \"weibull\",\"log_normal\",\"log_logistic\",\"gamma\" and \"gompertz\". The flexible “chen” and “additive chen weibull” are present as well. It should be noted that ParametricDiscreteTimePH is far more expensive to train than Cox, and cannot be blindly used in most situations. Be careful to set 'pytensor_mode' to 'NUMBA' or'FAST_COMPILE' when using multiprocessing, as JAX is not multiprocessing safe. ParametricDiscreteTimePH is using pymc under the hood, and we have access to both the coefficient and base hazards priors via the ‘coef_prior_normal_sigma’ and ‘base_harard_prior_exponential_lam’ parameters, respectively."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "b1b2ddf1-1188-4e03-8c84-971ffb703ee0",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
ParametricDiscreteTimePH()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.