{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"# Install all necessary libraries\n!pip3 install pillow pydicom tensorflow #python-varname","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Import all necessary libraries\n\nimport csv\nimport numpy as np\nimport pandas as pd\nimport os\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport torch\nimport torch.optim as optim\nfrom torch import nn\nimport sys\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset as TorchDataset\nfrom torch.utils.data import DataLoader\nfrom torchvision import datasets, transforms, models\nfrom PIL import Image\nimport PIL.ImageOps\nimport math\nimport sys\nimport random\nimport collections\nimport requests\nimport glob\nimport pydicom as pdm\n#import varname\nimport tensorflow as tf\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense, Dropout, Activation, Flatten\nfrom tensorflow.keras.layers import Conv2D, MaxPooling2D\nfrom keras.utils import to_categorical\nfrom keras.callbacks import ModelCheckpoint\nfrom tqdm.autonotebook import tqdm\nfrom sklearn.utils import shuffle\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.metrics import accuracy_score\nfrom sklearn.metrics import classification_report\nimport cv2\nimport gc\nfrom tensorflow.python.keras.utils.data_utils import Sequence\nfrom tensorflow.python.ops import array_ops","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.seterr(divide='ignore', invalid='ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Define Constants"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Input parameters\nINPUT_DIR = \"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/\"\n\n\n# Preprocessing Parameters\nHU_RESCALE = True\nWINDOW = True\nNORMALIZE = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Model Structure Parameters\nBATCH_SIZE = 16  # can increase with greater gpu strength\nVAL_BATCH_RATIO = .2  # ratio of (training entries):(validation entries)\ncriterion = nn.CrossEntropyLoss()\n#optimizer = torch.optim.Adam()\n#Activation Function = F.relu(): relu usually better than sigmoid, esp for CNNs\nINITIAL_LR = 4e-5  # initial learning rate may decay over time - this is max LR\nMODEL_PATH = '/kaggle/input/baseline-resnext50/resnext50_10.pth'\n\nPYTORCH_EPOCHS = 2\nKERAS_EPOCHS = 10\n\n\n# Model Analysis Parameters\nIMAGE_MODE = '2d'\nMANUAL_ANY = False\nHEMORRHAGE_TYPES = ('epidural', 'intraparenchymal', 'intraventricular', 'subarachnoid', 'subdural')\n\nN_DATA = 5000  # full data = 752803 entries\nVALID_SIZE = math.ceil(N_DATA*VAL_BATCH_RATIO) #validation data size\nN_THREADS = 4  # num parallell ops\nSHUFFLE = True\n\nSUBMISSION_SIZE = 100\n\nKERAS = True\n\n#Keras Image Ranges\nKERAS_VALID_IMAGE_START = 1000\nKERAS_VALID_IMAGE_END = 1200\nKERAS_TRAIN_IMAGE_START = 0\nKERAS_TRAIN_IMAGE_END = 1000","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Layer Parameters\nIN_CHANNELS = 1  # B&W = 1, Greyscale = 2, RGB = 3\n\nCONV1_CHANNELS = 512  # depth of conv layer 1 – one channel per pixel\nCONV1_KERN = 5  # kernal length and width\nCONV1_STRIDE = 1  # kernal must be able to move this much and cover all channels\nCONV1_PAD = 0  # padding for conv1\n\nCONV2_CHANNELS = 50\nCONV2_KERN = 3\nCONV2_STRIDE = 1\nCONV2_PAD = 0\n\nFC1_NEURONS = 500\n\nPOOL1_X = 2\nPOOL1_Y = 2\n\nPOOL2_X = 2\nPOOL2_Y = 2\n\nDR1_PROB = .2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Set device for training the model (If gpu available - GPU, otherwise - CPU)\n\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pandas Helper Class"},{"metadata":{"trusted":true},"cell_type":"code","source":"class PdFuncs():\n    \n    @staticmethod\n    def row_at_loc(df, val, search_col):\n        return df.loc[df[search_col] == val]\n\n    @staticmethod\n    def extract_row_value(row, idx):\n        try:\n            return row.values.flatten().tolist()[idx]\n        except AttributeError:\n            return row.flatten().tolist()[idx]\n\n    @staticmethod\n    def extract_col_value():\n        pass\n\n    @staticmethod\n    def extract_value(df, row, col):\n        try:\n            return df.iat(row, col)\n        except TypeError:\n            return df.at(row, df.columns[col])\n\n\n    @staticmethod\n    def match_extract(df, val, search_col, extract_col):\n        row = PdFuncs.row_at_loc(df, val, search_col)\n        try:\n            return PdFuncs.extract_row_value(row, extract_col)\n        except TypeError:\n            return PdFuncs.extract_row_value(row, df.columns.get_loc(extract_col))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Formula Calculations"},{"metadata":{"trusted":true},"cell_type":"code","source":"def output_size(inp, kern, stride, pad):  # verify dimensions, this is 3d (ret, ret, filters)\n    return ((inp - kern + 2*pad)//stride) + 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_any(): pass","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Processing"},{"metadata":{"trusted":true},"cell_type":"code","source":"# List and define paths constants\n\nTRAIN = INPUT_DIR + \"stage_2_train/\"\nTRAIN_CSV = INPUT_DIR + \"stage_2_train.csv\"\nTRAIN_EDIT_CSV = \"/kaggle/working/test.csv\"\nTEST = INPUT_DIR + \"stage_2_test/\"\nTEST_CSV = INPUT_DIR + \"stage_2_sample_submission.csv\"\nTEST_DIR = INPUT_DIR\nos.listdir(INPUT_DIR)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Delete corrupted image IDs\n\ncorrupted_images = ('ID_6431af929_epidural', 'ID_6431af929_intraparenchymal', 'ID_6431af929_intraventricular', 'ID_6431af929_subarachnoid', 'ID_6431af929_subdural', 'ID_6431af929_any')\n\nwith open(TRAIN_CSV, 'r') as inp, open(TRAIN_EDIT_CSV, 'w') as out:\n    writer = csv.writer(out)\n    for row in csv.reader(inp):\n        if row[0] not in corrupted_images:\n            writer.writerow(row)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess(dcm):\n    m = dcm.RescaleSlope\n    b = dcm.RescaleIntercept\n    c = dcm.WindowCenter\n    w = dcm.WindowWidth\n    \n    dcm = dcm.pixel_array\n    \n    # Convert to HU Units\n    if HU_RESCALE:\n        dcm = m*dcm + b\n    \n    # Window DCM\n    if WINDOW:\n        if isinstance(c, pdm.multival.MultiValue):\n            c = c[0]\n        if isinstance(w, pdm.multival.MultiValue):\n            w = w[0]\n    bounds = (c - w//2, c + w//2)\n    dcm = np.clip(dcm, bounds[0], bounds[1])\n    dcm = np.expand_dims(dcm, axis=2)\n#     dcm.resize((512, 512, 2))\n#     print(dcm.shape)\n    # Normalize Pixel Array\n    if NORMALIZE:\n        dcm = 255*(dcm - dcm.min()) / (dcm.max() - dcm.min())\n        \n    return dcm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_dcm(image_id, directory=TRAIN, mode=\"image\"):\n    dcm = pdm.dcmread(directory+f'ID_{image_id}.dcm')\n  \n    if mode == \"image\":\n        return preprocess(dcm)\n    if mode == \"patient\":\n        return dcm.PatientID\n    if mode == \"angle\":\n        return dcm.ImagePositionPatient  #https://stackoverflow.com/questions/30814720/dicom-and-the-image-position-patient\n    if mode == \"orientation\":\n        return dcm.ImageOrientationPatient","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_dcm(arr):\n    arr = np.squeeze(arr, axis=(2, ))\n    #print(arr.shape)\n    #image = Image.fromarray(np.asarray(arr))\n    plt.axis('off')\n    plt.imshow(arr, cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm1 = read_dcm('8d3864098')\ndisplay_dcm(dcm1)\n# print(dcm1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm2 = read_dcm('d3e935321')\ndisplay_dcm(dcm2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm1_pat = read_dcm('8d3864098', mode=\"patient\")\ndcm1_pat","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read & Display Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"class CSV_Sampler():\n    def __init__(self, path, n=BATCH_SIZE, mode=IMAGE_MODE, incl_label=True):\n        self.mode = mode\n        self.unclassified_fnames = pd.DataFrame()\n        self.classified_fnames = pd.DataFrame()\n        \n        self.IMGpath = \"\"\n        if path == TRAIN_CSV: \n            self.IMGpath = TRAIN\n        elif path == TEST_CSV:\n            self.IMGpath = TEST\n\n        self.incl_label=True\n        self.unclassified_fnames = self.pull_data(path)\n        self.classified_fnames = self.categorize_data(self.unclassified_fnames)\n        \n        self.IDs = self.sample_ids(self.classified_fnames, n, mode)\n        self.IMGs = self.sample_images(self.IDs, mode)\n\n    def __call__(self):\n        pass\n\n    def __getitem__(self, key):\n        return torch.from_numpy(self.IMGs[key])\n\n    def __len__(self):\n        return len(self.IMGs)\n\n    def __str__(self):\n        return str(self.IMGs)\n    \n    @staticmethod\n    def classify_image_category(x):\n        return \"_\".join(x.split(\"_\", 2)[2:])\n\n    @staticmethod\n    def retrieve_image_id(x):\n        return \"_\".join(x.split(\"_\", 2)[1:2])\n\n    @staticmethod\n    def retrieve_patient_id(x):\n        return read_dcm(x, directory=self.IMGpath, mode=\"patient\")\n\n\n    def pull_data(self, path):\n        df = pd.read_csv(path)\n        df.rename(columns={'Label':'Category'}, inplace=True)\n        if self.IMGpath == TRAIN:\n            df['Category'] = np.where(df.Category == 1, df['ID'].apply(self.classify_image_category), 'X')\n        elif self.IMGpath == TEST:\n            df['Category'] = np.where(df.Category == 0.5, df['ID'].apply(self.classify_image_category), \"\")\n        df['ID'] = df['ID'].apply(self.retrieve_image_id)\n        # df.insert(2, \"Patient\", df['ID'].apply(cls.retrieve_patient_id))  # applies to all here but causes performance issues\n\n        return df\n\n    def categorize_data(self, data):\n        df = self.unclassified_fnames\n        fnames = pd.crosstab(df['ID'], df['Category']).astype(int).rename_axis(index=None,columns=None)\n        fnames.reset_index(inplace=True)\n        fnames.rename(columns={'index':'ID'}, inplace=True)\n        if self.IMGpath == TRAIN:\n            del fnames['X']\n        elif self.IMGpath == TEST:\n            fnames.drop_duplicates('ID')\n            fnames = fnames.replace(True, 'Unknown')  # broken\n\n        return fnames\n\n    def label_image(self, id):\n        df = self.classified_fnames\n        label = []\n\n        if MANUAL_ANY:\n            pass\n        else:\n            for i in range(6):\n                label.append(int(PdFuncs.match_extract(df, id, 'ID', i+1)))\n      \n\n        # any = df.loc[df['ID'] == id][1]\n        # epidural = df.loc[df['ID'] == id][2]\n        # intraparenchymal = df.loc[df['ID'] == id][3]\n        # intraventricular = df.loc[df['ID'] == id][4]\n        # subarachnoid = df.loc[df['ID'] == id,'subarachnoid']\n        # subdural = df.loc[df['ID'] == id,'subdural']\n\n        # # data = {\n        # #     'any': any,\n        # #     'epidural': epidural,\n        # #     'intraparenchymal': intraparenchymal,\n        # #     'intraventricular': intraventricular,\n        # #     'subarachnoid': subarachnoid,\n        # #     'subdural': subdural,\n        # # }\n\n        # return [any, epidural, intraparenchymal, \n        #         intraventricular, subarachnoid, subdural]\n\n        return label\n    \n\n    def sample_ids(self, df, n, mode):\n        if mode == '2d':  # 1 image per patient\n            ids = set()\n            pats = set()\n\n        for _ in range(n):\n            while True:\n                row = df.sample()\n                imgID = PdFuncs.extract_row_value(row, 0)\n                if imgID in ids:\n                    continue\n                elif read_dcm(imgID, directory=self.IMGpath, mode=\"patient\") in pats:\n                    continue\n                else:\n                    ids.add(imgID)\n                    pats.add(read_dcm(imgID, directory=self.IMGpath, mode=\"patient\"))\n                    break\n\n            return ids\n\n        if mode == '3d': \n            pass\n\n    def sample_images(self, ids, mode):\n        lab_imgs = []\n\n        if mode == '2d':\n            for id in ids:\n                if self.incl_label:\n                    img = read_dcm(id, directory=self.IMGpath)\n                    lab = self.label_image(id)\n                    lab_imgs.append([img, lab])\n#                 else:\n           \n\n            return lab_imgs\n\n        if mode == '3d':\n            pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data = CSV_Sampler(TRAIN_CSV)\nvalidation_data = CSV_Sampler(TRAIN_CSV, n=math.ceil(VAL_BATCH_RATIO*BATCH_SIZE))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data.classified_fnames","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data.classified_fnames.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Multilabel logarithmic loss function\n\nclass FocalLoss(nn.Module):\n    def __init__(self, alpha=1, gamma=2, logits=False, reduce=True):\n        super(FocalLoss, self).__init__()\n        self.alpha = alpha\n        self.gamma = gamma\n        self.logits = logits\n        self.reduce = reduce\n\n    def forward(self, inputs, targets):\n        if self.logits:\n            BCE_loss = F.binary_cross_entropy_with_logits(inputs, targets, reduce=False)\n        else:\n            BCE_loss = F.binary_cross_entropy(inputs, targets, reduce=False)\n        pt = torch.exp(-BCE_loss)\n        F_loss = self.alpha * (1-pt)**self.gamma * BCE_loss\n\n        if self.reduce:\n            return torch.mean(F_loss)\n        else:\n            return F_loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# EDA (Exploratory Data Analysis)"},{"metadata":{},"cell_type":"markdown","source":"### Label Frequency"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(TRAIN_CSV)\nsns.countplot(df_train.Label).set_title(\"Label Frequency\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Hemmorhage Frequency"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['Sub_type'] = df_train['ID'].str.split(\"_\", n = 3, expand = True)[2]\ndf_train['PatientID'] = df_train['ID'].str.split(\"_\", n = 3, expand = True)[1]\n#df_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gbSub = df_train.groupby('Sub_type').sum()\ngbSub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"freq = sns.barplot(y=gbSub.index, x=gbSub.Label, palette=\"deep\")\nfreq.set_title('Hemmorhage Frequency')\nfreq.set_xlabel('Frequency')\nfreq.set_ylabel('Subtype')\nfreq","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Subtype Frequency"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig=plt.figure(figsize=(10, 8))\n\nfrq = sns.countplot(x=\"Sub_type\", hue=\"Label\", data=df_train)\nfrq.set_xlabel('Subtype')\nfrq.set_ylabel('Frequency')\n\nplt.title(\"Subtype Frequency\")\n\nfrq","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_image(img, window_center,window_width, intercept, slope, rescale=True):\n\n    img = (img*slope +intercept)\n    img_min = window_center - window_width//2\n    img_max = window_center + window_width//2\n    img[img<img_min] = img_min\n    img[img>img_max] = img_max\n    \n    if rescale:\n        # Extra rescaling to 0-1, not in the original notebook\n        img = (img - img_min) / (img_max - img_min)\n    \n    return img\n    \ndef first_in_field(x):\n    #get x[0] as in int is x is a 'pydicom.multival.MultiValue', otherwise get int(x)\n    if type(x) == pdm.multival.MultiValue:\n        return int(x[0])\n    else:\n        return int(x)\n\ndef get_windowing(data):\n    dicom_fields = [data[('0028','1050')].value, #window center\n                    data[('0028','1051')].value, #window width\n                    data[('0028','1052')].value, #intercept\n                    data[('0028','1053')].value] #slope\n    return [first_in_field(x) for x in dicom_fields]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def view_images(images, title = '', aug = None, path=TRAIN):\n    width = 5\n    height = 2\n    fig, axs = plt.subplots(height, width, figsize=(15,5))\n    \n    for im in range(0, height * width):\n        ''''\n        i = im // width\n        j = im % width\n        axs[i,j].imshow(image, cmap=plt.cm.bone) \n        axs[i,j].axis('off')'''''\n        \n        data = pdm.read_file(os.path.join(path,'ID_'+images[im]+ '.dcm'))\n        image = data.pixel_array\n        window_center , window_width, intercept, slope = get_windowing(data)\n        image_windowed = window_image(image, window_center, window_width, intercept, slope)\n\n\n        i = im // width\n        j = im % width\n        axs[i,j].imshow(image_windowed, cmap=plt.cm.bone) \n        axs[i,j].axis('off')\n        \n        \n    plt.suptitle(title)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'epidural') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Epidural Hemmorhages')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'intraparenchymal') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Intraparenchymal Hemmorhages')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'intraventricular') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Intraventricular Hemmorhages')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'subarachnoid') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Subarachnoid Hemmorhages')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'subdural') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Subdural Hemmorhages')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images(df_train[(df_train['Sub_type'] == 'any') & (df_train['Label'] == 1)][:10].PatientID.values, title = 'Hemmorhages of Any Type')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Pytorch implementation"},{"metadata":{},"cell_type":"markdown","source":"## Data Preparation"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train and validation dataframe preparation \n\ndf_train = pd.read_csv(TRAIN_CSV)\n# df_train = pd.read_csv(TRAIN_EDIT_CSV) # Use to train the whole data\ndf_train[['id', 'img', 'subtype']] = df_train['ID'].str.split('_', n=3, expand=True)\ndf_train['img'] = 'ID_' + df_train['img'] \n\ndf_train.drop_duplicates(inplace=True)\ndf_train = df_train.pivot(index='img', columns='subtype', values='Label').reset_index()\ndf_train['path'] = TRAIN + df_train['img'] + '.dcm'\n# print(len(df_train))\ndf_valid = df_train\ndf_train  = df_train.iloc[:N_DATA]  # train first N_DATA images\ndf_valid = df_valid.iloc[N_DATA:N_DATA + VALID_SIZE] # initialize the validation dataframe\ndf_valid.reset_index(inplace=True)\ndel df_valid['index']\n#df_train.rename(columns={'subtype':'idx'}, inplace=True)\ndf_train.head()\n# len(df_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Test dataframe preparation\n\ndf_test = pd.read_csv(TEST_CSV)\ndf_test[['id','img','subtype']] = df_test['ID'].str.split('_', expand=True)\ndf_test['img'] = 'ID_' + df_test['img']\ndf_test = df_test[['img', 'Label']]\ndf_test['path'] = TEST + df_test['img'] + '.dcm'\ndf_test.drop_duplicates(inplace=True)\n\ndf_test = df_test.reset_index(drop=True)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Dataset class"},{"metadata":{"trusted":true},"cell_type":"code","source":"class Dataset(TorchDataset):\n    \"\"\"Dataset preparation class. Contains a list of images in appropriate format\"\"\"\n    def __init__(self, df, labels):\n        self.data = df\n        self.labels = labels\n\n    def __len__(self):\n        \"\"\"The size of dataset\"\"\"\n        return len(self.data)\n\n    def __getitem__(self, index):\n        \"\"\"Get element by its index as a numpy array\"\"\"\n        \n        img_name = self.data.loc[index, 'path']   \n        \n        img_dcm = pdm.read_file(img_name)\n        img = Dataset.brain_window(img_dcm)\n        img = cv2.resize(img, (200,200))\n        \n        img = np.stack((img,)*3, axis=-1)  # add another dimension\n        img = np.transpose(img, (2, 1, 0))  # fix dimension sizes\n    \n                \n        if self.labels:        \n            labels = torch.tensor(\n                self.data.loc[index, ['epidural', 'intraparenchymal', 'intraventricular', 'subarachnoid', 'subdural', 'any']])\n            return {'image': img, 'labels': labels}   \n        else:\n            return {'image': img}\n  \n    @staticmethod      \n    def brain_window(img):\n        \"\"\"Process .dcm formatted images. All cells should are windowed between 0 and 80, then the image is normalized\"\"\"\n        window_min = 0\n        window_max = 80\n        _, _, intercept, slope = Dataset.get_windowing(img)\n        img = img.pixel_array.astype('float32')\n        img = img * slope + intercept\n        img = np.clip(img, window_min, window_max)  # Window\n        img = (img - np.min(img)) / 1e-5 + (np.max(img) - np.min(img))  # Normalize, 1e-5 so now div by 0\n        return img\n    \n    @staticmethod\n    def get_windowing(data):\n        \"\"\"Get data from .dcm image\"\"\"\n        dicom_fields = [data[('0028','1050')].value, #window center\n                        data[('0028','1051')].value, #window width\n                        data[('0028','1052')].value, #intercept\n                        data[('0028','1053')].value] #slope\n        return [Dataset.first_in_field(x) for x in dicom_fields]\n  \n    @staticmethod\n    def first_in_field(x):\n        if type(x) == pdm.multival.MultiValue:\n            return int(x[0])\n        else:\n            return int(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Set parameters for training\n\nparams = {'batch_size': BATCH_SIZE,\n          'shuffle': SHUFFLE,\n          'num_workers': N_THREADS}\n\n# Create Dataset based on our dataframe\n\ntrain_dataset = Dataset(df= df_train, labels=True)\ntest_dataset = Dataset(df= df_test, labels=False)\nvalid_dataset = Dataset(df = df_valid, labels=True)\n\n# Create DataLoader based on our dataset\n\ndata_train_generator = DataLoader(train_dataset, **params)\ndata_test_generator = DataLoader(test_dataset,**params)\ndata_valid_generator = DataLoader(valid_dataset, **params)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Example Images"},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_images(dataloader, num_images = 5):\n    batch = next(iter(dataloader))\n    fig, axs = plt.subplots(1, num_images, figsize=(15, 5))\n    \n    for i in np.arange(num_images):\n        axs[i].axis(\"off\")\n        axs[i].imshow(batch['image'][i][0].numpy(), cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Training Images"},{"metadata":{"trusted":true},"cell_type":"code","source":"display_images(data_train_generator)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test Images"},{"metadata":{"trusted":true},"cell_type":"code","source":"display_images(data_test_generator)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Set resnet pretrained model, add a linear layer as the last layer because we have 6 classes, not 1000\n\nmodel0 = models.resnext50_32x4d(pretrained=True)  # 1d batch size x 3d images\nmodel = torch.nn.Sequential(model0, torch.nn.Linear(1000, 6) )  # add linear layer for output of 6, 1000 is default resnet output\n\n# Set device (GPU or CPU) and loss function\n\nmodel = model.to(device)\ncriterion = torch.nn.BCEWithLogitsLoss()  # logarithmic loss BCE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def model_summary_pytorch(model):\n    layers = [module for module in model.modules() if type(module) != nn.Sequential]\n    print(layers)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_summary_pytorch(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Set number of epochs and define optimizer\n\nn_epochs = PYTORCH_EPOCHS\noptimizer = optim.Adam(model.parameters(), lr=INITIAL_LR)  # maybe implement changing LR\nepochs = list(range(1, n_epochs+1))\n\ntry:\n    model.load_state_dict(torch.load(MODEL_PATH))\n    torch.save(model.state_dict(), 'resnext50_0.pth') \nexcept Exception:  # figure out exception\n    print('Resnet loaded')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train Model"},{"metadata":{},"cell_type":"markdown","source":"### Performance Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_model_report(actual, predicted, mode):\n    print(mode + ' mode report :\\n')\n    actual = actual.cpu().numpy()\n    predicted = predicted.cpu().detach().numpy()\n    #print(actual)\n    #print(predicted)\n#     actual = np.add(np.argmax(actual, axis = 1), 1).tolist()\n    actual = np.argmax(actual, axis = 1).tolist()\n    predicted = np.argmax(predicted, axis = 1).tolist()\n    matrix = confusion_matrix(actual, predicted)\n    #report = classification_report(actual, predicted)\n    print('Confusion Matrix :')\n    print(matrix)\n    #print('Report :')\n    #print(report)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Training"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Standard pytorch training procedure\n\ntraining_losses = []\nval_losses = []\nfor epoch in range(1, n_epochs+1):\n    \n    print(f'Epoch {epoch}/{n_epochs}')\n    print('-' * 10)\n\n    model.train()    \n    tr_loss = 0\n    \n    tk_train = tqdm(data_train_generator, desc='Images Processed')  # data enumerator\n    \n    for step, batch in enumerate(tk_train):\n        \n        inputs = batch[\"image\"]\n        labels = batch[\"labels\"]\n\n        inputs = inputs.to(device, dtype=torch.float)\n        labels = labels.to(device, dtype=torch.float)\n\n        outputs = model(inputs)\n        loss = criterion(outputs, labels)\n                \n        loss.backward()\n\n        optimizer.step()\n        optimizer.zero_grad()\n     \n        tr_loss += loss.item()\n        \n        # Save scores for confusion matrix\n            \n#     torch.save(model.state_dict(), f'resnext50_{epoch}.pth') \n    if(epoch == n_epochs):\n        show_model_report(labels, outputs, mode='training')\n    \n    model.eval()\n    val_loss = 0\n    tk_valid = tqdm(data_valid_generator, desc='Images Processed')\n    \n    for step, batch in enumerate(tk_valid):\n        inputs = batch[\"image\"]\n        labels = batch[\"labels\"]\n\n        inputs = inputs.to(device, dtype=torch.float)\n        labels = labels.to(device, dtype=torch.float)\n\n        outputs = model(inputs)\n        loss = criterion(outputs, labels)\n        \n        optimizer.zero_grad()\n        \n        val_loss += loss.item()\n    \n    if(epoch == n_epochs):\n        show_model_report(labels, outputs, mode='validation')\n    \n    epoch_train_loss = tr_loss / len(data_train_generator)\n    epoch_val_loss = val_loss / len(data_valid_generator)\n    training_losses.append(epoch_train_loss)\n    val_losses.append(epoch_val_loss)\n    print('Training Loss: {:.7f}'.format(epoch_train_loss))\n    print('Validation Loss: {:.7f}'.format(epoch_val_loss))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = sns.lineplot(epochs, training_losses)\nfig.set_xlabel('Epoch')\nfig.set_ylabel('Training Loss')\nfig.set_xticks(list(range(1, PYTORCH_EPOCHS+1)))\nfig","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = sns.lineplot(epochs, val_losses)\nfig.set_xlabel('Epoch')\nfig.set_ylabel('Validation Loss')\nfig.set_xticks(list(range(1, PYTORCH_EPOCHS+1)))\nfig","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"torch.save(model.state_dict(), f'resnext50_improved.pth')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Predict a random image"},{"metadata":{"trusted":true},"cell_type":"code","source":"def predict_image(image_id):\n    df = df_test.loc[df_test['img'] == 'ID_'+image_id]\n    print(df)\n    dataset = Dataset(df = df, labels=False)\n    dataloader = DataLoader(dataset, **params)\n    for param in model.parameters():\n        param.requires_grad = False\n        model.eval()\n\n    test_pred = np.zeros((len(dataset) * 6, 1))\n\n    tk_test = tqdm(dataloader)\n\n    for i, batch in enumerate(tk_test):\n        batch = batch[\"image\"]\n        print(batch.shape)\n        batch = batch.to(device, dtype=torch.float)\n        with torch.no_grad():\n            pred = model(batch)\n            test_pred[(i * 1 * 6):((i + 1) * 1 * 6)] = torch.sigmoid(pred).detach().cpu().reshape((len(batch) * 6, 1))\n    dictionary = {'ID': [image_id + '_epidural', image_id + '_intraparenchymal', image_id + '_intraventricular', image_id + '_subarachnoid', image_id + '_subdural', image_id + '_any'], 'Label': [0, 0, 0, 0, 0, 0]}\n    submission = pd.DataFrame(dictionary) \n    submis = pd.concat([submission.drop(columns=['Label']), pd.DataFrame(test_pred)], axis=1)\n    submis.columns = ['ID', 'Label']\n    #print(submis)\n    return submis","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predict_image('0fbf6a978')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Predict Test Set"},{"metadata":{"trusted":true},"cell_type":"code","source":"for param in model.parameters():\n    param.requires_grad = False\n    \nmodel.eval()\n\ntest_pred = np.zeros((len(test_dataset) * 6, 1))\n\ntk_test = tqdm(data_test_generator)\n\nfor i, batch in enumerate(tk_test):\n    if i == SUBMISSION_SIZE//BATCH_SIZE:\n        break\n    \n    batch = batch[\"image\"]\n    batch = batch.to(device, dtype=torch.float)\n    #print(batch.shape)\n\n    with torch.no_grad():\n        \n        pred = model(batch)\n        #print(pred.shape)\n        test_pred[(i * BATCH_SIZE * 6):((i + 1) * BATCH_SIZE * 6)] = torch.sigmoid(\n            pred).detach().cpu().reshape((len(batch) * 6, 1))\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission =  pd.read_csv(TEST_CSV)\nsubmission = pd.concat([submission.drop(columns=['Label']), pd.DataFrame(test_pred)], axis=1)\nsubmission.columns = ['ID', 'Label']\n\nsubmission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SUBMIT_CSV = '/kaggle/working/submission.csv'\ndf_submit = pd.read_csv(SUBMIT_CSV)\ndf_submit[:50]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Keras implementation"},{"metadata":{"trusted":true},"cell_type":"code","source":"if not KERAS: raise Exception('Keras implementation not selected')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(TRAIN_CSV)\ntrain.shape\nread_data = train.copy()\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['filename'] = train['ID'].apply(lambda x: \"ID_\" + x.split('_')[1] + \".dcm\")\ntrain['type'] = train['ID'].apply(lambda x: x.split('_')[2])\ntrain.head()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = train[['Label', 'filename', 'type']].drop_duplicates().pivot(\n    index='filename', columns='type', values='Label').reset_index()\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = shuffle(train)\ntrain_sample=train\ntrain_sample2=train_sample.reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_sample2.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train_sample2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"yvals=pd.DataFrame(train_sample2,columns=['any','epidural','intraparenchymal','intraventricular','subarachnoid','subdural'])\nyvals.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xhead=pd.DataFrame(train_sample2,columns=['filename'])\nxhead.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def first_in_field(x):\n    #get x[0] as in int is x is a 'pydicom.multival.MultiValue', otherwise get int(x)\n    if type(x) == pdm.multival.MultiValue:\n        return int(x[0])\n    else:\n        return int(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_windowing(data):\n    dicom_fields = [data[('0028','1050')].value, #window center\n                    data[('0028','1051')].value, #window width\n                    data[('0028','1052')].value, #intercept\n                    data[('0028','1053')].value] #slope\n    return [first_in_field(x) for x in dicom_fields]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class DataGenerator(Sequence):\n    'Generates data for Keras'\n    def __init__(self, list_IDs_labels, batch_size=100, dim=(512,512)):\n        'Initialization'\n        self.dim = dim\n        self.batch_size = batch_size\n        self.list_IDs = list_IDs_labels\n        self.on_epoch_end()\n\n    def __len__(self):\n        'Denotes the number of batches per epoch'\n        return int(np.floor(len(self.list_IDs) / self.batch_size))\n\n    def __getitem__(self, index):\n        'Generate one batch of data'\n        # Generate indices of the batch\n        indexes = self.indexes[index*self.batch_size:(index+1)*self.batch_size]\n\n        # Find list of IDs\n        list_IDs_temp = [self.list_IDs['filename'][k] for k in indexes]\n        list_label_temp=[[int(self.list_IDs['any'][i]),int(self.list_IDs['epidural'][i]),int(self.list_IDs['intraparenchymal'][i]),int(self.list_IDs['intraventricular'][i]),int(self.list_IDs['subarachnoid'][i]),int(self.list_IDs['subdural'][i])] for i in indexes]\n        # Generate data\n        X, y = self.__data_generation(list_IDs_temp,list_label_temp)\n        X = tf.cast(X, tf.float32)\n        return X, y\n\n    def on_epoch_end(self):\n        'Updates indexes after each epoch'\n        self.indexes = np.arange(len(self.list_IDs))\n\n    def __data_generation(self, list_IDs_temp,list_label_temp):\n        'Generates data containing batch_size samples' # X : (n_samples, *dim, n_channels)\n        # Initialization\n        X = []\n        y = []\n\n        # Generate data\n        for i, ID in enumerate(list_IDs_temp):\n            # Store sample\n            ds=pdm.dcmread(TRAIN +list_IDs_temp[i] )\n            temp=ds.pixel_array\n            window_center , window_width, intercept, slope = get_windowing(ds)\n            img = window_image(temp, 50, 100, intercept, slope)\n            resized = cv2.resize(img, (200, 200))\n            X.append(resized)       \n        X=np.array(X).reshape(-1,200,200,1)\n        y_train=np.asarray(list_label_temp) \n        return X,y_train\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"valid=train_sample2[KERAS_VALID_IMAGE_START:KERAS_VALID_IMAGE_END]\nvalid=valid.reset_index(drop=True)\nvalid.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"traingen=DataGenerator(train_sample2[KERAS_TRAIN_IMAGE_START:KERAS_TRAIN_IMAGE_END])\nvalidgen=DataGenerator(valid)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def focal_loss(prediction_tensor, target_tensor, weights=None, alpha=0.25, gamma=2):\n    r\"\"\"Compute focal loss for predictions.\n        Multi-labels Focal loss formula:\n            FL = -alpha * (z-p)^gamma * log(p) -(1-alpha) * p^gamma * log(1-p)\n                 ,which alpha = 0.25, gamma = 2, p = sigmoid(x), z = target_tensor.\n    Args:\n     prediction_tensor: A float tensor of shape [batch_size, num_anchors,\n        num_classes] representing the predicted logits for each class\n     target_tensor: A float tensor of shape [batch_size, num_anchors,\n        num_classes] representing one-hot encoded classification targets\n     weights: A float tensor of shape [batch_size, num_anchors]\n     alpha: A scalar tensor for focal loss alpha hyper-parameter\n     gamma: A scalar tensor for focal loss gamma hyper-parameter\n    Returns:\n        loss: A (scalar) tensor representing the value of the loss function\n    \"\"\"\n    prediction_tensor = tf.dtypes.cast(prediction_tensor, dtype=tf.float32)\n    sigmoid_p = tf.nn.sigmoid(prediction_tensor)\n    zeros = array_ops.zeros_like(sigmoid_p, dtype=sigmoid_p.dtype)\n    \n    # For poitive prediction, only need consider front part loss, back part is 0;\n    # target_tensor > zeros <=> z=1, so poitive coefficient = z - p.\n    pos_p_sub = array_ops.where(target_tensor > zeros, target_tensor - sigmoid_p, zeros)\n    \n    # For negative prediction, only need consider back part loss, front part is 0;\n    # target_tensor > zeros <=> z=1, so negative coefficient = 0.\n    neg_p_sub = array_ops.where(target_tensor > zeros, zeros, sigmoid_p)\n    per_entry_cross_ent = - alpha * (pos_p_sub ** gamma) * tf.math.log(tf.clip_by_value(sigmoid_p, 1e-8, 1.0)) \\\n                          - (1 - alpha) * (neg_p_sub ** gamma) * tf.math.log(tf.clip_by_value(1.0 - sigmoid_p, 1e-8, 1.0))\n    return tf.reduce_sum(per_entry_cross_ent)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = Sequential()\n\nmodel.add(Conv2D(32, (3, 3), input_shape=(200, 200,1)))\nmodel.add(Activation('relu'))\nmodel.add(MaxPooling2D(pool_size=(2, 2)))\nmodel.add(Dropout(0.2))\nmodel.add(Conv2D(32,(3,3)))\nmodel.add(Activation('relu'))\nmodel.add(MaxPooling2D(pool_size=(5,5)))\n# model.add(Dropout(0.25))\n\nmodel.add(Flatten())\n\n\n# model.add(Dense(100))\n# model.add(Activation('relu'))\n\nmodel.add(Dense(50))\nmodel.add(Activation('relu'))\n\nmodel.add(Dense(6))\nmodel.add(Activation('sigmoid'))\n\n# model.compile(loss='categorical_crossentropy',\n#               optimizer='adam',\n#               metrics=['accuracy'])\nmodel.compile(loss=focal_loss,\n               optimizer='adam',\n               metrics=['accuracy'])\n\n#model.fit(X,y_train,batch_size=32,epochs=6,validation_split=0.5)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Layer maps from training"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_filters_keras(model):\n    for layer in model.layers:\n        if 'conv' in layer.name:\n            weights, bias= layer.get_weights()\n            print(layer.name, weights.shape)\n        \n            #normalize filter values between  0 and 1 for visualization\n            f_min, f_max = weights.min(), weights.max()\n            filters = (weights - f_min) / (f_max - f_min)  \n            print(filters.shape[3])\n            filter_cnt=1\n        \n            #plotting all the filters\n            for i in range(filters.shape[3]):\n                #get the filters\n                filt=filters[:,:,:, i]\n                #plotting each of the channel, color image RGB channels\n                for j in range(filters.shape[2]):\n                    ax= plt.subplot(filters.shape[3], filters.shape[2], filter_cnt)\n                    ax.set_xticks([])\n                    ax.set_yticks([])\n                    plt.imshow(filt[:,:, j])\n                    filter_cnt+=1\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_filters_keras(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_grid_keras(datagenerator, image_index=0):\n    data = datagenerator[image_index][0]\n    successive_outputs = [layer.output for layer in model.layers[1:]]\n\n    visualization_model = tf.keras.models.Model(inputs = model.input, outputs=successive_outputs)\n    successive_feature_maps = visualization_model.predict(data)\n    layer_names = [layer.name for layer in model.layers]\n    for layer_name, feature_map in zip(layer_names, successive_feature_maps):\n        print(feature_map.shape)\n        if(len(feature_map.shape) == 4):\n            n_features = feature_map.shape[-1]\n            size = feature_map.shape[1]\n            display_grid = np.zeros((size, size * n_features))\n            for i in range(n_features):\n                data = feature_map[0, :, :, i]\n                data -= data.mean()\n                data /= data.std()\n                data *= 64\n                data += 128\n                data = np.clip(data, 0, 255).astype('uint8')\n                display_grid[:, i * size : (i + 1) * size] = data\n            scale = 20. / n_features\n            plt.figure(figsize=(scale*n_features, scale))\n            plt.title(layer_name)\n            plt.grid(False)\n            plt.imshow(display_grid, aspect='auto', cmap='bone')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_grid_keras(datagenerator=traingen)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history=model.fit_generator(generator=traingen,validation_data=validgen,use_multiprocessing=False,\n                    workers=-1,epochs=KERAS_EPOCHS)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}