{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Implementation of EfficientNetB5 for the APTOS 2019 competition with Keras"},{"metadata":{},"cell_type":"markdown","source":"In this kernel we will implement EfficientNet for medical images ([APTOS 2019](https://www.kaggle.com/c/aptos2019-blindness-detection) competition). EfficientNet was released this June (2019) by Google AI and is the new state-of-the-art on ImageNet. It introduces a systematic way to scale CNN (Convolutional Neural Networks) in a nearly optimal way. For this kernel we will use the B5 version, but feel free to play with the larger models. This kernel provides weights for EfficientNetB0 through B5. Weights for EfficientNetB6 and B7 can be found in [Google AI's repository for EfficientNet](https://github.com/tensorflow/tpu/tree/master/models/official/efficientnet). I highly recommend you to read [the EfficientNet paper](https://arxiv.org/pdf/1905.11946.pdf) as it signifies a fundamental shift in how the Deep Learning community will approach model scaling!\n\nAlso, check out this [video on EfficientNet by Henry AI Labs](https://youtu.be/3svIm5UC94I) for a clear explanation!\n\nThis kernel is based heavily on Carlo Lepelaars' work, which can be found [here](https://www.kaggle.com/carlolepelaars/efficientnetb5-with-keras-aptos-2019). Thanks to Carlo for sharing his work and allowing everyone to learn from it. \n\nIf you like this Kaggle kernel, feel free to give an upvote and leave a comment!"},{"metadata":{},"cell_type":"markdown","source":"Image: an overview of model architectures and their performance on [ImageNet](http://www.image-net.org/). We can see that EfficientNet achieves state-of-the-art and uses a lot less parameters than most modern CNN architectures.\n\n[Source: EfficientNet Paper](https://arxiv.org/pdf/1905.11946.pdf)\n\n![](https://raw.githubusercontent.com/tensorflow/tpu/master/models/official/efficientnet/g3doc/params.png)"},{"metadata":{},"cell_type":"markdown","source":"## Table Of Contents"},{"metadata":{},"cell_type":"markdown","source":"- [Dependencies](#1)\n- [Preparation](#2)\n- [Metric (Quadratic Weighted Kappa)](#3)\n- [EDA (Exploratory Data Analysis)](#4)\n- [Preprocessing](#5)\n- [Modeling (EfficientNetB5)](#6)\n- [Evaluation](#7)\n- [Submission](#8)"},{"metadata":{},"cell_type":"markdown","source":"## Dependencies <a id=\"1\"></a>"},{"metadata":{},"cell_type":"markdown","source":"Special thanks to [qubvel](https://github.com/qubvel/efficientnet) for sharing an amazing wrapper to get the EfficientNet architecture in one line of code!"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import os\nimport sys\n# Repository source: https://github.com/qubvel/efficientnet\nsys.path.append(os.path.abspath('../input/efficientnet/efficientnet-master/efficientnet-master/'))\nfrom efficientnet import EfficientNetB5","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_kg_hide-input":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# Standard dependencies\nimport cv2\nimport time\nimport math\nimport scipy as sp\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom functools import partial\nimport matplotlib.pyplot as plt\n\n# Machine Learning\nimport tensorflow as tf\nimport keras\nfrom keras import initializers\nfrom keras import regularizers\nfrom keras import constraints\nfrom keras import backend as K\nfrom keras.activations import elu\nfrom keras.optimizers import Adam\nfrom keras.models import Sequential\nfrom keras.models import load_model\nfrom keras.engine import Layer, InputSpec\nfrom keras.utils.generic_utils import get_custom_objects\nfrom keras.callbacks import Callback, EarlyStopping, ReduceLROnPlateau\nfrom keras.layers import Dense, Conv2D, Flatten, GlobalAveragePooling2D, Dropout\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom sklearn.metrics import cohen_kappa_score\n\n# Path specifications\nKAGGLE_DIR = '../input/aptos2019-blindness-detection/'\nTRAIN_DF_PATH = KAGGLE_DIR + \"train.csv\"\nTEST_DF_PATH = KAGGLE_DIR + 'test.csv'\nTRAIN_IMG_PATH = KAGGLE_DIR + \"train_images/\"\nTEST_IMG_PATH = KAGGLE_DIR + 'test_images/'\nWEIGHTS_PATH = os.path.join(KAGGLE_DIR, '..','trained-models', 'effnet_modelB5.h5')\n\n# Specify title of our final model\nSAVED_MODEL_NAME = 'effnet_modelB5.h5'\n\n# Set seed for reproducability\nseed = 42\nnp.random.seed(seed)\ntf.set_random_seed(seed)\n\n# For keeping time. GPU limit for this competition is set to ± 9 hours.\nt_start = time.time()\n\n# File sizes and specifications\nprint('\\n# Files and file sizes')\nfor file in os.listdir(KAGGLE_DIR):\n    print('{}| {} MB'.format(file.ljust(30), \n                             str(round(os.path.getsize(KAGGLE_DIR + file) / 1000000, 2))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preparation <a id=\"2\"></a>"},{"metadata":{},"cell_type":"markdown","source":"By examining the data we can readily see that we do not have that much data (± 700 samples per class). It is probably a good idea to use data augmentation to increase robustness of our model (See the modeling section).\n\nWe could also try to use additional data from previous competitions to increase performance. Although I do not implement this in the kernel, feel free to experiment with adding data. Additional data can be found in [this Kaggle dataset](https://www.kaggle.com/benjaminwarner/resized-2015-2019-blindness-detection-images) (± 35000 additional images)."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Image IDs and Labels (TRAIN)\")\ntrain_df = pd.read_csv(TRAIN_DF_PATH)\n# Add extension to id_code\ntrain_df['id_code'] = train_df['id_code'] + \".png\"\nprint(f\"Training images: {train_df.shape[0]}\")\ndisplay(train_df.head())\nprint(\"Image IDs (TEST)\")\ntest_df = pd.read_csv(TEST_DF_PATH)\n# Add extension to id_code\ntest_df['id_code'] = test_df['id_code'] + \".png\"\nprint(f\"Testing Images: {test_df.shape[0]}\")\ndisplay(test_df.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The original image size from the [EfficientNet paper](https://arxiv.org/pdf/1905.11946.pdf) for EfficientNetB5 is 600x600x3. We are however not bound by this and can use a smaller size if we want. The original image sizes used for every version of EfficientNet are:\n\n- EfficientNetB0 - (224, 224, 3)\n- EfficientNetB1 - (240, 240, 3)\n- EfficientNetB2 - (260, 260, 3)\n- EfficientNetB3 - (300, 300, 3)\n- EfficientNetB4 - (380, 380, 3)\n- **EfficientNetB5 - (456, 456, 3)**\n- EfficientNetB6 - (528, 528, 3)\n- EfficientNetB7 - (600, 600, 3)\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Specify image size\nIMG_WIDTH = 456\nIMG_HEIGHT = 456\nIMAGE_SIZE = 456\nCHANNELS = 3","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Metric (Quadratic Weighted Kappa) <a id=\"3\"></a>"},{"metadata":{},"cell_type":"markdown","source":"The metric that is used for this competition is Quadratic Weighted Kappa (QWK) ([Kaggle's Explanation](https://www.kaggle.com/c/aptos2019-blindness-detection/overview/evaluation)) \n\nThe formula for weighted kappa is:\n\n![](https://wikimedia.org/api/rest_v1/media/math/render/svg/2a496e1cef7d812b83bdbb725d291748cf0183f5)\n\nIn this case we are going to optimize Mean Squared Error (MSE) (See Modeling section) since we are using regression and by optimizing MSE we are also optimizing QWK as long as we round predictions afterwards. Additionally we are going to same the model which achieves the best QWK score on the validation data through a custom Keras Callback.\n\nFor a more detailed and practical explanation of QWK I highly recommend [this Kaggle kernel](https://www.kaggle.com/aroraaman/quadratic-kappa-metric-explained-in-5-simple-steps)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_preds_and_labels(model, generator):\n    \"\"\"\n    Get predictions and labels from the generator\n    \"\"\"\n    preds = []\n    labels = []\n    for _ in range(int(np.ceil(generator.samples / BATCH_SIZE))):\n        x, y = next(generator)\n        preds.append(model.predict(x))\n        labels.append(y)\n    # Flatten list of numpy arrays\n    return np.concatenate(preds).ravel(), np.concatenate(labels).ravel()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Metrics(Callback):\n    \"\"\"\n    A custom Keras callback for saving the best model\n    according to the Quadratic Weighted Kappa (QWK) metric\n    \"\"\"\n    def on_train_begin(self, logs={}):\n        \"\"\"\n        Initialize list of QWK scores on validation data\n        \"\"\"\n        self.val_kappas = []\n\n    def on_epoch_end(self, epoch, logs={}):\n        \"\"\"\n        Gets QWK score on the validation data\n        \"\"\"\n        # Get predictions and convert to integers\n        y_pred, labels = get_preds_and_labels(model, val_generator)\n        y_pred = np.rint(y_pred).astype(np.uint8).clip(0, 4)\n        # We can use sklearns implementation of QWK straight out of the box\n        # as long as we specify weights as 'quadratic'\n        _val_kappa = cohen_kappa_score(labels, y_pred, weights='quadratic')\n        self.val_kappas.append(_val_kappa)\n        print(f\"val_kappa: {round(_val_kappa, 4)}\")\n        if _val_kappa == max(self.val_kappas):\n            print(\"Validation Kappa has improved. Saving model.\")\n            self.model.save(SAVED_MODEL_NAME)\n        return","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## EDA (Exploratory Data Analysis) <a id=\"4\"></a>"},{"metadata":{},"cell_type":"markdown","source":"For EDA on image datasets I think one should at least examine the label distribution, the images before preprocessing and the images after preprocessing. Through examining these three aspects we can get a good sense of the problem. Note that the distribution on the test set can still vary wildly from the training data."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Label distribution\ntrain_df['diagnosis'].value_counts().sort_index().plot(kind=\"bar\", \n                                                       figsize=(12,5), \n                                                       rot=0)\nplt.title(\"Label Distribution (Training Set)\", \n          weight='bold', \n          fontsize=18)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Label\", fontsize=17)\nplt.ylabel(\"Frequency\", fontsize=17);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will visualize a random image from every label to get a general sense of the distinctive features that seperate the classes. We will take this into account and try to enhance these features in our preprocessing. For these images there some to be increasingly more spots and stains on the retina as diabetic retinopathy worsens."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Example from every label\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\n\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\")\n    X = cv2.cvtColor(X, cv2.COLOR_BGR2RGB)\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preprocessing <a id=\"5\"></a>"},{"metadata":{},"cell_type":"markdown","source":"Here we will use the auto-cropping method with Ben's preprocessing as explained in [this kernel](https://www.kaggle.com/ratthachat/aptos-updatedv14-preprocessing-ben-s-cropping).\n\nIn addition, we will also include the circle cropping method to further standardize the images.From the [kernel here](https://www.kaggle.com/taindow/pre-processing-train-and-test-images): \"Here we resize the image after the first crop before drawing the circle. The result is that a larger portion of the zoomed in images is retained, though it will be somewhat stretched. I think a slightly stretch is better than losing so much information.\"\n\nThere are also 3 more potential preprocessing methods below. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Version 1 Preprocessing code (Ben's + Circular Cropping)\n\ndef crop_image_from_gray(img, tol=7):\n    \"\"\"\n    Applies masks to the orignal image and \n    returns the a preprocessed image with \n    3 channels\n    \"\"\"\n    # If for some reason we only have two channels\n    if img.ndim == 2:\n        mask = img > tol\n        return img[np.ix_(mask.any(1),mask.any(0))]\n    # If we have a normal RGB images\n    elif img.ndim == 3:\n        gray_img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)\n        mask = gray_img > tol\n        \n        check_shape = img[:,:,0][np.ix_(mask.any(1),mask.any(0))].shape[0]\n        if (check_shape == 0): # image is too dark so that we crop out everything,\n            return img # return original image\n        else:\n            img1=img[:,:,0][np.ix_(mask.any(1),mask.any(0))]\n            img2=img[:,:,1][np.ix_(mask.any(1),mask.any(0))]\n            img3=img[:,:,2][np.ix_(mask.any(1),mask.any(0))]\n            img = np.stack([img1,img2,img3],axis=-1)\n        return img\n    \ndef circle_crop_v2(img):\n    \"\"\"\n    Create circular crop around image centre\n    \"\"\"\n    #No need to cv2.imread the image as it has been done. \n    \n    height, width, depth = img.shape\n    largest_side = np.max((height, width))\n    img = cv2.resize(img, (largest_side, largest_side))\n\n    height, width, depth = img.shape\n\n    x = int(width / 2)\n    y = int(height / 2)\n    r = np.amin((x, y))\n\n    circle_img = np.zeros((height, width), np.uint8)\n    cv2.circle(circle_img, (x, y), int(r), 1, thickness=-1)\n    img = cv2.bitwise_and(img, img, mask=circle_img)\n    img = crop_image_from_gray(img)\n\n    return img\n\ndef preprocess_image(image, sigmaX=25):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Apply masks\n    3. Apply circular cropping to standardize images \n    4. Resize image to desired size\n    5. Remove local average colour to bring out features\n    \"\"\"\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    image = crop_image_from_gray(image)\n    image = circle_crop_v2(image)\n    image = cv2.resize(image, (IMG_WIDTH, IMG_HEIGHT))\n    image = cv2.addWeighted (image,4, cv2.GaussianBlur(image, (0,0) ,sigmaX), -4, 128)\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Version 2 Preprocessing Code (Image Centre Resizing + Gaussian Background Subtraction)\n\ndef info_image(im):\n    # Compute the center (cx, cy) and radius of the eye\n    cy = im.shape[0]//2\n    midline = im[cy,:]\n    midline = np.where(midline>midline.mean()/3)[0]\n    if len(midline)>im.shape[1]//2:\n        x_start, x_end = np.min(midline), np.max(midline)\n    else: # This actually rarely happens p~1/10000\n        x_start, x_end = im.shape[1]//10, 9*im.shape[1]//10\n    cx = (x_start + x_end)/2\n    r = (x_end - x_start)/2\n    return cx, cy, r\n\ndef resize_image(im, augmentation=False):\n    # Crops, resizes and potentially augments the image to IMAGE_SIZE\n    cx, cy, r = info_image(im)\n    scaling = IMAGE_SIZE/(2*r)\n    rotation = 0\n    if augmentation:\n        scaling *= 1 + 0.3 * (np.random.rand()-0.5)\n        rotation = 360 * np.random.rand()\n    M = cv2.getRotationMatrix2D((cx,cy), rotation, scaling)\n    M[0,2] -= cx - IMAGE_SIZE/2\n    M[1,2] -= cy - IMAGE_SIZE/2\n    return cv2.warpAffine(im,M,(IMAGE_SIZE,IMAGE_SIZE)) # This is the most important line\n\n#Error occurs with cv2.medianBlur\ndef subtract_median_bg_image(im):\n    #im = im / (np.max(im.shape))\n    #im = img_as_ubyte(im)\n    #k = np.max(im.shape)//20*2+1\n    k = 45\n    bg = cv2.medianBlur(im, k)\n    return cv2.addWeighted (im, 4, bg, -4, 128)\n\ndef subtract_gaussian_bg_image(im):\n    k = np.max(im.shape)/15\n    bg = cv2.GaussianBlur(im ,(0,0) ,k)\n    return cv2.addWeighted (im, 4, bg, -4, 128)\n\n#Constrast stretching not working well at this moment \ndef contrast_stretching(img):        \n    rr, gg, bb = cv2.split(img)    \n    imgray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)    \n    im = imgray    \n    ih, iw = imgray.shape    \n    (minVal, maxVal, minLoc, maxLoc) = cv2.minMaxLoc(imgray)    \n    for i in range(ih):        \n        for j in range(iw):            \n            im[i, j] = 255 * ((gg[i, j] - minVal) / (maxVal - minVal))        \n    limg = cv2.merge((rr, im, bb))    \n    return limg\n\n#Histogram normalizaiton not working well at this moment \ndef histogram_normalization(image):    \n    hist,bins = np.histogram(image.flatten(),256,[0,256])    \n    cdf = hist.cumsum()   \n    # cdf_normalized = cdf * hist.max()/ cdf.max()    \n    cdf_m = np.ma.masked_equal(cdf,0)    \n    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())    \n    cdf = np.ma.filled(cdf_m,0).astype('uint8')     \n    img2 = cdf[image]    \n    return img2\n\n# To remove irregularities along the circular boundary of the image\nPARAM = 96\ndef Radius_Reduction(img,PARAM):\n    h,w,c=img.shape\n    Frame=np.zeros((h,w,c),dtype=np.uint8)\n    cv2.circle(Frame,(int(math.floor(w/2)),int(math.floor(h/2))),int(math.floor((h*PARAM)/float(2*100))), (255,255,255), -1)\n    Frame1=cv2.cvtColor(Frame, cv2.COLOR_BGR2GRAY)\n    img1 =cv2.bitwise_and(img,img,mask=Frame1)\n    return img1\n\ndef preprocess_imagev2(image, sigmaX=10):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Find image centre and resize image\n    3. Subtract the Gaussian background (remove local average colour)  \n    4. Remove irregularities along the circular boundary of the image\n    \"\"\"\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    image = resize_image(image)\n    image = subtract_gaussian_bg_image(image)\n    image = Radius_Reduction(image, PARAM)\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Version 3 Preprocessing Code (Contrast Limited Adaptive Histogram Equalization)\n\ndef info_image(im):\n    # Compute the center (cx, cy) and radius of the eye\n    cy = im.shape[0]//2\n    midline = im[cy,:]\n    midline = np.where(midline>midline.mean()/3)[0]\n    if len(midline)>im.shape[1]//2:\n        x_start, x_end = np.min(midline), np.max(midline)\n    else: # This actually rarely happens p~1/10000\n        x_start, x_end = im.shape[1]//10, 9*im.shape[1]//10\n    cx = (x_start + x_end)/2\n    r = (x_end - x_start)/2\n    return cx, cy, r\n\ndef resize_image(im, augmentation=False):\n    # Crops, resizes and potentially augments the image to IMAGE_SIZE\n    cx, cy, r = info_image(im)\n    scaling = IMAGE_SIZE/(2*r)\n    rotation = 0\n    if augmentation:\n        scaling *= 1 + 0.3 * (np.random.rand()-0.5)\n        rotation = 360 * np.random.rand()\n    M = cv2.getRotationMatrix2D((cx,cy), rotation, scaling)\n    M[0,2] -= cx - IMAGE_SIZE/2\n    M[1,2] -= cy - IMAGE_SIZE/2\n    return cv2.warpAffine(im,M,(IMAGE_SIZE,IMAGE_SIZE)) # This is the most important line\n\ndef adaptive_histo_equalization(res_image):\n    lab= cv2.cvtColor(res_image, cv2.COLOR_BGR2LAB)\n    l, a, b = cv2.split(lab)\n    clahe = cv2.createCLAHE(clipLimit=4.0, tileGridSize=(8,8))\n    cl = clahe.apply(l)\n    limg = cv2.merge((cl,a,b))\n    final = cv2.cvtColor(limg, cv2.COLOR_LAB2BGR)\n    return final \n\ndef preprocess_imagev3(image, sigmaX=10):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Find image centre and resize image\n    3. Conduct CLAHE    \n    \"\"\"\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    image = resize_image(image)\n    image = adaptive_histo_equalization(image)\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Version 4 Preprocessing code (V1 + Kirsch Operator + Merging)\n\ndef Kirsch(crop):    \n    Input=crop[:,:,1]    \n    a,b=Input.shape    \n    Kernel=np.zeros((3,3,8)) #windows declearations(8 windows)    \n    Kernel[:,:,0]=np.array([[5,5,5],[-3,0,-3],[-3,-3,-3]])     \n    Kernel[:,:,1]=np.array([[-3,5,5],[-3,0,5],[-3,-3,-3]])    \n    Kernel[:,:,2]=np.array([[-3,-3,5],[-3,0,5],[-3,-3,5]])    \n    Kernel[:,:,3]=np.array([[-3,-3,-3],[-3,0,5],[-3,5,5]])    \n    Kernel[:,:,4]=np.array([[-3,-3,-3],[-3,0,-3],[5,5,5]])    \n    Kernel[:,:,5]=np.array([[-3,-3,-3],[5,0,-3],[5,5,-3]])    \n    Kernel[:,:,6]=np.array([[5,-3,-3],[5,0,-3],[5,-3,-3]])    \n    Kernel[:,:,7]=np.array([[5,5,-3],[5,0,-3],[-3,-3,-3]])    \n    #Kernel=(1/float(15))*Kernel    \n    #Convolution output    \n    dst=np.zeros((a,b,8))    \n    for x in range(0,8):        \n        dst[:,:,x] = cv2.filter2D(Input,-1,Kernel[:,:,x])    \n    Out=np.zeros((a,b))    \n    for y in range(0,a-1):        \n        for z in range(0,b-1):            \n            Out[y,z]=max(dst[y,z,:])    \n    Out=np.uint8(Out)            \n    return Out\n\ndef preprocess_imagev4(image):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Apply masks for intial cropping\n    3. Apply circular cropping to standardize images \n    4. Resize image to desired size\n    5. Subtract local background with larger SigmaX\n    \"\"\"\n    image = preprocess_image(image, sigmaX=10)\n    image_edge = Kirsch(image)\n    return image_edge","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Version 5 Preprocessing Code (V3 + V1)\n\ndef info_image(im):\n    # Compute the center (cx, cy) and radius of the eye\n    cy = im.shape[0]//2\n    midline = im[cy,:]\n    midline = np.where(midline>midline.mean()/3)[0]\n    if len(midline)>im.shape[1]//2:\n        x_start, x_end = np.min(midline), np.max(midline)\n    else: # This actually rarely happens p~1/10000\n        x_start, x_end = im.shape[1]//10, 9*im.shape[1]//10\n    cx = (x_start + x_end)/2\n    r = (x_end - x_start)/2\n    return cx, cy, r\n\ndef resize_image(im, augmentation=False):\n    # Crops, resizes and potentially augments the image to IMAGE_SIZE\n    cx, cy, r = info_image(im)\n    scaling = IMAGE_SIZE/(2*r)\n    rotation = 0\n    if augmentation:\n        scaling *= 1 + 0.3 * (np.random.rand()-0.5)\n        rotation = 360 * np.random.rand()\n    M = cv2.getRotationMatrix2D((cx,cy), rotation, scaling)\n    M[0,2] -= cx - IMAGE_SIZE/2\n    M[1,2] -= cy - IMAGE_SIZE/2\n    return cv2.warpAffine(im,M,(IMAGE_SIZE,IMAGE_SIZE)) # This is the most important line\n\ndef adaptive_histo_equalization(res_image):\n    lab= cv2.cvtColor(res_image, cv2.COLOR_BGR2LAB)\n    l, a, b = cv2.split(lab)\n    l = l.astype(np.uint8)\n    a = a.astype(np.uint8)\n    b = b.astype(np.uint8)\n    clahe = cv2.createCLAHE(clipLimit=4.0, tileGridSize=(8,8))\n    cl = clahe.apply(l)\n    limg = cv2.merge((cl,a,b))\n    final = cv2.cvtColor(limg, cv2.COLOR_LAB2BGR)\n    return final \n\ndef preprocess_imagev5(image, sigmaX=5):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Find image centre and resize image\n    3. Conduct CLAHE\n    4. Remove local average colour to bring out features \n    \"\"\"\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    image = crop_image_from_gray(image)\n    image = circle_crop_v2(image)\n    image = cv2.resize(image, (IMG_WIDTH, IMG_HEIGHT))\n    #image = resize_image(image)\n    image = adaptive_histo_equalization(image)\n    image = cv2.addWeighted (image,4, cv2.GaussianBlur(image, (0,0) ,sigmaX), -4, 128)\n    return image.astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After preprocessing we have managed to enhance the distinctive features in the images. This will increase performance when we train our EfficientNet model."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Example of preprocessed images from every label, using v1 preprocessing\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_image(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example of preprocessed images from every label, using v2 preprocessing\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_imagev2(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example of preprocessed images from every label, using v3 preprocessing\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_imagev3(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example of preprocessed images from every label, using v4 preprocessing\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_imagev4(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example of preprocessed images from every label, using v5 preprocessing\nfig, ax = plt.subplots(1, 5, figsize=(30, 12))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_imagev5(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=15)\n    ax[i].axis('off')\n    ax[i].imshow(X);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Labels for training data\ny_labels = train_df['diagnosis'].values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Modeling (EfficientNetB5) <a id=\"6\"></a>"},{"metadata":{},"cell_type":"markdown","source":"Since we want to optimize the Quadratic Weighted Kappa score we can formulate this challenge as a regression problem. In this way we are more flexible in our optimization and we can yield higher scores than solely optimizing for accuracy. We will optimize a pre-trained EfficientNetB5 with a few added layers. The metric that we try to optimize is the [Mean Squared Error](https://en.wikipedia.org/wiki/Mean_squared_error). This is the mean of squared differences between our predictions and labels, as showed in the formula below. By optimizing this metric we are also optimizing for Quadratic Weighted Kappa if we round the predictions afterwards.\n\n![](https://study.com/cimages/multimages/16/4e7cf150-0179-4d89-86f2-5cbb1f51c266_meansquarederrorformula.png)\n\nSince we are not provided with that much data (3662 images), we will augment the data to make the model more robust. We will rotate the data on any angle. Also, we will flip the data both horizontally and vertically. Lastly, we will divide the data by 128 for normalization."},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 4\n\n# Add Image augmentation to our generator, we are using version 1 preprocessing here. \ntrain_datagen = ImageDataGenerator(rotation_range=360,\n                                   horizontal_flip=True,\n                                   vertical_flip=True,\n                                   validation_split=0.10,\n                                   preprocessing_function=preprocess_image, \n                                   rescale=1 / 128.\n                                  )\n\n# Use the dataframe to define train and validation generators\ntrain_generator = train_datagen.flow_from_dataframe(train_df, \n                                                    x_col='id_code', \n                                                    y_col='diagnosis',\n                                                    directory = TRAIN_IMG_PATH,\n                                                    target_size=(IMG_WIDTH, IMG_HEIGHT),\n                                                    batch_size=BATCH_SIZE,\n                                                    class_mode='other', \n                                                    subset='training')\n\nval_generator = train_datagen.flow_from_dataframe(train_df, \n                                                  x_col='id_code', \n                                                  y_col='diagnosis',\n                                                  directory = TRAIN_IMG_PATH,\n                                                  target_size=(IMG_WIDTH, IMG_HEIGHT),\n                                                  batch_size=BATCH_SIZE,\n                                                  class_mode='other',\n                                                  subset='validation')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Thanks to the amazing wrapper by [qubvel](https://github.com/qubvel/efficientnet) we can load in a model like the Keras API. We specify the input shape and that we want the model without the top (the final Dense layer). Then we load in the weights which are provided in [this Kaggle dataset](https://www.kaggle.com/ratthachat/efficientnet-keras-weights-b0b5). Note that we will use the [RAdam optimizer](https://arxiv.org/pdf/1908.03265v1.pdf) since it often yields better convergence than Vanilla Adam. Thanks to CyberZHG who implemented [RAdam for Keras](https://github.com/CyberZHG/keras-radam/blob/master/keras_radam/optimizers.py)."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Code Source: https://github.com/CyberZHG/keras-radam/blob/master/keras_radam/optimizers.py\nclass RAdam(keras.optimizers.Optimizer):\n    \"\"\"RAdam optimizer.\n    # Arguments\n        lr: float >= 0. Learning rate.\n        beta_1: float, 0 < beta < 1. Generally close to 1.\n        beta_2: float, 0 < beta < 1. Generally close to 1.\n        epsilon: float >= 0. Fuzz factor. If `None`, defaults to `K.epsilon()`.\n        decay: float >= 0. Learning rate decay over each update.\n        weight_decay: float >= 0. Weight decay for each param.\n        amsgrad: boolean. Whether to apply the AMSGrad variant of this\n            algorithm from the paper \"On the Convergence of Adam and\n            Beyond\".\n        total_steps: int >= 0. Total number of training steps. Enable warmup by setting a positive value.\n        warmup_proportion: 0 < warmup_proportion < 1. The proportion of increasing steps.\n        min_lr: float >= 0. Minimum learning rate after warmup.\n    # References\n        - [Adam - A Method for Stochastic Optimization](https://arxiv.org/abs/1412.6980v8)\n        - [On the Convergence of Adam and Beyond](https://openreview.net/forum?id=ryQu7f-RZ)\n        - [On The Variance Of The Adaptive Learning Rate And Beyond](https://arxiv.org/pdf/1908.03265v1.pdf)\n    \"\"\"\n\n    def __init__(self, lr=0.001, beta_1=0.9, beta_2=0.999,\n                 epsilon=None, decay=0., weight_decay=0., amsgrad=False,\n                 total_steps=0, warmup_proportion=0.1, min_lr=0., **kwargs):\n        super(RAdam, self).__init__(**kwargs)\n        with K.name_scope(self.__class__.__name__):\n            self.iterations = K.variable(0, dtype='int64', name='iterations')\n            self.lr = K.variable(lr, name='lr')\n            self.beta_1 = K.variable(beta_1, name='beta_1')\n            self.beta_2 = K.variable(beta_2, name='beta_2')\n            self.decay = K.variable(decay, name='decay')\n            self.weight_decay = K.variable(weight_decay, name='weight_decay')\n            self.total_steps = K.variable(total_steps, name='total_steps')\n            self.warmup_proportion = K.variable(warmup_proportion, name='warmup_proportion')\n            self.min_lr = K.variable(lr, name='min_lr')\n        if epsilon is None:\n            epsilon = K.epsilon()\n        self.epsilon = epsilon\n        self.initial_decay = decay\n        self.initial_weight_decay = weight_decay\n        self.initial_total_steps = total_steps\n        self.amsgrad = amsgrad\n\n    def get_updates(self, loss, params):\n        grads = self.get_gradients(loss, params)\n        self.updates = [K.update_add(self.iterations, 1)]\n\n        lr = self.lr\n\n        if self.initial_decay > 0:\n            lr = lr * (1. / (1. + self.decay * K.cast(self.iterations, K.dtype(self.decay))))\n\n        t = K.cast(self.iterations, K.floatx()) + 1\n\n        if self.initial_total_steps > 0:\n            warmup_steps = self.total_steps * self.warmup_proportion\n            decay_steps = self.total_steps - warmup_steps\n            lr = K.switch(\n                t <= warmup_steps,\n                lr * (t / warmup_steps),\n                lr * (1.0 - K.minimum(t, decay_steps) / decay_steps),\n            )\n\n        ms = [K.zeros(K.int_shape(p), dtype=K.dtype(p), name='m_' + str(i)) for (i, p) in enumerate(params)]\n        vs = [K.zeros(K.int_shape(p), dtype=K.dtype(p), name='v_' + str(i)) for (i, p) in enumerate(params)]\n\n        if self.amsgrad:\n            vhats = [K.zeros(K.int_shape(p), dtype=K.dtype(p), name='vhat_' + str(i)) for (i, p) in enumerate(params)]\n        else:\n            vhats = [K.zeros(1, name='vhat_' + str(i)) for i in range(len(params))]\n\n        self.weights = [self.iterations] + ms + vs + vhats\n\n        beta_1_t = K.pow(self.beta_1, t)\n        beta_2_t = K.pow(self.beta_2, t)\n\n        sma_inf = 2.0 / (1.0 - self.beta_2) - 1.0\n        sma_t = sma_inf - 2.0 * t * beta_2_t / (1.0 - beta_2_t)\n\n        for p, g, m, v, vhat in zip(params, grads, ms, vs, vhats):\n            m_t = (self.beta_1 * m) + (1. - self.beta_1) * g\n            v_t = (self.beta_2 * v) + (1. - self.beta_2) * K.square(g)\n\n            m_corr_t = m_t / (1.0 - beta_1_t)\n            if self.amsgrad:\n                vhat_t = K.maximum(vhat, v_t)\n                v_corr_t = K.sqrt(vhat_t / (1.0 - beta_2_t) + self.epsilon)\n                self.updates.append(K.update(vhat, vhat_t))\n            else:\n                v_corr_t = K.sqrt(v_t / (1.0 - beta_2_t) + self.epsilon)\n\n            r_t = K.sqrt((sma_t - 4.0) / (sma_inf - 4.0) *\n                         (sma_t - 2.0) / (sma_inf - 2.0) *\n                         sma_inf / sma_t)\n\n            p_t = K.switch(sma_t > 5, r_t * m_corr_t / v_corr_t, m_corr_t)\n\n            if self.initial_weight_decay > 0:\n                p_t += self.weight_decay * p\n\n            p_t = p - lr * p_t\n\n            self.updates.append(K.update(m, m_t))\n            self.updates.append(K.update(v, v_t))\n            new_p = p_t\n\n            # Apply constraints.\n            if getattr(p, 'constraint', None) is not None:\n                new_p = p.constraint(new_p)\n\n            self.updates.append(K.update(p, new_p))\n        return self.updates\n\n    def get_config(self):\n        config = {\n            'lr': float(K.get_value(self.lr)),\n            'beta_1': float(K.get_value(self.beta_1)),\n            'beta_2': float(K.get_value(self.beta_2)),\n            'decay': float(K.get_value(self.decay)),\n            'weight_decay': float(K.get_value(self.weight_decay)),\n            'epsilon': self.epsilon,\n            'amsgrad': self.amsgrad,\n            'total_steps': float(K.get_value(self.total_steps)),\n            'warmup_proportion': float(K.get_value(self.warmup_proportion)),\n            'min_lr': float(K.get_value(self.min_lr)),\n        }\n        base_config = super(RAdam, self).get_config()\n        return dict(list(base_config.items()) + list(config.items()))\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"class GroupNormalization(Layer):\n    \"\"\"Group normalization layer\n    Group Normalization divides the channels into groups and computes within each group\n    the mean and variance for normalization. GN's computation is independent of batch sizes,\n    and its accuracy is stable in a wide range of batch sizes\n    # Arguments\n        groups: Integer, the number of groups for Group Normalization.\n        axis: Integer, the axis that should be normalized\n            (typically the features axis).\n            For instance, after a `Conv2D` layer with\n            `data_format=\"channels_first\"`,\n            set `axis=1` in `BatchNormalization`.\n        epsilon: Small float added to variance to avoid dividing by zero.\n        center: If True, add offset of `beta` to normalized tensor.\n            If False, `beta` is ignored.\n        scale: If True, multiply by `gamma`.\n            If False, `gamma` is not used.\n            When the next layer is linear (also e.g. `nn.relu`),\n            this can be disabled since the scaling\n            will be done by the next layer.\n        beta_initializer: Initializer for the beta weight.\n        gamma_initializer: Initializer for the gamma weight.\n        beta_regularizer: Optional regularizer for the beta weight.\n        gamma_regularizer: Optional regularizer for the gamma weight.\n        beta_constraint: Optional constraint for the beta weight.\n        gamma_constraint: Optional constraint for the gamma weight.\n    # Input shape\n        Arbitrary. Use the keyword argument `input_shape`\n        (tuple of integers, does not include the samples axis)\n        when using this layer as the first layer in a model.\n    # Output shape\n        Same shape as input.\n    # References\n        - [Group Normalization](https://arxiv.org/abs/1803.08494)\n    \"\"\"\n\n    def __init__(self,\n                 groups=32,\n                 axis=-1,\n                 epsilon=1e-5,\n                 center=True,\n                 scale=True,\n                 beta_initializer='zeros',\n                 gamma_initializer='ones',\n                 beta_regularizer=None,\n                 gamma_regularizer=None,\n                 beta_constraint=None,\n                 gamma_constraint=None,\n                 **kwargs):\n        super(GroupNormalization, self).__init__(**kwargs)\n        self.supports_masking = True\n        self.groups = groups\n        self.axis = axis\n        self.epsilon = epsilon\n        self.center = center\n        self.scale = scale\n        self.beta_initializer = initializers.get(beta_initializer)\n        self.gamma_initializer = initializers.get(gamma_initializer)\n        self.beta_regularizer = regularizers.get(beta_regularizer)\n        self.gamma_regularizer = regularizers.get(gamma_regularizer)\n        self.beta_constraint = constraints.get(beta_constraint)\n        self.gamma_constraint = constraints.get(gamma_constraint)\n\n    def build(self, input_shape):\n        dim = input_shape[self.axis]\n\n        if dim is None:\n            raise ValueError('Axis ' + str(self.axis) + ' of '\n                             'input tensor should have a defined dimension '\n                             'but the layer received an input with shape ' +\n                             str(input_shape) + '.')\n\n        if dim < self.groups:\n            raise ValueError('Number of groups (' + str(self.groups) + ') cannot be '\n                             'more than the number of channels (' +\n                             str(dim) + ').')\n\n        if dim % self.groups != 0:\n            raise ValueError('Number of groups (' + str(self.groups) + ') must be a '\n                             'multiple of the number of channels (' +\n                             str(dim) + ').')\n\n        self.input_spec = InputSpec(ndim=len(input_shape),\n                                    axes={self.axis: dim})\n        shape = (dim,)\n\n        if self.scale:\n            self.gamma = self.add_weight(shape=shape,\n                                         name='gamma',\n                                         initializer=self.gamma_initializer,\n                                         regularizer=self.gamma_regularizer,\n                                         constraint=self.gamma_constraint)\n        else:\n            self.gamma = None\n        if self.center:\n            self.beta = self.add_weight(shape=shape,\n                                        name='beta',\n                                        initializer=self.beta_initializer,\n                                        regularizer=self.beta_regularizer,\n                                        constraint=self.beta_constraint)\n        else:\n            self.beta = None\n        self.built = True\n\n    def call(self, inputs, **kwargs):\n        input_shape = K.int_shape(inputs)\n        tensor_input_shape = K.shape(inputs)\n\n        # Prepare broadcasting shape.\n        reduction_axes = list(range(len(input_shape)))\n        del reduction_axes[self.axis]\n        broadcast_shape = [1] * len(input_shape)\n        broadcast_shape[self.axis] = input_shape[self.axis] // self.groups\n        broadcast_shape.insert(1, self.groups)\n\n        reshape_group_shape = K.shape(inputs)\n        group_axes = [reshape_group_shape[i] for i in range(len(input_shape))]\n        group_axes[self.axis] = input_shape[self.axis] // self.groups\n        group_axes.insert(1, self.groups)\n\n        # reshape inputs to new group shape\n        group_shape = [group_axes[0], self.groups] + group_axes[2:]\n        group_shape = K.stack(group_shape)\n        inputs = K.reshape(inputs, group_shape)\n\n        group_reduction_axes = list(range(len(group_axes)))\n        group_reduction_axes = group_reduction_axes[2:]\n\n        mean = K.mean(inputs, axis=group_reduction_axes, keepdims=True)\n        variance = K.var(inputs, axis=group_reduction_axes, keepdims=True)\n\n        inputs = (inputs - mean) / (K.sqrt(variance + self.epsilon))\n\n        # prepare broadcast shape\n        inputs = K.reshape(inputs, group_shape)\n        outputs = inputs\n\n        # In this case we must explicitly broadcast all parameters.\n        if self.scale:\n            broadcast_gamma = K.reshape(self.gamma, broadcast_shape)\n            outputs = outputs * broadcast_gamma\n\n        if self.center:\n            broadcast_beta = K.reshape(self.beta, broadcast_shape)\n            outputs = outputs + broadcast_beta\n\n        outputs = K.reshape(outputs, tensor_input_shape)\n\n        return outputs\n\n    def get_config(self):\n        config = {\n            'groups': self.groups,\n            'axis': self.axis,\n            'epsilon': self.epsilon,\n            'center': self.center,\n            'scale': self.scale,\n            'beta_initializer': initializers.serialize(self.beta_initializer),\n            'gamma_initializer': initializers.serialize(self.gamma_initializer),\n            'beta_regularizer': regularizers.serialize(self.beta_regularizer),\n            'gamma_regularizer': regularizers.serialize(self.gamma_regularizer),\n            'beta_constraint': constraints.serialize(self.beta_constraint),\n            'gamma_constraint': constraints.serialize(self.gamma_constraint)\n        }\n        base_config = super(GroupNormalization, self).get_config()\n        return dict(list(base_config.items()) + list(config.items()))\n\n    def compute_output_shape(self, input_shape):\n        return input_shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load in EfficientNetB5\neffnet = EfficientNetB5(weights=None,\n                        include_top=False,\n                        input_shape=(IMG_WIDTH, IMG_HEIGHT, CHANNELS))\neffnet.load_weights('../input/efficientnet-keras-weights-b0b5/efficientnet-b5_imagenet_1000_notop.h5')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Batch Normalization becomes unstable with small batch sizes (<16) and that is why we use [Group Normalization ](https://arxiv.org/pdf/1803.08494.pdf) layers instead. Big thanks to [Somshubra Majumdar](https://github.com/titu1994) for building an implementation of Group Normalization for Keras.\n\nKeras makes it incredibly easy to replace layers. Just loop through the layers and replace each Batch Normalization layer with a Group Normalization layer."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Replace all Batch Normalization layers by Group Normalization layers\nfor i, layer in enumerate(effnet.layers):\n    if \"batch_normalization\" in layer.name:\n        effnet.layers[i] = GroupNormalization(groups=32, axis=-1, epsilon=0.00001)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_model():\n    \"\"\"\n    A custom implementation of EfficientNetB5\n    for the APTOS 2019 competition\n    (Regression)\n    \"\"\"\n    model = Sequential()\n    model.add(effnet)\n    model.add(GlobalAveragePooling2D())\n    model.add(Dropout(0.5))\n    model.add(Dense(5, activation=elu))\n    model.add(Dense(1, activation=\"linear\"))\n    model.compile(loss='mse',\n                  optimizer=RAdam(lr=0.00005), \n                  metrics=['mse', 'acc'])\n    print(model.summary())\n    return model\n\n# Initialize model\nmodel = build_model()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Setting class weights in order to emphasize the minority classes \n\nclass_weight = {0: 1.,\n                1: 5.,\n                2: 1.75,\n                3: 5.,\n                4: 5.}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We train all layers in the network. This is the traditional transfer learning approach were we can optimize and model for almost any image content. Since the pre-trained model was trained on [ImageNet](http://www.image-net.org/) and not on medical images, there are some limitations to this approach for our challenge.\n\nAfter each epoch we save the model if it is better than the previous one, according to the Quadratic Weighted Kappa score on the validation set. We also stop training if the MSE on the validation set doesn't go down for 4 epochs. This way we can counter overfitting.\n\nAnother option we could use is to directly use Quadratic Weighted Kappa as a loss function. Feel free to experiment with this. An implementation of a [QWK loss function for Tensorflow/Keras can be found in this Kaggle kernel](https://www.kaggle.com/christofhenkel/weighted-kappa-loss-for-keras-tensorflow)."},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# For tracking Quadratic Weighted Kappa score\nkappa_metrics = Metrics()\n# Monitor MSE to avoid overfitting and save best model\nes = EarlyStopping(monitor='val_loss', mode='auto', verbose=1, patience=12)\nrlr = ReduceLROnPlateau(monitor='val_loss', \n                        factor=0.5, \n                        patience=4, \n                        verbose=1, \n                        mode='auto', \n                        epsilon=0.0001)\n\n# Begin training\nmodel.fit_generator(train_generator,\n                    steps_per_epoch=train_generator.samples // BATCH_SIZE,\n                    epochs=25,\n                    validation_data=val_generator,\n                    validation_steps = val_generator.samples // BATCH_SIZE,\n                    class_weight = class_weight,\n                    callbacks=[kappa_metrics, es, rlr])","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Visualize mse\nhistory_df = pd.DataFrame(model.history.history)\nhistory_df[['loss', 'val_loss']].plot(figsize=(12,5))\nplt.title(\"Loss (MSE)\", fontsize=16, weight='bold')\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"Loss (MSE)\")\nhistory_df[['acc', 'val_acc']].plot(figsize=(12,5))\nplt.title(\"Accuracy\", fontsize=16, weight='bold')\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"% Accuracy\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Evaluation <a id=\"7\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load best weights according to MSE\nmodel.load_weights(SAVED_MODEL_NAME)\n#model.load_weights(WEIGHTS_PATH)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To evaluate our performance we predict values from the generator and round them of to the nearest integer to get valid predictions. After that we calculate the Quadratic Weighted Kappa score on the training set and the validation set."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Calculate QWK on train set\ny_train_preds, train_labels = get_preds_and_labels(model, train_generator)\ny_train_preds = np.rint(y_train_preds).astype(np.uint8).clip(0, 4)\n\n# Calculate score\ntrain_score = cohen_kappa_score(train_labels, y_train_preds, weights=\"quadratic\")\n\n# Calculate QWK on validation set\ny_val_preds, val_labels = get_preds_and_labels(model, val_generator)\ny_val_preds = np.rint(y_val_preds).astype(np.uint8).clip(0, 4)\n\n# Calculate score\nval_score = cohen_kappa_score(val_labels, y_val_preds, weights=\"quadratic\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f\"The Training Cohen Kappa Score is: {round(train_score, 5)}\")\nprint(f\"The Validation Cohen Kappa Score is: {round(val_score, 5)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can optimize the validation score by doing a [Grid Search](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html) over rounding thresholds instead of doing \"normal\" rounding. The \"OptimizedRounder\" class by [Abhishek Thakur](https://www.kaggle.com/abhishek) is a great way to do this. The original class can be found in [this Kaggle kernel](https://www.kaggle.com/abhishek/optimizer-for-quadratic-weighted-kappa)."},{"metadata":{"trusted":true},"cell_type":"code","source":"class OptimizedRounder(object):\n    \"\"\"\n    An optimizer for rounding thresholds\n    to maximize Quadratic Weighted Kappa score\n    \"\"\"\n    def __init__(self):\n        self.coef_ = 0\n\n    def _kappa_loss(self, coef, X, y):\n        \"\"\"\n        Get loss according to\n        using current coefficients\n        \"\"\"\n        X_p = np.copy(X)\n        for i, pred in enumerate(X_p):\n            if pred < coef[0]:\n                X_p[i] = 0\n            elif pred >= coef[0] and pred < coef[1]:\n                X_p[i] = 1\n            elif pred >= coef[1] and pred < coef[2]:\n                X_p[i] = 2\n            elif pred >= coef[2] and pred < coef[3]:\n                X_p[i] = 3\n            else:\n                X_p[i] = 4\n\n        ll = cohen_kappa_score(y, X_p, weights='quadratic')\n        return -ll\n\n    def fit(self, X, y):\n        \"\"\"\n        Optimize rounding thresholds\n        \"\"\"\n        loss_partial = partial(self._kappa_loss, X=X, y=y)\n        initial_coef = [0.5, 1.5, 2.5, 3.5]\n        self.coef_ = sp.optimize.minimize(loss_partial, initial_coef, method='nelder-mead')\n\n    def predict(self, X, coef):\n        \"\"\"\n        Make predictions with specified thresholds\n        \"\"\"\n        X_p = np.copy(X)\n        for i, pred in enumerate(X_p):\n            if pred < coef[0]:\n                X_p[i] = 0\n            elif pred >= coef[0] and pred < coef[1]:\n                X_p[i] = 1\n            elif pred >= coef[1] and pred < coef[2]:\n                X_p[i] = 2\n            elif pred >= coef[2] and pred < coef[3]:\n                X_p[i] = 3\n            else:\n                X_p[i] = 4\n        return X_p\n\n    def coefficients(self):\n        return self.coef_['x']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Optimize on validation data and evaluate again\ny_val_preds, val_labels = get_preds_and_labels(model, val_generator)\noptR = OptimizedRounder()\noptR.fit(y_val_preds, val_labels)\ncoefficients = optR.coefficients()\nopt_val_predictions = optR.predict(y_val_preds, coefficients)\nnew_val_score = cohen_kappa_score(val_labels, opt_val_predictions, weights=\"quadratic\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f\"Optimized Thresholds:\\n{coefficients}\\n\")\nprint(f\"The Validation Quadratic Weighted Kappa (QWK)\\n\\\nwith optimized rounding thresholds is: {round(new_val_score, 5)}\\n\")\nprint(f\"This is an improvement of {round(new_val_score - val_score, 5)}\\n\\\nover the unoptimized rounding\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission <a id=\"8\"></a>"},{"metadata":{},"cell_type":"markdown","source":"Since the test set is not that large we will not be using a generator for making the final predictions on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Place holder for diagnosis column\ntest_df['diagnosis'] = np.zeros(test_df.shape[0]) \n# For preprocessing test images, remember to change the preprocessing function to correct version \ntest_generator = ImageDataGenerator(preprocessing_function=preprocess_image, \n                                    rescale=1 / 128.).flow_from_dataframe(test_df, \n                                                                          x_col='id_code', \n                                                                          y_col='diagnosis',\n                                                                          directory=TEST_IMG_PATH,\n                                                                          target_size=(IMG_WIDTH, IMG_HEIGHT),\n                                                                          batch_size=BATCH_SIZE,\n                                                                          class_mode='other',\n                                                                          shuffle=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As mentioned earlier, we use custom thresholds to optimize our score. The same thresholds should be used when creating the final predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Make final predictions, round predictions and save to csv\ny_test, _ = get_preds_and_labels(model, test_generator)\ny_test = optR.predict(y_test, coefficients).astype(np.uint8)\ntest_df['diagnosis'] = y_test\n# Remove .png from ids\ntest_df['id_code'] = test_df['id_code'].str.replace(r'.png$', '')\ntest_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After creating the submission I always check the format and the distribution of the test predictions. Do they makes sense given the label distribution of the training data?"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Check submission\nprint(\"Submission File\")\ndisplay(test_df.head())","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Label distribution\ntrain_df['diagnosis'].value_counts().sort_index().plot(kind=\"bar\", \n                                                       figsize=(12,5), \n                                                       rot=0)\nplt.title(\"Label Distribution (Training Set)\", \n          weight='bold', \n          fontsize=18)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Label\", fontsize=17)\nplt.ylabel(\"Frequency\", fontsize=17);","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Distribution of predictions\ntest_df['diagnosis'].value_counts().sort_index().plot(kind=\"bar\", \n                                                      figsize=(12,5), \n                                                      rot=0)\nplt.title(\"Label Distribution (Predictions)\", \n          weight='bold', \n          fontsize=18)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Label\", fontsize=17)\nplt.ylabel(\"Frequency\", fontsize=17);","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Check kernels run-time. GPU limit for this competition is set to ± 9 hours.\nt_finish = time.time()\ntotal_time = round((t_finish-t_start) / 3600, 4)\nprint('Kernel runtime = {} hours ({} minutes)'.format(total_time, \n                                                      int(total_time*60)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"That's it! Feel free to experiment with this kernel and try a larger implementation of EfficientNet. This kernel provides weights for EfficientNetB0 through B5. Weights for EfficientNetB6 and B7 can be found in [Google AI's repository for EfficientNet](https://github.com/tensorflow/tpu/tree/master/models/official/efficientnet).\n\n**If you like this Kaggle kernel, feel free to give an upvote and leave a comment! I will try to implement your suggestions in this kernel!**"}],"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}