{"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":{"papermill":{"duration":28.282223,"end_time":"2021-10-13T15:16:52.388015","exception":false,"start_time":"2021-10-13T15:16:24.105792","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:40:08.760526Z","iopub.execute_input":"2021-10-13T20:40:08.761613Z","iopub.status.idle":"2021-10-13T20:40:18.454915Z","shell.execute_reply.started":"2021-10-13T20:40:08.76157Z","shell.execute_reply":"2021-10-13T20:40:18.453702Z"},"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\nfrom sklearn.metrics import roc_auc_score\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom torch.utils.data import Dataset, DataLoader, RandomSampler, SequentialSampler\nfrom monai.networks.nets import SegResNet, DenseNet","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":7.534199,"end_time":"2021-10-13T15:16:59.939612","exception":false,"start_time":"2021-10-13T15:16:52.405413","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T21:02:11.71627Z","iopub.execute_input":"2021-10-13T21:02:11.716879Z","iopub.status.idle":"2021-10-13T21:02:11.724359Z","shell.execute_reply.started":"2021-10-13T21:02:11.716841Z","shell.execute_reply":"2021-10-13T21:02:11.723404Z"},"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/brain-segresnet'\nPREDICTIONS_PATH = '../input/rsnamiccai-btrc-dataset/predictions'","metadata":{"papermill":{"duration":0.021073,"end_time":"2021-10-13T15:16:59.975608","exception":false,"start_time":"2021-10-13T15:16:59.954535","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:07.469936Z","iopub.execute_input":"2021-10-13T20:59:07.470277Z","iopub.status.idle":"2021-10-13T20:59:07.47621Z","shell.execute_reply.started":"2021-10-13T20:59:07.470247Z","shell.execute_reply":"2021-10-13T20:59:07.474498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train = pd.read_csv(f'{PROCESSED_DATA_PATH}/train_task2.csv', dtype={'BraTS21ID': 'object', 'MGMT_value': np.uint8, 'fold': np.uint8})\ndf_test = pd.read_csv(f'{RAW_DATA_PATH}/sample_submission.csv', usecols=['BraTS21ID'], dtype={'BraTS21ID': 'object'})\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":{"papermill":{"duration":0.048908,"end_time":"2021-10-13T15:17:00.038984","exception":false,"start_time":"2021-10-13T15:16:59.990076","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:07.832818Z","iopub.execute_input":"2021-10-13T20:59:07.833104Z","iopub.status.idle":"2021-10-13T20:59:07.85477Z","shell.execute_reply.started":"2021-10-13T20:59:07.833077Z","shell.execute_reply":"2021-10-13T20:59:07.853919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## DICOM Preprocessing","metadata":{"papermill":{"duration":0.014559,"end_time":"2021-10-13T15:17:00.069486","exception":false,"start_time":"2021-10-13T15:17:00.054927","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":0.055035,"end_time":"2021-10-13T15:17:00.139275","exception":false,"start_time":"2021-10-13T15:17:00.08424","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:09.099372Z","iopub.execute_input":"2021-10-13T20:59:09.099793Z","iopub.status.idle":"2021-10-13T20:59:09.143744Z","shell.execute_reply.started":"2021-10-13T20:59:09.099753Z","shell.execute_reply":"2021-10-13T20:59:09.142433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Tumor Segmentation","metadata":{"papermill":{"duration":0.01696,"end_time":"2021-10-13T15:17:00.174251","exception":false,"start_time":"2021-10-13T15:17:00.157291","status":"completed"},"tags":[]}},{"cell_type":"raw","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":{"papermill":{"duration":0.044897,"end_time":"2021-10-13T15:17:00.235805","exception":false,"start_time":"2021-10-13T15:17:00.190908","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:41:27.575141Z","iopub.execute_input":"2021-10-13T20:41:27.575514Z","iopub.status.idle":"2021-10-13T20:41:27.602742Z","shell.execute_reply.started":"2021-10-13T20:41:27.575473Z","shell.execute_reply":"2021-10-13T20:41:27.60134Z"}}},{"cell_type":"raw","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.execute_input":"2021-10-13T15:17:00.31463Z","iopub.status.busy":"2021-10-13T15:17:00.313919Z","iopub.status.idle":"2021-10-13T15:17:00.316156Z","shell.execute_reply":"2021-10-13T15:17:00.316562Z","shell.execute_reply.started":"2021-10-13T14:49:35.904231Z"},"papermill":{"duration":0.064702,"end_time":"2021-10-13T15:17:00.316697","exception":false,"start_time":"2021-10-13T15:17:00.251995","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def load_predictions_and_evaluate_model(df, model_name, mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'), sigmoid=False):\n\n    \"\"\"\n    Evaluate predictions of the given model for every mpMRI type\n\n    Parameters\n    ----------\n    df_train [pandas.DataFrame of shape (n_samples)]: DataFrame with MGMT_value column\n    model_name (str): Model name used in the predictions files\n    \"\"\"\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}_{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[predictions_column_name] = df_train_predictions[predictions_column_name].values\n\n        for 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_predictions.loc[val_idx, predictions_column_name])\n            print(f'Fold {fold} - ROC AUC Score: {fold_score:.6f}')\n        oof_score = roc_auc_score(df_train['MGMT_value'], df_train_predictions[predictions_column_name])\n        print(f'{\"-\" * 30}\\nOOF ROC AUC Score: {oof_score:.6}\\n{\"-\" * 30}\\n')\n","metadata":{"papermill":{"duration":0.025573,"end_time":"2021-10-13T15:17:00.357208","exception":false,"start_time":"2021-10-13T15:17:00.331635","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:15.589555Z","iopub.execute_input":"2021-10-13T20:59:15.590189Z","iopub.status.idle":"2021-10-13T20:59:15.600371Z","shell.execute_reply.started":"2021-10-13T20:59:15.590156Z","shell.execute_reply":"2021-10-13T20:59:15.599165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"load_predictions_and_evaluate_model(\n    df=df_train,\n    model_name='densenet',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"execution":{"iopub.execute_input":"2021-10-13T15:17:00.390495Z","iopub.status.busy":"2021-10-13T15:17:00.389938Z","iopub.status.idle":"2021-10-13T15:17:00.54427Z","shell.execute_reply":"2021-10-13T15:17:00.544759Z","shell.execute_reply.started":"2021-10-13T14:49:43.68283Z"},"papermill":{"duration":0.173015,"end_time":"2021-10-13T15:17:00.544931","exception":false,"start_time":"2021-10-13T15:17:00.371916","status":"completed"},"tags":[]}},{"cell_type":"raw","source":"df_train['densenet_blend_predictions'] = (df_train['densenet_FLAIR_predictions'] * 0.25) +\\\n                                         (df_train['densenet_T1w_predictions'] * 0.25) +\\\n                                         (df_train['densenet_T1wCE_predictions'] * 0.25) +\\\n                                         (df_train['densenet_T2w_predictions'] * 0.25)\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.execute_input":"2021-10-13T15:17:00.593668Z","iopub.status.busy":"2021-10-13T15:17:00.59049Z","iopub.status.idle":"2021-10-13T15:17:00.615519Z","shell.execute_reply":"2021-10-13T15:17:00.616166Z","shell.execute_reply.started":"2021-10-13T14:51:57.216489Z"},"papermill":{"duration":0.052905,"end_time":"2021-10-13T15:17:00.61634","exception":false,"start_time":"2021-10-13T15:17:00.563435","status":"completed"},"tags":[]}},{"cell_type":"code","source":"class MRIClassificationDataset(Dataset):\n\n    def __init__(self, cases, targets, mri_type, tta=0, transforms=None):\n\n        self.cases = cases\n        self.targets = targets\n        self.mri_type = mri_type\n        self.transforms = transforms\n        self.tta = tta\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        if self.tta==1:\n            mri = mri[::-1]\n        else:\n            pass\n            # Other tta here..\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":{"papermill":{"duration":0.029554,"end_time":"2021-10-13T15:17:00.662754","exception":false,"start_time":"2021-10-13T15:17:00.6332","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:19.298134Z","iopub.execute_input":"2021-10-13T20:59:19.298682Z","iopub.status.idle":"2021-10-13T20:59:19.316878Z","shell.execute_reply.started":"2021-10-13T20:59:19.298622Z","shell.execute_reply":"2021-10-13T20:59:19.315768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CustomModel(nn.Module):\n    def __init__(self):\n        super().__init__()\n\n        self.tumor = 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        \n        self.srn = 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        del self.srn.up_layers, self.srn.conv_final, self.srn.up_samples\n\n        dims = 120-16-8\n        self.tail = nn.Sequential(\n            nn.BatchNorm1d(dims),\n            nn.Linear(dims, dims, bias=False),\n            nn.GELU(),\n            \n            nn.Dropout(0.1),\n            nn.Linear(dims,1)\n        )\n        \n        # Freeze Stuff\n        self.srn.convInit.conv.weight.requires_grad = False\n        self.srn.down_layers[0][1].norm1.weight.requires_grad = False\n        self.srn.down_layers[0][1].norm1.bias.requires_grad = False\n        self.srn.down_layers[0][1].norm2.weight.requires_grad = False\n        self.srn.down_layers[0][1].norm2.bias.requires_grad = False\n        self.srn.down_layers[0][1].conv1.conv.weight.requires_grad = False\n        self.srn.down_layers[0][1].conv2.conv.weight.requires_grad = False\n        for p in self.tumor.parameters():\n            p.requires_grad = False\n    \n    def forward_path(self, x):\n        xs = self.tumor.convInit(x)\n        if self.tumor.dropout_prob is not None:\n            xs = self.tumor.dropout(xs)\n            \n        down_x = []\n        for down in self.tumor.down_layers:\n            xs = down(xs)\n            down_x.append(xs)\n        down_x.reverse()\n        for i, (up, upl) in enumerate(zip(self.tumor.up_samples, self.tumor.up_layers)):\n            xs = up(xs) + down_x[i + 1]\n            xs = upl(xs)\n        if self.tumor.use_conv_final:\n            xs = self.tumor.conv_final(xs)\n            \n        # probably should be max...\n        x *= xs.mean(axis=1, keepdim=True)#.sigmoid()\n        \n        x = self.srn.convInit(x)\n        if self.srn.dropout_prob is not None:\n            x = self.srn.dropout(x)\n\n        down_x = []\n        for down in self.srn.down_layers:\n            x = down(x)\n            down_x.append(x)\n\n        return down_x\n\n    def forward(self, x):\n        stack = self.forward_path(x)\n\n        x = torch.cat([\n            F.adaptive_max_pool3d(y, 1).flatten(1)\n            for y in stack[2:]\n        ], axis=1) + torch.cat([\n            F.adaptive_avg_pool3d(y, 1).flatten(1)\n            for y in stack[2:]\n        ], axis=1)\n        \n        # stack\n#         torch.Size([8, 8, 144, 144, 144])\n#         torch.Size([8, 16, 72, 72, 72])\n#         torch.Size([8, 32, 36, 36, 36])\n#         torch.Size([8, 64, 18, 18, 18])\n\n        return self.tail(x)","metadata":{"execution":{"iopub.status.busy":"2021-10-13T20:59:52.852115Z","iopub.execute_input":"2021-10-13T20:59:52.852472Z","iopub.status.idle":"2021-10-13T20:59:52.876941Z","shell.execute_reply.started":"2021-10-13T20:59:52.852405Z","shell.execute_reply":"2021-10-13T20:59:52.875839Z"},"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\n        return self.backbone(x)\n    \n\nmodel_configs = {\n    'densenet': {\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    'CustomModel': {},\n}\n\n\ndef get_model(model_name, **kwargs):\n    \n    if model_name == 'densenet':\n        return DenseNetModel(**kwargs)\n    \n    if model_name == 'CustomModel':\n        return CustomModel()\n","metadata":{"papermill":{"duration":0.025561,"end_time":"2021-10-13T15:17:00.704092","exception":false,"start_time":"2021-10-13T15:17:00.678531","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:53.067815Z","iopub.execute_input":"2021-10-13T20:59:53.069269Z","iopub.status.idle":"2021-10-13T20:59:53.083882Z","shell.execute_reply.started":"2021-10-13T20:59:53.069224Z","shell.execute_reply":"2021-10-13T20:59:53.082063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"folds_to_use = {\n    'densenet': {\n        'FLAIR': [1, 3, 4, 5],\n        'T1w': [2],\n        'T1wCE': [1, 3, 4],\n        'T2w': [1, 3, 4, 5]\n    },\n    'CustomModel':  {\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}","metadata":{"execution":{"iopub.status.busy":"2021-10-13T20:59:53.188419Z","iopub.execute_input":"2021-10-13T20:59:53.189161Z","iopub.status.idle":"2021-10-13T20:59:53.201085Z","shell.execute_reply.started":"2021-10-13T20:59:53.189116Z","shell.execute_reply":"2021-10-13T20:59:53.199906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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 = get_model(model_name, **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 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":{"papermill":{"duration":0.030508,"end_time":"2021-10-13T15:17:00.750675","exception":false,"start_time":"2021-10-13T15:17:00.720167","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:53.349061Z","iopub.execute_input":"2021-10-13T20:59:53.349828Z","iopub.status.idle":"2021-10-13T20:59:53.381778Z","shell.execute_reply.started":"2021-10-13T20:59:53.349786Z","shell.execute_reply":"2021-10-13T20:59:53.380555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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    device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')\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        test_dataset = MRIClassificationDataset(\n            cases=df_test['BraTS21ID'].values.tolist(),\n            targets=None,\n            mri_type=mri_type,\n            transforms=None,\n            tta=0 # TTA <--- we can wrap in for loop and try 1,2,3,4 different augs, time permitting\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        for run in range(2):\n            for fold in range(1,11):\n                # Actually 5 fold model, but snapshot ensemble of best score + best loss models.\n                model = get_model(model_name, **model_configs[model_name])\n                model.load_state_dict(torch.load(f'{MODELS_PATH}/{mri_type}_run{run}_{fold}.pth', map_location='cpu')['model'])\n                model.to(device)\n                model.eval()\n\n                predictions = []\n                with torch.no_grad():\n                    for mri in 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}. Fold {fold} Run {run} - ({mri_type})')\n\n        df_test[predictions_column_name] = df_test[predictions_column_name] / (5*2*2)\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-13T21:00:50.774257Z","iopub.execute_input":"2021-10-13T21:00:50.774824Z","iopub.status.idle":"2021-10-13T21:00:50.788721Z","shell.execute_reply.started":"2021-10-13T21:00:50.774777Z","shell.execute_reply":"2021-10-13T21:00:50.787301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inference(\n    df=df_test,\n    model_name='CustomModel',\n    mri_types=('FLAIR', 'T1w', 'T1wCE', 'T2w'),\n    sigmoid=True\n)","metadata":{"papermill":{"duration":972.122915,"end_time":"2021-10-13T15:33:12.890164","exception":false,"start_time":"2021-10-13T15:17:00.767249","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T21:02:18.689312Z","iopub.execute_input":"2021-10-13T21:02:18.689677Z","iopub.status.idle":"2021-10-13T21:24:19.659975Z","shell.execute_reply.started":"2021-10-13T21:02:18.689644Z","shell.execute_reply":"2021-10-13T21:24:19.657787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test","metadata":{"papermill":{"duration":0.064845,"end_time":"2021-10-13T15:33:12.975003","exception":false,"start_time":"2021-10-13T15:33:12.910158","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:02.703546Z","iopub.execute_input":"2021-10-13T20:59:02.703825Z","iopub.status.idle":"2021-10-13T20:59:02.727405Z","shell.execute_reply.started":"2021-10-13T20:59:02.703787Z","shell.execute_reply":"2021-10-13T20:59:02.72662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test['MGMT_value'] = (\n    df_test['CustomModel_FLAIR_predictions'].values +\n    df_test['CustomModel_T1w_predictions'].values +\n    df_test['CustomModel_T1wCE_predictions'].values +\n    df_test['CustomModel_T2w_predictions'].values\n)\ndf_test[['BraTS21ID', 'MGMT_value']].to_csv('submission.csv', index=False)","metadata":{"papermill":{"duration":0.037513,"end_time":"2021-10-13T15:33:13.032773","exception":false,"start_time":"2021-10-13T15:33:12.99526","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-10-13T20:59:02.860657Z","iopub.execute_input":"2021-10-13T20:59:02.861243Z","iopub.status.idle":"2021-10-13T20:59:02.939433Z","shell.execute_reply.started":"2021-10-13T20:59:02.861213Z","shell.execute_reply":"2021-10-13T20:59:02.937721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"papermill":{"duration":0.019566,"end_time":"2021-10-13T15:33:13.071998","exception":false,"start_time":"2021-10-13T15:33:13.052432","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]}]}