{"cells":[{"metadata":{},"cell_type":"markdown","source":"Thanks to Benjamin Warner for sharing the external dataset (2015 data)\n\nThe data can be found at: https://www.kaggle.com/benjaminwarner/resized-2015-2019-blindness-detection-images"},{"metadata":{},"cell_type":"markdown","source":"> 1 **. Imports:**\n\n> Import libraries that will be needed. Pytorch will be used as the main framework"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport torch\nimport cv2\n\nfrom torch import nn\nfrom torch import optim\nfrom torch.utils.data import Dataset, DataLoader, SubsetRandomSampler\nfrom torchvision import transforms, models\nfrom torch.optim import lr_scheduler\nfrom collections import OrderedDict\nfrom PIL import Image\nimport scipy as sp\nfrom functools import partial\nfrom sklearn import metrics\nfrom collections import Counter\nimport json\n\nfrom sklearn.model_selection import train_test_split\n\nimport seaborn as sns\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\ndevice = torch.device(\"cuda:0\")\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/resized-2015-2019-blindness-detection-images/labels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/aptos2019-blindness-detection","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Check if GPU is available for training:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_on_gpu = torch.cuda.is_available()\n\nif not train_on_gpu:\n    print('CUDA is not available.  Training on CPU ...')\nelse:\n    print('CUDA is available!  Training on GPU ...')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> 2 **. Exploratory Data Analysis**"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv(\"../input/aptos2019-blindness-detection/train.csv\")\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"type_percents = df[\"diagnosis\"].value_counts(normalize=True)\ntype_percents.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.barplot(x=type_percents.index, y=type_percents.values*100)\\\n          .set(xlabel=\"severity of diabetic retinopathy\", ylabel='Percent (%)')    \n\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Images are not well distributed among the severity of diabetic retinopathy."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(25, 25))\n# display 20 images\ntrain_imgs = os.listdir(\"../input/aptos2019-blindness-detection/train_images\")\nfor idx, img in enumerate(np.random.choice(train_imgs, 15)):\n    ax = fig.add_subplot(5, 15//5, idx+1, xticks=[], yticks=[])\n    im = Image.open(\"../input/aptos2019-blindness-detection/train_images/\" + img)\n    plt.imshow(im)\n    lab = df.loc[df['id_code'] == img.split('.')[0], 'diagnosis'].values[0]\n    ax.set_title(f'Label: {lab}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* We can see that the images differ in brightness intensity. It is probably a good idea to differ the brightness randomly to augment our training data later on, to increase robustness. \n\n* Looking at the images, we see a bright circle on the eye. This represents the optic disk with the main blood vessel portruding out of the optic disk.\n\nHow can we distinguish the severity based on these images? <br>\n* For stage 0, these are non contractors of Diabetic Retinopathy. According to Wikipedia (see [here](http://en.wikipedia.org/wiki/Diabetic_retinopathy)), subsequently, mild, moderate severe and proliferative. For non-proliferative, (mild and moderate in this case), microaneurysms can form. These are really small, to be precise, these have to be detected at the pixel level. For proliferative, abnormal new blood vessels (neovascularisation) form at the back of the eye. There may be a leakage of blood "},{"metadata":{},"cell_type":"markdown","source":"![](http://) 3 **. Load the data and perform preprocessing**[[](http://)](http://)\n\nThis involves:\n* Reading in the data\n* Preprocessing using Ben Graham's method and random cropping\n* Transforming the data\n\nSee Pytorch Documentation for Dataset: <br> https://pytorch.org/docs/stable/_modules/torch/utils/data/dataset.html#Dataset\n\nTutorial: <br>\nhttps://pytorch.org/tutorials/beginner/data_loading_tutorial.html\n\nPreprocessing kernel shared by Neuron Engineer:\nhttps://www.kaggle.com/ratthachat/aptos-simple-preprocessing-decoloring-cropping"},{"metadata":{"trusted":true},"cell_type":"code","source":"# def crop_image1(img,tol=7):\n#     # img is image data\n#     # tol  is tolerance\n        \n#     mask = img>tol\n#     return img[np.ix_(mask.any(1),mask.any(0))]\n\n\n# def crop_image(img,tol=7):\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#         h,w,_=img.shape\n# #         print(h,w)\n#         img1=cv2.resize(crop_image1(img[:,:,0]),(w,h))\n#         img2=cv2.resize(crop_image1(img[:,:,1]),(w,h))\n#         img3=cv2.resize(crop_image1(img[:,:,2]),(w,h))\n        \n# #         print(img1.shape,img2.shape,img3.shape)\n#         img[:,:,0]=img1\n#         img[:,:,1]=img2\n#         img[:,:,2]=img3\n#         return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AptosDrTrainDataset(Dataset):\n    def __init__(self, csv_file, transform=None):\n        self.data = pd.read_csv(csv_file)\n        self.transform = transform\n        self.targets = self.data.diagnosis\n        \n    def __len__(self):\n        return(len(self.data))\n    \n    def __getitem__(self, idx):\n        img_name = os.path.join(\"../input/aptos2019-blindness-detection/train_images\", \n                                self.data.loc[idx, 'id_code'] + '.png')\n#         IMG_SIZE = 224\n#         image = cv2.imread(img_name)\n#         image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)\n#         image = cv2.resize(image, (IMG_SIZE, IMG_SIZE))\n#         image = cv2.addWeighted ( image,4, cv2.GaussianBlur( image , (0,0) , IMG_SIZE/10) ,-4 ,128) # the trick is to add this line\n#         image = crop_image(image)\n        image = Image.open(img_name)\n        label = torch.tensor(self.data.loc[idx, 'diagnosis'])\n        if self.transform:\n            transformed_image = self.transform(image)\n            return transformed_image, label\n        \n        return image, label","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# train_transforms = transforms.Compose([transforms.ToPILImage(),\n#                                        transforms.ToTensor(),\n#                                        transforms.Normalize([0.485, 0.456, 0.406], \n#                                                             [0.229, 0.224, 0.225])])\n\n# train_transforms = transforms.Compose([transforms.ToPILImage(),\n#                                        transforms.Resize((224, 224)),\n#                                        transforms.RandomHorizontalFlip(),\n#                                        transforms.RandomVerticalFlip(),\n#                                        transforms.ToTensor(),\n#                                        transforms.Normalize([0.485, 0.456, 0.406], \n#                                                             [0.229, 0.224, 0.225])])\n\ntrain_transforms = transforms.Compose([transforms.Resize((224, 224)),\n                                       transforms.RandomHorizontalFlip(),\n                                       transforms.RandomVerticalFlip(),\n                                       transforms.ToTensor(),\n                                       transforms.Normalize([0.485, 0.456, 0.406], \n                                                            [0.229, 0.224, 0.225])])\n\ndataset = AptosDrTrainDataset(csv_file=\"../input/aptos2019-blindness-detection/train.csv\", transform=train_transforms)\n\nlabel = dataset.targets\n\ntrain_idx, valid_idx= train_test_split(\n    np.arange(len(label)), test_size=0.2, random_state=42, shuffle=True, stratify=label)\n\ntrain_sampler = SubsetRandomSampler(train_idx)\nvalid_sampler = SubsetRandomSampler(valid_idx)\n\ntrain_data_loader = torch.utils.data.DataLoader(dataset, batch_size=64, num_workers=0, \n                                                sampler=train_sampler)\n\nvalid_data_loader = torch.utils.data.DataLoader(dataset, batch_size=64, num_workers=0, \n                                                sampler=valid_sampler)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**![](http://) 4. Training a model** <br>\n\nThis involves:\n\n* Loading a pretrained model\n* Freezing most of the layers and retraining a few layers"},{"metadata":{"trusted":true},"cell_type":"code","source":"for image, label in train_data_loader:\n    print(image.shape)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = models.resnet152(pretrained=True)\nmodel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Replace last layer of model with our own"},{"metadata":{"trusted":true},"cell_type":"code","source":"model.fc = nn.Sequential(\n                          nn.BatchNorm1d(2048, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True),\n                          nn.Dropout(p=0.25),\n                          nn.Linear(in_features=2048, out_features=2048, bias=True),\n                          nn.ReLU(),\n                          nn.BatchNorm1d(2048, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True),\n                          nn.Dropout(p=0.5),\n                          nn.Linear(in_features=2048, out_features=1, bias=True),\n                         )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Find total parameters and trainable parameters\ntotal_params = sum(p.numel() for p in model.parameters())\nprint(f'{total_params:,} total parameters.')\ntotal_trainable_params = sum(\n    p.numel() for p in model.parameters() if p.requires_grad)\nprint(f'{total_trainable_params:,} training parameters.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.features[28]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Freeze all layers except for last few layers:"},{"metadata":{"trusted":true},"cell_type":"code","source":"for param in model.parameters():\n    param.requires_grad = False\n# for param in model.avgpool.parameters():\n#     param.requires_grad = True\n# for param in model.features[28].parameters():\n#      param.requires_grad = True\n# for param in model.features[29].parameters():\n#      param.requires_grad = True\n# for param in model.features[30].parameters():\n#      param.requires_grad = True\nfor param in model.avgpool.parameters():\n    param.requires_grad = True\nfor param in model.layer4[2].parameters():\n    param.requires_grad = True\nfor param in model.fc.parameters():\n    param.requires_grad = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Specifying loss and optimizer:"},{"metadata":{"trusted":true},"cell_type":"code","source":"criterion = nn.MSELoss()\noptimizer = optim.Adam(model.parameters(), lr=0.001)\nexp_lr_scheduler = lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Training Loop:"},{"metadata":{"trusted":true},"cell_type":"code","source":"if train_on_gpu:\n    model.cuda()\n\n# Number of epochs:\nn_epochs = 25\nvalid_loss_min = np.Inf\n  \nfor epoch in range(1, n_epochs + 1):\n  # keep track of training & validation loss\n  train_loss = 0.0\n  valid_loss = 0.0\n  \n  model.train()\n  exp_lr_scheduler.step()\n  \n  for data, target in train_data_loader:\n    # move tensors to GPU if CUDA is available\n    if train_on_gpu:\n        data, target = data.cuda(), target.cuda()\n    data = data.to(device, dtype=torch.float)\n    target = target.view(-1, 1)\n    target = target.to(device, dtype=torch.float)\n    # clear the gradients of all optimized variables\n    optimizer.zero_grad()\n    # forward pass: compute predicted outputs by passing inputs to the model\n    output = model(data)\n    output.cuda()\n    # calculate batch loss\n    loss = criterion(output, target)\n    # backward pass: compute gradient of loss wrt model parameters\n    loss.backward()\n    # perform single optimization step (parameter update)\n    optimizer.step()\n    # update training loss\n    train_loss += loss.item()*data.size(0)\n    \n   ### Validating the model ###\n  \n  model.eval()\n  \n#   class_correct = list(0. for i in range(5))\n#   class_total = list(0. for i in range(5))\n\n  for data, target in valid_data_loader:\n    # move tensors to GPU if CUDA is available\n    if train_on_gpu:\n        data, target = data.cuda(), target.cuda()\n    data = data.to(device, dtype=torch.float)\n    target = target.view(-1, 1)\n    target = target.to(device, dtype=torch.float)\n    # forward pass:\n    output = model(data)\n    # calculate batch loss\n    loss = criterion(output, target)\n    # update validation loss\n    valid_loss += loss.item()*data.size(0)\n    \n#     # convert output probabilities to predicted class\n#     _, pred = torch.max(output, 1)    \n#     # compare predictions to true label\n#     correct_tensor = pred.eq(target.data.view_as(pred))\n#     correct = np.squeeze(correct_tensor.numpy()) if not train_on_gpu else np.squeeze(correct_tensor.cpu().numpy())\n#     # calculate validation accuracy for each object class\n#     for i in range(target.data.size()[0]):\n#         label = target.data[i]\n#         class_correct[label] += correct[i].item()\n#         class_total[label] += 1\n    \n  # calculate average losses\n  train_loss = train_loss / len(train_data_loader.dataset)\n  valid_loss = valid_loss / len(valid_data_loader.dataset)\n  \n  # print training/validation statistics \n  print('Epoch: {} \\tTraining Loss: {:.6f} \\tValidation Loss: {:.6f}'.format(\n        epoch, train_loss, valid_loss))\n  \n#     print('\\nValidation Accuracy (Overall): %2d%% (%2d/%2d)' % (\n#     100. * np.sum(class_correct) / np.sum(class_total),\n#     np.sum(class_correct), np.sum(class_total)))\n    \n  # Save state_dict if validation loss decreased\n  if valid_loss <= valid_loss_min:\n    print('Validation loss decreased ({:.6f} --> {:.6f}).  Saving model ...'.format(\n        valid_loss_min,\n        valid_loss))\n    valid_loss_min = valid_loss\n    state = {\n    'epoch': epoch,\n    'state_dict': model.state_dict(),\n    'optimizer': optimizer.state_dict(),\n    'valid_loss_min': valid_loss_min\n    }\n    torch.save(state, 'aptos_resnet_152_9.pt')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load Model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"model.load_state_dict(torch.load(\"../input/aptos-resnet-152-8/aptos_resnet_152_8.pt\")[\"state_dict\"])\nmodel = model.to(device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Optimize using Quadratic Weighted Kappa:\n\nCredits goes to Abhishek for sharing"},{"metadata":{"trusted":true},"cell_type":"code","source":"for param in model.parameters():\n    param.requires_grad = False\n\nmodel.eval()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls /working","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AptosDrTestDataset(Dataset):\n    def __init__(self, csv_file, transform=None):\n        self.data = pd.read_csv(csv_file)\n        self.transform = transform\n        \n    def __len__(self):\n        return(len(self.data))\n    \n    def __getitem__(self, idx):\n        img_name = os.path.join('../input/resized-2015-2019-blindness-detection-images/resized test 15', \n                                self.data.loc[idx, 'image'] + '.jpg')\n        image = Image.open(img_name)\n        label = torch.tensor(self.data.loc[idx, 'level'])\n        if self.transform:\n            transformed_image = self.transform(image)\n            return transformed_image, label\n        return image, label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_transforms = transforms.Compose([transforms.RandomResizedCrop(224),\n                                       transforms.ToTensor(),\n                                       transforms.Normalize([0.485, 0.456, 0.406], \n                                                            [0.229, 0.224, 0.225])])\n\ntest_dataset = AptosDrTestDataset(csv_file='../input/resized-2015-2019-blindness-detection-images/labels/testLabels15.csv', transform=test_transforms)\ntest_data_loader = torch.utils.data.DataLoader(test_dataset, batch_size=32, shuffle=False, num_workers=4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = np.empty((0,1), int)\ntarget_lst = np.empty((0,1), int)\n\nfor data, target in valid_data_loader:\n    # move tensors to GPU if CUDA is available\n    if train_on_gpu:\n        data, target = data.cuda(), target.cuda()\n    # forward pass: compute predicted outputs by passing inputs to the model\n    output = model(data)\n    # convert output probabilities to predicted class\n#     _, pred = torch.max(output, 1) \n    preds = np.append(preds,  output.cpu().numpy(), axis=0)\n    target_ex =  np.expand_dims(target.cpu().numpy(), axis=1)\n    print(target_ex.shape)\n    target_lst = np.append(target_lst, target_ex, axis=0)\n    print(preds.shape)\n    print(target_lst.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class OptimizedRounder(object):\n    def __init__(self):\n        self.coef_ = 0\n\n    def _kappa_loss(self, coef, X, y):\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 = metrics.cohen_kappa_score(y, X_p, weights='quadratic')\n        return -ll\n\n    def fit(self, X, y):\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        print(-loss_partial(self.coef_['x']))\n\n    def predict(self, X, coef):\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":"preds.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds.T.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"target_lst.T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"optR = OptimizedRounder()\noptR.fit(preds, target_lst)\ncoefficients = optR.coefficients()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coefficients","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_df = pd.read_csv(\"../input/aptos2019-blindness-detection/sample_submission.csv\")\nsample_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coef = [0.5, 1.5, 2.5, 3.5]\n\nfor i, pred in enumerate(preds):\n    if pred < coef[0]:\n        preds[i] = 0\n    elif pred >= coef[0] and pred < coef[1]:\n        preds[i] = 1\n    elif pred >= coef[1] and pred < coef[2]:\n        preds[i] = 2\n    elif pred >= coef[2] and pred < coef[3]:\n        preds[i] = 3\n    else:\n        preds[i] = 4\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_df.diagnosis = preds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_df.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}