{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PMF Example \n",
    "\n",
    "In this notebook, we will present a simple example of the `PMF` method described in [this paper](https://arxiv.org/abs/1910.06724).\n",
    "\n",
    "For a more verbose introduction to `pycox` see [this notebook](https://nbviewer.jupyter.org/github/havakv/pycox/blob/master/examples/01_introduction.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# For preprocessing\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn_pandas import DataFrameMapper \n",
    "\n",
    "import torch # For building the networks \n",
    "import torchtuples as tt # Some useful functions\n",
    "\n",
    "from pycox.datasets import metabric\n",
    "from pycox.models import PMF\n",
    "from pycox.evaluation import EvalSurv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Uncomment to install `sklearn-pandas`\n",
    "# ! pip install sklearn-pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1234)\n",
    "_ = torch.manual_seed(123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset\n",
    "\n",
    "We load the METABRIC data set as a pandas DataFrame and split the data in in train, test and validation.\n",
    "\n",
    "The `duration` column gives the observed times and the `event` column contains indicators of whether the observation is an event (1) or a censored observation (0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = metabric.read_df()\n",
    "df_test = df_train.sample(frac=0.2)\n",
    "df_train = df_train.drop(df_test.index)\n",
    "df_val = df_train.sample(frac=0.2)\n",
    "df_train = df_train.drop(df_val.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>x0</th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>x4</th>\n",
       "      <th>x5</th>\n",
       "      <th>x6</th>\n",
       "      <th>x7</th>\n",
       "      <th>x8</th>\n",
       "      <th>duration</th>\n",
       "      <th>event</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>5.603834</td>\n",
       "      <td>7.811392</td>\n",
       "      <td>10.797988</td>\n",
       "      <td>5.967607</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>56.840000</td>\n",
       "      <td>99.333336</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>5.284882</td>\n",
       "      <td>9.581043</td>\n",
       "      <td>10.204620</td>\n",
       "      <td>5.664970</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>85.940002</td>\n",
       "      <td>95.733330</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>6.654017</td>\n",
       "      <td>5.341846</td>\n",
       "      <td>8.646379</td>\n",
       "      <td>5.655888</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>66.910004</td>\n",
       "      <td>239.300003</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5.456747</td>\n",
       "      <td>5.339741</td>\n",
       "      <td>10.555724</td>\n",
       "      <td>6.008429</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>67.849998</td>\n",
       "      <td>56.933334</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5</td>\n",
       "      <td>5.425826</td>\n",
       "      <td>6.331182</td>\n",
       "      <td>10.455145</td>\n",
       "      <td>5.749053</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>70.519997</td>\n",
       "      <td>123.533333</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x0        x1         x2        x3   x4   x5   x6   x7         x8  \\\n",
       "0  5.603834  7.811392  10.797988  5.967607  1.0  1.0  0.0  1.0  56.840000   \n",
       "1  5.284882  9.581043  10.204620  5.664970  1.0  0.0  0.0  1.0  85.940002   \n",
       "3  6.654017  5.341846   8.646379  5.655888  0.0  0.0  0.0  0.0  66.910004   \n",
       "4  5.456747  5.339741  10.555724  6.008429  1.0  0.0  0.0  1.0  67.849998   \n",
       "5  5.425826  6.331182  10.455145  5.749053  1.0  1.0  0.0  1.0  70.519997   \n",
       "\n",
       "     duration  event  \n",
       "0   99.333336      0  \n",
       "1   95.733330      1  \n",
       "3  239.300003      0  \n",
       "4   56.933334      1  \n",
       "5  123.533333      0  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature transforms\n",
    "\n",
    "The METABRIC dataset has  9 covariates: `x0, ..., x8`.\n",
    "We will standardize the 5 numerical covariates, and leave the binary covariates as is.\n",
    "Note that PyTorch require variables of type `'float32'`.\n",
    "\n",
    "We like using the `sklearn_pandas.DataFrameMapper` to make feature mappers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols_standardize = ['x0', 'x1', 'x2', 'x3', 'x8']\n",
    "cols_leave = ['x4', 'x5', 'x6', 'x7']\n",
    "\n",
    "standardize = [([col], StandardScaler()) for col in cols_standardize]\n",
    "leave = [(col, None) for col in cols_leave]\n",
    "\n",
    "x_mapper = DataFrameMapper(standardize + leave)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = x_mapper.fit_transform(df_train).astype('float32')\n",
    "x_val = x_mapper.transform(df_val).astype('float32')\n",
    "x_test = x_mapper.transform(df_test).astype('float32')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Label transforms\n",
    "\n",
    "The survival methods require individual label transforms, so we have included a proposed `label_transform` for each method.\n",
    "In this case `label_transform` is just a shorthand for the class `pycox.preprocessing.label_transforms.LabTransDiscreteTime`.\n",
    "\n",
    "The `PMF` is a discrete-time method, meaning it requires discretization of the event times to be applied to continuous-time data.\n",
    "We let `num_durations` define the size of this (equidistant) discretization grid, meaning our network will have `num_durations` output nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_durations = 10\n",
    "labtrans = PMF.label_transform(num_durations)\n",
    "get_target = lambda df: (df['duration'].values, df['event'].values)\n",
    "y_train = labtrans.fit_transform(*get_target(df_train))\n",
    "y_val = labtrans.transform(*get_target(df_val))\n",
    "\n",
    "train = (x_train, y_train)\n",
    "val = (x_val, y_val)\n",
    "\n",
    "# We don't need to transform the test labels\n",
    "durations_test, events_test = get_target(df_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pycox.preprocessing.label_transforms.LabTransDiscreteTime"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(labtrans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural net\n",
    "\n",
    "We make a neural net with `torch`.\n",
    "For simple network structures, we can use the `MLPVanilla` provided by `torchtuples`.\n",
    "For building more advanced network architectures, see for example [the tutorials by PyTroch](https://pytorch.org/tutorials/).\n",
    "\n",
    "The following net is an MLP with two hidden layers (with 32 nodes each), ReLU activations, and `out_features` output nodes.\n",
    "We also have batch normalization and dropout between the layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_features = x_train.shape[1]\n",
    "num_nodes = [32, 32]\n",
    "out_features = labtrans.out_features\n",
    "batch_norm = True\n",
    "dropout = 0.1\n",
    "\n",
    "net = tt.practical.MLPVanilla(in_features, num_nodes, out_features, batch_norm, dropout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you instead want to build this network with `torch` you can uncomment the following code.\n",
    "It is essentially equivalent to the `MLPVanilla`, but without the `torch.nn.init.kaiming_normal_` weight initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# net = torch.nn.Sequential(\n",
    "#     torch.nn.Linear(in_features, 32),\n",
    "#     torch.nn.ReLU(),\n",
    "#     torch.nn.BatchNorm1d(32),\n",
    "#     torch.nn.Dropout(0.1),\n",
    "    \n",
    "#     torch.nn.Linear(32, 32),\n",
    "#     torch.nn.ReLU(),\n",
    "#     torch.nn.BatchNorm1d(32),\n",
    "#     torch.nn.Dropout(0.1),\n",
    "    \n",
    "#     torch.nn.Linear(32, out_features)\n",
    "# )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the model\n",
    "\n",
    "To train the model we need to define an optimizer. You can choose any `torch.optim` optimizer, but here we instead use one from `tt.optim` as it has some added functionality.\n",
    "We use the `Adam` optimizer, but instead of choosing a learning rate, we will use the scheme proposed by [Smith 2017](https://arxiv.org/pdf/1506.01186.pdf) to find a suitable learning rate with `model.lr_finder`. See [this post](https://towardsdatascience.com/finding-good-learning-rate-and-the-one-cycle-policy-7159fe1db5d6) for an explanation.\n",
    "\n",
    "We also set `duration_index` which connects the output nodes of the network the the discretization times. This is only useful for prediction and does not affect the training procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PMF(net, tt.optim.Adam, duration_index=labtrans.cuts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3jUZdbw8e9JT0gIkBBaDAm9SoAgVWkWRFSsaxdFkV31cZu7uvvsuqvr6vPa1q6IiLqKq4gdLKsgKr0TmrRAQksgEJKQST3vHzPBgCkzyUwm5Xyuay4yv3pmRE7u8ju3qCrGGGOMOwL8HYAxxpjGw5KGMcYYt1nSMMYY4zZLGsYYY9xmScMYY4zbLGkYY4xxW5C/A/Cm2NhYTUxM9HcYxhjTqKxevfqwqrZ159gmlTQSExNZtWqVv8MwxphGRUT2uHusdU8ZY4xxmyUNY4wxbrOkYYwxxm1NakzDGNM0FRcXk5GRgcPh8HcojVpYWBjx8fEEBwfX+hqWNIwxDV5GRgZRUVEkJiYiIv4Op1FSVY4cOUJGRgZJSUm1vo51TxljGjyHw0FMTIwljDoQEWJiYurcWrOWRiVUlWW7sgkNDiAuKpS2UaGEBgX6OyxjmjVLGHXnje/QkkYlluw8wvUzl5+yrVtcJJ/ePYqwYEsexpjmy7qnKrEyLZsAgZdvHMyjl/fnmiFnsCMzjw0ZOf4OzRjjB8eOHeOFF17w+LyJEydy7Ngxj8+bMmUKc+fO9fi8+mBJoxJr9h6jR7soLujbnmvOSuAPE3oBzmRijGl+qkoapaWl1Z43f/58WrVq5auw/MK6p05TVqas23uUi87seHJbmxYhdIuLZJUlDWP87u+fbGLz/uNevWafji154OK+Ve6/77772LlzJ8nJyQQHBxMZGUmHDh1Yt24dmzdvZvLkyaSnp+NwOLjnnnuYNm0a8FNpo7y8PC688EJGjRrFkiVL6NSpEx999BHh4eE1xvb111/z+9//npKSEoYMGcKLL75IaGgo9913Hx9//DFBQUGcf/75PP7447z33nv8/e9/JzAwkOjoaBYvXuy176icJY3T7Dqcz3FHCQMTTv3tIKVza+ZvPEBZmRIQYANyxjQnjz76KKmpqaxbt45FixZx0UUXkZqaenLq6qxZs2jTpg0FBQUMGTKEK664gpiYmFOusX37dubMmcMrr7zC1Vdfzfvvv88NN9xQ7X0dDgdTpkzh66+/pkePHtx00028+OKL3HTTTXzwwQds3boVETnZBfbggw/yxRdf0KlTp1p1i7nDksZp1uw9CsCg05NGYhveWZnO9sw8eraP8kdoxhiotkVQX84666xTnnV45pln+OCDDwBIT09n+/btP0saSUlJJCcnAzB48GDS0tJqvM+2bdtISkqiR48eANx88808//zz3HXXXYSFhXHbbbdx0UUXMWnSJABGjhzJlClTuPrqq7n88su98VF/xsY0TrN27zFahgXRJTbylO1DElsDNq5hjIEWLVqc/HnRokX897//ZenSpaxfv56BAwdW+ixEaGjoyZ8DAwMpKSmp8T6qWun2oKAgVqxYwRVXXMGHH37IhAkTAHjppZf4xz/+QXp6OsnJyRw5csTTj1YjSxqnWbv3KMkJrX/WBZXQJoK2UaE2rmFMMxQVFUVubm6l+3JycmjdujURERFs3bqVZcuWee2+vXr1Ii0tjR07dgDw5ptvMnr0aPLy8sjJyWHixIn861//Yt26dQDs3LmToUOH8uCDDxIbG0t6errXYinXbLunFmw8QEJMBH07Rp/clldYwo+Hcrmgb/ufHS8iDElszao9R+szTGNMAxATE8PIkSPp168f4eHhtGvX7uS+CRMm8NJLL3HmmWfSs2dPhg0b5rX7hoWF8dprr3HVVVedHAifPn062dnZXHrppTgcDlSVp556CoB7772X7du3o6qMHz+eAQMGeC2WclJV86cxSklJ0ZoWYSotUx76dDOzl6TRLS6SL399zslWxZIdh7lu5nJev/UsRvf4+SJWs77fzYOfbmbp/ePoEF3zrAdjjHds2bKF3r17+zuMJqGy71JEVqtqijvn+7R7SkRmiUimiKRWsf9eEVnneqWKSKmItHHtayUic0Vkq4hsEZHhdY2noKiUX/57NbOXpDEksTU7MvP4dnvWyf3lg+DJ8ZXPq05xjWusSrPWhjGmefL1mMZsYEJVO1X1MVVNVtVk4H7gW1UtHzR4GvhcVXsBA4AtdQnkcF4h17yyjK+2HOKvk/rw9u3DaN8yjJnf7Tp5zNq9x+jatgXREZWXDe7ToSURIYE2rmGM8Yo777yT5OTkU16vvfaav8Oqlk/HNFR1sYgkunn4tcAcABFpCZwDTHFdpwgoqkssf5y7gW0Hj/PSDYNPjllMGZnIowu2snn/cXp3iGJt+jHG94qr8hpBgQEMTGjFSmtpGFPvVLXJFS18/vnn6/V+3hiOaBCzp0QkAmeL5H3Xpi5AFvCaiKwVkZki0qKKc6eJyCoRWZWVlVXZIaRnn+CbbZlMO7vLKYPc1w5JICIkkJnf72LPkRNk5xcxMKF1tbGmdG7D1oPHyXUUe/5BjTG1EhYWxpEjR7zyj15zVb6eRlhYWJ2u01BmT10M/FChayoIGATcrarLReRp4D7gL6efqKozgBngHAiv7OLvrNyLAL84K+GU7dERwVydcgZvLd9Dj3bOB/ZOfxL8dEMS21Cmzq6scyoZLDfGeF98fDwZGRlU9YuhcU/5yn0V7crK8+gaDSVpXIOra8olA8hQ1fL65HNxJg2PFZeW8Z+VGYzrFUenVj+f8XTryCTeWJrGv/77Iy1CAk8mj6okJ7QiMEBYlZZtScOYehIcHFyn1eZM1RakHvToeL93T4lINDAa+Kh8m6oeBNJFpKdr03hgc22u/9XmQxzOK+S6oQmV7k+IieCCvu1xFJcx4AxnQqhOZGgQvTtEsWyXDYYbYxo/R3H1lXpP5+spt3OApUBPEckQkakiMl1Eplc47DLgS1XNP+30u4G3RGQDkAz8szYxvLV8D51ahTO6R9UD3Led7fwNpqauqXIT+rZnRVo28zceqE1IxhjTYBQUeZY0fD176lo3jpmNc2ru6dvXAW49bFKVXVl5/LDjCL8/v0e1LYhBCa154fpBDE1q49Z17xjdla+2ZHLf+xsYcEarSru9jDGmMXCUNKCWhr/NWbGXoADh6pQzqj1ORJjYvwMxkaHVHlcuODCAZ65JprRM+c076ygtc39Gx/JdR5j49HfcPGsFf/5gIy99u5N16b4pYWyMMTUpKCrz6PgmmzQcxaXMXZ3BeX3aEdeyblPMKtM5pgUPTe7HirRsnl+4w+3zXvshjb3Zzum9n208wKMLtnLli0v4brvNCjHG1L8GNabhT/M3HuDoiWKuH9rZZ/e4fFA8k5M78vTX21m9p+aB8VxHMd9sy+TKwfF8cvco1v31fFb/77l0i4vkl/9ew6b9tga5MaZ+WdLA+RDLjMW76BYXyYiuMTWfUAcPTe5Hh+gw7p27gaKS6pt5X246RFFJGRcP+Gkp2ZjIUF67ZQhRYUHc8tpK9h0r8Gm8xhhTUYElDfhmayZbD+byqzFdfb40a1RYMA9d2o9dWfm8+v3uao/9ZMN+OrUK/9mqgB2iw5l9y1kUFJcyZdYKck7Y0+bGmPrR7FsaqspzC3fQqVX4Kb/R+9LYXnGc16cdz36znf1VtBSy84v4fvthLh7QsdL6OT3bR/HyjYNJO5LP3e+stXIJxph6UVDczAfCl+3KZu3eY0wf3YXgwPr7eH+d1IfSMuXhzyovxrsg9QAlZcrFAzpUeY0RXWP566Q+LP4xizeX7fE4hq+3HOKR+Vs8ms1ljGnemn1L44VFO4iNDOWqGqbZetsZbSK4c2w3Ptt4gO+3H/7Z/k/W76dr2xb06dCy2uvcMKwzY3q25Z/zt7Aj0/2aMMdOFPG799bz8uJdPPHlNo/jN8Y0T806aRQUlfLd9sPcdnYSYcGB9X7/aed0oXNMBH/9OPWUQfGDOQ6W786usmuqIhHh/11xJuHBgfz23XUUl7rXdHzyqx85XlDM+F5xvLBoJ5+s31+nz2KMaR6a9UB4Zm4hLcOCuL6KOlO+FhYcyAMX92FXVj7T3lzF7sPOyiifbTyAKm6PscS1DOORy/uzISOHZ7/eXuPxWw8e59/L9nDDsM68eMNgUjq35t6560ndZ1N4jTHVa9YtjeOOYqaMSCQqrPKV9+rDuF7t+MukPqzcnc35T33LPz7dzLw1GfTt2JKubSPdvs6Efh24YlA8zy3cwbJdR6o8TlX5+8ebaRkezG/P60FIUAAv3jCY1hEh3PHmag7nFXrjY/nN+vRjXP3yUp5fuKPKSQb1qbRMycotZMuB46zZe9TtlqAxDZGq4vBwIFya0iydmMTeun3Tetq0CPF3KGTmOnj8i228tzoDVbj/wl7cMbqrR9c47ihm8nM/sD+ngGevHcR5fdr97Jj5Gw/wq7fW8NClfblxeOLJ7RszcrjypSV0i4tk1pQhtPPBU/G+dtxRzMSnvyM7v4gTRaWIwIiuMVw/tDMT+1c9ocAXVu85yr3vrWf3kXwq/i8TGxnK5YM6cdXgeLq3i6K4tIx9RwtIO5JPfOsIusW5/4uCMfXNUVxKr798zp7/m7RaVd2q9dekkkZKSoquWrXK32GcInVfDvPW7ON/xnejVYTnyexwXiFTZ69k474cHr6sP9dWWEgqr7CEC55aTFRYEJ/ePYqg02aLLdyWyV1vrSEqLJiZN6fQr1N0nT9PfVFV7p6zlgWpB3n3juHERYXy/poM5q3Zx97sE1w/NIEHLu5LSJDvG8vvrkrnfz9IpX10GJOTOxIbFUpsZChlqny8bj/fbM2kpExp3zKMrLzCk7PXQoICmHlTiq27Yhqso/lFDHzoK0saTc2JohLufGsNC7dlccc5XWgZHsz32w+zes9RikrLePv2oYzoGlvpuVsOHGfq7JUcKyjm6WsGVtpaaYjeXZXOH+Zu4N4LenLn2G4nt5eWKY99sY2Xvt3J0KQ2vHD9ILcLTVaUcfQED326mRW7sxnZLZYJ/dozpmcckaE/FX4uKS3jn/O3MuuH3YzsFsPz1w2qNPEfzivkw7X7SN2XQ3zrCDrHRNCpVTgPfbaFXVl5vHrzEEZ1r/y/jzH+tP9YASMe/caSRlNUXFrGn+Zt5L3VGQD06dCSUd1jOa9PO4YkVl/SPfO4g9vfWMWGfTk8dXUykwd2qo+QK1VapizffYS1e4+xKyufnVl5pB3JJym2BZcP7MSkMzuSfaKISc98z8CEVrw5dWilZe0/XLuPP7y/gbioUO4Y3ZWMoyfYmZnH7sP5dI+L4vphCYzsGvuzigBFJWW8+v1unnFNMBjXK45lu45wJL+IkKAAuraNpPyM3MJi0rMLuGVkIn+e2PtnLbmaZOcXcd0ry9h9OJ9ZU4YwspslDtOw7MrKY9wT31rSaKpUldR9x+nYKszj364Likq58dXlbM/M4+vfjSa2ivN3ZeXx8fr9fLxuP4fzCrlsYCeuHZpAr/aVP19SUlpGxtECsk8UcWan6Er/YS0rU5bvzuazjfv5PPUgh/OKAGjXMpSubSPpHBPBmj3H2HYol+BAITo8mNIy5fNfn1PtWMy69GNMe2MVmbmFhAQGkBgbQUKbFqzZe5Ts/CKSYltwdcoZRIYGkplbSFZuISvSstmVlc8Ffdvx14v70qlVOKVlyuo9R/li00H2HMkHfko0F53ZnssGxlcZQ02O5BVy/czlpB3J59WbLXGYhmXT/hwueub7hpM0RGQWMAnIVNV+ley/F7je9TYI6A20VdVs1/5AYBWwT1Un1XS/pp406mpHZi4XPv0dF5/ZkSd/kXzKvtR9Odw/byMb9+UgAsOSYoiJDHEWWSwtY1BCKwYmtCbPUUJuYTG5jhIyjhaQnn2CElcf/vl92vHsdQMJDfrpGZmColJ+9dZqFm7LIiw4gPG92nHRmR0Y1T2WlhVmuakqmw8c54M1+1j0Yxb/e1FvxvSserXFcvmFJWTlFhLfOvxkwiosKWXBxoO8uWwPq/ccBSBAnMUh41uHc9fYbozvXX/ddEfyCrnuleXsPpzPU79I5qIzTx3En7/xAC9/u5MpIxOZnNypxmd5jPGW1XuOcsWLSxpU0jgHyAPeqCxpnHbsxcBvVHVchW2/xbl6X0tLGt7x+BfbeG7hDt6+bSgjXL/1bjuYyy9mLCU8OJCpo5KYdGZH2kc7f8PPzi9i3poM/rMynX3HCogKCyIqLJiosCA6RIeRGNOCxNgWHMpx8MRXPzK2Z1tevGEwYcGB5Jwo5tbXV7J271H+NLE31w1NICLEp4tF/sy+YwUEBwoxLUJrXP/dl3JOFHPbGytZtecof7u4LzePSKSgqJQHP93MnBV7iQoNIrewhHG94nj4sn50iLbVII3v/bDjMNfPXN5wkgaAiCQCn7qRNN4GFqrqK6738cDrwMPAby1peIejuJTzn1pMUICw4Ndns/+Yg6tfXooA700fTueYFrW+9tvL9/LnDzcysmssD1/WjzveXM2urHyeviaZC+t5imxD5Cgu5e45a/lq8yFuGt6ZpTuPsD0zj+mju/Lrc7vz1vK9PPbFVoIDArhrXDd6tIsiJjKENi1CiIsKq5eZYqZ5+e/mQ9z2xiqPkkb9/tpXBRGJACYAd1XY/C/gD0BUDedOA6YBJCT450nwxiQsOJCHJvfj5lkr+MenW/hmayalZcp/pg2rU8IAuG5oAiFBAfxh7nrGPfEtoUEBzJpiM4fKhQUH8uL1g/jLR5t4Y+keYiNDeXPqWZzd3Tkld+qoJM7tHcd972/kkQVbTzm3RUgg5/Roy7m92zGuVxytG8CzSKbx83R9cGggSQO4GPihwlhG+TjIahEZU92JqjoDmAHOloavA20KRvdoy6QzO/Dmsj1EhQYxZ9owurerNje77crB8YQGBfDy4p38Y3J/ks9oVfNJzUhQYAD/vKwf5/WJ48z4Vj+bkNA5pgVv3z6U9OwCsvIKyc4v4kheIRv35fDfLYdYkHqQAIHrh3bmb5f09WuXm2n8Cooab9K4BphT4f1I4BIRmQiEAS1F5N+qeoNfomuC/jqpD8WlZUw7p4vXH/q7eEDHelvLpDESEcb1qnogXkRIiIkgISbi5LZrgH9M7kfqvuO8s3Ivby7bQ15hCY9deabHU4GNKeeoYbXRyvg9aYhINDAaOJkQVPV+4H7X/jHA7y1heFdcyzBevtGtLkzTQIgI/eOj6R/fnw7RYTz+5Y+UlClPXj2gXteOMU2Ho6G1NERkDjAGiBWRDOABIBhAVV9yHXYZ8KWq5vsyFmOakrvGdSc4MIBHFmyluKSM64YmcOi4g8zcQopKypgyItHGPUyNPK1wCz5OGqp6rRvHzAZmV7N/EbDIWzEZ01TcMborwYEBPPjpZj7fdPCUfd9szeSt24ee8iyMMacrKC4lyMNxMb93Txljau/WUUmcldSGguJS2kWFEdcylKU7jzDtzVVMnb2S128965RnY3IKijlRVGLPgRjAmTTCPVywzpKGMY3c6RMZxvaK4+lrBnLX22u4483VvHJTCunZJ3htSRrz1mTgKC5jWJc2XDMkgQn92vtllUvTMDiKywj1RdIQkQBgANARKAA2qeohjyM0xtSLif078P+uHMDv31vP+Ce+Zd+xAkKCApic3JH41hHMXZ3Br/+zjpYfBXH3uO7cdnaSlS9phhzFpYSHeDaJotqkISJdgT8C5wLbgSycU2B7iMgJ4GXgdVW15cuMaWCuHBxPYUkpr36/m9+f34Nrz0o4WejyrrHdWLb7CDO/283D87eQlVfI/Rf2ssTRzDiKSwkL8m5L4x/Ai8Adelq9ERGJA64DbsRZ7sMY08BcP7Qz1w/t/LPtAQHCiK6xDEuK4W+fbGLG4l3kOkr4x+R+9sBgM1JQXEp4iBeTRnWzn1Q1E2epD2NMIxUQIPz9kr5EhQXx/MKd5BWW2HMfzUhBUanHY1o1dU9dXt1+VZ3n0d2MMQ2OiHDvBb2ICgvm0QVbSc8+wZNXD6BLW1vfvKlzlJQRHe7ZtOyauqcudv0ZB4wAvnG9H4vz2QlLGsY0EdNHdyW+dTh//iCVi575nj9N7MUNwzrbOEcT5igqpX1LzxZ0q6l76hYAEfkU6KOqB1zvOwDP1zJOY0wDNenMjqR0bsMf3t/AXz7axOebDpLSuQ3hIYGEBwfSsVU45/aOs0TSRDhKvNw9VUFiecJwOQT08OhOxphGoX10GK/fMoR/L9/LE19u44cdR07ZP+nMDjx25QCPB1BNw1NQ5LuH+xaJyBc4K9EqzqKbCz0LzxjTWIgINw7rzI3DOlNWpjhKSikoKuU/q9J57IttpB3JZ8aNKXRsZU+WN2aOYs9bGm5NkVDVu4CXcD7glwzMUNW7PY7QGNPoBAQIESFBxESG8qsx3Zh5Uwpph09wyXM/nFyD3TROjuIy3yQNlzXAZ6r6G+ALEfHOqj3GmEZlfO92fPCrEbQIDWTq6yvJLyzxd0imFkrLlKLSMo+7p9xKGiJyOzAX5xPgAJ2ADz26kzGmyejeLoonrx7AsRPFzF2d4e9wTC2Ul0UPC/bsmRx3j74T52p6xwFUdTvOabjGmGZqcOc2DExoxavf76a0zFZabmwKXEnD0wkN7iaNQlUtKn8jIkE4B8SNMc3Y7Wd3YW/2Cb7afLDmg02DcrKl4WHtKXeTxrci8icgXETOA94DPqnpJBGZJSKZIpJaxf57RWSd65UqIqUi0kZEzhCRhSKyRUQ2icg97n8kY0x9uaBve85oE84r3+32dyjGQyeTho9aGvfhrHC7EbgDmA/8rxvnzQYmVLVTVR9T1WRVTca5Jvi3qpoNlAC/U9XewDDgThHp42asxph6Ehgg3DoyidV7jrJmr82kakwKipzFyX0yEK6qZar6iqpepapXun6usXtKVRcD2W7Gci3O50BQ1QOqusb1cy6wBefguzGmgbk65QxahgUx87td/g7FeMBR4sOBcBEZKSJficiPIrJLRHaLiNf+hohIBM4WyfuV7EsEBgLLvXU/Y4z3tAgN4rqhnfk89SDp2Sf8HY5xU0GRayDcR89pvAo8CYwChgAprj+95WLgB1fX1EkiEokzkfxaVY9XdqKITBORVSKyKisry4shGWPcNWVEIgEivPjtTn+HYtz005Rb3ySNHFVdoKqZqnqk/OVhjNW5BlfXVDkRCcaZMN6qrgS7qs5Q1RRVTWnbtq0XQzLGuKt9dBg3Du/M28v38tVmWwm6MSjwRdIQkUEiMghYKCKPicjw8m2u7XUmItHAaOCjCtsEZ+tmi6o+6Y37GGN8674Le9G/UzS/e3eddVM1AoXFzoFwT8c0aipY+MRp71Mq/KzAuOpOFpE5wBggVkQygAeAYABVfcl12GXAl6qaX+HUkTiXkd0oIutc2/6kqvNriNcY4yehQYE8f90gLnr2O+56ew3vTh9OqIfPAJj6c/LhPm9WuVXVsQAi0kVVTxn4FpEuNV28uuViKxwzG+fU3IrbvgesYL8xjUxCTASPXTmA6f9ezSPzt/K3S/r6OyRTBV8/ET63km3veXQnY0yzMKFfe6aOSmL2kjQ+XLvP3+GYKtT2ifCa1gjvBfQFok9bL7wlEOZZiMaY5uKPE3qRui+H3723nrDgACb06+DvkMxpCopLCQkKICDAs06dmloaPYFJQCuc02LLX4OA22sRpzGmGQgJCuDVKUNIPqMVd729lv/ajKoGp7C4jLAgzwbBoeYxjY+Aj0RkuKourW1wxpjmJzI0iNduGcKNr67gV2+t4eWbBjO2pxXHbigKikprtWSvu2kmXUQ+cBUfPCQi74tIvMd3M8Y0Ky3DgnnjlrPo0T6SO95cbSv9NSAFxZ6vDw7uJ43XgI+BjjhrQH3i2maMMdWKjgjmzVuH0q5lKPe8s5bjjmJ/h2So3frg4H7SiFPV11S1xPWaDdjj18YYt7RuEcLT1wzkQI6DP3+Qihv1To2PFfg4aWSJyA0iEuh63QB4s4yIMaaJG5TQmt+c251P1u/n/TU2FdffCovLPH4aHNxPGrcCVwMHXa8rXduMMcZtvxzTjaFJbfjrR6nsPpxf8wnGZ3w6pqGqe1X1ElVt63pNVtU9Ht/NGNOsBQYIT/0imeDAAO55Zy3FpWX+DqnZ8umYhojE2+wpY4w3dGwVzv9d0Z8NGTm8tMhKqfuLzZ4yxjQaE/p14OIBHXnmm+1sO5jr73CaJUdxqcfrg4P7SaOtzZ4yxnjT3y7uQ8uwYP4wdz0l1k1V7xzFZR7XnQL3k8Zhmz1ljPGmmMhQ/n5pX9Zn5PDq97v9HU6zU1BcSnhI/cyeOoDNnjLGeMFF/TtwQd92PPHVj+zMyvN3OM1GcWkZpWXqu5bGabOn4mz2lDHGG0SEhyb3Izw4kN++u578whJ/h9Qs1HYtDXB/9lSSiDwpIvNE5OPylxvnzXLNuEqtYv+9IrLO9UoVkVIRaePaN0FEtonIDhG5z7OPZYxpLOKiwvi/K84kdV8ON766nJwCKzPia46i2q0PDu53T30IpAHP4lwCtvxVk9nAhKp2qupjqpqsqsnA/cC3qpotIoHA88CFQB/gWhHp42asxphGZkK/9jx/3SA27svh+pnLyM4v8ndITZrj5PrgvksaDlV9RlUXquq35a+aTlLVxUC2m/e4Fpjj+vksYIeq7lLVIuAd4FI3r2OMaYQm9GvPjJtS2H4oj2tmLCUz1+HvkJqs2q4PDu4njadF5AERGS4ig8pfHt+tCiISgbNF8r5rUycgvcIhGa5txpgmbGzPOF6bMoSMowX8+p11/g6nyTq51Gstak9VuwhTBf2BG4FxQPmEanW994aLgR9UtbxVUtn6g5WWxRSRacA0gISEBC+FY4zxlxHdYrlnfHceWbCVLQeO07tDS3+H1OTUR0vjMqCLqo5W1bGul7cSBsA1/NQ1Bc6WxRkV3scD+ys7UVVnqGqKqqa0bWvPGxrTFFwzJIHw4EBm/5Dm71CapPKWRqgPk8Z6nOuEe52IRAOjgY8qbF4JdHfN2grBmVRqnK1ljGkaoiOCuWxQJ7oKTCgAABnuSURBVD5ct88GxX3AUQ8tjXbAVhH5wsMpt3OApUBPEckQkakiMl1Eplc47DLgS1U9WSdZVUuAu4AvgC3Au6q6yd0PZYxp/KaMSKSwpIw5K/b6O5Qmpy7Pabg7pvGAx1cGVPVaN46ZjXNq7unb5wPza3NfY0zj16NdFKO6xfLvZXuYdk4XggM9H7Q1lftpyq3vyoisAr5zTbM9AEQDSzy+mzHGeGDKiEQO5Dj4ctMhf4fSpBQU+b57ajEQJiKdgK+BW6ikdWCMMd40rlccnWMimL3EChp6k6PE90+Ei6qeAC4HnlXVy4C+Ht/NGGM8EBAg3DQ8kZVpR9mQcczf4TQZ5WVEQoN81z0lIjIcuB74zLXN8xRljDEeuiolnpZhQVw7YxnPfbP9ZNeKqb2C4lLCggMQqeyRuOq5mzTuwVkb6gNV3SQiXYCFHt/NGGM81DIsmI/uGsWo7rE8/uWPjH18Ee+uSke10ud9jRscxWW1Gs8AN2dPuWpILa7wfhfwP7W6ozHGeCgptgUv35jCit3ZPDx/C3+Yu4Gs3ELuHNvN36E1SrVdHxxqaGmIyAwR6V/FvhYicquIXF+rOxtjjIfOSmrDh78awUVnduCpr35kY0aOv0NqlBzFpbUaBIeau6deAP4iIltE5D0RecG1RsZ3OKfcRgFza3VnY4ypBRHh4cn9iI0M5df/WWtjHLXgs6ShqutU9WpgCM71Lb7DWc7jNlUdoKpPq2phre5sjDG11CoihMevGsDOrHweXbDF3+E0Oo7islo92Afuj2nkAYtqdQdjjPGBUd1juWVkIq/9kMa43u0Y3cMKlrqroLi0ViVEwP3ZU8YY0+D8cUIvusdFcu9768mz9cXdVlDko4FwY4xpyMKCA3n4sv5k5hYyf8MBf4fTaDhKSmtVFh1qkTREJEBEbFUUY0yDMCSxNV1iWzB3dYa/Q2k0HL5uaYjI2yLSUkRaAJuBbSJyb63uaIwxXiQiXDE4nhVp2ew9csLf4TQKjpLaD4S7e1YfVT0OTMZZrjwB5/Kvxhjjd5cN7IQIvL/GWhvuqI8xjWARCcaZND5S1WKqWLPbGGPqW8dW4YzsGsv7azIoK7N/mqqjqq7aU75NGi8DaUALYLGIdAaO1+qOxhjjA1cM7kTG0QJWpGX7O5QGrbCkfAEmHyYNVX1GVTup6kR12gOMrek819PjmSKSWs0xY0RknYhsEpFvK2z/jWtbqojMEZEwtz6RMaZZuqBveyJDg3jfBsSrVZf1wcH9gfB7XAPhIiKvisgaYJwbp84GJlRz3VY4S5Vcoqp9gatc2zvhLIiYoqr9cJZhv8adWI0xzVNESBAT+7dn/sYDnCiyZzaq8tNSr77tnrrVNRB+PtAW58p9j9Z0kqs6bnVtxeuAeaq613V8ZoV9QUC4iAQBEcB+N2M1xjRTVwyKJ7+olM9TD/o7lAaroLylEeLb2VPlK3VMBF5T1fUVttVFD6C1iCwSkdUichOAqu4DHgf24lyTPEdVv6w0MJFpIrJKRFZlZWV5ISRjTGM1JLENCW0ieG+VdVFVpbx7KizIty2N1SLyJc6k8YWIRAFltbrjqYKAwcBFwAU4K+r2EJHWwKVAEtARaCEiN1R2AVWdoaopqprStq3VnjGmOQsIEK4fmsDSXUf4cO0+f4fTIJW3NMJ8XHtqKnAfMMS1VngIzi6qusoAPlfVfFU9jHOhpwHAucBuVc1yTe+dB4zwwv2MMU3c1FFJDElszZ8/2MiurDx/h9PglK8P7tOBcFUtA+KB/xWRx4ERqrqhVnc81UfA2SISJCIRwFBgC85uqWEiEiHORWzHu7YbY0y1ggIDeObagQQHBXDX22tPdscYp3xX0ojwZUtDRB7FuU74Ztfrf0TkETfOmwMsBXqKSIaITBWR6SIyHUBVtwCfAxuAFcBMVU1V1eU4F3daA2x0xTnD409njGmWOkSH8/iVA9h84DiPzLffNyvKdRQDEBUWXKvz3VpPA+dYRrKrxYGIvA6sBe6v7iRVvbamC6vqY8BjlWx/AHjAzfiMMeYU5/Zpx9RRSbz6/W6Gd41hQr8O/g6pQch1OKcjR4W5+8//qTyZc9Wqws/RtbqbMcbUoz9O6EWfDi15ZMFWVK28CHBy3RFfJ41HgLUiMtvVylgN/LNWdzTGmHoSEhTA1FFJ7DlyglV7jvo7nAbhuKOYkMAAQn055VZV5wDDcM5imgcMV9V3anVHY4ypRxf2b0+LkEDm2rMbgLN7qratDKghaYjIoPIX0AHnFNl0oKNrmzHGNGjO8iId+MzKiwB1Txo1nflENfsU9+pPGWOMX105OJ73VmfwxaaDXDYw3t/h+FWuo7jWM6eghqShqjVWsgUQkfNU9ataR2GMMT5UXl5k7uoMSxq+7J7ywP956TrGGON1AQHCFYPiWbLzCBlHm/eSsHkNJGl4o3ihMcb4zOWDOqEKH6xp3jWpch3FRIbWvnvKW0nDJkAbYxq0M9pEMLxLDHPXZDTrZzYaSveUMcY0eFcOjm/Wz2yUlSl5RSW0bABJI81L1zHGGJ+5sL9zSdiXv93p71D8Iq+oBNXa150C92tPISIjgMSK56jqG64/L691BMYYU08iQoK4e1w3Hlmwlf9uPsS5fdr5O6R6Vde6U+B+lds3ca6kNwoY4nql1PquxhjjJ7eOSqJ7XCR/+2QTBUXNq2x6XSvcgvstjRSgjzbn0SNjTJMQHBjAg5f249pXlvHCoh387vye/g6p3uS5WhqR9TCmkQq0r/VdjDGmARneNYbJyR15+dtdzWp1P593T4nIJyLyMRALbBaRL0Tk4/JXre9qjDF+9qeLehMaFMADH29qNlNwj7u6p+oye6qmMx+v9ZUBEZkFTAIyVbVfFceMAf4FBAOHVXW0a3srYCbQD+dzILeq6tK6xGOMMeXiosL43fk9+Nsnm5m3Zh9XDG765UV+amn4rvbUtwAikgQcUFWH63044M60g9nAc8Able10JYYXgAmquldE4irsfhr4XFWvFJEQIMKN+xljjNtuHJ7I/I0H+ctHqQw4oxXd4iL9HZJP1dvsKeA9oKzC+1LXtmqp6mIgu5pDrgPmqepe1/GZACLSEjgHeNW1vUhVj7kZqzHGuCUwQHjm2oGEBQdy51trmvxsqlxHMYEBQnhw7RZgAveTRpCqFpW/cf0cUuu7/qQH0FpEFonIahG5ybW9C5AFvCYia0Vkpoi0qOwCIjJNRFaJyKqsrCwvhGSMaU7aR4fx5NUD2HYol79/ssnf4fhUXqGzhIhI7csFups0skTkkvI3InIpcLjWd/1JEDAYuAi4APiLiPRwbR8EvKiqA4F84L7KLqCqM1Q1RVVT2rZt64WQjDHNzZiecfxqTFfeWZnOh2ubbkHDXEcJkaG175oC95PGdOBPIrJXRNKBPwJ31OnOThk4xy3yVfUwsBgY4NqeoarLXcfNxZlEjDHGJ357Xg+GJLbmTx9sZN+xAn+H4xN1XYAJ3F8jfKeqDgP64HzIb4Sq7qjTnZ0+As4WkSARiQCGAltU9SCQLiLlT92MBzZ74X7GGFOpoMAAnvpFMiWlynPfbPd3OD5xvI4VbsGz2lMXAX2BsPL+MFV9sIZz5gBjgFgRyQAewDm1FlV9SVW3iMjnwAacA+0zVTXVdfrdwFuumVO7gFs8+FzGGOOx+NYRXDc0gTeX7eGOc7qSGFvpUGqjlesooVOrsDpdw62kISIv4ZzyOhbnsxNXAitqOk9Vr3XjmMeAxyrZvg6rb2WMqWe/GtuVd1bu5emvt/PUL5L9HY5XObunoup0DXfHNEao6k3AUVX9OzAcOKNOdzbGmAYoLiqMm0ck8uG6fWw/lOvvcLyqrgswgftJo3xU6ISIdASKgaQ63dkYYxqo6ed0pUVIEE9+9aO/Q/EaVSWvsP5mT33qenr7MWANzkWX5tTpzsYY00C1bhHC1FFJLEg9SOq+HH+H4xUFxaWUlmm9zZ56SFWPqer7QGegl6r+tU53NsaYBmzq2UlEhwfzxJfb/B2KV3ijhAi4vwhTmIj8VkTmAW8Dt4pI3YbgjTGmAWsZFswdo7uwcFsW69IbfxWjnxZgqp/uqTdwTrd9FmcBwt7Am3W6szHGNHA3DU+kVURwk3hu47irpdGyjt1T7qacnqo6oML7hSKyvk53NsaYBi4yNIhbRybx5Fc/sml/Dn07Rvs7pFqr1+4pYK2IDCt/IyJDgR/qdGdjjGkEbh6RSFRoEM99440iGP7jjaVeoeaV+zaKyAac5T2WiEiaiOwGluIsXW6MMU1adHgwU0YmsiD1INsONt7nNn4a0/Dt7KlJwMXABJzPZYzGWRYkCWdlWmOMafJuHZlEi5BAnlvYeFsb9dI9pap7qnvV6c7GGNNItG4Rwg3DO/Pphv3szMrzdzi1kusoRgQiQ+pnTMMYY5q128/uQmhQAM830tbGcUcJkSFBBATUfgEmsKRhjDFuiY0M5dqzEvh43X4yjp7wdzge80bdKbCkYYwxbrvt7C4AzPxut58j8Zw3FmACSxrGGOO2Tq3CuSS5I/9Zmc7R/CJ/h+ORvMKSOk+3BUsaxhjjkemju1JQXMrrS9P8HYpHGkX3lIjMEpFMEUmt5pgxIrJORDaJyLen7QsUkbUi8qkv4zTGGHf1aBfF+F5xvL4kjRNFJf4Ox22NpXtqNs5nPCrlKrf+AnCJqvYFrjrtkHuALT6LzhhjamH6mK4cPVHMf1am+zsUtzWKloaqLgayqznkOmCequ51HZ9ZvkNE4nE+QDjTlzEaY4ynhiS2IaVza2Z+t5vi0jJ/h+OWRpE03NADaC0ii0RktYjcVGHfv4A/AI3jv4gxplmZPror+44V8Mn6/f4OpUaO4lKKSsvqXOEW/J80goDBOFsUFwB/EZEeIjIJyFTV1TVdQESmicgqEVmVlZXl43CNMcZpXK84eraL4oVFOyktU3+HU628Qlexwjou9Qr+TxoZwOeqmq+qh4HFwABgJHCJiKQB7wDjROTflV1AVWeoaoqqprRt27a+4jbGNHMBAcJd47qxIzOPBakH/B1OtbxVdwr8nzQ+As4WkSARicBZTXeLqt6vqvGqmghcA3yjqjf4M1BjjDndxP4d6Nq2Bc9+vYOyBtza8FaFW/D9lNs5OMuo9xSRDBGZKiLTRWQ6gKpuAT4HNgArgJmqWuX0XGOMaUgCA4S7x3Vn26Fcvtx80N/hVMmbLY26X6EaqnqtG8c8BjxWzf5FwCLvRWWMMd5z8YCOPPP1dp7+egcX9G2PSN0KAvqCt9YHB/93TxljTKMWGCDcObYbWw4c56vNh/wdTqW8tT44WNIwxpg6uzS5IwltInjmm+2oNryxjaY0EG6MMY1eUGAAd43tRuq+43y2seHNpCpfH7xFE5hya4wxTcJlgzrRv1M098/bSHp2w1pvI9dRTHhwIMGBdf8n35KGMcZ4QXBgAM9fNwiAO99eQ2FJqZ8j+om3SoiAJQ1jjPGahJgIHrtyABsycnhk/lZ/h3NSbmGxJQ1jjGmIJvRrz60jk5i9JI35DWR8w9nSqPvMKbCkYYwxXnffhb0YcEYr/jh3AwdyCvwdDsete8oYYxqukKAAnr1mIIUlZTz11Y/+Doc8h3VPGWNMg5YQE8FNwzszd3UGPx7K9WssuY4SokKte8oYYxq0O8d2o0VoEP/vc/8OitvsKWOMaQRatwjhl2O68t8tmazYXd0ipr5TXFpGQXGpDYQbY0xjcOvIJNq3DOORBVv8UmLk0HEHAC3DraVhjDENXlhwIL85rztr9x7ji031Xz793ZXpiMDYnnFeuZ5PS6MbY4yBKwbFM/O73dw3byNvr0inbWQobaNCuaBvOwYmtPbZfR3Fpfx7+V7G92pHYmwLr1zTWhrGGONjQYEBPPWLZFI6tyHnRBFLdx7m1e93cfOsFRx3rXXhCx+v2092fhG3jkz02jWtpWGMMfWgX6doZt6ccvJ96r4cJj37Pa//kMbd47t7/X6qyqwfdtOrfRTDu8Z47bq+Xu51lohkikiVS7iKyBgRWScim0TkW9e2M0RkoYhscW2/x5dxGmNMfevXKZpze8cx8/vdJ1fW86alu46w9WAut45M8upqgr7unpoNTKhqp4i0Al4ALlHVvsBVrl0lwO9UtTcwDLhTRPr4OFZjjKlX94zvQU5BMW8s3eP1a8/6Po02LUK4JLmjV6/r06ShqouB6iYnXwfMU9W9ruMzXX8eUNU1rp9zgS1AJ1/Gaowx9a1/fDTjesXxyne7yCss8dp10w7n8/XWQ9wwNIGw4ECvXRf8P6bRAwgWkUVAFPC0qr5R8QARSQQGAssru4CITAOmASQkJPgwVGOM8b57xnfn0ud/4PUladw5tlutrrEu/RgLt2YSHhJIeHAg3+84TFCAcMOwzl6O1v9JIwgYDIwHwoGlIrJMVX8EEJFI4H3g16p6vLILqOoMYAZASkpKw1uc1xhjqjHgjFaM6dmWmd/tYsqIRI+XZE3dl8O1M5ZRUHzqok9Xp8QT1zLMm6EC/k8aGcBhVc0H8kVkMTAA+FFEgnEmjLdUdZ4/gzTGGF+6Z3x3LnthCXfPWcuNwzszqlusW0uzZh53cPsbq2gVEczCX40hOjyYguJSCopLaRcV6pNY/Z00PgKeE5EgIAQYCjwlzqH+V4EtqvqkPwM0xhhfG5jQmv8Z353XftjNN1szaR0RzMT+HbhnfPcqWwsFRaXc/sYqcgqKeW/6cNpHO48LD/HuGMbpfJo0RGQOMAaIFZEM4AEgGEBVX1LVLSLyObABKANmqmqqiIwCbgQ2isg61+X+pKrzfRmvMcb4y2/P68GdY7uy+MfDfLJ+P3NXZ7B05xHemTbsZ4mjrEz5/Xvr2bAvhxk3ptC3Y3S9xSn+KKDlKykpKbpq1Sp/h2GMMXW2Ki2bm2atoEN0GO9MG05bV3dT5nEHf/tkE/M3HuRPE3sx7Zyudb6XiKxW1ZSaj7QyIsYY0yClJLbhtSlD2H/MwXWvLCMz18EbS9MY/8S3/HdLJn+Y0JPbz+5S73FZS8MYYxqwpTuPcMvsFZQpFJWUMapbLA9N7keSlwoQgmctDX8PhBtjjKnG8K4xvHrzEB77Yhu3jEzkkgEdvVoWxFOWNIwxpoEb2S2Wkd1i/R0GYGMaxhhjPGBJwxhjjNssaRhjjHGbJQ1jjDFus6RhjDHGbZY0jDHGuM2ShjHGGLdZ0jDGGOO2JlVGRERygW21PD0ayKnF/tO3V/e+sp8rbosFDnsUdfWx1bS/KcXuTrxV/dyQY6/43mKvPraa9tcU++nvm1PsnVW1bU3BA6CqTeYFrKrDuTNqs//07dW9r+zn07bVKn6L3b14q/kcDTb2ar5vi93LsVcTp8Ve4WXdUz/5pJb7T99e3fvKfq7pvu6w2H++zdOfa6M+Yq/43mJ37/zaxn76e4u9Ek2te2qVulmpsSFqzPFb7P5hsftHc469qbU0Zvg7gDpqzPFb7P5hsftHs429SbU0jDHG+FZTa2kYY4zxIUsaxhhj3GZJwxhjjNuaTdIQkbNF5CURmSkiS/wdjydEJEBEHhaRZ0XkZn/H4wkRGSMi37m++zH+jsdTItJCRFaLyCR/x+IJEent+s7nisgv/R2Pp0Rksoi8IiIficj5/o7HEyLSRUReFZG5/o7FHa6/46+7vu/razq+USQNEZklIpkiknra9gkisk1EdojIfdVdQ1W/U9XpwKfA676MtyJvxA5cCnQCioEMX8V6Oi/FrkAeEEbjix3gj8C7vomycl76+77F9ff9aqBep4Z6Kf4PVfV2YArwCx+Gewovxb5LVaf6NtLqefg5Lgfmur7vS2q8eF2eDKyvF3AOMAhIrbAtENgJdAFCgPVAH6A/zsRQ8RVX4bx3gZaNKXbgPuAO17lzG1nsAa7z2gFvNbLYzwWuwfkP16TGFLvrnEuAJcB19RW7N+N3nfcEMKiRxl5v/6/W8XPcDyS7jnm7pmsH0Qio6mIRSTxt81nADlXdBSAi7wCXquojQKVdCSKSAOSo6nEfhnsKb8QuIhlAkettqe+iPZW3vneXo0CoL+KsjJe+97FAC5z/YxWIyHxVLfNp4Hjve1fVj4GPReQz4G3fRfyz+3rjuxfgUWCBqq7xbcQ/8fLfeb/x5HPg7AGIB9bhRu9To0gaVegEpFd4nwEMreGcqcBrPovIfZ7GPg94VkTOBhb7MjA3eBS7iFwOXAC0Ap7zbWg18ih2Vf0zgIhMAQ7XR8Kohqff+xic3Q6hwHyfRuYeT//O342zpRctIt1U9SVfBlcDT7/7GOBhYKCI3O9KLg1BVZ/jGeA5EbkIN0qNNOakIZVsq/ZJRVV9wEexeMqj2FX1BM6E1xB4Gvs8nEmvIfD47wyAqs72fige8/R7XwQs8lUwteBp/M/g/MesIfA09iPAdN+FU2uVfg5VzQducfcijWIgvAoZwBkV3scD+/0Ui6csdv+w2P2nMcffmGOvyCufozEnjZVAdxFJEpEQnAOWH/s5JndZ7P5hsftPY46/McdekXc+h79G9z2cCTAHOMBPU06nurZPBH7EOSPgz/6O02JvOC+L3eJvbrHX1+ewgoXGGGPc1pi7p4wxxtQzSxrGGGPcZknDGGOM2yxpGGOMcZslDWOMMW6zpGGMMcZtljSM8SERyfN3DMZ4kyUNY+qZiAT6OwZjasuShjH1QJwrGC4UkbeBjf6Ox5jaasxVbo1pbM4C+qnqbn8HYkxtWUvDmPqzwhKGaewsaRhTf/L9HYAxdWVJwxhjjNssaRhjjHGblUY3xhjjNmtpGGOMcZslDWOMMW6zpGGMMcZtljSMMca4zZKGMcYYt1nSMMYY4zZLGsYYY9xmScMYY4zb/j8RN7Rdkb4EMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "batch_size = 256\n",
    "lr_finder = model.lr_finder(x_train, y_train, batch_size, tolerance=4)\n",
    "_ = lr_finder.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04229242874389523"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr_finder.get_best_lr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often, this learning rate is a little high, so we instead set it manually to 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.optimizer.set_lr(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We include the `EarlyStopping` callback to stop training when the validation loss stops improving. After training, this callback will also load the best performing model in terms of validation loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:\t[0s / 0s],\t\ttrain_loss: 1.6639,\tval_loss: 1.4607\n",
      "1:\t[0s / 0s],\t\ttrain_loss: 1.5303,\tval_loss: 1.4279\n",
      "2:\t[0s / 0s],\t\ttrain_loss: 1.4937,\tval_loss: 1.4188\n",
      "3:\t[0s / 0s],\t\ttrain_loss: 1.4411,\tval_loss: 1.4113\n",
      "4:\t[0s / 0s],\t\ttrain_loss: 1.4168,\tval_loss: 1.3912\n",
      "5:\t[0s / 0s],\t\ttrain_loss: 1.3792,\tval_loss: 1.3784\n",
      "6:\t[0s / 0s],\t\ttrain_loss: 1.3632,\tval_loss: 1.3795\n",
      "7:\t[0s / 0s],\t\ttrain_loss: 1.3575,\tval_loss: 1.3768\n",
      "8:\t[0s / 0s],\t\ttrain_loss: 1.3495,\tval_loss: 1.3730\n",
      "9:\t[0s / 0s],\t\ttrain_loss: 1.3305,\tval_loss: 1.3724\n",
      "10:\t[0s / 0s],\t\ttrain_loss: 1.3109,\tval_loss: 1.3740\n",
      "11:\t[0s / 0s],\t\ttrain_loss: 1.3149,\tval_loss: 1.3700\n",
      "12:\t[0s / 0s],\t\ttrain_loss: 1.3116,\tval_loss: 1.3695\n",
      "13:\t[0s / 0s],\t\ttrain_loss: 1.2873,\tval_loss: 1.3767\n",
      "14:\t[0s / 0s],\t\ttrain_loss: 1.2961,\tval_loss: 1.3831\n",
      "15:\t[0s / 0s],\t\ttrain_loss: 1.2934,\tval_loss: 1.3877\n",
      "16:\t[0s / 0s],\t\ttrain_loss: 1.2987,\tval_loss: 1.3941\n",
      "17:\t[0s / 0s],\t\ttrain_loss: 1.2882,\tval_loss: 1.3949\n",
      "18:\t[0s / 0s],\t\ttrain_loss: 1.2891,\tval_loss: 1.3841\n",
      "19:\t[0s / 0s],\t\ttrain_loss: 1.2712,\tval_loss: 1.3858\n",
      "20:\t[0s / 0s],\t\ttrain_loss: 1.2582,\tval_loss: 1.3860\n",
      "21:\t[0s / 0s],\t\ttrain_loss: 1.2768,\tval_loss: 1.3951\n",
      "22:\t[0s / 0s],\t\ttrain_loss: 1.2572,\tval_loss: 1.3902\n"
     ]
    }
   ],
   "source": [
    "epochs = 100\n",
    "callbacks = [tt.callbacks.EarlyStopping()]\n",
    "log = model.fit(x_train, y_train, batch_size, epochs, callbacks, val_data=val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVdrA8d+TZNILKZQQSKEEkA6hKArYEBDFwir2+rK2tezq6jYLq6+6xXV9Lei6LGvD3juiiIiUgPQSWiChpkAq6ef9497AEFImySQzSZ7v5zOfuXPvmXuf3EyeuTnn3HPEGINSSqn2y8fTASillGpZmuiVUqqd00SvlFLtnCZ6pZRq5zTRK6VUO+fn6QBqExMTYxITEz0dhlJKtRmrVq3KNsZ0rm2bVyb6xMREUlNTPR2GUkq1GSKyu65tWnWjlFLtnCZ6pZRq5zTRK6VUO+eVdfRKqfanvLyczMxMSkpKPB1KmxYYGEiPHj1wOBwuv0cTvVKqVWRmZhIWFkZiYiIi4ulw2iRjDDk5OWRmZpKUlOTy+7TqRinVKkpKSoiOjtYk3wwiQnR0dKP/K9JEr5RqNZrkm68p59ArE31ZRZWnQ1BKqXbDKxN9YWmFp0NQSql2wysTfVGZJnqllHsdOXKE559/vtHvmzp1KkeOHGn0+66//nrefffdRr+vJXhnoi+t9HQISql2pq5EX1lZf775/PPP6dSpU0uF1Sq8sntleWUV+/OOEhsR5OlQlFIt4JFPNrJpX75b93lK93AeumBgndsfeOABduzYwbBhw3A4HISGhhIbG8uaNWvYtGkTF110ERkZGZSUlHDXXXcxa9Ys4PjYW4WFhUyZMoXTTz+dpUuXEhcXx0cffURQUMN5auHChdx7771UVFQwatQoXnjhBQICAnjggQf4+OOP8fPzY9KkSfztb3/jnXfe4ZFHHsHX15eIiAgWL17c7HPTYKIXkbnANOCQMWZQHWUmAk8DDiDbGDPBXp8OFACVQIUxJsXVwFLTD3PBUE30Sin3eOKJJ9iwYQNr1qxh0aJFnH/++WzYsOFYf/S5c+cSFRXF0aNHGTVqFJdeeinR0dEn7GPbtm3Mnz+ff/3rX1x22WW89957XH311fUet6SkhOuvv56FCxeSnJzMtddeywsvvMC1117LBx98wJYtWxCRY9VDs2fP5quvviIuLq5JVUa1ceWKfh7wLPBKbRtFpBPwPDDZGLNHRLrUKHKmMSa7MUH5iJCanssFQ7s35m1KqTaivivv1jJ69OgTbjp65pln+OCDDwDIyMhg27ZtJyX6pKQkhg0bBsDIkSNJT09v8Dhbt24lKSmJ5ORkAK677jqee+457rjjDgIDA7n55ps5//zzmTZtGgDjxo3j+uuv57LLLuOSSy5xx4/acB29MWYxkFtPkSuB940xe+zyh5obVLC/LyvTDzd3N0opVaeQkJBjy4sWLeKbb77hp59+Yu3atQwfPrzWm5ICAgKOLfv6+lJR0XDHEWNMrev9/PxYsWIFl156KR9++CGTJ08GYM6cOTz66KNkZGQwbNgwcnJyGvujncQdjbHJQKSILBKRVSJyrdM2A3xtr59V305EZJaIpIpIKhWlbDmQT35JuRvCU0opCAsLo6CgoNZteXl5REZGEhwczJYtW1i2bJnbjtu/f3/S09PZvn07AK+++ioTJkygsLCQvLw8pk6dytNPP82aNWsA2LFjB2PGjGH27NnExMSQkZHR7Bjc0RjrB4wEzgaCgJ9EZJkxJg0YZ4zZZ1fnLBCRLfZ/CCcxxrwEvATQf/AwU2Lg5z1HmJBc64QpSinVKNHR0YwbN45BgwYRFBRE165dj22bPHkyc+bMYciQIfTr14+xY8e67biBgYH85z//4Re/+MWxxthbbrmF3Nxcpk+fTklJCcYY/vGPfwBw3333sW3bNowxnH322QwdOrTZMUhd/1acUEgkEfi0tsZYEXkACDTGPGy//jfwpTHmnRrlHgYKjTF/a+h4I0aONHnn/ZnbJvbmN5P6ufBjKKW83ebNmxkwYICnw2gXajuXIrKqrg4v7qi6+Qg4Q0T8RCQYGANsFpEQEQmzAwgBJgEbXNmhjwgDu4ezMr2+pgGllFKucKV75XxgIhAjIpnAQ1jdKDHGzDHGbBaRL4F1QBXwsjFmg4j0Aj6wB+DxA94wxnzpamApCVG8sWI3ZRVV+Pt55X1dSinF7bffzo8//njCurvuuosbbrjBQxGdrMFEb4y5woUyfwX+WmPdTqDJlUujEiOZ++MuNu7LY3h8ZFN3o5RSLeq5557zdAgN8tpL5ZGJVnJP1W6WSinVLF6b6LuEBZIYHaz19Eop1Uxem+gBUhKjSN19uM4bDpRSSjXMqxP9qMRIcovK2Jld5OlQlFKqzfLqRJ+SGAVAqlbfKKVaWWhoaJ3b0tPTGTSo1jEevZJXJ/peMSFEhfjruDdKKdUMXjkefTURISUhUq/olWpvvngADqx37z67DYYpT9S5+f777ychIYHbbrsNgIcffhgRYfHixRw+fJjy8nIeffRRpk+f3qjDlpSUcOutt5Kamoqfnx9PPfUUZ555Jhs3buSGG26grKyMqqoq3nvvPbp3785ll11GZmYmlZWV/OlPf+Lyyy9v1o/tCq9O9ACjEqP4etNBDhWU0CUs0NPhKKXaqJkzZ3L33XcfS/Rvv/02X375Jffccw/h4eFkZ2czduxYLrzwQuwbPV1S3Y9+/fr1bNmyhUmTJpGWlsacOXO46667uOqqqygrK6OyspLPP/+c7t2789lnnwHWYGqtwesTfYrdn35V+mGmDI71cDRKKbeo58q7pQwfPpxDhw6xb98+srKyiIyMJDY2lnvuuYfFixfj4+PD3r17OXjwIN26dXN5v0uWLOFXv/oVYI1UmZCQQFpaGqeeeiqPPfYYmZmZXHLJJfTt25fBgwdz7733cv/99zNt2jTOOOOMlvpxT+DVdfQAA7tHEOjw0Xp6pVSzzZgxg3fffZe33nqLmTNn8vrrr5OVlcWqVatYs2YNXbt2rXUc+vrU1f37yiuv5OOPPyYoKIjzzjuPb7/9luTkZFatWsXgwYP53e9+x+zZs93xYzXI66/o/f18GNazE6m7tZ5eKdU8M2fO5H/+53/Izs7m+++/5+2336ZLly44HA6+++47du/e3eh9jh8/ntdff52zzjqLtLQ09uzZQ79+/di5cye9evXizjvvZOfOnaxbt47+/fsTFRXF1VdfTWhoKPPmzXP/D1kLr0/0YNXTP79oB0WlFYQEtImQlVJeaODAgRQUFBAXF0dsbCxXXXUVF1xwASkpKQwbNoz+/fs3ep+33XYbt9xyC4MHD8bPz4958+YREBDAW2+9xWuvvYbD4aBbt248+OCDrFy5kvvuuw8fHx8cDgcvvPBCC/yUJ3NpPPrWlpKSYlJTU4+9/j4ti+vmruD1m8cwrk+MByNTSjWVjkfvPp4Yj77FjYjvhI+g494opVQTtIl6kLBAB/27hetIlkqpVrV+/XquueaaE9YFBASwfPlyD0XUNK5MPDIXmAYcqm0qQbvMROBprAlJso0xE+z1k4F/Ar5YE5I0uU/VqMRI3lmVSUVlFX6+beIfEaVUDcaYRvVR97TBgwcfm7TbWzSlut2VjDkPmFzXRhHpBDwPXGiMGQj8wl7vCzwHTAFOAa4QkVMaHaEtJTGK4rJKNu+vfRZ3pZR3CwwMJCcnR0ejbQZjDDk5OQQGNu7mUVdmmFpsTw5elyuB940xe+zyh+z1o4Ht9kxTiMibwHRgU6MitFXfOLUyPZfBPSKasgullAf16NGDzMxMsrKyPB1KmxYYGEiPHj0a9R531NEnAw4RWQSEAf80xrwCxAEZTuUysSYOr5WIzAJmAcTHx5+0PTYiiB6RQaTuzuXG05PcELZSqjU5HA6SkvRv1xPckej9gJHA2UAQ8JOILANqq4ir8382Y8xLwEtgda+srcyoxCiWbM9uc/V8SinlSe5o1cwEvjTGFBljsoHFWJOCZwI9ncr1APY150ApiZFkFZSyJ7e4ObtRSqkOxR2J/iPgDBHxE5FgrOqZzcBKoK+IJImIPzAT+Lg5BxplT0Si494opZTrXOleOR+YCMSISCbwEFY3Sowxc4wxm0XkS2AdUIXVjXKD/d47gK+wulfONcZsbE6wfTqHEhHkIDU9lxkjG9cYoZRSHZUrvW6ucKHMX4G/1rL+c+DzpoV2Mh8fayISvUNWKaVc1+buPEpJjGJHVhE5haWeDkUppdqENpfoR1VPRLJb6+mVUsoVbS7RD+4Rgb+fD6ma6JVSyiVtLtEH+PkytEeE1tMrpZSL2lyiB6uefsPePI6WVXo6FKWU8nptMtGPSoykvNKwNvOIp0NRSimv1yYT/ch468apVK2+UUqpBrXJRB8R7KBf1zC9Q1YppVzQJhM9WOPerN59mMoqHdtaKaXq02YT/ajEKApKK9h6QCciUUqp+rTZRF89EUnqbq2nV0qp+rTZRB/XKYjYiECtp1dKqQa02UQvIqQkRrFyV67OQamUUvVos4kerP70B/JL2HvkqKdDUUopr9WmE31KQnV/eq2+UUqpujSY6EVkrogcEpENdWyfKCJ5IrLGfjzotC1dRNbb61PdGThAv25hhAX46bg3SilVD1cmB58HPAu8Uk+ZH4wx0+rYdqY9l6zb+foIIxIi9YpeKaXq0eAVvTFmMeC1l8yjEiPZerCAvOJyT4eilFJeyV119KeKyFoR+UJEBjqtN8DXIrJKRGa56VgnSLEnDF+1x2u/i5RSyqPckehXAwnGmKHA/wEfOm0bZ4wZAUwBbheR8XXtRERmiUiqiKRmZWW5fPChPTrh8BXtT6+UUnVodqI3xuQbYwrt5c8Bh4jE2K/32c+HgA+A0fXs5yVjTIoxJqVz584uHz/I35dBcRE6kqVSStWh2YleRLqJiNjLo+195ohIiIiE2etDgElArT13mmtUYhRrM/IoKdeJSJRSqiZXulfOB34C+olIpojcJCK3iMgtdpEZwAYRWQs8A8w01q2qXYEl9voVwGfGmC9b4odISYikrLKKDXvzWmL3SinVpjXYvdIYc0UD25/F6n5Zc/1OYGjTQ3PdyARrgLOV6YePNc4qpZSytOk7Y6tFhwbQu3OI1tMrpVQt2kWiB6uePnX3Yap0IhKllDpBu0n0KYlR5B0tZ3tWoadDUUopr9JuEv0oeyKS5TtzPByJUkp5l3aT6OOjgunXNYynFqSRnl3k6XCUUsprtJtELyK8eM1IAG6ct5IjxWUejkgppbxDu0n0AIkxIbx0bQqZh4/yy1dXUVZR5emQlFLK49pVoger981ffzGE5btyeeD9dTrNoFKqw3NlPPo2Z/qwONKzi/nHN2kkRodw59l9PR2SUkp5TLtM9AB3nt2H3TlFPLUgjYToYKYPi/N0SEop5RHtNtGLCI9fOpjMI0e57511dO8UxCgdHkEp1QG1uzp6ZwF+vrx49UjiIoOY9UqqdrtUSnVI7TrRA0SG+DP3+lEYtNulUqpjaveJHiApJoSXrtFul0qpjqlDJHqA0UlR/GWGdrtUSnU87bYxtjYXDY8jPaeIp7/ZRlJ0CL/SbpdKqQ7AlRmm5orIIRGpdRpAEZkoInkissZ+POi0bbKIbBWR7SLygDsDb6q7zu7LxcPj+PuCND5as9fT4SilVItzpepmHjC5gTI/GGOG2Y/ZACLiCzwHTAFOAa4QkVOaE6w7iAhPXDqY0UlR3PfOOlbqZCVKqXauwURvjFkMNCUbjga2G2N2GmPKgDeB6U3Yj9tpt0ulVEfirsbYU0VkrYh8ISID7XVxQIZTmUx7Xa1EZJaIpIpIalZWlpvCqpt2u1RKdRTuSPSrgQRjzFDg/4AP7fVSS9k6u7oYY14yxqQYY1I6d+7shrAapt0ulVIdQbMTvTEm3xhTaC9/DjhEJAbrCr6nU9EewL7mHs/dRidF8eSMwSzflcvDn2z0dDhKKeV2zU70ItJNRMReHm3vMwdYCfQVkSQR8QdmAh8393gt4eLhPfifM5J4Y/keVu3WxlmlVPviSvfK+cBPQD8RyRSRm0TkFhG5xS4yA9ggImuBZ4CZxlIB3AF8BWwG3jbGeO0l893nJNMtPJAHP9pIZZXeTKWUaj/EG+8QTUlJMampqa1+3E/W7uNX83/m0YsGcfXYhFY/vlJKNZWIrDLGpNS2rcMMgeCKaUNiGdsrir9+tZXcIu2Fo5RqHzTROxERHrlwEIWlFfz1q62eDkcppdxCE30N/bqFcd2piby5cg/rMo94OhyllGo2TfS1uPvcvkSHBPDgRxup0oZZpVQbp4m+FuGBDn43pT9rMo7w7upMT4ejlFLNoom+DhcPj2NkQiRPfrGFvKPlng5HKaWaTBN9HXx8hEcuHMjh4jL+sSDN0+EopVSTaaKvx6C4CK4ak8ArP6WzeX++p8NRSqkm0UTfgN9MSiYiyMFDH23U6QeVUm2SJvoGdAr257eT+7MiPZeP13rdmGxKKdUgTfQuuCylJ0N6RPDYZ5spLK3wdDhKKdUomuhd4Gs3zB4qKOWZhds8HY5SSjWKdyb68mJPR3CS4fGRXJ7Sk7lLdrH9UIGnw1FKKZd5Z6LP2QEHN3k6ipP8dnI/gv19efjjTdowq5RqM7wz0YvAK9OthO9FokMD+M2kfizZns2XGw54OhyllHKJKxOPzBWRQyKyoYFyo0SkUkRmOK2rFJE19sP12aWi+4Cpgv9eCEf2uPy21nDVmHj6dwvjz59u4mhZpafDUUqpBrlyRT8PmFxfARHxBZ7Emk3K2VFjzDD7caHLUfkFwrUfQlkB/PcCyN/v8ltbmp+vD7OnD2JfXgnPL9ru6XCUUqpBDSZ6Y8xioKGJVH8FvAccckdQAHQbDFe/D0XZVjVOUbbbdt1co5OiuHh4HC9+v5P07CJPh6OUUvVyx+TgccDFwJxaNgeKSKqILBORixrYzyy7bGpWVpa1skcKXPm2VX3z6kVw9HBzw3Wb303pj8NXmP2p9zUaK6WUM3c0xj4N3G+Mqa3COt6ew/BK4GkR6V3XTowxLxljUowxKZ07dz6+IXEczHwNsrbCazOg1Du6NnYJD+Tuc5L5dsshFm4+6OlwlFKqTu5I9CnAmyKSDswAnq++ejfG7LOfdwKLgOFNOkKfc+AX82Dfz/DGTCjzjn72149LpE+XUB75ZBMl5dowq5TyTs1O9MaYJGNMojEmEXgXuM0Y86GIRIpIAICIxADjgKbXc/Q/Hy55CXb/CG9dDRWlzQ292Ry+Pjxy4UD25Bbzr8U7PR2OUkrVyq+hAiIyH5gIxIhIJvAQ4AAwxtRWL19tAPCiiFRhfaE8YYxpXoX24BlQfhQ+vgPevdG6yvd1NGuXzTWuTwznD47lnwu3ERLgxw3jEhERj8aklFLOGkz0xpgrXN2ZMeZ6p+WlwOCmhVWPEddYQyR88Vv48Fa4+EXw8XX7YRrj8UsHU1ZZxexPN/FzxhGeuGQwIQENnlqllGoV3nlnbEPG/BLOfgjWvwOf3g0eHo4gPNDBi1eP5LeT+/HZun1c9NyPbD9U6NGYlFKqWttM9ABn/BrG3werX4EvH/B4svfxEW6b2IdXbxpDblEZ059dwufrvedGL6VUx9V2Ez3AmX+AsbfD8jmwcLanowGsOvtP7zyd5G5h3Pb6ah77bBPllVWeDksp1YG17UQvAuc9BiNvgCVPweK/eToiAGIjgnhr1qlcd2oC//phF1f9azmHCko8HZZSqoNq24kerGR//lMwZCZ8+2er6+WaN6DQfaMxNIW/nw+PTB/E05cPY/3ePM5/ZgkrdjU0koRSSrmfeOO46ikpKSY1NbVxb6qsgIUPw7p3oNAeQrj7CEg+D/qeC7HDwccz32tbDuRz62ur2ZNbzO+m9Oem05O0C6ZSyq1EZJU9EsHJ29pNoq9mDBxYB2lfw7avIDMVMBDSxUr4fSdB7zMhMMKtMTckv6Sc+95Zy1cbD3L+4FienDGEUO2CqZRyk46V6GsqyoHt31hJf/tCKDkCPn4Qf6qV9PtOgs79rCqgFmaM4aXFO3nyyy0kxYTw4jUj6dMlrMWPq5Rq/zp2ondWWQGZK62kv20BHLTnUukUD33Pg95nQfxYCI5y/7Gd/LQjh1/NX01xWSV/mTGEaUO6t+jxlFLtnyb6uuRlwravraS/c9HxSck7D4CEUyH+NOs5oofbD30gr4Tb31jNqt2HuXFcEn84fwC+Plpvr5RqGk30rigvgX2rYfdS65GxwprhCiAi3k78p0LCaRCT7JaqnvLKKh77bDPzlqZz8+lJ/HHaKc3ep1KqY6ov0WtrYDVHoJXEE06zXldVwoH1sOcnK/Hv+BbWvWVtC462kn78qdYXQLeh4Nv4U+nw9eHhCwdijOHlJbvo2zWUy0fFu/GHUkopTfR18/GF7sOsx9hbrd48OTtgz1LY/ZP1vOVTq6wjxOrJc+5siK5zbpU6/WnaKezMLuKPH24gMTqEMb2i3fzDKKU6Mq26aY78/ccT/7q3oLIczvw9jL2t0Vf4eUfLufj5HzlcVMZHt59OfHRwCwWtlGqP6qu6aft3xnpSeCwMuhTO/xvcvty6ql/wJ3j5bNi/rlG7ighy8O/rRlFl4Kb/rqSgpLyFglZKdTQuJXoRmSsih0RkQwPlRolIpYjMcFp3nYhssx/XNTdgrxXeHWa+ATP+A/l74aWJ8M0jViOvi5JiQnjhqhHsyi7izvk/U1nlff9tKaXaHlev6OcBk+srICK+wJPAV07rorBmpBoDjAYeEpHIJkXaFojAoEvg9hUwdKY10NqccVZjrotO6xPDI9MH8t3WLB7/fHMLBquU6ihcSvTGmMVAQyNy/Qp4D3AeTew8YIExJtcYcxhYQANfGO1CcBRc9Dxc8wFUlsF/psCn90BJvktvv2pMAteflsjLS3bx1so9LRysUqq9c0sdvYjEARcDNeeQjQMynF5n2utq28csEUkVkdSsrCx3hOV5vc+C25ZZY+avmgfPjYGtX7j01j+eP4Az+sbwxw83sGxnTsvGqZRq19zVGPs0cL8xprLG+truKqq14tkY85IxJsUYk9K5c2c3heUF/ENg8v/CTd9AUCeYPxPeuQEK6/8y8/P14dkrR9AzKphbX1vFnpziVgpYKdXeuKsffQrwpj30bgwwVUQqsK7gJzqV6wEsctMx25YeI2HW9/DjP2HxX2Dnd3De41Zdfh132Vb3xLnouR+56b8ree+20wgPdLRy4Eo1Qf5+607zvL1W9WVlmdX9uLL0+HJFqb2uzF5vL1ev9w+B/ufDKRdBSDu+tyRnhz0Uy9eQvw+i+0BMX+sO/Jhk63VQp2YdwuV+9CKSCHxqjBnUQLl5drl37cbYVcAIe/NqYKQxpt76/jbTj76psrbCx3dCxjKremfa0xCZUGfxpTuyufbfKzi9bwz/vm6UjomjvEtxrpXU9/5sP68+PidETT5+4OsPvg77OcBp2R/8/I9vz98POdtAfK2uy4MutRJ/Kw8x7nblJbD7x+PJPXentT6mn5XUc3dYyb/KqYt1aFc78fc98Tm8x7F5Npo9BIKIzMe6Mo8RkUysnjQOAGNMzXr5Y4wxuSLyZ2ClvWp2Q0m+Q+jcD274AlL/Dd88bNXdD7kMRt0EsUNPKn5ab6snzh8+2MD/fr6ZP+mYOMpTSgtg/1ormVcn9SO7j2+P7gu9JliT/sSNgMgk8As4nrx9fF0/ljFwcCNseBc2vAcf3mp9MfQ910r6yZPBv43cWHgk4/gAiru+twZQ9AuEpPHWDZZ9z4XIxOPlKyus85qd5vTYBhvet4Zar+YXBDF9rKRfD70z1tPyMmHRE7D+Xag4CnEpkHKj1U3TEXRC0Yc/3si8pek8cclgZo7WMXFUCzIGinMgZ7t18191Us9O41gzW0Q8xA23knr34dZwIS11tW0M7F1l/Z1s/MD6j8ERAv2mwOAZ1n/GfgEtc+ymqCyHjOXHk/uhTdb6Tgn2rHeTIPH0k/7GG1T9e3FO/vay3L1OR6/0ekePwNo3rav87DQI7ATDrrKSfkwfACoqq7jxv6ks3Z7NazePYayOiaOa6+iR41UFOTvs5e2QsxNK846XC+l8/Cq9OrGHeqjTRFWldW/Khvdg00dwNNf6ghlwgXWlnzi+4SFIjLHaAkoLoDS/xnMBVFWAqbLKmSrrgXF6XXN99bpK64txx3fW+fNxWAMlVk9yFNO3xSY50mGK2xJjIH2JlfA3f2J94JImWNU6/aaSVwaXPP8juUVlfHj7OBKiQzwdsfJ2ZUUnJ/Hq5WLnrrsCET0hupdVVxzV2xqkr+tACI9rlVnYGq2y3JpLYsN7sPlTa2jxkM5WtY6v43jiLi2w7mM5IZm30DAjYbHHpy1NmgCB4S1znBo00bdVBQfh51etPvh5GRDaDUZcS0avy7jglV3EhAbwvvbEUc7Kiq3htff9bD9WW//eO/dqDou1k3iNhB6ZZA3X3VaVl1hVJRves3q1+fpDQJj9CHdarrmulm2+DkBAfOyHvXxsnRxfV3O9X6BHvhQ10bd1VZVWPV/qv61nEXLizuLeXSmYXmfy8vWj8fOt5ZYIY6CiBMqPHn9UVC8XW38YjiDrCii0CwRFHWvBV21ARak1Hea+n+0eLz9D1ma7OgGrp0b3EVbdeed+VkKP6gUBoZ6NW7UITfTtyeF06wp/9atQnE1GVWfyHTH0DIMw33KkvMRO4nZSbwzxhZAYK/FXJ/8TlrtY26vX++p/Eq2mstxq0Dt2pf4zHNx0vPohONqpUdR+hMd6NmbVqjTRt0cVpbD5Ew4seYW9WYfJLfcjMCiUPnExdIuOQvyDrK5XjiBwBFv/kjuCrdfH1gdaXwiFh6Ao6/jzseVD1h28dX1hBEdbdbqdelo9MDr1PP66U7zVoOyN9bptQXGu1WtjzzLree9q66YisBoenRN69xHWvMZ6rjs0TfTtXEVlFe+v3ss/F25j75GjjEyI5DfnJnNanxj3HKC08HjSL8o6vlywz+ofnJdhPdf8QvAPOzH5O38phHWz6lB9/Ky+1b4Oe9nveF2nKyrLrYa1skKr0bG00GqQKy201mkmsLUAABZTSURBVNV8XVYEQZFWn+VOCdaNap3iG9/NzZ2qZy/LWHY8sWenWdt8HFbVS88xdo+X4VZduiZ1VYMm+g6irKKKt1MzePbb7RzIL2Fsryh+fW4/RidFtfzBq/v3HtlzPPEfe95jPTvf6NGQ6qTvY99kc+y1/UVQbif16qtcV/bnH2rdVl+cY7VdOAvtZid9O/k7fxGExzXuRp+GVJRaNx1VJ/U9y6A429oW2MlK6vFjrUf34Z79ElJthib6DqakvJL5K/bw3Hc7yC4s5Yy+Mfz63GSGx3t4KoCSfOsGsbwMKDhg1S9XVVpdSKsflRUnvq6qtMs5v660kl9AqNVjwj/UWvYPtXpM1PbaL+D4VbAxUHgQDu+27j48vNtq+6hezs883qAJ1pdERA8r+QdFOX3p+Nb+JeT82rlMUfbJ1TBRvaDnWIgfYz3HJGuDuGoSTfQd1NGySl5btpsXvt9BblEZZ/Xvwq/PTWZQXBsfK6SlVZZbX0jOyb/6uSTPuinm2JeO85dSldNyLX20fRzWEBfVV+s9x1gN20q5gSb6Dq6otIJ5S9N5afFO8o6Wc97ArtxzbjL9u7XOjRwdVlWV05dChT1olxfdpq/aFU30CoD8knLmLtnFv3/YRUFpBecPieXylJ4MiA2nc5gmIKXaMk306gRHisv41w87+c+P6RSXWXPFxIT6079bOANiw+jfLZz+sWH06RJKgJ8bGyGVUi1GE72qVX5JORv25rF5fwFb9uez5UABWw8WUFZhNUT6+Qi9O4fS3yn5nxIbTpewAES79ynlVZo9Hr1qn8IDHZzWO4bTeh/vb19RWUV6TjGb9+ez5UA+W/YXkJp+mI/W7DtWJjLYQf9u4QzpEcHY3tGMTowiJEA/Skp5K72iVy7JKy63Ev+BArYcyGfT/gI278unrLIKPx9hWM9OnNY7mlN7xzA8vhOBDq3yUao1NavqRkTmAtOAQ7VNIygi04E/A1VABXC3MWaJva0SWG8X3WOMudCVgDXRtw1HyypZtfswS3dks3RHDusyj1BlIMDPh5TESE7rHcOpvaMZEhdR+6BrSim3aW6iHw8UAq/UkehDgSJjjBGRIcDbxpj+9rZCY0yjh8rTRN825ZeUs3JXLkt35LB0Rw6b9+cDEBrgx+ikKPuKP5oB3cLx0XlvlXKrZtXRG2MW2xOD17W90OllCCcMfK06kvBAB2cP6MrZA7oCkFNYyvJdufy4PZufduTw7ZZDAHQKdjAhuTO/PjdZJ05RqhW4pQVNRC4GHge6AOc7bQoUkVSsKp0njDEf1rOPWcAsgPh4nQ+1PYgODWDq4FimDraGy92fd5Sf7Kv9Lzcc4IsNB7h1Qm9undhb6/SVakEuNcbaV/Sf1lZ1U6PceOBBY8w59uvuxph9ItIL+BY42xizo6HjadVN+3cwv4RHP9vMJ2v3kRAdzMMXDuTMfjocgFJNVV/VjVtbyIwxi4HeIhJjv95nP+8EFgHD3Xk81XZ1DQ/k/64Yzus3j8HXR7jhPyv55aup7D3SyMlSlFINanaiF5E+Yt89IyIjAH8gR0QiRSTAXh8DjAM2Nfd4qn0Z1yeGL+8az33n9eP7tCzO+fv3vLBox7GbtpRSzddgHb2IzAcmAjEikgk8BDgAjDFzgEuBa0WkHDgKXG73wBkAvCgiVVhfKE8YYzTRq5P4+/lw+5l9mD6sO7M/2cSTX27hvdWZzJ4+8ISbuZRSTaM3TCmv8+2Wgzz08UYyco8yfVh3/jB1AF3CAz0dllJerdXq6JVyh7P6d2XBPRO486w+fLH+AGf//XvmLtlFRaVW5yjVFJrolVcKdPjy60n9+Oqe8QxPiGT2p5u44NkfWbX7sKdDU6rN0USvvFpSTAj/vWEUz181gsNFZVz6wlJ+++5acovKPB2aUm2GJnrl9USEqYNjWfibCfxyfC/eX72Xs/++iPdWZeKNbUxKeRtN9KrNCAnw43dTB/DpnaeTFBPCb95Zy1UvL2dnVmHDb1aqA9NEr9qc/t3CefeW03j0okGs35vH5H/+wDMLt1FaUenp0JTySproVZvk4yNcPTaBhb+ewKRTuvLUgjSm/vMHVuzK9XRoSnkdTfSqTesSHsizV47gPzeMorSiiste/In7313HkWJtrFWqmiZ61S6c2a8LX98znl9O6MW7qzM5++/f8+HPe7WxVik00at2JNjfj99NGcAnd5xOz6hg7n5rDdf8ewXp2UWeDk0pj9IhEFS7VFlleGP5bv7y5VZKK6u486w+zBrfG3+/+q9tjDEcLi4nI7eYPbnFZBwuJiP3KBn2Mlgjb8ZGBNItPJBuzs8RgXQODdBpE5VHNGsqQU/QRK/c5WB+CbM/2cRn6/fTt0so/3vJYAZ1jyDzsJ3Ic4vZk3vUTujWo6jsxN47USH+9IwKpmdkED4iHMgv4UCe9SirMSyDj0DnsAC6RQTRLTyA2IigY18Mg+LC6dMlrDV/fNWBaKJXHd7CzQd58KONtY53H+jwIT4qmJ6RwVZCt5N6fHQwPSKDCQ2ofZDX6qv//XlHOZhfwv68Eg7mWc/OXwYFpRWA9SXw28n9+eX4XtgjeyvlNprolQKKyyp45afdVFRW0TPKSuLxUcHEhPq3aOItLK3gQN5R/rFgG5+t38+kU7ryt8uGEh7oaLFjqo5HE71SXsAYw9wf03n88830iAzihatHMiA23NNhqXai2cMUi8hcETkkIhvq2D5dRNaJyBoRSRWR0522XSci2+zHdU37EZRq+0SEm05PYv6ssRSXVXLx8z/y/upMT4elOgBXuwfMAybXs30hMNQYMwy4EXgZQESisGakGgOMBh4SkcgmR6tUOzAqMYpP7zydoT068eu31/KHD9br8A2qRbmU6O1Jv+u8t9wYU2iO1wGFANXL5wELjDG5xpjDwALq/8JQqkPoEhbI6zeP4Zfje/H68j1c9uIynRhdtZgG54x1lYhcDDwOdAHOt1fHARlOxTLtdbW9fxYwCyA+Pt5dYSnltfx8ffjd1AEMj+/Eve+sY9ozP/DPmcMZn9zZ06E1Wk5hKWkHC9l2qIC0gwXW8sECAvx8mTyoG+cPiWVkfCQ+PtrbyBNcbowVkUTgU2PMoAbKjQceNMacIyL3AQHGmEftbX8Cio0xf69vH9oYqzqanVmF3PraatIOFXDPOcnccWafZiVFYwyb9xeweFsWK3bl4vAVokMDiAnxJzo0gOhQf6JDAogJtV53CnK4dLzDRWVWIj9kJfK0gwVsO1hIjtNEMGEBfvTtGkpy1zByi8pYlJZFWUUVXcMDmDIoVpN+C6mvMdZtV/TVjDGLRaS3iMRgXcFPdNrcA1jk7mMq1db16hzKB7efxu/fX89TC9L4ec9h/nH5MDoF+7u8j8NFZSzZns33aVksTsviUEEpAH26hOIjkJp+mNziMmq7tvMRiDqW+K0vgehQf6KC/ckpKmPboQK2Higku7D02HtC7YR+zoCu9O0aSt+uYSR3DaVbeOAJ3VULSytYuPkgn63bzxsr9jBvafqxpD91cCwpCZr0W5pbruhFpA+wwxhjRGQE8AlWUo8EVgEj7KKrgZHGmHrHktUretVRGWN4bfkeZn+yka7hgcy5eiSD4iJqLVtZZVibeYTvt2bxfVoW6zKPUGUgIsjBGX1jmJDcmfHJnekaHnjCew4Xl5FTWEZOYSk5Rcefs2uuKyyjoLSCEH9f+nQNI7mLdZVefbUeGxHY6PsPnJN+zSt9TfrN0+x+9CIyH+vKPAY4iNWTxgFgjJkjIvcD1wLlwFHgPmPMEvu9NwK/t3f1mDHmPw0dTxO96uh+3nOY215fTU5RGX+ePpDLR1ntVgfzS/g+zUrsS7Zlk3e0HBEY1rPTscQ+tEcnfN2ULEsrKnH4+LRI8q1O+p+v3893W+tO+sYYyiqrKKuwH07Lpfaj5npfH5jYrwuBDl+3x+2t9IYppdqgnMJS7npzDUu2ZzOxX2cO5JWw5UABAF3CAhif3JkJyZ05vU8MkSGuV/F4I+ekv2hrFqUVVQT4+WAMJ40n5KoekUH8fuoApgzq1iGGnNBEr1QbVVllePqbNOYtTWdQ9wgm9LOSe/9uYe02eRWWVvDtlkOszTiCw9cHfz8fAvx88LeX/Z2WA/yOr7PK+OLv58O+vKM8+cUWthwoYHRSFA9OO6XOKrD2QhO9UqrDqawyvLlyD3//Oo3DxWVcNrIn957Xj85hAZ4OrUU0ewgEpZRqa3x9hKvGJPDdvRO5aVwS763O5My/LWLO9zs63J3ImuiVUu1aRJCDP047ha/vGc/YXlE88cUWzn1qMV9uONBhpprURK+U6hB6dQ7l5etG8epNowl0+HDLa6u48l/L2bw/39OhAZCeXcSynTktsm+to1dKdTgVlVXMX7GHpxakkXe0nMtHxfObScnEhLZ+/X1RaQXPfredl3/YSXmlYfLAbjx84UC6RQQ2/GYn2hirlFK1yCsu5+mFabz6026CHL7ceXZfrjstscG5hd3BGMNn6/fz2Geb2Z9XwoyRPUiMDubZ77bj5+PDfef14+qxCS7fE6GJXiml6rH9UCGPfbaJ77ZmkRgdzO1n9uGCod1b7IartIMFPPTRRn7amcPA7uHMnj6QkQlRAOzJKeYPH67nh23ZDO3ZiccvHswp3RueoEYTvVJKuWDR1kM8/vkWth4sICLIwWUpPbhqTAKJMSFu2X9BSTlPf7ONeUvTCQ3w477z+nHF6PiTrtqNMXy8dh9//nQTh4vLufn0JO46py/B/nUPT6aJXimlXGSMYdnOXF5btpuvNh6gosowPrkz14xN4Kz+XZo0vIQxhg9+3sv/fr6FnKJSZo6K577z+hHVwB3NR4rLeOKLLby5MoMekUH8+aJBnNmvS61lNdErpVQTHMwv4c0VGbyxYjcH80uJ6xTElWPiuSylp8s3Xm3cl8dDH20kdfdhhvXsxOzpAxnSo1Oj4lixK5fff7Ce7YcKmTYklgcvOIUuYSc21mqiV0qpZqiorOKbzQd5ddluftyeg8NXmDIolmtOTSAlIbLW4Sjyisv5+4KtvLZsN52C/Xlgcn9mjOzR5AHiSisqefH7nTz73XYC/Hx4YEp/rhgVf2x/muiVUspNth8q5PXlu3l3VSYFJRX07xbG1WMTuGh4HKEBflRVGd5OzeAvX23lSHEZ14xN4Nfn9iMi2OGW4+/MKuQPH2zgp505jEyI5PFLBpPcNUwTvVJKuVtxWQUfr9nHq8t2s3FfPqEBfkwf1p0Ne/NYm5nHqMRIHrlwkEs9ZhrLGMN7q/fy2GebKCip4JcTevHbyQM00SulVEswxrAm4wivLtvNp+v2ExHk4PdT+3PRsLgWH2E0t6iMxz7bzHurM9n95LSmJ3oRmQtMAw7VMbvUVcD99stC4FZjzFp7WzpQAFQCFXUFUZMmeqVUW1RYWoHDVwjwa90JT5Zuz2Zc387NGr1yHjC5nu27gAnGmCHAn4GXamw/0xgzzNUkr5RSbVVogF+rJ3mA0/rE1Lu9wcnB7cm+E+vZvtTp5TKsuWKVUkp5CXcP6HAT8IXTawN8LSKrRGRWfW8UkVkikioiqVlZWW4OSymlOq4Gr+hdJSJnYiX6051WjzPG7BORLsACEdlijFlc2/uNMS9hV/ukpKR4XwuxUkq1UW65oheRIcDLwHRjzLEBlY0x++znQ8AHwGh3HE8ppZTrmp3oRSQeeB+4xhiT5rQ+RETCqpeBScCG5h5PKaVU4zRYdSMi84GJQIyIZAIPAQ4AY8wc4EEgGnje7jNa3Y2yK/CBvc4PeMMY82UL/AxKKaXq4Uqvmysa2H4zcHMt63cCQ5semlJKKXfQOWOVUqqd88ohEESkANjq6ThqEQNkezqIWmhcjaNxNY7G1TieiivBGNO5tg1u617pZlu98U5aEUnVuFyncTWOxtU4GpfrtOpGKaXaOU30SinVznlroq85MJq30LgaR+NqHI2rcTQuF3llY6xSSin38dYreqWUUm6iiV4ppdo5jyV6EZksIltFZLuIPFDL9gARecvevry+MfHdGFNPEflORDaLyEYRuauWMhNFJE9E1tiPB1s6Lqdjp4vIevu4J03BJZZn7HO2TkRGtEJM/ZzOxRoRyReRu2uUaZVzJiJzReSQiGxwWhclIgtEZJv9HFnHe6+zy2wTketaIa6/isgW+/f0gYh0quO99f7OWyCuh0Vkr9Pvamod763377cF4nrLKaZ0EVlTx3tb8nzVmh+84TPWIGNMqz8AX2AH0AvwB9YCp9Qocxswx16eCbzVCnHFAiPs5TAgrZa4JgKfeui8pQMx9WyfijUfgABjgeUe+L0ewLpxo9XPGTAeGAFscFr3F+ABe/kB4Mla3hcF7LSfI+3lyBaOaxLgZy8/WVtcrvzOWyCuh4F7Xfg91/v36+64amz/O/CgB85XrfnBGz5jDT08dUU/GthujNlpjCkD3gSm1ygzHfivvfwucLZIy860a4zZb4xZbS8XAJuBuJY8pptNB14xlmVAJxGJbcXjnw3sMMbsbsVjHmOsuQ5ya6x2/hz9F7iolreeBywwxuQaYw4DC6h/+sxmx2WM+doYU2G/9MjMbHWcL1e48vfbInHZOeAyYL67jueqevKDxz9jDfFUoo8DMpxeZ3JyQj1Wxv6DyMMaJbNV2FVFw4HltWw+VUTWisgXIjKwtWKi4Rm7XDmvLWkmdf8BeuqcdTXG7AfrDxXoUksZT5+3GzlxZjZnLs/S5kZ32FVKc+uohvDk+ToDOGiM2VbH9lY5XzXyg9d/xjyV6Gu7Mq/Zz9OVMi1CREKB94C7jTH5NTavxqqaGAr8H/Bha8RkG2eMGQFMAW4XkfE1tnvynPkDFwLv1LLZk+fMFZ48b38AKoDX6yjS0O/c3V4AegPDgP1Y1SQ1eex8AVdQ/9V8i5+vBvJDnW+rZV2r9W33VKLPBHo6ve4B7KurjIj4ARE07d/MRhERB9Yv8XVjzPs1txtj8o0xhfby54BDROqfgt1NTMMzdrlyXlvKFGC1MeZgzQ2ePGfAwerqK/v5UC1lPHLe7Aa5acBVxq7IrcmF37lbGWMOGmMqjTFVwL/qOJ6nzpcfcAnwVl1lWvp81ZEfvPYzVs1TiX4l0FdEkuwrwZnAxzXKfAxUt0zPAL6t64/BXez6v38Dm40xT9VRplt1W4GIjMY6hzm1lXVzbK7M2PUxcK1YxgJ51f9StoI6r7Q8dc5szp+j64CPainzFTBJRCLtqopJ9roWIyKTgfuBC40xxXWUafVZ2mq06Vxcx/Fc+fttCecAW4wxmbVtbOnzVU9+8MrP2Alaq9W3llboqVit1juAP9jrZmN98AECsaoBtgMrgF6tENPpWP9OrQPW2I+pwC3ALXaZO4CNWD0NlgGntdL56mUfc619/Opz5hybAM/Z53Q9kNJKsQVjJe4Ip3Wtfs6wvmj2A+VYV1A3YbXrLAS22c9RdtkU4GWn995of9a2Aze0Qlzbsepsqz9n1T3MugOf1/c7b+G4XrU/O+uwElhszbjs1yf9/bZkXPb6edWfKaeyrXm+6soPHv+MNfTQIRCUUqqd0ztjlVKqndNEr5RS7ZwmeqWUauc00SulVDuniV4ppdo5TfRKKdXOaaJXSql27v8BgFqBJa86DFAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = log.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction\n",
    "\n",
    "For evaluation we first need to obtain survival estimates for the test set.\n",
    "This can be done with `model.predict_surv` which returns an array of survival estimates, or with `model.predict_surv_df` which returns the survival estimates as a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "surv = model.predict_surv_df(x_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the survival estimates for the first 5 individuals.\n",
    "Note that the time scale is correct because we have set `model.duration_index` to be the grid points.\n",
    "We have, however, only defined the survival estimates at the 10 times in our discretization grid, so, the survival estimates is a step function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAczUlEQVR4nO3df3RV5Z3v8feXEEiXgBRIBAkWGEAhSvkRQKbqxFudQurAXBXFHxXRDnautnVpu6r1jjrOcrBWb0tvvVOpCtbeEbW9jlylWFvLTJdXRVAUBKMMohx+SAAREAIh+d4/ziYcshNyQs7O3if5vNY6K2fvs8/Zn2xCvnn2fvbzmLsjIiKSqUvcAUREJHlUHEREJETFQUREQlQcREQkRMVBRERCusYd4ET069fPBw8eHHcMEZG8sXLlyh3uXpzt9nlZHAYPHsyKFSvijiEikjfM7KPWbK/TSiIiEqLiICIiISoOIiISkpfXHERE4lBbW0sqlaKmpibuKM0qKiqitLSUwsLCNn2OioOISJZSqRQ9e/Zk8ODBmFnccULcnZ07d5JKpRgyZEibPivS00pm9piZbTezNc28bmb2MzNbb2bvmNm4KPOIiLRFTU0Nffv2TWRhADAz+vbtm5OWTdTXHBYCU47z+lRgePCYA/xLxHlERNokqYXhiFzli/S0krv/h5kNPs4m04FfeXrc8NfMrLeZDXD3rcf73B0bd/DLWb/MYdK2qy/5jBt+/L24Y4iI5ETcvZUGApsyllPBuhAzm2NmK8xsBe50oT4xj7rCARRs79UuB0xEZOnSpZx++ukMGzaM++67L5J9xH1Buqn2T5OzD7n7fGA+QPngk/36isVR5mqVR5f9DU1/KyIiuVVXV8eNN97ISy+9RGlpKRMmTGDatGmMGjUqp/uJuzikgEEZy6XAlhbf1W84zH4hqkytt+zhuBOISCexfPlyhg0bxtChQwGYOXMmzz33XIcrDouBm8xsETAJ+Kyl6w0iIknwj//3XdZu2ZPTzxx1ai/u+puy426zefNmBg06+jd1aWkpr7/+ek5zQMTFwcyeBCqAfmaWAu4CCgHc/RfAEqASWA/sB2ZHmUdEJN+l++8cK4oeVFH3VrqihdcduDHKDCIiUWjpL/yolJaWsmnT0X48qVSKU089Nef7ibu3koiItMKECRP44IMP+PDDDzl06BCLFi1i2rRpOd9P3NccRESkFbp27crPf/5zvva1r1FXV8d1111HWVnuWzEqDiIieaayspLKyspI96HTSiIiEqLiICIiITqtlCP1OLOXJqsnbuXQSmaMmBF3DBHJQyoOOVDg0K0WZv7Pd+OO0mD/4QN8NOkTuEfFQURaT8UhB7oWdQGvp/5g7vsan6jTqj+gcHl13DFEJE+pOORAYY+u1H6hnoXn/CDuKA1mPH0jBXX1cccQkTyl4pAD3Qq60K2gC0/dMDnuKA2e/Y1GiRXpiK677jqef/55SkpKWLOmyUk2c0K9lURE8si1117L0qVLI9+PioOISB4577zz6NOnT+T70WklEZET8bvbYNvq3H5m/7NgajQzu7WWWg4iIhKiloOIyIlIyF/4UVHLQUREQlQcRETyyBVXXMHkyZOpqqqitLSURx99NJL96LSSiEgeefLJJ9tlP2o5iIhIiIqDiIiEqDiIiEiIioOIiISoOIiISIiKg4iIhKg4iIjkkU2bNnH++eczcuRIysrKmDdvXiT70X0OIiJ5pGvXrjz44IOMGzeOvXv3Mn78eC688EJGjRqV0/2o5SAikkcGDBjAuHHjAOjZsycjR45k8+bNOd+PWg4iIifgR8t/xHu73svpZ57R5wx+MDH76YY3btzIW2+9xaRJk3KaA1QcOrR6DjJ76ey4YxyjcmglM0bMiDuGSN7bt28fl1xyCT/96U/p1atXzj9fxaGDKvCekLBppKt2VQGoOEiH0Jq/8HOttraWSy65hKuuuoqLL744kn2oOHRQhfTmUE0P9n80J+4oDeq6PcD2PQfjjiGS19yd66+/npEjR3LLLbdEth9dkO6g+vbozkndk1X79x88zI7PVRxE2uKVV17hiSee4OWXX2bMmDGMGTOGJUuW5Hw/kf72MLMpwDygAHjE3e9r9PppwONA72Cb29w9999lJ3RKz+6c0rM7T90wOe4oDSYtSFaxEslH55xzDu4e+X4iazmYWQHwEDAVGAVcYWaNO+L+d+Bpdx8LzAT+V1R5REQke1GeVpoIrHf3De5+CFgETG+0jQNHLrOfDGyJMI+IiGQpyuIwENiUsZwK1mW6G7jazFLAEuDbzX2Ymc0xsxVmtqK6ujrXWUVEJEOUJ4Gb6kjZ+ETZFcBCd3/QzCYDT5jZme5eH3qj+3xgPkB5eXn0J9xaacf+Yp598M24YzSoKZpCSer/wTeuiTtKg9u2beTVUSfHHUNEshBlcUgBgzKWSwmfNroemALg7q+aWRHQD9geYa6cG9FvHewA6B93lAZ7u5dA6V8yrGZp3FEanLa9Ju4IIpKlKIvDG8BwMxsCbCZ9wfnKRtt8DHwVWGhmI4EiIO/OGZWVrKasZDXMnhV3lAbpVkwfvnRr40MenzcvHB93BBHJUmTXHNz9MHAT8CKwjnSvpHfN7B4zmxZsdivwd2b2NvAkcK23Rx8tEZE8VVNTw8SJE/nyl79MWVkZd911VyT7ibTjeXDPwpJG6+7MeL4W+EqUGUREOpLu3bvz8ssv06NHD2praznnnHOYOnUqZ599dk73ozukRUTyiJnRo0cPID3GUm1tLWa5H0hNt6yKiJyAbf/8zxxcl9shu7uPPIP+P/xhi9vV1dUxfvx41q9fz4033hjJkN1qOYiI5JmCggJWrVpFKpVi+fLlrFmzJuf7UMtBROQEZPMXftR69+5NRUUFS5cu5cwzz8zpZ6vlICKSR6qrq9m9ezcABw4c4A9/+ANnnHFGzvejloOISB7ZunUrs2bNoq6ujvr6ei677DIuuuiinO9HxUFEJI+MHj2at956K/L96LSSiIiEqDiIiEiIioOIiISoOIiISIiKg4iIhKg4iIhIiIqDiEieqaurY+zYsZHc33CEioOISJ6ZN28eI0eOjHQfKg4iInkklUrxwgsv8M1vfjPS/egOaRGRE/Dnp99nx6Z9Of3MfoN6cO5lI467zc0338z999/P3r17c7rvxlQccmXbaljw9bhTHLXtMjipOO4UIpJDzz//PCUlJYwfP55ly5ZFui8Vh1w469K4E4Qd+jzuBCIdWkt/4UfhlVdeYfHixSxZsoSamhr27NnD1Vdfza9//euc70vFIRfKZ6cfSfL9x+NOICI5NnfuXObOnQvAsmXLeOCBByIpDKAL0iIi0gS1HERE8lBFRQUVFRWRfb5aDiIiEqKWg7Sreg4ye2myrs9UDq1kxogZcccQSRS1HKTdFHhPutA97hjHqNpVxZINS+KOIXnE3eOOcFy5yqeWg7SbQnpT6L1ZMGVB3FEaJK0VI8lWVFTEzp076du3L2YWd5wQd2fnzp0UFRW1+bNUHEREslRaWkoqlaK6ujruKM0qKiqitLS0zZ+j4iAikqXCwkKGDBkSd4x2oWsOIiISouIgIiIhKg4iIhIS6TUHM5sCzAMKgEfc/b4mtrkMuBtw4G13vzLKTJ3Jjv3FPPvgm3HHaLBr4MX02FsVdwwRyUJkxcHMCoCHgAuBFPCGmS1297UZ2wwHbge+4u6fmllJVHk6mxH91sEOgP5xR2lwqFsx+3rGnUJEshFly2EisN7dNwCY2SJgOrA2Y5u/Ax5y908B3H17hHk6lbKS1ZSVrIbZs+KO0uDRWb+KO4KIZCnKaw4DgU0Zy6lgXaYRwAgze8XMXgtOQ4mISMyibDk0dftg4/u6uwLDgQqgFPizmZ3p7rtDH2Y2B5gDcNppp+U2qYiIHCPKlkMKGJSxXApsaWKb59y91t0/BKpIF4sQd5/v7uXuXl5crOkvRUSiFGXL4Q1guJkNATYDM4HGPZH+DbgCWGhm/UifZtoQYSaRkKpdVYkcY0mjxUqcsioOZlYOnAucChwA1gB/cPddzb3H3Q+b2U3Ai6S7sj7m7u+a2T3ACndfHLz212a2FqgDvu/uO9v0HYm0QuXQyrgjNKlqV7rLr4qDxOW4xcHMrgW+A3wIrCR92qcIOAf4gZmtAf7B3T9u6v3uvgRY0mjdnRnPHbgleIi0uxkjZiTyF3ASWzLSubTUcjiJ9D0IB5p60czGkL5G0GRxEBGR/HTc4uDuDzX3mpl1c/dVuY8kHVldvXP5w6/GHeMY08cM5MpJ6gEnkimr3kpmtszMBmcsTyR9wVkka4UFXSjokqwJUtZu3cNzqzbHHUMkcbLtrTQXWGpmPyN9I9tUQCdFpVUKC7pQWNCFp26YHHeUBklrxYgkRVbFwd1fNLNvAS+RHrFnrLtvizSZdEj1+/fz0TeuiTtGg2u37mH16ZOA5BQskSTItivrPwCXAecBo4FlZnaru78QZTjpWLr27cthgJq4kxzVv1p9KUSaku1ppX7AxKDX0qtmthR4BFBxkKx1LS6ma3ExX7o1OaOyr73gb+OOIJJI2Z5W+m6j5Y9ID8UtIiIdkGaCExGREBUHEREJUXEQEZGQlsZWOi/Lz9nY3PhKIiKSf1q6IJ3tjW7PovGVREQ6jJbGVtJd0CIxSeI8E5pjovOIcrIfETlBSZxnQnNMdC4qDiIJlMR5JpLWipFoqbeSiIiEZDtk9xPZrBMRkY4h29NKZZkLZlYAjM99HMmpbathwdfjTnHUtsvgpOK4U4hIFlq6z+F24IfAF8xsz5HVwCFgfsTZpC3OujTuBGGHPo87QZP6V3+cqGHEj+h10UV88fLL4o4hnVRLXVnnAnPNbK67395OmSQXymenH0ny/cfjThCSnssB+sSco7Ga994DUHGQ2LTUchjs7hubKwxmZsBAd09Fkk4kYitHV7BydEWiZqcDEtmSkc6lpWsOPzazLsBzwEqgGigChgHnA18F7gJUHEREOpCWTivNMLNRwFXAdcAA4ACwjvREP/e6e4Lm9RIRkVxosbeSu68F7miHLCIikhDHvc/BzCaYWf+M5WvM7Dkz+5mZJe0anoiI5EhLLYeHgQugYfju+4BvA2NId2VNYH9JSbId+4t59sE3447RYPTWg2z/YkHcMUQSp6XiUODuu4LnlwPz3f23wG/NbFW00aSjGdFvHewA6N/Spu2mx4H6uCOIJFKLxcHMurr7YdI9k+a04r0ixygrWU1ZyWqYPSvuKA3u/d6yuCOIJFJLv+CfBP7dzHaQ7qX0ZwAzGwZ8FnE2ERGJSUtdWe81sz+S7sL6e3f34KUupK89iLROwsZ7+lLtxezp0jvuGCKJk01X1teaWPd+NHGkQ0vgeE9FXgP1u+OOIZI4kV43MLMpwDygAHjE3e9rZrtLgWeACe6+IspMEqMEjvdUs+rRuCOIJFJkk/0Ew3o/BEwFRgFXBHdbN96uJ/Ad4PWosoiISOtEORPcRGC9u29w90PAImB6E9v9E3A/oGE4REQSIsriMBDYlLGcCtY1MLOxwCB3f76lDzOzOWa2wsxWVFdX5zapiIgcI8riYE2s84YX06O9/gS4NZsPc/f57l7u7uXFxZpNTEQkSlEWhxQwKGO5FNiSsdwTOBNYZmYbgbOBxWZWHmEmERHJQpTF4Q1guJkNMbNuwExg8ZEX3f0zd+/n7oPdfTDwGjBNvZVEROIXWXEIhty4CXiR9PwPT7v7u2Z2j5lNi2q/IiLSdpHe5+DuS4Aljdbd2cy2FVFmERGR7EV5WklERPKUioOIiIRo2G3p9Iq8JlGDAQKwbQucpC7bEh8VB+nU9nTpncyB9w59HncC6eRUHKRT+7SgL58W9IXZL8Qd5Vj/OjZdIJLUorFP1JrpRFQcRJIoib+E1ZrpVFQcRJKoZ39qNu/mo5f7xp2kwcxtG1lXVht3DGknKg7S6e0/eJjLH3417hjHGN9rJBUDD1IUd5AMJdvrgcNxx5B2ouIgnVq/Ht3YEXeIJjxRPI6Voyt46obJcUdp8N6FZXFHkHak4iCdWknPIkp6FnHHDePijnKMpLVkpPPRTXAiIhKiloOIZM/rk9W9FuCsSxM3N3lHoJaDiGSnoBAsYb8ytq2G1b+JO0WHpJaDiGSnoFv6kaQbBpPWiulAEvZngIiIJIGKg4iIhOi0kkhCrd26J1FdWmceqqOwQH9PdhYqDiIJNH3MwLgjhNTVO/UcYPbSBPUMsk+o9JOYEXeODkjFQSSBrpx0GldOOi3uGMdY/Ewv6tgbd4xjVHEIDBWHCKg4iEhWCulNofdmwZQFcUdpMHthedwROiydQBQRkRAVBxERCVFxEBGREBUHEREJUXEQEZEQFQcREQlRcRARkRAVBxERCVFxEBGREN0hLZ3ejtQ+nn3wzbhjhIyYeApl5yZvjCXpHFQcpFMbMfGUuCM0aUdqH4CKg8RGxUE6tbJzBybyF3ASWzJJVcWhZI0UC1QOrWTGiPweDjDSaw5mNsXMqsxsvZnd1sTrt5jZWjN7x8z+aGZfijKPiHQslX4Sp9Mt7hjHqNpVxZINS+KO0WaRtRzMrAB4CLgQSAFvmNlid1+bsdlbQLm77zezvwfuBy6PKpOIdCwz6MEM7wFJGik2Ya2YExVly2EisN7dN7j7IWARMD1zA3f/k7vvDxZfA0ojzCMiIlmKsjgMBDZlLKeCdc25Hvhdcy+a2RwzW2FmK6qrq3MUUUREmhLlBWlrYp03uaHZ1UA58FfNfZi7zwfmA5SXlzf5OSISrf7VH/PRN66JO8ZR27bQa1QPvhh3jg4oyuKQAgZlLJcCWxpvZGYXAHcAf+XuByPMIyJtsPr0SQD0iTlHpprth4B9Kg4RiLI4vAEMN7MhwGZgJnBl5gZmNhZ4GJji7tsjzCIibbRydAVPFI9j1IBecUdpcNfHczhUVx93jA4psuLg7ofN7CbgRaAAeMzd3zWze4AV7r4Y+DHQA3jGzAA+dvdpUWUSkRM3fUzy7gepr3dqUXGIQqQ3wbn7EmBJo3V3Zjy/IMr9i0juXDnpNK6cdFrcMY6x9l+NIq+BBV+PO8pR9gmcVBx3ijbTHdIikrc+69Ib6nfHHeNYhz6PO0FOqDiISN7aXdCX3QV9GTX73+KOctTC8rgT5ISG7BYRkRAVBxERCVFxEBGREBUHEREJUXEQEZEQ9VYSkbyWtPGeZm47wLqy/P/Vmv/fgYh0Wkkc76lkez1wOO4YbabiIJJQO1L7Ejdd6IiJpyRqWtWVoytYObqCp26YHHeUBu9dWBZ3hJxQcRBJoBETT4k7QsiO1D6ARBUHiY6Kg0gClZ07MHG/hJPWipFoqbeSiIiEqDiIiEiITiuJiOTYfpzZS2fHHaNNVBxERHKoLwXgtbBtddxRGlRxqNXvUXEQkby2duseLn/41bhjNLj58MkMtT0s8OT0OJtdu4HWHiEVBxHJW0mcunTz4ZPZ3b0vw/N8jgkVBxHJW0mcuvR3z3SMX6sd47sQEUmQJI73tLCV71FxEBHJoeSO99Q6Kg4ikjWN99SyjjLek4qDiGRF4z11LioOIpIVjfeUvaR1r53l3ur3qDiIiORQErvX1nnrR0pScRARyaEkdq997je9Wv0eDbwnItLB7eo3vtXvUctBRPKaelC17OMrvwWLft6q96g4iEjeSmIPqi0f7GbLB7t5f/kncUdpMPoE3qPiICJ5K4k9qN798+ZEFYYTpeIgIpJDSSxYAHyvdZvrgrSIiIREWhzMbIqZVZnZejO7rYnXu5vZU8Hrr5vZ4CjziIhIdiIrDmZWADwETAVGAVeY2ahGm10PfOruw4CfAD+KKo+IiGQvypbDRGC9u29w90PAImB6o22mA48Hz38DfNXMLMJMIiKShSiLw0BgU8ZyKljX5Dbufhj4DOjb1IeZ2RwzW2FmK6qrqyOIKyIiR0RZHJpqATQe/SmbbdIr3ee7e7m7lxcXF7c5nIiINC/K4pACBmUslwJbmtvGzLoCJwO7IswkIiJZiLI4vAEMN7MhZtYNmAksbrTNYmBW8PxS4GX3ExhbVkREcsqi/F1sZpXAT4EC4DF3v9fM7gFWuPtiMysCngDGkm4xzHT3DVl87l6gKrLg0eoH7Ig7xAlS9ngoezw6WvYvuXvW5+QjLQ5RMbMV7l4ed44ToezxUPZ4KHs8cpFdd0iLiEiIioOIiITka3GYH3eANlD2eCh7PJQ9Hm3OnpfXHEREJFr52nIQEZEIqTiIiEhIXhWHloYATxoz22hmq81slZmtCNb1MbOXzOyD4OsX4855hJk9ZmbbzWxNxrom81raz4J/i3fMbFx8yZvNfreZbQ6O/6rgvpsjr90eZK8ys6/FkxrMbJCZ/cnM1pnZu2b23WB94o/7cbIn/rgHWYrMbLmZvR3k/8dg/ZBgCoEPgikFugXrEzPFwHGyLzSzDzOO/Zhgfet/btw9Lx6kb6T7T2Ao0A14GxgVd64WMm8E+jVadz9wW/D8NuBHcefMyHYeMA5Y01JeoBL4Henxsc4GXk9g9ruB7zWx7ajg56c7MCT4uSqIKfcAYFzwvCfwfpAv8cf9ONkTf9yDPAb0CJ4XAq8Hx/Rp0jfkAvwC+Pvg+X8DfhE8nwk8lcDsC4FLm9i+1T83+dRyyGYI8HyQOUz548DfxpjlGO7+H4THtmou73TgV572GtDbzAa0T9KwZrI3ZzqwyN0PuvuHwHrSP1/tzt23uvubwfO9wDrSoxUn/rgfJ3tzEnPcAYJjuC9YLAweDvwX0lMIQPjYJ2KKgeNkb06rf27yqThkMwR40jjwezNbaWZzgnWnuPtWSP/nAkpiS5ed5vLmy7/HTUEz+rGMU3iJzB6cphhL+q/AvDrujbJDnhx3Mysws1XAduAl0q2Z3Z6eQgCOzZj1FAPtoXF2dz9y7O8Njv1PzKx7sK7Vxz6fikPWw3snyFfcfRzp2fBuNLPz4g6UQ/nw7/EvwF8AY4CtwIPB+sRlN7MewG+Bm919z/E2bWJd0rLnzXF39zp3H0N61OiJwMimNgu+Jip/4+xmdiZwO3AGMAHoA/wg2LzV2fOpOGQzBHiiuPuW4Ot24FnSP3yfHGnOBV+3x5cwK83lTfy/h7t/EvwHqgd+ydFTGInKbmaFpH+5/m93/z/B6rw47k1lz5fjnsnddwPLSJ+P723pKQTg2IyJnGIgI/uU4FSfu/tBYAFtOPb5VByyGQI8MczsJDPreeQ58NfAGo4dpnwW8Fw8CbPWXN7FwDVBL4izgc+OnAZJikbnVP8r6eMP6ewzg94nQ4DhwPL2zgfpXiTAo8A6d/8fGS8l/rg3lz0fjjuAmRWbWe/g+ReAC0hfN/kT6SkEIHzsEzHFQDPZ38v4g8JIXyvJPPat+7mJ62r7iTxIX3F/n/R5wTviztNC1qGke2a8Dbx7JC/pc5R/BD4IvvaJO2tG5idJnwaoJf2XxvXN5SXdTH0o+LdYDZQnMPsTQbZ3gv8cAzK2vyPIXgVMjTH3OaSb9+8Aq4JHZT4c9+NkT/xxD7KMBt4Kcq4B7gzWDyVdtNYDzwDdg/VFwfL64PWhCcz+cnDs1wC/5miPplb/3Gj4DBERCcmn00oiItJOVBxERCRExUFEREJUHEREJETFQUREQrq2vIlI52VmR7qUAvQH6oDqYHm/u/9lLMFEIqaurCJZMrO7gX3u/kDcWUSiptNKIifIzPYFXyvM7N/N7Gkze9/M7jOzq4Lx9leb2V8E2xWb2W/N7I3g8ZV4vwOR5qk4iOTGl4HvAmcB3wBGuPtE4BHg28E284CfuPsE4JLgNZFE0jUHkdx4w4OxaszsP4HfB+tXA+cHzy8ARmVMAdDLzHp6ei4EkURRcRDJjYMZz+szlus5+v+sCzDZ3Q+0ZzCRE6HTSiLt5/fATUcWjszvK5JEKg4i7ec7QHkwS9da4FtxBxJpjrqyiohIiFoOIiISouIgIiIhKg4iIhKi4iAiIiEqDiIiEqLiICIiISoOIiIS8v8BsYezOPAkv9UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv.iloc[:, :5].plot(drawstyle='steps-post')\n",
    "plt.ylabel('S(t | x)')\n",
    "_ = plt.xlabel('Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is, therefore, often beneficial to interpolate the survival estimates, see [this paper](https://arxiv.org/abs/1910.06724) for a discussion.\n",
    "Linear interpolation (constant density interpolation) can be performed with the `interpolate` method. We also need to choose how many points we want to replace each grid point with. Her we will use 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "surv = model.interpolate(10).predict_surv_df(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5hU1Znv8e8rIG24yFVBGgSjclOC0kI0ypATVCSeUaOAmkRFZ0jmSCbXOerkTKLOJF5GY8wTnyTGiIkZbyTjkUPaS4zDTDQGaSLKRVsZAWmuDSgXsRGadf7YVdW7d9e1e++qXVW/z/P0Q9euTdWyBN5+17vWu8w5h4iIiN8RpR6AiIjEj4KDiIh0oOAgIiIdKDiIiEgHCg4iItJB91IPoDMGDRrkRo4cWephiIiUjeXLl+9wzg3O9/6yDA4jR46koaGh1MMQESkbZrahkPs1rSQiIh0oOIiISAcKDiIi0kFZ1hxERErh4MGDNDU10dLSUuqhZFRTU0NtbS09evTo0usoOIiI5KmpqYk+ffowcuRIzKzUw+nAOcfOnTtpampi1KhRXXqtSKeVzOxBM9tuZqsyPG9m9iMzW2tmr5vZ6VGOR0SkK1paWhg4cGAsAwOAmTFw4MBQMpuoaw4PATOyPH8BcFLiax7wk4jHIyLSJXENDElhjS/SaSXn3H+Z2cgst1wE/Mp5fcP/bGb9zGyoc25LttfdsX4HP7/65wAcPmY3X/rXb4U1ZBERofSrlYYBG32PmxLXOjCzeWbWYGYNOMcRHKa1x1C6be9blIGKiMTFM888w+jRoznxxBO5/fbbI3mPUgeHdPlP2tOHnHP3O+fqnHN1g/q8z3XTFtHjYBNHfnSYJ8+dlPp64LvzIx6yiEjptLa2cv311/P000+zZs0aHn30UdasWRP6+5R6tVITMNz3uBbYnPN3DToJ5v4Ot/TnHORw6vKYjfvh8T/w5J8mpa7tPOtM/uaWH4c3YhGREnrllVc48cQTOeGEEwC4/PLLeeqppxg3blyo71Pq4LAImG9mjwFTgN256g1+A/scCX3gkl8uB+CB785n4J9eTj0fDBYKFCISllv+32rWbN4T6muOO64v3/2f47Pes2nTJoYPb/uZura2lqVLl4Y6Dog4OJjZo8A0YJCZNQHfBXoAOOd+CtQDM4G1wH5gbqHvsWP/YJ68+y8ADOx7LSd/5ybGn+OVLfzBYvj2/bS+/AfmPtP2FjNPmMmsk2d19j9PRKTovPU77UWxgirq1UpX5HjeAdd39vVPHvQG7AAYAsCOpn0AqeDgzxIWTZ/GyK3N/PXdrwPQavt5cdwr1M+oT92jYCEi+cr1E35Uamtr2bixbR1PU1MTxx13XOjvU+pppS4Zf8xKxh9+BAacCsCTW2fD3vTtyntOv5zm55+hJhF0R25sZPy7h3lztYKFiJSPM844g7fffpt169YxbNgwHnvsMR555JHQ36esgwOnXtb+8UcfZLz1/Bu/DDd+OfX42dt/yuECgoUChYjEQffu3fnxj3/M+eefT2trK9deey3jx4efxVi6+au4q6urc+kO+3nyH37Jjv2DGXTCkNS1kycfm5pmyiYZLJIGb3+Hjcf2YNE3J9C4q5HRA0azYMaCcP4DRKQsvfHGG4wdO7bUw8gp3TjNbLlzri7f1yjvzCEgVw0im2Bm8fT0ixm+dQOz73qPFjvAn055jbmomC0i1aGigsP4Y1Yy/piVMPdqgNQqps44YvoMmhOZxJiN+xmzEd5cpfqEiFSHigoOYfJnEqpPiEi1qbzgsHUlLPhs4vvZ7DgwLJVB5Ft/CMq3mP3u26vZf+hDNkzZBrcqOIhI+aqs4BBYvXRyj+eB6cCAguoPuQSDxZ3X38YJr79ErwPdGdH8Nu7P72qznYiUtcoKDnVzva+E8Qs+y3iegLlXd6n+kEvtVZ9n4YppAMx67Aa6ubYltY27Gr3rCg4iUkYqKziUyJVTRnDllBEAPL1wAIO37mX2Xe8B0GIHWHnapuxHHomI5Onaa69l8eLFHHPMMaxalfaQzVBUVXDY0bSvXQbR2RpENv5VTgDDt7Vw+NXNmmYSkVBcc801zJ8/n6uuuirS96ma4HDy5GPbPQ6zBuEXrEcsmj4Nd3gPa7Z43RtbbCPb9xxQcBCRTpk6dSrr16+P/H2qJjiMP2dYu0AQZQ3Cb0jvY+mxbg83/1srAB+4Q/zp1DXM7atMQqSsPX2jtzoyTENOhQuiOdmtUKU+CS56yaWtCz4LDcVvf3HSFZcyYMIpjBval3FD+zJqxyHOWd0t9Xzjrkbq36nP8goiIsVX2ZmDf2lrMsLXFXxkRJf0nzOb/nNmpx6vmX4xrS2H2L9hHgCtR97FyuY3UjUJZREiZSImP+FHpbKDg39pa3JjnE8xCtRBA3v3ZEjzu1yz8A4APmjdxSsTPkbzBVr2KiLxUfnTShmcPPlYBtX2Tj3e0bSPt17ZFvn7BqeZPr5rH5P/8jH2b5hH64GhbN9zIPIxiEj5uuKKKzjzzDNpbGyktraWX/ziF5G8T2VnDlmUqkAdnGbadvEceu3zAsL+A4doso1a9ioiGT366KNFeZ+qDQ5xcWyfnhy9aR13vvgTXtuymxdGH8XLdVr2KiKlVV3Bwd+UD7yCta9AXYoaRN8LL0x9P2rXPs5vrGXLhG8BsMbdThNrlUmISNFVT3AIHikaWL1UrE1yQf5ppg1fvIpxwONfOhOA8x44i93uldS9KliLSLFUT3AINOULrl4qVQ0iqOXNN9nwRW9b/Ne37OHZIZ9gyzleY6bgsldQJiEi0aja1Upx1PfCC6kZMyb1uHZXE+dvfS31eP+uCXQ71BbAtIFORKJSPZlDJ/hrEMWoPwRXMgWnmeb8DNZs+CT7P+oLeJmElr6KSBQUHDLw1yCKVX9Ixz/N9M29B1gybCLLh04DOi591RSTSOXbuHEjV111FVu3buWII45g3rx5fPWrXw39fRQcMvDXIEpVf/CvZAI4etM65vTpyf/+0k2AV7De0vIn1mzZo2WvIlWie/fu3H333Zx++uns3buXSZMmce655zJu3Lhw3yfUVys3OZa2llq6aSa/a069gqdWTIWPtOxVpFoMHTqUoUOHAtCnTx/Gjh3Lpk2bFBxCk2Npa1z5p5k+BVxw4YX0nzNby15FiuyOV+7gzV1vhvqaYwaM4YbJN+R9//r163n11VeZMmVKqOOAag4OOZa2BpVig1xQcJqp5U3vD2b/ObPp3zqVLU0T2xWrtexVpHLt27ePSy+9lB/+8If07ds39Nev3uBQgFJtkAvKNs100cT2Y9m/awIfG9D2WJmESLgK+Qk/bAcPHuTSSy/l85//PJ/73OcieQ8FhzzEZYNcOslpJv8UE2jZq0ilcs5x3XXXMXbsWL7xjW9E9j4KDmXMP83kn2KCNJnEgUPsKN7QRCQiL730Eg8//DCnnnoqEydOBOD73/8+M2fODPV9Ig0OZjYDuBfoBjzgnLs98PwI4JdAv8Q9NzrnSrflt4DVS3GoQQT7MvldOWUEV04ZkXo8ZUF3WrQnQqTsnX322TjnIn+fyIKDmXUD7gPOBZqAZWa2yDm3xnfb/wGecM79xMzGAfXAyKjGlFUBq5fiUoMI8q9kAi+zSAaPo1sne+EX1R9EJLcoM4fJwFrn3DsAZvYYcBHgDw4OSJbZjwY2Rzie7ApYvRTHGkS2lUwA/Vun0r91KgtmnNluBZOISDpRBodhwEbf4yYguBj3ZuA5M/sK0AuYnunFzGweMA9gxIgRmW6rWrk2zAGs2bKHOT97mfVH7qG1+yYtcxWRjKIMDpbmWnCi7ArgIefc3WZ2JvCwmZ3inDvc4Tc6dz9wP0BdXV30E24FKnaTvnxk6sukZa4ikkuUwaEJGO57XEvHaaPrgBkAzrmXzawGGARsj3BcoYtLkz6/bH2Z5vwM+Og8Fszwur1qmklEgqIMDsuAk8xsFLAJuBy4MnDPu8BngIfMbCxQAzRHOKbC+FcvZVm5FIcmfUH5TDOJiGQSWXBwzh0ys/nAs3jrZB50zq02s1uBBufcIuCbwM/N7Ot4U07XuGKs0cqHf/VSmfRd6orGXY2qQYiUgZaWFqZOncqBAwc4dOgQl112Gbfcckvo7xPpPofEnoX6wLXv+L5fg9c/Ln78q5dy9F0KisMeiELMPKH95hnVIETiq2fPnrzwwgv07t2bgwcPcvbZZ3PBBRfwyU9+MtT30Q7pkMV1DwS0FaivSZxNPSf1TC0XTbwltWlONQiR+DIzevfuDXg9lg4ePIhZuvU/XaPgELI47oGA9gXq2l1NnA885K0FYM2WPQDtdlRrmkkku63f/z4H3gi3ZXfPsWMY8o//mPO+1tZWJk2axNq1a7n++uvVsrtcxWGaKdhqo/3Z1C+3u1fTTCLx1q1bN1asWMH777/PJZdcwqpVqzjllFNCfQ8Fh4jFeZopk1knz2oXCDTNJNJRPj/hR61fv35MmzaNZ555RsGhZDp5pGhcp5lEpDw1NzfTo0cP+vXrx4cffsjzzz/PDTeEf7aEgkM+yvRI0Wz8u6c7Fqi9lt+ZahCqP4iUzpYtW7j66qtpbW3l8OHDzJ49mwsDm17DoOCQjwKPFM2l1DWI4O7pXAVqfw1C9QeR0powYQKvvvpq5O+j4FBkcahBpNs9na1A7a9BqP4gUh0UHIqsEmoQWuYqUvkUHKQgWuYqUh0UHDqrk6uX0olDu+98C9Ra5ipSHRQcOiPE1UtxaPddaIFaRCqfgkNnhLh6KQ7tvgstUAepBiFSeRQcYqbUy1wLpRqESPG1trZSV1fHsGHDWLx4cSTvoeAQI3FY5loo1SBEiu/ee+9l7Nix7NmzJ7L3UHCIkTgtc81WoA7ung7SbmqR6DQ1NfG73/2Ob3/72/zgBz+I7H0UHKSDbAXqXMVp7aaWavHHJ95ix8Z9ob7moOG9OWf2yVnv+drXvsadd97J3r17Q33vIAWHsOR53nShSlGDyFagzlWc1m5qkegsXryYY445hkmTJrFkyZJI30vBIQwRnTddjjUIkWqR6yf8KLz00kssWrSI+vp6Wlpa2LNnD1/4whf49a9/Hfp7KTiEoQvnTWcTpxpEZ2mZq0h4brvtNm677TYAlixZwl133RVJYAAFh7JTqqWumc+fzlyg1jJXkfKl4BCFEFtr+JVqmqnQ86eTtMxVJDrTpk1j2rRpkb2+gkPYIjwYqFTTTIWcPy0ilUHBIWwhHwxUaVSDECkPCg5lLg4dXfMVrEE0bGugYVsD9e/Ut7tHwULizDmHmZV6GBk550J5HQWHYihCDaKYy1yDu6dXjp4CnJnz9wVrEAvfWtguMASDhQKFxE1NTQ07d+5k4MCBsQwQzjl27txJTU1Nl19LwSFqRapBFKv+ENw9PaT5XT44cKhd7SFXe42kbMFCK5skjmpra2lqaqK5ubnUQ8mopqaG2traLr+OgkPUKqwGEdw9ve3iOfTadyD1uCtnPwR3V6s+IXHTo0cPRo0aVephFIWCQwUpxR6IY/v05Ng+PUNfvaT6hEhpKTiUQgQ1iEprtZGrPqFpJ5FoKTgUW0Q1iFK22uhsgboQ6TbUqTW4SHQiDQ5mNgO4F+gGPOCcuz3NPbOBmwEHvOacuzLKMZVcEWsQxZhmylWgzrc4XSj/tJOmnETCF1lwMLNuwH3AuUATsMzMFjnn1vjuOQm4CfiUc+49MzsmqvFUm2JNM2UrUHelOJ2LP5PItSQWFCxEChVl5jAZWOucewfAzB4DLgLW+O75W+A+59x7AM657RGOJ74iOAuiVNNM/gJ1sVprqD4hEr4og8MwYKPvcRMwJXDPyQBm9hLe1NPNzrlnIhxT/ER0FkQ1y1WfAGUSIrlEGRzSbR8M7uvuDpwETANqgT+a2SnOufc7vJjZPGAewIgR4U9TlExEZ0GkU6xWG/723lEUpwuVa1msAoVIR1EGhyZguO9xLbA5zT1/ds4dBNaZWSNesFgWfDHn3P3A/QB1dXXhNA+pIsVqteEvUHdl93SYtBNbpHBRBodlwElmNgrYBFwOBFci/V/gCuAhMxuEN830ToRjir+I+jAVq9WGv0Ad5u7pMGkntkhueQUHM6sDzgGOAz4EVgHPO+d2Zfo9zrlDZjYfeBavnvCgc261md0KNDjnFiWeO8/M1gCtwD8453Z26b+onEXYhymoGMtco9o9HSbtxBZJL2twMLNrgL8H1gHLgUagBjgbuMHMVgH/5Jx7N93vd87VA/WBa9/xfe+AbyS+pEh7ICptN3VXaKWTSHq5ModeeHsQPkz3pJlNxKsRpA0OEk/FXOZajN3TYdJKJxFP1uDgnLsv03NmdqRzbkX4Q5JKEWZ771LRSiepVvnWHJYA1zjn1iceTwZ+DnwispGJJ6ICdVAUNYgo23sXi1Y6SbXKd7XSbcAzZvYjvM1tFwDaqRW1IhWoi1WDKIcCdS5a6STVIq/g4Jx71sy+DPwe2AGc5pzbGunIpGgFatUgOic45aRMQipJvtNK/wTMBqYCE4AlZvZN59zvohyclE4U00zpahDlTMVrqWT5TisNAiYnVi29bGbPAA8ACg7FVoQaRFTTTMEaxJrpFzN42waenn5x6toR02dw/o1f7tL7lIqK11JJ8p1W+mrg8Qa8VtxSTEWqQRRrmumI6TNofr6tz+LIpkZ4qJENq/8EeJmGP5jEnYrXUknM24dWXurq6lxDQ0Oph1F6yQxibrQJXDI4XPLN0yN9nzuvv41TG5cybmhfWt58k5oxYzj+4V9F+p7FkpxyGj1gdOqaMgkpJjNb7pyry/d+HRMqeSlGR9flE6axfMI0Hv/SmWz44lXtitdQfpmEn4rXUm4UHMpdBAcFBRWroyt4ex/m/OxlJvUdy7RhB6hJXN+/bBn7ly1jz+LFqXvLKVioeC3lJldvpal5vs76TP2VJEJFOiioWB1dL5rYFnAeHnx6KosAeO/xJ9oFhpY33wQom+AQpIZ/EndZaw5mtiDP13ky0WW1KFRzSKOI9YcdTfsYVNs7dS2KaabkBrlkcAhKTjvVjBkDlFcWkU66hn+jB4xmwYx8/wqKZBdqzcE5p13Q0k5cOrr690yU+5QT5J52UhYhxaaaQyUpwh6IYu6mzsa/ZyI45VQJwcI/7aTitZSCgkOlKOJBQcWSLE4nZergGtxcVwn1iWAPJ5FiU3CoFEXqw5ROFK02/MVpKKyDazBYVNqyWJFiyLe30sPOuS/muibVJ6oaxJVTRrQLBF3p4Brs6RScdlKgEOko38xhvP+BmXUDJoU/HAlVFdUgssk27VQuU07aEyHFlmufw03APwJHmdme5GXgI+D+iMcmXRGsQWx40fta+Zu25yOqRxRjN3VX+INFOUw5aXe1lEKupay3AbeZ2W3OuZuKNCYJQ7AG0bCgLTBEWKyOcjd1vgXqQgSnnOKYSaRb5ioStVyZw0jn3PpMgcHMDBjmnGuKZHQSHn+wWPDZyKacotpN3ZUCdTYqXoukl6vm8K9mdgTwFLAcaAZqgBOBTwOfAb4LKDiUk1xTTsl7YrQMNswCdTblkEmAahASvVzTSrPMbBzweeBaYCjwIfAG3kE/33POtUQ+SglXtiknCDVYRLHMNUrpMom4UQ1CiiHnaiXn3Brg20UYi5RKocEiz0ARdasNfw0ijPpDJnGbZlINQoohV83hDGCjc25r4vFVwKXABuBm59yu6IcoRZctWBSQVUS5zNVfgwir/pBOuUwziYQtV+bwM2A6pNp33w58BZiIt5T1ssy/VSqGP1gEs4oCVz6FNc3kr0FEVX+A8phmEolCruDQzZcdzAHud879Fvitma2IdmgSS+nadOS58ikuHV27yj/NVOoppiR1cJWw5QwOZtbdOXcIb2XSvAJ+r1SDAjbblcNu6lz800xxmWJSB1eJQq5/4B8F/tPMduCtUvojgJmdCOyOeGxSDkKqT3RVFBvk0gnuro4DdXCVKORayvo9M/sD3hLW51zbsXFH4NUeRNrLVp8IBouts6HX4C6/ZVQb5PIRt5VMImHJZynrn9Nceyua4UhFybVE9qMP2LF/cJcL1MXaIBcU15VM2iAnYYi0bmBmM4B7gW7AA8652zPcdxmwEDjDOafDoStVIFicfMeNsBnY2gzAjv2DYe/WsilQx3ElkzbISVgiCw6Jtt73AefitddYZmaLEpvq/Pf1Af4eWBrVWCSexn9mNON9mcSTfzkfdtK28ilmLTzyUeppJm2Qk7BEmTlMBtY6594BMLPHgIuANYH7/hm4E/hWhGOROApOO22ohw+aE993rZhdrAK1X1ynmUQ6I8rgMAzY6HvcBEzx32BmpwHDnXOLzSxrcDCzeSSW0o4YEX2hUUqgzxB27O7Nk7v+GdjKyUf9F+NJbLIrIFiUqkAdx2kmkc6KMjhYmmsu9aTX7fUe4Jp8Xsw5dz+JA4bq6upcjtulDLU7C2J3b+gzm/FzE2WqAnZml6pAnU4cNsypQC2dEWVwaAKG+x7X4pUfk/oApwBLvGMhGAIsMrO/VlG6OmU9CyLXzuwY1ifisGFOBWrprCiDwzLgJDMbBWwCLgeuTD7pnNsNDEo+NrMlwLcUGCQv/p3ZeUw5laIGEYcNcypQS2dFFhycc4fMbD7wLN5S1gedc6vN7FagwTm3KKr3lsqQtUlfAZvtfnSghaeOPos/4P0UXcxNcn6lXskkUohI9zk45+qB+sC172S4d1qUY5HyUlCTvhyb7Y7d1cA8GpjXzws0q4/czUv7Pw2cGfq4M4nTSibVICQfap4nsdSlJn05gsXIg++EMcSCxOWsatUgJF8KDlL5AsFi/ffP9gJECYvZpcokVIOQfCk4SNV56ahPs/+jVj62ZTfjP1rZsZhdiE4GFu2JkLhTcJCyEdYpcr3P+lvuWuFNr4zb8u9cduTLjO/MgNKtkipEILCoYC1xouAgZSHMU+TaHzEKt/I5Hp/bieJ0cJVUIQKBpe+APdDvYGpzX8v2j2DrSvrv/6V3f4RTXzpFTtJRcJCyEMtT5IKF70IEAkv/iX3pP7Fv6vGGR3z7RSM8NEmnyEkmCg4ipZArsLyQWM30wkDYO4W+x39I/+MTz+WaziogcOgUOclEwUHKlr8G0dn6A3TcPR1UjN3UQe1ab2x6H/qMof/cX3kXsk1nFfFoVqlsCg5Slto16etC/SHYwTVo6bpdLF23i6dWbGr3e0raeiNb1pHraNZ0fMFDG+QkydqOhS4fdXV1rqFBLZjEk8weLvnm6aG/9iNL320XGJau2wXAlFEDgOIEiuSGuZoxY1LX8l7JlKtovuFF79fjz2Yh+6jv/THoMwTwAsXoAaNZMGNBV4YvMWFmy51zdfner8xBJItg+29/sChWVtGlDXO5ahu+4DFrwwpmARx/NgBz7SPYu7UzQ5YKoOAgUgB/sEiXVUQRLCLdMJetgWHLbu8ruZM8SLWMiqbgIBUhrA1yhciWVUDHYBFmVhHJhrlglrHwAu/Y1nQzz1kOW5LKoOAgZS/MDXJdUawpqKL1ZeozhMaDu5k74JjUpVSBOlM2IRVDwUHKXiw3yBHdFFSx+jLl7ODqP4kPNM1UYRQcRIqg0CkoKCyziOKs6qwdXP0n8YGmmSqQgoNUpFLUIAoRZr2iJGdVpzvTWyqKgoNUnLjUIAqRLVjkOta0mGdVZ23Sp2mmiqLgIBUnrjWIQrTvHPtyhxYf2TKJqFp/Z23Sp2mmiqPgIBJzwRYf2TKJKFcyBZv0dWi18alr2oKFppnKntpnSMV78u6/sKNpH4Nqe6euxa0GUYhkJjFuaFuL70yZRDKDOP7hX4U6hoVvLaT+nfrU44Zt3t/HumMT3Rm2rmTm7veY1bet5YemmUpL7TNEAsqxBpFNIZkERDPNFFzJFAwWjd2Ao/szK/mzp6aZyo4yB6k6UTbqK4VgJuHPIt57/An2LF6cujfZwC/sTCIoOe00esBo78LWlcx0vZg19z8jfV/JTJmDSJXxZxLBLKJYG+aCghvoGuwADXaA+oe8f5tmHncOs867pyhjkc5RcBApc8GVTblEtZrJr8O003Nfp37zH4FEoNjyPPU6NyLWFBykKoV1ilwcZVv2WrS+TAGzzruH5D/9Cxf8FfX2Qeo5nV0dTwoOUnXCOkUujnIVq0s1zeQ3i97M2rIO3HbAOzeisXmlTqCLGQUHqTr+TXLluEEum+BO63w20EXRlymrwIa5mbvfg6P7px43bGugYVtDu9VPChbFp+AgUsFyZRJx6Ms0a8FnmbVlZSqTWMgA77jShGCwUKAoDgUHqXpxb9LXFbkzieFcNP9fuHLKiJJMMQEdMolZW9cxa8ipMMs7u9q/h0L1ieJRcJCqVmkb5HLJlUkUYyVTBzk6vOZs26FMIhLaBCfiU2mtNnLxb6Cb9PoSpm1awbF9egLF2zDXwYLPejuqh5zadi3ReiNn2w4ULDKJ1SY4M5sB3At0Ax5wzt0eeP4bwN8Ah4Bm4Frn3IYoxySSTTVnEg8PPp3lE6bx+JfOBEqzkgnI2uE1Z9sOTTuFJrLMwcy6AW8B5wJNwDLgCufcGt89nwaWOuf2m9nfAdOcc3NyvbYyBymWasokknUIf3BIZg9JRZlmCkpOM839Xc5bk9NNC2YsiHJEZSlOmcNkYK1z7h0AM3sMuAhIBQfn3H/47v8z8IUIxyNSsGrLJPzF6kl9xzJt2AFqEs8VbTVTOv6DhHJ0d816IJHkLcrgMAzY6HvcBEzJcv91wNOZnjSzecA8gBEj8jtXV6Sr0h0cVKm7q4PF6lhOM+Xo7pr1QCIpSJTTSrOA851zf5N4/EVgsnPuK2nu/QIwH/gr59yBXK+taSUpldV/3MRbr2wDSE03VUp316Bgt9drFt5B7a4m+o4fB5THFFO7zrBUdyYRp2mlJmC473EtsDl4k5lNB75NnoFBpJSCu6sreY9EMJN4dsgnOB8YB+xftoz9y5a1awdetGCR51nVwc6wyiQKE2Xm0B2vIP0ZYBNeQfpK59xq3z2nAb8BZjjn3s73tZU5SBz4swiorkwiuOx1/7JlAHzsjDNS90cSLBoWwMrftD1OLnlVJpFTbDIH59whM5sPPIu3lPVB59xqM7sVaBA7HLMAAAsOSURBVHDOLQL+FegNLDQzgHedc38d1ZhEwpSrHgGVm0kE6xHBQ4UiyyxybJjLRplEYbQJTiQk1ZRJ5DrHumgn0GXZMJdLtS17jU3mIFJt0mUSlaozrcEj6f6aZcNcPtSKIzNlDiIRCW6gq6QppqBCMolI6xPBTCJLFpFud/XoAaMrNpMoNHNQcBCJSDUte31k6bs8tWJT6nEyUCRrEn7p6hMQUrDwF6w3vOj9evzZbc9nCRbBgnWlZREKDiIxVE1tOCB3JuEXWbAocGVTsDV4pWURCg4iMVRNxWooLJMIyhUsOp1VFFC8rsRitYKDSBmo9kwiUxaRjj9YdGnVUzCTyDLtVIl7IrRaSaQMVFtDP//qpuDKplz8K5+Cq56ggEwiuEci3bRT4j7tiVDmIBIL1ZRJFFKPCIq0mJ1lpVMlZBLKHETKUDVlErn2SGQT3D+Ra2d2QYEiS/fXaswklDmIxJD2SORfk/ALbT9FjuJ1OS57VeYgUgH8mcTmt99n89vvt1vtVEnBoiuZRJA/s+hSv6fgzusNL3pfiRrFTPZB76OBys0ilDmIxFxwGezmt98H4LiT+qWuVVKwCDOT8OtSv6cseybKZdmrlrKKVLhcwaLcA0VX9kgUoktnZPumnebaNhq7wejBbVNQcZxm0rSSSIULNvgLtulI3lOurpwyol2WkMwkkmdbh5FFgBcI/AoqZvumnWbu3Ay9ekGrV8Ru5CPYuzV2waFQyhxEKkglLon1ZxJRZRHQhWJ2YMpp7sF3aDyyB6OHTEpdi0MmoWklkSpW6fWJqOoRQV3ZT7FwwV9Rf3A7HNkLgAbzTj+uO7bt3+VSBAsFBxFJqbRgEaxHLF23C4ApowakrhUjWGQtZgcyiYW7XqW+Vy+o8VY3BYNFsQKFgoOIZFRpxexcwSKKQAEdi9lZ6xNZgkUqULieqednHncOs867J/QxKziISN4q7cwJf7BIl1UElWSznS9YLGQf9fZB21MRBgsFBxHplEorZgeziqCwpqTC7Pe08LmvU7/5j6nHwWDRlUCh4CAinVJp9Ylc8qlf+IXZHNAvW+DwB4u0WYXrxSzagnm2k+4UHEQkFJVWn8glW6bRlSwjGCz8cgUOaAseubIKWnYz84MPmDXgtLbf7AsWCg4iEgl/sEiXVQRVUvAoVpYRlC14NH/YzIvjjNfO8vpwNWzz/k3MFCzs2noFBxGJVjCrCMoVPMo9cBSaZQRFkXU0f9jMzg93pp7fe3AvL447gt0Tj4KW3Tz0d2sVHESktLIFj2qrZQTlEzz8MgWSfLOOdz/eh70H9zLr928qOIhIfOVT+A6qpuDh15VAElxuO65RwUFEykihU1SVFChyCSuQzFjyCNc+9mMFBxGpHIUWwqMS96CUK5A88eWz1LJbRCqHv0V5riwjKulO44ubo4DL6Znx+ScKfD0FBxEpG8GzLIqlVEGplBQcRERyKFVQCtW3Crv9iGhGISIi5SzS4GBmM8ys0czWmtmNaZ7vaWaPJ55famYjoxyPiIjkJ7LgYGbdgPuAC4BxwBVmNi5w23XAe865E4F7gDuiGo+IiOQvysxhMrDWOfeOc+4j4DHgosA9FwG/THz/G+AzZmYRjklERPIQZXAYBmz0PW5KXEt7j3PuELAbGJjuxcxsnpk1mFlDc3NzBMMVEZGkKINDugwguOMun3u8i87d75yrc87VDR48uMuDExGRzKIMDk3AcN/jWmBzpnvMrDtwNLArwjGJiEgeogwOy4CTzGyUmR0JXA4sCtyzCLg68f1lwAuuHPt5iIhUmEh7K5nZTOCHQDfgQefc98zsVqDBObfIzGqAh4HT8DKGy51z7+TxunuBxsgGHq1BwI5SD6KTNPbS0NhLo9LGfrxzLu85+bJsvGdmDYU0kIoTjb00NPbS0NhLI4yxa4e0iIh0oOAgIiIdlGtwuL/UA+gCjb00NPbS0NhLo8tjL8uag4iIRKtcMwcREYmQgoOIiHRQVsEhVwvwuDGz9Wa20sxWmFlD4toAM/u9mb2d+LV/qceZZGYPmtl2M1vlu5Z2vOb5UeL/xetmdnrpRp5x7Deb2abE578ise8m+dxNibE3mtn5pRk1mNlwM/sPM3vDzFab2VcT12P/uWcZe+w/98RYaszsFTN7LTH+WxLXRyWOEHg7caTAkYnrsTliIMvYHzKzdb7PfmLieuF/bpxzZfGFt5Huv4ETgCOB14BxpR5XjjGvBwYFrt0J3Jj4/kbgjlKP0ze2qcDpwKpc4wVmAk/j9cf6JLA0hmO/GfhWmnvHJf789ARGJf5cdSvRuIcCpye+7wO8lRhf7D/3LGOP/eeeGI8BvRPf9wCWJj7TJ/A25AL8FPi7xPf/C/hp4vvLgcdjOPaHgMvS3F/wn5tyyhzyaQFeDvxtyn8JXFzCsbTjnPsvOva2yjTei4BfOc+fgX5mNrQ4I+0ow9gzuQh4zDl3wDm3DliL9+er6JxzW5xzf0l8vxd4A69bcew/9yxjzyQ2nztA4jPcl3jYI/HlgP+Bd4QAdPzsY3HEQJaxZ1Lwn5tyCg75tACPGwc8Z2bLzWxe4tqxzrkt4P3lAo4p2ejyk2m85fL/Y34ijX7QN4UXy7EnpilOw/spsKw+98DYoUw+dzPrZmYrgO3A7/Gymfedd4QAtB9j3kcMFENw7M655Gf/vcRnf4+Z9UxcK/izL6fgkHd77xj5lHPudLzT8K43s6mlHlCIyuH/x0+AjwMTgS3A3YnrsRu7mfUGfgt8zTm3J9utaa7Fbexl87k751qdcxPxukZPBsamuy3xa6zGHxy7mZ0C3ASMAc4ABgA3JG4veOzlFBzyaQEeK865zYlftwNP4v3h25ZM5xK/bi/dCPOSabyx///hnNuW+At0GPg5bVMYsRq7mfXA+8f135xz/564XBafe7qxl8vn7uecex9Ygjcf38+8IwSg/RhjecSAb+wzElN9zjl3AFhAFz77cgoO+bQAjw0z62VmfZLfA+cBq2jfpvxq4KnSjDBvmca7CLgqsQrik8Du5DRIXATmVC/B+/zBG/vlidUno4CTgFeKPT7wVpEAvwDecM79wPdU7D/3TGMvh88dwMwGm1m/xPdHAdPx6ib/gXeEAHT87GNxxECGsb/p+4HC8Gol/s++sD83paq2d+YLr+L+Ft684LdLPZ4cYz0Bb2XGa8Dq5Hjx5ij/ALyd+HVAqcfqG/OjeNMAB/F+0rgu03jx0tT7Ev8vVgJ1MRz7w4mxvZ74yzHUd/+3E2NvBC4o4bjPxkvvXwdWJL5mlsPnnmXssf/cE2OZALyaGOcq4DuJ6yfgBa21wEKgZ+J6TeLx2sTzJ8Rw7C8kPvtVwK9pW9FU8J8btc8QEZEOymlaSUREikTBQUREOlBwEBGRDhQcRESkAwUHERHpoHvuW0Sql5kll5QCDAFagebE4/3OubNKMjCRiGkpq0iezOxmYJ9z7q5Sj0UkappWEukkM9uX+HWamf2nmT1hZm+Z2e1m9vlEv/2VZvbxxH2Dzey3ZrYs8fWp0v4XiGSm4CASjk8AXwVOBb4InOycmww8AHwlcc+9wD3OuTOASxPPicSSag4i4VjmEr1qzOy/gecS11cCn058Px0Y5zsCoK+Z9XHeWQgisaLgIBKOA77vD/seH6bt79kRwJnOuQ+LOTCRztC0kkjxPAfMTz5Inu8rEkcKDiLF8/dAXeKUrjXAl0s9IJFMtJRVREQ6UOYgIiIdKDiIiEgHCg4iItKBgoOIiHSg4CAiIh0oOIiISAcKDiIi0sH/B9t++M+pH/H7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv.iloc[:, :5].plot(drawstyle='steps-post')\n",
    "plt.ylabel('S(t | x)')\n",
    "_ = plt.xlabel('Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation\n",
    "\n",
    "The `EvalSurv` class contains some useful evaluation criteria for time-to-event prediction.\n",
    "We set `censor_surv = 'km'` to state that we want to use Kaplan-Meier for estimating the censoring distribution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = EvalSurv(surv, durations_test, events_test, censor_surv='km')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Concordance\n",
    "\n",
    "We start with the event-time concordance by [Antolini et al. 2005](https://onlinelibrary.wiley.com/doi/10.1002/sim.2427)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6523591504514816"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev.concordance_td('antolini')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Brier Score\n",
    "\n",
    "We can plot the the [IPCW Brier score](https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0258%2819990915/30%2918%3A17/18%3C2529%3A%3AAID-SIM274%3E3.0.CO%3B2-5) for a given set of times.\n",
    "Here we just use 100 time-points between the min and max duration in the test set.\n",
    "Note that the score becomes unstable for the highest times. It is therefore common to disregard the rightmost part of the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xV9f348dc79+ZmkEFCmGGEEfYmIA5ABdyKtThaB1oVtdWfs1+1Q1s7tVpt3dZRq21ddaBiBWSIAyTslZCwkwBJgGwy7+f3xz3QEG7ITXJvTu7N+/l43EfuPePe9+GQ+85nizEGpZRSqqEwuwNQSinVPmmCUEop5ZUmCKWUUl5pglBKKeWVJgillFJeOe0OwF+SkpJMSkqK3WEopVRQWb16daExpqu3fSGTIFJSUkhPT7c7DKWUCioisruxfVrFpJRSyitNEEoppbzSBKGUUsorTRBKKaW80gShlFLKK00QSimlvNIEoZRSyquQGQehVHNV17rZmFtMQWklB8urKaqoYVL/RCamJNodmlLtgiYI1aGUVtawOCOfhVsOsCyzgNKq2hOOuXxCb352wTASOrlsiFCp9kMThAp5xhjW7S3i39/t4eP1+zhSU0dSjIsLRvXkrKHd6J0QRZcYF9EuJy8s287fvtzBFxn5/N+5Q7h0XDKR4Q67L0EpW0iorCiXlpZmdKoN1dDKHQf5zadb2JRbQrTLwcWje3F5Wm/G9U3AESZez8nYX8LP3t/Imj1FxEY4uWBUTy4bn0xaSmKj5ygVrERktTEmzes+TRAqFBWUVvGH+Vt5f20uyZ2juO3Mgcwa24vYyHCfzne7DSt2HuT9NbnM37iPiuo6EqLDOX1QElNTu3LqwC70TohCRBOGCm6aIFTIq6lzszmvhHV7DrNubxFfZORTWVPH3KkDuP2sVKJcLa8mqqiuZdHWfJZlFrA8q4D80ioAkmJcjOndmbF9OjOqdzyjkuPpEhPhr0tSqk2cLEFoG4QKCT/6+yqWZxUC0C02gmmDu3L3zMEM7BrT6veOdjm5ZEwvLhnTC2MMmQdKWbXrMOv2FLFu72G+yMg/dmyv+EiuntyPW6YOwOnQXuQquGmCUEFvY04xy7MKmTt1ANeflkLP+MiAVf2ICEN7xDG0RxzXTu4HQEllDZtzS9iUW8zX2wv50+eZLNhygCcuH82gbrHN/ozqWjcllTWUVdZSVlVLQicXyZ2j/H0pSjVJq5hU0Lv/vQ3MW5/Hyp9PJ87HNoZA+nTDPn7x4UbKq+u48Yz+pPVLYFjPuGOJq7rWTVlVLW5jiHCG4XKGcbi8hkVbD7BwywG+3X6Q6jr3ce85tEcs5wzvzvRh3RmZHK+N5cpvtIpJhaziiho+Wp/L98Ylt4vkAHDh6J5M6p/ILz/cxPNLtx/b3snloNZtqKp1N3puvy7RXHtqP/p1iSYmwkmnCCd7DlawcOsBnlmSzV8XZxMT4WRc385MTEnE6RDyio6QV1SJI0z4/vhkpg/rTrhWbyk/0AShgtp/1uRQWePmGqu6p73oGhvBC9dOoLSyhsz9pWzdV8L2gnIinGHERnq++B1hntJEVa2bCGcY0wZ3ZVC3GK/VYzdPHcCh8mqWZxWwatch0ncd5slF2zAGEqLD6RkfxaHyahZuOUDX2AiumtiHn5w1SMdwqFbRBKGCljGGN1fsZnzfzozoFW93OF7FRoaTlpJImh+m70js5GLW2GRmjU0GPKPCHWFCtMvza1znNizNzOefK/fw9OJsnGFh3DkjtdWfqzouLYeqoPXN9oPsKCzn2lPbV+mhrcRGhh9LDgCOMGH6sO68ev1Ezh7ajX98u4vKmjr7AlRBTxOEChrlVbVsyCliX/ERaurcvPHtbhI7uTh/ZE+7Q2t3bprSn4Pl1XywNtfuUFQQ0yom1a4s2nKA17/dxRmDkrg8rQ+JnVxU1dbx5oo9PLskm0Pl1QCIgDFw67SBWs/uxakDujAyOY6Xl+/gyrQ+hGmvJ9UCmiBUu2CM4bml23l8QSYJ0S6WZxXyxMJtnDO8O2v3FJFbdIQzBiVx1aQ+lFbWcqCkkuIjNdw0pb/dobdLIsLNUwZw51vrWJKZz/Rh3e0OSQUhTRCqTVTV1pFfUkWfxOgT9lXW1HH/fzbw0bo8Lh7Ti8e+P5o9hyp4c8VuPlyby4CunXhs9mhOH5RkQ+TB64JRPXn0swxe+nKHJgjVIpogVMB9k13ILz7cxM6D5fz4zIHcNWPwsX76W/JKuPfd9WzdV8J95wzmJ2cNQkQY0iOW31w6kkdmjdAJ8Voo3BHGj87oz28/3cqGnCJG9+5sd0gqyGgjtQqYwrIq7n57HT98eSW1bsNFo3vx7JLtXPnit+w+WM4zi7OY9exXFJRW8er1adx+duoJyUCTQ+tcObEPsRFOfv7BJjL2l9gdjgoyOtWG8jtjDO+tzuG3n26lorqWW6cNPDZo6+P1efzs/Y2UVddiDFw0uie/mTVSV28LoE825PHzDzZRWlnDNZP7cc/MwXSO1n9v5aHTfas2s+dgBT/7YCNfZRcyMSWBP1w26oQJ6/YequDxBZnMHN6di0b3sinSjuVweTV/XriNf67cjdMRRpdOLuKjwomLCmdqahJXTepLkk5V3iFpglABk1d0hLdX7WV7QRnbC8rJzi8lwung/vOHcvWkvtq9sp3Zuq+ED9bmcqi8mpIjNRworWL93iJcjjAuHN2TS8b2ol9iNL06R2n34Q5CJ+tTAWGM4bY3V7Mht5g+CdEM7NqJqalJzDkthV46PXW7NKxnHMN6xh23LTu/jDdX7Oa91TnHDazrkxjFC9dMaLfTmKjAC2gJQkTOA/4COICXjTF/bLD/HuAmoBYoAH5kjNlt7ZsD/MI69LfGmNdP9llagmh7SzLzueG1VfzhslH8YFJfu8NRrVRWVcum3GJyDx8hr+gI/1ixm4TocD6+4wwinFqaCFW2lCBExAE8C8wEcoBVIjLPGLOl3mFrgTRjTIWI3AY8BlwpIonAw0AaYIDV1rmHAxWvah5jDH9ZlEVy5yi+P7633eEoP4iJcDJ5QJdjr0cmx3PD31fx1KIs7j9vqI2RKbsEspvrJCDbGLPDGFMNvAXMqn+AMWaJMabCerkCOPpNcy6w0BhzyEoKC4HzAhiraqYvswpZt7eIH581EJdTe0uHorOGduPKtD68uGw7a/bo32YdUSB/s5OBvfVe51jbGnMj8FlzzhWRuSKSLiLpBQUFrQxX+cpTethGr/hILp/Qx+5wVAD94qJh9IyP4r531nOkWmeG7WgCmSC8dV/x2uAhItfgqU76U3PONca8ZIxJM8akde3atcWBqub5Ovsga/YUcdtZg7T0EOJiI8P50+zR7Cgs50+fZ9odjmpjgfztzgHq/3nZG8hreJCIzAB+DlxijKlqzrmq7Rlj+MsX2+gRF8kVadr20BGcNiiJq0/py9+/2cn6vUV2h6PaUCATxCogVUT6i4gLuAqYV/8AERkHvIgnOeTX2/U5cI6IJIhIAnCOtU3ZbNWuw6zadZhbpw3Qni0dyP3nDyUpJoIH3t9ITV3ja2qr0BKwBGGMqQVux/PFvhV4xxizWUQeEZFLrMP+BMQA74rIOhGZZ517CPgNniSzCnjE2qZs9sKy7SR2cnHlRO3W2pHERYbzyKwRbN1Xwqtf7bQ7HNVGAjpQzhgzH5jfYNtD9Z7POMm5rwKvBi461VwZ+0tYnJHPPTMHE+XS0kNHc+6IHswc3p0nF23j/JE9SU6IYlNuMSt3HqS8qg5jDG4DE/olcNbQbnaHq/xAR1Irn724bAfRLgfXddA1oDs6EeGRWSOY8cQyrn5lBaWVtRRV1NTb7+ld4jbw4PlDmTt1gM7GG+Q0QSif5ByuYN76PK4/LUVnAu3AesZH8fDFI3hmSTYzhnVnSmoSpw9KoksnFyJCVW0d9727gT98lsGBkip+ceEwnY8riGmCUD55eflOwgRd4lNxxcQ+XDHR+/iXCKeDv1w5lq4xEbz69U4Kyqp46sqxODRJBCXtxK6adKi8mrdW7WHW2GR6xuskfOrkwsKEX140jJ+eO4SP1+fxyQbtoR6sNEGoJj23JJuqWje3ThtgdygqSIgIt00bSGq3GJ5fup1QWVago9EEoU5qR0EZf/9mF1em9Tlh4R+lTiYsTLh12kAy9peyJDO/6RNUu6MJQp3U7+dvJTLcwb3nDLE7FBWELhnbi+TOUTy/dLvdoagW0AShGrU8q4BFW/O5/exBdI3V5ShV84U7wrh5Sn9rBL6OdQ02miCUV7V1bn7zyRb6JkZzw+kpdoejgtiVE/uS2MmlpYggpN1clVf/XrWXbQfKeOGaCTrnkmqVKJeDG05L4YmF23h5+Q6qat3sL66k1m0Y2LUTg7rFMKRHrPaQa4c0QagTVNe6eWZxFpP6J3LuiO52h6NCwHWnpvDS8h389tOtAMRHhSPCcSOxbzg9hZ9dMIxwh1ZstBeaINQJPl6fx4GSKh6bPUanSlB+ER8dzqJ7plFZU0e32Mhjc3kdLKsiO7+MTzfu47Wvd7Ept5hnfziebnGRFFVUM3/jfnIOV3DPzME4NXG0OU0Q6jjGGP62fAdDuscyNTXJ7nBUCOkeF3nCti4xEXSJieCUAV1IS0nk/vc2cOHTXzE6OZ4vswqoqfOMn+gWG8H1px8/in9f8RHKq2pDovv1wbIqYiKd7a46V1OyOs5X2YVk7C/lpin9tfSg2tQlY3rx4U9OJzHaxea8Em44vT+f3HEGU1KT+PPCbRSWVR07triihstf+JZrXv4Otzt4B+FtLyjj7rfXMfF3i3h2SftrxNcShDrOS1/uoGtsBJeM7WV3KKoDGtIjls/vnoox5tgfKA9fPILznvqSP/03k0dnj8btNtzzzjpyDh8BYFNeMaN7d7Yz7GbLOlDKc0u389G6XFzOMKJdTjL2ldgd1gm0BKGO2bqvhOVZhVx/Wkq7K+qqjqV+6XVQtxh+dEZ/3k7fy7q9Rbz45Q6+yMjnzumphAks2nLAxkh9Z4xh1a5D3PT6KmY++SWfbdrHTVMG8NX9Z5OWkkBu0RG7QzyBliDUMS8v30lUuIOrT9HV4lT7csfZg/hgbS53vrWWvYcquHB0T+6akco32wtZtDWfe4JgpP/v52/lb8t3khAdzl0zUrnu1BQSO3mmzk/uHNUu1/vWEoQCYPfBcuatz+WKtN663oNqd2Ijw/nZBUPZfbCClKROPPr90YgIM4Z1Z8u+knb513d9R6rr+NfKPZw3ogffPDCdu2YMPpYcAHp1juJwRQ0V1bU2RnkiTRCKOrfh3nfWExnu4NYzB9odjlJeXTo2md9cOpLXrp9ITISn8mPGcM84nS+2tu9qpiWZ+ZRX13Hdqf28LtfbO8EzSDD3cPtKdJogFC8v30H67sP8+pIROppVtVsiwrWT+9GvS6dj2wZ2jWFAUicW1muHqKqt44bXvuOx/2a0+TTjxhiqa90nbP94fR5JVndeb5I7e37vctpZSUgTRAeXub+UJxZs49wR3fneuGS7w1Gq2aYP68aKHQcprfSMyn7sv5ksySzguaXb+dPnmW0ay2tf72Li7xZxoKTy2LayqloWZ+Rz4ageja6sl2yVIPI0Qaj2orrWzd1vryMuysnvvzdKxz2ooDRjWHdq6gzLswpZmpnPK1/t5LpT+3H1KX15bul2nlmc1SZxVNbU8dzS7RQfqeHxeolp0ZYDVNW6uXhM413Hu8VG4gyTdlfFpL2YOrC/fpHFln0lvHTtBLrE6HTeKjhN6JdA5+hw3knfy6bcEoZ0j+VnFwzD5QjjSHUdjy/YhtMRxtwpAwgL4NrY76/JpbCsikkpiby3Joc5p6UwMjmej9fn0Ss+kvF9Exo91xEm9IiPbHeN7VqC6KBW7z7Ec0uzmT2hN+eM6GF3OEq1mNMRxtlDurE0s4CSyhr+8oOxRIY7CAsTHps9mgtH9eSPn2Vw4dNfsTjjQEDaJercnilqRiXH87c5aSREu/jdp1spqqjmy6wCLhzds8nklNw5SquYlP3Kqmq5++31JCdE8fDFw+0OR6lWO/pHzs/OH8rQHnHHtjsdYfz1B+N46sqxVFTX8qO/pzP7hW/9/pf6gs372VlYzi3TBhAfFc7dM1L5dsdB/u+9DdTUmZNWLx2V3Dmq3VUxaYLogB75eDM5hyt48oqxxEaG2x2OUq127ojufHLHGcw5LeWEfY4w4dJxySy6Zxq//94oMveX8sB/NvitJGGM4YVl2+nXJZrzR/YE4AeT+jKoWwwLthygX5doRiXHN/k+yQlR7C+ppKbuxF5QdtEE0cF8vnk/76TncNuZA0lLSbQ7HKX8QkQYmRx/0o4W4Y4wfnhKX3567hCWZxXyyYZ9fvnsb3ccZH1OMTdPGXCsl5LTEcbPLxwGwMWje/nUASS5cxRuA/uLK5s8tq1oguhACsuqePD9jYxMjuPO6YPtDkcpW1wzuR+je8fzyCdbKD5S0/QJTXh+6XaSYlzMntD7uO1nDu7Ka9dP5JZpA3x6n/bY1VUTRAfy0EebKKus5ckrxuJy6q1XHZMjTPjdpaM4WFbFEwtaN07iu52HWJ5VyI1nDCAy/PgR0iLCWUO7+VyN28saLNeeejLpt0QH8cmGPOZv3M9dM1NJ7R78C6wo1Rqjesdz3akpvLFiN+taOEmeMYbH/pvhWczIS9tHcx0dTd2eGqo1QXQAhWVVPPTRZsb0jmfuFN+Ku0qFunvPGUy32Ajuemst+SXNr/dfmllA+u7D3DE91ev8Ss0VGe4gKcbV6hJEZU0dW/00gaEOlAtxxhh++aGnaunxy8four5KWWIjw3nu6glc+8pKrn55JW/NnezzgFG32/DY55n0TYzmyrQ+fospuXNUi77YN+UW8/TiLDbl/i8xdI4OZ8FdU+nmZalXX+m3RQjLLTrC7f9ey2ebtGpJKW8m9EvglTkT2XOogmtf+Y7iCt8arT/ZuI+t+0q4Z+Zgv7bn9WpmgthfXMl9767n4me+4rudh5jUP5F7Zg7mj5eN4kh1HQ++v7FV3Xm1BBGCKmvqeHHZDp5flo0xcNeMVG6ZqtN4K+XNqQO78NJ1adz8ejrXvfYdb8+dfEKDc301dW7+vCCToT1iucSHAXDNkdw5iiWZ+cctuepNUUU1L325g9e+3kWd2zB3ygB+fNYg4qP+1yBeUV3HI59s4d30HK6Y2LJSTkBLECJynohkiki2iDzgZf9UEVkjIrUiMrvBvjoRWWc95gUyzlBSUFrFZc99w5OLtjF9aHe+uHcad80Y3OgskkopmDa4K09dNZb1e4t449vdJz328c8z2XWwgp+eO8TvczslJ0RRWePmUHm11/1lVbX89Ysspjy6hOeXbWfmcM/v+IMXDDsuOQBcf1oKkwck8sgnW9h7qKJF8QSsBCEiDuBZYCaQA6wSkXnGmC31DtsDXA/c5+UtjhhjxgYqvlC091AF17yykvySKl6Zk8b0Yd3tDkmpoHHBqJ5MG9yVpxdncXkjKyt+tnEfL365g2sm9w3I71f9rq4N20OqauuY/fw3ZOwvZebw7tx7zuDjphVpKCxM+NPsMZz/l+Xc++56bq7XQWViSoJPK0cGsgQxCcg2xuwwxlQDbwGz6h9gjNlljNkAtJ+x5UEqY38J33/+G4oqanjzplM0OSjVAg9eMJTSqlqeXZJ9wr7s/FLue3c9Y/t05pcXBWYOs5N1dX1qURYZ+0t54ZoJ/O26tJMmh6P6JEbz0MXD+W7nIW7+R/qxxxMLtvkUTyDbIJKBvfVe5wCnNOP8SBFJB2qBPxpjPvRncKFkzZ7DXP/qd0S7nLx766kM1sZopVpkaI84Zo/vzevf7Oa6U1PokxgNeKp2bnljNZHhDp6/ZjwRztZ3a/Xm2NKjDRqq1+45zIvLtnNlWh/OG9m82ZevSOvDKf0TKa30rHc99x/pPo8gD2QJwlvlXHOa0/saY9KAHwJPicgJrawiMldE0kUkvaCgoKVxBrWvswu55uWVJHZy8d5tmhyUaq17zxlCWBg8viCT6lo376zayyXPfMXOwnKe/uG4gC7LGx8VTieX47gEUVlTx73vrqdHXCQ/v2hYi963X5dOjEyOZ2RyPDGRTq/LonoTyBJEDlC/6bw3kOfrycaYPOvnDhFZCowDtjc45iXgJYC0tLS2XXy2HVi05QA//tcaBiR14h83TqJbbMv7OyulPHrER3LTGQN4Zkk2324/SH5pFcN7xvHynDROG5gU0M8WEU9X13pVTE8syGRHQTlv3DiJOD/MvhzhdFDt44yxgUwQq4BUEekP5AJX4SkNNElEEoAKY0yViCQBpwOPBSzSIPRVViG3vLmakcnxvH7DRJ8anJRSvrll2gDmb9xH97hIHr98DFNSk9psSd7khCiy8st4ZnEWC7ccYH1OMVef0pcpqV398v4uZ5j9JQhjTK2I3A58DjiAV40xm0XkESDdGDNPRCYCHwAJwMUi8mtjzAhgGPCiiLjxVIP9sUHvpw7N7Tb89tMt9EmI4p83nUJMhA5nUcqfYiPDWXzfmbZ8dt/EaJZmFvD4gm2M6dOZ/ztvCD86vb/f3t/laAcJAsAYMx+Y32DbQ/Wer8JT9dTwvG+AUYGMLZh9tmk/GftLeerKsZoclAoxt04byNg+nTl9UBLdWzFNRmNczjCK/NVILSLRIvJLEfmb9TpVRC5qZYyqherchicXbSO1W4xPyxgqpYJLr85RXDa+d0CSAzSvismXXkyvAVXAqdbrHOC3LQtNtdYnG/LIzi/T0dFKqRZxOcOoqq3z6VhfEsRAY8xjQA2AMeYI3ruwqgCrrXPz1KIshvaI5fxm9oVWSimAiGa0QfiSIKpFJAprDIM1HqGq5eGplvpgbS47C8u5Z+Zgv88Bo5TqGPzdi+lh4L9AHxH5J54up9e3ODrVIptyi/n9/K2M7h3PzOE6jYZSqmUinGH+GQchno6/GcBlwGQ8VUt3GmMKWxuk8t2aPYeZ8+p3xEWG8/QPxrVZf2ylVOjxWwnCGGNE5ENjzATgU38Ep5pnxY6D3Pj3VSTFRvCvmycfm8xLKaVawt9VTCtEZKI1ZkG1kez8Ml7/ZhfvpO+lT2I0/7zplIB1e1NKdRwuh4Nat6HObZrsCelLgjgLuEVEdgPleKqZjDFmdOtDVQ1tLyjjV/M2szyrEJcjjEvG9uLB84f6vFauUkqdzNElUqtr3US5Tj4rrS8J4nw/xKR89If5W1m3p4j7zhnMVZP6kqSJQSnlR35NEMaY3SIyBphibVpujFnf2iDViY5U17E8q5AfTOrL7Wen2h2OUioEHU0QVXV1wMlnh/Vlqo07gX8C3azHmyJyR6ujVCdYnlVAVa1bu7EqpQImol4Joim+VDHdCJxijCkHEJFHgW+Bp1seovJm4ZYDxEU6mdQ/0e5QlFIhqjkJwpeR1ALUn7ijDp1qw+/q3IbFGfmcNbQb4Y5ALvSnlOrIXNb3iy+D5XwpQbwGrBSRD6zXlwKvtDQ45d2aPYc5WF6t1UtKqYA61gZR44cEYYz5s7Xk5xl4Sg43GGPWti5E1dDCLQcIdwjTBvtn1SillPLmWC8mf5QgRGQysNkYs8Z6HSsipxhjVrYyTmUxxrBwywFOHZhErB/WnFVKqcYcq2LyUxvE80BZvdfl1jblJ9sLythZWK7VS0qpgHP5u5HaGGOOvjDGuAnwUqUdzYItBwCYOUwThFIqsCKcnsFxVX5KEDtE5P+JSLj1uBPY0boQVX0LtxxgdO94esTrXEtKqcBqThuELwniVuA0IBfPcqOnAHNbHp6qb39xJWv3FHGOVi8ppdqAXwfKGWPygataHZXyav7GfQBcMKqnzZEopTqCY91cfViX2pepNh4TkTireukLESkUkWtaH6YC+HTjPob3jGNA1xi7Q1FKdQD+7sV0jjGmBLgITxXTYOCnrYhPWfKKjrB692EuHK2lB6VU2/B3L6ajHfMvAP5tjDnU4sjUcY5WL12o1UtKqTbSnAThS3fVj0UkAzgC/FhEugKVrQlQeXy6cR8jesWRktTJ7lCUUh2EM0wQ8VMvJmPMA8CpQJoxpgaoAGa1OsoOLudwBWv3FGn1klKqTYkIET6uS+3TgDdjzOF6z8vxjKZWrfDZxv2AVi8ppdqeyxHmt4FyKgA+2biPUcnx9Oui1UtKqbblcjpanyDEo4/folIA7D1Uwfq9Wr2klLKHr1VMJ00Q1hxMH/orKOWZufWPn2UQ7hCtXlJK2cLlDPPbVBsrRGRi60NSAB+szeXTjfu4e+Zg+iRG2x2OUqoDcjnCqPZhJLUvjdRnAbeKyC48jdOCp3AxulURdkB7D1Xw8EebmZSSyC1TB9odjlKqg3L5sRfT+a0PR9W5Dfe+sx4DPHHFGBxhuqy3UsoeEf6qYjLG7Ab6AGdbzyt8OU/9jzGGvyzaxne7DvHrS0Zo1ZJSyla+liB8mazvYeB+4EFrUzjwpi9BiMh5IpIpItki8oCX/VNFZI2I1IrI7Ab75ohIlvWY48vntUdVtXXc/58N/HVxNpeNS+ay8cl2h6SU6uBcTt/GQfhSxfQ9YBywBsAYkycisU2dJCIO4FlgJp5J/laJyDxjzJZ6h+0Brgfua3BuIvAwkAYYYLV17mGCSH5pJbe9uYbVuw9zx9mDuHvGYES0akkpZS9PI7V/EkS1McaIiAEQEV9Hdk0Cso0xO6zz3sIzRcexBGGM2WXtaxjpucDCoxMDishC4Dzg3z5+tu3Kq2r53rPfcLC8imd/OF7HPCil2g2/VTEB74jIi0BnEbkZWAT8zYfzkoG99V7nWNt84dO5IjJXRNJFJL2goMDHt24by7MKyC06wnNXa3JQSrUvfqtiMsY8LiIzgRJgCPCQMWahDzF4q0sxPpzn87nGmJeAlwDS0tJ8fe828cXWfOIinUxJ7Wp3KEopdRxfezH5OlnfQsCXpFBfDp7eT0f1BvKace6ZDc5d2szPt43bbViSmc/UwV0Jd2iHL6VU+xLhdLSuiklEvrJ+lopISb1HqYiU+BDDKiBVRPqLiAvPutbzfIz/c+AcEUkQkQTgHGtbUNiYW0xhWTXTh3WzOxSllDpBqwfKGWPOsH422V6ygZsAABBpSURBVGOpkfNrReR2PF/sDuBVY8xmEXkESDfGzLOm8PgASAAuFpFfG2NGGGMOichv8CQZgEeCaSW7LzLyCROYNlgThFKq/fFM993KqTZEJAzYYIwZ2ZIgjDHzgfkNtj1U7/kqPNVH3s59FXi1JZ9rt8UZBxjXN4HETi67Q1FKqRO4nGG4DdQ20Q7R1GyubmC9iPT1Z3Ch7EBJJZtySzh7qJYelFLt07F1qZtIEL40UvcENovId9RbSc4Yc0kr4gtZSzLyAbT9QSnVbrmszjNNtUP4kiB+7Yd4OowvMvLpFR/JkO4tarpRSqmAO1aCaG2CMMYsO/pcRJKAg9ZCQqqBypo6vs4u5LLxyTqlhlKq3YqwEkRTg+VO1s11sogsFZH3RWSciGwCNgEHROQ8fwYbKlbuPERFdR3Th3a3OxSllGqUy8cEcbISxDPAz4B4YDFwvjFmhYgMxTMn0n/9EmmIqHMbnv4ii7hIJ6cO7GJ3OEop1agIH6uYTtaLyWmMWWCMeRfYb4xZAWCMyfBXkKHkpS93kL77ML+eNYLIcIfd4SilVKN87cV0sgRR/8wjDfZpG0Q9W/eV8OeFmZw/sgeXjtX1HpRS7ZvL4fkjtjWN1GOsKTUEiKo3vYYAkX6IMSRU1dZx99vriI9y8bvvjdLGaaVUu9fqXkzGGK0n8cGTC7PI2F/KK3PSdOS0Uioo/K+K6eTTbehUo62wbFsBL365nR9M6sP0YdpzSSkVHHwdKKcJooXyio5w11trGdI9locuGmF3OEop5bOI8FaOg1CNq6lzc/u/1lBd6+bZq8cT5dLaOKVU8DhagmjNOAjViEc/y2DNniKe+eE4BnaNsTscpZRqFn+Mg1BefLfzEC9/tZM5p/bjotG97A5HKaWazddeTJogmunZJdkkxbh48IJhdoeilFIt4o+BcqqBzXnFLNtWwA2n99fR0kqpoKW9mALghWU7iIlwcs3kfnaHopRSLeZ0hBEmmiD8ZldhOZ9uyOOayf2Ijwq3OxyllGqVCKejyXWpNUH46KXlO3A6wvjR6Sl2h6KUUq3mcoZpCcIf8ksqeS89h9kTetMtTqehUkoFP5czTBup/eGNFbupdbuZO2WA3aEopZRfuBxhOpLaH77ZfpDxfRNISepkdyhKKeUXEVrF1Hp1bsOWvBJG9Y63OxSllPIbbYPwg+0FZRypqWNUsiYIpVTo0DYIP9iUWwzASE0QSqkQEuEMo6pGE0SrbMwtJircoZPyKaVCipYg/GBTbjHDe8XhCNOlRJVSocPl0DaIVqlzGzbnlWj7g1Iq5GgjdSvtLCyjorpO2x+UUiHH5XRoFVNrbMotAWBkcpzNkSillH9pFVMrbcwtJjI8jEHaQK2UCjEup46kbpWNucUM6xmH06H/TEqp0BLhDNPZXFvKfXQEtbY/KKVCkO1TbYjIeSKSKSLZIvKAl/0RIvK2tX+liKRY21NE5IiIrLMeLwQyTm92HSynrKpWG6iVUiHJl3EQzkB9uIg4gGeBmUAOsEpE5hljttQ77EbgsDFmkIhcBTwKXGnt226MGRuo+Jqy8egI6l6aIJRSocflCMOYkx8TyBLEJCDbGLPDGFMNvAXManDMLOB16/l7wHQRaRcj0jblFuNyhpHaXRuolVKhx+Vs+us/kAkiGdhb73WOtc3rMcaYWqAY6GLt6y8ia0VkmYhM8fYBIjJXRNJFJL2goMCvwR9toA7XBmqlVAiyO0F4Kwk0LNA0dsw+oK8xZhxwD/AvETlhMIIx5iVjTJoxJq1r166tDvgo97ER1Dr+QSkVmuxOEDlAn3qvewN5jR0jIk4gHjhkjKkyxhwEMMasBrYDgwMY63E25hZTWllLWr/EtvpIpZRqUy4fakcCmSBWAaki0l9EXMBVwLwGx8wD5ljPZwOLjTFGRLpajdyIyAAgFdgRwFiPszSzABGYOth/pRKllGpPIsIdTR4TsF5MxphaEbkd+BxwAK8aYzaLyCNAujFmHvAK8IaIZAOH8CQRgKnAIyJSC9QBtxpjDgUq1oaWbctndO/OJHZytdVHKqVUm/KlBBGwBAFgjJkPzG+w7aF6zyuBy72c9x/gP4GMrTFFFdWs21vEHWen2vHxSinVJiJsboMISl9mFeI2MG2IVi8ppUKX3Y3UQWlpZj6do8MZ07uz3aEopVTAaIJoJrfb8OW2QqakdtUV5JRSIc3uXkxBZ8u+EgrLqjhTey8ppUKcliCaaWlmPqDdW5VSoU8bqZtpaWYBI5Pj6BobYXcoSikVUFqCaIbiihrW7DnMmYO72R2KUkoFnCaIZvh6u6d765navVUp1QFEOJoeSa0JwrJ2z2FczjBGa/dWpVQHoCWIZtiQU8zwnnE+/aMppVSw0wThozq3YVNuMaN76+pxSqmOwREmTY730gQB7Cwso7y6jlG6/rRSqgNparCcJgg81UuAtj8opTqUiHBNEE3akFNMVLiDQd10/WmlVMehJQgfbMgpYmRynM6/pJTqUJpqqO7wCaK2zm2tP63VS0qpjkUTRBOy8suoqnUzpo82UCulOhatYmrCRquBWnswKaU6mqYm7OvwCWJ9ThGxEU5SunSyOxSllGpTWsXUhI25xYxMjidMG6iVUh1MhPPk8zF16ARRVVvH1n0ljNb2B6VUB6QliJPYtr+MmjrDaO3BpJTqgLSR+iTW5xQB6BxMSqkOSUsQJ7Exp5iE6HB6J0TZHYpSSrU5TRCNWLatgI/W5zIxJRERbaBWSnU8miC8+HTDPm56fRUDkmL4/WWj7A5HKaVs0VQbhLON4rCdMYaiiho+2ZDHw/M2M6FfAi/PmUh8VLjdoSmllC2ams015BPEwi0HePS/GeQePsKRmjrAs+7081dPIMrV9JqsSikVqiI6egnir19kUV5Vyw9P6UuvzlH0TYzmzCFdCW/iH0YppULdOSN6nHR/SCeIjP0lbMwt5uGLh3PD6f3tDkcppdqVkU3MQRfSf0a/l55DuEOYNTbZ7lCUUirohGyCqKlz8+G6XKYP7U5iJ5fd4SilVNAJ2QSxJCOfwrJqZk/obXcoSikVlEI2Qby3OoekmAjOHNLV7lCUUiooBTRBiMh5IpIpItki8oCX/REi8ra1f6WIpNTb96C1PVNEzm3O5xaWVbE4I5/Lxifj1N5KSinVIgH79hQRB/AscD4wHPiBiAxvcNiNwGFjzCDgSeBR69zhwFXACOA84Dnr/Xzy0bo8at1Gq5eUUqoVAtnNdRKQbYzZASAibwGzgC31jpkF/Mp6/h7wjHgmRpoFvGWMqQJ2iki29X7fNvZh2w6UMvPPywDYV1zJmN7xDO4e698rUkqpDiSQCSIZ2FvvdQ5wSmPHGGNqRaQY6GJtX9Hg3BP6qorIXGAuQFyvAaR2jwFgcPdYrp7c1z9XoZRSHVQgE4S3KVKNj8f4ci7GmJeAlwDS0tLMc1dPaG6MSimlGhHIFtwcoE+9172BvMaOEREnEA8c8vFcpZRSARTIBLEKSBWR/iLiwtPoPK/BMfOAOdbz2cBiY4yxtl9l9XLqD6QC3wUwVqWUUg0ErIrJalO4HfgccACvGmM2i8gjQLoxZh7wCvCG1Qh9CE8SwTruHTwN2rXAT4wxdYGKVSml1InE8wd78EtLSzPp6el2h6GUUkFFRFYbY9K87dNRZEoppbzSBKGUUsorTRBKKaW80gShlFLKq5BppBaRUiDT7jj8JAkotDsIPwmVawmV6wC9lvbIzuvoZ4zxOu11KC05mtlYS3ywEZF0vZb2JVSuA/Ra2qP2eh1axaSUUsorTRBKKaW8CqUE8ZLdAfiRXkv7EyrXAXot7VG7vI6QaaRWSinlX6FUglBKKeVHmiCUUkp5FRIJQkTOE5FMEckWkQfsjqe5RGSXiGwUkXUikm5tSxSRhSKSZf1MsDvOhkTkVRHJF5FN9bZ5jVs8/mrdow0iMt6+yE/UyLX8SkRyrfuyTkQuqLfvQetaMkXkXHuiPpGI9BGRJSKyVUQ2i8id1vaguy8nuZZgvC+RIvKdiKy3ruXX1vb+IrLSui9vW0sjYC118LZ1LStFJMWWwI0xQf3AM5X4dmAA4ALWA8PtjquZ17ALSGqw7THgAev5A8CjdsfpJe6pwHhgU1NxAxcAn+FZLXAysNLu+H24ll8B93k5drj1/ywC6G/9/3PYfQ1WbD2B8dbzWGCbFW/Q3ZeTXEsw3hcBYqzn4cBK69/7HeAqa/sLwG3W8x8DL1jPrwLetiPuUChBTAKyjTE7jDHVwFvALJtj8odZwOvW89eBS22MxStjzJd41vGor7G4ZwH/MB4rgM4i0rNtIm1aI9fSmFnAW8aYKmPMTiAbz/9D2xlj9hlj1ljPS4GteNZzD7r7cpJraUx7vi/GGFNmvQy3HgY4G3jP2t7wvhy9X+8B00XE21LMARUKCSIZ2FvvdQ4n/0/UHhlggYisFpG51rbuxph94PlFAbrZFl3zNBZ3sN6n262ql1frVfMFxbVY1RLj8Py1GtT3pcG1QBDeFxFxiMg6IB9YiKeEU2SMqbUOqR/vsWux9hcDXdo24tBIEN6yarD13T3dGDMeOB/4iYhMtTugAAjG+/Q8MBAYC+wDnrC2t/trEZEY4D/AXcaYkpMd6mVbe7+WoLwvxpg6Y8xYoDeeks0wb4dZP9vFtYRCgsgB+tR73RvIsymWFjHG5Fk/84EP8PznOXC0qG/9zLcvwmZpLO6gu0/GmAPWL7Ub+Bv/q65o19ciIuF4vlD/aYx539oclPfF27UE6305yhhTBCzF0wbRWUSOzolXP95j12Ltj8f3KlC/CYUEsQpItXoDuPA06MyzOSafiUgnEYk9+hw4B9iE5xrmWIfNAT6yJ8JmayzuecB1Vq+ZyUDx0SqP9qpBXfz38NwX8FzLVVZPk/5AKvBdW8fnjVVP/Qqw1Rjz53q7gu6+NHYtQXpfuopIZ+t5FDADT5vKEmC2dVjD+3L0fs0GFhurxbpN2d26748Hnp4Y2/DU6f3c7niaGfsAPD0v1gObj8aPp77xCyDL+plod6xeYv83niJ+DZ6/eG5sLG48ReZnrXu0EUizO34fruUNK9YNeH5he9Y7/ufWtWQC59sdf724zsBTFbEBWGc9LgjG+3KSawnG+zIaWGvFvAl4yNo+AE8SywbeBSKs7ZHW62xr/wA74tapNpRSSnkVClVMSimlAkAThFJKKa80QSillPJKE4RSSimvNEEopZTyytn0IUqp+kTkaJdRgB5AHVBgva4wxpxmS2BK+Zl2c1WqFUTkV0CZMeZxu2NRyt+0ikkpPxKRMuvnmSKyTETeEZFtIvJHEbnaWhNgo4gMtI7rKiL/EZFV1uN0e69Aqf/RBKFU4IwB7gRGAdcCg40xk4CXgTusY/4CPGmMmQh839qnVLugbRBKBc4qY81rJCLbgQXW9o3AWdbzGcDwelP9x4lIrPGsf6CUrTRBKBU4VfWeu+u9dvO/370w4FRjzJG2DEwpX2gVk1L2WgDcfvSFiIy1MRaljqMJQil7/T8gzVodbQtwq90BKXWUdnNVSinllZYglFJKeaUJQimllFeaIJRSSnmlCUIppZRXmiCUUkp5pQlCKaWUV5oglFJKefX/ATniW5xYrShOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "time_grid = np.linspace(durations_test.min(), durations_test.max(), 100)\n",
    "ev.brier_score(time_grid).plot()\n",
    "plt.ylabel('Brier score')\n",
    "_ = plt.xlabel('Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Negative binomial log-likelihood\n",
    "\n",
    "In a similar manner, we can plot the the [IPCW negative binomial log-likelihood](https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0258%2819990915/30%2918%3A17/18%3C2529%3A%3AAID-SIM274%3E3.0.CO%3B2-5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXzU9ZnA8c8zM8kEkhAISThCAgmE+xIiIIcHpRXoClq1QtV6s62l9twttrvW1dquZ2tbut7Wo57VVqoo4FEqIkdAbkhIwpEAuYBc5E6++8dMYgiTi2Tym+N5v155MfO78vz4wTzzvcUYg1JKqeBlszoApZRS1tJEoJRSQU4TgVJKBTlNBEopFeQ0ESilVJBzWB1AZ8XExJhhw4ZZHYZSSvmVbdu2FRljYj3t87tEMGzYMNLS0qwOQyml/IqIHGltn1YNKaVUkNNEoJRSQU4TgVJKBTlNBEopFeQ0ESilVJDTRKCUUkFOE4FSSgU5vxtHoFRnVNfVc+BEGUdPVVBeXUdZVS11DYb54waSHBthdXhK+QRNBCqgFJRVsSn7FJuzT7Izt5j0vDJq689dc+OhD9K5eGQst8wcxiUjY7HZxIJolfINmgiU3ztWXMmbaTm8t+sEBwvKAYhwOpic0Jfb5yQzMT6K5NgI+vRyEBkWQkV1Ha9uyeEvm49wy5+3Eh0eykXJ/Zk5oj8zh8cwrH9vRDQxqOAh/rZCWWpqqtEpJoJbWVUtmQXlHMwvZ/WeE6zPKMQYmDm8PxePjOWi5P6MG9wHh73tJrCaugbW7svj4wMFbMw8SV5pFQBxkU6mJUUzPSmacfFRpMRFEBkW0hO3ppTXiMg2Y0yqx32aCJS/KCit4oZnN5ORX960bUAfJ99MTeCbqQkkRPc+72sbY8guOsPnWSfZevgUm7NPNSUGgPi+vUiKCWdgVBiDo8JIjo3gikmDsWuVkvITbSUCrRpSfuOpf2WTVXiGn35tJCMHRJIyIJLE6N7d8mEsIgyPjWB4bAQ3zBiKMYbc05Wk55WRnl9GRr6rwXnDwSIKyqpoMPDXbbn8fukFRIeHdsPdKWUdTQTKLxRX1PDKlqMsmjSY5XNTvP77RISE6N4kRPdm3tgBZ+2rq2/gre25/Pc7e7niDxt44oapTBgS1anrl1TUklFQRk1dAzX1DdTXGyYmRBEXGdadt6FUh2giUH7hhY1HqKip5zuXDLc6FBx2G9ddmMiYQX34zkvbuPqJjSwc7+qOmhwbTlSvEIrKqykoreZURQ1Ou43eTge9Q+0cKjrD5uxT7M8rpWWtrE1gRnJ/rpg0mEtHxTIgMuyc3kx19Q3YbaKN2apbaRuB8nkVNXXM+t+PmTq0H8/cdKHV4ZzlZHk19/5jH9sOn+J4SdU5+0PtNmobGpo+9J0OG1OH9mN6Un8mJkQRHuogxC40GMP6jCLe3Xmc7KIzTccmRPdmYJ8wTp2poaCsiqLyGgb2CWNOSgyzU2K4ZGQsfXtr1ZRqnzYWK7/27IZD3P/uPt767kymDu1ndTitqqyp5/DJM5RW1hIb6SQ20kmE01Xorqpt4ExNHZFhDpwOe6vXMMaw93gpX+QUc/TkGY6eqiC/tJr+4aHE9QkjNiKUrMIzbMgsoqSylv7hoaz90cX0j3D21G0qP6WNxcqv1DcYBLDZhJq6Bp75NJvpSdE+nQQAeoXaGTOoT6v7eoW2ngAaiQjj46MYH992m0N9g2FjVhE3P7+VR9dl8OurJpxXzEqBJgJlkWc3HOLAiVIuHzeQOSNjcDrsHC46w0ubjvBGWg6VNfXERDgJd9o5UVLFb76hH3TN2W3CnJRYbpwxlBc+P8z10xMZN7hzDdZKNdJEoHqUMYZH12bwx08yCXXYeHNbLpFOBykDIth+tBiHTZg/fiCJ0b0pKKumsKya1KHRXDLS45rbQe9H80byzo5j3PePfby2bIY2IqvzoolA9RhjDL95/wBP/SubpdMS+OUV49iUfZL3d+ex61gJP5yXwremJRLXR7tQdlRU7xB+8rVR/Nff97B6dx5fnzjI6pCUH/JqY7GIzAceB+zAM8aY//VwzDeBewED7DTGfKuta2pjse/amVPM858dYkCfML7/lZSmhlKAqtp6HnhvPy9tOsJNFw3ll1eM04neukl9g+Hrv/+Usqo61v7oYsKd+v1OncuSXkMiYgcygK8CucBWYKkxZl+zY1KAN4C5xpjTIhJnjClo67qaCHzP5uyT/PGTTD49WESk00F5TR2D+oRx/5XjuWRkLG9uy+XxDw+SV1rFHXOS+PnCMVqF0c0+zzrJ0qc3ER0eyrcvGsqNM4ZqTyJ1FqsSwUXAvcaYy93v7wYwxvym2TEPARnGmGc6el1NBL5jZ04xD605wGeZJ4mJcHLHnCSunzGU9Lwy7n57Fxn55cREhFJUXsOUxL785/zRzEjub3XYAWvr4VM8uT6LD/cX4HTYmBAfRV2Doba+AafDxoLxg1h8wWAdvRykrEoE1wDzjTG3u9/fCEw3xixvdszfcZUaZuGqPrrXGPOBh2stA5YBJCYmTj1y5IhXYlatO15cSWZBORU1dVTU1LN2bz4f7M0jOjyUOy8dzg0zhhIW8mX3yJq6Bp76VxZbDp/mxhlDmTcmTksBPSSzoIznPjvM4aIzhNhthNiFwrJqduaWYLcJl4yMZUJ8FP0jQokOD2VKYj8G9+1lddjKy6waR+Dpf33LrOMAUoBLgSHApyIy3hhTfNZJxjwFPAWuEkH3h6racqKkknmPraeipr5pW3ionR/OS+G22Ukep2gOddh6ZE4gda4RcZEexxVkFpTz1vZcVu04zscHvqyBjeoVwj9/ein9dPK8oOXNRJALJDR7PwQ47uGYTcaYWuCQiKTjSgxbvRiX6qQ/fpxJbX0Dz99yIXGRTnqHOoiLdGqjpJ8ZERfBz+aP5mfzR1Nb38Dpihoy8sr59nObefyjg9y7aJzVISqLeHPx+q1AiogkiUgosARY1eKYvwOXAYhIDDASyPZiTKqTck5V8PrWHJZcmMhlo+IYNziKpJhwTQJ+LsRuIy4yjNkpMSyZlsjLm46QVVje/okqIHktERhj6oDlwBpgP/CGMWaviNwnIovch60BTorIPuAT4D+MMSe9FZPqvN99eBC7TVg+d4TVoSgv+dG8kYSF2PnN6gNWh6Is4tWvdcaY1cDqFtvuafbaAD92/ygfk1lQzt++yOXWWUkM0EFeASs20smdlw3noQ/S2ZhVxMzhMVaHpHqYN6uGlJ/73YcZhIXY+c6l1q8BoLzr1llJxPftxa/e3U99g/bHCDaaCJRH+46X8u6uE9wyaxgxOjAp4IWF2PnZgtHsO1HKsxu0mS7YaCJQ5zDGcO+qvfTrHcKyOVoaCBZXTBzE/HEDeeiDdLYfPW11OKoHaSJQ53h7+zG2HD7FigWjiep97hgBFZhEhAevmcigvmF8/5UvKK6osTok1UM0EaizlFTW8pv39zMlsS/XTk1o/wQVUKJ6hfDHpVMoKKvip2/uwhhDVW09+46X8vGBfNZnFLIxs4gvjp7WtoQAop3B1VkeXZvOqTM1vHDrNJ0dNEhNSujL3QvGcN+7+5j94CfklVZ5/NCfOzqOP10/5aypRZR/0kSgmuzOLXFPEz1MV7sKcrfMGkbu6UqOFVdw9ZR4UgZEMqRfLxoM1NY38MVR14SDNz+/hWduuvCsKceV/9GnpwCoq2/gF3/fTUyEkx9/baTV4SiLiQj3XDG21f0zkvszKCqMn7y5k+uf3sSfb5mmcxX5MW0jUAA8/9lhduWW8MsrxtLHwyRySrV05QXxPHHDVPbnlXHHi2l4c5Er5V2aCBRHTp7h0XXpzBszgK9P0KUOVcd9dewA7l88jrQjp3l/T57V4ajzpIkgyBljuPvt3YTYbPzqyvG6ZoDqtGumJjBqQCQPfnCAmroGq8NR50ETQZB7Iy2HjVknWbFwNAOjdD4h1Xl2m7Bi4WiOnKzglc26aJQ/0kQQxIrKq/nVe/uZlhTN0gsTrQ5H+bFLR8Yyc3h/fv9xJqVVtVaHozpJE0EQ++PHmVTU1PPrq8brmAHVJSLC3QvGcOpMDU+uzwJc1Y7VdfXaiOwHtPtokDp6soK/bD7CN1MTGBEXaXU4KgBMGBLF4smDeWJ9Ni99foSKmnrqGgyRYQ5S4iJIiYvkwqRorrogHrt+8fApmgiC1KPr0rHbhB/O03WFVff5+cIxhDsdhNiEcKeDXiF28suqOJhfzrr9+byelsMLGw/zwFXjmTikLwD1DYZ9x0sZEOUkLlLbqaygiSAI7TlWwjs7jvO9y4brgjOqWw3oE8avr5rgcZ8xhn/sOsH97+5j8crPuHrKEMqqavk86ySlVXWkxEWw+gdzCLFrjXVP07/xIPTgBwfo2zuEf79Ep5hWPUdEWDRpMB/95BJuumgYb2/PZc+xUhaMH8T3547gYEE5L35+bq+jqtp6C6INLloiCDKfHizk04NF/NfXx+gIYmWJPmEh3LtoHHcvHE2o3YaIYIxhZ24Jv1uXwaJJg4mNdC2GtDGriNv+nMbPF47mxouGWRt4ANMSQRApKq/mp2/uZFj/3twwY6jV4agg53TYmwYwigi/vGIsVXX1PPTBAQC2Hz3N7S+kUVlbzytbcqwMNeBpIggSdfUN7sVGavnT9VN16mDlc4bHRnDr7CTe3JbLq1uOcvNzW4iNdPK9y4az/0QpGfllVocYsLyaCERkvoiki0imiKzwsP9mESkUkR3un9u9GU8we2xdBp9nn+SBqyYwdnAfq8NRyqPvz00hLtLJ3W/vJtzp4C+3T+fmmUnYBFbtOG51eF3S0GA4UVJpdRgeeS0RiIgdWAksAMYCS0XE07y2rxtjJrt/nvFWPMFs7d48/vTPLJZOS+SaqUOsDkepVkU4HTxw1QQmxEfx8u3TGdKvN7GRTmaNiOGdncf8dnDa51knWbRyAzP/92OyCsutDucc3iwRTAMyjTHZxpga4DVgsRd/n/LgYH4ZP3ljJxPio/hlG/PLK+Urvjp2AP/4/myGx0Y0bVs8OZ6cU5V8kVN81rGnz/j2uspZheXc8WIaS5/exLHTlRgDe4+XWh3WObyZCOKB5i08ue5tLV0tIrtE5K8i4nGRXBFZJiJpIpJWWFjojVgD0qkzNdz2QhrOEDtP3KjtAsp/XT5uAKEO21nVQys/yeSC+9fxP//YS229b816ml1Yzo9e38FXH1vPxswi/uPyUfzzp5chAlkFwVUi8DSGvGW57h/AMGPMROBD4AVPFzLGPGWMSTXGpMbGxnZzmIGppq6B77y8jbzSKp769lTi+/ayOiSlzltkWAjzxsTx7q4T1NU38Mrmozy8Jp2RAyJ4/rPDXP/0ZgrKqqwOk/oG17Tu8x5bz/t7TnD7nGTW/+dlfO+yEUT1DmFIv15kF52xOsxzeDMR5ALNv+EPAc5q7THGnDTGVLvfPg1M9WI8QcMYw3//fQ9bDp3i4WsmMiWxn9UhKdVliybFU1RezX3v7uO//r6bS0fF8t5dc3h8yWR2HSvmij9sIO3wKUtj3JBZxKtbjnLdhYl8+p9z+fnCMcREOJv2D4+NCLoSwVYgRUSSRCQUWAKsan6AiDRfDmsRsN+L8QSNZzcc4vW0HL4/dwSLJ3uqjVPK/1w6KpbIMAcvfn6EKYn9+L/rpxJit7F4cjx/u3MWToed657axOMfHqTOy1VFa/bmse3I6XO2v5GWQ9/eIdy7aGzToLjmkmMiOFR0hoYG32r09loiMMbUAcuBNbg+4N8wxuwVkftEZJH7sLtEZK+I7ATuAm72VjzB4sN9+Tywej8LJwzkR/N0EXoVOMJC7Nx00TAuHNaPZ2+6kF6hX7Z5jRnUh/fums2iSYP57YcZrsbZYu901cw5VcHyV7Zz16tfnLUiW3FFDev25nPl5HicDs/tccmx4VTW1pNXan01VnNeHUdgjFltjBlpjBlujHnAve0eY8wq9+u7jTHjjDGTjDGXGWMOeDOeQLfveCl3vfYFE+KjePTaybrGgAo4P718FG9+ZyZRvc+dHiUyLITfXjeZ3143if0nyvi333/KYS/Uxz/+0UHqGwzHiit5a3tu0/Z3dhynpr6Ba1Nb76Ld2BOqq11Ic05V8OAHB5j94Mf87Yvc9k9oh44sDhAFZVXc/sJW+oSF8PS3U8/6tqRUMLnqgiG8s3wWALe+sJWSiu5bMe1gfhlvb8/l1llJTE7oyx8/zmwqFbyRlsO4wX0YNziq1fOHx4YDkF14fgkqs6CMm57bwsUPf8KT67Moq6rjDx9ldrmqSROBnztcdIb7/rGPrzy6ntMVtTxzU6pOLa2C3vDYCJ68MZWcUxXc+cq2bute+ti6DHqF2LnzshH8YF5KU6lg7/ES9h4v5ZupHnvAN4mNdBLhdJB9HiWCz7NOctWfNrL7WAl3zU3hsxVzuW/xOLKLzvCvg13rVq+JwE8VlFVxx4tpXPrIP3nx88NcNiqOv373IsbHt/5tRKlgMi0pmt98YyKfZZ7knnf2dnlU8q7cYt7fk8ftc5KJDg/l0pGxTHKXCl7dcpRQu43Fkwe3eQ0RYXhsOFmdLBG8s+MY335uMwP6hLFq+Sx+9NWRDIrqxYLxg4iLdPLnjYe7cGc6DbVf2nbkNHf+ZRsllbX84CspXD89kTgtBSh1jmumDiGrsJz/+2cW8X3DWD73/Ffke3hNOv16h3D7nCTA9aH+w3kp3PL8Vl7edJR/mziIvr1D271OcmwEm7NPduh3NjQYVn6SyaPrMpiRHM2TN6Se1T4S6rBxw4yhPLYug+zCcpKbjcbuDC0R+BFjDC9tOsKSpz7H6bDztztd3ww0CSjVuv/42iiuuiCeR9Zm8Myn2ed1jfUZrnU87rx0BJHN1vFoLBUAXNtOtVCj4bHhHC+poqKmrs3jiitquP3FNB5dl8FVF8Tzwq3TPDaSL52WSKjd5nFRn47SEoGfqG8w/HLVHl7edJTLRsXyu+su8PiPQil1NptNePiaiVTX1fOr9/bjdNg6tchNRU0dv/jbbpJjw7nxorPX8RAR7ls0jje35TB7REyHrtf4rT278EyrVbm7cov57svbKSir4r7F47hxxtCmtRtaio108m+TBvFmWg4//trI81pwShOBH6iqreeuV79g7b58/v3iZH42f7R2DVWqExx2G48vuYCaum389zt7CXc6+MaUjs3E+7sPD5J7upLXl83wOF/XpIS+TaWCjkhu7DlU5DkRZOSXcc0TnxMb4eTN78xkcgeufcvMJN7efoyXPj/SNMNwmMPe4S+Lmgh8XHGFa+K47UdPc+8VY7l5VpLVISnll0LsNlZeP4Ubn93Cfe/uY/74gfQObfsjcM+xEp75NJul0xKYnty/W+IY1j+8zcnnHl2bTqjdxt+/N8vj6GRPJgyJYurQfjy8Jp2H16QDYBNY+6OLGREX2e752kbgwxoaDLf8eSu7j5Ww8ltTNAko1UVOh52fzR9FcUUtb2xte/nLuvoG7n57N/0jnKxYMKbbYggLsbc6+dzOnGLW7M3njjnJHU4CjX533WR+fdUEfn3VBP79kmQaDOSVVLd/IpoIfNo/dh3ni6PFPHDleBZOGNT+CUqpdk0dGk3q0H48s+FQm3MSPbE+i93HSrj3inFE9ere9rjkmAiPYwkeWZtOdHgot83p/Je+hOjefGt6It+ansjX3Z8X1XX1HTpXE4GPqq6r55G16YwZ1IerO1iXqZTqmGUXJ5N7upLVe/I87n8jLYdH1mZwxaTBLJwwsNt///DYCLILz558bmNWkbtn0nAinF2rtW+c66i6rmMD6TQR+KiXNx0l51Qldy/QhmGlutu8MQNIjg3nyfVZ5ww0+2BPHive2sWclBgeuXZiq711uqLl5HPGGB5Zk87APmHcMGNoO2e3z+lwfbRricCPlVTW8oePDzInJYaLR+pCPEp1N5tNWDYnmb3HS9mY9eXgrg0Hi7jr1S+YnNCXJ2+c2uosol3V2HMoq7Ccqtp6nt1wiO1Hi7nrKyndspKgM8SdCGo7ViLQXkM+6In1WZRU1vKz+aOtDkWpgHXlBfE8ui6DR9em86+MQtZnFHIgr4zRAyN5/uZp7fYo6ooR7rEErhHBZyiprOWCxL5tzlzaGaH2xhKBJgK/tOXQKZ7bcIgrJ8frvEFKeVFYiJ1bZyXx4AcH2H2shAuHRbNiwWiuS03w+mDN2EgnMRFO9hwrYf74QSydlsCMpP7dVg3sdJcqajQR+BdjDC9+foT7391HQnRvLQ0o1QPumJPE1KH9GDe4D+FdbKDtDBHhneWzCHPY6B/RuW6iHdHZNgJNBBYzxlBYXs2D76fz1vZc5o2J47HrJp/XMHGlVOc47DamJUVb8rvj+/by2rUdNsEmWjXk837xt91szDrJseLKpuLbD+elcNfcFO0lpJTqEhHB6bBrIvBlmQVl/GXzUaYNi+arYwcQ37cXkxL6dmhOEaWU6ohQh43qWq0a8lnv73YNYvn90gsYGKVTSCulup/TYdMBZb5s9Z48pg7tp0lAKeU1zhBbh3sNaSLoYYeKzrD/RCkLxnf/sHWllGrUmTYCryYCEZkvIukikikiK9o47hoRMSKS6s14fMH7e04AsEAnkVNKeZGrasjiKSZExA6sBBYAY4GlIjLWw3GRwF3AZm/F4kve353HpIS+Xu06ppRSvtJGMA3INMZkG2NqgNeAxR6Oux94CKjyYiw+IedUBbuPlbBQq4WUUl7m6jVkfSKIB5qv/JDr3tZERC4AEowx77Z1IRFZJiJpIpJWWFjY/ZH2kKZqofFaLaSU8i6nw051G+stNOfNROBpVFTTfK8iYgN+C/ykvQsZY54yxqQaY1JjY/13Ns7Vu/MYH9+HxP69rQ5FKRXgnJ0YR+DNRJALJDR7PwQ43ux9JDAe+KeIHAZmAKsCtcH4WHElO3KKtTSglOoRzhC7T3Qf3QqkiEiSiIQCS4BVjTuNMSXGmBhjzDBjzDBgE7DIGJPmxZgs8/cvjgHokpNKqR7hE43Fxpg6YDmwBtgPvGGM2Ssi94nIIm/9Xl9U32B4ZfNRLkruT1JMuNXhKKWCQGe6j3p1igljzGpgdYtt97Ry7KXejMVK/8oo5FhxJXcv1KmllVI9w1d6DSm3lzcdISbCydfGardRpVTP8JVeQwrIPV3Bx+kFLLkwgVCH/nUrpXqG0+Gaa8gY0+6x+snkZa9tcQ2lWDItoZ0jlVKq+zQtYN+BBmNNBF5UU9fAa1tzmDsqjiH9dOyAUqrnOB2udYs1EVhs3b58isqruWHGUKtDUUoFmdBOrFusicBLSipqeWxdOvF9e3HxSP8dDa2U8k9NC9h3oOfQeScCEbn6fM8NdFW19dzxUhpHT1Xw8LUTsesaxEqpHtaYCGo60HOoKyWC33bh3IDV0GD4yZs72XLoFI9+czIzh8dYHZJSKgg1tRF4s0SA50nlgpoxhgdW7+e9XSf4+cLRLJo02OqQlFJB6steQ+23EXRlZHH7nVODSEVNHXe/vZt3dhzn5pnDuGNOstUhKaWCWFMbQQd6DbWZCERkN54/8AXQYbJuR06e4d9f2kZ6fhn/cfkovnvJcES0wKSUsk63JQLg37ohnoC2K7eYG57ZjIjw51umcYn2EFJK+YAv2wi6WDVkjDnS2j4R+QyY1cnYAs7Tnx7CbhNWLZ9NQrQOGlNK+Yae6jWU2IVzA0JVbT0f789n/viBmgSUUj6lp3oNBX1j8casIs7U1HP5OG0uUUr5ls7MNdReY/E3WtsF9OpsYIHmgz15RIY5dKyAUsrnODsxxUR7jcVXtLHv3Y6HFHjq6htYty+feWMG6PTSSimfE9pdvYaMMbd0T0iBZ8uhU5yuqNVqIaWUTwq1uxuLu6H7KCJyCXDaGLNLRL4JXAxkAX8yxlR3LVT/9cHePMJCbNpdVCnlkxx2Gw6bdL1qSERWAhOBMBFJByKAD4CZwHPA9V0P1/80NBjW7M3j0pFx9Aq1Wx2OUkp55OzgusXtlQguM8aMFZEw4BgQZ4ypF5EngV3dEKdf2pFbTH5pNfPHa7WQUsp3OUPs3bIwTRWAMaYKOGKMqXe/N0BtexcXkfkiki4imSKywsP+74jIbhHZISIbRGRsuxH7gA/25BFiF+aOibM6FKWUalWo3dYtvYbiROTHuLqLNr7G/b7NynERsQMrga8CucBWEVlljNnX7LBXjDFPuI9fBDwGzG83agsdK67krW25zBoRQ5+wEKvDUUqpVjlDbN0y19DTQKSH1wDPtHPuNCDTGJMNICKvAYuBpkRgjCltdnw4Pj5IrbSqlluf30pNXQN3LxhjdThKKdUmp8PW9V5Dxpj/6UIM8UBOs/e5wPSWB4nI94AfA6HAXE8XEpFlwDKAxERrZraorW/gzpe3k1VYzp9vmcaogZHtn6SUUhZyOjrWRtBer6F72thtjDH3t3W6p3M8XGQlsFJEvgX8F3CTh2OeAp4CSE1N7fFSgzGGn7+9mw2ZRTx8zURmp+hIYqWU73M6OtZG0F5j8RkPPwC3AT9r59xcIKHZ+yHA8TaOfw24sp1rWuK1rTm8uS2Xu+aO4NrUhPZPUEopH+AM6Ybuo8aYRxtfi0gk8APgFlwf2o+2dp7bViBFRJJwdT1dAnyr+QEikmKMOeh++3XgID7mWHElD7y3n4uS+/PDeSOtDkcppTos1G6jtLKu3eM6MrI4Glcd/vXAC8AUY8zp9s4zxtSJyHJgDWAHnjPG7BWR+4A0Y8wqYLmIzMPVFfU0HqqFrGSMYcVbuzDG8NA1E7HZdNUxpZT/cLURdH1k8cPAN3DVz08wxpR3JghjzGpgdYtt9zR7/YPOXK+nvbolh08PFvGrK8fregNKKb/jDOlYr6H22gh+AgzG1Yh7XERK3T9lIlLazrl+Lfd0BQ+8t49ZI/pz/fSgX4NHKeWHXI3FXW8jCNr5le9/1zXc4cGrJ+pC9Eopv9TR7qNB+0Hflj3HSlizN59lFw9nSD+tElJK+SfXpHNd7z4alB7/6CB9whzcMnuY1aEopdR5C+1g1ZAmghb2HCth3b58bp+TrHMJKaX8mtNhp67BUN/Q9jhcTQQt/O5DV2ng5lnDrPPegpwAAA7ASURBVA5FKaW6pHEB+/Z6DmkiaGZ3bgkf7s/nDi0NKKUCQEcXsNdE0MzjH2UQ1SuEm7Q0oJQKAE6HawXF9toJNBG4lVXV8uH+Aq6fnqilAaVUQGgqEbQz35AmArfMAteg6ckJfS2ORCmlukeoVg11zkF3Ihg5QNcZUEoFhi/bCLRE0CEH88twOmw6p5BSKmA4Q7SNoFMOFpQzPDYCu84wqpQKENprqJMO5peTMiDC6jCUUqrbaNVQJ5RX13GsuJKUOE0ESqnAEaq9hjouy91QnKINxUqpAPLlOAKtGmpXRn4ZgJYIlFIBpbFqSKeY6IDMgnJCHTYStceQUiqANM41pG0EHXCwoJzkmHAcdv3rUEoFDp1iohMy8st0IJlSKuBo99EOqqipI/e09hhSSgWeULsP9BoSkfkiki4imSKywsP+H4vIPhHZJSIfichQb8bjSWZTjyFNBEqpwGKzCaH29lcp81oiEBE7sBJYAIwFlorI2BaHfQGkGmMmAn8FHvJWPK05mK9dR5VSgcvpsFnaa2gakGmMyTbG1ACvAYubH2CM+cQYU+F+uwkY4sV4PDpYUE6o3cZQ7TGklApAzhCbpW0E8UBOs/e57m2tuQ1434vxeHQwv4zkWO0xpJQKTE6Hvd2qIYcXf7+n2ds8rqAsIjcAqcAlrexfBiwDSExM7K74AFeJYOKQqG69plJK+Qqnw8I2AlwlgIRm74cAx1seJCLzgF8Ai4wx1Z4uZIx5yhiTaoxJjY2N7bYAK2vqyTldQUqctg8opQJTqMNGda11VUNbgRQRSRKRUGAJsKr5ASJyAfAkriRQ4MVYPMoqLMcYGKk9hpRSAcrpsFFTb1GJwBhTBywH1gD7gTeMMXtF5D4RWeQ+7GEgAnhTRHaIyKpWLucVB/LccwxpIlBKBSinw97uOAJvthFgjFkNrG6x7Z5mr+d58/e3Z1P2Sfr2DiE5RhOBUiowOUNsnKmua/OYoO0qY4xhY2YRM4f3x6arkimlApTVjcU+7fDJCo6XVDFzeIzVoSillNeEaiJo3YbMIgBmj9BEoJQKXK5xBDrpnEcbM4uI79uLof11RLFSKnBZPcWEz6pvMGzMOsnM4f0R0fYBpVTg0jaCVuw7XkpJZS2zU7RaSCkV2Jwh7XcfDcpE0Ng+cNHw/hZHopRS3uV02KjSNoJzbcwqYuSACOIiw6wORSmlvCrUbsN4nOXtS0GXCKpq69ly6BSztLeQUioINC5g35agSwTbj56muq6BWTp+QCkVBBoXsG9L0CWCjZknsduE6cnRVoeilFJe17iAfVuCLhFsyj7JhPgoIsNCrA5FKaW8TquGWjDGcCCvTBeiUUoFDa0aaiH3dCXl1XWMGqgL0SilgkNoB5bhDapEkO5ef2C0JgKlVJDQqqEWDuSVAjBygCYCpVRw0KqhFg7klTGkXy9tKFZKBQ3tNdRCel6ZVgsppYKKVg01U11XT3bRGUYP7GN1KEop1WO0sbiZzIJy6huM9hhSSgUVZ4i2ETTRHkNKqWCkbQTNpOeVEWq3MSwm3OpQlFKqx1ieCERkvoiki0imiKzwsP9iEdkuInUico03Y9mfV8aIuAhCOlBfppRSgcLS7qMiYgdWAguAscBSERnb4rCjwM3AK96Ko1F6XqlWCymlgk6IXWhvRV6HF3//NCDTGJMNICKvAYuBfY0HGGMOu/e1vY5aFxVX1JBfWq0NxUqpoCMi7fYc8mY9STyQ0+x9rntbp4nIMhFJE5G0wsLCTp9/oLGheJB2HVVKBZ/22gm8mQg8FUbaWTDNM2PMU8aYVGNMamxsbKfPP3DCNbWEVg0ppYJRe11IvZkIcoGEZu+HAMe9+PtalZ5fRt/eIcRFOq349UopZSkrSwRbgRQRSRKRUGAJsMqLv69VB9xTS0h7LSZKKRWA3v/BnDb3ey0RGGPqgOXAGmA/8IYxZq+I3CciiwBE5EIRyQWuBZ4Ukb3dHUdDg3HPMaTtA0qp4NTeRJve7DWEMWY1sLrFtnuavd6Kq8rIa/bnlVJRU6/tA0op1YqAH131h48yiXA6uHzcQKtDUUopnxTQiWBXbjEf7M3j9jlJ9AsPtTocpZTySQGdCB5Zm0G/3iHcNjvJ6lCUUspnBWwi2JR9kn9lFHLnpSN0RTKllGpDQCYCYwyPrElnQB8nN1401OpwlFLKpwVcIjDG8GZaLmlHTvP9uSmEdWBRBqWUCmZe7T7a03bmFPPA6v1sOXSKCfFRfDM1of2TlFIqyAVMIvjVu/t4ZsMh+oeH8qsrx7PkwgQcuvaAUkq1KyASwZnqOp777BALJwzkwasnauOwUkp1QkB8Zd5zrIQGA9dOTdAkoJRSnRQQiWBHTjEAE4dEWRyJUkr5n4BIBDtzi0mM7k3/CJ1mWimlOisgEsGOo8VMSuhrdRhKKeWX/D4RFJRWcbykismaCJRS6rz4fSJobB+YnKDtA0opdT78PhHszC3GYRPGDdZEoJRS58PvE8GOnGJGD4rUqSSUUuo8+XUiaGgw7MopYdIQbR9QSqnz5deJILuonLLqOm0oVkqpLvDrRLAjpwRAE4FSSnWBnyeC00Q4HQyPjbA6FKWU8lteTQQiMl9E0kUkU0RWeNjvFJHX3fs3i8iwzlx/Z04JE4dEYbNJd4WslFJBx2uJQETswEpgATAWWCoiY1scdhtw2hgzAvgt8GBHr19VW8/+E6VaLaSUUl3kzWmopwGZxphsABF5DVgM7Gt2zGLgXvfrvwJ/FBExxpjWLpqRX8a8x9ZTW99AXYPRqSWUUqqLvJkI4oGcZu9zgemtHWOMqROREqA/UNTaRcNC7IwaEAnA9KRoZo+I6c6YlVIq6HgzEXiquG/5Tb8jxyAiy4BlAImJiay8fkrXo1NKKQV4t7E4F2i+aPAQ4Hhrx4iIA4gCTrW8kDHmKWNMqjEmNTY21kvhKqVUcPJmItgKpIhIkoiEAkuAVS2OWQXc5H59DfBxW+0DSimlup/Xqobcdf7LgTWAHXjOGLNXRO4D0owxq4BngZdEJBNXSWCJt+JRSinlmVcXrzfGrAZWt9h2T7PXVcC13oxBKaVU2/x6ZLFSSqmu00SglFJBThOBUkoFOU0ESikV5MTfemuKSBmQbnUc3SSGNkZR+5FAuQ/Qe/FFgXIfYO29DDXGeByI5dVeQ16SboxJtTqI7iAiaYFwL4FyH6D34osC5T7Ad+9Fq4aUUirIaSJQSqkg54+J4CmrA+hGgXIvgXIfoPfiiwLlPsBH78XvGouVUkp1L38sESillOpGmgiUUirI+VUiEJH5IpLuXux+hdXxdIaIHBaR3SKyQ0TS3NuiRWSdiBx0/9nP6jg9EZHnRKRARPY02+YxdnH5vfsZ7RIRn1pFqJV7uVdEjrmfzQ4RWdhs393ue0kXkcutifpcIpIgIp+IyH4R2SsiP3Bv97vn0sa9+NVzEZEwEdkiIjvd9/E/7u1JIrLZ/Uxed0/Lj4g43e8z3fuHWRa8McYvfnBNZZ0FJAOhwE5grNVxdSL+w0BMi20PASvcr1cAD1odZyuxXwxMAfa0FzuwEHgf1+pzM4DNVsffgXu5F/iph2PHuv+dOYEk978/u9X34I5tEDDF/ToSyHDH63fPpY178avn4v67jXC/DgE2u/+u3wCWuLc/AXzX/fpO4An36yXA61bF7k8lgmlApjEm2xhTA7wGLLY4pq5aDLzgfv0CcKWFsbTKGPMvzl05rrXYFwMvGpdNQF8RGdQzkbavlXtpzWLgNWNMtTHmEJCJ69+h5YwxJ4wx292vy4D9uNYA97vn0sa9tMYnn4v777bc/TbE/WOAucBf3dtbPpPGZ/VX4Csi4mn5Xq/zp0TQtNC9Wy5t/2PxNQZYKyLb3GswAwwwxpwA138GIM6y6Dqvtdj99Tktd1eZPNesis4v7sVdpXABrm+gfv1cWtwL+NlzERG7iOwACoB1uEorxcaYOvchzWNtug/3/hKgf89G7OJPiaBDC937sFnGmCnAAuB7InKx1QF5iT8+p/8DhgOTgRPAo+7tPn8vIhIBvAX80BhT2tahHrb5+r343XMxxtQbYybjWqN9GjDG02HuP33mPvwpETQtdO82BDhuUSydZow57v6zAPgbrn8k+Y3Fc/efBdZF2Gmtxe53z8kYk+/+D9wAPM2X1Qw+fS8iEoLrg/Mvxpi33Zv98rl4uhd/fS4Axphi4J+42gj6ikjjvG7NY226D/f+KDpebdmt/CkRbAVS3C3wobgaV1ZZHFOHiEi4iEQ2vga+BuzBFf9N7sNuAt6xJsLz0lrsq4Bvu3upzABKGqsqfFWLuvKrcD0bcN3LEnfvjiQgBdjS0/F54q5LfhbYb4x5rNkuv3surd2Lvz0XEYkVkb7u172AebjaOz4BrnEf1vKZND6ra4CPjbvluMdZ3dLemR9cPR8ycNW7/cLqeDoRdzKuXg47gb2NseOqD/wIOOj+M9rqWFuJ/1VcRfNaXN9ibmstdlzF3ZXuZ7QbSLU6/g7cy0vuWHfh+s85qNnxv3DfSzqwwOr4m8U1G1c1wi5gh/tnoT8+lzbuxa+eCzAR+MId7x7gHvf2ZFyJKhN4E3C6t4e532e69ydbFbtOMaGUUkHOn6qGlFJKeYEmAqWUCnKaCJRSKshpIlBKqSCniUAppYKcPy5er1SPEJHGrpgAA4F6oND9vsIYM9OSwJTqZtp9VKkOEJF7gXJjzCNWx6JUd9OqIaXOg4iUu/+8VETWi8gbIpIhIv8rIte756XfLSLD3cfFishbIrLV/TPL2jtQ6kuaCJTquknAD4AJwI3ASGPMNOAZ4PvuYx4HfmuMuRC42r1PKZ+gbQRKdd1W4563R0SygLXu7buBy9yv5wFjm00330dEIo1r/n2lLKWJQKmuq272uqHZ+wa+/D9mAy4yxlT2ZGBKdYRWDSnVM9YCyxvfiMhkC2NR6iyaCJTqGXcBqe7VtvYB37E6IKUaafdRpZQKcloiUEqpIKeJQCmlgpwmAqWUCnKaCJRSKshpIlBKqSCniUAppYKcJgKllApy/w+UpHDFkEtimAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ev.nbll(time_grid).plot()\n",
    "plt.ylabel('NBLL')\n",
    "_ = plt.xlabel('Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Integrated scores\n",
    "\n",
    "The two time-dependent scores above can be integrated over time to produce a single score [Graf et al. 1999](https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0258%2819990915/30%2918%3A17/18%3C2529%3A%3AAID-SIM274%3E3.0.CO%3B2-5). In practice this is done by numerical integration over a defined `time_grid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16816426632487566"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev.integrated_brier_score(time_grid) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4969224495850561"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev.integrated_nbll(time_grid) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}