{"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":"# Peek\n\nFirst, let's see the structure of our files.","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport gc\nimport time\nimport matplotlib.pyplot as plt\nfrom IPython.display import clear_output\nfrom tensorflow.keras.models import load_model\nfrom tensorflow.keras.callbacks import ModelCheckpoint as MC\nfrom tensorflow.keras import backend as K\nfrom sklearn.metrics import confusion_matrix, roc_curve, auc, recall_score, accuracy_score, balanced_accuracy_score, precision_score\nimport pickle\nfrom IPython.display import FileLink\n\n\n\nroot = '/kaggle/input/rsna-str-pulmonary-embolism-detection'\nfor item in os.listdir(root):\n    path = os.path.join(root, item)\n    if os.path.isfile(path):\n        print(path)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-31T21:21:21.013138Z","iopub.execute_input":"2022-01-31T21:21:21.013549Z","iopub.status.idle":"2022-01-31T21:21:21.032195Z","shell.execute_reply.started":"2022-01-31T21:21:21.013508Z","shell.execute_reply":"2022-01-31T21:21:21.031517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data\n\nNext, we load all '.csv' files into memory and peek into their makeup.","metadata":{}},{"cell_type":"code","source":"print('Reading train data...')\ndata = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\nprint(data.shape)\ndata.head()","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","execution":{"iopub.status.busy":"2022-01-31T21:20:36.333616Z","iopub.execute_input":"2022-01-31T21:20:36.333959Z","iopub.status.idle":"2022-01-31T21:20:39.762678Z","shell.execute_reply.started":"2022-01-31T21:20:36.333919Z","shell.execute_reply":"2022-01-31T21:20:39.76198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Plot EDA\nfor col in data.columns[3:len(data.columns)]:\n    proportion=[sum(data[col]), len(data)-sum(data[col])]\n    proportion=pd.Series(proportion)\n    #colors\n    colors = ['#A9C4EB','#10739E']\n    #explsion\n    explode = (0.05,0.05)\n    proportion.plot(kind='pie',colors = colors, labels=[True,False], \n        autopct='%1.1f%%', startangle=90, pctdistance=0.85, \n        explode = explode)\n#     plt.pie(proportion, colors = colors, labels=[True,False], \n#             autopct='%1.1f%%', startangle=90, pctdistance=0.85, \n#             explode = explode)\n    plt.title(col)\n    #draw circle\n    centre_circle = plt.Circle((0,0),0.70,fc='white')\n    fig = plt.gcf()\n    fig.gca().add_artist(centre_circle)\n    # Equal aspect ratio ensures that pie is drawn as a circle\n#     ax1.axis('equal')  \n    plt.tight_layout()\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train/Test Split\nSince this project is for academic purposes, we will split the give train set into train+test set, to evaluate performance of the model better.\nWe will split by StudyInstanceID which is the exam ID, so that images from the same exam don't show up in train and test.","metadata":{}},{"cell_type":"code","source":"# Separate the labels from the test set and save for later\nkeys=['negative_exam_for_pe', 'qa_motion',\n       'qa_contrast', 'flow_artifact', 'rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1',\n       'leftsided_pe', 'chronic_pe', 'true_filling_defect_not_pe',\n       'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate']\n\ny=data['pe_present_on_image'].copy()\nx=data.drop('pe_present_on_image',axis=1)\n\nfor key in keys:\n    y=pd.concat([y,x[key]],axis=1)\n    x=x.drop(key, axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:35:57.15789Z","iopub.execute_input":"2022-01-31T16:35:57.1585Z","iopub.status.idle":"2022-01-31T16:35:59.678553Z","shell.execute_reply.started":"2022-01-31T16:35:57.158459Z","shell.execute_reply":"2022-01-31T16:35:59.67776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import GroupShuffleSplit,StratifiedShuffleSplit, GroupKFold\n\n#train,test=train_test_split(data, test_size=0.2, random_state=42, shuffle=False)\n\n#This split allows train and test sets to not have images from the same exam\nsplitter = GroupShuffleSplit(test_size=.20, random_state=42, n_splits=5)\n# splitter = StratifiedShuffleSplit(test_size=.20, random_state=42, n_splits=2) #this split would shuffle images, both were tried but no better results\nsplit = splitter.split(x,y, groups=data.StudyInstanceUID)\ntrain_inds, test_inds = next(split)\n\ntrain = data.iloc[train_inds]\ntest = data.iloc[test_inds]","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:35:59.679869Z","iopub.execute_input":"2022-01-31T16:35:59.680222Z","iopub.status.idle":"2022-01-31T16:36:02.502591Z","shell.execute_reply.started":"2022-01-31T16:35:59.680184Z","shell.execute_reply":"2022-01-31T16:36:02.501668Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()\n# print('%PE positive images:', round(sum(train.pe_present_on_image)/len(train)*100),'%')","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:02.506393Z","iopub.execute_input":"2022-01-31T16:36:02.506762Z","iopub.status.idle":"2022-01-31T16:36:02.522198Z","shell.execute_reply.started":"2022-01-31T16:36:02.506723Z","shell.execute_reply":"2022-01-31T16:36:02.521174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.head()\ny_test=test['pe_present_on_image'].copy()\ntest=test.drop('pe_present_on_image',axis=1)\n\nfor key in keys:\n    y_test=pd.concat([y_test,test[key]],axis=1)\n    test=test.drop(key, axis=1)\n    \n\nprint('%PE positive images:', round(sum(y_test.pe_present_on_image)/len(train)*100),'%\\n',\n      '# of PE positive images:', sum(y_test.pe_present_on_image))","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:02.524076Z","iopub.execute_input":"2022-01-31T16:36:02.524662Z","iopub.status.idle":"2022-01-31T16:36:03.090717Z","shell.execute_reply.started":"2022-01-31T16:36:02.524623Z","shell.execute_reply":"2022-01-31T16:36:03.089977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:03.092402Z","iopub.execute_input":"2022-01-31T16:36:03.092821Z","iopub.status.idle":"2022-01-31T16:36:03.107815Z","shell.execute_reply.started":"2022-01-31T16:36:03.092771Z","shell.execute_reply":"2022-01-31T16:36:03.105548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After checking the keys, we will now work on the fuction to get the image array from a DICOM image. We will be using a code snippet by [eladwar](https://www.kaggle.com/eladwar) from his notebook [here](https://www.kaggle.com/eladwar/20-seconds-or-less).","metadata":{}},{"cell_type":"code","source":"import vtk\nfrom vtk.util import numpy_support\nimport cv2\nimport albumentations as albu\n\nreader = vtk.vtkDICOMImageReader()\n\ntrain_transform = [albu.RandomBrightnessContrast(p=0.3),\n                   albu.VerticalFlip(p=0.5),\n                   albu.HorizontalFlip(p=0.5),\n                   albu.Downscale(p=1.0,scale_min=0.35,scale_max=0.75,),\n                   albu.Resize(512, 512)]\ntransform=albu.Compose(train_transform)\n\ndef normalize_image(image):\n    min = -1000\n    max = 400\n    image[image < min] = min\n    image[image > max] = max\n#     image = (image - min) / (max - min)\n#     image = image.astype(\"float64\")\n    return image\n\ndef get_img(path):\n    reader.SetFileName(path)\n    reader.Update()\n    _extent = reader.GetDataExtent()\n    ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]\n\n    ConstPixelSpacing = reader.GetPixelSpacing()\n    imageData = reader.GetOutput()\n    pointData = imageData.GetPointData()\n    arrayData = pointData.GetArray(0)\n    ArrayDicom = numpy_support.vtk_to_numpy(arrayData)\n    ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')\n    ArrayDicom = cv2.resize(ArrayDicom,(512,512))\n    \n    ArrayDicom = normalize_image(ArrayDicom)\n    return ArrayDicom","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:03.109031Z","iopub.execute_input":"2022-01-31T16:36:03.109356Z","iopub.status.idle":"2022-01-31T16:36:05.167683Z","shell.execute_reply.started":"2022-01-31T16:36:03.109331Z","shell.execute_reply":"2022-01-31T16:36:05.166906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After defining our image reader, we will test it with a sample DICOM image to load.","metadata":{}},{"cell_type":"code","source":"#test read a dcom file and view it\nfpath = \"../input/rsna-str-pulmonary-embolism-detection/train/6897fa9de148/2bfbb7fd2e8b/be0b7524ffb4.dcm\"\nds = get_img(fpath)\n\nimport matplotlib.pyplot as plt\n\n#Convert dcom file to 8bit color\nfunc = lambda x: int((2**15 + x)*(255/2**16))\nint16_to_uint8 = np.vectorize(func)\n\ndef show_dicom_images(dcom,transform_img=False):\n    f, ax = plt.subplots(1,2, figsize=(16,20))\n    if transform_img:\n        data_row_img = int16_to_uint8(dcom)\n        data_row_img=transform(image=data_row_img)['image']\n        title1='8-bit DICOM Transformed Image'\n    else:\n        data_row_img = int16_to_uint8(dcom)\n        title1='8-bit DICOM Image'\n        \n    ax[0].imshow(data_row_img, cmap=plt.cm.bone)\n    ax[1].imshow(ds, cmap=plt.cm.bone)\n    #print(data_row_img)\n    ax[0].axis('off')\n    ax[0].set_title(title1)\n    ax[1].axis('off')\n    ax[1].set_title('16-bit DICOM Image')\n    plt.show()\n\n\nshow_dicom_images(ds,False)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:05.168843Z","iopub.execute_input":"2022-01-31T16:36:05.169195Z","iopub.status.idle":"2022-01-31T16:36:05.569518Z","shell.execute_reply.started":"2022-01-31T16:36:05.169159Z","shell.execute_reply":"2022-01-31T16:36:05.568654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Creation","metadata":{}},{"cell_type":"markdown","source":"I decided to try and use a pre-trained Xception model as a feature-extractor. To simplify my coding, I did a multi-output model with each output having one node activated by a sigmoid. Also, since we are dealing with numbers between 1 and 0, I decided to use binary_crossentropy as the loss function.","metadata":{}},{"cell_type":"markdown","source":"**Embedding Model Creation**","metadata":{}},{"cell_type":"code","source":"# import tensorflow as tf\n# from tensorflow import keras\n# from tensorflow.keras.models import Model, Sequential\n# from tensorflow.keras.layers import Input, Dense, Dropout, Conv2D, MaxPooling2D, Flatten, GlobalAveragePooling2D, Embedding, Dot, Reshape, Multiply\n\n# inputs = Input((512, 512, 3))\n# #x = Conv2D(3, (1, 1), activation='relu')(inputs)\n# base_model = keras.applications.Xception(\n#     include_top=False,\n#     weights=\"imagenet\"\n# )\n\n# base_model.trainable = False\n\n# # First Branch of Neural Network - Classify presense of PE\n# base_outputs = base_model(inputs, training=False)\n# pool_outputs = keras.layers.GlobalAveragePooling2D()(base_outputs)\n# drop_outputs = Dropout(0.25)(pool_outputs)\n# outputs1_dense1 = Dense(1024, activation='relu')(drop_outputs)\n# outputs1_dense2 = Dense(256, activation='relu')(outputs1_dense1)\n# outputs1_dense3 = Dense(64, activation='relu')(outputs1_dense2)\n\n# #First Output\n# ppoi = Dense(1, activation='sigmoid', name='pe_present_on_image')(outputs1_dense3)\n# indt = Dense(1, activation='sigmoid', name='indeterminate')(outputs1_dense3)\n# # Second Branch of Neural Network - Classify severity and position of PE, based on the presence of PE\n\n# #Embedding labels of pe prediction into the input of this model - the embedding weights will be learned during\n# #backpropagation\n# pe_input = Input((None, 1))\n# embedding_size = 2048\n\n# #input_dim is the batch_size which in this case is 1000*0.9 (0.10 is for validation)\n# pe_embedding = Embedding(name='pe_embedding',input_dim = 900, output_dim = embedding_size, input_length=1)(ppoi)\n# # multi=Multiply()([outputs1_dense2,pe_embedding])\n# # merged=Dot(axes=1)([multi, multi])\n# reshape=Reshape(target_shape = [embedding_size])(pe_embedding)\n# merged=reshape*drop_outputs\n# #Merging the embeddings\n# #merged = Dot(name = 'dot_product', normalize = True, axes = 0)([drop_outputs,pe_embedding])\n# outputs2_dense1 = Dense(256, activation='relu')(merged)\n# outputs2_dense2 = Dense(128, activation='relu')(outputs2_dense1)\n# outputs2_dense3 = Dense(64, activation='relu')(outputs2_dense2)\n# #Outputs\n# lspe = Dense(1, activation='sigmoid', name='leftsided_pe')(outputs2_dense3)\n# rspe = Dense(1, activation='sigmoid', name='rightsided_pe')(outputs2_dense3)\n# cnpe = Dense(1, activation='sigmoid', name='central_pe')(outputs2_dense3)\n# cpe = Dense(1, activation='sigmoid', name='chronic_pe')(outputs2_dense3)\n# aacpe = Dense(1, activation='sigmoid', name='acute_and_chronic_pe')(outputs2_dense3)\n\n# model = Model(inputs=inputs, outputs={'pe_present_on_image':ppoi,\n#                                       'leftsided_pe':lspe,\n#                                       'chronic_pe':cpe,\n#                                       'rightsided_pe':rspe,\n#                                       'acute_and_chronic_pe':aacpe,\n#                                       'central_pe':cnpe,\n#                                       'indeterminate':indt})\n\n# opt = keras.optimizers.Adam(lr=0.00001)\n\n\n# model.compile(optimizer=opt,\n#               loss='binary_crossentropy',\n#               metrics=tf.keras.metrics.AUC())\n# model.summary()\n# model.save('pe_detection_model.h5')\n# del model\n# K.clear_session()\n# gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:05.570685Z","iopub.execute_input":"2022-01-31T16:36:05.571004Z","iopub.status.idle":"2022-01-31T16:36:10.319341Z","shell.execute_reply.started":"2022-01-31T16:36:05.570972Z","shell.execute_reply":"2022-01-31T16:36:10.318629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create an embedding of the PE prediction (present/not present) and the image itself, that will be inputs to our second model, which will predict location and severity","metadata":{}},{"cell_type":"markdown","source":"**VGG Model Creation**","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras.models import Model, Sequential\nfrom tensorflow.keras.layers import Input, Dense, Dropout, Conv2D, MaxPooling2D, Flatten, GlobalAveragePooling2D, Embedding, Dot, Reshape, Multiply\ninputs = Input((512, 512, 3))\nvgg=keras.applications.VGG16(\n    include_top=False,\n    weights=None)\n\n# First Branch of Neural Network - Classify presense of PE\nvgg_outputs = vgg(inputs, training=True)\npool_outputs = keras.layers.GlobalAveragePooling2D()(vgg_outputs)\ndrop_outputs = Dropout(0.25)(pool_outputs)\noutputs1_dense1 = Dense(1024, activation='relu')(drop_outputs)\noutputs1_dense2 = Dense(256, activation='relu')(outputs1_dense1)\noutputs1_dense3 = Dense(64, activation='relu')(outputs1_dense2)\n\n#Outputs\nppoi = Dense(1, activation='sigmoid', name='pe_present_on_image')(outputs1_dense3)\nindt = Dense(1, activation='sigmoid', name='indeterminate')(outputs1_dense3)\nlspe = Dense(1, activation='sigmoid', name='leftsided_pe')(outputs1_dense3)\nrspe = Dense(1, activation='sigmoid', name='rightsided_pe')(outputs1_dense3)\ncnpe = Dense(1, activation='sigmoid', name='central_pe')(outputs1_dense3)\ncpe = Dense(1, activation='sigmoid', name='chronic_pe')(outputs1_dense3)\naacpe = Dense(1, activation='sigmoid', name='acute_and_chronic_pe')(outputs1_dense3)\n\nmodel = Model(inputs=inputs, outputs={'pe_present_on_image':ppoi,\n                                      'leftsided_pe':lspe,\n                                      'chronic_pe':cpe,\n                                      'rightsided_pe':rspe,\n                                      'acute_and_chronic_pe':aacpe,\n                                      'central_pe':cnpe,\n                                      'indeterminate':indt})\n\nopt = keras.optimizers.Adam(lr=0.00001)\n\n\nmodel.compile(optimizer=opt,\n              loss='binary_crossentropy',\n              metrics=tf.keras.metrics.AUC())\nmodel.summary()\nmodel.save('pe_detection_model.h5')\ndel model\nK.clear_session()\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:10.320689Z","iopub.execute_input":"2022-01-31T16:36:10.321026Z","iopub.status.idle":"2022-01-31T16:36:10.326253Z","shell.execute_reply.started":"2022-01-31T16:36:10.320989Z","shell.execute_reply":"2022-01-31T16:36:10.325121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training Model","metadata":{}},{"cell_type":"markdown","source":"Before we can train our model, we would be needing an image generator. This is so that our training code would look much cleaner and nicer.","metadata":{}},{"cell_type":"markdown","source":"*custom_dcom_image_generator* is a function that constructs the batch on which we will train and test.\nIf train, it will yield (return a generator, iterable only once) with the X and Y, where X is a group of images\nIf test, it will yield the X on which we will predict.\n\nBy using yield, we only iterate through the batch once and don't save it in memory.","metadata":{}},{"cell_type":"code","source":"def convert_to_rgb(array):\n    array = array.reshape((512, 512, 1))\n    return np.stack([array, array, array], axis=2).reshape((512, 512, 3))\n    \ndef custom_dcom_image_generator(batch_size, dataset, test=False, debug=False):\n    \n    fnames = dataset[['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID']]\n    \n    if not test:\n#         Y=dataset[classes]\n        Y = dataset[['pe_present_on_image', 'leftsided_pe',\n                     'chronic_pe', 'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate'\n                    ]]\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/train'\n        \n    else:\n        prefix = 'input/rsna-str-pulmonary-embolism-detection/train'\n    X = []\n    batch = 0\n    for st, sr, so in fnames.values:\n        if debug:\n            print(f\"Current file: ../{prefix}/{st}/{sr}/{so}.dcm\")\n\n        dicom = get_img(f\"../{prefix}/{st}/{sr}/{so}.dcm\")\n        \n#         eightbit=int16_to_uint8(dicom)\n#         transformed_img=transform(image=dicom)['image']\n#         transformed_img=convert_to_rgb(transformed_img)\n        \n#         X.append(transformed_img)\n\n        image = convert_to_rgb(dicom)\n        X.append(image)\n            \n        del st, sr, so\n        \n        #If we reached the end of the batch\n        if len(X) == batch_size:\n            if test:\n                #yield is used to save memory\n                yield np.array(X)\n                del X\n            else:\n                yield np.array(X), Y[batch*batch_size:(batch+1)*batch_size].values\n                del X\n                \n            gc.collect()\n            X = []\n            batch += 1\n        \n    if test:\n        yield np.array(X)\n    else:\n        yield np.array(X), Y[batch*batch_size:(batch+1)*batch_size].values\n        del Y\n    del X\n    gc.collect()\n    return","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:10.328729Z","iopub.execute_input":"2022-01-31T16:36:10.328978Z","iopub.status.idle":"2022-01-31T16:36:10.345344Z","shell.execute_reply.started":"2022-01-31T16:36:10.328954Z","shell.execute_reply":"2022-01-31T16:36:10.34448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we will be training our model with the train data. We will be using a small batch, fitted with 3 epochs (arbitrary choice) of training to minimize our RAM usage. This is due to our large model, which will eat up quite a large portion of our RAM.\n\nWe will also use sampling for our train data. This will shuffle the data.","metadata":{}},{"cell_type":"code","source":"def train_model(model_path, train_data, train_size, batch_size, max_train_time, debug):\n    #Train loop\n    for n, (x, y) in enumerate(custom_dcom_image_generator(batch_size, train_data, False, debug)):\n        \n\n        if len(x) < 10: #Tries to filter out empty or short data\n            break\n\n        clear_output(wait=True)\n        print(\"Training batch: %i - %i\" %(batch_size*n, batch_size*(n+1)))\n        model = load_model(model_path)\n        hist = model.fit(\n            x[:train_size], \n            #Y values are in a dict as there's more than one target for training output\n            {'pe_present_on_image':y[:train_size, 0],\n             'leftsided_pe':y[:train_size, 1],\n             'chronic_pe':y[:train_size, 2],\n             'rightsided_pe':y[:train_size, 3],\n             'acute_and_chronic_pe':y[:train_size, 4],\n             'central_pe':y[:train_size, 5],\n             'indeterminate':y[:train_size, 6]},\n\n            callbacks = checkpoint,\n\n            validation_split=0.2,\n            epochs=3,\n            batch_size=8,\n            verbose=debug\n        )\n\n        print(\"Metrics for batch validation:\")\n        model.evaluate(x[train_size:],\n                       {'pe_present_on_image':y[train_size:, 0],\n                        'leftsided_pe':y[train_size:, 1],\n                        'chronic_pe':y[train_size:, 2],\n                        'rightsided_pe':y[train_size:, 3],\n                        'acute_and_chronic_pe':y[train_size:, 4],\n                        'central_pe':y[train_size:, 5],\n                        'indeterminate':y[train_size:, 6]\n                       }\n                      )\n\n        try:\n            for key in hist.history.keys():\n                history[key] = np.concatenate([history[key], hist.history[key]], axis=0)\n        except:\n            for key in hist.history.keys():\n                history[key] = hist.history[key]\n\n        #To make sure that our model don't train overtime\n        if time.time() - start >= max_train_time:\n            print(\"Time's up!\")\n            break\n\n        model.save('pe_detection_model.h5')\n        del model, x, y, hist\n        K.clear_session()\n        gc.collect()\n    \n    \n    return history, model","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:10.348313Z","iopub.execute_input":"2022-01-31T16:36:10.348569Z","iopub.status.idle":"2022-01-31T16:36:10.364769Z","shell.execute_reply.started":"2022-01-31T16:36:10.348543Z","shell.execute_reply":"2022-01-31T16:36:10.363963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = {}\nstart = time.time()\ntrain_data=train.sample(frac=1)\nbatch_size = 1000\ndebug = 0\n#90% for training, 10% for validation\ntrain_size = int(batch_size*0.9)\n\nmax_train_time = 3600 * 2 #hours to seconds of training\n\ncheckpoint = MC(filepath='../working/pe_detection_model.h5', monitor='val_loss', save_best_only=True, verbose=1)\n#Train loop\nhistory,trained_model = train_model('../working/pe_detection_model.h5', train_data, train_size, batch_size, max_train_time, debug)\ntrained_model.save('pe_detection_model_trained.h5')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now look at the history of the training for our data. Since the data is trained across several different batches, there should be some spikes among the values reflected here.","metadata":{}},{"cell_type":"code","source":"from IPython.display import FileLink\nimport csv, pickle\n\nFileLink(r'pe_detection_model_trained.h5')","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:35:37.327818Z","iopub.execute_input":"2022-01-31T16:35:37.328146Z","iopub.status.idle":"2022-01-31T16:35:37.336814Z","shell.execute_reply.started":"2022-01-31T16:35:37.328113Z","shell.execute_reply":"2022-01-31T16:35:37.335823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open('history.pkl', 'wb') as f:\n    pickle.dump(history, f)\n    \nFileLink(r'history.pkl')","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:36:10.366105Z","iopub.execute_input":"2022-01-31T16:36:10.366822Z","iopub.status.idle":"2022-01-31T16:36:10.66812Z","shell.execute_reply.started":"2022-01-31T16:36:10.366783Z","shell.execute_reply":"2022-01-31T16:36:10.666866Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Import history file if needed\nwith open('../input/models/hist_file_embed_2048_lr0.00001_norm.pkl', 'rb') as f:\n    history=pickle.load(f)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:50:22.451805Z","iopub.execute_input":"2022-01-31T16:50:22.45217Z","iopub.status.idle":"2022-01-31T16:50:22.465512Z","shell.execute_reply.started":"2022-01-31T16:50:22.452139Z","shell.execute_reply":"2022-01-31T16:50:22.464839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for key in history.keys():\n    if key.startswith('val'):\n        continue\n    else:\n        epoch = range(len(history[key]))\n        plt.plot(epoch, history[key]) #X=epoch, Y=value\n        plt.plot(epoch, history['val_'+key])\n        plt.title(key)\n        if 'accuracy' in key:\n            plt.axis([0, len(history[key]), -0.1, 1.1]) #Xmin, Xmax, Ymin, Ymax\n        plt.legend(['train', 'validation'], loc='upper right')\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-31T16:50:23.716863Z","iopub.execute_input":"2022-01-31T16:50:23.717188Z","iopub.status.idle":"2022-01-31T16:50:25.755536Z","shell.execute_reply.started":"2022-01-31T16:50:23.717159Z","shell.execute_reply":"2022-01-31T16:50:25.754627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# End of Part 1","metadata":{}},{"cell_type":"markdown","source":"# Start of Part 2","metadata":{}},{"cell_type":"markdown","source":"# Prediction\n\nWe will now proceed to predict our data. We will use our own test set, obtained from the original train set of the submission, in order to obtain the metrics we want (ROC_AUC, F1)","metadata":{}},{"cell_type":"code","source":"predictions = {}\nstopper = 3600 * 4 #8 hours limit for prediction\npred_start_time = time.time()\\\n\np, c = time.time(), time.time()\nbatch_size = 1000\n    \nl = 0\nn = test.shape[0]\n\nfor x in custom_dcom_image_generator(batch_size, test, True, False):\n    clear_output(wait=True)\n    model = load_model(\"../input/models/embed_2048_lr0.00001_norm.h5\")\n#     model = trained_model\n    preds = model.predict(x, batch_size=8, verbose=1)\n    \n    try:\n        for key in preds.keys():\n            predictions[key] += preds[key].flatten().tolist()\n            \n    except Exception as e:\n        print(e)\n        for key in preds.keys():\n            predictions[key] = preds[key].flatten().tolist()\n            \n    l = (l+batch_size)%n\n    print('Total predicted:', len(predictions['indeterminate']),'/', n)\n    p, c = c, time.time()\n    print(\"One batch time: %.2f seconds\" %(c-p))\n    print(\"ETA: %.2f\" %((n-l)*(c-p)/batch_size))\n    \n    if c - pred_start_time >= stopper:\n        print(\"Time's up!\")\n        break\n    \n    del model\n    K.clear_session()\n    \n    del x, preds\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-30T20:37:48.236372Z","iopub.execute_input":"2022-01-30T20:37:48.236785Z","iopub.status.idle":"2022-01-31T00:20:36.463345Z","shell.execute_reply.started":"2022-01-30T20:37:48.23675Z","shell.execute_reply":"2022-01-31T00:20:36.462492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now check the shape of each key in predictions to make sure we predicted everything.","metadata":{}},{"cell_type":"code","source":"pred_file = \"pred_file.pkl\"\nwith open(pred_file,'wb') as f:\n    pickle.dump(predictions, f)\nFileLink(pred_file)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T00:20:36.467099Z","iopub.execute_input":"2022-01-31T00:20:36.467381Z","iopub.status.idle":"2022-01-31T00:20:36.649234Z","shell.execute_reply.started":"2022-01-31T00:20:36.467352Z","shell.execute_reply":"2022-01-31T00:20:36.648403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Import predictions if necessary","metadata":{}},{"cell_type":"code","source":"import pickle \nwith open('../input/predictions/pred_embed_2048_lr0.00001_norm.pkl','rb') as f:\n    predictions=pickle.load(f)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.153222Z","iopub.execute_input":"2022-01-31T19:10:33.153586Z","iopub.status.idle":"2022-01-31T19:10:33.811392Z","shell.execute_reply.started":"2022-01-31T19:10:33.153554Z","shell.execute_reply":"2022-01-31T19:10:33.810596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for key in predictions.keys():\n    print(key, np.array(predictions[key]).shape)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.813106Z","iopub.execute_input":"2022-01-31T19:10:33.81346Z","iopub.status.idle":"2022-01-31T19:10:33.864898Z","shell.execute_reply.started":"2022-01-31T19:10:33.813432Z","shell.execute_reply":"2022-01-31T19:10:33.863139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we will convert predictions into a dataframe based on the test dataframe. We will also copy all unique **StudyInstanceUID** for predicting later.","metadata":{}},{"cell_type":"code","source":"for key in y_test.keys():\n    print(key, np.array(y_test[key]).shape)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.865799Z","iopub.status.idle":"2022-01-31T19:10:33.866189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate the image-level test performance ","metadata":{}},{"cell_type":"code","source":"def GetFinalScores(pred, y_test, agg=False):\n    test_scores={}\n    pred=round(pred)\n    if not agg:\n        keys=['pe_present_on_image', 'leftsided_pe','chronic_pe', 'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate']\n    else:\n        keys=['negative_exam_for_pe', 'leftsided_pe','chronic_pe', 'rightsided_pe', 'acute_and_chronic_pe', 'central_pe', 'indeterminate']\n    for key in keys:\n        fpr, tpr, thresholds = roc_curve(y_test[key], pred[key])\n        roc_auc = auc(fpr, tpr)\n        conf_matrix=confusion_matrix(y_test[key], pred[key])\n\n#         precision = precision_score(y_test[key], pred[key])\n\n        recall = recall_score(y_test[key], pred[key])\n\n        accuracy = balanced_accuracy_score(y_test[key], pred[key])\n        \n        test_scores[key]=[roc_auc,conf_matrix,recall,accuracy]\n\n    return test_scores","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.867482Z","iopub.status.idle":"2022-01-31T19:10:33.868141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Prediction size may be smaller than y_test if 4 hours weren't enough to predict entire test set\npred_df=pd.DataFrame.from_dict(predictions)\n\nif len(y_test)!= len(pred_df):\n    y_test_cropped=y_test[:len(pred_df)]\n    test_scores=GetFinalScores(pred_df,y_test_cropped)\nelse:\n    test_scores=GetFinalScores(pred_df,y_test)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.869482Z","iopub.status.idle":"2022-01-31T19:10:33.87012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_scores","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.871444Z","iopub.status.idle":"2022-01-31T19:10:33.872073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_ids = []\nfor v in test.StudyInstanceUID:\n    if v not in test_ids:\n        test_ids.append(v)\n\ntest_preds =pd.DataFrame(predictions)\ntest_preds.insert(0, \"StudyInstanceUID\", test.StudyInstanceUID.values)\ntest_preds","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.873406Z","iopub.status.idle":"2022-01-31T19:10:33.874041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Aggregate image-level predictions to exam-level","metadata":{}},{"cell_type":"markdown","source":"Now, we have image-level predictions.\nEach Series contains several images, and each Study/Exam contains several Series of Images.\nWe need the Exam level labels.\n\nFor that, we need to aggregate the image-level labels.\n\n- First task, define the label *negative_exam_for_pe*, which defines if the total exam is positive or negative for PE. To do this:\n    -  Check for any image level presence of PE of high confidence (threshold of 0.5). If any image has pe_presence>0.5, its positive, otherwise its negative.\n    -  If positive, set the study/exam level label of indeterminate to half the minimum probability of indeterminate image-level label and the *negative_exam_for_pe to 0*.\n    -  If negative, the study/exam level label can be indeterminate or negative. If any image level label of indeterminate is above 0.5 (confidently indeterminate), set the study/exam level label of indeterminate to the max probability of indeterminate image-level label and the *negative_exam_pe* to 1. Otherwise the exam is negative, the indeterminate study/exam level label is set to half the minimum and the *negative_exam_pe* is set to 1.\n    \n    \n- Second task, define the exam level label of severity: *chronic or acute_chronic_pe* and if none of these labels are confident, set to *acute*.\n    -  Check for any image level labels of acute_chronic and chronic of high confidence (threshold of 0.5).\n    -  If chronic, to agggregate the image level labels, set the exam-level label of chronic to 0.5+mean(image-level chronic label)/2 and the examl-level label of acute_chronic to mean(image-level acute_chronic label)/2. Similar is done if any image level label positive for acute_chronic with confidence (above 0.5)\n    - If none, we set the exam level label to acute, by setting the exam level label of chronic and acute_chronic to mean(image-level)/2\n\n- Third task, define the exam level label of position\n    -  Same process as before but with the labels, right, left and center.","metadata":{}},{"cell_type":"code","source":"from scipy.special import softmax\n\nlabel_agg = {key:[] for key in \n             ['id', 'negative_exam_for_pe',\n              'leftsided_pe', 'chronic_pe',\n              'rightsided_pe', 'acute_and_chronic_pe',\n              'central_pe', 'indeterminate']\n            }\n\n# label_agg = {key:[] for key in \n#              ['id', 'negative_exam_for_pe', 'rv_lv_ratio_gte_1',\n#               'rv_lv_ratio_lt_1', 'leftsided_pe', 'chronic_pe',\n#               'rightsided_pe', 'acute_and_chronic_pe',\n#               'central_pe', 'indeterminate']\n#             }\n\nfor uid in test_ids:\n    #to the label aggregation dataframe, we also add the studyinstance to which we refer in the field id\n    temp = test_preds.loc[test_preds.StudyInstanceUID == uid]\n    label_agg['id'].append(uid)\n    \n    n = temp.shape[0]\n    #Check for any image level presence of PE of high confidence\n    positive = any(temp.pe_present_on_image >= 0.5) #50% threshhold\n    \n    #Only one from positive, negative and indeterminate should have value>0.5\n    #per exam\n    if positive: \n        #if positive, set the study/exam level label of indeterminate to half the minimum \n        #and the negative_exam_for_pe to 0 (n)\n        label_agg['indeterminate'].append(temp.indeterminate.min()/2) \n        label_agg['negative_exam_for_pe'].append(0)\n    else:\n        #if negative i.e., no image in the study has more than 50% chance of presence of PE\n        if any(temp.indeterminate >= 0.5):\n            #it can be indeterminate, if any of the images has more than 50% chance of indeterminate\n            label_agg['indeterminate'].append(temp.indeterminate.max())\n            label_agg['negative_exam_for_pe'].append(1)\n        else:\n            #or it can be negative otherwise\n            label_agg['indeterminate'].append(temp.indeterminate.min()/2)\n            label_agg['negative_exam_for_pe'].append(1)\n    \n    #I decided that the total ratio should be equal to 1, so I used softmax\n    #We modify the weights by multiplying the bigger by 2 and dividing the smaller by 2\n#     a, b = temp[['rv_lv_ratio_gte_1', 'rv_lv_ratio_lt_1']].mean().values\n#     if a > b:\n#         a, b = a*2, b/2\n#     elif a < b:\n#         a, b = a/2, b*2\n#     a, b = softmax([a, b])\n#     if positive:\n#         label_agg['rv_lv_ratio_gte_1'].append(a)\n#         label_agg['rv_lv_ratio_lt_1'].append(b)\n#     else:\n#         label_agg['rv_lv_ratio_gte_1'].append(a/2)\n#         label_agg['rv_lv_ratio_lt_1'].append(b/2)\n    \n    #Next is for Chronic (C), Acute-Chronic (AC) and Acute (A) PE\n    #We need to see if we got a high confidence value from either C or AC\n    #If there is, we add it to a 50% based score for high confidence\n    #and half weight for low confidence score\n    if any(temp['acute_and_chronic_pe'] > 0.5): #50% confidence level\n        label_agg['acute_and_chronic_pe'].append(0.5 + temp['acute_and_chronic_pe'].mean()/2)\n        label_agg['chronic_pe'].append(temp['chronic_pe'].mean()/2)\n        \n    elif any(temp['chronic_pe'] > 0.5):\n        label_agg['acute_and_chronic_pe'].append(temp['acute_and_chronic_pe'].mean()/2)\n        label_agg['chronic_pe'].append(0.5 + temp['chronic_pe'].mean()/2)\n        \n    else: #Else, we set both to half values, as we declare the A as the value\n        label_agg['acute_and_chronic_pe'].append(temp['acute_and_chronic_pe'].mean()/2)\n        label_agg['chronic_pe'].append(temp['chronic_pe'].mean()/2)\n    \n    #for right, left, central, we use the same metric above\n    for key in ['leftsided_pe', 'rightsided_pe', 'central_pe']:\n        if positive:\n            label_agg[key].append(0.5 + temp[key].mean()/2)\n        else:\n            label_agg[key].append(temp[key].mean()/2)","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.875373Z","iopub.status.idle":"2022-01-31T19:10:33.875965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate the exam-level test scores","metadata":{}},{"cell_type":"markdown","source":"To do so, we assume that the true exam-level labels are the same as the image-level labels for the test set. If we take a look at the data, if one image of the exam has acute_chronic_pe=1, all the others have it too.","metadata":{}},{"cell_type":"code","source":"label_agg=pd.DataFrame(label_agg)\n\ny_test_agg=pd.DataFrame(y_test)\ny_test_agg.insert(0,'StudyInstanceUID',test.StudyInstanceUID.values)\ny_test_agg=y_test_agg.drop_duplicates(subset=\"StudyInstanceUID\")\ny_test_agg=y_test_agg.drop(columns=['pe_present_on_image','qa_motion',\n                                    'qa_contrast', 'flow_artifact', 'rv_lv_ratio_gte_1', \n                                    'rv_lv_ratio_lt_1','true_filling_defect_not_pe'])\n# for uid in y_test.StudyInstanceUID:\n#     if uid not in y_test_agg:\n#         y_test_agg.append(y_test.loc(\"S\"))\n\ntest_scores_exam=GetFinalScores(label_agg,y_test_agg,True)\n","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.877242Z","iopub.status.idle":"2022-01-31T19:10:33.877871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_scores_exam","metadata":{"execution":{"iopub.status.busy":"2022-01-31T19:10:33.879198Z","iopub.status.idle":"2022-01-31T19:10:33.879844Z"},"trusted":true},"execution_count":null,"outputs":[]}]}