{"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":"markdown","source":"# Introduction\nThis notebook provides solution to [Histopathologic Cancer Detection](https://www.kaggle.com/c/histopathologic-cancer-detection/overview) challenge on Kaggle. This is a perfect Computer Vision problem where we are tasked with the detection of cancer by identifying metastatic tissue in histopathologic scans of lymph nodes using Deep Learning.\n\n![Header Image](https://storage.googleapis.com/kaggle-competitions/kaggle/11848/logos/header.png?t=2018-11-15-01-52-19)\n\n\n### 1. Understanding the Problem:\nOur goal is to create an algorithm to identify metastatic cancer in small image patches taken from larger digital pathology scans. \n\nObviously I don't know biology to understand this problem right away, here is what I found online about histopathology.\n\n> Histopathology is the study of the signs of the disease using the microscopic examination of a biopsy or surgical specimen that is processed and fixed onto glass slides. To visualize different components of the tissue under a microscope, the sections are dyed with one or more stains.\n\n### Motivation:\nLymph nodes are small glands that filter the fluid in the lymphatic system and they are the first place a breast cancer is likely to spread. Histological assessment of lymph node metastases is part of determining the stage of breast cancer in TNM classification which is a globally recognized standard for classifying the extent of spread of cancer. \n> **The diagnostic procedure for pathologists is tedious and time-consuming as a large area of tissue has to be examined and small metastases can be easily missed.** \n\nThat makes using Machine Learning a great choice both in terms of accuracy and ease of usability. It could bring a great change altogether. \n\n### 2. Understanding the Data:\n\n**The train data we have here contains 220,025 images and the test set contains 57,468 images.** \n\nIt is important to take into account that this data is only a subset of the original [PCam dataset](https://github.com/basveeling/pcam) which in the end is derived from the [Camelyon16 Challenge dataset](https://camelyon16.grand-challenge.org/Data/), which contains 400 H&E stained whole slide images of sentinel lymph node sections that were acquired and digitized at 2 different centers using a 40x objective. The PCam's dataset including this one uses 10x undersampling to increase the field of view, which gives the resultant pixel resolution of 2.43 microns.\n\nHere's what Kaggle says,\n\n> The original PCam dataset contains duplicate images due to its probabilistic sampling, however, the version presented on Kaggle does not contain duplicates. We have otherwise maintained the same data and splits as the PCam benchmark.\n\nOur training data has a class distribution of 60:40 negative and positive samples which is not bad.\n\nI also found that these data were obtained as a result of routine clinical practices and similar to how a trained pathologist would examine similar images for identifying metastases. However, some relevant information about the surroundings might be left out with these small-sized image samples (I guess).\n\n### 3. Understanding the Images\n > You are predicting the labels for the images in the test folder. A positive label indicates that the center 32x32px region of a patch contains at least one pixel of tumor tissue. Tumor tissue in the outer region of the patch does not influence the label. This outer region is provided to enable fully-convolutional models that do not use zero-padding, to ensure consistent behavior when applied to a whole-slide image.\n \nThis from the competition's description means that the centers of the images are the ones that really matter.\n\nAs you might already know, **this is a binary classification problem**.\n\n### 4. Understanding the Evaluation Metric\nThe evaluation metric is the **Area Under ROC Curve** which is also called **AU-ROC/AOC Curve**. It is one of the most important evaluation metrics for checking any classification model’s performance.\n\nAUC - ROC curve is a performance measurement for classification problem at various thresholds settings. ROC is a probability curve and AUC represents degree or measure of separability. **It tells how much model is capable of distinguishing between classes.** Higher the AUC, better the model is at predicting 0s as 0s and 1s as 1s. By analogy, higher the AUC (close to 1), better the model is at distinguishing between patients with disease and no disease. The curve is plotted with True Positive Rates Vs the False Positive Rates along the x and y axes respectively.\n\n\nROC                        |  AUC \n:-------------------------:|:-------------------------:\n![ROC Curve](http://gim.unmc.edu/dxtests/roccomp.jpg)  |   ![AUC Curve](https://i.ibb.co/mBKh6ZB/roc.pnghttps://i.ibb.co/mBKh6ZB/roc.png)\n","metadata":{}},{"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)\n\nimport os\nbase_dir = '../input/'\nprint(os.listdir(base_dir))\n\n# Matplotlib for visualization\nimport matplotlib.pyplot as plt\nplt.style.use(\"ggplot\")\n\n# OpenCV Image Library\nimport cv2\n\n# Import PyTorch\nimport torchvision.transforms as transforms\nfrom torch.utils.data.sampler import SubsetRandomSampler\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import TensorDataset, DataLoader, Dataset\nimport torchvision\nimport torch.optim as optim\n\n# Import useful sklearn functions\nimport sklearn\nfrom sklearn.metrics import roc_auc_score, accuracy_score\nfrom PIL import Image","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading Data and EDA\nHaving a look at the data, just like any other image classification problem we have a csv file with image ids and labels. The directories train, test contain the actual images.","metadata":{}},{"cell_type":"code","source":"full_train_df = pd.read_csv(\"../input/train_labels.csv\")\nfull_train_df.head()","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Train Size: {}\".format(len(os.listdir('../input/train/'))))\nprint(\"Test Size: {}\".format(len(os.listdir('../input/test/'))))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels_count = full_train_df.label.value_counts()\n\n%matplotlib inline\nplt.pie(labels_count, labels=['No Cancer', 'Cancer'], startangle=180, \n        autopct='%1.1f', colors=['#00ff99','#FF96A7'], shadow=True)\nplt.figure(figsize=(16,16))\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing Images\nClassifying metastases is probably not an easy task for a trained pathologist and extremely difficult for an untrained eye when we take a look at the image.","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(30, 6))\n# display 20 images\ntrain_imgs = os.listdir(base_dir+\"train\")\nfor idx, img in enumerate(np.random.choice(train_imgs, 20)):\n    ax = fig.add_subplot(2, 20//2, idx+1, xticks=[], yticks=[])\n    im = Image.open(base_dir+\"train/\" + img)\n    plt.imshow(im)\n    lab = full_train_df.loc[full_train_df['id'] == img.split('.')[0], 'label'].values[0]\n    ax.set_title('Label: %s'%lab)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sampling\nSince the train dataset contains 220.025 images we can sample out a shuffled part of that, in this case 160000 samples and train on them to make predictions later. ","metadata":{}},{"cell_type":"code","source":"# Number of samples in each class\nSAMPLE_SIZE = 80000\n\n# Data paths\ntrain_path = '../input/train/'\ntest_path = '../input/test/'\n\n# Use 80000 positive and negative examples\ndf_negatives = full_train_df[full_train_df['label'] == 0].sample(SAMPLE_SIZE, random_state=42)\ndf_positives = full_train_df[full_train_df['label'] == 1].sample(SAMPLE_SIZE, random_state=42)\n\n# Concatenate the two dfs and shuffle them up\ntrain_df = sklearn.utils.shuffle(pd.concat([df_positives, df_negatives], axis=0).reset_index(drop=True))\n\ntrain_df.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Pre-processing for our PyTorch\nFirst we turn our data into PyTorch dataset then the data is sampled into train and validation sets. Data Augmentations are added for train data to improve performance.","metadata":{}},{"cell_type":"code","source":"# Our own custom class for datasets\nclass CreateDataset(Dataset):\n    def __init__(self, df_data, data_dir = './', transform=None):\n        super().__init__()\n        self.df = df_data.values\n        self.data_dir = data_dir\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.df)\n    \n    def __getitem__(self, index):\n        img_name,label = self.df[index]\n        img_path = os.path.join(self.data_dir, img_name+'.tif')\n        image = cv2.imread(img_path)\n        if self.transform is not None:\n            image = self.transform(image)\n        return image, label","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"transforms_train = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.RandomHorizontalFlip(p=0.4),\n    transforms.RandomVerticalFlip(p=0.4),\n    transforms.RandomRotation(20),\n    transforms.ToTensor(),\n    # We the get the following mean and std for the channels of all the images\n    #transforms.Normalize((0.70244707, 0.54624322, 0.69645334), (0.23889325, 0.28209431, 0.21625058))\n    transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))\n])\n\ntrain_data = CreateDataset(df_data=train_df, data_dir=train_path, transform=transforms_train)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set Batch Size\nbatch_size = 128\n\n# Percentage of training set to use as validation\nvalid_size = 0.1\n\n# obtain training indices that will be used for validation\nnum_train = len(train_data)\nindices = list(range(num_train))\n# np.random.shuffle(indices)\nsplit = int(np.floor(valid_size * num_train))\ntrain_idx, valid_idx = indices[split:], indices[:split]\n\n# Create Samplers\ntrain_sampler = SubsetRandomSampler(train_idx)\nvalid_sampler = SubsetRandomSampler(valid_idx)\n\n# prepare data loaders (combine dataset and sampler)\ntrain_loader = DataLoader(train_data, batch_size=batch_size, sampler=train_sampler)\nvalid_loader = DataLoader(train_data, batch_size=batch_size, sampler=valid_sampler)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"transforms_test = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.ToTensor(),\n    #transforms.Normalize((0.70244707, 0.54624322, 0.69645334), (0.23889325, 0.28209431, 0.21625058))\n    transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))\n])\n\n# creating test data\nsample_sub = pd.read_csv(\"../input/sample_submission.csv\")\ntest_data = CreateDataset(df_data=sample_sub, data_dir=test_path, transform=transforms_test)\n\n# prepare the test loader\ntest_loader = DataLoader(test_data, batch_size=batch_size, shuffle=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Defining Model Architecture\nI'm using a Deep Convolutional Neural Network for this task building which is fairly straight-forward in PyTorch if you understand how it works. This is one of many architectures I tried that gave better results.","metadata":{}},{"cell_type":"code","source":"class CNN(nn.Module):\n    def __init__(self):\n        super(CNN,self).__init__()\n        # Convolutional and Pooling Layers\n        self.conv1=nn.Sequential(\n                nn.Conv2d(in_channels=3,out_channels=32,kernel_size=3,stride=1,padding=0),\n                nn.BatchNorm2d(32),\n                nn.ReLU(inplace=True),\n                nn.MaxPool2d(2,2))\n        self.conv2=nn.Sequential(\n                nn.Conv2d(in_channels=32,out_channels=64,kernel_size=2,stride=1,padding=1),\n                nn.BatchNorm2d(64),\n                nn.ReLU(inplace=True),\n                nn.MaxPool2d(2,2))\n        self.conv3=nn.Sequential(\n                nn.Conv2d(in_channels=64,out_channels=128,kernel_size=3,stride=1,padding=1),\n                nn.BatchNorm2d(128),\n                nn.ReLU(inplace=True),\n                nn.MaxPool2d(2,2))\n        self.conv4=nn.Sequential(\n                nn.Conv2d(in_channels=128,out_channels=256,kernel_size=3,stride=1,padding=1),\n                nn.BatchNorm2d(256),\n                nn.ReLU(inplace=True),\n                nn.MaxPool2d(2,2))\n        self.conv5=nn.Sequential(\n                nn.Conv2d(in_channels=256, out_channels=512, kernel_size=3, stride=1, padding=1),\n                nn.BatchNorm2d(512),\n                nn.ReLU(inplace=True),\n                nn.MaxPool2d(2,2))\n        \n        self.dropout2d = nn.Dropout2d()\n        \n        \n        self.fc=nn.Sequential(\n                nn.Linear(512*3*3,1024),\n                nn.ReLU(inplace=True),\n                nn.Dropout(0.4),\n                nn.Linear(1024,512),\n                nn.Dropout(0.4),\n                nn.Linear(512, 1),\n                nn.Sigmoid())\n        \n    def forward(self,x):\n        \"\"\"Method for Forward Prop\"\"\"\n        x=self.conv1(x)\n        x=self.conv2(x)\n        x=self.conv3(x)\n        x=self.conv4(x)\n        x=self.conv5(x)\n        #print(x.shape) <-- Life saving debugging step :D\n        x=x.view(x.shape[0],-1)\n        x=self.fc(x)\n        return x","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check if CUDA is available\ntrain_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 ...')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training and Validation","metadata":{}},{"cell_type":"code","source":"# create a complete CNN\nmodel = CNN()\nprint(model)\n\n# Move model to GPU if available\nif train_on_gpu: model.cuda()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Trainable Parameters\npytorch_total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)\nprint(\"Number of trainable parameters: \\n{}\".format(pytorch_total_params))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# specify loss function (categorical cross-entropy loss)\ncriterion = nn.BCELoss()\n\n# specify optimizer\noptimizer = optim.Adam(model.parameters(), lr=0.00015)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# number of epochs to train the model\nn_epochs = 20\n\nvalid_loss_min = np.Inf\n\n# keeping track of losses as it happen\ntrain_losses = []\nvalid_losses = []\nval_auc = []\ntest_accuracies = []\nvalid_accuracies = []\nauc_epoch = []\n\nfor epoch in range(1, n_epochs+1):\n\n    # keep track of training and validation loss\n    train_loss = 0.0\n    valid_loss = 0.0\n    \n    ###################\n    # train the model #\n    ###################\n    model.train()\n    for data, target in train_loader:\n        # move tensors to GPU if CUDA is available\n        if train_on_gpu:\n            data, target = data.cuda(), target.cuda().float()\n        target = target.view(-1, 1)\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        # calculate the batch loss\n        loss = criterion(output, target)\n        # backward pass: compute gradient of the loss with respect to model parameters\n        loss.backward()\n        # perform a single optimization step (parameter update)\n        optimizer.step()\n        # Update Train loss and accuracies\n        train_loss += loss.item()*data.size(0)\n        \n    ######################    \n    # validate the model #\n    ######################\n    model.eval()\n    for data, target in valid_loader:\n        # move tensors to GPU if CUDA is available\n        if train_on_gpu:\n            data, target = data.cuda(), target.cuda().float()\n        # forward pass: compute predicted outputs by passing inputs to the model\n        target = target.view(-1, 1)\n        output = model(data)\n        # calculate the batch loss\n        loss = criterion(output, target)\n        # update average validation loss \n        valid_loss += loss.item()*data.size(0)\n        #output = output.topk()\n        y_actual = target.data.cpu().numpy()\n        y_pred = output[:,-1].detach().cpu().numpy()\n        val_auc.append(roc_auc_score(y_actual, y_pred))        \n    \n    # calculate average losses\n    train_loss = train_loss/len(train_loader.sampler)\n    valid_loss = valid_loss/len(valid_loader.sampler)\n    valid_auc = np.mean(val_auc)\n    auc_epoch.append(np.mean(val_auc))\n    train_losses.append(train_loss)\n    valid_losses.append(valid_loss)\n        \n    # print training/validation statistics \n    print('Epoch: {} | Training Loss: {:.6f} | Validation Loss: {:.6f} | Validation AUC: {:.4f}'.format(\n        epoch, train_loss, valid_loss, valid_auc))\n    \n    ##################\n    # Early Stopping #\n    ##################\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        torch.save(model.state_dict(), 'best_model.pt')\n        valid_loss_min = valid_loss","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\n%config InlineBackend.figure_format = 'retina'\n\nplt.plot(train_losses, label='Training loss')\nplt.plot(valid_losses, label='Validation loss')\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Loss\")\nplt.legend(frameon=False)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\n%config InlineBackend.figure_format = 'retina'\n\nplt.plot(auc_epoch, label='Validation AUC/Epochs')\nplt.legend(\"\")\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Area Under the Curve\")\nplt.legend(frameon=False)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Best parameters learned from training into our model to make predictions later\nmodel.load_state_dict(torch.load('best_model.pt'))","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predictions on Test set","metadata":{}},{"cell_type":"code","source":"# Turn off gradients\nmodel.eval()\n\npreds = []\nfor batch_i, (data, target) in enumerate(test_loader):\n    data, target = data.cuda(), target.cuda()\n    output = model(data)\n\n    pr = output.detach().cpu().numpy()\n    for i in pr:\n        preds.append(i)\n\n# Create Submission file        \nsample_sub['label'] = preds","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(len(sample_sub)):\n    sample_sub.label[i] = np.float(sample_sub.label[i]) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_sub.to_csv('submission.csv', index=False)\nsample_sub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing Preditions:","metadata":{}},{"cell_type":"code","source":"def imshow(img):\n    '''Helper function to un-normalize and display an image'''\n    # unnormalize\n    img = img / 2 + 0.5\n    # convert from Tensor image and display\n    plt.imshow(np.transpose(img, (1, 2, 0)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# obtain one batch of training images\ndataiter = iter(test_loader)\nimages, labels = dataiter.next()\nimages = images.numpy() # convert images to numpy for display\n\n# plot the images in the batch, along with the corresponding labels\nfig = plt.figure(figsize=(25, 4))\n# display 20 images\nfor idx in np.arange(20):\n    ax = fig.add_subplot(2, 20/2, idx+1, xticks=[], yticks=[])\n    imshow(images[idx])\n    prob = \"Cancer\" if(sample_sub.label[idx] >= 0.5) else \"Normal\" \n    ax.set_title('{}'.format(prob))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"How cool is that? Now this model can be used to predict Cancer, maybe even in real-world, the AUC score I was able to achieve with this model on test set is ~0.95 which shows the model is doing way better than just guessing, it might be very much reliable if a few tweaks are to be made to take it even closer to 1.   ","metadata":{}},{"cell_type":"markdown","source":"### Authored By,\n[Abhinand](http://kaggle.com/abhinand05)","metadata":{}}]}