{"cells":[{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Introduction to Problem and Notebook\n\nOver the last 50 years, several scoring systems have been proposed that allow pathologists to grade tumours based on their appearance. In addition to inter-pathologist variability, some systems which require quantitative analysis, are too time-consuming to use in a routine clinical setting. Automation of this process through machine learning and computer vision is therefore potentially a high-impact topic of research.\n\nFrom [Kaggle](\"https://www.kaggle.com/c/histopathologic-cancer-detection/overview\"): \"In this competition, you must create an algorithm to identify metastatic cancer in small image patches taken from larger digital pathology scans. The data for this competition is a slightly modified version of the PatchCamelyon (PCam) benchmark dataset.\"\n\nThe main influences on this notebook are certainly the [Joni Juvonen notebook](\"https://www.kaggle.com/qitvision/a-complete-ml-pipeline-fast-ai\") entry using Fast AI and the [Antonio de Perio blog post](\"https://towardsdatascience.com/histopathological-cancer-detection-with-deep-neural-networks-3399be879671\") on Towards Data Science. \n\nSeveral other excellent resources were also useful for direction. I shall cite these as we go and also summarise at the end of the notebook.\n\nAt a high-level this notebook is organised as follows:\n* Background on how we approach the challenge and information relevant to the problem\n* Exploratory data analysis\n* Modelling the classification problem\n* Model performance\n\n## Notes on Approaching the Challenge\n\nWe began by information gathering. This was certainly a useful step in order to confront some of the author's initial biases; namely that a CNN coded from scratch CNN would certainly be the only solution. In addition, researching the problem benefits us with context, clarity, and insights into the data. \n\nThe notebook is built on the following general approach:\n* Information gathering\n* Split the problem into smaller parts (see \"Plan Overview\" at the end of the notebook)\n* Execute\n* Refine plan as more information and obstacles become apparent\n\nThe review paper [Deep neural network models for computational histopathology: A survey](\"https://arxiv.org/pdf/1912.12378.pdf\") and article [Artificial Intelligence in Pathology](\"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6344799/\") were extremely helpful in understanding that there is a consensus on the machine learning approach to classification of PCam. In particular, they highlight patch-level CNNs and transfer-learning as the preffered approach. We also take this approach here:\n* \"Methods based on CNNs have consistently outperformed other handcrafted methods in a variety of deep learning\"\n* \"The most popular and widely adopted technique in digital pathology is the use of transfer learning approach\"\n\nWe discuss CNNs, transfer-learning, and the PCam data set in more detail later in the notebook. \n\nIn analysing model performance we pay close attention to false negatives, i.e., when metastastic cancer exists but is not detected. This seems to be an important statistic. According to [Saba 2020](\"https://reader.elsevier.com/reader/sd/pii/S1876034120305633?token=104840BB146B11F0239FAFE45FC879C93B6D4A7CDDED333EADDB2F5FA1BDE8BA672FD2C0C22DB5753011ED959877D131\") GLOBOCAN statistics show that 18.1m new cancer cases emerged in 2018 that caused 9.6m deaths. If we take this (approx. 50%) death rate and apply the false-negative rate 3% the trained model in [our Fast AI notebook](\"https://www.kaggle.com/qitvision/a-complete-ml-pipeline-fast-ai\") we are potentially missing 299k untreated deaths. Unfortunately [Shapley Values](\"https://towardsdatascience.com/making-sense-of-shapley-values-dc67a8e4c5e8\") are not available for deep-learning models but Fast AI has suggested visual attention maps to address interpretability.\n\n## Useful Terminology\n\nThe following definitions and terms helped provide clarity and make swifter modelling choices:\n* Metastasis is the spread of cancer from the point of onset to different parts of the body. Metastases most commonly develop when cancer cells break away from the main tumor and enter the bloodstream or lymphatic system. These systems carry fluids around the body\n* TIF files contain an image saved in the Tagged Image File Format (TIFF), a high-quality graphics format\n* Whole slide images (WSIs) are high quality digital scans of glass slides used in pathology\n* Patch-based learning splits the WSI into smaller (non-overlapping) patches which are used as training records for a machine-learning model. After the patch-level CNN is trained, another ML model is often developed for the whole image level decision\n\n## Introduction to Data\n\nPCam (PatchCamelyon) is a binary classification image dataset containing labeled low-resolution images of lymph node sections extracted from digital histopathological scans in tif format. A positive label indicates that the center 32x32px region of a patch contains at least one pixel of tumor tissue.\n\nPCam is derived from the Camelyon16 Challenge which contains 400 H&E stained WSIs of sentinel lymph node sections. The process of patch production is described in its [creator's post](\"http://basveeling.nl/posts/pcam/\") and the associated git repository. PCam was created specifically for machine learning reasearches to be able to approach the problem of cancer detection as easily as they would the MNIST or CIFAR10.\n\n## Initial Data Load\n\nWe will load data directly from the Kaggle directory. We will load only the training-data from the Kaggle competition as this is labelled and can be used for training our model and assessing performance."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport os\nimport gc\nimport cv2\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nimport random\nfrom fastai import *\nfrom fastai.vision import *\nfrom torchvision.models import resnet50\nfrom PIL import Image\nfrom sklearn.utils import shuffle\nfrom sklearn.decomposition import PCA\nfrom sklearn.model_selection import train_test_split\nfrom skimage.restoration import estimate_sigma","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = pd.read_csv('/kaggle/input/train_labels.csv')\ntrain_path = '/kaggle/input/train'\n# quick look at the label stats\nprint(\"Value counts of data labels:\")\nprint(data['label'].value_counts())\nprint(\"\\nSample data labels:\")\nprint(data.head())","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"We have 131k negative and 89k positive labels leading to an imbalance of 60%-40% indexed by the image ID. This does not seem imbalanced enough to warrant methods such as weighting our positive class, subsampling, or oversampling.\n\nAs it helps to visualise the data we will load 5 random samples of the positive and negative classes. We highlight the 32x32px region where we must search for evidence of metastastic cancer."},{"metadata":{"trusted":true,"_uuid":"90811cbf83e65ba2679ed1044798a71b2dd89818","_kg_hide-input":true},"cell_type":"code","source":"# As per qitvision reverse the cv2 colours from bgr to rgb - we will then see the same images as directly from PCam;\n# plus I wouldn't want this to affect a pretrained model's ability\n\ndef readImage(path):\n    # OpenCV reads the image in bgr format by default\n    bgr_img = cv2.imread(path)\n    # We flip it to rgb for visualization purposes\n    b,g,r = cv2.split(bgr_img)\n    rgb_img = cv2.merge([r,g,b])\n    return rgb_img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# random sampling\nshuffled_data = shuffle(data, random_state = 27)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"20e119981333d59b0ba4a52b46a9aa5b4bb2436f","_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,5, figsize=(20,8))\nfig.suptitle('Histopathologic scans of lymph node sections',fontsize=20)\n# Negatives\nfor i, idx in enumerate(shuffled_data[shuffled_data['label'] == 0]['id'][:5]):\n    path = os.path.join(train_path, idx)\n    ax[0,i].imshow(readImage(path + '.tif'))\n    # Create a Rectangle patch\n    box = patches.Rectangle((32,32),32,32,linewidth=4,edgecolor='b',facecolor='none', linestyle=':', capstyle='round')\n    ax[0,i].add_patch(box)\nax[0,0].set_ylabel('Negative samples', size='large')\n# Positives\nfor i, idx in enumerate(shuffled_data[shuffled_data['label'] == 1]['id'][:5]):\n    path = os.path.join(train_path, idx)\n    ax[1,i].imshow(readImage(path + '.tif'))\n    # Create a Rectangle patch\n    box = patches.Rectangle((32,32),32,32,linewidth=4,edgecolor='r',facecolor='none', linestyle=':', capstyle='round')\n    ax[1,i].add_patch(box)\nax[1,0].set_ylabel('Tumor tissue samples', size='large')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Exploratory Data Analysis\n\nWe will now begin to examine the data quality and the high-level properties of the data at the class level.\n\nThe following code:\n* Records all images which are dark (larges pixel value < 10)\n* Records all images which are light (lowest pixel value > 245)\n* Records a robust wavelet-based estimator of the (Gaussian) noise standard deviation\n* Checks if any images are loaded with the wrong data type (np.int8)\n\nDue to time constraints we only load 10k images. Testing has shown that the conclusions/insights scale to more data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# As we count the statistics, we can check if there are any completely black or white images\ndark_th = 10      # If no pixel reaches this threshold, image is considered too dark\nbright_th = 245   # If no pixel is under this threshold, image is considerd too bright\ntoo_dark_idx = []\ntoo_bright_idx = []\nbad_dtypes = []\nnoise_array_pos = np.array([])\nnoise_array_neg = np.array([])\nmeans_pos = np.array([])\nmeans_neg = np.array([])\nstds_pos = np.array([])\nstds_neg = np.array([])\n\nN = 10000 # max length of positive and negative samples\npositive_samples = []\nnegative_samples = []\n\niterable = shuffled_data[\"id\"].head(10000)\n#for i, idx in tqdm(enumerate(iterable), 'computing statistics...(220025 it total)'):\nfor i, idx in enumerate(iterable): \n    # What is the label\n    label = shuffled_data.loc[shuffled_data[\"id\"] == idx, \"label\"].values[0]\n    \n    # Read image\n    path = os.path.join(train_path, idx)\n    img = readImage(path + '.tif')\n    imagearray = img.reshape(-1,3)\n    \n    # Check for anomylous data types:\n    if img.dtype is not np.dtype('uint8'):\n        bad_dtypes.append(idx)\n        \n    # Build separate positive and negative samples for examination later\n    if label == 1 and len(positive_samples) < N:\n        positive_samples.append(img)\n        # Record \"noise\" of image\n        noise_array_pos = np.append(noise_array_pos, estimate_sigma(img, multichannel=True, average_sigmas=True))\n        # Record mean and std of each image\n        means_pos = np.append(means_pos, np.mean(imagearray))\n        stds_pos = np.append(stds_pos, np.std(imagearray))\n    if label == 0 and len(negative_samples) < N:\n        negative_samples.append(img)\n        # Record \"noise\" of image\n        noise_array_neg = np.append(noise_array_neg, estimate_sigma(img, multichannel=True, average_sigmas=True))\n        # Record mean and std of each image\n        means_neg = np.append(means_neg, np.mean(imagearray))\n        stds_neg = np.append(stds_neg, np.std(imagearray))\n        \n    # is this too dark\n    if(imagearray.max() < dark_th):\n        too_dark_idx.append(idx)\n        continue # do not include in statistics\n    \n    # is this too bright\n    if(imagearray.min() > bright_th):\n        too_bright_idx.append(idx)\n        continue # do not include in statistics","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"positive_samples = np.array(positive_samples)\nnegative_samples = np.array(negative_samples)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"# of samples which are too dark: \" + str(len(too_dark_idx)))\nprint(\"# of samples which are too light: \" + str(len(too_bright_idx)))\nprint(\"# of samples which have the wrong dtype: \" + str(len(bad_dtypes)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that there are no samples which do not have int dtypes. In addition we can only find a minority of samples which are either very light or very dark. As per [qitvision](\"https://www.kaggle.com/qitvision/a-complete-ml-pipeline-fast-ai\") we do not expect these anomolies to affect results.\n\nThe following code plots the distribution of RGB values for positive and negative classes"},{"metadata":{"trusted":true},"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()\n\nnr_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, positive samples\");\naxs[1].set_title(\"Mean brightness, negative 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\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In agreement with [gomezp](\"https://www.kaggle.com/gomezp/complete-beginner-s-guide-eda-keras-lb-0-93\") (although here we flip RGB accoring to OpenCV) we can see that:\n* Negative samples are bimodel in terms of brightness (one bright peak and one dark prak)\n* Negative samples have skewed to darker RGB channels in general\n* Positive samples have a lighter G channel than RB\n\nNext we plot the \"average\" positive and negative class pictures where positive images should appear less green and brighter in general."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot average images\nfig,axs = plt.subplots(1,2,sharey=True, sharex = True, figsize=(8,2),dpi=150)\naxs[0].imshow(Image.fromarray(np.mean(positive_samples,axis=(0)), 'RGB'))\naxs[1].imshow(Image.fromarray(np.mean(negative_samples,axis=(0)), 'RGB'))\naxs[0].set_title(\"Average image, positive samples\");\naxs[1].set_title(\"Average image, negative samples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We also examine the correlation between different features. Initially we have 96x96x3 = 27,648 features per image which is too large to process quickly. We settle for correlations of (i) the RGB means of each pixel for (ii) the center of each patch. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(1,2,sharey=True, sharex = True, figsize=(8,2),dpi=150)\naxs[0].set_title(\"Correlation of pixels, positives\");\naxs[1].set_title(\"Correlation of pixels, negatives\");\ntemp = np.mean(positive_samples[:,32:64,32:64,:], axis = 3).flatten().reshape(positive_samples.shape[0], 32*32)\ntemp = pd.DataFrame(data = temp, columns = [str(i) for i in range(temp.shape[1])])\nsns.heatmap(temp.corr(), ax = axs[0])\ntemp = np.mean(negative_samples[:,32:64,32:64,:], axis = 3).flatten().reshape(negative_samples.shape[0], 32*32)\ntemp = pd.DataFrame(data = temp, columns = [str(i) for i in range(temp.shape[1])])\nsns.heatmap(temp.corr(), ax = axs[1])\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here positive and negative samples look very similar to the eye.\n\nNext, we show the scatter plots of means vs standard deviations of pixel values (across RGB)."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(1,2,sharey=True, sharex = True, figsize=(8,2),dpi=150)\naxs[0].set_title(\"Mean vs Std, positives\");\naxs[1].set_title(\"Mean vs Std, negatives\");\naxs[0].scatter(means_pos, stds_pos)\naxs[1].scatter(means_neg, stds_neg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is a lot of overlap in the distributions. Positives seem to come from a smaller sample - which is natural given we have far less observations here. Negative samples seem to have an order of magnitude more outliers. One interesting non-overlapping section is the low/mid-mean/mid-std region that only positive samples seem to occupy. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Free memory\ntemp = None\naxs = None\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### A Note on Noise\n\nWe plot the Gaussian noise distribution estimated via skimage.restoration.estimate_sigma. We see a bimodal distribution for noise with the main peak close to zero. In general lower numbers mean more noise."},{"metadata":{"trusted":true},"cell_type":"code","source":"nr_of_bins = 100 #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(noise_array_pos.flatten(),bins=nr_of_bins,density=True);\naxs[1].hist(noise_array_neg.flatten(),bins=nr_of_bins,density=True);\naxs[0].set_title(\"Noise estimation of positives\");\naxs[1].set_title(\"Noise estimation of negatives\");\nfor i in [0, 1]:\n    axs[i].set_xlabel(\"Noise estimation\")\n    axs[i].set_ylabel(\"Relative frequency\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Gaussian noise distribution is estimated via skimage.restoration.estimate_sigma. We see a bimodal distribution for negative-class noise with the main peak close to zero. The low peak may warrent a noise cancelling method such as [wavelet denoising](\"https://scikit-image.org/docs/dev/auto_examples/filters/plot_denoise_wavelet.html#sphx-glr-auto-examples-filters-plot-denoise-wavelet-py\").\n\nWe can also study the PCA of the images. The below plots show us that total variance in pixel observations (RGB across 96x96 pixels) explained by the first 2000 principal components is around 90% for both positive and negative samples. A higher variance seems to be explained by fewer principal components for negative samples in general which ties in with the idea that negative samples are more \"noisy\"."},{"metadata":{"trusted":true},"cell_type":"code","source":"pca = PCA(700)\nfig, axs = plt.subplots(1,2,sharey=True, sharex = True, figsize=(8,2),dpi=150)\naxs[0].set_title(\"PCA variance, positive samples\");\naxs[0].set_xlabel(\"number of components\")\naxs[0].set_ylabel(\"cumulative explained variance\")\naxs[1].set_title(\"PCA variance, negative samples\");\naxs[1].set_xlabel(\"number of components\")\naxs[1].set_ylabel(\"cumulative explained variance\")\npca.fit(positive_samples.flatten().reshape(positive_samples.shape[0], 96*96*3))\naxs[0].plot(np.cumsum(pca.explained_variance_ratio_))\npca.fit(negative_samples.flatten().reshape(negative_samples.shape[0], 96*96*3))\naxs[1].plot(np.cumsum(pca.explained_variance_ratio_))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Free memory\npca = None\npositive_samples, negative_samples = None, None\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Processing\n\nFor image-data processing we use ImageDataBunch (PyTorch) directly as in [Antonio de Perio's notebook.](\"https://github.com/humanunsupervised/humanunsupervised.com/blob/master/pcam/pcam-cancer-detection.ipynb\"). This is over other methods such as ImageDataGenerator in Keras.\n\nScaling, renormalisation, centering, and even data augmentation and batching is done for free by ImageDataBunch through the function get_transforms. Normalisation is done using imagenet_stats."},{"metadata":{"trusted":true},"cell_type":"code","source":"tfms = get_transforms(do_flip=True)\n\nbs=64 # also the default batch size\ndata = ImageDataBunch.from_csv(\n    '/kaggle/input/', \n    ds_tfms=tfms, \n    size=224, \n    suffix=\".tif\",\n    folder=\"train\", \n    test=\"test\",\n    csv_labels=\"train_labels.csv\", \n    bs=bs)\n\ndata.normalize(imagenet_stats)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model Choise\n\nResults from previous sections indicate that (i) there are some distributional differences between positive and negative classes, (ii) many features are required to explain differences, (iii) noise (and therefore overfitting) may be an issue. \n\nWe will work with resnet50 due to the successful implementation [here](\"https://github.com/humanunsupervised/humanunsupervised.com/blob/master/pcam/pcam-cancer-detection.ipynb\") and the glowing reviews [here](\"https://arxiv.org/pdf/1912.12378.pdf\") and [here](\"https://towardsdatascience.com/an-overview-of-resnet-and-its-variants-5281e2f56035\").\n\nResnet50 is a 50-layer deep convolutional neural network (CNN). It can handle large amounts of data and features. In addition, using ImageDataBunch for data augmentation and regularisation we can combat overfitting.\n\n### Quick Notes on CNNs\n\nCNNs general start with several layers of data-processing inorder to reduce noise and increase the visibility of important features which are easier to learn by neural networks. The important components of a CNN can be summarised as follows:\n* Convolutional Layer; This layer applies some filters to an image to remove low and high level features\n* Stride; Is the value which indicates that the filter will shift the image in increments of \"Stride\" pixels (steps)\n* Padding; After convolution, to manage the size difference between the input and outputs\n* Activation Layer; Non-linearity through activation function\n* Pooling Layer; Reduces the offset between the new representation and the number of parameters within the neural network. Generally creates smaller outputs that contain enough information\n* Flattening Layer; Data to 1d for the Fully-Connected Layer\n* Fully-Connected Layer; The deep NN to learn the features and make predictions\n\n### Training\nWe will train via transfer learning. We will specifically use the one [cycle policy](https://arxiv.org/abs/1803.09820). As per [qitvision](\"https://www.kaggle.com/qitvision/a-complete-ml-pipeline-fast-ai\"), the fastai library has implemented a training function for one cycle policy that we can use with only a few lines of code.\n\nThe method relies on finding the optimal `learning rate` and `weight decay` (L2 penalty) values. The optimal lr is just before the base of the loss and before the start of divergence. It is important that the loss is still descending where we select the learning rate.\n\nIn the following code we use the qitvision optimal learning rates and weight decay. This is because of the time-cost of an optimal search and the assumption that the data set and the modelling framework is equivalent.\n\nFor the readers information we plot the evolution of the learning rate and the evolution of the loss function. The results show that the lr increases and then decreasese. In addition the loss functions of the training and validation sets seem to smoothly decrease and do not diverge."},{"metadata":{"trusted":true},"cell_type":"code","source":"def getLearner():\n    return create_cnn(data, models.resnet50, pretrained=True, path='.', metrics=error_rate, ps=0.5, callback_fns=ShowGraph)\n\nlearner = getLearner()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max_lr = 2e-2\nwd = 1e-4\n\n# 1cycle policy\nlearner.fit_one_cycle(cyc_len=5, max_lr=max_lr, wd=wd)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds, y, loss = learner.get_preds(with_loss=True)\n# get accuracy\nacc = accuracy(preds, y)\nprint('The accuracy is {0} %.'.format(acc*100))\n# plot learning rate of the one cycle\nlearner.recorder.plot_lr()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learner.save('./pcam_resnet50_frozen')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that the final accuracy is impressive. However, this is using resnet50 weights trained to another data set and we may acheive better performance by retraining.\n\n## Finetuning\n\nFinetuning involves unfreezing the pretrained model, lowering the learning rate, and continuing training. This has proven to improve results in computer vision. \n\nThe following code executes this. We can also see the loss function evolving over training. There is a divergence of training and validation performance indicating some overfitting to the training set."},{"metadata":{"trusted":true},"cell_type":"code","source":"# load the baseline model\nlearner.load('./pcam_resnet50_frozen')\n\n# unfreeze\nlearner.unfreeze()\n\n# Fit new model with lower learning rates\nlearner.fit_one_cycle(cyc_len=5, max_lr=slice(4e-5,4e-4))\nlearner.recorder.plot_losses()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"learner.save('./pcam_resnet50_finetuned')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model Performance\n\nThe following code gives us the final performance of the model in terms of accuracy, the confusion matrix, the area under curve, and the precision/recall curve. In addition, we study examples of false negatives using fastai's ClassificationInterpretation. Here areas of the images important to classification are highlighted as a heatmap."},{"metadata":{"trusted":true},"cell_type":"code","source":"preds,y, loss = learner.get_preds(with_loss=True)\n# get accuracy\nacc = accuracy(preds, y)\nprint('The accuracy is {0} %.'.format(acc))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Confusion matrix\ninterp = ClassificationInterpretation.from_learner(learner)\ninterp.plot_confusion_matrix(title='Confusion matrix')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import roc_curve, auc\n# probs from log preds\nprobs = np.exp(preds[:,1])\n# Compute ROC curve\nfpr, tpr, thresholds = roc_curve(y, probs, pos_label=1)\n\n# Compute ROC area\nroc_auc = auc(fpr, tpr)\nprint('ROC area is {0}'.format(roc_auc))\n\nplt.figure()\nplt.plot(fpr, tpr, color='darkorange', label='ROC curve (area = %0.2f)' % roc_auc)\nplt.plot([0, 1], [0, 1], color='navy', linestyle='--')\nplt.xlim([-0.01, 1.0])\nplt.ylim([0.0, 1.01])\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.title('Receiver operating characteristic')\nplt.legend(loc=\"lower right\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"losses,idxs = interp.top_losses()\ninterp.plot_top_losses(16, figsize=(16,16))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def pr_curve(preds, y):\n    \"\"\"\n    Function to create precision recall curve\n    \n    Inputs\n    ------\n    preds - {tf} [prob of negative class, prob of positive class] for each observation\n    y - {np.array} actual class value\n    \n    Outputs\n    -------\n    recall - {np.array} recall for different probability thresholds\n    precision - {np.array} precision for different probability thresholds\n    \"\"\"\n    temp_df = pd.DataFrame(data = np.array(preds), columns = [\"negative_prob\", \"positive_prob\"])\n    temp_df[\"actual\"] = y\n    \n    precision = np.array([])\n    recall = np.array([])\n    for threshold in range(0, 100, 10):\n        threshold /= 100.\n        temp_df[\"predicted_class\"] = temp_df[\"positive_prob\"] > threshold\n        temp_df[\"predicted_class\"] = temp_df[\"predicted_class\"].astype(int)\n        temp_df[\"true_positives\"] = temp_df.apply(lambda x: 1 if x[\"actual\"] == 1 and x[\"predicted_class\"] == 1 else 0, axis = 1)\n        temp_df[\"false_positives\"] = temp_df.apply(lambda x: 1 if x[\"actual\"] == 0 and x[\"predicted_class\"] == 1 else 0, axis = 1)\n        temp_df[\"true_negatives\"] = temp_df.apply(lambda x: 1 if x[\"actual\"] == 0 and x[\"predicted_class\"] == 0 else 0, axis = 1)\n        temp_df[\"false_negatives\"] = temp_df.apply(lambda x: 1 if x[\"actual\"] == 1 and x[\"predicted_class\"] == 0 else 0, axis = 1)\n        \n        p = temp_df[\"true_positives\"].sum() / float(temp_df[\"true_positives\"].sum() + temp_df[\"false_positives\"].sum())\n        r = temp_df[\"true_positives\"].sum() / float(temp_df[\"true_positives\"].sum() + temp_df[\"false_negatives\"].sum())\n        \n        precision = np.append(precision, p)\n        recall = np.append(recall, r)\n    \n    return recall, precision","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"recall, precision = pr_curve(preds, y)\n\nfig,axs = plt.subplots(1,1,sharey=True,figsize=(3,3),dpi=150)\naxs.set_title(\"P/R Curve for Positive Class\");\naxs.set_xlabel(\"Recall\")\naxs.set_ylabel(\"Precision\")\n\n#RGB channels\naxs.plot(recall, precision)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Interpretation\n* There may be overfitting after fine tuning\n* Fine tuning does improve results on the validation set\n* The AUC and precision/recall results look positive but there is room for improvement\n* False negatives with features highlighted show us that whitespace (potentially the brightness) may be leading to missclassifcation\n\nWe should most likely submit the transfer-learned model before fine tuning.\n\n## Potential Next Steps\n* We still do not know the model performance on unseen (testing data)\n* Noise reduction and more regularisation may help with overfitting\n* More methodical tuning over epochs and hyper-parameters would certainly improve performance\n* More features such as chemistry and phenotypic information could be used to help with the classification\n* Due to noise we may be susceptable to adversarial examples"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Unforeseen Obstacles\n\nSome challenges that I had to overcome in creating this notebook:\n* Enabling GPU for kaggle\n* Enabling an environment with fastai v1.0\n* Being unable to save data in \\input\\\n\nFor reference, here is the original plan to execute on the notebook:\n1. Load data\n* Plot samples and highlight 32x32 centre\n* Look for nulls and duplicates \n* Examine data balance --> over/under sampling is not needed\n* Examine RGB distribution\n* Plot mean vs sd for RGB --> look for outliers and anomalies\n* Correlation of pixles\n* PCA and number of important features --> may be needed for feature selection\n* PCA for positives vs negatives\n* Estimate noise --> may need adversarial examples if noise is higher (1/255 precision of images)\n* Normalise data\n* Create crops and augmented images\n* Choose model - pretrained NN - ResNet50\n* Image clustering on ResNet50 features\n* Shuffle date + split into training testing\n* Transfer learning --> use callback to save the best weights\n* Regularisation\n* Fine tunining\n* Examine loss function for overfitting\n* Examine precision/recall matrix\n* Investigate P/R curve\n* Examine permutation importance of false negatives --> how are these different from true negatives\n* Summarise references"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## References\n\n### Blogs\n1. https://towardsdatascience.com/histopathological-cancer-detection-with-deep-neural-networks-3399be879671\n    - Jupyter notebook uses fast.ai to augment PCam and transfer-learn/fine-tune a Resnet50 model. Further regularisation techniques are used.\n2. https://towardsdatascience.com/metastasis-detection-using-cnns-transfer-learning-and-data-augmentation-684761347b59\n    - Transfer learns with NasNet on augmented PCam. Predicts on augmented images and averages output.\n3. https://towardsdatascience.com/an-overview-of-resnet-and-its-variants-5281e2f56035\n    - The core idea of ResNet is introducing a so-called “identity shortcut connection” that skips one or more layers. \n4. Introduced serveral papers that propose interesting variants of ResNet or provide insightful interpretation of it.\n    - It is quite clear why removing a couple of layers in a ResNet architecture doesn’t compromise its performance too much — the architecture \n    - has many independent effective paths and the majority of them remain intact after we remove a couple of layers. \n5. http://basveeling.nl/posts/pcam/\n    - Datasets such as Camelyon require developing intricate dataloaders, that take care of balancing different types of tissue, \n    - and can efficiently work with the terabytes of data available.\n    - Packs the clinically-relevant task of metastasis detection into a straight-forward binary image classification task, akin to CIFAR-10 and MNIST.\n6. https://medium.com/bbm406f19/week-1-histopathological-cancer-detection-30649d8dd847\n    - 4 week project to classify PCam images using Vgg16. References Antony de Perio.\n7. https://towardsdatascience.com/feature-importance-with-neural-network-346eb6205743\n    - Permutation importance is calculated after a model has been fitted\n    - If I randomly shuffle a single feature in the data, leaving the target and all others in place, how would that affect the final prediction performances?\n8. https://geertlitjens.nl/post/getting-started-with-camelyon/\n    - Test set accuracy is 0.8563 with Vgg16\n\n### Kaggle\n1. https://www.kaggle.com/emackie/a-complete-ml-pipeline-fast-ai/edit\n2. https://www.kaggle.com/emackie/complete-beginner-s-guide-eda-keras-lb-0-93/edit\n    - Good comparison of positive an negative samples at EDA stage\n    - Builds and trains CNN ROC fast model\n3. https://www.kaggle.com/artgor/simple-eda-and-model-in-pytorch\n    - Defines their own 3-level CNN\n    - Then goes to ResNet\n4. https://www.kaggle.com/fmarazzi/baseline-keras-cnn-roc-fast-10min-0-925-lb\n    - Loads data in batches to avoid memory error\n    - Augments images with ImageDataGenerator\n    - Defines own 3-level CNN\n    - Does normalisation\n    - Predicts on testing set directly (not on augmented images + averaging)\n5. https://www.kaggle.com/cc786537662/initial-eda-image-processing\n    - Interesting EDA for image processing\n\n### Papers\n1. https://reader.elsevier.com/reader/sd/pii/S1876034120305633?token=ABA60CEBDF6831129168852D8363149C13690BB0BDBAFC23E032ED6FE278A1836BA034EEA21BE7726614DB65F2AC5EF1\n    - Systematic review of current techniques in diag-nosis and cure of several cancers affecting human body badly. Not including PCam\n    - According toGLOBOCAN statistics, about 18.1 million new cancer cases emergedin 2018 that caused 9.6 million cancer deaths\n    - Contains a brief description and comparison of currentreported techniques, results, datasets for cancer\n2. https://arxiv.org/pdf/1912.12378.pdf\n    - Survey of over 130 papers.\n3. https://arxiv.org/pdf/1912.03609v3.pdf\n    - Investigate to what extent alternative variants of Artificial Neural Networks (ANNs) are susceptible to adversarial attacks\n    - On PCam\n    - Stochastic ANNs are more robust but this is more prevalent in MNIST\n4. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6344799/\n    - The history of artificial intelligence in the medical domain, recent advances in artificial intelligence applied to pathology, and future prospects\n    - Normalization is one of the techniques used to reduce digital immaging\n5. https://arxiv.org/abs/1806.03962\n    - Propose a new model for digital pathology segmentation, based on the observation that histopathology images are inherently symmetric under rotation and reflection\n6. https://arxiv.org/pdf/1905.03151.pdf\n    - This paper advocates against permute-and-predict methods for interpreting black box functions\n    - Forces the original model to extrapolate to regions where there is little to no data\n    - Recommends explicitly removing features, conditional permutations, or model distillation methods\n7. https://www.nature.com/articles/nmeth.4397.pdf\n    - Image-based cell profiling is a high-throughput strategy for the quantification of phenotypic differences among a variety of cell populations\n    - Interesting EDA\n\n### Other references\n1. https://datascience.stackexchange.com/questions/29223/exploratory-data-analysis-with-image-datset\n    - Standard practices for EDA with image data\n2. https://keras.io/guides/transfer_learning/\n    - Transfer learning and fine tunining in Keras. Explanation of trainable = False.\n3. https://stackoverflow.com/questions/2440504/noise-estimation-noise-measurement-in-image\n    - Scikit Image has an noise estimate sigma function which works with colour images"},{"metadata":{"trusted":true,"_uuid":"c6d6eefdbeba2be86625abce612a52e7ac19375d"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}