{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":29653,"databundleVersionId":2420395,"sourceType":"competition"},{"sourceId":4696290,"sourceType":"datasetVersion","datasetId":1542467}],"dockerImageVersionId":30121,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Use stacked images (3D) and Densenet121 3D model\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\nThis notebook is based on the implementation of Densenet121 3D available here:\nhttps://www.kaggle.com/mikecho/monai-v060-deep-learning-in-healthcare-imaging\n\nIt builds 4 models with only one MRI type, then ensembles all of them computing average probabilities\n","metadata":{}},{"cell_type":"code","source":"import os\nimport sys \nimport json\nimport glob\nimport random\nimport re\nimport collections\nimport time\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\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":"2024-03-10T07:53:44.592101Z","iopub.execute_input":"2024-03-10T07:53:44.592563Z","iopub.status.idle":"2024-03-10T07:53:47.147368Z","shell.execute_reply.started":"2024-03-10T07:53:44.592464Z","shell.execute_reply":"2024-03-10T07:53:47.146603Z"},"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\"\nmonaipath = \"/kaggle/tmp/monai/\"","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:47.148629Z","iopub.execute_input":"2024-03-10T07:53:47.148908Z","iopub.status.idle":"2024-03-10T07:53:47.153539Z","shell.execute_reply.started":"2024-03-10T07:53:47.14888Z","shell.execute_reply":"2024-03-10T07:53:47.152696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir -p {monaipath}\n!cp -r {input_monaipath}/* {monaipath}","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:47.155234Z","iopub.execute_input":"2024-03-10T07:53:47.155538Z","iopub.status.idle":"2024-03-10T07:53:49.046836Z","shell.execute_reply.started":"2024-03-10T07:53:47.155513Z","shell.execute_reply":"2024-03-10T07:53:49.045951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mri_types = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\nSIZE = 256\nNUM_IMAGES = 64\nBATCH_SIZE = 4\nN_EPOCHS = 16\nSEED = 12345\nLEARNING_RATE = 0.0005\nLR_DECAY = 0.9\n\nsys.path.append(monaipath)\n\nfrom monai.networks.nets.densenet import DenseNet121","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":"2024-03-10T07:53:49.048792Z","iopub.execute_input":"2024-03-10T07:53:49.049185Z","iopub.status.idle":"2024-03-10T07:53:49.195408Z","shell.execute_reply.started":"2024-03-10T07:53:49.049142Z","shell.execute_reply":"2024-03-10T07:53:49.193723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to load images","metadata":{}},{"cell_type":"code","source":"def load_dicom_image(path, img_size=SIZE):\n    dicom = pydicom.read_file(path)\n    data = dicom.pixel_array\n    if np.min(data)==np.max(data):\n        data = np.zeros((img_size,img_size))\n        return data\n    \n    data = cv2.resize(data, (img_size, img_size))\n    return data\n\n\ndef natural_sort(l): \n    convert = lambda text: int(text) if text.isdigit() else text.lower()\n    alphanum_key = lambda key: [convert(c) for c in re.split('([0-9]+)', key)]\n    return sorted(l, key=alphanum_key)\n\n\ndef load_dicom_images_3d(scan_id, num_imgs=NUM_IMAGES, img_size=SIZE, mri_type=\"FLAIR\", split=\"train\"):\n    files = natural_sort(glob.glob(f\"{data_directory}/{split}/{scan_id}/{mri_type}/*.dcm\"))\n    \n    every_nth = len(files) / num_imgs\n    indexes = [min(int(round(i*every_nth)), len(files)-1) for i in range(0,num_imgs)]\n    \n    files_to_load = [files[i] for i in indexes]\n    \n    img3d = np.stack([load_dicom_image(f) for f in files_to_load]).T \n    \n    img3d = img3d - np.min(img3d)\n    if np.max(img3d) != 0:\n        img3d = img3d / np.max(img3d)\n    \n    return np.expand_dims(img3d,0)\n\n\nload_dicom_images_3d(\"00000\", mri_type=mri_types[0]).shape","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":"2024-03-10T07:53:49.196198Z","iopub.status.idle":"2024-03-10T07:53:49.196606Z"},"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":"2024-03-10T07:53:49.197665Z","iopub.status.idle":"2024-03-10T07:53:49.19818Z"},"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":"2024-03-10T07:53:49.199226Z","iopub.status.idle":"2024-03-10T07:53:49.199845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.tail()","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:49.20093Z","iopub.status.idle":"2024-03-10T07:53:49.201524Z"},"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, mri_type=None, split=\"train\"):\n        self.paths = paths\n        self.targets = targets\n        self.mri_type = mri_type\n        self.split = split\n          \n    def __len__(self):\n        return len(self.paths)\n    \n    def __getitem__(self, index):\n        scan_id = self.paths[index]\n        if self.targets is None:\n            data = load_dicom_images_3d(str(scan_id).zfill(5), mri_type=self.mri_type, split=self.split)\n        else:\n            data = load_dicom_images_3d(str(scan_id).zfill(5), mri_type=self.mri_type, split=\"train\")\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":"2024-03-10T07:53:49.202705Z","iopub.status.idle":"2024-03-10T07:53:49.203269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def build_model():\n    model = DenseNet121(spatial_dims=3, in_channels=1, out_channels=1)\n    return model    ","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:49.204302Z","iopub.status.idle":"2024-03-10T07:53:49.204861Z"},"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\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        \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                     \"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 = torch.tensor(batch[\"X\"]).float().to(self.device)\n            targets = batch[\"y\"].to(self.device)\n            self.optimizer.zero_grad()\n            outputs = self.model(X).squeeze(1)\n            loss = self.criterion(outputs, targets)\n                \n            loss.backward()\n\n            sum_loss += loss.detach().item()\n            \n            self.optimizer.step()\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\n                output = self.model(torch.tensor(batch[\"X\"]).float().to(self.device)).squeeze(1)\n                loss = self.criterion(output, targets)\n                sum_loss += loss.detach().item()\n                output = torch.sigmoid(output)\n\n                y_all.extend(batch[\"y\"].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, mri_type):\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":"2024-03-10T07:53:49.20642Z","iopub.status.idle":"2024-03-10T07:53:49.20701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train models","metadata":{}},{"cell_type":"code","source":"device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\ndef train_mri_type(df_train, df_valid, mri_type):\n    print(df_train.shape, df_valid.shape)\n    display(df_train.head())\n    display(df_valid.head())\n    \n    train_data_retriever = Dataset(\n        df_train[\"BraTS21ID\"].values, \n        df_train[\"MGMT_value\"].values, \n        mri_type\n    )\n\n    valid_data_retriever = Dataset(\n        df_valid[\"BraTS21ID\"].values, \n        df_valid[\"MGMT_value\"].values,\n        mri_type\n    )\n\n    train_loader = torch_data.DataLoader(\n        train_data_retriever,\n        batch_size=BATCH_SIZE,\n        shuffle=True,\n        num_workers=8,\n    )\n\n    valid_loader = torch_data.DataLoader(\n        valid_data_retriever, \n        batch_size=BATCH_SIZE,\n        shuffle=False,\n        num_workers=8,\n    )\n\n    model = build_model()\n    model.to(device)\n\n    optimizer = torch.optim.Adam(model.parameters(), lr=LEARNING_RATE)\n\n    criterion = torch_functional.binary_cross_entropy_with_logits\n\n    trainer = Trainer(\n        model, \n        device, \n        optimizer, \n        criterion\n    )\n\n    history = trainer.fit(\n        N_EPOCHS, \n        train_loader, \n        valid_loader, \n        f\"{mri_type}\", \n        N_EPOCHS,\n    )\n    \n    trainer.display_plots(mri_type)\n    \n    return trainer.lastmodel\n\nmodelfiles = None\n\nif not modelfiles:\n    modelfiles = [train_mri_type(df_train, df_valid, m) for m in mri_types]\n    print(modelfiles)","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":[],"execution":{"iopub.status.busy":"2024-03-10T07:53:49.208315Z","iopub.status.idle":"2024-03-10T07:53:49.208885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{}},{"cell_type":"code","source":"def predict(modelfile, df, mri_type, split):\n    print(\"Predict:\", modelfile, mri_type, df.shape)\n    data_retriever = Dataset(\n        df.index.values, \n        mri_type=mri_type,\n        split=split\n    )\n\n    data_loader = torch_data.DataLoader(\n        data_retriever,\n        batch_size=4,\n        shuffle=False,\n        num_workers=8,\n    )\n   \n    model = build_model()\n    model.to(device)\n    \n    checkpoint = torch.load(modelfile)\n    model.load_state_dict(checkpoint[\"model_state_dict\"])\n    model.eval()\n    \n    y_pred = []\n    ids = []\n\n    for e, batch in enumerate(data_loader,1):\n        print(f\"{e}/{len(data_loader)}\", end=\"\\r\")\n        with torch.no_grad():\n            tmp_pred = torch.sigmoid(model(torch.tensor(batch[\"X\"]).float().to(device)).squeeze(1)).cpu().numpy().squeeze()\n            if tmp_pred.size == 1:\n                y_pred.append(tmp_pred)\n            else:\n                y_pred.extend(tmp_pred.tolist())\n            ids.extend(batch[\"id\"].numpy().tolist())\n            \n    preddf = pd.DataFrame({\"BraTS21ID\": ids, \"MGMT_value\": y_pred}) \n    preddf = preddf.set_index(\"BraTS21ID\")\n    return preddf","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:49.210039Z","iopub.status.idle":"2024-03-10T07:53:49.210619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Validation","metadata":{}},{"cell_type":"code","source":"df_pred = df_valid.set_index(\"BraTS21ID\")\ndf_pred[\"MGMT_pred\"] = 0\nfor m, mtype in zip(modelfiles,  mri_types):\n    pred = predict(m, df_pred, mtype, \"train\")\n    df_pred[\"MGMT_pred\"] += pred[\"MGMT_value\"]\ndf_pred[\"MGMT_pred\"] /= len(modelfiles)\nauc = roc_auc_score(df_pred[\"MGMT_value\"], df_pred[\"MGMT_pred\"])\nprint(f\"Validation ensemble AUC: {auc:.4f}\")\nsns.displot(df_pred[\"MGMT_pred\"])","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:49.211643Z","iopub.status.idle":"2024-03-10T07:53:49.212209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv(f\"{data_directory}/sample_submission.csv\", index_col=\"BraTS21ID\")\n\nsubmission[\"MGMT_value\"] = 0\nfor m, mtype in zip(modelfiles, mri_types):\n    pred = predict(m, submission, mtype, split=\"test\")\n    submission[\"MGMT_value\"] += pred[\"MGMT_value\"]\n\nsubmission[\"MGMT_value\"] /= len(modelfiles)\nsubmission[\"MGMT_value\"].to_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-03-10T07:53:49.213388Z","iopub.status.idle":"2024-03-10T07:53:49.213966Z"},"trusted":true},"execution_count":null,"outputs":[]}]}