{"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":"!pip install ../input/rsnamiccai-btrc-dataset/packages/monai-0.7.0-202109240007-py3-none-any.whl","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:34:40.721306Z","iopub.execute_input":"2021-10-15T07:34:40.722399Z","iopub.status.idle":"2021-10-15T07:35:09.095287Z","shell.execute_reply.started":"2021-10-15T07:34:40.722357Z","shell.execute_reply":"2021-10-15T07:35:09.093779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nfrom glob import glob\nfrom tqdm.notebook import tqdm\n\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport scipy.ndimage\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.metrics import roc_auc_score, roc_curve, auc\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader, RandomSampler, SequentialSampler\nfrom monai.networks.nets import SegResNet, DenseNet","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-10-15T08:23:36.389091Z","iopub.execute_input":"2021-10-15T08:23:36.389862Z","iopub.status.idle":"2021-10-15T08:23:36.397339Z","shell.execute_reply.started":"2021-10-15T08:23:36.389831Z","shell.execute_reply":"2021-10-15T08:23:36.396087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"RAW_DATA_PATH = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\nPROCESSED_DATA_PATH = '../input/rsnamiccai-btrc-dataset'\nMODELS_PATH = '../input/rsnamiccai-btrc-dataset'\nPREDICTIONS_PATH = '../input/rsnamiccai-btrc-dataset/predictions'","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:37:29.131442Z","iopub.execute_input":"2021-10-15T07:37:29.132250Z","iopub.status.idle":"2021-10-15T07:37:29.140315Z","shell.execute_reply.started":"2021-10-15T07:37:29.132217Z","shell.execute_reply":"2021-10-15T07:37:29.139120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_test_dtypes = {\n    'BraTS21ID': 'object',\n    'MGMT_value': np.uint8,\n    'fold': np.uint8\n}\n\ndf_train = pd.read_csv(f'{PROCESSED_DATA_PATH}/train_task2.csv', dtype=train_test_dtypes)\ndf_test = pd.read_csv(f'{RAW_DATA_PATH}/sample_submission.csv', usecols=['BraTS21ID'], dtype=train_test_dtypes)\n\nprint(f'Training Set Shape: {df_train.shape} - Memory Usage: {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set Shape: {df_test.shape} - Memory Usage: {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:50:58.660003Z","iopub.execute_input":"2021-10-15T07:50:58.660706Z","iopub.status.idle":"2021-10-15T07:50:58.684897Z","shell.execute_reply.started":"2021-10-15T07:50:58.660674Z","shell.execute_reply":"2021-10-15T07:50:58.683892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## DICOM Preprocessing","metadata":{}},{"cell_type":"code","source":"def apply_manual_voi_lut(dicom_file, window_width, window_center):\n\n    \"\"\"\n    Create a VOI LUT with given window width and window center and apply it to dicom file's pixel array\n\n    Parameters\n    ----------\n    dicom_file (pydicom.dataset.FileDataset): Dicom file read into memory\n    window_width (int): Width of the modality pixel values\n    window_center (int): Center of the modality pixel values\n\n    Returns\n    -------\n    image [array-like of shape (width, height)]: Array of 2D image after manual VOI LUT applied to pixel array\n    \"\"\"\n\n    min_pixel_value = int(np.amin(dicom_file.pixel_array))\n    max_pixel_value = int(np.amax(dicom_file.pixel_array))\n\n    # Make an empty array for the LUT the size of the pixel 'width' in the raw pixel data\n    voi_lut = [0] * (max_pixel_value + 1)\n\n    # Invert pixel values and window level for MONOCHROME1 photometric interpretation\n    invert = False\n    if dicom_file.PhotometricInterpretation == 'MONOCHROME1':\n        invert = True\n    else:\n        window_center = (max_pixel_value - min_pixel_value) - window_center\n\n    # Loop through the pixels and calculate each LUT value\n    for pixel_value in range(min_pixel_value, max_pixel_value):\n        modality_lut_value = pixel_value * float(dicom_file.RescaleSlope) + float(dicom_file.RescaleIntercept)\n        voi_lut_value = (((modality_lut_value - window_center) / window_width + 0.5) * 255.0)\n        clamped_value = min(max(voi_lut_value, 0), 255)\n\n        if invert:\n            voi_lut[pixel_value] = round(255 - clamped_value)\n        else:\n            voi_lut[pixel_value] = round(clamped_value)\n\n    voi_lut = np.array(voi_lut)\n    return np.uint8(voi_lut[dicom_file.pixel_array])\n\n\ndef apply_auto_voi_lut(dicom_file):\n\n    \"\"\"\n    Apply VOI LUT if it exists in the dicom file, otherwise use window width and window center given in the dicom file\n\n    Parameters\n    ----------\n    dicom_file (pydicom.dataset.FileDataset): Dicom file read into memory\n\n    Returns\n    -------\n    image [array-like of shape (width, height)]: Array of 2D image after automatic VOI LUT applied to pixel array\n    \"\"\"\n\n    image = apply_voi_lut(dicom_file.pixel_array, dicom_file)\n\n    if dicom_file.PhotometricInterpretation == 'MONOCHROME1':\n        image = np.amax(image) - image\n\n    image = image - np.min(image)\n    image = image / np.max(image)\n    image = (image * 255).astype(np.uint8)\n\n    return image\n\n\ndef get_plane(dicom_file):\n\n    \"\"\"\n    Extract image plane from ImageOrientationPatient field of the dicom file\n\n    Parameters\n    ----------\n    dicom_file (pydicom.dataset.FileDataset): Dicom file read into memory\n\n    Returns\n    -------\n    plane (str): Image plane (Coronal, Sagittal or Axial)\n    \"\"\"\n\n    image_orientation_patient = dicom_file[0x0020, 0x0037]\n\n    row_x = round(image_orientation_patient[0])\n    row_y = round(image_orientation_patient[1])\n    col_x = round(image_orientation_patient[3])\n    col_y = round(image_orientation_patient[4])\n\n    plane = None\n\n    if row_x == 1 and row_y == 0 and col_x == 0 and col_y == 0:\n        plane = 'Coronal'\n\n    if row_x == 0 and row_y == 1 and col_x == 0 and col_y == 0:\n        plane = 'Sagittal'\n\n    if row_x == 1 and row_y == 0 and col_x == 0 and col_y == 1:\n        plane = 'Axial'\n\n    return plane\n\n\ndef change_spacing(mri, current_spacing, new_spacing):\n\n    \"\"\"\n    Change spacing of height, width and depth\n\n    Parameters\n    ----------\n    mri [array-like of shape (depth, width, height)]: Array of 3D mpMRI\n    current_spacing [array-like of shape (3)]: Array of current spacings in Z, X, Y direction (millimeters)\n    new_spacing [array-like of shape (3)]: Array of new spacings in Z, X, Y direction (millimeters)\n\n    Returns\n    -------\n    mri [array-like of shape (depth, width, height)]: Array of 3D mpMRI after change spacings\n    \"\"\"\n\n    resize_factor = current_spacing / new_spacing\n    normalized_shape = np.round(mri.shape * resize_factor)\n    resize_factor = normalized_shape / mri.shape\n    mri = scipy.ndimage.interpolation.zoom(mri, resize_factor, mode='nearest')\n\n    return mri\n\n\ndef load_mri(mri_path, window_width=None, window_center=None, voi_lut=None, new_spacing=None, reorder_plane=False, resize_shape=None, verbose=False):\n\n    \"\"\"\n    Read slices of mpMRI into memory and apply preprocessing steps\n\n    Parameters\n    ----------\n    mri_path (str): Directory of the mpMRI\n    window_width (int or None): Width of the modality pixel values\n    window_center (int or None): Center of the modality pixel values\n    voi_lut (str or None): Whether to use manual or auto VOI LUT (\"manual\" or \"auto\")\n    new_spacing [array-like of shape (3)] or None: Array of new spacings in Z, X, Y direction (millimeters)\n    reorder_plane (bool): Whether reorder planes or not\n    resize_shape (int or None): Resize shape of the planes\n    verbose (bool): Verbosity flag\n\n    Returns\n    -------\n    mri [np.ndarray of shape (depth, width, height)]: Array of 3D mpMRI\n    \"\"\"\n\n    slice_paths = sorted(glob(f'{mri_path}/*.dcm'), key=lambda x: int(str(x).split('-')[-1].split('.')[0]))\n    dicom_files = [pydicom.dcmread(slice_path) for slice_path in slice_paths]\n    slices = []\n\n    for i, dicom_file in enumerate(dicom_files):\n\n        if voi_lut == 'manual':\n            # Applying manually created voi lut\n            image = apply_manual_voi_lut(dicom_file=dicom_file, window_width=window_width, window_center=window_center)\n        elif voi_lut == 'auto':\n            # Applying voi lut of the dicom file\n            image = apply_auto_voi_lut(dicom_file=dicom_file)\n        else:\n            # Not applying voi lut\n            image = dicom_file.pixel_array\n\n            if dicom_file.PhotometricInterpretation == 'MONOCHROME1':\n                image = np.amax(image) - image\n\n            # Exclude empty slices\n            if np.all(image == np.min(image)):\n                continue\n\n        slices.append(image)\n\n    # Not processing all zero mpMRIs\n    if len(slices) == 0:\n        return None\n\n    mri = np.stack(slices)\n\n    # Change spacing if new spacing is given\n    current_spacing = np.array([float(dicom_files[0].SliceThickness)] + list(dicom_files[0].PixelSpacing), dtype=np.float32)\n    if new_spacing is not None:\n        if np.any(np.array(new_spacing) != current_spacing):\n            mri = change_spacing(mri=mri, current_spacing=current_spacing, new_spacing=new_spacing)\n\n    positions = [dicom_file.ImagePositionPatient for dicom_file in dicom_files]\n    plane = get_plane(dicom_file=dicom_files[0])\n\n    # Reorder plane if it is set to True\n    if reorder_plane:\n        if plane == 'Coronal':\n            if positions[0][1] < positions[-1][1]:\n                mri = mri[::-1]\n            mri = mri.transpose((1, 0, 2))\n        elif plane == 'Sagittal':\n            if positions[0][0] < positions[-1][0]:\n                mri = mri[::-1]\n            mri = mri.transpose((1, 2, 0))\n            mri = np.rot90(mri, 2, axes=(1, 2))\n        elif plane == 'Axial':\n            if positions[0][2] > positions[-1][2]:\n                mri = mri[::-1]\n            mri = np.rot90(mri, 2)\n\n    # Crop non-zero slices along Z-X, Z-Y and X-Y axes\n    mmin = np.array((mri > 0).nonzero()).min(axis=1)\n    mmax = np.array((mri > 0).nonzero()).max(axis=1)\n    mri = mri[\n        mmin[0]:mmax[0] + 1,\n        mmin[1]:mmax[1] + 1,\n        mmin[2]:mmax[2] + 1,\n    ]\n\n    # Resize sampled planes from longest axis if resize shape is given\n    if resize_shape is not None:\n\n        resized_mri = np.zeros((resize_shape, resize_shape, resize_shape), dtype=np.int16)\n\n        if np.argmax(mri.shape) == 0:\n            for i, s in enumerate(np.linspace(0, mri.shape[0] - 1, resize_shape)):\n                resized_mri[i] = cv2.resize(mri[int(s)], (resize_shape, resize_shape), interpolation=cv2.INTER_LANCZOS4)\n        elif np.argmax(mri.shape) == 1:\n            for i, s in enumerate(np.linspace(0, mri.shape[1] - 1, resize_shape)):\n                resized_mri[:, i] = cv2.resize(mri[:, int(s)], (resize_shape, resize_shape), interpolation=cv2.INTER_LANCZOS4)\n        elif np.argmax(mri.shape) == 2:\n            for i, s in enumerate(np.linspace(0, mri.shape[2] - 1, resize_shape)):\n                resized_mri[:, :, i] = cv2.resize(mri[:, :, int(s)], (resize_shape, resize_shape), interpolation=cv2.INTER_LANCZOS4)\n\n        mri = resized_mri\n\n    if verbose:\n        print(f'{mri_path} - MRI Shape: {mri.shape} - Mean: {np.mean(mri):.2f} - Std: {np.std(mri):.2f} - Min: {np.min(mri):.2f} - Max: {np.max(mri):.2f} - Type: {mri.dtype}')\n\n    return mri\n","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:37:39.493905Z","iopub.execute_input":"2021-10-15T07:37:39.494240Z","iopub.status.idle":"2021-10-15T07:37:39.535277Z","shell.execute_reply.started":"2021-10-15T07:37:39.494189Z","shell.execute_reply":"2021-10-15T07:37:39.534008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Tumor Segmentation","metadata":{}},{"cell_type":"code","source":"class SegmentationFeatureExtractor:\n\n    def __init__(self, df_train, df_test, train_features_path=None):\n\n        self.df_train = df_train\n        self.df_test = df_test\n        self.train_features_path = train_features_path\n        self.device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')\n\n    def load_segmentation_models(self):\n\n        models = {}\n\n        for mri_type in ['FLAIR', 'T1w']:\n            model = SegResNet(\n                spatial_dims=3,\n                in_channels=1,\n                out_channels=3,\n                init_filters=8,\n                dropout_prob=0.0,\n                blocks_down=(1, 2, 2, 4),\n                blocks_up=(1, 1, 1),\n                upsample_mode='nontrainable'\n            )\n            model.load_state_dict(torch.load(f'{MODELS_PATH}/segresnet/segresnet_{mri_type}.pt'))\n            model.to(self.device)\n            model.eval()\n            models[mri_type] = model\n            print(f'Loaded segresnet_{mri_type}.pt model')\n\n        return models\n\n    def extract_features(self):\n\n        models = self.load_segmentation_models()\n        \n        print('\\nCreating training set features')\n        if self.train_features_path:\n            # Load precomputed training set features\n            df_train_features = pd.read_csv(f'{PROCESSED_DATA_PATH}/train_features.csv')\n            self.df_train = pd.concat([self.df_train, df_train_features], axis=1)\n        else:\n            # Create training set features\n            for case in tqdm(self.df_train['BraTS21ID'].values):\n                case_mpmris = os.listdir(f'{RAW_DATA_PATH}/train/{case}')\n                for mri_type in ['FLAIR', 'T1w']:\n                    \n                    mri = load_mri(\n                        mri_path=f'{RAW_DATA_PATH}/train/{case}/{mri_type}',\n                        window_width=None,\n                        window_center=None,\n                        voi_lut=None,\n                        new_spacing=None,\n                        reorder_plane=True,\n                        resize_shape=144,\n                        verbose=False\n                    )\n                    \n                    mri = (mri - mri.mean()) / mri.std()\n                    mri = torch.as_tensor(mri, dtype=torch.float)\n                    mri = torch.unsqueeze(mri, 0)\n                    mri = torch.unsqueeze(mri, 0)\n                    mri = mri.to(self.device)\n\n                    mask = torch.sigmoid(models[mri_type](mri)).detach().cpu().numpy()\n                    mask = np.squeeze(mask, axis=0)\n                    mask = np.round(mask)\n\n                    self.df_train.loc[self.df_train['BraTS21ID'] == case, f'{mri_type}_Whole_Tumor_Area'] = np.sum(mask[0])\n                    self.df_train.loc[self.df_train['BraTS21ID'] == case, f'{mri_type}_Tumor_Core_Area'] = np.sum(mask[1])\n                    self.df_train.loc[self.df_train['BraTS21ID'] == case, f'{mri_type}_Enhancing_Tumor_Area'] = np.sum(mask[2])\n        \n        print('Creating test set features')\n        # Create test set features\n        for case in tqdm(self.df_test['BraTS21ID'].values):\n            case_mpmris = os.listdir(f'{RAW_DATA_PATH}/test/{case}')\n            for mri_type in ['FLAIR', 'T1w']:\n                \n                mri = load_mri(\n                    mri_path=f'{RAW_DATA_PATH}/test/{case}/{mri_type}',\n                    window_width=None,\n                    window_center=None,\n                    voi_lut=None,\n                    new_spacing=None,\n                    reorder_plane=True,\n                    resize_shape=144,\n                    verbose=False\n                )\n                \n                # This block is required for a failproof pipeline\n                # Some excluded problematic MRIs in training set caused this block to trigger\n                if mri is None:\n                    continue\n\n                mri = (mri - mri.mean()) / mri.std()\n                mri = torch.as_tensor(mri, dtype=torch.float)\n                mri = torch.unsqueeze(mri, 0)\n                mri = torch.unsqueeze(mri, 0)\n                mri = mri.to(self.device)\n\n                mask = torch.sigmoid(models[mri_type](mri)).detach().cpu().numpy()\n                mask = np.squeeze(mask, axis=0)\n                mask = np.round(mask)\n\n                self.df_test.loc[self.df_test['BraTS21ID'] == case, f'{mri_type}_Whole_Tumor_Area'] = np.sum(mask[0])\n                self.df_test.loc[self.df_test['BraTS21ID'] == case, f'{mri_type}_Tumor_Core_Area'] = np.sum(mask[1])\n                self.df_test.loc[self.df_test['BraTS21ID'] == case, f'{mri_type}_Enhancing_Tumor_Area'] = np.sum(mask[2])\n                \n        return self.df_train.copy(deep=True), self.df_test.copy(deep=True)\n","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:37:43.321947Z","iopub.execute_input":"2021-10-15T07:37:43.322244Z","iopub.status.idle":"2021-10-15T07:37:43.349656Z","shell.execute_reply.started":"2021-10-15T07:37:43.322215Z","shell.execute_reply":"2021-10-15T07:37:43.347965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"segmentation_feature_extractor = SegmentationFeatureExtractor(\n    df_train=df_train,\n    df_test=df_test,\n    train_features_path=None\n)\n\n#df_train, df_test = segmentation_feature_extractor.extract_features()","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:38:25.741013Z","iopub.execute_input":"2021-10-15T07:38:25.741986Z","iopub.status.idle":"2021-10-15T07:38:25.797917Z","shell.execute_reply.started":"2021-10-15T07:38:25.741952Z","shell.execute_reply":"2021-10-15T07:38:25.795932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3D Classification","metadata":{}},{"cell_type":"code","source":"def load_predictions_and_evaluate_model(df, model_name, mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'), sigmoid=False):\n\n    print(f'\\n{\"-\" * 30}\\nEvaluating {model_name}\\n{\"-\" * 30}\\n')\n    for mri_type in mri_types:\n\n        predictions_column_name = f'{model_name[:-3]}_{mri_type}_predictions'\n        train_column_name = f'{model_name}_{mri_type}_predictions'\n        print(f'{mri_type}\\n{\"-\" * len(mri_type)}')\n\n        df_train_predictions = pd.read_csv(f'{PREDICTIONS_PATH}/train_{model_name}_{mri_type}_predictions.csv')\n        if sigmoid:\n            df_train_predictions[predictions_column_name] = 1 / (1 + np.exp(-df_train_predictions[predictions_column_name].values))\n            \n        df[train_column_name] = df_train_predictions[predictions_column_name].values\n\n        for fold in sorted(df['fold'].unique()):\n            _, val_idx = df.loc[df['fold'] != fold].index, df.loc[df_train['fold'] == fold].index\n            fold_score = roc_auc_score(df.loc[val_idx, 'MGMT_value'], df.loc[val_idx, train_column_name])\n            print(f'Fold {fold} - ROC AUC Score: {fold_score:.6f}')\n        oof_score = roc_auc_score(df['MGMT_value'], df[train_column_name])\n        print(f'{\"-\" * 30}\\nOOF ROC AUC Score: {oof_score:.6}\\n{\"-\" * 30}\\n')\n","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:51:16.068357Z","iopub.execute_input":"2021-10-15T07:51:16.068652Z","iopub.status.idle":"2021-10-15T07:51:16.083903Z","shell.execute_reply.started":"2021-10-15T07:51:16.068624Z","shell.execute_reply":"2021-10-15T07:51:16.082698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"load_predictions_and_evaluate_model(\n    df=df_train,\n    model_name='densenet121',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:51:16.797207Z","iopub.execute_input":"2021-10-15T07:51:16.797529Z","iopub.status.idle":"2021-10-15T07:51:16.924925Z","shell.execute_reply.started":"2021-10-15T07:51:16.797486Z","shell.execute_reply":"2021-10-15T07:51:16.924035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"load_predictions_and_evaluate_model(\n    df=df_train,\n    model_name='densenet169',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:51:17.389131Z","iopub.execute_input":"2021-10-15T07:51:17.389540Z","iopub.status.idle":"2021-10-15T07:51:17.510853Z","shell.execute_reply.started":"2021-10-15T07:51:17.389510Z","shell.execute_reply":"2021-10-15T07:51:17.509960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Blending","metadata":{}},{"cell_type":"code","source":"prediction_columns = [col for col in df_train.columns if col.endswith('predictions')]\n\nfig = plt.figure(figsize=(16, 16), dpi=100)\nsns.heatmap(\n    df_train[prediction_columns + ['MGMT_value']].corr(),\n    annot=True,\n    square=True,\n    cmap='coolwarm',\n    annot_kws={'size': 12},\n    fmt='.4f'\n)\n\nplt.tick_params(axis='x', labelsize=10, rotation=90)\nplt.tick_params(axis='y', labelsize=10, rotation=0)\nplt.title('Prediction Correlations', size=20, pad=20)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-15T07:53:35.770766Z","iopub.execute_input":"2021-10-15T07:53:35.771054Z","iopub.status.idle":"2021-10-15T07:53:36.845179Z","shell.execute_reply.started":"2021-10-15T07:53:35.771025Z","shell.execute_reply":"2021-10-15T07:53:36.844122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train['densenet_blend_predictions'] = (\n    (df_train['densenet121_FLAIR_predictions'] * 0.125) +\n    (df_train['densenet121_T1w_predictions'] * 0.125) +\n    (df_train['densenet121_T1wCE_predictions'] * 0.125) +\n    (df_train['densenet121_T2w_predictions'] * 0.125) +\n    (df_train['densenet169_FLAIR_predictions'] * 0.125) +\n    (df_train['densenet169_T1w_predictions'] * 0.125) +\n    (df_train['densenet169_T1wCE_predictions'] * 0.125) +\n    (df_train['densenet169_T2w_predictions'] * 0.125)\n)\n\n\nprint(f'Blend\\n{\"-\" * 5}')\nfor fold in sorted(df_train['fold'].unique()):\n    _, val_idx = df_train.loc[df_train['fold'] != fold].index, df_train.loc[df_train['fold'] == fold].index\n    fold_score = roc_auc_score(df_train.loc[val_idx, 'MGMT_value'], df_train.loc[val_idx, 'densenet_blend_predictions'])\n    print(f'Fold {fold} - ROC AUC Score: {fold_score:.6f}')\noof_score = roc_auc_score(df_train['MGMT_value'], df_train['densenet_blend_predictions'])\nprint(f'{\"-\" * 30}\\nOOF ROC AUC Score: {oof_score:.6}\\n{\"-\" * 30}\\n')","metadata":{"execution":{"iopub.status.busy":"2021-10-15T08:35:14.804800Z","iopub.execute_input":"2021-10-15T08:35:14.805071Z","iopub.status.idle":"2021-10-15T08:35:14.851046Z","shell.execute_reply.started":"2021-10-15T08:35:14.805043Z","shell.execute_reply":"2021-10-15T08:35:14.850207Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"code","source":"class MRIClassificationDataset(Dataset):\n\n    def __init__(self, cases, targets, mri_type, transforms=None):\n\n        self.cases = cases\n        self.targets = targets\n        self.mri_type = mri_type\n        self.transforms = transforms\n\n    def __len__(self):\n        return len(self.cases)\n\n    def __getitem__(self, idx):\n\n        \"\"\"\n        Get the idxth element in the dataset\n\n        Parameters\n        ----------\n        idx (int): Index of the sample (0 <= idx < len(self.cases))\n\n        Returns\n        -------\n        mri [torch.FloatTensor of shape (channel, depth, height, width)]: Preprocessed 4D mpMRI\n        target [torch.FloatTensor of shape (1)]: MGMT value\n        \"\"\"\n\n        mri = load_mri(\n            mri_path=f'{RAW_DATA_PATH}/test/{self.cases[idx]}/{self.mri_type}',\n            window_width=None,\n            window_center=None,\n            voi_lut=None,\n            new_spacing=None,\n            reorder_plane=True,\n            resize_shape=144,\n            verbose=False\n        )\n        if self.transforms is not None:\n            mri = self.transforms(mri)\n\n        mri = (mri - mri.mean()) / mri.std()\n        mri = torch.as_tensor(mri, dtype=torch.float)\n        mri = torch.unsqueeze(mri, 0)\n\n        if self.targets is not None:\n            target = self.targets[idx]\n            target = torch.as_tensor(target, dtype=torch.float)\n            target = torch.unsqueeze(target, 0)\n            return mri, target\n        else:\n            return mri\n","metadata":{"execution":{"iopub.status.busy":"2021-10-15T08:36:31.633702Z","iopub.execute_input":"2021-10-15T08:36:31.633979Z","iopub.status.idle":"2021-10-15T08:36:31.644212Z","shell.execute_reply.started":"2021-10-15T08:36:31.633953Z","shell.execute_reply":"2021-10-15T08:36:31.643236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class DenseNetModel(nn.Module):\n\n    def __init__(self, init_features, growth_rate, block_config, bn_size, dropout_prob):\n\n        super(DenseNetModel, self).__init__()\n\n        self.backbone = DenseNet(\n            spatial_dims=3,\n            in_channels=1,\n            out_channels=1,\n            init_features=init_features,\n            growth_rate=growth_rate,\n            block_config=block_config,\n            bn_size=bn_size,\n            dropout_prob=dropout_prob\n        )\n\n    def forward(self, x):\n        return self.backbone(x)\n    \n\nmodel_configs = {\n    'densenet121': {\n        'init_features': 64,\n        'growth_rate': 32,\n        'block_config': (6, 12, 24, 16),\n        'bn_size': 4,\n        'dropout_prob': 0\n    },\n    'densenet169': {\n        'init_features': 64,\n        'growth_rate': 32,\n        'block_config': (6, 12, 48, 32),\n        'bn_size': 4,\n        'dropout_prob': 0\n    }\n}","metadata":{"execution":{"iopub.status.busy":"2021-10-15T08:38:24.185747Z","iopub.execute_input":"2021-10-15T08:38:24.186066Z","iopub.status.idle":"2021-10-15T08:38:24.194738Z","shell.execute_reply.started":"2021-10-15T08:38:24.186039Z","shell.execute_reply":"2021-10-15T08:38:24.193788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"folds_to_use = {\n    'densenet121': {\n        'FLAIR': [1, 2, 3, 4, 5],\n        'T1w': [1, 2, 3, 4, 5],\n        'T1wCE': [1, 2, 3, 4, 5],\n        'T2w': [1, 2, 3, 4, 5]\n    },\n    'densenet169': {\n        'FLAIR': [1, 2, 3, 4, 5],\n        'T1w': [1, 2, 3, 4, 5],\n        'T1wCE': [1, 2, 3, 4, 5],\n        'T2w': [1, 2, 3, 4, 5]\n    }\n}\n\n\ndef inference(df, model_name, mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'), sigmoid=True):\n\n    print(f'\\n{\"-\" * 30}\\nRunning {model_name} for Inference ({mri_types}\\n{\"-\" * 30}')\n    \n    for mri_type in mri_types:\n        \n        predictions_column_name = f'{model_name}_{mri_type}_predictions'\n        df_test[predictions_column_name] = 0\n\n        for fold in sorted(df_train['fold'].unique()):\n            \n            if fold not in folds_to_use[model_name][mri_type]:\n                continue\n\n            test_dataset = MRIClassificationDataset(\n                cases=df_test['BraTS21ID'].values.tolist(),\n                targets=None,\n                mri_type=mri_type,\n                transforms=None\n            )\n            test_loader = DataLoader(\n                test_dataset,\n                batch_size=16,\n                sampler=SequentialSampler(test_dataset),\n                pin_memory=True,\n                drop_last=False,\n                num_workers=4,\n            )\n\n            device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')\n            model = DenseNetModel(**model_configs[model_name])\n            model.load_state_dict(torch.load(f'{MODELS_PATH}/{model_name}/{model_name}_{mri_type}_fold{fold}.pt'))\n            model.to(device)\n            model.eval()\n\n            predictions = []\n            with torch.no_grad():\n                for mri in tqdm(test_loader):\n                    mri = mri.to(device)\n                    output = model(mri)\n                    output = output.detach().cpu().numpy().flatten().tolist()\n                    predictions += output\n\n            df_test[predictions_column_name] += predictions\n            print(f'Finished Inference for {model_name} Model Fold {fold} ({mri_type})')\n            \n        df_test[predictions_column_name] = df_test[predictions_column_name] / len(folds_to_use[model_name][mri_type])\n        if sigmoid:\n            df_test[predictions_column_name] = 1 / (1 + np.exp(-df_test[predictions_column_name].values))\n","metadata":{"execution":{"iopub.status.busy":"2021-10-15T08:48:28.909664Z","iopub.execute_input":"2021-10-15T08:48:28.909954Z","iopub.status.idle":"2021-10-15T08:48:28.924806Z","shell.execute_reply.started":"2021-10-15T08:48:28.909925Z","shell.execute_reply":"2021-10-15T08:48:28.923666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inference(\n    df=df_test,\n    model_name='densenet121',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-15T08:48:29.401636Z","iopub.execute_input":"2021-10-15T08:48:29.402580Z","iopub.status.idle":"2021-10-15T09:12:52.362185Z","shell.execute_reply.started":"2021-10-15T08:48:29.402548Z","shell.execute_reply":"2021-10-15T09:12:52.360023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inference(\n    df=df_test,\n    model_name='densenet169',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"execution":{"iopub.status.busy":"2021-10-15T09:17:09.683779Z","iopub.execute_input":"2021-10-15T09:17:09.684070Z","iopub.status.idle":"2021-10-15T09:42:47.510798Z","shell.execute_reply.started":"2021-10-15T09:17:09.684028Z","shell.execute_reply":"2021-10-15T09:42:47.508485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission","metadata":{}},{"cell_type":"code","source":"df_test['MGMT_value'] = (\n    (df_test['densenet121_FLAIR_predictions'] * 0.25) +\n    #(df_test['densenet121_T1w_predictions'] * 0.125) +\n    #(df_test['densenet121_T1wCE_predictions'] * 0.125) +\n    (df_test['densenet121_T2w_predictions'] * 0.25) +\n    (df_test['densenet169_FLAIR_predictions'] * 0.25) +\n    #(df_test['densenet169_T1w_predictions'] * 0.125) +\n    #(df_test['densenet169_T1wCE_predictions'] * 0.125) +\n    (df_test['densenet169_T2w_predictions'] * 0.25)\n)\ndf_test[['BraTS21ID', 'MGMT_value']].to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-10-15T09:44:09.135493Z","iopub.execute_input":"2021-10-15T09:44:09.136210Z","iopub.status.idle":"2021-10-15T09:44:09.152675Z","shell.execute_reply.started":"2021-10-15T09:44:09.136180Z","shell.execute_reply":"2021-10-15T09:44:09.151712Z"},"trusted":true},"execution_count":null,"outputs":[]}]}