{"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":"# Over view\nBaseとなるモデルを作成する\n- model:Efficientnet 3D ver  \nref:https://www.kaggle.com/rluethy/efficientnet3d-with-one-mri-type","metadata":{}},{"cell_type":"markdown","source":"### Contents\n工事中","metadata":{}},{"cell_type":"markdown","source":"# Import dependancies","metadata":{}},{"cell_type":"code","source":"import os\nimport sys \nimport json\nimport glob\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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-06T10:38:39.259299Z","iopub.execute_input":"2021-08-06T10:38:39.25966Z","iopub.status.idle":"2021-08-06T10:38:39.26747Z","shell.execute_reply.started":"2021-08-06T10:38:39.259623Z","shell.execute_reply":"2021-08-06T10:38:39.266513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"data_directory = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\npytorch3dpath = \"../input/efficientnetpyttorch3d/EfficientNet-PyTorch-3D\"\n\nmri_types = ['FLAIR','T1w','T1wCE','T2w']\nSIZE = 256\nNUM_IMAGES = 128 \n#DUBUGのため、一旦１２８→１０に\nsys.path.append(pytorch3dpath)\nfrom efficientnet_pytorch_3d import EfficientNet3D","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:39.26936Z","iopub.execute_input":"2021-08-06T10:38:39.269874Z","iopub.status.idle":"2021-08-06T10:38:39.282034Z","shell.execute_reply.started":"2021-08-06T10:38:39.269834Z","shell.execute_reply":"2021-08-06T10:38:39.281063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Functions","metadata":{}},{"cell_type":"markdown","source":"### 1. 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    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,\n                       cv2.ROTATE_90_CLOCKWISE, \n                       cv2.ROTATE_90_COUNTERCLOCKWISE,\n                       cv2.ROTATE_180]\n        data = cv2.rotate(data, rot_choices[rotate])\n    data = cv2.resize(data, (img_size,img_size))\n    # アウトプットは　H　* W \n    return data\n\ndef load_dicom_images_3d(scan_id, \n                         num_imgs=NUM_IMAGES, \n                         img_size=SIZE, \n                         mri_type='FLAIR',\n                         split='train',\n                         rotate=0):\n    # mri_type フォルダ内の画像pathを昇順？にソートして取得\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    #　数字の昇順にsortされると思うが、なぜこれでsortできるかわからない...\n    # 数字を返すべきところでfile pathを返しているように思う\n    \n    # ---------3D配列生成----------\n    middle = len(files)//2\n    middle_num_imgs = num_imgs//2\n    p1 = max(0,middle - middle_num_imgs)\n    p2 = min(len(files),middle + middle_num_imgs)\n    \n    # 総画像の中央から必要枚数分抜き出す\n    # D * H *W * Cから　H * W * D に\n    img3d = np.stack([load_dicom_image(f,rotate=rotate)for f in files[p1:p2]]).transpose(1,2,0)\n    \n    # 総画像が必要枚数（num_imgs)無かった場合,空の配列で補完\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    \n    #0-1で正規化    \n    img3d = img3d - np.min(img3d)\n    img3d = img3d / np.max(img3d)\n    \n    # チャネルの次元を加えてあげてreturn(モノクロ画像なので一次元)\n    return np.expand_dims(img3d,0)\n\n    \n#　症例00000に対してテスト\ntest_2d = load_dicom_image(glob.glob(f'{data_directory}/train/00000/FLAIR/*.dcm')[0])\nprint(test_2d.shape)\nplt.imshow(test_2d)\ntest_3d = load_dicom_images_3d('00000')\nprint(test_3d.shape)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:39.291471Z","iopub.execute_input":"2021-08-06T10:38:39.291885Z","iopub.status.idle":"2021-08-06T10:38:40.780951Z","shell.execute_reply.started":"2021-08-06T10:38:39.291845Z","shell.execute_reply":"2021-08-06T10:38:40.779989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.set seed","metadata":{}},{"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(42)","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:40.782563Z","iopub.execute_input":"2021-08-06T10:38:40.782931Z","iopub.status.idle":"2021-08-06T10:38:40.788542Z","shell.execute_reply.started":"2021-08-06T10:38:40.78289Z","shell.execute_reply":"2021-08-06T10:38:40.787436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{}},{"cell_type":"markdown","source":"## 1. train/test splits ","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv(f'{data_directory}/train_labels.csv')\n\ndf_train, df_valid = sk_model_selection.train_test_split(\n    train_df,\n    test_size = 0.2,\n    random_state = 42,\n    #陰性、陽性の比率がtrain と　testで等しくなる用にsplit\n    stratify = train_df['MGMT_value']\n)","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:40.790538Z","iopub.execute_input":"2021-08-06T10:38:40.790909Z","iopub.status.idle":"2021-08-06T10:38:40.804147Z","shell.execute_reply.started":"2021-08-06T10:38:40.790869Z","shell.execute_reply":"2021-08-06T10:38:40.803249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. define Dataset class ","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        # paths: pd.sereis(index(patientid),path)\n        # target: 正解ラベル\n        \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            # idは５桁なのでゼロバディングする\n            \n            data = load_dicom_images_3d(str(scan_id).zfill(5), 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\n            data = load_dicom_images_3d(str(scan_id).zfill(5), mri_type=self.mri_type[index], split=\"train\", rotate=rotation)\n\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}\n","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:40.805935Z","iopub.execute_input":"2021-08-06T10:38:40.806316Z","iopub.status.idle":"2021-08-06T10:38:40.816648Z","shell.execute_reply.started":"2021-08-06T10:38:40.806277Z","shell.execute_reply":"2021-08-06T10:38:40.815649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define model","metadata":{}},{"cell_type":"code","source":"class Model(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.net = EfficientNet3D.from_name(\n            'efficientnet-b0',\n            override_params = {'num_classes':2}, \n            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-06T10:38:40.818121Z","iopub.execute_input":"2021-08-06T10:38:40.818821Z","iopub.status.idle":"2021-08-06T10:38:40.828892Z","shell.execute_reply.started":"2021-08-06T10:38:40.818783Z","shell.execute_reply":"2021-08-06T10:38:40.828013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define Trainer:","metadata":{}},{"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.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            \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 True:\n            # if self.best_valid_score < valid_auc: \n            if self.best_valid_score > valid_loss: \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_loss, self.lastmodel\n                )\n                self.best_valid_score = valid_loss\n                self.n_patience = 0\n            else:\n                self.n_patience += 1\n            # patience = epoch数になっているのがよくわからない\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            outputs = self.model(X).squeeze(1)\n            \n            loss = self.criterion(outputs, targets)\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        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\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        \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)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-06T10:38:40.832057Z","iopub.execute_input":"2021-08-06T10:38:40.832369Z","iopub.status.idle":"2021-08-06T10:38:40.853528Z","shell.execute_reply.started":"2021-08-06T10:38:40.83234Z","shell.execute_reply":"2021-08-06T10:38:40.852742Z"},"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    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    else:\n        df_train.loc[:,\"MRI_Type\"] = mri_type\n        df_valid.loc[:,\"MRI_Type\"] = mri_type\n\n    print(df_train.shape, df_valid.shape)\n    display(df_train.head())\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    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        # DEBUGのため、一旦3で回す\n        100, \n        train_loader, \n        valid_loader, \n        f\"{mri_type}\", \n        3,\n    )\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-06T10:38:40.855031Z","iopub.execute_input":"2021-08-06T10:38:40.855493Z","iopub.status.idle":"2021-08-06T12:14:23.240831Z","shell.execute_reply.started":"2021-08-06T10:38:40.855454Z","shell.execute_reply":"2021-08-06T12:14:23.239444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predict","metadata":{}},{"cell_type":"code","source":"def predict(modelfile,df,mri_type,split):\n    print('Pridict:',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-06T12:15:01.074439Z","iopub.execute_input":"2021-08-06T12:15:01.074807Z","iopub.status.idle":"2021-08-06T12:15:01.087586Z","shell.execute_reply.started":"2021-08-06T12:15:01.074771Z","shell.execute_reply":"2021-08-06T12:15:01.08645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ensemble for validation","metadata":{}},{"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'Validtion ensemble AUC:{auc:.4f}')\nsns.displot(df_valid['MGMT_pred'])","metadata":{"execution":{"iopub.status.busy":"2021-08-06T12:15:09.490843Z","iopub.execute_input":"2021-08-06T12:15:09.491172Z","iopub.status.idle":"2021-08-06T12:19:22.674461Z","shell.execute_reply.started":"2021-08-06T12:15:09.49114Z","shell.execute_reply":"2021-08-06T12:19:22.67365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ensemble for submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv(f'{data_directory}/sample_submission.csv',index_col = 'BraTS21ID')\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-06T12:19:46.63398Z","iopub.execute_input":"2021-08-06T12:19:46.634324Z","iopub.status.idle":"2021-08-06T12:22:54.750374Z","shell.execute_reply.started":"2021-08-06T12:19:46.634291Z","shell.execute_reply":"2021-08-06T12:22:54.749427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2021-08-06T12:24:59.626949Z","iopub.execute_input":"2021-08-06T12:24:59.627304Z","iopub.status.idle":"2021-08-06T12:24:59.647022Z","shell.execute_reply.started":"2021-08-06T12:24:59.627269Z","shell.execute_reply":"2021-08-06T12:24:59.645922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import FileLink\nfor modelfile in modelfiles:\n    FileLink(modelfile)","metadata":{"execution":{"iopub.status.busy":"2021-08-06T12:28:52.010617Z","iopub.execute_input":"2021-08-06T12:28:52.010943Z","iopub.status.idle":"2021-08-06T12:28:52.014722Z","shell.execute_reply.started":"2021-08-06T12:28:52.010912Z","shell.execute_reply":"2021-08-06T12:28:52.013899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}