doc/user_guide/random-survival-forest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Random Survival Forests\n",
    "\n",
    "This notebook demonstrates how to use [Random Survival Forests](https://scikit-survival.readthedocs.io/en/latest/api/generated/sksurv.ensemble.RandomSurvivalForest.html#sksurv.ensemble.RandomSurvivalForest) introduced in [scikit-survival](https://github.com/sebp/scikit-survival) 0.11.\n",
    "\n",
    "As it's popular counterparts for classification and regression, a Random Survival Forest is an ensemble\n",
    "of tree-based learners. A Random Survival Forest ensures that individual trees are de-correlated by 1)\n",
    "building each tree on a different bootstrap sample of the original training data, and 2)\n",
    "at each node, only evaluate the split criterion for a randomly selected subset of\n",
    "features and thresholds. Predictions are formed by aggregating predictions of individual\n",
    "trees in the ensemble.\n",
    "\n",
    "To demonstrate Random Survival Forest, we are going to use data from the German Breast Cancer Study Group (GBSG-2) on the treatment of node-positive breast cancer patients. It contains data on 686 women\n",
    "and 8 prognostic factors:\n",
    "1. age,\n",
    "2. estrogen receptor (`estrec`),\n",
    "3. whether or not a hormonal therapy was administered (`horTh`),\n",
    "4. menopausal status (`menostat`),\n",
    "5. number of positive lymph nodes (`pnodes`),\n",
    "6. progesterone receptor (`progrec`),\n",
    "7. tumor size (`tsize`,\n",
    "8. tumor grade (`tgrade`).\n",
    "\n",
    "The goal is to predict recurrence-free survival time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from sklearn import set_config\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import OrdinalEncoder\n",
    "\n",
    "from sksurv.datasets import load_gbsg2\n",
    "from sksurv.preprocessing import OneHotEncoder\n",
    "from sksurv.ensemble import RandomSurvivalForest\n",
    "\n",
    "set_config(display=\"text\")  # displays text representation of estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to load the data and transform it into numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = load_gbsg2()\n",
    "\n",
    "grade_str = X.loc[:, \"tgrade\"].astype(object).values[:, np.newaxis]\n",
    "grade_num = OrdinalEncoder(categories=[[\"I\", \"II\", \"III\"]]).fit_transform(grade_str)\n",
    "\n",
    "X_no_grade = X.drop(\"tgrade\", axis=1)\n",
    "Xt = OneHotEncoder().fit_transform(X_no_grade)\n",
    "Xt.loc[:, \"tgrade\"] = grade_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the data is split into 75% for training and 25% for testing, so we can determine\n",
    "how well our model generalizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_state = 20\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(Xt, y, test_size=0.25, random_state=random_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training\n",
    "\n",
    "Several split criterion have been proposed in the past, but the most widespread one is based\n",
    "on the log-rank test, which you probably know from comparing survival curves among two or more\n",
    "groups. Using the training data, we fit a Random Survival Forest comprising 1000 trees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomSurvivalForest(min_samples_leaf=15, min_samples_split=10,\n",
       "                     n_estimators=1000, n_jobs=-1, random_state=20)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rsf = RandomSurvivalForest(\n",
    "    n_estimators=1000, min_samples_split=10, min_samples_leaf=15, n_jobs=-1, random_state=random_state\n",
    ")\n",
    "rsf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check how well the model performs by evaluating it on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6759696016771488"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rsf.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives a concordance index of 0.68, which is a good a value and matches the results\n",
    "reported in the [Random Survival Forests paper](https://projecteuclid.org/euclid.aoas/1223908043)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting\n",
    "\n",
    "For prediction, a sample is dropped down each tree in the forest until it reaches a terminal node.\n",
    "Data in each terminal is used to non-parametrically estimate the survival and cumulative hazard\n",
    "function using the Kaplan-Meier and Nelson-Aalen estimator, respectively. In addition, a risk score\n",
    "can be computed that represents the expected number of events for one particular terminal node.\n",
    "The ensemble prediction is simply the average across all trees in the forest.\n",
    "\n",
    "Let's first select a couple of patients from the test data\n",
    "according to the number of positive lymph nodes and age."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>estrec</th>\n",
       "      <th>horTh=yes</th>\n",
       "      <th>menostat=Post</th>\n",
       "      <th>pnodes</th>\n",
       "      <th>progrec</th>\n",
       "      <th>tsize</th>\n",
       "      <th>tgrade</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>119</th>\n",
       "      <td>33.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>574</th>\n",
       "      <td>34.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>421</th>\n",
       "      <td>36.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>76.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>65.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>80.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226</th>\n",
       "      <td>72.0</td>\n",
       "      <td>1091.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      age  estrec  horTh=yes  menostat=Post  pnodes  progrec  tsize  tgrade\n",
       "119  33.0     0.0        0.0            0.0     1.0     26.0   35.0     2.0\n",
       "574  34.0    37.0        0.0            0.0     1.0      0.0   40.0     2.0\n",
       "421  36.0    14.0        0.0            0.0     1.0     76.0   36.0     1.0\n",
       "24   65.0    64.0        0.0            1.0    26.0      2.0   70.0     2.0\n",
       "8    80.0    59.0        0.0            1.0    30.0      0.0   39.0     1.0\n",
       "226  72.0  1091.0        1.0            1.0    36.0      2.0   34.0     2.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_test_sorted = X_test.sort_values(by=[\"pnodes\", \"age\"])\n",
    "X_test_sel = pd.concat((X_test_sorted.head(3), X_test_sorted.tail(3)))\n",
    "\n",
    "X_test_sel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predicted risk scores indicate that risk for the last three patients is quite\n",
    "a bit higher than that of the first three patients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     91.477609\n",
       "1    102.897552\n",
       "2     75.883786\n",
       "3    170.502092\n",
       "4    171.210066\n",
       "5    148.691835\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(rsf.predict(X_test_sel))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a more detailed insight by considering the predicted survival function. It shows that the biggest difference occurs roughly within the first 750 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8iUlEQVR4nO3de3hU5bnw/+9NAkQgnEM4KYdigCCIGsBzQ6sVkf2jWgHFqsVW2r3rb1dtu/W1u1Jpt7V9q9VWW4tWaq0I0urGRkrF1lhF5KBSzgkpEAnhHIYQIJCE5/1jzUpWJnNYk5mVmczcn+vKRWbNmpVnlnHuPIf7fsQYg1JKqfTVIdENUEoplVgaCJRSKs1pIFBKqTSngUAppdKcBgKllEpzmYluQLT69u1rhg4d2uL4iRMn6Nq1a9s3KMnofbDofbDofWiS7vfio48+OmyMyQn2XLsLBEOHDmX9+vUtjhcXF1NYWNj2DUoyeh8seh8seh+apPu9EJHyUM/p0JBSSqU5DQRKKZXmNBAopVSaa3dzBEoplSh1dXVUVFRQW1ub6KaElJWVxeDBg+nYsaPr12ggUEoplyoqKsjOzmbo0KGISKKb04IxhiNHjlBRUcGwYcNcv86zoSEReUFEDorI5hDPi4j8QkTKRGSjiFzsVVuUUioeamtr6dOnT1IGAQARoU+fPlH3WLycI/gdMCXM89cD5/u/5gK/9rAtSikVF8kaBGytaZ9nQ0PGmH+IyNAwp0wHfm+sOtgfikhPERlgjNnnRXvuff5aKs8eanasukNPjmb0AaBv187069456GunDp/KjLwZXjRLKaUSLpFzBIOAPY7HFf5jLQKBiMzF6jWQm5tLcXFxi4vV1NQEPW77zDvnMbRhoOOIwZxTxweXbuWsAU5DzdHmrzndsSf7M0/h8/nIqQyakJd0It2HdKH3waL3oUk87kWPHj04fvx4fBrUSitXruSBBx6goaGBO++8k/vvv7/FObW1tVG910QGgmD9l6C75BhjFgALAAoKCkyw7MBIWYO7//AGx0/6Gh+fafDR6VQ2SzsN5sDxWg7XnG52/pgzmwCYOXAYDXK83WQkpnv2pE3vg0XvQ5N43Itt27aRnZ0dnwa1QkNDA9/97ndZuXIlgwcPZsKECcyYMYP8/Pxm52VlZXHRRRe5vm4iA0EFcK7j8WCg0qsf9pXnn2j2+Jez76Ku4ShLPh0LwOgrChl3TdOUxpqlj9Ntx+tknT1M7anDXjVLKaVcW7t2LSNGjGD48OEA3HLLLSxbtqxFIIhWIgPBG8A9IrIYmAQc82p+IJhzOg6lnkwO76nhzMn9nPCdaRYIJs34NvBtaheMZ09mHTMXjKdjRgc6dc9l6oV36ZyBUmnukT9vYWtldVyvmT+wO/P+bUzI5/fu3cu55zb9/Tx48GDWrFkT88/1LBCIyCtAIdBXRCqAeUBHAGPMs8ByYCpQBpwE5njVlmDyOnRgj0ykg3Rhn/yDo/vKeeZr9zQ75zOXXEFBz4lwfC1nzxo6nT1OyfHT8M8XNBAopdpcsD3m47GKyctVQ7dGeN4A3/Tq50cydkoeQ4qKoBY+PC4c7Naj2fO1x3ezpXg3Q/IvYMoVj/E//+rB0H2vQf//5Uz1gQS1WimVLML95e6VwYMHs2dP0xqbiooKBg4cGOYV7qRtZnGvWTPpNWum9eD2O4B6hjz/dOPzLz6wgKOVn1BZUsYJ3xmm33Ivy7iJBvMWezucZs7vCqwTu+boUJFSqk1MmDCBHTt2sGvXLgYNGsTixYtZtGhRzNdN20AQqHb7dspvv6Px8bDMPLoM+iJ79rzO0cpSxny8iiX/fiv/+atxdKjbBJ2A2mOsl9OsXz0fQIOBUspTmZmZPP3001x33XU0NDRw1113MWZM7D0TDQRA92nTWhzL3fwGQ+pLyb7yerYUv8y/PloF3MqBjP9gW2U1+QO68/lTy5l64rfM79uH+RoMlFJtYOrUqUydOjWu19RAQMAwkV/57XdQu307o4GyrEGNx6ePb/r+0QOXcmvGCR5mMfP79mH5zuUaCJRS7Y4GghDsXkLt9u2Y8y+gru4ISx55kAuuKGT2161lpovWfMpDr8P0Yx9Q0PUIHN+fyCYrpVSraCAIwe4llN9+B1lncjnRsWPjxLGdbzB70nkALHvjcuDPlJzYy5wV1ipYrU+klGovdIeyCLpPm8awg4fIPeezSEYOp46fafb87EnnsfO8GYyv7cVIOgGw/sB65q+ez9LSpYloslJKRUV7BBH0mjWT84uKOHfDk7x1/gXU1viaJZ595pIroMNQJlV14N9rKulkcllKb+Z3qNIJZKVUu6A9Ahe6T5tG1qhR9KtpoFOHno3Ha2v2sfX9d5lyTjeWNVzOjg5DAZhBNx72nQRg+c7lCWixUkq5pz0CFxpXFQUknr1w3/0crSylw5bV7Bw+g1v2XUf+me4A/EL+m4JTn+oEslIqru666y6Kioro168fmzcH3QAyahoIYlBwwxdY+VwpRys/4XPnXED/jl3YD6zZVcWTGRfCuZ+yzTGB7KSTyUqp1vjKV77CPffcwx133BH5ZJc0EETJmYHcA+iRPYBjx/dQW76WvM7jKMzqyueyerG1443kn/gHJzrUt7jG+gPrWX9gPaDzB0qp6Fx99dXs3r07rtfUQBCFwAzkk+vWcW7vbI6d24/6k2+TmVFG5TaoO93AxQMvokCy+I+anXTteNB6wdiboWAOS0uXMn/1fE1AU6o9+8uDsH9TfK/Zfyxc/1h8r+mCBoIoBGYgH13yKsybR6ehQzk4cnjj8cqSMmqObGRV/8mcPNNAl33HGFq3k5rjteQWzGFG3gydRFZKJQ0NBDHoNWsm1UVF9F+3jvHTb24MEvby0m6X383PNlg1Qb6z737yj26DhTdYPQOlVPuWgL/cvaKBIEbdp03j5Lp17J83j+qiIgDONmRypuEgvf0VSwEW/Pwaupx6hzHl70P5+zDmckrqjjFnxRydOFZKJZTmEcSo16yZ9H/kEbpMmNB4rF9NA4C/Yqnlb12mcsuZ77Ogx38CMPXTfzKgoUtjFvKcFXM0E1kpFdGtt97KZZddRklJCYMHD+a3v/1tzNfUHkEctKheevsdvNHQr9k5dtXSv2ENFc099gvyT9ey9eaHWb5zOSVVJYCuIlJKhffKK6/E/ZoaCNrI7EnnNRapg8vY8ug7gPXBPyNvBnNWzKGkqkSHipRSbU4DgUdMQz2mrs71+VOHWz0FO8cg1KoiDRJKqXjTOQIP2PkGp2v3svHtFSHPG3NmE2uWPg5YPYOFUxby8GUPU5BbEPR8rWqqlPKC9gg80GvWTLos386x2rWsf/Otxv0LnGrOvxG2bMJsWsqi82Y0DhvZQ0XBaCKaUsoL2iPwSF6HDkjmYI7uK+eF++5nySMPNusdTJrxbQ70LuDSDtuo+eA5V9eckTeDgtyCxrkE7RkopeJBA4FHhtSXcl6HjkhGDr6DJ6ksKWP9m281Oyf38i8D1gqiBT//PovWfBrxulOHT2Vk75G67FQpFTcaCDw0gWNc943vM+ziuRjpy9F95TzztXt45mv38MJ997PRN4A1Yx4GrGDgpmcQOJdQUlWi5SqUShN79uxh8uTJjB49mjFjxvDUU0/F5boaCDx0ct06Blau4sZvX0z+lZ8lq9sAAE6fqudoZSkrn3ua3VsPsfFcK8nsilPvuL62HRDs3oH2CpRKfZmZmTz++ONs27aNDz/8kGeeeYatW7fGft04tE0F4Sw9ATDl328FrHITW97by9svLCUzo4xDu3cBw8joNJYxZzbB+oVQ0HL/glCmDp/aOEwEkENOvN+KUipJDBgwgAEDrD8os7OzGT16NHv37iU/Pz+m62og8Iidabx/3rzGOkTdp02j16yZjLlqEKVrr6Zyxzi6ZFv1iVadM9kKBJv+GFUgsFcP2auJ7sy6M/5vRinVwk/W/oTtVdvjes1RvUfxwMQHXJ27e/duPvnkEyZNmhTzz9VA4CE7GFQXFVG7fXuzY3kTc6nc4cN30Nrb+G/9p3PFqXcY04qfo2WtlUovNTU1fOlLX+LJJ5+ke/fuMV9PA4HH7DpE5bff0Wx3s27AxIu/xHvr4WhlKTnZH3OyQwMnPv2Eml98ntzsrKaL+De0UUolD7d/ucdbXV0dX/rSl7jtttu46aab4nJNDQRtJHB3s9rt2+nHn+gz+EoO765g5IkyPhl8DQ1HV8LhExyuOQ3QbEObSEqqSniqw1O8uOJFLUWhVAoyxvDVr36V0aNHc//998ftuhoI2khghdLy2+/g5Lp1dPnibWRlD6Vv987Muu+HLFpzN8s27G08776993Fp1XrWLH2cSTO+HfL6dq0in8/XrF6RBgSlUseqVat46aWXGDt2LOPHjwfg0UcfZerUqTFdVwNBgtiris7s3oXp3FScrnmVUlizdAZsmY/ZtJRZVZcDVklr5znQVJqiuLiYQwMPsXzn8saAYD+vlGrfrrzySowxcb+u5hEkiL2hDYCpq6Ni6+agBersUhRjMz7l4SPf5Tv77o+YeOZMOgN0IlkpFZYGggTqNWsmHbK7c7bzOAC2rSoOel7u5V+m63kXMWZAD/Kl3HXimV2bSCmlwtFAkGCD6neS2XkcWdlDObR7V4vidIC1YmjOmzDnTXZ3HJ6YhiqlUpangUBEpohIiYiUiciDQZ7vISJ/FpF/isgWEUm7NZJD6kvp6SulS6dh5AwdRsXWzax87umw+xicPNPArN+sdlWkDtASFEqpsDwLBCKSATwDXA/kA7eKSGAe9DeBrcaYC4FC4HER6eRVm5KRvaw0q24As+Y9xrV33wOEHibq260zXTplsGZXFQ+9viliQLBXE+mGNkqpULzsEUwEyowxO40xZ4DFwPSAcwyQLSKClWNVBdR72KakY88T2MZdM4XB+ReEPD83O4sxA3rw6I1jmTSsN1v3VTdbbhpoRt6MxkljDQZKqWDEi6VIACJyMzDFGPM1/+PbgUnGmHsc52QDbwCjgGxgljHmzSDXmgvMBcjNzb1k8eLFLX5eTU0N3bp18+KteG7Pi59S3XkwQ3qW0W1KHiXLrPc3cvotLc4d/8n36HlsMyV5/8G+gdfx4zWnAPg/k84BQt+HVcdXsbjKuu4tvW/hiuwrvHo7SaE9/z7Ek96HJvG4Fz169GDEiBFxalH0amtrmTJlCmfOnKG+vp7p06fzve99r8V5ZWVlHDt2rNmxyZMnf2SMCbp6xMs8AglyLDDqXAdsAD4HfAZYKSLvGWOqm73ImAXAAoCCggJTWFjY4sLFxcUEO94erPtkGWtLoNw3gomfHKNnz55UbN1M7/ralttcdvsaFN3LyNJfMTIvj1/3HAVAYeFlQOj7UEgheaV5zF89n8VViynrWJbSyWbt+fchnvQ+NInHvdi2bRvZ2dnxaVArdOvWjXfffZdu3bpRV1fHlVdeyRe/+EUuvfTSZudlZWVx0UUXub6ul0NDFcC5jseDgcqAc+YArxlLGbALq3eQVibcN52JI48DULb1BKOvKARCzBMUzIFpT1rfb/pjVD/HHiYqyC1oLF2tQ0VKtR8i0tirqauro66uDmtkPTZe9gjWAeeLyDBgL3ALMDvgnE+BzwPviUguMBLY6WGbktaE+6ZTdvciqk0Ptiypokf2gMYksxa9goI5VhDYv4mHzXf545nLmPUbK+N4YISfY2cgLy1dyvzV8xvLVwMp3UNQKt72P/oop7fFtwx159Gj6P/QQ2HPaWho4JJLLqGsrIxvfvObcSlD7VmPwBhTD9wD/BXYBrxqjNkiIt8QkW/4T/shcLmIbAL+BjxgjDnsVZuS3Yj8rnSXY1SbHlBndaZCrR5i7M3Qfyznn93NzZ1WR5w0DuTsHQC65aVS7URGRgYbNmygoqKCtWvXsnnz5piv6WmtIWPMcmB5wLFnHd9XAl/wsg3tyYT7pjMBeOXuRVSdcyW9evtCn1wwBwrm0GnhDYzZv4nFnX7IqpOTsTpY7ti9A4A5K9IuhUOpmET6y91rPXv2pLCwkBUrVnDBBaFXGrqhmcVJaFC9NTp2supE5JP9PYOhdTsZd3QlxXvqIr9GKdUuHTp0CJ/PB8CpU6d4++23GTUq9mnViD0CEbnAGBN730O5NnZKHnteL+VAhoSeJ7D5ewY1v/g8HD7B77acYdtvVjc+HaxSqVKqfdq3bx933nknDQ0NnD17lpkzZzItYK+T1nAzNPSsP9v3d8AiY4wv5p+qwuo1aybnrvgRhxlFfX0FK597GiB0MMBKNMutWs+92X9nNTcDsHWftQrXbSAoqSppMUSkE8hKJY9x48bxySefxP26EYeGjDFXArdhLQVdLyKLROTauLdENTN2Sh59T2XRI2siEGbSuPEF1of/V7p8yJKvX8aSr19G/oDurNlV5aom0dThUxnZe2SzY/YS0zkr5jBnxRxdaqpUinI1WWyM2SEi/w2sB34BXOQvC/GQMeY1LxuYrnrNmkmHtxdxOiOPHh33RH6BvaTUP34I1rDQml1VLNuwN2KvwDlxbFtaurRxJVFJVUnjeUqp1BKxRyAi40Tk51hLQD8H/JsxZrT/+5973L60NiK/KwCnTjaELlEdxuxJ5zFpWO9W/3x7g5uFUxa26C0opVKHmx7B08BzWH/9n7IPGmMq/b0E5RE7yexg53F0zdhCxdbNVGy15u1DzRf0PLYZ1i+0eghxFmwOwaZzCUq1X26Wj75mjHnJGQRE5FsAxpiXPGuZAqxeQWbncXSVyRFLVNvzBBTdCwtvgIU38PmTy13PE4QTbA7BpnMJSrVvbnoEdwBPBhz7CvBUvBujWrJ7BWD1AratKg5beqKktJSRpzdaj/dvYnrXWh7lUlfzBOEEm0OwOecS1h9Yz/oD6xtfo5RKfiF7BCJyq4j8GRgmIm84vt4BjrRdE5VT2IJ0wL6B1zVua0n/seRmZzFpWO+49ApCcc4l2HsfaLkKpbzT0NDARRddFJccAgg/NPQB8Diw3f+v/fVtIPSCduWJsydPUn77HZx7pDrsxjXBTB8/CCCqWkStNSNvBgW5BY3zCTpMpFT8PfXUU4wePTpu1ws5NGSMKQfKgcvi9tNUq2T26UNVh968f+xycp99i/p8w/7DB8JnHNvK32f22L+xbNioxl6B15nG9vaY9jCRs3egk8pKxaaiooI333yT733vezzxxBNxuWbIQCAi7xtjrhSR4zTfUEYAY4zpHuKlKs7yb7iA0rUHqNzRCV/PPHKrXwes4aGwgWDszVD+Pmz6I9PH/7pxn2Nwn23cGs5S184g4AwMGhBUe/feq6Uc3lMT12v2PbcbV83MC3vOvffey09/+lOOHz8et58brkdwpf/fxG3HowAYc9Ugxlw1iC3v7aX45RLqsi5j8IiujbkFttFXFEJmVtML7SQzmj74H3p9Ew+9vollG/Z6XococILZDgyanKZU6xQVFdGvXz8uueQSiouL43bdcD2CsJlIxpiquLVCuTLmqkFs/P27QNOksc3OMTjvs9dCiO347A/9ZRv2smZXFWt2VTU77jU7MGjJa5UKIv3l7oVVq1bxxhtvsHz5cmpra6murubLX/4yf/jDH2K6brjlox9hDQmF2nt4eEw/WbVaVUZ/Tm86zqx5jzUe2/j2ClY+9zRVO7aFfe3sSecxe9J5LFrzabPeAbRtpVJncpoOEynlzo9//GN+/OMfA9YezD/72c9iDgIQfmhoWMxXV3E3Ir8ra0usvY0nOI7bOQY+R62hRuXvt8g2dvYOIPpKpbGwJ5NB5w2USgbh8ghG+f+9ONhX2zVROU24bzq9G/ZTldGfLe+5WA5qZxsH2eh+9qTzWlWpNFaBeQf2clPNPVDKvcLCQoqKiuJyrXBDQ/cDc7FyBwIZrKJzKgEG1e+kKqM/W9/czJirBjV7rqayovmyUnvCOEivwMmuVOocKgp2Trx7DDpvoFTihRsamuv/d3LbNUe5Ye9gVt/QfAXv6CsKqdi6ueVGNvYy0qJ7rcdBgkHgUFGgthg6Wn9gPUtLl+rwkFJtzM1WlVnAfwBXYvUE3gOeNcbUetw2FYK9V0GgcddMoaS0hE/fXcnK555m26piRl9RyLhr/B/8RfdGDAahPuhn/WY1W/dVM8u/DWa8ewdTh09tnCvQQKBU23JTdO73wHHgl/7HtwIvAfp/a4JVmx68cvciMvv0ITMnh7yJueTkX8jIvJGNxemaylY7gsGmP0ZdptouUwE0Lj0N7D3EEhxm5M3QOQKlEsRNIBhpjLnQ8fgdEfmnVw1S7ozI70rZ1mOcPV7NmePVVB011B86RM6Uzoy7ZgrjrpnSuKS0MQPZni/Yv8kqUz32ZtcBwdlbWLTm0xZBwBkc2nIZqlIqdm4CwScicqkx5kMAEZkErPK2WSqSCfdNZwJwdMmrVBcV8X5NJgdlOBllBgqtc4KWrbZXEZW/b31B1L2DYENIdnAI7C1oUFAq+YVbPrpJRDYCk4APRGS3iOwCVgNXt1UDVXi9Zs1kyEu/59xO+wE4vtHX7PkWZasL5lglqqc9aT12bGLD+oWtboe9FPXRG8c2bo+5dV91m1Q8VSqdDB06lLFjxzJ+/HgKCgrics1wPYL4FLpWbcJeSXT2nC7Njtu9ghbsXoCdX7B/U/PjreTsLdgTy9EI3A5Tk8yUaumdd96hb9++cbtepDLUjUSkH5AV4nSVYPZKorP19UGfD7qrWcGcpg/+hTc0zR04RTGPEIxzpZFTsCEjZ8YxBC9jHUgDhVKxc7N89P/DSiobCBwEhgDbgDHeNk21RnXnwWx5b2+zRDM7vyBs2Wp77sDJnkewew1RBgXnSiOnUDkJoaqVhqJVTFUivfO7BRws3xnXa/YbMpzJX5kb9hwR4Qtf+AIiwte//nXmzg1/vhtuJot/CFwKvG2MuUhEJmMtIVVJxs44Ln7Z+oC0g4E9PGSXrbZyC1rud9ziQ379wqYg0IqgECovITAnwRbYSwi3TzLAnBVzNAlNpZ1Vq1YxcOBADh48yLXXXsuoUaO4+urYpm3dBII6Y8wREekgIh2MMe+IyE9i+qnKE0PqSzmzezclI2dT/HIJpWsPkDcxlzFXDWqcNG6eWxBhdzNncAgXFCCq3kKwnkJrMpftJLT5q+c39hx0qEi1lUh/uXtl4MCBAPTr148bb7yRtWvXtkkg8IlIN6yM4pdF5CAQfCBaJVT3adMYNG8eAIdGfI5Du84CVs8gZG6BW6GCAkS9FDVYT2HWb1ZHvZWm/YFvBwHnnILP5+PFFS+6uo5Ng4hKdidOnODs2bNkZ2dz4sQJ3nrrLR5++OGYr+smEEwHaoF7gduAHsD8mH+yirtes2ZSWlrC+aU7GPTeD/l4/LeorM9rNmcQchVRNAKHkdYvjFi+IhK76N2yDXuj6hU4h48izSmEYwcR+5pKJaMDBw5w4403AlBfX8/s2bOZMiWKP+hCiBgIjDEnRKQ/MBGoAv5qjDkS809Wnjh11VUM+f73ObrkVXKffQtfz7ygVUqDriJqrQJH+Qq7hEWUE8uzJ50Xc86BMygUFxdTGGKntmCWli5tHGLSQKCS1fDhw/nnP+Nf2MHNqqGvAQ8Df8fareyXIjLfGPNC3Fuj4qbXrJlcAhyIpkppLJx5CcHmEMBVcAicRG6rzGStdaTSmZuhoe8CF9m9ABHpA3wAaCBIcnZugb2JjXN4CGDlc0/HPxgUzGk5hwCuVh0FTiInolxFYEKbTecPVCpzEwgqsKqP2o4De9xcXESmAE8BGcDzxpjHgpxTCDwJdAQOG2M+6+bayp1Qm9g4g0HUE8eRRFqKGiKLOXAS2VnczouKp4ECE9oam67baSoHYwwiwbZyTw7GmKhfEzIQiMj9/m/3AmtEZBnWfgTTgbWRLiwiGcAzwLVYwWSdiLxhjNnqOKcn8CtgijHmU3/2soqjUJvYQMv8AiB4jkE8tCKLOVLF03hvlhMqb8GehNYENpWVlcWRI0fo06dPUgYDYwxHjhwhKyu6IhDhegTZ/n//5f+yLXN57YlAmTFmJ4CILMYKIlsd58wGXjPGfApgjDno8trKpVCb2Njs/AKAQ7t3AXEaJgonWBZzhFpHoZactgXdTlPZBg8eTEVFBYcOHUp0U0LKyspi8ODBUb0mXK2hR5yPRSTbOmxqXF57EM2HkCqwKpk65QEdRaQYK/A8ZYz5vcvrqziw8wuAxl6B54INHS28IeK+yslAM5nTW8eOHRk2bFiimxF3blYNXYC1I1lv/+PDwB3GmC2RXhrkWODgVSZwCfB54BxgtYh8aIwpDWjDXGAuQG5uLsXFxS0uXFNTE/R4ugl2H+rr66nuPJiiBxfQbUpeyNf6fD5qKitY+qunyMm/MOR5XhjQeRwjeR/fe8+zocbd/2g+3ylKjp7lBy+tpPDcjs2e8+L3YUTdCNZjZTIv+tjqZRV0LeCK7Cvi+nPiSf+/aKL3IjQ3k8ULgPuNMe9A4+Tuc8DlEV5XAZzreDwYqAxyzmFjzAnghIj8A7gQaBYIjDEL/O2goKDABFsfHu268VQV7D6s+2QZa0ug3DeC3i9XMiK/KxPum97itb3ra1n53NOcPbSPwsJvtVGLbYWwcCM9wfV/x8pzPuWh1zex7WQ3flB4WbPnvPh9KKSQvNK8ZpnMZafLKOtYlrSTyPr/RRO9F6G5CQRd7SAAYIwpFpGuLl63DjhfRIZhTTjfgjUn4LQMeFpEMoFOWENHP3fVcuXahPumw8+XUbb1BNWmB2VbjzEhyHlBdzRra4GTyGFyD+wktGAF7Hy+U/y6JLo5BDcrkIJlMocql52swUGpQG4CwU4R+T7W8BDAl4FdkV5kjKkXkXuAv2ItH33BGLNFRL7hf/5ZY8w2EVkBbATOYi0x3dyaN6LCs7e2fOXulnkFTq5KVnslcBI5VGKaw/29P8cTETunkbVmBZIdFIKVttAVRqo9cRMI7gIeAV7zP/4H4Go2zxizHFgecOzZgMf/F/i/bq6nYhcqr8DmqmS1V4LVMAoTBNi/iUn9YcnXv93iKWsY4LIgLwquNUXvbMGWneoKI9WehA0E/lyApcaYa9qoPcpj4fIKbPaS0jZbThpKsNVFToF5CDGwi9499PqmuCWt6babqr0IGwiMMQ0iclJEehhjjrVVo5R3nGUn1v18WdBJY3tJ6ZJHHmyWbObUpj2FNmB/0AcGgdaWuQi37aYGBJVs3AwN1QKbRGQlcMI+aIz5T89apTw1Ir8ra0tgbUk2hAgG0DzZzCnhPQWPBEtac2Y0RzOPEGrbTZ07UMnITSB40/+lUoS9imhtSXbYYOBMNnNy9hQS3jPwOAnNGRziMY+gcwcqGbnZj+BFEekEjMJKCCsxxpzxvGXKU4HBoOLxjwEat7YMp1XbXnph7M1WILD3QLCPeRQUWrt5TiCtcKqSjZvM4qnAb7DqDQkwTES+boz5i9eNU96acN90Ts35EXszhwM9OVxhVQ+JFAhi3vYyXpx7IECz5abjfT7Y1TO660UIIsHyFqKdSA5V4VSHjFQiuRkaegKYbIwpAxCRz2ANFWkgSAFD6kvJ3fwGWbWj+CBrCvWH+rh+bVIkoIXbSzkaEYre2Zx7JsSSexBozoo5uspIJYybQHDQDgJ+OwGtEpoiuk+b1vj92ZMnORgm2SyYhCagBXIEhQ3RlhNwuRQ1cM4gXsKtMgp1vgYJFS9uAsEWEVkOvIo1RzADa2+BmwCMMa+Fe7FKbr1mzaTXrJkAlM/5ET6GU7r2gOtAYPcKUkICq5+GWmUUjA4jqXhzEwiygAOAvXPYIaxKpP+GFRg0EKSIIfWl7M0cTuUOouoVpIRgE8/O50IEB6/2WA41hASataziz82qIf2tSyM5ZX+nauRsil8uoXTtAVeriIDkWU7aWoETzzZ7Atp5jl8i91gOtvJIh4tUa7npEag00X3aNAbNm0enoUM5csFU16uIkmY5aaxC7bVcdK/1ZZ/j15o9liH2ABFs5ZEOF6lYaCBQjXrNmkl1URE5q19m7JQ8ihnB4YoaXvfnGEDwPIOkWU7qBfuDP0QwcIq0xzLEZ59lLXKn4k0DgWqm+7RpnFy3jv3z5tHnstuoz50AdAOI2ENIqYljpyiCgS1YuQpou32WlYpGyEAgIveHe6Ex5on4N0clmr2CqLqoiNzNbzCkvpQhj1nbSL/++MeNPYRwcwcJzSvwSmAwiCGT2Z5gjvccQrB5A5/Px4srXnR9DZ1nSE/hegTZbdYKlVTsJaXlt9/R7HjexFwAKnf4qNzhCzqZnFR5BfEWJpPZbUCwJ5jjPbEcKmM5GjrPkL5CBgJjzCNt2RCVnGq3b28MCN2Ai4E+mXkcyJ3QGBCgabgocGMbSLGS1cEymV1mJUPTkFGwiWX7+dYItdw0mn16g2U3R6I9iNTgptZQFvBVYAxWTgEAxpi7PGyXSgLOrGOnnNUvk8PLHLrsNjZ1vpzil62/JO1g4Cxfba8k2raqOLUCAjQFhYU3tNxrGSLut+ycWH7o9U0tNsXxcvlpMNH2KuzsZ9BeRHvnZrL4JWA7cB0wH7gN2OZlo1RycGYdOx1d8irVRUUMqS/lzO7dlIyc3WzrS2f56o1vr2jsIdjPpZzAvZYh6l4CNN8Uxzl01FYBIVwSWzBLS5cyf/V8lu9croGgnXMTCEYYY2aIyHR/SepFWBvSqzTlDBDdl7zKgTBbX4bb7SxlegjB8g+C9RJc9hCgaflpqHyEtu4tBDMjb0bjZjtuhpN0GCl5uQkEdf5/fSJyAbAfGOpZi1S74tz6MlxZisDdzpxDRs5zUiIwQMtegnNi2XlOhMAQLB8hWHBIVGBwO5wUqYheLD9fg0vs3ASCBSLSC/g+8AbWnOH3PW2ValcG1e+kKqN/i7kCp8DdzuwhI1vKBYbAXkJgiWyXgSHS9pkQn8nm1nI7nBSuiF5r6Sqn+HETCBYaYxqAd4HhHrdHtUNjp+Rx5tlFLeYKwnEbGOxz2z23gcE+N4xgw0j2RHOih4tCiXb+wQ3Npo4fN4Fgl4isAJYAfzfGGI/bpNqZXrNmcgmEnSuIJFhgSLlyFU7BAkNgsprN5c5pSrWWm0AwEqvk9DeBF0Tkz8BiY8z7nrZMtSvOuYJ1P19m7Ykcg2D5CIHa9dBRoEjVTyPsvPb5kxfzy2NXtihh4fOd4tcl7staJMMkdDTcTFTrPEJkbspQn8LalOZV/1zBU1jDRBket021MyPyu7K2BMq2nmBCHK4XOMHslJLLUUNVP420/Wb5+8zlfegBf6P1GcbxKIjXltxMVDsnqQPLbWiAaOKq6JyIfBaYBVwPrANaLi5XaW/CfdMpuzvyCiK3AoeLnEL1ElJOsOAQyD+sNLfnx8yd88NmT1mZxZe5+lHtrSCem3mHUJPUOtHcnJvM4l3ABqxewXeNMSe8bpRqv+wVRG4njVUcFMyJ3GtwKXDHNTeSeTjJGSyc5TZ0ork5Nz2CC40x1Z63RKWEsVPy2BPDpLFKnMAd19wItwFPsOsna8BId+HKUP+XMeanwP+ISIuVQsaY//S0ZapdsieNq00PXrl7EWDNHcQ6eaxcCFLvaLzPB7t6Wg9crD6K9oM61AY8gdrb/EO6CdcjsOsJrW+LhqjUMSK/K2VbjwFQbXpQtvVYXCaPA6XkvgetFazekZPL1UfR7q/gNnjM+s3qVg07RaK9jPgIV4b6z/5vNxpjPmmj9qgUMOG+6Y0f/K/c3dQ7iGfPwN73YOVzTwMptnqoNUJMKm+wx8XdrD6KolBetFoz7BRJpGEpDRLuuZkjeEJEBgBLsfIHtnjcJpVC7N5BvHsG9gd/SiedxZOb1UeBZbTjqDXDTpGEG5YKFiScORW7O1XTkLk37KRxOi0vdZNHMFlE+mMtGV0gIt2BJcaYH3neOtXu2b0De74gnlJ2j2TlSrjgEmnu4mTVOLr0Dn3tdFte6iqPwBizH/iFiLwD/BfwMKCBQEXl7MmTLba/DNR92rSgeyCEEyz72OfzceDdFc2OpVQmsleCbbBja8X+zIkSLEg4cypm/QY48wUWTgmeY5Fuy0vd5BGMxkomuxk4AiwGvu1xu1SKyezTh3qA2tDn1G7fDhBVIAiXfeyUckXsvBBuwtntZHOk67eTQJJuXFUfBV4BvmCMqYzm4iIyBaskRQbwvDHmsRDnTQA+BGYZY+KTGaOSSmZODr7T5/DxqHsbjwVufF9++x3N9kh20zsIlX0cuFdvyhexi4dw8whuJpvD8XAiWsUubCAQkQzgX8aYp6K9sP+1zwDXAhXAOhF5wxizNch5P0F3PUtpeRNzmz22N74vXXugMSA490g+uW4dJ9etA6LrIYSi8wkxcjPZHI6HE9EqdmEDgTGmQUT6iEgnY8yZKK89ESgzxuwEEJHFwHRga8B5/z/wJ/BkqblKEmOuGtTsr/8t7+2ldO2BZgEBRsD4ewGoH3SIXmuWwrx5VBcVAa2bPwhkzyfofEEChJp/0CGjhHMzNFQOrBKRN4DGOkPGmCcivG4QsMfxuAKY5DxBRAYBNwKfI0wgEJG5wFyA3NxciouLW5xTU1MT9Hi6aU/3odclYHoIx8oNPp+v2XO1JzpycuwN5HQ7yBmfj047dnBy3To+/cMfms6ZOIFTV10V9NrB7kOHnAF09Pka5wtWv/m/zZ7vff5ocvIvjMdbSxrJ8vswoPM4crN8EPDfueexzVD+Pr73nm/VdQ/kXs2+gde5Otd5L3y+UwAh7439+5gM964tuAkElf6vDkB2FNeWIMcCS1U8CTzg73mEvJAxZgGwAKCgoMA4x35tgWPC6ard3YfC4Idff/xjKndA5n/9mjFXDeLoklepLiqii//52u3b6Vm6gyHfD75ratD74H8cuBsaWJPJNZUVnD20r8W12nPvIXl+HwqDH/bPPfRszSX3b6Ln6Y2MLPyxq9Od98LOJwhVmdUuV50c9857bvIIHmnltSuAcx2PB2MFFKcCYLE/CPQFpopIvTHmf1v5M1WKyJuYS+UOH8Uvl7QYNgKozdrOoPqdDGnFtYNNMAcLDqCrjTwXy9zDwhvCL3eFsMNO4Upe7O5UTd+unVvXrnbIzfLRd2j5lzzGmM9FeOk64HwRGQbsBW4BZgdcY5jj5/wOKNIgoIDG+QQrCLRU3aG3yywYd0KtPrJXG9krjqB99xBSSqT6SmFWKkUqeXHydD2HW9uudsjN/0rfcXyfBXwJrCXh4Rhj6kXkHqzVQBnAC8aYLSLyDf/zz7aivSqNBE4wO3mRqRyM/YFvB4GU3BmtvYrUmwjoMTgrsc4GZl8SurcwaWEmtbIn6RLLRvUexQMTH4j7dd0MDX0UcGiViLzr5uLGmOXA8oBjQQOAMeYrbq6plO3s8WqOLnk1LstLw3H2FpY88mCzTGbtHSSxcD2GCHkNPRomptVmvG6GhpwVOToAlwD9PWuRUi5k9unDwYz+bFrxAVfParuf68xktucPtq0q1oCQjAJ6DBucE+cR8hp6NVxNr4arQ5agSDVuhoY+wpojEKwhoV3AV71slFKR5N9wAQdfLmHPqb5t0iuwOXsH9gSzDhelJi/2T4hV/sDuzPu3MXG/bodIJxhjhhljhvv/Pd8Y8wVjzPtxb4lSURhz1SD69azD1zOPTStKE9KGcddMYda8x8gZOizyyapdmT5+EPkD0me71XBbVU4A9vgrjyIid2BNFJcDPzDGVLVNE5UKzu4V7M0cnuimBK2AGkiHj9oPL/ZPSGbhhoZ+A1wDICJXA49hlYMYj5XcFWHtllLeGnPVIDb+/t1m5a3jUYYiWm4qoDrnE8JdRwOFSoRwgSDD8Vf/LGCBMeZPwJ9EZIPnLVPKBWd5a7tQXXVREb18Psp/+0KL870IFKFyEJxCJazZAgOFBgXVlsIGAhHJNMbUA5/HX+vHxeuUajPO8tb1gw5Rf+QIAPV968nMbP5r2q/iA0YUFbV5jwEiBwtnoNDJZ9XWwn2gvwK8KyKHgVPAewAiMgI41gZtUyoiZ3nrzJwcMnNyAKtoWLeePRufq9zho2rITVD+GoTYJS0Rw0q2wFwFpdpSyEBgjPkfEfkbMAB4yxhjl5nogDVXoFTChco+tgqMXdz4eMt7eyl+uYTtQ27iYMN+q05RfdNqI+ewUqBEBgil2kKk/Qg+DHIsMWv1lIqBs3ZR5Q6oyujPkQumNj5vDyuFChAQnw1ylEpGOtav0obde7A3xXHKzMnhoK9j0ABxZvcucp99i/PjuEGOUslEA4FKO6GGk8IFCF/PPDqd/oDczW8A3vcOAvMSdBWR8pIGAqX8wgWI4pdL2NT5cvaOH06/ig9aTDjHs5cQmJegOQjKaxoIlIrAOb9waNdZGHw5I2pXND5fu307EL9eQuBS00g5CLrcVMVKA4FSLti9hdcf/5jDFZ34eNS9jc/FsluaG5FyEJY88iAVWzez8e0VGgxUq2ggUCoKzrwFW7x3S4vW6CsKqdi6udkuavZxMrMS1i7VfmggUCoKweYR2mq3tFACd1GDpuGi3M9qD0FFpoFAqRQQOHxk76Tm8y2md32tDhmpsDQQKBUHdgXUZMkxsFceOYeMdGVRAMd+xkGNDb2ncarRQKBUjDL79KGqQ2/eP3Z5s8SzQG0ZJOwewtJfPcXZQ/tcLUGFNFqGGm4/Y4i4p3Gq0UCgVIzyb7jAX7qiE76eeRxq2N/inLPHq8MGiUDxCho5+RdSWPitiEtQwV2+QqB2GzgC9jNuIcKexqlGA4FSMWpZuqJni3Mqd/haBInAuka2eOclQHz2TAgUKXC02yCRhjQQKBUnoTKTgRZBonKHr0VdI5vXeQmhuAkWTuECRzS9i6QNGJHmEBLBo3kLDQRKtYHAIBGsrpEt0XkJboULHG57F60ZjnLyLIhEmkNIhPL3rX81ECiVGsL1HhKdlxAPbnsX0Q5HOXk6NBVpDiER/uLdhkUaCJRKQvZyVEjtstfRDkc5hQsiKVl/6frHPLu0BgKlkkxmnz7UA9TqxjjhhAsiut1ndDQQKJVk7L2Xh3x7NkeXvMr+efPYP29e0G00banca1De00CgVBKzP9zDBQEvlpuq9KKBQKkkdLiihtcf/9j/aAR5c38ccnK5PGCTHKWipYFAqSQTWOr6cEUNQMhAoFSsNBAolWQCl5Zam+E4ewjN1WZNSUgCmkodGgiUSnLBNsNxqu7Qm7NnTrYYIuo+bRrk9vOyaSpFaCBQKsmFSz4DWPrgX6k/2wVqm47ZE8h89S6PW6dSgQYCpdq5zJwcfKfPabGPcr+KD+j7+BOU//YFQJeYqtA0ECjVzgUbOjreuR8Mvpy++z8FmhLTwi1DhdQKFhVbN7Px7RWplV3sEU8DgYhMAZ4CMoDnjTGPBTx/G/CA/2EN8O/GmH962SalUk2woSNrgrkTH4+/l549e1I/6BD1R46EvU60eyYESqYgMvqKwsY6RBoIIvMsEIhIBvAMcC1QAawTkTeMMVsdp+0CPmuMOSoi1wMLgEletUmpdGH3Enw+H9CUrRxOsD0TnELtnwDuehxtvUNba4vZpSMvewQTgTJjzE4AEVkMTAcaA4Ex5gPH+R8Cgz1sj1Jpw+4lFBcXU1h4savXRNpYJ9T+CQD1gw6Re2BdqwJFMvUk0pUYY7y5sMjNwBRjzNf8j28HJhlj7glx/neAUfb5Ac/NBeYC5ObmXrJ48eIWr6+pqaFbt25xfAftk94Hi94HS7zuQ1WZ4Vh56M+Kk4esf7s4Oh09hgi9RwgA57z3Hllr17V4XacdOwA4c/75Mbcx0Nq+XTnbswcjp98C6O/E5MmTPzLGFAR7zssegQQ5FvQ3SUQmA18Frgz2vDFmAdawEQUFBaawsLDFOdZfPi2Ppxu9Dxa9D5a43YcIlwjcaKdyh4+ThwxyrAcAR7vcQN5/3dViLuPoklepLiqiS+wtbKZ2+3Y6dR1I1tAhje9ffydC8zIQVADnOh4PBioDTxKRccDzwPXGmPCzWUqppBRpB7bKHT4qd/iC7MoWvo5Sa5XffgfUVXNo967GktQ+n48D765I3q0xE8jLQLAOOF9EhgF7gVuA2c4TROQ84DXgdmNM8MFFpVS743ZrztABIja1WVPoymqyhnZtdty5q5kGhCaeBQJjTL2I3AP8FWv56AvGmC0i8g3/888CDwN9gF+JCEB9qDEspVT7FSo7OtzezbGo7tCbbnVjmDqqd+NEdHFxMb3ra9m2qjg1dzCLgad5BMaY5cDygGPPOr7/GtBiclgplR4ilc9oraUP/pUzx6x9HJwrkuxdzXQHs+Y0s1gplXIyc3I46OtI+enDIauyOucP2guvhrM6xP2KSimVYHZC3d7M4UGfH31FITlDh7Vlk2J2aPcuz5LktEeglEo5Y64axMbfvxvy+XAb3ycrL3sv2iNQSqk0p4FAKaXSnA4NKaVUO9BvSPD5jnjQQKCUSlnVpgev3L0IgPr6el55eVHjcyPyuzLhvumJalrUJn9lrmfX1qEhpVRKGpHfle5yLOhz1aYHZVtPtHGLkpf2CJRSKWnCfdOZ4HjsLDpn9xKURQOBUiotnT150ipO1454tXeDDg0ppdJOZp8+dOgS7+LX3qrdvj3intOtpT0CpVTayczJwXf6HD4edW+im+JabdZ2a7tQD66tgUAplXbsEhTtSVVGf896MRoIlFJpx6uqp15671XvtmzRQKCUUu3AVTPzPLu2ThYrpVSa00CglFJpTgOBUkqlOQ0ESimV5jQQKKVUmtNAoJRSaU4DgVJKpTkNBEoplebEGJPoNkRFRA4B5UGe6gscbuPmJCO9Dxa9Dxa9D03S/V4MMcbkBHui3QWCUERkvTGmINHtSDS9Dxa9Dxa9D030XoSmQ0NKKZXmNBAopVSaS6VAsCDRDUgSeh8seh8seh+a6L0IIWXmCJRSSrVOKvUIlFJKtYIGAqWUSnPtPhCIyBQRKRGRMhF5MNHt8ZqI7BaRTSKyQUTW+4/1FpGVIrLD/28vx/n/x39vSkTkusS1PHYi8oKIHBSRzY5jUb93EbnEfw/LROQXIiJt/V5iEeI+/EBE9vp/LzaIyFTHc6l6H84VkXdEZJuIbBGRb/mPp93vRMyMMe32C8gA/gUMBzoB/wTyE90uj9/zbqBvwLGfAg/6v38Q+In/+3z/PekMDPPfq4xEv4cY3vvVwMXA5ljeO7AWuAwQ4C/A9Yl+b3G4Dz8AvhPk3FS+DwOAi/3fZwOl/vebdr8TsX619x7BRKDMGLPTGHMGWAxMT3CbEmE68KL/+xeBLzqOLzbGnDbG7ALKsO5Zu2SM+QdQFXA4qvcuIgOA7saY1cb6BPi94zXtQoj7EEoq34d9xpiP/d8fB7YBg0jD34lYtfdAMAjY43hc4T+Wygzwloh8JCJz/cdyjTH7wPqfA+jnP54O9yfa9z7I/33g8VRwj4hs9A8d2cMhaXEfRGQocBGwBv2diFp7DwTBxvFSfT3sFcaYi4HrgW+KyNVhzk3H+2ML9d5T9Z78GvgMMB7YBzzuP57y90FEugF/Au41xlSHOzXIsZS6F63V3gNBBXCu4/FgoDJBbWkTxphK/78HgdexhnoO+Lu3+P896D89He5PtO+9wv994PF2zRhzwBjTYIw5CzxH0xBgSt8HEemIFQReNsa85j+svxNRau+BYB1wvogME5FOwC3AGwluk2dEpKuIZNvfA18ANmO95zv9p90JLPN//wZwi4h0FpFhwPlYk2KpJKr37h8qOC4il/pXhtzheE27ZX/w+d2I9XsBKXwf/O3+LbDNGPOE4yn9nYhWomerY/0CpmKtFvgX8L1Et8fj9zoca9XDP4Et9vsF+gB/A3b4/+3teM33/PemhHa+EgJ4BWvYow7rr7ivtua9AwVYH5T/Ap7Gn2HfXr5C3IeXgE3ARqwPvAFpcB+uxBrC2Qhs8H9NTcffiVi/tMSEUkqlufY+NKSUUipGGgiUUirNaSBQSqk0p4FAKaXSnAYCpZRKcxoIVEoRkT6OCpz7HRU5a0TkVx78vG+IyB0xvP53InJzPNukVLQyE90ApeLJGHMEq8wCIvIDoMYY8zMPf96zXl1bqbaiPQKVFkSkUESK/N//QEReFJG3xNrf4SYR+am/Hv0Kf9kCu0b9u/4Cf38NyN7Fca3v+L8vFpGfiMhaESkVkauCnC8i8rSIbBWRN2kqiIaIPCwi60Rks4gs8J/7GRH52HHO+SLykf/7x/zX2SgingU7lfo0EKh09RngBqzSxH8A3jHGjAVOATf4g8EvgZuNMZcALwD/4+K6mcaYicC9wLwgz98IjATGAncDlzuee9oYM8EYcwFwDjDNGPMv4JiIjPefMwf4nYj09l9rjDFmHPAj1+9cqQAaCFS6+osxpg6rLEMGsMJ/fBMwFOvD+gJgpYhsAP6b5oXJQrELn33kv06gq4FXjFUgrhL4u+O5ySKyRkQ2AZ8DxviPPw/MEZEMYBawCKgGaoHnReQm4KSLtikVlM4RqHR1GsAYc1ZE6kxTrZWzWP9fCLDFGHNZa64LNBD6/68WdV1EJAv4FVBgjNnjn9/I8j/9J6zexd+Bj/zzIIjIRODzWMUW78EKHkpFTXsESgVXAuSIyGVglTsWkTERXuPGP7AqYGb45xwm+4/bH/qH/fX1G1cSGWNqgb9i7Tmw0N+ebkAPY8xyrGGo8XFom0pT2iNQKghjzBn/ss5fiEgPrP9XnsSq+hqL17H+ct+EVTX3Xf/P84nIc/7ju7FKrDu9DNwEvOV/nA0s8/ckBLgvxnapNKbVR5VqB/wrk3oYY76f6Lao1KM9AqWSnIi8jrXKSecAlCe0R6CUUmlOJ4uVUirNaSBQSqk0p4FAKaXSnAYCpZRKcxoIlFIqzf0/EmPmEBJKptYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv = rsf.predict_survival_function(X_test_sel, return_array=True)\n",
    "\n",
    "for i, s in enumerate(surv):\n",
    "    plt.step(rsf.unique_times_, s, where=\"post\", label=str(i))\n",
    "plt.ylabel(\"Survival probability\")\n",
    "plt.xlabel(\"Time in days\")\n",
    "plt.legend()\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can also plot the predicted cumulative hazard function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7OklEQVR4nO3deXxV5bno8d+TBIgxiTEMYQglWBkVCxjAATG02iq19WK1Ra22dqD26Lm1p+2n0z1ae8+1np4OempbRevQVgVs69FDaRyORilVGRREUCCCSEDGEJIwZnjuH2vtsLKz9s5Ksld29t7P9/Phk73XlHe/xjx5p+cVVcUYY4yJlpXsAhhjjOmbLEAYY4zxZQHCGGOMLwsQxhhjfFmAMMYY4ysn2QVIpEGDBmlZWVm7Y4cOHeLkk09OToH6EKuHE6wuHFYPjkyvh9WrV+9T1cF+59IqQJSVlbFq1ap2x6qqqqioqEhOgfoQq4cTrC4cVg+OTK8HEdkW65x1MRljjPFlAcIYY4wvCxDGGGN8pdUYhB8RYevWrRw9ejTZRYkpNzeX0tJS+vXrl+yiGGNMm7QPECeffDIFBQWUlZUhIskuTgeqyv79+6mpqWH06NHJLo4xxrRJ+y6m7OxsBg4c2CeDAzgtnIEDB/bpFo4xJjOlfYAA+mxwiOjr5TPGZKaMCBDGGJOuli3exLLFm0J5tgWIXlBZWcm4ceM4/fTTufPOO5NdHGNMGtm3vZF92xtDeXZog9QiMhL4PTAUaAUWqOrdUdcIcDcwBzgMfFFVX3fPXeKeywYeUNWU/M3a0tLCTTfdxHPPPUdpaSnTpk3j05/+NBMnTkx20YwxUdYv28GmFbuTXYwu2VfTyKDS/FCeHWYLohn4lqpOAM4BbhKR6N+KlwJj3H/zgd8CiEg28Gv3/ETgap97U8KKFSs4/fTTOe200+jfvz/z5s3jqaeeSnaxjDE+Nq3Yzb6acP4aD8ug0nzGTi8J5dmhtSBU9QPgA/d1g4i8DYwANnguuxz4vTr7nr4qIkUiMgwoA6pVdQuAiCx0r/Xe22W3//d6Nuys78kjOpg4vJDbPnVGzPM7duxg5MiRbe9LS0t57bXXEloGY0ziDCrNZ+63pia7GIHtuuMOWAZc8IOEP7tX1kGISBkwBYj+zTgC2O55X+Me8zs+I8az5+O0PigpKaGqqqrd+cLCQhoaGgBoOt5ES0tLNz+Fv6bjTW3P93P48GGamk5cc+TIkXbvI44ePdqh7InU2NgY6vNTidWFw+rB4a2H+m0HyGpoYO2nfpTUMnVF/82bOT5mDO+E8N8y9AAhIvnAn4FbVDX6z3e/+Z0a53jHg6oLgAUA5eXlGp2V8Y033qCgoACAf/vM5C6UPDHGjBnDo48+2laG/fv3U1ZW1vY+Ijc3lylTpoRWjkzPWOlldeGwenB46+HxRx+j9fhxioqKklqmLpk2jaGXXcZHQvhvGWqAEJF+OMHhUVX9i88lNcBIz/tSYCfQP8bxlDNt2jQ2b97M1q1bGTFiBAsXLuSxxx5LdrGMMTFk5eUx6v7fJ7sYfUJog9TuDKXfAW+r6i9iXPY0cL04zgEOumMXK4ExIjJaRPoD89xrU05OTg733HMPn/jEJ5gwYQKf/exnOeOM2GMWxhjTV4TZgjgfuA5YJyJr3GM/AD4EoKr3AktxprhW40xzvcE91ywiNwPP4ExzfVBV14dY1lDNmTOHOXPmJLsYxhjTJWHOYvo7/mMJ3msUuCnGuaU4AcQYYxLmwKLF1C9ZwracsezIOY3m5mYef9Tp9q3XUyiUg0kuYd+R9tlcjTHGa13lJra3zqJuwGkAFDbXtJ0rlIOcPjFz96eOZgHCGJNRduScRmP/YoaPKWLs9BL2tmTZbK4YLBeTMSbjFLbWMvdbUznjghHJLkqfZgHCGGOMLwsQxhhjfFmA6AVf+tKXGDJkCGeeeWayi2KMMYFZgOgFX/ziF6msrEx2MYwxpkssQPSCWbNmUVxcnOxiGGNMl2TWNNe/fQ92rUvsM4dOgktTci8jY0waePHhBQDM/uL8hD87swKEMcakmT3btoT27MwKEPaXvjHGBGZjEMYYY3xZgOgFV199Neeeey4bN26ktLSU3/3ud8kukjHGdCqzupiS5PHHH092EYwxAbz5fCVvL69KdjG6ZO97WxlcNjqUZ1sLwhhjXG8vr2Lve1uTXYwuGVw2mgnnV4Ty7NBaECLyIHAZsEdVOywhFpHvANd6yjEBGKyqtSLyHtAAtADNqloeVjmNMcZrcNloPnebTWiBcFsQDwOXxDqpqv+hqpNVdTLwfeAlVa31XDLbPW/BwRhjkiDMHeVeFpGygJdfDVhHvTEmKSJjD2H256cicXb9DOnhToBY4tfF5LkmD6gBTo+0IERkK3AAUOA+VV0Q5/75wHyAkpKSsxcuXNjufGFhIWPGjOnhJwlfdXU1Bw+Gt9VhY2Mj+fn5oT0/lVhdODK1HrY/8j4AudMOULv5bVpaWjiy+wMA8oeXUjxmAoMnfiSZRexVs2fPXh2rp6YvzGL6FLA8qnvpfFXdKSJDgOdE5B1VfdnvZjd4LAAoLy/X6J2h3njjDQoKCsIpeQLl5uYyZcqU0J5fVVVlu2a5rC4cfb0e1i/bwaYVuxP+3IO6Bz2+geMv7QGcoFA68UwmnF/BWRfF7BXPSH0hQMwjqntJVXe6X/eIyJPAdMA3QPR127dv5/rrr2fXrl1kZWUxf/58vvGNbyS7WMb0GbECwc7NdQAMH1OU0O+nxzfQ0rK/LSjU5uT26UCZTEkNECJyCnAh8HnPsZOBLFVtcF9/HPhxkorYYzk5Ofz85z9n6tSpNDQ0cPbZZ3PxxRczceLEZBfNmD5h04rd7KtpZFBp++6uyJ7Rid4W9JF5tSDSNlOpqqoqoc9PJ2FOc30cqAAGiUgNcBvQD0BV73Uvmws8q6qHPLeWAE+KSKR8j6lqym6mMGzYMIYNGwZAQUEBEyZMYMeOHRYgjPEYVJrP3G9NTXYxTJQwZzFdHeCah3Gmw3qPbQFCGSH69xX/zju17yT0meOLx/Pd6d8NdO17773HG2+8wYwZMxJaBmNShV93kl/rwfQNtpK6lzQ2NvKZz3yGu+66i8LCwmQXx5ikiHQneQ0qzWfs9JIklcjE0xcGqXtN0L/0E62pqYnPfOYzXHvttVxxxRVJKYMxfUXQ7qQDixZTv2RJwr9/6+HDZOXlJfy56SijAkQyqCpf/vKXmTBhAv/yL/+S7OIY06uiu5SCdicdWLSYXbfdBkDetGkJLVNWXh45Awcm9JnpygJEyJYvX84f/vAHJk2axOTJkwG44447mDNnTnILZkxIvEEheqpq0YAjDHzrZbZdd1fcZxxeuRKAobffzqmf+2xCy5d7+/cS+rx0ZgEiZDNnziTM1erG9IauLFrzBoXoqarbrrueo++8A+PHx31G3rRpFF52WcKDg+kaCxDGmLjWL9tB1aMbgWCL1rxB4cCixdQvuIttbrKco++8Q+748Yz6w+9DLLFJFAsQxpi4Ii2HimvHdVi0FnMgeQ1sW3CiqygyjpA7fjyFl10WanlN4liAMCbDRLqL6upaObD69U6v31fTyJCiJvIXfL+tJRARHQCiWVdRarMAYUyGiaxFyAm4Nq1owBFOfe0JDn+wskMgsACQ3ixAGJOBigYcYeyaX1NUVNTptWHOKDJ9mwUIY9JUrJlHe7fWcvK+rfTfvBkCrDGwVkLmsgARsqNHjzJr1iyOHTtGc3MzV155Jbfffnuyi2UyQKwsqQXH9jB49yrqr72GCf/6r0kqnUkFFiBCNmDAAF544QXy8/Npampi5syZXHrppZxzzjnJLprJAH5pLbZddxeUwtYLLkhOoUzKsGR9IRORtm0dm5qaaGpqwk1lbowxfVpGtSB23XEHx95ObLrvARPGM/QHP4h7TUtLC2effTbV1dXcdNNNlu7b9IrmvXtp3r+/Q1qLyGI1YzoTWgtCRB4UkT0i8laM8xUiclBE1rj/bvWcu0RENopItYikfOKU7Oxs1qxZQ01NDStWrOCtt3yrxJiEat6/n9bDhzsct8VqJqgwWxAPA/cA8dbUL1PVdj+pIpIN/Bq4GKgBVorI06q6oacF6uwv/bAVFRVRUVFBZWUlZ555ZlLLYlJfZ/mR6rOKKcyDUffH+F/Qtto0nQitBaGqLwO13bh1OlCtqltU9TiwELg8oYXrRXv37qWurg6AI0eO8PzzzzPemvcmAfw23/EqbK1lRPOWXiyRSTfJHoM4V0TWAjuBb6vqemAEsN1zTQ0Qs9NeROYD8wFKSko6bEBeWFhIQ0NDgosdXHV1NTfeeCMtLS20trYyd+5cLrzwwg5lOnr0aKibpzc2Ntrm7K50qYu6ulZy8uHUs+t9z59atRAg5mdNl3roqsgfbJHPnqn1EEQyA8TrwChVbRSROcB/AWMAvyk+MfNlq+oCYAFAeXm5VlRUtDv/xhtvUFBQkKAid925557L2rVrO70uNzeXKVOmhFaOqqoqousmU6VyXXi7lZobnTUOFRX+u7Nt+92DAHwkxmdN5Xroid0vVQK0ffZMrYcgYgYIEYm7N6aq/qUn31hV6z2vl4rIb0RkEE6LYaTn0lKcFoYxGc+7+M32cjZhi9eC+JT7dQhwHvCC+342UAX0KECIyFBgt6qqiEzHGQ/ZD9QBY0RkNLADmAdc05PvZUw6CbqnszE9FTNAqOoNACKyBJioqh+474fhzDKKS0QeByqAQSJSA9wG9HOffS9wJfB1EWkGjgDz1Nl6rVlEbgaeAbKBB92xCWOMR8y9GFy23sH0VJAxiLJIcHDtBsZ2dpOqXt3J+XtwpsH6nVsKLA1QNmMyVv2SJXGDQKaud3jz+UreXl4V8/ze97YyuGx07xUohQUJEFUi8gzwOM5g8TzgxVBLZYwJJF237+zsl3w8NRuchailE/3XGg0uG82E8yu6WbLM0mmAUNWbRWQuMMs9tEBVnwy3WMYYP970GenSheQXDDr7JR9P6cQzmXB+BWdddEkiipfR4gYIEckC3lTVMwELCj3Q0tJCeXk5I0aMYEmcfmNj4mlLn5GVml1IQYOB/ZLvG+IGCFVtFZG1IvIhVX2/twqVju6++24mTJhAfb3/oiZjgsrKy4udPqOPiQ4IFgxSS5AxiGHAehFZARyKHFTVT4dWqjRTU1PDX//6V374wx/yi1/8ItnFMSZhOhsriA4IFgxSS5AAkTbbny1bvIl922PnrumOQSPzueCz8Sd13XLLLfz0pz9NasoPY7orXhDobKzAAkJqCzJI/VJvFCRdLVmyhCFDhnD22WdbvhfTJ3W1FeBlASC9dRogROQc4FfABKA/zuK1Q6paGHLZEq6zv/TDsHz5cp5++mmWLl3K0aNHqa+v5/Of/zx//OMfe70spu8LlMK7tTtJktt78/lKNv71v9j9UqW1AkxMQbqY7sFZ+/AEUA5cj5NUzwTwk5/8hJ/85CeAkxTsZz/7mQUHE5M315KfRKTwfvP5Sp6731mjWlRUZAHAxBQom6uqVotItqq2AA+JyD9CLpcxGcPbaogEh1i5lqK3D+2OSHfShy68mKv+6Rs9fp5JX0ECxGER6Q+sEZGfAh8AJ4dbrPRUUVFhaYUzUGfdRjs31wEwfEyRb4ZWb86lRC2OK514JoMnfqTHzzHpLUiAuA4n0+rNwDdxUnF/JsxCGZNOOus2GlLURMnulYxas8k5sAa2LThx/vDKlQDkTZvWpcVxsQafLRdRmvnb95yvl96Z8EcHCRCzgP9y92+4HUBELgOqE14aY9JIpOWwd2stBcf2MHVNpe91kQDAtGm+5/OmTaPwsss49XOf7dL3f3t5lW8wiOQi6vlQdwZa9RCs+1OyS9Hetr/DqJmhPDpIgPgV8C0RuVpV33aP/RiwfBHGxLHhr29Re0DJP/geg3evcra+8tHdABDE4LLRfO42/78sbdp1N6z7E+xaB0MnJbskJ4yaCZOuDOXRQQLEVuDLwJ9E5Eeq+gT+24IaYzya9+8n//BhZvb7B4U3hhMATBIMnQQ3/DXZpegVQQKEqurrInIh8LiIzMBZCxGXiDwIXAbscZP9RZ+/Fviu+7YR+LqqrnXPvQc0AC1As6qWB/kwxvQ1YeZNsn0PTNiyAlzzAYCq7gM+gbMnRJAcvA8D8SZWbwUuVNWzgP8LLIg6P1tVJ1twMMZfZIwhFtv3wPRUkFQbn/S8bgW+4/7r7L6XRaQsznnvWopXidlDm/rKysooKCggOzubnJwcVq1alewimRTlbTVEWgixxhhMN3Q2CN3Xxh9CFiTVxmCcrqCJQG7kuKp+NIHl+DLwN897BZ4VEQXuU9Xo1kXKefHFFxk0aFCyi2FSnHdmkrUQQtDZIPTQSaENCPdFQcYgHgUWAZ8EbgS+AOxNVAFEZDZOgPDO0zpfVXeKyBDgORF5R1VfjnH/fGA+QElJSYeZGYWFhUnPoqqqNDY2MmDAgJjXHD16NNRZJY2NjTZrxdVbddHc3AwkdrZQXV0d/YpOpeRCp/e2tgfPt58Jh7ceJtfVQe5I1oyO00nSCGRIvQUJEANV9Xci8g03s+tLIpKQDK8ichbwAHCpqu6PHFfVne7XPSLyJDAd8A0QbutiAUB5eblGr1R+4403KCgoAODFhxewZ1vP8thEGzLqNGZ/cX7ca7KysrjiiisQEb72ta8xf37H63Nzc5kyZUpCy+ZVVVVlq7hdvVUXjz/6GEBCv9fulyoT9kz7mXC0q4etRUBi/5ulsiABosn9+oGIfBLYSQLGC0TkQ8BfgOtUdZPn+MlAlqo2uK8/jrPuImUtX76c4cOHs2fPHi6++GLGjx/PrFmzOr/RpAy/dBqJyrxqTLIECRD/JiKnAN/CWTRXiJNyIy4ReRyoAAaJSA1wG9APQFXvBW4FBgK/ERE4MZ21BHjSPZYDPKaq/ktQu6izv/TDMnz4cACGDBnC3LlzWbFihQWINOOXTiMRmVeNSaYgs5giK6YPArODPlhVr+7k/FeAr/gc3wKkTRaxQ4cO0draSkFBAYcOHeLZZ5/l1ltvTXaxTAgiWVgjyfUSlVjPmGQJOovpq0CZ93pV/VJ4xUofu3fvZu7cuYAzaHnNNddwySWWdz/dNO/dS/P+/Wy77q52yfWCJtYzpi8K0sX0FLAMeB5nZbPpgtNOO421a9cmuxgmBN5xh9oDSv7hw5AVbm4lY3pTkACRp6rf7fwyYzKLd9yhsLWWEf13MeqhcNJqGJMMQQLEEhGZo6pLQy+NMX1MrM1+mvfupfaAUthay9Q1laGNN0TnW7L8SqY3xQwQItKAs6JZgB+IyDGcKa+Ck8CvsHeK2HOqijsrqk9S1WQXIaPF2/HNu9ubVyRT64j+uwC6tJFPtHhJ92o2vAU4O8CB5VcyvStmgFDVgt4sSFhaWlrYv38/AwcO7JNBQlXZv38/ubm5nV9sQhFvx7cOu725Ii2GIF1KnWVdjQ4CXqUTz2TC+RWcdZFNbDC9L0gXU0o7dOgQDQ0N7N2bsOwgCZebm0tpadrmKuzzmvfuJb9xv++Ob7F2ewvSYogEhngBIHLcgoDpi9I+QKgqo0dbn63xt37ZDvbU9aPInYEUrSczkiKJ9SwAmFSV9gHCZCbvuEJkjUJEc3NzW56k2uyhAIwMaQaSpeM2qSxQgBCRmcAYVX3IXTiXr6qxdyoxJski+0EXttbS2lAPQFZBx3kVxS27GNG8hUmXjO3tIhrT5wVZSX0bUA6MAx7Cyaf0R+D8cItmTPdFZhmdl/Uy9KNdN5FlMTUmmCAtiLnAFOB1cFJxi0hazHAy6S3M/aCj+c1UsjULJtUF2ZP6uDoT9RXa0nEbYzz89oe2NQsm1QVpQSwWkfuAIhH5KvAl4P5wi2VM6rEBaZNugqT7/pmIXAzU44xD3Kqqz4VeMmOMMUkVZJD6m8ATFhSMac877mDjDSlk1UOw7k9tbyfX1bVtNcqudTB0UlKK1RcF6WIqBJ4RkVpgIfAnVfVPXOMhIg8ClwF7VLXDElJx8l7cDcwBDgNfVNXX3XOXuOeygQdU1drtJpDI+oewt/t88/lKnrv/HsBZCW3jDSlk3Z9iB4Khk2DSlTFvfey193lqzY4QC9c9E4cXctunzkj4c4N0Md0O3C4iZwGfA14SkRpVvaiTWx8G7gFiTSO5FBjj/psB/BaYISLZwK+Bi4EaYKWIPK2qGwJ8HpPhIusf8htr2hLp9USsPEqR9BkXf/VmWyHd10S1EDqIBIcb/grAmi5Me3543eN8oP8gb0DfWmN8/NCHgZ8m/Lld+ZR7gF3AfmBIZxer6ssiUhbnksuB37szpF4VkSIRGYazc121u/UoIrLQvdYChOlUZP3DzKyXKbyka9lV/YJBrDxKlj6jD4vXQoBOWwnxHMxeQVbOB0wcmvi/1ntifPHAUJ4bZAzi6zgth8HAn4CvJuiv+RHAds/7GveY3/EZcco3H5gPUFJSQlVVVbvzjY2NHY5lokyph+bmZujfn61fcHfE9fnM0XWxd8Naaje/TePOGgDyh59InJg/vJTiMRMYPLHjNum1kNJ12pd/JobtfIaS3S936978xq005o9mzejvxL6okbafja7UQ3NzMzkM4wu5X+hW2UJzOJyfxSAtiFHALaq6JsHf2y/3tsY57ktVFwALAMrLyzW6qWirZh2ZUg+RHEvxPmt0XSx6qZKmugMZ1yroMz8Tfl1C2/7ufB01s+vPK5pC0aQrqSivCHR5V+ohZ+uvgPg/X+kk3oZBhapaj9uxJSLF3vOq2tMRwBpgpOd9KbAT6B/juDFA/A1+ujs4bWsYksivS2jUTKcbqPyGXi1KZ4PQh7W5z40/hCneJ30MZxbSajr+Za/AaT383k8DN7tjDDOAg6r6gYjsBcaIyGhgBzAPuKaH38ukifXLdlD16Eag4y5vgLM3dPOWwM978/lKaja8FXOvBtNLPIPGyfTUmh1s+KCeicP8N8zMG5DDoJMH9HKpkifejnKXuV+7NblbRB4HKoBBIlID3IaT6A9VvRdYijPFtRpnmusN7rlmEbkZeAZnmuuDqrq+O2Uw6SfScqi4dhxnXDCiw/lt193V4Vj04HNdXR27X3I2B4oMQtsU1QTrbCaRVx9bezBxWCGLvnau77kbKlNmp+WECDJI/T+q+rHOjkVT1as7Oa/ATTHOLcUJIMYAJ7qV9tU0MqSoifwF32fbgo7XRbYC9QYF29IzAbryCx+6NobQg1lFJlzxxiBygTycFsCpnOhiKgSG90LZTIbyG2PYubkOcPaIPvW1Jzj8wUryorYBhRNbgb7qJs8bXDa6QxDoM4OzqaSzqaPRkjSGYBIrXgvia8AtOMFgNScCRD3OQjZjQhFpKQwqzW87NnxMEWOnl5C/4Psc/mAlQ2+/Pf42oLe/bgPPidZHxglM74k3BnE3cLeI/LOq/qoXy2QyTHSLIRIc5n5ratuxA4sWU7/gLo6+8w5506Z1a49o002rHnK6jLoz5dSktCCpNn4lImcCE4Fcz/He2YnFpD3v9qAA+cDAt7a0G3A+vHIlAHnTplF4WddWSJseiow9ZMg4wYHsl7mh0meAC9hYu5FxxeN6uUTJE3TL0QqcALEUJ4fS34mdY8mYLmm3PWgMkcDQWcshMjht2VUTbNTMlBlP6GpCvbq6I/x24ysAbGh8Fgb9iQ92Q3lJeYdrxxWPY85pcxJW1r4uyIqPK4GPAG+o6g0iUgI8EG6xTKZJ1Pag3uBgU1eJO/uoXZrreEKehproDKmvbXVaojNGF3dyZUd5xW9yGLj13Fu5auxVCStTqgoSII6oaquINItIIU7Svp4ukjMZ7sCixdQvWQJAa+sssvLyEvZsG5z26OrsIz8hT0PtbHFaV80YXczlk0dwzYwPBbremdXmrHtw1jmUW3BwBQkQq0SkCGeb0dU4aa5WhFkok/7qlyxpW7OQlZdHzsBwslEaYs4+6kqa667qSqsgEhxiLU4zyRNkkPqf3Jf3ikglUKiqb4ZbLJPODixazOaa/uydfAu548fTWNPIoMH5nd9okqorv/S70s0zcVghl0/uuCreJF+8hXJT452L7P5mTFCRbqXNNf3ZOM5JrzUcGFSaz9jpJcktXDpK0PTUSGDoyi/9rnbzJNsTm55g6ZalGTdLqTPxWhA/j3NOgY8muCwmza2r3MT21lnUjXOGsGLlUzIJkqDpqZExglT7pR9LJBhE1NXVUb2tGnBmLmXSLKXOxFsoN7s3C2LS24FFi9l+ZBCNp5S2rYrubnCItQ0okFnTW4NsrTlqJo+1fIyn7nulw2nv9M54Um2MIDoARFu1exXQfhprJDDY4HR7QdZBXO933BbKma5YV7mJuqLzGFLU1G6FdHfEW+eQdtNb4wWBzhLiDZ3Ea/kf5QdPrgO6N+0TUm+MoLOuouhgYLm5Ygsyi8mbES0X+BjwOrZQznTBjhynW2niJ3u274J3/4aMmMoab5pqgIR4v7jvFaCWO+ZO6tA15J3emW7GFY/joUseSnYxUl6QWUz/7H0vIqcAfwitRCZtFbfs4owLgg1dxepGSuv9G/xaC5Hg0IMkeTNGF6fcuEFn3UTx2EBz4nRn77zDwJggF4rIJcDdOBv/PKCqd0ad/w5wracsE4DBqlorIu8BDUAL0KyqHde9m5RwYNFiWhvqySrwXwjlFwxi7eGQlvs3RAKDX5dRFxap+U1DTeQCtK7oyS948B8nCCrT0mGEKcgYxH/jzFoCyMLJybQ4wH3ZOGnBL8bZf3qliDytqhsi16jqfwD/4V7/KeCbUXtdz1bVfQE/i+mDDixazOp7n6Vu3DUMKWryvcZvTCEtAwH4txK8gaELeyhEBwS/aajJGj/o6ZRRGzTuG4K0IH7med0MbFPVmgD3TQeqVXULgLv39OXAhhjXXw08HuC5JoWsq9zUtuYh3vhDWqfH8AYFv1ZCNwNDdEBIxjTUWC2FSHCwcYDUJs7OnwEudPIwtQWUqL/0/a6/ErhEVb/ivr8OmKGqN/tcm4fTyjg98lwR2QocwGm93Keqvvl3RWQ+MB+gpKTk7IULF7Y739jYSH6+rdJNRj2ctGwZ1esHU1c0lmHlQvHp0nZu74a11G5+G4Aj+/Zy0qDBjLt8Xq+Uq7fqYtjOZyjZ/TJFB53usrpTnAC5u2QWHwz/RKBnVG1v4pWdze2ObTzQCsC4U7M4d3gOFSP7dat8XamH5Q3LWXVoVYfj1cec9QOnDzi9w7nyk8s5v+D8bpWtN2X674jZs2evjtWFH6SLaT7wf4EjQCvOznJK5wn7xOdYrGj0KWB5VNA5X1V3isgQ4DkReUdVO+SDdgPHAoDy8nKNnq5mU9gcyaiHbb97kGoGM6SoiSu+8omY+0QXFRU53Um9VL7Q68JvTGHSlRS5LYQiIEjHy2Ovvc/D6ztOUZ1RREJaCpF6CDJesKrWf0ygnNTvCrLfEbEF6WL6DnBGN8YCaoCRnvelwM4Y184jqntJVXe6X/eIyJM4XVaxNwwwSePNzOpVva+QulFjGT64iDefr+S5++8BnKCQtmMMcGJqajf2ZfaOK0S6kPymqAYV75d/XV0dj1Q+EmhA2MYEMlOQAPEuzsylrloJjBGR0cAOnCBwTfRF7rTZC4HPe46dDGSpaoP7+uPAj7tRBhOC6IDg3e0NYFvOWHbknEbtqKEAjJ1ewlsvOPMaLv7qzekXFKIHnrs5NfWx195vt6gtEWMKQQaL7Ze/iSVIgPg+8A8ReQ04Fjmoqv873k2q2iwiNwPP4ExzfVBV14vIje75e91L5wLPquohz+0lwJMiEinjY6paGfAzmZB5U3VD+93e1i/bwbpHNwK0S6nx1gtOyyHtggN0XMzWzf0TIi2HnrQY/MQaLLauFdOZIAHiPuAFYB3OGERgqroUZ5tS77F7o94/DDwcdWwLzi52po/KHT+eUX9wFtOvX7aD11fshp+/zs7NdUD7RHze1c9pJdJy6OFitki3UiQhXqotajPpK0iAaFbVfwm9JKbPi3QtVe8rZE/pebz+cyfjeyQoDB9T1CERn3fsIa1WP696CJbc4ryOjDUE4LeYzTtdtSdrFvzGG2xVsemJIAHiRXcm03/Tvosp7jRXk34iXUt7Jt9Cw4Ah5LrH42VnjcxaSumxh3iL2y67q8MgdLyNdfwWsyVq/YLfeIOtKjY9ESRARAaWv+85FmSaq0kjBxYt5vDKleRNm0bu+PHkQuCsrCk99hDdUoiIM0Mp3h7LiQgGtjjN9JYgyfoyJLm+iad+yRJ2DDufAyOuoq6mkUGlnS8sSouxh0jLwaelEE+i9k/wCwaxpqVaa8Ekmu0HYQLbe/pHaTx2UuAtQiPdSyk/9jBqZpeCQyL5dRvZtFTTW2w/CBPItpyx1GYPZXhpfqddS5EV03vf25q87qU4G+1MrquDrUXBnhNrLwYf3tlI3cmgGm+Q2bqNTDLYfhAmJu+CuO2ts2AAnbYcoldMJ631EG+jna4IsKbBL3led2Yj2SCz6WtC3Q/CpDbvgrisvDyGFDVxxgUj4u4JHcmx1CdmLcVYm7AmwQvEvGsYujMAHWk5WGvB9DWh7QdhUtuBRYvZXNOfvZNvIXf8eBprGhk02BmYjrcndFrnWPLwTmWNdCkFGZTubNDZWgumLwlzPwiTglb+8imqNxyitaGeOncfh+HAoNJ8TsrbyKLbF7cFh7TdvyFKZ4vburIpjw06m1QSM0CIyOlAiaq+FHX8AhEZoKrvhl4606tW/vIpVmwsgOwCigtgSFETEz95ZtsCOG9w6LMzk6LTX/RQdAK9iJ6sZ7BuJJMq4rUg7gJ+4HP8iHvuUyGUxyRR9YZDkF3A9HENTPtmh8S7QB/Z+S3ODKUO23fSsQVQV3eE3258JdC36k7K7Xgpti31hUkl8QJEmaq+GX1QVVeJSFl4RTLJVNyyyzc4JHXRW3RA8Nu2M8INDI+1fIynVu+A1a/4prcIqisthUhgiLe/gs1KMqkkXoDIjXPupEQXxCTX+mU7qM0eSnHLLt/zSV30Ft1lFGAjnqfue6Vt8Dj6l7yT5rrnq5wj/AKDjSmYdBAvQKwUka+q6v3egyLyZWB1uMUyvSUyKF2b7WzuM6J5S8xrk7LobdVDToth1Mwup9NOVLqLzkQGni0wmHQTL0DcgrNpz7WcCAjlQH+cTX46JSKXAHfjbBj0gKreGXW+AngK2Ooe+ouq/jjIvabn2g1Kt+xiRPMWJl0ytu28d71DrGmtofB2KUW6k7qxAU9vsoFnk45iBghV3Q2cJyKzgUjH819V9YUgDxaRbODXwMU4+1OvFJGnVXVD1KXLVPWybt5reqCzQWnveodem7kUnT21i/s69zTdhZ94g85gA88mfQVJtfEi8GI3nj0dqHZ3h0NEFgKXA0F+yffkXhPAgUWLaW2op7iADsHBm0upV2Yt+bUYYmRPjbfXAiRu8x2vzvZ1toFnk65EVTu/qjsPFrkSuERVv+K+vw6Yoao3e66pAP6M00rYCXzb3be603s9z5gPzAcoKSk5e+HChe3ONzY2kp/feWrqdBddD8d++yzVp15Efv+DjLriVAD2blhL7ea3adzprIPMH15K8ZgJDJ7Y891fh+18hpLdL/ueKzropOeoO8VpqO4umcUHwz/R7pqq7U28srOZjQecXW/HnZoV83udOzyHipH9Yp7vys/E8oblLKxdyOkDTucbQ78R6J5UYf9vODK9HmbPnr1aVTtOuaN7uZiCEp9j0dHodWCUqjaKyBzgv3DyPAW51zmougBYAFBeXq7ROXZsY3ZHdD08/uhOAMqvmn5iIdxLlTTVHUhsuoxI6yDe1NQipxupyG0xFAHev9Ufe+19Hl5/YrFakGmn8bqF6hrrKMopClb8Wmdm0jVTr6FibEWge1KF/b/hsHqILcwAUQOM9LwvxWkltFHVes/rpSLyGxEZFORe033e7qUzLvhou3MJ71KKTFHt4lhChHclc1cWq3XWLRSUzUwymSzMALESGCMio4EdwDxObF8KgIgMBXarqorIdJxkgPuBus7uNV1z0rJlbPvdgwBsrulP3bhrGFLUlJiZSvFWNkfWL3Rhiqp3nKE7K5kjYs0ssr8YjQkmtAChqs0icjPwDM5U1Qfd8YUb3fP3AlcCXxeRZpwUHvPUGRTxvTessqa7A4sW0/DCNqpLziOroJDacc6ah4JRtTx3/6OAs8Yh8EylrqxsDrCfgld07qMgXUrxNtoxxnRfmC0IVHUpsDTq2L2e1/cA9wS913RP/ZIl7C45j0ODTmfw6GKG42z889YLzwKd7N3g1zqIDgjdnIrqJ1aLId6Ygl9qC5tZZEzPhRogTN+wLWcsdSePZfjoonbbhb71QpzV0fEGl7s5nhARb51CrBZDvDEFGycwJhwWINLcgUWL2X5kUIftQn2T7/mtR+hhMIjo7gY74LQeVu1eRXlJua1WNqYXWYBIQ969pCMD0vn9D7ZNZ/XuG91uzMGbFC8BgcFvsLmrG+wAbV1L1mVkTO+yAJFmDixazOp7n2V31IC0Fr/Potu/B/jsGx29yU43ZxxF8waFnmywA043knUhGdO7LECkkUhw2BjZKnRMEcOBk/I2sr7qOcAZc2i3EC4691GMGUexAkG8vRZ6GhSgffeSMaZ3WYBII+sqN7UFh4prx7XbKhR8Zit5g0MnuY9iBYJEBIF4rHvJmOSxAJFGduScBrQPDpHB6Pzhpd0KDt41CWEGgojo6ayRfRase8mY3mcBIs0Ut+yi5dhxFt3+K+DEeEPxmAntL4zMVvIEh+hupJ6sYu6u6Omstp7BmOSxAJFGDh2v5kjTezx3/x6g/XhDbY7PDrKjZrYLDt7WQuRrb7QaotnmO8b0DRYg0kBk29CGphq0xT8ba1VVle+90WMMvdFaiLcq2lJkGNN3WIBIcd5tQ3No5qSCQZ1nY/Xs8xxZ1dwbrYVIYPBLjRFhXUrG9B0WIFKcd9vQLfVDYl43bOcz8NB/OG/cVdIL6qay4WDXVjV3lbe14A0MlhrDmL7PAkSKinQr7T+2A21+ni31Q+Km6y7Z/TIc3d62SnpB3VR+dXBml1c1R+tsv2ZvULDAYExqsQCRoqo3HKJeT0Gbn6eFOmBI7HTdqx5ytvUcNbNtlfT/3PcKE/Poccuhs415LCgYk7osQKSoQ8erOd70HjLgEMPKxsYed/Cud+jCvgwRnbUQIsHBZh0Zk35CDRAicglwN86mPw+o6p1R568Fvuu+bQS+rqpr3XPvAQ1AC9Aca1PtTLTyl09R37wTWg4womxc/E1+3PUOd/X7Cq+sHg+rXwGImW47WmctBBtUNiZ9hRYgRCQb+DVwMc4e0ytF5GlV3eC5bCtwoaoeEJFLgQXADM/52aq6L6wypppIltZNTedBNhTGm7HkJuA7vmMtr7dO4K6GjzJj0InTXRl7sBaCMZkpzBbEdKBaVbcAiMhC4HKgLUCo6j88178KlIZYnpS3rnIT21tncSC3Dj1SQ8HIM/0v9HQrvd46gadazuOLZ/TnR9d1Pt7gl+rC1iUYk5nCDBAjgO2e9zW0bx1E+zLwN897BZ4VEQXuU9UFiS9iaojMWNpDPi2sQo/UAHTsWoraBe525rNhxBVcPnkEw49s6fT7PLHpCX78yo+BE2sUrAvJmMwlqhrOg0WuAj6hql9x318HTFfVf/a5djbwG2Cmqu53jw1X1Z0iMgR4DvhnVX3Z5975wHyAkpKSsxcuXNjufGNjI/n5+Yn9cL3kpGXLqHvrCNWnXgRAy8E/0MIB8oYOpXjMBAZP/Ei76ye/8UPyG7fSmD+ahw+fw99yLuL7M04CYtfD8oblrDrkTEWtPlYNwLzieZxfcH6YHy2pUvlnIpGsHhyZXg+zZ89eHWuMN8wWRA0w0vO+FNgZfZGInAU8AFwaCQ4AqrrT/bpHRJ7E6bLqECDclsUCgPLycq2oqGh3vqqqiuhjqeLlR/7eFhychXAlQEnscYetRVA0haUTf8tdT65jxugiKiqcbqXoemhb1VzrWadAZkxJTeWfiUSyenBYPcQWZoBYCYwRkdHADmAecI33AhH5EPAX4DpV3eQ5fjKQpaoN7uuPAz8Osax9UiR994c/spste9bGXQgHsLvhKPsaj/GDjU7SvXiD0JHZSbZOwRgTS2gBQlWbReRm4Bmcaa4Pqup6EbnRPX8vcCswEPiNiMCJ6awlwJPusRzgMVWtDKusfdWh49W0Nj3P+qoT2VljTWl97YmfM6N2FVtbJwTOq2Szk4wx8YS6DkJVlwJLo47d63n9FeArPvdtAT4SfTxTRKazHj4uNFPvm50VnEysjf+4n/OPvMiM406rQSZdxaKrOs5WWt6wnEcqH2l7b7OTjDGdsZXUfcyBRYvZddttAMjUc+nfbwgtl/wT/2/NDnj3lbbrPnZ4KWcdeI5zst4GYH3/STSOmcuMq77V7nl+Yw1gs5OMMZ2zAJFkkdZCxOGVKwFo/NqdNK36G5KVxbqn7+Kb2f+gIPfEf64zjq+DLNhdXE7JeZ/nDJ8tQ+HEWMPpA07nmqnX2FiDMSYwCxBJ4A0KkYCQN21a29c9Uz/j7PEAoI38pN+jzuthMz1PmQmTrqQkRmCItBwiXUlfyP0CFWMrwvg4xpg0ZQGil3m7kPKmTWPX1EnsPLWAnCGDAThUd5y6Vc56weys/RRkNTg3evaO7kz0grc5p83xmWBsjDHxWYDoZfVLlvDm6CnsHVKA9OvH0YZa2HeY3GPOgrZjR5oBKDqliX5NjYzN28H6/pNidiF5Re/Yduu5t7Z1KVXtrArnAxlj0pYFiF62klOoKdwKR+vJ7VdGbkEZ+QPPonCIM3hcv3c7w7Ne4OL+jwPwausEGsfMjftMv608bW2DMaanLED0ogOLFrPz+EEAzqi4lku+fnWHa9bf8b8Zefxd1vefxPKTZrN5Yin79FXurYzdgrDAYIwJgwWIXhQZmB6QO6ItODz22vs8tWZH2zXfPt7C9v4f5owf/J0Nm56gMip5nh8LDMaYMFiACFH0FNZXDwpNJ+0nt19Z27HGf9zPtw8+T17/bADWFOxiadGp9K+8wXcswRhjeosFiBBEAsM7777DzlMLyCpwpqzuO8kZgD5lcDHr73CmrM53V0BHprD+TI6yNRvGYS0DY0xyWYBIAL/Fbu8XF/DWyCEA5OY7+yBlHW3mWP4Edhb/nmcLDnE0KxcYTctJgygc7Fy7sfaA5UgyxvQJFiASoH7JEo6+8w6548cDsH7aRWw7vhWAnLyLqCua7Fx3pImC3PW8W1DPptyTmDB0SodnWQoMY0xfYQGiByIth1cPCnvGnEnWgDwAjjY4wWHA0Eupy2niorzvUZCbw/N5x1mTe4CN/fsz4eQR1kowxvRpFiC6yNud9Oq+ZnYX5tB00n5ogVzKAMgtKCOrbArvH1nP1OKF3HfyyZB7CqvkGJBLeV4pcz7ypeR9CGOMCcACRBdFupN2n/lpdhTXoC17yc13AkLV4Ilt12Uf/CUDBr3Bj08aCEB5ySTKwQadjTEpwwJEHNGDzwDV+wrZM/kW9jTvQZtrOHX4WHLn3cKSl75Ns/4n2VkCwLqSFrytBQsKxphUE2qAEJFLgLtxdpR7QFXvjDov7vk5wGHgi6r6epB7w7Lyl09RveEQAI316zgCSPaJamoqbobDL6PNNQCsLnyNXW98ik0lrQCU6wD3aw5zhl/AVR//ZW8U2xhjEi60ACEi2cCvgYuBGmCliDytqhs8l10KjHH/zQB+C8wIeG8o1q5ZS0NTDTk0cxxnq89+WYPazudwlP7ZjezJP8zG4YdpGnaUbISzmvvxv0ZeaAHBGJM2wmxBTAeq3e1DEZGFwOWA95f85cDvVVWBV0WkSESGAWUB7k2Ye679Ms0tziK2Ft0PwL6iY2TTQl1JPTWlNR3u2davhdKcEv583athFMkYY5IuzAAxAtjueV+D00ro7JoRAe8FQETmA/MBSkpKqKqqane+sbGxw7Fo6nmdLQNpOamZDVPfAuBYvyKy84o73DMcmHpyeafP7iuC1EOmsLpwWD04rB5iCzNAiM8xDXhNkHudg6oLgAUA5eXlWlFR0e58VVUV0ceidXY+HQSph0xhdeGwenBYPcQWZoCoAUZ63pfScV+zWNf0D3CvMcaYEGWF+OyVwBgRGS0i/YF5wNNR1zwNXC+Oc4CDqvpBwHuNMcaEKLQWhKo2i8jNwDM4U1UfVNX1InKje/5eYCnOFNdqnGmuN8S7N6yyGmOM6SjUdRCquhQnCHiP3et5rcBNQe81xhjTe8LsYjLGGJPCLEAYY4zxZQHCGGOMLwsQxhhjfIkzTpweRGQvsC3q8CBgXxKK09dYPZxgdeGwenBkej2MUtXBfifSKkD4EZFVqlqe7HIkm9XDCVYXDqsHh9VDbNbFZIwxxpcFCGOMMb4yIUAsSHYB+girhxOsLhxWDw6rhxjSfgzCGGNM92RCC8IYY0w3WIAwxhjjK60DhIhcIiIbRaRaRL6X7PKETUTeE5F1IrJGRFa5x4pF5DkR2ex+PdVz/ffdutkoIp9IXsl7RkQeFJE9IvKW51iXP7eInO3WX7WI/KeI+G1c1WfFqIcficgO92dijYjM8ZxL13oYKSIvisjbIrJeRL7hHs+4n4keU9W0/IeTJvxd4DScDYjWAhOTXa6QP/N7wKCoYz8Fvue+/h7w7+7riW6dDABGu3WVnezP0M3PPQuYCrzVk88NrADOxdnR8G/Apcn+bAmohx8B3/a5Np3rYRgw1X1dAGxyP2/G/Uz09F86tyCmA9WqukVVjwMLgcuTXKZkuBx4xH39CPC/PMcXquoxVd2KsyfH9N4vXs+p6stAbdThLn1uERkGFKrqK+r8Zvi9556UEKMeYknnevhAVV93XzcAb+Psc59xPxM9lc4BYgSw3fO+xj2WzhR4VkRWi8h891iJOrv04X4d4h5P9/rp6uce4b6OPp4ObhaRN90uqEi3SkbUg4iUAVOA17CfiS5L5wDh11eY7nN6z1fVqcClwE0iMivOtZlYPxD7c6drffwW+DAwGfgA+Ll7PO3rQUTygT8Dt6hqfbxLfY6lVV10VzoHiBpgpOd9KbAzSWXpFaq60/26B3gSp8tot9tUxv26x7083eunq5+7xn0dfTylqepuVW1R1Vbgfk50I6Z1PYhIP5zg8Kiq/sU9bD8TXZTOAWIlMEZERotIf2Ae8HSSyxQaETlZRAoir4GPA2/hfOYvuJd9AXjKff00ME9EBojIaGAMzoBcuujS53a7HBpE5Bx3psr1nntSVuQXomsuzs8EpHE9uOX+HfC2qv7Cc8p+Jroq2aPkYf4D5uDMYHgX+GGyyxPyZz0NZybGWmB95PMCA4H/ATa7X4s99/zQrZuNpPDsDOBxnO6TJpy/+r7cnc8NlOP8An0XuAc300Cq/ItRD38A1gFv4vwiHJYB9TATpyvoTWCN+29OJv5M9PSfpdowxhjjK527mIwxxvSABQhjjDG+LEAYY4zxZQHCGGOMLwsQxhhjfFmAMBlBRAZ6Mpru8mQ4bRSR34Tw/W4Uket7cP/DInJlIstkTFflJLsAxvQGVd2Pk24CEfkR0KiqPwvx+90b1rON6S3WgjAZTUQqRGSJ+/pHIvKIiDwrzt4aV4jIT939ACrd9A2RPQJecpMiPhO1WhnPs77tvq4SkX8XkRUisklELvC5XkTkHhHZICJ/5UQiOUTkVhFZKSJvicgC99oPi8jrnmvGiMhq9/Wd7nPeFJHQgqBJfxYgjGnvw8AncVJA/xF4UVUnAUeAT7pB4lfAlap6NvAg8P8CPDdHVacDtwC3+ZyfC4wDJgFfBc7znLtHVaep6pnAScBlqvoucFBEJrvX3AA8LCLF7rPOUNWzgH8L/MmNiWIBwpj2/qaqTTjpKbKBSvf4OqAM55f4mcBzIrIG+D+0T+gWSyRh3Gr3OdFmAY+rk1hvJ/CC59xsEXlNRNYBHwXOcI8/ANwgItnA54DHgHrgKPCAiFwBHA5QNmN82RiEMe0dA1DVVhFp0hO5aFpx/n8RYL2qntud5wItxP7/rkPeGxHJBX4DlKvqdnf8JNc9/Wec1sgLwGp3nAURmQ58DCdB5c04QcWYLrMWhDFdsxEYLCLngpNWWkTO6OSeIF7GySia7Y5pzHaPR4LBPnd/g7aZTap6FHgGZ8+Hh9zy5AOnqOpSnO6syQkom8lQ1oIwpgtU9bg7/fQ/ReQUnP+H7sLJoNsTT+L8pb8OJwPxS+73qxOR+93j7+Gksfd6FLgCeNZ9XwA85bY8BPhmD8tlMphlczUmhbkzpU5R1X9NdllM+rEWhDEpSkSexJl1ZWMMJhTWgjDGGOPLBqmNMcb4sgBhjDHGlwUIY4wxvixAGGOM8WUBwhhjjK//D3Vu5+0sH6EzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv = rsf.predict_cumulative_hazard_function(X_test_sel, return_array=True)\n",
    "\n",
    "for i, s in enumerate(surv):\n",
    "    plt.step(rsf.unique_times_, s, where=\"post\", label=str(i))\n",
    "plt.ylabel(\"Cumulative hazard\")\n",
    "plt.xlabel(\"Time in days\")\n",
    "plt.legend()\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutation-based Feature Importance\n",
    "\n",
    "The implementation is based on scikit-learn's Random Forest implementation and inherits many\n",
    "features, such as building trees in parallel. What's currently missing is feature importances\n",
    "via the `feature_importance_` attribute.\n",
    "This is due to the way scikit-learn's implementation computes importances. It relies on\n",
    "a measure of *impurity* for each child node, and defines importance as the amount of\n",
    "decrease in impurity due to a split. For traditional regression, impurity would be measured by the variance, but for survival analysis\n",
    "there is no per-node impurity measure due to censoring. Instead, one could use the\n",
    "magnitude of the log-rank test statistic as an importance measure, but scikit-learn's\n",
    "implementation doesn't seem to allow this.\n",
    "\n",
    "Fortunately, this is not a big concern though, as scikit-learn's definition\n",
    "of feature importance is non-standard and differs from what Leo Breiman\n",
    "[proposed in the original Random Forest paper](https://github.com/scikit-learn/scikit-learn/pull/8027#issuecomment-327595859).\n",
    "Instead, we can use permutation to estimate feature importance, which is\n",
    "preferred over scikit-learn's definition. This is implemented in the\n",
    "[permutation_importance](https://scikit-learn.org/stable/modules/generated/sklearn.inspection.permutation_importance.html#sklearn.inspection.permutation_importance)\n",
    "function of scikit-learn,\n",
    "which is fully compatible with scikit-survival."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.inspection import permutation_importance\n",
    "\n",
    "result = permutation_importance(rsf, X_test, y_test, n_repeats=15, random_state=random_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>importances_mean</th>\n",
       "      <th>importances_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>pnodes</th>\n",
       "      <td>0.076616</td>\n",
       "      <td>0.019106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>0.016562</td>\n",
       "      <td>0.008774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>progrec</th>\n",
       "      <td>0.011513</td>\n",
       "      <td>0.013504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>horTh=yes</th>\n",
       "      <td>0.008220</td>\n",
       "      <td>0.004313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tgrade</th>\n",
       "      <td>0.004831</td>\n",
       "      <td>0.003773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tsize</th>\n",
       "      <td>0.002542</td>\n",
       "      <td>0.008665</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>menostat=Post</th>\n",
       "      <td>-0.000087</td>\n",
       "      <td>0.000978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>estrec</th>\n",
       "      <td>-0.007206</td>\n",
       "      <td>0.009461</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               importances_mean  importances_std\n",
       "pnodes                 0.076616         0.019106\n",
       "age                    0.016562         0.008774\n",
       "progrec                0.011513         0.013504\n",
       "horTh=yes              0.008220         0.004313\n",
       "tgrade                 0.004831         0.003773\n",
       "tsize                  0.002542         0.008665\n",
       "menostat=Post         -0.000087         0.000978\n",
       "estrec                -0.007206         0.009461"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(\n",
    "    {\n",
    "        k: result[k]\n",
    "        for k in (\n",
    "            \"importances_mean\",\n",
    "            \"importances_std\",\n",
    "        )\n",
    "    },\n",
    "    index=X_test.columns,\n",
    ").sort_values(by=\"importances_mean\", ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result shows that the number of positive lymph nodes (`pnodes`) is by far the most important\n",
    "feature. If its relationship to survival time is removed (by random shuffling),\n",
    "the concordance index on the test data drops on average by 0.076616 points.\n",
    "Again, this agrees with the results from the original\n",
    "[Random Survival Forests paper](https://projecteuclid.org/euclid.aoas/1223908043)."
   ]
  }
 ],
 "metadata": {
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}