{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nimport pydicom\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport plotly.express as px\nimport scipy\nfrom scipy import ndimage\n\nfrom skimage import measure\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LinearRegression\n\nimport numpy as np\n\nfrom glob import 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 *\n\nfrom scipy.ndimage import zoom\nfrom scipy.stats import kurtosis\nfrom scipy.stats import skew\n\nfrom tqdm import tqdm\n\nimport tensorflow as tf\nimport tensorflow_hub as hub\nfrom tensorflow.keras.layers.experimental import preprocessing\n\nimport keras\n\nos.getcwd()\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#Extract train.csv into a dataframe\ndata_path='../input/osic-pulmonary-fibrosis-progression'\ntraindf=pd.read_csv(f'{data_path}/train.csv')\n\ntraindf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf=pd.DataFrame({'pid': traindf['Patient'].unique()})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Converting traindf to coredf\n\nn_weeks=[]\nweeks=[]\nfvc=[]\npercent=[]\nage=[]\nsex=[]\nsmoking=[]\n\n\nfor pid in coredf['pid']:\n    tempdf=traindf[traindf['Patient']==pid]\n    \n    n_weeks.append(len(tempdf))\n    weeks.append(list(tempdf['Weeks'].values))\n    fvc.append(list(tempdf['FVC'].values))\n    percent.append(list(tempdf['Percent'].values))\n    age.append(tempdf['Age'].values[0])\n    sex.append(tempdf['Sex'].values[0])\n    smoking.append(tempdf['SmokingStatus'].values[0])\n    \n\ncoredf['n_weeks']=n_weeks\ncoredf['weeks']=weeks\ncoredf['fvc']=fvc\ncoredf['percent']=percent\ncoredf['age']=age\ncoredf['sex']=sex\ncoredf['smoking']=smoking\n\ndisplay(coredf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Adding dicom scan details to coredf\ndicoms=[]\nn_dicoms=[]\nbad_ids=[]\n\nfor pid in coredf['pid']:\n    try:\n        files=os.listdir(f\"{data_path}/train/{pid}\")\n\n        for i in range(len(files)):\n            files[i]=files[i][:-4]\n\n        dicoms.append(files)\n        n_dicoms.append(len(files))\n    \n    except:\n        bad_ids.append(pid)\n\ncoredf['dicoms']=dicoms\ncoredf['n_dicoms']=n_dicoms\n\ndisplay(coredf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf = coredf.sort_values('n_dicoms')\ncoredf=coredf.reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_scan(path,index):\n    dicom = dcmread(os.path.join(path,os.listdir(path)[0]))\n    \n    return dicom\n\ndef get_scans(path):\n    dicoms=[]\n\n    for filename in os.listdir(path):\n        dicoms.append(dcmread(os.path.join(path,filename)))\n\n    return dicoms\n\n\ndef get_images(scans):\n  temp=[]\n\n  for i in range(len(scans)):\n    try:\n        temp.append(scans[i].pixel_array)\n    except:\n        return None\n  \n  return np.asarray(temp,dtype=np.int16)\n\n'''def get_images(scans):\n  temp=[]\n\n  for i in range(len(scans)):\n    try:\n        images=np.array(scans[i].pixel_array,dtype=np.int16)\n    except:\n        return None\n\n    for image in images:\n        image[image<=-1000] = 0\n\n        intercept = scans[i].RescaleIntercept\n        slope = scans[i].RescaleSlope\n\n        if slope != 1:\n                image = slope * images.astype(np.int16)\n                image = image.astype(np.int16)\n\n        image += np.int16(intercept)\n\n\n        temp.append(image)\n        \n        \n    return np.asarray(temp,dtype=np.int16)'''\n\ndef crop_image(image):\n    if image.shape != (512,512):\n        left=((image.shape[0]-512)//2)\n        right=image.shape[0]-left\n        top=((image.shape[1]-512)//2)\n        bottom=image.shape[1]-top\n        \n        return image[left-1:right-1,top-1:bottom-1]\n    else:\n        return image\n    \ndef crop_images(images):\n    if images[0].shape != (512,512):\n        temp=[]\n        for image in images:\n            temp.append(crop_image(image))\n        \n        return np.asarray(temp)\n    else:\n        return images\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing=[scan[0].SliceThickness,scan[0].PixelSpacing[0],scan[0].PixelSpacing[1]]\n    spacing = np.array(list(spacing))\n\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)\n    \n    return image\n\n\ndef dimension_scaling(images):\n  '''\n  Returns images with channels. e.g. converts 512x512 images to 512x512x3\n\n  Inputs:\n  imageslices:  NumpyArray  Imageslices extracted from the DICOM files using fnc get_images_from_slices\n\n  '''\n\n  images=np.expand_dims(images,axis=-1)\n  images=np.concatenate((images,)*3, axis=-1)\n\n  return np.asarray(images)\n\ndef normalize(arr): \n    return (arr-arr.min())/(arr.max()-arr.min())\n    \n\ndef get_window_value(feature):\n    if type(feature) == pydicom.multival.MultiValue:\n        return np.int(feature[0])\n    else:\n        return np.int(feature)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Standardize the pixel values\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    # 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    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    #\n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    #\n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n    #\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    #\n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n        \n    return mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"m = tf.keras.Sequential([\n    hub.KerasLayer(\"https://tfhub.dev/tensorflow/efficientnet/b6/feature-vector/1\",\n                   trainable=False),  # Can be True, see below.\n    tf.keras.layers.Dense(7, activation='relu')\n])\nm.build([None, 512, 512, 3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coredf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"area=[]\npixelspacing_r = []\npixelspacing_c = []\nslice_thicknesses = []\nwindow_widths = []\nwindow_levels = []\neff_features=[]\nr_distance=[]\nc_distance=[]\narea_cm2=[]\nslice_volume_cm3=[]\n\nbad_ids=[]\n\n\nfor pid in coredf['pid'][:-1]:\n    images=get_images(get_scans(f'{data_path}/train/{pid}'))\n    \n    if type(images)!= type(None):\n        images=crop_images(images)\n        images=normalize(images)\n        \n        for i in tqdm(range(len(images)),desc=f\"{pid}\"):\n            try:\n                images[i]=images[i]*make_lungmask(images[i],False)\n            except:\n                images[i]=images[i]*np.zeros((512,512))\n        \n        images=dimension_scaling(images)\n        \n        features=np.mean(m.predict(images),axis=0)\n        \n        eff_features.append(features)\n\n    else:\n        bad_ids.append(pid)\n        eff_features.append(np.zeros(10,np.float32))\n\n    scan = get_scan(f'{data_path}/train/{pid}',0)\n    \n    area.append(scan.Columns*scan.Rows)\n    \n    window_widths.append(get_window_value(scan.WindowWidth))\n    window_levels.append(get_window_value(scan.WindowCenter))\n    \n    spacing = scan.PixelSpacing\n    pixelspacing_r.append(spacing[0])\n    pixelspacing_c.append(spacing[1])\n    \n    slice_thicknesses.append(scan.SliceThickness)\n    \n    r_distance.append(pixelspacing_r[-1] * scan.Rows)\n    c_distance.append(pixelspacing_c[-1] * scan.Columns)\n    area_cm2.append(0.1 * r_distance[-1] * 0.1 * c_distance[-1])\n    slice_volume_cm3.append(0.1 * slice_thicknesses[-1] * area_cm2[-1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Extracting model fit data\nage_list=[]\nsex_list=[]\nsmoking_list=[]\nbase_week=[]\nbase_fvc=[]\nbase_percent=[]\ncurrent_week=[]\nweeks_passed=[]\ncurrent_percent=[]\ncurrent_fvc=[]\n\nfor i in tqdm(range(len(coredf))):\n    weeks=coredf['weeks'][i]\n    fvc=coredf['fvc'][i]\n    percent=coredf['percent'][i]\n    age=coredf['age'][i]\n    \n    if coredf['sex'][i]=='Male':\n        sex=0\n    else:\n        sex=1\n\n    if coredf['smoking'][i] == 'Ex-smoker':\n        smoking=0\n    elif coredf['smoking'][i] == 'Never smoked':\n        smoking=1\n    else:\n        smoking=2\n    \n    for i in range(len(weeks)):            \n        for j in range(i+1,len(weeks)):\n            age_list.append(age)\n            sex_list.append(sex)\n            smoking_list.append(smoking)\n        \n            base_week.append(weeks[i])\n            base_percent.append(percent[i])\n            base_fvc.append(fvc[i])\n            \n            \n            weeks_passed.append(weeks[j]-weeks[i])\n            \n            \n            current_week.append(weeks[j])\n            current_percent.append(percent[j])\n            current_fvc.append(fvc[j])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"modeldf=pd.DataFrame({'age':normalize(np.asarray(age_list)),\n                      'sex':normalize(np.asarray(sex_list)),\n                      'smoking':normalize(np.asarray(smoking_list)),\n                        'base_week':normalize(np.asarray(base_week)),\n                      'base_fvc':normalize(np.asarray(base_fvc)),\n                      'base_percent': normalize(np.asarray(base_percent)),\n                      'current_week':normalize(np.asarray(current_week)),\n                      'weeks_passed':normalize(np.asarray(weeks_passed)),\n                      'current_percent': normalize(np.asarray(current_percent)),\n                      'current_fvc':np.asarray(current_fvc,dtype=np.float32)\n                     })","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"modeldf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"px.imshow(modeldf.corr())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"initializer = tf.keras.initializers.GlorotNormal()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = ((delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2))\n    \n    return keras.backend.mean(metric)\n#============================#\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 keras.backend.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#=================\n\ndef make_model():\n    ini_input=keras.Input(shape=(9,1))\n\n    y = keras.layers.LSTM(25,return_sequences=True)(ini_input)\n    y = keras.layers.LSTM(25,return_sequences=True)(y)\n    y = keras.layers.LSTM(25,return_sequences=True)(y)\n    y = keras.layers.Reshape([25*9])(y)\n    x = keras.layers.Dense(100,kernel_initializer=initializer, activation=\"relu\", name=\"d1\")(y)\n    x = keras.layers.Dense(100,kernel_initializer=initializer, activation=\"relu\", name=\"d2\")(x)\n    p1 = keras.layers.Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = keras.layers.Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = keras.layers.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n\n    model=keras.Model(inputs=ini_input,outputs=preds)\n\n    model.compile(loss=mloss(0.8),optimizer=tf.keras.optimizers.Adam(lr=0.2, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False),metrics=[score])\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\n\nkf = KFold(n_splits=5,shuffle=True)\n\ntrainx = modeldf[modeldf.columns[:-1]].values\ntrainy = modeldf[modeldf.columns[-1]].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"b_size=128  \ncount=1\nhistory={}\n\nfor train_ind,val_ind in kf.split(trainx):\n    print(f'Fold No. {count}')\n    \n    model=make_model()\n    \n    temp={}\n    \n    print('Training....')\n    \n    temp['train'] = model.fit(trainx[train_ind],trainy[train_ind],batch_size=b_size,epochs=500, verbose=0)\n    \n    print(f'Training score to beat: {laplace_log_likelihood(trainy[train_ind],np.mean(trainy[train_ind]),np.std(trainy[train_ind]))}')\n    \n    print(model.evaluate(trainx[train_ind],trainy[train_ind],batch_size=b_size))\n    \n    print(f'Validation score to beat: {laplace_log_likelihood(trainy[val_ind],np.mean(trainy[val_ind]),np.std(trainy[val_ind]))}')\n    temp['test'] = model.evaluate(trainx[val_ind],trainy[val_ind],batch_size=b_size,verbose=1)\n    \n    history[count]=temp\n    \n    count+=1\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x=model.predict(trainx[val_ind])\nx=x.mean(axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(x)\nplt.plot(trainy[val_ind])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import sklearn\nsklearn.metrics.r2_score(trainy[val_ind], x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"trainy[train_ind][0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred= np.array([[1990,2000,2010],\n                  [1980,2002,2010]])\n\n\nprint(y_pred[1] - y_pred[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.summary()","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}