{"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":"# Skipping Slices\nIn case of scans with very large number of slices, there isn't a lot of difference between the consecutive slices.\nFor eg. In a scan having 400 slices, there wouldn't be much difference between slice number 50,51 and 52. Therefore we can skip a few slices after every selected slice.\n\nHere I have calculated the number of slices to skip based on the shape of the scan and required stack size.\nUsing this approach along, I trained 4 different models on different MRI Sequences using different CV splits and then used ensembling.","metadata":{}},{"cell_type":"code","source":"!pip install monai","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport requests\nimport glob\nfrom tqdm import tqdm_notebook as tqdm\nimport random\nimport numpy as np\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.nn.functional as F\nfrom torchvision import transforms, utils\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\nfrom sklearn.metrics import classification_report, accuracy_score, roc_auc_score, confusion_matrix\nfrom sklearn.model_selection import StratifiedKFold\nimport albumentations as A\nfrom albumentations.pytorch import ToTensorV2\nimport gc\nfrom monai import transforms as T\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:41:34.158567Z","iopub.execute_input":"2021-10-08T04:41:34.159212Z","iopub.status.idle":"2021-10-08T04:41:38.437829Z","shell.execute_reply.started":"2021-10-08T04:41:34.159093Z","shell.execute_reply":"2021-10-08T04:41:38.437053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification'\ntrain_data = pd.read_csv(os.path.join(path, 'train_labels.csv'))\n\ntrain_data = train_data.dropna().reset_index(drop=True)\ntrain_data[\"kfold\"] = -1\n\ntrain_data = train_data.sample(frac=1, random_state=7).reset_index(drop=True)\n\nkf = StratifiedKFold(n_splits=4)\n\nfor fold, (trn_, val_) in enumerate(kf.split(X=train_data, y=train_data.MGMT_value.values)):\n    print(len(trn_), len(val_))\n    train_data.loc[val_, 'kfold'] = fold\n\nprint('Num of train samples:', len(train_data))\n# train_data.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:41:38.439423Z","iopub.execute_input":"2021-10-08T04:41:38.439689Z","iopub.status.idle":"2021-10-08T04:41:38.474525Z","shell.execute_reply.started":"2021-10-08T04:41:38.439653Z","shell.execute_reply":"2021-10-08T04:41:38.473821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_size = 256\nstack_size = 64\n\ndef dicom2array(path, voi_lut=True, fix_monochrome=True):\n    dicom = pydicom.read_file(path)\n    # VOI LUT (if available by DICOM device) is used to\n    # transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    data = cv2.resize(data, (img_size, img_size))\n    return data\n\ndef load_sequence(paths):\n    stack = []\n    \n    # load only non zero slices\n    for i, path in  enumerate(paths):\n        data = dicom2array(path)\n        if data.max() <= 0:\n            continue\n        else:\n            stack.append(data)\n    \n    # if all empty\n    if len(stack)==0:\n        return np.zeros((img_size,img_size,stack_size))\n    \n    stack = np.dstack(stack)\n    \n    # Skip slices(take every nth slice)\n    # Calculate value of n based on scan's shape and stack_size\n    n = stack.shape[2]//stack_size + 1 \n    # Select a random starting slice from first n slices\n    start = np.random.choice([i for i in range(n)]) \n    # Take every nth slice from the start slice\n    stack = stack[:,:,start::n] \n    \n    # If sequence is very small, repeat it multiple times. Better than leaving slices empty(ie. zero padding).\n    num_of_repetitions = stack_size//stack.shape[2]\n    stack = np.concatenate((stack,)*num_of_repetitions + (stack[:,:,:stack_size-stack.shape[2]*num_of_repetitions],), axis=2)\n    \n    return stack\n\ndef load_3d_dicom_images(scan_id, split = \"train\"):\n    \"\"\"\n    Sort the slices based on name \"correctly\". By default it comes like 1->10->100->101->102 (sorted by name in string format)\n    \"\"\"\n    # Flair\n    flair = sorted(glob.glob(f\"{path}/{split}/{scan_id}/FLAIR/*.dcm\"), key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n    flair_img = load_sequence(flair)\n    \n    # T1W\n#     t1w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1w/*.dcm\"), key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n#     t1w_img = load_sequence(t1w)\n    \n#     # T1WCE\n#     t1wce = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T1wCE/*.dcm\"), key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n#     t1wce_img = load_sequence(t1wce)\n    \n    \n#     # T2W\n#     t2w = sorted(glob.glob(f\"{path}/{split}/{scan_id}/T2w/*.dcm\"), key=lambda x: int(x.split('/')[-1].split('-')[-1].split('.')[0]))\n#     t2w_img = load_sequence(t2w)\n    return flair_img.astype(np.float32)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:09.762058Z","iopub.execute_input":"2021-10-08T04:42:09.762500Z","iopub.status.idle":"2021-10-08T04:42:09.775284Z","shell.execute_reply.started":"2021-10-08T04:42:09.762461Z","shell.execute_reply":"2021-10-08T04:42:09.774489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# let's write a simple pytorch dataloader\n\n\nclass BrainTumor(Dataset):\n    def __init__(self, df, path = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification', split = \"train\", validation_fold = 0):\n        \n        df.BraTS21ID = df.BraTS21ID.apply(lambda x: str(x).zfill(5))\n        self.labels = {}            \n        if split == \"val\":\n            self.split = 'train'\n            val_data = df[df.kfold==validation_fold]\n            brats = list(val_data[\"BraTS21ID\"])\n            mgmt = list(val_data[\"MGMT_value\"])\n            for b, m in zip(brats, mgmt):\n                self.labels[b] = m\n            \n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{self.split}/\" + \"/*\"))]\n            self.ids = [id for id in self.ids if id in val_data.BraTS21ID.values]\n        elif split == \"train\":\n            self.split = split\n            train_data = df[df.kfold!=validation_fold]\n            brats = list(train_data[\"BraTS21ID\"])\n            mgmt = list(train_data[\"MGMT_value\"])\n            for b, m in zip(brats, mgmt):\n                self.labels[b] = m\n            \n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{self.split}/\" + \"/*\"))]\n            self.ids = [id for id in self.ids if id in train_data.BraTS21ID.values]\n        else:\n            self.split = split\n            self.ids = [a.split(\"/\")[-1] for a in sorted(glob.glob(path + f\"/{self.split}/\" + \"/*\"))]\n            \n    \n    def __len__(self):\n        return len(self.ids)\n    \n    def preprocess(self, stack_dim):\n        return T.Compose(\n            [\n                T.AddChanneld(keys=['image'],), # Below steps won't work without adding channel\n                T.NormalizeIntensityd(keys=['image'],),\n#                 T.Spacingd(keys=['image'], pixdim=(1, 1, stack_dim)),\n#                 T.ResizeWithPadOrCropd(keys=['image'], mode='reflect', spatial_size=(img_size, img_size, stack_size)),\n            ]\n        )\n    \n    def get_transforms(self):\n        return A.Compose([\n                    A.OneOf([\n                        A.HorizontalFlip(),\n                        A.VerticalFlip(),\n                        A.IAAAffine(shear=(-15,15), mode='constant', cval=0),\n                        A.IAAAffine(scale=(0.9,1.2), mode='constant', cval=0),\n                        A.IAAAffine(translate_percent=(0,0.15), mode='constant', cval=0),\n                        A.IAAAffine(rotate=(-20,20), mode='constant', cval=0),\n                    ], p=0.7),\n                ])\n    \n    def __getitem__(self, idx):\n        imgs = load_3d_dicom_images(self.ids[idx], self.split)\n        \n        # Preprocess stack\n        stack_dim = imgs.shape[2]//stack_size*2 + 1\n        preprocess = self.preprocess(stack_dim)\n        # Extract out `3D` stacks to apply albumentation augmentations. \n        imgs = preprocess({'image':imgs})['image'][0,:,:,:] \n        \n        # Augment processed stack\n        transform = self.get_transforms()\n        imgs = transform(image=imgs)['image']\n        imgs = torch.unsqueeze(torch.from_numpy(imgs), 0) # Add channel dimension again.\n        \n        if self.split != \"test\":\n            label = self.labels[self.ids[idx]]\n            return imgs, torch.tensor([label], dtype = torch.float64)\n        else:\n            return imgs","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:10.578373Z","iopub.execute_input":"2021-10-08T04:42:10.578935Z","iopub.status.idle":"2021-10-08T04:42:10.598303Z","shell.execute_reply.started":"2021-10-08T04:42:10.578898Z","shell.execute_reply":"2021-10-08T04:42:10.597533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bs = 4\ntrain_dataset = BrainTumor(train_data, split='train', validation_fold=0)\ntrain_loader = DataLoader(train_dataset, batch_size=bs, shuffle=True, num_workers=4, pin_memory=True)\n\nval_dataset = BrainTumor(train_data, split='val', validation_fold=0)\nval_loader = DataLoader(val_dataset, batch_size=bs, shuffle=False, num_workers=4, pin_memory=True)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:11.258623Z","iopub.execute_input":"2021-10-08T04:42:11.259267Z","iopub.status.idle":"2021-10-08T04:42:11.310794Z","shell.execute_reply.started":"2021-10-08T04:42:11.259226Z","shell.execute_reply":"2021-10-08T04:42:11.310087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_imgs(imgs, cmap='gray'):\n    fig = plt.figure(figsize=(8,8))\n    for i in range(64):\n        img = imgs[:,:,i]\n        fig.add_subplot(8, 8, i+1)\n        plt.imshow(img, cmap=cmap)\n    plt.show()     \n\nfor i,(img, label) in enumerate(val_loader):\n    plot_imgs(img[0,0])\n    break","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:12.156932Z","iopub.execute_input":"2021-10-08T04:42:12.157697Z","iopub.status.idle":"2021-10-08T04:42:44.558080Z","shell.execute_reply.started":"2021-10-08T04:42:12.157652Z","shell.execute_reply":"2021-10-08T04:42:44.557273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Using DenseNet because the dataset size is small, dense connections can help in this case.\n\nsave_model_name = 'DenseNet-264-flair-fold-0.pt'\nfrom monai.networks.nets import DenseNet264\nmodel = DenseNet264(spatial_dims=3, in_channels=1, out_channels=1)\n# model.out = nn.Sequential(\n#                 nn.Dropout(p=0.5)\n#                 nn.Linear(1920, 1, bias=True)\n#             )\n\ncriterion = nn.BCEWithLogitsLoss()\n\nlr = 1e-4\nmin_lr = 3e-6\nn_epochs = 1 # Trial run\naccumulation_steps = 16\nearly_stop = 7\npatience = 2","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:49.258966Z","iopub.execute_input":"2021-10-08T04:42:49.259537Z","iopub.status.idle":"2021-10-08T04:42:49.826097Z","shell.execute_reply.started":"2021-10-08T04:42:49.259495Z","shell.execute_reply":"2021-10-08T04:42:49.825353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install torch_lr_finder\nfrom torch_lr_finder import LRFinder\noptimizer = torch.optim.AdamW(model.parameters(), lr=0.000001, weight_decay=0.001)\ngpu = torch.device(f\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nlr_finder = LRFinder(model, optimizer, criterion, device=gpu)\nlr_finder.range_test(train_loader, end_lr=0.001, num_iter=25)\nlr_finder.plot()\nlr_finder.reset()\n\n# Suggested LR: 2.74e-4","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Send updates on telegram. Really useful to keep track of training while commited.\ndef send_message(msg, chat_id, bot_token):\n    \"\"\"\n    params:\n    -------\n    msg: message you want to receive\n    chat_id: CHAT_ID\n    bot_token: API_KEY of your bot\n    \"\"\"\n\n    url  = f'https://api.telegram.org/bot{bot_token}/sendMessage'\n    data = {'chat_id': str(chat_id), 'text': f'{msg}'}\n    requests.post(url, data)\n\nchat_id = None\ntoken = ''\n# send_message('Training Started...', chat_id, token)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:53.776322Z","iopub.execute_input":"2021-10-08T04:42:53.777034Z","iopub.status.idle":"2021-10-08T04:42:53.782038Z","shell.execute_reply.started":"2021-10-08T04:42:53.776994Z","shell.execute_reply":"2021-10-08T04:42:53.781310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"optimizer = torch.optim.AdamW(model.parameters(),lr = lr, weight_decay=0.001, amsgrad=False)\nscheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, factor=1/3, patience=patience, min_lr=min_lr, verbose=True)\n\ngpu = torch.device(f\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nmodel.to(gpu)\n\nbest_loss = np.inf\n\nearly_stopping_counter = 0\nfor epoch in range(n_epochs):  # loop over the dataset multiple times\n\n    running_loss = 0.0\n    model.train()\n    train_bar = tqdm(enumerate(train_loader, 0), total=len(train_loader), desc = 'Training') \n    for i, data in train_bar: \n        x, y = data\n        \n        x = x.to(gpu)\n        y = y.to(gpu)\n\n        outputs = model(x)\n\n        loss = criterion(outputs, y)\n        running_loss += loss.item()\n        \n        loss = loss / accumulation_steps                \n        loss.backward()\n        \n        # Gradient Accumulation\n        if (i+1)%accumulation_steps == 0 or (i+1)==len(train_loader):             \n            optimizer.step()\n            optimizer.zero_grad()\n        train_bar.set_description(f\"Epoch: {epoch+1}, loss: {running_loss/(i+1)}\") \n\n    print(f\"epoch {epoch+1} train: {running_loss/len(train_loader)}\")\n#     send_message(f\"epoch {epoch+1} train: {running_loss/len(train_loader)}\", chat_id, token)\n\n    v_loss = 0.0\n    model.eval()\n    with torch.no_grad():\n        for i, data in tqdm(enumerate(val_loader, 0), total=len(val_loader)): \n\n            x, y = data\n\n            x = x.to(gpu)\n            y = y.to(gpu)\n\n            # forward\n            outputs = model(x)\n            loss = criterion(outputs, y)\n\n            # print statistics\n            v_loss += loss.item()\n        \n    v_loss = v_loss/len(val_loader)\n    print(f\"epoch {epoch+1} val: {v_loss}\")\n#     send_message(f\"epoch {epoch+1} val: {v_loss}\", chat_id, token)\n    scheduler.step(v_loss)\n    \n    if (v_loss) < best_loss:\n        best_loss = v_loss\n        torch.save(model,save_model_name)\n        print('Loss decreased, model saved.')\n        early_stopping_counter = 0\n    else:\n        early_stopping_counter += 1\n    if early_stopping_counter==early_stop:\n        print('Early Stopping!!!')\n        break\n    print('\\n')","metadata":{"execution":{"iopub.status.busy":"2021-10-08T04:42:56.960378Z","iopub.execute_input":"2021-10-08T04:42:56.960968Z"},"trusted":true},"execution_count":null,"outputs":[]}]}