{"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":"# Summary of the code!\n\n### I start by importing a pre-trained ResNet-50, freezing the weights, and adding a few extra fully connected convolutional layers with dropout to the end. The model takes individual scans of breasts as input and predicts whether the output is cancer or not. \n\n### Each scan is treated individually and not compared with any other view from the same breast or scans from the other breast (which is non-ideal given the amount of information shared between breast scans - ideally the network should be modified to accommodate two (or four) breast images). \n\n### For preprocessing I use histogram normalization, center cropping, image resizing, and conversion of single channel greyscale to triple channel RBG image. I am using a focal loss function with alpha (weight of positive samples to correct imbalanced dataset (only 2% positive samples)) of about 0.02 and gamma (weight of recall vs precision) of 0.5. (If I oversample the positive examples or do data augmentation to balance the dataset then I can modify this to gamma=0, which makes the focal loss equivalent to the binary cross entropy loss). \n\n### I set a learning rate of 0.001 and momentum of 0.9, which are default values for the stochastic gradient descent optimizer. \n\n### I split the data into train:validate 80:20 and set up a batch size of 32 for training. I am using only a percentage of the data (eg 10%) to train due to limited free GPU on Kaggle, so 8% train : 2% validate. I also am testing on a different 10% of the dataset.\n\n### I train for about 200 epochs with a patience of 20 epochs of no change in the loss for early stopping (ideally we train for >500 epochs given the loss seems to still have yet to flatline at 200). \n\n### Also, every 2 epochs we unfreeze one more layer from the ResNet-50 to fine tune the model to our dataset. \n\n### Lastly, I plot our loss and probabilistic F1-score on the validation data vs epochs, test on the held-out data that wasn't used for training in train_images, and create a submission on the hidden test dataset.","metadata":{}},{"cell_type":"markdown","source":"Todo: \n\n- make a submission by training on the full dataset -- first figure out how much (% data, epochs) you can train in 20 hours. Then train the model for 20 hours and DOWNLOAD AND SAVE THE MODEL FROM THE WORKING DIRECTORY. Then upload the trained model as a dataset to a new notebook which will be JUST for submissions. Cross your fingers and submit! (Ensure kaggle.json is in the location ~/.kaggle/kaggle.json to use the API.)\n- preprocessing method (center crop, print out images to ensure breast is cropped and correctly leveled, histogram normalization) - print out images!\n- exclude all images that are not CC or MLO\n- check VOI LUT or other examples\n- hyperparameter tuning experiment on 1% or less of the data, including learning rate, momentum, alpha and gamma of the focal loss function, adding/removing CNN output layers to the architecture\n- add other datasets!\n- try automated data-augmentation methods\n- try VGG16, VGG19, ResNet-34 on whole image\n- to handle multiple views, try making a double NN with two outputs, then one final set of FC layers connecting the two outputs initialized to an OR logic gate\n- to handle multiple breasts, try combining the multi-view output from the L breast with that of the R breast using a final set of FC layers initialized to some modified OR logic gate\n- try patch classifier\n- upscale GPU and train final model on full dataset no validation or test split >500 epochs with oversampled positives","metadata":{}},{"cell_type":"code","source":"!virtualenv kaggle_env\n!source kaggle_env/bin/activate\n# !pip install dicomsdl\n# !pip install gdcm\n# !pip install -U pylibjpeg pylibjpeg-openjpeg pylibjpeg-libjpeg\n!pip install /kaggle/input/rsna-2023-whl/dicomsdl-0.109.1-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl","metadata":{"execution":{"iopub.status.busy":"2023-04-18T02:52:25.997301Z","iopub.execute_input":"2023-04-18T02:52:25.998264Z","iopub.status.idle":"2023-04-18T02:52:42.375103Z","shell.execute_reply.started":"2023-04-18T02:52:25.99822Z","shell.execute_reply":"2023-04-18T02:52:42.373453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !pip install pydicom","metadata":{"execution":{"iopub.status.busy":"2023-04-18T02:52:42.379291Z","iopub.execute_input":"2023-04-18T02:52:42.3814Z","iopub.status.idle":"2023-04-18T02:52:42.387327Z","shell.execute_reply.started":"2023-04-18T02:52:42.381345Z","shell.execute_reply":"2023-04-18T02:52:42.386013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil\n\n!mkdir -p /root/.cache/torch/hub/checkpoints\n\nsrc = '../input/d/pytorch/resnet50/resnet50.pth'\ndst = '/root/.cache/torch/hub/checkpoints/resnet50-0676ba61.pth'\n\nshutil.copy(src, dst)","metadata":{"execution":{"iopub.status.busy":"2023-04-18T02:52:42.390126Z","iopub.execute_input":"2023-04-18T02:52:42.391294Z","iopub.status.idle":"2023-04-18T02:52:45.467494Z","shell.execute_reply.started":"2023-04-18T02:52:42.391252Z","shell.execute_reply":"2023-04-18T02:52:45.466063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Import packages and define functions for image preprocessing and training","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torchvision.models as models\nimport torchvision.transforms as transforms\nimport requests\nimport os\nimport dicomsdl\nfrom tqdm import tqdm\nimport time\n# import pydicom\n# from pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom PIL import Image\nfrom skimage import exposure\nfrom matplotlib import pyplot as plt\n%matplotlib inline \n\nstart_time = time.time()\n\ndef voi_lut(dicom):\n    # Extract the pixel data from the DICOM image\n    pixel_data = dicom.pixelData()\n\n    # Extract the Window Center and Window Width from the DICOM image\n    window_center = dicom.WindowCenter if type(dicom.WindowCenter) != list else dicom.WindowCenter[0]\n    window_width = dicom.WindowWidth if type(dicom.WindowWidth) != list else dicom.WindowWidth[0]\n    \n#     print(window_center,window_width)\n    # NOT REALLY SURE why the window center and width sometimes are lists insteada of floats\n    # but for now we'll just take the first item from each list... but this doesnt work sadly\n\n    # Compute the minimum and maximum values for the VOI LUT\n    min_value = window_center - 0.5 * window_width\n    max_value = window_center + 0.5 * window_width\n\n    # Normalize the pixel data using the VOI LUT\n    normalized_pixel_data = np.clip(pixel_data, min_value, max_value)\n    normalized_pixel_data = (normalized_pixel_data - min_value) / (max_value - min_value)\n    \n    return normalized_pixel_data\n\n# Image normalization code found at https://www.kaggle.com/code/raddar/popular-x-ray-image-normalization-techniques/notebook\ndef read_xray(path, normalize = True, fix_monochrome = True):\n    \n#     dicom = pydicom.read_file(path)\n#     print(path)\n    dicom = dicomsdl.open(path)\n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n    if normalize: # this is causing problems with pydicom so we are going to ignore it for now https://github.com/pydicom/pylibjpeg/issues/58\n        data = voi_lut(dicom) # also this function only works with pydicom\n    else:\n        data = dicom.pixelData() \n    # depending on this value, X-ray may look inverted - so we reinvert the image here:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    data = data - np.min(data)\n    return data\n\ndef flatten_folder(folder):\n    files = []\n    total_files = 54713\n    with tqdm(total=total_files, desc='Loading up data') as pbar:\n        for root, dirs, filenames in os.walk(folder):\n            for file in filenames:\n                files.append(os.path.join(root, file))\n                pbar.update(1)\n    return files\n\n# Define the Focal Loss function\nclass FocalLoss(nn.Module):\n    def __init__(self, alpha, gamma):\n        super(FocalLoss, self).__init__()\n        self.alpha = alpha\n        self.gamma = gamma\n\n    def forward(self, inputs, targets):\n        BCE_loss = nn.BCELoss()(inputs, targets)\n        pt = torch.exp(-BCE_loss)\n        loss = self.alpha * (1-pt)**self.gamma * BCE_loss\n        return loss\n    \ndef pfbeta(labels, predictions, beta=1):\n    y_true_count = 0\n    ctp = 0\n    cfp = 0\n\n    for idx in range(len(labels)):\n        prediction = min(max(predictions[idx], 0), 1)\n        if (labels[idx]):\n            y_true_count += 1\n            ctp += prediction\n        else:\n            cfp += prediction\n\n    beta_squared = beta * beta\n    c_precision = ctp / (ctp + cfp)\n    if y_true_count == 0:\n        c_recall = 1\n    else:\n        c_recall = ctp / y_true_count\n    if (c_precision > 0 and c_recall > 0):\n        result = (1 + beta_squared) * (c_precision * c_recall) / (beta_squared * c_precision + c_recall)\n        return result\n    else:\n        return 0 \n\n    \n# Define the DICOM dataset class\nclass DICOMDataset(torch.utils.data.Dataset):\n    def __init__(self, data_dir, csv_file, transform=None):\n        self.dir_name = data_dir\n        self.flat_data_dir = flatten_folder(data_dir)\n\n        \n        \n        \n        \n        \n        proportion_of_data_to_use = 0.05; total_positive_examples = 1158 # Currently using only a percentage of the data for training.\n#         proportion_of_data_to_use = 0.0005; total_positive_examples = 10 # Currently using only a percentage of the data for training.\n    \n    \n    \n    \n    \n    \n    \n    \n        test_proportion = 0.2\n        proportion_of_data_to_test = (proportion_of_data_to_use/(1-test_proportion))*test_proportion # 0.05 - Test on 20% of the data used (train:test = 80:20)\n        proportion_of_positive_examples_in_test = test_proportion * (1158/(53548+1158)) # About 0.2*0.02 # Oversample the train positive label proportions but keep the test set label distribution the same to avoid overestimated results\n        self.data_dir = []\n        self.test_idxs = []\n        self.labels_df = pd.read_csv(csv_file)\n        positive_example_idxs = []\n        \n        num_samples_in_test = round(proportion_of_data_to_use*54706*test_proportion*(1158/(53548+1158)))\n        num_positives_in_train = round(total_positive_examples - num_samples_in_test)\n        \n        # HERE we get 1-(proportion_of_positive_examples_in_test)% of the positive examples so we make sure we use those for training - we want to train on all the positive examples, test using the remaining 10% of samples\n        for idx in range(len(self.flat_data_dir)):\n            dicom_file = os.path.join(self.dir_name, self.flat_data_dir[idx])\n            # Get the label for the image from the labels DataFrame\n            image_id = int(os.path.splitext(os.path.basename(dicom_file))[0])\n            label = self.labels_df.loc[self.labels_df['image_id'] == image_id, 'cancer'].values[0]\n            if label == 1:\n                if len(positive_example_idxs) <= num_positives_in_train - 1:\n                    self.data_dir.append(self.flat_data_dir[idx])\n                else:\n                    self.test_idxs.append(idx)\n                positive_example_idxs.append(idx)\n                \n        print(f'{len(self.test_idxs)} is the number of positive cases in the test set')\n        \n        # HERE we add in the enough of the negative examples to use up proportion_of_data_to_use of the dataset for train\n        idx_ = 0\n        while len(self.data_dir) < proportion_of_data_to_use*len(self.flat_data_dir):\n            if idx_ not in positive_example_idxs:\n                self.data_dir.append(self.flat_data_dir[idx])\n            idx_ = idx_ + 1\n            \n        print(f'{len(self.data_dir)} length of data train dir')\n            \n        # HERE we add in the enough of the negative examples to use up proportion_of_data_to_use of the dataset for test\n        while len(self.test_idxs) < proportion_of_data_to_test*len(self.flat_data_dir):\n            if idx_ not in positive_example_idxs:\n                self.test_idxs.append(idx)\n            idx_ = idx_ + 1\n        \n        print(f'{len(self.test_idxs)} length of test dir')\n        \n        self.transform = transform\n        self.image_data = []\n        self.dicom_filenames = []\n        for idx in tqdm(range(len(self.data_dir)), desc='Preprocessing data'):\n            dicom_file = os.path.join(self.dir_name, self.data_dir[idx])\n            self.dicom_filenames.append(dicom_file)\n            # print(self.data_dir[idx]) # <<< HERE is where we want to cross check with train.csv to EXCLUDE the views that are not CC and MLO\n            # the above line will print out /kaggle/input/rsna-breast-cancer-detection/train_images/10706/763186195.dcm for example\n            image = read_xray(os.path.join(self.dir_name, dicom_file))\n            # !!! Use histogram normalization of the image !!!\n            image = exposure.equalize_hist(image)  \n            image = Image.fromarray(image)\n            if self.transform:\n                image = self.transform(image)\n            # Convert grayscale image to RGB\n            image = np.stack([image] * 3, axis=1)\n            self.image_data.append(image)\n\n    def __len__(self):\n        return len(self.image_data)\n\n    def __getitem__(self, idx):\n        dicom_file = self.dicom_filenames[idx]\n        # Get the label for the image from the labels DataFrame\n        image_id = int(os.path.splitext(os.path.basename(dicom_file))[0])\n        label = self.labels_df.loc[self.labels_df['image_id'] == image_id, 'cancer'].values[0]\n        label = torch.tensor(label, dtype=torch.float32)\n\n        return self.image_data[idx], label\n    ","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2023-04-18T02:52:45.471688Z","iopub.execute_input":"2023-04-18T02:52:45.472034Z","iopub.status.idle":"2023-04-18T02:52:48.603332Z","shell.execute_reply.started":"2023-04-18T02:52:45.471999Z","shell.execute_reply":"2023-04-18T02:52:48.602278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# Load in the pretrained ResNet50 model\n# This code is deprecated so we should import ResNet50_Weights and do model = models.resnet50(weights=ResNet50_Weights.DEFAULT)\nmodel = models.resnet50(pretrained=True)\n\n# Freeze all layers except the last one\nfor param in model.parameters():\n    param.requires_grad = False\nmodel.fc.requires_grad = True\n\n# Replace the last fully connected layer with a few new layer with dropout and 1 output node\nmodel.fc = nn.Sequential(\n    nn.Linear(2048, 512),\n    nn.ReLU(),\n    nn.Dropout(0.5),\n    nn.Linear(512, 128),\n    nn.ReLU(),\n    nn.Dropout(0.5),\n    nn.Linear(128, 1),\n    nn.Sigmoid()\n)\n\n# Load the DICOM data\n# FOR NOW, we only load in the images separately - consider combining pairs of images or even quadruplets, which will require changing the output layer\ndata_dir = '/kaggle/input/rsna-breast-cancer-detection/train_images'\ncsv_file = '/kaggle/input/rsna-breast-cancer-detection/train.csv'\ndataset = DICOMDataset(data_dir, csv_file, transform=transforms.Compose([\n    transforms.Resize(256),\n    transforms.CenterCrop(224), # need to check image to see if anything important is cropped out, if so add borders\n    transforms.ToTensor()\n#     ,transforms.Normalize([0.485, 0.456, 0.406],...\n    ]))\n\ntest_data_indexes = dataset.test_idxs\n\n# Define the device to use for training\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\n# Split the data into training and validation sets\ntrain_data, val_data = torch.utils.data.random_split(dataset, [int(0.875*len(dataset)), len(dataset)-int(0.875*len(dataset))])\n\n# Create the data loaders\nbatch_size = 32\ntrain_loader = torch.utils.data.DataLoader(train_data, batch_size=batch_size, shuffle=True)\nval_loader = torch.utils.data.DataLoader(val_data, batch_size=batch_size, shuffle=False)\n\n# Set loss function and optimizer\nnum_positive_samples=1158; num_negative_samples=53548 # Set the alpha of the focal loss function to be the ratio of positive to negative samples\ncriterion = FocalLoss(alpha=num_positive_samples/num_negative_samples, gamma=0.5) # set gamma to 0.5 for now, if we go higher then we may over prioritize recall over precision\noptimizer = optim.SGD(model.parameters(), lr=0.001, momentum=0.9) # defaults, can modify\n\nprint('Starting up model')\n# Move the model to the device\nmodel.to(device)\n\n# Specify the number of epochs to train for and the early stopping patience\n\n\n\n\n\n\n\n\n\n\n\n# num_epochs = 2 # At 200 epochs it looks like we should keep going\nnum_epochs = 200\n\n\n\n\n\n\n\n\n\npatience = num_epochs/10\n\n# Keep track of the best validation loss and the number of epochs without improvement\nbest_loss = float(\"inf\")\ncounter = 0\n\nimport pickle\n\n# # Saving train_loader\n# with open('train_loader.pkl', 'wb') as f:\n#     pickle.dump(train_loader, f)\n\n# # Saving val_loader\n# with open('val_loader.pkl', 'wb') as f:\n#     pickle.dump(val_loader, f)\n\n# # Loading train_loader\n# with open('train_loader.pkl', 'rb') as f:\n#     train_loader = pickle.load(f)\n\n# # Loading val_loader\n# with open('val_loader.pkl', 'rb') as f:\n#     val_loader = pickle.load(f)\n\nend_data_time = time.time()\ntotal_time = end_data_time - start_time\n\nprint(\"The data preprocessing took {:.0f}h {:.0f}m {:.0f}s\".format(total_time / 3600, (total_time % 3600) / 60, total_time % 60))","metadata":{"execution":{"iopub.status.busy":"2023-04-18T02:52:48.605769Z","iopub.execute_input":"2023-04-18T02:52:48.606777Z","iopub.status.idle":"2023-04-18T03:37:18.97036Z","shell.execute_reply.started":"2023-04-18T02:52:48.606732Z","shell.execute_reply":"2023-04-18T03:37:18.969329Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''In this example, the DICOM image is loaded using the pydicom.read_file function. \nThe pixel data is extracted from the DICOM image using the pixel_array property. \nThe Window Center and Window Width values are extracted from the DICOM image using the WindowCenter and WindowWidth properties. \nThe minimum and maximum values for the VOI LUT are computed using the Window Center and Window Width values. \nThe pixel data is then normalized using the VOI LUT by clipping the pixel values to the minimum and maximum values, \nand then scaling the pixel values to the range [0, 1]. Finally, the normalized image is displayed using matplotlib.pyplot.\n'''\n\n# # Display the normalized image using matplotlib\n# import matplotlib.pyplot as plt\n# plt.imshow(normalized_pixel_data, cmap='bone')\n# plt.show()\n\n\n# Code for undersampling negatives/oversampling positives\n\n# import pandas as pd\n# import numpy as np\n\n# # Load the data into a pandas DataFrame\n# df = pd.read_csv(\"your_dataset.csv\")\n\n# # Split the data into positive and negative examples\n# pos_examples = df[df[\"label\"] == 1]\n# neg_examples = df[df[\"label\"] == 0]\n\n# # Determine the size of the smaller and larger datasets\n# min_size = min(pos_examples.shape[0], neg_examples.shape[0])\n# max_size = max(pos_examples.shape[0], neg_examples.shape[0])\n\n# # Oversample the smaller dataset and undersample the larger dataset to create a balanced dataset\n# if pos_examples.shape[0] < neg_examples.shape[0]:\n#     pos_examples_sampled = pos_examples.sample(min_size, replace=True)\n#     neg_examples_sampled = neg_examples.sample(min_size, replace=False)\n# else:\n#     pos_examples_sampled = pos_examples.sample(max_size, replace=False)\n#     neg_examples_sampled = neg_examples.sample(max_size, replace=True)\n\n# # Combine the oversampled/undersampled datasets into a single balanced dataset\n# balanced_df = pd.concat([pos_examples_sampled, neg_examples_sampled])\n\n# # Shuffle the rows of the balanced dataset\n# balanced_df = balanced_df.sample(frac=1, random_state=42).reset_index(drop=True)\n\n# # Save the balanced dataset to a file\n# balanced_df.to_csv(\"balanced_dataset.csv\", index=False)\n\n# Better code for implementing the split into train:validate:test datasets\n\n# import torch\n# from torch.utils.data import Dataset, DataLoader\n# from torch.utils.data.sampler import Sampler\n\n# class StratifiedSampler(Sampler):\n#     def __init__(self, labels, num_samples):\n#         self.labels = labels\n#         self.num_samples = num_samples\n\n#     def __iter__(self):\n#         labels = torch.from_numpy(self.labels)\n#         positive_indices = (labels == 1).nonzero().squeeze(1)\n#         negative_indices = (labels == 0).nonzero().squeeze(1)\n        \n#         positive_sampler = torch.utils.data.sampler.SubsetRandomSampler(positive_indices, replacement=False, num_samples=self.num_samples // 2)\n#         negative_sampler = torch.utils.data.sampler.SubsetRandomSampler(negative_indices, replacement=False, num_samples=self.num_samples // 2)\n        \n#         return iter(torch.cat([i for i in positive_sampler], [i for i in negative_sampler]))\n\n#     def __len__(self):\n#         return self.num_samples\n\n# class CustomDataset(Dataset):\n#     def __init__(self, data, labels):\n#         self.data = data\n#         self.labels = labels\n\n#     def __len__(self):\n#         return len(self.data)\n\n#     def __getitem__(self, idx):\n#         return self.data[idx], self.labels[idx]\n\n# # Create the dataset\n# data = ... # Your data\n# labels = ... # Your labels\n# dataset = CustomDataset(data, labels)\n\n# # Split the dataset into training, validation, and test sets\n# num_samples = len(dataset)\n# train_size = int(0.8 * num_samples)\n# val_size = int(0.1 * num_samples)\n# test_size = num_samples - train_size - val_size\n\n# train_sampler = StratifiedSampler(labels, train_size)\n# val_sampler = StratifiedSampler(labels, val_size)\n# test_sampler = StratifiedSampler(labels, test_size)\n\n# train_loader = DataLoader(dataset, batch_size=batch_size, sampler=train_sampler)\n# val_loader = DataLoader(dataset, batch_size=batch_size, sampler=val_sampler)\n# test_loader = DataLoader(dataset, batch_size=batch_size, sampler=test_sampler)\n\n# Automated machine learning (AutoML) provides a way to automate the process of choosing the best data augmentation techniques for a specific dataset and model. The general process for using AutoML for data augmentation can be broken down into the following steps:\n\n# Prepare the data: Load the dataset into memory and perform any necessary preprocessing steps, such as normalizing the values.\n\n# Split the data: Split the data into a training set, validation set, and test set, with a stratified sampling approach to ensure a balanced distribution of positive and negative examples.\n\n# Define the model architecture: Define the model architecture, such as a Convolutional Neural Network (CNN) or Recurrent Neural Network (RNN), that you want to use for the task.\n\n# Define the search space: Define the search space for the AutoML method, which includes the types of data augmentation techniques that you want to consider, such as flipping, rotation, and scaling, as well as the range of values for each technique.\n\n# Train the model: Train the model using an AutoML framework, such as H2O.ai or TPOT, and let it search for the best data augmentation techniques for your specific dataset and model architecture.\n\n# Evaluate the results: Evaluate the performance of the model on the validation set, and select the best performing model based on the accuracy, F1 score, or other performance metric of interest.\n\n# Test the model: Test the model on the test set and use it to make predictions on new, unseen data.\n\n# It is important to note that while AutoML provides a convenient way to search for the optimal data augmentation techniques, it can be computationally expensive, especially for large datasets and complex model architectures. As such, it is important to carefully consider the trade-off between accuracy and computational cost when using AutoML for data augmentation.","metadata":{"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-04-18T03:37:18.972175Z","iopub.execute_input":"2023-04-18T03:37:18.972854Z","iopub.status.idle":"2023-04-18T03:37:18.987369Z","shell.execute_reply.started":"2023-04-18T03:37:18.972808Z","shell.execute_reply":"2023-04-18T03:37:18.986226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nprint('Training the model')\n\n# Lists to store the values of the metrics over the epochs\nval_pfbeta_list = []\nepoch_loss_list = []\nepoch_acc_list = []\n\n# Train the model for each epoch\nfor epoch in tqdm(range(num_epochs)):\n    running_loss = 0.0\n    running_corrects = 0\n    \n    # Early stopping code\n    # Evaluate the model on the validation data\n    with torch.no_grad():\n        val_loss = 0\n        val_pfbeta = 0\n        for inputs, labels in val_loader:\n            inputs = inputs.squeeze()\n            inputs = inputs.to(device)\n            labels = labels.to(device)\n            outputs = model(inputs).squeeze()\n            loss = criterion(outputs, labels)\n            val_loss += loss.item()\n            val_pfbeta += pfbeta(labels, outputs)\n        val_loss /= len(val_loader)\n        val_pfbeta /= len(val_loader)\n        val_pfbeta_list.append(val_pfbeta)\n        print('Validation Probabilistic F1-score: ',val_pfbeta)\n\n    # If the validation loss has not decreased for `patience` epochs, stop training\n    if val_loss > best_loss:\n        counter += 1\n        if counter >= patience:\n            break\n    else:\n        # Update the best validation loss and reset the counter\n        best_loss = val_loss\n        counter = 0\n\n    # Split the data into small batches\n    for inputs, labels in train_loader:\n        inputs=inputs.squeeze()\n        inputs = inputs.to(device)\n        labels = labels.to(device)\n\n        # Zero the gradients\n        optimizer.zero_grad()\n\n        # Forward pass\n        outputs = model(inputs)\n        _, preds = torch.max(outputs, 1)\n        loss = criterion(outputs.squeeze(), labels)\n\n        # Backward pass and optimization\n        loss.backward()\n        optimizer.step()\n\n        # Update running loss and accuracy\n        running_loss += loss.item() * inputs.size(0)\n        running_corrects += torch.sum(preds == labels.data)\n\n    # Calculate average loss and accuracy for the epoch\n    epoch_loss = running_loss / len(train_loader.dataset)\n    epoch_acc = running_corrects.double() / len(train_loader.dataset)\n    epoch_loss_list.append(epoch_loss)\n    epoch_acc_list.append(epoch_acc)\n\n    print(f'Epoch [{epoch + 1}/{num_epochs}], Loss: {epoch_loss:.4f}, Accuracy: {epoch_acc:.4f}')\n\n    # Save the model state after every 10% of total epochs\n    if epoch%(num_epochs/10) == (num_epochs/10)-1:\n        torch.save(model.state_dict(), f'resnet50_epoch_{epoch + 1}.pth')\n    \n    \n    \n    \n    ''' Try with unfreezing every 10 epochs instead of 2 '''\n    # Unfreeze one more layer after every 10 epochs\n    if (epoch+1) % 10 == 0:\n        for name, child in model.named_children():\n            if name == 'fc':\n                break\n            for param in child.parameters():\n                param.requires_grad = True\n\n                \n                \n                \n# Save the trained model\ntorch.save(model.state_dict(),'/kaggle/working/resnet50_breast_cancer.pth')\n\nend_training_time = time.time()\ntotal_time = end_training_time - end_data_time\n\nprint(\"The training took {:.0f}h {:.0f}m {:.0f}s\".format(total_time / 3600, (total_time % 3600) / 60, total_time % 60))","metadata":{"execution":{"iopub.status.busy":"2023-04-18T03:37:18.989146Z","iopub.execute_input":"2023-04-18T03:37:18.989713Z","iopub.status.idle":"2023-04-18T04:22:53.937224Z","shell.execute_reply.started":"2023-04-18T03:37:18.989669Z","shell.execute_reply":"2023-04-18T04:22:53.935953Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(6, 10))\n\n# Plot the training loss against the epochs\nax1 = plt.subplot(2, 1, 1)\nplt.plot(list(range(1, len(epoch_loss_list)+1)), epoch_loss_list)\nplt.title(\"Training Loss vs Epochs\")\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Training Loss\")\nplt.xlim(left=0)\nplt.ylim(bottom=0)\n\n# # Plot the accuracy against the epochs\n# if type(epoch_acc_list[0]) == torch.Tensor:\n#     epoch_acc_list=np.stack([t.cpu().numpy() for t in epoch_acc_list])\n# ax2 = plt.subplot(3, 1, 2)\n# plt.plot(list(range(1, len(epoch_acc_list)+1)), epoch_acc_list)\n# plt.title(\"Accuracy vs Epochs\")\n# plt.xlabel(\"Epochs\")\n# plt.ylabel(\"Accuracy\")\n\n# Plot the probabilistic F1-score against the epochs\nif type(val_pfbeta_list[0]) == torch.Tensor:\n    val_pfbeta_list=np.stack([t.cpu().numpy() for t in val_pfbeta_list])\nax2 = plt.subplot(2, 1, 2)\nplt.plot(list(range(1, len(val_pfbeta_list)+1)), val_pfbeta_list)\nplt.title(\"Probabilistic F1-Score vs Epochs\")\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Probabilistic F1-Score\")\nplt.xlim(left=0)\nplt.ylim(bottom=0)\n\n# Show the plot\nplt.tight_layout()\nplt.show()\n\nplt.savefig('/kaggle/working/training_loss_and_pf1_vs_epochs.png')","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:22:53.939054Z","iopub.execute_input":"2023-04-18T04:22:53.940624Z","iopub.status.idle":"2023-04-18T04:22:54.428018Z","shell.execute_reply.started":"2023-04-18T04:22:53.940582Z","shell.execute_reply":"2023-04-18T04:22:54.426968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import roc_curve, auc\nimport matplotlib.pyplot as plt\n\n\n\n\n\n\n\n\n# test_data_indexes = [0,1,2,3] # comment in for quick testing purposes only\n\n\n\n\n\n\n\n\n\n\n\n\n# Test on the held out test set - note that if you train on the entire dataset then you'll just be testing on training data \n# so this is only for when you are training on less than 90% of the data.\n\n# Define the DICOM dataset class\nclass DICOMTestingDataset(torch.utils.data.Dataset):\n    def __init__(self, data_dir, csv_file, test_data_indexes, transform=None):\n        self.dir_name = data_dir\n        self.flat_data_dir = flatten_folder(data_dir)\n        self.data_dir = [self.flat_data_dir[i] for i in test_data_indexes]\n        self.labels_df = pd.read_csv(csv_file)\n        self.transform = transform\n        self.image_data = []\n        self.dicom_filenames = []\n        progress_bar = tqdm(total=len(self.data_dir), desc='Preprocessing data')\n        for idx in range(len(self.data_dir)):\n            dicom_file = os.path.join(self.dir_name, self.data_dir[idx])\n            self.dicom_filenames.append(dicom_file)\n            # print(self.data_dir[idx]) # <<< HERE is where we want to cross check with train.csv to EXCLUDE the views that are not CC and MLO... but\n            # if we do that and we get some data in our test set that isn't either of those views then what? We can still just ignore them...\n            # the above line will print out /kaggle/input/rsna-breast-cancer-detection/train_images/10706/763186195.dcm for example\n            image = read_xray(os.path.join(self.dir_name, dicom_file))\n            # !!! Use histogram normalization of the image !!!\n            image = exposure.equalize_hist(image)  \n            image = Image.fromarray(image)\n            if self.transform:\n                image = self.transform(image)\n            # Convert grayscale image to RGB\n            image = np.stack([image] * 3, axis=1)\n            self.image_data.append(image)\n            progress_bar.update(1)\n        progress_bar.close()\n\n    def __len__(self):\n        return len(self.image_data)\n\n    def __getitem__(self, idx):\n        dicom_file = self.dicom_filenames[idx]\n        # Get the label for the image from the labels DataFrame\n        image_id = int(os.path.splitext(os.path.basename(dicom_file))[0])\n        label = self.labels_df.loc[self.labels_df['image_id'] == image_id, 'cancer'].values[0]\n        label = torch.tensor(label, dtype=torch.float32)\n\n        return self.image_data[idx], label\n    \ndf_train = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/train.csv')\n\n# Read the test data\ndata_dir = '/kaggle/input/rsna-breast-cancer-detection/train_images'\ncsv_file='/kaggle/input/rsna-breast-cancer-detection/train.csv'\ntest_data_indexes = test_data_indexes\n\n# Preprocess the test data\ntest_dataset = DICOMTestingDataset(data_dir, csv_file, test_data_indexes, transform=transforms.Compose([\n    transforms.Resize(256),\n    transforms.CenterCrop(224), # need to check image to see if anything important is cropped out, if so add borders\n    transforms.ToTensor()\n#     ,transforms.Normalize([0.485, 0.456, 0.406],...\n    ]))\n\ntest_loader = torch.utils.data.DataLoader(test_dataset, shuffle=False)\n\n# Use the model to make predictions\nmodel.eval()\n\nsensitivity, specificity = 0,0\n\n# Initialize true_labels and all_predictions lists\ntrue_labels = []\nall_predictions = []\n\nwith torch.no_grad():\n    predictions = []\n    true_positive = 0; false_positive = 0; false_negative = 0; true_negative = 0\n    \n    test_pfbeta = 0\n    for inputs, labels in tqdm(test_loader):\n        inputs, labels = inputs.to(device), labels.to(device)\n        output = model(inputs.squeeze(1)).squeeze()\n        \n        # Append true labels and predictions to the lists\n        true_labels.append(labels.cpu().numpy())\n        all_predictions.append(output.cpu().numpy())\n        \n        predictions.append(output.item())\n        if labels:\n            true_positive += output > 0.5\n            false_negative += output <= 0.5\n        else:\n            false_positive += output > 0.5\n            true_negative += output <= 0.5\n        outputs = model(inputs.squeeze(1))\n        test_pfbeta += pfbeta(labels, outputs)\n\n    if (true_positive + false_negative) == 0:\n        sensitivity = 1\n    else:\n        sensitivity = true_positive / (true_positive + false_negative)\n    if (true_negative + false_positive) == 0:\n        specificity = 1\n    else:\n        specificity = true_negative / (true_negative + false_positive)\n    test_pfbeta /= len(test_loader)\n\n    \nproportion_of_data_to_use = 0.0625\nif type(sensitivity) != float and type(sensitivity) != int:\n    sensitivity = sensitivity.item()\nif type(specificity) != float and type(specificity) != int:\n    specificity = specificity.item()\nif type(test_pfbeta) != float and type(test_pfbeta) != int:\n    test_pfbeta = test_pfbeta.item()\nprint(f'For the held-out test set using {proportion_of_data_to_use*100}% of the data, the sensitivity is {sensitivity*100:.1f}%, the specificity is {specificity*100:.1f}%, and the probabilistic F1-Score is {test_pfbeta:.2f}.')\n\n# Convert true_labels and all_predictions lists to NumPy arrays\ntrue_labels = np.array(true_labels)\nall_predictions = np.array(all_predictions)\n\n# Calculate true positive rates and false positive rates\nfpr, tpr, thresholds = roc_curve(true_labels, all_predictions)\n\n# Calculate the ROC-AUC score\nroc_auc = auc(fpr, tpr)\n\n# Plot the ROC curve\nplt.figure()\nlw = 2\nplt.plot(fpr, tpr, color='darkorange', lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)\nplt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\nplt.xlim([0.0, 1.0])\nplt.ylim([0.0, 1.05])\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.title('Receiver Operating Characteristic')\nplt.legend(loc=\"lower right\")\n\n# Save the ROC curve plot\nplt.savefig('/kaggle/working/roc_curve.png')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:22:54.430208Z","iopub.execute_input":"2023-04-18T04:22:54.430665Z","iopub.status.idle":"2023-04-18T04:30:18.084545Z","shell.execute_reply.started":"2023-04-18T04:22:54.430622Z","shell.execute_reply":"2023-04-18T04:30:18.083391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# Plot the ROC curve\nplt.figure(figsize=(8, 6))  # Increase the figure size\nlw = 2\nplt.plot(fpr, tpr, color='darkorange', lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)\nplt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\nplt.xlim([0.0, 1.0])\nplt.ylim([0.0, 1.05])\n\n# Set the font size and labels for the axes\nplt.xlabel('False Positive Rate', fontsize=14)\nplt.ylabel('True Positive Rate', fontsize=14)\n\n# Set the title and its font size\nplt.title('Receiver Operating Characteristic', fontsize=16)\n\n# Set the font size for the legend and place it in the lower right corner\nlegend = plt.legend(loc=\"lower right\", fontsize=12)\nplt.setp(legend.get_lines(), linewidth=2)  # Set the linewidth of the legend lines\n\n# Set the font size for the tick labels\nplt.xticks(fontsize=12)\nplt.yticks(fontsize=12)\n\n# Save the plot as a high-resolution image\nplt.savefig(\"ROC_curve_publication_ready.png\", dpi=300)\n\n# Display the plot\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:18.089272Z","iopub.execute_input":"2023-04-18T04:30:18.089615Z","iopub.status.idle":"2023-04-18T04:30:18.829402Z","shell.execute_reply.started":"2023-04-18T04:30:18.089565Z","shell.execute_reply":"2023-04-18T04:30:18.828084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\nimport seaborn as sns\n# Calculate confusion matrix\nthreshold = 0.5  # You can adjust this value based on your requirements\nall_predictions_binary = (all_predictions > threshold).astype(int)\ncm = confusion_matrix(true_labels, all_predictions_binary)\n\n# Plot the confusion matrix\nfig, ax = plt.subplots(figsize=(6, 6))\nsns.heatmap(cm, annot=True, fmt=\"d\", cmap=\"YlGnBu\", cbar=False, square=True,\n            xticklabels=[\"No cancer\", \"Cancer\"], yticklabels=[\"No cancer\", \"Cancer\"], ax=ax)\n\n# Customize the appearance of the plot\nax.set_xlabel(\"Predicted label\", fontsize=14, labelpad=10)\nax.set_ylabel(\"True label\", fontsize=14, labelpad=10)\nax.set_title(\"Confusion Matrix\", fontsize=18, pad=20)\nax.tick_params(axis='both', which='major', labelsize=12)\n\n# Annotate the cells with percentages\nnum_cells = cm.shape[0] * cm.shape[1]\nfor p, (i, j) in zip(ax.texts, [(i, j) for i in range(cm.shape[0]) for j in range(cm.shape[1])]):\n    percentage = cm[i, j] / cm.sum() * 100\n    p.set_text(f\"{cm[i, j]}\\n({percentage:.1f}%)\")\n\n# Save the confusion matrix plot\nplt.tight_layout()\nplt.savefig('/kaggle/working/confusion_matrix.png')\nplt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:18.830952Z","iopub.execute_input":"2023-04-18T04:30:18.831325Z","iopub.status.idle":"2023-04-18T04:30:19.32769Z","shell.execute_reply.started":"2023-04-18T04:30:18.831287Z","shell.execute_reply":"2023-04-18T04:30:19.325963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Test on the hidden test set for submission\n\n# Define the DICOM dataset class\nclass DICOMTestDataset(torch.utils.data.Dataset):\n    def __init__(self, data_dir, transform=None):\n        self.dir_name = data_dir\n        self.data_dir = flatten_folder(data_dir)\n        proportion_of_data_to_use = 1 # Currently using only a percentage of the data for training.\n        self.data_dir = self.data_dir[:int(proportion_of_data_to_use*len(self.data_dir))]\n        self.transform = transform\n        self.image_data = []\n        self.dicom_filenames = []\n        progress_bar = tqdm(total=len(self.data_dir), desc='Preprocessing data')\n        for idx in range(len(self.data_dir)):\n            dicom_file = os.path.join(self.dir_name, self.data_dir[idx])\n            self.dicom_filenames.append(dicom_file)\n            # print(self.data_dir[idx]) # <<< HERE is where we want to cross check with train.csv to EXCLUDE the views that are not CC and MLO... but\n            # if we do that and we get some data in our test set that isn't either of those views then what? We can still just ignore them...\n            # the above line will print out /kaggle/input/rsna-breast-cancer-detection/train_images/10706/763186195.dcm for example\n            image = read_xray(os.path.join(self.dir_name, dicom_file))\n            # !!! Use histogram normalization of the image !!!\n            image = exposure.equalize_hist(image)  \n            image = Image.fromarray(image)\n            if self.transform:\n                image = self.transform(image)\n            # Convert grayscale image to RGB\n            image = np.stack([image] * 3, axis=1)\n            self.image_data.append(image)\n            progress_bar.update(1)\n        progress_bar.close()\n\n    def __len__(self):\n        return len(self.image_data)\n\n    def __getitem__(self, idx):\n        dicom_file = self.dicom_filenames[idx]\n\n        return self.image_data[idx]\n    \ndf_test = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/test.csv')\n\n# Read the test data\ndata_dir = '/kaggle/input/rsna-breast-cancer-detection/test_images'\ncsv_file='/kaggle/input/rsna-breast-cancer-detection/test.csv'\n\n# Preprocess the test data\ntest_dataset = DICOMTestDataset(data_dir, transform=transforms.Compose([\n    transforms.Resize(256),\n    transforms.CenterCrop(224), # need to check image to see if anything important is cropped out, if so add borders\n    transforms.ToTensor()\n#     ,transforms.Normalize([0.485, 0.456, 0.406],...\n    ]))\n\ntest_loader = torch.utils.data.DataLoader(test_dataset, shuffle=False)\n\n# Use the model to make predictions\nmodel.eval()\n\nwith torch.no_grad():\n    cancer = []\n    for inputs in tqdm(test_loader):\n        inputs = inputs.to(device)\n        output = model(inputs.squeeze(1)).squeeze()\n#         output = torch.sigmoid(output)\n        cancer.append(output.item())\n\n\nsubmission = pd.DataFrame({'prediction_id': df_test['prediction_id'], 'cancer': cancer})\nsubmission = submission.groupby('prediction_id').mean().reset_index()\n# submission = submission.groupby('prediction_id').max().reset_index()\n\nsubmission.to_csv('submission.csv', index=False)\n\nend_testing_time = time.time()\ntotal_time = end_testing_time - end_training_time\n\nprint(\"The testing took {:.0f}h {:.0f}m {:.0f}s\".format(total_time / 3600, (total_time % 3600) / 60, total_time % 60))\n\nsubmission","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.335262Z","iopub.execute_input":"2023-04-18T04:30:19.339769Z","iopub.status.idle":"2023-04-18T04:30:19.356964Z","shell.execute_reply.started":"2023-04-18T04:30:19.339699Z","shell.execute_reply":"2023-04-18T04:30:19.355648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To make a test submission just comment out every other cell and submit this\n\n# import pandas as pd\n\n# # Test on the hidden test set for submission (here we just predict 0 for every point)\n\n# df_test = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/test.csv')\n\n# cancer = []\n# for item in df_test['prediction_id']:\n#     cancer.append(0)\n\n# submission = pd.DataFrame({'prediction_id': df_test['prediction_id'], 'cancer': cancer})\n# submission = submission.groupby('prediction_id').mean().reset_index()\n# # submission = submission.groupby('prediction_id').max().reset_index()\n\n# submission.to_csv('submission.csv', index=False)\n\n# import os\n# os.chdir(r'../working')\n# from IPython.display import FileLink\n# FileLink(r'submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.358799Z","iopub.execute_input":"2023-04-18T04:30:19.359247Z","iopub.status.idle":"2023-04-18T04:30:19.372777Z","shell.execute_reply.started":"2023-04-18T04:30:19.359201Z","shell.execute_reply":"2023-04-18T04:30:19.371839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import pickle\n# # Saving train_loader\n# with open('train_loader.pkl', 'wb') as f:\n#     pickle.dump(train_loader, f)\n\n# # Saving val_loader\n# with open('val_loader.pkl', 'wb') as f:\n#     pickle.dump(val_loader, f)","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.374358Z","iopub.execute_input":"2023-04-18T04:30:19.375057Z","iopub.status.idle":"2023-04-18T04:30:19.389005Z","shell.execute_reply.started":"2023-04-18T04:30:19.375017Z","shell.execute_reply":"2023-04-18T04:30:19.388269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Saving train_loader\n# with open('train_data.pkl', 'wb') as f:\n#     pickle.dump(train_data, f)\n\n# # Saving val_loader\n# with open('val_data.pkl', 'wb') as f:\n#     pickle.dump(val_data, f)","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.390378Z","iopub.execute_input":"2023-04-18T04:30:19.391031Z","iopub.status.idle":"2023-04-18T04:30:19.400933Z","shell.execute_reply.started":"2023-04-18T04:30:19.390995Z","shell.execute_reply":"2023-04-18T04:30:19.399933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nos.chdir(r'../working')\nfrom IPython.display import FileLink\nFileLink(r'resnet50_breast_cancer.pth')","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.402328Z","iopub.execute_input":"2023-04-18T04:30:19.402674Z","iopub.status.idle":"2023-04-18T04:30:19.415043Z","shell.execute_reply.started":"2023-04-18T04:30:19.402643Z","shell.execute_reply":"2023-04-18T04:30:19.414146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_pfbeta = pfbeta(true_labels, all_predictions)\nprint(f'For the held-out test set the Probabilistic F1 Score is {test_pfbeta}')","metadata":{"execution":{"iopub.status.busy":"2023-04-18T04:30:19.41663Z","iopub.execute_input":"2023-04-18T04:30:19.417402Z","iopub.status.idle":"2023-04-18T04:30:19.431473Z","shell.execute_reply.started":"2023-04-18T04:30:19.417364Z","shell.execute_reply":"2023-04-18T04:30:19.430342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}