{"cells":[{"metadata":{},"cell_type":"markdown","source":"
\n
\n

OSIC Pulmonary Fibrosis Progression

\n

Predict lung function decline

\n
"},{"metadata":{},"cell_type":"markdown","source":"#
1. Introduction
\n\n### 1.1 Pulmonary fibrosis[](http://)\n\n[The word \"**pulmonary**\" means lung and the word \"**fibrosis**\" means scar tissue— similar to scars that you may have on your skin from an old injury or surgery.](https://www.pulmonaryfibrosis.org/life-with-pf/about-pf) So, in its simplest sense, pulmonary fibrosis (PF) means scarring in the lungs. Over time, the scar tissue can destroy the normal lung and make it hard for oxygen to get into your blood. Low oxygen levels (and the stiff scar tissue itself) can cause you to feel short of breath, particularly when walking and exercising.\n\n\n\nImage Credits:- https://www.pulmonaryfibrosis.org/\n\n\n* **Prognosis** - Prognosis is a term for the predicted course of a disease. People commonly use the word to refer to an individual’s life expectancy, how long the person is likely to live. However, prognosis can also refer to the chance that a disease can be cured and the outlook for functional recovery, which includes the prospects of being able to return to work, engage in recreation, as well as the expected degree of help that will be necessary to accomplish activities of daily living.\n\n* **Forced vital capacity (FVC)** - FVC is the amount of air that can be forcibly exhaled from your lungs after taking the deepest breath possible, as measured by spirometry. FVC can also help doctors assess the progression of lung disease and evaluate the effectiveness of treatment.\n\nPatient's FVC volume can be compared with the standard FVC for similir age, sex, height, and weight. Patient FVC can also be compared with her/his own previous FVC values, if applicable, to determine whether pulmonary condition is progressing or if lung function is improving under treatment. FVC also may be expressed as a percentage of the predicted FVC.\n\nThe normal FVC range for an adult is between 3.0 and 5.0 L.\n\n\n### 1.2 Objective of Competition\n\nThe aim of this competition is to **predict a patient’s severity of decline in lung function** based on a CT scan of their lungs. Lung function is assessed based on output from a spirometer, which measures the **forced vital capacity (FVC)**, i.e. the volume of air exhaled. The challenge is to use machine learning techniques to make a prediction with the image, metadata, and baseline FVC as input.\n\n### 1.3 Evaluated of competition\n\nThis competition is evaluated on a modified version of the **Laplace Log Likelihood**. In medical applications, it is useful to evaluate a model's confidence in its decisions. Accordingly, the metric is designed to reflect both the accuracy and certainty of each prediction.\n\n\n\nImage Credits:- https://www.vosesoftware.com/riskwiki/Laplacedistribution.php\n\nFor each true FVC measurement, you will predict both an **FVC** and a **confidence measure** (standard deviation σ). The metric is computed as:\n\n\n**Confidence values smaller than 70 are clipped.**\n\n$$ \\large \\sigma_{clipped} = max(\\sigma, 70) $$\n\n\n**Errors greater than 1000 are also clipped in order to avoid large errors.**\n\n$$ \\large \\Delta = min ( |FVC_{true} - FVC_{predicted}|, 1000 ) $$\n\n\n**The metric is defined as:**\n\n$$ \\Large metric = - \\frac{\\sqrt{2} \\Delta}{\\sigma_{clipped}} - \\ln ( \\sqrt{2} \\sigma_{clipped} ) $$\n\n\nThe leaderboard is calculated with approximately 15% of the test data. The final results will be based on the other 85%, so the final standings may be different.\n\n"},{"metadata":{},"cell_type":"markdown","source":"## Contents\n\n* [Introduction](#Introduction)\n* [Importing libraries](#libraries)\n* [Load Data](#dataLoad)\n* [Exploratory Data Analysis](#EDA)\n* [Visualising Images : DECOM](#ImageVisuals)\n"},{"metadata":{},"cell_type":"markdown","source":"#
2. Importing libraries
"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Import Libraries\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport datetime\nimport random\nimport math\nimport matplotlib\nfrom termcolor import colored\nimport os\nfrom os import listdir\nfrom os.path import join, getsize\nimport glob\nimport cv2\n\n#plotly\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.offline import iplot\nimport cufflinks\ncufflinks.go_offline()\ncufflinks.set_config_file(world_readable=True, theme='pearl')\n\n# \nfrom skimage import measure\nfrom skimage.morphology import disk, opening, closing\n\nimport tensorflow as tf\nfrom tensorflow.keras import Model\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import train_test_split, KFold\n\nfrom tensorflow.keras.layers import (\n Dense, Dropout, Activation, Flatten, Input, GlobalAveragePooling2D, Add, Conv2D, AveragePooling2D, \n LeakyReLU, Concatenate \n)\n\n\n# Magic function to display In-Notebook display\n%matplotlib inline\n\n# Setting seabon style\nsns.set(style='darkgrid', palette='Set2')\n\n# Suppress warnings \nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Settings for pretty nice plots\nplt.style.use('fivethirtyeight')\nplt.show()\n\n# pydicom\nimport pydicom\n\n# Print versions of libraries\nprint(f\"Numpy version : Numpy {np.__version__}\")\nprint(f\"Pandas version : Pandas {pd.__version__}\")\nprint(f\"Matplotlib version : Matplotlib {matplotlib.__version__}\")\nprint(f\"Seaborn version : Seaborn {sns.__version__}\")\nprint(f\"Tensorflow version : Tensorflow {tf.__version__}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Install the EfficientNet Keras Library\n!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index\n!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import efficientnet.tfkeras as efn","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sets integer starting value "},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=100):\n random.seed(seed)\n os.environ['PYTHONHASHSEED'] = str(seed)\n np.random.seed(seed)\n tf.random.set_seed(seed)\n \nseed_everything(101)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#
2. Load Data
"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# List files available\nbase_dir = \"../input/osic-pulmonary-fibrosis-progression/\"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"list(os.listdir(base_dir))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Dataset Details \n* train.csv - the training set, contains full history of clinical information\n* test.csv - the test set, contains only the baseline measurement\n* train/ - contains the training patients' baseline CT scan in DICOM format\n* test/ - contains the test patients' baseline CT scan in DICOM format\n* sample_submission.csv - demonstrates the submission format"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train & Test set shape\ntrain_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv', encoding = 'latin-1')\ntest_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv', encoding = 'latin-1')\nsubmission_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv', encoding = 'latin-1')\n\nprint(colored('Training data set shape.......... : ','yellow'),train_df.shape)\nprint(colored('Test data set shape...............: ','red'),test_df.shape)\nprint(colored('Submission data set shape.........: ','blue'),submission_df.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# print top 5 rows of train set\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# print top 5 rows of test set\ntest_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Columns Details in train.csv and test.csv\n* **Patient** - a unique Id for each patient (also the name of the patient's DICOM folder)\n* **Weeks** - the relative number of weeks pre/post the baseline CT (may be negative)\n* **FVC** - the recorded lung capacity in ml\n* **Percent** - a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics\n* **Age** \n* **Sex** \n* **SmokingStatus** "},{"metadata":{},"cell_type":"markdown","source":"#
3. Exploratory Data Analysis
\n\nThe purpose of exploratory data analysis is to: Check for missing data and other mistakes. Gain maximum insight into the data set and its underlying structure. Uncover a parsimonious model, one which explains the data with a minimum number of predictor variables."},{"metadata":{},"cell_type":"markdown","source":"## Concise Summary"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Null values and Data types\nprint(colored('Train Set !!', 'yellow'))\nprint(colored('------------', 'yellow'))\nprint(train_df.info())\n\nprint('\\n')\n\nprint(colored('Test Set !!','red'))\nprint(colored('-----------','red'))\nprint(test_df.info())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is no missing values in train_df and test_df."},{"metadata":{},"cell_type":"markdown","source":"## Descriptive Statistics"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Null values and Data types\nprint(colored('Train Set !!', 'yellow'))\nprint(train_df.describe())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Patient Age ranges from 49 years to 88 years with 67 years average age and 7 years standard deviation.\n* The normal FVC range for an adult is between 3000ml to 5000ml. Dataset FVC ranges from 2690.48ml to 6399.00ml.\n* FVC are measured earliest at -5th week and latest by 133th week."},{"metadata":{},"cell_type":"markdown","source":"## Missing Values"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Total missing values for each feature\nprint(colored('Missing values in Train Set !!', 'yellow'))\nprint(train_df.isnull().sum())\n\nprint(\"\\n\")\n\nprint(colored('Missing values in Test Set !!', 'red'))\nprint(test_df.isnull().sum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.groupby( ['Sex','SmokingStatus'] )['FVC'].agg( ['mean','std','count'] ).sort_values(by=['Sex','count'],ascending=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is VERY weird: FVC and Percent are the highest for people that still smoke and the lowest for people that never smoked. HOWEVER, we need to keep in mind that the percentage of people that still smoke is very low. So, we CAN'T conclude that if a person smokes it's highly likely that will have a high FVC."},{"metadata":{},"cell_type":"markdown","source":"No missing value in either train/test datasets."},{"metadata":{},"cell_type":"markdown","source":"## Patients Counts"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Total number of Patient in the dataset(train+test)\n\nprint(colored(\"Total Patients in Train set... : \", 'yellow'),train_df['Patient'].count())\nprint(colored(\"Total Patients in Test set.... : \", 'red'),test_df['Patient'].count())\nprint(\"\\n\")\nprint(colored(\"Unique Patients in Train set...: \", 'yellow'),train_df['Patient'].nunique())\nprint(colored(\"Unique Patients in Test set....: \", 'red'),test_df['Patient'].nunique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"Few most repeated Patients in Train set: \", 'yellow'))\nprint(train_df['Patient'].value_counts().head())\n\nprint(\"\\n\")\n\nprint(colored(\"Few most repeated Patients in Test set: \", 'red'))\nprint(test_df['Patient'].value_counts().head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Unique patients\n\nLet's create a new data set having only unique patient details."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df_unique = train_df[['Patient', 'Age', 'Sex', 'SmokingStatus']].drop_duplicates().reset_index()\nprint(colored(\"Shape of unique patient data set : \",'yellow'),train_df_unique.shape)\ntrain_df_unique.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Frequency of a patient in Train set\n\nLet's count how many times a particular patient repeated in train set."},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_feq = train_df.groupby(['Patient'])['Patient'].count()\npatient_feq = pd.DataFrame({'Patient':patient_feq.index, 'Frequency':patient_feq.values})\n\n# Merge two dataframes based on patient's ids.\ntrain_df_unique = pd.merge(train_df_unique,patient_feq,how='inner',on='Patient')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df_unique.sort_values(by='Frequency', ascending=False).head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.bar(train_df_unique, x='Patient',y ='Frequency',color='Frequency')\nfig.update_layout(xaxis={'categoryorder':'total ascending'},title='Frequency of each patient')\nfig.update_xaxes(showticklabels=False)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Every patient is observed between 6 to 10 times however most of them have observed 9 times."},{"metadata":{},"cell_type":"markdown","source":"## Number of CT Scans for each patient in Train set\n\nwe are provided with a baseline chest CT scans at week intervats for each patients. So lets count how many CT Scan have each patient."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Creating unique patient lists \n# (here patient == dictory and files == CT Scan)\ntrain_dir = '../input/osic-pulmonary-fibrosis-progression/train/'\n\npatient_ids = os.listdir(train_dir)\npatient_ids = sorted(patient_ids)\n\n#Creating a new blank dataframe\nCtScan = pd.DataFrame(columns=['Patient','CtScanCount'])\n\n\nfor patient_id in patient_ids:\n # count number of images in each folder\n cnt = len(os.listdir(train_dir + patient_id))\n # insert patient id and ct scan count in dataframe\n CtScan.loc[len(CtScan)] = [patient_id,cnt]\n \n\n# Merge two dataframes based on patient's ids.\npatient_df = pd.merge(train_df_unique,CtScan,how='inner',on='Patient')\n\n# Reset index\n# https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.reset_index.html\npatient_df = patient_df.reset_index(drop=True)\n\n# Print new dataframe\npatient_df.head()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"CT Scans numbers in Train set \",\"yellow\"))\nprint(colored(\"Maximum number of CT Scans for a patient.... : \",\"blue\"),patient_df['CtScanCount'].max())\nprint(colored(\"Minimum number of CT Scans for a patient.... : \",\"blue\"),patient_df['CtScanCount'].min())\nprint(colored(\"Average number of CT Scans per patient...... : \",\"blue\"),round(patient_df['CtScanCount'].mean(),3))\nprint(colored(\"Total number of CT Scans of all patients.... : \",\"blue\"),patient_df['CtScanCount'].sum())\nprint(colored(\"Median of CT Scans counts................... : \",\"blue\"),patient_df['CtScanCount'].median())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Huge imbalance in the number of CT scans: half of the patients have less that 100 CT scans."},{"metadata":{},"cell_type":"markdown","source":"## Number of CT Scans for each patient in Test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Creating unique patient lists \n# (here patient == dictory and files == CT Scan)\ntest_dir = '../input/osic-pulmonary-fibrosis-progression/test/'\n\ntest_patient_ids = os.listdir(test_dir)\ntest_patient_ids = sorted(test_patient_ids)\n\n#Creating a new blank dataframe\nTestCtScan = pd.DataFrame(columns=['Patient','CtScanCount'])\n\nfor patient_id in test_patient_ids:\n # count number of images in each folder\n cnt = len(os.listdir(test_dir + patient_id))\n # insert patient id and ct scan count in dataframe\n TestCtScan.loc[len(TestCtScan)] = [patient_id,cnt]\n \n\n# Merge two dataframes based on patient's ids.\ntest_patient_df = pd.merge(test_df,TestCtScan,how='inner',on='Patient').reset_index()\n\n# Print new dataframe\ntest_patient_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"CT Scans numbers in Test set \",\"red\"))\nprint(colored(\"Maximum number of CT Scans for a patient... : \",\"green\"),test_patient_df['CtScanCount'].max())\nprint(colored(\"Minimum number of CT Scans for a patient... : \",\"green\"),test_patient_df['CtScanCount'].min())\nprint(colored(\"Average number of CT Scans per patient..... : \",\"green\"),test_patient_df['CtScanCount'].mean())\nprint(colored(\"Total number of CT Scans of all patients... : \",\"green\"),test_patient_df['CtScanCount'].sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Distribution of weeks"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['Weeks'].iplot(kind='hist',\n bins=100, xTitle='Weeks', yTitle='Frequency', \n linecolor='white',opacity=0.7,\n color='rgb(0, 200, 200)', theme='pearl',\n bargap=0.01, title='Distribution of Weeks')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of the patients CT scans done in between 4th to 20th week."},{"metadata":{},"cell_type":"markdown","source":"## Distribution of Patients age"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_df['Age'].iplot(kind='hist',\n bins=10, xTitle='Age', yTitle='Frequency', \n linecolor='white',opacity=0.7,\n color='rgb(0, 100, 200)', theme='pearl',\n bargap=0.01, title='Distribution of Age column')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Range of patients age is between 48-88 years where we have more records for patients in the age range 64-74 years."},{"metadata":{},"cell_type":"markdown","source":"### Distribution of Patient gender"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"Gender wise distribution of patients :\",\"blue\"))\nprint(patient_df['Sex'].value_counts())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sex_count = patient_df[\"Sex\"].value_counts()\nsex_labels = patient_df[\"Sex\"].unique()\n\nfig = px.pie(patient_df, values=sex_count, names=sex_labels, hover_name=sex_labels)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"More number of male patients than female patients."},{"metadata":{},"cell_type":"markdown","source":"### Distribution of Age vs Gender"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 6))\n\nsns.kdeplot(patient_df[patient_df['Sex'] == 'Male']['Age'], label = 'Male',shade=True)\nsns.kdeplot(patient_df[patient_df['Sex'] == 'Female']['Age'], label = 'Female',shade=True)\n\nplt.xlabel('Age (years)'); \nplt.ylabel('Density'); \nplt.title('Distribution of Ages');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Male and female records are almost distributed throughout the age range."},{"metadata":{},"cell_type":"markdown","source":"### Distribution of 'SmokingStatus' feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored('Total Smoking counts', 'red'))\nprint(patient_df['SmokingStatus'].value_counts())\n\nprint(\"\\n\")\nprint(colored(\"Male Smoking counts\",'blue'))\nprint(patient_df[patient_df['Sex']=='Male']['SmokingStatus'].value_counts())\n\nprint(\"\\n\")\nprint(colored(\"Female Smoking counts\",'green'))\nprint(patient_df[patient_df['Sex']=='Female']['SmokingStatus'].value_counts())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Distribution of Age vs SmokingStatus"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 6))\n\nsns.kdeplot(patient_df.loc[patient_df['SmokingStatus'] == 'Ex-smoker', 'Age'], label = 'Ex-smoker',shade=True)\nsns.kdeplot(patient_df.loc[patient_df['SmokingStatus'] == 'Never smoked', 'Age'], label = 'Never smoked',shade=True)\nsns.kdeplot(patient_df.loc[patient_df['SmokingStatus'] == 'Currently smokes', 'Age'], label = 'Currently smokes', shade=True)\n\n# Labeling of plot\nplt.xlabel('Age (years)'); \nplt.ylabel('Density'); \nplt.title('Distribution of Ages');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Gender wise smoking distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,8))\nsns.countplot(x='SmokingStatus', data=patient_df, hue='Sex')\nplt.title('Gender split by SmokingStatus', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Records with patient who have never smoked have almost equal distribution of male and female patients whereas majority of ex-smokers are males."},{"metadata":{},"cell_type":"markdown","source":"## FVC - The forced vital capacity\n\nLung function is assessed based on output from a spirometer, which measures the **forced vital capacity (FVC)**, i.e. the volume of air exhaled. FVC can also help doctors assess the progression of lung disease and evaluate the effectiveness of treatment.\n\nA person who has Diagnose obstructive lung diseases such as asthma and chronic obstructive pulmonary disease (COPD) has a lower FVC result than a healthy person.Decreases in the FVC value may mean the lung disease is getting worse.\n\n* Average normal values in healthy males aged 20-60 range from 5.5 to 4.75 liters.\n* Average normal values for females aged 20-60 range from 3.75 to 3.25 liters. \n* Percent- a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics. Percentage with normal test values falling between 80% and 120% of the average values.\n\nReferecne : https://www.nuvoair.com/blog/do-you-know-how-to-interpret-the-results-of-your-spirometry-test"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"Maximum value of FVC... :\",'blue'),colored(train_df['FVC'].max(),'blue'))\nprint(colored(\"Minimum value of FVC... :\",'green'),colored(train_df['FVC'].min(),'green'))\n\nprint(\"\\n\")\n\n# Distribution of FVC\nprint(colored(\"Distribution of FVC\",\"yellow\"))\nprint(colored(train_df['FVC'].value_counts(normalize=False, ascending=False, bins=62).head(),\"yellow\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC Distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['FVC'].iplot(kind='hist',\n xTitle='Lung Capacity(ml)', \n yTitle='Frequency', \n linecolor='black', \n bargap=0.2,\n title='Distribution of the FVC in the training set')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Smoking Status"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.violin(train_df, y='FVC', x='SmokingStatus', \n box=True, color='Sex', points=\"all\", hover_data=train_df.columns, title=\"FVC of various Smoking Status\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Age"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.scatter(train_df, x=\"Age\", y=\"FVC\", color='Sex', title='FVC values for Patient Age')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Males have higher FVC than females irrespective of age."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[train_df['FVC'] > 5000].sort_values(by='FVC', ascending=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Week"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = go.Figure()\nfig = px.scatter(train_df, x=\"Weeks\", y=\"FVC\", color='SmokingStatus')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of the FVC test is done in between 0 to 20 weeks. Also Ex-smoker males have more FVC than others."},{"metadata":{},"cell_type":"markdown","source":"### FVC of oldest and youngest person"},{"metadata":{"trusted":true},"cell_type":"code","source":"# patient = train_df[train_df['FVC'] == train_df['FVC'].max()]\npatient = train_df[(train_df['Age'] == train_df['Age'].max()) | (train_df['Age'] == train_df['Age'].min())]\nfig = px.line(patient, x=\"Weeks\", y=\"FVC\", color='Age',line_group=\"Sex\", hover_name=\"SmokingStatus\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Aging is associated with progressive decline in lung function as shown in predicted by above plot that FVC of older person is less than younger person."},{"metadata":{},"cell_type":"markdown","source":"## Percent\nPercent approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics. Percentage with normal test values falling between 80% and 120% of the average values."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored(\"Maximum value of Percent... :\",'blue'),colored(train_df['Percent'].max(),'blue'))\nprint(colored(\"Minimum value of Percent... :\",'green'),colored(train_df['Percent'].min(),'green'))\n\nprint(\"\\n\")\n\n# Distribution of Percent\nprint(colored(\"Distribution of Percent\",\"yellow\"))\nprint(colored(train_df['Percent'].value_counts(normalize=False, ascending=False, bins=62).head(),\"yellow\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Percent Distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['Percent'].iplot(kind='hist',\n xTitle='Percent', \n yTitle='Frequency', \n linecolor='black', \n bargap=0.2,\n title='Distribution of Percent in the training set')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Percent vs SmokingStatus"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.violin(train_df, y='Percent', x='SmokingStatus', \n box=True, color='Sex', points=\"all\", hover_data=train_df.columns, title=\"Percent of various Smoking Status\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Percent of oldest and youngest person"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient = train_df[(train_df['Age'] == train_df['Age'].max()) | (train_df['Age'] == train_df['Age'].min())]\nfig = px.line(patient, x=\"Weeks\", y=\"Percent\", color='Age',line_group=\"Sex\", hover_name=\"SmokingStatus\")\n\npatient = train_df[(train_df['Age'] == train_df['Age'].max()) | (train_df['Age'] == train_df['Age'].min())]\nfig = px.line(patient, x=\"Weeks\", y=\"Percent\", color='Age',line_group=\"Sex\", hover_name=\"SmokingStatus\")\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Percent vs Age "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.scatter(train_df, x=\"Age\", y=\"Percent\", color=\"SmokingStatus\", marginal_y=\"violin\",\n marginal_x=\"box\", trendline=\"ols\", template=\"simple_white\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC vs Percent"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.scatter(train_df, x=\"FVC\", y=\"Percent\", color='SmokingStatus', size='Age', \n hover_name='SmokingStatus',hover_data=['Weeks'])\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Patient Overlap"},{"metadata":{},"cell_type":"markdown","source":"## Correlation among varous features"},{"metadata":{"trusted":true},"cell_type":"code","source":"corrmat = train_df.corr() \nfig = px.imshow(corrmat, x=corrmat.columns, y=corrmat.columns)\nfig.update_xaxes(side=\"top\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* There is high correlation between FVC and Percent: when the volume of air increases, the Percent increases as well.\n\n* There is no correlation between FVC/Percent and Age, meaning that Age has no influence on the volume of exhaled air."},{"metadata":{},"cell_type":"markdown","source":"#
4. Visualising DICOM Images
"},{"metadata":{},"cell_type":"markdown","source":"### Digital Imaging and COmmunications in Medicine - DICOM\n\nDICOM(Digital Imaging and COmmunications in Medicine) is the de-facto standard that establishes rules that allow medical images(X-Ray, MRI, CT) and associated information to be exchanged between imaging equipment from different vendors, computers, and hospitals.\n\nDICOM files typically have a .dcm extension and provides a means of storing data in separate 'tags' such as patient information as well as image/pixel data. A DICOM file consists of a header and image data sets packed into a single file. The information within the header is organized as a constant and standardized series of tags.\n\nBy extracting data from these tags one can access important information regarding the patient demographics, study parameters, etc.\n\n### Pydicom\n\nPydicom is a python package for parsing DICOM files and makes it easy to covert DICOM files into pythonic structures for easier manipulation. Files are opened using pydicom.dcmread"},{"metadata":{},"cell_type":"markdown","source":"## Patients & their CT Scans in Training Images Folder"},{"metadata":{"trusted":true},"cell_type":"code","source":"## Patients & their CT Scans in Training Images Folder\n\nfile_len = folder_len = 0\nfiles = []\n\nfor dirpath, dirnames, filenames in os.walk(train_dir):\n file_len += len(filenames)\n folder_len += len(dirnames)\n files.append(len(filenames))\n\nprint(\"Training folder contains\", f'{file_len:,}', \"CT scans for all patients.\") \nprint('Training folder have only',f'{folder_len:,}', \"unique patients.\")\n\nprint(\"\\n\")\n\nprint('Each patient have', f'{round(np.mean(files)):,}', 'average number of CT scans.')\nprint('Maximum images per patient', f'{round(np.max(files)):,}')\nprint('Minimum images per patient', f'{round(np.min(files)):,}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Extracting DIOCOM files information in a dataframe"},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/schlerp/getting-to-know-dicom-and-the-data\n\ndef show_dcm_info(file_path):\n #print(colored(\"Filename.........:\",'yellow'),file_path)\n #print()\n print(colored(\"File Path...........:\",'blue'), file_path)\n \n dataset = pydicom.dcmread(file_path)\n\n pat_name = dataset.PatientName\n display_name = pat_name.family_name + \", \" + pat_name.given_name\n \n print(colored(\"Patient's name......:\",'blue'), display_name)\n print(colored(\"Patient id..........:\",'blue'), dataset.PatientID)\n print(colored(\"Patient's Sex.......:\",'blue'), dataset.PatientSex)\n print(colored(\"Modality............:\",'blue'), dataset.Modality)\n print(colored(\"Body Part Examined..:\",'blue'), dataset.BodyPartExamined)\n \n if 'PixelData' in dataset:\n rows = int(dataset.Rows)\n cols = int(dataset.Columns)\n print(colored(\"Image size..........:\",'blue'),\" {rows:d} x {cols:d}, {size:d} bytes\".format(\n rows=rows, cols=cols, size=len(dataset.PixelData)))\n if 'PixelSpacing' in dataset:\n print(colored(\"Pixel spacing.......:\",'blue'),dataset.PixelSpacing)\n dataset.PixelSpacing = [1, 1]\n plt.figure(figsize=(10, 10))\n plt.imshow(dataset.pixel_array, cmap='gray')\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for file_path in glob.glob(train_dir + '*/*.dcm'):\n show_dcm_info(file_path)\n break # Comment this out to see all","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_dcm_info(train_dir + 'ID00027637202179689871102/11.dcm')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_dir = train_dir + \"ID00123637202217151272140\"\n\nprint(\"total images for patient ID00123637202217151272140: \", len(os.listdir(patient_dir)))\n\n# view first (columns*rows) images in order\nfig=plt.figure(figsize=(16, 16))\ncolumns = 4\nrows = 5\nimglist = os.listdir(patient_dir)\nfor i in range(1, columns*rows +1):\n filename = patient_dir + \"/\" + str(i) + \".dcm\"\n ds = pydicom.dcmread(filename)\n fig.add_subplot(rows, columns, i)\n plt.imshow(ds.pixel_array, cmap='gray')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# view first (columns*rows) images in order\nfig=plt.figure(figsize=(16, 16))\ncolumns = 4\nrows = 5\nimglist = os.listdir(patient_dir)\nfor i in range(1, columns*rows +1):\n filename = patient_dir + \"/\" + str(i) + \".dcm\"\n ds = pydicom.dcmread(filename)\n fig.add_subplot(rows, columns, i)\n plt.imshow(ds.pixel_array, cmap='jet')\n #plt.imshow(cv2.cvtColor(ds.pixel_array, cv2.COLOR_BGR2RGB))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Loading DICOM files\n\nDicom files contain a lot of metadata (such as the pixel size, so how long one pixel is in every dimension in the real world).\n\nThis pixel size/coarseness of the scan differs from scan to scan (e.g. the distance between slices may differ), which can hurt performance of CNN approaches. \n\nBelow is code to load a scan, which consists of multiple slices, which we simply save in a Python list. Every folder in the dataset is one scan (so one patient). One metadata field is missing, the pixel size in the Z direction, which is the slice thickness. Fortunately we can infer this, and we add this to the metadata."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Ref : \n# https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial\n# https://www.kaggle.com/akh64bit/full-preprocessing-tutorial\n# https://www.researchgate.net/post/How_can_I_convert_pixel_intensity_values_to_housefield_CT_number\n\n# Load the scans in given folder path\ndef load_scan(path):\n slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n try:\n slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n except:\n slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n \n for s in slices:\n s.SliceThickness = slice_thickness\n \n return slices","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The unit of measurement in CT scans is the Hounsfield Unit (HU), which is a measure of radiodensity. CT scanners are carefully calibrated to accurately measure this. From Wikipedia:\n\n\n\nBy default however, the returned values are not in this unit. Let's fix this.\n\nSome scanners have cylindrical scanning bounds, but the output image is square. The pixels that fall outside of these bounds get the fixed value -2000. The first step is setting these values to 0, which currently corresponds to air. Next, let's go back to HU units, by multiplying with the rescale slope and adding the intercept (which are conveniently stored in the metadata of the scans!)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(slices):\n image = np.stack([s.pixel_array for s in slices])\n # Convert to int16 (from sometimes int16), \n # should be possible as values should always be low enough (<32k)\n image = image.astype(np.int16)\n\n # Set outside-of-scan pixels to 0\n # The intercept is usually -1024, so air is approximately 0\n image[image == -2000] = 0\n \n # Convert to Hounsfield units (HU)\n for slice_number in range(len(slices)):\n \n intercept = slices[slice_number].RescaleIntercept\n slope = slices[slice_number].RescaleSlope\n \n if slope != 1:\n image[slice_number] = slope * image[slice_number].astype(np.float64)\n image[slice_number] = image[slice_number].astype(np.int16)\n \n image[slice_number] += np.int16(intercept)\n \n return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at one of the patients."},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient = load_scan(train_dir + patient_ids[0])\nfirst_patient_pixels = get_pixels_hu(first_patient)\n\nplt.figure(figsize=(10, 10))\nplt.hist(first_patient_pixels.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\n# Show some slice in the middle\nplt.figure(figsize=(10, 10))\nplt.imshow(first_patient_pixels[15], cmap=plt.cm.gray)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking at the table from Wikipedia and this histogram, we can clearly see which pixels are air and which are tissue. We will use this for lung segmentation in a bit "},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at the first dicom file of a patient:"},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient_scan = load_scan(train_dir + patient_ids[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient_scan[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualization using gif"},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_lungwin(img, hu=[-1200., 600.]):\n lungwin = np.array(hu)\n newimg = (img-lungwin[0]) / (lungwin[1]-lungwin[0])\n newimg[newimg < 0] = 0\n newimg[newimg > 1] = 1\n newimg = (newimg * 255).astype('uint8')\n return newimg","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient_scan_array = set_lungwin(get_pixels_hu(first_patient_scan))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import imageio\nfrom IPython.display import Image\n\nimageio.mimsave(\"/tmp/gif.gif\", first_patient_scan_array, duration=0.00001)\nImage(filename=\"/tmp/gif.gif\", format='png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Transforming to Hounsfield Units \nBefore starting, let's plot the pixelarray distribution of some dicom files to get an impression of the raw data:\n\nRef : https://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nfor n in range(10):\n image = first_patient_scan[n].pixel_array.flatten()\n rescaled_image = image * first_patient_scan[n].RescaleSlope + first_patient_scan[n].RescaleIntercept\n sns.distplot(image.flatten(), ax=ax[0]);\n sns.distplot(rescaled_image.flatten(), ax=ax[1])\nax[0].set_title(\"Raw pixel array distributions for 10 examples\")\nax[1].set_title(\"HU unit distributions for 10 examples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are some raw values at -2000. They correspond to images with a circular boundary within the image. The \"outside\" of this circle value is often set to -2000 (or in other competitions I found also -3000) by default."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,4,figsize=(20,3))\nax[0].set_title(\"Original CT-scan\")\nax[0].imshow(first_patient_scan[0].pixel_array, cmap=\"bone\")\nax[1].set_title(\"Pixelarray distribution\");\nsns.distplot(first_patient_scan[0].pixel_array.flatten(), ax=ax[1]);\n\nax[2].set_title(\"CT-scan in HU\")\nax[2].imshow(first_patient_pixels[0], cmap=\"bone\")\nax[3].set_title(\"HU values distribution\");\nsns.distplot(first_patient_pixels[0].flatten(), ax=ax[3]);\n\nfor m in [0,2]:\n ax[m].grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The scan of our example patient had a circular boundary and now all raw values per slice are scaled to H-units."},{"metadata":{},"cell_type":"markdown","source":"## Tissue segmentation \n\nA scan may have a pixel spacing of [2.5, 0.5, 0.5], which means that the distance between slices is 2.5 millimeters. For a different scan this may be [1.5, 0.725, 0.725], this can be problematic for automatic analysis (e.g. using ConvNets)!\n\nIn order to reduce the problem space, we can segment the lungs (and usually some tissue around it). \n\nIt involves quite a few smart steps.\n\nThreshold the image (-320 HU is a good threshold, but it doesn't matter much for this approach).\nDo connected components, determine label of air around person, fill this with 1s in the binary image\nOptionally: For every axial slice in the scan, determine the largest solid connected component (the body+air around the person), and set others to 0. This fills the structures in the lungs in the mask.\nKeep only the largest air pocket (the human body has other pockets of air here and there).\n\nWith -320 we are separating between lungs (-700) /air (-1000) and tissue with values close to water (0)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_lung_mask(image):\n segmented = np.zeros(image.shape) \n \n for n in range(image.shape[0]):\n binary_image = np.array(image[n] > -320, dtype=np.int8)+1\n labels = measure.label(binary_image)\n \n background_label_1 = labels[0,0]\n background_label_2 = labels[0,-1]\n background_label_3 = labels[-1,0]\n background_label_4 = labels[-1,-1]\n \n #Fill the air around the person\n binary_image[background_label_1 == labels] = 2\n binary_image[background_label_2 == labels] = 2\n binary_image[background_label_3 == labels] = 2\n binary_image[background_label_4 == labels] = 2\n \n #We have a lot of remaining small signals outside of the lungs that need to be removed. \n #In our competition closing is superior to fill_lungs \n selem = disk(4)\n binary_image = closing(binary_image, selem)\n \n binary_image -= 1 #Make the image actual binary\n binary_image = 1-binary_image # Invert it, lungs are now 1\n \n segmented[n] = binary_image.copy() * image[n]\n \n return segmented","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented = segment_lung_mask(np.array([first_patient_pixels[20]]))\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(first_patient_pixels[20], cmap=\"Blues_r\")\nax[1].imshow(segmented[0], cmap=\"Blues_r\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_lungs = segment_lung_mask(first_patient_pixels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_lungs.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(6,5, figsize=(20,20))\nfor n in range(6):\n for m in range(5):\n ax[n,m].imshow(segmented_lungs[n*5+m], cmap=\"Blues_r\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Reset the index"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.reset_index(inplace = True , drop = True)\npatient_df.reset_index(inplace = True , drop = True)\ntest_df.reset_index(inplace = True , drop = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model Building"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv', encoding = 'latin-1')\n# test_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv', encoding = 'latin-1')\n# submission_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv', encoding = 'latin-1')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_sub = submission_df[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()\nprint(img_sub.sample(5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Dropout_model = 0.38559\nFVC_weight = 0.2\nConfidence_weight = 0.15\n\nconfig = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_efficientnet(model, shape):\n models_dict = {\n 'b0': efn.EfficientNetB0(input_shape=shape,weights=None,include_top=False),\n 'b1': efn.EfficientNetB1(input_shape=shape,weights=None,include_top=False),\n 'b2': efn.EfficientNetB2(input_shape=shape,weights=None,include_top=False),\n 'b3': efn.EfficientNetB3(input_shape=shape,weights=None,include_top=False),\n 'b4': efn.EfficientNetB4(input_shape=shape,weights=None,include_top=False),\n 'b5': efn.EfficientNetB5(input_shape=shape,weights=None,include_top=False),\n 'b6': efn.EfficientNetB6(input_shape=shape,weights=None,include_top=False),\n 'b7': efn.EfficientNetB7(input_shape=shape,weights=None,include_top=False)\n }\n return models_dict[model]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_model(shape=(512, 512, 1), model_class=None):\n inp = Input(shape=shape)\n base = get_efficientnet(model_class, shape)\n x = base(inp)\n x = GlobalAveragePooling2D()(x)\n inp2 = Input(shape=(5,)) # add the feature of MinFVC\n x2 = tf.keras.layers.GaussianNoise(0.2)(inp2)\n x = Concatenate()([x, x2]) \n x = Dropout(Dropout_model)(x)\n x = Dense(1)(x)\n model = Model([inp, inp2] , x)\n \n# weights = [w for w in os.listdir('../input/osic-model-weights') if model_class in w][0]\n# model.load_weights('../input/osic-model-weights/' + weights)\n return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_classes = ['b5'] #['b0','b1','b2','b3',b4','b5','b6','b7']\nmodels = [build_model(shape=(512, 512, 1), model_class=m) for m in model_classes]\nprint('Number of models: ' + str(len(models)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = models[0]\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE=128","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df['Patient'] = submission_df['Patient_Week'].apply(lambda x:x.split('_')[0])\nsubmission_df['Weeks'] = submission_df['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsubmission_df = submission_df[['Patient','Weeks','Confidence','Patient_Week']]\nsubmission_df = submission_df.merge(test_df.drop('Weeks', axis=1), on=\"Patient\")\n\nsubmission_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Merge Train, Test and Submission CSV files\ntrain_df['WHERE'] = 'train'\ntest_df['WHERE'] = 'val'\nsubmission_df['WHERE'] = 'test'\n\n# Merge train, test and submission dataset\ndata = train_df.append([test_df, submission_df])\n\nprint(colored('Train data set shape.......: ','yellow'), train_df.shape)\nprint(colored('Test data set shape........: ','green'), test_df.shape)\nprint(colored('Submission data set shape..: ','blue'), submission_df.shape)\nprint(colored('Comibined data set shape...: ','red'), data.shape)\n\nprint(\"\\n\")\n\nprint(colored('Unique Patient in Train data set shape.......: ','yellow'), train_df.Patient.nunique())\nprint(colored('Unique Patient in Test data set shape........: ','green'), test_df.Patient.nunique())\nprint(colored('Unique Patient in Submission data set shape..: ','blue'), submission_df.Patient.nunique())\nprint(colored('Unique Patient in Comibined data set shape...: ','red'), data.Patient.nunique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data.head(3).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\n\ndel base","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"COLS = ['Sex','SmokingStatus'] #,'Age'\nfeatures = []\n\nfor col in COLS:\n for mod in data[col].unique():\n features.append(mod)\n data[mod] = (data[col] == mod).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data.head(3).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#\ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\ndata['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min() )\ndata['week'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\ndata['percent'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )\ndata['FVC_Percent'] = data['FVC'] / data['Percent']\n\nfeatures += ['age','percent','week','BASE']\nprint(features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = data.loc[data.WHERE=='train']\ntest_df = data.loc[data.WHERE=='val']\nsubmission_df = data.loc[data.WHERE=='test']\n\ndel data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n\ndef score(y_true, y_pred):\n tf.dtypes.cast(y_true, tf.float32)\n tf.dtypes.cast(y_pred, tf.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 \n return K.mean(metric)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def qloss(y_true, y_pred):\n # Pinball loss for multiple quantiles\n qs = [0.2, 0.50, 0.8]\n q = tf.constant(np.array([qs]), dtype=tf.float32)\n e = y_true - y_pred\n v = tf.maximum(q*e, (q-1)*e)\n \n return K.mean(v)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_model(no_feature):\n z = L.Input((no_feature,), name=\"Patient\")\n x = L.Dense(100, activation=\"relu\", name=\"d1\")(z)\n x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), name=\"preds\")([p1, p2])\n \n model = M.Model(z, preds, name=\"CNN\")\n model.compile(loss=mloss(0.65), optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n \n return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(colored('Features for model building : ','yellow'),features)\n\ny = train_df['FVC'].values # train target\nX = train_df[features].values # fetures (1535, 9)\nze = submission_df[features].values # fetures of submission (730, 9) e: estimate\n\nprint(colored('Training data set shape for model building : ','yellow'),X.shape)\nprint(colored('Shape of submission.csv: ','yellow'),ze.shape) \n\nnh = X.shape[1] \nprint(colored('Number of features : ','yellow'),nh) # feature numbers (9,)\n\npe = np.zeros((ze.shape[0], 3)) #estimate of prediction\npred = np.zeros((X.shape[0], 3)) # prediction of truth ground","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = make_model(nh)\nprint(net.summary())\nprint(net.count_params())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"NFOLD = 5 # originally 5\nkf = KFold(n_splits=NFOLD)\n\n# %%time\ncnt = 0\nEPOCHS = 855\n\n\nfor tr_idx, val_idx in kf.split(X):\n cnt += 1\n print(\"\\n\")\n print(colored('Fold........... : ','red'),cnt)\n \n net = make_model(nh)\n net.fit(X[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, validation_data=(X[val_idx], y[val_idx]), verbose=0) #\n \n print(colored('Train.......... : ','yellow'),\"train\", net.evaluate(X[tr_idx], y[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n print(colored('Validation............ : ','yellow'),\"val\", net.evaluate(X[val_idx], y[val_idx], verbose=0, batch_size=BATCH_SIZE))\n \n print(colored('Predict Validation.. : ','yellow'))\n pred[val_idx] = net.predict(X[val_idx], batch_size=BATCH_SIZE, verbose=0)\n \n print(colored('Predict Test... : ','yellow'))\n pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\n\nprint(sigma_opt, sigma_mean)\nprint(unc.min(), unc.mean(), unc.max(), (unc>=0).mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.figure(figsize=(12,8))\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,6))\nplt.hist(unc)\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PREDICTION\nsubmission_df['FVC1'] = 1.*pe[:, 1]\nsubmission_df['Confidence1'] = pe[:, 2] - pe[:, 0]\nsubm = submission_df[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()\nsubm.loc[~subm.FVC1.isnull()].head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']\nif sigma_mean<70:\n subm['Confidence'] = sigma_opt\nelse:\n subm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"otest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nfor i in range(len(otest)):\n subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'FVC'] = otest.FVC[i]\n subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'Confidence'] = 0.1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_regression.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_sub = subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df1 = img_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)\ndf2 = reg_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df1[['Patient_Week']].copy()\ndf['FVC'] = FVC_weight*df1['FVC'] + (1-FVC_weight)*df2['FVC']\ndf['Confidence'] = Confidence_weight*df1['Confidence'] + (1-Confidence_weight)*df2['Confidence']\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Save submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Reference:\n* https://err.ersjournals.com/content/23/132/215\n* https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165727\n* https://www.kaggle.com/piantic/osic-pulmonary-fibrosis-progression-basic-eda\n* https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6738634/\n* https://www.semanticscholar.org/paper/Honeycomb-lung%3A-history-and-current-concepts.-Arakawa-Honma/9ea8579ddf8de97e308500ad73e680ba9b2c455d/figure/3\n* https://err.ersjournals.com/content/23/132/215\n* https://link.springer.com/article/10.1186/s12890-020-1061-x\n* https://www.pulmonologyadvisor.com/home/topics/restrictive-lung-disease/ct-honeycombing-in-interstitial-lung-disease-linked-to-higher-mortality-rates/\n* https://www.kaggle.com/thebigd8ta/osic-ensemble-iv/"},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}