{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Abstract\n\n**This notebook introduces our approach to the \"RSNA Intracranial Hemorrhage Detection\" Competition. In specific, we focus on pointing out how we tuned the ResNet50 pretrained neural net from the Keras library in order to achieve worthwile predictions of hemorrhages and their subtype.**\n"},{"metadata":{},"cell_type":"markdown","source":"### Acknowledgements\nWe benefited heavily from the public notebooks available for this competition, including amongst others:\n* [Akenesert](https://www.kaggle.com/akensert/inceptionv3-prev-resnet50-keras-baseline-model)'s notebook on the InceptionV3 implementation\n* [Ryan Epp](https://www.kaggle.com/reppic/gradient-sigmoid-windowing)'s notebook on windowing\n* [Marco E](https://www.kaggle.com/marcovasquez/basic-eda-data-visualization)'s EDA notebook\n* ...\n\n"},{"metadata":{},"cell_type":"markdown","source":"### Further Information\n\n* [Residual Networks in Keras 1](https://towardsdatascience.com/hitchhikers-guide-to-residual-networks-resnet-in-keras-385ec01ec8ff)/[Residual Networks in Keras 2](https://towardsdatascience.com/understanding-and-coding-a-resnet-in-keras-446d7ff84d33), and [ResNet50](https://www.kaggle.com/keras/resnet50) in specific \n* [Hyperparameter tuning with Keras Tuner](https://keras-team.github.io/keras-tuner/) (and a [simple introduction](https://www.mikulskibartosz.name/how-to-automaticallyselect-the-hyperparameters-of-a-resnet-neural-network/))"},{"metadata":{},"cell_type":"markdown","source":"# 1. Preparation"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nfrom math import ceil, floor, log\nfrom os import listdir\nfrom os.path import isfile, join\nimport matplotlib.pyplot as plt\n\nimport tensorflow as tf\nimport keras\n#from keras_applications.resnet import ResNet50\nfrom keras_applications.inception_v3 import InceptionV3\nfrom sklearn.model_selection import ShuffleSplit","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can set the directory right to the Kaggle backend. This makes handling the data thoroughout this notebook much easier:"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Set directories\nDATA_DIR = '../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection'\n\nTRAIN_IMAGES_DIR = DATA_DIR + '/stage_2_train/'\nTRAIN_CSV_DIR = DATA_DIR + '/stage_2_train.csv'\n\nTEST_IMAGES_DIR = DATA_DIR + '/stage_2_test/'\nTEST_CSV_DIR = DATA_DIR + '/stage_2_sample_submission.csv'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2. Data Cleaning and Exploration"},{"metadata":{},"cell_type":"markdown","source":"We have 4,416,818 unique training labels but there are duplicates we need to remove. So we define a function to load the train labels correctly and also a function to prepare the test labels."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check for duplicates in trainset csv\ntrain = pd.read_csv(TRAIN_CSV_DIR)\nprint(train.head(10))\nprint(len(train)) # 4,516,842 training labels\nprint(len(train.ID.unique())) # 4,516,818 unique training labels\n\n# Identify duplicate rows in trainset\nduplicates = train[train.duplicated(subset=None, keep='first')]\nprint(duplicates)\n\n# Define function to remove duplicates when loading trainset\ndef read_trainset(filename=TRAIN_CSV_DIR):\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        56346, 56347, 56348, 56349, 56350, 56351,\n        1171830, 1171831, 1171832, 1171833, 1171834,\n        1171835, 3705312, 3705313, 3705314, 3705315, \n        3705316, 3705317, 3842478, 3842479, 3842480,\n        3842481, 3842482, 3842483\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# Define function to load testset\ndef read_testset(filename= TEST_CSV_DIR):\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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then we load the .CSV data correctly for train as well as test."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the data as two dataframes\ndf = read_trainset()\ndf = df.sample(frac = .00001)\ntest_df = read_testset()\ntest_df = test_df.sample(frac = .00001)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(df.tail(7))\nprint(df.shape)  \n# 752,803 train images with 6 labels each = 4,516,818 labels \n# Amount of unique training labels from earlier equals amount of images labeled in csv","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Out of 752,803 train images we have 107,933 labeled as hemorrhage. The types of the hemorrhages identified vary heavily with epidural being very few and subdural being the most subtypes of hemorrhages identified:"},{"metadata":{"trusted":true},"cell_type":"code","source":"label_sum = df.sum()\nprint(label_sum)\nlabel_sum.plot.bar()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Image Preparation"},{"metadata":{},"cell_type":"markdown","source":"Before moving on with loading the images, let's first check how the images look like in the original data: "},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm = pydicom.dcmread(TRAIN_IMAGES_DIR + 'ID_ffffb670a' + '.dcm')\nplt.imshow(dcm.pixel_array, cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Throughout the Kaggle Competition, many challenges around DICOM images have been idetified and discussed. In the following, we use some of these insights as an important pre-processing step:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Shift everything up by 1000, then move the values larger than 2048 back to where they should have been. (JER)\ndef correct_dcm(dcm):\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n\n# Window the DICOM image\ndef window_image(dcm, window_center, window_width):\n    if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n        correct_dcm(dcm)\n    \n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n\n    return img\n\n# Save different images as R, G, B: in 3 dimensions\ndef bsb_window(dcm):\n    brain_img = window_image(dcm, 40, 80)\n    subdural_img = window_image(dcm, 80, 200)\n    soft_img = window_image(dcm, 40, 380)\n    \n    brain_img = (brain_img - 0) / 80\n    subdural_img = (subdural_img - (-20)) / 200\n    soft_img = (soft_img - (-150)) / 380\n    bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1,2,0)\n\n    return bsb_img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After the corrections and the windowing, the picture looks like this:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(bsb_window(dcm), cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We then create a helper function to load and store the images correctly. We also want 3-channel inputs for e.g. the ResNet. To do so, we use the bsb_window functino was defined earlier:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def data_read(path, desired_size):\n    dcm = pydicom.dcmread(path)\n    \n    try:\n        img = bsb_window(dcm)\n    except:\n        img = np.zeros(desired_size)\n    \n    img = cv2.resize(img, desired_size[:2], interpolation=cv2.INTER_LINEAR) # Make image smaller\n    \n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(data_read(TRAIN_IMAGES_DIR +'ID_ffffb670a'+'.dcm', (128, 128)), cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 4. Data Generator"},{"metadata":{},"cell_type":"markdown","source":"We then introduce a so-called Data Generator that enables us to load the data in parts and to feed it respectively into the neural net:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Inherits from keras.utils.Sequence object and thus should be safe for multiprocessing\nclass DataGenerator(keras.utils.Sequence):\n\n    def __init__(self, list_IDs, labels=None, batch_size=1, img_size=(512, 512, 1), \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(np.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.35, 1: 0.5})\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,] = data_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,] = data_read(self.img_dir+ID+\".dcm\", self.img_size)\n            \n            return X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"But before using the Data Generator we want to quickly check the amount of images: Does it match the amount of labels?"},{"metadata":{"trusted":true},"cell_type":"code","source":"#train_check = [f for f in listdir(TRAIN_IMAGES_DIR) if isfile(join(TRAIN_IMAGES_DIR, f))]\n#test_check = [f for f in listdir(TEST_IMAGES_DIR) if isfile(join(TEST_IMAGES_DIR, f))]\n\n#print('Amount train images:', len(train_check))\nprint('Amount train labels:', df.shape)\n\n#print('Amount test images:', len(test_check))\nprint('Amount test labels:', test_df.shape) \n\n# Looks good!","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5. Loss function"},{"metadata":{},"cell_type":"markdown","source":"We are ready to move on. Before continuing with the model, we need to specify the details for the evaluation of the predictive power:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras import backend as K\n\ndef 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\ndef _normalized_weighted_average(arr, weights=None):\n    \"\"\"\n    A simple Keras implementation that mimics that of \n    numpy.average(), specifically for 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\ndef 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\ndef 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()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 6. Model"},{"metadata":{},"cell_type":"markdown","source":"Model is divided into the following parts: \n\nThe input image is then passed through InceptionV3 (\"engine\"). InceptionV3 could be replaced by any of the available architectures in keras_application.\n\nFinally, the output from InceptionV3 goes through average pooling followed by two dense layers (including output layer)."},{"metadata":{"trusted":true},"cell_type":"code","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        # Commented out to save time\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 also save the predictions with np.save()\n\n\nclass 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,\n                             backend = keras.backend, layers = keras.layers,\n                             models = keras.models, utils = keras.utils)\n        \n        x = keras.layers.GlobalAveragePooling2D(name='avg_pool')(engine.output)\n#         x = keras.layers.Dropout(0.2)(x)\n#         x = keras.layers.Dense(keras.backend.int_shape(x)[1], activation=\"relu\", name=\"dense_hidden_1\")(x)\n#         x = keras.layers.Dropout(0.1)(x)\n        out = keras.layers.Dense(1, activation=\"sigmoid\", name='dense_output')(x)\n\n        self.model = keras.models.Model(inputs=engine.input, outputs=out)\n\n        self.model.compile(loss=\"binary_crossentropy\", optimizer=keras.optimizers.Adam(), 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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 7. Training and Prediction"},{"metadata":{},"cell_type":"markdown","source":"5. Train model and predict\n\nUsing train, validation and test set \n\nTraining for 5 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 5 epochs (not in this commit). fit_and_predict returns validation and test predictions for all epochs."},{"metadata":{"trusted":true},"cell_type":"code","source":"# train set (00%) and validation set (10%)\nss = ShuffleSplit(n_splits=10, test_size=0.2, random_state=42).split(df.index)\n\n# lets go for the first fold only\ntrain_idx, valid_idx = next(ss)\n\n# obtain model\nmodel = MyDeepModel(engine=InceptionV3, input_dims=(256, 256, 3), batch_size=32, learning_rate=5e-4,               \n                    num_epochs=2, decay_rate=0.8, decay_steps=1, weights=\"imagenet\", verbose=1)\n\n# obtain test + validation predictions (history.test_predictions, history.valid_predictions)\nhistory = model.fit_and_predict(df.iloc[train_idx], df.iloc[valid_idx], test_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 8. Hyperparameter Tuning"},{"metadata":{},"cell_type":"markdown","source":"# 9. Evaluation"},{"metadata":{},"cell_type":"markdown","source":"# 10. Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.iloc[:, :] = np.average(history.test_predictions, axis=0, weights=[0, 1, 2, 4, 6]) # let's do a weighted average for epochs (>1)\n\ntest_df = test_df.stack().reset_index()\n\ntest_df.insert(loc=0, column='ID', value=test_df['Image'].astype(str) + \"_\" + test_df['Diagnosis'])\n\ntest_df = test_df.drop([\"Image\", \"Diagnosis\"], axis=1)\n\ntest_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}