{"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":"# RSNA MICCAI Brain Tumor classifier -- Ordered MRI voxel data\n\n## Acknowledgement:\n\nVoxel Data Ordering\n* https://www.kaggle.com/davidbroberts/determining-dicom-image-order\n\nMobileNetV3 for 2D images\n* https://arxiv.org/abs/1905.02244\n* https://github.com/pytorch/vision/blob/v0.9.0/torchvision/models/mobilenetv3.py\n\nTroubleshooting\n* https://fullstackdeeplearning.com/spring2021/lecture-7/\n\nOther Works on 3D MRI voxel data\n* http://www.ajnr.org/content/ajnr/42/5/845.full.pdf","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport pydicom\nimport cv2 as cv\n\nimport torch\nfrom torch import nn\nfrom torch.nn import functional as F\nfrom torch.utils.data import DataLoader, Dataset, Subset\n# from torchvision import models\nimport torchvision\nimport kornia as K  # batch image augmentations with torch.Tensor\nfrom kornia.augmentation import AugmentationSequential\nfrom kornia.augmentation.base import AugmentationBase3D  # Subclassing this is too complicated.\nfrom kornia.enhance import invert\n\nfrom tqdm.notebook import tqdm\n\nfrom pathlib import Path\nfrom typing import Union, Tuple, List, Optional, Type, Dict, Iterable\nimport time\n\nDEBUG = False\nREPRODUCTIVE = True\nINFERENCE_ONLY = True\nUSE_CROSS_VALIDATION = True\n\nrandom_state = 42\nmodel_name = \"Net-3D\"\ndata_dir = Path(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification\")\nmodels_dir = Path(\"../input/model-weights-for-rsna-miccai-brain-tumor-dataset\")\n# models_dir = Path(\".\")  # If train model with local machine\n\ndevice = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n\ntime_begin = time.time()\n\nif REPRODUCTIVE:\n    np.random.seed(random_state)\n    torch.random.manual_seed(random_state)\ndisplay(list(data_dir.iterdir()), torch.__version__, torchvision.__version__)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data manipulation","metadata":{}},{"cell_type":"code","source":"mri_series = {0: \"FLAIR\", 1: \"T1w\", 2: \"T1wCE\", 3: \"T2w\"}\nmri_series_map = {v: k for k, v in mri_series.items()}\nplanes = {0: \"Unknown\", 1: \"Coronal\", 2: \"Sagittal\", 3: \"Axial\"}\nplanes_map = {v: k for k, v in planes.items()}","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels_train = pd.read_csv(data_dir / \"train_labels.csv\", dtype={\"BraTS21ID\": str})\nlabels_train","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def look_one_dcm(instance_id: str, img_dir: Path, mri_series=\"FLAIR\", verbose=False):\n    dcm_paths = list(img_dir.glob(\"./{}/{}/*.dcm\".format(instance_id.zfill(5), mri_series)))\n    print(\"Containing {} dicom files(including blank).\".format(len(dcm_paths)))\n    if dcm_paths:\n        dcm_mid = dcm_paths[(len(dcm_paths) - 1) // 2]\n        dcm_ds = pydicom.read_file(str(dcm_mid))\n        if verbose:\n            print(dir(dcm_ds))\n            print(dcm_ds)\n            print(type(dcm_ds[(\"0010\", \"0010\")].value))\n            print(dcm_ds[(\"0020\", \"0032\")].name, eval(str(dcm_ds[(\"0020\", \"0032\")].value)))\n            print(dir(dcm_ds[(\"0020\", \"0032\")]))\n            print(dcm_ds.pixel_array.dtype)\n        plt.imshow(dcm_ds.pixel_array, cmap=plt.cm.gray)\n        plt.show()\n\n\nlook_one_dcm(\"00000\", data_dir / \"train\", verbose=True)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_image_plane(loc):\n    row_x, row_y, row_z, col_x, col_y, col_z = [round(v) for v in loc]\n    if (row_x, row_y, col_x, col_y) == (1, 0, 0, 0): return planes[1]\n    if (row_x, row_y, col_x, col_y) == (0, 1, 0, 0): return planes[2]\n    if (row_x, row_y, col_x, col_y) == (1, 0, 0, 1): return planes[3]\n    return planes[0]\n\n\nclass DICOMMetaLoader(Dataset):\n    \n    def __init__(self, img_dir: Path, glob=None):\n        super(DICOMMetaLoader, self).__init__()\n        if glob is None:\n            glob = \"./*/*/*.dcm\"\n        self.dcm_paths = list(img_dir.glob(glob))\n    \n    def __len__(self): return len(self.dcm_paths)\n    \n    def __getitem__(self, idx) -> dict:\n        dcm_path = str(self.dcm_paths[idx])\n        dcm_obj = pydicom.read_file(dcm_path)\n        photometric = str(dcm_obj[0x28, 0x04])\n        array = dcm_obj.pixel_array\n        if photometric == \"MONOCHROME1\":\n            info_func = np.iinfo if np.issubdtype(array.dtype, np.integer) else np.finfo\n            array = info_func(array.dtype).max - array\n        image_mean, image_std = np.mean(array), np.std(array)\n        \n        impo_x, impo_y, impo_z = [float(v) for v in dcm_obj[0x20, 0x32]]\n        plane = get_image_plane(dcm_obj[0x20, 0x37])\n        \n        patient_id = str(dcm_obj[0x0010, 0x0020].value).strip().zfill(5)\n        series_desc = str(dcm_obj[0x0008, 0x103e].value).strip()\n        row = dict(dcm_path=dcm_path, BraTS21ID=patient_id, series_description=series_desc,\n                   image_mean=image_mean, image_std=image_std,\n                   plane=plane,\n                   image_position_x=impo_x, image_position_y=impo_y, image_position_z=impo_z)\n        return row\n\n\ndef get_meta_from_glob(img_dir: Path, glob=None) -> pd.DataFrame:\n    dcm_ds = DICOMMetaLoader(img_dir, glob)\n    dcm_dl = DataLoader(dcm_ds, batch_size=256, num_workers=6)\n    df = pd.DataFrame()\n    for item in tqdm(dcm_dl):\n        chunks = pd.DataFrame.from_dict({k:np.asarray(v) for k, v in item.items()})\n        df = pd.concat([df, chunks], ignore_index=True)\n    return df\n\n\ndf_train = get_meta_from_glob(data_dir / \"train\")","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To categorical data by mapping, \n\ndf_train.loc[:, \"plane\"] = df_train.loc[:, \"plane\"].map(planes_map)\ndf_train.loc[:, \"series_description\"] = df_train.loc[:, \"series_description\"].map(mri_series_map)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def keep_non_blank(df: pd.DataFrame):\n    \"\"\"\n    Keep data containing non blank image.\n    :params:\n        df: pd.DataFrame, requires \"image_std\" and \"image_mean\" in df.columns.\n    :returns:\n        pd.DataFrame: filtered DataFrame\n    \"\"\"\n    df = df.loc[(df[\"image_std\"] > 0) & (df[\"image_mean\"] > 0)]\n    return df\n\n\ndisplay(len(df_train))\ndf_train = keep_non_blank(df_train)\ndisplay(len(df_train))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def drop_by_id(df: pd.DataFrame, ids: List[Union[int, str]]):\n    ids = [str(s).zfill(5) for s in ids]\n    df = df.loc[~(df[\"BraTS21ID\"].isin(ids))].reset_index(drop=True)\n    return df\n\n\ndrop_ids = \"00109, 00123, 00709\".split(\", \")\ndf_train = drop_by_id(df_train, drop_ids)\nlabels_train = drop_by_id(labels_train, drop_ids)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def count_values(df: pd.DataFrame):\n    groupby = df.groupby([\"BraTS21ID\", \"series_description\"])\n    count = groupby.count()\n    display(count[\"dcm_path\"].describe())\n    display(count.loc[count[\"dcm_path\"] == count[\"dcm_path\"].min(), \"dcm_path\"])\n    display(count.loc[count[\"dcm_path\"] == count[\"dcm_path\"].max(), \"dcm_path\"])\n\n\ndisplay(df_train.describe())\ncount_values(df_train)\nlook_one_dcm(\"00571\", data_dir / \"train\", mri_series[0])\nlook_one_dcm(\"00818\", data_dir / \"train\", mri_series[0])\nlook_one_dcm(\"00012\", data_dir / \"train\", mri_series[3])","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def is_retrievable(df: pd.DataFrame,\n                   patient_id: str,\n                   series_desc_idx: int):\n    retrieved_idx = (df[\"BraTS21ID\"].eq(patient_id)) & (df[\"series_description\"].eq(series_desc_idx))\n    return True if retrieved_idx.sum() > 0 else False\n\n\ndef get_voxel_by_id_series(df: pd.DataFrame,\n                           patient_id: str,\n                           series_desc_idx: int = 0,\n                           size: Union[int, Tuple[int, int]] = 256) -> Tuple[np.ndarray, int]:\n    \"\"\"\n    :params:\n        :df: required columns: [dcm_path, BraTS21ID, series_description, plane,\n                                image_position_x, image_position_y, image_position_z]\n    \"\"\"\n    size = (int(size), int(size)) if isinstance(size, (int, float)) else size\n    retrieved_idx = (df[\"BraTS21ID\"].eq(patient_id)) & (df[\"series_description\"].eq(series_desc_idx))\n    assert retrieved_idx.sum() > 0, \"Nothing retrived.\"\n    retrieved_df = df.loc[retrieved_idx].copy()\n    plane = retrieved_df[\"plane\"].unique()\n    assert len(plane) == 1, \"Different plane in a folder.\"\n    img_pos_cols = [c for c in retrieved_df.columns if c.startswith(\"image_position_\")]\n    img_pos_stds = np.array([retrieved_df[c].std() for c in img_pos_cols])\n    img_pos_argsort = np.argsort(img_pos_stds)[::-1]\n    sorted_df = retrieved_df.sort_values([img_pos_cols[i] for i in img_pos_argsort], ascending=True, ignore_index=True)\n    voxel_stack = list()\n    for row in sorted_df.itertuples():\n        dcm_obj = pydicom.read_file(row.dcm_path)\n        array = dcm_obj.pixel_array\n        array = cv.resize(array, size)\n        dinfo = np.iinfo(array.dtype) if np.issubdtype(array.dtype, np.integer) else np.finfo(array.dtype)\n        array = (array / dinfo.max).astype(np.float32)  # like (a / 255) if a.dtype is uint8\n        if dcm_obj[0x0028, 0x0004] == \"MONOCHROME1\":\n            array = dinfo.max - array\n        voxel_stack.append(array)\n    voxel = np.stack(voxel_stack)\n    voxel = (voxel - np.min(voxel)) / max(np.max(voxel), 1e-8)  # min-max normalization\n    return voxel, plane[0]\n\n\ndef plot_voxel(voxel, max_n_plots=10, cols=10):\n    actual_n_plots = min(max_n_plots, len(voxel))\n    rows = int(np.ceil(actual_n_plots / cols))\n    fig, axes = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows), tight_layout=True)\n    for i in range(actual_n_plots):\n        axes[i // cols, i % cols].imshow(voxel[i, :, :], cmap=plt.cm.gray)\n        axes[i // cols, i % cols].set_axis_off()\n    plt.show()\n\n\nvox, plane = get_voxel_by_id_series(df_train, \"00571\", size=256)\nplot_voxel(vox, 14, 7)\nvox, plane = get_voxel_by_id_series(df_train, \"00571\", size=128)\nplot_voxel(vox, 14, 7)\nvox, plane = get_voxel_by_id_series(df_train, \"00571\", size=64)\nplot_voxel(vox, 14, 7)\nvox, plane = get_voxel_by_id_series(df_train, \"00571\", size=32)\nplot_voxel(vox, 14, 7)\ndisplay(plane)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class MRIVoxelDataset(Dataset):\n    \n    def __init__(self, meta_df: pd.DataFrame, label_df: Optional[pd.DataFrame] = None,\n                 voxel_size: Union[int, Tuple[int, int], Tuple[int, int, int]] = (64, 256, 256),\n                 including_series: np.ndarray = np.array(list(mri_series.keys()), dtype=np.int64)):\n        \"\"\"\n        :params:\n            :meta_df: required columns: [dcm_path, BraTS21ID, series_description, plane,\n                                         image_position_x, image_position_y, image_position_z]\n            :label_df(Optional): required columns: [BraTS21ID, MGMT_value]\n            :voxel_size: if int, the D, H, W will be set to the same;\n                         if (int, int), D by voxel_size[0], H, W by voxel_size[1];\n                         if (int, int, int), D, H, W will be set respectively.\n        \"\"\"\n        super(MRIVoxelDataset, self).__init__()\n        self.meta_df,self.label_df,self.voxel_size = meta_df,label_df,voxel_size\n        self.including_series = including_series\n        if isinstance(self.voxel_size, int):\n            self.voxel_size = tuple(self.voxel_size for _ in range(3))\n        elif isinstance(self.voxel_size, tuple):\n            if len(self.voxel_size) == 2:\n                self.voxel_size = (self.voxel_size[0], self.voxel_size[1], self.voxel_size[1])\n        self.meta_df = self.meta_df.loc[self.meta_df[\"series_description\"].isin(self.including_series)].copy()\n        if self.label_df is None:\n            self.label_df = pd.concat([pd.DataFrame.from_dict(\n                dict(BraTS21ID=self.meta_df[\"BraTS21ID\"].unique())\n            )], axis=1)\n            self.label_df.loc[:, \"BraTS21ID\"] = self.label_df[\"BraTS21ID\"].map(lambda i: str(i).zfill(5))\n            labels = np.full_like(self.label_df[\"BraTS21ID\"].values, np.nan, dtype=np.float64)\n            self.label_df.loc[:, \"MGMT_value\"] = labels\n\n        new_label_df = pd.DataFrame()\n        for v in self.meta_df[\"series_description\"].unique():\n            series_desc = pd.DataFrame({\n                    \"series_description\": np.full((len(self.label_df)), v, dtype=np.int64)\n                 })\n            df = self.label_df.reset_index(drop=True)\n            df = pd.concat([df, series_desc], axis=1)\n            new_label_df = pd.concat([new_label_df, df], axis=0)\n        self.label_df = new_label_df.reset_index(drop=True)\n\n        retrievables = list()\n        for i in range(len(self.label_df)):\n            row = self.label_df.iloc[i]\n            flag = is_retrievable(self.meta_df, row.BraTS21ID, row.series_description)\n            if not flag:\n                print(row.BraTS21ID, row.series_description)\n            retrievables.append(flag)\n        retrievables = np.asarray(retrievables)\n        self.label_df = self.label_df.iloc[retrievables]\n        print(f\"Got {len(self)} samples in dataset.\")\n\n    def __len__(self): return len(self.label_df)\n    \n    def __getitem__(self, idx):\n        row = self.label_df.iloc[idx]\n        voxel, plane = get_voxel_by_id_series(self.meta_df, row[\"BraTS21ID\"], row[\"series_description\"], self.voxel_size[1:])\n        voxel = torch.tensor(voxel, dtype=torch.float32).unsqueeze(0).unsqueeze(0)  # [N, C, D, H, W]\n        voxel = F.interpolate(voxel, self.voxel_size, mode=\"trilinear\", align_corners=False)\n        voxel = voxel.squeeze(0)\n        label = torch.tensor([row[\"MGMT_value\"]], dtype=torch.float32)\n        plane = torch.tensor(plane, dtype=torch.int64)\n        series_desc = torch.tensor(row[\"series_description\"], dtype=torch.int64)\n        return voxel, label, (series_desc, plane)\n\n\nif DEBUG:\n    ds_ = MRIVoxelDataset(df_train, labels_train, (64, 128), np.array([0], dtype=np.int64))\n    dl_ = DataLoader(ds_, batch_size=4, num_workers=4)\n    for voxel, label, (series_desc, plane) in dl_:\n        print(voxel.shape, label.shape, plane.shape, series_desc.shape)\n        print(voxel.dtype, label.dtype, plane.dtype, series_desc.dtype)\n        break","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The architecture of model","metadata":{}},{"cell_type":"code","source":"NormLayerClass = Type\nActivationLayerClass = Type\n\n\nclass SqueezeExcitation(nn.Module):\n    \n    def __init__(self, in_channels):\n        super(SqueezeExcitation, self).__init__()\n        self.in_channels = in_channels\n        self.squeeze_channels = self.in_channels // 4\n        \n        self.seq = nn.Sequential(\n            nn.AdaptiveAvgPool3d(1),\n            nn.Conv3d(self.in_channels, self.squeeze_channels, 1),\n            nn.ReLU(inplace=True),\n            nn.Conv3d(self.squeeze_channels, self.in_channels, 1),\n            nn.Hardsigmoid(inplace=True),\n        )\n    \n    def forward(self, x):\n        scale = self.seq(x)\n        out = scale * x\n        return out\n\n\nclass ConvBNActivation(nn.Module):\n    \n    def __init__(self, conv_config: dict,\n                 norm_layer_cls: NormLayerClass = nn.BatchNorm3d,\n                 activation_layer_cls: ActivationLayerClass = nn.ReLU,\n                 use_se: bool = False,\n        ) -> None:\n        super(ConvBNActivation, self).__init__()\n        layers = list()\n        layers.append(nn.Conv3d(**conv_config))\n        layers.append(norm_layer_cls(conv_config[\"out_channels\"]))\n        layers.append(activation_layer_cls(inplace=True))\n        if use_se:\n            layers.append(SqueezeExcitation(conv_config[\"out_channels\"]))\n        self.seq = nn.Sequential(*layers)\n    \n    def forward(self, x):\n        return self.seq(x)\n    \n    @staticmethod\n    def config(in_channels: int,\n               out_channels: int,\n               kernel_size: Union[int, Tuple[int, int, int]],\n               stride: Union[int, Tuple[int, int, int]] = 1,\n               padding: Union[int, Tuple[int, int, int]] = 0,\n               dilation: Union[int, Tuple[int, int, int]] = 1,\n               groups: int = 1,\n               bias: bool = True,\n               padding_mode: str = 'zeros',\n        ) -> dict:\n        return locals()\n\n\nclass BottleNeck(nn.Module):\n    \n    def __init__(self, residual_config: dict):\n        super(BottleNeck, self).__init__()\n        self.residual_config = residual_config\n        layers = list()\n        layers.append(ConvBNActivation(ConvBNActivation.config(\n            self.residual_config[\"in_channels\"],\n            self.residual_config[\"expand_channels\"],\n            1,\n            1,\n            0,\n        ), self.residual_config[\"norm_layer_cls\"], self.residual_config[\"activation_layer_cls\"]))\n        layers.append(ConvBNActivation(ConvBNActivation.config(\n            self.residual_config[\"expand_channels\"],\n            self.residual_config[\"expand_channels\"],\n            self.residual_config[\"kernel_size\"],\n            self.residual_config[\"stride\"],\n            self.residual_config[\"padding\"],\n            groups=self.residual_config[\"expand_channels\"],\n        ), self.residual_config[\"norm_layer_cls\"],\n           self.residual_config[\"activation_layer_cls\"],\n           self.residual_config[\"use_se\"]))\n        layers.append(ConvBNActivation(ConvBNActivation.config(\n            self.residual_config[\"expand_channels\"],\n            self.residual_config[\"out_channels\"],\n            1,\n            1,\n            0,\n        ), self.residual_config[\"norm_layer_cls\"], nn.Identity))\n        self.seq = nn.Sequential(*layers)\n        # The shortcut: Same as nn.Linear if channels at last dim.\n        self.shortcut = nn.Conv3d(self.residual_config[\"in_channels\"], self.residual_config[\"out_channels\"], 1)\n    \n    def forward(self, x):\n        post_seq = self.seq(x)\n        x = self.shortcut(x)\n        x = F.interpolate(x, post_seq.shape[-3:], mode=\"trilinear\", align_corners=False)\n        return x + post_seq\n    \n    @staticmethod\n    def config(in_channels: int,\n               out_channels: int,\n               expand_channels: int,\n               kernel_size: Union[int, Tuple[int, int, int]],\n               stride: Union[int, Tuple[int, int, int]] = 1,\n               padding: Union[int, Tuple[int, int, int]] = 0,\n               norm_layer_cls: NormLayerClass = nn.BatchNorm3d,\n               activation_layer_cls: ActivationLayerClass = nn.Hardswish,\n               use_se: bool = False,\n    ) -> dict:\n        return locals()\n\n\nclass NetFeatures(nn.Module):\n    \n    def __init__(self, in_channels, out_channels, residual_config_list: List[dict]):\n        super(NetFeatures, self).__init__()\n        self.in_channels,self.out_channels = in_channels,out_channels\n        self.residual_config_list = residual_config_list\n\n        first_conv_out_channels = self.residual_config_list[0][\"in_channels\"]\n        self.first_conv = ConvBNActivation(ConvBNActivation.config(\n            self.in_channels, first_conv_out_channels, 3, 2, 1), activation_layer_cls=nn.ReLU)\n        residual_layers = list()\n        for conf in self.residual_config_list:\n            residual_layers.append(BottleNeck(conf))\n        self.residual_block = nn.Sequential(*residual_layers)\n        last_conv_in_channels = self.residual_config_list[-1][\"out_channels\"]\n        self.last_conv = ConvBNActivation(ConvBNActivation.config(last_conv_in_channels, self.out_channels, 1),\n                                          nn.BatchNorm3d,\n                                          nn.Hardswish,\n                                          use_se=True)\n    \n    def forward(self, x):\n        x = self.first_conv(x)\n        x = self.residual_block(x)\n        x = self.last_conv(x)\n        return x\n\n\nclass ConcatEmbeddingLinear(nn.Module):\n    \n    def __init__(self, in_features: int, out_features: int, n_embeddings: int, embed_dim: Optional[int] = None):\n        super(ConcatEmbeddingLinear, self).__init__()\n        self.in_features,self.out_features = in_features,out_features\n        self.n_embeddings,self.embed_dim = n_embeddings,embed_dim\n        if self.embed_dim is None: self.embed_dim = self.in_features\n        \n        self.emb = nn.Embedding(self.n_embeddings, self.embed_dim)\n        self.fc = nn.Linear(self.in_features + self.embed_dim, self.out_features)\n    \n    def forward(self, x, idx_emb):\n        emb_out = self.emb(idx_emb)\n        concatenated = torch.cat([emb_out, x], dim=-1)\n        out = self.fc(concatenated)\n        return out\n\n\nclass Net(nn.Module):\n    \n    def __init__(self, in_channels, feature_out_channels, hidden_features, n_classes, n_series, n_planes,\n                 residual_config_list: List[dict]) -> None:\n        super(Net, self).__init__()\n        self.in_channels,self.feature_out_channels,self.n_classes = in_channels,feature_out_channels,n_classes\n        self.hidden_features = hidden_features\n        self.n_planes,self.n_series = n_planes,n_series\n        self.residual_config_list = residual_config_list\n        \n        self.features = NetFeatures(self.in_channels, self.feature_out_channels, self.residual_config_list)\n        self.pool_flat_linear = nn.Sequential(nn.AdaptiveAvgPool3d(1),\n            nn.Flatten(),\n            nn.Linear(self.features.out_channels, self.hidden_features),\n        )\n        self.emb_series = ConcatEmbeddingLinear(self.hidden_features, self.hidden_features, self.n_series)\n        self.emb_planes = ConcatEmbeddingLinear(self.hidden_features, self.hidden_features, self.n_planes)\n        self.classifier = nn.Sequential(\n            nn.Dropout(0.2),\n            nn.Linear(self.hidden_features, self.n_classes)\n        )\n\n    def forward(self, x, idx_series, idx_planes):\n        x = self.features(x)\n        x = self.pool_flat_linear(x)\n        x = self.emb_series(x, idx_series)\n        x = self.emb_planes(x, idx_planes)\n        out = self.classifier(x)\n        return out\n\n\ndef get_residual_config_backup():\n    # Like MobileNetV3 small, although it may be too deep.\n    # in_channels, out_channels, expand_channels, kernel_size, stride, padding, norm, activation, use_se\n    conf = list()\n    conf.append(BottleNeck.config(16, 16, 16, 3, 2, 1, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(16, 24, 72, 3, 2, 1, nn.BatchNorm3d, nn.ReLU, False))\n    conf.append(BottleNeck.config(24, 24, 88, 3, 1, 1, nn.BatchNorm3d, nn.ReLU, False))\n    conf.append(BottleNeck.config(24, 40, 96, 5, 2, 2, nn.BatchNorm3d, nn.ReLU, True))\n    conf.append(BottleNeck.config(40, 40, 240, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(40, 40, 240, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(40, 48, 120, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(48, 48, 144, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(48, 96, 288, 5, 2, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(96, 96, 576, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(96, 96, 576, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    return conf\n\n\ndef get_residual_config():\n    # in_channels, out_channels, expand_channels, kernel_size, stride, padding, norm, activation, use_se\n    conf = list()\n    conf.append(BottleNeck.config(16, 16, 16, 3, 2, 1, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(16, 24, 72, 3, 2, 1, nn.BatchNorm3d, nn.ReLU, False))\n    conf.append(BottleNeck.config(24, 24, 88, 3, 1, 1, nn.BatchNorm3d, nn.ReLU, False))\n    conf.append(BottleNeck.config(24, 40, 96, 5, 2, 2, nn.BatchNorm3d, nn.ReLU, True))\n    conf.append(BottleNeck.config(40, 40, 240, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(40, 80, 288, 5, 2, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    conf.append(BottleNeck.config(80, 96, 576, 5, 1, 2, nn.BatchNorm3d, nn.Hardswish, True))\n    return conf\n\n\nif DEBUG:\n    t_ = torch.ones(4, 1, 64, 256, 256, dtype=torch.float32)\n    l_ = torch.ones(4, 1, dtype=torch.float32)\n    s_ = torch.ones(4, dtype=torch.int64)\n    p_ = torch.ones(4, dtype=torch.int64)\n    config_ = get_residual_config()\n    net_ = Net(1, 512, 512, 1, 4, 4, config_).to(dtype=torch.float32)\n    print(net_)\n    with torch.no_grad():\n        o_ = net_(t_, s_, p_)\n        loss_ = F.binary_cross_entropy_with_logits(o_, l_)\n        print(loss_.item())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Voxel MRI augmentations","metadata":{}},{"cell_type":"code","source":"class RandomInvert3D(AugmentationBase3D):\n    \n    def __init__(\n        self,\n        max_val: Union[float, torch.Tensor] = torch.tensor(1.0),\n        return_transform: bool = False,\n        same_on_batch: bool = False,\n        p: float = 0.5,\n    ) -> None:\n        super(RandomInvert3D, self).__init__(\n            p=p, return_transform=return_transform, same_on_batch=same_on_batch, p_batch=1.0\n        )\n        self.max_val = max_val\n\n    def __repr__(self) -> str:\n        return self.__class__.__name__ + f\"({super().__repr__()})\"\n    \n    def generate_parameters(self, batch_shape: torch.Size):\n        return dict(max_val=torch.as_tensor(self.max_val), batch_shape=torch.as_tensor(batch_shape))\n    \n    def compute_transformation(self, input, params: Dict[str, torch.Tensor]):\n        return self.identity_matrix(input)\n\n    def apply_transform(\n        self, input: torch.Tensor,\n        params: Dict[str, torch.Tensor],\n        transform: Optional[torch.Tensor] = None\n    ) -> torch.Tensor:\n        max_val = params[\"max_val\"]\n        return invert(input, max_val)\n\n    \nNumeric = Union[int, float]\n\n\nclass RandomShift3D(nn.Module):\n    \n    def __init__(self,\n                 shift_limit: Union[Numeric, List[Numeric], Tuple[Numeric, Numeric]] = 0.125,\n                 p: float = 0.5):\n        super(RandomShift3D, self).__init__()\n        self.shift_limit,self.p = shift_limit,p\n        if isinstance(self.shift_limit, (float, int)):\n            self.shift_limit = np.array(((-abs(self.shift_limit), abs(self.shift_limit)),\n                                         (-abs(self.shift_limit), abs(self.shift_limit)),\n                                         (-abs(self.shift_limit), abs(self.shift_limit)),), dtype=np.float64)\n        elif isinstance(self.shift_limit, (tuple, list)):\n            self.shift_limit = np.array(self.shift_limit, dtype=np.float64)\n        else:\n            raise TypeError(\"shift_limit expects \")\n        self.shift_limit = np.clip(self.shift_limit, -1., 1.)\n        if self.shift_limit.shape[0] == 1:\n            self.shift_limit = np.concatenate([self.shift_limit, self.shift_limit, self.shift_limit])\n        assert self.shift_limit.shape == (3, 2), f\"\"\n    \n    def forward(self, tensor):\n        assert len(tensor.shape) == 5, f\"Requires 5 dims torch.Tensor[N, C, D, H, W], got {tensor.shape}\"\n        n, c, d, h, w = tensor.shape\n        apply_proba = np.random.uniform(size=(n,))\n        shift_size = np.random.uniform(low=self.shift_limit[:, 0], high=self.shift_limit[:, 1], size=(n, 3))\n        shift_d, shift_h, shift_w = (np.array(tensor.shape[2:])[np.newaxis, :] * shift_size).astype(np.int64).T\n        out = torch.zeros_like(tensor)\n        for i in range(n):\n            if apply_proba[i] <= self.p:\n                out[i, :,\n                    max(0, 0+shift_d[i]):min(d, d+shift_d[i]),\n                    max(0, 0+shift_h[i]):min(h, h+shift_h[i]),\n                    max(0, 0+shift_w[i]):min(w, w+shift_w[i]),\n                ] = tensor[i, :,\n                    max(0, 0-shift_d[i]):min(d, d-shift_d[i]),\n                    max(0, 0-shift_h[i]):min(h, h-shift_h[i]),\n                    max(0, 0-shift_w[i]):min(w, w-shift_w[i]),\n                ]\n            else:\n                out[i] = tensor[i]  # Unchanged.\n        return out\n\n\ndef get_augmentation(split=\"train\") -> nn.Sequential:\n    \"\"\"\n    Get Sequence of augmentations.\n    :return: nn.Sequential: requires input: torch.FloatTensor[N, C, D, H, W] in range[0., 1.]\n    \"\"\"\n    if split in (\"test\", \"val\"):\n        aug_list = nn.Sequential()\n    elif split == \"train\":\n        aug_list = nn.Sequential(\n            K.augmentation.RandomAffine3D(degrees=(5., 5., 90.), translate=(.05, .05, .05), scale=(.98, 1.02), p=.3),\n            K.augmentation.RandomHorizontalFlip3D(p=.3),\n#             K.augmentation.RandomVerticalFlip3D(p=.1),\n#             K.augmentation.RandomRotation3D((0., 0., 90.), p=1.0)\n            RandomShift3D(shift_limit=0.2, p=.3),\n            RandomInvert3D(p=.1),\n        )\n    else:\n        raise ValueError(f\"Argument `split` must in {{'train', 'val', 'test'}}, got {split}\")\n    aug_list.requires_grad_(False)\n    return aug_list\n\n\ndef plot_grid(t: torch.tensor) -> None:\n    \"\"\"\n    Plot image by middle index\n    :argument: t: torch.Tensor[N, C, D, H, W]\n    \"\"\"\n    from itertools import product\n    a = int(np.ceil(np.sqrt(len(t))))\n    fig, axes = plt.subplots(a, a, figsize=(14, 14))\n    for nth, (i, j) in zip(range(len(t)), product(range(a), range(a))):\n        nth_img = t[nth].squeeze(0).numpy()\n        nth_img_mid = nth_img[len(nth_img) // 2]\n        mean, std = np.mean(nth_img), np.std(nth_img)\n        axes[i, j].imshow(nth_img_mid, cmap=plt.cm.gray)\n        axes[i, j].set_title(f\"mean: {mean:.4f}, std: {std:.4f}\")\n        axes[i, j].set_axis_off()\n    plt.show()\n\n\n# Check the effect of augmentation.\nds_ = MRIVoxelDataset(df_train, labels_train, (64, 128))\ndl_ = DataLoader(ds_, batch_size=16, shuffle=True, num_workers=6)\naug_ = get_augmentation(split=\"train\")\nfor voxel, label, (_, _) in dl_:\n    voxel = aug_(voxel)\n    plot_grid(voxel)\n    break","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training and evaluating functions","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import roc_auc_score\n\n\n@torch.no_grad()\ndef evaluation(loader: DataLoader, model, aug_list, device):\n    if isinstance(loader.dataset, Subset):\n        df_copy: pd.DataFrame = (loader.dataset.dataset.label_df.iloc[loader.dataset.indices].copy())\n    else:\n        df_copy: pd.DataFrame = loader.dataset.label_df.copy()\n    df_copy[\"MGMT_value_pred\"] = np.full((len(df_copy),), fill_value=np.nan, dtype=np.float64)\n    df_copy = df_copy.reset_index(drop=True)\n    batch_size = loader.batch_size\n    model.to(device)\n    model.eval()\n    loss_val = 0.\n    for n, (voxel, label, (series_desc, plane)) in tqdm(enumerate(loader), desc=\"Evaluating with AUC\", total=len(loader)):\n        voxel, label, series_desc, plane = (aug_list(voxel).to(device), label.to(device),\n                                            series_desc.to(device), plane.to(device))\n        out = model(voxel, series_desc, plane)\n        loss = nn.functional.binary_cross_entropy_with_logits(out, label)\n        loss_val += loss.item()\n        pred_proba = torch.sigmoid(out.detach())[:, 0].cpu().numpy()\n        df_copy.iloc[n*batch_size:n*batch_size+len(voxel), df_copy.columns.get_loc(\"MGMT_value_pred\")] = pred_proba\n    loss_val /= max(1, len(loader))\n    auc = roc_auc_score(df_copy[\"MGMT_value\"].values.astype(\"int64\"), df_copy[\"MGMT_value_pred\"].values)\n    return loss_val, auc\n\n\ndef train_an_epoch(loader, model, aug_list, device, optimizer, epoch_idx=0):\n    model.train()\n    model.to(device)\n    loss_epoch = 0.\n    for voxel, label, (series_desc, plane) in tqdm(loader, desc=f\"Training epoch: {epoch_idx}\"):\n        voxel, label = aug_list(voxel).to(device), label.to(device)\n        series_desc, plane = series_desc.to(device), plane.to(device)\n        out = model(voxel, series_desc, plane)\n        loss = nn.functional.binary_cross_entropy_with_logits(out, label)\n        if np.isnan(loss.item()):\n            raise ValueError(\"loss is nan\")\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        loss_epoch += loss.item()\n    loss_epoch /= max(1, len(loader))\n    print(f\"Epoch: {e}, train loss: {loss_epoch}\")\n\n\nif DEBUG:\n    ds_ = MRIVoxelDataset(df_train, labels_train, (64, 128))\n    sub_ = Subset(ds_, list(range(len(ds_)))[:101])\n    dl_ = DataLoader(sub_, batch_size=4, num_workers=4)\n    aug_list_ = get_augmentation(split=\"val\")\n    conf_ = get_residual_config()\n    net_ = Net(1, 512, 512, 1, 4, 4, conf_).to(torch.float32)\n    optim_ = torch.optim.Adam(net_.parameters(), lr=3e-4)\n    for e in range(1):\n        train_an_epoch(dl_, net_, aug_list_, device, optim_, e)\n        loss_val, auc = evaluation(dl_, net_, aug_list_, device)\n        print(loss_val, auc)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train the model with data pipeline","metadata":{}},{"cell_type":"code","source":"voxel_size = (64, 64, 64)\nincluding_series = np.array([\n    mri_series_map[\"FLAIR\"],\n    mri_series_map[\"T1w\"],\n    mri_series_map[\"T1wCE\"],\n    mri_series_map[\"T2w\"],\n], dtype=np.int64)\nNumpyNDArray = Iterable\n\ndef get_dataset_in_pipeline(img_dir: Path, including_series: NumpyNDArray[np.int64],\n                            voxel_size: Union[int, Tuple[int, int], Tuple[int, int, int]] = (64, 256, 256),\n                            glob: str = None,\n                            df_labels: pd.DataFrame = None, drop_ids: List[str] = None):\n    df_meta = get_meta_from_glob(img_dir, glob)\n    df_meta.loc[:, \"plane\"] = df_meta.loc[:, \"plane\"].map(planes_map)\n    df_meta.loc[:, \"series_description\"] = df_meta.loc[:, \"series_description\"].map(mri_series_map)\n    df_meta = keep_non_blank(df_meta)\n    if df_labels is not None:\n        df_labels = drop_by_id(df_labels, drop_ids)\n        df_meta = drop_by_id(df_meta, drop_ids)\n    ds = MRIVoxelDataset(df_meta, df_labels, voxel_size, including_series)\n    return ds\n\n\nlabels_train = pd.read_csv(data_dir / \"train_labels.csv\", dtype={\"BraTS21ID\": str})\nds_train = get_dataset_in_pipeline(data_dir / \"train\", including_series, voxel_size,\n                                   None, labels_train, \"00109, 00123, 00709\".split(\", \"))\nds_test = get_dataset_in_pipeline(data_dir / \"test\", including_series, voxel_size)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train model by cross validation","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import StratifiedKFold\n\n# Parameters to construct Net\nin_channels = 1\nfeature_out_channels = 576\nhidden_features = 512\nn_classes = 1\nn_series = len(mri_series)\nn_planes = len(planes)\nresidual_config = get_residual_config()\n\n# Training Parameters\nn_splits = 5\nbatch_size = 16\nepochs = 15\nlr = 3e-4\nnum_workers = 6\nweight_decay = 1e-5\n\nsave_filename_template = \"{}-model-fold_{:03d}-best-state_dict.pt\"\nif USE_CROSS_VALIDATION and (not INFERENCE_ONLY):\n    splitter = StratifiedKFold(n_splits, shuffle=True, random_state=42)\n    df_label_train = ds_train.label_df\n    for nth_fold, (idx_trn, idx_val) in enumerate(splitter.split(df_label_train, df_label_train[\"series_description\"].values)):\n        ds_trn,ds_val = (Subset(ds_train, idx_trn), Subset(ds_train, idx_val))\n        dl_trn,dl_val = (DataLoader(ds_trn, batch_size=batch_size, shuffle=True, num_workers=num_workers),\n                         DataLoader(ds_val, batch_size=batch_size, num_workers=num_workers))\n        aug_trn,aug_val = get_augmentation(split=\"train\"),get_augmentation(split=\"val\")\n        model = Net(\n            in_channels,\n            feature_out_channels,\n            hidden_features,\n            n_classes,\n            n_series,\n            n_planes,\n            residual_config,\n        ).to(torch.float32)\n        optim = torch.optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)\n        scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(\n            optim, \"min\", factor=0.5, patience=3, cooldown=0, verbose=True\n        )\n        best_val_loss = np.inf\n        for e in range(epochs):\n            train_an_epoch(dl_trn, model, aug_trn, device, optim, e)\n            loss_val, auc_val = evaluation(dl_val, model, aug_val, device)\n            print(f\"fold {nth_fold}, epoch {e}, val loss: {loss_val:.6f}, val AUC: {auc_val:.6f}\")\n            if loss_val < best_val_loss:\n                print(f\"Best val loss changed from {best_val_loss:.6f} to {loss_val:.6f}\")\n                best_val_loss = loss_val\n                torch.save(model.state_dict(), save_filename_template.format(model_name, nth_fold))\n            scheduler.step(loss_val)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train model by full data\n\nTrain with augmented data, as well as on validation. (Overfit the random Augmented validation set)","metadata":{}},{"cell_type":"code","source":"# Parameters to construct Net\nin_channels = 1\nfeature_out_channels = 576\nhidden_features = 512\nn_classes = 1\nn_series = len(mri_series)\nn_planes = len(planes)\nresidual_config = get_residual_config()\n\n# Training Parameters\nbatch_size = 16\nepochs = 18\nlr = 3e-4\nnum_workers = 6\nweight_decay = 1e-5\n\nsave_filename_template = \"{}-model-whole-dataset-best-state_dict.pt\"\n\nif not INFERENCE_ONLY:\n    dl_trn,dl_val = (DataLoader(ds_train, batch_size=batch_size, shuffle=True, num_workers=num_workers),\n                     DataLoader(ds_train, batch_size=batch_size, num_workers=num_workers))\n    aug_trn,aug_val = get_augmentation(split=\"train\"),get_augmentation(split=\"train\")\n    model = Net(\n        in_channels,\n        feature_out_channels,\n        hidden_features,\n        n_classes,\n        n_series,\n        n_planes,\n        residual_config,\n    ).to(torch.float32)\n    optim = torch.optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)\n    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(\n        optim, \"min\", factor=0.3, patience=3, cooldown=0, verbose=True\n    )\n    best_val_loss = np.inf\n    for e in range(epochs):\n        train_an_epoch(dl_trn, model, aug_trn, device, optim, e)\n        loss_val, auc_val = evaluation(dl_val, model, aug_val, device)\n        print(f\"epoch {e}, val loss: {loss_val:.6f}, val AUC: {auc_val:.6f}\")\n        if loss_val < best_val_loss:\n            print(f\"Best val loss changed from {best_val_loss:.6f} to {loss_val:.6f}\")\n            best_val_loss = loss_val\n            torch.save(model.state_dict(), save_filename_template.format(model_name))\n        scheduler.step(loss_val)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Summary in train set","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\n\n\ndef load_model(path, *net_args, **net_kwargs):\n    net = Net(*net_args, **net_kwargs)\n    state_dict = torch.load(path)\n    net.load_state_dict(state_dict)\n    return net\n\n\n@torch.no_grad()\ndef train_performance(loader: DataLoader, models_list: List[Net], aug_list, device, n_pics=10):\n    \n    if isinstance(loader.dataset, Subset):\n        df_copy: pd.DataFrame = loader.dataset.dataset.label_df.iloc[loader.dataset.indices].copy()\n        meta_df = loader.dataset.dataset.meta_df.copy()\n        voxel_size = loader.dataset.dataset.voxel_size\n    else:\n        df_copy: pd.DataFrame = loader.dataset.label_df.copy()\n        meta_df = loader.dataset.meta_df.copy()\n        voxel_size = loader.dataset.voxel_size\n    print(df_copy[\"MGMT_value\"].value_counts())\n    batch_size = loader.batch_size\n    for i, model in enumerate(models_list):\n        df_copy.loc[:, f\"MGMT_value_{i}\"] = np.zeros_like(df_copy[\"BraTS21ID\"], dtype=np.float64)\n        model.to(device)\n        model.eval()\n        for n, (voxel, _, (series_desc, plane)) in tqdm(enumerate(loader),\n                                                        desc=f\"Inferencing with model idx {i}\", total=len(loader)):\n            voxel, series_desc, plane = aug_list(voxel).to(device), series_desc.to(device), plane.to(device)\n            out = model(voxel, series_desc, plane)\n            pred_proba = torch.sigmoid(out.detach())[:, 0].cpu().numpy()\n            df_copy.iloc[n*batch_size:n*batch_size+len(voxel), df_copy.columns.get_loc(f\"MGMT_value_{i}\")] = pred_proba\n    use_cols = [s for s in df_copy.columns if s.startswith(\"MGMT_value_\")]\n    df_copy[\"MGMT_value_pred\"] = df_copy.loc[:, use_cols].mean(axis=1)\n    \n    \n    auc = roc_auc_score(df_copy[\"MGMT_value\"].values.astype(\"int64\"), df_copy[\"MGMT_value_pred\"].values)\n    print(f\"AUC: {auc}\")\n\n    cm = confusion_matrix(df_copy[\"MGMT_value\"].values.astype(\"int64\"),\n                          df_copy[\"MGMT_value_pred\"].map(lambda v: 1 if v>=0.5 else 0).values.astype(\"int64\"))\n    fig, ax = plt.subplots(figsize=(8, 6))\n    sns.heatmap(cm, annot=True, fmt=\"d\", ax=ax)\n    ax.set_xlabel(\"Pred\")\n    ax.set_ylabel(\"True\")\n    plt.show()\n    \n    tp = (df_copy.loc[(df_copy[\"MGMT_value\"] == 1) & (df_copy[\"MGMT_value_pred\"] >= 0.5)]\n          .sort_values(\"MGMT_value_pred\", ascending=False))\n    fp = (df_copy.loc[(df_copy[\"MGMT_value\"] == 0) & (df_copy[\"MGMT_value_pred\"] >= 0.5)]\n          .sort_values(\"MGMT_value_pred\", ascending=False))\n    tn = (df_copy.loc[(df_copy[\"MGMT_value\"] == 0) & (df_copy[\"MGMT_value_pred\"] < 0.5)]\n          .sort_values(\"MGMT_value_pred\", ascending=True))\n    fn = (df_copy.loc[(df_copy[\"MGMT_value\"] == 1) & (df_copy[\"MGMT_value_pred\"] < 0.5)]\n          .sort_values(\"MGMT_value_pred\", ascending=True))\n    titles = [\"TP\", \"FP\", \"TN\", \"FN\"]\n    fig, axes = plt.subplots(n_pics, 4, figsize=(16, int(5*n_pics)))\n    for j, df in enumerate([tp, fp, tn, fn]):\n        for n, (_, row) in zip(range(n_pics), df.iterrows()):\n            voxel, plane = get_voxel_by_id_series(meta_df, row[\"BraTS21ID\"], row[\"series_description\"], voxel_size[1:])\n            img = voxel[(len(voxel)-1) // 2]\n            axes[n, j].set_title(f\"{titles[j]}, likelihood of positive: {row['MGMT_value_pred']:.4f}\\n\"\n                                 f\"type: {row['series_description']}, plane: {plane}.\")\n            axes[n, j].imshow(img, cmap=plt.cm.gray)\n            axes[n, j].set_axis_off()\n    plt.show()\n\n\naug_list = get_augmentation(split=\"val\")\nmodels_path = list(sorted(models_dir.glob(f\"{model_name}*whole*best-state_dict.pt\")))\nif USE_CROSS_VALIDATION:\n    models_path.extend(list(sorted(models_dir.glob(f\"{model_name}*fold*best-state_dict.pt\"))))\nprint(models_path)\nmodels_list = [load_model(path,\n                          in_channels,\n                          feature_out_channels,\n                          hidden_features,\n                          n_classes,\n                          n_series,\n                          n_planes,\n                          residual_config,\n) for path in models_path]\ndl_train = DataLoader(ds_train, batch_size=batch_size, num_workers=num_workers)\ntrain_performance(dl_train, models_list, aug_list, device)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@torch.no_grad()\ndef inference_by_models(loader: DataLoader, models_list: List[Net], aug_list, device):\n\n    df_copy: pd.DataFrame = loader.dataset.label_df.copy()\n    batch_size = loader.batch_size\n    for i, model in enumerate(models_list):\n        df_copy.loc[:, f\"MGMT_value_{i}\"] = np.full_like(df_copy[\"BraTS21ID\"], np.nan, dtype=np.float64)\n        model.to(device)\n        model.eval()\n        for n, (voxel, _, (series_desc, plane))  in tqdm(enumerate(loader),\n                                                         desc=f\"Inferencing with model idx {i}\", total=len(loader)):\n            voxel, series_desc, plane = aug_list(voxel).to(device), series_desc.to(device), plane.to(device)\n            out = model(voxel, series_desc, plane)\n            pred_proba = torch.sigmoid(out.detach())[:, 0]\n            df_copy.iloc[n*batch_size:n*batch_size+len(voxel),\n                         df_copy.columns.get_loc(f\"MGMT_value_{i}\")] = pred_proba.cpu().numpy()\n    df = df_copy.groupby(\"BraTS21ID\").mean()\n    use_cols = [s for s in df.columns if s.startswith(\"MGMT_value_\")]\n    df[\"MGMT_value\"] = df.loc[:, use_cols].mean(axis=1)\n    df = df.reset_index()\n    submission = df.loc[:, [\"BraTS21ID\", \"MGMT_value\"]].copy()\n    return submission\n\n\naug_list = get_augmentation(split=\"test\")\nmodels_path = list(sorted(models_dir.glob(f\"{model_name}*whole*best-state_dict.pt\")))\nif USE_CROSS_VALIDATION:\n    models_path.extend(list(sorted(models_dir.glob(f\"{model_name}*fold*best-state_dict.pt\"))))\nprint(models_path)\nmodels_list = [load_model(path,\n                          in_channels,\n                          feature_out_channels,\n                          hidden_features,\n                          n_classes,\n                          n_series,\n                          n_planes,\n                          residual_config,\n) for path in models_path]\ndl_test = DataLoader(ds_test, batch_size=batch_size, num_workers=num_workers)\nsubmission = inference_by_models(dl_test, models_list, aug_list, device)\nsubmission.to_csv(\"submission.csv\", index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.displot(submission[\"MGMT_value\"], kde=True, bins=50)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"time_end = time.time()\nduration = time_end - time_begin\nprint(f\"Kernel duration: {duration:.4f} seconds ({(duration / 60.):.4f} Min, {(duration / 3600.):.4f} Hour).\")","metadata":{},"execution_count":null,"outputs":[]}]}