{"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":"code","source":"import os\nimport json\nimport glob\nimport sys\nimport random\nimport collections\nimport time\nimport re\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\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\n\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import roc_auc_score","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:58.557774Z","iopub.execute_input":"2021-08-12T04:24:58.558179Z","iopub.status.idle":"2021-08-12T04:24:58.566387Z","shell.execute_reply.started":"2021-08-12T04:24:58.558141Z","shell.execute_reply":"2021-08-12T04:24:58.565194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if os.path.exists('../input/rsna-miccai-brain-tumor-radiogenomic-classification'):\n    data_directory = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\n    pytorch3dpath = '../input/efficientnetpyttorch3d/EfficientNet-PyTorch-3D'\n    \nelse:\n    data_directory = '/media/uakarsh121/data/kaggle/rsna-miccai-brain-tumor-radiogenomic-classification'\n    pytorch3dpath = \"EfficientNet-PyTorch-3D\"\n\nmri_types = ['FLAIR','T1w','T1wCE','T2w']\nSIZE = 256\nNUM_IMAGES = 64\n\nsys.path.append(pytorch3dpath)\nfrom efficientnet_pytorch_3d import EfficientNet3D","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:58.568458Z","iopub.execute_input":"2021-08-12T04:24:58.569085Z","iopub.status.idle":"2021-08-12T04:24:58.586603Z","shell.execute_reply.started":"2021-08-12T04:24:58.569045Z","shell.execute_reply":"2021-08-12T04:24:58.585415Z"},"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,voi_lut = True,rotate = 0):\n    dicom = pydicom.read_file(path)\n    data = dicom.pixel_array\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array,dicom)\n    else:\n        data =dicom.pixel_array\n    \n    if rotate>0:\n        rot_choices = [0, cv2.ROTATE_90_CLOCKWISE, cv2.ROTATE_90_COUNTERCLOCKWISE, cv2.ROTATE_180]\n        data = cv2.rotate(data,rot_choices[rotate])\n    data = cv2.resize(data,(img_size,img_size))\n    return data\n\ndef load_dicom_images_3d(scan_id,num_imgs = NUM_IMAGES,img_size = SIZE,mri_type = 'FLAIR',split='train',rotate = 0):\n    files = sorted(glob.glob(f'{data_directory}/{split}/{scan_id}/{mri_type}/*.dcm'),\n                  key =lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+',var)])\n    middle =len(files)//2\n    num_imgs2 = num_imgs//2\n    p1 = max(0,middle-num_imgs2)\n    p2 = min(len(files),middle+num_imgs2)\n    img3d = np.stack([load_dicom_image(f,rotate = rotate) for f in files[p1:p2]]).T\n    if img3d.shape[-1]<num_imgs:\n        n_zero = np.zeros((img_size,img_size,num_imgs-img3d.shape[-1]))\n        img3d = np.concatenate((img3d,n_zero),axis = -1)\n    if np.min(img3d)<np.max(img3d):\n        img3d = img3d-np.min(img3d)\n        img3d = img3d/np.max(img3d)\n    return np.expand_dims(img3d,0)\n\na = load_dicom_images_3d('00000')\nprint(a.shape)\nprint(np.min(a),np.max(a),np.mean(a),np.median(a))","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:58.588886Z","iopub.execute_input":"2021-08-12T04:24:58.589294Z","iopub.status.idle":"2021-08-12T04:24:59.065913Z","shell.execute_reply.started":"2021-08-12T04:24:58.589253Z","shell.execute_reply":"2021-08-12T04:24:59.064145Z"},"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(12)","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:59.069174Z","iopub.execute_input":"2021-08-12T04:24:59.06943Z","iopub.status.idle":"2021-08-12T04:24:59.077157Z","shell.execute_reply.started":"2021-08-12T04:24:59.069404Z","shell.execute_reply":"2021-08-12T04:24:59.076308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Train/Test splits**","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv(f'{data_directory}/train_labels.csv')\ndisplay(train_df)\n\ndf_train,df_valid = sk_model_selection.train_test_split(\ntrain_df,test_size = 0.2,random_state = 12,stratify = train_df['MGMT_value'])","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:59.078561Z","iopub.execute_input":"2021-08-12T04:24:59.07909Z","iopub.status.idle":"2021-08-12T04:24:59.107465Z","shell.execute_reply.started":"2021-08-12T04:24:59.079052Z","shell.execute_reply":"2021-08-12T04:24:59.10647Z"},"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,label_smoothing = 0.01,split='train',augment=False):\n        self.paths = paths\n        self.targets = targets\n        self.mri_type = mri_type\n        self.label_smoothing =label_smoothing\n        self.split=split\n        self.augment = augment\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),\n                                      mri_type = self.mri_type[index],split = self.split)\n        else:\n            if self.augment:\n                    rotation=np.random.randint(0,4)\n            else:\n                    rotation = 0\n            data = load_dicom_images_3d(str(scan_id).zfill(5),mri_type = self.mri_type[index],split='train',rotate = rotation)\n        if self.targets is None:\n            return {'X':torch.tensor(data).float(),'id':scan_id}\n        else:\n            y = torch.tensor(abs(self.targets[index]-self.label_smoothing),dtype = torch.float)\n            return {\"X\":torch.tensor(data).float(),'y':y}","metadata":{"execution":{"iopub.status.busy":"2021-08-12T08:23:10.178263Z","iopub.execute_input":"2021-08-12T08:23:10.178655Z","iopub.status.idle":"2021-08-12T08:23:10.191801Z","shell.execute_reply.started":"2021-08-12T08:23:10.178612Z","shell.execute_reply":"2021-08-12T08:23:10.190913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Model(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.net = EfficientNet3D.from_name('efficientnet-b0',override_params = {'num_classes':2},in_channels = 1)\n        n_features = self.net._fc.in_features\n        self.net._fc = nn.Linear(in_features  =n_features,out_features = 1,bias = True)\n        \n    def forward(self,x):\n        out =  self.net(x)\n        return out","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:59.12489Z","iopub.execute_input":"2021-08-12T04:24:59.125426Z","iopub.status.idle":"2021-08-12T04:24:59.135613Z","shell.execute_reply.started":"2021-08-12T04:24:59.125383Z","shell.execute_reply":"2021-08-12T04:24:59.134261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Trainer:\n    def __init__(self,model,device,optimizer,criterion):\n        self.model = model\n        self.device = device\n        self.optimizer = optimizer\n        self.criterion = criterion\n        \n        self.best_valid_score = np.inf\n        self.n_patience =0\n        self.lastmodel = None\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            self.info_message(\n                \"[Epoch Train: {}] loss: {:.4f}, time: {:.2f} s            \",\n                n_epoch, train_loss, train_time\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            \n            if self.best_valid_score>valid_loss:\n                self.save_model(n_epoch,save_path,valid_loss,valid_auc)\n                self.info_message(\"auc improved ffrom {:.4f} to {:.4f}. Saved model to '{}'\",\n                                 self.best_valid_score,valid_loss,self.lastmodel)\n                self.best_valid_score = valid_loss\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        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            outputs = self.model(X).squeeze(1)\n            loss = self.criterion(outputs,targets)\n            loss.backward()\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        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                X = batch['X'].to(self.device)\n                targets = batch['y'].to(self.device)\n                \n                outputs = self.model(X).squeeze(1)\n                loss = self.criterion(outputs,targets)\n                \n                sum_loss+=loss.detach().item()\n                y_all.extend(batch['y'].tolist())\n                outputs_all.extend(outputs.tolist())\n            message = 'Valid Step {}/{}, valid_loss: {:.4f}'\n            self.info_message(message, step, len(valid_loader), sum_loss/step, end=\"\\r\")\n            \n        y_all = [1 if x>0.5 else 0 for x in y_all]\n        auc = roc_auc_score(y_all,outputs_all)\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    @staticmethod\n    def info_message(message, *args, end=\"\\n\"):\n        print(message.format(*args), end=end)","metadata":{"execution":{"iopub.status.busy":"2021-08-12T04:24:59.137826Z","iopub.execute_input":"2021-08-12T04:24:59.138369Z","iopub.status.idle":"2021-08-12T04:24:59.16065Z","shell.execute_reply.started":"2021-08-12T04:24:59.138326Z","shell.execute_reply":"2021-08-12T04:24:59.159815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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    if mri_type =='all':\n        train_list = []\n        valid_list = []\n        for mri_type in mri_types:\n            df_train.loc[:,'MRI_type'] = mri_type\n            train_list.append(df_train.copy())\n            df_valid.loc[:,\"MRI_type\"] = mri_type\n            valid_list.append(df_valid.copy())\n            \n        df_train = pd.concat(train_list)\n        df_valid = pd.concat(valid_list)\n    \n    else:\n        df_train.loc[:,'MRI_type'] = mri_type\n        df_valid.loc[:,'MRI_type'] = mri_type\n    print(df_train.shape,df_valid.shape)\n    display(df_train.head())\n    \n    \n    train_data_retriever = Dataset(\n        df_train[\"BraTS21ID\"].values, \n        df_train[\"MGMT_value\"].values, \n        df_train[\"MRI_type\"].values,\n        augment=True\n    )\n\n    valid_data_retriever = Dataset(\n        df_valid[\"BraTS21ID\"].values, \n        df_valid[\"MGMT_value\"].values,\n        df_valid[\"MRI_type\"].values\n    )\n    \n    train_loader = torch_data.DataLoader(\n        train_data_retriever,\n        batch_size=4,\n        shuffle=True,\n        num_workers=8,\n    )\n\n    valid_loader = torch_data.DataLoader(\n        valid_data_retriever, \n        batch_size=4,\n        shuffle=False,\n        num_workers=8,\n    )\n\n    model = Model()\n    model.to(device)\n    \n    optimizer = torch.optim.Adam(model.parameters(), lr=0.001)\n    #optimizer = torch.optim.SGD(model.parameters(), lr=0.001, momentum=0.9)\n\n    criterion = torch_functional.binary_cross_entropy_with_logits\n\n    trainer = Trainer(\n        model, \n        device, \n        optimizer, \n        criterion\n    )\n    history = trainer.fit(\n        10, \n        train_loader, \n        valid_loader, \n        f\"{mri_type}\", \n        10,\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":{"execution":{"iopub.status.busy":"2021-08-12T04:24:59.162135Z","iopub.execute_input":"2021-08-12T04:24:59.162507Z","iopub.status.idle":"2021-08-12T08:15:05.302944Z","shell.execute_reply.started":"2021-08-12T04:24:59.162471Z","shell.execute_reply":"2021-08-12T08:15:05.300772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict Function","metadata":{}},{"cell_type":"code","source":"def predict(modelfile, df, mri_type, split):\n    print(\"Predict:\", modelfile, mri_type, df.shape)\n    df.loc[:,\"MRI_Type\"] = mri_type\n    data_retriever = Dataset(\n        df.index.values, \n        mri_type=df[\"MRI_Type\"].values,\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 = 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(batch[\"X\"].to(device))).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":"2021-08-12T08:15:05.309823Z","iopub.execute_input":"2021-08-12T08:15:05.31015Z","iopub.status.idle":"2021-08-12T08:15:05.325218Z","shell.execute_reply.started":"2021-08-12T08:15:05.310115Z","shell.execute_reply":"2021-08-12T08:15:05.324064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_valid","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_valid = df_valid.set_index(\"BraTS21ID\")\ndf_valid[\"MGMT_pred\"] = 0\nfor m, mtype in zip(modelfiles,  mri_types):\n    pred = predict(m, df_valid, mtype, \"train\")\n    df_valid[\"MGMT_pred\"] += pred[\"MGMT_value\"]\ndf_valid[\"MGMT_pred\"] /= len(modelfiles)\nauc = roc_auc_score(df_valid[\"MGMT_value\"], df_valid[\"MGMT_pred\"])\nprint(f\"Validation ensemble AUC: {auc:.4f}\")\nsns.displot(df_valid[\"MGMT_pred\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":"2021-08-12T08:23:22.372075Z","iopub.execute_input":"2021-08-12T08:23:22.372484Z","iopub.status.idle":"2021-08-12T08:25:31.02034Z","shell.execute_reply.started":"2021-08-12T08:23:22.372451Z","shell.execute_reply":"2021-08-12T08:25:31.019399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2021-08-12T08:27:58.215075Z","iopub.execute_input":"2021-08-12T08:27:58.2155Z","iopub.status.idle":"2021-08-12T08:27:58.290703Z","shell.execute_reply.started":"2021-08-12T08:27:58.215396Z","shell.execute_reply":"2021-08-12T08:27:58.289123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}