{"cells":[{"metadata":{},"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 if it ever becomes useful in your case.\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_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"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 import callbacks\nfrom keras.applications import DenseNet121\nfrom keras.applications.xception import Xception\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\nfrom tqdm import tqdm\nimport imgaug as ia\nfrom imgaug import augmenters as iaa\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAINING = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Loading & Exploration"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = pd.read_csv('../input/aptos2019-blindness-detection/test.csv')\nprint(test_df.shape)\nif TRAINING:\n    train_df = pd.read_csv('../input/aptos2019-blindness-detection/train.csv')\n    print(train_df.shape)\n    train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # test_df_19 = pd.read_csv('../input/aptos2019-blindness-detection/test.csv')\n# # print(test_df_19.shape)\n# if TRAINING:\n#     train_df_15 = pd.read_csv('../input/resized-2015-2019-blindness-detection-images/')\n#     print(train_df_19.shape)\n#     train_df_19.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!head -n5 ../input/resized-2015-2019-blindness-detection-images/labels/trainLabels15.csv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    train_df['diagnosis'].hist()\n    train_df['diagnosis'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Displaying some Sample Images"},{"metadata":{"_kg_hide-input":true,"trusted":true},"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    \nif TRAINING:\n    display_samples(train_df)","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"IMG_SIZE = 224\nNB_CHANNELS = 3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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 standardize(x):\n    x = x.astype(np.float32)\n    x = x / np.max(x)\n    return (x - np.mean(x)) / (np.std(x))\n\n\"\"\"\nnormalizes the image\n\"\"\"\ndef normalize(img):\n    img = ((img - np.min(img)) / (np.max(img) - np.min(img))) * 255\n    return img.astype(np.uint8)\n#     return img\n\n\"\"\"\ncrops black parts around the image (intensity is <= tol)\n\"\"\"\ndef crop_image(img, tol=10):\n    \n    # for one channel\n    def crop_image_1(img):\n        mask = img > tol\n        return img[np.ix_(mask.any(1), mask.any(0))]\n    \n    if img.ndim == 2:\n        return crop_image_1(img)\n    \n    elif img.ndim == 3:\n        try:\n            img_cpy = img.copy()\n            h, w, _ = img.shape\n            img1 = cv2.resize(crop_image_1(img[:, :, 0]), (w, h))\n            img2 = cv2.resize(crop_image_1(img[:, :, 1]), (w, h))\n            img3 = cv2.resize(crop_image_1(img[:, :, 2]), (w, h))\n\n            img[:,:,0] = img1\n            img[:,:,1] = img2\n            img[:,:,2] = img3\n            \n        except:\n            return img_cpy\n        \n        \n        \n        # add edges to img\n        \n#         sobelx = cv2.Sobel(img,cv2.CV_32F,1,0,ksize=5)\n#         sobely = cv2.Sobel(img,cv2.CV_32F,0,1,ksize=5)\n#         sobelx = cv2.cvtColor(sobelx, cv2.COLOR_RGB2GRAY)\n#         sobely = cv2.cvtColor(sobely, cv2.COLOR_RGB2GRAY)\n        \n#         img4 = sobelx ** 2 + sobely ** 2\n        \n#         img4 -= np.min(img4)\n#         img4 = img4 / np.max(img4)\n        \n#         img4 *= 255.\n#         img4.astype(np.uint8)\n        \n#         img[:,:,0] = (9.0 * img1 + img4) // 10\n#         img[:,:,1] = (9.0 * img2 + img4) // 10\n#         img[:,:,2] = (9.0 * img3 + img4) // 10\n        \n        \n        return img\n\n\n\"\"\"\ncrops black parts and enhances image (Ben Graham's method)\n\"\"\"\n# def preprocess_image(img):\n#     # Added this line\n#     img = normalize(img)\n    \n#     img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n#     img = crop_image(img)\n#     img = cv2.resize(img, (IMG_SIZE, IMG_SIZE))\n#     img = cv2.addWeighted(img, 4, cv2.GaussianBlur(img, (0, 0), IMG_SIZE/10), -4, 128)\n    \n#     return img.astype(np.float32)\n\n\ndef preprocess_image(im):\n#     im = Image.open(image_path)\n#     im = cv2.imread(image_path)\n    im = cv2.cvtColor(im, cv2.COLOR_BGR2RGB)\n    \n    im = normalize(im)\n    im = crop_image(im)\n#     im = im.resize((IMG_SIZE, ) * 2, resample=Image.LANCZOS)\n    im = cv2.resize(im, (IMG_SIZE, IMG_SIZE))\n    \n    im = cv2.addWeighted(im, 4, cv2.GaussianBlur(im, (0, 0), IMG_SIZE / 10), -4, 128)\n    \n#     im = standardize(im)\n#     im = normalize(im)\n    \n    return im.astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    fig, ax = plt.subplots(1, 4, figsize=[20, 5])\n    for i, image_id in enumerate(train_df['id_code']):\n    #     img = preprocess_image(f'../input/aptos2019-blindness-detection/train_images/{image_id}.png')\n        im = preprocess_image(cv2.imread(f'../input/aptos2019-blindness-detection/train_images/{image_id}.png'))\n#         im = crop_image(im)\n#         im = cv2.resize(im, (IMG_SIZE, IMG_SIZE))\n        ax[i].imshow(im)\n        ax[i].axis(\"off\")\n        if i == 3:\n            break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    N = train_df.shape[0]\n    x_train = np.empty((N, 224, 224, 3), dtype=np.uint8)\n\n    for i, image_id in enumerate(tqdm(train_df['id_code'])):\n        x_train[i, :, :, :] = preprocess_image(cv2.imread(\n            f'../input/aptos2019-blindness-detection/train_images/{image_id}.png'\n        ))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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(cv2.imread(\n        f'../input/aptos2019-blindness-detection/test_images/{image_id}.png'\n    ))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    y_train = pd.get_dummies(train_df['diagnosis']).values\n\n    print(x_train.shape)\n    print(y_train.shape)\nprint(x_test.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can split it into a training and validation set."},{"metadata":{"trusted":true},"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# )\n\nif TRAINING:\n    x_train, x_val, y_train, y_val = train_test_split(\n#         x_train, y_train_multi, \n        x_train, y_train, \n        test_size=0.2\n    )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Oversampling"},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     x_train_0 = x_train[[y_train == 1][0][:, 0]]\n#     x_train_1 = x_train[[y_train == 1][0][:, 1]]\n#     x_train_2 = x_train[[y_train == 1][0][:, 2]]\n#     x_train_3 = x_train[[y_train == 1][0][:, 3]]\n#     x_train_4 = x_train[[y_train == 1][0][:, 4]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# def oversample(x):\n#     indexes_asc = np.arange(len(x))\n#     indexes_ran = np.random.permutation(len(x))\n    \n#     return (0.5 * x[indexes_asc] + 0.5 * x[indexes_ran]).clip(0, 255).astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     x_train_1 = np.concatenate([x_train_1, oversample(x_train_1)])\n#     x_train_1 = np.concatenate([x_train_1, oversample(x_train_1)])\n\n#     x_train_2 = np.concatenate([x_train_2, oversample(x_train_2)])\n\n#     x_train_3 = np.concatenate([x_train_3, oversample(x_train_3)])\n#     x_train_3 = np.concatenate([x_train_3, oversample(x_train_3)])\n#     x_train_3 = np.concatenate([x_train_3, oversample(x_train_3)])\n\n#     x_train_4 = np.concatenate([x_train_4, oversample(x_train_4)])\n#     x_train_4 = np.concatenate([x_train_4, oversample(x_train_4)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     print(x_train_0.shape)\n#     print(x_train_1.shape)\n#     print(x_train_2.shape)\n#     print(x_train_3.shape)\n#     print(x_train_4.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     y_train_0 = np.zeros((len(x_train_0), 5), dtype=np.bool)\n#     y_train_0[:, 0] = True\n\n#     y_train_1 = np.zeros((len(x_train_1), 5), dtype=np.bool)\n#     y_train_1[:, 1] = True\n\n#     y_train_2 = np.zeros((len(x_train_2), 5), dtype=np.bool)\n#     y_train_2[:, 2] = True\n\n#     y_train_3 = np.zeros((len(x_train_3), 5), dtype=np.bool)\n#     y_train_3[:, 3] = True\n\n#     y_train_4 = np.zeros((len(x_train_4), 5), dtype=np.bool)\n#     y_train_4[:, 4] = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     x_train = np.concatenate([x_train_0, x_train_1, x_train_2, x_train_3, x_train_4])\n#     y_train = np.concatenate([y_train_0, y_train_1, y_train_2, y_train_3, y_train_4]).astype(np.uint8)\n\n#     del x_train_0\n#     del x_train_1\n#     del x_train_2\n#     del x_train_3\n#     del x_train_4\n\n#     del y_train_0\n#     del y_train_1\n#     del y_train_2\n#     del y_train_3\n#     del y_train_4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    indexes = np.random.permutation(len(x_train))\n\n    x_train = x_train[indexes]\n    y_train = y_train[indexes]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     print(x_train.dtype, x_train.shape)\n#     print(y_train.dtype, y_train.shape)","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n#     y_train_multi = np.empty(y_train.shape, dtype=y_train.dtype)\n#     y_train_multi[:, 4] = y_train[:, 4]\n\n    for i in range(3, -1, -1):\n#         y_train_multi[:, i] = np.logical_or(y_train[:, i], y_train_multi[:, i+1])\n        y_train[:, i] = np.logical_or(y_train[:, i], y_train[:, i + 1])\n        y_val[:, i] = np.logical_or(y_val[:, i], y_val[:, i + 1])\n\n#     print(\"Original y_train:\", y_train.sum(axis=0))\n#     print(\"Multilabel version:\", y_train_multi.sum(axis=0))\n    print(\"Multilabel version:\", y_train.sum(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    print(x_train.shape, x_train.dtype, np.min(x_train), np.max(x_train))\n    print(x_val.shape, x_val.dtype, np.min(x_val), np.max(x_val))\n    \nprint(x_test.shape, x_test.dtype, np.min(x_test), np.max(x_test))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Crop training and validation data as test data\nif TRAINING:\n    for i in tqdm(range(len(x_val))):\n        x_val[i] = cv2.resize(x_val[i, 20: -20, 20: -20, :], (IMG_SIZE, IMG_SIZE))\n        \n    for i in tqdm(range(len(x_train))):\n        x_train[i] = cv2.resize(x_train[i, 20: -20, 20: -20, :], (IMG_SIZE, IMG_SIZE))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    fig, ax = plt.subplots(5, 5, figsize=(10, 10))\n    indexes = np.random.randint(len(x_train), size=25)\n    for i in range(5):\n        for j in range(5):\n            ax[i][j].imshow(x_train[indexes[i * 5 + j]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    fig, ax = plt.subplots(5, 5, figsize=(10, 10))\n    indexes = np.random.randint(len(x_val), size=25)\n    for i in range(5):\n        for j in range(5):\n            ax[i][j].imshow(x_val[indexes[i * 5 + j]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(5, 5, figsize=(10, 10))\nindexes = np.random.randint(len(x_test), size=25)\nfor i in range(5):\n    for j in range(5):\n        ax[i][j].imshow(x_test[indexes[i * 5 + j]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    print(x_train.shape, x_train.dtype, np.min(x_train), np.max(x_train))\n    print(x_val.shape, x_val.dtype, np.min(x_val), np.max(x_val))\n    \nprint(x_test.shape, x_test.dtype, np.min(x_test), np.max(x_test))","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"class MixupGenerator():\n    def __init__(self, X_train, y_train, seq=None, 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        self.seq = seq\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 self.seq:\n            X = X.astype(np.float32)\n            X /= np.max(X)\n            X *= 255 \n            X = X.astype(np.uint8)\n            X = self.seq.augment_images(X)\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\n    \n    def get_sample(self):\n        return self.__data_generation(self.__get_exploration_order()[i * self.batch_size * 2:(i + 1) * self.batch_size * 2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sometimes = lambda aug: iaa.Sometimes(0.5, aug)\nseq = iaa.Sequential(\n        [\n            # apply the following augmenters to most images\n            iaa.Fliplr(0.5), # horizontally flip 50% of all images\n            iaa.Flipud(0.2), # vertically flip 20% of all images\n            sometimes(iaa.Affine(\n                scale={\"x\": (0.9, 1.1), \"y\": (0.9, 1.1)}, # scale images to 80-120% of their size, individually per axis\n                translate_percent={\"x\": (-0.1, 0.1), \"y\": (-0.1, 0.1)}, # translate by -20 to +20 percent (per axis)\n                rotate=(-10, 10), # rotate by -45 to +45 degrees\n                shear=(-5, 5), # shear by -16 to +16 degrees\n                order=[0, 1], # use nearest neighbour or bilinear interpolation (fast)\n                cval=(0, 255), # if mode is constant, use a cval between 0 and 255\n                mode=ia.ALL # use any of scikit-image's warping modes (see 2nd image from the top for examples)\n            )),\n            # execute 0 to 5 of the following (less important) augmenters per image\n            # don't execute all of them, as that would often be way too strong\n            iaa.SomeOf((0, 5),\n                [\n                    sometimes(iaa.Superpixels(p_replace=(0, 1.0), n_segments=(20, 200))), # convert images into their superpixel representation\n                    iaa.OneOf([\n                        iaa.GaussianBlur((0, 1.0)), # blur images with a sigma between 0 and 3.0\n                        iaa.AverageBlur(k=(3, 5)), # blur image using local means with kernel sizes between 2 and 7\n                        iaa.MedianBlur(k=(3, 5)), # blur image using local medians with kernel sizes between 2 and 7\n                    ]),\n                    iaa.Sharpen(alpha=(0, 1.0), lightness=(0.9, 1.1)), # sharpen images\n                    iaa.Emboss(alpha=(0, 1.0), strength=(0, 2.0)), # emboss images\n                    # search either for all edges or for directed edges,\n                    # blend the result with the original image using a blobby mask\n                    iaa.SimplexNoiseAlpha(iaa.OneOf([\n                        iaa.EdgeDetect(alpha=(0.5, 1.0)),\n                        iaa.DirectedEdgeDetect(alpha=(0.5, 1.0), direction=(0.0, 1.0)),\n                    ])),\n                    iaa.AdditiveGaussianNoise(loc=0, scale=(0.0, 0.01*255), per_channel=0.5), # add gaussian noise to images\n                    iaa.OneOf([\n                        iaa.Dropout((0.01, 0.05), per_channel=0.2), # randomly remove up to 10% of the pixels\n                        iaa.CoarseDropout((0.01, 0.03), size_percent=(0.01, 0.02), per_channel=0.2),\n                    ]),\n                    iaa.Invert(0.01, per_channel=True), # invert color channels\n                    iaa.Add((-2, 2), per_channel=0.5), # change brightness of images (by -10 to 10 of original value)\n                    iaa.AddToHueAndSaturation((-1, 1)), # change hue and saturation\n                    # either change the brightness of the whole image (sometimes\n                    # per channel) or change the brightness of subareas\n                    iaa.OneOf([\n                        iaa.Multiply((0.9, 1.1), per_channel=0.5),\n                        iaa.FrequencyNoiseAlpha(\n                            exponent=(-1, 0),\n                            first=iaa.Multiply((0.9, 1.1), per_channel=True),\n                            second=iaa.ContrastNormalization((0.9, 1.1))\n                        )\n                    ]),\n                    sometimes(iaa.ElasticTransformation(alpha=(0.5, 3.5), sigma=0.25)), # move pixels locally around (with random strengths)\n                    sometimes(iaa.PiecewiseAffine(scale=(0.01, 0.05))), # sometimes move parts of the image around\n                    sometimes(iaa.PerspectiveTransform(scale=(0.01, 0.1)))\n                ],\n                random_order=True\n            )\n        ],\n        random_order=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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    return ImageDataGenerator(\n        preprocessing_function=seq.augment_image\n    )\nif TRAINING:\n    # Using original generator\n    data_generator = create_datagen().flow(x_train, y_train, batch_size=BATCH_SIZE)\n    # Using Mixup\n    # mixup_generator = MixupGenerator(x_train, y_train, batch_size=BATCH_SIZE, alpha=0.2, datagen=create_datagen())()\n    mixup_generator = MixupGenerator(x_train, y_train, seq=seq, batch_size=BATCH_SIZE, alpha=0.2)()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    fig, ax = plt.subplots(5, 5, figsize=(15, 15))\n    indexes = np.random.randint(len(x_train), size=25)\n    for i in range(5):\n        for j in range(5):\n            ax[i][j].imshow(x_train[indexes[i * 5 + j]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Image and its augmentations\nif TRAINING:\n    index = np.random.randint(len(x_train), size=1)[0]\n\n    fig, ax = plt.subplots(1, 5, figsize=(15, 3))\n    for i in range(5):\n        if i == 0:\n            ax[i].imshow(x_train[index])\n            ax[i].axis(\"off\")\n        else:\n#             img = x_train[index] * 255\n#             img = img.astype(np.uint8)\n            img = seq.augment_image(x_train[index])\n            ax[i].imshow(img)\n            ax[i].axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"true_labels = np.array([1, 0, 1, 1, 0, 1])\npred_labels = np.array([1, 0, 0, 0, 0, 1])","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"accuracy_score(true_labels, pred_labels)","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"cohen_kappa_score(true_labels, pred_labels)","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"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).clip(0, 4)\n#         y_val = np.argmax(y_val, axis=1)\n        \n        y_pred = self.model.predict(X_val) > 0.5\n        y_pred = (y_pred.astype(int).sum(axis=1) - 1).clip(0, 4)\n#         y_pred = np.argmax(self.model.predict(X_val), axis=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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model: DenseNet-121"},{"metadata":{"trusted":true},"cell_type":"code","source":"base_model = Xception(\n    weights=None, #'imagenet',\n    include_top=False,\n    input_shape=(224,224,3)\n)\n\nbase_model.load_weights(\"../input/keras-pretrained-models/xception_weights_tf_dim_ordering_tf_kernels_notop.h5\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_model():\n    model = Sequential()\n    model.add(base_model)\n    model.add(layers.GlobalAveragePooling2D())\n    model.add(layers.Dropout(0.5))\n#     model.add(layers.Dense(1024, activation='sigmoid'))\n#     model.add(layers.Dropout(0.5))\n#     model.add(layers.BatchNormalization())\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = build_model()\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training & Evaluation"},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    kappa_metrics = Metrics()\n    rlr = callbacks.ReduceLROnPlateau(factor=0.5, patience=4, verbose=1)\n    es  = callbacks.EarlyStopping(patience=10, verbose=1, mode=\"min\")\n\n    history = model.fit_generator(\n        data_generator,\n        steps_per_epoch=x_train.shape[0] / BATCH_SIZE,\n        epochs=200,\n        validation_data=(x_val, y_val),\n        callbacks=[kappa_metrics, rlr, es]\n    )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    with open('history.json', 'w') as f:\n        json.dump(str(history.history), f)\n\n    history_df = pd.DataFrame(history.history)\n    history_df[['loss', 'val_loss']].plot()\n    history_df[['acc', 'val_acc']].plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    plt.plot(kappa_metrics.val_kappas)","execution_count":null,"outputs":[]},{"metadata":{},"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":{"trusted":true},"cell_type":"code","source":"# if TRAINING:\n#     model.load_weights('model.h5')\n#     y_val_pred = model.predict(x_val)\n\n#     def 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\n#     simplex = scipy.optimize.minimize(\n#         compute_score_inv, 0.5, method='nelder-mead'\n#     )\n\n#     best_threshold = simplex['x'][0]\n    \n#     print(best_threshold)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submit"},{"metadata":{"trusted":true},"cell_type":"code","source":"# if x_test.dtype == np.uint8:\n#     x_test = x_test.astype(np.float32)\n    \n# x_test = x_test / np.max(x_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if TRAINING:\n    weights_path = 'model.h5'\nelse:\n    weights_path = '../input/weights/model.h5'\nmodel.load_weights(weights_path)\n\ny_test = model.predict(x_test) > 0.5\ny_test = (y_test.astype(int).sum(axis=1) - 1).clip(0, 4)\n# y_test = np.argmax(model.predict(x_test), axis=1).astype(int)\n\ntest_df['diagnosis'] = y_test\ntest_df.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.hist()","execution_count":null,"outputs":[]}],"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"}},"nbformat":4,"nbformat_minor":1}