{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction\n\n## EfficientNet\n\nUsing the excellent pytorch-implementation from Luke Melas-Kyriazi:  \nhttps://github.com/lukemelas/EfficientNet-PyTorch\n\nEfficientNet-paper:  \nhttps://arxiv.org/abs/1905.11946\n\n## Pre-processing and data augmentation\n\nI use a pre-processing scheme similar to the one used by the [winner of the 2015 competition](http://blog.kaggle.com/2015/09/09/diabetic-retinopathy-winners-interview-1st-place-ben-graham).\n\nHowever, I made two changes:\n* added a center crop, such that the shape of the eye is not altered by resizing (during pre-processing)\n* increased the zoom-level of training images to match the zoom-level of test images (during data-augmentation)\n\nIn hindsight, that way probably not necessary but it's still interesting to think it through.\n\n**Motivation of avoiding shape-distorting resizing:**\n\nThe discrepancy between the local CV and the LB made me wary of a correlation between amount of black space around the eye and the diagnosis. (For instance, it would be conceivable that for more severe diagnoses, doctors take a closer look leading to a more zoomed-in image than for healthy eyes). If such a correlation existed in the training data but not in the public LB-data, that could explained the difference between validation- and LB-score.\n\nSince the black space around the eyes is mostly to the left and right of rectangular images, resizing those to square images alters the shape of the eyeball. If that correlates with the diagnosis, the model could learn features based on the shape as nicely explained in [this kernel](https://www.kaggle.com/taindow/be-careful-what-you-train-on).\n\nFor this reason, before resizing, I crop the largest possible square from the center of the image. This square image is then resized. This way, the shape of the eye is preserved exactly as in the original image.\n\n\n**Zoom difference between training and test-images:**\n \nA quick look at training- and test-images reveals that test-images tend to have a higher zoom which could also be a reason for the CV / LB - discrepancy. To remove this difference, I added a random-center-crop to the augmentation of the training-data but (but not the test-data). The random center-crop works just like the regular center-crop implemented in torchvision only randomly chooses the size of the cropped image within a user-specified range.\n\n## Runtime\n\nAbout half an hour."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%%capture\n\nimport warnings \nwarnings.filterwarnings(\"ignore\")\nimport os\nfrom os.path import join\nimport time\nfrom tqdm import tqdm\n\nimport numpy as np\nfrom numpy.random import choice\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.metrics import cohen_kappa_score\n\nimport PIL\nfrom PIL import Image\nimport cv2\n\nimport ipywidgets as widgets\nfrom ipywidgets import interactive\n\nimport torch\nfrom torch import nn\nfrom torch.nn import functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision import transforms, models as md\n\nimport sys\nsys.path.append('../input/efficientnet/efficientnet-pytorch/EfficientNet-PyTorch/')\nfrom efficientnet_pytorch import EfficientNet","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data pre-processing"},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = '../input/aptos2019-blindness-detection'\n\ntrain_dir = join(DATA_DIR, 'train_images')\nlabel_df  = pd.read_csv(join(DATA_DIR, 'train.csv'))\n\n\ndef train_validation_split(df, val_fraction=0.1):\n    val_ids  = np.random.choice(df.id_code, size=int(len(df) * val_fraction))\n    val_df   = df.query('id_code     in @val_ids')\n    train_df = df.query('id_code not in @val_ids')\n    return train_df, val_df\n\n\ntrain_df, val_df = train_validation_split(label_df)\nprint(train_df.shape, val_df.shape)\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_image_from_gray(img,tol=7):\n    \"\"\"\n    This function from:\n    https://www.kaggle.com/ratthachat/aptos-updatedv14-preprocessing-ben-s-cropping\n    \"\"\"\n    if img.ndim ==2:\n        mask = img>tol\n        return img[np.ix_(mask.any(1),mask.any(0))]\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    #         print(img1.shape,img2.shape,img3.shape)\n            img = np.stack([img1,img2,img3],axis=-1)\n    #         print(img.shape)\n        return img\n\n\ncv_to_pil = transforms.ToPILImage()\n\n    \ndef center_crop(image: PIL.Image):\n    \"\"\"\n    Only gets center square (of rectangular images) - no resizing\n    => diffently sized square images\n    \"\"\"\n    old_width, old_heigh = image.size\n    new_size = min(old_width, old_heigh)\n    \n    margin_x = (old_width - new_size) // 2\n    margin_y = (old_heigh - new_size) // 2\n    \n    left   = margin_x\n    right  = margin_x + new_size\n    top    = margin_y\n    bottom = margin_y + new_size\n    \n    return image.crop( (left, top, right, bottom) )\n\n\ndef process_image_ratio_invariant(cv2_image, size=256, do_center_crop=True):\n    \n    image = cv2.cvtColor(cv2_image, cv2.COLOR_BGR2RGB)\n    image = crop_image_from_gray(image)\n    #image = cv2.resize(image, (size, size))  # this would distort eyeball shape\n    \n    if do_center_crop is False:\n        return image\n    \n    # crop the largest possible square from the center\n    pil_img = cv_to_pil(image)\n    pil_img = center_crop(pil_img)\n    image   = np.array(pil_img).copy()\n    \n    # now we have quadratic, but differently sized images\n    # => resize without altering the shape of the eyeball\n    image = cv2.resize(image, (size, size))\n    \n    # add gaussian blur with sigma proportional to new size:\n    image = cv2.addWeighted (image, 4, cv2.GaussianBlur(image, (0, 0) , size/30) , -4 ,128)\n        \n    return cv_to_pil(image)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data + augmentation\n\n### Data augmentation\n\nBoth training-, and test-images are subjected to random horizontal flips and rotations. This inceases the diversity of the training data is hopefully decreases any systematic difference that might exist between training- and test-images.\n\nThe range of rotation between -20 and +20 degrees was chosen because the images show only a limited degree of rotation equivariance. (If an image is rotate by a few degrees, it looks just like another regular image. If it it's rotated by ~90 degrees it's doesn't look like a regular image any more).\n\n### In-memory dataset\n\nThe pre-processed data is stored in memory. While that takes a little while, it allows much faster training (less than a minute per epoch)."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n    \nclass Diabetic_Retionopathy_Data(Dataset):\n    \n    def __init__(self,\n                 image_dir: str,\n                 label_df: pd.DataFrame,\n                 train=True,\n                 transform=transforms.ToTensor(),\n                 sample_n=None,\n                 in_memory=False,\n                 write_images=False):\n        \"\"\"\n        @ image_dir:   path to directory with images\n        @ label_df:    df with image id (str) and label (0/1) - only for labeled test-set\n        @ transforms:  image transformation; by default no transformation\n        @ sample_n:    if not None, only use that many observations\n        \"\"\"\n        self.image_dir = image_dir\n        self.transform = transform\n        self.train     = train\n        self.in_memory = in_memory\n        \n        if sample_n:\n            label_df  = label_df.sample(n=min(sample_n, len(label_df)))\n            \n        ids            = set(label_df.id_code)\n        self.img_files = [f for f in os.listdir(image_dir) if f.split('.')[0] in ids]\n        label_df.index = label_df.id_code\n        self.label_df  = label_df.drop('id_code', axis=1)\n        \n        if in_memory:\n            \n            self.id2image = {}\n            for i, file_name in enumerate(self.img_files):\n                \n                if i and i % 500 == 0:\n                    print(f'{i} / {len(self.img_files)}')\n                \n                image = self._read_process_image(join(image_dir, file_name))\n                id_   = file_name.split('.')[0]\n                self.id2image[id_] = image\n                \n                if write_images:\n                    image.save(file_name)\n                    \n        print(f'Initialized datatset with {len(self.img_files)} images.\\n')\n        \n    @staticmethod\n    def _read_process_image(file_path: str, size=256):\n        image = cv2.imread(file_path)        \n        return process_image_ratio_invariant(image, size=size)        \n\n    def __getitem__(self, idx):\n\n        file_name = self.img_files[idx]\n        id_ = file_name.split('.')[0]\n        \n        if self.in_memory:\n            img = self.id2image[id_]\n        else:\n            img = self._read_process_image(join(self.image_dir, file_name))\n        \n        X   = self.transform(img)\n        \n        if self.train:\n            y = float(self.label_df.loc[id_].diagnosis)\n            return X, y, id_\n        else:\n            return X, id_\n    \n    def __len__(self):\n        return len(self.img_files)\n\n\nclass RandomCenterCrop(transforms.CenterCrop):\n    \"\"\"\n    Crops the PIL Image at the center.\n    :param: min_size, max_size: range of crop-size randomly within [min_size, max_size]\n    \"\"\"\n    def __init__(self, min_size: int, max_size: int):\n        self.min_size = min_size\n        self.max_size = max_size\n        \n    def __call__(self, img):\n        \"\"\"\n        Args:\n            img (PIL Image): Image to be cropped.\n        Returns:\n            PIL Image: Cropped image.\n        \"\"\"\n        size = np.random.randint(self.min_size, self.max_size + 1)\n        crop = transforms.CenterCrop( (size, size) )\n        return crop(img)\n\n    def __repr__(self):\n        return f'{self.__class__.__name__}: (min-size={self.min_size}, max-size={self.max_size})'\n\n\nbatchsize = 16\n\n# due to the large amount of data, random transformations might not be necessary...\ntrain_transform = transforms.Compose([\n    RandomCenterCrop(min_size=200, max_size=256),\n    transforms.Resize( (256, 256) ),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomRotation( (-20, 20) ),  \n    transforms.ToTensor(),\n    transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])\n])\n\ntrain = Diabetic_Retionopathy_Data(train_dir,\n                                   train_df,\n                                   transform=train_transform,\n                                   in_memory=True,\n                                   write_images=False)\nval   = Diabetic_Retionopathy_Data(train_dir,\n                                   val_df,\n                                   transform=train_transform,\n                                   in_memory=True,\n                                   write_images=False)\n\ntrain_loader = DataLoader(train, batch_size=batchsize, num_workers=4, shuffle=True)\nval_loader   = DataLoader(val,   batch_size=batchsize, num_workers=3, shuffle=False)\n\nX, y, _ = next(iter(val_loader))\nprint(f'batch-dimension:\\nX = {X.shape},\\ny = {y.shape}')\nprint(f'number of batches:\\ntrain: {len(train_loader)}\\nvalidation: {len(val_loader)}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Check pre-processing: compare raw vs. pre-processed images"},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_processed_images(image_dir, n=5, label_df=None, tf=None):\n    \n    sample_files = np.random.choice(os.listdir(image_dir), size=n)\n    \n    for file_name in sample_files:\n        \n        if label_df is not None:\n            id_ = file_name.split('.')[0]\n            diagnosis = label_df.query('id_code == @id_').diagnosis.item()\n        else:\n            diagnosis = 'unknown'\n        \n        image     = cv2.imread(join(image_dir, file_name))\n        raw_image = cv_to_pil(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))\n        \n        if tf is not None:\n            processed_image = tf(join(image_dir, file_name))\n        \n        fig, (ax1, ax2) = plt.subplots(1, 2)\n        fig.set_size_inches(10, 5)\n        \n        ax1.imshow(raw_image)\n        if tf is not None:\n            ax2.imshow(processed_image)\n        ax1.set_title('raw')\n        ax2.set_title('processed')\n        \n        fig.suptitle(f'diagnosis = {diagnosis}')            \n        plt.show()\n        \n    \nprint('TRAINING DATA:')\nshow_processed_images(join(DATA_DIR, 'train_images'),\n                      label_df=pd.concat([train_df, val_df]),\n                      tf=train._read_process_image)\nprint('TEST DATA:')\nshow_processed_images(join(DATA_DIR, 'test_images'),\n                      tf=train._read_process_image)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Define model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def count_parameters(model: nn.Module):\n    return sum([np.prod(x.shape) for x in model.parameters()])\n\n\ndef print_lr_schedule(lr: float, decay: float, num_epochs=20):\n    print('\\nlearning-rate schedule:')\n    for i in range(num_epochs):\n        if i % 2 == 0:\n            print(f'{i}\\t{lr:.6f}')\n        lr = lr* decay\n\n\nnet = EfficientNet.from_name('efficientnet-b0')\nnet.load_state_dict(torch.load('../input/efficientnet-pytorch/efficientnet-b0-08094119.pth'))\n\nnum_in_features = net._fc.in_features\nnet._fc = nn.Linear(num_in_features, 1)\n\nprint(f'number of parameters: {count_parameters(net)}')\nnet.train()\nnet.cuda()\n\nloss_function = nn.MSELoss()\nlr            = 0.0015\nlr_decay      = 0.97\n\noptimizer = torch.optim.Adam(net.parameters(), lr=lr, weight_decay=1e-5)\nscheduler = torch.optim.lr_scheduler.ExponentialLR(optimizer, gamma=lr_decay)\n\nprint_lr_schedule(lr, lr_decay)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train net"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nbest_epoch_score = np.inf\nprint('epoch\\ttrain-MSE\\tval-MSE\\tq-kappa\\tlr\\t\\ttime [min]')\nprint('------------------------------------------------------------------')\n\nfor epoch in range(25):\n    \n    start = time.time()\n    train_loss = []\n    \n    for i, (X, y, id_) in enumerate(train_loader):\n        \n        net.train()\n        optimizer.zero_grad()\n\n        out  = net(X.cuda())\n        loss = loss_function(out, y.float().cuda().view(-1, 1))    \n        \n        train_loss.append(loss.item())\n        \n        loss.backward()\n        optimizer.step()\n    \n    validation_loss = []\n    predictions     = np.array([])\n    truth           = np.array([])\n    \n    for  X, y, id_ in val_loader:\n\n        net.eval()\n\n        out  = net(X.cuda())\n        loss = loss_function(out, y.float().cuda().view(-1, 1))\n        \n        validation_loss.append(loss.item())\n        predictions = np.append(predictions, out.detach().cpu().numpy())\n        truth       = np.append(truth, y.detach().cpu().numpy().astype(int))\n\n    current_lr = optimizer.param_groups[0]['lr']\n    scheduler.step()\n    qk = cohen_kappa_score(predictions.round().astype(int), truth, weights='quadratic')\n    duration = (time.time() - start) / 60\n    print(f'{epoch}:\\t{np.mean(train_loss):.4f}\\t\\t{np.mean(validation_loss):.4f}\\t{qk:.4f}\\t{current_lr:.6f}\\t{duration:.2f}')\n    \n    if np.mean(validation_loss) < best_epoch_score:\n        torch.save(net.state_dict(), 'state_dict_best.pt')\n        best_epoch_score = np.mean(validation_loss)\n        best_epoch = epoch\n        \n\nprint(f'epoch with best validation-score: {best_epoch}')\n\nplt.hist(predictions, bins=5)\nplt.xlim(-1, 5)\nplt.title('distribution of predictions\\n(before rounding)')\nplt.show()\n\nplt.hist(train_df.diagnosis.values, bins=5)\nplt.xlim(-1, 5)\nplt.title('distribution of labels')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Make test-set predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_dir = join(DATA_DIR, 'test_images')\ntest_df  = pd.read_csv(join(DATA_DIR, 'test.csv'))\ntest_df.head(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Compare test- and training-images\n\nComparing test- and training-images side by side shows that, on average, test-images tend to have a higher zoom-level. The RandomCenterCrop applied to training-images only (as explained and implemented above) addresses this issue."},{"metadata":{"trusted":true},"cell_type":"code","source":"def sample_images(train_dir: str, test_dir: str, n=10):\n    \"\"\"\n    Show n train- and test-images side by side.\n    \"\"\"\n    train_files = choice(os.listdir(train_dir), size=n)\n    test_files  = choice(os.listdir(test_dir),  size=n)\n    images      = []\n    \n    for train_f, test_f in zip(train_files, test_files):\n        train_img = Image.open(join(train_dir, train_f))\n        test_img  = Image.open(join(test_dir,  test_f)) \n        images.append( (train_img, test_img) )\n        \n    def show_image(i):\n        train_image, test_image = images[i]\n        \n        fig, (ax1, ax2) = plt.subplots(1, 2)\n        fig.set_size_inches(20, 10)\n        ax1.imshow(train_image)\n        ax2.imshow(test_image)\n        ax1.set_title('train')\n        ax2.set_title('test')\n        plt.show()\n        \n    return interactive(show_image, i=(0, n-1))\n        \nsample_images(train_dir=train_dir,\n              test_dir=test_dir,\n              n=10)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"** Test-time augmentation **\n\nThe test-data is augmented in exactly the same way as the training data except for the RandomCenterCrop."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_transform = transforms.Compose([\n    transforms.Resize( (256, 256) ),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomRotation( (-20, 20) ),\n    transforms.ToTensor(),\n    transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])\n])\n\ntest_ds = Diabetic_Retionopathy_Data(test_dir,\n                                     test_df,\n                                     transform=test_transform,\n                                     train=False)\n\ntest_loader = DataLoader(test_ds, batch_size=batchsize)\n\nX, _ = next(iter(test_loader))\nprint(f'batch-dimension:\\nX = {X.shape}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Make predictions**"},{"metadata":{"trusted":true},"cell_type":"code","source":"net.load_state_dict(torch.load('state_dict_best.pt'))\nnet.eval()\nnet.cuda()\nid2prediction = {}\n\nfor i, (X, id_) in enumerate(test_loader):\n    out           = net(X.cuda())\n    preds         = out.detach().cpu().numpy().ravel()\n    id2prediction = {**id2prediction, **dict(zip( id_, preds.round().astype(int).tolist() ))}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df = pd.read_csv(join(DATA_DIR, 'sample_submission.csv'))\nsubmission_df.diagnosis = submission_df.id_code.map(id2prediction)\n\n# limit predictions to interval [0, 4] !!\nsubmission_df.diagnosis = submission_df.diagnosis.map(lambda p: max(p, 0))\nsubmission_df.diagnosis = submission_df.diagnosis.map(lambda p: min(p, 4))\n\nsubmission_df.to_csv('submission.csv', index=False)\ndisplay(submission_df.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Check the predictions**\n\nNote that the test-data seemst to have a different diagnosis-distribution than the training data. Or maybe the predictions all very off. Or maybe both :-)"},{"metadata":{"trusted":true},"cell_type":"code","source":"assert len(submission_df) == len(submission_df.dropna())\nassert set(submission_df.diagnosis) == {0, 1, 2, 3, 4}\nsubmission_df.diagnosis.hist()\nplt.title('distribution of test-set predictions');","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}