{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Overview & Remarks\n\n- In this notebook we estimate the linear decay in FVC for each patient over the weeks\n- The tissue within the segmented lung was extracted using horizontal and vertical perturbations\n- Several features like Lung Volume, Average tissue area for each patient etc. were extracted from the image\n- Average tissue area was found to show a high positive correlation with the decline in FVC over time, hence a couple of features related to tissue were used in the model.\n- EffiecientNet-B1 was trained with 30 epochs modified callbacks and training parameters\n- Mloss = 0.65"},{"metadata":{},"cell_type":"markdown","source":"# Acknowledgements\n\n- This notebook is primarily based on Wei Hao Khoong's work which can be found here: https://www.kaggle.com/khoongweihao/efficientnets-quantile-regression-inference\n\n- The lung segmentation technique used below is a modification of several common techniques present in various other kernels\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Imports"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index\n!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Regular Imports\nimport numpy as np\nimport pandas as pd\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm \nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import GroupKFold, KFold, StratifiedKFold\nimport seaborn as sns\nimport math\nimport cv2\nimport pydicom\nimport os\nimport glob\nimport pickle as pkl\nimport matplotlib.image as mpimg\nfrom tabulate import tabulate\nimport missingno as msno \nfrom IPython.display import display_html\nfrom PIL import Image\nimport gc\nfrom skimage.transform import resize\nimport copy\nimport re\nfrom scipy.stats import pearsonr\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n%matplotlib inline\n\n# Segmentation\nimport glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly.tools import FigureFactory as FF\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) \n\n# Model imports\nimport tensorflow as tf \nfrom tensorflow.keras.layers import (\n                                    Dense, Dropout, Activation, Flatten, Input, BatchNormalization, GlobalAveragePooling2D,\n                                    Add, Conv2D, AveragePooling2D, LeakyReLU, Concatenate , Lambda\n                                    )\nfrom tensorflow.keras import Model\nimport tensorflow.keras.models as M\nfrom tensorflow.keras.utils import Sequence\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.applications as tfa\nimport efficientnet.tfkeras as efn\nfrom sklearn.model_selection import train_test_split, KFold\nimport seaborn as sns\n\n\npd.set_option(\"display.max_columns\", 100)\ncustom_colors = ['#74a09e','#86c1b2','#98e2c6','#f3c969','#f2a553', '#d96548', '#c14953']\nsns.palplot(sns.color_palette(custom_colors))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"config = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"EPOCHS = 30\nBATCH_SIZE = 8\nNFOLD = 5\nLR = 0.003\nSAVE_BEST = True\nMODEL_CLASS = 'b1'\npath = '../input/osic-pulmonary-fibrosis-progression/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Helper Functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Areas with the same number of pixels on the edges are not required. Crop it.\ndef crop_image(img: np.ndarray):\n    edge_pixel_value = img[0, 0]\n    mask = img != edge_pixel_value\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\n# Load images, crop thick borders(if any) and resize\ndef load_image(path):\n    dataset = pydicom.dcmread(path)\n    img = dataset.pixel_array\n    img = crop_image(img)\n#     img = cv2.resize(img, (512,512))\n    return img\n\n# Get Nth percentile image\ndef get_img(perc, patient_id, data):\n    \n    l = glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id))\n    img_ids = []\n    for x in l:\n        y = x.split('/')[-1]\n        z = int(y.split('.')[0])\n        img_ids.append(z)\n\n    img_ids.sort()\n    \n    return img_ids[math.ceil(perc*(len(img_ids)))-1]\n\n# Get num of slices bw two percentiles\ndef num_img_bw_perc(p1, p2, patient_id, data):\n    \n    l = glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id))\n    img_ids = []\n    for x in l:\n        y = x.split('/')[-1]\n        z = int(y.split('.')[0])\n        img_ids.append(z)\n\n    img_ids.sort()\n    \n    return len(img_ids[math.ceil(p1*(len(img_ids)))-1:math.ceil(p2*(len(img_ids)))])-1\n    \n\n# Get number of images per patient\ndef get_num_images(patient_id, data):\n    \n    return len(glob.glob(path+'{0}/{1}/*.dcm'.format(data, patient_id)))\n\n# Get the lung area in the image slice\ndef lung_seg_pixel_ratio(img_array):\n    \n    c = 0\n    for i in range(img_array.shape[0]):\n        for j in range(img_array.shape[1]):\n            if img_array[i][j] != 0:\n                c+=1\n    \n    return c, round(c/(img_array.shape[0]*img_array.shape[1]),4)\n\n# https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\ndef make_lungmask(img, display=False):\n    \n    img = img.astype(float)\n    \n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    \n    # Find the average pixel value near the lungs\n    # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    \n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    \n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    \n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if ((B[2]-B[0]<row_size*0.9) and (B[3]-B[1]<col_size*0.9) and (B[2]-B[0]>row_size*0.20) \n            and (B[3]-B[1]>col_size*0.10) and (B[0]>row_size*0.03) and (B[2]<row_size*0.97) \n            and (B[1]>col_size*0.03) and (B[3]<col_size*0.97)):\n#         if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    \n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    # Compute Lung Area in the slice\n    lung_pixels, slice_lung_area = lung_seg_pixel_ratio(mask)\n    \n    # Tissue Mask\n    t_mask = tissue_mask(img, mask, shift_perc = 0.02)\n    \n    # Extract tissue features\n    num_t_pixels, tissue_by_total, tissue_by_lung = tissue_features(t_mask, img, thresh = 0.35)\n    \n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n        \n    return lung_pixels, slice_lung_area, num_t_pixels, tissue_by_total, tissue_by_lung\n\n# Get tissue mask\ndef tissue_mask(img, mask, shift_perc):\n    \n    r_dim, c_dim = img.shape[0], img.shape[1]\n    \n    # Move the image by shift_perc to the left\n    del_left_cols = int(shift_perc*c_dim)\n    \n    mask1, mask2 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask1 = mask1[:,del_left_cols:]\n    mask2[:,:c_dim-del_left_cols] = mask1\n    \n    # Move the image by shift_perc to the right\n    del_right_cols = int(shift_perc*c_dim)\n    \n    mask3, mask4 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask3 = mask3[:,:c_dim-del_right_cols]\n    mask4[:,del_right_cols:] = mask3\n    \n    # Move the image by shift_perc to the top\n    del_top_rows = int(shift_perc*c_dim)\n    \n    mask5, mask6 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask5 = mask5[del_top_rows:,:]\n    mask6[:r_dim-del_top_rows,:] = mask5\n    \n    # Move the image by shift_perc to the bottom\n    del_bottom_rows = int(shift_perc*r_dim)\n    \n    mask7, mask8 = mask.copy(), np.zeros((r_dim, c_dim)).astype(int)\n    mask7 = mask7[:r_dim-del_bottom_rows,:]\n    mask8[del_bottom_rows:,:] = mask7\n\n    #Obtain the final mask\n    final_mask = ((mask2==1) & (mask4==1) & (mask6==1) & (mask8==1)).astype(int)\n    \n    return final_mask\n    \n# Get tissue features\ndef tissue_features(tissue_mask, img, thresh = 0.35):\n    \n    final_img = tissue_mask*img\n\n    checker = np.zeros((final_img.shape[0], final_img.shape[1]))\n    counter, other_counter = 0, 0\n    for i in range(final_img.shape[0]):\n        for j in range(final_img.shape[1]):\n            if final_img[i][j]>=thresh:\n                checker[i][j] = 1\n                counter+=1\n            else:\n                checker[i][j] = 0\n                other_counter+=1\n                \n    tissue_by_total = counter/(final_img.shape[0]**2)\n    tissue_by_lung = counter/((tissue_mask==1).sum())\n                              \n#     if tissue_by_lung>1:\n#         tissue_by_total = (counter-(tissue_mask==0).sum())/(rescaled.shape[0]**2)\n#         tissue_by_lung = (counter-(tissue_mask==0).sum())/((tissue_mask==1).sum())\n#         counter = counter-(tissue_mask==0).sum()\n        \n    return counter, tissue_by_total, tissue_by_lung\n\n# Get dicom meta data\ndef get_dicom_meta(path):\n    \n    '''Get information from the .dcm files.\n    path: complete path to the .dcm file'''\n\n    image_data = pydicom.read_file(path)\n\n    # Dictionary to store the information from the image\n    observation_data = {\n                        \"SliceThickness\" : float(image_data.SliceThickness),\n                        \"PixelSpacing\" : float(image_data.PixelSpacing[0]),\n                        }\n    return observation_data\n\ndef seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    \n    \ndef get_tab(df):\n    vector = [(df.Age.values[0] - 30) / 30] \n    \n    if df.Sex.values[0].lower() == 'male':\n        vector.append(0)\n    else:\n        vector.append(1)\n    \n    if df.SmokingStatus.values[0] == 'Never smoked':\n        vector.extend([0,0])\n    elif df.SmokingStatus.values[0] == 'Ex-smoker':\n        vector.extend([1,1])\n    elif df.SmokingStatus.values[0] == 'Currently smokes':\n        vector.extend([0,1])\n    else:\n        vector.extend([1,0])\n        \n    if df.Avg_Tissue_30_60_Quartile.values[0] == 'Q1':\n        vector.extend([0])\n    elif df.Avg_Tissue_30_60_Quartile.values[0] == 'Q2':\n        vector.extend([1])\n    elif df.Avg_Tissue_30_60_Quartile.values[0] == 'Q3':\n        vector.extend([2])\n    elif df.Avg_Tissue_30_60_Quartile.values[0] == 'Q4':\n        vector.extend([3])\n    else:\n        vector.extend([4])\n        \n    vector.extend([df.Avg_Tissue_30_60.values[0]])\n    \n    return np.array(vector) \n\ndef get_img_1(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize((d.pixel_array - d.RescaleIntercept) / (d.RescaleSlope * 1000), (512, 512))\n\ndef score(fvc_true, fvc_pred, sigma):\n    sigma_clip = np.maximum(sigma, 70) # changed from 70, trie 66.7 too\n    delta = np.abs(fvc_true - fvc_pred)\n    delta = np.minimum(delta, 1000)\n    sq2 = np.sqrt(2)\n    metric = (delta / sigma_clip)*sq2 + np.log(sigma_clip* sq2)\n    return np.mean(metric)\n\nseed_everything(42)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Test Data And Train Data With Image Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = pd.read_csv('../input/osic-pf-train-meta-image/train_df_with_img_feat.csv')\n    \nts = pd.read_csv(f\"{path}/test.csv\")\n\nts['Where'] = 'test'\n\n# Fetch unique patient ids\nts_patient_ids = ts.Patient.unique().tolist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Image Feature Extraction"},{"metadata":{"trusted":true},"cell_type":"code","source":"percentile_range_1 = np.linspace(0.3,0.6,11)\npercentile_range_2 = range(30,63, 3)\npercentile_range = [(round(x,2),y) for x,y in zip(percentile_range_1,percentile_range_2)]\npercentile_range","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nimage_data = []\ncounter = 1\n\n\nfor data in [('test',ts_patient_ids)]:\n    \n    half_img_path = []\n    flag = -1\n\n    for i in data[1]:\n        \n        lung_list,num_t_pixels_list,tissue_by_total_list,tissue_by_lung_list = [], [], [], []\n        \n        for perc in percentile_range:\n            \n            img_no = get_img(perc[0], i, data[0])\n            \n            if img_no == flag:\n                print(perc[0], i, data[0])\n                continue\n            \n            flag = img_no\n            \n            try:\n                img= load_image(path+'{0}/{1}/{2}.dcm'.format(data[0],i,img_no))\n            except:\n                print('{0}/{1}'.format(i,img_no), 'Err')\n                continue\n                \n            try:\n                lung_pixels, area_ratio, num_t_pixels, tissue_by_total, tissue_by_lung = make_lungmask(img, display=False)\n                \n                if math.isnan(lung_pixels):\n                    print('nan lung_pixels',data[0],i,img_no)\n                else:\n                    lung_list.append((perc[0],lung_pixels))\n                \n                if math.isnan(num_t_pixels):\n                    print('nan num_t_pixels',data[0],i,img_no)\n                else:\n                    num_t_pixels_list.append(num_t_pixels)\n                    \n                if math.isnan(tissue_by_total):\n                    print('nan tissue_by_total',data[0],i,img_no)\n                else:\n                    tissue_by_total_list.append(tissue_by_total)\n                \n                if math.isnan(tissue_by_lung):\n                    print('nan tissue_by_lung',data[0],i,img_no)\n                else:\n                    tissue_by_lung_list.append(tissue_by_lung)    \n                \n#                 lung_list.append((perc[0],lung_pixels))\n#                 num_t_pixels_list.append(num_t_pixels)\n#                 tissue_by_total_list.append(tissue_by_total)\n#                 tissue_by_lung_list.append(tissue_by_lung)\n            except Exception as e:\n                print(data[0], i, img_no, e)\n                pass\n        \n        slice_thickness = get_dicom_meta(path+'{0}/{1}/{2}.dcm'.format(data[0], i, img_no))['SliceThickness']\n        pixel_spacing = get_dicom_meta(path+'{0}/{1}/{2}.dcm'.format(data[0], i, img_no))['PixelSpacing']\n        \n        try:\n            Avg_NumTissuePixel_30_60 = round(sum(num_t_pixels_list)/len(num_t_pixels_list),4)\n            Avg_Tissue_30_60 = round((sum(num_t_pixels_list)/len(num_t_pixels_list))*pixel_spacing,4)\n            Avg_Tissue_thickness_30_60 = round((sum(num_t_pixels_list)/len(num_t_pixels_list))*pixel_spacing*slice_thickness,4)\n            Avg_TissueByTotal_30_60 = round(sum(tissue_by_total_list)/len(tissue_by_total_list),4)\n            Avg_TissueByLung_30_60 = round(sum(tissue_by_lung_list)/len(tissue_by_lung_list),4)\n        except Exception as e:\n            print(data[0], i, img_no, e)\n            Avg_NumTissuePixel_30_60 = 0\n            Avg_Tissue_30_60 = 0\n            Avg_Tissue_thickness_30_60 = 0\n            Avg_TissueByTotal_30_60 = 0\n            Avg_TissueByLung_30_60 = 0\n        \n        num_im = num_img_bw_perc(percentile_range[0][0], percentile_range[1][0], i, data[0])\n        \n        if num_im == 0:\n            num_im = 1\n        \n        approx_vol = 0\n        for x in lung_list[:-1]:\n            approx_vol+=(x[1]*pixel_spacing*slice_thickness*num_im)\n            \n        patient_dict = {'Patient':i,\n                        'Data':data[0],\n                        'SliceThickness': slice_thickness,\n                        'PixelSpacing': pixel_spacing,\n                        'NumImgBw5Prec': num_im,\n                        'ApproxVol_30_60':round(approx_vol,4),\n                        'Avg_NumTissuePixel_30_60':Avg_NumTissuePixel_30_60,\n                        'Avg_Tissue_30_60':Avg_Tissue_30_60,\n                        'Avg_Tissue_thickness_30_60':Avg_Tissue_thickness_30_60,\n                        'Avg_TissueByTotal_30_60':Avg_TissueByTotal_30_60,\n                        'Avg_TissueByLung_30_60':Avg_TissueByLung_30_60}\n        \n        image_data.append(patient_dict)\n        \n        if counter%20 == 0:\n            print(counter)\n        counter+=1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_data_df = pd.DataFrame(image_data)\nimage_data_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Linear Decay - CNN"},{"metadata":{},"cell_type":"markdown","source":"## Merge Image and Meta "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Merge test data with image features\nts = ts.merge(image_data_df, left_on = ['Patient'], right_on = ['Patient'], how = 'left')\n\n# Merge train and test data\ndf = pd.concat([tr,ts])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, axes = plt.subplots(3, 2, figsize=(16, 12))\n\nsns.kdeplot(df.ApproxVol_30_60, shade = True, ax=axes[0, 0])\naxes[0,0].set_title('Average approximate lung volume')\nsns.kdeplot(df.Avg_NumTissuePixel_30_60, shade = True, ax=axes[0, 1])\naxes[0,1].set_title('Average number of tissue pixels')\nsns.kdeplot(df.Avg_Tissue_30_60, shade = True, ax=axes[1, 0])\naxes[1,0].set_title('Average tissue pixel area')\nsns.kdeplot(df.Avg_Tissue_thickness_30_60, shade = True, ax=axes[1, 1])\naxes[1,1].set_title('Average tissue pixel volume')\nsns.kdeplot(df.Avg_TissueByTotal_30_60, shade = True, ax=axes[2, 0])\naxes[2,0].set_title('Tissue area by Total area ratio')\nsns.kdeplot(df.Avg_TissueByLung_30_60, shade = True, ax=axes[2, 1])\naxes[2,1].set_title('Tissue area by Lung area ratio')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Clean the outliers"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Impute 0 values with median\nvariables = ['ApproxVol_30_60','Avg_NumTissuePixel_30_60','Avg_Tissue_30_60','Avg_Tissue_thickness_30_60','Avg_TissueByTotal_30_60','Avg_TissueByLung_30_60']\nfor var in variables:\n    median = df[var].quantile(q = 0.5, interpolation='linear')\n    df.loc[df[var]==0,var] = median\n\n# Impute extreme outliers in Avg_TissueByLung_20_60 by median\nmedian = df['Avg_TissueByLung_30_60'].quantile(q = 0.5, interpolation='linear')\ndf.loc[df['Avg_TissueByLung_30_60']>=0.3,'Avg_TissueByLung_30_60'] = median\n\n# Impute extreme outliers in Avg_TissueByLung_20_60 by median\nmedian = df['Avg_Tissue_thickness_30_60'].quantile(q = 0.5, interpolation='linear')\ndf.loc[df['Avg_Tissue_thickness_30_60']>=20000,'Avg_Tissue_thickness_30_60'] = median\n\n#Conver Avg_Tissue_30_60 to quartiles\ndf[\"Avg_Tissue_30_60_Quartile\"] = pd.qcut(df.Avg_Tissue_30_60, q = 4, labels = ['Q1','Q2','Q3','Q4'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, axes = plt.subplots(3, 2, figsize=(16, 12))\n\nsns.kdeplot(df.ApproxVol_30_60, shade = True, ax=axes[0, 0])\naxes[0,0].set_title('Average approximate lung volume')\nsns.kdeplot(df.Avg_NumTissuePixel_30_60, shade = True, ax=axes[0, 1])\naxes[0,1].set_title('Average number of tissue pixels')\nsns.kdeplot(df.Avg_Tissue_30_60, shade = True, ax=axes[1, 0])\naxes[1,0].set_title('Average tissue pixel area')\nsns.kdeplot(df.Avg_Tissue_thickness_30_60, shade = True, ax=axes[1, 1])\naxes[1,1].set_title('Average tissue pixel volume')\nsns.kdeplot(df.Avg_TissueByTotal_30_60, shade = True, ax=axes[2, 0])\naxes[2,0].set_title('Tissue area by Total area ratio')\nsns.kdeplot(df.Avg_TissueByLung_30_60, shade = True, ax=axes[2, 1])\naxes[2,1].set_title('Tissue area by Lung area ratio')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Normalize the data"},{"metadata":{"trusted":true},"cell_type":"code","source":"l = ['SliceThickness','PixelSpacing','ApproxVol_30_60','Avg_NumTissuePixel_30_60','Avg_Tissue_30_60',\n     'Avg_Tissue_thickness_30_60','Avg_TissueByTotal_30_60','Avg_TissueByLung_30_60']\n\nfor var in l:\n    df[var] = (df[var] - df[var].min() ) / ( df[var].max() - df[var].min() )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, axes = plt.subplots(3, 2, figsize=(16, 12))\n\nsns.kdeplot(df.ApproxVol_30_60, shade = True, ax=axes[0, 0])\naxes[0,0].set_title('Average approximate lung volume')\nsns.kdeplot(df.Avg_NumTissuePixel_30_60, shade = True, ax=axes[0, 1])\naxes[0,1].set_title('Average number of tissue pixels')\nsns.kdeplot(df.Avg_Tissue_30_60, shade = True, ax=axes[1, 0])\naxes[1,0].set_title('Average tissue pixel area')\nsns.kdeplot(df.Avg_Tissue_thickness_30_60, shade = True, ax=axes[1, 1])\naxes[1,1].set_title('Average tissue pixel volume')\nsns.kdeplot(df.Avg_TissueByTotal_30_60, shade = True, ax=axes[2, 0])\naxes[2,0].set_title('Tissue area by Total area ratio')\nsns.kdeplot(df.Avg_TissueByLung_30_60, shade = True, ax=axes[2, 1])\naxes[2,1].set_title('Tissue area by Lung area ratio')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = df[df.Where == 'train']\ntest = df[df.Where == 'test']\nprint(train.shape)\nprint(test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"quartile = ['Q1','Q2','Q3','Q4']\nfor q in quartile:\n    temp_df = train.loc[train.Avg_Tissue_30_60_Quartile==q,:]\n    patients = temp_df.Patient.unique()\n    fvc_diff_lst = []\n    print('Quartile: ', q)\n    for patient in patients:\n        fvc_lst = temp_df.loc[temp_df.Patient == patient, 'FVC'].tolist()\n        fvc_diff = fvc_lst[-1]-fvc_lst[0]\n        fvc_diff_lst.append(fvc_diff)\n    print('FVC Diff Mean: ', sum(fvc_diff_lst)/len(fvc_diff_lst))\n    print('FVC Diff Q1: ', np.quantile(fvc_diff_lst, .25))\n    print('FVC Diff Q2/Median: ', np.quantile(fvc_diff_lst, .5))\n    print('FVC Diff Q3: ', np.quantile(fvc_diff_lst, .75))\n    print('FVC Diff Q4: ', np.quantile(fvc_diff_lst, .90))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We observe that higher the presence of lung tissue, higher is the decline in FVC"},{"metadata":{},"cell_type":"markdown","source":"## Modeling"},{"metadata":{"trusted":true},"cell_type":"code","source":"A = {} \nTAB = {} \nP = [] \nfor i, p in tqdm(enumerate(train.Patient.unique())):\n    sub = train.loc[train.Patient == p, :] \n    fvc = sub.FVC.values\n    weeks = sub.Weeks.values\n    c = np.vstack([weeks, np.ones(len(weeks))]).T\n    a, b = np.linalg.lstsq(c, fvc)[0]\n    \n    A[p] = a\n    TAB[p] = get_tab(sub)\n    P.append(p)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## CNN for coeff prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.utils import Sequence\n\nclass IGenerator(Sequence):\n    BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n    def __init__(self, keys, a, tab, batch_size=32):\n        self.keys = [k for k in keys if k not in self.BAD_ID]\n        self.a = a\n        self.tab = tab\n        self.batch_size = batch_size\n        \n        self.train_data = {}\n        for p in train.Patient.values:\n            self.train_data[p] = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n    \n    def __len__(self):\n        return 1000\n    \n    def __getitem__(self, idx):\n        x = []\n        a, tab = [], [] \n        keys = np.random.choice(self.keys, size = self.batch_size)\n        for k in keys:\n            try:\n                i = np.random.choice(self.train_data[k], size=1)[0]\n                img = get_img_1(f'../input/osic-pulmonary-fibrosis-progression/train/{k}/{i}')\n                x.append(img)\n                a.append(self.a[k])\n                tab.append(self.tab[k])\n            except:\n                print(k, i)\n       \n        x,a,tab = np.array(x), np.array(a), np.array(tab)\n        x = np.expand_dims(x, axis=-1)\n        return [x, tab] , a","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.layers import (\n    Dense, Dropout, Activation, Flatten, Input, BatchNormalization, GlobalAveragePooling2D, Add, Conv2D, AveragePooling2D, \n    LeakyReLU, Concatenate \n)\nimport efficientnet.tfkeras as efn\n\ndef get_efficientnet(model, shape):\n    models_dict = {\n        'b0': efn.EfficientNetB0(input_shape=shape,weights=None,include_top=False),\n        'b1': efn.EfficientNetB1(input_shape=shape,weights=None,include_top=False),\n        'b2': efn.EfficientNetB2(input_shape=shape,weights=None,include_top=False),\n        'b3': efn.EfficientNetB3(input_shape=shape,weights=None,include_top=False),\n        'b4': efn.EfficientNetB4(input_shape=shape,weights=None,include_top=False),\n        'b5': efn.EfficientNetB5(input_shape=shape,weights=None,include_top=False),\n        'b6': efn.EfficientNetB6(input_shape=shape,weights=None,include_top=False),\n        'b7': efn.EfficientNetB7(input_shape=shape,weights=None,include_top=False)\n    }\n    return models_dict[model]\n\ndef build_model(shape=(512, 512, 1), model_class = None, fold=None):\n    inp = Input(shape=shape)\n    base = get_efficientnet(model_class, shape)\n    x = base(inp)\n    x = GlobalAveragePooling2D()(x)\n    inp2 = Input(shape=(6,))\n    x2 = tf.keras.layers.GaussianNoise(0.2)(inp2)\n    x = Concatenate()([x, x2]) \n    x = Dropout(0.5)(x) \n    x = Dense(1)(x)\n    model = Model([inp, inp2] , x)\n    \n    weights = [w for w in os.listdir('../input/lindecay-efficientnetb1-models') if str(fold) in w][0]\n    model.load_weights('../input/lindecay-efficientnetb1-models/' + weights)\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# from time import time\n# start_time = time()\n# for fold in range(5):\n#     model = build_model(shape=(512, 512, 1), model_class = 'b1', fold=fold)\n#     print(time()-start_time)\n#     start_time = time()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import fold-wise quantiles which give us best score\n# import pickle as pkl\n\n# with open('../input/lindecay-efficientnetb3-models/bst_quantile.pkl','rb') as f:\n#     bst_quantile = pkl.load(f)\n\nbst_quantile = [0.8, 0.5, 0.1, 0.1, 0.1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_class = 'b1'\nNFOLD = 5\n\nkf = KFold(n_splits=NFOLD, random_state=42,shuffle=False)\nP = np.array(P)\nsubs = []\nfolds_history = []\n\nfor fold, (tr_idx, val_idx) in enumerate(kf.split(P)):\n    print('#####################')\n    print('####### Fold %i ######'%fold)\n    print('#####################')\n    print('Predicting...')\n    \n    model = build_model(shape=(512, 512, 1), model_class = model_class, fold = fold)\n    q = bst_quantile[fold]\n    \n    sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \n    A_test, B_test, P_test,W, FVC= {}, {}, {}, {}, {} \n    STD, WEEK = {}, {} \n    for p in test.Patient.unique():\n        x = [] \n        tab = [] \n        ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/')\n        counter = 0\n        for i in ldir:\n            if int(i[:-4]) / len(ldir) < 0.7 and int(i[:-4]) / len(ldir) > 0.20:\n                x.append(get_img_1(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/{i}')) \n                tab.append(get_tab(test.loc[test.Patient == p, :])) \n                counter+=1\n            if counter == 20:\n                break\n        if len(x) <= 1:\n            continue\n        tab = np.array(tab) \n\n        x = np.expand_dims(x, axis=-1) \n        _a = model.predict([x, tab]) \n        a = np.quantile(_a, q)\n        A_test[p] = a\n        B_test[p] = test.FVC.values[test.Patient == p] - a*test.Weeks.values[test.Patient == p]\n        P_test[p] = test.Percent.values[test.Patient == p] \n        WEEK[p] = test.Weeks.values[test.Patient == p]\n\n    for k in sub.Patient_Week.values:\n        p, w = k.split('_')\n        w = int(w) \n\n        fvc = A_test[p] * w + B_test[p]\n        sub.loc[sub.Patient_Week == k, 'FVC'] = fvc\n        sub.loc[sub.Patient_Week == k, 'Confidence'] = (\n            P_test[p] - A_test[p] * abs(WEEK[p] - w) \n    ) \n\n    _sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()\n    subs.append(_sub)\n    \n    K.clear_session()\n    del model\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Averaging Predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"N = len(subs)\nsub = subs[0].copy() # ref\nsub[\"FVC\"] = 0\nsub[\"Confidence\"] = 0\nfor i in range(N):\n    sub[\"FVC\"] += subs[i][\"FVC\"] * (1/N)\n    sub[\"Confidence\"] += subs[i][\"Confidence\"] * (1/N)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_img.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Quantile Regression"},{"metadata":{},"cell_type":"markdown","source":"## Meta Data Processing"},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = pd.read_csv(f\"{path}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nchunk = pd.read_csv(f\"{path}/test.csv\")\n\n# Fetch unique patient ids\ntr_patient_ids = tr.Patient.unique().tolist()\nts_patient_ids = chunk.Patient.unique().tolist()\n\nsub = pd.read_csv(f\"{path}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\nmeta_data = tr.append([chunk, sub])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(tr.shape, chunk.shape, sub.shape, meta_data.shape)\nprint(tr.Patient.nunique(), chunk.Patient.nunique(), sub.Patient.nunique(), \n      meta_data.Patient.nunique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_data['min_week'] = meta_data['Weeks']\nmeta_data.loc[meta_data.WHERE=='test','min_week'] = np.nan","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_data['min_week'] = meta_data['Weeks']\nmeta_data.loc[meta_data.WHERE=='test','min_week'] = np.nan\nmeta_data['min_week'] = meta_data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = meta_data.loc[meta_data.Weeks == meta_data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)\nmeta_data = meta_data.merge(base, on='Patient', how='left')\n\nbase = meta_data.loc[meta_data.Weeks == meta_data.min_week]\nbase = base[['Patient','Percent']].copy()\nbase.columns = ['Patient','min_Percent']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)\nmeta_data = meta_data.merge(base, on='Patient', how='left')\n\nmeta_data['base_week'] = meta_data['Weeks'] - meta_data['min_week']\ndel base","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"COLS = ['Sex','SmokingStatus'] #,'Age'\nFE = []\nfor col in COLS:\n    for mod in meta_data[col].unique():\n        print(col, mod)\n        FE.append(mod)\n        meta_data[mod] = (meta_data[col] == mod).astype(int)\n#=================","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(meta_data.shape)\nmeta_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Prepare Image DF"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_data_df = pd.DataFrame(image_data)\nimage_data_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr1 = pd.read_csv('../input/osic-pf-train-meta-image/train_df_with_img_feat.csv')\ncols = image_data_df.columns\ntr1 = tr1[cols]\ntr1 = tr1.drop_duplicates(keep = 'first')\ntr1.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Merge train and test data\nimage_data_df = pd.concat([tr1,image_data_df])\n\nimage_data_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Clean the outliers"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Impute 0 values with median\nvariables = ['ApproxVol_30_60','Avg_NumTissuePixel_30_60','Avg_Tissue_30_60','Avg_Tissue_thickness_30_60','Avg_TissueByTotal_30_60','Avg_TissueByLung_30_60']\nfor var in variables:\n    median = image_data_df[var].quantile(q = 0.5, interpolation='linear')\n    image_data_df.loc[image_data_df[var]==0,var] = median\n\n# Impute extreme outliers in Avg_TissueByLung_20_60 by median\nmedian = image_data_df['Avg_TissueByLung_30_60'].quantile(q = 0.5, interpolation='linear')\nimage_data_df.loc[image_data_df['Avg_TissueByLung_30_60']>=0.3,'Avg_TissueByLung_30_60'] = median\n\n# Impute extreme outliers in Avg_TissueByLung_20_60 by median\nmedian = image_data_df['Avg_Tissue_thickness_30_60'].quantile(q = 0.5, interpolation='linear')\nimage_data_df.loc[image_data_df['Avg_Tissue_thickness_30_60']>=20000,'Avg_Tissue_thickness_30_60'] = median","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge Image and Meta"},{"metadata":{"trusted":true},"cell_type":"code","source":"final_data = meta_data.merge(image_data_df, left_on = ['Patient','WHERE'], right_on = ['Patient', 'Data'], how = 'left')\nprint(final_data.shape)\n# final_data = final_data.loc[final_data.Patient != 'ID00011637202177653955184',:]\nprint(final_data.shape)\nfinal_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Min-Max Normalization"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Min-Max normalization\n\nfinal_data['age'] = (final_data['Age'] - final_data['Age'].min() ) / ( final_data['Age'].max() - final_data['Age'].min() )\nfinal_data['BASE_FVC'] = (final_data['min_FVC'] - final_data['min_FVC'].min() ) / ( final_data['min_FVC'].max() - final_data['min_FVC'].min() )\nfinal_data['week'] = (final_data['base_week'] - final_data['base_week'].min() ) / ( final_data['base_week'].max() - final_data['base_week'].min() )\nfinal_data['BASE_percent'] = (final_data['min_Percent'] - final_data['min_Percent'].min() ) / ( final_data['min_Percent'].max() - final_data['min_Percent'].min() )\n\n\nl = ['SliceThickness','PixelSpacing','ApproxVol_30_60','Avg_NumTissuePixel_30_60','Avg_Tissue_30_60',\n     'Avg_Tissue_thickness_30_60','Avg_TissueByTotal_30_60','Avg_TissueByLung_30_60']\n\nfor var in l:\n    final_data[var] = (final_data[var] - final_data[var].min() ) / ( final_data[var].max() - final_data[var].min() )\n\nFE += ['age','week','BASE_FVC','BASE_percent']\nFE.extend(l)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = final_data.loc[final_data.WHERE=='train']\nchunk = final_data.loc[final_data.WHERE=='val']\nsub = final_data.loc[final_data.WHERE=='test']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"## Model Building"},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return K.mean(metric)\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    qs = [0.2, 0.50, 0.8]\n    q = tf.constant(np.array([qs]), dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return K.mean(v)\n#=============================#\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss\n#=================\ndef laplace_log_likelihood(actual_fvc, predicted_fvc, confidence, return_values = False):\n    \"\"\"\n    Calculates the modified Laplace Log Likelihood score for this competition.\n    \"\"\"\n    sd_clipped = np.maximum(confidence, 70)\n    delta = np.minimum(np.abs(actual_fvc - predicted_fvc), 1000)\n    metric = - np.sqrt(2) * delta / sd_clipped - np.log(np.sqrt(2) * sd_clipped)\n\n    if return_values:\n        return metric\n    else:\n        return np.mean(metric)\n#=================\ndef make_model(FE):\n    z = Input((len(FE),), name=\"Patient\")\n    x = Dense(100, activation=\"relu\", name=\"d1\")(z)\n#     x = L.Dropout(0.05)(x)\n    x = Dense(100, activation=\"relu\", name=\"d2\")(x)\n#     x = L.Dense(100, activation=\"relu\", name=\"d3\")(x)\n    p1 = Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n    \n    model = M.Model(z, preds, name=\"CNN\")\n    #model.compile(loss=qloss, optimizer=\"adam\", metrics=[score])\n    model.compile(loss=mloss(0.65), optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = make_model(FE)\nprint(net.summary())\nprint(net.count_params())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## GET TRAINING DATA AND TARGET VALUE\n\n# get target value\ny = tr['FVC'].values.astype(float)\n\n\n# get training & test data\nX_train = tr[FE].values\nX_test = sub[FE].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":{"trusted":true},"cell_type":"code","source":"%%time\n\nFE = ['Male', 'Female', 'Ex-smoker', 'Never smoked', 'Currently smokes', 'age', 'week', 'BASE_FVC', 'BASE_percent']\nimage_features = ['SliceThickness','PixelSpacing','ApproxVol_30_60','Avg_NumTissuePixel_30_60','Avg_Tissue_30_60',\n                  'Avg_TissueByTotal_30_60','Avg_TissueByLung_30_60']\nFE1 = FE\nFE2 = FE+['ApproxVol_30_60']\nFE3 = FE+['Avg_Tissue_thickness_30_60']\nFE4 = FE+['Avg_TissueByLung_30_60']\nFE5 = FE+['ApproxVol_30_60','Avg_Tissue_thickness_30_60','Avg_TissueByLung_30_60']\n\nmodel_cnt = 1\nfor FE in [FE1, FE2, FE3, FE4, FE5]:\n    \n    print('Model Count: ', model_cnt, '\\nFeatures Used: ', FE)\n    \n    # get target value\n    y = tr['FVC'].values.astype(float)\n\n    # get training & test data\n    X_train = tr[FE].values\n    X_test = sub[FE].values\n\n    # instantiate target arrays\n    globals()['train_preds_{}'.format(model_cnt)] = np.zeros((X_train.shape[0], 3))\n    globals()['test_preds_{}'.format(model_cnt)] = np.zeros((X_test.shape[0], 3))\n\n    NFOLD = 5\n    gkf = GroupKFold(n_splits=NFOLD)\n    groups = tr['Patient'].values\n\n    cnt = 0\n    EPOCHS = 800\n    BATCH_SIZE=128\n    for tr_idx, val_idx in gkf.split(X_train, y, groups):\n        cnt += 1\n        print(f\"FOLD {cnt}\")\n        net = make_model(FE)\n        net.fit(X_train[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, \n                validation_data=(X_train[val_idx], y[val_idx]), verbose=0) #\n        print(\"train\", net.evaluate(X_train[tr_idx], y[tr_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        print(\"predict val...\")\n        globals()['train_preds_{}'.format(model_cnt)][val_idx] = net.predict(X_train[val_idx], batch_size=BATCH_SIZE, verbose=0)\n        print(\"predict test...\")\n        globals()['test_preds_{}'.format(model_cnt)] += net.predict(X_test, batch_size=BATCH_SIZE, verbose=0) / NFOLD\n\n    predicted_fvc = globals()['train_preds_{}'.format(model_cnt)][:,1]\n    confidence = globals()['train_preds_{}'.format(model_cnt)][:,2]-globals()['train_preds_{}'.format(model_cnt)][:,0]\n    model_score = laplace_log_likelihood(actual_fvc = y, predicted_fvc = predicted_fvc, confidence = confidence,\n                           return_values = False)\n    print('Overall Score: ', model_score)\n    model_cnt+=1\n    #==============\n    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scores = []\n\nfor i in range(1,6):\n    globals()['predicted_fvc_{}'.format(i)] = globals()['train_preds_{}'.format(i)][:,1]\n    globals()['confidence_{}'.format(i)] = globals()['train_preds_{}'.format(i)][:,2]-globals()['train_preds_{}'.format(i)][:,0]\n    \nfor i1 in np.linspace(0,0.5,11):\n    for i2 in np.linspace(0,0.5,11):\n        for i3 in np.linspace(0,0.5,11):\n            for i4 in np.linspace(0,0.5,11):\n                for i5 in np.linspace(0,0.5,11):\n                    if i1+i2+i3+i4+i5 == 1:\n                        train_preds = (i1*train_preds_1+\n                                       i2*train_preds_2+\n                                       i3*train_preds_3+\n                                       i4*train_preds_4+\n                                       i5*train_preds_5)\n                        \n                        predicted_fvc = train_preds[:,1]\n                        confidence = train_preds[:,2]-train_preds[:,0]\n                        score = laplace_log_likelihood(actual_fvc = y, predicted_fvc = predicted_fvc, confidence = confidence, return_values = False)\n                        scores.append((i1,i2,i3,i4,i5,score))\n            \nscores = sorted(scores, key = lambda x: x[-1])\nscores[-10:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best_weights = scores[-1]\ntrain_preds = (best_weights[0]*train_preds_1+best_weights[1]*train_preds_2+best_weights[2]*train_preds_3+\n               best_weights[3]*train_preds_4+best_weights[4]*train_preds_5)\n\ntest_preds = (best_weights[0]*test_preds_1+best_weights[1]*test_preds_2+best_weights[2]*test_preds_3+\n               best_weights[3]*test_preds_4+best_weights[4]*test_preds_5)","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)\n\n# 161.8837622626985 237.76367962865177","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(train_preds[idxs, 0], label=\"q25\")\nplt.plot(train_preds[idxs, 1], label=\"q50\")\nplt.plot(train_preds[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","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":"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_meta.csv\", index = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_sub = submission[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Ensemble (Simple Blend)"},{"metadata":{"trusted":true},"cell_type":"code","source":"df1 = img_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)\ndf2 = reg_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df1[['Patient_Week']].copy()\ndf['FVC'] = 0.4*df1['FVC'] + 0.6*df2['FVC']\ndf['Confidence'] = 0.4*df1['Confidence'] + 0.6*df2['Confidence']\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}