{"cells":[{"metadata":{},"cell_type":"markdown","source":"This is a kernel under construction.\n\nThe goal is to use train data preprocessed in this kernel\n\nhttps://www.kaggle.com/spacelx/2020-pe-preprocessing-train-data\n\nand this kernel\n\nhttps://www.kaggle.com/spacelx/2020-pe-preprocessing-train-table\n\nto train basic convolutional neural networks in Pytorch to label pulmonary embolism in CT scans.\n\nThe train data has been preprocessed to extract lung features and has been resampled such that each study is associated with a 3D scan with 20 slices of 128x128 pixels each. The full preprocessed dataset can be found here\n\nhttps://www.kaggle.com/spacelx/2020pe-preprocessed-train-data\n\nIn the current state this kernel gets the data ready and sets up a CNN intended to only classify 3D scans into (negative / indeterminate / positive) PE exams. A second CNN of the same structure is then trained to classify PE-positive studies into (acute / chronic / acute & chronic) PE exams.\n\nExamining the confusion matrices from both stages, it's clear that none of that actually works at the moment. This is surely due to the preprocessing simplifying the train dataset quite a bit, and certainly also caused by the simplicity of the model used.\n\nIf anyone has any ideas for improvement, feel free to leave a comment! Would be nice if we could turn this into a somewhat-working-but-certainly-not-perfect kernel with a simple approach for fresh starters!"},{"metadata":{},"cell_type":"markdown","source":"# Setup\n\nLet's start by importing some useful packages."},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport time\n\n# path management\nfrom pathlib import Path\n\n# progress bars\nfrom tqdm import tqdm\n\n# plotting\nimport matplotlib.pyplot as plt\nimport seaborn as sn\n\n# pytorch\nimport torch\nfrom torch.utils.data import Dataset\nfrom torch.utils.data import DataLoader\n\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom torch.optim import SGD, Adam\n\n# zipfile management\nimport zipfile","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Saving paths of base data into variables"},{"metadata":{"trusted":true},"cell_type":"code","source":"comp_data_path = Path('../input/rsna-str-pulmonary-embolism-detection')\nprep_data_path = Path('../input/2020pe-preprocessed-train-data')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some settings - the preprocessed dataset is 20 slices of size 128x128 for each study, showing only the lung sections.\n\nThe preprocessing was done in \n\nhttps://www.kaggle.com/spacelx/2020-pe-preprocessing-train-data\n\nwhich is based on this kernel\n\nhttps://www.kaggle.com/allunia/pulmonary-dicom-preprocessing,\n\nand the complete preprocessed dataset can be found at\n\nhttps://www.kaggle.com/spacelx/2020pe-preprocessed-train-data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# set sizing\nNSCANS = 20\nNPX = 128\n\n# set device\ndevice = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we load an example study, just to show what our preprocessed data looks like."},{"metadata":{"trusted":true},"cell_type":"code","source":"allfiles = list((prep_data_path / f'proc_{NSCANS}_{NPX}_train').glob('*_data.npy'))\nsample_file = allfiles[0]\nsample_scans = np.load(str(sample_file), allow_pickle=True)\n\nfig, ax = plt.subplots(5, 4, figsize=(20,20))\nax = ax.flatten()\nfor m in range(NSCANS):\n    ax[m].imshow(sample_scans[m], cmap='Blues_r')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now reading in the preprocessed train data table, which includes one row per study (not per image as in the original version). Preprocessing was done in this kernel\n\nhttps://www.kaggle.com/spacelx/2020-pe-preprocessing-train-table.\n\nThere are 20 columns (one for each resampled slice) noting whether there was PE present in the base data or not."},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(prep_data_path / 'train_proc.csv', index_col=0)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dataset model\n\nJust quickly setting up a Pytorch dataset..."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class PE2020Dataset(Dataset):\n    def __init__(self, scans, labels, datapath):\n        self.scans = scans\n        self.labels = labels\n        self.datapath = Path(datapath)\n\n    def __len__(self):\n        return len(self.scans)\n\n    def __getitem__(self, i):\n        file = self.datapath / (self.scans[i] + '_data.npy')\n        x = torch.tensor(np.load(str(file)), dtype=torch.float).to(device)\n        y = torch.tensor(self.labels[i]).to(device)\n        return x, y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model\n\nWe define a basic CNN with two convolutional and two fully connected linear layers. We will use this several times later on.\nThe constructor parameter <code>NOUT</code> is used to set the number of output channels."},{"metadata":{"trusted":true},"cell_type":"code","source":"class PE2020Net(nn.Module):\n    def __init__(self, NOUT):\n        super(PE2020Net, self).__init__()\n        self.conv1 = nn.Conv2d(NSCANS, 32, 5)\n        self.conv2 = nn.Conv2d(32, 64, 5)\n        self.conv2_drop = nn.Dropout2d()\n        self.fc1 = nn.Linear(64*29*29, 500)\n        #self.fc2 = nn.Linear(5000, 5000)\n        self.fc3 = nn.Linear(500, 3)\n\n    def forward(self, x):\n        x = x.view(-1, NSCANS, NPX, NPX)\n        x = self.conv1(x)\n        x = F.max_pool2d(x, 2)\n        x = F.relu(x)\n        x = self.conv2(x)\n        x = self.conv2_drop(x)\n        x = F.max_pool2d(x, 2)\n        x = F.relu(x)\n        x = x.view(-1, 64*29*29)\n        x = self.fc1(x)\n        x = F.dropout(x, training=self.training)\n        x = F.relu(x)\n        #x = self.fc2(x)\n        #x = F.relu(x)\n        x = self.fc3(x)\n        x = F.log_softmax(x, dim=-1)\n        return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# First stage - classification into (PE / no PE / indeterminate)"},{"metadata":{},"cell_type":"markdown","source":"Set up a dataset with all scans and labels <code>positive_exam_for_pe</code>,  <code>negative_exam_for_pe</code> and  <code>indeterminate</code>, and create DataLoaders."},{"metadata":{"trusted":true},"cell_type":"code","source":"# separate scans from labels\nscans = train['dcmpath']\nall_labels = train.drop(labels='dcmpath', axis=1).astype(int)\nlabels = all_labels[['negative_exam_for_pe', 'indeterminate']].copy()\nlabels['positive_exam_for_pe'] = 1 - labels[['negative_exam_for_pe', 'indeterminate']].sum(axis=1)\n\n# keep label names for later\nlabel_names = labels.columns.tolist()\n# get label index for all studies\nlabels = np.where(labels.values)[1]\n\n# split into train and validation and create datasets\ntmp = train.sample(len(train)).index.values\nidx_train = tmp[:int(0.8*len(tmp))]\nidx_valid = tmp[int(0.8*len(tmp)):]\ntrainset = PE2020Dataset(scans.loc[idx_train].values, labels[idx_train], prep_data_path / f'proc_{NSCANS}_{NPX}_train')\nvalidset = PE2020Dataset(scans.loc[idx_valid].values, labels[idx_valid], prep_data_path / f'proc_{NSCANS}_{NPX}_train')\n\n# instantiate DataLoaders\nBATCHSIZE = 20\ntrainloader = DataLoader(trainset, batch_size=BATCHSIZE, shuffle=True)\nvalidloader = DataLoader(validset, batch_size=BATCHSIZE, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Instantiate model and train, validate after each epoch."},{"metadata":{"trusted":true},"cell_type":"code","source":"model = PE2020Net(3).to(device)\noptimiser = Adam(model.parameters(), lr=0.005)\nepochs = 20\naccuracy = {}\naccuracy['train'] = []\naccuracy['valid'] = []\n\nfor epoch in range(epochs):\n    \n    # training\n    correct = 0\n    losses = torch.tensor([])\n    for batch_data in tqdm(trainloader):\n        X, y = batch_data\n        # zero gradients\n        model.zero_grad()\n        # forward pass\n        output = model(X)\n        \n        # count accurate predictions\n        predicted = torch.max(output.detach(),1)[1]\n        correct += (predicted == y).sum()\n        \n        # calculate loss and keep it for later\n        loss = F.cross_entropy(output, y)\n        losses = torch.cat((losses, torch.tensor([loss.detach()])), 0)\n        \n        # backward pass\n        loss.backward()\n        # update parameters\n        optimiser.step()\n    \n    # calculate mean loss and accuracy and print\n    meanacc = float(correct) / (len(trainloader) * BATCHSIZE)\n    meanloss = float(losses.mean())\n    print('Epoch:', epoch, 'Loss:', meanloss, 'Accuracy:', meanacc)\n    accuracy['train'].append(meanacc)\n    # putting this in to keep the console clean\n    time.sleep(0.5)\n    \n    \n    \n    # validation\n    correct = 0\n    for batch_data in tqdm(validloader):\n        # forward pass\n        X, y = batch_data\n        with torch.no_grad():\n            output = model(X)\n        # get number of accurate predictions\n        predicted = torch.max(output,1)[1]\n        correct += (predicted == y).sum()\n    # calculate mean accuracy and print\n    meanacc = float(correct) / (len(validloader) * BATCHSIZE)\n    print('Validation epoch:', epoch, 'Accuracy:', meanacc)\n    accuracy['valid'].append(meanacc)\n    # putting this in to keep the console clean\n    time.sleep(0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plotting accuracy per epoch here - we don't see much improvement. Our model gets a bit over 60% right, which doesn't seem totally bad but also not really awesome. But that is to be expected; we boiled the dataset down by quite a bit and our model isn't really advanced either."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(accuracy['train'],label=\"Training Accuracy\")\nplt.plot(accuracy['valid'],label=\"Validation Accuracy\")\nplt.xlabel('No. of Epochs')\nplt.ylabel('Accuracy')\nplt.legend(frameon=False)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's quickly save the stage 1 model"},{"metadata":{"trusted":true},"cell_type":"code","source":"torch.save(model, 'stg1_model.pt')\nnp.save('stg1_label_names', label_names)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we'll check our predictions with a confusion matrix"},{"metadata":{"trusted":true},"cell_type":"code","source":"validset = iter(validloader)\nall_y = np.array([])\nall_y_pred = np.array([])\nfor x, y in validset:\n    with torch.no_grad():\n        y_pred = model(x)\n    all_y = np.append(all_y, y.tolist())\n    all_y_pred = np.append(all_y_pred, torch.max(y_pred,1)[1].tolist())\n    \ndf = pd.DataFrame(np.array([all_y, all_y_pred]).T, columns=['y','y_pred'])\nconfusion_matrix = pd.crosstab(df['y'], df['y_pred'], rownames=['y'], colnames=['y_pred'])\n\nsn.heatmap(confusion_matrix, annot=True)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A few things to unpack here. First of all, our model never predicts <code>indeterminate</code> - there are barely any examples of that kind in the set. Secondly, our model pretty much always predicts <code>negative_exam_for_pe</code>, and the 60-something percent accuracy we get simply reflects the 60% of PE-negative studies in the dataset.\nSo all in all we can safely say our model doesn't predict shit for now..."},{"metadata":{},"cell_type":"markdown","source":"# Second stage - for PE samples, classify into (acute / chronic / acute & chronic)\nLet's try something else - now we use the same model architecture to classify PE-positive studies by the PE type, so we only use PE-positive studies for training obviously. Apart from that, the steps are the same as in the previous stage."},{"metadata":{"trusted":true},"cell_type":"code","source":"# choose only PE-positive samples\ntrain_pe_positive = train[(train['negative_exam_for_pe'] == 0) &\n                          (train['indeterminate'] == 0)]\n\n# separate scans from labels\nscans = train_pe_positive['dcmpath']\nall_labels = train_pe_positive.drop(labels='dcmpath', axis=1).astype(int)\nlabels = all_labels[['acute_pe', 'chronic_pe', 'acute_and_chronic_pe']].copy()\n\n# keep label names for later\nlabel_names = labels.columns.tolist()\n# get label index for all studies\nlabels['label'] = np.where(labels.values)[1]\n\n# split into train and validation and create datasets\ntmp = train_pe_positive.sample(len(train_pe_positive)).index.values\nidx_train = tmp[:int(0.8*len(tmp))]\nidx_valid = tmp[int(0.8*len(tmp)):]\ntrainset = PE2020Dataset(scans.loc[idx_train].values, labels.loc[idx_train, 'label'].values,\n                         prep_data_path / f'proc_{NSCANS}_{NPX}_train')\nvalidset = PE2020Dataset(scans.loc[idx_valid].values, labels.loc[idx_train, 'label'].values,\n                         prep_data_path / f'proc_{NSCANS}_{NPX}_train')\n\n# instantiate DataLoaders\nBATCHSIZE = 20\ntrainloader = DataLoader(trainset, batch_size=BATCHSIZE, shuffle=True)\nvalidloader = DataLoader(validset, batch_size=BATCHSIZE, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We set up the same model as before and train it, this time on PE-positive studies only..."},{"metadata":{"trusted":true},"cell_type":"code","source":"model = PE2020Net(3).to(device)\noptimiser = Adam(model.parameters(), lr=0.005)\nepochs = 20\naccuracy = {}\naccuracy['train'] = []\naccuracy['valid'] = []\n\nfor epoch in range(epochs):\n    \n    # training\n    correct = 0\n    losses = torch.tensor([])\n    for batch_data in tqdm(trainloader):\n        X, y = batch_data\n        # zero gradients\n        model.zero_grad()\n        # forward pass\n        output = model(X)\n        \n        # count accurate predictions\n        predicted = torch.max(output.detach(),1)[1]\n        correct += (predicted == y).sum()\n        \n        # calculate loss and keep it for later\n        loss = F.cross_entropy(output, y)\n        losses = torch.cat((losses, torch.tensor([loss.detach()])), 0)\n        \n        # backward pass\n        loss.backward()\n        # update parameters\n        optimiser.step()\n    \n    # calculate mean loss and accuracy and print\n    meanacc = float(correct) / (len(trainloader) * BATCHSIZE)\n    meanloss = float(losses.mean())\n    print('Epoch:', epoch, 'Loss:', meanloss, 'Accuracy:', meanacc)\n    accuracy['train'].append(meanacc)\n    # putting this in to keep the console clean\n    time.sleep(0.5)\n    \n    \n    \n    # validation\n    correct = 0\n    for batch_data in tqdm(validloader):\n        # forward pass\n        X, y = batch_data\n        with torch.no_grad():\n            output = model(X)\n        # get number of accurate predictions\n        predicted = torch.max(output,1)[1]\n        correct += (predicted == y).sum()\n    # calculate mean accuracy and print\n    meanacc = float(correct) / (len(validloader) * BATCHSIZE)\n    print('Validation epoch:', epoch, 'Accuracy:', meanacc)\n    accuracy['valid'].append(meanacc)\n    # putting this in to keep the console clean\n    time.sleep(0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We plot the accuracy... looks like we're overfitting a bit, training accuracy increases while validation accuracy doesn't..."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(accuracy['train'],label=\"Training Accuracy\")\nplt.plot(accuracy['valid'],label=\"Validation Accuracy\")\nplt.xlabel('No. of Epochs')\nplt.ylabel('Accuracy')\nplt.legend(frameon=False)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"... and save the stage 2 model."},{"metadata":{"trusted":true},"cell_type":"code","source":"torch.save(model, 'stg2_model.pt')\nnp.save('stg2_label_names', label_names)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's plot a confusion matrix again"},{"metadata":{"trusted":true},"cell_type":"code","source":"validset = iter(validloader)\nall_y = np.array([])\nall_y_pred = np.array([])\nfor x, y in validset:\n    with torch.no_grad():\n        y_pred = model(x)\n    all_y = np.append(all_y, y.tolist())\n    all_y_pred = np.append(all_y_pred, torch.max(y_pred,1)[1].tolist())\n    \ndf = pd.DataFrame(np.array([all_y, all_y_pred]).T, columns=['y','y_pred'])\nconfusion_matrix = pd.crosstab(df['y'], df['y_pred'], rownames=['y'], colnames=['y_pred'])\n\nsn.heatmap(confusion_matrix, annot=True)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"label_names","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking at this, it is clear that nearly all PE-positive studies in the validation set (and so likely also in the training set) are acute, while nearly none are chronic or acute & chronic - we have a highly imbalanced dataset. No wonder our model has a hard time identifying chronic and acute & chronic studies. But then remembering stage 1 above, the model didn't fare much better even when the dataset was balanced..."}],"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}