{"cells":[{"metadata":{},"cell_type":"markdown","source":"## 1. Dependencies\n\nPackages are installed from a dataset because competition rules restrict internet usage. Packages are in [this](https://www.kaggle.com/gunesevitan/osic-pulmonary-fibrosis-datasets) dataset.\n\n* `gdcm==2.8.9`\n* `nrrd==0.4.2`\n\nTwo scans in training set requires `gdcm` to read .dcm files, but private test set scans can be read without it. `nrrd` is used for reading .nrrd files (masks)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# gdcm-2.8.9\n!cp ../input/osic-pulmonary-fibrosis-datasets/gdcm.tar .\n!tar -xvzf gdcm.tar\n!conda install --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2\n\n!pip install ../input/osic-pulmonary-fibrosis-datasets/pynrrd-master -f ./ --no-index # nrrd-0.4.2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* `efficientnet_pytorch==0.6.3`\n* `pretrainedmodels==0.7.4`\n* `timm==0.1.20`\n\n`timm-0.1.20` has to be installed initially because `segmentation_models_pytorch` requires that specific version. It is upgraded after installing `segmentation_models_pytorch`."},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install ../input/osic-pulmonary-fibrosis-datasets/efficientnet_pytorch-0.6.3 -f ./ --no-index # efficientnet_pytorch-0.6.3\n!pip install ../input/osic-pulmonary-fibrosis-datasets/pretrained-models.pytorch -f ./ --no-index # pretrainedmodels-0.7.4\n!pip install ../input/osic-pulmonary-fibrosis-datasets/timm-0.1.20-py3-none-any.whl # timm-0.1.20","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* `segmentation_models_pytorch==0.1.2`\n* `timm==0.2.1`\n\nA single line is changed in `segmentation_models_pytorch.encoders.timm_efficientnet.py` module because it was throwing an import error. Line 4 is changed from:\n\n`from timm.models.efficientnet import EfficientNet, Swish` \n\nto\n\n`\nfrom timm.models.layers.activations import Swish\nfrom timm.models.efficientnet import EfficientNet\n`\n\nAfter `segmentation_models_pytorch` is installed successfully, `timm` version is upgraded from `0.1.20` to `0.2.1`. Upgrading `timm` displays this error; `ERROR: segmentation-models-pytorch 0.1.2 has requirement timm==0.1.20, but you'll have timm 0.2.1 which is incompatible.`, but it is installed without any problem."},{"metadata":{"_kg_hide-input":false,"_kg_hide-output":false,"trusted":true},"cell_type":"code","source":"!pip install ../input/osic-pulmonary-fibrosis-datasets/segmentation_models_pytorch-0.1.2 -f ./ --no-index # segmentation_models_pytorch-0.1.2\n!pip install ../input/osic-pulmonary-fibrosis-datasets/timm-0.2.1-py3-none-any.whl # timm-0.2.1","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\n\nimport os\nimport pickle\nimport random\nimport gc\nfrom glob import glob\nfrom tqdm import tqdm\n\nimport numpy as np\nimport pandas as pd\npd.set_option('display.max_rows', 500)\npd.set_option('display.max_columns', 500)\npd.set_option('display.width', 1000)\nfrom scipy.stats import skew, mode, kurtosis\nfrom scipy.signal import find_peaks\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport cv2\nimport gdcm\nimport pydicom\nimport nrrd\n\nfrom sklearn.preprocessing import MinMaxScaler, StandardScaler\nfrom sklearn.linear_model import LinearRegression\n\nimport torch\nfrom torch import nn\nfrom torch.nn import functional as F\nfrom torch.utils.data import Dataset, DataLoader\nimport segmentation_models_pytorch as smp\n\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.layers import Input, Dense, Lambda, Dropout, BatchNormalization, GaussianDropout, Activation\nfrom tensorflow.keras.regularizers import l1, l2\nfrom tensorflow.keras.optimizers import Adam, Nadam, RMSprop\nfrom tensorflow.keras.callbacks import Callback\nfrom tensorflow.keras.utils import Sequence\n\nSEED = 1337\n\ndef seed_everything(seed):\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    random.seed(seed)    \n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    torch.manual_seed(seed)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ndf_test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\ndf_submission = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\n\nprint(f'Training Set Shape = {df_train.shape} - Patients = {df_train[\"Patient\"].nunique()}')\nprint(f'Training Set Memory Usage = {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set Shape = {df_test.shape} - Patients = {df_test[\"Patient\"].nunique()}')\nprint(f'Test Set Memory Usage = {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Sample Submission Shape = {df_submission.shape}')\nprint(f'Sample Submission Memory Usage = {df_submission.memory_usage().sum() / 1024 ** 2:.2f} MB')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Preprocessing"},{"metadata":{},"cell_type":"markdown","source":"### 2.1. Tabular Data\n\nDuplicate measurements on the same `[Patient, Weeks]` groups are averaged and duplicate rows are dropped.\n\nCategorical features; `Sex` and `SmokingStatus` are label encoded and one hot encoded.\n\n`Weeks_Passed` feature is created for capturing the effects of time on `FVC`, and `Age` feature is incremented by `1 / 52` for every passed week.\n\nFirst measurement of every patient is provided in test set and that feature is named `FVC_Baseline`. First `Percent` value of every patient is used as a feature, not the continous values. Scale can be set to `True` for scaling continuous features but scaled features yield worse results."},{"metadata":{},"cell_type":"markdown","source":"### 2.2. Cross-validation\n\nImplemented three different cv schemes because each of them aren't always correlated with the public leaderboard score. Voting is done between those cv schemes based on the score change. Those custom cv schemes are:\n\n* Double Stratified Shuffled KFold on `Patient` groups. Patients are stratified by `Sex` and `SmokingStatus` features and groups are shuffled.\n* Cluster Stratified Shuffled KFold on `Patient` groups. Patients are clustered by their last two `FVC` values and those clusters are used for stratification.\n* Regular Shuffled KFold on `Patient` groups. Patients are shuffled without any stratification.\n\nAll of the cv schemes above are splitting data on `Patient` level so there is no leak between folds. Data is split into 2 folds because of the low patient count (176)."},{"metadata":{"trusted":true},"cell_type":"code","source":"class TabularDataPreprocessor:\n    \n    def __init__(self, train, test, submission, n_folds, shuffle, ohe, scale):\n        \n        self.train = train.copy(deep=True)\n        self.train.sort_values(by=['Patient', 'Weeks'], inplace=True)        \n        self.test = test.copy(deep=True)\n        self.submission = submission.copy(deep=True)\n                \n        self.n_folds = n_folds\n        self.shuffle = shuffle\n        \n        self.ohe = ohe\n        self.scale = scale\n        \n    def drop_duplicates(self):\n        \n        \"\"\"        \n        Calculate mean FVC and mean Percent of [Patient, Weeks] groups and drop duplicate rows\n        This operation takes the mean of multiple measurements in a single week and uses it\n        \"\"\"\n        \n        self.train['FVC'] = self.train.groupby(['Patient', 'Weeks'])['FVC'].transform('mean')\n        self.train['Percent'] = self.train.groupby(['Patient', 'Weeks'])['Percent'].transform('mean')\n        self.train.drop_duplicates(inplace=True)\n        self.train.reset_index(drop=True, inplace=True)\n        \n    def label_encode(self):\n        \n        \"\"\"\n        Label Encode categorical features\n        \"\"\"\n                    \n        for df in [self.train, self.test]:\n            df['Sex'] = df['Sex'].map({'Male': 0, 'Female': 1})\n            df['Sex'] = df['Sex'].astype(np.uint8)\n            df['SmokingStatus'] = df['SmokingStatus'].map({'Never smoked': 0, 'Ex-smoker': 1, 'Currently smokes': 2})\n            df['SmokingStatus'] = df['SmokingStatus'].astype(np.uint8)\n            \n    def one_hot_encode(self):\n        \n        \"\"\"\n        One-hot Encode categorical features\n        \"\"\"\n        \n        for df in [self.train, self.test]:\n            df['Male'] = 0\n            df['Female'] = 0\n            df.loc[df['Sex'] == 0, 'Male'] = 1\n            df.loc[df['Sex'] == 1, 'Female'] = 1\n            \n            df['Never smoked'] = 0\n            df['Ex-smoker'] = 0\n            df['Currently smokes'] = 0\n            df.loc[df['SmokingStatus'] == 0, 'Never smoked'] = 1\n            df.loc[df['SmokingStatus'] == 1, 'Ex-smoker'] = 1\n            df.loc[df['SmokingStatus'] == 2, 'Currently smokes'] = 1\n            \n            for encoded_col in ['Male', 'Female', 'Never smoked', 'Ex-smoker', 'Currently smokes']:\n                df[encoded_col] = df[encoded_col].astype(np.uint8)\n\n    def create_folds(self):\n        \n        \"\"\"\n        Creates n number of folds for three different cross-validation schemes (n should be selected as 2 because of the low patient count)\n            \n        1. Double Stratified Shuffled Folds\n        -----------------------------------\n        Patients are stratified by Sex and SmokingStatus features and groups are shuffled\n        Patients listed below are split into n folds with all of their FVC measurements\n        \n        Male_Ex-smoker             106 / n\n        Male_Never smoked           26 / n\n        Female_Never smoked         23 / n\n        Female_Ex-smoker            12 / n\n        Male_Currently smokes        7 / n\n        Female_Currently smokes      2 / n\n        \n        2. Cluster Stratified Shuffled Folds\n        ------------------------------------\n        Patients are clustered by their last two FVC values and clusters are stratified\n        Patients listed below are split into n folds with all of their FVC measurements\n        \n        Cluster 2    167 / n\n        Cluster 3      7 / n\n        Cluster 1      2 / n  \n        \n        3. Regular Shuffled Folds\n        -------------------------\n        Patient groups are shuffled into n folds        \n        Patient Count 176 / n\n        \"\"\"\n        \n        # Double Stratified Shuffled Folds\n        self.train['Sex_SmokingStatus'] = self.train['Sex'].astype(str) + '_' + self.train['SmokingStatus'].astype(str)\n        for group in self.train['Sex_SmokingStatus'].unique():\n            patients = self.train[self.train['Sex_SmokingStatus'] == group]['Patient'].unique()\n            \n            if self.shuffle:\n                np.random.seed(SEED)\n                np.random.shuffle(patients)\n                \n            for fold, patient_group in enumerate(np.array_split(patients, self.n_folds), 1):\n                self.train.loc[self.train['Patient'].isin(patient_group), 'CV1_Fold'] = fold\n               \n        # Cluster Stratified Shuffled Folds\n        for patient_name in self.train['Patient'].unique():\n            z = (self.train[(self.train['Patient'] == patient_name)]['FVC'].values[-2:] - self.train[(self.train['Patient'] == patient_name)]['FVC'].values[-2:].mean()) / self.train[(self.train['Patient'] == patient_name)]['FVC'].values[-2:].std()\n            reg = LinearRegression(normalize=True).fit(self.train[(self.train['Patient'] == patient_name)]['Weeks'].values[-2:].reshape(-1, 1), z)\n\n            self.train.loc[self.train['Patient'] == patient_name, 'Intercept'] = reg.intercept_\n            self.train.loc[self.train['Patient'] == patient_name, 'Coef'] = reg.coef_[0]\n            \n        self.train.loc[self.train['Coef'] > 0.4, 'Cluster'] = 1\n        self.train.loc[(self.train['Coef'] < 0.4) & (self.train['Coef'] > -0.4), 'Cluster'] = 2\n        self.train.loc[self.train['Coef'] < -0.4, 'Cluster'] = 3\n\n        for group in self.train['Cluster'].unique():\n            patients = self.train[self.train['Cluster'] == group]['Patient'].unique()\n\n            if self.shuffle:\n                np.random.seed(SEED)\n                np.random.shuffle(patients)\n\n            for fold, patient_group in enumerate(np.array_split(patients, self.n_folds), 1):\n                self.train.loc[self.train['Patient'].isin(patient_group), 'CV2_Fold'] = fold\n               \n        # Regular Shuffled Folds\n        patients = self.train['Patient'].unique()\n        np.random.seed(SEED)\n        np.random.shuffle(patients)\n        \n        for fold, patient_group in enumerate(np.array_split(patients, self.n_folds), 1):\n            self.train.loc[self.train['Patient'].isin(patient_group), 'CV3_Fold'] = fold\n            \n        self.train.drop(columns=['Sex_SmokingStatus', 'Intercept', 'Coef', 'Cluster'], inplace=True)        \n    \n    def create_tabular_features(self):\n        \n        self.drop_duplicates()\n        self.create_folds()\n        self.label_encode()\n        if self.ohe:\n            self.one_hot_encode()\n            \n        # Training set preparation\n        self.train['Type'] = 'Train'\n        self.train['Weeks_Passed'] = self.train['Weeks'] - self.train.groupby('Patient')['Weeks'].transform('min')\n        self.train['FVC_Baseline'] = self.train.groupby('Patient')['FVC'].transform('first')\n        \n        # Use Baseline Percent or Percent for training set?\n        # Percent obviously increases CV score since it correlates well with FVC\n        # However, only first Percent value is given for test set\n        # In that case, using Percent by itself is wrong and doesn't make sense\n        # Even though it is wrong, it increases public leaderboard score as well\n        # Can we trust public leaderboard improvement and use Percent?\n        self.train['Percent'] = self.train.groupby('Patient')['Percent'].transform('first')\n        # Comment this line out for using Percent feature\n            \n        # Test set preparation\n        self.submission['Type'] = 'Test'\n        self.submission['Patient'] = self.submission['Patient_Week'].apply(lambda x: x.split('_')[0]).astype(str)\n        self.submission['Weeks'] = self.submission['Patient_Week'].apply(lambda x: x.split('_')[1]).astype(int)\n        self.submission.drop(columns=['Patient_Week', 'FVC', 'Confidence'], inplace=True)\n        self.test = self.submission.merge(self.test.rename(columns={'Weeks': 'Weeks_Baseline', 'FVC': 'FVC_Baseline'}), how='left', on='Patient')\n        self.test['Weeks_Passed'] = self.test['Weeks'] - self.test['Weeks_Baseline']\n        self.test.drop(columns=['Weeks_Baseline'], inplace=True)\n        \n        df_all = pd.concat([self.train, self.test], ignore_index=True, axis=0)        \n        df_all['Age'] += (df_all['Weeks_Passed'] / 52)\n        df_all['Age'] = df_all['Age'].astype(np.float32)\n        df_all['FVC_Baseline'] = df_all['FVC_Baseline'].astype(np.float32)\n        df_all['Percent'] = df_all['Percent'].astype(np.float32)\n        df_all['Weeks_Passed'] = df_all['Weeks_Passed'].astype(np.float32)        \n        df_all['Weeks'] = df_all['Weeks'].astype(np.int16)\n\n        df_all['FVC'] = df_all['FVC'].astype(np.float32)\n        \n        if self.scale:\n            scale_features = ['FVC_Baseline', 'Age', 'Percent', 'Weeks_Passed']\n            scaler = MinMaxScaler()\n            df_all.loc[:, scale_features] = scaler.fit_transform(df_all.loc[:, scale_features])\n        \n        df_train = df_all.loc[df_all['Type'] == 'Train', :].drop(columns=['Type'])\n        for i in range(1, 4):\n            df_train[f'CV{i}_Fold'] = df_train[f'CV{i}_Fold'].astype(np.uint8)\n        df_test = df_all.loc[df_all['Type'] == 'Test', :].drop(columns=['Type', 'FVC', 'CV1_Fold', 'CV2_Fold', 'CV3_Fold'])    \n        \n        return df_train.copy(deep=True), df_test.reset_index(drop=True).copy(deep=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tabular_data_preprocessor = TabularDataPreprocessor(train=df_train,\n                                                    test=df_test, \n                                                    submission=df_submission, \n                                                    n_folds=2,\n                                                    shuffle=True,\n                                                    ohe=True,\n                                                    scale=False)\n\ndf_train, df_test = tabular_data_preprocessor.create_tabular_features()\n\nprint(f'Training Set (Tabular Features) Shape = {df_train.shape} - Patients = {df_train[\"Patient\"].nunique()}')\nprint(f'Training Set (Tabular Features) Memory Usage = {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set (Tabular Features) Shape = {df_test.shape} - Patients = {df_test[\"Patient\"].nunique()}')\nprint(f'Test Set (Tabular Features) Memory Usage = {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2.3. Segmentation Models"},{"metadata":{"trusted":true},"cell_type":"code","source":"class CoordinateConvolution3D(nn.Module):\n    \n    def __init__(self, in_channels, with_radius=False, **kwargs):\n        \n        super(CoordinateConvolution3D, self).__init__()\n        self.params = nn.Parameter(torch.ones((2)))\n        self.conv = nn.Conv2d(in_channels + 1, **kwargs)\n\n    def forward(self, X, Z):\n        \n        batch_size, channels, dim_x, dim_y = X.size()\n        zz_channel = torch.ones([1, 1, dim_x, dim_y])\n        \n        out = torch.cat([X, zz_channel.to(X.device) * Z.reshape(-1, 1, 1, 1) * self.params[0] ** 2 + self.params[1]], dim=1)\n        return self.conv(out)\n    \n\nclass SegmenterModel(nn.Module):\n    \n    def __init__(self, checkpoint=None):\n        \n        super().__init__()\n        self.model = smp.Unet('timm-efficientnet-b0', classes=3, encoder_weights=None, decoder_attention_type='scse')\n        \n        self.final = CoordinateConvolution3D(\n            self.model.segmentation_head[0].in_channels,\n            out_channels = self.model.segmentation_head[0].out_channels,\n            kernel_size = self.model.segmentation_head[0].kernel_size,\n            stride = self.model.segmentation_head[0].stride,\n            padding = self.model.segmentation_head[0].padding,\n            bias = True\n        )\n        self.model.segmentation_head = nn.Identity()\n        \n        if checkpoint is not None:\n            checkpoint = torch.load(checkpoint)['model_state_dict']\n            for k in list(checkpoint.keys()):\n                checkpoint[k[len('module.'):]] = checkpoint[k]\n                del checkpoint[k]\n            self.load_state_dict(checkpoint)\n        \n    def forward(self, x, z):\n        return self.final(self.model(x), z)\n\n\nclass PatientDataset(Dataset):\n\n    def __init__(self, scans, data_info_path):\n        \n        super().__init__()\n        self.scans = scans\n        self.data_info = np.load(data_info_path)\n        self.n_slices = len(self)\n        \n    def __getitem__(self, index):\n        \n        if index == 0:\n            scan = np.concatenate([self.scans[index:index + 2], np.expand_dims(self.scans[index + 1], axis=0)], axis=0)\n        elif index == self.n_slices - 1:\n            scan = np.concatenate([np.expand_dims(self.scans[index - 1], axis=0), self.scans[index - 1:index + 1]], axis=0)\n        else:\n            scan = np.array(self.scans[index - 1:index + 2])\n\n        scan = torch.FloatTensor((scan - self.data_info[0]) / self.data_info[1])\n        z_pos = torch.FloatTensor([float(index) / self.n_slices]) - 0.5\n        \n        return scan, z_pos\n\n    def __len__(self):\n        return len(self.scans)\n    \n    \nclass InferenceEngine:\n    \n    def __init__(self, batch_size, num_workers):\n        \n        if torch.cuda.is_available():\n            self.device = torch.device('cuda')\n        else:\n            self.device = torch.device('cpu')\n            \n        self.batch_size = batch_size\n        self.num_workers = num_workers\n        \n    def load_models(self):\n                \n        # Distilled Noisy Lung Trachea Segmentation Model\n        data_info_201_path = '../input/osic-pulmonary-fibrosis-datasets/Segmentation Models/201_data_info.npy'\n        model_201_noisy_lung_trachea_path = '../input/osic-pulmonary-fibrosis-datasets/Segmentation Models/201_noisy_lung_trachea_model.bin'\n        model_201_noisy_lung_trachea = SegmenterModel(model_201_noisy_lung_trachea_path).to(self.device)\n        model_201_noisy_lung_trachea.eval()\n\n        # Distilled Smooth Lung Heart Segmentation Model\n        data_info_202_path = '../input/osic-pulmonary-fibrosis-datasets/Segmentation Models/202_data_info.npy'\n        model_202_smooth_lung_heart_path = '../input/osic-pulmonary-fibrosis-datasets/Segmentation Models/202_smooth_lung_heart_model.bin'\n        model_202_smooth_lung_heart = SegmenterModel(model_202_smooth_lung_heart_path).to(self.device)\n        model_202_smooth_lung_heart.eval()\n\n        return model_201_noisy_lung_trachea, data_info_201_path, model_202_smooth_lung_heart, data_info_202_path\n        \n    def run_inference(self, scan, model, data_info_path, softmax):\n        \n        loader = DataLoader(PatientDataset(scan, data_info_path), \n                            batch_size=self.batch_size,\n                            num_workers=self.num_workers,\n                            shuffle=False,\n                            pin_memory=True)\n        \n        model.eval()      \n        predictions = []\n        \n        for input_scans, input_zpos in loader:\n            with torch.no_grad():\n                input_scans, input_zpos = input_scans.to(self.device), input_zpos.to(self.device)\n                outputs = model(input_scans, input_zpos)\n                \n                if softmax:\n                    outputs = F.softmax(outputs, dim=1, _stacklevel=5)\n                else:\n                    outputs = F.sigmoid(outputs)\n                    \n                predictions.append(outputs.detach().cpu())\n                \n        if len(scan) > self.batch_size:\n            predictions = torch.cat(predictions, dim=0)\n        else:\n            predictions = predictions[0]\n                \n        torch.cuda.empty_cache()\n        return predictions.numpy()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2.4. Image Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ImageDataPreprocessor:\n    \n    def __init__(self, train, test, resize_shape, window_min, window_max, y_min, y_max):\n        \n        self.train = train.copy(deep=True)\n        self.test = test.copy(deep=True)\n        \n        self.resize_shape = resize_shape\n        \n        self.window_min = window_min\n        self.window_max = window_max\n        self.y_min = y_min\n        self.y_max = y_max\n    \n    def resize(self, s, interpolate=False):\n        \n        \"\"\"\n        Resize slices to given size with nearest neighbor interpolation or linear interpolation\n\n        Parameters\n        ----------\n        s : numpy array, shape = (Rows, Columns)\n        2D numpy array of slice\n        \n        interpolate: boolean\n        If True, use linear interpolation\n        If False, use nearest neighbor interpolation\n\n        Returns\n        -------\n        s : numpy array, shape = (resize_shape[0], resize_shape[1])\n        2D numpy array after resized to resize_shape\n        \"\"\"\n        \n        if s.shape[0] != self.resize_shape[0] and s.shape[1] != self.resize_shape[1]:\n            return cv2.resize(s, self.resize_shape, interpolation=cv2.INTER_LINEAR if interpolate else cv2.INTER_NEAREST)\n\n        return s\n    \n    def remove_padding(self, v, masks=None):\n        \n        \"\"\"\n        Remove junk padding from scan if present\n\n        Parameters\n        ----------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volumes with padding\n        \n        masks: dictionary of numpy arrays, shape = (Depth, Rows, Columns)\n        3D numpy arrays of masks with padding\n\n        Returns\n        -------\n        v : numpy array, shape = (Depth, (Rows - All Constant Rows), (Columns - All Constant Columns))\n        3D numpy array after all constant rows and columns are dropped\n        \n        masks : dictionary of numpy arrays, shape = (Depth, (Rows - All Constant Rows), (Columns - All Constant Columns))\n        3D numpy arrays of masks after all constant rows and columns are dropped\n        \"\"\"\n\n        if v.shape[1] == self.resize_shape[0] and v.shape[2] == self.resize_shape[1]:\n            return v, masks\n\n        # Constant value is taken from the first voxel\n        constant_value = v[0, 0, 0]\n\n        # Remove all constant horizontal lines\n        mid_slice = v[v.shape[0] // 2]\n        indexer = mid_slice == constant_value\n        v = v[:, ~np.all(indexer, axis=1)]\n        for key in list(masks.keys()):\n            if masks[key] is None:\n                continue\n            masks[key] = masks[key][:, ~np.all(indexer, axis=1)]            \n        del mid_slice, indexer\n\n        # Remove all constant vertical lines \n        mid_slice = v[v.shape[0] // 2]\n        indexer = mid_slice == constant_value\n        v = v[:, :, ~np.all(indexer, axis=0)]\n        for key in list(masks.keys()):\n            if masks[key] is None:\n                continue\n            masks[key] = masks[key][:, :, ~np.all(indexer, axis=0)]\n        del mid_slice, indexer\n\n        return v, masks\n    \n    def padded_crop(self, v, btop, bbottom, bleft, bright, bg):\n        \n        v_padded = np.ones((v.shape[0], bbottom - btop, bright - bleft), dtype=v.dtype) * bg\n        \n        pad_top = 0 if btop > 0 else -btop\n        pad_bottom = 0 if bbottom < v.shape[1] else bbottom - v.shape[1]\n        pad_left = 0 if bleft > 0 else -bleft\n        pad_right = 0 if bright < v.shape[2] else bright - v.shape[2]\n        \n        v_padded[:, pad_top:v_padded.shape[1] - pad_bottom, pad_left:v_padded.shape[2] - pad_right,] = v[:, max(0, btop):min(v.shape[1], bbottom), max(0, bleft):min(v.shape[2], bright)]\n        return v_padded\n        \n    def find_spikes(self, v):\n        \n        \"\"\"\n        Find the spikes in volumes\n\n        Parameters\n        ----------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volume\n\n        Returns\n        -------\n        tuple of spikes in v\n        \"\"\"\n        \n        mean = v.mean()\n        std = v.std() * 4.5\n        try:\n            # Take values 4.5 sigma away from mean\n            # It is inside a try/except block because a scan in private test set throws error here\n            hist = np.histogram(v[(v < mean + std) & (v > mean - std)].flatten(), bins=50, density=True)\n            peaks, peak_properties = find_peaks(hist[0], distance=3, prominence=0.00015, rel_height=0.5)\n\n            return (hist[1][peaks[-1] - 1], hist[1][peaks[-1]], hist[1][peaks[-1] + 1], hist[1][peaks[-1] + 2])\n        except:\n            return (None, None, None, None)\n\n    def fix_hu(self, v, slope, intercept):\n        \n        \"\"\"\n        Fix HU values in 3D volume\n\n        Parameters\n        ----------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volume\n        \n        slope: int\n        Rescale Slope from DICOM metadata\n        \n        intercept: int\n        Rescale Intercept from DICOM metadata\n\n        Returns\n        -------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volume with shifted hu values\n        \"\"\"\n                \n        v = (v * slope + intercept).clip(-32768, 32767).astype(np.int16)        \n        a, _, _, d = self.find_spikes(v=v)\n        if a and d:\n            mid = (a + d) / 2\n            if abs(mid) > 512:\n                v -= np.int16(round(mid / 1024) * 1024)\n        return v\n    \n    def window(self, v, window_min, window_max, y_min, y_max):\n        \n        \"\"\"\n        Window HU values in 3D volume\n\n        Parameters\n        ----------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volume\n        \n        window_min: int       \n        window_max: int\n        y_min: int\n        y_max: int\n\n        Returns\n        -------\n        v : numpy array, shape = (Depth, Rows, Columns)\n        3D numpy array of volume with windowed hu values\n        \"\"\"\n            \n        v = ((v - window_min) / (window_max - window_min) * (y_max - y_min) + y_min).clip(0, 255).astype(np.uint8)\n        return v\n        \n    def load_scan(self, scan_path, masks=None):\n        \n        \"\"\"\n        Preprocess and load the given patients scan\n\n        Parameters\n        ----------\n        scan_path : str Path of the CT scan\n\n        Returns\n        -------\n        scan : numpy array, shape = (n_slices, self.resize_shape[0], self.resize_shape[1])\n        numpy array after the cropped and resized slices are stacked\n        metadata : dict\n        dictionary of processed metadata        \n        \"\"\"\n\n        # Load scan from the given scan path\n        files_path = glob(f'{scan_path}/*.dcm')\n        patient_directory = [pydicom.dcmread(f) for f in files_path]\n        \n        if masks is not None:\n            pass\n        else:\n            masks = {\n                'noisy_lung': None,\n                'smooth_lung': None,\n                'heart': None,\n                'trachea': None\n            }\n        \n        try:\n            # Sort slices by ImagePositionPatient Z if the field exists\n            patient_directory.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n            slice_positions = np.round([s.ImagePositionPatient[2] for s in patient_directory], 4)\n            non_duplicate_idx = np.unique([np.where(slice_position == slice_positions)[0][0] for slice_position in slice_positions])\n            del slice_positions\n        except AttributeError:\n            # If ImagePositionPatient doesn't exist, sort them by InstanceNumber\n            # (Every CT scan has ImagePositionPatient in private test set)\n            patient_directory.sort(key=lambda s: int(s.InstanceNumber))\n            instance_numbers = np.array([int(s.InstanceNumber) for s in patient_directory])\n            non_duplicate_idx = np.unique([np.where(instance_number == instance_numbers)[0][0] for instance_number in instance_numbers])\n            del instance_numbers\n\n        keys = list(masks.keys())\n        for key in keys:\n            if not masks[key]:\n                # Explicitly set empty array to None\n                masks[key] = None\n                continue\n                \n            mask, _ = nrrd.read(masks[key][0])\n            mask = mask.transpose(2, 1, 0).astype(np.uint8)\n            \n            assert mask.shape[0] == len(patient_directory)\n            assert (mask.min() == 0) and (mask.max() == 1)\n            \n            masks[key] = (mask > 0).astype(np.uint8)\n            \n        # Order of operations\n        if masks['noisy_lung'] is not None and masks['smooth_lung'] is not None:\n            # All noisy lung must be contained in non-noisy lung\n            try:\n                masks['smooth_lung'][masks['noisy_lung'] > 0] = 1\n            except:\n                print(f\"{scan_path} smooth_lung mask != heart mask, {masks['smooth_lung'].shape}, {masks['heart'].shape}\")\n                return\n            \n        if masks['noisy_lung'] is not None and masks['trachea'] is not None:\n            # Noisy lung must have trachea removed\n            try:\n                masks['noisy_lung'][masks['trachea'] > 0] = 0\n            except:\n                print(f\"{scan_path} trachea mask != noisy lung mask, {masks['trachea'].shape}, {masks['noisy_lung'].shape}\")\n                return\n            \n        if masks['heart'] is not None and masks['smooth_lung'] is not None:\n            # Heart must have smooth lung removed\n            try:\n                masks['heart'][masks['smooth_lung'] > 0] = 0\n            except:\n                print(f\"{scan_path} smooth_lung mask != heart mask, {masks['smooth_lung'].shape}, {masks['heart'].shape}\")\n                return\n            \n        # Slice indices with the same ImagePositionPatient Z or InstanceNumber are removed from scan and masks\n        patient_directory = np.array(patient_directory)[non_duplicate_idx]\n        for key in keys:\n            if masks[key] is None:\n                continue\n            masks[key] = masks[key][non_duplicate_idx]\n        del non_duplicate_idx\n\n        # Slice indices of empty slices are removed from scan and masks\n        non_empty_idx = [i for i, s in enumerate(patient_directory) if not np.all(s.pixel_array == s.pixel_array[0, 0], axis=(-1, -2))]\n        patient_directory = patient_directory[non_empty_idx]\n        for key in keys:\n            if masks[key] is None:\n                continue\n            masks[key] = masks[key][non_empty_idx]\n        del non_empty_idx\n            \n        # DICOM Metadata\n        metadata = {}\n        \n        metadata['RescaleSlope'] = patient_directory[0].RescaleSlope\n        metadata['RescaleIntercept'] = patient_directory[0].RescaleIntercept\n        \n        pixel_spacings = np.zeros((len(patient_directory), 2))\n        slice_positions = np.zeros((len(patient_directory)))\n        slice_thickness = np.zeros((len(patient_directory)))\n        \n        for i, s in enumerate(patient_directory): \n            try:\n                pixel_spacings[i] = np.array(s.PixelSpacing)\n            except AttributeError:\n                pixel_spacings[i] = np.nan\n                \n            try:\n                slice_positions[i] = s.ImagePositionPatient[2]\n            except AttributeError:\n                slice_positions[i] = np.nan\n                \n            try:\n                slice_thickness[i] = float(s.SliceThickness)\n            except AttributeError:\n                slice_thickness[i] = np.nan\n        \n        # Fill missing pixel spacings and slice thinknesses with nanmean\n        pixel_spacings[np.isnan(pixel_spacings[:, 0]), 0] = np.nanmean(pixel_spacings[:, 0])\n        pixel_spacings[np.isnan(pixel_spacings[:, 1]), 1] = np.nanmean(pixel_spacings[:, 1])\n        slice_thickness[np.isnan(slice_thickness)] = np.nanmean(slice_thickness)\n        \n        # All PixelSpacing values are same in a CT scan but they have to be rounded and averaged along the axis 0 for floating point precision errors\n        metadata['PixelSpacing'] = list(np.round(pixel_spacings.mean(axis=0), 3))\n        del pixel_spacings\n        \n        # Two of the patients don't have their SpacingBetweenSlices or ImagePositionPatient fields so it is impossible to retrieve their SliceSpacings\n        # Their SliceSpacings are filled with the most common values of similar CT scans\n        # Most common rounded SliceSpacing value is selected for rest of the patients\n        patient_name = scan_path.split('/')[-1]\n        if patient_name == 'ID00128637202219474716089':\n            metadata['SliceSpacing'] = np.array([5.0] * len(patient_directory))\n        elif patient_name == 'ID00132637202222178761324':\n            metadata['SliceSpacing'] = np.array([0.7] * len(patient_directory))\n        else:\n            metadata['SliceSpacing'] = np.abs(np.diff(np.round(slice_positions, 3)))\n            metadata['SliceSpacing'] = np.concatenate([metadata['SliceSpacing'], [metadata['SliceSpacing'][-1]]])\n            \n            nan_mask = np.isnan(metadata['SliceSpacing'])\n            if metadata['SliceSpacing'].shape[0] == nan_mask.sum():\n                # All SliceSpacing values are NaNs\n                # Set SliceSpacing by filling values based on slice count\n                # This power function is regresses from known slice count - slice spacing pairs:\n                metadata['SliceSpacing'] = np.array([275.79 * (len(patient_directory) ** (-0.987))] * len(patient_directory))\n                \n                # Sandor's suggestions RE the above (override!):\n                # If slice count < 50,  slicespacing is 5 mm (99%)\n                # If slice  count > 150, slicespacing is 1 mm (99%)\n                if len(patient_directory) < 50:\n                    metadata['SliceSpacing'] = np.array([5.0] * len(patient_directory))\n                elif len(patient_directory) > 150:\n                    metadata['SliceSpacing'] = np.array([1.0] * len(patient_directory))\n                # Between 50 and 150 no good estimation can be given... i would go with 5 mm up to 75 and with 3 mm to 150\n                \n            elif nan_mask.sum() > 0:\n                # Fill intermittant NAs with the mean value\n                metadata['SliceSpacing'][nan_mask] = np.nanmean(metadata['SliceSpacing'])\n            del nan_mask\n        del slice_positions\n                \n        metadata['SliceThickness'] = slice_thickness\n        del slice_thickness\n        \n        # Remove padding\n        scan, masks = self.remove_padding(v=np.array([s.pixel_array for s in patient_directory], dtype=np.int16), masks=masks)\n        del patient_directory\n        \n        # Ensure scans and all available all have the same dimensions\n        for key in keys:\n            if masks[key] is not None:\n                assert masks[key].shape == scan.shape\n                \n        scan = self.fix_hu(v=scan, slope=metadata['RescaleSlope'], intercept=metadata['RescaleIntercept'])\n        \n        # Center crop 325 mm around the barycenter where values > the mean\n        # Cropped 335 mm (13.2 inches) in the x-y direction \n        center = np.argwhere(scan > scan.mean()).mean(axis=0).astype(np.int16)[1:]\n        bounding_box = (335 / np.array(metadata['PixelSpacing'])).astype(np.int16)\n        \n        metadata['original_scan_shape'] = scan.shape\n        metadata['btop'] = center[0] - bounding_box[0] // 2\n        metadata['bbottom'] = center[0] + bounding_box[0] // 2\n        metadata['bleft'] = center[1] - bounding_box[1] // 2\n        metadata['bright'] = center[1] + bounding_box[1] // 2\n        metadata['mm_crop'] = 335\n        del center, bounding_box\n        \n        scan = self.padded_crop(v=scan, btop=metadata['btop'], bbottom=metadata['bbottom'], bleft=metadata['bleft'], bright=metadata['bright'], bg=scan.min())\n        for key in keys:\n            if masks[key] is not None:                \n                masks[key] = self.padded_crop(v=masks[key], btop=metadata['btop'], bbot=metadata['bbottom'], bleft=metadata['bleft'], bright=metadata['bright'], bg=0)\n        \n        # Standardize the orientation of the patient\n        body = scan < -450\n        if body[:int(scan.shape[0] * 0.5)].sum() > body[int(scan.shape[0] * 0.5):].sum():\n            # Patient was put into the machine backwards\n            scan = scan[::-1]\n            for key in keys:\n                if masks[key] is None: continue\n                masks[key] = masks[key][::-1]\n        del body\n\n        # Bone at the 99th percentile\n        threshold = np.percentile(scan.flatten(), 99)\n        bone = (scan[2:-2,\n                     int(scan.shape[1] * 0.17):int(scan.shape[1] * 0.83):, # Cut the table out from top/bottom\n                     int(scan.shape[2] * 0.37):int(scan.shape[2] * 0.63):, # Limit the the center slices where the back bone should be\n                     ] >= threshold).mean(axis=0)\n        bone = bone.mean(axis=1)\n        if bone[:int(bone.shape[0] * 0.5)].sum() > bone[int(bone.shape[0] * 0.5):].sum():\n            # Patient is lying on their stomach\n            scan = scan[:, ::-1]\n            for key in keys:\n                if masks[key] is None: continue\n                masks[key] = masks[key][:, ::-1]                \n        del threshold, bone\n    \n        # Rescale each one to 512x512 or whatever\n        vscan = np.zeros((scan.shape[0], self.resize_shape[0], self.resize_shape[0]), dtype=np.int16)\n        vmasks = {\n            k: None if m is None else np.zeros((m.shape[0], self.resize_shape[0], self.resize_shape[0]), dtype=np.uint8)\n            for k,m in masks.items()\n        }\n        \n        for i, s in enumerate(scan):\n            vscan[i] = self.resize(s=s, interpolate=True)\n            for key in keys:\n                if vmasks[key] is None: continue\n                vmasks[key][i] = self.resize(s=masks[key][i], interpolate=False)\n\n        vscan = self.window(v=vscan, window_min=self.window_min, window_max=self.window_max, y_min=self.y_min, y_max=self.y_max)\n        return {'scan': vscan, **metadata, **vmasks}\n    \n    def find_largest_diameter(self, mask, mm_per_pixel):        \n        first = np.argmax(mask, axis=2)\n        last = mask.shape[2] - np.argmax(mask[:, :, ::-1], axis=2)\n        last[last == mask.shape[2]] = 0\n        first[first == 0] = mask.shape[2]\n        return (last - first).max() * mm_per_pixel\n    \n    def find_volume(self, mask, mm_per_pixel, mm_slice_spacing):\n        return (mask.sum(axis=(1, 2)) * mm_per_pixel * mm_per_pixel * mm_slice_spacing).sum() / 1000\n    \n    def find_surface_area(self, mask, mm_per_pixel, mm_slice_spacing):        \n        surface_area = 0\n        for idx, c in enumerate(mask):\n            depth = mm_slice_spacing[idx] * np.sqrt(2)\n            contours = cv2.findContours(c, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)[0]\n            surface_area += sum([cv2.arcLength(cnt, closed=True) * mm_per_pixel * depth for cnt in contours])\n\n        return surface_area\n    \n    def create_image_features(self):\n        \n        # Segmentation Models\n        # ------------------------\n        inference_engine = InferenceEngine(batch_size=32, num_workers=0)\n        model_201_noisy_lung_trachea, data_info_201, model_202_smooth_lung_heart, data_info_202 = inference_engine.load_models()\n        \n        # TRAINING SET\n        # ------------\n        # Set this to True for creating training set features on live kernel session\n        # Download features after creating them and upload them on a dataset\n        # Otherwise, GPU run-time exceeds 4 hours\n        if False:\n            for i, patient in enumerate(self.train['Patient'].unique(), 1):\n                \n                # Skipping 1018 slice CT scan because Kaggle kernel RAM is not enough for inference\n                # Created features are filled with mean for test purposes\n                # In final submissions, they should be switched with real values\n                if patient == 'ID00078637202199415319443':\n                    print(f'[{i}/{len(self.train[\"Patient\"].unique())}] Training Set - Patient: {patient} skipped\\n')\n                    continue\n                    \n                # --------------------\n                # Create features here\n                # --------------------\n            \n            features = ['smooth_lung_mean', 'smooth_lung_std', 'smooth_lung_skew',\n                        'chest_diameter', 'smooth_lung_volume', 'smooth_lung_surface_area', 'smooth_lung_lower_level_calcification',\n                        'noisy_lung_volume', 'noisy_lung_surface_area',\n                        'heart_diameter', 'heart_volume', 'heart_surface_area',\n                        'trachea_diameter', 'trachea_branching_angle']\n                            \n            self.train[(['Patient'] + features)].to_csv('train_features.csv', index=False)\n            print('Training features are created and saved\\n')\n            return\n\n        else:\n            # Merging precomputed training set features from dataset\n            df_train_features = pd.read_csv('../input/osic-pulmonary-fibrosis-datasets/train_features.csv').rename(columns={'patient_name': 'Patient'})\n            features = ['smooth_lung_mean', 'smooth_lung_std', 'smooth_lung_skew', 'smooth_lung_kurt',\n                        'chest_diameter', 'smooth_lung_volume', 'smooth_lung_surface_area', 'smooth_lung_lower_level_calcification',\n                        'noisy_lung_volume', 'noisy_lung_surface_area',\n                        'heart_diameter', 'heart_volume', 'heart_surface_area',\n                        'trachea_diameter', 'trachea_branching_angle']\n            \n            self.train = self.train.merge(df_train_features[(['Patient'] + features)], how='left', on='Patient')\n            print('Merged precomputed features to training set\\n')\n        \n        # TEST SET\n        # ------------\n        # Computing test set features on run-time\n        for i, patient in enumerate(self.test['Patient'].unique(), 1):\n            \n            results = self.load_scan(scan_path=f'../input/osic-pulmonary-fibrosis-progression/test/{patient}', masks=None)\n            mm_per_pixel = results['mm_crop'] / results['scan'].shape[1]\n            mm_slice_spacing = results['SliceSpacing']\n\n            mask_noisy_lung_trachea = inference_engine.run_inference(results['scan'], model_201_noisy_lung_trachea, data_info_201, softmax=True)\n            mask_smooth_lung_heart = inference_engine.run_inference(results['scan'], model_202_smooth_lung_heart, data_info_202, softmax=True)\n\n            mask_smooth_lung = (mask_smooth_lung_heart[:, 1] > 0.6).astype(np.uint8)\n            mask_heart = (mask_smooth_lung_heart[:, 2] > 0.6).astype(np.uint8)\n            mask_noisy_lung = (mask_noisy_lung_trachea[:, 1] > 0.6).astype(np.uint8)\n            mask_trachea = (mask_noisy_lung_trachea[:, 2] > 0.6).astype(np.uint8)\n            del mask_smooth_lung_heart, mask_noisy_lung_trachea\n\n            mask_smooth_lung[mask_trachea > 0] = 0 # Voxels predicted as trachea can't be smooth lung\n            mask_noisy_lung[mask_smooth_lung == 0] = 0 # Voxels predicted as smooth lung are overwritten to noisy lung\n            \n            print(f'[{i}/{len(self.test[\"Patient\"].unique())}] Test Set - Patient: {patient} - Scan Shape: {results[\"scan\"].shape} ({results[\"scan\"].nbytes >> 20}MB)')\n            print(f'{\" \" * 6}Smooth Lung Mask Shape: {mask_smooth_lung.shape} ({mask_smooth_lung.nbytes >> 20}MB) - Heart Mask Shape: {mask_heart.shape} ({mask_heart.nbytes >> 20}MB)')\n            print(f'{\" \" * 6}Noisy Lung  Mask Shape: {mask_noisy_lung.shape} ({mask_noisy_lung.nbytes >> 20}MB) - Trachea Mask Shape: {mask_trachea.shape} ({mask_trachea.nbytes >> 20}MB)\\n')\n            \n            # 1. Smooth Lung Statistical Features\n            # -----------------------------------\n            smooth_lung = results['scan'][mask_smooth_lung > 0]\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_mean'] = smooth_lung.mean()\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_std']  = smooth_lung.std()\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_skew'] = skew(smooth_lung)\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_kurt'] = kurtosis(smooth_lung)\n            del smooth_lung\n                \n            # 2. Smooth Lung Volume Features\n            # ------------------------------\n            if np.all(0 == mask_smooth_lung):\n                z_min, z_max = mask_smooth_lung.shape[0] // 20, mask_smooth_lung.shape[0] // 20 * 19\n                x_min, x_max = mask_smooth_lung.shape[1] // 20, mask_smooth_lung.shape[1] // 20 * 19\n                y_min, y_max = mask_smooth_lung.shape[2] // 20, mask_smooth_lung.shape[2] // 20 * 19\n            else:\n                zs, xs, ys = np.where(mask_smooth_lung != 0)\n                z_min, z_max = min(zs), max(zs) + 1\n                x_min, x_max = min(xs), max(xs) + 1\n                y_min, y_max = min(ys), max(ys) + 1\n                del zs, xs, ys\n                    \n            self.test.loc[self.test['Patient'] == patient, 'chest_diameter'] = self.find_largest_diameter(mask_smooth_lung, mm_per_pixel)\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_volume'] = self.find_volume(mask_smooth_lung, mm_per_pixel, mm_slice_spacing) \n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_surface_area'] = self.find_surface_area(mask_smooth_lung, mm_per_pixel, mm_slice_spacing)\n\n            # 3. Smooth Lung Lower Z-Level Calcification\n            # ------------------------------------------\n            last = z_min + max(1, (z_max - z_min) // 3)\n            scan_calcium = results['scan'][z_min:last]\n            mask_calcium = mask_smooth_lung[z_min:last].copy()\n            pixels = max(2, int(round(3 * mm_per_pixel)))\n\n            kernel_c = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (pixels * 2, pixels * 2))\n            for i, s in enumerate(scan_calcium):\n                # Errode mask a specific number of mm...    \n                mask_calcium[i] = cv2.erode(mask_calcium[i], kernel_c, iterations=1)\n\n            # Take only masked values > threshold\n            mask_calcium = ((scan_calcium > 800) * mask_calcium).astype(np.uint8)\n            del scan_calcium\n\n            # Denoise by getting rid of small calcification sprinkles\n            kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))\n\n            calcification_count = 0\n            for i in range(mask_calcium.shape[0]):\n                # Now, run connected components in 2d\n                calcification_count += cv2.connectedComponents(mask_calcium[i])[0] - 1\n\n            del mask_calcium, kernel, kernel_c, last, mask_smooth_lung\n            self.test.loc[self.test['Patient'] == patient, 'smooth_lung_lower_level_calcification'] = calcification_count / (z_max - z_min)\n\n            # 4. Noisy Lung Volume Features\n            # -----------------------------\n            self.test.loc[self.test['Patient'] == patient, 'noisy_lung_volume'] = self.find_volume(mask_noisy_lung, mm_per_pixel, mm_slice_spacing)\n            self.test.loc[self.test['Patient'] == patient, 'noisy_lung_surface_area'] = self.find_surface_area(mask_noisy_lung, mm_per_pixel, mm_slice_spacing)\n\n            # 5. Heart Volume Features\n            # ------------------------\n            mask_heart = mask_heart[z_min:z_max, x_min:x_max, y_min:y_max]\n            self.test.loc[self.test['Patient'] == patient, 'heart_diameter'] = self.find_largest_diameter(mask_heart, mm_per_pixel)\n            self.test.loc[self.test['Patient'] == patient, 'heart_volume'] = self.find_volume(mask_heart, mm_per_pixel, mm_slice_spacing[z_min:z_max])\n            self.test.loc[self.test['Patient'] == patient, 'heart_surface_area'] = self.find_surface_area(mask_heart, mm_per_pixel, mm_slice_spacing[z_min:z_max])\n            del mask_heart\n            \n            # 6. Trachea Volume Features\n            # --------------------------\n            # For trachea features, look within the smooth lung bounding box region\n            trachea = mask_trachea[z_min:z_max, x_min:x_max, y_min:y_max]\n            center_a, center_b = None, None\n            \n            for idx in range(trachea.shape[0] - 2, -1, -1):\n                num_labels, labels = cv2.connectedComponents(trachea[idx])\n                label_vals, counts = np.unique(labels, return_counts=True)\n\n                if num_labels < 3:\n                    continue\n\n                if (counts > 16).sum() < 3:\n                    continue\n\n                count_sort = np.argsort(counts)[::-1] # Descending\n                idx_bifurcation = idx\n                idx += 1\n\n                remove_16mm = mm_slice_spacing[z_min + idx:].cumsum() > 16\n                if remove_16mm.sum()==0:\n                    remove_16mm = remove_16mm.shape[0]\n                else:\n                    remove_16mm = np.nonzero(remove_16mm)[0][0] + 1\n\n                idx_single_trunk = idx + remove_16mm\n                if idx_single_trunk >= trachea.shape[0]:\n                    idx_single_trunk = trachea.shape[0] - 1\n\n                center_a = np.argwhere(labels==label_vals[count_sort[1]]).mean(axis=0)\n                center_b = np.argwhere(labels==label_vals[count_sort[2]]).mean(axis=0)\n                break\n                \n            if center_a is None:\n                self.test.loc[self.test['Patient'] == patient, 'trachea_diameter_normalized'] = np.nan\n                self.test.loc[self.test['Patient'] == patient, 'trachea_branching_angle'] = np.nan\n            else:\n                # -----------------------------------------------------------\n                # Feature: Trachea diameter before branching / chest diameter\n                largest_trachea = trachea[idx_single_trunk:]\n\n                # If we don't find a trachea in here, we need to do something...\n                # TODO: ...\n\n                # This will be normalized with chest diameter...\n                self.test.loc[self.test['Patient'] == patient, 'trachea_diameter'] = self.find_largest_diameter(largest_trachea, mm_per_pixel)\n                del largest_trachea\n\n                # --------------------------------    \n                # Feature: Trachea branching angle\n                center_c = np.argwhere(trachea[idx_single_trunk]==1).mean(axis=0)\n                v0 = np.array([center_c[0], center_c[1], idx_single_trunk])\n                v1 = np.array([center_a[0], center_a[1], idx_bifurcation]) - v0\n                v2 = np.array([center_b[0], center_b[1], idx_bifurcation]) - v0\n                v1_u = v1 / np.linalg.norm(v1)\n                v2_u = v2 / np.linalg.norm(v2)\n\n                self.test.loc[self.test['Patient'] == patient,'trachea_branching_angle'] = np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))\n                self.test.loc[self.test['Patient'] == patient,'trachea_branching_angle'] *= 180 / np.pi # Degrees\n\n        return self.train.copy(deep=True), self.test.copy(deep=True)\n            \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_data_preprocessor = ImageDataPreprocessor(\n    train=df_train,\n    test=df_test,\n    resize_shape=(512, 512),\n    window_min=-1024, window_max=250,\n    y_min=0, y_max=255\n)\n\ndf_train, df_test = image_data_preprocessor.create_image_features()\n\nprint(f'\\nTraining Set (Tabular Features + Image Features) Shape = {df_train.shape} - Patients = {df_train[\"Patient\"].nunique()}')\nprint(f'Training Set (Tabular Features + Image Features) Memory Usage = {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set (Tabular Features + Image Features) Shape = {df_test.shape} - Patients = {df_test[\"Patient\"].nunique()}')\nprint(f'Test Set (Tabular Features + Image Features) Memory Usage = {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Models"},{"metadata":{},"cell_type":"markdown","source":"### 2.1. MLP Quantile Regression"},{"metadata":{"trusted":true},"cell_type":"code","source":"class QuantileRegressorMLP:\n    \n    def __init__(self, model, cv, predictors, mlp_parameters, qr_parameters, hybrid_parameters):\n        \n        self.model = model\n        self.cv = cv\n        self.predictors = predictors\n\n        self.mlp_parameters = mlp_parameters\n        self.qr_parameters = qr_parameters\n        self.hybrid_parameters = hybrid_parameters\n        \n    def laplace_log_likelihood_metric(self, y_true, y_pred, sigma):\n        \n        sigma_clipped = np.maximum(sigma, 70)\n        delta_clipped = np.minimum(np.abs(y_true - y_pred), 1000)\n        score = - np.sqrt(2) * delta_clipped / sigma_clipped - np.log(np.sqrt(2) * sigma_clipped)\n\n        return np.mean(score)        \n            \n    def laplace_log_likelihood_loss(self, y_true, y_pred):\n                \n        K.cast(y_true, 'float32')\n        K.cast(y_pred, 'float32')\n    \n        sigma_lower_bound = K.constant(70, dtype='float32')\n        delta_upper_bound = K.constant(1000, dtype='float32')\n\n        sigma = y_pred[:, 1]\n        fvc_pred = y_pred[:, 0]\n\n        sigma_clipped = K.maximum(sigma, sigma_lower_bound)\n        delta = K.abs(y_true[:, 0] - fvc_pred)\n        delta_clipped = K.minimum(delta, delta_upper_bound)\n\n        score = (delta_clipped / sigma_clipped) * K.sqrt(K.cast(2, 'float32')) + K.log(sigma_clipped * K.sqrt(K.cast(2, 'float32')))\n        return K.mean(score)\n    \n    def quantile_loss(self, y_true, y_pred):\n        \n        quantiles = K.constant(np.array([self.qr_parameters['quantiles']]), dtype='float32')        \n        error = y_true - y_pred\n        return K.mean(K.maximum(quantiles * error, (quantiles - 1) * error))\n    \n    def hybrid_loss(self, w):\n        \n        def loss(y_true, y_pred):\n            return w * self.quantile_loss(y_true, y_pred) + (1 - w) * self.laplace_log_likelihood_loss(y_true, y_pred)\n        return loss\n    \n    def get_model(self, input_shape, m):\n        \n        model = None\n        \n        if m == 'MLP':            \n            input_layer = Input(shape=input_shape)\n            x = Dense(2 ** 7, activation='relu')(input_layer)\n            x = GaussianDropout(0.01)(x)\n            x = Dense(2 ** 7, activation='relu')(x)\n            x = GaussianDropout(0.01)(x)\n            output_layer = Dense(2, activation='relu')(x)\n\n            model = Model(input_layer, output_layer)\n            model.compile(loss=self.laplace_log_likelihood_loss, optimizer=tf.keras.optimizers.Adam(lr=self.mlp_parameters['lr']), metrics=[self.laplace_log_likelihood_loss])\n            \n        elif m == 'QR':            \n            input_layer = Input(shape=input_shape) \n            x = Dense(2 ** 6, activation='relu')(input_layer)   \n            x = GaussianDropout(0.01)(x)\n            x = Dense(2 ** 6, activation='relu')(x)\n            x = GaussianDropout(0.01)(x)\n            output_layer = Dense(3, activation='relu')(x)\n\n            model = Model(input_layer, output_layer)\n            model.compile(loss=self.quantile_loss, optimizer=Adam(lr=self.qr_parameters['lr']), metrics=[self.laplace_log_likelihood_loss])\n            \n        elif m == 'HYBRID':\n            input_layer = Input(shape=input_shape) \n            x = Dense(2 ** 6, activation='relu')(input_layer)   \n            x = GaussianDropout(0.01)(x)\n            x = Dense(2 ** 6, activation='relu')(x)\n            x = GaussianDropout(0.01)(x)\n            output_layer = Dense(3, activation='relu')(x)\n\n            model = Model(input_layer, output_layer)\n            model.compile(loss=self.hybrid_loss(0.5), optimizer=Adam(lr=self.qr_parameters['lr']), metrics=[self.laplace_log_likelihood_loss])\n\n        return model\n        \n    def train(self, X_train, y_train):\n        \n        self.mlp_scores_all = []\n        self.qr_scores_all = []\n        self.mlp_scores_last = []\n        self.qr_scores_last = []\n        self.hybrid_scores_all = []\n        self.hybrid_scores_last = []\n\n        self.mlp_oof = pd.DataFrame(np.zeros((len(y_train), 2)))\n        self.qr_oof = pd.DataFrame(np.zeros((len(y_train), len(self.qr_parameters['quantiles']))))\n        self.hybrid_oof = pd.DataFrame(np.zeros((len(y_train), len(self.hybrid_parameters['quantiles']))))\n            \n        self.mlp_models = {'CV1': [], 'CV2': [], 'CV3': []}\n        self.qr_models = {'CV1': [], 'CV2': [], 'CV3': []}\n        self.hybrid_models = {'CV1': [], 'CV2': [], 'CV3': []}\n        \n        models = [self.model] if self.model != 'Stack' else ['MLP', 'QR', 'HYBRID']        \n        for m in models:            \n            print(f'\\nRunning {m.upper()} Model\\n{(\"-\") * (14 + (len(m)))}')\n            \n            for cv in self.cv:               \n                for fold in sorted(X_train[f'CV{cv}_Fold'].unique()):\n                                        \n                    trn_idx, val_idx = X_train.loc[X_train[f'CV{cv}_Fold'] != fold].index, X_train.loc[X_train[f'CV{cv}_Fold'] == fold].index\n\n                    X_trn = X_train.loc[trn_idx, self.predictors]\n                    y_trn = y_train.loc[trn_idx]\n                    X_val = X_train.loc[val_idx, self.predictors]\n                    y_val = y_train.loc[val_idx]\n\n                    model = self.get_model(input_shape=X_trn.shape[1], m=m)                \n                    if m ==  'MLP':                    \n                        model.fit(X_trn, y_trn, epochs=self.mlp_parameters['epochs'], batch_size=self.mlp_parameters['batch_size'], verbose=0)\n                        self.mlp_models[f'CV{cv}'].append(model)                    \n                    elif m == 'QR':                    \n                        model.fit(X_trn, y_trn, epochs=self.qr_parameters['epochs'], batch_size=self.qr_parameters['batch_size'], verbose=0)\n                        self.qr_models[f'CV{cv}'].append(model)\n                    elif m == 'HYBRID':                    \n                        model.fit(X_trn, y_trn, epochs=self.hybrid_parameters['epochs'], batch_size=self.hybrid_parameters['batch_size'], verbose=0)\n                        self.hybrid_models[f'CV{cv}'].append(model)\n\n                    predictions = model.predict(X_val)\n                    \n                    if m == 'MLP':                    \n                        oof_predictions = predictions[:, 0]\n                        self.mlp_oof.iloc[val_idx, 0] = oof_predictions\n                        df_train.loc[val_idx, f'CV{cv}_MLP_FVC_Predictions'] = oof_predictions\n\n                        oof_confidence = predictions[:, 1]\n                        self.mlp_oof.iloc[val_idx, 1] = oof_confidence\n                        df_train.loc[val_idx, f'CV{cv}_MLP_Confidence_Predictions'] = oof_confidence\n                        \n                        fold_final_scores = []\n                        for df_patient in np.array_split(df_train.loc[val_idx].groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train.loc[val_idx, 'Patient'].nunique()):\n                            fold_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_MLP_FVC_Predictions'], df_patient[f'CV{cv}_MLP_Confidence_Predictions']))\n                        \n                    elif m == 'QR':                    \n                        oof_predictions = predictions[:, 1]\n                        oof_confidence = predictions[:, 2] - predictions[:, 0]\n                        for i, quantile in enumerate(self.qr_parameters['quantiles']):\n                            self.qr_oof.iloc[val_idx, i] = predictions[:, i]\n                            df_train.loc[val_idx, f'CV{cv}_QR_{quantile}_Predictions'] = predictions[:, i]\n                            \n                        fold_final_scores = []\n                        for df_patient in np.array_split(df_train.loc[val_idx].groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train.loc[val_idx, 'Patient'].nunique()):\n                            fold_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][1]}_Predictions'], (df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][2]}_Predictions'] - df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][0]}_Predictions'])))\n                            \n                    elif m == 'HYBRID':\n                        oof_predictions = predictions[:, 1]\n                        oof_confidence = predictions[:, 2] - predictions[:, 0]\n                        for i, quantile in enumerate(self.hybrid_parameters['quantiles']):\n                            self.hybrid_oof.iloc[val_idx, i] = predictions[:, i]\n                            df_train.loc[val_idx, f'CV{cv}_HYBRID_{quantile}_Predictions'] = predictions[:, i]\n                            \n                        fold_final_scores = []\n                        for df_patient in np.array_split(df_train.loc[val_idx].groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train.loc[val_idx, 'Patient'].nunique()):\n                            fold_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][1]}_Predictions'], (df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][2]}_Predictions'] - df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][0]}_Predictions'])))\n\n                    oof_score_all = self.laplace_log_likelihood_metric(y_val, oof_predictions, oof_confidence)\n                    print(f'CV {cv} {m} Fold {int(fold)} - X_train: {X_trn.shape} X_val: {X_val.shape} - All Measurement Score: {oof_score_all:.6} - Final 3 Measurement Score {np.mean(fold_final_scores):.6} [Std: {np.std(fold_final_scores):.6}]')\n\n                oof_final_scores = []\n                \n                if m == 'MLP':\n                    for df_patient in np.array_split(df_train.groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train['Patient'].nunique()):\n                        oof_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_MLP_FVC_Predictions'], df_patient[f'CV{cv}_MLP_Confidence_Predictions']))\n                    oof_all_score = self.laplace_log_likelihood_metric(y_train, self.mlp_oof.iloc[:, 0], self.mlp_oof.iloc[:, 1])                    \n                elif m == 'QR': \n                    for df_patient in np.array_split(df_train.groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train['Patient'].nunique()):\n                        oof_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][1]}_Predictions'], (df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][2]}_Predictions'] - df_patient[f'CV{cv}_QR_{self.qr_parameters[\"quantiles\"][0]}_Predictions'])))\n                    oof_all_score = self.laplace_log_likelihood_metric(y_train, self.qr_oof.iloc[:, 1], (self.qr_oof.iloc[:, 2] - self.qr_oof.iloc[:, 0]))\n                elif m == 'HYBRID': \n                    for df_patient in np.array_split(df_train.groupby('Patient').nth([-1, -2, -3]).reset_index(), df_train['Patient'].nunique()):\n                        oof_final_scores.append(self.laplace_log_likelihood_metric(df_patient['FVC'], df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][1]}_Predictions'], (df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][2]}_Predictions'] - df_patient[f'CV{cv}_HYBRID_{self.hybrid_parameters[\"quantiles\"][0]}_Predictions'])))\n                    oof_all_score = self.laplace_log_likelihood_metric(y_train, self.hybrid_oof.iloc[:, 1], (self.hybrid_oof.iloc[:, 2] - self.hybrid_oof.iloc[:, 0]))\n                    \n                print(f'{\"-\" * 30}\\nCV {cv} {m} All Measurement OOF Score {oof_all_score:.6} - Final 3 Measurement OOF Score {np.mean(oof_final_scores):.6} [Std: {np.std(oof_final_scores):.6}]\\n{\"-\" * 30}\\n')\n                    \n    def predict(self, X_test):\n        \n        for cv in self.cv:\n            mlp_predictions = np.zeros((len(X_test), 2))            \n            for model in self.mlp_models[f'CV{cv}']:\n                mlp_predictions += model.predict(X_test[self.predictors]) / len(self.mlp_models[f'CV{cv}'])\n    \n            X_test[f'CV{cv}_MLP_FVC_Predictions'] = mlp_predictions[:, 0]\n            X_test[f'CV{cv}_MLP_Confidence_Predictions'] = mlp_predictions[:, 1]\n            \n            qr_predictions = np.zeros((len(X_test), len(self.qr_parameters['quantiles'])))\n            for model in self.qr_models[f'CV{cv}']:\n                qr_predictions += model.predict(X_test[self.predictors]) / len(self.qr_models[f'CV{cv}'])\n\n            for i, quantile in enumerate(self.qr_parameters['quantiles']):    \n                X_test[f'CV{cv}_QR_{quantile}_Predictions'] = qr_predictions[:, i]\n                \n            hybrid_predictions = np.zeros((len(X_test), len(self.hybrid_parameters['quantiles'])))\n            for model in self.hybrid_models[f'CV{cv}']:\n                hybrid_predictions += model.predict(X_test[self.predictors]) / len(self.hybrid_models[f'CV{cv}'])\n\n            for i, quantile in enumerate(self.hybrid_parameters['quantiles']):    \n                X_test[f'CV{cv}_HYBRID_{quantile}_Predictions'] = hybrid_predictions[:, i]\n               \n    def plot_predictions(self, df, patient): \n              \n        mlp_prediction_columns = [f'CV{cv}_MLP_{target}_Predictions' for cv in self.cv for target in ['FVC', 'Confidence']]\n        if 'FVC' in df.columns:\n            mlp_scores = []\n            for cv in self.cv:            \n                score = self.laplace_log_likelihood_metric(df['FVC'], df[f'CV{cv}_MLP_FVC_Predictions'], df[f'CV{cv}_MLP_Confidence_Predictions'])\n                mlp_scores.append(round(score, 5))\n        \n        qr_prediction_columns = [f'CV{cv}_QR_{quantile}_Predictions' for cv in self.cv for quantile in self.qr_parameters['quantiles']]\n        if 'FVC' in df.columns:\n            qr_scores = []        \n            for i, cv in enumerate(self.cv):\n                score = self.laplace_log_likelihood_metric(df['FVC'], df[qr_prediction_columns[1 + (i * 3)]], (df[qr_prediction_columns[2 + (i * 3)]] - df[qr_prediction_columns[0 + (i * 3)]]))\n                qr_scores.append(round(score, 5))\n                \n        hybrid_prediction_columns = [f'CV{cv}_HYBRID_{quantile}_Predictions' for cv in self.cv for quantile in self.hybrid_parameters['quantiles']]\n        if 'FVC' in df.columns:\n            hybrid_scores = []        \n            for i, cv in enumerate(self.cv):\n                score = self.laplace_log_likelihood_metric(df['FVC'], df[hybrid_prediction_columns[1 + (i * 3)]], (df[hybrid_prediction_columns[2 + (i * 3)]] - df[hybrid_prediction_columns[0 + (i * 3)]]))\n                hybrid_scores.append(round(score, 5))\n\n        if 'FVC' in df.columns:\n            ax = df[(['Weeks', 'FVC'] + mlp_prediction_columns[::2] + qr_prediction_columns[1::3] + hybrid_prediction_columns[1::3])].set_index('Weeks').plot(figsize=(30, 6), style=['-b', 'r--', 'g--', 'b--', 'r:', 'g:', 'b:'])\n        else:\n            ax = df[(['Weeks'] + mlp_prediction_columns[::2] + qr_prediction_columns[1::3])].set_index('Weeks').plot(figsize=(30, 6), style=['r--', 'g--', 'b--', 'r:', 'g:', 'b:'])\n\n        ax.fill_between(df['Weeks'], df[qr_prediction_columns[2]], df[qr_prediction_columns[0]], alpha=0.05, label='CV1 QR Prediction Interval', color='red')\n        ax.fill_between(df['Weeks'], df[hybrid_prediction_columns[2]], df[hybrid_prediction_columns[0]], alpha=0.05, label='CV1 HYBRID Prediction Interval', color='blue')\n\n        ax.tick_params(axis='x', labelsize=20)\n        ax.tick_params(axis='y', labelsize=20)\n        ax.set_xlabel('')\n        ax.set_ylabel('')\n        if 'FVC' in df.columns:\n            ax.set_title(f'Patient: {patient} - MLP Scores: {mlp_scores} QR Scores: {qr_scores} HYBRID Scores: {hybrid_scores}', size=25, pad=25)\n        else:\n            ax.set_title(f'Patient: {patient}', size=25, pad=25)\n        ax.legend(prop={'size': 18})\n\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seed_everything(SEED)\n\nX_train = df_train.drop(columns=['FVC', 'Weeks'])\ny_train = df_train['FVC'].copy(deep=True)\n\n# Fill NaNs with mean\nfor feature in X_train.columns:\n    if X_train[feature].isnull().any():\n        X_train[feature] = X_train[feature].fillna(pd.concat([X_train[feature], df_test[feature]], axis=0).mean())\n        df_test[feature] = df_test[feature].fillna(pd.concat([X_train[feature], df_test[feature]], axis=0).mean())\n\n# Feature Scaling\nscaler = StandardScaler()\n# Scaling always hurts CV and LB score so decided not to scale any of the features\nscale_features = []\n\nif len(scale_features) > 0:\n    scaler.fit(pd.concat([X_train.loc[:, scale_features], df_test.loc[:, scale_features]], axis=0))\n    X_train.loc[:, scale_features] = scaler.transform(X_train.loc[:, scale_features])\n    df_test.loc[:, scale_features] = scaler.transform(df_test.loc[:, scale_features])\n\nmodel_parameters = {\n    'model': 'Stack',\n    'cv': [1],\n    'predictors': ['Age', 'Male', 'Female', 'Never smoked', 'Ex-smoker', 'Currently smokes', 'FVC_Baseline', 'Percent', 'Weeks_Passed',\n                   'smooth_lung_mean', 'smooth_lung_std', 'smooth_lung_skew', 'smooth_lung_kurt',\n                   'chest_diameter', 'smooth_lung_volume', 'smooth_lung_lower_level_calcification',\n                   'heart_diameter', 'heart_volume',\n                   'trachea_diameter', 'trachea_branching_angle'],\n    'mlp_parameters': {\n        'lr': 0.001,\n        'epochs': 25,\n        'batch_size': 2 ** 5\n    },\n    'qr_parameters': {\n        'quantiles': [0.25, 0.50, 0.75],\n        'lr': 0.001,\n        'epochs': 25,\n        'batch_size': 2 ** 5        \n    },\n    'hybrid_parameters': {\n    'quantiles': [0.25, 0.50, 0.75],\n    'lr': 0.001,\n    'epochs': 25,\n    'batch_size': 2 ** 5        \n    }\n}\n\nqr_mlp = QuantileRegressorMLP(**model_parameters)\nqr_mlp.train(X_train, y_train)\nqr_mlp.predict(df_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"for patient, df in list(df_train.groupby('Patient'))[:20]:        \n    qr_mlp.plot_predictions(df, patient)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for patient, df in list(df_test.groupby('Patient'))[:5]:        \n    qr_mlp.plot_predictions(df, patient)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Postprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"class SubmissionPipeline:\n    \n    def __init__(self, df_train, df_test):\n        \n        self.df_train = df_train\n        self.df_test = df_test\n        \n    def laplace_log_likelihood_metric(self, y_true, y_pred, sigma):\n        \n        sigma_clipped = np.maximum(sigma, 70)\n        delta_clipped = np.minimum(np.abs(y_true - y_pred), 1000)\n        score = - np.sqrt(2) * delta_clipped / sigma_clipped - np.log(np.sqrt(2) * sigma_clipped)\n\n        return np.mean(score)     \n        \n    def single_model(self, model):\n        \n        self.df_test['Patient_Week'] = self.df_test['Patient'].astype(str) + '_' + self.df_test['Weeks'].astype(str)\n        \n        prediction_cols = [col for col in self.df_train.columns if col.startswith(model)]\n        if model.split('_')[1] == 'MLP':\n            score = self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[prediction_cols[0]], self.df_train[prediction_cols[1]])\n            print(f'Single Model {model} Score: {score:.6}')\n            self.df_test['FVC'] = self.df_test[prediction_cols[0]]\n            self.df_test['Confidence'] = self.df_test[prediction_cols[1]]\n            \n        elif model.split('_')[1] == 'QR' or model.split('_')[1] == 'LGB':\n            score = self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[prediction_cols[1]], (self.df_train[prediction_cols[2]] - self.df_train[prediction_cols[0]]))\n            print(f'Single Model {model} Score: {score:.6}')\n            self.df_test['FVC'] = self.df_test[prediction_cols[1]]\n            self.df_test['Confidence'] = self.df_test[prediction_cols[2]] - self.df_test[prediction_cols[0]]\n            \n        print(f'\\n{self.df_test[[\"Patient_Week\", \"FVC\", \"Confidence\"]].describe()}')            \n        return self.df_test[['Patient_Week', 'FVC', 'Confidence']].copy(deep=True)\n        \n    def blend(self, by, model, cv):\n        \n        self.df_test['Patient_Week'] = self.df_test['Patient'].astype(str) + '_' + self.df_test['Weeks'].astype(str)\n        \n        if by == 'model':            \n            if model == 'MLP':\n                for df in [self.df_train, self.df_test]:\n                    df[f'{model}_FVC'] = (df['CV1_MLP_FVC_Predictions'] * 0.34) + (df['CV2_MLP_FVC_Predictions'] * 0.33) + (df['CV3_MLP_FVC_Predictions'] * 0.33)\n                    df[f'{model}_Confidence'] = (df['CV1_MLP_Confidence_Predictions'] * 0.34) + (df['CV2_MLP_Confidence_Predictions'] * 0.33) + (df['CV3_MLP_Confidence_Predictions'] * 0.33)\n                    \n                score = self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'{model}_FVC'], self.df_train[f'{model}_Confidence'])\n                single_model_scores = [round(self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'CV{i}_MLP_FVC_Predictions'], self.df_train[f'CV{i}_MLP_Confidence_Predictions']), 6) for i in range(1, 4)]\n                print(f'MLP Blend Score: {score:.6} - Single Model Scores: {single_model_scores}')\n                \n                self.df_test['FVC'] = self.df_test[f'{model}_FVC']\n                self.df_test['Confidence'] = self.df_test[f'{model}_Confidence']\n                \n            elif model == 'QR':\n                quantiles = [0.25, 0.50, 0.75]\n                for df in [self.df_train, self.df_test]:\n                    df[f'{model}_{quantiles[0]}_FVC'] = (df[f'CV1_QR_{quantiles[0]}_Predictions'] * 0.34) + (df[f'CV2_QR_{quantiles[0]}_Predictions'] * 0.33) + (df[f'CV3_QR_{quantiles[0]}_Predictions'] * 0.33)\n                    df[f'{model}_{quantiles[1]}_FVC'] = (df[f'CV1_QR_{quantiles[1]}_Predictions'] * 0.34) + (df[f'CV2_QR_{quantiles[1]}_Predictions'] * 0.33) + (df[f'CV3_QR_{quantiles[1]}_Predictions'] * 0.33)\n                    df[f'{model}_{quantiles[2]}_FVC'] = (df[f'CV1_QR_{quantiles[2]}_Predictions'] * 0.34) + (df[f'CV2_QR_{quantiles[2]}_Predictions'] * 0.33) + (df[f'CV3_QR_{quantiles[2]}_Predictions'] * 0.33)\n                    \n                score = self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'{model}_{quantiles[1]}_FVC'], (self.df_train[f'{model}_{quantiles[2]}_FVC'] - self.df_train[f'{model}_{quantiles[0]}_FVC']))\n                single_model_scores = [round(self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'CV{i}_QR_{quantiles[1]}_Predictions'], (self.df_train[f'CV{i}_QR_{quantiles[2]}_Predictions'] - self.df_train[f'CV{i}_QR_{quantiles[0]}_Predictions'])), 6) for i in range(1, 4)]\n                print(f'QR Blend Score: {score:.6} - Single Model Scores: {single_model_scores}')\n                \n                self.df_test['FVC'] = self.df_test[f'{model}_{quantiles[1]}_FVC']\n                self.df_test['Confidence'] = self.df_test[f'{model}_{quantiles[2]}_FVC'] - self.df_test[f'{model}_{quantiles[0]}_FVC']\n\n        elif by == 'cv':\n            quantiles = [0.25, 0.50, 0.75]\n            for df in [self.df_train, self.df_test]:\n                df[f'CV{cv}_FVC'] = (df[f'CV{cv}_MLP_FVC_Predictions'] * 0.5) + (df[f'CV{cv}_QR_{quantiles[1]}_Predictions'] * 0.5)\n                df[f'CV{cv}_Confidence'] = (df[f'CV{cv}_MLP_Confidence_Predictions'] * 0.5) + ((df[f'CV{cv}_QR_{quantiles[2]}_Predictions'] - df[f'CV{cv}_QR_{quantiles[0]}_Predictions']) * 0.5)\n                \n            score = self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'CV{cv}_FVC'], self.df_train[f'CV{cv}_Confidence'])\n            single_model_scores = [round(self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'CV{cv}_MLP_FVC_Predictions'], self.df_train[f'CV{cv}_MLP_Confidence_Predictions']), 6), round(self.laplace_log_likelihood_metric(self.df_train['FVC'], self.df_train[f'CV{cv}_QR_{quantiles[1]}_Predictions'], (self.df_train[f'CV{cv}_QR_{quantiles[2]}_Predictions'] - self.df_train[f'CV{cv}_QR_{quantiles[0]}_Predictions'])), 6)]\n            print(f'CV{cv} Blend Score: {score:.6} - Single Model Scores: {single_model_scores}')\n            \n            self.df_test['FVC'] = self.df_test[f'CV{cv}_FVC']\n            self.df_test['Confidence'] = self.df_test[f'CV{cv}_Confidence']\n            \n        print(f'\\n{self.df_test[[\"Patient_Week\", \"FVC\", \"Confidence\"]].describe()}') \n        return self.df_test[['Patient_Week', 'FVC', 'Confidence']].copy(deep=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = SubmissionPipeline(df_train, df_test)\ndf_submission = sub.blend(by='cv', model=None, cv=1)\ndf_submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}