{"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.6"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":14774,"databundleVersionId":875431,"sourceType":"competition"},{"sourceId":187731,"sourceType":"datasetVersion","datasetId":80814}],"dockerImageVersionId":28772,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# About this kernel\n\nIn this kernel, we will explore the complete workflow for the APTOS 2019 competition. We will go through:\n\n1. Loading & Exploration: A quick overview of the dataset\n2. Resize Images: We will resize both the training and test images to 224x224, so that it matches the ImageNet format.\n3. Mixup & Data Generator: We show how to create a data generator that will perform random transformation to our datasets (flip vertically/horizontally, rotation, zooming). This will help our model generalize better to the data, since it is fairly small (only ~3000 images).\n4. Quadratic Weighted Kappa: A thorough overview of the metric used for this competition, with an intuitive example. Check it out!\n5. Model: We will use a DenseNet-121 pre-trained on ImageNet. We will finetune it using Adam for 15 epochs, and evaluate it on an unseen validation set.\n6. Training & Evaluation: We take a look at the change in loss and QWK score through the epochs.\n\n### Unused Methods\n\nThroughout V15-V18 of this kernel, I ablated a few methods that I presented in this kernel. The highest LB score was achieved after I removed:\n* Mixup\n* Optimized Threshold\n\nI decided to keep them in the kernel if it ever becomes useful for you.\n\n### Citations & Resources\n\n* I had the idea of using mixup from [KeepLearning's ResNet50 baseline](https://www.kaggle.com/mathormad/aptos-resnet50-baseline). Since the implementation was in PyTorch, I instead used an [open-sourced keras implementation](https://github.com/yu4u/mixup-generator).\n* The transfer learning procedure is mostly inspired from my [previous kernel for iWildCam](https://www.kaggle.com/xhlulu/densenet-transfer-learning-iwildcam-2019). The workflow was however heavily modified since then.\n* Used similar [method as Abhishek](https://www.kaggle.com/abhishek/optimizer-for-quadratic-weighted-kappa) to find the optimal threshold.\n* [Lex's kernel](https://www.kaggle.com/lextoumbourou/blindness-detection-resnet34-ordinal-targets) prompted me to try using Multilabel instead of multiclass classification, which slightly improved the kappa score.","metadata":{}},{"cell_type":"code","source":"import 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 sklearn.metrics import confusion_matrix, classification_report\n\n%matplotlib inline","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:29.504138Z","iopub.execute_input":"2025-02-01T21:57:29.504491Z","iopub.status.idle":"2025-02-01T21:57:29.513759Z","shell.execute_reply.started":"2025-02-01T21:57:29.504447Z","shell.execute_reply":"2025-02-01T21:57:29.513013Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Set random seed for reproducibility.","metadata":{}},{"cell_type":"code","source":"np.random.seed(2019)\ntf.set_random_seed(2019)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:29.515491Z","iopub.execute_input":"2025-02-01T21:57:29.515911Z","iopub.status.idle":"2025-02-01T21:57:29.558218Z","shell.execute_reply.started":"2025-02-01T21:57:29.515850Z","shell.execute_reply":"2025-02-01T21:57:29.557727Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Loading & Exploration","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv('../input/aptos2019-blindness-detection/train.csv')\ntest_df = pd.read_csv('../input/aptos2019-blindness-detection/test.csv')\nprint(train_df.shape)\nprint(test_df.shape)\ntrain_df.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:29.559095Z","iopub.execute_input":"2025-02-01T21:57:29.559273Z","iopub.status.idle":"2025-02-01T21:57:29.599411Z","shell.execute_reply.started":"2025-02-01T21:57:29.559241Z","shell.execute_reply":"2025-02-01T21:57:29.598727Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_df['diagnosis'].hist()\ntrain_df['diagnosis'].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:29.600287Z","iopub.execute_input":"2025-02-01T21:57:29.600467Z","iopub.status.idle":"2025-02-01T21:57:29.786814Z","shell.execute_reply.started":"2025-02-01T21:57:29.600435Z","shell.execute_reply":"2025-02-01T21:57:29.786036Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Displaying some Sample Images","metadata":{}},{"cell_type":"code","source":"def display_samples(df, columns=4, rows=3):\n    fig=plt.figure(figsize=(5*columns, 4*rows))\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'../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_samples(train_df)","metadata":{"_kg_hide-input":true,"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:29.789999Z","iopub.execute_input":"2025-02-01T21:57:29.790383Z","iopub.status.idle":"2025-02-01T21:57:37.409066Z","shell.execute_reply.started":"2025-02-01T21:57:29.790234Z","shell.execute_reply":"2025-02-01T21:57:37.408051Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Resize Images\n\nWe will resize the images to 224x224, then create a single numpy array to hold the data.","metadata":{}},{"cell_type":"code","source":"def get_pad_width(im, new_shape, is_rgb=True):\n    pad_diff = new_shape - im.shape[0], new_shape - im.shape[1]\n    t, b = math.floor(pad_diff[0]/2), math.ceil(pad_diff[0]/2)\n    l, r = math.floor(pad_diff[1]/2), math.ceil(pad_diff[1]/2)\n    if is_rgb:\n        pad_width = ((t,b), (l,r), (0, 0))\n    else:\n        pad_width = ((t,b), (l,r))\n    return pad_width\n\ndef preprocess_image(image_path, desired_size=224):\n    im = Image.open(image_path)\n    im = im.resize((desired_size, )*2, resample=Image.LANCZOS)\n    \n    return im","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:37.411233Z","iopub.execute_input":"2025-02-01T21:57:37.411484Z","iopub.status.idle":"2025-02-01T21:57:37.418735Z","shell.execute_reply.started":"2025-02-01T21:57:37.411433Z","shell.execute_reply":"2025-02-01T21:57:37.418073Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"N = train_df.shape[0]\nx_train = np.empty((N, 224, 224, 3), dtype=np.uint8)\n\nfor i, image_id in enumerate(tqdm(train_df['id_code'])):\n    x_train[i, :, :, :] = preprocess_image(\n        f'../input/aptos2019-blindness-detection/train_images/{image_id}.png'\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-02-01T21:57:37.419978Z","iopub.execute_input":"2025-02-01T21:57:37.420219Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"N = test_df.shape[0]\nx_test = np.empty((N, 224, 224, 3), dtype=np.uint8)\n\nfor i, image_id in enumerate(tqdm(test_df['id_code'])):\n    x_test[i, :, :, :] = preprocess_image(\n        f'../input/aptos2019-blindness-detection/test_images/{image_id}.png'\n    )","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"y_train = pd.get_dummies(train_df['diagnosis']).values\n\nprint(x_train.shape)\nprint(y_train.shape)\nprint(x_test.shape)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Creating multilabels\n\nInstead of predicting a single label, we will change our target to be a multilabel problem; i.e., if the target is a certain class, then it encompasses all the classes before it. E.g. encoding a class 4 retinopathy would usually be `[0, 0, 0, 1]`, but in our case we will predict `[1, 1, 1, 1]`. For more details, please check out [Lex's kernel](https://www.kaggle.com/lextoumbourou/blindness-detection-resnet34-ordinal-targets).","metadata":{}},{"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))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now we can split it into a training and validation set.","metadata":{}},{"cell_type":"code","source":"x_train, x_val, y_train, y_val = train_test_split(\n    x_train, y_train_multi, \n    test_size=0.15, \n    random_state=2019\n)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Mixup & Data Generator\n\nPlease Note: Although I show how to construct Mixup, **it is currently unused**. Please see notice at the top of the kernel.","metadata":{}},{"cell_type":"code","source":"class MixupGenerator():\n    def __init__(self, X_train, y_train, batch_size=32, alpha=0.2, shuffle=True, datagen=None):\n        self.X_train = X_train\n        self.y_train = y_train\n        self.batch_size = batch_size\n        self.alpha = alpha\n        self.shuffle = shuffle\n        self.sample_num = len(X_train)\n        self.datagen = datagen\n\n    def __call__(self):\n        while True:\n            indexes = self.__get_exploration_order()\n            itr_num = int(len(indexes) // (self.batch_size * 2))\n\n            for i in range(itr_num):\n                batch_ids = indexes[i * self.batch_size * 2:(i + 1) * self.batch_size * 2]\n                X, y = self.__data_generation(batch_ids)\n\n                yield X, y\n\n    def __get_exploration_order(self):\n        indexes = np.arange(self.sample_num)\n\n        if self.shuffle:\n            np.random.shuffle(indexes)\n\n        return indexes\n\n    def __data_generation(self, batch_ids):\n        _, h, w, c = self.X_train.shape\n        l = np.random.beta(self.alpha, self.alpha, self.batch_size)\n        X_l = l.reshape(self.batch_size, 1, 1, 1)\n        y_l = l.reshape(self.batch_size, 1)\n\n        X1 = self.X_train[batch_ids[:self.batch_size]]\n        X2 = self.X_train[batch_ids[self.batch_size:]]\n        X = X1 * X_l + X2 * (1 - X_l)\n\n        if self.datagen:\n            for i in range(self.batch_size):\n                X[i] = self.datagen.random_transform(X[i])\n                X[i] = self.datagen.standardize(X[i])\n\n        if isinstance(self.y_train, list):\n            y = []\n\n            for y_train_ in self.y_train:\n                y1 = y_train_[batch_ids[:self.batch_size]]\n                y2 = y_train_[batch_ids[self.batch_size:]]\n                y.append(y1 * y_l + y2 * (1 - y_l))\n        else:\n            y1 = self.y_train[batch_ids[:self.batch_size]]\n            y2 = self.y_train[batch_ids[self.batch_size:]]\n            y = y1 * y_l + y2 * (1 - y_l)\n\n        return X, y","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"BATCH_SIZE = 32\n\ndef create_datagen():\n    return ImageDataGenerator(\n        zoom_range=0.15,  # set range for random zoom\n        # set mode for filling points outside the input boundaries\n        fill_mode='constant',\n        cval=0.,  # value used for fill_mode = \"constant\"\n        horizontal_flip=True,  # randomly flip images\n        vertical_flip=True,  # randomly flip images\n    )\n\n# Using original generator\ndata_generator = create_datagen().flow(x_train, y_train, batch_size=BATCH_SIZE, seed=2019)\n# Using Mixup\nmixup_generator = MixupGenerator(x_train, y_train, batch_size=BATCH_SIZE, alpha=0.2, datagen=create_datagen())()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Quadratic Weighted Kappa\n\nQuadratic Weighted Kappa (QWK, the greek letter $\\kappa$), 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\n### What is Cohen Kappa?\n\nAccording to the [wikipedia article](https://en.wikipedia.org/wiki/Cohen%27s_kappa), we have\n> The definition of $\\kappa$ is:\n> $$\\kappa \\equiv \\frac{p_o - p_e}{1 - p_e}$$\n> where $p_o$ is the relative observed agreement among raters (identical to accuracy), and $p_e$ is the hypothetical probability of chance agreement, using the observed data to calculate the probabilities of each observer randomly seeing each category.\n\n### How is it computed?\n\nLet's take the example of a binary classification problem. Say we have:","metadata":{}},{"cell_type":"code","source":"true_labels = np.array([1, 0, 1, 1, 0, 1])\npred_labels = np.array([1, 0, 0, 0, 0, 1])","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We can construct the following table:\n\n| true | pred | agreement      |\n|------|------|----------------|\n| 1    | 1    | true positive  |\n| 0    | 0    | true negative  |\n| 1    | 0    | false negative |\n| 1    | 0    | false negative |\n| 0    | 0    | true negative  |\n| 1    | 1    | true positive  |\n\n\nThen the \"observed proportionate agreement\" is calculated exactly the same way as accuracy:\n\n$$\np_o = acc = \\frac{tp + tn}{all} = {2 + 2}{6} = 0.66\n$$\n\nThis can be confirmed using scikit-learn:","metadata":{}},{"cell_type":"code","source":"accuracy_score(true_labels, pred_labels)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Additionally, we also need to compute `p_e`:\n\n$$p_{yes} = \\frac{tp + fp}{all} \\frac{tp + fn}{all} = \\frac{2}{6} \\frac{4}{6} = 0.222$$\n\n$$p_{no} = \\frac{fn + tn}{all} \\frac{fp + tn}{all} = \\frac{4}{6} \\frac{2}{6} = 0.222$$\n\n$$p_{e} = p_{yes} + p_{no} = 0.222 + 0.222 = 0.444$$\n\nFinally,\n\n$$\n\\kappa = \\frac{p_o - p_e}{1-p_e} = \\frac{0.666 - 0.444}{1 - 0.444} = 0.4\n$$\n\nLet's verify with scikit-learn:","metadata":{}},{"cell_type":"code","source":"cohen_kappa_score(true_labels, pred_labels)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### What is the weighted kappa?\n\nThe wikipedia page offer a very concise explanation: \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":{}},{"cell_type":"markdown","source":"### Creating keras callback for QWK","metadata":{}},{"cell_type":"code","source":"class Metrics(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\"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","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model: DenseNet-121","metadata":{}},{"cell_type":"code","source":"densenet = DenseNet121(\n    weights='../input/densenet-keras/DenseNet-BC-121-32-no-top.h5',\n    include_top=False,\n    input_shape=(224,224,3)\n)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def build_model():\n    model = Sequential()\n    model.add(densenet)\n    model.add(layers.GlobalAveragePooling2D())\n    model.add(layers.Dropout(0.5))\n    model.add(layers.Dense(5, activation='sigmoid'))\n    \n    model.compile(\n        loss='binary_crossentropy',\n        optimizer=Adam(lr=0.00005),\n        metrics=['accuracy']\n    )\n    \n    return model","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model = build_model()\nmodel.summary()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Training & Evaluation","metadata":{}},{"cell_type":"code","source":"kappa_metrics = Metrics()\n\nhistory = model.fit_generator(\n    data_generator,\n    steps_per_epoch=x_train.shape[0] / BATCH_SIZE,\n    epochs=15,\n    validation_data=(x_val, y_val),\n    callbacks=[kappa_metrics]\n)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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[['loss', 'val_loss']].plot()\nhistory_df[['acc', 'val_acc']].plot()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.plot(kappa_metrics.val_kappas)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Find best threshold\n\nPlease Note: Although I show how to construct a threshold optimizer, **it is currently unused**. Please see notice at the top of the kernel.","metadata":{}},{"cell_type":"code","source":"model.load_weights('model.h5')\ny_val_pred = model.predict(x_val)\n\ndef compute_score_inv(threshold):\n    y1 = y_val_pred > threshold\n    y1 = y1.astype(int).sum(axis=1) - 1\n    y2 = y_val.sum(axis=1) - 1\n    score = cohen_kappa_score(y1, y2, weights='quadratic')\n    \n    return 1 - score\n\nsimplex = scipy.optimize.minimize(\n    compute_score_inv, 0.5, method='nelder-mead'\n)\n\nbest_threshold = simplex['x'][0]","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load the best saved model (if not already loaded)\nmodel.load_weights('model.h5') \n\n# Get predictions for validation set\ny_pred_probs = model.predict(x_val)  # Predictions as probabilities\ny_pred_classes = (y_pred_probs > 0.5).astype(int) # Convert to classes (0 or 1)\ny_true_classes = y_val.argmax(axis=1) # Extract true class labels","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cm = confusion_matrix(y_true_classes, y_pred_classes.argmax(axis=1)) \n\nprint(\"Confusion Matrix:\")\nprint(cm)\n\n# Optional: Display the confusion matrix as a heatmap (visually appealing)\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nsns.heatmap(cm, annot=True, fmt='d', cmap='Blues') \nplt.xlabel('Predicted')\nplt.ylabel('Actual')\nplt.title('Confusion Matrix')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"report = classification_report(y_true_classes, y_pred_classes.argmax(axis=1)) \n\nprint(\"\\nClassification Report:\")\nprint(report)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Submit","metadata":{}},{"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)","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}