{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#  OSIC Pulmonary fibrosis progression \n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"<img src=\"https://media.giphy.com/media/WtUK5I9TbWiRcGrVZh/giphy.gif\">","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 1. Description ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Suppose you are diagonsed with pulmonary fibrosis ( Its a disorder  with no known cause  and no known cure created by scarring of lungs ) , its outcome can range from long term stability to rapid  deterioration and doctors aren’t easily able to tell where an individual may fall on that spectrum. So its our job to take the responsibilty as we are the data scientist hence we wont let anyone suffer from that disease anymore  ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 2 . So what is Pulmonary Fibrosis ?","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Pulmonary fibrosis is a lung disease that occurs when lung tissue becomes damaged and scarred. This thickened, stiff tissue makes it more difficult for your lungs to work properly. As pulmonary fibrosis worsens, you become progressively more short of breath.\n\n## The scarring associated with pulmonary fibrosis can be caused by a multitude of factors. But in most cases, doctors can't pinpoint what's causing the problem. When a cause can't be found, the condition is termed idiopathic pulmonary fibrosis.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"<img src = 'https://www.wikidoc.org/images/d/d3/Pulmonary_fibrosis.gif' > ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 3.  So what do we need to predict ?","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## We need to predict a patient’s severity of decline in lung function based on a CT scan of their  lungs ,  we need to determine the  lung function based on output from a spirometer, which measures the volume of air inhaled and exhaled. The challenge is to use machine learning techniques to make a prediction with the image, metadata, and baseline FVC as input.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":" # 4 . so what is FVC ?","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## It stands for forced vital capacity . It  is the amount of air that can be forcibly exhaled from your lungs after taking the deepest breath possible, as measured by spirometry.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 5. Spirometry ?? what is that !","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Spirometry is the most common of the pulmonary function tests. It measures lung function, specifically the amount and/or speed of air that can be inhaled and exhaled. Spirometry is helpful in assessing breathing patterns that identify conditions such as asthma, pulmonary fibrosis, cystic fibrosis, and COPD.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 1 :  Importing the libaraies ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf \nfrom tensorflow.keras.layers import *\nfrom tensorflow.keras.models import Sequential\nimport tensorflow.io as tfio\nfrom keras.preprocessing import image\nimport matplotlib.pyplot as plt \nimport glob as glob \nimport seaborn as sns \nimport plotly.express as px\nimport plotly.graph_objects as go\nimport pandas as pd \nimport numpy as np\nfrom skimage import morphology , segmentation , measure \nfrom sklearn.preprocessing import OneHotEncoder , LabelEncoder \nfrom sklearn.compose import ColumnTransformer\nimport os\nimport pydicom\n!pip install dicom\nimport dicom \nimport imageio\nfrom IPython.display import Image\nfrom timeit import timeit\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as Layers\nimport tensorflow.keras.models as Models\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom sklearn.model_selection import KFold, GroupKFold, StratifiedKFold\nfrom sklearn.metrics import mean_absolute_error","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2 :  Now its time for Exploratory Data Analysis","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# checking the no of rows and columns train data\ntrain_x = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\nprint('the no of rows is {} and the no of columns is {} '.format(train_x.shape[0] , train_x.shape[1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_x.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# checking the no of rows and columns in test data \ntest_x = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nprint('the no of rows is {} and the no of columns is {} '.format(test_x.shape[0] , test_x.shape[1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_x.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets check the no of males and females \nsns.countplot( x = 'Sex' , data = train_x )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## as we can clearly see that no of males are way higher then no of females hence it is a point to remember ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Distribution of images for each patient ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Thanks to PAB97 for the plot, couldn't figure out how to do it with Seaborn .\n\n\nnew_df = train_x.groupby(\n    [\n        train_x.Patient,\n        train_x.Age,train_x.Sex, \n        train_x.SmokingStatus\n    ]\n)['Patient'].count()\n\nnew_df.index = new_df.index.set_names(\n    [\n        'id',\n        'Age',\n        'Sex',\n        'SmokingStatus'\n    ]\n)\n\nnew_df = new_df.reset_index()\nnew_df.rename(columns = {'Patient': 'freq'},inplace = True)\n\nfig = px.bar(new_df, x='id',y ='freq',color='freq')\nfig.update_layout(\n    xaxis={'categoryorder':'total ascending'},\n    title='Distribution of images for each patient'\n)\nfig.update_xaxes(showticklabels=False)\nfig.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# distribution of age ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.histogram(\n    new_df, \n    x='Age',\n    nbins = 42\n)\n\nfig.update_traces(\n    marker_color='rgb(158,202,225)', \n    marker_line_color='rgb(8,48,107)',\n    marker_line_width=1.5, \n    opacity=0.6\n)\n\nfig.update_layout(\n    title = 'Distribution of Age'\n)\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets check the smoking status \nsns.countplot( x = 'SmokingStatus' , data = train_x )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# here we can see that no of ex smoker is way higher ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":" fig = px.histogram(\n    train_x, \n    x='Age',\n    color='SmokingStatus',\n    color_discrete_map=\n        {\n            'Never smoked':'yellow',\n            'Currently smokes':'cyan',\n            'Ex-smoker': 'green', \n        },\n    hover_data=train_x.columns\n)\n\nfig.update_layout(title='Distribution of Age w.r.t. SmokingStatus for unique patients')\n\nfig.update_traces(\n    marker_line_color='black',\n    marker_line_width=1.5, \n    opacity=0.85\n)\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# now analyzing the sex with respect to smoking status \nplt.figure(figsize = (5 , 5))\nsns.countplot(x = 'Sex' , hue = 'SmokingStatus' , data = train_x)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# as we can clearly see that male is dominating in case of ex smoker ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.histogram(\n    train_x, \n    x='Age',\n    color='Sex',\n    color_discrete_map=\n        {\n            'Male':'blue',\n            'Female':'mediumturquoise'\n        },\n    hover_data=train_x.columns\n)\n\nfig.update_layout(title='Distribution of Age w.r.t. sex for unique patients')\n\nfig.update_traces(\n    marker_line_color='black',\n    marker_line_width=1.5, \n    opacity=0.85\n)\n\nfig.show()\n50\n55\n60\n65\n70\n75\n80\n85\n0\n20\n40\n60\n80\n100\n120\n140\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# as we can clearly see that the no of males are too high in between age ( 64 - 74 )","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# now lets see the correlation between features using heatmap \nsns.heatmap(train_x.corr() , annot = True , cmap=plt.cm.cool)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# as we can see that the percent and FVC are having a good relationship","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets check the FVC distribution graph\na= sns.distplot(train_x['FVC'] , color = 'r' , )\na.set_title('Distribution plot of SVC ' , color = 'g'  , fontsize = 18)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# now lets check the percent distribution graph \nb = sns.distplot(train_x['Percent'] , color = 'g')\nb.set_title('Distribution plot of Percent' , color = 'r' , fontsize = 18)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# weeks \nimport plotly.express as px\ndata=px.bar(x=list(train_x['Weeks'].value_counts().keys()), y=list(train_x['Weeks'].value_counts().values) )\ndata","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets see the progression of FVC by sex \nfig = px.line(train_x, 'Weeks', 'FVC', line_group='Patient', color='Sex',\n             title='Pulmonary Condition Progression by Sex')\nfig.update_traces(mode='lines + markers')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# now lets check the pulmonary condition progression with respect to sex\nfig = px.line(train_x, 'Weeks', 'FVC', line_group='Patient', color='SmokingStatus',\n             title='Pulmonary Condition Progression by Smoking Status')\nfig.update_traces(mode='lines+markers')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Now friends the wait is over lets pre-process the DICOM files **","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets see the no of unique patient \nprint('The Number of Unique Patients in training data are : {}'.format(len(train_x['Patient'].unique()), \"\\n\"))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = '../input/osic-pulmonary-fibrosis-progression/train/'\n\noutput_path = '../input/output/'\ntrain_image_files = sorted(glob.glob(os.path.join(data_path, '*','*.dcm')))\npatients = os.listdir(data_path)\npatients.sort()\n\nprint('Some sample Patient ID''s :', len(train_image_files))\nprint(\"\\n\".join(train_image_files[:5]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Now lets create two helper functions \n### 1. load_scan will load all DICOM images from a folder into a list for manipulation.\n### 2. The voxel values in the images are raw. get_pixels_hu converts raw values into Houndsfeld units\n### 3. The transformation is linear. Therefore, so long as you have a slope and an intercept, you can rescale a voxel value to HU.\n### 4. Both the rescale intercept and rescale slope are stored in the DICOM header at the time of image acquisition (these values are scanner-dependent, so you will need external information).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    \"\"\"\n    Loads scans from a folder and into a list.\n    \n    Parameters: path (Folder path)\n    \n    Returns: slices (List of slices)\n    \"\"\"\n    \n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    \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    return slices\ndef get_pixels_hu(scans):\n    \"\"\"\n    Converts raw images to Hounsfield Units (HU).\n    \n    Parameters: scans (Raw images)\n    \n    Returns: image (NumPy array)\n    \"\"\"\n    \n    image = np.stack([s.pixel_array for s in scans])\n    image = image.astype(np.int16)\n\n    # Since the scanning equipment is cylindrical in nature and image output is square,\n    # we set the out-of-scan pixels to 0\n    image[image == -2000] = 0\n    \n    \n    # HU = m*P + b\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# we need to know what exactly is Housnfield unit .\n## The Hounsfield unit (HU) scale is a linear transformation of the original linear attenuation coefficient measurement into one in which the radiodensity of distilled water at standard pressure and temperature (STP) is defined as zero Hounsfield units (HU), while the radiodensity of air at STP is defined as -1000 HU.\n<img src =  'https://pbrainmd.files.wordpress.com/2015/10/hounsfield-2.jpg' >\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_patient_scans = load_scan(data_path + patients[2])\ntest_patient_images = get_pixels_hu(test_patient_scans)\n\n#We'll be taking a random slice to perform segmentation:\n\nfor imgs in range(len(test_patient_images[0:5])):\n    f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,15))\n    ax1.imshow(test_patient_images[imgs], cmap=plt.cm.bone)\n    ax1.set_title(\"Original Slice\")\n    \n    ax2.imshow(test_patient_images[imgs], cmap=plt.cm.bone)\n    ax2.set_title(\"Original Slice\")\n    \n    ax3.imshow(test_patient_images[imgs], cmap=plt.cm.bone)\n    ax3.set_title(\"Original Slice\")\n    plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# animated scan ","execution_count":null},{"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\n\n\nscans = load_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/')\nscan_array = set_lungwin(get_pixels_hu(scans))\n\nimageio.mimsave(\"/tmp/gif.gif\", scan_array, duration=0.00001)\nImage(filename=\"/tmp/gif.gif\", format='png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\ntrain_x.shape ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_x.shape ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# osic laplace function ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def eval_metric(FVC,FVC_Pred,sigma):\n    n = len(sigma)\n    a=np.empty(n)\n    a.fill(70)\n    sigma_clipped = np.maximum(sigma,a) \n    delta = np.minimum(np.abs(FVC,FVC_Pred),1000)\n    eval_metric = -np.sqrt(2)*delta/sigma_clipped - np.log(np.sqrt(2)*sigma_clipped)\n    return eval_metric","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# data wrangling and processing for tabular data ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## CHECK SUBMISSION FORMAT\nsub_df = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\n\nprint(f\"The sample submission contains: {sub_df.shape[0]} rows and {sub_df.shape[1]} columns.\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# split Patient_Week Column and re-arrage columns\nsub_df[['Patient','Weeks']] = sub_df.Patient_Week.str.split(\"_\",expand = True)\nsub_df =  sub_df[['Patient','Weeks','Confidence', 'Patient_Week']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df = sub_df.merge(test_x.drop('Weeks', axis = 1), on = \"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# introduce a column to indicate the source (train/test) for the data\ntrain_x['Source'] = 'train'\nsub_df['Source'] = 'test'\n\ndata_df = train_x.append([sub_df])\ndata_df.reset_index(inplace = True)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first big challenge is data wrangling: We could see that some patients take FVE measurements only after their baseline CT-Images, and some took measurements before that. So let's first find out what the actual baseline-week and baseline-FVC for each Patient is.\nWe start with the baseline week:\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_week(df):\n    # make a copy to not change original df    \n    _df = df.copy()\n    # ensure all Weeks values are INT and not accidentaly saved as string\n    _df['Weeks'] = _df['Weeks'].astype(int)\n    # as test data is containing all weeks, \n    _df.loc[_df.Source == 'test','min_week'] = np.nan\n    _df[\"min_week\"] = _df.groupby('Patient')['Weeks'].transform('min')\n    _df['baselined_week'] = _df['Weeks'] - _df['min_week']\n    \n    return _df   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = get_baseline_week(data_df)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What we can see here, is that the Patient with ID ending on \"430\" had his first FVC measure 4 weeks before the first (baseline) CT images ( = \"Weeks\" column -4) were taken. Then the patient took the next FVC measurement 9 weeks later. In the next step we need to baseline the FVC values. Note, that the BASELINE-FVC it not the minimum FVC, but the first measurement, meaning the measurement taken in the \"min_week\" or baselined_week = 0.\n\nFor getting the baselined FVC I first wrote the following straightforward function:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_FVC_old(df):\n    # copy the DF to not in-place change the original one\n    _df = df.copy()\n    # get only the rows containing the baseline (= min_weeks) and therefore the baseline FVC\n    baseline = _df.loc[_df.Weeks == _df.min_week]\n    baseline = baseline[['Patient','FVC']].copy()\n    baseline.columns = ['Patient','base_FVC']      \n    \n    # fill the df with the baseline FVC values\n    for idx in _df.index:\n        patient_id = _df.at[idx,'Patient']\n        _df.at[idx,'base_FVC'] = baseline.loc[baseline.Patient == patient_id, 'base_FVC'].iloc[0]\n    _df.drop(['min_week'], axis = 1)\n    \n    return _df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This apporach works fine, but as it contains a lot of look-ups, its slow and didn't feel right.\nBtw: there is an even worse approach: Using for row in df.iterrows() is roughly 8 times slower than using for idx in df.index.\nSo I looked up how other people solved it and I found a rough equivalent to the following function:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_FVC(df):\n    # same as above\n    _df = df.copy()\n    base = _df.loc[_df.Weeks == _df.min_week]\n    base = base[['Patient','FVC']].copy()\n    base.columns = ['Patient','base_FVC']\n    \n    # add a row which contains the cumulated sum of rows for each patient\n    base['nb'] = 1\n    base['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\n    \n    # drop all except the first row for each patient (=unique rows!), containing the min_week\n    base = base[base.nb == 1]\n    base.drop('nb', axis = 1, inplace = True)\n    \n    # merge the rows containing the base_FVC on the original _df\n    _df = _df.merge(base, on = 'Patient', how = 'left')    \n    _df.drop(['min_week'], axis = 1)\n    \n    return _df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The second apporach is using transform, which is not as known as apply, but faster for basic-operations not involving multiple columns of a dataframe. Here is an interesting post about it for those, who want to learn more: Apply vs transform.\n\nI wanted to know how much this speeds up the processing, you can find the results in the following:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def old_baseline_FVC():\n    return get_baseline_FVC_old(data_df)\n    pass\n\ndef new_baseline_FVC():\n    return get_baseline_FVC(data_df)\n    \n\nduration_old = timeit(old_baseline_FVC, number = 3)\nduration_new = timeit(new_baseline_FVC, number = 3)\n\nprint(f\"Taking the old, non-vectorized version took {duration_old / 3:.2f} sec, while the vectorized version only took {duration_new / 3:.3f} sec. That's {duration_old/duration_new:.0f} times faster!\" )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = get_baseline_FVC(data_df)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preparing the data for the Neural Network","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import OneHotEncoder , LabelEncoder\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler\nfrom sklearn.compose import ColumnTransformer\n\n# define which attributes shall not be transformed, are numeric or categorical\nno_transform_attribs = ['Patient', 'Weeks', 'min_week']\nnum_attribs = ['FVC', 'Percent', 'Age', 'baselined_week', 'base_FVC']\ncat_attribs = ['Sex', 'SmokingStatus']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def own_MinMaxColumnScaler(df, columns):\n    \"\"\"Adds columns with scaled numeric values to range [0, 1]\n    using the formula X_scld = (X - X.min) / (X.max - X.min)\"\"\"\n    for col in columns:\n        new_col_name = col + '_scld'\n        col_min = df[col].min()\n        col_max = df[col].max()        \n        df[new_col_name] = (df[col] - col_min) / ( col_max - col_min )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def own_OneHotColumnCreator(df, columns):\n    \"\"\"OneHot Encodes categorical features. Adds a column for each unique value per column\"\"\"\n    for col in cat_attribs:\n        for value in df[col].unique():\n            df[value] = (df[col] == value).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## APPLY DEFINED TRANSFORMATIONS\nown_MinMaxColumnScaler(data_df, num_attribs)\nown_OneHotColumnCreator(data_df, cat_attribs)\n\ndata_df[data_df.Source != \"train\"].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get back original data split\ntrain_df = data_df.loc[data_df.Source == 'train']\nsub = data_df.loc[data_df.Source == 'test']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay, so the second apporach (using our own implementation) was more straightforward and less code. Downside: if you want to replace the MinMaxScaler with another scaling method (RobustScaler, StdScaler), you need to implement it first.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Model & Loss\nIn this section we are going to define the loss & a first model. First we are taking care of the loss. We are trying to minimize the following:\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"######## CONFIG ########\n\n## Features\nfeatures_list = ['baselined_week_scld', 'Percent_scld', 'Age_scld', 'base_FVC_scld', 'Male', 'Female', 'Ex-smoker', 'Never smoked', 'Currently smokes']\n\n## Basics\nEPOCHS = 1000\nBATCH_SIZE = 128\n\n\n## LOSS; set tradeoff btw. Pinball-loss and adding score\n_lambda = 0.8 # 0.8 default\n\n\n## Optimizers\nADAM = tf.keras.optimizers.Adam(lr = 0.1,\n                                beta_1 = 0.9, \n                                beta_2 = 0.999,\n                                decay = 0.01)\nSGD = tf.keras.optimizers.SGD()\n\n# choose ADAM or SGD\noptimizer = ADAM\n\n\n## To-DO: Implement Callbacks for Learning Rate Schedulers\n\nlr_start   = 0.0001\nlr_max     = 0.0001 * BATCH_SIZE # higher batch size --> higher lr\nlr_min     = 0.00001\nlr_ramp_ep = EPOCHS * 0.3\nlr_sus_ep  = 0\nlr_decay   = 0.992\n\ndef test_the_scheduler(epoch):\n        if epoch < lr_ramp_ep:\n            lr = (lr_max - lr_start) / lr_ramp_ep * epoch + lr_start\n            \n        elif epoch < lr_ramp_ep + lr_sus_ep:\n            lr = lr_max\n            \n        else:\n            lr = (lr_max - lr_min) * lr_decay**(epoch - lr_ramp_ep - lr_sus_ep) + lr_min\n            \n        return lr\n\nrng = [i for i in range(EPOCHS)]\ny = [test_the_scheduler(x) for x in rng]\nplt.plot(rng, y)\nprint(\"Learning rate schedule: {:.3g} to {:.3g} to {:.3g}\".format(y[0], max(y), y[-1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create constants for the loss function\nC1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n\n# define competition metric\ndef score(y_true, y_pred):\n    \"\"\"Calculate the competition metric\"\"\"\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 = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype = tf.float32) )\n    metric = (delta / sigma_clip) * sq2 + tf.math.log(sigma_clip * sq2)\n    return K.mean(metric)\n\n# define pinball loss\ndef qloss(y_true, y_pred):\n    \"\"\"Calculate Pinball loss\"\"\"\n    # IMPORTANT: define quartiles, feel free to change here!\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    return K.mean(v)\n\n# combine competition metric and pinball loss to a joint loss function\ndef mloss(_lambda):\n    \"\"\"Combine Score and qloss\"\"\"\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda) * score(y_true, y_pred)\n    return loss\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Neural Network Model\nIn this section we build an initial neural Network. The code of this section is derived from Ulrich's great notebook, which also inspired me to change my loss to the above coded version. Please support the original Notebook creators! The chosen quartiles are simply derived by testing; using 0.25 and 0.75 leads to worse results.\n\nFor the architecture: It's good practice to use numbers of units following the schema 2^x, with x element of N (= resulting in 1,2,4,8,16,32,64,128,..).\nWe are going to use dropout for regularization and not a too broad and deep network, as the training data is very limited.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model():\n    \"Creates and returns a model\"\n    inp = Layers.Input((len(features_list),), name = \"Patient\")\n    x = Layers.Dense(128, activation = \"relu\", name = \"d1\")(inp)\n    x = Layers.Dropout(0.25)(x)\n    x = Layers.Dense(128, activation = \"relu\", name = \"d2\")(x)\n    x = Layers.Dropout(0.2)(x)\n    # predicting the \n    p1 = Layers.Dense(3, activation = \"relu\", name = \"p1\")(x)\n    # quantile adjusting p1 predictions\n    p2 = Layers.Dense(3, activation = \"relu\", name = \"p2\")(x)\n    preds = Layers.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis = 1), \n                     name = \"preds\")([p1, p2])\n    \n    model = Models.Model(inp, preds, name = \"NeuralNet\")\n    model.compile(loss = mloss(_lambda), optimizer = optimizer, metrics = [score])\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create neural Network\nneuralNet = get_model()\nneuralNet.summary()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## GET TRAINING DATA AND TARGET VALUE\n\n# get target value\ny = train_df['FVC'].values.astype(float)\n\n\n# get training & test data\nX_train = train_df[features_list].values\nX_test = sub[features_list].values\n\n# instantiate target arrays\ntrain_preds = np.zeros((X_train.shape[0], 3))\ntest_preds = np.zeros((X_test.shape[0], 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the following we want to create leak-free folds to get a robust cross-validation strategy in order to evaluate all our models & our training. The idea is to avoid having the same patient (= PatientID) in training- and in validation-Data, as this might lead to evaluate a higher CV-score for a model which is luckily learning/memorizing the data for a particular patientID which is also frequently occuring in the validation-data.\n\nThe idea on how to do that is coming from @PAB97 Pierre's great notebook (CHECK IT OUT!) Please note, that we still don't use propoer stratification based on 'Age', 'Sex', 'SmokingStatus'.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## Non-Stratified GroupKFold-split (can be further enhanced with stratification!)\n\"\"\"K-fold variant with non-overlapping groups.\nThe same group will not appear in two different folds: in this case we dont want to have overlapping patientIDs in TRAIN and VAL-Data!\nThe folds are approximately balanced in the sense that the number of distinct groups is approximately the same in each fold.\"\"\"\n\nNFOLDS = 10\ngkf = GroupKFold(n_splits = NFOLDS)\n# extract Patient IDs for ensuring \ngroups = train_df['Patient'].values\n\ncount = 0\nfor train_idx, val_idx in gkf.split(X_train, y, groups = groups):\n    count += 1\n    print(f\"FOLD {count}:\")\n    \n    # create and fit model\n    net = get_model()\n    net.fit(X_train[train_idx], y[train_idx], batch_size = BATCH_SIZE, epochs = EPOCHS, \n            validation_data = (X_train[val_idx], y[val_idx]), verbose = 0) \n    \n    # evaluate\n    print(\"Train:\", net.evaluate(X_train[train_idx], y[train_idx], verbose = 0, batch_size = BATCH_SIZE))\n    print(\"Val:\", net.evaluate(X_train[val_idx], y[val_idx], verbose = 0, batch_size = BATCH_SIZE))\n    \n    # generate predictions for the known train data and the unknown test data\n    train_preds[val_idx] = net.predict(X_train[val_idx], batch_size = BATCH_SIZE, verbose = 0)\n    \n    print(\"Predicting Test...\")\n    test_preds += net.predict(X_test, batch_size = BATCH_SIZE, verbose = 0) / NFOLDS","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the next section we are going to use the train_preds to calculate the optimized sigma, which is a measure for certainty or rather uncertainty. We can do that, as we have both: the model's estimate and the real data. We subtract the lower quartile from the upper quartile (defined in the loss function) and average it.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## FIND OPTIMIZED STANDARD-DEVIATION\nsigma_opt = mean_absolute_error(y, train_preds[:,1])\nsigma_uncertain = train_preds[:,2] - train_preds[:,0]\nsigma_mean = np.mean(sigma_uncertain)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## PREPARE SUBMISSION FILE WITH OUR PREDICTIONS\nsub['FVC1'] = test_preds[:, 1]\nsub['Confidence1'] = test_preds[:,2] - test_preds[:,0]\n\n# get rid of unused data and show some non-empty data\nsubmission = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()\nsubmission.loc[~submission.FVC1.isnull()].head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"org_test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nfor i in range(len(org_test)):\n    submission.loc[submission['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'FVC'] = org_test.FVC[i]\n    submission.loc[submission['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'Confidence'] = 70","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index = False)\n","execution_count":null,"outputs":[]},{"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}