{"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,"execution":{"iopub.status.busy":"2021-06-24T14:53:35.203574Z","iopub.execute_input":"2021-06-24T14:53:35.203873Z","iopub.status.idle":"2021-06-24T14:53:37.336941Z","shell.execute_reply.started":"2021-06-24T14:53:35.203823Z","shell.execute_reply":"2021-06-24T14:53:37.335902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"----------------------------------------------\n# Data understanding\n### What data do we have available?\n\n**220k training images and 57k evaluation images.** The dataset is a subset of the [PCam dataset](https://github.com/basveeling/pcam) and the only difference between these two is that all duplicate images have been removed. The PCam dataset 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\nAccording to the data description, there is a 50/50 balance between positive and negative examples in the training and test splits. However, **the training distribution seems to be 60/40 (negatives/positives)**. A positive label means that there is at least one pixel of tumor tissue in the center region (32 x 32px) of the image. **Tumor tissue in the outer region of the patch does not influence the label.** This means that a negatively labeled image could contain metastases in the outer region. Thus, it would be a good idea to crop the images to the center region.\n\n**Image file descriptors**\n\nDescription | \n:--------:|:-------:\nFormat | TIF\nSize | 96 x 96\nChannels | 3\nBits per channel | 8\nData type | Unsigned char\nCompression | Jpeg\n\n### Is the data relevant to the problem?\n\nThis dataset is a combination of two independent datasets collected in Radboud University Medical Center (Nijmegen, the Netherlands), and the University Medical Center Utrecht (Utrecht, the Netherlands). The slides are produced by routine clinical practices and 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.\n\n### Is it valid? Does it reflect our expectations?\n\nAccording to the data description, the dataset has been stripped of duplicates. However, this has not been confirmed by testing.\n\n> For the entire dataset, when the slide-level label was unclear during the inspection of the H&E-stained slide, an additional WSI with a consecutive tissue section, immunohistochemically stained for cytokeratin, was used to confirm the classification.\n- [1399 H&E-stained sentinel lymph node sections of breast cancer patients: the CAMELYON dataset](https://academic.oup.com/gigascience/article/7/6/giy065/5026175)\n\n### Is the data quality, quantity, recency sufficient?\n\n> All glass slides included in the CAMELYON dataset were part of routine clinical care and are thus of diagnostic quality. However, during the acquisition process, scanning can fail or result in out-of-focus images. As a quality-control measure, all slides were inspected manually after scanning. The inspection was performed by an experienced technician (Q.M. and N.S. for UMCU, M.H. or R.vd.L. for the other centers) to assess the quality of the scan; when in doubt, a pathologist was consulted on whether scanning issues might affect diagnosis.\n- [1399 H&E-stained sentinel lymph node sections of breast cancer patients: the CAMELYON dataset](https://academic.oup.com/gigascience/article/7/6/giy065/5026175)","metadata":{}},{"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","execution":{"iopub.status.busy":"2021-06-24T14:53:37.338618Z","iopub.execute_input":"2021-06-24T14:53:37.338887Z","iopub.status.idle":"2021-06-24T14:53:37.909546Z","shell.execute_reply.started":"2021-06-24T14:53:37.338839Z","shell.execute_reply":"2021-06-24T14:53:37.908769Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:37.911054Z","iopub.execute_input":"2021-06-24T14:53:37.911383Z","iopub.status.idle":"2021-06-24T14:53:42.882448Z","shell.execute_reply.started":"2021-06-24T14:53:37.911340Z","shell.execute_reply":"2021-06-24T14:53:42.881598Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:42.883820Z","iopub.execute_input":"2021-06-24T14:53:42.884334Z","iopub.status.idle":"2021-06-24T14:53:42.981422Z","shell.execute_reply.started":"2021-06-24T14:53:42.884282Z","shell.execute_reply":"2021-06-24T14:53:42.980595Z"},"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,"execution":{"iopub.status.busy":"2021-06-24T14:53:42.986945Z","iopub.execute_input":"2021-06-24T14:53:42.989226Z","iopub.status.idle":"2021-06-24T14:53:45.624702Z","shell.execute_reply.started":"2021-06-24T14:53:42.989169Z","shell.execute_reply":"2021-06-24T14:53:45.623905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Classifying metastases is probably not an easy task for a trained pathologist and extremely difficult for an untrained eye. According to [Libre Pathology](https://librepathology.org/wiki/Lymph_node_metastasis), lymph node metastases can have these features:\n\n> - Foreign cell population - key feature (Classic location: subcapsular sinuses)\n- Cells with cytologic features of malignancy\n    - Nuclear pleomorphism (variation in size, shape and staining).\n    - Nuclear atypia:\n        - **Nuclear enlargement**.\n        - **Irregular nuclear membrane**.\n        - **Irregular chromatin pattern, esp. asymmetry**.\n        - **Large or irregular nucleolus**.\n     - Abundant mitotic figures.\n- Cells in architectural arrangements seen in malignancy; highly variable - dependent on tumour type and differentiation.\n    - Gland formation.\n    - Single cells.\n    - Small clusters of cells.\n  \n**The takeaway from this is probably that irregular nuclear shapes, sizes or staining shades can indicate metastases.**\n\n### How is the data best transformed for modeling?\n\nWe know that the label of the image is influenced only by the center region (32 x 32px) so it would make sense to crop our data to that region only. However, some useful information about the surroundings could be lost if we crop too close.  This hypothesis could be confirmed by training models with varying crop sizes. My initial results with 32 x 32px size showed worse performance than with 48 x 48px but I haven't done a search for optimal size.\n\n### How may we increase the data quality?\n\nWe could inspect if the data contains bad data (too unfocused or corrupted) and remove those to increase the overall quality.\n\n### Preprocessing and augmentation\nThere are couple of ways we can use to avoid overfitting; more data, augmentation, regularization and less complex model architectures. Here we will define what image augmentations to use and add them directly to our image loader function. Note that if we apply augmentation here, augmentations will also be applied when we are predicting (inference). This is called test time augmentation (TTA) and it can improve our results if we run inference multiple times for each image and average out the predictions. \n\n**The augmentations we can use for this type of data:**\n- random rotation\n- random crop\n- random flip (horizontal and vertical both)\n- random lighting\n- random zoom (not implemented here)\n- Gaussian blur (not implemented here)\n\nWe will use OpenCV with image operations because in my experience, OpenCV is a lot faster than *PIL* or *scikit-image*.","metadata":{}},{"cell_type":"code","source":"from glob import glob \nimport numpy as np\nimport pandas as pd\nimport keras,cv2,os\n\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Flatten, BatchNormalization, Activation\nfrom keras.layers import Conv2D, MaxPool2D\n\nfrom tqdm import tqdm_notebook,trange\nimport matplotlib.pyplot as plt\n\nimport gc #garbage collection, we need to clean up and save memory\n#set paths to training and test data\npath = \"../input/\" #adapt this path, when running locally\ntrain_path = path + 'train/'\ntest_path = path + 'test/'\n\ndf = pd.DataFrame({'path': glob(os.path.join(train_path,'*.tif'))}) # load the filenames\ndf['id'] = df.path.map(lambda x: x.split('/')[3].split(\".\")[0]) # keep only the file names in 'id'\nlabels = pd.read_csv(path+\"train_labels.csv\") # read the provided labels\ndf = df.merge(labels, on = \"id\") # merge labels and filepaths\ndf.head(10) # print the first three entrys\ndef load_data(N,df):\n    \"\"\" This functions loads N images using the data df\n    \"\"\"\n    # allocate a numpy array for the images (N, 96x96px, 3 channels, values 0 - 255)\n    X = np.zeros([N,96,96,3],dtype=np.uint8) \n    #convert the labels to a numpy array too\n    y = np.squeeze(df.as_matrix(columns=['label']))[0:N]\n    #read images one by one, tdqm notebook displays a progress bar\n    for i, row in tqdm_notebook(df.iterrows(), total=N):\n        if i == N:\n            break\n        X[i] = cv2.imread(row['path'])\n          \n    return X,y\n\ndf.shape\n\n# Load 10k images\nN=10000\nX,y = load_data(N=N,df=df) ","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:07:16.049855Z","iopub.execute_input":"2021-06-24T15:07:16.050212Z","iopub.status.idle":"2021-06-24T15:08:30.749189Z","shell.execute_reply.started":"2021-06-24T15:07:16.050158Z","shell.execute_reply":"2021-06-24T15:08:30.748252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, for most people (like me) there is no easy way to discern which images contain cancer cells. There is a [variety of things](https://www.nature.com/articles/nmeth.4397.pdf) one can look at to get a deeper insight into the data.\n\n## Let's starting looking at the data distribution\n\nWe'll start by looking at how often the classes are represented. ","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(4, 2),dpi=150)\nplt.bar([1,0], [(y==0).sum(), (y==1).sum()]); #plot a bar chart of the label frequency\nplt.xticks([1,0],[\"Negative (N={})\".format((y==0).sum()),\"Positive (N={})\".format((y==1).sum())]);\nplt.ylabel(\"# of samples\")","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:30.754400Z","iopub.execute_input":"2021-06-24T15:08:30.756793Z","iopub.status.idle":"2021-06-24T15:08:30.978308Z","shell.execute_reply.started":"2021-06-24T15:08:30.756738Z","shell.execute_reply":"2021-06-24T15:08:30.977329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, we have about a 60 / 40 split of negative to positive samples in the data. This is important because it means that a trivial classifier, that just labels every sample as negative, would achieve an accuracy of 60%. Possible countermeasures to avoid a bias in the classifier and improve stability during training are, e.g., [over- and undersampling](https://en.wikipedia.org/wiki/Oversampling_and_undersampling_in_data_analysis).\n\n## Let's start looking at each class individually\n\nWe'll now split the data into positive and negative samples to get an idea what makes the classes unique. Such an analysis can often provide insight into possible [feature engineering](https://en.wikipedia.org/wiki/Feature_engineering) possibilities or which data transformations may be useful.","metadata":{}},{"cell_type":"code","source":"positive_samples = X[y == 1]\nnegative_samples = X[y == 0]","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:11:05.933377Z","iopub.execute_input":"2021-06-24T15:11:05.933689Z","iopub.status.idle":"2021-06-24T15:11:06.111171Z","shell.execute_reply.started":"2021-06-24T15:11:05.933637Z","shell.execute_reply":"2021-06-24T15:11:06.110179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we will compare the distribution of pixel values for each color channel (RGB) individually and all channels together","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:30.984439Z","iopub.execute_input":"2021-06-24T15:08:30.987257Z","iopub.status.idle":"2021-06-24T15:08:31.024528Z","shell.execute_reply.started":"2021-06-24T15:08:30.987191Z","shell.execute_reply":"2021-06-24T15:08:31.020866Z"}}},{"cell_type":"code","source":"nr_of_bins = 256 #each possible pixel value will get a bin in the following histograms\nfig,axs = plt.subplots(4,2,sharey=True,figsize=(8,8),dpi=150)\n\n#RGB channels\naxs[0,0].hist(positive_samples[:,:,:,0].flatten(),bins=nr_of_bins,density=True)\naxs[0,1].hist(negative_samples[:,:,:,0].flatten(),bins=nr_of_bins,density=True)\naxs[1,0].hist(positive_samples[:,:,:,1].flatten(),bins=nr_of_bins,density=True)\naxs[1,1].hist(negative_samples[:,:,:,1].flatten(),bins=nr_of_bins,density=True)\naxs[2,0].hist(positive_samples[:,:,:,2].flatten(),bins=nr_of_bins,density=True)\naxs[2,1].hist(negative_samples[:,:,:,2].flatten(),bins=nr_of_bins,density=True)\n\n#All channels\naxs[3,0].hist(positive_samples.flatten(),bins=nr_of_bins,density=True)\naxs[3,1].hist(negative_samples.flatten(),bins=nr_of_bins,density=True)\n\n#Set image labels\naxs[0,0].set_title(\"Positive samples (N =\" + str(positive_samples.shape[0]) + \")\");\naxs[0,1].set_title(\"Negative samples (N =\" + str(negative_samples.shape[0]) + \")\");\naxs[0,1].set_ylabel(\"Red\",rotation='horizontal',labelpad=35,fontsize=12)\naxs[1,1].set_ylabel(\"Green\",rotation='horizontal',labelpad=35,fontsize=12)\naxs[2,1].set_ylabel(\"Blue\",rotation='horizontal',labelpad=35,fontsize=12)\naxs[3,1].set_ylabel(\"RGB\",rotation='horizontal',labelpad=35,fontsize=12)\nfor i in range(4):\n    axs[i,0].set_ylabel(\"Relative frequency\")\naxs[3,0].set_xlabel(\"Pixel value\")\naxs[3,1].set_xlabel(\"Pixel value\")\nfig.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:11:10.086981Z","iopub.execute_input":"2021-06-24T15:11:10.087293Z","iopub.status.idle":"2021-06-24T15:11:25.123022Z","shell.execute_reply.started":"2021-06-24T15:11:10.087242Z","shell.execute_reply":"2021-06-24T15:11:25.122182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So we can already spot differences in the distributions of all channels individually and together between positive and negative samples: \n* Negative samples seem to have higher, i.e. brighter, pixel values in general and especially in the green color channel. \n* Interestingly, the positive samples have a darker green channel than red and blue while this is not true for the negative samples. However, very dark pixels are for both sample sets mostly only present in the green channel.\n* Furthermore, note the relatively high frequency of the pixel value 255. Looking at the data above we can see, that these can likely be attributed to the bright white image regions present in some images. They seem to be present in both positive and negative samples similarly frequently.\n\nNow, let's switch perspective and look at the distribution of mean image brightness, i.e. mean image pixel values. Note, previously we were averaging over all pixel values in the positive and negative samples. Now, we will take the mean of each individual image and look at that distribution.","metadata":{}},{"cell_type":"code","source":"nr_of_bins = 64 #we use a bit fewer bins to get a smoother image\nfig,axs = plt.subplots(1,2,sharey=True, sharex = True, figsize=(8,2),dpi=150)\naxs[0].hist(np.mean(positive_samples,axis=(1,2,3)),bins=nr_of_bins,density=True);\naxs[1].hist(np.mean(negative_samples,axis=(1,2,3)),bins=nr_of_bins,density=True);\naxs[0].set_title(\"Mean brightness, +ve samples\");\naxs[1].set_title(\"Mean brightness, -ve samples\");\naxs[0].set_xlabel(\"Image mean brightness\")\naxs[1].set_xlabel(\"Image mean brightness\")\naxs[0].set_ylabel(\"Relative frequency\")\naxs[1].set_ylabel(\"Relative frequency\");","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:12:03.912911Z","iopub.execute_input":"2021-06-24T15:12:03.913237Z","iopub.status.idle":"2021-06-24T15:12:04.819325Z","shell.execute_reply.started":"2021-06-24T15:12:03.913188Z","shell.execute_reply":"2021-06-24T15:12:04.818463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Once again, we can spot pretty obvious differenes between the positive and negative samples. The distribution of mean brightness for the positive samples looks almost like a normal distribution around a brightness of 150. The negative samples, however, seem to follow some bimodal distribution with peaks around 140 and 225. \n\n**Conclusions:** \n* There are some easily spotted differences in the distributions of pixel values and mean image brightness between positive and negative samples. This is good, because whatever model we will use can likely use this. \n* Some of the images seem to contain very bright regions, which are likely artifacts of the recording process. We might have to find a way to deal with them. They are almost equally distributed between positive and negative samples and, hence, probably not easily usable as a feature.\n* We have about 50% more negative than positive samples. This might require adjustments.","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2021-06-24T15:08:31.035672Z","iopub.status.idle":"2021-06-24T15:08:31.038191Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:45.712014Z","iopub.status.idle":"2021-06-24T14:53:45.712682Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:45.713960Z","iopub.status.idle":"2021-06-24T14:53:45.714692Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:45.715979Z","iopub.status.idle":"2021-06-24T14:53:45.716739Z"},"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":{"execution":{"iopub.status.busy":"2021-06-24T14:53:45.718014Z","iopub.status.idle":"2021-06-24T14:53:45.718742Z"},"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":{}}]}