main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%reset -f\n",
    "%matplotlib inline\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import rpy2.robjects as ro\n",
    "import pylab as pl\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "from sklearn.model_selection import train_test_split\n",
    "from tick.preprocessing.features_binarizer import FeaturesBinarizer\n",
    "from lifelines.utils import concordance_index\n",
    "from lifelines import CoxPHFitter\n",
    "from QNEM.inference import QNEM\n",
    "from load_Data_forCox import load_Data_forCox\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "pd.set_option('display.max_columns', 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data\n",
    "dataFrame index on OfficialID, drop=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>gender</th>\n",
       "      <th>age</th>\n",
       "      <th>t_stage_mod_6is4_7is1</th>\n",
       "      <th>n_stage</th>\n",
       "      <th>UICC_TNM</th>\n",
       "      <th>postop_chemo</th>\n",
       "      <th>nlnode</th>\n",
       "      <th>plnode</th>\n",
       "      <th>sidedness</th>\n",
       "      <th>ascending</th>\n",
       "      <th>transverse</th>\n",
       "      <th>descending</th>\n",
       "      <th>sigmoid</th>\n",
       "      <th>differentiation</th>\n",
       "      <th>mucinous_colloide</th>\n",
       "      <th>occlusion</th>\n",
       "      <th>perforation</th>\n",
       "      <th>venous_emboli</th>\n",
       "      <th>lymphatic_invasion</th>\n",
       "      <th>perineural_invasion</th>\n",
       "      <th>msi_gen</th>\n",
       "      <th>p53_status</th>\n",
       "      <th>kras_status</th>\n",
       "      <th>apc_status</th>\n",
       "      <th>braf_status</th>\n",
       "      <th>pi3k_status</th>\n",
       "      <th>Immunoscore</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OfficialID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PAR170</th>\n",
       "      <td>0</td>\n",
       "      <td>34.08</td>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>III</td>\n",
       "      <td>1</td>\n",
       "      <td>43.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PAR171</th>\n",
       "      <td>0</td>\n",
       "      <td>84.48</td>\n",
       "      <td>4</td>\n",
       "      <td>0.0</td>\n",
       "      <td>II</td>\n",
       "      <td>0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PAR172</th>\n",
       "      <td>1</td>\n",
       "      <td>61.05</td>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "      <td>III</td>\n",
       "      <td>0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PAR173</th>\n",
       "      <td>0</td>\n",
       "      <td>75.51</td>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>III</td>\n",
       "      <td>0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PAR174</th>\n",
       "      <td>0</td>\n",
       "      <td>94.07</td>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>II</td>\n",
       "      <td>0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            gender    age  t_stage_mod_6is4_7is1  n_stage UICC_TNM  \\\n",
       "OfficialID                                                           \n",
       "PAR170           0  34.08                      3      1.0      III   \n",
       "PAR171           0  84.48                      4      0.0       II   \n",
       "PAR172           1  61.05                      4      2.0      III   \n",
       "PAR173           0  75.51                      3      1.0      III   \n",
       "PAR174           0  94.07                      3      0.0       II   \n",
       "\n",
       "            postop_chemo  nlnode  plnode  sidedness  ascending  transverse  \\\n",
       "OfficialID                                                                   \n",
       "PAR170                 1    43.0     3.0        1.0        1.0         0.0   \n",
       "PAR171                 0    13.0     0.0        1.0        1.0         0.0   \n",
       "PAR172                 0    21.0     4.0        2.0        0.0         0.0   \n",
       "PAR173                 0     7.0     1.0        2.0        0.0         0.0   \n",
       "PAR174                 0    15.0     0.0        1.0        1.0         0.0   \n",
       "\n",
       "            descending  sigmoid  differentiation  mucinous_colloide  \\\n",
       "OfficialID                                                            \n",
       "PAR170             0.0      0.0              1.0                0.0   \n",
       "PAR171             0.0      0.0              1.0                1.0   \n",
       "PAR172             0.0      1.0              1.0                0.0   \n",
       "PAR173             0.0      1.0              2.0                0.0   \n",
       "PAR174             0.0      0.0              2.0                0.0   \n",
       "\n",
       "            occlusion  perforation  venous_emboli  lymphatic_invasion  \\\n",
       "OfficialID                                                              \n",
       "PAR170              0            0            0.0                 0.0   \n",
       "PAR171              0            0            0.0                 0.0   \n",
       "PAR172              0            0            1.0                 0.0   \n",
       "PAR173              0            0            0.0                 0.0   \n",
       "PAR174              0            0            0.0                 0.0   \n",
       "\n",
       "            perineural_invasion  msi_gen  p53_status  kras_status  apc_status  \\\n",
       "OfficialID                                                                      \n",
       "PAR170                      0.0      0.0         0.0          0.0         0.0   \n",
       "PAR171                      0.0      0.0         0.0          0.0         0.0   \n",
       "PAR172                      0.0      0.0         0.0          0.0         0.0   \n",
       "PAR173                      0.0      0.0         0.0          1.0         0.0   \n",
       "PAR174                      0.0      0.0         0.0          0.0         0.0   \n",
       "\n",
       "            braf_status  pi3k_status  Immunoscore  \n",
       "OfficialID                                         \n",
       "PAR170              0.0          0.0            1  \n",
       "PAR171              0.0          0.0            0  \n",
       "PAR172              0.0          0.0            1  \n",
       "PAR173              0.0          0.0            0  \n",
       "PAR174              0.0          0.0            2  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X, Y, delta = load_Data_forCox()\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dividing in train and test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "70% for training, 30% for testing.\n"
     ]
    }
   ],
   "source": [
    "X.index = np.arange(X.shape[0])\n",
    "delta.index = np.arange(X.shape[0])\n",
    "Y.index = np.arange(X.shape[0])\n",
    "\n",
    "## Split data into training and test sets ##\n",
    "test_size = .3  # proportion of data used for testing\n",
    "rs = ShuffleSplit(n_splits=1, test_size=test_size, random_state=0)\n",
    "\n",
    "for train_index, test_index in rs.split(X):\n",
    "    X_test = X.loc[test_index]\n",
    "    delta_test = delta.loc[test_index]\n",
    "    Y_test = Y.loc[test_index]\n",
    "\n",
    "    X = X.loc[train_index]\n",
    "    Y = Y.loc[train_index]\n",
    "    delta = delta.loc[train_index]  \n",
    "    \n",
    "print(\"%d%% for training, %d%% for testing.\" \n",
    "      % ((1 - test_size) * 100, test_size * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrouver la p_value du test de Wald avec l'immunoscore univarié"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Immunoscore       int8\n",
       "Y              float64\n",
       "delta            int64\n",
       "dtype: object"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=X[['Immunoscore']].astype(np.int8)\n",
    "df[\"Y\"] = Y\n",
    "df[\"delta\"] = delta\n",
    "df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=207, number of events=42\n",
      "\n",
      "               coef  exp(coef)  se(coef)       z      p  lower 0.95  upper 0.95     \n",
      "Immunoscore -0.9072     0.4036    0.2599 -3.4913 0.0005     -1.4165     -0.3979  ***\n",
      "---\n",
      "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \n",
      "\n",
      "Concordance = 0.646\n"
     ]
    }
   ],
   "source": [
    "# Using Cox Proportional Hazards model\n",
    "cph = CoxPHFitter()\n",
    "cph.fit(df, duration_col='Y', event_col='delta')\n",
    "cph.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C-index in prediction on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C-index on test set: 0.59\n"
     ]
    }
   ],
   "source": [
    "X_test_immunoscore = X_test[['Immunoscore']].astype(np.int8)\n",
    "hat_beta = cph.summary.coef\n",
    "marker_test = X_test_immunoscore.dot(hat_beta)\n",
    "c_index = concordance_index(np.ravel(Y_test), marker_test, np.ravel(delta_test))\n",
    "c_index = max(c_index, 1 - c_index)\n",
    "print(\"C-index on test set: %.2f\" % c_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modèles de survie multivarié et comparaison du C-index en prédiction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X, Y, delta = load_Data_forCox()\n",
    "\n",
    "X.index = np.arange(X.shape[0])\n",
    "delta.index = np.arange(X.shape[0])\n",
    "Y.index = np.arange(X.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_cont = X[[\"age\", \"nlnode\", \"plnode\"]]\n",
    "scaler = StandardScaler()\n",
    "X_cont = pd.DataFrame(scaler.fit_transform(X_cont), columns=X_cont.columns)\n",
    "X.drop([\"age\", \"nlnode\", \"plnode\"], axis=1, inplace=True)\n",
    "\n",
    "binarizer = FeaturesBinarizer(remove_first=True)\n",
    "X_bin = binarizer.fit_transform(X)\n",
    "X_bin = X_bin.todense()\n",
    "blocks_start = binarizer.blocks_start\n",
    "blocks_length = binarizer.blocks_length\n",
    "\n",
    "columns_names = []\n",
    "for i in range(X.shape[1]):\n",
    "    for j in range(1, blocks_length[i]):\n",
    "        columns_names += [X.columns[i] + \"=\" + str(j)]\n",
    "        \n",
    "X_bin = pd.DataFrame(X_bin, columns=columns_names)\n",
    "X = pd.concat([X_cont,X_bin], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO: \n",
    "Je te laisse creuser un peu le code de binarization, ça te permettra de bien le maitriser et de bien comprendre ce qui est fait. Ce qu'il faut c'est retrouver les vrais valeurs des catégories plutôt que de mettre `str(j)` dans les columns_names ds le code au dessus, de façon à pouvoir interpréter correctement ensuite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "70% for training, 30% for testing.\n"
     ]
    }
   ],
   "source": [
    "## Split data into training and test sets ##\n",
    "test_size = .3  # proportion of data used for testing\n",
    "rs = ShuffleSplit(n_splits=1, test_size=test_size, random_state=0)\n",
    "\n",
    "for train_index, test_index in rs.split(X):\n",
    "    X_test = X.loc[test_index]\n",
    "    delta_test = delta.loc[test_index]\n",
    "    Y_test = Y.loc[test_index]\n",
    "\n",
    "    X = X.loc[train_index]\n",
    "    Y = Y.loc[train_index]\n",
    "    delta = delta.loc[train_index]  \n",
    "    \n",
    "print(\"%d%% for training, %d%% for testing.\" \n",
    "      % ((1 - test_size) * 100, test_size * 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=207, number of events=42\n",
      "\n",
      "                            coef          exp(coef)  se(coef)       z      p  lower 0.95  upper 0.95     \n",
      "age                       0.3541             1.4249    0.2207  1.6044 0.1086     -0.0785      0.7867     \n",
      "nlnode                   -0.4729             0.6232    0.2893 -1.6345 0.1022     -1.0399      0.0942     \n",
      "plnode                    0.7123             2.0386    0.2200  3.2374 0.0012      0.2811      1.1435   **\n",
      "gender=1                 -0.2924             0.7465    0.3829 -0.7637 0.4451     -1.0428      0.4580     \n",
      "t_stage_mod_6is4_7is1=1  13.4409        687567.9480  474.0308  0.0284 0.9774   -915.6425    942.5243     \n",
      "t_stage_mod_6is4_7is1=2  28.0834 1572044958939.0122 1546.7711  0.0182 0.9855  -3003.5323   3059.6991     \n",
      "t_stage_mod_6is4_7is1=3  28.6492 2768277332243.5752 1546.7712  0.0185 0.9852  -3002.9665   3060.2650     \n",
      "n_stage=1                 6.8977           990.0000       nan     nan    nan         nan         nan     \n",
      "n_stage=2                 5.4878           241.7247       nan     nan    nan         nan         nan     \n",
      "UICC_TNM=1              -14.9659             0.0000 1472.3437 -0.0102 0.9919  -2900.7066   2870.7748     \n",
      "UICC_TNM=2              -20.9524             0.0000       nan     nan    nan         nan         nan     \n",
      "postop_chemo=1            0.3523             1.4223    0.5614  0.6275 0.5304     -0.7481      1.4527     \n",
      "postop_chemo=2            0.7230             2.0605    0.7479  0.9667 0.3337     -0.7428      2.1887     \n",
      "sidedness=1              -4.6066             0.0100       nan     nan    nan         nan         nan     \n",
      "ascending=1              10.2283         27676.4842       nan     nan    nan         nan         nan     \n",
      "transverse=1             10.2966         29633.0997       nan     nan    nan         nan         nan     \n",
      "descending=1             13.9575       1152598.9517 2361.6285  0.0059 0.9953  -4614.7494   4642.6644     \n",
      "sigmoid=1                14.8315       2762005.0734 2361.6284  0.0063 0.9950  -4613.8752   4643.5381     \n",
      "differentiation=1         0.1133             1.1199    0.5132  0.2207 0.8253     -0.8926      1.1191     \n",
      "differentiation=2         0.5613             1.7530    0.9207  0.6097 0.5421     -1.2432      2.3659     \n",
      "mucinous_colloide=1       0.1312             1.1402    0.4365  0.3006 0.7637     -0.7243      0.9867     \n",
      "occlusion=1              -0.2420             0.7851    0.5455 -0.4436 0.6574     -1.3111      0.8272     \n",
      "perforation=1             0.9518             2.5903    0.7728  1.2316 0.2181     -0.5629      2.4664     \n",
      "venous_emboli=1           0.6688             1.9519    0.5078  1.3172 0.1878     -0.3264      1.6640     \n",
      "lymphatic_invasion=1      0.3105             1.3642    0.4903  0.6333 0.5265     -0.6505      1.2716     \n",
      "perineural_invasion=1    -0.5611             0.5706    0.6431 -0.8725 0.3829     -1.8214      0.6993     \n",
      "msi_gen=1                 0.0099             1.0099    0.7635  0.0129 0.9897     -1.4866      1.5063     \n",
      "p53_status=1              1.1099             3.0342    0.7735  1.4349 0.1513     -0.4061      2.6260     \n",
      "kras_status=1            -1.8787             0.1528    0.9891 -1.8994 0.0575     -3.8174      0.0599    .\n",
      "apc_status=1              0.2884             1.3343    0.9495  0.3038 0.7613     -1.5726      2.1495     \n",
      "braf_status=1            -0.6233             0.5362    1.0263 -0.6073 0.5436     -2.6349      1.3883     \n",
      "pi3k_status=1           -15.4090             0.0000  776.3084 -0.0198 0.9842  -1536.9456   1506.1276     \n",
      "Immunoscore=1            -1.6008             0.2017    0.4129 -3.8774 0.0001     -2.4100     -0.7916  ***\n",
      "Immunoscore=2            -1.3813             0.2512    0.6513 -2.1208 0.0339     -2.6580     -0.1047    *\n",
      "---\n",
      "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \n",
      "\n",
      "Concordance = 0.820\n"
     ]
    }
   ],
   "source": [
    "df = X.copy()\n",
    "df[\"Y\"] = Y\n",
    "df[\"delta\"] = delta\n",
    "\n",
    "# Using Cox Proportional Hazards model\n",
    "cph = CoxPHFitter()\n",
    "cph.fit(df, duration_col='Y', event_col='delta')\n",
    "cph.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result in prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C-index on test set: 0.62\n"
     ]
    }
   ],
   "source": [
    "hat_beta = cph.summary.coef\n",
    "marker_test = X_test.dot(hat_beta)\n",
    "c_index = concordance_index(np.ravel(Y_test), marker_test, np.ravel(delta_test))\n",
    "c_index = max(c_index, 1 - c_index)\n",
    "print(\"C-index on test set: %.2f\" % c_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### donc simplement en prenant en compte toutes les variables cliniques que l'on a encodé, on améliore le C-index de prédiction en passant de 0.59 à 0.62.\n",
    "\n",
    "Reste maintenant à travailler sur les images avec les différentes méthodes d'extraction: c'est là dessus que le potentiel est énorme pour améliorer les performances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elastic-net Cox "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le code pour entraîner un cox avec une pénalité Elastic-Net, qui te sera très utile puisque la dimension du problème va significativement augmenter.\n",
    "\n",
    "On fait appel au package glmnet (https://web.stanford.edu/~hastie/Papers/v39i05.pdf) de R depuis ce notebook python grâce au package rpy2 qu'il te faut installer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO: \n",
    "mettre Y en jours (ce sera des entiers c'est mieux) et remplacer les Y=0 par Y=1 jour car Y=0 pose pb ds la plapart des modèles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# je remplace ici les Y=0 par Y=1 jours\n",
    "Y[Y==0] = 1/30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%R -i X,Y,delta\n",
    "%R library(survival)\n",
    "%R library(ggfortify)\n",
    "%R library(glmnet)\n",
    "%R Y = as.numeric(data.matrix(Y))\n",
    "%R delta = as.numeric(data.matrix(delta))\n",
    "%R cv = cv.glmnet(as.matrix(X),\\\n",
    "                  Surv(time=Y, event=delta, type=\"right\"),\\\n",
    "                  family=\"cox\",\\\n",
    "                  alpha=.9)\n",
    "#%R gammaChosen = cv$lambda.1se\n",
    "%R gammaChosen = cv$lambda.min\n",
    "%R fit <- glmnet(as.matrix(X),\\\n",
    "                 Surv(time=Y,event=delta,type=\"right\"),\\\n",
    "                 lambda=gammaChosen,\\\n",
    "                 family=\"cox\",\\\n",
    "                 alpha=.9)\n",
    "%R coeffs <- as.vector(fit$beta)\n",
    "%R -o coeffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAAE9CAYAAABOYwdYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xf8Y1Wd//HXmzoMSMcCOAxgYRUB\nERQBRdeyuyoIFhQHG239KYrgWhEpgthABBuDDWTYBRRXYUUpOvQ2QxUERaooCKiADJ3P749zMt9M\nJt/kJDc3mTDv5+ORR765uZ97br5pJ/d+zucoIjAzMzMzs+FaYtQ7YGZmZma2OHJH3MzMzMxsBNwR\nNzMzMzMbAXfEzczMzMxGwB1xMzMzM7MRcEfczMzMzGwE3BE3MzMzMxsBd8TNzMzMzEbAHXEzMzMz\nsxFYatQ7MCyrr756TJ8+fdS7YWZmZmZPYXPnzr0nItYoWXex6YhPnz6dOXPmjHo3zMzMzOwpTNKt\npes6NcXMzMzMbATcETczMzMzGwF3xM1GZNYsmD4dllgiXc+aNeo9MjMzs2FabHLEzRYls2bBHnvA\nvHnp9q23ptsAM2aMbr/MzMxseHxE3GwE9t13ohPeMG9eWm5mZmaLB3fEzUbgttt6W25mZmZPPe6I\nm43AtGm9LX+qcF68mZnZBHfEzUbgkENg6tQFl02dmpY/VTXy4m+9FSIm8uLdGTczs8VV7R1xSUtI\n2lvS9ZIelnS7pMMkLd/HtqZKuklSSPpGHftrNgwzZsDMmbDssun2Ouuk20/lgZrOizczM1vQMKqm\nfA34CPBT4DDgX/LtF0t6bUQ82cO2DgKKpgw1W9TNmAHHHJP+nj17pLsyFM6LNzMzW1CtR8QlvRD4\nMHBKRLwlIo6JiH2AfYBXA+/sYVubAh8F9q9lZ82sVotrXryZmdlk6k5N2QkQcETL8mOAecDOJRuR\ntGSO+SVwyiB30MyGY3HMizczM+uk7o745sCTwKXNCyPiYeDKfH+JvYENgD0HundmNjSLY168mZlZ\nJ3XniK8J3BMRj7S57w5gS0nLRMSjk21A0rrAgcBBEXGLpOmljUvaA9gDYJrPf5uN3OKWF29mZtZJ\n3UfEpwLtOuEADzet08l3gJuAw3ttPCJmRsRmEbHZGmt4jKeZmZmZLTrqPiI+D3j6JPdNaVqnLUk7\nA68DXhkRjw1438zMzMzMRqbuI+J/BlaXtGyb+9Yipa20TUvJMYcDvwDulPQcSc8B1smrrJSXrVzH\njpuZmZmZ1anujvhluY2XNi+UNAXYBJjTIXY5Us3wNwJ/aLrMzvfvnG/vNtA9NjMzMzMbgrpTU04E\nPkOq/31e0/LdSbnh8ye3lrQ+sHREXJ8XPQi8vc021wC+RSpl+D3g6sHvtpmZmZlZvWrtiEfENZK+\nCewp6RRSmkljZs1zgBOaVj+blHaiHPsY8OPWbTZVTfljRCx0v5mZmZnZOBjGFPcfBW4hlRF8I3AP\ncBTwuR6ntzczMzMze8qovSMeEU8Ah+VLp/WmF27vFvJRczMzMzOzcVX3YE0zMzMzM2vDHXEzMzMz\nsxFwR9zMzMzMbATcETczMzMzGwF3xM3MzMzMRsAdcTMzMzOzEXBH3MzMzMxsBNwRNzMzMzMbAXfE\nzczMzMxGwB1xMzMzM7MRcEfczMzMzGwE3BE3MzMzMxsBd8TNzMzMzEbAHXEzMzMzsxFwR9zMzMzM\nbARq74hLWkLS3pKul/SwpNslHSZp+YLY50k6SNLFku6W9ICkKyXtWxJvZmZmZraoGsYR8a8BhwPX\nAR8GTgY+ApwqqVv7uwB7A38EDgI+DtwAHAxcKGm5unbazMzMzKxOS9W5cUkvJHW+T4mItzYtvxk4\nEngncEKHTfwYODQi7mta9h1JfwD2BXYFvjHwHTczMzMzq1ndR8R3AgQc0bL8GGAesHOn4IiY09IJ\nbzgxX29YeQ/NxtCsWTB9OiyxRLqeNWvUe2RmZma9qvWIOLA58CRwafPCiHhY0pX5/n6sna/vqrBv\nZmNp1izYYw+YNy/dvvXWdBtgxozR7ZeZmZn1pu4j4msC90TEI23uuwNYXdIyvWxQ0pLAfsDjdE5r\nQdIekuZImnP33Xf30ozZImvffSc64Q3z5qXlZmZmNj7q7ohPBdp1wgEeblqnF0cALwc+FxE3dFox\nImZGxGYRsdkaa6zRYzM2TE61KHfbbb0tNzMzs0VT3R3xecCyk9w3pWmdIpI+D+wJzIyIQyvumy0i\nGqkWt94KEROpFu6MtzdtWm/LzczMbNFUd0f8z6T0k3ad8bVIaSuPlmxI0gHAZ4EfAB8Y2B7ayDnV\nojeHHAJTW84jTZ2alpuZmdn4qLsjfllu46XNCyVNATYB5pRsJHfC9weOBXaLiBjsbtooOdWiNzNm\nwMyZsGz+ebvOOum2B2qamZmNl7o74icCAXy0ZfnupNzw+ckHktaXtEHrBiR9jtQJ/xGwS0Q8Wd/u\n2ig41aJ3M2bAFlvANtvALbe4E25mZjaOau2IR8Q1wDeBt0g6RdJukg4jzbR5DgtWPTkb+F1zvKQP\nAQcCtwFnAe+StHPT5XV17r8Nh1MtzMz644HuZuOt7jrikI6G3wLsAbwRuAc4ilT1pNvR7Uad8Wmk\ntJRW5wBnDmY3bVQaR3N33RUeeSSlWhxyiI/ympl14jkFzMZf7R3xiHgCOCxfOq03vc2y9wHvq2O/\nbNEyYwYcc0z6e/bske6KmdlY6DTQ3R1xs/FQd464mZmZ1cAD3c3GnzviZmZmY8gD3c3GnzviZmZm\nY8gD3c3GnzviZmZmY8hzCpiNv2FUTTEzM7MaeKC72XjzEXEzMzMzsxFwR9zMzMzMbATcETczMzMz\nGwF3xM3M7CnB072b2bjxYE0zMxt7nu7dzMaRj4ibmdnY6zTdu5nZosodcTMzG3ue7t3MxpE74mZm\nNvY83buZjSN3xM3MbOx5unczG0fuiJuZ2djzdO9mNo6G0hGXtISkvSVdL+lhSbdLOkzS8sOINzOz\np74ZM2CLLWCbbeCWW9wJN7NF37COiH8NOBy4DvgwcDLwEeBUSSX7UDXezMzMzGyRUnsnVtILSZ3n\nUyLiLRFxTETsA+wDvBp4Z53x1psqE2KMcjKNUe33uE4gMo6PuWq74/gaGeXzNK773a9x/X9VMa6P\neRxjBxHfr3H9nlpsREStF+BgIIBXtCyfAjwI/KLO+MblJS95SVhnxx8fMXVqBExcpk5Ny+uMbdhm\nm3QZl/0e5WOuEjuOj7lqu+P4Ghnl8zSu+x2xeL2+GhanxzyOsYOI79eo2l3cAXOisJ+stH59JP0K\neC0wNSIeabnvAuB5EbFGXfENm222WcyZM6efh7DYmD49zUb3qTW+wAZTrp+/fMoU2OJlnWMvvgQe\nfnjh5SWxDVdema432aRs/UG0ParYhn4fc5XYcXzMVdsdx9fIKJ+ncd1vWLxeXw2L02Mex9hBxAPc\neGO6fs5zytZvbff6hzfgi3d/BkiDmW+5pXw71htJcyNis5J1hzHF/ZrAPa2d6OwOYEtJy0TEo4OO\nl7QHsAfANBeT7WqyiS/afXiUrlMS27DCCuXrDqrtUcU29PuYq8SO42Ou2u44vkZG+TyN637D4vX6\nalicHvM4xg4iHuCf/yxft9v2PdHVIqT00Hm/F+CPwG2T3HccKe1k5briGxenpnS3zjoLnr5qXNZZ\np97Yqka136N8zFWM42Ou2u44vkZG+TyN6373a1z/X1WM62Mex9hBxEf0l4I0rt9T444eUlOKVqpy\nAa4B7prkvpNyR3qZuuIbF3fEuxt1nuO47fe45t6N42N2Putwn6dx3e9+jev/q4pxfczjGDuI+IjR\njAWw/ixqHfFfAU8Ay7a57wLg7jrjGxd3xMscf3z6pSyl617erFViqxrVfo/yMVcxjo+5arvj+BoZ\n5fM0rvvdr3H9f1Uxro95HGMHEV9lgP44fk+Ns1464sMYrHkwsC/wyog4r2n5FOBe4NyI+I+64hs8\nWNPMzMzG1atela5nzx7lXliJXgZrDmMynBNJ6SMfbVm+OzAVmF/RUtL6kjboN97MzMzMbFzUXjUl\nIq6R9E1gT0mnAL8A/oU0M+Y5wAlNq58NrAOoz3gzMzMzs7EwjPKFkI5m30IqJfhG4B7gKOBzEfHk\nEOLNzMzMzBYpQ+mIR8QTwGH50mm96VXizczMzMzGxTByxM3MzMzMrIU74mZmZmZmI+COuJmZmZnZ\nCLgjbmZmZmY2Au6Im5mZmZmNgDviZmZmZmYj4I64mZmZmdkIuCNuZmZmZjYC7oibmZmZmY2AO+Jm\nZmZmZiPgjriZmZmZ2Qi4I25mZmZmNgLuiJuZmZmZjYA74mZmZmZmI+COuJmZmZnZCLgjbmZmZmY2\nAkPpiEt6j6QrJD0k6S5J35W0RmHsWpI+LekcSX+R9KCkayV9RdJqde+7mZmZmVkdau+IS9obOBa4\nD9gLOBp4JzBb0vIFm9gWOAC4F/gK8FHgwnx9paRn1rDbZmZmZma1WqrOjUtaHTgYuAx4TUQ8kZdf\nBvyc1DH/QpfNnAesExF3Ni07RtIlwDHAf+WLmZmZmdnYqPuI+PbAVOCoRiccICJOBW4Cdu62gYi4\ntqUT3nBivt5wEDtqZmZmZjZMdXfEN8/XF7W572JgA0kr9LnttfP1XX3Gm5mZmZmNTN0d8TXz9R1t\n7rsDUNM6vTowXx872QqS9pA0R9Kcu+++u89mzMzMzMwGryhHXNLKpMGRpY6MiL+R0lIAHmmzzsP5\nemqb+7rtz8eAtwMzI+LXk60XETOBmQCbbbZZ9NqOmZmZmVldSgdrrgzs38N2jwf+BszLt5cFHmpZ\nZ0q+nkcPJO1Gqp7yf8CevcSamZmZmS0qijriEXELKY2kV3/O12sBN7bctxYQTet0JWkX0hHuM4C3\nRsRjfeyTmZmZmdnI1Z0jflm+fnmb+7YAboiIf5ZsKHfCvwucBWwfEe3SXczMzMzMxkLdHfGfkVJS\n9pS0ZGOhpG2B9YBZzStLmiZpA0lLtyx/H6lm+K+BN0fEw5iZmZmZjbFaJ/SJiLsl7Qd8FThL0n+T\nUlI+BlwPHNESchywDbAucAuApO2A7wH3k2qHv1VaIEvmnxHxvzU+DDMzMzOzgau1Iw4QEYdJuhfY\nGziS1KE+CfhUYVrKpqQj9yuTK6C0uBVwR9zMzMzMxkrtHXGAiPgh8MOC9V7VZtkBwAED3iUzMzMz\ns5GqO0fczMzMzMzacEfczMzMzGwE3BE3MzMzMxsBd8TNzMzMzEbAHXEzMzMzsxFwR9zMzMzMbATc\nETczMzMzGwF3xM3MzMzMRsAdcTMzMzOzEXBH3MzMzMxsBNwRNzMzMzMbAXfEzczMzMxGwB1xMzMz\nM7MRcEfczMzMzGwE3BE3MzMzMxuBoXTEJb1H0hWSHpJ0l6TvSlqjz20tIekiSSHptEHvq5mZmZnZ\nMNTeEZe0N3AscB+wF3A08E5gtqTl+9jkB4ENB7eHZmZmZmbDV2tHXNLqwMHAZcBrImJmRHwO2Al4\nAalj3sv21ga+AOw/6H01MzMzMxumuo+Ibw9MBY6KiCcaCyPiVOAmYOcet/fNHPf1ge2hmZmZmdkI\n1N0R3zxfX9TmvouBDSStULIhSW8DtgU+0NypNzMzMzMbR3V3xNfM13e0ue8OQE3rTErSSsCRwNER\ncXFp45L2kDRH0py77767NMzMzMzMrHZLlawkaWXgoz1s98iI+BspLQXgkTbrPJyvp7a5r9WXST8a\nPt3DPhARM4GZAJtttln0EmtmZmZmVqeijjiwMr0NkDwe+BswL99eFnioZZ0p+XoeHUh6BbA78O6I\n+EcP+2BmZmZmtsgq6ohHxC2kNJJe/TlfrwXc2HLfWkA0rTOZbwBXAZdIek7LfVPzsn9ExD197J+Z\nmZmZ2UjUnSN+Wb5+eZv7tgBuiIh/dtnGOsAmwB9aLgCvzn8fUHlPzczMzMyGqDQ1pV8/Iw2y3FPS\nCY1qJ5K2BdYD9mteWdI0Us74HyPisbz4PcAybbZ9MjAX+CILH203MzMzM1uk1doRj4i7Je0HfBU4\nS9J/k1JSPgZcDxzREnIcsA2wLnBL3sbP221bEsCdEfHjWnbezMzMzKxGdR8RJyIOk3QvsDfp6Pj9\nwEnApwrSUszMzMzMnpJq74gDRMQPgR8WrPeqHrbZz+BRMzMzM7NFQt2DNc3MzMzMrA13xM3MzMzM\nRsAdcTMzMzOzEXBH3MzMzMxsBNwRNzMzMzMbAXfEzczMzMxGwB1xMzMzM7MRcEfczMzMzGwE3BE3\nMzMzMxsBd8TNzMzMzEbAHXEzMzMzsxFwR9zMzMxsETZrFlx8MZxzDkyfnm7bU4M74mZmZmaLqFmz\nYI894JFH0u1bb0233Rl/anBH3MzMzGwRte++MG/egsvmzUvLbfy5I25mZma2iLrttt6W23gZSkdc\n0nskXSHpIUl3SfqupDV63MZUSZ+TdG3ezt8kXSRph7r228zMzGyUpk3rbbmNl9o74pL2Bo4F7gP2\nAo4G3gnMlrR84TZWAc4HPg6cDewJHABcCawz+L02MzMzG71DDoGpUxdcNnVqWm7jb6k6Ny5pdeBg\n4DLgNRHxRF5+GfBzUsf8CwWbOhJ4DvCyiLiupt01MzMzW6TMmJGu9903paNMm5Y64Y3lNt7qPiK+\nPTAVOKrRCQeIiFOBm4Cdu21A0nTgXcAxEXGdpCUlrVDP7pqZmZktWmbMgFtugSefTNfuhD911N0R\n3zxfX9TmvouBDQo61f9O2s/rJP0ImAc8IOlPOe3FzMzMzGzs1N0RXzNf39HmvjsANa0zmefn60OB\nlwAfIB1Jvxk4XNKBkwVK2kPSHElz7r777p523MzMzMysTkU54pJWBj7aw3aPjIi/kdJSAB5ps87D\n+Xpqm/uaPS1fLwO8IiLuzft0EnAd8AlJR0TE31sDI2ImMBNgs802ix7238zMzMysVqWDNVcG9u9h\nu8cDfyOlkQAsCzzUss6UfN1Spn4hjbjTGp1wgIh4TNIJwOeALYDTO21k7ty590i6tWDfLVkduGfM\nYkfZ9jjGjrJtP+bxiB1l2+MYO8q2/ZjHI3aUbVfdbytXXNGvqCMeEbeQ0kh69ed8vRZwY8t9awHR\ntM5k/pSv72xz31/y9SrddiQieqpbvriTNCciNhun2FG2PY6xo2zbj3k8YkfZ9jjGjrJtP+bxiB1l\n21X32+pRd474Zfn65W3u2wK4ISL+2WUbl+brtdvc11j21z72zczMzMxsZOruiP+MlFqyp6QlGwsl\nbQusB8xqXlnSNEkbSFq6afG5wK3AtpLWalp3eeA9wD9oX5XFzMzMzGyRVWtHPCLuBvYDXgqclauY\nHAj8N3A9cERLyHHA70hpK41tPAF8kJRnfpGkT0nah1T+8NnAxyLiwTofx2Jq5hjGjrLtcYwdZdt+\nzOMRO8q2xzF2lG37MY9H7CjbrrrfVgNF1F9MRNL7gL1JpQjvB04DPhURf21ZbzawDbBuzktvvm8b\n0oDRzYElgSuAL+bJgczMzMzMxspQOuJmZmZmZragunPEzczMzMysDXfEzczMzMxGwB1xW4ikbrOd\nmpmZmVlF7ojbfJK2lHQdqaINkjaW9K3CWEnaWdLn8u1pkl7aQ9tbS3p//nsNSev28RAa23pdhdgV\n+o2t0Oaqw27zqWIUz1dVkjboM67SY228N/uM7Wufc+zYPUdQ7f+V46t8DlVqe1Qan+HDjq2iaruS\n+q5EUmespCUl/aekz0vaquW+z/bbrg2eO+LW7GvAvwH3AkTEVcArC2O/RZq4aad8+wHgmyWBkvYH\nPgl8Oi9aGji+sN12vlch9rpuK0h6kaSLJd0uaaakVZruu7RL7FaSfifpWkkvk3QmcFneVruJr4pJ\nuqbL/c+W9D+SzpP0meZ6/ZL+d9zazTo+X1Weqw7bPL2fuCZn9BnX9bXZxW4VYvvdZ6j5PVWw7X6f\nryr/L6j2OVSp7W7vybpigQNHEVvxPdm1XUmrTnJZDXjDohgLHE2qQncvcKSkw5vue0u3x2zDUzTF\nvS0+IuJ2Sc2LnigMfVlEbCrpirydv0tapjB2B+DFwOU59s+SntYpQNLPJ7sLWK1L7D4dYkuO3n0b\nOIBUy3434HxJ20XEH0k/Ijr5GrBjbuf/gO0j4nxJmwJHAVt1CpY02QeogGd2afv7wE/yfu8KnCNp\n24i4F1hnUWw3t13l+errucrPx2RtbtKlTSQd2SF+5Q5xlV6bku7vEL9cl9i+9jnHjvI91ffzVeX/\nleOrfA5Vbbvv92TF2Ks7xD6jxti+35NV2s3uJk0s2PzlGPn20xfR2JdGxEYAkr4BfEvSKaSDZeoY\naUPljrg1u13SlkDko5Z7kSZYKvGY0uypASm9BHiyMPbRiAhJjdjlC2JeAewM/LNluUgTSHXyBeAr\nwONt7is5S/S0iPhl/vurkuYCv5T0bvLj72DpiLgGQNLdEXE+QERcLqnrly9wImlG2nbtTOkSu0ZE\nfCf//WFJOwPnStquYL9H1S5Ue776fa4uA86h/RdWx05p9n7gY8Ajbe7bqc2yhqqvzX8Am0fEXa13\nSLq9S2y/+wyjfU9B/89Xlf8XVPscqtp2lfdkldhnkM6c/r1luYALa4yt8p6s0i7ATcBrIuK21jsK\nnqtRxc4/EBYRjwN7KKU8/ZqyH8c2LBHhiy9EBMDqpA/nu4C/ktJDViuMnQH8HPgTcAhwA/D2wtj/\nIp1GuwnYHbgI+HCXmNOBV09y37ldYi8EXjLJfbcX7O9VwEotyzYC/gDc2y226e/tW+77bUHbc4EN\n+9l34FpgSsuy1wI3An9ZFNut+nz1+1wBvwWeW+E18mtgy0nuu7nG1+bBpCNh7e77Uh37PKD97vs9\nVeX5qvL/yutU+Ryq2naV92SV2O8BW09y3wk1xvb9nqzSbl7nQ8DGk9zX7btqVLHHA//eZvluwGPd\nHrMvw7t4Qh8bGKUBXa8hHWU4OyJKj6Y3Bja9Psf+KiLOrGcvQdLzSV/u97S57xnR5ghVyzrvAm6K\niItblk8D9ouI3TvEbgecFRHzWpavD7w1Ir7cpe1XALdG+yMkm0XEnA6xewOXR8Q5LctfDHw5IiYd\nXDaqdvN6fT9f/T5Xkt4GXBMRN7S5b/uI6JjbrjQA9+HW57mbqq/NKvrd5xw7svdUXq/S8zWOKr4n\n+44dlcXxObbFgzviNt8kOaL3AXMi4meTxHSs+BERfxvEvplZbyStGBH3N65HvT+LulH+v/xcjQ9J\nG0TE9Y3rp3qs1c9VU6zZFNKglz/ky0bA2sCuko6YJGYuMCdf3w38PsfenZdNStIDku6f7FK605JO\nar7uhaRPNF8PM34Abe/cfP1Uj81xQ/9/S/rX5us+2j2i+bqHuEqvD2B2y3Uvbfe1zzlmZO+pHNfv\n8zW75brXdvv+HBpA26N6P5/ffD3E2L7fk1XazU5ouX6qx1rN3BG3ZhuR8h2PioijSHm8G5Cqmry+\nXUBErBsR6wFnAdtGxOoRsRrwJrqUPIuIp0XEisDXgU8Ba5E6/p8EeukAPCdfP7eHmIZ3tlwPM75q\n2/u0XD/VY2E0/++vtlz3qlECdJse46q+Phr6qZDQ7z7DaN9TUP356reiRJXPoaptj+o92Zj8rWSA\n/SBjqzzHVdptVqXyyDjGWk3cEbdmq7DgaOrlgVUj4gnaV1JotkVE/KJxIyJOB7YsbHe7iPhWRDwQ\nEfdHxLeBN/ey4wNQ9QNqlB+O4/ihPo6PeVRfYuP65TnK53gQ8eNoceukjeM+my3A5Qut2ZeBKyXN\nJn3AvRL4glI5wbO6xP5ZabauxkQ8M4A/F7b7oKQZwP+QymntBDzY476bmZmZjRUfEbf5IuJ7pAll\nrgdOAT4L/D4iHoyIj3cJ3wlYA/hpvjyd7rWHG95FmuSmUTbx7XmZmVXnEfm9GeX/y8/V+KjyXI1j\nrNXER8RtPkm7kSbxWRu4EtiCVNO764CYXB1lL6UZMSMiWie46BR7C9VSUXx60mxharm2zqr+vwaR\nFjJuz9U4psIMKuWpn+2MY6zVzEfErdlewOak+rKvJk07/4+SQEkvUpre/rfAtZLmStqwMHZtST+V\n9Nd8+YmktXvY76+0XPdidr7+TR+xVeOrtv37fL1QXd2naCyM5v/d+FH5QB9twkSlglk9xs3O1/2+\nPt7Rct2LfvcZRvuegv6fryr/L6j2OVS17VG9J/duuR5WbJX3ZJV2Ic2m2nz9VI+1mrmOuM0n6bKI\n2FzSlcDLIuIRSddGxAsLYi8E9o2I3+TbrwK+EBFdB2xKOpP0xf+jvGhnYEbBRC9LAh+JiK91a8PM\nzBY9kpYAVhhB7fae283jpR6KiCclPY9UVez0iHisztgcvw5pZtGzJC0HLBURRT9EqsRa/dwRt/kk\n/RR4P/BRUjrK34GlI+INBbFXRcTG3ZZNEntlRGzSbdkksZdGxEu7rVdK0uuqzOpZJb6XWEkiDYhd\nLyIOyjMQPjMiLq0jVtKmnbYZEZfXEduynaWB/8dEab1zgO8Ufgn2Fas04+mf8o/SV5FKfB4XEaVn\nitqVhLsPmBsRV3aImwp8DJgWEbtLei7w/Ig4rbDda1g4H/Q+Us3/gyPi3sJ9nS8iDi9se0tgOk2p\njxFxXEnsgOK3JnU6fiBpDVJn6+YuMVsBBwDr5HaVmo31CttcA9i9zX7v0iHmpIjYsc1z1Wh7o5K2\n87bWatr3Rtvn1hmbO5P/xcKPuWsqo6QTgA8ATwCXASsCX4+IrmcTqrwnq7Sb4+eSjiqvAlyQt/Fo\nRMyoOXZ3YA9SFbP18+fBdyLiNXXG2nC4I25tSdoGWAn4ZUQ8WrD+T4HLWfCo9ksiYoeC2LOBHwD/\nnRftBLy/8EPma8DSwIk0VVop7dy12d5tETGtn9iq8b3ESvo28CTwrxHxL5JWAc6IiM3riJX0JCnt\nqDGFeXOuYXT68q0S27Kd75Ke62PzoncDT0TEbnXF5rNDm5E6G78Afga8sOTHaY4/Icefmhe9Cbg6\nb+/kiPjyJHEnkibEek9EbJg75heW/DjN8V8mdTYaKSbvJNVPvhPYOiK2bROzf/7z+aQUtZ/n29sC\nl0ZE1wlfJP0IWJ80xuSJvDjkrGALAAAgAElEQVQi4iOF+101fn/S//v5EfE8SWuS/s9bdYm7npSq\nMLepXdr9YJkk/kLgvDbxP+kQ86yI+Es+WrmQiLi1sO0vkVJarmPB/9l2NcdeBXyHhR9zx4nccuyV\nEbFJrpa1KWkeibklPz6qvCertJvjL4+ITSV9GFguIr7cw0GjKrFXAi8FLomIF+dl10TEi+qMteHw\nYE1rKyLO6TFkF+BAUrUVSF9Kkx4NahN7FPA10pGhC0lH5ks0PsQOaloWdBhgKunnk90FrNatwSrx\nVdtu8rL8oX4FQET8XdIyNcbuA7wNeIhUZvKnPQzIrRLbbPOWMyy/zp2BOmOfjIjHJe0AHBURRzX+\nb4XWBjZtPN7cUfw/0pH5uaSSoe2sHxHvkLQTQETMy2cySr02IprPRFzT1BFo26GOiAPzPp6b9/mB\nfPuAvM8lNgNeEP0f4akavwNpbMvlABHxZ6UB5N3cF2nug35NjYhP9hIQEX/J10Ud7g62J/3w6DbX\nw6BjH48050M/ls5nqbYHvhERj0kqfc6rvCertAvphOLLSWcUd83LlhxC7CMR8WjjI0DSUpRXQKkS\na0PgjrgNRET8HSg6atUm9lag6xGYSWJf3UfYK0hH7Fs7gyIdOagzvmrbDY8p5cgHzD81/mRdsRFx\nBHCEpPVIR1fPlnQraRzApCkWVWNbPCFp/Yj4Y97v9Wg6EldT7GO5M/xe0pFhSEfWSz2dBSfDegx4\nRkQ8JKlT5+dRpVzOxnO0Pt0n1Wq2pKSXRk43krQ5E1/6j3eJfQbQfBbs0bysxG+BZwJ/6WFfBxn/\naEREo3OllJdb4jeSvkI6kDD//9zDmbXTJL0hmiY160bSA7RJSWlcR5p1uMRNpNdkP53pKrGnSvog\nqVxt8//sbwWxRwO3AFcB5+azAqW52lXek1XahZS2+WnSwYRr8+dI6aDiKrHnSPoMsJyk1wEfZOIs\nW52xNgTuiNtAVMwX7Dm/sin2GcAXgDUj4j8kvQB4eaSa6JO5GJjX7qi/pJLKAVXiq7bdcCS5Xruk\nQ0hHnD9bd2xE3CTpZ8BypPSO55HSCGqNzT5O6jDdROqsrEP5mZN+Y99Pyik9JCJulrQuE+lXJWYB\nl+THDanjcELuIF7XIW5/4JfAsyXNItX3f18P7e4GfF/SCqTHez+wW2730C6xxwGXKqWbQTp6+MPC\ndlcHrpN0KQt2zkp/aFeNP0nS0cDKSrmxuwDHFMS9LF9v1rSs45m1FnsBn5H0KOnHFnTpTEdEyZH6\nEvNIE7GdzYL/s5IDI1Vi35uvm+eYCKBrXn1EHEn6HGq4VVLpQZW+35MV222cKT6n6fZNFB6AqhJL\nSqHZFbgG+E9SSs53hxBrQ+AccRuIivmCPedXNsWeTsov3zciNs6n3a5YHPLfJG0AvIbU0To7In5X\nV2zT0ew3A7eTUkz+LyIeKmir79g221qWlMMMcEMvp9SrxOb4VYBnR8TVPcZtDjSqB10QEXMK41Yj\n1fIXcHFE3NMlpN02VgKIiPt6jNuUiVJn50ZE0al/pbElCylNdasan7fxOuD1pP/br6LC4OthkbQx\nC/6/i19jkt7bbnlEHNtu+aBiq5D0uUnaPajd8kWlXUm/oU1aR+EBp75i89nL46JgUOcgY2143BG3\ngZA0NyJe0mds0YCVSWIbJRevaBqIUjoApmo5qVGWslq1zeIHCtvuOVZpwOXVpIFR99PyhRIdKmpU\niW3ZzpLAG1n4zEnX+H5jJc0mpU0tRfqh+FdSZ7pjhZE2bT+jpd3busRsBVwZEQ8q5XRvSqruUDqA\nb1ngrSz8eCftcEhaMSLun+T1UZpy0DhL1Rj4e2lE/LUkblDx/cg/WPZnwao6B/XyA0bSdk3xs6O8\nws1epDOCjfE1OwAzI+KoHtpehnSGCdKPzKLPkSqxSgOI9yFV9tlDPVT2kfSxpptTSIOYf1d4FvRm\n2ndoux6Jr9Jujm/+jptCeo89HhGfqDn2fNLg+q6FEwYZa8Ph1BQblCr5gj3nVzZ5MB85bOSEbkEq\n01biXOAV+UjnGaRyUu8gDaapO75q25cDzyaVmBSwMnCnpLuA3buciegn9iAmvvxWKNzHQcQ2OxV4\nmHSKtTQfvmrsSrlzuhvpyNL+kno5WvlhUgfvLtLZnkYecLcqDd8GNs5HSvcBvkdKGWl7xLiNn5HL\nJFKe/3sCqWMyl4lc5YailANJO5ImtJmd44+S9PGI+HHJDvQbr4XzrRdQkG/9fVJ++o759rtJZ9re\nUrjfXyT9eGhMgLSXpK0i4tMF4buSBlA/mLf1JdKMxkUdcaUSfseScp9FSmd6b5SVIOw7lvT/mcvE\n2Z47gJOBrh3xiDisZT++CvyqoE1YMH1oCvB2oO2PxwG32+4M7wU5jarWWFIu/wVKg/2bq4OVHMSo\nEmtD4I64DUrf+YJM5Fc+Qsqv7GWw0j6kMmvrS7oAWIOU81xCkapR7Ap8K3I5qcLYqvFV2z4T+HFE\n/ApA0utJR1h+AHyLiZzXgcRGxAE97NvAYlusHT3UVh5Q7FKSnkXqoO3bR/xepKOERWXwmjweESHp\nzcA3I+J7+bVSau2I+PdeGoyIN+XrdXuJa7EvqULNX2H++I+zgKKOeL/xjXxrSZ8nDfT8EelzZAbw\nrIJ214+ItzbdPrDH9+MbgE0i4sm8H8cCV5AG53UjFhw43PjBVuow4PURcUNu+3mkUrAlZyirxFat\n7NNsKqnCUFdt3ktHKNXobpt2Mqh2YaGziUuQ/k8r1R0L/DFflgB6HVtQJdaGwB1xG4gqX95RYdBS\nRFye80qfT/ry6uW0rNR/Oamq8VXb3iIidm/ciIgzJH01Iv4zpyUMNHay3MqJTcTnJ7uzSmyL0yW9\nPiLOKFx/ELEHkY6YXRARlynlu/+hh/jbKT9D0+wBSZ8mVdh5pdIsgL1Ua7lQ0osi4po+2u47zQJY\noiWV5F5SB6BU1fjtYsEyld9WGr/SrZP2kKStI+J8mJ8a1OsYhpWBxhnA0g4WpB/AlygNjhVpLEWn\nweatlm50pAEi4vdKJfrqju27so8WnMRoSdIBlKLPAS04QdgSpCPkRX2ZKu1mzWeLHgduZuLzu7bY\nmCgtukK+XVz+tUqsDYc74laJpI6nbiPilMnuU7XZGidr93mSOrbbpEo5qarxVdv+i6RPkgY+Qkpr\nuSvnI3dLvegn9sE2y5YnfZGsRucvsyqxzS4Gfpo7pb2eOekrNiJOJp1ub9y+iXT2oNRNwGxJ/8eC\nKVvdTgu/A3gXsGtE3Kk0+2nR7H/Z1sD7cj7tI0w83pIJU9qlWWwZEZ8paPeXkn7FxORc7yBVaShV\nNf5Bpcla/ofU6dmJ9q+/Vv8PODbniovUoX5fD+0eClyhNCBPpB8xnyoJjIjDlcYibJ33+f1RODg2\nm6M0YdXx+fYM0gyqdcdWqezzpqa/HwfuiohuZTUbmtNLGh3aHSdZd5DtAvxLRDzcvKDgwEflWEkb\nks7yrJpv30Oa7OvaOmNtODxY0yqR9IOWRY0XVOOLv9MUz80dz3ZTPHealKfR7tNJOYq/zrdfTZqB\n8E1tA58iJK1O+iLcOi+6gDSh0n2kwVM31hGb459GSrnYFTgJOCwKB9RVjL2ZdLTwmujxg6vf2Hyq\n/tuk2t8bStqIdNT14ML4/dstbxylmiRmSeCs6K9GfmMbfc/WqJQD35xmsSSpElHp7INvJXXKAM6L\niJ92Wn+Q8ZKmA1/P8UF6bX80Im4pjF8RICJ6qS3diH0WCw4yvbOH2EaVmidJZ1+KZwbOnbkPMfF+\nPo+U7tb16HSV2BzfV2UfST+KiHd3WzZJ7Hr5B3HzsnUj4uY6283rXh4LTpTVdlkNsReSKoP9Jt9+\nFWkehi07BlaMteFwR9wGQlJjFPh0Js60RBSUhcqnNz/IxBGh84Bvtx49mCT2DOC9kWepy1+GP4yI\nf+sQc0REfFTSqbQffd+xZnGV+Kptj1LOcdyHdNTsWFIVj7/XHdu0jXOBVzU6iMOIlXQOadzD0TFR\nlee3EbFhr/vQY7tnA2+J3ssOVq58kjvir2qsm7c1u7QjPk4k7RwRx0tqWwWn25kLSRtExPWTnd0r\n6VDn1K23Az8hdWi3B04u/bE3KpM85vuAW7sdZW7tgCqVnb06Il5Q0G67Dm1R1a5+25X0TGAt0pmD\ndzGRw78i8J2I2KCO2KZtXNWSctV22aBjbTicmmKD8r/AP0gVORod6NJfeceSyto1Jlp4F6lCRMnp\nxmc3OuHZXcC0LjGNyR++Wrh/g4yv2jYwfxDbJ4AXkioHAD1NoNRTrNKsg28BZgIv6iXPsEpsi0aa\nx+n0luZRJXZqRFyqBcegdT2VPYAfXP8kTUt/JgtWOug2AUjlyif0kWYh6fyI2FqTzBbZLQVoAPGf\niDTg+Sja/78n+781Zt7sd5zKPsAeLJguMb9ZyiYEmgFs3DjwkFODrgQ6dsQlnRQRO2rBvOeJxjv8\ncKoS2+RbpLKaV5Oepw2Ba4GVJP2/aDMeQ2ncQ2OWx8ZZB5Fmb53ZqTGluQ9emLffnJq4Ik2fY5PE\n9t1u9m+ktJu1gebPjAfyduuKbbhJ0n5MfH/sTPpMqzvWhsBHxG0gqhwllHRd6xGJdssmif0G8FwW\nzCm9MSI+3M++jIt8JuBE0mymHyBVrbk7Ij5ZR6xSLfBHSJ3QnjpKVWJbttNzmkfV2Nxx35N0hHJT\nSW8j5W3/R5e4l0TEXPU5QY1GNNFKU/t9p1mMgqRtI+LUUf3fJE1pPYPXbtkksb8BdoiIf+TbKwOn\ndPtRLelZEfGXftKQqsQ2beMUYL/IucZKsxofRPqRf0p0mMtB0qFRVtqxOebNpLMF25EqZTU8APxP\nRFxYsI2e222Jf2sUTDRXQ+wqpPTB5rPGB5acVawSa8PhjrgNhKSZwFHRR5UGSccD34iIi/PtlwEf\nioj3FMa/hQVnpeuYUyrpctLkGf8dEX/sY3/7jq/adtN25kbESyRd3Th6pTy5UR2xkpaOHiYJGVTs\nJNubGhHzhhGrNIh2Jmkcwt9JA8N2jsKc4ypyyta0aKpqURDT9wDolu1sxMKTAXUdAK3qObiV4pti\neqoQIenLpCPQD5EGIG4E7B0Rx3cMnIjvOf+36ej9NNKPnjPz7deRfvyU1jAfycRi7Q6+NJapYFK1\n3EF8LguelSupff7yiLio23qDbrcp/o0sfDaxdGbOvmPtqcupKTYofVdpINVTvVBSY7bBacANjdOm\n3baROwglVVIaViGVGvuNpDtJR9NPjIg/DyG+atsNjS/Kv+QP9z9TOKlFn7GXkE5D96NK7HxK5R6/\nR5oUaJrSZDf/GREfrCs20qCw1+YOyxIR8UDhvrY95d+03Y6vaUnbktKXlgHWlbQJaabHbiktjRSJ\nKaSybleR3osbkaphvLxg37+f17+WiSo6Qdl77IUt21qKsprUA4nXghUiJOluyipEvD4iPiFpB9Lk\nNm8hTbrVsSOuifzf5SS9mAXzf6d2abNRnWQuaSK0htld4lqNamKxayV9mwWrL12nNAC0Y0deaYKs\nvUjpGleSBnxeRFkqzxWSPsTCHdqSWTmrtIuk75Ce11cD3yXNWVE0KU/F2DOBtzedNVmFdBZg0rFQ\ng4i1IYkIX3ypfAHWaXepEluyDdIX5h9Ig4TuJ52mvL9LzOVNf7+ClOt4J6l84B4F+9t3fNW2m2Lf\nRKpVvGGOnQtsW1csqWpGv6+NvmNbtnMJaUbQK5qW/bbOWGBZ0piFz5BqUX8O+FwPr+kv58uL8uVL\nwBcL4ufm56jnx5rXPYWUj9+4vSFpEqeS2Ov6eG4+nd97j+f3YeO9eC9waN3xTdu5EHh10+1Xkaoo\ndYv7bb7+LvDv+e+rCuLem99DD+TrxuXnpMG2lV/3Bftweb7+MPCJ/PeVQ4hdDvgY6UfET0mpblNJ\ntb1X6BJ7DakTfWW+vQEpnaWk3ZNJJU//mP//Z5AGf5fE9t1uXv/qlusVSJV96o5d6DO03bJBx/oy\nnIuPiNtAREFOYR2xpE7OthHxuz7bPg84T2kq8teRjuqUDN6pHF+x7b9HqqhxH+kIC0qTkNQVu4Ym\nqSwBXQc9VoltXfd2LThw8onJ1h1QbD9Txc9/TUt6XeRqK9knc3pStxrTj0XEfS3720vFl+dHU5pY\nRPxW0r8Uxl4k6QURcV1pYxFxKHBovzm4VeObLB+5TFve7ux8NqOb0yRdT0pN+X9KA5q75ndHyj0/\ntmL+75tIHct1SGepexo7kTYx/InFIuIh0hmYdgNV/ynpJ7HgbKXNHo6IhyUhadlIlWeeX7jPz4mI\nt0t6c0QcK+kEUt5ziSrtwsQkT/MkrUn6oVgyc2vV2CclTYuI2wCUcvtL84qrxNoQuCNu4+6uPjrh\nv29dEBFPkHJDf1lzfNW2G45i4XSPdssGFbsk6QhOP1NYV4ltdrukLYFQmv1vL6D0ue83tuep4ltI\n0lYRcUG+sSVlM0VeK+ldwJKSngt8hHS0t9TVWniilqsLY48jdcbvpMc0s4j4dJUc3Krx9FkhIiI+\nlfPE74uIJyQ9SKo7XyQiflIh//cI0pm9nuvjZ3vR/+RgVWK76VSh509Kg1L/FzhT0t+B0gMyjbSX\nf+RUpDtJ80mUqNIupB9sK5Mm17qc1KH97hBi9wXOVyqnKtLZ1D2GEGtD4MGaNtYkfR14JumDtbks\nXS8545Nt+71RodpClfjJYvPRqy1JM3N+remuFUmVFyatDVsxtmjiiUHHtmxnddJkLa8lfaGcAewV\nEffWFasKg5Bz/EuA75PSTEQa8LlLdBk0KWkq6Qv09XnRr4CDo6AKR46fQpotsjFN/bmU1+a/kVSW\n7xqajsKXnLmaLAc3CspqDii+uUIEpCOlB0SXChGS2g4Mj4jjCtttm/8bEV2nMFeqmvKa6KM+/qKs\n9H2vVFloJcoHie5Gqrn+IuCHpB/5+0XE0T3uX0/t5phlI092lHPhp5COshdNntRvbI5ZnfR+gB4m\nT6oaa/VzR9zGmhae2RO6zOjZw7YrdSDr6LzmL49XkcoOfqfprgeAUyPiDx22WSX2ipYUi8nWW6W1\n01MldtQkXQc8h1QtpddByM3bWYkUWDRBj6RNu3XW6yLpoojoOqhzkthrSBVALo6ITZRqP38hyiuA\nVIrvl1IFk4YpwGtI+dNvK4y/OiI2arpegdTBe0VB7Oak1JRz6L0+PhrynAKlOn3+qdrMmgvNotlu\n2aDbzeuOambNrUh57Q9K2pl09vLrhT+O+4614XBqio21iHh/jZuvmkpRJb5tbKT60+dI+mFTHnJj\ncFTHabmrxJI6JiXOZuEUlyqx8+VOw+4sXFavpFpCv7Ed64UXtLssTTPOKud8F6QsHKZUkePHpKo6\nv+2x3a2AA5jIO260WzKhzxU57/ZUej/LVDUHt694TTJxUtO+d6w2Ey3zDuQUgv+ZZPV2Gmca+sn/\nPYQ0gdMUUpWcXs0izQvwJprmBRhCbDedPv9aq+MsSXl1nJ+w8OfEjwvj+2pXFarjVIlt8m1gY6Vq\nT/uQKkAdB7Sdp2CAsTYE7ojbWFOqfftt4BmR6tduBGwXg5keuurpoirx3WIPlfQB0oDDy4AVJX09\nIr5SsO2eY6NgavRsoS/fKrEtfkZKNTiLHgZp9hOriSnii8oVdmm3n8Ger85f4DsCR0takdQhL31d\nfw/YO7fb6/9qubyvr29aVlq+sGoObr/xlWaqbeNBYN0e1j+1Tf7vMYWxa0afk6Flq0XE9yTt1fRj\n+7IhxHaz0ARhemrMrHkYE59Vvc6s2Wtsw+MREUqTGn0zP2dd054GEGvDEItA6RZffOn3Qjqd+1L6\nLPPWZduVSjxVie8Wy0T5rRmkD/elyWWxCrbdd2zBti+vK5bCsmqDiCWlotyUr58A7iEd5XwCuLmH\n7VR+LZJyYX8EPNpDzCWDeD4r7vc2pFkQlxlFfI9tnUoqOfhz4LT83HctM5ljlwC2bLq9LLBSD21/\nmVTHvN99vzhf/wp4I/Bi4I91xQIn5etrSAOAG5drevgMKi5J2RTzZuAH+X34g6bLkc3//0G32xL/\n1hHFnkMaVPsH0pioJUiDe2uN9WU4Fx8Rt3E3NSIu1YJl3h4f0LYvGGF8t9illap/bE+alfQxSaVH\n4KvEjtJpkt4QEb+oOzYi1gWQdAyposQv8u3/IP3fSl0o6UXR42BPpVKD7yCltdxLSh/4WA+b+I2k\nr5COYjenl3TNO5d0ZJvF9wFzIuJnk8SsGBH3N51JgNQxgzSYruNZkQHETzaBUmlOf/MR9ceBWyPi\nT11iIG38SUnfJHViiTT4rvjsB2lQ7X9JepR0hLbX8oUH5zEIHyNVP1qRdDakrti98vWbCtto5zRJ\ny0cPecv5tfczVZtZs+d2W6ydz049QDrjsSnwqYg4o+bYd5DmM9glIu6UNI109qVElVgbAg/WtLEm\n6XRgT+DkiNhU0tuAXSOip9xeSb+OCgOUqsT3EyvpI6RTv1eRjmRNA46PssFhfccWbLtoYGYvsZIe\nIHWyBCxP6uQ8RkGHpUpsjr8mIl7UbVmH+L4Ge0q6iJSjfHL0PutqoxJHqyh5nSlVitmANHEKpB8D\nNwOrATdFxEfbxJwWEW9Smlm38f+efx1dctMHEL9Op/t76Gj1RdJXSTM0nhJD/lKVtEZE9JXXXTF2\neeCh/EPkeaTXTGnlk6uBjUkzuP6QVGlmx4jomresVGbyYFJd7l/mbewdER1nQa3abo6/KiI2lvRv\npJz6zwI/irIBl33H5vhnkAYyQ6rI89eSuKqxVj93xG2sKdW9nUkqy/d3UodhRqcv3vxhvMAi4HnA\nDVA0/Xjf8VXb7rBdAUtGxOP59nujsHRiSaxSZYfVI+L0luVvINVyn5tvrxotOeFVYkdN0q9IueXN\n9bhfGYXTQ0/WQeylY6hUku/ZEVFaB7wSSRcDW0Wqb4/SNPPnkUoCXhMRLxjGfvSrl05H0w+1he6i\nh6PSeTvLk1KXHuolPr//ZgDrRsTnJT0beFZElE5//nvgFtJZk1Oih8pDFWPnkmpSr0I6g3cZKX1q\nRkHs5fnAyeeAOyLlLZdWELkyUkWdHUhH5fcBzo0O5VcH0W6Ob1TF+TowOyJ+WnrwoWLsjqSj2LNh\nfi3wj0fEj+uMtSHplLfiiy+L+oX0IbwPqebyfvnvXYFNOsT8nNSx2oBUVWI6cHv+e52CNvuOr9p2\nD/+XgeZqA79ut395v3/dZXt9x7asvwNNubfAysD2dcYCq5Lqj19BGoR3BLBqj//PrYH357/XIHW4\nusXMJqUJrEr6cXkJcHgPbT6DNGDz9Hz7BaQzRSWxN7T8r1YCbsh/dxu7cHbJshrjdyQN7jyWVBni\nZuBtvTxfw76QBpt/E/hdvr0KcFmP23gpcDgpt/00YOe6YxufE8CHgU/kv4vGYjCRt/x7es95vjZf\nfxf49/z3VXW3m+N/QJqD4A+kiidPA+YOIfYq4OlNt9fo4TH3HevLcC4j3wFffKlyAU7IH6pfJQ08\nvIF0Sv2yxpfDJHE7kCY52S7fvqnHdvuOr9p2YRsDHSjaqWNAlwFaVWJb1l3oS770cVaJzesu3+f/\ncn/SQMDf59trAheUPgfAbsCBffyvTid1Sq/Kt5fqoaOzK6kD+wPS6fub8n4sD3xlkpgppB8NV5E6\nkqvmy3Tg+oI2K8U3baevTgdpspOnNd1+GvCyHp/r7fLn0FeBN/UQ1+jQNg8476ujBKxO+gHyRN2x\npB+nLwcuBl6Yl5W+xp5JOmjyinx7GvCewtgvAtfn9pfOz3HR4OQq7eb1lyDldq+cb68GbNR0/wtr\nir2mzbZK/9d9x/oynMvId8AXX6pcSB3aFZpur0A66rEccF2X2OVJR4J+Bvypj7b7jq/adsH2B31E\n/MYO6096X9XYlnUX6oj28GXUVywp5ek64LZ8e2PgWz3s85Wk08FXdNqXdvtGqkN9BrB5aVxT/GX5\nurnd4soxue0358uaBevvxUQefKPazM2kjvGedcdP9pyWdjpInTq1xBW/h0idw7OBXfLlTAordJDO\ndizJRId8DXr7kbgiqf736aSDEl8CXjKE2Ffmz6/G0fB1SQMfi/a7y7Yv6nL/qqR0usZn6TOb7ntd\nXe0WxNdSNYqUWvIrUhnE9+Xn60uF2+071pfhXFw1xcbd01mwQsFjpJriD0nqWLkgIh4E9lGa6KDn\nmQSrxFdtu8CgJxM6S9IhwGcjf7rn3NYDSaknnVSJbTZH0uGk0/gAHyLVya4z9mukOsA/B4iIqyS9\nsnPIAh6NiGhUpckD3EocRPryPD8iLstjISad+bSNByWtRs5/lrQFqfLJpCRtEGkCnUa+7O35+pmS\nnhkdKq5ExNeBr0v6cEQcNdl6dcU3+WXO6//vfPsdpI5HN2q8NvP+PJlz40u9gZQO9ySApGNJnftP\nF8QeCfwUeHp+n7yNNJCv1FWkuusHRe/VRKrEzgOeBHbKFUgaA2wHoWNd8GgaS5I/Sx9suvtLpB9C\nA2+3wMAncQOIiI9LeiuwVV40MyJ+WrLRKrE2HO6I27ibBVwiqVFWbVvghNzhua5ToCam/r1K0ka5\no1ZcyqpK/ADaXmDGxsbymJixcdLyh33GfoyUk3mjpCvzso2BOaS0hU6qxDb7MGkcwImkL/wzSR3q\nWmMj4nYtWB6zlwlyTpJ0NLCypN1JR0u7TvQSESczUbWEiLiJ9JyV2of042F9SReQjrJ2m659H2AP\nUorX/Kab/u5acSUijpK0JQu/to4r2ekBxH9cabKXrfOi0k7HTbma0Lfz7Q+Sjsz3YmUmyiyuVBoU\nEbPywMfXkDpj20fE7xr3S1olOg+iXK/5R0SPqsTOAv4L+C2pQz5IVTr0VTrDVX9IVInvGBsRPyHN\nKtr7hivEWv1cNcXGnqTNmPi1f0FEzCmMay5l9QPS4LZeSln1HT+Atn/JxIyN8zuGEXHYpEGDiV2P\niWmir80dxCJVYgu3f6+408gAACAASURBVFS0TFVeNVbSj0kpRN8AXkZKodgsIt7Zw7ZfR5qlUsCv\nIqLr0TpJawC7s3CHdJce2l0KeH5u94YoKCuX43YEfhmprvd+pLzWz3c6It4U+yNgfVJKTuO1FRHx\nkcK2q8avC/wlIh7Ot5cjnSG7pUvc00lHpv+V1CE6G/hoFJZ5k7QTKT3lN6T/9ytJNaJPLInvsu2O\nVT2USgf+Fwu/VkpKVVaJPT8itu62Xj96qWSyqMTW2Xb+cfkl0hlg0VtVnr5jbTjcEbfF1gBKWVUp\nwVW17d9Gn9Ni9xOrVPv2adFS8kqpbvt9nTqXVWJ73MeBfwlKWp1UNeW1pC+wM4C9IuLewu322zG8\nkFQysPXHUtFRLUkfAmZFxD/y7VWAnSLiWwWxjTJrWwOfJw0+/FxEvKwg9nfAC/o9yjqA+DmkWRYf\nzbeXIf0437xzZH8kbRURF+SzTKuyYNnEOwfURscSd5KuAr7Dwq+VrqlXFWNfA+xE+tHSPGnUKd1i\nC7ZdZT6CKp8Dk81lMP95jjRh02TxF0fEFn22PWmspBuBbZvPlPSw3b5jbUhiEUhU98WXUVyoXsqq\nSgmuqm3PBF7U5+PuOZaUrrJGm+Wr031gVd+xPe5jLQOlKu7THJqmaAeWoaA0HT0MrCyNp7zCTKNi\ny6HAu3qMPZlUA7vf/a4a3+5xl1RNOZZczSLfXgX4fkHc3DpfPyXbprAEXg2xx+fX97FMTDff9X/W\nZjsLlQMFNqywX6d0uG9dYErT7eWA6d3aHcTzTMrDb769JOnHckls10pLdcT6MpyLc8RtcdaY+nfX\n6G/q3yrxVdveGnif0kyExTM2VohdNtrMwBcR9xQMQKwSOxKSjqJDzmYUpkoAS0U+OpvjHs1Habs5\nTdIbIuIXhe20WlLS/AGIkpYk/QgocUfOa38d8KV8tHeJwtjVgeskXcqCR0m3G1L83ZK2i4ifA0h6\nM3BPQdxGkc8e5Pb+LqnkiOxjSjORri3pyNY7e3idVHGqpA+SBnw2/89KJseqErt5RDy/lx3NY2O+\nS8op34U0Q+Z6+T2xY+QBoxHx2w7beDspdeoBSZ8lpU4dHDl1KiLe0mEXTiZVQmp4Ii/bvEu7jed5\nrQrP87MlfToiDs3vqZNIA3pLzJF0Imlgba9nH6rE2hC4I26LrUinjg9vun0bqY5u7fFV2wb+o4d1\nBxG7oqSlIs++2SBpadJRpbpiezHIigWNcQZbkSbDaeT6vp0ug4Bb9Nsx3Av4jFLln8foPa/zl8CJ\nuUMN8J95WYkdgX8HvhoR/5D0LODjhbEHFK5XV/wHgFmSvpFv/wl4d0HcEs2DIiWtStn345tIaUv/\nRnkFn151e12/N183P0cBrFew7SqxF0p6QUT08n74Gun1tQLwf6SBqecrVeo5iomxPp3sFxEn59Sp\n15IOYHybNIajm35/GA/ied6F9Nr8NPBq4BcRcURh7IqkKjWvb1oWQElnukqsDYFzxG2x0xhkpIWn\nty7q7FSJr9r2qEj6Imm2xj0jlQtD0gqk/Ol7IuKTdcROsr0VSf+rB1qWvy8ifjjIWKXp3rdu/IjI\nPx7Oi8I8UEnrk6pLrEl6jm8nTSByY0l8vyQtQaqA8tq86Ezgu5Gnra+57XWA50bEWZKmkuo9P9At\nblDxeRsrAETEP1uWvzcijm2z/ntIs/OeRHqe3gYcEhE/KmhrSdK4gcO7rTtJ/PqkuQQekfQq0gDu\n42Iiv3/VwiPUQ5Xz+ddnov571zNrzTnYkn4XEf/SdF/p+JorIuLFkg4lpfOdUJpTLulM4KiWH8Yf\niYjXdIvN628cEVeVrNsU0/yYlgaOJqXrfQ8gCgZB21ObO+Jm1pVSBY6DSeUGbyV96T6b9GWyX3So\nyFEltmU7mwPfJ816KOAfwC5RNrCsr1hJNwAvb3SE8qDHi/s4Jd+2Y9glZhXguTTVNo6Ic3tpN29n\nVWDtiLi619g+2tqd9ANg1YhYX9Jzge/00NGpFF+w/U6VKbYENiP9QJ4TPdTVlnRpRLy0z326Mrc7\nHfgFaZKcF0bEG7rE/WtE/FqpKsZCOqUeVIlt2sY6k8ROWoJV0lURsXH+e/uI+N+m+4oGkUs6DbiD\nlDq1KfAQaXDsxgWxjR/Ga+VFPf0wljSFNPPsC1nwfTlpNSNJv+mwyYiyCjU/oE2qXKd2BxFrw+HU\nFFvsSXoOqZTgddHfqPS+46u2PSz5iPCnJB0IPCcvvjEiHqoztsX3gA9GxHkA+dT0D0hHEOuK/SJw\nRf4ybZSlO6B0hyXtldt5ADgmHx37VESc0SVuN1J6ytqkUn5bABdRUMs7x88mTbm+FOlU+l8lXRgR\ne5fue58+BLyUNFskEfEHpdKAw4rvpm2aR36ediOdrhdwtKRjonxyoQtyOsyJNE0uU3i088mIeFzS\nDqSjtUdJKskd3oY0Ida2be7rlnpQJTatVDjnQYv9JE2NiHktnfD1KU/N6zt1KiL+CGzRzw/j7EfA\n9aQUlYOAGUDHz+2IeHWPbbRzWtPfU4AdgD8PIdaGwEfEbbGTO1VvjzRY8N2kiV7OJeUYzuz25Vsl\nvmrbozLZkbOGLkff+o5t2c5Cp597PZ3dZ+wzmcg/vSR6KEvXOAKoVMLxA6QZE3/UrV1J15AGkF0c\nEZtI2gD4QnQeiNYc3zh9vxvw7IjYX7ksYem+90PSJRHxsqb2lyJVmihqt2p8wfYnK1N5NenMRyN1\nanlSRZ/S/W531LP0aOclwBGk1JhtI+Lm0qPDhfvWNh2n7ti65R9ozUelbyuIWQ3YnzRgPYDzSdVM\nSsuRNl6XjRKfXVPVJIk0tiSAH5N+TL+Z1IE/OvJsrL3IqWfnR8SWXVceYKzVw0fEbXG0RsT/b+/O\nwySryjuOf38DKNswA1ExqKgoQgYFQRQUxBAVJTwuUZCAoMElRIyOWxLAJOISiPtClKCigCgGFwTE\nQUBZdUYFBAYUYmRUSIIkMMIEZRt++eOcsquX6rpVp+6t7rnv53n66bm3673nTg9Nnzr3Pe/rzoa5\nN5N+Ad+e81FXkDYN1RVfOva4vIjpjzc7q4v9VtBKYrtdorT58PQcdyBwcScHs8/q40Cx6t3ufStJ\nWw2Q19n5e/4pKe/3+vyLuZ97bN8jCaXaxTdIGiQdZv28UvgK0gSvKZdIOhrYSKmR0RHAOQ3G99Pr\ney8md0xdO8trpylc9TyM9Cbtn/Ik/PGklddRWUoqMdh07DSSduykSOVJ7N+RnoBcR6p88tsK13gx\nqfvrVsBtwNakVeodZovLvkxa+Oh0qX0l6SnG83pGTNZJo/uNpCcDt5Ia5czmk/k1DyFNwB9K6nq7\nH7A96Xs8qG0rjFtHbKhBTMRDG90v6VG2/xP4PyYeJd9Lqu1aZ3zp2OMytayXgf8hraysqjG2WycH\n9F1Tzu+crznb6uOgscXt3rMrJZ1Pql98lKSFVGsHfoukxaSSYxdIWk3Kr6/qPcC3STWEf6TU1fRn\nA8QP60hSDu1KUqWWb5HK1TUV38/3epz/PPADSWfm45eSN9NVIWlL4FhgK9v7SlpCepPd9xpOVUfe\n3HW8itQJcVRGWU2o1MmkvG5IaV9/QPoZeympsdCrKlzjvaRUrQvz6vTewCEVx/9D2+/tOn6fpAMr\nxgJ8Wmnvxt+TJtObkp5qzubZtp+S33jcmu/hPkmnA5Xe0Gtic7/y51tJb2JqjQ3NiNSU0DpKlQk+\nCXyN1A1vF9KkZU9SC/IP1RVfOva4SJo6gYV0/y8AjrH95Tpix00F7d5z/ALgqaRqCQ8l1cl+1CAp\nSJKeAyzK93Ffv9ePi1L1kFNtv3Ic8fka3Tn5nyW90eqbk59jdyH9HEJKN6ha4xlJy/K478ypSOuT\nmiA9pULsKmbeTFelhGCVextby/cZrtddNeVqUi3y+/NTomuqpAJJusL2rkpdQXe2/aC6NoH2if0I\n8ENSdRxI1XGeYfsdFWIXAPvbPqPfa6fEdf+dz7P9wq6vXW37qYNcL6x7YiIeWknSIlJDnSeRngzd\nApxl+4a640vHnkuUKnJcOMwv60Fj8/ftXaQNk5C6k77H9p11xaqg3XuOn3HT5Wy5w3lCer3t7auM\n0eMaTyLVVt7S9pMl7Qi82Pb7hr1mxXEvB/5k2DcMI4jvzsk/nLRa2Tcnv5SkH9l++tSJZpVJVs5b\n7tiQlE+8he1/HNG9lbSLHzq2x/VuAt5OahD1Pk8uX1h1Mn0haQX9ONIb29tIE/q+Oc95dXgTJtKQ\n1mPiqaTdv3TtFbZ37TfOlJhlpH1BU0tpPhI42xWr7eSf4cfRlcng6vtrho4N9YvUlNBKeQJ2wmyv\nkXS87TeNOr507LnE9h0Vc55HEfs5UprLK/LxoaRVyCobGIeN7fzC3g/4jO1zJQ0ymV3KxKbLvZU3\nXc4WYHutpBslbV1lA1oPnyFVkjgxX/NaSV8ilZGs002kCiJnM7l6SNUa26Xx3Tn5XxggJ7/U3XlC\n3elkujvQ9w0iwAwbBT8m6UpgJBNxeqfj1B07k0tI1XwAlkva0vav86S0SqMrSHnWvwPeSsrxXkRK\nxerL9sK8ADC1LOglFce+UNI7mF4dp2edd9u9GqitITUK6kvS50gVnq5nIrWt0v6aktjQjJiIh9Bb\nlS5vdcWXjt2InJ+5uqHYJ9h+edfxu/Pj7TpjS9q9w/CbLjcHrldq9d79C79qq/eNbf9wyhz0gV4v\nHqGf548FpJrtTccPm5Nf6m2knOFtJH0PeDgp7aEvTW74soBUU7zy72ZJxwIf8EQDoM2Bt9v+ewDb\nfz1DzNtmu2bnjc9MsSVsH5bH34i0EfcESZ3qJfv1i89Pi76ZN8c+yIAbSXs8ofo+ULVO/YGkSewR\nU873TCPS9I3fU91WYdzdbS+peI+jjA0NiIl4CKEvpXJ6U/PYtiDVo511g1VJ7BS/k7Sn7cvzdfcg\nrYzVGVvS7h2G33S5IZNXy8RgG/j+V6k2c2eFdn/gvweIH4rtd+fxNq5SAWPU8aSNnk8FbrL927xK\nfdgQ1xnUT4AzSa3E15D+vf+9Ymz3huAHgF8w8eSmin1tH905sL1a0p+SNhT2MsybnFE6BbiLiSpR\nB5PqiM/6985Pix6UtKhKStoMBn5CNcUS0iS8U/7wMtIm09l0b/ye1k2Zahu/l0takjf2DqokNjQg\ncsRD6KF0o9Jc2iRVStO76Bm43bnucl2xU66zE+mX9aJ8ajXwalfoGFkSOyoaYNPlTP/+GqAOuFKV\nlE8DzyL9XVcBr/RwTVgqk/RMUrWRTW1vnb/vh9ueuoJYV/xeM533EB1JByHpDNLE8ov51MHAYtsH\n1DluHvtaUo70vfl4I1Jn0Crl/MZC0k+mrtLOdK5H7FmkTbgXMPlp0Zt7Bk3EdnL5rwZ2s32vpOur\nfq96/Dsvst33jVPXU4DuSfwJtu+pEPsc0hOXW0kVtkTKaa+yuXXo2NCMWBEPobfS3NK5VDasSMkE\nboSTv7vyRrzN8nXvUqq5XHfsSFTJQ5X0BtIv623yBKtjIRXydaekHHwLuIiU7nA3qXZy1VzrYX2M\nVA3nbADb1/SaHNcU3/20YkNSjeorqV5uclhPnjKJvEhSpRVIFWxCzr4IfEeplTmkJwCVUjY0RMv2\nEblK0u62V+T72A24omLsZUzPb666wl9aFnTof2cmngJ8Ih9XegqQnUTa17KSwVOtSmJDA2IiHkJv\nHx9jfOnY66KvAbvYvqvr3FeBp9Uc26QvActIFSGO7Dq/ZrYNYV06E5LtSI/gzyK9qTuUVLatdrZv\nnpKbvrbXa0cdb3tSy3ZJjyFN7utWMrEs2YSM7ffnN22dPOf32v52xbEHbtk+Ik8Dvi+psxl5a+DG\nThpbn9Xag4FX2b4OQNJBpO/Zv/Qb1Paf5T8eo9QNdRFw3gD3XfLvXDKJ/x/bZw9wn6OKDQ2IiXho\nHUnnMEPd3o7OhjjbJ486vnTsNsp5nDsAiyR1T042o2sVb9Sx45BXQe8EDhoyvpNjfSnpjceafHwM\ncO6IbnM2N0t6FmClBiZLGWxiVxo/1S3AH/V9VbmSiWXJJmRIAywjvYEb1BNtHyDpJbZPUaqsc9kQ\n1xnUC/u/pKf9ga9KOhh4NmmfyT6DXmSASind+1w2YOLf2cBjSW9kqiiZxP84/9ucQ0ovASqXICyJ\nDQ2IiXhoo07TnJcBjwROy8cHAb+uOb507DbajrRxcTHQveK5Bnh9jbHz2ZZAdx76fflc3f6K9DTn\nUcB/AucDb2wqXtLxTLzR7TRTqtR8qVDJxHKojcSSLre9pyY6J/7+S1SoiZ0N07K9WGGq202S/pyU\nXvIrYB/bVTdtD6tSmcE+St6sbUSaRHe/4ahagrAkNjQgNmuG1tIMzRlmOldHfOnYbSTpmbaXNx07\nH0l6JynVobtl+7/ZPq6m8d5v++8kHWD7K03Hd13n1V2HDwC/sD3qWtgjNWUjsYA7gL+wfU0DY7+O\nlLb1FFL7+U2Bf7B9Yt1jD0rTqy89gvT06F6Aub75cIZN65PUvZE6zF0xEQ+tJemnwH62b8rHjwe+\n5a5ub3XFl47dRpI+QGpI8ztSXueOwFttnzZrYGHsfKVUt/jZ+fBSD9CyfYixVpK+p1cOU+2nNH7K\ntR5C6loLcKPt+2d7/VzRvZF4wLgv2D6037kZ4oZq2T4ubZ7I5t8Pb2J6d8y+fQVKYkMzIjUltNlb\ngYuV2i6LlO/3lw3Fl47dRvvY/ltJf0aqtfwy4FIm0nvqip2XbF9FM2kZkN7crAY2lXQXEzWSq6ZJ\nlMYDIOmPSdUpfpFjHyPp1XWXLyyh1CTq5eSJUmejqu1K3SJJeyC6r7c+FTYh235Q0t8C82Iivi5P\ntCv4Bqn6yTkMXvmkJDY0IFbEQ6vlX4Lb58MbOrV4m4gvHbttOvV+JX0W+Krt8yRdY3unOmNDf0pd\nQ++VdJbtlzQd33WdK4GDbd+Yj58EnG57rlXH+T1J55FSLK6kq0KM7Q/3DEpxRwFHk3KAO82PRNoP\n8GnbR1UY+59JreUrt2wPzZP0A9u7NR0bmhEr4qG1clWGw5mo33uxpBOrPsouiS8du6XOkXQDKb3k\nDZIeDvRthjGC2NDfcmAXUp3kccR3bNCZhAPY/vf8szaXPdr2wJs9c77/cZKOqzLp7uHA/Ll7Q6yZ\npWV7GIuPS3oXafNyd+WTKk+8SmJDA2JFPLRWXh3dgInmF4cCa22/ru740rHbStIWwJ1Ora43Bjaz\nfWvdsWF2kq4jtQp/L5Ob6gD9S6WVxndd53Okx++dlKNXAuu5/gY1Q5P0aeB42ysLrrE5sC2Tm/LM\n2XScMBhJx5F+R/ycifQS2+7bqKokNjQjJuKhtWZKTRgkXaEkvnTsNsorm29gcgfCfx3gCcRQsaE/\nSXuSJr2vIHfF7OJ+E+HS+K7rPJS0urtnPnUZ8Km5nPal1NTlicAqhmhBniufLAUeDVwN7A4srzhJ\ne9VM522fWu3uQxMk/QewxPZ9fV88wtjQjEhNCW22VtITbP8cQNI2DNYFsCS+dOw2OoH0FOFT+fjQ\nfK7KU4SS2NCHUw3sy3MJzpOaju+6zr3AR/LHfLFvYfxSUhfVFbb3VmpidWzF2Kd3/XlDUnfOq0jl\nFMPccR2pF8JtDceGBsREPLTZ35DaDHdXLjmsofjSsdvo6VOeGHxXUtVayyWxoQ9Jf2L7u8BqTe5g\nClRKTSmK77rOHsAxpJ+n7lJtczbn2fYv8xOBbW1/Pu9f2HSAS9xj+x5JnU2vN0jaruLYb+o+lrQY\n+PIAY4dmLAZukPQjJud5VylBWBIbGhAT8dBatr8jaVtS90VINYcrP8IuiS8du6XiCcTctRfwXVL3\n0t+XHez63G8iXRrfcRKpNOikCiRzWd5Ityvp/wWfJz25OQ3Yo+IlbskT6G8AF0haDQxb6u9u4PFD\nxob6vGtMsaEBkSMeWkvSesB+TG90UOmxdkl86dhtJOm5pInKTfnU44DDbF9UZ2zoT9LbmT6BJv+5\n73/XpfFd15l3pdokXQ3sDFxle+d87tqqOeJTrvUcUofO86rkBEs6h4lulQuAJcAZto8cdOwQwnBi\nRTy02TmkEnYrGa7RQUl86dht9D3gRFIe62+Ab5PK3tUdG/rrpFJsR8o7Pos0mX4R8MMG4jsukvRB\n0gr6fCnVdp9tSzKApE0GvUB+Y78lacMnwCOBX1UI/VDXnx8Afmn7lkHHD/WQtIaJN0qTvkSfRlcl\nsaFZsSIeWmvYVadRxJeO3UaSziDVmf5iPnUwsNj2AXXGhuokXQrsZ3tNPl4InGt7r9kjRxY/0xOO\nOV2qTdI7SKUHnw8cB7wG+JLt4yvGv4mUfvBrJpenq1p15bGk/PQLJW0ErN/5/ocQ6hcr4qHNlkna\nx/b5Y4gvHbuNnmx7SdfxRbn0W92xobotSZ0dO+7L5xqJt733AGPNFQYuJ71RfBLwj7YvGCB+KbCd\n7dsHHVjS64G/BLYAnkAqgfivpCdHIYQGxEQ8tNkK4ExJC4D7GfyRXUl86dhtdJWk3W2vAJC0G3BF\nA7GhulOBH0o6Mx+/FDi57nhJh9g+TdLbZvr6HN97sSlpFfwOUqv5aweMvxm4c8ix3wg8A/gBgO2f\nSXrEkNcKIQwhUlNCa0laBbwEWOkhfhBK4kvHbiNJPyXlEHdyX7cGbiTlts76KL4kNgxG0i7As/Ph\npbZ/XHe8pMNtn5grkExj+92D3MM4SNqR1HL+5cAttp/X5/WdNx07kP7bPpfJefFVNo3/wPZukn5s\ne2dJ65M2jcbPQwgNiRXx0GY3A9cVTIRL4kvHbqMXjik2DCBvjBx6c+Qw8bZPzJ/n/IR7FrcBtwK3\nA1VWpRfmz7/KHw/JH4O4RNLRwEaSng8cQdpIHkJoSKyIh9aSdDKwDbCMAVeSSuNLxw4hTCfpFGCp\n7d/k482BD9t+zXjvrDdJRwCvAB4OfIVUPnDg/QuSNiM93am80TKnxr0W2IeUHvdt4LOxQBBCc2JF\nPLTZqvwxzEpSaXzp2CGE6XbsTMIBbK+WtPM4b6iCxwBvsX31MMGSdiXVyF+Yj+8EXmP7yn6xth8E\nPpM/QghjECviIfQg6fipLaCbii8dO4Q2knQN8Me2V+fjLYBLbD9lvHdWH0nXAm+0fVk+3hP4VJU8\nb0l7AMcAjyUtzHU2jW9T3x2HELrFingIvVVtMV1HfOnYIbTRh4Hlkr5CmlTuD/zTeG+pdms7k3AA\n25dLeqBi7EnAW4ErgbV13FwIYXYxEQ8hhLBOsH2qpCuATgOflw2Tbz3PXCLpROB0Uk3yA4GLc/WZ\nfl1F77S9rIF7DCH0EKkpIfQg6Srbu4wjvnTsENpI0tYznbddpd37vNSjm2jHrF1FJf0zsB7wdSZv\nGh+66k0IYTCxIh5CbxpjfOnYIbTRuaRVYYCNgMeT6sXvMLY7qllhN9Hd8udduy/JxBOFEELNYiIe\nWk/SxrZ/O8OXPl53fOnYIYQJUzdl5vSMI8Z0O42QtBh4FfA4un6n235zv9jCSXwIYQQiNSW0lqRn\nAZ8FNrW9taSdgMNtV/rFXRJfOnYIoRpJK9fxqinfB1YAK4EHO+dtnzJLzCG2T+vqzjlJ9DMIoTmx\nIh7a7KPAC4CzAWxfI2mvhuJLxw4hTDFlYrkAeBrwX2O6naZsaHvGCfUsNsmfF876qhBC7WIiHlrN\n9s3SpHTsgUp4lcSXjh1CmGYhEzniD5DatX9tfLfTiC9Iej3wTSZvuLyjV4DtEyWtB9xl+6MN3GMI\noYeYiIc2uzmniFjSBsBS4KcNxZeOHUKY7lvA0UzOlz4S6NvcZh67D/gg8E4m3oQYmLUpj+21kg4i\nPZ0LIYxJ5IiH1pL0MNKmyOeRqpScDyy1fXvd8aVjhxCmk3Qj8A7gOibnS/9ybDdVM0k3Ac+w/b9D\nxH4U2AD4N+DuzvkoXxhCc2IiHkIIYZ0g6XLbe477Ppok6XzgpT2qL/WLnakG+ay1x0MIoxUT8dBa\nkj4xw+k7gStsn1VnfOnYIYTpJD0XOAj4DpPzpb8+tpuqmaQzSXXSL2Ly37lv+cIQwvhFjnhosw2B\n7YGv5OOXA6uAnSTtbfstNcaXjh1CmO4w0s/VBkykppjUOXJd9Y38MTBJWwLHAlvZ3lfSEuCZtk8a\n5Q2GEHqLFfHQWpJWAHvYXpuP1wcuA/YEVtpeUld86dghhOkk3Wh7u3HfR5MkvQg41/aDfV88PXYZ\n8HngnbZ3yv8f+vG6XHc9hLlmwbhvIIQx2hzYtOt4E2CLPDm+d+aQkcWXjh1CmO77eVW3TQ4Efibp\nA5K2HzD2YbbPID89sP0AUUY1hEZFakposw8AV0u6mFS5ZC/gWEmbABfWHF86dghhut1JP1erSG9o\nRdp8uM6WL7R9iKTNSLnxJ0syaZX7dNtr+oTfLekPyGUPJe1O2qsSQmhIpKaEVpP0h8Az8uGPbA/U\nha8kvnTsEMJkkh470/l1uXxhR55QHwq8hdST4InAJ2wfP0vMLsDxpM2e1wMPB/a3fW39dxxCgJiI\nh5aTtDmwLWnzJAC2L20ivnTsEEKQ9GLSJtUnAqcCp9i+TdLGwE9sP26W2A2BvwZeAKwBlgPH276n\n9hsPIQAxEQ8tJul1pI6WjwauJj3WXl61hm5JfOnYIYQAIOkU4KSZ3sRLeq7t78wSewZwF/DFfOpg\nYLHtA2q52RDCNDERD60laSXwdGCF7afmjU7H2n5Z3fGlY4cQQilJP5laoWmmcyGE+sRmzdBm99i+\nRxKSHmr7BkmDlD4riS8dO4TQYpLWkDdZTv0SaYPqZhUuc5Wk3W2vyNfcDbhihLcZQugjJuKhzW6R\ntJjUDOMCSauBQTZ1lcSXjh1CaDHbC0dwmaeRSj7+Kh9vDdyYn9it09VmQpgrIjUlBEDSc4BFwDLb\n9zcZXzp2CCEMo1eVmY42VJsJYdxiIh5aS9IXbB/a71wd8aVjhxBCCGH+i86aoc126D6QtB7pUW0T\n8aVjhxBCCGGe/D5M0gAAAaZJREFUi4l4aB1JR+WNTjtKuit/rAFuA86qM7507BBCCCGsOyI1JbSW\npONsHzWO+NKxQwghhDD/xYp4aLNvStoEQNIhkj7Sb/PSCONLxw4hhBDCPBcT8dBmJwC/lbQT8Hbg\n56QW0U3El44dQgghhHkuJuKhzR5wys16CfAvtj8JDFKbtyS+dOwQQgghzHPR0Ce02RpJRwGHAHtJ\nWgBs0FB86dghhBBCmOdiRTy02YHAvcBrbd8KPBr4YEPxpWOHEEIIYZ6Lqikh9CBpue1njiO+dOwQ\nQgghzH2xIh5CbxuOMb507BBCCCHMcTERD6G30sdFJfHxqCqEEEJYx8VEPIQQQgghhDGIiXhoLUnv\n73NOdcWXjh1CCCGE+S8m4qHNnj/DuX27/nxojfGlY4cQQghhnos64qF1JL0BOALYRtK1XV9aCHyv\nc2D7ulHHl44dQgghhHVHlC8MrSNpEbA5cBxwZNeX1ti+o8740rFDCCGEsO6IiXgIIYQQQghjEDni\nIYQQQgghjEFMxEMIIYQQQhiDmIiHEEIIIYQwBjERDyGEEEIIYQz+HxlthRCUD9GhAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hat_beta_cox_EL = np.array(coeffs)\n",
    "\n",
    "fig = pl.figure(figsize=(12, 3))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.stem(hat_beta_cox_EL, 'b', markerfmt='bo')\n",
    "ax.set_xlim([-5, len(hat_beta_cox_EL) + 5])\n",
    "pl.yticks(fontsize = 18)\n",
    "pl.xticks(np.arange(len(hat_beta_cox_EL)), X.columns, rotation=90)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TODO: \n",
    "essayer de faire tourner un Cox classique (CoxPHFitter) sur le support du $\\beta$ Elastic-Net (i.e. en gardant que les variables dont le coeff est non nul) et voir si ça améliore les prédictions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# C-mix model \n",
    "Il faut que tu installes le package en suivant les indications ici https://github.com/SimonBussy/C-mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Launching C-mix...\n",
      "\n",
      "Testing l_elastic_net=5.71e-05 on fold  0 1 2 3 4: avg_score=5.92e-01\n",
      "Testing l_elastic_net=7.85e-05 on fold  0 1 2 3 4: avg_score=6.14e-01\n",
      "Testing l_elastic_net=1.08e-04 on fold  0 1 2 3 4: avg_score=6.59e-01\n",
      "Testing l_elastic_net=1.48e-04 on fold  0 1 2 3 4: avg_score=6.14e-01\n",
      "Testing l_elastic_net=2.03e-04 on fold  0 1 2 3 4: avg_score=7.01e-01\n",
      "Testing l_elastic_net=2.80e-04 on fold  0 1 2 3 4: avg_score=6.30e-01\n",
      "Testing l_elastic_net=3.84e-04 on fold  0 1 2 3 4: avg_score=6.72e-01\n",
      "Testing l_elastic_net=5.28e-04 on fold  0 1 2 3 4: avg_score=5.94e-01\n",
      "Testing l_elastic_net=7.25e-04 on fold  0 1 2 3 4: avg_score=6.60e-01\n",
      "Testing l_elastic_net=9.96e-04 on fold  0 1 2 3 4: avg_score=6.13e-01\n",
      "Testing l_elastic_net=1.37e-03 on fold  0 1 2 3 4: avg_score=6.85e-01\n",
      "Testing l_elastic_net=1.88e-03 on fold  0 1 2 3 4: avg_score=7.09e-01\n",
      "Testing l_elastic_net=2.58e-03 on fold  0 1 2 3 4: avg_score=6.58e-01\n",
      "Testing l_elastic_net=3.55e-03 on fold  0 1 2 3 4: avg_score=7.09e-01\n",
      "Testing l_elastic_net=4.87e-03 on fold  0 1 2 3 4: avg_score=6.33e-01\n",
      "Testing l_elastic_net=6.70e-03 on fold  0 1 2 3 4: avg_score=7.28e-01\n",
      "Testing l_elastic_net=9.20e-03 on fold  0 1 2 3 4: avg_score=6.99e-01\n",
      "Testing l_elastic_net=1.26e-02 on fold  0 1 2 3 4: avg_score=6.82e-01\n",
      "Testing l_elastic_net=1.74e-02 on fold  0 1 2 3 4: avg_score=6.73e-01\n",
      "Testing l_elastic_net=2.38e-02 on fold  0 1 2 3 4: avg_score=7.30e-01\n",
      "Testing l_elastic_net=3.28e-02 on fold  0 1 2 3 4: avg_score=6.81e-01\n",
      "Testing l_elastic_net=4.50e-02 on fold  0 1 2 3 4: avg_score=6.87e-01\n",
      "Testing l_elastic_net=6.18e-02 on fold  0 1 2 3 4: avg_score=6.40e-01\n",
      "Testing l_elastic_net=8.50e-02 on fold  0 1 2 3 4: avg_score=6.80e-01\n",
      "Testing l_elastic_net=1.17e-01 on fold  0 1 2 3 4: avg_score=6.50e-01\n",
      "Testing l_elastic_net=1.60e-01 on fold  0 1 2 3 4: avg_score=4.97e-01\n",
      "Testing l_elastic_net=2.20e-01 on fold  0 1 2 3 4: avg_score=5.00e-01\n",
      "Testing l_elastic_net=3.03e-01 on fold  0 1 2 3 4: avg_score=5.00e-01\n",
      "Testing l_elastic_net=4.16e-01 on fold  0 1 2 3 4: avg_score=5.00e-01\n",
      "Testing l_elastic_net=5.71e-01 on fold  0 1 2 3 4: avg_score=5.00e-01\n",
      "Launching the solver QNEM...\n",
      "init: p0=0.00017446482120139393\n",
      "init: p1=0.04846762254402027\n",
      " n_iter  |   obj    | rel_obj \n",
      "       0 |  5.68608 |        1\n",
      "       1 |  5.57062 | 0.0203061\n",
      "       2 |  5.56575 | 0.000874007\n",
      "       3 |  5.56438 | 0.000245178\n",
      "       4 |  5.56369 | 0.000124105\n",
      "       5 |  5.56331 | 6.94415e-05\n",
      "       6 |   5.5631 | 3.74641e-05\n",
      "       7 |  5.56299 | 1.92554e-05\n",
      "       8 |  5.56294 | 9.60961e-06\n",
      "       9 |  5.56291 | 4.72578e-06\n",
      "      10 |   5.5629 | 2.30763e-06\n",
      "      11 |  5.56289 | 1.12186e-06\n",
      "      12 |  5.56289 | 5.46025e-07\n",
      "At the end: p0=0.00017015324074667003\n",
      "At the end: p1=0.048835172655645985\n",
      "Done solving using QNEM in 4.60e-02 seconds\n"
     ]
    }
   ],
   "source": [
    "n_samples, n_features = X.shape\n",
    "\n",
    "X = np.array(X)\n",
    "Y = np.ravel(Y)\n",
    "delta = np.ravel(delta)\n",
    "X_test = np.array(X_test)\n",
    "Y_test = np.ravel(Y_test)\n",
    "delta_test = np.ravel(delta_test)\n",
    "\n",
    "tol = 1e-6            # tolerance for the convergence stopping criterion \n",
    "eta = 0.3             # parameter controlling the trade-off between l1 \n",
    "                      # and l2 regularization in the elasticNet\n",
    "fit_intercept = True  # whether or not an intercept term is fitted\n",
    "gamma_chosen = '1se'  # way to select l_elasticNet_chosen: '1se' or 'min'\n",
    "warm_start = True     # at each L-BGFS-B iteration, reset beta to 0 or take \n",
    "                      # the previous value \n",
    "grid_size = 30        # grid size for the cross validation procedure\n",
    "metric = 'C-index'    # cross-validation metric: 'log_lik' or 'C-index'\n",
    "verbose = True \n",
    "\n",
    "## Choose between C-mix or CURE model ##\n",
    "model = \"C-mix\"  # \"C-mix\", \"CURE\"       \n",
    "\n",
    "if verbose:\n",
    "    print(\"\\nLaunching %s...\\n\" % model)\n",
    "\n",
    "learner = QNEM(l_elastic_net=0., eta=eta, max_iter=100, tol=tol, \n",
    "               warm_start=warm_start, verbose=verbose, model=model, \n",
    "               fit_intercept=fit_intercept)\n",
    "learner.n_features = n_features         \n",
    "\n",
    "## Cross-validation ##\n",
    "learner.cross_validate(X, Y, delta, n_folds=5, verbose=False, eta=eta, \n",
    "                       grid_size=grid_size, metric=metric)\n",
    "avg_scores = learner.scores.mean(axis=1)\n",
    "l_elastic_net_best = learner.l_elastic_net_best\n",
    "if gamma_chosen == '1se':\n",
    "    l_elastic_net_chosen = learner.l_elastic_net_chosen\n",
    "if gamma_chosen == 'min':\n",
    "    l_elastic_net_chosen = l_elastic_net_best\n",
    "    \n",
    "grid_elastic_net = learner.grid_elastic_net # get the cross-validation grid \n",
    "                                            # to plot learning curves\n",
    "\n",
    "## Run selected model with l_elasticNet_chosen ##\n",
    "learner = QNEM(l_elastic_net=l_elastic_net_chosen, eta=eta, tol=tol,\n",
    "               warm_start=warm_start, verbose=verbose, model=model,\n",
    "               fit_intercept=fit_intercept)\n",
    "learner.n_features = n_features\n",
    "learner.fit(X, Y, delta) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C-index on test set: 0.65\n"
     ]
    }
   ],
   "source": [
    "## Obtain the marker vector on test set ##\n",
    "coeffs = learner.coeffs\n",
    "marker_test = QNEM.predict_proba(X_test, fit_intercept, coeffs) \n",
    "c_index = concordance_index(Y_test, marker_test, delta_test)\n",
    "c_index = max(c_index, 1 - c_index)\n",
    "print(\"C-index on test set: %.2f\" % c_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Donc le C-mix améliore encore les perfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}