{"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\n\nAs an improvement upon the weighted mean baseline ([notebook](https://www.kaggle.com/code/pankajpansari/baseline-0-weighted-mean-probabilities)), we're going to use a 2D CNN to make predictions.\n\nA CT scan series is composed of many slices, that is, 2D images. Hence, each series is a sequence of 2D images. Let us ignore this sequential information, and think about making predictions on single 2D images. Two questions arise:\n\n1. How do we make the training set?\n    The labels are assigned to patients. Each patient can have multiple CT scans. Each CT scan series has many images. Let's do the simplest thing and assign to each image in each series, the label of the respective patient for all injury types.\n    \n2. Given predictions at image level, how do we combine them into predictions for patients?\n    The first idea is to use the mean. However, it may be that only a few images in the CT scan indicate some type of injury and that may be enough for the injury to be present. We could use max to aggregrate, but to avoid the noise/inaccuracies in the predictions, let's use 75% quantile value.\n\nSo in this notebook, we are going to fine-tune a pretrained CNN model using the CT scan images from the competition. Specifically, we are going to take a ResNet-18 model pre-trained on ImageNet and fine tune it on our CT scan data. We'll make use of the [fastai](https://docs.fast.ai) library.\n\nHere, we are not going to make use of segmentation data, DICOM tags, or meta data. To make our task simpler, we'll work with 128x128 images. We obtained this image data by taking the PNG dataset of [notebook](https://www.kaggle.com/code/theoviel/get-started-quicker-dicom-png-conversion) and resizing them beforehand using [this script](https://www.kaggle.com/code/pankajpansari/rsna-2023-abdominal-trauma-converting-dicom-to-png).","metadata":{}},{"cell_type":"markdown","source":"## Code\n\n### EDA\nWe import the necessary modules.","metadata":{}},{"cell_type":"code","source":"!pip install kornia pydicom","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os, random\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nfrom fastai.vision.all import *\nfrom fastai.basics import *\nfrom fastai.callback.all import *\nfrom fastai.medical.imaging import *\nimport shutil\nimport pydicom\nimport cv2\nimport glob\nimport time\nimport seaborn as sns\nimport torchvision\nfrom sklearn.model_selection import train_test_split\nfrom joblib import Parallel, delayed\n\nfrom tqdm.notebook import tqdm\nfrom joblib import Parallel, delayed\n\nrandom.seed(42)","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BASE_DIR = '.'","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us look at the labels provided in *train.csv*. We note that the labels are assigned to each patient and not to each CT scan image. There are many CT images associated with each patient; some of them may reflect the injury and some may not, depending upon their z-coordinate.\n\nHowever, for simplification, for each CT scan image we'll simply assign the label given to the respective patient.","metadata":{}},{"cell_type":"code","source":"patient_labels = pd.read_csv(os.path.join(BASE_DIR, 'train.csv'))\npatient_labels.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We do a sanity check that the *healthy* and *injury* (or *low injury* and *high injury*) probabilities add up to 1 for each type of injury.","metadata":{}},{"cell_type":"code","source":"#Check data has no NANs\nprint(patient_labels.isnull().any().any())\n#Check healthy and injury labels are complementary\nprint((patient_labels['bowel_healthy'] == np.abs(1 - patient_labels['bowel_injury'])).all())\nprint((patient_labels['extravasation_healthy'] == np.abs(1 - patient_labels['extravasation_injury'])).all())\nprint((patient_labels['kidney_healthy'] == np.abs(1 - patient_labels['kidney_low'] - patient_labels['kidney_high'])).all())\nprint((patient_labels['liver_healthy'] == np.abs(1 - patient_labels['liver_low'] - patient_labels['liver_high'])).all())\nprint((patient_labels['spleen_healthy'] == np.abs(1 - patient_labels['spleen_low'] - patient_labels['spleen_high'])).all())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The *any_injury* label seems a bit redundant since it can be derived from the other injury labels. It simply means whether at least one type of bowel/extravsation/liver/spleen/kidney injury is present.","metadata":{}},{"cell_type":"code","source":"#Check consistency of any_injury label with respect to other labels\nprint((patient_labels['any_injury'] == 1 - np.min(patient_labels[['bowel_healthy', 'extravasation_healthy', \n                                                          'kidney_healthy', 'liver_healthy', 'spleen_healthy']], axis = 1)).all())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Checking Data Imabalance and Correlations\n\nLet us find the percentage of samples having different forms of injuries. This will help us know the distribution of the labels and find out whether there is imbalance in positive/negative samples.","metadata":{}},{"cell_type":"code","source":"def get_pos_percent(df, label):\n    num_entries = df.shape[0]\n    return (df[label] == 1).sum()*100/num_entries\n\ninjury_labels = ['bowel_injury', 'extravasation_injury' , 'kidney_low' , 'kidney_high', 'liver_low',\n                 'liver_high', 'spleen_low', 'spleen_high']\n\n\nfor label in injury_labels:    \n    print(f'% of {label} samples = {get_pos_percent(patient_labels, label): .3f}')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Indeed there is a lot of imbalance between number of positive and negative samples for each type of injury.\n\nNow, let us check how much the injuries are correlated among themselves. We'll plot a heatmap of the Pearson correlation coefficients (between -1 and 1).","metadata":{}},{"cell_type":"code","source":"sns.heatmap(patient_labels[injury_labels].corr(), vmin = -1, vmax = 1, annot = True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The correlations are very weak, with most values close to 0 and maximum as 0.2.","metadata":{}},{"cell_type":"markdown","source":"### Data Preparation\n\n#### Splitting dataset\n\nLet us first form training and validation sets. Because of the imbalance, we would like to ensure that the class/injury distribution is similar in the three sets.","metadata":{}},{"cell_type":"code","source":"train_patients, validation_patients = train_test_split(patient_labels, test_size = 0.2, random_state = 42) ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Number of samples')\nprint(f'Training set: {train_patients.shape[0]}   Validation set: {validation_patients.shape[0]}')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'% of positive samples - Training, Validation')\nfor label in injury_labels:    \n    print(f'{label}: {get_pos_percent(train_patients, label): .3f}, {get_pos_percent(validation_patients, label): .3f}')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The distribution of the labels in not exactly equal in the training and validation sets, but it is more or less similar. For now, this will do for our purposes.","metadata":{}},{"cell_type":"markdown","source":"#### Multi-label Classification Data\n\nWe are going to consider our problem as multi-label classification. This is because a patient can have one of five different types of injuries, and the injuries can potentially co-exist.\n\nWe're going to use the patient labels in _train.csv_ to derive labels for each image.","metadata":{}},{"cell_type":"code","source":"filename_labels = pd.DataFrame()\nf_list = os.listdir('/storage/rsna-atd-128-png-pt1')\nfilename_labels['fname'] = pd.Series(f_list)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From _train.csv_ , we construct a dictionary _patient_dict_ which has patient_id as key and the list of injuries as labels.","metadata":{}},{"cell_type":"code","source":"train_patients['is_valid'] = False\nvalidation_patients['is_valid'] = True\ntrain_val_df = pd.concat([train_patients, validation_patients])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_val_df['is_valid'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_dict = {}\n\ntarget_labels = ['bowel_injury', 'extravasation_injury' , 'kidney_healthy' , 'kidney_low' , 'kidney_high', 'liver_healthy', 'liver_low',\n                 'liver_high', 'spleen_healthy', 'spleen_low', 'spleen_high']\nreduced_target_labels = ['bowel_injury', 'extravasation_injury' , 'kidney_injury', 'liver_injury', 'spleen_injury']\n\n\nfor idx, patient_id in enumerate(train_val_df['patient_id']):\n    entry = train_val_df.iloc[idx][target_labels]\n    is_valid = train_val_df.iloc[idx]['is_valid'] \n    patient_dict[patient_id] = (entry, is_valid)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Having constructed our _patient_dict_ dictionary, we can loop through all the PNG filenames and for each look up the appropriate entry from the dictionary using _patient_id_ as the key. We encode bowel and extravasation injury as binary variable. We encode liver, kidney and spleen as 0/1/2 depending on whether there is no injury (healthy), low injury or high injury.\n\nIt was important to construct a dictionary first, so that this lookup can be fast using hash table.","metadata":{}},{"cell_type":"code","source":"label_list = []\nis_valid_list = []\nreduced_target_labels = ['bowel_injury', 'extravasation_injury' , 'kidney_injury', 'liver_injury', 'spleen_injury']\n\nstart = time.time()\nfor scan_name in tqdm(filename_labels['fname']):\n    patient_id = int(scan_name.split('_')[0])\n    entry = patient_dict[patient_id][0]\n    \n    if_bowel = entry['bowel_injury']\n    if_extravasation = entry['extravasation_injury']\n    if_kidney = max(1*entry['kidney_healthy'], 2*entry['kidney_low'], 3*entry['kidney_high']) - 1\n    if_liver = max(1*entry['liver_healthy'], 2*entry['liver_low'], 3*entry['liver_high']) - 1\n    if_spleen = max(1*entry['spleen_healthy'], 2*entry['spleen_low'], 3*entry['spleen_high']) - 1\n\n    labels = [if_bowel, if_extravasation, if_kidney, if_liver, if_spleen]\n    label_list.append(labels)\n    \n    is_valid = patient_dict[patient_id][1]\n    is_valid_list.append(is_valid)\n\nend = time.time()\nprint(end - start)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filename_labels.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filename_labels[reduced_target_labels] = pd.DataFrame(label_list)\nfilename_labels['is_valid'] = pd.Series(is_valid_list)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"assert(filename_labels['bowel_injury'].all() in [0, 1])\nassert(filename_labels['extravasation_injury'].all() in [0, 1])\nassert(filename_labels['kidney_injury'].all() in [0, 1, 2])\nassert(filename_labels['liver_injury'].all() in [0, 1, 2])\nassert(filename_labels['spleen_injury'].all() in [0, 1, 2])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have the data in the format we wanted. Now we can construct the dataloader and train a CNN model.","metadata":{}},{"cell_type":"code","source":"print(filename_labels.shape)\nfilename_labels.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training","metadata":{}},{"cell_type":"markdown","source":"We're going to take a ResNet-18 model pretrained on ImageNet (directly available via fastai library) and fine-tune it on our small data.\n\nFirst, we construct our dataloader. Here, we do not perform any data augmentation.","metadata":{}},{"cell_type":"code","source":"SIZE = 128\n\nbowel_vocab = filename_labels['bowel_injury'].unique()\nextravasation_vocab = filename_labels['extravasation_injury'].unique()\nkidney_vocab = filename_labels['kidney_injury'].unique()\nliver_vocab = filename_labels['liver_injury'].unique()\nspleen_vocab = filename_labels['spleen_injury'].unique()\n\nblocks = (ImageBlock(cls=PILImageBW), \n          CategoryBlock(vocab = bowel_vocab),\n          CategoryBlock(vocab = extravasation_vocab),\n          CategoryBlock(vocab = kidney_vocab),\n          CategoryBlock(vocab = liver_vocab),\n          CategoryBlock(vocab = spleen_vocab))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kidney_vocab","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"getters = (ColReader('fname', pref = '/storage/rsna-atd-128-png-pt1/'), \n           ColReader('bowel_injury'), ColReader('extravasation_injury'),\n           ColReader('kidney_injury'), ColReader('liver_injury'), ColReader('spleen_injury'))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_block = DataBlock(blocks = blocks,\n                       getters = getters,\n                       splitter = ColSplitter('is_valid'),\n                       n_inp = 1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dls = data_block.dataloaders(filename_labels, bs = 128)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dls.c","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let  us visualize one batch. The labels are of the form (y1, y2, y3, y4, y5) where y1 represents bowel_injury and can take 0 or 1. y5 represents spleen_injury and can take 0/1/2.","metadata":{}},{"cell_type":"code","source":"dls.show_batch()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We need to change the head of the ResNet-18 model because we're dealing with a multi-label problem instead of a multi-class problem (for which the network was pretrained on ImageNet. Essentially, we're going to remove the output layer and replace it by one output node for bowel (binary), one for extravasation (binary), three for liver (ternary), three for kidney, and three for spleen. Hence, we have 5 new heads which we are going to connect to the *body* of ResNet18 - body is the full-network minus the layer output layer. ","metadata":{}},{"cell_type":"code","source":"class MultiHeadModel(Module):\n    \n    def __init__(self, body):\n    \n        self.body = body\n        nf = num_features_model(nn.Sequential(*self.body.children()))\n\n        self.bowel = create_head(nf, 1)\n        self.extravasation = create_head(nf, 1)\n        self.kidney = create_head(nf, 3)\n        self.liver = create_head(nf, 3)\n        self.spleen = create_head(nf, 3)\n        \n    def forward(self, x):\n        \n        y = self.body(x)\n        bowel = self.bowel(y)\n        extravasation = self.extravasation(y)\n        kidney = self.kidney(y)\n        liver = self.liver(y)\n        spleen = self.spleen(y)\n        return [bowel, extravasation, kidney, liver, spleen]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_model = create_vision_model(models.resnet18, 10, True, n_in = 1)\nbody = create_body(base_model, pretrained=True)\nnet = MultiHeadModel(body)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our loss is now combination of binary cross-entropy losses for the bowel and extravasation labels, and multi-class cross-entropy loss for kidney, liver, and spleen. Given the sample weights in the competition, we introduce weights for each loss type; hence mistakes in some categories (such as extravasation) are penalized much more than others (such as bowel). For ternary categories, the weight is the average of the weights of the low and high levels ( (2+4)/2 = 3).","metadata":{}},{"cell_type":"code","source":"class CombinationLoss(Module):\n    \"Cross entropy loss on multiple targets\"\n    def __init__(self, weights = [2, 6, 3, 3, 3]):\n        self.w = weights\n        \n    def forward(self, xs, *ys, reduction = 'mean'):\n        loss = 0\n        for i, w, x, y in zip(range(len(xs)), self.w, xs, ys):\n            if i < 2:\n                #loss += w*F.binary_cross_entropy_with_logits(x, y.unsqueeze(1).float(), reduction = reduction)\n                loss += w*torchvision.ops.sigmoid_focal_loss(x, y.unsqueeze(1).float(), reduction = reduction)\n            else:\n                loss += w*F.cross_entropy(x, y, reduction = reduction)\n        return loss","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Because we observed that there is severe imabalance in the data, accuracy is not a good metric to judge the quality of our model. Instead, we are going to look at the per-category recall and the average recall.","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import recall_score\n\nclass RecallPartial(Metric):\n    \"Stores predictions and targets on CPU in accumulate to perform final calculations with `func`.\"\n    def __init__(self, a=0, **kwargs):\n        self.func = partial(recall_score, average='macro', zero_division=0)\n        self.a = a\n\n    def reset(self): self.targs,self.preds = [],[]\n\n    def accumulate(self, learn):\n        if self.a < 2:\n            pred = learn.pred[self.a]>0\n        else:\n            pred = learn.pred[self.a].argmax(-1)\n        targ = learn.y[self.a]\n        pred,targ = to_detach(pred),to_detach(targ)\n        pred,targ = flatten_check(pred,targ)\n        self.preds.append(pred)\n        self.targs.append(targ)\n\n    @property\n    def value(self):\n        if len(self.preds) == 0: return\n        preds,targs = torch.cat(self.preds),torch.cat(self.targs)\n        return self.func(targs, preds)\n\n    @property\n    def name(self): return 'recall_' + filename_labels.columns[self.a+1].split('_')[0]\n    \nclass RecallCombine(Metric):\n    \n    def accumulate(self, learn):\n        scores = [learn.metrics[i].value for i in range(3)]\n        self.combine = np.average(scores, weights=[2,1,1])\n\n    @property\n    def value(self):\n        return self.combine","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We finally create a fast.ai learner object that encapsulates the dataset, model, loss function, and metrics.","metadata":{}},{"cell_type":"code","source":"learn = Learner(dls, net, loss_func = CombinationLoss(), metrics=[RecallPartial(a=i) for i in range(len(dls.c))] + [RecallCombine()],\n               default_cbs = False, cbs = [TrainEvalCallback, Recorder(train_metrics = False, valid_metrics = True), ProgressCallback])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at the default optimization settings used by fast.ai","metadata":{}},{"cell_type":"code","source":"learn.opt_func","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We do mixed-precision training to speed up our training process.","metadata":{}},{"cell_type":"code","source":"learn.to_fp16()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"fastai provides a handy function to let us find a good learning rate.","metadata":{}},{"cell_type":"code","source":"learn.lr_find()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's now fine tune our model for 4 epochs using the recommended learning rate.","metadata":{}},{"cell_type":"code","source":"learn.fine_tune(4, 1.5e-3)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As the epochs progress, the training error decreases to close to 0 (perfect recall for all categories for all images in the training set). However, the validation error is large and does not improve. The recall metrics on the validation set also indicate not very good performance.\n\nPerhaps our idea of assigning the same labels to all images in the CT scan series is not a good one. For instance, any organ is only present in a subset of images and absent in others. Besides, we've discarded the temporal information. Because of these, there is no meaningful pattern to learn and the network simply overfits to the training data.","metadata":{}},{"cell_type":"markdown","source":"### Save Model","metadata":{}},{"cell_type":"markdown","source":"We need to save our model, so that our [inference notebook](https://www.kaggle.com/code/pankajpansari/rsna-2023-atd-baseline-1-inference) can import it to make predictions on the test set.","metadata":{}},{"cell_type":"code","source":"learn.export('model_2.pt')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_dict = {'dataset': 'rsna-atd-128-png-pt1',\n     'model_name': 'model_2',\n     'architecture': 'resnet18',\n     'bs': 128,\n     'opt': 'adam with fastai defaults',\n     'lr': 1.5e-3,\n     'epochs': 4\n}","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import json\nwith open('saved_models.json', 'a') as outfile:\n    json.dump(a_dict, outfile)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#import gc; gc.collect()\n#torch.cuda.empty_cache()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!watch -n 1 nvidia-smi\n#!nvidia-smi --query-gpu=timestamp,pstate,temperature.gpu,utilization.gpu,utilization.memory,memory.total,memory.free,memory.used --format=csv -l 1","metadata":{},"execution_count":null,"outputs":[]}]}