{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport torch\nimport torchvision\nimport torchvision.transforms as transforms\nimport torch.nn as nn\n\n\n\nimport seaborn as sns\nimport cv2\nfrom skimage import io\n\nfrom sklearn.model_selection import train_test_split\nimport tensorflow as tf\nfrom tensorflow.python.keras import Sequential\nfrom tensorflow.keras import layers, optimizers\nfrom tensorflow.keras.models import Model, load_model\nfrom tensorflow.keras.initializers import glorot_uniform\nfrom tensorflow.keras.utils import plot_model\nfrom tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint, LearningRateScheduler\nimport tensorflow.keras.backend as K\n\nfrom tensorflow.keras.layers import Conv2D, BatchNormalization, Activation, MaxPool2D, Conv2DTranspose, Concatenate, Input\nfrom tensorflow.keras.applications import VGG19\nfrom tensorflow.keras.layers import (Dense, Dropout, Activation, Flatten, Input, Add,\n                                    BatchNormalization, LeakyReLU, Concatenate, GlobalAveragePooling2D,Conv2D, AveragePooling2D)\n\nfrom warnings import filterwarnings\nfilterwarnings('ignore')\n\nimport random\n\nimport glob\nfrom IPython.display import display\n\nfrom pathlib import Path\n\n# Set Color Palettes for the notebook\ncustom_colors = ['#74a09e','#86c1b2','#98e2c6','#f3c969','#f2a553', '#d96548', '#c14953']\nsns.palplot(sns.color_palette(custom_colors))\n\n# Set Style\nsns.set_style(\"whitegrid\")\nsns.despine(left=True, bottom=True)\n\nfrom scipy.stats import pearsonr\nimport pydicom\nimport re\nfrom sklearn.cluster import KMeans\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom tensorflow.keras.utils import Sequence","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:47.870048Z","iopub.execute_input":"2022-06-26T14:07:47.870416Z","iopub.status.idle":"2022-06-26T14:07:47.963178Z","shell.execute_reply.started":"2022-06-26T14:07:47.870389Z","shell.execute_reply":"2022-06-26T14:07:47.962071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"config = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:48.655969Z","iopub.execute_input":"2022-06-26T14:07:48.656381Z","iopub.status.idle":"2022-06-26T14:07:48.682569Z","shell.execute_reply.started":"2022-06-26T14:07:48.656348Z","shell.execute_reply":"2022-06-26T14:07:48.680562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ROOT = Path('../input/osic-pulmonary-fibrosis-progression/')\n\ntrain = pd.read_csv(ROOT / 'train.csv')\ntest = pd.read_csv(ROOT / 'test.csv')\nsub = pd.read_csv(ROOT / 'sample_submission.csv')\n\ntrain.info()","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:49.375892Z","iopub.execute_input":"2022-06-26T14:07:49.376809Z","iopub.status.idle":"2022-06-26T14:07:49.432702Z","shell.execute_reply.started":"2022-06-26T14:07:49.376772Z","shell.execute_reply":"2022-06-26T14:07:49.431818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:49.777872Z","iopub.execute_input":"2022-06-26T14:07:49.778310Z","iopub.status.idle":"2022-06-26T14:07:49.797949Z","shell.execute_reply.started":"2022-06-26T14:07:49.778270Z","shell.execute_reply":"2022-06-26T14:07:49.797121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train.loc[:,\"Patient\"]","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:50.042891Z","iopub.execute_input":"2022-06-26T14:07:50.043328Z","iopub.status.idle":"2022-06-26T14:07:50.047300Z","shell.execute_reply.started":"2022-06-26T14:07:50.043294Z","shell.execute_reply":"2022-06-26T14:07:50.046600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:50.308397Z","iopub.execute_input":"2022-06-26T14:07:50.308972Z","iopub.status.idle":"2022-06-26T14:07:50.313261Z","shell.execute_reply.started":"2022-06-26T14:07:50.308929Z","shell.execute_reply":"2022-06-26T14:07:50.312212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#data.iloc[<row s <column selection>]\n# train.iloc[2,1]","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:50.575521Z","iopub.execute_input":"2022-06-26T14:07:50.576071Z","iopub.status.idle":"2022-06-26T14:07:50.580010Z","shell.execute_reply.started":"2022-06-26T14:07:50.576035Z","shell.execute_reply":"2022-06-26T14:07:50.579128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Q: Are there any missing values?\", \"\\n\" +\n      \"A: {}\".format(train.isnull().values.any()))","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:50.819958Z","iopub.execute_input":"2022-06-26T14:07:50.820488Z","iopub.status.idle":"2022-06-26T14:07:50.827677Z","shell.execute_reply.started":"2022-06-26T14:07:50.820457Z","shell.execute_reply":"2022-06-26T14:07:50.826765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# **EDA**","metadata":{}},{"cell_type":"code","source":"print(\"There are {} unique patients in Train Data.\".format(len(train[\"Patient\"].unique())), \"\\n\")\n\n# Recordings per Patient\ndata = train.groupby(by=\"Patient\")[\"Weeks\"].count().reset_index(drop=False)\n# Sort by Weeks\ndata = data.sort_values(['Weeks']).reset_index(drop=True)\nprint(\"Minimum number of entries are: {}\".format(data[\"Weeks\"].min()), \"\\n\" +\n      \"Maximum number of entries are: {}\".format(data[\"Weeks\"].max()))\n\n# Plot\nplt.figure(figsize = (16, 6))\np = sns.barplot(data[\"Patient\"], data[\"Weeks\"], color=custom_colors[2])\n\nplt.title(\"Number of Entries per Patient\", fontsize = 17)\nplt.xlabel('Patient', fontsize=14)\nplt.ylabel('Frequency', fontsize=14)\n\np.axes.get_xaxis().set_visible(False);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:51.295402Z","iopub.execute_input":"2022-06-26T14:07:51.296324Z","iopub.status.idle":"2022-06-26T14:07:52.428683Z","shell.execute_reply.started":"2022-06-26T14:07:51.296280Z","shell.execute_reply":"2022-06-26T14:07:52.427834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Select unique bio info for the patients\ndata = train.groupby(by=\"Patient\")[[\"Patient\", \"Age\", \"Sex\", \"SmokingStatus\"]].first().reset_index(drop=True)\n\n# Figure\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (16, 6))\n\na = sns.distplot(data[\"Age\"], ax=ax1, color=custom_colors[1], hist=False, kde_kws=dict(lw=6, ls=\"--\"))\nb = sns.countplot(data[\"Sex\"], ax=ax2, palette=custom_colors[2:4])\nc = sns.countplot(data[\"SmokingStatus\"], ax=ax3, palette = custom_colors[4:7])\n\na.set_title(\"Patient Age Distribution\", fontsize=16)\nb.set_title(\"Sex Frequency\", fontsize=16)\nc.set_title(\"Smoking Status\", fontsize=16);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:52.430263Z","iopub.execute_input":"2022-06-26T14:07:52.430611Z","iopub.status.idle":"2022-06-26T14:07:52.911789Z","shell.execute_reply.started":"2022-06-26T14:07:52.430579Z","shell.execute_reply":"2022-06-26T14:07:52.911100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Min FVC value: {:,}\".format(train[\"FVC\"].min()), \"\\n\" +\n      \"Max FVC value: {:,}\".format(train[\"FVC\"].max()), \"\\n\" +\n      \"\\n\" +\n      \"Min Percent value: {:.4}%\".format(train[\"Percent\"].min()), \"\\n\" +\n      \"Max Percent value: {:.4}%\".format(train[\"Percent\"].max()))\n\n# Figure\nf, (ax1, ax2) = plt.subplots(1, 2, figsize = (16, 6))\n\na = sns.distplot(train[\"FVC\"], ax=ax1, color=custom_colors[6], hist=False, kde_kws=dict(lw=6, ls=\"--\"))\nb = sns.distplot(train[\"Percent\"], ax=ax2, color=custom_colors[4], hist=False, kde_kws=dict(lw=6, ls=\"-.\"))\n\na.set_title(\"FVC Distribution\", fontsize=16)\nb.set_title(\"Percent Distribution\", fontsize=16);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:52.912869Z","iopub.execute_input":"2022-06-26T14:07:52.913302Z","iopub.status.idle":"2022-06-26T14:07:53.384823Z","shell.execute_reply.started":"2022-06-26T14:07:52.913273Z","shell.execute_reply":"2022-06-26T14:07:53.383886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Minimum no. weeks before CT: {}\".format(train['Weeks'].min()), \"\\n\" +\n      \"Maximum no. weeks after CT: {}\".format(train['Weeks'].max()))\n\nplt.figure(figsize = (16, 6))\n\na = sns.distplot(train['Weeks'], color=custom_colors[3], hist=False, kde_kws=dict(lw=8, ls=\"--\"))\nplt.title(\"Number of weeks before/after the CT scan\", fontsize = 16)\nplt.xlabel(\"Weeks\", fontsize=14);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:53.387362Z","iopub.execute_input":"2022-06-26T14:07:53.387731Z","iopub.status.idle":"2022-06-26T14:07:53.680121Z","shell.execute_reply.started":"2022-06-26T14:07:53.387699Z","shell.execute_reply":"2022-06-26T14:07:53.679213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute Correlation\ncorr1, _ = pearsonr(train[\"FVC\"], train[\"Percent\"])\ncorr2, _ = pearsonr(train[\"FVC\"], train[\"Age\"])\ncorr3, _ = pearsonr(train[\"Percent\"], train[\"Age\"])\nprint(\"Pearson Corr FVC x Percent: {:.4}\".format(corr1), \"\\n\" +\n      \"Pearson Corr FVC x Age: {:.0}\".format(corr2), \"\\n\" +\n      \"Pearson Corr Percent x Age: {:.2}\".format(corr3))\n\n# Figure\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (16, 6))\n\na = sns.scatterplot(x = train[\"FVC\"], y = train[\"Percent\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax1)\n\nb = sns.scatterplot(x = train[\"FVC\"], y = train[\"Age\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax2)\n\nc = sns.scatterplot(x = train[\"Percent\"], y = train[\"Age\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax3)\n\na.set_title(\"Correlation between FVC and Percent\", fontsize = 16)\na.set_xlabel(\"FVC\", fontsize = 14)\na.set_ylabel(\"Percent\", fontsize = 14)\n\nb.set_title(\"Correlation between FVC and Age\", fontsize = 16)\nb.set_xlabel(\"FVC\", fontsize = 14)\nb.set_ylabel(\"Age\", fontsize = 14)\n\nc.set_title(\"Correlation between Percent and Age\", fontsize = 16)\nc.set_xlabel(\"Percent\", fontsize = 14)\nc.set_ylabel(\"Age\", fontsize = 14);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:53.681324Z","iopub.execute_input":"2022-06-26T14:07:53.681655Z","iopub.status.idle":"2022-06-26T14:07:54.614786Z","shell.execute_reply.started":"2022-06-26T14:07:53.681626Z","shell.execute_reply":"2022-06-26T14:07:54.614059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Figure\nf, (ax1, ax2) = plt.subplots(1,2, figsize = (16, 6))\n\na = sns.barplot(x = train[\"SmokingStatus\"], y = train[\"FVC\"], ax=ax1, palette=custom_colors[0:4])\nb = sns.barplot(x = train[\"SmokingStatus\"], y = train[\"Percent\"], ax=ax2, palette=custom_colors[4:7])\n\na.set_title(\"Mean FVC per Smoking Status\", fontsize=16)\nb.set_title(\"Mean Perc per Smoking Status\", fontsize=16);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:54.615710Z","iopub.execute_input":"2022-06-26T14:07:54.616511Z","iopub.status.idle":"2022-06-26T14:07:55.106084Z","shell.execute_reply.started":"2022-06-26T14:07:54.616477Z","shell.execute_reply":"2022-06-26T14:07:55.105183Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='Sex', y='FVC', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.107398Z","iopub.execute_input":"2022-06-26T14:07:55.107806Z","iopub.status.idle":"2022-06-26T14:07:55.247102Z","shell.execute_reply.started":"2022-06-26T14:07:55.107766Z","shell.execute_reply":"2022-06-26T14:07:55.246181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='Sex', y='Age', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.248598Z","iopub.execute_input":"2022-06-26T14:07:55.249498Z","iopub.status.idle":"2022-06-26T14:07:55.392480Z","shell.execute_reply.started":"2022-06-26T14:07:55.249450Z","shell.execute_reply":"2022-06-26T14:07:55.390995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='SmokingStatus', y='FVC', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.393962Z","iopub.execute_input":"2022-06-26T14:07:55.394560Z","iopub.status.idle":"2022-06-26T14:07:55.557547Z","shell.execute_reply.started":"2022-06-26T14:07:55.394518Z","shell.execute_reply":"2022-06-26T14:07:55.556652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='Sex', y='Percent', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.560793Z","iopub.execute_input":"2022-06-26T14:07:55.561418Z","iopub.status.idle":"2022-06-26T14:07:55.699344Z","shell.execute_reply.started":"2022-06-26T14:07:55.561375Z","shell.execute_reply":"2022-06-26T14:07:55.698269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='Sex', y='Weeks', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.700734Z","iopub.execute_input":"2022-06-26T14:07:55.701165Z","iopub.status.idle":"2022-06-26T14:07:55.845966Z","shell.execute_reply.started":"2022-06-26T14:07:55.701122Z","shell.execute_reply":"2022-06-26T14:07:55.844812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(x='SmokingStatus', y='Age', data=train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:55.847499Z","iopub.execute_input":"2022-06-26T14:07:55.848039Z","iopub.status.idle":"2022-06-26T14:07:56.019793Z","shell.execute_reply.started":"2022-06-26T14:07:55.847994Z","shell.execute_reply":"2022-06-26T14:07:56.018799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create Time variable to count in ascending order the times the Patient has done a check in FVC\ndata_time = train.groupby(by=\"Patient\")[\"Weeks\"].count().reset_index()\n# print(data_time)\ntrain[\"Time\"] = 0\n#print(train)\n\nfor patient, times in zip(data_time[\"Patient\"], data_time[\"Weeks\"]):\n    train.loc[train[\"Patient\"] == patient, 'Time'] = range(1, times+1)\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:56.021220Z","iopub.execute_input":"2022-06-26T14:07:56.021735Z","iopub.status.idle":"2022-06-26T14:07:56.176570Z","shell.execute_reply.started":"2022-06-26T14:07:56.021689Z","shell.execute_reply":"2022-06-26T14:07:56.175484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train)","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:56.178119Z","iopub.execute_input":"2022-06-26T14:07:56.179029Z","iopub.status.idle":"2022-06-26T14:07:56.193683Z","shell.execute_reply.started":"2022-06-26T14:07:56.178965Z","shell.execute_reply":"2022-06-26T14:07:56.192890Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# For graph purposes, keep only Patients that had a big difference in FVC between Time 1 and last Time\nmin_fvc = train[train[\"Time\"] == 1][[\"Patient\", \"FVC\"]].reset_index(drop=True)\n\nidx = train.groupby([\"Patient\"])[\"Weeks\"].transform(max) == train[\"Weeks\"]\nmax_fvc = train[idx][[\"Patient\", \"FVC\"]].reset_index(drop=True)\n# print(max_fvc)\n\n# Compute difference and select only top patients with biggest difference\ndata = pd.merge(min_fvc, max_fvc, how=\"inner\", on=\"Patient\")\ndata[\"Dif\"] = data[\"FVC_x\"] - data[\"FVC_y\"]\n\n# Select only top n\nl = list(data.sort_values(\"Dif\", ascending=False).head(10)[\"Patient\"])\nx = train[train[\"Patient\"].isin(l)]","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:56.194734Z","iopub.execute_input":"2022-06-26T14:07:56.195246Z","iopub.status.idle":"2022-06-26T14:07:56.221052Z","shell.execute_reply.started":"2022-06-26T14:07:56.195212Z","shell.execute_reply":"2022-06-26T14:07:56.220029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (16, 6))\n\na = sns.lineplot(x = x[\"Time\"], y = x[\"FVC\"], hue = x[\"Patient\"], legend=False,\n                 palette=sns.color_palette(\"GnBu_d\", 10), size=1)\n\nplt.title(\"Patient FVC decrease on Weeks\", fontsize = 16)\nplt.xlabel(\"Weeks\", fontsize=14)\nplt.ylabel(\"FVC\", fontsize=14);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:56.223069Z","iopub.execute_input":"2022-06-26T14:07:56.224093Z","iopub.status.idle":"2022-06-26T14:07:56.515486Z","shell.execute_reply.started":"2022-06-26T14:07:56.224049Z","shell.execute_reply":"2022-06-26T14:07:56.514820Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create base director for Train .dcm files\ndirector = \"../input/osic-pulmonary-fibrosis-progression/train\"\n\n# Create path column with the path to each patient's CT\ntrain[\"Path\"] = director + \"/\" + train[\"Patient\"]\n\n# Create variable that shows how many CT scans each patient has\ntrain[\"CT_number\"] = 0\n\nfor k, path in enumerate(train[\"Path\"]):\n    train[\"CT_number\"][k] = len(os.listdir(path))","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:07:56.516673Z","iopub.execute_input":"2022-06-26T14:07:56.517163Z","iopub.status.idle":"2022-06-26T14:08:02.246534Z","shell.execute_reply.started":"2022-06-26T14:07:56.517129Z","shell.execute_reply":"2022-06-26T14:08:02.245580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Minimum number of CT scans: {}\".format(train[\"CT_number\"].min()), \"\\n\" +\n      \"Maximum number of CT scans: {:,}\".format(train[\"CT_number\"].max()))\n\n# Scans per Patient\ndata = train.groupby(by=\"Patient\")[\"CT_number\"].first().reset_index(drop=False)\n# Sort by Weeks\ndata = data.sort_values(['CT_number']).reset_index(drop=True)\n\n\n# Plot\nplt.figure(figsize = (16, 6))\np = sns.barplot(data[\"Patient\"], data[\"CT_number\"], color=custom_colors[5])\nplt.axvline(x=85, color=custom_colors[2], linestyle='--', lw=3)\n\nplt.title(\"Number of CT Scans per Patient\", fontsize = 17)\nplt.xlabel('Patient', fontsize=14)\nplt.ylabel('Frequency', fontsize=14)\n\nplt.text(86, 850, \"Median=94\", fontsize=13)\n\np.axes.get_xaxis().set_visible(False);","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:08:02.247695Z","iopub.execute_input":"2022-06-26T14:08:02.248360Z","iopub.status.idle":"2022-06-26T14:08:03.470117Z","shell.execute_reply.started":"2022-06-26T14:08:02.248326Z","shell.execute_reply":"2022-06-26T14:08:03.469035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(data.Patient.iloc[-3])","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:08:03.471900Z","iopub.execute_input":"2022-06-26T14:08:03.472343Z","iopub.status.idle":"2022-06-26T14:08:03.478729Z","shell.execute_reply.started":"2022-06-26T14:08:03.472303Z","shell.execute_reply":"2022-06-26T14:08:03.477815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data_time=","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:08:41.328372Z","iopub.execute_input":"2022-06-26T14:08:41.328768Z","iopub.status.idle":"2022-06-26T14:08:41.332560Z","shell.execute_reply.started":"2022-06-26T14:08:41.328734Z","shell.execute_reply":"2022-06-26T14:08:41.331932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/19.dcm\"\ndataset = pydicom.dcmread(path)\n\nprint(\"Patient id.......:\", dataset.PatientID, \"\\n\" +\n      \"Modality.........:\", dataset.Modality, \"\\n\" +\n      \"Rows.............:\", dataset.Rows, \"\\n\" +\n      \"Columns..........:\", dataset.Columns)\n\nplt.figure(figsize = (7, 7))\nplt.imshow(dataset.pixel_array, cmap=\"plasma\")\nplt.axis('off');","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:08:41.600652Z","iopub.execute_input":"2022-06-26T14:08:41.601244Z","iopub.status.idle":"2022-06-26T14:08:41.874928Z","shell.execute_reply.started":"2022-06-26T14:08:41.601211Z","shell.execute_reply":"2022-06-26T14:08:41.874289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/ID00173637202238329754031\"\ndatasets = []\n\n# First Order the files in the dataset\nfiles = []\nfor dcm in list(os.listdir(patient_dir)):\n    files.append(dcm) \nfiles.sort(key=lambda f: int(re.sub('\\D', '', f)))\nprint(files)\n\n# Read in the Dataset\nfor dcm in files:\n    path = patient_dir + \"/\" + dcm\n    datasets.append(pydicom.dcmread(path))\n\n# Plot the images\nfig=plt.figure(figsize=(50,50))\ncolumns = 15\nrows = 40\n\nfor i in range(1, columns*rows +1):\n    img = datasets[i-1].pixel_array\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=\"plasma\")\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","metadata":{"execution":{"iopub.status.busy":"2022-06-26T14:08:41.876287Z","iopub.execute_input":"2022-06-26T14:08:41.876670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.loc[train['Patient'] == 'ID00078637202199415319443']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data[\"Patient\"]==\"ID00078637202199415319443\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from PIL import Image\nfrom IPython.display import Image as show_gif\nimport scipy.misc\nimport matplotlib","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_gif(number_of_CT = 87):\n    \"\"\"Picks a patient at random and creates a GIF with their CT scans.\"\"\"\n    \n    # Select one of the patients\n    # patient = \"ID00007637202177411956430\"\n    patient = train[train[\"CT_number\"] == number_of_CT].sample(random_state=1)[\"Patient\"].values[0]\n    \n    # === READ IN .dcm FILES ===\n    patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/\" + patient\n    datasets = []\n\n    # First Order the files in the dataset\n    files = []\n    for dcm in list(os.listdir(patient_dir)):\n        files.append(dcm) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Read in the Dataset from the Patient path\n    for dcm in files:\n        path = patient_dir + \"/\" + dcm\n        datasets.append(pydicom.dcmread(path))\n        \n        \n    # === SAVE AS .png ===\n    # Create directory to save the png files\n    if os.path.isdir(f\"png_{patient}\") == False:\n        os.mkdir(f\"png_{patient}\")\n\n    # Save images to PNG\n    for i in range(len(datasets)):\n        img = datasets[i].pixel_array\n        matplotlib.image.imsave(f'png_{patient}/img_{i}.png', img)\n        \n        \n    # === CREATE GIF ===\n    # First Order the files in the dataset (again)\n    files = []\n    for png in list(os.listdir(f\"../working/png_{patient}\")):\n        files.append(png) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Create the frames\n    frames = []\n\n    # Create frames\n    for file in files:\n    #     print(\"../working/png_images/\" + name)\n        new_frame = Image.open(f\"../working/png_{patient}/\" + file)\n        frames.append(new_frame)\n\n    # Save into a GIF file that loops forever\n    frames[0].save(f'gif_{patient}.gif', format='GIF',\n                   append_images=frames[1:],\n                   save_all=True,\n                   duration=200, loop=0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"create_gif(number_of_CT=12)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_gif(filename=\"./gif_ID00165637202237320314458.gif\", format='png', width=400, height=400)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    \n    # Find the average pixel value near the lungs\n        # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    \n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    \n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    \n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    \n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n    return mask*img","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Select a sample\npath = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/19.dcm\"\ndataset = pydicom.dcmread(path)\nimg = dataset.pixel_array\n\n# Masked image\nmask_img = make_lungmask(img, display=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# masking image for one patient","metadata":{}},{"cell_type":"code","source":"patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430\"\ndatasets = []\n\n# First Order the files in the dataset\nfiles = []\nfor dcm in list(os.listdir(patient_dir)):\n    files.append(dcm) \nfiles.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n# Read in the Dataset\nfor dcm in files:\n    path = patient_dir + \"/\" + dcm\n    datasets.append(pydicom.dcmread(path))\n    \nimgs = []\nfor data in datasets:\n    img = data.pixel_array\n    imgs.append(img)\n    \n    \n# Show masks\nfig=plt.figure(figsize=(16, 6))\ncolumns = 10\nrows = 3\n\nfor i in range(1, columns*rows +1):\n    img = make_lungmask(datasets[i-1].pixel_array)\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=\"gray\")\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# data preparation","metadata":{}},{"cell_type":"code","source":"def get_tab(df):\n    ''' \n    This function gives an array wrt each patient containing\n    feature like age, gender and smoking status\n    '''\n    vector = [(df.Age.values[0]-30)/30]\n    \n    if df.Sex.values[0].lower() == 'male':\n        vector.append(0)\n    else:\n        vector.append(1)\n        \n    if df.SmokingStatus.values[0] == 'Never smoked':\n        vector.extend([0,0])\n    elif df.SmokingStatus.values[0] == 'Ex-smoker':\n        vector.extend([1,1])\n    elif df.SmokingStatus.values[0] == 'Currently smokes':\n        vector.extend([0,1])\n    else:\n        vector.extend([1,0])\n        \n    return np.array(vector)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"A = {} #Stores slope value for each of the patient\nTAB = {} #Stores training data wrt each patient\nP = [] #Stores all unique patient id's\n\nfor i,p in enumerate(train.Patient.unique()):\n    sub = train.loc[train.Patient == p, :]\n    fvc = sub.FVC.values\n    week = sub.Weeks.values\n    #print(week)\n    c = np.vstack([week, np.ones(len(week))]).T\n    a, b = np.linalg.lstsq(c,fvc)[0]\n    #print(b)\n    \n    A[p] = a # Contains slope\n    TAB[p] = get_tab(sub) #Contains gender and smoking feature\n    P.append(p) #contains unique id\n   # print(TAB)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Creating CNN architecture for coeficient prediction:","metadata":{}},{"cell_type":"code","source":"def get_img(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize(d.pixel_array/2**11 ,(512,512))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class IGenerator(Sequence):\n    \n    ''' \n    This is the generator class, which generates an input of batch size 32\n    i.e 32 patient's 2 dicom image, and features from tabular data is generated. As output \n    from his generator x and y contains pixel_data of a dicom image, tab conatins patient's meta\n    information, and 'a' is the coeffiecient wrt each patient. \n    '''\n    BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n    def __init__(self, keys, a, tab, batch_size=16):\n        self.keys = [k for k in keys if k not in self.BAD_ID]\n        self.a = a\n        self.tab = tab\n        self.batch_size = batch_size\n        \n        self.train_data = {}\n        for p in train.Patient.values:\n            self.train_data[p] = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n            #print(p)\n    def __len__(self):\n        return 1000\n    \n    def __getitem__(self, idx):\n        x, y = [], []\n        a, tab = [], [] \n        keys = np.random.choice(self.keys, size = self.batch_size)\n        \n        for k in keys:\n            try:\n                i = np.random.choice(self.train_data[k], size=1)[0]\n                j = np.random.choice(self.train_data[k], size=1)[0]\n                img1 = get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{k}/{i}')\n                img2 = get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{k}/{j}')\n                \n                x.append(img1)\n                y.append(img2)\n                \n                a.append(self.a[k])\n                tab.append(self.tab[k])\n            except:\n                print(k, i)\n        \n        \n        x,y,a,tab = np.array(x),np.array(y), np.array(a), np.array(tab)\n        x = np.expand_dims(x, axis=-1)\n        y = np.expand_dims(y, axis=-1)\n        return [x,y, tab] , a\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def model_architechture(shape=(512,512,1)):\n    '''Architecture used here is inspired by this kaggle notebook \n    https://www.kaggle.com/miklgr500/linear-decay-based-on-resnet-cnn/notebook'''\n    \n    def res_block(x, filter_number):\n        _x = x\n        x = Conv2D(filter_number, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n        x = BatchNormalization()(x)\n        x = LeakyReLU(0.05)(x)\n        x = Conv2D(filter_number, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n        x = LeakyReLU(0.05)(x)\n        \n        x = Add()([_x, x])\n        return x\n    \n    #two input branch for images\n    input1 = Input(shape=shape, name= 'dicom_image_1')\n    input2 = Input(shape=shape, name= 'dicom_image_2')\n    \n    #image input branch 1 begins\n    x = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(input1)\n    x = BatchNormalization()(x)\n    x = LeakyReLU(0.05)(x)\n    \n    x = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU(0.05)(x)\n    \n    x = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x)\n    \n    #image input branch 2 begins\n    y = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(input2)\n    y = BatchNormalization()(y)\n    y = LeakyReLU(0.05)(y)\n    \n    y = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(y)\n    y = BatchNormalization()(y)\n    y = LeakyReLU(0.05)(y)\n    \n    y = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(y)\n    \n    #Concatinating image inputs\n    x_and_y = Concatenate()([x, y])\n    \n    x_and_y = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(x_and_y)\n    x_and_y = BatchNormalization()(x_and_y)\n    x_and_y = LeakyReLU(0.05)(x_and_y)\n    \n    x_and_y = Conv2D(16, kernel_size=(3, 3), strides=(1, 1), padding='same')(x_and_y)\n    for _ in range(2):\n        x_and_y = res_block(x_and_y, 16)\n    x_and_y = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x_and_y)\n    \n    x_and_y = Conv2D(64, kernel_size=(3, 3), strides=(1, 1), padding='same')(x_and_y)\n    for _ in range(3):\n        x_and_y = res_block(x_and_y, 64)\n    x_and_y = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x_and_y)    \n    \n    x_and_y = Conv2D(128, kernel_size=(3, 3), strides=(1, 1), padding='same')(x_and_y)\n    for _ in range(1):\n        x_and_y = res_block(x_and_y, 128)\n        \n   \n    x_and_y = GlobalAveragePooling2D()(x_and_y)\n    \n    #Patient tabular data input\n    input3 = Input(shape=(4,))\n    z = tf.keras.layers.GaussianNoise(0.2)(input3)\n    xyz = Concatenate()([x_and_y, z])\n    xyz = Dropout(0.6)(xyz) \n    xyz = Dense(1)(xyz)\n    return Model([input1, input2, input3] , xyz)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = model_architechture()\n\nmodel.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=0.001), loss='mae') \n\ntr_p, vl_p = train_test_split(P, shuffle=True, train_size= 0.8)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow import keras\nkeras.utils.plot_model(model,'img.png', show_shapes=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"er = tf.keras.callbacks.EarlyStopping(\n    monitor=\"val_loss\",\n    min_delta=1e-3,\n    patience=5,\n    verbose=0,\n    mode=\"auto\",\n    baseline=None,\n    restore_best_weights=True,\n)\n\nmodel.fit_generator(IGenerator(keys=tr_p, \n                               a = A, \n                               tab = TAB), \n                    steps_per_epoch = 20,\n                    validation_data=IGenerator(keys=vl_p, \n                               a = A, \n                               tab = TAB),\n                    validation_steps = 20, \n                    callbacks = [er], \n                    epochs=30)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.save('best_model.h5')\n            ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def score(fvc_true, fvc_pred, sigma):\n    sigma_clip = np.maximum(sigma,70)\n    delta = np.abs(fvc_true - fvc_pred)\n    delta = np.minimum(delta,1000)\n    sqrt = np.sqrt(2)\n    metric = (delta/sigma_clip)*sqrt + np.log(sigma_clip*sqrt)\n    return np.mean(metric)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm.notebook import tqdm\n\nmetric = []\nfor q in tqdm(range(1, 10)):\n    m = []\n    for p in vl_p:\n        x, y = [], []\n        tab = [] \n        \n        if p in ['ID00011637202177653955184', 'ID00052637202186188008618']:\n            continue\n            \n        img_set = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n        img_set = np.random.choice(img_set, size=20)\n        for i in img_set:\n            x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/{i}')) \n            y.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/{i}'))\n            tab.append(get_tab(train.loc[train.Patient == p, :])) \n        tab = np.array(tab) \n    \n        x = np.expand_dims(x, axis=-1)\n        y = np.expand_dims(y, axis=-1)\n        _a = model.predict([x,y, tab]) \n        a = np.quantile(_a, q / 10)\n        \n        percent_true = train.Percent.values[train.Patient == p]\n        fvc_true = train.FVC.values[train.Patient == p]\n        weeks_true = train.Weeks.values[train.Patient == p]\n        \n        fvc = a * (weeks_true - weeks_true[0]) + fvc_true[0]\n        percent = percent_true[0] - a * abs(weeks_true - weeks_true[0])\n        m.append(score(fvc_true, fvc, percent))\n    print(np.mean(m))\n    metric.append(np.mean(m))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"q = (np.argmin(metric) + 1)/ 10\nq","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \nsub.head() ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv') \ntest.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"A_test, B_test, P_test,W, FVC= {}, {}, {},{},{} \nSTD, WEEK = {}, {} \nfor p in test.Patient.unique():\n    x,y = [],[]\n    tab = [] \n    img_set = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/')\n    img_set = np.random.choice(img_set, size=20)\n    for i in img_set:\n        x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/{i}')) \n        y.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/{i}'))\n        tab.append(get_tab(test.loc[test.Patient == p, :])) \n    tab = np.array(tab) \n            \n    x = np.expand_dims(x, axis=-1) \n    y = np.expand_dims(y, axis=-1) \n    _a = model.predict([x,y, tab]) \n    a = np.quantile(_a, q)\n    A_test[p] = a\n    B_test[p] = test.FVC.values[test.Patient == p] - a*test.Weeks.values[test.Patient == p]\n    P_test[p] = test.Percent.values[test.Patient == p] \n    WEEK[p] = test.Weeks.values[test.Patient == p]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for k in sub.Patient_Week.values:\n    p, w = k.split('_')\n    w = int(w) \n    \n    fvc = A_test[p] * w + B_test[p]\n    sub.loc[sub.Patient_Week == k, 'FVC'] = fvc\n    sub.loc[sub.Patient_Week == k, 'Confidence'] = (\n        P_test[p] - A_test[p] * abs(WEEK[p] - w) \n) \nsub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}