{"cells":[{"metadata":{},"cell_type":"markdown","source":"For EDA part of this notebook, visit https://www.kaggle.com/gunawanmarbun/osic-pulmonary-eda/.\n\nFor inference notebook, visit https://www.kaggle.com/gunawanmarbun/osic-pulmonary-finetune. It use only middle slice + tabular data.\n\nThis notebook has 2 parts:\n* LGBM Modelling part, all ideas are from the osic-lgb-baseline notebook\n* CNN Modelling part, 4 model are developed in total. Pretrained weight of efficientnet does not seem to be helpful at all, I guess this is since we can have more that 255 unique values on our image data (depending on the window used of course)."},{"metadata":{},"cell_type":"markdown","source":"# Import Packages"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index\n!pip install ../input/osic-manual-bbox/efficientnet-1.1.0/efficientnet-1.1.0/ -f ./ --no-index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport gc\nimport glob\nimport math\nimport random\nfrom functools import partial, reduce\nfrom tqdm.auto import tqdm\nimport warnings\n\nfrom skimage.transform import resize\nimport pydicom\n\nimport scipy\nimport numpy as np\nimport pandas as pd\n\nimport lightgbm as lgb\nimport tensorflow as tf\nprint(f\"TF version: {tf.__version__}\")\nimport efficientnet.tfkeras as efn\n\nfrom sklearn.preprocessing import OneHotEncoder, OrdinalEncoder\nfrom sklearn.model_selection import GroupKFold, KFold, train_test_split\nfrom sklearn.metrics import mean_squared_error\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_slices_data(slices_data, n_cols=10, cmap='gray', **kwargs):\n    n_rows = math.ceil(slices_data.shape[0] / n_cols)\n    fig, axes = plt.subplots(n_rows, n_cols, figsize=(14, n_rows*1.5))\n    for img, ax in tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=slices_data.shape[0]):\n        ax.imshow(img, cmap=cmap, **kwargs)\n        ax.axis('off')\n    \n    missing_image_cnt = (n_rows * n_cols) - slices_data.shape[0]\n    if missing_image_cnt > 0:\n        for ax in axes.reshape(-1)[::-1][:-missing_image_cnt]:\n            ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utilities\n\n* `seed_everything`: to ensure reproducible results from model"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def seed_everything(seed=42):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n#     torch.manual_seed(seed)\n#     torch.cuda.manual_seed(seed)\n#     torch.backends.cudnn.deterministic = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Dataset Class Wrapper\n\nWrapper for convienient operation on dataset.\n* `OSICTrainDataset`: wrapper for training dataset\n* `OSICTestDataset`: wrapper for test dataset\n* `DICOMImages`: wrapper for DICOM images"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"class OSICTrainDataset:\n    def __init__(self, df):\n        self.df = df\n        self.create_base_df()\n        self._clean_dataset()\n        self._add_base_features()\n        self._add_col_id()\n        self.__sort_by_id()\n        \n    def create_base_df(self):\n        temp_dff = self.df.copy()\n        temp_dff['rank'] = temp_dff.groupby('Patient')['Weeks'].rank(method='min')\n        temp_dff = temp_dff[temp_dff['rank'] == 1]\n        temp_dff = temp_dff.drop_duplicates(subset='Patient')\n        self.base_df = temp_dff\n\n    def _clean_dataset(self):\n        \"\"\"\n        Preprocessing steps:\n            1. Drop duplicate Patient-Weeks combination\n        \"\"\"\n        self.__drop_duplicates()\n        \n    def __drop_duplicates(self):\n        before = self.df.shape[0]\n        self.df = self.df.drop_duplicates(subset=['Patient', 'Weeks'], keep='first').reset_index(drop=True)\n        after = self.df.shape[0]\n        print(f\"Dropped {before-after} rows of duplicate 'Patient-Weeks' values.\")\n        \n    def _add_base_features(self):\n        before = self.df.shape\n        temp_dff = self.df.copy()\n        temp_dff['rank'] = temp_dff.groupby('Patient')['Weeks'].rank(method='min')\n        \n        # Treate every FVC measurement as a baseline\n        all_dfs = []\n        for rank in sorted(temp_dff['rank'].unique()):\n            all_dfs.append(self.__get_ranked_base_features(temp_dff, rank))\n        self.df = pd.concat(all_dfs)\n        self.df = self.df.reset_index(drop=True)\n        after = self.df.shape\n        print(f\"Before-After shape adding base features: {before} {after}\")\n        \n    def __get_ranked_base_features(self, temp_dff, rank):\n        temp_df = temp_dff[temp_dff['rank'] == rank].reset_index(drop=True)\n        temp_df = temp_df.drop(['Sex', 'SmokingStatus', 'rank'], axis=1)\n        temp_df = temp_df.rename(columns={\n            \"FVC\": \"FVC_base\",\n            \"Percent\": \"Percent_base\",\n            \"Age\": \"Age_base\",\n            \"Weeks\": \"Weeks_base\"\n        })\n        temp_df = self.df[['Patient', 'Weeks', 'FVC', 'Sex', 'SmokingStatus']].merge(\n            temp_df,\n            how='inner',\n            on='Patient',\n        )\n        temp_df['Weeks_passed'] = temp_df['Weeks'] - temp_df['Weeks_base']\n        temp_df = temp_df[temp_df['Weeks_passed'] != 0]  # drop base observation\n        return temp_df\n    \n    def _add_col_id(self):\n        col_id = 'Patient_Week'\n        self.df[col_id] = self.df['Patient'] + '_' + self.df['Weeks'].astype(str)\n        print(f\"ID column '{col_id}' added. After adding shape: {self.df.shape}\")\n        \n    def __sort_by_id(self):\n        self.df = self.df.sort_values(by='Patient').reset_index(drop=True)\n        \n        \nclass OSICTestDataset:\n    def __init__(self, test_df, submission_df):\n        self.df = test_df\n        self.submission_df = submission_df\n        self._prepare_test_df()\n        self.__sort_by_id()\n    \n    def _prepare_test_df(self):\n        before = self.df.shape\n        self.submission_df[['Patient', 'Weeks']] = self.submission_df['Patient_Week'].str.split('_', expand=True)\n        self.df = self.submission_df.drop(['FVC', 'Confidence'], axis=1).merge(\n            self.df.rename(columns={\n                \"FVC\": \"FVC_base\",\n                \"Percent\": \"Percent_base\",\n                \"Age\": \"Age_base\",\n                \"Weeks\": \"Weeks_base\"\n            }),\n            how='left',\n            on='Patient'\n            )\n        self.df['Weeks'] = self.df['Weeks'].astype(int)\n        self.df['Weeks_passed'] = self.df['Weeks'] - self.df['Weeks_base']\n        self.df = self.df.reset_index(drop=True)\n        after = self.df.shape\n        print(f\"Before-After shape adding base features: {before} {after}\")\n        \n    def __sort_by_id(self):\n        self.df = self.df.sort_values(by='Patient').reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class DICOMImages:\n    DOUBLE_IDS = ['ID00078637202199415319443']\n    \"\"\"Wrapper for multiple slices of a patient CT-Scan results.\"\"\"\n    def __init__(self, id, dirpath='../input/osic-pulmonary-fibrosis-progression/train/'):\n        self.id = id\n        self.basepath = os.path.join(dirpath, self.id)\n        self.filepaths = glob.glob(os.path.join(self.basepath, \"*.dcm\"))\n        if self.id in self.DOUBLE_IDS:\n            self.filepaths = self.filepaths[:len(self.filepaths)//2]\n        sort_nicely(self.filepaths)\n        \n    def __iter__(self):\n        for filepath in self.filepaths:\n            yield pydicom.dcmread(filepath)\n\n    def __len__(self):\n        return len(self.filepaths)\n    \n    @property\n    def image_type(self):\n        \"\"\"\n        Infer dicom image type by its first slice metadata.\n        Categories:\n            - 'zero' : Rescale Intercept value is 0\n            - 'not-zero': Rescale Intercept value is either -1000 or -1024\n        \"\"\"\n        mapper = {0: 'zero'}\n        rescale_intercept = self.get_dicom_metadata(self.get_slice(index=0))['Rescale Intercept']\n        return {\n            'name': mapper.get(rescale_intercept, 'not-zero'),\n            'rescale_intercept': rescale_intercept\n        }\n        \n    @property\n    def slices(self):\n        return list(self.__iter__())\n    \n    def get_slice(self, index):\n        return pydicom.dcmread(self.filepaths[index])\n    \n    @property\n    def df(self):\n        return pd.DataFrame(\n            [self.get_dicom_metadata(slice) for slice in self.__iter__()]\n        )\n    \n    @staticmethod\n    def get_dicom_metadata(slice):\n        dict_res = {}\n        for x in slice.values():\n            if isinstance(x, pydicom.dataelem.RawDataElement):\n                metadata = pydicom.dataelem.DataElement_from_raw(x)\n            else:\n                metadata = x\n            if metadata.name == 'Pixel Data':\n                continue\n            dict_res.update({\n                f\"{metadata.name}\": metadata.value\n            })\n        return dict_res\n    \n    @property\n    def slices_data(self):\n        return np.stack([self._to_HU(slice) for slice in self.__iter__()])\n    \n    @property\n    def middle_filepath(self):\n        return self.filepaths[(len(self.filepaths)-1) // 2]\n\n    @property\n    def middle_slice_data(self):\n        mid_slice_index = (len(self.filepaths)-1) // 2\n        return self._to_HU(pydicom.dcmread(self.filepaths[mid_slice_index]))\n        \n    def sampled_slices_data(self, n_samples=30, ret_paths=False):\n        if len(self.filepaths) < n_samples:\n            msg = f\"Total slices is less than number of samples: {len(self.filepaths)} < {n_samples}.\"\n            msg += \" Number of samples default to total slices.\"\n            warnings.warn(msg, UserWarning)\n            n_samples = len(self.filepaths)\n        sample_indexes = np.linspace(0, len(self.slices)-1, n_samples).astype(int)\n        sampled_slices = np.array(self.slices)[sample_indexes]\n        if ret_paths:\n            sample_filepaths = np.array(self.filepaths)[sample_indexes]\n            return np.stack([self._to_HU(slice) for slice in sampled_slices]), sample_filepaths\n        else:\n            return np.stack([self._to_HU(slice) for slice in sampled_slices])\n\n    @staticmethod\n    def _to_HU(slice):\n        intercept, slope = slice.RescaleIntercept, slice.RescaleSlope\n        \n        slice_data = slice.pixel_array.astype(np.int16)\n        slice_data[slice_data <= -1000] = 0\n        \n        if slope != 1:\n            slice_data = slope * slice_data.astype(np.float64)\n            slice_data = slice_data.astype(np.int16)\n            \n        slice_data += np.int16(intercept)\n        return slice_data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utilities\n\n* `sort_nicely`: sort filepaths numerically"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import re\n\ndef tryint(s):\n    try:\n        return int(s)\n    except ValueError:\n        return s\n    \ndef alphanum_key(s):\n    \"\"\" Turn a string into a list of string and number chunks.\n        \"z23a\" -> [\"z\", 23, \"a\"]\n    \"\"\"\n    return [ tryint(c) for c in re.split('([0-9]+)', s) ]\n\ndef sort_nicely(l):\n    \"\"\" Sort the given list in the way that humans expect.\n    \"\"\"\n    l.sort(key=alphanum_key)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# References\n\n## LightGBM Baseline by yasufuminakama\nBasic lgbm for patient's tabular features. https://www.kaggle.com/yasufuminakama/osic-lgb-baseline"},{"metadata":{},"cell_type":"markdown","source":"# Data Preprocessing\n\nPreprocess both tabular and image data before modelling."},{"metadata":{"trusted":true},"cell_type":"code","source":"basepath = \"../input/osic-pulmonary-fibrosis-progression/\"\ntrain_df = pd.read_csv(f\"{basepath}train.csv\")\ntest_df = pd.read_csv(f\"{basepath}test.csv\")\nsubmission_df = pd.read_csv(f\"{basepath}sample_submission.csv\")\nprint(train_df.shape, test_df.shape, submission_df.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Patient's Data Preprocessing\n\nTrain Dataset:\n* Remove duplicate `Patient` and `Week` rows (total: 7)\n* Treate every FVC measurement as baseline. Unpack base column feature for each baseline.\n* Generate ID column `Patient_Week`\n\nTest Dataset:\n* Combine `submission_df` and `test_df`, treating `test_df` FVC measurement as baseline."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dataset = OSICTrainDataset(train_df)\ntest_dataset = OSICTestDataset(test_df, submission_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# DICOM Images Preprocessing\n\nFilter out bad ID that can't be read with `pydicom`. It might be fixed by installing `gdcm` but it's quite complicated now since we don't have access to internet. Also filter out those failed image ids."},{"metadata":{"trusted":true},"cell_type":"code","source":"bbox_map = pd.read_csv('../input/osic-manual-bbox/threshold_all.csv')\nprint(bbox_map.shape)\ncode_filter = (bbox_map['x'] == 0) & (bbox_map['y'] == 0) & (bbox_map['width'] == 1) & (bbox_map['height'] == 1)\nbbox_map = bbox_map[~code_filter]\nprint(bbox_map.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Ensure our test patient is in our filtered dataset\nfor x in test_df.Patient:\n    print(x in bbox_map.patient.tolist())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def check_bad_images(patient_ids):\n    bad_ids = []\n    exceptions = []\n    pbar = tqdm(range(100), leave=False)  # dummy range\n    for patient_id in tqdm(patient_ids, leave=False):\n        DICOMImage = DICOMImages(patient_id)\n        pbar.reset(total=len(DICOMImage))\n        for dicom_image in DICOMImage:\n            try:\n                _ = dicom_image.pixel_array\n            except Exception as e:\n                bad_ids.append(patient_id)\n                exceptions.append(e)\n                break\n            finally:\n                pbar.update()\n        pbar.refresh()\n    return bad_ids, exceptions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"bad_ids, exceptions = check_bad_images(train_df.Patient.unique())\nfor bad_id, exception in zip(bad_ids, exceptions):\n    print(bad_id, exception)\n# bad_ids = ['ID00011637202177653955184', 'ID00052637202186188008618']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create image data generator for modelling."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def load_image(filename, label=None):\n    image = preprocess_dicom(filename)\n    if label is None:\n        return image\n    else:\n        return image, label\n    \ndef preprocess_dicom(patient_id, img_width=512, img_height=512):\n    middle_slice_data = DICOMImages(patient_id).middle_slice_data\n    middle_slice_data = rescale(middle_slice_data)\n    middle_slice_data = np.expand_dims(middle_slice_data, axis=-1)\n    resized_slice_data = tf.image.resize_with_crop_or_pad(middle_slice_data, img_width, img_height)\n    return resized_slice_data\n        \ndef rescale(slice_data):\n    min_, max_ = slice_data.min(), slice_data.max()\n    rescaled = (slice_data-min_) / (max_-min_)\n    total_pixel_count = reduce(lambda a, b: a*b, slice_data.shape)\n    assert (rescaled >= 0).sum() == total_pixel_count\n    assert (rescaled <= 1).sum() == total_pixel_count\n    return rescaled","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Segmentation preprocessing is here."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from skimage.filters import threshold_otsu, median\nfrom skimage.segmentation import clear_border\nfrom skimage import morphology\nfrom scipy.ndimage import binary_fill_holes\n\n\ndef lung_segment(img):\n    thresh = threshold_otsu(img)\n    binary = img <= thresh\n\n    lungs = median(clear_border(binary))\n    lungs = morphology.binary_closing(lungs, selem=morphology.disk(7))\n    lungs = binary_fill_holes(lungs)\n\n    final = lungs*img\n    final[final == 0] = np.min(img)\n\n    return final, lungs\n\ndef morphological_segmentation(img):\n    segmented_img, _ = lung_segment(img)\n    return segmented_img\n\ndef segment_lung(slice_data, image_type, segment_func):\n    if image_type == 'zero':\n        slice_data[slice_data == 0] = -1000\n    segmented_image = segment_func(threshold_slices_data(slice_data, low=-1000, high=-400))\n    return segmented_image\n\n\nclass BoundingBox:\n    \"\"\"Initiation of bbox follows matplotlib Rectangle patch\"\"\"\n    def __init__(self, xy, width, height):\n        self.x, self.y = xy\n        self.width = width\n        self.height = height\n        \n    @property\n    def attribute_list(self):\n        return [(self.x, self.y), self.width, self.height]\n    \n    def __repr__(self):\n        return f\"Bbox (bottom left width height): {self.x} {self.y} {self.width} {self.height}\"\n\n    \ndef crop_recenter(image, bbox, pad_value=-1000):\n    x, y, width, height = bbox.x, bbox.y, bbox.width, bbox.height\n    cropped_image = image[ y:y+height, x:x+width ]\n    out_height, out_width = image.shape\n    \n    padded_image = np.ones(image.shape, dtype=np.int16) * pad_value\n    x_start = (out_width - width) // 2\n    y_start = (out_height - height) // 2\n    padded_image[ y_start:y_start+height, x_start:x_start+width ] = cropped_image\n    return padded_image\n\n\ndef threshold_slices_data(slices_data, low=-1000, high=-400):\n    copy = slices_data.copy()\n    copy[copy < low] = low\n    copy[copy > high] = high\n    return copy","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Generate the data generator, we will use keras `Sequence` class."},{"metadata":{"trusted":true},"cell_type":"code","source":"import collections\n\n\ndef create_scaler(min_, max_):\n    def scalar_scaler(val):\n        return (val - min_) / (max_ - min_)\n    return scalar_scaler\n\ndef osic_cat_encoder(cat):\n    mapper = {\n        'Male': 0,\n        'Female': 1,\n        'Never smoked': [0, 0],\n        'Ex-smoker': [0, 1],\n        'Currently smokes': [1, 0],\n    }\n    return mapper.get(cat, [1, 1])\n\ndef flatten(l):\n    for el in l:\n        if isinstance(el, collections.Iterable) and not isinstance(el, (str, bytes)):\n            yield from flatten(el)\n        else:\n            yield el","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (Chosen) CNN-MLP Dataset Generator\n\nThis is on `train_dataset.df` dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"base_df = train_dataset.df.copy().reset_index(drop=True)\nbase_test_df = test_dataset.df.copy().reset_index(drop=True)\nprint(base_df.shape, base_test_df.shape)\nbase_df = base_df[base_df.Patient.isin(bbox_map.patient)]\nbase_test_df = base_test_df[base_test_df.Patient.isin(bbox_map.patient)]\nprint(base_df.shape, base_test_df.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_pred_num = ['FVC_base', 'Percent_base', 'Age_base', 'Weeks_passed']\ncols_pred_cat = ['Sex', 'SmokingStatus']\ncols_pred = cols_pred_num + cols_pred_cat\n\n# Apply scaler and cat encoder, look for min-max range\n# in both train and test dataset\nfor col in cols_pred_num:\n    min_ = min(base_df[col].min(), test_dataset.df[col].min())\n    max_ = max(base_df[col].max(), test_dataset.df[col].max())\n    scaler = create_scaler(min_, max_)\n    base_df[col] = base_df[col].apply(scaler)\n    base_test_df[col] = base_test_df[col].apply(scaler)\nfor col in cols_pred_cat:\n    base_df[col] = base_df[col].apply(osic_cat_encoder)\n    base_test_df[col] = base_test_df[col].apply(osic_cat_encoder)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_X_y_cnnmlp(df):\n    base_df = df.copy()\n\n    mapper = {}\n    for patient_id in tqdm(base_df.Patient.unique(), leave=False):\n        dicom = DICOMImages(patient_id)\n        mapper.update({\n            f'{patient_id}': {\n                'filepath': dicom.middle_filepath,\n                'image_type': dicom.image_type['name'],\n            }\n        })\n\n    base_df['id'] = base_df['Patient']\n    base_df['filepath'] = base_df['Patient'].apply(lambda id: mapper[id]['filepath'])\n    base_df['image_type'] = base_df['Patient'].apply(lambda id: mapper[id]['image_type'])\n    \n    X = base_df[['id', 'filepath', 'image_type']].to_dict(orient='records')\n    for x, vector in tqdm(zip(X, base_df[cols_pred].values), leave=False, total=len(X)):\n        x.update({\n            'vector': list(flatten(vector.tolist()))\n        })\n    X = np.array(X)\n    if 'FVC' in base_df.columns.values:\n        y = base_df['FVC'].values\n        return X, y\n    else:\n        return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Take only middle slice data\nmapper = {}\nfor patient_id in tqdm(base_df.Patient.unique(), leave=False):\n    dicom = DICOMImages(patient_id)\n    mapper.update({\n        f'{patient_id}': {\n            'filepath': dicom.middle_filepath,\n            'image_type': dicom.image_type['name'],\n        }\n    })","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base_df['id'] = base_df['Patient']\nbase_df['filepath'] = base_df['Patient'].apply(lambda id: mapper[id]['filepath'])\nbase_df['image_type'] = base_df['Patient'].apply(lambda id: mapper[id]['image_type'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = base_df[['id', 'filepath', 'image_type']].to_dict(orient='records')\nfor x, vector in tqdm(zip(X, base_df[cols_pred].values), leave=False, total=len(X)):\n    x.update({\n        'vector': list(flatten(vector.tolist()))\n    })\nX = np.array(X)\ny = base_df['FVC'].values\nprint(X.shape, y.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## CNN+ Dataset Generator\n\nOptionally introduce noise to FVC measurement. Pertubate the FVC according to the error rate of measurement. Based from this [source](http://cpsa.ca/wp-content/uploads/2015/03/Spirometry_Flow_Volume_Measurement_Guidelines.pdf) and [ATS standard](https://www.atsjournals.org/doi/full/10.1164/rccm.201908-1590ST#_i12) the acceptable measurement error is around +-3%"},{"metadata":{"trusted":true},"cell_type":"code","source":"# weeks_scaler = create_scaler(train_df['Weeks'].min(), train_df['Weeks'].max())\n# percent_scaler = create_scaler(train_df['Percent'].min(), train_df['Percent'].max())\n# age_scaler = create_scaler(train_df['Age'].min(), train_df['Age'].max())\n# # fvc_scaler = create_scaler(train_df['FVC'].min(), train_df['FVC'].max())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# base_df = train_dataset.base_df.reset_index(drop=True)\n# base_df = base_df[base_df.Patient.isin(bbox_map.patient)]\n# base_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# base_df['Weeks'] = base_df['Weeks'].apply(weeks_scaler)\n# base_df['Percent'] = base_df['Percent'].apply(percent_scaler)\n# base_df['Age'] = base_df['Age'].apply(age_scaler)\n# base_df['Sex'] = base_df['Sex'].apply(osic_cat_encoder)\n# base_df['SmokingStatus'] = base_df['SmokingStatus'].apply(osic_cat_encoder)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# all_dicoms = [DICOMImages(patient_id) for patient_id in base_df.Patient.values]\n# X = []\n# y = []\n# for dicom in tqdm(all_dicoms, leave=False):\n#     patient_id = dicom.id\n#     image_type = dicom.image_type['name']\n#     base_data = base_df[base_df['Patient'] == patient_id].values\n#     base_FVC = base_data[0, 2]\n#     base_PCT = base_data[0, 3]\n    \n#     sampled_slices_data, sampled_filepaths = dicom.sampled_slices_data(100, ret_paths=True)  # sample or all\n#     ratio_from_middle = 0.25 * len(sampled_filepaths)\n#     mid_slice_index = len(sampled_filepaths) / 2\n#     left = int(mid_slice_index - ratio_from_middle)\n#     right = int(mid_slice_index + ratio_from_middle)\n#     sampled_filepaths = sampled_filepaths[left:right]\n#     total = len(sampled_filepaths)\n#     base_data = np.tile(base_data, (total, 1))\n\n#     # Assume 99.7% -- 3 standard deviation. For PCT, since it\n#     # is ratio of FVC then the error rate is sqrt(2) of fvc error\n#     # assuming error rate of both measurement is 3%\n#     FVC_noise = np.random.normal(0, 0.03/3, total) * base_FVC\n#     PCT_noise = np.random.normal(0, math.sqrt(2)*0.03/3, total) * base_PCT\n# #     FVC_noise = 0\n# #     PCT_noise = 0\n#     base_data[:, 2] = base_data[:, 2] + FVC_noise\n#     base_data[:, 3] = base_data[:, 3] + PCT_noise\n    \n#     # Change ID to filepaths\n#     base_data[:, 0] = sampled_filepaths\n    \n#     y_index = 2\n#     id_index = 0\n#     for data in base_data:\n#         x_index = list(set(range(7)) - set([y_index, id_index]))\n#         x = {\n#             'id': patient_id,\n#             'filepath': data[0],\n#             'image_type': image_type,\n#             'vector': list(flatten(data[x_index].tolist())),\n#         }\n#         X.append(x)\n#         y.append(data[y_index])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create Dataset Generator"},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\nfrom tensorflow.keras.utils import Sequence\n\n\nclass OSIC_CNNMLP_ImageGenerator(Sequence):\n    def __init__(self, vectors, labels, bbox_map, output_shape, batch_size=4, num_batch=0,\n                 segmented=True, shuffle=True, debug=False):\n        self.vectors = vectors\n        self.labels = labels\n        self.bbox_map = bbox_map\n        self.output_shape = output_shape\n        self.batch_size = batch_size\n        self.num_batch = num_batch\n        self.segmented = segmented\n        self.shuffle = shuffle\n        self.debug = debug\n        self._n = len(self.vectors)\n        self.on_epoch_end()\n\n    def __len__(self):\n        ct = len(self.vectors) // self.batch_size\n        ct += int((len(self.vectors) % self.batch_size)!=0)\n        self.num_batch = ct\n        return ct\n    \n    def __getitem__(self, batch_index):\n        indexes_in_batch = self.indexes[\n            batch_index * self.batch_size:(batch_index + 1) * self.batch_size\n        ]\n\n        selected_vectors = self.vectors[indexes_in_batch]\n        selected_labels = self.labels[indexes_in_batch]\n        out_height, out_width = self.output_shape\n        X = np.stack([self.__preprocess_vector(vector, out_width, out_height) for vector in selected_vectors])\n        X_meta = np.stack([vector['vector'] for vector in selected_vectors]).astype(np.float32)\n        y = np.expand_dims(selected_labels.astype(np.float32), axis=-1)\n        if self.debug:\n            selected_ids = [v['id'] for v in selected_vectors]\n            return [X, X_meta], y, selected_ids\n        else:\n            return [X, X_meta], y\n    \n    def on_epoch_end(self):\n        self.indexes = np.arange(self._n)\n        if self.shuffle == True:\n            np.random.shuffle(self.indexes)\n\n    def __preprocess_vector(self, vector, img_width=512, img_height=512):\n        slice_data = self.__to_HU(pydicom.dcmread(vector['filepath']))\n        image_type = vector['image_type']\n        patient_id = vector['id']\n        \n        if self.segmented:\n            try:\n                slice_data = segment_lung(slice_data, image_type, morphological_segmentation)\n            except:\n                print(\"Segmentation failed, returning original image.\")\n        x, y, width, height = self.bbox_map.loc[bbox_map.patient == patient_id, ['x', 'y', 'width', 'height']].values[0]\n        bbox = BoundingBox((x, y), width, height)\n        \n        slice_data = crop_recenter(slice_data, bbox)\n        slice_data = self.__rescale(slice_data)\n        resized_slice_data = resize(slice_data, (img_height, img_width), anti_aliasing=True)\n        resized_slice_data = np.expand_dims(resized_slice_data, axis=-1)\n        resized_slice_data = resized_slice_data.astype(np.float32)\n        return resized_slice_data\n    \n    def __rescale(self, slice_data):\n        min_, max_ = slice_data.min(), slice_data.max()\n        rescaled = (slice_data-min_) / (max_-min_)\n        \n        if np.isfinite(rescaled).all():\n            # total_pixel_count = reduce(lambda a, b: a*b, slice_data.shape)\n            # assert (rescaled >= 0).sum() == total_pixel_count\n            # assert (rescaled <= 1).sum() == total_pixel_count\n            return rescaled\n        else:\n            print(\"Rescaling failed, returning np.zeros() with original shape.\")\n            return np.zeros(slice_data.shape)\n\n    \n    def __to_HU(self, slice):\n        intercept, slope = slice.RescaleIntercept, slice.RescaleSlope\n        \n        slice_data = slice.pixel_array.astype(np.int16)\n        slice_data[slice_data <= -1000] = 0\n        \n        if slope != 1:\n            slice_data = slope * slice_data.astype(np.float64)\n            slice_data = slice_data.astype(np.int16)\n            \n        slice_data += np.int16(intercept)\n        return slice_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, random_state=42)\nprint(len(X_train), len(X_valid), len(y_train), len(y_valid))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"output_shape = (456, 456)\ntrain_image_dataset = OSIC_CNNMLP_ImageGenerator(\n    np.array(X_train), np.array(y_train),\n    bbox_map,\n    output_shape=output_shape, batch_size=8, segmented=True, shuffle=True, debug=True\n)\nvalid_image_dataset = OSIC_CNNMLP_ImageGenerator(\n    np.array(X_valid), np.array(y_valid),\n    bbox_map,\n    output_shape=output_shape, batch_size=8, segmented=True, shuffle=False\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Debug our dataset generator by seeing the output."},{"metadata":{"trusted":true},"cell_type":"code","source":"for data, labels, ids in tqdm(train_image_dataset, leave=False, total=len(X_train)):\n    try:\n        images, vector = data\n        print(vector.shape)\n        plot_slices_data(images.reshape(images.shape[0], 456, 456), cmap='Blues_r', n_cols=8)\n    except:\n        print(ids)\n    break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Modelling\n\nUsing the knowledge we have built before, now it's time to do the exciting part! Seed everything for reproducible results."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"class ModelExtractionCallback(object):\n    \"\"\"Callback class for retrieving trained model from lightgbm.cv()\n    NOTE: This class depends on '_CVBooster' which is hidden class,\n    so it might doesn't work if the specification is changed.\n    \"\"\"\n\n    def __init__(self):\n        self._model = None\n\n    def __call__(self, env):\n        # Saving _CVBooster object.\n        self._model = env.model\n\n    def _assert_called_cb(self):\n        if self._model is None:\n            # Throw exception if the callback class is not called.\n            raise RuntimeError('callback has not called yet')\n\n    @property\n    def boosters_proxy(self):\n        self._assert_called_cb()\n        # return Booster object\n        return self._model\n\n    @property\n    def raw_boosters(self):\n        self._assert_called_cb()\n        # return list of Booster\n        return self._model.boosters\n\n    @property\n    def best_iteration(self):\n        self._assert_called_cb()\n        # return boosting round when early stopping.\n        return self._model.best_iteration\n\ndef get_proxy_boosters_best_iter(extraction_cb):\n    proxy = extraction_cb.boosters_proxy\n    boosters = extraction_cb.raw_boosters\n    best_iteration = extraction_cb.best_iteration\n    return proxy, boosters, best_iteration\n\ndef loss_func(y_true, y_pred, weight):\n    confidence = weight\n    sigma_clipped = max(confidence, 70)\n    diff = abs(y_true - y_pred)\n    delta = min(diff, 1000)\n    score = -math.sqrt(2)*delta/sigma_clipped - np.log(math.sqrt(2)*sigma_clipped)\n    return -score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"GLOBAL CONFIGS\"\"\"\nSEED = 42\nseed_everything(SEED)\n\ncols_num = ['FVC_base', 'Percent_base', 'Age_base', 'Weeks_passed']\ncols_cat = ['Sex', 'SmokingStatus']\ncols_cat_oe = [c+'_oe' for c in cols_cat]\ncols_pred = cols_num + cols_cat_oe\n\ncol_target = 'FVC'\ncol_target_2 = 'Confidence'\ncol_score = 'Score'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model 1 (LGBM, Tabular)\n\n* Predict `FVC_pred`\n* Optimize `Confidence`\n* Predict `Confidence`, adding `FVC` and `FVC_pred` as features\n* Use all folds booster prediction by averaging"},{"metadata":{"trusted":true},"cell_type":"code","source":"lgbm_param = {\n    'objective': 'regression',\n    'boosting': 'gbdt',\n    'metric': 'rmse',\n    'learning_rate': 0.01,\n    'num_leaves': 31,\n    'max_depth': 2,\n    'colsample_bytree': 0.8,\n    'subsample': 0.8,\n    'subsample_freq': 1,\n    'num_threads': os.cpu_count() - 1,\n    'seed': 42\n}\n\nlgbm_train_param = {\n    \"verbose_eval\": 100,\n    \"num_boost_round\": 100000,\n    \"early_stopping_rounds\": 100,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_oe = OrdinalEncoder()\ntrain_dataset.df[cols_cat_oe] = train_oe.fit_transform(train_dataset.df[cols_cat])\ntest_dataset.df[cols_cat_oe] = train_oe.fit_transform(test_dataset.df[cols_cat])\n\ntrain_set = lgb.Dataset(train_dataset.df[cols_pred],\n                        label=train_dataset.df[col_target],\n                        group=train_dataset.df['Patient'].value_counts().sort_index())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"extraction_cb = ModelExtractionCallback()\n\nbst = lgb.cv(\n    lgbm_param,\n    train_set,\n    **lgbm_train_param,\n    folds = GroupKFold(n_splits=5),\n    seed=SEED,\n    callbacks=[extraction_cb]\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"proxy, boosters, best_iteration = get_proxy_boosters_best_iter(extraction_cb)\npredictions = proxy.predict(train_dataset.df[cols_pred], num_iteration=best_iteration)\nfor i, preds in enumerate(predictions):\n    rmse = np.sqrt(mean_squared_error(train_dataset.df[col_target], preds))\n    print(f\"Fold {i} rmse: {rmse}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"booster_id = 0\nfig, ax = plt.subplots(1, 2, figsize=(14, 4))\nlgb.plot_importance(boosters[booster_id], importance_type='gain', ax=ax[0])\nlgb.plot_importance(boosters[booster_id], importance_type='split', ax=ax[1])\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score_dataset(df):\n    scores = []\n    rows = df[['FVC', 'FVC_pred', 'Confidence']].values\n    for y_true, y_pred, weight in rows:\n        score = loss_func(y_true, y_pred, weight)\n        scores.append(score)\n    df[col_score] = scores\n    return -np.mean(scores)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Predict FVC\ntrain_dataset.df['FVC_pred'] = np.array(predictions).mean(axis=0)\ntest_dataset.df['FVC_pred'] = np.array(proxy.predict(test_dataset.df[cols_pred], num_iteration=best_iteration)).mean(axis=0)\n\n# Optimize score\ntrain_dataset.df['Confidence'] = 100\nnon_optimized_score = score_dataset(train_dataset.df)\n\nresults = []\nweight_arr = [100]\nfor y_true, y_pred in tqdm(train_dataset.df[['FVC', 'FVC_pred']].values, leave=False):\n    loss_partial = partial(loss_func, y_true, y_pred)\n    result = scipy.optimize.minimize(loss_partial, weight_arr, method='SLSQP')\n    x = result['x']\n    results.append(x[0])\n\ntrain_dataset.df['Confidence'] = results\noptimized_score = score_dataset(train_dataset.df)\nprint(f\"Non optimized score: {non_optimized_score}\")\nprint(f\"Optimized score: {optimized_score}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_pred_2 = list(set(cols_pred) - set(['FVC_base']))\ntrain_set = lgb.Dataset(train_dataset.df[cols_pred_2],\n                        label=train_dataset.df[col_target_2],\n                        group=train_dataset.df['Patient'].value_counts().sort_index())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"extraction_cb = ModelExtractionCallback()\n\nbst = lgb.cv(\n    lgbm_param,\n    train_set,\n    **lgbm_train_param,\n    folds = GroupKFold(n_splits=5),\n    seed=SEED,\n    callbacks=[extraction_cb]\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"proxy, boosters, best_iteration = get_proxy_boosters_best_iter(extraction_cb)\npredictions = proxy.predict(train_dataset.df[cols_pred_2], num_iteration=best_iteration)\nfor i, preds in enumerate(predictions):\n    rmse = np.sqrt(mean_squared_error(train_dataset.df[col_target_2], preds))\n    print(f\"Fold {i} rmse: {rmse}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Predict Confidence\ntrain_dataset.df[col_target_2] = np.array(predictions).mean(axis=0)\ntest_dataset.df[col_target_2] = np.array(proxy.predict(test_dataset.df[cols_pred_2], num_iteration=best_iteration)).mean(axis=0)\nprint(f\"Train score: {score_dataset(train_dataset.df)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model 2 (CNN, Ratio of Slices + Tabular Data)\n\nI ask myself on wether a radiologist scan through all slices or just sample some **big enough segmented lung** slices and start observing. Also I worry is that a model that is meticulously trained on all slices will more likely to learn the noise pattern instead of the real signal, I might be wrong, though.\n* Take ratio of all slices (or middle slice only) and train `FVC` measurement on that slice\n* Image input size is (height, width, 1) being a grayscaled image. If using pretrained model, will use that model native input shape. For efficientnet the image is simply concated together.\n* Image was simply resized to the desired input image\n* Qantile regression, hence multiple outputs from the model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow.keras.backend as K\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.layers import Dense, Input, Flatten, Conv2D, MaxPooling2D, GlobalAveragePooling2D, Concatenate\nfrom tensorflow.keras.optimizers import Adam, RMSprop, SGD\nfrom tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, LearningRateScheduler\n\n\nquantiles = [0.2, 0.5, 0.8]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score(y_true, y_pred):\n    y_true = tf.dtypes.cast(y_true, tf.float32)\n    y_pred = tf.dtypes.cast(y_pred, tf.float32)\n    C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return -K.mean(metric)\n\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss\n\ndef qloss(y_true, y_pred):\n    q = tf.constant(quantiles, dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return K.mean(v)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Model Zoo\n* (chosen) `effnet_mlp_model`, middle slice data only + 7 length vector of tabular data\n* `efficientnet_plus_model`, middle slice data only + 6 length vector of tabular data, no MLP for tabular data\n* `efficientnet_model`, middle slice data only, around -8.2x score\n* `dense_plus_model1`, middle slice data only + 6 length vector of tabular data, no MLP for tabular data\n* `dense_model1`, middle slice data only, around -8.0x score"},{"metadata":{"trusted":true},"cell_type":"code","source":"def effnet_mlp_model(input_shape, output_shape):\n    inp = Input(shape=input_shape)\n    base_model = efn.EfficientNetB5(\n        input_shape=input_shape,\n        weights=None,\n        include_top=False,\n    )\n    x1 = base_model(inp)\n    x1 = GlobalAveragePooling2D()(x1)\n    \n    inp2 = Input(shape=(7,))\n    x2 = Dense(100, activation='relu')(inp2)\n    x2 = Dense(100, activation='relu')(x2)\n    \n    x = Concatenate()([x1, x2])\n    output = Dense(output_shape, activation='linear', name='output')(x)\n    \n    model = Model(inputs=[inp, inp2], outputs=output, name='cnn_mlp_only_mid')\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def efficientnet_plus_model(input_shape, output_shape):\n    inp = Input(shape=input_shape)\n#     inp_concat = Concatenate()([inp, inp, inp])\n    base_model = efn.EfficientNetB5(\n        input_shape=input_shape,\n        weights=None,\n        include_top=False,\n    )\n#     base_model.load_weights(\"../input/efficientnet/efficientnet-b5_weights_tf_dim_ordering_tf_kernels_autoaugment_notop.h5\")\n    x1 = base_model(inp)\n    x1 = GlobalAveragePooling2D()(x1)\n    \n    inp2 = Input(shape=(6,))\n    \n    x = Concatenate()([x1, inp2])\n    output = Dense(output_shape, activation='linear', name='output')(x)\n    model = Model(inputs=[inp, inp2],\n                  outputs=output,\n                  name='effnetb5_plus_osic')\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def efficientnet_model(input_shape, output_shape):\n    \"\"\"\n    Although input shape can be modified, I prefer to use\n    efficientnet native input shape from model.\n    \"\"\"\n    inp = Input(shape=input_shape)\n#     inp_concat = Concatenate()([inp, inp, inp])\n    base_model = efn.EfficientNetB5(\n        input_shape=input_shape,\n        weights=None, \n        include_top=False,\n    )\n#     base_model.load_weights(\"../input/efficientnet/efficientnet-b5_weights_tf_dim_ordering_tf_kernels_autoaugment_notop.h5\")\n    x = base_model(inp)\n    x = GlobalAveragePooling2D()(x)\n    output = Dense(output_shape, activation='linear', name='output')(x)\n    model = Model(inputs=inp, outputs=output, name='effnetb5_osic')\n    \n    # On-off trainable layers\n#     model.layers[0].trainable = False\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def dense_plus_model1(input_shape, output_shape):\n    \"\"\"\n    Multiple quantile regression basic dense only model for image input.\n    \n    Parameters\n    ----------\n    input_shape : tuple\n        except input to be in grayscale\n    output_shape : int\n        how many quantile to fit\n        \n    Returns\n    -------\n    model : tf.keras.models.Model\n    \"\"\"\n    input = Input(shape=input_shape, name=\"input\")\n    x = Conv2D(16, (3, 3), activation='relu', padding='same')(input)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(32, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(64, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(128, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Flatten()(x)\n    x = Dense(512, activation='relu', name='dense1')(x)\n    x = Dense(128, activation='relu', name='dense2')(x)\n    x = Dense(64, activation='relu', name='dense3')(x)\n    x = Dense(16, activation='relu', name='dense4')(x)\n    \n    input2 = Input(shape=(6,))\n    x = Concatenate()([x, input2])\n    output = Dense(output_shape, activation='linear', name='output')(x)\n    \n    model = Model(inputs=[input, input2], outputs=output, name='dense_plus_model1')\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def dense_model1(input_shape, output_shape):\n    \"\"\"\n    Multiple quantile regression basic dense only model for image input.\n    \n    Parameters\n    ----------\n    input_shape : tuple\n        except input to be in grayscale\n    output_shape : int\n        how many quantile to fit\n        \n    Returns\n    -------\n    model : tf.keras.models.Model\n    \"\"\"\n    input = Input(shape=input_shape, name=\"input\")\n    x = Conv2D(16, (3, 3), activation='relu', padding='same')(input)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(32, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(64, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Conv2D(128, (3, 3), activation='relu', padding='same')(x)\n    x = MaxPooling2D((2, 2))(x)\n    x = Flatten()(x)\n    x = Dense(512, activation='relu', name='dense1')(x)\n    x = Dense(128, activation='relu', name='dense2')(x)\n    x = Dense(64, activation='relu', name='dense3')(x)\n    x = Dense(16, activation='relu', name='dense4')(x)\n    output = Dense(output_shape, activation='linear', name='output')(x)\n    \n    model = Model(inputs=input, outputs=output, name='dense_model1')\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Data Generator\noutput_shape = (512, 512)\n# train_image_dataset = OSICImageGenerator(X_train, y_train, bbox_map, output_shape=output_shape, segmented=True, shuffle=True, debug=False)\n# valid_image_dataset = OSICImageGenerator(X_valid, y_valid, bbox_map, output_shape=output_shape, segmented=True, shuffle=False)\n\ntrain_image_dataset = OSIC_CNNMLP_ImageGenerator(\n    np.array(X_train), np.array(y_train),\n    bbox_map,\n    output_shape=output_shape, batch_size=4, segmented=True, shuffle=True, debug=False\n)\nvalid_image_dataset = OSIC_CNNMLP_ImageGenerator(\n    np.array(X_valid), np.array(y_valid),\n    bbox_map,\n    output_shape=output_shape, batch_size=4, segmented=True, shuffle=False\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Learning Rate Scheduler\nLR_START = 1e-3\nLR_MAX = 0.03\nLR_RAMPUP_EPOCHS = 5\nLR_SUSTAIN_EPOCHS = 0\nLR_STEP_DECAY = 0.75\n\ndef lrfn(epoch):\n    if epoch < LR_RAMPUP_EPOCHS:\n        lr = (LR_MAX - LR_START) / LR_RAMPUP_EPOCHS * epoch + LR_START\n    elif epoch < LR_RAMPUP_EPOCHS + LR_SUSTAIN_EPOCHS:\n        lr = LR_MAX\n    else:\n        lr = LR_MAX * LR_STEP_DECAY**((epoch - LR_RAMPUP_EPOCHS - LR_SUSTAIN_EPOCHS)//10)\n    return lr\n    \nlr2 = LearningRateScheduler(lrfn, verbose = True)\n\nrng = [i for i in range(100)]\ny = [lrfn(x) for x in rng]\nplt.plot(rng, y); \nplt.xlabel('epoch',size=14); plt.ylabel('learning rate',size=14)\nplt.title('Training Schedule',size=16); plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# model = dense_model1(input_shape=(512, 512, 1), output_shape=len(quantiles))\n# model = dense_plus_model1(input_shape=(456, 456, 1), output_shape=len(quantiles))\n# model = efficientnet_model(input_shape=(456, 456, 1), output_shape=len(quantiles))\n# model = efficientnet_plus_model(input_shape=(456, 456, 1), output_shape=len(quantiles))\nmodel = effnet_mlp_model(input_shape=(512, 512, 1), output_shape=len(quantiles))\nmodel.compile(loss=mloss(0.65),\n              optimizer=Adam(learning_rate=0.03),\n              metrics=[score])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# model_path = \"dense_model1.h5\"\n# model_path = \"dense_plus_model1.h5\"\n# model_path = \"effnetb5.h5\"\n# model_path = \"effnetb5_plus.h5\"\nmodel_path = \"cnn_mlp_only_mid.h5\"\n\ncheckpoint_path = os.path.abspath(model_path)\n\ncp_callback = ModelCheckpoint(filepath=checkpoint_path,\n                              monitor='val_loss',\n                              save_best_only=True,\n                              mode='min',\n                              save_weights_only=True),\nes_callback = EarlyStopping(monitor='val_loss',\n                            min_delta=1e-4,\n                            patience=4,\n                            mode='min'),\nlr_callback = lr2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = model.fit(\n    train_image_dataset,\n    validation_data=valid_image_dataset,\n    initial_epoch=0,\n    callbacks=[\n        cp_callback,\n        es_callback,\n        lr_callback\n    ],\n    epochs=1  # estimated ~30 min per epoch, so maximum runtime of 6 hours equals 12 epoch, here I just use 50%\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_score = score(np.expand_dims(y_valid, axis=-1), model.predict(valid_image_dataset))\nprint(f\"Validation score: {valid_score}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create Submission"},{"metadata":{},"cell_type":"markdown","source":"## Model 1 (LGBM, Tabular)"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df = pd.read_csv(f\"{basepath}sample_submission.csv\")\nsubmission_df = submission_df[['Patient_Week']].merge(\n    test_dataset.df[['Patient_Week', 'FVC_pred', 'Confidence']]\\\n        .rename(columns={\"FVC_pred\": \"FVC\"}),\n    how='inner',\n    on='Patient_Week'\n)\nprint(submission_df.head(5))\nsubmission_df.to_csv(\"submission_lgbm.csv\", header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model 2 (CNN-MLP, Image+Tabular)"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_pred_num = ['FVC_base', 'Percent_base', 'Age_base', 'Weeks_passed']\ncols_pred_cat = ['Sex', 'SmokingStatus']\ncols_pred = cols_pred_num + cols_pred_cat","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = get_X_y_cnnmlp(base_test_df)\ntest_image_dataset = OSIC_CNNMLP_ImageGenerator(\n    np.array(X_test), np.array([2000]*len(X_test)),  # dummy target data\n    bbox_map,\n    output_shape=output_shape, batch_size=4, segmented=True, shuffle=False\n)\n\ny_pred = model.predict(test_image_dataset)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base_test_df['FVC'] = y_pred[:, 1]\nbase_test_df['Confidence'] = y_pred[:, 2] - y_pred[:, 0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df = pd.read_csv(f\"{basepath}sample_submission.csv\")\nsubmission_df = submission_df[['Patient_Week']].merge(\n    base_test_df[['Patient_Week', 'FVC', 'Confidence']],\n    how='inner',\n    on='Patient_Week'\n)\nprint(submission_df.head(5))\nsubmission_df.to_csv(\"submission_cnn_mlp.csv\", header=True, 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}