{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport json\nimport csv\nimport random\nimport pickle\nimport cv2\nimport numpy as np\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nimport torchvision.transforms as transforms\n\nfrom PIL import Image\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy.ndimage.measurements import label\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import roc_auc_score, roc_curve\n","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:06.143886Z","iopub.execute_input":"2021-05-23T09:04:06.14423Z","iopub.status.idle":"2021-05-23T09:04:08.909471Z","shell.execute_reply.started":"2021-05-23T09:04:06.144185Z","shell.execute_reply":"2021-05-23T09:04:08.908665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset class","metadata":{}},{"cell_type":"code","source":"class RefugeDataset(Dataset):\n\n def __init__(self, root_dir, split='train', output_size=(256,256)):\n # Define attributes\n self.output_size = output_size\n self.root_dir = root_dir\n self.split = split\n \n # Load data index\n with open(os.path.join(self.root_dir, self.split, 'index.json')) as f:\n self.index = json.load(f)\n \n self.images = []\n for k in range(len(self.index)):\n print('Loading {} image {}/{}...'.format(split, k, len(self.index)), end='\\r')\n img_name = os.path.join(self.root_dir, self.split, 'images', self.index[str(k)]['ImgName'])\n img = np.array(Image.open(img_name).convert('RGB'))\n img = transforms.functional.to_tensor(img)\n img = transforms.functional.resize(img, self.output_size, interpolation=Image.BILINEAR)\n self.images.append(img)\n \n # Load ground truth for 'train' and 'val' sets\n if split != 'test':\n self.segs = []\n for k in range(len(self.index)):\n print('Loading {} segmentation {}/{}...'.format(split, k, len(self.index)), end='\\r')\n seg_name = os.path.join(self.root_dir, self.split, 'gts', self.index[str(k)]['ImgName'].split('.')[0]+'.bmp')\n seg = np.array(Image.open(seg_name)).copy()\n seg = 255. - seg\n od = (seg>=127.).astype(np.float32)\n oc = (seg>=250.).astype(np.float32)\n od = torch.from_numpy(od[None,:,:])\n oc = torch.from_numpy(oc[None,:,:])\n od = transforms.functional.resize(od, self.output_size, interpolation=Image.NEAREST)\n oc = transforms.functional.resize(oc, self.output_size, interpolation=Image.NEAREST)\n seg = torch.cat([od, oc], dim=0)\n self.segs.append(seg)\n \n print('Succesfully loaded {} dataset.'.format(split) + ' '*50)\n \n \n def __len__(self):\n return len(self.index)\n\n def __getitem__(self, idx):\n # Image\n img = self.images[idx]\n \n # Return only images for 'test' set\n if self.split == 'test':\n return img\n \n # Else, images and ground truth\n else:\n # Label\n lab = torch.tensor(self.index[str(idx)]['Label'], dtype=torch.float32)\n\n # Segmentation masks\n seg = self.segs[idx]\n\n # Fovea localization\n f_x = self.index[str(idx)]['Fovea_X']\n f_y = self.index[str(idx)]['Fovea_Y']\n fov = torch.FloatTensor([f_x, f_y])\n \n return img, lab, seg, fov, self.index[str(idx)]['ImgName']","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:08.911041Z","iopub.execute_input":"2021-05-23T09:04:08.911361Z","iopub.status.idle":"2021-05-23T09:04:08.930079Z","shell.execute_reply.started":"2021-05-23T09:04:08.911334Z","shell.execute_reply":"2021-05-23T09:04:08.929069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Metrics","metadata":{}},{"cell_type":"code","source":"EPS = 1e-7\n\ndef compute_dice_coef(input, target):\n '''\n Compute dice score metric.\n '''\n batch_size = input.shape[0]\n return sum([dice_coef_sample(input[k,:,:], target[k,:,:]) for k in range(batch_size)])/batch_size\n\ndef dice_coef_sample(input, target):\n iflat = input.contiguous().view(-1)\n tflat = target.contiguous().view(-1)\n intersection = (iflat * tflat).sum()\n return (2. * intersection) / (iflat.sum() + tflat.sum())\n\n\ndef vertical_diameter(binary_segmentation):\n '''\n Get the vertical diameter from a binary segmentation.\n The vertical diameter is defined as the \"fattest\" area of the binary_segmentation parameter.\n '''\n\n # get the sum of the pixels in the vertical axis\n vertical_axis_diameter = np.sum(binary_segmentation, axis=1)\n\n # pick the maximum value\n diameter = np.max(vertical_axis_diameter, axis=1)\n\n # return it\n return diameter\n\n\n\ndef vertical_cup_to_disc_ratio(od, oc):\n '''\n Compute the vertical cup-to-disc ratio from a given labelling map.\n '''\n # compute the cup diameter\n cup_diameter = vertical_diameter(oc)\n # compute the disc diameter\n disc_diameter = vertical_diameter(od)\n\n return cup_diameter / (disc_diameter + EPS)\n\ndef compute_vCDR_error(pred_od, pred_oc, gt_od, gt_oc):\n '''\n Compute vCDR prediction error, along with predicted vCDR and ground truth vCDR.\n '''\n pred_vCDR = vertical_cup_to_disc_ratio(pred_od, pred_oc)\n gt_vCDR = vertical_cup_to_disc_ratio(gt_od, gt_oc)\n vCDR_err = np.mean(np.abs(gt_vCDR - pred_vCDR))\n return vCDR_err, pred_vCDR, gt_vCDR\n\n\ndef classif_eval(classif_preds, classif_gts):\n '''\n Compute AUC classification score.\n '''\n auc = roc_auc_score(classif_gts, classif_preds)\n return auc\n\n\ndef fov_error(pred_fov, gt_fov):\n '''\n Fovea localization error metric (mean root squared error).\n '''\n err = np.sqrt(np.sum((gt_fov-pred_fov)**2, axis=1)).mean()\n return err","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:15.981976Z","iopub.execute_input":"2021-05-23T09:04:15.982284Z","iopub.status.idle":"2021-05-23T09:04:15.995185Z","shell.execute_reply.started":"2021-05-23T09:04:15.982254Z","shell.execute_reply":"2021-05-23T09:04:15.994326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Post-processing functions","metadata":{}},{"cell_type":"code","source":"def refine_seg(pred):\n '''\n Only retain the biggest connected component of a segmentation map.\n '''\n np_pred = pred.numpy()\n \n largest_ccs = []\n for i in range(np_pred.shape[0]):\n labeled, ncomponents = label(np_pred[i,:,:])\n bincounts = np.bincount(labeled.flat)[1:]\n if len(bincounts) == 0:\n largest_cc = labeled == 0\n else:\n largest_cc = labeled == np.argmax(bincounts)+1\n largest_cc = torch.tensor(largest_cc, dtype=torch.float32)\n largest_ccs.append(largest_cc)\n largest_ccs = torch.stack(largest_ccs)\n \n return largest_ccs","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:19.485344Z","iopub.execute_input":"2021-05-23T09:04:19.485693Z","iopub.status.idle":"2021-05-23T09:04:19.492713Z","shell.execute_reply.started":"2021-05-23T09:04:19.485663Z","shell.execute_reply":"2021-05-23T09:04:19.491797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Network","metadata":{}},{"cell_type":"code","source":"class UNet(nn.Module):\n def __init__(self, n_channels=3, n_classes=2):\n super(UNet, self).__init__()\n self.n_channels = n_channels\n self.n_classes = n_classes\n self.epoch = 0\n\n self.inc = DoubleConv(n_channels, 64)\n self.down1 = Down(64, 128)\n self.down2 = Down(128, 256)\n self.down3 = Down(256, 512)\n factor = 2 \n self.down4 = Down(512, 1024 // factor)\n self.up1 = Up(1024, 512 // factor)\n self.up2 = Up(512, 256 // factor)\n self.up3 = Up(256, 128 // factor)\n self.up4 = Up(128, 64)\n self.output_layer = OutConv(64, n_classes)\n\n def forward(self, x):\n x1 = self.inc(x)\n x2 = self.down1(x1)\n x3 = self.down2(x2)\n x4 = self.down3(x3)\n x5 = self.down4(x4)\n out = self.up1(x5, x4)\n out = self.up2(out, x3)\n out = self.up3(out, x2)\n out = self.up4(out, x1)\n out = self.output_layer(out)\n out = torch.sigmoid(out)\n return out\n\n \nclass DoubleConv(nn.Module):\n \"\"\"(convolution => [BN] => ReLU) * 2\"\"\"\n\n def __init__(self, in_channels, out_channels, mid_channels=None):\n super().__init__()\n if not mid_channels:\n mid_channels = out_channels\n self.double_conv = nn.Sequential(\n nn.Conv2d(in_channels, mid_channels, kernel_size=3, padding=1),\n nn.BatchNorm2d(mid_channels),\n nn.ReLU(inplace=True),\n nn.Conv2d(mid_channels, out_channels, kernel_size=3, padding=1),\n nn.BatchNorm2d(out_channels),\n nn.ReLU(inplace=True)\n )\n\n def forward(self, x):\n return self.double_conv(x)\n\n\nclass Down(nn.Module):\n \"\"\"Downscaling with maxpool then double conv\"\"\"\n\n def __init__(self, in_channels, out_channels):\n super().__init__()\n self.maxpool_conv = nn.Sequential(\n nn.MaxPool2d(2),\n DoubleConv(in_channels, out_channels)\n )\n\n def forward(self, x):\n return self.maxpool_conv(x)\n\n\nclass Up(nn.Module):\n \"\"\"Upscaling then double conv\"\"\"\n\n def __init__(self, in_channels, out_channels):\n super().__init__()\n\n # Use the normal convolutions to reduce the number of channels\n self.up = nn.Upsample(scale_factor=2, mode='bilinear', align_corners=True)\n self.conv = DoubleConv(in_channels, out_channels, in_channels // 2)\n\n\n def forward(self, x1, x2):\n x1 = self.up(x1)\n # input is CHW\n diffY = x2.size()[2] - x1.size()[2]\n diffX = x2.size()[3] - x1.size()[3]\n\n x1 = F.pad(x1, [diffX // 2, diffX - diffX // 2,\n diffY // 2, diffY - diffY // 2])\n x = torch.cat([x2, x1], dim=1)\n return self.conv(x)\n\n\nclass OutConv(nn.Module):\n '''\n Simple convolution.\n '''\n def __init__(self, in_channels, out_channels):\n super(OutConv, self).__init__()\n self.conv = nn.Conv2d(in_channels, out_channels, kernel_size=1)\n\n def forward(self, x):\n return self.conv(x)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:21.890756Z","iopub.execute_input":"2021-05-23T09:04:21.891174Z","iopub.status.idle":"2021-05-23T09:04:21.913818Z","shell.execute_reply.started":"2021-05-23T09:04:21.891134Z","shell.execute_reply":"2021-05-23T09:04:21.912881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Settings","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 8\nnum_workers = 8\ntotal_epoch = 100","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:25.238906Z","iopub.execute_input":"2021-05-23T09:04:25.239232Z","iopub.status.idle":"2021-05-23T09:04:25.245654Z","shell.execute_reply.started":"2021-05-23T09:04:25.239201Z","shell.execute_reply":"2021-05-23T09:04:25.244782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create datasets and data loaders\nAll image files are loaded in RAM in order to speed up the pipeline. Therefore, each dataset creation should take a few minutes.","metadata":{}},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset(root_dir, \n split='train')\nval_set = RefugeDataset(root_dir, \n split='val')\ntest_set = RefugeDataset(root_dir, \n split='test')\n\n# Dataloaders\ntrain_loader = DataLoader(train_set, \n batch_size=batch_size, \n shuffle=True, \n num_workers=num_workers,\n pin_memory=True,\n )\nval_loader = DataLoader(val_set, \n batch_size=batch_size, \n shuffle=False, \n num_workers=num_workers,\n pin_memory=True,\n )\ntest_loader = DataLoader(test_set, \n batch_size=batch_size, \n shuffle=False, \n num_workers=num_workers,\n pin_memory=True)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:04:27.462329Z","iopub.execute_input":"2021-05-23T09:04:27.462678Z","iopub.status.idle":"2021-05-23T09:09:14.503224Z","shell.execute_reply.started":"2021-05-23T09:04:27.462649Z","shell.execute_reply":"2021-05-23T09:09:14.502243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Observation\n### Observing segmentations","metadata":{}},{"cell_type":"code","source":"image, label, mask, fov, img_name = train_set.__getitem__(0)\n\nR_img = image[0,] \nG_img = image[1,]\nB_img = image[2,]\n\nfig, axs = plt.subplots(2, 1, figsize=(16, 16))\n\naxs[0].imshow(image[0,])\naxs[1].imshow(mask[0,].squeeze())\naxs[1].imshow(mask[1,].squeeze(), alpha=0.5)\n\naxs[0].title.set_text(\"Retinal fundus image\")\naxs[1].title.set_text(\"Its associated mask\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Device, model, loss and optimizer","metadata":{}},{"cell_type":"code","source":"# Device\ndevice = torch.device(\"cuda:0\")\n\n# Network\nmodel = UNet(n_channels=3, n_classes=2).to(device)\n\n# Loss\nseg_loss = torch.nn.BCELoss(reduction='mean')","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:09:14.505002Z","iopub.execute_input":"2021-05-23T09:09:14.50554Z","iopub.status.idle":"2021-05-23T09:09:18.931985Z","shell.execute_reply.started":"2021-05-23T09:09:14.505498Z","shell.execute_reply":"2021-05-23T09:09:18.931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train for OC/OD segmentation","metadata":{}},{"cell_type":"code","source":"def train_model(model, lr, total_epoch):\n # Define parameters\n nb_train_batches = len(train_loader)\n nb_val_batches = len(val_loader)\n nb_iter = 0\n best_val_auc = 0.\n \n # Optimizer\n optimizer = optim.Adam(model.parameters(), lr=lr)\n while model.epoch < total_epoch:\n # Accumulators\n train_vCDRs, val_vCDRs = [], []\n train_classif_gts, val_classif_gts = [], []\n train_loss, val_loss = 0., 0.\n train_dsc_od, val_dsc_od = 0., 0.\n train_dsc_oc, val_dsc_oc = 0., 0.\n train_vCDR_error, val_vCDR_error = 0., 0.\n\n ############\n # TRAINING #\n ############\n model.train()\n train_data = iter(train_loader)\n for k in range(nb_train_batches):\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n # Forward pass\n logits = model(imgs)\n loss = seg_loss(logits, seg_gts)\n\n # Backward pass\n optimizer.zero_grad()\n loss.backward()\n optimizer.step()\n train_loss += loss.item() / nb_train_batches\n\n with torch.no_grad():\n # Compute segmentation metric\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n train_dsc_od += dsc_od.item()/nb_train_batches\n train_dsc_oc += dsc_oc.item()/nb_train_batches\n\n\n # Compute and store vCDRs\n vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n train_vCDRs += pred_vCDR.tolist()\n train_vCDR_error += vCDR_error / nb_train_batches\n train_classif_gts += classif_gts.cpu().numpy().tolist()\n\n # Increase iterations\n nb_iter += 1\n\n # Std out\n print('Epoch {}, iter {}/{}, loss {:.6f}'.format(model.epoch+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n end='\\r')\n\n # Train a logistic regression on vCDRs\n train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n train_classif_gts = np.array(train_classif_gts)\n clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\n train_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n train_auc = classif_eval(train_classif_preds, train_classif_gts)\n\n ##############\n # VALIDATION #\n ##############\n model.eval()\n with torch.no_grad():\n val_data = iter(val_loader)\n for k in range(nb_val_batches):\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n # Forward pass\n logits = model(imgs)\n val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n # Std out\n print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n end='\\r')\n\n # Compute segmentation metric\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n val_dsc_od += dsc_od.item()/nb_val_batches\n val_dsc_oc += dsc_oc.item()/nb_val_batches\n\n # Compute and store vCDRs\n vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n val_vCDRs += pred_vCDR.tolist()\n val_vCDR_error += vCDR_error / nb_val_batches\n val_classif_gts += classif_gts.cpu().numpy().tolist()\n\n\n # Glaucoma predictions from vCDRs\n val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n val_classif_gts = np.array(val_classif_gts)\n val_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\n val_auc = classif_eval(val_classif_preds, val_classif_gts)\n\n # Validation results\n print('VALIDATION epoch {}'.format(model.epoch+1)+' '*50)\n print('LOSSES: {:.4f} (train), {:.4f} (val)'.format(train_loss, val_loss))\n print('OD segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_od, val_dsc_od))\n print('OC segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_oc, val_dsc_oc))\n print('vCDR error: {:.4f} (train), {:.4f} (val)'.format(train_vCDR_error, val_vCDR_error))\n print('Classification (AUC): {:.4f} (train), {:.4f} (val)'.format(train_auc, val_auc))\n\n # Save model if best validation AUC is reached\n if val_auc > best_val_auc:\n torch.save(model.state_dict(), '/kaggle/working/best_AUC_weights.pth')\n with open('/kaggle/working/best_AUC_classifier.pkl', 'wb') as clf_file:\n pickle.dump(clf, clf_file)\n best_val_auc = val_auc\n print('Best validation AUC reached. Saved model weights and classifier.')\n print('_'*50)\n\n # End of epoch\n model.epoch += 1\n \n return best_val_auc","metadata":{"execution":{"iopub.status.busy":"2021-05-23T09:09:18.936655Z","iopub.execute_input":"2021-05-23T09:09:18.937002Z","iopub.status.idle":"2021-05-23T09:09:18.979909Z","shell.execute_reply.started":"2021-05-23T09:09:18.936967Z","shell.execute_reply":"2021-05-23T09:09:18.978937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MODEL 1: REPLACING THE LOGISTIC REGRESSION WITH A SVM or a XGBOOST","metadata":{}},{"cell_type":"code","source":"# Define parameters\nnb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)\nnb_iter = 0\nbest_val_auc = 0.\n\nwhile model.epoch < total_epoch:\n # Accumulators\n train_vCDRs, val_vCDRs = [], []\n train_classif_gts, val_classif_gts = [], []\n train_loss, val_loss = 0., 0.\n train_dsc_od, val_dsc_od = 0., 0.\n train_dsc_oc, val_dsc_oc = 0., 0.\n train_vCDR_error, val_vCDR_error = 0., 0.\n \n ############\n # TRAINING #\n ############\n model.train()\n train_data = iter(train_loader)\n for k in range(nb_train_batches):\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n # Forward pass\n logits = model(imgs)\n loss = seg_loss(logits, seg_gts)\n \n # Backward pass\n optimizer.zero_grad()\n loss.backward()\n optimizer.step()\n train_loss += loss.item() / nb_train_batches\n \n with torch.no_grad():\n # Compute segmentation metric\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n train_dsc_od += dsc_od.item()/nb_train_batches\n train_dsc_oc += dsc_oc.item()/nb_train_batches\n\n\n # Compute and store vCDRs\n vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n train_vCDRs += pred_vCDR.tolist()\n train_vCDR_error += vCDR_error / nb_train_batches\n train_classif_gts += classif_gts.cpu().numpy().tolist()\n \n # Increase iterations\n nb_iter += 1\n \n # Std out\n print('Epoch {}, iter {}/{}, loss {:.6f}'.format(model.epoch+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n end='\\r')\n \n # Train a logistic regression on vCDRs\n train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n train_classif_gts = np.array(train_classif_gts)\n #clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\n # XGBOOST for classification\n # Using early stopping. If the predictions don't improve (using AUC as metric), we build 10 more trees.\n # If none of the 10 trees improve the predictions, we stop.\n # clf = xgb.XGBClassifier(objective='binary:logistic', random_state=1).fit(train_vCDRs, \n # train_classif_gts, \n # verbose=False,\n # early_stopping_rounds=10,\n # eval_metric='aucpr',\n # eval_set=[(val_vCDRs, val_classif_gts)])\n clf = SVC(1, 'poly', 2, gamma=4.0, probability=True).fit(train_vCDRs, train_classif_gts)\n train_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n train_auc = classif_eval(train_classif_preds, train_classif_gts)\n \n ##############\n # VALIDATION #\n ##############\n model.eval()\n with torch.no_grad():\n val_data = iter(val_loader)\n for k in range(nb_val_batches):\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n # Forward pass\n logits = model(imgs)\n val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n # Std out\n print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n end='\\r')\n \n # Compute segmentation metric\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n val_dsc_od += dsc_od.item()/nb_val_batches\n val_dsc_oc += dsc_oc.item()/nb_val_batches\n \n # Compute and store vCDRs\n vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n val_vCDRs += pred_vCDR.tolist()\n val_vCDR_error += vCDR_error / nb_val_batches\n val_classif_gts += classif_gts.cpu().numpy().tolist()\n \n\n # Glaucoma predictions from vCDRs\n val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n val_classif_gts = np.array(val_classif_gts)\n val_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\n val_auc = classif_eval(val_classif_preds, val_classif_gts)\n \n # Validation results\n print('VALIDATION epoch {}'.format(model.epoch+1)+' '*50)\n print('LOSSES: {:.4f} (train), {:.4f} (val)'.format(train_loss, val_loss))\n print('OD segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_od, val_dsc_od))\n print('OC segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_oc, val_dsc_oc))\n print('vCDR error: {:.4f} (train), {:.4f} (val)'.format(train_vCDR_error, val_vCDR_error))\n print('Classification (AUC): {:.4f} (train), {:.4f} (val)'.format(train_auc, val_auc))\n \n # Save model if best validation AUC is reached\n if val_auc > best_val_auc:\n torch.save(model.state_dict(), '/kaggle/working/best_AUC_weights.pth')\n with open('/kaggle/working/best_AUC_classifier.pkl', 'wb') as clf_file:\n pickle.dump(clf, clf_file)\n best_val_auc = val_auc\n print('Best validation AUC reached. Saved model weights and classifier.')\n print('_'*50)\n \n # End of epoch\n model.epoch += 1\n \n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Gridsearch","metadata":{}},{"cell_type":"markdown","source":"

To find the best hyperparameters, we do a gridSearch on our model, evaluating the parameters on the best AUC

","metadata":{}},{"cell_type":"code","source":"param_dict = {'lr':[1e-3, 1e-4, 1e-5], 'total_epoch': [100, 150, 200]}","metadata":{"execution":{"iopub.status.busy":"2021-05-22T15:58:30.814295Z","iopub.execute_input":"2021-05-22T15:58:30.814611Z","iopub.status.idle":"2021-05-22T15:58:30.818817Z","shell.execute_reply.started":"2021-05-22T15:58:30.814582Z","shell.execute_reply":"2021-05-22T15:58:30.817929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# custom gridsearch\nresults = np.zeros([3,3])\ni = 0\nj = 0\nfor learning_rate in param_dict['lr']:\n for total_epoch in param_dict['total_epoch']:\n print('LEARNING RATE: {}, TOTAL_EPOCH: {}'.format(learning_rate, total_epoch))\n # Device\n device = torch.device(\"cuda:0\")\n # Network\n model = UNet(n_channels=3, n_classes=2).to(device)\n # Loss\n seg_loss = torch.nn.BCELoss(reduction='mean')\n \n best_auc = train_model(model, lr, total_epoch)\n results[i][j] = best_auc\n j+=1\n j = 0\n i+=1\nresults","metadata":{"execution":{"iopub.status.busy":"2021-05-22T15:58:32.555384Z","iopub.execute_input":"2021-05-22T15:58:32.555781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load best model + classifier","metadata":{}},{"cell_type":"code","source":"# Load model and classifier\nmodel = UNet(n_channels=3, n_classes=2).to(device)\nmodel.load_state_dict(torch.load('/kaggle/working/best_AUC_weights.pth'))\nwith open('/kaggle/working/best_AUC_classifier.pkl', 'rb') as clf_file:\n clf = pickle.load(clf_file)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:00:35.005026Z","iopub.execute_input":"2021-05-23T14:00:35.00537Z","iopub.status.idle":"2021-05-23T14:00:35.209897Z","shell.execute_reply.started":"2021-05-23T14:00:35.005337Z","shell.execute_reply":"2021-05-23T14:00:35.209065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Check performance is maintained on validation","metadata":{}},{"cell_type":"code","source":"model.eval()\nval_vCDRs = []\nval_classif_gts = []\nval_loss = 0.\nval_dsc_od = 0.\nval_dsc_oc = 0.\nnb_val_batches = len(val_loader)\nval_vCDR_error = 0.\nwith torch.no_grad():\n val_data = iter(val_loader)\n for k in range(nb_val_batches):\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n # Forward pass\n logits = model(imgs)\n val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n # Std out\n print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n end='\\r')\n\n # Compute segmentation metric\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n val_dsc_od += dsc_od.item()/nb_val_batches\n val_dsc_oc += dsc_oc.item()/nb_val_batches\n\n # Compute and store vCDRs\n vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n val_vCDRs += pred_vCDR.tolist()\n val_vCDR_error += vCDR_error / nb_val_batches\n val_classif_gts += classif_gts.cpu().numpy().tolist()\n\n\n# Glaucoma predictions from vCDRs\nval_vCDRs = np.array(val_vCDRs).reshape(-1,1)\nval_classif_gts = np.array(val_classif_gts)\nval_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\nval_auc = classif_eval(val_classif_preds, val_classif_gts)\n\n# Validation results\nprint('VALIDATION '+' '*50)\nprint('LOSSES: {:.4f} (val)'.format(val_loss))\nprint('OD segmentation (Dice Score): {:.4f} (val)'.format(val_dsc_od))\nprint('OC segmentation (Dice Score): {:.4f} (val)'.format(val_dsc_oc))\nprint('vCDR error: {:.4f} (val)'.format(val_vCDR_error))\nprint('Classification (AUC): {:.4f} (val)'.format(val_auc))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:00:37.823208Z","iopub.execute_input":"2021-05-23T14:00:37.823548Z","iopub.status.idle":"2021-05-23T14:00:43.215586Z","shell.execute_reply.started":"2021-05-23T14:00:37.823514Z","shell.execute_reply":"2021-05-23T14:00:43.214204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\nimport matplotlib.pyplot as plt\nimport seaborn as sn\nval_classif_preds_01 =[0 if x < 0.5 else 1 for x in val_classif_preds]\nconfusion_mtx = confusion_matrix(val_classif_gts, val_classif_preds_01)\n\nax = plt.axes()\nsn.heatmap(confusion_mtx, annot=True,annot_kws={\"size\": 25}, cmap=\"Reds\", ax = ax)\nax.set_title('Validation Accuracy Inception', size=14)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:00:50.186251Z","iopub.execute_input":"2021-05-23T14:00:50.186654Z","iopub.status.idle":"2021-05-23T14:00:50.414374Z","shell.execute_reply.started":"2021-05-23T14:00:50.186615Z","shell.execute_reply":"2021-05-23T14:00:50.413491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.eval()\nval_vCDRs = []\nval_classif_gts = []\nval_loss = 0.\nval_dsc_od = 0.\nval_dsc_oc = 0.\nval_vCDR_error = 0.\nwith torch.no_grad():\n val_data = iter(val_loader)\n # Loads data\n imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n # Forward pass\n logits_initial = model(imgs)\n val_loss += seg_loss(logits_initial, seg_gts).item() / nb_val_batches\n\n # Compute segmentation metric\n pred_od = refine_seg((logits_initial[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits_initial[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n gt_od = seg_gts[:,0,:,:].type(torch.int8)\n gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n dsc_od = compute_dice_coef(pred_od, gt_od)\n dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n val_dsc_od += dsc_od.item()/nb_val_batches\n val_dsc_oc += dsc_oc.item()/nb_val_batches\n \nfig, axs = plt.subplots(5, 1, figsize=(32, 32))\nlogits_initial = logits_initial.cpu()\nimgs = imgs.cpu()\ngt_od = gt_od.cpu()\ngt_oc = gt_oc.cpu()\n\naxs[0].imshow(imgs[0][0])\naxs[1].imshow(logits_initial[0][0])\naxs[2].imshow(logits_initial[0][1])\naxs[3].imshow(gt_od[0])\naxs[4].imshow(gt_oc[0])\n\n# axs[0].title.set_text(\"Retinal fundus image\")\n# axs[1].title.set_text(\"Its associated mask\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predictions on test set","metadata":{}},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs = []\nwith torch.no_grad():\n test_data = iter(test_loader)\n for k in range(nb_test_batches):\n # Loads data\n imgs = test_data.next()\n imgs = imgs.to(device)\n\n # Forward pass\n logits = model(imgs)\n\n # Std out\n print('Test iter {}/{}'.format(k+1, nb_test_batches) + ' '*50, \n end='\\r')\n \n # Compute segmentation\n pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n \n # Compute and store vCDRs\n pred_vCDR = vertical_cup_to_disc_ratio(pred_od.cpu().numpy(), pred_oc.cpu().numpy())\n test_vCDRs += pred_vCDR.tolist()\n \n\n # Glaucoma predictions from vCDRs\n test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n test_classif_preds = clf.predict_proba(test_vCDRs)[:,1]\n \n# Prepare and save .csv file\ndef create_submission_csv(prediction, submission_filename='/kaggle/working/submission2.csv'):\n \"\"\"Create a sumbission file in the appropriate format for evaluation.\n\n :param\n prediction: list of predictions (ex: [0.12720, 0.89289, ..., 0.29829])\n \"\"\"\n \n with open(submission_filename, mode='w') as csv_file:\n fieldnames = ['Id', 'Predicted']\n writer = csv.DictWriter(csv_file, fieldnames=fieldnames)\n writer.writeheader()\n\n for i, p in enumerate(prediction):\n writer.writerow({'Id': \"T{:04d}\".format(i+1), 'Predicted': '{:f}'.format(p)})\n\ncreate_submission_csv(test_classif_preds)\n\n# The submission.csv file is under /kaggle/working/submission.csv.\n# If you want to submit it, you should download it before closing the current kernel.","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:01:34.372844Z","iopub.execute_input":"2021-05-23T14:01:34.373225Z","iopub.status.idle":"2021-05-23T14:01:39.419527Z","shell.execute_reply.started":"2021-05-23T14:01:34.373186Z","shell.execute_reply":"2021-05-23T14:01:39.417732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MODEL 2 : BASIC CNN","metadata":{}},{"cell_type":"code","source":"import os\nimport json\nimport csv\nimport random\nimport pickle\nimport cv2\nimport numpy as np\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nimport torchvision.transforms as transforms\n\nfrom PIL import Image\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy.ndimage.measurements import label\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import roc_auc_score, roc_curve","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# importing the libraries\nimport pandas as pd\nimport numpy as np\n\n# for reading and displaying images\nfrom skimage.io import imread\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\n# for creating validation set\nfrom sklearn.model_selection import train_test_split\n\n# for evaluating the model\nfrom sklearn.metrics import accuracy_score\nfrom tqdm import tqdm\n\n# PyTorch libraries and modules\nimport torch\nfrom torch.autograd import Variable\nfrom torch.nn import Linear, ReLU, CrossEntropyLoss, Sequential, Conv2d, MaxPool2d, Module, Softmax, BatchNorm2d, Dropout\nfrom torch.optim import Adam, SGD","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Dataset: Loading images**","metadata":{}},{"cell_type":"code","source":"# loading training images\ntrain_img = []\ntrain_y=[]\n\nimport cv2\nimport os\n\ndef load_images_from_folder(folder):\n images = []\n for filename in os.listdir(folder):\n img = cv2.imread(os.path.join(folder,filename))\n # converting the type of pixel to float 32\n img = cv2.resize(img,(256,256))\n img = img.astype('float32')\n # normalizing the pixel values\n img /= 255.0\n \n if(filename[0]=='g'):\n train_y.append(1)\n else:\n train_y.append(0)\n if img is not None:\n images.append(img)\n \n return images\ntrain_img=load_images_from_folder('../input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/images')\n# converting the list to numpy array\ntrain_x = np.array(train_img)\ntrain_y=np.array(train_y)\n\ntrain_x.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def classif_eval(classif_preds, classif_gts):\n '''\n Compute AUC classification score.\n '''\n auc = roc_auc_score(classif_gts, classif_preds)\n return auc\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Train and Validation**","metadata":{}},{"cell_type":"code","source":"train_x, val_x, train_y, val_y = train_test_split(train_x, train_y, test_size = 0.1)\n(train_x.shape, train_y.shape), (val_x.shape, val_y.shape)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# converting training images into torch format\ntrain_x = train_x.reshape(360, 3,256, 256)\ntrain_x = torch.from_numpy(train_x)\n\n# converting the target into torch format\ntrain_y = train_y.astype(int);\ntrain_y = torch.from_numpy(train_y)\n\n# shape of training data\ntrain_x.shape, train_y.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# converting validation images into torch format\nval_x = val_x.reshape(40, 3, 256,256)\nval_x = torch.from_numpy(val_x)\n\n# converting the target into torch format\nval_y = val_y.astype(int);\nval_y = torch.from_numpy(val_y)\n\n# shape of validation data\nval_x.shape, val_y.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Net(Module): \n def __init__(self):\n super(Net, self).__init__()\n\n self.cnn_layers = Sequential(\n # Defining a 2D convolution layer\n Conv2d(3, 4, kernel_size=3, stride=1, padding=1),\n BatchNorm2d(4),\n ReLU(inplace=True),\n MaxPool2d(kernel_size=2, stride=2),\n # Defining another 2D convolution layer\n Conv2d(4, 4, kernel_size=3, stride=1, padding=1),\n BatchNorm2d(4),\n ReLU(inplace=True),\n MaxPool2d(kernel_size=2, stride=2),\n )\n\n self.linear_layers = Sequential(\n Linear(16384, 10)\n )\n\n # Defining the forward pass \n def forward(self, x):\n x = self.cnn_layers(x)\n x = x.view(x.size(0), -1)\n x = self.linear_layers(x)\n return x","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# defining the model\nmodel = Net()\n# defining the optimizer\noptimizer = Adam(model.parameters(), lr=0.07)\n# defining the loss function\ncriterion = CrossEntropyLoss()\n# checking if GPU is available\nif torch.cuda.is_available():\n model = model.cuda()\n criterion = criterion.cuda()\n \nprint(model)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train(epoch):\n model.train()\n tr_loss = 0\n # getting the training set\n x_train, y_train = Variable(train_x), Variable(train_y)\n # getting the validation set\n x_val, y_val = Variable(val_x), Variable(val_y)\n # converting the data into GPU format\n if torch.cuda.is_available():\n x_train = x_train.cuda()\n y_train = y_train.cuda()\n x_val = x_val.cuda()\n y_val = y_val.cuda()\n\n # clearing the Gradients of the model parameters\n optimizer.zero_grad()\n \n # prediction for training and validation set\n output_train = model(x_train)\n output_val = model(x_val)\n\n # computing the training and validation loss\n loss_train = criterion(output_train, y_train)\n loss_val = criterion(output_val, y_val)\n aux_train=classif_eval(output_train.cpu().detach().numpy()[:,1].tolist(), y_train.cpu().detach().numpy().tolist())\n auc_val=classif_eval(output_val.cpu().detach().numpy()[:,1].tolist(), y_val.cpu().detach().numpy().tolist())\n train_losses.append(loss_train)\n val_losses.append(loss_val)\n\n # computing the updated weights of all the model parameters\n loss_train.backward()\n optimizer.step()\n tr_loss = loss_train.item()\n if epoch%2 == 0:\n # printing the validation loss\n print('Epoch : ',epoch+1, '\\t', 'loss :', loss_val,'AUC :',auc_val)\n ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# defining the number of epochs\nn_epochs = 25\n# empty list to store training losses\ntrain_losses = []\n# empty list to store validation losses\nval_losses = []\n# training the model\nfor epoch in range(n_epochs):\n train(epoch)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(train_losses, label='Training loss')\nplt.plot(val_losses, label='Validation loss')\nplt.legend()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prediction for training set\nwith torch.no_grad():\n output = model(train_x.cuda())\n \nsoftmax = torch.exp(output).cpu()\nprob = list(softmax.numpy())\npredictions = np.argmax(prob, axis=1)\n\n# accuracy on training set\naccuracy_score(train_y, predictions)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**TEST**","metadata":{}},{"cell_type":"code","source":"\n# loading test images\ntest_img = []\ndef test_load_images_from_folder(folder):\n images = []\n for filename in os.listdir(folder):\n img = cv2.imread(os.path.join(folder,filename))\n # converting the type of pixel to float 32\n img = cv2.resize(img,(256,256))\n img = img.astype('float32')\n # normalizing the pixel values\n img /= 255.0\n images.append(img)\n \n \n return images\ntest_img=test_load_images_from_folder('../input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/test/images')\n# converting the list to numpy array\n\ntest_x = np.array(test_img)\ntest_x.shape\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# converting training images into torch format\ntest_x = test_x.reshape(400, 3, 256, 256)\ntest_x = torch.from_numpy(test_x)\ntest_x.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# generating predictions for test set\nwith torch.no_grad():\n output = model(test_x.cuda())\n\nsoftmax = torch.exp(output).cpu()\nprob = list(softmax.numpy())\npredictions = np.argmax(prob, axis=1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MODEL 3 :VGG-16","metadata":{}},{"cell_type":"code","source":"from IPython.display import clear_output\n!pip install imutils","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set everything up\nimport os\nimport cv2\nimport imutils as imutils\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport tensorflow as tf # machine learning\nfrom tqdm import tqdm # make your loops show a smart progress meter \nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nfrom sklearn.metrics import accuracy_score, confusion_matrix\nimport seaborn as sn\n\nRANDOM_SEED = 1\nIMG_SIZE = (256, 256) # size of vgg16 input","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Creation of the Dataset

","metadata":{}},{"cell_type":"markdown","source":"

This code is not runnable and was executed on a local notebook. The Dataset was then put on the kaggle plateform. ","metadata":{}},{"cell_type":"markdown","source":"

validation preparation

","metadata":{}},{"cell_type":"code","source":"#To prepare our data for their utilisation in the classification,we had to rearrange the directories.\n#To do so, we used another notebook and used this code : \n#validation_root = \"/home/iversenc/Documents/Telecom/Cours/eurecom-aml-2021-challenge-2/refuge_data/val/\"
\n#json = pandas.read_json(validation_root+\"index.json\")
\n#os.makedirs(validation_root + \"img_glaucomia\", exist_ok=True)
\n#os.makedirs(validation_root + \"img_normal\", exist_ok=True)
\n#for img in os.listdir(validation_root):
\n# print(img)
\n#print(json[i])
\n#for i in range(400):
\n# name = json[i].ImgName
\n# label = json[i].Label
\n# if(label==0):
\n# shutil.move(validation_root+\"images/\"+name, validation_root+\"img_normal/\"+name)
\n# else:
\n# shutil.move(validation_root+\"images/\"+name, validation_root+\"img_glaucomia/\"+name)
\n# \n#

This allowed us to put 2 directories for the sick and not sick images. This was done thank's to the index.

","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

training preparation

","metadata":{}},{"cell_type":"code","source":"#pip install tensorflowpip install tensorflow\n#!pip install keras\n#from keras.preprocessing.image import ImageDataGenerator, load_img, img_to_array\n#import os\n#import shutil","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train_root = \"/home/iversenc/Documents/Telecom/Cours/AML/eurecom/refuge_data/train2/\" \n#os.makedirs(train_root + \"img_glaucomia\", exist_ok=True)\n#os.makedirs(train_root + \"img_normal\", exist_ok=True)\n#for img in os.listdir(train_root):\n# print(img)\n#for img in os.listdir(train_root+\"images/\"):\n# name = list(img)\n# if(name[0]=='g'):\n# shutil.move(train_root+\"images/\"+img, train_root+\"img_glaucomia/\"+img)\n# else:\n# shutil.move(train_root+\"images/\"+img, train_root+\"img_normal/\"+img)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#for file in os.listdir(train_root+\"img_glaucomia\"):\n# name=file.split('.')\n# print(name)\n# if(len(name[0])==5):\n# for i in range(8):\n# shutil.copyfile(train_root+\"img_glaucomia/\"+file, train_root+\"img_glaucomia/\"+name[0]+str(i)+'6'+\".jpg\")\n#print(len(os.listdir(train_root+\"img_normal\")))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#datagen = ImageDataGenerator(\n # rotation_range=40, horizontal_flip=True, vertical_flip=True,\n # fill_mode='nearest')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Actual data_augmentation\n#image_mul = 6 # How many images we generate from 1 image\n#for img in os.listdir(train_root+\"img_glaucomia\"):\n # file_save_prefix = img.split(\".\")[0]\n # pic = load_img(train_root + 'img_glaucomia/' + img)\n # pic_array = img_to_array(pic)\n #pic_array = pic_array.reshape((1,) + pic_array.shape)\n #count = 0\n #for batch in datagen.flow(pic_array, batch_size=1, save_to_dir=train_root+\"img_glaucomia/\", save_prefix=file_save_prefix+str(count), save_format='jpg'):\n # count += 1\n # if count == image_mul:\n # break","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Definition of the paths

\n","metadata":{}},{"cell_type":"markdown","source":"

We saw that our data was quite limited. Thus we decided to create more images based on the ones from the dataset.
\n This also allowed us to have more data samples where the glaucomia is present. (more than 40)","metadata":{}},{"cell_type":"code","source":"train_root = \"../input/training-data-augmented/augmented_images_v3/augmented_images_v3/train_imgs/\"\ntest_root = \"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/test/images/\"\nvalidation_root = \"/kaggle/input/val-classified/val_classified/\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def repartition(path):\n print(path)\n for value in os.listdir(path):\n print(\"There is \"+ str(len(os.listdir(path + value))) + \" images for \" + value)\n print('\\n')\n \nprint(\"let's look at the repartition of the images in the training data\")\nrepartition(train_root)\nprint(\"let's look at the repartition of the images in the validation data\")\nrepartition(validation_root)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Let's visualize some of our images.
","metadata":{}},{"cell_type":"markdown","source":"

We already augmented the data so that we did not have to do it everytime we run our notebook.

","metadata":{}},{"cell_type":"code","source":"train_gen = tf.keras.preprocessing.image.ImageDataGenerator(\n preprocessing_function=tf.keras.applications.vgg16.preprocess_input\n )\n\ntest_gen = tf.keras.preprocessing.image.ImageDataGenerator(\n preprocessing_function=tf.keras.applications.vgg16.preprocess_input\n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_generator = train_gen.flow_from_directory(\n train_root,\n color_mode='rgb',\n target_size=IMG_SIZE,\n batch_size=10,\n shuffle = True,\n class_mode='binary',\n seed=RANDOM_SEED\n)\n\nvalidation_generator = test_gen.flow_from_directory(\n validation_root,\n color_mode='rgb',\n shuffle = True,\n target_size=IMG_SIZE,\n batch_size=20,\n class_mode='binary'\n )","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vgg16_weight_path = '../input/vgg16/vgg16_weights_tf_dim_ordering_tf_kernels_notop.h5'\nbase_model = tf.keras.applications.VGG16(\n weights=\"imagenet\",\n include_top=False,\n input_shape=(256,256,3)\n)\nmodel = tf.keras.models.Sequential()\nmodel.add(base_model)\nmodel.add(tf.keras.layers.Flatten())\nmodel.add(tf.keras.layers.Dropout(0.5))\nmodel.add(tf.keras.layers.Dense(1, activation='sigmoid'))\n\nmodel.layers[0].trainable = False\n\nmodel.compile(\n loss='binary_crossentropy',\n optimizer=tf.keras.optimizers.Adam(),\n metrics=[tf.keras.metrics.AUC(),'accuracy']\n)\n\nmodel.summary()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EPOCHS = 20 \nearly_stopping = tf.keras.callbacks.EarlyStopping(\n monitor='val_auc',\n mode='max',\n patience=6\n)\ncallbacks2 = keras.callbacks.ModelCheckpoint(\n # Path where to save the model\n # The two parameters below mean that we will overwrite\n # the current checkpoint if and only if\n # the `val_loss` score has improved.\n # The saved model name will include the current epoch.\n filepath=\"mymodelvgg_{epoch}\",\n save_best_only=True, # Only save a model if `val_loss` has improved.\n monitor=\"val_auc\",\n verbose=1,\n mode='max'\n )\n\n\nhistory = model.fit_generator(\n train_generator,\n steps_per_epoch=20,\n shuffle=True,\n epochs=EPOCHS,\n validation_data=validation_generator,\n validation_steps=15,\n callbacks=[early_stopping]\n)\n\nprint(\"Training Done\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.evaluate(validation_generator)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Model efficiency

<\\center>","metadata":{}},{"cell_type":"code","source":"def preprocess_imgs(path, img_size):\n set_new = []\n reality_array=[]\n for value in os.listdir(path):\n for img in os.listdir(path + value):\n img = cv2.imread(path + value + \"/\" + img)\n img = cv2.resize(\n img,\n dsize=img_size,\n interpolation=cv2.INTER_CUBIC\n )\n if(value==\"img_glaucomia\"):\n reality_array.append(1)\n else:\n reality_array.append(0)\n set_new.append(tf.keras.applications.vgg16.preprocess_input(img))\n \n return np.array(set_new),reality_array\n\nvalidation_data ,reality= preprocess_imgs(validation_root, img_size=IMG_SIZE)\npredictions = model.predict(validation_data)\npredictions = [0 if x > 0.5 else 1 for x in predictions]\naccuracy = accuracy_score(reality, predictions)\nprint(\"Validation Accuracy for VGG16:\", accuracy)\n\n\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc = history.history['auc']\nloss = history.history['loss']\nval_auc = history.history['val_auc']\nval_loss = history.history['val_loss']\nepochs_range = range(1, len(history.epoch) + 1)\n\n\nplt.figure(figsize=(10,5))\n\nplt.plot(epochs_range, acc, label='Training')\nplt.plot(epochs_range, val_auc, label='Validation')\nplt.legend(loc=\"best\")\nplt.xlabel('Epochs')\nplt.ylabel('Accuracy')\nplt.title('Model Accuracy')\nplt.grid(b=True, which='major', color='#666666', linestyle='-')\nplt.tight_layout()\nplt.show()\n\nplt.figure(figsize=(10,5))\n\nplt.plot(epochs_range, loss, label='Training')\nplt.plot(epochs_range, val_loss, label='Validation')\nplt.legend(loc=\"best\")\nplt.xlabel('Epochs')\nplt.ylabel('Loss')\nplt.title('Model Loss')\nplt.grid(b=True, which='major', color='#666666', linestyle='-')\nplt.tight_layout()\nplt.show()\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"confusion_mtx = confusion_matrix(reality, predictions)\n\nax = plt.axes()\nsn.heatmap(confusion_mtx, annot=True,annot_kws={\"size\": 25}, cmap=\"Blues\", ax = ax)\nax.set_title('Validation Accuracy VGG16', size=14)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model 4 : Inception v3","metadata":{}},{"cell_type":"markdown","source":"

Let's try a new model of CNN classifier to see if there is an evolution

\n

Indeed, google implemented a pretrained CNN called Inception

","metadata":{}},{"cell_type":"code","source":"train_generator_inception = train_gen.flow_from_directory(\n train_root,\n color_mode='rgb',\n target_size=IMG_SIZE,\n batch_size=10,\n shuffle = True,\n class_mode='binary',\n seed=RANDOM_SEED\n)\n\nvalidation_generator_inception = test_gen.flow_from_directory(\n validation_root,\n color_mode='rgb',\n shuffle = True,\n target_size=IMG_SIZE,\n batch_size=32,\n class_mode='binary'\n )","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras import layers\nfrom tensorflow.keras.applications.inception_v3 import InceptionV3\nbase_model_inception = InceptionV3(input_shape = (256, 256, 3), include_top = False, weights = 'imagenet')\n\nfrom tensorflow.keras.optimizers import RMSprop\nimport keras as keras\nx = layers.Flatten()(base_model_inception.output)\nx = layers.Dense(1024, activation='relu')(x)\nx = layers.Dropout(0.2)(x)\nx = layers.Dense(1, activation='sigmoid')(x)\n\nmodel_inception = tf.keras.models.Model(base_model_inception.input, x)\n\nmodel_inception.compile(optimizer = RMSprop(lr=0.0002), loss = 'binary_crossentropy', metrics = [tf.keras.metrics.AUC(),'accuracy'])\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"early_stopping_inception = tf.keras.callbacks.EarlyStopping(\n monitor='val_auc_2',\n mode='max',\n patience=6\n)\ncallbacks2 = keras.callbacks.ModelCheckpoint(\n # Path where to save the model\n # The two parameters below mean that we will overwrite\n # the current checkpoint if and only if\n # the `val_loss` score has improved.\n # The saved model name will include the current epoch.\n filepath=\"mymodel3_{epoch}\",\n save_best_only=True, # Only save a model if `val_loss` has improved.\n monitor=\"val_auc_2\",\n verbose=1,\n mode='max'\n )\n\ninc_history = model_inception.fit_generator(\n train_generator_inception, \n validation_data = validation_generator_inception, \n shuffle = True,\n steps_per_epoch = 100, \n epochs = 20,\n callbacks=[early_stopping_inception,callbacks2])\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import keras\nbest_model = keras.models.load_model(\"../input/model-087/model_inception1.h5\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# at this point, the top layers are well trained and we can start fine-tuning\n# convolutional layers from inception V3. We will freeze the bottom N layers\n# and train the remaining top layers.\n\n# let's visualize layer names and layer indices to see how many layers\n# we should freeze:\nfor i, layer in enumerate(best_model.layers):\n print(i, layer.name)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we chose to train the top 2 inception blocks, i.e. we will freeze\n# the first 249 layers and unfreeze the rest:\n# in other examples found it was 172 insted 249. \n# I took 249 according to https://keras.io/applications/#inceptionv3\nfor layer in best_model.layers[:249]:\n layer.trainable = False\nfor layer in best_model.layers[249:]:\n layer.trainable = True","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we need to recompile the model for these modifications to take effect\n# we use SGD with a low learning rate\nfrom keras.optimizers import SGD\n\nbest_model.compile(optimizer=SGD(lr=0.0001, momentum=0.9), loss='binary_crossentropy', metrics=[tf.keras.metrics.AUC(),'accuracy'])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inc_history = best_model.fit_generator(\n train_generator_inception, \n validation_data = validation_generator_inception, \n shuffle = True,\n steps_per_epoch = 50, \n epochs = 20)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Model efficiency

<\\center>","metadata":{}},{"cell_type":"markdown","source":"

Let's predict the validation data and evaluate the model on data on which he did not learn.

","metadata":{}},{"cell_type":"code","source":"validation_data ,reality= preprocess_imgs(validation_root, img_size=IMG_SIZE)\n\npredictions_inception = model_inception.predict(validation_data)\npredictions_inception = [0 if x > 0.5 else 1 for x in predictions_inception]\naccuracy_inception = accuracy_score(reality,predictions_inception)\n\nbest_predict = best_model.predict(validation_data)\nbest_predict = [0 if x > 0.5 else 1 for x in best_predict]\naccuracy_best_predict = accuracy_score(reality,best_predict)\n\nprint(\"Validation Accuracy for inception:\", accuracy_inception)\nprint(\"Validation Accuracy for inception best model:\", accuracy_best_predict)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_model.evaluate(validation_generator_inception)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Nevertheless, we figured that since our data was not equilibrated, thus the accuracy was not to take into consideration.

\n

A confusion matrix will be more appropriate to evaluate our model.

","metadata":{}},{"cell_type":"code","source":"auc = inc_history.history['auc_1']\nloss = inc_history.history['loss']\nval_loss = inc_history.history['val_loss']\nval_auc = inc_history.history['val_auc_1']\n\nepochs_range = range(1, len(inc_history.epoch) + 1)\n\n\nplt.figure(figsize=(10,5))\n\nplt.plot(epochs_range, auc, label='Training')\nplt.plot(epochs_range, val_auc, label='Validation')\nplt.legend(loc=\"best\")\nplt.xlabel('Epochs')\nplt.ylabel('Accuracy')\nplt.title('Model Accuracy')\nplt.grid(b=True, which='major', color='#666666', linestyle='-')\nplt.tight_layout()\nplt.show()\n\nplt.figure(figsize=(10,5))\n\nplt.plot(epochs_range, loss, label='Training')\nplt.plot(epochs_range, val_loss, label='Validation')\nplt.legend(loc=\"best\")\nplt.xlabel('Epochs')\nplt.ylabel('Loss')\nplt.title('Model Loss')\nplt.grid(b=True, which='major', color='#666666', linestyle='-')\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"confusion_mtx = confusion_matrix(reality, predictions_inception)\n\nax = plt.axes()\nsn.heatmap(confusion_mtx, annot=True,annot_kws={\"size\": 25}, cmap=\"Reds\", ax = ax)\nax.set_title('Validation Accuracy Inception', size=14)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"confusion_mtx = confusion_matrix(reality,best_predict)\n\nax = plt.axes()\nsn.heatmap(confusion_mtx, annot=True,annot_kws={\"size\": 25}, cmap=\"Reds\", ax = ax)\nax.set_title('Validation Accuracy Inception', size=14)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's save our model","metadata":{}},{"cell_type":"code","source":"tf.keras.models.save_model(\n model_inception, \"./model_inception2.h5\", overwrite=False, save_format=None,\n signatures=None, options=None, save_traces=True\n)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Submit the test data

","metadata":{}},{"cell_type":"code","source":"import cv2\ndef preprocess_imgs_test(path, img_size):\n set_new = []\n names=[]\n for img in os.listdir(path):\n \n names.append(img.split('.')[0])\n img = cv2.imread(path +img)\n img = cv2.resize(img,\n img_size,\n interpolation = cv2.INTER_CUBIC)\n set_new.append(img)\n return np.array(set_new), names\n\ntest_images,names = preprocess_imgs_test(test_root,IMG_SIZE)\npredictions = best_model.predict(test_images)\n#predictionss = model_inception.predict(test_images)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport csv\n# Prepare and save .csv file\ndef create_submission_csv(predictions, names, submission_filename='/kaggle/working/submission5.csv'):\n with open(submission_filename, mode='w') as csv_file:\n fieldnames = ['Id', 'Predicted']\n writer = csv.DictWriter(csv_file, fieldnames=fieldnames)\n writer.writeheader()\n for i in range(len(predictions)):\n \n writer.writerow({'Id':names[i], 'Predicted': 1-float(predictions[i])})\ncreate_submission_csv(predictions,names)","metadata":{},"execution_count":null,"outputs":[]}]}