{"cells":[{"metadata":{"_kg_hide-output":true,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"!pip install --upgrade pip","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"!pip install -q tensorflow-io","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import os\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport random\nimport re\n\n#color\nfrom colorama import Fore, Back, Style\n\nimport tensorflow as tf\n# import tensorflow_io as tfio\n\n#read dicom\nimport pydicom\n\n#visialisation\nfrom PIL import Image\nfrom IPython.display import Image as show_gif\nimport scipy.misc\nimport matplotlib\n\n#Pandas Profiling\nimport pandas_profiling as pdp\n\n#For segmentation and 3d plotting\nimport skimage\nfrom skimage import measure, morphology\nimport skimage.measure as measure\nfrom skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing\nfrom skimage.measure import label,regionprops, perimeter\nfrom skimage.morphology import binary_dilation, binary_opening\nfrom skimage.filters import roberts, sobel\nfrom skimage import measure, feature\nfrom skimage.segmentation import clear_border\nfrom skimage import data\nfrom scipy import ndimage as ndi\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.misc\n\n# Settings for pretty nice plots\nplt.style.use('seaborn-bright')\nplt.show()\n\n# https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\nfrom timeit import timeit\nfrom tqdm import tqdm\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, GroupKFold, StratifiedKFold\n\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as Layers\nimport tensorflow.keras.models as Models\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# seed function\n# https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\n\ndef seed_everything(seed): \n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h1 align=\"center\"> <span style='font-size:30px;'>&#9672;</span> OSIC Pulmonary Fibrosis Competition </h1>"},{"metadata":{},"cell_type":"markdown","source":"#### What is Idiopathic Pulmonary Fibrosis? \n\n<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/zxI4Hv6lNpA\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>\n\n\n### Purpose of this project: \n- <span style='font-size:15px;'> 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. Specifically: </span>\n    - FVC - final 3 values for each patient (only these will be used for the final score)\n    - Confidence Value - a thread about what it is here. It is measured in ml and it's \"how confident\" are you about the the estimated FVC. High value in confidence means that you are off a lot from the actual FVC, while a very low (or even 0) confidence means you're very sure on the FVC. (if somebody thinks this is not correct, please text me, I might be wrong about this)\n    \n### Additional info:\n- <span style='font-size:15px;'> **FVC** - the recorded lung capacity in ml </span>\n\n- <span style='font-size:15px;'> **Percent** - a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics </span>\n- <span style='font-size:15px;'> **Weeks** - the relative number of weeks pre/post the baseline CT (may be negative) </span>\n\n"},{"metadata":{},"cell_type":"markdown","source":"<h3 align=\"center\"> <span style='font-size:30px;'>&#9672;</span> <strong> EDA </strong> <span style='font-size:30px;'>&#9672;</span> </h3>"},{"metadata":{"trusted":true},"cell_type":"code","source":"dataDir = '../input/osic-pulmonary-fibrosis-progression'\nworkDir = '../input/output/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = pd.read_csv(os.path.join(dataDir, 'train.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\ntrain_unique_df = train_data.drop_duplicates(subset = ['Patient'], keep = 'first')\ntrain_unique_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\n# CHECK FOR DUPLICATES & DEAL WITH THEM\n# keep = False: All duplicates will be shown\ndupRows_df = train_data[train_data.duplicated(subset = ['Patient', 'Weeks'], keep = False )]\ndupRows_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'There are {dupRows_df.shape[0]} (= {dupRows_df.shape[0] / train_data.shape[0] * 100:.2f}%) duplicates.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\n# this one decides to drop the duplicate, but I can also get the average value of FVC.\n\ntrain_data.drop_duplicates(subset=['Patient','Weeks'], keep = False, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Function to plot characteristics of patients\ndef plot_columns_data(df, col, bins):\n    plt.figure(figsize = (8, 4))\n    plt.hist(df[col], bins = bins)\n    plt.xlabel('{}'.format(col))\n    plt.ylabel('Count')\n    try:\n        plt.text(50, 450, r'$\\mu={},\\ \\sigma={}$'.format(df[col].unique().mean().round(2), df[col].unique().std().round(2)), fontsize = 12)\n    except:\n        pass\n    plt.show() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(Fore.YELLOW + 'People\\'s ages, sorted: ', Style.RESET_ALL, sorted(train_data.Age.unique()))\nprint(Fore.YELLOW + 'Youngest person\\'s age:', Style.RESET_ALL, train_data.Age.unique().min())\nprint(Fore.YELLOW + 'Oldest person\\'s age  :', Style.RESET_ALL, train_data.Age.unique().max())\nprint(Fore.YELLOW + 'Mean age of people in the dataset:', Style.RESET_ALL, train_data.Age.unique().mean().round(2))\nprint(Fore.YELLOW + 'Person\\'s sex  : ', Style.RESET_ALL, train_data.Sex.unique())\nprint(Fore.YELLOW + 'Smoking status: ', Style.RESET_ALL, train_data.SmokingStatus.unique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_columns_data(train_data, 'Age', 7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_columns_data(train_data, 'Sex', 2)\ntrain_data.groupby(['Sex']).count()['Patient'].to_frame()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_columns_data(train_data, 'SmokingStatus', 3)\ntrain_data.groupby(['SmokingStatus']).count()['Patient'].to_frame()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('There are', Fore.YELLOW + '{}'.format(len(train_data[\"Patient\"].unique())), Style.RESET_ALL, 'unique patients in Train Data.', \"\\n\")\n\n# Recordings per Patient\ndata = train_data.groupby(\"Patient\")[\"Weeks\"].count().reset_index(drop=False)\n# Sort by Weeks\ndata = data.sort_values(['Weeks']).reset_index(drop=True)\nprint('Minimum recordings per patient:', Fore.YELLOW + '{}'.format(data[\"Weeks\"].min()), Style.RESET_ALL, \"\\n\" +\n      'Maximum recordings per patient:', Fore.YELLOW + '{}'.format(data[\"Weeks\"].max()), Style.RESET_ALL)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Min FVC value:', Fore.YELLOW + '{}'.format(train_data.FVC.unique().min()), Style.RESET_ALL)\nprint('Max FVC value:', Style.RESET_ALL, Fore.YELLOW + '{}'.format(train_data.FVC.unique().max()), Style.RESET_ALL)\nprint('Min Percent value:', Fore.MAGENTA + '{} %'.format(train_data.Percent.unique().min().round(2)), Style.RESET_ALL)\nprint('Max Percent value:', Fore.MAGENTA + '{} %'.format(train_data.Percent.unique().max().round(2)), Style.RESET_ALL)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Figure\nfig, (fvc, pct) = plt.subplots(1, 2, figsize = (16, 6))\n\nfvc = sns.distplot(train_data[\"FVC\"], ax=fvc, color = 'orange')\npct = sns.distplot(train_data[\"Percent\"], ax=pct, color = 'blue')\n\nfvc.set_title(\"FVC Distribution\", fontsize=16)\npct.set_title(\"Percent Distribution\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Min no. weeks pre(negative number)/post base CT scan:', Fore.YELLOW + '{}'.format(train_data['Weeks'].min()), Style.RESET_ALL)\nprint('Max no. weeks pre(negative number)/post base CT scan:', Fore.YELLOW + '{}'.format(train_data['Weeks'].max()), Style.RESET_ALL)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (16, 6))\n\nw = sns.distplot(train_data['Weeks'])\nplt.title(\"Number of weeks before/after the CT scan\", fontsize = 16)\nplt.xlabel(\"Weeks\", fontsize=14);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Let's check the correlations\nplt.figure(figsize=(16,6))\nsns.heatmap(train_data.corr(), \n            annot=True, \n            linewidths = .5,\n            square=True,\n            annot_kws={'size':12, 'weight': 'bold'},\n            center = 0,\n            cmap=sns.color_palette('bright'))\nplt.yticks(rotation = 0)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since the Percent field approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics it is normal that we see higher correlation between these two features. As a whole the correlation between the features is close to 0."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data = pd.read_csv(os.path.join(dataDir, 'test.csv'))\n# test_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_sub = pd.read_csv(os.path.join(dataDir, 'sample_submission.csv'))\n#sample_sub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Explore the DICOM Files"},{"metadata":{"trusted":true},"cell_type":"code","source":"trainDicomPath = '../input/osic-pulmonary-fibrosis-progression/train'\ntestDicomPath = '../input/osic-pulmonary-fibrosis-progression/test'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def getListOfFiles(dirName):\n    ''' For the given path, get the List of \n        all files in the directory tree \n    '''\n    # create a list of file and sub directories \n    # names in the given directory \n    listOfFile = os.listdir(dirName)\n    allFiles = list()\n    # Iterate over all the entries\n    for entry in listOfFile:\n        # Create full path\n        fullPath = os.path.join(dirName, entry)\n        # If entry is a directory then get the list of files in this directory \n        if os.path.isdir(fullPath):\n            allFiles = allFiles + getListOfFiles(fullPath)\n        else:\n            allFiles.append(fullPath)\n                \n    return allFiles   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get the list of all train files in directory tree at given path\nlistOfTrainDcmFiles = getListOfFiles(trainDicomPath)\nlen(listOfTrainDcmFiles)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get the list of all test files in directory tree at given path\nlistOfTestDcmFiles = getListOfFiles(testDicomPath)\nlen(listOfTestDcmFiles)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Show random scans\nrandom.seed(22)\nfor file in random.sample(listOfTrainDcmFiles, 3):\n    dataset = pydicom.dcmread(file)\n    \n    print(Fore.MAGENTA + \"Patient id.......:\", Style.RESET_ALL, dataset.PatientID, \"\\n\" +\n          Fore.MAGENTA + \"Modality.........:\", Style.RESET_ALL, dataset.Modality, \"\\n\" +\n          Fore.MAGENTA + \"Rows.............:\", Style.RESET_ALL, dataset.Rows, \"\\n\" +\n          Fore.MAGENTA + \"Columns..........:\", Style.RESET_ALL, dataset.Columns)\n    \n    plt.imshow(dataset.pixel_array, cmap=plt.cm.bone)\n    plt.show()\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### #A breath-in + hold for a Patient. \n#### *Credit: https://www.kaggle.com/andradaolteanu/pulmonary-fibrosis-competition-eda-dicom-prep*"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430\"\ndatasets = []\n\n# First Order the files in the dataset\nfiles = []\nfor dcm in list(os.listdir(patient_dir)):\n    files.append(dcm) \nfiles.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n# Read in the Dataset\nfor dcm in files:\n    path = patient_dir + \"/\" + dcm\n    datasets.append(pydicom.dcmread(path))\n\n# Plot the images\nfig=plt.figure(figsize=(16, 6))\ncolumns = 10\nrows = 3\n\nfor i in range(1, columns*rows +1):\n    img = datasets[i-1].pixel_array\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=plt.cm.bone)\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create base director for Train .dcm files\ndirector = \"../input/osic-pulmonary-fibrosis-progression/train\"\n\n# Create path column with the path to each patient's CT\ntrain_data[\"Path\"] = director + \"/\" + train_data[\"Patient\"]\n\n# Create variable that shows how many CT scans each patient has\ntrain_data[\"CT_number\"] = 0\n\nfor k, path in enumerate(train_data[\"Path\"]):\n    train_data[\"CT_number\"][k] = len(os.listdir(path))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_gif(number_of_CT = 87):\n    \"\"\"Picks a patient at random and creates a GIF with their CT scans.\"\"\"\n    \n    # Select one of the patients\n    # patient = \"ID00007637202177411956430\"\n    patient = train_data[train_data[\"CT_number\"] == number_of_CT].sample(random_state=1)[\"Patient\"].values[0]\n    \n    # === READ IN .dcm FILES ===\n    patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/\" + patient\n    datasets = []\n\n    # First Order the files in the dataset\n    files = []\n    for dcm in list(os.listdir(patient_dir)):\n        files.append(dcm) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Read in the Dataset from the Patient path\n    for dcm in files:\n        path = patient_dir + \"/\" + dcm\n        datasets.append(pydicom.dcmread(path))\n        \n        \n    # === SAVE AS .png ===\n    # Create directory to save the png files\n    if os.path.isdir(f\"png_{patient}\") == False:\n        os.mkdir(f\"png_{patient}\")\n\n    # Save images to PNG\n    for i in range(len(datasets)):\n        img = datasets[i].pixel_array\n        matplotlib.image.imsave(f'png_{patient}/img_{i}.png', img)\n        \n        \n    # === CREATE GIF ===\n    # First Order the files in the dataset (again)\n    files = []\n    for png in list(os.listdir(f\"../working/png_{patient}\")):\n        files.append(png) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Create the frames\n    frames = []\n\n    # Create frames\n    for file in files:\n    #     print(\"../working/png_images/\" + name)\n        new_frame = Image.open(f\"../working/png_{patient}/\" + file)\n        frames.append(new_frame)\n\n    # Save into a GIF file that loops forever\n    frames[0].save(f'gif_{patient}.gif', format='GIF',\n                   append_images=frames[1:],\n                   save_all=True,\n                   duration=200, loop=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create_gif(number_of_CT=12)\n# create_gif(number_of_CT=87)\n\n# print(\"First file len:\", len(os.listdir(\"../working/png_ID00165637202237320314458\")), \"\\n\" +\n      \n#       \"Third file len:\", len(os.listdir(\"../working/png_ID00340637202287399835821\")))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# show_gif(filename=\"./gif_ID00165637202237320314458.gif\", format='png', width=400, height=400)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# show_gif(filename=\"./gif_ID00340637202287399835821.gif\", format='png', width=400, height=400)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pandas Profiling"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"profile_train_df = pdp.ProfileReport(train_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"profile_train_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lung segmentation and 3d plot\n#### *Reference for this section: https://www.kaggle.com/arturscussel/lung-segmentation-and-candidate-points-generation*"},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_ct_scan(folder_name):\n        # Read the slices from the dicom file\n        slices = [pydicom.dcmread(folder_name + filename) for filename in os.listdir(folder_name)]\n        \n        # Sort the dicom slices in their respective order\n        slices.sort(key=lambda x: int(x.InstanceNumber))\n        \n        # Get the pixel values for all the slices\n        slices = np.stack([s.pixel_array for s in slices])\n        slices[slices == -2000] = 0\n        return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ct_scan_1 = read_ct_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/') \nct_scan_2 = read_ct_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00047637202184938901501/')\nct_scan_3 = read_ct_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00012637202177665765362/')\nct_scan_4 = read_ct_scan('../input/osic-pulmonary-fibrosis-progression/train/ID00062637202188654068490/')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ct_scan_2.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_ct_scan(scan):\n    f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(10, 10))\n    for i in range(0, scan.shape[0], 5):\n        plots[int(i / 20), int((i % 20) / 5)].axis('off')\n        plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_ct_scan(ct_scan_2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_segmented_lungs(im, plot=False):\n    \n    '''\n    This funtion segments the lungs from the given 2D slice.\n    '''\n    if plot == True:\n        f, plots = plt.subplots(8, 1, figsize=(5, 40))\n    '''\n    Step 1: Convert into a binary image. \n    '''\n    binary = im < 604\n    if plot == True:\n        plots[0].axis('off')\n        plots[0].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 2: Remove the blobs connected to the border of the image.\n    '''\n    cleared = clear_border(binary)\n    if plot == True:\n        plots[1].axis('off')\n        plots[1].imshow(cleared, cmap=plt.cm.bone) \n    '''\n    Step 3: Label the image.\n    '''\n    label_image = label(cleared)\n    if plot == True:\n        plots[2].axis('off')\n        plots[2].imshow(label_image, cmap=plt.cm.bone) \n    '''\n    Step 4: Keep the labels with 2 largest areas.\n    '''\n    areas = [r.area for r in regionprops(label_image)]\n    areas.sort()\n    if len(areas) > 2:\n        for region in regionprops(label_image):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       label_image[coordinates[0], coordinates[1]] = 0\n    binary = label_image > 0\n    if plot == True:\n        plots[3].axis('off')\n        plots[3].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 5: Erosion operation with a disk of radius 2. This operation is \n    seperate the lung nodules attached to the blood vessels.\n    '''\n    selem = disk(2)\n    binary = binary_erosion(binary, selem)\n    if plot == True:\n        plots[4].axis('off')\n        plots[4].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 6: Closure operation with a disk of radius 10. This operation is \n    to keep nodules attached to the lung wall.\n    '''\n    selem = disk(10)\n    binary = binary_closing(binary, selem)\n    if plot == True:\n        plots[5].axis('off')\n        plots[5].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 7: Fill in the small holes inside the binary mask of lungs.\n    '''\n    edges = roberts(binary)\n    binary = ndi.binary_fill_holes(edges)\n    if plot == True:\n        plots[6].axis('off')\n        plots[6].imshow(binary, cmap=plt.cm.bone) \n    '''\n    Step 8: Superimpose the binary mask on the input image.\n    '''\n    get_high_vals = binary == 0\n    im[get_high_vals] = 0\n    if plot == True:\n        plots[7].axis('off')\n        plots[7].imshow(im, cmap=plt.cm.gray) \n        \n    return im\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_lung_from_ct_scan(ct_scan):\n    return np.asarray([get_segmented_lungs(slice) for slice in ct_scan])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def further_remove_noise(segmented_ct_scan):\n    selem = ball(2)\n    binary = binary_closing(segmented_ct_scan, selem)\n\n    label_scan = label(binary)\n\n    areas = [r.area for r in regionprops(label_scan)]\n    areas.sort()\n\n    for r in regionprops(label_scan):\n        max_x, max_y, max_z = 0, 0, 0\n        min_x, min_y, min_z = 1000, 1000, 1000\n    \n        for c in r.coords:\n            max_z = max(c[0], max_z)\n            max_y = max(c[1], max_y)\n            max_x = max(c[2], max_x)\n        \n            min_z = min(c[0], min_z)\n            min_y = min(c[1], min_y)\n            min_x = min(c[2], min_x)\n        if (min_z == max_z or min_y == max_y or min_x == max_x or r.area > areas[-3]):\n            for c in r.coords:\n                segmented_ct_scan[c[0], c[1], c[2]] = 0\n        else:\n            index = (max((max_x - min_x), (max_y - min_y), (max_z - min_z))) / (min((max_x - min_x), (max_y - min_y) , (max_z - min_z)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_ct_scan_1 = segment_lung_from_ct_scan(ct_scan_1)\nsegmented_ct_scan_1[segmented_ct_scan_1 < 604] = 0\nfurther_remove_noise(segmented_ct_scan_1)\n# plot_ct_scan(segmented_ct_scan_1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_ct_scan_2 = segment_lung_from_ct_scan(ct_scan_2)\nsegmented_ct_scan_2[segmented_ct_scan_2 < 604] = 0\nfurther_remove_noise(segmented_ct_scan_2)\n# plot_ct_scan(segmented_ct_scan_2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_ct_scan_3 = segment_lung_from_ct_scan(ct_scan_3)\nsegmented_ct_scan_3[segmented_ct_scan_3 < 604] = 0\nfurther_remove_noise(segmented_ct_scan_3)\n#plot_ct_scan(segmented_ct_scan_3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_ct_scan_4 = segment_lung_from_ct_scan(ct_scan_4)\nsegmented_ct_scan_4[segmented_ct_scan_4 < 604] = 0\nfurther_remove_noise(segmented_ct_scan_4)\n#plot_ct_scan(segmented_ct_scan_4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image):\n    \n    # Position the scan upright, \n    # so the head of the patient would be at the top facing the camera\n    p = image.transpose(2,1,0)\n    p = p[:,:,::-1]\n    \n    verts, faces, normals, values = measure.marching_cubes_lewiner(p)\n    \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    #     ax.set_xlim(0, p.shape[0])\n    #     ax.set_ylim(0, p.shape[1])\n    #     ax.set_zlim(0, p.shape[2])\n    \n    ax.set_xlim(np.min(verts[:,0]), np.max(verts[:,0]))\n    ax.set_ylim(np.min(verts[:,1]), np.max(verts[:,1])) \n    ax.set_zlim(np.min(verts[:,2]), np.max(verts[:,2]))\n    \n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], alpha=0.1)\n    \n    face_color = [0.45, 0.45, 0.75]    \n    mesh.set_facecolor(face_color)\n    \n    mesh.set_edgecolor('k')\n    \n    ax.add_collection3d(mesh)\n    \n    plt.tight_layout()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **The normal FVC range for an adult is between 3.0 and 5.0 L.** (https://www.verywellhealth.com/forced-expiratory-capacity-measurement-914900)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_ct_scan_1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.loc[train_data.Patient == 'ID00007637202177411956430'].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_ct_scan_2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.loc[train_data.Patient == 'ID00047637202184938901501'].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_ct_scan_3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.loc[train_data.Patient == 'ID00012637202177665765362'].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_ct_scan_4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.loc[train_data.Patient == 'ID00062637202188654068490'].mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A lung nodule (or mass) is a small abnormal area that is sometimes found during a CT scan of the chest. These scans are done for many reasons, such as part of lung cancer screening, or to check the lungs if you have symptoms. Most lung nodules seen on CT scans are not cancer. They are more often the result of old infections, scar tissue, or other causes. But tests are often needed to be sure a nodule is not cancer. [Reference](https://www.cancer.org/cancer/lung-cancer/detection-diagnosis-staging/lung-nodules.html#:~:text=A%20lung%20nodule%20(or%20mass,CT%20scans%20are%20not%20cancer.)) \n\nOn its [web page](https://www.uwhealth.org/thoracic-surgery/lung-nodules/11643#:~:text=Lung%20nodules%20are%20masses%20which,%2C%20diagnosed%2C%20and%20treated%20appropriately.) UW Health mentions that *Lung nodules ...can be early indicators of cancer, infection, or other lung diseases such as Pulmonary Fibrosis (PF)*.\n\n[This paper](https://www.journalpulmonology.org/en-novelties-in-imaging-in-pulmonary-articulo-S2531043719301795) stresses out the importance of CT scans for detecting Pulmonary fibrosis and precise measuring and characterizing of nodules. But it does not link between PF and nodules. I would appreciate if someone can give more information about it.\nIt can be seen that patients with FVC below 3L have many nodules, compared to the ones that have FVC in the normal range of 3 to 5L.\nI cannot claim this is representative result since I have tested it only for four randomly chosed IDs but it is worth checking if it is true."},{"metadata":{"trusted":true},"cell_type":"markdown","source":"TODO: Check correlations\n\nTODO: Convert categories to numbers\n\nTODO: Age wise smokers correlations"},{"metadata":{},"cell_type":"markdown","source":"## Preprocessing of DICOM files to prepare for model\nBelow we use the preprocessing tools from this great kernel: https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial"},{"metadata":{"trusted":true},"cell_type":"code","source":"from pydicom import dcmread\nfrom scipy import ndimage\nfrom skimage import morphology","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the scans in given folder path\ndef load_scan(path):\n    slices = [dcmread(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":{"trusted":true},"cell_type":"code","source":"listOfTrainDcmFiles","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"madical_images = pydicom.read_file(listOfTrainDcmFiles[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"madical_images","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = madical_images.pixel_array","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(image.min())\nprint(image.max())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def transform_to_hu(medical_image, image):\n    intercept = medical_image.RescaleIntercept\n    slope = medical_image.RescaleSlope\n    hu_image = image * slope + intercept\n\n    return hu_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"hu_image = transform_to_hu(madical_images, image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(hu_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_image(image, window_center, window_width):\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    window_image = image.copy()\n    window_image[window_image < img_min] = img_min\n    window_image[window_image > img_max] = img_max\n    \n    return window_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def remove_noise(file_path, display=False):\n    medical_image = pydicom.read_file(file_path)\n    image = medical_image.pixel_array\n    \n    hu_image = transform_to_hu(medical_image, image)\n    brain_image = window_image(hu_image, 40, 80)\n\n    # morphology.dilation creates a segmentation of the image\n    # If one pixel is between the origin and the edge of a square of size\n    # 5x5, the pixel belongs to the same class\n    \n    # We can instead use a circule using: morphology.disk(2)\n    # In this case the pixel belongs to the same class if it's between the origin\n    # and the radius\n    \n    segmentation = morphology.dilation(brain_image, np.ones((5, 5)))\n    labels, label_nb = ndimage.label(segmentation)\n    \n    label_count = np.bincount(labels.ravel().astype(np.int))\n    # The size of label_count is the number of classes/segmentations found\n    \n    # We don't use the first class since it's the background\n    label_count[0] = 0\n    \n    # We create a mask with the class with more pixels\n    # In this case should be the brain\n    mask = labels == label_count.argmax()\n    \n    # Improve the brain mask\n    mask = morphology.dilation(mask, np.ones((5, 5)))\n    mask = ndimage.morphology.binary_fill_holes(mask)\n    mask = morphology.dilation(mask, np.ones((3, 3)))\n    \n    # Since the the pixels in the mask are zero's and one's\n    # We can multiple the original image to only keep the brain region\n    masked_image = mask * brain_image\n\n    if display:\n        plt.figure(figsize=(20, 5))\n        plt.subplot(141)\n        plt.imshow(brain_image)\n        plt.title('Original Image')\n        plt.axis('off')\n        \n        plt.subplot(142)\n        plt.imshow(mask)\n        plt.title('Mask')\n        plt.axis('off')\n\n        plt.subplot(143)\n        plt.imshow(masked_image)\n        plt.title('Final Image')\n        plt.axis('off')\n    \n    return masked_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"remove_noise(listOfTrainDcmFiles[11], display = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_and_plot_image(file_path, save=False):\n    medical_image = pydicom.read_file(file_path)\n    image = medical_image.pixel_array\n    \n    print(image.shape)\n    \n    hu_image = transform_to_hu(medical_image, image)\n    brain_image = window_image(hu_image, 40, 80)\n    bone_image = window_image(hu_image, 400, 1000)\n    \n    plt.figure(figsize=(30, 15))\n    plt.style.use('grayscale')\n\n    plt.subplot(151)\n    plt.imshow(image)\n    plt.title('Original')\n    plt.axis('off')\n\n    plt.subplot(152)\n    plt.imshow(hu_image)\n    plt.title('Hu image')\n    plt.axis('off')\n\n    plt.subplot(153)\n    plt.imshow(brain_image)\n    plt.title('brain image')\n    plt.axis('off')\n\n    plt.subplot(154)\n    plt.imshow(bone_image)\n    plt.title('bone image')\n    plt.axis('off')\n    \n    if save:\n        mpimg.imsave(os.path.join(output_path, f'{file_path[:-4]}-original.png'), image)\n        mpimg.imsave(os.path.join(output_path, f'{file_path[:-4]}-hu_image.png'), hu_image)\n        mpimg.imsave(os.path.join(output_path, f'{file_path[:-4]}-brain_image.png'), brain_image)\n        mpimg.imsave(os.path.join(output_path, f'{file_path[:-4]}-bone_image.png'), bone_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"load_and_plot_image(listOfTrainDcmFiles[11])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preprocessing for the model\n#### Credit: https://www.kaggle.com/ChristianDenich/quantile-reg-lr-schedulers-checkpoints\n##### TODO: Do some corrections later\n##### TODO: Write down every step that you plan to take. Much easier and don't need to try to remember what you had to do.\n##### This code only uses the tabular data but not the image data. Work on that. Will probably need to concatenate models. Check out other solutions."},{"metadata":{"trusted":true},"cell_type":"code","source":"# check submission data\nsub_data = pd.read_csv(os.path.join(dataDir, 'sample_submission.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# split Patient_Week Column and re-arrage columns\nsub_data[['Patient','Weeks']] = sub_data.Patient_Week.str.split(\"_\",expand = True)\nsub_data = sub_data[['Patient','Weeks','Confidence', 'Patient_Week']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_data = sub_data.merge(test_data.drop('Weeks', axis = 1), on = \"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.sample()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data.sample()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# introduce a column to indicate the source (train/test) for the data\ntrain_data['Source'] = 'train'\nsub_data['Source'] = 'test'\n\ndata_df = train_data.append([sub_data])\ndata_df.reset_index(inplace = True)\ndata_df.drop('Path', axis=1, 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:"},{"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    _df['min_week'] = _df['Weeks']\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. \n\nNote, that the *BASELINE-FVC* is not the minimum FVC, but the first measurement, meaning *the measurement taken in the \"min_week\" or baselined_week = 0.*\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get the baselined FVC\n# https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/179033\n# there is one more aproach, but slower. The comparison can be seen in the notebook here: https://www.kaggle.com/chrisden/6-82-quantile-reg-lr-schedulers-checkpoints\n\ndef 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\n\n# The nicest and fastest one\n\ndef get_baseline_FVC_new(df):\n    base = (\n        df\n        .loc[df.Weeks == df.min_week][['Patient','FVC']]\n        .rename({'FVC': 'base_FVC'}, axis=1)\n        .groupby('Patient')\n        .first()\n        .reset_index()\n    )\n    \n    # merge the rows containing the base_FVC on the original _df\n    _df = df.copy()\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":{"trusted":true},"cell_type":"code","source":"data_df = get_baseline_FVC_new(data_df)\ndata_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Standardization and normalization. Preparing for Neural network\n**chrisden** offers two approaches. One with creating classes and using sckit learn and another where I create my one scaling and one-hot encoding. I will use the first approach since I would like to try to work with classes more. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# import the necessary Encoders & Transformers\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler\nfrom sklearn.compose import ColumnTransformer","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define which attributes shall not be transformed, are numeric or categorical\nno_transform_attribs = ['Patient', 'Weeks', 'min_week', 'Source']\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":{},"cell_type":"markdown","source":"To get a ColumnTransformer who outputs the whole DataFrame compatible format, we need a class that takes attributes which we dont want to change, and simply passes them through. This class needs a fit and a transform method, so that the ColumnTransformer itself can use fit_transform like for the numerical and categorical attributes."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# class NoTransformer(BaseEstimator, TransformerMixin):\n#     \"\"\"Passes through data without any change and is compatible with ColumnTransformer class\"\"\"\n#     def fit(self, X, y=None):\n#         return self\n\n#     def transform(self, X):\n#         assert isinstance(X, pd.DataFrame)\n#         return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ## GET TRANSFORMED DATAFRAME\n\n# # create an instance of the ColumnTransformer\n# datawrangler = ColumnTransformer(([\n#      # the No-Transformer does not change the data and is applied to all no_transform_attribs \n#      ('original', NoTransformer(), no_transform_attribs),\n#      # Apply MinMax to the numerical attributes, here you can change to e.g. StdScaler()   \n#      ('MinMax', MinMaxScaler(), num_attribs),\n#      # OneHotEncoder all categorical attributes.   \n#      ('cat_encoder', OneHotEncoder(), cat_attribs),\n#     ]))\n\n# transformed_data_series = []\n# transformed_data_series = datawrangler.fit_transform(data_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Okay, now we encoded all the data and want to have a look at it. But wait, we only get back a series element? How to put that in a Dataframe again? Sadly, not as easy as I had hoped. We need to use pd.DataFrame(data, columns=column_names) function, for which we need the column names. But we used OneHot-Encoding. So the number of columns now depends on how many different values/categories a categorical value has, because for each unique value we get a separate column: e.g. If for the column \"SmokingStatus\" we only have the values \"Smoker\" and \"Never-Smoked\", we have two resulting columns if there are additional possible values like \"Ex_Smoker\" we get more columns. And we also should not get things wrong and mix the columns up. The following code is getting our data back to a Dataframe and preserving the correct order."},{"metadata":{"trusted":true},"cell_type":"code","source":"# # get column names for non-categorical data\n# new_col_names = no_transform_attribs + num_attribs\n\n# # extract possible values from the fitted transformer\n# categorical_values = [s for s in datawrangler.named_transformers_[\"cat_encoder\"].get_feature_names()]\n# new_col_names += categorical_values\n\n# # create Dataframe based on the extracted Column-Names\n# train_sklearn_df = pd.DataFrame(transformed_data_series, columns = new_col_names)\n# train_sklearn_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# new_col_names","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":{"trusted":true},"cell_type":"code","source":"train_df.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loss - see the formula in the credited notebook\n\nIn this section you can configure the following:\n\n* Features used for training\n* Basic training setup: BATCH_SIZE and EPOCHS,\n* Configuration for the loss function\n* Optimizers, Learning-Rate-Schedulers incl. Learning Rate start- & endpoint\n* Custom Logging Callback\n* Checkpoint-Saving Callback\nThe Learning-Rate scheduler below is inspired by Chris great [Melanoma-detection notebook](https://www.kaggle.com/cdeotte/triple-stratified-kfold-with-tfrecords).\nFeel free to experiment with the scheduler and it's max/min and decay values.\n\nEver wondered why lr_max is scaled by BATCH_SIZE and therefore bigger for larger batches? The reason for this is the following: the larger the BATCH_SIZE, the more averaged & smoothened a step of gradient decent is and the bigger our confidence in the direction of the step is. As there is less \"randomness\" in a huge averaged batch (compared with for example Stochastic Gradient Decent (=SGD) with batch size = 1) and our confidence in the direction is higher, the learning rate can be bigger to advance fast to the optimum."},{"metadata":{},"cell_type":"markdown","source":"######## CONFIG ########"},{"metadata":{},"cell_type":"raw","source":"# be careful, the resulsts are very SEED-DEPENDEND!\nseed_everything(1989)"},{"metadata":{"trusted":true},"cell_type":"code","source":"### Features: choose which features you want to use\n# you can exclude and include features by extending this feature list\nfeatures_list = ['baselined_week_scld', 'Age_scld', 'base_FVC_scld', 'Male', 'Female', 'Ex-smoker', 'Never smoked', 'Currently smokes']\n\n### Basics for training:\nEPOCHS = 1500\nBATCH_SIZE = 256\n\n\n### LOSS; set tradeoff btw. Pinball-loss and adding score\n_lambda = 0.8 # 0.8 default\n\n\n### Optimizers\n# choose ADAM or SGD\noptimizer = 'SGD'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Learning Rate Scheduler\ndef get_lr_callback(batch_size = 64, plot = False):\n    \"\"\"Returns a lr_scheduler callback which is used for training.\n    Feel free to change the values below!\n    \"\"\"\n    LR_START   = 0.001\n    LR_MAX     = 0.0001 * BATCH_SIZE # higher batch size --> higher lr\n    LR_MIN     = 0.000001\n    # 30% of all epochs are used for ramping up the LR and then declining starts\n    LR_RAMP_EP = EPOCHS * 0.3\n    # how many epochs shall LR_MAX be sustained\n    LR_SUS_EP  = 0\n    # rate of decay\n    LR_DECAY   = 0.993\n\n    def lr_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    \n    if plot == False:\n        # get the Keras-required callback with our LR for training\n        lr_callback = tf.keras.callbacks.LearningRateScheduler(lr_scheduler,verbose = False)\n        return lr_callback \n    \n    else: \n        return lr_scheduler","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot & check the LR-Scheulder for sanity-check\nlr_scheduler_plot = get_lr_callback(batch_size = 64, plot = True)\n\nrng = [i for i in range(EPOCHS)]\n\ny = [lr_scheduler_plot(x) for x in rng]\n\nplt.plot(rng, y)\n\nprint(f\"Learning rate schedule: {y[0]:.3f} to {max(y):.3f} to {y[-1]:.3f}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# logging & saving\nLOGGING = True\n\n# defining custom callbacks\nclass LogPrintingCallback(tf.keras.callbacks.Callback):\n    \n    def on_train_begin(self, logs = None):\n        print(\"Training started\")\n        # self.val_loss = [] not used for now\n        self.val_score = []        \n        \n    def on_epoch_end(self, epoch, logs = None):\n        # self.val_loss.append(logs['val_loss']) not used for now\n        self.val_score.append(logs['val_score'])\n        if epoch % 250 == 0 or epoch == (EPOCHS -1 ):\n            print(f\"The average val-loss for epoch {epoch} is {logs['val_loss']:.2f}\"\n                  f\" and the val-score is {logs['val_score']}\")\n            \n    def on_train_end(self, lowest_val_loss, logs = None):\n        # get index of best epoch\n        best_epoch = np.argmin(self.val_score)\n        # get score in best epoch\n        best_score = self.val_score[best_epoch]\n        print(f\"Stop training, best model was found and saved in epoch {best_epoch + 1} with val-score: {best_score}.\"\n              f\" Final results in this fold (last epoch):\") ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_checkpont_saver_callback(fold):\n    checkpt_saver = tf.keras.callbacks.ModelCheckpoint(\n        'fold-%i.h5'%fold,\n        monitor = 'val_score',\n        verbose = 0,\n        save_best_only = True,\n        save_weights_only = True,\n        mode = 'min',\n        save_freq = 'epoch')\n    \n    return checkpt_saver","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"####### loss function #######"},{"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    # Python is automatically broadcasting y_true with shape (1,0) to \n    # shape (3,0) in order to make this subtraction work\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.20, 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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Neural Network Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow_addons as tfa","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model(optimizer = 'ADAM', lr = 0.01):\n    \"Creates and returns a model\"\n    # instantiate optimizer\n    optimizer = tf.keras.optimizers.Adam(lr = LR) if optimizer == 'ADAM' else tf.keras.optimizers.SGD()\n    \n    # create model    \n    inp = Layers.Input((len(features_list),), name = \"Patient\")\n    x = Layers.BatchNormalization()(inp)\n    x = tfa.layers.WeightNormalization(Layers.Dense(160, activation = \"elu\", name = \"d1\"))(x)\n    x = Layers.BatchNormalization()(x)\n    x = Layers.Dropout(0.3)(x)\n    x = tfa.layers.WeightNormalization(Layers.Dense(128, activation = \"elu\", name = \"d2\"))(x)\n    x = Layers.BatchNormalization()(x)\n    x = Layers.Dropout(0.25)(x)\n    # predicting the 3 quantiles\n    q1 = Layers.Dense(3, activation = \"relu\", name = \"p1\")(x)\n    # generating another output for quantile adjusting the quantile predictions\n    q_adjust = Layers.Dense(3, activation = \"relu\", name = \"p2\")(x)\n    \n    # adding the tf.cumsum of q_adjust to the output q1\n    # to ensure increasing values [a < b < c]\n    # tf.cumsum([a, b, c]) --> [a, a + b, a + b + c]\n    preds = Layers.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis = 1), \n                     name = \"preds\")([q1, q_adjust])\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(optimizer, lr = 0.01)\nneuralNet.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train_df.drop(['Source'], axis = 1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","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# get training & test data\nX_train = train_df[features_list].values.astype(float)\nX_test = sub[features_list].values.astype(float)\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":{"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 = 7\ngkf = GroupKFold(n_splits = NFOLDS)\n# extract Patient IDs for ensuring \ngroups = train_df['Patient'].values\n\nOOF_val_score = []\nfold = 0\n\nfor train_idx, val_idx in gkf.split(X_train, y, groups = groups):\n    fold += 1\n    print(f\"FOLD {fold}:\")\n    \n    # callbacks: logging & model saving with checkpoints each fold\n    # callbacks = [get_lr_callback(BATCH_SIZE)]  # un-comment for using LRScheduler\n    reduce_lr_loss = tf.keras.callbacks.ReduceLROnPlateau(monitor = 'val_loss',\n                                                          factor = 0.6,\n                                                          patience = 150,\n                                                          verbose = 1,\n                                                          epsilon = 1e-4,\n                                                          mode = 'min')\n    \n    callbacks = [reduce_lr_loss]\n    \n    if LOGGING == True:\n        callbacks +=  [get_checkpont_saver_callback(fold),                     \n                     LogPrintingCallback()]\n\n    # build and train model\n    model = get_model(optimizer, lr = 0.05)\n    history = model.fit(X_train[train_idx], y[train_idx], \n              batch_size = BATCH_SIZE, \n              epochs = EPOCHS, \n              validation_data = (X_train[val_idx], y[val_idx]), \n              callbacks = callbacks,\n              verbose = 0) \n    \n    # evaluate\n    print(\"Train:\", model.evaluate(X_train[train_idx], y[train_idx], verbose = 0, batch_size = BATCH_SIZE, return_dict = True))\n    print(\"Val:\", model.evaluate(X_train[val_idx], y[val_idx], verbose = 0, batch_size = BATCH_SIZE, return_dict = True))\n    \n    ## Load best model to make pred\n    model.load_weights('fold-%i.h5'%fold)\n    train_preds[val_idx] = model.predict(X_train[val_idx],\n                                         batch_size = BATCH_SIZE,\n                                         verbose = 0)\n    \n    # append OOF evaluation to calculate OFF_Score\n    OOF_val_score.append(model.evaluate(X_train[val_idx], y[val_idx], verbose = 0, batch_size = BATCH_SIZE, return_dict = True)['score'])\n    \n    # predict on test set and average the predictions over all folds\n    print(\"Predicting Test...\")\n    test_preds += model.predict(X_test, batch_size = BATCH_SIZE, verbose = 0) / NFOLDS","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## PLOT results\n# fetch results from history\nscore = history.history['score']\nval_score = history.history['val_score']\n\nloss = history.history['loss']\nval_loss = history.history['val_loss']\n\nepochs_range = range(EPOCHS)\n\n# create plots\nplt.figure(figsize = (12,5))\nplt.plot(epochs_range, score, label = 'Training Score')\nplt.plot(epochs_range, val_score, label = 'Validation Score')\n# limit y-values for better zoom-scale. Remember that roughly -4.5 is the best possible score\n# plt.ylim(0.8 * np.mean(val_score), 1.2 * np.mean(val_score))\nplt.legend(loc = 'lower right')\nplt.title('Training and Validation Score')\nplt.show()\n\nplt.figure(figsize = (12,5))\nplt.plot(epochs_range, loss, label = 'Training Loss')\nplt.plot(epochs_range, val_loss, label = 'Validation Loss')\n# limit y-values for beter zoom-scale\n# plt.ylim(0.3 * np.mean(val_loss), 1.8 * np.mean(val_loss))\n\nplt.legend(loc = 'upper right')\nplt.title('Training and Validation Loss')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# OOF evaluation\n\nnp.mean(OOF_val_score)","execution_count":null,"outputs":[]},{"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.loc[~submission.FVC1.isnull(),'FVC'] = submission.loc[~submission.FVC1.isnull(),'FVC1']\n\nif sigma_mean < 70:\n    submission['Confidence'] = sigma_opt\nelse:\n    submission.loc[~submission.FVC1.isnull(),'Confidence'] = submission.loc[~submission.FVC1.isnull(),'Confidence1']","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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}