resnet5_keras_baseline_mode_v1_0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ResNet50 Keras baseline model\n",
    "\n",
    "This notebook takes you through some important steps in building a deep convnet in Keras for multilabel classification of brain CT scans. \n",
    "\n",
    "\n",
    "*Update (1):*\n",
    "* *training for 4 epochs instead of 3.*\n",
    "* *batch size lowered to 16 from 32.*\n",
    "* *training without learning rate decay.*\n",
    "* *Weighted BCE instead of \"plain\" BCE*\n",
    "* *training data lowered to 80% from 90%.*\n",
    "\n",
    "\n",
    "*Update (2):*\n",
    "* *adding competition metric for training*\n",
    "* *using custom Callback for validation and test sets instead of the `run()` function and 'global epochs'*\n",
    "* *training with \"plain\" BCE again*\n",
    "* *merging TestDataGenerator and TrainDataGenerator into one*\n",
    "* *adding undersampling (see inside `on_epoch_end`), will now run 6 epochs*\n",
    "\n",
    "*Update (3):*\n",
    "* *skipping/removing windowing (value clipping), but the transformation to Hounsfield Units is kept*\n",
    "* *removing initial layer (doing np.stack((img,)&ast;3, axis=-1)) instead*\n",
    "* *reducing learning rate to 5e-4 and add decay*\n",
    "* *increasing batch size to 32 from 16*\n",
    "* *Increasing training set to 90% of the data (10% for validation)*\n",
    "* *slight increase in undersampling*\n",
    "* *fixed some hardcoding for input dims/sizes*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pydicom\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import collections\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "from datetime import datetime\n",
    "\n",
    "from math import ceil, floor\n",
    "import cv2\n",
    "\n",
    "import tensorflow as tf\n",
    "import keras\n",
    "\n",
    "import sys\n",
    "\n",
    "from keras_applications.resnet import ResNet50\n",
    "\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "\n",
    "test_images_dir = '../input/rsna-intracranial-hemorrhage-detection/stage_1_test_images/'\n",
    "train_images_dir = '../input/rsna-intracranial-hemorrhage-detection/stage_1_train_images/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Helper functions\n",
    "\n",
    "* normalizing the pixel values between -1 and 1 \n",
    "* read and transform dcms to 3-channel inputs for e.g. ResNet50.\n",
    "\n",
    "\\* Source for windowing (although now partly removed from this kernel): https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_cell_guid": "79c7e3d0-c299-4dcb-8224-4455121ee9b0",
    "_uuid": "d629ff2d2480ee46fbb7e2d37f6b5fab8052498a",
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def _normalize(img):\n",
    "    if img.max() == img.min():\n",
    "        return np.zeros(img.shape)-1\n",
    "    return 2 * (img - img.min())/(img.max() - img.min()) - 1\n",
    "\n",
    "def _read(path, desired_size):\n",
    "    \"\"\"Will be used in DataGenerator\"\"\"\n",
    "    dcm = pydicom.dcmread(path)\n",
    "    slope, intercept = dcm.RescaleSlope, dcm.RescaleIntercept\n",
    "    try:\n",
    "        img = (dcm.pixel_array * slope + intercept)\n",
    "    except:\n",
    "        img = np.zeros(desired_size[:2])-1\n",
    "    if img.shape != desired_size[:2]:\n",
    "        img = cv2.resize(img, desired_size[:2], interpolation=cv2.INTER_LINEAR)\n",
    "    img = _normalize(img)\n",
    "    return np.stack((img,)*3, axis=-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Data generators\n",
    "\n",
    "Inherits from keras.utils.Sequence object and thus should be safe for multiprocessing.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class DataGenerator(keras.utils.Sequence):\n",
    "\n",
    "    def __init__(self, list_IDs, labels=None, batch_size=1, img_size=(512, 512, 3), \n",
    "                 img_dir=train_images_dir, *args, **kwargs):\n",
    "\n",
    "        self.list_IDs = list_IDs\n",
    "        self.labels = labels\n",
    "        self.batch_size = batch_size\n",
    "        self.img_size = img_size\n",
    "        self.img_dir = img_dir\n",
    "        self.on_epoch_end()\n",
    "\n",
    "    def __len__(self):\n",
    "        return int(ceil(len(self.indices) / self.batch_size))\n",
    "\n",
    "    def __getitem__(self, index):\n",
    "        indices = self.indices[index*self.batch_size:(index+1)*self.batch_size]\n",
    "        list_IDs_temp = [self.list_IDs[k] for k in indices]\n",
    "        \n",
    "        if self.labels is not None:\n",
    "            X, Y = self.__data_generation(list_IDs_temp)\n",
    "            return X, Y\n",
    "        else:\n",
    "            X = self.__data_generation(list_IDs_temp)\n",
    "            return X\n",
    "        \n",
    "    def on_epoch_end(self):\n",
    "        \n",
    "        \n",
    "        if self.labels is not None: # for training phase we undersample and shuffle\n",
    "            # keep probability of any=0 and any=1\n",
    "            keep_prob = self.labels.iloc[:, 0].map({0: 0.4, 1: 0.8})\n",
    "            keep = (keep_prob > np.random.rand(len(keep_prob)))\n",
    "            self.indices = np.arange(len(self.list_IDs))[keep]\n",
    "            np.random.shuffle(self.indices)\n",
    "        else:\n",
    "            self.indices = np.arange(len(self.list_IDs))\n",
    "\n",
    "    def __data_generation(self, list_IDs_temp):\n",
    "        X = np.empty((self.batch_size, *self.img_size))\n",
    "        \n",
    "        if self.labels is not None: # training phase\n",
    "            Y = np.empty((self.batch_size, 6), dtype=np.float32)\n",
    "        \n",
    "            for i, ID in enumerate(list_IDs_temp):\n",
    "                X[i,] = _read(self.img_dir+ID+\".dcm\", self.img_size)\n",
    "                Y[i,] = self.labels.loc[ID].values\n",
    "        \n",
    "            return X, Y\n",
    "        \n",
    "        else: # test phase\n",
    "            for i, ID in enumerate(list_IDs_temp):\n",
    "                X[i,] = _read(self.img_dir+ID+\".dcm\", self.img_size)\n",
    "            \n",
    "            return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. loss function and metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from keras import backend as K\n",
    "\n",
    "def weighted_log_loss(y_true, y_pred):\n",
    "    \"\"\"\n",
    "    Can be used as the loss function in model.compile()\n",
    "    ---------------------------------------------------\n",
    "    \"\"\"\n",
    "    \n",
    "    class_weights = np.array([2., 1., 1., 1., 1., 1.])\n",
    "    \n",
    "    eps = K.epsilon()\n",
    "    \n",
    "    y_pred = K.clip(y_pred, eps, 1.0-eps)\n",
    "\n",
    "    out = -(         y_true  * K.log(      y_pred) * class_weights\n",
    "            + (1.0 - y_true) * K.log(1.0 - y_pred) * class_weights)\n",
    "    \n",
    "    return K.mean(out, axis=-1)\n",
    "\n",
    "\n",
    "def _normalized_weighted_average(arr, weights=None):\n",
    "    \"\"\"\n",
    "    A simple Keras implementation that mimics that of \n",
    "    numpy.average(), specifically for the this competition\n",
    "    \"\"\"\n",
    "    \n",
    "    if weights is not None:\n",
    "        scl = K.sum(weights)\n",
    "        weights = K.expand_dims(weights, axis=1)\n",
    "        return K.sum(K.dot(arr, weights), axis=1) / scl\n",
    "    return K.mean(arr, axis=1)\n",
    "\n",
    "\n",
    "def weighted_loss(y_true, y_pred):\n",
    "    \"\"\"\n",
    "    Will be used as the metric in model.compile()\n",
    "    ---------------------------------------------\n",
    "    \n",
    "    Similar to the custom loss function 'weighted_log_loss()' above\n",
    "    but with normalized weights, which should be very similar \n",
    "    to the official competition metric:\n",
    "        https://www.kaggle.com/kambarakun/lb-probe-weights-n-of-positives-scoring\n",
    "    and hence:\n",
    "        sklearn.metrics.log_loss with sample weights\n",
    "    \"\"\"\n",
    "    \n",
    "    class_weights = K.variable([2., 1., 1., 1., 1., 1.])\n",
    "    \n",
    "    eps = K.epsilon()\n",
    "    \n",
    "    y_pred = K.clip(y_pred, eps, 1.0-eps)\n",
    "\n",
    "    loss = -(        y_true  * K.log(      y_pred)\n",
    "            + (1.0 - y_true) * K.log(1.0 - y_pred))\n",
    "    \n",
    "    loss_samples = _normalized_weighted_average(loss, class_weights)\n",
    "    \n",
    "    return K.mean(loss_samples)\n",
    "\n",
    "\n",
    "def weighted_log_loss_metric(trues, preds):\n",
    "    \"\"\"\n",
    "    Will be used to calculate the log loss \n",
    "    of the validation set in PredictionCheckpoint()\n",
    "    ------------------------------------------\n",
    "    \"\"\"\n",
    "    class_weights = [2., 1., 1., 1., 1., 1.]\n",
    "    \n",
    "    epsilon = 1e-7\n",
    "    \n",
    "    preds = np.clip(preds, epsilon, 1-epsilon)\n",
    "    loss = trues * np.log(preds) + (1 - trues) * np.log(1 - preds)\n",
    "    loss_samples = np.average(loss, axis=1, weights=class_weights)\n",
    "\n",
    "    return - loss_samples.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b. Model\n",
    "\n",
    "Model is divided into three parts: <br> \n",
    "\n",
    "* (REMOVED) The initial layer, which will transform/map input image of shape (\\_, \\_, 1) to another \"image\" of shape (\\_, \\_, 3).\n",
    "\n",
    "* The new input image is then passed through ResNet50 (which I named \"engine\"). ResNet50 could be replaced by any of the available architectures in keras_application.\n",
    "\n",
    "* Finally, the output from ResNet50 goes through average pooling followed by a dense output layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class PredictionCheckpoint(keras.callbacks.Callback):\n",
    "    \n",
    "    def __init__(self, test_df, valid_df, \n",
    "                 test_images_dir=test_images_dir, \n",
    "                 valid_images_dir=train_images_dir, \n",
    "                 batch_size=32, input_size=(224, 224, 3)):\n",
    "        \n",
    "        self.test_df = test_df\n",
    "        self.valid_df = valid_df\n",
    "        self.test_images_dir = test_images_dir\n",
    "        self.valid_images_dir = valid_images_dir\n",
    "        self.batch_size = batch_size\n",
    "        self.input_size = input_size\n",
    "        \n",
    "    def on_train_begin(self, logs={}):\n",
    "        self.test_predictions = []\n",
    "        self.valid_predictions = []\n",
    "        \n",
    "    def on_epoch_end(self,batch, logs={}):\n",
    "        self.test_predictions.append(\n",
    "            self.model.predict_generator(\n",
    "                DataGenerator(self.test_df.index, None, self.batch_size, self.input_size, self.test_images_dir), verbose=2)[:len(self.test_df)])\n",
    "        \n",
    "        self.valid_predictions.append(\n",
    "            self.model.predict_generator(\n",
    "                DataGenerator(self.valid_df.index, None, self.batch_size, self.input_size, self.valid_images_dir), verbose=2)[:len(self.valid_df)])\n",
    "        \n",
    "        print(\"validation loss: %.4f\" %\n",
    "              weighted_log_loss_metric(self.valid_df.values, \n",
    "                                   np.average(self.valid_predictions, axis=0, \n",
    "                                              weights=[2**i for i in range(len(self.valid_predictions))])))\n",
    "        \n",
    "        # here you could save the predictions with np.save()\n",
    "\n",
    "\n",
    "class MyDeepModel:\n",
    "    \n",
    "    def __init__(self, engine, input_dims, batch_size=5, num_epochs=4, learning_rate=1e-3, \n",
    "                 decay_rate=1.0, decay_steps=1, weights=\"imagenet\", verbose=1):\n",
    "        \n",
    "        self.engine = engine\n",
    "        self.input_dims = input_dims\n",
    "        self.batch_size = batch_size\n",
    "        self.num_epochs = num_epochs\n",
    "        self.learning_rate = learning_rate\n",
    "        self.decay_rate = decay_rate\n",
    "        self.decay_steps = decay_steps\n",
    "        self.weights = weights\n",
    "        self.verbose = verbose\n",
    "        self._build()\n",
    "\n",
    "    def _build(self):\n",
    "        \n",
    "    \n",
    "        engine = self.engine(include_top=False, weights=self.weights, input_shape=(*self.input_dims[:2], 3),\n",
    "                             backend = keras.backend, layers = keras.layers,\n",
    "                             models = keras.models, utils = keras.utils)\n",
    "        \n",
    "\n",
    "        x = keras.layers.GlobalAveragePooling2D(name='avg_pool')(engine.output)\n",
    "\n",
    "        out = keras.layers.Dense(6, activation=\"sigmoid\", name='dense_output')(x)\n",
    "\n",
    "        self.model = keras.models.Model(inputs=engine.input, outputs=out)\n",
    "\n",
    "        self.model.compile(loss=weighted_log_loss, optimizer=keras.optimizers.Adam(0.0), metrics=[weighted_loss])\n",
    "    \n",
    "\n",
    "    def fit_and_predict(self, train_df, valid_df, test_df):\n",
    "        \n",
    "        # callbacks\n",
    "        pred_history = PredictionCheckpoint(test_df, valid_df, input_size=self.input_dims)\n",
    "        #checkpointer = keras.callbacks.ModelCheckpoint(filepath='%s-{epoch:02d}.hdf5' % self.engine.__name__, verbose=1, save_weights_only=True, save_best_only=False)\n",
    "        scheduler = keras.callbacks.LearningRateScheduler(lambda epoch: self.learning_rate * pow(self.decay_rate, floor(epoch / self.decay_steps)))\n",
    "        \n",
    "        self.model.fit_generator(\n",
    "            DataGenerator(\n",
    "                train_df.index, \n",
    "                train_df, \n",
    "                self.batch_size, \n",
    "                self.input_dims, \n",
    "                train_images_dir\n",
    "            ),\n",
    "            epochs=self.num_epochs,\n",
    "            verbose=self.verbose,\n",
    "            use_multiprocessing=True,\n",
    "            workers=4,\n",
    "            callbacks=[pred_history, scheduler]\n",
    "        )\n",
    "        \n",
    "        return pred_history\n",
    "    \n",
    "    def save(self, path):\n",
    "        self.model.save_weights(path)\n",
    "    \n",
    "    def load(self, path):\n",
    "        self.model.load_weights(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Read csv files\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def read_testset(filename=\"../input/rsna-intracranial-hemorrhage-detection/stage_1_sample_submission.csv\"):\n",
    "    df = pd.read_csv(filename)\n",
    "    df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n",
    "    df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n",
    "    \n",
    "    df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n",
    "    df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n",
    "    \n",
    "    return df\n",
    "\n",
    "def read_trainset(filename=\"../input/rsna-intracranial-hemorrhage-detection/stage_1_train.csv\"):\n",
    "    df = pd.read_csv(filename)\n",
    "    df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n",
    "    df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n",
    "    \n",
    "    duplicates_to_remove = [\n",
    "        1598538, 1598539, 1598540, 1598541, 1598542, 1598543,\n",
    "        312468,  312469,  312470,  312471,  312472,  312473,\n",
    "        2708700, 2708701, 2708702, 2708703, 2708704, 2708705,\n",
    "        3032994, 3032995, 3032996, 3032997, 3032998, 3032999\n",
    "    ]\n",
    "    \n",
    "    df = df.drop(index=duplicates_to_remove)\n",
    "    df = df.reset_index(drop=True)\n",
    "    \n",
    "    df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n",
    "    df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n",
    "    \n",
    "    return df\n",
    "\n",
    "    \n",
    "test_df = read_testset()\n",
    "df = read_trainset()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test_df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Train model and predict\n",
    "\n",
    "*Using train, validation and test set* <br>\n",
    "\n",
    "Training for 6 epochs with Adam optimizer, with a learning rate of 0.0005 and decay rate of 0.8. The validation predictions are \\[exponentially weighted\\] averaged over all 6 epochs (same goes for the test set submission later). `fit_and_predict` returns validation and test predictions for all epochs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# train set (80%) and validation set (20%)\n",
    "ss = ShuffleSplit(n_splits=10, test_size=0.1, random_state=42).split(df.index)\n",
    "\n",
    "# lets go for the first fold only\n",
    "train_idx, valid_idx = next(ss)\n",
    "\n",
    "# obtain model\n",
    "model = MyDeepModel(engine=ResNet50, input_dims=(224, 224, 3), batch_size=32, learning_rate=5e-5,\n",
    "                    num_epochs=6, decay_rate=0.8, decay_steps=1, weights=\"imagenet\", verbose=2)\n",
    "\n",
    "# obtain test + validation predictions (history.test_predictions, history.valid_predictions)\n",
    "history = model.fit_and_predict(df.iloc[train_idx], df.iloc[valid_idx], test_df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Submit test predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test_df.iloc[:, :] = np.average(history.test_predictions, axis=0, weights=[2**i for i in range(len(history.test_predictions))])\n",
    "\n",
    "test_df = test_df.stack().reset_index()\n",
    "\n",
    "test_df.insert(loc=0, column='ID', value=test_df['Image'].astype(str) + \"_\" + test_df['Diagnosis'])\n",
    "\n",
    "test_df = test_df.drop([\"Image\", \"Diagnosis\"], axis=1)\n",
    "\n",
    "test_df.to_csv('submission.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Improvements\n",
    "\n",
    "Some improvements that could possibly be made:<br>\n",
    "* Image augmentation (which can be put in `_read()`)\n",
    "* Different learning rate and learning rate schedule\n",
    "* Increased input size\n",
    "* Train longer\n",
    "* Add dense layer and regularization (e.g. `keras.layers.Dropout()` before the output layer)\n",
    "* Adding some optimal windowing\n",
    "<br>\n",
    "<br>\n",
    "*Feel free to comment!*\n"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}