{"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":"## Using multimodal 2d CNN and independent augmentations for each of the MRI types\n\nThis notebook stacks number of slices per each modality and performs independent data augmentations with the use of Monai framework for each MRI type. The used architecture is 2d Densenet121 initialized with the standard pre-trained weights.\n\nAcknowledgements:\n\n- https://www.kaggle.com/rluethy/efficientnet3d-with-one-mri-type\n- https://www.kaggle.com/davidbroberts/determining-dicom-image-order\n- https://www.kaggle.com/ihelon/brain-tumor-eda-with-animations-and-modeling\n- https://www.kaggle.com/furcifer/torch-efficientnet3d-for-mri-no-train\n- https://github.com/shijianjian/EfficientNet-PyTorch-3D\n- https://www.kaggle.com/mikecho/rsna-miccai-monai-ensemble\n\nThis notebook is based on the framework Monai available here:\nhttps://www.kaggle.com/mikecho/monai-v060-deep-learning-in-healthcare-imaging\n\nNotebook to preprocess tensors:\nhttps://www.kaggle.com/mikecho/rsna-miccai-preprocessing-160x160x160/\n","metadata":{}},{"cell_type":"code","source":"import os\nimport sys \nimport json\nimport glob\nimport random\nimport re\nimport collections\nimport time\nimport math\nimport logging\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport torch\nfrom torch import nn\nfrom torch.utils import data as torch_data\nfrom sklearn import model_selection as sk_model_selection\nfrom torch.nn import functional as torch_functional\nimport torch.nn.functional as F\nfrom torch.cuda.amp import GradScaler, autocast\n\nfrom skimage.color import rgb2gray\nfrom skimage import data\nfrom skimage.filters import gaussian\n\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import roc_auc_score","metadata":{"papermill":{"duration":1.048295,"end_time":"2021-07-14T20:26:46.309722","exception":false,"start_time":"2021-07-14T20:26:45.261427","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:07:54.065959Z","iopub.execute_input":"2021-10-27T12:07:54.066756Z","iopub.status.idle":"2021-10-27T12:07:57.261349Z","shell.execute_reply.started":"2021-10-27T12:07:54.066664Z","shell.execute_reply":"2021-10-27T12:07:57.26061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir -p /root/.cache/torch/hub/checkpoints/\n!cp ../input/densenet121pretrained/densenet121-a639ec97.pth /root/.cache/torch/hub/checkpoints/","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:07:57.262866Z","iopub.execute_input":"2021-10-27T12:07:57.263131Z","iopub.status.idle":"2021-10-27T12:07:59.909998Z","shell.execute_reply.started":"2021-10-27T12:07:57.263095Z","shell.execute_reply":"2021-10-27T12:07:59.909083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pip install --upgrade ../input/mclahenumpy/mclahe-numpy","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:07:59.913103Z","iopub.execute_input":"2021-10-27T12:07:59.913651Z","iopub.status.idle":"2021-10-27T12:08:29.813703Z","shell.execute_reply.started":"2021-10-27T12:07:59.913612Z","shell.execute_reply":"2021-10-27T12:08:29.812867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from mclahe import mclahe","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:08:29.816395Z","iopub.execute_input":"2021-10-27T12:08:29.81711Z","iopub.status.idle":"2021-10-27T12:08:29.826819Z","shell.execute_reply.started":"2021-10-27T12:08:29.817065Z","shell.execute_reply":"2021-10-27T12:08:29.82602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_directory = '/kaggle/input/rsna-miccai-brain-tumor-radiogenomic-classification'\ninput_monaipath = \"/kaggle/input/monai-v060-deep-learning-in-healthcare-imaging/MONAI-0.7.0\"","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:08:29.828507Z","iopub.execute_input":"2021-10-27T12:08:29.829105Z","iopub.status.idle":"2021-10-27T12:08:29.993577Z","shell.execute_reply.started":"2021-10-27T12:08:29.829039Z","shell.execute_reply":"2021-10-27T12:08:29.992478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"parameter NUM_SLIZES indicates how many slices should be used for each modality","metadata":{}},{"cell_type":"code","source":"mri_types = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\nBATCH_SIZE = 32\nN_EPOCHS = 20\nSEED = 42\nBETA = 2.\nLEARNING_RATE = 0.0002\nLR_DECAY = 0.97\nNUM_SLICES = 20\nCHANNELS = NUM_SLICES * len(mri_types)\n\nsys.path.append(input_monaipath)\n\nfrom monai.networks.nets.densenet import DenseNet121\n\nfrom monai.data import TestTimeAugmentation\n\nfrom monai.transforms import (\n    Compose,\n    GibbsNoise,\n    RandAffine,\n    RandGibbsNoise,\n    RandKSpaceSpikeNoise,\n    RandRotate,\n    apply_transform,\n    RandAffined,\n    ConcatItemsd,\n    RandGibbsNoised,\n    RandGaussianSmoothd\n)","metadata":{"lines_to_end_of_cell_marker":2,"lines_to_next_cell":2,"papermill":{"duration":0.05565,"end_time":"2021-07-14T20:26:46.486521","exception":false,"start_time":"2021-07-14T20:26:46.430871","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:08:29.995283Z","iopub.execute_input":"2021-10-27T12:08:29.995938Z","iopub.status.idle":"2021-10-27T12:08:33.729866Z","shell.execute_reply.started":"2021-10-27T12:08:29.995893Z","shell.execute_reply":"2021-10-27T12:08:33.729097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to load images\n\nThe parameter \"randomized\" is used for additional data augmentation to provide maximum variability between epochs.\n\nFurthermore, with a probability 0.15 one of the modalities is removed (done twice).","metadata":{}},{"cell_type":"code","source":"def load_2d_tensor(study_id, randomized=False):\n    ds_path = f'../input/rsnamiccai-tensors-160x160x160-{1 if int(study_id) < 425 else 2}/kaggle/tmp/dataset/train'\n    tensors_dict = {}\n    \n    mri_types_randomized = mri_types.copy()\n    \n    if randomized:\n        if np.random.uniform() < 0.15:\n            del mri_types_randomized[np.random.randint(0, len(mri_types_randomized))]\n\n        if np.random.uniform() < 0.15:\n            del mri_types_randomized[np.random.randint(0, len(mri_types_randomized))]\n    \n    for mri_type in mri_types:\n        tensor_3d = np.load(f'{ds_path}/{study_id}/{mri_type}.npy')\n        if mri_type not in mri_types_randomized:\n            tensors_dict[mri_type] = np.zeros((NUM_SLICES, tensor_3d.shape[1], tensor_3d.shape[2]))\n        else:\n            tensors_list = []\n            step = tensor_3d.shape[-1] / NUM_SLICES\n            \n            first_index = 0\n            if randomized:\n                first_index = np.random.randint(int(math.ceil(step)))\n\n            for i in np.arange(first_index, tensor_3d.shape[-1], step):\n                tensors_list.append(tensor_3d[0, :, :, int(i)])\n\n            if len(tensors_list) > NUM_SLICES:\n                tensors_list = tensors_list[:NUM_SLICES]\n            elif len(tensors_list) < NUM_SLICES:\n                for i in range(NUM_SLICES - len(tensors_list)):\n                    tensors_list.append(tensors_list[-1])\n\n            tensors_dict[mri_type] = np.stack(tensors_list, axis=0)\n        \n    return tensors_dict","metadata":{"papermill":{"duration":0.035761,"end_time":"2021-07-14T20:26:46.726756","exception":false,"start_time":"2021-07-14T20:26:46.690995","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:20:30.904409Z","iopub.execute_input":"2021-10-27T12:20:30.905159Z","iopub.status.idle":"2021-10-27T12:20:30.928489Z","shell.execute_reply.started":"2021-10-27T12:20:30.905115Z","shell.execute_reply":"2021-10-27T12:20:30.927531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def set_seed(seed):\n    random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    if torch.cuda.is_available():\n        torch.cuda.manual_seed_all(seed)\n        torch.backends.cudnn.deterministic = True\n\nset_seed(SEED)","metadata":{"papermill":{"duration":0.668331,"end_time":"2021-07-14T20:27:48.114522","exception":false,"start_time":"2021-07-14T20:27:47.446191","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:08:33.892374Z","iopub.execute_input":"2021-10-27T12:08:33.892777Z","iopub.status.idle":"2021-10-27T12:08:33.953773Z","shell.execute_reply.started":"2021-10-27T12:08:33.89274Z","shell.execute_reply":"2021-10-27T12:08:33.953087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## train / test splits","metadata":{}},{"cell_type":"code","source":"samples_to_exclude = [109, 123, 709]\n\ntrain_df = pd.read_csv(f\"{data_directory}/train_labels.csv\")\nprint(\"original shape\", train_df.shape)\ntrain_df = train_df[~train_df.BraTS21ID.isin(samples_to_exclude)]\nprint(\"new shape\", train_df.shape)\ndisplay(train_df)\n\ndf_train, df_valid = sk_model_selection.train_test_split(\n    train_df, \n    test_size=0.2, \n    random_state=SEED, \n    stratify=train_df[\"MGMT_value\"],\n)\n","metadata":{"papermill":{"duration":0.633753,"end_time":"2021-07-14T20:27:49.350524","exception":false,"start_time":"2021-07-14T20:27:48.716771","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:08:33.956995Z","iopub.execute_input":"2021-10-27T12:08:33.9572Z","iopub.status.idle":"2021-10-27T12:08:33.998446Z","shell.execute_reply.started":"2021-10-27T12:08:33.957175Z","shell.execute_reply":"2021-10-27T12:08:33.997665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.tail()","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:08:34.001211Z","iopub.execute_input":"2021-10-27T12:08:34.001478Z","iopub.status.idle":"2021-10-27T12:08:34.010682Z","shell.execute_reply.started":"2021-10-27T12:08:34.001445Z","shell.execute_reply":"2021-10-27T12:08:34.009806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model and training classes","metadata":{}},{"cell_type":"code","source":"class Dataset(torch_data.Dataset):\n    def __init__(self, paths, targets=None, transforms=None, split=\"train\"):\n        self.paths = paths\n        self.targets = targets\n        self.split = split\n        self.transforms = transforms\n        \n        if transforms is not None:\n            self.transforms.set_random_state(seed=SEED)\n          \n    def __len__(self):\n        return len(self.paths)  #BATCH_SIZE * 2 # len(self.paths)\n    \n    def __getitem__(self, index):\n        scan_id = self.paths[index]\n        data = load_2d_tensor(str(scan_id).zfill(5), randomized=(self.split=='train'))\n        \n        if self.transforms is not None:\n            data = apply_transform(self.transforms, data, map_items=False)\n            data = data[\"img\"]\n            \n        if not isinstance(data, torch.Tensor):\n            data = torch.tensor(data, dtype=torch.float)\n            \n        if self.targets is None:\n            return {\"X\": data, \"id\": scan_id}\n        else:\n            return {\"X\": data, \"y\": torch.tensor(self.targets[index], dtype=torch.float)}\n","metadata":{"papermill":{"duration":0.634322,"end_time":"2021-07-14T20:27:50.594701","exception":false,"start_time":"2021-07-14T20:27:49.960379","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:44:46.299355Z","iopub.execute_input":"2021-10-27T12:44:46.299647Z","iopub.status.idle":"2021-10-27T12:44:46.309095Z","shell.execute_reply.started":"2021-10-27T12:44:46.299616Z","shell.execute_reply":"2021-10-27T12:44:46.308389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def build_model():\n    model = DenseNet121(spatial_dims=2, in_channels=CHANNELS, out_channels=1, pretrained=True)\n    return model    ","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:08:34.027255Z","iopub.execute_input":"2021-10-27T12:08:34.027545Z","iopub.status.idle":"2021-10-27T12:08:34.03411Z","shell.execute_reply.started":"2021-10-27T12:08:34.027516Z","shell.execute_reply":"2021-10-27T12:08:34.033356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Trainer:\n    def __init__(\n        self, \n        model, \n        device, \n        optimizer, \n        criterion\n    ):\n        self.model = model\n        self.device = device\n        self.optimizer = optimizer\n        self.lr_scheduler = torch.optim.lr_scheduler.ExponentialLR(self.optimizer, gamma=LR_DECAY)\n        self.criterion = criterion\n        self.scaler = GradScaler()\n\n        self.best_valid_score = .0\n        self.n_patience = 0\n        self.lastmodel = None\n        \n        self.val_losses = []\n        self.train_losses = []\n        self.val_auc = []\n        self.val_auc_tta = []\n        \n    def fit(self, epochs, train_loader, valid_loader, save_path, patience):      \n        for n_epoch in range(1, epochs + 1):\n            self.info_message(\"EPOCH: {}\", n_epoch)\n            \n            train_loss, train_time = self.train_epoch(train_loader)\n            valid_loss, valid_auc, valid_time = self.valid_epoch(valid_loader)\n            \n            self.train_losses.append(train_loss)\n            self.val_losses.append(valid_loss)\n            self.val_auc.append(valid_auc)\n            \n            self.info_message(\n                \"[Epoch Train: {}] loss: {:.4f}, time: {:.2f} s            \",\n                n_epoch, train_loss, train_time\n            )\n            \n            self.info_message(\n                \"[Epoch Valid: {}] loss: {:.4f}, auc: {:.4f}, time: {:.2f} s\",\n                n_epoch, valid_loss, valid_auc, valid_time\n            )\n\n            if self.best_valid_score < valid_auc: \n                self.save_model(n_epoch, save_path, valid_loss, valid_auc)\n                self.info_message(\n                     \"val_auc improved from {:.4f} to {:.4f}. Saved model to '{}'\", \n                    self.best_valid_score, valid_auc, self.lastmodel\n                )\n                self.best_valid_score = valid_auc\n                self.n_patience = 0\n            else:\n                self.n_patience += 1\n            \n            if self.n_patience >= patience:\n                self.info_message(\"\\nValid auc didn't improve last {} epochs.\", patience)\n                break\n            \n    def train_epoch(self, train_loader):\n        self.model.train()\n        t = time.time()\n        sum_loss = 0\n\n        for step, batch in enumerate(train_loader, 1):\n            X = batch[\"X\"].to(self.device)\n            targets = batch[\"y\"].to(self.device)\n            self.optimizer.zero_grad()\n            with autocast():\n                outputs = self.model(X).squeeze(1)\n                loss = self.criterion(outputs, targets)\n                \n            self.scaler.scale(loss).backward()\n\n            sum_loss += loss.detach().item()\n\n            self.scaler.step(self.optimizer)\n            \n            self.scaler.update()\n            \n            message = 'Train Step {}/{}, train_loss: {:.4f}'\n            self.info_message(message, step, len(train_loader), sum_loss/step, end=\"\\r\")\n            \n        self.lr_scheduler.step()\n        \n        return sum_loss/len(train_loader), int(time.time() - t)\n    \n    def valid_epoch(self, valid_loader):\n        self.model.eval()\n        t = time.time()\n        sum_loss = 0\n        y_all = []\n        outputs_all = []\n\n        for step, batch in enumerate(valid_loader, 1):\n            with torch.no_grad():\n                targets = batch[\"y\"].to(self.device)\n                output = self.model(batch[\"X\"].to(self.device)).squeeze(1)\n                loss = self.criterion(output, targets)\n                output = torch.sigmoid(output)\n                sum_loss += loss.detach().item()\n\n                y_all.extend(targets.tolist())\n                outputs_all.extend(output.tolist())\n\n            message = 'Valid Step {}/{}, valid_loss: {:.4f}'\n            self.info_message(message, step, len(valid_loader), sum_loss/step, end=\"\\r\")\n            \n        auc = roc_auc_score(y_all, outputs_all)\n        \n        return sum_loss/len(valid_loader), auc, int(time.time() - t)\n    \n    def save_model(self, n_epoch, save_path, loss, auc):\n        self.lastmodel = f\"{save_path}-e{n_epoch}-loss{loss:.3f}-auc{auc:.3f}.pth\"\n        torch.save(\n            {\n                \"model_state_dict\": self.model.state_dict(),\n                \"optimizer_state_dict\": self.optimizer.state_dict(),\n                \"best_valid_score\": self.best_valid_score,\n                \"n_epoch\": n_epoch,\n            },\n            self.lastmodel,\n        )\n        \n    def display_plots(self):\n        plt.figure(figsize=(10,5))\n        plt.title(\"Training and Validation Loss\")\n        plt.plot(self.val_losses,label=\"val\")\n        plt.plot(self.train_losses,label=\"train\")\n        plt.xlabel(\"iterations\")\n        plt.ylabel(\"Loss\")\n        plt.legend()\n        plt.show()\n        plt.close()\n        \n        plt.figure(figsize=(10,5))\n        plt.title(\"Validation AUC-ROC\")\n        plt.plot(self.val_auc,label=\"val\")\n        plt.xlabel(\"iterations\")\n        plt.ylabel(\"AUC\")\n        plt.legend()\n        plt.show()\n        plt.close()\n    \n    @staticmethod\n    def info_message(message, *args, end=\"\\n\"):\n        print(message.format(*args), end=end)","metadata":{"papermill":{"duration":0.637077,"end_time":"2021-07-14T20:27:58.09407","exception":false,"start_time":"2021-07-14T20:27:57.456993","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-27T12:08:34.037894Z","iopub.execute_input":"2021-10-27T12:08:34.038135Z","iopub.status.idle":"2021-10-27T12:08:34.06615Z","shell.execute_reply.started":"2021-10-27T12:08:34.038098Z","shell.execute_reply":"2021-10-27T12:08:34.065372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Loss function with DCA penalty for better model calibration","metadata":{}},{"cell_type":"code","source":"# https://github.com/GB-TonyLiang/DCA\ndef cross_entropy_with_dca_loss(logits, labels, pos_weights=None, alpha=1., beta=BETA):        \n    ce = F.binary_cross_entropy_with_logits(logits, labels, pos_weight=None)\n\n    confidences = torch.sigmoid(logits)\n    predictions = torch.round(confidences)\n    accuracies = predictions.eq(labels)\n    mean_conf = confidences.float().mean()\n    acc = accuracies.float().sum() / len(accuracies)\n    dca = torch.abs(mean_conf - acc)\n    loss = alpha * ce + beta * dca\n    \n    return loss","metadata":{"execution":{"iopub.status.busy":"2021-10-27T12:08:34.067664Z","iopub.execute_input":"2021-10-27T12:08:34.067958Z","iopub.status.idle":"2021-10-27T12:08:34.07838Z","shell.execute_reply.started":"2021-10-27T12:08:34.067921Z","shell.execute_reply":"2021-10-27T12:08:34.07749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## train models\n\nTrain models with augmented data where each modality is independently transformed with the use of Monai dictionary transformations","metadata":{}},{"cell_type":"code","source":"device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nprint(df_train.shape, df_valid.shape)\ndisplay(df_train.head())\ndisplay(df_valid.head())\n\ntrain_transforms = Compose([\n    RandAffined(mri_types, prob=0.5, rotate_range=np.pi, shear_range=0.05, translate_range=0.05, scale_range=0.05),\n    RandGibbsNoised(mri_types),\n    RandGaussianSmoothd(mri_types),\n    ConcatItemsd(mri_types, \"img\")\n])\n\nvalid_transforms = Compose([\n    ConcatItemsd(mri_types, \"img\")\n])\n\ntrain_data_retriever = Dataset(\n    df_train[\"BraTS21ID\"].values, \n    df_train[\"MGMT_value\"].values, \n    transforms=train_transforms\n)\n\nvalid_data_retriever = Dataset(\n    df_valid[\"BraTS21ID\"].values, \n    df_valid[\"MGMT_value\"].values,\n    transforms=valid_transforms\n)\n\ntrain_loader = torch_data.DataLoader(\n    train_data_retriever,\n    batch_size=BATCH_SIZE,\n    shuffle=True,\n    num_workers=2,\n)\n\nvalid_loader = torch_data.DataLoader(\n    valid_data_retriever, \n    batch_size=BATCH_SIZE,\n    shuffle=False,\n    num_workers=2,\n)\n\nmodel = build_model()\nmodel.to(device)\n\noptimizer = torch.optim.Adam(model.parameters(), lr=LEARNING_RATE)\n\ncriterion = cross_entropy_with_dca_loss\n\ntrainer = Trainer(\n    model, \n    device, \n    optimizer, \n    criterion\n)\n\nhistory = trainer.fit(\n    N_EPOCHS, \n    train_loader, \n    valid_loader, \n    f\"model2d\", \n    N_EPOCHS,\n)\n\ntrainer.display_plots()","metadata":{"lines_to_next_cell":2,"papermill":{"duration":447.387602,"end_time":"2021-07-14T20:35:26.110421","exception":false,"start_time":"2021-07-14T20:27:58.722819","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-10-27T12:44:52.28097Z","iopub.execute_input":"2021-10-27T12:44:52.281837Z"},"trusted":true},"execution_count":null,"outputs":[]}]}