{"cells":[{"metadata":{},"cell_type":"markdown","source":"# APTOS 2019 Blindness Detection\n## Detect diabetic retinopathy to stop blindness before it's too late\n\n### Submission: Aaron Balson\n### Date: August 2019\n\n\nProject Design:\n1.\tSetup ENV and download/import required packages.\n2.\tLoad given datasets to understand its format.\n3.\tExploratory data analysis on csv and images.\n4.\tPerform pre-processing, resizing/rescaling to match that of ImageNet (224x224) as we intend to use DenseNet (which was trained on ImageNet) for transfer learning.\n5.\tUse Multilabel instead of multiclass encoding. Eg: represent class 4 as [1,1,1,1,0] instead of traditional [0,0,0,1,0]\n6.\tPerform data augmentation using Data Generator since training set is small.\n7.\tAdopt DenseNet-121 Architecture and append our dense and output layers at end their conv layers.\n8.\tTrain and evaluate model based on accuracy and kappa score.\n9.  Plot training vs validation (accuracy & loss) to inspect overfitting.\n    (Re-run the training with hyperparameters tuning and optimizers based on accuracy, kappa score (if needed))\n10.\tTest the model and submit to Kaggle for public score.\n"},{"metadata":{},"cell_type":"markdown","source":"## **1. Setup ENV, download/import dependencies**"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n\nimport json\nimport math\nimport os\n\nimport cv2\nfrom PIL import Image\nimport numpy as np\nfrom keras import layers\nfrom keras.applications import DenseNet121\nfrom keras.callbacks import Callback, ModelCheckpoint\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.models import Sequential\nfrom keras.optimizers import Adam\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import cohen_kappa_score, accuracy_score\nimport scipy\nimport tensorflow as tf\nfrom tqdm import tqdm\nfrom keras.applications.densenet import DenseNet121\nimport seaborn as sns\nsns.set()\n\n\nfrom IPython.display import display\n\n\n%matplotlib inline\n\nEPOCHS = 50\nBATCH_SIZE = 16\nSEED = 20031976\nLRATE = 0.00005\nVERBOSE=0\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n#import os\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"## **2. Load Datasets**"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(SEED)\ntf.set_random_seed(SEED)\n\ntrain_df = pd.read_csv('/kaggle/input/aptos2019-blindness-detection/train.csv')\ntest_df = pd.read_csv('/kaggle/input/aptos2019-blindness-detection/test.csv')\nprint(\"Datasets loaded..\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **3. Exploratory Data Analysis**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Data\ndisplay(train_df.head(2))\ndisplay(test_df.head(2))\n\n# Shape of data\nprint(\"train_df shape = \",train_df.shape)\nprint(\"test_df shape = \",test_df.shape)\n\n# Distribution of data\ndisplay(train_df['diagnosis'].value_counts())\nsns.countplot(train_df['diagnosis'], color='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_image(df, rows, columns):\n    fig=plt.figure(figsize=(10, 10))\n\n    for i in range(columns*rows):\n        image_path = df.loc[i,'id_code']\n        image_id = df.loc[i,'diagnosis']\n        img = cv2.imread(f'/kaggle/input/aptos2019-blindness-detection/train_images/{image_path}.png')\n        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n        \n        fig.add_subplot(rows, columns, i+1)\n        plt.title(image_id)\n        plt.imshow(img)\n    \n    plt.tight_layout()\n\ndisplay_image(train_df, 4, 4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **4. Perform Pre-processing, Resampling**"},{"metadata":{},"cell_type":"markdown","source":"**Observation** The dataset is highly imbalanced, with many samples for level 0, and very little for level 3,4. It also has very less samples (3662) compared to other similar retinopathy competitions (~30k). We will use under sampling and over sampling to balance the dataset with each class having 1000 samples.\n\nImages are to be resized/rescaled to match that of ImageNet (224x224) as we intend to use ResNet50 or DenseNet50 (which were trained on ImageNet) for transfer learning."},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_image(image_path, desired_size=224):\n    im = Image.open(image_path)\n    im = im.resize((desired_size, )*2, resample=Image.LANCZOS)\n    return im\n\n# Trail-1 Under sampling by deleting oversized classes (Class:0)\ndef under_sample_make_all_same(df, categories, max_per_category):\n    df = pd.concat([df[df['diagnosis'] == c][:max_per_category] for c in categories])\n    df = df.sample(n=(max_per_category)*len(categories), replace=False, random_state=20031976)\n    df.index = np.arange(len(df))\n    return df\n#train_df = under_sample_make_all_same(train_df,[0,1,2,3,4], 193 ) \n#Under-sample class-0 (1805-805=1000) and Over-sample other classes so each class has 1000 entries\ntrain_df = train_df.drop(train_df[train_df['diagnosis'] == 0].sample(n=805, replace=False).index)\n\nN = train_df.shape[0]\nx_train = np.empty((N, 224, 224, 3), dtype=np.uint8)\n#tqdm\nfor i, image_id in enumerate((train_df['id_code'])):\n    x_train[i, :, :, :] = preprocess_image(\n        f'/kaggle/input/aptos2019-blindness-detection/train_images/{image_id}.png'\n    )\n    \nN = test_df.shape[0]\nx_test = np.empty((N, 224, 224, 3), dtype=np.uint8)\n\nfor i, image_id in enumerate((test_df['id_code'])):\n    x_test[i, :, :, :] = preprocess_image(\n        f'/kaggle/input/aptos2019-blindness-detection/test_images/{image_id}.png'\n    )\n    \ny_train = pd.get_dummies(train_df['diagnosis']).values\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"x_train.shape=\",x_train.shape)\nprint(\"y_train.shape=\",y_train.shape)\nprint(\"x_test.shape=\",x_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Trail-2 Over sampling by increasing undersized classes\nfrom imblearn.over_sampling import SMOTE, ADASYN\nx_resampled, y_resampled = SMOTE(random_state=SEED).fit_sample(x_train.reshape(x_train.shape[0], -1), train_df['diagnosis'].ravel())\n\nprint(\"x_resampled.shape=\",x_resampled.shape)\nprint(\"y_resampled.shape=\",y_resampled.shape)\n\nx_train = x_resampled.reshape(x_resampled.shape[0], 224, 224, 3)\ny_train = pd.get_dummies(y_resampled).values\n\n# Trail-3 No sampling\n\n# Each class should have 1000 samples now (5 x 1000 = 5000)\nprint(\"x_train.shape=\",x_train.shape)\nprint(\"y_train.shape=\",y_train.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5.Creating Multilabels\nEncoding a class by encompassing all the classes before it seems to work out well for kappa scores. Instead of predicting a single label, we will change our target to be a multilabel problem; i.e., encoding a class 4 retinopathy would usually be [0, 0, 0, 1, 0], but in our case we will predict [1, 1, 1, 1, 0]. For more details, please check out this paper on [Ordinal Regression](https://arxiv.org/abs/0704.1028), also implemented by [Lex's](https://www.kaggle.com/lextoumbourou/blindness-detection-resnet34-ordinal-targets) kernel."},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train_multi = np.empty(y_train.shape, dtype=y_train.dtype)\ny_train_multi[:, 4] = y_train[:, 4]\n\nfor i in range(3, -1, -1):\n    y_train_multi[:, i] = np.logical_or(y_train[:, i], y_train_multi[:, i+1])\n\nprint(\"Original y_train:\", y_train.sum(axis=0))\nprint(\"Multilabel version:\", y_train_multi.sum(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Split 85-15 training-validation sets\nx_sptrain, x_spval, y_sptrain, y_spval = train_test_split(\n    x_train, y_train_multi, \n    test_size=0.10, \n    random_state=SEED\n)\nprint(\"train-validation splitted ...\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 6. Image Data Augmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_datagen():\n    return ImageDataGenerator(\n        zoom_range=0.10,        # set range for random zoom\n        fill_mode='constant',   # set mode for filling points outside the input boundaries\n        cval=0.,                # value used for fill_mode = \"constant\"\n        horizontal_flip=True,   # randomly flip images\n        vertical_flip=True,     # randomly flip images\n        #rotation_range=20       # Degree range for random rotations\n    )\n\n# Using original generator\ndata_generator = create_datagen().flow(x_sptrain, y_sptrain, batch_size=BATCH_SIZE, seed=SEED)\nprint(\"Image data augmentated ...\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 7. CNNetwork Architecture (DenseNet Transfer Learning)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define evaluation metrics\n\nimport keras.backend as K\n\ndef precision(y_true, y_pred):\n    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))\n    predicted_positives = K.sum(K.round(K.clip(y_pred, 0, 1)))\n    precision = true_positives / (predicted_positives + K.epsilon())\n    return precision\n\ndef recall(y_true, y_pred):\n    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))\n    possible_positives = K.sum(K.round(K.clip(y_true, 0, 1)))\n    recall = true_positives / (possible_positives + K.epsilon())\n    return recall\n\ndef fbeta_score(y_true, y_pred, beta=1):\n    if beta < 0:\n        raise ValueError('The lowest choosable beta is zero (only precision).')\n\n    # If there are no true positives, fix the F score at 0 like sklearn.\n    if K.sum(K.round(K.clip(y_true, 0, 1))) == 0:\n        return 0\n\n    p = precision(y_true, y_pred)\n    r = recall(y_true, y_pred)\n    bb = beta ** 2\n    fbeta_score = (1 + bb) * (p * r) / (bb * p + r + K.epsilon())\n    return fbeta_score\n\ndef fmeasure(y_true, y_pred):\n    return fbeta_score(y_true, y_pred, beta=1)\n\ndef mean_pred(y_true, y_pred):\n    return K.mean(y_pred)\n\ndef f1_score(y_true, y_pred):\n    p = precision(y_true, y_pred)\n    r = recall(y_true, y_pred)\n    return 2*(p*r) / (p+r+K.epsilon())\n\nprint(\"Evaluation metrics defined ...\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Transfer Learning\ndensenet = DenseNet121(\n    weights='/kaggle/input/densenet121/DenseNet-BC-121-32-no-top.h5',\n    include_top=False,\n    input_shape=(224,224,3)\n)\n\nmodel = Sequential()\nmodel.add(densenet)\nmodel.add(layers.GlobalAveragePooling2D())\nmodel.add(layers.Dropout(0.2))\nmodel.add(layers.Dense(5, activation='sigmoid'))\n\nmodel.compile(\n    loss='binary_crossentropy',\n    optimizer=Adam(lr=LRATE),\n    metrics=['accuracy',mean_pred, precision, recall, f1_score, fbeta_score, fmeasure]\n)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 8. Training Model"},{"metadata":{},"cell_type":"markdown","source":"** Quadratic Weighted Kappa **\nQuadratic Weighted Kappa (QWK, the greek letter  κ ), also known as Cohen's Kappa, is the official evaluation metric. For our kernel, we will use a custom callback to monitor the score, and plot it at the end.\n\nWhat is the weighted kappa? The wikipedia page offer a very concise explanation:\n\n> The weighted kappa allows disagreements to be weighted differently and is especially useful when codes are ordered. Three matrices are involved, the matrix of observed scores, the matrix of expected scores based on chance agreement, and the weight matrix. Weight matrix cells located on the diagonal (upper-left to bottom-right) represent agreement and thus contain zeros. Off-diagonal cells contain weights indicating the seriousness of that disagreement.\n> \nSimply put, if two scores disagree, then the penalty will depend on how far they are apart. That means that our score will be higher if (a) the real value is 4 but the model predicts a 3, and the score will be lower if (b) the model instead predicts a 0. This metric makes sense for this competition, since the labels 0-4 indicates how severe the illness is. Intuitively, a model that predicts a severe retinopathy (3) when it is in reality a proliferative retinopathy (4) is probably better than a model that predicts a mild retinopathy (1)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# callback to keep track of kappa score during training\nclass KappaMetrics(Callback):\n    def on_train_begin(self, logs={}):\n        self.val_kappas = []\n        \n    def on_epoch_end(self, epoch, logs={}):\n        X_val, y_val = self.validation_data[:2]\n        y_val = y_val.sum(axis=1) - 1\n        \n        y_pred = self.model.predict(X_val) > 0.5\n        y_pred = y_pred.astype(int).sum(axis=1) - 1\n\n        _val_kappa = cohen_kappa_score(\n            y_val,\n            y_pred, \n            weights='quadratic'\n        )\n\n        self.val_kappas.append(_val_kappa)\n\n        print(f\"Epoch: {epoch+1} val_kappa: {_val_kappa:.4f}\")\n        \n        if _val_kappa == max(self.val_kappas):\n            print(\"Validation Kappa has improved. Saving model.\")\n            self.model.save('model.h5')\n\n        return\n    \nkappa_score = KappaMetrics()\n\n\n\nhistory = model.fit_generator(\n    data_generator,\n    steps_per_epoch=x_train.shape[0] / BATCH_SIZE,\n    epochs=EPOCHS,\n    validation_data=(x_spval, y_spval),\n    callbacks=[kappa_score],\n    verbose=VERBOSE\n)    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 9. Evaluate model using metrics"},{"metadata":{"trusted":true},"cell_type":"code","source":"with open('history.json', 'w') as f:\n    json.dump(history.history, f)\n\nhistory_df = pd.DataFrame(history.history)\nhistory_df.head(EPOCHS)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#history_df[['loss', 'val_loss']].plot()\n#history_df[['acc', 'val_acc']].plot()\n#history_df[['acc', 'precision', 'recall', 'f1_score', 'fbeta_score', 'fmeasure']].plot()\n\nf1, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(24, 4))\nt1 = f1.suptitle('CNN Performance', fontsize=12)\nf1.subplots_adjust(top=0.85, wspace=0.3)\n\nepoch_list = list(range(1,EPOCHS + 1))\nax1.plot(epoch_list, history.history['acc'], label='Train Accuracy')\nax1.plot(epoch_list, history.history['val_acc'], label='Validation Accuracy')\nax1.set_xticks(np.arange(0, EPOCHS + 1, 5))\nax1.set_ylabel('Accuracy %')\nax1.set_xlabel('Epoch')\nax1.set_title('Accuracy')\nl1 = ax1.legend(loc=\"best\")\n\nax2.plot(epoch_list, history.history['loss'], label='Train Loss')\nax2.plot(epoch_list, history.history['val_loss'], label='Validation Loss')\nax2.set_xticks(np.arange(0, EPOCHS + 1, 5))\nax2.set_ylabel('Loss %')\nax2.set_xlabel('Epoch')\nax2.set_title('Loss')\nl2 = ax2.legend(loc=\"best\")\n\nax3.plot(epoch_list, history.history['acc'], label='Accuracy')\nax3.plot(epoch_list, history.history['precision'], label='Precision')\nax3.plot(epoch_list, history.history['recall'], label='Recall')\nax3.plot(epoch_list, history.history['f1_score'], label='F1 score')\nax3.plot(epoch_list, history.history['fbeta_score'], label='Fbeta score')\nax3.plot(epoch_list, history.history['fmeasure'], label='FMeasure')\nax3.set_xticks(np.arange(0, EPOCHS + 1, 5))\nax3.set_ylabel('Score')\nax3.set_xlabel('Epoch')\nax3.set_title('Performance')\nl3 = ax3.legend(loc=\"best\")\n\nax4.plot(epoch_list, kappa_score.val_kappas, label='Kappa score')\nax4.set_xticks(np.arange(0, EPOCHS + 1, 5))\nax4.set_ylabel('Score')\nax4.set_xlabel('Epoch')\nax4.set_title('Kappa Metrics')\nl4 = ax4.legend(loc=\"best\")\n\ndisplay(\"Maximum Kappa Score: %s\" %max(kappa_score.val_kappas))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"** Observation: ** Model has shown good learning capability. It seems to generalize well since training vs validation accuracy are within 4% and loss within 0.1%. We can increase few more epoch to further train to acheive better kappa score without starting to overfit."},{"metadata":{},"cell_type":"markdown","source":"# 10. Predict & Submit for results\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_test = model.predict(x_test) > 0.5\ny_test = y_test.astype(int).sum(axis=1) - 1\n\ntest_df['diagnosis'] = y_test\ntest_df.to_csv('submission.csv',index=False)\ndisplay(test_df.head(5))\n\nimport datetime\nprint(\"Ran at UTC : \", datetime.datetime.utcnow())\n","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}