{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PC-Hazard Example \n",
    "\n",
    "In this notebook, we will present a simple example of the `PCHazard` 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 PCHazard\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 `PCHazard` is a continuous-time method, but requires defined intervals in which the hazard is constant. Hence, we need to perform an operation similar to a discretization of the time scale.\n",
    "We let `num_durations` define the size of this (equidistant) interval grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_durations = 10\n",
    "labtrans = PCHazard.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.LabTransPCHazard"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(labtrans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `y_train` now consist of three labels: the interval index, the event indicator, and the proportion of the interval before the event/censoring occur (i.e, $\\rho(t_i)$ in the [paper](https://arxiv.org/abs/1910.06724))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2, 2, 6, ..., 1, 5, 3]),\n",
       " array([0., 1., 0., ..., 1., 0., 0.], dtype=float32),\n",
       " array([0.7965465 , 0.69519496, 0.7370492 , ..., 0.06606603, 0.5865238 ,\n",
       "        0.96302533], dtype=float32))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_train"
   ]
  },
  {
   "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 `num_nodes` output nodes.\n",
    "We also have batch normalization and dropout between the layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PCHazard(net, tt.optim.Adam, duration_index=labtrans.cuts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEKCAYAAADw2zkCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV9dnH8c+Vk71DEgIkIQHCnrKHIi5EwFFF3IjaUh+txW5r26fjaeuo1daqpY6qbV1VVHDgQqaI7B02gYSVkED2CSfJ7/njnETAJGfknNwZ1/v1Oi9yzrnHlxDOld+4f7cYY1BKKdWxBVkdQCmllPW0GCillNJioJRSSouBUkoptBgopZRCi4FSSikg2OoAvkhKSjKZmZlWx1BKqTZl/fr1J4wxyQ291yaLQWZmJuvWrbM6hlJKtSkicrCx97SbSCmllBYDpZRSWgyUUkrRRscMlFLtj8PhIC8vD7vdbnWUNi88PJy0tDRCQkI83keLgVKqVcjLyyMmJobMzExExOo4bZYxhsLCQvLy8ujRo4fH+2k3kVKqVbDb7SQmJmohaCYRITEx0esWlhaDZjLGcKiwAl0KXKnm00LgHw19H1fuOdHkPgEtBiKSLiJLRCRbRLaLyNwGtokTkfdEZLNrmzsCmcmf7I4afvTfzUz80xIeWrRTC4JSqtX63wXbmnw/0C2DauBHxpj+wFjgXhEZcM429wI7jDFDgUnAn0UkNMC5mu14iZ0bnl3N2xsPMyozgWeX7+d37+8IWEHYdriY7UeKA3Jsd2pqDe9tPsKMv6/i1a8OWZJBqUA7deoUzzzzjNf7TZ06lVOnTnm93+zZs3nrrbe83s9XdkdNk+8HdADZGHMUOOr6ulREsoFUYMeZmwEx4mzXRANFOIuIXxVXOvjTxzu5sE9nLhuQ0qxjbco9xZx/raOsqpp5t47g8oEp/O79Hbz4RQ41tYbfXjWQgrIqlu0qYOVeZ9MsJTaclNhwuneK5KK+yQTbzq7D1TW1LN9TQJ+UGNISIutfN8bwwsoDPLRoJ7XGcOOodH5yeT86RflWL40xFJafxlFTS3WNobrW0C0+nLBg2ze2rak1fLD1KE8u3sPe/DJiwoN58J2tFJVXce9FWX5p0h8vsfPW+jxqaw2RYcFEhdrokRTFmJ6JzTpu5ekaNhw6ybieiQQFadeDcq+uGNxzzz1nvV5TU4PN9s3/H3U+/PDDQEfzi6rq2ibfb7HZRCKSCZwHfHXOW08BC4EjQAxwgzHmG6lFZA4wB6B79+5enftUxWlue2ENWw8X85/Vh7hsQAq/uWogqfERXv89th0u5ubnVpMYHcq/7hpPvy6xAPzv9AGE2IJ4dvl+luzKJ7eoEoDkmDDCQ4I4XlLFadc/RlbnaB6c2o+L+nZGRFi+u4D/e38He/LLCLEJN43uzr0XZREdFszP5m/h/S1HuXxgCmkJkby0KodF247x48l9mTKoC4lRoR5/KG/NK+YX725lS97ZLYzIUBsTspK4qG9nhqbHsSWvmJV7T/DlvkKKyk/TJyWap24+j8kDuvDA/C089sluTlY4+MXU/o1+0JbaHdiChIgQW4P5iisc/H3ZPl784kCDP6RTBnbhN1cNpEtcuEd/tzMVllVx58vr2Jx7iol9kvnz9UNJjgnz+jjKOr99bzs7jpT49ZgDusXy6ysHNvr+Aw88wL59+xg2bBghISFER0fTtWtXNm3axI4dO7jmmmvIzc3Fbrczd+5c5syZA3y9PE5ZWRlXXHEF559/PqtWrSI1NZUFCxYQEeH+c2bx4sX8+Mc/prq6mlGjRvH3v/+dsLAwHnjgARYuXEhwcDCTJ0/mscce48033+S3v/0tNpuNuLg4li9f7tHf310xkJbo5xaRaGAZ8AdjzNvnvDcDmAD8EOgFfAoMNcY0+pMwcuRI4+naRIVlVdz6whr2FZTx1E3nceBEOX/5bA8A917UiysGd6VnUpRHH6i5RRVc+/dVhNqCePue8aTEnv1BZYzh78v2sWL3CSZkJTKpb2cGdotFRDDGUFzp4Mt9hTz68S4OnChnfK9EIkJsLN6ZT0ZiJD+4tA9rcor479pcbEFC59gwDp+s5CeX9+PuC3siIuw6VsqvFmxjzYEiAGLDg+mZHE3flBgm9U3mgj7JRIedXeNL7Q7+/Mlu/vVlDonRYXz7/B7ERoQQ7Pog35x3iiU7Czh8qrJ+n5TYMCb0SmLywBQmD+hS/6FfW2v43fs7eGlVDlcP68Z3J/aif9eY+u/flrxT/GP5fhZtPUqtgdDgIOIjQoiPDCEuwvmICgtmyc58Squq+dawVH5wWR+6xIVTcbqG8qpq3t10mL9+tocQWxA/ntyH28ZlYjun6Kzae4LfvLed/NIq7p2UxW3jMggPsXGwsJzb/7mGo8V2bh2bwX9WHyQmPITHZw5lYp8G1+dSrUR2djb9+/cHrCkGOTk5TJ8+nW3btrF06VKmTZvGtm3b6qdnFhUV0alTJyorKxk1ahTLli0jMTHxrGKQlZXFunXrGDZsGDNnzuSqq67i1ltvbfB8s2fPZvr06UyfPp3evXuzePFi+vTpw6xZsxg+fDizZs1i3Lhx7Ny5ExHh1KlTxMfHM3jwYD766CNSU1PrX2vImd9PgN6/+JC9f5y23hgzsqHtA94yEJEQYD7wyrmFwOUO4GHjrEp7ReQA0A9Y09xzF5RWccvzqzlYWMHzs0bWfxhMG9KV3yzczmOf7OaxT3bTJTacCVlJTOyTxKQ+nYmL/OaFGsUVDu54aS12Rw2vfnvMNwqB6+/KPZOyuGdSVoPvxUeGcsXgrlw6IIVXvzrEXxfv4XR1LQ9c0Y87JmQSFmzjmvNS+e7Envz1sz2sySniX3eO4fzeSfXH6dslhjfmjOXL/YXsPFrK/hNlHDhRzqJtR3ljXS4hNmFsz0RSYsM5VeGguPI0+wrKOVlxmlvHZPDjy/sSF3H23+/6kekYY9ibX8b2IyUMSo2jV3LDBTIoSPj1lQNIiAzlL4t3s2DTEVJiw7iwTzK5RZV8ub+QmLBg7pzQg6SYME5WnKa4wuH8s9LB4VN2SiodjOmZyI8m96lvWQHERQQRFxHCPZOymD64G79csI3fvLeDp5fu47IBKUwekEKPpCge/WgXH2w9SnqnCAZ1i+MPH2bz4hcHmDU+k+eW76fWGF79zlhGZCQwc2Q69722gVn/XMOl/VPoFh9Op6hQEqPD6JkURe+UaJKjwzxuXeWX2NlxtAQDiOvftXfnaLr50MpUjWvqQ7uljB49+qx5+k8++STvvPMOALm5uezZs4fExLO7M3v06MGwYcMAGDFiBDk5OW7Ps2vXLnr06EGfPn0AuP3223n66af53ve+R3h4ON/+9reZNm0a06dPB2DChAnMnj2bmTNncu2113r0d6mpNThqmv7FP6DFwDUO8AKQbYx5vJHNDgGXACtEJAXoC+z3x/l/Nn8LuUWVvDh7FOOzvv5ATUuI5PnbR3GwsJwv9hbyxd4TfL7zOPM35GELEkZlJnBR385kJEaSEBlKQlQov3p3G4cKK3j5ztH0TolpVq4QWxC3j8/khlHpGAMRoWf3R2YkRvH4DcMa3V9EGN8rifG9vv47OWpqWX/wJJ/vzGfJznz25ZcRFxlKfEQIE7KSuOv8HgxLb/g3iLpj9k6J8ejvJiLMvbQ3N41OZ+nuApbtKmDRtmNEhQbz4NR+3DS6OzHhnl/52JDuiZG8fMcoPt1xnAWbjrBg4+H6wevwkCB+eFkf5kzsSXiIjVV7T/DwRzt5eNFO0jtF8PIdo+mZHA04i+eCe8/n0Y93smx3AWtziiiudJx1roTIEFITIqioqqHEXk2J3UFydBhjenRibM9EBnSLZW1OEYu2HmPtwSIaakwPTY9n6qAuTB7YhYxOkWd1n1Wcruar/UWsPlBIVnI0Vw7tRnhI433QqnWIioqq/3rp0qV89tlnfPnll0RGRjJp0qQG5/GHhX3dHWmz2aisrPzGNudqrHcmODiYNWvWsHjxYl5//XWeeuopPv/8c+bNm8dXX33FBx98wLBhw9i0adM3itK5qqqbHjyGwLcMJgC3AVtFZJPrtQeB7gDGmHnA/wEvichWnL9s/cwY0/SEWA+U2B2s2FPAnRN6nFUIzpSRGEVGYhQ3j+lOba1hU94pFmcfZ3F2Pg8t2vmN7f964zDG9WrewOaZ/PmBEGILYmzPRMb2TOTBqf3d7+AHnWPDmTkynZkj06mtNYj4d564iDB5oPMD1u6oYdW+E2w/XMK3hqeeNcg+PiuJBfdOYNW+Qvp3jf3G4HpEqO2s3zQdNbUUlp1mf0EZu46Xsvt4KUeL7UQnBhMbEUJMWDB5JytZtruAtzcert+vX5cY7r+kD+N6JRJiEwzO37jW5hTx0bZjPLRoJw8t2klYcBBpCRGkd4rE7qhh/cGTOGoMQQK1Bv74YTY3jOrOdcNTcdQY8kvtFJRWceSUnYOF5eQUlnO02M7E3sl8/9LePo1tKe/FxMRQWlra4HvFxcUkJCQQGRnJzp07Wb16td/O269fP3Jycti7dy9ZWVn8+9//5sILL6SsrIyKigqmTp3K2LFjycpy9jjs27ePMWPGMGbMGN577z1yc3PdFwNH0+MFEPjZRCtxfsA3tc0RYLK/z718dwGOGsOlHs4cCgoShndPYHj3BH5yeT8KSqvIL7VzstxBUcVpusSGM7pHJ3/HbDcCPWMnPMTGxf1SuLhfw/+eIsKERor+uUJsQXSJC6dLXHijvygA9V1n244UMyQtnl6u1sa5RmV24p5JWeQWVbBsdwEHC8vJLaok92QFInDnhB5c0DuZkZkJbDh0kn+tOsizy/cxb9m+bxyrW1w43RMjGZYezzsbD/POxsPcMrY790zK0kHwAEtMTGTChAkMGjSIiIgIUlK+/lmbMmUK8+bNY8iQIfTt25exY8f67bzh4eG8+OKLXH/99fUDyHfffTdFRUVcffXV2O12jDE88cQTAPzkJz9hz549GGO45JJLGDp0qNtzuBs8hhYaQPY3TwaQf/DGJpbuymfdLy/7xuCjUlY7fKqSFbsLiIsIITkmjM4x4XSODTurtXj4VCV/W7yHN9fnYRNheEY852clMSEriSFp8e3u5/rcAU/VPGd+P3NOlDPpsaUcfGS6dQPIVqiuqeXznflc0r9zu/sPo9qH1PgIbhzd9BTp1PgIHr5uCHMm9uSNtbms3HuiftJDvy4xPH3L8EZbK0qdyZOWQbssBusPnqS40sFl/Zt3cZlSrUHP5Gh+7hoHKiyrYsmuAv74YTZX/W0lf7x2MFcPS7U4oWrKvffeyxdffHHWa3PnzuWOO1pu5Z3WMIBsic+yjxNqC+ICnVeu2pnE6DBmjEhjQlYi9726kbmvb2LNgSJ+NX1Au5ihZIxpd4vVPf300y1+znO7/z1pGbTLVUsXZ+cztlfiNy6+Uqq96BoXwWtzxvLdC3vyyleHmPKX5Xy5r9DqWM0SHh5OYWGhLvjYTHX3MwgP//paKMtnE1lhX0EZ+0+UM3tCptVRlAqoEFsQP7+iPxf2TuaBt7dy03OruWl0Og9c0f8bFxa2BWlpaeTl5VFQUGB1lDav7k5nddwtUgftsBgszj4OwCU6XqA6iPFZSXx8/0Se+Gw3z6/Yz8JNR+jTJYbenaPp3TmGKYO6kN4p0v2BLBYSEuLVnbmU5zpkN9FnO/Lp3zVWL9RRHUpEqI0Hp/Znwb3nc92INMKDbXy+s4A/fJjNt55Zxf6CMqsjKgt1uAHkk+WnWXewiO9d9M21gZTqCAanxTE4La7++a5jpdz83Gpuff4r/nv3uLOu3FYdR4drGazce4JaAxdrF5FSgHNtpn/dNZqyqmpuef4r8ku8uy+uah88GTNoV8Vgb34ZIs41ZJRSTgO7xfHSnaNdq/h+xdJd+fX31lAdQ4drGRwqqqBrbHi7mG+tlD8N757A87NGcqzEzuwX1zLi958y9/WNLNutM3c6gg43tTSnsJyMxCj3GyrVAY3PSmLtLy7li70n+Hj7sfrlwWeNy+AX0/o3eOtT1T5UVdcQYmv6Yr52VQwOFlZw+UAdL1CqMeEhNi7pn8Il/VM4XV3Lox/t5PmVB9ice4qnbh7eJqagKu/ZHbVui3276SYqsTsoKj9N907aMlDKE6HBQfxy+gDm3TqC/SfKmfbkCj7ZfszqWCoAqqprCAtu+uO+3RSDQ4UVAGQm6m82SnljyqAuvH/f+XRPjGTOv9fzy3e3ejT7RLUdVdW1HacY5BSWA+iYgVI+yEiMYv7/jOc7F/TgP6sPceXfVpJ91L83pFfWqaqudTuxpt0Ug4OulkGGtgyU8klYsI1fTBvAv+4czckKB1c//QUfa7dRu1DlqCG0o7QMDhaWkxQdRpSuVKpUs0zsk8xH91/AgK6x3PPKBhZsOux+J9Wq2atrCetILQMdL1DKP5Kiw/jPt8cwMiOB+9/YxH/X5lodSTVDlaMDDSAfLKzQ8QKl/Cg6LJiX7hjNBb2T+en8Lfxj2T5qa/VeA21RhxkzsDtqOFZi1/ECpfwsItTGc7NGcMWgLjy0aCfXzVvFzmM6sNzWdJjZRIeKdPBYqUAJC7bxzC3DeXzmUA4WVjDtyZU8tChbp5+2IR2mmyjnhHNaaaZ2EykVECLCtcPTWPzDC5kxPI1/LNvPT9/aYnUs5SFny6ADdBPptFKlWkZCVCiPzBjCDy/rw8LNR3SmURtRVV1DeEg7bBkUlZ8+6/nBonLiIkKIjwy1KJFSHcs9k3pxXvd4fvXuNo4WV1odR7lR1V7XJjpabKe8qrr+uU4rVaplBduCeGLmMBw1hh+/uVlnGbVy9uoawtpjy6DWGD7YerT+eU5hOd11vECpFpWZFMWvpg/gi72FvPxljtVxVCNqag2OGtM+B5DDgoN4c53zIpjT1bUcPlmpLQOlLHDT6HQu7teZhxbtZOOhk1bHUQ2ou6tdu7zOICEylLU5J9lfUMbhU5XUGl2gTikriAh/mjGErnHh3PXyOvYXlFkdSZ2jqto5BbhdtgwSIkOxBQlvrs87Y7VSbRkoZYXE6DBevmM0Atz+4hryS+1WR1JnqLv/cbscQA62CRf1TWb++jz2F2gxUMpqmUlRvDB7FCdKT3PnS2spO2OCh7JW3cWB7bJlAHD9yHTyS6t4ZfVBIkNtJEeHWR1JqQ5tWHo8z9wynOyjpdz/+iaM0RlGrUFVex4zALi4X2eSokPZf6Kc7p0iEWn6Zs9KqcC7qF9nfn5FPz7LPs7CzUesjqNwXmMA7bhlEGIL4trhaYAuQ6FUa3LHhB4MS4/nNwu3U1hWZXWcDq9+ANnK6wxEJF1ElohItohsF5G5DWzzExHZ5HpsE5EaEenkyfFnjnQWg4wkHS9QqrWwBQmPzhhCWVU1v31vh9VxOjy7o3UMIFcDPzLG9AfGAveKyIAzNzDG/MkYM8wYMwz4ObDMGFPkycGzOsfw9M3DmT0+09+5lVLN0Cclhu9d1JuFm4/w2Y7jVsfp0Pw6tVREgkTkPBGZJiIXi0iKJ/sZY44aYza4vi4FsoHUJna5CXjNk2PXmTakK13jIrzZRSnVAv5nUi/6psTwy3e3UWJ3WB2nw/LLALKI9BKRZ4G9wMM4P6zvAT4VkdUicoeIeFpQMoHzgK8aeT8SmALMb+T9OSKyTkTWFRQUeHJKpZSFQoODeGTGEPJL7dz7yob631BVy/JXy+D3wH+AXsaYy40xtxpjZhhjhgBXAXHAbe7CiEg0zg/5+40xjd0m6Urgi8a6iIwxzxpjRhpjRiYnJ7s7pVKqFRiWHs+jM4ayYs8Jvv/aRqpraq2O1OHUjxm4GUAObupNY8xNTbyXD/zFXRARCcFZCF4xxrzdxKY34mUXkVKq9ZsxIo0yu4PfvLeDn761hceuH0pQkE4FbylV9RedNd1N1GQxEJFrm3rfzYc74pz8/wKQbYx5vInt4oALgVubOp5Sqm2aPaEHpfZq/vzpbqLDg/nd1YOsjtRhfD1m0IyWAc6uG4DOwHjgc9fzi4ClQJPFAJiAsxtpq4hscr32INAdwBgzz/Xat4BPjDHlbo6nlGqjvndxFicrHPzziwNMHdyVsT0TrY7UIdQVg1Bb87qJ7gAQkfeBAcaYo67nXYGn3YUwxqwE3LYHjTEvAS+5204p1XaJCD+d0pf3txzh8U9288Z3x+rKAS2gqrqG4CAh2E0x8PQ6g8y6QuByHOjjazilVMcUHmLjvouzWJNTxIo9J6yO0yHYHbVuZxKB58VgqYh8LCKzReR24ANgSXMCKqU6ppmj0kmNj+DPn+zSxexaQFV1jdtrDMDDYmCM+R4wDxgKDAOeNcbc16yESqkOKSzYxtxLerM5r5hP9erkgKvyc8sAYAPwgTHmB8DHIhLjazilVMd27fBUeiRF8finu6mt1dZBIFVV1xLmr5aBiHwHeAv4h+ulVOBdn9MppTq0YFsQ91/am53HSnl/61H3Oyif2R01fm0Z3ItzmmgJgDFmD87ppkop5ZMrh3SjX5cYHv1oJ5WndamKQPFrywCoMsacrnsiIsGAtu2UUj4LChJ+feVA8k5W8szSvVbHabeqqv3bMlgmIg8CESJyGfAm8F4z8imlFON6JfKt81L5x7L97C8oszpOu1RV7d8B5AeAAmAr8F3gQ+CXPqdTSimXn0/tR1hwEL9euF2nmgaA8zoD/00trTXGPGeMud61aulzRv/VlFJ+0DkmnB9f3pcVe07w4dZjVsdpd5zXGfipZSAiE0TkUxHZLSL7ReSAiOxvdkqllAJuHZvBwG6x/O797ZRVVVsdp12p8mfLAOfKo48D5wOjgJGuP5VSqtlsQcLvrh7E8ZIq/rs21+o47YpzNpH/xgyKjTGLjDH5xpjCukfzIiql1NdGZCQwJC2ON9fnWR2lXfHLbCIRGS4iw4ElIvInERlX95rrdaWU8pvrR6SRfbSEbYeLrY7SblQ5aj1am8jd/Qz+fM7zkWd8bYCLvcyllFKNumpoKv/3fjZvrc9jUGqc1XHavNpaw+kaz6aWurufwUUAItLTGHPWgLGI9GxWSqWUOkdcZAiXDUzh3U2HXVNO3f9Gqxp32nXPaX8OIL/VwGtveh5JKaU8c/2INE5VOPg8O9/qKG1elaOuGDSzZSAi/YCBQNw590OOBcJ9j6iUUg27oHcyXWLDeXN9HlcM7mp1nDbNXu1c88mT2UTuxgz6AtOBeL6+HzJAKfAd3+IppVTjbEHCtcNTmbdsH/kldjrH6u+dvqprGYR70E3kbsxgAbBARMYZY770SzqllHJjxog0nlm6j7c3HubuC3tZHafNqvKiZeDpmEGuiLwjIvkiclxE5otIWjMyKqVUo3omRzMyI4E31+XqekXNUFXt/wHkF4GFQDecN7Z5z/WaUkoFxA2j0tlXUM6X+/T6Vl/ZHa6WgR9XLe1sjHnRGFPterwEJPucUCml3LhyaDc6RYXy0qocq6O0WXUtA08uOvO0GBSIyK0iYnM9bgW0XCulAiY8xMaNo9L5LPs4eScrrI7TJtWPGfixZXAnMBM45nrMcL2mlFIBc+vYDESEf68+aHWUNqn+OgN/DSAbYw4ZY64yxiS7HtcYY/RfRykVUN3iI7h8YApvrM3V+yT7wO8DyCKSprOJlFJWuH1cJqcqHCzYdNjqKG1O3QCy325ug84mUkpZZHSPTvTrEsNLq3J0mqmXAjG1NFlnEymlrCAizB6fyc5jpaw5UGR1nDYlEAPIJ3Q2kVLKKlcPSyUuIoTX1hyyOkqb4s1Cdb7MJjqKziZSSrWgiFAbUwZ24bPs/Pp+cOWevbqG4CAh2BaY2USddTaRUqqlTR3SlbKqalbsOWF1lDajyuHZjW3A/aqlAIhID+A+IPPMfYwxV/mQTymlvDa+VyJxESEs2nqUywakWB2nTaiqriXMg6uPwcNiALwLvIBzFlGtj7mUUspnIbYgJg9I4aNtx1w3ede7oLnj/D551jLwdMzAbox50hizxBizrO7he0SllPLe1CFdKa2qZqV2FXnE7qj1aF0i8LwY/FVEfi0i40RkeN3D3U4iki4iS0QkW0S2i8jcRrabJCKbXNtokVFKNWhCryRiw4P5YOtRq6O0Cd60DDztJhoM3AZczNfdRMb1vCnVwI+MMRtEJAZYLyKfGmN21G0gIvHAM8AUY8whEensYSalVAcTGhzEZQO68MkO7SryRFW15wPInrYMvgX0NMZcaIy5yPVwVwgwxhw1xmxwfV0KZOO8gvlMNwNvG2MOubbTu2ArpRo1bUgXSu3VrNqrlzq545xN5N9uos0474PsMxHJBM4DvjrnrT5AgogsFZH1IjKrkf3niMg6EVlXUFDQnChKqTbs/KxkYrSryCNV1TUerVgKnncTpQA7RWQtUFX3oqdTS0UkGpgP3G+MKWkgwwjgEiAC+FJEVhtjdp+5kTHmWeBZgJEjR+oCJUp1UM6uohQ+2X6M098aTKiH3SAdkd1RS6co/04t/bWvYUQkBGcheMUY83YDm+QBJ4wx5UC5iCwHhgK7G9hWKaWYOqgrb284zBd7T3BRPx1mbIw3LQNPS+o6YIVrOulRIA5Y5W4nERGc1ydkG2Meb2SzBcAFIhIsIpHAGJxjC0op1aCJfZKJiwhh4eYjVkdp1bwZQPa0ZbAc5wd2ArAYZ3G4AbjFzX4TcM5C2ioim1yvPQh0BzDGzDPGZIvIR8AWnDOVnjfGbPMwl1KqAwoNDmLq4C4s2HSEytM1RITqrKKGOIuBf7uJxBhTISJ3AX8zxjx6xod7o4wxKwHxYLs/AX/yMItSSnHl0G68tiaXz7KPc+XQblbHaZXsjhqPbmwDnncTiYiMw9kS+MD1mpZipZRlxvRIJCU2TLuKmuBNy8DTYjAX+DnwjjFmu4j0BJb4mE8ppZrNFiRcOaQbS3flU1zhsDpOq2OM4bS/Lzozxix3LWH9iOv5fmPM95uRUymlmu2qYd1w1Bg+2q7XHJyr/paX/ugmEpFnRWRwI+9FicidIuJuEFkppQJicGocPZKiWLBJu4rOVXeXs3A/DSA/A/zKVRC2AQVAONAbiAX+CbziY1allGoWEeGqod148vM9HC+xk+S2l0IAABSXSURBVBIbbnWkVqP+/sf+aBkYYzYZY2YCo4CngRXAQuDbxpihxpi/GmOqmjqGUkoF0lXDumEMvL9Fu4rOVN9N5M+ppcaYMmCpr6GUUipQeiVHMyg1loWbDnPX+T2sjtNq1LcM/LxqqVJKtVpXDunG5rxicosqrI7Satjrxgz8fHMbpZRqtaYO7gqgK5meIeAtAxEJEpFYb/dTSqlASe8UydC0OD7UYlCvbjaRX4uBiLwqIrEiEgXsAHaJyE98DamUUv42bUhXtuQVc6hQu4rgzOsM/NtNNMB1H4JrgA9xLjR3mw/5lFIqIK4Y5Owq+nCbtg7g624if69NFOK6L8E1wAJjjAPnPZCVUqpVSO8UydD0eD7QKabA1wPI/l6b6B9ADhAFLBeRDODcO5YppZSlpg3uwtbD2lUEARpANsY8aYxJNcZMNU4HgYt8TqmUUgGgs4q+9vVFZ/4dQJ7rGkAWEXlBRDYAF/ucUimlAiAtwdVVtFXXKqoK0HUGd7oGkCcDycAdwMM+5FNKqYCaPrgr2w6XcLCw3Ooolqp0BOY6g7q7lU0FXjTGbMaDO5gppVRLu2JwFwDe6+A3vSm1O4gIsRFs828xWC8in+AsBh+LSAzO+xUrpVSrkpYQybieiby+Npea2o476bGkspq4iBCPt/e0GNwFPACMMsZUAKE4u4qUUqrVuXVsBnknK1m2O9/qKJYpsTuIjfD0NveezyaqBdKAX4rIY8B4Y8wW3yIqpVRgTR6YQnJMGP9ZfcjqKJYprnQQG+7nloGIPIzzPsg7XI/vi8hDPiVUSqkAC7EFcdOodJbsyu+wK5k6Wwb+7yaaClxmjPmnMeafwBRgmg/5lFKqRdw4ujsCvLqmY7YOAjVmABB/xtdxXuynlFItrlt8BJf2T+GNtbn1V+N2JCV2B7Hhfh4zAB4CNorISyLyMrAe+KMP+ZRSqsXcOjaDovLTfLTtmNVRWlRtraGkMgDdRMaY14CxwNuuxzhjzOs+pVRKqRZyflYSmYmR/Gf1QaujtKjy09XUGvw3gCwiw+seQFcgD8gFurleU0qpVisoSLhlTAZrc06y8dBJq+O0mBJ7NYBXYwbuOpT+3MR7Bl2fSCnVyt00pjv/WL6Phxbt5I05YxFp/4snlFQ6ALy6zqDJLY0xHq1MKiKXGWM+9fisSinVQqLDgpl7aR9+9e42PsvO57IBKVZHCrjiumLg7+sMPPCIn46jlFJ+d+OodHomR/Hwomyqa9r/Sjpftwxavhi0/3aXUqrNCrEF8bMp/dhXUM5/1+VZHSfgfBkz8Fcx6LirQSml2oTJA1IYmZHAE5/tpryq2uo4AWVlN5FSSrVqIsLPp/anoLSK51bstzpOQNV1E0UH4KIzd3L8dByllAqYERkJXDGoC88u309+id3qOAFTYncQEx6MLcjzHnyPi4GIjBeRm0VkVt2j7j1jzLWN7JMuIktEJFtEtovI3Aa2mSQixSKyyfX4X4/TK6WUl342pR+nq2t54rPdVkcJmJLKaq+6iMD9dQYAiMi/gV7AJqBukQ8D/MvNrtXAj4wxG1w3xFkvIp8aY3acs90KY8x0L3IrpZRPMpOiuG1cBi+vymH2+B707RJjdSS/K/ZyKQrwsBgAI4EBxhivBoqNMUeBo66vS0UkG0jFuQy2UkpZ4vsX92b++jweWpTNS3eMtjqO33m7SB143k20DejidaIziEgmcB7wVQNvjxORzSKySEQGNuc8SinlTkJUKPdd3JuluwpYsafA6jh+V1Lp8GpaKbhfm+g9EVkIJAE7RORjEVlY9/D0JCISDcwH7jfGlJzz9gYgwxgzFPgb8G4jx5gjIutEZF1BQfv7x1NKtaxZ4zNIS4jgDx9kt7t7JZfaq/3eTfSY73GcRCQEZyF4xRjz9rnvn1kcjDEfisgzIpJkjDlxznbPAs8CjBw5sn39yymlWlxYsI2fTenHfa9tZP76PGaOSrc6kt94e8tLcNMyMMYsM8YsAw4BX53xfA3gdk1Yca4I9QKQbYx5vJFturi2Q0RGuzIVevW3UEopH0wf0pXzusfzp092UdZOLkSrrqmlrKraq0XqwPMxgzeBMxf0qHG95s4E4Dbg4jOmjk4VkbtF5G7XNjOAbSKyGXgSuNHbgWqllPKFiPCr6QMoKK1i3tJ9Vsfxi7qi5u2YgaelI9gYc7ruiTHmtIiEutvJGLMSN+sWGWOeAp7yMIdSSvnV8O4JXD2sG8+t2M+No9NJS4i0OlKzlFQ6i4Ffu4nOUCAiV9U9EZGrgRNNbK+UUm3GT6f0A+CRj3ZZnKT5in1YsRQ8LwZ3Aw+KyCERyQV+BnzXqzMppVQrlRofwZyJPXlv8xHWHyyyOk6zlNjrFqkLwJiBMWafMWYsMADnxWfjjTF7vQ2plFKt1d0X9qJzTBi/ez+btjxsWbdIXVxkYMYMEJFpwEAgvO62ccaY33l1NqWUaqWiwoK5/9I+PPjOVtYdPMmozE5WR/KJL8tXg4ctAxGZB9wA3IdzQPh6IMOrMymlVCt39bBuRIbamL++7d4Ap76bKEBjBuONMbOAk8aY3wLjgPZzhYZSSuFsHVwxqCvvbzlK5eka9zu0QiWV1diChKhQm1f7eVoMKl1/VohIN8AB9PDqTEop1QZcNyKVsqpqPtlxzOooPqlbpK6uO99TnhaD90UkHvgTzrWEcoDXvDqTUkq1AWN7JJIaH8FbbbSryJflq8Hz2UT/Z4w5ZYyZj3OsoJ8xRm9Co5Rqd4KChOtGpLFy7wmOFle636GVKfFhXSLwfAA5XER+KCJvA68Cd4pIuNdnU0qpNuC64akYA29vOGx1FK+V2L1flwg87yb6F85ppX/DuXREf+DfXp9NKaXagIzEKEZndmL+hrw2d82BL/cyAM+LQV9jzF3GmCWuxxygj9dnU0qpNuK6EansLyhnY+4pq6N4xZflq8HzYrBRRMbWPRGRMcAXXp9NKaXaiKmDuxIeEsSb69rWQHKJPQADyCKyVUS2AGOAVSKSIyIHgC+BiT4lVUqpNiAmPIQrh3TjnY15nCw/7X6HVqCquga7o9brdYnA/XIU032LpJRSbd93JvbkzfV5/Gf1Qe67pLfVcdwqtft2LwNwUwyMMW7vZqaUUu1Vn5QYJvVN5uUvc/jOxJ6Eh3h3VW9L83X5avB8zEAppTqkORf05ETZad7d2PqnmZb4uEgdaDFQSqkmjeuVyMBusTy3Yj+1ta17mmmJq5sokNcZKKVUhyQizJnYk30F5SzZlW91nCbV38tAu4mUUsr/pg7uSre4cP6xfL/VUZrk670MQIuBUkq5FWIL4s7ze7DmQBHrclrvbTF9vZcBaDFQSimP3Di6O11iw/np/C2t9l4HJZXVhNqCCAv2/qNdi4FSSnkgOiyYP88cyv6Cch5elG11nAbVLV/t7b0MQIuBUkp5bEJWEndO6MHLXx5k2e4Cq+N8g3MpCu9nEoEWA6WU8spPp/Sld+dofvLm5la3TIWv9zIALQZKKeWV8BAbT9wwjJMVp/nlu9usjnMW570MtBgopVSLGJQaxz2Tsvhg61F2Hy+1Ok49X+9lAFoMlFLKJ7PGZRBiE15fk2t1lHrObiIdM1BKqRaTGB3G5AFdeHtjHnaH9VNNjTE+38sAtBgopZTPbhydzqkKBx9vP2Z1FOyOWhw1RgeQlVKqpU3olUR6p4hW0VV04EQ54Nu6RKDFQCmlfBYUJNw4qjtf7i8kx/VhbIWC0iru/s96OkWFcmHfZJ+OocVAKaWa4foRadiChNfXWtM6KK+q5q6X15JfaueF20eSGh/h03G0GCilVDN0jg3n4n6deWt9Ho6a2hY9t6Omlntf3cC2w8U8ddNwzuue4POxtBgopVQz3TQ6nRNlVSzOPt6i531k0U6W7irg99cM5tIBKc06VkCLgYiki8gSEckWke0iMreJbUeJSI2IzAhkJqWU8rcL+3QmNT6Cecv2Y0zL3A3NGMO7m44wbXBXbh7TvdnHC3TLoBr4kTGmPzAWuFdEBpy7kYjYgEeAjwOcRyml/M4WJMy9tDebck/x3pajLXLOQ0UVnCirYnxWol+OF9BiYIw5aozZ4Pq6FMgGUhvY9D5gPtC67ymnlFKNuG54Gv27xvLIop0tchHa+oMnARiR4fs4wZlabMxARDKB84Cvznk9FfgWMM/N/nNEZJ2IrCsoaH1LxyqlOjZbkPDLaf05fKqSl1blBPx86w+eJCYsmN6dY/xyvBYpBiISjfM3//uNMSXnvP0X4GfGmCZLqTHmWWPMSGPMyORk3+bRKqVUIE3ISuLifp15+vO9FJZVBfRc6w+e5LyMBGxB3t/IpiEBLwYiEoKzELxijHm7gU1GAq+LSA4wA3hGRK4JdC6llAqEB6f2o8JRw18X7wnYOUrsDnYdL2VEM6aSnivQs4kEeAHINsY83tA2xpgexphMY0wm8BZwjzHm3UDmUkqpQMnqHMNNo9N55atDLN0VmGHQTYdOYYz/xgsg8C2DCcBtwMUissn1mCoid4vI3QE+t1JKWeKHl/UlKzmaO15ay+Of7KKm1r/TTdcfPEmQwLDu8X47pm8LX3vIGLMS8LhDyxgzO3BplFKqZXSKCuXdeyfwqwXbePLzvaw/dJK/3ngeSdFhfjn+hkMn6dcllugw/32E6xXISikVABGhNh67fiiPXjeEdTknmfH3VX6ZclpTa9h46JRfu4hAi4FSSgXUzFHpPDdrJDmFFfzziwPNPt6uY6WUVVVrMVBKqbZmYp9kLu2fwjNL9lFQ2rwpp+sP+fdiszpaDJRSqgU8OLUfdkcNj3+6q1nHWZ9TROeYMNISfFuqujFaDJRSqgX0TI5m1rhM3libS/bRc6+99dz6QycZkZGAc+a+/2gxUEqpFjL3kt7ERoTw+w92+LS6aX6JndyiSr93EYEWA6WUajFxkSH84NI+fLG3kN9/kM3y3QWU2B0e778mpwjw/3gBBPg6A6WUUme7eUx3luzK559fHOCFlQcQgcGpcTxzy3DSEiIb3c8Yw/MrDpAaH8Gg1Di/59JioJRSLSjEFsRLd4ym1O5gc24xGw6d5Lnl+/neqxv573fHERrccIfN0t0FbMo9xUPXDibE5v9OHe0mUkopC8SEh3B+7yS+f0lvHpkxhE25p3j0o50NbmuM4YlPd5PeKYIZI9ICkkeLgVJKWWzq4K7cPi6D51ce4JPtx77x/uLsfLbkFXPfRb0D0ioALQZKKdUqPDitP4NSY/nxm5vJLaqof90YwxOf7SYjMZJvDW/oRpH+ocVAKaVagbBgG0/fPBxj4Pp5X/L4p7s5WFjOJzuOs/1ICd+/OHCtAgDxZa6r1UaOHGnWrVtndQyllPK7tTlFPLl4Dyv3nsAYiAy1kRIbzqc/mEhwM4uBiKw3xoxs6D2dTaSUUq3IqMxO/PuuMRwtruSdjYf5ePtx5l6S1exC4I62DJRSqoNoqmWgYwZKKaW0GCillNJioJRSCi0GSiml0GKglFIKLQZKKaXQYqCUUgotBkoppWijF52JSCngy12l44BiL99395q7r+v+TAJOaGbNrJk1s4WZM4wxyQ3uYYxpcw9gnY/7Pevt++5ec/f1GX9qZs2smTVzq8jc0KOjdRO958P77l5z97W7c7qjmT37WjN79r5m9l57zPwNbbWbaJ1pZH2N1koztwzN3DI0c8toycxttWXwrNUBfKCZW4ZmbhmauWW0WOY22TJQSinlX221ZaCUUsqPtBgopZTSYqCUUqodFgMRuUBE5onI8yKyyuo8nhCRIBH5g4j8TURutzqPJ0RkkoiscH2vJ1mdx1MiEiUi60VkutVZPCEi/V3f47dE5H+szuMJEblGRJ4TkQUiMtnqPJ4QkZ4i8oKIvGV1lqa4fn5fdn1/b/HnsVtVMRCRf4pIvohsO+f1KSKyS0T2isgDTR3DGLPCGHM38D7wciDzurI1OzNwNZAKOIC8QGU9I5s/MhugDAin7WQG+Bnw38CkPJuffp6zXT/PM4GATzH0U+Z3jTHfAWYDNwQwbl02f2Teb4y5K7BJG+Zl/muBt1zf36v8GsSXq9sC9QAmAsOBbWe8ZgP2AT2BUGAzMAAYjPMD/8xH5zP2+y8Q2xYyAw8A33Xt+1YbyRzk2i8FeKWNZL4UuBHnh9T0tpDZtc9VwCrg5raS2bXfn4HhbSxzwP//NTP/z4Fhrm1e9WeOYFoRY8xyEck85+XRwF5jzH4AEXkduNoY8xDQYFNfRLoDxcaYkgDGBfyTWUTygNOupzWBS+vkr++zy0kgLBA5z+Sn7/NFQBTO/1SVIvKhMaa2NWd2HWchsFBEPgBeDVRe17n88X0W4GFgkTFmQyDzgt9/nlucN/lxtsLTgE34uWenVRWDRqQCuWc8zwPGuNnnLuDFgCVyz9vMbwN/E5ELgOWBDNYErzKLyLXA5UA88FRgozXKq8zGmF8AiMhs4EQgC0ETvP0+T8LZNRAGfBjQZI3z9uf5PpytsDgRyTLGzAtkuEZ4+31OBP4AnCciP3cVDSs1lv9J4CkRmUbzl6w4S1soBtLAa01eKWeM+XWAsnjKq8zGmAqcBcxK3mZ+G2cRs5LXPxsAxpiX/B/FY95+n5cCSwMVxkPeZn4S54eWlbzNXAjcHbg4XmswvzGmHLgjECdsVQPIjcgD0s94ngYcsSiLpzRzy9DMLUMzt7wWz98WisFaoLeI9BCRUJwDgAstzuSOZm4ZmrllaOaW1/L5W3rk3M2o+mvAUb6eYnmX6/WpwG6co+u/sDqnZtbMmtn6rG05c2vMrwvVKaWUahPdREoppQJMi4FSSiktBkoppbQYKKWUQouBUkoptBgopZRCi4FSfiEiZVZnUKo5tBgoFSAiYrM6g1Ke0mKglB+J8w5wS0TkVWCr1XmU8lRbWLVUqbZmNDDIGHPA6iBKeUpbBkr53xotBKqt0WKglP+VWx1AKW9pMVBKKaXFQCmlFLqEtVJKKW0ZKKWUQouBUkoptBgopZRCi4FSSim0GCillEKLgVJKKbQYKKWUQouBUkop4P8BDdxPavG8WZYAAAAASUVORK5CYII=\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=8)\n",
    "_ = lr_finder.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10722672220103299"
      ]
     },
     "execution_count": 15,
     "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": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:\t[0s / 0s],\t\ttrain_loss: 2.6413,\tval_loss: 2.5174\n",
      "1:\t[0s / 0s],\t\ttrain_loss: 2.3406,\tval_loss: 2.3192\n",
      "2:\t[0s / 0s],\t\ttrain_loss: 2.1494,\tval_loss: 2.0976\n",
      "3:\t[0s / 0s],\t\ttrain_loss: 1.9278,\tval_loss: 1.8521\n",
      "4:\t[0s / 0s],\t\ttrain_loss: 1.7077,\tval_loss: 1.6468\n",
      "5:\t[0s / 0s],\t\ttrain_loss: 1.5645,\tval_loss: 1.5185\n",
      "6:\t[0s / 0s],\t\ttrain_loss: 1.4957,\tval_loss: 1.4789\n",
      "7:\t[0s / 0s],\t\ttrain_loss: 1.4887,\tval_loss: 1.4919\n",
      "8:\t[0s / 0s],\t\ttrain_loss: 1.4821,\tval_loss: 1.4874\n",
      "9:\t[0s / 0s],\t\ttrain_loss: 1.4499,\tval_loss: 1.4807\n",
      "10:\t[0s / 0s],\t\ttrain_loss: 1.4514,\tval_loss: 1.4747\n",
      "11:\t[0s / 0s],\t\ttrain_loss: 1.4425,\tval_loss: 1.4793\n",
      "12:\t[0s / 0s],\t\ttrain_loss: 1.4182,\tval_loss: 1.4848\n",
      "13:\t[0s / 0s],\t\ttrain_loss: 1.4268,\tval_loss: 1.4878\n",
      "14:\t[0s / 0s],\t\ttrain_loss: 1.4193,\tval_loss: 1.4861\n",
      "15:\t[0s / 0s],\t\ttrain_loss: 1.4099,\tval_loss: 1.4978\n",
      "16:\t[0s / 0s],\t\ttrain_loss: 1.3999,\tval_loss: 1.5006\n",
      "17:\t[0s / 0s],\t\ttrain_loss: 1.4148,\tval_loss: 1.4996\n",
      "18:\t[0s / 0s],\t\ttrain_loss: 1.4060,\tval_loss: 1.5080\n",
      "19:\t[0s / 0s],\t\ttrain_loss: 1.3819,\tval_loss: 1.5106\n",
      "20:\t[0s / 0s],\t\ttrain_loss: 1.3740,\tval_loss: 1.5111\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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXwV5bnA8d9zck72lSRsCYR9kS1ARBRBAeuCXnGr+0brpShabLVXrba1Xr1tr631WlHUuraIKCpa61JFFGVtgmENu4EEkITsC9nf+8ecHJKYkISck0k4z/fzmc/MmXnPmYfh5Jn3vPPOO2KMQSmllH9x2B2AUkqpzqfJXyml/JAmf6WU8kOa/JVSyg9p8ldKKT/ktGvHcXFxZsCAAXbtXimluqW0tLSjxpj4jn6Obcl/wIABpKam2rV7pZTqlkRkvzc+R5t9lFLKD2nyV0opP6TJXyml/JBtbf5KqVNPdXU12dnZVFRU2B1KtxccHExiYiIul8snn6/JXynlNdnZ2URERDBgwABExO5wui1jDHl5eWRnZzNw4ECf7EObfZRSXlNRUUFsbKwm/g4SEWJjY336C0qTv1LKqzTxe4evj6NtyT+npNKuXSullN+zLfkfKa4gp1gvCimllB1sbfZ5f9MhO3evlDrFFBYW8swzz7T7fbNmzaKwsLDd77v11ltZtmxZu9/XFdiW/ENcASxPP2jX7pVSp6CWkn9tbe0J3/fhhx8SHR3tq7C6JNu6ekaHuth6sJg9OaUM6RluVxhKKR/57T+2sf1QsVc/87S+kfzmP0a1uP3+++9n7969JCcn43K5CA8Pp0+fPqSnp7N9+3Yuu+wysrKyqKioYMGCBcydOxc4PtZYaWkpF110EWeffTZr1qwhISGB9957j5CQkFZjW7FiBffeey81NTWcfvrpPPvsswQFBXH//ffz/vvv43Q6Of/88/njH//IW2+9xW9/+1sCAgKIiopi1apVXjtGbdVqzV9E+onIShHJEJFtIrKghXLniki6u8yXrX1udGggDoH3tPavlPKS3//+9wwePJj09HQef/xxNmzYwGOPPcb27dsBeOmll0hLSyM1NZWnnnqKvLy8733G7t27mT9/Ptu2bSM6Opq333671f1WVFRw6623snTpUrZs2UJNTQ3PPvss+fn5vPvuu2zbto3Nmzfz0EMPAfDII4/wySefsGnTJt5//33vHoQ2akvNvwa4xxizUUQigDQR+dQYs72+gIhEA88AFxpjDohIz1Z37BCmDIljefpBfv6DYdo9TKlTzIlq6J1l0qRJjW6Seuqpp3j33XcByMrKYvfu3cTGxjZ6z8CBA0lOTgZg4sSJZGZmtrqfnTt3MnDgQIYNGwbALbfcwsKFC7nzzjsJDg7mtttu4+KLL+aSSy4BYMqUKdx6661cffXVXHHFFd74p7ZbqzV/Y8xhY8xG93IJkAEkNCl2PfCOMeaAu1xOW3Z+WXICWfnH2HigoH1RK6VUG4SFhXmWv/jiCz777DPWrl3Lpk2bGD9+fLM3UQUFBXmWAwICqKmpaXU/xphm1zudTjZs2MCVV17J8uXLufDCCwFYtGgRjz76KFlZWSQnJzf7C8TX2nXBV0QGAOOB9U02DQNiROQLEUkTkZtbeP9cEUkVkdTc3FwuGN2bYJeDd7/Rph+lVMdFRERQUlLS7LaioiJiYmIIDQ1lx44drFu3zmv7HTFiBJmZmezZsweAv/3tb5xzzjmUlpZSVFTErFmzePLJJ0lPTwdg7969nHHGGTzyyCPExcWRlZXltVjaqs0XfEUkHHgbuNsY0/QqjhOYCMwEQoC1IrLOGLOrYSFjzPPA8wApKSkmPMjJD07rzT83H+Y3/zEKV4DecKyUOnmxsbFMmTKF0aNHExISQq9evTzbLrzwQhYtWsTYsWMZPnw4kydP9tp+g4ODefnll/nhD3/oueA7b9488vPzmT17NhUVFRhj+POf/wzAL37xC3bv3o0xhpkzZzJu3DivxdJW0tLPlUaFRFzAB8Anxpgnmtl+PxBsjHnY/fpF4GNjzFstfWZKSopJTU1lRcYRfvxqKi/eksLMkb1aKq6U6gYyMjIYOXKk3WGcMpo7niKSZoxJ6ehnt6W3jwAvAhnNJX6394CpIuIUkVDgDKxrA62aNiyemFCXNv0opVQnakuzzxTgJmCLiKS71/0S6A9gjFlkjMkQkY+BzUAd8FdjzNa2BOAKcHDJ2L68lZZFaWUN4UE6yrRSqmuZP38+q1evbrRuwYIFzJkzx6aIOq7VTGuM+RpotR+mMeZx4PGTCeKy8Qn8bd1+Ptn6HVdOTDyZj1BKKZ9ZuHCh3SF4XZe4wjqhfzT9e4TqcA9KKdVJukTyFxFmJ/dl9Z6jOtKnUkp1gi6R/AFmJydQZ3SkT6WU6gxdJvkP6RnOmIQo3kvX5K+UUr7WZZI/wOzkvmw5WMSenFK7Q1FK+Ynw8JZHFc7MzGT06NGdGE3nsS/5V5d/b9Wl4/rqSJ9KKdUJ7OtUn7cHyvMhtIdnVc/IYKYMieO99EM60qdS3d1H98N3W7z7mb3HwEW/P2GR++67j6SkJO644w4AHn74YUSEVatWUVBQQHV1NY8++iizZ89u164rKiq4/fbbSU1Nxel08sQTTzB9+nS2bdvGnDlzqKqqoq6ujrfffpu+ffty9dVXk52dTW1tLb/61a+45pprTvqf7Qv21fzr6uCrP31v9WXJCRzIL2fjgfY/Uk0ppa699lqWLl3qef3mm28yZ84c3n33XTZu3MjKlSu55557WhyJsyX1ff23bNnCkiVLuOWWW6ioqGDRokUsWLCA9PR0UlNTSUxM5OOPP6Zv375s2rSJrVu3ekbz7Ersq/mH9oANz8OkuRCT5Fl9wejePLh8C8u/OcjEpBjbwlNKdVArNXRfGT9+PDk5ORw6dIjc3FxiYmLo06cPP/vZz1i1ahUOh4ODBw9y5MgRevfu3ebP/frrr7nrrrsAaxTPpKQkdu3axZlnnsljjz1GdnY2V1xxBUOHDmXMmDHce++93HfffVxyySVMnTrVV//ck2ZfzT+iD4gDVj7WaHV4kJPzRvbig82HqK6tsyk4pVR3dtVVV7Fs2TKWLl3Ktddey+LFi8nNzSUtLY309HR69erV7Fj+J9LSL4Xrr7+e999/n5CQEC644AI+//xzhg0bRlpaGmPGjOGBBx7gkUce8cY/y6vsS/4BLjhjHmx+Ew5vbrTp8vEJFJRXs2pXrk3BKaW6s2uvvZY33niDZcuWcdVVV1FUVETPnj1xuVysXLmS/fv3t/szp02bxuLFiwHYtWsXBw4cYPjw4ezbt49Bgwbx05/+lEsvvZTNmzdz6NAhQkNDufHGG7n33nvZuHGjt/+JHWZvV8+z74bgKPjs4Uar60f6XK59/pVSJ2HUqFGUlJSQkJBAnz59uOGGG0hNTSUlJYXFixczYsSIdn/mHXfcQW1tLWPGjOGaa67hlVdeISgoiKVLlzJ69GiSk5PZsWMHN998M1u2bGHSpEkkJyfz2GOPeZ7d25W0aTx/X6gfz581f4F/PQQ3vweDzvVsf2j5FpalZZP60A90pE+lugkdz9+7bB3P3+dO/0+I6gef/sbqAeR2+fgEKqrr+GTrdzYGp5RSpyb7k78rGKY/CIfTYds7ntUT+sfQr0eIjvSplPK5LVu2kJyc3Gg644wz7A7Lp7pGe8rYq2Ht0/D5f8PIS8EZiIhwWXICC1fuIaekgp4RwXZHqZRqA2NMt7tBc8yYMZ6Hq3cVvm6St7/mD+AIgPMehoJMSHvZs7p+pM9/bDpsV2RKqXYIDg4mLy/P54nrVGeMIS8vj+Bg31V6u0bNH2DIeTBgKnz5Bxh3HQRHMqRnOKMTInkv/SA/Pnug3REqpVqRmJhIdnY2ubnaTbujgoODSUz03ZMNu07yF4Ef/BZemGH1AJrxIGAN9/DoPzPYm1vK4PiWR99TStnP5XIxcKBW1LqDrtHsUy9hIoy63Gr/L7F6+XhG+vxGL/wqpZS3dK3kDzDjV1BbZTX/cHykz+Xph7QdUSmlvKTrJf/YwTBxDqS9Ckd3A9aFXx3pUymlvKfV5C8i/URkpYhkiMg2EVlwgrKni0itiFzVoajOuQ9cIbDitwBcMKoXQU6HPuRFKaW8pC01/xrgHmPMSGAyMF9ETmtaSEQCgD8An3Q4qvB4OOsuyPgHZP2biGAXPzitFx9sPqwjfSqllBe0mvyNMYeNMRvdyyVABpDQTNG7gLeBHK9EduadENYTPv01GMNlyQnkl1Xx1W7tQqaUUh3VrjZ/ERkAjAfWN1mfAFwOLGrl/XNFJFVEUlvtBxwUDufeBwfWwK6PmTYsnuhQF8u/0ZE+lVKqo9qc/EUkHKtmf7cxprjJ5ieB+4wxtSf6DGPM88aYFGNMSnx8fOs7nXAL9BgMnz1MoMNwydg+/Gv7d5RW1rQ1bKWUUs1oU/IXERdW4l9sjHmnmSIpwBsikglcBTwjIpd1OLoAF8z8NeTugPTXuSzZGunzX9t0pE+llOqItvT2EeBFIMMY80RzZYwxA40xA4wxA4BlwB3GmOVeifC02dbNXyv/h4l9g0mMCeFdveFLKaU6pC01/ynATcAMEUl3T7NEZJ6IzPNxfO5hHx6BkkPIhue4LDmB1XuOklPSvudvKqWUOq7VsX2MMV8DbR6f1Rhza0cCataAs2HoBfDVn7nihit5eiV8sOkwP9LB3pRS6qR0vTt8W3Leb6CymEEZzzE6IVIf8qKUUh3QfZJ/r1GQfD1seJ4bhwubs4vYl1tqd1RKKdUtdZ/kDzD9l4Awu+AVRGC5XvhVSqmT0r2Sf1QinPETQjKWcfPAEl5du5/C8iq7o1JKqW6neyV/gKk/h+BIfhGwhJKKap5ascfuiJRSqtvpfsk/JAam3kN41hc8MDKH19Zmatu/Ukq1U/dL/gCTfgKRidx67G8EOYXff7TD7oiUUqpb6Z7J3xUMU3+G63Aaj43L51/bj7BuX57dUSmlVLfRPZM/QPKNEN6b/yh+nb5RwTz6z+3U1eljHpVSqi26b/J3BcOUnxKw/2t+f3oZWw8W65g/SinVRt03+QNMvBVCY5l6+FXGJUbx+Cc7OVZ1wlGllVJK0d2Tf2AYnDkf2fsZ/zO5hu+KK3jhq312R6WUUl1e907+AKf/JwRHMWrPC1w0ujeLvtxLTrGO+KmUUifS/ZN/cCScMQ92fMCvTjdU19bxp3/tsjsqpZTq0rp/8gcr+QeG03fLM9xy5gDeTMti+6GmT5pUSilV79RI/qE94PTbYOs7LBgHUSEuHvtwO8Zo10+llGrOqZH8Ac68E5zBRKQ9zd0zh7J6Tx4rd+bYHZVSSnVJp07yD4+3un5ueoMbRgiD4sJ47J8ZVNfW2R2ZUkp1OadO8gc46y5wBOBa+388MGske3PLWLLhgN1RKaVUl3NqJf+oBEi+Ab75O+cl1HDmoFie/Gw3Rceq7Y5MKaW6lFMr+QOcfTfU1SJrn+bBi0dSUF7FMyt1zH+llGro1Ev+MQNg7DWQ+jKjo6q4ckIiL6/OJCu/3O7IlFKqy2g1+YtIPxFZKSIZIrJNRBY0U+YGEdnsntaIyDjfhNtGU38ONRWwbiH3nj+cAIfw+491zH+llKrXlpp/DXCPMWYkMBmYLyKnNSnzLXCOMWYs8N/A894Ns53ihsKoy2HDC/R2lfOTcwbxz82HSdufb2tYSinVVbSa/I0xh40xG93LJUAGkNCkzBpjTIH75Tog0duBttvUe6CqFDY8z9xpg+gVGcQjH2TomP9KKUU72/xFZAAwHlh/gmI/Bj5q4f1zRSRVRFJzc3Pbs+v26z0ahl8M654ltK6ce88fzqasQv6x+ZBv96uUUt1Am5O/iIQDbwN3G2OaHThHRKZjJf/7mttujHneGJNijEmJj48/mXjbZ9o9UFEIqS9y5YRERvWN5H8/3klFtY75r5Tyb21K/iLiwkr8i40x77RQZizwV2C2MaZrPFA3YSIMnglrnsZRc4wHLx7JwcJjvLT6W7sjU0opW7Wlt48ALwIZxpgnWijTH3gHuMkY07XGU572Cyg/Chtf5azBcZw3shfPrNzL0dJKuyNTSinbtKXmPwW4CZghIunuaZaIzBORee4yvwZigWfc21N9FXC7JZ0JSWfD6v+DmkoemDWCiupa/vxp1zpHKaVUZ3K2VsAY8zUgrZS5DbjNW0F53bR74W+XQfpiBqf8iBsnJ/Ha2kxuOWsAw3pF2B2dUkp1ulPvDt/mDDoXElLg6z9DbTULZg4lPMjJY//MsDsypZSyhX8kfxGr7b/wAGx5i5iwQH46cyhf7srlq90+7nKqlFJdkH8kf4BhF0CvMfDVn6CulpvPHEBCdAhPfrZbn/illPI7/pP8Ray2/7w9sP09Ap0O5p0ziLT9Bazbp8M+KKX8i/8kf4CRl0LccFj1R6ir44cp/YiPCOLplbvtjkwppTqVfyV/h8Ma8ydnG+z6mGBXAD+ZNojVe/JI21/Q+vuVUuoU4V/JH2D0ldaY/6seB2O4/oz+xIS6WKgPfFFK+RH/S/4BTjj753BoI+z9nNBAJ7dNHcTnO3LYerDI7uiUUqpT+F/yBxh3HUQmWG3/wE1nJhER7NTav1LKb/hn8ncGwpS74cAayPyayGAXc84awEdbv2PXkRK7o1NKKZ/zz+QPMOEmCI2DtQsBmDNlIKGBAfqwd6WUX/Df5O8KgZQfwc6PIG8vMWGB3DQ5ifc3HSLzaJnd0SmllE/5b/IHOP3H4HDCBuuRwz+eOhBXgINnv9hrc2BKKeVb/p38I3pbXT+/+TtUFNEzIpjrJvXn7Y3ZHCw8Znd0SinlM/6d/AEmz7Me9P7N3wGYO20QIvDcl1r7V0qdujT59x0P/c+C9Yugrpa+0SFcNTGRN/6dRU5xhd3RKaWUT2jyB5h8uzXc884PAbj9nCHU1hle+GqfzYEppZRvaPIHGHExRPeHdc8C0D82lNnj+vL3dQfIL6uyOTillPI+Tf4AjgCYNBf2r4ZD6QDcMX0wFTW1vPT1tzYHp5RS3qfJv974m8AVZrX9A0N6RjBrdB9eXZNJ0bFqm4NTSinv0uRfLyQaxt8AW9+GkiMAzJ8+hJLKGl5bk2lvbEop5WWa/Bs6Yx7UVkHqSwCc1jeS80b25MXV31JWWWNzcEop5T2tJn8R6SciK0UkQ0S2iciCZsqIiDwlIntEZLOITPBNuD4WOxiGXQipL0K11c1z/vQhFJZXs3j9fpuDU0op72lLzb8GuMcYMxKYDMwXkdOalLkIGOqe5gLPejXKzjT5dijLtZp/gPH9Yzh7SBzPr/qWiupam4NTSinvaDX5G2MOG2M2updLgAwgoUmx2cBrxrIOiBaRPl6PtjMMPAd6nmZ1+zQGgDtnDOFoaSVL/51lc3BKKeUd7WrzF5EBwHhgfZNNCUDDzJjN908QiMhcEUkVkdTc3Nz2RdpZRKza/5EtkPk1AGcM7MHpA2JY9OVeqmrqbA5QKaU6rs3JX0TCgbeBu40xxU03N/MW870VxjxvjEkxxqTEx8e3L9LONOaHEBrruelLRLhzxlAOF1XwzsZsm4NTSqmOa1PyFxEXVuJfbIx5p5ki2UC/Bq8TgUMdD88mrhCYOMca7iHfGuJh2tA4xiZG8cwXe6mp1dq/Uqp7a0tvHwFeBDKMMU+0UOx94GZ3r5/JQJEx5rAX4+x8p99m3fm73hrrX0S4c/oQDuSX84/N3fe8ppRS0Laa/xTgJmCGiKS7p1kiMk9E5rnLfAjsA/YALwB3+CbcThTZB0Zd4R7r32rlOm9kL0b0juDpz/dQV/e9Vi2llOo22tLb52tjjBhjxhpjkt3Th8aYRcaYRe4yxhgz3xgz2BgzxhiT6vvQO8Hk26GqxDPWv8MhzJ8+hL25ZXy87Tubg1NKqZOnd/ieSMIE6DcZNjwHdVYf/1lj+jAoLoy/fL4HY7T2r5TqnjT5t2by7VCQCbs+BiDAIdwxfQgZh4tZuTPH3tiUUuokafJvzYhLIKqfp9snwOzkviTGhPDUCq39K6W6J03+rQlwWmP9Z34FhzcD4ApwcPu5g0nPKmTN3jybA1RKqfbT5N8WExqP9Q9w1cREekUG8ZfPd9sYmFJKnRxN/m0REgPJ18OWt6DUaucPcgbwk2mDWbcvn39n5tscoFJKtY8m/7Y64yeNxvoHuG5Sf+LCA3lqhdb+lVLdiyb/toobCkPPh3//FWoqAQgJDOC2qYP4avdR0rMKbQ5QKaXaTpN/ezQZ6x/gxslJRIe6+IvW/pVS3Ygm//YYNB3iR8C6Zzxj/YcHOfnxlIGs2JHD1oNFNgeolFJto8m/PerH+v9uC+xf7Vl981kDiAhysnDlHhuDU0qpttPk315jr4GQHo1u+ooKcXHrlAF8tPU7dh0psTE4pZRqG03+7eUKgZQ5sOOfkP+tZ/WPpgwkLDCApz/X2r9SquvT5H8y6sf63/CCZ1VMWCA3npnEB5sPsS+31MbglFKqdZr8T0ZkXxh1OWx8zTPWP8B/Th1EoNPBwpV7bQxOKaVap8n/ZJ3hHus//XXPqrjwIK6flMTy9INk5ZfbGJxSSp2YJv+TlTgREidZ4/24x/oHmDttEAEiPPOF1v6VUl2XJv+OOPMOKPgWtr3rWdU7KpirT09kWVoWhwqP2RicUkq1TJN/R4ycDT1HwcrHoLbas3reOYMxBp77Umv/SqmuSZN/RzgcMOMhyN/XqO0/MSaUKycksuTfWeQUV9gYoFJKNU+Tf0cNvwgSUuDLP0D18UR/x/TB1NTW8fyqfTYGp5RSzdPk31EiMPPXUHyw0XDPSbFhzE5OYPH6A+SVVtoYoFJKfV+ryV9EXhKRHBHZ2sL2KBH5h4hsEpFtIjLH+2F2cYPOgYHnwFd/gsrjwzvMnz6EippaXvz62xO8WSmlOl9bav6vABeeYPt8YLsxZhxwLvAnEQnseGjdzMxfQ/lRWHf8UY9DeoYza0wfXlu7n8LyKhuDU0qpxlpN/saYVcCJnlNogAgRESDcXbbGO+F1I4kpMHwWrPkLlB8/XHfNGEJpZQ0vr860LzallGrCG23+TwMjgUPAFmCBMaauuYIiMldEUkUkNTc31wu77mKmPwiVxbDmKc+qEb0jOf+0Xry8+ltKKqpP8GallOo83kj+FwDpQF8gGXhaRCKbK2iMed4Yk2KMSYmPj/fCrruY3qNhzFVW00/JEc/qu2YMpbiihtfW7rcxOKWUOs4byX8O8I6x7AG+BUZ44XO7p3MfsB70/tUfPavGJEZx7vB4/vrVPsoq/a9FTCnV9Xgj+R8AZgKISC9gOOC/ndtjB8OEmyD1ZSg4XtO/a8ZQCsqrWbxea/9KKfu1pavnEmAtMFxEskXkxyIyT0TmuYv8N3CWiGwBVgD3GWOO+i7kbmDaf4E4rBu/3CYmxTBlSCzPr/qWiuraE7xZKaV8ry29fa4zxvQxxriMMYnGmBeNMYuMMYvc2w8ZY843xowxxow2xvzd92F3cVEJ1gNfNi2B3J2e1XfNGMrR0kre2HDAxuCUUkrv8PWdqT8HVyis/B/PqsmDYpk0oAeLvtxHZY3W/pVS9tHk7ythcTD5Dti+HA6le1bfNXMI3xVXsCwt28bglFL+TpO/L511JwRHw+ePeladPSSOcf2iefaLvVTXNns7hFJK+Zwmf18KjoKzfwZ7PoX9awAQEX46YwjZBcd495uDNgeolPJXmvx9bdJcCO8FKx4BYwCYMaIno/pG8szKPdTWGZsDVEr5I03+vhYYCtN+AQfWwp4VgFX7v2vGEDLzyvlg8yGbA1RK+SNN/p1hwi0Q3R8+P177P/+03gzrFc7Tn++hTmv/SqlOpsm/MzgDrWEfDm+CjPcBcDiE+dOHsDunlI+3fWdzgEopf6PJv7OMvQbihls9f+qsPv6XjO3LoLgwnvxsl/b7V0p1Kk3+ncURADMehKO7YPNSAAIcwi9njWTXkVJ+9+EOmwNUSvkTTf6daeSl0CcZVv4Oaqzn+p53Wi9+NGUgr6zJ5KMth20OUCnlLzT5dyYRmPkrKDoAG1/zrL7/ohGMS4ziv5Zt5kBeuY0BKqX8hSb/zjZ4JiRNgS//F6rKAAh0Onj6+gmIwJ1LNmr7v1LK5zT5dzYRmPErKMuBDc97VvfrEcrjPxzH5uwibf9XSvmcJn87JJ0JQ34AXz8JFUWe1ReM6u1p//94q7b/K6V8R5O/XWY8BBWFsObpRqvr2/9/oe3/Sikf0uRvl77JcNplsHYhlOZ6Vte3/4O2/yulfEeTv52mPwg1x+DrJxqt7tcjlMevstr/f/+Rtv8rpbxPk7+d4odB8g2w/jnY+3mjTReO7s2cKQN4eXUmH2/V4R+UUt6lyd9uF/4O4ofDm7dC7q5Gmx64aKS7/X8TWfna/q+U8h5N/nYLioDr3oAAF7x+NZTnezY1av9/fSNVNfrkL6WUd2jy7wpikuDa16H4ILx5M9RUeTbVt/9vyi7idx9l2BikUupUosm/q+h/Blz6F8j8Cj681zPuP1jt/7eepe3/SinvaTX5i8hLIpIjIltPUOZcEUkXkW0i8qV3Q/Qj466Fs38OG1+F9YsabXpg1gjGJkbxX9r+r5TygrbU/F8BLmxpo4hEA88AlxpjRgE/9E5ofmrGr2DEJfDJL2H3p57VQc4Anr5uAgZt/1dKdVyryd8YswrIP0GR64F3jDEH3OVzvBSbf3I44PLnoNcoeGsO5Bxv5+8fG8rjV41lk/b/V0p1kDfa/IcBMSLyhYikicjNLRUUkbkikioiqbm5uS0VU0HhVg+gwFB4/RooO+rZdOHoPtx61gBeWv0tn+jjH5VSJ8kbyd8JTAQuBi4AfiUiw5oraIx53hiTYoxJiY+P98KuT2FRiXDtEig9Aktv9Dz8Baz2/zEJUfziLW3/V0qdHG8k/2zgY2NMmTHmKLAKGOeFz1WJE2H2QjiwFj74uacHUJAzgIXXT8AYuHPJN9r+r5RqN28k//eAqSLiFJFQ4AxAO6R7y5ir4Jz7IP3vsOYvntX9Y0P536vGsimrkD98rO3/Sqn2cbZWQESWAOcCcSKSDfwGcAEYYxYZYzJE5DcUMCAAABIWSURBVGNgM1AH/NUY02K3UHUSzrkfcnfCp7+GuKEw/CIALhpjtf+/+PW3jOobyRUTEm0OVCnVXYhpcDNRZ0pJSTGpqam27LtbqiqHV2bB0d3wo0+g92gAKmtquf6F9aTtL+DisX14+D9GER8RZHOwSilfEZE0Y0xKRz9H7/DtLgJDrQvAQRGw5FootXrUBjkDWPKfk7n3/GF8uu0I5z3xJW+mZmHXSV0p1T1o8u9OIvvAdUusrp9v3ADVFYA1ANydM4by4YKpDO8VwX8t28yNL65nf16ZzQErpboqTf7dTd/xcPkiyN4A//hpozGAhvQM5425k3n0stFsyirigidX8dyXe6mp1d5ASqnGNPl3R6Mug+kPweal33sKmMMh3Dg5ic9+fg5Th8bzu492MHvharYeLGrhw5RS/kiTf3c17V4Y80NY8Qhsf/97m3tHBfP8TRN59oYJ5JRUMnvhan73UQbHqvSZwEopTf7dlwhc+jQkpMC7P4F9XzZqArKKCBeN6cNnPzuHqyYk8tyX+7jw/1axZs/RFj5UKeUvtKtnd1dyBP56HhQdgLjh1rDQY6+BqITvFV2z9yi/fGcLmXnlXJ2SyIOzTiMq1GVD0Eqd4urqoLbKPVVDbaV7uQYw7oqaAVPXYLm1OYBB+p3ula6emvxPBRXFsO1d2LTEGgoCgUHnwLjrYeQlEBh2vGh1LU9+tpsXvtpHTGggv710FLPG9EZE7ItfKbtVV0B5XpMpv/HrYwXWGFu1Ve5kXm0t11Q1SPTuqa7GZ6HKb4s1+atm5O+DTUutE0HhfggMh9Nmw7jrIGmKNWQ0sPVgEfe/s5mtB4s5b2Qv/vuyUfSJCrE5eHXKqSiGslx3kqw8nhzrl2sqj9eMG22varCuGkwtSACIw2rydNQv108NXzfdHmCtqyxtOcFXn6BbdHA0hMZCSAy4QiAg0D25wBl0fDkgyD0PdK93LzdcH+ACxIpH5PhyO+Yy/AJN/uoE6uogax2kvw7blkNVCUT1h3HXWCeC2MHU1Nbx8upM/vTpTgJEGNEnkqgQl2eKDHERGexstC4q1L0t2EVoYID1i6G2xnr+cNlR60sfGAou9xQYZv0hqlOLMVZNuPggFB9yzw83WD5kTVUlHdiJHE+u4rD2aercU+3x5bpa3G0irQuKhNAeVjJvNDW3LtZK/AGtjoLTqbx1h68mf39QVQ47P7ROBPtWWn8wiZMg+ToYdTkHyoN46vPdHCo8RtGxas9UUmH9dA2iigQ5+r2pnxwl0XGUnuQTQMv3EtQ6AqlzhmJcoeAKQQLDcARZk3hOEu55UAQERx2fgiK//9pxCvdTMAaqy61aanV5M7Xkyga14vp5RTPr3DVmAIfTOgE7nFbNs+Frz9T0tfP4Sbs093hCLzl8fLmmonHs4oDw3hDZ1z0lWPPwnu4k7k7kzga14fr1zsDjNer6hO9wumu7bTxunpNDgxND/cnB1Fm/gp2B3vu/sokmf3Vyig/DljchfQnkZlh/eMMvsrqNigOKsqDwABRlYQqtyVHe+ME7dRJAWVBPCly9yXP25DvpyUETx3c1EdRUV2CqyqG6nMC6CkKpJESqCKHhciWhVBIqlYQ5qgiTSkKoJJgKgk1lC4HXkyYnhCYnh/oTRFCEtS0oAoKi3HP3OleY904gdbVQVWZN1eUN5qVWAvfMS5q8LoXKkuZfmw7elCcBx5sdECvGuhr3VH1yn+lwWXeY1yf0hsk9MgEi+kB4ry5XSz4VafJXHWMMHE6HTW/Alresds96zmDrYTJR/SC6n9VcFN3v+OuIvm36I6+orqW4wvoFUXysmuKKGkoqqik+VuNe33C5hqJj1ew9UkRdRTFxARVM6efi7EQnKb0DiA04ZrUfVxSdeKospvUmAGnmBBFxfF1QhFVzbpTM3Qm9frm6zJrXtnayasDhcn9+OATWz8MbvG6yzhV2vJZcXxtuNA9qsL1BuRM1s9XXjj0ng5omJ4cm60wdhMVDaNyp/YurG9Hkr7ynpgqy1lvNLtH9rD92m3r/VNfWkZpZwIqMI6zYkcO3R60LcSN6R3DeyF7MHNmTcYnROBwtxFdX565ll1gni0r3cmVRk3XueUVR49eV7hq6M9BKvoENrl14mqfCrNee5WbKBDZN7OFWYlaqgzT5K7+wN7eUFRlH+Cwjh7T9BdTWGeLCA5k+vCczR/Zi6tA4woK0qUH5D03+yu8Ullfx5a5cPsvI4YudOZRU1BDodHDmoFjOG9mTGSN7kRCt3VXVqU2Tv/Jr1bV1/DsznxUZOazIOEJmnvUg+xG9I+jfI5Q6Y6ipM9Q2nYw1r6k1njJ1dY3LGgyRwS56hAUSGx5IbFiQZ7lHmDXFhVvrYkIDCWipCUopH9Dkr1QD9c1DK3fkUlBeRYBDcDoER/1cBGeAe+4QAhpNDgIEa+4AQSg6Vk1+WRV5ZZXkl1VRUN58LxkRiA6pP1EEEes+OfSODGZcv2iS+0cTGaxDaCjv8Vby18ZSdUoYHB/O4Phw5k4b7JPPr6mto6D8+Akhr7TKvVxFvvt1XlkVu3NKyS+ztoF1chjeK4IJSTGkJMUwMSmG/j1CdTgNZTtN/kq1gTPAQXxEkPv5yBGtli+pqGZTVhGp+/NJ21/AP9IP8fr6AwDEhQcyob91IpiYFMPohCiCXd69C9oYQ3lVLc4AIcipd1ir79Pkr5QPRAS7OHtoHGcPjQOgts6wO6eEtP0FpO0vYOP+Av61/QgArgBhdEIUExucEHpGBjf6vNo6Q2F5ledXRX5ZFfnlVeS7f3EUuLfllVrLeWVVVNVYN4uFBgYQHeIiKjSQqBAn0SGBRIdaQ3XUL0eHNHkd6iLEFaC/UE5h2uavlE2OllaycX8BaQesk8Gm7CJPwk6MCaFPVLAn0Rceq276uAaPiCAnMe5rDbFhgcQ0mNfU1lFYXk1h/bAd5dUUHquy1pVXU3WCR3y6AoSokEDiwq0L3LFN5sfXW9c6vP3rxVsqqmvJLijH6XDQIzyQiCBntz6pdVqbv4i8BFwC5BhjRp+g3OnAOuAaY8yyjgam1KkuLjyI80f15vxRvQGoqqlj26Ei65fBgQLySqsY3jvC6mEU6u5pFB7kWY4Nt2rpJ9usY4yhorqu0cmgqH75WLV7nfUr4mhpJQcOlHO0tJLyFp4GFxHkJC7COhE0OklEBBEfHmg1m4UHExcRSGigdxsdiiuqOZBXTmZeGfvzytnvmZfzXXHjMYgCAxzEhLnoEWbFWt+Ly7pY37hnV2xYIJHBrpZvKuzGWq35i8g0oBR4raXkLyIBwKdABfBSW5K/1vyV6p7Kq2rIK7VOCEdLq8grrSSvrIrcEmt+tKSSvDJrW0F5VbO/WMICAzzXUOLCg9wnhvoTRZBnW2x4IEHOAIwx5JVVeZJ6Zl45B/LKrHl+uecCe7248CCSYkOtqUcY/WNDqK3DujhfZjWXHb9gb02llc2PwR/gEGJCrZPX6L6RTEyKYUJSDEPiw205KXRazd8Ys0pEBrRS7C7gbeD0jgaklOraQgOdhPZw0q9HaKtla2rryC+rItd9osgtqfRMR0ut+e6cUtbuy6Owhe60USEuautMo+QsAn2jQkiKDeWCUb1Jig1lQGwo/XuE0T82lPCTuOu7orrWul5SevyEUN+bK7+sisNFFXyWcYS30rIBiAh2Mr5/DBP7xzAhKZrkftFEdKNuvR3+7SUiCcDlwAxaSf4iMheYC9C/f/+O7lop1cU5Axz0jAz+3gXs5lTW1JLnPkHUnxhySyrJLa3EIeJO8FZyT4wJ8XovpmBXAH2iQk74UCNjDJl55Z6muY37C3hyxS6MOd6td7z7wv2E/tEMjAvrstcX2nTB113z/6C5Zh8ReQv4kzFmnYi84i6nzT5KKb9QXFHNpqxCNu4vJO1AAd8cKPA8C6NHWCDj+0UzISmGCf1jGJ0QSXgHLzh3pZu8UoA33P+YOGCWiNQYY5Z74bOVUqpLiwx2MXVoPFOHxgNQV2fYk1tq9eRy/0JYsSPHUz4wwEF0qIuY0ECiQl3EuJejQwPd613WcoiLmLD6rriBBDq9O6R2h5O/MWZg/XKDmr8mfqWUX3I4hGG9IhjWK4JrJ1nN2wVlVaRnFbLzSAkF5VUUllldbgvKq/n2aBkbywspLK+iurbllpiwwACiQ733JLK2dPVcApwLxIlINvAbwAVgjFnktUiUUuoUFRMWyPQRPZk+omeLZervyi4oP971tqDcusej0D2+VGF5FWu8FJPe5KWUUt2It9r89blsSinlhzT5K6WUH9Lkr5RSfkiTv1JK+SFN/kop5Yc0+SullB/S5K+UUn5Ik79SSvkh227yEpESYKctO2+fOOCo3UG0gcbpXd0hzu4QI2ic3jbcGNP6g6RbYeczfHd64y41XxORVI3TezRO7+kOMYLG6W0i4pWhEbTZRyml/JAmf6WU8kN2Jv/nbdx3e2ic3qVxek93iBE0Tm/zSpy2XfBVSillH232UUopP6TJXyml/JDPk7+IXCgiO0Vkj4jc38z2IBFZ6t6+3v2w+E4lIv1EZKWIZIjINhFZ0EyZc0WkSETS3dOvOztOdxyZIrLFHcP3unyJ5Sn38dwsIhNsiHF4g+OULiLFInJ3kzK2HE8ReUlEckRka4N1PUTkUxHZ7Z7HtPDeW9xldovILZ0c4+MissP9f/quiES38N4Tfj86Ic6HReRgg//XWS2894R5oRPiXNogxkwRSW/hvZ15PJvNQz77fhpjfDYBAcBeYBAQCGwCTmtS5g5gkXv5WmCpL2NqIc4+wAT3cgSwq5k4z8V6PnGnxtZMrJlA3Am2zwI+AgSYDKy3Od4A4DsgqSscT2AaMAHY2mDd/wL3u5fvB/7QzPt6APvc8xj3ckwnxng+4HQv/6G5GNvy/eiEOB8G7m3Dd+KEecHXcTbZ/ifg113geDabh3z1/fR1zX8SsMcYs88YUwW8AcxuUmY28Kp7eRkwU0TEx3E1Yow5bIzZ6F4uATKAhM6MwYtmA68ZyzogWkT62BjPTGCvMWa/jTF4GGNWAflNVjf8Dr4KXNbMWy8APjXG5BtjCoBPgQs7K0ZjzL+MMTXul+uARF/suz1aOJZt0Za84DUnitOda64Glvhq/211gjzkk++nr5N/ApDV4HU230+qnjLuL3cREOvjuFrkbnYaD6xvZvOZIrJJRD4SkVGdGthxBviXiKSJyNxmtrflmHema2n5D6srHE+AXsaYw2D9AQLNPWW7Kx3XH2H9umtOa9+PznCnu3nqpRaaKLrSsZwKHDHG7G5huy3Hs0ke8sn309fJv7kafNO+pW0p0ylEJBx4G7jbGFPcZPNGrKaLccBfgOWdHZ/bFGPMBOAiYL6ITGuyvSsdz0DgUuCtZjZ3lePZVl3iuIrIg0ANsLiFIq19P3ztWWAwkAwcxmpSaapLHEu36zhxrb/Tj2creajFtzWz7oTH1NfJPxvo1+B1InCopTIi4gSiOLmfkh0iIi6sA77YGPNO0+3GmGJjTKl7+UPAJSJxnRwmxphD7nkO8C7WT+iG2nLMO8tFwEZjzJGmG7rK8XQ7Ut805p7nNFPG9uPqvoh3CXCDcTf0NtWG74dPGWOOGGNqjTF1wAst7N/2YwmefHMFsLSlMp19PFvIQz75fvo6+f8bGCoiA921wGuB95uUeR+ovzJ9FfB5S19sX3G3+70IZBhjnmihTO/6axEiMgnr2OV1XpQgImEiElG/jHURcGuTYu8DN4tlMlBU/5PRBi3WqrrC8Wyg4XfwFuC9Zsp8ApwvIjHupozz3es6hYhcCNwHXGqMKW+hTFu+Hz7V5PrS5S3svy15oTOcB+wwxmQ3t7Gzj+cJ8pBvvp+dcAV7FtZV673Ag+51j2B9iQGCsZoF9gAbgEG+jqmZGM/G+om0GUh3T7OAecA8d5k7gW1YPRPWAWfZEOcg9/43uWOpP54N4xRgoft4bwFSOjtOdxyhWMk8qsE6248n1snoMFCNVVv6MdY1phXAbve8h7tsCvDXBu/9kft7ugeY08kx7sFq063/ftb3kOsLfHii70cnx/k39/duM1bS6tM0Tvfr7+WFzozTvf6V+u9jg7J2Hs+W8pBPvp86vINSSvkhvcNXKaX8kCZ/pZTyQ5r8lVLKD2nyV0opP6TJXyml/JAmf6WU8kOa/JVSyg/9Pxda6Mgu94+jAAAAAElFTkSuQmCC\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.\n",
    "\n",
    "However, we need to define at how many points we want to get the predictions.\n",
    "The default (`model.sub = 1`) is just to use the interval knots, but by increasing the `model.sub` argument, we replace the knots with and equidistant number of points in each interval.\n",
    "This is very similar to the interpolation of the discrete methods such as `LogisticHazard` and `PMF`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "surv = model.predict_surv_df(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAd00lEQVR4nO3de3SV9Z3v8feXhBAkIEhAlIBAlXJRihBBp+qg1SlGBmYqN++iHTvnSKc91rOq9RyrzvJSO56WrrpaqXc9lTp2HDlKsVpl2uVSJCregCiDIJuLBFABIYQk3/PHfqJ7shNy209+e+98XmtlsZ9nP3v/PjwJfPN7nt/+/czdERERSdUjdAAREck+Kg4iIpJGxUFERNKoOIiISBoVBxERSVMYOkBHlJaW+ogRI0LHEBHJGa+//vpOdx/U1uNzsjiMGDGCysrK0DFERHKGmW1qz/G6rCQiImlUHEREJI2Kg4iIpMnJew4iIiEcOnSIRCJBTU1N6CgtKi4upqysjJ49e3bqfVQcRETaKJFI0LdvX0aMGIGZhY6Txt3ZtWsXiUSCkSNHduq9Yr2sZGYPmNkOM3u3hefNzH5hZuvN7G0zmxRnHhGRzqipqWHgwIFZWRgAzIyBAwdmpGcT9z2Hh4Dph3n+POCE6Otq4Fcx5xER6ZRsLQyNMpUv1stK7v5nMxtxmENmAY94ct7wV82sv5kd4+7bDve+OzdWc//l92YwafsV9N/EFYtuD5pBRCQuoUcrDQU2p2wnon1pzOxqM6s0s0oCL0FxqOdQ6j89LmwIEem2li9fzle/+lWOP/547rzzzljaCH1Durn+T7P/9bv7YmAxQHl5uV/18HfizHVYoXstItJ91dfXc8011/D8889TVlbGKaecwsyZMxk3blxG2wndc0gAw1K2y4CtgbKIiGS91157jeOPP55Ro0ZRVFTE/PnzefrppzPeTuiew1JgoZktAaYCn7V2v0FEJBvc8v/eY83WPRl9z3HH9uPHfzv+sMds2bKFYcO+/J26rKyMlStXZjQHxFwczOxxYBpQamYJ4MdATwB3/zWwDKgA1gP7gQVx5hERyXXJ8Tv/VRwjqOIerXRhK887cE2cGURE4tDab/hxKSsrY/PmL8fxJBIJjj322Iy3E/qeg4iItMMpp5zCBx98wIcffkhtbS1Llixh5syZGW8n9D0HERFph8LCQn75y1/yzW9+k/r6eq688krGj898L0bFQUQkx1RUVFBRURFrGyoOHVRU62y69LLQMeg3YwYD5s0NHUNE8ozuOXTA572htij8/Co169ax55lnQscQkTyknkMH7C0x9pbAcf/ySNAc2dBzEZH8pJ6DiIikUXEQEZE0uqzUQTXmLFge9gPd83evY2DvgWh+WBHJNPUcOuDI+h4Ue/gb0vvrDrDrwK7QMUSkC1155ZUMHjyYE088MdZ21HPogP4NRmHNUHjh8qA5dg75K0o+Xx00g4h0rSuuuIKFCxdy2WXxDkhRz6EDvnLUWkqKtoeOwaGiY9nXZ2LoGCLShc4880yOOuqo2NtRz6EDThu+ltOGr4UFzwbN8ZvLfxO0fZFu7Q/Xw/Z3MvueQ06C8+JZ2a291HMQEZE06jmIiHRElvyGHxf1HEREJI2Kg4hIDrnwwgs57bTTqKqqoqysjPvvvz+WdnRZqaO2vwMPnh80Qg9m4oT/vIWIdJ3HH3+8S9pRceiIk2aHThBRaRCReKg4dET5guRXaCvuDZ1ARPKU7jmIiEgaFQcREUmj4iAiImlUHEREJI2Kg4hIDtm8eTNnnXUWY8eOZfz48SxatCiWdjRaSUQkhxQWFnL33XczadIk9u7dy+TJkzn33HMZN25cRttRz0FEJIccc8wxTJo0CYC+ffsyduxYtmzZkvF21HMQEemAn7z2E9btXpfR9xxz1Bh+OOWHbT5+48aNvPnmm0ydOjWjOUA9BxGRnLRv3z4uuOACfv7zn9OvX7+Mv796DiIiHdCe3/Az7dChQ1xwwQVcfPHFfOtb34qlDfUcRERyiLtz1VVXMXbsWK699trY2om152Bm04FFQAFwn7vf2eT54cDDQP/omOvdfVmcmfJNA86C5WHneaoYVcGc0XOCZhDpLl5++WUeffRRTjrpJCZOTK4hf/vtt1NRUZHRdmIrDmZWANwDnAskgFVmttTd16Qc9r+AJ9z9V2Y2DlgGjIgrU74pcMDCzstatbsKQMVBpIucfvrpuHvs7cTZc5gCrHf3DQBmtgSYBaQWBwca76QcCWyNMU/e6enJrwenPxgsQ+hei4jEI857DkOBzSnbiWhfqpuBS8wsQbLX8N2W3szMrjazSjOrrK6uznRWERFJEWdxaO56R9O+0IXAQ+5eBlQAj5pZs5ncfbG7l7t7+aBBgzIcVUREUsVZHBLAsJTtMtIvG10FPAHg7q8AxUBpjJlERKQN4iwOq4ATzGykmRUB84GlTY75CPgGgJmNJVkcdM1IRCSw2IqDu9cBC4HngLUkRyW9Z2a3mtnM6LAfAP9gZm8BjwNXeFfchhcRkcOK9XMO0WcWljXZd1PK4zXA1+PMkPfqYNOllwVrfv7udaydXArTg0UQ6VZqamo488wzOXjwIHV1dcyePZtbbrkl4+1o+owc5r17YAcagmYYvGU/sDNoBpHupFevXrz44ouUlJRw6NAhTj/9dM477zxOPfXUjLaj4pDLevfAe/fguEWPBIuw7vwpwdoW6Y7MjJKSEiA5x9KhQ4ewGD4Mq+IgItIB22+/nYNrMztld6+xYxjyox+1elx9fT2TJ09m/fr1XHPNNZqyW0REoKCggNWrV5NIJHjttdd49913M96Geg4iIh3Qlt/w49a/f3+mTZvG8uXLOfHEEzP63uo5iIjkkOrqaj799FMADhw4wAsvvMCYMWMy3o56DiIiOWTbtm1cfvnl1NfX09DQwNy5c5kxY0bG21FxyHH7aofw1N1vBGt/96Cr6bt/dbD2RbqbCRMm8Oabb8bejopDDju6z9vRo2OCZagtOoa9wVoXkbioOOSwY0te59iS1xn/g6uCZbj/ivuCtS0i8dENaRERSaPiICIiaXRZSTqtweuDLxdaMapC61iLZJB6DtIpPQsK6WEFQTNU7a5i2YZlrR8oIm2mnoN0Ss8eRfTsUcSD0x8MliF0r0Wkq9XX11NeXs7QoUN55plnYmlDPQcRkRyzaNEixo4dG2sbKg4iIjkkkUjw7LPP8u1vfzvWdnRZSUSkA/7yxPvs3Lwvo+9ZOqyEM+aOPuwx3//+97nrrrvYuzfej5+q5yAikiOeeeYZBg8ezOTJk2NvSz2HHDfi0AZ48PxwAWpnQUFRuPZFAmntN/w4vPzyyyxdupRly5ZRU1PDnj17uOSSS3jssccy3pZ6Djns5d5nsbHnqLAhGhqgvjZsBpFu4o477iCRSLBx40aWLFnC2WefHUthAPUcctqfjqjgT0dU8LsFp4ULsWJxuLZFJDYqDiIiOWjatGlMmzYttvfXZSUREUmjnoN0WkOts+nSy4K1P3/3OtZOLoXpwSKI5B0VB+mUwiMKqKMeGsJlGLxlP7AzXADpVtwdMwsdo0XunpH3UXGQTiksKaSwpJDjfvpIsAzrzp8SrG3pXoqLi9m1axcDBw7MygLh7uzatYvi4uJOv5eKg4hIG5WVlZFIJKiurg4dpUXFxcWUlZV1+n1UHERE2qhnz56MHDkydIwuoeKQ49Zs28O8e18J1v7f1tbRs0CD3kTyTazFwcymA4uAAuA+d7+zmWPmAjcDDrzl7hfFmSmfzJo4NHQEGhqcQyHvRkf21x3IinUdtCKd5IvYioOZFQD3AOcCCWCVmS119zUpx5wA3AB83d0/MbPBceXJRxdNHc5FU4cHzfDY99a0flDMBvYeCAd2hY5B1e4qABUHyQtx9hymAOvdfQOAmS0BZgGp/5v8A3CPu38C4O47YswjeWpQ70EM6j0o6Gp0oBXpJL/EebF4KLA5ZTsR7Us1GhhtZi+b2avRZahmmdnVZlZpZpXZPFJARCQfxNlzaG4QcNNPZxQCJwDTgDLgL2Z2ort/mvZC98XAYoDy8vLMfMpDMqLYa8JOG759K/QZFK59kTwUZ3FIAMNStsuArc0c86q7HwI+NLMqksViVYy5JIP29OgPDWm1vGvVfh62fZE8FGdxWAWcYGYjgS3AfKDpSKR/By4EHjKzUpKXmTbEmEky7JOCgXxSMBAWPBsuxG9PDte2SJ6K7Z6Du9cBC4HngLXAE+7+npndamYzo8OeA3aZ2RrgJeB/unv4YSciIt1crJ9zcPdlwLIm+25KeezAtdGXiIhkiTYVBzMrB84AjgUOAO8CL7j77hiziYhIIIe9rGRmV5jZGyQ/qNYbqAJ2AKcDz5vZw2YW9lNYIiKSca31HPqQ/PTygeaeNLOJJEcXfZTpYCIiEs5hi4O739PSc2ZW5O6rMx9JRERCa+s9hxXAFe6+MdqeAvwG+FpsyUTaoWZHbdClSkHLlUp+aetopTuA5Wb2C5JTYJwHaCIZyQr9xpVA7cew/Z2gOQbvqGH/oc+Dz7GkmWElE9pUHNz9OTP7R+B5kov1nuzu22NNJtJGAy69kgETnwwdg+pH13OEh13bQjPDSqa09bLS/wbmAmcCE4AVZvYDdw/4sVjJFiUHGnjq7jcCJvgao6f8DePPCLu+xaDfnswgCDo7bOhei+SPtl5WKgWmRKOWXjGz5cB9gIpDN7djQEHoCOxM7AMIXhxE8klbLyt9r8n2JpKL+Eg3t720kO2lhdz4nUnBMoTttYjkJy3+KyIiaVQcREQkTawT74l0O7Wfh134yD7WwkeSEYctDmZ2ZhvfZ6O7awoN6d6y4T9lLXwkGdJaz6Gt4+KeQvMrSXfXd0jya8Ej4TI8VB6ubckrrc2tpEHTIiLdkG5Ii4hIGhUHERFJo+IgIiJp2lQczOzRtuwTEZH80NbPOYxP3TCzAmBy5uOI5LaadeuCrisxf/sBXhpvWTEBn6YOz22trSF9g5ntBSaY2Z7oay/JdaSf7pKEIjmi34wZFI8ZEzTD0B3OWe950AyQnDp82YZloWNIJ7Q2lPUO4A4zu8Pdb+iiTCI5acC8uQyYNzdohk3nnswYwk4bDpo6PB+01nMYAdBSYbCksszHEhGRkFq75/BTM+tB8hLS60A1UAwcD5wFfAP4MZCIM6SIiHSt1i4rzTGzccDFwJXAMcABYC3JhX5uc/ea2FOKiEiXanW0kruvAW7sgiwikgmhZ4YFzQ6bB1q753CKmQ1J2b7MzJ42s1+Y2VHxxxORdukzCIr6hE6RLFCfV4dOIZ3QWs/hXuAc+GL67juB7wITgcXA7FjTSU5Ys20P8+59JVj7E7YdpLSkKFj7WSUbZoYFzQ6bB1orDgXuvjt6PA9Y7O6/B35vZqvjjSa5YNbEoaEjsP9gHTtDhxDJM60WBzMrdPc6kiOTrm7Ha6UbuGjqcC6aOjxohtuuWxG0fZF81NrcSo8D/2FmT5McpfQXADM7HvistTc3s+lmVmVm683s+sMcN9vM3MzUFxURyQKtDWW9zcz+RHII6x/dvfFz+T1I3ntoUTT/0j3AuSQ/B7HKzJZGo59Sj+sL/BOwsmN/BRERybS2DGV9tZl977fhvacA6919A4CZLQFmAWuaHPfPwF3AdW14TxHJFdkwpPak2VCuqTw6Is71HIYCm1O2E9G+L5jZycAwd3+mtTczs6vNrNLMKqurNUROJKtlw5Da7e/AO0+GzZDD4rypbM3s+2K6yGhajp8BV7Tlzdx9Mcnhs5SXl4efdlJEWtY4pDbkBIChey05Ls7ikACGpWyXAVtTtvsCJwIrzAxgCLDUzGa6e2WMuSQPlRxo4Km73wgdg9FTjmb8GeGH94p0VpzFYRVwgpmNBLYA84GLGp9098+A0sZtM1sBXKfCIO21Y0BB6AgA7EzsA1BxkLwQW3Fw9zozWwg8BxQAD7j7e2Z2K1Dp7kvjalu6l+2lhWwvLeTG70wKmiMbei4imRLrB9ncfRmwrMm+m1o4dlqcWUREpO3iHK0kIiI5SsVBRETSaH4kEYlF1e6qsGtJ28dUeB/mhEuQ01QcRCTjKkZVhI5AFbVgqDh0kIqDiGTcnNFzmDM67H/LC7SmRKfonoOIiKRRz0Ekz9SsW8emSy8LHYN+M2YwYN7c0DGkg1QcRPJIvxkzQkcAkgUKUHHIYSoOInlkwLy5WfEfcjb0XKRzdM9BRETSqDiIiEgaFQcREUmj4iAiIml0Q1rywppte5h37ytBM0zYdpDSkqKgGUQyRcVBct6sidmxuM7+g3XsDB1CJENUHCTnXTR1OBdNHR46BrddtyJ0BJGM0T0HERFJo56DiOStKmrDThseqRhVEXwiwvZScRCRvFThfcBCp0iuawGoOIiIZIM5lDDHS2D6g0FzZEPPpSN0z0FERNKoOIiISBoVBxERSaN7DiISi+CLDm3fSr9xJQwIlyCnqTiISMZlw6JDNTtqgX0qDh2k4iAiGZcNiw5tOvfkoO3nOt1zEBGRNCoOIiKSRsVBRETSqDiIiEiaWG9Im9l0YBFQANzn7nc2ef5a4NtAHVANXOnum+LMJCLdSO3n8OD5YTPYx9BnUNgMHRBbz8HMCoB7gPOAccCFZjauyWFvAuXuPgF4Ergrrjwi0s30GQRFfUKnSBaoz6tDp2i3OHsOU4D17r4BwMyWALOANY0HuPtLKce/ClwSYx6R2JUcaOCpu98ImmH0lKMZf0Z2rI4XVN8hya8Fj4TN8VB52PY7KM57DkOBzSnbiWhfS64C/tDSk2Z2tZlVmllldXXuVWHJfzsGFLCvd9jbeDsT+3j/tY+DZpD8EGfPobmZ1L3ZA80uAcqBv27pzdx9MbAYoLy8vNn3EQlpe2kh20sLufE7k4JlCN1rkfwRZ3FIAMNStsuArU0PMrNzgBuBv3b3gzHmERGRNoqzD7wKOMHMRppZETAfWJp6gJmdDNwLzHT3HTFmERGRdoit5+DudWa2EHiO5FDWB9z9PTO7Fah096XAT4ES4F/NDOAjd58ZVyaRuK3Ztod5974SrP0J2w5SWlIUrH3JH7F+zsHdlwHLmuy7KeXxOXG2L9KVZk0MP0Jo/8E6doYOIWmqqM255UI1K6tIhlw0dTgXTR0eNMNt160I2r6kq/A+zQ/P6UJVu6va/RoVBxGRGM2hhDleAtMfDJZhwfIFvEL7LndqbiUREUmj4iAiImlUHEREJI2Kg4iIpFFxEBGRNCoOIiKSRsVBRETSqDiIiEgafQhORPJWzbp1bLr0srAhtm+lX9lnDCDgcqXW/jU+VBxEJC/1mzEjdAQAanY5cCQDQoao/bzdL1FxEMkz+w/WBZ0ZttGsiUODzjU1YN5cBsybG6z9Rl/0XEIuV9qBpUp1z0Ekj5SWFHFEr/C/863ZtoenV28JHUM6IfxPkYhkzOC+xQzuWxx0qVIgK3ou0jnqOYiISBoVBxERSaPiICIiaVQcREQkjYqDiIik0WglkTyzM7GPp+5+I2iGCdsOsmNAQdAMkqKoT7tfouIgkkdGTzk6dAQASg40hI4gqY4aBfy5XS9RcRDJI+PPGMr4M4aGjsFt160IHUE6ScVBRCRmoScAnL97HQ+18zUqDiIiMcqGCQAHb9nf7teoOIiIxCgbJgD86O+mwVvte42Kg4jEIhtmhw09M2y2GN63/edAn3MQkYzLhtlhNTNs56jnICIZlw2zw4buteQ69RxERPJcr7Fj2v0a9RxEJG+t2bYnK3oQoe99DPnRj+DGG9v1mliLg5lNBxYBBcB97n5nk+d7AY8Ak4FdwDx33xhnJhHpGqGn8Th7rzOk5xFsD5YgaeWHu1n54e6cu/8RW3EwswLgHuBcIAGsMrOl7r4m5bCrgE/c/Xgzmw/8BJgXVyYR6RrZMI1Hj8/qmFZWwt8HXhXvtys/Cl4YVn64u92vibPnMAVY7+4bAMxsCTALSC0Os4Cbo8dPAr80M3N3jzGXiMQsG6bxeOruN4L3XgB6A/PpFTTDaT368T/a+Zo4b0gPBTanbCeifc0e4+51wGfAwObezMyuNrNKM6usrq6OIa6I5JPRU46mtKwkdIyscNzA7JqV1ZrZ17RH0JZjkjvdFwOLAcrLy9WzEJHDyobeS1a5rn2Hx9lzSADDUrbLgK0tHWNmhcCRQPsvjomISEbFWRxWASeY2UgzKwLmA0ubHLMUuDx6PBt4UfcbRETCi+2ykrvXmdlC4DmSQ1kfcPf3zOxWoNLdlwL3A4+a2XqSPYb5ceUREZG2i/VzDu6+DFjWZN9NKY9rgDlxZhARkfbT9BkiIpJGxUFERNKoOIiISBoVBxERSWO5OHLUzPYCVaFzdFApsDN0iA5S9jCUPYx8y36cuw9q6xvk6pTdVe5eHjpER5hZpbJ3PWUPQ9nDyER2XVYSEZE0Kg4iIpImV4vD4tABOkHZw1D2MJQ9jE5nz8kb0iIiEq9c7TmIiEiMVBxERCRNThUHM5tuZlVmtt7Mrg+dpzVmttHM3jGz1WZWGe07ysyeN7MPoj8HhM7ZyMweMLMdZvZuyr5m81rSL6LvxdtmFnSh3hay32xmW6Lzv9rMKlKeuyHKXmVm3wyTGsxsmJm9ZGZrzew9M/tetD/rz/thsmf9eY+yFJvZa2b2VpT/lmj/SDNbGZ3730VLDmBmvaLt9dHzI7Iw+0Nm9mHKuZ8Y7W//z42758QXyWm//xMYBRQBbwHjQudqJfNGoLTJvruA66PH1wM/CZ0zJduZwCTg3dbyAhXAH0iu5ncqsDILs98MXNfMseOin59ewMjo56ogUO5jgEnR477A+1G+rD/vh8me9ec9ymNASfS4J7AyOqdPAPOj/b8G/lv0+L8Dv44ezwd+l4XZHwJmN3N8u39ucqnnMAVY7+4b3L0WWALMCpypI2YBD0ePHwb+LmCW/8Ld/0z6Snwt5Z0FPOJJrwL9zeyYrkmaroXsLZkFLHH3g+7+IbCe5M9Xl3P3be7+RvR4L7CW5NrqWX/eD5O9JVlz3gGic7gv2uwZfTlwNvBktL/puW/8njwJfMPMmlvqOHaHyd6Sdv/c5FJxGApsTtlOcPgfxGzgwB/N7HUzuzrad7S7b4PkPy5gcLB0bdNS3lz5fiyMutEPpFzCy8rs0WWKk0n+FphT571JdsiR825mBWa2GtgBPE+yN/Opu9dFh6Rm/CJ/9PxnwMCuTfylptndvfHc3xad+5+ZWa9oX7vPfS4Vh+YqdLaPw/26u08CzgOuMbMzQwfKoFz4fvwK+AowEdgG3B3tz7rsZlYC/B74vrvvOdyhzezLtuw5c97dvd7dJ5Jc434KMLa5w6I/syp/0+xmdiJwAzAGOAU4CvhhdHi7s+dScUgAw1K2y4CtgbK0ibtvjf7cATxF8ofv48buXPTnjnAJ26SlvFn//XD3j6N/QA3Ab/jyEkZWZTezniT/c/2/7v5v0e6cOO/NZc+V857K3T8FVpC8Ht/fzBrnnUvN+EX+6PkjafulzNikZJ8eXepzdz8IPEgnzn0uFYdVwAnRSIIikjeElgbO1CIz62NmfRsfA38DvEsy8+XRYZcDT4dJ2GYt5V0KXBaNgjgV+KzxMki2aHJN9e9Jnn9IZp8fjT4ZCZwAvNbV+SA5ioTkWupr3f3/pDyV9ee9pey5cN4BzGyQmfWPHvcGziF53+QlYHZ0WNNz3/g9mQ286NHd3q7WQvZ1Kb9QGMl7Jannvn0/N6Hutnfki+Qd9/dJXhe8MXSeVrKOIjky4y3gvca8JK9R/gn4IPrzqNBZUzI/TvIywCGSv2lc1VJekt3Ue6LvxTtAeRZmfzTK9nb0j+OYlONvjLJXAecFzH06ye7928Dq6KsiF877YbJn/XmPskwA3oxyvgvcFO0fRbJorQf+FegV7S+OttdHz4/KwuwvRuf+XeAxvhzR1O6fG02fISIiaXLpspKIiHQRFQcREUmj4iAiImlUHEREJI2Kg4iIpCls/RCR7svMGoeUAgwB6oHqaHu/u/9VkGAiMdNQVpE2MrObgX3u/i+hs4jETZeVRDrIzPZFf04zs/8wsyfM7H0zu9PMLo7m23/HzL4SHTfIzH5vZquir6+H/RuItEzFQSQzvgZ8DzgJuBQY7e5TgPuA70bHLAJ+5u6nABdEz4lkJd1zEMmMVR7NVWNm/wn8Mdr/DnBW9PgcYFzKEgD9zKyvJ9dCEMkqKg4imXEw5XFDynYDX/476wGc5u4HujKYSEfospJI1/kjsLBxo3F9X5FspOIg0nX+CSiPVulaA/xj6EAiLdFQVhERSaOeg4iIpFFxEBGRNCoOIiKSRsVBRETSqDiIiEgaFQcREUmj4iAiImn+P8cTvwJaGRxkAAAAAElFTkSuQmCC\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": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.sub = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZjVdZ3/8ecbUAYYUGEUcUYEV+7D2wHCzNXSQmzjZ4KS5Q26i7s/Ldtyr7XcNW2vTNusta2rdE1M21KxXyuLlJpmmVch400it5JADDjIjdyMw3AzfH5/nJv5zHfOOXPOme/33L4e18U15+Y753w8jLzn/X5/bsw5h4iIiK9PsQcgIiKlR8FBRES6UXAQEZFuFBxERKQbBQcREemmX7EHkI+6ujo3atSoYg9DRKRsvPLKK9udc8dme31ZBodRo0bR1NRU7GGIiJQNM9uYy/UqK4mISDcKDiIi0o2Cg4iIdFOWPQcRkWI4ePAgzc3NtLe3F3soadXU1NDQ0MARRxzRq9dRcBARyVJzczODBw9m1KhRmFmxh9ONc44dO3bQ3NzM6NGje/VakZaVzOxBM3vXzN5M87yZ2XfNbJ2ZvWFmZ0Y5HhGR3mhvb2fYsGElGRgAzIxhw4aFktlE3XN4CJiR4fmLgDHxP/OBH0Q8HhGRXinVwJAQ1vgiLSs5535nZqMyXDILeNjF9g3/o5kdbWYjnHPvZHrd7Ru28aOr7wOg79EbuebeO8MasoiIUPzZSvXAJu9+c/yxbsxsvpk1mVkT8SMoDh1RzxFbR7LxyqvYeOVVvPfY45EPWESk2H71q18xbtw4TjnlFO66665I3qPYDelU+U/K04ecc/cD9wM0Nja66358Pd/9+/voGNzAMwc+zMD9HQz/4TOMWbw4+T1DPvEJjrn8skgGLiJSDB0dHdxwww08++yzNDQ0MGXKFD75yU8yceLEUN+n2MGhGTjRu98AbMn2m/s3rGDnLsdh+jC4tYH3+0PH9gcBGPnnvbQtW8aeeLBQoBCRSvDyyy9zyimncPLJJwMwd+5cnnzyyYoLDouAG83sUWAasLunfoPv+hlnwPInAPjhGxfT1uckfn/c5wA4auQ+6rf9ibOJBQkFChEJ0x3/u4KVW/aE+poTTxjCV/9mUsZrNm/ezIkndv5O3dDQwNKlS0MdB0QcHMzsZ8B5QJ2ZNQNfBY4AcM79EFgCzATWAW3AvJzeoHFe7A9w2h03s3mHo+5ADQC73SnsrDuFBwZv5ujzz2fUzhU08m63QAEKFiJSPmLzd7qKYgZV1LOVPt3D8w64IYz3mn7xpGQWAfDQhmPZ/P7ZwACsz7G8etIYFn9kB6eNHs05Kx0D49cpqxCRfPT0G35UGhoa2LSpcx5Pc3MzJ5xwQujvU+yyUni8LALgmgUXA0/CvKe44+b/YVjrSPo8N5Bt7OeBD23mh//4VQDee+zxZGBQoBCRUjdlyhTeeust1q9fT319PY8++ig//elPQ3+fygkOqbQshwUX87HBE3itYzo17kQGtO2DLTDvV/FAchTM/Ne/Yc7YORkDBShYiEjx9evXj+9973t8/OMfp6Ojg2uvvZZJk8LPYixV/arUNTY2uh4P+2la0FlmalkOx09OZhGD2/rz3qDtABy2Nt6qe4VBpx4EYObJM7sFCogFC4CBU6YAChQi1WjVqlVMmDCh2MPoUapxmtkrzrnGbF+jcjMHv8y04OKUWQTA0XsPc8KeMbRs28xh9vPcqpXMGQvHXH5Zl3/8VX4SkWpSucHBN3l28ub0Ix9n+qmrYN5TACx4eDlb3tiRsuSUyCKga7DwA0X76tXJ50VEKkV1BIdgFuGZd9Xk5O07bv4fhr1fT5/nPsJha2PhG8+y5NQlyecTwcIPFBuvvCr68YuIFFh1BIegeIkpafJsaJzHyFP/KplFJMpNe3e3ANB2aB8vb3iLOWOLNGYRkQKqvuDglZiAWKAAaJzXJYu46esvUtuyn4H7Y/OHj9m3j8Nvu5Qlp/bVq7tkEOpBiEi5q77gEFgP4TergWQWMe2Ck3jy9c3Jy05euZNjW+uxQMnptNHbOWffMGri16kHISKVoPqCQ5CfSXhZxBXTRnLFtJHJp7o0rt8/QB83kGae5/Hxu/jT2eNYMGMBoB6EiETr2muvZfHixRx33HG8+WbKQzZDoeCQZsorkMwioGvj+us3v8AxrSdQ8+urabBNbDvh7czn3YmIhOSaa67hxhtv5Kqrov1FVMHBlyaLCDrh1GFseWMHAMe8X8fhLfuTvYi5O1dT33IwmUGo/yAiYTr33HPZsGFD5O+j4ODLMOXVF5z+2of+yfu/Gd/B+RzBeNR/EKlov7yl85fIsBw/GS6K5mS3XCk4ZJJmyquvH0fTzx2d7DnMYx53nLWGcUP7Mvc/Ydi+bZxUyDGLiIRAwSGdDFNeg2r3HeYX97wKwPR9c6mre4XtQ9+i7dA+2Lcj6pGKSDGUyG/4UVFwSCfLKa/vHtO3y7fZjgGcPWAGl8z4Ck//51Tqmlt5+uKpAPT9+Plc8Pm7C/VfICKStz7FHkDZmDw7Vg+EWJCI7/jaUtePN8b055IvncklXzqTuoba5Lf0/fj5bI/fr2tupePp3xR82CJSWT796U8zffp01qxZQ0NDAz/60Y8ieR9lDtlKM+X1th27eeLAdC6/L/bUqe/sp672SIBYlvD52OOJ7EFEpDd+9rOfFeR9FBzy4fUjxhzewOwj4Wt8CoC2/Ydo93oQY6cOZ9KH6wG6lJhAZSYRKV0qK+WjcV5sy+95T3Fk/WlMGnEUj10/nceun07r8f1pHRD7WLc3t7L25a1A1xITqMwkIqVNmUPIWur60VLXj1uvPzOZPUDXEhPEykxth/al3MhPRKTYlDkUybABwxjYbwAAa3auYcnbS3r4DhGRwlHmEAZviuttO3bz0oDzgelArLSUqv9w7IBjGbxxB7f/dwerdx7iN+PfYB7KIkSkNCg49FZgsdyog28nb4+dOjx5e3tzK0AyOAz5xCeSz9W3HOR8juDRC2NZBKDgICJFpeDQW4HFchvuPIe2Ax1cft8fAJh1Tj1XTBvZpf8AdDtqdDywYMYC5v1qHmt2rkn2IkCZhIh02rRpE1dddRUtLS306dOH+fPnc9NNN4X+PgoOIaur7c/21v0ArHxnD0CXcyF6MvPkmV3uK5MQEV+/fv245557OPPMM9m7dy9nnXUWF154IRMnTgz3fUJ9NWH44BqGD67hsXnTk9lDgt9/gK49iMRRo1OBC7xtvoOZhLIIkeo2YsQIRowYAcDgwYOZMGECmzdvVnAoC97q6URz2u8/QNcehN9/CG7z7WcSyiJESsfdL9/N6p2rQ33N8UPH889T/znr6zds2MBrr73GtGnTQh0HKDiEz2tQ+83pSR+uT2YJQJcMIth/8M0ZOycZDNSPEJGE1tZWLr30Uv7jP/6DIUOGhP76Cg5h8xrUG+48J6+XSJSYEhKnyakfIVI6cvkNP2wHDx7k0ksv5TOf+Qyf+tSnInkPBYeIdZm5dHp9l+Z0qjUQfokJupaZ/CwC6JJBiEh1cM5x3XXXMWHCBL74xS9G9j6RBgczmwHcC/QFHnDO3RV4fiTwY+Do+DW3OOcqZqlwXW1/at9bxW07/om2Ax281nYBTPs3IP0aCL/EBN3LTEFqVotUl5deeolHHnmEyZMnc/rppwNw5513MnPmzB6+MzeRBQcz6wt8H7gQaAaWmdki59xK77J/AR53zv3AzCYCS4BRUY2p0Iaf/VlY/gSTgPf/8hoD93VutOf3IIJrIIL8MtMQbyaTmtUi1eecc87BORf5+0SZOUwF1jnn3gYws0eBWYAfHByQ6KQcBWyJcDyFl0P/Id02G5lmMmVqViuLEJHeiDI41AObvPvNQHC+1e3AM2b2OWAQcEG6FzOz+cB8gJEjs19UVg4ybbORaSaTT1mEiIQpyuBgKR4L5kKfBh5yzt1jZtOBR8zsA865w92+0bn7gfsBGhsbo8+pIjDq4NudZ1BD8hzqXEpM6WjKq4iEKcrg0Ayc6N1voHvZ6DpgBoBz7g9mVgPUAe9GOK6ieGnA+bQd6GDgO7uBWKBo3dvO8MbcZhylm+bq05RXEemtKIPDMmCMmY0GNgNzgSsC1/wF+CjwkJlNAGqAbRGOqWhqz/47vvV65z/aN7/zRQa27md4imvTbbORaZqrL9WUV/UjRCQXkQUH59whM7sReJrYNNUHnXMrzOxrQJNzbhHwJeC/zOwfiZWcrnGFaMMXwRXTRnZZ47Dizr4pr8u0zUau01wT1I8QkVxFus4hvmZhSeCx27zbK4EPRTmGUtalB5Gi/wD5T3P1aVaTSOVob2/n3HPPZf/+/Rw6dIjZs2dzxx13hP4+WiFdJLEN+WASxDbqgy7nQvjymeaajrIIkfLWv39/nn/+eWprazl48CDnnHMOF110ER/84AdDfR8FhyJ5buBMnhs4k8fmTe86gykgjGmuvmAWISLlxcyora0FYnssHTx4ELNUk0N7R8GhiFa+s4fL7/sDt+3YTV1t/5TN6TBWUmeiKa8i+Wm58072rwp3y+7+E8Zz/Fe+0uN1HR0dnHXWWaxbt44bbrhBW3ZXklmnd/YV2g50sD3NzKVs9bbEBCoziZSLvn378vrrr7Nr1y4uueQS3nzzTT7wgQ+E+h5WjpODGhsbXVNTU7GHEZoVd57DqINvM2jkGbEH4s3poF/c8yrbm1upa6hNPub3IBI2XnkV7atXUzN+fPKxbDKJRLN63NBxgLIIkaBVq1YxYcKEYg+jizvuuINBgwZx8803Jx9LNU4ze8U515jt6/YJb4iSr5cGnM+GI06O3WlZDsufSHnd2KnDuwSG7c2trH15a7frhnziE10CQ/vq1exZvLjHccw8eWYyMKzZuYYlb1fMBrkiFWPbtm3s2rULgH379vHrX/+a8d7/72FRWakEZNucznaaa77rITTlVaT0vfPOO1x99dV0dHRw+PBhLrvsMj4RWCAbBgWHMpdummtQrs1qTXkVKU2nnnoqr732WuTvo+BQIvyZS2MOb+DIFBv0BWWa5urLp1mtjfxEqpuCQwnwZy49cWA6s4+ML46DjAvksp3mms96CJ9mNYlUHwWHEuDvu3T5ffA1PhXrP0DGHkS+ci0xaSM/keqj4FAOWpZ324MpX/mUmILUjxCpfAoOpW7y7M7bee7B5AuWmLI5HyJIs5pEKp+CQwlKNKdjxjPr9B/Eyk557sGUTrbnQ2SiLEKkMik4lBi/OQ2xQAF0OQsilXyOGs13PYRPs5pECq+jo4PGxkbq6+tZnMUC13woOJSY4KFAnRlEnN9/gF73IILy2bwvQbOaRArj3nvvZcKECezZsyey91BwKCd+/wEy9iDSHTWaSW+b1ZrVJBK95uZmnnrqKW699Va+/e1vR/Y+Cg7lpHFe10CQpgeR6ajRTDI1q3PNIqBrJtG0tYmmrU3J/ZoUKKTcvfj4WrZvag31NetOrOXDl43NeM0XvvAFvvnNb7J3795Q3ztIwaEC5XrUaCphTHn1M4mFaxcmA4PKTSL5Wbx4MccddxxnnXUWL7zwQqTvpeBQBvzZS7NOr+/anA5xDYQvjCmvPjWupdL09Bt+FF566SUWLVrEkiVLaG9vZ8+ePXz2s5/lJz/5SejvpeBQ4vzZS91mLvk9iI2/j/1JbPcdCBTZbtCXShhTXn1qXIvk5xvf+Abf+MY3AHjhhRf41re+FUlgAAWHktd1a43AzCW/B9G0oDMwBBrV+ayB8IUx5dWnxrVI6VNwqBR+oFhwcZdy0yRg0l/HMol8+g+p9LZZ7dNCOpHcnXfeeZx33nmRvb6CQyXKYcprPsJoVvu0HYdI6VFwKDNdt9ZI0aCG1FNeE5lEy2Vs31+fd/8Bwp/y6lMWIVIaFBzKSL5ba/iZxFj7X3DnQstRbG87Fva25BwcfIXMIkCZhBSfcw4zK/Yw0nLOhfI6FtYLFVJjY6Nramoq9jCKLpFBPHb99Oy/yWtc/+L1i+HIQVzy71eHMp5EFlHjHXbem0zCXxsBsYV0AI3DGwEFCim89evXM3jwYIYNG1aSAcI5x44dO9i7dy+jR4/u8pyZveKca8z2tZQ5VBu/5PRPP4YD74e2V1PYU16Ds5q0kE6KraGhgebmZrZt21bsoaRVU1NDQ0NDr19HmUMZyytz8Pzi9iVs3wZ1A+M/6O27GTvgd0yasC92v5eL6oKZRG/7Eb5EyWnc0HGAsgiRnihzqDIZV0/3YOxHT4OXtwLHA7B9405wRzGJx0OZ4RR2P8KnxrVItCLNHMxsBnAv0Bd4wDl3V4prLgNuBxzwJ+fcFT29rjKHmJ8u/QtPvr4ZiAWJiSOG5J1FQOceTJd86czOGU7HT449WUZZBCiTEAkqmczBzPoC3wcuBJqBZWa2yDm30rtmDPBl4EPOuffM7LioxlOJMq6e7q1MW3Mkns8hWBQqiwDtACsShijLSlOBdc65twHM7FFgFrDSu+bvgO87594DcM69G+F4JAudezCdxtipH4tNc/W35oC8Sk5hb+TnU+NaJHxRBod6YJN3vxmYFrhmLICZvUSs9HS7c+5XqV7MzOYD8wFGjsy+rl5NetN/gAx7MGVaVAe9yiIg/ExCayVEei/K4JBqEnCwwdEPGAOcBzQAL5rZB5xzu7p9o3P3A/dDrOcQ7lDLX8bdW7OU9TnUOewGm0qqjfzCXGXtU8lJJD9RBodm4ETvfgOwJcU1f3TOHQTWm9kaYsFiWYTjqkiR9h+C0u0Gm2dvIsp+hEpOIvmJbLaSmfUD1gIfBTYT+wf/CufcCu+aGcCnnXNXm1kd8BpwunNuR6bX1mylzHq7/gFimcP25lbqGmqBLPdgCvYmNv4+9vWkc2JfswgUUc5qCtJaCakmJTNbyTl3yMxuBJ4m1k940Dm3wsy+BjQ55xbFn/uYma0EOoB/6ikwSGHkdQZEsDeRKatIEyiizCKCMp1xnXhewUKqlVZIV6DL7/tDct1DQj4N6oQu6x/yFQwU0JlRQMpgEfZeTZloHyepdJFkDmbWCHwYOAHYB7wJ/No5tzOvUUqk8t69NQP/mFHIY6vvdH0KSJtVRD2ryafehEhXGTMHM7sG+DywHngFeBeoITYF9UPEgsS/Ouf+EvlIPcocctPbHsSKFzez9uWtyfuJXkSvMglfpqzCyygK2Y/wqTchlSDszGEQsdXL+9K82enEZhcVNDhIYflTXKGHaa75yHL205Che6D+aCD6foRPvQmpRnn3HMzsSOfcgZDHkxVlDrkJ9iB603+AkHoQ2cgw+2njT7fQvsNRM+lUoHBZhHoTUq6i6jm8AFzjnNsQvz8V+C/gtDzGKAUWxgK5IL8Hkc9Ro1nJMPtpyPAWaB8ALctp29RO27Jl7Fm8OHlpVMFCvQmpFlllDmb2cWK7q36X2LYYFwF/65wLub6QHWUO+QtjDYTfgwi9/5AtL1C899yr7Nk4AGqOAqD93QPUjK7npP95pqBDUm9CSlkkmYNz7mkz+3vgWWA7cIZzriXPMUqZy3qbjSh5WcUxkxdwjFd+2vjgKtr/vJGNF54BwJCJtRxzeue03t5uP56OehNSSbItK/0rcBlwLnAq8IKZfck591SUg5No9HaDvqCClJgyCZSfhmy7CZ79HRDLIqC1MzjksRdUtvySU7A3oZKTlJtsy0r3ArckZi2Z2UnEDu+5MOLxpaSyUv7CPiCoJEpMGXRbSLe3hSEn7YsFizy298iXSk5SbFGVlW4K3N9I7BAfKTNhb9BXEiWmDIIL6dpWb6JtNezZORr2TosFipNIvWmgr5eBI1PJSYFCSpHOkJaKFtwe/L3HHk/OamrfvAsGj+eYeQ93nzbrC6EUla7kpHKTlCrtrVTFotiDyd/JFYrUg8hS1iuue9oXypdj4ND511IoJbMrq5S+sPdg8ndyhRx2cy0Sv+TUtmxZ+rUSmfaF8uWRYegwIilVPe2tdG6Wr7OhkPsrKXOIRhhrIHwFW0kdAr/cBCQzipMeeTj7F/EDR8tyOH4yzMttQl+w5DRu6DgWzFiQ02uIpJJr5tBTcMj2p/IX8fMZCkLBIRpRbLOR84FBJaLXm/wlztk+fnLnY70sOSmLkN4ItazknItmXp+UpLC32cjrwKASkanklFWg8M/Zhl6XnLSoTgpNDWlJqZpLTEFdZjjlU26CrLclTyfVojqVnCQXoZaVSpWCQ/SiCA7lWmLyhXI6XQi9CZWcJFearSShCXObjXIuMfm6LarLp+Tkz35K9CYWXNz5fC9LTgoUEoZs91Z6xDl3ZU+PSeUIu/9Q6iups5VxUV0+BxDl2ZvQojqJWrZ7K73qnDvTu98XWO6cmxjl4NJRWamwoi4xQfmWmXyhHGPay5KTFtVJOqGWlczsy8BXgAFmtifxMHAAuD/vUUpVK7fFctnyS055H2OaqeSUx6I6ZRKSr2wzh284575cgPFkRZlDYYW9zUZQOc9kSmfjlVcB5D6rydfTth1ZBAs1riUh7MxhlHNuQ7rAYGYG1DvnmnMcp5SRsLfZSKXoZ0JEoH316mSQgDzKTJm27ciyN+FnEsoiJBc9NaT/3cz6AE8CrwDbgBrgFOB84KPAVwEFhwrmb/MN4Wz17auUmUy+4KymvMtMCRnO06Zleec1AX7jet6vtKZVstdjWcnMJgKfAT4EjAD2AauAp4AnnHPtUQ8ySGWl4gp7mw1fpayHCAqlzJROcKuONFmESkzVLfR1Ds65lcCtvRqVVJSwp7n6KjGLSPDLTHnNZErHnw6b6tCieLBQiUly0dPGe1OATc65lvj9q4BLgY3A7c65nQUZZYAyh9IR9jRXXyU1qkPZgiMbqXoT0G2rDk15rT5h78r6KnCBc25nfPvuR4HPAacDE5xzs9N+c4QUHEpHIUtMUBllpkhLTEFpZjwtpJUl9n7ysjUcYNygehbM+WX0Y5KiCLus1NfLDi4H7nfO/Rz4uZm9nu8gpXIUqsQElVVmiqzEFJRmxtMcapnjOoPuvP1roX1dTmsqpLL1GBzMrJ9z7hCxmUnzc/heqQL+TKawZzH5W25AeW+74QtlsVw+gjOefAsvYs37m5nHVmjfzczf3sqcFH0LqR49/QP/M+C3Zrad2CylFwHM7BRgd08vbmYzgHuBvsADzrm70lw3G1gITHHOqV5UxvzN+iDcMhNUxnoIf38mfx1EMc087VpI7M+0bTnUHMWcRMU5VZM7QUGjYvV02M/Xzew5YlNYn3GdDYo+xHoPacX3X/o+cCGxdRDLzGxRfPaTf91g4PPA0vz+E6RURL1YrlJnMvV6sVwIgush1uxcw7yhx8WerD2bma1tdGtVZ5gZJeUvm6msf0zx2NosXnsqsM459zaAmT0KzAJWBq77N+CbwM1ZvKaUsKgXywV3dq2ELCL0xXIh6LY/08HdMGIcc4IHC/W0atunoFF2IjvsJ14qmuGc+9v4/SuBac65G71rzgD+xTl3qZm9ANycrqxkZvOJ9zxGjhx51saNGyMZt4QnyplMK17czNqXtwIkZzRVwpTXUHZ2DVnWi+eCwSIhzwONJFyldNiPpXgsGYni23J8B7gmmxdzzt1PfCfYxsbG8ju+rgpFOZOpUs6HCCpaszqDrBfPpWt4+wcZSdmIMjg0Ayd69xuALd79wcAHgBdi+/dxPLDIzD6ppnRliHImU5BfYoLyLTOVYrM6ZT8ivk+TFs5VriiDwzJgjJmNBjYDc4ErEk8653YDdYn7PZWVpPyFeeyoT+shCkdbcFSPyIKDc+6Qmd0IPE1sKuuDzrkVZvY1oMk5tyiq95bSU6gSE1ROmakUS0yZsgjIkEnkcU62FFekC9mcc0uAJYHHbktz7XlRjkWKK1hiino9RCUIlphKYcqrL+tT54LnZGfYYlxKh1Y5S8Hp8KDcleKUVz+LgAz9iGCjWg3qsqDgIAWnw4Ny52cRUDrNap/6EZVFwUFKQpjN6kqd5hpUas1qzWqqLAoOUnRRNqsrVSk2q309ZhF+g1rN6ZIU2QrpKOk8h8oV9qrqSjowKJ1SXFXt63aw0N6W+F5NtVo9XUCltEJaJGdRZBGVskAunXLKIiCwV5Oa0yVLmYOUrDCOIPX3YILK2ocplWAWASWeSbQsZ6YbxJx5vy32sCqeMgepKL1dD1GpC+TSKcUpr0F+JtFk+2my/Sx5KPZv1swTPsycj32nWEMTj4KDlKxCrIeoNKmmvJbyrKaFz/wjS7a8CMAa1w5bXux+boQUhcpKUjbCaFYnzoGoa4idn1xp/Yeg9x57nD2LFwMky00nPfJwkUeV2ryHGlnDAcYdf1byMU2BDY/KSlKxwmhWV+ICuUwybcFRClmEb6Yb1GWj/6atTTRtbWJJ/PhSBYrCUnCQspFpf6Zss4hqWSCXSqnPappDLXPeWQ/uXQAWMpQltQMBBYpiUHCQsqSFc7kr9Y38ghv0zWlZz5zjJ8OcBSxcuzAZGLQ1R2Go5yBlL9iLgOwyiWrrP/j8XgSUaD8isQYisECu26I6lElkQz0HqTr5zmqqtv6DrxxmNaUTXFSnklM0lDlIxcln8Vw1bLORiZ9JtC1bBsDAKVOSzxclWCy4uHN7DUi7B5NfcmraGvt3oXF4fN2EAkWSMgcRojuStFL5mUSqklPimoLyexAZDgjqsm5CvYnQKHOQivPTpX/hydc3AyR7ET1lEcH+A1RXDyKTktjYL03/IRP1JrpS5iBVLzjlNRt+/wGqrweRiT8Ftm3ZMtqWLUtmFgUNFDlu8531MaaSkjIHqWi9mckE1duDSKdoK66bFsDyJ2K389zmO5hJVFsWocxBxNOb/Zkq7RzqMBRtrYR/DnWe23zrGNPcKDhIRcv3vOpqnuaareAOsEUtOWUheIypZKbgIFUnm5lM1bzNRraCayWCJafENZHw+w+Q11GjOuM6MwUHqSr5bruhElPPCrbJX2CbjUzTXNNRialnakhL1cp2C3D/NLlKP0kuLAVtXOcxzdVXLVNe1ZAWyVK2WYRKTLkLZhGRy3Gaq09TXlNTcJCqlc96COhaYgKVmbIR6aymLFdSp+M3qqEzk6j2foSCg0hcNudVa7Fc7iI/1zqEafrczZAAAA5ASURBVK4+9SNiFBxEyH49hF9iApWZspFqB9hI9aLEBN2nvPpZBFRPJqHgIELq9RDZbt6nmUwlpJclpqBq7kdEGhzMbAZwL9AXeMA5d1fg+S8CfwscArYB1zrnNkY5JpFsZNus1mK5/EQ2zTVYYurleohU/YhqEVlwMLO+wPeBC4FmYJmZLXLOrfQuew1odM61mdk/AN8ELo9qTCLZyva8as1kyl3BzrIOYT1EKtXSrI4yc5gKrHPOvQ1gZo8Cs4BkcHDO/ca7/o/AZyMcj0heclk4pxJTzwo2zdXPIiD0ZnWln0AXZXCoBzZ595uBaRmuvw74ZbonzWw+MB9g5Egd3CKFkymLgM5MQiWm/BT0eNIQm9WVfrBQlMHBUjyWcjm2mX0WaAT+Ot2LOefuB+6H2ArpMAYokqtMs5pUYspdwUpMEHqzutJnNUUZHJqBE737DcCW4EVmdgFwK/DXzrn9EY5HpNd6M6tJuivoSupMzeo8sghfJc5qijI4LAPGmNloYDMwF7jCv8DMzgDuA2Y4596NcCwikcjUj9BK6twV7HyICLMIqIxV1pEFB+fcITO7EXia2FTWB51zK8zsa0CTc24R8O9ALbDQzAD+4pz7ZFRjEglbun7E8e2HmHhU5/9e6kH0LPKV1L6Qp7wGVcIq60jXOTjnlgBLAo/d5t2+IMr3FykkP4t4/mAbLSOG8Nj1sd1b1YPoWcFXUidEMOU1Uz+iXLIIrZAWCUmmWU2nvrOfutojizm8slSQmUypprxG1I8opyxCwUEkAsFZTXv2HaSmtYOv3/wCACecOox5V00uwsjKR0FnMvkintVULnTYj0gBLHh4OVve2AHAgPc72DeoL7d+67ziDqqMJE6Wqxk/Hijg+dSJLOJ4L5D3IpNIBIcFMxaEMbqc6LAfkRLkZwlfv/kF2vYf0vTXHJREFgGhZBLlsh5CwUGkwOpqj2R7/PbS9TtZun4nT76+Ofm8gkV3BT9ZLiHkfkQ5rYdQcBApsOMG19Bn9yHmtvbnIzXHsPKIDlriz/W0d5PEFGw9RJCfSWz8fezP8ic6n+shUJTTLq8KDiIF5u/B1Gf3Ic5rqOWS+JRXrbjuWUHXQwT5mUTTgs7AENKOr6VEDWmRIkqsf7jkS7Hg8NOlf0mWmJau3wnAtNFDAQWKdIrWrPbl2bhOrIEYN3QcEG3/QQ1pkTLjb7MxALhtamwTv2CgUG8itaI1q315Nq5LeQ2EMgeRIlrx4mbWvrw1eX97cyt1DbXJTCLBDxQQ601MHDGEx66fXrCxloOSyCKgeyaRZRYB0U1zzTVzUHAQKSG/uOfVZICA9Jv1JXoTE0cMAZRFJLz32OPsWbwYgLZlywAYOGVK8vmCBQu/H7Hx97GvJ50T+5omUERdYlJZSaSMZXtgkL8CO1hyquZA4U959QMFlEjjOjjDCZLBotRKTMocREpUsFmdTqYmNlR3sPCVRMnJDxSQNquIosSkzEGkgmRzJrW/4V+q3kTimmpXEo3r4KK6dFmFbWVNX4q6klqZg0iJ8pvV6RrVPVFvIrWSyCKCvECxcOdrLBk0CGqOAmANBxg3qJ4Fc36Z98urIS1SgYKNasjuZDk/k9AMp04l07hOJ1B+mrd/LQAL+o+NPZDH5n8KDiIVKNspr5koi0gtVeO6Zvx4Tnrk4SKOqqt5Cy9izfubGceR0L6bme+/z5yhZ6S+OE3gUHAQqQLZTnn1KYvITimWnBauXciSt2OHaq7ZtpxxHbDADe9+YYZpswoOIlWgt/2IYBYByiQSSr3klHE9RIb1FXbtEgUHkWrS2ywClEmkEyw5BYNFMQJFlywiHiRSTnkNBAq7Y4+Cg0g18bOILW/tAuCEMUcnn88mWKgfkZ1MWUUxAkUwi4A0U15/eQs2824FB5FqlW/jWgvpclcK5Sc/iwBo2hr7d7FxeCwG+IFCPQcRSVLJqTBKpfyUqeSk4CAiSVpIVxylUH4KlpweuughBQcR6S6fLAJ0AFFv9VR+8oUZOPwsomlrE29e86aCg4h0F0bjWr2J3gmWn3xRZhh3v3w3t0y7RcFBRDILY8V1sDeRKlgkKGj0LOoMQz0HEclZviUnXzBYJGQKGqDAkUouGUZQusCh4CAiOeup5OTLNXCkCxrQc+DwKYjE5Bs4Rv3kEQUHEclfsOTkCwaOfDIMX6bA4VMQyU6mwKHgICKRiTLDyCSKIJJJJQaYkiormdkM4F6gL/CAc+6uwPP9gYeBs4AdwOXOuQ09va6Cg0jx5ZJhZFKMIJJJWAEmnWIFnpI5JtTM+gLfBy4EmoFlZrbIObfSu+w64D3n3ClmNhe4G7g8qjGJSHgmfbg+7T/qmQKHb8tbu9jy1q6srs3GAGAu/ZP3w2ysh2Hp+p0sXb8zstcPU5RnSE8F1jnn3gYws0eBWYAfHGYBt8dvPwF8z8zMlWOtS0SSMgUOX7ZBJB/5Bp5ggAnTR2qOYXvrAXhrfySvn8njOV4fZXCoBzZ595uBaemucc4dMrPdwDBge/DFzGw+MB9g5MjKqgWKVKtsg0g+ogw8+TpucA3HDa4p9jCyEmVwsBSPBTOCbK6JPejc/cD9EOs59G5oIlLpogw8Zenm3C7vE80ogFimcKJ3vwHYku4aM+sHHAXsjHBMIiKShSiDwzJgjJmNNrMjgbnAosA1i4Cr47dnA8+r3yAiUnyRlZXiPYQbgaeJTWV90Dm3wsy+BjQ55xYBPwIeMbN1xDKGuVGNR0REshdlzwHn3BJgSeCx27zb7cCc4PeJiEhxRVlWEhGRMqXgICIi3Sg4iIhINwoOIiLSTVnuympme4E1xR5HnupIsQK8TGjsxaGxF0eljf0k59yx2b5ApLOVIrQml90FS4mZNWnshaexF4fGXhxhjF1lJRER6UbBQUREuinX4HB/sQfQCxp7cWjsxaGxF0evx16WDWkREYlWuWYOIiISIQUHERHppqyCg5nNMLM1ZrbOzG4p9nh6YmYbzGy5mb1uZk3xx4aa2bNm9lb86zHFHmeCmT1oZu+a2ZveYynHazHfjf9dvGFmZxZv5GnHfruZbY5//q+b2UzvuS/Hx77GzD5enFGDmZ1oZr8xs1VmtsLMboo/XvKfe4axl/znHh9LjZm9bGZ/io//jvjjo81safyzfyx+5ABm1j9+f138+VElOPaHzGy999mfHn88958b51xZ/CG27fefgZOBI4E/AROLPa4exrwBqAs89k3glvjtW4C7iz1Ob2znAmcCb/Y0XmAm8Etip/l9EFhagmO/Hbg5xbUT4z8//YHR8Z+rvkUa9wjgzPjtwcDa+PhK/nPPMPaS/9zj4zGgNn77CGBp/DN9HJgbf/yHwD/Eb/9f4Ifx23OBx0pw7A8Bs1Ncn/PPTTllDlOBdc65t51zB4BHgVlFHlM+ZgE/jt/+MfB/ijiWLpxzv6P7SXzpxjsLeNjF/BE42sxGFGak3aUZezqzgEedc/udc+uBdcR+vgrOOfeOc+7V+O29wCpiZ6uX/OeeYezplMznDhD/DFvjd4+I/3HAR4An4o8HP/vE38kTwEfNLNVRx5HLMPZ0cv65KafgUA9s8u43k/kHsRQ44Bkze8XM5scfG+6cewdi/3MBxxVtdNlJN95y+fu4MZ5GP+iV8Epy7PEyxRnEfgssq889MHYok8/dzPqa2evAu8CzxLKZXc65Q/FL/DEmxx9/fjcwrLAj7hQcu3Mu8dl/Pf7Zf8fM+scfy/mzL6fgkCpCl/o83A85584ELgJuMLNziz2gEJXD38cPgL8CTgfeAe6JP15yYzezWuDnwBecc3syXZrisVIbe9l87s65Dufc6cTOuJ8KTEh1WfxrSY0/OHYz+wDwZWA8MAUYCvxz/PKcx15OwaEZONG73wBsKdJYsuKc2xL/+i7wC2I/fFsT6Vz867vFG2FW0o235P8+nHNb4/8DHQb+i84SRkmN3cyOIPaP63875/5f/OGy+NxTjb1cPnefc24X8AKxevzRZpbYd84fY3L88eePIvtSZmS8sc+Il/qcc24/sIBefPblFByWAWPiMwmOJNYQWlTkMaVlZoPMbHDiNvAx4E1iY746ftnVwJPFGWHW0o13EXBVfBbEB4HdiTJIqQjUVC8h9vlDbOxz47NPRgNjgJcLPT6IzSIhdpb6Kufct72nSv5zTzf2cvjcAczsWDM7On57AHABsb7Jb4DZ8cuCn33i72Q28LyLd3sLLc3YV3u/UBixXon/2ef2c1Osbns+f4h13NcSqwveWuzx9DDWk4nNzPgTsCIxXmI1yueAt+JfhxZ7rN6Yf0asDHCQ2G8a16UbL7E09fvxv4vlQGMJjv2R+NjeiP/PMcK7/tb42NcAFxVx3OcQS+/fAF6P/5lZDp97hrGX/OceH8upwGvxcb4J3BZ//GRiQWsdsBDoH3+8Jn5/Xfz5k0tw7M/HP/s3gZ/QOaMp558bbZ8hIiLdlFNZSURECkTBQUREulFwEBGRbhQcRESkGwUHERHppl/Pl4hULzNLTCkFOB7oALbF77c5584uysBEIqaprCJZMrPbgVbn3LeKPRaRqKmsJJInM2uNfz3PzH5rZo+b2Vozu8vMPhPfb3+5mf1V/LpjzeznZrYs/udDxf0vEElPwUEkHKcBNwGTgSuBsc65qcADwOfi19wLfMc5NwW4NP6cSElSz0EkHMtcfK8aM/sz8Ez88eXA+fHbFwATvSMAhpjZYBc7C0GkpCg4iIRjv3f7sHf/MJ3/n/UBpjvn9hVyYCL5UFlJpHCeAW5M3Emc7ytSihQcRArn80Bj/JSulcDfF3tAIuloKquIiHSjzEFERLpRcBARkW4UHEREpBsFBxER6UbBQUREulFwEBGRbhQcRESkm/8PKaasmNwOb2wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "surv = model.predict_surv_df(x_test)\n",
    "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": 23,
   "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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6573402857263979"
      ]
     },
     "execution_count": 24,
     "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": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU1dn48e+dTDIhCwkQEiAJSyDsyBZZFBEsKFqrVm3Frdbdtla7vnV527r091q7vbVVW3HXVtwXXuuGsrggqyxhC4QA2SCQBBKyZzLn98c8YAiTZJLM5JlM7s915crMs2Tuhwlz55z7OeeIMQallFKquTC7A1BKKRWcNEEopZTyShOEUkoprzRBKKWU8koThFJKKa8cdgfgL4mJiWbo0KF2h6GUUt3Khg0bSowx/b3tC5kEMXToUNavX293GEop1a2IyP6W9mkXk1JKKa80QSillPJKE4RSSimvNEEopZTyShOEUkoprzRBKKWU8koThFJKKa9CZhyEUh1VXtPAtqJydh44xuTBCUwe3MfukJQKCpogVI/jdhs25B3hvawDLNt5iP2l1Sf2hQncfk4Gd5wzAke4NrBVz6YJQvUY+WXV/GvNft7eWEhxRR2RjjDOGpHIFaenMX5QPMMSY3jkk9387ZPdfL77MI8snExa32i7w1bKNhIqK8plZmYanWpDebM6t5SnPtvLJzuLCRPhnNFJXHjaQM4ZnURcVMQpxy/ZXMS9b2ZR3+jmmxMGcuX0wWQO6YOI2BC9UoElIhuMMZne9mkLQoUsV6ObP36UzRMrc+kXE8mP5ozgqumDGZTQq9XzLpo4iCmDE3hiZS5vbyzkzY2FjEiK5Xszh3DZlFRinPrfRvUM2oJQIamsqp47Fm/k85wSrp4+mF9fOJaoiPB2/5zqehfvbjnAv9fksTn/KHFRDq7ITOOaGUMYmhgTgMiV6lqttSA0QaiQsLekinX7yiitrKe0so73tx7kcGUdv7t4PN89Pc0vr/FV3hGe/WIf72cdwOU2TEpL4JJJg5gzKomSyjr2llRReLSGszISmTqkr19eU6lA0wShQprbbZj9x+UUHKkBICoijGGJsTx82QROS03w++sVV9Ty9sZC3t5UxI4DFV6PuWr6YH61YDTxvU6tcSgVTLQGoULaZzklFByp4feXTuCiSYOIjgzsr3Vy7yhuPXs4t549nF3Fx1i3r4xB8b0YmhhD35hIHl22m6c/38vS7cX813mjOHfsAOKj204UbrehpLKO/CPVFBypIb/M8722oZE+MZH0iY4kKc7JmSMS9e4q1SW0BaG6vdte3MDafWV8efc5OB3trzMEQlZBOXe9uYVtRRWECUwZ3IezMvozIN5JXFQEcVEOyqrq2VV8jOyDleQerqTgaA31LvdJPycx1kmvyDCOVjVwrM51YvvI5FjOGZ1M/zgnlbUuKusaiIuKYOG0NJLiorr6clU3pl1MKmQdPlbHzIc+4YZZw7jngjF2h3OSRrdhU/5RVmYfYnn2YbIKy085xhEmDEuMYURSLIP7RpPapxepfaJJ69uLlIRoekV+nfDqXW4KjlSzPPswy3YWsya3DJfb8/+3V0Q4ta5GIsLD+M7UVG6Znc6QflpEV23TBKFC1j9W7OHhD3byyc/PZnj/WLvDaVVlnYuj1fUcq3VxrNZFfK8IhiXGEOno2Ijt6noXDS5DjDMcR3gYe0uqWPRpLm9sKKDRGF64YRpnjkj081WoUKM1CBWS3G7Dy+vymD6sb9AnB4BYp4NYP46hiI50QOTXz4clxvDQpRP46bwMLn7sCx5bnqMJQnWKTjajuq3VuaXsL63mymmD7Q4lqCT1juLamUNYtaeUnQe932WllC80Qaigsr+0irvfzGLp9mJcjV8XbLcWlnPH4o1c+PfP+PNH2WwrKueltXnE94pgwfgBNkYcnK48fTBREWE898U+u0NR3Zh2MamgcfhYHdc8vYb8shoWr80jubeTSyansL2ogs92lxDrdDBqQByPLc/h78tyALj+zKEdGiEd6vrERPLtyam8+VUB/7VgNH1jIts+SalmNEGooFBZ5+L659ZScqyeN35wBqWVdby8Lp8nP80lMdbJXeeP5qrpg+kdFUFJZR1Ltxezbm8ZN5+VbnfoQev6M4eyeG0ei9fm8aO5I+wOR3VDeheTsl29y82Nz69j1Z5Snrouk7mjkk7sK69uoFdkeIfv9OnprnlqDTmHKvnsV3OJ0PUtlBet3cWkvzHKVnml1dz4/Do+213C7y+dcFJyAIiPjtDk0Ak3zBrKwYpa3tlURKj8Mai6jnYxqYArr27gjx/t5IucUmYO78e5Y5OZPLgPT3++l3+u3IMjTPjdJeP5TqZ/JtVTX5szMonh/WP4xWubefiDnUxKS2BSWgLjU+IZN6g3ibFOu0NUQSygXUwisgB4BAgHnjLG/L7Z/p8BNwEu4DBwgzFmv7XvOuC/rUN/Z4x5vrXX0i6m4GOM4a2NhfzPezsoq6pnRno/NuUfpbq+8cQxF00cxD0XjGFAvE4PESjFFbV8uO0gG/OOsin/KHtLqk7sG9A7irvOH80lk1NsjFDZyZaR1CISDuwC5gMFwDrgSmPM9ibHzAXWGGOqReQHwBxjzBUi0hdYD2QCBtgATDXGHGnp9TRBBJeNeUd46L2drN1XxqS0BH53yXjGp8RT29DIl3tKWZ1bypxRScwc3s/uUHuc8uoGth0oZ3tRBe9uOUBWYTlPXDOVeWOT7Q5N2cCuBDETuM8Yc571/G4AY8xDLRw/GXjUGHOmiFyJJ1ncau17AlhhjFnc0utpgggOe0uq+OOHO3kv6yD9YiL5+bmjWHh6GmFhulxnMKqsc3HVk6vJPniMl26erutY9EB2TbWRAuQ3eV4ATG/l+BuB91s595Q2sIjcAtwCMHiwjqa1Q3FFLS+vzWfnwQp2HKhgf1k1vSLCueMbGdwyO92vU0so/4t1Onj2+6dz+T+/5Ppn1/HcDdMY0jcaR3gYURFhQTM7rrJHIP/3evuT0WtzRUSuwdOddHZ7zjXGLAIWgacF0bEwVWf8+KWNrNtfxtB+MYwZ2JvLp6by3dN1yunupF+skxdumMZl/1jFpY+vOrE9OjKcp67L5IzhOp9TTxXIBFEANL0tJRUoan6QiMwD7gXONsbUNTl3TrNzVwQkStVhX+UdYe2+Mn5z4VhumDXM7nBUJ6T1jeatH53JyuzDuNxu6l1uXly9n7veyOLDn8w+adpx1XMEMkGsAzJEZBhQCCwErmp6gFV3eAJYYIw51GTXh8D/iEgf6/m5wN0BjFV1wKKVucT3iuAKP635rOyVktCLq6Z/3VU7PiWehYtW89ePd3F3kK21obpGwEYgGWNcwO14Pux3AK8aY7aJyAMicpF12B+BWOA1EdkkIkusc8uAB/EkmXXAA9Y2FST2llTx4faDXDNjMDFaZwhJM9L7ceW0NJ78LJetXhY7UqFPp9pQHXLvW1m8tr6Az++aq/WGEFZe08C8v6ykf6yTd24/U6frCEE61Ybyq5LKOl7fUMBlU1M0OYS4+F4RPHjxOLYfqOCJlXvsDkd1MU0Qqt1eWLWP+kY3N+lMqj3CgvED+dbEQfxl6S4+313S4nHGGJ77Yi/PfbG3C6NTgaQJQrVLZZ2LF1bvZ96Y5G6xzKfyj99fOoGMpDhuX/wV+WXVp+x3uw33/9927rO+Pt112IYolb9pglDt8uznezla3cDtur5AjxLjdPDEtVNpdBtufXEDNU3m03I1uvnF65t5btU+vn/GUDKSYvnl65s5Wl1vY8TKH/T2E+Wz8uoGFn2Wy/yxyUxMS7A7HNXFhibG8MjCSdz4/HpueXE94wbF0+h2s7Wwgi9zS/n5/JHcfs4IthVV8O3Hv+Det7by6FWTEdFpVrorTRDKZ09+lktlnYufzR9pdyjKJueMTuae88fwx4+yWbO3DEeYEBURzoMXj+PamUMBz/iJn8wbyR8/zGb+pmSdKbYb0wShfFJSWcczX+zlwtMGMWZgb7vDUTa6eXY6N89u/QaF284ezvKdh/j121uZlJbA0MSYLopO+ZPWIJRP/rliD7UNjfxkXobdoahuIDxM+N8rJhEeLtz0wnoqahvsDkl1gCYI1aaD5bW8uHo/l05J1TuXlM/S+kbz+NVT2FdSxZ2LN9LoDo1BuT2JJgjVKmMM//12Fga48xvaelDtc8bwRO6/eBzLsw/z8Ac77Q5HtZPWIFSrFq/N5+Mdh/j1hWNJ6xttdziqG7p6+hB2HTzGok9zyT1cxegBcYxIimVIv2j6xTjpExNBrNOhdzsFIU0QqkV7S6p48N3tzBqRyPVnDLU7HNWN/frCsbgNfJFTwvLsQ6d0N8VEhnPnvAxumpWuqw8GEU0QyquGRjc/eWUTkY4w/vSdifqfVnWKIzyMBy8ZD0C9y83+0iryj1RTWllPWVU9a/aW8T/v7WTVnlL+/J2J9It12hyxAk0QqgWPL9/D5vyjPHbVFAbE64R8yn8iHWFkJMeRkRx3Ytsts9P51+r9PPifHZz/yGdcNX0wgtBoDP1jI7lmxhDtgrKBJgh1imO1DTz5WS4XTBjAN08baHc4qgcQEa6dOZSpQ/py58sb+evHu0/an9o3mrmjkmyKrufSBKFO8er6AirrXNx29nC7Q1E9zNhBvfnop7NxuQ1hIjS6DXP/tIJHPt7NnJH9T2pFuN0Gg2fMRXdW29DIvtIqhvaLISoiuJZ21QShTtLoNjy/ah+ZQ/pwWqrOt6S6nogQEe750A8PE344dzj3vrWVz3aXMHtkf8Dze3rT8+uoqHXxxg/OsDNcnxljKDhSw67iY2QXH2PXwWNsP1DBnsNVnkkQz07n7vODa2lXTRDqJJ/sKCavrJq7zh9tdyhKAXD51FQeXZbDI5/s5qyMRESEvy/bzfJsz5Ti+WXV3eIW7HveymLx2vwTzwfGRzF2YG/OGzeAdzYVsb2owsbovNMEoU7yzBd7SUnoxbljk+0ORSkAnI5wfjhnOL9+Zxtf7iml0Rge+WQ3s0Yk8nlOCR/vKOb6M4fZHWablu88zMz0fvzivJFkJMfROyrixL68smrW7ztiY3Te6UhqdcK2onJW55bxvZlDcOjawyqIfCczjeTeTh56fyc/eXkTGUmxLPreVEYkxfLxjmK7w2vTgfIaDlbUcu64ZKYO6XtScgBIT4ylqLzmpHU2goF+CqgTnvtiH70iwll4+mC7Q1HqJFER4dx29nCyCsupbWjk8aunEh3pYN6YZNbkllFeE9yTAW7MOwrA5MF9vO5P7x+DMZ7BqcFEE4QCoOhoDe9sKuKyqSnER0e0fYJSXezKaYM5b1wyf104mRFJnkkj549NwuU2rAzyJU435R8l0hHG2Bamyk/v75kOPbeksivDapMmCGVNyLeV8DDh1tl6a6sKTlER4TxxbSbzm9THJqX1oV9MJB9vD+5upo15Rxg3qDeRDu8fucOs9TJyD2sLQgWZJZuLWLbzEL84b1S3uBtEqePCw4RzRiexPPsQDY3uE9sPVdRSWeeyMbKvNTS62VJQzuQ0791LANGRDgbFR5F7WFsQKoiUVtZx35JtTEpL4Ps6IZ/qhuaNTeZYrYt1e8sAWJVTwpw/rWDWw8v4x4o9VNfbmyiyDx6jzuVm8uDWxxWl948lV2sQKpg88O52Kutc/OHy07r9iFTVM52VkUikI4ylO4pZvvMQ1z+3jtQ+vZiUlsDDH+xk9h+W89KaPNvi25jnuX11UlpbCSKG3MNVGBM8CyvpOIge7MNtB3lnUxE/mZfByCYTpynVnURHOpg1IpE3vyrkX6v3M2pAHC/cMJ2+MZGs31fGHz7I5p63sugbE8GC8V0/t9jGvKMkxjpJ7dOr1ePSE2OorHNx+FgdSb2DY4JMbUH0UFsKjvLTVzYxPqU3P5wzwu5wlOqUeWOSKa9p4LTUBF66eQZ9YyIByBzal3/dNJ3TUuP51RtZHCiv6fLYNuUfZVJaQpuz0Q637szaE0SFak0QPVBeaTU3PLeOPtGRPHPd6S3eWaFUd3HZ1BQevmwCL9ww7ZRBaJGOMB5ZOJmGRjc/fWVTl66NfaSqntySqjbrD+CpQUBw3eqqnww9TFlVPdc9uxaX2/D8DdOCpimrVGc4HeFccfpgYpzee82HJcZw/0XjWJ1bxj9X7umyuDYVHB8g13aCGNg7iqiIsKC61VVrED1InauRm19YT9HRGv590/QTg42U6gkun5rKyl2H+cvSXeSVVnNaWjwTUxMYPSAuYFPLbMo7igg+zYwcFiYMS4xlTxDd6qoJoocwxvCbt7exYf8RHr96CplD+9odklJdSkT4f9+eQKPb8OH2g7yy3jOz6rwxSTx13ekBec2N+UcZlRxHbAstm+bS+8eQVVAekFg6QhNED/GvNXm8sj6f2+eO4IIJukqc6pnie0Xwj2umYowhv6yG57/cx9Of72XD/jKmDvHvH01ut2Fz/lEumDDA53OGJ8bwftYB6lyNOB32Lx6kNYgeYO3eMu5fso1zRifx0/kj7Q5HKduJCIP7RfPzc0fSJzqCR5fl+P01Ptl5iPKaBmak9/P5nPT+sbgN7C+t9ns8HaEJIkRV1rlYtaeEx1fk8MN/byCtbzT/e8UkHQynVBPRkQ5unDWM5dmH2Vrov66dRrfhjx/uJD0xhm+2o8V+YtK+IKlDBDRBiMgCEckWkRwRucvL/tki8pWIuETk8mb7GkVkk/W1JJBxhhJXo5tfvraZCfd9yFVPruEPH2TTNyaSJ783lfheOkurUs1dO3MocU4Hj6/wXyvi7Y2F7Cqu5GfnjmxXAfz4pH3tHQvhanRTXt1AwZFqdh6sYN2+MoqOdn7MR8BqECISDjwGzAcKgHUissQYs73JYXnA94FfePkRNcaYSYGKLxQ1NLr5ySub+M+WA1w7YwjnjEliUmoCfaxBQ0qpU8X3iuC6M4by2Ioccg4dY0RS52YVqHM18peluxif0psL2jlyOy4qguTeTq+3umYVlPPS2v3sKq6kqs5FdX0jVXUuKutc1Lncpxw/LDGG5b+Y09HLAAJbpJ4G5BhjcgFE5GXgYuBEgjDG7LP2nXp1ql0aGt3csXgj7289yD0XjOYWnbZbKZ/dMGsYT3++l8eX7+EvV3Tu79KX1uRReLSGhy6dQFgHunTTE2PZebCCLQVHKa9pIK+smlfX5bO5oJxeEeFMHpxA35hoYp0OoiPDiXU6iLG+4qIcxDkdfJ5Twr/X5FFWVX9iVHlHBDJBpAD5TZ4XANPbcX6UiKwHXMDvjTFvNz9ARG4BbgEYPLjnroJWU9/InS9v5KPtxfz6wrHcOCv41+dVKpj0jYnk6umDeXbVPn527khS+3Rs2vvKOhePLsthZno/zspI7NDPyEiO5YUv93PRo198vS0plvsvGse3p6ScMlLcm4ToSP69Jo+swnLOHtm/Q3FAYBOEt9TZnjHug40xRSKSDiwTkSxjzElDII0xi4BFAJmZmcEzBWIXyi+r5tYXN7DjYAX3XzSO63TKbqU65IZZw3h21T5e/HI/d18wpt3n1zY08svXNlNaVc8vF4xqc+6lltx+zggmpibQu1cEfaIj6BsTybDEmHb9vPEpnpXrtuQfDdoEUQCkNXmeChT5erIxpsj6nisiK4DJQNeNke8GVuWU8KOXvsLlNjzz/dOZOyrJ7pCU6rYGJfRiwbgBLF6bx53zMoiO9P3j8UhVPTe/sJ71+49w7wVjmNLC2tO+SIqL4rKpqR0+Hzy1jPT+MWzp5J1ZgbyLaR2QISLDRCQSWAj4dDeSiPQREaf1OBE4kya1CwXLdhZz7TNrSYx1suT2WZoclPKD688cSkWtize/KvT5nLzSai77xyq2FJTz9ysnc/Ps9ABG6LvTUuI7PSo7YAnCGOMCbgc+BHYArxpjtonIAyJyEYCInC4iBcB3gCdEZJt1+hhgvYhsBpbjqUFogrA0ug3/895OhiXG8NaPzjxxa5xSqnOmDunD+JTePLdq30kL97yyLo8H392Ou9lMsMdqG7jyydWUVtXzr5um862Jg7o65BZNSE3gYEUthypqO/wzAjrVhjHmPeC9Ztt+0+TxOjxdT83PWwVMCGRs3dn/bS4i51Alj189xec5XpRSbRMRrj9jGD9/bTOf55RwVkZ/XlqTxz1vZQEwoHfUSS2E3727gwPlNbx22xlMHdLxbqVAOC01HoCswnK+0cFZm3UkdTfjanTzyCe7GT0gjgXjfJ/jRSnlmwsnDiQxNpJnv9jHWxsLuPftLM4ZncR545J5+IOdbM73TOG9bGcxr6zP57azhwddcgAYO7A3YQJbOtHNpAmim3lrYyF7S6r42fyRHbrHWinVOqcjnKumD2HZzkP84rUtzEzvx+NXT+EPl00kuXcUP168kfyyan71RhajB8Rx57wMu0P2KsbpYERSLFmdKFS3mSBEJFpEfi0iT1rPM0Tkwg6/ouqwhkY3f1u2mwkp8cwfm2x3OEqFrGtmDMbpCGNSWgJPfi+TqIhw4qMj+NuVkyg8WsMFj3zGkap6/vzdiUEx62pLJqQksKXg6En1lPbwpQXxLFAHzLSeFwC/69CrqU55fUMB+WU1/Gz+yA7fY62UaltSXBSf/PxsXrp5+kmr1E0d0pefzR/JsToXd3wjg3GD4m2Msm0T0+IpqaznQHnHCtW+VDiHG2OuEJErAYwxNaKfTl2utLKOP3+0i0lpCcwZ1fGBL0op37Q0mvoHZw9n1ohEJqQEd3IATsS4paCcQQm9aHQbHnx3OzPS+7LAh3mifGlB1ItIL6xR0CIyHE+LQnURYwy/eiOLipoGHrp0grYelLJRWJgwMS2hW9QAxwzsjSNMyCr0FNb/8MFOnlu1j4+2Fft0vi8tiN8CHwBpIvJvPIPWvt+xcFVHLF6bz8c7ivnvb45hzMDedoejlOomoiLCGZkcx5aCct7YUMATn+YCeJ391ZtWE4TVlbQTuBSYgWd+pTuNMSWdilr5LPdwJQ++u51ZIxK54UydhE8p1T6npcazZHMRa3LLmJnejyPV9dS5Gn06t9UuJuMpfb9tjCk1xvzHGPOuJoeuc+hYLXe+vAlnRBh//u7EbtGkVUoFlwmp8VTXNzIgPorHr55CjNNBbYMfWhCW1SJyujXqWXWBnEOVPPVZLm9+VYjL7eYf10wluYMjIZVSPdvcUUmcMbwf9100jj4xkTgdYT63IHxJEHOBW0VkP1CFp5vJGGNO63jIqiWL1+Zx95tZOB1hfPf0VG6alc5QnWtJKdVBgxJ68dLNM048dzrCqKxz+XSuLwni/A7GpTrg+VX7GDeoN8/fMI3EWKfd4SilQozTEU6dj11Mbd7maozZDyQA37K+Eqxtys/yy6rZefAYl0xK0eSglAoIZ4TvXUy+TLVxJ/BvIMn6+peI/LhTESqvlm733Jus02gopQIlyhHu1yL1jcB0Y0wVgIg8DHwJ/L3DESqvlm4vJiMpVmsOSqmA8WsLAk9RuulPa8T7etOqE45W17N2X5m2HpRSAeW5i8l/LYhngTUi8pb1/BLg6Q7GplqwIvswjW6jCUIpFVBOR7j/EoQx5i8isgKYhaflcL0xZmOnIlSnWLq9mP5xTiamJtgdilIqhDkdYTS6Da5GN47w1juR2kwQIjID2GaM+cp6Hici040xa/wTrqpzNbIi+xAXTRqko6WVUgEVFeFZv6LW5Sa2jQThSw3iH0Blk+dV1jblJ1/uKaWqvlG7l5RSAeeM8Hzs1zW0Xaj2qUhtmixHZIxx41vtQvlo6fZioiPDOWN4ot2hKKVCnNNhJQgf6hC+JIhcEblDRCKsrzuB3M6FqI5zuw0f7yhmdkb/E00/pZQKlONLpPorQdwGnAEU4lludDpwS8fDU01tzD9KcUUd543X7iWlVOB93YJou4vJl7uYDgELOx2V8ur9rANEhofxjTGaIJRSgfd1DcIPLQgR+YOI9La6lz4RkRIRuabzYSpjDO9vPchZGYn0joqwOxylVA8QZXUx1fqpSH2uMaYCuBBPF9NI4JediE9ZNheUU3i0hvMntL14uFJK+cOJFoSfahDH/7S9AFhsjCnrcGTqJO9nHSAiXJiv3UtKqS7SniK1L7er/p+I7ARqgB+KSH+gtjMBKk/30ntbD3DmiETio7V7SSnVNdpTpPZlPYi7gJlApjGmAagGLu5ciGpbUQX5ZTVcMF67l5RSXedEC8KHIrVPA96MMUeaPK7CM5padcJ7WQcIDxMdPa2U6lJRVg2i1h8tCOV/xhjeyzrAGcP70Scm0u5wlFI9SHtaEK0mCPFI809Y6rhtRRXsK63mfO1eUkp1Mb/dxWTNwfS2X6JSgKf18ND7O4hzOlgwfoDd4SilepjIcD8WqYHVInJ6J2NSlje+KuSLnFL+6/zR9NXuJaVUFwsLEyLDfVtVzpci9VzgNhHZh6c4LXgaF6d1KsoeqKSyjt/9ZzuZQ/pw9bTBdoejlOqhnI4wv42kPh9IB84BvoVnRPW3fAlCRBaISLaI5IjIXV72zxaRr0TEJSKXN9t3nYjstr6u8+X1gt2D726nqs7FQ5dO0IWBlFK2cUb4tuyoL+Mg9gNpwDnW42pfzhORcOAxPAlmLHCliIxtdlge8H3gpWbn9gV+i2fm2GnAb0WkT1uvGcw+3l7MO5uK+OGcEWQkx9kdjlKqB3M6wvwzDkJEfgtkAqOAZ/FMvfEv4Mw2Tp0G5Bhjcq2f8zKeAXbbjx9gjNln7Wse6XnA0uPTeojIUmABsLjNKwoyew5X8tePd/PuliIykmL54dzhdoeklOrhnBFh/pnuG/g2MBn4CsAYUyQivvwJnALkN3l+fC0JX3g7N6X5QSJyC9baFIMHB1effqPb8Ot3tvLy2jycjnBuO3s4t85OP3EPslJK2cXp8K2LyZcEUW+MMSJiAEQkxscYvHWyGy/bOnyuMWYRsAggMzPT15/dJTbmHeGlNXlckZnGLxeMIjHWaXdISikFWF1MfprN9VUReQJIEJGbgY+BJ304rwBP7eK4VKDIh/M6e25QWLnrMOFhwj3fHKPJQSkVVKIi/HQXkzHmT8DrwBt46hC/Mcb83YcY1gEZIjJMRCLxrKMsHGsAABDOSURBVEq3xIfzAD4EzhWRPlZx+lxrW7exctdhJqclEN9LZ2pVSgUXf3YxYYxZCixtTwDGGJeI3I7ngz0ceMYYs01EHgDWG2OWWAPw3gL6AN8SkfuNMeOMMWUi8iCeJAPwQHdah6Kkso4tBeX8fP5Iu0NRSqlTeO5i6kSRWkQ+N8bMEpFjnNz/f3ygXO+2frgx5j3gvWbbftPk8To83Ufezn0GeKat1whGn+8uAWD2yP42R6KUUqdyRoRT35kWhDFmlvVdb9pvp5W7DtM3JpIJKfF2h6KUUqfwS5FaRMJEZKvfouoB3G7Dp7sOc1ZGoo6WVkoFJb9MtWGMcQObRSS4BhkEsW1FFZRW1XO2di8ppYJUlI9TbfhSpB4IbBORtTRZSc4Yc1HHwwtdK3cdAuCsDE0QSqng5Oli8s9I6vs7H07PsXLXYcan9KZ/nI59UEoFJ6cjnIZGQ6O79fHFbSYIY8zK449FJBEotRYSUs2U1zTwVd5Rbjs73e5QlFKqRcdXlWvrTqYWaxAiMkNEVojImyIy2SpWbwWKRWSBP4MNFV/klNDoNszW7iWlVBBzOjwf/W0VqltrQTwK3APEA8uA840xq0VkNJ5ZVT/wS6QhoqHRzd+X5ZAU52TKkG49M7lSKsQdnzS0rUJ1a3cxOYwxHxljXgMOGmNWAxhjdvoryFDy1Gd72XGggvsvGkdEuC9TXCmllD2iInxbl7q1T7KmqaWm2T6tQTSxv7SKv368i/ljk1kwfoDd4SilVKt8bUG01sU0UUQq8Eyt0ct6jPU8yg8xhgRjDPe8lUVEeBgPXjweER0cp5QKbsdrEG2tKtfaVBu6so0P3viqkC9ySnnw4nEMiNe8qZQKfk4/dDGpNmQfPMb9S7YxdUgfrp4+xO5wlFLKJ8e7mGrbaEFoguigQxW13PDcOqIiw/nblZN13iWlVLfhjyK1akF1vYsbn19PWVU9z1x3OikJvewOSSmlfOaP21yVF41uw50vb2JbUTl/v3IyE1J1Sm+lVPdyokitLQj/WrK5kKXbi/nvb45l3thku8NRSql2O1Gk1hqE/xhjWPTpXkYkxfL9M4baHY5SSnXI10VqbUH4zao9pew4UMFNs4ZpUVop1W193cWkLQi/efKzXBJjI7lkcordoSilVIdpgvCz3cXHWJF9mGtnDCUqQscQKqW6L0d4GI4w0SK1vzz12V6cjjCunakD4pRS3Z/TEaZFan84dKyWtzYWcvnUVPrGRNodjlJKdZrTh3WpNUH44KU1eTS43dw4a5jdoSillF84HWF6F5M/rNpTysTUBNL7x9odilJK+UWUtiA6z+02bC+q4DQdMa2UCiFOR5gWqTtrX2kVlXUuxg/SBKGUCh2eBKEtiE7ZWuRZJ2lcSm+bI1FKKf9xOsL1LqbO2lZYTmR4GBlJcXaHopRSfuOMCKNWu5g6Z2tROaMGxBHp0H8qpVTo0HEQnWSMYWthBeO1e0kpFWI84yC0BdFhBUdqKK9pYJwWqJVSIUaL1J20ragcgAkpmiCUUqHF6dBxEJ2ytbCC8DBh1AAtUCulQouOpO6krUXlZCTF6uytSqmQ44ywuYtJRBaISLaI5IjIXV72O0XkFWv/GhEZam0fKiI1IrLJ+vpnIOP0xlOgLme8di8ppUKQ0xFOfRsJwhGoFxeRcOAxYD5QAKwTkSXGmO1NDrsROGKMGSEiC4GHgSusfXuMMZMCFV9bDh2ro6SynvGD9A4mpVToiYpou30QyBbENCDHGJNrjKkHXgYubnbMxcDz1uPXgW+ISFCs5bm10FOg1haEUioUHV+XujWBTBApQH6T5wXWNq/HGGNcQDnQz9o3TEQ2ishKETnL2wuIyC0isl5E1h8+fNivwWcVliMCYwZqC0IpFXqcPgz+DWSC8NYSMD4ecwAYbIyZDPwMeElETvmkNsYsMsZkGmMy+/fv3+mAm9paWMHw/rHEOAPWC6eUUraxO0EUAGlNnqcCRS0dIyIOIB4oM8bUGWNKAYwxG4A9wMgAxnqKbUXlWn9QSoUspw93ZwYyQawDMkRkmIhEAguBJc2OWQJcZz2+HFhmjDEi0t8qciMi6UAGkBvAWE+SX1bNgfJaJqYldNVLKqVUl4ryoQURsP4TY4xLRG4HPgTCgWeMMdtE5AFgvTFmCfA08KKI5ABleJIIwGzgARFxAY3AbcaYskDF2tyK7EMAnD3Sv91WSikVLHxpQQS0g90Y8x7wXrNtv2nyuBb4jpfz3gDeCGRsrVmefZgh/aIZlhhjVwhKKRVQdtcguqXahkZW7Slh7qgkguSOW6WU8jtNEB2wOreU2gY3c0Zp95JSKnTZPQ6iW1qRfZioiDBmpPdr+2CllOqmnDaPpO52jDEs23mIM4Yn6gR9SqmQ5stnnCaIJvaWVJFXVs1c7V5SSoU4rUG00/Jsz3Qdc0Yl2RyJUkoFliaIdlqRfYgRSbGk9Y22OxSllAooLVK3Q1WdizW5ZczRwXFKqR4gIlxo605+TRCWL/eUUt/oZu5o7V5SSoU+ESGqjVaEJgjLhrwjOMKEqUP62B2KUkp1ibZuddUEYdlaWM7I5Di9vVUp1WO0VajWBIFn/ENWYTkTdPU4pVQP0lahWhMEUHCkhqPVDUxI1QShlOo5tAXhgyxr/WltQSilehKtQfggq7AcR5gwakCc3aEopVSX0buYfKAFaqVUT6QtiDYcL1CfpvUHpVQPo0XqNhwvUI/X+oNSqofRInUbtECtlOqpNEG0IauwnIhwYfRALVArpXqWtuquPT5BHC9Q+zKzoVJKhRJtQbTCGMOWAh1BrZTqmZzagmhZwZEaymu0QK2U6pm0BdGK4wVqvcVVKdUTzRzer9X9PT5BRITrCGqlVM90xvDEVvf32ARRVlXPx9uLtUCtlFIt6JEJYseBCi569HP2l1Vzxzcy7A5HKaWCksPuALqC222obmikstbF2n1l3PXGFuKiHLx660wmpSXYHZ5SSgWlkE8QL3y5j/uWbMNtvt42MS2BRddOJbl3lG1xKaVUsAvpBGGM4ZnP95KRFMflU1OJjXLQJzqCOaOSdOZWpZRqQ0gniK/yjrCvtJo/fWcil09NtTscpZTqVkK6SP36hkJ6RYSzYPwAu0NRSqluJ2QTRG1DI+9uKeL88QOIdYZ0Q0kppQIiZBPExzuKOVbr4jLtWlJKqQ4J2QTxxoYCBsZHMSO99aHkSimlvAtoghCRBSKSLSI5InKXl/1OEXnF2r9GRIY22Xe3tT1bRM5rz+seOlbLp7tL+PbkFMLDpPMXopRSPVDAEoSIhAOPAecDY4ErRWRss8NuBI4YY0YA/ws8bJ07FlgIjAMWAI9bP88n72wsotFtuHSKdi8ppVRHBbJ6Ow3IMcbkAojIy8DFwPYmx1wM3Gc9fh14VETE2v6yMaYO2CsiOdbP+7KlF9tVfIz5f1kJwIHyWiamJTAiKda/V6SUUj1IIBNECpDf5HkBML2lY4wxLhEpB/pZ21c3Ozel+QuIyC3ALQC9B6WTkexJCCOT47h6xmD/XIVSSvVQgUwQ3jr/jY/H+HIuxphFwCKAzMxM8/jVU9sbo1JKqRYEskhdAKQ1eZ4KFLV0jIg4gHigzMdzlVJKBVAgE8Q6IENEholIJJ6i85JmxywBrrMeXw4sM8YYa/tC6y6nYUAGsDaAsSqllGomYF1MVk3hduBDIBx4xhizTUQeANYbY5YATwMvWkXoMjxJBOu4V/EUtF3Aj4wxjYGKVSml1KnE8wd795eZmWnWr19vdxhKKdWtiMgGY0ymt30hO5JaKaVU52iCUEop5ZUmCKWUUl5pglBKKeVVyBSpReQYkG13HH6SCJTYHYSfhMq1hMp1gF5LMLLzOoYYY/p72xFKK+lkt1SJ725EZL1eS3AJlesAvZZgFKzXoV1MSimlvNIEoZRSyqtQShCL7A7Aj/Ragk+oXAfotQSjoLyOkClSK6WU8q9QakEopZTyI00QSimlvAqJBCEiC0QkW0RyROQuu+NpLxHZJyJZIrJJRNZb2/qKyFIR2W1972N3nM2JyDMickhEtjbZ5jVu8fib9R5tEZEp9kV+qhau5T4RKbTel00ickGTfXdb15ItIufZE/WpRCRNRJaLyA4R2SYid1rbu9370sq1dMf3JUpE1orIZuta7re2DxORNdb78oq1NALWUgevWNeyRkSG2hK4MaZbf+GZSnwPkA5EApuBsXbH1c5r2AckNtv2B+Au6/FdwMN2x+kl7tnAFGBrW3EDFwDv41ktcAawxu74fbiW+4BfeDl2rPV75gSGWb9/4XZfgxXbQGCK9TgO2GXF2+3el1aupTu+LwLEWo8jgDXWv/erwEJr+z+BH1iPfwj803q8EHjFjrhDoQUxDcgxxuQaY+qBl4GLbY7JHy4GnrcePw9cYmMsXhljPsWzjkdTLcV9MfCC8VgNJIjIwK6JtG0tXEtLLgZeNsbUGWP2Ajl4fg9tZ4w5YIz5ynp8DNiBZz33bve+tHItLQnm98UYYyqtpxHWlwHOAV63tjd/X46/X68D3xARb0sxB1QoJIgUIL/J8wJa/yUKRgb4SEQ2iMgt1rZkY8wB8PxHAZJsi659Woq7u75Pt1tdL8806ebrFtdidUtMxvPXard+X5pdC3TD90VEwkVkE3AIWIqnhXPUGOOyDmka74lrsfaXA/26NuLQSBDesmp3u3f3TGPMFOB84EciMtvugAKgO75P/wCGA5OAA8Cfre1Bfy0iEgu8AfzEGFPR2qFetgX7tXTL98UY02iMmQSk4mnZjPF2mPU9KK4lFBJEAZDW5HkqUGRTLB1ijCmyvh8C3sLzy1N8vKlvfT9kX4Tt0lLc3e59MsYUW/+p3cCTfN1dEdTXIiIReD5Q/22MedPa3C3fF2/X0l3fl+OMMUeBFXhqEAkicnxOvKbxnrgWa388vneB+k0oJIh1QIZ1N0AknoLOEptj8pmIxIhI3PHHwLnAVjzXcJ112HXAO/ZE2G4txb0E+J5118wMoPx4l0ewatYX/2087wt4rmWhdafJMCADWNvV8Xlj9VM/Dewwxvylya5u9760dC3d9H3pLyIJ1uNewDw8NZXlwOXWYc3fl+Pv1+XAMmNVrLuU3dV9f3zhuRNjF54+vXvtjqedsafjufNiM7DtePx4+hs/AXZb3/vaHauX2BfjaeI34PmL58aW4sbTZH7Meo+ygEy74/fhWl60Yt2C5z/swCbH32tdSzZwvt3xN4lrFp6uiC3AJuvrgu74vrRyLd3xfTkN2GjFvBX4jbU9HU8SywFeA5zW9ijreY61P92OuHWqDaWUUl6FQheTUkqpANAEoZRSyitNEEoppbzSBKGUUsorTRBKKaW8crR9iFKqKRE5fssowACgEThsPa82xpxhS2BK+Zne5qpUJ4jIfUClMeZPdseilL9pF5NSfiQildb3OSKyUkReFZFdIvJ7EbnaWhMgS0SGW8f1F5E3RGSd9XWmvVeg1Nc0QSgVOBOBO4EJwLXASGPMNOAp4MfWMY8A/2uMOR24zNqnVFDQGoRSgbPOWPMaicge4CNrexYw13o8DxjbZKr/3iISZzzrHyhlK00QSgVOXZPH7ibP3Xz9fy8MmGmMqenKwJTyhXYxKWWvj4Dbjz8RkUk2xqLUSTRBKGWvO4BMa3W07cBtdgek1HF6m6tSSimvtAWhlFLKK00QSimlvNIEoZRSyitNEEoppbzSBKGUUsorTRBKKaW80gShlFLKq/8PcwNKPba2IhEAAAAASUVORK5CYII=\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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5bnA8d+TySSTHbIQQhIIS9h3AiooVUu9uGLrhtpWrb3cWulq+6m99Vpr721tvdXqLa1Va7Va9xWVitUibiib7JCVJQskISH7OjPv/WMmGEI2kgxnluf7+eSTOee8c/IcTsgz73LeV4wxKKWUCl1hVgeglFLKWpoIlFIqxGkiUEqpEKeJQCmlQpwmAqWUCnHhVgdwqpKTk01WVpbVYSilVEDZsmXLUWNMSnfHAi4RZGVlsXnzZqvDUEqpgCIiB3s6pk1DSikV4jQRKKVUiNNEoJRSIU4TgVJKhThNBEopFeI0ESilVIjTRKCUUiEu4J4jUGogGludFFY2UFjZwNH6Ni6dNYqRCQ6rw1LKL2giUEGpuLqJT4qq2HSgmk0HjrH/aOMJx//37VxuWjSWW74wnoRou0VRKuUfNBGooGCMYc/hOtbuLuft3UfYd6QegIQoO/OzhnPF3HQmjIhjwogYbGFhPPhuPn9+v5BnNh7imvmZnD95BPPGDMdu09ZSFXok0FYoy8nJMTrFhOpsT1kd/7NmDx8VVCEC88ckcsG0VBZPTGFCSixhYdLj++77Zx7r8ypodxniHOGcO2kEV8xN55zsFGw9vE+pQCQiW4wxOd0e00SgAkldSzu1Te20udw0t7l46pODPLe5mIQoOyvPm8Dlc9JJjo08pXPWt7TzUcFR/rWvgrf3lFPT1M7IeAdfmZvO2dnJzEhPIM7haT5yutzsP9pIY5uLWRkJiGiyUIFBE4EKCnnl9Vz84Ae0uz7/nQ0PE25YmMV3z88ekrb+VqeLd/dW8MLmYtbnVeI2IALjkmOICLdRWNFAm8sNwIrF4/jphZM1GaiA0Fsi0D4CFTAeeb8IW5jwy2XTcdhtRISHMSM9gczE6CH7GZHhNi6akcZFM9I41tjG9pIadpTUsqOkhnaXYXF2MpNGxrHl4DEefr+I5jYXv7hsWo/NT0oFAk0EKiBU1Lfw2rYyrpmfyfIFo0/LzxweE8G5k0Zw7qQRJx378px0Yh3h/Hl9Ec3tLn5zxcw++xTcbkNDm5OWNhdNbS5qm9sprWmm5FgTZTUtxETaGBnvIDXeQWZiNONTYokI185r5XuaCFRAeGrDQdrdbm5alGV1KACICLcvnUy0PZz738nj7d1HmJ6ewIz0BMYkxQBgMLS2u8mvaGDP4Tpyj9TR0u7u9nyxkeE0t7twuU9s9hqXEkNWUgy1ze1U1rdS2dDKqIQoFk5IYtH4ZM4cn0RspP43VoOjfQTK77W0uzjr1+8yb0wij97QbROnpd7adYT38yvZWVLLviN1J/RhgGcI65S0OKamJTBqmIOoCBtRdhtxDjvpw6JIHx5FQpQdl9tQ1dDKkboW9h9tJPdIPblH6ik+1sSwqAhS4iNJjomg6GgjG/dX0+p0kxofydrvL2ZYdIRFV68ChfYRqID28tZSjjW1881zxlodSreWTh/J0ukjAU9nc1VDG2EihAnYwoTEmIh+dSjbwoQR8Q5GxDuYmTGs17It7S7ez6vkW09t4ffv5HPXZdOG5FpUaNIGSGWJosoGPsivpN11YlNJfUs7r20rZV1uBc1tLtxuw18+LGJGegJnjE20KNr+iwy3MWpYFCMTPH/Qk2IjfTKqyGG3ccG0kSxfMJonPzlIfnn9kP8MFTp8WiMQkaXAA4ANeNQYc083Za4G7gIMsN0Yc50vY1LWO1jVyJUPbaC6sY3h0XYunJFGzpjh/GtfBf/cU06r05McImxhTEmLo7CykQeWz9Zhmt247UsTeX17Gb98cy9P3DRf/43UgPgsEYiIDVgFfAkoATaJyGpjzJ5OZbKBnwKLjDHHROTk4RkqqBxrbOOmv27CGMO9V87k/fyjvLK1lKc/PcTwaDtX52Ry+ZxRNLa6+CC/kvfzjjI1LZ6LZqRZHbpfSoqN5HtfzOa/39zLutwKzp+canVIKgD5rLNYRM4C7jLG/Jt3+6cAxphfdyrzWyDPGPNof8+rncWBq9Xp4muPbmRbcQ1///czmJ/laeppanOSV97A1LR4HS45AG1ON0sfeB9j4NVvL9JJ9FS3rOosTgeKO22XAGd0KTMRQEQ+wtN8dJcx5q2uJxKRFcAKgNGjT88YcjUwO0tqEYFpo+JPaKaoqGvhF2/sYeOBah68ds7xJAAQHRHO7MzeO0dVzyLCw/ivS6Zy0183MevutxmdGM20UfGMTHBgEzneCX3DWWMI10n1VDd8mQi6a6zsWv0IB7KBc4EM4AMRmW6MqTnhTcY8DDwMnhrB0IeqBsMYw4aiKh58N59PiqoByEqK5tJZoxifEsvq7WWsz6vE5Tb8ZOlkLps1yuKIg895k0bwwrfOYuP+anaX1bK7rI4P84/iMga3MbS0uymsbOB/Lp+u/QjqJL5MBCVAZqftDKCsmzKfGGPagf0ikosnMWzyYVxqiDS3uVi7+whPfXKQzQePkRofyZ2XTCU6wsbrO8pYta4At4GR8Q7+Y/E4rpiXwfiUWKvDDlrzsxJPqGl19pu39vGn9woZmxTDvy8ed5ojU/7Ol4lgE5AtImOBUmA50HVE0KvAtcDjIpKMp6moyIcxqQE6UttCfkU99S1O6lva+exQDW/sOExDq5PMxCjuXjaNq3MycdhtACxfMJrK+lYOVTcxO3OYTulssR9fMIlD1U386h97yUyMYul07XxXn/NZIjDGOEVkJbAWT/v/Y8aY3SJyN7DZGLPae+wCEdkDuIAfG2OqfBWTGpiK+haW3Leehlbn8X3REZ7J2a6cl8GCrMRuJ11LiYskJe7UpoRWvhEWJvzuqlkcrmnm+89tY9m+Sk+zkdswMyOBGxf558N66vTQKSZUn361Zi+PflDEI1/PIX14FHEOO0kxEcc//avAUdXQyn88uYXiY02Eh4XR0u7iWFMbH91+PmkJUVaHp3xIp5hQA1bd2MZTnxzk0lmj+OIUHaMe6JJiI3nxloXHtw9VNbH43nU8v6mE7y3JtjAyZSUdS6Z69deP9tPU5uLW8yZYHYrygdFJ0ZyTnczzm4tPmPlUhRZNBKpHdS3tPP7xAZZOG8nE1Dirw1E+cu2C0ZTWNPN+fqXVoSiLaCJQPXpyw0HqW5ysPF9rA8FsyZRUkmMjeObTQ1aHoiyiiUB1q6nNyaMfFHHepBSmpydYHY7yoYjwMK6cl8m7+yoor2uxOhxlAU0Eqlt/XFfIsaZ2rQ2EiOXzM3G5DS9sLu67sAo6mgjUSXaU1PCn9YVcMTeDeWP8fw0ANXhZyTEsmpDEMxuLcWunccjRRKBO0Op08aMXtpMcG8Gdl061Ohx1GnV0Gj+ntYKQo4lAneCBd/LJK2/gnitmkhCl0xmHkgunp3FOdjI/X72bnSW1vZbNPVJ/0upyKnBpIlDHbSuu4aH1hVydk8F5k3SNoFBjCxMeWD6HlNhIvvXUFo41tp1UprnNxe0v7eDffv8+P35hO4E2M4HqniYCBXjmE/rOM1tJjXdwxyXaJBSqEmMi+OP1c6msb+V7z2074SGzvPJ6lq36kOc2F3PWuCRe3VbGox/stzBaNVR0iglFQ6uTbzy+iaP1bTyz4kziHdokFMpmZQ7jrsum8Z+v7GTJfeuJibRhEyG3vJ7YyHD+9o0FnD0hmVuf3sqv/7GXSSPjWDwxxeqw1SBoIghx7S43tzy1hb2H63nk6/N0pTAFwLULMmlobWfzgWO43Aan23DxjFH85MJJjIhzAHDvlbMoqmzkO898xuqVixiTFGNx1GqgdPbREGaM4bYXtvPy1lJ+c8UMrpmvy4CqU1Nc3cSlf/iQlNhIXvr2Qq1N+rHeZh/VPoIQ9q99Fby8tZTvfjFbk4AakMzEaP54/Vz2H21k5dOf4dSRRAFJE0GIMsZw3z/zGJ0YzXf06WE1CAvHJ/Pfl0/n/bxKfvnGHqvDUQOgfQQh6u095ewuq+N/r5qF3aafB9TgLF8wmqKjjTz8fhFjk2O4YWEWIro8aaDQRBCC3G7D79/JZ2xyDJfPHmV1OCpI/GTpZIoqG7nr9T3c/04+41NiGJ8SS8bwaFLjI0mNdzBxZBzpw3QlNH+jiSAErd19hL2H6/j9NbMJ19qAGiK2MOH/rp3Di1uK2XeknsLKBt7Lq6SyvvWEMtfMz+T7S7KPjz5S1tNEEGLcbsP97+QxPiWGS2dpbUANragIG187K+uEfW1ON5UNrRypbeH17WU89clBXv2slJsWZTEmMYY2lxuny82Z45OYPDLemsBDnCaCEPP6jjLyyhv4v2vnYAvTNlzlexHhYaQPiyJ9WBTzxgznxoVZ3Pt2LqvWFZ5QbmxyDO/88Av6e2kBnyYCEVkKPADYgEeNMfd0OX4jcC9Q6t31B2PMo76MKZRV1LVw1+rdzMxI4OIZaVaHo0JUVnIMq66by88vbaHdZbDbhA/yjnLbC9t5e/cRLtTfzdPOZw3EImIDVgEXAlOBa0Wku0lsnjPGzPZ+aRLwEWMMP3pxB83tLu67ejZh+qlLWWxEnIP0YVGMiHNw+Zx0xiRF89D6Qp3IzgK+7ClcABQYY4qMMW3As8AyH/481Yu/bTjI+3mV/OyiKUwYEWt1OEqdwBYm/Ps549heUssnRdXH9ztdbv624QCHa5utCy4E+DIRpAOdV7go8e7r6goR2SEiL4pIpg/jCVn55fX8as1ezp2UwlfPHGN1OEp168p5GSTHRvDQek/fgcvtmQLlztd2c+/aXIujC26+TATdtT10rfO9DmQZY2YC7wBPdHsikRUisllENldWVg5xmMGtqqGVlU9/RkxkOL+9cqY+5KP8lsNu48aFWazPq2RXaS0/emE7r20rY3RiNP/YeYSGVqfVIQYtXyaCEqDzJ/wMoKxzAWNMlTGmY5DxI8C87k5kjHnYGJNjjMlJSdHpbvurvK6Fax7+hANVjTy4fI6O21Z+72tnZhETYeP6Rz/llc9K+dEFE7n/mtk0t7tYs+Ow1eEFLV8mgk1AtoiMFZEIYDmwunMBEek8POAyYK8P4wkpxdVNXPXQBg7XNPPENxZwdnay1SEp1aeEaDvXLhhNbXM7P1gykZXnZzN39DDGJcfw4pYSq8MLWj4bPmqMcYrISmAtnuGjjxljdovI3cBmY8xq4LsichngBKqBG30VTyg5cLSRax/5hMZWJ0998wzmjB5udUhK9duPl07iwhkjmTcmEQAR4Yp5Gdy7NpcDRxvJStZ1D4aarkcQZMpqmrnqoQ00tTn5+zfPZOoofVJTBb4jtS0svOddVp43gR9eMAnwDIkurWkmNd7hlxMnut2GstpmCisbyTtST255PbXN7dx39SziLFi3obf1CPTJ4iBSWd/KVx/9lLrmdp5ZoUlABY+RCQ7Ozk7hpa2lfH/JRFqcLn768k5e21ZGlN3GnNHDWDA2kWvmZ5KWYN2kdk6Xm5++vJNtxTUcrG6izfn5+gwJUXZqm9vZUFjFBdNGWhZjdzQRBAGX23CkroWbH9/E4doWnrx5AdPTE6wOS6khdeW8DL77zGc8vfEQT244SF5FPSsWj6PN6WbTgWoeeDefNTsP8/p3ziYy3GZJjFsOHuOFLSWcOS6R8yaPICsphnEpMUxKjSMqwsb0n69le0mNJgI1NIoqG/jJSzsoqGigprkdYyDCFsZfbswhJyvR6vCUGnIXTE0lzhHOHa/uYni0nSduWsDiiZ+PIlyXW8FNf93E/f/M5/YLJ1sS47rcSsLDhIe/ntPtsp0TU+PYUVJrQWS900QQgD7Ir+TWv28l3BbGxTPTSIyJJCkmgnljhmtNQAUth93Gt74wng2FVdxzxQwyhkefcPy8SSO4JieTh98v5IJpqcy1YJDEun0VzM9K7HHt5lmZw3hzRxnGGL96pkcTQQAxxvD4xwf47zf3kj0ilke+nkNmYnTfb1QqSNx63gRuPa/npVXvuGQKHxYc5UfPb+fN755DVMTpayIqrWkmt7ye/7yo59rIrIwEntl4iANVTYz1o9FP/tfVrrrV0Orkh89v5xev7+H8ySN48ZaFmgSU6iLOYee3V86k6Ggjd7y6i437qyk51oTT5e77zYO0bl8FAOdPHtFjmZkZwwDYUVLj83hOhdYIAsDOklq+88xWDlU38YMlE/nO+RN09lClerBoQjLfWDSWxz7az0tbPQ+hOexhPPL1HM7J9t3MBO/lVpCZGMX4lJ4ndZyYGovDHsb24lqWze5u6jVraCLwc89uPMR/vbaL5NhInl1xFgvGakewUn35r0umcP2Zoyk91kxpTTOr1hXwm7f2cfaEZJ+0zbe0u/iooIqrcjJ6PX+4LYzpoxLYrjUC1V8lx5q487XdLBibyB+uncvwmAirQ1IqIIgI41Nij386DxP4yUs7WZ9XybmTem66GahP91fT3O7ivF6ahTrMzBjG0xsP4nS5/WbNcP+IQnXrvrfzEIF7r5ylSUCpQfjynAxGJThYta7AJ+dft68Chz2Ms8Yl9Vl2VmYCLe1u8sobfBLLQGgi8FN7yup4ZVspNy7KYtQw656UVCoYRISHsWLxODYdOManRVVDem5jDP/aV8HC8ck47H2PUprl7TD2p+YhTQR+6p639hHvsPPtL/Q8VE4p1X/LF4wmOTaCPwxxraDoaCOHqpv61SwEMCYpmoQo+4BHDrW73OSV17N6exkPvptPcXXTgM7TmfYR+KEP848eX1YyIfr0T06lVDBy2G3cfPY4fvPWPrYX1zArc9iQnPf5zZ6FGM+b1L8RSSLCzIwEthf3/oRxbVM7B6oaOVDVSGFlI4UVDRRUNLD/aCNtnYbD1re087OLu1sOvv80EfgZt9twz1t7SR8WxdfO0mUllRpKXz1zNH96r4A/vVfIQ1/rdh2sU7JxfzWPvF/EVfMyTnrSuTczMxJ4aH0RzW2uEx56c7sNz24q5oF38yivaz2+P0wgMzGaCSmxnDs5hckj45g8Mp7bnt/O3sP1g74OTQR+5o/vFbCrtI77r5nVr/ZGpVT/xTnsfO2sMfzpvUKKq5sG9VBmbXM7P3huG6MTo7nrsmmn9N5ZGcNwuQ17DtceX3dhV2ktd7y6i23FNczPGs7NZ48lKymGrOQYRidGd/v3YNqoeNblVgz4GjpoIvAj63Ir+N0/81g2exSX+9HDJkoFk6+eOYaH1hfx5CcH+c+LpgzoHMYY7nh1F0fqWnjploXERJ7an9KOZqnf/COX+Khwyuta2V1Wy/DoCO67ehZfnpPer+cdpqTF88KWEirqWwa1FK12FvuJg1WNfO+Zz5g8Mp57vqKLzCvlK2kJUSydPpJnNx6iqc05oHO8tLWU17eX8YMl2cweQF9DaryDnDHDKaxsoLSmhcSYCFYsHs+/bjuXr8zt/aG0zianxQGwb5DNQ1oj8ANNbU7+48ktiAh//uq80zpRllKh6BuLsnhzx2Fe+ayU68/of19cY6uTe9fm8sSGAyzISuSWcwc+qu/FWxYO+L0dpoz0LD6193DdCVNynypNBBbZuL+aNTsPs7O0lt1ltbQ63Tx+0wJGJ+lEckr52tzRw5mRnsDjHx3gugWjERHqW9p58N18lkxJ5YxuHgxbn1fJf768k7LaZr5+5hh+vHQyNovn/BoeE8HIeAd7D9cN6jyaCCxQ09TG9Y9+QnhYGNPT47luwRi+OGUEiyYkWx2aUiFBRLhxYRa3vbCdjwqqSIqN4Nt/38r+o408s7GYF285i8kjP1/q9a8f7ecXr+9hfEoML37rrOMdvP5gSloc+45o01DAeXt3Oe0uw0u3nHl8Wlql1Ol1yaw0fv2Pvfx89S5Ka5qJd9hZdd1c7n5jNzc/vplXbl3IiDgHf9twgF+8vod/m5bKA8vn+N1ovilp8XyQf5RWp2vAS3RqZ7EF3tx5mMzEKGboamJKWSYy3MZ1Z4yhsLKROZnDefO753DxzDT+csN8qhvb+OYTm/nLh/u587XdfGlqKv937Vy/SwIAk9PicboNhRWNAz6HTxOBiCwVkVwRKRCR23spd6WIGBHJ8WU8/uBYYxsfFRzlohlpOjJIKYt9+9zxPPL1HJ68eQEpcZEATE9P4MFr57CztJZfvrGHJVNGsOq6uUSE++fn5qnekUOD6SfwWdOQiNiAVcCXgBJgk4isNsbs6VIuDvgu8KmvYvEnb+85gtNtuGTGKKtDUSrkOew2vjQ19aT9X5qaym+vmMm24hruvHSq3yYBgKykGCLDwwaVCHx5dQuAAmNMkTGmDXgWWNZNuV8CvwVafBiL33hz5xFGJ0YzPT2+78JKKctclZPJ/3x5xoDb3U+XcFsYE1MH12Hsy0SQDhR32i7x7jtOROYAmcaYN3o7kYisEJHNIrK5srJy6CM9TTqahS6eqc1CSqmhMyUtjr2H6zDGHN9XXtf/z9a+TATd/aU7HqWIhAH3A7f1dSJjzMPGmBxjTE5Kiu/WHPW1tbuP4HIbLp6RZnUoSqkgMiUtnqrGNiobPBPV/fWj/Zzxq3fZU9a/5iJfDh8tATI7bWcAZZ2244DpwHveT8cjgdUicpkxZrMP47LMmzsPMyYpmmmjtFlIKTV0Jh9/wriegooG/vvNvQDHE0NffJkINgHZIjIWKAWWA9d1HDTG1ALHn6ASkfeAHwVrEqhubOPjwir+Y/E4bRZSSg2pqWmeRPDu3nLe2HGY2MhwapvbaW5z9ev9PmsaMsY4gZXAWmAv8LwxZreI3C0il/nq5/qrNTsP43IbLtJmIaXUEEuItjMqwcHfNhyk3enmd1fNAqDV2b9E4NMni40xa4A1Xfbd2UPZc30Zi9Ve+ayU7BGx2iyklPKJKWnxHK5r4ffLZzPV+3fG8hqB+tzBqka2HDx2StPLKqXUqbjtgkk89NV5fHFKKg7vkNfmdj+oESiPl7eWIgKXz9GHyJRSvjF1VPzxmkDHVPYt7e7e3nKc1gh8zBjDq9tKWTg+ibSEKKvDUUqFgEjvk9D9rRFoIvCxrYeOcbCqiS/PybA6FKVUiBARHPYwWjQR+IeXtpbisIexdPpIq0NRSoWQKLtNE4E/aHW6eHPHYZZOG0nsKS5urZRSgxFlt+moIX+wbl8Ftc3tfHmuNgsppU4vh92mfQRWa2l38dD6IlLiIlk0/uT1T5VSypccdpvvRw2JyBUDfW+wc7sNP3x+m2cu80umEm7TfKuUOr2iIk5PH8H9g3hv0DLGcPcbe1iz8wh3XDyFS2fpswNKqdPPYQ87LU1D+ohsFy3tLv7wrwIe//gAN589lm+eM87qkJRSIepURg0NZiiL6btI8KttaueP7xWw8UA1u0praXcZLpmZxs8ummJ1aEqpEHYqncW9JgIR2Un3f/AFz/oBIe/+d/L424YDzBk9nG+cPZb5YxI5d1IKYWFaYVJKWcdht9HSz+GjfdUILhl8OMFtXW4FX5iYwl9vWmB1KEopdVyU3UaLs3+jhnpNBMaYgz0dE5GPgEWnFlpw2X+0kYNVTdx89lirQ1FKqRNERZyeB8pGD+K9QWHdvgoAzp04wuJIlFLqRI5wz6ihzgva92QwiSDkO4vX5VYwLiWG0UnRVoeilFIncHinom7tR/NQX53FX+npEBDScyo3tTn5dH81XztzjNWhKKXUSaLsHWsSuHB4X/ekr87iS3s59saphRVcNhRW0eZ0c94kbRZSSvmfjj/+ze0uhvVRtq/O4puGKqhg815uJdERNuaPHW51KEopdZLPawSDbBoCEJEvAMeMMTtE5GpgMVAI/NEY0zq4UAOTMYZ1uRUsHJ9MZHjvVS6llLKCw+5dpawfI4f66iNYBcwEHCKSC8QCbwELgceA6wcZa0AqrGyg5Fgzt5w73upQlFKqW52bhvrSV43gPGPMVBFxAKXACGOMS0T+DOzo6+QishR4ALABjxpj7uly/FvArYALaABWGGP29Bm1xd7LrQTgXO0fUEr5qY6modZ+JIK+ho+2ABhjWoCDxhiXd9sA7b29UURswCrgQmAqcK2ITO1S7GljzAxjzGzgt8B9fUbsB/65p5yJqbGkDwvpgVNKKT82lDWCESLyQzzDRTte491O6eO9C4ACY0wRgIg8CywDjn/iN8bUdSofQwA8m/DWriN8ur+anyydbHUoSinVo6iIoUsEjwBx3bwGeLSP96YDxZ22S4AzuhYSkVuBHwIRwPndnUhEVgArAEaPtu6B5mONbdzx6k6mjYrnm+fotBJKKf81ZKOGjDG/GEQc3U2/edInfmPMKmCViFwH3AHc0E2Zh4GHAXJyciyrNfx89W5qm9t58uYzsOuqY0opPxbZMWposDUCEbmzl8PGGPPLXo6XAJmdtjOAsl7KPwv8qbd4rPTWriOs3l7GD780kSlp8VaHo5RSvTpeI+jH8NG+PtY2dvMFcDPwkz7euwnIFpGxIhIBLAdWdy4gItmdNi8G8vuM2AJVDa3c8eoupo2K1yGjSqmA4Og0xURf+moa+l3HaxGJA74H3ITn0/vvenqf971OEVkJrMUzfPQxY8xuEbkb2GyMWQ2sFJEleEYgHaObZiGrOV1uVj79GXUt7Tx58wJtElJKBQS7LYzwMBmSzmJEJBFPZ+71wBPAXGPMsf4EYoxZA6zpsu/OTq+/15/zWOlXa/axoaiK3101S5uElFIBJaqfy1X21UdwL/AVPB21M4wxDUMTXmB45bMSHvtoPzcuzOKKeRlWh6OUUqfEEWHr16ihvto5bgNG4RnNUyYidd6vehGp6+O9AW1XaS23v7STM8Ym8rOLdSF6pVTgcdjDhqSPIGQbxO9+Yw/Dou2sun6u9gsopQJSlL1/y1XqX7hu7CmrY+P+ar559jiSYyOtDkcppQbEs4C9JoIBeeLjA0TZbVydk9l3YaWU8lORWiMYmGONbby6rZQvz00nIdpudThKKTVgUXZbv/oINBF08eymYlqdbm44K8vqUJRSalA8iWDwo4ZCitPl5skNB1g4PolJI+P6LK+UUv7MYaf6wz0AAA5uSURBVA/r13MEmgg6eWdvOWW1LdywMMvqUJRSatCiIvr3QJkmgk7++tEB0odFsWRKqtWhKKXUoDm0j+DU1LW08+n+aq7KycAW1t0M2kopFVg0EZyi/HLP7BnTRyVYHIlSSg2NKLuNdpfB6eq9w1gTgVd+eT0AE1O1k1gpFRwc3sVpWpyaCPolv6IBhz2MjOG6IL1SKjh0LE7T10Nlmgi88srrmTAiljDtH1BKBYn+Lk6jicArv7yBiSO0WUgpFTw0EZyCupZ2jtS1kK39A0qpIHK8aUgTQd86Rgxlj4i1OBKllBo6UREdNQLtLO6TjhhSSgWjjlFDWiPoBx0xpJQKRg4dNdR/OmJIKRWMOvoIWvtYnManiUBElopIrogUiMjt3Rz/oYjsEZEdIvKuiIzxZTw9KajQEUNKqeBjeY1ARGzAKuBCYCpwrYhM7VLsMyDHGDMTeBH4ra/i6UldSzuHa1uYkKodxUqp4OIPo4YWAAXGmCJjTBvwLLCscwFjzDpjTJN38xMgw4fxdKtjxJDWCJRSwcYfRg2lA8Wdtku8+3pyM/APH8bTrYIKHTGklApOkeH9GzUU7sMYuut5Nd0WFPkqkAN8oYfjK4AVAKNHjx6q+ADIK9cRQ0qp4CQiOOxhlj5ZXAJkdtrOAMq6FhKRJcDPgMuMMa3dncgY87AxJscYk5OSkjKkQeqIIaVUMOvPAva+TASbgGwRGSsiEcByYHXnAiIyB/gzniRQ4cNYeqQjhpRSwcxht1k3asgY4wRWAmuBvcDzxpjdInK3iFzmLXYvEAu8ICLbRGR1D6fzCR0xpJQKdlH2vtct9mUfAcaYNcCaLvvu7PR6iS9/fl8KKnTEkFIquHmWq9S5hnq093AdoCOGlFLBy+rOYr/3cWEVqfGRZCbqiCGlVHCKiui7aShkE4HbbdhQWMWiCcmI6IghpVRwsnrUkF/bd6Se6sY2Fo1PtjoUpZTymch+dBaHbCL4uPAoAIsmaCJQSgWvKLuNFp2GunsfFhxlXEoMIxMcVoeilFI+47CH0eLUUUMnaXO62bi/WpuFlFJBL8rKB8r82faSGpraXCyakGR1KEop5VP9eaAsJBPBRwVHEYGzxmmNQCkV3CK9axL0JmQTwYz0BBKi7VaHopRSPhWlieBkja1OPjtUw0LtH1BKhYCOxWl6E3KJYOOBapxuo/0DSqmQ4LD3/Wc+5BLBxwVHibCFkTMm0epQlFLK57RpqBubDhxj9uhh/aouKaVUoHNoIjiRMYb88nqmpsVbHYpSSp0Wmgi6KK1pprHNRbYuRKOUChHaNNRFfsdCNLr+gFIqROiooS7yy+sByB6hNQKlVGhwhGsiOEF+eQMpcZEMi46wOhSllDotHBE6fPQEeRUNWhtQSoUU7SPoxBhDQXm99g8opUKKjhrqpKy2hcY2FxO0RqCUCiF2WxjhYb0vx+vTRCAiS0UkV0QKROT2bo4vFpGtIuIUkSt9GUtHR7HWCJRSoaav5iGfJQIRsQGrgAuBqcC1IjK1S7FDwI3A076Ko0N+uWfoqPYRKKVCTV9TUYf78GcvAAqMMUUAIvIssAzY01HAGHPAe6z3ddSGQF55PcmxkQyP0RFDSqnQktTH3z1fNg2lA8Wdtku8+06ZiKwQkc0isrmysnJAweTriCGlVIha+4PFvR73ZSLornfCDORExpiHjTE5xpiclJSUgbyfgooGJurUEkopdRJfJoISILPTdgZQ5sOf16PDtS00tDqZoB3FSil1El8mgk1AtoiMFZEIYDmw2oc/r0d5HSOGtGlIKaVO4rNEYIxxAiuBtcBe4HljzG4RuVtELgMQkfkiUgJcBfxZRHb7IpYC72Rz2VojUEqpk/hy1BDGmDXAmi777uz0ehOeJiOf8owYiiBRRwwppdRJQuLJ4vyKBn2iWCmlehD0icAzx1AD2SO0WUgppboT9Ikgv6KB+lanDh1VSqkeBHUicLsN//XqLuIiw7lg2kirw1FKKb8U1Ing7xsP8en+au64ZAqp8Q6rw1FKKb8UtImg5FgT96zZy9kTkrk6J7PvNyilVIgKykRgjOGnL+/EAL/+ygxEep+LWymlQplPnyM4ndqcbnaX1fLZoRo+Lqzig/yj/HLZNDITo60OTSml/FrQJILL/vAh+454ppIYleDgpkVZXH/GGIujUkop/xcUiaCiroV9R+q5cWEW3/rCeEYmaMewUkr1V1D0EWwvqQXg0llpmgSUUuoUBUUi2FFSgy1MmJqWYHUoSikVcIIiEWwrrmFiahxREb2vy6mUUupkAZ8IjDHsKKlldqbWBpRSaiACPhEcrGqitrmdmRnDrA5FKaUCUsAngu0lNQDMzNAagVJKDUTAJ4IdJbU47GFM1NXHlFJqQAI+EWwvrmHaqATstoC/FKWUskRA//V0utzsKqtllvYPKKXUgAV0Isgrb6Cl3c0sHTGklFIDFtCJYMfxjmKtESil1EAFdCLYXlJLvCOcrCSdYVQppQbKp4lARJaKSK6IFIjI7d0cjxSR57zHPxWRrFM5//biGmZlDtP1BpRSahB8lghExAasAi4EpgLXisjULsVuBo4ZYyYA9wO/6e/5W9pd5JbX6/MDSik1SL6chnoBUGCMKQIQkWeBZcCeTmWWAXd5X78I/EFExBhjejppXnk9X7pvPe0uNy630RFDSik1SL5MBOlAcaftEuCMnsoYY5wiUgskAUc7FxKRFcAKgPhR48hOjQVgwdhEFk1I9knwSikVKnyZCLpruO/6Sb8/ZTDGPAw8DJCTk2P+eP28wUenlFIK8G1ncQmQ2Wk7AyjrqYyIhAMJQLUPY1JKKdWFLxPBJiBbRMaKSASwHFjdpcxq4Abv6yuBf/XWP6CUUmro+axpyNvmvxJYC9iAx4wxu0XkbmCzMWY18BfgSREpwFMTWO6reJRSSnXPp4vXG2PWAGu67Luz0+sW4CpfxqCUUqp3Af1ksVJKqcHTRKCUUiFOE4FSSoU4TQRKKRXiJNBGa4pIPZBrdRxDJJkuT1EHqGC5DtBr8UfBch1g7bWMMcakdHfAp6OGfCTXGJNjdRBDQUQ2B8O1BMt1gF6LPwqW6wD/vRZtGlJKqRCniUAppUJcICaCh60OYAgFy7UEy3WAXos/CpbrAD+9loDrLFZKKTW0ArFGoJRSaghpIlBKqRAXUIlARJaKSK53sfvbrY7nVIjIARHZKSLbRGSzd1+iiPxTRPK934dbHWd3ROQxEakQkV2d9nUbu3g86L1HO0RkrnWRn6yHa7lLREq992abiFzU6dhPvdeSKyL/Zk3UJxORTBFZJyJ7RWS3iHzPuz/g7ksv1xJQ90VEHCKyUUS2e6/jF979Y0XkU+89ec47LT8iEundLvAez7IseGNMQHzhmcq6EBgHRADbgalWx3UK8R8Akrvs+y1wu/f17cBvrI6zh9gXA3OBXX3FDlwE/APP6nNnAp9aHX8/ruUu4EfdlJ3q/T2LBMZ6f/9sVl+DN7Y0YK73dRyQ54034O5LL9cSUPfF+28b631tBz71/ls/Dyz37n8IuMX7+tvAQ97Xy4HnrIo9kGoEC4ACY0yRMaYNeBZYZnFMg7UMeML7+gngcgtj6ZEx5n1OXjmup9iXAX8zHp8Aw0Qk7fRE2rcerqUny4BnjTGtxpj9QAGe30PLGWMOG2O2el/XA3vxrAEecPell2vpiV/eF++/bYN30+79MsD5wIve/V3vSce9ehH4ooh0t3yvzwVSIji+0L1XCb3/svgbA7wtIltEZIV3X6ox5jB4/jMAIyyL7tT1FHug3qeV3iaTxzo10QXEtXibFObg+QQa0Pely7VAgN0XEbGJyDagAvgnntpKjTHG6S3SOdbj1+E9Xgsknd6IPQIpEfRroXs/tsgYMxe4ELhVRBZbHZCPBOJ9+hMwHpgNHAZ+593v99ciIrHAS8D3jTF1vRXtZp+/X0vA3RdjjMsYMxvPGu0LgCndFfN+95vrCKREcHyhe68MoMyiWE6ZMabM+70CeAXPL0l5R/Xc+73CughPWU+xB9x9MsaUe/8Du4FH+LyZwa+vRUTseP5w/t0Y87J3d0Del+6uJVDvC4AxpgZ4D08fwTAR6ZjXrXOsx6/DezyB/jdbDqlASgSbgGxvD3wEns6V1RbH1C8iEiMicR2vgQuAXXjiv8Fb7AbgNWsiHJCeYl8NfN07SuVMoLajqcJfdWkr/zKeewOea1nuHd0xFsgGNp7u+LrjbUv+C7DXGHNfp0MBd196upZAuy8ikiIiw7yvo4AlePo71gFXeot1vScd9+pK4F/G23N82lnd034qX3hGPuThaXf7mdXxnELc4/CMctgO7O6IHU974LtAvvd7otWx9hD/M3iq5u14PsXc3FPseKq7q7z3aCeQY3X8/biWJ72x7sDznzOtU/mfea8lF7jQ6vg7xXU2nmaEHcA279dFgXhfermWgLovwEzgM2+8u4A7vfvH4UlUBcALQKR3v8O7XeA9Ps6q2HWKCaWUCnGB1DSklFLKBzQRKKVUiNNEoJRSIU4TgVJKhThNBEopFeICcfF6pU4LEekYigkwEnABld7tJmPMQksCU2qI6fBRpfpBRO4CGowx/2t1LEoNNW0aUmoARKTB+/1cEVkvIs+LSJ6I3CMi13vnpd8pIuO95VJE5CUR2eT9WmTtFSj1OU0ESg3eLOB7wAzga8BEY8wC4FHgO94yDwD3G2PmA1d4jynlF7SPQKnB22S88/aISCHwtnf/TuA87+slwNRO083Hi0ic8cy/r5SlNBEoNXitnV67O227+fz/WBhwljGm+XQGplR/aNOQUqfH28DKjg0RmW1hLEqdQBOBUqfHd4Ec72pbe4BvWR2QUh10+KhSSoU4rREopVSI00SglFIhThOBUkqFOE0ESikV4jQRKKVUiNNEoJRSIU4TgVJKhbj/B89BMMDxIa/4AAAAAElFTkSuQmCC\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": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16404163283537904"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev.integrated_brier_score(time_grid) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4844648690929724"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ev.integrated_nbll(time_grid) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}