psteinb/random-survival-forest-shap.ipynb

## random-survival-forest-shap.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "First, we need to load the data and transform it into numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "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 = np.column_stack((Xt.values, grade_num))\n",
    "\n",
    "feature_names = X_no_grade.columns.tolist() + [\"tgrade\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "random_state = 20\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    Xt, y, test_size=0.25, random_state=random_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomSurvivalForest(max_features='sqrt', min_samples_leaf=15,\n",
       "                     min_samples_split=10, n_estimators=1000, n_jobs=-1,\n",
       "                     random_state=20)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rsf = RandomSurvivalForest(n_estimators=1000,\n",
    "                           min_samples_split=10,\n",
    "                           min_samples_leaf=15,\n",
    "                           max_features=\"sqrt\",\n",
    "                           n_jobs=-1,\n",
    "                           random_state=random_state)\n",
    "rsf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We can check how well the model performs by evaluating it on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false,
    "scrolled": true
   },
   "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</th>\n",
       "      <th>menostat</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>0</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>1</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>2</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>3</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>4</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>5</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  menostat  pnodes  progrec  tsize  tgrade\n",
       "0  33.0     0.0    0.0       0.0     1.0     26.0   35.0     2.0\n",
       "1  34.0    37.0    0.0       0.0     1.0      0.0   40.0     2.0\n",
       "2  36.0    14.0    0.0       0.0     1.0     76.0   36.0     1.0\n",
       "3  65.0    64.0    0.0       1.0    26.0      2.0   70.0     2.0\n",
       "4  80.0    59.0    0.0       1.0    30.0      0.0   39.0     1.0\n",
       "5  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": [
    "a = np.empty(X_test.shape[0], dtype=[(\"age\", float), (\"pnodes\", float)])\n",
    "a[\"age\"] = X_test[:, 0]\n",
    "a[\"pnodes\"] = X_test[:, 4]\n",
    "\n",
    "sort_idx = np.argsort(a, order=[\"pnodes\", \"age\"])\n",
    "X_test_sel = pd.DataFrame(\n",
    "    X_test[np.concatenate((sort_idx[:3], sort_idx[-3:]))],\n",
    "    columns=feature_names)\n",
    "\n",
    "X_test_sel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7+ElEQVR4nO3de3xU1bnw8d9DAkQg3CFgQC6FCEEoaAAVtaH1QpGWagUU6wVb8T2nvq2X9tSjp1I5Pdb2rba2tla0UmtFkLYcaKRUbI0VpFxElGsCcpFwDYQhCSFAkvX+sWeHnclcdjKzMzOZ5/v55ENmz549K5skT9Z61nqWGGNQSimVutrEuwFKKaXiSwOBUkqlOA0ESimV4jQQKKVUitNAoJRSKS493g1oqp49e5qBAwc2On7q1Ck6duzY8g1KMHofLHofLHofLHof4IMPPjhmjOkV7LmkCwQDBw5kw4YNjY4XFhaSn5/f8g1KMHofLHofLHofLHofQET2hXpOh4aUUirFaSBQSqkUp4FAKaVSXNLlCJRSKl7OnTtHSUkJ1dXV8W5KSBkZGfTr14+2bdu6fo0GAqWUcqmkpITMzEwGDhyIiMS7OY0YYzh+/DglJSUMGjTI9es8GxoSkZdF5KiIbAnxvIjIL0Rkl4h8LCKXetUWpZSKherqanr06JGQQQBAROjRo0eTeyxe5gh+B0wK8/wXgaH+j9nA8x62RSmlYiJRg4CtOe3zbGjIGPNPERkY5pSpwO+NVQf7XyLSVUT6GmMOedGeB166joN1pQ2Olbfpyom0HgD07Nie3p3bN3rd5MGTmZYzzYsmKaVUQohnjiAb2O94XOI/1igQiMhsrF4DWVlZFBYWNrpYZWVl0OO2z7xzEQNrL3QcMZgLzvH+5duoM8AZqDzR8DU7MurYcGQDxUXFTMic4PLLiq9I9yFV6H2w6H2wxOo+dOnShYqKiugbFIWVK1fyve99j9raWu666y4eeuihRudUV1c36etNimSxMWYeMA8gLy/PBFshGGnl4N4/LKOiylf/+Gytj3anM1ncrh9HKqo5Vnmmwfkjzm5mcWZH5vbswUbfOzz2pcdi8aV4TldQWvQ+WPQ+WGJ1H7Zv305mZmb0DWqm2tpavvvd77Jy5Ur69evH2LFjmTZtGrm5uQ3Oy8jIYMyYMa6vG89AcADo73jcz3/ME3e/9EyDx7+ceQ/nak+w6NORAAyfkM+oa8+nNNYufprcnUsYc/oY1W2OedUspZRybd26dQwZMoTBgwcDcOutt7J06dJGgaCp4hkIlgH3i8hCYDxw0qv8QDAXtB1IDekc21/J2arDnPKdbRAIxk97GHiY6nmj2Z9+junzRtM2rQ3tOmcx+bP3aN5AqRT3xF+2su1geUyvmXthZ+Z8aUTI5w8cOED//uf/fu7Xrx9r166N+n09CwQi8jqQD/QUkRJgDtAWwBjzG2A5MBnYBVQBs7xqSzA5bdqwX8bRRjpwSP7JiUP7+NU37m9wzmcum0Be13FQsY66OkO7ugqKKs7ARy9rIFBKtRpezhq6LcLzBvimV+8fychJOQwoKIBq+FeFcLRTlwbPV1fsZWvhXibdez//cds8ZrywhoGH/gx9/pez5Ufi1GqlVKII95e7V7Kzs9m///wcm5KSErKzs6O+blIki73QbcZ0us2Ybj24406ghgEvPVf//Cvfm8exvctY+aJ1bOroXJZyM7XmLQ60OcOs3+VZJ3bspUNFSqkWMXbsWHbu3MmePXvIzs5m4cKFLFiwIOrrpmwgCFS9Ywf77riz/vGg9Bx8Ha6lpuptVr74HCPyb2fRv93Gt349ijbnNkM7oPokG+QMG9bMBdBgoJTyVHp6Os899xw33HADtbW13HPPPYwYEX3PRKuPAp2nTCFj2LAGx7K2LGMMlfQc+GUAPvlgNQBH0v6d7Qd/TNWZX3LF6Vt5/NhxAJZ/9HLLNloplZImT55McXExn3zyCY89Fptp7dojIGCYyG/fHXfSa81rjHziCRYcH1h/fOro8+NxTx65nNvSTpHX8S9gGq5aVkqpZKE9ghA6T5kCwOE5czDnznG26jCLnniESyq2sei+K1h03xU8edNIXq/9AuV04GxtXZxbrJRSzaM9ghDsHsLhOXNoX92dmrZwsGhXg/UGM8dfBMCfPoRdbc4xa8X5GbBao0gplSw0EIRhB4NBv3mL0iGf50j1e5yuONvgnJnjL6J0TXs2cQIqDkNmHzYc2cCGIxsATSArpRKfBoIIus2YztCCAvpv+jlvDb2E6kpfg4Vnn7lsAheYG5h/+Bdw+CgMuIrFdGdumzLm6mwipVQS0ByBC/asot6VtbRr07X+eHXlIbatepetVdczr8u3YMBVAEyjE4/7qgBYvnt5PJqslFKuaY/AhfpZRQELz15+8CFOHCym7/6N/LLPVfy9w+T61/xC/ou8059aw0VKKRUj99xzDwUFBfTu3ZstW4JuANlkGgiikHfj9ax8sZj2ldu5tfKzHD1Ww+Ge6azdU8bP0z4L/T9l+6kDDZLINk0mK6Wa4+677+b+++/nzjvvjHyySxoImsi5ArkL0CWzLycrSuh0fBNdK0eRn9GRz2d0Y1vbm8g99U9OtalpdI2isiJAcwdKqaa75ppr2Lt3b0yvqYGgCey1Bbaq9evp3z2Tk/17U1P1Nulpuzi4Hc6dqeXSC8eQJxn8e+VuOrY9ar1g5C2QNytoD0EplWT++ggc3hzba/YZCV98KrbXdEEDQRMErkA+segNmDOHdgMHcvTiwfXHDxbtovL4x6zuM5Gqs7V0OHSSged2U1lRTVaeBgGlVGLRQBCFbjOmU15QQJ/16xk99Zb6IGFPL+105b38dJOVQP7OoYfIPbEd5t8Ineogs0/c2q2UioE4/OXuFQ0EUeo8ZQpV69dzeM4cygsKAKirTeds7VG6b1zNon+ztmWY97Nr6XD6HUbsWwV9elNUXVo/RKSJY6VUPOk6gih1mzGdPk88QYexY+uP9a6sBc5XLAX4e4fJ3Hr2+8zr8i0mnzrFoLNWbaINRzYwd81cZq2YxeLixS3beKVU0rntttu44oorKCoqol+/fvz2t7+N+praI4iBRtVL77iTZbW9G5xjVy39O5P5zqG3efRMGiO+Pp/FxYtZvnu5ziRSSrny+uuvx/ya2iNoITPHX1RftbRDu7T649NypjF/0nwu7n4xRWVF2jNQSrU47RF4xNTWYM6dc33+5MFWUtkuWBeqNIXmE5RSsaY9Ag/Y6w3OVB/g47dXBD2n/9lP2PrkVaxd/DRwvmfw+BWPk5eVF/Q1dj5BewxKqVjSHoEHus2YToflOzhZvY4Nb75Vv3+BrXLoTezfuYT+Zz5h2+bFLLhoWv3eBtNypoX8i39x8WLmrpnL8t3LtVeglIoZ7RF4JKdNGyS9HycO7ePlBx9q0DMYP+1hRjy6itM9crm8zXYq33/R1TWn5UwjLytPcwlKqZjSQOCRATXFXNSmLZLWixMHi1n54nOseL5htj/ryq8BMObk28x4YQ0L1n4a8bqTB0/m4u4X67RTpVTMaCDw0FhOcsP/+T49B34ZgK2Fr/Grb9zPr75xvxUU8mZxpHseubKP7xx6yFXPIDCXUFRWpHseKJUi9u/fz8SJE8nNzWXEiBE8++yzMbmu5gg8VLV+PRdOWc1dP57Niucz6xeYVVfsZWvhXk4e/Qgq+jO8cwa5fEiH0++4vradS9ACdkqljvT0dJ5++mkuvfRSKioquOyyy7juuuvIzc2N7roxap8K4Cw9ATDp324DrHITr3xvHsdLNnLKd5Yq32nIHEla25PNfi87ZzB58GR60SsWzVdKJaC+ffvSt29fADIzMxk+fDgHDhzQQJCo7JXGh+fMqQ8G9rG8KTdS+NoQTp+GDl0Lonofe/2BvTL5roy7orqeUsqdH6/7MTvKdsT0msO6D+N7477n6ty9e/fy4YcfMn78+KjfV3MEHrLrEIEVEE4segOAEVdnk3/7xQD4jlZxyne22e/hXJmslEoNlZWVfPWrX+XnP/85nTt3jvp62iPwWGDPwK5QeuGUKeTfPoG//QZOV5yFXucXmfXs1J6szIzzF/FvaKOUShxu/3KPtXPnzvHVr36V22+/nZtvvjkm19RA0ALsYGAHgeodVndyxKvTKXzF+i+wF5lVVNdQUV3DscozAE3a0KaorIhn2zzLKyte0VIUSrVCxhi+/vWvM3z4cB566KGYXVcDQQtxVijdd8ed9XsfmzYXIG3bMn7aw8DDLFj7KUs3Hah/3YMHHuTysg2sXfy0/5zg7FyBz+drUK9IA4JSrcfq1at59dVXGTlyJKNHjwbgySefZPLkyVFdVwNBHNi1iKrWr6fus2MbJGpmjr+ovtwEwNrF02DrXMzmxcwouxKwSlo7z4Hz00kLCwspvbCU5buX1wcE+3mlVHK76qqrMMbE/LqaLI6DbjOmM+DV39cnks9WH2XRE48ELVA3ftrDHOmex8i0T3n8+HddLTxzLjoDdMGZUiosDQRx1G3GdDpkDKFtWjdK9+5h++rCoOdlXfk1Ol40hhF9u5Ar+5jgcuGZXZtIKaXC0UAQZx3bDaFnx2vpNXAQpXv3BO8Z5M2CWW/CrDfZ23ZwfBqqlGq1PA0EIjJJRIpEZJeIPBLk+YtE5B0R+VBEPhaR6DIeSWz4hHx6DRxEybYtrHzxuZD7GABUna11XaQOrH0MtDCdUioUzwKBiKQBvwK+COQCt4lI4Dro/wLeMMaMAW4Ffu1VexLdqGsnMWPOU1x37/0AIYeJenZqT4d2aazdU8ajSzZHDAj2bCLd0EYpFYqXPYJxwC5jzG5jzFlgITA14BwD2MviugAHPWxPUhh17ST65V4S8vmszAxG9O3CkzeNZPyg7mw7VN5gummgaTnTNGmslApLvJiKBCAitwCTjDHf8D++AxhvjLnfcU5f4C2gG9ARuNYY80GQa80GZgNkZWVdtnDhwkbvV1lZSadOnbz4Ujy1/5VPOZXWg461x+mWVU2nSTkULbW+voun3tro/NEfPkbXk1soyvl3Dl14Az9aexqA/xx/ARD6Pjx7+FkOnD1Adrts8jrmMSFzgodfVfwl6/dDrOl9sMTqPnTp0oUhQ4bEoEXNU11dzaRJkzh79iw1NTVMnTqVxx57rNF5u3bt4uTJhoUsJ06c+IExJujskXivI7gN+J0x5mkRuQJ4VUQuMcbUOU8yxswD5gHk5eWZ/Pz8RhcqLCwk2PFEt/7DpezadpLytB6kHz/JlPx8jry7gpJtW+heU91om0s6fQMKHuDi4l9zcU4Oz3cdBkB+/hVA6PtQWnx+bcGuM7vIuTinVa8tSNbvh1jT+2CJ1X3Yvn07mZmZ0TeomTp16sS7775Lp06dOHfuHFdddRVf+cpXuPzyyxucl5GRwZgxY1xf18uhoQNAf8fjfv5jTl8H3gAwxqwBMoCeHrYp4Yx9cCq3vTiTznI+eg+fkA+EyBPkzYIpP7c+3/xH1+8TuLZAdzdTKvmISH3P5ty5c5w7dw4Rifq6XvYI1gNDRWQQVgC4FZgZcM6nwBeA34nIcKxAUOphmxJauenC6/cuAKBLZl9Ktm3h47dXNO4V5M2ygsDhzTxuvssfz17BjBesFccXRngPuxfgXHls5w60HIVS7h1+8knObI9tGer2w4fR59FHw55TW1vLZZddxq5du/jmN7+Z2GWojTE1wP3A34DtWLODtorIXBH5sv+0h4F7ReQj4HXgbuNV0iLBDcntWN8rKDdd4JzVmQo1e4iRt0CfkQyt28st7dZETBo7BW53CeiWl0olibS0NDZt2kRJSQnr1q1jy5YtUV/T0xyBMWY5sDzg2OOOz7cBrTtr6dLYB6cy1v/56/cuoOyCq+jW3Rf6BXmzIG8W7ebfyIjDm1nY7r9ZXTURq4Pljl2fCNAtL5Vqokh/uXuta9euTJw4kRUrVnDJJaFnGrqhK4sTUHbNbgCqyk5FPtnfMxh4bjejTqykcP85j1unlIqX0tJSfD4fAKdPn2blypUMGzYs6utG7BGIyEhjzOao30m5NnJSDvuXFHMkTULnCWz+nkHlL74Ax07xu61n2f7Cmvqng1UqVUolp0OHDnHXXXdRW1tLXV0d06dPZ4q/mnE03AwN/VpE2gO/A14zxjR/l3XlSrcZ0+m/4occYxg1NSWsfPE5gNDBAGuhWVbZBh7I/AdruAWAbYfKAVwHgqKyokZDRJpAVipxjBo1ig8//DDm1404NGSMuRq4HWsq6AciskBErot5S1QDIyfl0PN0Bl0yxgFhksb1L7B++d/d4V8suu8KFt13Bbl9O7N2T5mrmkSTB09utO/xhiMbdJqpUinAVbLYGLNTRP4L2AD8Ahgj1uTVR40xf/aygamq24zptHl7AWfScujSdn/kF9hTSv3jh2ANC63dU8bSTQci9gqciWPb4uLFLN+9nKKyovpzlFKtT8QegYiMEpGfYU0B/TzwJWPMcP/nP/O4fSltSG5HAE5X1dbnCppi5viLGD+oe7Pf355mGthTUEq1Lm56BL8EXsL66/+0fdAYc9DfS1AeGfvgVHbdu4Cj7UdB1dusfPE5tq8uZPiE/LD5Ai8Eyx84aS5BqeTlJhAsMca86jwgIt82xjwbeFzF3pDcjpQVjaJjXTmZg6so3bsHCJ047lS5B+bfWP/4C1WX8uSRy1mw9tNmzx6yS1mHoiuUlUpubgLBncDPA47dDTwb68aoxuxeAe2GMGPOTBY98UjoKaUjb6HS56Or/fjwZqZ2rOZJLneVJwglWP7Ayc4lwPmgYL9OKZX4QuYIROQ2EfkLMEhEljk+3gHKWq6JyilSQbpNY/6nfltL+owkKzOD8YO6u5491Bx2LsFZ1E7LVSjlndraWsaMGROTNQQQPln8PvA0sMP/r/3xMHBDTN5duVZXUc6JRW9E3LgmmKmjswFc1yKKxrScaeRl5dXnFHTaqVKx9+yzzzJ8+PCYXS/k0JAxZh+wD7giZu+mmiW9Rw+OpvXhg98sYGhBATWd23L42JHwK45t+1Yxc+TfWTpoWH2vwOuVxnZOITB3YD+nQ0ZKNV9JSQlvvvkmjz32GM8880xMrhkyEIjIKmPMVSJSgbWlZP1TgDHGdA7xUhVjuTdewtHXiii6eCZHKnbTvvwjwBoeChsIRt4C+1bB5j8ydfTz9fscg/vVxs1h5xScuQNoGBg0IKhk994bxRzbXxnTa/bs34mrp+eEPeeBBx7gJz/5CRUVFTF733A9gqv8/8ZvOx4FwIirraGd4nVHOFrcjzZ1Heg3pCOle/ew6IlH6s8bPiEf0jPOv9BeZMb5X/yPLtncIsEAGieZdYGaUtEpKCigd+/eXHbZZRQWFsbsuuF6BGFXIhljNGHcgkZcnc2Iq7PrN66xk8a2km1bKNm2hYs+dx2E2JIvMBgs3XSgRYvS2YFBS16r1iDSX+5eWL16NcuWLWP58uVUV1dTXl7O1772Nf7whz9Edd1w00c/wBoSCrYPmgEGR/XOqtnK0vpwZnMFM+Y8VX/s47dXsPLF5yjbuT3sa+1f+ks3HWDtnrL6EhTQspVKnQvUdJhIKXd+9KMf8aMf/Qiw9mH+6U9/GnUQgPBDQ4OivrqKuSG5HVlXBLu2narfyAasBWbbVxfW1ypvYN8q2DDfGirCCgYzx1/EgrWf1geBplYqjYZzgZrmDZSKv3DrCIb5/7002EfLNVE5jX1wKt1rD1OW1oet77mYDuqvShpso/uZ4y9qVqXSaAWuO7Cnm+raA6Xcy8/Pp6CgICbXCjc09BAwG2vtQCCDVXROxUF2zW7K0vqw7c0t9YnkkOyEcUCvIJBdqdTOHQR73oveguYNlIq/kD0CY8xs/78Tg3xoEIijkZNy6Oorpub48UbPnT5WyqInHmlYqTRMr8A2c/xFPHnTyKDVSrcdKm+RxWgbjmzQBWhKxYGbrSozgH8HrsLqCbwH/MYYU+1x21QI9l4FgYZPyMfn89XPIDpfqdTfKzi82SpIN/KWoD0DO3cQaMYLa9h2qJwZHm6BOXnw5PpcgeYJlGpZborO/R6owCpHDTATeBXQn9Y4KzddeP3eBaT36EF6r17kjBvJxVMz6F5TzfbVhfUBAWCU3Ss47N9+OsQQUTB2iQpb4Gwj+5xoAsO0nGmaI1AqTtwEgkuMMbmOx++IyDavGqTcGZLbkV3bTlJXUc7ZinLKThhqSkvpNak9o66dxKhrJ9VPKd2+upBRc56yfvnPvzFiviBQYE/BOdsIGgeGlpyGqpSKnptAsFFELjfG/AtARMZjbVmp4mjsg1MZC5xY9AblBQWsqkznqAwmbZeBfOsce0ppg7LVdtmJgges4aIQw0ThhAsMGhSUSj7hVhZvxsoJtAXeF5FP/Y8HYFUkVQmg24zpdJsxnX2zfoiPwVR87Gvw/PAJ+fX5glHXTjr/S9+eSeSvRQQ0KyhAw8AQr7UJSqWKgQMHkpmZSVpaGunp6WzYEP3f5eF6BLEpdK1axMhJOexfUkzdBR0aHLd7BQ3kzbI+Nsw/HwSakTsIxhkUnMllt3TFsVKRvfPOO/Ts2TNm14tUhrqeiPQGMkKcruLMnklUV1MT9Pmgu5rZAQGs3IE9q8ipmb0EW+BsI1uwIaNQK47D0WChVPTcTB/9MtaisguBo1hDQ9uBEd42TcVKo+GhYOxZRU7OoaNmBITA2Ua2UENGzmqlgSWsg9Eqpiqe3vndPI7u2x3Ta/YeMJiJd88Oe46IcP311yMi3HfffcyeHf58N9wki/8buBx42xgzRkQmAl+L+p2VJ0617cWSpzeSMy6rftWxPTxkl6221hYEBARn78BmDx0F5hLAVWBoyroEaNhLiLRPMsCsFbPqF6FpMFCpYtWqVWRnZ3P06FGuu+46hg0bxjXXXBPVNd0EgnPGmOMi0kZE2hhj3hGRn0f1rsoT2TW7qeEiDu5sz8GdPorXHakPCHbZ6gZrCyLtbhYslwDNDgy2YD2F5iSW7UVoc9fMre896FCRaimR/nL3Sna29fPTu3dvbrrpJtatW9cigcAnIp2wVhS/JiJHgVNRvavyxICaYnps+TPHR9/MgfTBlO6pA6y9DIKuLYgUCGyBvYXAwNDERHOwnsKMF9Y0eStN+xe+HQSceQWfz8crK15xdR3QAKKSw6lTp6irqyMzM5NTp07x1ltv8fjjj0d9XTeBYCpQDTwA3A50AeZG/c4q5jpPmYLP52NATTG91rzGxtHf5mBNDlvfO9BomCgqgYHBmWhuZnLZLnq3dNOBJvUKmppXCEZLYatkceTIEW666SYAampqmDlzJpMmufyDLoyIgcAYc0pE+gDjgDLgb8aYxtXOVNx1mzGdj7J689n8fE4seoOs37yFr2tO0CqlQWcRNZedaI4iuTxz/EVRF7ZzBoXCwkLyQ+zUFki30FTJYvDgwXz00Ucxv66bWUPfAB4H/oG1W9kvRWSuMeblmLdGxUy3GdO5DDiypJia2s4NnrNnEa188TnARa4gksBcQrAcArgKDl4XtwtGS2GrVOdmaOi7wBi7FyAiPYD3AQ0ECc5eW2BvYuMcHgJY+eJzsQsGEDq5DK5WMUcqbtcSQcG5oM2mw0WqtXMTCI5jVR+1VfiPRSQik4BngTTgJWPMU0HOmQ78AKt8xUfGmJlurq3cCbWJjTMYNClx7Ea4qagQMrnclBpGtlgGB+eCtvpmByxs06CgjDGIBNvKPTEYY5r8mnC1hh7yf7oLWCsiS7F+WU8FPo50YRFJA34FXAeUAOtFZJkxZpvjnKHAfwITjDEn/KuXVQzZpScCh4eg8foCIPgag1hoxirmUDWMbLGuZRRs7YIzAa05BJWRkcHx48fp0aNHQgYDYwzHjx8nI6NpRSDC9Qgy/f9+4v+wLXV57XHALmPMbgARWYgVRJwlrO8FfmWMOQFgjDnq8trKpVCb2Njs9QUApXv3ADEaJgon2CrmCFNQQ0059ZozOGgOQfXr14+SkhJKS0vj3ZSQMjIy6NevX5NeE67W0BPOx/61BBhjKl1eOxvY73hcAowPOCfHf+3VWMNHPzDGrEC1GHt9AVDfK/BcsKGjZuyTEA+6kjm1tW3blkGDBsW7GTHnZtbQJVg7knX3Pz4G3GmM2Rqj9x+KVUG/H/BPERlpjPEFtGE2MBsgKyuLwsLCRheqrKwMejzVBLsPNTU1lLfvR8Ej8+g0KSfka30+H5UHS1j862fplftZj1vaUN/2o7iYVfjee4lNle5+0Hy+0xSdqOMHr64kv3/bBs958f0w5NwQNmCtZF6w0epl5XXMY0LmhJi+Tyzpz4VF70N4bpLF84CHjDHvAIhIPvAicGWE1x0A+jse9/MfcyoB1hpjzgF7RKQYKzCsd55kjJnnbwd5eXkm2Pzwpswbb82C3Yf1Hy5lXRHs8w1h5CurGDkph24zpjd6bfeaala++Bx1pYfIz/92C7XYlg/zP6YruP5/PHjBpzy6ZDPbqzrxg/wrGjznxfdDPvnkFOc0WMm868wuci7OSdgegv5cWPQ+hOcmEHS0gwCAMaZQRDq6eN16YKiIDMIKALdi7Xfs9L/AbcB8EemJNVQU23J+irEPToWfLWVdUSab218JK97nmhmNz4vJquNoBSaRw6w9sBehBStg5/Od5vki9zkEt7OPAlcyz10zt0GtIyedYaSShZtAsFtEvo81PARW5dGIv6yNMTUicj/wN6zx/5eNMVtFZC6wwRizzP/c9f49kGuB7+qqZW+MfXAqHd47QOFrRWxufyU9HOsKEkZgEjnUwjSHh7p/nmcidk7Da+7so8BaR046w0glEzeB4B7gCeDPWNNH3/Mfi8gYsxxYHnDsccfnBnjI/6E8NuLqbI6/PJ/N7a8MWnbC1iLTSYOJVNwu0OHNjO8Di+57uNFT1lDAFUFe1JizLHZT1yWEKpetM4xUMgkbCPxrAf5sjJnYQu1RHgu3rgDiNJ00lGCzi5wC1yE0k72iOdaL1gJXKetQkUpUYQOBMaZWROpEpIsx5mRLNUp5J9K6gsDppM7egVOL9hQ8Zq9RCLZorbllLgJXKWuFU5XI3AwNVQKbRWQljn0IjDHf8qxVynPlpguv37uAIbkdrWRyEM7egVPcewoeCbZozRkcmpJLCBwy0gqnKpG5CQR/9n+oVmJIbkd2bTtJWVof1hUBP1saNBg4ewdOzp5C3HsGHi9CcwaH5myeY9MKpyqRudmP4BURaQcMw0oWFxljznreMuWZsQ9OZSyw3j+ldF1RJiVPbwRosNdxKM3a9tILI2+xAkHBA2GrmsZKczfPCRSswiloDkHFj5uVxZOBF7DqDQkwSETuM8b81evGKW+NfXAqp2f9kAPpg4GuHCuxqodECgRRbXsZS/YvfDsIOKabjvb5YE9X99dyEUCCrVtoaiI5WIVT0OmmKr7cDA09A0w0xuwCEJHPAG8CGghagQE1xWRtWUZG9TDez5hETWkP16+1F6DFdLezpnLOLIo03TSUJuy57NwzoTnrD8JNN9VZRipe3ASCCjsI+O2m4f4EKol1njKl/vO6qiqOBmxiE4m921ncegVOjqCwqSklBZpQ8C4wZxAr4WYZhTpfg4SKFTeBYIOILAfewMoRTMPaW+BmAGOMJpKTWLcZ0+vrDu2b9UN8DKZ43RHXgSAhylJEy5lrgLhUPw01yygYHUZSseYmEGQAR4DP+R+XAhcAX8IKDBoIWokBNcX1+YKUYv/iL3igYeLZFiZ/4NUey6GGkEBXLavYczNrSL/rUkhdRTmle8pY0oRZREDiTCdtrsDEs81OQDvP8YvnHsu6t7KKJTc9ApUiOk+ZQtZv3qI0szPQ3fUsooSZThqtUHst2z0F+xy/5uyxDNEHiGAzj3S4SEVDA4Gq123GdIYWFDC0egUDHp7Jkqc3cqyksr53AMF7CAkzndQLgcNGzmMBIu2xDLHZZznYsJEOF6loaCBQjVTv2MG+O+6kR3oONVljgU4AEXsIrSJxHEwTgoEtWLkKaJl9lpVqqpCBQETCloY2xjwT++aoeHNOJ83asowBNcUMeOr3AA16COFyB3FdV+CVUAnlZqxkbm7J60iC5Q18Ph+vrHjF9TU0z5CawvUIMlusFSphNJhOesedDZ7LGZcFwMGdPg7u9FG87kijgJBQ6wpiLdRKZudzEYQqee1F3qCpNM+QukIGAmPMEy3ZEJWY7GEisAaILgV6pOdwJGtsfUCA88NF9vBQ0s8iCiVwJbPdO3AZCIKVvLaDgv18c4SabtqUvXqDrW6ORHsQrYObWkMZwNeBEVhrCgAwxrjapUwlL+cwkVOvNa/Ri9coveJ2Nre/ksLXrL8k7WAQOIvIzhu0yqCw+Y+N91m2Rdhv2ZlYfnTJZh5dsrlBgtnL6afBNLVXYa9+Bu1FJDs3yeJXgR3ADcBc4HZgu5eNUonBOUzkdGLRG5QXFDCgppize/dSdPHMBltfOmcR2UGgte5h0GifZVsT6hfZv+ydQcA5dNRSASHcIrZgFhcvZu6auSzfvVwDQZJzEwiGGGOmichUf0nqBVj7FqsU5QwQnRe9wZEQW19G2u2sVfQQQm2nOf/Gxj0Flz0EOD/9NNbbZ8bStJxp9ZvtuBlO0mGkxOUmEJzz/+sTkUuAw0Bv75qkkkmkrS9tgbudBQ4b2eckfWCwBfYUHCWyG5wTITC42T4T4hcY3A4nRSqiF837a3CJnptAME9EugHfB5Zh5Qy/72mrVNIpb9M97LTSwN3OnMNG0AoDQ2BPIbBEtsvAEGn7TIhNsrm53A4nhSui11w6yyl23ASC+caYWuBdYLDH7VFJKLtmN3VnqzhaDDWlpa5qE7kNDPa5SS9SYGhiTiFwGMlONMd7uCiUpuYf3NDV1LHjJhDsEZEVwCLgH8YY43GbVJIZOSmHAQUFrKq8hpq6Ds26RrDA0OrKVTgFBoZgOQVbhEVr9s5pSjWXm0AwDJgCfBN4WUT+Aiw0xqzytGUqadjJ4/fvXUBZWh/W/2wpYx+cGtU1A9cjBJPUQ0eBQs0+CjaEFMQXqi7llyevalTCwuc7zfNF7staJEISuincJKonD55ML3q1UIuSk5sy1FVYm9K84c8VPIs1TJTmcdtUkhmS25F1RbBr2ynGxuB6gQlmp1Y3HTXU7CM322/uW8VsVkEX+DvNX2Eci4J4LclNotpOUg9pP6RBqQ1NMjfkquiciHwOmAFMAjYAjSeXq5Q39sGp7PL3Cpqy3WUogcNFTqF6Ca1OqADh5F/hPLvrRmbP+u8GT1kri69w9VbJVhDPTd7BTlL7fL76Y7oQrjE3K4v3Ah9i9Qq+a4w55XWjVPLKrtlNWVqfBgvMlMfsFc4xELjjmhuJPJxkBwtnqQ1dCNeYmx7BKGNMuectUa3CyEk57A+xwEwltsAd19wItwFPsOsnQsCwF8Kp88KVof4PY8xPgP8RkUYzhYwx3/K0ZSop2QvMyk0XXr93AUNyO0adOFYuBZl1NNrngz1drQcuZh819Rd1qA14AiVb/iHVhOsR2PWENrREQ1TrMSS3I7u2naTcdGHXtpMxSRwH0yr3PWiuULOObC5nHzV1fwW3wWPGC2uaNewUSaL0MpJduDLUf/F/utkYszHUeUoFGvvgVMYCr997vmcAxLR3YO97sPLF54BWNHuouUIklTfZY+NuZh81YVFbUzVn2CmSSMNSGiTcc5MjeFpE+gB/BBYZY7Z43CbVStg9AyDmvQP7F3+rXnQWS25mHwVbzBYjzRl2iiTcsFRgkAhcT7G3XTm16QfCrkFIpSmmbtYRTPQHgunACyLSGSsg/NDz1qmkZvcMgPpeQSy12j2SlSvhgkuk3EVV2Sg6dA997VSrY+RqHYEx5jDwCxF5B/gP4HFAA4FqkrqqqkbbXwbqPGVK0D0Qwgm2+tjn83Hk3RUNjrWqlcheCVXmApq1P3O8BAaJwPUUM14Azl7P/EnB11ikWh0jN+sIhmMtJvsqcByr5tDDHrdLtTLpPXpQA1Ad+pzqHTsAmhQIwq0+dmp1Rey8EC7h7DbZHOn6SRJIUo2bHsHLwELgBmPMwaZcXEQmYZWkSANeMsY8FeK8r2LlIMYaY3SWUiuU3qsXvjMXsHHYA/XHAktW77vjzgZ7JLvpHYRafRy4V2+rL2IXC+HyCG6SzeF4mIhW0QsbCEQkDdhjjHm2qRf2v/ZXwHVACbBeRJYZY7YFnJcJfBtY29T3UMkjZ1xWg8fBNr537pFctX49VevXA03rIYSi+YQouUk2h+NhIlpFL2wgMMbUikh/EWlnjDnbxGuPA3YZY3YDiMhCYCqwLeC8/wZ+DHy3iddXSWTE1dkN/vrf+t4BCl8rovC1IorXHfEfHQKjHwCgJruUbmsXw5w5lBcUAM3LHwSy8wmaL4iDZpbZVt5ztR8BsFpElgH1dYaMMc9EeF02sN/xuAQY7zxBRC4F+htj3hSRkIFARGYDswGysrIoLCxsdE5lZWXQ46kmme5D3zzh5D7ToCCYrfpUW6pG3kivTkc56/PRbudOqtav59M//OH8OePGcvrqq4NeO9h9aNOrL219vvp8wZo3/7fB892HDqdX7mej/bISSqJ8P/RtP4qsDB8E/F93PbkF9q3C995LzbrukaxrOHThDRHPC7wPPt9pgJD3xv6eTIR71xLcBIJP/B9tgMxYvbGItAGeAe6OdK4xZh4wDyAvL884x35tgWPCqSqp7kN+6KeWPL2Rgzsh/T+eZ8TV2ZxY9AblBQXY295U79hB1+KdDPh+8F1Tg94H/+PA3dDASiZXHiyhrvRQo2slc+8hcb4f8oMf9uceujbnkoc30/XMx1yc/6OIpwbeB3tNQajKrHbJ6sS4d95zs47giWZe+wDQ3/G4n/+YLRO4BCgUEYA+wDIR+bImjFXOuCwO7vQ5ho7ODxsBVGfsILtmNwOace1gCeZgwQF0tpHnosk9hNvVDSIOOYUrebG3XTk9O7ZvXruSkJvpo+8AwYrOfT7CS9cDQ0VkEFYAuBWY6Xj9SaCn430Kge9oEFBwPoF8Pn/QUHmb7i5XwbgTavaRPdvInnEEyd1DaFXCTXeNMEspUsmLqjM1HGtuu5KQmx+l7zg+z8BaT1AT6UXGmBoRuR/4G9b00ZeNMVtFZC6wwRizrDkNVqkjMMHs5MVK5WDsX/h2EGh1O6Mls3C9iYDeQoMqrMDMkbcw877QvYXx89Oplv0Jt7BsWPdhfG/c92J+XTdDQx8EHFotIuvcXNwYsxxYHnDs8RDn5ru5plK2uopyTix6IybTS8Nx9hYWPfFIg5XM2jtIUFH0FgC61I5Lqc143QwNOStytAEuA7p41iKlXEjv0YOjaX3YvOJ9rpnRcu/rXMms+YMEFtBb2ORMFrtY09Ct9hq61V4TsgRFa+NmaOgDrByBYA0J7QG+7mWjlIok98ZLOPpaEftP92yRXoHN2TvQ1cqtmxf7J0Qr98LOzPnSiJhf183Q0KCYv6tSURpxdTbb3tzCUXJavFdg09XKrZcX+ycksnBbVY4F9vsrjyIid2IlivcBPzDGlLVME5UKzu4VHEgfHNd2BKt+GozmE5KHF/snJLJwPYIXgGsBROQa4Cng/wKjsRZ3RdgbTylvjbg6m49//26D8taxKEPRFE2tfhquB6GBQsVLuECQ5virfwYwzxjzJ+BPIrLJ85Yp5YKzvLVdqK68oIBuPh/7fvtyo/NjHShCrT8IFGrBmi0wUGhQUC0pbCAQkXRjTA3wBfy1fly8TqkW4yxvXZNdSs3x4wDU9KwhPb3ht2nvkvcZUlDQoj0GW6SA4QwUulahBexbZZW30GJ3QPhf6K8D74rIMeA08B6AiAwBTrZA25SKyFneOr1XL9J79QKsomGdunatf+7gTh9lA26GfX+GELuktfSwklPgWgXloZG3nN9kRwMBECYQGGP+R0T+DvQF3jLG2GUm2mDlCpSKu1Crj60iY5fWP7bLXu8YcDNHaw9bdYpqiuufdw4rBYpngFAeyJsV3SY7rVCk/Qj+FeRYcbBzlUpkztpFB3dCWVofjl8yuf55e1gpVICA2GyQo1Qi0rF+lTLs3sPW9w40KmaX3qsXR31tgwaIs3v3kPWbtxgaww1ylEokGghUygk1nBQuQPi65tDuzPtkbbFqJXodCLSekWpJGgiU8gsXIApfK2Jz+ys5MHowvUveb5BwjnUPIVg9o0grmDVYqGhoIFAqAmd+oXRPHfS7kiHVKwBrpzSIbQ8hsJ5RpCCg001VtDQQKOWC3VtY8vRGjpW0Y+OwB4Dodkpzw82CtUVPPELJti18/PYKDQaqWTQQKNUEznULEPud0ppj+IR8SrZtabCLmn2c9Iy4tUslDw0ESjVBYB6hpXZKCydwFzU4P1yU9TntIajINBAo1QoEDiHZO6n5fAvpXlOtQ0YqLA0ESkUpntVPQ7FnHjmHjHRmUQDHnsZBjbwlZUpQaCBQKgrpPXpQ1qY779dNoq6ivMHCs0AtGSTsHsLiXz9LXekhnYYaKNyexuBqX+PWRAOBUlHIvfES/yK03hzc6cPXNYfS2sONzosUJALFKmj0yv0s+fnfdjUN1W2wcErawBGwp3EjLvY1bk00ECgVBWfy+PzK5K6NzgsWJALrGtm8XpsQiptg4RQpcCRtkEhBGgiUipFQK5OhcZA4uNPXqK6Rzeu1CaG43WTHFi5wNKV3kbABI1IOIR48yltoIFCqBQQGiWB1jWyJsDbBjXCBw23vojnDUU6eBZFIOYR42LfK+lcDgVKtQ7jeQyKsTYhWrLbwDMfToalIOYR4+Kt3GxZpIFAqASXilFQvNHU4yilcEGmV9Ze++JRnl9ZAoFSCSe/RgxqAam8Sx61FuCCi2302jQYCpRKMvffygIdnsu+OO6nesaO+dxBKa+41KO9pIFAqgXWeMiXiOdprUNHSQKBUAjpWUsmSpzcCQ8iZ/aOQiWUgYm9BqUg0ECiVYJylro+VVAKEDQRKRUsDgVIJxjm11NoIx+4dBFedMSkuC9BU66GBQKkEFrgRTjDlbbpTd7aq0RBR5ylTIKu3V01TrYgGAqUSWLiFZ7bFj/yNmroOUH3+mJ1A5uv3eNg61VpoIFAqyaX36oXvzAX1+yiDVa+od8n79Hz6Gfb99mVAp5iq0DQQKJXkgg0fVbTvDf2upOfhTwGoWr+eqvXrKY9QBluDRWryNBCIyCTgWSANeMkY81TA8w8B3wBqgFLgHmPMPi/bpFRrE2z4yEoyt2Pj6Afo2rUrNdml1Bw/HvY6Td0zIZAGkeTlWSAQkTTgV8B1QAmwXkSWGWO2OU77EMgzxlSJyL8BPwFmeNUmpVKF3Uvw+XzA+dXK4YTbWAdC758A7nocGigSl5c9gnHALmPMbgARWQhMBeoDgTHmHcf5/wK+5mF7lEoZdi+hsLCQ/PxLXb0m0sY6ofZPAKjJLiXryPpmBQoNEPEnxhhvLixyCzDJGPMN/+M7gPHGmPtDnP8ccNgY88Mgz80GZgNkZWVdtnDhwkavr6yspFOnTjH8CpKT3geL3gdLrO5D2S7DyX2hf1dUlVr/dnB0OroMELoPEQAueO89Mtatb/S6djt3AnB26NCo2+i0oX0ddZmZDL3rPkC/HwAmTpz4gTEmL9hzCZEsFpGvAXnA54I9b4yZB8wDyMvLM/n5+Y3Osf7yaXw81eh9sOh9sMTsPkS4ROBGOwd3+qgqNcjJLgCc6HAjOf9xT6NcxolFb1BeUECH6FvYQJuTh0k/fbr+a9fvh/C8DAQHgP6Ox/38xxoQkWuBx4DPGWPOeNgepZRHIu3AdnCnj4M7fUF2ZYtcS6k52tz6FU6a2vpy1D6fjyPvrkjcbTHjzMtAsB4YKiKDsALArcBM5wkiMgZ4AWsI6aiHbVFKtSC3W3OGDhDROdf+UjLObm9wzLmjmQaEhjwLBMaYGhG5H/gb1vTRl40xW0VkLrDBGLMM+H9AJ2CxiAB8aoz5sldtUkrFR6gV0uH2bo5GbcY4Orcbwow51t+ehYWFdK+pZvvqwta5e1mUPM0RGGOWA8sDjj3u+PxaL99fKZXY3JTQaI5g+z7bO5rp7mWNJUSyWCmlWlLp3j1JGRC8GtLSQKCUSinDJ+THuwnN4uWQlgYCpVRKCbfpfSLzsgfTxrMrK6WUSgoaCJRSKsXp0JBSSiWB3gMGe3ZtDQRKqVap3HSpn0ZaU1PD66+dn1I6JLcjYx+cGq+mNcvEu2d7dm0dGlJKtTpDcjvSWU4Gfa7cdGHXtlMt3KLEpj0CpVSrM/bBqYx1PHYWnQu22CzVaSBQSqWcuqoq9t1xZ7yb0WRe7d2gQ0NKqZSS3qMHbTrEuvC196p37Ii453RzaY9AKZVS0nv1wnfmAjYOeyDeTWmS6owd1nahHlxbA4FSKqXY+zknm7K0Pp71ZDQQKKVSilcVT7323hvB94OOBQ0ESimVBK6enuPZtTVZrJRSKU4DgVJKpTgNBEopleI0ECilVIrTQKCUUilOA4FSSqU4DQRKKZXiNBAopVSKE2NMvNvQJCJSCuwL8lRP4FgLNycR6X2w6H2w6H2w6H2AAcaYXsGeSLpAEIqIbDDG5MW7HfGm98Gi98Gi98Gi9yE8HRpSSqkUp4FAKaVSXGsKBPPi3YAEoffBovfBovfBovchjFaTI1BKKdU8ralHoJRSqhk0ECilVIpL+kAgIpNEpEhEdonII/Fuj9dEZK+IbBaRTSKywX+su4isFJGd/n+7+Y+LiPzCf28+FpFL49v66IjIyyJyVES2OI41+WsXkbv85+8Ukbvi8bVEI8R9+IGIHPB/X2wSkcmO5/7Tfx+KROQGx/Gk/tkRkf4i8o6IbBORrSLybf/xlPueiJoxJmk/gDTgE2Aw0A74CMiNd7s8/pr3Aj0Djv0EeMT/+SPAj/2fTwb+CghwObA23u2P8mu/BrgU2NLcrx3oDuz2/9vN/3m3eH9tMbgPPwC+E+TcXP/PRXtgkP/nJa01/OwAfYFL/Z9nAsX+rzflviei/Uj2HsE4YJcxZrcx5iywEJga5zbFw1TgFf/nrwBfcRz/vbH8C+gqIn3j0L6YMMb8EygLONzUr/0GYKUxpswYcwJYCUzyvPExFOI+hDIVWGiMOWOM2QPswvq5SfqfHWPMIWPMRv/nFcB2IJsU/J6IVrIHgmxgv+Nxif9Ya2aAt0TkAxGZ7T+WZYw55P/8MJDl/zwV7k9Tv/bWfE/u9w95vGwPh5Ai90FEBgJjgLXo90STJXsgSEVXGWMuBb4IfFNErnE+aay+bkrOCU7lrx14HvgMMBo4BDwd19a0IBHpBPwJeMAYU+58LsW/J1xL9kBwAOjveNzPf6zVMsYc8P97FFiC1cU/Yg/5+P896j89Fe5PU7/2VnlPjDFHjDG1xpg64EWs7wto5fdBRNpiBYHXjDF/9h/W74kmSvZAsB4YKiKDRKQdcCuwLM5t8oyIdBSRTPtz4HpgC9bXbM90uAtY6v98GXCnf7bE5cBJR5e5tWjq1/434HoR6eYfPrnefyypBeR+bsL6vgDrPtwqIu1FZBAwFFhHK/jZEREBfgtsN8Y843hKvyeaKt7Z6mg/sGYCFGPNgHgs3u3x+GsdjDW74yNgq/31Aj2AvwM7gbeB7v7jAvzKf282A3nx/hqi/Ppfxxr2OIc1jvv15nztwD1YSdNdwKx4f10xug+v+r/Oj7F+4fV1nP+Y/z4UAV90HE/qnx3gKqxhn4+BTf6Pyan4PRHth5aYUEqpFJfsQ0NKKaWipIFAKaVSnAYCpZRKcRoIlFIqxWkgUEqpFKeBQLUqItLDUYHzsKMiZ6WI/NqD9/s/InJnFK//nYjcEss2KdVU6fFugFKxZIw5jlVmARH5AVBpjPmph+/3G6+urVRL0R6BSgkiki8iBf7PfyAir4jIeyKyT0RuFpGfiLXPwwp/2QJE5DIReddf4O9vwSq3+q/1Hf/nhSLyYxFZJyLFInJ1kPNFRJ7z7wPwNtDb8dzjIrJeRLaIyDz/uZ8RkY2Oc4baj0XkKbFq8X8sIp4FO9X6aSBQqeozwOeBLwN/AN4xxowETgM3+oPBL4FbjDGXAS8D/+PiuunGmHHAA8CcIM/fBFyMVTf/TuBKx3PPGWPGGmMuAS4AphhjPgFOisho/zmzgPki0sN/rRHGmFHAD11/5UoF0ECgUtVfjTHnsEoNpAEr/Mc3AwOxfllfAqwUkU3Af2EVI4vELnz2gf86ga4BXjdWgbiDwD8cz00UkbUishkrSI3wH38JmCUiacAMYAFwEqgGfisiNwNVLtqmVFCaI1Cp6gyAMaZORM6Z87VW6rB+LgTYaoy5ojnXBWppws+XiGQAv8aqf7Pfn9/I8D/9J6zexT+AD/x5EERkHPAF4BbgfqzgoVSTaY9AqeCKgF4icgVY5Y5FZESE17jxT2CGiKT5cw4T/cftX/rH/PX162cSGWOqsaphPg/M97enE9DFGLMceBD4bAzaplKU9giUCsIYc9Y/rfMXItIF62fl51hVX6OxBOsv923Ap8Aa//v5RORFrPLRh7HKRDu9hpUTeMv/OBNY6u9JCPBQlO1SKUyrjyqVBPwzk7oYY74f77ao1kd7BEolOBFZwvlZTkrFnPYIlFIqxWmyWCmlUpwGAqWUSnEaCJRSKsVpIFBKqRSngUAppVLc/wdLW4LiS4VSjQAAAABJRU5ErkJggg==\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.event_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": {
    "collapsed": false
   },
   "source": [
    "Alternatively, we can also plot the predicted cumulative hazard function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA700lEQVR4nO3de3xU9Zn48c+TBIgxiSFcwiWUS7kLLWCAqoihVYvU6g+1FbXeaku7q7u1l31tbXe1ul1r77q1W8WK2la52K2VpRTF1Sj1xkVQBAVSEEmAcElCCNdcnt8f50w4mcxMZpI5mczM8369eGXmnDNnvvMlyZPv7fmKqmKMMcYEy0h0AYwxxnRPFiCMMcaEZAHCGGNMSBYgjDHGhGQBwhhjTEhZiS5APPXt21eHDRvW6tjRo0c588wzE1OgbsTq4TSrC4fVgyPd62H9+vUHVbVfqHMpFSCGDRvGunXrWh0rKyujtLQ0MQXqRqweTrO6cFg9ONK9HkRkV7hz1sVkjDEmJAsQxhhjQrIAYYwxJqSUGoMIRUTYuXMnJ06cSHRRwsrOzqa4uJgePXokuijGGNMi5QPEmWeeSV5eHsOGDUNEEl2cNlSVQ4cOUVFRwfDhwxNdHGOMaZHyXUyZmZn06dOnWwYHcFo4ffr06dYtHGNMekr5AAF02+AQ0N3LZ4xJT2kRIIwxJlWtXrqN1Uu3+XJvCxBdYOXKlYwZM4aRI0dy//33J7o4xpgUcnB3PQd31/tyb98GqUVkCPA7oAhQYIGqPhh0jQAPAnOAY8DNqvq2e+4m4N/cS3+oqk/6VVY/NTU1cdttt7Fq1SqKi4uZOnUql19+OePHj0900YwxQTavrmTbmqpEFyMmByvq6Vuc68u9/WxBNALfVtXxwKeA20Qk+LfipcAo99984DcAIlII3A1MB6YBd4tIbx/L6ps1a9YwcuRIRowYQc+ePZk3bx7PPfdcootljAlh25oqDlb489e4X/oW5zJ6WpEv9/atBaGqe4G97uMjIvI+MBjY4rnsCuB36ux7+qaIFIjIQKAUWKWq1QAisgqYDSzqTJnu+d/NbNlT15lbtDF+UD53f/7ssOcrKysZMmRIy/Pi4mLeeuutuJbBGBM/fYtzmfvtKYkuRtT23XcfrAYu+F7c790l6yBEZBgwGQj+zTgY2O15XuEeC3c81L3n47Q+KCoqoqysrNX5/Px8jhw5AkDDqQaampo6+ClCazjV0HL/UI4fP05Dw+lrjh8/zqlTp9q85sSJE23KHk/19fW+3j+ZWF04rB4c3nqo21VDxpEjvPP5HyS0TLHouX07p0aN4gMf/i99DxAikgv8D3CHqsb3z3dAVRcACwBKSko0OCvjhg0byMvLA+CHV02K99u3a+TIkfzhD39oKcOhQ4cYPnx4y/OA7OxsJk+e7Fs50j1jpZfVhcPqweGth0VPPU3zqVMUFBQktEwxmTqVAZddxid9+L/0NUCISA+c4PCUqv4pxCWVwBDP82L3WCVON5P3eJk/pfTX1KlT2b59Ozt37mTw4MEsXryYp59+OtHFMsaEkZGTw9BHf5foYnQLvg1SuzOUHgPeV9VfhLlsGXCjOD4FHHbHLp4HLhGR3u7g9CXusaSTlZXFQw89xGc/+1nGjRvHF7/4Rc4+O/yYhTHGdBd+tiDOB24ANonIRvfY94CPAajqw8AKnCmu5TjTXG9xz1WLyH8Aa93X3RsYsE5Gc+bMYc6cOYkuhjHGxMTPWUx/AyLmkHBnL90W5txCYKEPRTPGpLGaJUupW76cXVmjqcwaQWNjI4uecrp96/Qs8uVwgkvYfaR8NldjjPHatHIbu5tnUttrBAD5jRUt5/LlMCPHp+/+1MEsQBhj0kpl1gjqexYyaFQBo6cVcaApw2ZzhWG5mIwxaSe/uZq5357C2ReEXF5lXBYgjDHGhGQBwhhjTEgWILrAl7/8Zfr378+ECRMSXRRjjImaBYgucPPNN7Ny5cpEF8MYY2JiAaILzJw5k8LCwkQXwxhjYpJe01z/+l3Ytym+9xwwES61XeKMMYnx8hMLAJh18/y43zu9AoQxxqSY/bt2+Hbv9AoQ9pe+McZEzcYgjDHGhGQBogtce+21nHvuuWzdupXi4mIee+yxRBfJGGPalV5dTAmyaFGnttI2xnSRd19cyfuvlSW6GDE58OFO+g0b7su9rQVhjDGu918r48CHOxNdjJj0GzacceeX+nJv31oQIrIQuAzYr6ptlhCLyL8A13vKMQ7o524W9CFwBGgCGlW1xK9yGmOMV79hw7nmbpvQAv62IJ4AZoc7qao/VdVJqjoJuBN4JWjXuFnueQsOxhiTAH7uKPeqiAyL8vJrAeuoN8YkRGDswc/+/GQkzq6fPt3cCRDLQ3Uxea7JASqAkYEWhIjsBGoABR5R1QURXj8fmA9QVFR0zuLFi1udz8/PZ9SoUZ38JP4rLy/n8GH/tjqsr68nNzfXt/snE6sLR7rWw+4nPwIge2oN1dvfp6mpieNVewHIHVRM4ahx9Bv/yUQWsUvNmjVrfbiemu4wi+nzwGtB3UszVLVSRPoDq0TkA1V9NdSL3eCxAKCkpESDd4basGEDeXl5/pQ8jrKzs5k8ebJv9y8rK7Nds1xWF47uXg+bV1eybU1V3O97WPejp7Zw6pX9gBMUisdPYNz5pXziorC94mmpOwSIeQR1L6lqpft1v4g8C0wDQgaI7m737t3ceOONVFVVISLMnz+fb3zjG4kuljHdRrhAsGd7LQCDRhXE9f301Baamg61BIXqrOxuHSgTKaEBQkTOAi4EvuQ5diaQoapH3MeXAPcmqIidlpWVxc9//nOmTJnCkSNHOOecc7j44osZP358ootmTLewbU0VByvq6VvcursrsGd0vLcFfXJeNYi0zFQqKyuL6/1TiZ/TXBcBpUBfEakA7gZ6AKjqw+5lc4EXVPWo56VFwLMiEijf06qatJspDBw4kIEDBwKQl5fHuHHjqKystABhjEff4lzmfntKoothgvg5i+naKK55Amc6rPfYDsCXEaIfr/kxH1R/ENd7ji0cy79O+9eorv3www/ZsGED06dPj2sZjEkWobqTQrUeTPdgK6m7SH19PVdddRUPPPAA+fn5iS6OMQkR6E7y6lucy+hpRQkqkYmkOwxSd5lo/9KPt4aGBq666iquv/56rrzyyoSUwZjuItrupJolS6lbvjzu79987BgZOTlxv28qSqsAkQiqyq233sq4ceP41re+lejiGNOlgruUou1OqlmylH133w1AztSpcS1TRk4OWX36xPWeqcoChM9ee+01fv/73zNx4kQmTZoEwH333cecOXMSWzBjfOINCsFTVQt6HafPe6+y64YHIt7j2Nq1AAy45x56X/PFuJYv+57vxvV+qcwChM9mzJiBn6vVjekKsSxa8waF4Kmqu264kRMffABjx0a8R87UqeRfdlncg4OJjQUIY0xEm1dXUvbUViC6RWveoFCzZCl1Cx5gl5ss58QHH5A9dixDf/87H0ts4sUChDEmokDLofT6MW0WrYUdSN4Iuxac7ioKjCNkjx1L/mWX+VpeEz8WIIxJM4HuotraZmrWv93u9Qcr6ulf0EDugjtbWgIBwQEgmHUVJTcLEMakmcBahKwo16YV9DpO77ee4djetW0CgQWA1GYBwpg0VNDrOKM3/pqCgoJ2r/VzRpHp3ixAGJOiws08OrCzmjMP7qTn9u0QxRoDayWkLwsQPjtx4gQzZ87k5MmTNDY2cvXVV3PPPfckulgmDYTLkpp3cj/9qtZRd/11jPv3f09Q6UwysADhs169evHSSy+Rm5tLQ0MDM2bM4NJLL+VTn/pUootm0kCotBa7bngAimHnBRckplAmaViyPp+JSMu2jg0NDTQ0NOCmMjfGmG4trVoQ++67j5Pvxzfdd69xYxnwve9FvKapqYlzzjmH8vJybrvtNkv3bbpE44EDNB461CatRWCxmjHt8a0FISILRWS/iLwX5nypiBwWkY3uv7s852aLyFYRKReRpE+ckpmZycaNG6moqGDNmjW8917IKjEmrhoPHaL52LE2x22xmomWny2IJ4CHgEhr6leraqvvVBHJBH4NXAxUAGtFZJmqbulsgdr7S99vBQUFzJo1i5UrVzJhwoSElsUkv/byI9VlFJKfA0MfDfMjaFttmnb41oJQ1VeB6g68dBpQrqo7VPUUsBi4Iq6F60IHDhygtrYWgOPHj7Nq1SrGWvPexEGozXe88purGdy4owtLZFJNoscgzhWRd4A9wHdUdTMwGNjtuaYCCNtpLyLzgfkARUVFbTYgz8/P58iRI3EudvTKy8v5+te/TlNTE83NzcydO5cLL7ywTZlOnDjh6+bp9fX1tjm7K1Xqora2maxc6H1OXcjzvcsWA4T9rKlSD7EK/MEW+OzpWg/RSGSAeBsYqqr1IjIH+DMwKtabqOoCYAFASUmJlpaWtjq/YcMG8vLyOl3Yjjr33HN555132r0uOzubyZMn+1aOsrIygusmXSVzXXi7lRrrnTUOpaWhd2fb9dhCAD4Z5rMmcz10RtUrKwFaPnu61kM0wgYIEYm4N6aq/qkzb6yqdZ7HK0Tkv0WkL1AJDPFcWuweMybteRe/2V7Oxm+RWhCfd7/2B84DXnKfzwJeBzoVIERkAFClqioi03DGQw4BtcAoERmOExjmAdd15r2MSSXR7ulsTGeFDRCqeguAiLwAjFfVve7zgTgzlCISkUVAKdBXRCqAu4Ee7r0fBq4G/kFEGoHjwDx1tl5rFJHbgeeBTGChOzZhjPEIuxeDy9Y7mM6KZgxiSCA4uKqAj7X3IlW9tp3zD+FMgw11bgWwIoqyGZO26pYvjxgE0nW9w7svruT918rCnj/w4U76DRvedQVKYtEEiP8TkeeBRe7za4AX/SuSMSZaqbp9Z3u/5COp2OIsRC0eH3qtUb9hwxl3fmkHS5Ze2g0Qqnq7iMwFZrqHFqjqs/4WyxgTijd9Rqp0IYUKBu39ko+kePwExp1fyicumh2P4qW1iAHCXdW8WVXHAhYUOqGpqYmSkhIGDx7M8gj9xsZE0pI+IyM5u5CiDQb2S757iBggVLXJzYn0MVX9qKsKlYoefPBBxo0bR11d6EVNxkQrIycnfPqMbiY4IFgwSC7RjEH0BjaLyBrgaOCgql7uW6lSTEVFBX/5y1/4/ve/zy9+8YtEF8eYuGlvrCA4IFgwSC7RBIiU2XJq9dJtHNwdPndNR/QdkssFXxwd8Zo77riDn/zkJwlN+WFMR0UKAu2NFVhASG7RDFK/0hUFSVXLly+nf//+nHPOOZbvxXRLsbYCvCwApLZ2A4SIfAr4FTAO6ImzeO2oqub7XLa4a+8vfT+89tprLFu2jBUrVnDixAnq6ur40pe+xB/+8IcuL4vp/qJK4d3ckSTJrb374kq2/uXPVL2y0loBJqxoupgewkl38QxQAtwIdP1v2iT1ox/9iB/96EeAkxTsZz/7mQUHE5Y311Io8Ujh/e6LK1n1qLNGtaCgwAKACSuqbK6qWi4imaraBDwuIhuAO/0tmjHpwdtqCASHcLmWgrcP7YhAd9LHLryYL/zjNzp9P5O6ogkQx0SkJ7BRRH4C7MXHjYZSWWlpqaUVTkPtdRvt2V4LwKBRBSEztHpzLsVrcVzx+An0G//JTt/HpLZoAsQNOAHhduCbOKm4r/KzUMakkva6jfoXNFBUtZahG7c5BzbCrgWnzx9buxaAnKlTY1ocF27w2XIRpZi/ftf5eun9cb91NAFiJvBnd/+GewBE5DKgPO6lMSaFBFoOB3ZWk3dyP1M2rgx5XSAAMHVqyPM5U6eSf9ll9L7mizG9//uvlYUMBoFcRJ0f6k5D6x6HTX9MdCla2/U3GDrDl1tHEyB+BXxbRK5V1ffdY/cCli/CmAi2/OU9qmuU3MMf0q9qnbP1VQgdDQDR6DdsONfcHfovS5t23QGb/gj7NsGAiYkuyWlDZ8DEq325dTQBYidwK/BHEfmBqj4DiC+lMSaFNB46RO6xY8zo8Tr5X/cnAJgEGDARbvlLokvRJaIJEKqqb4vIhcAiEZmOsxYiIhFZCFwG7FfVNhOsReR64F9xgs0R4B9U9R333IfusSagUVVLovw8xnQrfuZNsn0PjN+imY20F0BVDwKfBRSIJgfvE0CkidU7gQtVdSLwH8CCoPOzVHWSBQdjQguMMYRj+x6Yzoom1cbnPI+bgX9x/7X3uldFZFiE8697nr5J2B7a5Dds2DDy8vLIzMwkKyuLdevWJbpIJkl5Ww2BFkK4MQbTAe0NQne38QefRZNqox9OV9B4IDtwXFU/Hcdy3Ar81fNcgRdERIFHVDW4dZF0Xn75Zfr27ZvoYpgk552ZZC0EH7Q3CD1gom8Dwt1RNGMQTwFLgM8BXwduAg7EqwAiMgsnQHjnac1Q1UoR6Q+sEpEPVPXVMK+fD8wHKCoqajMzIz8/P+FZVFWV+vp6evXqFfaaEydO+DqrpL6+3matuLqqLhobG4H4zhaqra2lR0Fvii50em+rO3F/+55weOthUm0tZA9h4/AInST1QJrUWzQBoo+qPiYi33Azu74iImvj8eYi8gngt8ClqnoocFxVK92v+0XkWWAaEDJAuK2LBQAlJSUavFJ5w4YN5OXlAfDyEwvYv6tzeWyC9R86glk3z494TUZGBldeeSUiwte+9jXmz297fXZ2NpMnT45r2bzKyspsFberq+pi0VNPA8T1vapeWRm3e9r3hKNVPewsAOL7f5bMogkQDe7XvSLyOWAPUNjZNxaRjwF/Am5Q1W2e42cCGap6xH18Cc66i6T1t7/9jcGDB7N//34uvvhixo4dy8yZM9t/oUkaodJpxCvzqjGJEk2A+KGInAV8G2fRXD5Oyo2IRGQRUAr0FZEK4G6gB4CqPgzcBfQB/ltE4PR01iLgWfdYFvC0qoZeghqj9v7S98vgwYMB6N+/P3PnzmXNmjUWIFJMqHQa8ci8akwiRTOLKbBi+jAwK9obq+q17Zz/CvCVEMd3ACmTRezo0aM0NzeTl5fH0aNHeeGFF7jrrrsSXSzjg0AW1kByvXgl1jMmUaKdxfRVYJj3elX9sn/FSh1VVVXMnTsXcAYtr7vuOmbPtrz7qabxwAEaDx1i1w0PtEquF21iPWO6o2i6mJ4DVgMv4qxsNjEYMWIE77zzTqKLYXzgHXeorlFyjx2DDH9zKxnTlaIJEDmq+q++l8SYJOMdd8hvrmZwz30MfdyftBrGJEI0AWK5iMxR1RW+l8aYbibcZj+NBw5QXaPkN1czZeNK38YbgvMtWX4l05XCBggROYKzolmA74nISZwpr4KTwC+/a4rYeaqKOyuqW1LVRBchrUXa8c2725tXIFPr4J77AGLayCdYpKR7FVveA5wd4MDyK5muFTZAqGpeVxbEL01NTRw6dIg+ffp0yyChqhw6dIjs7Oz2Lza+iLTjW5vd3lyBFkM0XUrtZV0NDgJexeMnMO78Uj5xkU1sMF0vmi6mpHb06FGOHDnCgQNxyw4Sd9nZ2RQXp2yuwm6v8cABcusPhdzxLdxub9G0GAKBIVIACBy3IGC6o5QPEKrK8OHWZ2tC27y6kv21PShwZyAF68yMpEBiPQsAJlmlfIAw6ck7rhBYoxDQ2NjYkiepOnMAAEN8moFk6bhNMosqQIjIDGCUqj7uLpzLVdXwO5UYk2CB/aDzm6tpPlIHQEZe23kVhU37GNy4g4mzR3d1EY3p9qJZSX03UAKMAR7Hyaf0B+B8f4tmTMcFZhmdl/Eq9KBVN5FlMTUmOtG0IOYCk4G3AVR1j4ikxAwnk9r83A86WKiZSrZmwSS7aPakPqXORH2FlnTcxhiPUPtD25oFk+yiaUEsFZFHgAIR+SrwZeBRf4tlTPKxAWmTaqJJ9/0zEbkYqMMZh7hLVVf5XjJjjDEJFc0g9beAJRYUjGnNO+5g4w1JZN3jsOmPLU8n1da2bDXKvk0wYGJCitUdRdPFlAe8ICLVwBLgGVUNnbgmiIgsBC4D9qtqm2Wk4uS+eBCYAxwDblbVt91zNwH/5l76Q1V9Mpr3NOktsP7B7+0+331xJasefQhwVkLbeEMS2fTH8IFgwESYeHXYlz791kc8t7HSx8J1zPhB+dz9+bPjft9oupjuAe4RkU8A1wCviEiFql4Uxf2fAB4Cwk0luRQY5f6bDvwGmC4ihThblJbgDI6vF5FlqloTxXuaNBZY/5BbX9GSSK8zwuVRCqTPuPirt9sK6e4mqIXQRiA43PIXADbGMO35iU2L2Kuvk9Ore60xPnX048BP4n7fWD7lfmAfcAjoH80LVPVVERkW4ZIrgN+5s6TeFJECERmIs5f1KlWtBhCRVcBsYFEM5TVpKLD+YUbGq+TPji27aqhgEC6PkqXP6MYitRCg3VZCJIcz15CRtZfxA+L/13pnjC3s48t9oxmD+Efgi0A/4Bngq6q6JU7vPxjY7Xle4R4LdzxU+eYD8wGKioooKytrdb6+vr7NsXSULvXQ2NgIPXuy8yZ3R9wQnzm4Lg5seYfq7e9Tv6cCgNxBpxMn5g4qpnDUOPqNb7tNejUkdZ125++JgXuep6jq1Q69Nrd+J/W5w9k4/F/CX1RPy/dGLPXQ2NhIFgO5KfumDpXNN8f8+V6MpgUxBLhDVTfG/d3jQFUXAAsASkpKNLipaKtmHelSD4EcS5E+a3BdLHllJQ21NWnXKug23xOhuoR2/c35OnRG7PcrmEzBxKspLSmN6vJY6iFr56+AyN9fqSTShkH5qloH/NR9Xug9H+j+6aRKnAAUUOweq8TpZvIeL4vD+5kUEGmDn44OTtsahgQK1SU0dIbTDVRyS5cWpb1B6GPa2O3GH/wU6ZM+jTMDaT2nd5YLUGBEHN5/GXC7iCzGGaQ+rKp7ReR54D4R6e1edwlwZxzezyS5zasrKXtqK9B2lzfA2Ru6cUfU93v3xZVUbHkv7F4Npot4Bo0T6bmNlWzZW8f4gaE3zMzplUXfM3t1cakSJ9KOcpe5Xzs8uVtEFuG0BPqKSAXOzKQe7n0fBlbgTHEtx5nmeot7rlpE/gNwd2vh3ji1WEySC7QcSq8fw9kXtB2W2nXDA22OBQ8+19bWUvWKszlQYBDapqjGWXsziby62dqD8QPzWfK1c0Oeu2Vl0uy0HBfRDFL/n6p+pr1joajqte2cV+C2MOcWAgvbew+THgLdSgcr6ulf0EDugjvZtaDtdYGtQL1Bwbb0jINYfuFDbGMInZhVZPwVaQwiG8jB+eu/N6e7mPIJM6PImHgINcawZ3st4OwR3futZzi2dy05QduAwumtQN90k+f1Gza8TRDoNoOzyaS9qaPBEjSGYOIrUgvia8AdwCCccYhAgKjDWfxmjC8CLYW+xbktxwaNKmD0tCJyF9zJsb1rGXDPPZG3Ab3nbRt4jrduMk5guk6kMYgHgQdF5J9U9VddWCaTZoJbDIHgMPfbU1qO1SxZSt2CBzjxwQfkTJ3aoT2iTQete9zpMurIlFOT1KJJtfErEZkAjAeyPce7ZicWk/K824MC5AJ93tvRasD52FpnvkLO1KnkXxbbCmnTSYGxhzQZJ6jJfJVbVoYY4AK2Vm9lTOGYLi5R4kS75WgpToBYgZM/6W+Ez69kTExabQ8aRiAwtNdyCAxOW3bVOBs6I2nGE2JNqFdbe5zfbH0DgC31L0DfP7K3CkqKStpcO6ZwDHNGzIlbWbu7aFZ8XA18EtigqreISBHOntTGxE28tgf1BgebukrE2Uet0lxH4vM01HhnSH1rp9MSnT68sJ0r28opfJdjwF3n3sUXRn8hbmVKVtEEiOOq2iwijSKSj5O0b0h7LzImkpolS6lbvhyA5uaZZOTkxO3eNjjtEevso1B8noba3uK0WE0fXsgVkwZz3fSPRXW9M6vNWffgrHMoseDgiiZArBORApxtRtfjpLl6w89CmdRXt3x5y5qFjJwcsvr4k43SEHb2USxprmMVS6sgEBzCLU4ziRPNIPU/ug8fFpGVQL6qvutvsUwqq1mylO0VPTkw6Q6yx46lvqKevv1y23+hSahYfunH0s0zfmA+V0yypVXdUaSFclMinQvs/GZMtALdStsrerJ1zHWAs8imb3Euo6cVJbZwqShO01MDgSGWX/qxdvMk2jPbnmHFjhVpN0upPZFaED+PcE6BT8e5LCbFbVq5jd3NM6kd4+R5DJdPycRJnKanBsYIku2XfjiBYBBQW1tL+a5ywJm5lE6zlNoTaaHcrK4siEltNUuWsvt4X+rPKm5ZFd3R4BBuG1Agvaa3RrO15tAZPN30GZ57pO2woXd6ZyTJNkYQHACCrataB7SexhoIDDY43Vo06yBuDHXcFsqZWGxauY3agvPoX9DQaoV0R0Ra55By01sjBYH2EuINmMhbuZ/me89uAjo27ROSb4ygva6i4GBgubnCi2YWkzcjWjbwGeBtbKGciUFlltOtNP5zndt3wbt/Q1pMZY00TTWKhHi/eOQNoJr75k5s0zXknd6ZasYUjuHx2Y8nuhhJL5pZTP/kfe5OeV3sV4FM6ips2sfZF0Q3dBWuGyml928I1VoIBIdOJMmbPrww6cYN2usmisQGmuOnI3vnHQWi6uQVkdnAg0Am8FtVvT/o/C+BwFhHDtBfVQvcc03AJvfcR6p6eQfKarqBmiVLaT5SR0Ze6IVQoYJBuD0cUnL/hkBgCNVlFMMitVDTUOO5AC0WnfkFD6HHCaKVbukw/BTNGMT/4sxaAsjAycm0NIrXZQK/Bi4GKoC1IrJMVbcErlHVb3qu/ydgsucWx1V1UhSfwXRjNUuWsv7hF6gdcx39CxpCXhNqTCElAwGEbiV4A0MMeygEB4RQ01ATNX7Q2SmjNmjcPUTTgviZ53EjsEtVK6J43TSgXFV3ALj7Tl8BbAlz/bU4W5KaFLJp5baWNQ+Rxh9SOj2GNyiEaiV0MDAEB4RETEMN11IIBAcbB0hu4uz6GcWFTh6mloDS3h7RInI1MFtVv+I+vwGYrqq3h7h2KPAmUKyqTe6xRmAjTlC6X1X/HOZ95gPzAYqKis5ZvLj18Eh9fT25ubZKNxH1cMbq1ZRv7kdtwWgGlgiFI6Xl3IEt71C9/X0Ajh88wBl9+zHminldUq6uqouBe56nqOpVCg473WW1ZzkBsqpoJnsHfTaqe5TtbuCNPY2tjm2taQZgTO8Mzh2URemQHh0qXyz18NqR11h3dF2b4+UnnfUDI3uNbHOu5MwSzs87v0Nl60rp/jti1qxZ61U1ZF9eNF1M84F7gRNAM87OcgqMiGMZ5wF/DAQH11BVrRSREcBLIrJJVf8e/EJVXQAsACgpKdHg6Wo2hc2RiHrY9dhCyulH/4IGrvzKZ8PuE11QUOB0J3VR+Xyvi1BjChOvpsBtIRQA0XS8PP3WRzyxue0U1ekFxKWlEKiHaMYL1lWHHhMoIfm7gux3RHjRdDH9CzBBVQ/GeO9KWmd9LXaPhTIPuM17QFUr3a87RKQMZ3yiTYAwiefNzOpVfjCf2qGjGdSvgHdfXMmqR52daovHT0jdMQY4PTW1A/sye8cVAl1IoaaoRivSL//a2lqeXPlkVAPCNiaQnqIJEH8HjnXg3muBUSIyHCcwzAOuC75IRMYCvfFkiBWR3sAxVT0pIn2B84GfdKAMxgfBAcG72xvArqzRVGaNoHroAABGTyvivZeceQ0Xf/X21AsKwQPPHZya+vRbH7Va1BaPMYVoBovtl78JJ5oAcSfwuoi8BZwMHFTVf470IlVtFJHbgedxprkuVNXNInIvsE5Vl7mXzgMWa+vBkHHAIyLSjDNz6n7v7CeTWN5U3dB6t7fNqyvZ9NRWgFYpNd57yWk5pFxwgLaL2Tq4f0Kg5dCZFkMo4QaLrWvFtCeaAPEI8BLOmoTmWG6uqitwtin1Hrsr6PkPQrzudcC/LaxMp2WPHcvQ3zuL6TevruTtNVXw87fZs70WaJ2Iz7v6OaUEWg6dXMwW6FYKJMRLtkVtJnVFEyB6qOq3fC+J6fYCXUvlB/PZX3web//cyfgeCAqDRhW0ScTnHXtIqdXP6x6H5Xc4jwNjDVEItZjNO121M2sWQo032Kpi0xnRBIi/ujOZ/pfWXUwRp7ma1BPoWto/6Q6O9OpPtns8UnbWwKylpB57iLS47bIH2gxCR9pYJ9RitnitXwg13mCrik1nRBMgrnW/3uk5Fu9prqabq1mylGNr15IzdSrZY8eSDVFnZU3qsYfglkJAhBlKkfZYjkcwsMVppqtEk6wvTZLrm0jqli+ncuD51Az+ArUV9fQtbn9hUUqMPQRaDiFaCpHEa/+EUMEg3LRUay2YeLP9IEzUDoz8NPUnz4h6i9BA91LSjz0MnRFTcIinUN1GNi3VdBXbD8JEZVfWaKozBzCoOLfdrqXAiukDH+5MXPdShI12JtXWws6C6O4Tbi+GELyzkTqSQTXSILN1G5lEsP0gTFjeBXG7m2dCL9ptOQSvmE5Y6yHSRjuxiGJNQ6jkeR2ZjWSDzKa78XU/CJPcvAviMnJy6F/QwNkXDI64J3Qgx1K3mLUUZm3CxjgvEPOuYejIAHSg5WCtBdPd+LYfhEluNUuWsr2iJwcm3UH22LHUV9TTt58zMB1pT+iUzrHk4Z3KGuhSimZQur1BZ2stmO7Ez/0gTBJa+8vnKN9ylOYjddS6+zgMAvoW53JGzlaW3LO0JTik7P4NQdpb3BbLpjw26GySSdgAISIjgSJVfSXo+Pki0itU6m2T3Nb+8jnWbM2DzDwK86B/QQPjPzehZQGcNzh025lJwekvOik4gV5AZ9YzWDeSSRaRWhAP0HpxXECde+7zPpTHJFD5lqOQmce0MUeY+s02iXeBbrLzW4QZSm2276RtC6C29ji/2fpGqFe30ZGU25FSbFvqC5NMIgWIIlXdFHxQVTeJyDD/imQSqbBpX8jgkNBFb8EBIdS2nQFuYHi66TM8t74S1r8RMr1FtGJpKQQCQ6T9FWxWkkkmkQJEQYRzZ8S5HCbBNq+upDpzAIVN+0KeT+iit+Auoyg24nnukTdaBo+Df8k7aa47v8o5IFRgsDEFkwoiBYh1IvJVVX3Ue1BEvgKs97dYpqsEBqWrM53NfQY37gh7bUIWva173GkxDJ0RczrteKW7aE9g4NkCg0k1kQLEHcCzInI9pwNCCdATmBvNzUVkNvAgzoZBv1XV+4PO3wz8lNNbkT6kqr91z90E/Jt7/Ieq+mQ072mi12pQumkfgxt3MHH26Jbz3vUO4aa1+sLbpRToTurABjxdyQaeTSoKGyBUtQo4T0RmAYGO57+o6kvR3FhEMoFfAxcDFcBaEVkWYme4Jap6e9BrC4G7cQKSAuvd19ZE894mOu0NSnvXO3TZzKXg7Kkx7uvc2XQXoUQadAYbeDapK5pUGy8DL3fg3tOAclXdASAii4ErgGi2Dv0ssCqw54SIrAJmA4s6UA4TQs2SpTQfqaMwjzbBwZtLqUtmLYVqMYTJnhpprwWI3+Y7Xu3t62wDzyZVdSTVRrQGA7s9zyuA6SGuu0pEZgLbgG+q6u4wrw350+5uZjQfoKioiLKyslbn6+vr2xxLR8H1cPJPG6ntfRG5PQ+3HD+w5R2qt79P/R5nHWTuoGIy+g2MS/0N3PM8RVWvhjxXcNhJz1F71gQ4awJVRTPZWz8cPO9btruBN/Y0srXG2fV2TO+MkPca0zuDcwdlUTrkJBzfQVlZ2zGVWL4nXjvyGuuq1zGy10huyr4p/IV7oGxPdPfsLuxnw2H1EJ6fASIa/wssUtWTIvI14Eng07HcQFUXAAsASkpKNDjHjm3M7giuh0VP7QGg5AvTTi+Ee2UlDbU18U2XEWgdRJqaWuB0IxW4LYYCwPu3+tNvfcQTm08vVotm2mmkbqHa+loKsgqiK361MzPpuinXUTq6NKrXJAv72XBYPYTnZ4CoBIZ4nhdzejAaAFU95Hn6W+AnnteWBr22LO4lTFPe7qWzL2gdj+PepRSYohrjWEKAdyVzLIvV2usWipbNTDLpzM8AsRYYJSLDcX7hzwNadXaLyEBV3es+vRx43338PHCfiPR2n19C6FXdJkpnrF7NrscWArC9oie1Y66jf0FDfGYqRVrZHFi/EMMUVe84Q0dWMgeEm1lkfzEaEx3fAoSqNorI7Ti/7DOBhaq6WUTuBdap6jLgn0XkcpwkgNXAze5rq0XkP3CCDMC9gQFrE7uaJUs58tIuyovOIyMvn+oxzpqHvKHVrHr0KcBZ4xD1TKVYVjZHsZ+CV3Duo2i6lCJttGOM6ThfxyBUdQWwIujYXZ7HdxKmZaCqC4GFfpYvXdQtX05V0Xkc7TuSfsMLGYSz8c97L70AtLN3Q6jWQXBA6OBU1FDCtRgijSmESm1hM4uM6bxED1KbLrArazS1Z45m0PCCVtuFvvdShNXRkQaXOzieEBBpnUK4FkOkMQUbJzDGHxYgUlzNkqXsPt63zXahIZPvhVqP0MlgENDRDXbAaT2sq1pHSVGJrVY2pgtZgEhB3r2kAwPSuT0Pt0xn9e4b3WrMwZsULw6BIdRgc6wb7AAtXUvWZWRM17IAkWJqlixl/cMvUBU0IK2FH7Hknu8CIfaNDt5kp4MzjoJ5g0JnNtgBpxvJupCM6VoWIFJIIDhsDWwVOqqAQcAZOVvZXLYKcMYcWi2EC859FGbGUbhAEGmvhc4GBWjdvWSM6VoWIFLIppXbWoJD6fVjWm0VCiFmK3mDQzu5j8IFgngEgUise8mYxLEAkUIqs0YArYNDYDA6d1Bxh4KDd02Cn4EgIHg6a2CfBeteMqbrWYBIMYVN+2g6eYol9/wKOD3eUDhqXOsLA7OVPMEhuBupM6uYOyp4OqutZzAmcSxApJCjp8o53vAhqx7dD7Qeb6jOym77gqEzWgUHb2sh8LUrWg3BbPMdY7oHCxApILBt6JGGCrQpdDbWcOmMg8cYuqK1EGlVtKXIMKb7sACR5LzbhmbRyBl5fdvPxurZ5zmwqrkrWguBwBAqNUaAdSkZ031YgEhy3m1Dd9T1D3vdwD3Pw+M/dZ64q6QX1E5hy+HYVjXHytta8AYGS41hTPdnASJJBbqVDp2sRBtfZEdd/4jpuouqXoUTu1tWSS+oncKvDs+IeVVzsPb2a/YGBQsMxiQXCxBJqnzLUer0LLTxRZqoBfqHT9e97nFnW8+hM1pWSf/fI28wPodOtxza25jHgoIxycsCRJI6eqqcUw0fIr2OMnDY6PDjDt71DjHsyxDQXgshEBxs1pExqcfXACEis4EHcTYM+q2q3h90/lvAV3A2DDoAfFlVd7nnmoBN7qUfqerlfpY1maz95XPUNe6BphoGDxsTeZMfd73DAz2+whvrx8L6NwDCptsO1l4LwQaVjUldvgUIEckEfg1cDFQAa0Vkmapu8Vy2AShR1WMi8g84e1Jf4547rqqT/CpfMgpkad3WcB5kQn6kGUtuAr5Tle/wdvM4Hjjyaab3PX06lrEHayEYk578bEFMA8pVdQeAiCwGrgBaAoSqvuy5/k3gSz6WJ+ltWrmN3c0zqcmuRY9XkDdkQugLPd1KbzeP47mm87j57J784Ib2xxtCpbqwdQnGpCc/A8RgYLfneQUwPcL1twJ/9TzPFpF1ON1P96vqn+NewiQRmLG0n1yaWIcerwBo27UUtAvcPcxny+AruWLSYAYd39Hu+zyz7RnufeNe4PQaBetCMiZ9iar6c2ORq4HZqvoV9/kNwHRVvT3EtV8CbgcuVNWT7rHBqlopIiOAl4DPqOrfQ7x2PjAfoKio6JzFixe3Ol9fX09ubm58P1wXOWP1amrfO05574sAaDr8e5qoIWfAAApHjaPf+E+2un7Shu+TW7+T+tzhPHHsU/w16yLunH4GEL4eXjvyGuuOOlNRy0+WAzCvcB7n553v50dLqGT+nognqwdHutfDrFmz1qtqyHz6frYgKoEhnufF7rFWROQi4Pt4ggOAqla6X3eISBkwGWgTIFR1AbAAoKSkREtLS1udLysrI/hYsnj1yb+1BAdnIVwRUBR+3GFnARRMZsX43/DAs5uYPryA0lKnWym4HlpWNVd71imQHlNSk/l7Ip6sHhxWD+H5GSDWAqNEZDhOYJgHXOe9QEQmA4/gtDT2e473Bo6p6kkR6QucjzOAnVYC6bs//skqdux/J+JCOICqIyc4WH+S7211Jn9FGoQOzE6ydQrGmHB8CxCq2igitwPP40xzXaiqm0XkXmCdqi4DfgrkAs+ICJyezjoOeEREmoEMnDGILSHfKIUdPVVOc8OLbC47nZ013JTWt575OdOr17GzeVzUeZVsdpIxJhJf10Go6gpgRdCxuzyPLwrzuteBiX6WrTsLTGc9dkpopC5kdlZwMrHWv/4o5x9/memnnFaDTPwCS77QdrbSa0de48mVT7Y8t9lJxpj22ErqbqZmyVL23X03ADLlXHr26E/T7H/kPzdWwt/faLnuM8dW8ImaVXwq430ANvecSP2ouUz/wrdb3S/UWAPY7CRjTPssQCRYoLUQcGztWgDqv3Y/Dev+imRksGnZA3wz83Xysk//d519ahNkQFVhCUXnfYmzQ2wZCqfHGkb2Gsl1U66zsQZjTNQsQCSANygEAkLO1KktX/dPucrZ4wFA6/lRj6ecxwNneO4yAyZeTVGYwBBoOQS6km7KvonS0aV+fBxjTIqyANHFvF1IOVOnsm/KRPb0ziOrfz8Ajtaeonads14wM+MQeRlHnBd69o5uT/CCtzkj5sCe+H4OY0zqswDRxeqWL+fd4ZM50D8P6dGDE0eq4eAxsk86C9pOHm8EoOCsBno01DM6p5LNPSeG7ULyCt6x7a5z72rpUirbU+bPBzLGpCwLEF1sLWdRkb8TTtSR3WMY2XnDyO3zCfL7O4PHdQd2MyjjJS7uuQiAN5vHUT9qbsR7htrK09Y2GGM6ywJEF6pZspQ9pw4DcHbp9cz+h2vbXLP5vn9myKm/s7nnRF47YxbbxxdzUN/k4ZXhWxAWGIwxfrAA0YUCA9O9sge3BIen3/qI5zaezkDynVNN7O75cc7+3t/Ysu0ZVgYlzwvFAoMxxg8WIHwUPIX1zcNCwxmHyO4xrOVY/euP8p3DL5LTMxOAjXn7WFHQm54rbwk5lmCMMV3FAoQPAoHhg79/wJ7eeWTkOVNWD57hDECf1a+Qzfc5U1bnuyugA1NYfyYn2JkJY7CWgTEmsSxAxEGoxW4fFebx3pD+AGTnFgOQcaKRk7nj2FP4O17IO8qJjGxgOE1n9CW/n3Pt1uoay5FkjOkWLEDEQd3y5Zz44AOyx44FYPPUi9h1aicAWTkXUVswybnueAN52Zv5e14d27LPYNyAyW3uZSkwjDHdhQWITgi0HN48LOwfNYGMXjkAnDjiBIdeAy6lNquBi3K+S152Fi/mnGJjdg1be/Zk3JmDrZVgjOnWLEDEyNud9ObBRqrys2g44xA0QTbDAMjOG0bGsMl8dHwzUwoX88iZZ0L2WayTk0A2JTnFzPnklxP3IYwxJgoWIGIU6E6qmnA5lYUVaNMBsnOdgFDWb3zLdZmHf0mvvhu494w+AJQUTaQEbNDZGJM0LEBEEDz4DFB+MJ/9k+5gf+N+tLGC3oNGkz3vDpa/8h0a9b/IzBAANhU14W0tWFAwxiQbXwOEiMwGHsTZUe63qnp/0PlewO+Ac4BDwDWq+qF77k7gVqAJ+GdVfd7Psgas/eVzlG85CkB93SaOA5J5upoaChvh2KtoYwUA6/PfYt+Gz7OtqBmAEu3lfs1izqAL+MIlv+yKYhtjTNz5FiBEJBP4NXAxUAGsFZFlQVuH3grUqOpIEZkH/Bi4RkTG4+xhfTYwCHhRREarapNf5Q14Z+M7HGmoIItGTuFs9dkjo2/L+SxO0DOznv25x9g66BgNA0+QifCJxh78vyEXWkAwxqQMP1sQ04ByVd0BICKLgSsAb4C4AviB+/iPwEPibE59BbBYVU8CO0Wk3L3fG/jgoetvpbHJWcTWpIcAOFhwkkyaqC2qo6K4os1rdvVoojiriP+54U0/imSMMQnnZ4AYDOz2PK8Apoe7RlUbReQw0Mc9/mbQaweHehMRmQ/MBygqKqKsrKzV+fr6+jbHgqnncab0oemMRrZMeQ+Akz0KyMwpbPOaQcCUM0vavXd3EU09pAurC4fVg8PqIbykH6RW1QXAAoCSkhItLS1tdb6srIzgY8HaO58KoqmHdGF14bB6cFg9hJfh470rgSGe58XusZDXiEgWcBbOYHU0rzXGGOMjPwPEWmCUiAwXkZ44g87Lgq5ZBtzkPr4aeElV1T0+T0R6ichwYBSwxseyGmOMCeJbF5M7pnA78DzONNeFqrpZRO4F1qnqMuAx4PfuIHQ1ThDBvW4pzoB2I3BbV8xgMsYYc5qvYxCqugJYEXTsLs/jE0DIFWSq+p/Af/pZPmOMMeH52cVkjDEmiVmAMMYYE5IFCGOMMSFZgDDGGBOSOLNKU4OIHAB2BR3uCxxMQHG6G6uH06wuHFYPjnSvh6Gq2i/UiZQKEKGIyDpVLUl0ORLN6uE0qwuH1YPD6iE862IyxhgTkgUIY4wxIaVDgFiQ6AJ0E1YPp1ldOKweHFYPYaT8GIQxxpiOSYcWhDHGmA6wAGGMMSaklA4QIjJbRLaKSLmIfDfR5fGbiHwoIptEZKOIrHOPFYrIKhHZ7n7t7R4XEfkvt27eFZEpiS19x4nIQhHZLyLveY7F/LlF5Cb3+u0iclOo9+rOwtTDD0Sk0v2e2Cgiczzn7nTrYauIfNZzPKl/bkRkiIi8LCJbRGSziHzDPZ523xOdpqop+Q8nxfjfgRFAT+AdYHyiy+XzZ/4Q6Bt07CfAd93H3wV+7D6eA/wVEOBTwFuJLn8nPvdMYArwXkc/N1AI7HC/9nYf9070Z4tDPfwA+E6Ia8e7PxO9gOHuz0pmKvzcAAOBKe7jPGCb+3nT7nuis/9SuQUxDShX1R2qegpYDFyR4DIlwhXAk+7jJ4H/5zn+O3W8CRSIyMAElK/TVPVVnP1EvGL93J8FVqlqtarWAKuA2b4XPo7C1EM4VwCLVfWkqu4EynF+ZpL+50ZV96rq2+7jI8D7OHvap933RGelcoAYDOz2PK9wj6UyBV4QkfUiMt89VqSqe93H+4Ai93Gq10+snzuV6+N2t+tkYaBbhTSpBxEZBkwG3sK+J2KWygEiHc1Q1SnApcBtIjLTe1KddnPazWtO18/t+g3wcWASsBf4eUJL04VEJBf4H+AOVa3znkvz74mopXKAqASGeJ4Xu8dSlqpWul/3A8/idBdUBbqO3K/73ctTvX5i/dwpWR+qWqWqTaraDDyK8z0BKV4PItIDJzg8pap/cg/b90SMUjlArAVGichwEemJs9/1sgSXyTcicqaI5AUeA5cA7+F85sDsi5uA59zHy4Ab3RkcnwIOe5rfqSDWz/08cImI9Ha7YS5xjyW1oHGluTjfE+DUwzwR6SUiw4FRwBpS4OdGRARnv/v3VfUXnlP2PRGrRI+S+/kPZ3bCNpxZGd9PdHl8/qwjcGacvANsDnxeoA/wf8B24EWg0D0uwK/dutkElCT6M3Tisy/C6T5pwOknvrUjnxv4Ms5gbTlwS6I/V5zq4ffu53wX5xfhQM/133frYStwqed4Uv/cADNwuo/eBTa6/+ak4/dEZ/9Zqg1jjDEhpXIXkzHGmE6wAGGMMSYkCxDGGGNCsgBhjDEmJAsQxhhjQrIAYdKCiPTxZDTd58lwWi8i/+3D+31dRG7sxOufEJGr41kmY2KVlegCGNMVVPUQTroJROQHQL2q/szH93vYr3sb01WsBWHSmoiUishy9/EPRORJEVktIrtE5EoR+Yk4e2ysdNM3ICLniMgrblLE50NlwXXv9R33cZmI/FhE1ojINhG5IMT1IiIPufswvAj095y7S0TWish7IrLAvfbjIvK255pRgecicr84eyG8KyK+BUGT+ixAGNPax4FPA5cDfwBeVtWJwHHgc26Q+BVwtaqeAywE/jOK+2ap6jTgDuDuEOfnAmNw9i24ETjPc+4hVZ2qqhOAM4DLVPXvwGERmeRecwvwuIj0ce91tqp+Avhh1J/cmCAWIIxp7a+q2oCTciETWOke3wQMw/klPgFYJSIbgX/DSeLWnkDCuPXufYLNBBapk1hvD/CS59wsEXlLRDbhBK+z3eO/BW4RkUzgGuBp4DBwAnhMRK4EjkVRNmNCsjEIY1o7CaCqzSLSoKdz0TTj/LwIsFlVz+3IfYEmYvi5E5Fs4L9x8gPtdsdPst3T/4PTGnkJWO+OsyAi04DPAFcDt+MEFWNiZi0IY2KzFegnIueCk1ZaRM5u5zXReBW4RkQy3TGNWe7xQDA46O5v0DKzSVVP4GQX/Q3wuFueXOAsVV0BfBP4ZBzKZtKUtSCMiYGqnnKnn/6XiJyF8zP0AE4G3c54Fucv/S3AR8Ab7vvVisijOGm69+Gk4/Z6CmfM4QX3eR7wnNvyEOBbnSyXSWOWzdWYJObOlDpLVf890WUxqcdaEMYkKRF5ltOzroyJO2tBGGOMCckGqY0xxoRkAcIYY0xIFiCMMcaEZAHCGGNMSBYgjDHGhPT/AQNi3SyELLY5AAAAAElFTkSuQmCC\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.event_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": {
    "collapsed": false
   },
   "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",
    "[ELI5](https://eli5.readthedocs.io/en/latest/overview.html) library,\n",
    "which is fully compatible with scikit-survival."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <style>\n",
       "    table.eli5-weights tr:hover {\n",
       "        filter: brightness(85%);\n",
       "    }\n",
       "</style>\n",
       "\n",
       "\n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "        <table class=\"eli5-weights eli5-feature-importances\" style=\"border-collapse: collapse; border: none; margin-top: 0em; table-layout: auto;\">\n",
       "    <thead>\n",
       "    <tr style=\"border: none;\">\n",
       "        <th style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">Weight</th>\n",
       "        <th style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">Feature</th>\n",
       "    </tr>\n",
       "    </thead>\n",
       "    <tbody>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 80.00%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0676\n",
       "                \n",
       "                    &plusmn; 0.0229\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                pnodes\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 91.29%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0206\n",
       "                \n",
       "                    &plusmn; 0.0139\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                age\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 92.19%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0177\n",
       "                \n",
       "                    &plusmn; 0.0468\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                progrec\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 95.29%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0086\n",
       "                \n",
       "                    &plusmn; 0.0098\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                horTh\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 97.61%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0032\n",
       "                \n",
       "                    &plusmn; 0.0198\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                tsize\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(120, 100.00%, 97.63%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                0.0032\n",
       "                \n",
       "                    &plusmn; 0.0060\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                tgrade\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(0, 100.00%, 99.21%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                -0.0007\n",
       "                \n",
       "                    &plusmn; 0.0018\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                menostat\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "        <tr style=\"background-color: hsl(0, 100.00%, 96.19%); border: none;\">\n",
       "            <td style=\"padding: 0 1em 0 0.5em; text-align: right; border: none;\">\n",
       "                -0.0063\n",
       "                \n",
       "                    &plusmn; 0.0207\n",
       "                \n",
       "            </td>\n",
       "            <td style=\"padding: 0 0.5em 0 0.5em; text-align: left; border: none;\">\n",
       "                estrec\n",
       "            </td>\n",
       "        </tr>\n",
       "    \n",
       "    \n",
       "    </tbody>\n",
       "</table>\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "    \n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import eli5\n",
    "from eli5.sklearn import PermutationImportance\n",
    "\n",
    "perm = PermutationImportance(rsf, n_iter=15, random_state=random_state)\n",
    "perm.fit(X_test, y_test)\n",
    "eli5.show_weights(perm, feature_names=feature_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "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.0676 points.\n",
    "Again, this agrees with the results from the original\n",
    "[Random Survival Forests paper](https://projecteuclid.org/euclid.aoas/1223908043)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import shap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "The passed model is not callable and cannot be analyzed directly with the given masker! Model: RandomSurvivalForest(max_features='sqrt', min_samples_leaf=15,\n                     min_samples_split=10, n_estimators=1000, n_jobs=-1,\n                     random_state=20)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_52797/2067042318.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mexpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshap\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mExplainer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrsf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/foo/venv/bar/lib64/python3.9/site-packages/shap/explainers/_explainer.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, masker, link, algorithm, output_names, feature_names, **kwargs)\u001b[0m\n\u001b[1;32m    166\u001b[0m                 \u001b[0;31m# if we get here then we don't know how to handle what was given to us\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    167\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 168\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"The passed model is not callable and cannot be analyzed directly with the given masker! Model: \"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    169\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    170\u001b[0m             \u001b[0;31m# build the right subclass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: The passed model is not callable and cannot be analyzed directly with the given masker! Model: RandomSurvivalForest(max_features='sqrt', min_samples_leaf=15,\n                     min_samples_split=10, n_estimators=1000, n_jobs=-1,\n                     random_state=20)"
     ]
    }
   ],
   "source": [
    "expl = shap.Explainer(rsf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "Model type not yet supported by TreeExplainer: <class 'sksurv.ensemble.forest.RandomSurvivalForest'>",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_52797/1379721347.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mexpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshap\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTreeExplainer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrsf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/foo/venv/bar/lib64/python3.9/site-packages/shap/explainers/_tree.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, data, model_output, feature_perturbation, feature_names, **deprecated_options)\u001b[0m\n\u001b[1;32m    145\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfeature_perturbation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfeature_perturbation\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    146\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpected_value\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 147\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTreeEnsemble\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata_missing\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel_output\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    148\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel_output\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel_output\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    149\u001b[0m         \u001b[0;31m#self.model_output = self.model.model_output # this allows the TreeEnsemble to translate model outputs types by how it loads the model\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/foo/venv/bar/lib64/python3.9/site-packages/shap/explainers/_tree.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, model, data, data_missing, model_output)\u001b[0m\n\u001b[1;32m    975\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbase_offset\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minit_params\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mparam_idx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    976\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 977\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Model type not yet supported by TreeExplainer: \"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    978\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    979\u001b[0m         \u001b[0;31m# build a dense numpy version of all the tree objects\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: Model type not yet supported by TreeExplainer: <class 'sksurv.ensemble.forest.RandomSurvivalForest'>"
     ]
    }
   ],
   "source": [
    "expl = shap.TreeExplainer(rsf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "argv": [
    "python",
    "-m",
    "ipykernel_launcher",
    "-f",
    "{connection_file}"
   ],
   "display_name": "Python 3 (ipykernel)",
   "env": null,
   "interrupt_mode": "signal",
   "language": "python",
   "metadata": {
    "debugger": true
   },
   "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.8.5"
  },
  "name": "random-survival-forest-shap.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}